ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/freemol/trunk/src/mengine/src/pcm7.c
Revision: 114
Committed: Wed Apr 29 19:07:31 2009 UTC (11 years, 3 months ago) by gilbertke
File size: 27635 byte(s)
Log Message:
added code for amber energy functions
Line User Rev File contents
1 gilbertke 103 #define EXTERN extern
2    
3     #include "pcwin.h"
4     #include "pcmod.h"
5    
6     #include "angles.h"
7     #include "attached.h"
8     #include "energies.h"
9     #include "bonds_ff.h"
10     #include "derivs.h"
11     #include "hess.h"
12     #include "utility.h"
13     #include "cutoffs.h"
14     #include "nonbond.h"
15     #include "torsions.h"
16 gilbertke 110 #include "fix.h"
17 gilbertke 103
18     EXTERN struct t_strbnd {
19 gilbertke 105 int nstrbnd, **isb;
20     float *ksb1, *ksb2;
21 gilbertke 103 } strbnd;
22 gilbertke 114 EXTERN struct t_improp {
23     int nimptors, **iiprop;
24     float *v1, *v2, *v3;
25     int *ph1, *ph2, *ph3;
26     } improp;
27 gilbertke 103
28     EXTERN struct t_minim_control {
29     int type, method, field, added_const;
30     char added_path[256],added_name[256];
31     } minim_control;
32    
33     EXTERN struct t_minim_values {
34     int iprint, ndc, nconst;
35     float dielc;
36     } minim_values;
37     struct t_solvent {
38     int type;
39     double EPSin, EPSsolv;
40     double doffset, p1,p2,p3,p4,p5;
41     double *shct,*asolv,*rsolv,*vsolv,*gpol,*rborn;
42     } solvent;
43    
44     int Missing_constants;
45     EXTERN double hesscut;
46    
47     void message_alert(char *,char *);
48     void pcm7_min(void);
49     int initialise(void);
50     void attach(void);
51     void get_bonds(void);
52     void get_angles(void);
53 gilbertke 110 void get_torsions(int natom,int *type,int **iat,int **bo);
54 gilbertke 103 void set_field(int);
55     double energy(void);
56     double get_total_energy(void);
57     double get_total_deriv_x(int i);
58     double get_total_deriv_y(int i);
59     double get_total_deriv_z(int i);
60     double print_energy(void);
61 gilbertke 110 void set_active(int natom,int *use,int natom_fix,int *katm_fix);
62 gilbertke 103 int setup_calculation(void);
63     void end_calculation(void);
64    
65 gilbertke 114 void ebond(int nbnd,int **,int *use,double *x,double *y,double *z,double *bl,double *bk,double *estr);
66     void ebond1(int natom,int nbnd,int **,int *use,double *x,double *y,double *z,double *bl,double *bk,double *estr,double **deb);
67 gilbertke 105 void eangle(int nang,int **,int *use,double *x,double *y,double *z,int *angin,int *angtype,float *anat,float *acon,double *ebend);
68     void eangle1(int nang,int natom,int **,int *use,double *x,double *y,double *z,int *angin,int *angtype,float *anat,float *acon,double *ebend,double **dea);
69 gilbertke 103 void ebufcharge(int natom,int *use,int **skip,double *x,double *y,double *z,double *charge,double chrgcut,double *eu);
70     void ebufcharge1(int natom,int *use,int **skip,double *x,double *y,double *z,double *charge,double chrgcut,double *eu,double **deqq);
71     void ehal(int natom,int *type,int *use,double *x,double *y,double *z,double vdwcut,int **skip,double **vrad,double **veps,double *evdw,double *e14);
72     void ehal1(int natom,int *type,int *use,double *x,double *y,double *z,double vdwcut,int **skip,double **vrad,double **veps,double *evdw,double *e14,
73     double **devdw,double **de14);
74     void etorsion(int ntor,int **i14,int *use,int *type,double *x,double *y,double *z,int *phin1,int *phin2,int *phin3,int *phin4,int *phin5,int *phin6,float *vin1,float *vin2,
75     float *vin3,float *vin4,float *vin5,float *vin6,double *etor);
76     void etorsion1(int natom,int ntor,int **i14,int *use,double *x,double *y,double *z,int *phin1,int *phin2,int *phin3,int *phin4,int *phin5,int *phin6,float *vin1,float *vin2,
77     float *vin3,float *vin4,float *vin5,float *vin6,double *etor,double **detor);
78 gilbertke 105 void eopbend_wilson(int nang,int **,int *use,double *x,double *y,double *z,int *angin,float *copb,double *eopb);
79     void eopbend_wilson1(int natom,int nang,int **,int *use,double *x,double *y,double *z,int *angin,float *copb,double *eopb,double **deopb);
80     void estrbnd(int nstrbnd,int **,int **,int *use,double *x,double *y,double *z,float *anat,double *bl,double *bk,float *ksb1,
81 gilbertke 103 float *ksb2,double *estrbnd);
82 gilbertke 105 void estrbnd1(int natom,int nstrbnd,int **,int **,int *use,double *x,double *y,double *z,float *anat,double *bl,double *bk,float *ksb1,
83 gilbertke 103 float *ksb2,double *estrbnd,double **destb);
84     void egeom(int natom,int *use,double *x,double *y,double *z,double *egeom);
85     void egeom1(int natom,int *use,double *x,double *y,double *z,double *egeom,double **degeom);
86 gilbertke 110 void esolv(int natom,double chrgcut,int **skip,int *use,double *charge,double *x,double *y,double *z,double *esolv);
87     void esolv1(int natom,double chrgcut,int **skip,int *use,double *charge,double *x,double *y,double *z,double *esolv,double **desolv,double *drb);
88 gilbertke 103
89     void hessian(int, double *,int *, int *, int *,double *);
90 gilbertke 105 void ebond2(int ia,int **,int **,int **,double *x,double *y,double *z,double *bl,double *bk,float **hessx,float **hessy,float **hessz);
91     void eangle2(int i,int nang,int **,int *angin, double *x,double *y,double *z,float *acon,float *anat,double **dea,float **hessx,float **hessy,float **hessz);
92 gilbertke 103 void etorsion2(int iatom,int ntor,int **i14,int *use,double *x,double *y,double *z,int *phin1,int *phin2,int *phin3,int *phin4,int *phin5,int *phin6,float *vin1,float *vin2,
93     float *vin3,float *vin4,float *vin5,float *vin6,float **hessx,float **hessy,float **hessz);
94     void ehal2(int iatom,int natom,int *type,int *use,double *x,double *y,double *z,double vdwcut,int **skip,double **vrad,double **veps,
95     float **hessx,float **hessy,float **hessz);
96     void ebufcharge2(int i,int natom,int *use,int **skip,double *x,double *y,double *z,double *charge,double chrgcut,float **hessx,float **hessy,float **hessz);
97 gilbertke 105 void eopbend_wilson2(int i,int nang,int **,int *use,double *x,double *y,double *z,int *angin,float *copb,
98 gilbertke 103 float **hessx,float **hessy,float **hessz);
99 gilbertke 105 void estrbnd2(int iatom,int nstrbnd,int **,int **,int *use,double *x,double *y,double *z,float *anat,double *bl,double *bk,float *ksb1,
100 gilbertke 103 float *ksb2,float **hessx,float **hessy,float **hessz);
101 gilbertke 110 void esolv2(int ia,int natom,double chrgcut,double *charge,double *x,double *y,double *z,float **hessx,float **hessy,float **hessz);
102 gilbertke 103 void egeom2(int ,int natom,int *use,double *x,double *y,double *z,float **hessx,float **hessy,float **hessz);
103 gilbertke 114 // gaff and amber
104     void eimptors(int nimptors,int **iiprop,int *use,int *type,double *x,double *y,double *z,float *v1in,float *v2in,float *v3in,int *ph1in,int *ph2in,int *ph3in, double *e_imp);
105     void eimptors1(int nimptors,int natom,int **iiprop,int *use,int *type,double *x,double *y,double *z,float *v1in,float *v2in,float *v3in,int *ph1in,int *ph2in,int *ph3in,
106     double *e_imp,double **deimprop);
107     void eimptors2(int iatom,int nimptors,int **iiprop,int *use,int *type,double *x,double *y,double *z,float *v1in,float *v2in,float *v3in,int *ph1in,int *ph2in,int *ph3in,
108     float **hessx,float **hessy,float **hessz);
109     void elj(int natom,int *type,int *use,double *x,double *y,double *z,double vdwcut,int **skip,double **vrad,double **veps,double *evdw,double *e14);
110     void elj1(int natom,int *type,int *use,double *x,double *y,double *z,double vdwcut,int **skip,double **vrad,double **veps,double *evdw,double *e14,
111     double **devdw,double **de14);
112     void elj2(int iatom,int natom,int *type,int *use,double *x,double *y,double *z,double vdwcut,int **skip,double **vrad,double **veps,
113     float **hessx,float **hessy,float **hessz);
114     void echarge(void);
115     void echarge1(void);
116     void echarge2(int);
117     //
118 gilbertke 103
119 gilbertke 111 int kbond(void);
120 gilbertke 103 void kangle(void);
121     void kstrbnd(void);
122 gilbertke 111 void kcharge(int natom,int *type,int *atomnum,long int *flags,int **iat,int **bo,double *charge,double *sigma_charge,double *formal_charge);
123     void kopbend(int *type,int *tclass);
124     void ktorsion(int *type,int *tclass,int **iat,int **bo);
125     void kvdw(int natom,int *type,int *atomnum);
126     void ksolv(int natom,int *atomnum,int **iat,int **bo);
127 gilbertke 103
128 gilbertke 114 // gaff
129     void kimptors(void);
130    
131 gilbertke 103 void get_memory(void);
132     void free_memory(void);
133     void gradient(void);
134     void read_datafiles(char *);
135     void get_added_const(void);
136     void hdel(int);
137     void set_atomtypes(int);
138     void type(void);
139     void zero_data(void);
140     void potoff(void);
141     int use_bond(void);
142     int use_angle(void);
143     int use_strbnd(void);
144     int use_opbend_wilson(void);
145     int use_tors(void);
146     int use_strtor(void);
147     int use_hal(void);
148     int use_charge(void);
149     int use_bufcharge(void);
150     int use_geom(void);
151     int use_solv(void);
152 gilbertke 114 int use_imptor(void);
153     int use_lj(void);
154 gilbertke 103
155     // ==================================================================
156     int setup_calculation()
157     {
158     int i,j,nRet;
159     char string[30];
160     double etot;
161    
162     potoff(); // turn off all force field terms before resetting
163     if (minim_control.field == MMFF94)
164     {
165     set_field(MMFF94);
166     set_atomtypes(MMFF94);
167     if (minim_control.added_const)
168     get_added_const();
169 gilbertke 114 } else if (minim_control.field == GAFF)
170     {
171     set_field(GAFF);
172     set_atomtypes(GAFF);
173     } else
174 gilbertke 103 {
175     if (minim_control.added_const)
176     {
177     strcpy(string,"mmxconst.prm");
178     zero_data();
179     read_datafiles(string);
180     get_added_const();
181     }
182     type();
183     set_field(MMX);
184     }
185    
186    
187     // allocate memeory for derivatives
188     get_memory();
189    
190     // setup bonds, angles, torsions, improper torsions and angles, allenes
191     get_bonds();
192     get_angles();
193 gilbertke 110 get_torsions(natom,atom.type,atom.iat,atom.bo);
194 gilbertke 103
195     attach();
196     // need allene
197    
198     // set active atoms
199 gilbertke 110 set_active(natom,atom.use,fx_atom.natom_fix,fx_atom.katom_fix);
200 gilbertke 103 // setup non_bonded list of atoms to skip
201     for (i=1; i <= natom; i++)
202     {
203     for (j=0; j < MAXIAT; j++)
204     if (atom.iat[i][j] != 0 && atom.bo[i][j] != 9)
205     {
206     skip[i][atom.iat[i][j]] = i;
207     skip[atom.iat[i][j]][i] = atom.iat[i][j];
208     }
209     for (j=0; j < attached.n13[i]; j++)
210     {
211     skip[i][attached.i13[j][i]] = i;
212     skip[attached.i13[j][i]][i] = attached.i13[j][i];
213     }
214     for(j=0; j < attached.n14[i]; j++)
215     {
216     skip[i][attached.i14[j][i]] = -i;
217     skip[attached.i14[j][i]][i] = -attached.i14[j][i];
218     }
219     }
220     /* assign local geometry potential function parameters */
221     Missing_constants = FALSE;
222     if (use_bond() || use_strbnd()) nRet = kbond();
223     if (nRet == FALSE) // use end_calc to free memory
224     {
225     energies.total = 9999.;
226     return FALSE;
227     }
228     if (use_angle() || use_strbnd()) kangle();
229 gilbertke 111 if (use_angle() || use_opbend_wilson()) kopbend(atom.type,atom.tclass);
230     if (use_tors() || use_strtor()) ktorsion(atom.type,atom.tclass,atom.iat,atom.bo);
231 gilbertke 103 if (use_strbnd()) kstrbnd();
232 gilbertke 114 if (use_imptor()) kimptors();
233 gilbertke 103
234 gilbertke 111 if (use_hal()) kvdw(natom,atom.type,atom.atomnum);
235     if (use_charge() || use_bufcharge()) kcharge(natom,atom.type,atom.atomnum,atom.flags,atom.iat,atom.bo,
236     atom.charge,atom.sigma_charge,atom.formal_charge);
237     if (use_solv()) ksolv(natom,atom.atomnum,atom.iat,atom.bo);
238 gilbertke 103
239     if (Missing_constants == TRUE)
240     {
241     energies.total = 9999.0;
242     return (FALSE);
243     }
244     if (minim_values.iprint == TRUE)
245     etot = print_energy();
246     else
247     etot = energy();
248     return TRUE;
249     }
250     // =================================================================
251     void end_calculation()
252     {
253    
254     free_memory();
255     }
256     // ==================================================================
257     void gradient()
258     {
259     int i, j;
260    
261     energies.estr = 0.0;
262     energies.ebend = 0.0;
263     energies.estrbnd = 0.0;
264     energies.e14 = 0.0;
265     energies.evdw = 0.0;
266     energies.etor = 0.0;
267     energies.eu = 0.0;
268     energies.eopb = 0.0;
269     energies.eangang = 0.0;
270     energies.estrtor = 0.0;
271     energies.ehbond = 0.0;
272     energies.efix = 0.0;
273     energies.eimprop = 0.0;
274     energies.eimptors = 0.0;
275     energies.eurey = 0.0;
276     energies.esolv = 0.0;
277     energies.egeom =0.0;
278    
279     for (i=1; i <= natom; i++)
280     {
281     for (j=0; j < 3; j++)
282     {
283     deriv.deb[i][j] = 0.0;
284     deriv.dea[i][j] = 0.0;
285     deriv.destb[i][j] = 0.0;
286     deriv.deopb[i][j] = 0.0;
287     deriv.detor[i][j] = 0.0;
288     deriv.de14[i][j] = 0.0;
289     deriv.devdw[i][j]= 0.0;
290     deriv.deqq[i][j] = 0.0;
291     deriv.deaa[i][j] = 0.0;
292     deriv.destor[i][j] = 0.0;
293     deriv.deimprop[i][j] = 0.0;
294     deriv.dehb[i][j] = 0.0;
295     deriv.desolv[i][j] = 0.0;
296     deriv.degeom[i][j] = 0.0;
297     }
298     }
299    
300 gilbertke 114 if (use_bond())ebond1(natom,bonds_ff.nbnd,bonds_ff.i12,atom.use,atom.x,atom.y,atom.z,bonds_ff.bl,bonds_ff.bk,&energies.estr,deriv.deb);
301 gilbertke 103 if (use_angle())eangle1(angles.nang,natom, angles.i13,atom.use,atom.x,atom.y,atom.z,angles.angin,angles.angtype,angles.anat,angles.acon,&energies.ebend,deriv.dea);
302     if (use_tors())etorsion1(natom,torsions.ntor,torsions.i14,atom.use,atom.x,atom.y,atom.z,torsions.ph1,torsions.ph2,torsions.ph3,torsions.ph4,torsions.ph5,torsions.ph6,
303     torsions.v1,torsions.v2,torsions.v3,torsions.v4,torsions.v5,torsions.v6,&energies.etor,deriv.detor);
304 gilbertke 114 // mmff
305 gilbertke 103 if (use_strbnd())estrbnd1(natom,strbnd.nstrbnd,strbnd.isb,angles.i13,atom.use,atom.x,atom.y,atom.z,angles.anat,bonds_ff.bl,bonds_ff.bk,strbnd.ksb1,strbnd.ksb2,
306     &energies.estrbnd,deriv.destb);
307 gilbertke 114 if (use_opbend_wilson()) eopbend_wilson1(natom,angles.nang,angles.i13,atom.use,atom.x,atom.y,atom.z,angles.angin,angles.copb,&energies.eopb,deriv.deopb);
308 gilbertke 103 if (use_hal()) ehal1(natom,atom.type,atom.use,atom.x,atom.y,atom.z,cutoffs.vdwcut,skip,nonbond.vrad,nonbond.veps,&energies.evdw,&energies.e14,deriv.devdw,deriv.de14);
309     if (use_bufcharge()) ebufcharge1(natom,atom.use,skip,atom.x,atom.y,atom.z,atom.charge,cutoffs.chrgcut,&energies.eu,deriv.deqq);
310 gilbertke 114 // gaff
311     if (use_imptor()) eimptors1(improp.nimptors,natom,improp.iiprop,atom.use,atom.type,atom.x,atom.y,atom.z,improp.v1,improp.v2,improp.v3,improp.ph1,improp.ph2,improp.ph3,
312     &energies.eimptors,deriv.deimprop);
313     if (use_lj()) elj1(natom,atom.type,atom.use,atom.x,atom.y,atom.z,cutoffs.vdwcut,skip,nonbond.vrad,nonbond.veps,&energies.evdw,&energies.e14,deriv.devdw,deriv.de14);
314     // if (use_charge())
315 gilbertke 103
316     if (use_geom()) egeom1(natom,atom.use,atom.x,atom.y,atom.z,&energies.egeom,deriv.degeom);
317 gilbertke 110 if (use_solv()) esolv1(natom,cutoffs.chrgcut,skip,atom.use,atom.charge,atom.x,atom.y,atom.z,&energies.esolv,deriv.desolv,deriv.drb);
318 gilbertke 103
319     energies.total = energies.estr + energies.ebend + energies.etor + energies.estrbnd + energies.e14+
320     energies.evdw + energies.eu + energies.ehbond + energies.eangang + energies.estrtor +
321     energies.eimprop + energies.eimptors + energies.eopb + energies.eurey + energies.esolv + energies.egeom;
322     for (i=1; i <= natom; i++)
323     {
324     for (j=0; j < 3; j++)
325     {
326     deriv.d1[i][j] = deriv.deb[i][j] + deriv.dea[i][j] + deriv.deaa[i][j] +
327     deriv.destb[i][j] + deriv.detor[i][j] + deriv.deopb[i][j] + deriv.dehb[i][j] +
328     deriv.destor[i][j] + deriv.deqq[i][j] + deriv.devdw[i][j] + deriv.de14[i][j] +
329     deriv.deimprop[i][j] + deriv.deub[i][j] + deriv.desolv[i][j] + deriv.degeom[i][j];
330    
331     }
332     }
333     }
334     // ==================================================================
335     double get_total_deriv_x(int i)
336     {
337     return (deriv.d1[i][0]);
338     }
339     double get_total_deriv_y(int i)
340     {
341     return (deriv.d1[i][1]);
342     }
343     double get_total_deriv_z(int i)
344     {
345     return (deriv.d1[i][2]);
346     }
347     // ============================================
348     void print_energy_totals()
349     {
350     printf("\nEnergy : %f\n",energies.total);
351     printf("Str: %f Bnd: %f Tor: %f\n", energies.estr, energies.ebend, energies.etor);
352     printf("StrBnd: %f VDW: %f QQ: %f\n",energies.estrbnd,energies.evdw+energies.e14+energies.ehbond, energies.eu);
353     printf("OOP: %f AA: %f Strtor: %f\n",energies.eopb,energies.eangang, energies.estrtor);
354     if (use_solv()) printf("Solv: %f \n",energies.esolv);
355     }
356     // ==================================================================
357     double get_total_energy()
358     {
359     return energies.total;
360     }
361     // ==================================================================
362     double energy()
363     {
364     double etot;
365     int i;
366    
367     energies.total = 0.0;
368     energies.estr = 0.0;
369     energies.ebend = 0.0;
370     energies.etor = 0.0;
371     energies.estrbnd = 0.0;
372     energies.evdw = 0.0;
373     energies.e14 = 0.0;
374     energies.ehbond = 0.0;
375     energies.eu = 0.0;
376     energies.eangang = 0.0;
377     energies.estrtor = 0.0;
378     energies.eimprop = 0.0;
379     energies.eimptors = 0.0;
380     energies.eopb = 0.0;
381     energies.eurey = 0.0;
382     energies.esolv = 0.0;
383     energies.egeom = 0.0;
384    
385    
386 gilbertke 114 if (use_bond())ebond(bonds_ff.nbnd,bonds_ff.i12,atom.use,atom.x,atom.y,atom.z,bonds_ff.bl,bonds_ff.bk,&energies.estr);
387 gilbertke 103 if (use_angle())eangle(angles.nang,angles.i13,atom.use,atom.x,atom.y,atom.z,angles.angin,angles.angtype,angles.anat,angles.acon,&energies.ebend);
388     if (use_tors())etorsion(torsions.ntor,torsions.i14,atom.use,atom.type,atom.x,atom.y,atom.z,torsions.ph1,torsions.ph2,torsions.ph3,torsions.ph4,torsions.ph5,
389     torsions.ph6,torsions.v1,torsions.v2,torsions.v3,torsions.v4,torsions.v5,torsions.v6,&energies.etor);
390 gilbertke 114 // mmff
391 gilbertke 103 if (use_strbnd())estrbnd(strbnd.nstrbnd,strbnd.isb,angles.i13,atom.use,atom.x,atom.y,atom.z,angles.anat,bonds_ff.bl,bonds_ff.bk,strbnd.ksb1,strbnd.ksb2,
392     &energies.estrbnd);
393 gilbertke 114 if (use_opbend_wilson()) eopbend_wilson(angles.nang,angles.i13,atom.use,atom.x,atom.y,atom.z,angles.angin,angles.copb,&energies.eopb);
394 gilbertke 103 if (use_hal()) ehal(natom,atom.type,atom.use,atom.x,atom.y,atom.z,cutoffs.vdwcut,skip,nonbond.vrad,nonbond.veps,&energies.evdw,&energies.e14);
395     if (use_bufcharge()) ebufcharge(natom,atom.use,skip,atom.x,atom.y,atom.z,atom.charge,cutoffs.chrgcut,&energies.eu);
396    
397 gilbertke 114 // gaff
398     if (use_imptor()) eimptors(improp.nimptors,improp.iiprop,atom.use,atom.type,atom.x,atom.y,atom.z,improp.v1,improp.v2,improp.v3,improp.ph1,improp.ph2,improp.ph3,&energies.eimptors);
399     if (use_lj()) elj(natom,atom.type,atom.use,atom.x,atom.y,atom.z,cutoffs.vdwcut,skip,nonbond.vrad,nonbond.veps,&energies.evdw,&energies.e14);
400     // if (use_charge())
401    
402 gilbertke 103 if (use_geom()) egeom(natom,atom.use,atom.x,atom.y,atom.z,&energies.egeom);
403 gilbertke 110 if (use_solv()) esolv(natom,cutoffs.chrgcut,skip,atom.use,atom.charge,atom.x,atom.y,atom.z,&energies.esolv);
404 gilbertke 103
405     energies.total = energies.estr + energies.ebend + energies.etor + energies.estrbnd
406     + energies.evdw + energies.e14 + energies.ehbond + energies.eu + energies.eangang +
407     energies.estrtor + energies.eimprop + energies.eimptors + energies.eopb + energies.eurey + energies.esolv + energies.egeom;
408     etot = energies.total;
409     return (etot);
410     }
411     /* =========================================================================== */
412     double print_energy()
413     {
414     double etot;
415     int i;
416    
417     energies.total = 0.0;
418     energies.estr = 0.0;
419     energies.ebend = 0.0;
420     energies.etor = 0.0;
421     energies.estrbnd = 0.0;
422     energies.evdw = 0.0;
423     energies.e14 = 0.0;
424     energies.ehbond = 0.0;
425     energies.eu = 0.0;
426     energies.eangang = 0.0;
427     energies.estrtor = 0.0;
428     energies.eimprop = 0.0;
429     energies.eimptors = 0.0;
430     energies.eopb = 0.0;
431     energies.eurey = 0.0;
432     energies.esolv = 0.0;
433     energies.egeom = 0.0;
434    
435    
436 gilbertke 114 if (use_bond())ebond(bonds_ff.nbnd,bonds_ff.i12,atom.use,atom.x,atom.y,atom.z,bonds_ff.bl,bonds_ff.bk,&energies.estr);
437 gilbertke 103 if (use_angle())eangle(angles.nang,angles.i13,atom.use,atom.x,atom.y,atom.z,angles.angin,angles.angtype,angles.anat,angles.acon,&energies.ebend);
438     if (use_tors())etorsion(torsions.ntor,torsions.i14,atom.use,atom.type,atom.x,atom.y,atom.z,torsions.ph1,torsions.ph2,torsions.ph3,torsions.ph4,torsions.ph5,
439     torsions.ph6,torsions.v1,torsions.v2,torsions.v3,torsions.v4,torsions.v5,torsions.v6,&energies.etor);
440 gilbertke 114 // mmff
441     if (use_opbend_wilson()) eopbend_wilson(angles.nang,angles.i13,atom.use,atom.x,atom.y,atom.z,angles.angin,angles.copb,&energies.eopb);
442 gilbertke 103 if (use_strbnd())estrbnd(strbnd.nstrbnd,strbnd.isb,angles.i13,atom.use,atom.x,atom.y,atom.z,angles.anat,bonds_ff.bl,bonds_ff.bk,strbnd.ksb1,strbnd.ksb2,
443     &energies.estrbnd);
444     if (use_hal()) ehal(natom,atom.type,atom.use,atom.x,atom.y,atom.z,cutoffs.vdwcut,skip,nonbond.vrad,nonbond.veps,&energies.evdw,&energies.e14);
445     if (use_bufcharge()) ebufcharge(natom,atom.use,skip,atom.x,atom.y,atom.z,atom.charge,cutoffs.chrgcut,&energies.eu);
446 gilbertke 114
447     // gaff
448     if (use_imptor()) eimptors(improp.nimptors,improp.iiprop,atom.use,atom.type,atom.x,atom.y,atom.z,improp.v1,improp.v2,improp.v3,improp.ph1,improp.ph2,improp.ph3,&energies.eimptors);
449     if (use_lj()) elj(natom,atom.type,atom.use,atom.x,atom.y,atom.z,cutoffs.vdwcut,skip,nonbond.vrad,nonbond.veps,&energies.evdw,&energies.e14);
450     // if (use_charge())
451 gilbertke 103
452     if (use_geom()) egeom(natom,atom.use,atom.x,atom.y,atom.z,&energies.egeom);
453 gilbertke 110 if (use_solv()) esolv(natom,cutoffs.chrgcut,skip,atom.use,atom.charge,atom.x,atom.y,atom.z,&energies.esolv);
454 gilbertke 103
455     energies.total = energies.estr + energies.ebend + energies.etor + energies.estrbnd
456     + energies.evdw + energies.e14 + energies.ehbond + energies.eu + energies.eangang +
457     energies.estrtor + energies.eimprop + energies.eimptors + energies.eopb + energies.eurey + energies.esolv + energies.egeom;
458     etot = energies.total;
459    
460     return (etot);
461     }
462     // ====================== hessian =======================
463     void hessian(int maxhess, double *h, int *hinit, int *hstop,int *hindex, double *hdiag)
464     {
465     int i,j,k,nhess;
466     double hmax;
467     char hatext[60];
468    
469     int *keep; // keep[MAXATOM];
470    
471     keep = ivector(0,natom);
472    
473     nhess = 0;
474     for (i=1; i <= natom; i++)
475     {
476    
477     hinit[nhess] = 0;
478     hstop[nhess] = 0;
479     hdiag[nhess] = 0.0;
480     nhess++;
481     hinit[nhess] = 0;
482     hstop[nhess] = 0;
483     hdiag[nhess] = 0.0;
484     nhess++;
485     hinit[nhess] = 0;
486     hstop[nhess] = 0;
487     hdiag[nhess] = 0.0;
488     nhess++;
489     }
490     /* periodic boundary conditions set here */
491    
492     // if (use_picalc) piseq(0,0);
493    
494     nhess= 0;
495     for (i=1; i <= natom; i++)
496     {
497     for(k=1; k <= natom; k++)
498     {
499     for(j=0; j < 3; j++)
500     {
501     hess.hessx[k][j] = 0.0;
502     hess.hessy[k][j] = 0.0;
503     hess.hessz[k][j] = 0.0;
504     }
505     }
506    
507     if (atom.use[i])
508     {
509     if (use_bond())ebond2(i,bonds_ff.i12,atom.iat,atom.bo,atom.x,atom.y,atom.z,bonds_ff.bl,bonds_ff.bk,hess.hessx,hess.hessy,hess.hessz);
510     if (use_angle())eangle2(i,angles.nang,angles.i13,angles.angin,atom.x,atom.y,atom.z,angles.acon,angles.anat,deriv.dea,hess.hessx,hess.hessy,hess.hessz);
511 gilbertke 114 if (use_tors())etorsion2(i,torsions.ntor,torsions.i14,atom.use,atom.x,atom.y,atom.z,torsions.ph1,torsions.ph2,torsions.ph3,torsions.ph4,torsions.ph5,
512     torsions.ph6,torsions.v1,torsions.v2,torsions.v3,torsions.v4,torsions.v5,torsions.v6,hess.hessx,hess.hessy,hess.hessz);
513    
514     // mmff
515 gilbertke 103 if (use_strbnd())estrbnd2(i,strbnd.nstrbnd,strbnd.isb,angles.i13,atom.use,atom.x,atom.y,atom.z,angles.anat,bonds_ff.bl,bonds_ff.bk,strbnd.ksb1,strbnd.ksb2,
516     hess.hessx,hess.hessy,hess.hessz);
517     if (use_opbend_wilson()) eopbend_wilson2(i,angles.nang,angles.i13,atom.use,atom.x,atom.y,atom.z,angles.angin,angles.copb,hess.hessx,hess.hessy,hess.hessz);
518     if (use_hal())ehal2(i,natom,atom.type,atom.use,atom.x,atom.y,atom.z,cutoffs.vdwcut,skip,nonbond.vrad,nonbond.veps,hess.hessx,hess.hessy,hess.hessz);
519     if (use_bufcharge())ebufcharge2(i,natom,atom.use,skip,atom.x,atom.y,atom.z,atom.charge,cutoffs.chrgcut,hess.hessx,hess.hessy,hess.hessz);
520    
521 gilbertke 114 // gaff
522     if (use_imptor()) eimptors2(i,improp.nimptors,improp.iiprop,atom.use,atom.type,atom.x,atom.y,atom.z,improp.v1,improp.v2,improp.v3,improp.ph1,improp.ph2,improp.ph3,hess.hessx,
523     hess.hessy,hess.hessz);
524     if (use_lj()) elj2(i,natom,atom.type,atom.use,atom.x,atom.y,atom.z,cutoffs.vdwcut,skip,nonbond.vrad,nonbond.veps,hess.hessx,hess.hessy,hess.hessz);
525     //if (use_charge())
526    
527 gilbertke 103 if (use_geom()) egeom2(i,natom,atom.use,atom.x,atom.y,atom.z,hess.hessx,hess.hessy,hess.hessz);
528 gilbertke 110 if (use_solv()) esolv2(i,natom,cutoffs.chrgcut,atom.charge,atom.x,atom.y,atom.z,hess.hessx,hess.hessy,hess.hessz);
529 gilbertke 103
530     hdiag[(i-1)*3] += hess.hessx[i][0];
531     hdiag[(i-1)*3+1] += hess.hessy[i][1];
532     hdiag[(i-1)*3+2] += hess.hessz[i][2];
533    
534     for (k=i+1; k <= natom; k++)
535     {
536     keep[k] = FALSE;
537     if (atom.use[k])
538     {
539     hmax = fabs(hess.hessx[k][0]);
540     for (j=0; j < 3;j++)
541     {
542     if (fabs(hess.hessx[k][j]) > hmax)
543     hmax = fabs(hess.hessx[k][j]);
544     }
545     for (j=0; j < 3;j++)
546     {
547     if (fabs(hess.hessy[k][j]) > hmax)
548     hmax = fabs(hess.hessy[k][j]);
549     }
550     for (j=0; j < 3;j++)
551     {
552     if (fabs(hess.hessz[k][j]) > hmax)
553     hmax = fabs(hess.hessz[k][j]);
554     }
555     if (hmax >= hesscut) keep[k] = TRUE;
556     }
557     }
558    
559     hinit[(i-1)*3] = nhess; // x component of new atom
560    
561     hindex[nhess] = 3*(i-1) + 1;
562     h[nhess] = hess.hessx[i][1];
563     nhess++;
564     hindex[nhess] = 3*(i-1) + 2;
565     h[nhess] = hess.hessx[i][2];
566     nhess++;
567    
568     for (k=i+1; k <= natom; k++)
569     {
570     if (keep[k])
571     {
572     for (j=0; j < 3; j++)
573     {
574     hindex[nhess] = 3*(k-1) +j;
575     h[nhess] = hess.hessx[k][j];
576     nhess++;
577     }
578     }
579     }
580     hstop[(i-1)*3] = nhess;
581    
582     hinit[(i-1)*3+1] = nhess; // y component of atom i
583     hindex[nhess] = 3*(i-1) + 2;
584     h[nhess] = hess.hessy[i][2];
585     nhess++;
586     for (k=i+1; k <= natom; k++)
587     {
588     if (keep[k])
589     {
590     for (j=0; j < 3; j++)
591     {
592     hindex[nhess] = 3*(k-1) +j;
593     h[nhess] = hess.hessy[k][j];
594     nhess++;
595     }
596     }
597     }
598     hstop[(i-1)*3+1] = nhess;
599     hinit[(i-1)*3+2] = nhess;
600     for (k=i+1; k <= natom; k++)
601     {
602     if (keep[k])
603     {
604     for (j=0; j < 3; j++)
605     {
606     hindex[nhess] = 3*(k-1) +j;
607     h[nhess] = hess.hessz[k][j];
608     nhess++;
609     }
610     }
611     }
612     hstop[(i-1)*3+2] = nhess;
613    
614     if (nhess > maxhess)
615     {
616     sprintf(hatext,"Too many hessian elements: %d elements %d max\n",nhess,maxhess);
617     message_alert(hatext,"Error");
618     fprintf(pcmlogfile,"Too many hessian elements !\n");
619     strcpy(Openbox.fname,"Hesserr.pcm");
620     exit(0);
621     }
622     }
623     }
624     free_ivector(keep, 0,natom);
625     }
626