ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/freemol/trunk/src/mengine/src/read_sdf.c
Revision: 48
Committed: Thu Jul 31 16:34:52 2008 UTC (13 years, 2 months ago) by tjod
File size: 17096 byte(s)
Log Message:
Add modification notice to files modified for Freemol
per mengine/smi23d license requirements

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