From b1181dd9005f4dff1563d7ff87f9f0fc2a8e9ea3 Mon Sep 17 00:00:00 2001 From: sjplimp Date: Mon, 14 Apr 2008 20:50:48 +0000 Subject: [PATCH] git-svn-id: svn://svn.icms.temple.edu/lammps-ro/trunk@1743 f3b2605a-c512-4ea7-a41b-209d697bcdaa --- src/finish.cpp | 46 +++++---- src/min.cpp | 22 +---- src/min.h | 4 +- src/min_cg.cpp | 237 +++++++++++++--------------------------------- src/min_cg.h | 9 +- src/min_cg_fr.cpp | 108 --------------------- src/min_cg_fr.h | 29 ------ src/min_sd.cpp | 43 +++++---- src/min_sd.h | 2 +- src/minimize.cpp | 12 ++- src/style.h | 2 - src/update.h | 2 +- 12 files changed, 140 insertions(+), 376 deletions(-) delete mode 100644 src/min_cg_fr.cpp delete mode 100644 src/min_cg_fr.h diff --git a/src/finish.cpp b/src/finish.cpp index 1240c7e85a..269d53b501 100644 --- a/src/finish.cpp +++ b/src/finish.cpp @@ -88,29 +88,43 @@ void Finish::end(int flag) if (me == 0) { if (screen) { fprintf(screen,"Minimization stats:\n"); - fprintf(screen," E initial, next-to-last, final = %g %g %g\n", + fprintf(screen," Stopping criterion = %s\n", + update->minimize->stopstr); + fprintf(screen," Energy initial, next-to-last, final = \n" + " %18.12g %18.12g %18.12g\n", update->minimize->einitial,update->minimize->eprevious, update->minimize->efinal); - fprintf(screen," Gradient 2-norm init/final= %g %g\n", - update->minimize->gnorm2_init,update->minimize->gnorm2_final); - fprintf(screen," Gradient inf-norm init/final= %g %g\n", - update->minimize->gnorminf_init, - update->minimize->gnorminf_final); - fprintf(screen," Iterations = %d\n",update->minimize->niter); - fprintf(screen," Force evaluations = %d\n",update->minimize->neval); + fprintf(screen," Force two-norm initial, final = %g %g\n", + update->minimize->fnorm2_init,update->minimize->fnorm2_final); + fprintf(screen," Force max component initial, final = %g %g\n", + update->minimize->fnorminf_init, + update->minimize->fnorminf_final); + fprintf(screen," Final line search alpha, max atom move = %g %g\n", + update->minimize->alpha_final, + update->minimize->alpha_final* + update->minimize->fnorminf_final); + fprintf(screen," Iterations, force evaluations = %d %d\n", + update->minimize->niter,update->minimize->neval); } if (logfile) { fprintf(logfile,"Minimization stats:\n"); - fprintf(logfile," E initial, next-to-last, final = %g %g %g\n", + fprintf(logfile," Stopping criterion = %s\n", + update->minimize->stopstr); + fprintf(logfile," Energy initial, next-to-last, final = \n" + " %18.12g %18.12g %18.12g\n", update->minimize->einitial,update->minimize->eprevious, update->minimize->efinal); - fprintf(logfile," Gradient 2-norm init/final= %g %g\n", - update->minimize->gnorm2_init,update->minimize->gnorm2_final); - fprintf(logfile," Gradient inf-norm init/final= %g %g\n", - update->minimize->gnorminf_init, - update->minimize->gnorminf_final); - fprintf(logfile," Iterations = %d\n",update->minimize->niter); - fprintf(logfile," Force evaluations = %d\n",update->minimize->neval); + fprintf(logfile," Force two-norm initial, final = %g %g\n", + update->minimize->fnorm2_init,update->minimize->fnorm2_final); + fprintf(logfile," Force max component initial, final = %g %g\n", + update->minimize->fnorminf_init, + update->minimize->fnorminf_final); + fprintf(logfile," Final line search alpha, max atom move = %g %g\n", + update->minimize->alpha_final, + update->minimize->alpha_final* + update->minimize->fnorminf_final); + fprintf(logfile," Iterations, force evaluations = %d %d\n", + update->minimize->niter,update->minimize->neval); } } if (me == 0) { diff --git a/src/min.cpp b/src/min.cpp index 08e406429a..bf141fd8f7 100644 --- a/src/min.cpp +++ b/src/min.cpp @@ -20,17 +20,11 @@ using namespace LAMMPS_NS; -#define SCAN 0 // same as in min_cg.cpp -#define SECANT 1 - /* ---------------------------------------------------------------------- */ Min::Min(LAMMPS *lmp) : Pointers(lmp) { - linestyle = SECANT; - dmin = 1.0e-5; dmax = 0.1; - lineiter = 10; elist_atom = NULL; vlist_global = vlist_atom = NULL; @@ -53,24 +47,10 @@ void Min::modify_params(int narg, char **arg) int iarg = 0; while (iarg < narg) { - if (strcmp(arg[iarg],"linestyle") == 0) { - if (iarg+2 > narg) error->all("Illegal min_modify command"); - if (strcmp(arg[iarg+1],"scan") == 0) linestyle = SCAN; - else if (strcmp(arg[iarg+1],"secant") == 0) linestyle = SECANT; - else error->all("Illegal min_modify command"); - iarg += 2; - } else if (strcmp(arg[iarg],"dmin") == 0) { - if (iarg+2 > narg) error->all("Illegal min_modify command"); - dmin = atof(arg[iarg+1]); - iarg += 2; - } else if (strcmp(arg[iarg],"dmax") == 0) { + if (strcmp(arg[iarg],"dmax") == 0) { if (iarg+2 > narg) error->all("Illegal min_modify command"); dmax = atof(arg[iarg+1]); iarg += 2; - } else if (strcmp(arg[iarg],"lineiter") == 0) { - if (iarg+2 > narg) error->all("Illegal min_modify command"); - lineiter = atoi(arg[iarg+1]); - iarg += 2; } else error->all("Illegal min_modify command"); } } diff --git a/src/min.h b/src/min.h index 2fddc3e0d5..35ad2d1753 100644 --- a/src/min.h +++ b/src/min.h @@ -21,8 +21,10 @@ namespace LAMMPS_NS { class Min : protected Pointers { public: double einitial,efinal,eprevious; - double gnorm2_init,gnorminf_init,gnorm2_final,gnorminf_final; + double fnorm2_init,fnorminf_init,fnorm2_final,fnorminf_final; + double alpha_final; int niter,neval; + char *stopstr; double dmin,dmax; int linestyle,lineiter; diff --git a/src/min_cg.cpp b/src/min_cg.cpp index 7723729068..f73c9745be 100644 --- a/src/min_cg.cpp +++ b/src/min_cg.cpp @@ -48,12 +48,15 @@ using namespace LAMMPS_NS; #define MIN(A,B) ((A) < (B)) ? (A) : (B) #define MAX(A,B) ((A) > (B)) ? (A) : (B) -#define EPS 1.0e-6 -#define SCAN_FACTOR 2.0 -#define SECANT_EPS 1.0e-3 +#define EPS_ENERGY 1.0e-8 +#define BACKTRACK_SLOPE 0.5 +#define ALPHA_REDUCE 0.5 -#define SCAN 0 // same as in min.cpp -#define SECANT 1 +enum{FAIL,MAXITER,MAXEVAL,ETOL,FTOL}; // same as in other min classes + +char *stopstrings[] = {"failed linesearch","max iterations", + "max force evaluations","energy tolerance", + "force tolerance"}; /* ---------------------------------------------------------------------- */ @@ -117,8 +120,7 @@ void MinCG::init() // set ptr to linemin function - if (linestyle == SCAN) linemin = &MinCG::linemin_scan; - else if (linestyle == SECANT) linemin = &MinCG::linemin_secant; + linemin = &MinCG::linemin_backtrack; } /* ---------------------------------------------------------------------- @@ -149,17 +151,18 @@ void MinCG::run() f = atom->f[0]; tmp = 0.0; for (int i = 0; i < ndof; i++) tmp += f[i]*f[i]; - MPI_Allreduce(&tmp,&gnorm2_init,1,MPI_DOUBLE,MPI_SUM,world); - gnorm2_init = sqrt(gnorm2_init); + MPI_Allreduce(&tmp,&fnorm2_init,1,MPI_DOUBLE,MPI_SUM,world); + fnorm2_init = sqrt(fnorm2_init); tmp = 0.0; for (int i = 0; i < ndof; i++) tmp = MAX(fabs(f[i]),tmp); - MPI_Allreduce(&tmp,&gnorminf_init,1,MPI_DOUBLE,MPI_MAX,world); + MPI_Allreduce(&tmp,&fnorminf_init,1,MPI_DOUBLE,MPI_MAX,world); // minimizer iterations timer->barrier_start(TIME_LOOP); - iterate(update->nsteps); + int stop_condition = iterate(update->nsteps); + stopstr = stopstrings[stop_condition]; // account for early exit from iterate loop due to convergence // set niter/nsteps for Finish stats to print @@ -206,12 +209,12 @@ void MinCG::run() f = atom->f[0]; tmp = 0.0; for (int i = 0; i < ndof; i++) tmp += f[i]*f[i]; - MPI_Allreduce(&tmp,&gnorm2_final,1,MPI_DOUBLE,MPI_SUM,world); - gnorm2_final = sqrt(gnorm2_final); + MPI_Allreduce(&tmp,&fnorm2_final,1,MPI_DOUBLE,MPI_SUM,world); + fnorm2_final = sqrt(fnorm2_final); tmp = 0.0; for (int i = 0; i < ndof; i++) tmp = MAX(fabs(f[i]),tmp); - MPI_Allreduce(&tmp,&gnorminf_final,1,MPI_DOUBLE,MPI_MAX,world); + MPI_Allreduce(&tmp,&fnorminf_final,1,MPI_DOUBLE,MPI_MAX,world); } /* ---------------------------------------------------------------------- @@ -269,10 +272,10 @@ void MinCG::setup() Polak-Ribiere formulation ------------------------------------------------------------------------- */ -void MinCG::iterate(int n) +int MinCG::iterate(int n) { - int i,gradsearch,fail,ntimestep; - double alpha,beta,gg,dot[2],dotall[2]; + int i,fail,ntimestep; + double beta,gg,dot[2],dotall[2]; double *f; f = atom->f[0]; @@ -283,7 +286,6 @@ void MinCG::iterate(int n) MPI_Allreduce(dot,&gg,1,MPI_DOUBLE,MPI_SUM,world); neval = 0; - gradsearch = 1; for (niter = 0; niter < n; niter++) { @@ -292,26 +294,21 @@ void MinCG::iterate(int n) // line minimization along direction h from current atom->x eprevious = ecurrent; - fail = (this->*linemin)(ndof,atom->x[0],h,ecurrent,dmin,dmax,alpha,neval); + fail = (this->*linemin)(ndof,atom->x[0],h,ecurrent,dmin,dmax, + alpha_final,neval); + if (fail) return FAIL; - // if max_eval exceeded, all done - // if linemin failed or energy did not decrease sufficiently: - // if searched in grad direction, then all done - // else force next search to be in grad direction (CG restart) + // function evaluation criterion - if (neval >= update->max_eval) break; + if (neval >= update->max_eval) return MAXEVAL; - if (fail || fabs(ecurrent-eprevious) <= - update->tolerance * 0.5*(fabs(ecurrent) + fabs(eprevious) + EPS)) { - if (gradsearch == 1) break; - gradsearch = -1; - } + // energy tolerance criterion - // update h from new f = -Grad(x) and old g - // old g,h must have migrated with atoms to do this correctly - // done if size sq of grad vector < EPS - // force new search dir to be grad dir if need to restart CG - // set gradsearch to 1 if will search in grad dir on next iteration + if (fabs(ecurrent-eprevious) <= + update->etol * 0.5*(fabs(ecurrent) + fabs(eprevious) + EPS_ENERGY)) + return ETOL; + + // force tolerance criterion f = atom->f[0]; dot[0] = dot[1] = 0.0; @@ -321,13 +318,12 @@ void MinCG::iterate(int n) } MPI_Allreduce(dot,dotall,2,MPI_DOUBLE,MPI_SUM,world); + if (dotall[0] < update->ftol * update->ftol) return FTOL; + + // update h from new f = -Grad(x) and old g + beta = MAX(0.0,(dotall[0] - dotall[1])/gg); gg = dotall[0]; - if (gg < EPS) break; - - if (gradsearch == -1) beta = 0.0; - if (beta == 0.0) gradsearch = 1; - else gradsearch = 0; for (i = 0; i < ndof; i++) { g[i] = f[i]; @@ -342,6 +338,8 @@ void MinCG::iterate(int n) timer->stamp(TIME_OUTPUT); } } + + return MAXITER; } /* ---------------------------------------------------------------------- @@ -481,11 +479,11 @@ void MinCG::force_clear() dir = search direction as vector eng = current energy at initial x min/max dist = min/max distance to move any atom coord - output: return 0 if successful move, set alpha - return 1 if failed, no move, no need to set alpha - alpha = distance moved along dir to set x to min-eng config + output: return 0 if successful move, non-zero alpha alpha + return 1 if failed, alpha = 0.0 + alpha = distance moved along dir to set x to minimun eng config caller has several quantities set via last call to eng_force() - INSURE last call to eng_force() is consistent with returns + must insure last call to eng_force() is consistent with returns if fail, eng_force() of original x if succeed, eng_force() at x + alpha*dir atom->x = coords at new configuration @@ -493,24 +491,30 @@ void MinCG::force_clear() ecurrent = energy of new configuration NOTE: when call eng_force: n,x,dir,eng may change due to atom migration updated values are returned by eng_force() - these routines CANNOT store atom-based quantities b/c of migration + b/c of migration, linemin routines CANNOT store atom-based quantities ------------------------------------------------------------------------- */ /* ---------------------------------------------------------------------- - linemin: scan forward by larger and larger steps (SCAN_FACTOR) + linemin: backtracking line search (Alg 3.1, p 41 in Nocedal and Wright) uses no gradient info, but should be very robust - start at mindist, continue until maxdist - quit as soon as energy starts to rise + start at maxdist, backtrack until energy decrease is sufficient ------------------------------------------------------------------------- */ -int MinCG::linemin_scan(int n, double *x, double *dir, double eng, - double mindist, double maxdist, - double &alpha, int &nfunc) +int MinCG::linemin_backtrack(int n, double *x, double *dir, double eng, + double mindist, double maxdist, + double &alpha, int &nfunc) { - int i; - double fmax,fme,elowest,alphamin,alphamax,alphalast; + int i,ilist,m; + double fdotdirme,fdotdirall,fme,fmax,eoriginal,alphamax; + + // stopping criterion + + double *f = atom->f[0]; + fdotdirme = 0.0; + for (i = 0; i < n; i++) fdotdirme += f[i]*dir[i]; + MPI_Allreduce(&fdotdirme,&fdotdirall,1,MPI_DOUBLE,MPI_SUM,world); + if (output->thermo->normflag) fdotdirall /= atom->natoms; - // alphamin = step that moves some atom coord by mindist // alphamax = step that moves some atom coord by maxdist fme = 0.0; @@ -518,133 +522,26 @@ int MinCG::linemin_scan(int n, double *x, double *dir, double eng, MPI_Allreduce(&fme,&fmax,1,MPI_DOUBLE,MPI_MAX,world); if (fmax == 0.0) return 1; - alphamin = mindist/fmax; alphamax = maxdist/fmax; - // if minstep is already uphill, fail - // if eng increases, stop and return previous alpha - // if alphamax, stop and return alphamax + // reduce alpha by factor of ALPHA_REDUCE until energy decrease is sufficient - elowest = eng; - alpha = alphamin; + eoriginal = eng; + alpha = alphamax; while (1) { for (i = 0; i < n; i++) x[i] += alpha*dir[i]; eng_force(&n,&x,&dir,&eng); nfunc++; - if (alpha == alphamin && eng >= elowest) { - for (i = 0; i < n; i++) x[i] -= alpha*dir[i]; + if (eng <= eoriginal - BACKTRACK_SLOPE*alpha*fdotdirall) return 0; + + for (i = 0; i < n; i++) x[i] -= alpha*dir[i]; + + alpha *= ALPHA_REDUCE; + if (alpha == 0.0) { eng_force(&n,&x,&dir,&eng); - nfunc++; return 1; } - if (eng > elowest) { - for (i = 0; i < n; i++) x[i] += (alphalast-alpha)*dir[i]; - eng_force(&n,&x,&dir,&eng); - nfunc++; - alpha = alphalast; - return 0; - } - if (alpha == alphamax) return 0; - - elowest = eng; - alphalast = alpha; - alpha *= SCAN_FACTOR; - if (alpha > alphamax) alpha = alphamax; - } -} - -/* ---------------------------------------------------------------------- - linemin: use secant approximation to estimate parabola minimum at each step - should converge more quickly/accurately than "scan", but can be less robust -------------------------------------------------------------------------- */ - -int MinCG::linemin_secant(int n, double *x, double *dir, double eng, - double mindist, double maxdist, - double &alpha, int &nfunc) -{ - int i,iter; - double eta,eta_prev,alphamin,alphamax,alphadelta,fme,fmax,dsq,e0,tmp; - double *f; - double epssq = SECANT_EPS * SECANT_EPS; - - // stopping criterion for secant iterations - - fme = 0.0; - for (i = 0; i < n; i++) fme += dir[i]*dir[i]; - MPI_Allreduce(&fme,&dsq,1,MPI_DOUBLE,MPI_SUM,world); - - // alphamin = smallest allowed step of mindist - // alphamax = largest allowed step (in single iteration) of maxdist - - fme = 0.0; - for (i = 0; i < n; i++) fme = MAX(fme,fabs(dir[i])); - MPI_Allreduce(&fme,&fmax,1,MPI_DOUBLE,MPI_MAX,world); - if (fmax == 0.0) return 1; - - alphamin = mindist/fmax; - alphamax = maxdist/fmax; - - // eval func at alphamin - // exit if minstep is already uphill - - e0 = eng; - for (i = 0; i < n; i++) x[i] += alphamin*dir[i]; - eng_force(&n,&x,&dir,&eng); - nfunc++; - - if (eng >= e0) { - for (i = 0; i < n; i++) x[i] -= alphamin*dir[i]; - eng_force(&n,&x,&dir,&eng); - nfunc++; - return 1; - } - - // secant iterations - // alphadelta = new increment to move, alpha = accumulated move - // first step is alpha = 0, first previous step is at mindist - // prevent func evals for alpha outside mindist to maxdist - // if happens on 1st iteration and alpha < mindist - // secant approx is likely searching - // for a maximum (negative alpha), so reevaluate at alphamin - // if happens on 1st iteration and alpha > maxdist - // wants to take big step, so reevaluate at alphamax - - f = atom->f[0]; - tmp = 0.0; - for (i = 0; i < n; i++) tmp -= f[i]*dir[i]; - MPI_Allreduce(&tmp,&eta_prev,1,MPI_DOUBLE,MPI_SUM,world); - - alpha = alphamin; - alphadelta = -alphamin; - - for (iter = 0; iter < lineiter; iter++) { - alpha += alphadelta; - for (i = 0; i < n; i++) x[i] += alphadelta*dir[i]; - eng_force(&n,&x,&dir,&eng); - nfunc++; - - f = atom->f[0]; - tmp = 0.0; - for (i = 0; i < n; i++) tmp -= f[i]*dir[i]; - MPI_Allreduce(&tmp,&eta,1,MPI_DOUBLE,MPI_SUM,world); - - alphadelta *= eta / (eta_prev - eta); - eta_prev = eta; - if (alphadelta*alphadelta*dsq <= epssq) break; - - if (alpha+alphadelta < alphamin || alpha+alphadelta > alphamax) { - if (iter == 0) { - if (alpha+alphadelta < alphamin) alpha = alphamin; - else alpha = alphamax; - for (i = 0; i < n; i++) x[i] += alpha*dir[i]; - eng_force(&n,&x,&dir,&eng); - nfunc++; - } - break; - } - } - - return 0; + } } diff --git a/src/min_cg.h b/src/min_cg.h index 586f13362b..8db0ef262b 100644 --- a/src/min_cg.h +++ b/src/min_cg.h @@ -24,7 +24,7 @@ class MinCG : public Min { virtual ~MinCG() {} void init(); void run(); - virtual void iterate(int); + virtual int iterate(int); protected: int pairflag,torqueflag; @@ -42,11 +42,8 @@ class MinCG : public Min { typedef int (MinCG::*FnPtr)(int, double *, double *, double, double, double, double &, int &); FnPtr linemin; // ptr to linemin functions - - int linemin_scan(int, double *, double *, double, - double, double, double &, int &); - int linemin_secant(int, double *, double *, double, - double, double, double &, int &); + int linemin_backtrack(int, double *, double *, double, + double, double, double &, int &); void setup(); void setup_vectors(); diff --git a/src/min_cg_fr.cpp b/src/min_cg_fr.cpp deleted file mode 100644 index fbac308558..0000000000 --- a/src/min_cg_fr.cpp +++ /dev/null @@ -1,108 +0,0 @@ -/* ---------------------------------------------------------------------- - 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. -------------------------------------------------------------------------- */ - -#include "math.h" -#include "mpi.h" -#include "min_cg_fr.h" -#include "atom.h" -#include "update.h" -#include "output.h" -#include "timer.h" - -using namespace LAMMPS_NS; - -#define MIN(A,B) ((A) < (B)) ? (A) : (B) -#define MAX(A,B) ((A) > (B)) ? (A) : (B) - -#define EPS 1.0e-6 - -/* ---------------------------------------------------------------------- */ - -MinCGFR::MinCGFR(LAMMPS *lmp) : MinCG(lmp) {} - -/* ---------------------------------------------------------------------- - minimization via conjugate gradient iterations - Fletcher-Reeves formulation -------------------------------------------------------------------------- */ - -void MinCGFR::iterate(int n) -{ - int i,gradsearch,fail; - double alpha,beta,gg,dot,dotall; - double *f; - - f = atom->f[0]; - for (int i = 0; i < ndof; i++) h[i] = g[i] = f[i]; - - dot = 0.0; - for (i = 0; i < ndof; i++) dot += f[i]*f[i]; - MPI_Allreduce(&dot,&gg,1,MPI_DOUBLE,MPI_SUM,world); - - neval = 0; - gradsearch = 1; - - for (niter = 0; niter < n; niter++) { - - update->ntimestep++; - - // line minimization along direction h from current atom->x - - eprevious = ecurrent; - fail = (this->*linemin)(ndof,atom->x[0],h,ecurrent,dmin,dmax,alpha,neval); - - // if max_eval exceeded, all done - // if linemin failed or energy did not decrease sufficiently: - // all done if searched in grad direction - // else force next search to be in grad direction (CG restart) - - if (neval >= update->max_eval) break; - - if (fail || fabs(ecurrent-eprevious) <= - update->tolerance * 0.5*(fabs(ecurrent) + fabs(eprevious) + EPS)) { - if (gradsearch == 1) break; - gradsearch = -1; - } - - // update h from new f = -Grad(x) and old g - // old g,h must have migrated with atoms to do this correctly - // done if size sq of grad vector < EPS - // force new search dir to be grad dir if need to restart CG - // set gradsesarch to 1 if will search in grad dir on next iteration - - f = atom->f[0]; - dot = 0.0; - for (i = 0; i < ndof; i++) dot += f[i]*f[i]; - MPI_Allreduce(&dot,&dotall,1,MPI_DOUBLE,MPI_SUM,world); - - beta = dotall/gg; - gg = dotall; - if (gg < EPS) break; - - if (gradsearch == -1) beta = 0.0; - if (beta == 0.0) gradsearch = 1; - else gradsearch = 0; - - for (i = 0; i < ndof; i++) { - g[i] = f[i]; - h[i] = g[i] + beta*h[i]; - } - - // output for thermo, dump, restart files - - if (output->next == update->ntimestep) { - timer->stamp(); - output->write(update->ntimestep); - timer->stamp(TIME_OUTPUT); - } - } -} diff --git a/src/min_cg_fr.h b/src/min_cg_fr.h deleted file mode 100644 index d82dfeae7e..0000000000 --- a/src/min_cg_fr.h +++ /dev/null @@ -1,29 +0,0 @@ -/* ---------------------------------------------------------------------- - 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_CG_FR_H -#define MIN_CG_FR_H - -#include "min_cg.h" - -namespace LAMMPS_NS { - -class MinCGFR : public MinCG { - public: - MinCGFR(class LAMMPS *lmp); - void iterate(int); -}; - -} - -#endif diff --git a/src/min_sd.cpp b/src/min_sd.cpp index 83a678af10..07975049dc 100644 --- a/src/min_sd.cpp +++ b/src/min_sd.cpp @@ -21,7 +21,9 @@ using namespace LAMMPS_NS; -#define EPS 1.0e-6 +#define EPS 1.0e-8 + +enum{FAIL,MAXITER,MAXEVAL,ETOL,FTOL}; // same as in other min classes /* ---------------------------------------------------------------------- */ @@ -31,10 +33,10 @@ MinSD::MinSD(LAMMPS *lmp) : MinCG(lmp) {} minimization via steepest descent ------------------------------------------------------------------------- */ -void MinSD::iterate(int n) +int MinSD::iterate(int n) { - int i,fail; - double alpha,dot,dotall; + int i,fail,ntimestep; + double dot,dotall; double *f; f = atom->f[0]; @@ -44,39 +46,46 @@ void MinSD::iterate(int n) for (niter = 0; niter < n; niter++) { - update->ntimestep++; + ntimestep = ++update->ntimestep; // line minimization along direction h from current atom->x eprevious = ecurrent; - fail = (this->*linemin)(ndof,atom->x[0],h,ecurrent,dmin,dmax,alpha,neval); + fail = (this->*linemin)(ndof,atom->x[0],h,ecurrent,dmin,dmax, + alpha_final,neval); + if (fail) return FAIL; - // if max_eval exceeded, all done - // if linemin failed or energy did not decrease sufficiently, all done + // function evaluation criterion - if (neval >= update->max_eval) break; + if (neval >= update->max_eval) return MAXEVAL; - if (fail || fabs(ecurrent-eprevious) <= - update->tolerance * 0.5*(fabs(ecurrent) + fabs(eprevious) + EPS)) - break; + // energy tolerance criterion - // set h to new f = -Grad(x) - // done if size sq of grad vector < EPS + if (fabs(ecurrent-eprevious) <= + update->etol * 0.5*(fabs(ecurrent) + fabs(eprevious) + EPS)) + return ETOL; + + // force tolerance criterion f = atom->f[0]; dot = 0.0; for (i = 0; i < ndof; i++) dot += f[i]*f[i]; MPI_Allreduce(&dot,&dotall,1,MPI_DOUBLE,MPI_SUM,world); - if (dotall < EPS) break; + + if (dotall < update->ftol * update->ftol) return FTOL; + + // set h to new f = -Grad(x) for (i = 0; i < ndof; i++) h[i] = f[i]; // output for thermo, dump, restart files - if (output->next == update->ntimestep) { + if (output->next == ntimestep) { timer->stamp(); - output->write(update->ntimestep); + output->write(ntimestep); timer->stamp(TIME_OUTPUT); } } + + return MAXITER; } diff --git a/src/min_sd.h b/src/min_sd.h index 70f7add92c..be502162d4 100644 --- a/src/min_sd.h +++ b/src/min_sd.h @@ -21,7 +21,7 @@ namespace LAMMPS_NS { class MinSD : public MinCG { public: MinSD(class LAMMPS *); - void iterate(int); + int iterate(int); }; } diff --git a/src/minimize.cpp b/src/minimize.cpp index 944f5831aa..a2a1b82acf 100644 --- a/src/minimize.cpp +++ b/src/minimize.cpp @@ -29,14 +29,18 @@ Minimize::Minimize(LAMMPS *lmp) : Pointers(lmp) {} void Minimize::command(int narg, char **arg) { - if (narg != 3) error->all("Illegal minimize command"); + if (narg != 4) error->all("Illegal minimize command"); if (domain->box_exist == 0) error->all("Minimize command before simulation box is defined"); - update->tolerance = atof(arg[0]); - update->nsteps = atoi(arg[1]); - update->max_eval = atoi(arg[2]); + update->etol = atof(arg[0]); + update->ftol = atof(arg[1]); + update->nsteps = atoi(arg[2]); + update->max_eval = atoi(arg[3]); + + if (update->etol < 0.0 || update->ftol < 0.0) + error->all("Illegal minimize command"); update->beginstep = update->firststep = update->ntimestep; update->endstep = update->laststep = update->firststep + update->nsteps; diff --git a/src/style.h b/src/style.h index a255fdef8c..f65c316da7 100644 --- a/src/style.h +++ b/src/style.h @@ -279,13 +279,11 @@ IntegrateStyle(verlet,Verlet) #ifdef MinimizeInclude #include "min_cg.h" -#include "min_cg_fr.h" #include "min_sd.h" #endif #ifdef MinimizeClass MinimizeStyle(cg,MinCG) -MinimizeStyle(cg/fr,MinCGFR) MinimizeStyle(sd,MinSD) # endif diff --git a/src/update.h b/src/update.h index 4f0a8556ad..1e5847cbe8 100644 --- a/src/update.h +++ b/src/update.h @@ -21,7 +21,7 @@ namespace LAMMPS_NS { class Update : protected Pointers { public: double dt; // timestep - double tolerance; // minimizer tolerance + double etol,ftol; // minimizer tolerances on energy/force int ntimestep; // current step (dynamics or min iter) int nsteps; // # of steps to run (dynamics or min iter) int whichflag; // 0 for time integration, 1 for minimization