diff --git a/src/min.cpp b/src/min.cpp index afbb44cb05..027ca3412a 100644 --- a/src/min.cpp +++ b/src/min.cpp @@ -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; + } + } } /* ---------------------------------------------------------------------- */ diff --git a/src/min.h b/src/min.h index e01999f5bc..65008c2652 100644 --- a/src/min.h +++ b/src/min.h @@ -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); }; diff --git a/src/min_cg.cpp b/src/min_cg.cpp index add6e7fc43..53a80d01d0 100644 --- a/src/min_cg.cpp +++ b/src/min_cg.cpp @@ -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; - } - } -} diff --git a/src/min_cg.h b/src/min_cg.h index 95d66f81dc..6a6480bffe 100644 --- a/src/min_cg.h +++ b/src/min_cg.h @@ -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); }; } diff --git a/src/min_sd.cpp b/src/min_sd.cpp index a2f13e7bad..7342f8d99a 100644 --- a/src/min_sd.cpp +++ b/src/min_sd.cpp @@ -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 diff --git a/src/min_sd.h b/src/min_sd.h index be502162d4..fa06f0e94a 100644 --- a/src/min_sd.h +++ b/src/min_sd.h @@ -18,7 +18,7 @@ namespace LAMMPS_NS { -class MinSD : public MinCG { +class MinSD : public Min { public: MinSD(class LAMMPS *); int iterate(int);