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