Commit JT 051319

- added a cubic anisotropy in fix_precession_spin
- added associated doc and examples
- corrected neb/spin commands in doc/src/
- added tools/spin/ description
This commit is contained in:
julient31 2019-05-13 16:59:39 -06:00
parent a7c9560dc1
commit 11f223416c
13 changed files with 392 additions and 55 deletions

View File

@ -83,7 +83,7 @@ An alphabetic list of all general LAMMPS commands.
"molecule"_molecule.html,
"ndx2group"_group2ndx.html,
"neb"_neb.html,
"neb_spin"_neb_spin.html,
"neb/spin"_neb_spin.html,
"neigh_modify"_neigh_modify.html,
"neighbor"_neighbor.html,
"newton"_newton.html,

Binary file not shown.

After

Width:  |  Height:  |  Size: 25 KiB

View File

@ -0,0 +1,21 @@
\documentclass[preview]{standalone}
\usepackage{varwidth}
\usepackage[utf8x]{inputenc}
\usepackage{amsmath,amssymb,amsthm,bm}
\begin{document}
\begin{varwidth}{50in}
\begin{equation}
\bm{H}_{cubic} = -\sum_{{ i}=1}^{N} K_{1}
\Big[
\left(\vec{s}_{i} \cdot \vec{n1} \right)^2
\left(\vec{s}_{i} \cdot \vec{n2} \right)^2 +
\left(\vec{s}_{i} \cdot \vec{n2} \right)^2
\left(\vec{s}_{i} \cdot \vec{n3} \right)^2 +
\left(\vec{s}_{i} \cdot \vec{n1} \right)^2
\left(\vec{s}_{i} \cdot \vec{n3} \right)^2 \Big]
+K_{2}^{(c)} \left(\vec{s}_{i} \cdot \vec{n1} \right)^2
\left(\vec{s}_{i} \cdot \vec{n2} \right)^2
\left(\vec{s}_{i} \cdot \vec{n3} \right)^2 \nonumber
\end{equation}
\end{varwidth}
\end{document}

View File

