ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/freemol/trunk/smi23d/src/mengine/ehal.c
Revision: 3
Committed: Mon Jun 9 21:38:26 2008 UTC (11 years, 4 months ago) by tjod
File size: 11286 byte(s)
Log Message:
test

Line File contents
1 #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 }