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