diff --git a/doc/src/min_style.txt b/doc/src/min_style.txt index b1a9da997d..f8c05d5483 100644 --- a/doc/src/min_style.txt +++ b/doc/src/min_style.txt @@ -11,7 +11,7 @@ min_style command :h3 min_style style :pre -style = {cg} or {hftn} or {sd} or {quickmin} or {fire} or {spinmin} :ul +style = {cg} or {hftn} or {sd} or {quickmin} or {fire} or {spin} :ul [Examples:] @@ -62,18 +62,14 @@ the velocity non-parallel to the current force vector. The velocity of each atom is initialized to 0.0 by this style, at the beginning of a minimization. -Style {spinmin} is a damped spin dynamics with a variable +Style {spin} is a damped spin dynamics with a variable timestep as described in "(Tranchida)"_#Tranchida. -The value of the fictitious Gilbert damping and of the dividing -factor for the adaptive timestep can be modified by the -{alpha_damp} and {discret_factor} options respectively. -Those options can be defined using the "min_modify"_min_modify.html -command. +See the "min/spin"_min_spin.html doc page for more information. Either the {quickmin} and {fire} styles are useful in the context of nudged elastic band (NEB) calculations via the "neb"_neb.html command. -The {spinmin} style is useful in the context of geodesic nudged +The {spin} style is useful in the context of geodesic nudged elastic band (GNEB) calculations via the "neb/spin"_neb_spin.html command. diff --git a/examples/SPIN/gneb/in.gneb.iron b/examples/SPIN/gneb/in.gneb.iron index 7aab0c04c3..a8028392a1 100644 --- a/examples/SPIN/gneb/in.gneb.iron +++ b/examples/SPIN/gneb/in.gneb.iron @@ -55,7 +55,7 @@ variable u universe 1 2 3 4 #dump 1 all custom 100 dump.lammpstrj type x y z c_outsp[1] c_outsp[2] c_outsp[3] c_outsp[4] c_outsp[5] c_outsp[6] c_outsp[7] dump 1 all custom 200 dump.$u type x y z c_outsp[1] c_outsp[2] c_outsp[3] -min_style spinmin +min_style spin min_modify alpha_damp 1.0 discret_factor 10.0 neb/spin 1.0e-12 1.0e-12 50000 50000 10 final final.iron_spin #neb/spin 1.0e-6 1.0e-6 1000 10 10 final final.iron_spin diff --git a/examples/SPIN/spinmin/in.spinmin.bfo b/examples/SPIN/spinmin/in.spinmin.bfo index a00af8833c..5b678c8a4d 100644 --- a/examples/SPIN/spinmin/in.spinmin.bfo +++ b/examples/SPIN/spinmin/in.spinmin.bfo @@ -52,6 +52,6 @@ thermo_modify format float %20.15g compute outsp all property/atom spx spy spz sp fmx fmy fmz dump 1 all custom 50 dump.lammpstrj type x y z c_outsp[1] c_outsp[2] c_outsp[3] c_outsp[4] c_outsp[5] c_outsp[6] c_outsp[7] -min_style spinmin +min_style spin min_modify alpha_damp 1.0 discret_factor 10.0 minimize 0.0 0.0 10000 1000 diff --git a/examples/SPIN/spinmin/in.spinmin.iron b/examples/SPIN/spinmin/in.spinmin.iron index 5a15082122..b87a811cc7 100644 --- a/examples/SPIN/spinmin/in.spinmin.iron +++ b/examples/SPIN/spinmin/in.spinmin.iron @@ -52,6 +52,6 @@ thermo_modify format float %20.15g compute outsp all property/atom spx spy spz sp fmx fmy fmz dump 1 all custom 100 dump.lammpstrj type x y z c_outsp[1] c_outsp[2] c_outsp[3] c_outsp[4] c_outsp[5] c_outsp[6] c_outsp[7] -min_style spinmin +min_style spin min_modify alpha_damp 1.0 discret_factor 10.0 minimize 1.0e-10 1.0e-10 100000 1000 diff --git a/src/SPIN/README b/src/SPIN/README index e371e39767..c3c02b1445 100644 --- a/src/SPIN/README +++ b/src/SPIN/README @@ -8,13 +8,16 @@ atom in the system * integrating the equations of motion for the coupled spin-lattice system * implementing magnetic pair interactions and magnetic forces * thermostating and applying a transverse damping to the magnetic spins +* minimizing spin configurations with an adaptive timestep scheme +* performing geodesic NEB calculations * computing and outputing magnetic quantities The different options provided by this package are explained in the LAMMPS documentation. Once you have successfully built LAMMPS with this package, you can test -it using one of the input files provided from the examples/SPIN dir: +it using one of the input files provided from the examples/SPIN dir. +For example: ./lmp_serial < lammps/examples/SPIN/cobalt_hcp/in.spin.cobalt_hcp diff --git a/src/REPLICA/fix_neb_spin.cpp b/src/SPIN/fix_neb_spin.cpp similarity index 100% rename from src/REPLICA/fix_neb_spin.cpp rename to src/SPIN/fix_neb_spin.cpp diff --git a/src/REPLICA/fix_neb_spin.h b/src/SPIN/fix_neb_spin.h similarity index 100% rename from src/REPLICA/fix_neb_spin.h rename to src/SPIN/fix_neb_spin.h diff --git a/src/SPIN/fix_precession_spin.cpp b/src/SPIN/fix_precession_spin.cpp index 433a260e83..6ccb692033 100644 --- a/src/SPIN/fix_precession_spin.cpp +++ b/src/SPIN/fix_precession_spin.cpp @@ -230,8 +230,8 @@ void FixPrecessionSpin::compute_single_precession(int i, double spi[3], double f void FixPrecessionSpin::compute_zeeman(int i, double fmi[3]) { double **sp = atom->sp; - fmi[0] += sp[i][0]*hx; - fmi[1] += sp[i][1]*hy; + fmi[0] += sp[i][3]*hx; + fmi[1] += sp[i][3]*hy; fmi[2] += sp[i][3]*hz; } diff --git a/src/SPIN/min_spinmin.cpp b/src/SPIN/min_spin.cpp similarity index 90% rename from src/SPIN/min_spinmin.cpp rename to src/SPIN/min_spin.cpp index 808a5359bc..ac8f22186e 100644 --- a/src/SPIN/min_spinmin.cpp +++ b/src/SPIN/min_spin.cpp @@ -19,7 +19,7 @@ #include #include -#include "min_spinmin.h" +#include "min_spin.h" #include "universe.h" #include "atom.h" #include "force.h" @@ -27,7 +27,6 @@ #include "output.h" #include "timer.h" #include "error.h" - #include #include #include "modify.h" @@ -46,11 +45,11 @@ using namespace MathConst; /* ---------------------------------------------------------------------- */ -MinSpinMin::MinSpinMin(LAMMPS *lmp) : Min(lmp) {} +MinSpin::MinSpin(LAMMPS *lmp) : Min(lmp) {} /* ---------------------------------------------------------------------- */ -void MinSpinMin::init() +void MinSpin::init() { alpha_damp = 1.0; discret_factor = 10.0; @@ -63,21 +62,43 @@ void MinSpinMin::init() /* ---------------------------------------------------------------------- */ -void MinSpinMin::setup_style() +void MinSpin::setup_style() { double **v = atom->v; int nlocal = atom->nlocal; + // check if the atom/spin style is defined + + if (!atom->sp_flag) + error->all(FLERR,"min/spin requires atom/spin style"); + for (int i = 0; i < nlocal; i++) v[i][0] = v[i][1] = v[i][2] = 0.0; } +/* ---------------------------------------------------------------------- */ + +int MinSpin::modify_param(int narg, char **arg) +{ + if (strcmp(arg[0],"alpha_damp") == 0) { + if (narg < 2) error->all(FLERR,"Illegal fix_modify command"); + alpha_damp = force->numeric(FLERR,arg[1]); + return 2; + } + if (strcmp(arg[0],"discret_factor") == 0) { + if (narg < 2) error->all(FLERR,"Illegal fix_modify command"); + discret_factor = force->numeric(FLERR,arg[1]); + return 2; + } + return 0; +} + /* ---------------------------------------------------------------------- set current vector lengths and pointers called after atoms have migrated ------------------------------------------------------------------------- */ -void MinSpinMin::reset_vectors() +void MinSpin::reset_vectors() { // atomic dof @@ -96,7 +117,7 @@ void MinSpinMin::reset_vectors() minimization via damped spin dynamics ------------------------------------------------------------------------- */ -int MinSpinMin::iterate(int maxiter) +int MinSpin::iterate(int maxiter) { bigint ntimestep; double fmdotfm,fmdotfmall; @@ -172,7 +193,7 @@ int MinSpinMin::iterate(int maxiter) evaluate max timestep ---------------------------------------------------------------------- */ -double MinSpinMin::evaluate_dt() +double MinSpin::evaluate_dt() { double dtmax; double fmsq; @@ -218,7 +239,7 @@ double MinSpinMin::evaluate_dt() geometric damped advance of spins ---------------------------------------------------------------------- */ -void MinSpinMin::advance_spins(double dts) +void MinSpin::advance_spins(double dts) { int nlocal = atom->nlocal; int *mask = atom->mask; @@ -282,7 +303,7 @@ void MinSpinMin::advance_spins(double dts) compute and return ||mag. torque||_2^2 ------------------------------------------------------------------------- */ -double MinSpinMin::fmnorm_sqr() +double MinSpin::fmnorm_sqr() { int i,n; double *fmatom; diff --git a/src/SPIN/min_spinmin.h b/src/SPIN/min_spin.h similarity index 78% rename from src/SPIN/min_spinmin.h rename to src/SPIN/min_spin.h index abc532a3d5..569bcbaab2 100644 --- a/src/SPIN/min_spinmin.h +++ b/src/SPIN/min_spin.h @@ -13,23 +13,24 @@ #ifdef MINIMIZE_CLASS -MinimizeStyle(spinmin,MinSpinMin) +MinimizeStyle(spin,MinSpin) #else -#ifndef LMP_MIN_SPINMIN_H -#define LMP_MIN_SPINMIN_H +#ifndef LMP_MIN_SPIN_H +#define LMP_MIN_SPIN_H #include "min.h" namespace LAMMPS_NS { -class MinSpinMin : public Min { +class MinSpin : public Min { public: - MinSpinMin(class LAMMPS *); - ~MinSpinMin() {} + MinSpin(class LAMMPS *); + ~MinSpin() {} void init(); void setup_style(); + int modify_param(int, char **); void reset_vectors(); int iterate(int); double evaluate_dt(); @@ -43,6 +44,9 @@ class MinSpinMin : public Min { double dt; double dts; + double alpha_damp; // damping for spin minimization + double discret_factor; // factor for spin timestep evaluation + double *spvec; // variables for atomic dof, as 1d vector double *fmvec; // variables for atomic dof, as 1d vector diff --git a/src/REPLICA/neb_spin.cpp b/src/SPIN/neb_spin.cpp similarity index 93% rename from src/REPLICA/neb_spin.cpp rename to src/SPIN/neb_spin.cpp index 38f5b530da..f5d9a75020 100644 --- a/src/REPLICA/neb_spin.cpp +++ b/src/SPIN/neb_spin.cpp @@ -109,7 +109,7 @@ NEB_spin::NEB_spin(LAMMPS *lmp, double etol_in, double ftol_in, int n1steps_in, sp[i][0] = spfinal[0]; sp[i][1] = spfinal[1]; sp[i][2] = spfinal[2]; - + ii += 3; } } @@ -228,8 +228,8 @@ void NEB_spin::run() if (update->minimize->searchflag) error->all(FLERR,"NEB_spin requires damped dynamics minimizer"); - if (strcmp(update->minimize_style,"spinmin") != 0) - error->all(FLERR,"NEB_spin requires spinmin minimizer"); + if (strcmp(update->minimize_style,"spin") != 0) + error->all(FLERR,"NEB_spin requires spin minimizer"); // setup regular NEB_spin minimization @@ -530,10 +530,14 @@ void NEB_spin::readfile(char *file, int flag) spfinal[0] = spx; spfinal[1] = spy; spfinal[2] = spz; + + printf("test spinit[0]:%g \n",sp[m][0]); // interpolate intermediate spin states initial_rotation(spinit,spfinal,fraction); + + printf("test spfinal[0]:%g \n",spfinal[0]); sp[m][0] = spfinal[0]; sp[m][1] = spfinal[1]; @@ -556,6 +560,8 @@ void NEB_spin::readfile(char *file, int flag) nread += nchunk; } + printf("test sp[1][2]:%g \n",sp[1][2]); + // check that all atom IDs in file were found by a proc if (flag == 0) { @@ -613,11 +619,15 @@ void NEB_spin::initial_rotation(double *spi, double *sploc, double fraction) spfy = sploc[1]; spfz = sploc[2]; + iphi = itheta = fphi = ftheta = 0.0; + iphi = acos(spiz); - itheta = acos(spix/sin(iphi)); + if (sin(iphi) != 0.0) + itheta = acos(spix/sin(iphi)); fphi = acos(spfz); - ftheta = acos(spfx/sin(fphi)); + if (sin(fphi) != 0.0) + ftheta = acos(spfx/sin(fphi)); kphi = iphi + fraction*(fphi-iphi); ktheta = itheta + fraction*(ftheta-itheta); @@ -626,12 +636,69 @@ void NEB_spin::initial_rotation(double *spi, double *sploc, double fraction) spky = sin(ktheta)*sin(kphi); spkz = cos(kphi); - iknorm = spkx*spkx+spky*spky+spkz*spkz; + double knormsq = spkx*spkx+spky*spky+spkz*spkz; + if (knormsq != 0.0) + iknorm = 1.0/sqrt(knormsq); spkx *= iknorm; spky *= iknorm; spkz *= iknorm; + //sploc[0] = spkx; + //sploc[1] = spky; + //sploc[2] = spkz; + + //double kx,ky,kz; + //double spix,spiy,spiz,spfx,spfy,spfz; + //double kcrossx,kcrossy,kcrossz,knormsq; + //double spkx,spky,spkz; + //double sdot,omega,iknorm; + + //spix = spi[0]; + //spiy = spi[1]; + //spiz = spi[2]; + + //spfx = sploc[0]; + //spfy = sploc[1]; + //spfz = sploc[2]; + // + //kx = spiy*spfz - spiz*spfy; + //ky = spiz*spfx - spix*spfz; + //kz = spix*spfy - spiy*spfx; + + //knormsq = kx*kx+ky*ky+kz*kz; + // + //if (knormsq != 0.0) { + // iknorm = 1.0/sqrt(knormsq); + // kx *= iknorm; + // ky *= iknorm; + // kz *= iknorm; + //} + // + //kcrossx = ky*spiz - kz*spiy; + //kcrossy = kz*spix - kx*spiz; + //kcrossz = kx*spiy - ky*spix; + + //sdot = spix*spfx + spiy*spfy + spiz*spfz; + + //omega = acos(sdot); + //omega *= fraction; + + //spkx = spix*cos(omega) + kcrossx*sin(omega); + //spky = spiy*cos(omega) + kcrossy*sin(omega); + //spkz = spiz*cos(omega) + kcrossz*sin(omega); + // + //iknorm = 1.0/sqrt(spkx*spkx+spky*spky+spkz*spkz); + //if (iknorm == 0.0) + // error->all(FLERR,"Incorrect rotation operation"); + + //spkx *= iknorm; + //spky *= iknorm; + //spkz *= iknorm; + + printf("init: %g %g %g \n",spix,spiy,spiz); + printf("fina: %g %g %g \n",spkx,spky,spkz); + sploc[0] = spkx; sploc[1] = spky; sploc[2] = spkz; diff --git a/src/REPLICA/neb_spin.h b/src/SPIN/neb_spin.h similarity index 100% rename from src/REPLICA/neb_spin.h rename to src/SPIN/neb_spin.h diff --git a/src/min.cpp b/src/min.cpp index c75db6e2b0..2a42a444a0 100644 --- a/src/min.cpp +++ b/src/min.cpp @@ -655,15 +655,11 @@ void Min::modify_params(int narg, char **arg) else if (strcmp(arg[iarg+1],"forcezero") == 0) linestyle = 2; else error->all(FLERR,"Illegal min_modify command"); iarg += 2; - } else if (strcmp(arg[iarg],"alpha_damp") == 0) { - if (iarg+2 > narg) error->all(FLERR,"Illegal min_modify command"); - alpha_damp = force->numeric(FLERR,arg[iarg+1]); - iarg += 2; - } else if (strcmp(arg[iarg],"discret_factor") == 0) { - if (iarg+2 > narg) error->all(FLERR,"Illegal min_modify command"); - discret_factor = force->numeric(FLERR,arg[iarg+1]); - iarg += 2; - } else error->all(FLERR,"Illegal min_modify command"); + } else { + int n = modify_param(narg-iarg,&arg[iarg]); + if (n == 0) error->all(FLERR,"Illegal fix_modify command"); + iarg += n; + } } } diff --git a/src/min.h b/src/min.h index ba1885671e..a63254231c 100644 --- a/src/min.h +++ b/src/min.h @@ -38,6 +38,7 @@ class Min : protected Pointers { int request(class Pair *, int, double); virtual bigint memory_usage() {return 0;} void modify_params(int, char **); + virtual int modify_param(int, char **) {return 0;} double fnorm_sqr(); double fnorm_inf(); @@ -58,11 +59,6 @@ class Min : protected Pointers { double dmax; // max dist to move any atom in one step int linestyle; // 0 = backtrack, 1 = quadratic, 2 = forcezero - // spinmin quantities - - double alpha_damp; // damping for spin minimization - double discret_factor; // factor for spin timestep evaluation - int nelist_global,nelist_atom; // # of PE,virial computes to check int nvlist_global,nvlist_atom; class Compute **elist_global; // lists of PE,virial Computes