ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/freemol/trunk/src/mengine/src/read_sdf.c
Revision: 58
Committed: Mon Dec 1 06:44:59 2008 UTC (13 years ago) by wdelano
File size: 13323 byte(s)
Log Message:
code merge 20081130
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 "atom_k.h"
7     #include "energies.h"
8 wdelano 58 #include "fix.h"
9 tjod 3
10 wdelano 58 #define TABULATOR_INCLUDE_IMPLEMENTATION
11     #include "tabulator.h"
12    
13 tjod 3 EXTERN struct t_files {
14     int nfiles, append, batch, icurrent, ibatno;
15     } files;
16     EXTERN struct ElementType { char symbol[3];
17     int atomnum;
18     float weight, covradius, vdwradius;
19     int s,p,d,f,type ;
20     } Elements[];
21     EXTERN struct t_dipolemom {
22     double total, xdipole, ydipole, zdipole;
23     } dipolemom;
24     EXTERN struct t_logp {
25     float logp;
26     } logp_calc;
27     EXTERN struct t_vibdata {
28     char ptgrp[4];
29     float mom_ix,mom_iy,mom_iz;
30     float etot,htot,stot,gtot,cptot;
31     } vibdata;
32    
33 wdelano 58 // ============================
34    
35 tjod 3 int make_atom(int, float, float, float,char *);
36     void make_bond(int, int, int);
37     int read_sdf(int,int);
38     int rd_sdf(FILE *);
39     void write_sdf(int);
40     void message_alert(char *, char *);
41     void hdel(int);
42     FILE * fopen_path ( char * , char * , char * ) ;
43     int FetchRecord(FILE *, char *);
44     void avgleg(void);
45 wdelano 58 void get_tag(char *,char *);
46     // ==================================
47 tjod 3 /*
48     * flags - indicates whether dipole, xlogp or vibrational calculations were done
49     * and if so write them to the output file. Look at the definitions in pcmod.h
50     *
51     */
52     void write_sdf(int flags)
53     {
54     int i,j,lptest,nbond,junk;
55 wdelano 58 FILE *wfile = pcmoutfile;
56 tjod 3
57     lptest = 0;
58     for( i = 1; i <= natom; i++ )
59     {
60     if( atom[i].mmx_type == 20 )
61     {
62     lptest = 1;
63     hdel( lptest );
64     break;
65     }
66     }
67     nbond = 0;
68     /* ** calculate the number of bonds in the molecule ** */
69     for( j = 1; j <= natom; j++ )
70     {
71     for( i = 0; i < MAXIAT; i++ )
72     {
73     if( atom[j].iat[i] != 0 )
74     {
75     if( atom[j].iat[i] < j )
76     nbond = nbond + 1;
77     }
78     }
79     }
80     /* now write the concord file */
81 tjod 15 /*
82 tjod 3 if (files.append)
83     wfile = fopen_path(Savebox.path,Savebox.fname,"a");
84     else
85     wfile = fopen_path(Savebox.path,Savebox.fname,"w");
86 tjod 15 */
87 tjod 3
88     j = strlen(Struct_Title);
89     for (i=0; i < j; i++)
90     {
91     if (Struct_Title[i] == '\n')
92     {
93     Struct_Title[i] = '\0';
94     break;
95     }
96     }
97     fprintf(wfile,"%s\n",Struct_Title);
98     fprintf(wfile," PCMODEL v9.1 1.00000 0.00000\n");
99     fprintf(wfile,"\n");
100    
101     fprintf(wfile,"%3d%3d 0 0 0 0 1 V2000\n",natom,nbond);
102    
103     for (i=1; i <= natom; i++)
104     {
105     junk = 0;
106     if (atom[i].mmx_type == 41) junk = 3;
107     else if (atom[i].mmx_type == 42) junk = 5;
108     else if (atom[i].mmx_type == 66)
109     {
110     if (atom[i].bo[0] == 1)
111     junk = 5;
112     }
113     fprintf(wfile,"%10.4f%10.4f%10.4f %-3s 0 %d 0 0 0 0 \n",atom[i].x, atom[i].y,
114     atom[i].z,Elements[atom[i].atomnum-1].symbol,junk);
115     }
116     for (i=1; i <= natom; i++)
117     {
118     for (j=0; j < MAXIAT; j++)
119     {
120     if (atom[i].iat[j] != 0 && i < atom[i].iat[j])
121     fprintf(wfile,"%3d%3d%3d 0\n",i, atom[i].iat[j], atom[i].bo[j]);
122     }
123     }
124     fprintf(wfile,"M END\n");
125     fprintf(wfile,"> <title>\n");
126     fprintf(wfile,"%s\n",Struct_Title);
127    
128     fprintf(wfile,"> <MMFF94 energy>\n");
129     fprintf(wfile,"%f\n",energies.total);
130    
131     if (flags & DO_DIPOLE) {
132     fprintf(wfile,"> <dipole moment>\n");
133     fprintf(wfile,"%f\n",dipolemom.total);
134     }
135    
136     if (flags & DO_XLOGP) {
137     fprintf(wfile,"> <xLogP>\n");
138     fprintf(wfile,"%f\n",logp_calc.logp);
139     }
140    
141     if (flags & DO_VIBRATION) {
142     fprintf(wfile,"> <Point Group>\n");
143     fprintf(wfile,"%s\n",vibdata.ptgrp);
144    
145     fprintf(wfile,"> <Moments of Inertia>\n");
146     fprintf(wfile,"%f8.3 %f8.3 %f8.3 \n",vibdata.mom_ix,vibdata.mom_iy,vibdata.mom_iz);
147    
148     fprintf(wfile,"> <thermodynamics dE, dH, S, dG, CP>\n");
149     fprintf(wfile,"%f8.3 %f8.3 %f8.3 %f8.3 %f8.3 \n",vibdata.etot, vibdata.htot,
150     vibdata.stot,vibdata.gtot,vibdata.cptot);
151     }
152    
153     fprintf(wfile,"\n");
154     fprintf(wfile,"$$$$\n");
155 tjod 15 // fclose(wfile);
156 tjod 3 }
157 wdelano 58
158 tjod 3 /* =================================== */
159     // fast read - assume file is open and positioned
160     // read structure and down to end $$$$
161     // and return
162     int rd_sdf(FILE *rfile)
163     {
164 wdelano 58 int i, j, niatom, nibond, ia1, ia2, ia3, ia4, ibond, itype, newatom;
165     int ncount,junk,junk1,got_title,istereo,iz,nvalue;
166 tjod 3 int jji, jj1, jj2, jj3, jj4;
167 wdelano 58 int n_row;
168 tjod 3 int has_Aromatic, xPlus,yPlus,zPlus, planar;
169     int icount,itemp[4];
170     int Ret_Val;
171     float xtmp, ytmp, ztmp, dz;
172 wdelano 58 float fconst,min,max;
173 tjod 3 char c1[4],c2[4];
174     char atomchar[3];
175     char inputline[150];
176 wdelano 58 char tag[30];
177     char ***tab;
178 tjod 3
179     Ret_Val = TRUE;
180     got_title = FALSE;
181     xPlus = yPlus = zPlus = FALSE;
182 tjod 15 if ( 0 == FetchRecord(rfile,inputline) )return -1;
183     // sscanf(inputline,"SDF %s",Struct_Title);
184     strncpy(Struct_Title, inputline, sizeof(Struct_Title));
185 tjod 3 got_title = TRUE;
186     /* if (strlen(inputline) > 4)
187     {
188     iz = strlen(inputline);
189     if (iz > 60) iz = 59;
190     for (i=4; i < iz; i++)
191     Struct_Title[i-4] = inputline[i];
192     Struct_Title[i] = '\0';
193     got_title = TRUE;
194     } */
195     FetchRecord(rfile,inputline); // blank line
196     FetchRecord(rfile,inputline); // blank line
197    
198     FetchRecord(rfile,inputline); // natom and nbond
199     for (i=0; i <4; i++)
200     {
201     c1[i] = inputline[i];
202     c2[i] = inputline[i+3];
203     }
204     c1[3] = '\0';c2[3] = '\0';
205     niatom = atoi(c1);
206     nibond = atoi(c2);
207     if (niatom == 0)
208     return FALSE;
209    
210     for (i=0; i < niatom; i++)
211     {
212     FetchRecord(rfile,inputline);
213     sscanf(inputline,"%f %f %f %s %d %d",&xtmp, &ytmp, &ztmp, atomchar,&junk,&junk1);
214    
215     itype = 0;
216     if (xtmp != 0.0) xPlus= TRUE;
217     if (ytmp != 0.0) yPlus= TRUE;
218     if (ztmp != 0.0) zPlus= TRUE;
219    
220     iz = strlen(atomchar);
221     if (itype == 0 && (atomchar[0] == 'N' || atomchar[0] == 'O') && iz == 1) // nitro fix
222     {
223     if (atomchar[0] == 'N') itype = 8;
224     if (atomchar[0] == 'O') itype = 6;
225    
226     if (junk1 != 0)
227     {
228     if (junk1 == 3 && itype == 8)
229     itype = 41; // N+
230     if (junk1 == 5 && itype == 6)
231     itype = 42; // O-
232     }
233     }
234    
235     newatom = make_atom(itype,xtmp,ytmp,ztmp,atomchar);
236     if (itype != 0)
237     atom[newatom].mmx_type = itype;
238     if (newatom == -1)
239     {
240 wdelano 58 fprintf(pcmlogfile, "Atom %d %s not recognized structure skipped\n",i,atomchar);
241 tjod 3 Ret_Val = FALSE;
242     }
243     }
244     has_Aromatic = FALSE;
245    
246     planar = TRUE;
247     if (xPlus && yPlus && zPlus) planar = FALSE;
248    
249     for (i=0; i < nibond; i++)
250     {
251     FetchRecord(rfile,inputline);
252     for (j=0; j <4; j++)
253     {
254     c1[j] = inputline[j];
255     c2[j] = inputline[j+3];
256     }
257     c1[3] = '\0';c2[3] = '\0';
258     ia1 = atoi(c1);
259     ia2 = atoi(c2);
260     istereo = 0;
261     sscanf(inputline,"%3d%3d%3d%3d",&junk, &junk1, &ibond,&istereo);
262     if (ibond >= 4)
263     {
264     ibond = 1;
265     has_Aromatic = TRUE;
266 tjod 31 // TJO
267     Ret_Val = FALSE;
268 tjod 3 }
269     make_bond(ia1, ia2, ibond);
270     if (istereo == 1 && planar)
271     {
272     atom[ia2].z += 0.7;
273     if (atom[ia2].atomnum == 1)
274     atom[ia2].type = 60;
275     }
276     if (istereo == 6 && planar)
277     {
278     atom[ia2].z -= 0.7;
279     if (atom[ia2].atomnum == 1)
280     atom[ia2].type = 60;
281     }
282     }
283     // read to end of structure
284 wdelano 58 // parse any tags here
285     // agreed tags
286     // <FIXED_ATOMS>
287     // <RESTRAINED_ATOMS>
288     // <RESTRAINED_DISTANCES>
289     // <RESTRAINED_ANGLES>
290     // <RESTRAINED_DIHEDRALS>
291    
292 tjod 3 while (FetchRecord(rfile,inputline))
293     {
294     if (strncasecmp(inputline,"$$$$",4) == 0)
295     {
296     return Ret_Val;
297 wdelano 58 } else if (inputline[0] == '>')
298     {
299     get_tag(inputline,tag);
300     if (strcasecmp(tag,"fixed_atoms") == 0)
301     {
302     tab = tabulator_new_from_file_using_header(rfile,"ATOM",0);
303     if (tab)
304     {
305     fprintf(pcmlogfile,"got table\n");
306     n_row = tabulator_height(tab);
307     for (i=0; i < n_row; i++)
308     {
309     ia1 = atoi(tab[i][0]);
310     if (ia1 > 0 && ia1 <= natom)
311     {
312     fx_atom.katom_fix[fx_atom.natom_fix] = ia1;
313     fx_atom.natom_fix++;
314     }
315     }
316     }
317     } else if (strcasecmp(tag,"restrained_atoms") == 0)
318     {
319     tab = tabulator_new_from_file_using_header(rfile,"ATOM1 MAX F_CONST X Y Z",0);
320     if (tab)
321     {
322     n_row = tabulator_height(tab);
323     for (i=0; i < n_row; i++)
324     {
325     ia1 = atoi(tab[i][0]);
326     max = atof(tab[i][1]);
327     fconst = atof(tab[i][2]);
328     if (ia1 > 0 && ia1 <= natom)
329     {
330     restrain_atom.katom_restrain[restrain_atom.natom_restrain] = ia1;
331     restrain_atom.restrain_const[restrain_atom.natom_restrain] = fconst;
332     restrain_atom.restrain_max[restrain_atom.natom_restrain] = max;
333     restrain_atom.restrain_position[restrain_atom.natom_restrain][0] = atom[ia1].x; // default positions
334     restrain_atom.restrain_position[restrain_atom.natom_restrain][1] = atom[ia1].y;
335     restrain_atom.restrain_position[restrain_atom.natom_restrain][2] = atom[ia1].z;
336     if (strlen(tab[i][3]) > 0) // x postion
337     {
338     xtmp = atof(tab[i][3]);
339     restrain_atom.restrain_position[restrain_atom.natom_restrain][0] = xtmp;
340     }
341     if (strlen(tab[i][4]) > 0) // y postion
342     {
343     xtmp = atof(tab[i][4]);
344     restrain_atom.restrain_position[restrain_atom.natom_restrain][1] = xtmp;
345     }
346     if (strlen(tab[i][5]) > 0) // y postion
347     {
348     xtmp = atof(tab[i][3]);
349     restrain_atom.restrain_position[restrain_atom.natom_restrain][2] = xtmp;
350     }
351     restrain_atom.natom_restrain++;
352     }
353     }
354     }
355     } else if (strcasecmp(tag,"restrained_distances") == 0)
356     {
357     tab = tabulator_new_from_file_using_header(rfile,"ATOM1 ATOM2 MIN MAX F_CONST",0);
358     if (tab)
359     {
360     n_row = tabulator_height(tab);
361     for (i=0; i < n_row; i++)
362     {
363     ia1 = atoi(tab[i][0]);
364     ia2 = atoi(tab[i][1]);
365     min = atof(tab[i][2]);
366     max = atof(tab[i][3]);
367     fconst = atof(tab[i][4]);
368     if ((ia1 > 0 && ia1 <= natom) && (ia2 > 0 && ia2 < natom))
369     {
370     fx_dist.kdfix[fx_dist.ndfix][0] = ia1;
371     fx_dist.kdfix[fx_dist.ndfix][1] = ia2;
372     fx_dist.fdconst[fx_dist.ndfix] = fconst;
373     fx_dist.min_dist[fx_dist.ndfix] = min;
374     fx_dist.max_dist[fx_dist.ndfix] = max;
375     fx_dist.ndfix++;
376     }
377     }
378     }
379     } else if (strcasecmp(tag,"restrained_angles") == 0)
380     {
381     tab = tabulator_new_from_file_using_header(rfile,"ATOM1 ATOM2 ATOM3 MIN MAX F_CONST",0);
382     if (tab)
383     {
384     n_row = tabulator_height(tab);
385     for (i=0; i < n_row; i++)
386     {
387     ia1 = atoi(tab[i][0]);
388     ia2 = atoi(tab[i][1]);
389     ia3 = atoi(tab[i][2]);
390     min = atof(tab[i][3]);
391     max = atof(tab[i][4]);
392     fconst = atof(tab[i][5]);
393     if ((ia1 > 0 && ia1 <= natom) && (ia2 > 0 && ia2 < natom) && (ia3 > 0 && ia3 <= natom))
394     {
395     fx_angle.kafix[fx_angle.nafix][0] = ia1;
396     fx_angle.kafix[fx_angle.nafix][1] = ia2;
397     fx_angle.kafix[fx_angle.nafix][2] = ia3;
398     fx_angle.faconst[fx_angle.nafix] = fconst;
399     fx_angle.min_ang[fx_angle.nafix] = min;
400     fx_angle.max_ang[fx_angle.nafix] = max;
401     fx_angle.nafix++;
402     }
403     }
404     }
405     } else if (strcasecmp(tag,"restrained_dihedrals") == 0)
406     {
407     tab = tabulator_new_from_file_using_header(rfile,"ATOM1 ATOM2 ATOM3 ATOM4 MIN MAX F_CONST",0);
408     if (tab)
409     {
410     n_row = tabulator_height(tab);
411     for (i=0; i < n_row; i++)
412     {
413     ia1 = atoi(tab[i][0]);
414     ia2 = atoi(tab[i][1]);
415     ia3 = atoi(tab[i][2]);
416     ia4 = atoi(tab[i][3]);
417     min = atof(tab[i][4]);
418     max = atof(tab[i][5]);
419     fconst = atof(tab[i][6]);
420     if ((ia1 > 0 && ia1 <= natom) && (ia2 > 0 && ia2 < natom) && (ia3 > 0 && ia3 <= natom) && (ia4 > 0 && ia4 <= natom))
421     {
422     fx_torsion.ktfix[fx_torsion.ntfix][0] = ia1;
423     fx_torsion.ktfix[fx_torsion.ntfix][1] = ia2;
424     fx_torsion.ktfix[fx_torsion.ntfix][2] = ia3;
425     fx_torsion.ktfix[fx_torsion.ntfix][3] = ia4;
426     fx_torsion.ftconst[fx_torsion.ntfix] = fconst;
427     fx_torsion.min_tor[fx_torsion.ntfix] = min;
428     fx_torsion.max_tor[fx_torsion.ntfix] = max;
429     fx_torsion.ntfix++;
430     }
431     }
432     }
433     }
434 tjod 3 }
435     }
436     return Ret_Val;
437     }
438 wdelano 58 // =====================
439     void get_tag(char *line,char *tag)
440 tjod 3 {
441 wdelano 58 int i,j,icount;
442     icount = 0;
443     strcpy(tag,"");
444     for (i=0; i < strlen(line); i++)
445 tjod 3 {
446 wdelano 58 if (line[i] == '<')
447     {
448     for (j=i+1; j < strlen(line); j++)
449     {
450     if (line[j] == '>')
451     {
452     tag[icount] = '\0';
453     return;
454     } else
455     {
456     tag[icount] = line[j];
457     icount++;
458     }
459     }
460     }
461 tjod 3 }
462     }