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

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

Browse Source
This commit is contained in:
sjplimp
2009-08-14 20:54:10 +00:00
parent 5234b15cd5
commit 33ff40dd7e
19 changed files with 1116 additions and 628 deletions
Download Patch File Download Diff File Expand all files Collapse all files

2
src/compute_pressure.cpp
View File

@ -236,10 +236,12 @@ void ComputePressure::compute_vector()
vector[0] = (ke_tensor[0] + virial[0]) * inv_volume * nktv2p;
vector[1] = (ke_tensor[1] + virial[1]) * inv_volume * nktv2p;
vector[3] = (ke_tensor[3] + virial[3]) * inv_volume * nktv2p;
vector[2] = vector[4] = vector[5] = 0.0;
} else {
vector[0] = virial[0] * inv_volume * nktv2p;
vector[1] = virial[1] * inv_volume * nktv2p;
vector[3] = virial[3] * inv_volume * nktv2p;
vector[2] = vector[4] = vector[5] = 0.0;
}
}
}

1
src/fix.h
View File

@ -105,6 +105,7 @@ class Fix : protected Pointers {
virtual void final_integrate_respa(int) {}
virtual void min_post_force(int) {}
virtual double min_energy(double *) {return 0.0;}
virtual void min_store() {}
virtual void min_step(double, double *) {}

24
src/fix_box_relax.cpp
View File

@ -318,34 +318,34 @@ void FixBoxRelax::min_store()
}
/* ----------------------------------------------------------------------
change the box dimensions by fraction ds = alpha*fextra
change the box dimensions by fraction ds = alpha*hextra
------------------------------------------------------------------------- */
void FixBoxRelax::min_step(double alpha, double *fextra)
void FixBoxRelax::min_step(double alpha, double *hextra)
{
if (press_couple == XYZ) {
ds[0] = ds[1] = ds[2] = alpha*fextra[0];
ds[0] = ds[1] = ds[2] = alpha*hextra[0];
} else {
if (p_flag[0]) ds[0] = alpha*fextra[0];
if (p_flag[1]) ds[1] = alpha*fextra[1];
if (p_flag[2]) ds[2] = alpha*fextra[2];
if (p_flag[0]) ds[0] = alpha*hextra[0];
if (p_flag[1]) ds[1] = alpha*hextra[1];
if (p_flag[2]) ds[2] = alpha*hextra[2];
}
remap();
}
/* ----------------------------------------------------------------------
max allowed step size along fextra
max allowed step size along hextra
------------------------------------------------------------------------- */
double FixBoxRelax::max_alpha(double *fextra)
double FixBoxRelax::max_alpha(double *hextra)
{
double alpha = 0.0;
if (press_couple == XYZ) alpha = vmax/fabs(fextra[0]);
if (press_couple == XYZ) alpha = vmax/fabs(hextra[0]);
else {
alpha = vmax/fabs(fextra[0]);
alpha = MIN(alpha,vmax/fabs(fextra[1]));
alpha = MIN(alpha,vmax/fabs(fextra[2]));
alpha = vmax/fabs(hextra[0]);
alpha = MIN(alpha,vmax/fabs(hextra[1]));
alpha = MIN(alpha,vmax/fabs(hextra[2]));
}
return alpha;
}

2
src/fix_box_relax.h
View File

