diff --git a/src/min.cpp b/src/min.cpp index 4bd825050..4c00c13f8 100644 --- a/src/min.cpp +++ b/src/min.cpp @@ -1,692 +1,692 @@ /* ---------------------------------------------------------------------- LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator http://lammps.sandia.gov, Sandia National Laboratories Steve Plimpton, sjplimp@sandia.gov Copyright (2003) Sandia Corporation. Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains certain rights in this software. This software is distributed under the GNU General Public License. See the README file in the top-level LAMMPS directory. ------------------------------------------------------------------------- */ /* ---------------------------------------------------------------------- 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 "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) /* ---------------------------------------------------------------------- */ Min::Min(LAMMPS *lmp) : Pointers(lmp) { dmax = 0.1; linestyle = 0; elist_atom = NULL; vlist_global = vlist_atom = NULL; nextra_global = 0; fextra = NULL; nextra_atom = 0; xextra_atom = fextra_atom = NULL; extra_peratom = extra_nlen = NULL; extra_max = NULL; requestor = NULL; } /* ---------------------------------------------------------------------- */ Min::~Min() { delete [] elist_atom; delete [] vlist_global; delete [] vlist_atom; delete [] fextra; memory->sfree(xextra_atom); memory->sfree(fextra_atom); memory->sfree(extra_peratom); memory->sfree(extra_nlen); memory->sfree(extra_max); memory->sfree(requestor); } /* ---------------------------------------------------------------------- */ void Min::init() { // create fix needed for storing atom-based quantities // 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]; // clear out extra global and per-atom dof // will receive requests for new per-atom dof during pair init() // can then add vectors to fix_minimize in setup() nextra_global = 0; delete [] fextra; fextra = NULL; nextra_atom = 0; memory->sfree(xextra_atom); memory->sfree(fextra_atom); memory->sfree(extra_peratom); memory->sfree(extra_nlen); memory->sfree(extra_max); memory->sfree(requestor); xextra_atom = fextra_atom = NULL; extra_peratom = extra_nlen = NULL; extra_max = NULL; requestor = NULL; // 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; // style-specific initialization init_style(); } /* ---------------------------------------------------------------------- setup before run ------------------------------------------------------------------------- */ void Min::setup() { if (comm->me == 0 && screen) fprintf(screen,"Setting up minimization ...\n"); // setup extra global dof due to fixes // cannot be done in init() b/c update init() is before modify init() nextra_global = modify->min_dof(); if (nextra_global) fextra = new double[nextra_global]; // compute for potential energy int id = modify->find_compute("thermo_pe"); if (id < 0) error->all("Minimization could not find thermo_pe compute"); pe_compute = modify->compute[id]; // style-specific setup does two tasks // setup extra global dof vectors // setup extra per-atom dof vectors due to requests from Pair classes // cannot be done in init() b/c update init() is before modify/pair init() setup_style(); // ndoftotal = total dof for entire minimization problem // dof for atoms, extra per-atom, extra global double ndofme = 3.0*atom->nlocal; for (int m = 0; m < nextra_atom; m++) ndofme += extra_peratom[m]*atom->nlocal; MPI_Allreduce(&ndofme,&ndoftotal,1,MPI_DOUBLE,MPI_SUM,world); ndoftotal += nextra_global; // setup domain, communication and neighboring // acquire ghosts // build neighbor lists atom->setup(); if (triclinic) domain->x2lamda(atom->nlocal); domain->pbc(); domain->reset_box(); comm->setup(); if (neighbor->style) neighbor->setup_bins(); comm->exchange(); if (atom->sortfreq > 0) atom->sort(); comm->borders(); if (triclinic) domain->lamda2x(atom->nlocal+atom->nghost); neighbor->build(); neighbor->ncalls = 0; // 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_comm(); modify->setup(vflag); output->setup(1); // atoms may have migrated in comm->exchange() reset_vectors(); } /* ---------------------------------------------------------------------- setup without output or one-time post-init setup flag = 0 = just force calculation flag = 1 = reneighbor and force calculation ------------------------------------------------------------------------- */ void Min::setup_minimal(int flag) { // setup domain, communication and neighboring // acquire ghosts // build neighbor lists if (flag) { 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; } // 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_comm(); modify->setup(vflag); // atoms may have migrated in comm->exchange() reset_vectors(); } /* ---------------------------------------------------------------------- perform minimization, calling iterate() for nsteps ------------------------------------------------------------------------- */ void Min::run(int nsteps) { // possible stop conditions char *stopstrings[] = {"max iterations", "max force evaluations", "energy tolerance", "force tolerance", "search direction is not downhill", "linesearch alpha is zero", "forces are zero", "quadratic factors are zero", "trust region too small", "HFTN minimizer error"}; // stats for Finish to print ecurrent = pe_compute->compute_scalar(); if (nextra_global) ecurrent += modify->min_energy(fextra); if (output->thermo->normflag) ecurrent /= atom->natoms; einitial = ecurrent; fnorm2_init = sqrt(fnorm_sqr()); fnorminf_init = fnorm_inf(); // minimizer iterations int stop_condition = iterate(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 engergy_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 < nsteps) { niter++; update->nsteps = niter; if (update->restrict_output == 0) { 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); ecurrent = energy_force(0); output->write(update->ntimestep); } // stats for Finish to print efinal = ecurrent; fnorm2_final = sqrt(fnorm_sqr()); fnorminf_final = fnorm_inf(); } /* ---------------------------------------------------------------------- */ void Min::cleanup() { // reset reneighboring criteria neighbor->every = neigh_every; neighbor->delay = neigh_delay; neighbor->dist_check = neigh_dist_check; // delete fix at end of run, so its atom arrays won't persist modify->delete_fix("MINIMIZE"); } /* ---------------------------------------------------------------------- evaluate potential energy and forces may migrate atoms due to reneighboring return new energy, which should include nextra_global dof return negative gradient stored in atom->f return negative gradient for nextra_global dof in fextra ------------------------------------------------------------------------- */ double Min::energy_force(int resetflag) { // check for reneighboring // always communicate since minimizer moved atoms int nflag = neighbor->decide(); if (nflag == 0) { timer->stamp(); comm->forward_comm(); timer->stamp(TIME_COMM); } else { if (modify->n_min_pre_exchange) modify->min_pre_exchange(); 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(); if (atom->sortfreq > 0 && update->ntimestep >= atom->nextsort) atom->sort(); comm->borders(); if (triclinic) domain->lamda2x(atom->nlocal+atom->nghost); timer->stamp(TIME_COMM); neighbor->build(); timer->stamp(TIME_NEIGHBOR); } 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_comm(); 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 double energy = pe_compute->compute_scalar(); if (nextra_global) energy += modify->min_energy(fextra); if (output->thermo->normflag) energy /= atom->natoms; // if reneighbored, atoms migrated // if resetflag = 1, update x0 of atoms crossing PBC // reset vectors used by lo-level minimizer if (nflag) { if (resetflag) fix_minimize->reset_coords(); reset_vectors(); } return energy; } /* ---------------------------------------------------------------------- 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; } } } /* ---------------------------------------------------------------------- clear force on own & ghost atoms setup and clear other arrays as needed ------------------------------------------------------------------------- */ void Min::request(Pair *pair, int peratom, double maxvalue) { int n = nextra_atom + 1; xextra_atom = (double **) memory->srealloc(xextra_atom,n*sizeof(double *), "min:xextra_atom"); fextra_atom = (double **) memory->srealloc(fextra_atom,n*sizeof(double *), "min:fextra_atom"); extra_peratom = (int *) memory->srealloc(extra_peratom,n*sizeof(int), "min:extra_peratom"); extra_nlen = (int *) memory->srealloc(extra_nlen,n*sizeof(int), "min:extra_nlen"); extra_max = (double *) memory->srealloc(extra_max,n*sizeof(double), "min:extra_max"); requestor = (Pair **) memory->srealloc(requestor,n*sizeof(Pair *), "min:requestor"); requestor[nextra_atom] = pair; extra_peratom[nextra_atom] = peratom; extra_max[nextra_atom] = maxvalue; nextra_atom++; } /* ---------------------------------------------------------------------- */ void Min::modify_params(int narg, char **arg) { if (narg == 0) error->all("Illegal min_modify command"); int iarg = 0; while (iarg < narg) { if (strcmp(arg[iarg],"dmax") == 0) { if (iarg+2 > narg) error->all("Illegal min_modify command"); dmax = atof(arg[iarg+1]); iarg += 2; } else if (strcmp(arg[iarg],"line") == 0) { if (iarg+2 > narg) error->all("Illegal min_modify command"); if (strcmp(arg[iarg+1],"backtrack") == 0) linestyle = 0; else if (strcmp(arg[iarg+1],"quadratic") == 0) linestyle = 1; else error->all("Illegal min_modify command"); iarg += 2; } else error->all("Illegal min_modify command"); } } /* ---------------------------------------------------------------------- setup lists of computes for global and per-atom PE and pressure ------------------------------------------------------------------------- */ void Min::ev_setup() { delete [] elist_atom; delete [] vlist_global; delete [] vlist_atom; elist_atom = NULL; vlist_global = vlist_atom = NULL; nelist_atom = 0; nvlist_global = nvlist_atom = 0; for (int i = 0; i < modify->ncompute; i++) { if (modify->compute[i]->peatomflag) nelist_atom++; if (modify->compute[i]->pressflag) nvlist_global++; if (modify->compute[i]->pressatomflag) nvlist_atom++; } if (nelist_atom) elist_atom = new Compute*[nelist_atom]; if (nvlist_global) vlist_global = new Compute*[nvlist_global]; if (nvlist_atom) vlist_atom = new Compute*[nvlist_atom]; nelist_atom = 0; nvlist_global = nvlist_atom = 0; for (int i = 0; i < modify->ncompute; i++) { if (modify->compute[i]->peatomflag) elist_atom[nelist_atom++] = modify->compute[i]; if (modify->compute[i]->pressflag) vlist_global[nvlist_global++] = modify->compute[i]; if (modify->compute[i]->pressatomflag) vlist_atom[nvlist_atom++] = modify->compute[i]; } } /* ---------------------------------------------------------------------- set eflag,vflag for current iteration based on computes that need info on this ntimestep always set eflag_global = 1, since need energy every iteration eflag = 0 = no energy computation eflag = 1 = global energy only eflag = 2 = per-atom energy only eflag = 3 = both global and per-atom energy vflag = 0 = no virial computation (pressure) vflag = 1 = global virial with pair portion via sum of pairwise interactions vflag = 2 = global virial with pair portion via F dot r including ghosts vflag = 4 = per-atom virial only vflag = 5 or 6 = both global and per-atom virial ------------------------------------------------------------------------- */ void Min::ev_set(int ntimestep) { int i; int eflag_global = 1; int eflag_atom = 0; for (i = 0; i < nelist_atom; i++) if (elist_atom[i]->matchstep(ntimestep)) break; if (i < nelist_atom) eflag_atom = 2; if (eflag_global) update->eflag_global = update->ntimestep; if (eflag_atom) update->eflag_atom = update->ntimestep; eflag = eflag_global + eflag_atom; int vflag_global = 0; for (i = 0; i < nvlist_global; i++) if (vlist_global[i]->matchstep(ntimestep)) break; if (i < nvlist_global) vflag_global = virial_style; int vflag_atom = 0; for (i = 0; i < nvlist_atom; i++) if (vlist_atom[i]->matchstep(ntimestep)) break; if (i < nvlist_atom) vflag_atom = 4; if (vflag_global) update->vflag_global = update->ntimestep; if (vflag_atom) update->vflag_atom = update->ntimestep; vflag = vflag_global + vflag_atom; } /* ---------------------------------------------------------------------- compute and return ||force||_2^2 ------------------------------------------------------------------------- */ double Min::fnorm_sqr() { int i,n; double *fatom; double local_norm2_sqr = 0.0; - for (i = 0; i < n3; i++) local_norm2_sqr += f[i]*f[i]; + for (i = 0; i < nvec; i++) local_norm2_sqr += fvec[i]*fvec[i]; if (nextra_atom) { for (int m = 0; m < nextra_atom; m++) { fatom = fextra_atom[m]; n = extra_nlen[m]; for (i = 0; i < n; i++) local_norm2_sqr += fatom[i]*fatom[i]; } } double norm2_sqr = 0.0; MPI_Allreduce(&local_norm2_sqr,&norm2_sqr,1,MPI_DOUBLE,MPI_SUM,world); if (nextra_global) for (i = 0; i < nextra_global; i++) norm2_sqr += fextra[i]*fextra[i]; return norm2_sqr; } /* ---------------------------------------------------------------------- compute and return ||force||_inf ------------------------------------------------------------------------- */ double Min::fnorm_inf() { int i,n; double *fatom; double local_norm_inf = 0.0; - for (i = 0; i < n3; i++) - local_norm_inf = MAX(fabs(f[i]),local_norm_inf); + for (i = 0; i < nvec; i++) + local_norm_inf = MAX(fabs(fvec[i]),local_norm_inf); if (nextra_atom) { for (int m = 0; m < nextra_atom; m++) { fatom = fextra_atom[m]; n = extra_nlen[m]; for (i = 0; i < n; i++) local_norm_inf = MAX(fabs(fatom[i]),local_norm_inf); } } double norm_inf = 0.0; MPI_Allreduce(&local_norm_inf,&norm_inf,1,MPI_DOUBLE,MPI_MAX,world); if (nextra_global) for (i = 0; i < nextra_global; i++) norm_inf = MAX(fabs(fextra[i]),norm_inf); return norm_inf; } diff --git a/src/min.h b/src/min.h index 358880839..b73c4aa9c 100644 --- a/src/min.h +++ b/src/min.h @@ -1,102 +1,102 @@ /* ---------------------------------------------------------------------- LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator http://lammps.sandia.gov, Sandia National Laboratories Steve Plimpton, sjplimp@sandia.gov Copyright (2003) Sandia Corporation. Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains certain rights in this software. This software is distributed under the GNU General Public License. See the README file in the top-level LAMMPS directory. ------------------------------------------------------------------------- */ #ifndef LMP_MIN_H #define LMP_MIN_H #include "pointers.h" namespace LAMMPS_NS { class Min : protected Pointers { public: double einitial,efinal,eprevious; double fnorm2_init,fnorminf_init,fnorm2_final,fnorminf_final; double alpha_final; int niter,neval; char *stopstr; Min(class LAMMPS *); virtual ~Min(); void init(); void setup(); void setup_minimal(int); void run(int); void cleanup(); void request(class Pair *, int, double); double memory_usage() {return 0.0;} void modify_params(int, char **); virtual void init_style() {} virtual void setup_style() = 0; virtual void reset_vectors() = 0; virtual int iterate(int) = 0; protected: int eflag,vflag; // flags for energy/virial computation int virial_style; // compute virial explicitly or implicitly double dmax; // max dist to move any atom in one step int linestyle; // 0 = backtrack, 1 = quadratic int nelist_atom; // # of PE,virial computes to check int nvlist_global,nvlist_atom; class Compute **elist_atom; // list of PE,virial Computes class Compute **vlist_global; class Compute **vlist_atom; int pairflag,torqueflag; int triclinic; // 0 if domain is orthog, 1 if triclinic int narray; // # of arrays stored by fix_minimize class FixMinimize *fix_minimize; // fix that stores auxiliary data class Compute *pe_compute; // compute for potential energy double ecurrent; // current potential energy double ndoftotal; // total dof for entire problem - int n3; // local atomic dof - double *x; // variables for atomic dof, as 1d vector - double *f; // force vector for atomic dof, as 1d vector + int nvec; // local atomic dof = length of xvec + double *xvec; // variables for atomic dof, as 1d vector + double *fvec; // force vector for atomic dof, as 1d vector int nextra_global; // # of extra global dof due to fixes double *fextra; // force vector for extra global dof // xextra is stored by fix int nextra_atom; // # of sets of extra per-atom dof double **xextra_atom; // variables for extra per-atom dof sets double **fextra_atom; // force vectors for extra per-atom dof sets int *extra_peratom; // # of per-atom values in each set int *extra_nlen; // total local length of each set, e.g 3*nlocal double *extra_max; // max allowed change in one iter for each set class Pair **requestor; // Pair that requested each set int neigh_every,neigh_delay,neigh_dist_check; // neighboring params double energy_force(int); void force_clear(); double compute_force_norm_sqr(); double compute_force_norm_inf(); void ev_setup(); void ev_set(int); double fnorm_sqr(); double fnorm_inf(); }; } #endif diff --git a/src/min_cg.cpp b/src/min_cg.cpp index 327f35ede..252f19fd1 100644 --- a/src/min_cg.cpp +++ b/src/min_cg.cpp @@ -1,189 +1,189 @@ /* ---------------------------------------------------------------------- LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator http://lammps.sandia.gov, Sandia National Laboratories Steve Plimpton, sjplimp@sandia.gov Copyright (2003) Sandia Corporation. Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains certain rights in this software. This software is distributed under the GNU General Public License. See the README file in the top-level LAMMPS directory. ------------------------------------------------------------------------- */ #include "math.h" #include "string.h" #include "mpi.h" #include "min_cg.h" #include "atom.h" #include "update.h" #include "output.h" #include "timer.h" using namespace LAMMPS_NS; #define MAXATOMS 0x7FFFFFFF // EPS_ENERGY = minimum normalization for energy tolerance #define EPS_ENERGY 1.0e-8 // 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) /* ---------------------------------------------------------------------- */ MinCG::MinCG(LAMMPS *lmp) : MinLineSearch(lmp) {} /* ---------------------------------------------------------------------- minimization via conjugate gradient iterations ------------------------------------------------------------------------- */ int MinCG::iterate(int niter_max) { int i,m,n,fail,ntimestep; double beta,gg,dot[2],dotall[2]; double *fatom,*gatom,*hatom; // nlimit = max # of CG iterations before restarting // set to ndoftotal unless too big int nlimit = static_cast (MIN(MAXATOMS,ndoftotal)); // initialize working vectors - for (i = 0; i < n3; i++) h[i] = g[i] = f[i]; + for (i = 0; i < nvec; i++) h[i] = g[i] = fvec[i]; if (nextra_atom) for (m = 0; m < nextra_atom; m++) { fatom = fextra_atom[m]; gatom = gextra_atom[m]; hatom = hextra_atom[m]; n = extra_nlen[m]; for (i = 0; i < n; i++) hatom[i] = gatom[i] = fatom[i]; } if (nextra_global) for (i = 0; i < nextra_global; i++) hextra[i] = gextra[i] = fextra[i]; gg = fnorm_sqr(); neval = 0; for (niter = 0; niter < niter_max; niter++) { ntimestep = ++update->ntimestep; // line minimization along direction h from current atom->x eprevious = ecurrent; fail = (this->*linemin)(ecurrent,alpha_final,neval); if (fail) return fail; // function evaluation criterion if (neval >= update->max_eval) return MAXEVAL; // energy tolerance criterion if (fabs(ecurrent-eprevious) < update->etol * 0.5*(fabs(ecurrent) + fabs(eprevious) + EPS_ENERGY)) return ETOL; // force tolerance criterion dot[0] = dot[1] = 0.0; - for (i = 0; i < n3; i++) { - dot[0] += f[i]*f[i]; - dot[1] += f[i]*g[i]; + for (i = 0; i < nvec; i++) { + dot[0] += fvec[i]*fvec[i]; + dot[1] += fvec[i]*g[i]; } if (nextra_atom) for (m = 0; m < nextra_atom; m++) { fatom = fextra_atom[m]; gatom = gextra_atom[m]; n = extra_nlen[m]; for (i = 0; i < n; i++) { dot[0] += fatom[i]*fatom[i]; dot[1] += fatom[i]*gatom[i]; } } MPI_Allreduce(dot,dotall,2,MPI_DOUBLE,MPI_SUM,world); if (nextra_global) for (i = 0; i < nextra_global; i++) { dotall[0] += fextra[i]*fextra[i]; dotall[1] += fextra[i]*gextra[i]; } if (dotall[0] < update->ftol*update->ftol) return FTOL; // update new search direction h from new f = -Grad(x) and old g // this is Polak-Ribieri formulation // beta = dotall[0]/gg would be Fletcher-Reeves // reinitialize CG every ndof iterations by setting beta = 0.0 beta = MAX(0.0,(dotall[0] - dotall[1])/gg); if ((niter+1) % nlimit == 0) beta = 0.0; gg = dotall[0]; - for (i = 0; i < n3; i++) { - g[i] = f[i]; + for (i = 0; i < nvec; i++) { + g[i] = fvec[i]; h[i] = g[i] + beta*h[i]; } if (nextra_atom) for (m = 0; m < nextra_atom; m++) { fatom = fextra_atom[m]; gatom = gextra_atom[m]; hatom = hextra_atom[m]; n = extra_nlen[m]; for (i = 0; i < n; i++) { gatom[i] = fatom[i]; hatom[i] = gatom[i] + beta*hatom[i]; } } if (nextra_global) for (i = 0; i < nextra_global; i++) { gextra[i] = fextra[i]; hextra[i] = gextra[i] + beta*hextra[i]; } // reinitialize CG if new search direction h is not downhill dot[0] = 0.0; - for (i = 0; i < n3; i++) dot[0] += g[i]*h[i]; + for (i = 0; i < nvec; i++) dot[0] += g[i]*h[i]; if (nextra_atom) for (m = 0; m < nextra_atom; m++) { gatom = gextra_atom[m]; hatom = hextra_atom[m]; n = extra_nlen[m]; for (i = 0; i < n; i++) dot[0] += gatom[i]*hatom[i]; } MPI_Allreduce(dot,dotall,1,MPI_DOUBLE,MPI_SUM,world); if (nextra_global) for (i = 0; i < nextra_global; i++) dotall[0] += gextra[i]*hextra[i]; if (dotall[0] <= 0.0) { - for (i = 0; i < n3; i++) h[i] = g[i]; + for (i = 0; i < nvec; i++) h[i] = g[i]; if (nextra_atom) for (m = 0; m < nextra_atom; m++) { gatom = gextra_atom[m]; hatom = hextra_atom[m]; n = extra_nlen[m]; for (i = 0; i < n; i++) hatom[i] = gatom[i]; } if (nextra_global) for (i = 0; i < nextra_global; i++) hextra[i] = gextra[i]; } // output for thermo, dump, restart files if (output->next == ntimestep) { timer->stamp(); output->write(ntimestep); timer->stamp(TIME_OUTPUT); } } return MAXITER; } diff --git a/src/min_hftn.cpp b/src/min_hftn.cpp index a8be7d019..f0b70d4b6 100644 --- a/src/min_hftn.cpp +++ b/src/min_hftn.cpp @@ -1,1680 +1,1679 @@ /* ---------------------------------------------------------------------- LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator http://lammps.sandia.gov, Sandia National Laboratories Steve Plimpton, sjplimp@sandia.gov Copyright (2003) Sandia Corporation. Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains certain rights in this software. This software is distributed under the GNU General Public License. See the README file in the top-level LAMMPS directory. ------------------------------------------------------------------------- */ /* ---------------------------------------------------------------------- Author: Todd Plantenga (SNL) Sources: "Numerical Optimization", Nocedal and Wright, 2nd Ed, p170 "Parallel Unconstrained Min", Plantenga, SAND98-8201 ------------------------------------------------------------------------- */ #include "math.h" #include "string.h" #include "atom.h" #include "fix_minimize.h" #include "min_hftn.h" #include "modify.h" #include "output.h" #include "pair.h" #include "update.h" #include "timer.h" #define MIN(A,B) (((A) < (B)) ? (A) : (B)) #define MAX(A,B) (((A) > (B)) ? (A) : (B)) using namespace LAMMPS_NS; /* ---------------------------------------------------------------------- * This class performs Hessian-free truncated Newton minimization on an * unconstrained molecular potential. The algorithm avoids computing the * Hessian matrix, but obtains a near-quadratic rate of convergence. ------------------------------------------------------------------------- */ /* ---------------------------------------------------------------------- File local data ------------------------------------------------------------------------- */ //---- CONSTANTS MAP TO stopstrings DECLARED IN Min.run (min.cpp). static const int STOP_MAX_ITER = 0; //-- MAX ITERATIONS EXCEEDED static const int STOP_MAX_FORCE_EVALS = 1; //-- MAX FORCE EVALUATIONS EXCEEDED static const int STOP_ENERGY_TOL = 2; //-- STEP DID NOT CHANGE ENERGY static const int STOP_FORCE_TOL = 3; //-- CONVERGED TO DESIRED FORCE TOL static const int STOP_TR_TOO_SMALL = 8; //-- TRUST REGION TOO SMALL static const int STOP_ERROR = 9; //-- INTERNAL ERROR static const int NO_CGSTEP_BECAUSE_F_TOL_SATISFIED = 0; static const int CGSTEP_NEWTON = 1; static const int CGSTEP_TO_TR = 2; static const int CGSTEP_TO_DMAX = 3; static const int CGSTEP_NEGATIVE_CURVATURE = 4; static const int CGSTEP_MAX_INNER_ITERS = 5; static const int CGSTEP_UNDETERMINED = 6; //---- WHEN TESTING ENERGY_TOL, THE ENERGY MAGNITUDE MUST BE AT LEAST THIS BIG. static const double MIN_ETOL_MAG = 1.0e-8; //---- MACHINE PRECISION IS SOMETIMES DEFINED BY THE C RUNTIME. #ifdef DBL_EPSILON #define MACHINE_EPS DBL_EPSILON #else #define MACHINE_EPS 2.220446049250313e-16 #endif /* ---------------------------------------------------------------------- Constructor ------------------------------------------------------------------------- */ MinHFTN::MinHFTN(LAMMPS *lmp) : Min(lmp) { for (int i = 1; i < NUM_HFTN_ATOM_BASED_VECTORS; i++) _daExtraGlobal[i] = NULL; for (int i = 0; i < NUM_HFTN_ATOM_BASED_VECTORS; i++) _daExtraAtom[i] = NULL; _fpPrint = NULL; return; } /* ---------------------------------------------------------------------- Destructor ------------------------------------------------------------------------- */ MinHFTN::~MinHFTN (void) { for (int i = 1; i < NUM_HFTN_ATOM_BASED_VECTORS; i++) if (_daExtraGlobal[i] != NULL) delete [] _daExtraGlobal[i]; for (int i = 0; i < NUM_HFTN_ATOM_BASED_VECTORS; i++) if (_daExtraAtom[i] != NULL) delete [] _daExtraAtom[i]; return; } /* ---------------------------------------------------------------------- Public method init_style ------------------------------------------------------------------------- */ void MinHFTN::init_style() { for (int i = 1; i < NUM_HFTN_ATOM_BASED_VECTORS; i++) { if (_daExtraGlobal[i] != NULL) delete [] _daExtraGlobal[i]; _daExtraGlobal[i] = NULL; } for (int i = 0; i < NUM_HFTN_ATOM_BASED_VECTORS; i++) { if (_daExtraAtom[i] != NULL) delete [] _daExtraAtom[i]; _daExtraAtom[i] = NULL; } return; } /* ---------------------------------------------------------------------- Public method setup_style ------------------------------------------------------------------------- */ void MinHFTN::setup_style() { //---- ALLOCATE MEMORY FOR ATOMIC DEGREES OF FREEDOM. for (int i = 0; i < NUM_HFTN_ATOM_BASED_VECTORS; i++) fix_minimize->add_vector(3); //---- ALLOCATE MEMORY FOR EXTRA GLOBAL DEGREES OF FREEDOM. //---- THE FIX MODULE TAKES CARE OF THE FIRST VECTOR, X0 (XK). if (nextra_global) { for (int i = 1; i < NUM_HFTN_ATOM_BASED_VECTORS; i++) _daExtraGlobal[i] = new double[nextra_global]; } //---- ALLOCATE MEMORY FOR EXTRA PER-ATOM DEGREES OF FREEDOM. if (nextra_atom) { for (int i = 0; i < NUM_HFTN_ATOM_BASED_VECTORS; i++) _daExtraAtom[i] = new double*[nextra_atom]; for (int m = 0; m < nextra_atom; m++) { for (int i = 0; i < NUM_HFTN_ATOM_BASED_VECTORS; i++) fix_minimize->add_vector (extra_peratom[m]); } } return; } /* ---------------------------------------------------------------------- Public method reset_vectors After an energy/force calculation, atoms may migrate from one processor to another. Any local vector correlated with atom positions or forces must also be migrated. This is accomplished by a subclass of Fix. This method updates local pointers to the latest Fix copies. ------------------------------------------------------------------------- */ void MinHFTN::reset_vectors() { - n3 = 3 * atom->nlocal; + nvec = 3 * atom->nlocal; //---- ATOMIC DEGREES OF FREEDOM. - if (n3 > 0) { - x = atom->x[0]; - f = atom->f[0]; + if (nvec > 0) { + xvec = atom->x[0]; + fvec = atom->f[0]; } for (int i = 0; i < NUM_HFTN_ATOM_BASED_VECTORS; i++) _daAVectors[i] = fix_minimize->request_vector (i); //---- EXTRA PER-ATOM DEGREES OF FREEDOM. if (nextra_atom) { int n = NUM_HFTN_ATOM_BASED_VECTORS; for (int m = 0; m < nextra_atom; m++) { extra_nlen[m] = extra_peratom[m] * atom->nlocal; requestor[m]->min_pointers (&xextra_atom[m], &fextra_atom[m]); for (int i = 0; i < NUM_HFTN_ATOM_BASED_VECTORS; i++) _daExtraAtom[i][m] = fix_minimize->request_vector (n++); } } return; } /* ---------------------------------------------------------------------- Public method iterate Upon entry, Min::setup() and Min::run have executed, and energy has already been evaluated at the initial point. Return an integer code that maps to a stop condition in min.cpp. ------------------------------------------------------------------------- */ int MinHFTN::iterate(int) { //---- TURN THIS ON TO GENERATE AN OPTIMIZATION PROGRESS FILE. bool bPrintProgress = false; if (bPrintProgress) open_hftn_print_file_(); double dFinalEnergy = 0.0; double dFinalFnorm2 = 0.0; modify->min_clearstore(); int nStopCode = execute_hftn_ (bPrintProgress, einitial, fnorm2_init, dFinalEnergy, dFinalFnorm2); modify->min_clearstore(); if (bPrintProgress) close_hftn_print_file_(); return( nStopCode ); } /* ---------------------------------------------------------------------- Private method execute_hftn_ @param[in] bPrintProgress - if true then print progress to a file @param[in] dInitialEnergy - energy at input x @param[in] dInitialForce2 - |F|_2 at input x @param[out] dFinalEnergy - energy at output x @param[out] dFinalForce2 - |F|_2 at output x Return stop code described in the enumeration at the top of this file, and the following: atom->x - positions at output x atom->f - forces evaluated at output x ------------------------------------------------------------------------- */ int MinHFTN::execute_hftn_(const bool bPrintProgress, const double dInitialEnergy, const double dInitialForce2, double & dFinalEnergy, double & dFinalForce2) { //---- DEFINE OUTPUTS PRINTED BY "Finish". eprevious = dInitialEnergy; alpha_final = 0.0; neval = 0; niter = 0; dFinalEnergy = dInitialEnergy; dFinalForce2 = dInitialForce2; if (dInitialForce2 < update->ftol) return( STOP_FORCE_TOL ); //---- SAVE ATOM POSITIONS BEFORE AN ITERATION. fix_minimize->store_box(); - for (int i = 0; i < n3; i++) - _daAVectors[VEC_XK][i] = atom->x[0][i]; + for (int i = 0; i < nvec; i++) + _daAVectors[VEC_XK][i] = xvec[i]; if (nextra_atom) { for (int m = 0; m < nextra_atom; m++) { double * xatom = xextra_atom[m]; double * xkAtom = _daExtraAtom[VEC_XK][m]; int n = extra_nlen[m]; for (int i = 0; i < n; i++) xkAtom[i] = xatom[i]; } } if (nextra_global) modify->min_store(); double dXInf = calc_xinf_using_mpi_(); //---- FIND THE NUMBER OF UNKNOWNS. - int nLocalNumUnknowns = n3 + nextra_atom; + int nLocalNumUnknowns = nvec + nextra_atom; MPI_Allreduce (&nLocalNumUnknowns, &_nNumUnknowns, 1, MPI_INT, MPI_SUM, world); //---- INITIALIZE THE TRUST RADIUS BASED ON THE GRADIENT. double dTrustRadius = 1.5 * dInitialForce2; //---- TRUST RADIUS MUST KEEP STEPS FROM LETTING ATOMS MOVE SO FAR THEY //---- VIOLATE PHYSICS OR JUMP BEYOND A PARALLEL PROCESSING DOMAIN. //---- LINE SEARCH METHODS DO THIS BY RESTRICTING THE LARGEST CHANGE //---- OF ANY ATOM'S COMPONENT TO dmax. AN EXACT CHECK IS MADE LATER, //---- BUT THIS GUIDES DETERMINATION OF A MAX TRUST RADIUS. double dMaxTrustRadius = dmax * sqrt (_nNumUnknowns); dTrustRadius = MIN (dTrustRadius, dMaxTrustRadius); double dLastNewtonStep2 = dMaxTrustRadius; if (bPrintProgress) hftn_print_line_ (false, -1, neval, dInitialEnergy, dInitialForce2, -1, dTrustRadius, 0.0, 0.0, 0.0); bool bHaveEvaluatedAtX = true; double dCurrentEnergy = dInitialEnergy; double dCurrentForce2 = dInitialForce2; for (niter = 0; niter < update->nsteps; niter++) { (update->ntimestep)++; //---- CALL THE INNER LOOP TO GET THE NEXT TRUST REGION STEP. double dCgForce2StopTol = MIN ((dCurrentForce2 / 2.0), 0.1 / (niter+1)); dCgForce2StopTol = MAX (dCgForce2StopTol, update->ftol); double dNewEnergy; double dNewForce2; int nStepType; double dStepLength2; double dStepLengthInf; if (compute_inner_cg_step_ (dTrustRadius, dCgForce2StopTol, update->max_eval, bHaveEvaluatedAtX, dCurrentEnergy, dCurrentForce2, dNewEnergy, dNewForce2, nStepType, dStepLength2, dStepLengthInf) == false) { //---- THERE WAS AN ERROR. RESTORE TO LAST ACCEPTED STEP. if (nextra_global) modify->min_step (0.0, _daExtraGlobal[VEC_CG_P]); - for (int i = 0; i < n3; i++) - atom->x[0][i] = _daAVectors[VEC_XK][i]; + for (int i = 0; i < nvec; i++) + xvec[i] = _daAVectors[VEC_XK][i]; if (nextra_atom) { for (int m = 0; m < nextra_atom; m++) { double * xatom = xextra_atom[m]; double * xkAtom = _daExtraAtom[VEC_XK][m]; int n = extra_nlen[m]; for (int i = 0; i < n; i++) xatom[i] = xkAtom[i]; } } dFinalEnergy = energy_force (0); neval++; dFinalForce2 = sqrt (fnorm_sqr()); return( STOP_ERROR ); } //---- STOP IF THE CURRENT POSITION WAS FOUND TO BE ALREADY GOOD ENOUGH. //---- IN THIS CASE THE ENERGY AND FORCES ARE ALREADY COMPUTED. if (nStepType == NO_CGSTEP_BECAUSE_F_TOL_SATISFIED) { if (bPrintProgress) hftn_print_line_ (true, niter+1, neval, dNewEnergy, dNewForce2, nStepType, dTrustRadius, dStepLength2, 0.0, 0.0); dFinalEnergy = dNewEnergy; dFinalForce2 = dNewForce2; return( STOP_FORCE_TOL ); } //---- COMPUTE THE DIRECTIONAL DERIVATIVE H(x_k) p. bool bUseForwardDiffs = (dCurrentForce2 > 1000.0 * sqrt (MACHINE_EPS)); evaluate_dir_der_ (bUseForwardDiffs, VEC_CG_P, VEC_CG_HD, true, dCurrentEnergy); //---- COMPUTE p^T grad(x_k) AND SAVE IT FOR PRED. double dGradDotP = calc_grad_dot_v_using_mpi_ (VEC_CG_P); //---- MOVE TO THE NEW POINT AND EVALUATE ENERGY AND FORCES. //---- THIS IS THE PLACE WHERE energy_force IS ALLOWED TO RESET. - for (int i = 0; i < n3; i++) - atom->x[0][i] = _daAVectors[VEC_XK][i] + _daAVectors[VEC_CG_P][i]; + for (int i = 0; i < nvec; i++) + xvec[i] = _daAVectors[VEC_XK][i] + _daAVectors[VEC_CG_P][i]; if (nextra_atom) { for (int m = 0; m < nextra_atom; m++) { double * xatom = xextra_atom[m]; double * xkAtom = _daExtraAtom[VEC_XK][m]; double * pAtom = _daExtraAtom[VEC_CG_P][m]; int n = extra_nlen[m]; for (int i = 0; i < n; i++) xatom[i] = xkAtom[i] + pAtom[i]; } } if (nextra_global) modify->min_step (1.0, _daExtraGlobal[VEC_CG_P]); dNewEnergy = energy_force (1); neval++; dNewForce2 = sqrt (fnorm_sqr()); double dAred = dCurrentEnergy - dNewEnergy; //---- STOP IF THE FORCE TOLERANCE IS MET. if (dNewForce2 < update->ftol) { if (bPrintProgress) hftn_print_line_ (true, niter+1, neval, dNewEnergy, dNewForce2, nStepType, dTrustRadius, dStepLength2, dAred, -1.0); //---- (IMPLICITLY ACCEPT THE LAST STEP TO THE NEW POINT.) dFinalEnergy = dNewEnergy; dFinalForce2 = dNewForce2; return( STOP_FORCE_TOL ); } //---- STOP IF THE ACTUAL ENERGY REDUCTION IS TINY. if (nStepType != CGSTEP_TO_DMAX) { double dMag = 0.5 * (fabs (dCurrentEnergy) + fabs (dNewEnergy)); dMag = MAX (dMag, MIN_ETOL_MAG); if ( (fabs (dAred) < (update->etol * dMag)) || (dStepLengthInf == 0.0) ) { if (bPrintProgress) hftn_print_line_ (true, niter+1, neval, dNewEnergy, dNewForce2, nStepType, dTrustRadius, dStepLength2, dAred, -1.0); //---- (IMPLICITLY ACCEPT THE LAST STEP TO THE NEW POINT.) dFinalEnergy = dNewEnergy; dFinalForce2 = dNewForce2; return( STOP_ENERGY_TOL ); } } //---- COMPUTE THE PREDICTED REDUCTION - p^T grad - 0.5 p^T Hp double dPHP = calc_dot_prod_using_mpi_ (VEC_CG_P, VEC_CG_HD); double dPred = - dGradDotP - (0.5 * dPHP); //---- ACCEPT OR REJECT THE STEP PROPOSED BY THE INNER CG LOOP. //---- WHEN NEAR A SOLUTION, THE FORCE NORM IS PROBABLY MORE ACCURATE, //---- SO DON'T ACCEPT A STEP THAT REDUCES ENERGY SOME TINY AMOUNT //---- WHILE INCREASING THE FORCE NORM. bool bStepAccepted = (dAred > 0.0) && ( (dNewForce2 < dCurrentForce2) || (dCurrentForce2 > 1.0e-6)); if (bStepAccepted) { //---- THE STEP IS ACCEPTED. if (bPrintProgress) hftn_print_line_ (true, niter+1, neval, dNewEnergy, dNewForce2, nStepType, dTrustRadius, dStepLength2, dAred, dPred); fix_minimize->store_box(); modify->min_clearstore(); - for (int i = 0; i < n3; i++) - _daAVectors[VEC_XK][i] = atom->x[0][i]; + for (int i = 0; i < nvec; i++) + _daAVectors[VEC_XK][i] = xvec[i]; if (nextra_atom) { for (int m = 0; m < nextra_atom; m++) { double * xatom = xextra_atom[m]; double * xkAtom = _daExtraAtom[VEC_XK][m]; int n = extra_nlen[m]; for (int i = 0; i < n; i++) xkAtom[i] = xatom[i]; } } if (nextra_global) modify->min_store(); if (niter > 0) eprevious = dCurrentEnergy; dCurrentEnergy = dNewEnergy; dCurrentForce2 = dNewForce2; bHaveEvaluatedAtX = true; if (nStepType == CGSTEP_NEWTON) dLastNewtonStep2 = dStepLength2; //---- UPDATE THE TRUST REGION BASED ON AGREEMENT BETWEEN //---- THE ACTUAL REDUCTION AND THE PREDICTED MODEL REDUCTION. if ((dAred > 0.75 * dPred) && (dStepLength2 >= 0.99 * dTrustRadius)) dTrustRadius = 2.0 * dTrustRadius; dTrustRadius = MIN (dTrustRadius, dMaxTrustRadius); //---- DMAX VIOLATIONS TRUNCATE THE CG STEP WITHOUT COMPARISONS; //---- BETTER TO ADJUST THE TRUST REGION SO DMAX STOPS HAPPENING. if (nStepType == CGSTEP_TO_DMAX) { if (dStepLength2 <= MACHINE_EPS) dTrustRadius = 0.1 * dTrustRadius; else dTrustRadius = MIN (dTrustRadius, 2.0 * dStepLength2); } } else { //---- THE STEP IS REJECTED. if (bPrintProgress) hftn_print_line_ (false, niter+1, neval, dCurrentEnergy, dCurrentForce2, nStepType, dTrustRadius, dStepLength2, dAred, dPred); //---- RESTORE THE LAST X_K POSITION. if (nextra_global) modify->min_step (0.0, _daExtraGlobal[VEC_CG_P]); - for (int i = 0; i < n3; i++) - atom->x[0][i] = _daAVectors[VEC_XK][i]; + for (int i = 0; i < nvec; i++) + xvec[i] = _daAVectors[VEC_XK][i]; if (nextra_atom) { for (int m = 0; m < nextra_atom; m++) { double * xatom = xextra_atom[m]; double * xkAtom = _daExtraAtom[VEC_XK][m]; int n = extra_nlen[m]; for (int i = 0; i < n; i++) xatom[i] = xkAtom[i]; } } bHaveEvaluatedAtX = false; //---- UPDATE THE TRUST REGION. //---- EXPERIMENTS INDICATE NEGATIVE CURVATURE CAN TAKE A BAD //---- STEP A LONG WAY, SO BE MORE AGGRESSIVE IN THIS CASE. //---- ALSO, IF NEAR A SOLUTION AND DONE WITH NEWTON STEPS, //---- THEN REDUCE TO SOMETHING NEAR THE LAST GOOD NEWTON STEP. if ((nStepType == CGSTEP_NEGATIVE_CURVATURE) && (-dAred > dPred)) dTrustRadius = 0.10 * MIN (dTrustRadius, dStepLength2); else if ( (nStepType == CGSTEP_TO_DMAX) && (dStepLength2 <= MACHINE_EPS)) dTrustRadius = 0.10 * dTrustRadius; else if (-dAred > dPred) dTrustRadius = 0.20 * MIN (dTrustRadius, dStepLength2); else dTrustRadius = 0.25 * MIN (dTrustRadius, dStepLength2); if ( (nStepType != CGSTEP_NEWTON) && (dCurrentForce2 < sqrt (MACHINE_EPS))) dTrustRadius = MIN (dTrustRadius, 2.0 * dLastNewtonStep2); dLastNewtonStep2 = dMaxTrustRadius; //---- STOP IF THE TRUST RADIUS IS TOO SMALL TO CONTINUE. if ( (dTrustRadius <= 0.0) || (dTrustRadius <= MACHINE_EPS * MAX (1.0, dXInf))) { dFinalEnergy = dCurrentEnergy; dFinalForce2 = dCurrentForce2; return( STOP_TR_TOO_SMALL ); } } //---- OUTPUT FOR thermo, dump, restart FILES. if (output->next == update->ntimestep) { //---- IF THE LAST STEP WAS REJECTED, THEN REEVALUATE ENERGY AND //---- FORCES AT THE OLD POINT SO THE OUTPUT DOES NOT DISPLAY //---- THE INCREASED ENERGY OF THE REJECTED STEP. if (bStepAccepted == false) { dCurrentEnergy = energy_force (1); neval++; } timer->stamp(); output->write (update->ntimestep); timer->stamp (TIME_OUTPUT); } //---- RETURN IF NUMBER OF EVALUATIONS EXCEEDED. if (neval >= update->max_eval) { dFinalEnergy = dCurrentEnergy; dFinalForce2 = dCurrentForce2; return( STOP_MAX_FORCE_EVALS ); } } //-- END for LOOP OVER niter dFinalEnergy = dCurrentEnergy; dFinalForce2 = dCurrentForce2; return( STOP_MAX_ITER ); } /* ---------------------------------------------------------------------- Private method compute_inner_cg_step_ Execute CG using Hessian-vector products approximated by finite difference directional derivatives. On input these must be defined: atom->x - positions at x atom->f - ignored VEC_XK - positions at x On output these are defined: atom->x - unchanged atom->f - forces evaluated at x, but only if nStepType == NO_CGSTEP VEC_XK - unchanged VEC_CG_P - step from VEC_XK to new positions During processing these are modified: VEC_CG_D - conjugate gradient inner loop step VEC_CG_HD - Hessian-vector product VEC_CG_R - residual of inner loop step VEC_DIF1 - temp storage VEC_DIF2 - temp storage @param[in] dTrustRadius - trust region radius for this subiteration @param[in] dForceTol - stop tolerance on |F|_2 for this subiteration @param[in] nMaxEvals - total energy/force evaluations allowed @param[in] bHaveEvalAtXin - true if forces are valid at input x @param[in] dEnergyAtXin - energy at input x, if bHaveEvalAtXin is true @param[in] dForce2AtXin - |F|_2 at input x, if bHaveEvalAtXin is true @param[out] dEnergyAtXout - energy at output x, if NO_CGSTEP (see below) @param[out] dForce2AtXout - |F|_2 at output x, if NO_CGSTEP (see below) @param[out] nStepType - step type for hftn_print_line_() @param[out] dStepLength2 - |step|_2 @param[out] dStepLengthInf - |step|_inf Return false if there was a fatal error. If nStepType equals NO_CGSTEP_BECAUSE_F_TOL_SATISFIED, then the energy and forces are evaluated and returned in dEnergyAtXout, dForce2AtXout; else energy and forces are not evaluated. ------------------------------------------------------------------------- */ bool MinHFTN::compute_inner_cg_step_(const double dTrustRadius, const double dForceTol, const int nMaxEvals, const bool bHaveEvalAtXin, const double dEnergyAtXin, const double dForce2AtXin, double & dEnergyAtXout, double & dForce2AtXout, int & nStepType, double & dStepLength2, double & dStepLengthInf) { //---- SET p_0 = 0. if (nextra_global) { for (int i = 0; i < nextra_global; i++) _daExtraGlobal[VEC_CG_P][i] = 0.0; } - for (int i = 0; i < n3; i++) + for (int i = 0; i < nvec; i++) _daAVectors[VEC_CG_P][i] = 0.0; if (nextra_atom) { for (int m = 0; m < nextra_atom; m++) { double * pAtom = _daExtraAtom[VEC_CG_P][m]; int n = extra_nlen[m]; for (int i = 0; i < n; i++) pAtom[i] = 0.0; } } double dPP = 0.0; //---- OBTAIN THE ENERGY AND FORCES AT THE INPUT POSITION. double dEnergyAtX = dEnergyAtXin; double dForce2AtX = dForce2AtXin; if (bHaveEvalAtXin == false) { dEnergyAtX = energy_force (0); neval++; dForce2AtX = sqrt (fnorm_sqr()); } //---- RETURN IMMEDIATELY IF THE FORCE TOLERANCE IS ALREADY MET. //---- THE STEP TYPE INFORMS THE CALLER THAT ENERGY AND FORCES HAVE //---- BEEN EVALUATED. if (dForce2AtX <= dForceTol) { dEnergyAtXout = dEnergyAtX; dForce2AtXout = dForce2AtX; nStepType = NO_CGSTEP_BECAUSE_F_TOL_SATISFIED; dStepLength2 = 0.0; dStepLengthInf = 0.0; return( true ); } //---- r_0 = -grad (FIRST SEARCH DIRECTION IS STEEPEST DESCENT) //---- d_0 = r_0 //---- REMEMBER THAT FORCES = -GRADIENT. if (nextra_global) { for (int i = 0; i < nextra_global; i++) { _daExtraGlobal[VEC_CG_R][i] = fextra[i]; _daExtraGlobal[VEC_CG_D][i] = fextra[i]; } } - for (int i = 0; i < n3; i++) { - _daAVectors[VEC_CG_R][i] = atom->f[0][i]; - _daAVectors[VEC_CG_D][i] = atom->f[0][i]; + for (int i = 0; i < nvec; i++) { + _daAVectors[VEC_CG_R][i] = fvec[i]; + _daAVectors[VEC_CG_D][i] = fvec[i]; } if (nextra_atom) { for (int m = 0; m < nextra_atom; m++) { double * rAtom = _daExtraAtom[VEC_CG_R][m]; double * dAtom = _daExtraAtom[VEC_CG_D][m]; int n = extra_nlen[m]; for (int i = 0; i < n; i++) { rAtom[i] = fextra_atom[m][i]; dAtom[i] = fextra_atom[m][i]; } } } double dRR = dForce2AtX * dForce2AtX; double dR0norm2 = sqrt (dRR); //---- LIMIT THE NUMBER OF INNER CG ITERATIONS. //---- BASE IT ON THE NUMBER OF UNKNOWNS, OR MAXIMUM EVALUATIONS ASSUMING //---- FORWARD DIFFERENCES ARE USED. //---- NOTE THAT SETTING MAX=1 GIVES STEEPEST DESCENT. int nLimit1 = _nNumUnknowns / 5; if (nLimit1 < 100) nLimit1 = MIN (_nNumUnknowns, 100); int nLimit2 = (nMaxEvals - neval) / 2; int nMaxInnerIters = MIN (nLimit1, nLimit2); //---- FURTHER LIMIT ITERATIONS IF NEAR MACHINE ROUNDOFF. //---- THE METHOD CAN WASTE A LOT EVALUATIONS WITH LITTLE PAYOFF PROSPECT. if (dForce2AtX < (sqrt (MACHINE_EPS) * MAX (1.0, fabs (dEnergyAtX))) ) nMaxInnerIters = MIN (nMaxInnerIters, _nNumUnknowns / 20); bool bUseForwardDiffs = (dForce2AtX > 1000.0 * sqrt (MACHINE_EPS)); //---- MAIN CG LOOP. for (int nInnerIter = 0; nInnerIter < nMaxInnerIters; nInnerIter++) { //---- COMPUTE HESSIAN-VECTOR PRODUCT: H(x_k) d_i. double dDummyEnergy; evaluate_dir_der_ (bUseForwardDiffs, VEC_CG_D, VEC_CG_HD, false, dDummyEnergy); //---- CALCULATE d_i^T H d_i AND d_i^T d_i. double dDHD; double dDD; calc_dhd_dd_using_mpi_ (dDHD, dDD); //---- HANDLE NEGATIVE CURVATURE. if (dDHD <= (MACHINE_EPS * dDD)) { //---- PROJECT BOTH DIRECTIONS TO THE TRUST RADIUS AND DECIDE //---- WHICH MAKES A BETTER PREDICTED REDUCTION. //---- p_i^T H(x_k) d_i AND grad_i^T d_i. double dPdotD = calc_dot_prod_using_mpi_ (VEC_CG_P, VEC_CG_D); double dPdotHD = calc_dot_prod_using_mpi_ (VEC_CG_P, VEC_CG_HD); //---- MOVE TO X_K AND COMPUTE ENERGY AND FORCES. if (nextra_global) modify->min_step (0.0, _daExtraGlobal[VEC_CG_P]); - for (int i = 0; i < n3; i++) - atom->x[0][i] = _daAVectors[VEC_XK][i]; + for (int i = 0; i < nvec; i++) + xvec[i] = _daAVectors[VEC_XK][i]; if (nextra_atom) { for (int m = 0; m < nextra_atom; m++) { double * xatom = xextra_atom[m]; double * xkAtom = _daExtraAtom[VEC_XK][m]; int n = extra_nlen[m]; for (int i = 0; i < n; i++) xatom[i] = xkAtom[i]; } } dEnergyAtX = energy_force (0); neval++; double dGradDotD = calc_grad_dot_v_using_mpi_ (VEC_CG_D); double tau = compute_to_tr_ (dPP, dPdotD, dDD, dTrustRadius, true, dDHD, dPdotHD, dGradDotD); //---- MOVE THE POINT. if (nextra_global) { double * pGlobal = _daExtraGlobal[VEC_CG_P]; double * dGlobal = _daExtraGlobal[VEC_CG_D]; for (int i = 0; i < nextra_global; i++) { pGlobal[i] += tau * dGlobal[i]; } } - for (int i = 0; i < n3; i++) + for (int i = 0; i < nvec; i++) _daAVectors[VEC_CG_P][i] += tau * _daAVectors[VEC_CG_D][i]; if (nextra_atom) { for (int m = 0; m < nextra_atom; m++) { double * pAtom = _daExtraAtom[VEC_CG_P][m]; double * dAtom = _daExtraAtom[VEC_CG_D][m]; int n = extra_nlen[m]; for (int i = 0; i < n; i++) pAtom[i] += tau * dAtom[i]; } } nStepType = CGSTEP_NEGATIVE_CURVATURE; calc_plengths_using_mpi_ (dStepLength2, dStepLengthInf); return( true ); } //---- COMPUTE THE OPTIMAL STEP LENGTH BASED ON THE QUADRATIC CG MODEL. double dAlpha = dRR / dDHD; //---- MIGHT WANT TO ENABLE THIS TO DEBUG INTERNAL CG STEPS. //fprintf (_fpPrint, " alpha = %11.8f neval=%4d\n", dAlpha, neval); //---- p_i+1 = p_i + alpha_i d_i //---- (SAVE THE CURRENT p_i IN CASE THE STEP HAS TO BE SHORTENED.) if (nextra_global) { double * pGlobal = _daExtraGlobal[VEC_CG_P]; double * dGlobal = _daExtraGlobal[VEC_CG_D]; double * d1Global = _daExtraGlobal[VEC_DIF1]; for (int i = 0; i < nextra_global; i++) { d1Global[i] = pGlobal[i]; pGlobal[i] += dAlpha * dGlobal[i]; } } - for (int i = 0; i < n3; i++) { + for (int i = 0; i < nvec; i++) { _daAVectors[VEC_DIF1][i] = _daAVectors[VEC_CG_P][i]; _daAVectors[VEC_CG_P][i] += dAlpha * _daAVectors[VEC_CG_D][i]; } if (nextra_atom) { for (int m = 0; m < nextra_atom; m++) { double * pAtom = _daExtraAtom[VEC_CG_P][m]; double * dAtom = _daExtraAtom[VEC_CG_D][m]; double * d1Atom = _daExtraAtom[VEC_DIF1][m]; int n = extra_nlen[m]; for (int i = 0; i < n; i++) { d1Atom[i] = pAtom[i]; pAtom[i] += dAlpha * dAtom[i]; } } } //---- COMPUTE VECTOR PRODUCTS p_i+1^T p_i+1 AND p_i^T d_i. double dPnewDotPnew; double dPoldDotD; calc_ppnew_pdold_using_mpi_ (dPnewDotPnew, dPoldDotD); nStepType = CGSTEP_UNDETERMINED; //---- IF STEP LENGTH IS TOO LARGE, THEN REDUCE IT AND RETURN. double tau; if (step_exceeds_TR_ (dTrustRadius, dPP, dPoldDotD, dDD, tau)) { adjust_step_to_tau_ (tau); nStepType = CGSTEP_TO_TR; } if (step_exceeds_DMAX_()) { adjust_step_to_tau_ (0.0); nStepType = CGSTEP_TO_DMAX; } if ((nStepType == CGSTEP_TO_TR) || (nStepType == CGSTEP_TO_DMAX)) { calc_plengths_using_mpi_ (dStepLength2, dStepLengthInf); return( true ); } dStepLength2 = sqrt (dPnewDotPnew); //---- r_i+1 = r_i - alpha * H d_i if (nextra_global) { double * rGlobal = _daExtraGlobal[VEC_CG_R]; double * hdGlobal = _daExtraGlobal[VEC_CG_HD]; for (int i = 0; i < nextra_global; i++) rGlobal[i] -= dAlpha * hdGlobal[i]; } - for (int i = 0; i < n3; i++) + for (int i = 0; i < nvec; i++) _daAVectors[VEC_CG_R][i] -= dAlpha * _daAVectors[VEC_CG_HD][i]; if (nextra_atom) { for (int m = 0; m < nextra_atom; m++) { double * rAtom = _daExtraAtom[VEC_CG_R][m]; double * hdAtom = _daExtraAtom[VEC_CG_HD][m]; int n = extra_nlen[m]; for (int i = 0; i < n; i++) rAtom[i] -= dAlpha * hdAtom[i]; } } double dRnewDotRnew = calc_dot_prod_using_mpi_ (VEC_CG_R, VEC_CG_R); //---- IF RESIDUAL IS SMALL ENOUGH, THEN RETURN THE CURRENT STEP. if (sqrt (dRnewDotRnew) < dForceTol * dR0norm2) { nStepType = CGSTEP_NEWTON; calc_plengths_using_mpi_ (dStepLength2, dStepLengthInf); return( true ); } //---- beta = r_i+1^T r_i+1 / r_i^T r_i //---- d_i+1 = r_i+1 + beta d_i double dBeta = dRnewDotRnew / dRR; if (nextra_global) { double * rGlobal = _daExtraGlobal[VEC_CG_R]; double * dGlobal = _daExtraGlobal[VEC_CG_D]; for (int i = 0; i < nextra_global; i++) dGlobal[i] = rGlobal[i] + dBeta * dGlobal[i]; } - for (int i = 0; i < n3; i++) + for (int i = 0; i < nvec; i++) _daAVectors[VEC_CG_D][i] = _daAVectors[VEC_CG_R][i] + dBeta * _daAVectors[VEC_CG_D][i]; if (nextra_atom) { for (int m = 0; m < nextra_atom; m++) { double * rAtom = _daExtraAtom[VEC_CG_R][m]; double * dAtom = _daExtraAtom[VEC_CG_D][m]; int n = extra_nlen[m]; for (int i = 0; i < n; i++) dAtom[i] = rAtom[i] + dBeta * dAtom[i]; } } //---- CONTINUE THE LOOP. dRR = dRnewDotRnew; dPP = dPnewDotPnew; } nStepType = CGSTEP_MAX_INNER_ITERS; calc_plengths_using_mpi_ (dStepLength2, dStepLengthInf); return( true ); } /* ---------------------------------------------------------------------- Private method calc_xinf_using_mpi_ ------------------------------------------------------------------------- */ double MinHFTN::calc_xinf_using_mpi_(void) const { double dXInfLocal = 0.0; - for (int i = 0; i < n3; i++) - dXInfLocal = MAX (dXInfLocal, fabs (atom->x[0][i])); + for (int i = 0; i < nvec; i++) + dXInfLocal = MAX(dXInfLocal,fabs(xvec[i])); double dXInf; MPI_Allreduce (&dXInfLocal, &dXInf, 1, MPI_DOUBLE, MPI_MAX, world); if (nextra_atom) { for (int m = 0; m < nextra_atom; m++) { double * xatom = xextra_atom[m]; int n = extra_nlen[m]; double dXInfLocalExtra = 0.0; for (int i = 0; i < n; i++) dXInfLocalExtra = MAX (dXInfLocalExtra, fabs (xatom[i])); double dXInfExtra; MPI_Allreduce (&dXInfLocalExtra, &dXInfExtra, 1, MPI_DOUBLE, MPI_MAX, world); dXInf = MAX (dXInf, dXInfExtra); } } return( dXInf ); } /* ---------------------------------------------------------------------- Private method calc_dot_prod_using_mpi_ ------------------------------------------------------------------------- */ double MinHFTN::calc_dot_prod_using_mpi_(const int nIx1, const int nIx2) const { double dDotLocal = 0.0; - for (int i = 0; i < n3; i++) + for (int i = 0; i < nvec; i++) dDotLocal += _daAVectors[nIx1][i] * _daAVectors[nIx2][i]; if (nextra_atom) { for (int m = 0; m < nextra_atom; m++) { double * i1Atom = _daExtraAtom[nIx1][m]; double * i2Atom = _daExtraAtom[nIx2][m]; int n = extra_nlen[m]; for (int i = 0; i < n; i++) dDotLocal += i1Atom[i] * i2Atom[i]; } } double dDot; MPI_Allreduce (&dDotLocal, &dDot, 1, MPI_DOUBLE, MPI_SUM, world); if (nextra_global) { for (int i = 0; i < nextra_global; i++) { double * i1Global = _daExtraGlobal[nIx1]; double * i2Global = _daExtraGlobal[nIx2]; dDot += i1Global[i] * i2Global[i]; } } return( dDot ); } /* ---------------------------------------------------------------------- Private method calc_grad_dot_v_using_mpi_ ------------------------------------------------------------------------- */ double MinHFTN::calc_grad_dot_v_using_mpi_(const int nIx) const { //---- ASSUME THAT FORCES HAVE BEEN EVALUATED AT DESIRED ATOM POSITIONS. //---- REMEMBER THAT FORCES = -GRADIENT. double dGradDotVLocal = 0.0; - for (int i = 0; i < n3; i++) - dGradDotVLocal += - _daAVectors[nIx][i] * atom->f[0][i]; + for (int i = 0; i < nvec; i++) + dGradDotVLocal += - _daAVectors[nIx][i] * fvec[i]; if (nextra_atom) { - for (int m = 0; m < nextra_atom; m++) { + for (int m = 0; m < nextra_atom; m++) { double * iAtom = _daExtraAtom[nIx][m]; int n = extra_nlen[m]; - for (int i = 0; i < n; i++) + for (int i = 0; i < n; i++) dGradDotVLocal += - iAtom[i] * fextra_atom[m][i]; } } double dGradDotV; MPI_Allreduce (&dGradDotVLocal, &dGradDotV, 1, MPI_DOUBLE, MPI_SUM, world); if (nextra_global) { for (int i = 0; i < nextra_global; i++) { double * iGlobal = _daExtraGlobal[nIx]; dGradDotV += - iGlobal[i] * fextra[i]; } } return( dGradDotV ); } /* ---------------------------------------------------------------------- Private method calc_dhd_dd_using_mpi_ ------------------------------------------------------------------------- */ void MinHFTN::calc_dhd_dd_using_mpi_(double & dDHD, double & dDD) const { double dDHDLocal = 0.0; double dDDLocal = 0.0; - for (int i = 0; i < n3; i++) { + for (int i = 0; i < nvec; i++) { dDHDLocal += _daAVectors[VEC_CG_D][i] * _daAVectors[VEC_CG_HD][i]; dDDLocal += _daAVectors[VEC_CG_D][i] * _daAVectors[VEC_CG_D][i]; } if (nextra_atom) { for (int m = 0; m < nextra_atom; m++) { double * dAtom = _daExtraAtom[VEC_CG_D][m]; double * hdAtom = _daExtraAtom[VEC_CG_HD][m]; int n = extra_nlen[m]; for (int i = 0; i < n; i++) { dDHDLocal += dAtom[i] * hdAtom[i]; dDDLocal += dAtom[i] * dAtom[i]; } } } double daDotsLocal[2]; daDotsLocal[0] = dDHDLocal; daDotsLocal[1] = dDDLocal; double daDots[2]; MPI_Allreduce (daDotsLocal, daDots, 2, MPI_DOUBLE, MPI_SUM, world); if (nextra_global) { double * dGlobal = _daExtraGlobal[VEC_CG_D]; double * hdGlobal = _daExtraGlobal[VEC_CG_HD]; for (int i = 0; i < nextra_global; i++) { daDots[0] += dGlobal[i] * hdGlobal[i]; daDots[1] += dGlobal[i] * dGlobal[i]; } } dDHD = daDots[0]; dDD = daDots[1]; return; } /* ---------------------------------------------------------------------- Private method calc_ppnew_pdold_using_mpi_ ------------------------------------------------------------------------- */ void MinHFTN::calc_ppnew_pdold_using_mpi_(double & dPnewDotPnew, double & dPoldDotD) const { double dPnewDotPnewLocal = 0.0; double dPoldDotDLocal = 0.0; - for (int i = 0; i < n3; i++) { + for (int i = 0; i < nvec; i++) { dPnewDotPnewLocal += _daAVectors[VEC_CG_P][i] * _daAVectors[VEC_CG_P][i]; dPoldDotDLocal += _daAVectors[VEC_DIF1][i] * _daAVectors[VEC_CG_D][i]; } if (nextra_atom) { for (int m = 0; m < nextra_atom; m++) { double * dAtom = _daExtraAtom[VEC_CG_D][m]; double * pAtom = _daExtraAtom[VEC_CG_P][m]; double * d1Atom = _daExtraAtom[VEC_DIF1][m]; int n = extra_nlen[m]; for (int i = 0; i < n; i++) { dPnewDotPnewLocal += pAtom[i] * pAtom[i]; dPoldDotDLocal += d1Atom[i] * dAtom[i]; } } } double daDotsLocal[2]; daDotsLocal[0] = dPnewDotPnewLocal; daDotsLocal[1] = dPoldDotDLocal; double daDots[2]; MPI_Allreduce (daDotsLocal, daDots, 2, MPI_DOUBLE, MPI_SUM, world); if (nextra_global) { for (int i = 0; i < nextra_global; i++) { double * dGlobal = _daExtraGlobal[VEC_CG_D]; double * pGlobal = _daExtraGlobal[VEC_CG_P]; double * d1Global = _daExtraGlobal[VEC_DIF1]; daDots[0] += pGlobal[i] * pGlobal[i]; daDots[1] += d1Global[i] * dGlobal[i]; } } dPnewDotPnew = daDots[0]; dPoldDotD = daDots[1]; return; } /* ---------------------------------------------------------------------- Private method calc_plengths_using_mpi_ ------------------------------------------------------------------------- */ void MinHFTN::calc_plengths_using_mpi_(double & dStepLength2, double & dStepLengthInf) const { double dPPLocal = 0.0; double dPInfLocal = 0.0; - for (int i = 0; i < n3; i++) { + for (int i = 0; i < nvec; i++) { dPPLocal += _daAVectors[VEC_CG_P][i] * _daAVectors[VEC_CG_P][i]; dPInfLocal = MAX (dPInfLocal, fabs (_daAVectors[VEC_CG_P][i])); } if (nextra_atom) { for (int m = 0; m < nextra_atom; m++) { double * pAtom = _daExtraAtom[VEC_CG_P][m]; int n = extra_nlen[m]; for (int i = 0; i < n; i++) { dPPLocal += pAtom[i] * pAtom[i]; dPInfLocal = MAX (dPInfLocal, fabs (pAtom[i])); } } } double dPP; MPI_Allreduce (&dPPLocal, &dPP, 1, MPI_DOUBLE, MPI_SUM, world); double dPInf; MPI_Allreduce (&dPInfLocal, &dPInf, 1, MPI_DOUBLE, MPI_MAX, world); if (nextra_global) { for (int i = 0; i < nextra_global; i++) { double * pGlobal = _daExtraGlobal[VEC_CG_P]; dPP += pGlobal[i] * pGlobal[i]; dPInf = MAX (dPInf, fabs (pGlobal[i])); } } dStepLength2 = sqrt (dPP); dStepLengthInf = dPInf; return; } /* ---------------------------------------------------------------------- Private method step_exceeds_TR_ ------------------------------------------------------------------------- */ bool MinHFTN::step_exceeds_TR_(const double dTrustRadius, const double dPP, const double dPD, const double dDD, double & dTau) const { double dPnewNorm2; double dPnewNormInf; calc_plengths_using_mpi_ (dPnewNorm2, dPnewNormInf); if (dPnewNorm2 > dTrustRadius) { dTau = compute_to_tr_ (dPP, dPD, dDD, dTrustRadius, false, 0.0, 0.0, 0.0); return( true ); } //---- STEP LENGTH IS NOT TOO LONG. dTau = 0.0; return( false ); } /* ---------------------------------------------------------------------- Private method step_exceeds_DMAX_ Check that atoms do not move too far: for atom coordinates: limit is dmax for extra per-atom DOF: limit is extra_max[] for extra global DOF: limit is set by modify->max_alpha, which calls fix_box_relax->max_alpha ------------------------------------------------------------------------- */ bool MinHFTN::step_exceeds_DMAX_(void) const { double dAlpha = dmax * sqrt (_nNumUnknowns); double dPInfLocal = 0.0; - for (int i = 0; i < n3; i++) + for (int i = 0; i < nvec; i++) dPInfLocal = MAX (dPInfLocal, fabs (_daAVectors[VEC_CG_P][i])); double dPInf; MPI_Allreduce (&dPInfLocal, &dPInf, 1, MPI_DOUBLE, MPI_MAX, world); if (dPInf > dmax) return( true ); if (dPInf > MACHINE_EPS) dAlpha = MIN (dAlpha, dmax / dPInf); if (nextra_atom) { for (int m = 0; m < nextra_atom; m++) { double * pAtom = _daExtraAtom[VEC_CG_P][m]; dPInfLocal = 0.0; int n = extra_nlen[m]; for (int i = 0; i < n; i++) dPInfLocal = MAX (dPInfLocal, fabs (pAtom[i])); MPI_Allreduce (&dPInfLocal, &dPInf, 1, MPI_DOUBLE, MPI_MAX, world); if (dPInf > extra_max[m]) return( true ); if (dPInf > MACHINE_EPS) dAlpha = MIN (dAlpha, extra_max[m] / dPInf); } } if (nextra_global) { //---- IF THE MAXIMUM DISTANCE THAT THE GLOBAL BOX CONSTRAINT WILL //---- ALLOW IS SMALLER THAN THE PROPOSED DISTANCE, THEN THE STEP //---- IS TOO LONG. PROPOSED DISTANCE IS ESTIMATED BY |P|_INF. double dAlphaExtra = modify->max_alpha (_daExtraGlobal[VEC_CG_P]); if (dAlphaExtra < dAlpha) return( true ); } //---- STEP LENGTH IS NOT TOO LONG. return( false ); } /* ---------------------------------------------------------------------- Private method adjust_step_to_tau_ Adjust the step so that VEC_CG_P = VEC_DIF1 + tau * VEC_CG_D. ------------------------------------------------------------------------- */ void MinHFTN::adjust_step_to_tau_(const double tau) { if (nextra_global) { double * pGlobal = _daExtraGlobal[VEC_CG_P]; double * dGlobal = _daExtraGlobal[VEC_CG_D]; double * d1Global = _daExtraGlobal[VEC_DIF1]; for (int i = 0; i < nextra_global; i++) pGlobal[i] = d1Global[i] + (tau * dGlobal[i]); } - for (int i = 0; i < n3; i++) { + for (int i = 0; i < nvec; i++) { _daAVectors[VEC_CG_P][i] = _daAVectors[VEC_DIF1][i] + (tau * _daAVectors[VEC_CG_D][i]); } if (nextra_atom) { for (int m = 0; m < nextra_atom; m++) { double * pAtom = _daExtraAtom[VEC_CG_P][m]; double * dAtom = _daExtraAtom[VEC_CG_D][m]; double * d1Atom = _daExtraAtom[VEC_DIF1][m]; int n = extra_nlen[m]; for (int i = 0; i < n; i++) pAtom[i] = d1Atom[i] + (tau * dAtom[i]); } } return; } /* ---------------------------------------------------------------------- Private method compute_to_tr_ Return the value tau that solves || p + tau*d ||_2 = dTrustRadius If both roots are considered, the TR method chooses the one that minimizes grad^T (p + tau*d) + 0.5 (p + tau*d)^T H (p + tau*d) @param[in] dPP - p^T p @param[in] dPD - p^T d @param[in] dDD - d^T d @param[in] dTrustRadius - distance to match @param[in] bConsiderBothRoots - evaluate both roots, or return the positive @param[in] dDHD - d^T H d @param[in] dPdotHD - p^T H d @param[in] dGradDotD - grad(x_k)^T d ------------------------------------------------------------------------- */ double MinHFTN::compute_to_tr_(const double dPP, const double dPD, const double dDD, const double dTrustRadius, const bool bConsiderBothRoots, const double dDHD, const double dPdotHD, const double dGradDotD) const { //---- SOLVE A QUADRATIC EQUATION FOR TAU. //---- THE COEFFICIENTS ARE SUCH THAT THERE ARE ALWAYS TWO REAL ROOTS, //---- ONE POSITIVE AND ONE NEGATIVE. //---- CHECK FOR ERRONEOUS DATA. if ( (dDD <= 0.0) || (dPP < 0.0) || (dTrustRadius < 0.0) || (dTrustRadius * dTrustRadius < dPP) ) { printf ("HFTN internal error - bad data given to compute_to_tr_()\n"); return( 0.0 ); } double dTRsqrd = dTrustRadius * dTrustRadius; double dDiscr = (dPD * dPD) - (dDD * (dPP - dTRsqrd)); dDiscr = MAX (0.0, dDiscr); //-- SHOULD NEVER BE NEGATIVE dDiscr = sqrt (dDiscr); double dRootPos = (-dPD + dDiscr) / dDD; double dRootNeg = (-dPD - dDiscr) / dDD; if (bConsiderBothRoots == false) return( dRootPos ); //---- EVALUATE THE CG OBJECTIVE FUNCTION FOR EACH ROOT. double dTmpTerm = dGradDotD + dPdotHD; double dCgRedPos = (dRootPos * dTmpTerm) + (0.5 * dRootPos*dRootPos * dDHD); double dCgRedNeg = (dRootNeg * dTmpTerm) + (0.5 * dRootNeg*dRootNeg * dDHD); if ((-dCgRedPos) > (-dCgRedNeg)) return( dRootPos ); else return( dRootNeg ); } /* ---------------------------------------------------------------------- Private method evaluate_dir_der_ Compute the directional derivative using a finite difference approximation. This is equivalent to the Hessian times direction vector p. As a side effect, the method computes the energy and forces at x. On input these must be defined: atom->x - positions at x atom->f - ignored nIxDir - VEC_ index of the direction p nIxResult - ignored On output these are defined: atom->x - unchanged atom->f - forces evaluated at x, only if bEvaluateAtX is true nIxDir - unchanged nIxResult - directional derivative Hp During processing these are modified: VEC_DIF1 VEC_DIF2 @param[in] bUseForwardDiffs - if true use forward difference approximation, else use central difference @param[in] nIxDir - VEC_ index of the direction @param[in] nIxResult - VEC_ index to place the result (it is acceptable for nIxDir = nIxResult) @param[in] bEvaluateAtX - if true, then evaluate at x before returning @param[out] dNewEnergy - energy at x, if bEvaluateAtX is true @param[out] dNewForce2 - |F|_2 at x, if bEvaluateAtX is true ------------------------------------------------------------------------- */ void MinHFTN::evaluate_dir_der_(const bool bUseForwardDiffs, const int nIxDir, const int nIxResult, const bool bEvaluateAtX, double & dNewEnergy) { //---- COMPUTE THE MAGNITUDE OF THE DIRECTION VECTOR: |p|_2. double dDirNorm2SqrdLocal = 0.0; - for (int i = 0; i < n3; i++) + for (int i = 0; i < nvec; i++) dDirNorm2SqrdLocal += _daAVectors[nIxDir][i] * _daAVectors[nIxDir][i]; if (nextra_atom) { for (int m = 0; m < nextra_atom; m++) { double * iAtom = _daExtraAtom[nIxDir][m]; int n = extra_nlen[m]; for (int i = 0; i < n; i++) dDirNorm2SqrdLocal += iAtom[i] * iAtom[i]; } } double dDirNorm2Sqrd = 0.0; MPI_Allreduce (&dDirNorm2SqrdLocal, &dDirNorm2Sqrd, 1, MPI_DOUBLE, MPI_SUM, world); if (nextra_global) { for (int i = 0; i < nextra_global; i++) { double * iGlobal = _daExtraGlobal[nIxDir]; dDirNorm2Sqrd += iGlobal[i] * iGlobal[i]; } } double dDirNorm2 = sqrt (dDirNorm2Sqrd); //---- IF THE STEP IS TOO SMALL, RETURN ZERO FOR THE DERIVATIVE. if (dDirNorm2 == 0.0) { - for (int i = 0; i < n3; i++) + for (int i = 0; i < nvec; i++) _daAVectors[nIxResult][i] = 0.0; if (nextra_atom) { for (int m = 0; m < nextra_atom; m++) { double * iAtom = _daExtraAtom[nIxDir][m]; int n = extra_nlen[m]; for (int i = 0; i < n; i++) iAtom[i] = 0; } } if (nextra_global) { for (int i = 0; i < nextra_global; i++) _daExtraGlobal[nIxDir][i] = 0.0; } if (bEvaluateAtX) { dNewEnergy = energy_force (0); neval++; } return; } //---- FORWARD DIFFERENCES ARE LESS ACCURATE THAN CENTRAL DIFFERENCES, //---- BUT REQUIRE ONLY 2 ENERGY+FORCE EVALUATIONS VERSUS 3 EVALUATIONS. //---- STORAGE REQUIREMENTS ARE THE SAME. if (bUseForwardDiffs) { //---- EQUATION IS FROM THE OLD LAMMPS VERSION, SAND98-8201. double dEps = 2.0 * sqrt (1000.0 * MACHINE_EPS) / dDirNorm2; //---- SAVE A COPY OF x. fix_minimize->store_box(); - for (int i = 0; i < n3; i++) - _daAVectors[VEC_DIF1][i] = atom->x[0][i]; + for (int i = 0; i < nvec; i++) + _daAVectors[VEC_DIF1][i] = xvec[i]; if (nextra_atom) { for (int m = 0; m < nextra_atom; m++) { double * xatom = xextra_atom[m]; double * d1Atom = _daExtraAtom[VEC_DIF1][m]; int n = extra_nlen[m]; for (int i = 0; i < n; i++) d1Atom[i] = xatom[i]; } } if (nextra_global) { modify->min_pushstore(); modify->min_store(); } //---- EVALUATE FORCES AT x + eps*p. if (nextra_global) modify->min_step (dEps, _daExtraGlobal[nIxDir]); - for (int i = 0; i < n3; i++) - atom->x[0][i] += dEps * _daAVectors[nIxDir][i]; + for (int i = 0; i < nvec; i++) + xvec[i] += dEps * _daAVectors[nIxDir][i]; if (nextra_atom) { for (int m = 0; m < nextra_atom; m++) { double * xatom = xextra_atom[m]; double * iAtom = _daExtraAtom[nIxDir][m]; int n = extra_nlen[m]; for (int i = 0; i < n; i++) xatom[i] += dEps * iAtom[i]; } } energy_force (0); neval++; //---- STORE THE FORCE IN DIF2. if (nextra_global) { for (int i = 0; i < nextra_global; i++) _daExtraGlobal[VEC_DIF2][i] = fextra[i]; } - for (int i = 0; i < n3; i++) - _daAVectors[VEC_DIF2][i] = atom->f[0][i]; + for (int i = 0; i < nvec; i++) + _daAVectors[VEC_DIF2][i] = fvec[i]; if (nextra_atom) { for (int m = 0; m < nextra_atom; m++) { double * d2Atom = _daExtraAtom[VEC_DIF2][m]; int n = extra_nlen[m]; for (int i = 0; i < n; i++) d2Atom[i] = fextra_atom[m][i]; } } //---- MOVE BACK TO x AND EVALUATE FORCES. if (nextra_global) { modify->min_step (0.0, _daExtraGlobal[VEC_DIF1]); modify->min_popstore(); } - for (int i = 0; i < n3; i++) - atom->x[0][i] = _daAVectors[VEC_DIF1][i]; + for (int i = 0; i < nvec; i++) + xvec[i] = _daAVectors[VEC_DIF1][i]; if (nextra_atom) { for (int m = 0; m < nextra_atom; m++) { double * xatom = xextra_atom[m]; double * d1Atom = _daExtraAtom[VEC_DIF1][m]; int n = extra_nlen[m]; for (int i = 0; i < n; i++) xatom[i] += d1Atom[i]; } } dNewEnergy = energy_force (0); neval++; //---- COMPUTE THE DIFFERENCE VECTOR: [grad(x + eps*p) - grad(x)] / eps. //---- REMEMBER THAT FORCES = -GRADIENT. - for (int i = 0; i < n3; i++) - _daAVectors[nIxResult][i] - = (atom->f[0][i] - _daAVectors[VEC_DIF2][i]) / dEps; + for (int i = 0; i < nvec; i++) + _daAVectors[nIxResult][i] = (fvec[i] - _daAVectors[VEC_DIF2][i]) / dEps; if (nextra_atom) { for (int m = 0; m < nextra_atom; m++) { double * iAtom = _daExtraAtom[nIxResult][m]; double * d2Atom = _daExtraAtom[VEC_DIF2][m]; int n = extra_nlen[m]; for (int i = 0; i < n; i++) iAtom[i] = (fextra_atom[m][i] - d2Atom[i]) / dEps; } } if (nextra_global) { for (int i = 0; i < nextra_global; i++) _daExtraGlobal[nIxResult][i] = (fextra[i] - _daExtraGlobal[VEC_DIF2][i]) / dEps; } } else { //-- bUseForwardDiffs == false //---- EQUATION IS FROM THE OLD LAMMPS VERSION, SAND98-8201. double dEps = pow (3000.0 * MACHINE_EPS, 0.33333333) / dDirNorm2; //---- SAVE A COPY OF x. fix_minimize->store_box(); - for (int i = 0; i < n3; i++) - _daAVectors[VEC_DIF1][i] = atom->x[0][i]; + for (int i = 0; i < nvec; i++) + _daAVectors[VEC_DIF1][i] = xvec[i]; if (nextra_atom) { for (int m = 0; m < nextra_atom; m++) { double * xatom = xextra_atom[m]; double * d1Atom = _daExtraAtom[VEC_DIF1][m]; int n = extra_nlen[m]; for (int i = 0; i < n; i++) d1Atom[i] = xatom[i]; } } if (nextra_global) { modify->min_pushstore(); modify->min_store(); } //---- EVALUATE FORCES AT x + eps*p. if (nextra_global) modify->min_step (dEps, _daExtraGlobal[nIxDir]); - for (int i = 0; i < n3; i++) - atom->x[0][i] += dEps * _daAVectors[nIxDir][i]; + for (int i = 0; i < nvec; i++) + xvec[i] += dEps * _daAVectors[nIxDir][i]; if (nextra_atom) { for (int m = 0; m < nextra_atom; m++) { double * xatom = xextra_atom[m]; double * iAtom = _daExtraAtom[nIxDir][m]; int n = extra_nlen[m]; for (int i = 0; i < n; i++) xatom[i] += dEps * iAtom[i]; } } energy_force (0); neval++; //---- STORE THE FORCE IN DIF2. if (nextra_global) { for (int i = 0; i < nextra_global; i++) _daExtraGlobal[VEC_DIF2][i] = fextra[i]; } - for (int i = 0; i < n3; i++) - _daAVectors[VEC_DIF2][i] = atom->f[0][i]; + for (int i = 0; i < nvec; i++) + _daAVectors[VEC_DIF2][i] = fvec[i]; if (nextra_atom) { for (int m = 0; m < nextra_atom; m++) { double * d2Atom = _daExtraAtom[VEC_DIF2][m]; int n = extra_nlen[m]; for (int i = 0; i < n; i++) d2Atom[i] = fextra_atom[m][i]; } } //---- EVALUATE FORCES AT x - eps*p. if (nextra_global) modify->min_step (-dEps, _daExtraGlobal[nIxDir]); - for (int i = 0; i < n3; i++) - atom->x[0][i] = _daAVectors[VEC_DIF1][i] + for (int i = 0; i < nvec; i++) + xvec[i] = _daAVectors[VEC_DIF1][i] - dEps * _daAVectors[nIxDir][i]; if (nextra_atom) { for (int m = 0; m < nextra_atom; m++) { double * xatom = xextra_atom[m]; double * iAtom = _daExtraAtom[nIxDir][m]; double * d1Atom = _daExtraAtom[VEC_DIF1][m]; int n = extra_nlen[m]; for (int i = 0; i < n; i++) xatom[i] = d1Atom[i] - dEps * iAtom[i]; } } energy_force (0); neval++; //---- COMPUTE THE DIFFERENCE VECTOR: //---- [grad(x + eps*p) - grad(x - eps*p)] / 2*eps. //---- REMEMBER THAT FORCES = -GRADIENT. if (nextra_global) { double * iGlobal = _daExtraGlobal[nIxResult]; double * d2Global = _daExtraGlobal[VEC_DIF2]; for (int i = 0; i < nextra_global; i++) iGlobal[i] = (fextra[i] - d2Global[i]) / (2.0 + dEps); } - for (int i = 0; i < n3; i++) - _daAVectors[nIxResult][i] - = (atom->f[0][i] - _daAVectors[VEC_DIF2][i]) / (2.0 * dEps); + for (int i = 0; i < nvec; i++) + _daAVectors[nIxResult][i] = + (fvec[i] - _daAVectors[VEC_DIF2][i]) / (2.0 * dEps); if (nextra_atom) { for (int m = 0; m < nextra_atom; m++) { double * iAtom = _daExtraAtom[nIxResult][m]; double * d2Atom = _daExtraAtom[VEC_DIF2][m]; int n = extra_nlen[m]; for (int i = 0; i < n; i++) iAtom[i] = (fextra_atom[m][i] - d2Atom[i]) / (2.0 + dEps); } } if (bEvaluateAtX) { //---- EVALUATE FORCES AT x. if (nextra_global) { modify->min_step (0.0, _daExtraGlobal[VEC_DIF1]); modify->min_popstore(); } - for (int i = 0; i < n3; i++) - atom->x[0][i] = _daAVectors[VEC_DIF1][i]; + for (int i = 0; i < nvec; i++) + xvec[i] = _daAVectors[VEC_DIF1][i]; if (nextra_atom) { for (int m = 0; m < nextra_atom; m++) { double * xatom = xextra_atom[m]; double * d1Atom = _daExtraAtom[VEC_DIF1][m]; int n = extra_nlen[m]; for (int i = 0; i < n; i++) xatom[i] = d1Atom[i]; } } dNewEnergy = energy_force (0); neval++; } } return; } /* ---------------------------------------------------------------------- Private method open_hftn_print_file_ ------------------------------------------------------------------------- */ void MinHFTN::open_hftn_print_file_(void) { int nMyRank; MPI_Comm_rank (world, &nMyRank); char szTmp[50]; sprintf (szTmp, "progress_MinHFTN_%d.txt", nMyRank); _fpPrint = fopen (szTmp, "w"); if (_fpPrint == NULL) { printf ("*** MinHFTN cannot open file '%s'\n", szTmp); printf ("*** continuing...\n"); return; } fprintf (_fpPrint, " Iter Evals Energy |F|_2" " Step TR used |step|_2 ared pred\n"); return; } /* ---------------------------------------------------------------------- Private method hftn_print_line_ Step types: 1 - Nw (inner iteration converged like a Newton step) 2 - TR (inner iteration reached the trust region boundary) 3 - Neg (inner iteration ended with negative curvature) ------------------------------------------------------------------------- */ void MinHFTN::hftn_print_line_(const bool bIsStepAccepted, const int nIteration, const int nTotalEvals, const double dEnergy, const double dForce2, const int nStepType, const double dTrustRadius, const double dStepLength2, const double dActualRed, const double dPredictedRed) const { const char sFormat1[] = " %4d %5d %14.8f %11.5e\n"; const char sFormatA[] = " %4d %5d %14.8f %11.5e %3s %9.3e %8.2e %10.3e %10.3e\n"; const char sFormatR[] = "r %4d %5d %14.8f %11.5e %3s %9.3e %8.2e %10.3e %10.3e\n"; if (_fpPrint == NULL) return; char sStepType[4]; if (nStepType == NO_CGSTEP_BECAUSE_F_TOL_SATISFIED) strcpy (sStepType, " - "); else if (nStepType == CGSTEP_NEWTON) strcpy (sStepType, "Nw "); else if (nStepType == CGSTEP_TO_TR) strcpy (sStepType, "TR "); else if (nStepType == CGSTEP_TO_DMAX) strcpy (sStepType, "dmx"); else if (nStepType == CGSTEP_NEGATIVE_CURVATURE) strcpy (sStepType, "Neg"); else if (nStepType == CGSTEP_MAX_INNER_ITERS) strcpy (sStepType, "its"); else strcpy (sStepType, "???"); if (nIteration == -1) { fprintf (_fpPrint, sFormat1, 0, nTotalEvals, dEnergy, dForce2); } else { if (bIsStepAccepted) fprintf (_fpPrint, sFormatA, nIteration, nTotalEvals, dEnergy, dForce2, sStepType, dTrustRadius, dStepLength2, dActualRed, dPredictedRed); else fprintf (_fpPrint, sFormatR, nIteration, nTotalEvals, dEnergy, dForce2, sStepType, dTrustRadius, dStepLength2, dActualRed, dPredictedRed); } fflush (_fpPrint); return; } /* ---------------------------------------------------------------------- Private method close_hftn_print_file_ ------------------------------------------------------------------------- */ void MinHFTN::close_hftn_print_file_(void) { if (_fpPrint != NULL) fclose (_fpPrint); return; } diff --git a/src/min_linesearch.cpp b/src/min_linesearch.cpp index 6eb96cbec..933e606e6 100644 --- a/src/min_linesearch.cpp +++ b/src/min_linesearch.cpp @@ -1,546 +1,546 @@ /* ---------------------------------------------------------------------- LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator http://lammps.sandia.gov, Sandia National Laboratories Steve Plimpton, sjplimp@sandia.gov Copyright (2003) Sandia Corporation. Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains certain rights in this software. This software is distributed under the GNU General Public License. See the README file in the top-level LAMMPS directory. ------------------------------------------------------------------------- */ /* ---------------------------------------------------------------------- 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 "min_linesearch.h" #include "atom.h" #include "update.h" #include "neighbor.h" #include "domain.h" #include "modify.h" #include "fix_minimize.h" #include "pair.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 IDEAL_TOL 1.0e-8 #define QUADRATIC_TOL 0.1 #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) /* ---------------------------------------------------------------------- */ MinLineSearch::MinLineSearch(LAMMPS *lmp) : Min(lmp) { gextra = hextra = NULL; x0extra_atom = gextra_atom = hextra_atom = NULL; } /* ---------------------------------------------------------------------- */ MinLineSearch::~MinLineSearch() { delete [] gextra; delete [] hextra; delete [] x0extra_atom; delete [] gextra_atom; delete [] hextra_atom; } /* ---------------------------------------------------------------------- */ void MinLineSearch::init_style() { if (linestyle == 0) linemin = &MinLineSearch::linemin_backtrack; else if (linestyle == 1) linemin = &MinLineSearch::linemin_quadratic; delete [] gextra; delete [] hextra; gextra = hextra = NULL; delete [] x0extra_atom; delete [] gextra_atom; delete [] hextra_atom; x0extra_atom = gextra_atom = hextra_atom = NULL; } /* ---------------------------------------------------------------------- */ void MinLineSearch::setup_style() { // memory for x0,g,h for atomic dof fix_minimize->add_vector(3); fix_minimize->add_vector(3); fix_minimize->add_vector(3); // memory for g,h for extra global dof, fix stores x0 if (nextra_global) { gextra = new double[nextra_global]; hextra = new double[nextra_global]; } // memory for x0,g,h for extra per-atom dof if (nextra_atom) { x0extra_atom = new double*[nextra_atom]; gextra_atom = new double*[nextra_atom]; hextra_atom = new double*[nextra_atom]; for (int m = 0; m < nextra_atom; m++) { fix_minimize->add_vector(extra_peratom[m]); fix_minimize->add_vector(extra_peratom[m]); fix_minimize->add_vector(extra_peratom[m]); } } } /* ---------------------------------------------------------------------- set current vector lengths and pointers called after atoms have migrated ------------------------------------------------------------------------- */ void MinLineSearch::reset_vectors() { // atomic dof - n3 = 3 * atom->nlocal; - if (n3) x = atom->x[0]; - if (n3) f = atom->f[0]; + nvec = 3 * atom->nlocal; + if (nvec) xvec = atom->x[0]; + if (nvec) fvec = atom->f[0]; x0 = fix_minimize->request_vector(0); g = fix_minimize->request_vector(1); h = fix_minimize->request_vector(2); // extra per-atom dof if (nextra_atom) { int n = 3; for (int m = 0; m < nextra_atom; m++) { extra_nlen[m] = extra_peratom[m] * atom->nlocal; requestor[m]->min_pointers(&xextra_atom[m],&fextra_atom[m]); x0extra_atom[m] = fix_minimize->request_vector(n++); gextra_atom[m] = fix_minimize->request_vector(n++); hextra_atom[m] = fix_minimize->request_vector(n++); } } } /* ---------------------------------------------------------------------- line minimization methods find minimum-energy starting at x along h direction input args: eoriginal = energy at initial x input extra: n,x,x0,f,h for atomic, extra global, extra per-atom dof output args: return 0 if successful move, non-zero alpha return non-zero if failed alpha = distance moved along h for x at min eng config nfunc = updated counter of eng/force function evals output extra: if fail, energy_force() of original x if succeed, energy_force() at x + alpha*h atom->x = coords at new configuration atom->f = force at new configuration ecurrent = energy of new configuration ------------------------------------------------------------------------- */ /* ---------------------------------------------------------------------- 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 MinLineSearch::linemin_backtrack(double eoriginal, double &alpha, int &nfunc) { int i,m,n; double fdothall,fdothme,hme,hmax,hmaxall; double de_ideal,de; double *xatom,*x0atom,*fatom,*hatom; // fdothall = projection of search dir along downhill gradient // if search direction is not downhill, exit with error fdothme = 0.0; - for (i = 0; i < n3; i++) fdothme += f[i]*h[i]; + for (i = 0; i < nvec; i++) fdothme += fvec[i]*h[i]; if (nextra_atom) for (m = 0; m < nextra_atom; m++) { fatom = fextra_atom[m]; hatom = hextra_atom[m]; n = extra_nlen[m]; for (i = 0; i < n; i++) fdothme += fatom[i]*hatom[i]; } MPI_Allreduce(&fdothme,&fdothall,1,MPI_DOUBLE,MPI_SUM,world); if (nextra_global) for (i = 0; i < nextra_global; i++) fdothall += fextra[i]*hextra[i]; if (output->thermo->normflag) fdothall /= atom->natoms; if (fdothall <= 0.0) return DOWNHILL; // set alpha so no dof is changed by more than max allowed amount // for atom coords, max amount = dmax // for extra per-atom dof, max amount = extra_max[] // for extra global dof, max amount is set by fix // also insure alpha <= ALPHA_MAX // else will have to backtrack from huge value when forces are tiny // if all search dir components are already 0.0, exit with error hme = 0.0; - for (i = 0; i < n3; i++) hme = MAX(hme,fabs(h[i])); + for (i = 0; i < nvec; i++) hme = MAX(hme,fabs(h[i])); MPI_Allreduce(&hme,&hmaxall,1,MPI_DOUBLE,MPI_MAX,world); alpha = MIN(ALPHA_MAX,dmax/hmaxall); if (nextra_atom) for (m = 0; m < nextra_atom; m++) { hme = 0.0; fatom = fextra_atom[m]; n = extra_nlen[m]; for (i = 0; i < n; i++) hme = MAX(hme,fabs(hatom[i])); MPI_Allreduce(&hme,&hmax,1,MPI_DOUBLE,MPI_MAX,world); alpha = MIN(alpha,extra_max[m]/hmax); hmaxall = MAX(hmaxall,hmax); } if (nextra_global) { double alpha_extra = modify->max_alpha(hextra); alpha = MIN(alpha,alpha_extra); for (i = 0; i < nextra_global; i++) hmaxall = MAX(hmaxall,fabs(hextra[i])); } if (hmaxall == 0.0) return ZEROFORCE; // store box and values of all dof at start of linesearch fix_minimize->store_box(); - for (i = 0; i < n3; i++) x0[i] = x[i]; + for (i = 0; i < nvec; i++) x0[i] = xvec[i]; if (nextra_atom) for (m = 0; m < nextra_atom; m++) { xatom = xextra_atom[m]; x0atom = x0extra_atom[m]; n = extra_nlen[m]; for (i = 0; i < n; i++) x0atom[i] = xatom[i]; } if (nextra_global) modify->min_store(); // Important diagnostic: test the gradient against energy // double etmp; // double alphatmp = alphamax*1.0e-4; // etmp = alpha_step(alphatmp,1,nfunc); // printf("alpha = %g dele = %g dele_force = %g err = %g\n", // alphatmp,etmp-eoriginal,-alphatmp*fdothall, // etmp-eoriginal+alphatmp*fdothall); // alpha_step(0.0,1,nfunc); // backtrack with alpha until energy decrease is sufficient while (1) { ecurrent = alpha_step(alpha,1,nfunc); // if energy change is better than ideal, exit with success de_ideal = -BACKTRACK_SLOPE*alpha*fdothall; de = ecurrent - eoriginal; if (de <= de_ideal) { if (nextra_global) { int itmp = modify->min_reset_ref(); if (itmp) ecurrent = energy_force(1); } 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) { ecurrent = alpha_step(0.0,0,nfunc); return ZEROALPHA; } } } /* ---------------------------------------------------------------------- // linemin: quadratic line search (adapted from Dennis and Schnabel) // The objective function is approximated by a quadratic // function in alpha, for sufficiently small alpha. // This idea is the same as that used in the well-known secant // method. However, since the change in the objective function // (difference of two finite numbers) is not known as accurately // as the gradient (which is close to zero), all the expressions // are written in terms of gradients. In this way, we can converge // the LAMMPS forces much closer to zero. // // We know E,Eprev,fh,fhprev. The Taylor series about alpha_prev // truncated at the quadratic term is: // // E = Eprev - del_alpha*fhprev + (1/2)del_alpha^2*Hprev // // and // // fh = fhprev - del_alpha*Hprev // // where del_alpha = alpha-alpha_prev // // We solve these two equations for Hprev and E=Esolve, giving: // // Esolve = Eprev - del_alpha*(f+fprev)/2 // // We define relerr to be: // // relerr = |(Esolve-E)/Eprev| // = |1.0 - (0.5*del_alpha*(f+fprev)+E)/Eprev| // // If this is accurate to within a reasonable tolerance, then // we go ahead and use a secant step to fh = 0: // // alpha0 = alpha - (alpha-alphaprev)*fh/delfh; // ------------------------------------------------------------------------- */ int MinLineSearch::linemin_quadratic(double eoriginal, double &alpha, int &nfunc) { int i,m,n; double fdothall,fdothme,hme,hmax,hmaxall; double de_ideal,de; double delfh,engprev,relerr,alphaprev,fhprev,ff,fh,alpha0,fh0,ff0; double dot[2],dotall[2]; double *xatom,*x0atom,*fatom,*hatom; double alphamax; // fdothall = projection of search dir along downhill gradient // if search direction is not downhill, exit with error fdothme = 0.0; - for (i = 0; i < n3; i++) fdothme += f[i]*h[i]; + for (i = 0; i < nvec; i++) fdothme += fvec[i]*h[i]; if (nextra_atom) for (m = 0; m < nextra_atom; m++) { fatom = fextra_atom[m]; hatom = hextra_atom[m]; n = extra_nlen[m]; for (i = 0; i < n; i++) fdothme += fatom[i]*hatom[i]; } MPI_Allreduce(&fdothme,&fdothall,1,MPI_DOUBLE,MPI_SUM,world); if (nextra_global) for (i = 0; i < nextra_global; i++) fdothall += fextra[i]*hextra[i]; if (output->thermo->normflag) fdothall /= atom->natoms; if (fdothall <= 0.0) return DOWNHILL; // set alphamax so no dof is changed by more than max allowed amount // for atom coords, max amount = dmax // for extra per-atom dof, max amount = extra_max[] // for extra global dof, max amount is set by fix // also insure alphamax <= ALPHA_MAX // else will have to backtrack from huge value when forces are tiny // if all search dir components are already 0.0, exit with error hme = 0.0; - for (i = 0; i < n3; i++) hme = MAX(hme,fabs(h[i])); + for (i = 0; i < nvec; i++) hme = MAX(hme,fabs(h[i])); MPI_Allreduce(&hme,&hmaxall,1,MPI_DOUBLE,MPI_MAX,world); alphamax = MIN(ALPHA_MAX,dmax/hmaxall); if (nextra_atom) for (m = 0; m < nextra_atom; m++) { hme = 0.0; fatom = fextra_atom[m]; n = extra_nlen[m]; for (i = 0; i < n; i++) hme = MAX(hme,fabs(hatom[i])); MPI_Allreduce(&hme,&hmax,1,MPI_DOUBLE,MPI_MAX,world); alphamax = MIN(alphamax,extra_max[m]/hmax); hmaxall = MAX(hmaxall,hmax); } if (nextra_global) { double alpha_extra = modify->max_alpha(hextra); alphamax = MIN(alphamax,alpha_extra); for (i = 0; i < nextra_global; i++) hmaxall = MAX(hmaxall,fabs(hextra[i])); } if (hmaxall == 0.0) return ZEROFORCE; // store box and values of all dof at start of linesearch fix_minimize->store_box(); - for (i = 0; i < n3; i++) x0[i] = x[i]; + for (i = 0; i < nvec; i++) x0[i] = xvec[i]; if (nextra_atom) for (m = 0; m < nextra_atom; m++) { xatom = xextra_atom[m]; x0atom = x0extra_atom[m]; n = extra_nlen[m]; for (i = 0; i < n; i++) x0atom[i] = xatom[i]; } if (nextra_global) modify->min_store(); // backtrack with alpha until energy decrease is sufficient // or until get to small energy change, then perform quadratic projection alpha = alphamax; fhprev = fdothall; engprev = eoriginal; alphaprev = 0.0; // Important diagnostic: test the gradient against energy // double etmp; // double alphatmp = alphamax*1.0e-4; // etmp = alpha_step(alphatmp,1,nfunc); // printf("alpha = %g dele = %g dele_force = %g err = %g\n", // alphatmp,etmp-eoriginal,-alphatmp*fdothall, // etmp-eoriginal+alphatmp*fdothall); // alpha_step(0.0,1,nfunc); while (1) { ecurrent = alpha_step(alpha,1,nfunc); // compute new fh, alpha, delfh dot[0] = dot[1] = 0.0; - for (i = 0; i < n3; i++) { - dot[0] += f[i]*f[i]; - dot[1] += f[i]*h[i]; + for (i = 0; i < nvec; i++) { + dot[0] += fvec[i]*fvec[i]; + dot[1] += fvec[i]*h[i]; } if (nextra_atom) for (m = 0; m < nextra_atom; m++) { xatom = xextra_atom[m]; hatom = hextra_atom[m]; n = extra_nlen[m]; for (i = 0; i < n; i++) { dot[0] += fatom[i]*fatom[i]; dot[1] += fatom[i]*hatom[i]; } } MPI_Allreduce(dot,dotall,2,MPI_DOUBLE,MPI_SUM,world); if (nextra_global) { for (i = 0; i < nextra_global; i++) { dotall[0] += fextra[i]*fextra[i]; dotall[1] += fextra[i]*hextra[i]; } } ff = dotall[0]; fh = dotall[1]; if (output->thermo->normflag) { ff /= atom->natoms; fh /= atom->natoms; } 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) { ecurrent = alpha_step(0.0,0,nfunc); return ZEROQUAD; } // Check if ready for quadratic projection, equivalent to secant method // alpha0 = projected alpha relerr = fabs(1.0-(0.5*(alpha-alphaprev)*(fh+fhprev)+ecurrent)/engprev); alpha0 = alpha - (alpha-alphaprev)*fh/delfh; if (relerr <= QUADRATIC_TOL && alpha0 > 0.0 && alpha0 < alphamax) { ecurrent = alpha_step(alpha0,1,nfunc); if (ecurrent < eoriginal) { if (nextra_global) { int itmp = modify->min_reset_ref(); if (itmp) ecurrent = energy_force(1); } return 0; } } // if backtracking energy change is better than ideal, exit with success de_ideal = -BACKTRACK_SLOPE*alpha*fdothall; de = ecurrent - eoriginal; if (de <= de_ideal) { if (nextra_global) { int itmp = modify->min_reset_ref(); if (itmp) ecurrent = energy_force(1); } return 0; } // save previous state fhprev = fh; engprev = ecurrent; 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) { ecurrent = alpha_step(0.0,0,nfunc); return ZEROALPHA; } } } /* ---------------------------------------------------------------------- */ double MinLineSearch::alpha_step(double alpha, int resetflag, int &nfunc) { int i,n,m; double *xatom,*x0atom,*hatom; // reset to starting point if (nextra_global) modify->min_step(0.0,hextra); - for (i = 0; i < n3; i++) x[i] = x0[i]; + for (i = 0; i < nvec; i++) xvec[i] = x0[i]; if (nextra_atom) for (m = 0; m < nextra_atom; m++) { xatom = xextra_atom[m]; x0atom = x0extra_atom[m]; n = extra_nlen[m]; for (i = 0; i < n; i++) xatom[i] = x0atom[i]; } // step forward along h if (alpha > 0.0) { if (nextra_global) modify->min_step(alpha,hextra); - for (i = 0; i < n3; i++) x[i] += alpha*h[i]; + for (i = 0; i < nvec; i++) xvec[i] += alpha*h[i]; if (nextra_atom) for (m = 0; m < nextra_atom; m++) { xatom = xextra_atom[m]; hatom = hextra_atom[m]; n = extra_nlen[m]; for (i = 0; i < n; i++) xatom[i] += alpha*hatom[i]; } } // compute and return new energy nfunc++; return energy_force(resetflag); } diff --git a/src/min_sd.cpp b/src/min_sd.cpp index 99d31cab5..3d0e2c0cb 100644 --- a/src/min_sd.cpp +++ b/src/min_sd.cpp @@ -1,109 +1,109 @@ /* ---------------------------------------------------------------------- LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator http://lammps.sandia.gov, Sandia National Laboratories Steve Plimpton, sjplimp@sandia.gov Copyright (2003) Sandia Corporation. Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains certain rights in this software. This software is distributed under the GNU General Public License. See the README file in the top-level LAMMPS directory. ------------------------------------------------------------------------- */ #include "math.h" #include "mpi.h" #include "min_sd.h" #include "atom.h" #include "update.h" #include "output.h" #include "timer.h" using namespace LAMMPS_NS; // EPS_ENERGY = minimum normalization for energy tolerance #define EPS_ENERGY 1.0e-8 // same as in other min classes enum{MAXITER,MAXEVAL,ETOL,FTOL,DOWNHILL,ZEROALPHA,ZEROFORCE,ZEROQUAD}; /* ---------------------------------------------------------------------- */ MinSD::MinSD(LAMMPS *lmp) : MinLineSearch(lmp) {} /* ---------------------------------------------------------------------- minimization via steepest descent ------------------------------------------------------------------------- */ int MinSD::iterate(int niter_max) { int i,m,n,fail,ntimestep; double fdotf; double *fatom,*hatom; // initialize working vectors - for (i = 0; i < n3; i++) h[i] = f[i]; + for (i = 0; i < nvec; i++) h[i] = fvec[i]; if (nextra_atom) for (m = 0; m < nextra_atom; m++) { fatom = fextra_atom[m]; hatom = hextra_atom[m]; n = extra_nlen[m]; for (i = 0; i < n; i++) hatom[i] = fatom[i]; } if (nextra_global) for (i = 0; i < nextra_global; i++) hextra[i] = fextra[i]; neval = 0; for (niter = 0; niter < niter_max; niter++) { ntimestep = ++update->ntimestep; // line minimization along h from current position x // h = downhill gradient direction eprevious = ecurrent; fail = (this->*linemin)(ecurrent,alpha_final,neval); if (fail) return fail; // function evaluation criterion if (neval >= update->max_eval) return MAXEVAL; // energy tolerance criterion if (fabs(ecurrent-eprevious) < update->etol * 0.5*(fabs(ecurrent) + fabs(eprevious) + EPS_ENERGY)) return ETOL; // force tolerance criterion fdotf = fnorm_sqr(); if (fdotf < update->ftol*update->ftol) return FTOL; // set new search direction h to f = -Grad(x) - for (i = 0; i < n3; i++) h[i] = f[i]; + for (i = 0; i < nvec; i++) h[i] = fvec[i]; if (nextra_atom) for (m = 0; m < nextra_atom; m++) { fatom = fextra_atom[m]; hatom = hextra_atom[m]; n = extra_nlen[m]; for (i = 0; i < n; i++) hatom[i] = fatom[i]; } if (nextra_global) for (i = 0; i < nextra_global; i++) hextra[i] = fextra[i]; // output for thermo, dump, restart files if (output->next == ntimestep) { timer->stamp(); output->write(ntimestep); timer->stamp(TIME_OUTPUT); } } return MAXITER; } diff --git a/src/run.cpp b/src/run.cpp index 41ff1dd81..63de1ec2a 100644 --- a/src/run.cpp +++ b/src/run.cpp @@ -1,230 +1,231 @@ /* ---------------------------------------------------------------------- LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator http://lammps.sandia.gov, Sandia National Laboratories Steve Plimpton, sjplimp@sandia.gov Copyright (2003) Sandia Corporation. Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains certain rights in this software. This software is distributed under the GNU General Public License. See the README file in the top-level LAMMPS directory. ------------------------------------------------------------------------- */ #include "stdlib.h" #include "string.h" #include "run.h" #include "domain.h" #include "update.h" #include "integrate.h" #include "modify.h" #include "output.h" #include "finish.h" #include "input.h" #include "timer.h" #include "error.h" using namespace LAMMPS_NS; #define MIN(a,b) ((a) < (b) ? (a) : (b)) #define MAX(a,b) ((a) > (b) ? (a) : (b)) #define MAXLINE 2048 /* ---------------------------------------------------------------------- */ Run::Run(LAMMPS *lmp) : Pointers(lmp) {} /* ---------------------------------------------------------------------- */ void Run::command(int narg, char **arg) { if (narg < 1) error->all("Illegal run command"); if (domain->box_exist == 0) error->all("Run command before simulation box is defined"); int nsteps = atoi(arg[0]); // parse optional args int uptoflag = 0; int startflag = 0; int stopflag = 0; int start,stop; int preflag = 1; int postflag = 1; int nevery = 0; - int first,last,ncommands; + int ncommands = 0; + int first,last; int iarg = 1; while (iarg < narg) { if (strcmp(arg[iarg],"upto") == 0) { if (iarg+1 > narg) error->all("Illegal run command"); uptoflag = 1; iarg += 1; } else if (strcmp(arg[iarg],"start") == 0) { if (iarg+2 > narg) error->all("Illegal run command"); startflag = 1; start = atoi(arg[iarg+1]); iarg += 2; } else if (strcmp(arg[iarg],"stop") == 0) { if (iarg+2 > narg) error->all("Illegal run command"); stopflag = 1; stop = atoi(arg[iarg+1]); iarg += 2; } else if (strcmp(arg[iarg],"pre") == 0) { if (iarg+2 > narg) error->all("Illegal run command"); if (strcmp(arg[iarg+1],"no") == 0) preflag = 0; else if (strcmp(arg[iarg+1],"yes") == 0) preflag = 1; else error->all("Illegal run command"); iarg += 2; } else if (strcmp(arg[iarg],"post") == 0) { if (iarg+2 > narg) error->all("Illegal run command"); if (strcmp(arg[iarg+1],"no") == 0) postflag = 0; else if (strcmp(arg[iarg+1],"yes") == 0) postflag = 1; else error->all("Illegal run command"); iarg += 2; // all remaining args are commands // first,last = arg index of first/last commands // set ncommands = 0 if single command and it is NULL } else if (strcmp(arg[iarg],"every") == 0) { if (iarg+3 > narg) error->all("Illegal run command"); nevery = atoi(arg[iarg+1]); if (nevery <= 0) error->all("Illegal run command"); first = iarg+2; last = narg-1; ncommands = last-first + 1; if (ncommands == 1 && strcmp(arg[first],"NULL") == 0) ncommands = 0; iarg = narg; } else error->all("Illegal run command"); } // adjust nsteps if upto was specified if (uptoflag) nsteps -= update->ntimestep; // if nevery, make copies of arg strings that are commands // required because re-parsing commands via input->one() will wipe out args char **commands = NULL; if (nevery && ncommands > 0) { commands = new char*[ncommands]; ncommands = 0; for (int i = first; i <= last; i++) { int n = strlen(arg[i]) + 1; commands[ncommands] = new char[n]; strcpy(commands[ncommands],arg[i]); ncommands++; } } // error check if (uptoflag && nsteps < 0) error->all("Run command upto value is before current timestep"); if (startflag && start > update->ntimestep) error->all("Run command start value is after start of run"); if (stopflag && stop < update->ntimestep + nsteps) error->all("Run command stop value is before end of run"); // perform a single run // use start/stop to set begin/end step // if pre or 1st run, do System init/setup, // else just init timer and setup output // if post, do full Finish, else just print time update->whichflag = 1; if (nevery == 0) { update->nsteps = nsteps; update->firststep = update->ntimestep; update->laststep = update->ntimestep + nsteps; if (startflag) update->beginstep = start; else update->beginstep = update->firststep; if (stopflag) update->endstep = stop; else update->endstep = update->laststep; if (preflag || update->first_update == 0) { lmp->init(); update->integrate->setup(); } else { timer->init(); output->setup(0); } timer->barrier_start(TIME_LOOP); update->integrate->run(nsteps); timer->barrier_stop(TIME_LOOP); update->integrate->cleanup(); Finish finish(lmp); finish.end(postflag); // perform multiple runs optionally interleaved with invocation command(s) // use start/stop to set begin/end step // if pre or 1st iteration of multiple runs, do System init/setup, // else just init timer and setup output // if post or last iteration, do full Finish, else just print time } else { int iter = 0; int nleft = nsteps; while (nleft > 0 || iter == 0) { nsteps = MIN(nleft,nevery); update->nsteps = nsteps; update->firststep = update->ntimestep; update->laststep = update->ntimestep + nsteps; if (startflag) update->beginstep = start; else update->beginstep = update->firststep; if (stopflag) update->endstep = stop; else update->endstep = update->laststep; if (preflag || iter == 0) { lmp->init(); update->integrate->setup(); } else { timer->init(); output->setup(0); } timer->barrier_start(TIME_LOOP); update->integrate->run(nsteps); timer->barrier_stop(TIME_LOOP); update->integrate->cleanup(); Finish finish(lmp); if (postflag || nleft <= nsteps) finish.end(1); else finish.end(0); // wrap command invocation with clearstep/addstep // since a command may invoke computes via variables if (ncommands) { modify->clearstep_compute(); for (int i = 0; i < ncommands; i++) char *command = input->one(commands[i]); modify->addstep_compute(update->ntimestep + nevery); } nleft -= nsteps; iter++; } } update->whichflag = 0; update->firststep = update->laststep = 0; update->beginstep = update->endstep = 0; if (commands) { for (int i = 0; i < ncommands; i++) delete [] commands[i]; delete [] commands; } }