--- trunk/src/brains/SimInfo.cpp 2005/05/05 14:47:35 523 +++ trunk/src/brains/SimInfo.cpp 2006/04/25 02:09:01 945 @@ -48,15 +48,23 @@ #include #include +#include #include "brains/SimInfo.hpp" #include "math/Vector3.hpp" #include "primitives/Molecule.hpp" +#include "UseTheForce/fCutoffPolicy.h" +#include "UseTheForce/DarkSide/fElectrostaticSummationMethod.h" +#include "UseTheForce/DarkSide/fElectrostaticScreeningMethod.h" +#include "UseTheForce/DarkSide/fSwitchingFunctionType.h" #include "UseTheForce/doForces_interface.h" -#include "UseTheForce/notifyCutoffs_interface.h" +#include "UseTheForce/DarkSide/electrostatic_interface.h" +#include "UseTheForce/DarkSide/switcheroo_interface.h" #include "utils/MemoryUtils.hpp" #include "utils/simError.h" #include "selection/SelectionManager.hpp" +#include "io/ForceFieldOptions.hpp" +#include "UseTheForce/ForceField.hpp" #ifdef IS_MPI #include "UseTheForce/mpiComponentPlan.h" @@ -64,47 +72,54 @@ namespace oopse { #endif namespace oopse { + std::set getRigidSet(int index, std::map >& container) { + std::map >::iterator i = container.find(index); + std::set result; + if (i != container.end()) { + result = i->second; + } - 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), + return result; + } + + SimInfo::SimInfo(ForceField* ff, Globals* simParams) : + forceField_(ff), simParams_(simParams), + ndf_(0), fdf_local(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 + 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; + std::vector components = simParams->getComponents(); + + for (std::vector::iterator i = components.begin(); i !=components.end(); ++i) { + molStamp = (*i)->getMoleculeStamp(); + nMolWithSameStamp = (*i)->getNMol(); addMoleculeStamp(molStamp, nMolWithSameStamp); //calculate atoms in molecules nGlobalAtoms_ += molStamp->getNAtoms() *nMolWithSameStamp; - //calculate atoms in cutoff groups int nAtomsInGroups = 0; int nCutoffGroupsInStamp = molStamp->getNCutoffGroups(); for (int j=0; j < nCutoffGroupsInStamp; j++) { - cgStamp = molStamp->getCutoffGroup(j); + cgStamp = molStamp->getCutoffGroupStamp(j); nAtomsInGroups += cgStamp->getNMembers(); } nGroups += nCutoffGroupsInStamp * nMolWithSameStamp; + nCutoffAtoms += nAtomsInGroups * nMolWithSameStamp; //calculate atoms in rigid bodies @@ -112,7 +127,7 @@ namespace oopse { int nRigidBodiesInStamp = molStamp->getNRigidBodies(); for (int j=0; j < nRigidBodiesInStamp; j++) { - rbStamp = molStamp->getRigidBody(j); + rbStamp = molStamp->getRigidBodyStamp(j); nAtomsInRigidBodies += rbStamp->getNMembers(); } @@ -121,18 +136,21 @@ namespace oopse { } - //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 + //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(); #ifdef IS_MPI @@ -148,7 +166,6 @@ namespace oopse { } molecules_.clear(); - delete stamps_; delete sman_; delete simParams_; delete forceField_; @@ -255,8 +272,8 @@ namespace oopse { } } - }//end for (integrableObject) - }// end for (mol) + } + } // n_constraints is local, so subtract them on each processor ndf_local -= nConstraints_; @@ -273,6 +290,15 @@ namespace oopse { } + int SimInfo::getFdf() { +#ifdef IS_MPI + MPI_Allreduce(&fdf_local,&fdf_,1,MPI_INT,MPI_SUM, MPI_COMM_WORLD); +#else + fdf_ = fdf_local; +#endif + return fdf_; + } + void SimInfo::calcNdfRaw() { int ndfRaw_local; @@ -335,7 +361,36 @@ namespace oopse { int b; int c; int d; + + std::map > atomGroups; + + Molecule::RigidBodyIterator rbIter; + RigidBody* rb; + Molecule::IntegrableObjectIterator ii; + StuntDouble* integrableObject; + + for (integrableObject = mol->beginIntegrableObject(ii); integrableObject != NULL; + integrableObject = mol->nextIntegrableObject(ii)) { + + if (integrableObject->isRigidBody()) { + rb = static_cast(integrableObject); + std::vector atoms = rb->getAtoms(); + std::set rigidAtoms; + for (int i = 0; i < atoms.size(); ++i) { + rigidAtoms.insert(atoms[i]->getGlobalIndex()); + } + for (int i = 0; i < atoms.size(); ++i) { + atomGroups.insert(std::map >::value_type(atoms[i]->getGlobalIndex(), rigidAtoms)); + } + } else { + std::set oneAtomSet; + oneAtomSet.insert(integrableObject->getGlobalIndex()); + atomGroups.insert(std::map >::value_type(integrableObject->getGlobalIndex(), oneAtomSet)); + } + } + + for (bond= mol->beginBond(bondIter); bond != NULL; bond = mol->nextBond(bondIter)) { a = bond->getAtomA()->getGlobalIndex(); b = bond->getAtomB()->getGlobalIndex(); @@ -346,10 +401,17 @@ namespace oopse { a = bend->getAtomA()->getGlobalIndex(); b = bend->getAtomB()->getGlobalIndex(); c = bend->getAtomC()->getGlobalIndex(); + std::set rigidSetA = getRigidSet(a, atomGroups); + std::set rigidSetB = getRigidSet(b, atomGroups); + std::set rigidSetC = getRigidSet(c, atomGroups); - exclude_.addPair(a, b); - exclude_.addPair(a, c); - exclude_.addPair(b, c); + exclude_.addPairs(rigidSetA, rigidSetB); + exclude_.addPairs(rigidSetA, rigidSetC); + exclude_.addPairs(rigidSetB, rigidSetC); + + //exclude_.addPair(a, b); + //exclude_.addPair(a, c); + //exclude_.addPair(b, c); } for (torsion= mol->beginTorsion(torsionIter); torsion != NULL; torsion = mol->nextTorsion(torsionIter)) { @@ -357,17 +419,36 @@ namespace oopse { b = torsion->getAtomB()->getGlobalIndex(); c = torsion->getAtomC()->getGlobalIndex(); d = torsion->getAtomD()->getGlobalIndex(); + std::set rigidSetA = getRigidSet(a, atomGroups); + std::set rigidSetB = getRigidSet(b, atomGroups); + std::set rigidSetC = getRigidSet(c, atomGroups); + std::set rigidSetD = getRigidSet(d, atomGroups); + exclude_.addPairs(rigidSetA, rigidSetB); + exclude_.addPairs(rigidSetA, rigidSetC); + exclude_.addPairs(rigidSetA, rigidSetD); + exclude_.addPairs(rigidSetB, rigidSetC); + exclude_.addPairs(rigidSetB, rigidSetD); + exclude_.addPairs(rigidSetC, rigidSetD); + + /* + exclude_.addPairs(rigidSetA.begin(), rigidSetA.end(), rigidSetB.begin(), rigidSetB.end()); + exclude_.addPairs(rigidSetA.begin(), rigidSetA.end(), rigidSetC.begin(), rigidSetC.end()); + exclude_.addPairs(rigidSetA.begin(), rigidSetA.end(), rigidSetD.begin(), rigidSetD.end()); + exclude_.addPairs(rigidSetB.begin(), rigidSetB.end(), rigidSetC.begin(), rigidSetC.end()); + exclude_.addPairs(rigidSetB.begin(), rigidSetB.end(), rigidSetD.begin(), rigidSetD.end()); + exclude_.addPairs(rigidSetC.begin(), rigidSetC.end(), rigidSetD.begin(), rigidSetD.end()); + + 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) { @@ -392,7 +473,35 @@ namespace oopse { int b; int c; int d; + + std::map > atomGroups; + + Molecule::RigidBodyIterator rbIter; + RigidBody* rb; + Molecule::IntegrableObjectIterator ii; + StuntDouble* integrableObject; + for (integrableObject = mol->beginIntegrableObject(ii); integrableObject != NULL; + integrableObject = mol->nextIntegrableObject(ii)) { + + if (integrableObject->isRigidBody()) { + rb = static_cast(integrableObject); + std::vector atoms = rb->getAtoms(); + std::set rigidAtoms; + for (int i = 0; i < atoms.size(); ++i) { + rigidAtoms.insert(atoms[i]->getGlobalIndex()); + } + for (int i = 0; i < atoms.size(); ++i) { + atomGroups.insert(std::map >::value_type(atoms[i]->getGlobalIndex(), rigidAtoms)); + } + } else { + std::set oneAtomSet; + oneAtomSet.insert(integrableObject->getGlobalIndex()); + atomGroups.insert(std::map >::value_type(integrableObject->getGlobalIndex(), oneAtomSet)); + } + } + + for (bond= mol->beginBond(bondIter); bond != NULL; bond = mol->nextBond(bondIter)) { a = bond->getAtomA()->getGlobalIndex(); b = bond->getAtomB()->getGlobalIndex(); @@ -404,9 +513,17 @@ namespace oopse { b = bend->getAtomB()->getGlobalIndex(); c = bend->getAtomC()->getGlobalIndex(); - exclude_.removePair(a, b); - exclude_.removePair(a, c); - exclude_.removePair(b, c); + std::set rigidSetA = getRigidSet(a, atomGroups); + std::set rigidSetB = getRigidSet(b, atomGroups); + std::set rigidSetC = getRigidSet(c, atomGroups); + + exclude_.removePairs(rigidSetA, rigidSetB); + exclude_.removePairs(rigidSetA, rigidSetC); + exclude_.removePairs(rigidSetB, rigidSetC); + + //exclude_.removePair(a, b); + //exclude_.removePair(a, c); + //exclude_.removePair(b, c); } for (torsion= mol->beginTorsion(torsionIter); torsion != NULL; torsion = mol->nextTorsion(torsionIter)) { @@ -415,16 +532,36 @@ namespace oopse { c = torsion->getAtomC()->getGlobalIndex(); d = torsion->getAtomD()->getGlobalIndex(); + std::set rigidSetA = getRigidSet(a, atomGroups); + std::set rigidSetB = getRigidSet(b, atomGroups); + std::set rigidSetC = getRigidSet(c, atomGroups); + std::set rigidSetD = getRigidSet(d, atomGroups); + + exclude_.removePairs(rigidSetA, rigidSetB); + exclude_.removePairs(rigidSetA, rigidSetC); + exclude_.removePairs(rigidSetA, rigidSetD); + exclude_.removePairs(rigidSetB, rigidSetC); + exclude_.removePairs(rigidSetB, rigidSetD); + exclude_.removePairs(rigidSetC, rigidSetD); + + /* + exclude_.removePairs(rigidSetA.begin(), rigidSetA.end(), rigidSetB.begin(), rigidSetB.end()); + exclude_.removePairs(rigidSetA.begin(), rigidSetA.end(), rigidSetC.begin(), rigidSetC.end()); + exclude_.removePairs(rigidSetA.begin(), rigidSetA.end(), rigidSetD.begin(), rigidSetD.end()); + exclude_.removePairs(rigidSetB.begin(), rigidSetB.end(), rigidSetC.begin(), rigidSetC.end()); + exclude_.removePairs(rigidSetB.begin(), rigidSetB.end(), rigidSetD.begin(), rigidSetD.end()); + exclude_.removePairs(rigidSetC.begin(), rigidSetC.end(), rigidSetD.begin(), rigidSetD.end()); + + 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) { @@ -462,7 +599,10 @@ namespace oopse { //setup fortran force field /** @deprecate */ int isError = 0; - initFortranFF( &fInfo_.SIM_uses_RF , &isError ); + + setupElectrostaticSummationMethod( isError ); + setupSwitchingFunction(); + if(isError){ sprintf( painCave.errMsg, "ForceField error: There was an error initializing the forceField in fortran.\n" ); @@ -506,6 +646,7 @@ namespace oopse { int useLennardJones = 0; int useElectrostatic = 0; int useEAM = 0; + int useSC = 0; int useCharge = 0; int useDirectional = 0; int useDipole = 0; @@ -517,14 +658,34 @@ namespace oopse { int useDirectionalAtom = 0; int useElectrostatics = 0; //usePBC and useRF are from simParams - int usePBC = simParams_->getPBC(); - int useRF = simParams_->getUseRF(); + int usePBC = simParams_->getUsePeriodicBoundaryConditions(); + int useRF; + int useSF; + std::string myMethod; + // set the useRF logical + useRF = 0; + useSF = 0; + + + if (simParams_->haveElectrostaticSummationMethod()) { + std::string myMethod = simParams_->getElectrostaticSummationMethod(); + toUpper(myMethod); + if (myMethod == "REACTION_FIELD") { + useRF=1; + } else { + if (myMethod == "SHIFTED_FORCE") { + useSF = 1; + } + } + } + //loop over all of the atom types for (i = atomTypes.begin(); i != atomTypes.end(); ++i) { useLennardJones |= (*i)->isLennardJones(); useElectrostatic |= (*i)->isElectrostatic(); useEAM |= (*i)->isEAM(); + useSC |= (*i)->isSC(); useCharge |= (*i)->isCharge(); useDirectional |= (*i)->isDirectional(); useDipole |= (*i)->isDipole(); @@ -575,6 +736,9 @@ namespace oopse { temp = useEAM; MPI_Allreduce(&temp, &useEAM, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD); + temp = useSC; + MPI_Allreduce(&temp, &useSC, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD); + temp = useShape; MPI_Allreduce(&temp, &useShape, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD); @@ -583,7 +747,10 @@ namespace oopse { temp = useRF; MPI_Allreduce(&temp, &useRF, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD); - + + temp = useSF; + MPI_Allreduce(&temp, &useSF, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD); + #endif fInfo_.SIM_uses_PBC = usePBC; @@ -596,12 +763,14 @@ namespace oopse { fInfo_.SIM_uses_StickyPower = useStickyPower; fInfo_.SIM_uses_GayBerne = useGayBerne; fInfo_.SIM_uses_EAM = useEAM; + fInfo_.SIM_uses_SC = useSC; fInfo_.SIM_uses_Shapes = useShape; fInfo_.SIM_uses_FLARB = useFLARB; fInfo_.SIM_uses_RF = useRF; + fInfo_.SIM_uses_SF = useSF; - if( fInfo_.SIM_uses_Dipoles && fInfo_.SIM_uses_RF) { - + if( myMethod == "REACTION_FIELD") { + if (simParams_->haveDielectric()) { fInfo_.dielect = simParams_->getDielectric(); } else { @@ -611,10 +780,7 @@ namespace oopse { "\tsetting a dielectric constant!\n"); painCave.isFatal = 1; simError(); - } - - } else { - fInfo_.dielect = 0.0; + } } } @@ -650,7 +816,11 @@ namespace oopse { totalMass = cg->getMass(); for(atom = cg->beginAtom(ai); atom != NULL; atom = cg->nextAtom(ai)) { - mfact.push_back(atom->getMass()/totalMass); + // Check for massless groups - set mfact to 1 if true + if (totalMass != 0) + mfact.push_back(atom->getMass()/totalMass); + else + mfact.push_back( 1.0 ); } } @@ -759,83 +929,236 @@ namespace oopse { #endif - double SimInfo::calcMaxCutoffRadius() { + void SimInfo::setupCutoff() { + + ForceFieldOptions& forceFieldOptions_ = forceField_->getForceFieldOptions(); + // Check the cutoff policy + int cp = TRADITIONAL_CUTOFF_POLICY; // Set to traditional by default - 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)); + std::string myPolicy; + if (forceFieldOptions_.haveCutoffPolicy()){ + myPolicy = forceFieldOptions_.getCutoffPolicy(); + }else if (simParams_->haveCutoffPolicy()) { + myPolicy = simParams_->getCutoffPolicy(); } - double maxCutoffRadius = *(std::max_element(cutoffRadius.begin(), cutoffRadius.end())); -#ifdef IS_MPI - //pick the max cutoff radius among the processors -#endif + if (!myPolicy.empty()){ + toUpper(myPolicy); + 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(); + } + } + } + } + notifyFortranCutoffPolicy(&cp); - return maxCutoffRadius; - } - - void SimInfo::getCutoff(double& rcut, double& rsw) { - - if (fInfo_.SIM_uses_Charges | fInfo_.SIM_uses_Dipoles | fInfo_.SIM_uses_RF) { + // Check the Skin Thickness for neighborlists + double skin; + if (simParams_->haveSkinThickness()) { + skin = simParams_->getSkinThickness(); + notifyFortranSkinThickness(&skin); + } - if (!simParams_->haveRcut()){ - sprintf(painCave.errMsg, + // Check if the cutoff was set explicitly: + if (simParams_->haveCutoffRadius()) { + rcut_ = simParams_->getCutoffRadius(); + if (simParams_->haveSwitchingRadius()) { + rsw_ = simParams_->getSwitchingRadius(); + } else { + if (fInfo_.SIM_uses_Charges | + fInfo_.SIM_uses_Dipoles | + fInfo_.SIM_uses_RF) { + + rsw_ = 0.85 * rcut_; + sprintf(painCave.errMsg, + "SimCreator Warning: No value was set for the switchingRadius.\n" + "\tOOPSE will use a default value of 85 percent of the cutoffRadius.\n" + "\tswitchingRadius = %f. for this simulation\n", rsw_); + painCave.isFatal = 0; + simError(); + } else { + rsw_ = rcut_; + sprintf(painCave.errMsg, + "SimCreator Warning: No value was set for the switchingRadius.\n" + "\tOOPSE will use the same value as the cutoffRadius.\n" + "\tswitchingRadius = %f. for this simulation\n", rsw_); + painCave.isFatal = 0; + simError(); + } + } + + notifyFortranCutoffs(&rcut_, &rsw_); + + } else { + + // For electrostatic atoms, we'll assume a large safe value: + if (fInfo_.SIM_uses_Charges | fInfo_.SIM_uses_Dipoles | fInfo_.SIM_uses_RF) { + 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; + painCave.isFatal = 0; simError(); - rcut = 15.0; - } else{ - rcut = simParams_->getRcut(); - } + rcut_ = 15.0; + + if (simParams_->haveElectrostaticSummationMethod()) { + std::string myMethod = simParams_->getElectrostaticSummationMethod(); + toUpper(myMethod); + if (myMethod == "SHIFTED_POTENTIAL" || myMethod == "SHIFTED_FORCE") { + if (simParams_->haveSwitchingRadius()){ + sprintf(painCave.errMsg, + "SimInfo Warning: A value was set for the switchingRadius\n" + "\teven though the electrostaticSummationMethod was\n" + "\tset to %s\n", myMethod.c_str()); + painCave.isFatal = 1; + simError(); + } + } + } + + if (simParams_->haveSwitchingRadius()){ + rsw_ = simParams_->getSwitchingRadius(); + } else { + sprintf(painCave.errMsg, + "SimCreator Warning: No value was set for switchingRadius.\n" + "\tOOPSE will use a default value of\n" + "\t0.85 * cutoffRadius for the switchingRadius\n"); + painCave.isFatal = 0; + simError(); + rsw_ = 0.85 * rcut_; + } + notifyFortranCutoffs(&rcut_, &rsw_); + } else { + // We didn't set rcut explicitly, and we don't have electrostatic atoms, so + // We'll punt and let fortran figure out the cutoffs later. + + notifyFortranYouAreOnYourOwn(); - 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(); - } + void SimInfo::setupElectrostaticSummationMethod( int isError ) { + + int errorOut; + int esm = NONE; + int sm = UNDAMPED; + double alphaVal; + double dielectric; - if (simParams_->haveRsw()) { - rsw = simParams_->getRsw(); + errorOut = isError; + alphaVal = simParams_->getDampingAlpha(); + dielectric = simParams_->getDielectric(); + + if (simParams_->haveElectrostaticSummationMethod()) { + std::string myMethod = simParams_->getElectrostaticSummationMethod(); + toUpper(myMethod); + if (myMethod == "NONE") { + esm = NONE; } else { - rsw = rcut; + if (myMethod == "SWITCHING_FUNCTION") { + esm = SWITCHING_FUNCTION; + } else { + if (myMethod == "SHIFTED_POTENTIAL") { + esm = SHIFTED_POTENTIAL; + } else { + if (myMethod == "SHIFTED_FORCE") { + esm = SHIFTED_FORCE; + } else { + if (myMethod == "REACTION_FIELD") { + esm = REACTION_FIELD; + } else { + // throw error + sprintf( painCave.errMsg, + "SimInfo error: Unknown electrostaticSummationMethod.\n" + "\t(Input file specified %s .)\n" + "\telectrostaticSummationMethod must be one of: \"none\",\n" + "\t\"shifted_potential\", \"shifted_force\", or \n" + "\t\"reaction_field\".\n", myMethod.c_str() ); + painCave.isFatal = 1; + simError(); + } + } + } + } } + } + if (simParams_->haveElectrostaticScreeningMethod()) { + std::string myScreen = simParams_->getElectrostaticScreeningMethod(); + toUpper(myScreen); + if (myScreen == "UNDAMPED") { + sm = UNDAMPED; + } else { + if (myScreen == "DAMPED") { + sm = DAMPED; + if (!simParams_->haveDampingAlpha()) { + //throw error + sprintf( painCave.errMsg, + "SimInfo warning: dampingAlpha was not specified in the input file.\n" + "\tA default value of %f (1/ang) will be used.\n", alphaVal); + painCave.isFatal = 0; + simError(); + } + } else { + // throw error + sprintf( painCave.errMsg, + "SimInfo error: Unknown electrostaticScreeningMethod.\n" + "\t(Input file specified %s .)\n" + "\telectrostaticScreeningMethod must be one of: \"undamped\"\n" + "or \"damped\".\n", myScreen.c_str() ); + painCave.isFatal = 1; + simError(); + } + } } + + // let's pass some summation method variables to fortran + setElectrostaticSummationMethod( &esm ); + setFortranElectrostaticMethod( &esm ); + setScreeningMethod( &sm ); + setDampingAlpha( &alphaVal ); + setReactionFieldDielectric( &dielectric ); + initFortranFF( &errorOut ); } - void SimInfo::setupCutoff() { - getCutoff(rcut_, rsw_); - double rnblist = rcut_ + 1; // skin of neighbor list + void SimInfo::setupSwitchingFunction() { + int ft = CUBIC; + + if (simParams_->haveSwitchingFunctionType()) { + std::string funcType = simParams_->getSwitchingFunctionType(); + toUpper(funcType); + if (funcType == "CUBIC") { + ft = CUBIC; + } else { + if (funcType == "FIFTH_ORDER_POLYNOMIAL") { + ft = FIFTH_ORDER_POLY; + } else { + // throw error + sprintf( painCave.errMsg, + "SimInfo error: Unknown switchingFunctionType. (Input file specified %s .)\n\tswitchingFunctionType must be one of: \"cubic\" or \"fifth_order_polynomial\".", funcType.c_str() ); + painCave.isFatal = 1; + simError(); + } + } + } - //Pass these cutoff radius etc. to fortran. This function should be called once and only once - notifyFortranCutoffs(&rcut_, &rsw_, &rnblist); + // send switching function notification to switcheroo + setFunctionType(&ft); + } void SimInfo::addProperty(GenericData* genData) { @@ -945,6 +1268,148 @@ namespace oopse { 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