diff -Naur lammps-16Aug09/src/MAKE/Makefile.linux lammps-17Aug09/src/MAKE/Makefile.linux --- lammps-16Aug09/src/MAKE/Makefile.linux 2009-08-13 10:57:57.000000000 -0600 +++ lammps-17Aug09/src/MAKE/Makefile.linux 2009-08-14 14:35:57.000000000 -0600 @@ -14,7 +14,7 @@ DEPFLAGS = -M LINK = icc LINKFORT = -L/opt/intel/fce/10.0.023/lib -#LINKGPU = -L/usr/local/cuda/lib +#LINKGPU = -L/usr/local/cuda/lib64// LINKFLAGS = -O $(PKGPATH) $(LINKFORT) $(LINKGPU) USRLIB = -lfftw -lmpich $(PKGLIB) FORTLIB = -lifcore -lsvml -lompstub -limf diff -Naur lammps-16Aug09/src/Makefile.package.empty lammps-17Aug09/src/Makefile.package.empty --- lammps-16Aug09/src/Makefile.package.empty 1969-12-31 17:00:00.000000000 -0700 +++ lammps-17Aug09/src/Makefile.package.empty 2009-08-13 16:24:56.000000000 -0600 @@ -0,0 +1,6 @@ +# Settings for libraries used by specific LAMMPS packages +# this file is auto-edited when those packages are included/excluded + +PKGINC = +PKGPATH = +PKGLIB = diff -Naur lammps-16Aug09/src/compute_pressure.cpp lammps-17Aug09/src/compute_pressure.cpp --- lammps-16Aug09/src/compute_pressure.cpp 2009-07-02 14:53:24.000000000 -0600 +++ lammps-17Aug09/src/compute_pressure.cpp 2009-08-14 14:54:10.000000000 -0600 @@ -236,10 +236,12 @@ 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; } } } diff -Naur lammps-16Aug09/src/fix.h lammps-17Aug09/src/fix.h --- lammps-16Aug09/src/fix.h 2009-08-06 16:21:38.000000000 -0600 +++ lammps-17Aug09/src/fix.h 2009-08-14 14:54:10.000000000 -0600 @@ -105,6 +105,7 @@ 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 *) {} diff -Naur lammps-16Aug09/src/fix_box_relax.cpp lammps-17Aug09/src/fix_box_relax.cpp --- lammps-16Aug09/src/fix_box_relax.cpp 2009-07-02 13:47:06.000000000 -0600 +++ lammps-17Aug09/src/fix_box_relax.cpp 2009-08-14 14:54:10.000000000 -0600 @@ -318,34 +318,34 @@ } /* ---------------------------------------------------------------------- - 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; } diff -Naur lammps-16Aug09/src/fix_box_relax.h lammps-17Aug09/src/fix_box_relax.h --- lammps-16Aug09/src/fix_box_relax.h 2009-03-25 08:59:54.000000000 -0600 +++ lammps-17Aug09/src/fix_box_relax.h 2009-08-14 14:54:10.000000000 -0600 @@ -24,11 +24,13 @@ ~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: diff -Naur lammps-16Aug09/src/fix_minimize.cpp lammps-17Aug09/src/fix_minimize.cpp --- lammps-16Aug09/src/fix_minimize.cpp 2009-02-25 16:58:04.000000000 -0700 +++ lammps-17Aug09/src/fix_minimize.cpp 2009-08-14 14:54:10.000000000 -0600 @@ -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 @@ 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 @@ 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); } /* ---------------------------------------------------------------------- */ @@ -56,12 +58,119 @@ } /* ---------------------------------------------------------------------- + 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 @@ 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::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 @@ 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::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; } diff -Naur lammps-16Aug09/src/fix_minimize.h lammps-17Aug09/src/fix_minimize.h --- lammps-16Aug09/src/fix_minimize.h 2009-02-25 16:58:04.000000000 -0700 +++ lammps-17Aug09/src/fix_minimize.h 2009-08-14 14:54:10.000000000 -0600 @@ -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 @@ 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(); }; } diff -Naur lammps-16Aug09/src/lammps.cpp lammps-17Aug09/src/lammps.cpp --- lammps-16Aug09/src/lammps.cpp 2008-12-12 17:47:23.000000000 -0700 +++ lammps-17Aug09/src/lammps.cpp 2009-08-14 14:54:10.000000000 -0600 @@ -304,7 +304,7 @@ 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 diff -Naur lammps-16Aug09/src/min.cpp lammps-17Aug09/src/min.cpp --- lammps-16Aug09/src/min.cpp 2009-06-04 11:53:33.000000000 -0600 +++ lammps-17Aug09/src/min.cpp 2009-08-14 14:54:10.000000000 -0600 @@ -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 @@ 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 @@ 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 @@ 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 @@ 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 @@ 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 @@ { 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 @@ if (triclinic) domain->lamda2x(atom->nlocal+atom->nghost); neighbor->build(); neighbor->ncalls = 0; - setup_vectors(); // compute all forces @@ -202,6 +234,10 @@ modify->setup(vflag); output->setup(1); + + // atoms may have migrated in comm->exchange() + + reset_vectors(); } /* ---------------------------------------------------------------------- @@ -210,9 +246,6 @@ void Min::run() { - int i; - double tmp,*f; - // possible stop conditions char *stopstrings[] = {"max iterations","max force evaluations", @@ -221,57 +254,21 @@ "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]; + // stats for Finish to print + ecurrent = pe_compute->compute_scalar(); - if (nextra) ecurrent += modify->min_energy(fextra); + if (nextra_global) ecurrent += modify->min_energy(fextra); if (output->thermo->normflag) ecurrent /= atom->natoms; - - // stats for Finish to print einitial = ecurrent; - - f = NULL; - if (ndof) f = atom->f[0]; - tmp = 0.0; - for (i = 0; i < ndof; i++) tmp += f[i]*f[i]; - MPI_Allreduce(&tmp,&fnorm2_init,1,MPI_DOUBLE,MPI_SUM,world); - if (nextra) - for (i = 0; i < nextra; i++) fnorm2_init += fextra[i]*fextra[i]; - fnorm2_init = sqrt(fnorm2_init); - - tmp = 0.0; - for (i = 0; i < ndof; i++) tmp = MAX(fabs(f[i]),tmp); - MPI_Allreduce(&tmp,&fnorminf_init,1,MPI_DOUBLE,MPI_MAX,world); - if (nextra) - for (i = 0; i < nextra; i++) - fnorminf_init = MAX(fabs(fextra[i]),fnorminf_init); + fnorm2_init = sqrt(fnorm_sqr()); + fnorminf_init = fnorm_inf(); // minimizer iterations @@ -298,16 +295,13 @@ 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,40 +314,23 @@ // 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(); if (nflag == 0) { @@ -375,34 +352,6 @@ 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 @@ // 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; -} - -/* ---------------------------------------------------------------------- - set ndof and vector pointers after atoms have migrated -------------------------------------------------------------------------- */ + if (nflag) { + if (resetflag) fix_minimize->reset_coords(); + reset_vectors(); + } -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 @@ } /* ---------------------------------------------------------------------- - line minimization methods - find minimum-energy starting at x along dir direction - input: n = # of degrees of freedom on this proc - x = ptr to atom->x[0] as vector - dir = search direction as vector - x0 = ptr to fix->x0[0] as vector, for storing initial coords - eoriginal = energy at initial x - maxdist = max distance to move any atom coord - output: return 0 if successful move, non-zero alpha - return non-zero if failed - alpha = distance moved along dir to set x to minimun eng config - caller has several quantities set via last call to eng_force() - must insure last call to eng_force() is consistent with returns - if fail, eng_force() of original x - if succeed, eng_force() at x + alpha*dir - atom->x = coords at new configuration - atom->f = force (-Grad) is evaulated at new configuration - ecurrent = energy of new configuration - NOTE: when call eng_force: n,x,dir,x0,eng may change due to atom migration - updated values are returned by eng_force() - b/c of migration, linemin routines CANNOT store atom-based quantities -------------------------------------------------------------------------- */ - -/* ---------------------------------------------------------------------- - linemin: backtracking line search (Proc 3.1, p 41 in Nocedal and Wright) - uses no gradient info, but should be very robust - start at maxdist, backtrack until energy decrease is sufficient -------------------------------------------------------------------------- */ - -int Min::linemin_backtrack(int n, double *x, double *dir, - double *x0, double eoriginal, double maxdist, - double &alpha, int &nfunc) -{ - int i,m; - double fdotdirall,fdotdirme,hmax,hme,alpha_extra; - double eng,de_ideal,de; - - double *f = NULL; - if (n) f = atom->f[0]; - - // fdotdirall = projection of search dir along downhill gradient - // if search direction is not downhill, exit with error - - fdotdirme = 0.0; - for (i = 0; i < n; i++) fdotdirme += f[i]*dir[i]; - MPI_Allreduce(&fdotdirme,&fdotdirall,1,MPI_DOUBLE,MPI_SUM,world); - if (nextra) - for (i = 0; i < nextra; i++) fdotdirall += fextra[i]*hextra[i]; - if (output->thermo->normflag) fdotdirall /= atom->natoms; - if (fdotdirall <= 0.0) return DOWNHILL; - - // initial alpha = stepsize to change any atom coord by maxdist - // alpha <= ALPHA_MAX, else backtrack from huge value when forces are tiny - // if all search dir components are already 0.0, exit with error - - hme = 0.0; - for (i = 0; i < n; i++) hme = MAX(hme,fabs(dir[i])); - MPI_Allreduce(&hme,&hmax,1,MPI_DOUBLE,MPI_MAX,world); - alpha = MIN(ALPHA_MAX,maxdist/hmax); - if (nextra) { - double alpha_extra = modify->max_alpha(hextra); - alpha = MIN(alpha,alpha_extra); - for (i = 0; i < nextra; i++) - hmax = MAX(hmax,fabs(hextra[i])); - } - if (hmax == 0.0) return ZEROFORCE; - - // store coords and other dof at start of linesearch - - 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 + clear force on own & ghost atoms + setup and clear other arrays as needed ------------------------------------------------------------------------- */ -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; - } - } +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++; } /* ---------------------------------------------------------------------- */ @@ -882,45 +566,62 @@ } /* ---------------------------------------------------------------------- - 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]; - - boxhi0[0] = domain->boxhi[0]; - boxhi0[1] = domain->boxhi[1]; - boxhi0[2] = domain->boxhi[2]; + int i,n; + double *fatom; + + 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; + + 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); + } + } + + 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); - 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; - - 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; + return norm_inf; } diff -Naur lammps-16Aug09/src/min.h lammps-17Aug09/src/min.h --- lammps-16Aug09/src/min.h 2009-06-04 11:53:33.000000000 -0600 +++ lammps-17Aug09/src/min.h 2009-08-14 14:54:10.000000000 -0600 @@ -29,17 +29,23 @@ Min(class LAMMPS *); virtual ~Min(); void init(); + void setup(); void run(); - 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; + void request(class Pair *, int, double); + double memory_usage() {return 0.0;} + void modify_params(int, char **); + 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 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 - int nextra; // extra dof due to fixes - double *fextra; // vectors for extra dof - double *gextra; - double *hextra; + double ndoftotal; // total dof for entire problem - double boxlo0[3]; // box size at start of linesearch - double boxhi0[3]; + 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 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 - // ptr to linemin functions + int neigh_every,neigh_delay,neigh_dist_check; // neighboring params - void setup(); - void eng_force(int *, double **, double **, double **, double *, int); - void setup_vectors(); + 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(); }; } diff -Naur lammps-16Aug09/src/min_cg.cpp lammps-17Aug09/src/min_cg.cpp --- lammps-16Aug09/src/min_cg.cpp 2009-04-03 13:46:03.000000000 -0600 +++ lammps-17Aug09/src/min_cg.cpp 2009-08-14 14:54:10.000000000 -0600 @@ -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 @@ /* ---------------------------------------------------------------------- */ -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; - - 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]; - - 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]; + int nlimit = static_cast (MIN(MAXATOMS,ndoftotal)); - neval = 0; + // initialize working vectors + + 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(); - for (niter = 0; niter < n; niter++) { + neval = 0; + 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 @@ // 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 @@ // 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 @@ // 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 @@ return MAXITER; } - diff -Naur lammps-16Aug09/src/min_cg.h lammps-17Aug09/src/min_cg.h --- lammps-16Aug09/src/min_cg.h 2009-04-03 13:46:03.000000000 -0600 +++ lammps-17Aug09/src/min_cg.h 2009-08-14 14:54:10.000000000 -0600 @@ -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); diff -Naur lammps-16Aug09/src/min_linesearch.cpp lammps-17Aug09/src/min_linesearch.cpp --- lammps-16Aug09/src/min_linesearch.cpp 1969-12-31 17:00:00.000000000 -0700 +++ lammps-17Aug09/src/min_linesearch.cpp 2009-08-14 14:54:10.000000000 -0600 @@ -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; + } + } +} diff -Naur lammps-16Aug09/src/min_linesearch.h lammps-17Aug09/src/min_linesearch.h --- lammps-16Aug09/src/min_linesearch.h 1969-12-31 17:00:00.000000000 -0700 +++ lammps-17Aug09/src/min_linesearch.h 2009-08-14 14:54:10.000000000 -0600 @@ -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 diff -Naur lammps-16Aug09/src/min_sd.cpp lammps-17Aug09/src/min_sd.cpp --- lammps-16Aug09/src/min_sd.cpp 2009-04-03 13:46:03.000000000 -0600 +++ lammps-17Aug09/src/min_sd.cpp 2009-08-14 14:54:10.000000000 -0600 @@ -31,37 +31,41 @@ /* ---------------------------------------------------------------------- */ -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; - - double *x = NULL; - double *f = NULL; - - 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]; + int i,m,n,fail,ntimestep; + double fdotf; + double *fatom,*hatom; + + // initialize working vectors + + 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 @@ // 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 diff -Naur lammps-16Aug09/src/min_sd.h lammps-17Aug09/src/min_sd.h --- lammps-16Aug09/src/min_sd.h 2009-04-03 13:46:03.000000000 -0600 +++ lammps-17Aug09/src/min_sd.h 2009-08-14 14:54:10.000000000 -0600 @@ -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); diff -Naur lammps-16Aug09/src/minimize.cpp lammps-17Aug09/src/minimize.cpp --- lammps-16Aug09/src/minimize.cpp 2008-04-14 14:50:48.000000000 -0600 +++ lammps-17Aug09/src/minimize.cpp 2009-08-14 14:54:10.000000000 -0600 @@ -48,6 +48,7 @@ update->whichflag = 1; lmp->init(); + update->minimize->setup(); update->minimize->run(); Finish finish(lmp); diff -Naur lammps-16Aug09/src/modify.cpp lammps-17Aug09/src/modify.cpp --- lammps-16Aug09/src/modify.cpp 2009-08-06 16:21:38.000000000 -0600 +++ lammps-17Aug09/src/modify.cpp 2009-08-14 14:54:10.000000000 -0600 @@ -428,26 +428,26 @@ } /* ---------------------------------------------------------------------- - 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 @@ 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(); } diff -Naur lammps-16Aug09/src/modify.h lammps-17Aug09/src/modify.h --- lammps-16Aug09/src/modify.h 2009-08-06 16:21:38.000000000 -0600 +++ lammps-17Aug09/src/modify.h 2009-08-14 14:54:10.000000000 -0600 @@ -60,6 +60,7 @@ void final_integrate_respa(int); void min_post_force(int); + double min_energy(double *); void min_store(); void min_step(double, double *); diff -Naur lammps-16Aug09/src/pair.h lammps-17Aug09/src/pair.h --- lammps-16Aug09/src/pair.h 2009-06-08 12:29:00.000000000 -0600 +++ lammps-17Aug09/src/pair.h 2009-08-14 14:54:10.000000000 -0600 @@ -96,6 +96,7 @@ 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