ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/freemol/branches/sync4pymol12/src/mengine/src/estrbnd.c
Revision: 63
Committed: Wed Dec 3 03:46:32 2008 UTC (11 years, 1 month ago) by gilbertke
Original Path: trunk/src/mengine/src/estrbnd.c
File size: 24046 byte(s)
Log Message:
updated read_sdf type_mmx and first pass at using best practices
Line File contents
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 virial.virx += xab*dedxia + xcb*dedxic;
218 virial.viry += yab*dedyia + ycb*dedyic;
219 virial.virz += zab*dedzia + zcb*dedzic;
220 }
221 }
222 }
223 }
224
225 void estrbnd2(int iatom)
226 {
227 int i,j,k,ij;
228 int ia,ib,ic;
229 double dt,dr,angle,dot,cosine;
230 double xia,yia,zia;
231 double xib,yib,zib;
232 double xic,yic,zic;
233 double xab,yab,zab,rab;
234 double xcb,ycb,zcb,rcb;
235 double xp,yp,zp,rp,rp2;
236 double terma,termc;
237 double xrab,yrab,zrab,rab2;
238 double xrcb,yrcb,zrcb,rcb2;
239 double xabp,yabp,zabp;
240 double xcbp,ycbp,zcbp;
241 double ddtdxia,ddtdyia,ddtdzia;
242 double ddtdxib,ddtdyib,ddtdzib;
243 double ddtdxic,ddtdyic,ddtdzic;
244 double ddrdxia,ddrdyia,ddrdzia;
245 double ddrdxib,ddrdyib,ddrdzib;
246 double ddrdxic,ddrdyic,ddrdzic;
247 double dtxiaxia,dtxiayia,dtxiazia;
248 double dtxibxib,dtxibyib,dtxibzib;
249 double dtxicxic,dtxicyic,dtxiczic;
250 double dtyiayia,dtyiazia,dtziazia;
251 double dtyibyib,dtyibzib,dtzibzib;
252 double dtyicyic,dtyiczic,dtziczic;
253 double dtxibxia,dtxibyia,dtxibzia;
254 double dtyibxia,dtyibyia,dtyibzia;
255 double dtzibxia,dtzibyia,dtzibzia;
256 double dtxibxic,dtxibyic,dtxibzic;
257 double dtyibxic,dtyibyic,dtyibzic;
258 double dtzibxic,dtzibyic,dtzibzic;
259 double dtxiaxic,dtxiayic,dtxiazic;
260 double dtyiaxic,dtyiayic,dtyiazic;
261 double dtziaxic,dtziayic,dtziazic;
262 double drxiaxia,drxiayia,drxiazia;
263 double drxibxib,drxibyib,drxibzib;
264 double drxicxic,drxicyic,drxiczic;
265 double dryiayia,dryiazia,drziazia;
266 double dryibyib,dryibzib,drzibzib;
267 double dryicyic,dryiczic,drziczic;
268
269 for (i=0; i < strbnd.nstrbnd; i++)
270 {
271 ij = strbnd.isb[i][0];
272 ia = angles.i13[ij][0];
273 ib = angles.i13[ij][1];
274 ic = angles.i13[ij][2];
275 if (iatom == ia || iatom == ib || iatom == ic)
276 {
277 xia = atom[ia].x;
278 yia = atom[ia].y;
279 zia = atom[ia].z;
280 xib = atom[ib].x;
281 yib = atom[ib].y;
282 zib = atom[ib].z;
283 xic = atom[ic].x;
284 yic = atom[ic].y;
285 zic = atom[ic].z;
286 xab = xia - xib;
287 yab = yia - yib;
288 zab = zia - zib;
289 rab = sqrt(xab*xab + yab*yab + zab*zab);
290 xcb = xic - xib;
291 ycb = yic - yib;
292 zcb = zic - zib;
293 rcb = sqrt(xcb*xcb + ycb*ycb + zcb*zcb);
294 xp = ycb*zab - zcb*yab;
295 yp = zcb*xab - xcb*zab;
296 zp = xcb*yab - ycb*xab;
297 rp = sqrt(xp*xp + yp*yp + zp*zp);
298 if (rp != 0.0)
299 {
300 dot = xab*xcb + yab*ycb + zab*zcb;
301 cosine = dot / (rab*rcb);
302 if (cosine < -1.0)
303 cosine = -1.0;
304 if (cosine > 1.0)
305 cosine = 1.0;
306 angle = radian * acos(cosine);
307
308 dt = angle - angles.anat[ij];
309 terma = -radian / (rab*rab*rp);
310 termc = radian / (rcb*rcb*rp);
311 ddtdxia = terma * (yab*zp-zab*yp);
312 ddtdyia = terma * (zab*xp-xab*zp);
313 ddtdzia = terma * (xab*yp-yab*xp);
314 ddtdxic = termc * (ycb*zp-zcb*yp);
315 ddtdyic = termc * (zcb*xp-xcb*zp);
316 ddtdzic = termc * (xcb*yp-ycb*xp);
317 ddtdxib = -ddtdxia - ddtdxic;
318 ddtdyib = -ddtdyia - ddtdyic;
319 ddtdzib = -ddtdzia - ddtdzic;
320
321 rab2 = 2.0 / (rab*rab);
322 xrab = xab * rab2;
323 yrab = yab * rab2;
324 zrab = zab * rab2;
325 rcb2 = 2.0 / (rcb*rcb);
326 xrcb = xcb * rcb2;
327 yrcb = ycb * rcb2;
328 zrcb = zcb * rcb2;
329 rp2 = 1.0 / (rp*rp);
330 xabp = (yab*zp-zab*yp) * rp2;
331 yabp = (zab*xp-xab*zp) * rp2;
332 zabp = (xab*yp-yab*xp) * rp2;
333 xcbp = (ycb*zp-zcb*yp) * rp2;
334 ycbp = (zcb*xp-xcb*zp) * rp2;
335 zcbp = (xcb*yp-ycb*xp) * rp2;
336
337 dtxiaxia = terma*(xab*xcb-dot) + ddtdxia*(xcbp-xrab);
338 dtxiayia = terma*(zp+yab*xcb) + ddtdxia*(ycbp-yrab);
339 dtxiazia = terma*(zab*xcb-yp) + ddtdxia*(zcbp-zrab);
340 dtyiayia = terma*(yab*ycb-dot) + ddtdyia*(ycbp-yrab);
341 dtyiazia = terma*(xp+zab*ycb) + ddtdyia*(zcbp-zrab);
342 dtziazia = terma*(zab*zcb-dot) + ddtdzia*(zcbp-zrab);
343 dtxicxic = termc*(dot-xab*xcb) - ddtdxic*(xabp+xrcb);
344 dtxicyic = termc*(zp-ycb*xab) - ddtdxic*(yabp+yrcb);
345 dtxiczic = -termc*(yp+zcb*xab) - ddtdxic*(zabp+zrcb);
346 dtyicyic = termc*(dot-yab*ycb) - ddtdyic*(yabp+yrcb);
347 dtyiczic = termc*(xp-zcb*yab) - ddtdyic*(zabp+zrcb);
348 dtziczic = termc*(dot-zab*zcb) - ddtdzic*(zabp+zrcb);
349 dtxiaxic = terma*(yab*yab+zab*zab) - ddtdxia*xabp;
350 dtxiayic = -terma*xab*yab - ddtdxia*yabp;
351 dtxiazic = -terma*xab*zab - ddtdxia*zabp;
352 dtyiaxic = -terma*xab*yab - ddtdyia*xabp;
353 dtyiayic = terma*(xab*xab+zab*zab) - ddtdyia*yabp;
354 dtyiazic = -terma*yab*zab - ddtdyia*zabp;
355 dtziaxic = -terma*xab*zab - ddtdzia*xabp;
356 dtziayic = -terma*yab*zab - ddtdzia*yabp;
357 dtziazic = terma*(xab*xab+yab*yab) - ddtdzia*zabp;
358
359 dtxibxia = -dtxiaxia - dtxiaxic;
360 dtxibyia = -dtxiayia - dtyiaxic;
361 dtxibzia = -dtxiazia - dtziaxic;
362 dtyibxia = -dtxiayia - dtxiayic;
363 dtyibyia = -dtyiayia - dtyiayic;
364 dtyibzia = -dtyiazia - dtziayic;
365 dtzibxia = -dtxiazia - dtxiazic;
366 dtzibyia = -dtyiazia - dtyiazic;
367 dtzibzia = -dtziazia - dtziazic;
368 dtxibxic = -dtxicxic - dtxiaxic;
369 dtxibyic = -dtxicyic - dtxiayic;
370 dtxibzic = -dtxiczic - dtxiazic;
371 dtyibxic = -dtxicyic - dtyiaxic;
372 dtyibyic = -dtyicyic - dtyiayic;
373 dtyibzic = -dtyiczic - dtyiazic;
374 dtzibxic = -dtxiczic - dtziaxic;
375 dtzibyic = -dtyiczic - dtziayic;
376 dtzibzic = -dtziczic - dtziazic;
377 dtxibxib = -dtxibxia - dtxibxic;
378 dtxibyib = -dtxibyia - dtxibyic;
379 dtxibzib = -dtxibzia - dtxibzic;
380 dtyibyib = -dtyibyia - dtyibyic;
381 dtyibzib = -dtyibzia - dtyibzic;
382 dtzibzib = -dtzibzia - dtzibzic;
383
384 j = strbnd.isb[i][1];
385 k = strbnd.isb[i][2];
386 dr = 0.0;
387 terma = 0.0;
388 termc = 0.0;
389 if (j > -1)
390 {
391 terma = units.stbnunit* strbnd.ksb1[i] ;
392 dr = terma * (rab-bonds_ff.bl[j]);
393 terma = terma / rab;
394 }
395 if (k > -1)
396 {
397 termc = units.stbnunit* strbnd.ksb2[i];
398 dr = dr + termc*(rcb-bonds_ff.bl[k]);
399 termc = termc / rcb;
400 }
401 ddrdxia = terma * xab;
402 ddrdyia = terma * yab;
403 ddrdzia = terma * zab;
404 ddrdxic = termc * xcb;
405 ddrdyic = termc * ycb;
406 ddrdzic = termc * zcb;
407 ddrdxib = -ddrdxia - ddrdxic;
408 ddrdyib = -ddrdyia - ddrdyic;
409 ddrdzib = -ddrdzia - ddrdzic;
410
411 xab = xab / rab;
412 yab = yab / rab;
413 zab = zab / rab;
414 xcb = xcb / rcb;
415 ycb = ycb / rcb;
416 zcb = zcb / rcb;
417
418 drxiaxia = terma * (1.0-xab*xab);
419 drxiayia = -terma * xab*yab;
420 drxiazia = -terma * xab*zab;
421 dryiayia = terma * (1.0-yab*yab);
422 dryiazia = -terma * yab*zab;
423 drziazia = terma * (1.0-zab*zab);
424 drxicxic = termc * (1.0-xcb*xcb);
425 drxicyic = -termc * xcb*ycb;
426 drxiczic = -termc * xcb*zcb;
427 dryicyic = termc * (1.0-ycb*ycb);
428 dryiczic = -termc * ycb*zcb;
429 drziczic = termc * (1.0-zcb*zcb);
430 drxibxib = drxiaxia + drxicxic;
431 drxibyib = drxiayia + drxicyic;
432 drxibzib = drxiazia + drxiczic;
433 dryibyib = dryiayia + dryicyic;
434 dryibzib = dryiazia + dryiczic;
435 drzibzib = drziazia + drziczic;
436
437 if (ia == iatom)
438 {
439 hess.hessx[ia][0] += dt*drxiaxia + dr*dtxiaxia + 2.0*ddtdxia*ddrdxia;
440 hess.hessx[ia][1] += dt*drxiayia + dr*dtxiayia + ddtdxia*ddrdyia + ddtdyia*ddrdxia;
441 hess.hessx[ia][2] += dt*drxiazia + dr*dtxiazia + ddtdxia*ddrdzia + ddtdzia*ddrdxia;
442 hess.hessy[ia][0] += dt*drxiayia + dr*dtxiayia + ddtdyia*ddrdxia + ddtdxia*ddrdyia;
443 hess.hessy[ia][1] += dt*dryiayia + dr*dtyiayia + 2.0*ddtdyia*ddrdyia;
444 hess.hessy[ia][2] += dt*dryiazia + dr*dtyiazia + ddtdyia*ddrdzia + ddtdzia*ddrdyia;
445 hess.hessz[ia][0] += dt*drxiazia + dr*dtxiazia + ddtdzia*ddrdxia + ddtdxia*ddrdzia;
446 hess.hessz[ia][1] += dt*dryiazia + dr*dtyiazia + ddtdzia*ddrdyia + ddtdyia*ddrdzia;
447 hess.hessz[ia][2] += dt*drziazia + dr*dtziazia + 2.0*ddtdzia*ddrdzia;
448 hess.hessx[ib][0] += -dt*drxiaxia + dr*dtxibxia + ddtdxia*ddrdxib + ddtdxib*ddrdxia;
449 hess.hessx[ib][1] += -dt*drxiayia + dr*dtxibyia + ddtdxia*ddrdyib + ddtdyib*ddrdxia;
450 hess.hessx[ib][2] += -dt*drxiazia + dr*dtxibzia + ddtdxia*ddrdzib + ddtdzib*ddrdxia;
451 hess.hessy[ib][0] += -dt*drxiayia + dr*dtyibxia + ddtdyia*ddrdxib + ddtdxib*ddrdyia;
452 hess.hessy[ib][1] += -dt*dryiayia + dr*dtyibyia + ddtdyia*ddrdyib + ddtdyib*ddrdyia;
453 hess.hessy[ib][2] += -dt*dryiazia + dr*dtyibzia + ddtdyia*ddrdzib + ddtdzib*ddrdyia;
454 hess.hessz[ib][0] += -dt*drxiazia + dr*dtzibxia + ddtdzia*ddrdxib + ddtdxib*ddrdzia;
455 hess.hessz[ib][1] += -dt*dryiazia + dr*dtzibyia + ddtdzia*ddrdyib + ddtdyib*ddrdzia;
456 hess.hessz[ib][2] += -dt*drziazia + dr*dtzibzia + ddtdzia*ddrdzib + ddtdzib*ddrdzia;
457 hess.hessx[ic][0] += dr*dtxiaxic + ddtdxia*ddrdxic + ddtdxic*ddrdxia;
458 hess.hessx[ic][1] += dr*dtxiayic + ddtdxia*ddrdyic + ddtdyic*ddrdxia;
459 hess.hessx[ic][2] += dr*dtxiazic + ddtdxia*ddrdzic + ddtdzic*ddrdxia;
460 hess.hessy[ic][0] += dr*dtyiaxic + ddtdyia*ddrdxic + ddtdxic*ddrdyia;
461 hess.hessy[ic][1] += dr*dtyiayic + ddtdyia*ddrdyic + ddtdyic*ddrdyia;
462 hess.hessy[ic][2] += dr*dtyiazic + ddtdyia*ddrdzic + ddtdzic*ddrdyia;
463 hess.hessz[ic][0] += dr*dtziaxic + ddtdzia*ddrdxic + ddtdxic*ddrdzia;
464 hess.hessz[ic][1] += dr*dtziayic + ddtdzia*ddrdyic + ddtdyic*ddrdzia;
465 hess.hessz[ic][2] += dr*dtziazic + ddtdzia*ddrdzic + ddtdzic*ddrdzia;
466 } else if (ib == iatom)
467 {
468 hess.hessx[ib][0] += dt*drxibxib + dr*dtxibxib + 2.0*ddtdxib*ddrdxib;
469 hess.hessx[ib][1] += dt*drxibyib + dr*dtxibyib + ddtdxib*ddrdyib + ddtdyib*ddrdxib;
470 hess.hessx[ib][2] += dt*drxibzib + dr*dtxibzib + ddtdxib*ddrdzib + ddtdzib*ddrdxib;
471 hess.hessy[ib][0] += dt*drxibyib + dr*dtxibyib + ddtdyib*ddrdxib + ddtdxib*ddrdyib;
472 hess.hessy[ib][1] += dt*dryibyib + dr*dtyibyib + 2.0*ddtdyib*ddrdyib;
473 hess.hessy[ib][2] += dt*dryibzib + dr*dtyibzib + ddtdyib*ddrdzib + ddtdzib*ddrdyib;
474 hess.hessz[ib][0] += dt*drxibzib + dr*dtxibzib + ddtdzib*ddrdxib + ddtdxib*ddrdzib;
475 hess.hessz[ib][1] += dt*dryibzib + dr*dtyibzib + ddtdzib*ddrdyib + ddtdyib*ddrdzib;
476 hess.hessz[ib][2] += dt*drzibzib + dr*dtzibzib + 2.0*ddtdzib*ddrdzib;
477 hess.hessx[ia][0] += -dt*drxiaxia + dr*dtxibxia + ddtdxib*ddrdxia + ddtdxia*ddrdxib;
478 hess.hessx[ia][1] += -dt*drxiayia + dr*dtxibyia + ddtdxib*ddrdyia + ddtdyia*ddrdxib;
479 hess.hessx[ia][2] += -dt*drxiazia + dr*dtxibzia + ddtdxib*ddrdzia + ddtdzia*ddrdxib;
480 hess.hessy[ia][0] += -dt*drxiayia + dr*dtyibxia + ddtdyib*ddrdxia + ddtdxia*ddrdyib;
481 hess.hessy[ia][1] += -dt*dryiayia + dr*dtyibyia + ddtdyib*ddrdyia + ddtdyia*ddrdyib;
482 hess.hessy[ia][2] += -dt*dryiazia + dr*dtyibzia + ddtdyib*ddrdzia + ddtdzia*ddrdyib;
483 hess.hessz[ia][0] += -dt*drxiazia + dr*dtzibxia + ddtdzib*ddrdxia + ddtdxia*ddrdzib;
484 hess.hessz[ia][1] += -dt*dryiazia + dr*dtzibyia + ddtdzib*ddrdyia + ddtdyia*ddrdzib;
485 hess.hessz[ia][2] += -dt*drziazia + dr*dtzibzia + ddtdzib*ddrdzia + ddtdzia*ddrdzib;
486 hess.hessx[ic][0] += -dt*drxicxic + dr*dtxibxic + ddtdxib*ddrdxic + ddtdxic*ddrdxib;
487 hess.hessx[ic][1] += -dt*drxicyic + dr*dtxibyic + ddtdxib*ddrdyic + ddtdyic*ddrdxib;
488 hess.hessx[ic][2] += -dt*drxiczic + dr*dtxibzic + ddtdxib*ddrdzic + ddtdzic*ddrdxib;
489 hess.hessy[ic][0] += -dt*drxicyic + dr*dtyibxic + ddtdyib*ddrdxic + ddtdxic*ddrdyib;
490 hess.hessy[ic][1] += -dt*dryicyic + dr*dtyibyic + ddtdyib*ddrdyic + ddtdyic*ddrdyib;
491 hess.hessy[ic][2] += -dt*dryiczic + dr*dtyibzic + ddtdyib*ddrdzic + ddtdzic*ddrdyib;
492 hess.hessz[ic][0] += -dt*drxiczic + dr*dtzibxic + ddtdzib*ddrdxic + ddtdxic*ddrdzib;
493 hess.hessz[ic][1] += -dt*dryiczic + dr*dtzibyic + ddtdzib*ddrdyic + ddtdyic*ddrdzib;
494 hess.hessz[ic][2] += -dt*drziczic + dr*dtzibzic + ddtdzib*ddrdzic + ddtdzic*ddrdzib;
495 }else if (ic == iatom)
496 {
497 hess.hessx[ic][0] += dt*drxicxic + dr*dtxicxic + 2.0*ddtdxic*ddrdxic;
498 hess.hessx[ic][1] += dt*drxicyic + dr*dtxicyic + ddtdxic*ddrdyic + ddtdyic*ddrdxic;
499 hess.hessx[ic][2] += dt*drxiczic + dr*dtxiczic + ddtdxic*ddrdzic + ddtdzic*ddrdxic;
500 hess.hessy[ic][0] += dt*drxicyic + dr*dtxicyic + ddtdyic*ddrdxic + ddtdxic*ddrdyic;
501 hess.hessy[ic][1] += dt*dryicyic + dr*dtyicyic + 2.0*ddtdyic*ddrdyic;
502 hess.hessy[ic][2] += dt*dryiczic + dr*dtyiczic + ddtdyic*ddrdzic + ddtdzic*ddrdyic;
503 hess.hessz[ic][0] += dt*drxiczic + dr*dtxiczic + ddtdzic*ddrdxic + ddtdxic*ddrdzic;
504 hess.hessz[ic][1] += dt*dryiczic + dr*dtyiczic + ddtdzic*ddrdyic + ddtdyic*ddrdzic;
505 hess.hessz[ic][2] += dt*drziczic + dr*dtziczic + 2.0*ddtdzic*ddrdzic;
506 hess.hessx[ib][0] += -dt*drxicxic + dr*dtxibxic + ddtdxic*ddrdxib + ddtdxib*ddrdxic;
507 hess.hessx[ib][1] += -dt*drxicyic + dr*dtxibyic + ddtdxic*ddrdyib + ddtdyib*ddrdxic;
508 hess.hessx[ib][2] += -dt*drxiczic + dr*dtxibzic + ddtdxic*ddrdzib + ddtdzib*ddrdxic;
509 hess.hessy[ib][0] += -dt*drxicyic + dr*dtyibxic + ddtdyic*ddrdxib + ddtdxib*ddrdyic;
510 hess.hessy[ib][1] += -dt*dryicyic + dr*dtyibyic + ddtdyic*ddrdyib + ddtdyib*ddrdyic;
511 hess.hessy[ib][2] += -dt*dryiczic + dr*dtyibzic + ddtdyic*ddrdzib + ddtdzib*ddrdyic;
512 hess.hessz[ib][0] += -dt*drxiczic + dr*dtzibxic + ddtdzic*ddrdxib + ddtdxib*ddrdzic;
513 hess.hessz[ib][1] += -dt*dryiczic + dr*dtzibyic + ddtdzic*ddrdyib + ddtdyib*ddrdzic;
514 hess.hessz[ib][2] += -dt*drziczic + dr*dtzibzic + ddtdzic*ddrdzib + ddtdzib*ddrdzic;
515 hess.hessx[ia][0] += dr*dtxiaxic + ddtdxic*ddrdxia + ddtdxia*ddrdxic;
516 hess.hessx[ia][1] += dr*dtyiaxic + ddtdxic*ddrdyia + ddtdyia*ddrdxic;
517 hess.hessx[ia][2] += dr*dtziaxic + ddtdxic*ddrdzia + ddtdzia*ddrdxic;
518 hess.hessy[ia][0] += dr*dtxiayic + ddtdyic*ddrdxia + ddtdxia*ddrdyic;
519 hess.hessy[ia][1] += dr*dtyiayic + ddtdyic*ddrdyia + ddtdyia*ddrdyic;
520 hess.hessy[ia][2] += dr*dtziayic + ddtdyic*ddrdzia + ddtdzia*ddrdyic;
521 hess.hessz[ia][0] += dr*dtxiazic + ddtdzic*ddrdxia + ddtdxia*ddrdzic;
522 hess.hessz[ia][1] += dr*dtyiazic + ddtdzic*ddrdyia + ddtdyia*ddrdzic;
523 hess.hessz[ia][2] += dr*dtziazic + ddtdzic*ddrdzia + ddtdzia*ddrdzic;
524 }
525 }
526 }
527 }
528 }
529