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  //fepb BKR
00047  SimParameters *const simParams = Node::Object()->simParameters;
00048  const int step = tuples[0].p[0]->p->flags.step;
00049  const BigReal alchLambda = simParams->getCurrentLambda(step);
00050  const BigReal alchLambda2 = simParams->alchLambda2;
00051  const BigReal bond_lambda_1 = simParams->getBondLambda(alchLambda);
00052  const BigReal bond_lambda_2 = simParams->getBondLambda(1-alchLambda);
00053  const BigReal bond_lambda_12 = simParams->getBondLambda(alchLambda2);
00054  const BigReal bond_lambda_22 = simParams->getBondLambda(1-alchLambda2);
00055  Molecule *const mol = Node::Object()->molecule;
00056  //fepe
00057 
00058  for ( int ituple=0; ituple<ntuple; ++ituple ) {
00059   const AnisoElem &tup = tuples[ituple];
00060   enum { size = 4 };
00061   const AtomID (&atomID)[size](tup.atomID);
00062   const int    (&localIndex)[size](tup.localIndex);
00063   TuplePatchElem * const(&p)[size](tup.p);
00064   const Real (&scale)(tup.scale);
00065   const AnisoValue * const(&value)(tup.value);
00066 
00067   DebugM(3, "::computeForce() localIndex = " << localIndex[0] << " "
00068                << localIndex[1] << " " << localIndex[2] << " "
00069                << localIndex[3] << std::endl);
00070 
00071 #ifdef CALCULATE_ANISO
00072   // used some comments from Ed Harder's implementation in CHARMM
00073 
00074   const BigReal kpar0  = 2*value->k11;  // force constants
00075   const BigReal kperp0 = 2*value->k22;
00076   const BigReal kiso0  = 2*value->k33;
00077 
00078   const Position & ri = p[0]->x[localIndex[0]].position;    // atom I
00079   const Position & rj = p[0]->x[localIndex[0]+1].position;  // atom I's Drude
00080   const Position & rl = p[1]->x[localIndex[1]].position;    // atom L
00081   const Position & rm = p[2]->x[localIndex[2]].position;    // atom M
00082   const Position & rn = p[3]->x[localIndex[3]].position;    // atom N
00083 
00084   // calculate parallel and perpendicular displacement vectors
00085   Vector r_il = lattice.delta(ri,rl);  // shortest vector image:  ri - rl
00086   Vector r_mn = lattice.delta(rm,rn);  // shortest vector image:  rm - rn
00087 
00088   BigReal r_il_invlen = r_il.rlength();  // need recip lengths of r_il, r_mn
00089   BigReal r_mn_invlen = r_mn.rlength();
00090 
00091   Vector u1 = r_il * r_il_invlen;  // normalize r_il, r_mn
00092   Vector u2 = r_mn * r_mn_invlen;
00093 
00094   Vector dr = rj - ri;  // Drude displacement vector (ri, rj are in same patch)
00095 
00096   BigReal dpar  = dr * u1;  // parallel displacement
00097   BigReal dperp = dr * u2;  // perpendicular displacement
00098 
00099   // aniso spring energy
00100   // kpar reduces response along carbonyl vector
00101   // kperp reduces response perp to bond vector
00102   //   (reg in and out of plane response)
00103   BigReal eaniso;
00104   eaniso = 0.5*kpar0*dpar*dpar + 0.5*kperp0*dperp*dperp + 0.5*kiso0*(dr*dr);
00105 
00106   // calculate force vectors in one direction only:
00107   // fi = -(fj + fl),  fn = -fm
00108 
00109   // force on atom j
00110   Vector fj = -kiso0 * dr;
00111   fj -= kpar0 * dpar * u1;
00112   fj -= kperp0 * dperp * u2;
00113 
00114   // force on atom l
00115   Vector fl = kpar0 * dpar * r_il_invlen * dr;
00116   fl -= kpar0 * dpar * dpar * r_il_invlen * u1;
00117 
00118   // force on atom m
00119   Vector fm = kperp0 * dperp * dperp * r_mn_invlen * u2;
00120   fm -= kperp0 * dperp * r_mn_invlen * dr;
00121 
00122   //fepb - BKR scaling of alchemical bonded terms
00123   //       NB: TI derivative is the _unscaled_ energy.
00124   if ( simParams->alchOn ) {
00125     switch ( mol->get_fep_bonded_type(atomID, 4) ) {
00126     case 1:
00127       reduction[anisoEnergyIndex_ti_1] += eaniso;
00128       reduction[anisoEnergyIndex_f] += (bond_lambda_12 - bond_lambda_1)*eaniso;
00129       eaniso *= bond_lambda_1;
00130       fj *= bond_lambda_1;
00131       fl *= bond_lambda_1;
00132       fm *= bond_lambda_1;
00133       break;
00134     case 2:
00135       reduction[anisoEnergyIndex_ti_2] += eaniso;
00136       reduction[anisoEnergyIndex_f] += (bond_lambda_22 - bond_lambda_2)*eaniso;
00137       eaniso *= bond_lambda_2;
00138       fj *= bond_lambda_2;
00139       fl *= bond_lambda_2;
00140       fm *= bond_lambda_2; 
00141       break;
00142     }
00143   }
00144   //fepe
00145 
00146   // accumulate forces
00147   p[0]->f[localIndex[0]] -= (fj + fl);
00148   p[0]->f[localIndex[0]+1] += fj;
00149   p[1]->f[localIndex[1]] += fl;
00150   p[2]->f[localIndex[2]] += fm;
00151   p[3]->f[localIndex[3]] -= fm;
00152 
00153   // update potential
00154   reduction[anisoEnergyIndex] += eaniso;
00155 
00156   // update virial
00157   reduction[virialIndex_XX] += fj.x * dr.x - fl.x * r_il.x + fm.x * r_mn.x;
00158   reduction[virialIndex_XY] += fj.x * dr.y - fl.x * r_il.y + fm.x * r_mn.y;
00159   reduction[virialIndex_XZ] += fj.x * dr.z - fl.x * r_il.z + fm.x * r_mn.z;
00160   reduction[virialIndex_YX] += fj.y * dr.x - fl.y * r_il.x + fm.y * r_mn.x;
00161   reduction[virialIndex_YY] += fj.y * dr.y - fl.y * r_il.y + fm.y * r_mn.y;
00162   reduction[virialIndex_YZ] += fj.y * dr.z - fl.y * r_il.z + fm.y * r_mn.z;
00163   reduction[virialIndex_ZX] += fj.z * dr.x - fl.z * r_il.x + fm.z * r_mn.x;
00164   reduction[virialIndex_ZY] += fj.z * dr.y - fl.z * r_il.y + fm.z * r_mn.y;
00165   reduction[virialIndex_ZZ] += fj.z * dr.z - fl.z * r_il.z + fm.z * r_mn.z;
00166 
00167   // update pressure profile data
00168   if (pressureProfileData) {
00169     BigReal zi = p[0]->x[localIndex[0]].position.z;
00170     BigReal zj = p[0]->x[localIndex[0]+1].position.z;
00171     BigReal zl = p[1]->x[localIndex[1]].position.z;
00172     BigReal zm = p[2]->x[localIndex[2]].position.z;
00173     BigReal zn = p[3]->x[localIndex[3]].position.z;
00174     int ni = (int)floor((zi-pressureProfileMin)/pressureProfileThickness);
00175     int nj = (int)floor((zj-pressureProfileMin)/pressureProfileThickness);
00176     int nl = (int)floor((zl-pressureProfileMin)/pressureProfileThickness);
00177     int nm = (int)floor((zm-pressureProfileMin)/pressureProfileThickness);
00178     int nn = (int)floor((zn-pressureProfileMin)/pressureProfileThickness);
00179     pp_clamp(ni, pressureProfileSlabs);
00180     pp_clamp(nj, pressureProfileSlabs);
00181     pp_clamp(nl, pressureProfileSlabs);
00182     pp_clamp(nm, pressureProfileSlabs);
00183     pp_clamp(nn, pressureProfileSlabs);
00184     int pi = p[0]->x[localIndex[0]].partition;
00185     int pj = p[0]->x[localIndex[0]+1].partition;
00186     int pl = p[1]->x[localIndex[1]].partition;
00187     int pm = p[2]->x[localIndex[2]].partition;
00188     int pn = p[3]->x[localIndex[3]].partition;
00189     int pt = pressureProfileAtomTypes;
00190     pp_reduction(pressureProfileSlabs, nj, ni,
00191         pj, pi, pt, fj.x * dr.x, fj.y * dr.y, fj.z * dr.z,
00192         pressureProfileData);
00193     pp_reduction(pressureProfileSlabs, ni, nl,
00194         pi, pl, pt, -fl.x * r_il.x, -fl.y * r_il.y, -fl.z * r_il.z,
00195         pressureProfileData);
00196     pp_reduction(pressureProfileSlabs, nm, nn,
00197         pm, pn, pt, fm.x * r_mn.x, fm.y * r_mn.y, fm.z * r_mn.z,
00198         pressureProfileData);
00199   }
00200 #endif
00201 
00202  }
00203 }
00204 
00205 
00206 void AnisoElem::submitReductionData(BigReal *data, SubmitReduction *reduction)
00207 {
00208   reduction->item(REDUCTION_BOND_ENERGY) += data[anisoEnergyIndex];
00209   reduction->item(REDUCTION_BONDED_ENERGY_F) += data[anisoEnergyIndex_f];
00210   reduction->item(REDUCTION_BONDED_ENERGY_TI_1) += data[anisoEnergyIndex_ti_1];
00211   reduction->item(REDUCTION_BONDED_ENERGY_TI_2) += data[anisoEnergyIndex_ti_2];
00212   ADD_TENSOR(reduction,REDUCTION_VIRIAL_NORMAL,data,virialIndex);
00213 }
00214 

Generated on Tue May 22 01:17:12 2018 for NAMD by  doxygen 1.4.7