ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/freemol/trunk/src/mengine/src/read_sdf.c
Revision: 27
Committed: Tue Jul 8 19:10:25 2008 UTC (13 years, 5 months ago) by tjod
File size: 16985 byte(s)
Log Message:
move mengine src to reflect new freemol directory structure

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     }