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