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