ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/freemol/trunk/src/mengine/src/pcm7.c
Revision: 63
Committed: Wed Dec 3 03:46:32 2008 UTC (13 years ago) by gilbertke
File size: 17549 byte(s)
Log Message:
updated read_sdf type_mmx and first pass at using best practices
Line File contents
1 #define EXTERN extern
2
3 #include "pcwin.h"
4 #include "pcmod.h"
5
6 #include "angles.h"
7
8 #include "energies.h"
9 #include "rings.h"
10 #include "torsions.h"
11 #include "bonds_ff.h"
12 #include "derivs.h"
13 #include "hess.h"
14 #include "pot.h"
15 #include "gmmx.h"
16 #include "utility.h"
17 #include "solv.h"
18 #include "field.h"
19 #include "dipmom.h"
20
21 struct t_cbonds {
22 int nbnd, icbonds[MAXCBND][3], lpd[MAXCBND], ia1[MAXCBND],ia2[MAXCBND];
23 double vrad[MAXCBND],veps[MAXCBND];
24 } cbonds;
25 EXTERN struct t_high_coord {
26 int ncoord, i13[400][3];
27 } high_coord;
28
29 EXTERN struct t_minim_control {
30 int type, method, field, added_const;
31 char added_path[256],added_name[256];
32 } minim_control;
33
34 EXTERN struct t_minim_values {
35 int iprint, ndc, nconst;
36 float dielc;
37 } minim_values;
38 struct t_solvent {
39 int type;
40 double doffset, p1,p2,p3,p4,p5;
41 double *shct,*asolv,*rsolv,*vsolv,*gpol,*rborn;
42 } solvent;
43
44 int Missing_constants;
45 FILE *errfile;
46
47 FILE * fopen_path ( char * , char * , char * ) ;
48 void message_alert(char *,char *);
49 void get_added_const(void);
50 void assign_gast_charge();
51 void lattice(void);
52 void bounds(void);
53 void pcm7_min(void);
54 int iscoord_bond(int, int);
55 void get_coordbonds(void);
56 int isbond(int, int);
57 int ishbond(int, int);
58 int istorsion(int, int);
59 void get_hbond(void);
60 void mark_hbond(void);
61 void reset_atom_data(void);
62 int initialise(void);
63 void get_bonds(void);
64 void get_torsions(void);
65 void get_rings(void);
66 void set_field(void);
67 double energy(void);
68 double print_energy(void);
69 void pcmfout(int);
70 void orbital(void);
71 void free_orbit(void);
72 void set_active(void);
73 void estrtor(void);
74 void echarge(void);
75 void ewald(void);
76 void eurey(void);
77 void ebuck(void);
78 void ebuckmm3(void);
79 void estrbnd(void);
80 void edipole(void);
81 void etorsion(void);
82 void eangang(void);
83 void eangle(void);
84 void ebond(void);
85 void eopbend(void);
86 void ehbond(void);
87 void eimprop(void);
88 void eimptors(void);
89 void elj(void);
90 void ecoord_bond(void);
91 void elj_qq(void);
92 void ebuck_charge(void);
93 void ehigh_coord(void);
94
95 void estrtor1(void);
96 void echarge1(void);
97 void ebuck1(void);
98 void ebuckmm31(void);
99 void estrbnd1(void);
100 void edipole1(void);
101 void etorsion1(void);
102 void eangang1(void);
103 void eangle1(void);
104 void ebond1(void);
105 void eopbend1(void);
106 void ehbond1(void);
107 void eimprop1(void);
108 void eimptors1(void);
109 void elj1(void);
110 void ecoord_bond1(void);
111 void ewald1(void);
112 void eurey1(void);
113 void elj_qq1(void);
114 void ebuck_charge1(void);
115 void ehigh_coord1(void);
116
117 void kangle(void);
118 void korbit(void);
119 void ktorsion(void);
120 void kdipole(void);
121 void kstrbnd(void);
122 void kcharge(void);
123 void ksolv(void);
124 void kangang(void);
125 void kopbend(void);
126 void kstrtor(void);
127 void khbond(void);
128 void piseq(int, int);
129 void kimprop(void);
130 void kimptors(void);
131 void kcoord_bonds(void);
132 void kvdw(void);
133 int kbond(void);
134 void kewald(void);
135 void kurey(void);
136
137 void get_memory(void);
138 void free_memory(void);
139 void gradient(void);
140 void minimize(void);
141 void dynamics(void);
142 void read_datafiles(char *);
143 void attach(void);
144 void eheat(void);
145 void pirite(void);
146 void piden(void);
147 void hdel(int);
148 void set_atomtypes(int);
149 void type(void);
150 void generate_bonds(void);
151 void zero_data(void);
152 void ehal(void);
153 void ehal1(void);
154 void ebufcharge(void);
155 void ebufcharge1(void);
156 void eopbend_wilson(void);
157 void eopbend_wilson1(void);
158 void hbondreset(void);
159 void vibrate(void);
160 void pireset(void);
161 void adjust_mmfftypes(void);
162 void build_neighbor_list(int);
163 int setup_calculation(void);
164 void end_calculation(void);
165 void egeom1(void);
166 void esolv1(void);
167 void egeom(void);
168 void esolv(void);
169 void xlogp(float *);
170
171 void mmxsub(int ia)
172 {
173 int nret;
174
175 minim_control.type = ia; // 0 to minimize, 1 for single point
176 nret = setup_calculation();
177 if (nret == FALSE)
178 {
179 end_calculation();
180 return;
181 }
182
183 pcm7_min();
184 end_calculation();
185 }
186 // ================================================================
187 void pcm7_min()
188 {
189 int i, print;
190 double etot;
191
192 print = FALSE;
193 if (minim_values.iprint == TRUE)
194 {
195 print = TRUE;
196 minim_values.iprint = FALSE;
197 }
198
199 if (minim_control.type == 1) // single point calculation
200 return;
201
202 minimize();
203
204 if (print == TRUE)
205 minim_values.iprint = TRUE;
206
207 if (minim_values.iprint == TRUE)
208 etot = print_energy();
209 else
210 etot = energy();
211
212 // compute dipole moment here
213 dipolemom.xdipole = 0.0;
214 dipolemom.ydipole = 0.0;
215 dipolemom.zdipole = 0.0;
216 if (pot.use_charge || pot.use_bufcharge) charge_dipole();
217 if (pot.use_dipole) dipole_dipole();
218 dipolemom.total = sqrt(dipolemom.xdipole*dipolemom.xdipole +
219 dipolemom.ydipole*dipolemom.ydipole +
220 dipolemom.zdipole*dipolemom.zdipole);
221 }
222 // ==================================================================
223 int setup_calculation()
224 {
225 int nRet;
226 char string[30];
227 double etot;
228 if (minim_control.field == MM3)
229 {
230 // copy mm3 types into type file and retype
231 pot.use_hbond = FALSE;
232 hbondreset();
233 hdel(1);
234 strcpy(string,"mm3.prm");
235 zero_data();
236 read_datafiles(string);
237 set_field();
238 default_intype = MM3;
239 default_outtype = MM3;
240 set_atomtypes(MM3);
241 generate_bonds();
242 if (minim_control.added_const)
243 get_added_const();
244 } else if (minim_control.field == MMFF94)
245 {
246 // copy mm3 types into type file and retype
247 pot.use_hbond = FALSE;
248 pot.use_picalc = FALSE;
249 hbondreset();
250 pireset();
251 hdel(1);
252 set_field();
253 default_intype = MMFF94;
254 default_outtype = MMFF94;
255 set_atomtypes(MMFF94);
256 generate_bonds();
257 if (minim_control.added_const)
258 get_added_const();
259 } else if (minim_control.field == AMBER)
260 {
261 pot.use_hbond = FALSE;
262 pot.use_picalc = FALSE;
263 hbondreset();
264 pireset();
265 hdel(2);
266 set_field();
267 generate_bonds();
268 if (minim_control.added_const)
269 get_added_const();
270 } else if (minim_control.field == OPLSAA)
271 {
272 pot.use_hbond = FALSE;
273 pot.use_picalc = FALSE;
274 hbondreset();
275 pireset();
276 hdel(1);
277 set_field();
278 generate_bonds();
279 if (minim_control.added_const)
280 get_added_const();
281 } else
282 {
283 if (minim_control.added_const)
284 {
285 strcpy(string,"mmxconst.prm");
286 zero_data();
287 read_datafiles(string);
288 get_added_const();
289 }
290 type();
291 set_field();
292 generate_bonds();
293 }
294
295
296 // allocate memeory for derivatives
297 get_memory();
298
299 // setup bonds, angles, torsions, improper torsions and angles, allenes
300 get_bonds();
301 get_angles();
302 get_torsions();
303 get_coordbonds();
304
305 if (high_coord.ncoord > 0)
306 pot.use_highcoord = TRUE;
307
308 attach();
309
310
311 get_rings();
312 // need allene
313
314 // set active atoms
315 set_active();
316
317 /* assign local geometry potential function parameters */
318 Missing_constants = FALSE;
319 // errfile = fopen_path(pcwindir,"pcmod.err","w");
320 if (pot.use_bond || pot.use_strbnd) nRet = kbond();
321 if (nRet == FALSE) // use end_calc to free memory
322 {
323 // message_alert("Parameters missing. See pcmod.err for information","Error");
324 // fclose(errfile);
325 energies.total = 9999.;
326 return FALSE;
327 }
328 if (pot.use_angle || pot.use_strbnd || pot.use_angang) kangle();
329 if (pot.use_angle || pot.use_opbend || pot.use_opbend_wilson) kopbend();
330 if (pot.use_tors || pot.use_strtor || pot.use_tortor) ktorsion();
331 if (pot.use_strbnd) kstrbnd();
332
333 if (pot.use_buck || pot.use_hal || pot.use_lj) kvdw();
334 if ((pot.use_charge || pot.use_bufcharge || pot.use_ewald)) kcharge();
335
336 if (Missing_constants == TRUE)
337 {
338 // message_alert("Parameters missing. See pcmod.err for information","Error");
339 // fclose(errfile);
340 energies.total = 9999.0;
341 return (FALSE);
342 } else
343 {
344 // fclose(errfile);
345 // remove("pcmod.err");
346 }
347 if (minim_values.iprint == TRUE)
348 etot = print_energy();
349 else
350 etot = energy();
351 return TRUE;
352 }
353 // =================================================================
354 void end_calculation()
355 {
356
357 free_memory();
358 }
359 // ==================================================================
360 void gradient()
361 {
362 int i, j;
363
364 energies.estr = 0.0;
365 energies.ebend = 0.0;
366 energies.estrbnd = 0.0;
367 energies.e14 = 0.0;
368 energies.evdw = 0.0;
369 energies.etor = 0.0;
370 energies.eu = 0.0;
371 energies.eopb = 0.0;
372 energies.eangang = 0.0;
373 energies.estrtor = 0.0;
374 energies.ehbond = 0.0;
375 energies.efix = 0.0;
376 energies.eimprop = 0.0;
377 energies.eimptors = 0.0;
378 energies.eurey = 0.0;
379 energies.esolv = 0.0;
380 energies.egeom =0.0;
381
382 virial.virx = 0.0;
383 virial.viry = 0.0;
384 virial.virz = 0.0;
385
386 for (i=1; i <= natom; i++)
387 atom[i].energy = 0.0;
388
389 for (i=1; i <= natom; i++)
390 {
391 for (j=0; j < 3; j++)
392 {
393 deriv.deb[i][j] = 0.0;
394 deriv.dea[i][j] = 0.0;
395 deriv.destb[i][j] = 0.0;
396 deriv.deopb[i][j] = 0.0;
397 deriv.detor[i][j] = 0.0;
398 deriv.de14[i][j] = 0.0;
399 deriv.devdw[i][j]= 0.0;
400 deriv.deqq[i][j] = 0.0;
401 deriv.deaa[i][j] = 0.0;
402 deriv.destor[i][j] = 0.0;
403 deriv.deimprop[i][j] = 0.0;
404 deriv.dehb[i][j] = 0.0;
405 deriv.desolv[i][j] = 0.0;
406 deriv.degeom[i][j] = 0.0;
407 }
408 }
409
410 if (pot.use_bond)ebond1();
411 if (pot.use_angle)eangle1();
412 if (pot.use_opbend_wilson) eopbend_wilson1();
413 if (pot.use_tors)etorsion1();
414 if (pot.use_strbnd)estrbnd1();
415
416 // mmff
417 if (pot.use_hal) ehal1();
418 if (pot.use_bufcharge) ebufcharge1();
419
420 if (pot.use_geom) egeom1();
421
422 energies.total = energies.estr + energies.ebend + energies.etor + energies.estrbnd + energies.e14+
423 energies.evdw + energies.eu + energies.ehbond + energies.eangang + energies.estrtor +
424 energies.eimprop + energies.eimptors + energies.eopb + energies.eurey + energies.esolv + energies.egeom;
425 for (i=1; i <= natom; i++)
426 {
427 for (j=0; j < 3; j++)
428 {
429 deriv.d1[i][j] = deriv.deb[i][j] + deriv.dea[i][j] + deriv.deaa[i][j] +
430 deriv.destb[i][j] + deriv.detor[i][j] + deriv.deopb[i][j] + deriv.dehb[i][j] +
431 deriv.destor[i][j] + deriv.deqq[i][j] + deriv.devdw[i][j] + deriv.de14[i][j] +
432 deriv.deimprop[i][j] + deriv.deub[i][j] + deriv.desolv[i][j] + deriv.degeom[i][j];
433
434 }
435 }
436 }
437 // ==================================================================
438 double energy()
439 {
440 double etot;
441 int i;
442 // struct timeb start,end;
443
444 energies.total = 0.0;
445 energies.estr = 0.0;
446 energies.ebend = 0.0;
447 energies.etor = 0.0;
448 energies.estrbnd = 0.0;
449 energies.evdw = 0.0;
450 energies.e14 = 0.0;
451 energies.ehbond = 0.0;
452 energies.eu = 0.0;
453 energies.eangang = 0.0;
454 energies.estrtor = 0.0;
455 energies.eimprop = 0.0;
456 energies.eimptors = 0.0;
457 energies.eopb = 0.0;
458 energies.eurey = 0.0;
459 energies.esolv = 0.0;
460 energies.egeom = 0.0;
461
462 for (i=1; i <= natom; i++)
463 atom[i].energy = 0.0;
464
465 if (pot.use_bond)ebond();
466 if (pot.use_angle)eangle();
467 if (pot.use_opbend_wilson) eopbend_wilson();
468 if (pot.use_tors)etorsion();
469 if (pot.use_strbnd)estrbnd();
470 // mmff
471 if (pot.use_hal) ehal();
472 if (pot.use_bufcharge) ebufcharge();
473
474 if (pot.use_geom) egeom();
475
476 energies.total = energies.estr + energies.ebend + energies.etor + energies.estrbnd
477 + energies.evdw + energies.e14 + energies.ehbond + energies.eu + energies.eangang +
478 energies.estrtor + energies.eimprop + energies.eimptors + energies.eopb + energies.eurey + energies.esolv + energies.egeom;
479 etot = energies.total;
480 return (etot);
481 }
482 /* =========================================================================== */
483 double print_energy()
484 {
485 double etot;
486 int i;
487
488 energies.total = 0.0;
489 energies.estr = 0.0;
490 energies.ebend = 0.0;
491 energies.etor = 0.0;
492 energies.estrbnd = 0.0;
493 energies.evdw = 0.0;
494 energies.e14 = 0.0;
495 energies.ehbond = 0.0;
496 energies.eu = 0.0;
497 energies.eangang = 0.0;
498 energies.estrtor = 0.0;
499 energies.eimprop = 0.0;
500 energies.eimptors = 0.0;
501 energies.eopb = 0.0;
502 energies.eurey = 0.0;
503 energies.esolv = 0.0;
504 energies.egeom = 0.0;
505
506 for (i=1; i <= natom; i++)
507 atom[i].energy = 0.0;
508
509 if (pot.use_bond)ebond();
510 if (pot.use_angle)eangle();
511 if (pot.use_opbend_wilson) eopbend_wilson();
512 if (pot.use_tors)etorsion();
513 if (pot.use_strbnd)estrbnd();
514 // mmff
515 if (pot.use_hal) ehal();
516 if (pot.use_bufcharge) ebufcharge();
517
518 if (pot.use_geom) egeom();
519
520 energies.total = energies.estr + energies.ebend + energies.etor + energies.estrbnd
521 + energies.evdw + energies.e14 + energies.ehbond + energies.eu + energies.eangang +
522 energies.estrtor + energies.eimprop + energies.eimptors + energies.eopb + energies.eurey + energies.esolv + energies.egeom;
523 etot = energies.total;
524
525 return (etot);
526 }
527
528 /* ================================================================== */
529
530 int iscoord_bond(int i, int j)
531 {
532 int k;
533 long int mask;
534
535 mask = 1 << METCOORD_MASK ;
536
537 for (k=0; k < MAXIAT; k++)
538 {
539 if (atom[i].iat[k] == j && atom[i].bo[k] == 9)
540 {
541 if (atom[i].type >= 300)
542 {
543 if (atom[j].type == 2 || atom[j].type == 4 || atom[j].type == 29 ||
544 atom[j].type == 30 || atom[j].type == 40 || atom[j].type == 48 )
545 return TRUE;
546 }
547 if (atom[j].type >= 300)
548 {
549 if (atom[i].type == 2 || atom[i].type == 4 || atom[i].type == 29 ||
550 atom[i].type == 30 || atom[i].type == 40 || atom[i].type == 48 )
551 return TRUE;
552 }
553 }
554 }
555 // throughout Metal-lone pairs if coordinated bond
556 // metal to donor atom
557 if ( (atom[i].type >= 300) && atom[j].flags & mask)
558 return TRUE;
559 if ( (atom[j].type >= 300) && atom[i].flags & mask)
560 return TRUE;
561
562 return FALSE;
563 }
564
565 void get_coordbonds()
566 {
567 int i, j, k, iatype, jatm;
568
569 cbonds.nbnd = 0;
570 pot.use_coordb = FALSE;
571
572 for (i=1; i <= natom; i++)
573 {
574 if (atom[i].mmx_type >= 300)
575 {
576 for (j=0; j <= MAXIAT; j++)
577 {
578 if (atom[i].bo[j] == 9) // coord bond
579 {
580 iatype = atom[atom[i].iat[j]].mmx_type;
581 if (iatype == 2 || iatype == 4 || iatype == 29 ||
582 iatype == 30 || iatype == 40 || iatype == 48 )
583 {
584 cbonds.icbonds[cbonds.nbnd][0] = i;
585 cbonds.icbonds[cbonds.nbnd][1] = atom[i].iat[j];
586 cbonds.icbonds[cbonds.nbnd][2] = 0;
587 cbonds.nbnd++;
588 } else
589 {
590 jatm = atom[i].iat[j];
591 for (k = 0; k < MAXIAT; k++)
592 {
593 if (atom[atom[jatm].iat[k]].mmx_type == 20) // lp
594 {
595 cbonds.icbonds[cbonds.nbnd][0] = i;
596 cbonds.icbonds[cbonds.nbnd][1] = jatm;
597 cbonds.icbonds[cbonds.nbnd][2] = atom[jatm].iat[k];
598 cbonds.nbnd++;
599 }
600 }
601 }
602 if (cbonds.nbnd > MAXCBND)
603 {
604 fprintf(pcmlogfile,"Error - Too many coordinated bonds!\nProgram will now quit.\n");
605 exit(0);
606 }
607 }
608 }
609 }
610 }
611 if (cbonds.nbnd > 0)
612 pot.use_coordb = TRUE;
613 }
614 int isbond(int i, int j)
615 {
616 int k;
617 for (k=0; k < MAXIAT; k++)
618 {
619 if (atom[i].iat[k] == j && atom[i].bo[k] != 9)
620 return TRUE;
621 }
622 return FALSE;
623 }
624
625 void get_bonds()
626 {
627 int i, j;
628 long int mask;
629
630 bonds_ff.nbnd = 0;
631 mask = 1L << 0; // flag for pi atoms
632
633 for (i=1; i <= natom; i++)
634 {
635 for (j=i+1; j <= natom; j++)
636 {
637 if (isbond(i,j))
638 {
639 bonds_ff.i12[bonds_ff.nbnd][0] = i;
640 bonds_ff.i12[bonds_ff.nbnd][1] = j;
641 bonds_ff.nbnd++;
642 if (bonds_ff.nbnd > MAXBND)
643 {
644 fprintf(pcmlogfile,"Error - Too many bonds!\nProgram will now quit.\n");
645 exit(0);
646 }
647 }
648 }
649 }
650 }