@ -77,6 +77,7 @@ Post-processing tools :h3
"python"_#pythontools,
"reax"_#reax_tool,
"smd"_#smd,
"spin"_#spin,
"xmgrace"_#xmgrace :tb(c=6,ea=c,a=l)
Miscellaneous tools :h3
@ -511,6 +512,20 @@ Ernst Mach Institute in Germany (georg.ganzenmueller at emi.fhg.de).
:line
spin tool :h4,link(spin)
The spin sub-directory contains a C file interpolate.c which can
be compiled and used to perform a cubic polynomial interpolation of
the MEP following a GNEB calculation.
See the README file in tools/spin/interpolate_gneb for more details.
This tool was written by the SPIN package author, Julien
Tranchida at Sandia National Labs (jtranch at sandia.gov, and by Aleksei
Ivanov, at University of Iceland (ali5 at hi.is).
:line
vim tool :h4,link(vim)
The files in the tools/vim directory are add-ons to the VIM editor

View File

@ -62,12 +62,12 @@ Commands :h1
mass
message
min_modify
min_spin
min/spin
min_style
minimize
molecule
neb
neb_spin
neb/spin
neigh_modify
neighbor
newton

View File

@ -14,19 +14,23 @@ fix ID group precession/spin style args :pre
ID, group are documented in "fix"_fix.html command :ulb,l
precession/spin = style name of this fix command :l
style = {zeeman} or {anisotropy} :l
style = {zeeman} or {anisotropy} or {cubic} :l
{zeeman} args = H x y z
H = intensity of the magnetic field (in Tesla)
x y z = vector direction of the field
{anisotropy} args = K x y z
K = intensity of the magnetic anisotropy (in eV)
x y z = vector direction of the anisotropy :pre
{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 :pre
:ule
[Examples:]
fix 1 all precession/spin zeeman 0.1 0.0 0.0 1.0
fix 1 all precession/spin anisotropy 0.001 0.0 0.0 1.0
fix 1 3 precession/spin anisotropy 0.001 0.0 0.0 1.0
fix 1 iron precession/spin cubic 0.001 0.0005 1.0 0.0 0.0 0.0 1.0 0.0 0.0 0.0 1.0
fix 1 all precession/spin zeeman 0.1 0.0 0.0 1.0 anisotropy 0.001 0.0 0.0 1.0 :pre
[Description:]
@ -50,10 +54,29 @@ for the magnetic spins in the defined group:
with n defining the direction of the anisotropy, and K (in eV) its intensity.
If K>0, an easy axis is defined, and if K<0, an easy plane is defined.
In both cases, the choice of (x y z) imposes the vector direction for the force.
Only the direction of the vector is important; it's length is ignored.
Style {cubic} is used to simulate a cubic anisotropy, with three
possible easy axis for the magnetic spins in the defined group:
Both styles can be combined within one single command line.
:c,image(Eqs/fix_spin_cubic.jpg)
with K1 and K2c (in eV) the intensity coefficients and
n1, n2 and n3 defining the three anisotropic directions
defined by the command (from n1x to n3z).
For n1 = (100), n2 = (010), and n3 = (001), K1 < 0 defines an
iron type anisotropy (easy axis along the (001)-type cube
edges), and K1 > 0 defines a nickel type anisotropy (easy axis
along the (111)-type cube diagonals).
K2^c > 0 also defines easy axis along the (111)-type cube
diagonals.
See chapter 2 of "(Skomski)"_#Skomski1 for more details on cubic
anisotropies.
In all cases, the choice of (x y z) only imposes the vector
directions for the forces. Only the direction of the vector is
important; it's length is ignored (the entered vectors are
normalized).
Those styles can be combined within one single command line.
:line
@ -85,3 +108,9 @@ package"_Build_package.html doc page for more info.
"atom_style spin"_atom_style.html
[Default:] none
:line
:link(Skomski1)
[(Skomski)] Skomski, R. (2008). Simple models of magnetism.
Oxford University Press.

View File

@ -6,7 +6,7 @@
:line
neb command :h3
neb/spin command :h3
[Syntax:]

View File

@ -0,0 +1,60 @@
# bcc iron in a 3d periodic box
clear
units metal
atom_style spin
dimension 3
boundary p p p
# necessary for the serial algorithm (sametag)
atom_modify map array
lattice bcc 2.8665
region box block 0.0 5.0 0.0 5.0 0.0 5.0
create_box 1 box
create_atoms 1 box
# setting mass, mag. moments, and interactions for bcc iron
mass 1 55.845
set group all spin 2.2 -1.0 0.0 0.0
velocity all create 100 4928459 rot yes dist gaussian
pair_style hybrid/overlay eam/alloy spin/exchange 3.5
pair_coeff * * eam/alloy Fe_Mishin2006.eam.alloy Fe
pair_coeff * * spin/exchange exchange 3.4 0.02726 0.2171 1.841
neighbor 0.1 bin
neigh_modify every 10 check yes delay 20
fix 1 all precession/spin cubic 0.001 0.0005 1.0 0.0 0.0 0.0 1.0 0.0 0.0 0.0 1.0
fix_modify 1 energy yes
fix 2 all langevin/spin 0.0 0.0 21
fix 3 all nve/spin lattice yes
timestep 0.0001
# compute and output options
compute out_mag all spin
compute out_pe all pe
compute out_ke all ke
compute out_temp all temp
variable magx equal c_out_mag[1]
variable magy equal c_out_mag[2]
variable magz equal c_out_mag[3]
variable magnorm equal c_out_mag[4]
variable emag equal c_out_mag[5]
variable tmag equal c_out_mag[6]
thermo_style custom step time v_magx v_magy v_magz v_magnorm v_tmag v_emag pe etotal
thermo 50
compute outsp all property/atom spx spy spz sp fmx fmy fmz
dump 100 all custom 1 dump_iron.lammpstrj type x y z c_outsp[1] c_outsp[2] c_outsp[3]
run 2000
# min_style spin
# min_modify alpha_damp 1.0 discrete_factor 10
# minimize 1.0e-16 1.0e-16 10000 10000

View File

@ -0,0 +1,60 @@
# fcc nickel in a 3d periodic box
clear
units metal
atom_style spin
dimension 3
boundary p p p
# necessary for the serial algorithm (sametag)
atom_modify map array
lattice fcc 3.524
region box block 0.0 5.0 0.0 5.0 0.0 5.0
create_box 1 box
create_atoms 1 box
# setting mass, mag. moments, and interactions for cobalt
mass 1 58.69
set group all spin/random 31 0.63
#set group all spin 0.63 0.0 0.0 1.0
velocity all create 100 4928459 rot yes dist gaussian
pair_style hybrid/overlay eam/alloy spin/exchange 4.0
pair_coeff * * eam/alloy Ni99.eam.alloy Ni
pair_coeff * * spin/exchange exchange 4.0 0.50 0.2280246862 1.229983475
neighbor 0.1 bin
neigh_modify every 10 check yes delay 20
fix 1 all precession/spin cubic -0.0001 0.0 1.0 0.0 0.0 0.0 1.0 0.0 0.0 0.0 1.0 &
zeeman 0.0 0.0 0.0 1.0
fix_modify 1 energy yes
fix 2 all langevin/spin 0.0 0.0 21
fix 3 all nve/spin lattice yes
timestep 0.0001
# compute and output options
compute out_mag all spin
compute out_pe all pe
compute out_ke all ke
compute out_temp all temp
variable magz equal c_out_mag[3]
variable magnorm equal c_out_mag[4]
variable emag equal c_out_mag[5]
variable tmag equal c_out_mag[6]
thermo_style custom step time v_magnorm v_emag temp v_tmag etotal
thermo 50
compute outsp all property/atom spx spy spz sp fmx fmy fmz
dump 50 all custom 1 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]
run 2000

