--- trunk/OOPSE/libmdtools/SimSetup.cpp 2003/10/27 22:07:48 823 +++ trunk/OOPSE/libmdtools/SimSetup.cpp 2004/02/06 21:37:59 1035 @@ -1,7 +1,7 @@ #include -#include +#include #include -#include +#include #include #include #include "SimSetup.hpp" @@ -9,6 +9,7 @@ #include "parse_me.h" #include "Integrator.hpp" #include "simError.h" +#include "ConjugateMinimizer.hpp" #ifdef IS_MPI #include "mpiBASS.h" @@ -24,12 +25,37 @@ #define NPTxyz_ENS 4 -#define FF_DUFF 0 -#define FF_LJ 1 -#define FF_EAM 2 +#define FF_DUFF 0 +#define FF_LJ 1 +#define FF_EAM 2 +#define FF_H2O 3 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; @@ -120,10 +146,13 @@ void SimSetup::createSim(void){ makeOutNames(); - // make the integrator - - makeIntegrator(); - + if (globals->haveMinimizer()) + // make minimizer + makeMinimizer(); + else + // make the integrator + makeIntegrator(); + #ifdef IS_MPI mpiSim->mpiRefresh(); #endif @@ -150,7 +179,6 @@ void SimSetup::makeMolecules(void){ bend_set* theBends; torsion_set* theTorsions; - //init the forceField paramters the_ff->readParams(); @@ -158,6 +186,9 @@ void SimSetup::makeMolecules(void){ // init the atoms + double phi, theta, psi; + double sux, suy, suz; + double Axx, Axy, Axz, Ayx, Ayy, Ayz, Azx, Azy, Azz; double ux, uy, uz, u, uSqr; for (k = 0; k < nInfo; k++){ @@ -194,9 +225,33 @@ void SimSetup::makeMolecules(void){ info[k].n_oriented++; molInfo.myAtoms[j] = dAtom; - ux = currentAtom->getOrntX(); - uy = currentAtom->getOrntY(); - uz = currentAtom->getOrntZ(); + // Directional Atoms have standard unit vectors which are oriented + // in space using the three Euler angles. We assume the standard + // unit vector was originally along the z axis below. + + phi = currentAtom->getEulerPhi() * M_PI / 180.0; + theta = currentAtom->getEulerTheta() * M_PI / 180.0; + psi = currentAtom->getEulerPsi()* M_PI / 180.0; + + Axx = (cos(phi) * cos(psi)) - (sin(phi) * cos(theta) * sin(psi)); + Axy = (sin(phi) * cos(psi)) + (cos(phi) * cos(theta) * sin(psi)); + Axz = sin(theta) * sin(psi); + + Ayx = -(cos(phi) * sin(psi)) - (sin(phi) * cos(theta) * cos(psi)); + Ayy = -(sin(phi) * sin(psi)) + (cos(phi) * cos(theta) * cos(psi)); + Ayz = sin(theta) * cos(psi); + + Azx = sin(phi) * sin(theta); + Azy = -cos(phi) * sin(theta); + Azz = cos(theta); + + sux = 0.0; + suy = 0.0; + suz = 1.0; + + ux = (Axx * sux) + (Ayx * suy) + (Azx * suz); + uy = (Axy * sux) + (Ayy * suy) + (Azy * suz); + uz = (Axz * sux) + (Ayz * suy) + (Azz * suz); uSqr = (ux * ux) + (uy * uy) + (uz * uz); @@ -584,6 +639,9 @@ void SimSetup::gatherInfo(void){ } else if (!strcasecmp(force_field, "EAM")){ ffCase = FF_EAM; + } + else if (!strcasecmp(force_field, "WATER")){ + ffCase = FF_H2O; } else{ sprintf(painCave.errMsg, "SimSetup Error. Unrecognized force field -> %s\n", @@ -613,8 +671,8 @@ void SimSetup::gatherInfo(void){ } 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(); @@ -646,8 +704,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(); } @@ -665,8 +723,49 @@ void SimSetup::gatherInfo(void){ " Please give nMol in the components.\n"); 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 for (i = 0; i < nInfo; i++){ @@ -696,58 +795,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; @@ -787,9 +845,9 @@ void SimSetup::gatherInfo(void){ for (int i = 0; i < nInfo; i++){ info[i].setSeed(seedValue); } - + #ifdef IS_MPI - strcpy(checkPointMsg, "Succesfully gathered all information from Bass\n"); + strcpy(checkPointMsg, "Successfully gathered all information from Bass\n"); MPIcheckPoint(); #endif // is_mpi } @@ -822,18 +880,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(); @@ -841,8 +893,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; @@ -851,12 +905,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(); } @@ -866,28 +921,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; @@ -895,18 +945,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; @@ -934,21 +983,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"); + "Cannot intialize a simulation without an initial configuration file.\n"); painCave.isFatal = 1;; simError(); - -#else - - initFromBass(); - - -#endif + } #ifdef IS_MPI @@ -1100,6 +1142,10 @@ void SimSetup::createFF(void){ case FF_EAM: the_ff = new EAM_FF(); + break; + + case FF_H2O: + the_ff = new WATER(); break; default: @@ -1244,8 +1290,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(); @@ -1387,7 +1433,7 @@ 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(); } @@ -1410,7 +1456,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(); } @@ -1420,7 +1466,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(); } @@ -1430,7 +1476,7 @@ 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(); } @@ -1453,27 +1499,29 @@ 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(); } @@ -1496,7 +1544,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(); } @@ -1506,7 +1554,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(); } @@ -1516,7 +1564,7 @@ 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(); } @@ -1569,8 +1617,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(); } @@ -1585,8 +1633,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(); @@ -1604,8 +1653,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"); @@ -1649,7 +1699,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(); } @@ -1660,3 +1710,80 @@ void SimSetup::setupZConstraint(SimInfo& theInfo){ //push data into siminfo, therefore, we can retrieve later theInfo.addProperty(zconsParaData); } + +void SimSetup::makeMinimizer(){ + + OOPSEMinimizerBase* myOOPSEMinimizerBase; + ObjFunctor1 * objFunc; + OutputFunctor* outputFunc; + ConcreteNLModel1* nlp; + MinimizerParameterSet* param; + ConjugateMinimizerBase* minimizer; + int dim; + + for (int i = 0; i < nInfo; i++){ + //creat + myOOPSEMinimizerBase = new OOPSEMinimizerBase(&(info[i]), the_ff); + + info[i].the_integrator = myOOPSEMinimizerBase; + //creat the object functor; + objFunc = (ObjFunctor1*) new ClassMemObjFunctor1 + (myOOPSEMinimizerBase, &OOPSEMinimizerBase::calcGradient); + + //creat output functor; + outputFunc = new ClassMemOutputFunctor + (myOOPSEMinimizerBase, &OOPSEMinimizerBase::output); + + //creat nonlinear model + dim = myOOPSEMinimizerBase->getDim(); + nlp = new ConcreteNLModel1(dim, objFunc); + + nlp->setX(myOOPSEMinimizerBase->getCoor()); + + //prepare parameter set for minimizer + param = new MinimizerParameterSet(); + param->setDefaultParameter(); + + if (globals->haveMinimizer()){ + param->setFTol(globals->getMinFTol()); + } + + if (globals->haveMinGTol()){ + param->setGTol(globals->getMinGTol()); + } + + if (globals->haveMinMaxIter()){ + param->setMaxIteration(globals->getMinMaxIter()); + } + + if (globals->haveMinWriteFrq()){ + param->setMaxIteration(globals->getMinMaxIter()); + } + + if (globals->haveMinWriteFrq()){ + param->setWriteFrq(globals->getMinWriteFrq()); + } + + if (globals->haveMinResetFrq()){ + param->setResetFrq(globals->getMinResetFrq()); + } + + if (globals->haveMinLSMaxIter()){ + param->setLineSearchMaxIteration(globals->getMinLSMaxIter()); + } + + if (globals->haveMinLSTol()){ + param->setLineSearchTol(globals->getMinLSTol()); + } + + //creat the minimizer + minimizer = new PRCGMinimizer(nlp, param); + minimizer->setLineSearchStrategy(nlp, GoldenSection); + minimizer->setOutputFunctor(outputFunc); + + //store the minimizer into simInfo + info[i].the_minimizer = minimizer; + info[i].has_minimizer = true; + } + +}