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