ComputeNonbondedCUDAExcl.C

Go to the documentation of this file.
00001 
00007 #include "InfoStream.h"
00008 #include "ComputeNonbondedCUDAExcl.h"
00009 #include "Molecule.h"
00010 #include "Parameters.h"
00011 #include "LJTable.h"
00012 #include "Node.h"
00013 #include "ReductionMgr.h"
00014 #include "Lattice.h"
00015 #include "PressureProfile.h"
00016 #include "Debug.h"
00017 
00018 
00019 // static initialization
00020 int ExclElem::pressureProfileSlabs = 0;
00021 int ExclElem::pressureProfileAtomTypes = 1;
00022 BigReal ExclElem::pressureProfileThickness = 0;
00023 BigReal ExclElem::pressureProfileMin = 0;
00024 
00025 void ExclElem::getMoleculePointers
00026     (Molecule* mol, int* count, int32*** byatom, Exclusion** structarray)
00027 {
00028 #ifdef MEM_OPT_VERSION
00029   NAMD_die("Should not be called in ExclElem::getMoleculePointers in memory optimized version!");
00030 #else
00031   *count = mol->numExclusions;
00032   *byatom = mol->exclusionsByAtom;
00033   *structarray = mol->exclusions;
00034 #endif
00035 }
00036 
00037 void ExclElem::getParameterPointers(Parameters *p, const int **v) {
00038   *v = 0;
00039 }
00040 
00041 void ExclElem::computeForce(ExclElem *tuples, int ntuple, BigReal *reduction, 
00042                             BigReal *pressureProfileData)
00043 {
00044   const Lattice & lattice = tuples[0].p[0]->p->lattice;
00045   const Flags &flags = tuples[0].p[0]->p->flags;
00046   if ( ! flags.doNonbonded ) return;
00047   const int doFull = flags.doFullElectrostatics;
00048   const int doEnergy = flags.doEnergy;
00049 
00050   for ( int ituple=0; ituple<ntuple; ++ituple ) {
00051     const ExclElem &tup = tuples[ituple];
00052     enum { size = 2 };
00053     const AtomID (&atomID)[size](tup.atomID);
00054     const int    (&localIndex)[size](tup.localIndex);
00055     TuplePatchElem * const(&p)[size](tup.p);
00056     const Real (&scale)(tup.scale);
00057     const int (&modified)(tup.modified);
00058 
00059     const CompAtom &p_i = p[0]->x[localIndex[0]];
00060     const CompAtom &p_j = p[1]->x[localIndex[1]];
00061 
00062     // compute vectors between atoms and their distances
00063     const Vector r12 = lattice.delta(p_i.position, p_j.position);
00064     BigReal r2 = r12.length2();
00065 
00066     if ( r2 > cutoff2 ) continue;
00067 
00068     if ( modified && r2 < 1.0 ) r2 = 1.0;  // match CUDA interpolation
00069 
00070     r2 += r2_delta;
00071 
00072     union { double f; int64 i; } r2i;
00073     r2i.f = r2;
00074     const int r2_delta_expc = 64 * (r2_delta_exp - 1023);
00075     int table_i = (r2i.i >> (32+14)) + r2_delta_expc;  // table_i >= 0
00076 
00077     const BigReal* const table_four_i = table_noshort + 16*table_i;
00078 
00079     BigReal diffa = r2 - r2_table[table_i];
00080 
00081     BigReal fast_a = 0., fast_b = 0., fast_c = 0., fast_d = 0.;
00082     BigReal slow_a, slow_b, slow_c, slow_d;
00083 
00084   if ( modified ) {
00085 
00086     const LJTable::TableEntry * lj_pars =
00087             ljTable->table_row(p_i.vdwType) + 2 * p_j.vdwType;
00088 
00089     // modified - normal = correction
00090     const BigReal A = scaling * ( (lj_pars+1)->A - lj_pars->A );
00091     const BigReal B = scaling * ( (lj_pars+1)->B - lj_pars->B );
00092 
00093     BigReal vdw_d = A * table_four_i[0] - B * table_four_i[4];
00094     BigReal vdw_c = A * table_four_i[1] - B * table_four_i[5];
00095     BigReal vdw_b = A * table_four_i[2] - B * table_four_i[6];
00096     BigReal vdw_a = A * table_four_i[3] - B * table_four_i[7];
00097 
00098     const BigReal kqq = (1.0 - scale14) *
00099             COULOMB * p_i.charge * p_j.charge * scaling * dielectric_1;
00100 
00101     fast_a =      kqq * fast_table[4*table_i+0];  // not used!
00102     fast_b = 2. * kqq * fast_table[4*table_i+1];
00103     fast_c = 4. * kqq * fast_table[4*table_i+2];
00104     fast_d = 6. * kqq * fast_table[4*table_i+3];
00105 
00106     if ( doFull ) {
00107       slow_a =      kqq * slow_table[4*table_i+3];  // not used!
00108       slow_b = 2. * kqq * slow_table[4*table_i+2];
00109       slow_c = 4. * kqq * slow_table[4*table_i+1];
00110       slow_d = 6. * kqq * slow_table[4*table_i+0];
00111     }
00112 
00113     if ( doEnergy ) {
00114       reduction[vdwEnergyIndex] -=
00115                 ( ( diffa * (1./6.)*vdw_d + 0.25*vdw_c ) * diffa
00116                         + 0.5*vdw_b ) * diffa + vdw_a;
00117       reduction[electEnergyIndex] -=
00118                 ( ( diffa * (1./6.)*fast_d + 0.25*fast_c ) * diffa
00119                         + 0.5*fast_b ) * diffa + fast_a;
00120       if ( doFull ) {
00121         reduction[fullElectEnergyIndex] -=
00122                 ( ( diffa * (1./6.)*slow_d + 0.25*slow_c ) * diffa
00123                         + 0.5*slow_b ) * diffa + slow_a;
00124       }
00125     }
00126 
00127     fast_a += vdw_a;
00128     fast_b += vdw_b;
00129     fast_c += vdw_c;
00130     fast_d += vdw_d;
00131 
00132   } else if ( doFull ) {  // full exclusion
00133 
00134     const BigReal kqq = 
00135             COULOMB * p_i.charge * p_j.charge * scaling * dielectric_1;
00136 
00137     slow_d = kqq * ( table_four_i[8]  - table_four_i[12] );
00138     slow_c = kqq * ( table_four_i[9]  - table_four_i[13] );
00139     slow_b = kqq * ( table_four_i[10] - table_four_i[14] );
00140     slow_a = kqq * ( table_four_i[11] - table_four_i[15] );  // not used!
00141 
00142     if ( doEnergy ) {
00143       reduction[fullElectEnergyIndex] -=
00144                 ( ( diffa * (1./6.)*slow_d + 0.25*slow_c ) * diffa
00145                         + 0.5*slow_b ) * diffa + slow_a;
00146     }
00147   }
00148 
00149   register BigReal fast_dir =
00150                   (diffa * fast_d + fast_c) * diffa + fast_b;
00151 
00152   const Force f12 = fast_dir * r12;
00153 
00154   //  Now add the forces to each force vector
00155   p[0]->r->f[Results::nbond][localIndex[0]] += f12;
00156   p[1]->r->f[Results::nbond][localIndex[1]] -= f12;
00157 
00158   // reduction[nonbondedEnergyIndex] += energy;
00159   reduction[virialIndex_XX] += f12.x * r12.x;
00160   reduction[virialIndex_XY] += f12.x * r12.y;
00161   reduction[virialIndex_XZ] += f12.x * r12.z;
00162   reduction[virialIndex_YX] += f12.y * r12.x;
00163   reduction[virialIndex_YY] += f12.y * r12.y;
00164   reduction[virialIndex_YZ] += f12.y * r12.z;
00165   reduction[virialIndex_ZX] += f12.z * r12.x;
00166   reduction[virialIndex_ZY] += f12.z * r12.y;
00167   reduction[virialIndex_ZZ] += f12.z * r12.z;
00168 
00169   if ( doFull ) {
00170     register BigReal slow_dir =
00171                   (diffa * slow_d + slow_c) * diffa + slow_b;
00172 
00173     const Force slow_f12 = slow_dir * r12;
00174 
00175     p[0]->r->f[Results::slow][localIndex[0]] += slow_f12;
00176     p[1]->r->f[Results::slow][localIndex[1]] -= slow_f12;
00177 
00178     // reduction[nonbondedEnergyIndex] += energy;
00179     reduction[slowVirialIndex_XX] += slow_f12.x * r12.x;
00180     reduction[slowVirialIndex_XY] += slow_f12.x * r12.y;
00181     reduction[slowVirialIndex_XZ] += slow_f12.x * r12.z;
00182     reduction[slowVirialIndex_YX] += slow_f12.y * r12.x;
00183     reduction[slowVirialIndex_YY] += slow_f12.y * r12.y;
00184     reduction[slowVirialIndex_YZ] += slow_f12.y * r12.z;
00185     reduction[slowVirialIndex_ZX] += slow_f12.z * r12.x;
00186     reduction[slowVirialIndex_ZY] += slow_f12.z * r12.y;
00187     reduction[slowVirialIndex_ZZ] += slow_f12.z * r12.z;
00188   }
00189 
00190   }
00191 }
00192 
00193 void ExclElem::submitReductionData(BigReal *data, SubmitReduction *reduction)
00194 {
00195   reduction->item(REDUCTION_ELECT_ENERGY) += data[electEnergyIndex];
00196   reduction->item(REDUCTION_LJ_ENERGY) += data[vdwEnergyIndex];
00197   reduction->item(REDUCTION_ELECT_ENERGY_SLOW) += data[fullElectEnergyIndex];
00198   ADD_TENSOR(reduction,REDUCTION_VIRIAL_NBOND,data,virialIndex);
00199   ADD_TENSOR(reduction,REDUCTION_VIRIAL_SLOW,data,slowVirialIndex);
00200 }
00201 

Generated on Thu Sep 21 01:17:11 2017 for NAMD by  doxygen 1.4.7