
static __device__ void SUB_METHOD_NAME( calculatePmeDirectElectrostaticPairIxnT2, _kernel)( 
                                        PmeDirectElectrostaticParticle& atomI, PmeDirectElectrostaticParticle& atomJ,
                                        const float4 delta, const float4 bn
#ifdef APPLY_SCALE
                                        , const float* scalingFactors
#endif
                                        ){

    float xr                    = delta.x;
    float yr                    = delta.y;
    float zr                    = delta.z;
    float rr1                   = delta.w;

    // set the permanent multipole and induced dipole values;

    float di1                   = atomI.labFrameDipole[0];
    float di2                   = atomI.labFrameDipole[1];
    float di3                   = atomI.labFrameDipole[2];

    float qi1                   = atomI.labFrameQuadrupole[0];
    float qi2                   = atomI.labFrameQuadrupole[1];
    float qi3                   = atomI.labFrameQuadrupole[2];
    float qi5                   = atomI.labFrameQuadrupole[3];
    float qi6                   = atomI.labFrameQuadrupole[4];
    //float qi9                   = atomI.labFrameQuadrupole[5];
    float qi9                   = -(atomI.labFrameQuadrupole[0] + atomI.labFrameQuadrupole[3]);

    float bn1                   = bn.x;
    float bn2                   = bn.y;
    float bn3                   = bn.z;

    // apply Thole polarization damping to scale factors

    float scale3                = 1.0f;
    float scale5                = 1.0f;
    float scale7                = 1.0f;

    float damp                  = atomI.damp*atomJ.damp;
    if( damp != 0.0f ){
        float pgamma  = atomI.thole < atomJ.thole ? atomI.thole : atomJ.thole;
        float ratio   = 1.0f/(rr1*damp);
            damp      = -pgamma*ratio*ratio*ratio;
        if( damp > -50.0f ){
            float expdamp    = exp(damp);
            scale3           = 1.0f - expdamp;
            scale5           = 1.0f - (1.0f-damp)*expdamp;
            scale7           = 1.0f - (1.0f-damp+0.6f*damp*damp)*expdamp;
        }
    }


    float rr3                   = rr1*rr1*rr1;
#ifdef APPLY_SCALE
    float dsc3                  = rr3*(1.0f - scale3*scalingFactors[DScaleIndex]);
    float dsc5                  = (3.0f*rr3*rr1*rr1)* (1.0f - scale5*scalingFactors[DScaleIndex]);
    float dsc7                  = (15.0f*rr3*rr3*rr1)*(1.0f - scale7*scalingFactors[DScaleIndex]);

    float psc3                  = rr3*(1.0f - scale3*scalingFactors[PScaleIndex]);
    float psc5                  = (3.0f*rr3*rr1*rr1)*(1.0f - scale5*scalingFactors[PScaleIndex]);
    float psc7                  = (15.0f*rr3*rr3*rr1)*(1.0f - scale7*scalingFactors[PScaleIndex]);
#else
    float psc3                  = rr3*(1.0f - scale3);
    float psc5                  = (3.0f*rr3*rr1*rr1)*(1.0f - scale5);
    float psc7                  = (15.0f*rr3*rr3*rr1)*(1.0f - scale7);
#endif

    float prefactor1            = 0.5f*(psc3 - bn1);
#ifdef APPLY_SCALE
    float prefactor2            = 0.5f*(dsc3 - bn1);
#endif

    float dixuk1                = di2*atomJ.inducedDipole[2]  - di3*atomJ.inducedDipole[1];
    float dixukp1               = di2*atomJ.inducedDipoleP[2] - di3*atomJ.inducedDipoleP[1];

#ifdef APPLY_SCALE
    float ttm2i1                = prefactor1*dixuk1 + prefactor2*dixukp1;
#else
    float ttm2i1                = prefactor1*(dixuk1 + dixukp1);
#endif

    float dixuk2                = di3*atomJ.inducedDipole[0]  - di1*atomJ.inducedDipole[2];
    float dixukp2               = di3*atomJ.inducedDipoleP[0] - di1*atomJ.inducedDipoleP[2];

#ifdef APPLY_SCALE
    float ttm2i2                = prefactor1*dixuk2 + prefactor2*dixukp2;
#else
    float ttm2i2                = prefactor1*(dixuk2 + dixukp2);
#endif

    float dixuk3                = di1*atomJ.inducedDipole[1]  - di2*atomJ.inducedDipole[0];
    float dixukp3               = di1*atomJ.inducedDipoleP[1] - di2*atomJ.inducedDipoleP[0];
#ifdef APPLY_SCALE
    float ttm2i3                = prefactor1*dixuk3 + prefactor2*dixukp3;
#else
    float ttm2i3                = prefactor1*(dixuk3 + dixukp3);
#endif

    float sci4                  = atomJ.inducedDipole[0]*xr  + atomJ.inducedDipole[1]*yr  + atomJ.inducedDipole[2]*zr;
    float scip4                 = atomJ.inducedDipoleP[0]*xr + atomJ.inducedDipoleP[1]*yr + atomJ.inducedDipoleP[2]*zr;
    float gti2                  = bn2*(sci4+scip4);
#ifdef APPLY_SCALE
    float gtri2                 = (sci4*psc5+scip4*dsc5);
#else
    float gtri2                 = psc5*(sci4+scip4);
#endif
    prefactor1                  = 0.5f*(gti2 - gtri2);

    ttm2i1                     += prefactor1*( di2*zr - di3*yr );
    ttm2i2                     += prefactor1*( di3*xr - di1*zr );
    ttm2i3                     += prefactor1*( di1*yr - di2*xr );

    float qir1                  = qi1*xr + qi2*yr + qi3*zr;
    float qir2                  = qi2*xr + qi5*yr + qi6*zr;
    float qir3                  = qi3*xr + qi6*yr + qi9*zr;

#ifdef APPLY_SCALE
    prefactor1                  = sci4*psc7 + scip4*dsc7 - bn3*(sci4+scip4);
#else
    prefactor1                  = psc7*(sci4+scip4) - bn3*(sci4+scip4);
#endif
    ttm2i1                     += prefactor1*( yr*qir3 - zr*qir2 );
    ttm2i2                     += prefactor1*( zr*qir1 - xr*qir3 );
    ttm2i3                     += prefactor1*( xr*qir2 - yr*qir1 );

    float qiuk1                 = qi1*atomJ.inducedDipole[0]  + qi2*atomJ.inducedDipole[1]  + qi3*atomJ.inducedDipole[2];
    float qiuk2                 = qi2*atomJ.inducedDipole[0]  + qi5*atomJ.inducedDipole[1]  + qi6*atomJ.inducedDipole[2];
    float qiuk3                 = qi3*atomJ.inducedDipole[0]  + qi6*atomJ.inducedDipole[1]  + qi9*atomJ.inducedDipole[2];

    float qiukp1                = qi1*atomJ.inducedDipoleP[0] + qi2*atomJ.inducedDipoleP[1] + qi3*atomJ.inducedDipoleP[2];
    float qiukp2                = qi2*atomJ.inducedDipoleP[0] + qi5*atomJ.inducedDipoleP[1] + qi6*atomJ.inducedDipoleP[2];
    float qiukp3                = qi3*atomJ.inducedDipoleP[0] + qi6*atomJ.inducedDipoleP[1] + qi9*atomJ.inducedDipoleP[2];

    prefactor1                  = (bn2 - psc5);
#ifdef APPLY_SCALE
    prefactor2                  = (bn2 - dsc5);
#endif
    float ukxqir1               = atomJ.inducedDipole[1]*qir3  - atomJ.inducedDipole[2]*qir2;
    float ukxqirp1              = atomJ.inducedDipoleP[1]*qir3 - atomJ.inducedDipoleP[2]*qir2;
    float rxqiuk1               = yr*qiuk3  - zr*qiuk2;
    float rxqiukp1              = yr*qiukp3 - zr*qiukp2;

#ifdef APPLY_SCALE
    ttm2i1                     += prefactor1*(ukxqir1 + rxqiuk1) + prefactor2*(ukxqirp1 + rxqiukp1);
#else
    ttm2i1                     += prefactor1*( ukxqir1 + rxqiuk1 + ukxqirp1 + rxqiukp1 );
#endif

    float ukxqir2               = atomJ.inducedDipole[2]*qir1  - atomJ.inducedDipole[0]*qir3;
    float ukxqirp2              = atomJ.inducedDipoleP[2]*qir1 - atomJ.inducedDipoleP[0]*qir3;
    float rxqiuk2               = zr*qiuk1  - xr*qiuk3;
    float rxqiukp2              = zr*qiukp1 - xr*qiukp3;
#ifdef APPLY_SCALE
    ttm2i2                     += prefactor1*(ukxqir2 + rxqiuk2) + prefactor2*(ukxqirp2 + rxqiukp2);
#else
    ttm2i2                     += prefactor1*( ukxqir2 + rxqiuk2 + ukxqirp2 + rxqiukp2 );
#endif

    float ukxqir3               = atomJ.inducedDipole[0]*qir2  - atomJ.inducedDipole[1]*qir1;
    float ukxqirp3              = atomJ.inducedDipoleP[0]*qir2 - atomJ.inducedDipoleP[1]*qir1;
    float rxqiuk3               = xr*qiuk2  - yr*qiuk1;
    float rxqiukp3              = xr*qiukp2 - yr*qiukp1;
#ifdef APPLY_SCALE
    ttm2i3                     += prefactor1*(ukxqir3 + rxqiuk3) + prefactor2*(ukxqirp3 + rxqiukp3);
#else
    ttm2i3                     += prefactor1*(ukxqir3 + rxqiuk3 + ukxqirp3 + rxqiukp3 );
#endif

    atomI.torque[0]            += ttm2i1;
    atomI.torque[1]            += ttm2i2;
    atomI.torque[2]            += ttm2i3;

/*
    torque[0]            += ttm2i1;
    torque[1]            += ttm2i2;
    torque[2]            += ttm2i3;
*/

    return;

}
