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