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 User Rev File contents
1 tjod 3 #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 gilbertke 103 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 tjod 3
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 gilbertke 103 int i, j, k,field;
40 tjod 3 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 gilbertke 103 field = get_field();
48 tjod 3 // mm3 setup
49 gilbertke 103 if (field == MM3)
50 tjod 3 {
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 gilbertke 103 ita = atom.type[i];
81     ianum = atom.atomnum[i];
82 tjod 3 for (j= i+1; j <= natom; j++)
83     {
84 gilbertke 103 itb = atom.type[j];
85     ibnum = atom.atomnum[j];
86 tjod 3 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 gilbertke 103 if (field == MMFF94)
98 tjod 3 {
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 gilbertke 103 if ((atom.atomnum[i] == 1 && ita != 5) || (atom.atomnum[j] ==1 && itb != 5) ) // polar hydrogens
103 tjod 3 {
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 gilbertke 103 if (strcmp(get_radiustype(),"SIGMA") == 0)
144 tjod 3 {
145     rad_ia *= 1.122462048;
146     rad_ib *= 1.122462048;
147     }
148 gilbertke 103 if (strcmp(get_radiussize(),"DIAMETER") == 0)
149 tjod 3 {
150     rad_ia /= 2.0;
151     rad_ib /= 2.0;
152     }
153     eps1 = vdw1.eps[ita];
154     eps2 = vdw1.eps[itb];
155    
156 gilbertke 103 if (strcmp(get_radiusrule(),"ARITHMETIC") == 0)
157 tjod 3 {
158     nonbond.vrad[ita][itb] = rad_ia + rad_ib;
159     nonbond.vrad[itb][ita] = rad_ia + rad_ib;
160 gilbertke 103 } else if (strcmp(get_radiusrule(),"GEOMETRIC") == 0)
161 tjod 3 {
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 gilbertke 103 } else if (strcmp(get_radiusrule(),"CUBIC-MEAN") == 0)
165 tjod 3 {
166 gilbertke 103 if ((atom.atomnum[i] == 1 && ita != 5) || (atom.atomnum[j] ==1 && itb != 5) ) // polar hydrogens
167 tjod 3 {
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 gilbertke 103 if (strcmp(get_epsrule(),"ARITHMETIC") == 0)
182 tjod 3 {
183     nonbond.veps[ita][itb] = 0.5*(eps1 + eps2);
184     nonbond.veps[itb][ita] = 0.5*(eps1 + eps2);
185 gilbertke 103 } else if (strcmp(get_epsrule(),"GEOMETRIC") == 0)
186 tjod 3 {
187     nonbond.veps[ita][itb] = sqrt( eps1*eps2 );
188     nonbond.veps[itb][ita] = sqrt( eps1*eps2 );
189 gilbertke 103 } else if (strcmp(get_epsrule(),"HARMONIC") == 0)
190 tjod 3 {
191     nonbond.veps[ita][itb] = 2.0* (eps1*eps2)/(eps1+eps2);
192     nonbond.veps[itb][ita] = 2.0* (eps1*eps2)/(eps1+eps2);
193 gilbertke 103 } else if (strcmp(get_epsrule(),"HHG") == 0)
194 tjod 3 {
195 gilbertke 103 if ((atom.atomnum[i] == 1 && ita != 5) || (atom.atomnum[j] ==1 && itb != 5) ) // polar hydrogens
196 tjod 3 {
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 gilbertke 103 if (field == MM3)
213 tjod 3 {
214     nonbond.ipif[ita][itb] = 1;
215     nonbond.ipif[itb][ita] = 1;
216 gilbertke 103 } else if (field == MMX)
217 tjod 3 {
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 gilbertke 103 numeral(atom.type[i],pa,3);
234     numeral(atom.type[j],pb,3);
235 tjod 3
236 gilbertke 103 if (atom.type[i] < atom.type[j])
237 tjod 3 {
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 gilbertke 103 ita = atom.type[i];
281 tjod 3 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     }