ComputeThole.C

Go to the documentation of this file.
00001 
00007 #include "InfoStream.h"
00008 #include "ComputeThole.h"
00009 #include "Molecule.h"
00010 #include "Parameters.h"
00011 #include "Node.h"
00012 #include "ReductionMgr.h"
00013 #include "Lattice.h"
00014 #include "PressureProfile.h"
00015 #include "Debug.h"
00016 
00017 #define CALCULATE_THOLE_CORRECTION
00018 
00019 // static initialization
00020 int TholeElem::pressureProfileSlabs = 0;
00021 int TholeElem::pressureProfileAtomTypes = 1;
00022 BigReal TholeElem::pressureProfileThickness = 0;
00023 BigReal TholeElem::pressureProfileMin = 0;
00024 
00025 void TholeElem::getMoleculePointers
00026     (Molecule* mol, int* count, int32*** byatom, Thole** structarray)
00027 {
00028 #ifdef MEM_OPT_VERSION
00029   NAMD_die("Should not be called in TholeElem::getMoleculePointers in memory optimized version!");
00030 #else
00031   *count = mol->numTholes;
00032   *byatom = mol->tholesByAtom;
00033   *structarray = mol->tholes;
00034 #endif
00035 }
00036 
00037 void TholeElem::getParameterPointers(Parameters *p, const TholeValue **v) {
00038   *v = NULL;  // parameters are stored in the structure
00039 }
00040 
00041 void TholeElem::computeForce(TholeElem *tuples, int ntuple, BigReal *reduction, 
00042                                 BigReal *pressureProfileData)
00043 {
00044  const Lattice & lattice = tuples[0].p[0]->p->lattice;
00045 
00046  for ( int ituple=0; ituple<ntuple; ++ituple ) {
00047   const TholeElem &tup = tuples[ituple];
00048   enum { size = 4 };
00049   const AtomID (&atomID)[size](tup.atomID);
00050   const int    (&localIndex)[size](tup.localIndex);
00051   TuplePatchElem * const(&p)[size](tup.p);
00052   const Real (&scale)(tup.scale);
00053   const TholeValue * const(&value)(tup.value);
00054 
00055   DebugM(3, "::computeForce() localIndex = " << localIndex[0] << " "
00056                << localIndex[1] << " " << localIndex[2] << " "
00057                << localIndex[3] << std::endl);
00058 
00059 #ifdef CALCULATE_THOLE_CORRECTION
00060   const BigReal aa = value->aa;
00061   const BigReal qq = value->qq;
00062 
00063   //  Calculate the vectors between atoms
00064   const Position & rai = p[0]->x[localIndex[0]].position;  // atom i
00065   const Position & rdi = p[1]->x[localIndex[1]].position;  // atom i's Drude
00066   const Position & raj = p[2]->x[localIndex[2]].position;  // atom j
00067   const Position & rdj = p[3]->x[localIndex[3]].position;  // atom j's Drude
00068 
00069   // r_ij = r_i - r_j
00070   Vector raa = lattice.delta(rai,raj);  // shortest vector image:  rai - raj
00071   Vector rad = lattice.delta(rai,rdj);  // shortest vector image:  rai - rdj
00072   Vector rda = lattice.delta(rdi,raj);  // shortest vector image:  rdi - raj
00073   Vector rdd = lattice.delta(rdi,rdj);  // shortest vector image:  rdi - rdj
00074 
00075   // 1/r, r = |r_ij|
00076   BigReal raa_invlen = raa.rlength();  // reciprocal of length
00077   BigReal rad_invlen = rad.rlength();
00078   BigReal rda_invlen = rda.rlength();
00079   BigReal rdd_invlen = rdd.rlength();
00080 
00081   // ar
00082   BigReal auaa = aa / raa_invlen;
00083   BigReal auad = aa / rad_invlen;
00084   BigReal auda = aa / rda_invlen;
00085   BigReal audd = aa / rdd_invlen;
00086 
00087   // exp(-ar)
00088   BigReal expauaa = exp(-auaa);
00089   BigReal expauad = exp(-auad);
00090   BigReal expauda = exp(-auda);
00091   BigReal expaudd = exp(-audd);
00092 
00093   // (1 + ar/2)
00094   BigReal polyauaa = 1 + 0.5*auaa;
00095   BigReal polyauad = 1 + 0.5*auad;
00096   BigReal polyauda = 1 + 0.5*auda;
00097   BigReal polyaudd = 1 + 0.5*audd;
00098 
00099   // U(r) = qq/r (1 - (1 + ar/2) exp(-ar))
00100   BigReal ethole = 0;
00101   ethole += qq * raa_invlen * (1 - polyauaa * expauaa);
00102   ethole += -qq * rad_invlen * (1 - polyauad * expauad);
00103   ethole += -qq * rda_invlen * (1 - polyauda * expauda);
00104   ethole += qq * rdd_invlen * (1 - polyaudd * expaudd);
00105 
00106   polyauaa = 1 + auaa*polyauaa;
00107   polyauad = 1 + auad*polyauad;
00108   polyauda = 1 + auda*polyauda;
00109   polyaudd = 1 + audd*polyaudd;
00110 
00111   BigReal raa_invlen3 = raa_invlen * raa_invlen * raa_invlen;
00112   BigReal rad_invlen3 = rad_invlen * rad_invlen * rad_invlen;
00113   BigReal rda_invlen3 = rda_invlen * rda_invlen * rda_invlen;
00114   BigReal rdd_invlen3 = rdd_invlen * rdd_invlen * rdd_invlen;
00115 
00116   // df = (1/r) (dU/dr)
00117   BigReal dfaa = qq * raa_invlen3 * (polyauaa*expauaa - 1);
00118   BigReal dfad = -qq * rad_invlen3 * (polyauad*expauad - 1);
00119   BigReal dfda = -qq * rda_invlen3 * (polyauda*expauda - 1);
00120   BigReal dfdd = qq * rdd_invlen3 * (polyaudd*expaudd - 1);
00121 
00122   Vector faa = -dfaa * raa;
00123   Vector fad = -dfad * rad;
00124   Vector fda = -dfda * rda;
00125   Vector fdd = -dfdd * rdd;
00126 
00127   p[0]->f[localIndex[0]] += faa + fad;
00128   p[1]->f[localIndex[1]] += fda + fdd;
00129   p[2]->f[localIndex[2]] -= faa + fda;
00130   p[3]->f[localIndex[3]] -= fad + fdd;
00131 
00132   DebugM(3, "::computeForce() -- ending with delta energy " << ethole
00133       << std::endl);
00134   reduction[tholeEnergyIndex] += ethole;
00135 
00136   reduction[virialIndex_XX] += faa.x * raa.x + fad.x * rad.x
00137     + fda.x * rda.x + fdd.x * rdd.x;
00138   reduction[virialIndex_XY] += faa.x * raa.y + fad.x * rad.y
00139     + fda.x * rda.y + fdd.x * rdd.y;
00140   reduction[virialIndex_XZ] += faa.x * raa.z + fad.x * rad.z
00141     + fda.x * rda.z + fdd.x * rdd.z;
00142   reduction[virialIndex_YX] += faa.y * raa.x + fad.y * rad.x
00143     + fda.y * rda.x + fdd.y * rdd.x;
00144   reduction[virialIndex_YY] += faa.y * raa.y + fad.y * rad.y
00145     + fda.y * rda.y + fdd.y * rdd.y;
00146   reduction[virialIndex_YZ] += faa.y * raa.z + fad.y * rad.z
00147     + fda.y * rda.z + fdd.y * rdd.z;
00148   reduction[virialIndex_ZX] += faa.z * raa.x + fad.z * rad.x
00149     + fda.z * rda.x + fdd.z * rdd.x;
00150   reduction[virialIndex_ZY] += faa.z * raa.y + fad.z * rad.y
00151     + fda.z * rda.y + fdd.z * rdd.y;
00152   reduction[virialIndex_ZZ] += faa.z * raa.z + fad.z * rad.z
00153     + fda.z * rda.z + fdd.z * rdd.z;
00154 
00155   if (pressureProfileData) {
00156     BigReal zai = p[0]->x[localIndex[0]].position.z;
00157     BigReal zdi = p[1]->x[localIndex[1]].position.z;
00158     BigReal zaj = p[2]->x[localIndex[2]].position.z;
00159     BigReal zdj = p[3]->x[localIndex[3]].position.z;
00160     int nai = (int)floor((zai-pressureProfileMin)/pressureProfileThickness);
00161     int ndi = (int)floor((zdi-pressureProfileMin)/pressureProfileThickness);
00162     int naj = (int)floor((zaj-pressureProfileMin)/pressureProfileThickness);
00163     int ndj = (int)floor((zdj-pressureProfileMin)/pressureProfileThickness);
00164     pp_clamp(nai, pressureProfileSlabs);
00165     pp_clamp(ndi, pressureProfileSlabs);
00166     pp_clamp(naj, pressureProfileSlabs);
00167     pp_clamp(ndj, pressureProfileSlabs);
00168     int pai = p[0]->x[localIndex[0]].partition;
00169     int pdi = p[1]->x[localIndex[1]].partition;
00170     int paj = p[2]->x[localIndex[2]].partition;
00171     int pdj = p[3]->x[localIndex[3]].partition;
00172     int pn = pressureProfileAtomTypes;
00173     pp_reduction(pressureProfileSlabs, nai, naj,
00174         pai, paj, pn, faa.x * raa.x, faa.y * raa.y, faa.z * raa.z,
00175         pressureProfileData);
00176     pp_reduction(pressureProfileSlabs, nai, ndj,
00177         pai, pdj, pn, fad.x * rad.x, fad.y * rad.y, fad.z * rad.z,
00178         pressureProfileData);
00179     pp_reduction(pressureProfileSlabs, ndi, naj,
00180         pdi, paj, pn, fda.x * rda.x, fda.y * rda.y, fda.z * rda.z,
00181         pressureProfileData);
00182     pp_reduction(pressureProfileSlabs, ndi, ndj,
00183         pdi, pdj, pn, fdd.x * rdd.x, fdd.y * rdd.y, fdd.z * rdd.z,
00184         pressureProfileData);
00185   }
00186 #endif
00187 
00188  }
00189 }
00190 
00191 
00192 // The energy from the screened Coulomb correction of Thole is 
00193 // accumulated into the electrostatic potential energy.
00194 void TholeElem::submitReductionData(BigReal *data, SubmitReduction *reduction)
00195 {
00196   reduction->item(REDUCTION_ELECT_ENERGY) += data[tholeEnergyIndex];
00197   ADD_TENSOR(reduction,REDUCTION_VIRIAL_NORMAL,data,virialIndex);
00198 }
00199 

Generated on Thu Nov 23 01:17:12 2017 for NAMD by  doxygen 1.4.7