1 |
#define EXTERN extern |
2 |
|
3 |
#include "pcwin.h" |
4 |
#include "pcmod.h" |
5 |
|
6 |
#include "energies.h" |
7 |
#include "angles.h" |
8 |
#include "bonds_ff.h" |
9 |
#include "hess.h" |
10 |
#include "derivs.h" |
11 |
|
12 |
EXTERN struct t_strbnd { |
13 |
int nstrbnd, isb[MAXANG][3]; |
14 |
float ksb1[MAXANG],ksb2[MAXANG]; |
15 |
} strbnd; |
16 |
|
17 |
EXTERN struct t_minim_values { |
18 |
int iprint, ndc, nconst; |
19 |
float dielc; |
20 |
} minim_values; |
21 |
|
22 |
void estrbnd() |
23 |
{ |
24 |
int i,ij,j,k,ia,ib,ic; |
25 |
double e,dr,dt,angle,dot,cosine,dr1; |
26 |
double xia,yia,zia,xib,yib,zib; |
27 |
double xic,yic,zic; |
28 |
double rab,xab,yab,zab; |
29 |
double rcb,xcb,ycb,zcb; |
30 |
|
31 |
energies.estrbnd = 0; |
32 |
if (minim_values.iprint && strbnd.nstrbnd > 0) |
33 |
{ |
34 |
fprintf(pcmlogfile,"\nStretch-Bend Terms\n"); |
35 |
fprintf(pcmlogfile," At1 At2 At3 Angle Anat StbnCst1 StbnCst2 Estb\n"); |
36 |
} |
37 |
|
38 |
for (i=0; i < strbnd.nstrbnd; i++) |
39 |
{ |
40 |
ij = strbnd.isb[i][0]; |
41 |
ia = angles.i13[ij][0]; |
42 |
ib = angles.i13[ij][1]; |
43 |
ic = angles.i13[ij][2]; |
44 |
if (atom[ia].use || atom[ib].use || atom[ic].use) |
45 |
{ |
46 |
xia = atom[ia].x; |
47 |
yia = atom[ia].y; |
48 |
zia = atom[ia].z; |
49 |
xib = atom[ib].x; |
50 |
yib = atom[ib].y; |
51 |
zib = atom[ib].z; |
52 |
xic = atom[ic].x; |
53 |
yic = atom[ic].y; |
54 |
zic = atom[ic].z; |
55 |
|
56 |
xab = xia - xib; |
57 |
yab = yia - yib; |
58 |
zab = zia - zib; |
59 |
rab = sqrt(xab*xab + yab*yab + zab*zab); |
60 |
xcb = xic - xib; |
61 |
ycb = yic - yib; |
62 |
zcb = zic - zib; |
63 |
rcb = sqrt(xcb*xcb + ycb*ycb + zcb*zcb); |
64 |
if (rab*rcb != 0.0) |
65 |
{ |
66 |
dot = xab*xcb + yab*ycb + zab*zcb; |
67 |
cosine = dot / (rab*rcb); |
68 |
if (cosine < -1.0) |
69 |
cosine = -1.0; |
70 |
if (cosine > 1.0) |
71 |
cosine = 1.0; |
72 |
angle = radian * acos(cosine); |
73 |
dt = angle - angles.anat[ij]; |
74 |
j = strbnd.isb[i][1]; |
75 |
k = strbnd.isb[i][2]; |
76 |
dr = 0.0; |
77 |
dr1 = 0.0; |
78 |
if (j > -1) dr = (rab - bonds_ff.bl[j])*strbnd.ksb1[i]; |
79 |
if (k > -1) dr1 = (rcb - bonds_ff.bl[k])*strbnd.ksb2[i]; |
80 |
e = units.stbnunit*dt*(dr+dr1); |
81 |
energies.estrbnd += e; |
82 |
atom[ia].energy += e; |
83 |
atom[ib].energy += e; |
84 |
atom[ic].energy += e; |
85 |
|
86 |
if (minim_values.iprint) |
87 |
fprintf(pcmlogfile,"StrBnd: %2s(%-3d)-%2s(%-3d)-%2s(%-3d) : %-8.3f %-8.3f %-8.3f %-8.3f = %-8.4f\n", |
88 |
atom[ia].name, ia,atom[ib].name, ib,atom[ic].name, ic,angle,angles.anat[ij], |
89 |
strbnd.ksb1[i],strbnd.ksb2[i], e); |
90 |
} |
91 |
} |
92 |
} |
93 |
} |
94 |
|
95 |
void estrbnd1() |
96 |
{ |
97 |
int i,ij,j,k,ia,ib,ic; |
98 |
double e,dr,dr1,dt,angle,dot,cosine; |
99 |
double xia,yia,zia,xib,yib,zib; |
100 |
double xic,yic,zic; |
101 |
double rab,xab,yab,zab; |
102 |
double rcb,xcb,ycb,zcb; |
103 |
double xp,yp,zp,rp; |
104 |
double dedxia,dedyia,dedzia; |
105 |
double dedxib,dedyib,dedzib; |
106 |
double dedxic,dedyic,dedzic; |
107 |
double term,terma,termc, rab2, rcb2; |
108 |
double ksb1, ksb2; |
109 |
double dtemp; |
110 |
|
111 |
energies.estrbnd = 0; |
112 |
for (i=0; i <= natom; i++) |
113 |
{ |
114 |
deriv.destb[i][0] = 0.0; |
115 |
deriv.destb[i][1] = 0.0; |
116 |
deriv.destb[i][2] = 0.0; |
117 |
} |
118 |
|
119 |
for (i=0; i < strbnd.nstrbnd; i++) |
120 |
{ |
121 |
ij = strbnd.isb[i][0]; |
122 |
ia = angles.i13[ij][0]; |
123 |
ib = angles.i13[ij][1]; |
124 |
ic = angles.i13[ij][2]; |
125 |
if (atom[ia].use || atom[ib].use || atom[ic].use) |
126 |
{ |
127 |
xia = atom[ia].x; |
128 |
yia = atom[ia].y; |
129 |
zia = atom[ia].z; |
130 |
xib = atom[ib].x; |
131 |
yib = atom[ib].y; |
132 |
zib = atom[ib].z; |
133 |
xic = atom[ic].x; |
134 |
yic = atom[ic].y; |
135 |
zic = atom[ic].z; |
136 |
|
137 |
xab = xia - xib; |
138 |
yab = yia - yib; |
139 |
zab = zia - zib; |
140 |
rab2 = xab*xab + yab*yab + zab*zab; |
141 |
rab = sqrt(xab*xab + yab*yab + zab*zab); |
142 |
xcb = xic - xib; |
143 |
ycb = yic - yib; |
144 |
zcb = zic - zib; |
145 |
rcb2 = xcb*xcb + ycb*ycb + zcb*zcb; |
146 |
rcb = sqrt(xcb*xcb + ycb*ycb + zcb*zcb); |
147 |
xp = ycb*zab - zcb*yab; |
148 |
yp = zcb*xab - xcb*zab; |
149 |
zp = xcb*yab - ycb*xab; |
150 |
rp = sqrt(xp*xp + yp*yp + zp*zp); |
151 |
if (rp != 0.0) |
152 |
{ |
153 |
dot = xab*xcb + yab*ycb + zab*zcb; |
154 |
cosine = dot / (rab*rcb); |
155 |
if (cosine < -1.0) |
156 |
cosine = -1.0; |
157 |
if (cosine > 1.0) |
158 |
cosine = 1.0; |
159 |
angle = radian * acos(cosine); |
160 |
dt = angle - angles.anat[ij]; |
161 |
|
162 |
dtemp = dot*dot/(rab2*rcb2); |
163 |
if (fabs(dtemp) >= 1.0) |
164 |
dtemp = 0.9999*dtemp/fabs(dtemp); |
165 |
termc = sqrt(1.00 - dtemp); |
166 |
terma = rab*rcb; |
167 |
j = strbnd.isb[i][1]; |
168 |
k = strbnd.isb[i][2]; |
169 |
term = -units.stbnunit*radian; |
170 |
if (j < 0) |
171 |
{ |
172 |
dr = 0.0; |
173 |
ksb1 = 0.0; |
174 |
} else |
175 |
{ |
176 |
dr = strbnd.ksb1[i]*(rab - bonds_ff.bl[j]); |
177 |
ksb1 = strbnd.ksb1[i]; |
178 |
} |
179 |
if (k < 0) |
180 |
{ |
181 |
dr1 = 0.0; |
182 |
ksb2 = 0.0; |
183 |
} else |
184 |
{ |
185 |
dr1 = strbnd.ksb2[i]*(rcb - bonds_ff.bl[k]); |
186 |
ksb2 = strbnd.ksb2[i]; |
187 |
} |
188 |
e = units.stbnunit*dt*(dr+dr1); |
189 |
dedxia = (term *(xcb/terma - (xab*dot)/(rab2*terma))*(dr+dr1))/ termc + (units.stbnunit*dt*ksb1*xab)/rab; |
190 |
dedyia = (term *(ycb/terma - (yab*dot)/(rab2*terma))*(dr+dr1))/ termc + (units.stbnunit*dt*ksb1*yab)/rab; |
191 |
dedzia = (term *(zcb/terma - (zab*dot)/(rab2*terma))*(dr+dr1))/ termc + (units.stbnunit*dt*ksb1*zab)/rab; |
192 |
|
193 |
dedxic = (term *(xab/terma - (xcb*dot)/(rcb2*terma))*(dr+dr1))/ termc + (units.stbnunit*dt*ksb2*xcb)/rcb; |
194 |
dedyic = (term *(yab/terma - (ycb*dot)/(rcb2*terma))*(dr+dr1))/ termc + (units.stbnunit*dt*ksb2*ycb)/rcb; |
195 |
dedzic = (term *(zab/terma - (zcb*dot)/(rcb2*terma))*(dr+dr1))/ termc + (units.stbnunit*dt*ksb2*zcb)/rcb; |
196 |
|
197 |
dedxib = ( (term /(rab*rcb))*( (xcb*dot)/rcb2 + (xab*dot)/rab2 + (2.0*xib-xia-xic))*(dr+dr1))/termc + |
198 |
units.stbnunit*dt*( -(ksb1*xab)/rab - ksb2*xcb/rcb); |
199 |
dedyib = ( (term /(rab*rcb))*( (ycb*dot)/rcb2 + (yab*dot)/rab2 + (2.0*yib-yia-yic))*(dr+dr1))/termc + |
200 |
units.stbnunit*dt*( -(ksb1*yab)/rab - ksb2*ycb/rcb); |
201 |
dedzib = ( (term /(rab*rcb))*( (zcb*dot)/rcb2 + (zab*dot)/rab2 + (2.0*zib-zia-zic))*(dr+dr1))/termc + |
202 |
units.stbnunit*dt*( -(ksb1*zab)/rab - ksb2*zcb/rcb); |
203 |
|
204 |
energies.estrbnd += e; |
205 |
deriv.destb[ia][0] += dedxia; |
206 |
deriv.destb[ia][1] += dedyia; |
207 |
deriv.destb[ia][2] += dedzia; |
208 |
|
209 |
deriv.destb[ib][0] += dedxib; |
210 |
deriv.destb[ib][1] += dedyib; |
211 |
deriv.destb[ib][2] += dedzib; |
212 |
|
213 |
deriv.destb[ic][0] += dedxic; |
214 |
deriv.destb[ic][1] += dedyic; |
215 |
deriv.destb[ic][2] += dedzic; |
216 |
|
217 |
} |
218 |
} |
219 |
} |
220 |
} |
221 |
|
222 |
void estrbnd2(int iatom) |
223 |
{ |
224 |
int i,j,k,ij; |
225 |
int ia,ib,ic; |
226 |
double dt,dr,angle,dot,cosine; |
227 |
double xia,yia,zia; |
228 |
double xib,yib,zib; |
229 |
double xic,yic,zic; |
230 |
double xab,yab,zab,rab; |
231 |
double xcb,ycb,zcb,rcb; |
232 |
double xp,yp,zp,rp,rp2; |
233 |
double terma,termc; |
234 |
double xrab,yrab,zrab,rab2; |
235 |
double xrcb,yrcb,zrcb,rcb2; |
236 |
double xabp,yabp,zabp; |
237 |
double xcbp,ycbp,zcbp; |
238 |
double ddtdxia,ddtdyia,ddtdzia; |
239 |
double ddtdxib,ddtdyib,ddtdzib; |
240 |
double ddtdxic,ddtdyic,ddtdzic; |
241 |
double ddrdxia,ddrdyia,ddrdzia; |
242 |
double ddrdxib,ddrdyib,ddrdzib; |
243 |
double ddrdxic,ddrdyic,ddrdzic; |
244 |
double dtxiaxia,dtxiayia,dtxiazia; |
245 |
double dtxibxib,dtxibyib,dtxibzib; |
246 |
double dtxicxic,dtxicyic,dtxiczic; |
247 |
double dtyiayia,dtyiazia,dtziazia; |
248 |
double dtyibyib,dtyibzib,dtzibzib; |
249 |
double dtyicyic,dtyiczic,dtziczic; |
250 |
double dtxibxia,dtxibyia,dtxibzia; |
251 |
double dtyibxia,dtyibyia,dtyibzia; |
252 |
double dtzibxia,dtzibyia,dtzibzia; |
253 |
double dtxibxic,dtxibyic,dtxibzic; |
254 |
double dtyibxic,dtyibyic,dtyibzic; |
255 |
double dtzibxic,dtzibyic,dtzibzic; |
256 |
double dtxiaxic,dtxiayic,dtxiazic; |
257 |
double dtyiaxic,dtyiayic,dtyiazic; |
258 |
double dtziaxic,dtziayic,dtziazic; |
259 |
double drxiaxia,drxiayia,drxiazia; |
260 |
double drxibxib,drxibyib,drxibzib; |
261 |
double drxicxic,drxicyic,drxiczic; |
262 |
double dryiayia,dryiazia,drziazia; |
263 |
double dryibyib,dryibzib,drzibzib; |
264 |
double dryicyic,dryiczic,drziczic; |
265 |
|
266 |
for (i=0; i < strbnd.nstrbnd; i++) |
267 |
{ |
268 |
ij = strbnd.isb[i][0]; |
269 |
ia = angles.i13[ij][0]; |
270 |
ib = angles.i13[ij][1]; |
271 |
ic = angles.i13[ij][2]; |
272 |
if (iatom == ia || iatom == ib || iatom == ic) |
273 |
{ |
274 |
xia = atom[ia].x; |
275 |
yia = atom[ia].y; |
276 |
zia = atom[ia].z; |
277 |
xib = atom[ib].x; |
278 |
yib = atom[ib].y; |
279 |
zib = atom[ib].z; |
280 |
xic = atom[ic].x; |
281 |
yic = atom[ic].y; |
282 |
zic = atom[ic].z; |
283 |
xab = xia - xib; |
284 |
yab = yia - yib; |
285 |
zab = zia - zib; |
286 |
rab = sqrt(xab*xab + yab*yab + zab*zab); |
287 |
xcb = xic - xib; |
288 |
ycb = yic - yib; |
289 |
zcb = zic - zib; |
290 |
rcb = sqrt(xcb*xcb + ycb*ycb + zcb*zcb); |
291 |
xp = ycb*zab - zcb*yab; |
292 |
yp = zcb*xab - xcb*zab; |
293 |
zp = xcb*yab - ycb*xab; |
294 |
rp = sqrt(xp*xp + yp*yp + zp*zp); |
295 |
if (rp != 0.0) |
296 |
{ |
297 |
dot = xab*xcb + yab*ycb + zab*zcb; |
298 |
cosine = dot / (rab*rcb); |
299 |
if (cosine < -1.0) |
300 |
cosine = -1.0; |
301 |
if (cosine > 1.0) |
302 |
cosine = 1.0; |
303 |
angle = radian * acos(cosine); |
304 |
|
305 |
dt = angle - angles.anat[ij]; |
306 |
terma = -radian / (rab*rab*rp); |
307 |
termc = radian / (rcb*rcb*rp); |
308 |
ddtdxia = terma * (yab*zp-zab*yp); |
309 |
ddtdyia = terma * (zab*xp-xab*zp); |
310 |
ddtdzia = terma * (xab*yp-yab*xp); |
311 |
ddtdxic = termc * (ycb*zp-zcb*yp); |
312 |
ddtdyic = termc * (zcb*xp-xcb*zp); |
313 |
ddtdzic = termc * (xcb*yp-ycb*xp); |
314 |
ddtdxib = -ddtdxia - ddtdxic; |
315 |
ddtdyib = -ddtdyia - ddtdyic; |
316 |
ddtdzib = -ddtdzia - ddtdzic; |
317 |
|
318 |
rab2 = 2.0 / (rab*rab); |
319 |
xrab = xab * rab2; |
320 |
yrab = yab * rab2; |
321 |
zrab = zab * rab2; |
322 |
rcb2 = 2.0 / (rcb*rcb); |
323 |
xrcb = xcb * rcb2; |
324 |
yrcb = ycb * rcb2; |
325 |
zrcb = zcb * rcb2; |
326 |
rp2 = 1.0 / (rp*rp); |
327 |
xabp = (yab*zp-zab*yp) * rp2; |
328 |
yabp = (zab*xp-xab*zp) * rp2; |
329 |
zabp = (xab*yp-yab*xp) * rp2; |
330 |
xcbp = (ycb*zp-zcb*yp) * rp2; |
331 |
ycbp = (zcb*xp-xcb*zp) * rp2; |
332 |
zcbp = (xcb*yp-ycb*xp) * rp2; |
333 |
|
334 |
dtxiaxia = terma*(xab*xcb-dot) + ddtdxia*(xcbp-xrab); |
335 |
dtxiayia = terma*(zp+yab*xcb) + ddtdxia*(ycbp-yrab); |
336 |
dtxiazia = terma*(zab*xcb-yp) + ddtdxia*(zcbp-zrab); |
337 |
dtyiayia = terma*(yab*ycb-dot) + ddtdyia*(ycbp-yrab); |
338 |
dtyiazia = terma*(xp+zab*ycb) + ddtdyia*(zcbp-zrab); |
339 |
dtziazia = terma*(zab*zcb-dot) + ddtdzia*(zcbp-zrab); |
340 |
dtxicxic = termc*(dot-xab*xcb) - ddtdxic*(xabp+xrcb); |
341 |
dtxicyic = termc*(zp-ycb*xab) - ddtdxic*(yabp+yrcb); |
342 |
dtxiczic = -termc*(yp+zcb*xab) - ddtdxic*(zabp+zrcb); |
343 |
dtyicyic = termc*(dot-yab*ycb) - ddtdyic*(yabp+yrcb); |
344 |
dtyiczic = termc*(xp-zcb*yab) - ddtdyic*(zabp+zrcb); |
345 |
dtziczic = termc*(dot-zab*zcb) - ddtdzic*(zabp+zrcb); |
346 |
dtxiaxic = terma*(yab*yab+zab*zab) - ddtdxia*xabp; |
347 |
dtxiayic = -terma*xab*yab - ddtdxia*yabp; |
348 |
dtxiazic = -terma*xab*zab - ddtdxia*zabp; |
349 |
dtyiaxic = -terma*xab*yab - ddtdyia*xabp; |
350 |
dtyiayic = terma*(xab*xab+zab*zab) - ddtdyia*yabp; |
351 |
dtyiazic = -terma*yab*zab - ddtdyia*zabp; |
352 |
dtziaxic = -terma*xab*zab - ddtdzia*xabp; |
353 |
dtziayic = -terma*yab*zab - ddtdzia*yabp; |
354 |
dtziazic = terma*(xab*xab+yab*yab) - ddtdzia*zabp; |
355 |
|
356 |
dtxibxia = -dtxiaxia - dtxiaxic; |
357 |
dtxibyia = -dtxiayia - dtyiaxic; |
358 |
dtxibzia = -dtxiazia - dtziaxic; |
359 |
dtyibxia = -dtxiayia - dtxiayic; |
360 |
dtyibyia = -dtyiayia - dtyiayic; |
361 |
dtyibzia = -dtyiazia - dtziayic; |
362 |
dtzibxia = -dtxiazia - dtxiazic; |
363 |
dtzibyia = -dtyiazia - dtyiazic; |
364 |
dtzibzia = -dtziazia - dtziazic; |
365 |
dtxibxic = -dtxicxic - dtxiaxic; |
366 |
dtxibyic = -dtxicyic - dtxiayic; |
367 |
dtxibzic = -dtxiczic - dtxiazic; |
368 |
dtyibxic = -dtxicyic - dtyiaxic; |
369 |
dtyibyic = -dtyicyic - dtyiayic; |
370 |
dtyibzic = -dtyiczic - dtyiazic; |
371 |
dtzibxic = -dtxiczic - dtziaxic; |
372 |
dtzibyic = -dtyiczic - dtziayic; |
373 |
dtzibzic = -dtziczic - dtziazic; |
374 |
dtxibxib = -dtxibxia - dtxibxic; |
375 |
dtxibyib = -dtxibyia - dtxibyic; |
376 |
dtxibzib = -dtxibzia - dtxibzic; |
377 |
dtyibyib = -dtyibyia - dtyibyic; |
378 |
dtyibzib = -dtyibzia - dtyibzic; |
379 |
dtzibzib = -dtzibzia - dtzibzic; |
380 |
|
381 |
j = strbnd.isb[i][1]; |
382 |
k = strbnd.isb[i][2]; |
383 |
dr = 0.0; |
384 |
terma = 0.0; |
385 |
termc = 0.0; |
386 |
if (j > -1) |
387 |
{ |
388 |
terma = units.stbnunit* strbnd.ksb1[i] ; |
389 |
dr = terma * (rab-bonds_ff.bl[j]); |
390 |
terma = terma / rab; |
391 |
} |
392 |
if (k > -1) |
393 |
{ |
394 |
termc = units.stbnunit* strbnd.ksb2[i]; |
395 |
dr = dr + termc*(rcb-bonds_ff.bl[k]); |
396 |
termc = termc / rcb; |
397 |
} |
398 |
ddrdxia = terma * xab; |
399 |
ddrdyia = terma * yab; |
400 |
ddrdzia = terma * zab; |
401 |
ddrdxic = termc * xcb; |
402 |
ddrdyic = termc * ycb; |
403 |
ddrdzic = termc * zcb; |
404 |
ddrdxib = -ddrdxia - ddrdxic; |
405 |
ddrdyib = -ddrdyia - ddrdyic; |
406 |
ddrdzib = -ddrdzia - ddrdzic; |
407 |
|
408 |
xab = xab / rab; |
409 |
yab = yab / rab; |
410 |
zab = zab / rab; |
411 |
xcb = xcb / rcb; |
412 |
ycb = ycb / rcb; |
413 |
zcb = zcb / rcb; |
414 |
|
415 |
drxiaxia = terma * (1.0-xab*xab); |
416 |
drxiayia = -terma * xab*yab; |
417 |
drxiazia = -terma * xab*zab; |
418 |
dryiayia = terma * (1.0-yab*yab); |
419 |
dryiazia = -terma * yab*zab; |
420 |
drziazia = terma * (1.0-zab*zab); |
421 |
drxicxic = termc * (1.0-xcb*xcb); |
422 |
drxicyic = -termc * xcb*ycb; |
423 |
drxiczic = -termc * xcb*zcb; |
424 |
dryicyic = termc * (1.0-ycb*ycb); |
425 |
dryiczic = -termc * ycb*zcb; |
426 |
drziczic = termc * (1.0-zcb*zcb); |
427 |
drxibxib = drxiaxia + drxicxic; |
428 |
drxibyib = drxiayia + drxicyic; |
429 |
drxibzib = drxiazia + drxiczic; |
430 |
dryibyib = dryiayia + dryicyic; |
431 |
dryibzib = dryiazia + dryiczic; |
432 |
drzibzib = drziazia + drziczic; |
433 |
|
434 |
if (ia == iatom) |
435 |
{ |
436 |
hess.hessx[ia][0] += dt*drxiaxia + dr*dtxiaxia + 2.0*ddtdxia*ddrdxia; |
437 |
hess.hessx[ia][1] += dt*drxiayia + dr*dtxiayia + ddtdxia*ddrdyia + ddtdyia*ddrdxia; |
438 |
hess.hessx[ia][2] += dt*drxiazia + dr*dtxiazia + ddtdxia*ddrdzia + ddtdzia*ddrdxia; |
439 |
hess.hessy[ia][0] += dt*drxiayia + dr*dtxiayia + ddtdyia*ddrdxia + ddtdxia*ddrdyia; |
440 |
hess.hessy[ia][1] += dt*dryiayia + dr*dtyiayia + 2.0*ddtdyia*ddrdyia; |
441 |
hess.hessy[ia][2] += dt*dryiazia + dr*dtyiazia + ddtdyia*ddrdzia + ddtdzia*ddrdyia; |
442 |
hess.hessz[ia][0] += dt*drxiazia + dr*dtxiazia + ddtdzia*ddrdxia + ddtdxia*ddrdzia; |
443 |
hess.hessz[ia][1] += dt*dryiazia + dr*dtyiazia + ddtdzia*ddrdyia + ddtdyia*ddrdzia; |
444 |
hess.hessz[ia][2] += dt*drziazia + dr*dtziazia + 2.0*ddtdzia*ddrdzia; |
445 |
hess.hessx[ib][0] += -dt*drxiaxia + dr*dtxibxia + ddtdxia*ddrdxib + ddtdxib*ddrdxia; |
446 |
hess.hessx[ib][1] += -dt*drxiayia + dr*dtxibyia + ddtdxia*ddrdyib + ddtdyib*ddrdxia; |
447 |
hess.hessx[ib][2] += -dt*drxiazia + dr*dtxibzia + ddtdxia*ddrdzib + ddtdzib*ddrdxia; |
448 |
hess.hessy[ib][0] += -dt*drxiayia + dr*dtyibxia + ddtdyia*ddrdxib + ddtdxib*ddrdyia; |
449 |
hess.hessy[ib][1] += -dt*dryiayia + dr*dtyibyia + ddtdyia*ddrdyib + ddtdyib*ddrdyia; |
450 |
hess.hessy[ib][2] += -dt*dryiazia + dr*dtyibzia + ddtdyia*ddrdzib + ddtdzib*ddrdyia; |
451 |
hess.hessz[ib][0] += -dt*drxiazia + dr*dtzibxia + ddtdzia*ddrdxib + ddtdxib*ddrdzia; |
452 |
hess.hessz[ib][1] += -dt*dryiazia + dr*dtzibyia + ddtdzia*ddrdyib + ddtdyib*ddrdzia; |
453 |
hess.hessz[ib][2] += -dt*drziazia + dr*dtzibzia + ddtdzia*ddrdzib + ddtdzib*ddrdzia; |
454 |
hess.hessx[ic][0] += dr*dtxiaxic + ddtdxia*ddrdxic + ddtdxic*ddrdxia; |
455 |
hess.hessx[ic][1] += dr*dtxiayic + ddtdxia*ddrdyic + ddtdyic*ddrdxia; |
456 |
hess.hessx[ic][2] += dr*dtxiazic + ddtdxia*ddrdzic + ddtdzic*ddrdxia; |
457 |
hess.hessy[ic][0] += dr*dtyiaxic + ddtdyia*ddrdxic + ddtdxic*ddrdyia; |
458 |
hess.hessy[ic][1] += dr*dtyiayic + ddtdyia*ddrdyic + ddtdyic*ddrdyia; |
459 |
hess.hessy[ic][2] += dr*dtyiazic + ddtdyia*ddrdzic + ddtdzic*ddrdyia; |
460 |
hess.hessz[ic][0] += dr*dtziaxic + ddtdzia*ddrdxic + ddtdxic*ddrdzia; |
461 |
hess.hessz[ic][1] += dr*dtziayic + ddtdzia*ddrdyic + ddtdyic*ddrdzia; |
462 |
hess.hessz[ic][2] += dr*dtziazic + ddtdzia*ddrdzic + ddtdzic*ddrdzia; |
463 |
} else if (ib == iatom) |
464 |
{ |
465 |
hess.hessx[ib][0] += dt*drxibxib + dr*dtxibxib + 2.0*ddtdxib*ddrdxib; |
466 |
hess.hessx[ib][1] += dt*drxibyib + dr*dtxibyib + ddtdxib*ddrdyib + ddtdyib*ddrdxib; |
467 |
hess.hessx[ib][2] += dt*drxibzib + dr*dtxibzib + ddtdxib*ddrdzib + ddtdzib*ddrdxib; |
468 |
hess.hessy[ib][0] += dt*drxibyib + dr*dtxibyib + ddtdyib*ddrdxib + ddtdxib*ddrdyib; |
469 |
hess.hessy[ib][1] += dt*dryibyib + dr*dtyibyib + 2.0*ddtdyib*ddrdyib; |
470 |
hess.hessy[ib][2] += dt*dryibzib + dr*dtyibzib + ddtdyib*ddrdzib + ddtdzib*ddrdyib; |
471 |
hess.hessz[ib][0] += dt*drxibzib + dr*dtxibzib + ddtdzib*ddrdxib + ddtdxib*ddrdzib; |
472 |
hess.hessz[ib][1] += dt*dryibzib + dr*dtyibzib + ddtdzib*ddrdyib + ddtdyib*ddrdzib; |
473 |
hess.hessz[ib][2] += dt*drzibzib + dr*dtzibzib + 2.0*ddtdzib*ddrdzib; |
474 |
hess.hessx[ia][0] += -dt*drxiaxia + dr*dtxibxia + ddtdxib*ddrdxia + ddtdxia*ddrdxib; |
475 |
hess.hessx[ia][1] += -dt*drxiayia + dr*dtxibyia + ddtdxib*ddrdyia + ddtdyia*ddrdxib; |
476 |
hess.hessx[ia][2] += -dt*drxiazia + dr*dtxibzia + ddtdxib*ddrdzia + ddtdzia*ddrdxib; |
477 |
hess.hessy[ia][0] += -dt*drxiayia + dr*dtyibxia + ddtdyib*ddrdxia + ddtdxia*ddrdyib; |
478 |
hess.hessy[ia][1] += -dt*dryiayia + dr*dtyibyia + ddtdyib*ddrdyia + ddtdyia*ddrdyib; |
479 |
hess.hessy[ia][2] += -dt*dryiazia + dr*dtyibzia + ddtdyib*ddrdzia + ddtdzia*ddrdyib; |
480 |
hess.hessz[ia][0] += -dt*drxiazia + dr*dtzibxia + ddtdzib*ddrdxia + ddtdxia*ddrdzib; |
481 |
hess.hessz[ia][1] += -dt*dryiazia + dr*dtzibyia + ddtdzib*ddrdyia + ddtdyia*ddrdzib; |
482 |
hess.hessz[ia][2] += -dt*drziazia + dr*dtzibzia + ddtdzib*ddrdzia + ddtdzia*ddrdzib; |
483 |
hess.hessx[ic][0] += -dt*drxicxic + dr*dtxibxic + ddtdxib*ddrdxic + ddtdxic*ddrdxib; |
484 |
hess.hessx[ic][1] += -dt*drxicyic + dr*dtxibyic + ddtdxib*ddrdyic + ddtdyic*ddrdxib; |
485 |
hess.hessx[ic][2] += -dt*drxiczic + dr*dtxibzic + ddtdxib*ddrdzic + ddtdzic*ddrdxib; |
486 |
hess.hessy[ic][0] += -dt*drxicyic + dr*dtyibxic + ddtdyib*ddrdxic + ddtdxic*ddrdyib; |
487 |
hess.hessy[ic][1] += -dt*dryicyic + dr*dtyibyic + ddtdyib*ddrdyic + ddtdyic*ddrdyib; |
488 |
hess.hessy[ic][2] += -dt*dryiczic + dr*dtyibzic + ddtdyib*ddrdzic + ddtdzic*ddrdyib; |
489 |
hess.hessz[ic][0] += -dt*drxiczic + dr*dtzibxic + ddtdzib*ddrdxic + ddtdxic*ddrdzib; |
490 |
hess.hessz[ic][1] += -dt*dryiczic + dr*dtzibyic + ddtdzib*ddrdyic + ddtdyic*ddrdzib; |
491 |
hess.hessz[ic][2] += -dt*drziczic + dr*dtzibzic + ddtdzib*ddrdzic + ddtdzic*ddrdzib; |
492 |
}else if (ic == iatom) |
493 |
{ |
494 |
hess.hessx[ic][0] += dt*drxicxic + dr*dtxicxic + 2.0*ddtdxic*ddrdxic; |
495 |
hess.hessx[ic][1] += dt*drxicyic + dr*dtxicyic + ddtdxic*ddrdyic + ddtdyic*ddrdxic; |
496 |
hess.hessx[ic][2] += dt*drxiczic + dr*dtxiczic + ddtdxic*ddrdzic + ddtdzic*ddrdxic; |
497 |
hess.hessy[ic][0] += dt*drxicyic + dr*dtxicyic + ddtdyic*ddrdxic + ddtdxic*ddrdyic; |
498 |
hess.hessy[ic][1] += dt*dryicyic + dr*dtyicyic + 2.0*ddtdyic*ddrdyic; |
499 |
hess.hessy[ic][2] += dt*dryiczic + dr*dtyiczic + ddtdyic*ddrdzic + ddtdzic*ddrdyic; |
500 |
hess.hessz[ic][0] += dt*drxiczic + dr*dtxiczic + ddtdzic*ddrdxic + ddtdxic*ddrdzic; |
501 |
hess.hessz[ic][1] += dt*dryiczic + dr*dtyiczic + ddtdzic*ddrdyic + ddtdyic*ddrdzic; |
502 |
hess.hessz[ic][2] += dt*drziczic + dr*dtziczic + 2.0*ddtdzic*ddrdzic; |
503 |
hess.hessx[ib][0] += -dt*drxicxic + dr*dtxibxic + ddtdxic*ddrdxib + ddtdxib*ddrdxic; |
504 |
hess.hessx[ib][1] += -dt*drxicyic + dr*dtxibyic + ddtdxic*ddrdyib + ddtdyib*ddrdxic; |
505 |
hess.hessx[ib][2] += -dt*drxiczic + dr*dtxibzic + ddtdxic*ddrdzib + ddtdzib*ddrdxic; |
506 |
hess.hessy[ib][0] += -dt*drxicyic + dr*dtyibxic + ddtdyic*ddrdxib + ddtdxib*ddrdyic; |
507 |
hess.hessy[ib][1] += -dt*dryicyic + dr*dtyibyic + ddtdyic*ddrdyib + ddtdyib*ddrdyic; |
508 |
hess.hessy[ib][2] += -dt*dryiczic + dr*dtyibzic + ddtdyic*ddrdzib + ddtdzib*ddrdyic; |
509 |
hess.hessz[ib][0] += -dt*drxiczic + dr*dtzibxic + ddtdzic*ddrdxib + ddtdxib*ddrdzic; |
510 |
hess.hessz[ib][1] += -dt*dryiczic + dr*dtzibyic + ddtdzic*ddrdyib + ddtdyib*ddrdzic; |
511 |
hess.hessz[ib][2] += -dt*drziczic + dr*dtzibzic + ddtdzic*ddrdzib + ddtdzib*ddrdzic; |
512 |
hess.hessx[ia][0] += dr*dtxiaxic + ddtdxic*ddrdxia + ddtdxia*ddrdxic; |
513 |
hess.hessx[ia][1] += dr*dtyiaxic + ddtdxic*ddrdyia + ddtdyia*ddrdxic; |
514 |
hess.hessx[ia][2] += dr*dtziaxic + ddtdxic*ddrdzia + ddtdzia*ddrdxic; |
515 |
hess.hessy[ia][0] += dr*dtxiayic + ddtdyic*ddrdxia + ddtdxia*ddrdyic; |
516 |
hess.hessy[ia][1] += dr*dtyiayic + ddtdyic*ddrdyia + ddtdyia*ddrdyic; |
517 |
hess.hessy[ia][2] += dr*dtziayic + ddtdyic*ddrdzia + ddtdzia*ddrdyic; |
518 |
hess.hessz[ia][0] += dr*dtxiazic + ddtdzic*ddrdxia + ddtdxia*ddrdzic; |
519 |
hess.hessz[ia][1] += dr*dtyiazic + ddtdzic*ddrdyia + ddtdyia*ddrdzic; |
520 |
hess.hessz[ia][2] += dr*dtziazic + ddtdzic*ddrdzia + ddtdzia*ddrdzic; |
521 |
} |
522 |
} |
523 |
} |
524 |
} |
525 |
} |
526 |
|