ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/freemol/branches/sync4pymol12/src/mengine/src/read_sdf.c
Revision: 99
Committed: Mon Jan 19 04:31:29 2009 UTC (10 years, 7 months ago) by wdelano
File size: 18838 byte(s)
Log Message:
synchronized with trunk, less openmp lib/include
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 wdelano 99 #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     // ==================================
41     /*
42     * flags - indicates whether dipole, xlogp or vibrational calculations were done
43     * and if so write them to the output file. Look at the definitions in pcmod.h
44     *
45     */
46     void write_sdf(int flags)
47     {
48     int i,j,lptest,nbond,junk;
49     FILE *wfile = pcmoutfile;
50    
51     lptest = 0;
52     for( i = 1; i <= natom; i++ )
53     {
54     if( atom[i].mmx_type == 20 )
55     {
56     lptest = 1;
57     hdel( lptest );
58     break;
59     }
60     }
61     nbond = 0;
62     /* ** calculate the number of bonds in the molecule ** */
63     for( j = 1; j <= natom; j++ )
64     {
65     for( i = 0; i < MAXIAT; i++ )
66     {
67     if( atom[j].iat[i] != 0 )
68     {
69     if( atom[j].iat[i] < j )
70     nbond = nbond + 1;
71     }
72     }
73     }
74     /* now write the concord file */
75     /*
76     if (files.append)
77     wfile = fopen_path(Savebox.path,Savebox.fname,"a");
78     else
79     wfile = fopen_path(Savebox.path,Savebox.fname,"w");
80     */
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," MENGINE v1.0 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     fprintf(wfile,"> <MMFF94_Atomtypes>\n");
126     fprintf(wfile,"+ Atom Atomtype\n");
127     for (i=1; i <=natom; i++)
128     fprintf(wfile,"| %d %d\n",i,atom[i].type);
129     fprintf(wfile,"\n");
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     // fclose(wfile);
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 wdelano 99 int numbonds, *ib1, *ib2;
171 gilbertke 67 float xtmp, ytmp, ztmp, dz;
172     float fconst,min,max;
173     char c1[4],c2[4];
174     char atomchar[3];
175     char inputline[150];
176     char tag[30];
177     char ***tab;
178    
179     Ret_Val = TRUE;
180     got_title = FALSE;
181     xPlus = yPlus = zPlus = FALSE;
182     if ( 0 == FetchRecord(rfile,inputline) )return -1;
183     // sscanf(inputline,"SDF %s",Struct_Title);
184     strncpy(Struct_Title, inputline, sizeof(Struct_Title));
185     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 wdelano 99 // allocate space for aromatic bonds
210     numbonds = 0;
211     ib1 = ivector(0,nibond);
212     ib2 = ivector(0,nibond);
213     //
214 gilbertke 67 for (i=0; i < niatom; i++)
215     {
216     FetchRecord(rfile,inputline);
217     sscanf(inputline,"%f %f %f %s %d %d",&xtmp, &ytmp, &ztmp, atomchar,&junk,&junk1);
218    
219     itype = 0;
220     if (xtmp != 0.0) xPlus= TRUE;
221     if (ytmp != 0.0) yPlus= TRUE;
222     if (ztmp != 0.0) zPlus= TRUE;
223    
224     iz = strlen(atomchar);
225     if (itype == 0 && (atomchar[0] == 'N' || atomchar[0] == 'O') && iz == 1) // nitro fix
226     {
227     if (atomchar[0] == 'N') itype = 8;
228     if (atomchar[0] == 'O') itype = 6;
229    
230     if (junk1 != 0)
231     {
232     if (junk1 == 3 && itype == 8)
233     itype = 41; // N+
234     if (junk1 == 5 && itype == 6)
235     itype = 42; // O-
236     }
237     }
238     newatom = make_atom(itype,xtmp,ytmp,ztmp,atomchar);
239     if (itype != 0)
240     atom[newatom].mmx_type = itype;
241     if (newatom == -1)
242     {
243     fprintf(pcmlogfile, "Atom %d %s not recognized structure skipped\n",i,atomchar);
244     Ret_Val = FALSE;
245     }
246     }
247     has_Aromatic = FALSE;
248    
249     planar = TRUE;
250     if (xPlus && yPlus && zPlus) planar = FALSE;
251    
252     for (i=0; i < nibond; i++)
253     {
254     FetchRecord(rfile,inputline);
255     for (j=0; j <4; j++)
256     {
257     c1[j] = inputline[j];
258     c2[j] = inputline[j+3];
259     }
260     c1[3] = '\0';c2[3] = '\0';
261     ia1 = atoi(c1);
262     ia2 = atoi(c2);
263     istereo = 0;
264     sscanf(inputline,"%3d%3d%3d%3d",&junk, &junk1, &ibond,&istereo);
265     if (ibond >= 4)
266     {
267     // use bonds array as temporary storage for aromatic bonds
268     make_bond(ia1,ia2,1);
269 wdelano 99 ib1[numbonds] = ia1;
270     ib2[numbonds] = ia2;
271     numbonds++;
272 gilbertke 67 has_Aromatic = TRUE;
273     } else
274     make_bond(ia1, ia2, ibond);
275     if (istereo == 1 && planar)
276     {
277     atom[ia2].z += 0.7;
278     if (atom[ia2].atomnum == 1)
279     atom[ia2].type = 60;
280     }
281     if (istereo == 6 && planar)
282     {
283     atom[ia2].z -= 0.7;
284     if (atom[ia2].atomnum == 1)
285     atom[ia2].type = 60;
286     }
287     }
288     // read to end of structure
289     // parse any tags here
290     // agreed tags
291     // <FIXED_ATOMS>
292     // <RESTRAINED_ATOMS>
293     // <RESTRAINED_DISTANCES>
294     // <RESTRAINED_ANGLES>
295     // <RESTRAINED_DIHEDRALS>
296    
297     while (FetchRecord(rfile,inputline))
298     {
299     if (strncasecmp(inputline,"$$$$",4) == 0)
300     {
301     goto L_DONE;
302     } else if (inputline[0] == '>')
303     {
304     get_tag(inputline,tag);
305     if (strcasecmp(tag,"fixed_atoms") == 0)
306     {
307     tab = tabulator_new_from_file_using_header(rfile,"ATOM",0);
308     if (tab)
309     {
310     fprintf(pcmlogfile,"got table\n");
311     n_row = tabulator_height(tab);
312     for (i=0; i < n_row; i++)
313     {
314     ia1 = atoi(tab[i][0]);
315     if (ia1 > 0 && ia1 <= natom)
316     {
317     fx_atom.katom_fix[fx_atom.natom_fix] = ia1;
318     fx_atom.natom_fix++;
319     }
320     }
321     }
322     } else if (strcasecmp(tag,"restrained_atoms") == 0)
323     {
324     tab = tabulator_new_from_file_using_header(rfile,"ATOM MAX F_CONST X Y Z",0);
325     if (tab)
326     {
327     n_row = tabulator_height(tab);
328     for (i=0; i < n_row; i++)
329     {
330     ia1 = atoi(tab[i][0]);
331     max = atof(tab[i][1]);
332     fconst = atof(tab[i][2]);
333     if (ia1 > 0 && ia1 <= natom)
334     {
335     restrain_atom.katom_restrain[restrain_atom.natom_restrain] = ia1;
336     restrain_atom.restrain_const[restrain_atom.natom_restrain] = fconst;
337     restrain_atom.restrain_max[restrain_atom.natom_restrain] = max;
338     restrain_atom.restrain_position[restrain_atom.natom_restrain][0] = atom[ia1].x; // default positions
339     restrain_atom.restrain_position[restrain_atom.natom_restrain][1] = atom[ia1].y;
340     restrain_atom.restrain_position[restrain_atom.natom_restrain][2] = atom[ia1].z;
341     if (strlen(tab[i][3]) > 0) // x postion
342     {
343     xtmp = atof(tab[i][3]);
344     restrain_atom.restrain_position[restrain_atom.natom_restrain][0] = xtmp;
345     }
346     if (strlen(tab[i][4]) > 0) // y postion
347     {
348     xtmp = atof(tab[i][4]);
349     restrain_atom.restrain_position[restrain_atom.natom_restrain][1] = xtmp;
350     }
351     if (strlen(tab[i][5]) > 0) // y postion
352     {
353     xtmp = atof(tab[i][5]);
354     restrain_atom.restrain_position[restrain_atom.natom_restrain][2] = xtmp;
355     }
356     restrain_atom.natom_restrain++;
357     }
358     }
359     }
360     } else if (strcasecmp(tag,"restrained_distances") == 0)
361     {
362     tab = tabulator_new_from_file_using_header(rfile,"ATOM1 ATOM2 MIN MAX F_CONST",0);
363     if (tab)
364     {
365     n_row = tabulator_height(tab);
366     for (i=0; i < n_row; i++)
367     {
368     ia1 = atoi(tab[i][0]);
369     ia2 = atoi(tab[i][1]);
370     min = atof(tab[i][2]);
371     max = atof(tab[i][3]);
372     fconst = atof(tab[i][4]);
373     if ((ia1 > 0 && ia1 <= natom) && (ia2 > 0 && ia2 < natom))
374     {
375     fx_dist.kdfix[fx_dist.ndfix][0] = ia1;
376     fx_dist.kdfix[fx_dist.ndfix][1] = ia2;
377     fx_dist.fdconst[fx_dist.ndfix] = fconst;
378     fx_dist.min_dist[fx_dist.ndfix] = min;
379     fx_dist.max_dist[fx_dist.ndfix] = max;
380     fx_dist.ndfix++;
381     }
382     }
383     }
384     } else if (strcasecmp(tag,"restrained_angles") == 0)
385     {
386     tab = tabulator_new_from_file_using_header(rfile,"ATOM1 ATOM2 ATOM3 MIN MAX F_CONST",0);
387     if (tab)
388     {
389     n_row = tabulator_height(tab);
390     for (i=0; i < n_row; i++)
391     {
392     ia1 = atoi(tab[i][0]);
393     ia2 = atoi(tab[i][1]);
394     ia3 = atoi(tab[i][2]);
395     min = atof(tab[i][3]);
396     max = atof(tab[i][4]);
397     fconst = atof(tab[i][5]);
398     if ((ia1 > 0 && ia1 <= natom) && (ia2 > 0 && ia2 < natom) && (ia3 > 0 && ia3 <= natom))
399     {
400     fx_angle.kafix[fx_angle.nafix][0] = ia1;
401     fx_angle.kafix[fx_angle.nafix][1] = ia2;
402     fx_angle.kafix[fx_angle.nafix][2] = ia3;
403     fx_angle.faconst[fx_angle.nafix] = fconst;
404     fx_angle.min_ang[fx_angle.nafix] = min;
405     fx_angle.max_ang[fx_angle.nafix] = max;
406     fx_angle.nafix++;
407     }
408     }
409     }
410     } else if (strcasecmp(tag,"restrained_dihedrals") == 0)
411     {
412     tab = tabulator_new_from_file_using_header(rfile,"ATOM1 ATOM2 ATOM3 ATOM4 MIN MAX F_CONST",0);
413     if (tab)
414     {
415     n_row = tabulator_height(tab);
416     for (i=0; i < n_row; i++)
417     {
418     ia1 = atoi(tab[i][0]);
419     ia2 = atoi(tab[i][1]);
420     ia3 = atoi(tab[i][2]);
421     ia4 = atoi(tab[i][3]);
422     min = atof(tab[i][4]);
423     max = atof(tab[i][5]);
424     fconst = atof(tab[i][6]);
425     if ((ia1 > 0 && ia1 <= natom) && (ia2 > 0 && ia2 < natom) && (ia3 > 0 && ia3 <= natom) && (ia4 > 0 && ia4 <= natom))
426     {
427     fx_torsion.ktfix[fx_torsion.ntfix][0] = ia1;
428     fx_torsion.ktfix[fx_torsion.ntfix][1] = ia2;
429     fx_torsion.ktfix[fx_torsion.ntfix][2] = ia3;
430     fx_torsion.ktfix[fx_torsion.ntfix][3] = ia4;
431     fx_torsion.ftconst[fx_torsion.ntfix] = fconst;
432     fx_torsion.min_tor[fx_torsion.ntfix] = min;
433     fx_torsion.max_tor[fx_torsion.ntfix] = max;
434     fx_torsion.ntfix++;
435     }
436     }
437     }
438     }
439     }
440     }
441     // if (has_aromatic)
442     // need to deal with aromatic bonds
443     L_DONE:
444     quick_type();
445     if (has_Aromatic)
446     {
447 wdelano 99 for (i=0; i < numbonds; i++)
448     mopaco(ib1[i],ib2[i]);
449 gilbertke 67 }
450 wdelano 99 free_ivector(ib1,0,nibond);
451     free_ivector(ib2,0,nibond);
452 gilbertke 67 return Ret_Val;
453     }
454     // =====================
455     void get_tag(char *line,char *tag)
456     {
457     int i,j,icount;
458     icount = 0;
459     strcpy(tag,"");
460     for (i=0; i < strlen(line); i++)
461     {
462     if (line[i] == '<')
463     {
464     for (j=i+1; j < strlen(line); j++)
465     {
466     if (line[j] == '>')
467     {
468     tag[icount] = '\0';
469     return;
470     } else
471     {
472     tag[icount] = line[j];
473     icount++;
474     }
475     }
476     }
477     }
478     }
479     /* ------------------------------------*/
480     // new mopaco version - test valency of each atom and add bond if possible
481     void mopaco(int ia, int ib)
482     {
483     int ibi, ibj, it, jt, kk;
484    
485     it = atom[ia].atomnum;
486     jt = atom[ib].atomnum;
487    
488     if( it == 6 || it == 7 || it == 8 || it == 15 || it == 16 )
489     {
490     if( jt == 6 || jt == 7 || jt == 8 || jt == 15 || jt == 16 )
491     {
492     /* c=n */
493     if( (it == 7 && jt == 6) || (it == 6 && jt == 7) )
494     {
495     ibi = 0;
496     for( kk = 0; kk < MAXIAT; kk++ )
497     {
498     if( atom[ia].bo[kk] != 9 )
499     ibi = ibi + atom[ia].bo[kk];
500     }
501     ibj = 0;
502     for( kk = 0; kk < MAXIAT; kk++ )
503     {
504     if( atom[ib].bo[kk] != 9 )
505     ibj = ibj + atom[ib].bo[kk];
506     }
507     if( it == 7 )
508     {
509     if( ibi <= 2 && ibj <= 3 )
510     {
511     make_bond( ia, ib, 1 );
512     }
513     }else if( jt == 7 )
514     {
515     if( ibj <= 2 && ibi <= 3 )
516     {
517     make_bond( ia, ib, 1 );
518     }
519     }
520     /* c=c bond */
521     } else if( (it == 6 && jt == 6) )
522     {
523     ibi = 0;
524     for( kk = 0; kk < MAXIAT; kk++ )
525     {
526     if( atom[ia].bo[kk] != 9 )
527     ibi = ibi + atom[ia].bo[kk];
528     }
529     ibj = 0;
530     for( kk = 0; kk < MAXIAT; kk++ )
531     {
532     if( atom[ib].bo[kk] != 9 )
533     ibj = ibj + atom[ib].bo[kk];
534     }
535     if( ibi <= 3 && ibj <= 3 )
536     {
537     make_bond( ia, ib, 1 );
538     }
539     }
540     }
541     }
542     }
543