git-svn-id: svn://svn.icms.temple.edu/lammps-ro/trunk@3199 f3b2605a-c512-4ea7-a41b-209d697bcdaa

src/ASPHERE/fix_nve_asphere.cpp

@@ -27,9 +27,6 @@
 #include "memory.h"
 #include "error.h"
 
-#define TOLERANCE 1.0e-6
-#define EPSILON 1.0e-7
-
 using namespace LAMMPS_NS;
 
 /* ---------------------------------------------------------------------- */

src/fix.h

@@ -109,6 +109,9 @@ class Fix : protected Pointers {
 
   virtual double min_energy(double *) {return 0.0;}
   virtual void min_store() {}
+  virtual void min_clearstore() {}
+  virtual void min_pushstore() {}
+  virtual void min_popstore() {}
   virtual void min_step(double, double *) {}
   virtual double max_alpha(double *) {return 0.0;}
   virtual int min_dof() {return 0;}

src/fix_box_relax.cpp

@@ -179,6 +179,8 @@ FixBoxRelax::FixBoxRelax(LAMMPS *lmp, int narg, char **arg) :
   dimension = domain->dimension;
   nrigid = 0;
   rfix = 0;
+
+  current_lifo = 0;
 }
 
 /* ---------------------------------------------------------------------- */
