Magnetism with Questaal (DFT and QSGW)
Introduction
The bare spin susceptibility contains only Stoner (spin-flip) excitations; a Dyson equation relates the true response function to bare spin susceptibility:
The bare spin susceptibility is easily calculated via Alder-Wiser (“sum-over-states”) approach.
Projection of (bare) RPA transverse susceptibility on muffin-tin ( ) magnetisation gives the site-indexed susceptibility:
and we can define (isotropic) exchange interactions and with such that the spin waves , for local , which is a generally good approximation.
Because we project the susceptibility onto sites before inversion, is very easy to evaluate, but because the real part of the susceptibility is obtained using Kramers-Kronig from the imaginary part, it is necessary to calculate the imaginary part of the susceptibility over a broad frequency range.
At present, only the isotropic part of is available because the susceptibility maker does not currently support spin-orbit coupling.
Iron
DFT ground state calculation
Setup the crystal structure using a “.cif” file using ASE, requesting a primitive (not conventional cell):
cp /mnt/ceph-training/course_materials/cif/Fe.cif .
ase convert Fe.cif fe.poscar --read-args primitive_cell=True subtrans_included=False
poscar2init fe.poscar >init.fe
The “init” file is used to construct the real input file using blm (input generator). Additionally, we specify k-meshes and G-cut-off.
blm --nk=16 --gmax=9.9 --simple --gw~nk=6 --mag init.fe
The options specify:
- 161616 mesh of points in the Brillouin zone (tetrahedron integration is used by default in lmf)
- a g-cut-off of 9.9 a.u. defines the “smooth mesh” used for evalating various matrix elements in the full-potential method
- a “simple” formatting of the input file is requested
- 666 mesh of points in the Brillouin zone for evaluating the GW self-energy
- spin polarisation is enabled
The main input file is setup ctrl.ext .
The LMTO basis is setup and written to file basp.ext .
We cause a starting density to be spin polarised to break the spin symmetry; the moment is specified (by hand) as a list by , i.e. s,p,d,f starting moments.
sed -i -e 's/atom=Fe z= 26 r= 2.345586 lmx=3 lmxa=5/atom=Fe z= 26 r= 2.345586 lmx=3 lmxa=5 mmom=0 0 2/' ctrl.fe
The calculation now proceeds by:
- Solving the free atom problem using lmfa (atom solver). This provides the starting density.
lmfa ctrl.fe | tee llmfa
- Solving the band problem for the crystal using lmf (full-potential LMTO code).
mpirun lmf ctrl.fe | tee llmf_lda
At each iteration of the calculation, the total spin moment and spin moments within each atomic (augmentation) sphere are printed. In lmf the charges are represented by a sum of local quantities (expanded in spherical harmonics) and smooth-mesh quantities, the latter extending also into the interstitial.
For LDA:
mkrout: Qtrue sm,loc local true mm smooth mm local mm ovlp
1 7.040021 2.794699 4.245322 2.238593 0.387740 1.850853 2.153371
- Plotting bands along typical symmetry directions.
lmchk is a tool for analysing the structure and calculation setup, showing neighbour shells, bond angles, etc, and also sets up standard BZ paths for plotting bands.
lmchk --syml ctrl.fe | tee llmchk
lmf invoked with “–band~fn=syml” will calculate the eigenvalues along the specified path for plotting.
mpirun lmf --band~fn=syml ctrl.fe | tee llmf_band_lda
mv bnds.fe bnds_lda.fe
band-plot syml.fe bnds_lda.fe -e -8,8 -o fe_lda.pdf
Tutorials for partial DOS, Mulliken analysis, and core-level spectroscopy, can be found on this page.
Calculation of the Heisenberg exchange interaction parameters
Once the DFT ground state has been obtained, the Heisenberg interactions can be calculated. By default, interactions between sites with moments greater than 0.2 are calculated.
The calculation depends on evaluating the bare (Kohn-Sham) transverse susceptibility, on a mesh of -points spanning the irreducible Brillouin zone and a mesh of frequencies. These two meshes determine the quality of the calculation.
The transverse susceptibility is calculated using the GW machinery, invoked with lmgw.sh. The flags to control these meshes are: “GW_NKABC” (larger mesh is better) and “GW_DELRE” (smaller is better).
The calculation starts with a converged DFT calculation ( rst.ext as well as basp.ext and ctrl.ext files). It is easiest to start in a new directory.
mkdir jij-lda
cp {basp,ctrl,rst}.fe jij-lda/
cd jij-lda
sed -i -e 's/nkabc= 6/nkabc= 12 delre=0.001/' ctrl.fe
lmgw.sh --chipm --mpirun-1='OMP_NUM_THREADS=30 mpirun -n 1' --mpirun-n='OMP_NUM_THREADS=1 mpirun -n 30' ctrl.fe
The Heisenberg $J$ can be calculated from the static part of the bare transverse susceptibility: ; this is made efficient by projecting the susceptibility onto atomic sites.
lmf --x0j ctrl.fe
This produces a real-space tabulation of the Heisenberg to a range determined by the -mesh used.
The inverse of the site projected bare susceptibility is written to rsj-xinv.ext .
The convention used for the tabulation of the Heisenberg integrals is: , i.e. positive means ferromagnetic, are unit spin.
(An approximation to this inverse is which is the conventional “Lichtenstein formula” and , the inverse of the on-site susceptibility.)
ib jb vector J
1 1 0.0000000 0.0000000 0.0000000 -91.759933
1 1 1.4332500 1.4332500 -1.4332500 0.969301
1 1 1.4332500 -1.4332500 1.4332500 0.969301
1 1 -1.4332500 1.4332500 1.4332500 0.969301
1 1 1.4332500 1.4332500 1.4332500 0.969301
1 1 -1.4332500 -1.4332500 1.4332500 0.969301
1 1 -1.4332500 1.4332500 -1.4332500 0.969301
1 1 1.4332500 -1.4332500 -1.4332500 0.969301
1 1 -1.4332500 -1.4332500 -1.4332500 0.969301
1 1 2.8665000 0.0000000 0.0000000 0.750299
1 1 0.0000000 2.8665000 0.0000000 0.750299
1 1 0.0000000 0.0000000 2.8665000 0.750299
1 1 -2.8665000 0.0000000 0.0000000 0.750299
1 1 0.0000000 -2.8665000 0.0000000 0.750299
1 1 0.0000000 0.0000000 -2.8665000 0.750299
QSGW
QSGW calculations are performed using lmgw.sh (GW and BSE driver script).
The calculation proceeds by setting up a basis of wave-function products consisting of local parts expanded in spherical harmonics and interstitial parts in plane waves, calculating bare Coulomb integrals, the polarisability and finally the GW self-energy.
QSGW involves forming a new single-particle starting point from the interacting GW self-energy, and using this as a new starting point for GW. The converged QSGW’s eigenvalues/functions provide an ideal replacing the LDA (or other Vxc) starting point. The self-energy is stored in the file sigm.ext, which, together with the restart file ( rst.ext ) is enough to restart/use the converged calculation.
Copy the converged DFT solution to a new directory.
mkdir fe-qsgw
cp {ctrl,rst,basp,syml}.fe fe-qsgw/
cd fe-qsgw/
Invoke the QSGW driver, providing options (MPI launcher commands and options) to control the parallel execution of different parts of the calculation.
lmgw.sh --mpirun-1='OMP_NUM_THREADS=30 mpirun -n 1' --mpirun-n='OMP_NUM_THREADS=1 mpirun -n 30' ctrl.fe
This calculation requires 7 iterations for convergence. The QSGW moment:
mkrout: Qtrue sm,loc local true mm smooth mm local mm ovlp
1 6.994343 2.819829 4.174514 2.305116 0.383774 1.921342 2.137876
Repeat the band plotting using the QSGW potential
mpirun lmf --band~fn=syml ctrl.fe >llmf_band_qsgw
mv bnds.fe bnds_qsgw.fe
band-plot syml.fe bnds_qsgw.fe -e -8,8 -o fe_qsgw.pdf
band-plot syml.fe bnds_qsgw.fe:c=orange bnds_lda.fe:c=grey -e -8,8 -o fe_comparison.pdf
Repeat the calculation of the Heisenberg using the QSGW potential, taking care to remember to improve the k-mesh and frequency meshes. The GW_NKABC token determines the susceptibility mesh for the $J_{ij}$ calculation, which may be different to the mesh used for constructing the self-energy. Typically the susceptibility must be calculated on a finer BZ mesh, and convergence of the should be tested.
mkdir jij-qsgw
cp {basp,ctrl,rst,sigm}.fe jij-qsgw/
cd jij-qsgw
sed -i -e 's/nkabc= 6/nkabc= 12 delre=0.001/' ctrl.fe
lmgw.sh --chipm --mpirun-1='OMP_NUM_THREADS=30 mpirun -n 1' --mpirun-n='OMP_NUM_THREADS=30 mpirun -n 1' ctrl.fe
lmf --x0j ctrl.fe
This gives a different result because the sigm.ext file (and corresponding rst.ext file) provide a QSGW potential that, when present, automatically replace the DFT Vxc.
The QSGW interactions, from rsj_xinv.ext are changed from their LDA values, NN interactions becoming smaller while NNN grow slightly.
ib jb vector J
1 1 1.4332500 1.4332500 -1.4332500 0.355066
1 1 1.4332500 -1.4332500 1.4332500 0.355066
1 1 -1.4332500 1.4332500 1.4332500 0.355066
1 1 1.4332500 1.4332500 1.4332500 0.355066
1 1 -1.4332500 -1.4332500 1.4332500 0.355066
1 1 -1.4332500 1.4332500 -1.4332500 0.355066
1 1 1.4332500 -1.4332500 -1.4332500 0.355066
1 1 -1.4332500 -1.4332500 -1.4332500 0.355066
1 1 2.8665000 0.0000000 0.0000000 0.995604
1 1 0.0000000 2.8665000 0.0000000 0.995604
1 1 0.0000000 0.0000000 2.8665000 0.995604
1 1 -2.8665000 0.0000000 0.0000000 0.995604
1 1 0.0000000 -2.8665000 0.0000000 0.995604
1 1 0.0000000 0.0000000 -2.8665000 0.995604
Nickel oxide
The antiferromagnetic structure requires a supercell. This can be built from the primitive cell using lmscell which interactively prompts for the desired lattice vectors of the target supercell.
lmscell can be caused to write the new structure to a “site.ext” file. The control file can be edited to use the site file, replacing the structure or site blocks. Without the “–simple” option, blm will use “site.ext” by default.
We add the flags for the screened method in the blm and lmgw.sh calls to obtain more rapidly converging QSGW results with k-mesh.
Setup
Generate the primitive cell and supercell structure with AF configuration
cp /mnt/ceph-training/course_materials/cif/NiO.cif .
ase convert NiO.cif nio.poscar --read-args primitive_cell=True subtrans_included=False
poscar2init nio.poscar >init.nio
blm --nk=6 --gw~nk=3 --mag --gmax=11.2 --simple --rsham init.nio
lmchk ctrl.nio > llmchk_primitive
echo "2.089 2.089 4.178 2.089 4.178 2.089 4.178 2.089 2.089" | lmscell --wsite~fn=site ctrl.nio
See here for documentation on the supercell maker.
Edit the ctrl file to setup an AF configuration and to use the “site.ext” file for the structural information
# Autogenerated from init.nio using:
# blm --nk=6 --gw~nk=3 --mag --gmax=11.2 --simple --rsham init.nio
#symgrp i*r3(1,1,-1) r4x
#symgrp spglib
ham
gmax= 11.2
autobas[pnu=1 loc=1 mto=4 gw=1 pfloat=2,1]
forces= 0
nspin= 2
so= 0
xcfun="lm_lda_bh"
sigp[emax=3]
str
rmaxa=4.73698
env[mode=1 nel=2 el=-0.5 -1.2]
iter
nit= 99
mix= b2,b=0.3,k=7
conv= 1e-5
convc= 3e-5
bz
nkabc= 6
gw
nkabc= 3
gcutb=3.51
gcutx=2.86
mixbeta=1
struc
file=site
spec
atom=Ni z= 28 r= 2.192965 lmx=3 lmxa=5 idmod=0,0,1 mmom=0 0 2
atom=NiX z= 28 r= 2.192965 lmx=3 lmxa=5 idmod=0,0,1 mmom=0 0 -2
atom=O z= 8 r= 1.754673 lmx=3 lmxa=5
site
file=site
Swap one site from “up” to “down” by marking it “NiX” in the site.nio file:
sed -i -e 's/Ni 4.1780000 4.1780000 4.1780000/NiX 4.1780000 4.1780000 4.1780000/' site.nio
lmchk ctrl.nio > llmchk_afm
Add the new basis parameters and initial density for the new “NiX” species
echo " NiX RSMH= 1.462 1.462 0.961 0.961 EH= -0.5 -0.5 -0.5 -0.5 RSMH2= 1.462 1.462 0.961 EH2= -1.2 -1.2 -1.2 P= 4.743 4.547 3.895 4.12 5.102 6.078 PZ= 0 0 4.4" >> basp.nio
lmfa ctrl.nio
Run the QSGW calculation
lmgw.sh --mpirun-n='OMP_NUM_THREADS=1 mpirun -n 30' --cmd-gd='OMP_NUM_THREADS=2 mpirun -n 8' --lm-args='--v8 --tbeh' ctrl.nio
Optionally, you can explicitly add the AF symmetry by specifying in the control file “symgrp find afm:e::(1/2,1/2,1/2)”.