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