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, 7 months ago) by gilbertke
File size: 12641 byte(s)
Log Message:
further cleanup and localization of atom data
Line User Rev File contents
1 wdelano 58 #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 gilbertke 103 double EPSin, EPSsolv;
11 wdelano 58 double doffset, p1,p2,p3,p4,p5;
12     double *shct,*asolv,*rsolv,*vsolv,*gpol,*rborn;
13     } solvent;
14    
15 gilbertke 110 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 gilbertke 103 // ===============================
21 gilbertke 110 void esolv(int natom,double chrgcut,int **skip,int *use, double *charge,double *x,double *y,double *z,double *esolv)
22 wdelano 58 {
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 gilbertke 103 cutoff = chrgcut*chrgcut;
34 wdelano 58
35 gilbertke 110 born(natom,skip,x,y,z);
36 wdelano 58 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 gilbertke 103 if (rb > 0.001)
43 wdelano 58 {
44     term1 = (ri/rb)*(ri/rb);
45     e = ai*term*(ri+probe)*(ri+probe)*term1*term1*term1;
46 gilbertke 103 *esolv += e;
47 wdelano 58 }
48 gilbertke 110
49 wdelano 58 }
50     // polarization term
51 gilbertke 103 f = -cc*(1.0/solvent.EPSin - 1.0/solvent.EPSsolv);
52 wdelano 58
53 gilbertke 103 for (i=1; i <= natom; i++)
54 wdelano 58 {
55 gilbertke 110 if (charge[i] != 0.0)
56 wdelano 58 {
57     for(j=i; j <= natom; j++)
58     {
59 gilbertke 103 if (use[i] || use[j])
60 wdelano 58 {
61 gilbertke 110 if (charge[j] != 0.0)
62 wdelano 58 {
63 gilbertke 103 xi = x[i];
64     yi = y[i];
65     zi = z[i];
66 wdelano 58
67 gilbertke 103 xr = xi - x[j];
68     yr = yi - y[j];
69     zr = zi - z[j];
70 wdelano 58 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 gilbertke 103 e = f*charge[i]*charge[j]/fgb;
76 wdelano 58 if (i == j) e *= 0.5;
77 gilbertke 103 *esolv += e;
78 wdelano 58 }
79     }
80     }
81     }
82     }
83     }
84     }
85 gilbertke 103 // =======================
86 gilbertke 110 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 wdelano 58 {
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 gilbertke 103 cutoff = chrgcut*chrgcut;
100 wdelano 58 for (i=1; i <= natom; i++)
101     {
102 gilbertke 103 drb[i] = 0.0;
103     desolv[i][0] = 0.0;
104     desolv[i][1] = 0.0;
105     desolv[i][2] = 0.0;
106 wdelano 58 }
107    
108 gilbertke 110 born(natom,skip,x,y,z);
109 wdelano 58 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 gilbertke 103 if (rb > 0.01)
116 wdelano 58 {
117     term1 = (ri/rb)*(ri/rb);
118     e = ai*term*(ri+probe)*(ri+probe)*term1*term1*term1;
119 gilbertke 103 *esolv += e;
120     drb[i] -= 6.0*e/rb;
121 wdelano 58 }
122     }
123     // polarization term
124 gilbertke 103 f = -cc*(1.0/solvent.EPSin - 1.0/solvent.EPSsolv);
125 wdelano 58
126 gilbertke 103 for (i=1; i <= natom; i++)
127 wdelano 58 {
128 gilbertke 110 if (fabs(charge[i]) > 0.001)
129 gilbertke 103 {
130     xi = x[i];
131     yi = y[i];
132     zi = z[i];
133 wdelano 58 rbi = solvent.rborn[i];
134     for(j=i; j <= natom; j++)
135     {
136 gilbertke 103 if (use[i] || use[j])
137 wdelano 58 {
138 gilbertke 110 if (fabs(charge[j]) > 0.001)
139 gilbertke 103 {
140     xr = xi - x[j];
141     yr = yi - y[j];
142     zr = zi - z[j];
143 wdelano 58 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 gilbertke 103 e = f*charge[i]*charge[j]/fgb;
153 wdelano 58 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 gilbertke 103 *esolv += e;
159 wdelano 58 drbi = derb*rbk;
160 gilbertke 103 drb[i] += drbi;
161 wdelano 58 } else
162     {
163 gilbertke 103 *esolv += e;
164 wdelano 58 de /= r;
165     dedx = de*xr;
166     dedy = de*yr;
167     dedz = de*zr;
168 gilbertke 103 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 wdelano 58
175     drbi = derb*rbk;
176     drbk = derb*rbi;
177 gilbertke 103 drb[i] += drbi;
178     drb[j] += drbk;
179 wdelano 58 }
180     }
181     }
182     }
183     }
184     }
185     }
186     //
187 gilbertke 110 born1(natom,skip,x,y,z,desolv,drb);
188 wdelano 58 }
189     /* ===================================================== */
190 gilbertke 110 void born1(int natom,int **skip,double *x,double *y,double *z,double **desolv,double *drb)
191 wdelano 58 {
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 gilbertke 89 skip[i][i] = i;
216 gilbertke 103 xi = x[i];
217     yi = y[i];
218     zi = z[i];
219 wdelano 58 gpi = 2.0*solvent.rborn[i]*solvent.rborn[i]/cc;
220     for (k = 1; k <= natom; k++)
221     {
222 gilbertke 110 if (skip[i][k] != i)
223 wdelano 58 {
224 gilbertke 103 xr = x[k] - xi;
225     yr = y[k] - yi;
226     zr = z[k] - zi;
227 wdelano 58 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 gilbertke 103 de = drb[i]*solvent.p4*gpi*vk*(4.0*ccf-dccf)/r6;
247 wdelano 58 dedx = de*xr;
248     dedy = de*yr;
249     dedz = de*zr;
250 gilbertke 103 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 wdelano 58
257     }
258     }
259 gilbertke 110
260 wdelano 58 }
261     } else if (solvent.type == HCT)
262     {
263     for (i=1; i <= natom; i++)
264     {
265 gilbertke 103 xi = x[i];
266     yi = y[i];
267     zi = z[i];
268 wdelano 58 ri = solvent.rsolv[i] + solvent.doffset;
269     rb2 = solvent.rborn[i]*solvent.rborn[i];
270     for (k=1; k <= natom; k++)
271     {
272 gilbertke 110 if (i != k)
273 wdelano 58 {
274 gilbertke 103 xr = x[k] - xi;
275     yr = y[k] - yi;
276     zr = z[k] - zi;
277 wdelano 58 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 gilbertke 103 de = drb[i] * rb2 * (dlik*t1+duik*t2+t3) / r;
297 wdelano 58 dedx = de*xr;
298     dedy = de*yr;
299     dedz = de*zr;
300 gilbertke 103 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 wdelano 58 }
307     }
308     }
309 gilbertke 110
310 wdelano 58 }
311     }
312     }
313     /* ============================================== */
314 gilbertke 110 void esolv2(int ia,int natom,double chrgcut,double *charge,double *x,double *y,double *z,float **hessx,float **hessy,float **hessz)
315 wdelano 58 {
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 gilbertke 110 if (fabs(charge[ia]) > 0.001)
332 wdelano 58 return;
333    
334 gilbertke 103 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 wdelano 58
340     for (k = 1; k <= natom; k++)
341     {
342 gilbertke 110 if (ia != k && charge[k] != 0.0)
343 wdelano 58 {
344 gilbertke 103 xr = xi - x[k];
345     yr = yi - y[k];
346     zr = zi - z[k];
347 wdelano 58 r2 = xr*xr + yr*yr + zr*zr;
348     if (r2 < cutoff)
349     {
350     r = sqrt(r2);
351 gilbertke 103 fik = fi*charge[k];
352 wdelano 58 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 gilbertke 103 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 wdelano 58 }
384     }
385     }
386     }
387     }
388    
389    
390    
391    
392    
393