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 User Rev File contents
1 tjod 3 #define EXTERN extern
2    
3 wdelano 106 #ifdef KEG_OPENMP
4 gilbertke 104 #include <omp.h>
5 wdelano 106 #endif
6 gilbertke 89
7 tjod 3 #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 gilbertke 103 //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 tjod 3
20 gilbertke 103 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 tjod 3 {
31 gilbertke 113 int i,ia, ib, ita, itb, j;
32 tjod 3 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 gilbertke 103 cutoff = vdwcut*vdwcut;
41     *evdw = 0.0;
42     *e14 = 0.0;
43 tjod 3 sigma = 1.12;
44     kappa = 1.07;
45     kappa7 = pow(kappa,7.0);
46    
47     if (minim_values.iprint)
48     {
49 wdelano 58 fprintf(pcmlogfile,"\nVDW Terms - MMFF Buffered 14-7 Potential\n");
50     fprintf(pcmlogfile," At1 At2 Rik Radius Eps Evdw\n");
51 tjod 3 }
52     for (i=1; i < natom; i++)
53     {
54    
55     for(j=i+1; j <= natom; j++)
56     {
57 gilbertke 89 if (skip[i][j] != i)
58 tjod 3 {
59     ia = i;
60     ib = j;
61 gilbertke 103 if (use[ia] || use[ib])
62 tjod 3 {
63 gilbertke 103 ita = type[ia];
64     itb = type[ib];
65 tjod 3
66 gilbertke 103 xi = x[ia];
67     yi = y[ia];
68     zi = z[ia];
69     xk = x[ib];
70     yk = y[ib];
71     zk = z[ib];
72 tjod 3
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 gilbertke 103 rv = vrad[ita][itb];
80     eps = veps[ita][itb];
81 gilbertke 89 if (skip[i][j] == -i)
82 tjod 3 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 gilbertke 89 if (skip[i][j] == -i)
91 gilbertke 103 *e14 += e;
92 tjod 3 else
93 gilbertke 103 *evdw += e;
94 tjod 3 if (minim_values.iprint)
95 wdelano 58 fprintf(pcmlogfile,"VDW: %-4d - %-4d %-8.3f %-8.3f %-8.3f = %8.4f\n",ia, ib, sqrt(rik2),rv,eps,e);
96 tjod 3 }
97     }
98     }
99     }
100     }
101     }
102 gilbertke 103 // ==============================
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 tjod 3 {
106 gilbertke 113 int i,ia, ib, j, ita, itb,skipij;
107 tjod 3 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 gilbertke 86 double rv7,rho, tau,tau7,tau8, rv14;
114 tjod 3 double dedx,dedy,dedz,de;
115 gilbertke 103 double cutoff,v14scale;
116     double sum_evdw,sum_e14;
117     double **ldevdw,**lde14;
118 tjod 3
119 gilbertke 103 icount++;
120     cutoff = vdwcut*vdwcut;
121     v14scale = units.v14scale;
122 tjod 3
123 gilbertke 103 *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 tjod 3
130     sigma = 1.12;
131     kappa = 1.07;
132 gilbertke 86 kappa7 = kappa*kappa*kappa*kappa*kappa*kappa*kappa;
133 tjod 3
134     for (i=1; i <= natom; i++)
135     {
136 gilbertke 103 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 tjod 3 }
143    
144 gilbertke 113 /*#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 tjod 3 for (i=1; i < natom; i++)
148     {
149     ia = i;
150 gilbertke 103 ita = type[ia];
151     xi = x[ia];
152     yi = y[ia];
153     zi = z[ia];
154 tjod 3
155     for(j=i+1; j <= natom; j++)
156     {
157 gilbertke 103 if (use[i] || use[j])
158 tjod 3 {
159 gilbertke 89 if (skip[i][j] != i)
160 tjod 3 {
161 gilbertke 103 skipij = skip[i][j];
162 tjod 3 ib = j;
163 gilbertke 103 itb = type[j];
164     xk = x[ib];
165     yk = y[ib];
166     zk = z[ib];
167 tjod 3
168     xr = xi - xk;
169     yr = yi - yk;
170     zr = zi - zk;
171     rik2 = xr*xr + yr*yr + zr*zr;
172 gilbertke 103 if (rik2 < cutoff)
173 tjod 3 {
174 gilbertke 103 rv = vrad[ita][itb];
175     eps = veps[ita][itb];
176     if (skipij == -i)
177     eps /= v14scale;
178 tjod 3 rik = sqrt(rik2);
179 gilbertke 87 rv7 = rv*rv*rv*rv*rv*rv*rv; // pow(rv,7.0);
180 tjod 3 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 gilbertke 87 tau7 = tau*tau*tau*tau*tau*tau*tau;
186     tau8 = tau7*tau;
187 gilbertke 86 e = eps*kappa7*(rv7/tau7)*(sigma*rv7/rho-2.00);
188 tjod 3
189 gilbertke 86 de = -7.00*eps*kappa7*(rv7/tau8)*(sigma*rv7/rho-2.00)
190     -7.00*eps*kappa7*sigma*rv14*rik6/(rho*rho*tau7);
191 tjod 3
192     de /= rik;
193     dedx = de*xr;
194     dedy = de*yr;
195     dedz = de*zr;
196 gilbertke 103 if (skipij == -i)
197 tjod 3 {
198 gilbertke 103 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 tjod 3 } else
206     {
207 gilbertke 103 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 tjod 3 }
215 gilbertke 103 }
216 tjod 3 }
217     }
218     }
219     }
220 gilbertke 103 *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 tjod 3 }
234 gilbertke 103 // =============================
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 tjod 3 {
238 gilbertke 124 int j,ia, ib, k, ita, itb;
239 tjod 3 double xi,yi,zi,xr,yr,zr;
240     double xk,yk,zk;
241 gilbertke 86 double rv,eps;
242 tjod 3 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 gilbertke 103 cutoff = vdwcut*vdwcut;
251 tjod 3
252     sigma = 1.12;
253     kappa = 1.07;
254     kappa7 = pow(kappa,7.0);
255    
256 gilbertke 103 ita = type[iatom];
257 tjod 3
258     ia = iatom;
259 gilbertke 103 ita = type[ia];
260     xi = x[ia];
261     yi = y[ia];
262     zi = z[ia];
263     skip[iatom][iatom] = iatom;
264 tjod 3 for (j=1; j <= natom; j++)
265     {
266     ib = j;
267 gilbertke 89 if (skip[iatom][j] != iatom)
268 tjod 3 {
269 gilbertke 103 itb = type[ib];
270 tjod 3
271 gilbertke 103 xk = x[ib];
272     yk = y[ib];
273     zk = z[ib];
274 tjod 3
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 gilbertke 103 rv = vrad[ita][itb];
282     eps = veps[ita][itb];
283 gilbertke 89 if (skip[iatom][j] == -iatom)
284 tjod 3 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 gilbertke 103 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 tjod 3 }
331     } else
332     {
333     for (k=0; k < 3; k++)
334     {
335 gilbertke 103 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 tjod 3 }
342     }
343     }
344     }
345     }
346     }