ComputeBonds.C

Go to the documentation of this file.
00001 
00007 #include "InfoStream.h"
00008 #include "ComputeBonds.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 
00018 // static initialization
00019 int BondElem::pressureProfileSlabs = 0;
00020 int BondElem::pressureProfileAtomTypes = 1;
00021 BigReal BondElem::pressureProfileThickness = 0;
00022 BigReal BondElem::pressureProfileMin = 0;
00023 
00024 void BondElem::getMoleculePointers
00025     (Molecule* mol, int* count, int32*** byatom, Bond** structarray)
00026 {
00027 #ifdef MEM_OPT_VERSION
00028   NAMD_die("Should not be called in BondElem::getMoleculePointers in memory optimized version!");
00029 #else
00030   *count = mol->numBonds;
00031   *byatom = mol->bondsByAtom;
00032   *structarray = mol->bonds;
00033 #endif
00034 }
00035 
00036 void BondElem::getParameterPointers(Parameters *p, const BondValue **v) {
00037   *v = p->bond_array;
00038 }
00039 
00040 void BondElem::computeForce(BondElem *tuples, int ntuple, BigReal *reduction, 
00041                             BigReal *pressureProfileData)
00042 {
00043  SimParameters *const simParams = Node::Object()->simParameters;
00044  const Lattice & lattice = tuples[0].p[0]->p->lattice;
00045 
00046  //fepb BKR
00047  const int step = tuples[0].p[0]->p->flags.step;
00048  const BigReal alchLambda = simParams->getCurrentLambda(step);
00049  const BigReal alchLambda2 = simParams->alchLambda2;
00050  const BigReal bond_lambda_1 = simParams->getBondLambda(alchLambda);
00051  const BigReal bond_lambda_2 = simParams->getBondLambda(1-alchLambda);
00052  const BigReal bond_lambda_12 = simParams->getBondLambda(alchLambda2);
00053  const BigReal bond_lambda_22 = simParams->getBondLambda(1-alchLambda2);
00054  Molecule *const mol = Node::Object()->molecule;
00055  //fepe
00056 
00057  for ( int ituple=0; ituple<ntuple; ++ituple ) {
00058   const BondElem &tup = tuples[ituple];
00059   enum { size = 2 };
00060   const AtomID (&atomID)[size](tup.atomID);
00061   const int    (&localIndex)[size](tup.localIndex);
00062   TuplePatchElem * const(&p)[size](tup.p);
00063   const Real (&scale)(tup.scale);
00064   const BondValue * const(&value)(tup.value);
00065 
00066   DebugM(1, "::computeForce() localIndex = " << localIndex[0] << " "
00067                << localIndex[1] << std::endl);
00068 
00069   // skip Lonepair bonds (other k=0. bonds have been filtered out)
00070   if (0. == value->k) continue;
00071 
00072   BigReal scal;         // force scaling
00073   BigReal energy;       // energy from the bond
00074 
00075   //DebugM(3, "::computeForce() -- starting with bond type " << bondType << std::endl);
00076 
00077   // get the bond information
00078   Real k = value->k * scale;
00079   Real x0 = value->x0;
00080 
00081   // compute vectors between atoms and their distances
00082   const Vector r12 = lattice.delta(p[0]->x[localIndex[0]].position,
00083                                         p[1]->x[localIndex[1]].position);
00084 
00085   if (0. == x0) {
00086     BigReal r2 = r12.length2();
00087     scal = -2.0*k;    // bond interaction for equilibrium length 0
00088     energy = k*r2;
00089     // TODO: Instead flag by mass for Drude particles, otherwise mixing and 
00090     // matching Drude and alchemical "twin" particles will likely not work as 
00091     // expected.
00092     if(simParams->drudeOn) { // for Drude bonds
00093       BigReal drudeBondLen = simParams->drudeBondLen;
00094       BigReal drudeBondConst = simParams->drudeBondConst;
00095       if ( (drudeBondConst > 0) && ( ! simParams->drudeHardWallOn ) &&
00096           (r2 > drudeBondLen*drudeBondLen) ) {
00097         // add a quartic restraining potential to keep Drude bond short
00098         BigReal r = sqrt(r2);
00099         BigReal diff = r - drudeBondLen;
00100         BigReal diff2 = diff*diff;
00101 
00102         scal += -4*drudeBondConst * diff2 * diff / r;
00103         energy += drudeBondConst * diff2 * diff2;
00104       }
00105     }
00106   }
00107   else {
00108     BigReal r = r12.length();  // Distance between atoms
00109     BigReal diff = r - x0;     // Compare it to the rest bond
00110 
00111     //  Add the energy from this bond to the total energy
00112     energy = k*diff*diff;
00113 
00114     //  Determine the magnitude of the force
00115     diff *= -2.0*k;
00116 
00117     //  Scale the force vector accordingly
00118     scal = (diff/r);
00119   }
00120 
00121   //fepb - BKR scaling of alchemical bonded terms
00122   //       NB: TI derivative is the _unscaled_ energy.
00123   if ( simParams->alchOn ) {
00124     switch ( mol->get_fep_bonded_type(atomID, 2) ) {
00125     case 1:
00126       reduction[bondEnergyIndex_ti_1] += energy;
00127       reduction[bondEnergyIndex_f] += (bond_lambda_12 - bond_lambda_1)*energy;
00128       energy *= bond_lambda_1;
00129       scal *= bond_lambda_1;
00130       break;
00131     case 2:
00132       reduction[bondEnergyIndex_ti_2] += energy;
00133       reduction[bondEnergyIndex_f] += (bond_lambda_22 - bond_lambda_2)*energy;
00134       energy *= bond_lambda_2;
00135       scal *= bond_lambda_2;
00136       break;
00137     }
00138   }
00139   //fepe
00140   const Force f12 = scal * r12;
00141 
00142 
00143   //  Now add the forces to each force vector
00144   p[0]->f[localIndex[0]] += f12;
00145   p[1]->f[localIndex[1]] -= f12;
00146 
00147   DebugM(3, "::computeForce() -- ending with delta energy " << energy << std::endl);
00148   reduction[bondEnergyIndex] += 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 (pressureProfileData) {
00160     BigReal z1 = p[0]->x[localIndex[0]].position.z;
00161     BigReal z2 = p[1]->x[localIndex[1]].position.z;
00162     int n1 = (int)floor((z1-pressureProfileMin)/pressureProfileThickness);
00163     int n2 = (int)floor((z2-pressureProfileMin)/pressureProfileThickness);
00164     pp_clamp(n1, pressureProfileSlabs);
00165     pp_clamp(n2, pressureProfileSlabs);
00166     int p1 = p[0]->x[localIndex[0]].partition;
00167     int p2 = p[1]->x[localIndex[1]].partition;
00168     int pn = pressureProfileAtomTypes;
00169     pp_reduction(pressureProfileSlabs,
00170                 n1, n2, 
00171                 p1, p2, pn, 
00172                 f12.x * r12.x, f12.y * r12.y, f12.z * r12.z,
00173                 pressureProfileData);
00174   } 
00175 
00176  }
00177 }
00178 
00179 void BondElem::submitReductionData(BigReal *data, SubmitReduction *reduction)
00180 {
00181   reduction->item(REDUCTION_BOND_ENERGY) += data[bondEnergyIndex];
00182   reduction->item(REDUCTION_BONDED_ENERGY_F) += data[bondEnergyIndex_f];
00183   reduction->item(REDUCTION_BONDED_ENERGY_TI_1) += data[bondEnergyIndex_ti_1];
00184   reduction->item(REDUCTION_BONDED_ENERGY_TI_2) += data[bondEnergyIndex_ti_2];
00185   ADD_TENSOR(reduction,REDUCTION_VIRIAL_NORMAL,data,virialIndex);
00186 }
00187 

Generated on Tue Sep 19 01:17:10 2017 for NAMD by  doxygen 1.4.7