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__ static void loadFixedFieldShared( struct FixedFieldParticle* sA, unsigned int atomI 
#ifdef GK
    , float* bornR
#endif
)
{
    // 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_X         = cAmoebaSim.pLabFrameDipole[atomI*3];
    sA->labFrameDipole_Y         = cAmoebaSim.pLabFrameDipole[atomI*3+1];
    sA->labFrameDipole_Z         = 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->damp                     = cAmoebaSim.pDampingFactorAndThole[atomI].x;
    sA->thole                    = cAmoebaSim.pDampingFactorAndThole[atomI].y;
#ifdef GK
    sA->bornR                    = bornR[atomI];
#endif

}

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

__device__ static 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__ static 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__ static void calculateFixedEFieldPairIxn_kernel( FixedFieldParticle& atomI, FixedFieldParticle& atomJ,
                                                           float field[2][3]
#ifdef AMOEBA_DEBUG
                                                    , float4 debugArray[12]
#endif
)
{

 
    // ---------------------------------------------------------------------------------------
 
    // get deltaR and r between 2 atoms
 
    float deltaR[3];
    deltaR[0]           = atomJ.x - atomI.x;
    deltaR[1]           = atomJ.y - atomI.y;
    deltaR[2]           = atomJ.z - atomI.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          = atomI.damp*atomJ.damp;
    float dampExp;
    if( damp != 0.0f && r < cAmoebaSim.scalingDistanceCutoff ){

        // get scaling factors
      
        float ratio     = r/damp;
        float pGamma    = atomJ.thole > atomI.thole ? atomI.thole : atomJ.thole; 
        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;
 
    float  qDotDelta[3];
    qDotDelta[0]        = deltaR[0]*atomJ.labFrameQuadrupole_XX + deltaR[1]*atomJ.labFrameQuadrupole_XY + deltaR[2]*atomJ.labFrameQuadrupole_XZ;
    qDotDelta[1]        = deltaR[0]*atomJ.labFrameQuadrupole_XY + deltaR[1]*atomJ.labFrameQuadrupole_YY + deltaR[2]*atomJ.labFrameQuadrupole_YZ;
    qDotDelta[2]        = deltaR[0]*atomJ.labFrameQuadrupole_XZ + deltaR[1]*atomJ.labFrameQuadrupole_YZ + deltaR[2]*atomJ.labFrameQuadrupole_ZZ;
 
    float dotdd         = deltaR[0]*atomJ.labFrameDipole_X      + deltaR[1]*atomJ.labFrameDipole_Y      + deltaR[2]*atomJ.labFrameDipole_Z;
    float dotqd         = deltaR[0]*qDotDelta[0]                + deltaR[1]*qDotDelta[1]                + deltaR[2]*qDotDelta[2];

    float factor        = -rr3*atomJ.q + rr5*dotdd - rr7*dotqd;
 
    field[0][0]         = deltaR[0]*factor - rr3*atomJ.labFrameDipole_X + rr5_2*qDotDelta[0];
    field[0][1]         = deltaR[1]*factor - rr3*atomJ.labFrameDipole_Y + rr5_2*qDotDelta[1];
    field[0][2]         = deltaR[2]*factor - rr3*atomJ.labFrameDipole_Z + rr5_2*qDotDelta[2];
 
    qDotDelta[0]        = deltaR[0]*atomI.labFrameQuadrupole_XX + deltaR[1]*atomI.labFrameQuadrupole_XY + deltaR[2]*atomI.labFrameQuadrupole_XZ;
    qDotDelta[1]        = deltaR[0]*atomI.labFrameQuadrupole_XY + deltaR[1]*atomI.labFrameQuadrupole_YY + deltaR[2]*atomI.labFrameQuadrupole_YZ;
    qDotDelta[2]        = deltaR[0]*atomI.labFrameQuadrupole_XZ + deltaR[1]*atomI.labFrameQuadrupole_YZ + deltaR[2]*atomI.labFrameQuadrupole_ZZ;
 
    dotdd               = deltaR[0]*atomI.labFrameDipole_X    + deltaR[1]*atomI.labFrameDipole_Y    + deltaR[2]*atomI.labFrameDipole_Z;
    dotqd               = deltaR[0]*qDotDelta[0] + deltaR[1]*qDotDelta[1] + deltaR[2]*qDotDelta[2];
    factor              = rr3*atomI.q + rr5*dotdd + rr7*dotqd;
 
    field[1][0]         = deltaR[0]*factor - rr3*atomI.labFrameDipole_X - rr5_2*qDotDelta[0];
    field[1][1]         = deltaR[1]*factor - rr3*atomI.labFrameDipole_Y - rr5_2*qDotDelta[1];
    field[1][2]         = deltaR[2]*factor - rr3*atomI.labFrameDipole_Z - rr5_2*qDotDelta[2];

}