--- trunk/src/brains/SimInfo.cpp 2005/11/21 22:59:21 764 +++ trunk/src/brains/SimInfo.cpp 2007/04/20 18:15:48 1129 @@ -53,17 +53,22 @@ #include "brains/SimInfo.hpp" #include "math/Vector3.hpp" #include "primitives/Molecule.hpp" +#include "primitives/StuntDouble.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/DarkSide/neighborLists_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" #include "UseTheForce/DarkSide/simParallel_interface.h" @@ -80,18 +85,16 @@ namespace oopse { return result; } - 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), + 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) { + sman_(NULL), fortranInitialized_(false), calcBoxDipole_(false), + useAtomicVirial_(true) { - - std::vector >::iterator i; MoleculeStamp* molStamp; int nMolWithSameStamp; int nCutoffAtoms = 0; // number of atoms belong to cutoff groups @@ -99,23 +102,23 @@ namespace oopse { 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(); } @@ -128,7 +131,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(); } @@ -167,7 +170,6 @@ namespace oopse { } molecules_.clear(); - delete stamps_; delete sman_; delete simParams_; delete forceField_; @@ -274,8 +276,8 @@ namespace oopse { } } - }//end for (integrableObject) - }// end for (mol) + } + } // n_constraints is local, so subtract them on each processor ndf_local -= nConstraints_; @@ -292,6 +294,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; @@ -593,8 +604,11 @@ namespace oopse { /** @deprecate */ int isError = 0; + setupCutoff(); + setupElectrostaticSummationMethod( isError ); setupSwitchingFunction(); + setupAccumulateBoxDipole(); if(isError){ sprintf( painCave.errMsg, @@ -602,9 +616,6 @@ namespace oopse { painCave.isFatal = 1; simError(); } - - - setupCutoff(); calcNdf(); calcNdfRaw(); @@ -654,25 +665,35 @@ namespace oopse { int usePBC = simParams_->getUsePeriodicBoundaryConditions(); int useRF; int useSF; + int useSP; + int useBoxDipole; + std::string myMethod; // set the useRF logical useRF = 0; useSF = 0; + useSP = 0; if (simParams_->haveElectrostaticSummationMethod()) { std::string myMethod = simParams_->getElectrostaticSummationMethod(); toUpper(myMethod); - if (myMethod == "REACTION_FIELD") { - useRF=1; - } else { - if (myMethod == "SHIFTED_FORCE") { - useSF = 1; - } + if (myMethod == "REACTION_FIELD"){ + useRF = 1; + } else if (myMethod == "SHIFTED_FORCE"){ + useSF = 1; + } else if (myMethod == "SHIFTED_POTENTIAL"){ + useSP = 1; } } + + if (simParams_->haveAccumulateBoxDipole()) + if (simParams_->getAccumulateBoxDipole()) + useBoxDipole = 1; + useAtomicVirial_ = simParams_->getUseAtomicVirial(); + //loop over all of the atom types for (i = atomTypes.begin(); i != atomTypes.end(); ++i) { useLennardJones |= (*i)->isLennardJones(); @@ -742,8 +763,17 @@ namespace oopse { 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); + MPI_Allreduce(&temp, &useSF, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD); + temp = useSP; + MPI_Allreduce(&temp, &useSP, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD); + + temp = useBoxDipole; + MPI_Allreduce(&temp, &useBoxDipole, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD); + + temp = useAtomicVirial_; + MPI_Allreduce(&temp, &useAtomicVirial_, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD); + #endif fInfo_.SIM_uses_PBC = usePBC; @@ -761,21 +791,9 @@ namespace oopse { fInfo_.SIM_uses_FLARB = useFLARB; fInfo_.SIM_uses_RF = useRF; fInfo_.SIM_uses_SF = useSF; - - if( myMethod == "REACTION_FIELD") { - - 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(); - } - } - + fInfo_.SIM_uses_SP = useSP; + fInfo_.SIM_uses_BoxDipole = useBoxDipole; + fInfo_.SIM_uses_AtomicVirial = useAtomicVirial_; } void SimInfo::setupFortranSim() { @@ -792,14 +810,14 @@ namespace oopse { } //calculate mass ratio of cutoff group - std::vector mfact; + std::vector mfact; SimInfo::MoleculeIterator mi; Molecule* mol; Molecule::CutoffGroupIterator ci; CutoffGroup* cg; Molecule::AtomIterator ai; Atom* atom; - double totalMass; + RealType totalMass; //to avoid memory reallocation, reserve enough space for mfact mfact.reserve(getNCutoffGroups()); @@ -860,6 +878,14 @@ namespace oopse { "succesfully sent the simulation information to fortran.\n"); MPIcheckPoint(); #endif // is_mpi + + // Setup number of neighbors in neighbor list if present + if (simParams_->haveNeighborListNeighbors()) { + int nlistNeighbors = simParams_->getNeighborListNeighbors(); + setNeighbors(&nlistNeighbors); + } + + } @@ -924,10 +950,23 @@ namespace oopse { void SimInfo::setupCutoff() { + ForceFieldOptions& forceFieldOptions_ = forceField_->getForceFieldOptions(); + // Check the cutoff policy - int cp = TRADITIONAL_CUTOFF_POLICY; - if (simParams_->haveCutoffPolicy()) { - std::string myPolicy = simParams_->getCutoffPolicy(); + int cp = TRADITIONAL_CUTOFF_POLICY; // Set to traditional by default + + // Set LJ shifting bools to false + ljsp_ = false; + ljsf_ = false; + + std::string myPolicy; + if (forceFieldOptions_.haveCutoffPolicy()){ + myPolicy = forceFieldOptions_.getCutoffPolicy(); + }else if (simParams_->haveCutoffPolicy()) { + myPolicy = simParams_->getCutoffPolicy(); + } + + if (!myPolicy.empty()){ toUpper(myPolicy); if (myPolicy == "MIX") { cp = MIX_CUTOFF_POLICY; @@ -950,7 +989,7 @@ namespace oopse { notifyFortranCutoffPolicy(&cp); // Check the Skin Thickness for neighborlists - double skin; + RealType skin; if (simParams_->haveSkinThickness()) { skin = simParams_->getSkinThickness(); notifyFortranSkinThickness(&skin); @@ -962,9 +1001,39 @@ namespace oopse { if (simParams_->haveSwitchingRadius()) { rsw_ = simParams_->getSwitchingRadius(); } else { - rsw_ = rcut_; + 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_); + + if (simParams_->haveElectrostaticSummationMethod()) { + std::string myMethod = simParams_->getElectrostaticSummationMethod(); + toUpper(myMethod); + + if (myMethod == "SHIFTED_POTENTIAL") { + ljsp_ = true; + } else if (myMethod == "SHIFTED_FORCE") { + ljsf_ = true; + } + } + notifyFortranCutoffs(&rcut_, &rsw_, &ljsp_, &ljsf_); } else { @@ -981,7 +1050,15 @@ namespace oopse { if (simParams_->haveElectrostaticSummationMethod()) { std::string myMethod = simParams_->getElectrostaticSummationMethod(); toUpper(myMethod); - if (myMethod == "SHIFTED_POTENTIAL" || myMethod == "SHIFTED_FORCE") { + + // For the time being, we're tethering the LJ shifted behavior to the + // electrostaticSummationMethod keyword options + if (myMethod == "SHIFTED_POTENTIAL") { + ljsp_ = true; + } else if (myMethod == "SHIFTED_FORCE") { + ljsf_ = true; + } + if (myMethod == "SHIFTED_POTENTIAL" || myMethod == "SHIFTED_FORCE") { if (simParams_->haveSwitchingRadius()){ sprintf(painCave.errMsg, "SimInfo Warning: A value was set for the switchingRadius\n" @@ -1004,7 +1081,9 @@ namespace oopse { simError(); rsw_ = 0.85 * rcut_; } - notifyFortranCutoffs(&rcut_, &rsw_); + + notifyFortranCutoffs(&rcut_, &rsw_, &ljsp_, &ljsf_); + } 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. @@ -1020,12 +1099,10 @@ namespace oopse { int errorOut; int esm = NONE; int sm = UNDAMPED; - double alphaVal; - double dielectric; - + RealType alphaVal; + RealType dielectric; + errorOut = isError; - alphaVal = simParams_->getDampingAlpha(); - dielectric = simParams_->getDielectric(); if (simParams_->haveElectrostaticSummationMethod()) { std::string myMethod = simParams_->getElectrostaticSummationMethod(); @@ -1042,8 +1119,17 @@ namespace oopse { if (myMethod == "SHIFTED_FORCE") { esm = SHIFTED_FORCE; } else { - if (myMethod == "REACTION_FIELD") { + if (myMethod == "REACTION_FIELD") { esm = REACTION_FIELD; + dielectric = simParams_->getDielectric(); + if (!simParams_->haveDielectric()) { + // throw warning + sprintf( painCave.errMsg, + "SimInfo warning: dielectric was not specified in the input file\n\tfor the reaction field correction method.\n" + "\tA default value of %f will be used for the dielectric.\n", dielectric); + painCave.isFatal = 0; + simError(); + } } else { // throw error sprintf( painCave.errMsg, @@ -1070,13 +1156,22 @@ namespace oopse { if (myScreen == "DAMPED") { sm = DAMPED; if (!simParams_->haveDampingAlpha()) { - //throw error + // first set a cutoff dependent alpha value + // we assume alpha depends linearly with rcut from 0 to 20.5 ang + alphaVal = 0.5125 - rcut_* 0.025; + // for values rcut > 20.5, alpha is zero + if (alphaVal < 0) alphaVal = 0; + + // throw warning 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); + "\tA default value of %f (1/ang) will be used for the cutoff of\n\t%f (ang).\n", alphaVal, rcut_); painCave.isFatal = 0; simError(); + } else { + alphaVal = simParams_->getDampingAlpha(); } + } else { // throw error sprintf( painCave.errMsg, @@ -1092,7 +1187,7 @@ namespace oopse { // let's pass some summation method variables to fortran setElectrostaticSummationMethod( &esm ); - notifyFortranElectrostaticMethod( &esm ); + setFortranElectrostaticMethod( &esm ); setScreeningMethod( &sm ); setDampingAlpha( &alphaVal ); setReactionFieldDielectric( &dielectric ); @@ -1125,6 +1220,17 @@ namespace oopse { } + void SimInfo::setupAccumulateBoxDipole() { + + // we only call setAccumulateBoxDipole if the accumulateBoxDipole parameter is true + if ( simParams_->haveAccumulateBoxDipole() ) + if ( simParams_->getAccumulateBoxDipole() ) { + setAccumulateBoxDipole(); + calcBoxDipole_ = true; + } + + } + void SimInfo::addProperty(GenericData* genData) { properties_.addProperty(genData); } @@ -1181,20 +1287,20 @@ namespace oopse { Molecule* mol; Vector3d comVel(0.0); - double totalMass = 0.0; + RealType totalMass = 0.0; for (mol = beginMolecule(i); mol != NULL; mol = nextMolecule(i)) { - double mass = mol->getMass(); + RealType mass = mol->getMass(); totalMass += mass; comVel += mass * mol->getComVel(); } #ifdef IS_MPI - double tmpMass = totalMass; + RealType tmpMass = totalMass; Vector3d tmpComVel(comVel); - MPI_Allreduce(&tmpMass,&totalMass,1,MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD); - MPI_Allreduce(tmpComVel.getArrayPointer(), comVel.getArrayPointer(),3,MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD); + MPI_Allreduce(&tmpMass,&totalMass,1,MPI_REALTYPE,MPI_SUM, MPI_COMM_WORLD); + MPI_Allreduce(tmpComVel.getArrayPointer(), comVel.getArrayPointer(),3,MPI_REALTYPE,MPI_SUM, MPI_COMM_WORLD); #endif comVel /= totalMass; @@ -1207,19 +1313,19 @@ namespace oopse { Molecule* mol; Vector3d com(0.0); - double totalMass = 0.0; + RealType totalMass = 0.0; for (mol = beginMolecule(i); mol != NULL; mol = nextMolecule(i)) { - double mass = mol->getMass(); + RealType mass = mol->getMass(); totalMass += mass; com += mass * mol->getCom(); } #ifdef IS_MPI - double tmpMass = totalMass; + RealType tmpMass = totalMass; Vector3d tmpCom(com); - 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(&tmpMass,&totalMass,1,MPI_REALTYPE,MPI_SUM, MPI_COMM_WORLD); + MPI_Allreduce(tmpCom.getArrayPointer(), com.getArrayPointer(),3,MPI_REALTYPE,MPI_SUM, MPI_COMM_WORLD); #endif com /= totalMass; @@ -1243,23 +1349,23 @@ namespace oopse { Molecule* mol; - double totalMass = 0.0; + RealType totalMass = 0.0; for (mol = beginMolecule(i); mol != NULL; mol = nextMolecule(i)) { - double mass = mol->getMass(); + RealType mass = mol->getMass(); totalMass += mass; com += mass * mol->getCom(); comVel += mass * mol->getComVel(); } #ifdef IS_MPI - double tmpMass = totalMass; + RealType 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); + MPI_Allreduce(&tmpMass,&totalMass,1,MPI_REALTYPE,MPI_SUM, MPI_COMM_WORLD); + MPI_Allreduce(tmpCom.getArrayPointer(), com.getArrayPointer(),3,MPI_REALTYPE,MPI_SUM, MPI_COMM_WORLD); + MPI_Allreduce(tmpComVel.getArrayPointer(), comVel.getArrayPointer(),3,MPI_REALTYPE,MPI_SUM, MPI_COMM_WORLD); #endif com /= totalMass; @@ -1278,12 +1384,12 @@ namespace oopse { 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; + RealType xx = 0.0; + RealType yy = 0.0; + RealType zz = 0.0; + RealType xy = 0.0; + RealType xz = 0.0; + RealType yz = 0.0; Vector3d com(0.0); Vector3d comVel(0.0); @@ -1295,7 +1401,7 @@ namespace oopse { Vector3d thisq(0.0); Vector3d thisv(0.0); - double thisMass = 0.0; + RealType thisMass = 0.0; @@ -1333,8 +1439,8 @@ namespace oopse { #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); + MPI_Allreduce(tmpI.getArrayPointer(), inertiaTensor.getArrayPointer(),9,MPI_REALTYPE,MPI_SUM, MPI_COMM_WORLD); + MPI_Allreduce(tmpAngMom.getArrayPointer(), angularMomentum.getArrayPointer(),3,MPI_REALTYPE,MPI_SUM, MPI_COMM_WORLD); #endif return; @@ -1355,7 +1461,7 @@ namespace oopse { Vector3d thisr(0.0); Vector3d thisp(0.0); - double thisMass; + RealType thisMass; for (mol = beginMolecule(i); mol != NULL; mol = nextMolecule(i)) { thisMass = mol->getMass(); @@ -1368,12 +1474,67 @@ namespace oopse { #ifdef IS_MPI Vector3d tmpAngMom; - MPI_Allreduce(tmpAngMom.getArrayPointer(), angularMomentum.getArrayPointer(),3,MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD); + MPI_Allreduce(tmpAngMom.getArrayPointer(), angularMomentum.getArrayPointer(),3,MPI_REALTYPE,MPI_SUM, MPI_COMM_WORLD); #endif return angularMomentum; } - + StuntDouble* SimInfo::getIOIndexToIntegrableObject(int index) { + return IOIndexToIntegrableObject.at(index); + } + + void SimInfo::setIOIndexToIntegrableObject(const std::vector& v) { + IOIndexToIntegrableObject= v; + } + + /* Returns the Volume of the simulation based on a ellipsoid with semi-axes + based on the radius of gyration V=4/3*Pi*R_1*R_2*R_3 + where R_i are related to the principle inertia moments R_i = sqrt(C*I_i/N), this reduces to + V = 4/3*Pi*(C/N)^3/2*sqrt(det(I)). See S.E. Baltazar et. al. Comp. Mat. Sci. 37 (2006) 526-536. + */ + void SimInfo::getGyrationalVolume(RealType &volume){ + Mat3x3d intTensor; + RealType det; + Vector3d dummyAngMom; + RealType sysconstants; + RealType geomCnst; + + geomCnst = 3.0/2.0; + /* Get the inertial tensor and angular momentum for free*/ + getInertiaTensor(intTensor,dummyAngMom); + + det = intTensor.determinant(); + sysconstants = geomCnst/(RealType)nGlobalIntegrableObjects_; + volume = 4.0/3.0*NumericConstant::PI*pow(sysconstants,3.0/2.0)*sqrt(det); + return; + } + + void SimInfo::getGyrationalVolume(RealType &volume, RealType &detI){ + Mat3x3d intTensor; + Vector3d dummyAngMom; + RealType sysconstants; + RealType geomCnst; + + geomCnst = 3.0/2.0; + /* Get the inertial tensor and angular momentum for free*/ + getInertiaTensor(intTensor,dummyAngMom); + + detI = intTensor.determinant(); + sysconstants = geomCnst/(RealType)nGlobalIntegrableObjects_; + volume = 4.0/3.0*NumericConstant::PI*pow(sysconstants,3.0/2.0)*sqrt(detI); + return; + } +/* + void SimInfo::setStuntDoubleFromGlobalIndex(std::vector v) { + assert( v.size() == nAtoms_ + nRigidBodies_); + sdByGlobalIndex_ = v; + } + + StuntDouble* SimInfo::getStuntDoubleFromGlobalIndex(int index) { + //assert(index < nAtoms_ + nRigidBodies_); + return sdByGlobalIndex_.at(index); + } +*/ }//end namespace oopse