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

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

Browse Source
This commit is contained in:
sjplimp
2009-04-03 19:46:03 +00:00
parent 5265bf5e87
commit 90dacbb8a5
6 changed files with 775 additions and 765 deletions
Download Patch File Download Diff File Expand all files Collapse all files

723
src/min.cpp
View File

@ -11,16 +11,61 @@
See the README file in the top-level LAMMPS directory.
------------------------------------------------------------------------- */
/* ----------------------------------------------------------------------
Contributing author: Aidan Thompson (SNL)
improved CG and backtrack ls, added quadratic 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
------------------------------------------------------------------------- */
#include "math.h"
#include "stdlib.h"
#include "string.h"
#include "min.h"
#include "atom.h"
#include "domain.h"
#include "comm.h"
#include "update.h"
#include "modify.h"
#include "fix_minimize.h"
#include "compute.h"
#include "neighbor.h"
#include "force.h"
#include "pair.h"
#include "bond.h"
#include "angle.h"
#include "dihedral.h"
#include "improper.h"
#include "kspace.h"
#include "output.h"
#include "thermo.h"
#include "timer.h"
#include "error.h"
using namespace LAMMPS_NS;
// ALPHA_MAX = max alpha allowed to avoid long backtracks
// ALPHA_REDUCE = reduction ratio, should be in range [0.5,1)
// BACKTRACK_SLOPE, should be in range (0,0.5]
// QUADRATIC_TOL = tolerance on alpha0, should be in range [0.1,1)
// IDEAL_TOL = ideal energy tolerance for backtracking
// EPS_QUAD = tolerance for quadratic projection
#define ALPHA_MAX 1.0
#define ALPHA_REDUCE 0.5
#define BACKTRACK_SLOPE 0.4
#define QUADRATIC_TOL 0.1
#define IDEAL_TOL 1.0e-8
#define EPS_QUAD 1.0e-28
// same as in other min classes
enum{MAXITER,MAXEVAL,ETOL,FTOL,DOWNHILL,ZEROALPHA,ZEROFORCE,ZEROQUAD};
#define MIN(A,B) ((A) < (B)) ? (A) : (B)
#define MAX(A,B) ((A) > (B)) ? (A) : (B)
/* ---------------------------------------------------------------------- */
Min::Min(LAMMPS *lmp) : Pointers(lmp)
@ -30,6 +75,8 @@ Min::Min(LAMMPS *lmp) : Pointers(lmp)
elist_atom = NULL;
vlist_global = vlist_atom = NULL;
fextra = gextra = hextra = NULL;
}
/* ---------------------------------------------------------------------- */
@ -39,6 +86,682 @@ Min::~Min()
delete [] elist_atom;
delete [] vlist_global;
delete [] vlist_atom;
delete [] fextra;
delete [] gextra;
delete [] hextra;
}
/* ---------------------------------------------------------------------- */
void Min::init()
{
// create fix needed for storing atom-based gradient vectors
// will delete it at end of run
char **fixarg = new char*[3];
fixarg[0] = (char *) "MINIMIZE";
fixarg[1] = (char *) "all";
fixarg[2] = (char *) "MINIMIZE";
modify->add_fix(3,fixarg);
delete [] fixarg;
fix_minimize = (FixMinimize *) modify->fix[modify->nfix-1];
// zero gradient vectors before first atom exchange
setup_vectors();
for (int i = 0; i < ndof; i++) h[i] = g[i] = 0.0;
// virial_style:
// 1 if computed explicitly by pair->compute via sum over pair interactions
// 2 if computed implicitly by pair->virial_compute via sum over ghost atoms
if (force->newton_pair) virial_style = 2;
else virial_style = 1;
// setup lists of computes for global and per-atom PE and pressure
ev_setup();
// set flags for what arrays to clear in force_clear()
// clear torques if array exists
torqueflag = 0;
if (atom->torque) torqueflag = 1;
// orthogonal vs triclinic simulation box
triclinic = domain->triclinic;
// reset reneighboring criteria if necessary
neigh_every = neighbor->every;
neigh_delay = neighbor->delay;
neigh_dist_check = neighbor->dist_check;
if (neigh_every != 1 || neigh_delay != 0 || neigh_dist_check != 1) {
if (comm->me == 0)
error->warning("Resetting reneighboring criteria during minimization");
}
neighbor->every = 1;
neighbor->delay = 0;
neighbor->dist_check = 1;
// set ptr to linemin function
if (linestyle == 0) linemin = &Min::linemin_backtrack;
else if (linestyle == 1) linemin = &Min::linemin_quadratic;
}
/* ----------------------------------------------------------------------
setup before run
------------------------------------------------------------------------- */
void Min::setup()
{
if (comm->me == 0 && screen) fprintf(screen,"Setting up minimization ...\n");
// setup domain, communication and neighboring
// acquire ghosts
// build neighbor lists
// reset gradient vector ptrs
if (triclinic) domain->x2lamda(atom->nlocal);
domain->pbc();
domain->reset_box();
comm->setup();
if (neighbor->style) neighbor->setup_bins();
comm->exchange();
comm->borders();
if (triclinic) domain->lamda2x(atom->nlocal+atom->nghost);
neighbor->build();
neighbor->ncalls = 0;
setup_vectors();
// compute all forces
ev_set(update->ntimestep);
force_clear();
if (force->pair) force->pair->compute(eflag,vflag);
if (atom->molecular) {
if (force->bond) force->bond->compute(eflag,vflag);
if (force->angle) force->angle->compute(eflag,vflag);
if (force->dihedral) force->dihedral->compute(eflag,vflag);
if (force->improper) force->improper->compute(eflag,vflag);
}
if (force->kspace) {
force->kspace->setup();
force->kspace->compute(eflag,vflag);
}
if (force->newton) comm->reverse_communicate();
modify->setup(vflag);
output->setup(1);
}
/* ----------------------------------------------------------------------
perform minimization, with setup first
------------------------------------------------------------------------- */
void Min::run()
{
int i;
double tmp,*f;
// 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"};
// set initial force & energy
setup();
// setup any extra dof due to fixes
// can't be done until now b/c update init() comes before modify init()
delete [] fextra;
delete [] gextra;
delete [] hextra;
nextra = modify->min_dof();
if (nextra) {
fextra = new double[nextra];
gextra = new double[nextra];
hextra = new double[nextra];
}
// compute potential energy of system
// normalize energy if thermo PE does
int id = modify->find_compute("thermo_pe");
if (id < 0) error->all("Minimization could not find thermo_pe compute");
pe_compute = modify->compute[id];
ecurrent = pe_compute->compute_scalar();
if (nextra) ecurrent += modify->min_energy(fextra);
if (output->thermo->normflag) ecurrent /= atom->natoms;
// stats for Finish to print
einitial = ecurrent;
f = NULL;
if (ndof) f = atom->f[0];
tmp = 0.0;
for (i = 0; i < ndof; i++) tmp += f[i]*f[i];
MPI_Allreduce(&tmp,&fnorm2_init,1,MPI_DOUBLE,MPI_SUM,world);
if (nextra)
for (i = 0; i < nextra; i++) fnorm2_init += fextra[i]*fextra[i];
fnorm2_init = sqrt(fnorm2_init);
tmp = 0.0;
for (i = 0; i < ndof; i++) tmp = MAX(fabs(f[i]),tmp);
MPI_Allreduce(&tmp,&fnorminf_init,1,MPI_DOUBLE,MPI_MAX,world);
if (nextra)
for (i = 0; i < nextra; i++)
fnorminf_init = MAX(fabs(fextra[i]),fnorminf_init);
// minimizer iterations
timer->barrier_start(TIME_LOOP);
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
// set output->next values to this timestep
// call eng_force to insure vflag is set when forces computed
// output->write does final output for thermo, dump, restart files
// add ntimestep to all computes that store invocation times
// since are hardwireing call to thermo/dumps and computes may not be ready
if (niter < update->nsteps) {
niter++;
update->nsteps = niter;
for (int idump = 0; idump < output->ndump; idump++)
output->next_dump[idump] = update->ntimestep;
output->next_dump_any = update->ntimestep;
if (output->restart_every) output->next_restart = update->ntimestep;
output->next_thermo = update->ntimestep;
modify->addstep_compute_all(update->ntimestep);
int ntmp;
double *xtmp,*htmp,*x0tmp,etmp;
eng_force(&ntmp,&xtmp,&htmp,&x0tmp,&etmp,0);
output->write(update->ntimestep);
}
timer->barrier_stop(TIME_LOOP);
// delete fix at end of run, so its atom arrays won't persist
modify->delete_fix("MINIMIZE");
// reset reneighboring criteria
neighbor->every = neigh_every;
neighbor->delay = neigh_delay;
neighbor->dist_check = neigh_dist_check;
// stats for Finish to print
efinal = ecurrent;
f = NULL;
if (ndof) f = atom->f[0];
tmp = 0.0;
for (i = 0; i < ndof; i++) tmp += f[i]*f[i];
MPI_Allreduce(&tmp,&fnorm2_final,1,MPI_DOUBLE,MPI_SUM,world);
if (nextra)
for (i = 0; i < nextra; i++)
fnorm2_final += fextra[i]*fextra[i];
fnorm2_final = sqrt(fnorm2_final);
tmp = 0.0;
for (i = 0; i < ndof; i++) tmp = MAX(fabs(f[i]),tmp);
MPI_Allreduce(&tmp,&fnorminf_final,1,MPI_DOUBLE,MPI_MAX,world);
if (nextra)
for (i = 0; i < nextra; i++)
fnorminf_final = MAX(fabs(fextra[i]),fnorminf_final);
}
/* ----------------------------------------------------------------------
evaluate potential energy and forces
may migrate atoms
if resetflag = 1, update x0 by PBC for atoms that migrate
new energy stored in ecurrent and returned (in case caller not in class)
negative gradient will be stored in atom->f
------------------------------------------------------------------------- */
void Min::eng_force(int *pndof, double **px, double **ph, double **px0,
double *peng, int resetflag)
{
// check for reneighboring
// always communicate since minimizer moved atoms
// if reneighbor, have to setup_vectors() since atoms migrated
int nflag = neighbor->decide();
if (nflag == 0) {
timer->stamp();
comm->communicate();
timer->stamp(TIME_COMM);
} else {
if (triclinic) domain->x2lamda(atom->nlocal);
domain->pbc();
if (domain->box_change) {
domain->reset_box();
comm->setup();
if (neighbor->style) neighbor->setup_bins();
}
timer->stamp();
comm->exchange();
comm->borders();
if (triclinic) domain->lamda2x(atom->nlocal+atom->nghost);
timer->stamp(TIME_COMM);
neighbor->build();
timer->stamp(TIME_NEIGHBOR);
setup_vectors();
if (resetflag) {
double **x = atom->x;
double **x0 = fix_minimize->x0;
int nlocal = atom->nlocal;
double dx,dy,dz,dx0,dy0,dz0;
for (int i = 0; i < nlocal; i++) {
dx = dx0 = x[i][0] - x0[i][0];
dy = dy0 = x[i][1] - x0[i][1];
dz = dz0 = x[i][2] - x0[i][2];
domain->minimum_image(dx,dy,dz);
if (dx != dx0) x0[i][0] = x[i][0] - dx;
if (dy != dy0) x0[i][1] = x[i][1] - dy;
if (dz != dz0) x0[i][2] = x[i][2] - dz;
}
}
}
ev_set(update->ntimestep);
force_clear();
timer->stamp();
if (force->pair) {
force->pair->compute(eflag,vflag);
timer->stamp(TIME_PAIR);
}
if (atom->molecular) {
if (force->bond) force->bond->compute(eflag,vflag);
if (force->angle) force->angle->compute(eflag,vflag);
if (force->dihedral) force->dihedral->compute(eflag,vflag);
if (force->improper) force->improper->compute(eflag,vflag);
timer->stamp(TIME_BOND);
}
if (force->kspace) {
force->kspace->compute(eflag,vflag);
timer->stamp(TIME_KSPACE);
}
if (force->newton) {
comm->reverse_communicate();
timer->stamp(TIME_COMM);
}
// fixes that affect minimization
if (modify->n_min_post_force) modify->min_post_force(vflag);
// compute potential energy of system
// normalize if thermo PE does
ecurrent = pe_compute->compute_scalar();
if (nextra) ecurrent += modify->min_energy(fextra);
if (output->thermo->normflag) ecurrent /= atom->natoms;
// return updated ptrs to caller since atoms may have migrated
*pndof = ndof;
if (ndof) *px = atom->x[0];
else *px = NULL;
*ph = h;
*px0 = x0;
*peng = ecurrent;
}
/* ----------------------------------------------------------------------
set ndof and vector pointers after atoms have migrated
------------------------------------------------------------------------- */
void Min::setup_vectors()
{
ndof = 3 * atom->nlocal;
if (ndof) g = fix_minimize->gradient[0];
else g = NULL;
if (ndof) h = fix_minimize->searchdir[0];
else h = NULL;
if (ndof) x0 = fix_minimize->x0[0];
else x0 = NULL;
}
/* ----------------------------------------------------------------------
clear force on own & ghost atoms
setup and clear other arrays as needed
------------------------------------------------------------------------- */
void Min::force_clear()
{
// clear global force array
// nall includes ghosts only if either newton flag is set
int nall;
if (force->newton) nall = atom->nlocal + atom->nghost;
else nall = atom->nlocal;
double **f = atom->f;
for (int i = 0; i < nall; i++) {
f[i][0] = 0.0;
f[i][1] = 0.0;
f[i][2] = 0.0;
}
if (torqueflag) {
double **torque = atom->torque;
for (int i = 0; i < nall; i++) {
torque[i][0] = 0.0;
torque[i][1] = 0.0;
torque[i][2] = 0.0;
}
}
}
/* ----------------------------------------------------------------------
line minimization methods
find minimum-energy starting at x along dir direction
input: n = # of degrees of freedom on this proc
x = ptr to atom->x[0] as vector
dir = search direction as vector
x0 = ptr to fix->x0[0] as vector, for storing initial coords
eoriginal = energy at initial x
maxdist = max distance to move any atom coord
output: return 0 if successful move, non-zero alpha
return non-zero if failed
alpha = distance moved along dir to set x to minimun eng config
caller has several quantities set via last call to eng_force()
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
atom->f = force (-Grad) is evaulated at new configuration
ecurrent = energy of new configuration
NOTE: when call eng_force: n,x,dir,x0,eng may change due to atom migration
updated values are returned by eng_force()
b/c of migration, linemin routines CANNOT store atom-based quantities
------------------------------------------------------------------------- */
/* ----------------------------------------------------------------------
linemin: backtracking line search (Proc 3.1, p 41 in Nocedal and Wright)
uses no gradient info, but should be very robust
start at maxdist, backtrack until energy decrease is sufficient
------------------------------------------------------------------------- */
int Min::linemin_backtrack(int n, double *x, double *dir,
double *x0, double eoriginal, double maxdist,
double &alpha, int &nfunc)
{
int i,m;
double fdotdirall,fdotdirme,hmax,hme,alpha_extra;
double eng,de_ideal,de;
double *f = NULL;
if (n) f = atom->f[0];
// fdotdirall = projection of search dir along downhill gradient
// if search direction is not downhill, exit with error
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 (nextra)
for (i = 0; i < nextra; i++) fdotdirall += fextra[i]*hextra[i];
if (output->thermo->normflag) fdotdirall /= atom->natoms;
if (fdotdirall <= 0.0) return DOWNHILL;
// initial alpha = stepsize to change any atom coord by maxdist
// alpha <= ALPHA_MAX, else 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 < n; i++) hme = MAX(hme,fabs(dir[i]));
MPI_Allreduce(&hme,&hmax,1,MPI_DOUBLE,MPI_MAX,world);
alpha = MIN(ALPHA_MAX,maxdist/hmax);
if (nextra) {
double alpha_extra = modify->max_alpha(hextra);
alpha = MIN(alpha,alpha_extra);
for (i = 0; i < nextra; i++)
hmax = MAX(hmax,fabs(hextra[i]));
}
if (hmax == 0.0) return ZEROFORCE;
// store coords and other dof at start of linesearch
for (i = 0; i < n; i++) x0[i] = x[i];
if (nextra) modify->min_store();
// backtrack with alpha until energy decrease is sufficient
while (1) {
if (nextra) modify->min_step(0.0,hextra);
for (i = 0; i < n; i++) x[i] = x0[i];
if (nextra) modify->min_step(alpha,hextra);
for (i = 0; i < n; i++) x[i] += alpha*dir[i];
eng_force(&n,&x,&dir,&x0,&eng,1);
nfunc++;
// if energy change is better than ideal, exit with success
de_ideal = -BACKTRACK_SLOPE*alpha*fdotdirall;
de = eng - eoriginal;
if (de <= de_ideal) return 0;
// reduce alpha
alpha *= ALPHA_REDUCE;
// backtracked all the way to 0.0
// reset to starting point, exit with error
if (alpha <= 0.0 || de_ideal >= -IDEAL_TOL) {
if (nextra) modify->min_step(0.0,hextra);
for (i = 0; i < n; i++) x[i] = x0[i];
eng_force(&n,&x,&dir,&x0,&eng,0);
nfunc++;
return ZEROALPHA;
}
}
}
/* ----------------------------------------------------------------------
linemin: quadratic line search (adapted from Dennis and Schnabel)
basic idea is to backtrack until change in energy is sufficiently small
based on ENERGY_QUADRATIC, then use a quadratic approximation
using forces at two alpha values to project to minimum
use forces rather than energy change to do projection
this is b/c the forces are going to zero and can become very small
unlike energy differences which are the difference of two finite
values and are thus limited by machine precision
two changes that were critical to making this method work:
a) limit maximum step to alpha <= 1
b) ignore energy criterion if delE <= ENERGY_QUADRATIC
several other ideas also seemed to help:
c) making each step from starting point (alpha = 0), not previous alpha
d) quadratic model based on forces, not energy
e) exiting immediately if f.dir <= 0 (search direction not downhill)
so that CG can restart
a,c,e were also adopted for the backtracking linemin function
------------------------------------------------------------------------- */
int Min::linemin_quadratic(int n, double *x, double *dir,
double *x0, double eoriginal, double maxdist,
double &alpha, int &nfunc)
{
int i,m;
double fdotdirall,fdotdirme,hmax,hme,alphamax,alpha_extra;
double eng,de_ideal,de;
double delfh,engprev,relerr,alphaprev,fhprev,ff,fh,alpha0,fh0,ff0;
double dot[2],dotall[2];
double *f = atom->f[0];
// fdotdirall = projection of search dir along downhill gradient
// if search direction is not downhill, exit with error
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 (nextra)
for (i = 0; i < nextra; i++) fdotdirall += fextra[i]*hextra[i];
if (output->thermo->normflag) fdotdirall /= atom->natoms;
if (fdotdirall <= 0.0) return DOWNHILL;
// initial alpha = stepsize to change any atom coord by maxdist
// alpha <= ALPHA_MAX, else 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 < n; i++) hme = MAX(hme,fabs(dir[i]));
MPI_Allreduce(&hme,&hmax,1,MPI_DOUBLE,MPI_MAX,world);
alpha = MIN(ALPHA_MAX,maxdist/hmax);
if (nextra) {
double alpha_extra = modify->max_alpha(hextra);
alpha = MIN(alpha,alpha_extra);
for (i = 0; i < nextra; i++)
hmax = MAX(hmax,fabs(hextra[i]));
}
if (hmax == 0.0) return ZEROFORCE;
// store coords and other dof at start of linesearch
for (i = 0; i < n; i++) x0[i] = x[i];
if (nextra) modify->min_store();
// backtrack with alpha until energy decrease is sufficient
// or until get to small energy change, then perform quadratic projection
fhprev = fdotdirall;
engprev = eoriginal;
alphaprev = 0.0;
while (1) {
if (nextra) modify->min_step(0.0,hextra);
for (i = 0; i < n; i++) x[i] = x0[i];
if (nextra) modify->min_step(alpha,hextra);
for (i = 0; i < n; i++) x[i] += alpha*dir[i];
eng_force(&n,&x,&dir,&x0,&eng,1);
nfunc++;
// compute new fh, alpha, delfh
dot[0] = dot[1] = 0.0;
for (i = 0; i < ndof; i++) {
dot[0] += f[i]*f[i];
dot[1] += f[i]*dir[i];
}
MPI_Allreduce(dot,dotall,2,MPI_DOUBLE,MPI_SUM,world);
if (nextra) {
for (i = 0; i < nextra; 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;
}
delfh = fh - fhprev;
// if fh or delfh is epsilon, reset to starting point, exit with error
if (fabs(fh) < EPS_QUAD || fabs(delfh) < EPS_QUAD) {
if (nextra) modify->min_step(0.0,hextra);
for (i = 0; i < n; i++) x[i] = x0[i];
eng_force(&n,&x,&dir,&x0,&eng,0);
nfunc++;
return ZEROQUAD;
}
// check if ready for quadratic projection, equivalent to secant method
// alpha0 = projected alpha
relerr = fabs(1.0+(0.5*alpha*(alpha-alphaprev)*(fh+fhprev)-eng)/engprev);
alpha0 = alpha - (alpha-alphaprev)*fh/delfh;
if (relerr <= QUADRATIC_TOL && alpha0 > 0.0) {
if (nextra) modify->min_step(0.0,hextra);
for (i = 0; i < n; i++) x[i] = x0[i];
if (nextra) modify->min_step(alpha0,hextra);
for (i = 0; i < n; i++) x[i] += alpha0*dir[i];
eng_force(&n,&x,&dir,&x0,&eng,1);
nfunc++;
// if backtracking energy change is better than ideal, exit with success
de_ideal = -BACKTRACK_SLOPE*alpha0*fdotdirall;
de = eng - eoriginal;
if (de <= de_ideal || de_ideal >= -IDEAL_TOL) return 0;
// drop back from alpha0 to alpha
if (nextra) modify->min_step(0.0,hextra);
for (i = 0; i < n; i++) x[i] = x0[i];
if (nextra) modify->min_step(alpha,hextra);
for (i = 0; i < n; i++) x[i] += alpha*dir[i];
eng_force(&n,&x,&dir,&x0,&eng,1);
nfunc++;
}
// if backtracking energy change is better than ideal, exit with success
de_ideal = -BACKTRACK_SLOPE*alpha*fdotdirall;
de = eng - eoriginal;
if (de <= de_ideal) return 0;
// save previous state
fhprev = fh;
engprev = eng;
alphaprev = alpha;
// reduce alpha
alpha *= ALPHA_REDUCE;
// backtracked all the way to 0.0
// reset to starting point, exit with error
if (alpha <= 0.0 || de_ideal >= -IDEAL_TOL) {
if (nextra) modify->min_step(0.0,hextra);
for (i = 0; i < n; i++) x[i] = x0[i];
eng_force(&n,&x,&dir,&x0,&eng,0);
nfunc++;
return ZEROALPHA;
}
}
}
/* ---------------------------------------------------------------------- */

43
src/min.h
View File

@ -28,13 +28,13 @@ class Min : protected Pointers {
Min(class LAMMPS *);
virtual ~Min();
virtual void init() = 0;
virtual void run() = 0;
virtual int iterate(int) = 0;
virtual double memory_usage() {return 0.0;}
void init();
void run();
double memory_usage() {return 0.0;}
void modify_params(int, char **);
virtual int iterate(int) = 0;
protected:
int eflag,vflag; // flags for energy/virial computation
int virial_style; // compute virial explicitly or implicitly
@ -48,6 +48,39 @@ class Min : protected Pointers {
class Compute **vlist_global;
class Compute **vlist_atom;
int pairflag,torqueflag;
int neigh_every,neigh_delay,neigh_dist_check; // copies of reneigh criteria
int triclinic; // 0 if domain is orthog, 1 if triclinic
class FixMinimize *fix_minimize; // fix that stores gradient vecs
class Compute *pe_compute; // compute for potential energy
double ecurrent; // current potential energy
double mindist,maxdist; // min/max dist for coord delta in line search
int ndof; // # of degrees-of-freedom on this proc
double *g,*h; // local portion of gradient, searchdir vectors
double *x0; // coords at start of linesearch
int nextra; // extra dof due to fixes
double *fextra; // vectors for extra dof
double *gextra;
double *hextra;
// ptr to linemin functions
void setup();
void eng_force(int *, double **, double **, double **, double *, int);
void setup_vectors();
void force_clear();
typedef int (Min::*FnPtr)(int, double *, double *, double *, double,
double, double &, int &);
FnPtr linemin;
int linemin_backtrack(int, double *, double *, double *, double,
double, double &, int &);
int linemin_quadratic(int, double *, double *, double *, double,
double, double &, int &);
void ev_setup();
void ev_set(int);
};

729
src/min_cg.cpp
View File

@ -11,325 +11,31 @@
See the README file in the top-level LAMMPS directory.
------------------------------------------------------------------------- */
/* ----------------------------------------------------------------------
Contributing author: Aidan Thompson (SNL)
improved CG and backtrack ls, added quadratic 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
------------------------------------------------------------------------- */
#include "math.h"
#include "string.h"
#include "mpi.h"
#include "min_cg.h"
#include "neighbor.h"
#include "domain.h"
#include "comm.h"
#include "atom.h"
#include "force.h"
#include "pair.h"
#include "bond.h"
#include "angle.h"
#include "dihedral.h"
#include "improper.h"
#include "kspace.h"
#include "output.h"
#include "thermo.h"
#include "update.h"
#include "modify.h"
#include "compute.h"
#include "fix_minimize.h"
#include "thermo.h"
#include "output.h"
#include "timer.h"
#include "memory.h"
#include "error.h"
using namespace LAMMPS_NS;
#define MIN(A,B) ((A) < (B)) ? (A) : (B)
#define MAX(A,B) ((A) > (B)) ? (A) : (B)
// ALPHA_MAX = max alpha allowed to avoid long backtracks
// ALPHA_REDUCE = reduction ratio, should be in range [0.5,1)
// BACKTRACK_SLOPE, should be in range (0,0.5]
// QUADRATIC_TOL = tolerance on alpha0, should be in range [0.1,1)
// EPS_ENERGY = minimum normalization for energy tolerance
// IDEAL_TOL = ideal energy tolerance for backtracking
// EPS_QUAD = tolerance for quadratic projection
#define ALPHA_MAX 1.0
#define ALPHA_REDUCE 0.5
#define BACKTRACK_SLOPE 0.4
#define QUADRATIC_TOL 0.1
#define EPS_ENERGY 1.0e-8
#define IDEAL_TOL 1.0e-8
#define EPS_QUAD 1.0e-28
// same as in other min classes
enum{MAXITER,MAXEVAL,ETOL,FTOL,DOWNHILL,ZEROALPHA,ZEROFORCE,ZEROQUAD};
/* ---------------------------------------------------------------------- */
MinCG::MinCG(LAMMPS *lmp) : Min(lmp)
{
fextra = gextra = hextra = NULL;
}
#define MIN(A,B) ((A) < (B)) ? (A) : (B)
#define MAX(A,B) ((A) > (B)) ? (A) : (B)
/* ---------------------------------------------------------------------- */
MinCG::~MinCG()
{
delete [] fextra;
delete [] gextra;
delete [] hextra;
}
/* ---------------------------------------------------------------------- */
void MinCG::init()
{
// create fix needed for storing atom-based gradient vectors
// will delete it at end of run
char **fixarg = new char*[3];
fixarg[0] = (char *) "MINIMIZE";
fixarg[1] = (char *) "all";
fixarg[2] = (char *) "MINIMIZE";
modify->add_fix(3,fixarg);
delete [] fixarg;
fix_minimize = (FixMinimize *) modify->fix[modify->nfix-1];
// zero gradient vectors before first atom exchange
setup_vectors();
for (int i = 0; i < ndof; i++) h[i] = g[i] = 0.0;
// virial_style:
// 1 if computed explicitly by pair->compute via sum over pair interactions
// 2 if computed implicitly by pair->virial_compute via sum over ghost atoms
if (force->newton_pair) virial_style = 2;
else virial_style = 1;
// setup lists of computes for global and per-atom PE and pressure
ev_setup();
// set flags for what arrays to clear in force_clear()
// clear torques if array exists
torqueflag = 0;
if (atom->torque) torqueflag = 1;
// orthogonal vs triclinic simulation box
triclinic = domain->triclinic;
// reset reneighboring criteria if necessary
neigh_every = neighbor->every;
neigh_delay = neighbor->delay;
neigh_dist_check = neighbor->dist_check;
if (neigh_every != 1 || neigh_delay != 0 || neigh_dist_check != 1) {
if (comm->me == 0)
error->warning("Resetting reneighboring criteria during minimization");
}
neighbor->every = 1;
neighbor->delay = 0;
neighbor->dist_check = 1;
// set ptr to linemin function
if (linestyle == 0) linemin = &MinCG::linemin_backtrack;
else if (linestyle == 1) linemin = &MinCG::linemin_quadratic;
}
/* ----------------------------------------------------------------------
perform minimization, with setup first
------------------------------------------------------------------------- */
void MinCG::run()
{
int i;
double tmp,*f;
// 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"};
// set initial force & energy
// normalize energy if thermo PE does
setup();
// setup any extra dof due to fixes
// can't be done until now b/c update init() comes before modify init()
delete [] fextra;
delete [] gextra;
delete [] hextra;
nextra = modify->min_dof();
if (nextra) {
fextra = new double[nextra];
gextra = new double[nextra];
hextra = new double[nextra];
}
// compute potential energy of system
// normalize if thermo PE does
int id = modify->find_compute("thermo_pe");
if (id < 0) error->all("Minimization could not find thermo_pe compute");
pe_compute = modify->compute[id];
ecurrent = pe_compute->compute_scalar();
if (nextra) ecurrent += modify->min_energy(fextra);
if (output->thermo->normflag) ecurrent /= atom->natoms;
// stats for Finish to print
einitial = ecurrent;
f = NULL;
if (ndof) f = atom->f[0];
tmp = 0.0;
for (i = 0; i < ndof; i++) tmp += f[i]*f[i];
MPI_Allreduce(&tmp,&fnorm2_init,1,MPI_DOUBLE,MPI_SUM,world);
if (nextra)
for (i = 0; i < nextra; i++) fnorm2_init += fextra[i]*fextra[i];
fnorm2_init = sqrt(fnorm2_init);
tmp = 0.0;
for (i = 0; i < ndof; i++) tmp = MAX(fabs(f[i]),tmp);
MPI_Allreduce(&tmp,&fnorminf_init,1,MPI_DOUBLE,MPI_MAX,world);
if (nextra)
for (i = 0; i < nextra; i++)
fnorminf_init = MAX(fabs(fextra[i]),fnorminf_init);
// minimizer iterations
timer->barrier_start(TIME_LOOP);
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
// set output->next values to this timestep
// call eng_force to insure vflag is set when forces computed
// output->write does final output for thermo, dump, restart files
// add ntimestep to all computes that store invocation times
// since are hardwireing call to thermo/dumps and computes may not be ready
if (niter < update->nsteps) {
niter++;
update->nsteps = niter;
for (int idump = 0; idump < output->ndump; idump++)
output->next_dump[idump] = update->ntimestep;
output->next_dump_any = update->ntimestep;
if (output->restart_every) output->next_restart = update->ntimestep;
output->next_thermo = update->ntimestep;
modify->addstep_compute_all(update->ntimestep);
int ntmp;
double *xtmp,*htmp,*x0tmp,etmp;
eng_force(&ntmp,&xtmp,&htmp,&x0tmp,&etmp);
output->write(update->ntimestep);
}
timer->barrier_stop(TIME_LOOP);
// delete fix at end of run, so its atom arrays won't persist
modify->delete_fix("MINIMIZE");
// reset reneighboring criteria
neighbor->every = neigh_every;
neighbor->delay = neigh_delay;
neighbor->dist_check = neigh_dist_check;
// stats for Finish to print
efinal = ecurrent;
f = NULL;
if (ndof) f = atom->f[0];
tmp = 0.0;
for (i = 0; i < ndof; i++) tmp += f[i]*f[i];
MPI_Allreduce(&tmp,&fnorm2_final,1,MPI_DOUBLE,MPI_SUM,world);
if (nextra)
for (i = 0; i < nextra; i++)
fnorm2_final += fextra[i]*fextra[i];
fnorm2_final = sqrt(fnorm2_final);
tmp = 0.0;
for (i = 0; i < ndof; i++) tmp = MAX(fabs(f[i]),tmp);
MPI_Allreduce(&tmp,&fnorminf_final,1,MPI_DOUBLE,MPI_MAX,world);
if (nextra)
for (i = 0; i < nextra; i++)
fnorminf_final = MAX(fabs(fextra[i]),fnorminf_final);
}
/* ----------------------------------------------------------------------
setup before run
------------------------------------------------------------------------- */
void MinCG::setup()
{
if (comm->me == 0 && screen) fprintf(screen,"Setting up minimization ...\n");
// setup domain, communication and neighboring
// acquire ghosts
// build neighbor lists
// reset gradient vector ptrs
if (triclinic) domain->x2lamda(atom->nlocal);
domain->pbc();
domain->reset_box();
comm->setup();
if (neighbor->style) neighbor->setup_bins();
comm->exchange();
comm->borders();
if (triclinic) domain->lamda2x(atom->nlocal+atom->nghost);
neighbor->build();
neighbor->ncalls = 0;
setup_vectors();
// compute all forces
ev_set(update->ntimestep);
force_clear();
if (force->pair) force->pair->compute(eflag,vflag);
if (atom->molecular) {
if (force->bond) force->bond->compute(eflag,vflag);
if (force->angle) force->angle->compute(eflag,vflag);
if (force->dihedral) force->dihedral->compute(eflag,vflag);
if (force->improper) force->improper->compute(eflag,vflag);
}
if (force->kspace) {
force->kspace->setup();
force->kspace->compute(eflag,vflag);
}
if (force->newton) comm->reverse_communicate();
modify->setup(vflag);
output->setup(1);
}
MinCG::MinCG(LAMMPS *lmp) : Min(lmp) {}
/* ----------------------------------------------------------------------
minimization via conjugate gradient iterations
@ -343,6 +49,10 @@ int MinCG::iterate(int n)
double *x = NULL;
double *f = NULL;
int ndoftotal;
MPI_Allreduce(&ndof,&ndoftotal,1,MPI_INT,MPI_SUM,world);
ndoftotal += nextra;
if (ndof) f = atom->f[0];
for (i = 0; i < ndof; i++) h[i] = g[i] = f[i];
if (nextra)
@ -355,7 +65,6 @@ int MinCG::iterate(int n)
if (nextra)
for (i = 0; i < nextra; i++) gg += fextra[i]*fextra[i];
int ndoftotal = ndof + nextra;
neval = 0;
for (niter = 0; niter < n; niter++) {
@ -420,6 +129,7 @@ int MinCG::iterate(int n)
dot[0] = 0.0;
for (i = 0; i < ndof; i++) dot[0] += g[i]*h[i];
MPI_Allreduce(dot,dotall,1,MPI_DOUBLE,MPI_SUM,world);
if (nextra)
for (i = 0; i < nextra; i++)
dotall[0] += gextra[i]*hextra[i];
@ -443,424 +153,3 @@ int MinCG::iterate(int n)
return MAXITER;
}
/* ----------------------------------------------------------------------
set ndof and vector pointers after atoms have migrated
------------------------------------------------------------------------- */
void MinCG::setup_vectors()
{
ndof = 3 * atom->nlocal;
if (ndof) g = fix_minimize->gradient[0];
else g = NULL;
if (ndof) h = fix_minimize->searchdir[0];
else h = NULL;
if (ndof) x0 = fix_minimize->x0[0];
else x0 = NULL;
}
/* ----------------------------------------------------------------------
evaluate potential energy and forces
may migrate atoms
new energy stored in ecurrent and returned (in case caller not in class)
negative gradient will be stored in atom->f
------------------------------------------------------------------------- */
void MinCG::eng_force(int *pndof, double **px, double **ph, double **px0,
double *peng)
{
// check for reneighboring
// always communicate since minimizer moved atoms
// if reneighbor, have to setup_vectors() since atoms migrated
int nflag = neighbor->decide();
if (nflag == 0) {
timer->stamp();
comm->communicate();
timer->stamp(TIME_COMM);
} else {
if (triclinic) domain->x2lamda(atom->nlocal);
domain->pbc();
if (domain->box_change) {
domain->reset_box();
comm->setup();
if (neighbor->style) neighbor->setup_bins();
}
timer->stamp();
comm->exchange();
comm->borders();
if (triclinic) domain->lamda2x(atom->nlocal+atom->nghost);
timer->stamp(TIME_COMM);
neighbor->build();
timer->stamp(TIME_NEIGHBOR);
setup_vectors();
}
ev_set(update->ntimestep);
force_clear();
timer->stamp();
if (force->pair) {
force->pair->compute(eflag,vflag);
timer->stamp(TIME_PAIR);
}
if (atom->molecular) {
if (force->bond) force->bond->compute(eflag,vflag);
if (force->angle) force->angle->compute(eflag,vflag);
if (force->dihedral) force->dihedral->compute(eflag,vflag);
if (force->improper) force->improper->compute(eflag,vflag);
timer->stamp(TIME_BOND);
}
if (force->kspace) {
force->kspace->compute(eflag,vflag);
timer->stamp(TIME_KSPACE);
}
if (force->newton) {
comm->reverse_communicate();
timer->stamp(TIME_COMM);
}
// fixes that affect minimization
if (modify->n_min_post_force) modify->min_post_force(vflag);
// compute potential energy of system
// normalize if thermo PE does
ecurrent = pe_compute->compute_scalar();
if (nextra) ecurrent += modify->min_energy(fextra);
if (output->thermo->normflag) ecurrent /= atom->natoms;
// return updated ptrs to caller since atoms may have migrated
*pndof = ndof;
if (ndof) *px = atom->x[0];
else *px = NULL;
*ph = h;
*px0 = x0;
*peng = ecurrent;
}
/* ----------------------------------------------------------------------
clear force on own & ghost atoms
setup and clear other arrays as needed
------------------------------------------------------------------------- */
void MinCG::force_clear()
{
// clear global force array
// nall includes ghosts only if either newton flag is set
int nall;
if (force->newton) nall = atom->nlocal + atom->nghost;
else nall = atom->nlocal;
double **f = atom->f;
for (int i = 0; i < nall; i++) {
f[i][0] = 0.0;
f[i][1] = 0.0;
f[i][2] = 0.0;
}
if (torqueflag) {
double **torque = atom->torque;
for (int i = 0; i < nall; i++) {
torque[i][0] = 0.0;
torque[i][1] = 0.0;
torque[i][2] = 0.0;
}
}
}
/* ----------------------------------------------------------------------
line minimization methods
find minimum-energy starting at x along dir direction
input: n = # of degrees of freedom on this proc
x = ptr to atom->x[0] as vector
dir = search direction as vector
x0 = ptr to fix->x0[0] as vector, for storing initial coords
eng = current energy at initial x
maxdist = max distance to move any atom coord
output: return 0 if successful move, non-zero alpha
return non-zero if failed
alpha = distance moved along dir to set x to minimun eng config
caller has several quantities set via last call to eng_force()
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
atom->f = force (-Grad) is evaulated at new configuration
ecurrent = energy of new configuration
NOTE: when call eng_force: n,x,dir,x0,eng may change due to atom migration
updated values are returned by eng_force()
b/c of migration, linemin routines CANNOT store atom-based quantities
------------------------------------------------------------------------- */
/* ----------------------------------------------------------------------
linemin: backtracking line search (Proc 3.1, p 41 in Nocedal and Wright)
uses no gradient info, but should be very robust
start at maxdist, backtrack until energy decrease is sufficient
------------------------------------------------------------------------- */
int MinCG::linemin_backtrack(int n, double *x, double *dir,
double *x0, double eng, double maxdist,
double &alpha, int &nfunc)
{
int i,m;
double fdotdirall,fdotdirme,hmax,hme,alpha_extra,eoriginal;
double de_ideal,de;
double *f = NULL;
if (n) f = atom->f[0];
// fdotdirall = projection of search dir along downhill gradient
// if search direction is not downhill, exit with error
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 (nextra)
for (i = 0; i < nextra; i++) fdotdirall += fextra[i]*hextra[i];
if (output->thermo->normflag) fdotdirall /= atom->natoms;
if (fdotdirall <= 0.0) return DOWNHILL;
// initial alpha = stepsize to change any atom coord by maxdist
// alpha <= ALPHA_MAX, else 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 < n; i++) hme = MAX(hme,fabs(dir[i]));
MPI_Allreduce(&hme,&hmax,1,MPI_DOUBLE,MPI_MAX,world);
alpha = MIN(ALPHA_MAX,maxdist/hmax);
if (nextra) {
double alpha_extra = modify->max_alpha(hextra);
alpha = MIN(alpha,alpha_extra);
for (i = 0; i < nextra; i++)
hmax = MAX(hmax,fabs(hextra[i]));
}
if (hmax == 0.0) return ZEROFORCE;
// store coords and other dof at start of linesearch
for (i = 0; i < n; i++) x0[i] = x[i];
if (nextra) modify->min_store();
// eoriginal = energy at start of linesearch
eng_force(&n,&x,&dir,&x0,&eng);
nfunc++;
eoriginal = eng;
// backtrack with alpha until energy decrease is sufficient
while (1) {
if (nextra) modify->min_step(0.0,hextra);
for (i = 0; i < n; i++) x[i] = x0[i];
if (nextra) modify->min_step(alpha,hextra);
for (i = 0; i < n; i++) x[i] += alpha*dir[i];
eng_force(&n,&x,&dir,&x0,&eng);
nfunc++;
// if energy change is better than ideal, exit with success
de_ideal = -BACKTRACK_SLOPE*alpha*fdotdirall;
de = eng - eoriginal;
if (de <= de_ideal) return 0;
// reduce alpha
alpha *= ALPHA_REDUCE;
// backtracked all the way to 0.0
// reset to starting point, exit with error
if (alpha <= 0.0 || de_ideal >= -IDEAL_TOL) {
if (nextra) modify->min_step(0.0,hextra);
for (i = 0; i < n; i++) x[i] = x0[i];
eng_force(&n,&x,&dir,&x0,&eng);
nfunc++;
return ZEROALPHA;
}
}
}
/* ----------------------------------------------------------------------
linemin: quadratic line search (adapted from Dennis and Schnabel)
basic idea is to backtrack until change in energy is sufficiently small
based on ENERGY_QUADRATIC, then use a quadratic approximation
using forces at two alpha values to project to minimum
use forces rather than energy change to do projection
this is b/c the forces are going to zero and can become very small
unlike energy differences which are the difference of two finite
values and are thus limited by machine precision
two changes that were critical to making this method work:
a) limit maximum step to alpha <= 1
b) ignore energy criterion if delE <= ENERGY_QUADRATIC
several other ideas also seemed to help:
c) making each step from starting point (alpha = 0), not previous alpha
d) quadratic model based on forces, not energy
e) exiting immediately if f.dir <= 0 (search direction not downhill)
so that CG can restart
a,c,e were also adopted for the backtracking linemin function
------------------------------------------------------------------------- */
int MinCG::linemin_quadratic(int n, double *x, double *dir,
double *x0, double eng, double maxdist,
double &alpha, int &nfunc)
{
int i,m;
double fdotdirall,fdotdirme,hmax,hme,alphamax,alpha_extra,eoriginal;
double de_ideal,de;
double delfh,engprev,relerr,alphaprev,fhprev,ff,fh,alpha0,fh0,ff0;
double dot[2],dotall[2];
double *f = atom->f[0];
// fdotdirall = projection of search dir along downhill gradient
// if search direction is not downhill, exit with error
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 (nextra)
for (i = 0; i < nextra; i++) fdotdirall += fextra[i]*hextra[i];
if (output->thermo->normflag) fdotdirall /= atom->natoms;
if (fdotdirall <= 0.0) return DOWNHILL;
// initial alpha = stepsize to change any atom coord by maxdist
// alpha <= ALPHA_MAX, else 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 < n; i++) hme = MAX(hme,fabs(dir[i]));
MPI_Allreduce(&hme,&hmax,1,MPI_DOUBLE,MPI_MAX,world);
alpha = MIN(ALPHA_MAX,maxdist/hmax);
if (nextra) {
double alpha_extra = modify->max_alpha(hextra);
alpha = MIN(alpha,alpha_extra);
for (i = 0; i < nextra; i++)
hmax = MAX(hmax,fabs(hextra[i]));
}
if (hmax == 0.0) return ZEROFORCE;
// store coords and other dof at start of linesearch
for (i = 0; i < n; i++) x0[i] = x[i];
if (nextra) modify->min_store();
// eoriginal = energy at start of linesearch
eng_force(&n,&x,&dir,&x0,&eng);
nfunc++;
eoriginal = eng;
// backtrack with alpha until energy decrease is sufficient
// or until get to small energy change, then perform quadratic projection
fhprev = fdotdirall;
engprev = eoriginal;
alphaprev = 0.0;
while (1) {
if (nextra) modify->min_step(0.0,hextra);
for (i = 0; i < n; i++) x[i] = x0[i];
if (nextra) modify->min_step(alpha,hextra);
for (i = 0; i < n; i++) x[i] += alpha*dir[i];
eng_force(&n,&x,&dir,&x0,&eng);
nfunc++;
// compute new fh, alpha, delfh
dot[0] = dot[1] = 0.0;
for (i = 0; i < ndof; i++) {
dot[0] += f[i]*f[i];
dot[1] += f[i]*dir[i];
}
MPI_Allreduce(dot,dotall,2,MPI_DOUBLE,MPI_SUM,world);
if (nextra) {
for (i = 0; i < nextra; 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;
}
delfh = fh - fhprev;
// if fh or delfh is epsilon, reset to starting point, exit with error
if (fabs(fh) < EPS_QUAD || fabs(delfh) < EPS_QUAD) {
if (nextra) modify->min_step(0.0,hextra);
for (i = 0; i < n; i++) x[i] = x0[i];
eng_force(&n,&x,&dir,&x0,&eng);
nfunc++;
return ZEROQUAD;
}
// check if ready for quadratic projection, equivalent to secant method
// alpha0 = projected alpha
relerr = fabs(1.0+(0.5*alpha*(alpha-alphaprev)*(fh+fhprev)-eng)/engprev);
alpha0 = alpha - (alpha-alphaprev)*fh/delfh;
if (relerr <= QUADRATIC_TOL && alpha0 > 0.0) {
if (nextra) modify->min_step(0.0,hextra);
for (i = 0; i < n; i++) x[i] = x0[i];
if (nextra) modify->min_step(alpha0,hextra);
for (i = 0; i < n; i++) x[i] += alpha0*dir[i];
eng_force(&n,&x,&dir,&x0,&eng);
nfunc++;
// if backtracking energy change is better than ideal, exit with success
de_ideal = -BACKTRACK_SLOPE*alpha0*fdotdirall;
de = eng - eoriginal;
if (de <= de_ideal || de_ideal >= -IDEAL_TOL) return 0;
// drop back from alpha0 to alpha
if (nextra) modify->min_step(0.0,hextra);
for (i = 0; i < n; i++) x[i] = x0[i];
if (nextra) modify->min_step(alpha,hextra);
for (i = 0; i < n; i++) x[i] += alpha*dir[i];
eng_force(&n,&x,&dir,&x0,&eng);
nfunc++;
}
// if backtracking energy change is better than ideal, exit with success
de_ideal = -BACKTRACK_SLOPE*alpha*fdotdirall;
de = eng - eoriginal;
if (de <= de_ideal) return 0;
// save previous state
fhprev = fh;
engprev = eng;
alphaprev = alpha;
// reduce alpha
alpha *= ALPHA_REDUCE;
// backtracked all the way to 0.0
// reset to starting point, exit with error
if (alpha <= 0.0 || de_ideal >= -IDEAL_TOL) {
if (nextra) modify->min_step(0.0,hextra);
for (i = 0; i < n; i++) x[i] = x0[i];
eng_force(&n,&x,&dir,&x0,&eng);
nfunc++;
return ZEROALPHA;
}
}
}

39
src/min_cg.h
View File

@ -21,44 +21,7 @@ namespace LAMMPS_NS {
class MinCG : public Min {
public:
MinCG(class LAMMPS *);
virtual ~MinCG();
void init();
void run();
virtual int iterate(int);
protected:
int pairflag,torqueflag;
int neigh_every,neigh_delay,neigh_dist_check; // copies of reneigh criteria
int triclinic; // 0 if domain is orthog, 1 if triclinic
class FixMinimize *fix_minimize; // fix that stores gradient vecs
class Compute *pe_compute; // compute for potential energy
double ecurrent; // current potential energy
double mindist,maxdist; // min/max dist for coord delta in line search
int ndof; // # of degrees-of-freedom on this proc
double *g,*h; // local portion of gradient, searchdir vectors
double *x0; // coords at start of linesearch
int nextra; // extra dof due to fixes
double *fextra; // vectors for extra dof
double *gextra;
double *hextra;
// ptr to linemin functions
typedef int (MinCG::*FnPtr)(int, double *, double *, double *, double,
double, double &, int &);
FnPtr linemin;
int linemin_backtrack(int, double *, double *, double *, double,
double, double &, int &);
int linemin_quadratic(int, double *, double *, double *, double,
double, double &, int &);
void setup();
void setup_vectors();
void eng_force(int *, double **, double **, double **, double *);
void force_clear();
int iterate(int);
};
}

4
src/min_sd.cpp
View File

@ -21,6 +21,8 @@
using namespace LAMMPS_NS;
// EPS_ENERGY = minimum normalization for energy tolerance
#define EPS_ENERGY 1.0e-8
// same as in other min classes
@ -29,7 +31,7 @@ enum{MAXITER,MAXEVAL,ETOL,FTOL,DOWNHILL,ZEROALPHA,ZEROFORCE,ZEROQUAD};
/* ---------------------------------------------------------------------- */
MinSD::MinSD(LAMMPS *lmp) : MinCG(lmp) {}
MinSD::MinSD(LAMMPS *lmp) : Min(lmp) {}
/* ----------------------------------------------------------------------
minimization via steepest descent

2
src/min_sd.h
View File

@ -18,7 +18,7 @@
namespace LAMMPS_NS {
class MinSD : public MinCG {
class MinSD : public Min {
public:
MinSD(class LAMMPS *);
int iterate(int);