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, 9 months ago) by gilbertke
File size: 12545 byte(s)
Log Message:
updated read_sdf type_mmx and first pass at using best practices
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     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     }