--- trunk/src/brains/SimInfo.cpp 2006/04/25 02:09:01 945 +++ trunk/src/brains/SimInfo.cpp 2006/09/21 18:25:17 1045 @@ -53,6 +53,7 @@ #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" @@ -89,7 +90,7 @@ 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) { + sman_(NULL), fortranInitialized_(false), calcBoxDipole_(false) { MoleculeStamp* molStamp; int nMolWithSameStamp; @@ -599,9 +600,12 @@ namespace oopse { //setup fortran force field /** @deprecate */ int isError = 0; + + setupCutoff(); setupElectrostaticSummationMethod( isError ); setupSwitchingFunction(); + setupAccumulateBoxDipole(); if(isError){ sprintf( painCave.errMsg, @@ -609,9 +613,6 @@ namespace oopse { painCave.isFatal = 1; simError(); } - - - setupCutoff(); calcNdf(); calcNdfRaw(); @@ -661,6 +662,8 @@ namespace oopse { int usePBC = simParams_->getUsePeriodicBoundaryConditions(); int useRF; int useSF; + int useSP; + int useBoxDipole; std::string myMethod; // set the useRF logical @@ -671,14 +674,18 @@ namespace oopse { if (simParams_->haveElectrostaticSummationMethod()) { std::string myMethod = simParams_->getElectrostaticSummationMethod(); toUpper(myMethod); - if (myMethod == "REACTION_FIELD") { + if (myMethod == "REACTION_FIELD"){ useRF=1; - } else { - if (myMethod == "SHIFTED_FORCE") { - useSF = 1; - } + } else if (myMethod == "SHIFTED_FORCE"){ + useSF = 1; + } else if (myMethod == "SHIFTED_POTENTIAL"){ + useSP = 1; } } + + if (simParams_->haveAccumulateBoxDipole()) + if (simParams_->getAccumulateBoxDipole()) + useBoxDipole = 1; //loop over all of the atom types for (i = atomTypes.begin(); i != atomTypes.end(); ++i) { @@ -749,8 +756,14 @@ 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); + #endif fInfo_.SIM_uses_PBC = usePBC; @@ -768,6 +781,8 @@ namespace oopse { fInfo_.SIM_uses_FLARB = useFLARB; fInfo_.SIM_uses_RF = useRF; fInfo_.SIM_uses_SF = useSF; + fInfo_.SIM_uses_SP = useSP; + fInfo_.SIM_uses_BoxDipole = useBoxDipole; if( myMethod == "REACTION_FIELD") { @@ -799,14 +814,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()); @@ -966,7 +981,7 @@ namespace oopse { notifyFortranCutoffPolicy(&cp); // Check the Skin Thickness for neighborlists - double skin; + RealType skin; if (simParams_->haveSkinThickness()) { skin = simParams_->getSkinThickness(); notifyFortranSkinThickness(&skin); @@ -1056,11 +1071,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()) { @@ -1106,10 +1120,16 @@ 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(); } @@ -1161,6 +1181,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); } @@ -1217,20 +1248,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; @@ -1243,19 +1274,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; @@ -1279,23 +1310,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; @@ -1314,12 +1345,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); @@ -1331,7 +1362,7 @@ namespace oopse { Vector3d thisq(0.0); Vector3d thisv(0.0); - double thisMass = 0.0; + RealType thisMass = 0.0; @@ -1369,8 +1400,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; @@ -1391,7 +1422,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(); @@ -1404,12 +1435,30 @@ 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; + } + +/* + 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