ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/freemol/trunk/src/mengine/src/read_sdf.c
Revision: 64
Committed: Thu Dec 4 02:03:23 2008 UTC (12 years, 9 months ago) by gilbertke
File size: 20816 byte(s)
Log Message:
fixed bug in routine to recognize delocalized bonds in kbond.c
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 #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
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
28 // ============================
29
30 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 void get_tag(char *,char *);
38 static void mopaco(int start_atom,int end_atom);
39 // ==================================
40 /*
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 FILE *wfile = pcmoutfile;
49
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 /*
75 if (files.append)
76 wfile = fopen_path(Savebox.path,Savebox.fname,"a");
77 else
78 wfile = fopen_path(Savebox.path,Savebox.fname,"w");
79 */
80
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 fprintf(wfile," MENGINE v1.0 1.00000 0.00000\n");
92 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 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 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 // fclose(wfile);
155 }
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 float xtmp, ytmp, ztmp, dz;
171 float fconst,min,max;
172 char c1[4],c2[4];
173 char atomchar[3];
174 char inputline[150];
175 char tag[30];
176 char ***tab;
177
178 Ret_Val = TRUE;
179 got_title = FALSE;
180 xPlus = yPlus = zPlus = FALSE;
181 if ( 0 == FetchRecord(rfile,inputline) )return -1;
182 // sscanf(inputline,"SDF %s",Struct_Title);
183 strncpy(Struct_Title, inputline, sizeof(Struct_Title));
184 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 fprintf(pcmlogfile, "Atom %d %s not recognized structure skipped\n",i,atomchar);
239 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 // TJO
265 Ret_Val = FALSE;
266 }
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 // parse any tags here
283 // agreed tags
284 // <FIXED_ATOMS>
285 // <RESTRAINED_ATOMS>
286 // <RESTRAINED_DISTANCES>
287 // <RESTRAINED_ANGLES>
288 // <RESTRAINED_DIHEDRALS>
289
290 while (FetchRecord(rfile,inputline))
291 {
292 if (strncasecmp(inputline,"$$$$",4) == 0)
293 {
294 goto L_DONE;
295 } 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 tab = tabulator_new_from_file_using_header(rfile,"ATOM MAX F_CONST X Y Z",0);
318 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 xtmp = atof(tab[i][5]);
347 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 }
433 }
434 // if (has_aromatic)
435 // need to deal with aromatic bonds
436 L_DONE:
437 if (has_Aromatic)
438 {
439 mopaco(1,natom);
440 }
441 return Ret_Val;
442 }
443 // =====================
444 void get_tag(char *line,char *tag)
445 {
446 int i,j,icount;
447 icount = 0;
448 strcpy(tag,"");
449 for (i=0; i < strlen(line); i++)
450 {
451 if (line[i] == '<')
452 {
453 for (j=i+1; j < strlen(line); j++)
454 {
455 if (line[j] == '>')
456 {
457 tag[icount] = '\0';
458 return;
459 } else
460 {
461 tag[icount] = line[j];
462 icount++;
463 }
464 }
465 }
466 }
467 }
468 /* ------------------------------------*/
469 // currently based on distance tests
470 // this should be changed to only look at bonds passed in with bond order >= 4
471 // but we need a mechanism to mark this -- Kevin
472 void mopaco(int istart,int ifinish)
473 {
474 int i, ibi, ibj, it, j, jt, kk;
475 float disij, dx, dy, dz, r, xii, yii, zii,rbdis[4];
476
477 rbdis[0] = 1.42;
478 rbdis[1] = 1.25;
479 rbdis[2] = 1.37;
480 rbdis[3] = 1.25;
481
482 /* find all double and triple bonds - search only types for
483 * carbon, oxygen, nitrogen, sulfur, phosphorous */
484 for( i = istart; i <= (ifinish - 1); i++ )
485 {
486 it = atom[i].atomnum;
487 if( it == 6 || it == 7 || it == 8 || it == 15 || it == 16 )
488 {
489 for( j = 0; j < MAXIAT; j++ )
490 {
491 if( atom[i].iat[j] != 0 )
492 {
493 jt = atom[atom[i].iat[j]].atomnum;
494 if( jt == 6 || jt == 7 || jt == 8 || jt == 15 || jt == 16 )
495 {
496 xii = atom[i].x - atom[atom[i].iat[j]].x;
497 yii = atom[i].y - atom[atom[i].iat[j]].y;
498 zii = atom[i].z - atom[atom[i].iat[j]].z;
499 if( fabs( xii ) <= 2.3 )
500 {
501 if( fabs( yii ) <= 2.3 )
502 {
503 if( fabs( zii ) <= 2.3 )
504 {
505 disij = sqrt( xii*xii + yii*yii + zii*zii );
506 /* c=o bond */
507 if( ((it == 8 && jt == 6) || (it == 6 && jt == 8)) && disij <= rbdis[3] + 0.05 )
508 {
509 if( i < atom[i].iat[j] )
510 make_bond( i, atom[i].iat[j],1 );
511 /* alkyne */
512 }else if( (it == 6 && jt == 6) && disij <= rbdis[1] )
513 {
514 if( i < atom[i].iat[j] )
515 make_bond( i, atom[i].iat[j],2 );
516 /* c=n */
517 }else if( ((it == 7 && jt == 6) || (it == 6 && jt == 7))
518 && disij <= rbdis[2] )
519 {
520 ibi = 0;
521 for( kk = 0; kk < MAXIAT; kk++ )
522 {
523 if( atom[i].bo[kk] != 9 )
524 ibi = ibi + atom[i].bo[kk];
525 }
526 ibj = 0;
527 for( kk = 0; kk < MAXIAT; kk++ )
528 {
529 if( atom[atom[i].iat[j]].bo[kk] != 9 )
530 ibj = ibj + atom[atom[i].iat[j]].bo[kk];
531 }
532 if( it == 7 )
533 {
534 if( ibi <= 2 && ibj <= 3 )
535 {
536 if( i < atom[i].iat[j] )
537 make_bond( i, atom[i].iat[j], 1 );
538 break;
539 }
540 }else if( jt == 7 )
541 {
542 if( ibj <= 2 && ibi <= 3 )
543 {
544 if( i < atom[i].iat[j] )
545 make_bond( i, atom[i].iat[j], 1 );
546 break;
547 }
548 }
549 /*nitrile */ }else if( ((it == 7 && jt == 6) || (it == 6 && jt == 7)) && disij <= 1.18 )
550 {
551 if( i < atom[i].iat[j] )
552 make_bond( i, atom[i].iat[j], 2 );
553 /* N=O */ }else if ( ((it == 7 && jt == 8) || (it == 8 && jt == 7)) && disij <= 1.26)
554 {
555 ibi = 0;
556 if (it == 7)
557 {
558 for (kk = 0; kk < MAXIAT; kk++)
559 {
560 if (atom[i].bo[kk] != 9)
561 ibi += atom[i].bo[kk];
562 }
563 } else
564 {
565 for (kk = 0; kk < MAXIAT; kk++)
566 {
567 if (atom[atom[i].iat[j]].bo[kk] != 9)
568 ibi += atom[atom[i].iat[j]].bo[kk];
569 }
570 }
571 if (ibi <= 3)
572 {
573 if ( i < atom[i].iat[j])
574 make_bond(i,atom[i].iat[j],1);
575 break;
576 }
577 /* c=c bond */ }else if( (it == 6 && jt == 6) && disij <= rbdis[0] )
578 {
579 ibi = 0;
580 for( kk = 0; kk < MAXIAT; kk++ )
581 {
582 if( atom[i].bo[kk] != 9 )
583 ibi = ibi + atom[i].bo[kk];
584 }
585 ibj = 0;
586 for( kk = 0; kk < MAXIAT; kk++ )
587 {
588 if( atom[atom[i].iat[j]].bo[kk] != 9 )
589 ibj = ibj + atom[atom[i].iat[j]].bo[kk];
590 }
591 if( ibi <= 3 && ibj <= 3 )
592 {
593 if( i < atom[i].iat[j] )
594 make_bond( i, atom[i].iat[j], 1 );
595 }
596 }
597 }
598 }
599 }
600 }
601 }
602 }
603 }
604 }
605 // type();
606 return;
607 }