ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/freemol/trunk/smi23d/src/mengine/read_sdf.c
Revision: 3
Committed: Mon Jun 9 21:38:26 2008 UTC (11 years, 4 months ago) by tjod
File size: 16878 byte(s)
Log Message:
test

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