ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/freemol/branches/sync4pymol12/src/mengine/src/ksolv.c
Revision: 99
Committed: Mon Jan 19 04:31:29 2009 UTC (10 years, 9 months ago) by wdelano
File size: 11040 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 "attached.h"
9    
10     #define STILL 1
11     #define HCT 2
12    
13     EXTERN struct t_solvent {
14     int type;
15     double EPSin, EPSsolv;
16     double doffset, p1,p2,p3,p4,p5;
17     double *shct,*asolv,*rsolv,*vsolv,*gpol,*rborn;
18     } solvent;
19    
20    
21     int find_bond(int,int);
22     void born(void);
23     void ksolv(void);
24     void message_alert(char *, char *);
25    
26     // choices are gbsa models: Analytical Still and HCT
27     // allocate memory in get_mem
28    
29     void ksolv()
30     {
31     int i,j,k, jji,it,kt;
32     int ia,ib,ic,id;
33     double ri,ri2,rk,cc;
34     double r,r2,r4,rab,rbc;
35     double cosine,factor;
36     double h,ratio,term;
37    
38     solvent.doffset = -0.09;
39     solvent.p1 = 0.073;
40     solvent.p2 = 0.921;
41     solvent.p3 = 6.211;
42     solvent.p4 = 15.236;
43     solvent.p5 = 1.254;
44    
45     if (solvent.type == HCT)
46     {
47     for (i=1; i <= natom; i++)
48     {
49     solvent.asolv[i] = 0.0054;
50     solvent.shct[i] = 0.0;
51     if (atom[i].atomnum == 1) solvent.shct[i] = 0.85;
52     if (atom[i].mmx_type == 20) solvent.shct[i] = 0.85;
53     if (atom[i].atomnum == 6) solvent.shct[i] = 0.72;
54     if (atom[i].atomnum == 7) solvent.shct[i] = 0.79;
55     if (atom[i].atomnum == 8) solvent.shct[i] = 0.85;
56     if (atom[i].atomnum == 9) solvent.shct[i] = 0.88;
57     if (atom[i].atomnum == 15) solvent.shct[i] = 0.86;
58     if (atom[i].atomnum == 16) solvent.shct[i] = 0.96;
59     if (atom[i].atomnum == 26) solvent.shct[i] = 0.88;
60     }
61     } else if (solvent.type == STILL)
62     {
63     for (i=1; i <= natom; i++)
64     {
65     solvent.asolv[i] = 0.0049;
66     }
67     }
68     // assign standard radii
69     for (i=1; i <= natom; i++)
70     {
71     solvent.rsolv[i] = 0.0;
72     if (atom[i].atomnum == 1)
73     {
74     solvent.rsolv[i] = 1.25;
75     if (atom[atom[i].iat[0]].atomnum == 7) solvent.rsolv[i] = 1.15;
76     if (atom[atom[i].iat[0]].atomnum == 8) solvent.rsolv[i] = 1.05;
77     } else if (atom[i].mmx_type == 20)
78     {
79     solvent.rsolv[i] = 0.8;
80     } else if (atom[i].atomnum == 3)
81     {
82     solvent.rsolv[i] = 1.432;
83     } else if (atom[i].atomnum == 6)
84     {
85     solvent.rsolv[i] = 1.90;
86     jji = 0;
87     for (j=0; j < 4; j++)
88     {
89     if (atom[i].iat[j] != 0)
90     jji++;
91     }
92     if (jji == 3) solvent.rsolv[i] = 1.875;
93     if (jji == 2) solvent.rsolv[i] = 1.825;
94     } else if (atom[i].atomnum == 7)
95     {
96     solvent.rsolv[i] = 1.7063;
97     jji = 0;
98     for (j=0; j < 4; j++)
99     {
100     if (atom[i].iat[j] != 0)
101     jji++;
102     }
103     if (jji == 4) solvent.rsolv[i] = 1.625;
104     if (jji == 1) solvent.rsolv[i] = 1.60;
105     } else if (atom[i].atomnum == 8)
106     {
107     solvent.rsolv[i] = 1.535;
108     jji = 0;
109     for (j=0; j < 4; j++)
110     {
111     if (atom[i].iat[j] != 0)
112     jji++;
113     }
114     if (jji == 1) solvent.rsolv[i] = 1.48;
115     } else if (atom[i].atomnum == 9)
116     {
117     solvent.rsolv[i] = 1.47;
118     } else if (atom[i].atomnum == 10)
119     {
120     solvent.rsolv[i] = 1.39;
121     } else if (atom[i].atomnum == 11)
122     {
123     solvent.rsolv[i] = 1.992;
124     } else if (atom[i].atomnum == 12)
125     {
126     solvent.rsolv[i] = 1.70;
127     } else if (atom[i].atomnum == 14)
128     {
129     solvent.rsolv[i] = 1.80;
130     } else if (atom[i].atomnum == 15)
131     {
132     solvent.rsolv[i] = 1.87;
133     } else if (atom[i].atomnum == 16)
134     {
135     solvent.rsolv[i] = 1.775;
136     } else if (atom[i].atomnum == 17)
137     {
138     solvent.rsolv[i] = 1.735;
139     } else if (atom[i].atomnum == 18)
140     {
141     solvent.rsolv[i] = 1.70;
142     } else if (atom[i].atomnum == 19)
143     {
144     solvent.rsolv[i] = 2.123;
145     } else if (atom[i].atomnum == 20)
146     {
147     solvent.rsolv[i] = 1.817;
148     } else if (atom[i].atomnum == 35)
149     {
150     solvent.rsolv[i] = 1.90;
151     } else if (atom[i].atomnum == 36)
152     {
153     solvent.rsolv[i] = 1.812;
154     } else if (atom[i].atomnum == 37)
155     {
156     solvent.rsolv[i] = 2.26;
157     } else if (atom[i].atomnum == 53)
158     {
159     solvent.rsolv[i] = 2.10;
160     } else if (atom[i].atomnum == 54)
161     {
162     solvent.rsolv[i] = 1.967;
163     } else if (atom[i].atomnum == 55)
164     {
165     solvent.rsolv[i] = 2.507;
166     } else if (atom[i].atomnum == 56)
167     {
168     solvent.rsolv[i] = 2.188;
169     }else
170     {
171     message_alert("Error - Atom not available for GB/SA","Error");
172     break;
173     }
174     }
175     // atomic volumes for Still
176     if (solvent.type == STILL)
177     {
178     for (i=1; i <= natom; i++)
179     {
180     if (atom[i].mmx_type != 20)
181     {
182     solvent.vsolv[i] = (4.0*PI/3.0)*solvent.rsolv[i]*solvent.rsolv[i]*solvent.rsolv[i];
183     ri = solvent.rsolv[i];
184     ri2 = ri*ri;
185     for (j=0; j < MAXIAT; j++)
186     {
187     if (atom[i].iat[j] != 0 && atom[i].bo[j] != 9 && atom[atom[i].iat[j]].mmx_type != 20)
188     {
189     k = find_bond(i,atom[i].iat[j]);
190     rk = solvent.rsolv[atom[i].iat[j]];
191     r = 1.01 * bonds_ff.bl[k];
192     ratio = (rk*rk - ri2 - r*r)/(2.0*ri*r);
193     h = ri*(1.0+ratio);
194     term = (PI/3.0)*h*h*(3.0*ri-h);
195     solvent.vsolv[i] -= term;
196     }
197     }
198     }
199     }
200    
201     // 1,2 and 1,3 polarization
202     cc = 4.80298*4.80298*14.39418;
203    
204     for (i=1; i <= natom; i++)
205     {
206     if (atom[i].mmx_type != 20)
207     solvent.gpol[i] = -0.5*cc/(solvent.rsolv[i]+solvent.doffset+solvent.p1);
208     }
209    
210     for (i=0; i < bonds_ff.nbnd; i++)
211     {
212     it = bonds_ff.i12[i][0];
213     kt = bonds_ff.i12[i][1];
214     if (atom[it].mmx_type != 20 && atom[kt].mmx_type != 20)
215     {
216     r = bonds_ff.bl[i];
217     r4 = r*r*r*r;
218     solvent.gpol[it] += solvent.p2*solvent.vsolv[kt]/r4;
219     solvent.gpol[kt] += solvent.p2*solvent.vsolv[it]/r4;
220     }
221     }
222    
223     for (i=0; i < angles.nang; i++)
224     {
225     ia = angles.i13[i][0];
226     ib = angles.i13[i][1];
227     ic = angles.i13[i][2];
228     if (atom[ia].mmx_type != 20 && atom[ib].mmx_type != 20 && atom[ic].mmx_type != 20)
229     {
230     factor = 1.0;
231     for (j=0; j < MAXIAT; j++)
232     {
233     if (atom[ia].iat[j] != 0)
234     {
235     id = atom[ia].iat[j];
236     if (id == ic)
237     factor = 0.0;
238     else if (id != ib)
239     {
240     for (k=0; k < MAXIAT; k++)
241     {
242     if (atom[ic].iat[k] != 0)
243     {
244     if (atom[ic].iat[k] == id)
245     factor = 0.5;
246     }
247     }
248     }
249     }
250     }
251    
252     id = find_bond(ia,ib);
253     rab = bonds_ff.bl[id];
254     id = find_bond(ib,ic);
255     rbc = bonds_ff.bl[id];
256     cosine = cos(angles.anat[i]/radian);
257     r2 = rab*rab + rbc*rbc - 2.0*rab*rbc*cosine;
258     r4 = r2*r2;
259     solvent.gpol[ia] += factor*solvent.p3*solvent.vsolv[ic]/r4;
260     solvent.gpol[ic] += factor*solvent.p3*solvent.vsolv[ia]/r4;
261     }
262     }
263     }
264     }
265     /* =================================================== */
266     void born()
267     {
268     int i,k;
269     double cc;
270     double ratio;
271     double xi,yi,zi,ri;
272     double rk,sk,sk2,sum;
273     double lik,lik2,uik,uik2;
274     double xr,yr,zr,rvdw;
275     double r,r2,r4;
276     double gpi,pip5,p5inv;
277     double theta,term,ccf;
278    
279     cc = 4.80298*4.80298*14.39418;
280     if (solvent.type == STILL)
281     {
282     p5inv = 1.0/solvent.p5;
283     pip5 = PI*solvent.p5;
284     for (i=1; i <= natom; i++)
285     {
286     if (atom[i].mmx_type != 20)
287     {
288     xi = atom[i].x;
289     yi = atom[i].y;
290     zi = atom[i].z;
291     gpi = solvent.gpol[i];
292 wdelano 99 skip[i][i] = i;
293 wdelano 58 for (k = 1; k <= natom; k++)
294     {
295 wdelano 99 if (skip[i][k] != i && atom[k].mmx_type != 20)
296 wdelano 58 {
297     xr = atom[k].x - xi;
298     yr = atom[k].y - yi;
299     zr = atom[k].z - zi;
300     r2 = xr*xr + yr*yr + zr*zr;
301     r4 = r2*r2;
302     rvdw = solvent.rsolv[i] + solvent.rsolv[k];
303     ratio = r2/(rvdw*rvdw);
304     if (ratio > p5inv)
305     {
306     ccf = 1.0;
307     } else
308     {
309     theta = ratio*pip5;
310     term = 0.5*(1.0-cos(theta));
311     ccf = term*term;
312     }
313     gpi += solvent.p4*ccf*solvent.vsolv[k]/r4;
314     }
315     }
316     solvent.rborn[i] = -0.5*cc/gpi;
317     }
318     }
319     }else if (solvent.type == HCT)
320     {
321     for (i=1; i <= natom; i++)
322     {
323     if (atom[i].mmx_type != 20)
324     {
325     xi = atom[i].x;
326     yi = atom[i].y;
327     zi = atom[i].z;
328     ri = solvent.rsolv[i] + solvent.doffset;
329     sum = 1.0/ri;
330     for (k=1; k <= natom; k++)
331     {
332     if (i != k && atom[k].mmx_type != 20)
333     {
334     xr = atom[k].x - xi;
335     yr = atom[k].y - yi;
336     zr = atom[k].z - zi;
337     r2 = xr*xr + yr*yr + zr*zr;
338     r = sqrt(r2);
339     rk = solvent.rsolv[k]+solvent.doffset;
340     sk = rk * solvent.shct[k];
341     sk2 = sk*sk;
342     if (ri < (r+sk))
343     {
344     lik = 1.0/ MaxFun(ri,r-sk);
345     uik = 1.0/(r+sk);
346     lik2 = lik*lik;
347     uik2 = uik*uik;
348     term = lik - uik + 0.25*r*(uik2-lik2) + (0.5/r)*log(uik/lik)
349     + (0.25*sk2/r)*(lik2-uik2);
350     sum -= 0.5*term;
351     }
352     }
353     }
354     theta = 1.0/sum;
355     solvent.rborn[i] = MaxFun(ri,theta);
356     }
357     }
358     }
359     }
360    
361    
362    
363    
364