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

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