/* -------------------------------------------------------------------------- *
 *                                   OpenMM                                   *
 * -------------------------------------------------------------------------- *
 * This is part of the OpenMM molecular simulation toolkit originating from   *
 * Simbios, the NIH National Center for Physics-Based Simulation of           *
 * Biological Structures at Stanford, funded under the NIH Roadmap for        *
 * Medical Research, grant U54 GM072970. See https://simtk.org.               *
 *                                                                            *
 * Portions copyright (c) 2009 Stanford University and the Authors.           *
 * Authors: Scott Le Grand, Peter Eastman                                     *
 * Contributors:                                                              *
 *                                                                            *
 * This program is free software: you can redistribute it and/or modify       *
 * it under the terms of the GNU Lesser General Public License as published   *
 * by the Free Software Foundation, either version 3 of the License, or       *
 * (at your option) any later version.                                        *
 *                                                                            *
 * This program is distributed in the hope that it will be useful,            *
 * but WITHOUT ANY WARRANTY; without even the implied warranty of             *
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the              *
 * GNU Lesser General Public License for more details.                        *
 *                                                                            *
 * You should have received a copy of the GNU Lesser General Public License   *
 * along with this program.  If not, see <http://www.gnu.org/licenses/>.      *
 * -------------------------------------------------------------------------- */

#include "cudaKernels.h"
#include "amoebaCudaKernels.h"
#include "kCalculateAmoebaCudaUtilities.h"
#include "kCalculateAmoebaCudaKirkwoodParticle.h"
extern void kCalculateObcGbsaForces2(gpuContext gpu);

static __constant__ cudaGmxSimulation cSim;
static __constant__ cudaAmoebaGmxSimulation cAmoebaSim;

void SetCalculateAmoebaKirkwoodSim(amoebaGpuContext amoebaGpu)
{
    cudaError_t status;
    gpuContext gpu = amoebaGpu->gpuContext;
    status         = cudaMemcpyToSymbol(cSim, &gpu->sim, sizeof(cudaGmxSimulation));
    RTERROR(status, "SetCalculateAmoebaKirkwoodSim: cudaMemcpyToSymbol: SetSim copy to cSim failed");
    status         = cudaMemcpyToSymbol(cAmoebaSim, &amoebaGpu->amoebaSim, sizeof(cudaAmoebaGmxSimulation));
    RTERROR(status, "SetCalculateAmoebaKirkwoodSim: cudaMemcpyToSymbol: SetSim copy to cAmoebaSim failed");
}

void GetCalculateAmoebaKirkwoodSim(amoebaGpuContext amoebaGpu)
{
    cudaError_t status;
    gpuContext gpu = amoebaGpu->gpuContext;
    status = cudaMemcpyFromSymbol(&gpu->sim, cSim, sizeof(cudaGmxSimulation));
    RTERROR(status, "GetCalculateAmoebaKirkwoodSim: cudaMemcpyFromSymbol: SetSim copy from cSim failed");
    status = cudaMemcpyFromSymbol(&amoebaGpu->amoebaSim, cAmoebaSim, sizeof(cudaAmoebaGmxSimulation));
    RTERROR(status, "GetCalculateAmoebaKirkwoodSim: cudaMemcpyFromSymbol: SetSim copy from cAmoebaSim failed");
}

__device__ void loadKirkwoodShared( struct KirkwoodParticle* sA, unsigned int atomI )
{
    // coordinates & charge

    sA->x                        = cSim.pPosq[atomI].x;
    sA->y                        = cSim.pPosq[atomI].y;
    sA->z                        = cSim.pPosq[atomI].z;
    sA->q                        = cSim.pPosq[atomI].w;

    // lab dipole

    sA->labFrameDipole[0]        = cAmoebaSim.pLabFrameDipole[atomI*3];
    sA->labFrameDipole[1]        = cAmoebaSim.pLabFrameDipole[atomI*3+1];
    sA->labFrameDipole[2]        = cAmoebaSim.pLabFrameDipole[atomI*3+2];

    // lab quadrupole

    sA->labFrameQuadrupole_XX    = cAmoebaSim.pLabFrameQuadrupole[atomI*9];
    sA->labFrameQuadrupole_XY    = cAmoebaSim.pLabFrameQuadrupole[atomI*9+1];
    sA->labFrameQuadrupole_XZ    = cAmoebaSim.pLabFrameQuadrupole[atomI*9+2];
    sA->labFrameQuadrupole_YY    = cAmoebaSim.pLabFrameQuadrupole[atomI*9+4];
    sA->labFrameQuadrupole_YZ    = cAmoebaSim.pLabFrameQuadrupole[atomI*9+5];
    //sA->labFrameQuadrupole_ZZ    = cAmoebaSim.pLabFrameQuadrupole[atomI*9+8];
    sA->labFrameQuadrupole_ZZ    = -(sA->labFrameQuadrupole_XX + sA->labFrameQuadrupole_YY);

    // induced dipole

    sA->inducedDipole[0]         = cAmoebaSim.pInducedDipoleS[atomI*3];
    sA->inducedDipole[1]         = cAmoebaSim.pInducedDipoleS[atomI*3+1];
    sA->inducedDipole[2]         = cAmoebaSim.pInducedDipoleS[atomI*3+2];

    // induced dipole polar

    sA->inducedDipoleP[0]        = cAmoebaSim.pInducedDipolePolarS[atomI*3];
    sA->inducedDipoleP[1]        = cAmoebaSim.pInducedDipolePolarS[atomI*3+1];
    sA->inducedDipoleP[2]        = cAmoebaSim.pInducedDipolePolarS[atomI*3+2];

    sA->bornRadius               = cSim.pBornRadii[atomI];

}

