--- trunk/src/brains/SimInfo.cpp 2007/02/26 04:45:42 1121 +++ branches/development/src/brains/SimInfo.cpp 2010/10/02 19:54:41 1503 @@ -6,19 +6,10 @@ * redistribute this software in source and binary code form, provided * that the following conditions are met: * - * 1. Acknowledgement of the program authors must be made in any - * publication of scientific results based in part on use of the - * program. An acceptable form of acknowledgement is citation of - * the article in which the program was described (Matthew - * A. Meineke, Charles F. Vardeman II, Teng Lin, Christopher - * J. Fennell and J. Daniel Gezelter, "OOPSE: An Object-Oriented - * Parallel Simulation Engine for Molecular Dynamics," - * J. Comput. Chem. 26, pp. 252-271 (2005)) - * - * 2. Redistributions of source code must retain the above copyright + * 1. Redistributions of source code must retain the above copyright * notice, this list of conditions and the following disclaimer. * - * 3. Redistributions in binary form must reproduce the above copyright + * 2. Redistributions in binary form must reproduce the above copyright * notice, this list of conditions and the following disclaimer in the * documentation and/or other materials provided with the * distribution. @@ -37,6 +28,15 @@ * arising out of the use of or inability to use software, even if the * University of Notre Dame has been advised of the possibility of * such damages. + * + * SUPPORT OPEN SCIENCE! If you use OpenMD or its source code in your + * research, please cite the appropriate papers when you publish your + * work. Good starting points are: + * + * [1] Meineke, et al., J. Comp. Chem. 26, 252-271 (2005). + * [2] Fennell & Gezelter, J. Chem. Phys. 124, 234104 (2006). + * [3] Sun, Lin & Gezelter, J. Chem. Phys. 128, 24107 (2008). + * [4] Vardeman & Gezelter, in progress (2009). */ /** @@ -55,12 +55,9 @@ #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" @@ -74,7 +71,7 @@ #include "UseTheForce/DarkSide/simParallel_interface.h" #endif -namespace oopse { +namespace OpenMD { std::set getRigidSet(int index, std::map >& container) { std::map >::iterator i = container.find(index); std::set result; @@ -90,10 +87,12 @@ namespace oopse { 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) { + 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 @@ -101,6 +100,7 @@ namespace oopse { CutoffGroupStamp* cgStamp; RigidBodyStamp* rbStamp; int nRigidAtoms = 0; + std::vector components = simParams->getComponents(); for (std::vector::iterator i = components.begin(); i !=components.end(); ++i) { @@ -155,11 +155,7 @@ namespace oopse { + nGlobalRigidBodies_; nGlobalMols_ = molStampIds_.size(); - -#ifdef IS_MPI molToProcMap_.resize(nGlobalMols_); -#endif - } SimInfo::~SimInfo() { @@ -197,13 +193,14 @@ namespace oopse { nBonds_ += mol->getNBonds(); nBends_ += mol->getNBends(); nTorsions_ += mol->getNTorsions(); + nInversions_ += mol->getNInversions(); nRigidBodies_ += mol->getNRigidBodies(); nIntegrableObjects_ += mol->getNIntegrableObjects(); nCutoffGroups_ += mol->getNCutoffGroups(); nConstraints_ += mol->getNConstraintPairs(); - addExcludePairs(mol); - + addInteractionPairs(mol); + return true; } else { return false; @@ -222,12 +219,13 @@ namespace oopse { nBonds_ -= mol->getNBonds(); nBends_ -= mol->getNBends(); nTorsions_ -= mol->getNTorsions(); + nInversions_ -= mol->getNInversions(); nRigidBodies_ -= mol->getNRigidBodies(); nIntegrableObjects_ -= mol->getNIntegrableObjects(); nCutoffGroups_ -= mol->getNCutoffGroups(); nConstraints_ -= mol->getNConstraintPairs(); - removeExcludePairs(mol); + removeInteractionPairs(mol); molecules_.erase(mol->getGlobalIndex()); delete mol; @@ -353,150 +351,202 @@ namespace oopse { } - void SimInfo::addExcludePairs(Molecule* mol) { + void SimInfo::addInteractionPairs(Molecule* mol) { + ForceFieldOptions& options_ = forceField_->getForceFieldOptions(); std::vector::iterator bondIter; std::vector::iterator bendIter; std::vector::iterator torsionIter; + std::vector::iterator inversionIter; Bond* bond; Bend* bend; Torsion* torsion; + Inversion* inversion; int a; int b; int c; int d; - 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. + std::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); + std::vector atoms = rb->getAtoms(); + std::set rigidAtoms; + for (int i = 0; i < static_cast(atoms.size()); ++i) { + rigidAtoms.insert(atoms[i]->getGlobalIndex()); + } + for (int i = 0; i < static_cast(atoms.size()); ++i) { + atomGroups.insert(std::map >::value_type(atoms[i]->getGlobalIndex(), rigidAtoms)); + } } else { std::set oneAtomSet; oneAtomSet.insert(integrableObject->getGlobalIndex()); atomGroups.insert(std::map >::value_type(integrableObject->getGlobalIndex(), oneAtomSet)); } } + + for (bond= mol->beginBond(bondIter); bond != NULL; + bond = mol->nextBond(bondIter)) { - - - 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)) { + for (inversion= mol->beginInversion(inversionIter); inversion != NULL; + inversion = mol->nextInversion(inversionIter)) { + + a = inversion->getAtomA()->getGlobalIndex(); + b = inversion->getAtomB()->getGlobalIndex(); + c = inversion->getAtomC()->getGlobalIndex(); + d = inversion->getAtomD()->getGlobalIndex(); + + if (options_.havevdw12scale() || options_.haveelectrostatic12scale()) { + oneTwoInteractions_.addPair(a, b); + oneTwoInteractions_.addPair(a, c); + oneTwoInteractions_.addPair(a, d); + } else { + excludedInteractions_.addPair(a, b); + excludedInteractions_.addPair(a, c); + excludedInteractions_.addPair(a, d); + } + + if (options_.havevdw13scale() || options_.haveelectrostatic13scale()) { + oneThreeInteractions_.addPair(b, c); + oneThreeInteractions_.addPair(b, d); + oneThreeInteractions_.addPair(c, d); + } else { + excludedInteractions_.addPair(b, c); + excludedInteractions_.addPair(b, d); + excludedInteractions_.addPair(c, d); + } + } + + for (rb = mol->beginRigidBody(rbIter); rb != NULL; + rb = mol->nextRigidBody(rbIter)) { std::vector atoms = rb->getAtoms(); - for (int i = 0; i < atoms.size() -1 ; ++i) { - for (int j = i + 1; j < atoms.size(); ++j) { + for (int i = 0; i < static_cast(atoms.size()) -1 ; ++i) { + for (int j = i + 1; j < static_cast(atoms.size()); ++j) { a = atoms[i]->getGlobalIndex(); b = atoms[j]->getGlobalIndex(); - exclude_.addPair(a, b); + excludedInteractions_.addPair(a, b); } } } } - void SimInfo::removeExcludePairs(Molecule* mol) { + void SimInfo::removeInteractionPairs(Molecule* mol) { + ForceFieldOptions& options_ = forceField_->getForceFieldOptions(); std::vector::iterator bondIter; std::vector::iterator bendIter; std::vector::iterator torsionIter; + std::vector::iterator inversionIter; Bond* bond; Bend* bend; Torsion* torsion; + Inversion* inversion; int a; int b; int c; int d; std::map > atomGroups; - Molecule::RigidBodyIterator rbIter; RigidBody* rb; Molecule::IntegrableObjectIterator ii; StuntDouble* integrableObject; - for (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); + std::vector atoms = rb->getAtoms(); + std::set rigidAtoms; + for (int i = 0; i < static_cast(atoms.size()); ++i) { + rigidAtoms.insert(atoms[i]->getGlobalIndex()); + } + for (int i = 0; i < static_cast(atoms.size()); ++i) { + atomGroups.insert(std::map >::value_type(atoms[i]->getGlobalIndex(), rigidAtoms)); + } } else { std::set oneAtomSet; oneAtomSet.insert(integrableObject->getGlobalIndex()); @@ -504,84 +554,121 @@ namespace oopse { } } - - 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)) { + for (rb = mol->beginRigidBody(rbIter); rb != NULL; + rb = mol->nextRigidBody(rbIter)) { std::vector atoms = rb->getAtoms(); - for (int i = 0; i < atoms.size() -1 ; ++i) { - for (int j = i + 1; j < atoms.size(); ++j) { + for (int i = 0; i < static_cast(atoms.size()) -1 ; ++i) { + for (int j = i + 1; j < static_cast(atoms.size()); ++j) { a = atoms[i]->getGlobalIndex(); b = atoms[j]->getGlobalIndex(); - exclude_.removePair(a, b); + excludedInteractions_.removePair(a, b); } } } - + } - - + + void SimInfo::addMoleculeStamp(MoleculeStamp* molStamp, int nmol) { int curStampId; - + //index from 0 curStampId = moleculeStamps_.size(); @@ -666,12 +753,14 @@ namespace oopse { int useSF; int useSP; int useBoxDipole; + std::string myMethod; // set the useRF logical useRF = 0; useSF = 0; useSP = 0; + useBoxDipole = 0; if (simParams_->haveElectrostaticSummationMethod()) { @@ -690,6 +779,8 @@ namespace oopse { if (simParams_->getAccumulateBoxDipole()) useBoxDipole = 1; + useAtomicVirial_ = simParams_->getUseAtomicVirial(); + //loop over all of the atom types for (i = atomTypes.begin(); i != atomTypes.end(); ++i) { useLennardJones |= (*i)->isLennardJones(); @@ -767,6 +858,9 @@ namespace oopse { temp = useBoxDipole; MPI_Allreduce(&temp, &useBoxDipole, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD); + temp = useAtomicVirial_; + MPI_Allreduce(&temp, &useAtomicVirial_, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD); + #endif fInfo_.SIM_uses_PBC = usePBC; @@ -786,14 +880,14 @@ namespace oopse { fInfo_.SIM_uses_SF = useSF; fInfo_.SIM_uses_SP = useSP; fInfo_.SIM_uses_BoxDipole = useBoxDipole; + fInfo_.SIM_uses_AtomicVirial = useAtomicVirial_; } void SimInfo::setupFortranSim() { int isError; - int nExclude; + int nExclude, nOneTwo, nOneThree, nOneFour; std::vector fortranGlobalGroupMembership; - nExclude = exclude_.getSize(); isError = 0; //globalGroupMembership_ is filled by SimCreator @@ -825,7 +919,6 @@ namespace oopse { else mfact.push_back( 1.0 ); } - } } @@ -849,28 +942,40 @@ namespace oopse { } //setup fortran simulation - int nGlobalExcludes = 0; - int* globalExcludes = NULL; - int* excludeList = exclude_.getExcludeList(); - setFortranSim( &fInfo_, &nGlobalAtoms_, &nAtoms_, &identArray[0], &nExclude, excludeList , - &nGlobalExcludes, globalExcludes, &molMembershipArray[0], - &mfact[0], &nCutoffGroups_, &fortranGlobalGroupMembership[0], &isError); - if( isError ){ + nExclude = excludedInteractions_.getSize(); + nOneTwo = oneTwoInteractions_.getSize(); + nOneThree = oneThreeInteractions_.getSize(); + nOneFour = oneFourInteractions_.getSize(); + int* excludeList = excludedInteractions_.getPairList(); + int* oneTwoList = oneTwoInteractions_.getPairList(); + int* oneThreeList = oneThreeInteractions_.getPairList(); + int* oneFourList = oneFourInteractions_.getPairList(); + + setFortranSim( &fInfo_, &nGlobalAtoms_, &nAtoms_, &identArray[0], + &nExclude, excludeList, + &nOneTwo, oneTwoList, + &nOneThree, oneThreeList, + &nOneFour, oneFourList, + &molMembershipArray[0], &mfact[0], &nCutoffGroups_, + &fortranGlobalGroupMembership[0], &isError); + + if( isError ){ + sprintf( painCave.errMsg, "There was an error setting the simulation information in fortran.\n" ); painCave.isFatal = 1; - painCave.severity = OOPSE_ERROR; + painCave.severity = OPENMD_ERROR; simError(); } - -#ifdef IS_MPI + + sprintf( checkPointMsg, "succesfully sent the simulation information to fortran.\n"); - MPIcheckPoint(); -#endif // is_mpi - + + errorCheckPoint(); + // Setup number of neighbors in neighbor list if present if (simParams_->haveNeighborListNeighbors()) { int nlistNeighbors = simParams_->getNeighborListNeighbors(); @@ -881,9 +986,8 @@ 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; @@ -933,12 +1037,10 @@ namespace oopse { } sprintf(checkPointMsg, " mpiRefresh successful.\n"); - MPIcheckPoint(); + errorCheckPoint(); - - } - #endif + } void SimInfo::setupCutoff() { @@ -947,6 +1049,10 @@ namespace oopse { // Check the cutoff policy int cp = TRADITIONAL_CUTOFF_POLICY; // Set to traditional by default + // Set LJ shifting bools to false + ljsp_ = 0; + ljsf_ = 0; + std::string myPolicy; if (forceFieldOptions_.haveCutoffPolicy()){ myPolicy = forceFieldOptions_.getCutoffPolicy(); @@ -996,7 +1102,7 @@ namespace oopse { 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" + "\tOpenMD will use a default value of 85 percent of the cutoffRadius.\n" "\tswitchingRadius = %f. for this simulation\n", rsw_); painCave.isFatal = 0; simError(); @@ -1004,22 +1110,33 @@ namespace oopse { 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" + "\tOpenMD will use the same value as the cutoffRadius.\n" "\tswitchingRadius = %f. for this simulation\n", rsw_); painCave.isFatal = 0; simError(); } } + + if (simParams_->haveElectrostaticSummationMethod()) { + std::string myMethod = simParams_->getElectrostaticSummationMethod(); + toUpper(myMethod); + + if (myMethod == "SHIFTED_POTENTIAL") { + ljsp_ = 1; + } else if (myMethod == "SHIFTED_FORCE") { + ljsf_ = 1; + } + } + + notifyFortranCutoffs(&rcut_, &rsw_, &ljsp_, &ljsf_); - 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" + "\tOpenMD will use a default value of 15.0 angstroms" "\tfor the cutoffRadius.\n"); painCave.isFatal = 0; simError(); @@ -1028,7 +1145,15 @@ namespace oopse { if (simParams_->haveElectrostaticSummationMethod()) { std::string myMethod = simParams_->getElectrostaticSummationMethod(); toUpper(myMethod); - if (myMethod == "SHIFTED_POTENTIAL" || myMethod == "SHIFTED_FORCE") { + + // For the time being, we're tethering the LJ shifted behavior to the + // electrostaticSummationMethod keyword options + if (myMethod == "SHIFTED_POTENTIAL") { + ljsp_ = 1; + } else if (myMethod == "SHIFTED_FORCE") { + ljsf_ = 1; + } + if (myMethod == "SHIFTED_POTENTIAL" || myMethod == "SHIFTED_FORCE") { if (simParams_->haveSwitchingRadius()){ sprintf(painCave.errMsg, "SimInfo Warning: A value was set for the switchingRadius\n" @@ -1045,13 +1170,16 @@ namespace oopse { } else { sprintf(painCave.errMsg, "SimCreator Warning: No value was set for switchingRadius.\n" - "\tOOPSE will use a default value of\n" + "\tOpenMD 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_); + + Electrostatic::setElectrostaticCutoffRadius(rcut_, rsw_); + notifyFortranCutoffs(&rcut_, &rsw_, &ljsp_, &ljsf_); + } else { // We didn't set rcut explicitly, and we don't have electrostatic atoms, so // We'll punt and let fortran figure out the cutoffs later. @@ -1065,8 +1193,8 @@ namespace oopse { void SimInfo::setupElectrostaticSummationMethod( int isError ) { int errorOut; - int esm = NONE; - int sm = UNDAMPED; + ElectrostaticSummationMethod esm = NONE; + ElectrostaticScreeningMethod sm = UNDAMPED; RealType alphaVal; RealType dielectric; @@ -1153,12 +1281,11 @@ namespace oopse { } } - // let's pass some summation method variables to fortran - setElectrostaticSummationMethod( &esm ); - setFortranElectrostaticMethod( &esm ); - setScreeningMethod( &sm ); - setDampingAlpha( &alphaVal ); - setReactionFieldDielectric( &dielectric ); + + Electrostatic::setElectrostaticSummationMethod( esm ); + Electrostatic::setElectrostaticScreeningMethod( sm ); + Electrostatic::setDampingAlpha( alphaVal ); + Electrostatic::setReactionFieldDielectric( dielectric ); initFortranFF( &errorOut ); } @@ -1504,5 +1631,5 @@ namespace oopse { return sdByGlobalIndex_.at(index); } */ -}//end namespace oopse +}//end namespace OpenMD