--- trunk/src/brains/Thermo.cpp 2006/05/17 21:51:42 963 +++ branches/development/src/brains/Thermo.cpp 2012/06/05 17:48:40 1733 @@ -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 @@ -49,9 +50,10 @@ #include "brains/Thermo.hpp" #include "primitives/Molecule.hpp" #include "utils/simError.h" -#include "utils/OOPSEConstant.hpp" +#include "utils/PhysicalConstants.hpp" +#include "types/MultipoleAdapter.hpp" -namespace oopse { +namespace OpenMD { RealType Thermo::getKinetic() { SimInfo::MoleculeIterator miter; @@ -64,6 +66,7 @@ namespace oopse { int i; int j; int k; + RealType mass; RealType kinetic = 0.0; RealType kinetic_global = 0.0; @@ -71,8 +74,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 +105,7 @@ namespace oopse { #endif //is_mpi - kinetic = kinetic * 0.5 / OOPSEConstant::energyConvert; + kinetic = kinetic * 0.5 / PhysicalConstants::energyConvert; return kinetic; } @@ -138,10 +141,46 @@ 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; } + RealType Thermo::getElectronicTemperature() { + SimInfo::MoleculeIterator miter; + std::vector::iterator iiter; + Molecule* mol; + Atom* atom; + RealType cvel; + RealType cmass; + RealType kinetic = 0.0; + RealType kinetic_global = 0.0; + + for (mol = info_->beginMolecule(miter); mol != NULL; mol = info_->nextMolecule(miter)) { + for (atom = mol->beginFluctuatingCharge(iiter); atom != NULL; + atom = mol->nextFluctuatingCharge(iiter)) { + cmass = atom->getChargeMass(); + cvel = atom->getFlucQVel(); + + kinetic += cmass * cvel * cvel; + + } + } + +#ifdef IS_MPI + + MPI_Allreduce(&kinetic, &kinetic_global, 1, MPI_REALTYPE, MPI_SUM, + MPI_COMM_WORLD); + kinetic = kinetic_global; + +#endif //is_mpi + + kinetic = kinetic * 0.5; + return ( 2.0 * kinetic) / (info_->getNFluctuatingCharges()* PhysicalConstants::kb ); + } + + + + RealType Thermo::getVolume() { Snapshot* curSnapshot = info_->getSnapshotManager()->getCurrentSnapshot(); return curSnapshot->getVolume(); @@ -157,7 +196,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 +211,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: @@ -209,13 +246,15 @@ namespace oopse { RealType volume = this->getVolume(); Snapshot* curSnapshot = info_->getSnapshotManager()->getCurrentSnapshot(); - Mat3x3d tau = curSnapshot->statData.getTau(); + Mat3x3d stressTensor = curSnapshot->getStressTensor(); - pressureTensor = (p_global + OOPSEConstant::energyConvert* tau)/volume; - + pressureTensor = (p_global + + PhysicalConstants::energyConvert * stressTensor)/volume; + return pressureTensor; } + void Thermo::saveStat(){ Snapshot* currSnapshot = info_->getSnapshotManager()->getCurrentSnapshot(); Stats& stat = currSnapshot->statData; @@ -228,14 +267,303 @@ 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(); + // grab the heat flux if desired + if (simParams->havePrintHeatFlux()) { + if (simParams->getPrintHeatFlux()){ + Vector3d heatFlux = getHeatFlux(); + stat[Stats::HEATFLUX_X] = heatFlux(0); + stat[Stats::HEATFLUX_Y] = heatFlux(1); + stat[Stats::HEATFLUX_Z] = heatFlux(2); + } + } + + 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++; + } + } + } + + MultipoleAdapter ma = MultipoleAdapter(atom->getAtomType()); + if (ma.isDipole() ) { + Vector3d u_i = atom->getElectroFrame().getColumn(2); + moment = ma.getDipoleMoment(); + 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; + } + + // Returns the Heat Flux Vector for the system + Vector3d Thermo::getHeatFlux(){ + Snapshot* currSnapshot = info_->getSnapshotManager()->getCurrentSnapshot(); + SimInfo::MoleculeIterator miter; + std::vector::iterator iiter; + Molecule* mol; + StuntDouble* integrableObject; + RigidBody::AtomIterator ai; + Atom* atom; + Vector3d vel; + Vector3d angMom; + Mat3x3d I; + int i; + int j; + int k; + RealType mass; + + Vector3d x_a; + RealType kinetic; + RealType potential; + RealType eatom; + RealType AvgE_a_ = 0; + // Convective portion of the heat flux + Vector3d heatFluxJc = V3Zero; + + /* Calculate convective portion of the heat flux */ + for (mol = info_->beginMolecule(miter); mol != NULL; + mol = info_->nextMolecule(miter)) { + + for (integrableObject = mol->beginIntegrableObject(iiter); + integrableObject != NULL; + integrableObject = mol->nextIntegrableObject(iiter)) { + + mass = integrableObject->getMass(); + vel = integrableObject->getVel(); + + kinetic = mass * (vel[0]*vel[0] + vel[1]*vel[1] + vel[2]*vel[2]); + + if (integrableObject->isDirectional()) { + angMom = integrableObject->getJ(); + I = integrableObject->getI(); + + if (integrableObject->isLinear()) { + i = integrableObject->linearAxis(); + j = (i + 1) % 3; + k = (i + 2) % 3; + kinetic += angMom[j] * angMom[j] / I(j, j) + angMom[k] * angMom[k] / I(k, k); + } else { + kinetic += angMom[0]*angMom[0]/I(0, 0) + angMom[1]*angMom[1]/I(1, 1) + + angMom[2]*angMom[2]/I(2, 2); + } + } + + potential = 0.0; + + if (integrableObject->isRigidBody()) { + RigidBody* rb = dynamic_cast(integrableObject); + for (atom = rb->beginAtom(ai); atom != NULL; + atom = rb->nextAtom(ai)) { + potential += atom->getParticlePot(); + } + } else { + potential = integrableObject->getParticlePot(); + cerr << "ppot = " << potential << "\n"; + } + + potential *= PhysicalConstants::energyConvert; // amu A^2/fs^2 + // The potential may not be a 1/2 factor + eatom = (kinetic + potential)/2.0; // amu A^2/fs^2 + heatFluxJc[0] += eatom*vel[0]; // amu A^3/fs^3 + heatFluxJc[1] += eatom*vel[1]; // amu A^3/fs^3 + heatFluxJc[2] += eatom*vel[2]; // amu A^3/fs^3 + } + } + + std::cerr << "Heat flux heatFluxJc is: " << heatFluxJc << std::endl; + + /* The J_v vector is reduced in fortan so everyone has the global + * Jv. Jc is computed over the local atoms and must be reduced + * among all processors. + */ +#ifdef IS_MPI + MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &heatFluxJc[0], 3, MPI::REALTYPE, + MPI::SUM); +#endif + + // (kcal/mol * A/fs) * conversion => (amu A^3)/fs^3 + + Vector3d heatFluxJv = currSnapshot->getConductiveHeatFlux() * + PhysicalConstants::energyConvert; + + std::cerr << "Heat flux Jc is: " << heatFluxJc << std::endl; + std::cerr << "Heat flux Jv is: " << heatFluxJv << std::endl; + + // Correct for the fact the flux is 1/V (Jc + Jv) + return (heatFluxJv + heatFluxJc) / this->getVolume(); // amu / fs^3 + } +} //end namespace OpenMD