From cb90e03944e20f6bc415184dbb4d67c8e3bafa9b Mon Sep 17 00:00:00 2001 From: sjplimp Date: Thu, 5 Mar 2009 20:49:19 +0000 Subject: [PATCH] git-svn-id: svn://svn.icms.temple.edu/lammps-ro/trunk@2622 f3b2605a-c512-4ea7-a41b-209d697bcdaa --- src/compute_pressure.cpp | 9 ++- src/fix.h | 1 + src/fix_box_relax.cpp | 117 +++++++++++++++++++-------------------- src/fix_box_relax.h | 3 +- src/min_cg.cpp | 22 +++++--- src/modify.cpp | 24 ++++++++ src/modify.h | 1 + 7 files changed, 106 insertions(+), 71 deletions(-) diff --git a/src/compute_pressure.cpp b/src/compute_pressure.cpp index f40ac80ad3..05ed1ac232 100644 --- a/src/compute_pressure.cpp +++ b/src/compute_pressure.cpp @@ -64,12 +64,12 @@ ComputePressure::ComputePressure(LAMMPS *lmp, int narg, char **arg) : keflag = 1; pairflag = 1; bondflag = angleflag = dihedralflag = improperflag = 1; - fixflag = kspaceflag = 1; + kspaceflag = fixflag = 1; } else { keflag = 0; pairflag = 0; bondflag = angleflag = dihedralflag = improperflag = 0; - fixflag = kspaceflag = 0; + kspaceflag = fixflag = 0; int iarg = 4; while (iarg < narg) { if (strcmp(arg[iarg],"ke") == 0) keflag = 1; @@ -80,6 +80,11 @@ ComputePressure::ComputePressure(LAMMPS *lmp, int narg, char **arg) : else if (strcmp(arg[iarg],"improper") == 0) improperflag = 1; else if (strcmp(arg[iarg],"kspace") == 0) kspaceflag = 1; else if (strcmp(arg[iarg],"fix") == 0) fixflag = 1; + else if (strcmp(arg[iarg],"virial") == 0) { + pairflag = 1; + bondflag = angleflag = dihedralflag = improperflag = 1; + kspaceflag = fixflag = 1; + } else error->all("Illegal compute pressure command"); iarg++; } diff --git a/src/fix.h b/src/fix.h index 35388650c4..7bf7f144bc 100644 --- a/src/fix.h +++ b/src/fix.h @@ -105,6 +105,7 @@ class Fix : protected Pointers { virtual double min_energy(double *) {return 0.0;} virtual void min_store() {} virtual void min_step(double, double *) {} + virtual double max_alpha(double *) {return 0.0;} virtual int min_dof() {return 0;} virtual int pack_comm(int, int *, double *, int, int *) {return 0;} diff --git a/src/fix_box_relax.cpp b/src/fix_box_relax.cpp index 886a1a8740..03f84e547a 100644 --- a/src/fix_box_relax.cpp +++ b/src/fix_box_relax.cpp @@ -29,6 +29,9 @@ using namespace LAMMPS_NS; +#define MIN(A,B) ((A) < (B)) ? (A) : (B) +#define MAX(A,B) ((A) > (B)) ? (A) : (B) + enum{XYZ,XY,YZ,XZ,ANISO}; /* ---------------------------------------------------------------------- */ @@ -60,18 +63,24 @@ FixBoxRelax::FixBoxRelax(LAMMPS *lmp, int narg, char **arg) : (press_couple == XY || press_couple == YZ || press_couple == XZ)) error->all("Invalid fix box/relax command for a 2d simulation"); - if (strcmp(arg[4],"NULL") == 0) p_flag[0] = 0; - else { + if (strcmp(arg[4],"NULL") == 0) { + p_target[0] = 0.0; + p_flag[0] = 0; + } else { p_target[0] = atof(arg[4]); p_flag[0] = 1; } - if (strcmp(arg[5],"NULL") == 0) p_flag[1] = 0; - else { + if (strcmp(arg[5],"NULL") == 0) { + p_target[1] = 0.0; + p_flag[1] = 0; + } else { p_target[1] = atof(arg[5]); p_flag[1] = 1; } - if (strcmp(arg[6],"NULL") == 0) p_flag[2] = 0; - else { + if (strcmp(arg[6],"NULL") == 0) { + p_target[2] = 0.0; + p_flag[2] = 0; + } else { if (domain->dimension == 2) error->all("Invalid fix box/relax command for a 2d simulation"); p_target[2] = atof(arg[6]); @@ -79,9 +88,12 @@ FixBoxRelax::FixBoxRelax(LAMMPS *lmp, int narg, char **arg) : } } + pflagsum = p_flag[0] + p_flag[1] + p_flag[2]; + // process extra keywords allremap = 0; + smax = 0.0001; int iarg; if (press_couple == XYZ) iarg = 5; @@ -94,6 +106,10 @@ FixBoxRelax::FixBoxRelax(LAMMPS *lmp, int narg, char **arg) : else if (strcmp(arg[iarg+1],"partial") == 0) allremap = 1; else error->all("Illegal fix box/relax command"); iarg += 2; + } else if (strcmp(arg[iarg],"smax") == 0) { + if (iarg+2 > narg) error->all("Illegal fix box/relax command"); + smax = atof(arg[iarg+1]); + iarg += 2; } else error->all("Illegal fix box/relax command"); } @@ -138,7 +154,7 @@ FixBoxRelax::FixBoxRelax(LAMMPS *lmp, int narg, char **arg) : delete [] newarg; tflag = 1; - // create a new compute pressure style + // create a new compute pressure style (virial only) // id = fix-ID + press, compute group = all // pass id_temp as 4th arg to pressure constructor @@ -147,12 +163,13 @@ FixBoxRelax::FixBoxRelax(LAMMPS *lmp, int narg, char **arg) : strcpy(id_press,id); strcat(id_press,"_press"); - newarg = new char*[4]; + newarg = new char*[5]; newarg[0] = id_press; newarg[1] = (char *) "all"; newarg[2] = (char *) "pressure"; newarg[3] = id_temp; - modify->add_compute(4,newarg); + newarg[4] = (char *) "virial"; + modify->add_compute(5,newarg); delete [] newarg; pflag = 1; @@ -261,47 +278,20 @@ double FixBoxRelax::min_energy(double *fextra) fextra[0] = pv2e * (p_current[0] - p_target[0])*2.0*scale*volinit; } - } else if (press_couple == XY) { - scalex = domain->xprd/xprdinit; - scaley = domain->yprd/yprdinit; - scalez = domain->zprd/zprdinit; - eng = pv2e * (p_target[0] + p_flag[2]*p_target[2])/3.0 * - (scalex*scaley*scalez-1.0)*volinit; - scalex = scaley = domain->xprd/xprdinit; - eng = pv2e * p_target[0] * (scalex*scaley-1.0)*volinit; - fextra[0] = pv2e * (p_current[0] - p_target[0])*2.0*scale*volinit; - if (p_flag[2]) { - scalez = domain->zprd/zprdinit; - eng += pv2e * p_target[2] * (scalez-1.0)*volinit; - fextra[1] = pv2e * (p_current[2] - p_target[2])*scalex*scaley*volinit; - } - - } else if (press_couple == YZ) { - - } else if (press_couple == XZ) { - - } else if (press_couple == ANISO) { - scalex = domain->xprd/xprdinit; - scaley = domain->yprd/yprdinit; - scalez = domain->zprd/zprdinit; - if (dimension == 3) { - eng = pv2e * (p_flag[0]*p_target[0] + p_flag[1]*p_target[1] + - p_flag[2]*p_target[2])/3.0 * - (scalex*scaley*scalez-1.0)*volinit; - if (p_flag[0]) - fextra[0] = pv2e * (p_current[0] - p_target[0])*scaley*scalez*volinit; - if (p_flag[1]) - fextra[1] = pv2e * (p_current[1] - p_target[1])*scalex*scalez*volinit; - if (p_flag[2]) - fextra[2] = pv2e * (p_current[2] - p_target[2])*scalex*scaley*volinit; - } else { - eng = pv2e * (p_flag[0]*p_target[0] + p_flag[1]*p_target[1])/2.0 * - (scalex*scaley-1.0)*volinit; - if (p_flag[0]) - fextra[0] = pv2e * (p_current[0] - p_target[0])*scaley*volinit; - if (p_flag[1]) - fextra[1] = pv2e * (p_current[1] - p_target[1])*scalex*volinit; - } + } else { + scalex = scaley = scalez = 1.0; + if (p_flag[0]) scalex = domain->xprd/xprdinit; + if (p_flag[1]) scaley = domain->yprd/yprdinit; + if (p_flag[2]) scalez = domain->zprd/zprdinit; + eng = pv2e * (p_flag[0]*p_target[0] + p_flag[1]*p_target[1] + + p_flag[2]*p_target[2])/pflagsum * + (scalex*scaley*scalez-1.0)*volinit; + if (p_flag[0]) + fextra[0] = pv2e * (p_current[0] - p_target[0])*scaley*scalez*volinit; + if (p_flag[1]) + fextra[1] = pv2e * (p_current[1] - p_target[1])*scalex*scalez*volinit; + if (p_flag[2]) + fextra[2] = pv2e * (p_current[2] - p_target[2])*scalex*scaley*volinit; } return eng; @@ -332,16 +322,7 @@ void FixBoxRelax::min_step(double alpha, double *fextra) { if (press_couple == XYZ) { ds[0] = ds[1] = ds[2] = alpha*fextra[0]; - } else if (press_couple == XY) { - ds[0] = ds[1] = alpha*fextra[0]; - if (p_flag[2]) ds[2] = alpha*fextra[1]; - } else if (press_couple == YZ) { - ds[1] = ds[2] = alpha*fextra[0]; - if (p_flag[0]) ds[0] = alpha*fextra[1]; - } else if (press_couple == XZ) { - ds[0] = ds[2] = alpha*fextra[0]; - if (p_flag[1]) ds[1] = alpha*fextra[1]; - } else if (press_couple == ANISO) { + } 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]; @@ -350,6 +331,22 @@ void FixBoxRelax::min_step(double alpha, double *fextra) remap(); } +/* ---------------------------------------------------------------------- + max allowed step size along fextra +------------------------------------------------------------------------- */ + +double FixBoxRelax::max_alpha(double *fextra) +{ + double alpha = 0.0; + if (press_couple == XYZ) alpha = smax/fabs(fextra[0]); + else { + alpha = smax/fabs(fextra[0]); + alpha = MIN(alpha,smax/fabs(fextra[1])); + alpha = MIN(alpha,smax/fabs(fextra[2])); + } + return alpha; +} + /* ---------------------------------------------------------------------- return number of degrees of freedom added by this fix ------------------------------------------------------------------------- */ diff --git a/src/fix_box_relax.h b/src/fix_box_relax.h index 58367c1640..9fe9caa2f7 100644 --- a/src/fix_box_relax.h +++ b/src/fix_box_relax.h @@ -27,6 +27,7 @@ class FixBoxRelax : public Fix { double min_energy(double *); void min_store(); void min_step(double, double *); + double max_alpha(double *); int min_dof(); private: @@ -36,7 +37,7 @@ class FixBoxRelax : public Fix { double p_target[3],p_current[3]; double dilation[3]; double volinit,xprdinit,yprdinit,zprdinit; - double pv2e; + double smax,pv2e,pflagsum; double boxlo0[3],boxhi0[3]; // box bounds at start of line search double s0[6]; // scale matrix in Voigt notation diff --git a/src/min_cg.cpp b/src/min_cg.cpp index 41982fca91..0607e9f90e 100644 --- a/src/min_cg.cpp +++ b/src/min_cg.cpp @@ -612,7 +612,7 @@ int MinCG::linemin_backtrack(int n, double *x, double *dir, double &alpha, int &nfunc) { int i,m; - double fdotdirall,fdotdirme,hmax,hme,eoriginal; + double fdotdirall,fdotdirme,hmax,hme,alpha_extra,eoriginal; double de_ideal,de; double *f = NULL; @@ -630,17 +630,20 @@ int MinCG::linemin_backtrack(int n, double *x, double *dir, if (fdotdirall <= 0.0) return DOWNHILL; // initial alpha = stepsize to change any atom coord by maxdist - // limit alpha <= 1.0 else backtrack from huge value when forces are tiny + // 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); - if (nextra) + 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; - alpha = MIN(ALPHA_MAX,maxdist/hmax); // store coords and other dof at start of linesearch @@ -711,7 +714,7 @@ int MinCG::linemin_quadratic(int n, double *x, double *dir, double &alpha, int &nfunc) { int i,m; - double fdotdirall,fdotdirme,hmax,hme,alphamax,eoriginal; + double fdotdirall,fdotdirme,hmax,hme,alphamax,alpha_extra,eoriginal; double de_ideal,de; double delfh,engprev,relerr,alphaprev,fhprev,ff,fh,alpha0; double dot[2],dotall[2]; @@ -729,17 +732,20 @@ int MinCG::linemin_quadratic(int n, double *x, double *dir, if (fdotdirall <= 0.0) return DOWNHILL; // initial alpha = stepsize to change any atom coord by maxdist - // limit alpha <= 1.0 else backtrack from huge value when forces are tiny + // 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); - if (nextra) + 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; - alpha = MIN(ALPHA_MAX,maxdist/hmax); // store coords and other dof at start of linesearch diff --git a/src/modify.cpp b/src/modify.cpp index eb07d360be..4d89e8f826 100644 --- a/src/modify.cpp +++ b/src/modify.cpp @@ -53,6 +53,11 @@ using namespace LAMMPS_NS; #define MIN_POST_FORCE 16384 #define MIN_ENERGY 32768 +#define MIN(A,B) ((A) < (B)) ? (A) : (B) +#define MAX(A,B) ((A) > (B)) ? (A) : (B) + +#define BIG 1.0e20 + /* ---------------------------------------------------------------------- */ Modify::Modify(LAMMPS *lmp) : Pointers(lmp) @@ -435,6 +440,25 @@ void Modify::min_step(double alpha, double *fextra) } } +/* ---------------------------------------------------------------------- + compute max allowed step size along vector fextra, only for relevant fixes +------------------------------------------------------------------------- */ + +double Modify::max_alpha(double *fextra) +{ + int ifix,index; + + double alpha = BIG; + 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]); + alpha = MIN(alpha,alpha_one); + index += fix[ifix]->min_dof(); + } + return alpha; +} + /* ---------------------------------------------------------------------- extract extra dof for minimization, only for relevant fixes ------------------------------------------------------------------------- */ diff --git a/src/modify.h b/src/modify.h index 2e4a44a18c..a749f010dc 100644 --- a/src/modify.h +++ b/src/modify.h @@ -63,6 +63,7 @@ class Modify : protected Pointers { double min_energy(double *); void min_store(); void min_step(double, double *); + double max_alpha(double *); int min_dof(); void add_fix(int, char **);