@ -24,11 +24,13 @@ class FixBoxRelax : public Fix {
~FixBoxRelax();
int setmask();
void init();
double min_energy(double *);
void min_store();
void min_step(double, double *);
double max_alpha(double *);
int min_dof();
int modify_param(int, char **);
private:

195
src/fix_minimize.cpp
View File

@ -11,8 +11,10 @@
See the README file in the top-level LAMMPS directory.
------------------------------------------------------------------------- */
#include "stdlib.h"
#include "fix_minimize.h"
#include "atom.h"
#include "domain.h"
#include "memory.h"
#include "error.h"
@ -23,13 +25,13 @@ using namespace LAMMPS_NS;
FixMinimize::FixMinimize(LAMMPS *lmp, int narg, char **arg) :
Fix(lmp, narg, arg)
{
// perform initial allocation of atom-based arrays
// register with Atom class
nvector = 0;
peratom = NULL;
vectors = NULL;
// register callback to this fix from Atom class
// don't perform initial allocation here, must wait until add_vector()
gradient = NULL;
searchdir = NULL;
x0 = NULL;
grow_arrays(atom->nmax);
atom->add_callback(0);
}
@ -41,11 +43,11 @@ FixMinimize::~FixMinimize()
atom->delete_callback(id,0);
// delete locally stored arrays
// delete locally stored data
memory->destroy_2d_double_array(gradient);
memory->destroy_2d_double_array(searchdir);
memory->destroy_2d_double_array(x0);
memory->sfree(peratom);
for (int m = 0; m < nvector; m++) memory->sfree(vectors[m]);
memory->sfree(vectors);
}
/* ---------------------------------------------------------------------- */
@ -55,13 +57,120 @@ int FixMinimize::setmask()
return 0;
}
/* ----------------------------------------------------------------------
allocate/initialize memory for a new vector with N elements per atom
------------------------------------------------------------------------- */
void FixMinimize::add_vector(int n)
{
peratom = (int *)
memory->srealloc(peratom,(nvector+1)*sizeof(int),"minimize:peratom");
peratom[nvector] = n;
vectors = (double **)
memory->srealloc(vectors,(nvector+1)*sizeof(double *),"minimize:vectors");
vectors[nvector] = (double *)
memory->smalloc(atom->nmax*n*sizeof(double),"minimize:vector");
int ntotal = n*atom->nlocal;
for (int i = 0; i < ntotal; i++) vectors[nvector][i] = 0.0;
nvector++;
}
/* ----------------------------------------------------------------------
return a pointer to the Mth vector
------------------------------------------------------------------------- */
double *FixMinimize::request_vector(int m)
{
return vectors[m];
}
/* ----------------------------------------------------------------------
store box size at beginning of line search
------------------------------------------------------------------------- */
void FixMinimize::store_box()
{
boxlo[0] = domain->boxlo[0];
boxlo[1] = domain->boxlo[1];
boxlo[2] = domain->boxlo[2];
boxhi[0] = domain->boxhi[0];
boxhi[1] = domain->boxhi[1];
boxhi[2] = domain->boxhi[2];
}
/* ----------------------------------------------------------------------
reset x0 for atoms that moved across PBC via reneighboring in line search
x0 = 1st vector
must do minimum_image using original box stored at beginning of line search
swap & set_global_box() change to original box, then restore current box
------------------------------------------------------------------------- */
void FixMinimize::reset_coords()
{
box_swap();
domain->set_global_box();
double **x = atom->x;
double *x0 = vectors[0];
int nlocal = atom->nlocal;
double dx,dy,dz,dx0,dy0,dz0;
int n = 0;
for (int i = 0; i < nlocal; i++) {
dx = dx0 = x[i][0] - x0[n];
dy = dy0 = x[i][1] - x0[n+1];
dz = dz0 = x[i][2] - x0[n+2];
domain->minimum_image(dx,dy,dz);
if (dx != dx0) x0[n] = x[i][0] - dx;
if (dy != dy0) x0[n+1] = x[i][1] - dy;
if (dz != dz0) x0[n+2] = x[i][2] - dz;
n += 3;
}
box_swap();
domain->set_global_box();
}
/* ----------------------------------------------------------------------
swap current box size with stored box size
------------------------------------------------------------------------- */
void FixMinimize::box_swap()
{
double tmp;
tmp = boxlo[0];
boxlo[0] = domain->boxlo[0];
domain->boxlo[0] = tmp;
tmp = boxlo[1];
boxlo[1] = domain->boxlo[1];
domain->boxlo[1] = tmp;
tmp = boxlo[2];
boxlo[2] = domain->boxlo[2];
domain->boxlo[2] = tmp;
tmp = boxhi[0];
boxhi[0] = domain->boxhi[0];
domain->boxhi[0] = tmp;
tmp = boxhi[1];
boxhi[1] = domain->boxhi[1];
domain->boxhi[1] = tmp;
tmp = boxhi[2];
boxhi[2] = domain->boxhi[2];
domain->boxhi[2] = tmp;
}
/* ----------------------------------------------------------------------
memory usage of local atom-based arrays
------------------------------------------------------------------------- */
double FixMinimize::memory_usage()
{
double bytes = 3 * atom->nmax*3 * sizeof(double);
double bytes = 0.0;
for (int m = 0; m < nvector; m++)
bytes += atom->nmax*peratom[m]*sizeof(double);
return bytes;
}
@ -71,12 +180,9 @@ double FixMinimize::memory_usage()
void FixMinimize::grow_arrays(int nmax)
{
gradient =
memory->grow_2d_double_array(gradient,nmax,3,"fix_minimize:gradient");
searchdir =
memory->grow_2d_double_array(searchdir,nmax,3,"fix_minimize:searchdir");
x0 =
memory->grow_2d_double_array(x0,nmax,3,"fix_minimize:x0");
for (int m = 0; m < nvector; m++)
vectors[m] = (double *) memory->srealloc(vectors[m],peratom[m]*nmax,
"minimize:vector");
}
/* ----------------------------------------------------------------------
@ -85,15 +191,14 @@ void FixMinimize::grow_arrays(int nmax)
void FixMinimize::copy_arrays(int i, int j)
{
gradient[j][0] = gradient[i][0];
gradient[j][1] = gradient[i][1];
gradient[j][2] = gradient[i][2];
searchdir[j][0] = searchdir[i][0];
searchdir[j][1] = searchdir[i][1];
searchdir[j][2] = searchdir[i][2];
x0[j][0] = x0[i][0];
x0[j][1] = x0[i][1];
x0[j][2] = x0[i][2];
int m,iper,nper,ni,nj;
for (m = 0; m < nvector; m++) {
nper = peratom[m];
ni = nper*i;
nj = nper*j;
for (iper = 0; iper < nper; iper++) vectors[m][nj++] = vectors[m][ni++];
}
}
/* ----------------------------------------------------------------------
@ -102,16 +207,15 @@ void FixMinimize::copy_arrays(int i, int j)
int FixMinimize::pack_exchange(int i, double *buf)
{
buf[0] = gradient[i][0];
buf[1] = gradient[i][1];
buf[2] = gradient[i][2];
buf[3] = searchdir[i][0];
buf[4] = searchdir[i][1];
buf[5] = searchdir[i][2];
buf[6] = x0[i][0];
buf[7] = x0[i][1];
buf[8] = x0[i][2];
return 9;
int m,iper,nper,ni;
int n = 0;
for (m = 0; m < nvector; m++) {
nper = peratom[m];
ni = nper*i;
for (iper = 0; iper < nper; iper++) buf[n++] = vectors[m][ni++];
}
return n;
}
/* ----------------------------------------------------------------------
@ -120,14 +224,13 @@ int FixMinimize::pack_exchange(int i, double *buf)
int FixMinimize::unpack_exchange(int nlocal, double *buf)
{
gradient[nlocal][0] = buf[0];
gradient[nlocal][1] = buf[1];
gradient[nlocal][2] = buf[2];
searchdir[nlocal][0] = buf[3];
searchdir[nlocal][1] = buf[4];
searchdir[nlocal][2] = buf[5];
x0[nlocal][0] = buf[6];
x0[nlocal][1] = buf[7];
x0[nlocal][2] = buf[8];
return 9;
int m,iper,nper,ni;
int n = 0;
for (m = 0; m < nvector; m++) {
int nper = peratom[m];
ni = nper*nlocal;
for (iper = 0; iper < nper; iper++) vectors[m][ni++] = buf[n++];
}
return n;
}

18
src/fix_minimize.h
View File

@ -19,10 +19,9 @@
namespace LAMMPS_NS {
class FixMinimize : public Fix {
public:
double **gradient,**searchdir; // gradient vectors
double **x0; // initial coords at start of linesearch
friend class MinLineSearch;
public:
FixMinimize(class LAMMPS *, int, char **);
~FixMinimize();
int setmask();
@ -33,6 +32,19 @@ class FixMinimize : public Fix {
void copy_arrays(int, int);
int pack_exchange(int, double *);
int unpack_exchange(int, double *);
void add_vector(int);
double *request_vector(int);
void store_box();
void reset_coords();
private:
int nvector;
int *peratom;
double **vectors;
double boxlo[3],boxhi[3];
void box_swap();
};
}

2
src/lammps.cpp
View File

@ -304,7 +304,7 @@ void LAMMPS::create()
void LAMMPS::init()
{
update->init();
force->init();
force->init(); // pair must come after update due to minimizer
domain->init();
atom->init(); // atom must come after force:
// atom deletes extra array

615
src/min.cpp
View File

@ -41,28 +41,11 @@
#include "output.h"
#include "thermo.h"
#include "timer.h"
#include "memory.h"
#include "error.h"
using namespace LAMMPS_NS;
// ALPHA_MAX = max alpha allowed to avoid long backtracks
// ALPHA_REDUCE = reduction ratio, should be in range [0.5,1)
// BACKTRACK_SLOPE, should be in range (0,0.5]
// QUADRATIC_TOL = tolerance on alpha0, should be in range [0.1,1)
// IDEAL_TOL = ideal energy tolerance for backtracking
// EPS_QUAD = tolerance for quadratic projection
#define ALPHA_MAX 1.0
#define ALPHA_REDUCE 0.5
#define BACKTRACK_SLOPE 0.4
#define QUADRATIC_TOL 0.1
#define IDEAL_TOL 1.0e-8
#define EPS_QUAD 1.0e-28
// same as in other min classes
enum{MAXITER,MAXEVAL,ETOL,FTOL,DOWNHILL,ZEROALPHA,ZEROFORCE,ZEROQUAD};
#define MIN(A,B) ((A) < (B)) ? (A) : (B)
#define MAX(A,B) ((A) > (B)) ? (A) : (B)
@ -76,7 +59,14 @@ Min::Min(LAMMPS *lmp) : Pointers(lmp)
elist_atom = NULL;
vlist_global = vlist_atom = NULL;
fextra = gextra = hextra = NULL;
nextra_global = 0;
fextra = NULL;
nextra_atom = 0;
xextra_atom = fextra_atom = NULL;
extra_peratom = extra_nlen = NULL;
extra_max = NULL;
requestor = NULL;
}
/* ---------------------------------------------------------------------- */
@ -88,15 +78,20 @@ Min::~Min()
delete [] vlist_atom;
delete [] fextra;
delete [] gextra;
delete [] hextra;
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 gradient vectors
// create fix needed for storing atom-based quantities
// will delete it at end of run
char **fixarg = new char*[3];
@ -107,10 +102,25 @@ void Min::init()
delete [] fixarg;
fix_minimize = (FixMinimize *) modify->fix[modify->nfix-1];
// zero gradient vectors before first atom exchange
// 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()
setup_vectors();
for (int i = 0; i < ndof; i++) h[i] = g[i] = 0.0;
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
@ -139,6 +149,8 @@ void Min::init()
neigh_delay = neighbor->delay;
neigh_dist_check = neighbor->dist_check;
// reset reneighboring criteria if necessary
if (neigh_every != 1 || neigh_delay != 0 || neigh_dist_check != 1) {
if (comm->me == 0)
error->warning("Resetting reneighboring criteria during minimization");
@ -148,10 +160,9 @@ void Min::init()
neighbor->delay = 0;
neighbor->dist_check = 1;
// set ptr to linemin function
// style-specific initialization
if (linestyle == 0) linemin = &Min::linemin_backtrack;
else if (linestyle == 1) linemin = &Min::linemin_quadratic;
init_style();
}
/* ----------------------------------------------------------------------
@ -162,6 +173,28 @@ 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];
// 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
@ -177,7 +210,6 @@ void Min::setup()
if (triclinic) domain->lamda2x(atom->nlocal+atom->nghost);
neighbor->build();
neighbor->ncalls = 0;
setup_vectors();
// compute all forces
@ -202,6 +234,10 @@ void Min::setup()
modify->setup(vflag);
output->setup(1);
// atoms may have migrated in comm->exchange()
reset_vectors();
}
/* ----------------------------------------------------------------------
@ -210,9 +246,6 @@ void Min::setup()
void Min::run()
{
int i;
double tmp,*f;
// possible stop conditions
char *stopstrings[] = {"max iterations","max force evaluations",
@ -221,57 +254,21 @@ void Min::run()
"linesearch alpha is zero",
"forces are zero","quadratic factors are zero"};
// set initial force & energy
setup();
// setup any extra dof due to fixes
// can't be done until now b/c update init() comes before modify init()
delete [] fextra;
delete [] gextra;
delete [] hextra;
fextra = NULL;
gextra = NULL;
hextra = NULL;
nextra = modify->min_dof();
if (nextra) {
fextra = new double[nextra];
gextra = new double[nextra];
hextra = new double[nextra];
}
// compute potential energy of system
// normalize energy if thermo PE does
// 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];
ecurrent = pe_compute->compute_scalar();
if (nextra) ecurrent += modify->min_energy(fextra);
if (output->thermo->normflag) ecurrent /= atom->natoms;
// 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;
f = NULL;
if (ndof) f = atom->f[0];
tmp = 0.0;
for (i = 0; i < ndof; i++) tmp += f[i]*f[i];
MPI_Allreduce(&tmp,&fnorm2_init,1,MPI_DOUBLE,MPI_SUM,world);
if (nextra)
for (i = 0; i < nextra; i++) fnorm2_init += fextra[i]*fextra[i];
fnorm2_init = sqrt(fnorm2_init);
tmp = 0.0;
for (i = 0; i < ndof; i++) tmp = MAX(fabs(f[i]),tmp);
MPI_Allreduce(&tmp,&fnorminf_init,1,MPI_DOUBLE,MPI_MAX,world);
if (nextra)
for (i = 0; i < nextra; i++)
fnorminf_init = MAX(fabs(fextra[i]),fnorminf_init);
fnorm2_init = sqrt(fnorm_sqr());
fnorminf_init = fnorm_inf();
// minimizer iterations
@ -298,16 +295,13 @@ void Min::run()
output->next_thermo = update->ntimestep;
modify->addstep_compute_all(update->ntimestep);
int ntmp;
double *xtmp,*htmp,*x0tmp,etmp;
eng_force(&ntmp,&xtmp,&htmp,&x0tmp,&etmp,0);
ecurrent = energy_force(0);
output->write(update->ntimestep);
}
timer->barrier_stop(TIME_LOOP);
// delete fix at end of run, so its atom arrays won't persist
// delete fix_minimize at end of run
modify->delete_fix("MINIMIZE");
@ -320,39 +314,22 @@ void Min::run()
// stats for Finish to print
efinal = ecurrent;
f = NULL;
if (ndof) f = atom->f[0];
tmp = 0.0;
for (i = 0; i < ndof; i++) tmp += f[i]*f[i];
MPI_Allreduce(&tmp,&fnorm2_final,1,MPI_DOUBLE,MPI_SUM,world);
if (nextra)
for (i = 0; i < nextra; i++)
fnorm2_final += fextra[i]*fextra[i];
fnorm2_final = sqrt(fnorm2_final);
tmp = 0.0;
for (i = 0; i < ndof; i++) tmp = MAX(fabs(f[i]),tmp);
MPI_Allreduce(&tmp,&fnorminf_final,1,MPI_DOUBLE,MPI_MAX,world);
if (nextra)
for (i = 0; i < nextra; i++)
fnorminf_final = MAX(fabs(fextra[i]),fnorminf_final);
fnorm2_final = sqrt(fnorm_sqr());
fnorminf_final = fnorm_inf();
}
/* ----------------------------------------------------------------------
evaluate potential energy and forces
may migrate atoms
if resetflag = 1, update x0 by PBC for atoms that migrate
new energy stored in ecurrent and returned (in case caller not in class)
negative gradient will be stored in atom->f
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
------------------------------------------------------------------------- */
void Min::eng_force(int *pndof, double **px, double **ph, double **px0,
double *peng, int resetflag)
double Min::energy_force(int resetflag)
{
// check for reneighboring
// always communicate since minimizer moved atoms
// if reneighbor, have to setup_vectors() since atoms migrated
int nflag = neighbor->decide();
@ -375,34 +352,6 @@ void Min::eng_force(int *pndof, double **px, double **ph, double **px0,
timer->stamp(TIME_COMM);
neighbor->build();
timer->stamp(TIME_NEIGHBOR);
setup_vectors();
// update x0 for atoms that migrated
// must do minimum_image on box size when x0 was stored
// domain->set_global_box() changes to x0 box, then restores current box
if (resetflag) {
box_swap();
domain->set_global_box();
double **x = atom->x;
double **x0 = fix_minimize->x0;
int nlocal = atom->nlocal;
double dx,dy,dz,dx0,dy0,dz0;
for (int i = 0; i < nlocal; i++) {
dx = dx0 = x[i][0] - x0[i][0];
dy = dy0 = x[i][1] - x0[i][1];
dz = dz0 = x[i][2] - x0[i][2];
domain->minimum_image(dx,dy,dz);
if (dx != dx0) x0[i][0] = x[i][0] - dx;
if (dy != dy0) x0[i][1] = x[i][1] - dy;
if (dz != dz0) x0[i][2] = x[i][2] - dz;
}
box_swap();
domain->set_global_box();
}
}
ev_set(update->ntimestep);
@ -440,33 +389,20 @@ void Min::eng_force(int *pndof, double **px, double **ph, double **px0,
// compute potential energy of system
// normalize if thermo PE does
ecurrent = pe_compute->compute_scalar();
if (nextra) ecurrent += modify->min_energy(fextra);
if (output->thermo->normflag) ecurrent /= atom->natoms;
double energy = pe_compute->compute_scalar();
if (nextra_global) energy += modify->min_energy(fextra);
if (output->thermo->normflag) energy /= atom->natoms;
// return updated ptrs to caller since atoms may have migrated
// if reneighbored, atoms migrated
// if resetflag = 1, update x0 of atoms crossing PBC
// reset vectors used by lo-level minimizer
*pndof = ndof;
if (ndof) *px = atom->x[0];
else *px = NULL;
*ph = h;
*px0 = x0;
*peng = ecurrent;
}
if (nflag) {
if (resetflag) fix_minimize->reset_coords();
reset_vectors();
}
/* ----------------------------------------------------------------------
set ndof and vector pointers after atoms have migrated
------------------------------------------------------------------------- */
void Min::setup_vectors()
{
ndof = 3 * atom->nlocal;
if (ndof) g = fix_minimize->gradient[0];
else g = NULL;
if (ndof) h = fix_minimize->searchdir[0];
else h = NULL;
if (ndof) x0 = fix_minimize->x0[0];
else x0 = NULL;
return energy;
}
/* ----------------------------------------------------------------------
@ -501,282 +437,30 @@ void Min::force_clear()
}
/* ----------------------------------------------------------------------
line minimization methods
find minimum-energy starting at x along dir direction
input: n = # of degrees of freedom on this proc
x = ptr to atom->x[0] as vector
dir = search direction as vector
x0 = ptr to fix->x0[0] as vector, for storing initial coords
eoriginal = energy at initial x
maxdist = max distance to move any atom coord
output: return 0 if successful move, non-zero alpha
return non-zero if failed
alpha = distance moved along dir to set x to minimun eng config
caller has several quantities set via last call to eng_force()
must insure last call to eng_force() is consistent with returns
if fail, eng_force() of original x
if succeed, eng_force() at x + alpha*dir
atom->x = coords at new configuration
atom->f = force (-Grad) is evaulated at new configuration
ecurrent = energy of new configuration
NOTE: when call eng_force: n,x,dir,x0,eng may change due to atom migration
updated values are returned by eng_force()
b/c of migration, linemin routines CANNOT store atom-based quantities
clear force on own & ghost atoms
setup and clear other arrays as needed
------------------------------------------------------------------------- */
/* ----------------------------------------------------------------------
linemin: backtracking line search (Proc 3.1, p 41 in Nocedal and Wright)
uses no gradient info, but should be very robust
start at maxdist, backtrack until energy decrease is sufficient
------------------------------------------------------------------------- */
int Min::linemin_backtrack(int n, double *x, double *dir,
double *x0, double eoriginal, double maxdist,
double &alpha, int &nfunc)
void Min::request(Pair *pair, int peratom, double maxvalue)
{
int i,m;
double fdotdirall,fdotdirme,hmax,hme,alpha_extra;
double eng,de_ideal,de;
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");
double *f = NULL;
if (n) f = atom->f[0];
// fdotdirall = projection of search dir along downhill gradient
// if search direction is not downhill, exit with error
fdotdirme = 0.0;
for (i = 0; i < n; i++) fdotdirme += f[i]*dir[i];
MPI_Allreduce(&fdotdirme,&fdotdirall,1,MPI_DOUBLE,MPI_SUM,world);
if (nextra)
for (i = 0; i < nextra; i++) fdotdirall += fextra[i]*hextra[i];
if (output->thermo->normflag) fdotdirall /= atom->natoms;
if (fdotdirall <= 0.0) return DOWNHILL;
// initial alpha = stepsize to change any atom coord by maxdist
// alpha <= ALPHA_MAX, else backtrack from huge value when forces are tiny
// if all search dir components are already 0.0, exit with error
hme = 0.0;
for (i = 0; i < n; i++) hme = MAX(hme,fabs(dir[i]));
MPI_Allreduce(&hme,&hmax,1,MPI_DOUBLE,MPI_MAX,world);
alpha = MIN(ALPHA_MAX,maxdist/hmax);
if (nextra) {
double alpha_extra = modify->max_alpha(hextra);
alpha = MIN(alpha,alpha_extra);
for (i = 0; i < nextra; i++)
hmax = MAX(hmax,fabs(hextra[i]));
}
if (hmax == 0.0) return ZEROFORCE;
// store coords and other dof at start of linesearch
box_store();
for (i = 0; i < n; i++) x0[i] = x[i];
if (nextra) modify->min_store();
// backtrack with alpha until energy decrease is sufficient
while (1) {
if (nextra) modify->min_step(0.0,hextra);
for (i = 0; i < n; i++) x[i] = x0[i];
if (nextra) modify->min_step(alpha,hextra);
for (i = 0; i < n; i++) x[i] += alpha*dir[i];
eng_force(&n,&x,&dir,&x0,&eng,1);
nfunc++;
// if energy change is better than ideal, exit with success
de_ideal = -BACKTRACK_SLOPE*alpha*fdotdirall;
de = eng - eoriginal;
if (de <= de_ideal) return 0;
// reduce alpha
alpha *= ALPHA_REDUCE;
// backtracked all the way to 0.0
// reset to starting point, exit with error
if (alpha <= 0.0 || de_ideal >= -IDEAL_TOL) {
if (nextra) modify->min_step(0.0,hextra);
for (i = 0; i < n; i++) x[i] = x0[i];
eng_force(&n,&x,&dir,&x0,&eng,0);
nfunc++;
return ZEROALPHA;
}
}
}
/* ----------------------------------------------------------------------
linemin: quadratic line search (adapted from Dennis and Schnabel)
basic idea is to backtrack until change in energy is sufficiently small
based on ENERGY_QUADRATIC, then use a quadratic approximation
using forces at two alpha values to project to minimum
use forces rather than energy change to do projection
this is b/c the forces are going to zero and can become very small
unlike energy differences which are the difference of two finite
values and are thus limited by machine precision
two changes that were critical to making this method work:
a) limit maximum step to alpha <= 1
b) ignore energy criterion if delE <= ENERGY_QUADRATIC
several other ideas also seemed to help:
c) making each step from starting point (alpha = 0), not previous alpha
d) quadratic model based on forces, not energy
e) exiting immediately if f.dir <= 0 (search direction not downhill)
so that CG can restart
a,c,e were also adopted for the backtracking linemin function
------------------------------------------------------------------------- */
int Min::linemin_quadratic(int n, double *x, double *dir,
double *x0, double eoriginal, double maxdist,
double &alpha, int &nfunc)
{
int i,m;
double fdotdirall,fdotdirme,hmax,hme,alphamax,alpha_extra;
double eng,de_ideal,de;
double delfh,engprev,relerr,alphaprev,fhprev,ff,fh,alpha0,fh0,ff0;
double dot[2],dotall[2];
double *f = atom->f[0];
// fdotdirall = projection of search dir along downhill gradient
// if search direction is not downhill, exit with error
fdotdirme = 0.0;
for (i = 0; i < n; i++) fdotdirme += f[i]*dir[i];
MPI_Allreduce(&fdotdirme,&fdotdirall,1,MPI_DOUBLE,MPI_SUM,world);
if (nextra)
for (i = 0; i < nextra; i++) fdotdirall += fextra[i]*hextra[i];
if (output->thermo->normflag) fdotdirall /= atom->natoms;
if (fdotdirall <= 0.0) return DOWNHILL;
// initial alpha = stepsize to change any atom coord by maxdist
// alpha <= ALPHA_MAX, else backtrack from huge value when forces are tiny
// if all search dir components are already 0.0, exit with error
hme = 0.0;
for (i = 0; i < n; i++) hme = MAX(hme,fabs(dir[i]));
MPI_Allreduce(&hme,&hmax,1,MPI_DOUBLE,MPI_MAX,world);
alpha = MIN(ALPHA_MAX,maxdist/hmax);
if (nextra) {
double alpha_extra = modify->max_alpha(hextra);
alpha = MIN(alpha,alpha_extra);
for (i = 0; i < nextra; i++)
hmax = MAX(hmax,fabs(hextra[i]));
}
if (hmax == 0.0) return ZEROFORCE;
// store coords and other dof at start of linesearch
box_store();
for (i = 0; i < n; i++) x0[i] = x[i];
if (nextra) modify->min_store();
// backtrack with alpha until energy decrease is sufficient
// or until get to small energy change, then perform quadratic projection
fhprev = fdotdirall;
engprev = eoriginal;
alphaprev = 0.0;
while (1) {
if (nextra) modify->min_step(0.0,hextra);
for (i = 0; i < n; i++) x[i] = x0[i];
if (nextra) modify->min_step(alpha,hextra);
for (i = 0; i < n; i++) x[i] += alpha*dir[i];
eng_force(&n,&x,&dir,&x0,&eng,1);
nfunc++;
// compute new fh, alpha, delfh
dot[0] = dot[1] = 0.0;
for (i = 0; i < ndof; i++) {
dot[0] += f[i]*f[i];
dot[1] += f[i]*dir[i];
}
MPI_Allreduce(dot,dotall,2,MPI_DOUBLE,MPI_SUM,world);
if (nextra) {
for (i = 0; i < nextra; i++) {
dotall[0] += fextra[i]*fextra[i];
dotall[1] += fextra[i]*hextra[i];
}
}
ff = dotall[0];
fh = dotall[1];
if (output->thermo->normflag) {
ff /= atom->natoms;
fh /= atom->natoms;
}
delfh = fh - fhprev;
// if fh or delfh is epsilon, reset to starting point, exit with error
if (fabs(fh) < EPS_QUAD || fabs(delfh) < EPS_QUAD) {
if (nextra) modify->min_step(0.0,hextra);
for (i = 0; i < n; i++) x[i] = x0[i];
eng_force(&n,&x,&dir,&x0,&eng,0);
nfunc++;
return ZEROQUAD;
}
// check if ready for quadratic projection, equivalent to secant method
// alpha0 = projected alpha
relerr = fabs(1.0+(0.5*alpha*(alpha-alphaprev)*(fh+fhprev)-eng)/engprev);
alpha0 = alpha - (alpha-alphaprev)*fh/delfh;
if (relerr <= QUADRATIC_TOL && alpha0 > 0.0) {
if (nextra) modify->min_step(0.0,hextra);
for (i = 0; i < n; i++) x[i] = x0[i];
if (nextra) modify->min_step(alpha0,hextra);
for (i = 0; i < n; i++) x[i] += alpha0*dir[i];
eng_force(&n,&x,&dir,&x0,&eng,1);
nfunc++;
// if backtracking energy change is better than ideal, exit with success
de_ideal = -BACKTRACK_SLOPE*alpha0*fdotdirall;
de = eng - eoriginal;
if (de <= de_ideal || de_ideal >= -IDEAL_TOL) return 0;
// drop back from alpha0 to alpha
if (nextra) modify->min_step(0.0,hextra);
for (i = 0; i < n; i++) x[i] = x0[i];
if (nextra) modify->min_step(alpha,hextra);
for (i = 0; i < n; i++) x[i] += alpha*dir[i];
eng_force(&n,&x,&dir,&x0,&eng,1);
nfunc++;
}
// if backtracking energy change is better than ideal, exit with success
de_ideal = -BACKTRACK_SLOPE*alpha*fdotdirall;
de = eng - eoriginal;
if (de <= de_ideal) return 0;
// save previous state
fhprev = fh;
engprev = eng;
alphaprev = alpha;
// reduce alpha
alpha *= ALPHA_REDUCE;
// backtracked all the way to 0.0
// reset to starting point, exit with error
if (alpha <= 0.0 || de_ideal >= -IDEAL_TOL) {
if (nextra) modify->min_step(0.0,hextra);
for (i = 0; i < n; i++) x[i] = x0[i];
eng_force(&n,&x,&dir,&x0,&eng,0);
nfunc++;
return ZEROALPHA;
}
}
requestor[nextra_atom] = pair;
extra_peratom[nextra_atom] = peratom;
extra_max[nextra_atom] = maxvalue;
nextra_atom++;
}
/* ---------------------------------------------------------------------- */
@ -882,45 +566,62 @@ void Min::ev_set(int ntimestep)
}
/* ----------------------------------------------------------------------
store box size at beginning of line search
compute and return ||force||_2^2
------------------------------------------------------------------------- */
void Min::box_store()
double Min::fnorm_sqr()
{
boxlo0[0] = domain->boxlo[0];
boxlo0[1] = domain->boxlo[1];
boxlo0[2] = domain->boxlo[2];
int i,n;
double *fatom;
boxhi0[0] = domain->boxhi[0];
boxhi0[1] = domain->boxhi[1];
boxhi0[2] = domain->boxhi[2];
double local_norm2_sqr = 0.0;
for (i = 0; i < n3; i++) local_norm2_sqr += f[i]*f[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;
}
/* ----------------------------------------------------------------------
swap current box size with stored box size
compute and return ||force||_inf
------------------------------------------------------------------------- */
void Min::box_swap()
double Min::fnorm_inf()
{
double tmp;
int i,n;
double *fatom;
tmp = boxlo0[0];
boxlo0[0] = domain->boxlo[0];
domain->boxlo[0] = tmp;
tmp = boxlo0[1];
boxlo0[1] = domain->boxlo[1];
domain->boxlo[1] = tmp;
tmp = boxlo0[2];
boxlo0[2] = domain->boxlo[2];
domain->boxlo[2] = tmp;
double local_norm_inf = 0.0;
for (i = 0; i < n3; i++)
local_norm_inf = MAX(fabs(f[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);
}
}
tmp = boxhi0[0];
boxhi0[0] = domain->boxhi[0];
domain->boxhi[0] = tmp;
tmp = boxhi0[1];
boxhi0[1] = domain->boxhi[1];
domain->boxhi[1] = tmp;
tmp = boxhi0[2];
boxhi0[2] = domain->boxhi[2];
domain->boxhi[2] = tmp;
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;
}

62
src/min.h
View File

@ -29,17 +29,23 @@ class Min : protected Pointers {
Min(class LAMMPS *);
virtual ~Min();
void init();
void setup();
void run();
virtual void init_style() {}
virtual void setup_style() = 0;
virtual void reset_vectors() = 0;
virtual int iterate(int) = 0;
void request(class Pair *, int, double);
double memory_usage() {return 0.0;}
void modify_params(int, char **);
virtual int iterate(int) = 0;
protected:
int eflag,vflag; // flags for energy/virial computation
int virial_style; // compute virial explicitly or implicitly
double dmax; // max dist to move any atom in one linesearch
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
@ -49,45 +55,45 @@ class Min : protected Pointers {
class Compute **vlist_atom;
int pairflag,torqueflag;
int neigh_every,neigh_delay,neigh_dist_check; // copies of reneigh criteria
int triclinic; // 0 if domain is orthog, 1 if triclinic
class FixMinimize *fix_minimize; // fix that stores gradient vecs
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 mindist,maxdist; // min/max dist for coord delta in line search
int ndof; // # of degrees-of-freedom on this proc
double *g,*h; // local portion of gradient, searchdir vectors
double *x0; // coords at start of linesearch
double ndoftotal; // total dof for entire problem
int nextra; // extra dof due to fixes
double *fextra; // vectors for extra dof
double *gextra;
double *hextra;
int n3; // local atomic dof
double *x; // variables for atomic dof, as 1d vector
double *f; // force vector for atomic dof, as 1d vector
double boxlo0[3]; // box size at start of linesearch
double boxhi0[3];
int nextra_global; // # of extra global dof due to fixes
double *fextra; // force vector for extra global dof
// xextra is stored by fix
// ptr to linemin functions
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
void setup();
void eng_force(int *, double **, double **, double **, double *, int);
void setup_vectors();
int neigh_every,neigh_delay,neigh_dist_check; // neighboring params
double energy_force(int);
void force_clear();
typedef int (Min::*FnPtr)(int, double *, double *, double *, double,
double, double &, int &);
FnPtr linemin;
int linemin_backtrack(int, double *, double *, double *, double,
double, double &, int &);
int linemin_quadratic(int, double *, double *, double *, double,
double, double &, int &);
double compute_force_norm_sqr();
double compute_force_norm_inf();
void ev_setup();
void ev_set(int);
void box_store();
void box_swap();
double fnorm_sqr();
double fnorm_inf();
};
}

114
src/min_cg.cpp
View File

@ -22,6 +22,8 @@
using namespace LAMMPS_NS;
#define MAXATOMS 0x7FFFFFFF
// EPS_ENERGY = minimum normalization for energy tolerance
#define EPS_ENERGY 1.0e-8
@ -35,47 +37,48 @@ enum{MAXITER,MAXEVAL,ETOL,FTOL,DOWNHILL,ZEROALPHA,ZEROFORCE,ZEROQUAD};
/* ---------------------------------------------------------------------- */
MinCG::MinCG(LAMMPS *lmp) : Min(lmp) {}
MinCG::MinCG(LAMMPS *lmp) : MinLineSearch(lmp) {}
/* ----------------------------------------------------------------------
minimization via conjugate gradient iterations
------------------------------------------------------------------------- */
int MinCG::iterate(int n)
int MinCG::iterate(int niter_max)
{
int i,fail,ntimestep;
int i,m,n,fail,ntimestep;
double beta,gg,dot[2],dotall[2];
double *fatom,*gatom,*hatom;
double *x = NULL;
double *f = NULL;
// nlimit = max # of CG iterations before restarting
// set to ndoftotal unless too big
int ndoftotal;
MPI_Allreduce(&ndof,&ndoftotal,1,MPI_INT,MPI_SUM,world);
ndoftotal += nextra;
int nlimit = static_cast<int> (MIN(MAXATOMS,ndoftotal));
if (ndof) f = atom->f[0];
for (i = 0; i < ndof; i++) h[i] = g[i] = f[i];
if (nextra)
for (i = 0; i < nextra; i++)
hextra[i] = gextra[i] = fextra[i];
// initialize working vectors
dot[0] = 0.0;
for (i = 0; i < ndof; i++) dot[0] += f[i]*f[i];
MPI_Allreduce(dot,&gg,1,MPI_DOUBLE,MPI_SUM,world);
if (nextra)
for (i = 0; i < nextra; i++) gg += fextra[i]*fextra[i];
for (i = 0; i < n3; i++) h[i] = g[i] = f[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 < n; niter++) {
for (niter = 0; niter < niter_max; niter++) {
ntimestep = ++update->ntimestep;
// line minimization along direction h from current atom->x
eprevious = ecurrent;
if (ndof) x = atom->x[0];
fail = (this->*linemin)(ndof,x,h,x0,ecurrent,dmax,alpha_final,neval);
fail = (this->*linemin)(ecurrent,alpha_final,neval);
if (fail) return fail;
// function evaluation criterion
@ -90,20 +93,29 @@ int MinCG::iterate(int n)
// force tolerance criterion
if (ndof) f = atom->f[0];
dot[0] = dot[1] = 0.0;
for (i = 0; i < ndof; i++) {
for (i = 0; i < n3; i++) {
dot[0] += f[i]*f[i];
dot[1] += f[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)
for (i = 0; i < nextra; i++) {
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;
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
@ -111,15 +123,26 @@ int MinCG::iterate(int n)
// reinitialize CG every ndof iterations by setting beta = 0.0
beta = MAX(0.0,(dotall[0] - dotall[1])/gg);
if ((niter+1) % ndoftotal == 0) beta = 0.0;
if ((niter+1) % nlimit == 0) beta = 0.0;
gg = dotall[0];
for (i = 0; i < ndof; i++) {
for (i = 0; i < n3; i++) {
g[i] = f[i];
h[i] = g[i] + beta*h[i];
}
if (nextra)
for (i = 0; i < nextra; 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];
}
@ -127,18 +150,30 @@ int MinCG::iterate(int n)
// reinitialize CG if new search direction h is not downhill
dot[0] = 0.0;
for (i = 0; i < ndof; i++) dot[0] += g[i]*h[i];
for (i = 0; i < n3; 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)
for (i = 0; i < nextra; i++)
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 < ndof; i++) h[i] = g[i];
if (nextra)
for (i = 0; i < nextra; i++)
hextra[i] = gextra[i];
for (i = 0; i < n3; 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
@ -152,4 +187,3 @@ int MinCG::iterate(int n)
return MAXITER;
}

4
src/min_cg.h
View File

@ -14,11 +14,11 @@
#ifndef MIN_CG_H
#define MIN_CG_H
#include "min.h"
#include "min_linesearch.h"
namespace LAMMPS_NS {
class MinCG : public Min {
class MinCG : public MinLineSearch {
public:
MinCG(class LAMMPS *);
int iterate(int);

568
src/min_linesearch.cpp Normal file
View File

@ -0,0 +1,568 @@
/* ----------------------------------------------------------------------
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 QUADRATIC_TOL 0.1
#define IDEAL_TOL 1.0e-8
#define EPS_QUAD 1.0e-28
// same as in other min classes
enum{MAXITER,MAXEVAL,ETOL,FTOL,DOWNHILL,ZEROALPHA,ZEROFORCE,ZEROQUAD};
#define MIN(A,B) ((A) < (B)) ? (A) : (B)
#define MAX(A,B) ((A) > (B)) ? (A) : (B)
/* ---------------------------------------------------------------------- */
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];
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];
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]));
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];
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
while (1) {
if (nextra_global) modify->min_step(0.0,hextra);
for (i = 0; i < n3; i++) x[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];
}
if (nextra_global) modify->min_step(alpha,hextra);
for (i = 0; i < n3; i++) x[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];
}
ecurrent = energy_force(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) return 0;
// reduce alpha
alpha *= ALPHA_REDUCE;
// backtracked all the way to 0.0
// reset to starting point, exit with error
if (alpha <= 0.0 || de_ideal >= -IDEAL_TOL) {
if (nextra_global) modify->min_step(0.0,hextra);
for (i = 0; i < n3; i++) x[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];
}
ecurrent = energy_force(0);
nfunc++;
return ZEROALPHA;
}
}
}
/* ----------------------------------------------------------------------
linemin: quadratic line search (adapted from Dennis and Schnabel)
basic idea is to backtrack until change in energy is sufficiently small
based on ENERGY_QUADRATIC, then use a quadratic approximation
using forces at two alpha values to project to minimum
use forces rather than energy change to do projection
this is b/c the forces are going to zero and can become very small
unlike energy differences which are the difference of two finite
values and are thus limited by machine precision
two changes that were critical to making this method work:
a) limit maximum step to alpha <= 1
b) ignore energy criterion if delE <= ENERGY_QUADRATIC
several other ideas also seemed to help:
c) making each step from starting point (alpha = 0), not previous alpha
d) quadratic model based on forces, not energy
e) exiting immediately if f.dir <= 0 (search direction not downhill)
so that CG can restart
a,c,e were also adopted for the backtracking linemin function
------------------------------------------------------------------------- */
int 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;
// 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];
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]));
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];
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
fhprev = fdothall;
engprev = eoriginal;
alphaprev = 0.0;
while (1) {
if (nextra_global) modify->min_step(0.0,hextra);
for (i = 0; i < n3; i++) x[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];
}
if (nextra_global) modify->min_step(alpha,hextra);
for (i = 0; i < n3; i++) x[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];
}
ecurrent = energy_force(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];
}
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) {
if (nextra_global) modify->min_step(0.0,hextra);
for (i = 0; i < n3; i++) x[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];
}
ecurrent = energy_force(0);
nfunc++;
return ZEROQUAD;
}
// check if ready for quadratic projection, equivalent to secant method
// alpha0 = projected alpha
relerr = fabs(1.0+(0.5*alpha*(alpha-alphaprev)*
(fh+fhprev)-ecurrent)/engprev);
alpha0 = alpha - (alpha-alphaprev)*fh/delfh;
if (relerr <= QUADRATIC_TOL && alpha0 > 0.0) {
if (nextra_global) modify->min_step(0.0,hextra);
for (i = 0; i < n3; i++) x[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];
}
if (nextra_global) modify->min_step(alpha0,hextra);
for (i = 0; i < n3; i++) x[i] += alpha0*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] += alpha0*hatom[i];
}
ecurrent = energy_force(1);
nfunc++;
// if backtracking energy change is better than ideal, exit with success
de_ideal = -BACKTRACK_SLOPE*alpha0*fdothall;
de = ecurrent - eoriginal;
if (de <= de_ideal || de_ideal >= -IDEAL_TOL) {
alpha = alpha0;
return 0;
}
// drop back from alpha0 to alpha
if (nextra_global) modify->min_step(0.0,hextra);
for (i = 0; i < n3; i++) x[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];
}
if (nextra_global) modify->min_step(alpha,hextra);
for (i = 0; i < n3; i++) x[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];
}
ecurrent = energy_force(1);
nfunc++;
}
// if backtracking energy change is better than ideal, exit with success
de_ideal = -BACKTRACK_SLOPE*alpha*fdothall;
de = ecurrent - eoriginal;
if (de <= de_ideal) 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) {
if (nextra_global) modify->min_step(0.0,hextra);
for (i = 0; i < n3; i++) x[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];
}
ecurrent = energy_force(0);
nfunc++;
return ZEROALPHA;
}
}
}

