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

Browse Source
This commit is contained in:
sjplimp
2010-06-02 15:41:31 +00:00
parent c00c60c86d
commit b3132ab259
13 changed files with 134 additions and 81 deletions
Download Patch File Download Diff File Expand all files Collapse all files

4
src/Makefile.lib
View File

File diff suppressed because one or more lines are too long

7
src/fix_ave_time.cpp
View File

@ -124,7 +124,7 @@ FixAveTime::FixAveTime(LAMMPS *lmp, int narg, char **arg) :
ids[nvalues+icol] = new char[n];
strcpy(ids[nvalues+icol],ids[nvalues]);
}
nvalues += ncols;
nvalues += ncols-1;
}
} else if (mode == VECTOR && which[nvalues] == FIX &&
@ -144,10 +144,11 @@ FixAveTime::FixAveTime(LAMMPS *lmp, int narg, char **arg) :
ids[nvalues+icol] = new char[n];
strcpy(ids[nvalues+icol],ids[nvalues]);
}
nvalues += ncols;
nvalues += ncols-1;
}
}
} else nvalues++;
nvalues++;
iarg++;
} else break;
}

2
src/fix_box_relax.cpp
View File

@ -892,7 +892,7 @@ void FixBoxRelax::compute_press_target()
if (pflagsum) p_hydro /= pflagsum;
for (int i = 0; i < 3; i++) {
if (p_flag[i] && fabs(p_hydro - p_target[i] > 1.0e-6)) deviatoric_flag = 1;
if (p_flag[i] && fabs(p_hydro - p_target[i]) > 1.0e-6) deviatoric_flag = 1;
}
if (pstyle == TRICLINIC) {

16
src/fix_viscous.cpp
View File

@ -63,6 +63,7 @@ int FixViscous::setmask()
int mask = 0;
mask |= POST_FORCE;
mask |= POST_FORCE_RESPA;
mask |= MIN_POST_FORCE;
return mask;
}
@ -89,6 +90,13 @@ void FixViscous::setup(int vflag)
/* ---------------------------------------------------------------------- */
void FixViscous::min_setup(int vflag)
{
post_force(vflag);
}
/* ---------------------------------------------------------------------- */
void FixViscous::post_force(int vflag)
{
// apply drag force to atoms in group
@ -118,3 +126,11 @@ void FixViscous::post_force_respa(int vflag, int ilevel, int iloop)
{
if (ilevel == nlevels_respa-1) post_force(vflag);
}
/* ---------------------------------------------------------------------- */
void FixViscous::min_post_force(int vflag)
{
post_force(vflag);
}

2
src/fix_viscous.h
View File

@ -31,8 +31,10 @@ class FixViscous : public Fix {
int setmask();
void init();
void setup(int);
void min_setup(int);
void post_force(int);
void post_force_respa(int, int, int);
void min_post_force(int);
private:
double *gamma;

118
src/min.cpp
View File

@ -158,6 +158,8 @@ void Min::init()
neighbor->delay = 0;
neighbor->dist_check = 1;
niter = neval = 0;
// style-specific initialization
init_style();
@ -216,6 +218,10 @@ void Min::setup()
neighbor->build();
neighbor->ncalls = 0;
// atoms may have migrated in comm->exchange()
reset_vectors();
// compute all forces
ev_set(update->ntimestep);
@ -240,9 +246,15 @@ void Min::setup()
modify->setup(vflag);
output->setup(1);
// atoms may have migrated in comm->exchange()
// stats for Finish to print
reset_vectors();
ecurrent = pe_compute->compute_scalar();
if (nextra_global) ecurrent += modify->min_energy(fextra);
if (output->thermo->normflag) ecurrent /= atom->natoms;
einitial = ecurrent;
fnorm2_init = sqrt(fnorm_sqr());
fnorminf_init = fnorm_inf();
}
/* ----------------------------------------------------------------------
@ -270,6 +282,10 @@ void Min::setup_minimal(int flag)
neighbor->ncalls = 0;
}
// atoms may have migrated in comm->exchange()
reset_vectors();
// compute all forces
ev_set(update->ntimestep);
@ -293,30 +309,6 @@ void Min::setup_minimal(int flag)
modify->setup(vflag);
// atoms may have migrated in comm->exchange()
reset_vectors();
}
/* ----------------------------------------------------------------------
perform minimization, calling iterate() for nsteps
------------------------------------------------------------------------- */
void Min::run(int nsteps)
{
// possible stop conditions
char *stopstrings[] = {"max iterations",
"max force evaluations",
"energy tolerance",
"force tolerance",
"search direction is not downhill",
"linesearch alpha is zero",
"forces are zero",
"quadratic factors are zero",
"trust region too small",
"HFTN minimizer error"};
// stats for Finish to print
ecurrent = pe_compute->compute_scalar();
@ -326,22 +318,49 @@ void Min::run(int nsteps)
einitial = ecurrent;
fnorm2_init = sqrt(fnorm_sqr());
fnorminf_init = fnorm_inf();
}
/* ----------------------------------------------------------------------
perform minimization, calling iterate() for N steps
complete = 1 if caller wants to force N iterations (in every replica)
minimize command and PRD call with complete = 0
NEB calls with complete = 1
------------------------------------------------------------------------- */
void Min::run(int n, int complete)
{
// minimizer iterations
int stop_condition = iterate(nsteps);
stopstr = stopstrings[stop_condition];
int iter_start = niter;
stop_condition = iterate(n);
stopstr = stopstrings(stop_condition);
// account for early exit from iterate loop due to convergence
// set niter/nsteps for Finish stats to print
// if early exit from iterate loop and complete flag set
// perform remaining dummy iterations
if (stop_condition && complete) {
int ntimestep;
while (niter - iter_start < n) {
ntimestep = ++update->ntimestep;
niter++;
ecurrent = energy_force(0);
if (output->next == ntimestep) {
timer->stamp();
output->write(ntimestep);
timer->stamp(TIME_OUTPUT);
}
}
}
// if early exit from iterate loop and complete flag not set
// set update->nsteps to niter for Finish stats to print
// set output->next values to this timestep
// call engergy_force() to insure vflag is set when forces computed
// call energy_force() to insure vflag is set when forces computed
// output->write does final output for thermo, dump, restart files
// add ntimestep to all computes that store invocation times
// since are hardwireing call to thermo/dumps and computes may not be ready
// since are hardwiring call to thermo/dumps and computes may not be ready
if (niter < nsteps) {
niter++;
if (stop_condition && !complete) {
update->nsteps = niter;
if (update->restrict_output == 0) {
@ -356,18 +375,18 @@ void Min::run(int nsteps)
ecurrent = energy_force(0);
output->write(update->ntimestep);
}
// stats for Finish to print
efinal = ecurrent;
fnorm2_final = sqrt(fnorm_sqr());
fnorminf_final = fnorm_inf();
}
/* ---------------------------------------------------------------------- */
void Min::cleanup()
{
// stats for Finish to print
efinal = ecurrent;
fnorm2_final = sqrt(fnorm_sqr());
fnorminf_final = fnorm_inf();
// reset reneighboring criteria
neighbor->every = neigh_every;
@ -690,3 +709,22 @@ double Min::fnorm_inf()
return norm_inf;
}
/* ----------------------------------------------------------------------
possible stop conditions
------------------------------------------------------------------------- */
char *Min::stopstrings(int n)
{
char *strings[] = {"max iterations",
"max force evaluations",
"energy tolerance",
"force tolerance",
"search direction is not downhill",
"linesearch alpha is zero",
"forces are zero",
"quadratic factors are zero",
"trust region too small",
"HFTN minimizer error"};
return strings[n];
}

4
src/min.h
View File

@ -24,6 +24,7 @@ class Min : protected Pointers {
double fnorm2_init,fnorminf_init,fnorm2_final,fnorminf_final;
double alpha_final;
int niter,neval;
int stop_condition;
char *stopstr;
Min(class LAMMPS *);
@ -31,7 +32,7 @@ class Min : protected Pointers {
void init();
void setup();
void setup_minimal(int);
void run(int);
void run(int,int);
void cleanup();
void request(class Pair *, int, double);
double memory_usage() {return 0.0;}
@ -95,6 +96,7 @@ class Min : protected Pointers {
double fnorm_sqr();
double fnorm_inf();
char *stopstrings(int);
};
}

9
src/min_cg.cpp
View File

@ -43,7 +43,7 @@ MinCG::MinCG(LAMMPS *lmp) : MinLineSearch(lmp) {}
minimization via conjugate gradient iterations
------------------------------------------------------------------------- */
int MinCG::iterate(int niter_max)
int MinCG::iterate(int maxiter)
{
int i,m,n,fail,ntimestep;
double beta,gg,dot[2],dotall[2];
@ -70,15 +70,14 @@ int MinCG::iterate(int niter_max)
gg = fnorm_sqr();
neval = 0;
for (niter = 0; niter < niter_max; niter++) {
for (int iter = 0; iter < maxiter; iter++) {
ntimestep = ++update->ntimestep;
niter++;
// line minimization along direction h from current atom->x
eprevious = ecurrent;
fail = (this->*linemin)(ecurrent,alpha_final,neval);
fail = (this->*linemin)(ecurrent,alpha_final);
if (fail) return fail;
// function evaluation criterion

2
src/min_hftn.cpp
View File

@ -241,8 +241,6 @@ int MinHFTN::execute_hftn_(const bool bPrintProgress,
//---- DEFINE OUTPUTS PRINTED BY "Finish".
eprevious = dInitialEnergy;
alpha_final = 0.0;
neval = 0;
niter = 0;
dFinalEnergy = dInitialEnergy;
dFinalForce2 = dInitialForce2;

32
src/min_linesearch.cpp
View File

@ -162,7 +162,7 @@ void MinLineSearch::reset_vectors()
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
update neval counter of eng/force function evaluations
output extra: if fail, energy_force() of original x
if succeed, energy_force() at x + alpha*h
atom->x = coords at new configuration
@ -176,8 +176,7 @@ void MinLineSearch::reset_vectors()
start at maxdist, backtrack until energy decrease is sufficient
------------------------------------------------------------------------- */
int MinLineSearch::linemin_backtrack(double eoriginal, double &alpha,
int &nfunc)
int MinLineSearch::linemin_backtrack(double eoriginal, double &alpha)
{
int i,m,n;
double fdothall,fdothme,hme,hmax,hmaxall;
@ -248,16 +247,16 @@ int MinLineSearch::linemin_backtrack(double eoriginal, double &alpha,
// Important diagnostic: test the gradient against energy
// double etmp;
// double alphatmp = alphamax*1.0e-4;
// etmp = alpha_step(alphatmp,1,nfunc);
// etmp = alpha_step(alphatmp,1);
// printf("alpha = %g dele = %g dele_force = %g err = %g\n",
// alphatmp,etmp-eoriginal,-alphatmp*fdothall,
// etmp-eoriginal+alphatmp*fdothall);
// alpha_step(0.0,1,nfunc);
// alpha_step(0.0,1);
// backtrack with alpha until energy decrease is sufficient
while (1) {
ecurrent = alpha_step(alpha,1,nfunc);
ecurrent = alpha_step(alpha,1);
// if energy change is better than ideal, exit with success
@ -279,7 +278,7 @@ int MinLineSearch::linemin_backtrack(double eoriginal, double &alpha,
// reset to starting point, exit with error
if (alpha <= 0.0 || de_ideal >= -IDEAL_TOL) {
ecurrent = alpha_step(0.0,0,nfunc);
ecurrent = alpha_step(0.0,0);
return ZEROALPHA;
}
}
@ -323,8 +322,7 @@ int MinLineSearch::linemin_backtrack(double eoriginal, double &alpha,
//
------------------------------------------------------------------------- */
int MinLineSearch::linemin_quadratic(double eoriginal, double &alpha,
int &nfunc)
int MinLineSearch::linemin_quadratic(double eoriginal, double &alpha)
{
int i,m,n;
double fdothall,fdothme,hme,hmax,hmaxall;
@ -407,14 +405,14 @@ int MinLineSearch::linemin_quadratic(double eoriginal, double &alpha,
// Important diagnostic: test the gradient against energy
// double etmp;
// double alphatmp = alphamax*1.0e-4;
// etmp = alpha_step(alphatmp,1,nfunc);
// etmp = alpha_step(alphatmp,1);
// printf("alpha = %g dele = %g dele_force = %g err = %g\n",
// alphatmp,etmp-eoriginal,-alphatmp*fdothall,
// etmp-eoriginal+alphatmp*fdothall);
// alpha_step(0.0,1,nfunc);
// alpha_step(0.0,1);
while (1) {
ecurrent = alpha_step(alpha,1,nfunc);
ecurrent = alpha_step(alpha,1);
// compute new fh, alpha, delfh
@ -452,7 +450,7 @@ int MinLineSearch::linemin_quadratic(double eoriginal, double &alpha,
// if fh or delfh is epsilon, reset to starting point, exit with error
if (fabs(fh) < EPS_QUAD || fabs(delfh) < EPS_QUAD) {
ecurrent = alpha_step(0.0,0,nfunc);
ecurrent = alpha_step(0.0,0);
return ZEROQUAD;
}
@ -463,7 +461,7 @@ int MinLineSearch::linemin_quadratic(double eoriginal, double &alpha,
alpha0 = alpha - (alpha-alphaprev)*fh/delfh;
if (relerr <= QUADRATIC_TOL && alpha0 > 0.0 && alpha0 < alphamax) {
ecurrent = alpha_step(alpha0,1,nfunc);
ecurrent = alpha_step(alpha0,1);
if (ecurrent < eoriginal) {
if (nextra_global) {
int itmp = modify->min_reset_ref();
@ -500,7 +498,7 @@ int MinLineSearch::linemin_quadratic(double eoriginal, double &alpha,
// reset to starting point, exit with error
if (alpha <= 0.0 || de_ideal >= -IDEAL_TOL) {
ecurrent = alpha_step(0.0,0,nfunc);
ecurrent = alpha_step(0.0,0);
return ZEROALPHA;
}
}
@ -508,7 +506,7 @@ int MinLineSearch::linemin_quadratic(double eoriginal, double &alpha,
/* ---------------------------------------------------------------------- */
double MinLineSearch::alpha_step(double alpha, int resetflag, int &nfunc)
double MinLineSearch::alpha_step(double alpha, int resetflag)
{
int i,n,m;
double *xatom,*x0atom,*hatom;
@ -541,6 +539,6 @@ double MinLineSearch::alpha_step(double alpha, int resetflag, int &nfunc)
// compute and return new energy
nfunc++;
neval++;
return energy_force(resetflag);
}

8
src/min_linesearch.h
View File

@ -42,12 +42,12 @@ class MinLineSearch : public Min {
double **gextra_atom;
double **hextra_atom;
typedef int (MinLineSearch::*FnPtr)(double, double &, int &);
typedef int (MinLineSearch::*FnPtr)(double, double &);
FnPtr linemin;
int linemin_backtrack(double, double &, int &);
int linemin_quadratic(double, double &, int &);
int linemin_backtrack(double, double &);
int linemin_quadratic(double, double &);
double alpha_step(double, int, int &);
double alpha_step(double, int);
};
}

9
src/min_sd.cpp
View File

@ -37,7 +37,7 @@ MinSD::MinSD(LAMMPS *lmp) : MinLineSearch(lmp) {}
minimization via steepest descent
------------------------------------------------------------------------- */
int MinSD::iterate(int niter_max)
int MinSD::iterate(int maxiter)
{
int i,m,n,fail,ntimestep;
double fdotf;
@ -56,16 +56,15 @@ int MinSD::iterate(int niter_max)
if (nextra_global)
for (i = 0; i < nextra_global; i++) hextra[i] = fextra[i];
neval = 0;
for (niter = 0; niter < niter_max; niter++) {
for (int iter = 0; iter < maxiter; iter++) {
ntimestep = ++update->ntimestep;
niter++;
// line minimization along h from current position x
// h = downhill gradient direction
eprevious = ecurrent;
fail = (this->*linemin)(ecurrent,alpha_final,neval);
fail = (this->*linemin)(ecurrent,alpha_final);
if (fail) return fail;
// function evaluation criterion

2
src/minimize.cpp
View File

@ -51,7 +51,7 @@ void Minimize::command(int narg, char **arg)
update->minimize->setup();
timer->barrier_start(TIME_LOOP);
update->minimize->run(update->nsteps);
update->minimize->run(update->nsteps,0);
timer->barrier_stop(TIME_LOOP);
update->minimize->cleanup();