From 75ddde438c13cf32564a9271c03d6867c5db92bc Mon Sep 17 00:00:00 2001 From: julient31 Date: Tue, 12 Mar 2019 14:38:49 -0600 Subject: [PATCH] Commit JT 031219 - correct errors in fix_prec_spin - clean version of spinmin --- src/REPLICA/fix_neb_spin.cpp | 3 +- src/SPIN/fix_precession_spin.cpp | 10 +- src/SPIN/min_spinmin.cpp | 242 +++++++++---------------------- src/SPIN/min_spinmin.h | 4 +- src/min.cpp | 8 + src/min.h | 5 + 6 files changed, 86 insertions(+), 186 deletions(-) diff --git a/src/REPLICA/fix_neb_spin.cpp b/src/REPLICA/fix_neb_spin.cpp index 42450c2f0f..7a72936200 100644 --- a/src/REPLICA/fix_neb_spin.cpp +++ b/src/REPLICA/fix_neb_spin.cpp @@ -337,7 +337,8 @@ void FixNEB_spin::min_post_force(int /*vflag*/) //printf("test veng: %g / %g / %g \n",veng,vprev,vnext); //error->universe_all(FLERR,"End test"); - if (FreeEndFinal && ireplica == nreplica-1 && (update->ntimestep == 0)) EFinalIni = veng; + if (FreeEndFinal && ireplica == nreplica-1 && (update->ntimestep == 0)) + EFinalIni = veng; if (ireplica == 0) vIni=veng; diff --git a/src/SPIN/fix_precession_spin.cpp b/src/SPIN/fix_precession_spin.cpp index 08e4fd5a63..433a260e83 100644 --- a/src/SPIN/fix_precession_spin.cpp +++ b/src/SPIN/fix_precession_spin.cpp @@ -173,8 +173,6 @@ void FixPrecessionSpin::setup(int vflag) void FixPrecessionSpin::post_force(int /*vflag*/) { - printf("test inside post force (precession) \n"); - // update mag field with time (potential improvement) if (varflag != CONSTANT) { @@ -204,7 +202,7 @@ void FixPrecessionSpin::post_force(int /*vflag*/) if (aniso_flag) { // compute magnetic anisotropy compute_anisotropy(spi,fmi); - emag -= (spi[0]*fmi[0] + spi[1]*fmi[1] + spi[2]*fmi[2]); + emag -= 0.5*(spi[0]*fmi[0] + spi[1]*fmi[1] + spi[2]*fmi[2]); } fm[i][0] += fmi[0]; @@ -232,9 +230,9 @@ 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][3]*hx; - fmi[1] -= sp[i][3]*hy; - fmi[2] -= sp[i][3]*hz; + fmi[0] += sp[i][0]*hx; + fmi[1] += sp[i][1]*hy; + fmi[2] += sp[i][3]*hz; } /* ---------------------------------------------------------------------- */ diff --git a/src/SPIN/min_spinmin.cpp b/src/SPIN/min_spinmin.cpp index 2e03086bf0..cdd0d45287 100644 --- a/src/SPIN/min_spinmin.cpp +++ b/src/SPIN/min_spinmin.cpp @@ -46,14 +46,13 @@ MinSpinMin::MinSpinMin(LAMMPS *lmp) : Min(lmp) {} void MinSpinMin::init() { + alpha_damp = 1.0; + discret_factor = 10.0; + Min::init(); - dt = update->dt; + dts = dt = update->dt; last_negative = update->ntimestep; - - // test dts - dts = dt; - } /* ---------------------------------------------------------------------- */ @@ -94,8 +93,7 @@ void MinSpinMin::reset_vectors() int MinSpinMin::iterate(int maxiter) { bigint ntimestep; - //double vmax,vdotf,vdotfall,fdotf,fdotfall,scale; - //double dtvone,dtv,dtf,dtfm; + double fmdotfm,fmdotfmall; int flag,flagall; for (int iter = 0; iter < maxiter; iter++) { @@ -114,106 +112,6 @@ int MinSpinMin::iterate(int maxiter) advance_spins(dts); - - //// zero velocity if anti-parallel to force - //// else project velocity in direction of force - - //double **v = atom->v; - //double **f = atom->f; - //int nlocal = atom->nlocal; - - //vdotf = 0.0; - //for (int i = 0; i < nlocal; i++) - // vdotf += v[i][0]*f[i][0] + v[i][1]*f[i][1] + v[i][2]*f[i][2]; - //MPI_Allreduce(&vdotf,&vdotfall,1,MPI_DOUBLE,MPI_SUM,world); - - // sum vdotf over replicas, if necessary - // this communicator would be invalid for multiprocess replicas - - //if (update->multireplica == 1) { - // vdotf = vdotfall; - // MPI_Allreduce(&vdotf,&vdotfall,1,MPI_DOUBLE,MPI_SUM,universe->uworld); - //} - - //if (vdotfall < 0.0) { - // last_negative = ntimestep; - // for (int i = 0; i < nlocal; i++) - // v[i][0] = v[i][1] = v[i][2] = 0.0; - - //} else { - // fdotf = 0.0; - // for (int i = 0; i < nlocal; i++) - // fdotf += f[i][0]*f[i][0] + f[i][1]*f[i][1] + f[i][2]*f[i][2]; - // MPI_Allreduce(&fdotf,&fdotfall,1,MPI_DOUBLE,MPI_SUM,world); - - // // sum fdotf over replicas, if necessary - // // this communicator would be invalid for multiprocess replicas - - // if (update->multireplica == 1) { - // fdotf = fdotfall; - // MPI_Allreduce(&fdotf,&fdotfall,1,MPI_DOUBLE,MPI_SUM,universe->uworld); - // } - - // if (fdotfall == 0.0) scale = 0.0; - // else scale = vdotfall/fdotfall; - // for (int i = 0; i < nlocal; i++) { - // v[i][0] = scale * f[i][0]; - // v[i][1] = scale * f[i][1]; - // v[i][2] = scale * f[i][2]; - // } - //} - - //// limit timestep so no particle moves further than dmax - - //double *rmass = atom->rmass; - //double *mass = atom->mass; - //int *type = atom->type; - - //dtvone = dt; - - //for (int i = 0; i < nlocal; i++) { - // vmax = MAX(fabs(v[i][0]),fabs(v[i][1])); - // vmax = MAX(vmax,fabs(v[i][2])); - // if (dtvone*vmax > dmax) dtvone = dmax/vmax; - //} - //MPI_Allreduce(&dtvone,&dtv,1,MPI_DOUBLE,MPI_MIN,world); - - //// min dtv over replicas, if necessary - //// this communicator would be invalid for multiprocess replicas - - //if (update->multireplica == 1) { - // dtvone = dtv; - // MPI_Allreduce(&dtvone,&dtv,1,MPI_DOUBLE,MPI_MIN,universe->uworld); - //} - - //dtf = dtv * force->ftm2v; - - //// Euler integration step - - //double **x = atom->x; - - //if (rmass) { - // for (int i = 0; i < nlocal; i++) { - // dtfm = dtf / rmass[i]; - // x[i][0] += dtv * v[i][0]; - // x[i][1] += dtv * v[i][1]; - // x[i][2] += dtv * v[i][2]; - // v[i][0] += dtfm * f[i][0]; - // v[i][1] += dtfm * f[i][1]; - // v[i][2] += dtfm * f[i][2]; - // } - //} else { - // for (int i = 0; i < nlocal; i++) { - // dtfm = dtf / mass[type[i]]; - // x[i][0] += dtv * v[i][0]; - // x[i][1] += dtv * v[i][1]; - // x[i][2] += dtv * v[i][2]; - // v[i][0] += dtfm * f[i][0]; - // v[i][1] += dtfm * f[i][1]; - // v[i][2] += dtfm * f[i][2]; - // } - //} - eprevious = ecurrent; ecurrent = energy_force(0); neval++; @@ -237,37 +135,22 @@ int MinSpinMin::iterate(int maxiter) } } - //// magnetic force tolerance criterion - //// sync across replicas if running multi-replica minimization + // magnetic torque tolerance criterion + // sync across replicas if running multi-replica minimization - //if (update->fmtol > 0.0) { - // fmdotfm = fmnorm_sqr(); - // if (update->multireplica == 0) { - // if (fmdotfm < update->fmtol*update->fmtol) return FTOL; - // } else { - // if (fmdotfm < update->fmtol*update->fmtol) flag = 0; - // else flag = 1; - // MPI_Allreduce(&fmlag,&fmlagall,1,MPI_INT,MPI_SUM,universe->uworld); - // if (fmlagall == 0) return FTOL; - // } - //} - - //// force tolerance criterion - //// sync across replicas if running multi-replica minimization + if (update->ftol > 0.0) { + fmdotfm = fmnorm_sqr(); + if (update->multireplica == 0) { + if (fmdotfm < update->ftol*update->ftol) return FTOL; + } else { + if (fmdotfm < update->ftol*update->ftol) flag = 0; + else flag = 1; + MPI_Allreduce(&flag,&flagall,1,MPI_INT,MPI_SUM,universe->uworld); + if (flagall == 0) return FTOL; + } + } - //if (update->ftol > 0.0) { - // fdotf = fnorm_sqr(); - // if (update->multireplica == 0) { - // if (fdotf < update->ftol*update->ftol) return FTOL; - // } else { - // if (fdotf < update->ftol*update->ftol) flag = 0; - // else flag = 1; - // MPI_Allreduce(&flag,&flagall,1,MPI_INT,MPI_SUM,universe->uworld); - // if (flagall == 0) return FTOL; - // } - //} - - //// output for thermo, dump, restart files + // output for thermo, dump, restart files if (output->next == ntimestep) { timer->stamp(); @@ -299,12 +182,6 @@ double MinSpinMin::evaluate_dt() fmsq = fm[i][0]*fm[i][0]+fm[i][1]*fm[i][1]+fm[i][2]*fm[i][2]; fmaxsqone = MAX(fmaxsqone,fmsq); } - //printf("test inside evaluate dt, fmaxsqone = %g \n",fmaxsqone); - //for (int i = 0; i < nlocal; i++) - // if (mask[i] & groupbit) { - // fmsq = fm[i][0]*fm[i][0]+fm[i][1]*fm[i][1]+fm[i][2]*fm[i][2]; - // fmaxsqone = MAX(fmaxsqone,fmsq); - // } // finding max fm on this replica @@ -320,15 +197,13 @@ double MinSpinMin::evaluate_dt() MPI_Allreduce(&fmaxsqloc,&fmaxsqall,1,MPI_DOUBLE,MPI_MAX,universe->uworld); } - //if (fmaxsqall < fmaxsqloc) - // error->all(FLERR,"Incorrect fmaxall calc."); + if (fmaxsqall == 0.0) + error->all(FLERR,"Incorrect fmaxsqall calculation"); - // define max timestep - // dividing by 10 the inverse of max frequency - //printf("test inside evaluate dt, fmaxsqall = %g \n",fmaxsqall); + // define max timestep by dividing by the + // inverse of max frequency by discret_factor - dtmax = MY_2PI/(10.0*sqrt(fmaxsqall)); - //printf("test inside evaluate dt, dtmax = %g \n",dtmax); + dtmax = MY_2PI/(discret_factor*sqrt(fmaxsqall)); return dtmax; } @@ -345,28 +220,19 @@ void MinSpinMin::advance_spins(double dts) double **fm = atom->fm; double tdampx,tdampy,tdampz; double msq,scale,fm2,energy,dts2; - double alpha; double cp[3],g[3]; dts2 = dts*dts; - // fictitious Gilbert damping of 1 - - alpha = 1.0; - - // loop on all spins on proc. - //if (ireplica != nreplica-1 && ireplica != 0) - // for (int i = 0; i < nlocal; i++) - // if (mask[i] & groupbit) { for (int i = 0; i < nlocal; i++) { // calc. damping torque - tdampx = -alpha*(fm[i][1]*sp[i][2] - fm[i][2]*sp[i][1]); - tdampy = -alpha*(fm[i][2]*sp[i][0] - fm[i][0]*sp[i][2]); - tdampz = -alpha*(fm[i][0]*sp[i][1] - fm[i][1]*sp[i][0]); + tdampx = -alpha_damp*(fm[i][1]*sp[i][2] - fm[i][2]*sp[i][1]); + tdampy = -alpha_damp*(fm[i][2]*sp[i][0] - fm[i][0]*sp[i][2]); + tdampz = -alpha_damp*(fm[i][0]*sp[i][1] - fm[i][1]*sp[i][0]); // apply advance algorithm (geometric, norm preserving) @@ -401,22 +267,44 @@ void MinSpinMin::advance_spins(double dts) sp[i][1] *= scale; sp[i][2] *= scale; - // comm. sp[i] to atoms with same tag (for serial algo) - - // no need for simplecticity - //if (sector_flag == 0) { - // if (sametag[i] >= 0) { - // j = sametag[i]; - // while (j >= 0) { - // sp[j][0] = sp[i][0]; - // sp[j][1] = sp[i][1]; - // sp[j][2] = sp[i][2]; - // j = sametag[j]; - // } - // } - //} - // + // no comm. to atoms with same tag + // because no need for simplecticity } - +} + +/* ---------------------------------------------------------------------- + compute and return ||mag. torque||_2^2 +------------------------------------------------------------------------- */ + +double MinSpinMin::fmnorm_sqr() +{ + int i,n; + double *fmatom; + + 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; } diff --git a/src/SPIN/min_spinmin.h b/src/SPIN/min_spinmin.h index 943b2d2749..abc532a3d5 100644 --- a/src/SPIN/min_spinmin.h +++ b/src/SPIN/min_spinmin.h @@ -34,12 +34,12 @@ class MinSpinMin : public Min { int iterate(int); double evaluate_dt(); void advance_spins(double); - - class FixNEB_spin *fneb; + double fmnorm_sqr(); private: // global and spin timesteps + double dt; double dts; diff --git a/src/min.cpp b/src/min.cpp index cd9253f8d3..c75db6e2b0 100644 --- a/src/min.cpp +++ b/src/min.cpp @@ -655,6 +655,14 @@ 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"); } } diff --git a/src/min.h b/src/min.h index 92da97c664..ba1885671e 100644 --- a/src/min.h +++ b/src/min.h @@ -58,6 +58,11 @@ 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