--- branches/new_design/OOPSE-2.0/src/brains/SimInfo.cpp 2004/11/12 06:19:04 1733 +++ branches/new_design/OOPSE-2.0/src/brains/SimInfo.cpp 2005/01/04 22:18:36 1901 @@ -31,58 +31,90 @@ */ #include +#include #include "brains/SimInfo.hpp" +#include "math/Vector3.hpp" +#include "primitives/Molecule.hpp" +#include "UseTheForce/doForces_interface.h" +#include "UseTheForce/notifyCutoffs_interface.h" #include "utils/MemoryUtils.hpp" +#include "utils/simError.h" +#ifdef IS_MPI +#include "UseTheForce/mpiComponentPlan.h" +#include "UseTheForce/DarkSide/simParallel_interface.h" +#endif + namespace oopse { -SimInfo::SimInfo(const std::vector >& molStampPairs, - ForceField* ff, Globals* globals) : - forceField_(ff), globals_(globals), nAtoms_(0), nBonds_(0), - nBends_(0), nTorsions_(0), nRigidBodies_(0), nIntegrableObjects_(0), - nCutoffGroups_(0), nConstraints_(0), nZConstraint_(0), sman_(NULL) { +SimInfo::SimInfo(std::vector >& molStampPairs, + ForceField* ff, Globals* simParams) : + forceField_(ff), simParams_(simParams), + ndf_(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; - int nCutoffAtoms; // number of atoms belong to cutoff groups - int ngroups; //total cutoff groups defined in meta-data file MoleculeStamp* molStamp; int nMolWithSameStamp; - CutoffGroupStamp* cgStamp; - int nAtomsIngroups; - int nCutoffGroupsInStamp; - - nGlobalAtoms_ = 0; - ngroups = 0; + 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; 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(); - nAtomsIngroups = 0; - nCutoffGroupsInStamp = molStamp->getNCutoffGroups(); - for (int j=0; j < nCutoffGroupsInStamp; j++) { cgStamp = molStamp->getCutoffGroup(j); - nAtomsIngroups += cgStamp->getNMembers(); + nAtomsInGroups += cgStamp->getNMembers(); } - ngroups += *nMolWithSameStamp; - nCutoffAtoms += nAtomsIngroups * nMolWithSameStamp; + nGroups += nCutoffGroupsInStamp * nMolWithSameStamp; + nCutoffAtoms += nAtomsInGroups * nMolWithSameStamp; + + //calculate atoms in rigid bodies + int nAtomsInRigidBodies = 0; + int nRigidBodiesInStamp = molStamp->getNCutoffGroups(); + + for (int j=0; j < nRigidBodiesInStamp; j++) { + rbStamp = molStamp->getRigidBody(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; + nGlobalCutoffGroups_ = nGlobalAtoms_ - nCutoffAtoms + nGroups; - //initialize globalGroupMembership_, every element of this array will be 0 - globalGroupMembership_.insert(globalGroupMembership_.end(), nGlobalAtoms_, 0); + //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(); @@ -98,7 +130,7 @@ SimInfo::~SimInfo() { MemoryUtils::deleteVectorOfPointer(moleculeStamps_); delete sman_; - delete globals_; + delete simParams_; delete forceField_; } @@ -108,9 +140,9 @@ bool SimInfo::addMolecule(Molecule* mol) { MoleculeIterator i; i = molecules_.find(mol->getGlobalIndex()); - if (i != molecules_.end() ) { + if (i == molecules_.end() ) { - molecules_.insert(make_pair(mol->getGlobalIndex(), mol)); + molecules_.insert(std::make_pair(mol->getGlobalIndex(), mol)); nAtoms_ += mol->getNAtoms(); nBonds_ += mol->getNBonds(); @@ -119,8 +151,10 @@ bool SimInfo::addMolecule(Molecule* mol) { nRigidBodies_ += mol->getNRigidBodies(); nIntegrableObjects_ += mol->getNIntegrableObjects(); nCutoffGroups_ += mol->getNCutoffGroups(); - nConstraints_ += mol->getNConstraints(); + nConstraints_ += mol->getNConstraintPairs(); + addExcludePairs(mol); + return true; } else { return false; @@ -142,8 +176,9 @@ bool SimInfo::removeMolecule(Molecule* mol) { nRigidBodies_ -= mol->getNRigidBodies(); nIntegrableObjects_ -= mol->getNIntegrableObjects(); nCutoffGroups_ -= mol->getNCutoffGroups(); - nConstraints_ -= mol->getNConstraints(); + nConstraints_ -= mol->getNConstraintPairs(); + removeExcludePairs(mol); molecules_.erase(mol->getGlobalIndex()); delete mol; @@ -159,12 +194,12 @@ Molecule* SimInfo::beginMolecule(MoleculeIterator& i) Molecule* SimInfo::beginMolecule(MoleculeIterator& i) { i = molecules_.begin(); - return i == molecules_.end() ? NULL : *i; + return i == molecules_.end() ? NULL : i->second; } Molecule* SimInfo::nextMolecule(MoleculeIterator& i) { ++i; - return i == molecules_.end() ? NULL : *i; + return i == molecules_.end() ? NULL : i->second; } @@ -205,7 +240,7 @@ void SimInfo::calcNdf() { // nZconstraints_ is global, as are the 3 COM translations for the // entire system: - ndf_ = ndf_ - 3 - nZconstraints_; + ndf_ = ndf_ - 3 - nZconstraint_; } @@ -256,7 +291,7 @@ void SimInfo::calcNdfTrans() { ndfTrans_ = ndfTrans_local; #endif - ndfTrans_ = ndfTrans_ - 3 - nZconstraints_; + ndfTrans_ = ndfTrans_ - 3 - nZconstraint_; } @@ -288,7 +323,7 @@ void SimInfo::addExcludePairs(Molecule* mol) { exclude_.addPair(b, c); } - for (torsion= mol->beginTorsion(torsionIter); torsion != NULL; torsion = mol->nextBond(torsionIter)) { + for (torsion= mol->beginTorsion(torsionIter); torsion != NULL; torsion = mol->nextTorsion(torsionIter)) { a = torsion->getAtomA()->getGlobalIndex(); b = torsion->getAtomB()->getGlobalIndex(); c = torsion->getAtomC()->getGlobalIndex(); @@ -333,7 +368,7 @@ void SimInfo::removeExcludePairs(Molecule* mol) { exclude_.removePair(b, c); } - for (torsion= mol->beginTorsion(torsionIter); torsion != NULL; torsion = mol->nextBond(torsionIter)) { + for (torsion= mol->beginTorsion(torsionIter); torsion != NULL; torsion = mol->nextTorsion(torsionIter)) { a = torsion->getAtomA()->getGlobalIndex(); b = torsion->getAtomB()->getGlobalIndex(); c = torsion->getAtomC()->getGlobalIndex(); @@ -354,29 +389,47 @@ void SimInfo::addMoleculeStamp(MoleculeStamp* molStamp int curStampId; //index from 0 - curStampId = molStampIds_.size(); + curStampId = moleculeStamps_.size(); moleculeStamps_.push_back(molStamp); - molStampIds_.insert(molStampIds_.end(), nmol, curStampId) + molStampIds_.insert(molStampIds_.end(), nmol, curStampId); } void SimInfo::update() { + setupSimType(); + #ifdef IS_MPI setupFortranParallel(); #endif setupFortranSim(); + + //setup fortran force field + /** @deprecate */ + int isError = 0; + initFortranFF( &fInfo_.SIM_uses_RF , &isError ); + if(isError){ + sprintf( painCave.errMsg, + "ForceField error: There was an error initializing the forceField in fortran.\n" ); + painCave.isFatal = 1; + simError(); + } + + + setupCutoff(); calcNdf(); calcNdfRaw(); calcNdfTrans(); + + fortranInitialized_ = true; } std::set SimInfo::getUniqueAtomTypes() { - typename SimInfo::MoleculeIterator mi; + SimInfo::MoleculeIterator mi; Molecule* mol; - typename Molecule::AtomIterator ai; + Molecule::AtomIterator ai; Atom* atom; std::set atomTypes; @@ -408,21 +461,21 @@ void SimInfo::setupSimType() { int useFLARB = 0; //it is not in AtomType yet int useDirectionalAtom = 0; int useElectrostatics = 0; - //usePBC and useRF are from globals - bool usePBC = globals_->getPBC(); - bool useRF = globals_->getUseRF(); + //usePBC and useRF are from simParams + int usePBC = simParams_->getPBC(); + int useRF = simParams_->getUseRF(); //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(); - useShape |= i->isShape(); + useLennardJones |= (*i)->isLennardJones(); + useElectrostatic |= (*i)->isElectrostatic(); + useEAM |= (*i)->isEAM(); + useCharge |= (*i)->isCharge(); + useDirectional |= (*i)->isDirectional(); + useDipole |= (*i)->isDipole(); + useGayBerne |= (*i)->isGayBerne(); + useSticky |= (*i)->isSticky(); + useShape |= (*i)->isShape(); } if (useSticky || useDipole || useGayBerne || useShape) { @@ -488,7 +541,18 @@ void SimInfo::setupSimType() { fInfo_.SIM_uses_RF = useRF; if( fInfo_.SIM_uses_Dipoles && fInfo_.SIM_uses_RF) { - fInfo_.dielect = dielectric; + + 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(); + } + } else { fInfo_.dielect = 0.0; } @@ -510,11 +574,11 @@ void SimInfo::setupFortranSim() { //calculate mass ratio of cutoff group std::vector mfact; - typename SimInfo::MoleculeIterator mi; + SimInfo::MoleculeIterator mi; Molecule* mol; - typename Molecule::CutoffGroupIterator ci; + Molecule::CutoffGroupIterator ci; CutoffGroup* cg; - typename Molecule::AtomIterator ai; + Molecule::AtomIterator ai; Atom* atom; double totalMass; @@ -544,13 +608,22 @@ void SimInfo::setupFortranSim() { } } + //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 //gloalExcludes and molMembershipArray should go away (They are never used) //why the hell fortran need to know molecule? //OOPSE = Object-Obfuscated Parallel Simulation Engine - - setFortranSim( &fInfo_, &nGlobalAtoms_, &nAtoms_, &identArray[0], &nExclude, exclude_->getExcludeList(), - &nGlobalExcludes, globalExcludes, molMembershipArray, + 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 ){ @@ -576,9 +649,9 @@ void SimInfo::setupFortranParallel() { //SimInfo is responsible for creating localToGlobalAtomIndex and localToGlobalGroupIndex std::vector localToGlobalAtomIndex(getNAtoms(), 0); std::vector localToGlobalCutoffGroupIndex; - typename SimInfo::MoleculeIterator mi; - typename Molecule::AtomIterator ai; - typename Molecule::CutoffGroupIterator ci; + SimInfo::MoleculeIterator mi; + Molecule::AtomIterator ai; + Molecule::CutoffGroupIterator ci; Molecule* mol; Atom* atom; CutoffGroup* cg; @@ -607,11 +680,11 @@ void SimInfo::setupFortranParallel() { parallelData.nGroupsGlobal = getNGlobalCutoffGroups(); parallelData.nGroupsLocal = getNCutoffGroups(); parallelData.myNode = worldRank; - MPI_Comm_size(MPI_COMM_WORLD, &(parallelData->nProcessors)); + MPI_Comm_size(MPI_COMM_WORLD, &(parallelData.nProcessors)); //pass mpiSimData struct and index arrays to fortran - setFsimParallel(parallelData, &(parallelData->nAtomsLocal), - &localToGlobalAtomIndex[0], &(parallelData->nGroupsLocal), + setFsimParallel(¶llelData, &(parallelData.nAtomsLocal), + &localToGlobalAtomIndex[0], &(parallelData.nGroupsLocal), &localToGlobalCutoffGroupIndex[0], &isError); if (isError) { @@ -629,6 +702,84 @@ void SimInfo::addProperty(GenericData* genData) { #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::setupCutoff() { + double rcut_; //cutoff radius + double rsw_; //switching radius + + if (fInfo_.SIM_uses_Charges | fInfo_.SIM_uses_Dipoles | fInfo_.SIM_uses_RF) { + + if (!simParams_->haveRcut()){ + 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_->getRcut(); + } + + if (!simParams_->haveRsw()){ + sprintf(painCave.errMsg, + "SimCreator Warning: No value was set for switchingRadius.\n" + "\tOOPSE will use a default value of\n" + "\t0.95 * cutoffRadius for the switchingRadius\n"); + painCave.isFatal = 0; + simError(); + rsw_ = 0.95 * rcut_; + } else{ + rsw_ = simParams_->getRsw(); + } + + } 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_->haveRcut()) { + rcut_ = simParams_->getRcut(); + } else { + //set cutoff radius to the maximum cutoff radius based on atom types in the whole system + rcut_ = calcMaxCutoffRadius(); + } + + if (simParams_->haveRsw()) { + rsw_ = simParams_->getRsw(); + } else { + rsw_ = rcut_; + } + + } + + double rnblist = rcut_ + 1; // skin of neighbor list + + //Pass these cutoff radius etc. to fortran. This function should be called once and only once + notifyFortranCutoffs(&rcut_, &rsw_, &rnblist); +} + void SimInfo::addProperty(GenericData* genData) { properties_.addProperty(genData); } @@ -653,10 +804,85 @@ GenericData* SimInfo::getPropertyByName(const std::str return properties_.getPropertyByName(propName); } +void SimInfo::setSnapshotManager(SnapshotManager* sman) { + sman_ = sman; -std::ostream& operator <<(ostream& o, SimInfo& info) { + Molecule* mol; + RigidBody* rb; + Atom* atom; + SimInfo::MoleculeIterator mi; + Molecule::RigidBodyIterator rbIter; + Molecule::AtomIterator atomIter;; + + for (mol = beginMolecule(mi); mol != NULL; mol = nextMolecule(mi)) { + + 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)) { + rb->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(); + } + +#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; } }//end namespace oopse +