ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/freemol/trunk/src/mengine/src/ehal.c
Revision: 86
Committed: Mon Jan 12 16:21:45 2009 UTC (10 years, 10 months ago) by gilbertke
File size: 11116 byte(s)
Log Message:
rewrote ehal to remove pow(x,y) function call
Line File contents
1 #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 fprintf(pcmlogfile,"\nVDW Terms - MMFF Buffered 14-7 Potential\n");
51 fprintf(pcmlogfile," At1 At2 Rik Radius Eps Evdw\n");
52 }
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 fprintf(pcmlogfile,"VDW: %-4d - %-4d %-8.3f %-8.3f %-8.3f = %8.4f\n",ia, ib, sqrt(rik2),rv,eps,e);
112 }
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 double rv7,rho, tau,tau7,tau8, rv14;
129 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 kappa7 = kappa*kappa*kappa*kappa*kappa*kappa*kappa;
140
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 rv7 = rv7*rv7*rv7*rv7*rv7*rv7*rv7; // pow(rv,7.0);
198 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 tau7 = tau*tau*tau*tau*tau*tau*tau;
204 tau8 = tau7*tau;
205 e = eps*kappa7*(rv7/tau7)*(sigma*rv7/rho-2.00);
206
207 de = -7.00*eps*kappa7*(rv7/tau8)*(sigma*rv7/rho-2.00)
208 -7.00*eps*kappa7*sigma*rv14*rik6/(rho*rho*tau7);
209
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 int j,ia, ib, k, ita, itb, it;
247 double xi,yi,zi,xr,yr,zr;
248 double xk,yk,zk;
249 double rv,eps;
250 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 }