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, 3 months ago) by gilbertke
File size: 10342 byte(s)
Log Message:
further cleanup and localization of atom data
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
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 void born(int natom,int **skip, double *x,double *y,double *z);
21 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 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 }
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 if (atom.atomnum[i] == 1)
70 {
71 solvent.rsolv[i] = 1.25;
72 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 {
76 solvent.rsolv[i] = 1.432;
77 } else if (atom.atomnum[i] == 6)
78 {
79 solvent.rsolv[i] = 1.90;
80 jji = 0;
81 for (j=0; j < 4; j++)
82 {
83 if (atom.iat[i][j] != 0)
84 jji++;
85 }
86 if (jji == 3) solvent.rsolv[i] = 1.875;
87 if (jji == 2) solvent.rsolv[i] = 1.825;
88 } else if (atom.atomnum[i] == 7)
89 {
90 solvent.rsolv[i] = 1.7063;
91 jji = 0;
92 for (j=0; j < 4; j++)
93 {
94 if (atom.iat[i][j] != 0)
95 jji++;
96 }
97 if (jji == 4) solvent.rsolv[i] = 1.625;
98 if (jji == 1) solvent.rsolv[i] = 1.60;
99 } else if (atom.atomnum[i] == 8)
100 {
101 solvent.rsolv[i] = 1.535;
102 jji = 0;
103 for (j=0; j < 4; j++)
104 {
105 if (atom.iat[i][j] != 0)
106 jji++;
107 }
108 if (jji == 1) solvent.rsolv[i] = 1.48;
109 } else if (atom.atomnum[i] == 9)
110 {
111 solvent.rsolv[i] = 1.47;
112 } else if (atom.atomnum[i] == 10)
113 {
114 solvent.rsolv[i] = 1.39;
115 } else if (atom.atomnum[i] == 11)
116 {
117 solvent.rsolv[i] = 1.992;
118 } else if (atom.atomnum[i] == 12)
119 {
120 solvent.rsolv[i] = 1.70;
121 } else if (atom.atomnum[i] == 14)
122 {
123 solvent.rsolv[i] = 1.80;
124 } else if (atom.atomnum[i] == 15)
125 {
126 solvent.rsolv[i] = 1.87;
127 } else if (atom.atomnum[i] == 16)
128 {
129 solvent.rsolv[i] = 1.775;
130 } else if (atom.atomnum[i] == 17)
131 {
132 solvent.rsolv[i] = 1.735;
133 } else if (atom.atomnum[i] == 18)
134 {
135 solvent.rsolv[i] = 1.70;
136 } else if (atom.atomnum[i] == 19)
137 {
138 solvent.rsolv[i] = 2.123;
139 } else if (atom.atomnum[i] == 20)
140 {
141 solvent.rsolv[i] = 1.817;
142 } else if (atom.atomnum[i] == 35)
143 {
144 solvent.rsolv[i] = 1.90;
145 } else if (atom.atomnum[i] == 36)
146 {
147 solvent.rsolv[i] = 1.812;
148 } else if (atom.atomnum[i] == 37)
149 {
150 solvent.rsolv[i] = 2.26;
151 } else if (atom.atomnum[i] == 53)
152 {
153 solvent.rsolv[i] = 2.10;
154 } else if (atom.atomnum[i] == 54)
155 {
156 solvent.rsolv[i] = 1.967;
157 } else if (atom.atomnum[i] == 55)
158 {
159 solvent.rsolv[i] = 2.507;
160 } else if (atom.atomnum[i] == 56)
161 {
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 if (atom.iat[i][j] != 0 && atom.bo[i][j] != 9)
180 {
181 k = find_bond(i,atom.iat[i][j]);
182 rk = solvent.rsolv[atom.iat[i][j]];
183 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
191 }
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 solvent.gpol[i] = -0.5*cc/(solvent.rsolv[i]+solvent.doffset+solvent.p1);
198
199
200 for (i=0; i < bonds_ff.nbnd; i++)
201 {
202 it = bonds_ff.i12[i][0];
203 kt = bonds_ff.i12[i][1];
204 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 }
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 if (atom.iat[ia][j] != 0)
220 {
221 id = atom.iat[ia][j];
222 if (id == ic)
223 factor = 0.0;
224 else if (id != ib)
225 {
226 for (k=0; k < MAXIAT; k++)
227 {
228 if (atom.iat[ic][k] != 0)
229 {
230 if (atom.iat[ic][k] == id)
231 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
248 }
249 }
250 }
251 /* =================================================== */
252 void born(int natom,int **skip, double *x,double *y,double *z)
253 {
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 xi = x[i];
273 yi = y[i];
274 zi = z[i];
275 gpi = solvent.gpol[i];
276 skip[i][i] = i;
277 for (k = 1; k <= natom; k++)
278 {
279 if (skip[i][k] != i)
280 {
281 xr = x[k] - xi;
282 yr = y[k] - yi;
283 zr = z[k] - zi;
284 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
302 }
303 }else if (solvent.type == HCT)
304 {
305 for (i=1; i <= natom; i++)
306 {
307 xi = x[i];
308 yi = y[i];
309 zi = z[i];
310 ri = solvent.rsolv[i] + solvent.doffset;
311 sum = 1.0/ri;
312 for (k=1; k <= natom; k++)
313 {
314 if (i != k)
315 {
316 xr = x[k] - xi;
317 yr = y[k] - yi;
318 zr = z[k] - zi;
319 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
339 }
340 }
341 }
342
343
344
345
346