__device__ void calculateKirkwoodPairIxnOrig_kernel( KirkwoodParticle& atomI,       KirkwoodParticle& atomJ,
                                                 float*  outputForce,           float outputTorque[2][3],
                                                 float*  outputBorn,            float*  outputBornPolar,
                                                 float* outputEnergy

 ){

    float e,ei;
    float xr,yr,zr;
    float xr2,yr2,zr2;
    float sxi,syi,szi;
    float sxk,syk,szk;
    float r2,rb2;
    float dedx,dedy,dedz;
    float drbi;
    float drbk;
    float dpdx,dpdy,dpdz;
    float dpbi;
    float dpbk;
    float fc,fd,fq;
    float expterm;
    float gf,gf2,gf3,gf5;
    float gf7,gf9,gf11;
    float expc,dexpc;
    float expc1,expcdexpc;
    float expcr,dexpcr;
    float dgfdr;
    float esym,ewi,ewk;
    float desymdx,dewidx,dewkdx;
    float desymdy,dewidy,dewkdy;
    float desymdz,dewidz,dewkdz;
    float dsumdr,desymdr;
    float dewidr,dewkdr;
    float dsymdr;
    float esymi,ewii,ewki;
    float dpwidx,dpwkdx;
    float dpsymdy,dpwidy,dpwkdy;
    float dpsymdz,dpwidz,dpwkdz;
    float dwipdr,dwkpdr;
    float duvdr;

    // set the bulk dielectric constant to the water value

    fc           = cAmoebaSim.electric * cAmoebaSim.fc;
    fd           = cAmoebaSim.electric * cAmoebaSim.fd;
    fq           = cAmoebaSim.electric * cAmoebaSim.fq;

    sxi          = atomI.inducedDipole[0] + atomI.inducedDipoleP[0];
    syi          = atomI.inducedDipole[1] + atomI.inducedDipoleP[1];
    szi          = atomI.inducedDipole[2] + atomI.inducedDipoleP[2];

    // decide whether to compute the current interaction;

    xr           = atomJ.x - atomI.x;
    yr           = atomJ.y - atomI.y;
    zr           = atomJ.z - atomI.z;

    xr2          = xr*xr;
    yr2          = yr*yr;
    zr2          = zr*zr;
    r2           = xr2 + yr2 + zr2;

    if( r2 > cAmoebaSim.scalingDistanceCutoff ){
    }

    sxk          = atomJ.inducedDipole[0] + atomJ.inducedDipoleP[0];
    syk          = atomJ.inducedDipole[1] + atomJ.inducedDipoleP[1];
    szk          = atomJ.inducedDipole[2] + atomJ.inducedDipoleP[2];

    rb2          = atomI.bornRadius*atomJ.bornRadius;

    expterm      = expf(-r2/(cAmoebaSim.gkc*rb2));
    expc         = expterm/cAmoebaSim.gkc;
    expcr        = r2*expterm/(cAmoebaSim.gkc*cAmoebaSim.gkc*rb2*rb2);
    dexpc        = -2.0f / (cAmoebaSim.gkc*rb2);
    dexpcr       = 2.0f / (cAmoebaSim.gkc*rb2*rb2);
    dgfdr        = 0.5f * expterm * (1.0f+r2/(rb2*cAmoebaSim.gkc));
    gf2          = 1.0f / (r2+rb2*expterm);
    gf           = sqrt(gf2);
    gf3          = gf2 * gf;
    gf5          = gf3 * gf2;
    gf7          = gf5 * gf2;
    gf9          = gf7 * gf2;
    gf11         = gf9 * gf2;

    // reaction potential auxiliary terms;

    float a00      =            gf;
    float a10      =          -gf3;
    float a20      =    3.0f * gf5;
    float a30      =  -15.0f * gf7;
    float a40      =  105.0f * gf9;
    float a50      = -945.0f * gf11;

    // Born radii derivatives of reaction potential auxiliary terms;

    float b00      = dgfdr * a10;
    float b10      = dgfdr * a20;
    float b20      = dgfdr * a30;
    float b30      = dgfdr * a40;
    float b40      = dgfdr * a50;

    // reaction potential gradient auxiliary terms;

    expc1          = 1.0f - expc;
    float a01      = expc1 * a10;
    float a11      = expc1 * a20;
    float a21      = expc1 * a30;
    float a31      = expc1 * a40;
    float a41      = expc1 * a50;

    // Born radii derivs of reaction potential gradient auxiliary terms;

    float b01      = b10 - expcr*a10 - expc*b10;
    float b11      = b20 - expcr*a20 - expc*b20;
    float b21      = b30 - expcr*a30 - expc*b30;
    float b31      = b40 - expcr*a40 - expc*b40;

    // 2nd reaction potential gradient auxiliary terms;

    expcdexpc      = -expc * dexpc;
    float a02      = expc1*a11 + expcdexpc*a10;
    float a12      = expc1*a21 + expcdexpc*a20;
    float a22      = expc1*a31 + expcdexpc*a30;
    float a32      = expc1*a41 + expcdexpc*a40;

    // Born radii derivatives of the 2nd reaction potential
    // gradient auxiliary terms

     float b02     = b11 - (expcr*(a11 + dexpc*a10)
                           + expc*(b11 + dexpcr*a10+dexpc*b10));

     float b12     = b21 - (expcr*(a21 + dexpc*a20)
                           +   expc*(b21 + dexpcr*a20+dexpc*b20));

     float b22     = b31 - (expcr*(a31 + dexpc*a30)
                           +   expc*(b31 + dexpcr*a30+dexpc*b30));

    // 3rd reaction potential gradient auxiliary terms

    expcdexpc      = 2.0f * expcdexpc;
    float a03      = expc1*a12 + expcdexpc*a11;
    float a13      = expc1*a22 + expcdexpc*a21;
    float a23      = expc1*a32 + expcdexpc*a31;

    expcdexpc      = -expc*dexpc*dexpc;
    a03            = a03 + expcdexpc*a10;
    a13            = a13 + expcdexpc*a20;
    a23            = a23 + expcdexpc*a30;

    // multiply the auxillary terms by their dieletric functions;

    a00            *= fc;
    a01            *= fc;
    a02            *= fc;
    a03            *= fc;

    b00            *= fc;
    b01            *= fc;
    b02            *= fc;

    a10            *= fd;
    a11            *= fd;
    a12            *= fd;
    a13            *= fd;

    b10            *= fd;
    b11            *= fd;
    b12            *= fd;

    a20            *= fq;
    a21            *= fq;
    a22            *= fq;
    a23            *= fq;

    b20            *= fq;
    b21            *= fq;
    b22            *= fq;

    // unweighted reaction potential tensor

    float gc1        = a00;

    float gux1       = xr * a10;
    float guy1       = yr * a10;
    float guz1       = zr * a10;

    float gqxx1      = xr2 * a20;
    float gqyy1      = yr2 * a20;
    float gqzz1      = zr2 * a20;

    float gqxy1      = xr * yr * a20;
    float gqxz1      = xr * zr * a20;
    float gqyz1      = yr * zr * a20;

    // Born radii derivs of unweighted reaction potential tensor

    float gc21       = b00;

    float gux21      = xr * b10;
    float guy21      = yr * b10;
    float guz21      = zr * b10;

    float gqxx21     = xr2 * b20;
    float gqyy21     = yr2 * b20;
    float gqzz21     = zr2 * b20;

    float gqxy21     = xr * yr * b20;
    float gqxz21     = xr * zr * b20;
    float gqyz21     = yr * zr * b20;

    // unweighted reaction potential gradient tensor;

    float gc2        = xr * a01;
    float gc3        = yr * a01;
    float gc4        = zr * a01;

    float gux2       = a10 + xr2*a11;
    float gux3       = xr * yr * a11;
    float gux4       = xr * zr * a11;

    float guy2       = gux3;
    float guy3       = a10 + yr2*a11;
    float guy4       = yr * zr * a11;
    float guz2       = gux4;
    float guz3       = guy4;
    float guz4       = a10 + zr2*a11;
    float gqxx2      = xr * (2.0f*a20+xr2*a21);
    float gqxx3      = yr * xr2 * a21;
    float gqxx4      = zr * xr2 * a21;
    float gqyy2      = xr * yr2 * a21;
    float gqyy3      = yr * (2.0f*a20+yr2*a21);
    float gqyy4      = zr * yr2 * a21;
    float gqzz2      = xr * zr2 * a21;
    float gqzz3      = yr * zr2 * a21;
    float gqzz4      = zr * (2.0f*a20+zr2*a21);
    float gqxy2      = yr * (a20+xr2*a21);
    float gqxy3      = xr * (a20+yr2*a21);
    float gqxy4      = zr * xr * yr * a21;
    float gqxz2      = zr * (a20+xr2*a21);
    float gqxz3      = gqxy4;
    float gqxz4      = xr * (a20+zr2*a21);
    float gqyz2      = gqxy4;
    float gqyz3      = zr * (a20+yr2*a21);
    float gqyz4      = yr * (a20+zr2*a21);

    // Born derivs of the unweighted reaction potential gradient tensor

    float gc22       = xr * b01;
    float gc23       = yr * b01;
    float gc24       = zr * b01;
    float gux22      = b10 + xr2*b11;
    float gux23      = xr * yr * b11;
    float gux24      = xr * zr * b11;
    float guy22      = gux23;
    float guy23      = b10 + yr2*b11;
    float guy24      = yr * zr * b11;
    float guz22      = gux24;
    float guz23      = guy24;
    float guz24      = b10 + zr2*b11;
    float gqxx22     = xr * (2.0f*b20+xr2*b21);
    float gqxx23     = yr * xr2 * b21;
    float gqxx24     = zr * xr2 * b21;
    float gqyy22     = xr * yr2 * b21;
    float gqyy23     = yr * (2.0f*b20+yr2*b21);
    float gqyy24     = zr * yr2 * b21;
    float gqzz22     = xr * zr2 * b21;
    float gqzz23     = yr * zr2 * b21;
    float gqzz24     = zr * (2.0f*b20 + zr2*b21);
    float gqxy22     = yr * (b20+xr2*b21);
    float gqxy23     = xr * (b20+yr2*b21);
    float gqxy24     = zr * xr * yr * b21;
    float gqxz22     = zr * (b20+xr2*b21);
    float gqxz23     = gqxy24;
    float gqxz24     = xr * (b20+zr2*b21);
    float gqyz22     = gqxy24;
    float gqyz23     = zr * (b20+yr2*b21);
    float gqyz24     = yr * (b20+zr2*b21);

    // unweighted 2nd reaction potential gradient tensor;

    float gc5        = a01 + xr2*a02;
    float gc6        = xr * yr * a02;
    float gc7        = xr * zr * a02;
    float gc8        = a01 + yr2*a02;
    float gc9        = yr * zr * a02;
    float gc10       = a01 + zr2*a02;
    float gux5       = xr * (3.0f*a11+xr2*a12);
    float gux6       = yr * (a11+xr2*a12);
    float gux7       = zr * (a11+xr2*a12);
    float gux8       = xr * (a11+yr2*a12);
    float gux9       = zr * xr * yr * a12;
    float gux10      = xr * (a11+zr2*a12);
    float guy5       = yr * (a11+xr2*a12);
    float guy6       = xr * (a11+yr2*a12);
    float guy7       = gux9;
    float guy8       = yr * (3.0f*a11+yr2*a12);
    float guy9       = zr * (a11+yr2*a12);
    float guy10      = yr * (a11+zr2*a12);
    float guz5       = zr * (a11+xr2*a12);
    float guz6       = gux9;
    float guz7       = xr * (a11+zr2*a12);
    float guz8       = zr * (a11+yr2*a12);
    float guz9       = yr * (a11+zr2*a12);
    float guz10      = zr * (3.0f*a11+zr2*a12);
    float gqxx5      = 2.0f*a20 + xr2*(5.0f*a21+xr2*a22);
    float gqxx6      = yr * xr * (2.0f*a21+xr2*a22);
    float gqxx7      = zr * xr * (2.0f*a21+xr2*a22);
    float gqxx8      = xr2 * (a21+yr2*a22);
    float gqxx9      = zr * yr * xr2 * a22;
    float gqxx10     = xr2 * (a21+zr2*a22);
    float gqyy5      = yr2 * (a21+xr2*a22);
    float gqyy6      = xr * yr * (2.0f*a21+yr2*a22);
    float gqyy7      = xr * zr * yr2 * a22;
    float gqyy8      = 2.0f*a20 + yr2*(5.0f*a21+yr2*a22);
    float gqyy9      = yr * zr * (2.0f*a21+yr2*a22);
    float gqyy10     = yr2 * (a21+zr2*a22);
    float gqzz5      = zr2 * (a21+xr2*a22);
    float gqzz6      = xr * yr * zr2 * a22;
    float gqzz7      = xr * zr * (2.0f*a21+zr2*a22);
    float gqzz8      = zr2 * (a21+yr2*a22);
    float gqzz9      = yr * zr * (2.0f*a21+zr2*a22);
    float gqzz10     = 2.0f*a20 + zr2*(5.0f*a21+zr2*a22);
    float gqxy5      = xr * yr * (3.0f*a21+xr2*a22);
    float gqxy6      = a20 + (xr2+yr2)*a21 + xr2*yr2*a22;
    float gqxy7      = zr * yr * (a21+xr2*a22);
    float gqxy8      = xr * yr * (3.0f*a21+yr2*a22);
    float gqxy9      = zr * xr * (a21+yr2*a22);
    float gqxy10     = xr * yr * (a21+zr2*a22);
    float gqxz5      = xr * zr * (3.0f*a21+xr2*a22);
    float gqxz6      = yr * zr * (a21+xr2*a22);
    float gqxz7      = a20 + (xr2+zr2)*a21 + xr2*zr2*a22;
    float gqxz8      = xr * zr * (a21+yr2*a22);
    float gqxz9      = xr * yr * (a21+zr2*a22);
    float gqxz10     = xr * zr * (3.0f*a21+zr2*a22);
    float gqyz5      = zr * yr * (a21+xr2*a22);
    float gqyz6      = xr * zr * (a21+yr2*a22);
    float gqyz7      = xr * yr * (a21+zr2*a22);
    float gqyz8      = yr * zr * (3.0f*a21+yr2*a22);
    float gqyz9      = a20 + (yr2+zr2)*a21 + yr2*zr2*a22;
    float gqyz10     = yr * zr * (3.0f*a21+zr2*a22);

    // Born radii derivatives of the unweighted 2nd reaction;
    // potential gradient tensor;

    float gc25       = b01 + xr2*b02;
    float gc26       = xr * yr * b02;
    float gc27       = xr * zr * b02;
    float gc28       = b01 + yr2*b02;
    float gc29       = yr * zr * b02;
    float gc30       = b01 + zr2*b02;
    float gux25      = xr * (3.0f*b11+xr2*b12);
    float gux26      = yr * (b11+xr2*b12);
    float gux27      = zr * (b11+xr2*b12);
    float gux28      = xr * (b11+yr2*b12);
    float gux29      = zr * xr * yr * b12;
    float gux30      = xr * (b11+zr2*b12);
    float guy25      = yr * (b11+xr2*b12);
    float guy26      = xr * (b11+yr2*b12);
    float guy27      = gux29;
    float guy28      = yr * (3.0f*b11+yr2*b12);
    float guy29      = zr * (b11+yr2*b12);
    float guy30      = yr * (b11+zr2*b12);
    float guz25      = zr * (b11+xr2*b12);
    float guz26      = gux29;
    float guz27      = xr * (b11+zr2*b12);
    float guz28      = zr * (b11+yr2*b12);
    float guz29      = yr * (b11+zr2*b12);
    float guz30      = zr * (3.0f*b11+zr2*b12);
    float gqxx25     = 2.0f*b20 + xr2*(5.0f*b21+xr2*b22);
    float gqxx26     = yr * xr * (2.0f*b21+xr2*b22);
    float gqxx27     = zr * xr * (2.0f*b21+xr2*b22);
    float gqxx28     = xr2 * (b21+yr2*b22);
    float gqxx29     = zr * yr * xr2 * b22;
    float gqxx30     = xr2 * (b21+zr2*b22);
    float gqyy25     = yr2 * (b21+xr2*b22);
    float gqyy26     = xr * yr * (2.0f*b21+yr2*b22);
    float gqyy27     = xr * zr * yr2 * b22;
    float gqyy28     = 2.0f*b20 + yr2*(5.0f*b21+yr2*b22);
    float gqyy29     = yr * zr * (2.0f*b21+yr2*b22);
    float gqyy30     = yr2 * (b21+zr2*b22);
    float gqzz25     = zr2 * (b21+xr2*b22);
    float gqzz26     = xr * yr * zr2 * b22;
    float gqzz27     = xr * zr * (2.0f*b21+zr2*b22);
    float gqzz28     = zr2 * (b21+yr2*b22);
    float gqzz29     = yr * zr * (2.0f*b21+zr2*b22);
    float gqzz30     = 2.0f*b20 + zr2*(5.0f*b21+zr2*b22);
    float gqxy25     = xr * yr * (3.0f*b21 + xr2*b22);
    float gqxy26     = b20 + (xr2+yr2)*b21 + xr2*yr2*b22;
    float gqxy27     = zr * yr * (b21+xr2*b22);
    float gqxy28     = xr * yr * (3.0f*b21+yr2*b22);
    float gqxy29     = zr * xr * (b21+yr2*b22);
    float gqxy30     = xr * yr * (b21+zr2*b22);
    float gqxz25     = xr * zr * (3.0f*b21+xr2*b22);
    float gqxz26     = yr * zr * (b21+xr2*b22);
    float gqxz27     = b20 + (xr2+zr2)*b21 + xr2*zr2*b22;
    float gqxz28     = xr * zr * (b21+yr2*b22);
    float gqxz29     = xr * yr * (b21+zr2*b22);
    float gqxz30     = xr * zr * (3.0f*b21+zr2*b22);
    float gqyz25     = zr * yr * (b21+xr2*b22);
    float gqyz26     = xr * zr * (b21+yr2*b22);
    float gqyz27     = xr * yr * (b21+zr2*b22);
    float gqyz28     = yr * zr * (3.0f*b21+yr2*b22);
    float gqyz29     = b20 + (yr2+zr2)*b21 + yr2*zr2*b22;
    float gqyz30     = yr * zr * (3.0f*b21+zr2*b22);

    // unweighted 3rd reaction potential gradient tensor;

    float gc11       = xr * (3.0f*a02+xr2*a03);
    float gc12       = yr * (a02+xr2*a03);
    float gc13       = zr * (a02+xr2*a03);
    float gc14       = xr * (a02+yr2*a03);
    float gc15       = xr * yr * zr * a03;
    float gc16       = xr * (a02+zr2*a03);
    float gc17       = yr * (3.0f*a02+yr2*a03);
    float gc18       = zr * (a02+yr2*a03);
    float gc19       = yr * (a02+zr2*a03);
    float gc20       = zr * (3.0f*a02+zr2*a03);
    float gux11      = 3.0f*a11 + xr2*(6.0f*a12+xr2*a13);
    float gux12      = xr * yr * (3.0f*a12+xr2*a13);
    float gux13      = xr * zr * (3.0f*a12+xr2*a13);
    float gux14      = a11 + (xr2+yr2)*a12 + xr2*yr2*a13;
    float gux15      = yr * zr * (a12+xr2*a13);
    float gux16      = a11 + (xr2+zr2)*a12 + xr2*zr2*a13;
    float gux17      = xr * yr * (3.0f*a12+yr2*a13);
    float gux18      = xr * zr * (a12+yr2*a13);
    float gux19      = xr * yr * (a12+zr2*a13);
    float gux20      = xr * zr * (3.0f*a12+zr2*a13);
    float guy11      = gux12;
    float guy12      = gux14;
    float guy13      = gux15;
    float guy14      = gux17;
    float guy15      = gux18;
    float guy16      = gux19;
    float guy17      = 3.0f*a11 + yr2*(6.0f*a12+yr2*a13);
    float guy18      = yr * zr * (3.0f*a12+yr2*a13);
    float guy19      = a11 + (yr2+zr2)*a12 + yr2*zr2*a13;
    float guy20      = yr * zr * (3.0f*a12+zr2*a13);
    float guz11      = gux13;
    float guz12      = gux15;
    float guz13      = gux16;
    float guz14      = gux18;
    float guz15      = gux19;
    float guz16      = gux20;
    float guz17      = guy18;
    float guz18      = guy19;
    float guz19      = guy20;
    float guz20      = 3.0f*a11 + zr2*(6.0f*a12+zr2*a13);

    float gqxx11     = xr * (12.0f*a21+xr2*(9.0f*a22 + xr2*a23));
    float gqxx12     = yr * (2.0f*a21+xr2*(5.0f*a22  + xr2*a23));
    float gqxx13     = zr * (2.0f*a21+xr2*(5.0f*a22  + xr2*a23));
    float gqxx14     = xr * (2.0f*a21+yr2*2.0f*a22   +xr2*(a22+yr2*a23));
    float gqxx15     = xr * yr * zr * (2.0f*a22+xr2*a23);
    float gqxx16     = xr * (2.0f*a21+zr2*2.0f*a22 +xr2*(a22+zr2*a23));
    float gqxx17     = yr * xr2 * (3.0f*a22+yr2*a23);
    float gqxx18     = zr * xr2 * (a22+yr2*a23);
    float gqxx19     = yr * xr2 * (a22+zr2*a23);
    float gqxx20     = zr * xr2 * (3.0f*a22+zr2*a23);

    float gqxy11     = yr * (3.0f*a21+xr2*(6.0f*a22 +xr2*a23));
    float gqxy12     = xr * (3.0f*(a21+yr2*a22) +xr2*(a22+yr2*a23));
    float gqxy13     = xr * yr * zr * (3.0f*a22+xr2*a23);
    float gqxy14     = yr * (3.0f*(a21+xr2*a22) +yr2*(a22+xr2*a23));
    float gqxy15     = zr * (a21+(yr2+xr2)*a22 +yr2*xr2*a23);
    float gqxy16     = yr * (a21+(xr2+zr2)*a22 +xr2*zr2*a23);
    float gqxy17     = xr * (3.0f*(a21+yr2*a22) +yr2*(3.0f*a22+yr2*a23));
    float gqxy18     = xr * yr * zr * (3.0f*a22+yr2*a23);
    float gqxy19     = xr * (a21+(yr2+zr2)*a22 +yr2*zr2*a23);
    float gqxy20     = xr * yr * zr * (3.0f*a22+zr2*a23);
    float gqxz11     = zr * (3.0f*a21+xr2*(6.0f*a22 +xr2*a23));
    float gqxz12     = xr * yr * zr * (3.0f*a22+xr2*a23);
    float gqxz13     = xr * (3.0f*(a21+zr2*a22) +xr2*(a22+zr2*a23));
    float gqxz14     = zr * (a21+(xr2+yr2)*a22 +xr2*yr2*a23);
    float gqxz15     = yr * (a21+(xr2+zr2)*a22 +zr2*xr2*a23);
    float gqxz16     = zr * (3.0f*(a21+xr2*a22) +zr2*(a22+xr2*a23));
    float gqxz17     = xr * yr * zr * (3.0f*a22+yr2*a23);
    float gqxz18     = xr * (a21+(zr2+yr2)*a22 +zr2*yr2*a23);
    float gqxz19     = xr * yr * zr * (3.0f*a22+zr2*a23);
    float gqxz20     = xr * (3.0f*a21+zr2*(6.0f*a22 +zr2*a23));
    float gqyy11     = xr * yr2 * (3.0f*a22+xr2*a23);
    float gqyy12     = yr * (2.0f*a21+xr2*2.0f*a22 +yr2*(a22+xr2*a23));
    float gqyy13     = zr * yr2 * (a22+xr2*a23);
    float gqyy14     = xr * (2.0f*a21+yr2*(5.0f*a22 +yr2*a23));
    float gqyy15     = xr * yr * zr * (2.0f*a22+yr2*a23);
    float gqyy16     = xr * yr2 * (a22+zr2*a23);
    float gqyy17     = yr * (12.0f*a21+yr2*(9.0f*a22 +yr2*a23));
    float gqyy18     = zr * (2.0f*a21+yr2*(5.0f*a22 +yr2*a23));
    float gqyy19     = yr * (2.0f*a21+zr2*2.0f*a22 +yr2*(a22+zr2*a23));
    float gqyy20     = zr * yr2 * (3.0f*a22+zr2*a23);
    float gqyz11     = xr * yr * zr * (3.0f*a22+xr2*a23);
    float gqyz12     = zr * (a21+(xr2+yr2)*a22 +xr2*yr2*a23);
    float gqyz13     = yr * (a21+(xr2+zr2)*a22 +xr2*zr2*a23);
    float gqyz14     = xr * yr * zr * (3.0f*a22+yr2*a23);
    float gqyz15     = xr * (a21+(yr2+zr2)*a22 +yr2*zr2*a23);
    float gqyz16     = xr * yr * zr * (3.0f*a22+zr2*a23);
    float gqyz17     = zr * (3.0f*a21+yr2*(6.0f*a22 +yr2*a23));
    float gqyz18     = yr * (3.0f*(a21+zr2*a22) +yr2*(a22+zr2*a23));
    float gqyz19     = zr * (3.0f*(a21+yr2*a22) +zr2*(a22+yr2*a23));
    float gqyz20     = yr * (3.0f*a21+zr2*(6.0f*a22 +zr2*a23));
    float gqzz11     = xr * zr2 * (3.0f*a22+xr2*a23);
    float gqzz12     = yr * (zr2*a22+xr2*(zr2*a23));
    float gqzz13     = zr * (2.0f*a21+xr2*2.0f*a22 +zr2*(a22+xr2*a23));
    float gqzz14     = xr * zr2 * (a22+yr2*a23);
    float gqzz15     = xr * yr * zr * (2.0f*a22+zr2*a23);
    float gqzz16     = xr * (2.0f*a21+zr2*(5.0f*a22 +zr2*a23));
    float gqzz17     = yr * zr2 * (3.0f*a22+yr2*a23);
    float gqzz18     = zr * (2.0f*a21+yr2*2.0f*a22 +zr2*(a22+yr2*a23));
    float gqzz19     = yr * (2.0f*a21+zr2*(5.0f*a22 +zr2*a23));
    float gqzz20     = zr * (12.0f*a21+zr2*(9.0f*a22 +zr2*a23));

    // electrostatic solvation energy of the permanent multipoles
    // in their own GK reaction potential

    esym = atomI.q * atomJ.q * gc1 - (atomI.labFrameDipole[0]*(atomJ.labFrameDipole[0]*gux2+atomJ.labFrameDipole[1]*guy2+atomJ.labFrameDipole[2]*guz2)
                                   +  atomI.labFrameDipole[1]*(atomJ.labFrameDipole[0]*gux3+atomJ.labFrameDipole[1]*guy3+atomJ.labFrameDipole[2]*guz3)
                                   +  atomI.labFrameDipole[2]*(atomJ.labFrameDipole[0]*gux4+atomJ.labFrameDipole[1]*guy4+atomJ.labFrameDipole[2]*guz4));

    ewi =  atomI.q*(atomJ.labFrameDipole[0]*gc2+atomJ.labFrameDipole[1]*gc3+atomJ.labFrameDipole[2]*gc4)
          -atomJ.q*(atomI.labFrameDipole[0]*gux1+atomI.labFrameDipole[1]*guy1+atomI.labFrameDipole[2]*guz1)
           +atomI.q*(atomJ.labFrameQuadrupole_XX*gc5+atomJ.labFrameQuadrupole_YY*gc8+atomJ.labFrameQuadrupole_ZZ*gc10
              +2.0f*(atomJ.labFrameQuadrupole_XY*gc6+atomJ.labFrameQuadrupole_XZ*gc7+atomJ.labFrameQuadrupole_YZ*gc9))
                 +atomJ.q*(atomI.labFrameQuadrupole_XX*gqxx1+atomI.labFrameQuadrupole_YY*gqyy1+atomI.labFrameQuadrupole_ZZ*gqzz1
              +2.0f*(atomI.labFrameQuadrupole_XY*gqxy1+atomI.labFrameQuadrupole_XZ*gqxz1+atomI.labFrameQuadrupole_YZ*gqyz1))
               - atomI.labFrameDipole[0]*(atomJ.labFrameQuadrupole_XX*gux5+atomJ.labFrameQuadrupole_YY*gux8+atomJ.labFrameQuadrupole_ZZ*gux10
              +2.0f*(atomJ.labFrameQuadrupole_XY*gux6+atomJ.labFrameQuadrupole_XZ*gux7+atomJ.labFrameQuadrupole_YZ*gux9))
               - atomI.labFrameDipole[1]*(atomJ.labFrameQuadrupole_XX*guy5+atomJ.labFrameQuadrupole_YY*guy8+atomJ.labFrameQuadrupole_ZZ*guy10
              +2.0f*(atomJ.labFrameQuadrupole_XY*guy6+atomJ.labFrameQuadrupole_XZ*guy7+atomJ.labFrameQuadrupole_YZ*guy9))
               - atomI.labFrameDipole[2]*(atomJ.labFrameQuadrupole_XX*guz5+atomJ.labFrameQuadrupole_YY*guz8+atomJ.labFrameQuadrupole_ZZ*guz10
              +2.0f*(atomJ.labFrameQuadrupole_XY*guz6+atomJ.labFrameQuadrupole_XZ*guz7+atomJ.labFrameQuadrupole_YZ*guz9))
               + atomJ.labFrameDipole[0]*(atomI.labFrameQuadrupole_XX*gqxx2+atomI.labFrameQuadrupole_YY*gqyy2+atomI.labFrameQuadrupole_ZZ*gqzz2
              +2.0f*(atomI.labFrameQuadrupole_XY*gqxy2+atomI.labFrameQuadrupole_XZ*gqxz2+atomI.labFrameQuadrupole_YZ*gqyz2))
               + atomJ.labFrameDipole[1]*(atomI.labFrameQuadrupole_XX*gqxx3+atomI.labFrameQuadrupole_YY*gqyy3+atomI.labFrameQuadrupole_ZZ*gqzz3
              +2.0f*(atomI.labFrameQuadrupole_XY*gqxy3+atomI.labFrameQuadrupole_XZ*gqxz3+atomI.labFrameQuadrupole_YZ*gqyz3))
               + atomJ.labFrameDipole[2]*(atomI.labFrameQuadrupole_XX*gqxx4+atomI.labFrameQuadrupole_YY*gqyy4+atomI.labFrameQuadrupole_ZZ*gqzz4
              +2.0f*(atomI.labFrameQuadrupole_XY*gqxy4+atomI.labFrameQuadrupole_XZ*gqxz4+atomI.labFrameQuadrupole_YZ*gqyz4))
              + atomI.labFrameQuadrupole_XX*(atomJ.labFrameQuadrupole_XX*gqxx5+atomJ.labFrameQuadrupole_YY*gqxx8+atomJ.labFrameQuadrupole_ZZ*gqxx10
              +2.0f*(atomJ.labFrameQuadrupole_XY*gqxx6+atomJ.labFrameQuadrupole_XZ*gqxx7+atomJ.labFrameQuadrupole_YZ*gqxx9))
              + atomI.labFrameQuadrupole_YY*(atomJ.labFrameQuadrupole_XX*gqyy5+atomJ.labFrameQuadrupole_YY*gqyy8+atomJ.labFrameQuadrupole_ZZ*gqyy10
              +2.0f*(atomJ.labFrameQuadrupole_XY*gqyy6+atomJ.labFrameQuadrupole_XZ*gqyy7+atomJ.labFrameQuadrupole_YZ*gqyy9))
              + atomI.labFrameQuadrupole_ZZ*(atomJ.labFrameQuadrupole_XX*gqzz5+atomJ.labFrameQuadrupole_YY*gqzz8+atomJ.labFrameQuadrupole_ZZ*gqzz10
              +2.0f*(atomJ.labFrameQuadrupole_XY*gqzz6+atomJ.labFrameQuadrupole_XZ*gqzz7+atomJ.labFrameQuadrupole_YZ*gqzz9))
              + 2.0f*(atomI.labFrameQuadrupole_XY*(atomJ.labFrameQuadrupole_XX*gqxy5+atomJ.labFrameQuadrupole_YY*gqxy8+atomJ.labFrameQuadrupole_ZZ*gqxy10
              +2.0f*(atomJ.labFrameQuadrupole_XY*gqxy6+atomJ.labFrameQuadrupole_XZ*gqxy7+atomJ.labFrameQuadrupole_YZ*gqxy9))
              + atomI.labFrameQuadrupole_XZ*(atomJ.labFrameQuadrupole_XX*gqxz5+atomJ.labFrameQuadrupole_YY*gqxz8+atomJ.labFrameQuadrupole_ZZ*gqxz10
              +2.0f*(atomJ.labFrameQuadrupole_XY*gqxz6+atomJ.labFrameQuadrupole_XZ*gqxz7+atomJ.labFrameQuadrupole_YZ*gqxz9))
              + atomI.labFrameQuadrupole_YZ*(atomJ.labFrameQuadrupole_XX*gqyz5+atomJ.labFrameQuadrupole_YY*gqyz8+atomJ.labFrameQuadrupole_ZZ*gqyz10
              +2.0f*(atomJ.labFrameQuadrupole_XY*gqyz6+atomJ.labFrameQuadrupole_XZ*gqyz7+atomJ.labFrameQuadrupole_YZ*gqyz9)));

    ewk = atomI.q*(atomJ.labFrameDipole[0]*gux1+atomJ.labFrameDipole[1]*guy1+atomJ.labFrameDipole[2]*guz1)
                      -atomJ.q*(atomI.labFrameDipole[0]*gc2+atomI.labFrameDipole[1]*gc3+atomI.labFrameDipole[2]*gc4)
                 +atomI.q*(atomJ.labFrameQuadrupole_XX*gqxx1+atomJ.labFrameQuadrupole_YY*gqyy1+atomJ.labFrameQuadrupole_ZZ*gqzz1
              +2.0f*(atomJ.labFrameQuadrupole_XY*gqxy1+atomJ.labFrameQuadrupole_XZ*gqxz1+atomJ.labFrameQuadrupole_YZ*gqyz1))
                 +atomJ.q*(atomI.labFrameQuadrupole_XX*gc5+atomI.labFrameQuadrupole_YY*gc8+atomI.labFrameQuadrupole_ZZ*gc10
              +2.0f*(atomI.labFrameQuadrupole_XY*gc6+atomI.labFrameQuadrupole_XZ*gc7+atomI.labFrameQuadrupole_YZ*gc9))
               - atomI.labFrameDipole[0]*(atomJ.labFrameQuadrupole_XX*gqxx2+atomJ.labFrameQuadrupole_YY*gqyy2+atomJ.labFrameQuadrupole_ZZ*gqzz2
              +2.0f*(atomJ.labFrameQuadrupole_XY*gqxy2+atomJ.labFrameQuadrupole_XZ*gqxz2+atomJ.labFrameQuadrupole_YZ*gqyz2))
               - atomI.labFrameDipole[1]*(atomJ.labFrameQuadrupole_XX*gqxx3+atomJ.labFrameQuadrupole_YY*gqyy3+atomJ.labFrameQuadrupole_ZZ*gqzz3
              +2.0f*(atomJ.labFrameQuadrupole_XY*gqxy3+atomJ.labFrameQuadrupole_XZ*gqxz3+atomJ.labFrameQuadrupole_YZ*gqyz3))
               - atomI.labFrameDipole[2]*(atomJ.labFrameQuadrupole_XX*gqxx4+atomJ.labFrameQuadrupole_YY*gqyy4+atomJ.labFrameQuadrupole_ZZ*gqzz4
              +2.0f*(atomJ.labFrameQuadrupole_XY*gqxy4+atomJ.labFrameQuadrupole_XZ*gqxz4+atomJ.labFrameQuadrupole_YZ*gqyz4))
               + atomJ.labFrameDipole[0]*(atomI.labFrameQuadrupole_XX*gux5+atomI.labFrameQuadrupole_YY*gux8+atomI.labFrameQuadrupole_ZZ*gux10
              +2.0f*(atomI.labFrameQuadrupole_XY*gux6+atomI.labFrameQuadrupole_XZ*gux7+atomI.labFrameQuadrupole_YZ*gux9))
               + atomJ.labFrameDipole[1]*(atomI.labFrameQuadrupole_XX*guy5+atomI.labFrameQuadrupole_YY*guy8+atomI.labFrameQuadrupole_ZZ*guy10
              +2.0f*(atomI.labFrameQuadrupole_XY*guy6+atomI.labFrameQuadrupole_XZ*guy7+atomI.labFrameQuadrupole_YZ*guy9))
               + atomJ.labFrameDipole[2]*(atomI.labFrameQuadrupole_XX*guz5+atomI.labFrameQuadrupole_YY*guz8+atomI.labFrameQuadrupole_ZZ*guz10
              +2.0f*(atomI.labFrameQuadrupole_XY*guz6+atomI.labFrameQuadrupole_XZ*guz7+atomI.labFrameQuadrupole_YZ*guz9))
              + atomI.labFrameQuadrupole_XX*(atomJ.labFrameQuadrupole_XX*gqxx5+atomJ.labFrameQuadrupole_YY*gqyy5+atomJ.labFrameQuadrupole_ZZ*gqzz5
              +2.0f*(atomJ.labFrameQuadrupole_XY*gqxy5+atomJ.labFrameQuadrupole_XZ*gqxz5+atomJ.labFrameQuadrupole_YZ*gqyz5))
              + atomI.labFrameQuadrupole_YY*(atomJ.labFrameQuadrupole_XX*gqxx8+atomJ.labFrameQuadrupole_YY*gqyy8+atomJ.labFrameQuadrupole_ZZ*gqzz8
              +2.0f*(atomJ.labFrameQuadrupole_XY*gqxy8+atomJ.labFrameQuadrupole_XZ*gqxz8+atomJ.labFrameQuadrupole_YZ*gqyz8))
              + atomI.labFrameQuadrupole_ZZ*(atomJ.labFrameQuadrupole_XX*gqxx10+atomJ.labFrameQuadrupole_YY*gqyy10+atomJ.labFrameQuadrupole_ZZ*gqzz10
              +2.0f*(atomJ.labFrameQuadrupole_XY*gqxy10+atomJ.labFrameQuadrupole_XZ*gqxz10+atomJ.labFrameQuadrupole_YZ*gqyz10))
       + 2.0f*(atomI.labFrameQuadrupole_XY*(atomJ.labFrameQuadrupole_XX*gqxx6+atomJ.labFrameQuadrupole_YY*gqyy6+atomJ.labFrameQuadrupole_ZZ*gqzz6
              +2.0f*(atomJ.labFrameQuadrupole_XY*gqxy6+atomJ.labFrameQuadrupole_XZ*gqxz6+atomJ.labFrameQuadrupole_YZ*gqyz6))
              + atomI.labFrameQuadrupole_XZ*(atomJ.labFrameQuadrupole_XX*gqxx7+atomJ.labFrameQuadrupole_YY*gqyy7+atomJ.labFrameQuadrupole_ZZ*gqzz7
              +2.0f*(atomJ.labFrameQuadrupole_XY*gqxy7+atomJ.labFrameQuadrupole_XZ*gqxz7+atomJ.labFrameQuadrupole_YZ*gqyz7))
              + atomI.labFrameQuadrupole_YZ*(atomJ.labFrameQuadrupole_XX*gqxx9+atomJ.labFrameQuadrupole_YY*gqyy9+atomJ.labFrameQuadrupole_ZZ*gqzz9
              +2.0f*(atomJ.labFrameQuadrupole_XY*gqxy9+atomJ.labFrameQuadrupole_XZ*gqxz9+atomJ.labFrameQuadrupole_YZ*gqyz9)));

    desymdx = atomI.q * atomJ.q * gc2 - (atomI.labFrameDipole[0]*(atomJ.labFrameDipole[0]*gux5+atomJ.labFrameDipole[1]*guy5+atomJ.labFrameDipole[2]*guz5)
                              +  atomI.labFrameDipole[1]*(atomJ.labFrameDipole[0]*gux6+atomJ.labFrameDipole[1]*guy6+atomJ.labFrameDipole[2]*guz6)
                              +  atomI.labFrameDipole[2]*(atomJ.labFrameDipole[0]*gux7+atomJ.labFrameDipole[1]*guy7+atomJ.labFrameDipole[2]*guz7));

    dewidx = atomI.q*(atomJ.labFrameDipole[0]*gc5+atomJ.labFrameDipole[1]*gc6+atomJ.labFrameDipole[2]*gc7)
                      -atomJ.q*(atomI.labFrameDipole[0]*gux2+atomI.labFrameDipole[1]*guy2+atomI.labFrameDipole[2]*guz2)
                 +atomI.q*(atomJ.labFrameQuadrupole_XX*gc11+atomJ.labFrameQuadrupole_YY*gc14+atomJ.labFrameQuadrupole_ZZ*gc16
              +2.0f*(atomJ.labFrameQuadrupole_XY*gc12+atomJ.labFrameQuadrupole_XZ*gc13+atomJ.labFrameQuadrupole_YZ*gc15))
                 +atomJ.q*(atomI.labFrameQuadrupole_XX*gqxx2+atomI.labFrameQuadrupole_YY*gqyy2+atomI.labFrameQuadrupole_ZZ*gqzz2
              +2.0f*(atomI.labFrameQuadrupole_XY*gqxy2+atomI.labFrameQuadrupole_XZ*gqxz2+atomI.labFrameQuadrupole_YZ*gqyz2))
               - atomI.labFrameDipole[0]*(atomJ.labFrameQuadrupole_XX*gux11+atomJ.labFrameQuadrupole_YY*gux14+atomJ.labFrameQuadrupole_ZZ*gux16
              +2.0f*(atomJ.labFrameQuadrupole_XY*gux12+atomJ.labFrameQuadrupole_XZ*gux13+atomJ.labFrameQuadrupole_YZ*gux15))
               - atomI.labFrameDipole[1]*(atomJ.labFrameQuadrupole_XX*guy11+atomJ.labFrameQuadrupole_YY*guy14+atomJ.labFrameQuadrupole_ZZ*guy16
              +2.0f*(atomJ.labFrameQuadrupole_XY*guy12+atomJ.labFrameQuadrupole_XZ*guy13+atomJ.labFrameQuadrupole_YZ*guy15))
               - atomI.labFrameDipole[2]*(atomJ.labFrameQuadrupole_XX*guz11+atomJ.labFrameQuadrupole_YY*guz14+atomJ.labFrameQuadrupole_ZZ*guz16
              +2.0f*(atomJ.labFrameQuadrupole_XY*guz12+atomJ.labFrameQuadrupole_XZ*guz13+atomJ.labFrameQuadrupole_YZ*guz15))
               + atomJ.labFrameDipole[0]*(atomI.labFrameQuadrupole_XX*gqxx5+atomI.labFrameQuadrupole_YY*gqyy5+atomI.labFrameQuadrupole_ZZ*gqzz5
              +2.0f*(atomI.labFrameQuadrupole_XY*gqxy5+atomI.labFrameQuadrupole_XZ*gqxz5+atomI.labFrameQuadrupole_YZ*gqyz5))
               + atomJ.labFrameDipole[1]*(atomI.labFrameQuadrupole_XX*gqxx6+atomI.labFrameQuadrupole_YY*gqyy6+atomI.labFrameQuadrupole_ZZ*gqzz6
              +2.0f*(atomI.labFrameQuadrupole_XY*gqxy6+atomI.labFrameQuadrupole_XZ*gqxz6+atomI.labFrameQuadrupole_YZ*gqyz6))
               + atomJ.labFrameDipole[2]*(atomI.labFrameQuadrupole_XX*gqxx7+atomI.labFrameQuadrupole_YY*gqyy7+atomI.labFrameQuadrupole_ZZ*gqzz7
              +2.0f*(atomI.labFrameQuadrupole_XY*gqxy7+atomI.labFrameQuadrupole_XZ*gqxz7+atomI.labFrameQuadrupole_YZ*gqyz7))
              + atomI.labFrameQuadrupole_XX*(atomJ.labFrameQuadrupole_XX*gqxx11+atomJ.labFrameQuadrupole_YY*gqxx14+atomJ.labFrameQuadrupole_ZZ*gqxx16
              +2.0f*(atomJ.labFrameQuadrupole_XY*gqxx12+atomJ.labFrameQuadrupole_XZ*gqxx13+atomJ.labFrameQuadrupole_YZ*gqxx15))
              + atomI.labFrameQuadrupole_YY*(atomJ.labFrameQuadrupole_XX*gqyy11+atomJ.labFrameQuadrupole_YY*gqyy14+atomJ.labFrameQuadrupole_ZZ*gqyy16
              +2.0f*(atomJ.labFrameQuadrupole_XY*gqyy12+atomJ.labFrameQuadrupole_XZ*gqyy13+atomJ.labFrameQuadrupole_YZ*gqyy15))
              + atomI.labFrameQuadrupole_ZZ*(atomJ.labFrameQuadrupole_XX*gqzz11+atomJ.labFrameQuadrupole_YY*gqzz14+atomJ.labFrameQuadrupole_ZZ*gqzz16
              +2.0f*(atomJ.labFrameQuadrupole_XY*gqzz12+atomJ.labFrameQuadrupole_XZ*gqzz13+atomJ.labFrameQuadrupole_YZ*gqzz15))
       + 2.0f*(atomI.labFrameQuadrupole_XY*(atomJ.labFrameQuadrupole_XX*gqxy11+atomJ.labFrameQuadrupole_YY*gqxy14+atomJ.labFrameQuadrupole_ZZ*gqxy16
              +2.0f*(atomJ.labFrameQuadrupole_XY*gqxy12+atomJ.labFrameQuadrupole_XZ*gqxy13+atomJ.labFrameQuadrupole_YZ*gqxy15))
              + atomI.labFrameQuadrupole_XZ*(atomJ.labFrameQuadrupole_XX*gqxz11+atomJ.labFrameQuadrupole_YY*gqxz14+atomJ.labFrameQuadrupole_ZZ*gqxz16
              +2.0f*(atomJ.labFrameQuadrupole_XY*gqxz12+atomJ.labFrameQuadrupole_XZ*gqxz13+atomJ.labFrameQuadrupole_YZ*gqxz15))
              + atomI.labFrameQuadrupole_YZ*(atomJ.labFrameQuadrupole_XX*gqyz11+atomJ.labFrameQuadrupole_YY*gqyz14+atomJ.labFrameQuadrupole_ZZ*gqyz16
              +2.0f*(atomJ.labFrameQuadrupole_XY*gqyz12+atomJ.labFrameQuadrupole_XZ*gqyz13+atomJ.labFrameQuadrupole_YZ*gqyz15)));

    dewkdx = atomI.q*(atomJ.labFrameDipole[0]*gux2+atomJ.labFrameDipole[1]*guy2+atomJ.labFrameDipole[2]*guz2)
                      -atomJ.q*(atomI.labFrameDipole[0]*gc5+atomI.labFrameDipole[1]*gc6+atomI.labFrameDipole[2]*gc7)
                 +atomI.q*(atomJ.labFrameQuadrupole_XX*gqxx2+atomJ.labFrameQuadrupole_YY*gqyy2+atomJ.labFrameQuadrupole_ZZ*gqzz2
              +2.0f*(atomJ.labFrameQuadrupole_XY*gqxy2+atomJ.labFrameQuadrupole_XZ*gqxz2+atomJ.labFrameQuadrupole_YZ*gqyz2))
                 +atomJ.q*(atomI.labFrameQuadrupole_XX*gc11+atomI.labFrameQuadrupole_YY*gc14+atomI.labFrameQuadrupole_ZZ*gc16
              +2.0f*(atomI.labFrameQuadrupole_XY*gc12+atomI.labFrameQuadrupole_XZ*gc13+atomI.labFrameQuadrupole_YZ*gc15))
               - atomI.labFrameDipole[0]*(atomJ.labFrameQuadrupole_XX*gqxx5+atomJ.labFrameQuadrupole_YY*gqyy5+atomJ.labFrameQuadrupole_ZZ*gqzz5
              +2.0f*(atomJ.labFrameQuadrupole_XY*gqxy5+atomJ.labFrameQuadrupole_XZ*gqxz5+atomJ.labFrameQuadrupole_YZ*gqyz5))
               - atomI.labFrameDipole[1]*(atomJ.labFrameQuadrupole_XX*gqxx6+atomJ.labFrameQuadrupole_YY*gqyy6+atomJ.labFrameQuadrupole_ZZ*gqzz6
              +2.0f*(atomJ.labFrameQuadrupole_XY*gqxy6+atomJ.labFrameQuadrupole_XZ*gqxz6+atomJ.labFrameQuadrupole_YZ*gqyz6))
               - atomI.labFrameDipole[2]*(atomJ.labFrameQuadrupole_XX*gqxx7+atomJ.labFrameQuadrupole_YY*gqyy7+atomJ.labFrameQuadrupole_ZZ*gqzz7
              +2.0f*(atomJ.labFrameQuadrupole_XY*gqxy7+atomJ.labFrameQuadrupole_XZ*gqxz7+atomJ.labFrameQuadrupole_YZ*gqyz7))
               + atomJ.labFrameDipole[0]*(atomI.labFrameQuadrupole_XX*gux11+atomI.labFrameQuadrupole_YY*gux14+atomI.labFrameQuadrupole_ZZ*gux16
              +2.0f*(atomI.labFrameQuadrupole_XY*gux12+atomI.labFrameQuadrupole_XZ*gux13+atomI.labFrameQuadrupole_YZ*gux15))
               + atomJ.labFrameDipole[1]*(atomI.labFrameQuadrupole_XX*guy11+atomI.labFrameQuadrupole_YY*guy14+atomI.labFrameQuadrupole_ZZ*guy16
              +2.0f*(atomI.labFrameQuadrupole_XY*guy12+atomI.labFrameQuadrupole_XZ*guy13+atomI.labFrameQuadrupole_YZ*guy15))
               + atomJ.labFrameDipole[2]*(atomI.labFrameQuadrupole_XX*guz11+atomI.labFrameQuadrupole_YY*guz14+atomI.labFrameQuadrupole_ZZ*guz16
              +2.0f*(atomI.labFrameQuadrupole_XY*guz12+atomI.labFrameQuadrupole_XZ*guz13+atomI.labFrameQuadrupole_YZ*guz15))
              + atomI.labFrameQuadrupole_XX*(atomJ.labFrameQuadrupole_XX*gqxx11+atomJ.labFrameQuadrupole_YY*gqyy11+atomJ.labFrameQuadrupole_ZZ*gqzz11
              +2.0f*(atomJ.labFrameQuadrupole_XY*gqxy11+atomJ.labFrameQuadrupole_XZ*gqxz11+atomJ.labFrameQuadrupole_YZ*gqyz11))
              + atomI.labFrameQuadrupole_YY*(atomJ.labFrameQuadrupole_XX*gqxx14+atomJ.labFrameQuadrupole_YY*gqyy14+atomJ.labFrameQuadrupole_ZZ*gqzz14
              +2.0f*(atomJ.labFrameQuadrupole_XY*gqxy14+atomJ.labFrameQuadrupole_XZ*gqxz14+atomJ.labFrameQuadrupole_YZ*gqyz14))
              + atomI.labFrameQuadrupole_ZZ*(atomJ.labFrameQuadrupole_XX*gqxx16+atomJ.labFrameQuadrupole_YY*gqyy16+atomJ.labFrameQuadrupole_ZZ*gqzz16
              +2.0f*(atomJ.labFrameQuadrupole_XY*gqxy16+atomJ.labFrameQuadrupole_XZ*gqxz16+atomJ.labFrameQuadrupole_YZ*gqyz16))
       + 2.0f*(atomI.labFrameQuadrupole_XY*(atomJ.labFrameQuadrupole_XX*gqxx12+atomJ.labFrameQuadrupole_YY*gqyy12+atomJ.labFrameQuadrupole_ZZ*gqzz12
              +2.0f*(atomJ.labFrameQuadrupole_XY*gqxy12+atomJ.labFrameQuadrupole_XZ*gqxz12+atomJ.labFrameQuadrupole_YZ*gqyz12))
              + atomI.labFrameQuadrupole_XZ*(atomJ.labFrameQuadrupole_XX*gqxx13+atomJ.labFrameQuadrupole_YY*gqyy13+atomJ.labFrameQuadrupole_ZZ*gqzz13
              +2.0f*(atomJ.labFrameQuadrupole_XY*gqxy13+atomJ.labFrameQuadrupole_XZ*gqxz13+atomJ.labFrameQuadrupole_YZ*gqyz13))
              + atomI.labFrameQuadrupole_YZ*(atomJ.labFrameQuadrupole_XX*gqxx15+atomJ.labFrameQuadrupole_YY*gqyy15+atomJ.labFrameQuadrupole_ZZ*gqzz15
              +2.0f*(atomJ.labFrameQuadrupole_XY*gqxy15+atomJ.labFrameQuadrupole_XZ*gqxz15+atomJ.labFrameQuadrupole_YZ*gqyz15)));

    dedx = desymdx + 0.5f*(dewidx + dewkdx);

    desymdy = atomI.q * atomJ.q * gc3
                           - (atomI.labFrameDipole[0]*(atomJ.labFrameDipole[0]*gux6+atomJ.labFrameDipole[1]*guy6+atomJ.labFrameDipole[2]*guz6)
                             +atomI.labFrameDipole[1]*(atomJ.labFrameDipole[0]*gux8+atomJ.labFrameDipole[1]*guy8+atomJ.labFrameDipole[2]*guz8)
                             +atomI.labFrameDipole[2]*(atomJ.labFrameDipole[0]*gux9+atomJ.labFrameDipole[1]*guy9+atomJ.labFrameDipole[2]*guz9));

    dewidy = atomI.q*(atomJ.labFrameDipole[0]*gc6+atomJ.labFrameDipole[1]*gc8+atomJ.labFrameDipole[2]*gc9)
                      -atomJ.q*(atomI.labFrameDipole[0]*gux3+atomI.labFrameDipole[1]*guy3+atomI.labFrameDipole[2]*guz3)
                 +atomI.q*(atomJ.labFrameQuadrupole_XX*gc12+atomJ.labFrameQuadrupole_YY*gc17+atomJ.labFrameQuadrupole_ZZ*gc19
              +2.0f*(atomJ.labFrameQuadrupole_XY*gc14+atomJ.labFrameQuadrupole_XZ*gc15+atomJ.labFrameQuadrupole_YZ*gc18))
                 +atomJ.q*(atomI.labFrameQuadrupole_XX*gqxx3+atomI.labFrameQuadrupole_YY*gqyy3+atomI.labFrameQuadrupole_ZZ*gqzz3
              +2.0f*(atomI.labFrameQuadrupole_XY*gqxy3+atomI.labFrameQuadrupole_XZ*gqxz3+atomI.labFrameQuadrupole_YZ*gqyz3))
               - atomI.labFrameDipole[0]*(atomJ.labFrameQuadrupole_XX*gux12+atomJ.labFrameQuadrupole_YY*gux17+atomJ.labFrameQuadrupole_ZZ*gux19
              +2.0f*(atomJ.labFrameQuadrupole_XY*gux14+atomJ.labFrameQuadrupole_XZ*gux15+atomJ.labFrameQuadrupole_YZ*gux18))
               - atomI.labFrameDipole[1]*(atomJ.labFrameQuadrupole_XX*guy12+atomJ.labFrameQuadrupole_YY*guy17+atomJ.labFrameQuadrupole_ZZ*guy19
              +2.0f*(atomJ.labFrameQuadrupole_XY*guy14+atomJ.labFrameQuadrupole_XZ*guy15+atomJ.labFrameQuadrupole_YZ*guy18))
               - atomI.labFrameDipole[2]*(atomJ.labFrameQuadrupole_XX*guz12+atomJ.labFrameQuadrupole_YY*guz17+atomJ.labFrameQuadrupole_ZZ*guz19
              +2.0f*(atomJ.labFrameQuadrupole_XY*guz14+atomJ.labFrameQuadrupole_XZ*guz15+atomJ.labFrameQuadrupole_YZ*guz18))
               + atomJ.labFrameDipole[0]*(atomI.labFrameQuadrupole_XX*gqxx6+atomI.labFrameQuadrupole_YY*gqyy6+atomI.labFrameQuadrupole_ZZ*gqzz6
              +2.0f*(atomI.labFrameQuadrupole_XY*gqxy6+atomI.labFrameQuadrupole_XZ*gqxz6+atomI.labFrameQuadrupole_YZ*gqyz6))
               + atomJ.labFrameDipole[1]*(atomI.labFrameQuadrupole_XX*gqxx8+atomI.labFrameQuadrupole_YY*gqyy8+atomI.labFrameQuadrupole_ZZ*gqzz8
              +2.0f*(atomI.labFrameQuadrupole_XY*gqxy8+atomI.labFrameQuadrupole_XZ*gqxz8+atomI.labFrameQuadrupole_YZ*gqyz8))
               + atomJ.labFrameDipole[2]*(atomI.labFrameQuadrupole_XX*gqxx9+atomI.labFrameQuadrupole_YY*gqyy9+atomI.labFrameQuadrupole_ZZ*gqzz9
              +2.0f*(atomI.labFrameQuadrupole_XY*gqxy9+atomI.labFrameQuadrupole_XZ*gqxz9+atomI.labFrameQuadrupole_YZ*gqyz9))
              + atomI.labFrameQuadrupole_XX*(atomJ.labFrameQuadrupole_XX*gqxx12+atomJ.labFrameQuadrupole_YY*gqxx17+atomJ.labFrameQuadrupole_ZZ*gqxx19
              +2.0f*(atomJ.labFrameQuadrupole_XY*gqxx14+atomJ.labFrameQuadrupole_XZ*gqxx15+atomJ.labFrameQuadrupole_YZ*gqxx18))
              + atomI.labFrameQuadrupole_YY*(atomJ.labFrameQuadrupole_XX*gqyy12+atomJ.labFrameQuadrupole_YY*gqyy17+atomJ.labFrameQuadrupole_ZZ*gqyy19
              +2.0f*(atomJ.labFrameQuadrupole_XY*gqyy14+atomJ.labFrameQuadrupole_XZ*gqyy15+atomJ.labFrameQuadrupole_YZ*gqyy18))
              + atomI.labFrameQuadrupole_ZZ*(atomJ.labFrameQuadrupole_XX*gqzz12+atomJ.labFrameQuadrupole_YY*gqzz17+atomJ.labFrameQuadrupole_ZZ*gqzz19
              +2.0f*(atomJ.labFrameQuadrupole_XY*gqzz14+atomJ.labFrameQuadrupole_XZ*gqzz15+atomJ.labFrameQuadrupole_YZ*gqzz18))
       + 2.0f*(atomI.labFrameQuadrupole_XY*(atomJ.labFrameQuadrupole_XX*gqxy12+atomJ.labFrameQuadrupole_YY*gqxy17+atomJ.labFrameQuadrupole_ZZ*gqxy19
              +2.0f*(atomJ.labFrameQuadrupole_XY*gqxy14+atomJ.labFrameQuadrupole_XZ*gqxy15+atomJ.labFrameQuadrupole_YZ*gqxy18))
              + atomI.labFrameQuadrupole_XZ*(atomJ.labFrameQuadrupole_XX*gqxz12+atomJ.labFrameQuadrupole_YY*gqxz17+atomJ.labFrameQuadrupole_ZZ*gqxz19
              +2.0f*(atomJ.labFrameQuadrupole_XY*gqxz14+atomJ.labFrameQuadrupole_XZ*gqxz15+atomJ.labFrameQuadrupole_YZ*gqxz18))
              + atomI.labFrameQuadrupole_YZ*(atomJ.labFrameQuadrupole_XX*gqyz12+atomJ.labFrameQuadrupole_YY*gqyz17+atomJ.labFrameQuadrupole_ZZ*gqyz19
              +2.0f*(atomJ.labFrameQuadrupole_XY*gqyz14+atomJ.labFrameQuadrupole_XZ*gqyz15+atomJ.labFrameQuadrupole_YZ*gqyz18)));

    dewkdy = atomI.q*(atomJ.labFrameDipole[0]*gux3+atomJ.labFrameDipole[1]*guy3+atomJ.labFrameDipole[2]*guz3)
                      -atomJ.q*(atomI.labFrameDipole[0]*gc6+atomI.labFrameDipole[1]*gc8+atomI.labFrameDipole[2]*gc9)
                 +atomI.q*(atomJ.labFrameQuadrupole_XX*gqxx3+atomJ.labFrameQuadrupole_YY*gqyy3+atomJ.labFrameQuadrupole_ZZ*gqzz3
              +2.0f*(atomJ.labFrameQuadrupole_XY*gqxy3+atomJ.labFrameQuadrupole_XZ*gqxz3+atomJ.labFrameQuadrupole_YZ*gqyz3))
                 +atomJ.q*(atomI.labFrameQuadrupole_XX*gc12+atomI.labFrameQuadrupole_YY*gc17+atomI.labFrameQuadrupole_ZZ*gc19
              +2.0f*(atomI.labFrameQuadrupole_XY*gc14+atomI.labFrameQuadrupole_XZ*gc15+atomI.labFrameQuadrupole_YZ*gc18))
               - atomI.labFrameDipole[0]*(atomJ.labFrameQuadrupole_XX*gqxx6+atomJ.labFrameQuadrupole_YY*gqyy6+atomJ.labFrameQuadrupole_ZZ*gqzz6
              +2.0f*(atomJ.labFrameQuadrupole_XY*gqxy6+atomJ.labFrameQuadrupole_XZ*gqxz6+atomJ.labFrameQuadrupole_YZ*gqyz6))
               - atomI.labFrameDipole[1]*(atomJ.labFrameQuadrupole_XX*gqxx8+atomJ.labFrameQuadrupole_YY*gqyy8+atomJ.labFrameQuadrupole_ZZ*gqzz8
              +2.0f*(atomJ.labFrameQuadrupole_XY*gqxy8+atomJ.labFrameQuadrupole_XZ*gqxz8+atomJ.labFrameQuadrupole_YZ*gqyz8))
               - atomI.labFrameDipole[2]*(atomJ.labFrameQuadrupole_XX*gqxx9+atomJ.labFrameQuadrupole_YY*gqyy9+atomJ.labFrameQuadrupole_ZZ*gqzz9
              +2.0f*(atomJ.labFrameQuadrupole_XY*gqxy9+atomJ.labFrameQuadrupole_XZ*gqxz9+atomJ.labFrameQuadrupole_YZ*gqyz9))
               + atomJ.labFrameDipole[0]*(atomI.labFrameQuadrupole_XX*gux12+atomI.labFrameQuadrupole_YY*gux17+atomI.labFrameQuadrupole_ZZ*gux19
              +2.0f*(atomI.labFrameQuadrupole_XY*gux14+atomI.labFrameQuadrupole_XZ*gux15+atomI.labFrameQuadrupole_YZ*gux18))
               + atomJ.labFrameDipole[1]*(atomI.labFrameQuadrupole_XX*guy12+atomI.labFrameQuadrupole_YY*guy17+atomI.labFrameQuadrupole_ZZ*guy19
              +2.0f*(atomI.labFrameQuadrupole_XY*guy14+atomI.labFrameQuadrupole_XZ*guy15+atomI.labFrameQuadrupole_YZ*guy18))
               + atomJ.labFrameDipole[2]*(atomI.labFrameQuadrupole_XX*guz12+atomI.labFrameQuadrupole_YY*guz17+atomI.labFrameQuadrupole_ZZ*guz19
              +2.0f*(atomI.labFrameQuadrupole_XY*guz14+atomI.labFrameQuadrupole_XZ*guz15+atomI.labFrameQuadrupole_YZ*guz18))
              + atomI.labFrameQuadrupole_XX*(atomJ.labFrameQuadrupole_XX*gqxx12+atomJ.labFrameQuadrupole_YY*gqyy12+atomJ.labFrameQuadrupole_ZZ*gqzz12
              +2.0f*(atomJ.labFrameQuadrupole_XY*gqxy12+atomJ.labFrameQuadrupole_XZ*gqxz12+atomJ.labFrameQuadrupole_YZ*gqyz12))
              + atomI.labFrameQuadrupole_YY*(atomJ.labFrameQuadrupole_XX*gqxx17+atomJ.labFrameQuadrupole_YY*gqyy17+atomJ.labFrameQuadrupole_ZZ*gqzz17
              +2.0f*(atomJ.labFrameQuadrupole_XY*gqxy17+atomJ.labFrameQuadrupole_XZ*gqxz17+atomJ.labFrameQuadrupole_YZ*gqyz17))
              + atomI.labFrameQuadrupole_ZZ*(atomJ.labFrameQuadrupole_XX*gqxx19+atomJ.labFrameQuadrupole_YY*gqyy19+atomJ.labFrameQuadrupole_ZZ*gqzz19
              +2.0f*(atomJ.labFrameQuadrupole_XY*gqxy19+atomJ.labFrameQuadrupole_XZ*gqxz19+atomJ.labFrameQuadrupole_YZ*gqyz19))
       + 2.0f*(atomI.labFrameQuadrupole_XY*(atomJ.labFrameQuadrupole_XX*gqxx14+atomJ.labFrameQuadrupole_YY*gqyy14+atomJ.labFrameQuadrupole_ZZ*gqzz14
              +2.0f*(atomJ.labFrameQuadrupole_XY*gqxy14+atomJ.labFrameQuadrupole_XZ*gqxz14+atomJ.labFrameQuadrupole_YZ*gqyz14))
              + atomI.labFrameQuadrupole_XZ*(atomJ.labFrameQuadrupole_XX*gqxx15+atomJ.labFrameQuadrupole_YY*gqyy15+atomJ.labFrameQuadrupole_ZZ*gqzz15
              +2.0f*(atomJ.labFrameQuadrupole_XY*gqxy15+atomJ.labFrameQuadrupole_XZ*gqxz15+atomJ.labFrameQuadrupole_YZ*gqyz15))
              + atomI.labFrameQuadrupole_YZ*(atomJ.labFrameQuadrupole_XX*gqxx18+atomJ.labFrameQuadrupole_YY*gqyy18+atomJ.labFrameQuadrupole_ZZ*gqzz18
              +2.0f*(atomJ.labFrameQuadrupole_XY*gqxy18+atomJ.labFrameQuadrupole_XZ*gqxz18+atomJ.labFrameQuadrupole_YZ*gqyz18)));

    dedy = desymdy + 0.5f*(dewidy + dewkdy);

    desymdz = atomI.q * atomJ.q * gc4
                           - (atomI.labFrameDipole[0]*(atomJ.labFrameDipole[0]*gux7+atomJ.labFrameDipole[1]*guy7+atomJ.labFrameDipole[2]*guz7)
                             +atomI.labFrameDipole[1]*(atomJ.labFrameDipole[0]*gux9+atomJ.labFrameDipole[1]*guy9+atomJ.labFrameDipole[2]*guz9)
                             +atomI.labFrameDipole[2]*(atomJ.labFrameDipole[0]*gux10+atomJ.labFrameDipole[1]*guy10+atomJ.labFrameDipole[2]*guz10));

    dewidz = atomI.q*(atomJ.labFrameDipole[0]*gc7+atomJ.labFrameDipole[1]*gc9+atomJ.labFrameDipole[2]*gc10)
                      -atomJ.q*(atomI.labFrameDipole[0]*gux4+atomI.labFrameDipole[1]*guy4+atomI.labFrameDipole[2]*guz4)
                 +atomI.q*(atomJ.labFrameQuadrupole_XX*gc13+atomJ.labFrameQuadrupole_YY*gc18+atomJ.labFrameQuadrupole_ZZ*gc20
              +2.0f*(atomJ.labFrameQuadrupole_XY*gc15+atomJ.labFrameQuadrupole_XZ*gc16+atomJ.labFrameQuadrupole_YZ*gc19))
                 +atomJ.q*(atomI.labFrameQuadrupole_XX*gqxx4+atomI.labFrameQuadrupole_YY*gqyy4+atomI.labFrameQuadrupole_ZZ*gqzz4
              +2.0f*(atomI.labFrameQuadrupole_XY*gqxy4+atomI.labFrameQuadrupole_XZ*gqxz4+atomI.labFrameQuadrupole_YZ*gqyz4))
               - atomI.labFrameDipole[0]*(atomJ.labFrameQuadrupole_XX*gux13+atomJ.labFrameQuadrupole_YY*gux18+atomJ.labFrameQuadrupole_ZZ*gux20
              +2.0f*(atomJ.labFrameQuadrupole_XY*gux15+atomJ.labFrameQuadrupole_XZ*gux16+atomJ.labFrameQuadrupole_YZ*gux19))
               - atomI.labFrameDipole[1]*(atomJ.labFrameQuadrupole_XX*guy13+atomJ.labFrameQuadrupole_YY*guy18+atomJ.labFrameQuadrupole_ZZ*guy20
              +2.0f*(atomJ.labFrameQuadrupole_XY*guy15+atomJ.labFrameQuadrupole_XZ*guy16+atomJ.labFrameQuadrupole_YZ*guy19))
               - atomI.labFrameDipole[2]*(atomJ.labFrameQuadrupole_XX*guz13+atomJ.labFrameQuadrupole_YY*guz18+atomJ.labFrameQuadrupole_ZZ*guz20
              +2.0f*(atomJ.labFrameQuadrupole_XY*guz15+atomJ.labFrameQuadrupole_XZ*guz16+atomJ.labFrameQuadrupole_YZ*guz19))
               + atomJ.labFrameDipole[0]*(atomI.labFrameQuadrupole_XX*gqxx7+atomI.labFrameQuadrupole_YY*gqyy7+atomI.labFrameQuadrupole_ZZ*gqzz7
              +2.0f*(atomI.labFrameQuadrupole_XY*gqxy7+atomI.labFrameQuadrupole_XZ*gqxz7+atomI.labFrameQuadrupole_YZ*gqyz7))
               + atomJ.labFrameDipole[1]*(atomI.labFrameQuadrupole_XX*gqxx9+atomI.labFrameQuadrupole_YY*gqyy9+atomI.labFrameQuadrupole_ZZ*gqzz9
              +2.0f*(atomI.labFrameQuadrupole_XY*gqxy9+atomI.labFrameQuadrupole_XZ*gqxz9+atomI.labFrameQuadrupole_YZ*gqyz9))
               + atomJ.labFrameDipole[2]*(atomI.labFrameQuadrupole_XX*gqxx10+atomI.labFrameQuadrupole_YY*gqyy10+atomI.labFrameQuadrupole_ZZ*gqzz10
              +2.0f*(atomI.labFrameQuadrupole_XY*gqxy10+atomI.labFrameQuadrupole_XZ*gqxz10+atomI.labFrameQuadrupole_YZ*gqyz10))
              + atomI.labFrameQuadrupole_XX*(atomJ.labFrameQuadrupole_XX*gqxx13+atomJ.labFrameQuadrupole_YY*gqxx18+atomJ.labFrameQuadrupole_ZZ*gqxx20
              +2.0f*(atomJ.labFrameQuadrupole_XY*gqxx15+atomJ.labFrameQuadrupole_XZ*gqxx16+atomJ.labFrameQuadrupole_YZ*gqxx19))
              + atomI.labFrameQuadrupole_YY*(atomJ.labFrameQuadrupole_XX*gqyy13+atomJ.labFrameQuadrupole_YY*gqyy18+atomJ.labFrameQuadrupole_ZZ*gqyy20
              +2.0f*(atomJ.labFrameQuadrupole_XY*gqyy15+atomJ.labFrameQuadrupole_XZ*gqyy16+atomJ.labFrameQuadrupole_YZ*gqyy19))
              + atomI.labFrameQuadrupole_ZZ*(atomJ.labFrameQuadrupole_XX*gqzz13+atomJ.labFrameQuadrupole_YY*gqzz18+atomJ.labFrameQuadrupole_ZZ*gqzz20
              +2.0f*(atomJ.labFrameQuadrupole_XY*gqzz15+atomJ.labFrameQuadrupole_XZ*gqzz16+atomJ.labFrameQuadrupole_YZ*gqzz19))
       + 2.0f*(atomI.labFrameQuadrupole_XY*(atomJ.labFrameQuadrupole_XX*gqxy13+atomJ.labFrameQuadrupole_YY*gqxy18+atomJ.labFrameQuadrupole_ZZ*gqxy20
              +2.0f*(atomJ.labFrameQuadrupole_XY*gqxy15+atomJ.labFrameQuadrupole_XZ*gqxy16+atomJ.labFrameQuadrupole_YZ*gqxy19))
              + atomI.labFrameQuadrupole_XZ*(atomJ.labFrameQuadrupole_XX*gqxz13+atomJ.labFrameQuadrupole_YY*gqxz18+atomJ.labFrameQuadrupole_ZZ*gqxz20
              +2.0f*(atomJ.labFrameQuadrupole_XY*gqxz15+atomJ.labFrameQuadrupole_XZ*gqxz16+atomJ.labFrameQuadrupole_YZ*gqxz19))
              + atomI.labFrameQuadrupole_YZ*(atomJ.labFrameQuadrupole_XX*gqyz13+atomJ.labFrameQuadrupole_YY*gqyz18+atomJ.labFrameQuadrupole_ZZ*gqyz20
              +2.0f*(atomJ.labFrameQuadrupole_XY*gqyz15+atomJ.labFrameQuadrupole_XZ*gqyz16+atomJ.labFrameQuadrupole_YZ*gqyz19)));

    dewkdz = atomI.q*(atomJ.labFrameDipole[0]*gux4+atomJ.labFrameDipole[1]*guy4+atomJ.labFrameDipole[2]*guz4)
                      -atomJ.q*(atomI.labFrameDipole[0]*gc7+atomI.labFrameDipole[1]*gc9+atomI.labFrameDipole[2]*gc10)
                 +atomI.q*(atomJ.labFrameQuadrupole_XX*gqxx4+atomJ.labFrameQuadrupole_YY*gqyy4+atomJ.labFrameQuadrupole_ZZ*gqzz4
              +2.0f*(atomJ.labFrameQuadrupole_XY*gqxy4+atomJ.labFrameQuadrupole_XZ*gqxz4+atomJ.labFrameQuadrupole_YZ*gqyz4))
                 +atomJ.q*(atomI.labFrameQuadrupole_XX*gc13+atomI.labFrameQuadrupole_YY*gc18+atomI.labFrameQuadrupole_ZZ*gc20
              +2.0f*(atomI.labFrameQuadrupole_XY*gc15+atomI.labFrameQuadrupole_XZ*gc16+atomI.labFrameQuadrupole_YZ*gc19))
               - atomI.labFrameDipole[0]*(atomJ.labFrameQuadrupole_XX*gqxx7+atomJ.labFrameQuadrupole_YY*gqyy7+atomJ.labFrameQuadrupole_ZZ*gqzz7
              +2.0f*(atomJ.labFrameQuadrupole_XY*gqxy7+atomJ.labFrameQuadrupole_XZ*gqxz7+atomJ.labFrameQuadrupole_YZ*gqyz7))
               - atomI.labFrameDipole[1]*(atomJ.labFrameQuadrupole_XX*gqxx9+atomJ.labFrameQuadrupole_YY*gqyy9+atomJ.labFrameQuadrupole_ZZ*gqzz9
              +2.0f*(atomJ.labFrameQuadrupole_XY*gqxy9+atomJ.labFrameQuadrupole_XZ*gqxz9+atomJ.labFrameQuadrupole_YZ*gqyz9))
               - atomI.labFrameDipole[2]*(atomJ.labFrameQuadrupole_XX*gqxx10+atomJ.labFrameQuadrupole_YY*gqyy10+atomJ.labFrameQuadrupole_ZZ*gqzz10
              +2.0f*(atomJ.labFrameQuadrupole_XY*gqxy10+atomJ.labFrameQuadrupole_XZ*gqxz10+atomJ.labFrameQuadrupole_YZ*gqyz10))
               + atomJ.labFrameDipole[0]*(atomI.labFrameQuadrupole_XX*gux13+atomI.labFrameQuadrupole_YY*gux18+atomI.labFrameQuadrupole_ZZ*gux20
              +2.0f*(atomI.labFrameQuadrupole_XY*gux15+atomI.labFrameQuadrupole_XZ*gux16+atomI.labFrameQuadrupole_YZ*gux19))
               + atomJ.labFrameDipole[1]*(atomI.labFrameQuadrupole_XX*guy13+atomI.labFrameQuadrupole_YY*guy18+atomI.labFrameQuadrupole_ZZ*guy20
              +2.0f*(atomI.labFrameQuadrupole_XY*guy15+atomI.labFrameQuadrupole_XZ*guy16+atomI.labFrameQuadrupole_YZ*guy19))
               + atomJ.labFrameDipole[2]*(atomI.labFrameQuadrupole_XX*guz13+atomI.labFrameQuadrupole_YY*guz18+atomI.labFrameQuadrupole_ZZ*guz20
              +2.0f*(atomI.labFrameQuadrupole_XY*guz15+atomI.labFrameQuadrupole_XZ*guz16+atomI.labFrameQuadrupole_YZ*guz19))
              + atomI.labFrameQuadrupole_XX*(atomJ.labFrameQuadrupole_XX*gqxx13+atomJ.labFrameQuadrupole_YY*gqyy13+atomJ.labFrameQuadrupole_ZZ*gqzz13
              +2.0f*(atomJ.labFrameQuadrupole_XY*gqxy13+atomJ.labFrameQuadrupole_XZ*gqxz13+atomJ.labFrameQuadrupole_YZ*gqyz13))
              + atomI.labFrameQuadrupole_YY*(atomJ.labFrameQuadrupole_XX*gqxx18+atomJ.labFrameQuadrupole_YY*gqyy18+atomJ.labFrameQuadrupole_ZZ*gqzz18
              +2.0f*(atomJ.labFrameQuadrupole_XY*gqxy18+atomJ.labFrameQuadrupole_XZ*gqxz18+atomJ.labFrameQuadrupole_YZ*gqyz18))
              + atomI.labFrameQuadrupole_ZZ*(atomJ.labFrameQuadrupole_XX*gqxx20+atomJ.labFrameQuadrupole_YY*gqyy20+atomJ.labFrameQuadrupole_ZZ*gqzz20
              +2.0f*(atomJ.labFrameQuadrupole_XY*gqxy20+atomJ.labFrameQuadrupole_XZ*gqxz20+atomJ.labFrameQuadrupole_YZ*gqyz20))
       + 2.0f*(atomI.labFrameQuadrupole_XY*(atomJ.labFrameQuadrupole_XX*gqxx15+atomJ.labFrameQuadrupole_YY*gqyy15+atomJ.labFrameQuadrupole_ZZ*gqzz15
              +2.0f*(atomJ.labFrameQuadrupole_XY*gqxy15+atomJ.labFrameQuadrupole_XZ*gqxz15+atomJ.labFrameQuadrupole_YZ*gqyz15))
              + atomI.labFrameQuadrupole_XZ*(atomJ.labFrameQuadrupole_XX*gqxx16+atomJ.labFrameQuadrupole_YY*gqyy16+atomJ.labFrameQuadrupole_ZZ*gqzz16
              +2.0f*(atomJ.labFrameQuadrupole_XY*gqxy16+atomJ.labFrameQuadrupole_XZ*gqxz16+atomJ.labFrameQuadrupole_YZ*gqyz16))
              + atomI.labFrameQuadrupole_YZ*(atomJ.labFrameQuadrupole_XX*gqxx19+atomJ.labFrameQuadrupole_YY*gqyy19+atomJ.labFrameQuadrupole_ZZ*gqzz19
              +2.0f*(atomJ.labFrameQuadrupole_XY*gqxy19+atomJ.labFrameQuadrupole_XZ*gqxz19+atomJ.labFrameQuadrupole_YZ*gqyz19)));

    dedz = desymdz + 0.5f*(dewidz + dewkdz);

    desymdr = atomI.q * atomJ.q * gc21
                           - (atomI.labFrameDipole[0]*(atomJ.labFrameDipole[0]*gux22+atomJ.labFrameDipole[1]*guy22+atomJ.labFrameDipole[2]*guz22)
                             +atomI.labFrameDipole[1]*(atomJ.labFrameDipole[0]*gux23+atomJ.labFrameDipole[1]*guy23+atomJ.labFrameDipole[2]*guz23)
                             +atomI.labFrameDipole[2]*(atomJ.labFrameDipole[0]*gux24+atomJ.labFrameDipole[1]*guy24+atomJ.labFrameDipole[2]*guz24));

    dewidr = atomI.q*(atomJ.labFrameDipole[0]*gc22+atomJ.labFrameDipole[1]*gc23+atomJ.labFrameDipole[2]*gc24)
                      -atomJ.q*(atomI.labFrameDipole[0]*gux21+atomI.labFrameDipole[1]*guy21+atomI.labFrameDipole[2]*guz21)
                 +atomI.q*(atomJ.labFrameQuadrupole_XX*gc25+atomJ.labFrameQuadrupole_YY*gc28+atomJ.labFrameQuadrupole_ZZ*gc30
              +2.0f*(atomJ.labFrameQuadrupole_XY*gc26+atomJ.labFrameQuadrupole_XZ*gc27+atomJ.labFrameQuadrupole_YZ*gc29))
                 +atomJ.q*(atomI.labFrameQuadrupole_XX*gqxx21+atomI.labFrameQuadrupole_YY*gqyy21+atomI.labFrameQuadrupole_ZZ*gqzz21
              +2.0f*(atomI.labFrameQuadrupole_XY*gqxy21+atomI.labFrameQuadrupole_XZ*gqxz21+atomI.labFrameQuadrupole_YZ*gqyz21))
               - atomI.labFrameDipole[0]*(atomJ.labFrameQuadrupole_XX*gux25+atomJ.labFrameQuadrupole_YY*gux28+atomJ.labFrameQuadrupole_ZZ*gux30
              +2.0f*(atomJ.labFrameQuadrupole_XY*gux26+atomJ.labFrameQuadrupole_XZ*gux27+atomJ.labFrameQuadrupole_YZ*gux29))
               - atomI.labFrameDipole[1]*(atomJ.labFrameQuadrupole_XX*guy25+atomJ.labFrameQuadrupole_YY*guy28+atomJ.labFrameQuadrupole_ZZ*guy30
              +2.0f*(atomJ.labFrameQuadrupole_XY*guy26+atomJ.labFrameQuadrupole_XZ*guy27+atomJ.labFrameQuadrupole_YZ*guy29))
               - atomI.labFrameDipole[2]*(atomJ.labFrameQuadrupole_XX*guz25+atomJ.labFrameQuadrupole_YY*guz28+atomJ.labFrameQuadrupole_ZZ*guz30
              +2.0f*(atomJ.labFrameQuadrupole_XY*guz26+atomJ.labFrameQuadrupole_XZ*guz27+atomJ.labFrameQuadrupole_YZ*guz29))
               + atomJ.labFrameDipole[0]*(atomI.labFrameQuadrupole_XX*gqxx22+atomI.labFrameQuadrupole_YY*gqyy22+atomI.labFrameQuadrupole_ZZ*gqzz22
              +2.0f*(atomI.labFrameQuadrupole_XY*gqxy22+atomI.labFrameQuadrupole_XZ*gqxz22+atomI.labFrameQuadrupole_YZ*gqyz22))
               + atomJ.labFrameDipole[1]*(atomI.labFrameQuadrupole_XX*gqxx23+atomI.labFrameQuadrupole_YY*gqyy23+atomI.labFrameQuadrupole_ZZ*gqzz23
              +2.0f*(atomI.labFrameQuadrupole_XY*gqxy23+atomI.labFrameQuadrupole_XZ*gqxz23+atomI.labFrameQuadrupole_YZ*gqyz23))
               + atomJ.labFrameDipole[2]*(atomI.labFrameQuadrupole_XX*gqxx24+atomI.labFrameQuadrupole_YY*gqyy24+atomI.labFrameQuadrupole_ZZ*gqzz24
              +2.0f*(atomI.labFrameQuadrupole_XY*gqxy24+atomI.labFrameQuadrupole_XZ*gqxz24+atomI.labFrameQuadrupole_YZ*gqyz24))
              + atomI.labFrameQuadrupole_XX*(atomJ.labFrameQuadrupole_XX*gqxx25+atomJ.labFrameQuadrupole_YY*gqxx28+atomJ.labFrameQuadrupole_ZZ*gqxx30
              +2.0f*(atomJ.labFrameQuadrupole_XY*gqxx26+atomJ.labFrameQuadrupole_XZ*gqxx27+atomJ.labFrameQuadrupole_YZ*gqxx29))
              + atomI.labFrameQuadrupole_YY*(atomJ.labFrameQuadrupole_XX*gqyy25+atomJ.labFrameQuadrupole_YY*gqyy28+atomJ.labFrameQuadrupole_ZZ*gqyy30
              +2.0f*(atomJ.labFrameQuadrupole_XY*gqyy26+atomJ.labFrameQuadrupole_XZ*gqyy27+atomJ.labFrameQuadrupole_YZ*gqyy29))
              + atomI.labFrameQuadrupole_ZZ*(atomJ.labFrameQuadrupole_XX*gqzz25+atomJ.labFrameQuadrupole_YY*gqzz28+atomJ.labFrameQuadrupole_ZZ*gqzz30
              +2.0f*(atomJ.labFrameQuadrupole_XY*gqzz26+atomJ.labFrameQuadrupole_XZ*gqzz27+atomJ.labFrameQuadrupole_YZ*gqzz29))
              + 2.0f*(atomI.labFrameQuadrupole_XY*(atomJ.labFrameQuadrupole_XX*gqxy25+atomJ.labFrameQuadrupole_YY*gqxy28+atomJ.labFrameQuadrupole_ZZ*gqxy30
              +2.0f*(atomJ.labFrameQuadrupole_XY*gqxy26+atomJ.labFrameQuadrupole_XZ*gqxy27+atomJ.labFrameQuadrupole_YZ*gqxy29))
              + atomI.labFrameQuadrupole_XZ*(atomJ.labFrameQuadrupole_XX*gqxz25+atomJ.labFrameQuadrupole_YY*gqxz28+atomJ.labFrameQuadrupole_ZZ*gqxz30
              +2.0f*(atomJ.labFrameQuadrupole_XY*gqxz26+atomJ.labFrameQuadrupole_XZ*gqxz27+atomJ.labFrameQuadrupole_YZ*gqxz29))
              + atomI.labFrameQuadrupole_YZ*(atomJ.labFrameQuadrupole_XX*gqyz25+atomJ.labFrameQuadrupole_YY*gqyz28+atomJ.labFrameQuadrupole_ZZ*gqyz30
              +2.0f*(atomJ.labFrameQuadrupole_XY*gqyz26+atomJ.labFrameQuadrupole_XZ*gqyz27+atomJ.labFrameQuadrupole_YZ*gqyz29)));

    dewkdr = atomI.q*(atomJ.labFrameDipole[0]*gux21+atomJ.labFrameDipole[1]*guy21+atomJ.labFrameDipole[2]*guz21)
                      -atomJ.q*(atomI.labFrameDipole[0]*gc22+atomI.labFrameDipole[1]*gc23+atomI.labFrameDipole[2]*gc24)
                 +atomI.q*(atomJ.labFrameQuadrupole_XX*gqxx21+atomJ.labFrameQuadrupole_YY*gqyy21+atomJ.labFrameQuadrupole_ZZ*gqzz21
              +2.0f*(atomJ.labFrameQuadrupole_XY*gqxy21+atomJ.labFrameQuadrupole_XZ*gqxz21+atomJ.labFrameQuadrupole_YZ*gqyz21))
                 +atomJ.q*(atomI.labFrameQuadrupole_XX*gc25+atomI.labFrameQuadrupole_YY*gc28+atomI.labFrameQuadrupole_ZZ*gc30
              +2.0f*(atomI.labFrameQuadrupole_XY*gc26+atomI.labFrameQuadrupole_XZ*gc27+atomI.labFrameQuadrupole_YZ*gc29))
               - atomI.labFrameDipole[0]*(atomJ.labFrameQuadrupole_XX*gqxx22+atomJ.labFrameQuadrupole_YY*gqyy22+atomJ.labFrameQuadrupole_ZZ*gqzz22
              +2.0f*(atomJ.labFrameQuadrupole_XY*gqxy22+atomJ.labFrameQuadrupole_XZ*gqxz22+atomJ.labFrameQuadrupole_YZ*gqyz22))
               - atomI.labFrameDipole[1]*(atomJ.labFrameQuadrupole_XX*gqxx23+atomJ.labFrameQuadrupole_YY*gqyy23+atomJ.labFrameQuadrupole_ZZ*gqzz23
              +2.0f*(atomJ.labFrameQuadrupole_XY*gqxy23+atomJ.labFrameQuadrupole_XZ*gqxz23+atomJ.labFrameQuadrupole_YZ*gqyz23))
               - atomI.labFrameDipole[2]*(atomJ.labFrameQuadrupole_XX*gqxx24+atomJ.labFrameQuadrupole_YY*gqyy24+atomJ.labFrameQuadrupole_ZZ*gqzz24
              +2.0f*(atomJ.labFrameQuadrupole_XY*gqxy24+atomJ.labFrameQuadrupole_XZ*gqxz24+atomJ.labFrameQuadrupole_YZ*gqyz24))
               + atomJ.labFrameDipole[0]*(atomI.labFrameQuadrupole_XX*gux25+atomI.labFrameQuadrupole_YY*gux28+atomI.labFrameQuadrupole_ZZ*gux30
              +2.0f*(atomI.labFrameQuadrupole_XY*gux26+atomI.labFrameQuadrupole_XZ*gux27+atomI.labFrameQuadrupole_YZ*gux29))
               + atomJ.labFrameDipole[1]*(atomI.labFrameQuadrupole_XX*guy25+atomI.labFrameQuadrupole_YY*guy28+atomI.labFrameQuadrupole_ZZ*guy30
              +2.0f*(atomI.labFrameQuadrupole_XY*guy26+atomI.labFrameQuadrupole_XZ*guy27+atomI.labFrameQuadrupole_YZ*guy29))
               + atomJ.labFrameDipole[2]*(atomI.labFrameQuadrupole_XX*guz25+atomI.labFrameQuadrupole_YY*guz28+atomI.labFrameQuadrupole_ZZ*guz30
              +2.0f*(atomI.labFrameQuadrupole_XY*guz26+atomI.labFrameQuadrupole_XZ*guz27+atomI.labFrameQuadrupole_YZ*guz29))
              + atomI.labFrameQuadrupole_XX*(atomJ.labFrameQuadrupole_XX*gqxx25+atomJ.labFrameQuadrupole_YY*gqyy25+atomJ.labFrameQuadrupole_ZZ*gqzz25
              +2.0f*(atomJ.labFrameQuadrupole_XY*gqxy25+atomJ.labFrameQuadrupole_XZ*gqxz25+atomJ.labFrameQuadrupole_YZ*gqyz25))
              + atomI.labFrameQuadrupole_YY*(atomJ.labFrameQuadrupole_XX*gqxx28+atomJ.labFrameQuadrupole_YY*gqyy28+atomJ.labFrameQuadrupole_ZZ*gqzz28
              +2.0f*(atomJ.labFrameQuadrupole_XY*gqxy28+atomJ.labFrameQuadrupole_XZ*gqxz28+atomJ.labFrameQuadrupole_YZ*gqyz28))
              + atomI.labFrameQuadrupole_ZZ*(atomJ.labFrameQuadrupole_XX*gqxx30+atomJ.labFrameQuadrupole_YY*gqyy30+atomJ.labFrameQuadrupole_ZZ*gqzz30
              +2.0f*(atomJ.labFrameQuadrupole_XY*gqxy30+atomJ.labFrameQuadrupole_XZ*gqxz30+atomJ.labFrameQuadrupole_YZ*gqyz30))
              + 2.0f*(atomI.labFrameQuadrupole_XY*(atomJ.labFrameQuadrupole_XX*gqxx26+atomJ.labFrameQuadrupole_YY*gqyy26+atomJ.labFrameQuadrupole_ZZ*gqzz26
              +2.0f*(atomJ.labFrameQuadrupole_XY*gqxy26+atomJ.labFrameQuadrupole_XZ*gqxz26+atomJ.labFrameQuadrupole_YZ*gqyz26))
              + atomI.labFrameQuadrupole_XZ*(atomJ.labFrameQuadrupole_XX*gqxx27+atomJ.labFrameQuadrupole_YY*gqyy27+atomJ.labFrameQuadrupole_ZZ*gqzz27
              +2.0f*(atomJ.labFrameQuadrupole_XY*gqxy27+atomJ.labFrameQuadrupole_XZ*gqxz27+atomJ.labFrameQuadrupole_YZ*gqyz27))
              + atomI.labFrameQuadrupole_YZ*(atomJ.labFrameQuadrupole_XX*gqxx29+atomJ.labFrameQuadrupole_YY*gqyy29+atomJ.labFrameQuadrupole_ZZ*gqzz29
              +2.0f*(atomJ.labFrameQuadrupole_XY*gqxy29+atomJ.labFrameQuadrupole_XZ*gqxz29+atomJ.labFrameQuadrupole_YZ*gqyz29)));

    dsumdr = desymdr + 0.5f*(dewidr + dewkdr);
    drbi = atomJ.bornRadius*dsumdr;
    drbk = atomI.bornRadius*dsumdr;

    // torque on permanent dipoles due to permanent reaction field

    float trq1   = 0.0f;
    float trq2   = 0.0f;
    float trq3   = 0.0f;

    float trq_k1 = 0.0f;
    float trq_k2 = 0.0f;
    float trq_k3 = 0.0f;

    if ( xr != 0.0f || yr != 0.0f || zr != 0.0f )
    {

        float fid1 = atomJ.labFrameDipole[0]*gux2 + atomJ.labFrameDipole[1]*gux3 + atomJ.labFrameDipole[2]*gux4
                + 0.5f*(atomJ.q*gux1+atomJ.labFrameQuadrupole_XX*gux5+atomJ.labFrameQuadrupole_YY*gux8+atomJ.labFrameQuadrupole_ZZ*gux10
                      +2.0f*(atomJ.labFrameQuadrupole_XY*gux6+atomJ.labFrameQuadrupole_XZ*gux7+atomJ.labFrameQuadrupole_YZ*gux9)
                      +atomJ.q*gc2+atomJ.labFrameQuadrupole_XX*gqxx2+atomJ.labFrameQuadrupole_YY*gqyy2+atomJ.labFrameQuadrupole_ZZ*gqzz2
                      +2.0f*(atomJ.labFrameQuadrupole_XY*gqxy2+atomJ.labFrameQuadrupole_XZ*gqxz2+atomJ.labFrameQuadrupole_YZ*gqyz2));

        float fid2 = atomJ.labFrameDipole[0]*guy2 + atomJ.labFrameDipole[1]*guy3 + atomJ.labFrameDipole[2]*guy4
                + 0.5f*(atomJ.q*guy1+atomJ.labFrameQuadrupole_XX*guy5+atomJ.labFrameQuadrupole_YY*guy8+atomJ.labFrameQuadrupole_ZZ*guy10
                      +2.0f*(atomJ.labFrameQuadrupole_XY*guy6+atomJ.labFrameQuadrupole_XZ*guy7+atomJ.labFrameQuadrupole_YZ*guy9)
                      +atomJ.q*gc3+atomJ.labFrameQuadrupole_XX*gqxx3+atomJ.labFrameQuadrupole_YY*gqyy3+atomJ.labFrameQuadrupole_ZZ*gqzz3
                      +2.0f*(atomJ.labFrameQuadrupole_XY*gqxy3+atomJ.labFrameQuadrupole_XZ*gqxz3+atomJ.labFrameQuadrupole_YZ*gqyz3));

        float fid3 = atomJ.labFrameDipole[0]*guz2 + atomJ.labFrameDipole[1]*guz3 + atomJ.labFrameDipole[2]*guz4
                + 0.5f*(atomJ.q*guz1+atomJ.labFrameQuadrupole_XX*guz5+atomJ.labFrameQuadrupole_YY*guz8+atomJ.labFrameQuadrupole_ZZ*guz10
                      +2.0f*(atomJ.labFrameQuadrupole_XY*guz6+atomJ.labFrameQuadrupole_XZ*guz7+atomJ.labFrameQuadrupole_YZ*guz9)
                      +atomJ.q*gc4+atomJ.labFrameQuadrupole_XX*gqxx4+atomJ.labFrameQuadrupole_YY*gqyy4+atomJ.labFrameQuadrupole_ZZ*gqzz4
                      +2.0f*(atomJ.labFrameQuadrupole_XY*gqxy4+atomJ.labFrameQuadrupole_XZ*gqxz4+atomJ.labFrameQuadrupole_YZ*gqyz4));

        float fkd1 = atomI.labFrameDipole[0]*gux2 + atomI.labFrameDipole[1]*gux3 + atomI.labFrameDipole[2]*gux4
                - 0.5f*(atomI.q*gux1+atomI.labFrameQuadrupole_XX*gux5+atomI.labFrameQuadrupole_YY*gux8+atomI.labFrameQuadrupole_ZZ*gux10
                      +2.0f*(atomI.labFrameQuadrupole_XY*gux6+atomI.labFrameQuadrupole_XZ*gux7+atomI.labFrameQuadrupole_YZ*gux9)
                      +atomI.q*gc2+atomI.labFrameQuadrupole_XX*gqxx2+atomI.labFrameQuadrupole_YY*gqyy2+atomI.labFrameQuadrupole_ZZ*gqzz2
                      +2.0f*(atomI.labFrameQuadrupole_XY*gqxy2+atomI.labFrameQuadrupole_XZ*gqxz2+atomI.labFrameQuadrupole_YZ*gqyz2));

        float fkd2 = atomI.labFrameDipole[0]*guy2 + atomI.labFrameDipole[1]*guy3 + atomI.labFrameDipole[2]*guy4
                - 0.5f*(atomI.q*guy1+atomI.labFrameQuadrupole_XX*guy5+atomI.labFrameQuadrupole_YY*guy8+atomI.labFrameQuadrupole_ZZ*guy10
                      +2.0f*(atomI.labFrameQuadrupole_XY*guy6+atomI.labFrameQuadrupole_XZ*guy7+atomI.labFrameQuadrupole_YZ*guy9)
                      +atomI.q*gc3+atomI.labFrameQuadrupole_XX*gqxx3+atomI.labFrameQuadrupole_YY*gqyy3+atomI.labFrameQuadrupole_ZZ*gqzz3
                      +2.0f*(atomI.labFrameQuadrupole_XY*gqxy3+atomI.labFrameQuadrupole_XZ*gqxz3+atomI.labFrameQuadrupole_YZ*gqyz3));

        float fkd3 = atomI.labFrameDipole[0]*guz2 + atomI.labFrameDipole[1]*guz3 + atomI.labFrameDipole[2]*guz4
                - 0.5f*(atomI.q*guz1+atomI.labFrameQuadrupole_XX*guz5+atomI.labFrameQuadrupole_YY*guz8+atomI.labFrameQuadrupole_ZZ*guz10
                      +2.0f*(atomI.labFrameQuadrupole_XY*guz6+atomI.labFrameQuadrupole_XZ*guz7+atomI.labFrameQuadrupole_YZ*guz9)
                      +atomI.q*gc4+atomI.labFrameQuadrupole_XX*gqxx4+atomI.labFrameQuadrupole_YY*gqyy4+atomI.labFrameQuadrupole_ZZ*gqzz4
                      +2.0f*(atomI.labFrameQuadrupole_XY*gqxy4+atomI.labFrameQuadrupole_XZ*gqxz4+atomI.labFrameQuadrupole_YZ*gqyz4));

        trq1    = atomI.labFrameDipole[1]*fid3 - atomI.labFrameDipole[2]*fid2;
        trq2    = atomI.labFrameDipole[2]*fid1 - atomI.labFrameDipole[0]*fid3;
        trq3    = atomI.labFrameDipole[0]*fid2 - atomI.labFrameDipole[1]*fid1;

        trq_k1  = atomJ.labFrameDipole[1]*fkd3 - atomJ.labFrameDipole[2]*fkd2;
        trq_k2  = atomJ.labFrameDipole[2]*fkd1 - atomJ.labFrameDipole[0]*fkd3;
        trq_k3  = atomJ.labFrameDipole[0]*fkd2 - atomJ.labFrameDipole[1]*fkd1;

        // torque on quadrupoles due to permanent reaction field gradient

        float fidg11 =
                - 0.5f*(atomJ.q*gqxx1+atomJ.labFrameDipole[0]*gqxx2+atomJ.labFrameDipole[1]*gqxx3+atomJ.labFrameDipole[2]*gqxx4
                      +atomJ.labFrameQuadrupole_XX*gqxx5+atomJ.labFrameQuadrupole_YY*gqxx8+atomJ.labFrameQuadrupole_ZZ*gqxx10
                      +2.0f*(atomJ.labFrameQuadrupole_XY*gqxx6+atomJ.labFrameQuadrupole_XZ*gqxx7+atomJ.labFrameQuadrupole_YZ*gqxx9)
                      +atomJ.q*gc5+atomJ.labFrameDipole[0]*gux5+atomJ.labFrameDipole[1]*guy5+atomJ.labFrameDipole[2]*guz5
                      +atomJ.labFrameQuadrupole_XX*gqxx5+atomJ.labFrameQuadrupole_YY*gqyy5+atomJ.labFrameQuadrupole_ZZ*gqzz5
                      +2.0f*(atomJ.labFrameQuadrupole_XY*gqxy5+atomJ.labFrameQuadrupole_XZ*gqxz5+atomJ.labFrameQuadrupole_YZ*gqyz5));

        float fidg12 =
                - 0.5f*(atomJ.q*gqxy1+atomJ.labFrameDipole[0]*gqxy2+atomJ.labFrameDipole[1]*gqxy3+atomJ.labFrameDipole[2]*gqxy4
                      +atomJ.labFrameQuadrupole_XX*gqxy5+atomJ.labFrameQuadrupole_YY*gqxy8+atomJ.labFrameQuadrupole_ZZ*gqxy10
                      +2.0f*(atomJ.labFrameQuadrupole_XY*gqxy6+atomJ.labFrameQuadrupole_XZ*gqxy7+atomJ.labFrameQuadrupole_YZ*gqxy9)
                      +atomJ.q*gc6+atomJ.labFrameDipole[0]*gux6+atomJ.labFrameDipole[1]*guy6+atomJ.labFrameDipole[2]*guz6
                      +atomJ.labFrameQuadrupole_XX*gqxx6+atomJ.labFrameQuadrupole_YY*gqyy6+atomJ.labFrameQuadrupole_ZZ*gqzz6
                      +2.0f*(atomJ.labFrameQuadrupole_XY*gqxy6+atomJ.labFrameQuadrupole_XZ*gqxz6+atomJ.labFrameQuadrupole_YZ*gqyz6));

        float fidg13 =
                - 0.5f*(atomJ.q*gqxz1+atomJ.labFrameDipole[0]*gqxz2+atomJ.labFrameDipole[1]*gqxz3+atomJ.labFrameDipole[2]*gqxz4
                      +atomJ.labFrameQuadrupole_XX*gqxz5+atomJ.labFrameQuadrupole_YY*gqxz8+atomJ.labFrameQuadrupole_ZZ*gqxz10
                      +2.0f*(atomJ.labFrameQuadrupole_XY*gqxz6+atomJ.labFrameQuadrupole_XZ*gqxz7+atomJ.labFrameQuadrupole_YZ*gqxz9)
                      +atomJ.q*gc7+atomJ.labFrameDipole[0]*gux7+atomJ.labFrameDipole[1]*guy7+atomJ.labFrameDipole[2]*guz7
                      +atomJ.labFrameQuadrupole_XX*gqxx7+atomJ.labFrameQuadrupole_YY*gqyy7+atomJ.labFrameQuadrupole_ZZ*gqzz7
                      +2.0f*(atomJ.labFrameQuadrupole_XY*gqxy7+atomJ.labFrameQuadrupole_XZ*gqxz7+atomJ.labFrameQuadrupole_YZ*gqyz7));

        float fidg22 =
                - 0.5f*(atomJ.q*gqyy1+atomJ.labFrameDipole[0]*gqyy2+atomJ.labFrameDipole[1]*gqyy3+atomJ.labFrameDipole[2]*gqyy4
                      +atomJ.labFrameQuadrupole_XX*gqyy5+atomJ.labFrameQuadrupole_YY*gqyy8+atomJ.labFrameQuadrupole_ZZ*gqyy10
                      +2.0f*(atomJ.labFrameQuadrupole_XY*gqyy6+atomJ.labFrameQuadrupole_XZ*gqyy7+atomJ.labFrameQuadrupole_YZ*gqyy9)
                      +atomJ.q*gc8+atomJ.labFrameDipole[0]*gux8+atomJ.labFrameDipole[1]*guy8+atomJ.labFrameDipole[2]*guz8
                      +atomJ.labFrameQuadrupole_XX*gqxx8+atomJ.labFrameQuadrupole_YY*gqyy8+atomJ.labFrameQuadrupole_ZZ*gqzz8
                      +2.0f*(atomJ.labFrameQuadrupole_XY*gqxy8+atomJ.labFrameQuadrupole_XZ*gqxz8+atomJ.labFrameQuadrupole_YZ*gqyz8));

        float fidg23 =
                - 0.5f*(atomJ.q*gqyz1+atomJ.labFrameDipole[0]*gqyz2+atomJ.labFrameDipole[1]*gqyz3+atomJ.labFrameDipole[2]*gqyz4
                      +atomJ.labFrameQuadrupole_XX*gqyz5+atomJ.labFrameQuadrupole_YY*gqyz8+atomJ.labFrameQuadrupole_ZZ*gqyz10
                      +2.0f*(atomJ.labFrameQuadrupole_XY*gqyz6+atomJ.labFrameQuadrupole_XZ*gqyz7+atomJ.labFrameQuadrupole_YZ*gqyz9)
                      +atomJ.q*gc9+atomJ.labFrameDipole[0]*gux9+atomJ.labFrameDipole[1]*guy9+atomJ.labFrameDipole[2]*guz9
                      +atomJ.labFrameQuadrupole_XX*gqxx9+atomJ.labFrameQuadrupole_YY*gqyy9+atomJ.labFrameQuadrupole_ZZ*gqzz9
                      +2.0f*(atomJ.labFrameQuadrupole_XY*gqxy9+atomJ.labFrameQuadrupole_XZ*gqxz9+atomJ.labFrameQuadrupole_YZ*gqyz9));

        float fidg33 =
                - 0.5f*(atomJ.q*gqzz1+atomJ.labFrameDipole[0]*gqzz2+atomJ.labFrameDipole[1]*gqzz3+atomJ.labFrameDipole[2]*gqzz4
                      +atomJ.labFrameQuadrupole_XX*gqzz5+atomJ.labFrameQuadrupole_YY*gqzz8+atomJ.labFrameQuadrupole_ZZ*gqzz10
                      +2.0f*(atomJ.labFrameQuadrupole_XY*gqzz6+atomJ.labFrameQuadrupole_XZ*gqzz7+atomJ.labFrameQuadrupole_YZ*gqzz9)
                      +atomJ.q*gc10+atomJ.labFrameDipole[0]*gux10+atomJ.labFrameDipole[1]*guy10+atomJ.labFrameDipole[2]*guz10
                      +atomJ.labFrameQuadrupole_XX*gqxx10+atomJ.labFrameQuadrupole_YY*gqyy10+atomJ.labFrameQuadrupole_ZZ*gqzz10
                   +2.0f*(atomJ.labFrameQuadrupole_XY*gqxy10+atomJ.labFrameQuadrupole_XZ*gqxz10+atomJ.labFrameQuadrupole_YZ*gqyz10));

        float fidg21 = fidg12;
        float fidg31 = fidg13;
        float fidg32 = fidg23;

        float fkdg11 =
                - 0.5f*(atomI.q*gqxx1-atomI.labFrameDipole[0]*gqxx2-atomI.labFrameDipole[1]*gqxx3-atomI.labFrameDipole[2] *gqxx4
                      +atomI.labFrameQuadrupole_XX*gqxx5+atomI.labFrameQuadrupole_YY*gqxx8+atomI.labFrameQuadrupole_ZZ*gqxx10
                      +2.0f*(atomI.labFrameQuadrupole_XY*gqxx6+atomI.labFrameQuadrupole_XZ*gqxx7+atomI.labFrameQuadrupole_YZ*gqxx9)
                      +atomI.q*gc5-atomI.labFrameDipole[0]*gux5-atomI.labFrameDipole[1]*guy5-atomI.labFrameDipole[2]*guz5
                      +atomI.labFrameQuadrupole_XX*gqxx5+atomI.labFrameQuadrupole_YY*gqyy5+atomI.labFrameQuadrupole_ZZ*gqzz5
                      +2.0f*(atomI.labFrameQuadrupole_XY*gqxy5+atomI.labFrameQuadrupole_XZ*gqxz5+atomI.labFrameQuadrupole_YZ*gqyz5));

        float fkdg12 =
                - 0.5f*(atomI.q*gqxy1-atomI.labFrameDipole[0]*gqxy2-atomI.labFrameDipole[1]*gqxy3-atomI.labFrameDipole[2]*gqxy4
                      +atomI.labFrameQuadrupole_XX*gqxy5+atomI.labFrameQuadrupole_YY*gqxy8+atomI.labFrameQuadrupole_ZZ*gqxy10
                      +2.0f*(atomI.labFrameQuadrupole_XY*gqxy6+atomI.labFrameQuadrupole_XZ*gqxy7+atomI.labFrameQuadrupole_YZ*gqxy9)
                      +atomI.q*gc6-atomI.labFrameDipole[0]*gux6-atomI.labFrameDipole[1]*guy6-atomI.labFrameDipole[2]*guz6
                      +atomI.labFrameQuadrupole_XX*gqxx6+atomI.labFrameQuadrupole_YY*gqyy6+atomI.labFrameQuadrupole_ZZ*gqzz6
                      +2.0f*(atomI.labFrameQuadrupole_XY*gqxy6+atomI.labFrameQuadrupole_XZ*gqxz6+atomI.labFrameQuadrupole_YZ*gqyz6));

        float fkdg13 =
                - 0.5f*(atomI.q*gqxz1-atomI.labFrameDipole[0]*gqxz2-atomI.labFrameDipole[1]*gqxz3-atomI.labFrameDipole[2]*gqxz4
                      +atomI.labFrameQuadrupole_XX*gqxz5+atomI.labFrameQuadrupole_YY*gqxz8+atomI.labFrameQuadrupole_ZZ*gqxz10
                      +2.0f*(atomI.labFrameQuadrupole_XY*gqxz6+atomI.labFrameQuadrupole_XZ*gqxz7+atomI.labFrameQuadrupole_YZ*gqxz9)
                      +atomI.q*gc7-atomI.labFrameDipole[0]*gux7-atomI.labFrameDipole[1]*guy7-atomI.labFrameDipole[2]*guz7
                      +atomI.labFrameQuadrupole_XX*gqxx7+atomI.labFrameQuadrupole_YY*gqyy7+atomI.labFrameQuadrupole_ZZ*gqzz7
                      +2.0f*(atomI.labFrameQuadrupole_XY*gqxy7+atomI.labFrameQuadrupole_XZ*gqxz7+atomI.labFrameQuadrupole_YZ*gqyz7));

        float fkdg22 =
                - 0.5f*(atomI.q*gqyy1-atomI.labFrameDipole[0]*gqyy2-atomI.labFrameDipole[1]*gqyy3-atomI.labFrameDipole[2]*gqyy4
                      +atomI.labFrameQuadrupole_XX*gqyy5+atomI.labFrameQuadrupole_YY*gqyy8+atomI.labFrameQuadrupole_ZZ*gqyy10
                      +2.0f*(atomI.labFrameQuadrupole_XY*gqyy6+atomI.labFrameQuadrupole_XZ*gqyy7+atomI.labFrameQuadrupole_YZ*gqyy9)
                      +atomI.q*gc8-atomI.labFrameDipole[0]*gux8-atomI.labFrameDipole[1]*guy8-atomI.labFrameDipole[2]*guz8
                      +atomI.labFrameQuadrupole_XX*gqxx8+atomI.labFrameQuadrupole_YY*gqyy8+atomI.labFrameQuadrupole_ZZ*gqzz8
                      +2.0f*(atomI.labFrameQuadrupole_XY*gqxy8+atomI.labFrameQuadrupole_XZ*gqxz8+atomI.labFrameQuadrupole_YZ*gqyz8));

        float fkdg23 =
                - 0.5f*(atomI.q*gqyz1-atomI.labFrameDipole[0]*gqyz2-atomI.labFrameDipole[1]*gqyz3-atomI.labFrameDipole[2]*gqyz4
                      +atomI.labFrameQuadrupole_XX*gqyz5+atomI.labFrameQuadrupole_YY*gqyz8+atomI.labFrameQuadrupole_ZZ*gqyz10
                      +2.0f*(atomI.labFrameQuadrupole_XY*gqyz6+atomI.labFrameQuadrupole_XZ*gqyz7+atomI.labFrameQuadrupole_YZ*gqyz9)
                      +atomI.q*gc9-atomI.labFrameDipole[0]*gux9-atomI.labFrameDipole[1]*guy9-atomI.labFrameDipole[2]*guz9
                      +atomI.labFrameQuadrupole_XX*gqxx9+atomI.labFrameQuadrupole_YY*gqyy9+atomI.labFrameQuadrupole_ZZ*gqzz9
                      +2.0f*(atomI.labFrameQuadrupole_XY*gqxy9+atomI.labFrameQuadrupole_XZ*gqxz9+atomI.labFrameQuadrupole_YZ*gqyz9));
        float fkdg33 =
                - 0.5f*(atomI.q*gqzz1-atomI.labFrameDipole[0]*gqzz2-atomI.labFrameDipole[1]*gqzz3-atomI.labFrameDipole[2]*gqzz4
                      +atomI.labFrameQuadrupole_XX*gqzz5+atomI.labFrameQuadrupole_YY*gqzz8+atomI.labFrameQuadrupole_ZZ*gqzz10
                      +2.0f*(atomI.labFrameQuadrupole_XY*gqzz6+atomI.labFrameQuadrupole_XZ*gqzz7+atomI.labFrameQuadrupole_YZ*gqzz9)
                      +atomI.q*gc10-atomI.labFrameDipole[0]*gux10-atomI.labFrameDipole[1]*guy10-atomI.labFrameDipole[2]*guz10
                      +atomI.labFrameQuadrupole_XX*gqxx10+atomI.labFrameQuadrupole_YY*gqyy10+atomI.labFrameQuadrupole_ZZ*gqzz10
                    +2.0f*(atomI.labFrameQuadrupole_XY*gqxy10+atomI.labFrameQuadrupole_XZ*gqxz10+atomI.labFrameQuadrupole_YZ*gqyz10));

        float fkdg21 = fkdg12;
        float fkdg31 = fkdg13;
        float fkdg32 = fkdg23;

        trq1   += 2.0f* (atomI.labFrameQuadrupole_XY*fidg13+atomI.labFrameQuadrupole_YY*fidg23+atomI.labFrameQuadrupole_YZ*fidg33
                           -atomI.labFrameQuadrupole_XZ*fidg12-atomI.labFrameQuadrupole_YZ*fidg22-atomI.labFrameQuadrupole_ZZ*fidg32);

        trq2   += 2.0f*(atomI.labFrameQuadrupole_XZ*fidg11+atomI.labFrameQuadrupole_YZ*fidg21+atomI.labFrameQuadrupole_ZZ*fidg31
                         -atomI.labFrameQuadrupole_XX*fidg13-atomI.labFrameQuadrupole_XY*fidg23-atomI.labFrameQuadrupole_XZ*fidg33);

        trq3   += 2.0f*(atomI.labFrameQuadrupole_XX*fidg12+atomI.labFrameQuadrupole_XY*fidg22+atomI.labFrameQuadrupole_XZ*fidg32
                         -atomI.labFrameQuadrupole_XY*fidg11-atomI.labFrameQuadrupole_YY*fidg21-atomI.labFrameQuadrupole_YZ*fidg31);

        trq_k1 += 2.0f*
                          (atomJ.labFrameQuadrupole_XY*fkdg13+atomJ.labFrameQuadrupole_YY*fkdg23+atomJ.labFrameQuadrupole_YZ*fkdg33
                          -atomJ.labFrameQuadrupole_XZ*fkdg12-atomJ.labFrameQuadrupole_YZ*fkdg22-atomJ.labFrameQuadrupole_ZZ*fkdg32);

        trq_k2 += 2.0f*
                          (atomJ.labFrameQuadrupole_XZ*fkdg11+atomJ.labFrameQuadrupole_YZ*fkdg21+atomJ.labFrameQuadrupole_ZZ*fkdg31
                          -atomJ.labFrameQuadrupole_XX*fkdg13-atomJ.labFrameQuadrupole_XY*fkdg23-atomJ.labFrameQuadrupole_XZ*fkdg33);

        trq_k3 += 2.0f*
                          (atomJ.labFrameQuadrupole_XX*fkdg12+atomJ.labFrameQuadrupole_XY*fkdg22+atomJ.labFrameQuadrupole_XZ*fkdg32
                          -atomJ.labFrameQuadrupole_XY*fkdg11-atomJ.labFrameQuadrupole_YY*fkdg21-atomJ.labFrameQuadrupole_YZ*fkdg31);
    }

    // electrostatic solvation energy of the permanent multipoles in
    // the GK reaction potential of the induced dipoles

    esymi =              -atomI.labFrameDipole[0]*(atomJ.inducedDipole[0]*gux2+atomJ.inducedDipole[1]*guy2+atomJ.inducedDipole[2]*guz2)
                        - atomI.labFrameDipole[1]*(atomJ.inducedDipole[0]*gux3+atomJ.inducedDipole[1]*guy3+atomJ.inducedDipole[2]*guz3)
                        - atomI.labFrameDipole[2]*(atomJ.inducedDipole[0]*gux4+atomJ.inducedDipole[1]*guy4+atomJ.inducedDipole[2]*guz4)
                        - atomJ.labFrameDipole[0]*(atomI.inducedDipole[0]*gux2+atomI.inducedDipole[1]*guy2+atomI.inducedDipole[2]*guz2)
                        - atomJ.labFrameDipole[1]*(atomI.inducedDipole[0]*gux3+atomI.inducedDipole[1]*guy3+atomI.inducedDipole[2]*guz3)
                        - atomJ.labFrameDipole[2]*(atomI.inducedDipole[0]*gux4+atomI.inducedDipole[1]*guy4+atomI.inducedDipole[2]*guz4);

    ewii = atomI.q*(atomJ.inducedDipole[0]*gc2+atomJ.inducedDipole[1]*gc3+atomJ.inducedDipole[2]*gc4)
                      - atomJ.q*(atomI.inducedDipole[0]*gux1+atomI.inducedDipole[1]*guy1+atomI.inducedDipole[2]*guz1)
                      - atomI.inducedDipole[0]*(atomJ.labFrameQuadrupole_XX*gux5+atomJ.labFrameQuadrupole_YY*gux8+atomJ.labFrameQuadrupole_ZZ*gux10
                     +2.0f*(atomJ.labFrameQuadrupole_XY*gux6+atomJ.labFrameQuadrupole_XZ*gux7+atomJ.labFrameQuadrupole_YZ*gux9))
                      - atomI.inducedDipole[1]*(atomJ.labFrameQuadrupole_XX*guy5+atomJ.labFrameQuadrupole_YY*guy8+atomJ.labFrameQuadrupole_ZZ*guy10
                     +2.0f*(atomJ.labFrameQuadrupole_XY*guy6+atomJ.labFrameQuadrupole_XZ*guy7+atomJ.labFrameQuadrupole_YZ*guy9))
                      - atomI.inducedDipole[2]*(atomJ.labFrameQuadrupole_XX*guz5+atomJ.labFrameQuadrupole_YY*guz8+atomJ.labFrameQuadrupole_ZZ*guz10
                     +2.0f*(atomJ.labFrameQuadrupole_XY*guz6+atomJ.labFrameQuadrupole_XZ*guz7+atomJ.labFrameQuadrupole_YZ*guz9))
                      + atomJ.inducedDipole[0]*(atomI.labFrameQuadrupole_XX*gqxx2+atomI.labFrameQuadrupole_YY*gqyy2+atomI.labFrameQuadrupole_ZZ*gqzz2
                     +2.0f*(atomI.labFrameQuadrupole_XY*gqxy2+atomI.labFrameQuadrupole_XZ*gqxz2+atomI.labFrameQuadrupole_YZ*gqyz2))
                      + atomJ.inducedDipole[1]*(atomI.labFrameQuadrupole_XX*gqxx3+atomI.labFrameQuadrupole_YY*gqyy3+atomI.labFrameQuadrupole_ZZ*gqzz3
                     +2.0f*(atomI.labFrameQuadrupole_XY*gqxy3+atomI.labFrameQuadrupole_XZ*gqxz3+atomI.labFrameQuadrupole_YZ*gqyz3))
                      + atomJ.inducedDipole[2]*(atomI.labFrameQuadrupole_XX*gqxx4+atomI.labFrameQuadrupole_YY*gqyy4+atomI.labFrameQuadrupole_ZZ*gqzz4
                     +2.0f*(atomI.labFrameQuadrupole_XY*gqxy4+atomI.labFrameQuadrupole_XZ*gqxz4+atomI.labFrameQuadrupole_YZ*gqyz4));

    ewki = atomI.q*(atomJ.inducedDipole[0]*gux1+atomJ.inducedDipole[1]*guy1+atomJ.inducedDipole[2]*guz1)
                      - atomJ.q*(atomI.inducedDipole[0]*gc2+atomI.inducedDipole[1]*gc3+atomI.inducedDipole[2]*gc4)
                      - atomI.inducedDipole[0]*(atomJ.labFrameQuadrupole_XX*gqxx2+atomJ.labFrameQuadrupole_YY*gqyy2+atomJ.labFrameQuadrupole_ZZ*gqzz2
                     +2.0f*(atomJ.labFrameQuadrupole_XY*gqxy2+atomJ.labFrameQuadrupole_XZ*gqxz2+atomJ.labFrameQuadrupole_YZ*gqyz2))
                      - atomI.inducedDipole[1]*(atomJ.labFrameQuadrupole_XX*gqxx3+atomJ.labFrameQuadrupole_YY*gqyy3+atomJ.labFrameQuadrupole_ZZ*gqzz3
                     +2.0f*(atomJ.labFrameQuadrupole_XY*gqxy3+atomJ.labFrameQuadrupole_XZ*gqxz3+atomJ.labFrameQuadrupole_YZ*gqyz3))
                      - atomI.inducedDipole[2]*(atomJ.labFrameQuadrupole_XX*gqxx4+atomJ.labFrameQuadrupole_YY*gqyy4+atomJ.labFrameQuadrupole_ZZ*gqzz4
                     +2.0f*(atomJ.labFrameQuadrupole_XY*gqxy4+atomJ.labFrameQuadrupole_XZ*gqxz4+atomJ.labFrameQuadrupole_YZ*gqyz4))
                      + atomJ.inducedDipole[0]*(atomI.labFrameQuadrupole_XX*gux5+atomI.labFrameQuadrupole_YY*gux8+atomI.labFrameQuadrupole_ZZ*gux10
                     +2.0f*(atomI.labFrameQuadrupole_XY*gux6+atomI.labFrameQuadrupole_XZ*gux7+atomI.labFrameQuadrupole_YZ*gux9))
                      + atomJ.inducedDipole[1]*(atomI.labFrameQuadrupole_XX*guy5+atomI.labFrameQuadrupole_YY*guy8+atomI.labFrameQuadrupole_ZZ*guy10
                     +2.0f*(atomI.labFrameQuadrupole_XY*guy6+atomI.labFrameQuadrupole_XZ*guy7+atomI.labFrameQuadrupole_YZ*guy9))
                      + atomJ.inducedDipole[2]*(atomI.labFrameQuadrupole_XX*guz5+atomI.labFrameQuadrupole_YY*guz8+atomI.labFrameQuadrupole_ZZ*guz10
                     +2.0f*(atomI.labFrameQuadrupole_XY*guz6+atomI.labFrameQuadrupole_XZ*guz7+atomI.labFrameQuadrupole_YZ*guz9));

    // electrostatic solvation free energy gradient of the permanent
    // multipoles in the reaction potential of the induced dipoles

    float dpsymdx =       - atomI.labFrameDipole[0]*(sxk*gux5+syk*guy5+szk*guz5)
                          - atomI.labFrameDipole[1]*(sxk*gux6+syk*guy6+szk*guz6)
                          - atomI.labFrameDipole[2]*(sxk*gux7+syk*guy7+szk*guz7)

                          - atomJ.labFrameDipole[0]*(sxi*gux5+syi*guy5+szi*guz5)
                          - atomJ.labFrameDipole[1]*(sxi*gux6+syi*guy6+szi*guz6)
                          - atomJ.labFrameDipole[2]*(sxi*gux7+syi*guy7+szi*guz7);

    dpwidx = atomI.q*(sxk*gc5+syk*gc6+szk*gc7)
                        - atomJ.q*(sxi*gux2+syi*guy2+szi*guz2)
                      - sxi*(atomJ.labFrameQuadrupole_XX*gux11+atomJ.labFrameQuadrupole_YY*gux14+atomJ.labFrameQuadrupole_ZZ*gux16
                     +2.0f*(atomJ.labFrameQuadrupole_XY*gux12+atomJ.labFrameQuadrupole_XZ*gux13+atomJ.labFrameQuadrupole_YZ*gux15))
                      - syi*(atomJ.labFrameQuadrupole_XX*guy11+atomJ.labFrameQuadrupole_YY*guy14+atomJ.labFrameQuadrupole_ZZ*guy16
                     +2.0f*(atomJ.labFrameQuadrupole_XY*guy12+atomJ.labFrameQuadrupole_XZ*guy13+atomJ.labFrameQuadrupole_YZ*guy15))
                      - szi*(atomJ.labFrameQuadrupole_XX*guz11+atomJ.labFrameQuadrupole_YY*guz14+atomJ.labFrameQuadrupole_ZZ*guz16
                     +2.0f*(atomJ.labFrameQuadrupole_XY*guz12+atomJ.labFrameQuadrupole_XZ*guz13+atomJ.labFrameQuadrupole_YZ*guz15))
                      + sxk*(atomI.labFrameQuadrupole_XX*gqxx5+atomI.labFrameQuadrupole_YY*gqyy5+atomI.labFrameQuadrupole_ZZ*gqzz5
                     +2.0f*(atomI.labFrameQuadrupole_XY*gqxy5+atomI.labFrameQuadrupole_XZ*gqxz5+atomI.labFrameQuadrupole_YZ*gqyz5))
                      + syk*(atomI.labFrameQuadrupole_XX*gqxx6+atomI.labFrameQuadrupole_YY*gqyy6+atomI.labFrameQuadrupole_ZZ*gqzz6
                     +2.0f*(atomI.labFrameQuadrupole_XY*gqxy6+atomI.labFrameQuadrupole_XZ*gqxz6+atomI.labFrameQuadrupole_YZ*gqyz6))
                      + szk*(atomI.labFrameQuadrupole_XX*gqxx7+atomI.labFrameQuadrupole_YY*gqyy7+atomI.labFrameQuadrupole_ZZ*gqzz7
                     +2.0f*(atomI.labFrameQuadrupole_XY*gqxy7+atomI.labFrameQuadrupole_XZ*gqxz7+atomI.labFrameQuadrupole_YZ*gqyz7));

    dpwkdx = atomI.q*(sxk*gux2+syk*guy2+szk*guz2)
                        - atomJ.q*(sxi*gc5+syi*gc6+szi*gc7)
                      - sxi*(atomJ.labFrameQuadrupole_XX*gqxx5+atomJ.labFrameQuadrupole_YY*gqyy5+atomJ.labFrameQuadrupole_ZZ*gqzz5
                     +2.0f*(atomJ.labFrameQuadrupole_XY*gqxy5+atomJ.labFrameQuadrupole_XZ*gqxz5+atomJ.labFrameQuadrupole_YZ*gqyz5))
                      - syi*(atomJ.labFrameQuadrupole_XX*gqxx6+atomJ.labFrameQuadrupole_YY*gqyy6+atomJ.labFrameQuadrupole_ZZ*gqzz6
                     +2.0f*(atomJ.labFrameQuadrupole_XY*gqxy6+atomJ.labFrameQuadrupole_XZ*gqxz6+atomJ.labFrameQuadrupole_YZ*gqyz6))
                      - szi*(atomJ.labFrameQuadrupole_XX*gqxx7+atomJ.labFrameQuadrupole_YY*gqyy7+atomJ.labFrameQuadrupole_ZZ*gqzz7
                     +2.0f*(atomJ.labFrameQuadrupole_XY*gqxy7+atomJ.labFrameQuadrupole_XZ*gqxz7+atomJ.labFrameQuadrupole_YZ*gqyz7))
                      + sxk*(atomI.labFrameQuadrupole_XX*gux11+atomI.labFrameQuadrupole_YY*gux14+atomI.labFrameQuadrupole_ZZ*gux16
                     +2.0f*(atomI.labFrameQuadrupole_XY*gux12+atomI.labFrameQuadrupole_XZ*gux13+atomI.labFrameQuadrupole_YZ*gux15))
                      + syk*(atomI.labFrameQuadrupole_XX*guy11+atomI.labFrameQuadrupole_YY*guy14+atomI.labFrameQuadrupole_ZZ*guy16
                     +2.0f*(atomI.labFrameQuadrupole_XY*guy12+atomI.labFrameQuadrupole_XZ*guy13+atomI.labFrameQuadrupole_YZ*guy15))
                      + szk*(atomI.labFrameQuadrupole_XX*guz11+atomI.labFrameQuadrupole_YY*guz14+atomI.labFrameQuadrupole_ZZ*guz16
                     +2.0f*(atomI.labFrameQuadrupole_XY*guz12+atomI.labFrameQuadrupole_XZ*guz13+atomI.labFrameQuadrupole_YZ*guz15));

    dpdx = 0.5f * (dpsymdx + 0.5f*(dpwidx + dpwkdx));

    dpsymdy = -atomI.labFrameDipole[0]*(sxk*gux6+syk*guy6+szk*guz6)
                          - atomI.labFrameDipole[1]*(sxk*gux8+syk*guy8+szk*guz8)
                          - atomI.labFrameDipole[2]*(sxk*gux9+syk*guy9+szk*guz9)
                          - atomJ.labFrameDipole[0]*(sxi*gux6+syi*guy6+szi*guz6)
                          - atomJ.labFrameDipole[1]*(sxi*gux8+syi*guy8+szi*guz8)
                          - atomJ.labFrameDipole[2]*(sxi*gux9+syi*guy9+szi*guz9);

    dpwidy = atomI.q*(sxk*gc6+syk*gc8+szk*gc9)
                        - atomJ.q*(sxi*gux3+syi*guy3+szi*guz3)
                         - sxi*(atomJ.labFrameQuadrupole_XX*gux12+atomJ.labFrameQuadrupole_YY*gux17+atomJ.labFrameQuadrupole_ZZ*gux19
                        +2.0f*(atomJ.labFrameQuadrupole_XY*gux14+atomJ.labFrameQuadrupole_XZ*gux15+atomJ.labFrameQuadrupole_YZ*gux18))
                         - syi*(atomJ.labFrameQuadrupole_XX*guy12+atomJ.labFrameQuadrupole_YY*guy17+atomJ.labFrameQuadrupole_ZZ*guy19
                        +2.0f*(atomJ.labFrameQuadrupole_XY*guy14+atomJ.labFrameQuadrupole_XZ*guy15+atomJ.labFrameQuadrupole_YZ*guy18))
                         - szi*(atomJ.labFrameQuadrupole_XX*guz12+atomJ.labFrameQuadrupole_YY*guz17+atomJ.labFrameQuadrupole_ZZ*guz19
                        +2.0f*(atomJ.labFrameQuadrupole_XY*guz14+atomJ.labFrameQuadrupole_XZ*guz15+atomJ.labFrameQuadrupole_YZ*guz18))
                         + sxk*(atomI.labFrameQuadrupole_XX*gqxx6+atomI.labFrameQuadrupole_YY*gqyy6+atomI.labFrameQuadrupole_ZZ*gqzz6
                        +2.0f*(atomI.labFrameQuadrupole_XY*gqxy6+atomI.labFrameQuadrupole_XZ*gqxz6+atomI.labFrameQuadrupole_YZ*gqyz6))
                         + syk*(atomI.labFrameQuadrupole_XX*gqxx8+atomI.labFrameQuadrupole_YY*gqyy8+atomI.labFrameQuadrupole_ZZ*gqzz8
                        +2.0f*(atomI.labFrameQuadrupole_XY*gqxy8+atomI.labFrameQuadrupole_XZ*gqxz8+atomI.labFrameQuadrupole_YZ*gqyz8))
                         + szk*(atomI.labFrameQuadrupole_XX*gqxx9+atomI.labFrameQuadrupole_YY*gqyy9+atomI.labFrameQuadrupole_ZZ*gqzz9
                        +2.0f*(atomI.labFrameQuadrupole_XY*gqxy9+atomI.labFrameQuadrupole_XZ*gqxz9+atomI.labFrameQuadrupole_YZ*gqyz9));

    dpwkdy = atomI.q*(sxk*gux3+syk*guy3+szk*guz3)
                        - atomJ.q*(sxi*gc6+syi*gc8+szi*gc9)
                      - sxi*(atomJ.labFrameQuadrupole_XX*gqxx6+atomJ.labFrameQuadrupole_YY*gqyy6+atomJ.labFrameQuadrupole_ZZ*gqzz6
                     +2.0f*(atomJ.labFrameQuadrupole_XY*gqxy6+atomJ.labFrameQuadrupole_XZ*gqxz6+atomJ.labFrameQuadrupole_YZ*gqyz6))
                      - syi*(atomJ.labFrameQuadrupole_XX*gqxx8+atomJ.labFrameQuadrupole_YY*gqyy8+atomJ.labFrameQuadrupole_ZZ*gqzz8
                     +2.0f*(atomJ.labFrameQuadrupole_XY*gqxy8+atomJ.labFrameQuadrupole_XZ*gqxz8+atomJ.labFrameQuadrupole_YZ*gqyz8))
                      - szi*(atomJ.labFrameQuadrupole_XX*gqxx9+atomJ.labFrameQuadrupole_YY*gqyy9+atomJ.labFrameQuadrupole_ZZ*gqzz9
                     +2.0f*(atomJ.labFrameQuadrupole_XY*gqxy9+atomJ.labFrameQuadrupole_XZ*gqxz9+atomJ.labFrameQuadrupole_YZ*gqyz9))
                      + sxk*(atomI.labFrameQuadrupole_XX*gux12+atomI.labFrameQuadrupole_YY*gux17+atomI.labFrameQuadrupole_ZZ*gux19
                     +2.0f*(atomI.labFrameQuadrupole_XY*gux14+atomI.labFrameQuadrupole_XZ*gux15+atomI.labFrameQuadrupole_YZ*gux18))
                      + syk*(atomI.labFrameQuadrupole_XX*guy12+atomI.labFrameQuadrupole_YY*guy17+atomI.labFrameQuadrupole_ZZ*guy19
                     +2.0f*(atomI.labFrameQuadrupole_XY*guy14+atomI.labFrameQuadrupole_XZ*guy15+atomI.labFrameQuadrupole_YZ*guy18))
                      + szk*(atomI.labFrameQuadrupole_XX*guz12+atomI.labFrameQuadrupole_YY*guz17+atomI.labFrameQuadrupole_ZZ*guz19
                     +2.0f*(atomI.labFrameQuadrupole_XY*guz14+atomI.labFrameQuadrupole_XZ*guz15+atomI.labFrameQuadrupole_YZ*guz18));

    dpdy    = 0.5f * (dpsymdy + 0.5f*(dpwidy + dpwkdy));

    dpsymdz = -atomI.labFrameDipole[0]*(sxk*gux7+syk*guy7+szk*guz7)
                          - atomI.labFrameDipole[1]*(sxk*gux9+syk*guy9+szk*guz9)
                          - atomI.labFrameDipole[2]*(sxk*gux10+syk*guy10+szk*guz10)
                          - atomJ.labFrameDipole[0]*(sxi*gux7+syi*guy7+szi*guz7)
                          - atomJ.labFrameDipole[1]*(sxi*gux9+syi*guy9+szi*guz9)
                          - atomJ.labFrameDipole[2]*(sxi*gux10+syi*guy10+szi*guz10);

    dpwidz = atomI.q*(sxk*gc7+syk*gc9+szk*gc10)
                        - atomJ.q*(sxi*gux4+syi*guy4+szi*guz4)
                      - sxi*(atomJ.labFrameQuadrupole_XX*gux13+atomJ.labFrameQuadrupole_YY*gux18+atomJ.labFrameQuadrupole_ZZ*gux20
                     +2.0f*(atomJ.labFrameQuadrupole_XY*gux15+atomJ.labFrameQuadrupole_XZ*gux16+atomJ.labFrameQuadrupole_YZ*gux19))
                      - syi*(atomJ.labFrameQuadrupole_XX*guy13+atomJ.labFrameQuadrupole_YY*guy18+atomJ.labFrameQuadrupole_ZZ*guy20
                     +2.0f*(atomJ.labFrameQuadrupole_XY*guy15+atomJ.labFrameQuadrupole_XZ*guy16+atomJ.labFrameQuadrupole_YZ*guy19))
                      - szi*(atomJ.labFrameQuadrupole_XX*guz13+atomJ.labFrameQuadrupole_YY*guz18+atomJ.labFrameQuadrupole_ZZ*guz20
                     +2.0f*(atomJ.labFrameQuadrupole_XY*guz15+atomJ.labFrameQuadrupole_XZ*guz16+atomJ.labFrameQuadrupole_YZ*guz19))
                      + sxk*(atomI.labFrameQuadrupole_XX*gqxx7+atomI.labFrameQuadrupole_YY*gqyy7+atomI.labFrameQuadrupole_ZZ*gqzz7
                     +2.0f*(atomI.labFrameQuadrupole_XY*gqxy7+atomI.labFrameQuadrupole_XZ*gqxz7+atomI.labFrameQuadrupole_YZ*gqyz7))
                      + syk*(atomI.labFrameQuadrupole_XX*gqxx9+atomI.labFrameQuadrupole_YY*gqyy9+atomI.labFrameQuadrupole_ZZ*gqzz9
                     +2.0f*(atomI.labFrameQuadrupole_XY*gqxy9+atomI.labFrameQuadrupole_XZ*gqxz9+atomI.labFrameQuadrupole_YZ*gqyz9))
                      + szk*(atomI.labFrameQuadrupole_XX*gqxx10+atomI.labFrameQuadrupole_YY*gqyy10+atomI.labFrameQuadrupole_ZZ*gqzz10
                     +2.0f*(atomI.labFrameQuadrupole_XY*gqxy10+atomI.labFrameQuadrupole_XZ*gqxz10+atomI.labFrameQuadrupole_YZ*gqyz10));

    dpwkdz = atomI.q*(sxk*gux4+syk*guy4+szk*guz4)
                        - atomJ.q*(sxi*gc7+syi*gc9+szi*gc10)
                      - sxi*(atomJ.labFrameQuadrupole_XX*gqxx7+atomJ.labFrameQuadrupole_YY*gqyy7+atomJ.labFrameQuadrupole_ZZ*gqzz7
                     +2.0f*(atomJ.labFrameQuadrupole_XY*gqxy7+atomJ.labFrameQuadrupole_XZ*gqxz7+atomJ.labFrameQuadrupole_YZ*gqyz7))
                      - syi*(atomJ.labFrameQuadrupole_XX*gqxx9+atomJ.labFrameQuadrupole_YY*gqyy9+atomJ.labFrameQuadrupole_ZZ*gqzz9
                     +2.0f*(atomJ.labFrameQuadrupole_XY*gqxy9+atomJ.labFrameQuadrupole_XZ*gqxz9+atomJ.labFrameQuadrupole_YZ*gqyz9))
                      - szi*(atomJ.labFrameQuadrupole_XX*gqxx10+atomJ.labFrameQuadrupole_YY*gqyy10+atomJ.labFrameQuadrupole_ZZ*gqzz10
                     +2.0f*(atomJ.labFrameQuadrupole_XY*gqxy10+atomJ.labFrameQuadrupole_XZ*gqxz10+atomJ.labFrameQuadrupole_YZ*gqyz10))
                      + sxk*(atomI.labFrameQuadrupole_XX*gux13+atomI.labFrameQuadrupole_YY*gux18+atomI.labFrameQuadrupole_ZZ*gux20
                     +2.0f*(atomI.labFrameQuadrupole_XY*gux15+atomI.labFrameQuadrupole_XZ*gux16+atomI.labFrameQuadrupole_YZ*gux19))
                      + syk*(atomI.labFrameQuadrupole_XX*guy13+atomI.labFrameQuadrupole_YY*guy18+atomI.labFrameQuadrupole_ZZ*guy20
                     +2.0f*(atomI.labFrameQuadrupole_XY*guy15+atomI.labFrameQuadrupole_XZ*guy16+atomI.labFrameQuadrupole_YZ*guy19))
                      + szk*(atomI.labFrameQuadrupole_XX*guz13+atomI.labFrameQuadrupole_YY*guz18+atomI.labFrameQuadrupole_ZZ*guz20
                     +2.0f*(atomI.labFrameQuadrupole_XY*guz15+atomI.labFrameQuadrupole_XZ*guz16+atomI.labFrameQuadrupole_YZ*guz19));

    dpdz = 0.5f * (dpsymdz + 0.5f*(dpwidz + dpwkdz));

    // effective radii chain rule terms for the
    // electrostatic solvation free energy gradient of the permanent
    // multipoles in the reaction potential of the induced dipoles

    dsymdr = -atomI.labFrameDipole[0]*(sxk*gux22+syk*guy22+szk*guz22)
                          - atomI.labFrameDipole[1]*(sxk*gux23+syk*guy23+szk*guz23)
                          - atomI.labFrameDipole[2]*(sxk*gux24+syk*guy24+szk*guz24)
                          - atomJ.labFrameDipole[0]*(sxi*gux22+syi*guy22+szi*guz22)
                          - atomJ.labFrameDipole[1]*(sxi*gux23+syi*guy23+szi*guz23)
                          - atomJ.labFrameDipole[2]*(sxi*gux24+syi*guy24+szi*guz24);

    dwipdr = atomI.q*(sxk*gc22+syk*gc23+szk*gc24)
                         - atomJ.q*(sxi*gux21+syi*guy21+szi*guz21)
                      - sxi*(atomJ.labFrameQuadrupole_XX*gux25+atomJ.labFrameQuadrupole_YY*gux28+atomJ.labFrameQuadrupole_ZZ*gux30
                     +2.0f*(atomJ.labFrameQuadrupole_XY*gux26+atomJ.labFrameQuadrupole_XZ*gux27+atomJ.labFrameQuadrupole_YZ*gux29))
                      - syi*(atomJ.labFrameQuadrupole_XX*guy25+atomJ.labFrameQuadrupole_YY*guy28+atomJ.labFrameQuadrupole_ZZ*guy30
                     +2.0f*(atomJ.labFrameQuadrupole_XY*guy26+atomJ.labFrameQuadrupole_XZ*guy27+atomJ.labFrameQuadrupole_YZ*guy29))
                      - szi*(atomJ.labFrameQuadrupole_XX*guz25+atomJ.labFrameQuadrupole_YY*guz28+atomJ.labFrameQuadrupole_ZZ*guz30
                     +2.0f*(atomJ.labFrameQuadrupole_XY*guz26+atomJ.labFrameQuadrupole_XZ*guz27+atomJ.labFrameQuadrupole_YZ*guz29))
                      + sxk*(atomI.labFrameQuadrupole_XX*gqxx22+atomI.labFrameQuadrupole_YY*gqyy22+atomI.labFrameQuadrupole_ZZ*gqzz22
                     +2.0f*(atomI.labFrameQuadrupole_XY*gqxy22+atomI.labFrameQuadrupole_XZ*gqxz22+atomI.labFrameQuadrupole_YZ*gqyz22))
                      + syk*(atomI.labFrameQuadrupole_XX*gqxx23+atomI.labFrameQuadrupole_YY*gqyy23+atomI.labFrameQuadrupole_ZZ*gqzz23
                     +2.0f*(atomI.labFrameQuadrupole_XY*gqxy23+atomI.labFrameQuadrupole_XZ*gqxz23+atomI.labFrameQuadrupole_YZ*gqyz23))
                      + szk*(atomI.labFrameQuadrupole_XX*gqxx24+atomI.labFrameQuadrupole_YY*gqyy24+atomI.labFrameQuadrupole_ZZ*gqzz24
                     +2.0f*(atomI.labFrameQuadrupole_XY*gqxy24+atomI.labFrameQuadrupole_XZ*gqxz24+atomI.labFrameQuadrupole_YZ*gqyz24));

    dwkpdr = atomI.q*(sxk*gux21+syk*guy21+szk*guz21)
                         - atomJ.q*(sxi*gc22+syi*gc23+szi*gc24)
                      - sxi*(atomJ.labFrameQuadrupole_XX*gqxx22+atomJ.labFrameQuadrupole_YY*gqyy22+atomJ.labFrameQuadrupole_ZZ*gqzz22
                     +2.0f*(atomJ.labFrameQuadrupole_XY*gqxy22+atomJ.labFrameQuadrupole_XZ*gqxz22+atomJ.labFrameQuadrupole_YZ*gqyz22))
                      - syi*(atomJ.labFrameQuadrupole_XX*gqxx23+atomJ.labFrameQuadrupole_YY*gqyy23+atomJ.labFrameQuadrupole_ZZ*gqzz23
                     +2.0f*(atomJ.labFrameQuadrupole_XY*gqxy23+atomJ.labFrameQuadrupole_XZ*gqxz23+atomJ.labFrameQuadrupole_YZ*gqyz23))
                      - szi*(atomJ.labFrameQuadrupole_XX*gqxx24+atomJ.labFrameQuadrupole_YY*gqyy24+atomJ.labFrameQuadrupole_ZZ*gqzz24
                     +2.0f*(atomJ.labFrameQuadrupole_XY*gqxy24+atomJ.labFrameQuadrupole_XZ*gqxz24+atomJ.labFrameQuadrupole_YZ*gqyz24))
                      + sxk*(atomI.labFrameQuadrupole_XX*gux25+atomI.labFrameQuadrupole_YY*gux28+atomI.labFrameQuadrupole_ZZ*gux30
                     +2.0f*(atomI.labFrameQuadrupole_XY*gux26+atomI.labFrameQuadrupole_XZ*gux27+atomI.labFrameQuadrupole_YZ*gux29))
                      + syk*(atomI.labFrameQuadrupole_XX*guy25+atomI.labFrameQuadrupole_YY*guy28+atomI.labFrameQuadrupole_ZZ*guy30
                     +2.0f*(atomI.labFrameQuadrupole_XY*guy26+atomI.labFrameQuadrupole_XZ*guy27+atomI.labFrameQuadrupole_YZ*guy29))
                      + szk*(atomI.labFrameQuadrupole_XX*guz25+atomI.labFrameQuadrupole_YY*guz28+atomI.labFrameQuadrupole_ZZ*guz30
                     +2.0f*(atomI.labFrameQuadrupole_XY*guz26+atomI.labFrameQuadrupole_XZ*guz27+atomI.labFrameQuadrupole_YZ*guz29));

    dsumdr = dsymdr + 0.5f*(dwipdr + dwkpdr);
    dpbi = 0.5f*atomJ.bornRadius*dsumdr;
    dpbk = 0.5f*atomI.bornRadius*dsumdr;

    // mutual polarization electrostatic solvation free energy gradient

    if( cAmoebaSim.polarizationType == 0 ){

        dpdx = dpdx - 0.5f *
                           (atomI.inducedDipole[0]*(atomJ.inducedDipoleP[0]*gux5+atomJ.inducedDipoleP[1]*gux6+atomJ.inducedDipoleP[2]*gux7)
                           +atomI.inducedDipole[1]*(atomJ.inducedDipoleP[0]*guy5+atomJ.inducedDipoleP[1]*guy6+atomJ.inducedDipoleP[2]*guy7)
                           +atomI.inducedDipole[2]*(atomJ.inducedDipoleP[0]*guz5+atomJ.inducedDipoleP[1]*guz6+atomJ.inducedDipoleP[2]*guz7)
                           +atomJ.inducedDipole[0]*(atomI.inducedDipoleP[0]*gux5+atomI.inducedDipoleP[1]*gux6+atomI.inducedDipoleP[2]*gux7)
                           +atomJ.inducedDipole[1]*(atomI.inducedDipoleP[0]*guy5+atomI.inducedDipoleP[1]*guy6+atomI.inducedDipoleP[2]*guy7)
                           +atomJ.inducedDipole[2]*(atomI.inducedDipoleP[0]*guz5+atomI.inducedDipoleP[1]*guz6+atomI.inducedDipoleP[2]*guz7));

        dpdy = dpdy - 0.5f *
                           (atomI.inducedDipole[0]*(atomJ.inducedDipoleP[0]*gux6+atomJ.inducedDipoleP[1]*gux8+atomJ.inducedDipoleP[2]*gux9)
                           +atomI.inducedDipole[1]*(atomJ.inducedDipoleP[0]*guy6+atomJ.inducedDipoleP[1]*guy8+atomJ.inducedDipoleP[2]*guy9)
                           +atomI.inducedDipole[2]*(atomJ.inducedDipoleP[0]*guz6+atomJ.inducedDipoleP[1]*guz8+atomJ.inducedDipoleP[2]*guz9)
                           +atomJ.inducedDipole[0]*(atomI.inducedDipoleP[0]*gux6+atomI.inducedDipoleP[1]*gux8+atomI.inducedDipoleP[2]*gux9)
                           +atomJ.inducedDipole[1]*(atomI.inducedDipoleP[0]*guy6+atomI.inducedDipoleP[1]*guy8+atomI.inducedDipoleP[2]*guy9)
                           +atomJ.inducedDipole[2]*(atomI.inducedDipoleP[0]*guz6+atomI.inducedDipoleP[1]*guz8+atomI.inducedDipoleP[2]*guz9));

        dpdz = dpdz - 0.5f *
                           (atomI.inducedDipole[0]*(atomJ.inducedDipoleP[0]*gux7+atomJ.inducedDipoleP[1]*gux9+atomJ.inducedDipoleP[2]*gux10)
                           +atomI.inducedDipole[1]*(atomJ.inducedDipoleP[0]*guy7+atomJ.inducedDipoleP[1]*guy9+atomJ.inducedDipoleP[2]*guy10)
                           +atomI.inducedDipole[2]*(atomJ.inducedDipoleP[0]*guz7+atomJ.inducedDipoleP[1]*guz9+atomJ.inducedDipoleP[2]*guz10)
                           +atomJ.inducedDipole[0]*(atomI.inducedDipoleP[0]*gux7+atomI.inducedDipoleP[1]*gux9+atomI.inducedDipoleP[2]*gux10)
                           +atomJ.inducedDipole[1]*(atomI.inducedDipoleP[0]*guy7+atomI.inducedDipoleP[1]*guy9+atomI.inducedDipoleP[2]*guy10)
                           +atomJ.inducedDipole[2]*(atomI.inducedDipoleP[0]*guz7+atomI.inducedDipoleP[1]*guz9+atomI.inducedDipoleP[2]*guz10));

        duvdr = atomI.inducedDipole[0]*(atomJ.inducedDipoleP[0]*gux22+atomJ.inducedDipoleP[1]*gux23+atomJ.inducedDipoleP[2]*gux24)
                            + atomI.inducedDipole[1]*(atomJ.inducedDipoleP[0]*guy22+atomJ.inducedDipoleP[1]*guy23+atomJ.inducedDipoleP[2]*guy24)
                            + atomI.inducedDipole[2]*(atomJ.inducedDipoleP[0]*guz22+atomJ.inducedDipoleP[1]*guz23+atomJ.inducedDipoleP[2]*guz24)
                            + atomJ.inducedDipole[0]*(atomI.inducedDipoleP[0]*gux22+atomI.inducedDipoleP[1]*gux23+atomI.inducedDipoleP[2]*gux24)
                            + atomJ.inducedDipole[1]*(atomI.inducedDipoleP[0]*guy22+atomI.inducedDipoleP[1]*guy23+atomI.inducedDipoleP[2]*guy24)
                            + atomJ.inducedDipole[2]*(atomI.inducedDipoleP[0]*guz22+atomI.inducedDipoleP[1]*guz23+atomI.inducedDipoleP[2]*guz24);

        dpbi = dpbi - 0.5f*atomJ.bornRadius*duvdr;
        dpbk = dpbk - 0.5f*atomI.bornRadius*duvdr;
    }

    // torque due to induced reaction field on permanent dipoles

    float fid1 = 0.5f * (sxk*gux2+syk*guy2+szk*guz2);
    float fid2 = 0.5f * (sxk*gux3+syk*guy3+szk*guz3);
    float fid3 = 0.5f * (sxk*gux4+syk*guy4+szk*guz4);
    float fkd1 = 0.5f * (sxi*gux2+syi*guy2+szi*guz2);
    float fkd2 = 0.5f * (sxi*gux3+syi*guy3+szi*guz3);
    float fkd3 = 0.5f * (sxi*gux4+syi*guy4+szi*guz4);


    float trqi1   = atomI.labFrameDipole[1]*fid3 - atomI.labFrameDipole[2]*fid2;
    float trqi2   = atomI.labFrameDipole[2]*fid1 - atomI.labFrameDipole[0]*fid3;
    float trqi3   = atomI.labFrameDipole[0]*fid2 - atomI.labFrameDipole[1]*fid1;

    float trqi_k1 = atomJ.labFrameDipole[1]*fkd3 - atomJ.labFrameDipole[2]*fkd2;
    float trqi_k2 = atomJ.labFrameDipole[2]*fkd1 - atomJ.labFrameDipole[0]*fkd3;
    float trqi_k3 = atomJ.labFrameDipole[0]*fkd2 - atomJ.labFrameDipole[1]*fkd1;


    // torque due to induced reaction field gradient on quadrupoles;

    float fidg11 = -0.25f *
                              ( (sxk*gqxx2+syk*gqxx3+szk*gqxx4)
                              + (sxk*gux5+syk*guy5+szk*guz5));

    float fidg12 = -0.25f *
                              ( (sxk*gqxy2+syk*gqxy3+szk*gqxy4)
                              + (sxk*gux6+syk*guy6+szk*guz6));

    float fidg13 = -0.25f *
                              ( (sxk*gqxz2+syk*gqxz3+szk*gqxz4)
                              + (sxk*gux7+syk*guy7+szk*guz7));

    float fidg22 = -0.25f *
                              ( (sxk*gqyy2+syk*gqyy3+szk*gqyy4)
                              + (sxk*gux8+syk*guy8+szk*guz8));

    float fidg23 = -0.25f *
                              ( (sxk*gqyz2+syk*gqyz3+szk*gqyz4)
                              + (sxk*gux9+syk*guy9+szk*guz9));

    float fidg33 = -0.25f *
                              ( (sxk*gqzz2+syk*gqzz3+szk*gqzz4)
                              + (sxk*gux10+syk*guy10+szk*guz10));

    float fidg21 = fidg12;
    float fidg31 = fidg13;
    float fidg32 = fidg23;

    float fkdg11 = 0.25f *
                              ( (sxi*gqxx2+syi*gqxx3+szi*gqxx4)
                              + (sxi*gux5+syi*guy5+szi*guz5));

    float fkdg12 = 0.25f *
                              ( (sxi*gqxy2+syi*gqxy3+szi*gqxy4)
                              + (sxi*gux6+syi*guy6+szi*guz6));
    float fkdg13 = 0.25f *
                              ( (sxi*gqxz2+syi*gqxz3+szi*gqxz4)
                              + (sxi*gux7+syi*guy7+szi*guz7));
    float fkdg22 = 0.25f *
                              ( (sxi*gqyy2+syi*gqyy3+szi*gqyy4)
                              + (sxi*gux8+syi*guy8+szi*guz8));
    float fkdg23 = 0.25f *
                              ( (sxi*gqyz2+syi*gqyz3+szi*gqyz4)
                              + (sxi*gux9+syi*guy9+szi*guz9));
    float fkdg33 = 0.25f *
                              ( (sxi*gqzz2+syi*gqzz3+szi*gqzz4)
                              + (sxi*gux10+syi*guy10+szi*guz10));
    float fkdg21 = fkdg12;
    float fkdg31 = fkdg13;
    float fkdg32 = fkdg23;

    trqi1 += 2.0f*(atomI.labFrameQuadrupole_XY*fidg13+atomI.labFrameQuadrupole_YY*fidg23+atomI.labFrameQuadrupole_YZ*fidg33
                        -atomI.labFrameQuadrupole_XZ*fidg12-atomI.labFrameQuadrupole_YZ*fidg22-atomI.labFrameQuadrupole_ZZ*fidg32);
    trqi2 += 2.0f*(atomI.labFrameQuadrupole_XZ*fidg11+atomI.labFrameQuadrupole_YZ*fidg21+atomI.labFrameQuadrupole_ZZ*fidg31
                        -atomI.labFrameQuadrupole_XX*fidg13-atomI.labFrameQuadrupole_XY*fidg23-atomI.labFrameQuadrupole_XZ*fidg33);

    trqi3 += 2.0f*(atomI.labFrameQuadrupole_XX*fidg12+atomI.labFrameQuadrupole_XY*fidg22+atomI.labFrameQuadrupole_XZ*fidg32
                        -atomI.labFrameQuadrupole_XY*fidg11-atomI.labFrameQuadrupole_YY*fidg21-atomI.labFrameQuadrupole_YZ*fidg31);

    trqi_k1 += 2.0f*
                        (atomJ.labFrameQuadrupole_XY*fkdg13+atomJ.labFrameQuadrupole_YY*fkdg23+atomJ.labFrameQuadrupole_YZ*fkdg33
                        -atomJ.labFrameQuadrupole_XZ*fkdg12-atomJ.labFrameQuadrupole_YZ*fkdg22-atomJ.labFrameQuadrupole_ZZ*fkdg32);

    trqi_k2 += 2.0f*
                        (atomJ.labFrameQuadrupole_XZ*fkdg11+atomJ.labFrameQuadrupole_YZ*fkdg21+atomJ.labFrameQuadrupole_ZZ*fkdg31
                        -atomJ.labFrameQuadrupole_XX*fkdg13-atomJ.labFrameQuadrupole_XY*fkdg23-atomJ.labFrameQuadrupole_XZ*fkdg33);

    trqi_k3 += 2.0f*
                        (atomJ.labFrameQuadrupole_XX*fkdg12+atomJ.labFrameQuadrupole_XY*fkdg22+atomJ.labFrameQuadrupole_XZ*fkdg32
                        -atomJ.labFrameQuadrupole_XY*fkdg11-atomJ.labFrameQuadrupole_YY*fkdg21-atomJ.labFrameQuadrupole_YZ*fkdg31);

    // total permanent and induced energies for this interaction

    e                        = esym + 0.5f*(ewi+ewk);
    ei                       = 0.5f * (esymi + 0.5f*(ewii+ewki));

    outputForce[0]           = (dedx + dpdx);
    outputForce[1]           = (dedy + dpdy);
    outputForce[2]           = (dedz + dpdz);

    outputTorque[0][0]       = (trq1 + trqi1);
    outputTorque[0][1]       = (trq2 + trqi2);
    outputTorque[0][2]       = (trq3 + trqi3);

    outputTorque[1][0]       = (trq_k1 + trqi_k1);
    outputTorque[1][1]       = (trq_k2 + trqi_k2);
    outputTorque[1][2]       = (trq_k3 + trqi_k3);

    outputBorn[0]            = drbi;
    outputBorn[1]            = drbk;

    outputBornPolar[0]       = dpbi;
    outputBornPolar[1]       = dpbk;

    *outputEnergy            = (e + ei);

}

