ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/freemol/trunk/src/mengine/src/read_sdf.c
Revision: 65
Committed: Fri Dec 5 02:21:07 2008 UTC (13 years, 7 months ago) by gilbertke
File size: 20857 byte(s)
Log Message:
added call to type in rd_sdf, fixed problem in type_mmx
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 gilbertke 63 #include "dipmom.h"
10     #include "draw.h"
11     #include "vibrate.h"
12 tjod 3
13 wdelano 58 #define TABULATOR_INCLUDE_IMPLEMENTATION
14     #include "tabulator.h"
15    
16 tjod 3 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 wdelano 58 // ============================
29 gilbertke 65 void quick_type(void);
30 tjod 3 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 wdelano 58 void get_tag(char *,char *);
38 gilbertke 63 static void mopaco(int start_atom,int end_atom);
39 wdelano 58 // ==================================
40 tjod 3 /*
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 wdelano 58 FILE *wfile = pcmoutfile;
49 tjod 3
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 tjod 15 /*
75 tjod 3 if (files.append)
76     wfile = fopen_path(Savebox.path,Savebox.fname,"a");
77     else
78     wfile = fopen_path(Savebox.path,Savebox.fname,"w");
79 tjod 15 */
80 tjod 3
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 gilbertke 63 fprintf(wfile," MENGINE v1.0 1.00000 0.00000\n");
92 tjod 3 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 gilbertke 63 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 tjod 3 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 tjod 15 // fclose(wfile);
155 tjod 3 }
156 wdelano 58
157 tjod 3 /* =================================== */
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 wdelano 58 int i, j, niatom, nibond, ia1, ia2, ia3, ia4, ibond, itype, newatom;
164     int ncount,junk,junk1,got_title,istereo,iz,nvalue;
165 tjod 3 int jji, jj1, jj2, jj3, jj4;
166 wdelano 58 int n_row;
167 tjod 3 int has_Aromatic, xPlus,yPlus,zPlus, planar;
168     int icount,itemp[4];
169     int Ret_Val;
170     float xtmp, ytmp, ztmp, dz;
171 wdelano 58 float fconst,min,max;
172 tjod 3 char c1[4],c2[4];
173     char atomchar[3];
174     char inputline[150];
175 wdelano 58 char tag[30];
176     char ***tab;
177 tjod 3
178     Ret_Val = TRUE;
179     got_title = FALSE;
180     xPlus = yPlus = zPlus = FALSE;
181 tjod 15 if ( 0 == FetchRecord(rfile,inputline) )return -1;
182     // sscanf(inputline,"SDF %s",Struct_Title);
183     strncpy(Struct_Title, inputline, sizeof(Struct_Title));
184 tjod 3 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 wdelano 58 fprintf(pcmlogfile, "Atom %d %s not recognized structure skipped\n",i,atomchar);
239 tjod 3 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     ibond = 1;
263     has_Aromatic = TRUE;
264 tjod 31 // TJO
265     Ret_Val = FALSE;
266 tjod 3 }
267     make_bond(ia1, ia2, ibond);
268     if (istereo == 1 && planar)
269     {
270     atom[ia2].z += 0.7;
271     if (atom[ia2].atomnum == 1)
272     atom[ia2].type = 60;
273     }
274     if (istereo == 6 && planar)
275     {
276     atom[ia2].z -= 0.7;
277     if (atom[ia2].atomnum == 1)
278     atom[ia2].type = 60;
279     }
280     }
281     // read to end of structure
282 wdelano 58 // parse any tags here
283     // agreed tags
284     // <FIXED_ATOMS>
285     // <RESTRAINED_ATOMS>
286     // <RESTRAINED_DISTANCES>
287     // <RESTRAINED_ANGLES>
288     // <RESTRAINED_DIHEDRALS>
289    
290 tjod 3 while (FetchRecord(rfile,inputline))
291     {
292     if (strncasecmp(inputline,"$$$$",4) == 0)
293     {
294 gilbertke 64 goto L_DONE;
295 wdelano 58 } else if (inputline[0] == '>')
296     {
297     get_tag(inputline,tag);
298     if (strcasecmp(tag,"fixed_atoms") == 0)
299     {
300     tab = tabulator_new_from_file_using_header(rfile,"ATOM",0);
301     if (tab)
302     {
303     fprintf(pcmlogfile,"got table\n");
304     n_row = tabulator_height(tab);
305     for (i=0; i < n_row; i++)
306     {
307     ia1 = atoi(tab[i][0]);
308     if (ia1 > 0 && ia1 <= natom)
309     {
310     fx_atom.katom_fix[fx_atom.natom_fix] = ia1;
311     fx_atom.natom_fix++;
312     }
313     }
314     }
315     } else if (strcasecmp(tag,"restrained_atoms") == 0)
316     {
317 wdelano 62 tab = tabulator_new_from_file_using_header(rfile,"ATOM MAX F_CONST X Y Z",0);
318 wdelano 58 if (tab)
319     {
320     n_row = tabulator_height(tab);
321     for (i=0; i < n_row; i++)
322     {
323     ia1 = atoi(tab[i][0]);
324     max = atof(tab[i][1]);
325     fconst = atof(tab[i][2]);
326     if (ia1 > 0 && ia1 <= natom)
327     {
328     restrain_atom.katom_restrain[restrain_atom.natom_restrain] = ia1;
329     restrain_atom.restrain_const[restrain_atom.natom_restrain] = fconst;
330     restrain_atom.restrain_max[restrain_atom.natom_restrain] = max;
331     restrain_atom.restrain_position[restrain_atom.natom_restrain][0] = atom[ia1].x; // default positions
332     restrain_atom.restrain_position[restrain_atom.natom_restrain][1] = atom[ia1].y;
333     restrain_atom.restrain_position[restrain_atom.natom_restrain][2] = atom[ia1].z;
334     if (strlen(tab[i][3]) > 0) // x postion
335     {
336     xtmp = atof(tab[i][3]);
337     restrain_atom.restrain_position[restrain_atom.natom_restrain][0] = xtmp;
338     }
339     if (strlen(tab[i][4]) > 0) // y postion
340     {
341     xtmp = atof(tab[i][4]);
342     restrain_atom.restrain_position[restrain_atom.natom_restrain][1] = xtmp;
343     }
344     if (strlen(tab[i][5]) > 0) // y postion
345     {
346 wdelano 62 xtmp = atof(tab[i][5]);
347 wdelano 58 restrain_atom.restrain_position[restrain_atom.natom_restrain][2] = xtmp;
348     }
349     restrain_atom.natom_restrain++;
350     }
351     }
352     }
353     } else if (strcasecmp(tag,"restrained_distances") == 0)
354     {
355     tab = tabulator_new_from_file_using_header(rfile,"ATOM1 ATOM2 MIN MAX F_CONST",0);
356     if (tab)
357     {
358     n_row = tabulator_height(tab);
359     for (i=0; i < n_row; i++)
360     {
361     ia1 = atoi(tab[i][0]);
362     ia2 = atoi(tab[i][1]);
363     min = atof(tab[i][2]);
364     max = atof(tab[i][3]);
365     fconst = atof(tab[i][4]);
366     if ((ia1 > 0 && ia1 <= natom) && (ia2 > 0 && ia2 < natom))
367     {
368     fx_dist.kdfix[fx_dist.ndfix][0] = ia1;
369     fx_dist.kdfix[fx_dist.ndfix][1] = ia2;
370     fx_dist.fdconst[fx_dist.ndfix] = fconst;
371     fx_dist.min_dist[fx_dist.ndfix] = min;
372     fx_dist.max_dist[fx_dist.ndfix] = max;
373     fx_dist.ndfix++;
374     }
375     }
376     }
377     } else if (strcasecmp(tag,"restrained_angles") == 0)
378     {
379     tab = tabulator_new_from_file_using_header(rfile,"ATOM1 ATOM2 ATOM3 MIN MAX F_CONST",0);
380     if (tab)
381     {
382     n_row = tabulator_height(tab);
383     for (i=0; i < n_row; i++)
384     {
385     ia1 = atoi(tab[i][0]);
386     ia2 = atoi(tab[i][1]);
387     ia3 = atoi(tab[i][2]);
388     min = atof(tab[i][3]);
389     max = atof(tab[i][4]);
390     fconst = atof(tab[i][5]);
391     if ((ia1 > 0 && ia1 <= natom) && (ia2 > 0 && ia2 < natom) && (ia3 > 0 && ia3 <= natom))
392     {
393     fx_angle.kafix[fx_angle.nafix][0] = ia1;
394     fx_angle.kafix[fx_angle.nafix][1] = ia2;
395     fx_angle.kafix[fx_angle.nafix][2] = ia3;
396     fx_angle.faconst[fx_angle.nafix] = fconst;
397     fx_angle.min_ang[fx_angle.nafix] = min;
398     fx_angle.max_ang[fx_angle.nafix] = max;
399     fx_angle.nafix++;
400     }
401     }
402     }
403     } else if (strcasecmp(tag,"restrained_dihedrals") == 0)
404     {
405     tab = tabulator_new_from_file_using_header(rfile,"ATOM1 ATOM2 ATOM3 ATOM4 MIN MAX F_CONST",0);
406     if (tab)
407     {
408     n_row = tabulator_height(tab);
409     for (i=0; i < n_row; i++)
410     {
411     ia1 = atoi(tab[i][0]);
412     ia2 = atoi(tab[i][1]);
413     ia3 = atoi(tab[i][2]);
414     ia4 = atoi(tab[i][3]);
415     min = atof(tab[i][4]);
416     max = atof(tab[i][5]);
417     fconst = atof(tab[i][6]);
418     if ((ia1 > 0 && ia1 <= natom) && (ia2 > 0 && ia2 < natom) && (ia3 > 0 && ia3 <= natom) && (ia4 > 0 && ia4 <= natom))
419     {
420     fx_torsion.ktfix[fx_torsion.ntfix][0] = ia1;
421     fx_torsion.ktfix[fx_torsion.ntfix][1] = ia2;
422     fx_torsion.ktfix[fx_torsion.ntfix][2] = ia3;
423     fx_torsion.ktfix[fx_torsion.ntfix][3] = ia4;
424     fx_torsion.ftconst[fx_torsion.ntfix] = fconst;
425     fx_torsion.min_tor[fx_torsion.ntfix] = min;
426     fx_torsion.max_tor[fx_torsion.ntfix] = max;
427     fx_torsion.ntfix++;
428     }
429     }
430     }
431     }
432 tjod 3 }
433     }
434 gilbertke 63 // if (has_aromatic)
435     // need to deal with aromatic bonds
436 gilbertke 64 L_DONE:
437 gilbertke 65 quick_type();
438 gilbertke 63 if (has_Aromatic)
439 gilbertke 64 {
440 gilbertke 63 mopaco(1,natom);
441 gilbertke 64 }
442 tjod 3 return Ret_Val;
443     }
444 wdelano 58 // =====================
445     void get_tag(char *line,char *tag)
446 tjod 3 {
447 wdelano 58 int i,j,icount;
448     icount = 0;
449     strcpy(tag,"");
450     for (i=0; i < strlen(line); i++)
451 tjod 3 {
452 wdelano 58 if (line[i] == '<')
453     {
454     for (j=i+1; j < strlen(line); j++)
455     {
456     if (line[j] == '>')
457     {
458     tag[icount] = '\0';
459     return;
460     } else
461     {
462     tag[icount] = line[j];
463     icount++;
464     }
465     }
466     }
467 tjod 3 }
468     }
469 gilbertke 63 /* ------------------------------------*/
470     // currently based on distance tests
471     // this should be changed to only look at bonds passed in with bond order >= 4
472     // but we need a mechanism to mark this -- Kevin
473     void mopaco(int istart,int ifinish)
474     {
475     int i, ibi, ibj, it, j, jt, kk;
476     float disij, dx, dy, dz, r, xii, yii, zii,rbdis[4];
477    
478     rbdis[0] = 1.42;
479     rbdis[1] = 1.25;
480     rbdis[2] = 1.37;
481     rbdis[3] = 1.25;
482    
483     /* find all double and triple bonds - search only types for
484     * carbon, oxygen, nitrogen, sulfur, phosphorous */
485     for( i = istart; i <= (ifinish - 1); i++ )
486     {
487     it = atom[i].atomnum;
488     if( it == 6 || it == 7 || it == 8 || it == 15 || it == 16 )
489     {
490     for( j = 0; j < MAXIAT; j++ )
491     {
492     if( atom[i].iat[j] != 0 )
493     {
494     jt = atom[atom[i].iat[j]].atomnum;
495     if( jt == 6 || jt == 7 || jt == 8 || jt == 15 || jt == 16 )
496     {
497     xii = atom[i].x - atom[atom[i].iat[j]].x;
498     yii = atom[i].y - atom[atom[i].iat[j]].y;
499     zii = atom[i].z - atom[atom[i].iat[j]].z;
500     if( fabs( xii ) <= 2.3 )
501     {
502     if( fabs( yii ) <= 2.3 )
503     {
504     if( fabs( zii ) <= 2.3 )
505     {
506     disij = sqrt( xii*xii + yii*yii + zii*zii );
507     /* c=o bond */
508     if( ((it == 8 && jt == 6) || (it == 6 && jt == 8)) && disij <= rbdis[3] + 0.05 )
509     {
510     if( i < atom[i].iat[j] )
511     make_bond( i, atom[i].iat[j],1 );
512     /* alkyne */
513     }else if( (it == 6 && jt == 6) && disij <= rbdis[1] )
514     {
515     if( i < atom[i].iat[j] )
516     make_bond( i, atom[i].iat[j],2 );
517     /* c=n */
518     }else if( ((it == 7 && jt == 6) || (it == 6 && jt == 7))
519     && disij <= rbdis[2] )
520     {
521     ibi = 0;
522     for( kk = 0; kk < MAXIAT; kk++ )
523     {
524     if( atom[i].bo[kk] != 9 )
525     ibi = ibi + atom[i].bo[kk];
526     }
527     ibj = 0;
528     for( kk = 0; kk < MAXIAT; kk++ )
529     {
530     if( atom[atom[i].iat[j]].bo[kk] != 9 )
531     ibj = ibj + atom[atom[i].iat[j]].bo[kk];
532     }
533     if( it == 7 )
534     {
535     if( ibi <= 2 && ibj <= 3 )
536     {
537     if( i < atom[i].iat[j] )
538     make_bond( i, atom[i].iat[j], 1 );
539     break;
540     }
541     }else if( jt == 7 )
542     {
543     if( ibj <= 2 && ibi <= 3 )
544     {
545     if( i < atom[i].iat[j] )
546     make_bond( i, atom[i].iat[j], 1 );
547     break;
548     }
549     }
550     /*nitrile */ }else if( ((it == 7 && jt == 6) || (it == 6 && jt == 7)) && disij <= 1.18 )
551     {
552     if( i < atom[i].iat[j] )
553     make_bond( i, atom[i].iat[j], 2 );
554     /* N=O */ }else if ( ((it == 7 && jt == 8) || (it == 8 && jt == 7)) && disij <= 1.26)
555     {
556     ibi = 0;
557     if (it == 7)
558     {
559     for (kk = 0; kk < MAXIAT; kk++)
560     {
561     if (atom[i].bo[kk] != 9)
562     ibi += atom[i].bo[kk];
563     }
564     } else
565     {
566     for (kk = 0; kk < MAXIAT; kk++)
567     {
568     if (atom[atom[i].iat[j]].bo[kk] != 9)
569     ibi += atom[atom[i].iat[j]].bo[kk];
570     }
571     }
572     if (ibi <= 3)
573     {
574     if ( i < atom[i].iat[j])
575     make_bond(i,atom[i].iat[j],1);
576     break;
577     }
578     /* c=c bond */ }else if( (it == 6 && jt == 6) && disij <= rbdis[0] )
579     {
580     ibi = 0;
581     for( kk = 0; kk < MAXIAT; kk++ )
582     {
583     if( atom[i].bo[kk] != 9 )
584     ibi = ibi + atom[i].bo[kk];
585     }
586     ibj = 0;
587     for( kk = 0; kk < MAXIAT; kk++ )
588     {
589     if( atom[atom[i].iat[j]].bo[kk] != 9 )
590     ibj = ibj + atom[atom[i].iat[j]].bo[kk];
591     }
592     if( ibi <= 3 && ibj <= 3 )
593     {
594     if( i < atom[i].iat[j] )
595     make_bond( i, atom[i].iat[j], 1 );
596     }
597     }
598     }
599     }
600     }
601     }
602     }
603     }
604     }
605     }
606     // type();
607     return;
608     }