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