--- trunk/src/brains/SimInfo.cpp 2005/03/09 17:30:29 413 +++ trunk/src/brains/SimInfo.cpp 2005/09/16 21:07:45 608 @@ -1,4 +1,4 @@ - /* +/* * Copyright (c) 2005 The University of Notre Dame. All Rights Reserved. * * The University of Notre Dame grants you ("Licensee") a @@ -52,6 +52,8 @@ #include "brains/SimInfo.hpp" #include "math/Vector3.hpp" #include "primitives/Molecule.hpp" +#include "UseTheForce/fCutoffPolicy.h" +#include "UseTheForce/DarkSide/fElectrostaticSummationMethod.h" #include "UseTheForce/doForces_interface.h" #include "UseTheForce/notifyCutoffs_interface.h" #include "utils/MemoryUtils.hpp" @@ -65,27 +67,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) { + 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 +102,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,50 +114,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 -} + } -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_; -} + } -int SimInfo::getNGlobalConstraints() { + int SimInfo::getNGlobalConstraints() { int nGlobalConstraints; #ifdef IS_MPI MPI_Allreduce(&nConstraints_, &nGlobalConstraints, 1, MPI_INT, MPI_SUM, @@ -164,76 +165,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(); + 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); + addExcludePairs(mol); - return true; + return true; } else { - return false; + return false; } -} + } -bool SimInfo::removeMolecule(Molecule* mol) { + 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(); - - removeExcludePairs(mol); - molecules_.erase(mol->getGlobalIndex()); + nAtoms_ -= mol->getNAtoms(); + nBonds_ -= mol->getNBonds(); + nBends_ -= mol->getNBends(); + nTorsions_ -= mol->getNTorsions(); + nRigidBodies_ -= mol->getNRigidBodies(); + nIntegrableObjects_ -= mol->getNIntegrableObjects(); + nCutoffGroups_ -= mol->getNCutoffGroups(); + nConstraints_ -= mol->getNConstraintPairs(); - delete mol; + removeExcludePairs(mol); + molecules_.erase(mol->getGlobalIndex()); + + 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; @@ -243,20 +244,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 @@ -272,9 +273,9 @@ void SimInfo::calcNdf() { // entire system: ndf_ = ndf_ - 3 - nZconstraint_; -} + } -void SimInfo::calcNdfRaw() { + void SimInfo::calcNdfRaw() { int ndfRaw_local; MoleculeIterator i; @@ -286,20 +287,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 @@ -307,9 +308,9 @@ void SimInfo::calcNdfRaw() { #else ndfRaw_ = ndfRaw_local; #endif -} + } -void SimInfo::calcNdfTrans() { + void SimInfo::calcNdfTrans() { int ndfTrans_local; ndfTrans_local = 3 * nIntegrableObjects_ - nConstraints_; @@ -323,9 +324,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; @@ -338,39 +339,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; @@ -383,39 +396,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 @@ -423,9 +449,9 @@ void SimInfo::addMoleculeStamp(MoleculeStamp* molStamp moleculeStamps_.push_back(molStamp); molStampIds_.insert(molStampIds_.end(), nmol, curStampId); -} + } -void SimInfo::update() { + void SimInfo::update() { setupSimType(); @@ -438,12 +464,14 @@ void SimInfo::update() { //setup fortran force field /** @deprecate */ int isError = 0; - initFortranFF( &fInfo_.SIM_uses_RF , &isError ); + + setupElectrostaticSummationMethod( 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(); } @@ -454,9 +482,9 @@ void SimInfo::update() { calcNdfTrans(); fortranInitialized_ = true; -} + } -std::set SimInfo::getUniqueAtomTypes() { + std::set SimInfo::getUniqueAtomTypes() { SimInfo::MoleculeIterator mi; Molecule* mol; Molecule::AtomIterator ai; @@ -465,16 +493,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(); @@ -487,33 +515,34 @@ 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; int useElectrostatics = 0; //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(); + 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 @@ -540,6 +569,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); @@ -552,9 +584,6 @@ void SimInfo::setupSimType() { temp = useFLARB; MPI_Allreduce(&temp, &useFLARB, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD); - temp = useRF; - MPI_Allreduce(&temp, &useRF, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD); - #endif fInfo_.SIM_uses_PBC = usePBC; @@ -564,32 +593,32 @@ 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; fInfo_.SIM_uses_FLARB = useFLARB; - fInfo_.SIM_uses_RF = useRF; 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; @@ -599,7 +628,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 @@ -616,14 +645,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 !!!) @@ -633,48 +662,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); @@ -690,15 +716,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); + } } @@ -718,21 +744,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; @@ -744,7 +770,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())); @@ -753,94 +779,158 @@ 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); -} + + int cp = TRADITIONAL_CUTOFF_POLICY; + if (simParams_->haveCutoffPolicy()) { + std::string myPolicy = simParams_->getCutoffPolicy(); + if (myPolicy == "MIX") { + cp = MIX_CUTOFF_POLICY; + } else { + if (myPolicy == "MAX") { + cp = MAX_CUTOFF_POLICY; + } else { + if (myPolicy == "TRADITIONAL") { + cp = TRADITIONAL_CUTOFF_POLICY; + } else { + // throw error + sprintf( painCave.errMsg, + "SimInfo error: Unknown cutoffPolicy. (Input file specified %s .)\n\tcutoffPolicy must be one of: \"Mix\", \"Max\", or \"Traditional\".", myPolicy.c_str() ); + painCave.isFatal = 1; + simError(); + } + } + } + } + notifyFortranCutoffs(&rcut_, &rsw_, &rnblist, &cp); + } -void SimInfo::addProperty(GenericData* genData) { + void SimInfo::setupElectrostaticSummationMethod( int isError ) { + + int errorOut; + int esm = NONE; + double alphaVal; + + errorOut = isError; + + if (simParams_->haveElectrostaticSummationMethod()) { + std::string myMethod = simParams_->getElectrostaticSummationMethod(); + if (myMethod == "NONE") { + esm = NONE; + } else { + if (myMethod == "UNDAMPED_WOLF") { + esm = UNDAMPED_WOLF; + } else { + if (myMethod == "DAMPED_WOLF") { + esm = DAMPED_WOLF; + if (!simParams_->haveDampingAlpha()) { + //throw error + sprintf( painCave.errMsg, + "SimInfo warning: dampingAlpha was not specified in the input file. A default value of %f (1/ang) will be used for the Damped Wolf Method.", simParams_->getDampingAlpha()); + painCave.isFatal = 0; + simError(); + } + alphaVal = simParams_->getDampingAlpha(); + } else { + if (myMethod == "REACTION_FIELD") { + esm = REACTION_FIELD; + } else { + // throw error + sprintf( painCave.errMsg, + "SimInfo error: Unknown electrostaticSummationMethod. (Input file specified %s .)\n\telectrostaticSummationMethod must be one of: \"none\", \"undamped_wolf\", \"damped_wolf\", or \"reaction_field\".", myMethod.c_str() ); + painCave.isFatal = 1; + simError(); + } + } + } + } + } + initFortranFF( &esm, &alphaVal, &errorOut ); + } + + 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; @@ -852,18 +942,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; @@ -872,9 +962,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 @@ -887,9 +977,9 @@ Vector3d SimInfo::getComVel(){ comVel /= totalMass; return comVel; -} + } -Vector3d SimInfo::getCom(){ + Vector3d SimInfo::getCom(){ SimInfo::MoleculeIterator i; Molecule* mol; @@ -897,9 +987,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 @@ -913,12 +1003,154 @@ Vector3d SimInfo::getCom(){ return com; -} + } -std::ostream& operator <<(std::ostream& o, SimInfo& info) { + std::ostream& operator <<(std::ostream& o, SimInfo& info) { return o; -} + } + + + /* + Returns center of mass and center of mass velocity in one function call. + */ + + void SimInfo::getComAll(Vector3d &com, Vector3d &comVel){ + SimInfo::MoleculeIterator i; + Molecule* mol; + + + double totalMass = 0.0; + + for (mol = beginMolecule(i); mol != NULL; mol = nextMolecule(i)) { + double mass = mol->getMass(); + totalMass += mass; + com += mass * mol->getCom(); + comVel += mass * mol->getComVel(); + } + +#ifdef IS_MPI + double tmpMass = totalMass; + Vector3d tmpCom(com); + Vector3d tmpComVel(comVel); + 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); + MPI_Allreduce(tmpComVel.getArrayPointer(), comVel.getArrayPointer(),3,MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD); +#endif + + com /= totalMass; + comVel /= totalMass; + } + + /* + Return intertia tensor for entire system and angular momentum Vector. + + + [ Ixx -Ixy -Ixz ] + J =| -Iyx Iyy -Iyz | + [ -Izx -Iyz Izz ] + */ + + void SimInfo::getInertiaTensor(Mat3x3d &inertiaTensor, Vector3d &angularMomentum){ + + + double xx = 0.0; + double yy = 0.0; + double zz = 0.0; + double xy = 0.0; + double xz = 0.0; + double yz = 0.0; + Vector3d com(0.0); + Vector3d comVel(0.0); + + getComAll(com, comVel); + + SimInfo::MoleculeIterator i; + Molecule* mol; + + Vector3d thisq(0.0); + Vector3d thisv(0.0); + + double thisMass = 0.0; + + + + + for (mol = beginMolecule(i); mol != NULL; mol = nextMolecule(i)) { + + thisq = mol->getCom()-com; + thisv = mol->getComVel()-comVel; + thisMass = mol->getMass(); + // Compute moment of intertia coefficients. + xx += thisq[0]*thisq[0]*thisMass; + yy += thisq[1]*thisq[1]*thisMass; + zz += thisq[2]*thisq[2]*thisMass; + + // compute products of intertia + xy += thisq[0]*thisq[1]*thisMass; + xz += thisq[0]*thisq[2]*thisMass; + yz += thisq[1]*thisq[2]*thisMass; + + angularMomentum += cross( thisq, thisv ) * thisMass; + + } + + + inertiaTensor(0,0) = yy + zz; + inertiaTensor(0,1) = -xy; + inertiaTensor(0,2) = -xz; + inertiaTensor(1,0) = -xy; + inertiaTensor(1,1) = xx + zz; + inertiaTensor(1,2) = -yz; + inertiaTensor(2,0) = -xz; + inertiaTensor(2,1) = -yz; + inertiaTensor(2,2) = xx + yy; + +#ifdef IS_MPI + Mat3x3d tmpI(inertiaTensor); + Vector3d tmpAngMom; + MPI_Allreduce(tmpI.getArrayPointer(), inertiaTensor.getArrayPointer(),9,MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD); + MPI_Allreduce(tmpAngMom.getArrayPointer(), angularMomentum.getArrayPointer(),3,MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD); +#endif + + return; + } + + //Returns the angular momentum of the system + Vector3d SimInfo::getAngularMomentum(){ + + Vector3d com(0.0); + Vector3d comVel(0.0); + Vector3d angularMomentum(0.0); + + getComAll(com,comVel); + + SimInfo::MoleculeIterator i; + Molecule* mol; + + Vector3d thisr(0.0); + Vector3d thisp(0.0); + + double thisMass; + + for (mol = beginMolecule(i); mol != NULL; mol = nextMolecule(i)) { + thisMass = mol->getMass(); + thisr = mol->getCom()-com; + thisp = (mol->getComVel()-comVel)*thisMass; + + angularMomentum += cross( thisr, thisp ); + + } + +#ifdef IS_MPI + Vector3d tmpAngMom; + MPI_Allreduce(tmpAngMom.getArrayPointer(), angularMomentum.getArrayPointer(),3,MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD); +#endif + + return angularMomentum; + } + + }//end namespace oopse