ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/freemol/trunk/src/mengine/src/read_sdf.c
Revision: 31
Committed: Fri Jul 18 18:47:27 2008 UTC (11 years, 10 months ago) by tjod
File size: 17021 byte(s)
Log Message:
When aromatic bonds are encountered generate error on stderr and do not output sdf file.

Line File contents
1 #define EXTERN extern
2
3 #include "pcwin.h"
4 #include "pcmod.h"
5 #include "substr.h"
6 #include "atom_k.h"
7 #include "energies.h"
8
9 EXTERN struct t_files {
10 int nfiles, append, batch, icurrent, ibatno;
11 } files;
12 EXTERN struct ElementType { char symbol[3];
13 int atomnum;
14 float weight, covradius, vdwradius;
15 int s,p,d,f,type ;
16 } Elements[];
17 EXTERN struct t_dipolemom {
18 double total, xdipole, ydipole, zdipole;
19 } dipolemom;
20 EXTERN struct t_logp {
21 float logp;
22 } logp_calc;
23 EXTERN struct t_vibdata {
24 char ptgrp[4];
25 float mom_ix,mom_iy,mom_iz;
26 float etot,htot,stot,gtot,cptot;
27 } vibdata;
28
29 int make_atom(int, float, float, float,char *);
30 void make_bond(int, int, int);
31 int read_sdf(int,int);
32 int rd_sdf(FILE *);
33 void write_sdf(int);
34 void message_alert(char *, char *);
35 void hdel(int);
36 FILE * fopen_path ( char * , char * , char * ) ;
37 void read_schakal(int,int);
38 void write_schakal(void);
39 void mopaco(int,int);
40 int FetchRecord(FILE *, char *);
41 void avgleg(void);
42
43 /*
44 * flags - indicates whether dipole, xlogp or vibrational calculations were done
45 * and if so write them to the output file. Look at the definitions in pcmod.h
46 *
47 */
48 void write_sdf(int flags)
49 {
50 int i,j,lptest,nbond,junk;
51 FILE *wfile;
52
53 lptest = 0;
54 for( i = 1; i <= natom; i++ )
55 {
56 if( atom[i].mmx_type == 20 )
57 {
58 lptest = 1;
59 hdel( lptest );
60 break;
61 }
62 }
63 nbond = 0;
64 /* ** calculate the number of bonds in the molecule ** */
65 for( j = 1; j <= natom; j++ )
66 {
67 for( i = 0; i < MAXIAT; i++ )
68 {
69 if( atom[j].iat[i] != 0 )
70 {
71 if( atom[j].iat[i] < j )
72 nbond = nbond + 1;
73 }
74 }
75 }
76 /* now write the concord file */
77 /*
78 if (files.append)
79 wfile = fopen_path(Savebox.path,Savebox.fname,"a");
80 else
81 wfile = fopen_path(Savebox.path,Savebox.fname,"w");
82 */
83 wfile = stdout;
84
85 j = strlen(Struct_Title);
86 for (i=0; i < j; i++)
87 {
88 if (Struct_Title[i] == '\n')
89 {
90 Struct_Title[i] = '\0';
91 break;
92 }
93 }
94 fprintf(wfile,"%s\n",Struct_Title);
95 fprintf(wfile," PCMODEL v9.1 1.00000 0.00000\n");
96 fprintf(wfile,"\n");
97
98 fprintf(wfile,"%3d%3d 0 0 0 0 1 V2000\n",natom,nbond);
99
100 for (i=1; i <= natom; i++)
101 {
102 junk = 0;
103 if (atom[i].mmx_type == 41) junk = 3;
104 else if (atom[i].mmx_type == 42) junk = 5;
105 else if (atom[i].mmx_type == 66)
106 {
107 if (atom[i].bo[0] == 1)
108 junk = 5;
109 }
110 fprintf(wfile,"%10.4f%10.4f%10.4f %-3s 0 %d 0 0 0 0 \n",atom[i].x, atom[i].y,
111 atom[i].z,Elements[atom[i].atomnum-1].symbol,junk);
112 }
113 for (i=1; i <= natom; i++)
114 {
115 for (j=0; j < MAXIAT; j++)
116 {
117 if (atom[i].iat[j] != 0 && i < atom[i].iat[j])
118 fprintf(wfile,"%3d%3d%3d 0\n",i, atom[i].iat[j], atom[i].bo[j]);
119 }
120 }
121 fprintf(wfile,"M END\n");
122 fprintf(wfile,"> <title>\n");
123 fprintf(wfile,"%s\n",Struct_Title);
124
125 fprintf(wfile,"> <MMFF94 energy>\n");
126 fprintf(wfile,"%f\n",energies.total);
127
128 if (flags & DO_DIPOLE) {
129 fprintf(wfile,"> <dipole moment>\n");
130 fprintf(wfile,"%f\n",dipolemom.total);
131 }
132
133 if (flags & DO_XLOGP) {
134 fprintf(wfile,"> <xLogP>\n");
135 fprintf(wfile,"%f\n",logp_calc.logp);
136 }
137
138 if (flags & DO_VIBRATION) {
139 fprintf(wfile,"> <Point Group>\n");
140 fprintf(wfile,"%s\n",vibdata.ptgrp);
141
142 fprintf(wfile,"> <Moments of Inertia>\n");
143 fprintf(wfile,"%f8.3 %f8.3 %f8.3 \n",vibdata.mom_ix,vibdata.mom_iy,vibdata.mom_iz);
144
145 fprintf(wfile,"> <thermodynamics dE, dH, S, dG, CP>\n");
146 fprintf(wfile,"%f8.3 %f8.3 %f8.3 %f8.3 %f8.3 \n",vibdata.etot, vibdata.htot,
147 vibdata.stot,vibdata.gtot,vibdata.cptot);
148 }
149
150 fprintf(wfile,"\n");
151 fprintf(wfile,"$$$$\n");
152 // fclose(wfile);
153 }
154 /* =================================== */
155 // fast read - assume file is open and positioned
156 // read structure and down to end $$$$
157 // and return
158 int rd_sdf(FILE *rfile)
159 {
160 int i, j, niatom, nibond, ia1, ia2, ibond, itype, newatom;
161 int ncount,junk,junk1,got_title,istereo,iz;
162 int jji, jj1, jj2, jj3, jj4;
163 int has_Aromatic, xPlus,yPlus,zPlus, planar;
164 int icount,itemp[4];
165 int Ret_Val;
166 float xtmp, ytmp, ztmp, dz;
167 char c1[4],c2[4];
168 char atomchar[3];
169 char inputline[150];
170
171 Ret_Val = TRUE;
172 got_title = FALSE;
173 xPlus = yPlus = zPlus = FALSE;
174 if ( 0 == FetchRecord(rfile,inputline) )return -1;
175 // sscanf(inputline,"SDF %s",Struct_Title);
176 strncpy(Struct_Title, inputline, sizeof(Struct_Title));
177 got_title = TRUE;
178 /* if (strlen(inputline) > 4)
179 {
180 iz = strlen(inputline);
181 if (iz > 60) iz = 59;
182 for (i=4; i < iz; i++)
183 Struct_Title[i-4] = inputline[i];
184 Struct_Title[i] = '\0';
185 got_title = TRUE;
186 } */
187 FetchRecord(rfile,inputline); // blank line
188 FetchRecord(rfile,inputline); // blank line
189
190 FetchRecord(rfile,inputline); // natom and nbond
191 for (i=0; i <4; i++)
192 {
193 c1[i] = inputline[i];
194 c2[i] = inputline[i+3];
195 }
196 c1[3] = '\0';c2[3] = '\0';
197 niatom = atoi(c1);
198 nibond = atoi(c2);
199 if (niatom == 0)
200 return FALSE;
201
202 for (i=0; i < niatom; i++)
203 {
204 FetchRecord(rfile,inputline);
205 sscanf(inputline,"%f %f %f %s %d %d",&xtmp, &ytmp, &ztmp, atomchar,&junk,&junk1);
206
207 itype = 0;
208 if (xtmp != 0.0) xPlus= TRUE;
209 if (ytmp != 0.0) yPlus= TRUE;
210 if (ztmp != 0.0) zPlus= TRUE;
211
212 iz = strlen(atomchar);
213 if (itype == 0 && (atomchar[0] == 'N' || atomchar[0] == 'O') && iz == 1) // nitro fix
214 {
215 if (atomchar[0] == 'N') itype = 8;
216 if (atomchar[0] == 'O') itype = 6;
217
218 if (junk1 != 0)
219 {
220 if (junk1 == 3 && itype == 8)
221 itype = 41; // N+
222 if (junk1 == 5 && itype == 6)
223 itype = 42; // O-
224 }
225 }
226
227 newatom = make_atom(itype,xtmp,ytmp,ztmp,atomchar);
228 if (itype != 0)
229 atom[newatom].mmx_type = itype;
230 if (newatom == -1)
231 {
232 printf("Atom %d %s not recognized structure skipped\n",i,atomchar);
233 Ret_Val = FALSE;
234 }
235 }
236 has_Aromatic = FALSE;
237
238 planar = TRUE;
239 if (xPlus && yPlus && zPlus) planar = FALSE;
240
241 for (i=0; i < nibond; i++)
242 {
243 FetchRecord(rfile,inputline);
244 for (j=0; j <4; j++)
245 {
246 c1[j] = inputline[j];
247 c2[j] = inputline[j+3];
248 }
249 c1[3] = '\0';c2[3] = '\0';
250 ia1 = atoi(c1);
251 ia2 = atoi(c2);
252 istereo = 0;
253 sscanf(inputline,"%3d%3d%3d%3d",&junk, &junk1, &ibond,&istereo);
254 if (ibond >= 4)
255 {
256 ibond = 1;
257 has_Aromatic = TRUE;
258 // TJO
259 Ret_Val = FALSE;
260 }
261 make_bond(ia1, ia2, ibond);
262 if (istereo == 1 && planar)
263 {
264 atom[ia2].z += 0.7;
265 if (atom[ia2].atomnum == 1)
266 atom[ia2].type = 60;
267 }
268 if (istereo == 6 && planar)
269 {
270 atom[ia2].z -= 0.7;
271 if (atom[ia2].atomnum == 1)
272 atom[ia2].type = 60;
273 }
274 }
275 // read to end of structure
276 while (FetchRecord(rfile,inputline))
277 {
278 if (strncasecmp(inputline,"$$$$",4) == 0)
279 {
280 return Ret_Val;
281 }
282 }
283 return Ret_Val;
284 }
285 /* =================================== */
286 int read_sdf(int istruct, int isubred)
287 {
288 int i, j, niatom, nibond, ia1, ia2, ibond, itype, newatom;
289 int ncount, ibotptr,junk,junk1,got_title,istereo,iz;
290 int jji, jj1, jj2, jj3, jj4;
291 int has_Aromatic, xPlus,yPlus,zPlus, planar;
292 int icount,itemp[4];
293 float xtmp, ytmp, ztmp, dz;
294 char c1[4],c2[4];
295 char atomchar[3];
296 char inputline[150];
297 FILE *infile;
298
299 infile = fopen_path(Openbox.path,Openbox.fname,"r");
300
301 if (infile == NULL)
302 {
303 message_alert("Error opening Concord file","Concord Setup");
304 return FALSE;
305 }
306
307 if (isubred == 0)
308 ibotptr = 0;
309 else
310 {
311 ibotptr = natom;
312 substr.istract[substr.icurstr] = TRUE;
313 }
314
315 if (istruct > 1)
316 {
317 ncount = 0;
318 while (FetchRecord(infile,inputline))
319 {
320 if (strncasecmp(inputline,"$$$$",4) == 0)
321 {
322 ncount++;
323 if (ncount+1 == istruct)
324 goto L_1;
325 }
326 }
327 }
328 // start file read here
329 L_1:
330 got_title = FALSE;
331 xPlus = yPlus = zPlus = FALSE;
332 FetchRecord(infile,inputline);
333 sscanf(inputline,"SDF %s",Struct_Title);
334 got_title = TRUE;
335 /* if (strlen(inputline) > 4)
336 {
337 iz = strlen(inputline);
338 if (iz > 60) iz = 59;
339 for (i=4; i < iz; i++)
340 Struct_Title[i-4] = inputline[i];
341 Struct_Title[i] = '\0';
342 got_title = TRUE;
343 } */
344 FetchRecord(infile,inputline); // blank line
345 FetchRecord(infile,inputline); // blank line
346
347 FetchRecord(infile,inputline); // natom and nbond
348 for (i=0; i <4; i++)
349 {
350 c1[i] = inputline[i];
351 c2[i] = inputline[i+3];
352 }
353 c1[3] = '\0';c2[3] = '\0';
354 niatom = atoi(c1);
355 nibond = atoi(c2);
356 if (niatom == 0)
357 return FALSE;
358
359 for (i=0; i < niatom; i++)
360 {
361 FetchRecord(infile,inputline);
362 sscanf(inputline,"%f %f %f %s %d %d",&xtmp, &ytmp, &ztmp, atomchar,&junk,&junk1);
363
364 itype = 0;
365 if (xtmp != 0.0) xPlus= TRUE;
366 if (ytmp != 0.0) yPlus= TRUE;
367 if (ztmp != 0.0) zPlus= TRUE;
368
369 iz = strlen(atomchar);
370 if (itype == 0 && (atomchar[0] == 'N' || atomchar[0] == 'O') && iz == 1) // nitro fix
371 {
372 if (atomchar[0] == 'N') itype = 8;
373 if (atomchar[0] == 'O') itype = 6;
374
375 if (junk1 != 0)
376 {
377 if (junk1 == 3 && itype == 8)
378 itype = 41; // N+
379 if (junk1 == 5 && itype == 6)
380 itype = 42; // O-
381 }
382 }
383
384 newatom = make_atom(itype,xtmp,ytmp,ztmp,atomchar);
385 if (newatom == -1)
386 {
387 printf("Atom %d %s not recognized structure skipped\n",i,atomchar);
388 return FALSE;
389 }
390 if (isubred == 1)
391 {
392 atom[newatom].substr[0] = 0;
393 atom[newatom].substr[0] |= (1L << substr.icurstr);
394 }
395 }
396 has_Aromatic = FALSE;
397
398 planar = TRUE;
399 if (xPlus && yPlus && zPlus) planar = FALSE;
400
401 for (i=0; i < nibond; i++)
402 {
403 FetchRecord(infile,inputline);
404 for (j=0; j <4; j++)
405 {
406 c1[j] = inputline[j];
407 c2[j] = inputline[j+3];
408 }
409 c1[3] = '\0';c2[3] = '\0';
410 ia1 = atoi(c1);
411 ia2 = atoi(c2);
412 istereo = 0;
413 sscanf(inputline,"%3d%3d%3d%3d",&junk, &junk1, &ibond,&istereo);
414 if (ibond >= 4)
415 {
416 ibond = 1;
417 has_Aromatic = TRUE;
418 }
419 make_bond(ia1+ibotptr, ia2+ibotptr, ibond);
420 if (istereo == 1 && planar)
421 {
422 atom[ia2+ibotptr].z += 0.7;
423 if (atom[ia2+ibotptr].atomnum == 1)
424 atom[ia2+ibotptr].type = 60;
425 }
426 if (istereo == 6 && planar)
427 {
428 atom[ia2+ibotptr].z -= 0.7;
429 if (atom[ia2+ibotptr].atomnum == 1)
430 atom[ia2+ibotptr].type = 60;
431 }
432 }
433 FetchRecord(infile,inputline); // M END line
434 FetchRecord(infile,inputline);
435 if (got_title == FALSE)
436 strcpy(Struct_Title,inputline);
437 fclose(infile);
438 ncount = strlen(Struct_Title);
439 for (i=0; i < ncount; i++)
440 {
441 if (Struct_Title[i] == '\n')
442 {
443 Struct_Title[i] = '\0';
444 break;
445 }
446 }
447 // search for quaternary centers - check if flat
448 dz = 0.0;
449 for (i=1; i <= natom; i++)
450 {
451 dz += atom[i].z;
452 jji = jj1 = jj2 = jj3 = jj4 = 0;
453 if (atom[i].type != 0)
454 {
455 for (j=0; j < 6; j++)
456 {
457 if (atom[i].iat[j] != 0)
458 jji++;
459 }
460 }
461 if (jji == 4)
462 {
463 if (atom[atom[i].iat[0]].z == 0.0 && atom[atom[i].iat[0]].z == 0.0 &&
464 atom[atom[i].iat[0]].z == 0.0 && atom[atom[i].iat[0]].z == 0.0) // flat center
465 {
466 for (j=0; j < 4; j++)
467 {
468 itemp[j] = 0;
469 if (atom[atom[i].iat[0]].iat[j] != 0) jj1++;
470 if (atom[atom[i].iat[1]].iat[j] != 0) jj2++;
471 if (atom[atom[i].iat[2]].iat[j] != 0) jj3++;
472 if (atom[atom[i].iat[3]].iat[j] != 0) jj4++;
473 }
474 icount = 0;
475 if (jj1 == 1)
476 {
477 itemp[icount] = atom[i].iat[0];
478 icount++;
479 }
480 if (jj2 == 1)
481 {
482 itemp[icount] = atom[i].iat[1];
483 icount++;
484 }
485 if (jj3 == 1)
486 {
487 itemp[icount] = atom[i].iat[2];
488 icount++;
489 }
490 if (jj3 == 1)
491 {
492 itemp[icount] = atom[i].iat[3];
493 icount++;
494 }
495 if (icount >= 2)
496 {
497 atom[itemp[0]].z += 0.5;
498 atom[itemp[1]].z -= 0.5;
499 }
500 }
501 }
502 }
503 for (i=1; i < natom; i++)
504 {
505 for (j=i+1; j <= natom; j++)
506 {
507 xtmp = atom[i].x - atom[j].x;
508 ytmp = atom[i].y - atom[j].y;
509 ztmp = atom[i].z - atom[j].z;
510 if ( (xtmp + ytmp + ztmp) < 0.1)
511 {
512 atom[i].x += 0.1;
513 atom[i].y += 0.1;
514 atom[i].z += 0.1;
515 atom[j].x -= 0.1;
516 atom[j].y -= 0.1;
517 atom[j].z -= 0.1;
518 }
519 }
520 }
521
522 if (has_Aromatic == TRUE)
523 mopaco(1,natom);
524 return TRUE;
525 }
526 // ==================================================
527 void read_schakal(int isel,int isubred)
528 {
529 int i, itype, newatom, iz;
530 float xtmp, ytmp, ztmp;
531 char dummy[30];
532 char atomchar[6];
533 char inputline[150];
534 FILE *infile;
535
536 infile = fopen_path(Openbox.path,Openbox.fname,"r");
537
538 if (infile == NULL)
539 {
540 message_alert("Error opening Schakal file","Schakal Setup");
541 return;
542 }
543 while ( FetchRecord(infile,inputline))
544 {
545 sscanf(inputline,"%s",dummy);
546 if (strcmp(dummy,"TITLE") == 0)
547 {
548 sscanf(inputline,"%s %s",dummy,Struct_Title);
549 } else if (strcmp(dummy,"ATOM") == 0)
550 {
551 sscanf(inputline,"%s %s %f %f %f",dummy,atomchar,&xtmp,&ytmp,&ztmp);
552 // strip off numbers
553 iz = strlen(atomchar);
554 for (i=0; i < iz; i++)
555 {
556 if (isdigit(atomchar[i]))
557 {
558 atomchar[i] = '\0';
559 break;
560 }
561 }
562 itype = 0;
563 newatom = make_atom(itype,xtmp,ytmp,ztmp,atomchar);
564 } else if (strcmp(dummy,"END") == 0)
565 {
566 break;
567 }
568 }
569 fclose(infile);
570 mopaco(1,natom);
571 //
572 }
573 // ==================================================
574 void write_schakal()
575 {
576 int i,j,lptest;
577 char atomchar[6];
578 FILE *wfile;
579
580 lptest = 0;
581 for( i = 1; i <= natom; i++ )
582 {
583 if( atom[i].mmx_type == 20 )
584 {
585 lptest = 1;
586 hdel( lptest );
587 break;
588 }
589 }
590
591 /* now write the schakal file */
592 if (files.append)
593 wfile = fopen_path(Savebox.path,Savebox.fname,"a");
594 else
595 wfile = fopen_path(Savebox.path,Savebox.fname,"w");
596
597 j = strlen(Struct_Title);
598 for (i=0; i < j; i++)
599 {
600 if (Struct_Title[i] == '\n')
601 {
602 Struct_Title[i] = '\0';
603 break;
604 }
605 }
606 fprintf(wfile,"TITLE %s\n",Struct_Title);
607 fprintf(wfile,"CELL \n");
608 for (i=1; i <= natom; i++)
609 {
610 sprintf(atomchar,"%s%d",Elements[atom[i].atomnum-1].symbol,i);
611
612 fprintf(wfile,"ATOM %-5s %12.7f%12.7f%12.7f\n",atomchar,atom[i].x,
613 atom[i].y, atom[i].z);
614 }
615 fprintf(wfile,"END \n");
616 fclose(wfile);
617 }