kCalculateAmoebaCudaFixedFieldParticle.h

struct FixedFieldParticle {

    // coordinates charge

    float x;
    float y;
    float z;
    float q;

    // lab frame dipole

    float labFrameDipole_X;
    float labFrameDipole_Y;
    float labFrameDipole_Z;

    // lab frame quadrupole

    float labFrameQuadrupole_XX;
    float labFrameQuadrupole_XY;
    float labFrameQuadrupole_XZ;
    float labFrameQuadrupole_YY;
    float labFrameQuadrupole_YZ;
    float labFrameQuadrupole_ZZ;

    // scaling factor

    float thole;
    float damp;

    // field accumulators

    float eField[3];
    float eFieldP[3];

#ifdef GK

    // Born radius

    float bornR;

    // GK field

    float gkField[3];

#endif
};

__device__ void loadFixedFieldShared( struct FixedFieldParticle* sA, unsigned int atomI,
                                      float4* atomCoord, float* labDipole, float* labQuadrupole,
                                      float2* dampingFactorAndThole
#ifdef GK
    , float* bornR
#endif
)
{
    // coordinates & charge

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

    // lab dipole

    sA->labFrameDipole_X         = labDipole[atomI*3];
    sA->labFrameDipole_Y         = labDipole[atomI*3+1];
    sA->labFrameDipole_Z         = labDipole[atomI*3+2];

    // lab quadrupole

    sA->labFrameQuadrupole_XX    = labQuadrupole[atomI*9];
    sA->labFrameQuadrupole_XY    = labQuadrupole[atomI*9+1];
    sA->labFrameQuadrupole_XZ    = labQuadrupole[atomI*9+2];
    sA->labFrameQuadrupole_YY    = labQuadrupole[atomI*9+4];
    sA->labFrameQuadrupole_YZ    = labQuadrupole[atomI*9+5];
    sA->labFrameQuadrupole_ZZ    = labQuadrupole[atomI*9+8];

    sA->damp                     = dampingFactorAndThole[atomI].x;
    sA->thole                    = dampingFactorAndThole[atomI].y;
#ifdef GK
    sA->bornR                    = bornR[atomI];
#endif

}

// load struct and arrays w/ shared data in sA

__device__ void loadFixedFieldParticleData( struct FixedFieldParticle* sA, 
                                            float4* jCoord, float* jDipole, float* jQuadrupole
#ifdef GK
, float* bornR
#endif
)
{

    // load coords, charge, ...

    jCoord->x               = sA->x;
    jCoord->y               = sA->y;
    jCoord->z               = sA->z;
    jCoord->w               = sA->q;

    jDipole[0]              = sA->labFrameDipole_X;
    jDipole[1]              = sA->labFrameDipole_Y;
    jDipole[2]              = sA->labFrameDipole_Z;

    jQuadrupole[0]          = sA->labFrameQuadrupole_XX;
    jQuadrupole[1]          = sA->labFrameQuadrupole_XY;
    jQuadrupole[2]          = sA->labFrameQuadrupole_XZ;

    jQuadrupole[3]          = sA->labFrameQuadrupole_XY;
    jQuadrupole[4]          = sA->labFrameQuadrupole_YY;
    jQuadrupole[5]          = sA->labFrameQuadrupole_YZ;

    jQuadrupole[6]          = sA->labFrameQuadrupole_XZ;
    jQuadrupole[7]          = sA->labFrameQuadrupole_YZ;
    jQuadrupole[8]          = sA->labFrameQuadrupole_ZZ;
 
#ifdef GK
    *bornR                  = sA->bornR;
#endif
}

// zero fields

__device__ void zeroFixedFieldParticleSharedField( struct FixedFieldParticle* sA ) 
{

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

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

#ifdef GK
    sA->gkField[0]   = 0.0f;
    sA->gkField[1]   = 0.0f;
    sA->gkField[2]   = 0.0f;
#endif
}

// body of fixed E-field calculation

