--- trunk/src/brains/Thermo.cpp 2004/09/24 16:27:58 3 +++ branches/development/src/brains/Thermo.cpp 2013/05/15 15:09:35 1874 @@ -1,450 +1,970 @@ +/* + * Copyright (c) 2005 The University of Notre Dame. All Rights Reserved. + * + * The University of Notre Dame grants you ("Licensee") a + * non-exclusive, royalty free, license to use, modify and + * redistribute this software in source and binary code form, provided + * that the following conditions are met: + * + * 1. Redistributions of source code must retain the above copyright + * notice, this list of conditions and the following disclaimer. + * + * 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. + * + * This software is provided "AS IS," without a warranty of any + * kind. All express or implied conditions, representations and + * warranties, including any implied warranty of merchantability, + * fitness for a particular purpose or non-infringement, are hereby + * excluded. The University of Notre Dame and its licensors shall not + * be liable for any damages suffered by licensee as a result of + * using, modifying or distributing the software or its + * derivatives. In no event will the University of Notre Dame or its + * licensors be liable for any lost revenue, profit or data, or for + * direct, indirect, special, consequential, incidental or punitive + * damages, however caused and regardless of the theory of liability, + * 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, 234107 (2008). + * [4] Kuang & Gezelter, J. Chem. Phys. 133, 164101 (2010). + * [5] Vardeman, Stocker & Gezelter, J. Chem. Theory Comput. 7, 834 (2011). + */ + #include #include -using namespace std; #ifdef IS_MPI #include #endif //is_mpi #include "brains/Thermo.hpp" -#include "primitives/SRI.hpp" -#include "integrators/Integrator.hpp" +#include "primitives/Molecule.hpp" #include "utils/simError.h" -#include "math/MatVec3.h" +#include "utils/PhysicalConstants.hpp" +#include "types/FixedChargeAdapter.hpp" +#include "types/FluctuatingChargeAdapter.hpp" +#include "types/MultipoleAdapter.hpp" +#ifdef HAVE_QHULL +#include "math/ConvexHull.hpp" +#include "math/AlphaHull.hpp" +#endif -#ifdef IS_MPI -#define __C -#include "brains/mpiSimulation.hpp" -#endif // is_mpi +using namespace std; +namespace OpenMD { -inline double roundMe( double x ){ - return ( x >= 0 ) ? floor( x + 0.5 ) : ceil( x - 0.5 ); -} + RealType Thermo::getTranslationalKinetic() { + Snapshot* snap = info_->getSnapshotManager()->getCurrentSnapshot(); -Thermo::Thermo( SimInfo* the_info ) { - info = the_info; - int baseSeed = the_info->getSeed(); - - gaussStream = new gaussianSPRNG( baseSeed ); -} + if (!snap->hasTranslationalKineticEnergy) { + SimInfo::MoleculeIterator miter; + vector::iterator iiter; + Molecule* mol; + StuntDouble* sd; + Vector3d vel; + RealType mass; + RealType kinetic(0.0); + + for (mol = info_->beginMolecule(miter); mol != NULL; + mol = info_->nextMolecule(miter)) { + + for (sd = mol->beginIntegrableObject(iiter); sd != NULL; + sd = mol->nextIntegrableObject(iiter)) { + + mass = sd->getMass(); + vel = sd->getVel(); + + kinetic += mass * (vel[0]*vel[0] + vel[1]*vel[1] + vel[2]*vel[2]); + + } + } + +#ifdef IS_MPI + MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &kinetic, 1, MPI::REALTYPE, + MPI::SUM); +#endif + + kinetic = kinetic * 0.5 / PhysicalConstants::energyConvert; + + + snap->setTranslationalKineticEnergy(kinetic); + } + return snap->getTranslationalKineticEnergy(); + } -Thermo::~Thermo(){ - delete gaussStream; -} + RealType Thermo::getRotationalKinetic() { + Snapshot* snap = info_->getSnapshotManager()->getCurrentSnapshot(); -double Thermo::getKinetic(){ + if (!snap->hasRotationalKineticEnergy) { + SimInfo::MoleculeIterator miter; + vector::iterator iiter; + Molecule* mol; + StuntDouble* sd; + Vector3d angMom; + Mat3x3d I; + int i, j, k; + RealType kinetic(0.0); + + for (mol = info_->beginMolecule(miter); mol != NULL; + mol = info_->nextMolecule(miter)) { + + for (sd = mol->beginIntegrableObject(iiter); sd != NULL; + sd = mol->nextIntegrableObject(iiter)) { + + if (sd->isDirectional()) { + angMom = sd->getJ(); + I = sd->getI(); + + if (sd->isLinear()) { + i = sd->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); + } + } + } + } + +#ifdef IS_MPI + MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &kinetic, 1, MPI::REALTYPE, + MPI::SUM); +#endif + + kinetic = kinetic * 0.5 / PhysicalConstants::energyConvert; + + snap->setRotationalKineticEnergy(kinetic); + } + return snap->getRotationalKineticEnergy(); + } - const double e_convert = 4.184E-4; // convert kcal/mol -> (amu A^2)/fs^2 - double kinetic; - double amass; - double aVel[3], aJ[3], I[3][3]; - int i, j, k, kl; + - double kinetic_global; - vector integrableObjects = info->integrableObjects; - - kinetic = 0.0; - kinetic_global = 0.0; + RealType Thermo::getKinetic() { + Snapshot* snap = info_->getSnapshotManager()->getCurrentSnapshot(); - for (kl=0; klgetVel(aVel); - amass = integrableObjects[kl]->getMass(); + if (!snap->hasKineticEnergy) { + RealType ke = getTranslationalKinetic() + getRotationalKinetic(); + snap->setKineticEnergy(ke); + } + return snap->getKineticEnergy(); + } - for(j=0; j<3; j++) - kinetic += amass*aVel[j]*aVel[j]; + RealType Thermo::getPotential() { - if (integrableObjects[kl]->isDirectional()){ - - integrableObjects[kl]->getJ( aJ ); - integrableObjects[kl]->getI( I ); + // ForceManager computes the potential and stores it in the + // Snapshot. All we have to do is report it. - if (integrableObjects[kl]->isLinear()) { - i = integrableObjects[kl]->linearAxis(); - j = (i+1)%3; - k = (i+2)%3; - kinetic += aJ[j]*aJ[j]/I[j][j] + aJ[k]*aJ[k]/I[k][k]; - } else { - for (j=0; j<3; j++) - kinetic += aJ[j]*aJ[j] / I[j][j]; - } - } + Snapshot* snap = info_->getSnapshotManager()->getCurrentSnapshot(); + return snap->getPotentialEnergy(); } -#ifdef IS_MPI - MPI_Allreduce(&kinetic,&kinetic_global,1,MPI_DOUBLE, - MPI_SUM, MPI_COMM_WORLD); - kinetic = kinetic_global; -#endif //is_mpi - - kinetic = kinetic * 0.5 / e_convert; - return kinetic; -} + RealType Thermo::getTotalEnergy() { -double Thermo::getPotential(){ - - double potential_local; - double potential; - int el, nSRI; - Molecule* molecules; + Snapshot* snap = info_->getSnapshotManager()->getCurrentSnapshot(); - molecules = info->molecules; - nSRI = info->n_SRI; + if (!snap->hasTotalEnergy) { + snap->setTotalEnergy(this->getKinetic() + this->getPotential()); + } - potential_local = 0.0; - potential = 0.0; - potential_local += info->lrPot; - - for( el=0; eln_mol; el++ ){ - potential_local += molecules[el].getPotential(); + return snap->getTotalEnergy(); } - // Get total potential for entire system from MPI. -#ifdef IS_MPI - MPI_Allreduce(&potential_local,&potential,1,MPI_DOUBLE, - MPI_SUM, MPI_COMM_WORLD); -#else - potential = potential_local; -#endif // is_mpi + RealType Thermo::getTemperature() { - return potential; -} + Snapshot* snap = info_->getSnapshotManager()->getCurrentSnapshot(); -double Thermo::getTotalE(){ + if (!snap->hasTemperature) { - double total; + RealType temperature = ( 2.0 * this->getKinetic() ) + / (info_->getNdf()* PhysicalConstants::kb ); - total = this->getKinetic() + this->getPotential(); - return total; -} + snap->setTemperature(temperature); + } + + return snap->getTemperature(); + } -double Thermo::getTemperature(){ + RealType Thermo::getElectronicTemperature() { + Snapshot* snap = info_->getSnapshotManager()->getCurrentSnapshot(); - const double kb = 1.9872156E-3; // boltzman's constant in kcal/(mol K) - double temperature; + if (!snap->hasElectronicTemperature) { + + SimInfo::MoleculeIterator miter; + vector::iterator iiter; + Molecule* mol; + Atom* atom; + RealType cvel; + RealType cmass; + RealType kinetic(0.0); + RealType eTemp; + + 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::COMM_WORLD.Allreduce(MPI::IN_PLACE, &kinetic, 1, MPI::REALTYPE, + MPI::SUM); +#endif - temperature = ( 2.0 * this->getKinetic() ) / ((double)info->ndf * kb ); - return temperature; -} + kinetic *= 0.5; + eTemp = (2.0 * kinetic) / + (info_->getNFluctuatingCharges() * PhysicalConstants::kb ); + + snap->setElectronicTemperature(eTemp); + } -double Thermo::getVolume() { + return snap->getElectronicTemperature(); + } - return info->boxVol; -} -double Thermo::getPressure() { + RealType Thermo::getVolume() { + Snapshot* snap = info_->getSnapshotManager()->getCurrentSnapshot(); + return snap->getVolume(); + } - // Relies on the calculation of the full molecular pressure tensor - - const double p_convert = 1.63882576e8; - double press[3][3]; - double pressure; + RealType Thermo::getPressure() { + Snapshot* snap = info_->getSnapshotManager()->getCurrentSnapshot(); - this->getPressureTensor(press); + if (!snap->hasPressure) { + // Relies on the calculation of the full molecular pressure tensor + + Mat3x3d tensor; + RealType pressure; + + tensor = getPressureTensor(); + + pressure = PhysicalConstants::pressureConvert * + (tensor(0, 0) + tensor(1, 1) + tensor(2, 2)) / 3.0; + + snap->setPressure(pressure); + } + + return snap->getPressure(); + } - pressure = p_convert * (press[0][0] + press[1][1] + press[2][2]) / 3.0; + Mat3x3d Thermo::getPressureTensor() { + // returns pressure tensor in units amu*fs^-2*Ang^-1 + // routine derived via viral theorem description in: + // Paci, E. and Marchi, M. J.Phys.Chem. 1996, 100, 4314-4322 + Snapshot* snap = info_->getSnapshotManager()->getCurrentSnapshot(); - return pressure; -} + if (!snap->hasPressureTensor) { + + Mat3x3d pressureTensor; + Mat3x3d p_tens(0.0); + RealType mass; + Vector3d vcom; + + SimInfo::MoleculeIterator i; + vector::iterator j; + Molecule* mol; + StuntDouble* sd; + for (mol = info_->beginMolecule(i); mol != NULL; + mol = info_->nextMolecule(i)) { + + for (sd = mol->beginIntegrableObject(j); sd != NULL; + sd = mol->nextIntegrableObject(j)) { + + mass = sd->getMass(); + vcom = sd->getVel(); + p_tens += mass * outProduct(vcom, vcom); + } + } + +#ifdef IS_MPI + MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, p_tens.getArrayPointer(), 9, + MPI::REALTYPE, MPI::SUM); +#endif + + RealType volume = this->getVolume(); + Mat3x3d stressTensor = snap->getStressTensor(); + + pressureTensor = (p_tens + + PhysicalConstants::energyConvert * stressTensor)/volume; + + snap->setPressureTensor(pressureTensor); + } + return snap->getPressureTensor(); + } -double Thermo::getPressureX() { - // Relies on the calculation of the full molecular pressure tensor - - const double p_convert = 1.63882576e8; - double press[3][3]; - double pressureX; - this->getPressureTensor(press); - pressureX = p_convert * press[0][0]; + Vector3d Thermo::getSystemDipole() { + Snapshot* snap = info_->getSnapshotManager()->getCurrentSnapshot(); - return pressureX; -} + if (!snap->hasSystemDipole) { + SimInfo::MoleculeIterator miter; + vector::iterator aiter; + Molecule* mol; + Atom* atom; + RealType charge; + 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)) { + + charge = 0.0; + + FixedChargeAdapter fca = FixedChargeAdapter(atom->getAtomType()); + if ( fca.isFixedCharge() ) { + charge = fca.getCharge(); + } + + FluctuatingChargeAdapter fqa = FluctuatingChargeAdapter(atom->getAtomType()); + if ( fqa.isFluctuatingCharge() ) { + charge += atom->getFlucQPos(); + } + + charge *= chargeToC; + + ri = atom->getPos(); + snap->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()) { + dipoleVector += atom->getDipole() * debyeToCm; + } + } + } + + +#ifdef IS_MPI + MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &pChg, 1, MPI::REALTYPE, + MPI::SUM); + MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &nChg, 1, MPI::REALTYPE, + MPI::SUM); -double Thermo::getPressureY() { + MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &pCount, 1, MPI::INTEGER, + MPI::SUM); + MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &nCount, 1, MPI::INTEGER, + MPI::SUM); - // Relies on the calculation of the full molecular pressure tensor - - const double p_convert = 1.63882576e8; - double press[3][3]; - double pressureY; + MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, pPos.getArrayPointer(), 3, + MPI::REALTYPE, MPI::SUM); + MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, nPos.getArrayPointer(), 3, + MPI::REALTYPE, MPI::SUM); - this->getPressureTensor(press); + MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, dipoleVector.getArrayPointer(), + 3, MPI::REALTYPE, MPI::SUM); +#endif + + // 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; + snap->setSystemDipole(boxDipole); + } - pressureY = p_convert * press[1][1]; + return snap->getSystemDipole(); + } - return pressureY; -} + // Returns the Heat Flux Vector for the system + Vector3d Thermo::getHeatFlux(){ + Snapshot* currSnapshot = info_->getSnapshotManager()->getCurrentSnapshot(); + SimInfo::MoleculeIterator miter; + vector::iterator iiter; + Molecule* mol; + StuntDouble* sd; + RigidBody::AtomIterator ai; + Atom* atom; + Vector3d vel; + Vector3d angMom; + Mat3x3d I; + int i; + int j; + int k; + RealType mass; -double Thermo::getPressureZ() { + Vector3d x_a; + RealType kinetic; + RealType potential; + RealType eatom; + // Convective portion of the heat flux + Vector3d heatFluxJc = V3Zero; - // Relies on the calculation of the full molecular pressure tensor - - const double p_convert = 1.63882576e8; - double press[3][3]; - double pressureZ; + /* Calculate convective portion of the heat flux */ + for (mol = info_->beginMolecule(miter); mol != NULL; + mol = info_->nextMolecule(miter)) { + + for (sd = mol->beginIntegrableObject(iiter); + sd != NULL; + sd = mol->nextIntegrableObject(iiter)) { + + mass = sd->getMass(); + vel = sd->getVel(); - this->getPressureTensor(press); + kinetic = mass * (vel[0]*vel[0] + vel[1]*vel[1] + vel[2]*vel[2]); + + if (sd->isDirectional()) { + angMom = sd->getJ(); + I = sd->getI(); - pressureZ = p_convert * press[2][2]; + if (sd->isLinear()) { + i = sd->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); + } + } - return pressureZ; -} + potential = 0.0; + if (sd->isRigidBody()) { + RigidBody* rb = dynamic_cast(sd); + for (atom = rb->beginAtom(ai); atom != NULL; + atom = rb->nextAtom(ai)) { + potential += atom->getParticlePot(); + } + } else { + potential = sd->getParticlePot(); + } -void Thermo::getPressureTensor(double press[3][3]){ - // returns pressure tensor in units amu*fs^-2*Ang^-1 - // routine derived via viral theorem description in: - // Paci, E. and Marchi, M. J.Phys.Chem. 1996, 100, 4314-4322 + 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 + } + } - const double e_convert = 4.184e-4; - - double molmass, volume; - double vcom[3]; - double p_local[9], p_global[9]; - int i, j, k; - - for (i=0; i < 9; i++) { - p_local[i] = 0.0; - p_global[i] = 0.0; - } - - // use velocities of integrableObjects and their masses: - - for (i=0; i < info->integrableObjects.size(); i++) { - - molmass = info->integrableObjects[i]->getMass(); + /* The J_v vector is reduced in the forceManager 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 - info->integrableObjects[i]->getVel(vcom); - - p_local[0] += molmass * (vcom[0] * vcom[0]); - p_local[1] += molmass * (vcom[0] * vcom[1]); - p_local[2] += molmass * (vcom[0] * vcom[2]); - p_local[3] += molmass * (vcom[1] * vcom[0]); - p_local[4] += molmass * (vcom[1] * vcom[1]); - p_local[5] += molmass * (vcom[1] * vcom[2]); - p_local[6] += molmass * (vcom[2] * vcom[0]); - p_local[7] += molmass * (vcom[2] * vcom[1]); - p_local[8] += molmass * (vcom[2] * vcom[2]); + // (kcal/mol * A/fs) * conversion => (amu A^3)/fs^3 + Vector3d heatFluxJv = currSnapshot->getConductiveHeatFlux() * + PhysicalConstants::energyConvert; + + // Correct for the fact the flux is 1/V (Jc + Jv) + return (heatFluxJv + heatFluxJc) / this->getVolume(); // amu / fs^3 } - // Get total for entire system from MPI. - -#ifdef IS_MPI - MPI_Allreduce(p_local,p_global,9,MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD); -#else - for (i=0; i<9; i++) { - p_global[i] = p_local[i]; - } -#endif // is_mpi - volume = this->getVolume(); + Vector3d Thermo::getComVel(){ + Snapshot* snap = info_->getSnapshotManager()->getCurrentSnapshot(); + if (!snap->hasCOMvel) { - - for(i = 0; i < 3; i++) { - for (j = 0; j < 3; j++) { - k = 3*i + j; - press[i][j] = (p_global[k] + info->tau[k]*e_convert) / volume; + SimInfo::MoleculeIterator i; + Molecule* mol; + + Vector3d comVel(0.0); + RealType totalMass(0.0); + + for (mol = info_->beginMolecule(i); mol != NULL; + mol = info_->nextMolecule(i)) { + RealType mass = mol->getMass(); + totalMass += mass; + comVel += mass * mol->getComVel(); + } + +#ifdef IS_MPI + MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &totalMass, 1, MPI::REALTYPE, + MPI::SUM); + MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, comVel.getArrayPointer(), 3, + MPI::REALTYPE, MPI::SUM); +#endif + + comVel /= totalMass; + snap->setCOMvel(comVel); } + return snap->getCOMvel(); } -} -void Thermo::velocitize() { - - double aVel[3], aJ[3], I[3][3]; - int i, j, l, m, n, vr, vd; // velocity randomizer loop counters - double vdrift[3]; - double vbar; - const double kb = 8.31451e-7; // kb in amu, angstroms, fs, etc. - double av2; - double kebar; - double temperature; - int nobj; + Vector3d Thermo::getCom(){ + Snapshot* snap = info_->getSnapshotManager()->getCurrentSnapshot(); - if (!info->have_target_temp) { - sprintf( painCave.errMsg, - "You can't resample the velocities without a targetTemp!\n" - ); - painCave.isFatal = 1; - painCave.severity = OOPSE_ERROR; - simError(); - return; - } + if (!snap->hasCOM) { + + SimInfo::MoleculeIterator i; + Molecule* mol; + + Vector3d com(0.0); + RealType totalMass(0.0); + + for (mol = info_->beginMolecule(i); mol != NULL; + mol = info_->nextMolecule(i)) { + RealType mass = mol->getMass(); + totalMass += mass; + com += mass * mol->getCom(); + } + +#ifdef IS_MPI + MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &totalMass, 1, MPI::REALTYPE, + MPI::SUM); + MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, com.getArrayPointer(), 3, + MPI::REALTYPE, MPI::SUM); +#endif + + com /= totalMass; + snap->setCOM(com); + } + return snap->getCOM(); + } - nobj = info->integrableObjects.size(); - - temperature = info->target_temp; - - kebar = kb * temperature * (double)info->ndfRaw / - ( 2.0 * (double)info->ndf ); - - for(vr = 0; vr < nobj; vr++){ - - // uses equipartition theory to solve for vbar in angstrom/fs + /** + * Returns center of mass and center of mass velocity in one + * function call. + */ + void Thermo::getComAll(Vector3d &com, Vector3d &comVel){ + Snapshot* snap = info_->getSnapshotManager()->getCurrentSnapshot(); - av2 = 2.0 * kebar / info->integrableObjects[vr]->getMass(); - vbar = sqrt( av2 ); + if (!(snap->hasCOM && snap->hasCOMvel)) { - // picks random velocities from a gaussian distribution - // centered on vbar + SimInfo::MoleculeIterator i; + Molecule* mol; + + RealType totalMass(0.0); + + com = 0.0; + comVel = 0.0; + + for (mol = info_->beginMolecule(i); mol != NULL; + mol = info_->nextMolecule(i)) { + RealType mass = mol->getMass(); + totalMass += mass; + com += mass * mol->getCom(); + comVel += mass * mol->getComVel(); + } + +#ifdef IS_MPI + MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &totalMass, 1, MPI::REALTYPE, + MPI::SUM); + MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, com.getArrayPointer(), 3, + MPI::REALTYPE, MPI::SUM); + MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, comVel.getArrayPointer(), 3, + MPI::REALTYPE, MPI::SUM); +#endif + + com /= totalMass; + comVel /= totalMass; + snap->setCOM(com); + snap->setCOMvel(comVel); + } + com = snap->getCOM(); + comVel = snap->getCOMvel(); + return; + } + + /** + * \brief Return inertia tensor for entire system and angular momentum + * Vector. + * + * + * + * [ Ixx -Ixy -Ixz ] + * I =| -Iyx Iyy -Iyz | + * [ -Izx -Iyz Izz ] + */ + void Thermo::getInertiaTensor(Mat3x3d &inertiaTensor, + Vector3d &angularMomentum){ - for (j=0; j<3; j++) - aVel[j] = vbar * gaussStream->getGaussian(); + Snapshot* snap = info_->getSnapshotManager()->getCurrentSnapshot(); - info->integrableObjects[vr]->setVel( aVel ); + if (!(snap->hasInertiaTensor && snap->hasCOMw)) { + + RealType xx = 0.0; + RealType yy = 0.0; + RealType zz = 0.0; + RealType xy = 0.0; + RealType xz = 0.0; + RealType yz = 0.0; + Vector3d com(0.0); + Vector3d comVel(0.0); + + getComAll(com, comVel); + + SimInfo::MoleculeIterator i; + Molecule* mol; + + Vector3d thisq(0.0); + Vector3d thisv(0.0); + + RealType thisMass = 0.0; + + for (mol = info_->beginMolecule(i); mol != NULL; + mol = info_->nextMolecule(i)) { + + thisq = mol->getCom()-com; + thisv = mol->getComVel()-comVel; + thisMass = mol->getMass(); + // Compute moment of intertia coefficients. + xx += thisq[0]*thisq[0]*thisMass; + yy += thisq[1]*thisq[1]*thisMass; + zz += thisq[2]*thisq[2]*thisMass; + + // compute products of intertia + xy += thisq[0]*thisq[1]*thisMass; + xz += thisq[0]*thisq[2]*thisMass; + yz += thisq[1]*thisq[2]*thisMass; + + angularMomentum += cross( thisq, thisv ) * thisMass; + } + + inertiaTensor(0,0) = yy + zz; + inertiaTensor(0,1) = -xy; + inertiaTensor(0,2) = -xz; + inertiaTensor(1,0) = -xy; + inertiaTensor(1,1) = xx + zz; + inertiaTensor(1,2) = -yz; + inertiaTensor(2,0) = -xz; + inertiaTensor(2,1) = -yz; + inertiaTensor(2,2) = xx + yy; + +#ifdef IS_MPI + MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, inertiaTensor.getArrayPointer(), + 9, MPI::REALTYPE, MPI::SUM); + MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, + angularMomentum.getArrayPointer(), 3, + MPI::REALTYPE, MPI::SUM); +#endif + + snap->setCOMw(angularMomentum); + snap->setInertiaTensor(inertiaTensor); + } - if(info->integrableObjects[vr]->isDirectional()){ + angularMomentum = snap->getCOMw(); + inertiaTensor = snap->getInertiaTensor(); + + return; + } - info->integrableObjects[vr]->getI( I ); - if (info->integrableObjects[vr]->isLinear()) { - - l= info->integrableObjects[vr]->linearAxis(); - m = (l+1)%3; - n = (l+2)%3; - - aJ[l] = 0.0; - vbar = sqrt( 2.0 * kebar * I[m][m] ); - aJ[m] = vbar * gaussStream->getGaussian(); - vbar = sqrt( 2.0 * kebar * I[n][n] ); - aJ[n] = vbar * gaussStream->getGaussian(); + Mat3x3d Thermo::getBoundingBox(){ + + Snapshot* snap = info_->getSnapshotManager()->getCurrentSnapshot(); + + if (!(snap->hasBoundingBox)) { + + SimInfo::MoleculeIterator i; + Molecule::RigidBodyIterator ri; + Molecule::AtomIterator ai; + Molecule* mol; + RigidBody* rb; + Atom* atom; + Vector3d pos, bMax, bMin; + int index = 0; + + for (mol = info_->beginMolecule(i); mol != NULL; + mol = info_->nextMolecule(i)) { - } else { - for (j = 0 ; j < 3; j++) { - vbar = sqrt( 2.0 * kebar * I[j][j] ); - aJ[j] = vbar * gaussStream->getGaussian(); - } - } // else isLinear + //change the positions of atoms which belong to the rigidbodies + for (rb = mol->beginRigidBody(ri); rb != NULL; + rb = mol->nextRigidBody(ri)) { + rb->updateAtoms(); + } + + for(atom = mol->beginAtom(ai); atom != NULL; + atom = mol->nextAtom(ai)) { + + pos = atom->getPos(); - info->integrableObjects[vr]->setJ( aJ ); + if (index == 0) { + bMax = pos; + bMin = pos; + } else { + for (int i = 0; i < 3; i++) { + bMax[i] = max(bMax[i], pos[i]); + bMin[i] = min(bMin[i], pos[i]); + } + } + index++; + } + } - }//isDirectional +#ifdef IS_MPI + MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &bMax[0], 3, MPI::REALTYPE, + MPI::MAX); + MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &bMin[0], 3, MPI::REALTYPE, + MPI::MIN); +#endif + Mat3x3d bBox = Mat3x3d(0.0); + for (int i = 0; i < 3; i++) { + bBox(i,i) = bMax[i] - bMin[i]; + } + snap->setBoundingBox(bBox); + } + + return snap->getBoundingBox(); } - - // Get the Center of Mass drift velocity. - - getCOMVel(vdrift); - // Corrects for the center of mass drift. - // sums all the momentum and divides by total mass. - - for(vd = 0; vd < nobj; vd++){ + + // Returns the angular momentum of the system + Vector3d Thermo::getAngularMomentum(){ + Snapshot* snap = info_->getSnapshotManager()->getCurrentSnapshot(); - info->integrableObjects[vd]->getVel(aVel); - - for (j=0; j < 3; j++) - aVel[j] -= vdrift[j]; + if (!snap->hasCOMw) { + + Vector3d com(0.0); + Vector3d comVel(0.0); + Vector3d angularMomentum(0.0); + + getComAll(com, comVel); + + SimInfo::MoleculeIterator i; + Molecule* mol; + + Vector3d thisr(0.0); + Vector3d thisp(0.0); + + RealType thisMass; + + for (mol = info_->beginMolecule(i); mol != NULL; + mol = info_->nextMolecule(i)) { + thisMass = mol->getMass(); + thisr = mol->getCom() - com; + thisp = (mol->getComVel() - comVel) * thisMass; - info->integrableObjects[vd]->setVel( aVel ); + angularMomentum += cross( thisr, thisp ); + } + +#ifdef IS_MPI + MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, + angularMomentum.getArrayPointer(), 3, + MPI::REALTYPE, MPI::SUM); +#endif + + snap->setCOMw(angularMomentum); + } + + return snap->getCOMw(); } - -} - -void Thermo::getCOMVel(double vdrift[3]){ - - double mtot, mtot_local; - double aVel[3], amass; - double vdrift_local[3]; - int vd, j; - int nobj; - - nobj = info->integrableObjects.size(); - - mtot_local = 0.0; - vdrift_local[0] = 0.0; - vdrift_local[1] = 0.0; - vdrift_local[2] = 0.0; - for(vd = 0; vd < nobj; vd++){ + + /** + * Returns the Volume of the system based on a ellipsoid with + * semi-axes based on the radius of gyration V=4/3*Pi*R_1*R_2*R_3 + * where R_i are related to the principle inertia moments + * R_i = sqrt(C*I_i/N), this reduces to + * V = 4/3*Pi*(C/N)^3/2*sqrt(det(I)). + * See S.E. Baltazar et. al. Comp. Mat. Sci. 37 (2006) 526-536. + */ + RealType Thermo::getGyrationalVolume(){ + Snapshot* snap = info_->getSnapshotManager()->getCurrentSnapshot(); - amass = info->integrableObjects[vd]->getMass(); - info->integrableObjects[vd]->getVel( aVel ); + if (!snap->hasGyrationalVolume) { + + Mat3x3d intTensor; + RealType det; + Vector3d dummyAngMom; + RealType sysconstants; + RealType geomCnst; + RealType volume; + + geomCnst = 3.0/2.0; + /* Get the inertial tensor and angular momentum for free*/ + getInertiaTensor(intTensor, dummyAngMom); + + det = intTensor.determinant(); + sysconstants = geomCnst / (RealType)(info_->getNGlobalIntegrableObjects()); + volume = 4.0/3.0*NumericConstant::PI*pow(sysconstants,geomCnst)*sqrt(det); - for(j = 0; j < 3; j++) - vdrift_local[j] += aVel[j] * amass; - - mtot_local += amass; + snap->setGyrationalVolume(volume); + } + return snap->getGyrationalVolume(); } + + void Thermo::getGyrationalVolume(RealType &volume, RealType &detI){ + Snapshot* snap = info_->getSnapshotManager()->getCurrentSnapshot(); -#ifdef IS_MPI - MPI_Allreduce(&mtot_local,&mtot,1,MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD); - MPI_Allreduce(vdrift_local,vdrift,3,MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD); -#else - mtot = mtot_local; - for(vd = 0; vd < 3; vd++) { - vdrift[vd] = vdrift_local[vd]; - } -#endif + if (!(snap->hasInertiaTensor && snap->hasGyrationalVolume)) { - for (vd = 0; vd < 3; vd++) { - vdrift[vd] = vdrift[vd] / mtot; + Mat3x3d intTensor; + Vector3d dummyAngMom; + RealType sysconstants; + RealType geomCnst; + + geomCnst = 3.0/2.0; + /* Get the inertia tensor and angular momentum for free*/ + this->getInertiaTensor(intTensor, dummyAngMom); + + detI = intTensor.determinant(); + sysconstants = geomCnst/(RealType)(info_->getNGlobalIntegrableObjects()); + volume = 4.0/3.0*NumericConstant::PI*pow(sysconstants,geomCnst)*sqrt(detI); + snap->setGyrationalVolume(volume); + } else { + volume = snap->getGyrationalVolume(); + detI = snap->getInertiaTensor().determinant(); + } + return; } -} + RealType Thermo::getTaggedAtomPairDistance(){ + Snapshot* currSnapshot = info_->getSnapshotManager()->getCurrentSnapshot(); + Globals* simParams = info_->getSimParams(); + + if (simParams->haveTaggedAtomPair() && + simParams->havePrintTaggedPairDistance()) { + if ( simParams->getPrintTaggedPairDistance()) { + + pair tap = simParams->getTaggedAtomPair(); + Vector3d pos1, pos2, rab; + +#ifdef IS_MPI + int mol1 = info_->getGlobalMolMembership(tap.first); + int mol2 = info_->getGlobalMolMembership(tap.second); -void Thermo::getCOM(double COM[3]){ + int proc1 = info_->getMolToProc(mol1); + int proc2 = info_->getMolToProc(mol2); - double mtot, mtot_local; - double aPos[3], amass; - double COM_local[3]; - int i, j; - int nobj; + RealType data[3]; + if (proc1 == worldRank) { + StuntDouble* sd1 = info_->getIOIndexToIntegrableObject(tap.first); + pos1 = sd1->getPos(); + data[0] = pos1.x(); + data[1] = pos1.y(); + data[2] = pos1.z(); + MPI::COMM_WORLD.Bcast(data, 3, MPI::REALTYPE, proc1); + } else { + MPI::COMM_WORLD.Bcast(data, 3, MPI::REALTYPE, proc1); + pos1 = Vector3d(data); + } - mtot_local = 0.0; - COM_local[0] = 0.0; - COM_local[1] = 0.0; - COM_local[2] = 0.0; - - nobj = info->integrableObjects.size(); - for(i = 0; i < nobj; i++){ - - amass = info->integrableObjects[i]->getMass(); - info->integrableObjects[i]->getPos( aPos ); - - for(j = 0; j < 3; j++) - COM_local[j] += aPos[j] * amass; - - mtot_local += amass; + if (proc2 == worldRank) { + StuntDouble* sd2 = info_->getIOIndexToIntegrableObject(tap.second); + pos2 = sd2->getPos(); + data[0] = pos2.x(); + data[1] = pos2.y(); + data[2] = pos2.z(); + MPI::COMM_WORLD.Bcast(data, 3, MPI::REALTYPE, proc2); + } else { + MPI::COMM_WORLD.Bcast(data, 3, MPI::REALTYPE, proc2); + 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); + return rab.length(); + } + return 0.0; + } + return 0.0; } -#ifdef IS_MPI - MPI_Allreduce(&mtot_local,&mtot,1,MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD); - MPI_Allreduce(COM_local,COM,3,MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD); + RealType Thermo::getHullVolume(){ +#ifdef HAVE_QHULL + Snapshot* snap = info_->getSnapshotManager()->getCurrentSnapshot(); + if (!snap->hasHullVolume) { + Hull* surfaceMesh_; + + Globals* simParams = info_->getSimParams(); + const std::string ht = simParams->getHULL_Method(); + + if (ht == "Convex") { + surfaceMesh_ = new ConvexHull(); + } else if (ht == "AlphaShape") { + surfaceMesh_ = new AlphaHull(simParams->getAlpha()); + } else { + return 0.0; + } + + // Build a vector of stunt doubles to determine if they are + // surface atoms + std::vector localSites_; + Molecule* mol; + StuntDouble* sd; + SimInfo::MoleculeIterator i; + Molecule::IntegrableObjectIterator j; + + for (mol = info_->beginMolecule(i); mol != NULL; + mol = info_->nextMolecule(i)) { + for (sd = mol->beginIntegrableObject(j); + sd != NULL; + sd = mol->nextIntegrableObject(j)) { + localSites_.push_back(sd); + } + } + + // Compute surface Mesh + surfaceMesh_->computeHull(localSites_); + snap->setHullVolume(surfaceMesh_->getVolume()); + + delete surfaceMesh_; + } + + return snap->getHullVolume(); #else - mtot = mtot_local; - for(i = 0; i < 3; i++) { - COM[i] = COM_local[i]; - } + return 0.0; #endif - - for (i = 0; i < 3; i++) { - COM[i] = COM[i] / mtot; } -} -void Thermo::removeCOMdrift() { - double vdrift[3], aVel[3]; - int vd, j, nobj; - nobj = info->integrableObjects.size(); - - // Get the Center of Mass drift velocity. - - getCOMVel(vdrift); - - // Corrects for the center of mass drift. - // sums all the momentum and divides by total mass. - - for(vd = 0; vd < nobj; vd++){ - - info->integrableObjects[vd]->getVel(aVel); - - for (j=0; j < 3; j++) - aVel[j] -= vdrift[j]; - - info->integrableObjects[vd]->setVel( aVel ); - } }