From 09d0df43e215e6ce646334b5b60c3e8e53ff85f0 Mon Sep 17 00:00:00 2001 From: julient31 Date: Fri, 21 Feb 2020 17:53:14 -0700 Subject: [PATCH] Commit JT 022120 - added message for only one precession/spin (+doc) - added a per pair/spin class emag table --- doc/src/fix_precession_spin.rst | 7 +++++- src/SPIN/compute_spin.cpp | 18 +++++++++++++- src/SPIN/fix_precession_spin.cpp | 9 +++++++ src/SPIN/pair_spin.cpp | 6 +++++ src/SPIN/pair_spin.h | 4 +++ src/SPIN/pair_spin_dipole_cut.cpp | 10 ++++++++ src/SPIN/pair_spin_dipole_long.cpp | 24 ++++++++++++++---- src/SPIN/pair_spin_dmi.cpp | 15 ++++++++++- src/SPIN/pair_spin_exchange.cpp | 40 ++++++++++++++++++++++++++++++ src/SPIN/pair_spin_exchange.h | 2 ++ src/SPIN/pair_spin_magelec.cpp | 14 +++++++++++ src/SPIN/pair_spin_neel.cpp | 14 +++++++++++ 12 files changed, 155 insertions(+), 8 deletions(-) diff --git a/doc/src/fix_precession_spin.rst b/doc/src/fix_precession_spin.rst index 9cd15119bd..2d23ed1037 100644 --- a/doc/src/fix_precession_spin.rst +++ b/doc/src/fix_precession_spin.rst @@ -82,6 +82,7 @@ function for the same parameters. .. image:: JPG/zeeman_langevin.jpg :align: center + :width: 600 The temperature effects are accounted for by connecting the spin :math:`i` to a thermal bath using a Langevin thermostat (see @@ -159,11 +160,15 @@ No information about this fix is written to :doc:`binary restart files Restrictions """""""""""" - The *precession/spin* style is part of the SPIN package. This style is only enabled if LAMMPS was built with this package, and if the atom\_style "spin" was declared. See the :doc:`Build package ` doc page for more info. +The *precession/spin* style can only be declared once. If more +than one precession type (for example combining an anisotropy and a Zeeman interactions) +has to be declared, they have to be chained in the same command +line (as shown in the examples above). + Related commands """""""""""""""" diff --git a/src/SPIN/compute_spin.cpp b/src/SPIN/compute_spin.cpp index 0612e5720e..8a71be019b 100644 --- a/src/SPIN/compute_spin.cpp +++ b/src/SPIN/compute_spin.cpp @@ -69,6 +69,22 @@ void ComputeSpin::init() hbar = force->hplanck/MY_2PI; kb = force->boltz; + // loop 1: obtain # of Pairs, and # of Pair/Spin styles + + if (force->pair_match("spin",0,0)) { // only one Pair/Spin style + pair = force->pair_match("spin",0,0); + npairs = pair->instance_total; + npairspin = 1; + } else if (force->pair_match("spin",0,1)) { // more than one Pair/Spin style + pair = force->pair_match("spin",0,1); + npairs = pair->instance_total; + for (int i = 0; ipair_match("spin",0,i)) { + npairspin ++; + } + } + } + // init length of vector of ptrs to Pair/Spin styles if (npairspin > 0) { @@ -166,7 +182,7 @@ void ComputeSpin::compute_vector() if (pair_spin_flag) { for (int k = 0; k < npairspin; k++) { - // spin_pairs[k]->compute_single_pair(i,fmi); + magenergy += spin_pairs[k]->emag[i]; } } diff --git a/src/SPIN/fix_precession_spin.cpp b/src/SPIN/fix_precession_spin.cpp index f9307d7ad0..57e4549718 100644 --- a/src/SPIN/fix_precession_spin.cpp +++ b/src/SPIN/fix_precession_spin.cpp @@ -197,6 +197,15 @@ void FixPrecessionSpin::init() error->all(FLERR,"Illegal precession/spin command"); } + // check that fix precession/spin is only declared once + + int iprec = 0; + for (int iforce = 0; iforce < modify->nfix; iforce++) + if (strstr(modify->fix[iforce]->style,"precession/spin")) iprec++; + if (iprec > 1) + error->all(FLERR,"precession/spin command can only be declared once"); + + varflag = CONSTANT; if (magfieldstyle != CONSTANT) varflag = EQUAL; diff --git a/src/SPIN/pair_spin.cpp b/src/SPIN/pair_spin.cpp index f167e3455c..01b8775eab 100644 --- a/src/SPIN/pair_spin.cpp +++ b/src/SPIN/pair_spin.cpp @@ -29,6 +29,7 @@ #include "fix.h" #include "force.h" #include "math_const.h" +#include "memory.h" #include "modify.h" #include "neighbor.h" #include "neigh_request.h" @@ -98,4 +99,9 @@ void PairSpin::init_style() if (ifix >=0) lattice_flag = ((FixNVESpin *) modify->fix[ifix])->lattice_flag; + // test emag list storing mag energies + // init. size of energy stacking lists + + nlocal_max = atom->nlocal; + memory->grow(emag,nlocal_max,"pair/spin:emag"); } diff --git a/src/SPIN/pair_spin.h b/src/SPIN/pair_spin.h index 34f12d8d59..c8292236a3 100644 --- a/src/SPIN/pair_spin.h +++ b/src/SPIN/pair_spin.h @@ -31,6 +31,10 @@ friend class FixNVESpin; virtual void compute(int, int) {} virtual void compute_single_pair(int, double *) {} + + // test emag list storing mag energies + int nlocal_max; // max value of nlocal (for size of lists) + double *emag; // energy list protected: double hbar; // Planck constant (eV.ps.rad-1) diff --git a/src/SPIN/pair_spin_dipole_cut.cpp b/src/SPIN/pair_spin_dipole_cut.cpp index a7372b480d..6029f8bdbb 100644 --- a/src/SPIN/pair_spin_dipole_cut.cpp +++ b/src/SPIN/pair_spin_dipole_cut.cpp @@ -64,6 +64,9 @@ PairSpinDipoleCut::~PairSpinDipoleCut() memory->destroy(setflag); memory->destroy(cut_spin_long); memory->destroy(cutsq); + + // test emag list storing mag energies + memory->destroy(emag); } } @@ -185,6 +188,13 @@ void PairSpinDipoleCut::compute(int eflag, int vflag) numneigh = list->numneigh; firstneigh = list->firstneigh; + // test emag list storing mag energies + // checking size of emag + if (nlocal_max < nlocal) { // grow emag lists if necessary + nlocal_max = nlocal; + memory->grow(emag,nlocal_max,"pair/spin:emag"); + } + // computation of the exchange interaction // loop over atoms and their neighbors diff --git a/src/SPIN/pair_spin_dipole_long.cpp b/src/SPIN/pair_spin_dipole_long.cpp index 124522a9b9..7856035159 100644 --- a/src/SPIN/pair_spin_dipole_long.cpp +++ b/src/SPIN/pair_spin_dipole_long.cpp @@ -69,6 +69,9 @@ PairSpinDipoleLong::~PairSpinDipoleLong() memory->destroy(setflag); memory->destroy(cut_spin_long); memory->destroy(cutsq); + + // test emag list storing mag energies + memory->destroy(emag); } } @@ -212,6 +215,13 @@ void PairSpinDipoleLong::compute(int eflag, int vflag) numneigh = list->numneigh; firstneigh = list->firstneigh; + // test emag list storing mag energies + // checking size of emag + if (nlocal_max < nlocal) { // grow emag lists if necessary + nlocal_max = nlocal; + memory->grow(emag,nlocal_max,"pair/spin:emag"); + } + pre1 = 2.0 * g_ewald / MY_PIS; pre2 = 4.0 * pow(g_ewald,3.0) / MY_PIS; pre3 = 8.0 * pow(g_ewald,5.0) / MY_PIS; @@ -221,16 +231,20 @@ void PairSpinDipoleLong::compute(int eflag, int vflag) for (ii = 0; ii < inum; ii++) { i = ilist[ii]; + itype = type[i]; + + jlist = firstneigh[i]; + jnum = numneigh[i]; xi[0] = x[i][0]; xi[1] = x[i][1]; xi[2] = x[i][2]; - jlist = firstneigh[i]; - jnum = numneigh[i]; spi[0] = sp[i][0]; spi[1] = sp[i][1]; spi[2] = sp[i][2]; spi[3] = sp[i][3]; - itype = type[i]; + + // test emag list storing mag energies + emag[i] = 0.0; for (jj = 0; jj < jnum; jj++) { j = jlist[jj]; @@ -294,9 +308,9 @@ void PairSpinDipoleLong::compute(int eflag, int vflag) if (eflag) { if (rsq <= local_cut2) { - evdwl -= spi[0]*fmi[0] + spi[1]*fmi[1] + - spi[2]*fmi[2]; + evdwl -= (spi[0]*fmi[0] + spi[1]*fmi[1] + spi[2]*fmi[2]); evdwl *= 0.5*hbar; + emag[i] += evdwl; } } else evdwl = 0.0; diff --git a/src/SPIN/pair_spin_dmi.cpp b/src/SPIN/pair_spin_dmi.cpp index 04c2dc408d..266bc05da4 100644 --- a/src/SPIN/pair_spin_dmi.cpp +++ b/src/SPIN/pair_spin_dmi.cpp @@ -53,6 +53,9 @@ PairSpinDmi::~PairSpinDmi() memory->destroy(vmech_dmy); memory->destroy(vmech_dmz); memory->destroy(cutsq); + + // test emag list storing mag energies + memory->destroy(emag); } } @@ -191,6 +194,13 @@ void PairSpinDmi::compute(int eflag, int vflag) numneigh = list->numneigh; firstneigh = list->firstneigh; + // test emag list storing mag energies + // checking size of emag + if (nlocal_max < nlocal) { // grow emag lists if necessary + nlocal_max = nlocal; + memory->grow(emag,nlocal_max,"pair/spin:emag"); + } + // dmi computation // loop over all atoms @@ -206,7 +216,9 @@ void PairSpinDmi::compute(int eflag, int vflag) spi[0] = sp[i][0]; spi[1] = sp[i][1]; spi[2] = sp[i][2]; - + + // test emag list storing mag energies + emag[i] = 0.0; // loop on neighbors @@ -260,6 +272,7 @@ void PairSpinDmi::compute(int eflag, int vflag) if (eflag) { evdwl -= (spi[0]*fmi[0] + spi[1]*fmi[1] + spi[2]*fmi[2]); evdwl *= 0.5*hbar; + emag[i] += evdwl; } else evdwl = 0.0; if (evflag) ev_tally_xyz(i,j,nlocal,newton_pair, diff --git a/src/SPIN/pair_spin_exchange.cpp b/src/SPIN/pair_spin_exchange.cpp index 6eacb04ee3..d645515506 100644 --- a/src/SPIN/pair_spin_exchange.cpp +++ b/src/SPIN/pair_spin_exchange.cpp @@ -50,6 +50,9 @@ PairSpinExchange::~PairSpinExchange() memory->destroy(J2); memory->destroy(J3); memory->destroy(cutsq); // to be implemented + + // test emag list storing mag energies + memory->destroy(emag); } } @@ -176,6 +179,13 @@ void PairSpinExchange::compute(int eflag, int vflag) numneigh = list->numneigh; firstneigh = list->firstneigh; + // test emag list storing mag energies + // checking size of emag + if (nlocal_max < nlocal) { // grow emag lists if necessary + nlocal_max = nlocal; + memory->grow(emag,nlocal_max,"pair/spin:emag"); + } + // computation of the exchange interaction // loop over atoms and their neighbors @@ -191,6 +201,9 @@ void PairSpinExchange::compute(int eflag, int vflag) spi[0] = sp[i][0]; spi[1] = sp[i][1]; spi[2] = sp[i][2]; + + // test emag list storing mag energies + emag[i] = 0.0; // loop on neighbors @@ -243,6 +256,10 @@ void PairSpinExchange::compute(int eflag, int vflag) if (eflag) { evdwl -= (spi[0]*fmi[0] + spi[1]*fmi[1] + spi[2]*fmi[2]); evdwl *= 0.5*hbar; + // printf("test ex energy: %g \n",evdwl); + // evdwl = -0.5*compute_energy(i,j,rsq,spi,spj); + // printf("test ex energy: %g \n",evdwl); + emag[i] += evdwl; // evdwl *= hbar; } else evdwl = 0.0; @@ -385,6 +402,29 @@ void PairSpinExchange::compute_exchange_mech(int i, int j, double rsq, double ei fi[2] -= Jex_mech*eij[2]; } +/* ---------------------------------------------------------------------- + compute energy of spin pair i and j +------------------------------------------------------------------------- */ + +// double PairSpinExchange::compute_energy(int i, int j, double rsq, double spi[3], double spj[3]) +// { +// int *type = atom->type; +// int itype, jtype; +// double Jex, ra; +// double energy = 0.0; +// itype = type[i]; +// jtype = type[j]; +// +// Jex = J1_mech[itype][jtype]; +// ra = rsq/J3[itype][jtype]/J3[itype][jtype]; +// Jex = 4.0*Jex*ra; +// Jex *= (1.0-J2[itype][jtype]*ra); +// Jex *= exp(-ra); +// +// energy = Jex*(spi[0]*spj[0]+spi[1]*spj[1]+spi[2]*spj[2]); +// return energy; +// } + /* ---------------------------------------------------------------------- allocate all arrays ------------------------------------------------------------------------- */ diff --git a/src/SPIN/pair_spin_exchange.h b/src/SPIN/pair_spin_exchange.h index 19eafeb5ca..4e9e6bfac8 100644 --- a/src/SPIN/pair_spin_exchange.h +++ b/src/SPIN/pair_spin_exchange.h @@ -39,6 +39,8 @@ class PairSpinExchange : public PairSpin { void compute_exchange(int, int, double, double *, double *); void compute_exchange_mech(int, int, double, double *, double *, double *, double *); + // double compute_energy(int , int , double , double *, double *); + void write_restart(FILE *); void read_restart(FILE *); void write_restart_settings(FILE *); diff --git a/src/SPIN/pair_spin_magelec.cpp b/src/SPIN/pair_spin_magelec.cpp index fabad4ae4d..ef91ab764a 100644 --- a/src/SPIN/pair_spin_magelec.cpp +++ b/src/SPIN/pair_spin_magelec.cpp @@ -51,6 +51,9 @@ PairSpinMagelec::~PairSpinMagelec() memory->destroy(v_mey); memory->destroy(v_mez); memory->destroy(cutsq); // to be deteled + + // test emag list storing mag energies + memory->destroy(emag); } } @@ -185,6 +188,13 @@ void PairSpinMagelec::compute(int eflag, int vflag) numneigh = list->numneigh; firstneigh = list->firstneigh; + // test emag list storing mag energies + // checking size of emag + if (nlocal_max < nlocal) { // grow emag lists if necessary + nlocal_max = nlocal; + memory->grow(emag,nlocal_max,"pair/spin:emag"); + } + // magneto-electric computation // loop over atoms and their neighbors @@ -200,6 +210,9 @@ void PairSpinMagelec::compute(int eflag, int vflag) spi[0] = sp[i][0]; spi[1] = sp[i][1]; spi[2] = sp[i][2]; + + // test emag list storing mag energies + emag[i] = 0.0; // loop on neighbors @@ -252,6 +265,7 @@ void PairSpinMagelec::compute(int eflag, int vflag) if (eflag) { evdwl -= (spi[0]*fmi[0] + spi[1]*fmi[1] + spi[2]*fmi[2]); evdwl *= 0.5*hbar; + emag[i] += evdwl; } else evdwl = 0.0; if (evflag) ev_tally_xyz(i,j,nlocal,newton_pair, diff --git a/src/SPIN/pair_spin_neel.cpp b/src/SPIN/pair_spin_neel.cpp index 6cf5d4843f..e158906b75 100644 --- a/src/SPIN/pair_spin_neel.cpp +++ b/src/SPIN/pair_spin_neel.cpp @@ -54,6 +54,9 @@ PairSpinNeel::~PairSpinNeel() memory->destroy(q2); memory->destroy(q3); memory->destroy(cutsq); // to be deleted + + // test emag list storing mag energies + memory->destroy(emag); } } @@ -190,6 +193,13 @@ void PairSpinNeel::compute(int eflag, int vflag) numneigh = list->numneigh; firstneigh = list->firstneigh; + // test emag list storing mag energies + // checking size of emag + if (nlocal_max < nlocal) { // grow emag lists if necessary + nlocal_max = nlocal; + memory->grow(emag,nlocal_max,"pair/spin:emag"); + } + // computation of the neel interaction // loop over atoms and their neighbors @@ -206,6 +216,9 @@ void PairSpinNeel::compute(int eflag, int vflag) spi[1] = sp[i][1]; spi[2] = sp[i][2]; + // test emag list storing mag energies + emag[i] = 0.0; + // loop on neighbors for (jj = 0; jj < jnum; jj++) { @@ -262,6 +275,7 @@ void PairSpinNeel::compute(int eflag, int vflag) // evdwl = (spi[0]*fmi[0] + spi[1]*fmi[1] + spi[2]*fmi[2]); evdwl = compute_neel_energy(i,j,rsq,eij,spi,spj); evdwl *= 0.5*hbar; + emag[i] += evdwl; } else evdwl = 0.0; if (evflag) ev_tally_xyz(i,j,nlocal,newton_pair,