--- trunk/OOPSE/libmdtools/SimSetup.cpp 2004/01/19 21:17:39 965 +++ trunk/OOPSE/libmdtools/SimSetup.cpp 2004/03/16 19:22:56 1091 @@ -9,6 +9,8 @@ #include "parse_me.h" #include "Integrator.hpp" #include "simError.h" +//#include "ConjugateMinimizer.hpp" +#include "OOPSEMinimizer.hpp" #ifdef IS_MPI #include "mpiBASS.h" @@ -24,9 +26,10 @@ #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; @@ -144,10 +147,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 @@ -174,7 +180,6 @@ void SimSetup::makeMolecules(void){ bend_set* theBends; torsion_set* theTorsions; - //init the forceField paramters the_ff->readParams(); @@ -182,6 +187,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++){ @@ -218,10 +226,34 @@ 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); u = sqrt(uSqr); @@ -609,6 +641,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", force_field); @@ -811,9 +846,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 } @@ -1110,6 +1145,10 @@ void SimSetup::createFF(void){ the_ff = new EAM_FF(); break; + case FF_H2O: + the_ff = new WATER(); + break; + default: sprintf(painCave.errMsg, "SimSetup Error. Unrecognized force field in case statement.\n"); @@ -1624,7 +1663,26 @@ void SimSetup::setupZConstraint(SimInfo& theInfo){ } theInfo.addProperty(zconsForcePolicy); + + //set zcons gap + DoubleData* zconsGap = new DoubleData(); + zconsGap->setID(ZCONSGAP_ID); + + if (globals->haveZConsGap()){ + zconsGap->setData(globals->getZconsGap()); + theInfo.addProperty(zconsGap); + } + + //set zcons fixtime + DoubleData* zconsFixtime = new DoubleData(); + zconsFixtime->setID(ZCONSFIXTIME_ID); + if (globals->haveZConsFixTime()){ + zconsFixtime->setData(globals->getZconsFixtime()); + theInfo.addProperty(zconsFixtime); + } + + //Determine the name of ouput file and add it into SimInfo's property list //Be careful, do not use inFileName, since it is a pointer which //point to a string at master node, and slave nodes do not contain that string @@ -1672,3 +1730,73 @@ void SimSetup::setupZConstraint(SimInfo& theInfo){ //push data into siminfo, therefore, we can retrieve later theInfo.addProperty(zconsParaData); } + +void SimSetup::makeMinimizer(){ + + OOPSEMinimizer* myOOPSEMinimizer; + MinimizerParameterSet* param; + char minimizerName[100]; + + for (int i = 0; i < nInfo; i++){ + + //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->haveMinStepSize()){ + param->setStepSize(globals->getMinStepSize()); + } + + if (globals->haveMinLSMaxIter()){ + param->setLineSearchMaxIteration(globals->getMinLSMaxIter()); + } + + if (globals->haveMinLSTol()){ + param->setLineSearchTol(globals->getMinLSTol()); + } + + strcpy(minimizerName, globals->getMinimizer()); + + if (!strcasecmp(minimizerName, "CG")){ + myOOPSEMinimizer = new PRCGMinimizer(&(info[i]), the_ff, param); + } + else if (!strcasecmp(minimizerName, "SD")){ + //myOOPSEMinimizer = MinimizerFactory.creatMinimizer("", &(info[i]), the_ff, param); + myOOPSEMinimizer = new SDMinimizer(&(info[i]), the_ff, param); + } + else{ + sprintf(painCave.errMsg, + "SimSetup error: Unrecognized Minimizer, use Conjugate Gradient \n"); + painCave.isFatal = 0; + simError(); + + myOOPSEMinimizer = new PRCGMinimizer(&(info[i]), the_ff, param); + } + info[i].the_integrator = myOOPSEMinimizer; + + //store the minimizer into simInfo + info[i].the_minimizer = myOOPSEMinimizer; + info[i].has_minimizer = true; + } + +}