diff --git a/src/SPIN/min_spin.cpp b/src/SPIN/min_spin.cpp index f56c9f0d96..d229927c29 100644 --- a/src/SPIN/min_spin.cpp +++ b/src/SPIN/min_spin.cpp @@ -119,7 +119,7 @@ void MinSpin::reset_vectors() int MinSpin::iterate(int maxiter) { bigint ntimestep; - double fmdotfm; + double fmdotfm,fmsq,fmsqall; int flag,flagall; for (int iter = 0; iter < maxiter; iter++) { @@ -166,8 +166,20 @@ int MinSpin::iterate(int maxiter) // magnetic torque tolerance criterion // sync across replicas if running multi-replica minimization + fmdotfm = fmsq = fmsqall = 0.0; if (update->ftol > 0.0) { - fmdotfm = max_torque(); + if (normstyle == 1) { // max torque norm + fmsq = max_torque(); + fmsqall = fmsq; + if (update->multireplica == 0) + MPI_Allreduce(&fmsq,&fmsqall,1,MPI_INT,MPI_MAX,universe->uworld); + } else { // Euclidean torque norm + fmsq = total_torque(); + fmsqall = fmsq; + if (update->multireplica == 0) + MPI_Allreduce(&fmsq,&fmsqall,1,MPI_INT,MPI_SUM,universe->uworld); + } + fmdotfm = fmsqall*fmsqall; if (update->multireplica == 0) { if (fmdotfm < update->ftol*update->ftol) return FTOL; } else { @@ -297,77 +309,3 @@ void MinSpin::advance_spins(double dts) // because no need for simplecticity } } - -/* ---------------------------------------------------------------------- - compute and return ||mag. torque||_2^2 -------------------------------------------------------------------------- */ - -double MinSpin::fmnorm_sqr() -{ - int nlocal = atom->nlocal; - double tx,ty,tz; - double **sp = atom->sp; - double **fm = atom->fm; - - // calc. magnetic torques - - double local_norm2_sqr = 0.0; - for (int i = 0; i < nlocal; i++) { - tx = (fm[i][1]*sp[i][2] - fm[i][2]*sp[i][1]); - ty = (fm[i][2]*sp[i][0] - fm[i][0]*sp[i][2]); - tz = (fm[i][0]*sp[i][1] - fm[i][1]*sp[i][0]); - - local_norm2_sqr += tx*tx + ty*ty + tz*tz; - } - - // no extra atom calc. for spins - - if (nextra_atom) - error->all(FLERR,"extra atom option not available yet"); - - double norm2_sqr = 0.0; - MPI_Allreduce(&local_norm2_sqr,&norm2_sqr,1,MPI_DOUBLE,MPI_SUM,world); - - return norm2_sqr; -} - -/* ---------------------------------------------------------------------- - compute and return max_i||mag. torque_i||_2 -------------------------------------------------------------------------- */ - -double MinSpin::max_torque() -{ - double fmsq,fmaxsqone,fmaxsqloc,fmaxsqall; - int nlocal = atom->nlocal; - double hbar = force->hplanck/MY_2PI; - double tx,ty,tz; - double **sp = atom->sp; - double **fm = atom->fm; - - fmsq = fmaxsqone = fmaxsqloc = fmaxsqall = 0.0; - for (int i = 0; i < nlocal; i++) { - tx = fm[i][1] * sp[i][2] - fm[i][2] * sp[i][1]; - ty = fm[i][2] * sp[i][0] - fm[i][0] * sp[i][2]; - tz = fm[i][0] * sp[i][1] - fm[i][1] * sp[i][0]; - fmsq = tx * tx + ty * ty + tz * tz; - fmaxsqone = MAX(fmaxsqone,fmsq); - } - - // finding max fm on this replica - - fmaxsqloc = fmaxsqone; - MPI_Allreduce(&fmaxsqone,&fmaxsqloc,1,MPI_DOUBLE,MPI_MAX,world); - - // finding max fm over all replicas, if necessary - // this communicator would be invalid for multiprocess replicas - - fmaxsqall = fmaxsqloc; - if (update->multireplica == 1) { - fmaxsqall = fmaxsqloc; - MPI_Allreduce(&fmaxsqloc,&fmaxsqall,1,MPI_DOUBLE,MPI_MAX,universe->uworld); - } - - // multiply it by hbar so that units are in eV - - return sqrt(fmaxsqall) * hbar; -} diff --git a/src/SPIN/min_spin.h b/src/SPIN/min_spin.h index d6d49203d5..f2df81e58c 100644 --- a/src/SPIN/min_spin.h +++ b/src/SPIN/min_spin.h @@ -35,8 +35,6 @@ class MinSpin : public Min { int iterate(int); double evaluate_dt(); void advance_spins(double); - double fmnorm_sqr(); - double max_torque(); private: diff --git a/src/SPIN/min_spin_oso_cg.cpp b/src/SPIN/min_spin_oso_cg.cpp index 1c91fa1500..16a95c5c02 100644 --- a/src/SPIN/min_spin_oso_cg.cpp +++ b/src/SPIN/min_spin_oso_cg.cpp @@ -29,6 +29,7 @@ #include "universe.h" #include "atom.h" #include "citeme.h" +#include "comm.h" #include "force.h" #include "update.h" #include "output.h" @@ -99,6 +100,13 @@ void MinSpinOSO_CG::init() Min::init(); + // warning if line_search combined to gneb + + if ((nreplica >= 1) && (linestyle != 4) && (comm->me == 0)) + error->warning(FLERR,"Line search incompatible gneb"); + + // set back use_line_search to 0 if more than one replica + if (linestyle == 3 && nreplica == 1){ use_line_search = 1; } @@ -175,7 +183,7 @@ int MinSpinOSO_CG::iterate(int maxiter) { int nlocal = atom->nlocal; bigint ntimestep; - double fmdotfm; + double fmdotfm,fmsq,fmsqall; int flag, flagall; double **sp = atom->sp; double der_e_cur_tmp = 0.0; @@ -261,8 +269,20 @@ int MinSpinOSO_CG::iterate(int maxiter) // magnetic torque tolerance criterion // sync across replicas if running multi-replica minimization + fmdotfm = fmsq = fmsqall = 0.0; if (update->ftol > 0.0) { - fmdotfm = max_torque(); + if (normstyle == 1) { // max torque norm + fmsq = max_torque(); + fmsqall = fmsq; + if (update->multireplica == 0) + MPI_Allreduce(&fmsq,&fmsqall,1,MPI_INT,MPI_MAX,universe->uworld); + } else { // Euclidean torque norm + fmsq = total_torque(); + fmsqall = fmsq; + if (update->multireplica == 0) + MPI_Allreduce(&fmsq,&fmsqall,1,MPI_INT,MPI_SUM,universe->uworld); + } + fmdotfm = fmsqall*fmsqall; if (update->multireplica == 0) { if (fmdotfm < update->ftol*update->ftol) return FTOL; } else { @@ -353,6 +373,7 @@ void MinSpinOSO_CG::calc_search_direction() MPI_Allreduce(&g2old,&g2old_global,1,MPI_DOUBLE,MPI_SUM,world); // Sum over all replicas. Good for GNEB. + if (nreplica > 1) { g2 = g2_global * factor; g2old = g2old_global * factor; @@ -361,7 +382,9 @@ void MinSpinOSO_CG::calc_search_direction() } if (fabs(g2_global) < 1.0e-60) beta = 0.0; else beta = g2_global / g2old_global; + // calculate conjugate direction + for (int i = 0; i < 3 * nlocal; i++) { p_s[i] = (beta * p_s[i] - g_cur[i]) * factor; g_old[i] = g_cur[i] * factor; @@ -379,7 +402,7 @@ void MinSpinOSO_CG::advance_spins() { int nlocal = atom->nlocal; double **sp = atom->sp; - double rot_mat[9]; // exponential of matrix made of search direction + double rot_mat[9]; // exponential of matrix made of search direction double s_new[3]; // loop on all spins on proc. @@ -394,47 +417,6 @@ void MinSpinOSO_CG::advance_spins() } } -/* ---------------------------------------------------------------------- - compute and return max_i||mag. torque_i||_2 -------------------------------------------------------------------------- */ - -double MinSpinOSO_CG::max_torque() -{ - double fmsq,fmaxsqone,fmaxsqloc,fmaxsqall; - int nlocal = atom->nlocal; - double factor; - double hbar = force->hplanck/MY_2PI; - - if (use_line_search) factor = 1.0; - else factor = hbar; - - // finding max fm on this proc. - - fmsq = fmaxsqone = fmaxsqloc = fmaxsqall = 0.0; - for (int i = 0; i < nlocal; i++) { - fmsq = 0.0; - for (int j = 0; j < 3; j++) - fmsq += g_cur[3 * i + j] * g_cur[3 * i + j]; - fmaxsqone = MAX(fmaxsqone,fmsq); - } - - // finding max fm on this replica - - fmaxsqloc = fmaxsqone; - MPI_Allreduce(&fmaxsqone,&fmaxsqloc,1,MPI_DOUBLE,MPI_MAX,world); - - // finding max fm over all replicas, if necessary - // this communicator would be invalid for multiprocess replicas - - fmaxsqall = fmaxsqloc; - if (update->multireplica == 1) { - fmaxsqall = fmaxsqloc; - MPI_Allreduce(&fmaxsqloc,&fmaxsqall,1,MPI_DOUBLE,MPI_MAX,universe->uworld); - } - - return sqrt(fmaxsqall) * factor; -} - /* ---------------------------------------------------------------------- calculate 3x3 matrix exponential using Rodrigues' formula (R. Murray, Z. Li, and S. Shankar Sastry, @@ -456,15 +438,14 @@ void MinSpinOSO_CG::rodrigues_rotation(const double *upp_tr, double *out) fabs(upp_tr[1]) < 1.0e-40 && fabs(upp_tr[2]) < 1.0e-40){ - // if upp_tr is zero, return unity matrix - for(int k = 0; k < 3; k++){ - for(int m = 0; m < 3; m++){ - if (m == k) - out[3 * k + m] = 1.0; - else - out[3 * k + m] = 0.0; + // if upp_tr is zero, return unity matrix + + for(int k = 0; k < 3; k++){ + for(int m = 0; m < 3; m++){ + if (m == k) out[3 * k + m] = 1.0; + else out[3 * k + m] = 0.0; + } } - } return; } @@ -512,13 +493,14 @@ void MinSpinOSO_CG::rodrigues_rotation(const double *upp_tr, double *out) void MinSpinOSO_CG::vm3(const double *m, const double *v, double *out) { for(int i = 0; i < 3; i++){ - //out[i] *= 0.0; out[i] = 0.0; - for(int j = 0; j < 3; j++) - out[i] += *(m + 3 * j + i) * v[j]; + for(int j = 0; j < 3; j++) out[i] += *(m + 3 * j + i) * v[j]; } } +/* ---------------------------------------------------------------------- + advance spins +------------------------------------------------------------------------- */ void MinSpinOSO_CG::make_step(double c, double *energy_and_der) { @@ -586,7 +568,7 @@ int MinSpinOSO_CG::calc_and_make_step(double a, double b, int index) } return 1; } - else{ + else { double r,f0,f1,df0,df1; r = b - a; f0 = eprevious; diff --git a/src/SPIN/min_spin_oso_cg.h b/src/SPIN/min_spin_oso_cg.h index 41253f440f..30d9adf066 100644 --- a/src/SPIN/min_spin_oso_cg.h +++ b/src/SPIN/min_spin_oso_cg.h @@ -25,44 +25,44 @@ MinimizeStyle(spin_oso_cg, MinSpinOSO_CG) namespace LAMMPS_NS { class MinSpinOSO_CG: public Min { - public: - MinSpinOSO_CG(class LAMMPS *); - virtual ~MinSpinOSO_CG(); - void init(); - void setup_style(); - int modify_param(int, char **); - void reset_vectors(); - int iterate(int); - private: - double dt; // global timestep - double dts; // spin timestep - int ireplica,nreplica; // for neb - double *spvec; // variables for atomic dof, as 1d vector - double *fmvec; // variables for atomic dof, as 1d vector - double *g_old; // gradient vector at previous step - double *g_cur; // current gradient vector - double *p_s; // search direction vector - double **sp_copy; // copy of the spins - int local_iter; // for neb - int nlocal_max; // max value of nlocal (for size of lists) - double discrete_factor; // factor for spin timestep evaluation + public: + MinSpinOSO_CG(class LAMMPS *); + virtual ~MinSpinOSO_CG(); + void init(); + void setup_style(); + void reset_vectors(); + int modify_param(int, char **); + int iterate(int); - double evaluate_dt(); - void advance_spins(); - void calc_gradient(); - void calc_search_direction(); - double maximum_rotation(double *); - void vm3(const double *, const double *, double *); - void rodrigues_rotation(const double *, double *); - int calc_and_make_step(double, double, int); - int awc(double, double, double, double); - void make_step(double, double *); - double max_torque(); - double der_e_cur; // current derivative along search dir. - double der_e_pr; // previous derivative along search dir. - int use_line_search; // use line search or not. + private: + int local_iter; // for neb + int nlocal_max; // max value of nlocal (for size of lists) + int use_line_search; // use line search or not. + int ireplica,nreplica; // for neb + double dt; // global timestep + double dts; // spin timestep + double discrete_factor; // factor for spin timestep evaluation + double der_e_cur; // current derivative along search dir. + double der_e_pr; // previous derivative along search dir. + double *spvec; // variables for atomic dof, as 1d vector + double *fmvec; // variables for atomic dof, as 1d vector + double *g_old; // gradient vector at previous step + double *g_cur; // current gradient vector + double *p_s; // search direction vector + double **sp_copy; // copy of the spins - bigint last_negative; + void advance_spins(); + void calc_gradient(); + void calc_search_direction(); + void vm3(const double *, const double *, double *); + void rodrigues_rotation(const double *, double *); + void make_step(double, double *); + int calc_and_make_step(double, double, int); + int awc(double, double, double, double); + double evaluate_dt(); + double maximum_rotation(double *); + + bigint last_negative; }; } diff --git a/src/SPIN/min_spin_oso_lbfgs.cpp b/src/SPIN/min_spin_oso_lbfgs.cpp index b9315d706e..2913ef4101 100644 --- a/src/SPIN/min_spin_oso_lbfgs.cpp +++ b/src/SPIN/min_spin_oso_lbfgs.cpp @@ -26,9 +26,9 @@ #include #include #include "min_spin_oso_lbfgs.h" -#include "universe.h" #include "atom.h" #include "citeme.h" +#include "comm.h" #include "force.h" #include "update.h" #include "output.h" @@ -107,6 +107,13 @@ void MinSpinOSO_LBFGS::init() Min::init(); + // warning if line_search combined to gneb + + if ((nreplica >= 1) && (linestyle != 4) && (comm->me == 0)) + error->warning(FLERR,"Line search incompatible gneb"); + + // set back use_line_search to 0 if more than one replica + if (linestyle != 4 && nreplica == 1){ use_line_search = 1; } @@ -188,7 +195,7 @@ int MinSpinOSO_LBFGS::iterate(int maxiter) { int nlocal = atom->nlocal; bigint ntimestep; - double fmdotfm; + double fmdotfm,fmsq,fmsqall; int flag, flagall; double **sp = atom->sp; double der_e_cur_tmp = 0.0; @@ -280,8 +287,20 @@ int MinSpinOSO_LBFGS::iterate(int maxiter) // magnetic torque tolerance criterion // sync across replicas if running multi-replica minimization + fmdotfm = fmsq = fmsqall = 0.0; if (update->ftol > 0.0) { - fmdotfm = max_torque(); + if (normstyle == 1) { // max torque norm + fmsq = max_torque(); + fmsqall = fmsq; + if (update->multireplica == 0) + MPI_Allreduce(&fmsq,&fmsqall,1,MPI_INT,MPI_MAX,universe->uworld); + } else { // Euclidean torque norm + fmsq = total_torque(); + fmsqall = fmsq; + if (update->multireplica == 0) + MPI_Allreduce(&fmsq,&fmsqall,1,MPI_INT,MPI_SUM,universe->uworld); + } + fmdotfm = fmsqall*fmsqall; if (update->multireplica == 0) { if (fmdotfm < update->ftol*update->ftol) return FTOL; } else { @@ -534,42 +553,6 @@ void MinSpinOSO_LBFGS::advance_spins() } } -/* ---------------------------------------------------------------------- - compute and return max_i||mag. torque_i||_2 -------------------------------------------------------------------------- */ - -double MinSpinOSO_LBFGS::max_torque() -{ - double fmsq,fmaxsqone,fmaxsqloc,fmaxsqall; - int nlocal = atom->nlocal; - - // finding max fm on this proc. - - fmsq = fmaxsqone = fmaxsqloc = fmaxsqall = 0.0; - for (int i = 0; i < nlocal; i++) { - fmsq = 0.0; - for (int j = 0; j < 3; j++) - fmsq += g_cur[3 * i + j] * g_cur[3 * i + j]; - fmaxsqone = MAX(fmaxsqone,fmsq); - } - - // finding max fm on this replica - - fmaxsqloc = fmaxsqone; - MPI_Allreduce(&fmaxsqone,&fmaxsqloc,1,MPI_DOUBLE,MPI_MAX,world); - - // finding max fm over all replicas, if necessary - // this communicator would be invalid for multiprocess replicas - - fmaxsqall = fmaxsqloc; - if (update->multireplica == 1) { - fmaxsqall = fmaxsqloc; - MPI_Allreduce(&fmaxsqloc,&fmaxsqall,1,MPI_DOUBLE,MPI_MAX,universe->uworld); - } - - return sqrt(fmaxsqall); -} - /* ---------------------------------------------------------------------- calculate 3x3 matrix exponential using Rodrigues' formula (R. Murray, Z. Li, and S. Shankar Sastry, diff --git a/src/SPIN/min_spin_oso_lbfgs.h b/src/SPIN/min_spin_oso_lbfgs.h index 204f6bf058..9bd36afa8b 100644 --- a/src/SPIN/min_spin_oso_lbfgs.h +++ b/src/SPIN/min_spin_oso_lbfgs.h @@ -25,45 +25,45 @@ MinimizeStyle(spin_oso_lbfgs, MinSpinOSO_LBFGS) namespace LAMMPS_NS { class MinSpinOSO_LBFGS: public Min { - public: - MinSpinOSO_LBFGS(class LAMMPS *); - virtual ~MinSpinOSO_LBFGS(); - void init(); - void setup_style(); - int modify_param(int, char **); - void reset_vectors(); - int iterate(int); - private: - int ireplica,nreplica; // for neb - double *spvec; // variables for atomic dof, as 1d vector - double *fmvec; // variables for atomic dof, as 1d vector - double *g_old; // gradient vector at previous step - double *g_cur; // current gradient vector - double *p_s; // search direction vector - int local_iter; // for neb - int nlocal_max; // max value of nlocal (for size of lists) + public: + MinSpinOSO_LBFGS(class LAMMPS *); + virtual ~MinSpinOSO_LBFGS(); + void init(); + void setup_style(); + int modify_param(int, char **); + void reset_vectors(); + int iterate(int); - void advance_spins(); - void calc_gradient(); - void calc_search_direction(); - double maximum_rotation(double *); - void vm3(const double *, const double *, double *); - void rodrigues_rotation(const double *, double *); - int calc_and_make_step(double, double, int); - int awc(double, double, double, double); - void make_step(double, double *); - double max_torque(); - double der_e_cur; // current derivative along search dir. - double der_e_pr; // previous derivative along search dir. - int use_line_search; // use line search or not. - double maxepsrot; + private: + int local_iter; // for neb + int use_line_search; // use line search or not. + int nlocal_max; // max value of nlocal (for size of lists) + int ireplica,nreplica; // for neb + double der_e_cur; // current derivative along search dir. + double der_e_pr; // previous derivative along search dir. + double maxepsrot; + double *spvec; // variables for atomic dof, as 1d vector + double *fmvec; // variables for atomic dof, as 1d vector + double *g_old; // gradient vector at previous step + double *g_cur; // current gradient vector + double *p_s; // search direction vector - double *rho; // estimation of curvature - double **ds; // change in rotation matrix between two iterations, da - double **dy; // change in gradients between two iterations, dg - double **sp_copy; // copy of the spins - int num_mem; // number of stored steps - bigint last_negative; + void advance_spins(); + void calc_gradient(); + void calc_search_direction(); + void vm3(const double *, const double *, double *); + void rodrigues_rotation(const double *, double *); + void make_step(double, double *); + int calc_and_make_step(double, double, int); + int awc(double, double, double, double); + double maximum_rotation(double *); + + double *rho; // estimation of curvature + double **ds; // change in rotation matrix between two iterations, da + double **dy; // change in gradients between two iterations, dg + double **sp_copy; // copy of the spins + int num_mem; // number of stored steps + bigint last_negative; }; } diff --git a/src/min.cpp b/src/min.cpp index 2a42a444a0..e476b1abc8 100644 --- a/src/min.cpp +++ b/src/min.cpp @@ -42,10 +42,12 @@ #include "output.h" #include "thermo.h" #include "timer.h" +#include "math_const.h" #include "memory.h" #include "error.h" using namespace LAMMPS_NS; +using namespace MathConst; /* ---------------------------------------------------------------------- */ @@ -54,6 +56,7 @@ Min::Min(LAMMPS *lmp) : Pointers(lmp) dmax = 0.1; searchflag = 0; linestyle = 1; + normstyle = 0; elist_global = elist_atom = NULL; vlist_global = vlist_atom = NULL; @@ -653,6 +656,14 @@ void Min::modify_params(int narg, char **arg) if (strcmp(arg[iarg+1],"backtrack") == 0) linestyle = 0; else if (strcmp(arg[iarg+1],"quadratic") == 0) linestyle = 1; else if (strcmp(arg[iarg+1],"forcezero") == 0) linestyle = 2; + else if (strcmp(arg[iarg+1],"spin_cubic") == 0) linestyle = 3; + else if (strcmp(arg[iarg+1],"spin_none") == 0) linestyle = 4; + else error->all(FLERR,"Illegal min_modify command"); + iarg += 2; + } else if (strcmp(arg[iarg],"norm") == 0) { + if (iarg+2 > narg) error->all(FLERR,"Illegal min_modify command"); + if (strcmp(arg[iarg+1],"euclidean") == 0) normstyle = 0; + else if (strcmp(arg[iarg+1],"max") == 0) normstyle = 1; else error->all(FLERR,"Illegal min_modify command"); iarg += 2; } else { @@ -816,6 +827,69 @@ double Min::fnorm_inf() return norm_inf; } +/* ---------------------------------------------------------------------- + compute and return sum_i||mag. torque_i||_2 (in eV) +------------------------------------------------------------------------- */ + +double Min::total_torque() +{ + double fmsq,ftotsqone,ftotsqall; + int nlocal = atom->nlocal; + double hbar = force->hplanck/MY_2PI; + double tx,ty,tz; + double **sp = atom->sp; + double **fm = atom->fm; + + fmsq = ftotsqone = ftotsqall = 0.0; + for (int i = 0; i < nlocal; i++) { + tx = fm[i][1] * sp[i][2] - fm[i][2] * sp[i][1]; + ty = fm[i][2] * sp[i][0] - fm[i][0] * sp[i][2]; + tz = fm[i][0] * sp[i][1] - fm[i][1] * sp[i][0]; + fmsq = tx * tx + ty * ty + tz * tz; + ftotsqone += fmsq; + } + + // summing all fmsqtot on this replica + + MPI_Allreduce(&ftotsqone,&ftotsqall,1,MPI_DOUBLE,MPI_SUM,world); + + // multiply it by hbar so that units are in eV + + return sqrt(ftotsqall) * hbar; +} + +/* ---------------------------------------------------------------------- + compute and return max_i ||mag. torque_i|| (in eV) +------------------------------------------------------------------------- */ + +double Min::max_torque() +{ + double fmsq,fmaxsqone,fmaxsqall; + int nlocal = atom->nlocal; + double hbar = force->hplanck/MY_2PI; + double tx,ty,tz; + double **sp = atom->sp; + double **fm = atom->fm; + + fmsq = fmaxsqone = fmaxsqall = 0.0; + for (int i = 0; i < nlocal; i++) { + tx = fm[i][1] * sp[i][2] - fm[i][2] * sp[i][1]; + ty = fm[i][2] * sp[i][0] - fm[i][0] * sp[i][2]; + tz = fm[i][0] * sp[i][1] - fm[i][1] * sp[i][0]; + fmsq = tx * tx + ty * ty + tz * tz; + fmaxsqone = MAX(fmaxsqone,fmsq); + } + + // finding max fm on this replica + + fmaxsqall = fmaxsqone; + MPI_Allreduce(&fmaxsqone,&fmaxsqall,1,MPI_DOUBLE,MPI_MAX,world); + + // multiply it by hbar so that units are in eV + + return sqrt(fmaxsqall) * hbar; +} + /* ---------------------------------------------------------------------- possible stop conditions ------------------------------------------------------------------------- */ diff --git a/src/min.h b/src/min.h index a63254231c..e18d0dd677 100644 --- a/src/min.h +++ b/src/min.h @@ -42,6 +42,10 @@ class Min : protected Pointers { double fnorm_sqr(); double fnorm_inf(); + // methods for spin minimizers + double max_torque(); + double total_torque(); + virtual void init_style() {} virtual void setup_style() = 0; virtual void reset_vectors() = 0; @@ -56,8 +60,11 @@ class Min : protected Pointers { int virial_style; // compute virial explicitly or implicitly int external_force_clear; // clear forces locally or externally - double dmax; // max dist to move any atom in one step - int linestyle; // 0 = backtrack, 1 = quadratic, 2 = forcezero + double dmax; // max dist to move any atom in one step + int linestyle; // 0 = backtrack, 1 = quadratic, 2 = forcezero + // 3 = spin_cubic, 4 = spin_none + + int normstyle; // 0 = Euclidean norm, 1 = inf. norm int nelist_global,nelist_atom; // # of PE,virial computes to check int nvlist_global,nvlist_atom; @@ -102,9 +109,6 @@ class Min : protected Pointers { double energy_force(int); void force_clear(); - double compute_force_norm_sqr(); - double compute_force_norm_inf(); - void ev_setup(); void ev_set(bigint); diff --git a/src/min_cg.cpp b/src/min_cg.cpp index 20e8cc30dd..9801e57f4d 100644 --- a/src/min_cg.cpp +++ b/src/min_cg.cpp @@ -37,7 +37,7 @@ MinCG::MinCG(LAMMPS *lmp) : MinLineSearch(lmp) {} int MinCG::iterate(int maxiter) { int i,m,n,fail,ntimestep; - double beta,gg,dot[2],dotall[2]; + double beta,gg,dot[2],dotall[2],fmax,fmaxall; double *fatom,*gatom,*hatom; // nlimit = max # of CG iterations before restarting @@ -87,10 +87,12 @@ int MinCG::iterate(int maxiter) // force tolerance criterion + fmax = fmaxall = 0.0; dot[0] = dot[1] = 0.0; for (i = 0; i < nvec; i++) { dot[0] += fvec[i]*fvec[i]; dot[1] += fvec[i]*g[i]; + fmax = MAX(fmax,fvec[i]*fvec[i]); } if (nextra_atom) for (m = 0; m < nextra_atom; m++) { @@ -100,16 +102,22 @@ int MinCG::iterate(int maxiter) for (i = 0; i < n; i++) { dot[0] += fatom[i]*fatom[i]; dot[1] += fatom[i]*gatom[i]; + fmax = MAX(fmax,fatom[i]*fatom[i]); } } MPI_Allreduce(dot,dotall,2,MPI_DOUBLE,MPI_SUM,world); + MPI_Allreduce(&fmax,&fmaxall,2,MPI_DOUBLE,MPI_MAX,world); 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 (normstyle == 1) { // max force norm + if (fmax < update->ftol*update->ftol) return FTOL; + } else { // Euclidean force norm + 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 diff --git a/src/min_fire.cpp b/src/min_fire.cpp index a50071d562..a0a3bce8ba 100644 --- a/src/min_fire.cpp +++ b/src/min_fire.cpp @@ -80,7 +80,7 @@ void MinFire::reset_vectors() int MinFire::iterate(int maxiter) { bigint ntimestep; - double vmax,vdotf,vdotfall,vdotv,vdotvall,fdotf,fdotfall; + double vmax,vdotf,vdotfall,vdotv,vdotvall,fdotf,fdotfloc,fdotfall; double scale1,scale2; double dtvone,dtv,dtf,dtfm; int flag,flagall; @@ -250,7 +250,15 @@ int MinFire::iterate(int maxiter) // sync across replicas if running multi-replica minimization if (update->ftol > 0.0) { - fdotf = fnorm_sqr(); + if (normstyle == 1) { // max force norm + fdotf = fnorm_inf(); + fdotfloc = fdotf; + MPI_Allreduce(&fdotfloc,&fdotf,1,MPI_INT,MPI_MAX,universe->uworld); + } else { // Euclidean force norm + fdotf = fnorm_sqr(); + fdotfloc = fdotf; + MPI_Allreduce(&fdotfloc,&fdotf,1,MPI_INT,MPI_SUM,universe->uworld); + } if (update->multireplica == 0) { if (fdotf < update->ftol*update->ftol) return FTOL; } else { diff --git a/src/min_hftn.cpp b/src/min_hftn.cpp index 0c834fbeb4..9f8695f151 100644 --- a/src/min_hftn.cpp +++ b/src/min_hftn.cpp @@ -20,6 +20,7 @@ #include #include #include "atom.h" +#include "error.h" #include "fix_minimize.h" #include "min_hftn.h" #include "modify.h" @@ -111,6 +112,9 @@ void MinHFTN::init() { Min::init(); + if (normstyle == 1) + error->all(FLERR,"Incorrect min_modify option"); + for (int i = 1; i < NUM_HFTN_ATOM_BASED_VECTORS; i++) { if (_daExtraGlobal[i] != NULL) delete [] _daExtraGlobal[i]; diff --git a/src/min_quickmin.cpp b/src/min_quickmin.cpp index 8b48816355..d6507cfcde 100644 --- a/src/min_quickmin.cpp +++ b/src/min_quickmin.cpp @@ -76,7 +76,7 @@ void MinQuickMin::reset_vectors() int MinQuickMin::iterate(int maxiter) { bigint ntimestep; - double vmax,vdotf,vdotfall,fdotf,fdotfall,scale; + double vmax,vdotf,vdotfall,fdotf,fdotfloc,fdotfall,scale; double dtvone,dtv,dtf,dtfm; int flag,flagall; @@ -216,7 +216,15 @@ int MinQuickMin::iterate(int maxiter) // sync across replicas if running multi-replica minimization if (update->ftol > 0.0) { - fdotf = fnorm_sqr(); + if (normstyle == 1) { // max force norm + fdotf = fnorm_inf(); + fdotfloc = fdotf; + MPI_Allreduce(&fdotfloc,&fdotf,1,MPI_INT,MPI_MAX,universe->uworld); + } else { // Euclidean force norm + fdotf = fnorm_sqr(); + fdotfloc = fdotf; + MPI_Allreduce(&fdotfloc,&fdotf,1,MPI_INT,MPI_SUM,universe->uworld); + } if (update->multireplica == 0) { if (fdotf < update->ftol*update->ftol) return FTOL; } else { diff --git a/src/min_sd.cpp b/src/min_sd.cpp index 5d44437ca0..60386df82c 100644 --- a/src/min_sd.cpp +++ b/src/min_sd.cpp @@ -79,7 +79,8 @@ int MinSD::iterate(int maxiter) // force tolerance criterion - fdotf = fnorm_sqr(); + if (normstyle == 1) fdotf = fnorm_inf(); // max force norm + else fdotf = fnorm_sqr(); // Euclidean force norm if (fdotf < update->ftol*update->ftol) return FTOL; // set new search direction h to f = -Grad(x)