#undef SUB_METHOD_NAME
#undef F1
#define SUB_METHOD_NAME(a, b) a##F1##b
#define F1
#include "kCalculateAmoebaCudaKirkwood_b.h"
#undef F1
#undef SUB_METHOD_NAME

#undef F2
#define SUB_METHOD_NAME(a, b) a##F2##b
#define F2
#include "kCalculateAmoebaCudaKirkwood_b.h"
#undef F2
#undef SUB_METHOD_NAME
/*
#undef F1
#undef F2
#define SUB_METHOD_NAME(a, b) a##F1F2##b
#define F1
#define F2
#include "kCalculateAmoebaCudaKirkwood_b.h"
#undef F1
#undef F2
#undef SUB_METHOD_NAME
*/

#undef SUB_METHOD_NAME
#undef T1
#define SUB_METHOD_NAME(a, b) a##T1##b
#define T1
#include "kCalculateAmoebaCudaKirkwood_b.h"
#undef T1
#undef SUB_METHOD_NAME

#undef T2
#define SUB_METHOD_NAME(a, b) a##T2##b
#define T2
#include "kCalculateAmoebaCudaKirkwood_b.h"
#undef T2
#undef SUB_METHOD_NAME

#undef T1
#undef T2
#define SUB_METHOD_NAME(a, b) a##T1T2##b
#define T1
#define T2
#include "kCalculateAmoebaCudaKirkwood_b.h"
#undef T1
#undef T2
#undef SUB_METHOD_NAME

