ComputeAniso.C

Go to the documentation of this file.
00001 
00007 #include "InfoStream.h"
00008 #include "ComputeAniso.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_ANISO
00018 
00019 // static initialization
00020 int AnisoElem::pressureProfileSlabs = 0;
00021 int AnisoElem::pressureProfileAtomTypes = 1;
00022 BigReal AnisoElem::pressureProfileThickness = 0;
00023 BigReal AnisoElem::pressureProfileMin = 0;
00024 
00025 void AnisoElem::getMoleculePointers
00026     (Molecule* mol, int* count, int32*** byatom, Aniso** structarray)
00027 {
00028 #ifdef MEM_OPT_VERSION
00029   NAMD_die("Should not be called in AnisoElem::getMoleculePointers in memory optimized version!");
00030 #else
00031   *count = mol->numAnisos;
00032   *byatom = mol->anisosByAtom;
00033   *structarray = mol->anisos;
00034 #endif
00035 }
00036 
00037 void AnisoElem::getParameterPointers(Parameters *p, const AnisoValue **v) {
00038   *v = NULL;  // parameters are stored in the structure
00039 }
00040 
00041 void AnisoElem::computeForce(AnisoElem *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 AnisoElem &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 AnisoValue * 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_ANISO
00060   // used some comments from Ed Harder's implementation in CHARMM
00061 
00062   const BigReal kpar0  = 2*value->k11;  // force constants
00063   const BigReal kperp0 = 2*value->k22;
00064   const BigReal kiso0  = 2*value->k33;
00065 
00066   const Position & ri = p[0]->x[localIndex[0]].position;    // atom I
00067   const Position & rj = p[0]->x[localIndex[0]+1].position;  // atom I's Drude
00068   const Position & rl = p[1]->x[localIndex[1]].position;    // atom L
00069   const Position & rm = p[2]->x[localIndex[2]].position;    // atom M
00070   const Position & rn = p[3]->x[localIndex[3]].position;    // atom N
00071 
00072   // calculate parallel and perpendicular displacement vectors
00073   Vector r_il = lattice.delta(ri,rl);  // shortest vector image:  ri - rl
00074   Vector r_mn = lattice.delta(rm,rn);  // shortest vector image:  rm - rn
00075 
00076   BigReal r_il_invlen = r_il.rlength();  // need recip lengths of r_il, r_mn
00077   BigReal r_mn_invlen = r_mn.rlength();
00078 
00079   Vector u1 = r_il * r_il_invlen;  // normalize r_il, r_mn
00080   Vector u2 = r_mn * r_mn_invlen;
00081 
00082   Vector dr = rj - ri;  // Drude displacement vector (ri, rj are in same patch)
00083 
00084   BigReal dpar  = dr * u1;  // parallel displacement
00085   BigReal dperp = dr * u2;  // perpendicular displacement
00086 
00087   // aniso spring energy
00088   // kpar reduces response along carbonyl vector
00089   // kperp reduces response perp to bond vector
00090   //   (reg in and out of plane response)
00091   BigReal eaniso;
00092   eaniso = 0.5*kpar0*dpar*dpar + 0.5*kperp0*dperp*dperp + 0.5*kiso0*(dr*dr);
00093 
00094   // calculate force vectors in one direction only:
00095   // fi = -(fj + fl),  fn = -fm
00096 
00097   // force on atom j
00098   Vector fj = -kiso0 * dr;
00099   fj -= kpar0 * dpar * u1;
00100   fj -= kperp0 * dperp * u2;
00101 
00102   // force on atom l
00103   Vector fl = kpar0 * dpar * r_il_invlen * dr;
00104   fl -= kpar0 * dpar * dpar * r_il_invlen * u1;
00105 
00106   // force on atom m
00107   Vector fm = kperp0 * dperp * dperp * r_mn_invlen * u2;
00108   fm -= kperp0 * dperp * r_mn_invlen * dr;
00109 
00110   // accumulate forces
00111   p[0]->f[localIndex[0]] -= (fj + fl);
00112   p[0]->f[localIndex[0]+1] += fj;
00113   p[1]->f[localIndex[1]] += fl;
00114   p[2]->f[localIndex[2]] += fm;
00115   p[3]->f[localIndex[3]] -= fm;
00116 
00117   // update potential
00118   reduction[anisoEnergyIndex] += eaniso;
00119 
00120   // update virial
00121   reduction[virialIndex_XX] += fj.x * dr.x - fl.x * r_il.x + fm.x * r_mn.x;
00122   reduction[virialIndex_XY] += fj.x * dr.y - fl.x * r_il.y + fm.x * r_mn.y;
00123   reduction[virialIndex_XZ] += fj.x * dr.z - fl.x * r_il.z + fm.x * r_mn.z;
00124   reduction[virialIndex_YX] += fj.y * dr.x - fl.y * r_il.x + fm.y * r_mn.x;
00125   reduction[virialIndex_YY] += fj.y * dr.y - fl.y * r_il.y + fm.y * r_mn.y;
00126   reduction[virialIndex_YZ] += fj.y * dr.z - fl.y * r_il.z + fm.y * r_mn.z;
00127   reduction[virialIndex_ZX] += fj.z * dr.x - fl.z * r_il.x + fm.z * r_mn.x;
00128   reduction[virialIndex_ZY] += fj.z * dr.y - fl.z * r_il.y + fm.z * r_mn.y;
00129   reduction[virialIndex_ZZ] += fj.z * dr.z - fl.z * r_il.z + fm.z * r_mn.z;
00130 
00131   // update pressure profile data
00132   if (pressureProfileData) {
00133     BigReal zi = p[0]->x[localIndex[0]].position.z;
00134     BigReal zj = p[0]->x[localIndex[0]+1].position.z;
00135     BigReal zl = p[1]->x[localIndex[1]].position.z;
00136     BigReal zm = p[2]->x[localIndex[2]].position.z;
00137     BigReal zn = p[3]->x[localIndex[3]].position.z;
00138     int ni = (int)floor((zi-pressureProfileMin)/pressureProfileThickness);
00139     int nj = (int)floor((zj-pressureProfileMin)/pressureProfileThickness);
00140     int nl = (int)floor((zl-pressureProfileMin)/pressureProfileThickness);
00141     int nm = (int)floor((zm-pressureProfileMin)/pressureProfileThickness);
00142     int nn = (int)floor((zn-pressureProfileMin)/pressureProfileThickness);
00143     pp_clamp(ni, pressureProfileSlabs);
00144     pp_clamp(nj, pressureProfileSlabs);
00145     pp_clamp(nl, pressureProfileSlabs);
00146     pp_clamp(nm, pressureProfileSlabs);
00147     pp_clamp(nn, pressureProfileSlabs);
00148     int pi = p[0]->x[localIndex[0]].partition;
00149     int pj = p[0]->x[localIndex[0]+1].partition;
00150     int pl = p[1]->x[localIndex[1]].partition;
00151     int pm = p[2]->x[localIndex[2]].partition;
00152     int pn = p[3]->x[localIndex[3]].partition;
00153     int pt = pressureProfileAtomTypes;
00154     pp_reduction(pressureProfileSlabs, nj, ni,
00155         pj, pi, pt, fj.x * dr.x, fj.y * dr.y, fj.z * dr.z,
00156         pressureProfileData);
00157     pp_reduction(pressureProfileSlabs, ni, nl,
00158         pi, pl, pt, -fl.x * r_il.x, -fl.y * r_il.y, -fl.z * r_il.z,
00159         pressureProfileData);
00160     pp_reduction(pressureProfileSlabs, nm, nn,
00161         pm, pn, pt, fm.x * r_mn.x, fm.y * r_mn.y, fm.z * r_mn.z,
00162         pressureProfileData);
00163   }
00164 #endif
00165 
00166  }
00167 }
00168 
00169 
00170 void AnisoElem::submitReductionData(BigReal *data, SubmitReduction *reduction)
00171 {
00172   reduction->item(REDUCTION_BOND_ENERGY) += data[anisoEnergyIndex];
00173   ADD_TENSOR(reduction,REDUCTION_VIRIAL_NORMAL,data,virialIndex);
00174 }
00175 

Generated on Sat Nov 18 01:17:12 2017 for NAMD by  doxygen 1.4.7