ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/freemol/trunk/src/mengine/src/kvdw.c
Revision: 90
Committed: Fri Jan 16 21:55:34 2009 UTC (12 years, 6 months ago) by gilbertke
File size: 12523 byte(s)
Log Message:
code cleanup
Line User Rev File contents
1 tjod 3 #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     }