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