--- trunk/src/brains/SimInfo.cpp 2005/09/07 20:46:46 586 +++ trunk/src/brains/SimInfo.cpp 2005/10/19 19:24:40 691 @@ -53,7 +53,9 @@ #include "math/Vector3.hpp" #include "primitives/Molecule.hpp" #include "UseTheForce/fCutoffPolicy.h" +#include "UseTheForce/DarkSide/fElectrostaticSummationMethod.h" #include "UseTheForce/doForces_interface.h" +#include "UseTheForce/DarkSide/electrostatic_interface.h" #include "UseTheForce/notifyCutoffs_interface.h" #include "utils/MemoryUtils.hpp" #include "utils/simError.h" @@ -81,7 +83,7 @@ namespace oopse { 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; @@ -106,6 +108,7 @@ namespace oopse { } nGroups += nCutoffGroupsInStamp * nMolWithSameStamp; + nCutoffAtoms += nAtomsInGroups * nMolWithSameStamp; //calculate atoms in rigid bodies @@ -122,18 +125,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 @@ -463,8 +469,9 @@ namespace oopse { //setup fortran force field /** @deprecate */ int isError = 0; - initFortranFF( &fInfo_.SIM_uses_RF, &fInfo_.SIM_uses_UW, - &fInfo_.SIM_uses_DW, &isError ); + + setupElectrostaticSummationMethod( isError ); + if(isError){ sprintf( painCave.errMsg, "ForceField error: There was an error initializing the forceField in fortran.\n" ); @@ -519,11 +526,22 @@ namespace oopse { int useDirectionalAtom = 0; int useElectrostatics = 0; //usePBC and useRF are from simParams - int usePBC = simParams_->getPBC(); - int useRF = simParams_->getUseRF(); - int useUW = simParams_->getUseUndampedWolf(); - int useDW = simParams_->getUseDampedWolf(); + int usePBC = simParams_->getUsePeriodicBoundaryConditions(); + int useRF; + std::string myMethod; + // set the useRF logical + useRF = 0; + + + if (simParams_->haveElectrostaticSummationMethod()) { + std::string myMethod = simParams_->getElectrostaticSummationMethod(); + toUpper(myMethod); + if (myMethod == "REACTION_FIELD") { + useRF=1; + } + } + //loop over all of the atom types for (i = atomTypes.begin(); i != atomTypes.end(); ++i) { useLennardJones |= (*i)->isLennardJones(); @@ -588,12 +606,6 @@ namespace oopse { temp = useRF; MPI_Allreduce(&temp, &useRF, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD); - temp = useUW; - MPI_Allreduce(&temp, &useUW, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD); - - temp = useDW; - MPI_Allreduce(&temp, &useDW, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD); - #endif fInfo_.SIM_uses_PBC = usePBC; @@ -609,11 +621,9 @@ namespace oopse { fInfo_.SIM_uses_Shapes = useShape; fInfo_.SIM_uses_FLARB = useFLARB; fInfo_.SIM_uses_RF = useRF; - fInfo_.SIM_uses_UW = useUW; - fInfo_.SIM_uses_DW = useDW; - if( fInfo_.SIM_uses_Dipoles && fInfo_.SIM_uses_RF) { - + if( myMethod == "REACTION_FIELD") { + if (simParams_->haveDielectric()) { fInfo_.dielect = simParams_->getDielectric(); } else { @@ -623,12 +633,8 @@ namespace oopse { "\tsetting a dielectric constant!\n"); painCave.isFatal = 1; simError(); - } - - } else { - fInfo_.dielect = 0.0; + } } - } void SimInfo::setupFortranSim() { @@ -662,7 +668,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 ); } } @@ -798,7 +808,7 @@ namespace oopse { if (fInfo_.SIM_uses_Charges | fInfo_.SIM_uses_Dipoles | fInfo_.SIM_uses_RF) { - if (!simParams_->haveRcut()){ + if (!simParams_->haveCutoffRadius()){ sprintf(painCave.errMsg, "SimCreator Warning: No value was set for the cutoffRadius.\n" "\tOOPSE will use a default value of 15.0 angstroms" @@ -807,10 +817,10 @@ namespace oopse { simError(); rcut = 15.0; } else{ - rcut = simParams_->getRcut(); + rcut = simParams_->getCutoffRadius(); } - if (!simParams_->haveRsw()){ + if (!simParams_->haveSwitchingRadius()){ sprintf(painCave.errMsg, "SimCreator Warning: No value was set for switchingRadius.\n" "\tOOPSE will use a default value of\n" @@ -819,22 +829,22 @@ namespace oopse { simError(); rsw = 0.95 * rcut; } else{ - rsw = simParams_->getRsw(); + rsw = simParams_->getSwitchingRadius(); } } 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(); + if (simParams_->haveCutoffRadius()) { + rcut = simParams_->getCutoffRadius(); } else { //set cutoff radius to the maximum cutoff radius based on atom types in the whole system rcut = calcMaxCutoffRadius(); } - if (simParams_->haveRsw()) { - rsw = simParams_->getRsw(); + if (simParams_->haveSwitchingRadius()) { + rsw = simParams_->getSwitchingRadius(); } else { rsw = rcut; } @@ -851,6 +861,7 @@ namespace oopse { int cp = TRADITIONAL_CUTOFF_POLICY; if (simParams_->haveCutoffPolicy()) { std::string myPolicy = simParams_->getCutoffPolicy(); + toUpper(myPolicy); if (myPolicy == "MIX") { cp = MIX_CUTOFF_POLICY; } else { @@ -869,9 +880,67 @@ namespace oopse { } } } + + + if (simParams_->haveSkinThickness()) { + double skinThickness = simParams_->getSkinThickness(); + } + notifyFortranCutoffs(&rcut_, &rsw_, &rnblist, &cp); + // also send cutoff notification to electrostatics + setElectrostaticCutoffRadius(&rcut_, &rsw_); } + void SimInfo::setupElectrostaticSummationMethod( int isError ) { + + int errorOut; + int esm = NONE; + double alphaVal; + double dielectric; + + errorOut = isError; + alphaVal = simParams_->getDampingAlpha(); + dielectric = simParams_->getDielectric(); + + if (simParams_->haveElectrostaticSummationMethod()) { + std::string myMethod = simParams_->getElectrostaticSummationMethod(); + toUpper(myMethod); + if (myMethod == "NONE") { + esm = NONE; + } else { + if (myMethod == "UNDAMPED_WOLF") { + esm = UNDAMPED_WOLF; + } else { + if (myMethod == "DAMPED_WOLF") { + esm = DAMPED_WOLF; + if (!simParams_->haveDampingAlpha()) { + //throw error + sprintf( painCave.errMsg, + "SimInfo warning: dampingAlpha was not specified in the input file. A default value of %f (1/ang) will be used for the Damped Wolf Method.", alphaVal); + painCave.isFatal = 0; + simError(); + } + } else { + if (myMethod == "REACTION_FIELD") { + esm = REACTION_FIELD; + } else { + // throw error + sprintf( painCave.errMsg, + "SimInfo error: Unknown electrostaticSummationMethod. (Input file specified %s .)\n\telectrostaticSummationMethod must be one of: \"none\", \"undamped_wolf\", \"damped_wolf\", or \"reaction_field\".", myMethod.c_str() ); + painCave.isFatal = 1; + simError(); + } + } + } + } + } + // let's pass some summation method variables to fortran + setElectrostaticSummationMethod( &esm ); + setDampedWolfAlpha( &alphaVal ); + setReactionFieldDielectric( &dielectric ); + initFortranFF( &esm, &errorOut ); + } + void SimInfo::addProperty(GenericData* genData) { properties_.addProperty(genData); }