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