ICODE-ADC/lammps
4
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

Merge pull request #2714 from julient31/STT2

Browse Source 
Improvement of the SPIN package by removing a limitation
This commit is contained in:
Axel Kohlmeyer
2021-04-21 14:55:17 -04:00
committed by GitHub
parent 2e01bb8f62 6cad2ba829
commit a64a7d5096
11 changed files with 389 additions and 73 deletions
Download Patch File Download Diff File Expand all files Collapse all files

34
doc/src/fix_precession_spin.rst
View File

@ -12,7 +12,7 @@ Syntax
* ID, group are documented in :doc:`fix <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) <Skomski1>` 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) <Chirac1>` 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
<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.

11
examples/SPIN/test_problems/README
View File

@ -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.

25
examples/SPIN/test_problems/run_all.sh Executable file
View File

@ -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 ..

4
examples/SPIN/test_problems/validation_nvt/in.spin.nvt_lattice
View File

@ -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

9
examples/SPIN/test_problems/validation_nvt/in.spin.nvt_spin
View File

@ -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

20
src/SPIN/fix_langevin_spin.cpp
View File

@ -24,9 +24,11 @@
#include "fix_langevin_spin.h"
#include <cmath>
#include <cstring>
#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);
}
}

8
src/SPIN/fix_langevin_spin.h
View File

@ -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;

90
src/SPIN/fix_nve_spin.cpp
View File

@ -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)
{
@ -75,6 +76,9 @@ FixNVESpin::FixNVESpin(LAMMPS *lmp, int narg, char **arg) :
npairs = 0;
npairspin = 0;
// test nprec
nprecspin = nlangspin = nsetspin = 0;
// checking if map array or hash is defined
if (atom->map_style == Atom::MAP_NONE)
@ -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; i<npairs; i++) {
if (force->pair_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++;
}
}
// ptrs on the FixLangevinSpin class
// init length of vector of ptrs to precession/spin styles
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++;
}
}
if (maglangevin_flag) {
if (locklangevinspin->tdamp_flag == 1) tdamp_flag = 1;
if (locklangevinspin->temp_flag == 1) temp_flag = 1;
// init length of vector of ptrs to precession/spin styles
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];
}
/* ----------------------------------------------------------------------

17
src/SPIN/fix_nve_spin.h
View File

@ -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

217
src/SPIN/fix_precession_spin.cpp
View File

@ -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,41 +124,111 @@ 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;

25
src/SPIN/fix_precession_spin.h
View File

@ -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();
};