--- trunk/src/brains/SimInfo.cpp 2005/03/07 22:39:33 398 +++ trunk/src/brains/SimInfo.cpp 2005/05/05 14:47:35 523 @@ -1,4 +1,4 @@ - /* +/* * Copyright (c) 2005 The University of Notre Dame. All Rights Reserved. * * The University of Notre Dame grants you ("Licensee") a @@ -65,27 +65,27 @@ namespace oopse { namespace oopse { -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), selectMan_(NULL) { + SimInfo::SimInfo(MakeStamps* stamps, std::vector >& molStampPairs, + ForceField* ff, Globals* simParams) : + stamps_(stamps), 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; - MoleculeStamp* molStamp; - int nMolWithSameStamp; - 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; + std::vector >::iterator i; + MoleculeStamp* molStamp; + int nMolWithSameStamp; + 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) { + for (i = molStampPairs.begin(); i !=molStampPairs.end(); ++i) { molStamp = i->first; nMolWithSameStamp = i->second; @@ -100,8 +100,8 @@ SimInfo::SimInfo(std::vectorgetNCutoffGroups(); for (int j=0; j < nCutoffGroupsInStamp; j++) { - cgStamp = molStamp->getCutoffGroup(j); - nAtomsInGroups += cgStamp->getNMembers(); + cgStamp = molStamp->getCutoffGroup(j); + nAtomsInGroups += cgStamp->getNMembers(); } nGroups += nCutoffGroupsInStamp * nMolWithSameStamp; @@ -112,53 +112,49 @@ SimInfo::SimInfo(std::vectorgetNRigidBodies(); for (int j=0; j < nRigidBodiesInStamp; j++) { - rbStamp = molStamp->getRigidBody(j); - nAtomsInRigidBodies += rbStamp->getNMembers(); + 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; + //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; - //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_; + //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(); + nGlobalMols_ = molStampIds_.size(); #ifdef IS_MPI - molToProcMap_.resize(nGlobalMols_); + molToProcMap_.resize(nGlobalMols_); #endif - selectMan_ = new SelectionManager(this); - selectMan_->selectAll(); -} + } -SimInfo::~SimInfo() { + SimInfo::~SimInfo() { std::map::iterator i; for (i = molecules_.begin(); i != molecules_.end(); ++i) { - delete i->second; + delete i->second; } molecules_.clear(); - - MemoryUtils::deletePointers(moleculeStamps_); - + + delete stamps_; delete sman_; delete simParams_; delete forceField_; - delete selectMan_; -} + } -int SimInfo::getNGlobalConstraints() { + int SimInfo::getNGlobalConstraints() { int nGlobalConstraints; #ifdef IS_MPI MPI_Allreduce(&nConstraints_, &nGlobalConstraints, 1, MPI_INT, MPI_SUM, @@ -167,76 +163,76 @@ int SimInfo::getNGlobalConstraints() { nGlobalConstraints = nConstraints_; #endif return nGlobalConstraints; -} + } -bool SimInfo::addMolecule(Molecule* mol) { + bool SimInfo::addMolecule(Molecule* mol) { MoleculeIterator i; i = molecules_.find(mol->getGlobalIndex()); if (i == molecules_.end() ) { - molecules_.insert(std::make_pair(mol->getGlobalIndex(), mol)); + molecules_.insert(std::make_pair(mol->getGlobalIndex(), mol)); - nAtoms_ += mol->getNAtoms(); - nBonds_ += mol->getNBonds(); - nBends_ += mol->getNBends(); - nTorsions_ += mol->getNTorsions(); - nRigidBodies_ += mol->getNRigidBodies(); - nIntegrableObjects_ += mol->getNIntegrableObjects(); - nCutoffGroups_ += mol->getNCutoffGroups(); - nConstraints_ += mol->getNConstraintPairs(); - - addExcludePairs(mol); - - return true; - } else { - return false; - } -} + nAtoms_ += mol->getNAtoms(); + nBonds_ += mol->getNBonds(); + nBends_ += mol->getNBends(); + nTorsions_ += mol->getNTorsions(); + nRigidBodies_ += mol->getNRigidBodies(); + nIntegrableObjects_ += mol->getNIntegrableObjects(); + nCutoffGroups_ += mol->getNCutoffGroups(); + nConstraints_ += mol->getNConstraintPairs(); -bool SimInfo::removeMolecule(Molecule* mol) { + addExcludePairs(mol); + + return true; + } else { + return false; + } + } + + bool SimInfo::removeMolecule(Molecule* mol) { MoleculeIterator i; i = molecules_.find(mol->getGlobalIndex()); if (i != molecules_.end() ) { - assert(mol == i->second); + assert(mol == i->second); - nAtoms_ -= mol->getNAtoms(); - nBonds_ -= mol->getNBonds(); - nBends_ -= mol->getNBends(); - nTorsions_ -= mol->getNTorsions(); - nRigidBodies_ -= mol->getNRigidBodies(); - nIntegrableObjects_ -= mol->getNIntegrableObjects(); - nCutoffGroups_ -= mol->getNCutoffGroups(); - nConstraints_ -= mol->getNConstraintPairs(); + nAtoms_ -= mol->getNAtoms(); + nBonds_ -= mol->getNBonds(); + nBends_ -= mol->getNBends(); + nTorsions_ -= mol->getNTorsions(); + nRigidBodies_ -= mol->getNRigidBodies(); + nIntegrableObjects_ -= mol->getNIntegrableObjects(); + nCutoffGroups_ -= mol->getNCutoffGroups(); + nConstraints_ -= mol->getNConstraintPairs(); - removeExcludePairs(mol); - molecules_.erase(mol->getGlobalIndex()); + removeExcludePairs(mol); + molecules_.erase(mol->getGlobalIndex()); - delete mol; + delete mol; - return true; + return true; } else { - return false; + return false; } -} + } -Molecule* SimInfo::beginMolecule(MoleculeIterator& i) { + Molecule* SimInfo::beginMolecule(MoleculeIterator& i) { i = molecules_.begin(); return i == molecules_.end() ? NULL : i->second; -} + } -Molecule* SimInfo::nextMolecule(MoleculeIterator& i) { + Molecule* SimInfo::nextMolecule(MoleculeIterator& i) { ++i; return i == molecules_.end() ? NULL : i->second; -} + } -void SimInfo::calcNdf() { + void SimInfo::calcNdf() { int ndf_local; MoleculeIterator i; std::vector::iterator j; @@ -246,20 +242,20 @@ void SimInfo::calcNdf() { ndf_local = 0; for (mol = beginMolecule(i); mol != NULL; mol = nextMolecule(i)) { - for (integrableObject = mol->beginIntegrableObject(j); integrableObject != NULL; - integrableObject = mol->nextIntegrableObject(j)) { + for (integrableObject = mol->beginIntegrableObject(j); integrableObject != NULL; + integrableObject = mol->nextIntegrableObject(j)) { - ndf_local += 3; + ndf_local += 3; - if (integrableObject->isDirectional()) { - if (integrableObject->isLinear()) { - ndf_local += 2; - } else { - ndf_local += 3; - } - } + if (integrableObject->isDirectional()) { + if (integrableObject->isLinear()) { + ndf_local += 2; + } else { + ndf_local += 3; + } + } - }//end for (integrableObject) + }//end for (integrableObject) }// end for (mol) // n_constraints is local, so subtract them on each processor @@ -275,9 +271,9 @@ void SimInfo::calcNdf() { // entire system: ndf_ = ndf_ - 3 - nZconstraint_; -} + } -void SimInfo::calcNdfRaw() { + void SimInfo::calcNdfRaw() { int ndfRaw_local; MoleculeIterator i; @@ -289,20 +285,20 @@ void SimInfo::calcNdfRaw() { ndfRaw_local = 0; for (mol = beginMolecule(i); mol != NULL; mol = nextMolecule(i)) { - for (integrableObject = mol->beginIntegrableObject(j); integrableObject != NULL; - integrableObject = mol->nextIntegrableObject(j)) { + for (integrableObject = mol->beginIntegrableObject(j); integrableObject != NULL; + integrableObject = mol->nextIntegrableObject(j)) { - ndfRaw_local += 3; + ndfRaw_local += 3; - if (integrableObject->isDirectional()) { - if (integrableObject->isLinear()) { - ndfRaw_local += 2; - } else { - ndfRaw_local += 3; - } - } + if (integrableObject->isDirectional()) { + if (integrableObject->isLinear()) { + ndfRaw_local += 2; + } else { + ndfRaw_local += 3; + } + } - } + } } #ifdef IS_MPI @@ -310,9 +306,9 @@ void SimInfo::calcNdfRaw() { #else ndfRaw_ = ndfRaw_local; #endif -} + } -void SimInfo::calcNdfTrans() { + void SimInfo::calcNdfTrans() { int ndfTrans_local; ndfTrans_local = 3 * nIntegrableObjects_ - nConstraints_; @@ -326,9 +322,9 @@ void SimInfo::calcNdfTrans() { ndfTrans_ = ndfTrans_ - 3 - nZconstraint_; -} + } -void SimInfo::addExcludePairs(Molecule* mol) { + void SimInfo::addExcludePairs(Molecule* mol) { std::vector::iterator bondIter; std::vector::iterator bendIter; std::vector::iterator torsionIter; @@ -341,39 +337,51 @@ void SimInfo::addExcludePairs(Molecule* mol) { int d; for (bond= mol->beginBond(bondIter); bond != NULL; bond = mol->nextBond(bondIter)) { - a = bond->getAtomA()->getGlobalIndex(); - b = bond->getAtomB()->getGlobalIndex(); - exclude_.addPair(a, b); + a = bond->getAtomA()->getGlobalIndex(); + b = bond->getAtomB()->getGlobalIndex(); + exclude_.addPair(a, b); } for (bend= mol->beginBend(bendIter); bend != NULL; bend = mol->nextBend(bendIter)) { - a = bend->getAtomA()->getGlobalIndex(); - b = bend->getAtomB()->getGlobalIndex(); - c = bend->getAtomC()->getGlobalIndex(); + a = bend->getAtomA()->getGlobalIndex(); + b = bend->getAtomB()->getGlobalIndex(); + c = bend->getAtomC()->getGlobalIndex(); - exclude_.addPair(a, b); - exclude_.addPair(a, c); - exclude_.addPair(b, c); + exclude_.addPair(a, b); + exclude_.addPair(a, c); + exclude_.addPair(b, c); } for (torsion= mol->beginTorsion(torsionIter); torsion != NULL; torsion = mol->nextTorsion(torsionIter)) { - a = torsion->getAtomA()->getGlobalIndex(); - b = torsion->getAtomB()->getGlobalIndex(); - c = torsion->getAtomC()->getGlobalIndex(); - d = torsion->getAtomD()->getGlobalIndex(); + a = torsion->getAtomA()->getGlobalIndex(); + b = torsion->getAtomB()->getGlobalIndex(); + c = torsion->getAtomC()->getGlobalIndex(); + d = torsion->getAtomD()->getGlobalIndex(); - exclude_.addPair(a, b); - exclude_.addPair(a, c); - exclude_.addPair(a, d); - exclude_.addPair(b, c); - exclude_.addPair(b, d); - exclude_.addPair(c, d); + exclude_.addPair(a, b); + exclude_.addPair(a, c); + exclude_.addPair(a, d); + exclude_.addPair(b, c); + exclude_.addPair(b, d); + exclude_.addPair(c, d); } - -} + Molecule::RigidBodyIterator rbIter; + RigidBody* rb; + 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) { + for (int j = i + 1; j < atoms.size(); ++j) { + a = atoms[i]->getGlobalIndex(); + b = atoms[j]->getGlobalIndex(); + exclude_.addPair(a, b); + } + } + } -void SimInfo::removeExcludePairs(Molecule* mol) { + } + + void SimInfo::removeExcludePairs(Molecule* mol) { std::vector::iterator bondIter; std::vector::iterator bendIter; std::vector::iterator torsionIter; @@ -386,39 +394,52 @@ void SimInfo::removeExcludePairs(Molecule* mol) { int d; for (bond= mol->beginBond(bondIter); bond != NULL; bond = mol->nextBond(bondIter)) { - a = bond->getAtomA()->getGlobalIndex(); - b = bond->getAtomB()->getGlobalIndex(); - exclude_.removePair(a, b); + a = bond->getAtomA()->getGlobalIndex(); + b = bond->getAtomB()->getGlobalIndex(); + exclude_.removePair(a, b); } for (bend= mol->beginBend(bendIter); bend != NULL; bend = mol->nextBend(bendIter)) { - a = bend->getAtomA()->getGlobalIndex(); - b = bend->getAtomB()->getGlobalIndex(); - c = bend->getAtomC()->getGlobalIndex(); + a = bend->getAtomA()->getGlobalIndex(); + b = bend->getAtomB()->getGlobalIndex(); + c = bend->getAtomC()->getGlobalIndex(); - exclude_.removePair(a, b); - exclude_.removePair(a, c); - exclude_.removePair(b, c); + exclude_.removePair(a, b); + exclude_.removePair(a, c); + exclude_.removePair(b, c); } for (torsion= mol->beginTorsion(torsionIter); torsion != NULL; torsion = mol->nextTorsion(torsionIter)) { - a = torsion->getAtomA()->getGlobalIndex(); - b = torsion->getAtomB()->getGlobalIndex(); - c = torsion->getAtomC()->getGlobalIndex(); - d = torsion->getAtomD()->getGlobalIndex(); + a = torsion->getAtomA()->getGlobalIndex(); + b = torsion->getAtomB()->getGlobalIndex(); + c = torsion->getAtomC()->getGlobalIndex(); + d = torsion->getAtomD()->getGlobalIndex(); - exclude_.removePair(a, b); - exclude_.removePair(a, c); - exclude_.removePair(a, d); - exclude_.removePair(b, c); - exclude_.removePair(b, d); - exclude_.removePair(c, d); + exclude_.removePair(a, b); + exclude_.removePair(a, c); + exclude_.removePair(a, d); + exclude_.removePair(b, c); + exclude_.removePair(b, d); + exclude_.removePair(c, d); } -} + Molecule::RigidBodyIterator rbIter; + RigidBody* rb; + 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) { + for (int j = i + 1; j < atoms.size(); ++j) { + a = atoms[i]->getGlobalIndex(); + b = atoms[j]->getGlobalIndex(); + exclude_.removePair(a, b); + } + } + } + } -void SimInfo::addMoleculeStamp(MoleculeStamp* molStamp, int nmol) { + + void SimInfo::addMoleculeStamp(MoleculeStamp* molStamp, int nmol) { int curStampId; //index from 0 @@ -426,9 +447,9 @@ void SimInfo::addMoleculeStamp(MoleculeStamp* molStamp moleculeStamps_.push_back(molStamp); molStampIds_.insert(molStampIds_.end(), nmol, curStampId); -} + } -void SimInfo::update() { + void SimInfo::update() { setupSimType(); @@ -443,10 +464,10 @@ void SimInfo::update() { 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(); + sprintf( painCave.errMsg, + "ForceField error: There was an error initializing the forceField in fortran.\n" ); + painCave.isFatal = 1; + simError(); } @@ -457,9 +478,9 @@ void SimInfo::update() { calcNdfTrans(); fortranInitialized_ = true; -} + } -std::set SimInfo::getUniqueAtomTypes() { + std::set SimInfo::getUniqueAtomTypes() { SimInfo::MoleculeIterator mi; Molecule* mol; Molecule::AtomIterator ai; @@ -468,16 +489,16 @@ std::set SimInfo::getUniqueAtomTypes() { for(mol = beginMolecule(mi); mol != NULL; mol = nextMolecule(mi)) { - for(atom = mol->beginAtom(ai); atom != NULL; atom = mol->nextAtom(ai)) { - atomTypes.insert(atom->getAtomType()); - } + for(atom = mol->beginAtom(ai); atom != NULL; atom = mol->nextAtom(ai)) { + atomTypes.insert(atom->getAtomType()); + } } return atomTypes; -} + } -void SimInfo::setupSimType() { + void SimInfo::setupSimType() { std::set::iterator i; std::set atomTypes; atomTypes = getUniqueAtomTypes(); @@ -490,6 +511,7 @@ void SimInfo::setupSimType() { int useDipole = 0; int useGayBerne = 0; int useSticky = 0; + int useStickyPower = 0; int useShape = 0; int useFLARB = 0; //it is not in AtomType yet int useDirectionalAtom = 0; @@ -500,23 +522,24 @@ 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(); + useStickyPower |= (*i)->isStickyPower(); + useShape |= (*i)->isShape(); } - if (useSticky || useDipole || useGayBerne || useShape) { - useDirectionalAtom = 1; + if (useSticky || useStickyPower || useDipole || useGayBerne || useShape) { + useDirectionalAtom = 1; } if (useCharge || useDipole) { - useElectrostatics = 1; + useElectrostatics = 1; } #ifdef IS_MPI @@ -543,6 +566,9 @@ void SimInfo::setupSimType() { temp = useSticky; MPI_Allreduce(&temp, &useSticky, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD); + temp = useStickyPower; + MPI_Allreduce(&temp, &useStickyPower, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD); + temp = useGayBerne; MPI_Allreduce(&temp, &useGayBerne, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD); @@ -567,6 +593,7 @@ void SimInfo::setupSimType() { fInfo_.SIM_uses_Charges = useCharge; fInfo_.SIM_uses_Dipoles = useDipole; fInfo_.SIM_uses_Sticky = useSticky; + fInfo_.SIM_uses_StickyPower = useStickyPower; fInfo_.SIM_uses_GayBerne = useGayBerne; fInfo_.SIM_uses_EAM = useEAM; fInfo_.SIM_uses_Shapes = useShape; @@ -575,24 +602,24 @@ void SimInfo::setupSimType() { if( fInfo_.SIM_uses_Dipoles && fInfo_.SIM_uses_RF) { - 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(); - } + 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; + fInfo_.dielect = 0.0; } -} + } -void SimInfo::setupFortranSim() { + void SimInfo::setupFortranSim() { int isError; int nExclude; std::vector fortranGlobalGroupMembership; @@ -602,7 +629,7 @@ void SimInfo::setupFortranSim() { //globalGroupMembership_ is filled by SimCreator for (int i = 0; i < nGlobalAtoms_; i++) { - fortranGlobalGroupMembership.push_back(globalGroupMembership_[i] + 1); + fortranGlobalGroupMembership.push_back(globalGroupMembership_[i] + 1); } //calculate mass ratio of cutoff group @@ -619,14 +646,14 @@ void SimInfo::setupFortranSim() { mfact.reserve(getNCutoffGroups()); 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)) { - mfact.push_back(atom->getMass()/totalMass); - } + totalMass = cg->getMass(); + for(atom = cg->beginAtom(ai); atom != NULL; atom = cg->nextAtom(ai)) { + mfact.push_back(atom->getMass()/totalMass); + } - } + } } //fill ident array of local atoms (it is actually ident of AtomType, it is so confusing !!!) @@ -636,48 +663,45 @@ void SimInfo::setupFortranSim() { 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()); - } + for(atom = mol->beginAtom(ai); atom != NULL; atom = mol->nextAtom(ai)) { + identArray.push_back(atom->getIdent()); + } } //fill molMembershipArray //molMembershipArray is filled by SimCreator std::vector molMembershipArray(nGlobalAtoms_); for (int i = 0; i < nGlobalAtoms_; i++) { - molMembershipArray[i] = globalMolMembership_[i] + 1; + 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 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); + &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(); + 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"); + "succesfully sent the simulation information to fortran.\n"); MPIcheckPoint(); #endif // is_mpi -} + } #ifdef IS_MPI -void SimInfo::setupFortranParallel() { + void SimInfo::setupFortranParallel() { //SimInfo is responsible for creating localToGlobalAtomIndex and localToGlobalGroupIndex std::vector localToGlobalAtomIndex(getNAtoms(), 0); @@ -693,15 +717,15 @@ void SimInfo::setupFortranParallel() { 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(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); - } + //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); + } } @@ -721,21 +745,21 @@ void SimInfo::setupFortranParallel() { &localToGlobalCutoffGroupIndex[0], &isError); if (isError) { - sprintf(painCave.errMsg, - "mpiRefresh errror: fortran didn't like something we gave it.\n"); - painCave.isFatal = 1; - simError(); + sprintf(painCave.errMsg, + "mpiRefresh errror: fortran didn't like something we gave it.\n"); + painCave.isFatal = 1; + simError(); } sprintf(checkPointMsg, " mpiRefresh successful.\n"); MPIcheckPoint(); -} + } #endif -double SimInfo::calcMaxCutoffRadius() { + double SimInfo::calcMaxCutoffRadius() { std::set atomTypes; @@ -747,7 +771,7 @@ double SimInfo::calcMaxCutoffRadius() { //query the max cutoff radius among these atom types for (i = atomTypes.begin(); i != atomTypes.end(); ++i) { - cutoffRadius.push_back(forceField_->getRcutFromAtomType(*i)); + cutoffRadius.push_back(forceField_->getRcutFromAtomType(*i)); } double maxCutoffRadius = *(std::max_element(cutoffRadius.begin(), cutoffRadius.end())); @@ -756,94 +780,93 @@ double SimInfo::calcMaxCutoffRadius() { #endif return maxCutoffRadius; -} + } -void SimInfo::getCutoff(double& rcut, double& rsw) { + void SimInfo::getCutoff(double& rcut, double& rsw) { if (fInfo_.SIM_uses_Charges | fInfo_.SIM_uses_Dipoles | fInfo_.SIM_uses_RF) { - if (!simParams_->haveRcut()){ - sprintf(painCave.errMsg, + 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(); - } + painCave.isFatal = 0; + simError(); + rcut = 15.0; + } else{ + rcut = simParams_->getRcut(); + } - if (!simParams_->haveRsw()){ - sprintf(painCave.errMsg, + 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(); - } + 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 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_->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; - } + if (simParams_->haveRsw()) { + rsw = simParams_->getRsw(); + } else { + rsw = rcut; + } } -} + } -void SimInfo::setupCutoff() { + void SimInfo::setupCutoff() { getCutoff(rcut_, rsw_); 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) { + void SimInfo::addProperty(GenericData* genData) { properties_.addProperty(genData); -} + } -void SimInfo::removeProperty(const std::string& propName) { + void SimInfo::removeProperty(const std::string& propName) { properties_.removeProperty(propName); -} + } -void SimInfo::clearProperties() { + void SimInfo::clearProperties() { properties_.clearProperties(); -} + } -std::vector SimInfo::getPropertyNames() { + std::vector SimInfo::getPropertyNames() { return properties_.getPropertyNames(); -} + } -std::vector SimInfo::getProperties() { + std::vector SimInfo::getProperties() { return properties_.getProperties(); -} + } -GenericData* SimInfo::getPropertyByName(const std::string& propName) { + GenericData* SimInfo::getPropertyByName(const std::string& propName) { return properties_.getPropertyByName(propName); -} + } -void SimInfo::setSnapshotManager(SnapshotManager* sman) { - //if (sman_ == sman_) { - // return; - //} - - //delete sman_; + void SimInfo::setSnapshotManager(SnapshotManager* sman) { + if (sman_ == sman) { + return; + } + delete sman_; sman_ = sman; Molecule* mol; @@ -855,18 +878,18 @@ void SimInfo::setSnapshotManager(SnapshotManager* sman for (mol = beginMolecule(mi); mol != NULL; mol = nextMolecule(mi)) { - for (atom = mol->beginAtom(atomIter); atom != NULL; atom = mol->nextAtom(atomIter)) { - atom->setSnapshotManager(sman_); - } + 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_); - } + for (rb = mol->beginRigidBody(rbIter); rb != NULL; rb = mol->nextRigidBody(rbIter)) { + rb->setSnapshotManager(sman_); + } } -} + } -Vector3d SimInfo::getComVel(){ + Vector3d SimInfo::getComVel(){ SimInfo::MoleculeIterator i; Molecule* mol; @@ -875,9 +898,9 @@ Vector3d SimInfo::getComVel(){ for (mol = beginMolecule(i); mol != NULL; mol = nextMolecule(i)) { - double mass = mol->getMass(); - totalMass += mass; - comVel += mass * mol->getComVel(); + double mass = mol->getMass(); + totalMass += mass; + comVel += mass * mol->getComVel(); } #ifdef IS_MPI @@ -890,9 +913,9 @@ Vector3d SimInfo::getComVel(){ comVel /= totalMass; return comVel; -} + } -Vector3d SimInfo::getCom(){ + Vector3d SimInfo::getCom(){ SimInfo::MoleculeIterator i; Molecule* mol; @@ -900,9 +923,9 @@ Vector3d SimInfo::getCom(){ double totalMass = 0.0; for (mol = beginMolecule(i); mol != NULL; mol = nextMolecule(i)) { - double mass = mol->getMass(); - totalMass += mass; - com += mass * mol->getCom(); + double mass = mol->getMass(); + totalMass += mass; + com += mass * mol->getCom(); } #ifdef IS_MPI @@ -916,12 +939,12 @@ Vector3d SimInfo::getCom(){ return com; -} + } -std::ostream& operator <<(std::ostream& o, SimInfo& info) { + std::ostream& operator <<(std::ostream& o, SimInfo& info) { return o; -} + } }//end namespace oopse