ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/freemol/trunk/src/mengine/src/read_sdf.c
Revision: 89
Committed: Thu Jan 15 14:13:58 2009 UTC (13 years, 4 months ago) by gilbertke
File size: 18635 byte(s)
Log Message:
resolved problem in read_sdf, removed omp pragma in ehal
Line User Rev File contents
1 gilbertke 67 /* NOTICE: this source code file has been modified for use with FreeMOL */
2     #define EXTERN extern
3    
4     #include "pcwin.h"
5     #include "pcmod.h"
6     #include "atom_k.h"
7     #include "energies.h"
8     #include "fix.h"
9     #include "dipmom.h"
10     #include "draw.h"
11     #include "vibrate.h"
12    
13     #define TABULATOR_INCLUDE_IMPLEMENTATION
14     #include "tabulator.h"
15    
16     EXTERN struct t_files {
17     int nfiles, append, batch, icurrent, ibatno;
18     } files;
19     EXTERN struct ElementType { char symbol[3];
20     int atomnum;
21     float weight, covradius, vdwradius;
22     int s,p,d,f,type ;
23     } Elements[];
24     EXTERN struct t_logp {
25     float logp;
26     } logp_calc;
27    
28     // ============================
29     void quick_type(void);
30     int read_sdf(int,int);
31     int rd_sdf(FILE *);
32     void write_sdf(int);
33     void hdel(int);
34     FILE * fopen_path ( char * , char * , char * ) ;
35     int FetchRecord(FILE *, char *);
36     void avgleg(void);
37     void get_tag(char *,char *);
38     static void mopaco(int start_atom,int end_atom);
39     // ==================================
40     /*
41     * flags - indicates whether dipole, xlogp or vibrational calculations were done
42     * and if so write them to the output file. Look at the definitions in pcmod.h
43     *
44     */
45     void write_sdf(int flags)
46     {
47     int i,j,lptest,nbond,junk;
48     FILE *wfile = pcmoutfile;
49    
50     lptest = 0;
51     for( i = 1; i <= natom; i++ )
52     {
53     if( atom[i].mmx_type == 20 )
54     {
55     lptest = 1;
56     hdel( lptest );
57     break;
58     }
59     }
60     nbond = 0;
61     /* ** calculate the number of bonds in the molecule ** */
62     for( j = 1; j <= natom; j++ )
63     {
64     for( i = 0; i < MAXIAT; i++ )
65     {
66     if( atom[j].iat[i] != 0 )
67     {
68     if( atom[j].iat[i] < j )
69     nbond = nbond + 1;
70     }
71     }
72     }
73     /* now write the concord file */
74     /*
75     if (files.append)
76     wfile = fopen_path(Savebox.path,Savebox.fname,"a");
77     else
78     wfile = fopen_path(Savebox.path,Savebox.fname,"w");
79     */
80    
81     j = strlen(Struct_Title);
82     for (i=0; i < j; i++)
83     {
84     if (Struct_Title[i] == '\n')
85     {
86     Struct_Title[i] = '\0';
87     break;
88     }
89     }
90     fprintf(wfile,"%s\n",Struct_Title);
91     fprintf(wfile," MENGINE v1.0 1.00000 0.00000\n");
92     fprintf(wfile,"\n");
93    
94     fprintf(wfile,"%3d%3d 0 0 0 0 1 V2000\n",natom,nbond);
95    
96     for (i=1; i <= natom; i++)
97     {
98     junk = 0;
99     if (atom[i].mmx_type == 41) junk = 3;
100     else if (atom[i].mmx_type == 42) junk = 5;
101     else if (atom[i].mmx_type == 66)
102     {
103     if (atom[i].bo[0] == 1)
104     junk = 5;
105     }
106     fprintf(wfile,"%10.4f%10.4f%10.4f %-3s 0 %d 0 0 0 0 \n",atom[i].x, atom[i].y,
107     atom[i].z,Elements[atom[i].atomnum-1].symbol,junk);
108     }
109     for (i=1; i <= natom; i++)
110     {
111     for (j=0; j < MAXIAT; j++)
112     {
113     if (atom[i].iat[j] != 0 && i < atom[i].iat[j])
114     fprintf(wfile,"%3d%3d%3d 0\n",i, atom[i].iat[j], atom[i].bo[j]);
115     }
116     }
117     fprintf(wfile,"M END\n");
118     fprintf(wfile,"> <title>\n");
119     fprintf(wfile,"%s\n",Struct_Title);
120    
121     fprintf(wfile,"> <MMFF94 energy>\n");
122     fprintf(wfile,"%f\n",energies.total);
123    
124     fprintf(wfile,"> <MMFF94_Atomtypes>\n");
125     fprintf(wfile,"+ Atom Atomtype\n");
126     for (i=1; i <=natom; i++)
127     fprintf(wfile,"| %d %d\n",i,atom[i].type);
128     fprintf(wfile,"\n");
129    
130     if (flags & DO_DIPOLE) {
131     fprintf(wfile,"> <dipole moment>\n");
132     fprintf(wfile,"%f\n",dipolemom.total);
133     }
134    
135     if (flags & DO_XLOGP) {
136     fprintf(wfile,"> <xLogP>\n");
137     fprintf(wfile,"%f\n",logp_calc.logp);
138     }
139    
140     if (flags & DO_VIBRATION) {
141     fprintf(wfile,"> <Point Group>\n");
142     fprintf(wfile,"%s\n",vibdata.ptgrp);
143    
144     fprintf(wfile,"> <Moments of Inertia>\n");
145     fprintf(wfile,"%f8.3 %f8.3 %f8.3 \n",vibdata.mom_ix,vibdata.mom_iy,vibdata.mom_iz);
146    
147     fprintf(wfile,"> <thermodynamics dE, dH, S, dG, CP>\n");
148     fprintf(wfile,"%f8.3 %f8.3 %f8.3 %f8.3 %f8.3 \n",vibdata.etot, vibdata.htot,
149     vibdata.stot,vibdata.gtot,vibdata.cptot);
150     }
151    
152     fprintf(wfile,"\n");
153     fprintf(wfile,"$$$$\n");
154     // fclose(wfile);
155     }
156    
157     /* =================================== */
158     // fast read - assume file is open and positioned
159     // read structure and down to end $$$$
160     // and return
161     int rd_sdf(FILE *rfile)
162     {
163     int i, j, niatom, nibond, ia1, ia2, ia3, ia4, ibond, itype, newatom;
164     int ncount,junk,junk1,got_title,istereo,iz,nvalue;
165     int jji, jj1, jj2, jj3, jj4;
166     int n_row;
167     int has_Aromatic, xPlus,yPlus,zPlus, planar;
168     int icount,itemp[4];
169     int Ret_Val;
170     float xtmp, ytmp, ztmp, dz;
171     float fconst,min,max;
172     char c1[4],c2[4];
173     char atomchar[3];
174     char inputline[150];
175     char tag[30];
176     char ***tab;
177    
178     Ret_Val = TRUE;
179     got_title = FALSE;
180     xPlus = yPlus = zPlus = FALSE;
181     if ( 0 == FetchRecord(rfile,inputline) )return -1;
182     // sscanf(inputline,"SDF %s",Struct_Title);
183     strncpy(Struct_Title, inputline, sizeof(Struct_Title));
184     got_title = TRUE;
185     /* if (strlen(inputline) > 4)
186     {
187     iz = strlen(inputline);
188     if (iz > 60) iz = 59;
189     for (i=4; i < iz; i++)
190     Struct_Title[i-4] = inputline[i];
191     Struct_Title[i] = '\0';
192     got_title = TRUE;
193     } */
194     FetchRecord(rfile,inputline); // blank line
195     FetchRecord(rfile,inputline); // blank line
196    
197     FetchRecord(rfile,inputline); // natom and nbond
198     for (i=0; i <4; i++)
199     {
200     c1[i] = inputline[i];
201     c2[i] = inputline[i+3];
202     }
203     c1[3] = '\0';c2[3] = '\0';
204     niatom = atoi(c1);
205     nibond = atoi(c2);
206     if (niatom == 0)
207     return FALSE;
208    
209     for (i=0; i < niatom; i++)
210     {
211     FetchRecord(rfile,inputline);
212     sscanf(inputline,"%f %f %f %s %d %d",&xtmp, &ytmp, &ztmp, atomchar,&junk,&junk1);
213    
214     itype = 0;
215     if (xtmp != 0.0) xPlus= TRUE;
216     if (ytmp != 0.0) yPlus= TRUE;
217     if (ztmp != 0.0) zPlus= TRUE;
218    
219     iz = strlen(atomchar);
220     if (itype == 0 && (atomchar[0] == 'N' || atomchar[0] == 'O') && iz == 1) // nitro fix
221     {
222     if (atomchar[0] == 'N') itype = 8;
223     if (atomchar[0] == 'O') itype = 6;
224    
225     if (junk1 != 0)
226     {
227     if (junk1 == 3 && itype == 8)
228     itype = 41; // N+
229     if (junk1 == 5 && itype == 6)
230     itype = 42; // O-
231     }
232     }
233     newatom = make_atom(itype,xtmp,ytmp,ztmp,atomchar);
234     if (itype != 0)
235     atom[newatom].mmx_type = itype;
236     if (newatom == -1)
237     {
238     fprintf(pcmlogfile, "Atom %d %s not recognized structure skipped\n",i,atomchar);
239     Ret_Val = FALSE;
240     }
241     }
242     has_Aromatic = FALSE;
243    
244     planar = TRUE;
245     if (xPlus && yPlus && zPlus) planar = FALSE;
246    
247     for (i=0; i < nibond; i++)
248     {
249     FetchRecord(rfile,inputline);
250     for (j=0; j <4; j++)
251     {
252     c1[j] = inputline[j];
253     c2[j] = inputline[j+3];
254     }
255     c1[3] = '\0';c2[3] = '\0';
256     ia1 = atoi(c1);
257     ia2 = atoi(c2);
258     istereo = 0;
259     sscanf(inputline,"%3d%3d%3d%3d",&junk, &junk1, &ibond,&istereo);
260     if (ibond >= 4)
261     {
262     // use bonds array as temporary storage for aromatic bonds
263     make_bond(ia1,ia2,1);
264     bonds.ia1[bonds.numbonds] = ia1;
265     bonds.ia2[bonds.numbonds] = ia2;
266     bonds.numbonds++;
267     has_Aromatic = TRUE;
268     } else
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     // parse any tags here
285     // agreed tags
286     // <FIXED_ATOMS>
287     // <RESTRAINED_ATOMS>
288     // <RESTRAINED_DISTANCES>
289     // <RESTRAINED_ANGLES>
290     // <RESTRAINED_DIHEDRALS>
291    
292     while (FetchRecord(rfile,inputline))
293     {
294     if (strncasecmp(inputline,"$$$$",4) == 0)
295     {
296     goto L_DONE;
297     } 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,"ATOM 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][5]);
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     }
435     }
436     // if (has_aromatic)
437     // need to deal with aromatic bonds
438     L_DONE:
439     quick_type();
440     if (has_Aromatic)
441     {
442     for (i=0; i < bonds.numbonds; i++)
443     mopaco(bonds.ia1[i],bonds.ia2[i]);
444     }
445     return Ret_Val;
446     }
447     // =====================
448     void get_tag(char *line,char *tag)
449     {
450     int i,j,icount;
451     icount = 0;
452     strcpy(tag,"");
453     for (i=0; i < strlen(line); i++)
454     {
455     if (line[i] == '<')
456     {
457     for (j=i+1; j < strlen(line); j++)
458     {
459     if (line[j] == '>')
460     {
461     tag[icount] = '\0';
462     return;
463     } else
464     {
465     tag[icount] = line[j];
466     icount++;
467     }
468     }
469     }
470     }
471     }
472     /* ------------------------------------*/
473     // new mopaco version - test valency of each atom and add bond if possible
474     void mopaco(int ia, int ib)
475     {
476     int ibi, ibj, it, jt, kk;
477    
478     it = atom[ia].atomnum;
479     jt = atom[ib].atomnum;
480    
481     if( it == 6 || it == 7 || it == 8 || it == 15 || it == 16 )
482     {
483     if( jt == 6 || jt == 7 || jt == 8 || jt == 15 || jt == 16 )
484     {
485     /* c=n */
486     if( (it == 7 && jt == 6) || (it == 6 && jt == 7) )
487     {
488     ibi = 0;
489     for( kk = 0; kk < MAXIAT; kk++ )
490     {
491     if( atom[ia].bo[kk] != 9 )
492     ibi = ibi + atom[ia].bo[kk];
493     }
494     ibj = 0;
495     for( kk = 0; kk < MAXIAT; kk++ )
496     {
497     if( atom[ib].bo[kk] != 9 )
498     ibj = ibj + atom[ib].bo[kk];
499     }
500     if( it == 7 )
501     {
502     if( ibi <= 2 && ibj <= 3 )
503     {
504     make_bond( ia, ib, 1 );
505     }
506     }else if( jt == 7 )
507     {
508     if( ibj <= 2 && ibi <= 3 )
509     {
510     make_bond( ia, ib, 1 );
511     }
512     }
513     /* c=c bond */
514     } else if( (it == 6 && jt == 6) )
515     {
516     ibi = 0;
517     for( kk = 0; kk < MAXIAT; kk++ )
518     {
519     if( atom[ia].bo[kk] != 9 )
520     ibi = ibi + atom[ia].bo[kk];
521     }
522     ibj = 0;
523     for( kk = 0; kk < MAXIAT; kk++ )
524     {
525     if( atom[ib].bo[kk] != 9 )
526     ibj = ibj + atom[ib].bo[kk];
527     }
528     if( ibi <= 3 && ibj <= 3 )
529     {
530     make_bond( ia, ib, 1 );
531     }
532     }
533     }
534     }
535     }
536