ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/freemol/trunk/src/mengine/src/read_sdf.c
Revision: 27
Committed: Tue Jul 8 19:10:25 2008 UTC (13 years, 5 months ago) by tjod
File size: 16985 byte(s)
Log Message:
move mengine src to reflect new freemol directory structure

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 }
259 make_bond(ia1, ia2, ibond);
260 if (istereo == 1 && planar)
261 {
262 atom[ia2].z += 0.7;
263 if (atom[ia2].atomnum == 1)
264 atom[ia2].type = 60;
265 }
266 if (istereo == 6 && planar)
267 {
268 atom[ia2].z -= 0.7;
269 if (atom[ia2].atomnum == 1)
270 atom[ia2].type = 60;
271 }
272 }
273 // read to end of structure
274 while (FetchRecord(rfile,inputline))
275 {
276 if (strncasecmp(inputline,"$$$$",4) == 0)
277 {
278 return Ret_Val;
279 }
280 }
281 return Ret_Val;
282 }
283 /* =================================== */
284 int read_sdf(int istruct, int isubred)
285 {
286 int i, j, niatom, nibond, ia1, ia2, ibond, itype, newatom;
287 int ncount, ibotptr,junk,junk1,got_title,istereo,iz;
288 int jji, jj1, jj2, jj3, jj4;
289 int has_Aromatic, xPlus,yPlus,zPlus, planar;
290 int icount,itemp[4];
291 float xtmp, ytmp, ztmp, dz;
292 char c1[4],c2[4];
293 char atomchar[3];
294 char inputline[150];
295 FILE *infile;
296
297 infile = fopen_path(Openbox.path,Openbox.fname,"r");
298
299 if (infile == NULL)
300 {
301 message_alert("Error opening Concord file","Concord Setup");
302 return FALSE;
303 }
304
305 if (isubred == 0)
306 ibotptr = 0;
307 else
308 {
309 ibotptr = natom;
310 substr.istract[substr.icurstr] = TRUE;
311 }
312
313 if (istruct > 1)
314 {
315 ncount = 0;
316 while (FetchRecord(infile,inputline))
317 {
318 if (strncasecmp(inputline,"$$$$",4) == 0)
319 {
320 ncount++;
321 if (ncount+1 == istruct)
322 goto L_1;
323 }
324 }
325 }
326 // start file read here
327 L_1:
328 got_title = FALSE;
329 xPlus = yPlus = zPlus = FALSE;
330 FetchRecord(infile,inputline);
331 sscanf(inputline,"SDF %s",Struct_Title);
332 got_title = TRUE;
333 /* if (strlen(inputline) > 4)
334 {
335 iz = strlen(inputline);
336 if (iz > 60) iz = 59;
337 for (i=4; i < iz; i++)
338 Struct_Title[i-4] = inputline[i];
339 Struct_Title[i] = '\0';
340 got_title = TRUE;
341 } */
342 FetchRecord(infile,inputline); // blank line
343 FetchRecord(infile,inputline); // blank line
344
345 FetchRecord(infile,inputline); // natom and nbond
346 for (i=0; i <4; i++)
347 {
348 c1[i] = inputline[i];
349 c2[i] = inputline[i+3];
350 }
351 c1[3] = '\0';c2[3] = '\0';
352 niatom = atoi(c1);
353 nibond = atoi(c2);
354 if (niatom == 0)
355 return FALSE;
356
357 for (i=0; i < niatom; i++)
358 {
359 FetchRecord(infile,inputline);
360 sscanf(inputline,"%f %f %f %s %d %d",&xtmp, &ytmp, &ztmp, atomchar,&junk,&junk1);
361
362 itype = 0;
363 if (xtmp != 0.0) xPlus= TRUE;
364 if (ytmp != 0.0) yPlus= TRUE;
365 if (ztmp != 0.0) zPlus= TRUE;
366
367 iz = strlen(atomchar);
368 if (itype == 0 && (atomchar[0] == 'N' || atomchar[0] == 'O') && iz == 1) // nitro fix
369 {
370 if (atomchar[0] == 'N') itype = 8;
371 if (atomchar[0] == 'O') itype = 6;
372
373 if (junk1 != 0)
374 {
375 if (junk1 == 3 && itype == 8)
376 itype = 41; // N+
377 if (junk1 == 5 && itype == 6)
378 itype = 42; // O-
379 }
380 }
381
382 newatom = make_atom(itype,xtmp,ytmp,ztmp,atomchar);
383 if (newatom == -1)
384 {
385 printf("Atom %d %s not recognized structure skipped\n",i,atomchar);
386 return FALSE;
387 }
388 if (isubred == 1)
389 {
390 atom[newatom].substr[0] = 0;
391 atom[newatom].substr[0] |= (1L << substr.icurstr);
392 }
393 }
394 has_Aromatic = FALSE;
395
396 planar = TRUE;
397 if (xPlus && yPlus && zPlus) planar = FALSE;
398
399 for (i=0; i < nibond; i++)
400 {
401 FetchRecord(infile,inputline);
402 for (j=0; j <4; j++)
403 {
404 c1[j] = inputline[j];
405 c2[j] = inputline[j+3];
406 }
407 c1[3] = '\0';c2[3] = '\0';
408 ia1 = atoi(c1);
409 ia2 = atoi(c2);
410 istereo = 0;
411 sscanf(inputline,"%3d%3d%3d%3d",&junk, &junk1, &ibond,&istereo);
412 if (ibond >= 4)
413 {
414 ibond = 1;
415 has_Aromatic = TRUE;
416 }
417 make_bond(ia1+ibotptr, ia2+ibotptr, ibond);
418 if (istereo == 1 && planar)
419 {
420 atom[ia2+ibotptr].z += 0.7;
421 if (atom[ia2+ibotptr].atomnum == 1)
422 atom[ia2+ibotptr].type = 60;
423 }
424 if (istereo == 6 && planar)
425 {
426 atom[ia2+ibotptr].z -= 0.7;
427 if (atom[ia2+ibotptr].atomnum == 1)
428 atom[ia2+ibotptr].type = 60;
429 }
430 }
431 FetchRecord(infile,inputline); // M END line
432 FetchRecord(infile,inputline);
433 if (got_title == FALSE)
434 strcpy(Struct_Title,inputline);
435 fclose(infile);
436 ncount = strlen(Struct_Title);
437 for (i=0; i < ncount; i++)
438 {
439 if (Struct_Title[i] == '\n')
440 {
441 Struct_Title[i] = '\0';
442 break;
443 }
444 }
445 // search for quaternary centers - check if flat
446 dz = 0.0;
447 for (i=1; i <= natom; i++)
448 {
449 dz += atom[i].z;
450 jji = jj1 = jj2 = jj3 = jj4 = 0;
451 if (atom[i].type != 0)
452 {
453 for (j=0; j < 6; j++)
454 {
455 if (atom[i].iat[j] != 0)
456 jji++;
457 }
458 }
459 if (jji == 4)
460 {
461 if (atom[atom[i].iat[0]].z == 0.0 && atom[atom[i].iat[0]].z == 0.0 &&
462 atom[atom[i].iat[0]].z == 0.0 && atom[atom[i].iat[0]].z == 0.0) // flat center
463 {
464 for (j=0; j < 4; j++)
465 {
466 itemp[j] = 0;
467 if (atom[atom[i].iat[0]].iat[j] != 0) jj1++;
468 if (atom[atom[i].iat[1]].iat[j] != 0) jj2++;
469 if (atom[atom[i].iat[2]].iat[j] != 0) jj3++;
470 if (atom[atom[i].iat[3]].iat[j] != 0) jj4++;
471 }
472 icount = 0;
473 if (jj1 == 1)
474 {
475 itemp[icount] = atom[i].iat[0];
476 icount++;
477 }
478 if (jj2 == 1)
479 {
480 itemp[icount] = atom[i].iat[1];
481 icount++;
482 }
483 if (jj3 == 1)
484 {
485 itemp[icount] = atom[i].iat[2];
486 icount++;
487 }
488 if (jj3 == 1)
489 {
490 itemp[icount] = atom[i].iat[3];
491 icount++;
492 }
493 if (icount >= 2)
494 {
495 atom[itemp[0]].z += 0.5;
496 atom[itemp[1]].z -= 0.5;
497 }
498 }
499 }
500 }
501 for (i=1; i < natom; i++)
502 {
503 for (j=i+1; j <= natom; j++)
504 {
505 xtmp = atom[i].x - atom[j].x;
506 ytmp = atom[i].y - atom[j].y;
507 ztmp = atom[i].z - atom[j].z;
508 if ( (xtmp + ytmp + ztmp) < 0.1)
509 {
510 atom[i].x += 0.1;
511 atom[i].y += 0.1;
512 atom[i].z += 0.1;
513 atom[j].x -= 0.1;
514 atom[j].y -= 0.1;
515 atom[j].z -= 0.1;
516 }
517 }
518 }
519
520 if (has_Aromatic == TRUE)
521 mopaco(1,natom);
522 return TRUE;
523 }
524 // ==================================================
525 void read_schakal(int isel,int isubred)
526 {
527 int i, itype, newatom, iz;
528 float xtmp, ytmp, ztmp;
529 char dummy[30];
530 char atomchar[6];
531 char inputline[150];
532 FILE *infile;
533
534 infile = fopen_path(Openbox.path,Openbox.fname,"r");
535
536 if (infile == NULL)
537 {
538 message_alert("Error opening Schakal file","Schakal Setup");
539 return;
540 }
541 while ( FetchRecord(infile,inputline))
542 {
543 sscanf(inputline,"%s",dummy);
544 if (strcmp(dummy,"TITLE") == 0)
545 {
546 sscanf(inputline,"%s %s",dummy,Struct_Title);
547 } else if (strcmp(dummy,"ATOM") == 0)
548 {
549 sscanf(inputline,"%s %s %f %f %f",dummy,atomchar,&xtmp,&ytmp,&ztmp);
550 // strip off numbers
551 iz = strlen(atomchar);
552 for (i=0; i < iz; i++)
553 {
554 if (isdigit(atomchar[i]))
555 {
556 atomchar[i] = '\0';
557 break;
558 }
559 }
560 itype = 0;
561 newatom = make_atom(itype,xtmp,ytmp,ztmp,atomchar);
562 } else if (strcmp(dummy,"END") == 0)
563 {
564 break;
565 }
566 }
567 fclose(infile);
568 mopaco(1,natom);
569 //
570 }
571 // ==================================================
572 void write_schakal()
573 {
574 int i,j,lptest;
575 char atomchar[6];
576 FILE *wfile;
577
578 lptest = 0;
579 for( i = 1; i <= natom; i++ )
580 {
581 if( atom[i].mmx_type == 20 )
582 {
583 lptest = 1;
584 hdel( lptest );
585 break;
586 }
587 }
588
589 /* now write the schakal file */
590 if (files.append)
591 wfile = fopen_path(Savebox.path,Savebox.fname,"a");
592 else
593 wfile = fopen_path(Savebox.path,Savebox.fname,"w");
594
595 j = strlen(Struct_Title);
596 for (i=0; i < j; i++)
597 {
598 if (Struct_Title[i] == '\n')
599 {
600 Struct_Title[i] = '\0';
601 break;
602 }
603 }
604 fprintf(wfile,"TITLE %s\n",Struct_Title);
605 fprintf(wfile,"CELL \n");
606 for (i=1; i <= natom; i++)
607 {
608 sprintf(atomchar,"%s%d",Elements[atom[i].atomnum-1].symbol,i);
609
610 fprintf(wfile,"ATOM %-5s %12.7f%12.7f%12.7f\n",atomchar,atom[i].x,
611 atom[i].y, atom[i].z);
612 }
613 fprintf(wfile,"END \n");
614 fclose(wfile);
615 }