--- trunk/src/brains/SimInfo.cpp 2005/11/11 15:22:11 726 +++ branches/development/src/brains/SimInfo.cpp 2013/02/20 13:52:51 1849 @@ -6,19 +6,10 @@ * redistribute this software in source and binary code form, provided * that the following conditions are met: * - * 1. Acknowledgement of the program authors must be made in any - * publication of scientific results based in part on use of the - * program. An acceptable form of acknowledgement is citation of - * the article in which the program was described (Matthew - * A. Meineke, Charles F. Vardeman II, Teng Lin, Christopher - * J. Fennell and J. Daniel Gezelter, "OOPSE: An Object-Oriented - * Parallel Simulation Engine for Molecular Dynamics," - * J. Comput. Chem. 26, pp. 252-271 (2005)) - * - * 2. Redistributions of source code must retain the above copyright + * 1. Redistributions of source code must retain the above copyright * notice, this list of conditions and the following disclaimer. * - * 3. Redistributions in binary form must reproduce the above copyright + * 2. Redistributions in binary form must reproduce the above copyright * notice, this list of conditions and the following disclaimer in the * documentation and/or other materials provided with the * distribution. @@ -37,6 +28,16 @@ * arising out of the use of or inability to use software, even if the * University of Notre Dame has been advised of the possibility of * such damages. + * + * SUPPORT OPEN SCIENCE! If you use OpenMD or its source code in your + * research, please cite the appropriate papers when you publish your + * work. Good starting points are: + * + * [1] Meineke, et al., J. Comp. Chem. 26, 252-271 (2005). + * [2] Fennell & Gezelter, J. Chem. Phys. 124, 234104 (2006). + * [3] Sun, Lin & Gezelter, J. Chem. Phys. 128, 24107 (2008). + * [4] Kuang & Gezelter, J. Chem. Phys. 133, 164101 (2010). + * [5] Vardeman, Stocker & Gezelter, J. Chem. Theory Comput. 7, 834 (2011). */ /** @@ -48,158 +49,141 @@ #include #include +#include #include "brains/SimInfo.hpp" #include "math/Vector3.hpp" #include "primitives/Molecule.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/electrostatic_interface.h" -#include "UseTheForce/notifyCutoffs_interface.h" -#include "UseTheForce/DarkSide/switcheroo_interface.h" +#include "primitives/StuntDouble.hpp" #include "utils/MemoryUtils.hpp" #include "utils/simError.h" #include "selection/SelectionManager.hpp" - +#include "io/ForceFieldOptions.hpp" +#include "brains/ForceField.hpp" +#include "nonbonded/SwitchingFunction.hpp" #ifdef IS_MPI -#include "UseTheForce/mpiComponentPlan.h" -#include "UseTheForce/DarkSide/simParallel_interface.h" -#endif +#include +#endif -namespace oopse { - - 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), +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) { - - - 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 - CutoffGroupStamp* cgStamp; - RigidBodyStamp* rbStamp; - int nRigidAtoms = 0; + nGlobalIntegrableObjects_(0), nGlobalRigidBodies_(0), nGlobalFluctuatingCharges_(0), + nAtoms_(0), nBonds_(0), nBends_(0), nTorsions_(0), nInversions_(0), + nRigidBodies_(0), nIntegrableObjects_(0), nCutoffGroups_(0), + nConstraints_(0), nFluctuatingCharges_(0), sman_(NULL), topologyDone_(false), + calcBoxDipole_(false), useAtomicVirial_(true) { - for (i = molStampPairs.begin(); i !=molStampPairs.end(); ++i) { - molStamp = i->first; - nMolWithSameStamp = i->second; - - 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); - 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->getRigidBody(j); - nAtomsInRigidBodies += rbStamp->getNMembers(); - } - - nGlobalRigidBodies_ += nRigidBodiesInStamp * nMolWithSameStamp; - nRigidAtoms += nAtomsInRigidBodies * nMolWithSameStamp; - + 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(); + + 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; } molecules_.clear(); - delete stamps_; delete sman_; delete simParams_; 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()); @@ -212,12 +196,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; @@ -226,8 +211,6 @@ namespace oopse { } else { return false; } - - } @@ -243,38 +226,54 @@ namespace oopse { void SimInfo::calcNdf() { - int ndf_local; + int ndf_local, nfq_local; MoleculeIterator i; - std::vector::iterator j; + vector::iterator j; + vector::iterator k; + Molecule* mol; - StuntDouble* integrableObject; + StuntDouble* sd; + Atom* atom; ndf_local = 0; + nfq_local = 0; for (mol = beginMolecule(i); mol != NULL; mol = nextMolecule(i)) { - for (integrableObject = mol->beginIntegrableObject(j); integrableObject != NULL; - integrableObject = mol->nextIntegrableObject(j)) { + for (sd = mol->beginIntegrableObject(j); sd != NULL; + sd = mol->nextIntegrableObject(j)) { + ndf_local += 3; - if (integrableObject->isDirectional()) { - if (integrableObject->isLinear()) { + if (sd->isDirectional()) { + if (sd->isLinear()) { ndf_local += 2; } else { ndf_local += 3; } } - - }//end for (integrableObject) - }// end for (mol) + } + + for (atom = mol->beginFluctuatingCharge(k); atom != NULL; + atom = mol->nextFluctuatingCharge(k)) { + if (atom->isFluctuatingCharge()) { + nfq_local++; + } + } + } + ndfLocal_ = ndf_local; + // n_constraints is local, so subtract them on each processor ndf_local -= nConstraints_; #ifdef IS_MPI - MPI_Allreduce(&ndf_local,&ndf_,1,MPI_INT,MPI_SUM, MPI_COMM_WORLD); + MPI::COMM_WORLD.Allreduce(&ndf_local, &ndf_, 1, MPI::INT,MPI::SUM); + MPI::COMM_WORLD.Allreduce(&nfq_local, &nGlobalFluctuatingCharges_, 1, + MPI::INT, MPI::SUM); #else ndf_ = ndf_local; + nGlobalFluctuatingCharges_ = nfq_local; #endif // nZconstraints_ is global, as are the 3 COM translations for the @@ -283,25 +282,54 @@ namespace oopse { } + int SimInfo::getFdf() { +#ifdef IS_MPI + MPI::COMM_WORLD.Allreduce(&fdf_local, &fdf_, 1, MPI::INT, MPI::SUM); +#else + fdf_ = fdf_local; +#endif + return fdf_; + } + + unsigned int SimInfo::getNLocalCutoffGroups(){ + int nLocalCutoffAtoms = 0; + Molecule* mol; + MoleculeIterator mi; + CutoffGroup* cg; + Molecule::CutoffGroupIterator ci; + + for (mol = beginMolecule(mi); mol != NULL; mol = nextMolecule(mi)) { + + for (cg = mol->beginCutoffGroup(ci); cg != NULL; + cg = mol->nextCutoffGroup(ci)) { + nLocalCutoffAtoms += cg->getNumAtom(); + + } + } + + return nAtoms_ - nLocalCutoffAtoms + nCutoffGroups_; + } + void SimInfo::calcNdfRaw() { int ndfRaw_local; MoleculeIterator i; - std::vector::iterator j; + vector::iterator j; Molecule* mol; - StuntDouble* integrableObject; + StuntDouble* sd; // Raw degrees of freedom that we have to set ndfRaw_local = 0; for (mol = beginMolecule(i); mol != NULL; mol = nextMolecule(i)) { - for (integrableObject = mol->beginIntegrableObject(j); integrableObject != NULL; - integrableObject = mol->nextIntegrableObject(j)) { + for (sd = mol->beginIntegrableObject(j); sd != NULL; + sd = mol->nextIntegrableObject(j)) { + ndfRaw_local += 3; - if (integrableObject->isDirectional()) { - if (integrableObject->isLinear()) { + if (sd->isDirectional()) { + if (sd->isLinear()) { ndfRaw_local += 2; } else { ndfRaw_local += 3; @@ -312,7 +340,7 @@ namespace oopse { } #ifdef IS_MPI - MPI_Allreduce(&ndfRaw_local,&ndfRaw_,1,MPI_INT,MPI_SUM, MPI_COMM_WORLD); + MPI::COMM_WORLD.Allreduce(&ndfRaw_local, &ndfRaw_, 1, MPI::INT, MPI::SUM); #else ndfRaw_ = ndfRaw_local; #endif @@ -325,7 +353,8 @@ namespace oopse { #ifdef IS_MPI - MPI_Allreduce(&ndfTrans_local,&ndfTrans_,1,MPI_INT,MPI_SUM, MPI_COMM_WORLD); + MPI::COMM_WORLD.Allreduce(&ndfTrans_local, &ndfTrans_, 1, + MPI::INT, MPI::SUM); #else ndfTrans_ = ndfTrans_local; #endif @@ -334,124 +363,322 @@ 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; + + // 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* sd; - for (bond= mol->beginBond(bondIter); bond != NULL; bond = mol->nextBond(bondIter)) { + for (sd = mol->beginIntegrableObject(ii); sd != NULL; + sd = mol->nextIntegrableObject(ii)) { + + if (sd->isRigidBody()) { + rb = static_cast(sd); + 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 { + set oneAtomSet; + oneAtomSet.insert(sd->getGlobalIndex()); + atomGroups.insert(map >::value_type(sd->getGlobalIndex(), oneAtomSet)); + } + } + + 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(); + + 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)) { - 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; + + map > atomGroups; + Molecule::RigidBodyIterator rbIter; + RigidBody* rb; + Molecule::IntegrableObjectIterator ii; + StuntDouble* sd; - for (bond= mol->beginBond(bondIter); bond != NULL; bond = mol->nextBond(bondIter)) { + for (sd = mol->beginIntegrableObject(ii); sd != NULL; + sd = mol->nextIntegrableObject(ii)) { + + if (sd->isRigidBody()) { + rb = static_cast(sd); + 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 { + set oneAtomSet; + oneAtomSet.insert(sd->getGlobalIndex()); + atomGroups.insert(map >::value_type(sd->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)) { - 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_.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)) { + 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(); @@ -459,559 +686,276 @@ namespace oopse { molStampIds_.insert(molStampIds_.end(), nmol, curStampId); } - void SimInfo::update() { - setupSimType(); - -#ifdef IS_MPI - setupFortranParallel(); -#endif - - setupFortranSim(); - - //setup fortran force field - /** @deprecate */ - int isError = 0; - - setupElectrostaticSummationMethod( isError ); - setupSwitchingFunction(); - - if(isError){ - sprintf( painCave.errMsg, - "ForceField error: There was an error initializing the forceField in fortran.\n" ); - painCave.isFatal = 1; - simError(); - } - - - setupCutoff(); - + /** + * update + * + * Performs the global checks and variable settings after the + * objects have been created. + * + */ + void SimInfo::update() { + setupSimVariables(); calcNdf(); calcNdfRaw(); calcNdfTrans(); - - fortranInitialized_ = true; } - - std::set SimInfo::getUniqueAtomTypes() { + + /** + * getSimulatedAtomTypes + * + * Returns an STL set of AtomType* that are actually present in this + * simulation. Must query all processors to assemble this information. + * + */ + set SimInfo::getSimulatedAtomTypes() { SimInfo::MoleculeIterator mi; Molecule* mol; Molecule::AtomIterator ai; Atom* atom; - std::set atomTypes; - + set atomTypes; + for(mol = beginMolecule(mi); mol != NULL; mol = nextMolecule(mi)) { - - for(atom = mol->beginAtom(ai); atom != NULL; atom = mol->nextAtom(ai)) { + for(atom = mol->beginAtom(ai); atom != NULL; + atom = mol->nextAtom(ai)) { atomTypes.insert(atom->getAtomType()); - } - + } + } + +#ifdef IS_MPI + + // loop over the found atom types on this processor, and add their + // numerical idents to a vector: + + vector foundTypes; + set::iterator i; + for (i = atomTypes.begin(); i != atomTypes.end(); ++i) + foundTypes.push_back( (*i)->getIdent() ); + + // count_local holds the number of found types on this processor + int count_local = foundTypes.size(); + + int nproc = MPI::COMM_WORLD.Get_size(); + + // we need arrays to hold the counts and displacement vectors for + // all processors + vector counts(nproc, 0); + vector disps(nproc, 0); + + // fill the counts array + MPI::COMM_WORLD.Allgather(&count_local, 1, MPI::INT, &counts[0], + 1, MPI::INT); + + // use the processor counts to compute the displacement array + disps[0] = 0; + int totalCount = counts[0]; + for (int iproc = 1; iproc < nproc; iproc++) { + disps[iproc] = disps[iproc-1] + counts[iproc-1]; + totalCount += counts[iproc]; } + // we need a (possibly redundant) set of all found types: + vector ftGlobal(totalCount); + + // now spray out the foundTypes to all the other processors: + MPI::COMM_WORLD.Allgatherv(&foundTypes[0], count_local, MPI::INT, + &ftGlobal[0], &counts[0], &disps[0], + MPI::INT); + + vector::iterator j; + + // foundIdents is a stl set, so inserting an already found ident + // will have no effect. + set foundIdents; + + for (j = ftGlobal.begin(); j != ftGlobal.end(); ++j) + foundIdents.insert((*j)); + + // now iterate over the foundIdents and get the actual atom types + // that correspond to these: + set::iterator it; + for (it = foundIdents.begin(); it != foundIdents.end(); ++it) + atomTypes.insert( forceField_->getAtomType((*it)) ); + +#endif + return atomTypes; } - void SimInfo::setupSimType() { - std::set::iterator i; - std::set atomTypes; - atomTypes = getUniqueAtomTypes(); - - int useLennardJones = 0; - int useElectrostatic = 0; - int useEAM = 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; - std::string myMethod; - // set the useRF logical - useRF = 0; - useSF = 0; + int getGlobalCountOfType(AtomType* atype) { + /* + set atypes = getSimulatedAtomTypes(); + map counts_; + for(mol = beginMolecule(mi); mol != NULL; mol = nextMolecule(mi)) { + for(atom = mol->beginAtom(ai); atom != NULL; + atom = mol->nextAtom(ai)) { + atom->getAtomType(); + } + } + */ + return 0; + } - if (simParams_->haveElectrostaticSummationMethod()) { - std::string myMethod = simParams_->getElectrostaticSummationMethod(); - toUpper(myMethod); - if (myMethod == "REACTION_FIELD") { - useRF=1; - } else { - if (myMethod == "SHIFTED_FORCE") { - useSF = 1; - } + void SimInfo::setupSimVariables() { + useAtomicVirial_ = simParams_->getUseAtomicVirial(); + // we only call setAccumulateBoxDipole if the accumulateBoxDipole + // parameter is true + calcBoxDipole_ = false; + if ( simParams_->haveAccumulateBoxDipole() ) + if ( simParams_->getAccumulateBoxDipole() ) { + calcBoxDipole_ = true; } - } - + + set::iterator i; + set atomTypes; + atomTypes = getSimulatedAtomTypes(); + bool usesElectrostatic = false; + bool usesMetallic = false; + bool usesDirectional = false; + bool usesFluctuatingCharges = false; //loop over all of the atom types for (i = atomTypes.begin(); i != atomTypes.end(); ++i) { - useLennardJones |= (*i)->isLennardJones(); - useElectrostatic |= (*i)->isElectrostatic(); - useEAM |= (*i)->isEAM(); - 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(); + usesFluctuatingCharges |= (*i)->isFluctuatingCharge(); } - if (useSticky || useStickyPower || useDipole || useGayBerne || useShape) { - useDirectionalAtom = 1; - } +#ifdef IS_MPI + bool temp; + temp = usesDirectional; + MPI::COMM_WORLD.Allreduce(&temp, &usesDirectionalAtoms_, 1, MPI::BOOL, + MPI::LOR); + + temp = usesMetallic; + MPI::COMM_WORLD.Allreduce(&temp, &usesMetallicAtoms_, 1, MPI::BOOL, + MPI::LOR); + + temp = usesElectrostatic; + MPI::COMM_WORLD.Allreduce(&temp, &usesElectrostaticAtoms_, 1, MPI::BOOL, + MPI::LOR); - if (useCharge || useDipole) { - useElectrostatics = 1; - } + temp = usesFluctuatingCharges; + MPI::COMM_WORLD.Allreduce(&temp, &usesFluctuatingCharges_, 1, MPI::BOOL, + MPI::LOR); +#else -#ifdef IS_MPI - int temp; + usesDirectionalAtoms_ = usesDirectional; + usesMetallicAtoms_ = usesMetallic; + usesElectrostaticAtoms_ = usesElectrostatic; + usesFluctuatingCharges_ = usesFluctuatingCharges; - temp = usePBC; - MPI_Allreduce(&temp, &usePBC, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD); +#endif + + requiresPrepair_ = usesMetallicAtoms_ ? true : false; + requiresSkipCorrection_ = usesElectrostaticAtoms_ ? true : false; + requiresSelfCorrection_ = usesElectrostaticAtoms_ ? true : false; + } - 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); + vector SimInfo::getGlobalAtomIndices() { + SimInfo::MoleculeIterator mi; + Molecule* mol; + Molecule::AtomIterator ai; + Atom* atom; - 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); + vector GlobalAtomIndices(getNAtoms(), 0); - temp = useGayBerne; - MPI_Allreduce(&temp, &useGayBerne, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD); + for (mol = beginMolecule(mi); mol != NULL; mol = nextMolecule(mi)) { + + for (atom = mol->beginAtom(ai); atom != NULL; atom = mol->nextAtom(ai)) { + GlobalAtomIndices[atom->getLocalIndex()] = atom->getGlobalIndex(); + } + } + return GlobalAtomIndices; + } - temp = useEAM; - MPI_Allreduce(&temp, &useEAM, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD); - temp = useShape; - MPI_Allreduce(&temp, &useShape, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD); + vector SimInfo::getGlobalGroupIndices() { + SimInfo::MoleculeIterator mi; + Molecule* mol; + Molecule::CutoffGroupIterator ci; + CutoffGroup* cg; - 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); - -#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_Shapes = useShape; - fInfo_.SIM_uses_FLARB = useFLARB; - fInfo_.SIM_uses_RF = useRF; - fInfo_.SIM_uses_SF = useSF; - - if( myMethod == "REACTION_FIELD") { + vector GlobalGroupIndices; + + for (mol = beginMolecule(mi); mol != NULL; mol = nextMolecule(mi)) { - 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(); - } + //local index of cutoff group is trivial, it only depends on the + //order of travesing + for (cg = mol->beginCutoffGroup(ci); cg != NULL; + cg = mol->nextCutoffGroup(ci)) { + GlobalGroupIndices.push_back(cg->getGlobalIndex()); + } } - + return GlobalGroupIndices; } - void SimInfo::setupFortranSim() { - int isError; - int nExclude; - std::vector fortranGlobalGroupMembership; - - nExclude = exclude_.getSize(); - isError = 0; - //globalGroupMembership_ is filled by SimCreator - for (int i = 0; i < nGlobalAtoms_; i++) { - fortranGlobalGroupMembership.push_back(globalGroupMembership_[i] + 1); - } + void SimInfo::prepareTopology() { //calculate mass ratio of cutoff group - 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()); + /** + * The mass factor is the relative mass of an atom to the total + * mass of the cutoff group it belongs to. By default, all atoms + * are their own cutoff groups, and therefore have mass factors of + * 1. We need some special handling for massless atoms, which + * will be treated as carrying the entire mass of the cutoff + * group. + */ + massFactors_.clear(); + massFactors_.resize(getNAtoms(), 1.0); for(mol = beginMolecule(mi); mol != NULL; mol = nextMolecule(mi)) { - for (cg = mol->beginCutoffGroup(ci); cg != NULL; cg = mol->nextCutoffGroup(ci)) { + for (cg = mol->beginCutoffGroup(ci); cg != NULL; + cg = mol->nextCutoffGroup(ci)) { totalMass = cg->getMass(); for(atom = cg->beginAtom(ai); atom != NULL; atom = cg->nextAtom(ai)) { // Check for massless groups - set mfact to 1 if true - if (totalMass != 0) - mfact.push_back(atom->getMass()/totalMass); + if (totalMass != 0) + massFactors_[atom->getLocalIndex()] = atom->getMass()/totalMass; else - mfact.push_back( 1.0 ); + massFactors_[atom->getLocalIndex()] = 1.0; } - } } - //fill ident array of local atoms (it is actually ident of AtomType, it is so confusing !!!) - std::vector identArray; + // Build the identArray_ - //to avoid memory reallocation, reserve enough space identArray - identArray.reserve(getNAtoms()); - + identArray_.clear(); + identArray_.reserve(getNAtoms()); for(mol = beginMolecule(mi); mol != NULL; mol = nextMolecule(mi)) { for(atom = mol->beginAtom(ai); atom != NULL; atom = mol->nextAtom(ai)) { - identArray.push_back(atom->getIdent()); + identArray_.push_back(atom->getIdent()); } } - - //fill molMembershipArray - //molMembershipArray is filled by SimCreator - std::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); + //scan topology - if( isError ){ - - sprintf( painCave.errMsg, - "There was an error setting the simulation information in fortran.\n" ); - painCave.isFatal = 1; - painCave.severity = OOPSE_ERROR; - simError(); - } - -#ifdef IS_MPI - sprintf( checkPointMsg, - "succesfully sent the simulation information to fortran.\n"); - MPIcheckPoint(); -#endif // is_mpi - } - - -#ifdef IS_MPI - void SimInfo::setupFortranParallel() { - - //SimInfo is responsible for creating localToGlobalAtomIndex and localToGlobalGroupIndex - std::vector localToGlobalAtomIndex(getNAtoms(), 0); - std::vector localToGlobalCutoffGroupIndex; - SimInfo::MoleculeIterator mi; - Molecule::AtomIterator ai; - Molecule::CutoffGroupIterator ci; - Molecule* mol; - Atom* atom; - CutoffGroup* cg; - mpiSimData parallelData; - int isError; - - for (mol = beginMolecule(mi); mol != NULL; mol = nextMolecule(mi)) { - - //local index(index in DataStorge) of atom is important - for (atom = mol->beginAtom(ai); atom != NULL; atom = mol->nextAtom(ai)) { - localToGlobalAtomIndex[atom->getLocalIndex()] = atom->getGlobalIndex() + 1; - } - - //local index of cutoff group is trivial, it only depends on the order of travesing - for (cg = mol->beginCutoffGroup(ci); cg != NULL; cg = mol->nextCutoffGroup(ci)) { - localToGlobalCutoffGroupIndex.push_back(cg->getGlobalIndex() + 1); - } - - } - - //fill up mpiSimData struct - parallelData.nMolGlobal = getNGlobalMolecules(); - parallelData.nMolLocal = getNMolecules(); - parallelData.nAtomsGlobal = getNGlobalAtoms(); - parallelData.nAtomsLocal = getNAtoms(); - parallelData.nGroupsGlobal = getNGlobalCutoffGroups(); - parallelData.nGroupsLocal = getNCutoffGroups(); - parallelData.myNode = worldRank; - MPI_Comm_size(MPI_COMM_WORLD, &(parallelData.nProcessors)); - - //pass mpiSimData struct and index arrays to fortran - setFsimParallel(¶llelData, &(parallelData.nAtomsLocal), - &localToGlobalAtomIndex[0], &(parallelData.nGroupsLocal), - &localToGlobalCutoffGroupIndex[0], &isError); - - if (isError) { - sprintf(painCave.errMsg, - "mpiRefresh errror: fortran didn't like something we gave it.\n"); - painCave.isFatal = 1; - simError(); - } - - sprintf(checkPointMsg, " mpiRefresh successful.\n"); - MPIcheckPoint(); - - - } - -#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) { - - if (fInfo_.SIM_uses_Charges | fInfo_.SIM_uses_Dipoles | fInfo_.SIM_uses_RF) { - - if (!simParams_->haveCutoffRadius()){ - 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_->getCutoffRadius(); - } - - if (!simParams_->haveSwitchingRadius()){ - 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; - } else{ - rsw = simParams_->getSwitchingRadius(); - } - - } 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_->haveCutoffRadius()) { - rcut = simParams_->getCutoffRadius(); - } else { - //set cutoff radius to the maximum cutoff radius based on atom types in the whole system - rcut = calcMaxCutoffRadius(); - } - - if (simParams_->haveSwitchingRadius()) { - rsw = simParams_->getSwitchingRadius(); - } else { - rsw = rcut; - } - - } - } - - 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(); - 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(); - } - } - } - } + int* excludeList = excludedInteractions_.getPairList(); + int* oneTwoList = oneTwoInteractions_.getPairList(); + int* oneThreeList = oneThreeInteractions_.getPairList(); + int* oneFourList = oneFourInteractions_.getPairList(); - - if (simParams_->haveSkinThickness()) { - double skinThickness = simParams_->getSkinThickness(); - } - - notifyFortranCutoffs(&rcut_, &rsw_, &rnblist, &cp); - // also send cutoff notification to electrostatics - setElectrostaticCutoffRadius(&rcut_, &rsw_); + topologyDone_ = true; } - void SimInfo::setupElectrostaticSummationMethod( int isError ) { - - int errorOut; - int esm = NONE; - int sm = UNDAMPED; - double alphaVal; - double dielectric; - - errorOut = isError; - alphaVal = simParams_->getDampingAlpha(); - dielectric = simParams_->getDielectric(); - - 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; - } else { - // throw error - sprintf( painCave.errMsg, - "SimInfo error: Unknown electrostaticSummationMethod. (Input file specified %s .)\n\telectrostaticSummationMethod must be one of: \"none\", \"shifted_potential\", \"shifted_force\", or \"reaction_field\".", 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()) { - //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.", alphaVal); - painCave.isFatal = 0; - simError(); - } - } else { - // throw error - sprintf( painCave.errMsg, - "SimInfo error: Unknown electrostaticScreeningMethod. (Input file specified %s .)\n\telectrostaticScreeningMethod must be one of: \"undamped\" or \"damped\".", myScreen.c_str() ); - painCave.isFatal = 1; - simError(); - } - } - } - - // let's pass some summation method variables to fortran - setElectrostaticSummationMethod( &esm ); - setScreeningMethod( &sm ); - setDampingAlpha( &alphaVal ); - setReactionFieldDielectric( &dielectric ); - initFortranFF( &esm, &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::addProperty(GenericData* genData) { properties_.addProperty(genData); } - void SimInfo::removeProperty(const std::string& propName) { + void SimInfo::removeProperty(const string& propName) { properties_.removeProperty(propName); } @@ -1019,15 +963,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); } @@ -1041,221 +985,65 @@ namespace oopse { Molecule* mol; RigidBody* rb; Atom* atom; + CutoffGroup* cg; SimInfo::MoleculeIterator mi; Molecule::RigidBodyIterator rbIter; - Molecule::AtomIterator atomIter;; + Molecule::AtomIterator atomIter; + Molecule::CutoffGroupIterator cgIter; for (mol = beginMolecule(mi); mol != NULL; mol = nextMolecule(mi)) { - for (atom = mol->beginAtom(atomIter); atom != NULL; atom = mol->nextAtom(atomIter)) { + for (atom = mol->beginAtom(atomIter); atom != NULL; + atom = mol->nextAtom(atomIter)) { atom->setSnapshotManager(sman_); } - for (rb = mol->beginRigidBody(rbIter); rb != NULL; rb = mol->nextRigidBody(rbIter)) { + for (rb = mol->beginRigidBody(rbIter); rb != NULL; + rb = mol->nextRigidBody(rbIter)) { rb->setSnapshotManager(sman_); } + + for (cg = mol->beginCutoffGroup(cgIter); cg != NULL; + cg = mol->nextCutoffGroup(cgIter)) { + cg->setSnapshotManager(sman_); + } } } - Vector3d SimInfo::getComVel(){ - SimInfo::MoleculeIterator i; - Molecule* mol; - Vector3d comVel(0.0); - double totalMass = 0.0; - - - for (mol = beginMolecule(i); mol != NULL; mol = nextMolecule(i)) { - double mass = mol->getMass(); - totalMass += mass; - comVel += mass * mol->getComVel(); - } + ostream& operator <<(ostream& o, SimInfo& info) { -#ifdef IS_MPI - double 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); -#endif - - comVel /= totalMass; - - return comVel; - } - - Vector3d SimInfo::getCom(){ - SimInfo::MoleculeIterator i; - Molecule* mol; - - Vector3d com(0.0); - double totalMass = 0.0; - - for (mol = beginMolecule(i); mol != NULL; mol = nextMolecule(i)) { - double mass = mol->getMass(); - totalMass += mass; - com += mass * mol->getCom(); - } - -#ifdef IS_MPI - double 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); -#endif - - com /= totalMass; - - return com; - - } - - std::ostream& operator <<(std::ostream& o, SimInfo& info) { - return o; } - - /* - Returns center of mass and center of mass velocity in one function call. - */ - - void SimInfo::getComAll(Vector3d &com, Vector3d &comVel){ - SimInfo::MoleculeIterator i; - Molecule* mol; - - - double totalMass = 0.0; - + + StuntDouble* SimInfo::getIOIndexToIntegrableObject(int index) { + if (index >= int(IOIndexToIntegrableObject.size())) { + sprintf(painCave.errMsg, + "SimInfo::getIOIndexToIntegrableObject Error: Integrable Object\n" + "\tindex exceeds number of known objects!\n"); + painCave.isFatal = 1; + simError(); + return NULL; + } else + return IOIndexToIntegrableObject.at(index); + } + + void SimInfo::setIOIndexToIntegrableObject(const vector& v) { + IOIndexToIntegrableObject= v; + } - for (mol = beginMolecule(i); mol != NULL; mol = nextMolecule(i)) { - double mass = mol->getMass(); - totalMass += mass; - com += mass * mol->getCom(); - comVel += mass * mol->getComVel(); - } - + int SimInfo::getNGlobalConstraints() { + int nGlobalConstraints; #ifdef IS_MPI - double 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::COMM_WORLD.Allreduce(&nConstraints_, &nGlobalConstraints, 1, + MPI::INT, MPI::SUM); +#else + nGlobalConstraints = nConstraints_; #endif - - com /= totalMass; - comVel /= totalMass; - } - - /* - Return intertia tensor for entire system and angular momentum Vector. + return nGlobalConstraints; + } +}//end namespace OpenMD - [ Ixx -Ixy -Ixz ] - J =| -Iyx Iyy -Iyz | - [ -Izx -Iyz Izz ] - */ - - 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; - Vector3d com(0.0); - Vector3d comVel(0.0); - - getComAll(com, comVel); - - SimInfo::MoleculeIterator i; - Molecule* mol; - - Vector3d thisq(0.0); - Vector3d thisv(0.0); - - double thisMass = 0.0; - - - - - for (mol = beginMolecule(i); mol != NULL; mol = nextMolecule(i)) { - - thisq = mol->getCom()-com; - thisv = mol->getComVel()-comVel; - thisMass = mol->getMass(); - // Compute moment of intertia coefficients. - xx += thisq[0]*thisq[0]*thisMass; - yy += thisq[1]*thisq[1]*thisMass; - zz += thisq[2]*thisq[2]*thisMass; - - // compute products of intertia - xy += thisq[0]*thisq[1]*thisMass; - xz += thisq[0]*thisq[2]*thisMass; - yz += thisq[1]*thisq[2]*thisMass; - - angularMomentum += cross( thisq, thisv ) * thisMass; - - } - - - inertiaTensor(0,0) = yy + zz; - inertiaTensor(0,1) = -xy; - inertiaTensor(0,2) = -xz; - inertiaTensor(1,0) = -xy; - inertiaTensor(1,1) = xx + zz; - inertiaTensor(1,2) = -yz; - inertiaTensor(2,0) = -xz; - inertiaTensor(2,1) = -yz; - inertiaTensor(2,2) = xx + yy; - -#ifdef IS_MPI - 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); -#endif - - return; - } - - //Returns the angular momentum of the system - Vector3d SimInfo::getAngularMomentum(){ - - Vector3d com(0.0); - Vector3d comVel(0.0); - Vector3d angularMomentum(0.0); - - getComAll(com,comVel); - - SimInfo::MoleculeIterator i; - Molecule* mol; - - Vector3d thisr(0.0); - Vector3d thisp(0.0); - - double thisMass; - - for (mol = beginMolecule(i); mol != NULL; mol = nextMolecule(i)) { - thisMass = mol->getMass(); - thisr = mol->getCom()-com; - thisp = (mol->getComVel()-comVel)*thisMass; - - angularMomentum += cross( thisr, thisp ); - - } - -#ifdef IS_MPI - Vector3d tmpAngMom; - MPI_Allreduce(tmpAngMom.getArrayPointer(), angularMomentum.getArrayPointer(),3,MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD); -#endif - - return angularMomentum; - } - - -}//end namespace oopse -