53
src/min_linesearch.h Normal file
View File

@ -0,0 +1,53 @@
/* ----------------------------------------------------------------------
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
http://lammps.sandia.gov, Sandia National Laboratories
Steve Plimpton, sjplimp@sandia.gov
Copyright (2003) Sandia Corporation. Under the terms of Contract
DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
certain rights in this software. This software is distributed under
the GNU General Public License.
See the README file in the top-level LAMMPS directory.
------------------------------------------------------------------------- */
#ifndef MIN_LSRCH_H
#define MIN_LSRCH_H
#include "min.h"
namespace LAMMPS_NS {
class MinLineSearch : public Min {
public:
MinLineSearch(class LAMMPS *);
~MinLineSearch();
void init_style();
void setup_style();
void reset_vectors();
protected:
// vectors needed by linesearch minimizers
// allocated and stored by fix_minimize
// x,f are stored by parent or Atom class or Pair class
double *x0; // coords at start of linesearch
double *g; // old gradient vector
double *h; // search direction vector
double *gextra; // g,h for extra global dof, x0 is stored by fix
double *hextra;
double **x0extra_atom; // x0,g,h for extra per-atom dof
double **gextra_atom;
double **hextra_atom;
typedef int (MinLineSearch::*FnPtr)(double, double &, int &);
FnPtr linemin;
int linemin_backtrack(double, double &, int &);
int linemin_quadratic(double, double &, int &);
};
}
#endif

