// distance between center of positions (AB) and (CD)
void pc_dist_cop22_hess(void *, double *dVdXX, index_t nb, double const &Factor, double const *pXyz, index_t const *iAtom)
{
   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 = Bx + Ax;
   double r2 = Dx + Cx;
   double r3 = 0.5*r1 - 0.5*r2;
   double r4 = sqr(r3);
   double r5 = By + Ay;
   double r6 = Dy + Cy;
   double r7 = 0.5*r5 - 0.5*r6;
   double r8 = sqr(r7);
   double r9 = Bz + Az;
   double r10 = Dz + Cz;
   double r11 = 0.5*r9 - 0.5*r10;
   double r12 = sqr(r11);
   double r13 = sqrt(r12 + r8 + r4);
   double r14 = 1.0/pow(r13,3);
   double r15 = 1.0/r13;
   double r16 = 0.25*r15;
   double r17 = r16 - 0.25*r14*r4;
   double r18 = -0.25*r14*r3*r7;
   double r19 = -0.25*r11*r14*r3;
   double r20 = 0.5*r2 - 0.5*r1;
   double r21 = -0.25*r15;
   double r22 = r21 - 0.25*r14*r20*r3;
   double r23 = 0.5*r6 - 0.5*r5;
   double r24 = -0.25*r14*r23*r3;
   double r25 = 0.5*r10 - 0.5*r9;
   double r26 = -0.25*r14*r25*r3;
   double r27 = r16 - 0.25*r14*r8;
   double r28 = -0.25*r11*r14*r7;
   double r29 = -0.25*r14*r20*r7;
   double r30 = r21 - 0.25*r14*r23*r7;
   double r31 = -0.25*r14*r25*r7;
   double r32 = r16 - 0.25*r12*r14;
   double r33 = -0.25*r11*r14*r20;
   double r34 = -0.25*r11*r14*r23;
   double r35 = r21 - 0.25*r11*r14*r25;
   double r36 = r16 - 0.25*r14*sqr(r20);
   double r37 = -0.25*r14*r20*r23;
   double r38 = -0.25*r14*r20*r25;
   double r39 = r16 - 0.25*r14*sqr(r23);
   double r40 = -0.25*r14*r23*r25;
   double r41 = r16 - 0.25*r14*sqr(r25);
   dVdXX[iAx + nb*(iAx)] += Factor*(r17);
   dVdXX[iAx + nb*(iAy)] += Factor*(r18);
   dVdXX[iAx + nb*(iAz)] += Factor*(r19);
   dVdXX[iAx + nb*(iBx)] += Factor*(r17);
   dVdXX[iAx + nb*(iBy)] += Factor*(r18);
   dVdXX[iAx + nb*(iBz)] += Factor*(r19);
   dVdXX[iAx + nb*(iCx)] += Factor*(r22);
   dVdXX[iAx + nb*(iCy)] += Factor*(r24);
   dVdXX[iAx + nb*(iCz)] += Factor*(r26);
   dVdXX[iAx + nb*(iDx)] += Factor*(r22);
   dVdXX[iAx + nb*(iDy)] += Factor*(r24);
   dVdXX[iAx + nb*(iDz)] += Factor*(r26);
   dVdXX[iAy + nb*(iAx)] += Factor*(r18);
   dVdXX[iAy + nb*(iAy)] += Factor*(r27);
   dVdXX[iAy + nb*(iAz)] += Factor*(r28);
   dVdXX[iAy + nb*(iBx)] += Factor*(r18);
   dVdXX[iAy + nb*(iBy)] += Factor*(r27);
   dVdXX[iAy + nb*(iBz)] += Factor*(r28);
   dVdXX[iAy + nb*(iCx)] += Factor*(r29);
   dVdXX[iAy + nb*(iCy)] += Factor*(r30);
   dVdXX[iAy + nb*(iCz)] += Factor*(r31);
   dVdXX[iAy + nb*(iDx)] += Factor*(r29);
   dVdXX[iAy + nb*(iDy)] += Factor*(r30);
   dVdXX[iAy + nb*(iDz)] += Factor*(r31);
   dVdXX[iAz + nb*(iAx)] += Factor*(r19);
   dVdXX[iAz + nb*(iAy)] += Factor*(r28);
   dVdXX[iAz + nb*(iAz)] += Factor*(r32);
   dVdXX[iAz + nb*(iBx)] += Factor*(r19);
   dVdXX[iAz + nb*(iBy)] += Factor*(r28);
   dVdXX[iAz + nb*(iBz)] += Factor*(r32);
   dVdXX[iAz + nb*(iCx)] += Factor*(r33);
   dVdXX[iAz + nb*(iCy)] += Factor*(r34);
   dVdXX[iAz + nb*(iCz)] += Factor*(r35);
   dVdXX[iAz + nb*(iDx)] += Factor*(r33);
   dVdXX[iAz + nb*(iDy)] += Factor*(r34);
   dVdXX[iAz + nb*(iDz)] += Factor*(r35);
   dVdXX[iBx + nb*(iAx)] += Factor*(r17);
   dVdXX[iBx + nb*(iAy)] += Factor*(r18);
   dVdXX[iBx + nb*(iAz)] += Factor*(r19);
   dVdXX[iBx + nb*(iBx)] += Factor*(r17);
   dVdXX[iBx + nb*(iBy)] += Factor*(r18);
   dVdXX[iBx + nb*(iBz)] += Factor*(r19);
   dVdXX[iBx + nb*(iCx)] += Factor*(r22);
   dVdXX[iBx + nb*(iCy)] += Factor*(r24);
   dVdXX[iBx + nb*(iCz)] += Factor*(r26);
   dVdXX[iBx + nb*(iDx)] += Factor*(r22);
   dVdXX[iBx + nb*(iDy)] += Factor*(r24);
   dVdXX[iBx + nb*(iDz)] += Factor*(r26);
   dVdXX[iBy + nb*(iAx)] += Factor*(r18);
   dVdXX[iBy + nb*(iAy)] += Factor*(r27);
   dVdXX[iBy + nb*(iAz)] += Factor*(r28);
   dVdXX[iBy + nb*(iBx)] += Factor*(r18);
   dVdXX[iBy + nb*(iBy)] += Factor*(r27);
   dVdXX[iBy + nb*(iBz)] += Factor*(r28);
   dVdXX[iBy + nb*(iCx)] += Factor*(r29);
   dVdXX[iBy + nb*(iCy)] += Factor*(r30);
   dVdXX[iBy + nb*(iCz)] += Factor*(r31);
   dVdXX[iBy + nb*(iDx)] += Factor*(r29);
   dVdXX[iBy + nb*(iDy)] += Factor*(r30);
   dVdXX[iBy + nb*(iDz)] += Factor*(r31);
   dVdXX[iBz + nb*(iAx)] += Factor*(r19);
   dVdXX[iBz + nb*(iAy)] += Factor*(r28);
   dVdXX[iBz + nb*(iAz)] += Factor*(r32);
   dVdXX[iBz + nb*(iBx)] += Factor*(r19);
   dVdXX[iBz + nb*(iBy)] += Factor*(r28);
   dVdXX[iBz + nb*(iBz)] += Factor*(r32);
   dVdXX[iBz + nb*(iCx)] += Factor*(r33);
   dVdXX[iBz + nb*(iCy)] += Factor*(r34);
   dVdXX[iBz + nb*(iCz)] += Factor*(r35);
   dVdXX[iBz + nb*(iDx)] += Factor*(r33);
   dVdXX[iBz + nb*(iDy)] += Factor*(r34);
   dVdXX[iBz + nb*(iDz)] += Factor*(r35);
   dVdXX[iCx + nb*(iAx)] += Factor*(r22);
   dVdXX[iCx + nb*(iAy)] += Factor*(r29);
   dVdXX[iCx + nb*(iAz)] += Factor*(r33);
   dVdXX[iCx + nb*(iBx)] += Factor*(r22);
   dVdXX[iCx + nb*(iBy)] += Factor*(r29);
   dVdXX[iCx + nb*(iBz)] += Factor*(r33);
   dVdXX[iCx + nb*(iCx)] += Factor*(r36);
   dVdXX[iCx + nb*(iCy)] += Factor*(r37);
   dVdXX[iCx + nb*(iCz)] += Factor*(r38);
   dVdXX[iCx + nb*(iDx)] += Factor*(r36);
   dVdXX[iCx + nb*(iDy)] += Factor*(r37);
   dVdXX[iCx + nb*(iDz)] += Factor*(r38);
   dVdXX[iCy + nb*(iAx)] += Factor*(r24);
   dVdXX[iCy + nb*(iAy)] += Factor*(r30);
   dVdXX[iCy + nb*(iAz)] += Factor*(r34);
   dVdXX[iCy + nb*(iBx)] += Factor*(r24);
   dVdXX[iCy + nb*(iBy)] += Factor*(r30);
   dVdXX[iCy + nb*(iBz)] += Factor*(r34);
   dVdXX[iCy + nb*(iCx)] += Factor*(r37);
   dVdXX[iCy + nb*(iCy)] += Factor*(r39);
   dVdXX[iCy + nb*(iCz)] += Factor*(r40);
   dVdXX[iCy + nb*(iDx)] += Factor*(r37);
   dVdXX[iCy + nb*(iDy)] += Factor*(r39);
   dVdXX[iCy + nb*(iDz)] += Factor*(r40);
   dVdXX[iCz + nb*(iAx)] += Factor*(r26);
   dVdXX[iCz + nb*(iAy)] += Factor*(r31);
   dVdXX[iCz + nb*(iAz)] += Factor*(r35);
   dVdXX[iCz + nb*(iBx)] += Factor*(r26);
   dVdXX[iCz + nb*(iBy)] += Factor*(r31);
   dVdXX[iCz + nb*(iBz)] += Factor*(r35);
   dVdXX[iCz + nb*(iCx)] += Factor*(r38);
   dVdXX[iCz + nb*(iCy)] += Factor*(r40);
   dVdXX[iCz + nb*(iCz)] += Factor*(r41);
   dVdXX[iCz + nb*(iDx)] += Factor*(r38);
   dVdXX[iCz + nb*(iDy)] += Factor*(r40);
   dVdXX[iCz + nb*(iDz)] += Factor*(r41);
   dVdXX[iDx + nb*(iAx)] += Factor*(r22);
   dVdXX[iDx + nb*(iAy)] += Factor*(r29);
   dVdXX[iDx + nb*(iAz)] += Factor*(r33);
   dVdXX[iDx + nb*(iBx)] += Factor*(r22);
   dVdXX[iDx + nb*(iBy)] += Factor*(r29);
   dVdXX[iDx + nb*(iBz)] += Factor*(r33);
   dVdXX[iDx + nb*(iCx)] += Factor*(r36);
   dVdXX[iDx + nb*(iCy)] += Factor*(r37);
   dVdXX[iDx + nb*(iCz)] += Factor*(r38);
   dVdXX[iDx + nb*(iDx)] += Factor*(r36);
   dVdXX[iDx + nb*(iDy)] += Factor*(r37);
   dVdXX[iDx + nb*(iDz)] += Factor*(r38);
   dVdXX[iDy + nb*(iAx)] += Factor*(r24);
   dVdXX[iDy + nb*(iAy)] += Factor*(r30);
   dVdXX[iDy + nb*(iAz)] += Factor*(r34);
   dVdXX[iDy + nb*(iBx)] += Factor*(r24);
   dVdXX[iDy + nb*(iBy)] += Factor*(r30);
   dVdXX[iDy + nb*(iBz)] += Factor*(r34);
   dVdXX[iDy + nb*(iCx)] += Factor*(r37);
   dVdXX[iDy + nb*(iCy)] += Factor*(r39);
   dVdXX[iDy + nb*(iCz)] += Factor*(r40);
   dVdXX[iDy + nb*(iDx)] += Factor*(r37);
   dVdXX[iDy + nb*(iDy)] += Factor*(r39);
   dVdXX[iDy + nb*(iDz)] += Factor*(r40);
   dVdXX[iDz + nb*(iAx)] += Factor*(r26);
   dVdXX[iDz + nb*(iAy)] += Factor*(r31);
   dVdXX[iDz + nb*(iAz)] += Factor*(r35);
   dVdXX[iDz + nb*(iBx)] += Factor*(r26);
   dVdXX[iDz + nb*(iBy)] += Factor*(r31);
   dVdXX[iDz + nb*(iBz)] += Factor*(r35);
   dVdXX[iDz + nb*(iCx)] += Factor*(r38);
   dVdXX[iDz + nb*(iCy)] += Factor*(r40);
   dVdXX[iDz + nb*(iCz)] += Factor*(r41);
   dVdXX[iDz + nb*(iDx)] += Factor*(r38);
   dVdXX[iDz + nb*(iDy)] += Factor*(r40);
   dVdXX[iDz + nb*(iDz)] += Factor*(r41);
   // 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;
}

/* kate: syntax c++; */

