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