--- trunk/src/brains/Thermo.cpp 2012/08/31 17:29:35 1792 +++ trunk/src/brains/Thermo.cpp 2014/09/26 22:22:28 2022 @@ -35,17 +35,17 @@ * * [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). + * [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 #ifdef IS_MPI #include #endif //is_mpi + +#include +#include #include "brains/Thermo.hpp" #include "primitives/Molecule.hpp" @@ -89,8 +89,8 @@ namespace OpenMD { } #ifdef IS_MPI - MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &kinetic, 1, MPI::REALTYPE, - MPI::SUM); + MPI_Allreduce(MPI_IN_PLACE, &kinetic, 1, MPI_REALTYPE, + MPI_SUM, MPI_COMM_WORLD); #endif kinetic = kinetic * 0.5 / PhysicalConstants::energyConvert; @@ -140,8 +140,8 @@ namespace OpenMD { } #ifdef IS_MPI - MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &kinetic, 1, MPI::REALTYPE, - MPI::SUM); + MPI_Allreduce(MPI_IN_PLACE, &kinetic, 1, MPI_REALTYPE, + MPI_SUM, MPI_COMM_WORLD); #endif kinetic = kinetic * 0.5 / PhysicalConstants::energyConvert; @@ -227,13 +227,13 @@ namespace OpenMD { } #ifdef IS_MPI - MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &kinetic, 1, MPI::REALTYPE, - MPI::SUM); + MPI_Allreduce(MPI_IN_PLACE, &kinetic, 1, MPI_REALTYPE, + MPI_SUM, MPI_COMM_WORLD); #endif kinetic *= 0.5; eTemp = (2.0 * kinetic) / - (info_->getNFluctuatingCharges() * PhysicalConstants::kb ); + (info_->getNFluctuatingCharges() * PhysicalConstants::kb ); snap->setElectronicTemperature(eTemp); } @@ -297,8 +297,8 @@ namespace OpenMD { } #ifdef IS_MPI - MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, p_tens.getArrayPointer(), 9, - MPI::REALTYPE, MPI::SUM); + MPI_Allreduce(MPI_IN_PLACE, p_tens.getArrayPointer(), 9, + MPI_REALTYPE, MPI_SUM, MPI_COMM_WORLD); #endif RealType volume = this->getVolume(); @@ -324,7 +324,6 @@ namespace OpenMD { Molecule* mol; Atom* atom; RealType charge; - RealType moment(0.0); Vector3d ri(0.0); Vector3d dipoleVector(0.0); Vector3d nPos(0.0); @@ -372,35 +371,31 @@ namespace OpenMD { 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; + 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); + MPI_Allreduce(MPI_IN_PLACE, &pChg, 1, MPI_REALTYPE, + MPI_SUM, MPI_COMM_WORLD); + MPI_Allreduce(MPI_IN_PLACE, &nChg, 1, MPI_REALTYPE, + MPI_SUM, MPI_COMM_WORLD); + + MPI_Allreduce(MPI_IN_PLACE, &pCount, 1, MPI_INTEGER, + MPI_SUM, MPI_COMM_WORLD); + MPI_Allreduce(MPI_IN_PLACE, &nCount, 1, MPI_INTEGER, + MPI_SUM, MPI_COMM_WORLD); + + MPI_Allreduce(MPI_IN_PLACE, pPos.getArrayPointer(), 3, + MPI_REALTYPE, MPI_SUM, MPI_COMM_WORLD); + MPI_Allreduce(MPI_IN_PLACE, nPos.getArrayPointer(), 3, + MPI_REALTYPE, MPI_SUM, MPI_COMM_WORLD); - 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); - - 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); - - MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, dipoleVector.getArrayPointer(), - 3, MPI::REALTYPE, MPI::SUM); + MPI_Allreduce(MPI_IN_PLACE, dipoleVector.getArrayPointer(), + 3, MPI_REALTYPE, MPI_SUM, MPI_COMM_WORLD); #endif // first load the accumulated dipole moment (if dipoles were present) @@ -421,6 +416,74 @@ namespace OpenMD { } return snap->getSystemDipole(); + } + + + Mat3x3d Thermo::getSystemQuadrupole() { + Snapshot* snap = info_->getSnapshotManager()->getCurrentSnapshot(); + + if (!snap->hasSystemQuadrupole) { + SimInfo::MoleculeIterator miter; + vector::iterator aiter; + Molecule* mol; + Atom* atom; + RealType charge; + Vector3d ri(0.0); + Vector3d dipole(0.0); + Mat3x3d qpole(0.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)) { + + ri = atom->getPos(); + snap->wrapVector(ri); + ri *= angstromToM; + + 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; + + qpole += 0.5 * charge * outProduct(ri, ri); + + MultipoleAdapter ma = MultipoleAdapter(atom->getAtomType()); + + if ( ma.isDipole() ) { + dipole = atom->getDipole() * debyeToCm; + qpole += 0.5 * outProduct( dipole, ri ); + } + + if ( ma.isQuadrupole() ) { + qpole += atom->getQuadrupole() * debyeToCm * angstromToM; + } + } + } + +#ifdef IS_MPI + MPI_Allreduce(MPI_IN_PLACE, qpole.getArrayPointer(), + 9, MPI_REALTYPE, MPI_SUM, MPI_COMM_WORLD); +#endif + + snap->setSystemQuadrupole(qpole); + } + + return snap->getSystemQuadrupole(); } // Returns the Heat Flux Vector for the system @@ -503,8 +566,8 @@ namespace OpenMD { * reduced among all processors. */ #ifdef IS_MPI - MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, &heatFluxJc[0], 3, MPI::REALTYPE, - MPI::SUM); + MPI_Allreduce(MPI_IN_PLACE, &heatFluxJc[0], 3, MPI_REALTYPE, + MPI_SUM, MPI_COMM_WORLD); #endif // (kcal/mol * A/fs) * conversion => (amu A^3)/fs^3 @@ -536,10 +599,10 @@ namespace OpenMD { } #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); + MPI_Allreduce(MPI_IN_PLACE, &totalMass, 1, MPI_REALTYPE, + MPI_SUM, MPI_COMM_WORLD); + MPI_Allreduce(MPI_IN_PLACE, comVel.getArrayPointer(), 3, + MPI_REALTYPE, MPI_SUM, MPI_COMM_WORLD); #endif comVel /= totalMass; @@ -567,10 +630,10 @@ namespace OpenMD { } #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_Allreduce(MPI_IN_PLACE, &totalMass, 1, MPI_REALTYPE, + MPI_SUM, MPI_COMM_WORLD); + MPI_Allreduce(MPI_IN_PLACE, com.getArrayPointer(), 3, + MPI_REALTYPE, MPI_SUM, MPI_COMM_WORLD); #endif com /= totalMass; @@ -605,12 +668,12 @@ namespace OpenMD { } #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); + MPI_Allreduce(MPI_IN_PLACE, &totalMass, 1, MPI_REALTYPE, + MPI_SUM, MPI_COMM_WORLD); + MPI_Allreduce(MPI_IN_PLACE, com.getArrayPointer(), 3, + MPI_REALTYPE, MPI_SUM, MPI_COMM_WORLD); + MPI_Allreduce(MPI_IN_PLACE, comVel.getArrayPointer(), 3, + MPI_REALTYPE, MPI_SUM, MPI_COMM_WORLD); #endif com /= totalMass; @@ -624,8 +687,8 @@ namespace OpenMD { } /** - * Return intertia tensor for entire system and angular momentum - * Vector. + * \brief Return inertia tensor for entire system and angular momentum + * Vector. * * * @@ -689,11 +752,11 @@ namespace OpenMD { 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); + MPI_Allreduce(MPI_IN_PLACE, inertiaTensor.getArrayPointer(), + 9, MPI_REALTYPE, MPI_SUM, MPI_COMM_WORLD); + MPI_Allreduce(MPI_IN_PLACE, + angularMomentum.getArrayPointer(), 3, + MPI_REALTYPE, MPI_SUM, MPI_COMM_WORLD); #endif snap->setCOMw(angularMomentum); @@ -706,6 +769,67 @@ namespace OpenMD { return; } + + 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)) { + + //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(); + + 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++; + } + } + +#ifdef IS_MPI + MPI_Allreduce(MPI_IN_PLACE, &bMax[0], 3, MPI_REALTYPE, + MPI_MAX, MPI_COMM_WORLD); + + MPI_Allreduce(MPI_IN_PLACE, &bMin[0], 3, MPI_REALTYPE, + MPI_MIN, MPI_COMM_WORLD); +#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(); + } + + // Returns the angular momentum of the system Vector3d Thermo::getAngularMomentum(){ Snapshot* snap = info_->getSnapshotManager()->getCurrentSnapshot(); @@ -736,9 +860,9 @@ namespace OpenMD { } #ifdef IS_MPI - MPI::COMM_WORLD.Allreduce(MPI::IN_PLACE, - angularMomentum.getArrayPointer(), 3, - MPI::REALTYPE, MPI::SUM); + MPI_Allreduce(MPI_IN_PLACE, + angularMomentum.getArrayPointer(), 3, + MPI_REALTYPE, MPI_SUM, MPI_COMM_WORLD); #endif snap->setCOMw(angularMomentum); @@ -842,7 +966,7 @@ namespace OpenMD { pos2 = sd2->getPos(); data[0] = pos2.x(); data[1] = pos2.y(); - data[2] = pos2.z(); + 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); @@ -864,12 +988,11 @@ namespace OpenMD { } RealType Thermo::getHullVolume(){ - Snapshot* snap = info_->getSnapshotManager()->getCurrentSnapshot(); - #ifdef HAVE_QHULL + Snapshot* snap = info_->getSnapshotManager()->getCurrentSnapshot(); if (!snap->hasHullVolume) { Hull* surfaceMesh_; - + Globals* simParams = info_->getSimParams(); const std::string ht = simParams->getHULL_Method(); @@ -901,10 +1024,15 @@ namespace OpenMD { // Compute surface Mesh surfaceMesh_->computeHull(localSites_); snap->setHullVolume(surfaceMesh_->getVolume()); + + delete surfaceMesh_; } + return snap->getHullVolume(); #else return 0.0; #endif } + + }