ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/freemol/trunk/src/mengine/src/read_sdf.c
Revision: 31
Committed: Fri Jul 18 18:47:27 2008 UTC (13 years, 10 months ago) by tjod
File size: 17021 byte(s)
Log Message:
When aromatic bonds are encountered generate error on stderr and do not output sdf file.

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 tjod 31 // TJO
259     Ret_Val = FALSE;
260 tjod 3 }
261     make_bond(ia1, ia2, ibond);
262     if (istereo == 1 && planar)
263     {
264     atom[ia2].z += 0.7;
265     if (atom[ia2].atomnum == 1)
266     atom[ia2].type = 60;
267     }
268     if (istereo == 6 && planar)
269     {
270     atom[ia2].z -= 0.7;
271     if (atom[ia2].atomnum == 1)
272     atom[ia2].type = 60;
273     }
274     }
275     // read to end of structure
276     while (FetchRecord(rfile,inputline))
277     {
278     if (strncasecmp(inputline,"$$$$",4) == 0)
279     {
280     return Ret_Val;
281     }
282     }
283     return Ret_Val;
284     }
285     /* =================================== */
286     int read_sdf(int istruct, int isubred)
287     {
288     int i, j, niatom, nibond, ia1, ia2, ibond, itype, newatom;
289     int ncount, ibotptr,junk,junk1,got_title,istereo,iz;
290     int jji, jj1, jj2, jj3, jj4;
291     int has_Aromatic, xPlus,yPlus,zPlus, planar;
292     int icount,itemp[4];
293     float xtmp, ytmp, ztmp, dz;
294     char c1[4],c2[4];
295     char atomchar[3];
296     char inputline[150];
297     FILE *infile;
298    
299     infile = fopen_path(Openbox.path,Openbox.fname,"r");
300    
301     if (infile == NULL)
302     {
303     message_alert("Error opening Concord file","Concord Setup");
304     return FALSE;
305     }
306    
307     if (isubred == 0)
308     ibotptr = 0;
309     else
310     {
311     ibotptr = natom;
312     substr.istract[substr.icurstr] = TRUE;
313     }
314    
315     if (istruct > 1)
316     {
317     ncount = 0;
318     while (FetchRecord(infile,inputline))
319     {
320     if (strncasecmp(inputline,"$$$$",4) == 0)
321     {
322     ncount++;
323     if (ncount+1 == istruct)
324     goto L_1;
325     }
326     }
327     }
328     // start file read here
329     L_1:
330     got_title = FALSE;
331     xPlus = yPlus = zPlus = FALSE;
332     FetchRecord(infile,inputline);
333     sscanf(inputline,"SDF %s",Struct_Title);
334     got_title = TRUE;
335     /* if (strlen(inputline) > 4)
336     {
337     iz = strlen(inputline);
338     if (iz > 60) iz = 59;
339     for (i=4; i < iz; i++)
340     Struct_Title[i-4] = inputline[i];
341     Struct_Title[i] = '\0';
342     got_title = TRUE;
343     } */
344     FetchRecord(infile,inputline); // blank line
345     FetchRecord(infile,inputline); // blank line
346    
347     FetchRecord(infile,inputline); // natom and nbond
348     for (i=0; i <4; i++)
349     {
350     c1[i] = inputline[i];
351     c2[i] = inputline[i+3];
352     }
353     c1[3] = '\0';c2[3] = '\0';
354     niatom = atoi(c1);
355     nibond = atoi(c2);
356     if (niatom == 0)
357     return FALSE;
358    
359     for (i=0; i < niatom; i++)
360     {
361     FetchRecord(infile,inputline);
362     sscanf(inputline,"%f %f %f %s %d %d",&xtmp, &ytmp, &ztmp, atomchar,&junk,&junk1);
363    
364     itype = 0;
365     if (xtmp != 0.0) xPlus= TRUE;
366     if (ytmp != 0.0) yPlus= TRUE;
367     if (ztmp != 0.0) zPlus= TRUE;
368    
369     iz = strlen(atomchar);
370     if (itype == 0 && (atomchar[0] == 'N' || atomchar[0] == 'O') && iz == 1) // nitro fix
371     {
372     if (atomchar[0] == 'N') itype = 8;
373     if (atomchar[0] == 'O') itype = 6;
374    
375     if (junk1 != 0)
376     {
377     if (junk1 == 3 && itype == 8)
378     itype = 41; // N+
379     if (junk1 == 5 && itype == 6)
380     itype = 42; // O-
381     }
382     }
383    
384     newatom = make_atom(itype,xtmp,ytmp,ztmp,atomchar);
385     if (newatom == -1)
386     {
387     printf("Atom %d %s not recognized structure skipped\n",i,atomchar);
388     return FALSE;
389     }
390     if (isubred == 1)
391     {
392     atom[newatom].substr[0] = 0;
393     atom[newatom].substr[0] |= (1L << substr.icurstr);
394     }
395     }
396     has_Aromatic = FALSE;
397    
398     planar = TRUE;
399     if (xPlus && yPlus && zPlus) planar = FALSE;
400    
401     for (i=0; i < nibond; i++)
402     {
403     FetchRecord(infile,inputline);
404     for (j=0; j <4; j++)
405     {
406     c1[j] = inputline[j];
407     c2[j] = inputline[j+3];
408     }
409     c1[3] = '\0';c2[3] = '\0';
410     ia1 = atoi(c1);
411     ia2 = atoi(c2);
412     istereo = 0;
413     sscanf(inputline,"%3d%3d%3d%3d",&junk, &junk1, &ibond,&istereo);
414     if (ibond >= 4)
415     {
416     ibond = 1;
417     has_Aromatic = TRUE;
418     }
419     make_bond(ia1+ibotptr, ia2+ibotptr, ibond);
420     if (istereo == 1 && planar)
421     {
422     atom[ia2+ibotptr].z += 0.7;
423     if (atom[ia2+ibotptr].atomnum == 1)
424     atom[ia2+ibotptr].type = 60;
425     }
426     if (istereo == 6 && planar)
427     {
428     atom[ia2+ibotptr].z -= 0.7;
429     if (atom[ia2+ibotptr].atomnum == 1)
430     atom[ia2+ibotptr].type = 60;
431     }
432     }
433     FetchRecord(infile,inputline); // M END line
434     FetchRecord(infile,inputline);
435     if (got_title == FALSE)
436     strcpy(Struct_Title,inputline);
437     fclose(infile);
438     ncount = strlen(Struct_Title);
439     for (i=0; i < ncount; i++)
440     {
441     if (Struct_Title[i] == '\n')
442     {
443     Struct_Title[i] = '\0';
444     break;
445     }
446     }
447     // search for quaternary centers - check if flat
448     dz = 0.0;
449     for (i=1; i <= natom; i++)
450     {
451     dz += atom[i].z;
452     jji = jj1 = jj2 = jj3 = jj4 = 0;
453     if (atom[i].type != 0)
454     {
455     for (j=0; j < 6; j++)
456     {
457     if (atom[i].iat[j] != 0)
458     jji++;
459     }
460     }
461     if (jji == 4)
462     {
463     if (atom[atom[i].iat[0]].z == 0.0 && atom[atom[i].iat[0]].z == 0.0 &&
464     atom[atom[i].iat[0]].z == 0.0 && atom[atom[i].iat[0]].z == 0.0) // flat center
465     {
466     for (j=0; j < 4; j++)
467     {
468     itemp[j] = 0;
469     if (atom[atom[i].iat[0]].iat[j] != 0) jj1++;
470     if (atom[atom[i].iat[1]].iat[j] != 0) jj2++;
471     if (atom[atom[i].iat[2]].iat[j] != 0) jj3++;
472     if (atom[atom[i].iat[3]].iat[j] != 0) jj4++;
473     }
474     icount = 0;
475     if (jj1 == 1)
476     {
477     itemp[icount] = atom[i].iat[0];
478     icount++;
479     }
480     if (jj2 == 1)
481     {
482     itemp[icount] = atom[i].iat[1];
483     icount++;
484     }
485     if (jj3 == 1)
486     {
487     itemp[icount] = atom[i].iat[2];
488     icount++;
489     }
490     if (jj3 == 1)
491     {
492     itemp[icount] = atom[i].iat[3];
493     icount++;
494     }
495     if (icount >= 2)
496     {
497     atom[itemp[0]].z += 0.5;
498     atom[itemp[1]].z -= 0.5;
499     }
500     }
501     }
502     }
503     for (i=1; i < natom; i++)
504     {
505     for (j=i+1; j <= natom; j++)
506     {
507     xtmp = atom[i].x - atom[j].x;
508     ytmp = atom[i].y - atom[j].y;
509     ztmp = atom[i].z - atom[j].z;
510     if ( (xtmp + ytmp + ztmp) < 0.1)
511     {
512     atom[i].x += 0.1;
513     atom[i].y += 0.1;
514     atom[i].z += 0.1;
515     atom[j].x -= 0.1;
516     atom[j].y -= 0.1;
517     atom[j].z -= 0.1;
518     }
519     }
520     }
521    
522     if (has_Aromatic == TRUE)
523     mopaco(1,natom);
524     return TRUE;
525     }
526     // ==================================================
527     void read_schakal(int isel,int isubred)
528     {
529     int i, itype, newatom, iz;
530     float xtmp, ytmp, ztmp;
531     char dummy[30];
532     char atomchar[6];
533     char inputline[150];
534     FILE *infile;
535    
536     infile = fopen_path(Openbox.path,Openbox.fname,"r");
537    
538     if (infile == NULL)
539     {
540     message_alert("Error opening Schakal file","Schakal Setup");
541     return;
542     }
543     while ( FetchRecord(infile,inputline))
544     {
545     sscanf(inputline,"%s",dummy);
546     if (strcmp(dummy,"TITLE") == 0)
547     {
548     sscanf(inputline,"%s %s",dummy,Struct_Title);
549     } else if (strcmp(dummy,"ATOM") == 0)
550     {
551     sscanf(inputline,"%s %s %f %f %f",dummy,atomchar,&xtmp,&ytmp,&ztmp);
552     // strip off numbers
553     iz = strlen(atomchar);
554     for (i=0; i < iz; i++)
555     {
556     if (isdigit(atomchar[i]))
557     {
558     atomchar[i] = '\0';
559     break;
560     }
561     }
562     itype = 0;
563     newatom = make_atom(itype,xtmp,ytmp,ztmp,atomchar);
564     } else if (strcmp(dummy,"END") == 0)
565     {
566     break;
567     }
568     }
569     fclose(infile);
570     mopaco(1,natom);
571     //
572     }
573     // ==================================================
574     void write_schakal()
575     {
576     int i,j,lptest;
577     char atomchar[6];
578     FILE *wfile;
579    
580     lptest = 0;
581     for( i = 1; i <= natom; i++ )
582     {
583     if( atom[i].mmx_type == 20 )
584     {
585     lptest = 1;
586     hdel( lptest );
587     break;
588     }
589     }
590    
591     /* now write the schakal file */
592     if (files.append)
593     wfile = fopen_path(Savebox.path,Savebox.fname,"a");
594     else
595     wfile = fopen_path(Savebox.path,Savebox.fname,"w");
596    
597     j = strlen(Struct_Title);
598     for (i=0; i < j; i++)
599     {
600     if (Struct_Title[i] == '\n')
601     {
602     Struct_Title[i] = '\0';
603     break;
604     }
605     }
606     fprintf(wfile,"TITLE %s\n",Struct_Title);
607     fprintf(wfile,"CELL \n");
608     for (i=1; i <= natom; i++)
609     {
610     sprintf(atomchar,"%s%d",Elements[atom[i].atomnum-1].symbol,i);
611    
612     fprintf(wfile,"ATOM %-5s %12.7f%12.7f%12.7f\n",atomchar,atom[i].x,
613     atom[i].y, atom[i].z);
614     }
615     fprintf(wfile,"END \n");
616     fclose(wfile);
617     }