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@2622 f3b2605a-c512-4ea7-a41b-209d697bcdaa

Browse Source
This commit is contained in:
sjplimp
2009-03-05 20:49:19 +00:00
parent 12f25b8829
commit cb90e03944
7 changed files with 106 additions and 71 deletions
Download Patch File Download Diff File Expand all files Collapse all files

9
src/compute_pressure.cpp
View File

@ -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++;
}

1
src/fix.h
View File

@ -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;}

101
src/fix_box_relax.cpp
View File

@ -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,32 +278,13 @@ 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) {
} 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])/3.0 *
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;
@ -294,14 +292,6 @@ double FixBoxRelax::min_energy(double *fextra)
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;
}
}
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
------------------------------------------------------------------------- */

3
src/fix_box_relax.h
View File

@ -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

22
src/min_cg.cpp
View File

@ -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

24
src/modify.cpp
View File

@ -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
------------------------------------------------------------------------- */

1
src/modify.h
View File

@ -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 **);