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, 1 month ago) by wdelano
File size: 18317 byte(s)
Log Message:
code merge 20081130
Line User Rev File contents
1 tjod 3 #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 wdelano 58 if (pot.use_charge || pot.use_bufcharge) charge_dipole();
238     if (pot.use_dipole) dipole_dipole();
239 tjod 3 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 wdelano 58 if ((pot.use_charge || pot.use_bufcharge || pot.use_ewald)) kcharge();
357 tjod 3
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 wdelano 58 energies.egeom =0.0;
404 tjod 3
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 wdelano 58 deriv.degeom[i][j] = 0.0;
430 tjod 3 }
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 wdelano 58
443     if (pot.use_geom) egeom1();
444 tjod 3
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 wdelano 58 energies.eimprop + energies.eimptors + energies.eopb + energies.eurey + energies.esolv + energies.egeom;
448 tjod 3 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 wdelano 58 deriv.deimprop[i][j] + deriv.deub[i][j] + deriv.desolv[i][j] + deriv.degeom[i][j];
456 tjod 3
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 wdelano 58 energies.egeom = 0.0;
484 tjod 3
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 wdelano 58
497     if (pot.use_geom) egeom();
498 tjod 3
499     energies.total = energies.estr + energies.ebend + energies.etor + energies.estrbnd
500     + energies.evdw + energies.e14 + energies.ehbond + energies.eu + energies.eangang +
501 wdelano 58 energies.estrtor + energies.eimprop + energies.eimptors + energies.eopb + energies.eurey + energies.esolv + energies.egeom;
502 tjod 3 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 wdelano 58 energies.egeom = 0.0;
528 tjod 3
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 wdelano 58 if (pot.use_geom) egeom();
542    
543 tjod 3 energies.total = energies.estr + energies.ebend + energies.etor + energies.estrbnd
544     + energies.evdw + energies.e14 + energies.ehbond + energies.eu + energies.eangang +
545 wdelano 58 energies.estrtor + energies.eimprop + energies.eimptors + energies.eopb + energies.eurey + energies.esolv + energies.egeom;
546 tjod 3 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 wdelano 58 fprintf(pcmlogfile,"Error - Too many coordinated bonds!\nProgram will now quit.\n");
628 tjod 3 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 wdelano 58 fprintf(pcmlogfile,"Error - Too many bonds!\nProgram will now quit.\n");
668 tjod 3 exit(0);
669     }
670     }
671     }
672     }
673     }