--- trunk/src/integrators/Velocitizer.cpp 2006/10/19 15:57:07 1079 +++ trunk/src/integrators/Velocitizer.cpp 2012/08/22 02:28:28 1782 @@ -6,19 +6,10 @@ * redistribute this software in source and binary code form, provided * that the following conditions are met: * - * 1. Acknowledgement of the program authors must be made in any - * publication of scientific results based in part on use of the - * program. An acceptable form of acknowledgement is citation of - * the article in which the program was described (Matthew - * A. Meineke, Charles F. Vardeman II, Teng Lin, Christopher - * J. Fennell and J. Daniel Gezelter, "OOPSE: An Object-Oriented - * Parallel Simulation Engine for Molecular Dynamics," - * J. Comput. Chem. 26, pp. 252-271 (2005)) - * - * 2. Redistributions of source code must retain the above copyright + * 1. Redistributions of source code must retain the above copyright * notice, this list of conditions and the following disclaimer. * - * 3. Redistributions in binary form must reproduce the above copyright + * 2. Redistributions in binary form must reproduce the above copyright * notice, this list of conditions and the following disclaimer in the * documentation and/or other materials provided with the * distribution. @@ -37,6 +28,16 @@ * arising out of the use of or inability to use software, even if the * University of Notre Dame has been advised of the possibility of * such damages. + * + * SUPPORT OPEN SCIENCE! If you use OpenMD or its source code in your + * research, please cite the appropriate papers when you publish your + * work. Good starting points are: + * + * [1] Meineke, et al., J. Comp. Chem. 26, 252-271 (2005). + * [2] Fennell & Gezelter, J. Chem. Phys. 124, 234104 (2006). + * [3] Sun, Lin & Gezelter, J. Chem. Phys. 128, 24107 (2008). + * [4] Kuang & Gezelter, J. Chem. Phys. 133, 164101 (2010). + * [5] Vardeman, Stocker & Gezelter, J. Chem. Theory Comput. 7, 834 (2011). */ #include "integrators/Velocitizer.hpp" @@ -50,14 +51,9 @@ #include "math/ParallelRandNumGen.hpp" #endif -/* Remove me after testing*/ -#include -#include -/*End remove me*/ - -namespace oopse { +namespace OpenMD { - Velocitizer::Velocitizer(SimInfo* info) : info_(info) { + Velocitizer::Velocitizer(SimInfo* info) : info_(info), thermo(info) { int seedValue; Globals * simParams = info->getSimParams(); @@ -87,12 +83,10 @@ namespace oopse { Vector3d aVel; Vector3d aJ; Mat3x3d I; - int l; - int m; - int n; + int l, m, n; Vector3d vdrift; RealType vbar; - /**@todo refactory kb */ + /**@todo refactor kb */ const RealType kb = 8.31451e-7; // kb in amu, angstroms, fs, etc. RealType av2; RealType kebar; @@ -102,18 +96,19 @@ namespace oopse { SimInfo::MoleculeIterator i; Molecule::IntegrableObjectIterator j; Molecule * mol; - StuntDouble * integrableObject; - + StuntDouble * sd; + kebar = kb * temperature * info_->getNdfRaw() / (2.0 * info_->getNdf()); + for( mol = info_->beginMolecule(i); mol != NULL; mol = info_->nextMolecule(i) ) { - for( integrableObject = mol->beginIntegrableObject(j); - integrableObject != NULL; - integrableObject = mol->nextIntegrableObject(j) ) { + + for( sd = mol->beginIntegrableObject(j); sd != NULL; + sd = mol->nextIntegrableObject(j) ) { // uses equipartition theory to solve for vbar in angstrom/fs - av2 = 2.0 * kebar / integrableObject->getMass(); + av2 = 2.0 * kebar / sd->getMass(); vbar = sqrt(av2); // picks random velocities from a gaussian distribution @@ -122,13 +117,13 @@ namespace oopse { for( int k = 0; k < 3; k++ ) { aVel[k] = vbar * randNumGen_->randNorm(0.0, 1.0); } - integrableObject->setVel(aVel); + sd->setVel(aVel); - if (integrableObject->isDirectional()) { - I = integrableObject->getI(); + if (sd->isDirectional()) { + I = sd->getI(); - if (integrableObject->isLinear()) { - l = integrableObject->linearAxis(); + if (sd->isLinear()) { + l = sd->linearAxis(); m = (l + 1) % 3; n = (l + 2) % 3; @@ -142,102 +137,81 @@ namespace oopse { vbar = sqrt(2.0 * kebar * I(k, k)); aJ[k] = vbar *randNumGen_->randNorm(0.0, 1.0); } - } // else isLinear + } - integrableObject->setJ(aJ); - } //isDirectional + sd->setJ(aJ); + } } - } //end for (mol = beginMolecule(i); ...) - - - + } + removeComDrift(); - // Remove angular drift if we are not using periodic boundary conditions. - if(!simParams->getUsePeriodicBoundaryConditions()) removeAngularDrift(); - + + // Remove angular drift if we are not using periodic boundary + // conditions: + + if(!simParams->getUsePeriodicBoundaryConditions()) removeAngularDrift(); } - - - + void Velocitizer::removeComDrift() { // Get the Center of Mass drift velocity. - Vector3d vdrift = info_->getComVel(); + Vector3d vdrift = thermo.getComVel(); SimInfo::MoleculeIterator i; Molecule::IntegrableObjectIterator j; Molecule * mol; - StuntDouble * integrableObject; + StuntDouble * sd; // Corrects for the center of mass drift. // sums all the momentum and divides by total mass. for( mol = info_->beginMolecule(i); mol != NULL; - mol = info_->nextMolecule(i) ) { - for( integrableObject = mol->beginIntegrableObject(j); - integrableObject != NULL; - integrableObject = mol->nextIntegrableObject(j) ) { - integrableObject->setVel(integrableObject->getVel() - vdrift); + mol = info_->nextMolecule(i) ) { + + for( sd = mol->beginIntegrableObject(j); sd != NULL; + sd = mol->nextIntegrableObject(j) ) { + + sd->setVel(sd->getVel() - vdrift); + } - } - + } } - - + void Velocitizer::removeAngularDrift() { // Get the Center of Mass drift velocity. Vector3d vdrift; Vector3d com; - info_->getComAll(com,vdrift); + thermo.getComAll(com, vdrift); Mat3x3d inertiaTensor; Vector3d angularMomentum; Vector3d omega; - - - - info_->getInertiaTensor(inertiaTensor,angularMomentum); + + thermo.getInertiaTensor(inertiaTensor, angularMomentum); + // We now need the inverse of the inertia tensor. - /* - std::cerr << "Angular Momentum before is " - << angularMomentum << std::endl; - std::cerr << "Inertia Tensor before is " - << inertiaTensor << std::endl; - */ - inertiaTensor =inertiaTensor.inverse(); - /* - std::cerr << "Inertia Tensor after inverse is " - << inertiaTensor << std::endl; - */ - omega = inertiaTensor*angularMomentum; - + inertiaTensor = inertiaTensor.inverse(); + omega = inertiaTensor * angularMomentum; + SimInfo::MoleculeIterator i; Molecule::IntegrableObjectIterator j; - Molecule * mol; - StuntDouble * integrableObject; + Molecule* mol; + StuntDouble* sd; Vector3d tempComPos; - - // Corrects for the center of mass angular drift. - // sums all the angular momentum and divides by total mass. + + // Corrects for the center of mass angular drift by summing all + // the angular momentum and dividing by the total mass. + for( mol = info_->beginMolecule(i); mol != NULL; mol = info_->nextMolecule(i) ) { - for( integrableObject = mol->beginIntegrableObject(j); - integrableObject != NULL; - integrableObject = mol->nextIntegrableObject(j) ) { - tempComPos = integrableObject->getPos()-com; - integrableObject->setVel((integrableObject->getVel() - vdrift)-cross(omega,tempComPos)); + + for( sd = mol->beginIntegrableObject(j); sd != NULL; + sd = mol->nextIntegrableObject(j) ) { + + tempComPos = sd->getPos() - com; + sd->setVel((sd->getVel() - vdrift) - cross(omega, tempComPos)); + } - } - - angularMomentum = info_->getAngularMomentum(); - /* - std::cerr << "Angular Momentum after is " - << angularMomentum << std::endl; - */ - - } - - - - + } + } }