__device__ void calculateFixedEFieldPairIxn_kernel( float4 atomCoordinatesI, float4 atomCoordinatesJ,
                                                    float dampingFactorI,    float dampingFactorJ,
                                                    float tholeI,            float tholeJ,
                                                    float* labDipoleI, float* labDipoleJ,
                                                    float* labQuadrupoleI, float* labQuadrupoleJ,
                                                    float scalingDistanceCutoff,
                                                    float field[2][3]
#ifdef AMOEBA_DEBUG
                                                    , float4 debugArray[12]
#endif
)
{

 
    // ---------------------------------------------------------------------------------------
 
    // get deltaR and r between 2 atoms
 
    float deltaR[3];
    deltaR[0]           = atomCoordinatesJ.x - atomCoordinatesI.x;
    deltaR[1]           = atomCoordinatesJ.y - atomCoordinatesI.y;
    deltaR[2]           = atomCoordinatesJ.z - atomCoordinatesI.z;

    float r             =  SQRT( deltaR[0]*deltaR[0] + deltaR[1]*deltaR[1] + deltaR[2]*deltaR[2] );
    float rI            =  1.0f/r;
    float r2I           =  rI*rI;

    float rr3           =  rI*r2I;
    float rr5           =  3.0f*rr3*r2I;
    float rr7           =  5.0f*rr5*r2I;
 
    // get scaling factors, if needed
    
    float damp          = dampingFactorI*dampingFactorJ;
    float dampExp;
    if( damp != 0.0f && r < scalingDistanceCutoff ){

        // get scaling factors
      
        float ratio     = r/damp;
        float pGamma    = tholeJ > tholeI ? tholeI : tholeJ; 
        damp            = ratio*ratio*ratio*pGamma;
        dampExp         = EXP( -damp );
    } else {
        dampExp         = 0.0f;
    }
      
    rr3                *= 1.0f - dampExp;
    rr5                *= 1.0f - ( 1.0f + damp )*dampExp;
    rr7                *= 1.0f - ( 1.0f + damp + (0.6f*damp*damp))*dampExp;
      
    float rr5_2         = rr5*2.0f;
 
#ifdef AMOEBA_DEBUG
    int index           = 0;

    // 0-2
    debugArray[index].x   = r;
    debugArray[index].y   = rr3;
    debugArray[index].z   = rr5;
    index++;
#endif

    float* dipole       = labDipoleJ;
    float* quadrupole   = labQuadrupoleJ;
    float  qDotDelta[3];
    qDotDelta[0]        = deltaR[0]*quadrupole[0] + deltaR[1]*quadrupole[1] + deltaR[2]*quadrupole[2];
    qDotDelta[1]        = deltaR[0]*quadrupole[3] + deltaR[1]*quadrupole[4] + deltaR[2]*quadrupole[5];
    qDotDelta[2]        = deltaR[0]*quadrupole[6] + deltaR[1]*quadrupole[7] + deltaR[2]*quadrupole[8];
 
    float dotdd         = deltaR[0]*dipole[0]    + deltaR[1]*dipole[1]    + deltaR[2]*dipole[2];
    float dotqd         = deltaR[0]*qDotDelta[0] + deltaR[1]*qDotDelta[1] + deltaR[2]*qDotDelta[2];
    float factor        = -rr3*atomCoordinatesJ.w + rr5*dotdd - rr7*dotqd;
 
#ifdef AMOEBA_DEBUG
    // 3-5
    debugArray[index].x   = dotdd;
    debugArray[index].y   = dotqd;
    debugArray[index].z   = factor;
    index++;
#endif

    field[0][0]         = deltaR[0]*factor - rr3*dipole[0] + rr5_2*qDotDelta[0];
    field[0][1]         = deltaR[1]*factor - rr3*dipole[1] + rr5_2*qDotDelta[1];
    field[0][2]         = deltaR[2]*factor - rr3*dipole[2] + rr5_2*qDotDelta[2];
 
    dipole              = labDipoleI;
    quadrupole          = labQuadrupoleI;
    qDotDelta[0]        = deltaR[0]*quadrupole[0] + deltaR[1]*quadrupole[1] + deltaR[2]*quadrupole[2];
    qDotDelta[1]        = deltaR[0]*quadrupole[3] + deltaR[1]*quadrupole[4] + deltaR[2]*quadrupole[5];
    qDotDelta[2]        = deltaR[0]*quadrupole[6] + deltaR[1]*quadrupole[7] + deltaR[2]*quadrupole[8];
 
    dotdd               = deltaR[0]*dipole[0]    + deltaR[1]*dipole[1]    + deltaR[2]*dipole[2];
    dotqd               = deltaR[0]*qDotDelta[0] + deltaR[1]*qDotDelta[1] + deltaR[2]*qDotDelta[2];
    factor              = rr3*atomCoordinatesI.w + rr5*dotdd + rr7*dotqd;
 
#ifdef AMOEBA_DEBUG
    // 6-8
    debugArray[index].x = dotdd;
    debugArray[index].y = dotqd;
    debugArray[index].z = factor;
    index++;
#endif

    field[1][0]         = deltaR[0]*factor - rr3*dipole[0] - rr5_2*qDotDelta[0];
    field[1][1]         = deltaR[1]*factor - rr3*dipole[1] - rr5_2*qDotDelta[1];
    field[1][2]         = deltaR[2]*factor - rr3*dipole[2] - rr5_2*qDotDelta[2];

#if 0
    float testValue  = 1.0f;
    field[0][0]      = testValue;
    field[0][1]      = testValue;
    field[0][2]      = testValue;
    field[1][0]      = testValue;
    field[1][1]      = testValue;
    field[1][2]      = testValue;
#endif
 
}