#undef SUB_METHOD_NAME
#undef B1
#define SUB_METHOD_NAME(a, b) a##B1##b
#define B1
#include "kCalculateAmoebaCudaKirkwood_b.h"
#undef B1
#undef SUB_METHOD_NAME

#undef B2
#define SUB_METHOD_NAME(a, b) a##B2##b
#define B2
#include "kCalculateAmoebaCudaKirkwood_b.h"
#undef B2
#undef SUB_METHOD_NAME

#undef B2
#define SUB_METHOD_NAME(a, b) a##B1B2##b
#define B1
#define B2
#include "kCalculateAmoebaCudaKirkwood_b.h"
#undef B1
#undef B2
#undef SUB_METHOD_NAME

#ifdef INCLUDE_TORQUE
__device__ void calculateKirkwoodPairIxn_kernel( KirkwoodParticle& atomI, KirkwoodParticle& atomJ, float forceFactor, float* outputEnergy  ){

//#define Original
#define New
#ifdef Original
    float outputBornPolar[2];
    return calculateKirkwoodPairIxnOrig_kernel( atomI, atomJ,
                                                outputForce, outputTorque, outputBornPolar, outputEnergy );
#endif
#ifdef New


    float force[3];
    float energy;
    calculateKirkwoodPairIxnF1_kernel( atomI, atomJ, &energy, force);
    calculateKirkwoodPairIxnF2_kernel( atomI, atomJ, &energy, force);

    force[0]                *= 0.5f;
    force[1]                *= 0.5f;
    force[2]                *= 0.5f;

    atomI.force[0]          += force[0];
    atomI.force[1]          += force[1];
    atomI.force[2]          += force[2];

    *outputEnergy           += 0.5f*forceFactor*energy;
    if( forceFactor == 1.0f ){
        atomJ.force[0]      -= force[0];
        atomJ.force[1]      -= force[1];
        atomJ.force[2]      -= force[2];
    }    

    calculateKirkwoodPairIxnT1_kernel( atomI, atomJ );
    calculateKirkwoodPairIxnT2_kernel( atomI, atomJ );

    //calculateKirkwoodPairIxnT1T2_kernel( atomI, atomJ );
    if( forceFactor == 1.0f ){
        calculateKirkwoodPairIxnT1_kernel( atomJ, atomI );
        calculateKirkwoodPairIxnT2_kernel( atomJ, atomI );
        //calculateKirkwoodPairIxnT1T2_kernel( atomJ, atomI );
    }


    calculateKirkwoodPairIxnB1B2_kernel( atomI, atomJ );

    return;
#endif
}
#endif

