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, 6 months ago) by wdelano
File size: 11040 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 "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 skip[i][i] = i;
293 for (k = 1; k <= natom; k++)
294 {
295 if (skip[i][k] != i && atom[k].mmx_type != 20)
296 {
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