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