@@ -301,20 +303,60 @@ double FixBoxRelax::min_energy(double *fextra)
 }
 
 /* ----------------------------------------------------------------------
-   store extra dof values for linesearch starting point
+   store extra dof values for minimization linesearch starting point
    boxlo0,boxhi0 = box dimensions
    s0 = ratio of current boxsize to initial boxsize
+   box values are pushed onto a LIFO stack so nested calls can be made
+   values are popped by calling min_step(0.0)
 ------------------------------------------------------------------------- */
 
 void FixBoxRelax::min_store()
 {
   for (int i = 0; i < 3; i++) {
-    boxlo0[i] = domain->boxlo[i];
-    boxhi0[i] = domain->boxhi[i];
+    boxlo0[current_lifo][i] = domain->boxlo[i];
+    boxhi0[current_lifo][i] = domain->boxhi[i];
   }
-  s0[0] = (boxhi0[0]-boxlo0[0])/xprdinit;
-  s0[1] = (boxhi0[1]-boxlo0[1])/yprdinit;
-  s0[2] = (boxhi0[2]-boxlo0[2])/zprdinit;
+  s0[0] = (boxhi0[current_lifo][0]-boxlo0[current_lifo][0])/xprdinit;
+  s0[1] = (boxhi0[current_lifo][1]-boxlo0[current_lifo][1])/yprdinit;
+  s0[2] = (boxhi0[current_lifo][2]-boxlo0[current_lifo][2])/zprdinit;
+}
+
+/* ----------------------------------------------------------------------
+   clear the LIFO stack for min_store
+------------------------------------------------------------------------- */
+
+void FixBoxRelax::min_clearstore()
+{
+  current_lifo = 0;
+}
+
+/* ----------------------------------------------------------------------
+   push the LIFO stack for min_store
+------------------------------------------------------------------------- */
+
+void FixBoxRelax::min_pushstore()
+{
+  if (current_lifo >= MAX_LIFO_DEPTH) {
+    error->all("Attempt to push beyond stack limit <FixBoxRelax>");
+    return;
+  }
+
+  current_lifo++;
+}
+
+
+/* ----------------------------------------------------------------------
+   pop the LIFO stack for min_store
+------------------------------------------------------------------------- */
+
+void FixBoxRelax::min_popstore()
+{
+  if (current_lifo <= 0) {
+    error->all("Attempt to pop empty stack <FixBoxRelax>");
+    return;
+  }
+
+  current_lifo--;
 }
 
 /* ----------------------------------------------------------------------
@@ -394,9 +436,11 @@ void FixBoxRelax::remap()
 
   for (i = 0; i < 3; i++)
     if (p_flag[i]) {
-      ctr = 0.5 * (boxlo0[i] + boxhi0[i]);
-      domain->boxlo[i] = boxlo0[i] + (boxlo0[i]-ctr)*ds[i]/s0[i];
-      domain->boxhi[i] = boxhi0[i] + (boxhi0[i]-ctr)*ds[i]/s0[i];
+      double currentBoxLo0 = boxlo0[current_lifo][i];
+      double currentBoxHi0 = boxhi0[current_lifo][i];
+      ctr = 0.5 * (currentBoxLo0 + currentBoxHi0);
+      domain->boxlo[i] = currentBoxLo0 + (currentBoxLo0-ctr)*ds[i]/s0[i];
+      domain->boxhi[i] = currentBoxHi0 + (currentBoxHi0-ctr)*ds[i]/s0[i];
     }
 
   domain->set_global_box();

src/fix_box_relax.h

@@ -16,6 +16,8 @@
 
 #include "fix.h"
 
+const int  MAX_LIFO_DEPTH = 2;    // to meet the needs of min_hftn
+
 namespace LAMMPS_NS {
 
 class FixBoxRelax : public Fix {
@@ -27,6 +29,9 @@ class FixBoxRelax : public Fix {
 
   double min_energy(double *);
   void min_store();
+  void min_clearstore();
+  void min_pushstore();
+  void min_popstore();
   void min_step(double, double *);
   double max_alpha(double *);
   int min_dof();
@@ -42,9 +47,11 @@ class FixBoxRelax : public Fix {
   double volinit,xprdinit,yprdinit,zprdinit;
   double vmax,pv2e,pflagsum;
 
-  double boxlo0[3],boxhi0[3];     // box bounds at start of line search
-  double s0[6];	          	  // scale matrix in Voigt notation
-  double ds[6];                   // increment in scale matrix
+  int    current_lifo;                  // LIFO stack pointer
+  double boxlo0[MAX_LIFO_DEPTH][3];     // box bounds at start of line search
+  double boxhi0[MAX_LIFO_DEPTH][3];
+  double s0[6];                         // scale matrix in Voigt notation
+  double ds[6];                         // increment in scale matrix
 
   char *id_temp,*id_press;
   class Compute *temperature,*pressure;

src/min.cpp

@@ -304,11 +304,16 @@ void Min::run(int nsteps)
 {
   // possible stop conditions
 
-  char *stopstrings[] = {"max iterations","max force evaluations",
-                         "energy tolerance","force tolerance",
-			 "search direction is not downhill",
-			 "linesearch alpha is zero",
-			 "forces are zero","quadratic factors are zero"};
+  char *stopstrings[] = {"max iterations",
+                         "max force evaluations",
+                         "energy tolerance",
+                         "force tolerance",
+                         "search direction is not downhill",
+                         "linesearch alpha is zero",
+                         "forces are zero",
+                         "quadratic factors are zero",
+                         "trust region too small",
+                         "HFTN minimizer error"};
 
   // stats for Finish to print
 

src/min_hftn.h

@@ -0,0 +1,119 @@
+/* ----------------------------------------------------------------------
+   LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
+   http://lammps.sandia.gov, Sandia National Laboratories
+   Steve Plimpton, sjplimp@sandia.gov
+
+   Copyright (2003) Sandia Corporation.  Under the terms of Contract
+   DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
+   certain rights in this software.  This software is distributed under 
+   the GNU General Public License.
+
+   See the README file in the top-level LAMMPS directory.
+------------------------------------------------------------------------- */
+
+#ifndef MIN_HFTN_H
+#define MIN_HFTN_H
+
+#include "min.h"
+
+namespace LAMMPS_NS
+{
+  //---- THE ALGORITHM NEEDS TO STORE THIS MANY ATOM-BASED VECTORS,
+  //---- IN ADDITION TO ATOM POSITIONS AND THE FORCE VECTOR.
+
+  static const int  NUM_HFTN_ATOM_BASED_VECTORS = 7;
+  static const int  VEC_XK    = 0;   //-- ATOM POSITIONS AT SUBITER START
+  static const int  VEC_CG_P  = 1;   //-- STEP p IN CG SUBITER
+  static const int  VEC_CG_D  = 2;   //-- DIRECTION d IN CG SUBITER
+  static const int  VEC_CG_HD = 3;   //-- HESSIAN-VECTOR PRODUCT Hd
+  static const int  VEC_CG_R  = 4;   //-- RESIDUAL r IN CG SUBITER
+  static const int  VEC_DIF1  = 5;   //-- FOR FINITE DIFFERENCING
+  static const int  VEC_DIF2  = 6;   //-- FOR FINITE DIFFERENCING
+
+class MinHFTN : public Min
+{
+  public:
+  
+  MinHFTN (LAMMPS *);
+  ~MinHFTN (void);
+  
+  void init_style();
+  void setup_style();
+  void reset_vectors();
+  int iterate (int);
+
+  private:
+
+  //---- ATOM-BASED STORAGE VECTORS.
+  double *   _daAVectors[NUM_HFTN_ATOM_BASED_VECTORS];
+  double **  _daExtraAtom[NUM_HFTN_ATOM_BASED_VECTORS];
+  
+  //---- GLOBAL DOF STORAGE.  ELEMENT [0] IS X0 (XK), NOT USED IN THIS ARRAY.
+  double *   _daExtraGlobal[NUM_HFTN_ATOM_BASED_VECTORS];
+  
+  int     _nNumUnknowns;
+  FILE *  _fpPrint;
+  
+  int   execute_hftn_ (const bool      bPrintProgress,
+		       const double    dInitialEnergy,
+		       const double    dInitialForce2,
+		       double &  dFinalEnergy,
+		       double &  dFinalForce2);
+  bool    compute_inner_cg_step_ (const double    dTrustRadius,
+				  const double    dForceTol,
+				  const int       nMaxEvals,
+				  const bool      bHaveEvalAtXin,
+				  const double    dEnergyAtXin,
+				  const double    dForce2AtXin,
+				  double &  dEnergyAtXout,
+				  double &  dForce2AtXout,
+				  int    &  nStepType,
+				  double &  dStepLength2,
+				  double &  dStepLengthInf);
+  double  calc_xinf_using_mpi_ (void) const;
+  double  calc_dot_prod_using_mpi_ (const int  nIx1,
+				    const int  nIx2) const;
+  double  calc_grad_dot_v_using_mpi_ (const int  nIx) const;
+  void    calc_dhd_dd_using_mpi_ (double &  dDHD,
+				  double &  dDD) const;
+  void    calc_ppnew_pdold_using_mpi_ (double &  dPnewDotPnew,
+				       double &  dPoldDotD) const;
+  void    calc_plengths_using_mpi_ (double &  dStepLength2,
+				    double &  dStepLengthInf) const;
+  bool    step_exceeds_TR_ (const double    dTrustRadius,
+			    const double    dPP,
+			    const double    dPD,
+			    const double    dDD,
+			    double &  dTau) const;
+  bool    step_exceeds_DMAX_ (void) const;
+  void    adjust_step_to_tau_ (const double  tau);
+  double  compute_to_tr_ (const double  dPP,
+			  const double  dPD,
+			  const double  dDD,
+			  const double  dTrustRadius,
+			  const bool    bConsiderBothRoots,
+			  const double  dDHD,
+			  const double  dPdotHD,
+			  const double  dGradDotD) const;
+  void    evaluate_dir_der_ (const bool      bUseForwardDiffs,
+			     const int       nIxDir,
+			     const int       nIxResult,
+			     const bool      bEvaluateAtX,
+			     double &  dNewEnergy);
+  void  open_hftn_print_file_ (void);
+  void  hftn_print_line_ (const bool    bIsStepAccepted,
+			  const int     nIteration,
+			  const int     nTotalEvals,
+			  const double  dEnergy,
+			  const double  dForce2,
+			  const int     nStepType,
+			  const double  dTrustRadius,
+			  const double  dStepLength2,
+			  const double  dActualRed,
+			  const double  dPredictedRed) const;
+  void  close_hftn_print_file_ (void);
+};
+ 
+}
+
+#endif

