Main Page | Namespace List | Class Hierarchy | Alphabetical List | Class List | File List | Namespace Members | Class Members | File Members

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(BigReal *reduction, 
00042                             BigReal *pressureProfileData)
00043 {
00044     const Flags &flags = p[0]->p->flags;
00045     if ( ! flags.doNonbonded ) return;
00046 
00047     const CompAtom &p_i = p[0]->x[localIndex[0]];
00048     const CompAtom &p_j = p[1]->x[localIndex[1]];
00049 
00050     // compute vectors between atoms and their distances
00051     const Lattice & lattice = p[0]->p->lattice;
00052     const Vector r12 = lattice.delta(p_i.position, p_j.position);
00053     BigReal r2 = r12.length2();
00054 
00055     if ( r2 > cutoff2 ) return;
00056 
00057     r2 += r2_delta;
00058 
00059     union { double f; int64 i; } r2i;
00060     r2i.f = r2;
00061     const int r2_delta_expc = 64 * (r2_delta_exp - 1023);
00062     int table_i = (r2i.i >> (32+14)) + r2_delta_expc;  // table_i >= 0
00063 
00064     const BigReal* const table_four_i = table_noshort + 16*table_i;
00065 
00066     BigReal diffa = r2 - r2_table[table_i];
00067 
00068     const int doFull = flags.doFullElectrostatics;
00069     const int doEnergy = flags.doEnergy;
00070 
00071     BigReal fast_a = 0., fast_b = 0., fast_c = 0., fast_d = 0.;
00072     BigReal slow_a, slow_b, slow_c, slow_d;
00073 
00074   if ( modified ) {
00075 
00076     const LJTable::TableEntry * lj_pars =
00077             ljTable->table_row(p_i.vdwType) + 2 * p_j.vdwType;
00078 
00079     // modified - normal = correction
00080     const BigReal A = scaling * ( (lj_pars+1)->A - lj_pars->A );
00081     const BigReal B = scaling * ( (lj_pars+1)->B - lj_pars->B );
00082 
00083     BigReal vdw_d = A * table_four_i[0] - B * table_four_i[4];
00084     BigReal vdw_c = A * table_four_i[1] - B * table_four_i[5];
00085     BigReal vdw_b = A * table_four_i[2] - B * table_four_i[6];
00086     BigReal vdw_a = A * table_four_i[3] - B * table_four_i[7];
00087 
00088     const BigReal kqq = (1.0 - scale14) *
00089             COULOMB * p_i.charge * p_j.charge * scaling * dielectric_1;
00090 
00091     fast_a =      kqq * fast_table[4*table_i+0];  // not used!
00092     fast_b = 2. * kqq * fast_table[4*table_i+1];
00093     fast_c = 4. * kqq * fast_table[4*table_i+2];
00094     fast_d = 6. * kqq * fast_table[4*table_i+3];
00095 
00096     if ( doFull ) {
00097       slow_a =      kqq * slow_table[4*table_i+3];  // not used!
00098       slow_b = 2. * kqq * slow_table[4*table_i+2];
00099       slow_c = 4. * kqq * slow_table[4*table_i+1];
00100       slow_d = 6. * kqq * slow_table[4*table_i+0];
00101     }
00102 
00103     if ( doEnergy ) {
00104       reduction[vdwEnergyIndex] -=
00105                 ( ( diffa * (1./6.)*vdw_d + 0.25*vdw_c ) * diffa
00106                         + 0.5*vdw_b ) * diffa + vdw_a;
00107       reduction[electEnergyIndex] -=
00108                 ( ( diffa * (1./6.)*fast_d + 0.25*fast_c ) * diffa
00109                         + 0.5*fast_b ) * diffa + fast_a;
00110       if ( doFull ) {
00111         reduction[fullElectEnergyIndex] -=
00112                 ( ( diffa * (1./6.)*slow_d + 0.25*slow_c ) * diffa
00113                         + 0.5*slow_b ) * diffa + slow_a;
00114       }
00115     }
00116 
00117     fast_a += vdw_a;
00118     fast_b += vdw_b;
00119     fast_c += vdw_c;
00120     fast_d += vdw_d;
00121 
00122   } else if ( doFull ) {  // full exclusion
00123 
00124     const BigReal kqq = 
00125             COULOMB * p_i.charge * p_j.charge * scaling * dielectric_1;
00126 
00127     slow_d = kqq * ( table_four_i[8]  - table_four_i[12] );
00128     slow_c = kqq * ( table_four_i[9]  - table_four_i[13] );
00129     slow_b = kqq * ( table_four_i[10] - table_four_i[14] );
00130     slow_a = kqq * ( table_four_i[11] - table_four_i[15] );  // not used!
00131 
00132     if ( doEnergy ) {
00133       reduction[fullElectEnergyIndex] -=
00134                 ( ( diffa * (1./6.)*slow_d + 0.25*slow_c ) * diffa
00135                         + 0.5*slow_b ) * diffa + slow_a;
00136     }
00137   }
00138 
00139   register BigReal fast_dir =
00140                   (diffa * fast_d + fast_c) * diffa + fast_b;
00141 
00142   const Force f12 = fast_dir * r12;
00143 
00144   //  Now add the forces to each force vector
00145   p[0]->r->f[Results::nbond][localIndex[0]] += f12;
00146   p[1]->r->f[Results::nbond][localIndex[1]] -= f12;
00147 
00148   // reduction[nonbondedEnergyIndex] += energy;
00149   reduction[virialIndex_XX] += f12.x * r12.x;
00150   reduction[virialIndex_XY] += f12.x * r12.y;
00151   reduction[virialIndex_XZ] += f12.x * r12.z;
00152   reduction[virialIndex_YX] += f12.y * r12.x;
00153   reduction[virialIndex_YY] += f12.y * r12.y;
00154   reduction[virialIndex_YZ] += f12.y * r12.z;
00155   reduction[virialIndex_ZX] += f12.z * r12.x;
00156   reduction[virialIndex_ZY] += f12.z * r12.y;
00157   reduction[virialIndex_ZZ] += f12.z * r12.z;
00158 
00159   if ( doFull ) {
00160     register BigReal slow_dir =
00161                   (diffa * slow_d + slow_c) * diffa + slow_b;
00162 
00163     const Force slow_f12 = slow_dir * r12;
00164 
00165     p[0]->r->f[Results::slow][localIndex[0]] += slow_f12;
00166     p[1]->r->f[Results::slow][localIndex[1]] -= slow_f12;
00167 
00168     // reduction[nonbondedEnergyIndex] += energy;
00169     reduction[slowVirialIndex_XX] += slow_f12.x * r12.x;
00170     reduction[slowVirialIndex_XY] += slow_f12.x * r12.y;
00171     reduction[slowVirialIndex_XZ] += slow_f12.x * r12.z;
00172     reduction[slowVirialIndex_YX] += slow_f12.y * r12.x;
00173     reduction[slowVirialIndex_YY] += slow_f12.y * r12.y;
00174     reduction[slowVirialIndex_YZ] += slow_f12.y * r12.z;
00175     reduction[slowVirialIndex_ZX] += slow_f12.z * r12.x;
00176     reduction[slowVirialIndex_ZY] += slow_f12.z * r12.y;
00177     reduction[slowVirialIndex_ZZ] += slow_f12.z * r12.z;
00178   }
00179 
00180 }
00181 
00182 void ExclElem::submitReductionData(BigReal *data, SubmitReduction *reduction)
00183 {
00184   reduction->item(REDUCTION_ELECT_ENERGY) += data[electEnergyIndex];
00185   reduction->item(REDUCTION_LJ_ENERGY) += data[vdwEnergyIndex];
00186   reduction->item(REDUCTION_ELECT_ENERGY_SLOW) += data[fullElectEnergyIndex];
00187   ADD_TENSOR(reduction,REDUCTION_VIRIAL_NBOND,data,virialIndex);
00188   ADD_TENSOR(reduction,REDUCTION_VIRIAL_SLOW,data,slowVirialIndex);
00189 }
00190 

Generated on Wed May 22 04:07:15 2013 for NAMD by  doxygen 1.3.9.1