--- trunk/src/brains/SimInfo.cpp 2006/12/05 00:17:24 1095 +++ trunk/src/brains/SimInfo.cpp 2008/04/25 15:14:47 1241 @@ -92,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; @@ -155,11 +156,7 @@ namespace oopse { + nGlobalRigidBodies_; nGlobalMols_ = molStampIds_.size(); - -#ifdef IS_MPI molToProcMap_.resize(nGlobalMols_); -#endif - } SimInfo::~SimInfo() { @@ -666,6 +663,7 @@ namespace oopse { int useSF; int useSP; int useBoxDipole; + std::string myMethod; // set the useRF logical @@ -690,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(); @@ -767,6 +767,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 fInfo_.SIM_uses_PBC = usePBC; @@ -786,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() { @@ -852,37 +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()) { - setNeighbors(simParams_->getNeighborListNeighbors()); + 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; @@ -932,13 +937,11 @@ namespace oopse { } sprintf(checkPointMsg, " mpiRefresh successful.\n"); - MPIcheckPoint(); + errorCheckPoint(); - +#endif } -#endif - void SimInfo::setupCutoff() { ForceFieldOptions& forceFieldOptions_ = forceField_->getForceFieldOptions(); @@ -946,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(); @@ -1009,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 { @@ -1027,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" @@ -1050,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. @@ -1455,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_);