diff --git a/doc/src/fix_precession_spin.rst b/doc/src/fix_precession_spin.rst index 5e818374a0..183753db5f 100644 --- a/doc/src/fix_precession_spin.rst +++ b/doc/src/fix_precession_spin.rst @@ -12,7 +12,7 @@ Syntax * ID, group are documented in :doc:`fix ` command * precession/spin = style name of this fix command -* style = *zeeman* or *anisotropy* or *cubic* +* style = *zeeman* or *anisotropy* or *cubic* or *stt* .. parsed-literal:: @@ -22,12 +22,12 @@ Syntax *anisotropy* args = K x y z K = intensity of the magnetic anisotropy (in eV) x y z = vector direction of the anisotropy - - .. parsed-literal:: - *cubic* args = K1 K2c n1x n1y n1x n2x n2y n2z n3x n3y n3z K1 and K2c = intensity of the magnetic anisotropy (in eV) n1x to n3z = three direction vectors of the cubic anisotropy + *stt* args = J x y z + J = intensity of the spin-transfer torque field + x y z = vector direction of the field Examples """""""" @@ -125,6 +125,11 @@ axis along the :math:`(1 1 1)`-type cube diagonals). :math:`K_2^c > diagonals. See chapter 2 of :ref:`(Skomski) ` for more details on cubic anisotropies. +Style *stt* is used to simulate the interaction between the spins and +a spin-transfer torque. +See equation (7) of :ref:`(Chirac) ` for more details about the +implemented spin-transfer torque term. + In all cases, the choice of :math:`(x y z)` only imposes the vector directions for the forces. Only the direction of the vector is important; its length is ignored (the entered vectors are @@ -132,6 +137,16 @@ normalized). Those styles can be combined within one single command line. +.. note:: + + The norm of all vectors defined with the precession/spin command + have to be non-zero. For example, defining + "fix 1 all precession/spin zeeman 0.1 0.0 0.0 0.0" would result + in an error message. + Since those vector components are used to compute the inverse of the + field (or anisotropy) vector norm, setting a zero-vector would result + in a division by zero. + ---------- Restart, fix_modify, output, run start/stop, minimize info @@ -162,11 +177,6 @@ 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 """""""""""""""" @@ -184,3 +194,9 @@ none **(Skomski)** Skomski, R. (2008). Simple models of magnetism. Oxford University Press. + +.. _Chirac1: + +**(Chirac)** Chirac, Théophile, et al. Ultrafast antiferromagnetic +switching in NiO induced by spin transfer torques. +Physical Review B 102.13 (2020): 134415. diff --git a/examples/SPIN/test_problems/README b/examples/SPIN/test_problems/README index 0a1362ec9c..17e0935a35 100644 --- a/examples/SPIN/test_problems/README +++ b/examples/SPIN/test_problems/README @@ -45,3 +45,14 @@ directory. results (computed by the python script). Note: This example is a reworked version of a test problem provided by Martin Kroger (ETHZ). + +- validation_nve: + simulates a small assembly of magnetic atoms (54). The atoms are + coupled by an exchange interaction and a mechanical potential + (EAM here). + This example represents an NVE run: the total energy of the + system is preserved, whereas the spin and lattice energy + reservoirs are exchanging energy. + Run as: ./run-test-nve.sh + Output: res_lammps.dat contains the data. The results are displayed + by nve_spin_lattice.pdf. diff --git a/examples/SPIN/test_problems/run_all.sh b/examples/SPIN/test_problems/run_all.sh new file mode 100755 index 0000000000..fb682e3ef5 --- /dev/null +++ b/examples/SPIN/test_problems/run_all.sh @@ -0,0 +1,25 @@ +#!/bin/bash + +# test 1: damping and exchange +cd validation_damped_exchange/ +./run-test-exchange.sh +rm dump.data res_lammps.dat res_llg.dat +cd .. + +# test 2: damping and Zeeman +cd validation_damped_precession/ +./run-test-prec.sh +rm res_lammps.dat res_llg.dat +cd .. + +# test 3: langevin, damping and Zeeman +cd validation_langevin_precession/ +./run-test-prec.sh +rm average_spin test-prec-spin.in res_lammps.dat res_langevin.dat +cd .. + +# test 4: NVE run, test Etot preservation +cd validation_nve/ +./run-test-nve.sh +rm nve_spin_lattice.pdf res_lammps.dat +cd .. diff --git a/examples/SPIN/test_problems/validation_nvt/in.spin.nvt_lattice b/examples/SPIN/test_problems/validation_nvt/in.spin.nvt_lattice index cc7bcfa68d..a1f739fe54 100644 --- a/examples/SPIN/test_problems/validation_nvt/in.spin.nvt_lattice +++ b/examples/SPIN/test_problems/validation_nvt/in.spin.nvt_lattice @@ -31,8 +31,8 @@ neighbor 0.1 bin neigh_modify every 10 check yes delay 20 fix 1 all precession/spin zeeman 0.0 0.0 0.0 1.0 -fix 2 all langevin 200.0 200.0 1.0 48279 -fix 3 all langevin/spin 0.0 0.00001 321 +fix 2 all langevin 200.0 200.0 0.1 48279 +fix 3 all langevin/spin 0.0 0.0 321 fix 4 all nve/spin lattice moving timestep 0.001 diff --git a/examples/SPIN/test_problems/validation_nvt/in.spin.nvt_spin b/examples/SPIN/test_problems/validation_nvt/in.spin.nvt_spin index 7e8152f481..102f76a2a6 100644 --- a/examples/SPIN/test_problems/validation_nvt/in.spin.nvt_spin +++ b/examples/SPIN/test_problems/validation_nvt/in.spin.nvt_spin @@ -22,16 +22,19 @@ mass 1 55.845 set group all spin 2.2 0.0 0.0 1.0 velocity all create 0 4928459 rot yes dist gaussian -pair_style hybrid/overlay eam/alloy spin/exchange 3.5 +# pair_style hybrid/overlay eam/alloy spin/exchange 3.5 +pair_style hybrid/overlay eam/alloy spin/exchange 4.0 spin/neel 4.0 pair_coeff * * eam/alloy Fe_Mishin2006.eam.alloy Fe pair_coeff * * spin/exchange exchange 3.4 0.1 0.2171 1.841 +pair_coeff * * spin/neel neel 4.0 0.02 0.0 1.841 0.0 0.0 1.0 neighbor 0.1 bin neigh_modify every 10 check yes delay 20 fix 1 all precession/spin zeeman 0.0 0.0 0.0 1.0 -fix 2 all langevin/spin 200.0 0.01 321 -fix 3 all nve/spin lattice moving +fix 2 all langevin 0.0 0.0 0.0 48279 +fix 3 all langevin/spin 200.0 0.01 321 +fix 4 all nve/spin lattice moving timestep 0.001 # compute and output options diff --git a/src/SPIN/fix_langevin_spin.cpp b/src/SPIN/fix_langevin_spin.cpp index c9cb0fa2ec..949833caa7 100644 --- a/src/SPIN/fix_langevin_spin.cpp +++ b/src/SPIN/fix_langevin_spin.cpp @@ -24,9 +24,11 @@ #include "fix_langevin_spin.h" #include #include +#include "atom.h" #include "comm.h" #include "error.h" #include "force.h" +#include "group.h" #include "math_const.h" #include "memory.h" #include "modify.h" @@ -42,7 +44,7 @@ using namespace MathConst; /* ---------------------------------------------------------------------- */ FixLangevinSpin::FixLangevinSpin(LAMMPS *lmp, int narg, char **arg) : - Fix(lmp, narg, arg), id_temp(nullptr), random(nullptr) + Fix(lmp, narg, arg), random(nullptr) { if (narg != 6) error->all(FLERR,"Illegal langevin/spin command"); @@ -106,10 +108,8 @@ void FixLangevinSpin::init() double hbar = force->hplanck/MY_2PI; // eV/(rad.THz) double kb = force->boltz; // eV/K - // D = (MY_2PI*alpha_t*gil_factor*kb*temp); D = (alpha_t*gil_factor*kb*temp); - // D = (12.0/MY_2PI)*(MY_2PI*alpha_t*gil_factor*kb*temp); D /= (hbar*dts); sigma = sqrt(2.0*D); } @@ -142,7 +142,7 @@ void FixLangevinSpin::add_tdamping(double spi[3], double fmi[3]) /* ---------------------------------------------------------------------- */ -void FixLangevinSpin::add_temperature(double fmi[3]) +void FixLangevinSpin::add_temperature(int i, double spi[3], double fmi[3]) { // double rx = sigma*(2.0*random->uniform() - 1.0); // double ry = sigma*(2.0*random->uniform() - 1.0); @@ -150,6 +150,7 @@ void FixLangevinSpin::add_temperature(double fmi[3]) double rx = sigma*random->gaussian(); double ry = sigma*random->gaussian(); double rz = sigma*random->gaussian(); + double hbar = force->hplanck/MY_2PI; // adding the random field @@ -163,3 +164,14 @@ void FixLangevinSpin::add_temperature(double fmi[3]) fmi[1] *= gil_factor; fmi[2] *= gil_factor; } + +/* ---------------------------------------------------------------------- */ + +void FixLangevinSpin::compute_single_langevin(int i, double spi[3], double fmi[3]) +{ + int *mask = atom->mask; + if (mask[i] & groupbit) { + if (tdamp_flag) add_tdamping(spi,fmi); + if (temp_flag) add_temperature(i,spi,fmi); + } +} diff --git a/src/SPIN/fix_langevin_spin.h b/src/SPIN/fix_langevin_spin.h index c73b33353b..090e5b666a 100644 --- a/src/SPIN/fix_langevin_spin.h +++ b/src/SPIN/fix_langevin_spin.h @@ -26,7 +26,7 @@ namespace LAMMPS_NS { class FixLangevinSpin : public Fix { public: - int tdamp_flag,ldamp_flag,temp_flag; // damping and temperature flags + int tdamp_flag,temp_flag; // damping and temperature flags FixLangevinSpin(class LAMMPS *, int, char **); virtual ~FixLangevinSpin(); @@ -34,7 +34,8 @@ class FixLangevinSpin : public Fix { void init(); void setup(int); void add_tdamping(double *, double *); // add transverse damping - void add_temperature(double *); // add temperature + void add_temperature(int, double *, double *); + void compute_single_langevin(int, double *, double *); protected: double alpha_t; // transverse mag. damping @@ -43,9 +44,6 @@ class FixLangevinSpin : public Fix { double D,sigma; // bath intensity var. double gil_factor; // gilbert's prefactor - char *id_temp; - class Compute *temperature; - int nlevels_respa; class RanMars *random; int seed; diff --git a/src/SPIN/fix_nve_spin.cpp b/src/SPIN/fix_nve_spin.cpp index eb6d918ab0..5da7c32b6c 100644 --- a/src/SPIN/fix_nve_spin.cpp +++ b/src/SPIN/fix_nve_spin.cpp @@ -60,7 +60,8 @@ enum{NONE}; FixNVESpin::FixNVESpin(LAMMPS *lmp, int narg, char **arg) : Fix(lmp, narg, arg), - pair(nullptr), spin_pairs(nullptr), + pair(nullptr), spin_pairs(nullptr), locklangevinspin(nullptr), + locksetforcespin(nullptr), lockprecessionspin(nullptr), rsec(nullptr), stack_head(nullptr), stack_foot(nullptr), backward_stacks(nullptr), forward_stacks(nullptr) { @@ -74,6 +75,9 @@ FixNVESpin::FixNVESpin(LAMMPS *lmp, int narg, char **arg) : nlocal_max = 0; npairs = 0; npairspin = 0; + + // test nprec + nprecspin = nlangspin = nsetspin = 0; // checking if map array or hash is defined @@ -125,7 +129,6 @@ FixNVESpin::FixNVESpin(LAMMPS *lmp, int narg, char **arg) : maglangevin_flag = 0; tdamp_flag = temp_flag = 0; setforce_spin_flag = 0; - } /* ---------------------------------------------------------------------- */ @@ -138,6 +141,8 @@ FixNVESpin::~FixNVESpin() memory->destroy(forward_stacks); memory->destroy(backward_stacks); delete [] spin_pairs; + delete [] locklangevinspin; + delete [] lockprecessionspin; } /* ---------------------------------------------------------------------- */ @@ -189,20 +194,20 @@ void FixNVESpin::init() // loop 2: fill vector with ptrs to Pair/Spin styles - int count = 0; + int count1 = 0; if (npairspin == 1) { - count = 1; + count1 = 1; spin_pairs[0] = (PairSpin *) force->pair_match("spin",0,0); } else if (npairspin > 1) { for (int i = 0; ipair_match("spin",0,i)) { - spin_pairs[count] = (PairSpin *) force->pair_match("spin",0,i); - count++; + spin_pairs[count1] = (PairSpin *) force->pair_match("spin",0,i); + count1++; } } } - if (count != npairspin) + if (count1 != npairspin) error->all(FLERR,"Incorrect number of spin pairs"); // set pair/spin and long/spin flags @@ -215,30 +220,71 @@ void FixNVESpin::init() } } - // ptrs FixPrecessionSpin classes - + // set ptrs for fix precession/spin styles + + // loop 1: obtain # of fix precession/spin styles + int iforce; for (iforce = 0; iforce < modify->nfix; iforce++) { if (strstr(modify->fix[iforce]->style,"precession/spin")) { - precession_spin_flag = 1; - lockprecessionspin = (FixPrecessionSpin *) modify->fix[iforce]; + nprecspin++; } } + + // init length of vector of ptrs to precession/spin styles - // ptrs on the FixLangevinSpin class + if (nprecspin > 0) { + lockprecessionspin = new FixPrecessionSpin*[nprecspin]; + } + + // loop 2: fill vector with ptrs to precession/spin styles + int count2 = 0; + if (nprecspin > 0) { + for (iforce = 0; iforce < modify->nfix; iforce++) { + if (strstr(modify->fix[iforce]->style,"precession/spin")) { + precession_spin_flag = 1; + lockprecessionspin[count2] = (FixPrecessionSpin *) modify->fix[iforce]; + count2++; + } + } + } + + if (count2 != nprecspin) + error->all(FLERR,"Incorrect number of fix precession/spin"); + + // set ptrs for fix langevin/spin styles + + // loop 1: obtain # of fix langevin/spin styles + for (iforce = 0; iforce < modify->nfix; iforce++) { if (strstr(modify->fix[iforce]->style,"langevin/spin")) { - maglangevin_flag = 1; - locklangevinspin = (FixLangevinSpin *) modify->fix[iforce]; + nlangspin++; } } + + // init length of vector of ptrs to precession/spin styles - if (maglangevin_flag) { - if (locklangevinspin->tdamp_flag == 1) tdamp_flag = 1; - if (locklangevinspin->temp_flag == 1) temp_flag = 1; + if (nlangspin > 0) { + locklangevinspin = new FixLangevinSpin*[nprecspin]; } + + // loop 2: fill vector with ptrs to precession/spin styles + count2 = 0; + if (nlangspin > 0) { + for (iforce = 0; iforce < modify->nfix; iforce++) { + if (strstr(modify->fix[iforce]->style,"langevin/spin")) { + maglangevin_flag = 1; + locklangevinspin[count2] = (FixLangevinSpin *) modify->fix[iforce]; + count2++; + } + } + } + + if (count2 != nlangspin) + error->all(FLERR,"Incorrect number of fix precession/spin"); + // ptrs FixSetForceSpin classes for (iforce = 0; iforce < modify->nfix; iforce++) { @@ -471,17 +517,16 @@ void FixNVESpin::ComputeInteractionsSpin(int i) // update magnetic precession interactions if (precession_spin_flag) { - lockprecessionspin->compute_single_precession(i,spi,fmi); + for (int k = 0; k < nprecspin; k++) { + lockprecessionspin[k]->compute_single_precession(i,spi,fmi); + } } // update langevin damping and random force if (maglangevin_flag) { // mag. langevin - if (tdamp_flag) { // transverse damping - locklangevinspin->add_tdamping(spi,fmi); - } - if (temp_flag) { // spin temperature - locklangevinspin->add_temperature(fmi); + for (int k = 0; k < nlangspin; k++) { + locklangevinspin[k]->compute_single_langevin(i,spi,fmi); } } @@ -496,7 +541,6 @@ void FixNVESpin::ComputeInteractionsSpin(int i) fm[i][0] = fmi[0]; fm[i][1] = fmi[1]; fm[i][2] = fmi[2]; - } /* ---------------------------------------------------------------------- diff --git a/src/SPIN/fix_nve_spin.h b/src/SPIN/fix_nve_spin.h index 5aa6b8e4e4..ac5bc57b25 100644 --- a/src/SPIN/fix_nve_spin.h +++ b/src/SPIN/fix_nve_spin.h @@ -61,11 +61,20 @@ friend class PairSpin; int tdamp_flag, temp_flag; int setforce_spin_flag; - // pointers to magnetic fixes + // pointers to fix langevin/spin styles - class FixPrecessionSpin *lockprecessionspin; - class FixLangevinSpin *locklangevinspin; - class FixSetForceSpin *locksetforcespin; + int nlangspin; + class FixLangevinSpin **locklangevinspin; + + // pointers to fix setforce/spin styles + + int nsetspin; + class FixSetForceSpin *locksetforcespin; // to be done + + // pointers to fix precession/spin styles + + int nprecspin; + class FixPrecessionSpin **lockprecessionspin; // pointers to magnetic pair styles diff --git a/src/SPIN/fix_precession_spin.cpp b/src/SPIN/fix_precession_spin.cpp index 17b9d3eb22..389ab1838d 100644 --- a/src/SPIN/fix_precession_spin.cpp +++ b/src/SPIN/fix_precession_spin.cpp @@ -67,6 +67,9 @@ FixPrecessionSpin::FixPrecessionSpin(LAMMPS *lmp, int narg, char **arg) : H_field = 0.0; nhx = nhy = nhz = 0.0; hx = hy = hz = 0.0; + stt_field = 0.0; + nsttx = nstty = nsttz = 0.0; + sttx = stty = sttz = 0.0; Ka = 0.0; nax = nay = naz = 0.0; Kax = Kay = Kaz = 0.0; @@ -74,8 +77,11 @@ FixPrecessionSpin::FixPrecessionSpin(LAMMPS *lmp, int narg, char **arg) : nc1x = nc1y = nc1z = 0.0; nc2x = nc2y = nc2z = 0.0; nc3x = nc3y = nc3z = 0.0; + K6 = 0.0; + n6x = n6y = n6z = 0.0; + m6x = m6y = m6z = 0.0; - zeeman_flag = aniso_flag = cubic_flag = 0; + zeeman_flag = stt_flag = aniso_flag = cubic_flag = hexaniso_flag = 0; int iarg = 3; while (iarg < narg) { @@ -87,6 +93,14 @@ FixPrecessionSpin::FixPrecessionSpin(LAMMPS *lmp, int narg, char **arg) : nhy = utils::numeric(FLERR,arg[iarg+3],false,lmp); nhz = utils::numeric(FLERR,arg[iarg+4],false,lmp); iarg += 5; + } else if (strcmp(arg[iarg],"stt") == 0) { + if (iarg+4 > narg) error->all(FLERR,"Illegal fix precession/spin command"); + stt_flag = 1; + stt_field = utils::numeric(FLERR,arg[iarg+1],false,lmp); + nsttx = utils::numeric(FLERR,arg[iarg+2],false,lmp); + nstty = utils::numeric(FLERR,arg[iarg+3],false,lmp); + nsttz = utils::numeric(FLERR,arg[iarg+4],false,lmp); + iarg += 5; } else if (strcmp(arg[iarg],"anisotropy") == 0) { if (iarg+4 > narg) error->all(FLERR,"Illegal fix precession/spin command"); aniso_flag = 1; @@ -110,40 +124,110 @@ FixPrecessionSpin::FixPrecessionSpin(LAMMPS *lmp, int narg, char **arg) : nc3y = utils::numeric(FLERR,arg[iarg+10],false,lmp); nc3z = utils::numeric(FLERR,arg[iarg+11],false,lmp); iarg += 12; + } else if (strcmp(arg[iarg],"hexaniso") == 0) { + if (iarg+7 > narg) error->all(FLERR,"Illegal fix precession/spin command"); + hexaniso_flag = 1; + K6 = utils::numeric(FLERR,arg[iarg+1],false,lmp); + n6x = utils::numeric(FLERR,arg[iarg+2],false,lmp); + n6y = utils::numeric(FLERR,arg[iarg+3],false,lmp); + n6z = utils::numeric(FLERR,arg[iarg+4],false,lmp); + m6x = utils::numeric(FLERR,arg[iarg+5],false,lmp); + m6y = utils::numeric(FLERR,arg[iarg+6],false,lmp); + m6z = utils::numeric(FLERR,arg[iarg+7],false,lmp); + iarg += 8; } else error->all(FLERR,"Illegal precession/spin command"); } // normalize vectors - double inorm; + double norm2,inorm; if (zeeman_flag) { - inorm = 1.0/sqrt(nhx*nhx + nhy*nhy + nhz*nhz); + norm2 = nhx*nhx + nhy*nhy + nhz*nhz; + if (norm2 == 0.0) + error->all(FLERR,"Illegal precession/spin command"); + inorm = 1.0/sqrt(norm2); nhx *= inorm; nhy *= inorm; nhz *= inorm; } + if (stt_flag) { + norm2 = nsttx*nsttx + nstty*nstty + nsttz*nsttz; + if (norm2 == 0.0) + error->all(FLERR,"Illegal precession/spin command"); + inorm = 1.0/sqrt(norm2); + nsttx *= inorm; + nstty *= inorm; + nsttz *= inorm; + } + if (aniso_flag) { - inorm = 1.0/sqrt(nax*nax + nay*nay + naz*naz); + norm2 = nax*nax + nay*nay + naz*naz; + if (norm2 == 0.0) + error->all(FLERR,"Illegal precession/spin command"); + inorm = 1.0/sqrt(norm2); nax *= inorm; nay *= inorm; naz *= inorm; } if (cubic_flag) { - inorm = 1.0/sqrt(nc1x*nc1x + nc1y*nc1y + nc1z*nc1z); + norm2 = nc1x*nc1x + nc1y*nc1y + nc1z*nc1z; + if (norm2 == 0.0) + error->all(FLERR,"Illegal precession/spin command"); + inorm = 1.0/sqrt(norm2); nc1x *= inorm; nc1y *= inorm; nc1z *= inorm; - inorm = 1.0/sqrt(nc2x*nc2x + nc2y*nc2y + nc2z*nc2z); + + norm2 = nc2x*nc2x + nc2y*nc2y + nc2z*nc2z; + if (norm2 == 0.0) + error->all(FLERR,"Illegal precession/spin command"); + inorm = 1.0/sqrt(norm2); nc2x *= inorm; nc2y *= inorm; nc2z *= inorm; - inorm = 1.0/sqrt(nc3x*nc3x + nc3y*nc3y + nc3z*nc3z); + + norm2 = nc3x*nc3x + nc3y*nc3y + nc3z*nc3z; + if (norm2 == 0.0) + error->all(FLERR,"Illegal precession/spin command"); + inorm = 1.0/sqrt(norm2); nc3x *= inorm; nc3y *= inorm; nc3z *= inorm; } + + if (hexaniso_flag) { + norm2 = n6x*n6x + n6y*n6y + n6z*n6z; + if (norm2 == 0.0) + error->all(FLERR,"Illegal precession/spin command"); + inorm = 1.0/sqrt(norm2); + n6x *= inorm; + n6y *= inorm; + n6z *= inorm; + + norm2 = m6x*m6x + m6y*m6y + m6z*m6z; + if (norm2 == 0.0) + error->all(FLERR,"Illegal precession/spin command"); + inorm = 1.0/sqrt(norm2); + m6x *= inorm; + m6y *= inorm; + m6z *= inorm; + l6x = (n6z*m6y-n6y*m6z); + l6y = (n6x*m6z-n6z*m6x); + l6z = (n6y*m6x-n6x*m6y); + + norm2 = l6x*l6x + l6y*l6y + l6z*l6z; + if (norm2 == 0.0) + error->all(FLERR,"Illegal precession/spin command"); + inorm = 1.0/sqrt(norm2); + l6x *= inorm; + l6y *= inorm; + l6z *= inorm; + m6x = (l6z*n6y-l6y*n6z); + m6y = (l6x*n6z-l6z*n6x); + m6z = (l6y*n6x-l6x*n6y); + } degree2rad = MY_PI/180.0; time_origin = update->ntimestep; @@ -185,6 +269,7 @@ void FixPrecessionSpin::init() Kah = Ka/hbar; k1ch = k1c/hbar; k2ch = k2c/hbar; + K6h = K6/hbar; if (utils::strmatch(update->integrate_style,"^respa")) { ilevel_respa = ((Respa *) update->integrate)->nlevels-1; @@ -199,15 +284,6 @@ 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; @@ -283,16 +359,26 @@ void FixPrecessionSpin::post_force(int /* vflag */) epreci -= compute_zeeman_energy(spi); } + if (stt_flag) { // compute Spin Transfer Torque + compute_stt(spi,fmi); + epreci -= compute_stt_energy(spi); + } + if (aniso_flag) { // compute magnetic anisotropy compute_anisotropy(spi,fmi); epreci -= compute_anisotropy_energy(spi); } - if (cubic_flag) { // compute cubic anisotropy + if (cubic_flag) { // compute cubic anisotropy compute_cubic(spi,fmi); epreci -= compute_cubic_energy(spi); } + if (hexaniso_flag) { // compute hexagonal anisotropy + compute_hexaniso(spi,fmi); + epreci -= compute_hexaniso_energy(spi); + } + emag[i] += epreci; eprec += epreci; fm[i][0] += fmi[0]; @@ -309,12 +395,16 @@ void FixPrecessionSpin::compute_single_precession(int i, double spi[3], double f int *mask = atom->mask; if (mask[i] & groupbit) { if (zeeman_flag) compute_zeeman(i,fmi); + if (stt_flag) compute_stt(spi,fmi); if (aniso_flag) compute_anisotropy(spi,fmi); if (cubic_flag) compute_cubic(spi,fmi); + if (hexaniso_flag) compute_hexaniso(spi,fmi); } } -/* ---------------------------------------------------------------------- */ +/* ---------------------------------------------------------------------- + Zeeman +------------------------------------------------------------------------- */ void FixPrecessionSpin::compute_zeeman(int i, double fmi[3]) { @@ -334,8 +424,32 @@ double FixPrecessionSpin::compute_zeeman_energy(double spi[4]) return energy; } +/* ---------------------------------------------------------------------- + STT +------------------------------------------------------------------------- */ + +void FixPrecessionSpin::compute_stt(double spi[3], double fmi[3]) +{ + double sx = spi[0]; + double sy = spi[1]; + double sz = spi[2]; + fmi[0] += 1.0*stt_field*( sy*nsttz-sz*nstty); + fmi[1] += 1.0*stt_field*(-sx*nsttz+sz*nsttx); + fmi[2] += 1.0*stt_field*( sx*nstty-sy*nsttx); +} + /* ---------------------------------------------------------------------- */ +double FixPrecessionSpin::compute_stt_energy(double spi[3]) +{ + double energy = 0.0; // Non-conservative force + return energy; +} + +/* ---------------------------------------------------------------------- + compute uniaxial anisotropy interaction for spin i +------------------------------------------------------------------------- */ + void FixPrecessionSpin::compute_anisotropy(double spi[3], double fmi[3]) { double scalar = nax*spi[0] + nay*spi[1] + naz*spi[2]; @@ -393,9 +507,7 @@ void FixPrecessionSpin::compute_cubic(double spi[3], double fmi[3]) fmi[2] += (fourz + sixz); } -/* ---------------------------------------------------------------------- - compute cubic aniso energy of spin i -------------------------------------------------------------------------- */ +/* ---------------------------------------------------------------------- */ double FixPrecessionSpin::compute_cubic_energy(double spi[3]) { @@ -412,6 +524,62 @@ double FixPrecessionSpin::compute_cubic_energy(double spi[3]) return energy; } +/* ---------------------------------------------------------------------- + compute hexagonal anisotropy interaction for spin i +------------------------------------------------------------------------- */ + +void FixPrecessionSpin::compute_hexaniso(double spi[3], double fmi[3]) +{ + double s_x,s_y,s_z; + double pf, phi, ssint2; + + // changing to the axes' frame + + s_x = l6x*spi[0]+l6y*spi[1]+l6z*spi[2]; + s_y = m6x*spi[0]+m6y*spi[1]+m6z*spi[2]; + s_z = n6x*spi[0]+n6y*spi[1]+n6z*spi[2]; + + // hexagonal anisotropy in the axes' frame + + phi = atan2(s_y,s_x); + ssint2 = s_x*s_x + s_y*s_y; // s^2sin^2(theta) + pf = 6.0 * K6h * ssint2*ssint2*sqrt(ssint2); // 6*K_6*s^5*sin^5(theta) + double fm_x = pf*cos(5*phi); + double fm_y = -pf*sin(5*phi); + double fm_z = 0; + + // back to the lab's frame + + fmi[0] += fm_x*l6x+fm_y*m6x+fm_z*n6x; + fmi[1] += fm_x*l6y+fm_y*m6y+fm_z*n6y; + fmi[2] += fm_x*l6z+fm_y*m6z+fm_z*n6z; +} + +/* ---------------------------------------------------------------------- + compute hexagonal aniso energy of spin i +------------------------------------------------------------------------- */ + +double FixPrecessionSpin::compute_hexaniso_energy(double spi[3]) +{ + double energy = 0.0; + double s_x,s_y,s_z, phi,ssint2; + + // changing to the axes' frame + + s_x = l6x*spi[0]+l6y*spi[1]+l6z*spi[2]; + s_y = m6x*spi[0]+m6y*spi[1]+m6z*spi[2]; + s_z = n6x*spi[0]+n6y*spi[1]+n6z*spi[2]; + + // hexagonal anisotropy in the axes' frame + + phi = atan2(s_y,s_z); + ssint2 = s_x*s_x + s_y*s_y; + + energy = K6 * ssint2*ssint2*ssint2*cos(6*phi); + + return 2.0*energy; +} + /* ---------------------------------------------------------------------- */ void FixPrecessionSpin::set_magneticprecession() @@ -421,6 +589,13 @@ void FixPrecessionSpin::set_magneticprecession() hy = H_field*nhy; hz = H_field*nhz; } + + if (stt_flag) { + sttx = stt_field*nsttx; + stty = stt_field*nstty; + sttz = stt_field*nsttz; + } + if (aniso_flag) { Kax = 2.0*Kah*nax; Kay = 2.0*Kah*nay; diff --git a/src/SPIN/fix_precession_spin.h b/src/SPIN/fix_precession_spin.h index 9c3c616077..d06a04f8f1 100644 --- a/src/SPIN/fix_precession_spin.h +++ b/src/SPIN/fix_precession_spin.h @@ -39,7 +39,7 @@ class FixPrecessionSpin : public Fix { void min_post_force(int); double compute_scalar(); - int zeeman_flag, aniso_flag, cubic_flag; + int zeeman_flag, stt_flag, aniso_flag, cubic_flag, hexaniso_flag; void compute_single_precession(int, double *, double *); // zeeman calculations @@ -47,6 +47,11 @@ class FixPrecessionSpin : public Fix { void compute_zeeman(int, double *); double compute_zeeman_energy(double *); + // stt calculations + + void compute_stt(double *, double *); + double compute_stt_energy(double *); + // uniaxial aniso calculations void compute_anisotropy(double *, double *); @@ -57,6 +62,11 @@ class FixPrecessionSpin : public Fix { void compute_cubic(double *, double *); double compute_cubic_energy(double *); + // hexagonal aniso calculations + + void compute_hexaniso(double *, double *); + double compute_hexaniso_energy(double *); + // storing magnetic energies int nlocal_max; // max nlocal (for list size) @@ -83,6 +93,12 @@ class FixPrecessionSpin : public Fix { double nhx, nhy, nhz; double hx, hy, hz; // temp. force variables + // STT intensity and direction + + double stt_field; + double nsttx, nstty, nsttz; + double sttx, stty, sttz; + // magnetic anisotropy intensity and direction double Ka; // aniso const. in eV @@ -98,6 +114,13 @@ class FixPrecessionSpin : public Fix { double nc2x,nc2y,nc2z; double nc3x,nc3y,nc3z; + // hexagonal anisotropy + double K6; // hexagonal aniso const. in eV + double K6h; // hexagonal aniso const. in rad.THz + double n6x,n6y,n6z; // main axis + double m6x,m6y,m6z; // secondary (perpendicular) axis + double l6x,l6y,l6z; // =(m x n) + void set_magneticprecession(); };