// bond bend harmonic energy
void eh_bend_hess(void *, double *dVdXX, index_t nb, double const &Factor, double const *pXyz, index_t const *iAtom, double theta0, double kijk)
{
   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;
   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 r1 = -1.0*Bx;
   double r2 = r1 + Ax;
   double r3 = -1.0*By;
   double r4 = r3 + Ay;
   double r5 = -1.0*Bz;
   double r6 = r5 + Az;
   double r7 = sqr(r6) + sqr(r4) + sqr(r2);
   double r8 = Cx + r1;
   double r9 = Cy + r3;
   double r10 = Cz + r5;
   double r11 = r6*r10 + r4*r9 + r2*r8;
   double r12 = sqr(r10) + sqr(r9) + sqr(r8);
   double r13 = 1.0/(1.0 - (1.0*sqr(r11))/(r12*r7));
   double r14 = sqrt(r7);
   double r15 = 1.0/r14;
   double r16 = sqrt(r12);
   double r17 = 1.0/r16;
   double r18 = 1.0/pow(r14,3);
   double r19 = r15*r8*r17 - 1.0*r11*r17*r18*r2;
   double r20 = r15*r9*r17 - 1.0*r11*r17*r18*r4;
   double r21 = r13*r19*r20*kijk;
   double r22 = r15*r10*r17 - 1.0*r11*r17*r18*r6;
   double r23 = r13*r19*r22*kijk;
   double r24 = 1.0/pow(r16,3);
   double r25 = r15*r17*((-1.0*Cx) + 2.0*Bx - 1.0*Ax) + r15*r8*r11*r24 + r2*r18*r11*r17;
   double r26 = r13*r19*r25*kijk;
   double r27 = r15*r17*((-1.0*Cy) + 2.0*By - 1.0*Ay) + r15*r9*r11*r24 + r4*r18*r11*r17;
   double r28 = r13*r19*r27*kijk;
   double r29 = r15*r17*((-1.0*Cz) + 2.0*Bz - 1.0*Az) + r15*r10*r11*r24 + r18*r6*r11*r17;
   double r30 = r13*r19*r29*kijk;
   double r31 = r2*r15*r17 - 1.0*r11*r15*r24*r8;
   double r32 = r13*r31*r19*kijk;
   double r33 = r4*r15*r17 - 1.0*r11*r15*r24*r9;
   double r34 = r13*r33*r19*kijk;
   double r35 = r15*r6*r17 - 1.0*r10*r11*r15*r24;
   double r36 = r13*r35*r19*kijk;
   double r37 = r13*r20*r22*kijk;
   double r38 = r13*r25*r20*kijk;
   double r39 = r13*r20*r27*kijk;
   double r40 = r13*r20*r29*kijk;
   double r41 = r13*r31*r20*kijk;
   double r42 = r13*r33*r20*kijk;
   double r43 = r13*r35*r20*kijk;
   double r44 = r13*r25*r22*kijk;
   double r45 = r13*r27*r22*kijk;
   double r46 = r13*r22*r29*kijk;
   double r47 = r13*r31*r22*kijk;
   double r48 = r13*r33*r22*kijk;
   double r49 = r13*r35*r22*kijk;
   double r50 = r13*r25*r27*kijk;
   double r51 = r13*r25*r29*kijk;
   double r52 = r13*r31*r25*kijk;
   double r53 = r13*r33*r25*kijk;
   double r54 = r13*r35*r25*kijk;
   double r55 = r13*r27*r29*kijk;
   double r56 = r13*r31*r27*kijk;
   double r57 = r13*r33*r27*kijk;
   double r58 = r13*r35*r27*kijk;
   double r59 = r13*r31*r29*kijk;
   double r60 = r13*r33*r29*kijk;
   double r61 = r13*r35*r29*kijk;
   double r62 = r13*r31*r33*kijk;
   double r63 = r13*r31*r35*kijk;
   double r64 = r13*r33*r35*kijk;
   dVdXX[iAx + nb*(iAx)] += Factor*(r13*sqr(r19)*kijk);
   dVdXX[iAx + nb*(iAy)] += Factor*(r21);
   dVdXX[iAx + nb*(iAz)] += Factor*(r23);
   dVdXX[iAx + nb*(iBx)] += Factor*(r26);
   dVdXX[iAx + nb*(iBy)] += Factor*(r28);
   dVdXX[iAx + nb*(iBz)] += Factor*(r30);
   dVdXX[iAx + nb*(iCx)] += Factor*(r32);
   dVdXX[iAx + nb*(iCy)] += Factor*(r34);
   dVdXX[iAx + nb*(iCz)] += Factor*(r36);
   dVdXX[iAy + nb*(iAx)] += Factor*(r21);
   dVdXX[iAy + nb*(iAy)] += Factor*(r13*sqr(r20)*kijk);
   dVdXX[iAy + nb*(iAz)] += Factor*(r37);
   dVdXX[iAy + nb*(iBx)] += Factor*(r38);
   dVdXX[iAy + nb*(iBy)] += Factor*(r39);
   dVdXX[iAy + nb*(iBz)] += Factor*(r40);
   dVdXX[iAy + nb*(iCx)] += Factor*(r41);
   dVdXX[iAy + nb*(iCy)] += Factor*(r42);
   dVdXX[iAy + nb*(iCz)] += Factor*(r43);
   dVdXX[iAz + nb*(iAx)] += Factor*(r23);
   dVdXX[iAz + nb*(iAy)] += Factor*(r37);
   dVdXX[iAz + nb*(iAz)] += Factor*(r13*sqr(r22)*kijk);
   dVdXX[iAz + nb*(iBx)] += Factor*(r44);
   dVdXX[iAz + nb*(iBy)] += Factor*(r45);
   dVdXX[iAz + nb*(iBz)] += Factor*(r46);
   dVdXX[iAz + nb*(iCx)] += Factor*(r47);
   dVdXX[iAz + nb*(iCy)] += Factor*(r48);
   dVdXX[iAz + nb*(iCz)] += Factor*(r49);
   dVdXX[iBx + nb*(iAx)] += Factor*(r26);
   dVdXX[iBx + nb*(iAy)] += Factor*(r38);
   dVdXX[iBx + nb*(iAz)] += Factor*(r44);
   dVdXX[iBx + nb*(iBx)] += Factor*(r13*sqr(r25)*kijk);
   dVdXX[iBx + nb*(iBy)] += Factor*(r50);
   dVdXX[iBx + nb*(iBz)] += Factor*(r51);
   dVdXX[iBx + nb*(iCx)] += Factor*(r52);
   dVdXX[iBx + nb*(iCy)] += Factor*(r53);
   dVdXX[iBx + nb*(iCz)] += Factor*(r54);
   dVdXX[iBy + nb*(iAx)] += Factor*(r28);
   dVdXX[iBy + nb*(iAy)] += Factor*(r39);
   dVdXX[iBy + nb*(iAz)] += Factor*(r45);
   dVdXX[iBy + nb*(iBx)] += Factor*(r50);
   dVdXX[iBy + nb*(iBy)] += Factor*(r13*sqr(r27)*kijk);
   dVdXX[iBy + nb*(iBz)] += Factor*(r55);
   dVdXX[iBy + nb*(iCx)] += Factor*(r56);
   dVdXX[iBy + nb*(iCy)] += Factor*(r57);
   dVdXX[iBy + nb*(iCz)] += Factor*(r58);
   dVdXX[iBz + nb*(iAx)] += Factor*(r30);
   dVdXX[iBz + nb*(iAy)] += Factor*(r40);
   dVdXX[iBz + nb*(iAz)] += Factor*(r46);
   dVdXX[iBz + nb*(iBx)] += Factor*(r51);
   dVdXX[iBz + nb*(iBy)] += Factor*(r55);
   dVdXX[iBz + nb*(iBz)] += Factor*(r13*sqr(r29)*kijk);
   dVdXX[iBz + nb*(iCx)] += Factor*(r59);
   dVdXX[iBz + nb*(iCy)] += Factor*(r60);
   dVdXX[iBz + nb*(iCz)] += Factor*(r61);
   dVdXX[iCx + nb*(iAx)] += Factor*(r32);
   dVdXX[iCx + nb*(iAy)] += Factor*(r41);
   dVdXX[iCx + nb*(iAz)] += Factor*(r47);
   dVdXX[iCx + nb*(iBx)] += Factor*(r52);
   dVdXX[iCx + nb*(iBy)] += Factor*(r56);
   dVdXX[iCx + nb*(iBz)] += Factor*(r59);
   dVdXX[iCx + nb*(iCx)] += Factor*(r13*sqr(r31)*kijk);
   dVdXX[iCx + nb*(iCy)] += Factor*(r62);
   dVdXX[iCx + nb*(iCz)] += Factor*(r63);
   dVdXX[iCy + nb*(iAx)] += Factor*(r34);
   dVdXX[iCy + nb*(iAy)] += Factor*(r42);
   dVdXX[iCy + nb*(iAz)] += Factor*(r48);
   dVdXX[iCy + nb*(iBx)] += Factor*(r53);
   dVdXX[iCy + nb*(iBy)] += Factor*(r57);
   dVdXX[iCy + nb*(iBz)] += Factor*(r60);
   dVdXX[iCy + nb*(iCx)] += Factor*(r62);
   dVdXX[iCy + nb*(iCy)] += Factor*(r13*sqr(r33)*kijk);
   dVdXX[iCy + nb*(iCz)] += Factor*(r64);
   dVdXX[iCz + nb*(iAx)] += Factor*(r36);
   dVdXX[iCz + nb*(iAy)] += Factor*(r43);
   dVdXX[iCz + nb*(iAz)] += Factor*(r49);
   dVdXX[iCz + nb*(iBx)] += Factor*(r54);
   dVdXX[iCz + nb*(iBy)] += Factor*(r58);
   dVdXX[iCz + nb*(iBz)] += Factor*(r61);
   dVdXX[iCz + nb*(iCx)] += Factor*(r63);
   dVdXX[iCz + nb*(iCy)] += Factor*(r64);
   dVdXX[iCz + nb*(iCz)] += Factor*(r13*sqr(r35)*kijk);
   // suppress potential unused variable warnings.
   (void)Ax; (void)Ay; (void)Az; (void)Bx; (void)By; (void)Bz; (void)Cx; (void)Cy;
   (void)Cz; (void)theta0; (void)kijk;
}

/* kate: syntax c++; */

