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

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 || use_gast_chrg) charge_dipole();
238 if (pot.use_dipole && !use_gast_chrg) 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) && !use_external_chrg && !use_gast_chrg) 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
404 virial.virx = 0.0;
405 virial.viry = 0.0;
406 virial.virz = 0.0;
407
408 for (i=1; i <= natom; i++)
409 atom[i].energy = 0.0;
410
411 for (i=1; i <= natom; i++)
412 {
413 for (j=0; j < 3; j++)
414 {
415 deriv.deb[i][j] = 0.0;
416 deriv.dea[i][j] = 0.0;
417 deriv.destb[i][j] = 0.0;
418 deriv.deopb[i][j] = 0.0;
419 deriv.detor[i][j] = 0.0;
420 deriv.de14[i][j] = 0.0;
421 deriv.devdw[i][j]= 0.0;
422 deriv.deqq[i][j] = 0.0;
423 deriv.deaa[i][j] = 0.0;
424 deriv.destor[i][j] = 0.0;
425 deriv.deimprop[i][j] = 0.0;
426 deriv.dehb[i][j] = 0.0;
427 deriv.desolv[i][j] = 0.0;
428 }
429 }
430
431 if (pot.use_bond)ebond1();
432 if (pot.use_angle)eangle1();
433 if (pot.use_opbend_wilson) eopbend_wilson1();
434 if (pot.use_tors)etorsion1();
435 if (pot.use_strbnd)estrbnd1();
436
437 // mmff
438 if (pot.use_hal) ehal1();
439 if (pot.use_bufcharge) ebufcharge1();
440
441 energies.total = energies.estr + energies.ebend + energies.etor + energies.estrbnd + energies.e14+
442 energies.evdw + energies.eu + energies.ehbond + energies.eangang + energies.estrtor +
443 energies.eimprop + energies.eimptors + energies.eopb + energies.eurey + energies.esolv;
444 for (i=1; i <= natom; i++)
445 {
446 for (j=0; j < 3; j++)
447 {
448 deriv.d1[i][j] = deriv.deb[i][j] + deriv.dea[i][j] + deriv.deaa[i][j] +
449 deriv.destb[i][j] + deriv.detor[i][j] + deriv.deopb[i][j] + deriv.dehb[i][j] +
450 deriv.destor[i][j] + deriv.deqq[i][j] + deriv.devdw[i][j] + deriv.de14[i][j] +
451 deriv.deimprop[i][j] + deriv.deub[i][j] + deriv.desolv[i][j];
452
453 }
454 }
455 }
456 // ==================================================================
457 double energy()
458 {
459 double etot;
460 int i;
461 // struct timeb start,end;
462
463 energies.total = 0.0;
464 energies.estr = 0.0;
465 energies.ebend = 0.0;
466 energies.etor = 0.0;
467 energies.estrbnd = 0.0;
468 energies.evdw = 0.0;
469 energies.e14 = 0.0;
470 energies.ehbond = 0.0;
471 energies.eu = 0.0;
472 energies.eangang = 0.0;
473 energies.estrtor = 0.0;
474 energies.eimprop = 0.0;
475 energies.eimptors = 0.0;
476 energies.eopb = 0.0;
477 energies.eurey = 0.0;
478 energies.esolv = 0.0;
479
480 for (i=1; i <= natom; i++)
481 atom[i].energy = 0.0;
482
483 if (pot.use_bond)ebond();
484 if (pot.use_angle)eangle();
485 if (pot.use_opbend_wilson) eopbend_wilson();
486 if (pot.use_tors)etorsion();
487 if (pot.use_strbnd)estrbnd();
488 // mmff
489 if (pot.use_hal) ehal();
490 if (pot.use_bufcharge) ebufcharge();
491
492 energies.total = energies.estr + energies.ebend + energies.etor + energies.estrbnd
493 + energies.evdw + energies.e14 + energies.ehbond + energies.eu + energies.eangang +
494 energies.estrtor + energies.eimprop + energies.eimptors + energies.eopb + energies.eurey + energies.esolv;
495 etot = energies.total;
496 return (etot);
497 }
498 /* =========================================================================== */
499 double print_energy()
500 {
501 double etot;
502 int i;
503
504 energies.total = 0.0;
505 energies.estr = 0.0;
506 energies.ebend = 0.0;
507 energies.etor = 0.0;
508 energies.estrbnd = 0.0;
509 energies.evdw = 0.0;
510 energies.e14 = 0.0;
511 energies.ehbond = 0.0;
512 energies.eu = 0.0;
513 energies.eangang = 0.0;
514 energies.estrtor = 0.0;
515 energies.eimprop = 0.0;
516 energies.eimptors = 0.0;
517 energies.eopb = 0.0;
518 energies.eurey = 0.0;
519 energies.esolv = 0.0;
520
521 for (i=1; i <= natom; i++)
522 atom[i].energy = 0.0;
523
524 if (pot.use_bond)ebond();
525 if (pot.use_angle)eangle();
526 if (pot.use_opbend_wilson) eopbend_wilson();
527 if (pot.use_tors)etorsion();
528 if (pot.use_strbnd)estrbnd();
529 // mmff
530 if (pot.use_hal) ehal();
531 if (pot.use_bufcharge) ebufcharge();
532
533 energies.total = energies.estr + energies.ebend + energies.etor + energies.estrbnd
534 + energies.evdw + energies.e14 + energies.ehbond + energies.eu + energies.eangang +
535 energies.estrtor + energies.eimprop + energies.eimptors + energies.eopb + energies.eurey + energies.esolv;
536 etot = energies.total;
537
538 return (etot);
539 }
540
541 /* ================================================================== */
542
543 int iscoord_bond(int i, int j)
544 {
545 int k;
546 long int mask;
547
548 mask = 1 << METCOORD_MASK ;
549
550 for (k=0; k < MAXIAT; k++)
551 {
552 if (atom[i].iat[k] == j && atom[i].bo[k] == 9)
553 {
554 if (atom[i].type >= 300)
555 {
556 if (atom[j].type == 2 || atom[j].type == 4 || atom[j].type == 29 ||
557 atom[j].type == 30 || atom[j].type == 40 || atom[j].type == 48 )
558 return TRUE;
559 }
560 if (atom[j].type >= 300)
561 {
562 if (atom[i].type == 2 || atom[i].type == 4 || atom[i].type == 29 ||
563 atom[i].type == 30 || atom[i].type == 40 || atom[i].type == 48 )
564 return TRUE;
565 }
566 }
567 }
568 // throughout Metal-lone pairs if coordinated bond
569 // metal to donor atom
570 if ( (atom[i].type >= 300) && atom[j].flags & mask)
571 return TRUE;
572 if ( (atom[j].type >= 300) && atom[i].flags & mask)
573 return TRUE;
574
575 return FALSE;
576 }
577
578 void get_coordbonds()
579 {
580 int i, j, k, iatype, jatm;
581
582 cbonds.nbnd = 0;
583 pot.use_coordb = FALSE;
584
585 for (i=1; i <= natom; i++)
586 {
587 if (atom[i].mmx_type >= 300)
588 {
589 for (j=0; j <= MAXIAT; j++)
590 {
591 if (atom[i].bo[j] == 9) // coord bond
592 {
593 iatype = atom[atom[i].iat[j]].mmx_type;
594 if (iatype == 2 || iatype == 4 || iatype == 29 ||
595 iatype == 30 || iatype == 40 || iatype == 48 )
596 {
597 cbonds.icbonds[cbonds.nbnd][0] = i;
598 cbonds.icbonds[cbonds.nbnd][1] = atom[i].iat[j];
599 cbonds.icbonds[cbonds.nbnd][2] = 0;
600 cbonds.nbnd++;
601 } else
602 {
603 jatm = atom[i].iat[j];
604 for (k = 0; k < MAXIAT; k++)
605 {
606 if (atom[atom[jatm].iat[k]].mmx_type == 20) // lp
607 {
608 cbonds.icbonds[cbonds.nbnd][0] = i;
609 cbonds.icbonds[cbonds.nbnd][1] = jatm;
610 cbonds.icbonds[cbonds.nbnd][2] = atom[jatm].iat[k];
611 cbonds.nbnd++;
612 }
613 }
614 }
615 if (cbonds.nbnd > MAXCBND)
616 {
617 fprintf(pcmoutfile,"Error - Too many coordinated bonds!\nProgram will now quit.\n");
618 exit(0);
619 }
620 }
621 }
622 }
623 }
624 if (cbonds.nbnd > 0)
625 pot.use_coordb = TRUE;
626 }
627 int isbond(int i, int j)
628 {
629 int k;
630 for (k=0; k < MAXIAT; k++)
631 {
632 if (atom[i].iat[k] == j && atom[i].bo[k] != 9)
633 return TRUE;
634 }
635 return FALSE;
636 }
637
638 void get_bonds()
639 {
640 int i, j;
641 long int mask;
642
643 bonds_ff.nbnd = 0;
644 mask = 1L << 0; // flag for pi atoms
645
646 for (i=1; i <= natom; i++)
647 {
648 for (j=i+1; j <= natom; j++)
649 {
650 if (isbond(i,j))
651 {
652 bonds_ff.i12[bonds_ff.nbnd][0] = i;
653 bonds_ff.i12[bonds_ff.nbnd][1] = j;
654 bonds_ff.nbnd++;
655 if (bonds_ff.nbnd > MAXBND)
656 {
657 fprintf(pcmoutfile,"Error - Too many bonds!\nProgram will now quit.\n");
658 exit(0);
659 }
660 }
661 }
662 }
663 }