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