2
src/.gitignore vendored
View File

@ -156,6 +156,8 @@
/fix_nve_spin.h
/fix_precession_spin.cpp
/fix_precession_spin.h
/fix_setforce_spin.cpp
/fix_setforce_spin.h
/min_spin.cpp
/min_spin.h
/neb_spin.cpp

View File

@ -25,8 +25,8 @@
#include <cstdio>
#include <cstdlib>
#include <cstring>
#include "atom.h"
#include "error.h"
#include "domain.h"
#include "error.h"
#include "fix_precession_spin.h"
@ -35,6 +35,7 @@
#include "math_const.h"
#include "memory.h"
#include "modify.h"
#include "neigh_list.h"
#include "respa.h"
#include "update.h"
#include "variable.h"
@ -71,8 +72,12 @@ FixPrecessionSpin::FixPrecessionSpin(LAMMPS *lmp, int narg, char **arg) : Fix(lm
Ka = 0.0;
nax = nay = naz = 0.0;
Kax = Kay = Kaz = 0.0;
k1c = k2c = 0.0;
nc1x = nc1y = nc1z = 0.0;
nc2x = nc2y = nc2z = 0.0;
nc3x = nc3y = nc3z = 0.0;
zeeman_flag = aniso_flag = 0;
zeeman_flag = aniso_flag = cubic_flag = 0;
int iarg = 3;
while (iarg < narg) {
@ -92,14 +97,61 @@ FixPrecessionSpin::FixPrecessionSpin(LAMMPS *lmp, int narg, char **arg) : Fix(lm
nay = force->numeric(FLERR,arg[iarg+3]);
naz = force->numeric(FLERR,arg[iarg+4]);
iarg += 5;
} else if (strcmp(arg[iarg],"cubic") == 0) {
if (iarg+2 > narg) error->all(FLERR,"Illegal fix precession/spin command");
cubic_flag = 1;
k1c = force->numeric(FLERR,arg[iarg+1]);
k2c = force->numeric(FLERR,arg[iarg+2]);
nc1x = force->numeric(FLERR,arg[iarg+3]);
nc1y = force->numeric(FLERR,arg[iarg+4]);
nc1z = force->numeric(FLERR,arg[iarg+5]);
nc2x = force->numeric(FLERR,arg[iarg+6]);
nc2y = force->numeric(FLERR,arg[iarg+7]);
nc2z = force->numeric(FLERR,arg[iarg+8]);
nc3x = force->numeric(FLERR,arg[iarg+9]);
nc3y = force->numeric(FLERR,arg[iarg+10]);
nc3z = force->numeric(FLERR,arg[iarg+11]);
iarg += 12;
} else error->all(FLERR,"Illegal precession/spin command");
}
// normalize vectors
double inorm;
if (zeeman_flag) {
inorm = 1.0/sqrt(nhx*nhx + nhy*nhy + nhz*nhz);
nhx *= inorm;
nhy *= inorm;
nhz *= inorm;
}
if (aniso_flag) {
inorm = 1.0/sqrt(nax*nax + nay*nay + naz*naz);
nax *= inorm;
nay *= inorm;
naz *= inorm;
}
if (cubic_flag) {
inorm = 1.0/sqrt(nc1x*nc1x + nc1y*nc1y + nc1z*nc1z);
nc1x *= inorm;
nc1y *= inorm;
nc1z *= inorm;
inorm = 1.0/sqrt(nc2x*nc2x + nc2y*nc2y + nc2z*nc2z);
nc2x *= inorm;
nc2y *= inorm;
nc2z *= inorm;
inorm = 1.0/sqrt(nc3x*nc3x + nc3y*nc3y + nc3z*nc3z);
nc3x *= inorm;
nc3y *= inorm;
nc3z *= inorm;
}
degree2rad = MY_PI/180.0;
time_origin = update->ntimestep;
eflag = 0;
emag = 0.0;
eprec = 0.0;
}
/* ---------------------------------------------------------------------- */
@ -130,8 +182,12 @@ void FixPrecessionSpin::init()
const double mub = 5.78901e-5; // in eV/T
const double gyro = mub/hbar; // in rad.THz/T
H_field *= gyro; // in rad.THz
Ka /= hbar; // in rad.THz
// convert field quantities to rad.THz
H_field *= gyro;
Kah = Ka/hbar;
k1ch = k1c/hbar;
k2ch = k2c/hbar;
if (strstr(update->integrate_style,"respa")) {
ilevel_respa = ((Respa *) update->integrate)->nlevels-1;
@ -185,53 +241,60 @@ void FixPrecessionSpin::post_force(int /* vflag */)
if (varflag != CONSTANT) {
modify->clearstep_compute();
modify->addstep_compute(update->ntimestep + 1);
set_magneticprecession(); // update mag. field if time-dep.
set_magneticprecession(); // update mag. field if time-dep.
}
double **sp = atom->sp;
int *mask = atom->mask;
double **fm = atom->fm;
double spi[3], fmi[3];
double **sp = atom->sp;
const int nlocal = atom->nlocal;
double spi[3], fmi[3], epreci;
eflag = 0;
emag = 0.0;
eprec = 0.0;
for (int i = 0; i < nlocal; i++) {
spi[0] = sp[i][0];
spi[1] = sp[i][1];
spi[2] = sp[i][2];
fmi[0] = fmi[1] = fmi[2] = 0.0;
if (mask[i] & groupbit) {
epreci = 0.0;
spi[0] = sp[i][0];
spi[1] = sp[i][1];
spi[2] = sp[i][2];
fmi[0] = fmi[1] = fmi[2] = 0.0;
if (zeeman_flag) { // compute Zeeman interaction
compute_zeeman(i,fmi);
emag -= (spi[0]*fmi[0] + spi[1]*fmi[1] + spi[2]*fmi[2]);
if (zeeman_flag) { // compute Zeeman interaction
compute_zeeman(i,fmi);
epreci -= hbar*(spi[0]*fmi[0] + spi[1]*fmi[1] + spi[2]*fmi[2]);
}
if (aniso_flag) { // compute magnetic anisotropy
compute_anisotropy(spi,fmi);
epreci -= compute_anisotropy_energy(spi);
}
if (cubic_flag) { // compute cubic anisotropy
compute_cubic(spi,fmi);
epreci -= compute_cubic_energy(spi);
}
eprec += epreci;
fm[i][0] += fmi[0];
fm[i][1] += fmi[1];
fm[i][2] += fmi[2];
}
if (aniso_flag) { // compute magnetic anisotropy
compute_anisotropy(spi,fmi);
emag -= 0.5*(spi[0]*fmi[0] + spi[1]*fmi[1] + spi[2]*fmi[2]);
}
fm[i][0] += fmi[0];
fm[i][1] += fmi[1];
fm[i][2] += fmi[2];
}
emag *= hbar;
}
/* ---------------------------------------------------------------------- */
void FixPrecessionSpin::compute_single_precession(int i, double spi[3], double fmi[3])
{
if (zeeman_flag) {
compute_zeeman(i,fmi);
}
if (aniso_flag) {
compute_anisotropy(spi,fmi);
int *mask = atom->mask;
if (mask[i] & groupbit) {
if (zeeman_flag) compute_zeeman(i,fmi);
if (aniso_flag) compute_anisotropy(spi,fmi);
if (cubic_flag) compute_cubic(spi,fmi);
}
}
/* ---------------------------------------------------------------------- */
void FixPrecessionSpin::compute_zeeman(int i, double fmi[3])
@ -254,6 +317,16 @@ void FixPrecessionSpin::compute_anisotropy(double spi[3], double fmi[3])
/* ---------------------------------------------------------------------- */
double FixPrecessionSpin::compute_anisotropy_energy(double spi[3])
{
double energy = 0.0;
double scalar = nax*spi[0] + nay*spi[1] + naz*spi[2];
energy = Ka*scalar*scalar;
return energy;
}
/* ---------------------------------------------------------------------- */
void FixPrecessionSpin::post_force_respa(int vflag, int ilevel, int /*iloop*/)
{
if (ilevel == ilevel_respa) post_force(vflag);
@ -264,17 +337,75 @@ void FixPrecessionSpin::post_force_respa(int vflag, int ilevel, int /*iloop*/)
void FixPrecessionSpin::set_magneticprecession()
{
if (zeeman_flag) {
hx = H_field*nhx;
hy = H_field*nhy;
hz = H_field*nhz;
hx = H_field*nhx;
hy = H_field*nhy;
hz = H_field*nhz;
}
if (aniso_flag) {
Kax = 2.0*Ka*nax;
Kay = 2.0*Ka*nay;
Kaz = 2.0*Ka*naz;
Kax = 2.0*Kah*nax;
Kay = 2.0*Kah*nay;
Kaz = 2.0*Kah*naz;
}
}
/* ----------------------------------------------------------------------
compute cubic aniso energy of spin i
------------------------------------------------------------------------- */
double FixPrecessionSpin::compute_cubic_energy(double spi[3])
{
double energy = 0.0;
double skx,sky,skz;
skx = spi[0]*nc1x+spi[1]*nc1y+spi[2]*nc1z;
sky = spi[0]*nc2x+spi[1]*nc2y+spi[2]*nc2z;
skz = spi[0]*nc3x+spi[1]*nc3y+spi[2]*nc3z;
energy = k1c*(skx*skx*sky*sky + sky*sky*skz*skz + skx*skx*skz*skz);
energy += k2c*skx*skx*sky*sky*skz*skz;
return energy;
}
/* ----------------------------------------------------------------------
compute cubic anisotropy interaction for spin i
------------------------------------------------------------------------- */
void FixPrecessionSpin::compute_cubic(double spi[3], double fmi[3])
{
double skx,sky,skz,skx2,sky2,skz2;
double four1,four2,four3,fourx,foury,fourz;
double six1,six2,six3,sixx,sixy,sixz;
skx = spi[0]*nc1x+spi[1]*nc1y+spi[2]*nc1z;
sky = spi[0]*nc2x+spi[1]*nc2y+spi[2]*nc2z;
skz = spi[0]*nc3x+spi[1]*nc3y+spi[2]*nc3z;
skx2 = skx*skx;
sky2 = sky*sky;
skz2 = skz*skz;
four1 = 2.0*skx*(sky2+skz2);
four2 = 2.0*sky*(skx2+skz2);
four3 = 2.0*skz*(skx2+sky2);
fourx = k1ch*(nc1x*four1 + nc2x*four2 + nc3x*four3);
foury = k1ch*(nc1y*four1 + nc2y*four2 + nc3y*four3);
fourz = k1ch*(nc1z*four1 + nc2z*four2 + nc3z*four3);
six1 = 2.0*skx*sky2*skz2;
six2 = 2.0*sky*skx2*skz2;
six3 = 2.0*skz*skx2*sky2;
sixx = k2ch*(nc1x*six1 + nc2x*six2 + nc3x*six3);
sixy = k2ch*(nc1y*six1 + nc2y*six2 + nc3y*six3);
sixz = k2ch*(nc1z*six1 + nc2z*six2 + nc3z*six3);
fmi[0] += fourx + sixx;
fmi[1] += foury + sixy;
fmi[2] += fourz + sixz;
}
/* ----------------------------------------------------------------------
potential energy in magnetic field
------------------------------------------------------------------------- */
@ -284,10 +415,10 @@ double FixPrecessionSpin::compute_scalar()
// only sum across procs one time
if (eflag == 0) {
MPI_Allreduce(&emag,&emag_all,1,MPI_DOUBLE,MPI_SUM,world);
MPI_Allreduce(&eprec,&eprec_all,1,MPI_DOUBLE,MPI_SUM,world);
eflag = 1;
}
return emag_all;
return eprec_all;
}
/* ---------------------------------------------------------------------- */

View File

@ -39,11 +39,20 @@ class FixPrecessionSpin : public Fix {
void min_post_force(int);
double compute_scalar();
int zeeman_flag, aniso_flag;
int zeeman_flag, aniso_flag, cubic_flag;
void compute_single_precession(int, double *, double *);
void compute_zeeman(int, double *);
void compute_anisotropy(double *, double *);
// uniaxial aniso calculations
void compute_anisotropy(double *, double *);
double compute_anisotropy_energy(double *);
// cubic aniso calculations
void compute_cubic(double *, double *);
double compute_cubic_energy(double *);
protected:
int style; // style of the magnetic precession
@ -52,7 +61,7 @@ class FixPrecessionSpin : public Fix {
int ilevel_respa;
int time_origin;
int eflag;
double emag, emag_all;
double eprec, eprec_all;
int varflag;
int magfieldstyle;
@ -67,10 +76,19 @@ class FixPrecessionSpin : public Fix {
// magnetic anisotropy intensity and direction
double Ka;
double Ka; // aniso const. in eV
double Kah; // aniso const. in rad.THz
double nax, nay, naz;
double Kax, Kay, Kaz; // temp. force variables
// cubic anisotropy intensity
double k1c,k2c; // cubic const. in eV
double k1ch,k2ch; // cubic const. in rad.THz
double nc1x,nc1y,nc1z;
double nc2x,nc2y,nc2z;
double nc3x,nc3y,nc3z;
void set_magneticprecession();
};

View File

@ -39,6 +39,7 @@ pymol_asphere convert LAMMPS output of ellipsoids to PyMol format
python Python scripts for post-processing LAMMPS output
reax Tools for analyzing output of ReaxFF simulations
smd convert Smooth Mach Dynamics triangles to VTK
spin perform a cubic polynomial interpolation of a GNEB MEP
vim add-ons to VIM editor for editing LAMMPS input scripts
xmgrace a collection of scripts to generate xmgrace plots