ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/freemol/trunk/smi23d/src/mengine/read.c
Revision: 3
Committed: Mon Jun 9 21:38:26 2008 UTC (12 years, 5 months ago) by tjod
File size: 21581 byte(s)
Log Message:
test

Line User Rev File contents
1 tjod 3 #define EXTERN extern
2    
3     #include "pcwin.h"
4     #include "pcmod.h"
5     #include "substr.h"
6     #include "field.h"
7     #include "pot.h"
8     #include "pdb.h"
9    
10     #include <errno.h>
11    
12     #define MAIN 0L
13     #define SUBSTRUCTURE 1L
14     #define MOPAC 0L
15     #define MNDO 1L
16    
17     int FetchRecord(FILE *, char *);
18     FILE * fopen_path ( char * , char * , char * ) ;
19     int is_bond(int, int);
20     int get_bondorder(int,int);
21     void substr_scale(void);
22     void writefile(int);
23     void readpid(char *);
24     void readfile(int,int);
25     void check_numfile(int);
26     void inxyz(int);
27     void incsd(int,int);
28     void mopaco(int,int);
29     void find_termini(int *,int *);
30     void find_nucterm(int *,int *);
31     void find_cp(int *);
32     void type(void);
33     void hdel(int);
34     void hcoord(void);
35     void pcmfout(int);
36     void pcmfin(int,int);
37     void mmp22mod(int,int);
38     void mm32mod(int);
39     void mod2mmp2(int);
40     void readz(int, int);
41     void mdl2mod(int);
42     void mod2c3d(void);
43     void mod2tinker(void);
44     void c3d2mod(int,int *);
45     void tink2mod(int,int *);
46     void readz(int, int);
47     void read_arcfile(int, int);
48     void write_mopac(void);
49     void alc2mod(int,int *);
50     void mod2alc(void);
51     void mod2sybyl(void);
52     void mod2sybyl2(void);
53     void sybyl2mod(int ,int *);
54     void mod2mac(void);
55     void mac2mod(int, int);
56     void read_psgvbin(int);
57     void read_psgvbout(int);
58     void write_psgvbin(void);
59     void read_freq_log(int,int);
60     void read_games(int,int);
61     void read_hondo(int,int);
62     void initialize(void);
63     void write_mm2(void);
64     void write_mm3(void);
65     void mod2pdb(void);
66     void make_bond(int,int,int);
67     void resetsubstrmem(void);
68     void subconnect(int,int);
69     void mark_hbond(void);
70     void set_atomtypes(int);
71     void mm2gaus(void);
72     void mm2games(void);
73     void write_hondo(void);
74     int read_sdf(int,int);
75     void write_sdf(void);
76     void message_alert(char *, char *);
77     void write_eht(void);
78     void read_gausfchk(int);
79     void read_tm(int,int);
80     void mod2tm(void);
81     void read_adf(int,int);
82     void write_adf(void);
83     void read_smiles(int,int);
84     void write_smiles(void);
85     void build_coord(void);
86    
87     EXTERN struct t_files {
88     int nfiles, append, batch, icurrent;
89     int istrselect;
90     } files;
91    
92     int filetypes[10000];
93    
94     int oh2, oh3, oh6;
95     void clean_string(char *);
96     // ========================================
97     void clean_string(char *astring)
98     {
99     int iz,i;
100     iz = strlen(astring);
101     for (i=0; i < iz; i++)
102     {
103     if (isprint(astring[i]) == 0)
104     astring[i] = ' ';
105     }
106     for (i = iz -1; i >= 0; i--)
107     {
108     if (astring[i] != ' ')
109     {
110     astring[i+1] = '\0';
111     break;
112     }
113     }
114     }
115     /*======================================================================*/
116     void * malloc_filename ( char *path , char *name )
117     {
118     int ix, iz;
119     char *filename;
120    
121     ix = 0;
122     if (path != NULL)
123     ix = strlen(path);
124     iz = strlen(name);
125     filename = malloc (ix + 1 + iz + 1) ;
126     if ( ix == 0 )
127     {
128     strcpy (filename,name);
129     }
130     else
131     {
132     sprintf (filename,"%s\\%s",path,name) ;
133     }
134     return ( (void *) filename ) ;
135     }
136     /*======================================================================*/
137     FILE * fopen_path ( char *path , char *name , char *mode )
138     {
139     char *filename = malloc_filename(path,name) ;
140     FILE *rval = fopen (filename,mode) ;
141     if ( rval == NULL )
142     {
143     fprintf (pcmoutfile,"open of --%s-- failed\n",filename) ;
144     fprintf(pcmoutfile,"Error %d %s\n",errno,strerror(errno));
145     fclose(pcmoutfile);
146     exit(0);
147     }
148     free (filename) ;
149     return ( rval ) ;
150     }
151    
152     /* ============================================= */
153     void readfile(int isel, int isubred)
154     {
155     int iftype, i;
156     int oldfield;
157    
158     /* need to determine number of structures in file
159     which structure the user wants to read and then
160     call the appropriate read routine. Check to see if
161     we already have structure and erase if necessary */
162    
163    
164     if (isubred == 0)
165     initialize();
166    
167    
168     iftype = Openbox.ftype;
169     if (isubred == 1)
170     {
171     for (i=0; i < MAXSS; i++)
172     {
173     if ( !substr.istract[i] )
174     break;
175     }
176     substr.icurstr = i;
177     substr.istract[i] = TRUE;
178     }
179    
180     switch(Openbox.ftype)
181     {
182     case FTYPE_PCM:
183     oldfield = field.type;
184     if (field.type != MMX)
185     field.type = MMX;
186     pcmfin(isel,isubred);
187     field.type = oldfield;
188     break;
189     case FTYPE_SDF: // need to check number of structures in file
190     oldfield = field.type;
191     if (field.type != MMX)
192     field.type = MMX;
193     i = read_sdf(isel,isubred);
194     field.type = oldfield;
195     break;
196     }
197     /* move substructure to right of main structures */
198     if (isubred == 1)
199     substr_scale();
200    
201     // convert to correct types for current force field
202     type();
203     for(i=1; i <= natom; i++)
204     {
205     if (atom[i].mmx_type == 5)
206     {
207     flags.noh = TRUE;
208     break;
209     }
210     }
211     substr.icurstr = -1;
212     return;
213     }
214     /* ================================================== */
215     void writefile(int isubred)
216     {
217     int oldfield;
218    
219     oldfield = field.type;
220     if (Savebox.ftype == FTYPE_MM3 && field.type == MMX)
221     {
222     // convert to mm3 types
223     set_atomtypes(MM3);
224     field.type = MM3;
225     type();
226     } else if (field.type == MM3)
227     {
228     // convert to MMX types
229     set_atomtypes(MMX);
230     field.type = MMX;
231     type();
232     }
233     if (isubred == MAIN)
234     {
235     switch(Savebox.ftype)
236     {
237     case FTYPE_PCM:
238     pcmfout(1);
239     break;
240     case FTYPE_SDF:
241     write_sdf();
242     break;
243     }
244     }
245     }
246     /* ------------------------------------------------- */
247     void check_numfile(int ftype)
248     {
249     int igau, header, model;
250     FILE *infile;
251     char inputline[160];
252    
253     header = FALSE;
254     model = FALSE;
255     infile = fopen_path(Openbox.path,Openbox.fname,"r");
256     if (infile == NULL)
257     {
258     message_alert("Error trying to check number of structures in file. Probably a bad filename.","Check Numfile");
259     fprintf(pcmoutfile,"Error in check numfile:: %s %s\n",Openbox.path, Openbox.fname);
260     files.nfiles = 0;
261     return;
262     }
263     igau = 0;
264     while ( fgets(inputline,159,infile) != NULL )
265     {
266     if (strncasecmp(inputline,"{PCM",4) == 0)
267     {
268     filetypes[files.nfiles] = FTYPE_PCM;
269     files.nfiles++;
270     }else if (strncasecmp(inputline,"SDF ",4) == 0)
271     {
272     filetypes[files.nfiles] = FTYPE_SDF;
273     files.nfiles++;
274     }
275     }
276     fclose(infile);
277     }
278     /* ------------------------------------*/
279     int get_bondorder(int ia,int ib)
280     {
281     int i;
282     for (i=0; i < MAXIAT; i++)
283     {
284     if (atom[ia].iat[i] == ib)
285     return atom[ia].bo[i];
286     }
287     return 0;
288     }
289     /* ------------------------------------*/
290     int is_bond(int ia, int ib)
291     {
292     int i;
293     for (i=0; i < MAXIAT; i++)
294     {
295     if (atom[ia].iat[i] == ib)
296     return TRUE;
297     }
298     return FALSE;
299     }
300     /* ------------------------------------*/
301     void mopaco(int istart,int ifinish)
302     {
303     int i, ibi, ibj, it, j, jt, kk;
304     float disij, dx, dy, dz, r, xii, yii, zii,rbdis[4];
305     static float radii[110] =
306     { 0.77, 0.77, 0.77, 0.77, 0.50, 0.66, 0.66, 0.70, 0.70, 0.70,
307     0.64, 0.99, 1.14, 1.33, 1.04, 1.04, 1.04, 1.04, 1.17, 0.01,
308     0.50, 0.77, 0.50, 0.50, 1.10, 0.88, 0.88, 0.50, 0.77, 0.77,
309     1.22, 1.40, 1.20, 1.17, 1.37, 0.50, 0.70, 1.04, 1.17, 0.77,
310     0.70, 0.66, 0.77, 0.70, 0.50, 0.66, 1.10, 0.77, 0.77, 0.77,
311     0.77, 0.66, 0.70, 1.33, 0.66, 0.77, 0.77, 0.70, 0.01, 0.50,
312     0.77, 0.77, 0.77, 0.70, 0.66, 0.70, 1.10, 0.66, 0.70, 0.50,
313     0.77, 0.77, 0.77, 0.77, 0.50, 0.66, 0.66, 0.70, 0.70, 0.70,
314     0.77, 0.77, 0.77, 0.77, 0.50, 0.66, 0.66, 0.70, 0.70, 0.70,
315     0.77, 0.77, 0.77, 0.77, 0.50, 0.66, 0.66, 0.70, 0.70, 0.70,
316     0.77, 0.77, 0.77, 0.77, 0.50, 0.66, 0.66, 0.70, 0.70, 0.70,};
317    
318     /* SUBROUTINE SENDS BACK THE ACTUAL NUMBER OF ATOMS, AND
319     * DETERMINES THE BOND CONNECTIONS and bond orders TO ESTABLISH THE attached
320     * ATOM LIST. ATOM TYPES HAVE BEEN CHANGED TO MM2 TYPES. THOSE THAT
321     * AREN'T IN MM2 ARE GIVEN #S = AT # + 60
322     * BOND ORDERS AND THE IATMS AND IBNDS LISTS ARE GENERATED
323     * NATOMS IS THE # OF ENTRIES IN THE Z MATRIX */
324    
325     rbdis[0] = 1.42;
326     rbdis[1] = 1.25;
327     rbdis[2] = 1.37;
328     rbdis[3] = 1.25;
329    
330     /* find all single bonds based on sum of covalent radii */
331     for( i = istart; i <= (ifinish - 1L); i++ )
332     {
333     it = atom[i].mmx_type;
334     for( j = i + 1L; j <= ifinish; j++ )
335     {
336     if (!is_bond(i,j))
337     {
338     jt = atom[j].mmx_type;
339     dx = atom[i].x - atom[j].x;
340     dy = atom[i].y - atom[j].y;
341     dz = atom[i].z - atom[j].z;
342     r = sqrt( dx*dx + dy*dy + dz*dz );
343     if( it >= 80L && jt >= 80L )
344     {
345     if( r <= (atom[i].radius + atom[j].radius) )
346     make_bond( i, j, 1 );
347     }else if( it >= 80L && jt < 80L )
348     {
349     if( r <= (atom[i].radius + radii[atom[j].mmx_type-1]) )
350     make_bond( i, j, 1);
351    
352     }else if( it < 80L && jt >= 80L )
353     {
354     if( r <= (radii[atom[i].mmx_type-1] + atom[j].radius + .1) )
355     make_bond( i, j, 1 );
356     }else if( r <= (radii[atom[i].mmx_type-1] + radii[atom[j].mmx_type-1]+ .1) )
357     {
358     if( i < j )
359     make_bond( i, j, 1 );
360     }
361     }
362     }
363     }
364    
365     /* find all double and triple bonds - search only types for
366     * carbon, oxygen, nitrogen, sulfur, phosphorous */
367     for( i = istart; i <= (ifinish - 1); i++ )
368     {
369     it = atom[i].atomnum;
370     if( it == 6 || it == 7 || it == 8 || it == 15 || it == 16 )
371     {
372     for( j = 0; j < MAXIAT; j++ )
373     {
374     if( atom[i].iat[j] != 0 )
375     {
376     jt = atom[atom[i].iat[j]].atomnum;
377     if( jt == 6 || jt == 7 || jt == 8 || jt == 15 || jt == 16 )
378     {
379     xii = atom[i].x - atom[atom[i].iat[j]].x;
380     yii = atom[i].y - atom[atom[i].iat[j]].y;
381     zii = atom[i].z - atom[atom[i].iat[j]].z;
382     if( fabs( xii ) <= 2.3 )
383     {
384     if( fabs( yii ) <= 2.3 )
385     {
386     if( fabs( zii ) <= 2.3 )
387     {
388     disij = sqrt( xii*xii + yii*yii + zii*zii );
389     /* c=o bond */
390     if( ((it == 8 && jt == 6) || (it == 6 && jt == 8)) && disij <= rbdis[3] + 0.05 )
391     {
392     if( i < atom[i].iat[j] )
393     make_bond( i, atom[i].iat[j],1 );
394     /* alkyne */
395     }else if( (it == 6 && jt == 6) && disij <= rbdis[1] )
396     {
397     if( i < atom[i].iat[j] )
398     make_bond( i, atom[i].iat[j],2 );
399     /* c=n */
400     }else if( ((it == 7 && jt == 6) || (it == 6 && jt == 7))
401     && disij <= rbdis[2] )
402     {
403     ibi = 0;
404     for( kk = 0; kk < MAXIAT; kk++ )
405     {
406     if( atom[i].bo[kk] != 9 )
407     ibi = ibi + atom[i].bo[kk];
408     }
409     ibj = 0;
410     for( kk = 0; kk < MAXIAT; kk++ )
411     {
412     if( atom[atom[i].iat[j]].bo[kk] != 9 )
413     ibj = ibj + atom[atom[i].iat[j]].bo[kk];
414     }
415     if( it == 7 )
416     {
417     if( ibi <= 2 && ibj <= 3 )
418     {
419     if( i < atom[i].iat[j] )
420     make_bond( i, atom[i].iat[j], 1 );
421     break;
422     }
423     }else if( jt == 7 )
424     {
425     if( ibj <= 2 && ibi <= 3 )
426     {
427     if( i < atom[i].iat[j] )
428     make_bond( i, atom[i].iat[j], 1 );
429     break;
430     }
431     }
432     /*nitrile */ }else if( ((it == 7 && jt == 6) || (it == 6 && jt == 7)) && disij <= 1.18 )
433     {
434     if( i < atom[i].iat[j] )
435     make_bond( i, atom[i].iat[j], 2 );
436     /* N=O */ }else if ( ((it == 7 && jt == 8) || (it == 8 && jt == 7)) && disij <= 1.26)
437     {
438     ibi = 0;
439     if (it == 7)
440     {
441     for (kk = 0; kk < MAXIAT; kk++)
442     {
443     if (atom[i].bo[kk] != 9)
444     ibi += atom[i].bo[kk];
445     }
446     } else
447     {
448     for (kk = 0; kk < MAXIAT; kk++)
449     {
450     if (atom[atom[i].iat[j]].bo[kk] != 9)
451     ibi += atom[atom[i].iat[j]].bo[kk];
452     }
453     }
454     if (ibi <= 3)
455     {
456     if ( i < atom[i].iat[j])
457     make_bond(i,atom[i].iat[j],1);
458     break;
459     }
460     /* c=c bond */ }else if( (it == 6 && jt == 6) && disij <= rbdis[0] )
461     {
462     ibi = 0;
463     for( kk = 0; kk < MAXIAT; kk++ )
464     {
465     if( atom[i].bo[kk] != 9 )
466     ibi = ibi + atom[i].bo[kk];
467     }
468     ibj = 0;
469     for( kk = 0; kk < MAXIAT; kk++ )
470     {
471     if( atom[atom[i].iat[j]].bo[kk] != 9 )
472     ibj = ibj + atom[atom[i].iat[j]].bo[kk];
473     }
474     if( ibi <= 3 && ibj <= 3 )
475     {
476     if( i < atom[i].iat[j] )
477     make_bond( i, atom[i].iat[j], 1 );
478     }
479     }
480     }
481     }
482     }
483     }
484     }
485     }
486     }
487     }
488     type();
489     return;
490     }
491     void find_nucterm(int *i5p,int *i3p)
492     {
493     int i, icount;
494     long int mask1,mask2;
495     // assume find first 3p end and second 5p
496     mask1 = 1L << P5; /* 5 prime terminus */
497     mask2 = 1L << P3; /* 3 prime terminus */
498    
499     icount = 0;
500     for (i=1; i <= natom; i++)
501     {
502     if (atom[i].flags & mask1)
503     {
504     icount++;
505     if (icount == 2)
506     {
507     *i5p = i;
508     break;
509     }
510     }
511     }
512     for (i=1; i <= natom; i++)
513     {
514     if (atom[i].flags & mask2)
515     {
516     *i3p = i;
517     break;
518     }
519     }
520     }
521     void find_termini(int *iNT,int *iCT)
522     {
523     long int i, j, mask1, mask2;
524    
525     mask1 = 1L << NTERM; /* N terminus */
526     mask2 = 1L << CNTERM; /* C terminus */
527    
528     for (i = 1; i < natom; i++)
529     {
530     if (atom[i].flags & mask2)
531     {
532     *iCT = i;
533     for (j = i+1; j < natom; j++)
534     {
535     if (atom[j].flags & mask1)
536     {
537     *iNT = j;
538     break;
539     }
540     }
541     break;
542     }
543     }
544     return;
545     }
546     void find_cp(int *iCP)
547     {
548     long int i, mask1;
549    
550     mask1 = 1L << DUMMY; /* connect point */
551    
552     for (i = 1; i < natom; i++)
553     {
554     if (atom[i].flags & mask1)
555     {
556     *iCP = i;
557     break;
558     }
559     }
560     return;
561     }
562    
563     void substr_scale()
564     {
565     int isub,ii;
566     float xmain_min,xmain_max, sub_xmin, xdif, sub_center, main_center;
567     float ymain_min,ymain_max, sub_ymax, sub_ymin, ydif;
568    
569     xmain_min = 1000.0F;
570     ymain_min = 1000.0F;
571     sub_xmin = 1000.0F;
572     sub_ymin = 1000.0F;
573    
574     xmain_max = -1000.0F;
575     ymain_max = -1000.0F;
576     sub_ymax = -1000.0F;
577    
578     isub = (1L << substr.icurstr);
579    
580     for (ii = 1; ii <= natom; ii++)
581     {
582     if (atom[ii].substr[0] == isub)
583     {
584     if ( atom[ii].x < sub_xmin )
585     sub_xmin = atom[ii].x;
586     if ( atom[ii].y < sub_ymin )
587     sub_ymin = atom[ii].y;
588     if ( atom[ii].y > sub_ymax )
589     sub_ymax = atom[ii].y;
590     } else
591     {
592     if ( atom[ii].x < xmain_min )
593     xmain_min = atom[ii].x;
594     if (atom[ii].x > xmain_max)
595     xmain_max = atom[ii].x;
596     if ( atom[ii].y < ymain_min )
597     ymain_min = atom[ii].y;
598     if (atom[ii].y > ymain_max)
599     ymain_max = atom[ii].y;
600     }
601     }
602    
603     if (xmain_max > sub_xmin)
604     {
605     xdif = xmain_max - sub_xmin;
606     for (ii = 1; ii <= natom; ii++)
607     {
608     if (atom[ii].substr[0] == isub)
609     atom[ii].x += xdif;
610     }
611     }
612    
613     main_center = ((ymain_max - ymain_min)/2.0F) + ymain_min;
614     sub_center = ((sub_ymax - sub_ymin)/2.0F) + sub_ymin;
615    
616     ydif = main_center - sub_center;
617    
618     for (ii = 1; ii <= natom; ii++)
619     {
620     if (atom[ii].substr[0] == isub)
621     atom[ii].y += ydif;
622     }
623     }
624     /* ============================================= */
625     int FetchRecord(FILE *fp, char *buffer)
626     {
627     int ch;
628     char *ptr;
629    
630     ptr = buffer;
631     do
632     {
633     ch = getc(fp);
634     if (ch >= ' ')
635     *ptr++ = ch;
636     else if (ch == '\n')
637     {
638     *ptr = 0;
639     return TRUE;
640     } else if (ch == '\r')
641     {
642     ch = getc(fp);
643     if (ch != '\n')
644     ungetc(ch,fp);
645     *ptr = 0;
646     return TRUE;
647     } else if (ch == EOF)
648     {
649     *ptr = 0;
650     return (ptr != buffer);
651     } else
652     *ptr++ = ch;
653     } while (ptr < buffer+255);
654    
655     do
656     {
657     ch = getc(fp);
658     } while (ch !='\n' && ch != '\r' && ch != EOF);
659     if (ch == '\r')
660     {
661     ch = getc(fp);
662     if (ch != '\n')
663     ungetc(ch,fp);
664     }
665     *ptr = 0;
666     return TRUE;
667     }
668    
669    
670    
671