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