__device__ void zeroKirkwoodParticleSharedField( struct KirkwoodParticle* sA )
{
    // zero shared fields

    sA->force[0]              = 0.0f;
    sA->force[1]              = 0.0f;
    sA->force[2]              = 0.0f;

#ifdef INCLUDE_TORQUE
    sA->torque[0]             = 0.0f;
    sA->torque[1]             = 0.0f;
    sA->torque[2]             = 0.0f;
#endif

    sA->dBornRadius           = 0.0f;
    sA->dBornRadiusPolar      = 0.0f;

}

// Include versions of the kernels for N^2 calculations.

#undef USE_OUTPUT_BUFFER_PER_WARP
#define METHOD_NAME(a, b) a##N2##b
#include "kCalculateAmoebaCudaKirkwood.h"
#define USE_OUTPUT_BUFFER_PER_WARP
#undef METHOD_NAME
#define METHOD_NAME(a, b) a##N2ByWarp##b
#include "kCalculateAmoebaCudaKirkwood.h"

// reduce psWorkArray_1_1 -> dBorn
// reduce psWorkArray_1_2 -> dBornPolar

__global__
#if (__CUDA_ARCH__ >= 200)
__launch_bounds__(GF1XX_THREADS_PER_BLOCK, 1)
#elif (__CUDA_ARCH__ >= 120)
__launch_bounds__(GT2XX_THREADS_PER_BLOCK, 1)
#else
__launch_bounds__(G8X_THREADS_PER_BLOCK, 1)
#endif
 void kReduceToObcBornForcePrefactor_kernel( unsigned int fieldComponents, unsigned int outputBuffers, float* fieldIn1, float* fieldIn2,
                                             float* fieldOut )
{
    unsigned int pos = blockIdx.x * blockDim.x + threadIdx.x;

    // Reduce field

    while (pos < fieldComponents)
    {

        float totalField = 0.0f;

        float* pFt1      = fieldIn1 + pos;
        float* pFt2      = fieldIn2 + pos;

        float bornRadius = cSim.pBornRadii[pos];
        float obcChain   = cSim.pObcChain[pos];

        unsigned int i   = outputBuffers;

        while (i >= 4)
        {
            totalField += pFt1[0] + pFt1[fieldComponents] + pFt1[2*fieldComponents] + pFt1[3*fieldComponents];
            totalField += pFt2[0] + pFt2[fieldComponents] + pFt2[2*fieldComponents] + pFt2[3*fieldComponents];
            pFt1       += fieldComponents*4;
            pFt2       += fieldComponents*4;
            i          -= 4;
        }

        if (i >= 2)
        {
            totalField += pFt1[0] + pFt1[fieldComponents];
            totalField += pFt2[0] + pFt2[fieldComponents];
            pFt1       += fieldComponents*2;
            pFt2       += fieldComponents*2;
            i          -= 2;
        }

        if (i > 0)
        {
            totalField += pFt1[0];
            totalField += pFt2[0];
        }

        fieldOut[pos]   = totalField*bornRadius*bornRadius*obcChain;
        pos            += gridDim.x * blockDim.x;
    }
}

