ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/freemol/trunk/src/mengine/src/read_sdf.c
Revision: 48
Committed: Thu Jul 31 16:34:52 2008 UTC (11 years, 8 months ago) by tjod
File size: 17096 byte(s)
Log Message:
Add modification notice to files modified for Freemol
per mengine/smi23d license requirements

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