// bond torsion semi-harmonic energy: dihedral X=B-A,Y=C-B,Z=D-C.  plane defined by X,Y  and Y,Z
void eh_torsion_hess(void *, double *dVdXX, index_t nb, double const &Factor, double const *pXyz, index_t const *iAtom, double theta0, double k)
{
   index_t iAx = 3.0*iAtom[0] + 0;
   index_t iAy = 3.0*iAtom[0] + 1;
   index_t iAz = 3.0*iAtom[0] + 2;
   index_t iBx = 3.0*iAtom[1] + 0;
   index_t iBy = 3.0*iAtom[1] + 1;
   index_t iBz = 3.0*iAtom[1] + 2;
   index_t iCx = 3.0*iAtom[2] + 0;
   index_t iCy = 3.0*iAtom[2] + 1;
   index_t iCz = 3.0*iAtom[2] + 2;
   index_t iDx = 3.0*iAtom[3] + 0;
   index_t iDy = 3.0*iAtom[3] + 1;
   index_t iDz = 3.0*iAtom[3] + 2;
   double Ax = pXyz[iAx];
   double Ay = pXyz[iAy];
   double Az = pXyz[iAz];
   double Bx = pXyz[iBx];
   double By = pXyz[iBy];
   double Bz = pXyz[iBz];
   double Cx = pXyz[iCx];
   double Cy = pXyz[iCy];
   double Cz = pXyz[iCz];
   double Dx = pXyz[iDx];
   double Dy = pXyz[iDy];
   double Dz = pXyz[iDz];
   double r1 = -1.0*Bx;
   double r2 = Cx + r1;
   double r3 = -1.0*By;
   double r4 = Cy + r3;
   double r5 = -1.0*Bz;
   double r6 = Cz + r5;
   double r7 = sqr(r6) + sqr(r4) + sqr(r2);
   double r8 = sqrt(r7);
   double r9 = -1.0*Cy;
   double r10 = r9 + By;
   double r11 = -1.0*Cx;
   double r12 = Dx + r11;
   double r13 = -1.0*r12*r4;
   double r14 = Dy + r9;
   double r15 = r2*r14;
   double r16 = r15 + r13;
   double r17 = r6*r12;
   double r18 = -1.0*Cz;
   double r19 = Dz + r18;
   double r20 = -1.0*r19*r2;
   double r21 = r20 + r17;
   double r22 = -1.0*Az;
   double r23 = Bz + r22;
   double r24 = -1.0*Ay;
   double r25 = By + r24;
   double r26 = -1.0*Ax;
   double r27 = Bx + r26;
   double r28 = -1.0*r14*r6;
   double r29 = r4*r19;
   double r30 = r29 + r28;
   double r31 = r27*r30 + r25*r21 + r23*r16;
   double r32 = r27*r4;
   double r33 = r32 - 1.0*r2*r25;
   double r34 = r23*r2;
   double r35 = r34 - 1.0*r27*r6;
   double r36 = r25*r6;
   double r37 = r36 - 1.0*r23*r4;
   double r38 = r37*r30 + r35*r21 + r33*r16;
   double r39 = 1.0/(sqr(r38) + r7*sqr(r31));
   double r40 = r38*r39*(r6*r14 - 1.0*r19*r4)*r8 - 1.0*(r6*r21 + r10*r16)*r31*r39*r8;
   double r41 = r18 + Bz;
   double r42 = r38*r39*(r2*r19 - 1.0*r12*r6)*r8 - 1.0*(r41*r30 + r2*r16)*r31*r39*r8;
   double r43 = r42*r40*k;
   double r44 = r11 + Bx;
   double r45 = (r4*r12 - 1.0*r14*r2)*r38*r39*r8 - 1.0*(r4*r30 + r44*r21)*r31*r39*r8;
   double r46 = r45*r40*k;
   double r47 = -1.0*Dy;
   double r48 = r47 + Cy;
   double r49 = 1.0/r8;
   double r50 = r38*r39*(r8*(r29 + r25*r19 + r28 + r23*r48) - 1.0*r2*r31*r49) - 1.0*((r18 + Az)*r21 + r35*r19 + (Cy + r24)*r16 + r33*r48)*r31*r39*r8;
   double r51 = r40*r50*k;
   double r52 = -1.0*Dz;
   double r53 = r52 + Cz;
   double r54 = r38*r39*(r8*(r20 + r27*r53 + r17 + r23*r12) - 1.0*r31*r4*r49) - 1.0*((Cz + r22)*r30 + r37*r53 + (r11 + Ax)*r16 + r33*r12)*r31*r39*r8;
   double r55 = r40*r54*k;
   double r56 = -1.0*Dx;
   double r57 = r56 + Cx;
   double r58 = r38*r39*(r8*(r15 + r27*r14 + r13 + r25*r57) - 1.0*r31*r49*r6) - 1.0*((r9 + Ay)*r30 + (Cx + r26)*r21 + r37*r14 + r35*r57)*r31*r39*r8;
   double r59 = r40*r58*k;
   double r60 = Dy + r3;
   double r61 = r52 + Bz;
   double r62 = r38*(r2*r49*r31 + r8*(r25*r61 + r23*r60))*r39 - 1.0*(r23*r21 + r35*r61 + (r3 + Ay)*r16 + r33*r60)*r31*r39*r8;
   double r63 = r40*r62*k;
   double r64 = r56 + Bx;
   double r65 = Dz + r5;
   double r66 = r38*(r4*r49*r31 + r8*(r27*r65 + r23*r64))*r39 - 1.0*((r5 + Az)*r30 + r37*r65 + r27*r16 + r33*r64)*r31*r39*r8;
   double r67 = r40*r66*k;
   double r68 = Dx + r1;
   double r69 = r47 + By;
   double r70 = r38*(r6*r49*r31 + r8*(r27*r69 + r25*r68))*r39 - 1.0*(r25*r30 + (r1 + Ax)*r21 + r37*r69 + r35*r68)*r31*r39*r8;
   double r71 = r40*r70*k;
   double r72 = (r36 + r23*r10)*r8*r38*r39 - 1.0*r31*(r6*r35 + r10*r33)*r39*r8;
   double r73 = r72*r40*k;
   double r74 = (r27*r41 + r34)*r8*r38*r39 - 1.0*r31*(r41*r37 + r2*r33)*r39*r8;
   double r75 = r74*r40*k;
   double r76 = (r32 + r25*r44)*r8*r38*r39 - 1.0*r31*(r4*r37 + r44*r35)*r39*r8;
   double r77 = r76*r40*k;
   double r78 = r45*r42*k;
   double r79 = r42*r50*k;
   double r80 = r42*r54*k;
   double r81 = r42*r58*k;
   double r82 = r42*r62*k;
   double r83 = r42*r66*k;
   double r84 = r42*r70*k;
   double r85 = r72*r42*k;
   double r86 = r74*r42*k;
   double r87 = r76*r42*k;
   double r88 = r45*r50*k;
   double r89 = r45*r54*k;
   double r90 = r45*r58*k;
   double r91 = r45*r62*k;
   double r92 = r45*r66*k;
   double r93 = r45*r70*k;
   double r94 = r72*r45*k;
   double r95 = r74*r45*k;
   double r96 = r76*r45*k;
   double r97 = r50*r54*k;
   double r98 = r50*r58*k;
   double r99 = r50*r62*k;
   double r100 = r50*r66*k;
   double r101 = r50*r70*k;
   double r102 = r72*r50*k;
   double r103 = r74*r50*k;
   double r104 = r76*r50*k;
   double r105 = r54*r58*k;
   double r106 = r62*r54*k;
   double r107 = r54*r66*k;
   double r108 = r54*r70*k;
   double r109 = r72*r54*k;
   double r110 = r74*r54*k;
   double r111 = r76*r54*k;
   double r112 = r62*r58*k;
   double r113 = r66*r58*k;
   double r114 = r58*r70*k;
   double r115 = r72*r58*k;
   double r116 = r74*r58*k;
   double r117 = r76*r58*k;
   double r118 = r62*r66*k;
   double r119 = r62*r70*k;
   double r120 = r72*r62*k;
   double r121 = r74*r62*k;
   double r122 = r76*r62*k;
   double r123 = r66*r70*k;
   double r124 = r72*r66*k;
   double r125 = r74*r66*k;
   double r126 = r76*r66*k;
   double r127 = r72*r70*k;
   double r128 = r74*r70*k;
   double r129 = r76*r70*k;
   double r130 = r74*r72*k;
   double r131 = r76*r72*k;
   double r132 = r76*r74*k;
   dVdXX[iAx + nb*(iAx)] += Factor*(sqr(r40)*k);
   dVdXX[iAx + nb*(iAy)] += Factor*(r43);
   dVdXX[iAx + nb*(iAz)] += Factor*(r46);
   dVdXX[iAx + nb*(iBx)] += Factor*(r51);
   dVdXX[iAx + nb*(iBy)] += Factor*(r55);
   dVdXX[iAx + nb*(iBz)] += Factor*(r59);
   dVdXX[iAx + nb*(iCx)] += Factor*(r63);
   dVdXX[iAx + nb*(iCy)] += Factor*(r67);
   dVdXX[iAx + nb*(iCz)] += Factor*(r71);
   dVdXX[iAx + nb*(iDx)] += Factor*(r73);
   dVdXX[iAx + nb*(iDy)] += Factor*(r75);
   dVdXX[iAx + nb*(iDz)] += Factor*(r77);
   dVdXX[iAy + nb*(iAx)] += Factor*(r43);
   dVdXX[iAy + nb*(iAy)] += Factor*(sqr(r42)*k);
   dVdXX[iAy + nb*(iAz)] += Factor*(r78);
   dVdXX[iAy + nb*(iBx)] += Factor*(r79);
   dVdXX[iAy + nb*(iBy)] += Factor*(r80);
   dVdXX[iAy + nb*(iBz)] += Factor*(r81);
   dVdXX[iAy + nb*(iCx)] += Factor*(r82);
   dVdXX[iAy + nb*(iCy)] += Factor*(r83);
   dVdXX[iAy + nb*(iCz)] += Factor*(r84);
   dVdXX[iAy + nb*(iDx)] += Factor*(r85);
   dVdXX[iAy + nb*(iDy)] += Factor*(r86);
   dVdXX[iAy + nb*(iDz)] += Factor*(r87);
   dVdXX[iAz + nb*(iAx)] += Factor*(r46);
   dVdXX[iAz + nb*(iAy)] += Factor*(r78);
   dVdXX[iAz + nb*(iAz)] += Factor*(sqr(r45)*k);
   dVdXX[iAz + nb*(iBx)] += Factor*(r88);
   dVdXX[iAz + nb*(iBy)] += Factor*(r89);
   dVdXX[iAz + nb*(iBz)] += Factor*(r90);
   dVdXX[iAz + nb*(iCx)] += Factor*(r91);
   dVdXX[iAz + nb*(iCy)] += Factor*(r92);
   dVdXX[iAz + nb*(iCz)] += Factor*(r93);
   dVdXX[iAz + nb*(iDx)] += Factor*(r94);
   dVdXX[iAz + nb*(iDy)] += Factor*(r95);
   dVdXX[iAz + nb*(iDz)] += Factor*(r96);
   dVdXX[iBx + nb*(iAx)] += Factor*(r51);
   dVdXX[iBx + nb*(iAy)] += Factor*(r79);
   dVdXX[iBx + nb*(iAz)] += Factor*(r88);
   dVdXX[iBx + nb*(iBx)] += Factor*(sqr(r50)*k);
   dVdXX[iBx + nb*(iBy)] += Factor*(r97);
   dVdXX[iBx + nb*(iBz)] += Factor*(r98);
   dVdXX[iBx + nb*(iCx)] += Factor*(r99);
   dVdXX[iBx + nb*(iCy)] += Factor*(r100);
   dVdXX[iBx + nb*(iCz)] += Factor*(r101);
   dVdXX[iBx + nb*(iDx)] += Factor*(r102);
   dVdXX[iBx + nb*(iDy)] += Factor*(r103);
   dVdXX[iBx + nb*(iDz)] += Factor*(r104);
   dVdXX[iBy + nb*(iAx)] += Factor*(r55);
   dVdXX[iBy + nb*(iAy)] += Factor*(r80);
   dVdXX[iBy + nb*(iAz)] += Factor*(r89);
   dVdXX[iBy + nb*(iBx)] += Factor*(r97);
   dVdXX[iBy + nb*(iBy)] += Factor*(sqr(r54)*k);
   dVdXX[iBy + nb*(iBz)] += Factor*(r105);
   dVdXX[iBy + nb*(iCx)] += Factor*(r106);
   dVdXX[iBy + nb*(iCy)] += Factor*(r107);
   dVdXX[iBy + nb*(iCz)] += Factor*(r108);
   dVdXX[iBy + nb*(iDx)] += Factor*(r109);
   dVdXX[iBy + nb*(iDy)] += Factor*(r110);
   dVdXX[iBy + nb*(iDz)] += Factor*(r111);
   dVdXX[iBz + nb*(iAx)] += Factor*(r59);
   dVdXX[iBz + nb*(iAy)] += Factor*(r81);
   dVdXX[iBz + nb*(iAz)] += Factor*(r90);
   dVdXX[iBz + nb*(iBx)] += Factor*(r98);
   dVdXX[iBz + nb*(iBy)] += Factor*(r105);
   dVdXX[iBz + nb*(iBz)] += Factor*(sqr(r58)*k);
   dVdXX[iBz + nb*(iCx)] += Factor*(r112);
   dVdXX[iBz + nb*(iCy)] += Factor*(r113);
   dVdXX[iBz + nb*(iCz)] += Factor*(r114);
   dVdXX[iBz + nb*(iDx)] += Factor*(r115);
   dVdXX[iBz + nb*(iDy)] += Factor*(r116);
   dVdXX[iBz + nb*(iDz)] += Factor*(r117);
   dVdXX[iCx + nb*(iAx)] += Factor*(r63);
   dVdXX[iCx + nb*(iAy)] += Factor*(r82);
   dVdXX[iCx + nb*(iAz)] += Factor*(r91);
   dVdXX[iCx + nb*(iBx)] += Factor*(r99);
   dVdXX[iCx + nb*(iBy)] += Factor*(r106);
   dVdXX[iCx + nb*(iBz)] += Factor*(r112);
   dVdXX[iCx + nb*(iCx)] += Factor*(sqr(r62)*k);
   dVdXX[iCx + nb*(iCy)] += Factor*(r118);
   dVdXX[iCx + nb*(iCz)] += Factor*(r119);
   dVdXX[iCx + nb*(iDx)] += Factor*(r120);
   dVdXX[iCx + nb*(iDy)] += Factor*(r121);
   dVdXX[iCx + nb*(iDz)] += Factor*(r122);
   dVdXX[iCy + nb*(iAx)] += Factor*(r67);
   dVdXX[iCy + nb*(iAy)] += Factor*(r83);
   dVdXX[iCy + nb*(iAz)] += Factor*(r92);
   dVdXX[iCy + nb*(iBx)] += Factor*(r100);
   dVdXX[iCy + nb*(iBy)] += Factor*(r107);
   dVdXX[iCy + nb*(iBz)] += Factor*(r113);
   dVdXX[iCy + nb*(iCx)] += Factor*(r118);
   dVdXX[iCy + nb*(iCy)] += Factor*(sqr(r66)*k);
   dVdXX[iCy + nb*(iCz)] += Factor*(r123);
   dVdXX[iCy + nb*(iDx)] += Factor*(r124);
   dVdXX[iCy + nb*(iDy)] += Factor*(r125);
   dVdXX[iCy + nb*(iDz)] += Factor*(r126);
   dVdXX[iCz + nb*(iAx)] += Factor*(r71);
   dVdXX[iCz + nb*(iAy)] += Factor*(r84);
   dVdXX[iCz + nb*(iAz)] += Factor*(r93);
   dVdXX[iCz + nb*(iBx)] += Factor*(r101);
   dVdXX[iCz + nb*(iBy)] += Factor*(r108);
   dVdXX[iCz + nb*(iBz)] += Factor*(r114);
   dVdXX[iCz + nb*(iCx)] += Factor*(r119);
   dVdXX[iCz + nb*(iCy)] += Factor*(r123);
   dVdXX[iCz + nb*(iCz)] += Factor*(sqr(r70)*k);
   dVdXX[iCz + nb*(iDx)] += Factor*(r127);
   dVdXX[iCz + nb*(iDy)] += Factor*(r128);
   dVdXX[iCz + nb*(iDz)] += Factor*(r129);
   dVdXX[iDx + nb*(iAx)] += Factor*(r73);
   dVdXX[iDx + nb*(iAy)] += Factor*(r85);
   dVdXX[iDx + nb*(iAz)] += Factor*(r94);
   dVdXX[iDx + nb*(iBx)] += Factor*(r102);
   dVdXX[iDx + nb*(iBy)] += Factor*(r109);
   dVdXX[iDx + nb*(iBz)] += Factor*(r115);
   dVdXX[iDx + nb*(iCx)] += Factor*(r120);
   dVdXX[iDx + nb*(iCy)] += Factor*(r124);
   dVdXX[iDx + nb*(iCz)] += Factor*(r127);
   dVdXX[iDx + nb*(iDx)] += Factor*(sqr(r72)*k);
   dVdXX[iDx + nb*(iDy)] += Factor*(r130);
   dVdXX[iDx + nb*(iDz)] += Factor*(r131);
   dVdXX[iDy + nb*(iAx)] += Factor*(r75);
   dVdXX[iDy + nb*(iAy)] += Factor*(r86);
   dVdXX[iDy + nb*(iAz)] += Factor*(r95);
   dVdXX[iDy + nb*(iBx)] += Factor*(r103);
   dVdXX[iDy + nb*(iBy)] += Factor*(r110);
   dVdXX[iDy + nb*(iBz)] += Factor*(r116);
   dVdXX[iDy + nb*(iCx)] += Factor*(r121);
   dVdXX[iDy + nb*(iCy)] += Factor*(r125);
   dVdXX[iDy + nb*(iCz)] += Factor*(r128);
   dVdXX[iDy + nb*(iDx)] += Factor*(r130);
   dVdXX[iDy + nb*(iDy)] += Factor*(sqr(r74)*k);
   dVdXX[iDy + nb*(iDz)] += Factor*(r132);
   dVdXX[iDz + nb*(iAx)] += Factor*(r77);
   dVdXX[iDz + nb*(iAy)] += Factor*(r87);
   dVdXX[iDz + nb*(iAz)] += Factor*(r96);
   dVdXX[iDz + nb*(iBx)] += Factor*(r104);
   dVdXX[iDz + nb*(iBy)] += Factor*(r111);
   dVdXX[iDz + nb*(iBz)] += Factor*(r117);
   dVdXX[iDz + nb*(iCx)] += Factor*(r122);
   dVdXX[iDz + nb*(iCy)] += Factor*(r126);
   dVdXX[iDz + nb*(iCz)] += Factor*(r129);
   dVdXX[iDz + nb*(iDx)] += Factor*(r131);
   dVdXX[iDz + nb*(iDy)] += Factor*(r132);
   dVdXX[iDz + nb*(iDz)] += Factor*(sqr(r76)*k);
   // suppress potential unused variable warnings.
   (void)Ax; (void)Ay; (void)Az; (void)Bx; (void)By; (void)Bz; (void)Cx; (void)Cy;
   (void)Cz; (void)Dx; (void)Dy; (void)Dz; (void)theta0; (void)k;
}

/* kate: syntax c++; */

