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