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, 11 months ago) by tjod
File size: 12797 byte(s)
Log Message:
test

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     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     }