Part 2: Energy differences comparing collinear magnetic structures ¶
By the end of this tutorial, you will be able to:
- determine nearest-neighbor and next-nearest-neighbor coupling for a Heisenberg model for NiO
- scale the energy and k-mesh with the size of the unit cell
5.1 Task¶
Compute the exchange-interaction parameters for NiO for a Heisenberg model including nearest and next-nearest neighbors based on a LDA+U calculation with U=8eV, J=0.95eV.
A Heisenberg model for NiO that accounts for nearest neighbor $\langle i,j \rangle$ and next nearest neighbor $\langle\langle i,j \rangle\rangle$ interactions $J_1$ and $J_2$ can be written as
$$ H = -J_1 \sum_{\langle i,j \rangle} \vec{e}_i \cdot \vec{e}_j - J_2 \sum_{\langle\langle i,j \rangle\rangle} \vec{e}_i \cdot \vec{e}_j, $$
where $\vec{e}_i, \vec{e}_j$ are unit vectors along the direction of the on-site magnetic moment. We want to consider a ferromagentic (FM) order, and two antiferromagnetic (AFM) orders, one that orders the $(001)/(110)$ planes (AF1), and a second (AF2) that orders the $(111)$ planes. The energies per formula unit in terms of the Heisenberg model are
$$ E_{FM} = E_0 - 6J_1 - 3 J_2, $$ $$ E_{AF1} = E_0 + 2J_1 - 3 J_2, $$ $$ E_{AF2} = E_0 + 3 J_2. $$
This set of equations can be used to express the exchange-interaction parameters $J_1$ and $J_2$ in terms of the total energies of the FM, AF1 and AF2 orders:
$$ J_1 = \frac{1}{8} \left(E_{AF1}-E_{FM}\right), $$ $$ J_2 = \frac{1}{24} \left(4E_{AF2}-3E_{AF1}-E_{FM}\right). $$
Using a collinear magnetic ab-initio calculation, we can obtain these total energies and thus provide model parameters. This is particularly useful if the model Hamiltonian can be solved to access other properties like the critical temperature.
A detailed discussion can be found in Archer et. al., PRB 84, 115114 (2011).
5.2 Input¶
The input files to run this example are prepared at $TUTORIALS/magnetism/e05_NiO_Heisenberg. Check them out!
AFM2/POSCAR
AFM NiO
4.11908875913573
1.0000000000000000 0.5000000000000000 0.5000000000000000
0.5000000000000000 1.0000000000000000 0.5000000000000000
0.5000000000000000 0.5000000000000000 1.0000000000000000
Ni O
2 2
Direct
0.0000000000000000 0.0000000000000000 0.0000000000000000
0.5000000000000000 0.5000000000000000 0.5000000000000000
0.2500000000000000 0.2500000000000000 0.2500000000000000
0.7500000000000000 0.7500000000000000 0.7500000000000000
NiO
4.11908875913573
1.0000000000000000 0.0000000000000000 0.0000000000000000
0.0000000000000000 1.0000000000000000 0.0000000000000000
0.0000000000000000 0.0000000000000000 1.0000000000000000
Ni O
4 4
Direct
0.0000000000000000 0.0000000000000000 0.0000000000000000
0.5000000000000000 0.5000000000000000 0.0000000000000000
0.5000000000000000 0.0000000000000000 0.5000000000000000
0.0000000000000000 0.5000000000000000 0.5000000000000000
0.5000000000000000 0.0000000000000000 0.0000000000000000
0.0000000000000000 0.5000000000000000 0.0000000000000000
0.0000000000000000 0.0000000000000000 0.5000000000000000
0.5000000000000000 0.5000000000000000 0.5000000000000000
AFM2/INCAR
SYSTEM = NiO
ISTART = 0
# Magnetism AFM2
ISPIN = 2
MAGMOM = 5.0 -5.0 2*0
LASPH = T
# FFT grid
ENCUT = 450
PREC = Acc
# Convergence criterion
EDIFF = 1E-7
# Smearing/k integration
ISMEAR = -5
EFERMI = MIDGAP
# Algorithm and density mixing
ALGO = Normal
IMIX = 4
AMIX = 0.8
BMIX = 0.00001
AMIX_MAG = 0.8
BMIX_MAG = 0.00001
LMAXMIX = 4
# compute bandgap
BANDGAP = WEIGHT
# compute on-site moment
LORBIT = 11
# Exchange-correlation functional
XC = CA
LDAU = .TRUE.
LDAUTYPE = 2
LDAUL = 2 -1
LDAUU = 8.00 0.00
LDAUJ = 0.95 0.00
LDAUPRINT = 1
AFM1/INCAR
SYSTEM = NiO
ISTART = 0
# Magnetism AFM1
ISPIN = 2
MAGMOM = 2*-5.0 2*5.0 4*0
LASPH = T
# FFT grid
ENCUT = 450
PREC = Acc
# Convergence criterion
EDIFF = 1E-7
# Smearing/k integration
ISMEAR = -5
EFERMI = MIDGAP
# Algorithm and density mixing
ALGO = Normal
IMIX = 4
AMIX = 0.8
BMIX = 0.00001
AMIX_MAG = 0.8
BMIX_MAG = 0.00001
LMAXMIX = 4
# compute bandgap
BANDGAP = WEIGHT
# compute on-site moment
LORBIT = 11
# Exchange-correlation functional
XC = CA
LDAU = .TRUE.
LDAUTYPE = 2
LDAUL = 2 -1
LDAUU = 8.00 0.00
LDAUJ = 0.95 0.00
LDAUPRINT = 1
FM/INCAR
SYSTEM = NiO
ISTART = 0
# Magnetism FM
ISPIN = 2
MAGMOM = 4*5.0 4*0
LASPH = T
# FFT grid
ENCUT = 450
PREC = Acc
# Convergence criterion
EDIFF = 1E-7
# Smearing/k integration
ISMEAR = -5
EFERMI = MIDGAP
# Algorithm and density mixing
ALGO = Normal
IMIX = 4
AMIX = 0.8
BMIX = 0.00001
AMIX_MAG = 0.8
BMIX_MAG = 0.00001
LMAXMIX = 4
# compute bandgap
BANDGAP = WEIGHT
# compute on-site moment
LORBIT = 11
# Exchange-correlation functional
XC = CA
LDAU = .TRUE.
LDAUTYPE = 2
LDAUL = 2 -1
LDAUU = 8.00 0.00
LDAUJ = 0.95 0.00
LDAUPRINT = 1
AFM2/KPOINTS
Regular k-point mesh
0
Gamma
5 5 5
Regular k-point mesh
0
Gamma
4 4 4
Pseudopotential: Ni_sv_GW and O_s.
The INCAR files set the basic settings for a DFT+U calculation. The cutoff energies and k-mesh density are not well converged here to keep the tutorial computationally cheap. The usual convergence tests, should be done for a real application. In order to obtain precise total energies, we recommend using the tetrahedron method by setting ISMEAR = -5 and the convergence criterion EDIFF should be set relatively tight. The difference between the INCAR files is the initial guess for the magnetic configuration set by MAGMOM. MAGMOM is used to determine the symmetry operations for the final result, so it may not be surprising that the calculations in AFM1 and FM subdirectories converge to different final magnetic configurations. However, in practice, even without symmetry constraints the final magnetic configuration is often close to the initial magnetic configuration. Here, this is an advantage, but generally that means the magnetic ground state is very hard to determine. Also note that
MAGMOM = 4*5 4*0
is short hand for
MAGMOM = 5 5 5 5 0 0 0 0
The number of subdivisions to define the k mesh in the KPOINTS file is chosen inversely proportional to the length of the lattice vectors. This results in a comparible k mesh density. To obtain accurate total energies, it is more important to converge each calculation with respect to k-mesh density than to use a comparable k mesh. Consider an extreme example: if one magnetic structure is metallic and the other gapped, long range interactions and thus a finer k mesh are expected to be more relevant for the metallic system.
5.3 Calculation¶
Open a terminal and navigate to this example's directory and run it by entering the following:
cd $TUTORIALS/magnetism/e05_*
for d in FM AFM1 AFM2
do
cd $d
cp ../POTCAR .
mpirun -np 2 vasp_std
cd ..
done
Now, use the total energy to compute the exchange-interaction parameters:
import py4vasp
EFM = py4vasp.Calculation.from_path("e05_NiO_Heisenberg/FM").energy.read()["free energy TOTEN"]/4
EAFM1 = py4vasp.Calculation.from_path("e05_NiO_Heisenberg/AFM1").energy.read()["free energy TOTEN"]/4
EAFM2 = py4vasp.Calculation.from_path("e05_NiO_Heisenberg/AFM2").energy.read()["free energy TOTEN"]/2
J1 = (EAFM1-EFM)/8
J2 = (4*EAFM2-3*EAFM1-EFM)/24
print("J1 in meV: ",J1*1000)
print("J2 in meV: ",J2*1000)
Mind that the unit cells of the magnetic configurations differ. Thus the total energy read from file are scaled accordingly before they can be used to compute the exchange-coupling parameters.
5.4 Questions¶
- Is MAGMOM = 2*3 equivalent to MAGMOM = 6?
- Which setting for ISMEAR is best for accurate total energies?
- Is it important to use comparable k-mesh density, when comparing total energies?
By the end of this tutorial, you will be able to:
- switch on spin-orbit coupling
- decide when to use vasp_ncl instead of vasp_std
6.1 Task¶
Compute the band structure of BCC Fe inluding spin-orbit coupling along the high-symmetry path H$_3$-Γ-H.
In order to include spin-orbit coupling (SOC), we need to employ the framework of spin-density functional theory as is used to describe noncollinear magnetism, see Hobbs et al. PRB 62, 11556 (2000). Hence, the system is described in terms of a $2\times2$ spin-density matrix with entries
$$ n(\vec{r}) = \begin{pmatrix} n_{\uparrow \uparrow}(\vec{r}) & n_{\uparrow \downarrow}(\vec{r}) \\ n_{\downarrow \uparrow}(\vec{r}) & n_{\downarrow \downarrow}(\vec{r}) \end{pmatrix}. $$
The SOC term in the zeroth-order-regular approximation (ZORA) reads
$$ \hat{H}_{SO}= \frac{1}{(2c)^2} \frac{K^2(r)}{r} \frac{\text{d}V(r)}{\text{d}r} \vec{\sigma} \cdot \vec{L}, $$ $$ K(r) = \frac{1}{ 1 - \frac{V(r)}{2c^2}} $$
in Gaussian CGS units, see Lenthe et al. J Chem Phys 99, 4597 (1993), Steiner et al. PRB 93, 224425 (2016) and Eq. (19) in Speelman et al. arXiv:2504.17380.
The energy scale of SOC is quite small, but it can have substantial consequences. Here, for instance, the high-symmetry path Γ-H$_3$ and Γ-H in BCC Fe would be equivalent without SOC, while the calculation including SOC shows so-called anticrossings that are typical for band structures including SOC.
6.2 Input¶
The input files to run this example are prepared at $TUTORIALS/magnetism/e06_bcc_Fe_SOC. Check them out!
BCC Fe
1.4315177495
-1.0000000000000000 1.0000000000000000 1.0000000000000000
1.0000000000000005 -1.0000000000000000 1.0000000000000000
1.0000000000000000 1.0000000000000000 -1.0000000000000000
Fe
1
direct
0.0000000000000000 0.0000000000000000 0.0000000000000000 Fe
ENCUT = 300
EDIFF = 1e-8
ALGO = Normal
ISMEAR = 0
SIGMA = 0.05
LNONCOLLINEAR = T
LSORBIT = T
MAGMOM = 0.0 0.0 2.0
LASPH = T
LMAXMIX = 4
AMIX = 0.2
BMIX = 0.0001 ! almost zero, but 0 will crash some versions
AMIX_MAG = 0.8
BMIX_MAG = 0.0001 ! almost zero, but 0 will crash some versions
LORBIT = 11
0
Gamma
12 12 12
k points for band structure
141 ! intersections
line mode
reciprocal
0.5000000000 0.5000000000 -0.5000000000 H3
0.0000000000 0.0000000000 0.0000000000 Γ
0.0000000000 0.0000000000 0.0000000000 Γ
0.5000000000 -0.5000000000 0.5000000000 H
Pseudopotential: Fe.
Check out the tags LNONCOLLINEAR, LSORBIT and note how the format of MAGMOM has changed in accordance!
A uniform k mesh is provided in the KPOINTS file and the high-symmetry path is set in the KPOINTS_OPT file. Hence, we will perform a single calculation without restarting to compute the band structure.
6.3 Calculation¶
Open a terminal, navigate to this example's directory and run VASP by entering the following:
cd $TUTORIALS/magnetism/e06_*
mpirun -np 2 vasp_ncl
Note that here we use the executable vasp_ncl. ncl is an abbreviation of noncollinear and refers to the fact that the full $2\times2$ spin-density matrix is treated as required for the description of noncollinear magnetism, but also for SOC in general.
Plot the band structure using py4vasp! Mind that the Fermi energy needs to be computed on the uniform k mesh. Thus, the value needs to be provided as an argument, when the eigenvalues along a k path should be interpreted. This is done in the following Python code:
import py4vasp as pv
# Read calculation
calc = pv.Calculation.from_path("e06_bcc_Fe_SOC")
# Fermi energy
ef = calc.band.to_dict()['fermi_energy']
# Plot band structure
calc.band.plot("kpoints_opt", fermi_energy=ef)
Zoom in around the Fermi energy to see the anti-crossings. Compare the path between Γ-H$_3$ and Γ-H! The anti-crossings appear at different point. Hence, the crystallographically equivalent directions are now different in terms of their magnetic spacegroup symmetry. The inclusion of SOC introduces a coupling of spin degrees of freedom and lattice degrees of freedom, so it is not enough to perform a collinear spin-polarized calculation to obtain a difference between Γ-H$_3$ and Γ-H.
6.4 Questions¶
- Can spin-orbit coupling be included using vasp_std?
- What are anti-crossings in the band structure?
By the end of this tutorial, you will be able to:
- define the orientation of spinor space with respect to Cartesian space
- restart a noncollinear calculation from a collinear calculation
- compute the magnetic anisotropy energy of a collinear magnetic structure
7.1 Task¶
Compute the magnetic anisotropy energy of bulk FeO using LDA+U with U = 4eV for two angles $\phi=0$ and $\phi=90$ scanning different values of $\theta$, where $\theta=0$ corresponds to the cubic [1,1,1] direction and $\phi=0$ lies within the cubic $(1,\bar{1},0)$ plane.
Magnetic anisotropy is a relativistic effect that arises due to spin-orbit coupling (SOC). FeO has a rock-salt-like crystal structure, which is slightly distorted such that it becomes monoclinic instead. For the ground-state magnetic structure the $(1,1,1)$ planes are AFM ordered. The magnetic anisotropy energy (MAE) for MnO, FeO, CoO, and NiO have been discussed in detail by Schrön, Rödl, and Bechstedt, PRB 86, 115134 (2012). Here, we reproduce the result of Fig. 5 (a).
The workflow for a MAE calculation is
- SCF calculation using vasp_std with ISPIN=2, ISYM=-1.
- Restart noncollinear calculation with fixed collinear charge density (ICHARG=11) including SOC and rotate SAXIS along the MAE path.
This approach has the advantage that the direction of the magnetization can be varied without worrying that the calculation will converge to a different direction of the magnetization. But it relies on the fact that SOC is a small perturbation so the charge density and the magnetization do not change significantly due to SOC. This workflow presumes that the on-site orbital angular moment orders based on the crystal field instead of following the direction of the on-site magnetic moment.
7.2 Input¶
The resources to create the input and run this example are prepared at $TUTORIALS/magnetism/e07_FeO_bulk_MAE. We will discuss how to setup the calculation in detail.
The path of the MAE calculation can be characterized in terms of the polar angles $\theta$ and $\phi$. We chose the angles the same way as Schrön et al.:
- $\theta=0$ if the magnetization is parallel to the cubic [111] direction
- $\phi=0$ if the magnetization lies in the $(1\bar{1}0)$ plane.
It is a bit tricky to find the correct path, since FeO is slightly distorted such that it is monoclinic instead of cubic. Let us walk through the process of determining the MAE path step-by-step:
- Start from the conventional unit cell of the cubic structure and create a 2x2x2 supercell.
- Apply the $q=(1/2,1/2,1/2)$ spin wave and monoclinic distortion.
- Obtain the primitive magnetic unit cell.
- Identify the cubic [111] direction and define a rotation in the cubic $(1\bar{1}0)$ plane in the supercell of the distorted structure.
- Compute the list of SAXIS for a path with $\theta\in[0,\pi]$ for $\phi=0$ and $\phi=90$ .
7.2.1 cubic structure¶
You can download the Fm$\bar{3}$m FCC rock-salt structure of FeO from the materials project. The corresponding CIF is prepared in this example's directory: e07_FeO_bulk_MAE/FeO-rock-salt.cif. Read the file using the Python package ASE, build a supercell and plot it.
from ase.io import read, write
from ase.build import make_supercell
import numpy as np
import spglib
from ase import Atoms
from ase.visualize.plot import plot_atoms
# Load rock-salt structure
atoms = read("e07_FeO_bulk_MAE/FeO-rock-salt.cif")
# Build supercell that can hold the spin wave with q=(1/2,1/2,1/2)
q=np.array([0.5,0.5,0.5]) # AFM2
supercell = make_supercell(atoms, np.diag([1/q[0], 1/q[1], 1/q[2]]))
# Plot
plot_atoms(supercell,radii=0.3,rotation=('-20x,20y,-0z'))
7.2.2 Spin wave and distortion¶
Step 2 is to apply the $q=(1/2,1/2,1/2)$ spin wave and monoclinic distortion. Below you can see a function with different types of distortions that can be applyed to an FCC lattice. After the distortion, we can obtain the new space group using Spglib:
def apply_distortion(lattice, case):
#c.f. Schrön et al Table II
if case=="FeO":
# monoclinic distortion for FeO
r=0.0014; e=0.0253; t=0.0251
elif case=="MnO":
# rhombohedral distortion for MnO
e, t = 0, 0 ; r=-0.0063
elif case=="tetragonal":
# tetragonal distortion
r, e = 0, 0 ; t = 0.004
elif case=="orthorhombic":
# orthorhombic distortion
r = 0; e, t = 0.001, 0.001
elif case==None:
# No distortion
r, e, t = 0, 0, 0
strain_tensor=np.array([ [-0.5*t, 0.5*e + r, r], [0.5*e + r, -0.5*t, r], [r,r,t] ])
return (strain_tensor + np.diag([1, 1, 1])) @ lattice
# Apply distortion to lattice
lattice=supercell.get_cell()
distorted_lattice=apply_distortion(lattice, "FeO")
# Compute the volume for the cubic supercell and the monoclinic supercell
# Note that the volume is not perfectly preserved, but the change is small.
a1, a2, a3 = lattice
print("Supercell volume:", np.dot(a1,np.cross(a2,a3)))
a1, a2, a3 = distorted_lattice
print("Distorted supercell volume:", np.dot(a1,np.cross(a2,a3)))
# Create ASE structure for distorted supercell: distorted_supercell_atoms
distorted_supercell_atoms = Atoms(
cell=distorted_lattice,
scaled_positions=supercell.get_scaled_positions(),
numbers=supercell.get_atomic_numbers(),
pbc=True
)
# Determine space group of distorted structure
cell = (
distorted_supercell_atoms.get_cell(),
distorted_supercell_atoms.get_scaled_positions(),
distorted_supercell_atoms.get_atomic_numbers()
)
symmetry = spglib.get_spacegroup(cell, symprec=1e-5)
print("Space group:", symmetry)
# Plot
plot_atoms(distorted_supercell_atoms,radii=0.3,rotation=('-20x,20y,-0z'))
To apply the spin wave, we change the atomic type Fe to Mn at spin down sites as a place holder:
# Assign AFM moments (q = (1/2, 1/2, 1/2))
positions = distorted_supercell_atoms.get_scaled_positions()
numbers = []
for i, atom in enumerate(distorted_supercell_atoms):
if atom.symbol == "Fe":
R=2*positions[i][:]
sign = np.cos(2*np.pi*np.dot(R, q))
number = 25 if round(sign) < 0 else 26
else:
number = 8
numbers.append(number)
distorted_spinwave_atoms = Atoms(
cell=distorted_supercell_atoms.get_cell(),
scaled_positions=distorted_supercell_atoms.get_scaled_positions(),
numbers=numbers,
pbc=True
)
# Determine magnetic space group of distorted structure
cell = (
distorted_spinwave_atoms.get_cell(),
distorted_spinwave_atoms.get_scaled_positions(),
distorted_spinwave_atoms.get_atomic_numbers()
)
symmetry = spglib.get_spacegroup(cell, symprec=1e-5)
print("Space group:", symmetry)
plot_atoms(distorted_spinwave_atoms,radii=0.3,rotation=('-20x,20y,-0z'))
7.2.3 Primitive magnetic unit cell¶
While it is easy to vizualize the magnetic structure in the supercell of the conventional cell, it is computationlly cheaper to use the primitive cell. The transformation to obtain the primitive magnetic cell of the AFM struture is well-known and was used in Part 1 Example 4 for NiO for the cubic system. The same transformation is used for the distorted cell and the positions are obtained using: $$ A_{prim} = P \cdot A_{conv}. $$ Hence, positions transform via the inverse transpose: $$ \vec{v}_{prim}=(P^{-1})^T \cdot \vec{v}_{conv}. $$
The positions must be filtered, such that we get the list of atoms in the primitive unit cell. We can then check that the space group is still what we expect. We also revert Mn back to Fe and sort the atoms by atomic species in order to write a POSCAR file. Finally, we write the MAGMOM for the AFM structure corresponding to the POSCAR file and plot the structure.
lattice, positions, numbers = cell
def filter(pos,prec=8):
if np.all(np.round(pos,prec)>=0) and np.all(np.round(pos,prec)<1):
return pos
else:
pass
# New lattice vectors (c.f. Part 1 Example 4 for NiO)
a1 = 0.5 * (lattice[0] + 0.5 * (lattice[1] + lattice[2]))
a2 = 0.5 * (lattice[1] + 0.5 * (lattice[2] + lattice[0]))
a3 = 0.5 * (lattice[2] + 0.5 * (lattice[0] + lattice[1]))
new_lattice = np.array([a1, a2, a3])
print("Primitive lattice:")
print(new_lattice)
print("a:", np.linalg.norm(a1), "b:", np.linalg.norm(a2), "c:", np.linalg.norm(a3))
print("Volume:",np.dot(a1,np.cross(a2,a3)))
# New positions and numbers (i.e. apply inverse transpose of transformation and filter)
P = np.dot(new_lattice, np.linalg.inv(lattice))
P_inv_T = np.linalg.inv(P).T
new_positions = []
new_numbers = []
for i, pos_conv in enumerate(positions):
pos_prim = filter(P_inv_T @ pos_conv)
if pos_prim is None:
pass
else:
new_positions.append(pos_prim)
new_numbers.append(numbers[i])
print("Atom", i, pos_conv, filter(pos_prim), numbers[i])
# Make sure the space group of the primitive cell is still the desired space group
prim_cell = (
new_lattice,
new_positions,
new_numbers
)
symmetry = spglib.get_spacegroup(prim_cell, symprec=1e-5)
print("Space group:", symmetry)
# Sort atomic types, revert Mn to Fe and create corresponding MAGMOM
def revert_type(n):
if n==25:
n=26
return n
def get_m(n):
if n==8:
m=0
elif n==25:
m=4
elif n==26:
m=-4
return m
sorted_index = np.argsort(new_numbers)
new_positions = [new_positions[i] for i in sorted_index]
magmom = np.array([[0,0,get_m(new_numbers[i])] for i in sorted_index])
new_numbers = [revert_type(new_numbers[i]) for i in sorted_index]
print("MAGMOM = ", ' '.join( map(str,magmom[:,2] )) )
print("MAGMOM = ", ' '.join( map(str,magmom.flatten())) )
# ASE atoms structure for primitive magnetic cell
primitive_atoms = Atoms(
cell=new_lattice,
scaled_positions=new_positions,
numbers=new_numbers,
pbc=True
)
# Write POSCAR
primitive_atoms.write("e07_FeO_bulk_MAE/POSCAR.prim")
# Plot
plot_atoms(primitive_atoms,radii=0.3,rotation=('-20x,20y,-0z'))
To obtain a first understanding of the relative orientation of the primitive cell and the conventional supercell, plot them on top of each other! Additionally mark the cubic [1,1,1] and [1,1,0] direction.
import plotly.graph_objects as go
def add_cell_plot(fig,atoms,color,inv_scale):
# Plot lattice vectors as lines of distorted supercell
origin = np.zeros(3)
colors = ["red", "green", "blue"]
A = atoms.get_cell()
for i in range(3):
vec = A[i]
fig.add_trace(go.Scatter3d(
x=[origin[0], vec[0]],
y=[origin[1], vec[1]],
z=[origin[2], vec[2]],
mode='lines',
line=dict(color=colors[i], width=4),
name=f'a{i+1}'
))
# Plot atoms of distorted supercell
positions = [ np.dot(v,A) for v in atoms.get_scaled_positions()]
numbers = atoms.get_atomic_numbers()
for i, atom in enumerate(atoms):
fig.add_trace(go.Scatter3d(
x=[positions[i][0]], y=[positions[i][1]], z=[positions[i][2]],
mode='markers',
marker=dict(size=numbers[i]/inv_scale, color=color),
name=atom.symbol
))
def add_direction(fig, direction, name, color="black"):
fig.add_trace(go.Scatter3d(
x=[0, direction[0]],
y=[0, direction[1]],
z=[0, direction[2]],
mode='lines+text',
line=dict(color=color, width=5, dash='dash'),
name=name,
text=['', name],
textposition='top center'
))
# Plot
fig = go.Figure()
add_cell_plot(fig, distorted_supercell_atoms,"red",4)
add_cell_plot(fig, primitive_atoms,"blue",8)
add_direction(fig, distorted_supercell_atoms.get_positions()[6], "[1 1 1]")
add_direction(fig, distorted_supercell_atoms.get_positions()[3], "[1 1 0]")
# Layout and display
fig.update_layout(
scene=dict(
xaxis_title='X (Å)',
yaxis_title='Y (Å)',
zaxis_title='Z (Å)',
aspectmode='data'
),
title="Primitive cell and conventional supercell with lattice vectors and cubic crystallographic directions",
template="plotly_white"
)
fig.show()
7.2.4 Cubic [111] direction and (1-10) plane¶
Let us analyze how to identify the directions and planes you know in the conventional cell, also in the primitive cell.
The [1,1,1] direction refers to the direction that points from the origin to the point of fractional coordinates [h,k,l]. It is easy to identify that the [1,1,1] direction in the cubic supercell points from the origin to the 7th atom. Since the order was not changed when we applied the distortion, the [1,1,1] direction in the distorted supercell also points to the 7th atom.
The $(1\bar{1}0)$ plane refers to Miller Indices (hkl) in the conventional unit cell, which are inversely proportional to the intercepts with the crystal axes where $\bar{1}$ corresponds to $-1$. Therefore, the $(1\bar{1}0)$ plane:
- Cuts the x-axis at $1$ (in units of lattice constants),
- Cuts the y-axis at $−1$,
- Is parallel to the z-axis (since 0 means $\infty$ intercept).
Note that the [1,1,1] direction is in the $(1 \bar{1} 0)$ plane. Another direction in that plane is the [1,1,0] direction, where we can find the 4th atom in both the cubic and the distorted cell. The rotation axis to rotate a vector within the $(1 \bar{1} 0)$ plane is the [-1,1,0] direction. See the last plot in the previous section!
Determine the Cartesian coordinates of the [-1,1,0] direction in the distorted conventional supercell! How can the crystallographic direction be expressed in terms of the primitive lattice?
lattice_prim = primitive_atoms.get_cell()
lattice_distorted = distorted_spinwave_atoms.get_cell()
P = np.dot(lattice_prim, np.linalg.inv(lattice_distorted))
P_inv_T = np.linalg.inv(P).T
vec_conv = np.array([-1,1,0])
vec_prim = P_inv_T @ vec_conv
print( "Crystallographic", vec_conv, "direction in conventional cell: ")
print( " Cartesian coordinates ", np.dot(vec_conv,lattice_distorted))
print( "Crystallographic", np.round(vec_prim,4), "direction in primitiv cell: ")
print( " Cartesian coordinates ", np.dot(vec_prim,lattice_prim))
With this knowledge, we can define functions to obtain the Cartesian coordinates for a crystallographic direction of the cubic cell. And additional functions to obtain a unit vector and the angle between two vectors. Use these functions to obtain a unit vector along [1,1,1], [-1,1,0], [1,1,0] and [0,0,-1] direction and check the angles.
def cartesian_of(conventional_direction):
""" Returns the Cartesian coordinates of a crystallographic direction
defined in the conventional cell. """
return np.dot(conventional_direction,distorted_supercell_atoms.get_cell())
def unit_vector(vector):
""" Returns the unit vector of the vector. """
return vector / np.linalg.norm(vector)
def angle_between(v1, v2):
""" Returns the angle in radians between vectors 'v1' and 'v2'::
>>> angle_between((1, 0, 0), (0, 1, 0))
90.0
>>> angle_between((1, 0, 0), (1, 0, 0))
0.0
>>> angle_between((1, 0, 0), (-1, 0, 0))
180.0
"""
v1_u = unit_vector(v1)
v2_u = unit_vector(v2)
return np.arccos(np.clip(np.dot(v1_u, v2_u), -1.0, 1.0))/np.pi*180
e111 = unit_vector( cartesian_of( [1,1,1]) )
e_110 = unit_vector( cartesian_of([-1,1,0]) )
e110 = unit_vector( cartesian_of( [1,1,0]) )
e00_1 = unit_vector( cartesian_of( [0,0,-1]) )
print("Unit vector of [1,1,1] direction in Cartesian coordinates: ", e111 )
print("Unit vector of [-1,1,0] direction in Cartesian coordinates: ", e_110 )
print("Unit vector of [1,1,0] direction in Cartesian coordinates: ", e110 )
print("")
print("Check angle between [1,1,1] direction and [-1,1,0] direction is 90:", angle_between(e111, e_110))
print("Note: The angle theta that rotates SAXIS from [1,1,1] to [1,1,0] is", angle_between(e111, e110))
print("Note: The angle theta that rotates SAXIS from [1,1,1] to [0,0,-1] is", angle_between(e111, e00_1))
7.2.5 SAXIS along MAE path¶
The magnetization will be indirectly rotated by changing the spin quantization axis using the SAXIS tag. Read the documentation of SAXIS on the VASP Wiki!
For FeO, to increase the angle $\theta$ for $\phi=0$, we need to rotate e111 (unit vector along [1,1,1] direction) around e_110 (unit vector along [-1,1,0] direction). For $\phi=90$, the rotation axis e_110 first needs to be rotated by $\pi/2$ around e111 and subsequently $\theta$ can be increased by rotating about the new rotated_e_110. We use Scipy to create a list of values for SAXIS.
from scipy.spatial.transform import Rotation as R
print("Create saxis list for a rotation axis:", e_110)
# Define a rotation axis and step size for phi=0
rotvec = -e_110*5/180*np.pi
r = R.from_rotvec(rotvec)
# Create a list of 36 directions for the spin-quantization axis phi=0; theta in [0,180]
saxis_list = []
saxis = e111
for i in range(2*18+1):
#theta = round(angle_between(saxis,e111))
#print(theta, saxis)
saxis_list.append(saxis)
saxis = r.apply(saxis, inverse=True)
# Define and apply a rotation axis and step size of 90 to rotate to phi=90
rotvec1 = -e111*90/180*np.pi
r1 = R.from_rotvec(rotvec1)
rotated_e_110 = r1.apply(e_110, inverse=True)
print("Create saxis_list90 for a rotation axis:", rotated_e_110 )
# Define a second rotation axis and step size for phi=90
rotvec2 = -rotated_e_110*5/180*np.pi
r2 = R.from_rotvec(rotvec2)
# Create a list of 36 directions for the spin-quantization axis phi=90; theta in [0,180]
saxis_list90 = []
saxis = e111
for i in range(2*18+1):
saxis_list90.append(saxis)
saxis = r2.apply(saxis, inverse=True)
Plot the values for SAXIS in the unit cells! Look close to the origin to see the black and red planes corresponding to $\phi = 0 °$ and $\phi = 90 °$.
import plotly.graph_objects as go
# Plot
fig = go.Figure()
for i, saxis in enumerate(saxis_list):
add_direction(fig, saxis, " ", color="red")
add_direction(fig, saxis_list90[i], " ")
add_cell_plot(fig, distorted_supercell_atoms,"red",4)
add_cell_plot(fig, primitive_atoms,"blue",8)
# Layout and display
fig.update_layout(
scene=dict(
xaxis_title='X (Å)',
yaxis_title='Y (Å)',
zaxis_title='Z (Å)',
aspectmode='data'
),
title="Red: phi=0, black: phi=90",
template="plotly_white"
)
fig.show()
Finally, we print a bash array to write a run script that scans all values of SAXIS:
# uncomment for direct comparison
#for i, saxis in enumerate(saxis_list):
# theta = round(angle_between(saxis,e111))
# theta2 = round(angle_between(saxis_list90[i],e111))
# print(theta, saxis, theta2, saxis_list90[i])
def print_saxis_for_bash(saxis_list, str=""):
print( "declare -a saxis"+str+"=(")
for saxis in saxis_list:
theta = angle_between(saxis,e111)
print("\"SAXIS=", saxis[0], saxis[1], saxis[2] ," # theta= ", round(theta),"\"")
print(")")
print("# phi=0")
print_saxis_for_bash(saxis_list)
print("# phi=90")
print_saxis_for_bash(saxis_list90, "90")
7.2.6 VASP input files¶
Following the discussion above, we have a POSCAR file that contains the primitiv unit cell, MAGMOM for the $q=(1/2,1/2,1/2)$ AFM structure for both the collinear and the noncollinear calculation, and the MAE path in terms of a bash array that contains values for SAXIS. This is now combined to run a collinear LDA+U calculation and multiple ncl calculations along the path.
O Fe
1.0000000000000000
4.2032711140995955 2.2376441116844221 2.3233554199863908
2.2376441116844221 4.2032711140995955 2.3233554199863908
2.1592266873017891 2.1592266873017891 4.6099460038514311
O Fe
2 2
Cartesian
2.1500354782714526 2.1500354782714517 2.3141642109560534
6.4501064348143577 6.4501064348143577 6.9424926328681602
4.3000709565429034 4.3000709565429025 4.6283284219121050
0.0000000000000000 0.0000000000000000 0.0000000000000000
INCAR.ispin2
SYSTEM = FeO
ENCUT = 520
PREC = Acc
NBANDS = 32
ISPIN = 2
LNONCOLLINEAR = F
LSORBIT = F
MAGMOM = 0 0 4 -4
LASPH = T
GGA_COMPAT = F
EDIFF = 1E-8
ISMEAR = 0
SIGMA = 0.01
EFERMI = MIDGAP
ALGO = Normal
IMIX = 4
AMIX = 0.8
BMIX = 0.00001
AMIX_MAG = 0.8
BMIX_MAG = 0.00001
LMAXMIX = 4
XC = CA
LDAU = .TRUE.
LDAUTYPE = 2
LDAUL = -1 2
LDAUU = 0.00 4.00
LDAUJ = 0.00 0.00
INCAR.path
SYSTEM = FeO
ENCUT = 520
PREC = Acc
NBANDS = 64
ISPIN = 1
LNONCOLLINEAR = T
LSORBIT = T
MAGMOM = 0 0 0 0 0 0 0 0 4 0 0 -4
LASPH = T
GGA_COMPAT = F
EDIFF = 1E-8
ISMEAR = -5
EFERMI = MIDGAP
ALGO = Normal
IMIX = 4
AMIX = 0.8
BMIX = 0.00001
AMIX_MAG = 0.8
BMIX_MAG = 0.00001
LMAXMIX = 4
XC = CA
LDAU = .TRUE.
LDAUTYPE = 2
LDAUL = -1 2
LDAUU = 0.00 4.00
LDAUJ = 0.00 0.00
LORBMOM = T
LORBIT = 11
LCHARG = F
ICHARG = 11
k points
0
gamma
4 4 4
0 0 0
Pseudopotential of O and Fe.
Check out the tags that are set in the INCAR files. Restarting a noncollinear calculation from a collinear calculation can be done by reading the charge density and magnetization. This implies that the cutoff energy (set by ENCUT) must be the same, so that the FFT grid stays the same. The number of bands (set by NBANDS) in the collinear calculation is set for spin up and down. Thus, it should be twice the number of bands in the noncollinear run to span the same space. Technically it is possible to change the k-point density when restarting from a charge density, but here the charge density is fixed (ICHARG) so it is recommended to use the same k mesh.
7.3 Calculation¶
Open a terminal, navigate to this example's directory and run VASP by entering the following:
cd $TUTORIALS/magnetism/e07_*
bash run
where run contains
# spin-polarized collinear calculation
vasp_rm
cp INCAR.ispin2 INCAR
mpirun -np 4 vasp_std
cp CHGCAR CHGCAR.ispin2
cp OUTCAR OUTCAR.ispin2
rm WAVECAR
# phi=0
declare -a saxis=(
"SAXIS= 0.5626772625617759 0.5626772625617759 0.6056307426080456 # theta= 0 "
"SAXIS= 0.5978601694824203 0.5978601694824203 0.5339723171597044 # theta= 5 "
"SAXIS= 0.6284929995154631 0.6284929995154631 0.4582500399564773 # theta= 10 "
"SAXIS= 0.6543426183277444 0.6543426183277444 0.379040203250242 # theta= 15 "
"SAXIS= 0.6752122947116761 0.6752122947116761 0.29694564172654014 # theta= 20 "
"SAXIS= 0.6909431978285214 0.6909431978285214 0.21259114456861852 # theta= 25 "
"SAXIS= 0.701415606006982 0.701415606006982 0.12661870043448697 # theta= 30 "
"SAXIS= 0.706549817897407 0.706549817897407 0.039682611535587275 # theta= 35 "
"SAXIS= 0.7063067590471882 0.7063067590471882 -0.04755548599811423 # theta= 40 "
"SAXIS= 0.7006882792809389 0.7006882792809389 -0.13443165756858255 # theta= 45 "
"SAXIS= 0.689737138622211 0.689737138622211 -0.22028472305289984 # theta= 50 "
"SAXIS= 0.6735366818638969 0.6735366818638969 -0.3044612887832321 # theta= 55 "
"SAXIS= 0.6522102042640306 0.6522102042640306 -0.3863207202671715 # theta= 60 "
"SAXIS= 0.6259200131944265 0.6259200131944265 -0.46524001780304913 # theta= 65 "
"SAXIS= 0.5948661928835755 0.5948661928835755 -0.5406185578838422 # theta= 70 "
"SAXIS= 0.5592850816548526 0.5592850816548526 -0.6118826643047289 # theta= 75 "
"SAXIS= 0.5194474732491708 0.5194474732491708 -0.6784899741854025 # theta= 80 "
"SAXIS= 0.47565655592110295 0.47565655592110295 -0.7399335656790775 # theta= 85 "
"SAXIS= 0.4282456049931944 0.4282456049931944 -0.7957458159538323 # theta= 90 "
"SAXIS= 0.37757544642952146 0.37757544642952146 -0.8455019600847179 # theta= 95 "
"SAXIS= 0.32403171073223197 0.32403171073223197 -0.8888233237713169 # theta= 100 "
"SAXIS= 0.26802189806057375 0.26802189806057375 -0.9253802052778198 # theta= 105 "
"SAXIS= 0.2099722769086278 0.2099722769086278 -0.9548943846623136 # theta= 110 "
"SAXIS= 0.15032463994467993 0.15032463994467993 -0.9771412411985333 # theta= 115 "
"SAXIS= 0.08953294170225376 0.08953294170225376 -0.9919514628752192 # theta= 120 "
"SAXIS= 0.028059843712005143 0.028059843712005143 -0.9992123349627554 # theta= 125 "
"SAXIS= -0.03362680663188873 -0.03362680663188873 -0.9988685978403212 # theta= 130 "
"SAXIS= -0.09505753667289296 -0.09505753667289296 -0.9909228675549725 # theta= 135 "
"SAXIS= -0.15576482146250653 -0.15576482146250653 -0.9754356159119439 # theta= 140 "
"SAXIS= -0.21528664190741975 -0.21528664190741975 -0.9525247102476972 # theta= 145 "
"SAXIS= -0.27317000101378897 -0.27317000101378897 -0.9223645163883203 # theta= 150 "
"SAXIS= -0.32897437146788694 -0.32897437146788694 -0.885184571620306 # theta= 155 "
"SAXIS= -0.38227504831495773 -0.38227504831495773 -0.8412678377732034 # theta= 160 "
"SAXIS= -0.4326663812203666 -0.4326663812203666 -0.7909485477092377 # theta= 165 "
"SAXIS= -0.4797648617135846 -0.4797648617135846 -0.734609661609414 # theta= 170 "
"SAXIS= -0.5232120419192184 -0.5232120419192184 -0.6726799524153013 # theta= 175 "
"SAXIS= -0.5626772625617784 -0.5626772625617784 -0.6056307426080476 # theta= 180 "
)
# phi=90
declare -a saxis90=(
"SAXIS= 0.5626772625617759 0.5626772625617759 0.6056307426080456 # theta= 0 "
"SAXIS= 0.4989076889845989 0.6221645224170376 0.6033261347875017 # theta= 5 "
"SAXIS= 0.4313411266455501 0.6769167345834957 0.5964298507829445 # theta= 10 "
"SAXIS= 0.36049179788183533 0.726517201666274 0.5849943754797388 # theta= 15 "
"SAXIS= 0.28689890886534086 0.7705884341612915 0.5691067397498706 # theta= 20 "
"SAXIS= 0.21112254591808358 0.8087950233783235 0.5488878580944611 # theta= 25 "
"SAXIS= 0.13373941291711103 0.840846194103659 0.5244916084114024 # theta= 30 "
"SAXIS= 0.05533844222977383 0.8664980175750523 0.4961036608916411 # theta= 35 "
"SAXIS= -0.02348368741719698 0.8855552679268913 0.46394006495691376 # theta= 40 "
"SAXIS= -0.10212709202308479 0.8978729079769164 0.4282456049931939 # theta= 45 "
"SAXIS= -0.17999324779265274 0.9033571930467524 0.3892919373937097 # theta= 50 "
"SAXIS= -0.25648954626362414 0.901966384415516 0.3473755230897608 # theta= 55 "
"SAXIS= -0.331033804414907 0.8937110669766841 0.30281537130402303 # theta= 60 "
"SAXIS= -0.40305869543091627 0.8786540686806632 0.25595061169774136 # theta= 65 "
"SAXIS= -0.47201606640120186 0.85690998237615 0.20713791338923507 # theta= 70 "
"SAXIS= -0.537381110095081 0.8286442936893604 0.15674877048654495 # theta= 75 "
"SAXIS= -0.5986563590615508 0.79407212157849 0.10516667479295681 # theta= 80 "
"SAXIS= -0.6553754716569693 0.7534565811485512 0.05278419720281984 # theta= 85 "
"SAXIS= -0.7071067811865489 0.7071067811865495 -2.7755575615628914e-16 # theta= 90 "
"SAXIS= -0.7534565811485509 0.65537547165697 -0.05278419720282042 # theta= 95 "
"SAXIS= -0.7940721215784899 0.5986563590615516 -0.10516667479295741 # theta= 100 "
"SAXIS= -0.8286442936893604 0.5373811100950819 -0.1567487704865456 # theta= 105 "
"SAXIS= -0.8569099823761503 0.4720160664012028 -0.2071379133892358 # theta= 110 "
"SAXIS= -0.8786540686806638 0.4030586954309172 -0.2559506116977422 # theta= 115 "
"SAXIS= -0.8937110669766849 0.33103380441490793 -0.302815371304024 # theta= 120 "
"SAXIS= -0.901966384415517 0.256489546263625 -0.34737552308976183 # theta= 125 "
"SAXIS= -0.9033571930467538 0.17999324779265355 -0.38929193739371093 # theta= 130 "
"SAXIS= -0.8978729079769181 0.10212709202308544 -0.4282456049931953 # theta= 135 "
"SAXIS= -0.8855552679268932 0.02348368741719746 -0.4639400649569153 # theta= 140 "
"SAXIS= -0.8664980175750544 -0.05533844222977356 -0.4961036608916428 # theta= 145 "
"SAXIS= -0.8408461941036615 -0.13373941291711103 -0.5244916084114044 # theta= 150 "
"SAXIS= -0.8087950233783261 -0.2111225459180839 -0.5488878580944633 # theta= 155 "
"SAXIS= -0.7705884341612942 -0.28689890886534153 -0.5691067397498729 # theta= 160 "
"SAXIS= -0.7265172016662768 -0.3604917978818364 -0.5849943754797412 # theta= 165 "
"SAXIS= -0.6769167345834984 -0.4313411266455516 -0.5964298507829471 # theta= 170 "
"SAXIS= -0.6221645224170403 -0.4989076889846009 -0.6033261347875045 # theta= 175 "
"SAXIS= -0.5626772625617784 -0.5626772625617784 -0.6056307426080485 # theta= 180 "
)
# noncollinear calculations Phi=0 (restart from CHGCAR.ispin2)
echo "# theta/5, OSZICAR" > total_energy.dat
for i in 0 3 6 15 18 27 #{0..36}
do
cp INCAR.path INCAR.path_$i
echo """
${saxis[$i]} """ >> INCAR.path_$i
cp INCAR.path_$i INCAR
cp CHGCAR.ispin2 CHGCAR
mpirun -np 4 vasp_ncl
cp OSZICAR OSZICAR.path$i
cp OUTCAR OUTCAR.path$i
cp vaspout.h5 vaspout.path$i.h5
echo $i $(grep "F=" OSZICAR.path$i) >> total_energy.dat
done
cat total_energy.dat
# noncollinear calculations Phi=90 (restart from CHGCAR.ispin2)
echo "# theta/5, OSZICAR" > total_energy_90.dat
for i in 0 3 6 15 18 27 #{0..36}
do
cp INCAR.path INCAR.path90_$i
echo """
${saxis90[$i]} """ >> INCAR.path90_$i
cp INCAR.path90_$i INCAR
cp CHGCAR.ispin2 CHGCAR
mpirun -np 4 vasp_ncl
cp OSZICAR OSZICAR.path90-$i
cp OUTCAR OUTCAR.path90-$i
cp vaspout.h5 vaspout.path90-$i.h5
echo $i $(grep "F=" OSZICAR.path90-$i) >> total_energy_90.dat
done
cat total_energy_90.dat
import numpy as np
import py4vasp as pv
from scipy.optimize import curve_fit
filename='./e07_FeO_bulk_MAE/total_energy_90_ref.dat'
data=[]
with open(filename,"r") as f:
f.readline()
for line in f:
data.append([float(line.split()[0])*5,float(line.split()[3])])
data = np.array(data)
dgr, energies = data[:,0], data[:,1]-data[0,1]
def fit_func(x, K):
return K * np.sin(np.radians(x))**2
popt, _ = curve_fit(fit_func, dgr, energies)
K_fit = popt[0]
# Generate smooth curve
dgr_fit = np.linspace(min(dgr), max(dgr), 500)
energies_fit = fit_func(dgr_fit, K_fit)
p= pv.plot(dgr,energies, "phi=90")+pv.plot(dgr_fit,energies_fit, "phi=90 fit")
filename='./e07_FeO_bulk_MAE/total_energy_ref.dat'
data=[]
with open(filename,"r") as f:
f.readline()
for line in f:
data.append([float(line.split()[0])*5,float(line.split()[3])])
data = np.array(data)
dgr, energies = data[:,0], data[:,1]-data[0,1]
p= p+ pv.plot(dgr,energies, "phi=0", xlabel="θ in °", ylabel= "MAE in eV")
p
pa="e07_FeO_bulk_MAE"
p1=p
filename=pa+'/total_energy_90.dat'
data=[]
with open(filename,"r") as f:
f.readline()
for line in f:
data.append([float(line.split()[0])*5,float(line.split()[3])])
data = np.array(data)
dgr, energies = data[:,0], data[:,1]-data[0,1]
def fit_func(x, K):
return K * np.sin(np.radians(x))**2
popt, _ = curve_fit(fit_func, dgr, energies)
K_fit = popt[0]
# Generate smooth curve
dgr_fit = np.linspace(min(dgr), max(dgr), 500)
energies_fit = fit_func(dgr_fit, K_fit)
p1= p1+pv.plot(dgr,energies, "phi=90 calc")+pv.plot(dgr_fit,energies_fit, "phi=90 fit calc")
filename=pa+'/total_energy.dat'
data=[]
with open(filename,"r") as f:
f.readline()
for line in f:
data.append([float(line.split()[0])*5,float(line.split()[3])])
data = np.array(data)
dgr, energies = data[:,0], data[:,1]-data[0,1]
p1= p1+ pv.plot(dgr,energies, "phi=0 calc", xlabel="θ in °", ylabel= "MAE in eV")
p1
Click to see the answer!
No, because of the monoclinic distortion.
7.4 Questions¶
- Which coordinate system are the orbital moments defined in?
- Which tag controls the orientation of the spinor space with respect to Cartesian space?
- Which steps comprise the workflow for a MAE calculation?