src/modify.cpp

@@ -440,6 +440,28 @@ void Modify::min_store()
     fix[list_min_energy[i]]->min_store();
 }
 
+/* ----------------------------------------------------------------------
+   mange state of extra dof on a stack, only for relevant fixes
+------------------------------------------------------------------------- */
+
+void Modify::min_clearstore()
+{
+  for (int i = 0; i < n_min_energy; i++)
+    fix[list_min_energy[i]]->min_clearstore();
+}
+
+void Modify::min_pushstore()
+{
+  for (int i = 0; i < n_min_energy; i++)
+    fix[list_min_energy[i]]->min_pushstore();
+}
+
+void Modify::min_popstore()
+{
+  for (int i = 0; i < n_min_energy; i++)
+    fix[list_min_energy[i]]->min_popstore();
+}
+
 /* ----------------------------------------------------------------------
    displace extra dof along vector hextra, only for relevant fixes
 ------------------------------------------------------------------------- */

src/modify.h

@@ -65,6 +65,9 @@ class Modify : protected Pointers {
   double min_energy(double *);
   void min_store();
   void min_step(double, double *);
+  void min_clearstore();
+  void min_pushstore();
+  void min_popstore();
   double max_alpha(double *);
   int min_dof();
 

src/style.h

@@ -297,11 +297,13 @@ IntegrateStyle(verlet,Verlet)
 
 #ifdef MinimizeInclude
 #include "min_cg.h"
+#include "min_hftn.h"
 #include "min_sd.h"
 #endif
 
 #ifdef MinimizeClass
 MinimizeStyle(cg,MinCG)
+MinimizeStyle(hftn,MinHFTN)
 MinimizeStyle(sd,MinSD)
 # endif