ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/freemol/trunk/src/mengine/src/ehal.c
Revision: 27
Committed: Tue Jul 8 19:10:25 2008 UTC (11 years, 10 months ago) by tjod
File size: 11286 byte(s)
Log Message:
move mengine src to reflect new freemol directory structure

Line User Rev File contents
1 tjod 3 #define EXTERN extern
2    
3     #include "pcwin.h"
4     #include "pcmod.h"
5     #include "utility.h"
6    
7     #include "energies.h"
8     #include "derivs.h"
9     #include "hess.h"
10     #include "torsions.h"
11     #include "nonbond.h"
12     #include "field.h"
13     #include "cutoffs.h"
14     #include "attached.h"
15     #include "pot.h"
16     #include "solv.h"
17    
18     EXTERN struct t_units {
19     double bndunit, cbnd, qbnd;
20     double angunit, cang, qang, pang, sang, aaunit;
21     double stbnunit, ureyunit, torsunit, storunit, v14scale;
22     double aterm, bterm, cterm, dielec, chgscale;
23     } units;
24     EXTERN struct t_minim_values {
25     int iprint, ndc, nconst;
26     float dielc;
27     } minim_values;
28    
29     int ishbond(int, int);
30     int iscoord_bond(int, int);
31    
32     EXTERN int *skip;
33    
34     void image(double *, double *, double *, int);
35    
36     void ehal()
37     {
38     int i,ia, ib, ita, itb, it, j, k;
39     double xi,yi,zi,xr,yr,zr;
40     double xk,yk,zk;
41     double e,rv,eps;
42     double rik,rik2,rik7;
43     double sigma, kappa, kappa7;
44     double rv7,rho, tau;
45     double cutoff;
46    
47     cutoff = cutoffs.vdwcut*cutoffs.vdwcut;
48     energies.evdw = 0.0;
49     energies.e14 = 0.0;
50     sigma = 1.12;
51     kappa = 1.07;
52     kappa7 = pow(kappa,7.0);
53    
54     if (minim_values.iprint)
55     {
56     fprintf(pcmoutfile,"\nVDW Terms - MMFF Buffered 14-7 Potential\n");
57     fprintf(pcmoutfile," At1 At2 Rik Radius Eps Evdw\n");
58     }
59     for (i=1; i <= natom; i++)
60     skip[i] = 0;
61     for (i=1; i < natom; i++)
62     {
63     for (j=0; j < MAXIAT; j++)
64     if (atom[i].iat[j] != 0 && atom[i].bo[j] != 9)
65     skip[atom[i].iat[j]] = i;
66     for (j=0; j < attached.n13[i]; j++)
67     skip[attached.i13[j][i]] = i;
68     for(k=0; k < attached.n14[i]; k++)
69     skip[attached.i14[k][i]] = -i;
70    
71     for(j=i+1; j <= natom; j++)
72     {
73     if (skip[j] != i)
74     {
75     ia = i;
76     ib = j;
77     if (atom[ia].use || atom[ib].use)
78     {
79     ita = atom[ia].type;
80     itb = atom[ib].type;
81     if (ita <= itb)
82     it = 1000*ita + itb;
83     else
84     it = 1000*itb + ita;
85    
86     xi = atom[ia].x;
87     yi = atom[ia].y;
88     zi = atom[ia].z;
89     xk = atom[ib].x;
90     yk = atom[ib].y;
91     zk = atom[ib].z;
92    
93     xr = xi - xk;
94     yr = yi - yk;
95     zr = zi - zk;
96     rik2 = xr*xr + yr*yr + zr*zr;
97     if (rik2 < cutoff)
98     {
99     rv = nonbond.vrad[ita][itb];
100     eps = nonbond.veps[ita][itb];
101     if (skip[j] == -i)
102     eps /= units.v14scale;
103     rik = sqrt(rik2);
104     rv7 = pow(rv,7.0);
105     rik7 = pow(rik,7.0);
106     rho = rik7 + (sigma-1.00)*rv7;
107     tau = rik + (kappa-1.00)*rv;
108    
109     e = eps*kappa7*(rv7/pow(tau,7.0))*(sigma*rv7/rho-2.00);
110     if (skip[j] == -i)
111     energies.e14 += e;
112     else
113     energies.evdw += e;
114     atom[ia].energy += e;
115     atom[ib].energy += e;
116     if (minim_values.iprint)
117     fprintf(pcmoutfile,"VDW: %-4d - %-4d %-8.3f %-8.3f %-8.3f = %8.4f\n",ia, ib, sqrt(rik2),rv,eps,e);
118     }
119     }
120     }
121     }
122     }
123     }
124    
125     void ehal1()
126     {
127     int i,ia, ib, j, k, ita, itb, it;
128     double xi,yi,zi,xr,yr,zr;
129     double xk,yk,zk;
130     double e,rv,eps;
131     double rik,rik2;
132     double rik6,rik7;
133     double sigma, kappa, kappa7;
134     double rv7,rho, tau, rv14;
135     double dedx,dedy,dedz,de;
136     double cutoff;
137    
138     cutoff = cutoffs.vdwcut*cutoffs.vdwcut;
139    
140     energies.evdw = 0.0;
141     energies.e14 = 0.0;
142    
143     sigma = 1.12;
144     kappa = 1.07;
145     kappa7 = pow(kappa,7.0);
146    
147     for (i=1; i <= natom; i++)
148     {
149     deriv.devdw[i][0] = 0.0;
150     deriv.devdw[i][1] = 0.0;
151     deriv.devdw[i][2] = 0.0;
152     deriv.de14[i][0] = 0.0;
153     deriv.de14[i][1] = 0.0;
154     deriv.de14[i][2] = 0.0;
155     }
156    
157     for (i=1; i <= natom; i++)
158     skip[i] = 0;
159    
160     for (i=1; i < natom; i++)
161     {
162     for (j=0; j < MAXIAT; j++)
163     if (atom[i].iat[j] != 0 && atom[i].bo[j] != 9)
164     skip[atom[i].iat[j]] = i;
165     for (j=0; j < attached.n13[i]; j++)
166     skip[attached.i13[j][i]] = i;
167     for(k=0; k < attached.n14[i]; k++)
168     skip[attached.i14[k][i]] = -i;
169     ia = i;
170     ita = atom[i].type;
171    
172     for(j=i+1; j <= natom; j++)
173     {
174     if (atom[i].use || atom[j].use)
175     {
176     if (skip[j] != i)
177     {
178     ib = j;
179     itb = atom[j].type;
180     if (ita <= itb)
181     it = 1000*ita + itb;
182     else
183     it = 1000*itb + ita;
184    
185     xi = atom[ia].x;
186     yi = atom[ia].y;
187     zi = atom[ia].z;
188     xk = atom[ib].x;
189     yk = atom[ib].y;
190     zk = atom[ib].z;
191    
192     xr = xi - xk;
193     yr = yi - yk;
194     zr = zi - zk;
195     rik2 = xr*xr + yr*yr + zr*zr;
196     if (rik2 < cutoff)
197     {
198     rv = nonbond.vrad[ita][itb];
199     eps = nonbond.veps[ita][itb];
200     if (skip[j] == -i)
201     eps /= units.v14scale;
202     rik = sqrt(rik2);
203     rv7 = pow(rv,7.0);
204     rv14 = rv7*rv7;
205     rik6 = rik2*rik2*rik2;
206     rik7 = rik6*rik;
207     rho = rik7 + (sigma-1.00)*rv7;
208     tau = rik + (kappa-1.00)*rv;
209    
210     e = eps*kappa7*(rv7/pow(tau,7.0))*(sigma*rv7/rho-2.00);
211    
212     de = -7.00*eps*kappa7*(rv7/pow(tau,8.00))*(sigma*rv7/rho-2.00)
213     -7.00*eps*kappa7*sigma*rv14*rik6/(rho*rho*pow(tau,7.00));
214    
215     de /= rik;
216     dedx = de*xr;
217     dedy = de*yr;
218     dedz = de*zr;
219    
220     if (skip[j] == -i)
221     {
222     energies.e14 += e;
223     deriv.de14[ia][0] += dedx;
224     deriv.de14[ia][1] += dedy;
225     deriv.de14[ia][2] += dedz;
226     deriv.de14[ib][0] -= dedx;
227     deriv.de14[ib][1] -= dedy;
228     deriv.de14[ib][2] -= dedz;
229     } else
230     {
231     energies.evdw += e;
232     deriv.devdw[ia][0] += dedx;
233     deriv.devdw[ia][1] += dedy;
234     deriv.devdw[ia][2] += dedz;
235     deriv.devdw[ib][0] -= dedx;
236     deriv.devdw[ib][1] -= dedy;
237     deriv.devdw[ib][2] -= dedz;
238     }
239     virial.virx += xr*dedx;
240     virial.viry += yr*dedy;
241     virial.virz += zr*dedz;
242     }
243     }
244     }
245     }
246     }
247     }
248    
249     void ehal2(int iatom)
250     {
251     int j,ia, ib, k, ita, itb, it, kk;
252     double xi,yi,zi,xr,yr,zr;
253     double xk,yk,zk;
254     double e,rv,eps;
255     double rik,rik2, rik5;
256     double rik6,rik7, rik12;
257     double sigma, kappa, kappa7;
258     double rv7,rho, tau, rv14;
259     double de,d2e;
260     double d2edx,d2edy,d2edz,term[3][3];
261     double cutoff;
262    
263     cutoff = cutoffs.vdwcut*cutoffs.vdwcut;
264    
265     sigma = 1.12;
266     kappa = 1.07;
267     kappa7 = pow(kappa,7.0);
268    
269     for (j=1; j <= natom; j++)
270     skip[j] = 0;
271    
272     skip[iatom] = iatom;
273     ita = atom[iatom].type;
274     for (k=0; k < MAXIAT; k++)
275     {
276     if (atom[iatom].iat[k] != 0 && atom[iatom].bo[k] != 9)
277     skip[atom[iatom].iat[k]] = iatom;
278     }
279     for(k=0; k < attached.n13[iatom]; k++)
280     skip[attached.i13[k][iatom]] = iatom;
281     for(k=0; k < attached.n14[iatom]; k++)
282     skip[attached.i14[k][iatom]] = -iatom;
283    
284     ia = iatom;
285     ita = atom[ia].type;
286     xi = atom[ia].x;
287     yi = atom[ia].y;
288     zi = atom[ia].z;
289    
290     for (j=1; j <= natom; j++)
291     {
292     ib = j;
293     if (skip[j] != iatom)
294     {
295     itb = atom[ib].type;
296     if (ita <= itb)
297     it = 1000*ita + itb;
298     else
299     it = 1000*itb + ita;
300    
301     xk = atom[ib].x;
302     yk = atom[ib].y;
303     zk = atom[ib].z;
304    
305     xr = xi - xk;
306     yr = yi - yk;
307     zr = zi - zk;
308     rik2 = xr*xr + yr*yr + zr*zr;
309     if (rik2 < cutoff)
310     {
311     rv = nonbond.vrad[ita][itb];
312     eps = nonbond.veps[ita][itb];
313     if (skip[j] == -iatom)
314     eps /= units.v14scale;
315    
316     rik = sqrt(rik2);
317     rv7 = pow(rv,7.0);
318     rv14 = rv7*rv7;
319     rik6 = rik2*rik2*rik2;
320     rik5 = rik6/rik;
321     rik7 = rik6*rik;
322     rik12 = rik6*rik6;
323     rho = rik7 + (sigma-1.00)*rv7;
324     tau = rik + (kappa-1.00)*rv;
325    
326     de = -7.00*eps*kappa7*(rv7/pow(tau,8.00))*(sigma*rv7/rho-2.00)
327     -7.00*eps*kappa7*sigma*rv14*rik6/(rho*rho*pow(tau,7.00));
328    
329     d2e = 56.0*eps*kappa7*(rv7/pow(tau,9.00))*(sigma*rv7/rho-2.00)
330     + 98.0*eps*kappa7*sigma*rv14*rik6/(rho*rho*pow(tau,8.00))
331     + 98.0*eps*kappa7*sigma*rv14*rik12/(rho*rho*rho*pow(tau,7.00))
332     - 42.0*eps*kappa7*sigma*rv14*rik5/(rho*rho*pow(tau,7.00));
333    
334     de /= rik;
335     d2e = (d2e-de)/rik2;
336     d2edx = d2e*xr;
337     d2edy = d2e*yr;
338     d2edz = d2e*zr;
339    
340     term[0][0] = d2edx*xr + de;
341     term[1][0] = d2edx*yr;
342     term[2][0] = d2edx*zr;
343     term[0][1] = term[1][0];
344     term[1][1] = d2edy*yr + de;
345     term[2][1] = d2edy*zr;
346     term[0][2] = term[2][0];
347     term[1][2] = term[2][1];
348     term[2][2] = d2edz*zr + de;
349    
350     if (ia == iatom)
351     {
352     for (k=0; k < 3; k++)
353     {
354     hess.hessx[ia][k] += term[k][0];
355     hess.hessy[ia][k] += term[k][1];
356     hess.hessz[ia][k] += term[k][2];
357     hess.hessx[ib][k] -= term[k][0];
358     hess.hessy[ib][k] -= term[k][1];
359     hess.hessz[ib][k] -= term[k][2];
360     }
361     } else
362     {
363     for (k=0; k < 3; k++)
364     {
365     hess.hessx[ia][k] -= term[k][0];
366     hess.hessy[ia][k] -= term[k][1];
367     hess.hessz[ia][k] -= term[k][2];
368     hess.hessx[ib][k] += term[k][0];
369     hess.hessy[ib][k] += term[k][1];
370     hess.hessz[ib][k] += term[k][2];
371     }
372     }
373     }
374     }
375     }
376     }