ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/freemol/branches/sync4pymol12/src/mengine/src/etorsion.c
Revision: 99
Committed: Mon Jan 19 04:31:29 2009 UTC (10 years, 10 months ago) by wdelano
File size: 40311 byte(s)
Log Message:
synchronized with trunk, less openmp lib/include
Line File contents
1 #define EXTERN extern
2
3 #include "pcwin.h"
4 #include "pcmod.h"
5
6 #include "energies.h"
7 #include "torsions.h"
8 #include "derivs.h"
9 #include "hess.h"
10 #include "pot.h"
11
12 double dihdrl(int,int,int,int);
13 void etorsion(void);
14 void etorsion1(void);
15 void etorsion2(int);
16
17
18 EXTERN struct t_minim_values {
19 int iprint, ndc, nconst;
20 float dielc;
21 } minim_values;
22
23 void etorsion()
24 {
25 int i,ia, ib, ic, id;
26 double e, rt2, ru2, rtru;
27 double xt,yt,zt,xu,yu,zu;
28 double v1,v2,v3,v4,v5,v6;
29 double s1,s2,s3,s4,s5,s6;
30 double cosine,cosine2,cosine3;
31 double cosine4,cosine5,cosine6;
32 double phi1,phi2,phi3;
33 double phi4,phi5,phi6;
34 double xia,yia,zia,xib,yib,zib;
35 double xic,yic,zic,xid,yid,zid;
36 double xba,yba,zba,xcb,ycb,zcb;
37 double xdc,ydc,zdc;
38 double xtu,ytu,ztu,rcb, angle,rdihed;
39 double width,sine,dt2;
40 int curang;
41
42 energies.etor = 0.0;
43 if (minim_values.iprint)
44 {
45 fprintf(pcmlogfile,"\nTorsion Terms\n");
46 fprintf(pcmlogfile," At1 At2 At3 At4 Types Angle V1 V2 V3 Etor\n");
47 }
48
49 width = 0.0;
50
51 for (i=0; i < torsions.ntor; i++)
52 {
53 ia = torsions.i14[i][0];
54 ib = torsions.i14[i][1];
55 ic = torsions.i14[i][2];
56 id = torsions.i14[i][3];
57 if (atom[ia].use || atom[ib].use || atom[ic].use || atom[id].use )
58 {
59 xia = atom[ia].x;
60 yia = atom[ia].y;
61 zia = atom[ia].z;
62 xib = atom[ib].x;
63 yib = atom[ib].y;
64 zib = atom[ib].z;
65 xic = atom[ic].x;
66 yic = atom[ic].y;
67 zic = atom[ic].z;
68 xid = atom[id].x;
69 yid = atom[id].y;
70 zid = atom[id].z;
71 xba = xib - xia;
72 yba = yib - yia;
73 zba = zib - zia;
74 xcb = xic - xib;
75 ycb = yic - yib;
76 zcb = zic - zib;
77 xdc = xid - xic;
78 ydc = yid - yic;
79 zdc = zid - zic;
80 xt = yba*zcb - ycb*zba;
81 yt = zba*xcb - zcb*xba;
82 zt = xba*ycb - xcb*yba;
83 xu = ycb*zdc - ydc*zcb;
84 yu = zcb*xdc - zdc*xcb;
85 zu = xcb*ydc - xdc*ycb;
86 rt2 = xt*xt + yt*yt + zt*zt;
87 ru2 = xu*xu + yu*yu + zu*zu;
88 rtru = sqrt(rt2 * ru2);
89 if (rtru != 0.0)
90 {
91 cosine = (xt*xu + yt*yu + zt*zu) / rtru;
92
93 v1 = torsions.v1[i];
94 s1 = torsions.ph1[i];
95 v2 = torsions.v2[i];
96 s2 = torsions.ph2[i];
97 v3 = torsions.v3[i];
98 s3 = torsions.ph3[i];
99 v4 = torsions.v4[i];
100 s4 = torsions.ph4[i];
101 v5 = torsions.v5[i];
102 s5 = torsions.ph5[i];
103 v6 = torsions.v6[i];
104 s6 = torsions.ph6[i];
105
106 // compute the powers of the cosine of this angle
107
108 cosine2 = cosine * cosine;
109 cosine3 = cosine2 * cosine;
110 cosine4 = cosine2 * cosine2;
111 cosine5 = cosine3 * cosine2;
112 cosine6 = cosine3 * cosine3;
113
114 phi1 = 1.0 + s1*cosine;
115 phi2 = 1.0 + s2*(2.0*cosine2 - 1.0);
116 phi3 = 1.0 + s3*(4.0*cosine3 - 3.0*cosine);
117 phi4 = 1.0 + s4*(8.0*(cosine4-cosine2) + 1.0);
118 phi5 = 1.0 + s5*(16.0*cosine5 - 20.0*cosine3 + 5.0*cosine);
119 phi6 = 1.0 + s6*(32.0*cosine6 - 48.0*cosine4 + 18.0*cosine2 - 1.0);
120
121 if (pot.use_deform)
122 {
123 phi1 = phi1*exp(-width);
124 phi2 = phi2*exp(-4.0*width);
125 phi3 = phi3*exp(-9.0*width);
126 phi4 = phi4*exp(-16.0*width);
127 phi5 = phi5*exp(-25.0*width);
128 phi6 = phi6*exp(-36.0*width);
129 }
130 e = units.torsunit * (v1*phi1 + v2*phi2 + v3*phi3
131 + v4*phi4 + v5*phi5 + v6*phi6);
132 energies.etor += e;
133 atom[ia].energy += e;
134 atom[ib].energy += e;
135 atom[ic].energy += e;
136 atom[id].energy += e;
137 if(minim_values.iprint)
138 {
139 rdihed = dihdrl(ia,ib,ic,id);
140 fprintf(pcmlogfile,"Tor: %2s(%-3d)- %2s(%-3d)- %2s(%-3d)- %2s(%-3d) %d %d %d %d %8.2f %-8.3f %-8.3f %-8.3f = %-8.4f\n",
141 atom[ia].name,ia,atom[ib].name,ib,atom[ic].name,ic,atom[id].name,id, atom[ia].type, atom[ib].type,
142 atom[ic].type, atom[id].type, rdihed,v1,v2,v3,e);
143 }
144 }
145 }
146 }
147 }
148 // ===========================================================
149 void etorsion1()
150 {
151 int i,ia, ib, ic, id;
152 int curang;
153 double angle;
154 double e, rt2, ru2, rtru,dedphi,sine,dt2;
155 double xt,yt,zt,xu,yu,zu;
156 double v1,v2,v3,v4,v5,v6;
157 double s1,s2,s3,s4,s5,s6;
158 double cosine,cosine2,cosine3;
159 double cosine4,cosine5,cosine6;
160 double phi1,phi2,phi3;
161 double phi4,phi5,phi6;
162 double xia,yia,zia,xib,yib,zib;
163 double xic,yic,zic,xid,yid,zid;
164 double xba,yba,zba,xcb,ycb,zcb;
165 double xdc,ydc,zdc;
166 double xca,yca,zca,xdb,ydb,zdb;
167 double xtu,ytu,ztu,rcb;
168 double dphi1,dphi2,dphi3;
169 double dphi4,dphi5,dphi6;
170 double dphidxt,dphidyt,dphidzt;
171 double dphidxu,dphidyu,dphidzu;
172 double dphidxia,dphidyia,dphidzia;
173 double dphidxib,dphidyib,dphidzib;
174 double dphidxic,dphidyic,dphidzic;
175 double dphidxid,dphidyid,dphidzid;
176 double dedxia,dedyia,dedzia;
177 double dedxib,dedyib,dedzib;
178 double dedxic,dedyic,dedzic;
179 double dedxid,dedyid,dedzid;
180 double dedxt,dedyt,dedzt, dedxu,dedyu,dedzu;
181 double deform1,deform2,deform3;
182 double deform4,deform5,deform6, width;
183
184 energies.etor = 0.0;
185 for (i=0; i <= natom; i++)
186 {
187 deriv.detor[i][0] = 0.0;
188 deriv.detor[i][1] = 0.0;
189 deriv.detor[i][2] = 0.0;
190 }
191
192 width = 0.0;
193
194 for (i=0; i < torsions.ntor; i++)
195 {
196 ia = torsions.i14[i][0];
197 ib = torsions.i14[i][1];
198 ic = torsions.i14[i][2];
199 id = torsions.i14[i][3];
200 if (atom[ia].use || atom[ib].use || atom[ic].use || atom[id].use )
201 {
202 xia = atom[ia].x;
203 yia = atom[ia].y;
204 zia = atom[ia].z;
205 xib = atom[ib].x;
206 yib = atom[ib].y;
207 zib = atom[ib].z;
208 xic = atom[ic].x;
209 yic = atom[ic].y;
210 zic = atom[ic].z;
211 xid = atom[id].x;
212 yid = atom[id].y;
213 zid = atom[id].z;
214 xba = xib - xia;
215 yba = yib - yia;
216 zba = zib - zia;
217 xcb = xic - xib;
218 ycb = yic - yib;
219 zcb = zic - zib;
220 xdc = xid - xic;
221 ydc = yid - yic;
222 zdc = zid - zic;
223 xt = yba*zcb - ycb*zba;
224 yt = zba*xcb - zcb*xba;
225 zt = xba*ycb - xcb*yba;
226 xu = ycb*zdc - ydc*zcb;
227 yu = zcb*xdc - zdc*xcb;
228 zu = xcb*ydc - xdc*ycb;
229 xtu = yt*zu - yu*zt;
230 ytu = zt*xu - zu*xt;
231 ztu = xt*yu - xu*yt;
232 rt2 = xt*xt + yt*yt + zt*zt;
233 ru2 = xu*xu + yu*yu + zu*zu;
234 rtru = sqrt(rt2 * ru2);
235 if (rtru != 0.0)
236 {
237 rcb = sqrt(xcb*xcb + ycb*ycb + zcb*zcb);
238 cosine = (xt*xu + yt*yu + zt*zu) / rtru;
239 sine = (xcb*xtu + ycb*ytu + zcb*ztu) / (rcb*rtru);
240
241 v1 = torsions.v1[i];
242 s1 = torsions.ph1[i];
243 v2 = torsions.v2[i];
244 s2 = torsions.ph2[i];
245 v3 = torsions.v3[i];
246 s3 = torsions.ph3[i];
247 v4 = torsions.v4[i];
248 s4 = torsions.ph4[i];
249 v5 = torsions.v5[i];
250 s5 = torsions.ph5[i];
251 v6 = torsions.v6[i];
252 s6 = torsions.ph6[i];
253 // compute the powers of the cosine of this angle
254
255 cosine2 = cosine * cosine;
256 cosine3 = cosine2 * cosine;
257 cosine4 = cosine2 * cosine2;
258 cosine5 = cosine3 * cosine2;
259 cosine6 = cosine3 * cosine3;
260
261 phi1 = 1.0 + s1*cosine;
262 phi2 = 1.0 + s2*(2.0*cosine2 - 1.0);
263 phi3 = 1.0 + s3*(4.0*cosine3 - 3.0*cosine);
264 phi4 = 1.0 + s4*(8.0*(cosine4-cosine2) + 1.0);
265 phi5 = 1.0 + s5*(16.0*cosine5 - 20.0*cosine3 + 5.0*cosine);
266 phi6 = 1.0 + s6*(32.0*cosine6 - 48.0*cosine4 + 18.0*cosine2 - 1.0);
267
268 dphi1 = s1;
269 dphi2 = s2 * (4.0*cosine);
270 dphi3 = s3 * (12.0*cosine2 - 3.0);
271 dphi4 = s4 * (32.0*cosine3 - 16.0*cosine);
272 dphi5 = s5 * (80.0*cosine4 - 60.0*cosine2 + 5.0);
273 dphi6 = s6 * (192.0*(cosine5-cosine3) + 36.0*cosine);
274
275 if (pot.use_deform)
276 {
277 deform1 = exp(-width);
278 deform2 = exp(-4.0*width);
279 deform3 = exp(-9.0*width);
280 deform4 = exp(-16.0*width);
281 deform5 = exp(-25.0*width);
282 deform6 = exp(-36.0*width);
283 phi1 *= deform1;
284 phi2 *= deform2;
285 phi3 *= deform3;
286 phi4 *= deform4;
287 phi5 *= deform5;
288 phi6 *= deform6;
289 dphi1 *= deform1;
290 dphi2 *= deform2;
291 dphi3 *= deform3;
292 dphi4 *= deform4;
293 dphi5 *= deform5;
294 dphi6 *= deform6;
295 }
296
297 e = units.torsunit * (v1*phi1 + v2*phi2 + v3*phi3
298 + v4*phi4 + v5*phi5 + v6*phi6);
299 dedphi = -sine * units.torsunit
300 * (v1*dphi1+v2*dphi2+v3*dphi3+v4*dphi4+v5*dphi5+v6*dphi6);
301
302 xca = xic - xia;
303 yca = yic - yia;
304 zca = zic - zia;
305 xdb = xid - xib;
306 ydb = yid - yib;
307 zdb = zid - zib;
308 dphidxt = (yt*zcb - ycb*zt) / (rt2*rcb);
309 dphidyt = (zt*xcb - zcb*xt) / (rt2*rcb);
310 dphidzt = (xt*ycb - xcb*yt) / (rt2*rcb);
311 dphidxu = -(yu*zcb - ycb*zu) / (ru2*rcb);
312 dphidyu = -(zu*xcb - zcb*xu) / (ru2*rcb);
313 dphidzu = -(xu*ycb - xcb*yu) / (ru2*rcb);
314
315 dphidxia = zcb*dphidyt - ycb*dphidzt;
316 dphidyia = xcb*dphidzt - zcb*dphidxt;
317 dphidzia = ycb*dphidxt - xcb*dphidyt;
318 dphidxib = yca*dphidzt - zca*dphidyt + zdc*dphidyu - ydc*dphidzu;
319 dphidyib = zca*dphidxt - xca*dphidzt + xdc*dphidzu - zdc*dphidxu;
320 dphidzib = xca*dphidyt - yca*dphidxt + ydc*dphidxu - xdc*dphidyu;
321 dphidxic = zba*dphidyt - yba*dphidzt + ydb*dphidzu - zdb*dphidyu;
322 dphidyic = xba*dphidzt - zba*dphidxt + zdb*dphidxu - xdb*dphidzu;
323 dphidzic = yba*dphidxt - xba*dphidyt + xdb*dphidyu - ydb*dphidxu;
324 dphidxid = zcb*dphidyu - ycb*dphidzu;
325 dphidyid = xcb*dphidzu - zcb*dphidxu;
326 dphidzid = ycb*dphidxu - xcb*dphidyu;
327 dedxia = dedphi * dphidxia;
328 dedyia = dedphi * dphidyia;
329 dedzia = dedphi * dphidzia;
330 dedxib = dedphi * dphidxib;
331 dedyib = dedphi * dphidyib;
332 dedzib = dedphi * dphidzib;
333 dedxic = dedphi * dphidxic;
334 dedyic = dedphi * dphidyic;
335 dedzic = dedphi * dphidzic;
336 dedxid = dedphi * dphidxid;
337 dedyid = dedphi * dphidyid;
338 dedzid = dedphi * dphidzid;
339
340 energies.etor += e;
341 deriv.detor[ia][0] += dedxia;
342 deriv.detor[ia][1] += dedyia;
343 deriv.detor[ia][2] += dedzia;
344
345 deriv.detor[ib][0] += dedxib;
346 deriv.detor[ib][1] += dedyib;
347 deriv.detor[ib][2] += dedzib;
348
349 deriv.detor[ic][0] += dedxic;
350 deriv.detor[ic][1] += dedyic;
351 deriv.detor[ic][2] += dedzic;
352
353 deriv.detor[id][0] += dedxid;
354 deriv.detor[id][1] += dedyid;
355 deriv.detor[id][2] += dedzid;
356
357 }
358 }
359 }
360 }
361 // ===========================================================
362 void etorsion2(iatom)
363 {
364 int i,ia,ib,ic,id;
365 int curang;
366 double agl, cosm, sinm;
367 double dot, denom, denom2, denom3;
368 double dedphi,d2edphi2,sine;
369 double term1, term2, term3;
370 double v1,v2,v3,v4,v5,v6;
371 double s1,s2,s3,s4,s5,s6;
372 double cosine,cosine2,cosine3;
373 double cosine4,cosine5,cosine6;
374 double xia,yia,zia,xib,yib,zib;
375 double xic,yic,zic,xid,yid,zid;
376 double xba,yba,zba,xcb,ycb,zcb;
377 double xdc,ydc,zdc;
378 double xca,yca,zca,xdb,ydb,zdb;
379 double xt,yt,zt,xu,yu,zu,xtu,ytu,ztu;
380 double rt2,ru2,rtru,rcb, rt2ru2;
381 double phi1,phi2,phi3;
382 double phi4,phi5,phi6;
383 double dphi1,dphi2,dphi3;
384 double dphi4,dphi5,dphi6;
385 double d2phi1,d2phi2,d2phi3;
386 double d2phi4,d2phi5,d2phi6;
387 double dphidxt,dphidyt,dphidzt;
388 double dphidxu,dphidyu,dphidzu;
389 double dphidxt1,dphidyt1,dphidzt1;
390 double dphidxu1,dphidyu1,dphidzu1;
391 double dphidxia,dphidyia,dphidzia;
392 double dphidxib,dphidyib,dphidzib;
393 double dphidxic,dphidyic,dphidzic;
394 double dphidxid,dphidyid,dphidzid;
395 double xycb2,xzcb2,yzcb2;
396 double rcbxt,rcbyt,rcbzt,rcbt2;
397 double rcbxu,rcbyu,rcbzu,rcbu2;
398 double dphidxibt,dphidyibt,dphidzibt;
399 double dphidxibu,dphidyibu,dphidzibu;
400 double dphidxict,dphidyict,dphidzict;
401 double dphidxicu,dphidyicu,dphidzicu;
402 double dxiaxia,dyiayia,dziazia,dxibxib,dyibyib,dzibzib;
403 double dxicxic,dyicyic,dziczic,dxidxid,dyidyid,dzidzid;
404 double dxiayia,dxiazia,dyiazia,dxibyib,dxibzib,dyibzib;
405 double dxicyic,dxiczic,dyiczic,dxidyid,dxidzid,dyidzid;
406 double dxiaxib,dxiayib,dxiazib,dyiaxib,dyiayib,dyiazib;
407 double dziaxib,dziayib,dziazib,dxiaxic,dxiayic,dxiazic;
408 double dyiaxic,dyiayic,dyiazic,dziaxic,dziayic,dziazic;
409 double dxiaxid,dxiayid,dxiazid,dyiaxid,dyiayid,dyiazid;
410 double dziaxid,dziayid,dziazid,dxibxic,dxibyic,dxibzic;
411 double dyibxic,dyibyic,dyibzic,dzibxic,dzibyic,dzibzic;
412 double dxibxid,dxibyid,dxibzid,dyibxid,dyibyid,dyibzid;
413 double dzibxid,dzibyid,dzibzid,dxicxid,dxicyid,dxiczid;
414 double dyicxid,dyicyid,dyiczid,dzicxid,dzicyid,dziczid;
415 double yuzab, zuydb, ytzba, rtru32, rtru322;
416 double zdcyu, xdczu, ydcxu, ycazt, zcaxt, xcayt;
417 double deform1,deform2,deform3;
418 double deform4,deform5,deform6,width;
419
420 width = 0.0;
421
422 for (i=0; i < torsions.ntor; i++)
423 {
424 ia = torsions.i14[i][0];
425 ib = torsions.i14[i][1];
426 ic = torsions.i14[i][2];
427 id = torsions.i14[i][3];
428 if (ia == iatom || ib == iatom || ic == iatom || id == iatom)
429 {
430 xia = atom[ia].x;
431 yia = atom[ia].y;
432 zia = atom[ia].z;
433 xib = atom[ib].x;
434 yib = atom[ib].y;
435 zib = atom[ib].z;
436 xic = atom[ic].x;
437 yic = atom[ic].y;
438 zic = atom[ic].z;
439 xid = atom[id].x;
440 yid = atom[id].y;
441 zid = atom[id].z;
442 xba = xib - xia;
443 yba = yib - yia;
444 zba = zib - zia;
445 xcb = xic - xib;
446 ycb = yic - yib;
447 zcb = zic - zib;
448 xdc = xid - xic;
449 ydc = yid - yic;
450 zdc = zid - zic;
451 xt = yba*zcb - ycb*zba;
452 yt = zba*xcb - zcb*xba;
453 zt = xba*ycb - xcb*yba;
454 xu = ycb*zdc - ydc*zcb;
455 yu = zcb*xdc - zdc*xcb;
456 zu = xcb*ydc - xdc*ycb;
457 xtu = yt*zu - yu*zt;
458 ytu = zt*xu - zu*xt;
459 ztu = xt*yu - xu*yt;
460 rt2 = xt*xt + yt*yt + zt*zt;
461 ru2 = xu*xu + yu*yu + zu*zu;
462 rtru = sqrt(rt2 * ru2);
463 if (rtru != 0.0)
464 {
465 rcb = sqrt(xcb*xcb + ycb*ycb + zcb*zcb);
466 cosine = (xt*xu + yt*yu + zt*zu) / rtru;
467 sine = (xcb*xtu + ycb*ytu + zcb*ztu) / (rcb*rtru);
468
469 v1 = torsions.v1[i];
470 s1 = torsions.ph1[i];
471 v2 = torsions.v2[i];
472 s2 = torsions.ph2[i];
473 v3 = torsions.v3[i];
474 s3 = torsions.ph3[i];
475 v4 = torsions.v4[i];
476 s4 = torsions.ph4[i];
477 v5 = torsions.v5[i];
478 s5 = torsions.ph5[i];
479 v6 = torsions.v6[i];
480 s6 = torsions.ph6[i];
481
482 // compute the powers of the cosine of this angle
483
484 cosine2 = cosine * cosine;
485 cosine3 = cosine2 * cosine;
486 cosine4 = cosine2 * cosine2;
487 cosine5 = cosine3 * cosine2;
488 cosine6 = cosine3 * cosine3;
489
490 phi1 = 1.0 + s1*cosine;
491 phi2 = 1.0 + s2*(2.0*cosine2 - 1.0);
492 phi3 = 1.0 + s3*(4.0*cosine3 - 3.0*cosine);
493 phi4 = 1.0 + s4*(8.0*(cosine4-cosine2) + 1.0);
494 phi5 = 1.0 + s5*(16.0*cosine5 - 20.0*cosine3 + 5.0*cosine);
495 phi6 = 1.0 + s6*(32.0*cosine6 - 48.0*cosine4 + 18.0*cosine2 - 1.0);
496
497 dphi1 = s1;
498 dphi2 = s2 * (4.0*cosine);
499 dphi3 = s3 * (12.0*cosine2 - 3.0);
500 dphi4 = s4 * (32.0*cosine3 - 16.0*cosine);
501 dphi5 = s5 * (80.0*cosine4 - 60.0*cosine2 + 5.0);
502 dphi6 = s6 * (192.0*(cosine5-cosine3) + 36.0*cosine);
503
504 d2phi1 = -s1 * cosine;
505 d2phi2 = -s2 * (8.0*cosine2 - 4.0);
506 d2phi3 = -s3 * (36.0*cosine3 - 27.0*cosine);
507 d2phi4 = -s4 * (128.0*(cosine4-cosine2) + 16.0);
508 d2phi5 = -s5 * (400.0*cosine5 - 500.0*cosine + 125.0*cosine);
509 d2phi6 = -s6 * (1152.0*cosine6 - 1728.0*cosine2 + 648.0*cosine2);
510
511 if (pot.use_deform)
512 {
513 deform1 = exp(-width);
514 deform2 = exp(-4.0*width);
515 deform3 = exp(-9.0*width);
516 deform4 = exp(-16.0*width);
517 deform5 = exp(-25.0*width);
518 deform6 = exp(-36.0*width);
519 phi1 *= deform1;
520 phi2 *= deform2;
521 phi3 *= deform3;
522 phi4 *= deform4;
523 phi5 *= deform5;
524 phi6 *= deform6;
525 dphi1 *= deform1;
526 dphi2 *= deform2;
527 dphi3 *= deform3;
528 dphi4 *= deform4;
529 dphi5 *= deform5;
530 dphi6 *= deform6;
531 d2phi1 *= deform1;
532 d2phi2 *= deform2;
533 d2phi3 *= deform3;
534 d2phi4 *= deform4;
535 d2phi5 *= deform5;
536 d2phi6 *= deform6;
537 }
538
539
540 dedphi = -sine * units.torsunit * (v1*dphi1 + v2*dphi2 + v3*dphi3
541 + v4*dphi4 + v5*dphi5 + v6*dphi6);
542 d2edphi2 = units.torsunit * (v1*d2phi1 + v2*d2phi2 + v3*d2phi3
543 + v4*d2phi4 + v5*d2phi5 + v6*d2phi6);
544
545 xca = xic - xia;
546 yca = yic - yia;
547 zca = zic - zia;
548 xdb = xid - xib;
549 ydb = yid - yib;
550 zdb = zid - zib;
551 dphidxt = (yt*zcb - ycb*zt) / (rt2*rcb);
552 dphidyt = (zt*xcb - zcb*xt) / (rt2*rcb);
553 dphidzt = (xt*ycb - xcb*yt) / (rt2*rcb);
554 dphidxu = -(yu*zcb - ycb*zu) / (ru2*rcb);
555 dphidyu = -(zu*xcb - zcb*xu) / (ru2*rcb);
556 dphidzu = -(xu*ycb - xcb*yu) / (ru2*rcb);
557
558 xycb2 = xcb*xcb + ycb*ycb;
559 xzcb2 = xcb*xcb + zcb*zcb;
560 yzcb2 = ycb*ycb + zcb*zcb;
561 rcbxt = -2.0 * rcb * dphidxt;
562 rcbyt = -2.0 * rcb * dphidyt;
563 rcbzt = -2.0 * rcb * dphidzt;
564 rcbt2 = rcb * rt2;
565 rcbxu = 2.0 * rcb * dphidxu;
566 rcbyu = 2.0 * rcb * dphidyu;
567 rcbzu = 2.0 * rcb * dphidzu;
568 rcbu2 = rcb * ru2;
569 dphidxibt = yca*dphidzt - zca*dphidyt;
570 dphidxibu = zdc*dphidyu - ydc*dphidzu;
571 dphidyibt = zca*dphidxt - xca*dphidzt;
572 dphidyibu = xdc*dphidzu - zdc*dphidxu;
573 dphidzibt = xca*dphidyt - yca*dphidxt;
574 dphidzibu = ydc*dphidxu - xdc*dphidyu;
575 dphidxict = zba*dphidyt - yba*dphidzt;
576 dphidxicu = ydb*dphidzu - zdb*dphidyu;
577 dphidyict = xba*dphidzt - zba*dphidxt;
578 dphidyicu = zdb*dphidxu - xdb*dphidzu;
579 dphidzict = yba*dphidxt - xba*dphidyt;
580 dphidzicu = xdb*dphidyu - ydb*dphidxu;
581
582 dphidxia = zcb*dphidyt - ycb*dphidzt;
583 dphidyia = xcb*dphidzt - zcb*dphidxt;
584 dphidzia = ycb*dphidxt - xcb*dphidyt;
585 dphidxib = dphidxibt + dphidxibu;
586 dphidyib = dphidyibt + dphidyibu;
587 dphidzib = dphidzibt + dphidzibu;
588 dphidxic = dphidxict + dphidxicu;
589 dphidyic = dphidyict + dphidyicu;
590 dphidzic = dphidzict + dphidzicu;
591 dphidxid = zcb*dphidyu - ycb*dphidzu;
592 dphidyid = xcb*dphidzu - zcb*dphidxu;
593 dphidzid = ycb*dphidxu - xcb*dphidyu;
594
595 dxiaxia = rcbxt*dphidxia;
596 dxiayia = rcbxt*dphidyia - zcb*rcb/rt2;
597 dxiazia = rcbxt*dphidzia + ycb*rcb/rt2;
598 dxiaxib = rcbxt*dphidxibt + xcb*(zca*ycb-yca*zcb)/rcbt2;
599 dxiayib = rcbxt*dphidyibt + dphidzt + (xca*zcb*xcb+zca*yzcb2)/rcbt2;
600 dxiazib = rcbxt*dphidzibt - dphidyt - (xca*ycb*xcb+yca*yzcb2)/rcbt2;
601 dxiaxic = rcbxt*dphidxict + xcb*xt/rcbt2;
602 dxiayic = rcbxt*dphidyict - dphidzt - (xba*zcb*xcb+zba*yzcb2)/rcbt2;
603 dxiazic = rcbxt*dphidzict + dphidyt + (xba*ycb*xcb+yba*yzcb2)/rcbt2;
604 dxiaxid = 0.0;
605 dxiayid = 0.0;
606 dxiazid = 0.0;
607 dyiayia = rcbyt*dphidyia;
608 dyiazia = rcbyt*dphidzia - xcb*rcb/rt2;
609 dyiaxib = rcbyt*dphidxibt - dphidzt - (yca*zcb*ycb+zca*xzcb2)/rcbt2;
610 dyiayib = rcbyt*dphidyibt + ycb*(xca*zcb-zca*xcb)/rcbt2;
611 dyiazib = rcbyt*dphidzibt + dphidxt + (yca*xcb*ycb+xca*xzcb2)/rcbt2;
612 dyiaxic = rcbyt*dphidxict + dphidzt + (yba*zcb*ycb+zba*xzcb2)/rcbt2;
613 dyiayic = rcbyt*dphidyict + ycb*yt/rcbt2;
614 dyiazic = rcbyt*dphidzict - dphidxt - (yba*xcb*ycb+xba*xzcb2)/rcbt2;
615 dyiaxid = 0.0;
616 dyiayid = 0.0;
617 dyiazid = 0.0;
618 dziazia = rcbzt*dphidzia;
619 dziaxib = rcbzt*dphidxibt + dphidyt + (zca*ycb*zcb+yca*xycb2)/rcbt2;
620 dziayib = rcbzt*dphidyibt - dphidxt - (zca*xcb*zcb+xca*xycb2)/rcbt2;
621 dziazib = rcbzt*dphidzibt + zcb*(yca*xcb-xca*ycb)/rcbt2;
622 dziaxic = rcbzt*dphidxict - dphidyt - (zba*ycb*zcb+yba*xycb2)/rcbt2;
623 dziayic = rcbzt*dphidyict + dphidxt + (zba*xcb*zcb+xba*xycb2)/rcbt2;
624 dziazic = rcbzt*dphidzict + zcb*zt/rcbt2;
625 dziaxid = 0.0;
626 dziayid = 0.0;
627 dziazid = 0.0;
628 dxibxic = -xcb*dphidxib/(rcb*rcb) - (yca*(zba*xcb+yt)-zca*(yba*xcb-zt))/rcbt2
629 - 2.0*(yt*zba-yba*zt)*dphidxibt/rt2 - (zdc*(ydb*xcb+zu)-ydc*(zdb*xcb-yu))/rcbu2
630 + 2.0*(yu*zdb-ydb*zu)*dphidxibu/ru2;
631 dxibyic = -ycb*dphidxib/(rcb*rcb) + dphidzt + dphidzu - (yca*(zba*ycb-xt)+zca*(xba*xcb+zcb*zba))/rcbt2
632 - 2.0*(zt*xba-zba*xt)*dphidxibt/rt2 + (zdc*(xdb*xcb+zcb*zdb)+ydc*(zdb*ycb+xu))/rcbu2
633 + 2.0*(zu*xdb-zdb*xu)*dphidxibu/ru2;
634 dxibxid = rcbxu*dphidxibu + xcb*xu/rcbu2;
635 dxibyid = rcbyu*dphidxibu - dphidzu - (ydc*zcb*ycb+zdc*xzcb2)/rcbu2;
636 dxibzid = rcbzu*dphidxibu + dphidyu + (zdc*ycb*zcb+ydc*xycb2)/rcbu2;
637 dyibzib = ycb*dphidzib/(rcb*rcb) - (xca*(xca*xcb+zcb*zca)+yca*(ycb*xca+zt))/rcbt2
638 - 2.0*(xt*zca-xca*zt)*dphidzibt/rt2 + (ydc*(xdc*ycb-zu)+xdc*(xdc*xcb+zcb*zdc))/rcbu2
639 + 2.0*(xu*zdc-xdc*zu)*dphidzibu/ru2;
640 dyibxic = -xcb*dphidyib/(rcb*rcb) - dphidzt - dphidzu + (xca*(zba*xcb+yt)+zca*(zba*zcb+ycb*yba))/rcbt2
641 - 2.0*(yt*zba-yba*zt)*dphidyibt/rt2 - (zdc*(zdb*zcb+ycb*ydb)+xdc*(zdb*xcb-yu))/rcbu2
642 + 2.0*(yu*zdb-ydb*zu)*dphidyibu/ru2;
643 dyibyic = -ycb*dphidyib/(rcb*rcb)- (zca*(xba*ycb+zt)-xca*(zba*ycb-xt))/rcbt2
644 - 2.0*(zt*xba-zba*xt)*dphidyibt/rt2 - (xdc*(zdb*ycb+xu)-zdc*(xdb*ycb-zu))/rcbu2
645 + 2.0*(zu*xdb-zdb*xu)*dphidyibu/ru2;
646 dyibxid = rcbxu*dphidyibu + dphidzu + (xdc*zcb*xcb+zdc*yzcb2)/rcbu2;
647 dyibyid = rcbyu*dphidyibu + ycb*yu/rcbu2;
648 dyibzid = rcbzu*dphidyibu - dphidxu - (zdc*xcb*zcb+xdc*xycb2)/rcbu2;
649 dzibxic = -xcb*dphidzib/(rcb*rcb) + dphidyt + dphidyu - (xca*(yba*xcb-zt)+yca*(zba*zcb+ycb*yba))/rcbt2
650 - 2.0*(yt*zba-yba*zt)*dphidzibt/rt2 + (ydc*(zdb*zcb+ycb*ydb)+xdc*(ydb*xcb+zu))/rcbu2
651 + 2.0*(yu*zdb-ydb*zu)*dphidzibu/ru2;
652 dzibzic = -zcb*dphidzib/(rcb*rcb) - (xca*(yba*zcb+xt)-yca*(xba*zcb-yt))/rcbt2
653 - 2.0*(xt*yba-xba*yt)*dphidzibt/rt2 - (ydc*(xdb*zcb+yu)-xdc*(ydb*zcb-xu))/rcbu2
654 + 2.0*(xu*ydb-xdb*yu)*dphidzibu/ru2;
655 dzibxid = rcbxu*dphidzibu - dphidyu - (xdc*ycb*xcb+ydc*yzcb2)/rcbu2;
656 dzibyid = rcbyu*dphidzibu + dphidxu + (ydc*xcb*ycb+xdc*xzcb2)/rcbu2;
657 dzibzid = rcbzu*dphidzibu + zcb*zu/rcbu2;
658 dxicxid = rcbxu*dphidxicu - xcb*(zdb*ycb-ydb*zcb)/rcbu2;
659 dxicyid = rcbyu*dphidxicu + dphidzu + (ydb*zcb*ycb+zdb*xzcb2)/rcbu2;
660 dxiczid = rcbzu*dphidxicu - dphidyu - (zdb*ycb*zcb+ydb*xycb2)/rcbu2;
661 dyicxid = rcbxu*dphidyicu - dphidzu - (xdb*zcb*xcb+zdb*yzcb2)/rcbu2;
662 dyicyid = rcbyu*dphidyicu - ycb*(xdb*zcb-zdb*xcb)/rcbu2;
663 dyiczid = rcbzu*dphidyicu + dphidxu + (zdb*xcb*zcb+xdb*xycb2)/rcbu2;
664 dzicxid = rcbxu*dphidzicu + dphidyu + (xdb*ycb*xcb+ydb*yzcb2)/rcbu2;
665 dzicyid = rcbyu*dphidzicu - dphidxu - (ydb*xcb*ycb+xdb*xzcb2)/rcbu2;
666 dziczid = rcbzu*dphidzicu - zcb*(ydb*xcb-xdb*ycb)/rcbu2;
667 dxidxid = rcbxu*dphidxid;
668 dxidyid = rcbxu*dphidyid + zcb*rcb/ru2;
669 dxidzid = rcbxu*dphidzid - ycb*rcb/ru2;
670 dyidyid = rcbyu*dphidyid;
671 dyidzid = rcbyu*dphidzid + xcb*rcb/ru2;
672 dzidzid = rcbzu*dphidzid;
673
674 dxibxib = -dxiaxib - dxibxic - dxibxid;
675 dxibyib = -dyiaxib - dxibyic - dxibyid;
676 dxibzib = -dxiazib - dzibxic - dzibxid;
677 dxibzic = -dziaxib - dxibzib - dxibzid;
678 dyibyib = -dyiayib - dyibyic - dyibyid;
679 dyibzic = -dziayib - dyibzib - dyibzid;
680 dzibzib = -dziazib - dzibzic - dzibzid;
681 dzibyic = -dyiazib - dyibzib - dzibyid;
682 dxicxic = -dxiaxic - dxibxic - dxicxid;
683 dxicyic = -dyiaxic - dyibxic - dxicyid;
684 dxiczic = -dziaxic - dzibxic - dxiczid;
685 dyicyic = -dyiayic - dyibyic - dyicyid;
686 dyiczic = -dziayic - dzibyic - dyiczid;
687 dziczic = -dziazic - dzibzic - dziczid;
688
689 if (iatom == ia)
690 {
691 hess.hessx[ia][0] += dedphi*dxiaxia + d2edphi2*dphidxia*dphidxia;
692 hess.hessy[ia][0] += dedphi*dxiayia + d2edphi2*dphidxia*dphidyia;
693 hess.hessz[ia][0] += dedphi*dxiazia + d2edphi2*dphidxia*dphidzia;
694 hess.hessx[ia][1] += dedphi*dxiayia + d2edphi2*dphidxia*dphidyia;
695 hess.hessy[ia][1] += dedphi*dyiayia + d2edphi2*dphidyia*dphidyia;
696 hess.hessz[ia][1] += dedphi*dyiazia + d2edphi2*dphidyia*dphidzia;
697 hess.hessx[ia][2] += dedphi*dxiazia + d2edphi2*dphidxia*dphidzia;
698 hess.hessy[ia][2] += dedphi*dyiazia + d2edphi2*dphidyia*dphidzia;
699 hess.hessz[ia][2] += dedphi*dziazia + d2edphi2*dphidzia*dphidzia;
700 hess.hessx[ib][0] += dedphi*dxiaxib + d2edphi2*dphidxia*dphidxib;
701 hess.hessy[ib][0] += dedphi*dyiaxib + d2edphi2*dphidyia*dphidxib;
702 hess.hessz[ib][0] += dedphi*dziaxib + d2edphi2*dphidzia*dphidxib;
703 hess.hessx[ib][1] += dedphi*dxiayib + d2edphi2*dphidxia*dphidyib;
704 hess.hessy[ib][1] += dedphi*dyiayib + d2edphi2*dphidyia*dphidyib;
705 hess.hessz[ib][1] += dedphi*dziayib + d2edphi2*dphidzia*dphidyib;
706 hess.hessx[ib][2] += dedphi*dxiazib + d2edphi2*dphidxia*dphidzib;
707 hess.hessy[ib][2] += dedphi*dyiazib + d2edphi2*dphidyia*dphidzib;
708 hess.hessz[ib][2] += dedphi*dziazib + d2edphi2*dphidzia*dphidzib;
709 hess.hessx[ic][0] += dedphi*dxiaxic + d2edphi2*dphidxia*dphidxic;
710 hess.hessy[ic][0] += dedphi*dyiaxic + d2edphi2*dphidyia*dphidxic;
711 hess.hessz[ic][0] += dedphi*dziaxic + d2edphi2*dphidzia*dphidxic;
712 hess.hessx[ic][1] += dedphi*dxiayic + d2edphi2*dphidxia*dphidyic;
713 hess.hessy[ic][1] += dedphi*dyiayic + d2edphi2*dphidyia*dphidyic;
714 hess.hessz[ic][1] += dedphi*dziayic + d2edphi2*dphidzia*dphidyic;
715 hess.hessx[ic][2] += dedphi*dxiazic + d2edphi2*dphidxia*dphidzic;
716 hess.hessy[ic][2] += dedphi*dyiazic + d2edphi2*dphidyia*dphidzic;
717 hess.hessz[ic][2] += dedphi*dziazic + d2edphi2*dphidzia*dphidzic;
718 hess.hessx[id][0] += dedphi*dxiaxid + d2edphi2*dphidxia*dphidxid;
719 hess.hessy[id][0] += dedphi*dyiaxid + d2edphi2*dphidyia*dphidxid;
720 hess.hessz[id][0] += dedphi*dziaxid + d2edphi2*dphidzia*dphidxid;
721 hess.hessx[id][1] += dedphi*dxiayid + d2edphi2*dphidxia*dphidyid;
722 hess.hessy[id][1] += dedphi*dyiayid + d2edphi2*dphidyia*dphidyid;
723 hess.hessz[id][1] += dedphi*dziayid + d2edphi2*dphidzia*dphidyid;
724 hess.hessx[id][2] += dedphi*dxiazid + d2edphi2*dphidxia*dphidzid;
725 hess.hessy[id][2] += dedphi*dyiazid + d2edphi2*dphidyia*dphidzid;
726 hess.hessz[id][2] += dedphi*dziazid + d2edphi2*dphidzia*dphidzid;
727 }else if (iatom == ib)
728 {
729 hess.hessx[ib][0] += dedphi*dxibxib + d2edphi2*dphidxib*dphidxib;
730 hess.hessy[ib][0] += dedphi*dxibyib + d2edphi2*dphidxib*dphidyib;
731 hess.hessz[ib][0] += dedphi*dxibzib + d2edphi2*dphidxib*dphidzib;
732 hess.hessx[ib][1] += dedphi*dxibyib + d2edphi2*dphidxib*dphidyib;
733 hess.hessy[ib][1] += dedphi*dyibyib + d2edphi2*dphidyib*dphidyib;
734 hess.hessz[ib][1] += dedphi*dyibzib + d2edphi2*dphidyib*dphidzib;
735 hess.hessx[ib][2] += dedphi*dxibzib + d2edphi2*dphidxib*dphidzib;
736 hess.hessy[ib][2] += dedphi*dyibzib + d2edphi2*dphidyib*dphidzib;
737 hess.hessz[ib][2] += dedphi*dzibzib + d2edphi2*dphidzib*dphidzib;
738 hess.hessx[ia][0] += dedphi*dxiaxib + d2edphi2*dphidxib*dphidxia;
739 hess.hessy[ia][0] += dedphi*dxiayib + d2edphi2*dphidyib*dphidxia;
740 hess.hessz[ia][0] += dedphi*dxiazib + d2edphi2*dphidzib*dphidxia;
741 hess.hessx[ia][1] += dedphi*dyiaxib + d2edphi2*dphidxib*dphidyia;
742 hess.hessy[ia][1] += dedphi*dyiayib + d2edphi2*dphidyib*dphidyia;
743 hess.hessz[ia][1] += dedphi*dyiazib + d2edphi2*dphidzib*dphidyia;
744 hess.hessx[ia][2] += dedphi*dziaxib + d2edphi2*dphidxib*dphidzia;
745 hess.hessy[ia][2] += dedphi*dziayib + d2edphi2*dphidyib*dphidzia;
746 hess.hessz[ia][2] += dedphi*dziazib + d2edphi2*dphidzib*dphidzia;
747 hess.hessx[ic][0] += dedphi*dxibxic + d2edphi2*dphidxib*dphidxic;
748 hess.hessy[ic][0] += dedphi*dyibxic + d2edphi2*dphidyib*dphidxic;
749 hess.hessz[ic][0] += dedphi*dzibxic + d2edphi2*dphidzib*dphidxic;
750 hess.hessx[ic][1] += dedphi*dxibyic + d2edphi2*dphidxib*dphidyic;
751 hess.hessy[ic][1] += dedphi*dyibyic + d2edphi2*dphidyib*dphidyic;
752 hess.hessz[ic][1] += dedphi*dzibyic + d2edphi2*dphidzib*dphidyic;
753 hess.hessx[ic][2] += dedphi*dxibzic + d2edphi2*dphidxib*dphidzic;
754 hess.hessy[ic][2] += dedphi*dyibzic + d2edphi2*dphidyib*dphidzic;
755 hess.hessz[ic][2] += dedphi*dzibzic + d2edphi2*dphidzib*dphidzic;
756 hess.hessx[id][0] += dedphi*dxibxid + d2edphi2*dphidxib*dphidxid;
757 hess.hessy[id][0] += dedphi*dyibxid + d2edphi2*dphidyib*dphidxid;
758 hess.hessz[id][0] += dedphi*dzibxid + d2edphi2*dphidzib*dphidxid;
759 hess.hessx[id][1] += dedphi*dxibyid + d2edphi2*dphidxib*dphidyid;
760 hess.hessy[id][1] += dedphi*dyibyid + d2edphi2*dphidyib*dphidyid;
761 hess.hessz[id][1] += dedphi*dzibyid + d2edphi2*dphidzib*dphidyid;
762 hess.hessx[id][2] += dedphi*dxibzid + d2edphi2*dphidxib*dphidzid;
763 hess.hessy[id][2] += dedphi*dyibzid + d2edphi2*dphidyib*dphidzid;
764 hess.hessz[id][2] += dedphi*dzibzid + d2edphi2*dphidzib*dphidzid;
765 }else if (iatom == ic)
766 {
767 hess.hessx[ic][0] += dedphi*dxicxic + d2edphi2*dphidxic*dphidxic;
768 hess.hessy[ic][0] += dedphi*dxicyic + d2edphi2*dphidxic*dphidyic;
769 hess.hessz[ic][0] += dedphi*dxiczic + d2edphi2*dphidxic*dphidzic;
770 hess.hessx[ic][1] += dedphi*dxicyic + d2edphi2*dphidxic*dphidyic;
771 hess.hessy[ic][1] += dedphi*dyicyic + d2edphi2*dphidyic*dphidyic;
772 hess.hessz[ic][1] += dedphi*dyiczic + d2edphi2*dphidyic*dphidzic;
773 hess.hessx[ic][2] += dedphi*dxiczic + d2edphi2*dphidxic*dphidzic;
774 hess.hessy[ic][2] += dedphi*dyiczic + d2edphi2*dphidyic*dphidzic;
775 hess.hessz[ic][2] += dedphi*dziczic + d2edphi2*dphidzic*dphidzic;
776 hess.hessx[ia][0] += dedphi*dxiaxic + d2edphi2*dphidxic*dphidxia;
777 hess.hessy[ia][0] += dedphi*dxiayic + d2edphi2*dphidyic*dphidxia;
778 hess.hessz[ia][0] += dedphi*dxiazic + d2edphi2*dphidzic*dphidxia;
779 hess.hessx[ia][1] += dedphi*dyiaxic + d2edphi2*dphidxic*dphidyia;
780 hess.hessy[ia][1] += dedphi*dyiayic + d2edphi2*dphidyic*dphidyia;
781 hess.hessz[ia][1] += dedphi*dyiazic + d2edphi2*dphidzic*dphidyia;
782 hess.hessx[ia][2] += dedphi*dziaxic + d2edphi2*dphidxic*dphidzia;
783 hess.hessy[ia][2] += dedphi*dziayic + d2edphi2*dphidyic*dphidzia;
784 hess.hessz[ia][2] += dedphi*dziazic + d2edphi2*dphidzic*dphidzia;
785 hess.hessx[ib][0] += dedphi*dxibxic + d2edphi2*dphidxic*dphidxib;
786 hess.hessy[ib][0] += dedphi*dxibyic + d2edphi2*dphidyic*dphidxib;
787 hess.hessz[ib][0] += dedphi*dxibzic + d2edphi2*dphidzic*dphidxib;
788 hess.hessx[ib][1] += dedphi*dyibxic + d2edphi2*dphidxic*dphidyib;
789 hess.hessy[ib][1] += dedphi*dyibyic + d2edphi2*dphidyic*dphidyib;
790 hess.hessz[ib][1] += dedphi*dyibzic + d2edphi2*dphidzic*dphidyib;
791 hess.hessx[ib][2] += dedphi*dzibxic + d2edphi2*dphidxic*dphidzib;
792 hess.hessy[ib][2] += dedphi*dzibyic + d2edphi2*dphidyic*dphidzib;
793 hess.hessz[ib][2] += dedphi*dzibzic + d2edphi2*dphidzic*dphidzib;
794 hess.hessx[id][0] += dedphi*dxicxid + d2edphi2*dphidxic*dphidxid;
795 hess.hessy[id][0] += dedphi*dyicxid + d2edphi2*dphidyic*dphidxid;
796 hess.hessz[id][0] += dedphi*dzicxid + d2edphi2*dphidzic*dphidxid;
797 hess.hessx[id][1] += dedphi*dxicyid + d2edphi2*dphidxic*dphidyid;
798 hess.hessy[id][1] += dedphi*dyicyid + d2edphi2*dphidyic*dphidyid;
799 hess.hessz[id][1] += dedphi*dzicyid + d2edphi2*dphidzic*dphidyid;
800 hess.hessx[id][2] += dedphi*dxiczid + d2edphi2*dphidxic*dphidzid;
801 hess.hessy[id][2] += dedphi*dyiczid + d2edphi2*dphidyic*dphidzid;
802 hess.hessz[id][2] += dedphi*dziczid + d2edphi2*dphidzic*dphidzid;
803 }else if (iatom == id)
804 {
805 hess.hessx[id][0] += dedphi*dxidxid + d2edphi2*dphidxid*dphidxid;
806 hess.hessy[id][0] += dedphi*dxidyid + d2edphi2*dphidxid*dphidyid;
807 hess.hessz[id][0] += dedphi*dxidzid + d2edphi2*dphidxid*dphidzid;
808 hess.hessx[id][1] += dedphi*dxidyid + d2edphi2*dphidxid*dphidyid;
809 hess.hessy[id][1] += dedphi*dyidyid + d2edphi2*dphidyid*dphidyid;
810 hess.hessz[id][1] += dedphi*dyidzid + d2edphi2*dphidyid*dphidzid;
811 hess.hessx[id][2] += dedphi*dxidzid + d2edphi2*dphidxid*dphidzid;
812 hess.hessy[id][2] += dedphi*dyidzid + d2edphi2*dphidyid*dphidzid;
813 hess.hessz[id][2] += dedphi*dzidzid + d2edphi2*dphidzid*dphidzid;
814 hess.hessx[ia][0] += dedphi*dxiaxid + d2edphi2*dphidxid*dphidxia;
815 hess.hessy[ia][0] += dedphi*dxiayid + d2edphi2*dphidyid*dphidxia;
816 hess.hessz[ia][0] += dedphi*dxiazid + d2edphi2*dphidzid*dphidxia;
817 hess.hessx[ia][1] += dedphi*dyiaxid + d2edphi2*dphidxid*dphidyia;
818 hess.hessy[ia][1] += dedphi*dyiayid + d2edphi2*dphidyid*dphidyia;
819 hess.hessz[ia][1] += dedphi*dyiazid + d2edphi2*dphidzid*dphidyia;
820 hess.hessx[ia][2] += dedphi*dziaxid + d2edphi2*dphidxid*dphidzia;
821 hess.hessy[ia][2] += dedphi*dziayid + d2edphi2*dphidyid*dphidzia;
822 hess.hessz[ia][2] += dedphi*dziazid + d2edphi2*dphidzid*dphidzia;
823 hess.hessx[ib][0] += dedphi*dxibxid + d2edphi2*dphidxid*dphidxib;
824 hess.hessy[ib][0] += dedphi*dxibyid + d2edphi2*dphidyid*dphidxib;
825 hess.hessz[ib][0] += dedphi*dxibzid + d2edphi2*dphidzid*dphidxib;
826 hess.hessx[ib][1] += dedphi*dyibxid + d2edphi2*dphidxid*dphidyib;
827 hess.hessy[ib][1] += dedphi*dyibyid + d2edphi2*dphidyid*dphidyib;
828 hess.hessz[ib][1] += dedphi*dyibzid + d2edphi2*dphidzid*dphidyib;
829 hess.hessx[ib][2] += dedphi*dzibxid + d2edphi2*dphidxid*dphidzib;
830 hess.hessy[ib][2] += dedphi*dzibyid + d2edphi2*dphidyid*dphidzib;
831 hess.hessz[ib][2] += dedphi*dzibzid + d2edphi2*dphidzid*dphidzib;
832 hess.hessx[ic][0] += dedphi*dxicxid + d2edphi2*dphidxid*dphidxic;
833 hess.hessy[ic][0] += dedphi*dxicyid + d2edphi2*dphidyid*dphidxic;
834 hess.hessz[ic][0] += dedphi*dxiczid + d2edphi2*dphidzid*dphidxic;
835 hess.hessx[ic][1] += dedphi*dyicxid + d2edphi2*dphidxid*dphidyic;
836 hess.hessy[ic][1] += dedphi*dyicyid + d2edphi2*dphidyid*dphidyic;
837 hess.hessz[ic][1] += dedphi*dyiczid + d2edphi2*dphidzid*dphidyic;
838 hess.hessx[ic][2] += dedphi*dzicxid + d2edphi2*dphidxid*dphidzic;
839 hess.hessy[ic][2] += dedphi*dzicyid + d2edphi2*dphidyid*dphidzic;
840 hess.hessz[ic][2] += dedphi*dziczid + d2edphi2*dphidzid*dphidzic;
841 }
842 }
843 }
844 }
845 }