__global__
#if (__CUDA_ARCH__ >= 200)
__launch_bounds__(GF1XX_THREADS_PER_BLOCK, 1)
#elif (__CUDA_ARCH__ >= 120)
__launch_bounds__(GT2XX_THREADS_PER_BLOCK, 1)
#else
__launch_bounds__(G8X_THREADS_PER_BLOCK, 1)
#endif
void kReduceToObcBornForcePrefactorAndSASA_kernel( unsigned int fieldComponents, unsigned int outputBuffers, float* fieldIn1, float* fieldIn2,
                                                   float* fieldOut )
{
    unsigned int pos = blockIdx.x * blockDim.x + threadIdx.x;

    float energy     = 0.0f;

    // Reduce field

    while (pos < fieldComponents)
    {

        float totalForce = 0.0f;

        float* pFt1      = fieldIn1 + pos;
        float* pFt2      = fieldIn2 + pos;

        float bornRadius = cSim.pBornRadii[pos];
        float obcChain   = cSim.pObcChain[pos];
        float2 obcData   = cSim.pObcData[pos];

        unsigned int i   = outputBuffers;

        while (i >= 4)
        {
            totalForce += pFt1[0] + pFt1[fieldComponents] + pFt1[2*fieldComponents] + pFt1[3*fieldComponents];
            totalForce += pFt2[0] + pFt2[fieldComponents] + pFt2[2*fieldComponents] + pFt2[3*fieldComponents];
            pFt1       += fieldComponents*4;
            pFt2       += fieldComponents*4;
            i          -= 4;
        }

        if (i >= 2)
        {
            totalForce += pFt1[0] + pFt1[fieldComponents];
            totalForce += pFt2[0] + pFt2[fieldComponents];
            pFt1       += fieldComponents*2;
            pFt2       += fieldComponents*2;
            i          -= 2;
        }

        if (i > 0)
        {
            totalForce += pFt1[0];
            totalForce += pFt2[0];
        }

        float r        = (obcData.x + cSim.dielectricOffset + cSim.probeRadius);
        float ratio6   = ( (obcData.x + cSim.dielectricOffset) / bornRadius);
              ratio6   = ratio6*ratio6*ratio6;
              ratio6   = ratio6*ratio6;
        float saTerm   = cSim.surfaceAreaFactor * r * r * ratio6;

        totalForce    += saTerm / bornRadius;
        totalForce    *= bornRadius * bornRadius * obcChain;

        fieldOut[pos]  = totalForce;

        energy        += saTerm;
        pos           += gridDim.x * blockDim.x;
    }

    cSim.pEnergy[blockIdx.x * blockDim.x + threadIdx.x] += energy / -6.0f;
}


