--- trunk/src/brains/SimInfo.cpp 2006/09/22 22:19:59 1050 +++ trunk/src/brains/SimInfo.cpp 2008/04/25 15:14:47 1241 @@ -59,6 +59,7 @@ #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" @@ -67,6 +68,7 @@ #include "io/ForceFieldOptions.hpp" #include "UseTheForce/ForceField.hpp" + #ifdef IS_MPI #include "UseTheForce/mpiComponentPlan.h" #include "UseTheForce/DarkSide/simParallel_interface.h" @@ -90,7 +92,8 @@ namespace oopse { nGlobalIntegrableObjects_(0), nGlobalRigidBodies_(0), nAtoms_(0), nBonds_(0), nBends_(0), nTorsions_(0), nRigidBodies_(0), nIntegrableObjects_(0), nCutoffGroups_(0), nConstraints_(0), - sman_(NULL), fortranInitialized_(false), calcBoxDipole_(false) { + sman_(NULL), fortranInitialized_(false), calcBoxDipole_(false), + useAtomicVirial_(true) { MoleculeStamp* molStamp; int nMolWithSameStamp; @@ -153,11 +156,7 @@ namespace oopse { + nGlobalRigidBodies_; nGlobalMols_ = molStampIds_.size(); - -#ifdef IS_MPI molToProcMap_.resize(nGlobalMols_); -#endif - } SimInfo::~SimInfo() { @@ -664,18 +663,20 @@ namespace oopse { 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; + useRF = 1; } else if (myMethod == "SHIFTED_FORCE"){ useSF = 1; } else if (myMethod == "SHIFTED_POTENTIAL"){ @@ -687,6 +688,8 @@ namespace oopse { 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(); @@ -763,6 +766,9 @@ namespace oopse { 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 @@ -783,6 +789,7 @@ namespace oopse { fInfo_.SIM_uses_SF = useSF; fInfo_.SIM_uses_SP = useSP; fInfo_.SIM_uses_BoxDipole = useBoxDipole; + fInfo_.SIM_uses_AtomicVirial = useAtomicVirial_; } void SimInfo::setupFortranSim() { @@ -849,30 +856,38 @@ namespace oopse { 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); - + setFortranSim( &fInfo_, &nGlobalAtoms_, &nAtoms_, &identArray[0], + &nExclude, excludeList , &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(); } - -#ifdef IS_MPI + + sprintf( checkPointMsg, "succesfully sent the simulation information to fortran.\n"); - MPIcheckPoint(); -#endif // is_mpi + + errorCheckPoint(); + + // Setup number of neighbors in neighbor list if present + if (simParams_->haveNeighborListNeighbors()) { + int nlistNeighbors = simParams_->getNeighborListNeighbors(); + setNeighbors(&nlistNeighbors); + } + + } -#ifdef IS_MPI void SimInfo::setupFortranParallel() { - +#ifdef IS_MPI //SimInfo is responsible for creating localToGlobalAtomIndex and localToGlobalGroupIndex std::vector localToGlobalAtomIndex(getNAtoms(), 0); std::vector localToGlobalCutoffGroupIndex; @@ -922,13 +937,11 @@ namespace oopse { } sprintf(checkPointMsg, " mpiRefresh successful.\n"); - MPIcheckPoint(); + errorCheckPoint(); - +#endif } -#endif - void SimInfo::setupCutoff() { ForceFieldOptions& forceFieldOptions_ = forceField_->getForceFieldOptions(); @@ -936,6 +949,10 @@ namespace oopse { // Check the cutoff policy 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(); @@ -999,8 +1016,18 @@ namespace oopse { 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 { @@ -1017,7 +1044,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" @@ -1040,7 +1075,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. @@ -1125,7 +1162,10 @@ namespace oopse { "\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, @@ -1442,6 +1482,43 @@ namespace oopse { 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_);