61
src/min_sd.cpp
View File

@ -31,37 +31,41 @@ enum{MAXITER,MAXEVAL,ETOL,FTOL,DOWNHILL,ZEROALPHA,ZEROFORCE,ZEROQUAD};
/* ---------------------------------------------------------------------- */
MinSD::MinSD(LAMMPS *lmp) : Min(lmp) {}
MinSD::MinSD(LAMMPS *lmp) : MinLineSearch(lmp) {}
/* ----------------------------------------------------------------------
minimization via steepest descent
------------------------------------------------------------------------- */
int MinSD::iterate(int n)
int MinSD::iterate(int niter_max)
{
int i,fail,ntimestep;
double dot,dotall;
int i,m,n,fail,ntimestep;
double fdotf;
double *fatom,*hatom;
double *x = NULL;
double *f = NULL;
// initialize working vectors
if (ndof) f = atom->f[0];
for (i = 0; i < ndof; i++) h[i] = f[i];
if (nextra)
for (i = 0; i < nextra; i++)
hextra[i] = fextra[i];
for (i = 0; i < n3; i++) h[i] = f[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 < n; niter++) {
for (niter = 0; niter < niter_max; niter++) {
ntimestep = ++update->ntimestep;
// line minimization along direction h from current atom->x
// line minimization along h from current position x
// h = downhill gradient direction
eprevious = ecurrent;
if (ndof) x = atom->x[0];
fail = (this->*linemin)(ndof,x,h,x0,ecurrent,dmax,alpha_final,neval);
fail = (this->*linemin)(ecurrent,alpha_final,neval);
if (fail) return fail;
// function evaluation criterion
@ -76,22 +80,21 @@ int MinSD::iterate(int n)
// force tolerance criterion
if (ndof) f = atom->f[0];
dot = 0.0;
for (i = 0; i < ndof; i++) dot += f[i]*f[i];
MPI_Allreduce(&dot,&dotall,1,MPI_DOUBLE,MPI_SUM,world);
if (nextra)
for (i = 0; i < nextra; i++)
dotall += fextra[i]*fextra[i];
if (dotall < update->ftol * update->ftol) return FTOL;
fdotf = fnorm_sqr();
if (fdotf < update->ftol*update->ftol) return FTOL;
// set new search direction h to f = -Grad(x)
for (i = 0; i < ndof; i++) h[i] = f[i];
if (nextra)
for (i = 0; i < nextra; i++)
hextra[i] = fextra[i];
for (i = 0; i < n3; i++) h[i] = f[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

4
src/min_sd.h
View File

@ -14,11 +14,11 @@
#ifndef MIN_SD_H
#define MIN_SD_H
#include "min_cg.h"
#include "min_linesearch.h"
namespace LAMMPS_NS {
class MinSD : public Min {
class MinSD : public MinLineSearch {
public:
MinSD(class LAMMPS *);
int iterate(int);

1
src/minimize.cpp
View File

@ -48,6 +48,7 @@ void Minimize::command(int narg, char **arg)
update->whichflag = 1;
lmp->init();
update->minimize->setup();
update->minimize->run();
Finish finish(lmp);

12
src/modify.cpp
View File

@ -428,26 +428,26 @@ void Modify::min_store()
}
/* ----------------------------------------------------------------------
displace extra dof along vector fextra, only for relevant fixes
displace extra dof along vector hextra, only for relevant fixes
------------------------------------------------------------------------- */
void Modify::min_step(double alpha, double *fextra)
void Modify::min_step(double alpha, double *hextra)
{
int ifix,index;
index = 0;
for (int i = 0; i < n_min_energy; i++) {
ifix = list_min_energy[i];
fix[ifix]->min_step(alpha,&fextra[index]);
fix[ifix]->min_step(alpha,&hextra[index]);
index += fix[ifix]->min_dof();
}
}
/* ----------------------------------------------------------------------
compute max allowed step size along vector fextra, only for relevant fixes
compute max allowed step size along vector hextra, only for relevant fixes
------------------------------------------------------------------------- */
double Modify::max_alpha(double *fextra)
double Modify::max_alpha(double *hextra)
{
int ifix,index;
@ -455,7 +455,7 @@ double Modify::max_alpha(double *fextra)
index = 0;
for (int i = 0; i < n_min_energy; i++) {
ifix = list_min_energy[i];
double alpha_one = fix[ifix]->max_alpha(&fextra[index]);
double alpha_one = fix[ifix]->max_alpha(&hextra[index]);
alpha = MIN(alpha,alpha_one);
index += fix[ifix]->min_dof();
}

1
src/modify.h
View File

@ -60,6 +60,7 @@ class Modify : protected Pointers {
void final_integrate_respa(int);
void min_post_force(int);
double min_energy(double *);
void min_store();
void min_step(double, double *);

1
src/pair.h
View File

@ -96,6 +96,7 @@ class Pair : protected Pointers {
virtual void *extract(char *) {return NULL;}
virtual void swap_eam(double *, double **) {}
virtual void reset_dt() {}
virtual void min_pointers(double **, double **) {}
protected:
int allocated; // 0/1 = whether arrays are allocated