ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/freemol/trunk/src/mengine/src/read_sdf.c
Revision: 109
Committed: Fri Mar 6 21:18:36 2009 UTC (11 years, 4 months ago) by wdelano
File size: 20780 byte(s)
Log Message:
Line File contents
1 /* 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
7 #include "fix.h"
8 #include "draw.h"
9 #include "job_control.h"
10 #include "vibrate.h"
11 #include "utility.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 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
34 static char Struct_Title[100];
35
36 // ============================
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 void get_molecule_memory(int);
48 void get_rings(int natom,int **iat,long int *flags);
49 double get_dipole_moment(void);
50 double get_total_energy(void);
51 char *get_structure_title(void);
52 // ==================================
53 char *get_structure_title()
54 {
55 if (strlen(Struct_Title) > 0)
56 return Struct_Title;
57 else
58 return NULL;
59 }
60 // ==================================
61 /*
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 if( atom.mmx_type[i] == 20 )
75 {
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 if( atom.iat[j][i] != 0 )
88 {
89 if( atom.iat[j][i] < j )
90 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 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 {
124 if (atom.bo[i][0] == 1)
125 junk = 5;
126 }
127 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 }
130 for (i=1; i <= natom; i++)
131 {
132 for (j=0; j < MAXIAT; j++)
133 {
134 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 }
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 fprintf(wfile,"%f\n",get_total_energy());
144
145 /* fprintf(wfile,"> <MMFF94_Atomtypes>\n");
146 fprintf(wfile,"+ Atom Atomtype\n");
147 for (i=1; i <=natom; i++)
148 fprintf(wfile,"| %d %d\n",i,atom.type[i]);
149 fprintf(wfile,"\n"); */
150
151 if (flags & DO_DIPOLE) {
152 fprintf(wfile,"> <dipole moment>\n");
153 fprintf(wfile,"%f\n",get_dipole_moment());
154 }
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 int numbonds, *ib1, *ib2;
191 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 // get memory for molecule
230 get_molecule_memory(niatom);
231 // allocate space for aromatic bonds
232 numbonds = 0;
233 ib1 = ivector(0,nibond);
234 ib2 = ivector(0,nibond);
235 //
236 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 atom.mmx_type[newatom] = itype;
263 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 ib1[numbonds] = ia1;
292 ib2[numbonds] = ia2;
293 numbonds++;
294 has_Aromatic = TRUE;
295 } else
296 make_bond(ia1, ia2, ibond);
297 if (istereo == 1 && planar)
298 {
299 atom.z[ia2] += 0.7;
300 if (atom.atomnum[ia2] == 1)
301 atom.type[ia2] = 60;
302 }
303 if (istereo == 6 && planar)
304 {
305 atom.z[ia2] -= 0.7;
306 if (atom.atomnum[ia2] == 1)
307 atom.type[ia2] = 60;
308 }
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 // <ELECTROSTATICS>
319 // <PARTIAL_CHARGES>
320
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 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 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 } else if (strcasecmp(tag,"electrostatics") == 0)
462 {
463 tab = tabulator_new_from_file_using_header(rfile,"TREATMENT DIELECTRIC SCALE_FACTOR GBSA_INTERNAL GBSA_EXTERNAL",0);
464 if (tabulator_height(tab))
465 {
466 if (strcmp(tab[0][0],"NONE") == 0)
467 {
468 job_control.use_charge = TRUE;
469 job_control.scale = 0.0;
470 } else if (strcmp(tab[0][0],"NORMAL") == 0)
471 {
472 if(tab[0][1])
473 units.dielec = atof(tab[0][1]);
474 else
475 units.dielec = 1.0F;
476 } else if (strcmp(tab[0][0],"SCALED") == 0)
477 {
478 job_control.use_charge = TRUE;
479 if(tab[0][2])
480 job_control.scale = atof(tab[0][2]);
481 else
482 job_control.scale = 1.0F;
483 } else if (strcmp(tab[0][0],"GBSA") == 0)
484 {
485 job_control.use_gbsa = TRUE;
486 solvent.type = 1; // STILL
487 if(tab[0][3])
488 solvent.EPSin = atof(tab[0][3]);
489 else
490 solvent.EPSin = 1.0F;
491 if(tab[0][4])
492 solvent.EPSsolv = atof(tab[0][4]);
493 else
494 solvent.EPSsolv = 78.30F;
495 }
496 }
497 }
498 }
499 }
500 // if (has_aromatic)
501 // need to deal with aromatic bonds
502 L_DONE:
503 get_rings(natom,atom.iat,atom.flags);
504 quick_type();
505 if (has_Aromatic)
506 {
507 for (i=0; i < numbonds; i++)
508 mopaco(ib1[i],ib2[i]);
509 }
510 free_ivector(ib1,0,nibond);
511 free_ivector(ib2,0,nibond);
512 return Ret_Val;
513 }
514 // =====================
515 void get_tag(char *line,char *tag)
516 {
517 int i,j,icount;
518 icount = 0;
519 strcpy(tag,"");
520 for (i=0; i < strlen(line); i++)
521 {
522 if (line[i] == '<')
523 {
524 for (j=i+1; j < strlen(line); j++)
525 {
526 if (line[j] == '>')
527 {
528 tag[icount] = '\0';
529 return;
530 } else
531 {
532 tag[icount] = line[j];
533 icount++;
534 }
535 }
536 }
537 }
538 }
539 /* ------------------------------------*/
540 // new mopaco version - test valency of each atom and add bond if possible
541 void mopaco(int ia, int ib)
542 {
543 int ibi, ibj, it, jt, kk;
544
545 it = atom.atomnum[ia];
546 jt = atom.atomnum[ib];
547
548 if( it == 6 || it == 7 || it == 8 || it == 15 || it == 16 )
549 {
550 if( jt == 6 || jt == 7 || jt == 8 || jt == 15 || jt == 16 )
551 {
552 /* c=n */
553 if( (it == 7 && jt == 6) || (it == 6 && jt == 7) )
554 {
555 ibi = 0;
556 for( kk = 0; kk < MAXIAT; kk++ )
557 {
558 if( atom.bo[ia][kk] != 9 )
559 ibi = ibi + atom.bo[ia][kk];
560 }
561 ibj = 0;
562 for( kk = 0; kk < MAXIAT; kk++ )
563 {
564 if( atom.bo[ib][kk] != 9 )
565 ibj = ibj + atom.bo[ib][kk];
566 }
567 if( it == 7 )
568 {
569 if( ibi <= 2 && ibj <= 3 )
570 {
571 make_bond( ia, ib, 1 );
572 }
573 }else if( jt == 7 )
574 {
575 if( ibj <= 2 && ibi <= 3 )
576 {
577 make_bond( ia, ib, 1 );
578 }
579 }
580 /* c=c bond */
581 } else if( (it == 6 && jt == 6) )
582 {
583 ibi = 0;
584 for( kk = 0; kk < MAXIAT; kk++ )
585 {
586 if( atom.bo[ia][kk] != 9 )
587 ibi = ibi + atom.bo[ia][kk];
588 }
589 ibj = 0;
590 for( kk = 0; kk < MAXIAT; kk++ )
591 {
592 if( atom.bo[ib][kk] != 9 )
593 ibj = ibj + atom.bo[ib][kk];
594 }
595 if( ibi <= 3 && ibj <= 3 )
596 {
597 make_bond( ia, ib, 1 );
598 }
599 }
600 }
601 }
602 }
603