ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/freemol/trunk/src/mengine/src/ehal.c
Revision: 124
Committed: Wed Jun 17 01:42:14 2009 UTC (11 years, 1 month ago) by gilbertke
File size: 11161 byte(s)
Log Message:
updated code to remove compiler warnings
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;
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, 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 */
147 for (i=1; i < natom; i++)
148 {
149 ia = i;
150 ita = type[ia];
151 xi = x[ia];
152 yi = y[ia];
153 zi = z[ia];
154
155 for(j=i+1; j <= natom; j++)
156 {
157 if (use[i] || use[j])
158 {
159 if (skip[i][j] != i)
160 {
161 skipij = skip[i][j];
162 ib = j;
163 itb = type[j];
164 xk = x[ib];
165 yk = y[ib];
166 zk = z[ib];
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 = vrad[ita][itb];
175 eps = veps[ita][itb];
176 if (skipij == -i)
177 eps /= 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 if (skipij == -i)
197 {
198 sum_e14 += e;
199 lde14[ia][0] += dedx;
200 lde14[ia][1] += dedy;
201 lde14[ia][2] += dedz;
202 lde14[ib][0] -= dedx;
203 lde14[ib][1] -= dedy;
204 lde14[ib][2] -= dedz;
205 } else
206 {
207 sum_evdw += e;
208 ldevdw[ia][0] += dedx;
209 ldevdw[ia][1] += dedy;
210 ldevdw[ia][2] += dedz;
211 ldevdw[ib][0] -= dedx;
212 ldevdw[ib][1] -= dedy;
213 ldevdw[ib][2] -= dedz;
214 }
215 }
216 }
217 }
218 }
219 }
220 *evdw = sum_evdw;
221 *e14 = sum_e14;
222 for (i=1;i<= natom; i++)
223 {
224 devdw[i][0] = ldevdw[i][0];
225 devdw[i][1] = ldevdw[i][1];
226 devdw[i][2] = ldevdw[i][2];
227 de14[i][0] = lde14[i][0];
228 de14[i][1] = lde14[i][1];
229 de14[i][2] = lde14[i][2];
230 }
231 free_dmatrix(ldevdw, 0,natom+1,0,3);
232 free_dmatrix(lde14,0,natom+1,0,3);
233 }
234 // =============================
235 void ehal2(int iatom,int natom,int *type,int *use,double *x,double *y,double *z,double vdwcut,int **skip,double **vrad,double **veps,
236 float **hessx,float **hessy,float **hessz)
237 {
238 int j,ia, ib, k, ita, itb;
239 double xi,yi,zi,xr,yr,zr;
240 double xk,yk,zk;
241 double rv,eps;
242 double rik,rik2, rik5;
243 double rik6,rik7, rik12;
244 double sigma, kappa, kappa7;
245 double rv7,rho, tau, rv14;
246 double de,d2e;
247 double d2edx,d2edy,d2edz,term[3][3];
248 double cutoff;
249
250 cutoff = vdwcut*vdwcut;
251
252 sigma = 1.12;
253 kappa = 1.07;
254 kappa7 = pow(kappa,7.0);
255
256 ita = type[iatom];
257
258 ia = iatom;
259 ita = type[ia];
260 xi = x[ia];
261 yi = y[ia];
262 zi = z[ia];
263 skip[iatom][iatom] = iatom;
264 for (j=1; j <= natom; j++)
265 {
266 ib = j;
267 if (skip[iatom][j] != iatom)
268 {
269 itb = type[ib];
270
271 xk = x[ib];
272 yk = y[ib];
273 zk = z[ib];
274
275 xr = xi - xk;
276 yr = yi - yk;
277 zr = zi - zk;
278 rik2 = xr*xr + yr*yr + zr*zr;
279 if (rik2 < cutoff)
280 {
281 rv = vrad[ita][itb];
282 eps = veps[ita][itb];
283 if (skip[iatom][j] == -iatom)
284 eps /= units.v14scale;
285
286 rik = sqrt(rik2);
287 rv7 = pow(rv,7.0);
288 rv14 = rv7*rv7;
289 rik6 = rik2*rik2*rik2;
290 rik5 = rik6/rik;
291 rik7 = rik6*rik;
292 rik12 = rik6*rik6;
293 rho = rik7 + (sigma-1.00)*rv7;
294 tau = rik + (kappa-1.00)*rv;
295
296 de = -7.00*eps*kappa7*(rv7/pow(tau,8.00))*(sigma*rv7/rho-2.00)
297 -7.00*eps*kappa7*sigma*rv14*rik6/(rho*rho*pow(tau,7.00));
298
299 d2e = 56.0*eps*kappa7*(rv7/pow(tau,9.00))*(sigma*rv7/rho-2.00)
300 + 98.0*eps*kappa7*sigma*rv14*rik6/(rho*rho*pow(tau,8.00))
301 + 98.0*eps*kappa7*sigma*rv14*rik12/(rho*rho*rho*pow(tau,7.00))
302 - 42.0*eps*kappa7*sigma*rv14*rik5/(rho*rho*pow(tau,7.00));
303
304 de /= rik;
305 d2e = (d2e-de)/rik2;
306 d2edx = d2e*xr;
307 d2edy = d2e*yr;
308 d2edz = d2e*zr;
309
310 term[0][0] = d2edx*xr + de;
311 term[1][0] = d2edx*yr;
312 term[2][0] = d2edx*zr;
313 term[0][1] = term[1][0];
314 term[1][1] = d2edy*yr + de;
315 term[2][1] = d2edy*zr;
316 term[0][2] = term[2][0];
317 term[1][2] = term[2][1];
318 term[2][2] = d2edz*zr + de;
319
320 if (ia == iatom)
321 {
322 for (k=0; k < 3; k++)
323 {
324 hessx[ia][k] += term[k][0];
325 hessy[ia][k] += term[k][1];
326 hessz[ia][k] += term[k][2];
327 hessx[ib][k] -= term[k][0];
328 hessy[ib][k] -= term[k][1];
329 hessz[ib][k] -= term[k][2];
330 }
331 } else
332 {
333 for (k=0; k < 3; k++)
334 {
335 hessx[ia][k] -= term[k][0];
336 hessy[ia][k] -= term[k][1];
337 hessz[ia][k] -= term[k][2];
338 hessx[ib][k] += term[k][0];
339 hessy[ib][k] += term[k][1];
340 hessz[ib][k] += term[k][2];
341 }
342 }
343 }
344 }
345 }
346 }