ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/freemol/trunk/src/mengine/src/pcm7.c
(Generate patch)
# Line 1 | Line 1
1 < #define EXTERN extern
2 <
3 < #include "pcwin.h"
4 < #include "pcmod.h"
5 <
6 < #include "attached.h"
7 < #include "energies.h"
8 < #include "bonds_ff.h"
9 < #include "derivs.h"
10 < #include "pot.h"
11 < #include "utility.h"
12 < #include "solv.h"  
13 < #include "dipmom.h"
14 <                
15 < EXTERN struct t_minim_control {
16 <        int type, method, field, added_const;
17 <        char added_path[256],added_name[256];
18 <        } minim_control;
19 <        
20 < EXTERN struct t_minim_values {
21 <        int iprint, ndc, nconst;
22 <        float dielc;
23 <        } minim_values;
24 < struct t_solvent {
25 <    int type;
26 <    double doffset, p1,p2,p3,p4,p5;
27 <    double *shct,*asolv,*rsolv,*vsolv,*gpol,*rborn;
28 <    } solvent;
29 <
30 < int Missing_constants;
31 <    
32 < void message_alert(char *,char *);
33 < void pcm7_min(void);
34 < int isbond(int, int);
35 < int initialise(void);
36 < void get_bonds(void);
37 < void get_torsions(void);
38 < void set_field(void);
39 < double energy(void);
40 < double print_energy(void);
41 < void pcmfout(int);
42 < void set_active(void);
43 < int setup_calculation(void);
44 < void end_calculation(void);
45 <
46 < void egeom1(void);
47 < void esolv1(void);
48 < void egeom(void);
49 < void esolv(void);
50 <
51 < void ehal(void);
52 < void ehal1(void);
53 < void ebufcharge(void);
54 < void ebufcharge1(void);
55 < void eopbend_wilson(void);
56 < void eopbend_wilson1(void);
57 <
58 < void estrbnd(void);
59 < void etorsion(void);
60 < void eangang(void);
61 < void eangle(void);
62 < void ebond(void);
63 <
64 < void estrbnd1(void);
65 < void etorsion1(void);
66 < void eangle1(void);
67 < void ebond1(void);
68 <
69 < void kangle(void);
70 < void ktorsion(void);
71 < void kstrbnd(void);
72 < void kcharge(void);
73 < void ksolv(void);
74 < void kopbend(void);
75 < void kvdw(void);
76 < int  kbond(void);
77 <
78 < void get_memory(void);
79 < void free_memory(void);
80 < void gradient(void);
81 < void minimize(void);
82 < void read_datafiles(char *);
83 < void get_added_const(void);
84 < void attach(void);
85 < void hdel(int);
86 < void set_atomtypes(int);
87 < void type(void);
88 < void zero_data(void);
89 < void xlogp(float *);
90 <
91 < // ==================================================================
92 < int setup_calculation()
93 < {
94 <  int i,j,nRet;
95 <    char string[30];
96 <    double etot;
97 <    if (minim_control.field == MMFF94)
98 <    {
99 <        // copy mm3 types into type file and retype
100 <        pot.use_hbond = FALSE;
101 <        pot.use_picalc = FALSE;
102 <        hdel(1);
103 <        set_field();
104 <        set_atomtypes(MMFF94);
105 <        if (minim_control.added_const)
106 <           get_added_const();
107 <    } else
108 <    {
109 <        if (minim_control.added_const)
110 <        {
111 <           strcpy(string,"mmxconst.prm");
112 <           zero_data();
113 <           read_datafiles(string);
114 <           get_added_const();
115 <        }
116 <        type();
117 <        set_field();
118 <    }
119 <
120 <            
121 <  // allocate memeory for derivatives
122 <  get_memory();
123 <
124 <  // setup bonds, angles, torsions, improper torsions and angles, allenes
125 <  get_bonds();
126 <  get_angles();
127 <  get_torsions();
128 <
129 <  attach();
130 <  // need allene
131 <
132 <  // set active atoms
133 <  set_active();
134 <  // setup non_bonded list of atoms to skip
135 <  for (i=1; i <= natom; i++)
136 <    {
137 <         for (j=0; j < MAXIAT; j++)
138 <             if (atom[i].iat[j] != 0 && atom[i].bo[j] != 9)
139 <               {
140 <                    skip[i][atom[i].iat[j]] = i;
141 <                    skip[atom[i].iat[j]][i] = atom[i].iat[j];
142 <               }
143 <         for (j=0; j < attached.n13[i]; j++)
144 <           {
145 <               skip[i][attached.i13[j][i]] = i;
146 <               skip[attached.i13[j][i]][i] = attached.i13[j][i];
147 <           }
148 <         for(j=0; j < attached.n14[i]; j++)
149 <           {
150 <            skip[i][attached.i14[j][i]] = -i;
151 <            skip[attached.i14[j][i]][i] = -attached.i14[j][i];
152 <           }
153 <    }
154 < /*  assign local geometry potential function parameters  */
155 <    Missing_constants = FALSE;
156 <     if (pot.use_bond || pot.use_strbnd)  nRet = kbond();
157 <     if (nRet == FALSE) // use end_calc to free memory
158 <     {
159 <         energies.total = 9999.;
160 <         return FALSE;
161 <     }
162 <     if (pot.use_angle || pot.use_strbnd || pot.use_angang) kangle();
163 <     if (pot.use_angle || pot.use_opbend || pot.use_opbend_wilson) kopbend();
164 <     if (pot.use_tors  || pot.use_strtor || pot.use_tortor) ktorsion();
165 <     if (pot.use_strbnd) kstrbnd();
166 <      
167 <      if (pot.use_buck || pot.use_hal || pot.use_lj) kvdw();
168 <      if ((pot.use_charge || pot.use_bufcharge || pot.use_ewald)) kcharge();
169 <
170 <       if (Missing_constants == TRUE)
171 <       {
172 <           energies.total = 9999.0;
173 <           return (FALSE);
174 <       }    
175 <       if (minim_values.iprint == TRUE)
176 <          etot = print_energy();
177 <       else
178 <          etot = energy();
179 <       return TRUE;
180 < }
181 < // =================================================================
182 < void end_calculation()
183 < {
184 <
185 <    free_memory();
186 < }
187 < // ==================================================================
188 < void gradient()
189 < {
190 <    int i, j;
191 <
192 <      energies.estr = 0.0;
193 <      energies.ebend = 0.0;
194 <      energies.estrbnd = 0.0;
195 <      energies.e14 = 0.0;
196 <      energies.evdw = 0.0;
197 <      energies.etor = 0.0;
198 <      energies.eu = 0.0;
199 <      energies.eopb = 0.0;
200 <      energies.eangang = 0.0;
201 <      energies.estrtor = 0.0;
202 <      energies.ehbond = 0.0;
203 <      energies.efix = 0.0;
204 <      energies.eimprop = 0.0;
205 <      energies.eimptors = 0.0;
206 <      energies.eurey = 0.0;
207 <      energies.esolv = 0.0;
208 <      energies.egeom =0.0;
209 <      
210 <      for (i=1; i <= natom; i++)
211 <           atom[i].energy = 0.0;
212 <      
213 <      for (i=1; i <= natom; i++)
214 <      {
215 <        for (j=0; j < 3; j++)
216 <        {
217 <            deriv.deb[i][j] = 0.0;
218 <            deriv.dea[i][j] = 0.0;
219 <            deriv.destb[i][j] = 0.0;
220 <            deriv.deopb[i][j] = 0.0;
221 <            deriv.detor[i][j] = 0.0;
222 <            deriv.de14[i][j] = 0.0;
223 <            deriv.devdw[i][j]= 0.0;
224 <            deriv.deqq[i][j] = 0.0;
225 <            deriv.deaa[i][j] = 0.0;
226 <            deriv.destor[i][j] = 0.0;
227 <            deriv.deimprop[i][j] = 0.0;
228 <            deriv.dehb[i][j] = 0.0;
229 <            deriv.desolv[i][j] = 0.0;
230 <            deriv.degeom[i][j] = 0.0;
231 <        }
232 <      }
233 <  
234 <      if (pot.use_bond)ebond1();
235 <      if (pot.use_angle)eangle1();
236 <      if (pot.use_opbend_wilson) eopbend_wilson1();
237 <      if (pot.use_tors)etorsion1();
238 <      if (pot.use_strbnd)estrbnd1();
239 <    
240 <      // mmff
241 <      if (pot.use_hal) ehal1();
242 <      if (pot.use_bufcharge) ebufcharge1();
243 <
244 <      if (pot.use_geom) egeom1();
245 <    
246 <      energies.total =  energies.estr + energies.ebend + energies.etor + energies.estrbnd + energies.e14+
247 <                        energies.evdw + energies.eu + energies.ehbond + energies.eangang + energies.estrtor +
248 <                        energies.eimprop + energies.eimptors + energies.eopb + energies.eurey + energies.esolv + energies.egeom;
249 <      for (i=1; i <= natom; i++)
250 <      {
251 <          for (j=0; j < 3; j++)
252 <          {
253 <              deriv.d1[i][j] = deriv.deb[i][j] + deriv.dea[i][j] + deriv.deaa[i][j] +
254 <                               deriv.destb[i][j] + deriv.detor[i][j] + deriv.deopb[i][j] + deriv.dehb[i][j] +
255 <                               deriv.destor[i][j] + deriv.deqq[i][j] + deriv.devdw[i][j] + deriv.de14[i][j] +
256 <                               deriv.deimprop[i][j] + deriv.deub[i][j] + deriv.desolv[i][j] + deriv.degeom[i][j];
257 <
258 <         }
259 <      }
260 < }
261 < // ==================================================================
262 < double energy()
263 < {
264 <    double etot;
265 <    int i;
266 <    
267 <      energies.total = 0.0;
268 <      energies.estr = 0.0;
269 <      energies.ebend = 0.0;
270 <      energies.etor = 0.0;
271 <      energies.estrbnd = 0.0;
272 <      energies.evdw = 0.0;
273 <      energies.e14 = 0.0;
274 <      energies.ehbond = 0.0;
275 <      energies.eu = 0.0;
276 <      energies.eangang = 0.0;
277 <      energies.estrtor = 0.0;
278 <      energies.eimprop = 0.0;
279 <      energies.eimptors = 0.0;
280 <      energies.eopb = 0.0;
281 <      energies.eurey = 0.0;
282 <      energies.esolv = 0.0;
283 <      energies.egeom = 0.0;
284 <
285 <      for (i=1; i <= natom; i++)
286 <           atom[i].energy = 0.0;
287 <
288 <      if (pot.use_bond)ebond();
289 <      if (pot.use_angle)eangle();
290 <      if (pot.use_opbend_wilson) eopbend_wilson();
291 <      if (pot.use_tors)etorsion();
292 <      if (pot.use_strbnd)estrbnd();
293 <      // mmff
294 <      if (pot.use_hal) ehal();
295 <      if (pot.use_bufcharge) ebufcharge();
296 <
297 <      if (pot.use_geom) egeom();
298 <                              
299 <      energies.total =  energies.estr + energies.ebend + energies.etor + energies.estrbnd
300 <             + energies.evdw + energies.e14 + energies.ehbond + energies.eu + energies.eangang +
301 <               energies.estrtor + energies.eimprop + energies.eimptors + energies.eopb + energies.eurey + energies.esolv + energies.egeom;
302 <      etot = energies.total;
303 <      return (etot);
304 < }
305 < /* =========================================================================== */
306 < double print_energy()
307 < {
308 <    double etot;
309 <    int i;
310 <    
311 <      energies.total = 0.0;
312 <      energies.estr = 0.0;
313 <      energies.ebend = 0.0;
314 <      energies.etor = 0.0;
315 <      energies.estrbnd = 0.0;
316 <      energies.evdw = 0.0;
317 <      energies.e14 = 0.0;
318 <      energies.ehbond = 0.0;
319 <      energies.eu = 0.0;
320 <      energies.eangang = 0.0;
321 <      energies.estrtor = 0.0;
322 <      energies.eimprop = 0.0;
323 <      energies.eimptors = 0.0;
324 <      energies.eopb = 0.0;
325 <      energies.eurey = 0.0;
326 <      energies.esolv = 0.0;
327 <      energies.egeom = 0.0;
328 <
329 <      for (i=1; i <= natom; i++)
330 <           atom[i].energy = 0.0;
331 <
332 <      if (pot.use_bond)ebond();
333 <      if (pot.use_angle)eangle();
334 <      if (pot.use_opbend_wilson) eopbend_wilson();
335 <      if (pot.use_tors)etorsion();
336 <      if (pot.use_strbnd)estrbnd();
337 <      // mmff
338 <      if (pot.use_hal) ehal();
339 <      if (pot.use_bufcharge) ebufcharge();
340 <      
341 <      if (pot.use_geom) egeom();
342 <
343 <      energies.total =  energies.estr + energies.ebend + energies.etor + energies.estrbnd
344 <             + energies.evdw + energies.e14 + energies.ehbond + energies.eu + energies.eangang +
345 <               energies.estrtor + energies.eimprop + energies.eimptors + energies.eopb + energies.eurey + energies.esolv + energies.egeom;
346 <      etot = energies.total;
347 <      
348 <      return (etot);
349 < }    
350 <
351 < /* ================================================================== */
352 < int isbond(int i, int j)
353 < {
354 <    int k;
355 <    for (k=0; k < MAXIAT; k++)
356 <    {
357 <        if (atom[i].iat[k] == j && atom[i].bo[k] != 9)
358 <            return TRUE;
359 <    }
360 <    return FALSE;
361 < }
362 < // ==========================
363 < void get_bonds()
364 < {
365 <    int i, j;
366 <
367 <    bonds_ff.nbnd = 0;    
368 <    for (i=1; i <= natom; i++)
369 <    {
370 <        for (j=i+1; j <= natom; j++)
371 <        {
372 <           if (isbond(i,j))
373 <           {
374 <              bonds_ff.i12[bonds_ff.nbnd][0] = i;
375 <              bonds_ff.i12[bonds_ff.nbnd][1] = j;
376 <              bonds_ff.nbnd++;
377 <              if (bonds_ff.nbnd > MAXBND)
378 <              {
379 <                 fprintf(pcmlogfile,"Error - Too many bonds!\nProgram will now quit.\n");
380 <                 exit(0);
381 <               }
382 <           }
383 <        }
384 <    }
385 < }
386 <
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 (*)[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 >

Diff Legend

Removed lines
+ Added lines
< Changed lines
> Changed lines