ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/freemol/trunk/src/mengine/src/pcm7.c
Revision: 115
Committed: Thu May 7 14:35:26 2009 UTC (11 years, 2 months ago) by gilbertke
File size: 28363 byte(s)
Log Message:
added charge-charge routines
Line File contents
1 #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 #include "fix.h"
17
18 EXTERN struct t_strbnd {
19 int nstrbnd, **isb;
20 float *ksb1, *ksb2;
21 } strbnd;
22 EXTERN struct t_improp {
23 int nimptors, **iiprop;
24 float *v1, *v2, *v3;
25 int *ph1, *ph2, *ph3;
26 } improp;
27
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 void get_torsions(int natom,int *type,int **iat,int **bo);
54 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 void set_active(int natom,int *use,int natom_fix,int *katm_fix);
62 int setup_calculation(void);
63 void end_calculation(void);
64
65 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 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 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 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 float *ksb2,double *estrbnd);
82 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 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 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
89 void hessian(int, double *,int *, int *, int *,double *);
90 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 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 void eopbend_wilson2(int i,int nang,int **,int *use,double *x,double *y,double *z,int *angin,float *copb,
98 float **hessx,float **hessy,float **hessz);
99 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 float *ksb2,float **hessx,float **hessy,float **hessz);
101 void esolv2(int ia,int natom,double chrgcut,double *charge,double *x,double *y,double *z,float **hessx,float **hessy,float **hessz);
102 void egeom2(int ,int natom,int *use,double *x,double *y,double *z,float **hessx,float **hessy,float **hessz);
103 // 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(int natom,int *use, int **skip,double *x,double *y,double *z,double *charge,double chrgcut,double *eu);
115 void echarge1(int natom,int *use, int **skip,double *x,double *y,double *z,double *charge,double chrgcut,double *eu,double **deqq);
116 void echarge2(int i,int natom,int *use, int **skip,double *x,double *y,double *z,double *charge,double chrgcut,float **hessx,float **hessy,float **hessz);
117 //
118
119 int kbond(void);
120 void kangle(void);
121 void kstrbnd(void);
122 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
128 // gaff
129 void kimptors(void);
130
131 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 int use_imptor(void);
153 int use_lj(void);
154
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 } else if (minim_control.field == GAFF)
170 {
171 set_field(GAFF);
172 set_atomtypes(GAFF);
173 } else
174 {
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 get_torsions(natom,atom.type,atom.iat,atom.bo);
194
195 attach();
196 // need allene
197
198 // set active atoms
199 set_active(natom,atom.use,fx_atom.natom_fix,fx_atom.katom_fix);
200 // 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 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 if (use_strbnd()) kstrbnd();
232 if (use_imptor()) kimptors();
233
234 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
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 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 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 // mmff
305 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 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 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 // 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()) echarge1(natom,atom.use,skip,atom.x,atom.y,atom.z,atom.charge,cutoffs.chrgcut,&energies.eu,deriv.deqq);
315
316 if (use_geom()) egeom1(natom,atom.use,atom.x,atom.y,atom.z,&energies.egeom,deriv.degeom);
317 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
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 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 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 // mmff
391 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 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 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 // 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()) echarge(natom,atom.use,skip,atom.x,atom.y,atom.z,atom.charge,cutoffs.chrgcut,&energies.eu);
401
402 if (use_geom()) egeom(natom,atom.use,atom.x,atom.y,atom.z,&energies.egeom);
403 if (use_solv()) esolv(natom,cutoffs.chrgcut,skip,atom.use,atom.charge,atom.x,atom.y,atom.z,&energies.esolv);
404
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 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 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 // 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 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
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()) echarge(natom,atom.use,skip,atom.x,atom.y,atom.z,atom.charge,cutoffs.chrgcut,&energies.eu);
451
452 if (use_geom()) egeom(natom,atom.use,atom.x,atom.y,atom.z,&energies.egeom);
453 if (use_solv()) esolv(natom,cutoffs.chrgcut,skip,atom.use,atom.charge,atom.x,atom.y,atom.z,&energies.esolv);
454
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 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 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 // 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()) echarge2(i,natom,atom.use,skip,atom.x,atom.y,atom.z,atom.charge,cutoffs.chrgcut,hess.hessx,hess.hessy,hess.hessz);
526
527 if (use_geom()) egeom2(i,natom,atom.use,atom.x,atom.y,atom.z,hess.hessx,hess.hessy,hess.hessz);
528 if (use_solv()) esolv2(i,natom,cutoffs.chrgcut,atom.charge,atom.x,atom.y,atom.z,hess.hessx,hess.hessy,hess.hessz);
529
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