ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/freemol/trunk/src/mengine/src/ehal.c
Revision: 87
Committed: Mon Jan 12 17:11:11 2009 UTC (11 years, 10 months ago) by gilbertke
File size: 11118 byte(s)
Log Message:
fixed corruption in mengine.c and bug in ehal.c
Line User Rev File contents
1 tjod 3 #define EXTERN extern
2    
3     #include "pcwin.h"
4     #include "pcmod.h"
5     #include "utility.h"
6    
7     #include "energies.h"
8     #include "derivs.h"
9     #include "hess.h"
10     #include "torsions.h"
11     #include "nonbond.h"
12     #include "field.h"
13     #include "cutoffs.h"
14     #include "attached.h"
15     #include "pot.h"
16     #include "solv.h"
17    
18     EXTERN struct t_minim_values {
19     int iprint, ndc, nconst;
20     float dielc;
21     } minim_values;
22    
23     int ishbond(int, int);
24     int iscoord_bond(int, int);
25    
26     EXTERN int *skip;
27    
28     void image(double *, double *, double *, int);
29    
30     void ehal()
31     {
32     int i,ia, ib, ita, itb, it, j, k;
33     double xi,yi,zi,xr,yr,zr;
34     double xk,yk,zk;
35     double e,rv,eps;
36     double rik,rik2,rik7;
37     double sigma, kappa, kappa7;
38     double rv7,rho, tau;
39     double cutoff;
40    
41     cutoff = cutoffs.vdwcut*cutoffs.vdwcut;
42     energies.evdw = 0.0;
43     energies.e14 = 0.0;
44     sigma = 1.12;
45     kappa = 1.07;
46     kappa7 = pow(kappa,7.0);
47    
48     if (minim_values.iprint)
49     {
50 wdelano 58 fprintf(pcmlogfile,"\nVDW Terms - MMFF Buffered 14-7 Potential\n");
51     fprintf(pcmlogfile," At1 At2 Rik Radius Eps Evdw\n");
52 tjod 3 }
53     for (i=1; i <= natom; i++)
54     skip[i] = 0;
55     for (i=1; i < natom; i++)
56     {
57     for (j=0; j < MAXIAT; j++)
58     if (atom[i].iat[j] != 0 && atom[i].bo[j] != 9)
59     skip[atom[i].iat[j]] = i;
60     for (j=0; j < attached.n13[i]; j++)
61     skip[attached.i13[j][i]] = i;
62     for(k=0; k < attached.n14[i]; k++)
63     skip[attached.i14[k][i]] = -i;
64    
65     for(j=i+1; j <= natom; j++)
66     {
67     if (skip[j] != i)
68     {
69     ia = i;
70     ib = j;
71     if (atom[ia].use || atom[ib].use)
72     {
73     ita = atom[ia].type;
74     itb = atom[ib].type;
75     if (ita <= itb)
76     it = 1000*ita + itb;
77     else
78     it = 1000*itb + ita;
79    
80     xi = atom[ia].x;
81     yi = atom[ia].y;
82     zi = atom[ia].z;
83     xk = atom[ib].x;
84     yk = atom[ib].y;
85     zk = atom[ib].z;
86    
87     xr = xi - xk;
88     yr = yi - yk;
89     zr = zi - zk;
90     rik2 = xr*xr + yr*yr + zr*zr;
91     if (rik2 < cutoff)
92     {
93     rv = nonbond.vrad[ita][itb];
94     eps = nonbond.veps[ita][itb];
95     if (skip[j] == -i)
96     eps /= units.v14scale;
97     rik = sqrt(rik2);
98     rv7 = pow(rv,7.0);
99     rik7 = pow(rik,7.0);
100     rho = rik7 + (sigma-1.00)*rv7;
101     tau = rik + (kappa-1.00)*rv;
102    
103     e = eps*kappa7*(rv7/pow(tau,7.0))*(sigma*rv7/rho-2.00);
104     if (skip[j] == -i)
105     energies.e14 += e;
106     else
107     energies.evdw += e;
108     atom[ia].energy += e;
109     atom[ib].energy += e;
110     if (minim_values.iprint)
111 wdelano 58 fprintf(pcmlogfile,"VDW: %-4d - %-4d %-8.3f %-8.3f %-8.3f = %8.4f\n",ia, ib, sqrt(rik2),rv,eps,e);
112 tjod 3 }
113     }
114     }
115     }
116     }
117     }
118    
119     void ehal1()
120     {
121     int i,ia, ib, j, k, ita, itb, it;
122     double xi,yi,zi,xr,yr,zr;
123     double xk,yk,zk;
124     double e,rv,eps;
125     double rik,rik2;
126     double rik6,rik7;
127     double sigma, kappa, kappa7;
128 gilbertke 86 double rv7,rho, tau,tau7,tau8, rv14;
129 tjod 3 double dedx,dedy,dedz,de;
130     double cutoff;
131    
132     cutoff = cutoffs.vdwcut*cutoffs.vdwcut;
133    
134     energies.evdw = 0.0;
135     energies.e14 = 0.0;
136    
137     sigma = 1.12;
138     kappa = 1.07;
139 gilbertke 86 kappa7 = kappa*kappa*kappa*kappa*kappa*kappa*kappa;
140 tjod 3
141     for (i=1; i <= natom; i++)
142     {
143     deriv.devdw[i][0] = 0.0;
144     deriv.devdw[i][1] = 0.0;
145     deriv.devdw[i][2] = 0.0;
146     deriv.de14[i][0] = 0.0;
147     deriv.de14[i][1] = 0.0;
148     deriv.de14[i][2] = 0.0;
149     }
150    
151     for (i=1; i <= natom; i++)
152     skip[i] = 0;
153    
154     for (i=1; i < natom; i++)
155     {
156     for (j=0; j < MAXIAT; j++)
157     if (atom[i].iat[j] != 0 && atom[i].bo[j] != 9)
158     skip[atom[i].iat[j]] = i;
159     for (j=0; j < attached.n13[i]; j++)
160     skip[attached.i13[j][i]] = i;
161     for(k=0; k < attached.n14[i]; k++)
162     skip[attached.i14[k][i]] = -i;
163     ia = i;
164     ita = atom[i].type;
165    
166     for(j=i+1; j <= natom; j++)
167     {
168     if (atom[i].use || atom[j].use)
169     {
170     if (skip[j] != i)
171     {
172     ib = j;
173     itb = atom[j].type;
174     if (ita <= itb)
175     it = 1000*ita + itb;
176     else
177     it = 1000*itb + ita;
178    
179     xi = atom[ia].x;
180     yi = atom[ia].y;
181     zi = atom[ia].z;
182     xk = atom[ib].x;
183     yk = atom[ib].y;
184     zk = atom[ib].z;
185    
186     xr = xi - xk;
187     yr = yi - yk;
188     zr = zi - zk;
189     rik2 = xr*xr + yr*yr + zr*zr;
190     if (rik2 < cutoff)
191     {
192     rv = nonbond.vrad[ita][itb];
193     eps = nonbond.veps[ita][itb];
194     if (skip[j] == -i)
195     eps /= units.v14scale;
196     rik = sqrt(rik2);
197 gilbertke 87 rv7 = rv*rv*rv*rv*rv*rv*rv; // pow(rv,7.0);
198 tjod 3 rv14 = rv7*rv7;
199     rik6 = rik2*rik2*rik2;
200     rik7 = rik6*rik;
201     rho = rik7 + (sigma-1.00)*rv7;
202     tau = rik + (kappa-1.00)*rv;
203 gilbertke 87 tau7 = tau*tau*tau*tau*tau*tau*tau;
204     tau8 = tau7*tau;
205 gilbertke 86 e = eps*kappa7*(rv7/tau7)*(sigma*rv7/rho-2.00);
206 tjod 3
207 gilbertke 86 de = -7.00*eps*kappa7*(rv7/tau8)*(sigma*rv7/rho-2.00)
208     -7.00*eps*kappa7*sigma*rv14*rik6/(rho*rho*tau7);
209 tjod 3
210     de /= rik;
211     dedx = de*xr;
212     dedy = de*yr;
213     dedz = de*zr;
214    
215     if (skip[j] == -i)
216     {
217     energies.e14 += e;
218     deriv.de14[ia][0] += dedx;
219     deriv.de14[ia][1] += dedy;
220     deriv.de14[ia][2] += dedz;
221     deriv.de14[ib][0] -= dedx;
222     deriv.de14[ib][1] -= dedy;
223     deriv.de14[ib][2] -= dedz;
224     } else
225     {
226     energies.evdw += e;
227     deriv.devdw[ia][0] += dedx;
228     deriv.devdw[ia][1] += dedy;
229     deriv.devdw[ia][2] += dedz;
230     deriv.devdw[ib][0] -= dedx;
231     deriv.devdw[ib][1] -= dedy;
232     deriv.devdw[ib][2] -= dedz;
233     }
234     virial.virx += xr*dedx;
235     virial.viry += yr*dedy;
236     virial.virz += zr*dedz;
237     }
238     }
239     }
240     }
241     }
242     }
243    
244     void ehal2(int iatom)
245     {
246 gilbertke 86 int j,ia, ib, k, ita, itb, it;
247 tjod 3 double xi,yi,zi,xr,yr,zr;
248     double xk,yk,zk;
249 gilbertke 86 double rv,eps;
250 tjod 3 double rik,rik2, rik5;
251     double rik6,rik7, rik12;
252     double sigma, kappa, kappa7;
253     double rv7,rho, tau, rv14;
254     double de,d2e;
255     double d2edx,d2edy,d2edz,term[3][3];
256     double cutoff;
257    
258     cutoff = cutoffs.vdwcut*cutoffs.vdwcut;
259    
260     sigma = 1.12;
261     kappa = 1.07;
262     kappa7 = pow(kappa,7.0);
263    
264     for (j=1; j <= natom; j++)
265     skip[j] = 0;
266    
267     skip[iatom] = iatom;
268     ita = atom[iatom].type;
269     for (k=0; k < MAXIAT; k++)
270     {
271     if (atom[iatom].iat[k] != 0 && atom[iatom].bo[k] != 9)
272     skip[atom[iatom].iat[k]] = iatom;
273     }
274     for(k=0; k < attached.n13[iatom]; k++)
275     skip[attached.i13[k][iatom]] = iatom;
276     for(k=0; k < attached.n14[iatom]; k++)
277     skip[attached.i14[k][iatom]] = -iatom;
278    
279     ia = iatom;
280     ita = atom[ia].type;
281     xi = atom[ia].x;
282     yi = atom[ia].y;
283     zi = atom[ia].z;
284    
285     for (j=1; j <= natom; j++)
286     {
287     ib = j;
288     if (skip[j] != iatom)
289     {
290     itb = atom[ib].type;
291     if (ita <= itb)
292     it = 1000*ita + itb;
293     else
294     it = 1000*itb + ita;
295    
296     xk = atom[ib].x;
297     yk = atom[ib].y;
298     zk = atom[ib].z;
299    
300     xr = xi - xk;
301     yr = yi - yk;
302     zr = zi - zk;
303     rik2 = xr*xr + yr*yr + zr*zr;
304     if (rik2 < cutoff)
305     {
306     rv = nonbond.vrad[ita][itb];
307     eps = nonbond.veps[ita][itb];
308     if (skip[j] == -iatom)
309     eps /= units.v14scale;
310    
311     rik = sqrt(rik2);
312     rv7 = pow(rv,7.0);
313     rv14 = rv7*rv7;
314     rik6 = rik2*rik2*rik2;
315     rik5 = rik6/rik;
316     rik7 = rik6*rik;
317     rik12 = rik6*rik6;
318     rho = rik7 + (sigma-1.00)*rv7;
319     tau = rik + (kappa-1.00)*rv;
320    
321     de = -7.00*eps*kappa7*(rv7/pow(tau,8.00))*(sigma*rv7/rho-2.00)
322     -7.00*eps*kappa7*sigma*rv14*rik6/(rho*rho*pow(tau,7.00));
323    
324     d2e = 56.0*eps*kappa7*(rv7/pow(tau,9.00))*(sigma*rv7/rho-2.00)
325     + 98.0*eps*kappa7*sigma*rv14*rik6/(rho*rho*pow(tau,8.00))
326     + 98.0*eps*kappa7*sigma*rv14*rik12/(rho*rho*rho*pow(tau,7.00))
327     - 42.0*eps*kappa7*sigma*rv14*rik5/(rho*rho*pow(tau,7.00));
328    
329     de /= rik;
330     d2e = (d2e-de)/rik2;
331     d2edx = d2e*xr;
332     d2edy = d2e*yr;
333     d2edz = d2e*zr;
334    
335     term[0][0] = d2edx*xr + de;
336     term[1][0] = d2edx*yr;
337     term[2][0] = d2edx*zr;
338     term[0][1] = term[1][0];
339     term[1][1] = d2edy*yr + de;
340     term[2][1] = d2edy*zr;
341     term[0][2] = term[2][0];
342     term[1][2] = term[2][1];
343     term[2][2] = d2edz*zr + de;
344    
345     if (ia == iatom)
346     {
347     for (k=0; k < 3; k++)
348     {
349     hess.hessx[ia][k] += term[k][0];
350     hess.hessy[ia][k] += term[k][1];
351     hess.hessz[ia][k] += term[k][2];
352     hess.hessx[ib][k] -= term[k][0];
353     hess.hessy[ib][k] -= term[k][1];
354     hess.hessz[ib][k] -= term[k][2];
355     }
356     } else
357     {
358     for (k=0; k < 3; k++)
359     {
360     hess.hessx[ia][k] -= term[k][0];
361     hess.hessy[ia][k] -= term[k][1];
362     hess.hessz[ia][k] -= term[k][2];
363     hess.hessx[ib][k] += term[k][0];
364     hess.hessy[ib][k] += term[k][1];
365     hess.hessz[ib][k] += term[k][2];
366     }
367     }
368     }
369     }
370     }
371     }