ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/freemol/trunk/smi23d/src/mengine/read.c
Revision: 3
Committed: Mon Jun 9 21:38:26 2008 UTC (11 years, 4 months ago) by tjod
File size: 21581 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 "field.h"
7 #include "pot.h"
8 #include "pdb.h"
9
10 #include <errno.h>
11
12 #define MAIN 0L
13 #define SUBSTRUCTURE 1L
14 #define MOPAC 0L
15 #define MNDO 1L
16
17 int FetchRecord(FILE *, char *);
18 FILE * fopen_path ( char * , char * , char * ) ;
19 int is_bond(int, int);
20 int get_bondorder(int,int);
21 void substr_scale(void);
22 void writefile(int);
23 void readpid(char *);
24 void readfile(int,int);
25 void check_numfile(int);
26 void inxyz(int);
27 void incsd(int,int);
28 void mopaco(int,int);
29 void find_termini(int *,int *);
30 void find_nucterm(int *,int *);
31 void find_cp(int *);
32 void type(void);
33 void hdel(int);
34 void hcoord(void);
35 void pcmfout(int);
36 void pcmfin(int,int);
37 void mmp22mod(int,int);
38 void mm32mod(int);
39 void mod2mmp2(int);
40 void readz(int, int);
41 void mdl2mod(int);
42 void mod2c3d(void);
43 void mod2tinker(void);
44 void c3d2mod(int,int *);
45 void tink2mod(int,int *);
46 void readz(int, int);
47 void read_arcfile(int, int);
48 void write_mopac(void);
49 void alc2mod(int,int *);
50 void mod2alc(void);
51 void mod2sybyl(void);
52 void mod2sybyl2(void);
53 void sybyl2mod(int ,int *);
54 void mod2mac(void);
55 void mac2mod(int, int);
56 void read_psgvbin(int);
57 void read_psgvbout(int);
58 void write_psgvbin(void);
59 void read_freq_log(int,int);
60 void read_games(int,int);
61 void read_hondo(int,int);
62 void initialize(void);
63 void write_mm2(void);
64 void write_mm3(void);
65 void mod2pdb(void);
66 void make_bond(int,int,int);
67 void resetsubstrmem(void);
68 void subconnect(int,int);
69 void mark_hbond(void);
70 void set_atomtypes(int);
71 void mm2gaus(void);
72 void mm2games(void);
73 void write_hondo(void);
74 int read_sdf(int,int);
75 void write_sdf(void);
76 void message_alert(char *, char *);
77 void write_eht(void);
78 void read_gausfchk(int);
79 void read_tm(int,int);
80 void mod2tm(void);
81 void read_adf(int,int);
82 void write_adf(void);
83 void read_smiles(int,int);
84 void write_smiles(void);
85 void build_coord(void);
86
87 EXTERN struct t_files {
88 int nfiles, append, batch, icurrent;
89 int istrselect;
90 } files;
91
92 int filetypes[10000];
93
94 int oh2, oh3, oh6;
95 void clean_string(char *);
96 // ========================================
97 void clean_string(char *astring)
98 {
99 int iz,i;
100 iz = strlen(astring);
101 for (i=0; i < iz; i++)
102 {
103 if (isprint(astring[i]) == 0)
104 astring[i] = ' ';
105 }
106 for (i = iz -1; i >= 0; i--)
107 {
108 if (astring[i] != ' ')
109 {
110 astring[i+1] = '\0';
111 break;
112 }
113 }
114 }
115 /*======================================================================*/
116 void * malloc_filename ( char *path , char *name )
117 {
118 int ix, iz;
119 char *filename;
120
121 ix = 0;
122 if (path != NULL)
123 ix = strlen(path);
124 iz = strlen(name);
125 filename = malloc (ix + 1 + iz + 1) ;
126 if ( ix == 0 )
127 {
128 strcpy (filename,name);
129 }
130 else
131 {
132 sprintf (filename,"%s\\%s",path,name) ;
133 }
134 return ( (void *) filename ) ;
135 }
136 /*======================================================================*/
137 FILE * fopen_path ( char *path , char *name , char *mode )
138 {
139 char *filename = malloc_filename(path,name) ;
140 FILE *rval = fopen (filename,mode) ;
141 if ( rval == NULL )
142 {
143 fprintf (pcmoutfile,"open of --%s-- failed\n",filename) ;
144 fprintf(pcmoutfile,"Error %d %s\n",errno,strerror(errno));
145 fclose(pcmoutfile);
146 exit(0);
147 }
148 free (filename) ;
149 return ( rval ) ;
150 }
151
152 /* ============================================= */
153 void readfile(int isel, int isubred)
154 {
155 int iftype, i;
156 int oldfield;
157
158 /* need to determine number of structures in file
159 which structure the user wants to read and then
160 call the appropriate read routine. Check to see if
161 we already have structure and erase if necessary */
162
163
164 if (isubred == 0)
165 initialize();
166
167
168 iftype = Openbox.ftype;
169 if (isubred == 1)
170 {
171 for (i=0; i < MAXSS; i++)
172 {
173 if ( !substr.istract[i] )
174 break;
175 }
176 substr.icurstr = i;
177 substr.istract[i] = TRUE;
178 }
179
180 switch(Openbox.ftype)
181 {
182 case FTYPE_PCM:
183 oldfield = field.type;
184 if (field.type != MMX)
185 field.type = MMX;
186 pcmfin(isel,isubred);
187 field.type = oldfield;
188 break;
189 case FTYPE_SDF: // need to check number of structures in file
190 oldfield = field.type;
191 if (field.type != MMX)
192 field.type = MMX;
193 i = read_sdf(isel,isubred);
194 field.type = oldfield;
195 break;
196 }
197 /* move substructure to right of main structures */
198 if (isubred == 1)
199 substr_scale();
200
201 // convert to correct types for current force field
202 type();
203 for(i=1; i <= natom; i++)
204 {
205 if (atom[i].mmx_type == 5)
206 {
207 flags.noh = TRUE;
208 break;
209 }
210 }
211 substr.icurstr = -1;
212 return;
213 }
214 /* ================================================== */
215 void writefile(int isubred)
216 {
217 int oldfield;
218
219 oldfield = field.type;
220 if (Savebox.ftype == FTYPE_MM3 && field.type == MMX)
221 {
222 // convert to mm3 types
223 set_atomtypes(MM3);
224 field.type = MM3;
225 type();
226 } else if (field.type == MM3)
227 {
228 // convert to MMX types
229 set_atomtypes(MMX);
230 field.type = MMX;
231 type();
232 }
233 if (isubred == MAIN)
234 {
235 switch(Savebox.ftype)
236 {
237 case FTYPE_PCM:
238 pcmfout(1);
239 break;
240 case FTYPE_SDF:
241 write_sdf();
242 break;
243 }
244 }
245 }
246 /* ------------------------------------------------- */
247 void check_numfile(int ftype)
248 {
249 int igau, header, model;
250 FILE *infile;
251 char inputline[160];
252
253 header = FALSE;
254 model = FALSE;
255 infile = fopen_path(Openbox.path,Openbox.fname,"r");
256 if (infile == NULL)
257 {
258 message_alert("Error trying to check number of structures in file. Probably a bad filename.","Check Numfile");
259 fprintf(pcmoutfile,"Error in check numfile:: %s %s\n",Openbox.path, Openbox.fname);
260 files.nfiles = 0;
261 return;
262 }
263 igau = 0;
264 while ( fgets(inputline,159,infile) != NULL )
265 {
266 if (strncasecmp(inputline,"{PCM",4) == 0)
267 {
268 filetypes[files.nfiles] = FTYPE_PCM;
269 files.nfiles++;
270 }else if (strncasecmp(inputline,"SDF ",4) == 0)
271 {
272 filetypes[files.nfiles] = FTYPE_SDF;
273 files.nfiles++;
274 }
275 }
276 fclose(infile);
277 }
278 /* ------------------------------------*/
279 int get_bondorder(int ia,int ib)
280 {
281 int i;
282 for (i=0; i < MAXIAT; i++)
283 {
284 if (atom[ia].iat[i] == ib)
285 return atom[ia].bo[i];
286 }
287 return 0;
288 }
289 /* ------------------------------------*/
290 int is_bond(int ia, int ib)
291 {
292 int i;
293 for (i=0; i < MAXIAT; i++)
294 {
295 if (atom[ia].iat[i] == ib)
296 return TRUE;
297 }
298 return FALSE;
299 }
300 /* ------------------------------------*/
301 void mopaco(int istart,int ifinish)
302 {
303 int i, ibi, ibj, it, j, jt, kk;
304 float disij, dx, dy, dz, r, xii, yii, zii,rbdis[4];
305 static float radii[110] =
306 { 0.77, 0.77, 0.77, 0.77, 0.50, 0.66, 0.66, 0.70, 0.70, 0.70,
307 0.64, 0.99, 1.14, 1.33, 1.04, 1.04, 1.04, 1.04, 1.17, 0.01,
308 0.50, 0.77, 0.50, 0.50, 1.10, 0.88, 0.88, 0.50, 0.77, 0.77,
309 1.22, 1.40, 1.20, 1.17, 1.37, 0.50, 0.70, 1.04, 1.17, 0.77,
310 0.70, 0.66, 0.77, 0.70, 0.50, 0.66, 1.10, 0.77, 0.77, 0.77,
311 0.77, 0.66, 0.70, 1.33, 0.66, 0.77, 0.77, 0.70, 0.01, 0.50,
312 0.77, 0.77, 0.77, 0.70, 0.66, 0.70, 1.10, 0.66, 0.70, 0.50,
313 0.77, 0.77, 0.77, 0.77, 0.50, 0.66, 0.66, 0.70, 0.70, 0.70,
314 0.77, 0.77, 0.77, 0.77, 0.50, 0.66, 0.66, 0.70, 0.70, 0.70,
315 0.77, 0.77, 0.77, 0.77, 0.50, 0.66, 0.66, 0.70, 0.70, 0.70,
316 0.77, 0.77, 0.77, 0.77, 0.50, 0.66, 0.66, 0.70, 0.70, 0.70,};
317
318 /* SUBROUTINE SENDS BACK THE ACTUAL NUMBER OF ATOMS, AND
319 * DETERMINES THE BOND CONNECTIONS and bond orders TO ESTABLISH THE attached
320 * ATOM LIST. ATOM TYPES HAVE BEEN CHANGED TO MM2 TYPES. THOSE THAT
321 * AREN'T IN MM2 ARE GIVEN #S = AT # + 60
322 * BOND ORDERS AND THE IATMS AND IBNDS LISTS ARE GENERATED
323 * NATOMS IS THE # OF ENTRIES IN THE Z MATRIX */
324
325 rbdis[0] = 1.42;
326 rbdis[1] = 1.25;
327 rbdis[2] = 1.37;
328 rbdis[3] = 1.25;
329
330 /* find all single bonds based on sum of covalent radii */
331 for( i = istart; i <= (ifinish - 1L); i++ )
332 {
333 it = atom[i].mmx_type;
334 for( j = i + 1L; j <= ifinish; j++ )
335 {
336 if (!is_bond(i,j))
337 {
338 jt = atom[j].mmx_type;
339 dx = atom[i].x - atom[j].x;
340 dy = atom[i].y - atom[j].y;
341 dz = atom[i].z - atom[j].z;
342 r = sqrt( dx*dx + dy*dy + dz*dz );
343 if( it >= 80L && jt >= 80L )
344 {
345 if( r <= (atom[i].radius + atom[j].radius) )
346 make_bond( i, j, 1 );
347 }else if( it >= 80L && jt < 80L )
348 {
349 if( r <= (atom[i].radius + radii[atom[j].mmx_type-1]) )
350 make_bond( i, j, 1);
351
352 }else if( it < 80L && jt >= 80L )
353 {
354 if( r <= (radii[atom[i].mmx_type-1] + atom[j].radius + .1) )
355 make_bond( i, j, 1 );
356 }else if( r <= (radii[atom[i].mmx_type-1] + radii[atom[j].mmx_type-1]+ .1) )
357 {
358 if( i < j )
359 make_bond( i, j, 1 );
360 }
361 }
362 }
363 }
364
365 /* find all double and triple bonds - search only types for
366 * carbon, oxygen, nitrogen, sulfur, phosphorous */
367 for( i = istart; i <= (ifinish - 1); i++ )
368 {
369 it = atom[i].atomnum;
370 if( it == 6 || it == 7 || it == 8 || it == 15 || it == 16 )
371 {
372 for( j = 0; j < MAXIAT; j++ )
373 {
374 if( atom[i].iat[j] != 0 )
375 {
376 jt = atom[atom[i].iat[j]].atomnum;
377 if( jt == 6 || jt == 7 || jt == 8 || jt == 15 || jt == 16 )
378 {
379 xii = atom[i].x - atom[atom[i].iat[j]].x;
380 yii = atom[i].y - atom[atom[i].iat[j]].y;
381 zii = atom[i].z - atom[atom[i].iat[j]].z;
382 if( fabs( xii ) <= 2.3 )
383 {
384 if( fabs( yii ) <= 2.3 )
385 {
386 if( fabs( zii ) <= 2.3 )
387 {
388 disij = sqrt( xii*xii + yii*yii + zii*zii );
389 /* c=o bond */
390 if( ((it == 8 && jt == 6) || (it == 6 && jt == 8)) && disij <= rbdis[3] + 0.05 )
391 {
392 if( i < atom[i].iat[j] )
393 make_bond( i, atom[i].iat[j],1 );
394 /* alkyne */
395 }else if( (it == 6 && jt == 6) && disij <= rbdis[1] )
396 {
397 if( i < atom[i].iat[j] )
398 make_bond( i, atom[i].iat[j],2 );
399 /* c=n */
400 }else if( ((it == 7 && jt == 6) || (it == 6 && jt == 7))
401 && disij <= rbdis[2] )
402 {
403 ibi = 0;
404 for( kk = 0; kk < MAXIAT; kk++ )
405 {
406 if( atom[i].bo[kk] != 9 )
407 ibi = ibi + atom[i].bo[kk];
408 }
409 ibj = 0;
410 for( kk = 0; kk < MAXIAT; kk++ )
411 {
412 if( atom[atom[i].iat[j]].bo[kk] != 9 )
413 ibj = ibj + atom[atom[i].iat[j]].bo[kk];
414 }
415 if( it == 7 )
416 {
417 if( ibi <= 2 && ibj <= 3 )
418 {
419 if( i < atom[i].iat[j] )
420 make_bond( i, atom[i].iat[j], 1 );
421 break;
422 }
423 }else if( jt == 7 )
424 {
425 if( ibj <= 2 && ibi <= 3 )
426 {
427 if( i < atom[i].iat[j] )
428 make_bond( i, atom[i].iat[j], 1 );
429 break;
430 }
431 }
432 /*nitrile */ }else if( ((it == 7 && jt == 6) || (it == 6 && jt == 7)) && disij <= 1.18 )
433 {
434 if( i < atom[i].iat[j] )
435 make_bond( i, atom[i].iat[j], 2 );
436 /* N=O */ }else if ( ((it == 7 && jt == 8) || (it == 8 && jt == 7)) && disij <= 1.26)
437 {
438 ibi = 0;
439 if (it == 7)
440 {
441 for (kk = 0; kk < MAXIAT; kk++)
442 {
443 if (atom[i].bo[kk] != 9)
444 ibi += atom[i].bo[kk];
445 }
446 } else
447 {
448 for (kk = 0; kk < MAXIAT; kk++)
449 {
450 if (atom[atom[i].iat[j]].bo[kk] != 9)
451 ibi += atom[atom[i].iat[j]].bo[kk];
452 }
453 }
454 if (ibi <= 3)
455 {
456 if ( i < atom[i].iat[j])
457 make_bond(i,atom[i].iat[j],1);
458 break;
459 }
460 /* c=c bond */ }else if( (it == 6 && jt == 6) && disij <= rbdis[0] )
461 {
462 ibi = 0;
463 for( kk = 0; kk < MAXIAT; kk++ )
464 {
465 if( atom[i].bo[kk] != 9 )
466 ibi = ibi + atom[i].bo[kk];
467 }
468 ibj = 0;
469 for( kk = 0; kk < MAXIAT; kk++ )
470 {
471 if( atom[atom[i].iat[j]].bo[kk] != 9 )
472 ibj = ibj + atom[atom[i].iat[j]].bo[kk];
473 }
474 if( ibi <= 3 && ibj <= 3 )
475 {
476 if( i < atom[i].iat[j] )
477 make_bond( i, atom[i].iat[j], 1 );
478 }
479 }
480 }
481 }
482 }
483 }
484 }
485 }
486 }
487 }
488 type();
489 return;
490 }
491 void find_nucterm(int *i5p,int *i3p)
492 {
493 int i, icount;
494 long int mask1,mask2;
495 // assume find first 3p end and second 5p
496 mask1 = 1L << P5; /* 5 prime terminus */
497 mask2 = 1L << P3; /* 3 prime terminus */
498
499 icount = 0;
500 for (i=1; i <= natom; i++)
501 {
502 if (atom[i].flags & mask1)
503 {
504 icount++;
505 if (icount == 2)
506 {
507 *i5p = i;
508 break;
509 }
510 }
511 }
512 for (i=1; i <= natom; i++)
513 {
514 if (atom[i].flags & mask2)
515 {
516 *i3p = i;
517 break;
518 }
519 }
520 }
521 void find_termini(int *iNT,int *iCT)
522 {
523 long int i, j, mask1, mask2;
524
525 mask1 = 1L << NTERM; /* N terminus */
526 mask2 = 1L << CNTERM; /* C terminus */
527
528 for (i = 1; i < natom; i++)
529 {
530 if (atom[i].flags & mask2)
531 {
532 *iCT = i;
533 for (j = i+1; j < natom; j++)
534 {
535 if (atom[j].flags & mask1)
536 {
537 *iNT = j;
538 break;
539 }
540 }
541 break;
542 }
543 }
544 return;
545 }
546 void find_cp(int *iCP)
547 {
548 long int i, mask1;
549
550 mask1 = 1L << DUMMY; /* connect point */
551
552 for (i = 1; i < natom; i++)
553 {
554 if (atom[i].flags & mask1)
555 {
556 *iCP = i;
557 break;
558 }
559 }
560 return;
561 }
562
563 void substr_scale()
564 {
565 int isub,ii;
566 float xmain_min,xmain_max, sub_xmin, xdif, sub_center, main_center;
567 float ymain_min,ymain_max, sub_ymax, sub_ymin, ydif;
568
569 xmain_min = 1000.0F;
570 ymain_min = 1000.0F;
571 sub_xmin = 1000.0F;
572 sub_ymin = 1000.0F;
573
574 xmain_max = -1000.0F;
575 ymain_max = -1000.0F;
576 sub_ymax = -1000.0F;
577
578 isub = (1L << substr.icurstr);
579
580 for (ii = 1; ii <= natom; ii++)
581 {
582 if (atom[ii].substr[0] == isub)
583 {
584 if ( atom[ii].x < sub_xmin )
585 sub_xmin = atom[ii].x;
586 if ( atom[ii].y < sub_ymin )
587 sub_ymin = atom[ii].y;
588 if ( atom[ii].y > sub_ymax )
589 sub_ymax = atom[ii].y;
590 } else
591 {
592 if ( atom[ii].x < xmain_min )
593 xmain_min = atom[ii].x;
594 if (atom[ii].x > xmain_max)
595 xmain_max = atom[ii].x;
596 if ( atom[ii].y < ymain_min )
597 ymain_min = atom[ii].y;
598 if (atom[ii].y > ymain_max)
599 ymain_max = atom[ii].y;
600 }
601 }
602
603 if (xmain_max > sub_xmin)
604 {
605 xdif = xmain_max - sub_xmin;
606 for (ii = 1; ii <= natom; ii++)
607 {
608 if (atom[ii].substr[0] == isub)
609 atom[ii].x += xdif;
610 }
611 }
612
613 main_center = ((ymain_max - ymain_min)/2.0F) + ymain_min;
614 sub_center = ((sub_ymax - sub_ymin)/2.0F) + sub_ymin;
615
616 ydif = main_center - sub_center;
617
618 for (ii = 1; ii <= natom; ii++)
619 {
620 if (atom[ii].substr[0] == isub)
621 atom[ii].y += ydif;
622 }
623 }
624 /* ============================================= */
625 int FetchRecord(FILE *fp, char *buffer)
626 {
627 int ch;
628 char *ptr;
629
630 ptr = buffer;
631 do
632 {
633 ch = getc(fp);
634 if (ch >= ' ')
635 *ptr++ = ch;
636 else if (ch == '\n')
637 {
638 *ptr = 0;
639 return TRUE;
640 } else if (ch == '\r')
641 {
642 ch = getc(fp);
643 if (ch != '\n')
644 ungetc(ch,fp);
645 *ptr = 0;
646 return TRUE;
647 } else if (ch == EOF)
648 {
649 *ptr = 0;
650 return (ptr != buffer);
651 } else
652 *ptr++ = ch;
653 } while (ptr < buffer+255);
654
655 do
656 {
657 ch = getc(fp);
658 } while (ch !='\n' && ch != '\r' && ch != EOF);
659 if (ch == '\r')
660 {
661 ch = getc(fp);
662 if (ch != '\n')
663 ungetc(ch,fp);
664 }
665 *ptr = 0;
666 return TRUE;
667 }
668
669
670
671