--- branches/development/src/brains/SimInfo.cpp 2010/12/27 18:35:59 1529 +++ branches/development/src/brains/SimInfo.cpp 2011/09/13 22:05:04 1627 @@ -54,23 +54,15 @@ #include "math/Vector3.hpp" #include "primitives/Molecule.hpp" #include "primitives/StuntDouble.hpp" -#include "UseTheForce/fCutoffPolicy.h" -#include "UseTheForce/DarkSide/fSwitchingFunctionType.h" -#include "UseTheForce/doForces_interface.h" -#include "UseTheForce/DarkSide/neighborLists_interface.h" -#include "UseTheForce/DarkSide/switcheroo_interface.h" #include "utils/MemoryUtils.hpp" #include "utils/simError.h" #include "selection/SelectionManager.hpp" #include "io/ForceFieldOptions.hpp" #include "UseTheForce/ForceField.hpp" -#include "nonbonded/InteractionManager.hpp" - - +#include "nonbonded/SwitchingFunction.hpp" #ifdef IS_MPI -#include "UseTheForce/mpiComponentPlan.h" -#include "UseTheForce/DarkSide/simParallel_interface.h" -#endif +#include +#endif using namespace std; namespace OpenMD { @@ -82,7 +74,7 @@ namespace OpenMD { nGlobalIntegrableObjects_(0), nGlobalRigidBodies_(0), nAtoms_(0), nBonds_(0), nBends_(0), nTorsions_(0), nInversions_(0), nRigidBodies_(0), nIntegrableObjects_(0), nCutoffGroups_(0), - nConstraints_(0), sman_(NULL), fortranInitialized_(false), + nConstraints_(0), sman_(NULL), topologyDone_(false), calcBoxDipole_(false), useAtomicVirial_(true) { MoleculeStamp* molStamp; @@ -136,6 +128,7 @@ namespace OpenMD { //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 @@ -279,7 +272,26 @@ namespace OpenMD { #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; @@ -656,27 +668,28 @@ namespace OpenMD { molStampIds_.insert(molStampIds_.end(), nmol, curStampId); } - void SimInfo::update() { - setupSimType(); - setupCutoffRadius(); - setupSwitchingRadius(); - setupCutoffMethod(); - setupSkinThickness(); - setupSwitchingFunction(); - setupAccumulateBoxDipole(); - -#ifdef IS_MPI - setupFortranParallel(); -#endif - setupFortranSim(); - fortranInitialized_ = true; - + /** + * update + * + * Performs the global checks and variable settings after the + * objects have been created. + * + */ + void SimInfo::update() { + setupSimVariables(); calcNdf(); calcNdfRaw(); calcNdfTrans(); } + /** + * 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; @@ -684,115 +697,85 @@ namespace OpenMD { Atom* atom; set atomTypes; - for(mol = beginMolecule(mi); mol != NULL; mol = nextMolecule(mi)) { - for(atom = mol->beginAtom(ai); atom != NULL; atom = mol->nextAtom(ai)) { + for(mol = beginMolecule(mi); mol != NULL; mol = nextMolecule(mi)) { + for(atom = mol->beginAtom(ai); atom != NULL; + atom = mol->nextAtom(ai)) { atomTypes.insert(atom->getAtomType()); } } - return atomTypes; - } - - /** - * setupCutoffRadius - * - * If the cutoffRadius was explicitly set, use that value. - * If the cutoffRadius was not explicitly set: - * Are there electrostatic atoms? Use 12.0 Angstroms. - * No electrostatic atoms? Poll the atom types present in the - * simulation for suggested cutoff values (e.g. 2.5 * sigma). - * Use the maximum suggested value that was found. - */ - void SimInfo::setupCutoffRadius() { - if (simParams_->haveCutoffRadius()) { - cutoffRadius_ = simParams_->getCutoffRadius(); - } else { - if (usesElectrostaticAtoms_) { - sprintf(painCave.errMsg, - "SimInfo Warning: No value was set for the cutoffRadius.\n" - "\tOpenMD will use a default value of 12.0 angstroms" - "\tfor the cutoffRadius.\n"); - painCave.isFatal = 0; - simError(); - cutoffRadius_ = 12.0; - } else { - RealType thisCut; - set::iterator i; - set atomTypes; - atomTypes = getSimulatedAtomTypes(); - for (i = atomTypes.begin(); i != atomTypes.end(); ++i) { - thisCut = InteractionManager::Instance()->getSuggestedCutoffRadius((*i)); - cutoffRadius_ = max(thisCut, cutoffRadius_); - } - sprintf(painCave.errMsg, - "SimInfo Warning: No value was set for the cutoffRadius.\n" - "\tOpenMD will use %lf angstroms.\n", - cutoffRadius_); - painCave.isFatal = 0; - simError(); - } - } +#ifdef IS_MPI - InteractionManager::Instance()->setCutoffRadius(cutoffRadius_); - } - - /** - * setupSwitchingRadius - * - * If the switchingRadius was explicitly set, use that value (but check it) - * If the switchingRadius was not explicitly set: use 0.85 * cutoffRadius_ - */ - void SimInfo::setupSwitchingRadius() { + // loop over the found atom types on this processor, and add their + // numerical idents to a vector: - if (simParams_->haveSwitchingRadius()) { - switchingRadius_ = simParams_->getSwitchingRadius(); - if (switchingRadius_ > cutoffRadius_) { - sprintf(painCave.errMsg, - "SimInfo Error: switchingRadius (%f) is larger than cutoffRadius(%f)\n", - switchingRadius_, cutoffRadius_); - painCave.isFatal = 1; - simError(); + vector foundTypes; + set::iterator i; + for (i = atomTypes.begin(); i != atomTypes.end(); ++i) + foundTypes.push_back( (*i)->getIdent() ); - } - } else { - switchingRadius_ = 0.85 * cutoffRadius_; - sprintf(painCave.errMsg, - "SimInfo Warning: No value was set for the switchingRadius.\n" - "\tOpenMD will use a default value of 85 percent of the cutoffRadius.\n" - "\tswitchingRadius = %f. for this simulation\n", switchingRadius_); - painCave.isFatal = 0; - simError(); - } - InteractionManager::Instance()->setSwitchingRadius(switchingRadius_); - } + // count_local holds the number of found types on this processor + int count_local = foundTypes.size(); - /** - * setupSkinThickness - * - * If the skinThickness was explicitly set, use that value (but check it) - * If the skinThickness was not explicitly set: use 1.0 angstroms - */ - void SimInfo::setupSkinThickness() { - if (simParams_->haveSkinThickness()) { - skinThickness_ = simParams_->getSkinThickness(); - } else { - skinThickness_ = 1.0; - sprintf(painCave.errMsg, - "SimInfo Warning: No value was set for the skinThickness.\n" - "\tOpenMD will use a default value of %f Angstroms\n" - "\tfor this simulation\n", skinThickness_); - painCave.isFatal = 0; - simError(); - } + 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() { + 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(); - - useAtomicVirial_ = simParams_->getUseAtomicVirial(); - + atomTypes = getSimulatedAtomTypes(); int usesElectrostatic = 0; int usesMetallic = 0; int usesDirectional = 0; @@ -802,46 +785,74 @@ namespace OpenMD { usesMetallic |= (*i)->isMetal(); usesDirectional |= (*i)->isDirectional(); } - + #ifdef IS_MPI int temp; temp = usesDirectional; MPI_Allreduce(&temp, &usesDirectionalAtoms_, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD); - + temp = usesMetallic; MPI_Allreduce(&temp, &usesMetallicAtoms_, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD); - + temp = usesElectrostatic; MPI_Allreduce(&temp, &usesElectrostaticAtoms_, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD); +#else + + usesDirectionalAtoms_ = usesDirectional; + usesMetallicAtoms_ = usesMetallic; + usesElectrostaticAtoms_ = usesElectrostatic; + #endif - fInfo_.SIM_uses_PBC = usesPeriodicBoundaries_; - fInfo_.SIM_uses_DirectionalAtoms = usesDirectionalAtoms_; - fInfo_.SIM_uses_MetallicAtoms = usesMetallicAtoms_; - fInfo_.SIM_requires_SkipCorrection = usesElectrostaticAtoms_; - fInfo_.SIM_requires_SelfCorrection = usesElectrostaticAtoms_; - fInfo_.SIM_uses_AtomicVirial = usesAtomicVirial_; + + requiresPrepair_ = usesMetallicAtoms_ ? true : false; + requiresSkipCorrection_ = usesElectrostaticAtoms_ ? true : false; + requiresSelfCorrection_ = usesElectrostaticAtoms_ ? true : false; } - void SimInfo::setupFortranSim() { - int isError; - int nExclude, nOneTwo, nOneThree, nOneFour; - vector fortranGlobalGroupMembership; + + vector SimInfo::getGlobalAtomIndices() { + SimInfo::MoleculeIterator mi; + Molecule* mol; + Molecule::AtomIterator ai; + Atom* atom; + + vector GlobalAtomIndices(getNAtoms(), 0); - notifyFortranSkinThickness(&skinThickness_); + 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; + } - int ljsp = cutoffMethod_ == SHIFTED_POTENTIAL ? 1 : 0; - int ljsf = cutoffMethod_ == SHIFTED_FORCE ? 1 : 0; - notifyFortranCutoffs(&cutoffRadius_, &switchingRadius_, &ljsp, &ljsf); - isError = 0; + vector SimInfo::getGlobalGroupIndices() { + SimInfo::MoleculeIterator mi; + Molecule* mol; + Molecule::CutoffGroupIterator ci; + CutoffGroup* cg; - //globalGroupMembership_ is filled by SimCreator - for (int i = 0; i < nGlobalAtoms_; i++) { - fortranGlobalGroupMembership.push_back(globalGroupMembership_[i] + 1); + vector GlobalGroupIndices; + + for (mol = beginMolecule(mi); mol != NULL; mol = nextMolecule(mi)) { + + //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::prepareTopology() { + int nExclude, nOneTwo, nOneThree, nOneFour; + //calculate mass ratio of cutoff group - vector mfact; SimInfo::MoleculeIterator mi; Molecule* mol; Molecule::CutoffGroupIterator ci; @@ -850,43 +861,43 @@ namespace OpenMD { Atom* atom; 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 !!!) - 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 - vector molMembershipArray(nGlobalAtoms_); - for (int i = 0; i < nGlobalAtoms_; i++) { - molMembershipArray[i] = globalMolMembership_[i] + 1; - } - //setup fortran simulation + //scan topology nExclude = excludedInteractions_.getSize(); nOneTwo = oneTwoInteractions_.getSize(); @@ -898,134 +909,9 @@ namespace OpenMD { 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 = OPENMD_ERROR; - simError(); - } - - - sprintf( checkPointMsg, - "succesfully sent the simulation information to fortran.\n"); - - errorCheckPoint(); - - // Setup number of neighbors in neighbor list if present - if (simParams_->haveNeighborListNeighbors()) { - int nlistNeighbors = simParams_->getNeighborListNeighbors(); - setNeighbors(&nlistNeighbors); - } - - - } - - - void SimInfo::setupFortranParallel() { -#ifdef IS_MPI - //SimInfo is responsible for creating localToGlobalAtomIndex and localToGlobalGroupIndex - vector localToGlobalAtomIndex(getNAtoms(), 0); - 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"); - errorCheckPoint(); - -#endif - } - - - void SimInfo::setupSwitchingFunction() { - int ft = CUBIC; - - if (simParams_->haveSwitchingFunctionType()) { - string funcType = simParams_->getSwitchingFunctionType(); - toUpper(funcType); - if (funcType == "CUBIC") { - ft = CUBIC; - } else { - if (funcType == "FIFTH_ORDER_POLYNOMIAL") { - ft = FIFTH_ORDER_POLY; - } else { - // throw error - sprintf( painCave.errMsg, - "SimInfo 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); - + topologyDone_ = true; } - void SimInfo::setupAccumulateBoxDipole() { - - // we only call setAccumulateBoxDipole if the accumulateBoxDipole parameter is true - if ( simParams_->haveAccumulateBoxDipole() ) - if ( simParams_->getAccumulateBoxDipole() ) { - calcBoxDipole_ = true; - } - - } - void SimInfo::addProperty(GenericData* genData) { properties_.addProperty(genData); } @@ -1060,9 +946,11 @@ namespace OpenMD { 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)) { @@ -1072,6 +960,10 @@ namespace OpenMD { 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_); } }