ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/freemol/trunk/smi23d/src/mengine/read_sdf.c
Revision: 15
Committed: Thu Jun 19 00:26:16 2008 UTC (11 years, 9 months ago) by tjod
File size: 16985 byte(s)
Log Message:
Create new main program, mengine.c to replace gmmx.c
simplifying command line parsing (none).
Make mengine read stdin (sdf) and write stdout (sdf) instead
of files.  Makefile reflects change from gmmx to mengine.

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 }