--- branches/new_design/OOPSE-2.0/src/brains/SimInfo.cpp 2004/11/12 17:40:03 1735 +++ branches/new_design/OOPSE-2.0/src/brains/SimInfo.cpp 2004/11/30 22:43:51 1807 @@ -33,53 +33,89 @@ #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); @@ -159,12 +195,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 +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_; } @@ -256,7 +292,7 @@ void SimInfo::calcNdfTrans() { ndfTrans_ = ndfTrans_local; #endif - ndfTrans_ = ndfTrans_ - 3 - nZconstraints_; + ndfTrans_ = ndfTrans_ - 3 - nZconstraint_; } @@ -288,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(); @@ -333,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(); @@ -357,26 +393,44 @@ 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() { + 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; @@ -414,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) { @@ -488,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; } @@ -510,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; @@ -555,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 ){ @@ -583,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; @@ -639,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 @@ -651,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 @@ -659,6 +726,61 @@ void SimInfo::addProperty(GenericData* genData) { 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 (!globals_->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_ = globals_->getRcut(); + } + + if (!globals_->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_ = globals_->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 (globals_->haveRcut()) { + rcut_ = globals_->getRcut(); + } else { + //set cutoff radius to the maximum cutoff radius based on atom types in the whole system + rcut_ = calcMaxCutoffRadius(); + } + + if (globals_->haveRsw()) { + rsw_ = globals_->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); } @@ -683,7 +805,29 @@ GenericData* SimInfo::getPropertyByName(const std::str return properties_.getPropertyByName(propName); } +void SimInfo::setSnapshotManager(SnapshotManager* sman) { + sman_ = sman; + 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 <<(ostream& o, SimInfo& info) { return o;