--- trunk/src/brains/Thermo.cpp 2006/05/17 21:51:42 963 +++ branches/development/src/brains/Thermo.cpp 2010/10/02 19:54:41 1503 @@ -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,199 @@ 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); + // grab the simulation box dipole moment if specified + if (info_->getCalcBoxDipole()){ + Vector3d totalDipole = getBoxDipole(); + stat[Stats::BOX_DIPOLE_X] = totalDipole(0); + stat[Stats::BOX_DIPOLE_Y] = totalDipole(1); + stat[Stats::BOX_DIPOLE_Z] = totalDipole(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 + + Vector3d Thermo::getBoxDipole() { + Snapshot* currSnapshot = info_->getSnapshotManager()->getCurrentSnapshot(); + SimInfo::MoleculeIterator miter; + std::vector::iterator aiter; + Molecule* mol; + Atom* atom; + RealType charge; + RealType moment(0.0); + Vector3d ri(0.0); + Vector3d dipoleVector(0.0); + Vector3d nPos(0.0); + Vector3d pPos(0.0); + RealType nChg(0.0); + RealType pChg(0.0); + int nCount = 0; + int pCount = 0; + + RealType chargeToC = 1.60217733e-19; + RealType angstromToM = 1.0e-10; + RealType debyeToCm = 3.33564095198e-30; + + for (mol = info_->beginMolecule(miter); mol != NULL; + mol = info_->nextMolecule(miter)) { + + for (atom = mol->beginAtom(aiter); atom != NULL; + atom = mol->nextAtom(aiter)) { + + if (atom->isCharge() ) { + charge = 0.0; + GenericData* data = atom->getAtomType()->getPropertyByName("Charge"); + if (data != NULL) { + + charge = (dynamic_cast(data))->getData(); + charge *= chargeToC; + + ri = atom->getPos(); + currSnapshot->wrapVector(ri); + ri *= angstromToM; + + if (charge < 0.0) { + nPos += ri; + nChg -= charge; + nCount++; + } else if (charge > 0.0) { + pPos += ri; + pChg += charge; + pCount++; + } + } + } + + if (atom->isDipole() ) { + Vector3d u_i = atom->getElectroFrame().getColumn(2); + GenericData* data = dynamic_cast(atom->getAtomType())->getPropertyByName("Dipole"); + if (data != NULL) { + moment = (dynamic_cast(data))->getData(); + + moment *= debyeToCm; + dipoleVector += u_i * moment; + } + } + } + } + + +#ifdef IS_MPI + RealType pChg_global, nChg_global; + int pCount_global, nCount_global; + Vector3d pPos_global, nPos_global, dipVec_global; + + MPI_Allreduce(&pChg, &pChg_global, 1, MPI_REALTYPE, MPI_SUM, + MPI_COMM_WORLD); + pChg = pChg_global; + MPI_Allreduce(&nChg, &nChg_global, 1, MPI_REALTYPE, MPI_SUM, + MPI_COMM_WORLD); + nChg = nChg_global; + MPI_Allreduce(&pCount, &pCount_global, 1, MPI_INTEGER, MPI_SUM, + MPI_COMM_WORLD); + pCount = pCount_global; + MPI_Allreduce(&nCount, &nCount_global, 1, MPI_INTEGER, MPI_SUM, + MPI_COMM_WORLD); + nCount = nCount_global; + MPI_Allreduce(pPos.getArrayPointer(), pPos_global.getArrayPointer(), 3, + MPI_REALTYPE, MPI_SUM, MPI_COMM_WORLD); + pPos = pPos_global; + MPI_Allreduce(nPos.getArrayPointer(), nPos_global.getArrayPointer(), 3, + MPI_REALTYPE, MPI_SUM, MPI_COMM_WORLD); + nPos = nPos_global; + MPI_Allreduce(dipoleVector.getArrayPointer(), + dipVec_global.getArrayPointer(), 3, + MPI_REALTYPE, MPI_SUM, MPI_COMM_WORLD); + dipoleVector = dipVec_global; +#endif //is_mpi + + // first load the accumulated dipole moment (if dipoles were present) + Vector3d boxDipole = dipoleVector; + // now include the dipole moment due to charges + // use the lesser of the positive and negative charge totals + RealType chg_value = nChg <= pChg ? nChg : pChg; + + // find the average positions + if (pCount > 0 && nCount > 0 ) { + pPos /= pCount; + nPos /= nCount; + } + + // dipole is from the negative to the positive (physics notation) + boxDipole += (pPos - nPos) * chg_value; + + return boxDipole; + } +} //end namespace OpenMD