--- trunk/src/brains/Thermo.cpp 2006/05/17 21:51:42 963 +++ branches/development/src/brains/Thermo.cpp 2010/07/09 23:08:25 1465 @@ -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,15 @@ * 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] Vardeman & Gezelter, in progress (2009). */ #include @@ -49,9 +49,9 @@ #include "brains/Thermo.hpp" #include "primitives/Molecule.hpp" #include "utils/simError.h" -#include "utils/OOPSEConstant.hpp" +#include "utils/PhysicalConstants.hpp" -namespace oopse { +namespace OpenMD { RealType Thermo::getKinetic() { SimInfo::MoleculeIterator miter; @@ -64,6 +64,7 @@ namespace oopse { int i; int j; int k; + RealType mass; RealType kinetic = 0.0; RealType kinetic_global = 0.0; @@ -71,8 +72,8 @@ namespace oopse { for (integrableObject = mol->beginIntegrableObject(iiter); integrableObject != NULL; integrableObject = mol->nextIntegrableObject(iiter)) { - RealType mass = integrableObject->getMass(); - Vector3d vel = integrableObject->getVel(); + mass = integrableObject->getMass(); + vel = integrableObject->getVel(); kinetic += mass * (vel[0]*vel[0] + vel[1]*vel[1] + vel[2]*vel[2]); @@ -102,7 +103,7 @@ namespace oopse { #endif //is_mpi - kinetic = kinetic * 0.5 / OOPSEConstant::energyConvert; + kinetic = kinetic * 0.5 / PhysicalConstants::energyConvert; return kinetic; } @@ -138,7 +139,7 @@ namespace oopse { RealType Thermo::getTemperature() { - RealType temperature = ( 2.0 * this->getKinetic() ) / (info_->getNdf()* OOPSEConstant::kb ); + RealType temperature = ( 2.0 * this->getKinetic() ) / (info_->getNdf()* PhysicalConstants::kb ); return temperature; } @@ -157,7 +158,7 @@ namespace oopse { tensor = getPressureTensor(); - pressure = OOPSEConstant::pressureConvert * (tensor(0, 0) + tensor(1, 1) + tensor(2, 2)) / 3.0; + pressure = PhysicalConstants::pressureConvert * (tensor(0, 0) + tensor(1, 1) + tensor(2, 2)) / 3.0; return pressure; } @@ -172,13 +173,11 @@ namespace oopse { tensor = getPressureTensor(); - pressure = OOPSEConstant::pressureConvert * tensor(direction, direction); + pressure = PhysicalConstants::pressureConvert * tensor(direction, direction); return pressure; } - - Mat3x3d Thermo::getPressureTensor() { // returns pressure tensor in units amu*fs^-2*Ang^-1 // routine derived via viral theorem description in: @@ -211,11 +210,12 @@ namespace oopse { Snapshot* curSnapshot = info_->getSnapshotManager()->getCurrentSnapshot(); Mat3x3d tau = curSnapshot->statData.getTau(); - pressureTensor = (p_global + OOPSEConstant::energyConvert* tau)/volume; - + pressureTensor = (p_global + PhysicalConstants::energyConvert* tau)/volume; + return pressureTensor; } + void Thermo::saveStat(){ Snapshot* currSnapshot = info_->getSnapshotManager()->getCurrentSnapshot(); Stats& stat = currSnapshot->statData; @@ -228,14 +228,80 @@ namespace oopse { stat[Stats::VOLUME] = getVolume(); Mat3x3d tensor =getPressureTensor(); - stat[Stats::PRESSURE_TENSOR_X] = tensor(0, 0); - stat[Stats::PRESSURE_TENSOR_Y] = tensor(1, 1); - stat[Stats::PRESSURE_TENSOR_Z] = tensor(2, 2); + stat[Stats::PRESSURE_TENSOR_XX] = tensor(0, 0); + stat[Stats::PRESSURE_TENSOR_XY] = tensor(0, 1); + stat[Stats::PRESSURE_TENSOR_XZ] = tensor(0, 2); + stat[Stats::PRESSURE_TENSOR_YX] = tensor(1, 0); + stat[Stats::PRESSURE_TENSOR_YY] = tensor(1, 1); + stat[Stats::PRESSURE_TENSOR_YZ] = tensor(1, 2); + stat[Stats::PRESSURE_TENSOR_ZX] = tensor(2, 0); + stat[Stats::PRESSURE_TENSOR_ZY] = tensor(2, 1); + stat[Stats::PRESSURE_TENSOR_ZZ] = tensor(2, 2); + + + Globals* simParams = info_->getSimParams(); + if (simParams->haveTaggedAtomPair() && + simParams->havePrintTaggedPairDistance()) { + if ( simParams->getPrintTaggedPairDistance()) { + + std::pair tap = simParams->getTaggedAtomPair(); + Vector3d pos1, pos2, rab; +#ifdef IS_MPI + std::cerr << "tap = " << tap.first << " " << tap.second << std::endl; + + int mol1 = info_->getGlobalMolMembership(tap.first); + int mol2 = info_->getGlobalMolMembership(tap.second); + std::cerr << "mols = " << mol1 << " " << mol2 << std::endl; + + int proc1 = info_->getMolToProc(mol1); + int proc2 = info_->getMolToProc(mol2); + + std::cerr << " procs = " << proc1 << " " <getIOIndexToIntegrableObject(tap.first); + std::cerr << " on proc " << proc1 << ", sd1 has global index= " << sd1->getGlobalIndex() << std::endl; + pos1 = sd1->getPos(); + data[0] = pos1.x(); + data[1] = pos1.y(); + data[2] = pos1.z(); + MPI_Bcast(data, 3, MPI_REALTYPE, proc1, MPI_COMM_WORLD); + } else { + MPI_Bcast(data, 3, MPI_REALTYPE, proc1, MPI_COMM_WORLD); + pos1 = Vector3d(data); + } + + + if (proc2 == worldRank) { + StuntDouble* sd2 = info_->getIOIndexToIntegrableObject(tap.second); + std::cerr << " on proc " << proc2 << ", sd2 has global index= " << sd2->getGlobalIndex() << std::endl; + pos2 = sd2->getPos(); + data[0] = pos2.x(); + data[1] = pos2.y(); + data[2] = pos2.z(); + MPI_Bcast(data, 3, MPI_REALTYPE, proc2, MPI_COMM_WORLD); + } else { + MPI_Bcast(data, 3, MPI_REALTYPE, proc2, MPI_COMM_WORLD); + pos2 = Vector3d(data); + } +#else + StuntDouble* at1 = info_->getIOIndexToIntegrableObject(tap.first); + StuntDouble* at2 = info_->getIOIndexToIntegrableObject(tap.second); + pos1 = at1->getPos(); + pos2 = at2->getPos(); +#endif + rab = pos2 - pos1; + currSnapshot->wrapVector(rab); + stat[Stats::TAGGED_PAIR_DISTANCE] = rab.length(); + } + } + /**@todo need refactorying*/ //Conserved Quantity is set by integrator and time is set by setTime } -} //end namespace oopse +} //end namespace OpenMD