--- trunk/src/brains/SimInfo.cpp 2005/11/03 23:12:27 716 +++ trunk/src/brains/SimInfo.cpp 2005/11/16 23:10:02 749 @@ -48,6 +48,7 @@ #include #include +#include #include "brains/SimInfo.hpp" #include "math/Vector3.hpp" @@ -55,9 +56,11 @@ #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 "utils/MemoryUtils.hpp" #include "utils/simError.h" #include "selection/SelectionManager.hpp" @@ -68,7 +71,16 @@ namespace oopse { #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; + } + SimInfo::SimInfo(MakeStamps* stamps, std::vector >& molStampPairs, ForceField* ff, Globals* simParams) : stamps_(stamps), forceField_(ff), simParams_(simParams), @@ -343,7 +355,36 @@ namespace oopse { 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)) { + + 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)); + } + } else { + std::set oneAtomSet; + oneAtomSet.insert(integrableObject->getGlobalIndex()); + atomGroups.insert(std::map >::value_type(integrableObject->getGlobalIndex(), oneAtomSet)); + } + } + + + for (bond= mol->beginBond(bondIter); bond != NULL; bond = mol->nextBond(bondIter)) { a = bond->getAtomA()->getGlobalIndex(); b = bond->getAtomB()->getGlobalIndex(); @@ -354,10 +395,17 @@ namespace oopse { 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_.addPair(a, b); - exclude_.addPair(a, c); - exclude_.addPair(b, c); + exclude_.addPairs(rigidSetA, rigidSetB); + exclude_.addPairs(rigidSetA, rigidSetC); + exclude_.addPairs(rigidSetB, rigidSetC); + + //exclude_.addPair(a, b); + //exclude_.addPair(a, c); + //exclude_.addPair(b, c); } for (torsion= mol->beginTorsion(torsionIter); torsion != NULL; torsion = mol->nextTorsion(torsionIter)) { @@ -365,17 +413,36 @@ namespace oopse { 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); + exclude_.addPairs(rigidSetA, rigidSetB); + exclude_.addPairs(rigidSetA, rigidSetC); + exclude_.addPairs(rigidSetA, rigidSetD); + exclude_.addPairs(rigidSetB, rigidSetC); + exclude_.addPairs(rigidSetB, rigidSetD); + exclude_.addPairs(rigidSetC, rigidSetD); + + /* + 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); + */ } - 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) { @@ -400,6 +467,34 @@ namespace oopse { 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)) { + + 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)); + } + } else { + std::set oneAtomSet; + oneAtomSet.insert(integrableObject->getGlobalIndex()); + atomGroups.insert(std::map >::value_type(integrableObject->getGlobalIndex(), oneAtomSet)); + } + } + for (bond= mol->beginBond(bondIter); bond != NULL; bond = mol->nextBond(bondIter)) { a = bond->getAtomA()->getGlobalIndex(); @@ -412,9 +507,17 @@ namespace oopse { b = bend->getAtomB()->getGlobalIndex(); c = bend->getAtomC()->getGlobalIndex(); - exclude_.removePair(a, b); - exclude_.removePair(a, c); - exclude_.removePair(b, c); + 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); } for (torsion= mol->beginTorsion(torsionIter); torsion != NULL; torsion = mol->nextTorsion(torsionIter)) { @@ -423,16 +526,36 @@ namespace oopse { 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); + + exclude_.removePairs(rigidSetA, rigidSetB); + exclude_.removePairs(rigidSetA, rigidSetC); + exclude_.removePairs(rigidSetA, rigidSetD); + exclude_.removePairs(rigidSetB, rigidSetC); + exclude_.removePairs(rigidSetB, rigidSetD); + exclude_.removePairs(rigidSetC, rigidSetD); + + /* + 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()); + + exclude_.removePair(a, b); exclude_.removePair(a, c); exclude_.removePair(a, d); exclude_.removePair(b, c); exclude_.removePair(b, d); exclude_.removePair(c, 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) { @@ -472,6 +595,7 @@ namespace oopse { int isError = 0; setupElectrostaticSummationMethod( isError ); + setupSwitchingFunction(); if(isError){ sprintf( painCave.errMsg, @@ -516,6 +640,7 @@ namespace oopse { int useLennardJones = 0; int useElectrostatic = 0; int useEAM = 0; + int useSC = 0; int useCharge = 0; int useDirectional = 0; int useDipole = 0; @@ -529,12 +654,12 @@ namespace oopse { //usePBC and useRF are from simParams int usePBC = simParams_->getUsePeriodicBoundaryConditions(); int useRF; - int useDW; + int useSF; std::string myMethod; // set the useRF logical useRF = 0; - useDW = 0; + useSF = 0; if (simParams_->haveElectrostaticSummationMethod()) { @@ -543,8 +668,8 @@ namespace oopse { if (myMethod == "REACTION_FIELD") { useRF=1; } else { - if (myMethod == "SHIFTED_POTENTIAL") { - useDW = 1; + if (myMethod == "SHIFTED_FORCE") { + useSF = 1; } } } @@ -554,6 +679,7 @@ namespace oopse { useLennardJones |= (*i)->isLennardJones(); useElectrostatic |= (*i)->isElectrostatic(); useEAM |= (*i)->isEAM(); + useSC |= (*i)->isSC(); useCharge |= (*i)->isCharge(); useDirectional |= (*i)->isDirectional(); useDipole |= (*i)->isDipole(); @@ -604,6 +730,9 @@ namespace oopse { temp = useEAM; MPI_Allreduce(&temp, &useEAM, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD); + temp = useSC; + MPI_Allreduce(&temp, &useSC, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD); + temp = useShape; MPI_Allreduce(&temp, &useShape, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD); @@ -613,8 +742,8 @@ namespace oopse { temp = useRF; MPI_Allreduce(&temp, &useRF, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD); - temp = useDW; - MPI_Allreduce(&temp, &useDW, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD); + temp = useSF; + MPI_Allreduce(&temp, &useSF, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD); #endif @@ -628,10 +757,11 @@ namespace oopse { fInfo_.SIM_uses_StickyPower = useStickyPower; fInfo_.SIM_uses_GayBerne = useGayBerne; fInfo_.SIM_uses_EAM = useEAM; + fInfo_.SIM_uses_SC = useSC; fInfo_.SIM_uses_Shapes = useShape; fInfo_.SIM_uses_FLARB = useFLARB; fInfo_.SIM_uses_RF = useRF; - fInfo_.SIM_uses_DampedWolf = useDW; + fInfo_.SIM_uses_SF = useSF; if( myMethod == "REACTION_FIELD") { @@ -978,6 +1108,32 @@ namespace oopse { 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); }