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, 6 months ago) by wdelano
File size: 12658 byte(s)
Log Message:
synchronized with trunk, less openmp lib/include
Line File contents
1 #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 skip[i][i] = i;
231 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 if (skip[i][k] != i && atom[k].mmx_type != 20)
238 {
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