ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/freemol/trunk/smi23d/src/mengine/kvdw.c
Revision: 3
Committed: Mon Jun 9 21:38:26 2008 UTC (11 years, 4 months ago) by tjod
File size: 12797 byte(s)
Log Message:
test

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