--- trunk/src/brains/SimInfo.cpp 2007/02/26 04:45:42 1121 +++ trunk/src/brains/SimInfo.cpp 2008/07/14 12:35:58 1277 @@ -90,10 +90,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 +103,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 +158,7 @@ namespace oopse { + nGlobalRigidBodies_; nGlobalMols_ = molStampIds_.size(); - -#ifdef IS_MPI molToProcMap_.resize(nGlobalMols_); -#endif - } SimInfo::~SimInfo() { @@ -197,6 +196,7 @@ namespace oopse { nBonds_ += mol->getNBonds(); nBends_ += mol->getNBends(); nTorsions_ += mol->getNTorsions(); + nInversions_ += mol->getNInversions(); nRigidBodies_ += mol->getNRigidBodies(); nIntegrableObjects_ += mol->getNIntegrableObjects(); nCutoffGroups_ += mol->getNCutoffGroups(); @@ -222,6 +222,7 @@ namespace oopse { nBonds_ -= mol->getNBonds(); nBends_ -= mol->getNBends(); nTorsions_ -= mol->getNTorsions(); + nInversions_ -= mol->getNInversions(); nRigidBodies_ -= mol->getNRigidBodies(); nIntegrableObjects_ -= mol->getNIntegrableObjects(); nCutoffGroups_ -= mol->getNCutoffGroups(); @@ -357,9 +358,11 @@ namespace oopse { 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; @@ -452,6 +455,42 @@ namespace oopse { */ } + 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(); + 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); + */ + } + 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) { @@ -469,9 +508,11 @@ namespace oopse { 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; @@ -534,7 +575,43 @@ 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_.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); + */ + } + 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(); + std::set rigidSetA = getRigidSet(a, atomGroups); std::set rigidSetB = getRigidSet(b, atomGroups); std::set rigidSetC = getRigidSet(c, atomGroups); @@ -666,6 +743,7 @@ namespace oopse { int useSF; int useSP; int useBoxDipole; + std::string myMethod; // set the useRF logical @@ -689,6 +767,8 @@ namespace oopse { if (simParams_->haveAccumulateBoxDipole()) if (simParams_->getAccumulateBoxDipole()) useBoxDipole = 1; + + useAtomicVirial_ = simParams_->getUseAtomicVirial(); //loop over all of the atom types for (i = atomTypes.begin(); i != atomTypes.end(); ++i) { @@ -767,6 +847,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,6 +869,7 @@ 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() { @@ -852,25 +936,26 @@ namespace oopse { 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); - + setFortranSim( &fInfo_, &nGlobalAtoms_, &nAtoms_, &identArray[0], + &nExclude, excludeList , &nGlobalExcludes, globalExcludes, + &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; 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 +966,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,13 +1017,11 @@ namespace oopse { } sprintf(checkPointMsg, " mpiRefresh successful.\n"); - MPIcheckPoint(); + errorCheckPoint(); - +#endif } -#endif - void SimInfo::setupCutoff() { ForceFieldOptions& forceFieldOptions_ = forceField_->getForceFieldOptions(); @@ -947,6 +1029,10 @@ namespace oopse { // Check the cutoff policy int cp = TRADITIONAL_CUTOFF_POLICY; // Set to traditional by default + // Set LJ shifting bools to false + ljsp_ = false; + ljsf_ = false; + std::string myPolicy; if (forceFieldOptions_.haveCutoffPolicy()){ myPolicy = forceFieldOptions_.getCutoffPolicy(); @@ -1010,8 +1096,18 @@ namespace oopse { simError(); } } - - notifyFortranCutoffs(&rcut_, &rsw_); + + if (simParams_->haveElectrostaticSummationMethod()) { + std::string myMethod = simParams_->getElectrostaticSummationMethod(); + toUpper(myMethod); + + if (myMethod == "SHIFTED_POTENTIAL") { + ljsp_ = true; + } else if (myMethod == "SHIFTED_FORCE") { + ljsf_ = true; + } + } + notifyFortranCutoffs(&rcut_, &rsw_, &ljsp_, &ljsf_); } else { @@ -1028,7 +1124,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_ = true; + } else if (myMethod == "SHIFTED_FORCE") { + ljsf_ = true; + } + if (myMethod == "SHIFTED_POTENTIAL" || myMethod == "SHIFTED_FORCE") { if (simParams_->haveSwitchingRadius()){ sprintf(painCave.errMsg, "SimInfo Warning: A value was set for the switchingRadius\n" @@ -1051,7 +1155,9 @@ namespace oopse { simError(); rsw_ = 0.85 * rcut_; } - notifyFortranCutoffs(&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.