ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/freemol/trunk/src/mengine/src/read_sdf.c
Revision: 62
Committed: Tue Dec 2 07:28:56 2008 UTC (12 years ago) by wdelano
File size: 13322 byte(s)
Log Message:
fixed typos
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
10 #define TABULATOR_INCLUDE_IMPLEMENTATION
11 #include "tabulator.h"
12
13 EXTERN struct t_files {
14 int nfiles, append, batch, icurrent, ibatno;
15 } files;
16 EXTERN struct ElementType { char symbol[3];
17 int atomnum;
18 float weight, covradius, vdwradius;
19 int s,p,d,f,type ;
20 } Elements[];
21 EXTERN struct t_dipolemom {
22 double total, xdipole, ydipole, zdipole;
23 } dipolemom;
24 EXTERN struct t_logp {
25 float logp;
26 } logp_calc;
27 EXTERN struct t_vibdata {
28 char ptgrp[4];
29 float mom_ix,mom_iy,mom_iz;
30 float etot,htot,stot,gtot,cptot;
31 } vibdata;
32
33 // ============================
34
35 int make_atom(int, float, float, float,char *);
36 void make_bond(int, int, int);
37 int read_sdf(int,int);
38 int rd_sdf(FILE *);
39 void write_sdf(int);
40 void message_alert(char *, char *);
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 // ==================================
47 /*
48 * flags - indicates whether dipole, xlogp or vibrational calculations were done
49 * and if so write them to the output file. Look at the definitions in pcmod.h
50 *
51 */
52 void write_sdf(int flags)
53 {
54 int i,j,lptest,nbond,junk;
55 FILE *wfile = pcmoutfile;
56
57 lptest = 0;
58 for( i = 1; i <= natom; i++ )
59 {
60 if( atom[i].mmx_type == 20 )
61 {
62 lptest = 1;
63 hdel( lptest );
64 break;
65 }
66 }
67 nbond = 0;
68 /* ** calculate the number of bonds in the molecule ** */
69 for( j = 1; j <= natom; j++ )
70 {
71 for( i = 0; i < MAXIAT; i++ )
72 {
73 if( atom[j].iat[i] != 0 )
74 {
75 if( atom[j].iat[i] < j )
76 nbond = nbond + 1;
77 }
78 }
79 }
80 /* now write the concord file */
81 /*
82 if (files.append)
83 wfile = fopen_path(Savebox.path,Savebox.fname,"a");
84 else
85 wfile = fopen_path(Savebox.path,Savebox.fname,"w");
86 */
87
88 j = strlen(Struct_Title);
89 for (i=0; i < j; i++)
90 {
91 if (Struct_Title[i] == '\n')
92 {
93 Struct_Title[i] = '\0';
94 break;
95 }
96 }
97 fprintf(wfile,"%s\n",Struct_Title);
98 fprintf(wfile," PCMODEL v9.1 1.00000 0.00000\n");
99 fprintf(wfile,"\n");
100
101 fprintf(wfile,"%3d%3d 0 0 0 0 1 V2000\n",natom,nbond);
102
103 for (i=1; i <= natom; i++)
104 {
105 junk = 0;
106 if (atom[i].mmx_type == 41) junk = 3;
107 else if (atom[i].mmx_type == 42) junk = 5;
108 else if (atom[i].mmx_type == 66)
109 {
110 if (atom[i].bo[0] == 1)
111 junk = 5;
112 }
113 fprintf(wfile,"%10.4f%10.4f%10.4f %-3s 0 %d 0 0 0 0 \n",atom[i].x, atom[i].y,
114 atom[i].z,Elements[atom[i].atomnum-1].symbol,junk);
115 }
116 for (i=1; i <= natom; i++)
117 {
118 for (j=0; j < MAXIAT; j++)
119 {
120 if (atom[i].iat[j] != 0 && i < atom[i].iat[j])
121 fprintf(wfile,"%3d%3d%3d 0\n",i, atom[i].iat[j], atom[i].bo[j]);
122 }
123 }
124 fprintf(wfile,"M END\n");
125 fprintf(wfile,"> <title>\n");
126 fprintf(wfile,"%s\n",Struct_Title);
127
128 fprintf(wfile,"> <MMFF94 energy>\n");
129 fprintf(wfile,"%f\n",energies.total);
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 /* =================================== */
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 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
210 for (i=0; i < niatom; i++)
211 {
212 FetchRecord(rfile,inputline);
213 sscanf(inputline,"%f %f %f %s %d %d",&xtmp, &ytmp, &ztmp, atomchar,&junk,&junk1);
214
215 itype = 0;
216 if (xtmp != 0.0) xPlus= TRUE;
217 if (ytmp != 0.0) yPlus= TRUE;
218 if (ztmp != 0.0) zPlus= TRUE;
219
220 iz = strlen(atomchar);
221 if (itype == 0 && (atomchar[0] == 'N' || atomchar[0] == 'O') && iz == 1) // nitro fix
222 {
223 if (atomchar[0] == 'N') itype = 8;
224 if (atomchar[0] == 'O') itype = 6;
225
226 if (junk1 != 0)
227 {
228 if (junk1 == 3 && itype == 8)
229 itype = 41; // N+
230 if (junk1 == 5 && itype == 6)
231 itype = 42; // O-
232 }
233 }
234
235 newatom = make_atom(itype,xtmp,ytmp,ztmp,atomchar);
236 if (itype != 0)
237 atom[newatom].mmx_type = itype;
238 if (newatom == -1)
239 {
240 fprintf(pcmlogfile, "Atom %d %s not recognized structure skipped\n",i,atomchar);
241 Ret_Val = FALSE;
242 }
243 }
244 has_Aromatic = FALSE;
245
246 planar = TRUE;
247 if (xPlus && yPlus && zPlus) planar = FALSE;
248
249 for (i=0; i < nibond; i++)
250 {
251 FetchRecord(rfile,inputline);
252 for (j=0; j <4; j++)
253 {
254 c1[j] = inputline[j];
255 c2[j] = inputline[j+3];
256 }
257 c1[3] = '\0';c2[3] = '\0';
258 ia1 = atoi(c1);
259 ia2 = atoi(c2);
260 istereo = 0;
261 sscanf(inputline,"%3d%3d%3d%3d",&junk, &junk1, &ibond,&istereo);
262 if (ibond >= 4)
263 {
264 ibond = 1;
265 has_Aromatic = TRUE;
266 // TJO
267 Ret_Val = FALSE;
268 }
269 make_bond(ia1, ia2, ibond);
270 if (istereo == 1 && planar)
271 {
272 atom[ia2].z += 0.7;
273 if (atom[ia2].atomnum == 1)
274 atom[ia2].type = 60;
275 }
276 if (istereo == 6 && planar)
277 {
278 atom[ia2].z -= 0.7;
279 if (atom[ia2].atomnum == 1)
280 atom[ia2].type = 60;
281 }
282 }
283 // read to end of structure
284 // parse any tags here
285 // agreed tags
286 // <FIXED_ATOMS>
287 // <RESTRAINED_ATOMS>
288 // <RESTRAINED_DISTANCES>
289 // <RESTRAINED_ANGLES>
290 // <RESTRAINED_DIHEDRALS>
291
292 while (FetchRecord(rfile,inputline))
293 {
294 if (strncasecmp(inputline,"$$$$",4) == 0)
295 {
296 return Ret_Val;
297 } else if (inputline[0] == '>')
298 {
299 get_tag(inputline,tag);
300 if (strcasecmp(tag,"fixed_atoms") == 0)
301 {
302 tab = tabulator_new_from_file_using_header(rfile,"ATOM",0);
303 if (tab)
304 {
305 fprintf(pcmlogfile,"got table\n");
306 n_row = tabulator_height(tab);
307 for (i=0; i < n_row; i++)
308 {
309 ia1 = atoi(tab[i][0]);
310 if (ia1 > 0 && ia1 <= natom)
311 {
312 fx_atom.katom_fix[fx_atom.natom_fix] = ia1;
313 fx_atom.natom_fix++;
314 }
315 }
316 }
317 } else if (strcasecmp(tag,"restrained_atoms") == 0)
318 {
319 tab = tabulator_new_from_file_using_header(rfile,"ATOM MAX F_CONST X Y Z",0);
320 if (tab)
321 {
322 n_row = tabulator_height(tab);
323 for (i=0; i < n_row; i++)
324 {
325 ia1 = atoi(tab[i][0]);
326 max = atof(tab[i][1]);
327 fconst = atof(tab[i][2]);
328 if (ia1 > 0 && ia1 <= natom)
329 {
330 restrain_atom.katom_restrain[restrain_atom.natom_restrain] = ia1;
331 restrain_atom.restrain_const[restrain_atom.natom_restrain] = fconst;
332 restrain_atom.restrain_max[restrain_atom.natom_restrain] = max;
333 restrain_atom.restrain_position[restrain_atom.natom_restrain][0] = atom[ia1].x; // default positions
334 restrain_atom.restrain_position[restrain_atom.natom_restrain][1] = atom[ia1].y;
335 restrain_atom.restrain_position[restrain_atom.natom_restrain][2] = atom[ia1].z;
336 if (strlen(tab[i][3]) > 0) // x postion
337 {
338 xtmp = atof(tab[i][3]);
339 restrain_atom.restrain_position[restrain_atom.natom_restrain][0] = xtmp;
340 }
341 if (strlen(tab[i][4]) > 0) // y postion
342 {
343 xtmp = atof(tab[i][4]);
344 restrain_atom.restrain_position[restrain_atom.natom_restrain][1] = xtmp;
345 }
346 if (strlen(tab[i][5]) > 0) // y postion
347 {
348 xtmp = atof(tab[i][5]);
349 restrain_atom.restrain_position[restrain_atom.natom_restrain][2] = xtmp;
350 }
351 restrain_atom.natom_restrain++;
352 }
353 }
354 }
355 } else if (strcasecmp(tag,"restrained_distances") == 0)
356 {
357 tab = tabulator_new_from_file_using_header(rfile,"ATOM1 ATOM2 MIN MAX F_CONST",0);
358 if (tab)
359 {
360 n_row = tabulator_height(tab);
361 for (i=0; i < n_row; i++)
362 {
363 ia1 = atoi(tab[i][0]);
364 ia2 = atoi(tab[i][1]);
365 min = atof(tab[i][2]);
366 max = atof(tab[i][3]);
367 fconst = atof(tab[i][4]);
368 if ((ia1 > 0 && ia1 <= natom) && (ia2 > 0 && ia2 < natom))
369 {
370 fx_dist.kdfix[fx_dist.ndfix][0] = ia1;
371 fx_dist.kdfix[fx_dist.ndfix][1] = ia2;
372 fx_dist.fdconst[fx_dist.ndfix] = fconst;
373 fx_dist.min_dist[fx_dist.ndfix] = min;
374 fx_dist.max_dist[fx_dist.ndfix] = max;
375 fx_dist.ndfix++;
376 }
377 }
378 }
379 } else if (strcasecmp(tag,"restrained_angles") == 0)
380 {
381 tab = tabulator_new_from_file_using_header(rfile,"ATOM1 ATOM2 ATOM3 MIN MAX F_CONST",0);
382 if (tab)
383 {
384 n_row = tabulator_height(tab);
385 for (i=0; i < n_row; i++)
386 {
387 ia1 = atoi(tab[i][0]);
388 ia2 = atoi(tab[i][1]);
389 ia3 = atoi(tab[i][2]);
390 min = atof(tab[i][3]);
391 max = atof(tab[i][4]);
392 fconst = atof(tab[i][5]);
393 if ((ia1 > 0 && ia1 <= natom) && (ia2 > 0 && ia2 < natom) && (ia3 > 0 && ia3 <= natom))
394 {
395 fx_angle.kafix[fx_angle.nafix][0] = ia1;
396 fx_angle.kafix[fx_angle.nafix][1] = ia2;
397 fx_angle.kafix[fx_angle.nafix][2] = ia3;
398 fx_angle.faconst[fx_angle.nafix] = fconst;
399 fx_angle.min_ang[fx_angle.nafix] = min;
400 fx_angle.max_ang[fx_angle.nafix] = max;
401 fx_angle.nafix++;
402 }
403 }
404 }
405 } else if (strcasecmp(tag,"restrained_dihedrals") == 0)
406 {
407 tab = tabulator_new_from_file_using_header(rfile,"ATOM1 ATOM2 ATOM3 ATOM4 MIN MAX F_CONST",0);
408 if (tab)
409 {
410 n_row = tabulator_height(tab);
411 for (i=0; i < n_row; i++)
412 {
413 ia1 = atoi(tab[i][0]);
414 ia2 = atoi(tab[i][1]);
415 ia3 = atoi(tab[i][2]);
416 ia4 = atoi(tab[i][3]);
417 min = atof(tab[i][4]);
418 max = atof(tab[i][5]);
419 fconst = atof(tab[i][6]);
420 if ((ia1 > 0 && ia1 <= natom) && (ia2 > 0 && ia2 < natom) && (ia3 > 0 && ia3 <= natom) && (ia4 > 0 && ia4 <= natom))
421 {
422 fx_torsion.ktfix[fx_torsion.ntfix][0] = ia1;
423 fx_torsion.ktfix[fx_torsion.ntfix][1] = ia2;
424 fx_torsion.ktfix[fx_torsion.ntfix][2] = ia3;
425 fx_torsion.ktfix[fx_torsion.ntfix][3] = ia4;
426 fx_torsion.ftconst[fx_torsion.ntfix] = fconst;
427 fx_torsion.min_tor[fx_torsion.ntfix] = min;
428 fx_torsion.max_tor[fx_torsion.ntfix] = max;
429 fx_torsion.ntfix++;
430 }
431 }
432 }
433 }
434 }
435 }
436 return Ret_Val;
437 }
438 // =====================
439 void get_tag(char *line,char *tag)
440 {
441 int i,j,icount;
442 icount = 0;
443 strcpy(tag,"");
444 for (i=0; i < strlen(line); i++)
445 {
446 if (line[i] == '<')
447 {
448 for (j=i+1; j < strlen(line); j++)
449 {
450 if (line[j] == '>')
451 {
452 tag[icount] = '\0';
453 return;
454 } else
455 {
456 tag[icount] = line[j];
457 icount++;
458 }
459 }
460 }
461 }
462 }