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