ICODE-ADC/lammps
4
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

Added forcezero linesearch method

Browse Source 
git-svn-id: svn://svn.icms.temple.edu/lammps-ro/trunk@7675 f3b2605a-c512-4ea7-a41b-209d697bcdaa
This commit is contained in:
athomps
2012-02-02 23:19:33 +00:00
parent 90e6bf8af0
commit 9d88ab7ee1
4 changed files with 404 additions and 21 deletions
Download Patch File Download Diff File Expand all files Collapse all files

420
src/min_linesearch.cpp
View File

@ -14,6 +14,8 @@
/* ----------------------------------------------------------------------
Contributing author: Aidan Thompson (SNL)
improved CG and backtrack ls, added quadratic ls
Contributing author: Asad Hasan (CMU)
added forcezero ls
Sources: Numerical Recipes frprmn routine
"Conjugate Gradient Method Without the Agonizing Pain" by
JR Shewchuk, http://www-2.cs.cmu.edu/~jrs/jrspapers.html#cg
@ -45,8 +47,8 @@ using namespace LAMMPS_NS;
#define ALPHA_MAX 1.0
#define ALPHA_REDUCE 0.5
#define BACKTRACK_SLOPE 0.4
#define EMACH 1.0e-8
#define QUADRATIC_TOL 0.1
#define EMACH 1.0e-8
#define EPS_QUAD 1.0e-28
// same as in other min classes
@ -79,6 +81,7 @@ void MinLineSearch::init_style()
{
if (linestyle == 0) linemin = &MinLineSearch::linemin_backtrack;
else if (linestyle == 1) linemin = &MinLineSearch::linemin_quadratic;
else if (linestyle == 2) linemin = &MinLineSearch::linemin_forcezero;
delete [] gextra;
delete [] hextra;
@ -160,12 +163,12 @@ void MinLineSearch::reset_vectors()
output args: return 0 if successful move, non-zero alpha
return non-zero if failed
alpha = distance moved along h for x at min eng config
update neval counter of eng/force function evaluations
output extra: if fail, energy_force() of original x
if succeed, energy_force() at x + alpha*h
update neval counter of eng/force function evaluations
output extra: if fail, energy_force() of original x
if succeed, energy_force() at x + alpha*h
atom->x = coords at new configuration
atom->f = force at new configuration
ecurrent = energy of new configuration
atom->f = force at new configuration
ecurrent = energy of new configuration
------------------------------------------------------------------------- */
/* ----------------------------------------------------------------------
@ -247,7 +250,7 @@ int MinLineSearch::linemin_backtrack(double eoriginal, double &alpha)
// double alphatmp = alphamax*1.0e-4;
// etmp = alpha_step(alphatmp,1);
// printf("alpha = %g dele = %g dele_force = %g err = %g\n",
// alphatmp,etmp-eoriginal,-alphatmp*fdothall,
// alphatmp,etmp-eoriginal,-alphatmp*fdothall,
// etmp-eoriginal+alphatmp*fdothall);
// alpha_step(0.0,1);
@ -326,7 +329,7 @@ int MinLineSearch::linemin_quadratic(double eoriginal, double &alpha)
double fdothall,fdothme,hme,hmax,hmaxall;
double de_ideal,de;
double delfh,engprev,relerr,alphaprev,fhprev,ff,fh,alpha0;
double dot[2],dotall[2];
double dot[2],dotall[2];
double *xatom,*x0atom,*fatom,*hatom;
double alphamax;
@ -405,7 +408,7 @@ int MinLineSearch::linemin_quadratic(double eoriginal, double &alpha)
// double alphatmp = alphamax*1.0e-4;
// etmp = alpha_step(alphatmp,1);
// printf("alpha = %g dele = %g dele_force = %g err = %g\n",
// alphatmp,etmp-eoriginal,-alphatmp*fdothall,
// alphatmp,etmp-eoriginal,-alphatmp*fdothall,
// etmp-eoriginal+alphatmp*fdothall);
// alpha_step(0.0,1);
@ -458,7 +461,7 @@ int MinLineSearch::linemin_quadratic(double eoriginal, double &alpha)
relerr = fabs(1.0-(0.5*(alpha-alphaprev)*(fh+fhprev)+ecurrent)/engprev);
alpha0 = alpha - (alpha-alphaprev)*fh/delfh;
if (relerr <= QUADRATIC_TOL && alpha0 > 0.0 && alpha0 < alphamax) {
if (relerr <= QUADRATIC_TOL && alpha0 > 0.0 && alpha0 < alphamax) {
ecurrent = alpha_step(alpha0,1);
if (ecurrent - eoriginal < EMACH) {
if (nextra_global) {
@ -502,6 +505,335 @@ int MinLineSearch::linemin_quadratic(double eoriginal, double &alpha)
}
}
/* ----------------------------------------------------------------------
forcezero linesearch method - seeks a zero of force in a robust manner.
(motivated by a line minimization routine of f77 DYNAMO code)
central idea:
In each linesearch we attempt to converge to a zero of force
(usual case) or reduces forces (worst case).
Energy does not play any role in the search procedure,
except we ensure that it doesn't increase.
pseudo code:
i) Fix an alpha max:
// also account for nextra atom & global
alpha_max <= dmax/hmaxall
ii) Initialize:
fhCurr = current_force.dot.search_direction
fhoriginal = fhCurr
// try decreasing the energy to 1/10 of initial
alpha_init = 0.1*fabs(eoriginal)/fhCurr;
// initial alpha is smaller than alpha_max
alpha_del = MIN(alpha_init, 0.5*alpha_max);
alpha = 0.0
iii) Loop:
backtrack = false
alpha += alpha_del
if (alpha > alpha_max):
// we have done enough in the search space
EXIT with success
Step with the new alpha
Compute:
current energy and 'fhCurr'
de = ecurrent - eprev
// ZERO_ENERGY = 1e-12, is max allowed energy increase
if (de > ZERO_ENERGY):
bactrack = true
// GRAD_TOL = 0.1
if ( (not backtrack) && (fabs(fhCurr/fh0) <= GRAD_TOL) ):
// forces sufficiently reduced without energy increase
EXIT with success
// projected force changed sign but didn't become small enough
if ( fhCurr < 0):
backtrack = true
if (bactrack):
// forces along search direction changed sign
if (fhCurr < 0):
Get alpha_del by solving for zero
of force (1D Newton's Method)
else:
// force didn't change sign but only energy increased,
// we overshot a minimum which is very close to a
// maximum (or there is an inflection point)
// New alpha_del should be much smaller
// ALPHA_FACT = 0.1
alpha_del *= ALPHA_FACT
// Check to see if new 'alpha_del' isn't too small
if (alpha_del < MIN_ALPHA):
EXIT with failure("linesearch alpha is zero")
Undo the step of alpha.
// continue the loop with a new alpha_del
else:
Get new alpha_del by linearizing force and solving for its zero
---------------------------------------------------------------------- */
int MinLineSearch::linemin_forcezero(double eoriginal, double &alpha)
{
int i,m,n;
double fdothall,fdothme,hme,hmax,hmaxall;
double de_ideal,de;
double *xatom,*x0atom,*fatom,*hatom;
double alpha_max, alpha_init, alpha_del;
// projection of: force on itself, current force on search direction,
double ffCurr, fhCurr;
// previous force on search direction, initial force on search direction
double fhPrev, fhoriginal;
// current energy, previous energy
double engCurr, engPrev;
bool backtrack;
// hardcoded constants
// factor by which alpha is reduced when backtracking
double ALPHA_FACT = 0.1;
// maximum amount by which we'll permit energy increase
double ZERO_ENERGY = 1e-12;
// fraction to which we want to reduce the directional derivative
double GRAD_TOL = 0.1;
// largest alpha increment which will trigger a failed_linesearch
double MIN_ALPHA_FAC = 1e-20;
double LIMIT_BOOST = 4.0;
// fdothall = projection of search dir along downhill gradient
// if search direction is not downhill, exit with error
fdothme = 0.0;
for (i = 0; i < nvec; i++) fdothme += fvec[i]*h[i];
if (nextra_atom)
for (m = 0; m < nextra_atom; m++) {
fatom = fextra_atom[m];
hatom = hextra_atom[m];
n = extra_nlen[m];
for (i = 0; i < n; i++)
fdothme += fatom[i]*hatom[i];
}
MPI_Allreduce(&fdothme,&fdothall,1,MPI_DOUBLE,MPI_SUM,world);
if (nextra_global)
for (i = 0; i < nextra_global; i++)
fdothall += fextra[i]*hextra[i];
if (output->thermo->normflag) fdothall /= atom->natoms;
if (fdothall <= 0.0) return DOWNHILL;
// set alpha so no dof is changed by more than max allowed amount
// for atom coords, max amount = dmax
// for extra per-atom dof, max amount = extra_max[]
// for extra global dof, max amount is set by fix
// also insure alpha <= ALPHA_MAX else will have
// to backtrack from huge value when forces are tiny
// if all search dir components are already 0.0, exit with error
hme = 0.0;
for (i = 0; i < nvec; i++)
hme = MAX(hme,fabs(h[i]));
MPI_Allreduce(&hme,&hmaxall,1,MPI_DOUBLE,MPI_MAX,world);
alpha_max = dmax/hmaxall;
if (nextra_atom)
for (m = 0; m < nextra_atom; m++) {
hatom = hextra_atom[m];
n = extra_nlen[m];
hme = 0.0;
for (i = 0; i < n; i++) hme = MAX(hme,fabs(hatom[i]));
MPI_Allreduce(&hme,&hmax,1,MPI_DOUBLE,MPI_MAX,world);
alpha_max = MIN(alpha_max,extra_max[m]/hmax);
hmaxall = MAX(hmaxall,hmax);
}
if (nextra_global) {
double alpha_extra = modify->max_alpha(hextra);
alpha_max = MIN(alpha_max, alpha_extra);
for (i = 0; i < nextra_global; i++)
hmaxall = MAX(hmaxall,fabs(hextra[i]));
}
if (hmaxall == 0.0) return ZEROFORCE;
// store box and values of all dof at start of linesearch
fix_minimize->store_box();
for (i = 0; i < nvec; i++) x0[i] = xvec[i];
if (nextra_atom)
for (m = 0; m < nextra_atom; m++) {
xatom = xextra_atom[m];
x0atom = x0extra_atom[m];
n = extra_nlen[m];
for (i = 0; i < n; i++) x0atom[i] = xatom[i];
}
if (nextra_global) modify->min_store();
// initialize important variables before main linesearch loop
ffCurr = 0.0;
fhCurr = fdothall;
fhoriginal = fhCurr;
engCurr = eoriginal;
// stores energy difference due to the current move
de = 0.0;
// choosing the initial alpha that we'll use
// rough estimate that'll decrease energy to 1/10
alpha_init = 0.1*fabs(eoriginal)/fdothall;
// initialize aplha to 0.0
alpha = 0.0;
// compute increment to alpha, ensure that we
// don't take the largest allowed alpha
// first alpha that will actually apply
alpha_del = MIN(alpha_init,0.5*alpha_max);
// main linesearch loop
while (1) {
backtrack = false;
fhPrev = fhCurr;
engPrev = engCurr;
// apply the increment to alpha, but first
// check whether we are still in allowed search space
alpha += alpha_del;
if (alpha > alpha_max) {
// undo the increment
alpha -= alpha_del;
if (nextra_global) {
int itmp = modify->min_reset_ref();
if (itmp) ecurrent = energy_force(1);
}
// exit linesearch with success: have done
// enough in allowed search space
return 0;
}
// move the system
// '1' updates coordinates of atoms which cross PBC
engCurr = alpha_step(alpha,1);
ecurrent = engCurr;
// compute the new directional derivative and also f_dot_f
fhCurr = compute_dir_deriv(ffCurr);
// energy change
de = engCurr - engPrev;
// if the function value increases measurably,
// then we have to reduce alpha
if (de >= ZERO_ENERGY)
backtrack = true;
// check if the directional derivative has sufficiently decreased
// NOTE: the fabs is essential here
if ((!backtrack) && (fabs(fhCurr/fhoriginal) <= GRAD_TOL)) {
if (nextra_global) {
int itmp = modify->min_reset_ref();
if (itmp) ecurrent = energy_force(1);
}
// we are done
return 0;
}
// check if the directional derivative changed sign
// but it's not small: we overshot the minima -- BACKTRACK
if (fhCurr < 0.0)
backtrack = true;
// backtrack by undoing step and choosing a new alpha
if (backtrack) {
// move back
alpha -= alpha_del;
// choose new alpha
// if the force changed sign, linearize force and
// solve for new alpha_del
if (fhCurr < 0.0)
alpha_del *= fhPrev/(fhPrev - fhCurr);
else
// force didn't change sign but only energy increased,
// we overshot a minimum which is very close to a maxima
// (or there is an inflection point)
// new alpha_del should be much smaller
alpha_del *= ALPHA_FACT;
// since we moved back ...
engCurr = engPrev;
ecurrent = engCurr;
fhCurr = fhPrev;
// if new move is too small then we have failed;
// exit with 'failed_linesearch'
if (hmaxall*alpha_del <= MIN_ALPHA_FAC) {
// undo all line minization moves
engCurr = alpha_step(0.0,1);
ecurrent= engCurr;
return ZEROALPHA;
}
} else {
// get a new alpha by linearizing force and start over
double boostFactor = LIMIT_BOOST;
// avoids problems near an energy inflection point
if (fhPrev > fhCurr)
boostFactor = fhCurr/(fhPrev - fhCurr);
// don't want to boost too much
boostFactor = MIN(boostFactor, LIMIT_BOOST);
alpha_del *= boostFactor;
}
}
}
/* ---------------------------------------------------------------------- */
double MinLineSearch::alpha_step(double alpha, int resetflag)
@ -515,11 +847,11 @@ double MinLineSearch::alpha_step(double alpha, int resetflag)
for (i = 0; i < nvec; i++) xvec[i] = x0[i];
if (nextra_atom)
for (m = 0; m < nextra_atom; m++) {
xatom = xextra_atom[m];
x0atom = x0extra_atom[m];
n = extra_nlen[m];
for (i = 0; i < n; i++) xatom[i] = x0atom[i];
requestor[m]->min_x_set(m);
xatom = xextra_atom[m];
x0atom = x0extra_atom[m];
n = extra_nlen[m];
for (i = 0; i < n; i++) xatom[i] = x0atom[i];
requestor[m]->min_x_set(m);
}
// step forward along h
@ -529,11 +861,11 @@ double MinLineSearch::alpha_step(double alpha, int resetflag)
for (i = 0; i < nvec; i++) xvec[i] += alpha*h[i];
if (nextra_atom)
for (m = 0; m < nextra_atom; m++) {
xatom = xextra_atom[m];
hatom = hextra_atom[m];
n = extra_nlen[m];
for (i = 0; i < n; i++) xatom[i] += alpha*hatom[i];
requestor[m]->min_x_set(m);
xatom = xextra_atom[m];
hatom = hextra_atom[m];
n = extra_nlen[m];
for (i = 0; i < n; i++) xatom[i] += alpha*hatom[i];
requestor[m]->min_x_set(m);
}
}
@ -542,3 +874,51 @@ double MinLineSearch::alpha_step(double alpha, int resetflag)
neval++;
return energy_force(resetflag);
}
/* ---------------------------------------------------------------------- */
// compute projection of force on: itself and the search direction
double MinLineSearch::compute_dir_deriv(double &ff)
{
int i,m,n;
double *xatom,*hatom, *fatom;
double dot[2],dotall[2];
double fh;
// compute new fh, alpha, delfh
dot[0] = dot[1] = 0.0;
for (i = 0; i < nvec; i++) {
dot[0] += fvec[i]*fvec[i];
dot[1] += fvec[i]*h[i];
}
if (nextra_atom)
for (m = 0; m < nextra_atom; m++) {
fatom = fextra_atom[m];
hatom = hextra_atom[m];
n = extra_nlen[m];
for (i = 0; i < n; i++) {
dot[0] += fatom[i]*fatom[i];
dot[1] += fatom[i]*hatom[i];
}
}
MPI_Allreduce(dot,dotall,2,MPI_DOUBLE,MPI_SUM,world);
if (nextra_global) {
for (i = 0; i < nextra_global; i++) {
dotall[0] += fextra[i]*fextra[i];
dotall[1] += fextra[i]*hextra[i];
}
}
ff = dotall[0];
fh = dotall[1];
if (output->thermo->normflag) {
ff /= atom->natoms;
fh /= atom->natoms;
}
return fh;
}