ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/freemol/trunk/src/mengine/src/pcm7.c
Revision: 103
Committed: Thu Feb 19 01:37:38 2009 UTC (12 years, 7 months ago) by gilbertke
File size: 24129 byte(s)
Log Message:
major rewrite - removing global data, adding electrostatics tag to read_sdf
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    
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 (*)[MAXIAT],int (*)[MAXIAT],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