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

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(BigReal *reduction, 
00041                             BigReal *pressureProfileData)
00042 {
00043   DebugM(1, "::computeForce() localIndex = " << localIndex[0] << " "
00044                << localIndex[1] << std::endl);
00045 
00046   // skip Lonepair bonds (other k=0. bonds have been filtered out)
00047   if (0. == value->k) return;
00048 
00049   BigReal scal;         // force scaling
00050   BigReal energy;       // energy from the bond
00051 
00052   //DebugM(3, "::computeForce() -- starting with bond type " << bondType << std::endl);
00053 
00054   // get the bond information
00055   Real k = value->k * scale;
00056   Real x0 = value->x0;
00057 
00058   // compute vectors between atoms and their distances
00059   const Lattice & lattice = p[0]->p->lattice;
00060   const Vector r12 = lattice.delta(p[0]->x[localIndex[0]].position,
00061                                         p[1]->x[localIndex[1]].position);
00062 
00063   if (0. == x0) {  // for Drude bonds
00064     SimParameters *simParams = Node::Object()->simParameters;
00065     BigReal drudeBondLen = simParams->drudeBondLen;
00066     BigReal drudeBondConst = simParams->drudeBondConst;
00067 
00068     BigReal r2 = r12.length2();
00069 
00070     scal = -2.0*k;    // bond interaction for equilibrium length 0
00071     energy = k * r2;
00072 
00073     if ( (drudeBondConst > 0) && ( ! simParams->drudeHardWallOn ) &&
00074         (r2 > drudeBondLen*drudeBondLen) ) {
00075       // add a quartic restraining potential to keep Drude bond short
00076       BigReal r = sqrt(r2);
00077       BigReal diff = r - drudeBondLen;
00078       BigReal diff2 = diff*diff;
00079 
00080       scal += -4*drudeBondConst * diff2 * diff / r;
00081       energy += drudeBondConst * diff2 * diff2;
00082     }
00083   }
00084   else {
00085     BigReal r = r12.length();  // Distance between atoms
00086     BigReal diff = r - x0;     // Compare it to the rest bond
00087 
00088     //  Add the energy from this bond to the total energy
00089     energy = k*diff*diff;
00090 
00091     //  Determine the magnitude of the force
00092     diff *= -2.0*k;
00093 
00094     //  Scale the force vector accordingly
00095     scal = (diff/r);
00096   }
00097   const Force f12 = scal * r12;
00098 
00099 
00100   //  Now add the forces to each force vector
00101   p[0]->f[localIndex[0]] += f12;
00102   p[1]->f[localIndex[1]] -= f12;
00103 
00104   DebugM(3, "::computeForce() -- ending with delta energy " << energy << std::endl);
00105   reduction[bondEnergyIndex] += energy;
00106   reduction[virialIndex_XX] += f12.x * r12.x;
00107   reduction[virialIndex_XY] += f12.x * r12.y;
00108   reduction[virialIndex_XZ] += f12.x * r12.z;
00109   reduction[virialIndex_YX] += f12.y * r12.x;
00110   reduction[virialIndex_YY] += f12.y * r12.y;
00111   reduction[virialIndex_YZ] += f12.y * r12.z;
00112   reduction[virialIndex_ZX] += f12.z * r12.x;
00113   reduction[virialIndex_ZY] += f12.z * r12.y;
00114   reduction[virialIndex_ZZ] += f12.z * r12.z;
00115 
00116   if (pressureProfileData) {
00117     BigReal z1 = p[0]->x[localIndex[0]].position.z;
00118     BigReal z2 = p[1]->x[localIndex[1]].position.z;
00119     int n1 = (int)floor((z1-pressureProfileMin)/pressureProfileThickness);
00120     int n2 = (int)floor((z2-pressureProfileMin)/pressureProfileThickness);
00121     pp_clamp(n1, pressureProfileSlabs);
00122     pp_clamp(n2, pressureProfileSlabs);
00123     int p1 = p[0]->x[localIndex[0]].partition;
00124     int p2 = p[1]->x[localIndex[1]].partition;
00125     int pn = pressureProfileAtomTypes;
00126     pp_reduction(pressureProfileSlabs,
00127                 n1, n2, 
00128                 p1, p2, pn, 
00129                 f12.x * r12.x, f12.y * r12.y, f12.z * r12.z,
00130                 pressureProfileData);
00131   } 
00132 }
00133 
00134 void BondElem::submitReductionData(BigReal *data, SubmitReduction *reduction)
00135 {
00136   reduction->item(REDUCTION_BOND_ENERGY) += data[bondEnergyIndex];
00137   ADD_TENSOR(reduction,REDUCTION_VIRIAL_NORMAL,data,virialIndex);
00138 }
00139 

Generated on Tue May 21 04:07:14 2013 for NAMD by  doxygen 1.3.9.1