ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/freemol/trunk/src/mengine/src/read_sdf.c
Revision: 15
Committed: Thu Jun 19 00:26:16 2008 UTC (13 years, 5 months ago) by tjod
Original Path: trunk/smi23d/src/mengine/read_sdf.c
File size: 16985 byte(s)
Log Message:
Create new main program, mengine.c to replace gmmx.c
simplifying command line parsing (none).
Make mengine read stdin (sdf) and write stdout (sdf) instead
of files.  Makefile reflects change from gmmx to mengine.

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 "atom_k.h"
7     #include "energies.h"
8    
9     EXTERN struct t_files {
10     int nfiles, append, batch, icurrent, ibatno;
11     } files;
12     EXTERN struct ElementType { char symbol[3];
13     int atomnum;
14     float weight, covradius, vdwradius;
15     int s,p,d,f,type ;
16     } Elements[];
17     EXTERN struct t_dipolemom {
18     double total, xdipole, ydipole, zdipole;
19     } dipolemom;
20     EXTERN struct t_logp {
21     float logp;
22     } logp_calc;
23     EXTERN struct t_vibdata {
24     char ptgrp[4];
25     float mom_ix,mom_iy,mom_iz;
26     float etot,htot,stot,gtot,cptot;
27     } vibdata;
28    
29     int make_atom(int, float, float, float,char *);
30     void make_bond(int, int, int);
31     int read_sdf(int,int);
32     int rd_sdf(FILE *);
33     void write_sdf(int);
34     void message_alert(char *, char *);
35     void hdel(int);
36     FILE * fopen_path ( char * , char * , char * ) ;
37     void read_schakal(int,int);
38     void write_schakal(void);
39     void mopaco(int,int);
40     int FetchRecord(FILE *, char *);
41     void avgleg(void);
42    
43     /*
44     * flags - indicates whether dipole, xlogp or vibrational calculations were done
45     * and if so write them to the output file. Look at the definitions in pcmod.h
46     *
47     */
48     void write_sdf(int flags)
49     {
50     int i,j,lptest,nbond,junk;
51     FILE *wfile;
52    
53     lptest = 0;
54     for( i = 1; i <= natom; i++ )
55     {
56     if( atom[i].mmx_type == 20 )
57     {
58     lptest = 1;
59     hdel( lptest );
60     break;
61     }
62     }
63     nbond = 0;
64     /* ** calculate the number of bonds in the molecule ** */
65     for( j = 1; j <= natom; j++ )
66     {
67     for( i = 0; i < MAXIAT; i++ )
68     {
69     if( atom[j].iat[i] != 0 )
70     {
71     if( atom[j].iat[i] < j )
72     nbond = nbond + 1;
73     }
74     }
75     }
76     /* now write the concord file */
77 tjod 15 /*
78 tjod 3 if (files.append)
79     wfile = fopen_path(Savebox.path,Savebox.fname,"a");
80     else
81     wfile = fopen_path(Savebox.path,Savebox.fname,"w");
82 tjod 15 */
83     wfile = stdout;
84 tjod 3
85     j = strlen(Struct_Title);
86     for (i=0; i < j; i++)
87     {
88     if (Struct_Title[i] == '\n')
89     {
90     Struct_Title[i] = '\0';
91     break;
92     }
93     }
94     fprintf(wfile,"%s\n",Struct_Title);
95     fprintf(wfile," PCMODEL v9.1 1.00000 0.00000\n");
96     fprintf(wfile,"\n");
97    
98     fprintf(wfile,"%3d%3d 0 0 0 0 1 V2000\n",natom,nbond);
99    
100     for (i=1; i <= natom; i++)
101     {
102     junk = 0;
103     if (atom[i].mmx_type == 41) junk = 3;
104     else if (atom[i].mmx_type == 42) junk = 5;
105     else if (atom[i].mmx_type == 66)
106     {
107     if (atom[i].bo[0] == 1)
108     junk = 5;
109     }
110     fprintf(wfile,"%10.4f%10.4f%10.4f %-3s 0 %d 0 0 0 0 \n",atom[i].x, atom[i].y,
111     atom[i].z,Elements[atom[i].atomnum-1].symbol,junk);
112     }
113     for (i=1; i <= natom; i++)
114     {
115     for (j=0; j < MAXIAT; j++)
116     {
117     if (atom[i].iat[j] != 0 && i < atom[i].iat[j])
118     fprintf(wfile,"%3d%3d%3d 0\n",i, atom[i].iat[j], atom[i].bo[j]);
119     }
120     }
121     fprintf(wfile,"M END\n");
122     fprintf(wfile,"> <title>\n");
123     fprintf(wfile,"%s\n",Struct_Title);
124    
125     fprintf(wfile,"> <MMFF94 energy>\n");
126     fprintf(wfile,"%f\n",energies.total);
127    
128     if (flags & DO_DIPOLE) {
129     fprintf(wfile,"> <dipole moment>\n");
130     fprintf(wfile,"%f\n",dipolemom.total);
131     }
132    
133     if (flags & DO_XLOGP) {
134     fprintf(wfile,"> <xLogP>\n");
135     fprintf(wfile,"%f\n",logp_calc.logp);
136     }
137    
138     if (flags & DO_VIBRATION) {
139     fprintf(wfile,"> <Point Group>\n");
140     fprintf(wfile,"%s\n",vibdata.ptgrp);
141    
142     fprintf(wfile,"> <Moments of Inertia>\n");
143     fprintf(wfile,"%f8.3 %f8.3 %f8.3 \n",vibdata.mom_ix,vibdata.mom_iy,vibdata.mom_iz);
144    
145     fprintf(wfile,"> <thermodynamics dE, dH, S, dG, CP>\n");
146     fprintf(wfile,"%f8.3 %f8.3 %f8.3 %f8.3 %f8.3 \n",vibdata.etot, vibdata.htot,
147     vibdata.stot,vibdata.gtot,vibdata.cptot);
148     }
149    
150     fprintf(wfile,"\n");
151     fprintf(wfile,"$$$$\n");
152 tjod 15 // fclose(wfile);
153 tjod 3 }
154     /* =================================== */
155     // fast read - assume file is open and positioned
156     // read structure and down to end $$$$
157     // and return
158     int rd_sdf(FILE *rfile)
159     {
160     int i, j, niatom, nibond, ia1, ia2, ibond, itype, newatom;
161     int ncount,junk,junk1,got_title,istereo,iz;
162     int jji, jj1, jj2, jj3, jj4;
163     int has_Aromatic, xPlus,yPlus,zPlus, planar;
164     int icount,itemp[4];
165     int Ret_Val;
166     float xtmp, ytmp, ztmp, dz;
167     char c1[4],c2[4];
168     char atomchar[3];
169     char inputline[150];
170    
171     Ret_Val = TRUE;
172     got_title = FALSE;
173     xPlus = yPlus = zPlus = FALSE;
174 tjod 15 if ( 0 == FetchRecord(rfile,inputline) )return -1;
175     // sscanf(inputline,"SDF %s",Struct_Title);
176     strncpy(Struct_Title, inputline, sizeof(Struct_Title));
177 tjod 3 got_title = TRUE;
178     /* if (strlen(inputline) > 4)
179     {
180     iz = strlen(inputline);
181     if (iz > 60) iz = 59;
182     for (i=4; i < iz; i++)
183     Struct_Title[i-4] = inputline[i];
184     Struct_Title[i] = '\0';
185     got_title = TRUE;
186     } */
187     FetchRecord(rfile,inputline); // blank line
188     FetchRecord(rfile,inputline); // blank line
189    
190     FetchRecord(rfile,inputline); // natom and nbond
191     for (i=0; i <4; i++)
192     {
193     c1[i] = inputline[i];
194     c2[i] = inputline[i+3];
195     }
196     c1[3] = '\0';c2[3] = '\0';
197     niatom = atoi(c1);
198     nibond = atoi(c2);
199     if (niatom == 0)
200     return FALSE;
201    
202     for (i=0; i < niatom; i++)
203     {
204     FetchRecord(rfile,inputline);
205     sscanf(inputline,"%f %f %f %s %d %d",&xtmp, &ytmp, &ztmp, atomchar,&junk,&junk1);
206    
207     itype = 0;
208     if (xtmp != 0.0) xPlus= TRUE;
209     if (ytmp != 0.0) yPlus= TRUE;
210     if (ztmp != 0.0) zPlus= TRUE;
211    
212     iz = strlen(atomchar);
213     if (itype == 0 && (atomchar[0] == 'N' || atomchar[0] == 'O') && iz == 1) // nitro fix
214     {
215     if (atomchar[0] == 'N') itype = 8;
216     if (atomchar[0] == 'O') itype = 6;
217    
218     if (junk1 != 0)
219     {
220     if (junk1 == 3 && itype == 8)
221     itype = 41; // N+
222     if (junk1 == 5 && itype == 6)
223     itype = 42; // O-
224     }
225     }
226    
227     newatom = make_atom(itype,xtmp,ytmp,ztmp,atomchar);
228     if (itype != 0)
229     atom[newatom].mmx_type = itype;
230     if (newatom == -1)
231     {
232     printf("Atom %d %s not recognized structure skipped\n",i,atomchar);
233     Ret_Val = FALSE;
234     }
235     }
236     has_Aromatic = FALSE;
237    
238     planar = TRUE;
239     if (xPlus && yPlus && zPlus) planar = FALSE;
240    
241     for (i=0; i < nibond; i++)
242     {
243     FetchRecord(rfile,inputline);
244     for (j=0; j <4; j++)
245     {
246     c1[j] = inputline[j];
247     c2[j] = inputline[j+3];
248     }
249     c1[3] = '\0';c2[3] = '\0';
250     ia1 = atoi(c1);
251     ia2 = atoi(c2);
252     istereo = 0;
253     sscanf(inputline,"%3d%3d%3d%3d",&junk, &junk1, &ibond,&istereo);
254     if (ibond >= 4)
255     {
256     ibond = 1;
257     has_Aromatic = TRUE;
258     }
259     make_bond(ia1, ia2, ibond);
260     if (istereo == 1 && planar)
261     {
262     atom[ia2].z += 0.7;
263     if (atom[ia2].atomnum == 1)
264     atom[ia2].type = 60;
265     }
266     if (istereo == 6 && planar)
267     {
268     atom[ia2].z -= 0.7;
269     if (atom[ia2].atomnum == 1)
270     atom[ia2].type = 60;
271     }
272     }
273     // read to end of structure
274     while (FetchRecord(rfile,inputline))
275     {
276     if (strncasecmp(inputline,"$$$$",4) == 0)
277     {
278     return Ret_Val;
279     }
280     }
281     return Ret_Val;
282     }
283     /* =================================== */
284     int read_sdf(int istruct, int isubred)
285     {
286     int i, j, niatom, nibond, ia1, ia2, ibond, itype, newatom;
287     int ncount, ibotptr,junk,junk1,got_title,istereo,iz;
288     int jji, jj1, jj2, jj3, jj4;
289     int has_Aromatic, xPlus,yPlus,zPlus, planar;
290     int icount,itemp[4];
291     float xtmp, ytmp, ztmp, dz;
292     char c1[4],c2[4];
293     char atomchar[3];
294     char inputline[150];
295     FILE *infile;
296    
297     infile = fopen_path(Openbox.path,Openbox.fname,"r");
298    
299     if (infile == NULL)
300     {
301     message_alert("Error opening Concord file","Concord Setup");
302     return FALSE;
303     }
304    
305     if (isubred == 0)
306     ibotptr = 0;
307     else
308     {
309     ibotptr = natom;
310     substr.istract[substr.icurstr] = TRUE;
311     }
312    
313     if (istruct > 1)
314     {
315     ncount = 0;
316     while (FetchRecord(infile,inputline))
317     {
318     if (strncasecmp(inputline,"$$$$",4) == 0)
319     {
320     ncount++;
321     if (ncount+1 == istruct)
322     goto L_1;
323     }
324     }
325     }
326     // start file read here
327     L_1:
328     got_title = FALSE;
329     xPlus = yPlus = zPlus = FALSE;
330     FetchRecord(infile,inputline);
331     sscanf(inputline,"SDF %s",Struct_Title);
332     got_title = TRUE;
333     /* if (strlen(inputline) > 4)
334     {
335     iz = strlen(inputline);
336     if (iz > 60) iz = 59;
337     for (i=4; i < iz; i++)
338     Struct_Title[i-4] = inputline[i];
339     Struct_Title[i] = '\0';
340     got_title = TRUE;
341     } */
342     FetchRecord(infile,inputline); // blank line
343     FetchRecord(infile,inputline); // blank line
344    
345     FetchRecord(infile,inputline); // natom and nbond
346     for (i=0; i <4; i++)
347     {
348     c1[i] = inputline[i];
349     c2[i] = inputline[i+3];
350     }
351     c1[3] = '\0';c2[3] = '\0';
352     niatom = atoi(c1);
353     nibond = atoi(c2);
354     if (niatom == 0)
355     return FALSE;
356    
357     for (i=0; i < niatom; i++)
358     {
359     FetchRecord(infile,inputline);
360     sscanf(inputline,"%f %f %f %s %d %d",&xtmp, &ytmp, &ztmp, atomchar,&junk,&junk1);
361    
362     itype = 0;
363     if (xtmp != 0.0) xPlus= TRUE;
364     if (ytmp != 0.0) yPlus= TRUE;
365     if (ztmp != 0.0) zPlus= TRUE;
366    
367     iz = strlen(atomchar);
368     if (itype == 0 && (atomchar[0] == 'N' || atomchar[0] == 'O') && iz == 1) // nitro fix
369     {
370     if (atomchar[0] == 'N') itype = 8;
371     if (atomchar[0] == 'O') itype = 6;
372    
373     if (junk1 != 0)
374     {
375     if (junk1 == 3 && itype == 8)
376     itype = 41; // N+
377     if (junk1 == 5 && itype == 6)
378     itype = 42; // O-
379     }
380     }
381    
382     newatom = make_atom(itype,xtmp,ytmp,ztmp,atomchar);
383     if (newatom == -1)
384     {
385     printf("Atom %d %s not recognized structure skipped\n",i,atomchar);
386     return FALSE;
387     }
388     if (isubred == 1)
389     {
390     atom[newatom].substr[0] = 0;
391     atom[newatom].substr[0] |= (1L << substr.icurstr);
392     }
393     }
394     has_Aromatic = FALSE;
395    
396     planar = TRUE;
397     if (xPlus && yPlus && zPlus) planar = FALSE;
398    
399     for (i=0; i < nibond; i++)
400     {
401     FetchRecord(infile,inputline);
402     for (j=0; j <4; j++)
403     {
404     c1[j] = inputline[j];
405     c2[j] = inputline[j+3];
406     }
407     c1[3] = '\0';c2[3] = '\0';
408     ia1 = atoi(c1);
409     ia2 = atoi(c2);
410     istereo = 0;
411     sscanf(inputline,"%3d%3d%3d%3d",&junk, &junk1, &ibond,&istereo);
412     if (ibond >= 4)
413     {
414     ibond = 1;
415     has_Aromatic = TRUE;
416     }
417     make_bond(ia1+ibotptr, ia2+ibotptr, ibond);
418     if (istereo == 1 && planar)
419     {
420     atom[ia2+ibotptr].z += 0.7;
421     if (atom[ia2+ibotptr].atomnum == 1)
422     atom[ia2+ibotptr].type = 60;
423     }
424     if (istereo == 6 && planar)
425     {
426     atom[ia2+ibotptr].z -= 0.7;
427     if (atom[ia2+ibotptr].atomnum == 1)
428     atom[ia2+ibotptr].type = 60;
429     }
430     }
431     FetchRecord(infile,inputline); // M END line
432     FetchRecord(infile,inputline);
433     if (got_title == FALSE)
434     strcpy(Struct_Title,inputline);
435     fclose(infile);
436     ncount = strlen(Struct_Title);
437     for (i=0; i < ncount; i++)
438     {
439     if (Struct_Title[i] == '\n')
440     {
441     Struct_Title[i] = '\0';
442     break;
443     }
444     }
445     // search for quaternary centers - check if flat
446     dz = 0.0;
447     for (i=1; i <= natom; i++)
448     {
449     dz += atom[i].z;
450     jji = jj1 = jj2 = jj3 = jj4 = 0;
451     if (atom[i].type != 0)
452     {
453     for (j=0; j < 6; j++)
454     {
455     if (atom[i].iat[j] != 0)
456     jji++;
457     }
458     }
459     if (jji == 4)
460     {
461     if (atom[atom[i].iat[0]].z == 0.0 && atom[atom[i].iat[0]].z == 0.0 &&
462     atom[atom[i].iat[0]].z == 0.0 && atom[atom[i].iat[0]].z == 0.0) // flat center
463     {
464     for (j=0; j < 4; j++)
465     {
466     itemp[j] = 0;
467     if (atom[atom[i].iat[0]].iat[j] != 0) jj1++;
468     if (atom[atom[i].iat[1]].iat[j] != 0) jj2++;
469     if (atom[atom[i].iat[2]].iat[j] != 0) jj3++;
470     if (atom[atom[i].iat[3]].iat[j] != 0) jj4++;
471     }
472     icount = 0;
473     if (jj1 == 1)
474     {
475     itemp[icount] = atom[i].iat[0];
476     icount++;
477     }
478     if (jj2 == 1)
479     {
480     itemp[icount] = atom[i].iat[1];
481     icount++;
482     }
483     if (jj3 == 1)
484     {
485     itemp[icount] = atom[i].iat[2];
486     icount++;
487     }
488     if (jj3 == 1)
489     {
490     itemp[icount] = atom[i].iat[3];
491     icount++;
492     }
493     if (icount >= 2)
494     {
495     atom[itemp[0]].z += 0.5;
496     atom[itemp[1]].z -= 0.5;
497     }
498     }
499     }
500     }
501     for (i=1; i < natom; i++)
502     {
503     for (j=i+1; j <= natom; j++)
504     {
505     xtmp = atom[i].x - atom[j].x;
506     ytmp = atom[i].y - atom[j].y;
507     ztmp = atom[i].z - atom[j].z;
508     if ( (xtmp + ytmp + ztmp) < 0.1)
509     {
510     atom[i].x += 0.1;
511     atom[i].y += 0.1;
512     atom[i].z += 0.1;
513     atom[j].x -= 0.1;
514     atom[j].y -= 0.1;
515     atom[j].z -= 0.1;
516     }
517     }
518     }
519    
520     if (has_Aromatic == TRUE)
521     mopaco(1,natom);
522     return TRUE;
523     }
524     // ==================================================
525     void read_schakal(int isel,int isubred)
526     {
527     int i, itype, newatom, iz;
528     float xtmp, ytmp, ztmp;
529     char dummy[30];
530     char atomchar[6];
531     char inputline[150];
532     FILE *infile;
533    
534     infile = fopen_path(Openbox.path,Openbox.fname,"r");
535    
536     if (infile == NULL)
537     {
538     message_alert("Error opening Schakal file","Schakal Setup");
539     return;
540     }
541     while ( FetchRecord(infile,inputline))
542     {
543     sscanf(inputline,"%s",dummy);
544     if (strcmp(dummy,"TITLE") == 0)
545     {
546     sscanf(inputline,"%s %s",dummy,Struct_Title);
547     } else if (strcmp(dummy,"ATOM") == 0)
548     {
549     sscanf(inputline,"%s %s %f %f %f",dummy,atomchar,&xtmp,&ytmp,&ztmp);
550     // strip off numbers
551     iz = strlen(atomchar);
552     for (i=0; i < iz; i++)
553     {
554     if (isdigit(atomchar[i]))
555     {
556     atomchar[i] = '\0';
557     break;
558     }
559     }
560     itype = 0;
561     newatom = make_atom(itype,xtmp,ytmp,ztmp,atomchar);
562     } else if (strcmp(dummy,"END") == 0)
563     {
564     break;
565     }
566     }
567     fclose(infile);
568     mopaco(1,natom);
569     //
570     }
571     // ==================================================
572     void write_schakal()
573     {
574     int i,j,lptest;
575     char atomchar[6];
576     FILE *wfile;
577    
578     lptest = 0;
579     for( i = 1; i <= natom; i++ )
580     {
581     if( atom[i].mmx_type == 20 )
582     {
583     lptest = 1;
584     hdel( lptest );
585     break;
586     }
587     }
588    
589     /* now write the schakal file */
590     if (files.append)
591     wfile = fopen_path(Savebox.path,Savebox.fname,"a");
592     else
593     wfile = fopen_path(Savebox.path,Savebox.fname,"w");
594    
595     j = strlen(Struct_Title);
596     for (i=0; i < j; i++)
597     {
598     if (Struct_Title[i] == '\n')
599     {
600     Struct_Title[i] = '\0';
601     break;
602     }
603     }
604     fprintf(wfile,"TITLE %s\n",Struct_Title);
605     fprintf(wfile,"CELL \n");
606     for (i=1; i <= natom; i++)
607     {
608     sprintf(atomchar,"%s%d",Elements[atom[i].atomnum-1].symbol,i);
609    
610     fprintf(wfile,"ATOM %-5s %12.7f%12.7f%12.7f\n",atomchar,atom[i].x,
611     atom[i].y, atom[i].z);
612     }
613     fprintf(wfile,"END \n");
614     fclose(wfile);
615     }