ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/freemol/trunk/src/mengine/src/mengine.c
Revision: 85
Committed: Mon Jan 12 15:05:44 2009 UTC (10 years, 9 months ago) by gilbertke
File size: 19338 byte(s)
Log Message:
more cleanup and code removal
Line File contents
1 /* NOTICE: this source code file has been modified for use with FreeMOL */
2
3 #define notKEG_DEBUG
4
5 #define EXTERN
6 #include "pcwin.h"
7 #include "pcmod.h"
8 #include "energies.h"
9 #include "pot.h"
10 #include "angles.h"
11 #include "attached.h"
12 #include "rings.h"
13 #include "torsions.h"
14 #include "nonbond.h"
15 #include "bonds_ff.h"
16 #include "derivs.h"
17 #include "hess.h"
18 #include "field.h"
19 #include "atom_k.h"
20 #include "cutoffs.h"
21 #include "solv.h"
22 #include "fix.h"
23 #include "vibrate.h"
24 #include "dipmom.h"
25
26 #include <stdio.h>
27 #include <stdlib.h>
28 #include <unistd.h>
29 #include <string.h>
30 //#include <malloc.h>
31
32 struct t_minim_values {
33 int iprint, ndc, nconst;
34 float dielc;
35 } minim_values;
36
37 struct t_minim_control {
38 int type, method, field, added_const;
39 char added_path[256],added_name[256];
40 } minim_control;
41 EXTERN struct t_files {
42 int nfiles, append, batch, icurrent, ibatno;
43 } files;
44 char Savename[80];
45
46 /*
47 void initialize_pcmodel(char*);
48 void mmp22mod(int,int);
49 void pcmfin(int, int);
50 void initialize(void);
51 void mmxsub(int);
52 void check_numfile(int, char *);
53 void mac2mod(int, int);
54 void bbchk(void);
55 void type(void);
56 void initialize_gmmx(void);
57 void read_gmmxinp(char*, char*, char*, FILE*, int);
58 void run_gmmx(void);
59 void search_rings(int);
60 */
61
62
63 void usage(void) {
64 printf("\nUsage: mengine [-hvdxi] reads sdf on stdin, writes on stdout\n\n"
65 "-h\tThis help page\n"
66 "-v\tVerbose output\n"
67 "-a\tAdd hydrogen atoms; delete existing ones\n"
68 "-d\tEvaluate dipole moment\n"
69 "-x\tEvaluate XlogP\n"
70 "-i\tEvaluate vibrational data\n\n");
71 }
72
73 /* ==================================================== */
74 int main(int argc, char *argv[])
75 {
76 int i;
77 FILE *infile;
78 char *line = NULL;
79
80 char *infilename = NULL;
81 char *outfilename = NULL;
82 char *mmff94filename = NULL;
83 char *mmxfilename = NULL;
84 char *logfilename = NULL;
85
86 int c;
87 int flags = 0;
88
89 while ((c = getopt(argc, argv, "e:o:p:c:havdxi")) != -1) {
90 switch(c) {
91 case 'a':
92 flags = flags | DO_ADDH;
93 break;
94 case 'd':
95 flags = flags | DO_DIPOLE;
96 break;
97 case 'x':
98 flags = flags | DO_XLOGP;
99 break;
100 case 'i':
101 flags = flags | DO_VIBRATION;
102 break;
103 case 'v':
104 VERBOSE = 1;
105 break;
106 case 'h':
107 usage();
108 exit(-1);
109 case 'o':
110 outfilename = strdup(optarg);
111 break;
112 #if 0
113 // EXTERNAL PARAM FILES STILL NOT WORKING
114 // SOON THEY WILL
115 case 'p':
116 mmff94filename = strdup(optarg);
117 break;
118 case 'c':
119 mmxfilename = strdup(optarg);
120 break;
121 #endif
122 case 'e':
123 logfilename = strdup(optarg);
124 break;
125 case '?':
126 if (isprint (optopt))
127 fprintf (stderr, "Unknown option `-%c'.\n", optopt);
128 else
129 fprintf (stderr,
130 "Unknown option character `\\x%x'.\n",
131 optopt);
132 return 1;
133 default:
134 abort ();
135 }
136 }
137
138 argc -= optind;
139 argv += optind;
140
141 if (argc == 1)
142 infilename = strdup(argv[0]);
143
144 #ifdef KEG_DEBUG
145 if(!infilename)
146 infilename = strdup("caffeine.sdf");
147 if(!outfilename)
148 outfilename = strdup("output.sdf");
149 #endif
150
151 /* default behavior is to use the hardcoded parameter files */
152
153 if(!mmff94filename) mmff94filename = strdup("|mmff94"); /* hard-coded builtins */
154 if(!mmxfilename) mmxfilename = strdup("|mmxconst"); /* hard-coded builtins */
155
156 /* if no input file provide, then use standard input */
157
158 if(infilename) {
159 infile = fopen(infilename,"rb");
160 } else {
161 infile = stdin;
162 }
163
164 /* if no output file provided, then use standard output */
165
166 if(outfilename) {
167 pcmoutfile = fopen(outfilename,"wb");
168 } else {
169 pcmoutfile = stdout;
170 }
171
172 /* if no log file provided, then use standard error */
173
174 if(logfilename) {
175 pcmlogfile = fopen(logfilename,"wb");
176 } else {
177 pcmlogfile = stderr;
178 }
179
180 /* initialize the system using parameter file */
181
182 initialize_pcmodel(mmxfilename);
183 initialize();
184 minim_values.iprint = 0; //FALSE;
185 #define EXTERN extern
186
187 #include "pcwin.h"
188 #include "pcmod.h"
189
190 #include "angles.h"
191 #include "torsions.h"
192 #include "rings.h"
193
194
195 int isbond(int, int);
196 int is_ring31(int);
197 int is_ring41(int);
198 int is_ring42(int,int);
199 int is_ring43(int, int, int);
200 int is_ring44(int *);
201 int is_ring51(int);
202 int is_ring52(int, int);
203 int is_ring53(int, int, int);
204 int is_ring54(int, int, int, int);
205 int is_ring55(int *);
206 int is_ring61(int);
207 int is_ring62(int, int);
208 int is_ring63(int,int,int);
209 int is_ring66(int *);
210 void ksort(int, int *);
211 void message_alert(char *, char *);
212 int is_cyclo6(int, int *);
213
214 int is_ring31(int ia)
215 {
216 int i,j;
217
218 for (i=0; i < rings.nring3; i++)
219 {
220 for (j=0; j < 3; j++)
221 {
222 if (ia == rings.r13[i][j])
223 {
224 return(TRUE);
225 }
226 }
227 }
228 return FALSE;
229 }
230 /* ============================ */
231 int is_ring41(int ia)
232 {
233 int i,j, found;
234
235 found = FALSE;
236 for (i=0; i < rings.nring4; i++)
237 {
238 for (j=0; j < 4; j++)
239 {
240 if (ia == rings.r14[i][j])
241 {
242 found = TRUE;
243 return (found);
244 }
245 }
246 }
247 return (found);
248 }
249 /* -------------------------------------------------------- */
250 int is_ring42(int ia, int ib)
251 {
252 int i,j, k, itmp, jtmp;
253 if (ia < ib )
254 {
255 itmp = ia;
256 jtmp = ib;
257 }else
258 {
259 itmp = ib;
260 jtmp = ia;
261 }
262 for (i=0; i < rings.nring4; i++)
263 {
264 for (j = 0; j < 3; j++)
265 {
266 if ( itmp == rings.r14[i][j])
267 {
268 for (k=j+1; k < 4; k++)
269 {
270 if (jtmp == rings.r14[i][k])
271 return(TRUE);
272 }
273 }
274 }
275 }
276 return(FALSE);
277 }
278 /* -------------------------------------------------------- */
279 int is_ring43(int ia, int ib, int ic)
280 {
281 int i, found;
282 int array[5];
283
284 found = FALSE;
285 array[0] = ia;
286 array[1] = ib;
287 array[2] = ic;
288 ksort(3,array);
289 for (i=0; i < rings.nring4; i++)
290 {
291 if ( array[0] == rings.r14[i][0] && rings.r14[i][1] == array[1] && rings.r14[i][2] == array[2])
292 {
293 found = TRUE;
294 return(found);
295 }
296 if ( rings.r14[i][0] == array[0] && rings.r14[i][1] == array[1] && rings.r14[i][3] == array[2])
297 {
298 found = TRUE;
299 return(found);
300 }
301 if ( rings.r14[i][0] == array[0] && rings.r14[i][2] == array[1] && rings.r14[i][3] == array[2])
302 {
303 found = TRUE;
304 return(found);
305 }
306 if ( rings.r14[i][1] == array[0] && rings.r14[i][2] == array[1] && rings.r14[i][3] == array[2])
307 {
308 found = TRUE;
309 return(found);
310 }
311 }
312 return(found);
313 }
314 /* -------------------------------------------------------- */
315 int is_ring44(int *array)
316 {
317 int i;
318 for (i=0; i < rings.nring4; i++)
319 {
320 if (array[0] == rings.r14[i][0] && array[1] == rings.r14[i][1] &&
321 array[2] == rings.r14[i][2] && array[3] == rings.r14[i][3])
322 return(TRUE);
323 }
324 return FALSE;
325 }
326 /* -------------------------------------------------------- */
327 int is_ring51(int ia)
328 {
329 int i,j;
330
331 for(i=0; i < rings.nring5; i++)
332 {
333 for (j=0; j < 5; j++)
334 {
335 if (rings.r15[i][j] == ia)
336 {
337 return(TRUE);
338 }
339 }
340 }
341 return(FALSE);
342 }
343 /* -------------------------------------------------------- */
344 int is_ring52(int ia, int ib)
345 {
346 int i, j,k, itmp, jtmp;
347 if (ia < ib)
348 {
349 itmp = ia;
350 jtmp = ib;
351 }else
352 {
353 itmp = ib;
354 jtmp = ia;
355 }
356 for (i=0; i < rings.nring5; i++)
357 {
358 for (j=0; j < 4; j++)
359 {
360 if (itmp == rings.r15[i][j])
361 {
362 for (k=j+1; k < 5; k++)
363 {
364 if (jtmp == rings.r15[i][k])
365 return(TRUE);
366 }
367 }
368 }
369 }
370 return(FALSE);
371 }
372 /* -------------------------------------------------------- */
373 int is_ring53(int ia, int ib, int ic)
374 {
375 int i, j,k,l;
376 int array[5];
377
378 array[0] = ia;
379 array[1] = ib;
380 array[2] = ic;
381 ksort(3,array);
382 for (i=0; i < rings.nring5; i++)
383 {
384 for (j=0; j < 3; j++)
385 {
386 if (array[0] == rings.r15[i][j])
387 {
388 for (k=j+1; k < 4; k++)
389 {
390 if (array[1] == rings.r15[i][k])
391 {
392 for (l=j+1; l < 5; l++)
393 {
394 if (array[2] == rings.r15[i][l])
395 return(TRUE);
396 }
397 }
398 }
399 }
400 }
401 }
402 return(FALSE);
403 }
404 /* -------------------------------------------------------- */
405 int is_ring54(int ia, int ib, int ic, int id)
406 {
407 int i;
408 int array[5];
409 array[0] = ia;
410 array[1] = ib;
411 array[2] = ic;
412 array[3] = id;
413 ksort(4,array);
414 for (i=0; i < rings.nring5; i++)
415 {
416 if (array[0] == rings.r15[i][0])
417 {
418 if (array[1] == rings.r15[i][1])
419 {
420 if (array[2] == rings.r15[i][2])
421 {
422 if (array[3] == rings.r15[i][3])
423 return(TRUE);
424 else if (array[3] == rings.r15[i][4])
425 return(TRUE);
426 else
427 return(FALSE);
428 } else if (array[2] == rings.r15[i][3] && array[3] == rings.r15[i][4])
429 return(TRUE);
430 } else if (array[1] == rings.r15[i][2] &&
431 array[2] == rings.r15[i][3] && array[3] == rings.r15[i][4])
432 return(TRUE);
433
434 } else if (array[0] == rings.r15[i][1] && array[1] == rings.r15[i][2] &&
435 array[2] == rings.r15[i][3] && array[3] == rings.r15[i][4])
436 return(TRUE);
437 }
438 return (FALSE);
439 }
440 /* -------------------------------------------------------- */
441 int is_ring55(int *array)
442 {
443 int i;
444 for(i=0; i < rings.nring5; i++)
445 {
446 if ( array[0] == rings.r15[i][0] && array[1] == rings.r15[i][1] && array[2] == rings.r15[i][2] &&
447 array[3] == rings.r15[i][3] && array[4] == rings.r15[i][4])
448 return(TRUE);
449 }
450 return (FALSE);
451 }
452 /* -------------------------------------------------------- */
453 int is_ring61(int ia)
454 {
455 int i, j;
456 for(i=0; i < rings.nring6; i++)
457 {
458 for(j=0; j < 6; j++)
459 {
460 if (ia == rings.r16[i][j])
461 return (TRUE);
462 }
463 }
464 return (FALSE);
465 }
466 /* -------------------------------------------------------- */
467 int is_ring62(int ia, int ib)
468 {
469 int i, j,k, itmp, jtmp;
470 if (ia < ib)
471 {
472 itmp = ia;
473 jtmp = ib;
474 }else
475 {
476 itmp = ib;
477 jtmp = ia;
478 }
479 for (i=0; i < rings.nring6; i++)
480 {
481 for (j=0; j < 5; j++)
482 {
483 if (itmp == rings.r16[i][j])
484 {
485 for (k=j+1; k < 6; k++)
486 {
487 if (jtmp == rings.r16[i][k])
488 return(TRUE);
489 }
490 }
491 }
492 }
493 return(FALSE);
494 }
495 /* -------------------------------------------------------- */
496 int is_ring63(int ia, int ib, int ic)
497 {
498 int i, j, k, l;
499 int array[5];
500
501 array[0] = ia;
502 array[1] = ib;
503 array[2] = ic;
504 ksort(3,array);
505 for (i=0; i < rings.nring6; i++)
506 {
507 for (j=0; j < 4; j++)
508 {
509 if (array[0] == rings.r16[i][j])
510 {
511 for (k=j+1; k < 5; k++)
512 {
513 if (array[1] == rings.r16[i][k])
514 {
515 for (l=k+1; l < 6; l++)
516 {
517 if (array[2] == rings.r16[i][l])
518 return(TRUE);
519 }
520 }
521 }
522 }
523 }
524 }
525 return(FALSE);
526 }
527 /* -------------------------------------------------------- */
528 int is_ring66(int *array)
529 {
530 int i;
531 for(i=0; i < rings.nring6; i++)
532 {
533 if ( array[0] == rings.r16[i][0] && array[1] == rings.r16[i][1] && array[2] == rings.r16[i][2] &&
534 array[3] == rings.r16[i][3] && array[4] == rings.r16[i][4] && array[5] == rings.r16[i][5])
535 return(TRUE);
536 }
537 return (FALSE);
538 }
539 /* -------------------------------------------------------- */
540
541 void get_rings()
542 {
543 int i,j, add_ring, k, l;
544 int jatm,katm;
545 int array[6];
546
547 rings.nring3 = 0;
548 rings.nring4 = 0;
549 rings.nring5 = 0;
550 rings.nring6 = 0;
551
552 // get three membered rings
553 for (i=0; i < angles.nang; i++)
554 {
555 if (isbond(angles.i13[i][0], angles.i13[i][2]))
556 {
557 array[0] = angles.i13[i][0];
558 array[1] = angles.i13[i][1];
559 array[2] = angles.i13[i][2];
560 ksort(3, array);
561 add_ring = TRUE;
562 for (j = 0; j < rings.nring3; j++)
563 {
564 if (array[0] == rings.r13[j][0] && array[1] == rings.r13[j][1]
565 && array[2] == rings.r13[j][2])
566 {
567 add_ring = FALSE;
568 break;
569 }
570 }
571 if (add_ring == TRUE)
572 {
573 rings.r13[rings.nring3][0] = array[0];
574 rings.r13[rings.nring3][1] = array[1];
575 rings.r13[rings.nring3][2] = array[2];
576 rings.nring3++;
577 }
578 }
579 }
580 // get four membered rings
581 for(i=1; i <= natom; i++)
582 {
583 for(j=0; j < MAXIAT; j++)
584 {
585 if (atom[i].iat[j] != 0 && atom[i].bo[j] != 9)
586 {
587 jatm = atom[i].iat[j];
588 for (k=0; k < MAXIAT; k++)
589 {
590 if (atom[jatm].iat[k] != 0 && atom[jatm].bo[k] != 9 && atom[jatm].iat[k] != i)
591 {
592 katm = atom[jatm].iat[k];
593 if (!(isbond(katm,i)))
594 {
595 for(l=0; l < MAXIAT; l++)
596 {
597 if (atom[katm].iat[l] != 0 && atom[katm].bo[l] != 9 && atom[katm].iat[l] != jatm)
598 {
599 if (isbond(i,atom[katm].iat[l]) )
600 {
601 array[0] = i;
602 array[1] = jatm;
603 array[2] = katm;
604 array[3] = atom[katm].iat[l];
605 ksort(4, array);
606 if ( is_ring44(array) == FALSE )
607 {
608 rings.r14[rings.nring4][0] = array[0];
609 rings.r14[rings.nring4][1] = array[1];
610 rings.r14[rings.nring4][2] = array[2];
611 rings.r14[rings.nring4][3] = array[3];
612 rings.nring4++;
613 }
614 }
615 }
616 }
617 }
618 }
619 }
620 }
621 }
622 }
623 // get five membered rings
624 for(i=0; i < torsions.ntor; i++)
625 {
626 if (!isbond(torsions.i14[i][0],torsions.i14[i][3]))
627 {
628 for (j=0; j < MAXIAT; j++)
629 {
630 jatm = atom[torsions.i14[i][0]].iat[j];
631 if ( jatm != torsions.i14[i][1] && jatm != torsions.i14[i][2] && jatm != torsions.i14[i][3] )
632 {
633 if (isbond(jatm,torsions.i14[i][3]))
634 {
635 array[0] = jatm;
636 array[1] = torsions.i14[i][0];
637 array[2] = torsions.i14[i][1];
638 array[3] = torsions.i14[i][2];
639 array[4] = torsions.i14[i][3];
640 ksort(5,array);
641 if ( is_ring55(array) == FALSE )
642 {
643 rings.r15[rings.nring5][0] = array[0];
644 rings.r15[rings.nring5][1] = array[1];
645 rings.r15[rings.nring5][2] = array[2];
646 rings.r15[rings.nring5][3] = array[3];
647 rings.r15[rings.nring5][4] = array[4];
648 rings.nring5++;
649 if (rings.nring5 >= natom)
650 message_alert("Error. Too many 5 membered rings!","ERROR");
651
652 }
653 }
654 }
655 }
656 }
657 }
658 // get six membered rings
659 for(i=0; i <= natom; i++)
660 {
661 if (atom[i].type != 5 && atom[i].type != 20)
662 {
663 if (is_cyclo6(i, array))
664 {
665 ksort(6,array);
666 if (is_ring66(array) == FALSE)
667 {
668 rings.r16[rings.nring6][0] = array[0];
669 rings.r16[rings.nring6][1] = array[1];
670 rings.r16[rings.nring6][2] = array[2];
671 rings.r16[rings.nring6][3] = array[3];
672 rings.r16[rings.nring6][4] = array[4];
673 rings.r16[rings.nring6][5] = array[5];
674 rings.nring6++;
675 if (rings.nring6 >= natom)
676 message_alert("Error. Too many 6 membered rings!","ERROR");
677 }
678 }
679 }
680 }
681
682 /* if (rings.nring3 != 0)
683 fprintf(pcmlogfile,"Number of three membered rings: %d\n",rings.nring3);
684 if (rings.nring4 != 0)
685 fprintf(pcmlogfile,"Number of four membered rings: %d\n",rings.nring4);
686 if (rings.nring5 != 0)
687 fprintf(pcmlogfile,"Number of five membered rings: %d\n",rings.nring5);
688 if (rings.nring6 != 0)
689 fprintf(pcmlogfile,"Number of six membered rings: %d\n",rings.nring6); */
690
691 }
692 /* =================================================== */
693 void ksort(int num, int array[])
694 {
695 int i,temp;
696 int found;
697
698 L_1:
699 found = FALSE;
700 for (i=0; i < num-1; i++)
701 {
702 if (array[i+1] < array[i])
703 {
704 temp = array[i];
705 array[i] = array[i+1];
706 array[i+1] = temp;
707 found = TRUE;
708 }
709 }
710 if (found == TRUE)
711 goto L_1;
712 }
713
714
715 /* run the job */
716
717 if (infile) {
718 read_gmmxinp(mmff94filename, infile, flags);
719 }
720
721 /* release temporary strings */
722
723 if(infilename) free(infilename);
724 if(outfilename) free(outfilename);
725 if(mmff94filename) free(mmff94filename);
726 if(mmxfilename) free(mmxfilename);
727 if(logfilename) free(logfilename);
728
729 /* close file handles */
730
731 if(infile && (infile!=stdin)) fclose(infile);
732 if(pcmoutfile && (pcmoutfile!=stdout)) fclose(pcmoutfile);
733 if(pcmlogfile && (pcmlogfile!=stderr)) fclose(pcmlogfile);
734
735 exit(0);
736 }