ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/freemol/branches/sync4pymol12/src/mengine/src/kvdw.c
Revision: 99
Committed: Mon Jan 19 04:31:29 2009 UTC (10 years, 5 months ago) by wdelano
File size: 12523 byte(s)
Log Message:
synchronized with trunk, less openmp lib/include
Line File contents
1 #define EXTERN extern
2 #include "pcwin.h"
3 #include "pcmod.h"
4 #include "nonbond.h"
5 #include "pot.h"
6 #include "field.h"
7 #include "atom_k.h"
8
9 void numeral(int,char *,int);
10
11 EXTERN struct t_vdw1 {
12 int nvdw;
13 float rad[MAXVDWCONST], eps[MAXVDWCONST];
14 int lpd[MAXVDWCONST], ihtyp[MAXVDWCONST], ihdon[MAXVDWCONST];
15 float alpha[MAXVDWCONST],n[MAXVDWCONST],a[MAXVDWCONST],g[MAXVDWCONST];
16 char da[MAXVDWCONST][2];
17 } vdw1;
18
19 EXTERN struct t_vdwpr_k {
20 int nvdwpr;
21 int ia1[MAXBONDCONST],ia2[MAXBONDCONST];
22 char kv[MAXBONDCONST][7];
23 float radius[MAXBONDCONST], eps[MAXBONDCONST];
24 } vdwpr_k;
25
26 struct t_mm3hbond {
27 int npair;
28 int htype[MAXBONDCONST];
29 int otype[MAXBONDCONST];
30 } mm3hbond;
31
32 EXTERN int Missing_constants;
33
34 void kvdw()
35 {
36 int i, j, k;
37 int it, ita, itb, found, ianum,ibnum;
38 int dapair;
39 char pa[4],pb[4],pt[7];
40 long int hbmask;
41 float rad_ia, rad_ib;
42 float eps1, eps2, seps1, seps2;
43 float rii,rjj,gammaij, rij, epsij;
44
45 // mm3 setup
46 if (field.type == MM3)
47 {
48 mm3hbond.npair = 0;
49 for (i=0; i < vdwpr_k.nvdwpr; i++)
50 {
51 if (atom_k.number[vdwpr_k.ia1[i]] == 1) // looking for hydrogens
52 {
53 mm3hbond.htype[mm3hbond.npair] = vdwpr_k.ia1[i];
54 mm3hbond.otype[mm3hbond.npair] = vdwpr_k.ia2[i];
55 mm3hbond.npair++;
56 } else if (atom_k.number[vdwpr_k.ia2[i]] == 1)
57 {
58 if (!( (vdwpr_k.ia2[i] == 5 && vdwpr_k.ia1[i] == 1) || (vdwpr_k.ia2[i] == 36 && vdwpr_k.ia1[i] == 1) ))
59 {
60 mm3hbond.htype[mm3hbond.npair] = vdwpr_k.ia2[i];
61 mm3hbond.otype[mm3hbond.npair] = vdwpr_k.ia1[i];
62 mm3hbond.npair++;
63 }
64 }
65 }
66 }
67 //
68 nonbond.npair = 0;
69 hbmask = 1 << HBOND_MASK;
70 for (i=0; i <= nonbond.maxnbtype; i++)
71 {
72 for (j=0; j <= nonbond.maxnbtype; j++)
73 nonbond.iNBtype[j][i] = 0;
74 }
75
76 for (i=1; i < natom; i++)
77 {
78 ita = atom[i].type;
79 ianum = atom[i].atomnum;
80 for (j= i+1; j <= natom; j++)
81 {
82 itb = atom[j].type;
83 ibnum = atom[j].atomnum;
84 if (ita < itb)
85 it = 1000*ita + itb;
86 else
87 it = 1000*itb + ita;
88
89 found = FALSE;
90 if (nonbond.iNBtype[ita][itb] != 0)
91 found = TRUE;
92
93 if (found == FALSE)
94 {
95 if (field.type == MMFF94)
96 {
97 dapair = FALSE;
98 rii = vdw1.a[ita]*pow(vdw1.alpha[ita],0.25);
99 rjj = vdw1.a[itb]*pow(vdw1.alpha[itb],0.25);
100 if ((atom[i].atomnum == 1 && ita != 5) || (atom[j].atomnum ==1 && itb != 5) ) // polar hydrogens
101 {
102 rij = 0.5*(rii+rjj);
103 if ( (strcmp(vdw1.da[ita],"A") == 0) && (strcmp(vdw1.da[itb],"D") == 0) )
104 {
105 dapair = TRUE;
106 }
107 if ( (strcmp(vdw1.da[itb],"A") == 0) && (strcmp(vdw1.da[ita],"D") == 0) )
108 {
109 dapair = TRUE;
110 }
111 } else
112 {
113 gammaij = (rii-rjj)/(rii+rjj);
114 rij = 0.5*(rii+rjj)*(1.00+0.2*(1.00- exp(-12.0*gammaij*gammaij)));
115 }
116 epsij = ((181.16*vdw1.g[ita]*vdw1.g[itb]*vdw1.alpha[ita]*vdw1.alpha[itb])/
117 ( sqrt(vdw1.alpha[ita]/vdw1.n[ita]) + sqrt(vdw1.alpha[itb]/vdw1.n[itb])) ) / pow(rij,6.0);
118
119 if (dapair == TRUE)
120 {
121 rij *= 0.8;
122 epsij *= 0.5;
123 }
124
125 nonbond.vrad[ita][itb] = rij;
126 nonbond.veps[ita][itb] = epsij;
127 nonbond.ipif[ita][itb] = 1.0;
128 nonbond.vrad[itb][ita] = rij;
129 nonbond.veps[itb][ita] = epsij;
130 nonbond.ipif[itb][ita] = 1.0;
131
132 nonbond.vrad14[ita][itb] = nonbond.vrad[ita][itb];
133 nonbond.vrad14[itb][ita] = nonbond.vrad[itb][ita];
134 nonbond.veps14[ita][itb] = nonbond.veps[ita][itb];
135 nonbond.veps14[itb][ita] = nonbond.veps[itb][ita];
136
137 } else
138 {
139 rad_ia = vdw1.rad[ita];
140 rad_ib = vdw1.rad[itb];
141 if (strcmp(field.radiustype,"SIGMA") == 0)
142 {
143 rad_ia *= 1.122462048;
144 rad_ib *= 1.122462048;
145 }
146 if (strcmp(field.radiussize,"DIAMETER") == 0)
147 {
148 rad_ia /= 2.0;
149 rad_ib /= 2.0;
150 }
151 eps1 = vdw1.eps[ita];
152 eps2 = vdw1.eps[itb];
153
154 if (strcmp(field.radiusrule,"ARITHMETIC") == 0)
155 {
156 nonbond.vrad[ita][itb] = rad_ia + rad_ib;
157 nonbond.vrad[itb][ita] = rad_ia + rad_ib;
158 } else if (strcmp(field.radiusrule,"GEOMETRIC") == 0)
159 {
160 nonbond.vrad[ita][itb] = 2.0*( sqrt(rad_ia)*sqrt(rad_ib));
161 nonbond.vrad[itb][ita] = 2.0*( sqrt(rad_ia)*sqrt(rad_ib));
162 } else if (strcmp(field.radiusrule,"CUBIC-MEAN") == 0)
163 {
164 if ((atom[i].atomnum == 1 && ita != 5) || (atom[j].atomnum ==1 && itb != 5) ) // polar hydrogens
165 {
166 nonbond.vrad[ita][itb] = rad_ia + rad_ib;
167 nonbond.vrad[itb][ita] = rad_ia + rad_ib;
168 }else
169 {
170 nonbond.vrad[ita][itb] = 2.0*( ((rad_ia*rad_ia*rad_ia) + (rad_ib*rad_ib*rad_ib))/(rad_ia*rad_ia + rad_ib*rad_ib) );
171 nonbond.vrad[itb][ita] = 2.0*( ((rad_ia*rad_ia*rad_ia) + (rad_ib*rad_ib*rad_ib))/(rad_ia*rad_ia + rad_ib*rad_ib) );
172 }
173 } else
174 {
175 nonbond.vrad[ita][itb] = rad_ia + rad_ib;
176 nonbond.vrad[itb][ita] = rad_ia + rad_ib;
177 }
178
179 if (strcmp(field.epsrule,"ARITHMETIC") == 0)
180 {
181 nonbond.veps[ita][itb] = 0.5*(eps1 + eps2);
182 nonbond.veps[itb][ita] = 0.5*(eps1 + eps2);
183 } else if (strcmp(field.epsrule,"GEOMETRIC") == 0)
184 {
185 nonbond.veps[ita][itb] = sqrt( eps1*eps2 );
186 nonbond.veps[itb][ita] = sqrt( eps1*eps2 );
187 } else if (strcmp(field.epsrule,"HARMONIC") == 0)
188 {
189 nonbond.veps[ita][itb] = 2.0* (eps1*eps2)/(eps1+eps2);
190 nonbond.veps[itb][ita] = 2.0* (eps1*eps2)/(eps1+eps2);
191 } else if (strcmp(field.epsrule,"HHG") == 0)
192 {
193 if ((atom[i].atomnum == 1 && ita != 5) || (atom[j].atomnum ==1 && itb != 5) ) // polar hydrogens
194 {
195 nonbond.veps[ita][itb] = 0.5*(eps1 + eps2);
196 nonbond.veps[itb][ita] = 0.5*(eps1 + eps2);
197 } else
198 {
199 seps1 = sqrt(eps1);
200 seps2 = sqrt(eps2);
201 nonbond.veps[ita][itb] = 4.0*(eps1*eps2)/( (seps1+seps2)*(seps1+seps2));
202 nonbond.veps[itb][ita] = 4.0*(eps1*eps2)/( (seps1+seps2)*(seps1+seps2));
203 }
204 }else
205 {
206 nonbond.veps[ita][itb] = sqrt( eps1*eps2 );
207 nonbond.veps[itb][ita] = sqrt( eps1*eps2 );
208 }
209
210 if (field.type == MM3)
211 {
212 nonbond.ipif[ita][itb] = 1;
213 nonbond.ipif[itb][ita] = 1;
214 } else if (field.type == MMX)
215 {
216 nonbond.ipif[ita][itb] = 1;
217 nonbond.ipif[itb][ita] = 1;
218 if( (ita == 2 || ita == 4 || ita == 40 || ita == 48) &&
219 (itb == 2 || itb == 4 || itb == 40 || itb == 48) )
220 {
221 nonbond.ipif[ita][itb] = 0;
222 nonbond.ipif[itb][ita] = 0;
223 }
224 }
225 nonbond.vrad14[ita][itb] = nonbond.vrad[ita][itb];
226 nonbond.vrad14[itb][ita] = nonbond.vrad[itb][ita];
227 nonbond.veps14[ita][itb] = nonbond.veps[ita][itb];
228 nonbond.veps14[itb][ita] = nonbond.veps[itb][ita];
229 }
230 // look for any special vdw pairs
231 numeral(atom[i].type,pa,3);
232 numeral(atom[j].type,pb,3);
233
234 if (atom[i].type < atom[j].type)
235 {
236 strcpy(pt,pa);
237 strcat(pt,pb);
238 } else
239 {
240 strcpy(pt,pb);
241 strcat(pt,pa);
242 }
243
244 // look for vdw pairs
245 for (k= 0; k < vdwpr_k.nvdwpr ; k++)
246 {
247 if (strcmp(vdwpr_k.kv[k],pt) == 0)
248 {
249 nonbond.vrad[ita][itb] = vdwpr_k.radius[k];
250 nonbond.vrad[itb][ita] = vdwpr_k.radius[k];
251 nonbond.veps[ita][itb] = vdwpr_k.eps[k];
252 nonbond.veps[itb][ita] = vdwpr_k.eps[k];
253 if (ianum != 1 && ibnum != 1 )
254 {
255 nonbond.veps[ita][itb] /= units.dielec;
256 nonbond.veps[itb][ita] /= units.dielec;
257 }
258 if (ita == 1 && itb == 5)
259 {
260 nonbond.vrad14[ita][itb] = nonbond.vrad[ita][itb];
261 nonbond.vrad14[itb][ita] = nonbond.vrad[itb][ita];
262 nonbond.veps14[itb][ita] = nonbond.veps[itb][ita];
263 nonbond.veps14[ita][itb] = nonbond.veps[ita][itb];
264 }
265 break;
266 }
267 }
268
269 nonbond.iNBtype[ita][itb] = it;
270 nonbond.iNBtype[itb][ita] = it;
271
272 if (pot.use_hbond)
273 {
274 if ( atom[i].flags & hbmask || atom[i].mmx_type == 20)
275 {
276 for (k=0; k < MAXIAT; k++)
277 {
278 if (atom[j].iat[k] != 0 && atom[j].bo[k] != 9)
279 {
280 if (atom[atom[j].iat[k]].flags & hbmask || atom[atom[j].iat[k]].mmx_type == 20)
281 {
282 nonbond.vrad[ita][itb] -= 0.25;
283 nonbond.vrad[itb][ita] -= 0.25;
284 break;
285 }
286 }
287 }
288 } else if (atom[j].flags & hbmask || atom[j].mmx_type == 20)
289 {
290 for (k=0; k < MAXIAT; k++)
291 {
292 if (atom[i].iat[k] != 0 && atom[i].bo[k] != 9)
293 {
294 if (atom[atom[i].iat[k]].flags & hbmask || atom[atom[i].iat[k]].mmx_type == 20)
295 {
296 nonbond.vrad[ita][itb] -= 0.25;
297 nonbond.vrad[itb][ita] -= 0.25;
298 break;
299 }
300 }
301 }
302 }
303 }
304 nonbond.npair++;
305 }
306 }
307 }
308 for (i=1; i <= natom ; i++)
309 {
310 ita = atom[i].type;
311 if (vdw1.rad[ita] == 0)
312 {
313 rii = vdw1.a[ita]*pow(vdw1.alpha[ita],0.25);
314 vdw1.rad[ita] = 0.5*rii;
315 }
316 }
317 }