--- trunk/src/brains/SimInfo.cpp 2005/09/16 21:07:45 608 +++ branches/development/src/brains/SimInfo.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). */ /** @@ -48,65 +48,81 @@ #include #include +#include #include "brains/SimInfo.hpp" #include "math/Vector3.hpp" #include "primitives/Molecule.hpp" +#include "primitives/StuntDouble.hpp" #include "UseTheForce/fCutoffPolicy.h" -#include "UseTheForce/DarkSide/fElectrostaticSummationMethod.h" +#include "UseTheForce/DarkSide/fSwitchingFunctionType.h" #include "UseTheForce/doForces_interface.h" -#include "UseTheForce/notifyCutoffs_interface.h" +#include "UseTheForce/DarkSide/neighborLists_interface.h" +#include "UseTheForce/DarkSide/switcheroo_interface.h" #include "utils/MemoryUtils.hpp" #include "utils/simError.h" #include "selection/SelectionManager.hpp" +#include "io/ForceFieldOptions.hpp" +#include "UseTheForce/ForceField.hpp" + #ifdef IS_MPI #include "UseTheForce/mpiComponentPlan.h" #include "UseTheForce/DarkSide/simParallel_interface.h" #endif -namespace oopse { +namespace OpenMD { + std::set getRigidSet(int index, std::map >& container) { + std::map >::iterator i = container.find(index); + std::set result; + if (i != container.end()) { + result = i->second; + } - SimInfo::SimInfo(MakeStamps* stamps, std::vector >& molStampPairs, - ForceField* ff, Globals* simParams) : - stamps_(stamps), forceField_(ff), simParams_(simParams), - ndf_(0), ndfRaw_(0), ndfTrans_(0), nZconstraint_(0), + return result; + } + + SimInfo::SimInfo(ForceField* ff, Globals* simParams) : + forceField_(ff), simParams_(simParams), + ndf_(0), fdf_local(0), ndfRaw_(0), ndfTrans_(0), nZconstraint_(0), nGlobalMols_(0), nGlobalAtoms_(0), nGlobalCutoffGroups_(0), nGlobalIntegrableObjects_(0), nGlobalRigidBodies_(0), - nAtoms_(0), nBonds_(0), nBends_(0), nTorsions_(0), nRigidBodies_(0), - nIntegrableObjects_(0), nCutoffGroups_(0), nConstraints_(0), - sman_(NULL), fortranInitialized_(false) { + nAtoms_(0), nBonds_(0), nBends_(0), nTorsions_(0), nInversions_(0), + nRigidBodies_(0), nIntegrableObjects_(0), nCutoffGroups_(0), + nConstraints_(0), sman_(NULL), fortranInitialized_(false), + calcBoxDipole_(false), useAtomicVirial_(true) { - - std::vector >::iterator i; + MoleculeStamp* molStamp; int nMolWithSameStamp; int nCutoffAtoms = 0; // number of atoms belong to cutoff groups - int nGroups = 0; //total cutoff groups defined in meta-data file + int nGroups = 0; //total cutoff groups defined in meta-data file CutoffGroupStamp* cgStamp; RigidBodyStamp* rbStamp; int nRigidAtoms = 0; - - for (i = molStampPairs.begin(); i !=molStampPairs.end(); ++i) { - molStamp = i->first; - nMolWithSameStamp = i->second; + + std::vector components = simParams->getComponents(); + + for (std::vector::iterator i = components.begin(); i !=components.end(); ++i) { + molStamp = (*i)->getMoleculeStamp(); + nMolWithSameStamp = (*i)->getNMol(); addMoleculeStamp(molStamp, nMolWithSameStamp); //calculate atoms in molecules nGlobalAtoms_ += molStamp->getNAtoms() *nMolWithSameStamp; - //calculate atoms in cutoff groups int nAtomsInGroups = 0; int nCutoffGroupsInStamp = molStamp->getNCutoffGroups(); for (int j=0; j < nCutoffGroupsInStamp; j++) { - cgStamp = molStamp->getCutoffGroup(j); + cgStamp = molStamp->getCutoffGroupStamp(j); nAtomsInGroups += cgStamp->getNMembers(); } nGroups += nCutoffGroupsInStamp * nMolWithSameStamp; + nCutoffAtoms += nAtomsInGroups * nMolWithSameStamp; //calculate atoms in rigid bodies @@ -114,7 +130,7 @@ namespace oopse { int nRigidBodiesInStamp = molStamp->getNRigidBodies(); for (int j=0; j < nRigidBodiesInStamp; j++) { - rbStamp = molStamp->getRigidBody(j); + rbStamp = molStamp->getRigidBodyStamp(j); nAtomsInRigidBodies += rbStamp->getNMembers(); } @@ -123,24 +139,23 @@ namespace oopse { } - //every free atom (atom does not belong to cutoff groups) is a cutoff group - //therefore the total number of cutoff groups in the system is equal to - //the total number of atoms minus number of atoms belong to cutoff group defined in meta-data - //file plus the number of cutoff groups defined in meta-data file + //every free atom (atom does not belong to cutoff groups) is a cutoff + //group therefore the total number of cutoff groups in the system is + //equal to the total number of atoms minus number of atoms belong to + //cutoff group defined in meta-data file plus the number of cutoff + //groups defined in meta-data file nGlobalCutoffGroups_ = nGlobalAtoms_ - nCutoffAtoms + nGroups; - //every free atom (atom does not belong to rigid bodies) is an integrable object - //therefore the total number of integrable objects in the system is equal to - //the total number of atoms minus number of atoms belong to rigid body defined in meta-data - //file plus the number of rigid bodies defined in meta-data file - nGlobalIntegrableObjects_ = nGlobalAtoms_ - nRigidAtoms + nGlobalRigidBodies_; - + //every free atom (atom does not belong to rigid bodies) is an + //integrable object therefore the total number of integrable objects + //in the system is equal to the total number of atoms minus number of + //atoms belong to rigid body defined in meta-data file plus the number + //of rigid bodies defined in meta-data file + nGlobalIntegrableObjects_ = nGlobalAtoms_ - nRigidAtoms + + nGlobalRigidBodies_; + nGlobalMols_ = molStampIds_.size(); - -#ifdef IS_MPI molToProcMap_.resize(nGlobalMols_); -#endif - } SimInfo::~SimInfo() { @@ -150,7 +165,6 @@ namespace oopse { } molecules_.clear(); - delete stamps_; delete sman_; delete simParams_; delete forceField_; @@ -179,13 +193,14 @@ namespace oopse { nBonds_ += mol->getNBonds(); nBends_ += mol->getNBends(); nTorsions_ += mol->getNTorsions(); + nInversions_ += mol->getNInversions(); nRigidBodies_ += mol->getNRigidBodies(); nIntegrableObjects_ += mol->getNIntegrableObjects(); nCutoffGroups_ += mol->getNCutoffGroups(); nConstraints_ += mol->getNConstraintPairs(); - addExcludePairs(mol); - + addInteractionPairs(mol); + return true; } else { return false; @@ -204,12 +219,13 @@ namespace oopse { nBonds_ -= mol->getNBonds(); nBends_ -= mol->getNBends(); nTorsions_ -= mol->getNTorsions(); + nInversions_ -= mol->getNInversions(); nRigidBodies_ -= mol->getNRigidBodies(); nIntegrableObjects_ -= mol->getNIntegrableObjects(); nCutoffGroups_ -= mol->getNCutoffGroups(); nConstraints_ -= mol->getNConstraintPairs(); - removeExcludePairs(mol); + removeInteractionPairs(mol); molecules_.erase(mol->getGlobalIndex()); delete mol; @@ -257,8 +273,8 @@ namespace oopse { } } - }//end for (integrableObject) - }// end for (mol) + } + } // n_constraints is local, so subtract them on each processor ndf_local -= nConstraints_; @@ -275,6 +291,15 @@ namespace oopse { } + int SimInfo::getFdf() { +#ifdef IS_MPI + MPI_Allreduce(&fdf_local,&fdf_,1,MPI_INT,MPI_SUM, MPI_COMM_WORLD); +#else + fdf_ = fdf_local; +#endif + return fdf_; + } + void SimInfo::calcNdfRaw() { int ndfRaw_local; @@ -326,124 +351,324 @@ namespace oopse { } - void SimInfo::addExcludePairs(Molecule* mol) { + void SimInfo::addInteractionPairs(Molecule* mol) { + ForceFieldOptions& options_ = forceField_->getForceFieldOptions(); std::vector::iterator bondIter; std::vector::iterator bendIter; std::vector::iterator torsionIter; + std::vector::iterator inversionIter; Bond* bond; Bend* bend; Torsion* torsion; + Inversion* inversion; int a; int b; int c; int d; - - for (bond= mol->beginBond(bondIter); bond != NULL; bond = mol->nextBond(bondIter)) { - a = bond->getAtomA()->getGlobalIndex(); - b = bond->getAtomB()->getGlobalIndex(); - exclude_.addPair(a, b); + + // atomGroups can be used to add special interaction maps between + // groups of atoms that are in two separate rigid bodies. + // However, most site-site interactions between two rigid bodies + // are probably not special, just the ones between the physically + // bonded atoms. Interactions *within* a single rigid body should + // always be excluded. These are done at the bottom of this + // function. + + std::map > atomGroups; + Molecule::RigidBodyIterator rbIter; + RigidBody* rb; + Molecule::IntegrableObjectIterator ii; + StuntDouble* integrableObject; + + for (integrableObject = mol->beginIntegrableObject(ii); + integrableObject != NULL; + integrableObject = mol->nextIntegrableObject(ii)) { + + if (integrableObject->isRigidBody()) { + rb = static_cast(integrableObject); + std::vector atoms = rb->getAtoms(); + std::set rigidAtoms; + for (int i = 0; i < static_cast(atoms.size()); ++i) { + rigidAtoms.insert(atoms[i]->getGlobalIndex()); + } + for (int i = 0; i < static_cast(atoms.size()); ++i) { + atomGroups.insert(std::map >::value_type(atoms[i]->getGlobalIndex(), rigidAtoms)); + } + } else { + std::set oneAtomSet; + oneAtomSet.insert(integrableObject->getGlobalIndex()); + atomGroups.insert(std::map >::value_type(integrableObject->getGlobalIndex(), oneAtomSet)); + } + } + + for (bond= mol->beginBond(bondIter); bond != NULL; + bond = mol->nextBond(bondIter)) { + + a = bond->getAtomA()->getGlobalIndex(); + b = bond->getAtomB()->getGlobalIndex(); + + if (options_.havevdw12scale() || options_.haveelectrostatic12scale()) { + oneTwoInteractions_.addPair(a, b); + } else { + excludedInteractions_.addPair(a, b); + } } - for (bend= mol->beginBend(bendIter); bend != NULL; bend = mol->nextBend(bendIter)) { + for (bend= mol->beginBend(bendIter); bend != NULL; + bend = mol->nextBend(bendIter)) { + a = bend->getAtomA()->getGlobalIndex(); b = bend->getAtomB()->getGlobalIndex(); c = bend->getAtomC()->getGlobalIndex(); + + if (options_.havevdw12scale() || options_.haveelectrostatic12scale()) { + oneTwoInteractions_.addPair(a, b); + oneTwoInteractions_.addPair(b, c); + } else { + excludedInteractions_.addPair(a, b); + excludedInteractions_.addPair(b, c); + } - exclude_.addPair(a, b); - exclude_.addPair(a, c); - exclude_.addPair(b, c); + if (options_.havevdw13scale() || options_.haveelectrostatic13scale()) { + oneThreeInteractions_.addPair(a, c); + } else { + excludedInteractions_.addPair(a, c); + } } - for (torsion= mol->beginTorsion(torsionIter); torsion != NULL; torsion = mol->nextTorsion(torsionIter)) { + for (torsion= mol->beginTorsion(torsionIter); torsion != NULL; + torsion = mol->nextTorsion(torsionIter)) { + a = torsion->getAtomA()->getGlobalIndex(); b = torsion->getAtomB()->getGlobalIndex(); c = torsion->getAtomC()->getGlobalIndex(); - d = torsion->getAtomD()->getGlobalIndex(); + d = torsion->getAtomD()->getGlobalIndex(); - exclude_.addPair(a, b); - exclude_.addPair(a, c); - exclude_.addPair(a, d); - exclude_.addPair(b, c); - exclude_.addPair(b, d); - exclude_.addPair(c, d); + if (options_.havevdw12scale() || options_.haveelectrostatic12scale()) { + oneTwoInteractions_.addPair(a, b); + oneTwoInteractions_.addPair(b, c); + oneTwoInteractions_.addPair(c, d); + } else { + excludedInteractions_.addPair(a, b); + excludedInteractions_.addPair(b, c); + excludedInteractions_.addPair(c, d); + } + + if (options_.havevdw13scale() || options_.haveelectrostatic13scale()) { + oneThreeInteractions_.addPair(a, c); + oneThreeInteractions_.addPair(b, d); + } else { + excludedInteractions_.addPair(a, c); + excludedInteractions_.addPair(b, d); + } + + if (options_.havevdw14scale() || options_.haveelectrostatic14scale()) { + oneFourInteractions_.addPair(a, d); + } else { + excludedInteractions_.addPair(a, d); + } } - Molecule::RigidBodyIterator rbIter; - RigidBody* rb; - for (rb = mol->beginRigidBody(rbIter); rb != NULL; rb = mol->nextRigidBody(rbIter)) { + for (inversion= mol->beginInversion(inversionIter); inversion != NULL; + inversion = mol->nextInversion(inversionIter)) { + + a = inversion->getAtomA()->getGlobalIndex(); + b = inversion->getAtomB()->getGlobalIndex(); + c = inversion->getAtomC()->getGlobalIndex(); + d = inversion->getAtomD()->getGlobalIndex(); + + if (options_.havevdw12scale() || options_.haveelectrostatic12scale()) { + oneTwoInteractions_.addPair(a, b); + oneTwoInteractions_.addPair(a, c); + oneTwoInteractions_.addPair(a, d); + } else { + excludedInteractions_.addPair(a, b); + excludedInteractions_.addPair(a, c); + excludedInteractions_.addPair(a, d); + } + + if (options_.havevdw13scale() || options_.haveelectrostatic13scale()) { + oneThreeInteractions_.addPair(b, c); + oneThreeInteractions_.addPair(b, d); + oneThreeInteractions_.addPair(c, d); + } else { + excludedInteractions_.addPair(b, c); + excludedInteractions_.addPair(b, d); + excludedInteractions_.addPair(c, d); + } + } + + for (rb = mol->beginRigidBody(rbIter); rb != NULL; + rb = mol->nextRigidBody(rbIter)) { std::vector atoms = rb->getAtoms(); - for (int i = 0; i < atoms.size() -1 ; ++i) { - for (int j = i + 1; j < atoms.size(); ++j) { + for (int i = 0; i < static_cast(atoms.size()) -1 ; ++i) { + for (int j = i + 1; j < static_cast(atoms.size()); ++j) { a = atoms[i]->getGlobalIndex(); b = atoms[j]->getGlobalIndex(); - exclude_.addPair(a, b); + excludedInteractions_.addPair(a, b); } } } } - void SimInfo::removeExcludePairs(Molecule* mol) { + void SimInfo::removeInteractionPairs(Molecule* mol) { + ForceFieldOptions& options_ = forceField_->getForceFieldOptions(); std::vector::iterator bondIter; std::vector::iterator bendIter; std::vector::iterator torsionIter; + std::vector::iterator inversionIter; Bond* bond; Bend* bend; Torsion* torsion; + Inversion* inversion; int a; int b; int c; int d; + + std::map > atomGroups; + Molecule::RigidBodyIterator rbIter; + RigidBody* rb; + Molecule::IntegrableObjectIterator ii; + StuntDouble* integrableObject; - for (bond= mol->beginBond(bondIter); bond != NULL; bond = mol->nextBond(bondIter)) { + for (integrableObject = mol->beginIntegrableObject(ii); + integrableObject != NULL; + integrableObject = mol->nextIntegrableObject(ii)) { + + if (integrableObject->isRigidBody()) { + rb = static_cast(integrableObject); + std::vector atoms = rb->getAtoms(); + std::set rigidAtoms; + for (int i = 0; i < static_cast(atoms.size()); ++i) { + rigidAtoms.insert(atoms[i]->getGlobalIndex()); + } + for (int i = 0; i < static_cast(atoms.size()); ++i) { + atomGroups.insert(std::map >::value_type(atoms[i]->getGlobalIndex(), rigidAtoms)); + } + } else { + std::set oneAtomSet; + oneAtomSet.insert(integrableObject->getGlobalIndex()); + atomGroups.insert(std::map >::value_type(integrableObject->getGlobalIndex(), oneAtomSet)); + } + } + + for (bond= mol->beginBond(bondIter); bond != NULL; + bond = mol->nextBond(bondIter)) { + a = bond->getAtomA()->getGlobalIndex(); - b = bond->getAtomB()->getGlobalIndex(); - exclude_.removePair(a, b); + b = bond->getAtomB()->getGlobalIndex(); + + if (options_.havevdw12scale() || options_.haveelectrostatic12scale()) { + oneTwoInteractions_.removePair(a, b); + } else { + excludedInteractions_.removePair(a, b); + } } - for (bend= mol->beginBend(bendIter); bend != NULL; bend = mol->nextBend(bendIter)) { + for (bend= mol->beginBend(bendIter); bend != NULL; + bend = mol->nextBend(bendIter)) { + a = bend->getAtomA()->getGlobalIndex(); b = bend->getAtomB()->getGlobalIndex(); c = bend->getAtomC()->getGlobalIndex(); + + if (options_.havevdw12scale() || options_.haveelectrostatic12scale()) { + oneTwoInteractions_.removePair(a, b); + oneTwoInteractions_.removePair(b, c); + } else { + excludedInteractions_.removePair(a, b); + excludedInteractions_.removePair(b, c); + } - exclude_.removePair(a, b); - exclude_.removePair(a, c); - exclude_.removePair(b, c); + if (options_.havevdw13scale() || options_.haveelectrostatic13scale()) { + oneThreeInteractions_.removePair(a, c); + } else { + excludedInteractions_.removePair(a, c); + } } - for (torsion= mol->beginTorsion(torsionIter); torsion != NULL; torsion = mol->nextTorsion(torsionIter)) { + for (torsion= mol->beginTorsion(torsionIter); torsion != NULL; + torsion = mol->nextTorsion(torsionIter)) { + a = torsion->getAtomA()->getGlobalIndex(); b = torsion->getAtomB()->getGlobalIndex(); c = torsion->getAtomC()->getGlobalIndex(); - d = torsion->getAtomD()->getGlobalIndex(); + d = torsion->getAtomD()->getGlobalIndex(); + + if (options_.havevdw12scale() || options_.haveelectrostatic12scale()) { + oneTwoInteractions_.removePair(a, b); + oneTwoInteractions_.removePair(b, c); + oneTwoInteractions_.removePair(c, d); + } else { + excludedInteractions_.removePair(a, b); + excludedInteractions_.removePair(b, c); + excludedInteractions_.removePair(c, d); + } - exclude_.removePair(a, b); - exclude_.removePair(a, c); - exclude_.removePair(a, d); - exclude_.removePair(b, c); - exclude_.removePair(b, d); - exclude_.removePair(c, d); + if (options_.havevdw13scale() || options_.haveelectrostatic13scale()) { + oneThreeInteractions_.removePair(a, c); + oneThreeInteractions_.removePair(b, d); + } else { + excludedInteractions_.removePair(a, c); + excludedInteractions_.removePair(b, d); + } + + if (options_.havevdw14scale() || options_.haveelectrostatic14scale()) { + oneFourInteractions_.removePair(a, d); + } else { + excludedInteractions_.removePair(a, d); + } } - Molecule::RigidBodyIterator rbIter; - RigidBody* rb; - for (rb = mol->beginRigidBody(rbIter); rb != NULL; rb = mol->nextRigidBody(rbIter)) { + for (inversion= mol->beginInversion(inversionIter); inversion != NULL; + inversion = mol->nextInversion(inversionIter)) { + + a = inversion->getAtomA()->getGlobalIndex(); + b = inversion->getAtomB()->getGlobalIndex(); + c = inversion->getAtomC()->getGlobalIndex(); + d = inversion->getAtomD()->getGlobalIndex(); + + if (options_.havevdw12scale() || options_.haveelectrostatic12scale()) { + oneTwoInteractions_.removePair(a, b); + oneTwoInteractions_.removePair(a, c); + oneTwoInteractions_.removePair(a, d); + } else { + excludedInteractions_.removePair(a, b); + excludedInteractions_.removePair(a, c); + excludedInteractions_.removePair(a, d); + } + + if (options_.havevdw13scale() || options_.haveelectrostatic13scale()) { + oneThreeInteractions_.removePair(b, c); + oneThreeInteractions_.removePair(b, d); + oneThreeInteractions_.removePair(c, d); + } else { + excludedInteractions_.removePair(b, c); + excludedInteractions_.removePair(b, d); + excludedInteractions_.removePair(c, d); + } + } + + for (rb = mol->beginRigidBody(rbIter); rb != NULL; + rb = mol->nextRigidBody(rbIter)) { std::vector atoms = rb->getAtoms(); - for (int i = 0; i < atoms.size() -1 ; ++i) { - for (int j = i + 1; j < atoms.size(); ++j) { + for (int i = 0; i < static_cast(atoms.size()) -1 ; ++i) { + for (int j = i + 1; j < static_cast(atoms.size()); ++j) { a = atoms[i]->getGlobalIndex(); b = atoms[j]->getGlobalIndex(); - exclude_.removePair(a, b); + excludedInteractions_.removePair(a, b); } } } - + } - - + + void SimInfo::addMoleculeStamp(MoleculeStamp* molStamp, int nmol) { int curStampId; - + //index from 0 curStampId = moleculeStamps_.size(); @@ -465,7 +690,11 @@ namespace oopse { /** @deprecate */ int isError = 0; + setupCutoff(); + setupElectrostaticSummationMethod( isError ); + setupSwitchingFunction(); + setupAccumulateBoxDipole(); if(isError){ sprintf( painCave.errMsg, @@ -473,9 +702,6 @@ namespace oopse { painCave.isFatal = 1; simError(); } - - - setupCutoff(); calcNdf(); calcNdfRaw(); @@ -510,6 +736,7 @@ namespace oopse { int useLennardJones = 0; int useElectrostatic = 0; int useEAM = 0; + int useSC = 0; int useCharge = 0; int useDirectional = 0; int useDipole = 0; @@ -521,13 +748,45 @@ namespace oopse { int useDirectionalAtom = 0; int useElectrostatics = 0; //usePBC and useRF are from simParams - int usePBC = simParams_->getPBC(); + int usePBC = simParams_->getUsePeriodicBoundaryConditions(); + int useRF; + int useSF; + int useSP; + int useBoxDipole; + std::string myMethod; + + // set the useRF logical + useRF = 0; + useSF = 0; + useSP = 0; + useBoxDipole = 0; + + + if (simParams_->haveElectrostaticSummationMethod()) { + std::string myMethod = simParams_->getElectrostaticSummationMethod(); + toUpper(myMethod); + if (myMethod == "REACTION_FIELD"){ + useRF = 1; + } else if (myMethod == "SHIFTED_FORCE"){ + useSF = 1; + } else if (myMethod == "SHIFTED_POTENTIAL"){ + useSP = 1; + } + } + + if (simParams_->haveAccumulateBoxDipole()) + if (simParams_->getAccumulateBoxDipole()) + useBoxDipole = 1; + + useAtomicVirial_ = simParams_->getUseAtomicVirial(); + //loop over all of the atom types for (i = atomTypes.begin(); i != atomTypes.end(); ++i) { useLennardJones |= (*i)->isLennardJones(); useElectrostatic |= (*i)->isElectrostatic(); useEAM |= (*i)->isEAM(); + useSC |= (*i)->isSC(); useCharge |= (*i)->isCharge(); useDirectional |= (*i)->isDirectional(); useDipole |= (*i)->isDipole(); @@ -578,12 +837,30 @@ namespace oopse { temp = useEAM; MPI_Allreduce(&temp, &useEAM, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD); + temp = useSC; + MPI_Allreduce(&temp, &useSC, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD); + temp = useShape; MPI_Allreduce(&temp, &useShape, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD); temp = useFLARB; MPI_Allreduce(&temp, &useFLARB, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD); + temp = useRF; + MPI_Allreduce(&temp, &useRF, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD); + + temp = useSF; + MPI_Allreduce(&temp, &useSF, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD); + + temp = useSP; + MPI_Allreduce(&temp, &useSP, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD); + + temp = useBoxDipole; + MPI_Allreduce(&temp, &useBoxDipole, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD); + + temp = useAtomicVirial_; + MPI_Allreduce(&temp, &useAtomicVirial_, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD); + #endif fInfo_.SIM_uses_PBC = usePBC; @@ -596,34 +873,21 @@ namespace oopse { fInfo_.SIM_uses_StickyPower = useStickyPower; fInfo_.SIM_uses_GayBerne = useGayBerne; fInfo_.SIM_uses_EAM = useEAM; + fInfo_.SIM_uses_SC = useSC; fInfo_.SIM_uses_Shapes = useShape; fInfo_.SIM_uses_FLARB = useFLARB; - - if( fInfo_.SIM_uses_Dipoles && fInfo_.SIM_uses_RF) { - - if (simParams_->haveDielectric()) { - fInfo_.dielect = simParams_->getDielectric(); - } else { - sprintf(painCave.errMsg, - "SimSetup Error: No Dielectric constant was set.\n" - "\tYou are trying to use Reaction Field without" - "\tsetting a dielectric constant!\n"); - painCave.isFatal = 1; - simError(); - } - - } else { - fInfo_.dielect = 0.0; - } - + fInfo_.SIM_uses_RF = useRF; + fInfo_.SIM_uses_SF = useSF; + fInfo_.SIM_uses_SP = useSP; + fInfo_.SIM_uses_BoxDipole = useBoxDipole; + fInfo_.SIM_uses_AtomicVirial = useAtomicVirial_; } void SimInfo::setupFortranSim() { int isError; - int nExclude; + int nExclude, nOneTwo, nOneThree, nOneFour; std::vector fortranGlobalGroupMembership; - nExclude = exclude_.getSize(); isError = 0; //globalGroupMembership_ is filled by SimCreator @@ -632,14 +896,14 @@ namespace oopse { } //calculate mass ratio of cutoff group - std::vector mfact; + std::vector mfact; SimInfo::MoleculeIterator mi; Molecule* mol; Molecule::CutoffGroupIterator ci; CutoffGroup* cg; Molecule::AtomIterator ai; Atom* atom; - double totalMass; + RealType totalMass; //to avoid memory reallocation, reserve enough space for mfact mfact.reserve(getNCutoffGroups()); @@ -649,9 +913,12 @@ namespace oopse { totalMass = cg->getMass(); for(atom = cg->beginAtom(ai); atom != NULL; atom = cg->nextAtom(ai)) { - mfact.push_back(atom->getMass()/totalMass); + // Check for massless groups - set mfact to 1 if true + if (totalMass != 0) + mfact.push_back(atom->getMass()/totalMass); + else + mfact.push_back( 1.0 ); } - } } @@ -675,33 +942,52 @@ namespace oopse { } //setup fortran simulation - int nGlobalExcludes = 0; - int* globalExcludes = NULL; - int* excludeList = exclude_.getExcludeList(); - setFortranSim( &fInfo_, &nGlobalAtoms_, &nAtoms_, &identArray[0], &nExclude, excludeList , - &nGlobalExcludes, globalExcludes, &molMembershipArray[0], - &mfact[0], &nCutoffGroups_, &fortranGlobalGroupMembership[0], &isError); - if( isError ){ + nExclude = excludedInteractions_.getSize(); + nOneTwo = oneTwoInteractions_.getSize(); + nOneThree = oneThreeInteractions_.getSize(); + nOneFour = oneFourInteractions_.getSize(); + int* excludeList = excludedInteractions_.getPairList(); + int* oneTwoList = oneTwoInteractions_.getPairList(); + int* oneThreeList = oneThreeInteractions_.getPairList(); + int* oneFourList = oneFourInteractions_.getPairList(); + + setFortranSim( &fInfo_, &nGlobalAtoms_, &nAtoms_, &identArray[0], + &nExclude, excludeList, + &nOneTwo, oneTwoList, + &nOneThree, oneThreeList, + &nOneFour, oneFourList, + &molMembershipArray[0], &mfact[0], &nCutoffGroups_, + &fortranGlobalGroupMembership[0], &isError); + + if( isError ){ + sprintf( painCave.errMsg, "There was an error setting the simulation information in fortran.\n" ); painCave.isFatal = 1; - painCave.severity = OOPSE_ERROR; + painCave.severity = OPENMD_ERROR; simError(); } - -#ifdef IS_MPI + + sprintf( checkPointMsg, "succesfully sent the simulation information to fortran.\n"); - MPIcheckPoint(); -#endif // is_mpi + + errorCheckPoint(); + + // Setup number of neighbors in neighbor list if present + if (simParams_->haveNeighborListNeighbors()) { + int nlistNeighbors = simParams_->getNeighborListNeighbors(); + setNeighbors(&nlistNeighbors); + } + + } -#ifdef IS_MPI void SimInfo::setupFortranParallel() { - +#ifdef IS_MPI //SimInfo is responsible for creating localToGlobalAtomIndex and localToGlobalGroupIndex std::vector localToGlobalAtomIndex(getNAtoms(), 0); std::vector localToGlobalCutoffGroupIndex; @@ -751,93 +1037,31 @@ namespace oopse { } sprintf(checkPointMsg, " mpiRefresh successful.\n"); - MPIcheckPoint(); + errorCheckPoint(); - - } - #endif - - double SimInfo::calcMaxCutoffRadius() { - - - std::set atomTypes; - std::set::iterator i; - std::vector cutoffRadius; - - //get the unique atom types - atomTypes = getUniqueAtomTypes(); - - //query the max cutoff radius among these atom types - for (i = atomTypes.begin(); i != atomTypes.end(); ++i) { - cutoffRadius.push_back(forceField_->getRcutFromAtomType(*i)); - } - - double maxCutoffRadius = *(std::max_element(cutoffRadius.begin(), cutoffRadius.end())); -#ifdef IS_MPI - //pick the max cutoff radius among the processors -#endif - - return maxCutoffRadius; } - void SimInfo::getCutoff(double& rcut, double& rsw) { + void SimInfo::setupCutoff() { - if (fInfo_.SIM_uses_Charges | fInfo_.SIM_uses_Dipoles | fInfo_.SIM_uses_RF) { - - if (!simParams_->haveRcut()){ - sprintf(painCave.errMsg, - "SimCreator Warning: No value was set for the cutoffRadius.\n" - "\tOOPSE will use a default value of 15.0 angstroms" - "\tfor the cutoffRadius.\n"); - painCave.isFatal = 0; - simError(); - rcut = 15.0; - } else{ - rcut = simParams_->getRcut(); - } + ForceFieldOptions& forceFieldOptions_ = forceField_->getForceFieldOptions(); - if (!simParams_->haveRsw()){ - sprintf(painCave.errMsg, - "SimCreator Warning: No value was set for switchingRadius.\n" - "\tOOPSE will use a default value of\n" - "\t0.95 * cutoffRadius for the switchingRadius\n"); - painCave.isFatal = 0; - simError(); - rsw = 0.95 * rcut; - } else{ - rsw = simParams_->getRsw(); - } + // Check the cutoff policy + int cp = TRADITIONAL_CUTOFF_POLICY; // Set to traditional by default - } else { - // if charge, dipole or reaction field is not used and the cutofff radius is not specified in - //meta-data file, the maximum cutoff radius calculated from forcefiled will be used - - if (simParams_->haveRcut()) { - rcut = simParams_->getRcut(); - } else { - //set cutoff radius to the maximum cutoff radius based on atom types in the whole system - rcut = calcMaxCutoffRadius(); - } + // Set LJ shifting bools to false + ljsp_ = 0; + ljsf_ = 0; - if (simParams_->haveRsw()) { - rsw = simParams_->getRsw(); - } else { - rsw = rcut; - } - + std::string myPolicy; + if (forceFieldOptions_.haveCutoffPolicy()){ + myPolicy = forceFieldOptions_.getCutoffPolicy(); + }else if (simParams_->haveCutoffPolicy()) { + myPolicy = simParams_->getCutoffPolicy(); } - } - void SimInfo::setupCutoff() { - getCutoff(rcut_, rsw_); - double rnblist = rcut_ + 1; // skin of neighbor list - - //Pass these cutoff radius etc. to fortran. This function should be called once and only once - - int cp = TRADITIONAL_CUTOFF_POLICY; - if (simParams_->haveCutoffPolicy()) { - std::string myPolicy = simParams_->getCutoffPolicy(); + if (!myPolicy.empty()){ + toUpper(myPolicy); if (myPolicy == "MIX") { cp = MIX_CUTOFF_POLICY; } else { @@ -855,51 +1079,251 @@ namespace oopse { } } } + } + notifyFortranCutoffPolicy(&cp); + + // Check the Skin Thickness for neighborlists + RealType skin; + if (simParams_->haveSkinThickness()) { + skin = simParams_->getSkinThickness(); + notifyFortranSkinThickness(&skin); + } + + // Check if the cutoff was set explicitly: + if (simParams_->haveCutoffRadius()) { + rcut_ = simParams_->getCutoffRadius(); + if (simParams_->haveSwitchingRadius()) { + rsw_ = simParams_->getSwitchingRadius(); + } else { + if (fInfo_.SIM_uses_Charges | + fInfo_.SIM_uses_Dipoles | + fInfo_.SIM_uses_RF) { + + rsw_ = 0.85 * rcut_; + sprintf(painCave.errMsg, + "SimCreator Warning: No value was set for the switchingRadius.\n" + "\tOpenMD will use a default value of 85 percent of the cutoffRadius.\n" + "\tswitchingRadius = %f. for this simulation\n", rsw_); + painCave.isFatal = 0; + simError(); + } else { + rsw_ = rcut_; + sprintf(painCave.errMsg, + "SimCreator Warning: No value was set for the switchingRadius.\n" + "\tOpenMD will use the same value as the cutoffRadius.\n" + "\tswitchingRadius = %f. for this simulation\n", rsw_); + painCave.isFatal = 0; + simError(); + } + } + + if (simParams_->haveElectrostaticSummationMethod()) { + std::string myMethod = simParams_->getElectrostaticSummationMethod(); + toUpper(myMethod); + + if (myMethod == "SHIFTED_POTENTIAL") { + ljsp_ = 1; + } else if (myMethod == "SHIFTED_FORCE") { + ljsf_ = 1; + } + } + + notifyFortranCutoffs(&rcut_, &rsw_, &ljsp_, &ljsf_); + + } else { + + // For electrostatic atoms, we'll assume a large safe value: + if (fInfo_.SIM_uses_Charges | fInfo_.SIM_uses_Dipoles | fInfo_.SIM_uses_RF) { + sprintf(painCave.errMsg, + "SimCreator Warning: No value was set for the cutoffRadius.\n" + "\tOpenMD will use a default value of 15.0 angstroms" + "\tfor the cutoffRadius.\n"); + painCave.isFatal = 0; + simError(); + rcut_ = 15.0; + + if (simParams_->haveElectrostaticSummationMethod()) { + std::string myMethod = simParams_->getElectrostaticSummationMethod(); + toUpper(myMethod); + + // For the time being, we're tethering the LJ shifted behavior to the + // electrostaticSummationMethod keyword options + if (myMethod == "SHIFTED_POTENTIAL") { + ljsp_ = 1; + } else if (myMethod == "SHIFTED_FORCE") { + ljsf_ = 1; + } + if (myMethod == "SHIFTED_POTENTIAL" || myMethod == "SHIFTED_FORCE") { + if (simParams_->haveSwitchingRadius()){ + sprintf(painCave.errMsg, + "SimInfo Warning: A value was set for the switchingRadius\n" + "\teven though the electrostaticSummationMethod was\n" + "\tset to %s\n", myMethod.c_str()); + painCave.isFatal = 1; + simError(); + } + } + } + + if (simParams_->haveSwitchingRadius()){ + rsw_ = simParams_->getSwitchingRadius(); + } else { + sprintf(painCave.errMsg, + "SimCreator Warning: No value was set for switchingRadius.\n" + "\tOpenMD will use a default value of\n" + "\t0.85 * cutoffRadius for the switchingRadius\n"); + painCave.isFatal = 0; + simError(); + rsw_ = 0.85 * rcut_; + } + + Electrostatic::setElectrostaticCutoffRadius(rcut_, rsw_); + notifyFortranCutoffs(&rcut_, &rsw_, &ljsp_, &ljsf_); + + } else { + // We didn't set rcut explicitly, and we don't have electrostatic atoms, so + // We'll punt and let fortran figure out the cutoffs later. + + notifyFortranYouAreOnYourOwn(); + + } } - notifyFortranCutoffs(&rcut_, &rsw_, &rnblist, &cp); } void SimInfo::setupElectrostaticSummationMethod( int isError ) { int errorOut; - int esm = NONE; - double alphaVal; - + ElectrostaticSummationMethod esm = NONE; + ElectrostaticScreeningMethod sm = UNDAMPED; + RealType alphaVal; + RealType dielectric; + errorOut = isError; if (simParams_->haveElectrostaticSummationMethod()) { std::string myMethod = simParams_->getElectrostaticSummationMethod(); + toUpper(myMethod); if (myMethod == "NONE") { esm = NONE; } else { - if (myMethod == "UNDAMPED_WOLF") { - esm = UNDAMPED_WOLF; + if (myMethod == "SWITCHING_FUNCTION") { + esm = SWITCHING_FUNCTION; } else { - if (myMethod == "DAMPED_WOLF") { - esm = DAMPED_WOLF; - if (!simParams_->haveDampingAlpha()) { - //throw error - sprintf( painCave.errMsg, - "SimInfo warning: dampingAlpha was not specified in the input file. A default value of %f (1/ang) will be used for the Damped Wolf Method.", simParams_->getDampingAlpha()); - painCave.isFatal = 0; - simError(); - } - alphaVal = simParams_->getDampingAlpha(); - } else { - if (myMethod == "REACTION_FIELD") { - esm = REACTION_FIELD; + if (myMethod == "SHIFTED_POTENTIAL") { + esm = SHIFTED_POTENTIAL; + } else { + if (myMethod == "SHIFTED_FORCE") { + esm = SHIFTED_FORCE; } else { - // throw error - sprintf( painCave.errMsg, - "SimInfo error: Unknown electrostaticSummationMethod. (Input file specified %s .)\n\telectrostaticSummationMethod must be one of: \"none\", \"undamped_wolf\", \"damped_wolf\", or \"reaction_field\".", myMethod.c_str() ); - painCave.isFatal = 1; - simError(); - } - } + if (myMethod == "REACTION_FIELD") { + esm = REACTION_FIELD; + dielectric = simParams_->getDielectric(); + if (!simParams_->haveDielectric()) { + // throw warning + sprintf( painCave.errMsg, + "SimInfo warning: dielectric was not specified in the input file\n\tfor the reaction field correction method.\n" + "\tA default value of %f will be used for the dielectric.\n", dielectric); + painCave.isFatal = 0; + simError(); + } + } else { + // throw error + sprintf( painCave.errMsg, + "SimInfo error: Unknown electrostaticSummationMethod.\n" + "\t(Input file specified %s .)\n" + "\telectrostaticSummationMethod must be one of: \"none\",\n" + "\t\"shifted_potential\", \"shifted_force\", or \n" + "\t\"reaction_field\".\n", myMethod.c_str() ); + painCave.isFatal = 1; + simError(); + } + } + } } } } - initFortranFF( &esm, &alphaVal, &errorOut ); + + if (simParams_->haveElectrostaticScreeningMethod()) { + std::string myScreen = simParams_->getElectrostaticScreeningMethod(); + toUpper(myScreen); + if (myScreen == "UNDAMPED") { + sm = UNDAMPED; + } else { + if (myScreen == "DAMPED") { + sm = DAMPED; + if (!simParams_->haveDampingAlpha()) { + // first set a cutoff dependent alpha value + // we assume alpha depends linearly with rcut from 0 to 20.5 ang + alphaVal = 0.5125 - rcut_* 0.025; + // for values rcut > 20.5, alpha is zero + if (alphaVal < 0) alphaVal = 0; + + // throw warning + sprintf( painCave.errMsg, + "SimInfo warning: dampingAlpha was not specified in the input file.\n" + "\tA default value of %f (1/ang) will be used for the cutoff of\n\t%f (ang).\n", alphaVal, rcut_); + painCave.isFatal = 0; + simError(); + } else { + alphaVal = simParams_->getDampingAlpha(); + } + + } else { + // throw error + sprintf( painCave.errMsg, + "SimInfo error: Unknown electrostaticScreeningMethod.\n" + "\t(Input file specified %s .)\n" + "\telectrostaticScreeningMethod must be one of: \"undamped\"\n" + "or \"damped\".\n", myScreen.c_str() ); + painCave.isFatal = 1; + simError(); + } + } + } + + + Electrostatic::setElectrostaticSummationMethod( esm ); + Electrostatic::setElectrostaticScreeningMethod( sm ); + Electrostatic::setDampingAlpha( alphaVal ); + Electrostatic::setReactionFieldDielectric( dielectric ); + initFortranFF( &errorOut ); + } + + void SimInfo::setupSwitchingFunction() { + int ft = CUBIC; + + if (simParams_->haveSwitchingFunctionType()) { + std::string funcType = simParams_->getSwitchingFunctionType(); + toUpper(funcType); + if (funcType == "CUBIC") { + ft = CUBIC; + } else { + if (funcType == "FIFTH_ORDER_POLYNOMIAL") { + ft = FIFTH_ORDER_POLY; + } else { + // throw error + sprintf( painCave.errMsg, + "SimInfo error: Unknown switchingFunctionType. (Input file specified %s .)\n\tswitchingFunctionType must be one of: \"cubic\" or \"fifth_order_polynomial\".", funcType.c_str() ); + painCave.isFatal = 1; + simError(); + } + } + } + + // send switching function notification to switcheroo + setFunctionType(&ft); + + } + + void SimInfo::setupAccumulateBoxDipole() { + + // we only call setAccumulateBoxDipole if the accumulateBoxDipole parameter is true + if ( simParams_->haveAccumulateBoxDipole() ) + if ( simParams_->getAccumulateBoxDipole() ) { + setAccumulateBoxDipole(); + calcBoxDipole_ = true; + } + } void SimInfo::addProperty(GenericData* genData) { @@ -958,20 +1382,20 @@ namespace oopse { Molecule* mol; Vector3d comVel(0.0); - double totalMass = 0.0; + RealType totalMass = 0.0; for (mol = beginMolecule(i); mol != NULL; mol = nextMolecule(i)) { - double mass = mol->getMass(); + RealType mass = mol->getMass(); totalMass += mass; comVel += mass * mol->getComVel(); } #ifdef IS_MPI - double tmpMass = totalMass; + RealType tmpMass = totalMass; Vector3d tmpComVel(comVel); - MPI_Allreduce(&tmpMass,&totalMass,1,MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD); - MPI_Allreduce(tmpComVel.getArrayPointer(), comVel.getArrayPointer(),3,MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD); + MPI_Allreduce(&tmpMass,&totalMass,1,MPI_REALTYPE,MPI_SUM, MPI_COMM_WORLD); + MPI_Allreduce(tmpComVel.getArrayPointer(), comVel.getArrayPointer(),3,MPI_REALTYPE,MPI_SUM, MPI_COMM_WORLD); #endif comVel /= totalMass; @@ -984,19 +1408,19 @@ namespace oopse { Molecule* mol; Vector3d com(0.0); - double totalMass = 0.0; + RealType totalMass = 0.0; for (mol = beginMolecule(i); mol != NULL; mol = nextMolecule(i)) { - double mass = mol->getMass(); + RealType mass = mol->getMass(); totalMass += mass; com += mass * mol->getCom(); } #ifdef IS_MPI - double tmpMass = totalMass; + RealType tmpMass = totalMass; Vector3d tmpCom(com); - MPI_Allreduce(&tmpMass,&totalMass,1,MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD); - MPI_Allreduce(tmpCom.getArrayPointer(), com.getArrayPointer(),3,MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD); + MPI_Allreduce(&tmpMass,&totalMass,1,MPI_REALTYPE,MPI_SUM, MPI_COMM_WORLD); + MPI_Allreduce(tmpCom.getArrayPointer(), com.getArrayPointer(),3,MPI_REALTYPE,MPI_SUM, MPI_COMM_WORLD); #endif com /= totalMass; @@ -1020,23 +1444,23 @@ namespace oopse { Molecule* mol; - double totalMass = 0.0; + RealType totalMass = 0.0; for (mol = beginMolecule(i); mol != NULL; mol = nextMolecule(i)) { - double mass = mol->getMass(); + RealType mass = mol->getMass(); totalMass += mass; com += mass * mol->getCom(); comVel += mass * mol->getComVel(); } #ifdef IS_MPI - double tmpMass = totalMass; + RealType tmpMass = totalMass; Vector3d tmpCom(com); Vector3d tmpComVel(comVel); - MPI_Allreduce(&tmpMass,&totalMass,1,MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD); - MPI_Allreduce(tmpCom.getArrayPointer(), com.getArrayPointer(),3,MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD); - MPI_Allreduce(tmpComVel.getArrayPointer(), comVel.getArrayPointer(),3,MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD); + MPI_Allreduce(&tmpMass,&totalMass,1,MPI_REALTYPE,MPI_SUM, MPI_COMM_WORLD); + MPI_Allreduce(tmpCom.getArrayPointer(), com.getArrayPointer(),3,MPI_REALTYPE,MPI_SUM, MPI_COMM_WORLD); + MPI_Allreduce(tmpComVel.getArrayPointer(), comVel.getArrayPointer(),3,MPI_REALTYPE,MPI_SUM, MPI_COMM_WORLD); #endif com /= totalMass; @@ -1055,12 +1479,12 @@ namespace oopse { void SimInfo::getInertiaTensor(Mat3x3d &inertiaTensor, Vector3d &angularMomentum){ - double xx = 0.0; - double yy = 0.0; - double zz = 0.0; - double xy = 0.0; - double xz = 0.0; - double yz = 0.0; + 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); @@ -1072,7 +1496,7 @@ namespace oopse { Vector3d thisq(0.0); Vector3d thisv(0.0); - double thisMass = 0.0; + RealType thisMass = 0.0; @@ -1110,8 +1534,8 @@ namespace oopse { #ifdef IS_MPI Mat3x3d tmpI(inertiaTensor); Vector3d tmpAngMom; - MPI_Allreduce(tmpI.getArrayPointer(), inertiaTensor.getArrayPointer(),9,MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD); - MPI_Allreduce(tmpAngMom.getArrayPointer(), angularMomentum.getArrayPointer(),3,MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD); + MPI_Allreduce(tmpI.getArrayPointer(), inertiaTensor.getArrayPointer(),9,MPI_REALTYPE,MPI_SUM, MPI_COMM_WORLD); + MPI_Allreduce(tmpAngMom.getArrayPointer(), angularMomentum.getArrayPointer(),3,MPI_REALTYPE,MPI_SUM, MPI_COMM_WORLD); #endif return; @@ -1132,7 +1556,7 @@ namespace oopse { Vector3d thisr(0.0); Vector3d thisp(0.0); - double thisMass; + RealType thisMass; for (mol = beginMolecule(i); mol != NULL; mol = nextMolecule(i)) { thisMass = mol->getMass(); @@ -1145,12 +1569,67 @@ namespace oopse { #ifdef IS_MPI Vector3d tmpAngMom; - MPI_Allreduce(tmpAngMom.getArrayPointer(), angularMomentum.getArrayPointer(),3,MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD); + MPI_Allreduce(tmpAngMom.getArrayPointer(), angularMomentum.getArrayPointer(),3,MPI_REALTYPE,MPI_SUM, MPI_COMM_WORLD); #endif return angularMomentum; } - -}//end namespace oopse + StuntDouble* SimInfo::getIOIndexToIntegrableObject(int index) { + return IOIndexToIntegrableObject.at(index); + } + + void SimInfo::setIOIndexToIntegrableObject(const std::vector& v) { + IOIndexToIntegrableObject= v; + } + /* Returns the Volume of the simulation 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. + */ + void SimInfo::getGyrationalVolume(RealType &volume){ + Mat3x3d intTensor; + RealType det; + Vector3d dummyAngMom; + RealType sysconstants; + RealType geomCnst; + + geomCnst = 3.0/2.0; + /* Get the inertial tensor and angular momentum for free*/ + getInertiaTensor(intTensor,dummyAngMom); + + det = intTensor.determinant(); + sysconstants = geomCnst/(RealType)nGlobalIntegrableObjects_; + volume = 4.0/3.0*NumericConstant::PI*pow(sysconstants,3.0/2.0)*sqrt(det); + return; + } + + void SimInfo::getGyrationalVolume(RealType &volume, RealType &detI){ + Mat3x3d intTensor; + Vector3d dummyAngMom; + RealType sysconstants; + RealType geomCnst; + + geomCnst = 3.0/2.0; + /* Get the inertial tensor and angular momentum for free*/ + getInertiaTensor(intTensor,dummyAngMom); + + detI = intTensor.determinant(); + sysconstants = geomCnst/(RealType)nGlobalIntegrableObjects_; + volume = 4.0/3.0*NumericConstant::PI*pow(sysconstants,3.0/2.0)*sqrt(detI); + return; + } +/* + void SimInfo::setStuntDoubleFromGlobalIndex(std::vector v) { + assert( v.size() == nAtoms_ + nRigidBodies_); + sdByGlobalIndex_ = v; + } + + StuntDouble* SimInfo::getStuntDoubleFromGlobalIndex(int index) { + //assert(index < nAtoms_ + nRigidBodies_); + return sdByGlobalIndex_.at(index); + } +*/ +}//end namespace OpenMD +