--- trunk/src/brains/SimInfo.cpp 2007/02/26 04:45:42 1121 +++ branches/development/src/brains/SimInfo.cpp 2010/12/29 19:59:21 1532 @@ -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). */ /** @@ -54,19 +54,14 @@ #include "math/Vector3.hpp" #include "primitives/Molecule.hpp" #include "primitives/StuntDouble.hpp" -#include "UseTheForce/fCutoffPolicy.h" -#include "UseTheForce/DarkSide/fElectrostaticSummationMethod.h" -#include "UseTheForce/DarkSide/fElectrostaticScreeningMethod.h" -#include "UseTheForce/DarkSide/fSwitchingFunctionType.h" #include "UseTheForce/doForces_interface.h" #include "UseTheForce/DarkSide/neighborLists_interface.h" -#include "UseTheForce/DarkSide/electrostatic_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" +#include "nonbonded/SwitchingFunction.hpp" #ifdef IS_MPI @@ -74,96 +69,86 @@ #include "UseTheForce/DarkSide/simParallel_interface.h" #endif -namespace oopse { - 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; - } - - return result; - } +using namespace std; +namespace OpenMD { 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), calcBoxDipole_(false) { - - 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 - CutoffGroupStamp* cgStamp; - RigidBodyStamp* rbStamp; - int nRigidAtoms = 0; - std::vector components = simParams->getComponents(); + 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) { + + 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 + CutoffGroupStamp* cgStamp; + RigidBodyStamp* rbStamp; + int nRigidAtoms = 0; + + vector components = simParams->getComponents(); + + for (vector::iterator i = components.begin(); i !=components.end(); ++i) { + molStamp = (*i)->getMoleculeStamp(); + nMolWithSameStamp = (*i)->getNMol(); - 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->getCutoffGroupStamp(j); - nAtomsInGroups += cgStamp->getNMembers(); - } - - nGroups += nCutoffGroupsInStamp * nMolWithSameStamp; - - nCutoffAtoms += nAtomsInGroups * nMolWithSameStamp; - - //calculate atoms in rigid bodies - int nAtomsInRigidBodies = 0; - int nRigidBodiesInStamp = molStamp->getNRigidBodies(); - - for (int j=0; j < nRigidBodiesInStamp; j++) { - rbStamp = molStamp->getRigidBodyStamp(j); - nAtomsInRigidBodies += rbStamp->getNMembers(); - } - - nGlobalRigidBodies_ += nRigidBodiesInStamp * nMolWithSameStamp; - nRigidAtoms += nAtomsInRigidBodies * nMolWithSameStamp; - + 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->getCutoffGroupStamp(j); + nAtomsInGroups += cgStamp->getNMembers(); } - - //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_; - - nGlobalMols_ = molStampIds_.size(); - -#ifdef IS_MPI - molToProcMap_.resize(nGlobalMols_); -#endif - + + nGroups += nCutoffGroupsInStamp * nMolWithSameStamp; + + nCutoffAtoms += nAtomsInGroups * nMolWithSameStamp; + + //calculate atoms in rigid bodies + int nAtomsInRigidBodies = 0; + int nRigidBodiesInStamp = molStamp->getNRigidBodies(); + + for (int j=0; j < nRigidBodiesInStamp; j++) { + rbStamp = molStamp->getRigidBodyStamp(j); + nAtomsInRigidBodies += rbStamp->getNMembers(); + } + + nGlobalRigidBodies_ += nRigidBodiesInStamp * nMolWithSameStamp; + nRigidAtoms += nAtomsInRigidBodies * nMolWithSameStamp; + } - + + //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_; + + nGlobalMols_ = molStampIds_.size(); + molToProcMap_.resize(nGlobalMols_); + } + SimInfo::~SimInfo() { - std::map::iterator i; + map::iterator i; for (i = molecules_.begin(); i != molecules_.end(); ++i) { delete i->second; } @@ -174,42 +159,33 @@ namespace oopse { delete forceField_; } - int SimInfo::getNGlobalConstraints() { - int nGlobalConstraints; -#ifdef IS_MPI - MPI_Allreduce(&nConstraints_, &nGlobalConstraints, 1, MPI_INT, MPI_SUM, - MPI_COMM_WORLD); -#else - nGlobalConstraints = nConstraints_; -#endif - return nGlobalConstraints; - } bool SimInfo::addMolecule(Molecule* mol) { MoleculeIterator i; - + i = molecules_.find(mol->getGlobalIndex()); if (i == molecules_.end() ) { - - molecules_.insert(std::make_pair(mol->getGlobalIndex(), mol)); - + + molecules_.insert(make_pair(mol->getGlobalIndex(), mol)); + nAtoms_ += mol->getNAtoms(); 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; } } - + bool SimInfo::removeMolecule(Molecule* mol) { MoleculeIterator i; i = molecules_.find(mol->getGlobalIndex()); @@ -222,12 +198,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; @@ -236,8 +213,6 @@ namespace oopse { } else { return false; } - - } @@ -255,7 +230,7 @@ namespace oopse { void SimInfo::calcNdf() { int ndf_local; MoleculeIterator i; - std::vector::iterator j; + vector::iterator j; Molecule* mol; StuntDouble* integrableObject; @@ -306,7 +281,7 @@ namespace oopse { int ndfRaw_local; MoleculeIterator i; - std::vector::iterator j; + vector::iterator j; Molecule* mol; StuntDouble* integrableObject; @@ -353,235 +328,324 @@ namespace oopse { } - void SimInfo::addExcludePairs(Molecule* mol) { - std::vector::iterator bondIter; - std::vector::iterator bendIter; - std::vector::iterator torsionIter; + void SimInfo::addInteractionPairs(Molecule* mol) { + ForceFieldOptions& options_ = forceField_->getForceFieldOptions(); + vector::iterator bondIter; + vector::iterator bendIter; + vector::iterator torsionIter; + vector::iterator inversionIter; Bond* bond; Bend* bend; Torsion* torsion; + Inversion* inversion; int a; int b; int c; int d; - std::map > atomGroups; + // 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. + map > atomGroups; Molecule::RigidBodyIterator rbIter; RigidBody* rb; Molecule::IntegrableObjectIterator ii; StuntDouble* integrableObject; - for (integrableObject = mol->beginIntegrableObject(ii); integrableObject != NULL; - integrableObject = mol->nextIntegrableObject(ii)) { - + 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 < atoms.size(); ++i) { - rigidAtoms.insert(atoms[i]->getGlobalIndex()); - } - for (int i = 0; i < atoms.size(); ++i) { - atomGroups.insert(std::map >::value_type(atoms[i]->getGlobalIndex(), rigidAtoms)); - } + rb = static_cast(integrableObject); + vector atoms = rb->getAtoms(); + 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(map >::value_type(atoms[i]->getGlobalIndex(), rigidAtoms)); + } } else { - std::set oneAtomSet; + set oneAtomSet; oneAtomSet.insert(integrableObject->getGlobalIndex()); - atomGroups.insert(std::map >::value_type(integrableObject->getGlobalIndex(), oneAtomSet)); + atomGroups.insert(map >::value_type(integrableObject->getGlobalIndex(), oneAtomSet)); } } + + for (bond= mol->beginBond(bondIter); bond != NULL; + bond = mol->nextBond(bondIter)) { - - - for (bond= mol->beginBond(bondIter); bond != NULL; bond = mol->nextBond(bondIter)) { a = bond->getAtomA()->getGlobalIndex(); - b = bond->getAtomB()->getGlobalIndex(); - exclude_.addPair(a, b); + 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(); - std::set rigidSetA = getRigidSet(a, atomGroups); - std::set rigidSetB = getRigidSet(b, atomGroups); - std::set rigidSetC = getRigidSet(c, atomGroups); - - exclude_.addPairs(rigidSetA, rigidSetB); - exclude_.addPairs(rigidSetA, rigidSetC); - exclude_.addPairs(rigidSetB, rigidSetC); - //exclude_.addPair(a, b); - //exclude_.addPair(a, c); - //exclude_.addPair(b, c); + if (options_.havevdw12scale() || options_.haveelectrostatic12scale()) { + oneTwoInteractions_.addPair(a, b); + oneTwoInteractions_.addPair(b, c); + } else { + excludedInteractions_.addPair(a, b); + excludedInteractions_.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(); - std::set rigidSetA = getRigidSet(a, atomGroups); - std::set rigidSetB = getRigidSet(b, atomGroups); - std::set rigidSetC = getRigidSet(c, atomGroups); - std::set rigidSetD = getRigidSet(d, atomGroups); + d = torsion->getAtomD()->getGlobalIndex(); - exclude_.addPairs(rigidSetA, rigidSetB); - exclude_.addPairs(rigidSetA, rigidSetC); - exclude_.addPairs(rigidSetA, rigidSetD); - exclude_.addPairs(rigidSetB, rigidSetC); - exclude_.addPairs(rigidSetB, rigidSetD); - exclude_.addPairs(rigidSetC, rigidSetD); + 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); + } - /* - exclude_.addPairs(rigidSetA.begin(), rigidSetA.end(), rigidSetB.begin(), rigidSetB.end()); - exclude_.addPairs(rigidSetA.begin(), rigidSetA.end(), rigidSetC.begin(), rigidSetC.end()); - exclude_.addPairs(rigidSetA.begin(), rigidSetA.end(), rigidSetD.begin(), rigidSetD.end()); - exclude_.addPairs(rigidSetB.begin(), rigidSetB.end(), rigidSetC.begin(), rigidSetC.end()); - exclude_.addPairs(rigidSetB.begin(), rigidSetB.end(), rigidSetD.begin(), rigidSetD.end()); - exclude_.addPairs(rigidSetC.begin(), rigidSetC.end(), rigidSetD.begin(), rigidSetD.end()); - - - 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_.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); + } } - 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 (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)) { + vector atoms = rb->getAtoms(); + 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) { - std::vector::iterator bondIter; - std::vector::iterator bendIter; - std::vector::iterator torsionIter; + void SimInfo::removeInteractionPairs(Molecule* mol) { + ForceFieldOptions& options_ = forceField_->getForceFieldOptions(); + vector::iterator bondIter; + vector::iterator bendIter; + vector::iterator torsionIter; + vector::iterator inversionIter; Bond* bond; Bend* bend; Torsion* torsion; + Inversion* inversion; int a; int b; int c; int d; - std::map > atomGroups; - + map > atomGroups; Molecule::RigidBodyIterator rbIter; RigidBody* rb; Molecule::IntegrableObjectIterator ii; StuntDouble* integrableObject; - for (integrableObject = mol->beginIntegrableObject(ii); integrableObject != NULL; - integrableObject = mol->nextIntegrableObject(ii)) { - + 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 < atoms.size(); ++i) { - rigidAtoms.insert(atoms[i]->getGlobalIndex()); - } - for (int i = 0; i < atoms.size(); ++i) { - atomGroups.insert(std::map >::value_type(atoms[i]->getGlobalIndex(), rigidAtoms)); - } + rb = static_cast(integrableObject); + vector atoms = rb->getAtoms(); + 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(map >::value_type(atoms[i]->getGlobalIndex(), rigidAtoms)); + } } else { - std::set oneAtomSet; + set oneAtomSet; oneAtomSet.insert(integrableObject->getGlobalIndex()); - atomGroups.insert(std::map >::value_type(integrableObject->getGlobalIndex(), oneAtomSet)); + atomGroups.insert(map >::value_type(integrableObject->getGlobalIndex(), oneAtomSet)); } } - - for (bond= mol->beginBond(bondIter); bond != NULL; bond = mol->nextBond(bondIter)) { + 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(); - - std::set rigidSetA = getRigidSet(a, atomGroups); - std::set rigidSetB = getRigidSet(b, atomGroups); - std::set rigidSetC = getRigidSet(c, atomGroups); - - exclude_.removePairs(rigidSetA, rigidSetB); - exclude_.removePairs(rigidSetA, rigidSetC); - exclude_.removePairs(rigidSetB, rigidSetC); - //exclude_.removePair(a, b); - //exclude_.removePair(a, c); - //exclude_.removePair(b, c); + if (options_.havevdw12scale() || options_.haveelectrostatic12scale()) { + oneTwoInteractions_.removePair(a, b); + oneTwoInteractions_.removePair(b, c); + } else { + excludedInteractions_.removePair(a, b); + excludedInteractions_.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); + } - std::set rigidSetA = getRigidSet(a, atomGroups); - std::set rigidSetB = getRigidSet(b, atomGroups); - std::set rigidSetC = getRigidSet(c, atomGroups); - std::set rigidSetD = getRigidSet(d, atomGroups); + if (options_.havevdw13scale() || options_.haveelectrostatic13scale()) { + oneThreeInteractions_.removePair(a, c); + oneThreeInteractions_.removePair(b, d); + } else { + excludedInteractions_.removePair(a, c); + excludedInteractions_.removePair(b, d); + } - exclude_.removePairs(rigidSetA, rigidSetB); - exclude_.removePairs(rigidSetA, rigidSetC); - exclude_.removePairs(rigidSetA, rigidSetD); - exclude_.removePairs(rigidSetB, rigidSetC); - exclude_.removePairs(rigidSetB, rigidSetD); - exclude_.removePairs(rigidSetC, rigidSetD); + if (options_.havevdw14scale() || options_.haveelectrostatic14scale()) { + oneFourInteractions_.removePair(a, d); + } else { + excludedInteractions_.removePair(a, d); + } + } - /* - exclude_.removePairs(rigidSetA.begin(), rigidSetA.end(), rigidSetB.begin(), rigidSetB.end()); - exclude_.removePairs(rigidSetA.begin(), rigidSetA.end(), rigidSetC.begin(), rigidSetC.end()); - exclude_.removePairs(rigidSetA.begin(), rigidSetA.end(), rigidSetD.begin(), rigidSetD.end()); - exclude_.removePairs(rigidSetB.begin(), rigidSetB.end(), rigidSetC.begin(), rigidSetC.end()); - exclude_.removePairs(rigidSetB.begin(), rigidSetB.end(), rigidSetD.begin(), rigidSetD.end()); - exclude_.removePairs(rigidSetC.begin(), rigidSetC.end(), rigidSetD.begin(), rigidSetD.end()); + for (inversion= mol->beginInversion(inversionIter); inversion != NULL; + inversion = mol->nextInversion(inversionIter)) { - - exclude_.removePair(a, b); - exclude_.removePair(a, c); - exclude_.removePair(a, d); - exclude_.removePair(b, c); - exclude_.removePair(b, d); - exclude_.removePair(c, d); - */ + 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 (rb = mol->beginRigidBody(rbIter); rb != NULL; + rb = mol->nextRigidBody(rbIter)) { + vector atoms = rb->getAtoms(); + 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(); @@ -589,211 +653,261 @@ namespace oopse { molStampIds_.insert(molStampIds_.end(), nmol, curStampId); } + + /** + * update + * + * Performs the global checks and variable settings after the objects have been + * created. + * + */ void SimInfo::update() { + + setupSimVariables(); + setupCutoffs(); + setupSwitching(); + setupElectrostatics(); + setupNeighborlists(); - setupSimType(); - #ifdef IS_MPI setupFortranParallel(); #endif - setupFortranSim(); + fortranInitialized_ = true; - //setup fortran force field - /** @deprecate */ - int isError = 0; - - setupCutoff(); - - setupElectrostaticSummationMethod( isError ); - setupSwitchingFunction(); - setupAccumulateBoxDipole(); - - if(isError){ - sprintf( painCave.errMsg, - "ForceField error: There was an error initializing the forceField in fortran.\n" ); - painCave.isFatal = 1; - simError(); - } - calcNdf(); calcNdfRaw(); calcNdfTrans(); - - fortranInitialized_ = true; } - - std::set SimInfo::getUniqueAtomTypes() { + + set SimInfo::getSimulatedAtomTypes() { SimInfo::MoleculeIterator mi; Molecule* mol; Molecule::AtomIterator ai; Atom* atom; - std::set atomTypes; - - for(mol = beginMolecule(mi); mol != NULL; mol = nextMolecule(mi)) { - + set atomTypes; + + for(mol = beginMolecule(mi); mol != NULL; mol = nextMolecule(mi)) { for(atom = mol->beginAtom(ai); atom != NULL; atom = mol->nextAtom(ai)) { atomTypes.insert(atom->getAtomType()); + } + } + return atomTypes; + } + + /** + * setupCutoffs + * + * Sets the values of cutoffRadius and cutoffMethod + * + * cutoffRadius : realType + * If the cutoffRadius was explicitly set, use that value. + * If the cutoffRadius was not explicitly set: + * Are there electrostatic atoms? Use 12.0 Angstroms. + * No electrostatic atoms? Poll the atom types present in the + * simulation for suggested cutoff values (e.g. 2.5 * sigma). + * Use the maximum suggested value that was found. + * + * cutoffMethod : (one of HARD, SWITCHED, SHIFTED_FORCE, SHIFTED_POTENTIAL) + * If cutoffMethod was explicitly set, use that choice. + * If cutoffMethod was not explicitly set, use SHIFTED_FORCE + */ + void SimInfo::setupCutoffs() { + + if (simParams_->haveCutoffRadius()) { + cutoffRadius_ = simParams_->getCutoffRadius(); + } else { + if (usesElectrostaticAtoms_) { + sprintf(painCave.errMsg, + "SimInfo: No value was set for the cutoffRadius.\n" + "\tOpenMD will use a default value of 12.0 angstroms" + "\tfor the cutoffRadius.\n"); + painCave.isFatal = 0; + painCave.severity = OPENMD_INFO; + simError(); + cutoffRadius_ = 12.0; + } else { + RealType thisCut; + set::iterator i; + set atomTypes; + atomTypes = getSimulatedAtomTypes(); + for (i = atomTypes.begin(); i != atomTypes.end(); ++i) { + thisCut = InteractionManager::Instance()->getSuggestedCutoffRadius((*i)); + cutoffRadius_ = max(thisCut, cutoffRadius_); + } + sprintf(painCave.errMsg, + "SimInfo: No value was set for the cutoffRadius.\n" + "\tOpenMD will use %lf angstroms.\n", + cutoffRadius_); + painCave.isFatal = 0; + painCave.severity = OPENMD_INFO; + simError(); + } + } + + InteractionManager::Instance()->setCutoffRadius(cutoffRadius_); + + map stringToCutoffMethod; + stringToCutoffMethod["HARD"] = HARD; + stringToCutoffMethod["SWITCHING_FUNCTION"] = SWITCHING_FUNCTION; + stringToCutoffMethod["SHIFTED_POTENTIAL"] = SHIFTED_POTENTIAL; + stringToCutoffMethod["SHIFTED_FORCE"] = SHIFTED_FORCE; + + if (simParams_->haveCutoffMethod()) { + string cutMeth = toUpperCopy(simParams_->getCutoffMethod()); + map::iterator i; + i = stringToCutoffMethod.find(cutMeth); + if (i == stringToCutoffMethod.end()) { + sprintf(painCave.errMsg, + "SimInfo: Could not find chosen cutoffMethod %s\n" + "\tShould be one of: " + "HARD, SWITCHING_FUNCTION, SHIFTED_POTENTIAL, or SHIFTED_FORCE\n", + cutMeth.c_str()); + painCave.isFatal = 1; + painCave.severity = OPENMD_ERROR; + simError(); + } else { + cutoffMethod_ = i->second; } - + } else { + sprintf(painCave.errMsg, + "SimInfo: No value was set for the cutoffMethod.\n" + "\tOpenMD will use SHIFTED_FORCE.\n"); + painCave.isFatal = 0; + painCave.severity = OPENMD_INFO; + simError(); + cutoffMethod_ = SHIFTED_FORCE; } - return atomTypes; + InteractionManager::Instance()->setCutoffMethod(cutoffMethod_); } + + /** + * setupSwitching + * + * Sets the values of switchingRadius and + * If the switchingRadius was explicitly set, use that value (but check it) + * If the switchingRadius was not explicitly set: use 0.85 * cutoffRadius_ + */ + void SimInfo::setupSwitching() { + + if (simParams_->haveSwitchingRadius()) { + switchingRadius_ = simParams_->getSwitchingRadius(); + if (switchingRadius_ > cutoffRadius_) { + sprintf(painCave.errMsg, + "SimInfo: switchingRadius (%f) is larger than cutoffRadius(%f)\n", + switchingRadius_, cutoffRadius_); + painCave.isFatal = 1; + painCave.severity = OPENMD_ERROR; + simError(); + } + } else { + switchingRadius_ = 0.85 * cutoffRadius_; + sprintf(painCave.errMsg, + "SimInfo: 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", switchingRadius_); + painCave.isFatal = 0; + painCave.severity = OPENMD_WARNING; + simError(); + } + + InteractionManager::Instance()->setSwitchingRadius(switchingRadius_); - void SimInfo::setupSimType() { - std::set::iterator i; - std::set atomTypes; - atomTypes = getUniqueAtomTypes(); + SwitchingFunctionType ft; - int useLennardJones = 0; - int useElectrostatic = 0; - int useEAM = 0; - int useSC = 0; - int useCharge = 0; - int useDirectional = 0; - int useDipole = 0; - int useGayBerne = 0; - int useSticky = 0; - int useStickyPower = 0; - int useShape = 0; - int useFLARB = 0; //it is not in AtomType yet - int useDirectionalAtom = 0; - int useElectrostatics = 0; - //usePBC and useRF are from simParams - int usePBC = simParams_->getUsePeriodicBoundaryConditions(); - int useRF; - int useSF; - int useSP; - int useBoxDipole; - std::string myMethod; + if (simParams_->haveSwitchingFunctionType()) { + 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 : Unknown switchingFunctionType. (Input file specified %s .)\n" + "\tswitchingFunctionType must be one of: " + "\"cubic\" or \"fifth_order_polynomial\".", + funcType.c_str() ); + painCave.isFatal = 1; + painCave.severity = OPENMD_ERROR; + simError(); + } + } + } - // set the useRF logical - useRF = 0; - useSF = 0; - useSP = 0; + InteractionManager::Instance()->setSwitchingFunctionType(ft); + } + /** + * setupSkinThickness + * + * If the skinThickness was explicitly set, use that value (but check it) + * If the skinThickness was not explicitly set: use 1.0 angstroms + */ + void SimInfo::setupSkinThickness() { + if (simParams_->haveSkinThickness()) { + skinThickness_ = simParams_->getSkinThickness(); + } else { + skinThickness_ = 1.0; + sprintf(painCave.errMsg, + "SimInfo Warning: No value was set for the skinThickness.\n" + "\tOpenMD will use a default value of %f Angstroms\n" + "\tfor this simulation\n", skinThickness_); + painCave.isFatal = 0; + simError(); + } + } - 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; + void SimInfo::setupSimType() { + set::iterator i; + set atomTypes; + atomTypes = getSimulatedAtomTypes(); + useAtomicVirial_ = simParams_->getUseAtomicVirial(); + + int usesElectrostatic = 0; + int usesMetallic = 0; + int usesDirectional = 0; //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(); - useGayBerne |= (*i)->isGayBerne(); - useSticky |= (*i)->isSticky(); - useStickyPower |= (*i)->isStickyPower(); - useShape |= (*i)->isShape(); + usesElectrostatic |= (*i)->isElectrostatic(); + usesMetallic |= (*i)->isMetal(); + usesDirectional |= (*i)->isDirectional(); } - if (useSticky || useStickyPower || useDipole || useGayBerne || useShape) { - useDirectionalAtom = 1; - } - - if (useCharge || useDipole) { - useElectrostatics = 1; - } - #ifdef IS_MPI int temp; + temp = usesDirectional; + MPI_Allreduce(&temp, &usesDirectionalAtoms_, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD); - temp = usePBC; - MPI_Allreduce(&temp, &usePBC, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD); + temp = usesMetallic; + MPI_Allreduce(&temp, &usesMetallicAtoms_, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD); - temp = useDirectionalAtom; - MPI_Allreduce(&temp, &useDirectionalAtom, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD); - - temp = useLennardJones; - MPI_Allreduce(&temp, &useLennardJones, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD); - - temp = useElectrostatics; - MPI_Allreduce(&temp, &useElectrostatics, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD); - - temp = useCharge; - MPI_Allreduce(&temp, &useCharge, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD); - - temp = useDipole; - MPI_Allreduce(&temp, &useDipole, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD); - - temp = useSticky; - MPI_Allreduce(&temp, &useSticky, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD); - - temp = useStickyPower; - MPI_Allreduce(&temp, &useStickyPower, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD); - - temp = useGayBerne; - MPI_Allreduce(&temp, &useGayBerne, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD); - - 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 = usesElectrostatic; + MPI_Allreduce(&temp, &usesElectrostaticAtoms_, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD); #endif - - fInfo_.SIM_uses_PBC = usePBC; - fInfo_.SIM_uses_DirectionalAtoms = useDirectionalAtom; - fInfo_.SIM_uses_LennardJones = useLennardJones; - fInfo_.SIM_uses_Electrostatics = useElectrostatics; - fInfo_.SIM_uses_Charges = useCharge; - fInfo_.SIM_uses_Dipoles = useDipole; - fInfo_.SIM_uses_Sticky = useSticky; - 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; - fInfo_.SIM_uses_RF = useRF; - fInfo_.SIM_uses_SF = useSF; - fInfo_.SIM_uses_SP = useSP; - fInfo_.SIM_uses_BoxDipole = useBoxDipole; + fInfo_.SIM_uses_PBC = usesPeriodicBoundaries_; + fInfo_.SIM_uses_DirectionalAtoms = usesDirectionalAtoms_; + fInfo_.SIM_uses_MetallicAtoms = usesMetallicAtoms_; + fInfo_.SIM_requires_SkipCorrection = usesElectrostaticAtoms_; + fInfo_.SIM_requires_SelfCorrection = usesElectrostaticAtoms_; + fInfo_.SIM_uses_AtomicVirial = usesAtomicVirial_; } void SimInfo::setupFortranSim() { int isError; - int nExclude; - std::vector fortranGlobalGroupMembership; + int nExclude, nOneTwo, nOneThree, nOneFour; + vector fortranGlobalGroupMembership; - nExclude = exclude_.getSize(); + notifyFortranSkinThickness(&skinThickness_); + + int ljsp = cutoffMethod_ == SHIFTED_POTENTIAL ? 1 : 0; + int ljsf = cutoffMethod_ == SHIFTED_FORCE ? 1 : 0; + notifyFortranCutoffs(&cutoffRadius_, &switchingRadius_, &ljsp, &ljsf); + isError = 0; //globalGroupMembership_ is filled by SimCreator @@ -802,7 +916,7 @@ namespace oopse { } //calculate mass ratio of cutoff group - std::vector mfact; + vector mfact; SimInfo::MoleculeIterator mi; Molecule* mol; Molecule::CutoffGroupIterator ci; @@ -825,12 +939,11 @@ namespace oopse { else mfact.push_back( 1.0 ); } - } } //fill ident array of local atoms (it is actually ident of AtomType, it is so confusing !!!) - std::vector identArray; + vector identArray; //to avoid memory reallocation, reserve enough space identArray identArray.reserve(getNAtoms()); @@ -843,34 +956,46 @@ namespace oopse { //fill molMembershipArray //molMembershipArray is filled by SimCreator - std::vector molMembershipArray(nGlobalAtoms_); + vector molMembershipArray(nGlobalAtoms_); for (int i = 0; i < nGlobalAtoms_; i++) { molMembershipArray[i] = globalMolMembership_[i] + 1; } //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); + + nExclude = excludedInteractions_.getSize(); + nOneTwo = oneTwoInteractions_.getSize(); + nOneThree = oneThreeInteractions_.getSize(); + nOneFour = oneFourInteractions_.getSize(); - if( isError ){ + 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(); @@ -881,12 +1006,11 @@ namespace oopse { } -#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; + vector localToGlobalAtomIndex(getNAtoms(), 0); + vector localToGlobalCutoffGroupIndex; SimInfo::MoleculeIterator mi; Molecule::AtomIterator ai; Molecule::CutoffGroupIterator ci; @@ -933,259 +1057,14 @@ namespace oopse { } sprintf(checkPointMsg, " mpiRefresh successful.\n"); - MPIcheckPoint(); + errorCheckPoint(); - - } - #endif - - void SimInfo::setupCutoff() { - - ForceFieldOptions& forceFieldOptions_ = forceField_->getForceFieldOptions(); - - // Check the cutoff policy - int cp = TRADITIONAL_CUTOFF_POLICY; // Set to traditional by default - - std::string myPolicy; - if (forceFieldOptions_.haveCutoffPolicy()){ - myPolicy = forceFieldOptions_.getCutoffPolicy(); - }else if (simParams_->haveCutoffPolicy()) { - myPolicy = simParams_->getCutoffPolicy(); - } - - if (!myPolicy.empty()){ - toUpper(myPolicy); - if (myPolicy == "MIX") { - cp = MIX_CUTOFF_POLICY; - } else { - if (myPolicy == "MAX") { - cp = MAX_CUTOFF_POLICY; - } else { - if (myPolicy == "TRADITIONAL") { - cp = TRADITIONAL_CUTOFF_POLICY; - } else { - // throw error - sprintf( painCave.errMsg, - "SimInfo error: Unknown cutoffPolicy. (Input file specified %s .)\n\tcutoffPolicy must be one of: \"Mix\", \"Max\", or \"Traditional\".", myPolicy.c_str() ); - painCave.isFatal = 1; - simError(); - } - } - } - } - 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" - "\tOOPSE 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" - "\tOOPSE will use the same value as the cutoffRadius.\n" - "\tswitchingRadius = %f. for this simulation\n", rsw_); - painCave.isFatal = 0; - simError(); - } - } - - notifyFortranCutoffs(&rcut_, &rsw_); - - } 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" - "\tOOPSE 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); - 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" - "\tOOPSE will use a default value of\n" - "\t0.85 * cutoffRadius for the switchingRadius\n"); - painCave.isFatal = 0; - simError(); - rsw_ = 0.85 * rcut_; - } - notifyFortranCutoffs(&rcut_, &rsw_); - } 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(); - - } - } } - void SimInfo::setupElectrostaticSummationMethod( int isError ) { - - int errorOut; - int esm = NONE; - int 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 == "SWITCHING_FUNCTION") { - esm = SWITCHING_FUNCTION; - } else { - if (myMethod == "SHIFTED_POTENTIAL") { - esm = SHIFTED_POTENTIAL; - } else { - if (myMethod == "SHIFTED_FORCE") { - esm = SHIFTED_FORCE; - } else { - 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(); - } - } - } - } - } - } - - 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(); - } - } - } - - // let's pass some summation method variables to fortran - setElectrostaticSummationMethod( &esm ); - setFortranElectrostaticMethod( &esm ); - setScreeningMethod( &sm ); - setDampingAlpha( &alphaVal ); - 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() { @@ -1193,7 +1072,6 @@ namespace oopse { // we only call setAccumulateBoxDipole if the accumulateBoxDipole parameter is true if ( simParams_->haveAccumulateBoxDipole() ) if ( simParams_->getAccumulateBoxDipole() ) { - setAccumulateBoxDipole(); calcBoxDipole_ = true; } @@ -1203,7 +1081,7 @@ namespace oopse { properties_.addProperty(genData); } - void SimInfo::removeProperty(const std::string& propName) { + void SimInfo::removeProperty(const string& propName) { properties_.removeProperty(propName); } @@ -1211,15 +1089,15 @@ namespace oopse { properties_.clearProperties(); } - std::vector SimInfo::getPropertyNames() { + vector SimInfo::getPropertyNames() { return properties_.getPropertyNames(); } - std::vector SimInfo::getProperties() { + vector SimInfo::getProperties() { return properties_.getProperties(); } - GenericData* SimInfo::getPropertyByName(const std::string& propName) { + GenericData* SimInfo::getPropertyByName(const string& propName) { return properties_.getPropertyByName(propName); } @@ -1302,7 +1180,7 @@ namespace oopse { } - std::ostream& operator <<(std::ostream& o, SimInfo& info) { + ostream& operator <<(ostream& o, SimInfo& info) { return o; } @@ -1345,7 +1223,7 @@ namespace oopse { [ Ixx -Ixy -Ixz ] - J =| -Iyx Iyy -Iyz | + J =| -Iyx Iyy -Iyz | [ -Izx -Iyz Izz ] */ @@ -1452,7 +1330,7 @@ namespace oopse { return IOIndexToIntegrableObject.at(index); } - void SimInfo::setIOIndexToIntegrableObject(const std::vector& v) { + void SimInfo::setIOIndexToIntegrableObject(const vector& v) { IOIndexToIntegrableObject= v; } @@ -1494,7 +1372,7 @@ namespace oopse { return; } /* - void SimInfo::setStuntDoubleFromGlobalIndex(std::vector v) { + void SimInfo::setStuntDoubleFromGlobalIndex(vector v) { assert( v.size() == nAtoms_ + nRigidBodies_); sdByGlobalIndex_ = v; } @@ -1504,5 +1382,16 @@ namespace oopse { return sdByGlobalIndex_.at(index); } */ -}//end namespace oopse + int SimInfo::getNGlobalConstraints() { + int nGlobalConstraints; +#ifdef IS_MPI + MPI_Allreduce(&nConstraints_, &nGlobalConstraints, 1, MPI_INT, MPI_SUM, + MPI_COMM_WORLD); +#else + nGlobalConstraints = nConstraints_; +#endif + return nGlobalConstraints; + } +}//end namespace OpenMD +