__global__
#if (__CUDA_ARCH__ >= 200)
__launch_bounds__(GF1XX_THREADS_PER_BLOCK, 1)
#elif (__CUDA_ARCH__ >= 120)
__launch_bounds__(GT2XX_THREADS_PER_BLOCK, 1)
#else
__launch_bounds__(G8X_THREADS_PER_BLOCK, 1)
#endif
 void kReduceToBornForcePrefactor_kernel( unsigned int fieldComponents, unsigned int outputBuffers, float* fieldIn1, float* fieldIn2,
                                          float* fieldOut )
{
    unsigned int pos = blockIdx.x * blockDim.x + threadIdx.x;

    // Reduce field

    while (pos < fieldComponents)
    {

        float totalField = 0.0f;

        float* pFt1      = fieldIn1 + pos;
        float* pFt2      = fieldIn2 + pos;

        //float bornRadius = cSim.pBornRadii[pos];
        //float obcChain   = cSim.pObcChain[pos];

        unsigned int i   = outputBuffers;

        while (i >= 4)
        {
            totalField += pFt1[0] + pFt1[fieldComponents] + pFt1[2*fieldComponents] + pFt1[3*fieldComponents];
            totalField += pFt2[0] + pFt2[fieldComponents] + pFt2[2*fieldComponents] + pFt2[3*fieldComponents];
            pFt1       += fieldComponents*4;
            pFt2       += fieldComponents*4;
            i          -= 4;
        }

        if (i >= 2)
        {
            totalField += pFt1[0] + pFt1[fieldComponents];
            totalField += pFt2[0] + pFt2[fieldComponents];
            pFt1       += fieldComponents*2;
            pFt2       += fieldComponents*2;
            i          -= 2;
        }

        if (i > 0)
        {
            totalField += pFt1[0];
            totalField += pFt2[0];
        }

        //fieldOut[pos]   = totalField*bornRadius*bornRadius*obcChain;
        fieldOut[pos]   = totalField;
        pos            += gridDim.x * blockDim.x;
    }
}

