--- trunk/OOPSE/libmdtools/SimSetup.cpp 2003/09/22 21:23:25 780 +++ trunk/OOPSE/libmdtools/SimSetup.cpp 2004/01/19 21:17:39 965 @@ -1,10 +1,9 @@ #include -#include +#include #include -#include +#include #include #include - #include "SimSetup.hpp" #include "ReadWrite.hpp" #include "parse_me.h" @@ -22,14 +21,42 @@ #define NVT_ENS 1 #define NPTi_ENS 2 #define NPTf_ENS 3 +#define NPTxyz_ENS 4 + #define FF_DUFF 0 #define FF_LJ 1 #define FF_EAM 2 using namespace std; + +/** + * Check whether dividend is divisble by divisor or not + */ +bool isDivisible(double dividend, double divisor){ + double tolerance = 0.000001; + double quotient; + double diff; + int intQuotient; + + quotient = dividend / divisor; + + if (quotient < 0) + quotient = -quotient; + + intQuotient = int (quotient + tolerance); + + diff = fabs(fabs(dividend) - intQuotient * fabs(divisor)); + if (diff <= tolerance) + return true; + else + return false; +} + SimSetup::SimSetup(){ + + initSuspend = false; isInfoArray = 0; nInfo = 1; @@ -52,6 +79,7 @@ void SimSetup::setSimInfo(SimInfo* the_info, int theNi info = the_info; nInfo = theNinfo; isInfoArray = 1; + initSuspend = true; } @@ -90,7 +118,6 @@ void SimSetup::createSim(void){ #endif // is_mpi void SimSetup::createSim(void){ - int i, j, k, globalAtomIndex; // gather all of the information from the Bass file @@ -106,8 +133,11 @@ void SimSetup::createSim(void){ // initialize the system coordinates - if (!isInfoArray){ + if ( !initSuspend ){ initSystemCoords(); + + if( !(globals->getUseInitTime()) ) + info[0].currentTime = 0.0; } // make the output filenames @@ -129,7 +159,7 @@ void SimSetup::makeMolecules(void){ void SimSetup::makeMolecules(void){ - int k, l; + int k; int i, j, exI, exJ, tempEx, stampID, atomOffset, excludeOffset; molInit molInfo; DirectionalAtom* dAtom; @@ -551,7 +581,7 @@ void SimSetup::gatherInfo(void){ void SimSetup::gatherInfo(void){ - int i, j, k; + int i; ensembleCase = -1; ffCase = -1; @@ -602,10 +632,13 @@ void SimSetup::gatherInfo(void){ else if (!strcasecmp(ensemble, "NPTf")){ ensembleCase = NPTf_ENS; } + else if (!strcasecmp(ensemble, "NPTxyz")){ + ensembleCase = NPTxyz_ENS; + } else{ sprintf(painCave.errMsg, - "SimSetup Warning. Unrecognized Ensemble -> %s, " - "reverting to NVE for this simulation.\n", + "SimSetup Warning. Unrecognized Ensemble -> %s \n" + "\treverting to NVE for this simulation.\n", ensemble); painCave.isFatal = 0; simError(); @@ -637,8 +670,8 @@ void SimSetup::gatherInfo(void){ if (!the_components[i]->haveNMol()){ // we have a problem sprintf(painCave.errMsg, - "SimSetup Error. No global NMol or component NMol" - " given. Cannot calculate the number of atoms.\n"); + "SimSetup Error. No global NMol or component NMol given.\n" + "\tCannot calculate the number of atoms.\n"); painCave.isFatal = 1; simError(); } @@ -657,6 +690,47 @@ void SimSetup::gatherInfo(void){ painCave.isFatal = 1; simError(); } + + //check whether sample time, status time, thermal time and reset time are divisble by dt + if (!isDivisible(globals->getSampleTime(), globals->getDt())){ + sprintf(painCave.errMsg, + "Sample time is not divisible by dt.\n" + "\tThis will result in samples that are not uniformly\n" + "\tdistributed in time. If this is a problem, change\n" + "\tyour sampleTime variable.\n"); + painCave.isFatal = 0; + simError(); + } + + if (globals->haveStatusTime() && !isDivisible(globals->getSampleTime(), globals->getDt())){ + sprintf(painCave.errMsg, + "Status time is not divisible by dt.\n" + "\tThis will result in status reports that are not uniformly\n" + "\tdistributed in time. If this is a problem, change \n" + "\tyour statusTime variable.\n"); + painCave.isFatal = 0; + simError(); + } + + if (globals->haveThermalTime() && !isDivisible(globals->getThermalTime(), globals->getDt())){ + sprintf(painCave.errMsg, + "Thermal time is not divisible by dt.\n" + "\tThis will result in thermalizations that are not uniformly\n" + "\tdistributed in time. If this is a problem, change \n" + "\tyour thermalTime variable.\n"); + painCave.isFatal = 0; + simError(); + } + + if (globals->haveResetTime() && !isDivisible(globals->getResetTime(), globals->getDt())){ + sprintf(painCave.errMsg, + "Reset time is not divisible by dt.\n" + "\tThis will result in integrator resets that are not uniformly\n" + "\tdistributed in time. If this is a problem, change\n" + "\tyour resetTime variable.\n"); + painCave.isFatal = 0; + simError(); + } // set the status, sample, and thermal kick times @@ -687,58 +761,17 @@ void SimSetup::gatherInfo(void){ } // check for the temperature set flag - + if (globals->haveTempSet()) info[i].setTemp = globals->getTempSet(); - // get some of the tricky things that may still be in the globals + // check for the extended State init - double boxVector[3]; - if (globals->haveBox()){ - boxVector[0] = globals->getBox(); - boxVector[1] = globals->getBox(); - boxVector[2] = globals->getBox(); - - info[i].setBox(boxVector); - } - else if (globals->haveDensity()){ - double vol; - vol = (double) tot_nmol / globals->getDensity(); - boxVector[0] = pow(vol, (1.0 / 3.0)); - boxVector[1] = boxVector[0]; - boxVector[2] = boxVector[0]; - - info[i].setBox(boxVector); - } - else{ - if (!globals->haveBoxX()){ - sprintf(painCave.errMsg, - "SimSetup error, no periodic BoxX size given.\n"); - painCave.isFatal = 1; - simError(); - } - boxVector[0] = globals->getBoxX(); - - if (!globals->haveBoxY()){ - sprintf(painCave.errMsg, - "SimSetup error, no periodic BoxY size given.\n"); - painCave.isFatal = 1; - simError(); - } - boxVector[1] = globals->getBoxY(); - - if (!globals->haveBoxZ()){ - sprintf(painCave.errMsg, - "SimSetup error, no periodic BoxZ size given.\n"); - painCave.isFatal = 1; - simError(); - } - boxVector[2] = globals->getBoxZ(); - - info[i].setBox(boxVector); - } + info[i].useInitXSstate = globals->getUseInitXSstate(); + info[i].orthoTolerance = globals->getOrthoBoxTolerance(); + } - + //setup seed for random number generator int seedValue; @@ -813,18 +846,12 @@ void SimSetup::finalInfoCheck(void){ if (!globals->haveECR()){ sprintf(painCave.errMsg, - "SimSetup Warning: using default value of 1/2 the smallest " - "box length for the electrostaticCutoffRadius.\n" - "I hope you have a very fast processor!\n"); + "SimSetup Warning: No value was set for electrostaticCutoffRadius.\n" + "\tOOPSE will use a default value of 15.0 angstroms" + "\tfor the electrostaticCutoffRadius.\n"); painCave.isFatal = 0; simError(); - double smallest; - smallest = info[i].boxL[0]; - if (info[i].boxL[1] <= smallest) - smallest = info[i].boxL[1]; - if (info[i].boxL[2] <= smallest) - smallest = info[i].boxL[2]; - theEcr = 0.5 * smallest; + theEcr = 15.0; } else{ theEcr = globals->getECR(); @@ -832,8 +859,10 @@ void SimSetup::finalInfoCheck(void){ if (!globals->haveEST()){ sprintf(painCave.errMsg, - "SimSetup Warning: using default value of 0.05 * the " - "electrostaticCutoffRadius for the electrostaticSkinThickness\n"); + "SimSetup Warning: No value was set for electrostaticSkinThickness.\n" + "\tOOPSE will use a default value of\n" + "\t0.05 * electrostaticCutoffRadius\n" + "\tfor the electrostaticSkinThickness\n"); painCave.isFatal = 0; simError(); theEst = 0.05 * theEcr; @@ -842,12 +871,13 @@ void SimSetup::finalInfoCheck(void){ theEst = globals->getEST(); } - info[i].setEcr(theEcr, theEst); + info[i].setDefaultEcr(theEcr, theEst); if (!globals->haveDielectric()){ sprintf(painCave.errMsg, - "SimSetup Error: You are trying to use Reaction Field without" - "setting a dielectric constant!\n"); + "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(); } @@ -857,28 +887,23 @@ void SimSetup::finalInfoCheck(void){ if (usesDipoles){ if (!globals->haveECR()){ sprintf(painCave.errMsg, - "SimSetup Warning: using default value of 1/2 the smallest " - "box length for the electrostaticCutoffRadius.\n" - "I hope you have a very fast processor!\n"); - painCave.isFatal = 0; - simError(); - double smallest; - smallest = info[i].boxL[0]; - if (info[i].boxL[1] <= smallest) - smallest = info[i].boxL[1]; - if (info[i].boxL[2] <= smallest) - smallest = info[i].boxL[2]; - theEcr = 0.5 * smallest; + "SimSetup Warning: No value was set for electrostaticCutoffRadius.\n" + "\tOOPSE will use a default value of 15.0 angstroms" + "\tfor the electrostaticCutoffRadius.\n"); + painCave.isFatal = 0; + simError(); + theEcr = 15.0; } else{ theEcr = globals->getECR(); } - + if (!globals->haveEST()){ sprintf(painCave.errMsg, - "SimSetup Warning: using default value of 0.05 * the " - "electrostaticCutoffRadius for the " - "electrostaticSkinThickness\n"); + "SimSetup Warning: No value was set for electrostaticSkinThickness.\n" + "\tOOPSE will use a default value of\n" + "\t0.05 * electrostaticCutoffRadius\n" + "\tfor the electrostaticSkinThickness\n"); painCave.isFatal = 0; simError(); theEst = 0.05 * theEcr; @@ -886,18 +911,17 @@ void SimSetup::finalInfoCheck(void){ else{ theEst = globals->getEST(); } - - info[i].setEcr(theEcr, theEst); + + info[i].setDefaultEcr(theEcr, theEst); } } } - #ifdef IS_MPI strcpy(checkPointMsg, "post processing checks out"); MPIcheckPoint(); #endif // is_mpi } - + void SimSetup::initSystemCoords(void){ int i; @@ -914,7 +938,6 @@ void SimSetup::initSystemCoords(void){ if (worldRank == 0){ #endif //is_mpi inName = globals->getInitialConfig(); - double* tempDouble = new double[1000000]; fileInit = new InitializeFromFile(inName); #ifdef IS_MPI } @@ -926,21 +949,14 @@ void SimSetup::initSystemCoords(void){ delete fileInit; } else{ -#ifdef IS_MPI - + // no init from bass - + sprintf(painCave.errMsg, - "Cannot intialize a parallel simulation without an initial configuration file.\n"); - painCave.isFatal; + "Cannot intialize a simulation without an initial configuration file.\n"); + painCave.isFatal = 1;; simError(); - -#else - - initFromBass(); - - -#endif + } #ifdef IS_MPI @@ -1158,7 +1174,7 @@ void SimSetup::calcSysValues(void){ } void SimSetup::calcSysValues(void){ - int i, j, k; + int i; int* molMembershipArray; @@ -1236,8 +1252,8 @@ void SimSetup::mpiMolDivide(void){ if (local_atoms != info[0].n_atoms){ sprintf(painCave.errMsg, - "SimSetup error: mpiSim's localAtom (%d) and SimSetup's" - " localAtom (%d) are not equal.\n", + "SimSetup error: mpiSim's localAtom (%d) and SimSetup's\n" + "\tlocalAtom (%d) are not equal.\n", info[0].n_atoms, local_atoms); painCave.isFatal = 1; simError(); @@ -1257,7 +1273,11 @@ void SimSetup::makeSysArrays(void){ void SimSetup::makeSysArrays(void){ - int i, j, k, l; + +#ifndef IS_MPI + int k, j; +#endif // is_mpi + int i, l; Atom** the_atoms; Molecule* the_molecules; @@ -1340,19 +1360,24 @@ void SimSetup::makeIntegrator(void){ void SimSetup::makeIntegrator(void){ int k; + NVE* myNVE = NULL; NVT* myNVT = NULL; NPTi >* myNPTi = NULL; NPTf >* myNPTf = NULL; + NPTxyz >* myNPTxyz = NULL; for (k = 0; k < nInfo; k++){ switch (ensembleCase){ case NVE_ENS: if (globals->haveZconstraints()){ setupZConstraint(info[k]); - new ZConstraint >(&(info[k]), the_ff); + myNVE = new ZConstraint >(&(info[k]), the_ff); } - else - new NVE(&(info[k]), the_ff); + else{ + myNVE = new NVE(&(info[k]), the_ff); + } + + info->the_integrator = myNVE; break; case NVT_ENS: @@ -1370,10 +1395,12 @@ void SimSetup::makeIntegrator(void){ else{ sprintf(painCave.errMsg, "SimSetup error: If you use the NVT\n" - " ensemble, you must set tauThermostat.\n"); + "\tensemble, you must set tauThermostat.\n"); painCave.isFatal = 1; simError(); } + + info->the_integrator = myNVT; break; case NPTi_ENS: @@ -1391,7 +1418,7 @@ void SimSetup::makeIntegrator(void){ else{ sprintf(painCave.errMsg, "SimSetup error: If you use a constant pressure\n" - " ensemble, you must set targetPressure in the BASS file.\n"); + "\tensemble, you must set targetPressure in the BASS file.\n"); painCave.isFatal = 1; simError(); } @@ -1401,7 +1428,7 @@ void SimSetup::makeIntegrator(void){ else{ sprintf(painCave.errMsg, "SimSetup error: If you use an NPT\n" - " ensemble, you must set tauThermostat.\n"); + "\tensemble, you must set tauThermostat.\n"); painCave.isFatal = 1; simError(); } @@ -1411,10 +1438,12 @@ void SimSetup::makeIntegrator(void){ else{ sprintf(painCave.errMsg, "SimSetup error: If you use an NPT\n" - " ensemble, you must set tauBarostat.\n"); + "\tensemble, you must set tauBarostat.\n"); painCave.isFatal = 1; simError(); } + + info->the_integrator = myNPTi; break; case NPTf_ENS: @@ -1432,32 +1461,79 @@ void SimSetup::makeIntegrator(void){ else{ sprintf(painCave.errMsg, "SimSetup error: If you use a constant pressure\n" - " ensemble, you must set targetPressure in the BASS file.\n"); + "\tensemble, you must set targetPressure in the BASS file.\n"); painCave.isFatal = 1; simError(); } if (globals->haveTauThermostat()) myNPTf->setTauThermostat(globals->getTauThermostat()); + else{ sprintf(painCave.errMsg, "SimSetup error: If you use an NPT\n" - " ensemble, you must set tauThermostat.\n"); + "\tensemble, you must set tauThermostat.\n"); painCave.isFatal = 1; simError(); } if (globals->haveTauBarostat()) myNPTf->setTauBarostat(globals->getTauBarostat()); + else{ sprintf(painCave.errMsg, "SimSetup error: If you use an NPT\n" - " ensemble, you must set tauBarostat.\n"); + "\tensemble, you must set tauBarostat.\n"); painCave.isFatal = 1; simError(); } + + info->the_integrator = myNPTf; break; + case NPTxyz_ENS: + if (globals->haveZconstraints()){ + setupZConstraint(info[k]); + myNPTxyz = new ZConstraint > >(&(info[k]), the_ff); + } + else + myNPTxyz = new NPTxyz >(&(info[k]), the_ff); + + myNPTxyz->setTargetTemp(globals->getTargetTemp()); + + if (globals->haveTargetPressure()) + myNPTxyz->setTargetPressure(globals->getTargetPressure()); + else{ + sprintf(painCave.errMsg, + "SimSetup error: If you use a constant pressure\n" + "\tensemble, you must set targetPressure in the BASS file.\n"); + painCave.isFatal = 1; + simError(); + } + + if (globals->haveTauThermostat()) + myNPTxyz->setTauThermostat(globals->getTauThermostat()); + else{ + sprintf(painCave.errMsg, + "SimSetup error: If you use an NPT\n" + "\tensemble, you must set tauThermostat.\n"); + painCave.isFatal = 1; + simError(); + } + + if (globals->haveTauBarostat()) + myNPTxyz->setTauBarostat(globals->getTauBarostat()); + else{ + sprintf(painCave.errMsg, + "SimSetup error: If you use an NPT\n" + "\tensemble, you must set tauBarostat.\n"); + painCave.isFatal = 1; + simError(); + } + + info->the_integrator = myNPTxyz; + break; + default: sprintf(painCave.errMsg, "SimSetup Error. Unrecognized ensemble in case statement.\n"); @@ -1503,8 +1579,8 @@ void SimSetup::setupZConstraint(SimInfo& theInfo){ } else{ sprintf(painCave.errMsg, - "ZConstraint error: If you use an ZConstraint\n" - " , you must set sample time.\n"); + "ZConstraint error: If you use a ZConstraint,\n" + "\tyou must set zconsTime.\n"); painCave.isFatal = 1; simError(); } @@ -1519,8 +1595,9 @@ void SimSetup::setupZConstraint(SimInfo& theInfo){ else{ double defaultZConsTol = 0.01; sprintf(painCave.errMsg, - "ZConstraint Waring: Tolerance for z-constraint methodl is not specified\n" - " , default value %f is used.\n", + "ZConstraint Warning: Tolerance for z-constraint method is not specified.\n" + "\tOOPSE will use a default value of %f.\n" + "\tTo set the tolerance, use the zconsTol variable.\n", defaultZConsTol); painCave.isFatal = 0; simError(); @@ -1538,8 +1615,9 @@ void SimSetup::setupZConstraint(SimInfo& theInfo){ } else{ sprintf(painCave.errMsg, - "ZConstraint Warning: User does not set force Subtraction policy, " - "PolicyByMass is used\n"); + "ZConstraint Warning: No force subtraction policy was set.\n" + "\tOOPSE will use PolicyByMass.\n" + "\tTo set the policy, use the zconsForcePolicy variable.\n"); painCave.isFatal = 0; simError(); zconsForcePolicy->setData("BYMASS"); @@ -1583,7 +1661,7 @@ void SimSetup::setupZConstraint(SimInfo& theInfo){ //check the uniqueness of index if(!zconsParaData->isIndexUnique()){ sprintf(painCave.errMsg, - "ZConstraint Error: molIndex is not unique\n"); + "ZConstraint Error: molIndex is not unique!\n"); painCave.isFatal = 1; simError(); }