--- trunk/OOPSE/libmdtools/SimSetup.cpp	2004/01/19 21:17:39	965
+++ trunk/OOPSE/libmdtools/SimSetup.cpp	2004/02/06 21:37:59	1035
@@ -9,6 +9,7 @@
 #include "parse_me.h"
 #include "Integrator.hpp"
 #include "simError.h"
+#include "ConjugateMinimizer.hpp"
 
 #ifdef IS_MPI
 #include "mpiBASS.h"
@@ -24,9 +25,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 +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
@@ -174,7 +179,6 @@ void SimSetup::makeMolecules(void){
   bend_set* theBends;
   torsion_set* theTorsions;
 
-
   //init the forceField paramters
 
   the_ff->readParams();
@@ -182,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++){
@@ -218,10 +225,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 +640,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 +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
 }
@@ -1110,6 +1144,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");
@@ -1672,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<OOPSEMinimizerBase>
+                                              (myOOPSEMinimizerBase, &OOPSEMinimizerBase::calcGradient);
+
+    //creat output functor;
+    outputFunc =  new ClassMemOutputFunctor<OOPSEMinimizerBase>
+                               (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;
+  }
+
+}