__global__
#if (__CUDA_ARCH__ >= 200)
__launch_bounds__(GF1XX_THREADS_PER_BLOCK, 1)
#elif (__CUDA_ARCH__ >= 120)
__launch_bounds__(GT2XX_THREADS_PER_BLOCK, 1)
#else
__launch_bounds__(G8X_THREADS_PER_BLOCK, 1)
#endif
void kReduceToBornForcePrefactorAndSASA_kernel( unsigned int fieldComponents, unsigned int outputBuffers, float* fieldIn1, float* fieldIn2,
                                                float* fieldOut )
{
    unsigned int pos = blockIdx.x * blockDim.x + threadIdx.x;

    float energy     = 0.0f;

    // Reduce field

    while (pos < fieldComponents)
    {

        float totalForce = 0.0f;

        float* pFt1      = fieldIn1 + pos;
        float* pFt2      = fieldIn2 + pos;

        float bornRadius = cSim.pBornRadii[pos];
        //float obcChain   = cSim.pObcChain[pos];
        float2 obcData   = cSim.pObcData[pos];

        unsigned int i   = outputBuffers;

        while (i >= 4)
        {
            totalForce += pFt1[0] + pFt1[fieldComponents] + pFt1[2*fieldComponents] + pFt1[3*fieldComponents];
            totalForce += pFt2[0] + pFt2[fieldComponents] + pFt2[2*fieldComponents] + pFt2[3*fieldComponents];
            pFt1       += fieldComponents*4;
            pFt2       += fieldComponents*4;
            i          -= 4;
        }

        if (i >= 2)
        {
            totalForce += pFt1[0] + pFt1[fieldComponents];
            totalForce += pFt2[0] + pFt2[fieldComponents];
            pFt1       += fieldComponents*2;
            pFt2       += fieldComponents*2;
            i          -= 2;
        }

        if (i > 0)
        {
            totalForce += pFt1[0];
            totalForce += pFt2[0];
        }

        float r        = (obcData.x + cSim.dielectricOffset + cSim.probeRadius);
        float ratio6   = ( (obcData.x + cSim.dielectricOffset) / bornRadius);
              ratio6   = ratio6*ratio6*ratio6;
              ratio6   = ratio6*ratio6;
        float saTerm   = cSim.surfaceAreaFactor * r * r * ratio6;

        totalForce    += saTerm / bornRadius;
        //totalForce    *= bornRadius * bornRadius * obcChain;

        fieldOut[pos]  = totalForce;

        energy        += saTerm;
        pos           += gridDim.x * blockDim.x;
    }

    cSim.pEnergy[blockIdx.x * blockDim.x + threadIdx.x] += energy / -6.0f;
}

static void kReduceToBornForcePrefactor( amoebaGpuContext amoebaGpu )
{

    gpuContext gpu = amoebaGpu->gpuContext;
    if( amoebaGpu->includeObcCavityTerm ){

        kReduceToBornForcePrefactorAndSASA_kernel<<<gpu->sim.nonbond_blocks, gpu->sim.bsf_reduce_threads_per_block>>>(
                                                    gpu->sim.paddedNumberOfAtoms, gpu->sim.outputBuffers,
                                                    amoebaGpu->psWorkArray_1_1->_pDevData,
                                                    amoebaGpu->psWorkArray_1_2->_pDevData,
                                                    amoebaGpu->gpuContext->psBornForce->_pDevData );
    } else {

        kReduceToBornForcePrefactor_kernel<<<gpu->sim.nonbond_blocks, gpu->sim.bsf_reduce_threads_per_block>>>(
                                             gpu->sim.paddedNumberOfAtoms, gpu->sim.outputBuffers,
                                             amoebaGpu->psWorkArray_1_1->_pDevData,
                                             amoebaGpu->psWorkArray_1_2->_pDevData,
                                             amoebaGpu->gpuContext->psBornForce->_pDevData );
    }
    LAUNCHERROR("kReduceToBornForcePrefactor");
}

/**---------------------------------------------------------------------------------------

   Compute Amoeba electrostatic force & torque

   @param amoebaGpu        amoebaGpu context
   @param gpu              OpenMM gpu Cuda context

   --------------------------------------------------------------------------------------- */

void kCalculateAmoebaKirkwood( amoebaGpuContext amoebaGpu )
{

    // ---------------------------------------------------------------------------------------

    gpuContext gpu = amoebaGpu->gpuContext;

    // apparently debug array can take up nontrivial no. registers
    // on first pass, set threads/block and based on that setting the energy buffer array

    static unsigned int threadsPerBlock = 0;
    if( threadsPerBlock == 0 ){
        unsigned int maxThreads;
        if (gpu->sm_version >= SM_20)
            //maxThreads = 384;
            maxThreads = 512;
        else if (gpu->sm_version >= SM_12)
            maxThreads = 128;
        else
            maxThreads = 64;
        threadsPerBlock = std::min(getThreadsPerBlock(amoebaGpu, sizeof(KirkwoodParticle), gpu->sharedMemoryPerBlock ), maxThreads);

    }

    kClearFields_1( amoebaGpu );

    if (gpu->bOutputBufferPerWarp){
        kCalculateAmoebaCudaKirkwoodN2ByWarpForces_kernel<<<gpu->sim.nonbond_blocks, threadsPerBlock, sizeof(KirkwoodParticle)*threadsPerBlock>>>(
                                                                           gpu->psWorkUnit->_pDevData);
    } else {

        kCalculateAmoebaCudaKirkwoodN2Forces_kernel<<<gpu->sim.nonbond_blocks, threadsPerBlock, sizeof(KirkwoodParticle)*threadsPerBlock>>>(
                                                                           gpu->psWorkUnit->_pDevData);
    }
    LAUNCHERROR("kCalculateAmoebaCudaKirkwoodN2Forces");

    if( 0 ){
        gpuContext gpu                       = amoebaGpu->gpuContext;
        std::vector<int> fileId;
        fileId.push_back( 0 );
        VectorOfDoubleVectors outputVector;
        //cudaLoadCudaFloat4Array( gpu->natoms, 3, gpu->psPosq4,              outputVector, gpu->psAtomIndex->_pSysData, 1.0f );
        //cudaLoadCudaFloatArray( gpu->natoms,  3, amoebaGpu->psLabFrameDipole,     outputVector, gpu->psAtomIndex->_pSysData, 1.0f );
        CUDAStream<float>* temp  = new CUDAStream<float>(3*gpu->sim.paddedNumberOfAtoms, 1, "Temp1");

        reduceAndCopyCUDAStreamFloat( amoebaGpu->psWorkArray_3_1, temp, 1.0 );
        cudaLoadCudaFloatArray( gpu->natoms,  3, temp, outputVector, gpu->psAtomIndex->_pSysData, 1.0f ); 
        cudaWriteVectorOfDoubleVectorsToFile( "Torq", fileId, outputVector );
        delete temp;
    }   
    kReduceToBornForcePrefactor( amoebaGpu );

    // Tinker's Born1 && E-diff

    //kCalculateObcGbsaForces2( amoebaGpu->gpuContext );
    kCalculateGrycukGbsaForces2( amoebaGpu );
    kCalculateAmoebaKirkwoodEDiff( amoebaGpu );

   // ---------------------------------------------------------------------------------------
}