--- branches/new_design/OOPSE-2.0/src/brains/SimInfo.cpp 2004/11/13 05:08:12 1738 +++ branches/new_design/OOPSE-2.0/src/brains/SimInfo.cpp 2004/12/02 16:04:19 1832 @@ -31,57 +31,92 @@ */ #include - +#include #include "brains/SimInfo.hpp" +#include "primitives/Molecule.hpp" +#include "UseTheForce/doForces_interface.h" +#include "UseTheForce/notifyCutoffs_interface.h" #include "utils/MemoryUtils.hpp" +#include "utils/simError.h" namespace oopse { -SimInfo::SimInfo(const std::vector >& molStampPairs, +SimInfo::SimInfo(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), + nCutoffGroups_(0), nConstraints_(0), nZconstraint_(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; + int nCutoffAtoms; // number of atoms belong to cutoff groups + int nGroups; //total cutoff groups defined in meta-data file CutoffGroupStamp* cgStamp; - int nAtomsIngroups; - int nCutoffGroupsInStamp; - + int nAtomsInGroups; + int nCutoffGroupsInStamp; + + RigidBodyStamp* rbStamp; + int nAtomsInRigidBodies; + int nRigidBodiesInStamp; + int nRigidAtoms; + int nRigidBodies; + nGlobalAtoms_ = 0; - ngroups = 0; + nGroups = 0; + nCutoffAtoms = 0; + nRigidBodies = 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; - - nAtomsIngroups = 0; + + + //calculate atoms in cutoff groups + 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 + nAtomsInRigidBodies = 0; + nRigidBodiesInStamp = molStamp->getNCutoffGroups(); + + for (int j=0; j < nRigidBodiesInStamp; j++) { + rbStamp = molStamp->getRigidBody(j); + nRigidBodiesInStamp += rbStamp->getNMembers(); + } + + nRigidBodies += 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; + //every free atom (atom does not belong to rigid bodies) is a rigid body + //therefore the total number of cutoff groups 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 + nRigidBodies; + //initialize globalGroupMembership_, every element of this array will be 0 globalGroupMembership_.insert(globalGroupMembership_.end(), nGlobalAtoms_, 0); @@ -111,7 +146,7 @@ bool SimInfo::addMolecule(Molecule* mol) { i = molecules_.find(mol->getGlobalIndex()); 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(); @@ -206,7 +241,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_; } @@ -257,7 +292,7 @@ void SimInfo::calcNdfTrans() { ndfTrans_ = ndfTrans_local; #endif - ndfTrans_ = ndfTrans_ - 3 - nZconstraints_; + ndfTrans_ = ndfTrans_ - 3 - nZconstraint_; } @@ -289,7 +324,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(); @@ -334,7 +369,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(); @@ -358,7 +393,7 @@ void SimInfo::addMoleculeStamp(MoleculeStamp* molStamp curStampId = molStampIds_.size(); moleculeStamps_.push_back(molStamp); - molStampIds_.insert(molStampIds_.end(), nmol, curStampId) + molStampIds_.insert(molStampIds_.end(), nmol, curStampId); } void SimInfo::update() { @@ -371,19 +406,20 @@ void SimInfo::update() { 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(); - //notify fortran whether reaction field is used or not. It is deprecated now - //int isError = 0; - //initFortranFF( &useReactionField, &isError ); - - //if(isError){ - // sprintf( painCave.errMsg, - // "SimCreator::initFortran() error: There was an error initializing the forceField in fortran.\n" ); - // painCave.isFatal = 1; - // simError(); - //} - calcNdf(); calcNdfRaw(); calcNdfTrans(); @@ -392,9 +428,9 @@ std::set SimInfo::getUniqueAtomTypes() { } 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; @@ -432,15 +468,15 @@ void SimInfo::setupSimType() { //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) { @@ -506,7 +542,18 @@ void SimInfo::setupSimType() { fInfo_.SIM_uses_RF = useRF; if( fInfo_.SIM_uses_Dipoles && fInfo_.SIM_uses_RF) { - fInfo_.dielect = dielectric; + + if (globals_->haveDielectric()) { + fInfo_.dielect = globals_->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; } @@ -528,11 +575,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; @@ -573,9 +620,11 @@ void SimInfo::setupFortranSim() { //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 ){ @@ -601,9 +650,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; @@ -657,8 +706,8 @@ double SimInfo::calcMaxCutoffRadius() { double SimInfo::calcMaxCutoffRadius() { - std::vector atomTypes; - std::vector::iterator i; + std::set atomTypes; + std::set::iterator i; std::vector cutoffRadius; //get the unique atom types @@ -669,7 +718,7 @@ double SimInfo::calcMaxCutoffRadius() { cutoffRadius.push_back(forceField_->getRcutFromAtomType(*i)); } - double maxCutoffRadius = std::max_element(cutoffRadius.begin(), cutoffRadius.end()); + double maxCutoffRadius = *(std::max_element(cutoffRadius.begin(), cutoffRadius.end())); #ifdef IS_MPI //pick the max cutoff radius among the processors #endif @@ -719,7 +768,7 @@ void SimInfo::setupCutoff() { } if (globals_->haveRsw()) { - rsw_ = globals_->getRsw() + rsw_ = globals_->getRsw(); } else { rsw_ = rcut_; } @@ -756,9 +805,31 @@ 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_); + } + } + +} +std::ostream& operator <<(std::ostream& o, SimInfo& info) { + return o; }