ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/freemol/trunk/src/mengine/src/esolv.c
Revision: 110
Committed: Thu Mar 12 01:43:50 2009 UTC (12 years, 1 month ago) by gilbertke
File size: 12641 byte(s)
Log Message:
further cleanup and localization of atom data
Line File contents
1 #define EXTERN extern
2
3 #include "pcwin.h"
4
5 #define STILL 1
6 #define HCT 2
7
8 EXTERN struct t_solvent {
9 int type;
10 double EPSin, EPSsolv;
11 double doffset, p1,p2,p3,p4,p5;
12 double *shct,*asolv,*rsolv,*vsolv,*gpol,*rborn;
13 } solvent;
14
15 void esolv(int natom,double chrgcut,int **skip,int *use,double *charge,double *x,double *y,double *z,double *esolv);
16 void esolv1(int natom,double chrgcut,int **skip,int *use,double *charge,double *x,double *y,double *z,double *esolv,double **desolv,double *drb);
17 void esolv2(int ia,int natom,double chrgcut,double *charge,double *x,double *y,double *z,float **hessx,float **hessy,float **hessz);
18 void born1(int natom,int **skip,double *x,double *y,double *z,double **desolv,double *drb);
19 void born(int natom,int **skip, double *x,double *y,double *z);
20 // ===============================
21 void esolv(int natom,double chrgcut,int **skip,int *use, double *charge,double *x,double *y,double *z,double *esolv)
22 {
23 int i,j;
24 double e,ai,ri,rb;
25 double term,probe,term1;
26 double dwater, cc,f,cutoff;
27 double xi,yi,zi,xr,yr,zr,rik2;
28 double rb2,fgb;
29
30 probe = 1.40;
31 dwater = 78.3;
32 cc = 4.80298*4.80298*14.39418;
33 cutoff = chrgcut*chrgcut;
34
35 born(natom,skip,x,y,z);
36 term = 4.0*PI;
37 for (i=1; i <= natom; i++)
38 {
39 ai = solvent.asolv[i];
40 ri = solvent.rsolv[i];
41 rb = solvent.rborn[i];
42 if (rb > 0.001)
43 {
44 term1 = (ri/rb)*(ri/rb);
45 e = ai*term*(ri+probe)*(ri+probe)*term1*term1*term1;
46 *esolv += e;
47 }
48
49 }
50 // polarization term
51 f = -cc*(1.0/solvent.EPSin - 1.0/solvent.EPSsolv);
52
53 for (i=1; i <= natom; i++)
54 {
55 if (charge[i] != 0.0)
56 {
57 for(j=i; j <= natom; j++)
58 {
59 if (use[i] || use[j])
60 {
61 if (charge[j] != 0.0)
62 {
63 xi = x[i];
64 yi = y[i];
65 zi = z[i];
66
67 xr = xi - x[j];
68 yr = yi - y[j];
69 zr = zi - z[j];
70 rik2 = xr*xr + yr*yr + zr*zr;
71 if (rik2 < cutoff)
72 {
73 rb2 = solvent.rborn[i]*solvent.rborn[j];
74 fgb = sqrt(rik2 + rb2*exp(-0.25*rik2/rb2));
75 e = f*charge[i]*charge[j]/fgb;
76 if (i == j) e *= 0.5;
77 *esolv += e;
78 }
79 }
80 }
81 }
82 }
83 }
84 }
85 // =======================
86 void esolv1(int natom,double chrgcut,int **skip,int *use,double *charge,double *x,double *y,double *z,double *esolv,double **desolv,double *drb)
87 {
88 int i,j;
89 double e,de,derb,ai,ri,rb;
90 double term,probe,term1;
91 double dwater, cc,f,cutoff;
92 double xi,yi,zi,xr,yr,zr,rik2,r;
93 double rb2,fgb,fgb2,expterm;
94 double rbi,rbk,dedx,dedy,dedz,drbi,drbk;
95
96 probe = 1.40;
97 dwater = 78.3;
98 cc = 4.80298*4.80298*14.39418;
99 cutoff = chrgcut*chrgcut;
100 for (i=1; i <= natom; i++)
101 {
102 drb[i] = 0.0;
103 desolv[i][0] = 0.0;
104 desolv[i][1] = 0.0;
105 desolv[i][2] = 0.0;
106 }
107
108 born(natom,skip,x,y,z);
109 term = 4.0*PI;
110 for (i=1; i <= natom; i++)
111 {
112 ai = solvent.asolv[i];
113 ri = solvent.rsolv[i];
114 rb = solvent.rborn[i];
115 if (rb > 0.01)
116 {
117 term1 = (ri/rb)*(ri/rb);
118 e = ai*term*(ri+probe)*(ri+probe)*term1*term1*term1;
119 *esolv += e;
120 drb[i] -= 6.0*e/rb;
121 }
122 }
123 // polarization term
124 f = -cc*(1.0/solvent.EPSin - 1.0/solvent.EPSsolv);
125
126 for (i=1; i <= natom; i++)
127 {
128 if (fabs(charge[i]) > 0.001)
129 {
130 xi = x[i];
131 yi = y[i];
132 zi = z[i];
133 rbi = solvent.rborn[i];
134 for(j=i; j <= natom; j++)
135 {
136 if (use[i] || use[j])
137 {
138 if (fabs(charge[j]) > 0.001)
139 {
140 xr = xi - x[j];
141 yr = yi - y[j];
142 zr = zi - z[j];
143 rik2 = xr*xr + yr*yr + zr*zr;
144 if (rik2 < cutoff)
145 {
146 r = sqrt(rik2);
147 rbk = solvent.rborn[j];
148 rb2 = rbi*rbk;
149 expterm = exp(-0.25*rik2/rb2);
150 fgb2 = rik2 + rb2*expterm;
151 fgb = sqrt(fgb2);
152 e = f*charge[i]*charge[j]/fgb;
153 de = -e*(r-0.25*r*expterm)/fgb2;
154 derb = -e*expterm*(0.5+0.125*rik2/rb2)/fgb2;
155 if (i == j)
156 {
157 e *= 0.5;
158 *esolv += e;
159 drbi = derb*rbk;
160 drb[i] += drbi;
161 } else
162 {
163 *esolv += e;
164 de /= r;
165 dedx = de*xr;
166 dedy = de*yr;
167 dedz = de*zr;
168 desolv[i][0] += dedx;
169 desolv[i][1] += dedy;
170 desolv[i][2] += dedz;
171 desolv[j][0] -= dedx;
172 desolv[j][1] -= dedy;
173 desolv[j][2] -= dedz;
174
175 drbi = derb*rbk;
176 drbk = derb*rbi;
177 drb[i] += drbi;
178 drb[j] += drbk;
179 }
180 }
181 }
182 }
183 }
184 }
185 }
186 //
187 born1(natom,skip,x,y,z,desolv,drb);
188 }
189 /* ===================================================== */
190 void born1(int natom,int **skip,double *x,double *y,double *z,double **desolv,double *drb)
191 {
192 int i,k;
193 double xi,yi,zi,cc,rvdw;
194 double xr,yr,zr;
195 double de;
196 double r,r2,r6;
197 double p5inv,pip5;
198 double gpi,vk,ratio;
199 double ccf,cosq,dccf;
200 double sinq,term,theta;
201 double rb2,ri,rk,sk,sk2;
202 double lik,lik2,lik3;
203 double uik,uik2,uik3;
204 double dlik,duik;
205 double t1,t2,t3;
206 double dedx,dedy,dedz;
207
208 cc = 4.80298*4.80298*14.39418;
209 if (solvent.type == STILL)
210 {
211 p5inv = 1.0/solvent.p5;
212 pip5 = PI*solvent.p5;
213 for (i=1; i <= natom; i++)
214 {
215 skip[i][i] = i;
216 xi = x[i];
217 yi = y[i];
218 zi = z[i];
219 gpi = 2.0*solvent.rborn[i]*solvent.rborn[i]/cc;
220 for (k = 1; k <= natom; k++)
221 {
222 if (skip[i][k] != i)
223 {
224 xr = x[k] - xi;
225 yr = y[k] - yi;
226 zr = z[k] - zi;
227 vk = solvent.vsolv[k];
228 r2 = xr*xr + yr*yr + zr*zr;
229 r = sqrt(r2);
230 r6 = r2*r2*r2;
231 rvdw = solvent.rsolv[i] + solvent.rsolv[k];
232 ratio = r2/(rvdw*rvdw);
233 if (ratio > p5inv)
234 {
235 ccf = 1.0;
236 dccf = 0.0;
237 } else
238 {
239 theta = ratio*pip5;
240 cosq = cos(theta);
241 term = 0.5*(1.0-cosq);
242 ccf = term*term;
243 sinq = sin(theta);
244 dccf = 2.0*term*sinq*pip5*ratio;
245 }
246 de = drb[i]*solvent.p4*gpi*vk*(4.0*ccf-dccf)/r6;
247 dedx = de*xr;
248 dedy = de*yr;
249 dedz = de*zr;
250 desolv[i][0] += dedx;
251 desolv[i][1] += dedy;
252 desolv[i][2] += dedz;
253 desolv[k][0] -= dedx;
254 desolv[k][1] -= dedy;
255 desolv[k][2] -= dedz;
256
257 }
258 }
259
260 }
261 } else if (solvent.type == HCT)
262 {
263 for (i=1; i <= natom; i++)
264 {
265 xi = x[i];
266 yi = y[i];
267 zi = z[i];
268 ri = solvent.rsolv[i] + solvent.doffset;
269 rb2 = solvent.rborn[i]*solvent.rborn[i];
270 for (k=1; k <= natom; k++)
271 {
272 if (i != k)
273 {
274 xr = x[k] - xi;
275 yr = y[k] - yi;
276 zr = z[k] - zi;
277 r2 = xr*xr + yr*yr + zr*zr;
278 r = sqrt(r2);
279 rk = solvent.rsolv[k]+solvent.doffset;
280 sk = rk * solvent.shct[k];
281 sk2 = sk*sk;
282 if (ri < (r+sk) )
283 {
284 lik = 1.0/ MaxFun(ri,r-sk);
285 uik = 1.0/(r+sk);
286 lik2 = lik*lik;
287 uik2 = uik*uik;
288 lik3 = lik * lik2;
289 uik3 = uik * uik2;
290 dlik = 1.0;
291 if (ri >= (r-sk)) dlik = 0.0;
292 duik = 1.0;
293 t1 = 0.5*lik2 + 0.25*sk2*lik3/r - 0.25*((lik/r)+(lik3*r));
294 t2 = -0.5*uik2 - 0.25*sk2*(uik3/r) + 0.25*((uik/r)+uik3*r);
295 t3 = 0.125*(1.0+sk2/r2)*(lik2-uik2) + 0.25*log(uik/lik)/r2;
296 de = drb[i] * rb2 * (dlik*t1+duik*t2+t3) / r;
297 dedx = de*xr;
298 dedy = de*yr;
299 dedz = de*zr;
300 desolv[i][0] += dedx;
301 desolv[i][1] += dedy;
302 desolv[i][2] += dedz;
303 desolv[k][0] -= dedx;
304 desolv[k][1] -= dedy;
305 desolv[k][2] -= dedz;
306 }
307 }
308 }
309
310 }
311 }
312 }
313 /* ============================================== */
314 void esolv2(int ia,int natom,double chrgcut,double *charge,double *x,double *y,double *z,float **hessx,float **hessy,float **hessz)
315 {
316 int jj,k;
317 double probe,cc,cutoff;
318 double fi,fik,de,d2e;
319 double d2edx,d2edy,d2edz;
320 double xi,yi,zi,xr,yr,zr;
321 double r,r2;
322 double dwater,rb2;
323 double expterm;
324 double fgb,fgb2,dfgb,dfgb2,d2fgb;
325 double term[3][3];
326
327 probe = 1.40;
328 dwater = 78.3;
329 cc = 4.80298*4.80298*14.39418;
330
331 if (fabs(charge[ia]) > 0.001)
332 return;
333
334 cutoff = chrgcut*chrgcut;
335 fi = -cc*(1.0/solvent.EPSin - 1.0/solvent.EPSsolv)*charge[ia];
336 xi = x[ia];
337 yi = y[ia];
338 zi = z[ia];
339
340 for (k = 1; k <= natom; k++)
341 {
342 if (ia != k && charge[k] != 0.0)
343 {
344 xr = xi - x[k];
345 yr = yi - y[k];
346 zr = zi - z[k];
347 r2 = xr*xr + yr*yr + zr*zr;
348 if (r2 < cutoff)
349 {
350 r = sqrt(r2);
351 fik = fi*charge[k];
352 rb2 = solvent.rborn[ia]*solvent.rborn[k];
353 expterm = exp(-0.25*r2/rb2);
354 fgb2 = r2 + rb2*expterm;
355 fgb = sqrt(fgb2);
356 dfgb = (1.0-0.25*expterm)*r/fgb;
357 dfgb2 = dfgb*dfgb;
358 d2fgb = -dfgb2/fgb + dfgb/r + 0.125*(r2/rb2)*expterm/fgb;
359 de = -fik*dfgb/fgb2;
360 d2e = -fik*(d2fgb-2.0*dfgb2/fgb)/fgb2;
361 de /= r;
362 d2e = (d2e-de)/r;
363 d2edx = d2e*xr;
364 d2edy = d2e*yr;
365 d2edz = d2e*zr;
366 term[0][0] = d2edx*xr + de;
367 term[1][0] = d2edx*yr;
368 term[2][0] = d2edx*zr;
369 term[0][1] = term[1][0];
370 term[1][1] = d2edy*yr + de;
371 term[2][1] = d2edy*zr;
372 term[0][2] = term[2][0];
373 term[1][2] = term[2][1];
374 term[2][2] = d2edz*zr + de;
375 for (jj=0; jj < 3; jj++)
376 {
377 hessx[ia][jj] += term[jj][0];
378 hessy[ia][jj] += term[jj][1];
379 hessz[ia][jj] += term[jj][2];
380 hessx[k][jj] -= term[jj][0];
381 hessy[k][jj] -= term[jj][1];
382 hessz[k][jj] -= term[jj][2];
383 }
384 }
385 }
386 }
387 }
388
389
390
391
392
393