Transverse susceptibility and spin waves

This tutorial provides a basic demonstration of a spin wave calculation based on the transverse susceptibility; the main reference for the method and implementation can be found in reference1.

The site/magnetisation projected transverse (spin) susceptibility $\langle m(r_i)\|\chi^{+-}(q,\omega)\|m(r_j) \rangle$ provides a compact description of the magnetic response of the system that straightforwardly connects to atomistic models of magnetism and, combining the site-resolved magnetisation with the magnetic sum rules, the site projected susceptibility can be easily calculated from the non-interacting (bare) susceptibility using a TDDFT-like expression with a local kernel.

Questaal can calculate the site projected bare susceptibility efficiently using the $GW$ infrastructure; this can be post-processed to provide Heisenberg exchange integrals, $J^0_{ij}$ directly, using only the static part and identifying the $J=[\chi^{0+-}(q,\omega=0)]^{-1}$ . The inverse of the static part of the bare susceptibility can be shown to be approximately equivalent to the Lichtenstein formula implemented in various DFT codes.

The full susceptibility is related to the bare susceptibility by $\chi^{+-} = \chi^{0+-} - U$ . In TD-DFT, $U$ is the exchange-correlation kernel $f_{xc}$ and is static and local; an alternative is to define $U$ so as to ensure that the resulting $\chi^{+-}$ obeys the Goldstone mode (zero energy excitation at zero wavevector, assuming no spin-orbit interaction).

Enforcing that spin waves present in the susceptibility obey the Goldstone mode (zero frequency excitations with wave-vector zero corresponding to uniform rotation of the body in the absence of spin-orbit coupling) allows the kernel relating the full susceptibility to the bare susceptibility to be evaluated, assuming a rigid, spherical magnetisation at each site. The kernel used may be either static, full-filling the sum rule at $\omega=0$ , as is commonly done [2], or dynamic where the sum-rule is enforced at all frequencies, as derived by Kotani and van Schilfgaarde[1]. The latter agrees accurately with spin waves derived from the static Heisenberg exchange parameters.

Questaal has a tool (in lmf) that calculates the full susceptibility on a set of q-points provided in the same way as for plotting band structures. To do this, the rather smooth (in $q$ and $\omega$ ) bare susceptibility is approximated at each point in the path by explicit Fourier transform using the bare susceptibility provided on the set of irreducible Brillouin zone points. At each point on the path, the bare susceptibility is interpolated on a fine frequency mesh before transforming to the full susceptibility.

Overview of the calculation

The calculation proceeds identically in the case of DFT or QSGW (when a sigm.ext is present)

Generally, because the k-mesh is an important parameter that should be studied, it is sensible to start in a separate directory to any $GW$ calculation, etc, copying the appropriate sigm and rst files, etc. This is especially true because the process currently requires BZ_NKABC to be set to the same value as GW_NKABC and because the calculation typically uses a finer frequency mesh than is needed for QSGW.

Note that only sites with moments greater than $0.1\mu_B$ are considered magnetic and included in the calculation.

prepare a syml.ext file (perhaps by using lmchk –syml) specifying the q-path (Cartesian only) for plotting spin-waves
set the parameters GW_NKABC, the q-mesh for evaluating $\chi^{0+-}$ (bare susceptibility), and GW_DELRE its frequency binning
evaluate the bare susceptibility $\chi^{0+-}$ using lmgw.sh (adding –mpirun-n, etc, as usual)
edit the control file to ensure that BZ_NKABC matches GW_NKABC
post-process $\chi^{0+-}$ for $\chi^{RSA+-}$ using a special mode in lmf

    lmgw.sh --chipm ctrl.ext
    lmf --x0j~band~emax=0.5~ev ctrl.ext

Option band causes the calculation of the spin-waves from the static $J^0$ and the full susceptibility on the specified q-path for frequencies between 0 and emax Rydberg (default: 1 eV), with linear frequency mesh given by de (default: 0.1 meV). If the token ev is present, then emax and de are interpreted as electron-volt values. Options are lowercase. It is expedient to set mesh to cover only the area desired for plotting.

lmf –x0j also writes the real-space Heisenberg exchange, according to the convention: $H=-\sum_{i\neq j} J_{ij}e_i.e_j$ ; the range of interactions printed is determined by the initial q-point mesh for $\chi^0$ .

In general, as is the case for dielectric function calculations, the frequency mesh, specified by GW_DELRE may need to be made finer than the default used for QS $GW$ calculations and the results should be checked for convergence in this. GW_DELRE=0.001 is generally a reasonable choice. Similarly the k-sampling must be carefully checked; the necessary k-mesh will be larger than that used for the self-energy.

Options for –x0j mode

Options are separated using “~”

band calculate the RSA transverse susceptibility along a path specified in the syml.ext file (Cartesian q)

when band is set, the following options are used:

emax maximum frequency for calculating the RSA susceptibility, default is 0.073 Ryd (1 eV)
emin minimum frequency for calculating the RSA susceptibility, default is 0 Ryd
de frequency step for calculating the RSA transverse susceptibility, default is 7.3e-6 Ryd (0.1 meV)
ev switch to specify emax, emin, de in eV (lowercase only)
umod 1: dynamic RSA $U$ , 2: static RSA $U$

The Heisenberg parameters are always printed, but the spin-wave eigenvalues are only provided if band is specified.

Output files

file	produced by	contains
chipm.h5	lmgw.sh –chipm	bare susceptibility on IBZ q-points
spin_wave_lswt.ext	lmf –x0j~band	LSWT eigenvalues from $J^0$ on q path (syml.ext)
chipm_rsa.h5	lmf –x0j~band	RSA susceptibility on q path (syml.ext) and fine $\omega$ mesh
		RSA kernel $U(\omega)_i$
		lists of q and $\omega$ values
rsj-xinv.ext	lmf –x0j	real-space $J_{ij}$ table in simple format

Example 1: BCC Fe

Fe is the model itinerant, metallic ferromagnet. It shows sharp spin-waves at the centre of the Brillouin zone, which give rise to the spin-wave stiffness, but these become rapidly broadened as they reach the Stoner continuum. The resolution of the spin-wave spectrum requires a much finer k sampling for the susceptibility calculation than is typically necessary for insulators; similarly the resulting $J$ are much longer ranged. Spin wave excitations in iron extend over about 800 meV.

control file

% const nit=20
bz nkabc=20
gw delre=0.001 nkabc=10
ham xcfun='lm_lda_ca' nspin=2 autobas[mto=4 lmto=5]
site atom=Fe xpos=0 0 0
iter convc=1e-7 mix[b3,b=0.8] nit={nit}
spec atom=Fe z=26 lmx=4 lmxa=5 r/a={sqrt(3/4)/2} mmom=0 0 2
struc nbas=1 alat={2.867/0.529174} plat=-1/2 1/2 1/2 1/2 -1/2 1/2 1/2 1/2 -1/2

commands

First converge an LDA calculation:

lmfa ctrl.fe |tee llmfa
mv basp0.fe basp.fe 
lmf -vnit=20 ctrl.fe |tee llmf

The calculation should converge to a moment of $2.21\mu_B$ in 15 iterations.

 it 15 of 20       ehf=   -2541.1359472   ehk=   -2541.1359466
 From last iter    ehf=   -2541.1359472   ehk=   -2541.1359466
 diffe(q)=  0.0000000 (0.0000000)    tol= 0.0000100 (0.0000001)   more=F
c mmom=2.2130466 ehf=-2541.1359472 ehk=-2541.1359466

Prepare a new directory with the rst.fe file:

mkdir chipm
cp ctrl.fe basp.fe rst.fe chipm/

Change to the “chipm” directory and calculate the bare susceptibility. This is the most time consuming part and should be parallelised by invoking lmgw.sh using “mpirun-n” or “mpirun-pm” with many MPI processes.

lmgw.sh --chipm  ctrl.fe

For the post processing, obtain a suitable band plotting path and calculate the full susceptibility on fine frequency mesh upto 1 eV. You must edit the ctrl.fe file to set BZ_NKABC equal to GW_NKABC (10).

lmchk --syml ctrl.fe
lmf --x0j~band~emax=1.0~ev ctrl.fe

lmf outputs some information about the interpolation energy meshes, as well as the moment, the static limit of $U$ as well as the local $1/\chi^{0+-}_{ii}$ for each site.

 Rigid spin approximation spin waves utility (1 magnetic sites)
 chipm.h5 -> chipm_rsa.h5

 Dynamic sum-rule U

 data mesh: 4013 energy points (using 79), smallest positive value = 5.031e-4 eV
 interpolation mesh: 10001 energy points, h = 1.000e-4 eV, min = 0.000, max = 1.000 eV

 read qp lists from disk from file syml.ext
 383 qp on band path
 LSWT eigenvalues on qp path (from static J0) to file spin_wave_lwst.ext

 site=1 pos=0.0000 0.0000 0.0000 mom=2.246 U(0)=-0.9238 I=-0.9238 eV

 Heisenberg J = [X^{0,+-}(w=0)]^-1 to file rsj-xinv.ext
 writing RS J upto rmax=4.9237 alat (889 of 1000 vectors)
 Exit 0 exit x0j mode

(The sign of $U$ and $I$ depends on the definition of $\chi$ .)

plot

The chipm-rsa.h5 file contains the RSA $\chi^{+-}$ on the q-points specified by the syml file. It can easily be accessed using python. The imaginary part shows sharp (Lorenzian) peaks at the spin wave excitation frequencies; these are most easily illustrated if the heat map uses a logarithmic color scale; they become diffuse as they meet Stoner (spin-flip) excitations.

On top of the imaginary part of the susceptibility are plotted the spin waves given by the Heisenberg model (from the static limit of the bare susceptibility); they agree in where spin-wave remain well defined, in particular around the centre of the Brillouin zone. The fine details of the spectrum are not yet fully converged at nk=10.

Note that the damping at zero wavevector is almost identically zero.

#!/usr/bin/env python3
import h5py
import numpy as np
import matplotlib.pyplot as plt

f = h5py.File('chipm_rsa.h5', 'r')
swfile = 'spin_wave_lswt.fe'

ryd = 13.60569
try:
    ichi = np.transpose(np.log10(np.fabs(np.array(f['chirsa']).imag)))
    omega = np.array(f['realw']) * ryd
    iq = np.arange(0, np.array(f['nq'])[0])
    jij_sw_data = np.transpose(np.sort(np.loadtxt(swfile)))/1000.0  # sort eigenvalues to allow line plot, data is meV
    nsw = jij_sw_data.shape[0] - 3 # discard three columns holding q vector
except:
    pass
else:
    CS = plt.pcolormesh(iq, omega, np.trace(ichi), vmin=0.0, vmax=4, cmap='coolwarm', rasterized=True)
    colorbar = plt.colorbar(CS, label='log(trace(χ^{+-}))')

    for sw in range(nsw):
        plt.plot(iq,jij_sw_data[3+sw], 'g--', markersize=4, linewidth=2)

    plt.ylim([0,omega[-1]]) # only positive omega
    plt.xlabel('q index')
    plt.ylabel('ω (eV)')
    plt.show()
    #plt.savefig('chi-rsa.pdf')

Fe spin susceptibility (LDA, nk=10)

Unfortunately the special q point labels must be added manually at this time.

Example 2: MnO

MnO is a large moment antiferromagnetic insulator with well defined spin-waves throughout the Brillouin zone that have been well characterised by inelastic neutron scattering. Its spin wave spectrum is very similar to nickel oxide, but rather lower in energy, experimentally it is known that the maximum spin-wave energy is around 28 meV. We specify an antiferromagnetic configuration by distinguishing the two Mn sites with opposite initial moments and specify an antiferromagnetic symmetry:

For simplicity, the same k-mesh is used here for the DFT and susceptibility calculations.

control file

bz nkabc=4
gw delre=0.001
ham xcfun=1 nspin=2 autobas[mto=4 lmto=5]
iter convc=1e-6 mix[b3,b=0.9] nit=20
struc nbas=4 alat=8.40 plat=1/2 1/2 1 1/2 1 1/2 1 1/2 1/2
spec
 atom=Mn z=25 r=2.45 lmx=3 lmxa=4 mmom=0 0 1
 atom=Mnx z=25 r=2.45 lmx=3 lmxa=4 mmom=0 0 -1
 atom=o z=8 r=1.75 lmx=2 lmxa=4
site
 atom=Mn xpos=0 0 0
 atom=Mnx xpos=1/2 1/2 1/2
 atom=o xpos=1/4 1/4 1/4
 atom=o xpos=3/4 3/4 3/4
symgrp i*r3d r2(1,0,-1) afm:e::(1/2,1/2,1/2)

Note that GW_NKABC defaults to BZ_NKABC if not specified.

commands

It is appropriate to call lmgw.sh with the mpirun commands. Parallelisation of the susceptibility calculation will use mpirun-n or mpirun-pm.

Unfortunately, the LDA spin-wave size is of order $>$ 100 meV, so we choose the maximum energy for plotting to be 200 meV:

lmfa --usebasp ctrl.mno
lmchk --syml ctrl.mno
lmgw.sh --chipm  ctrl.mno
lmf --x0j~band~emax=0.2~ev ctrl.mno

plot

The MnO spin waves are well defined throughout the Brillouin zone and show antiferromagnetic (linear) dispersion around the zone centre. Here we plot the trace over the two sites.

#!/usr/bin/env python3
import h5py
import numpy as np
import matplotlib.pyplot as plt

f = h5py.File('chipm_rsa.h5', 'r')
swfile = 'spin_wave_lswt.mno'

ryd = 13.60569
try:
    ichi = np.transpose(np.log10(np.fabs(np.array(f['chirsa']).imag)))
    omega = np.array(f['realw']) * ryd
    iq = np.arange(0, np.array(f['nq'])[0])
    jij_sw_data = np.transpose(np.sort(np.loadtxt(swfile)))/1000.0  # sort eigenvalues to allow line plot, data is meV
    nsw = jij_sw_data.shape[0] - 3 # discard three columns holding q vector
except:
    pass
else:
    CS = plt.pcolormesh(iq, omega, np.trace(ichi), vmin=-20.0, vmax=10.0, shading='nearest', cmap='coolwarm', rasterized=True)
    colorbar = plt.colorbar(CS, label='log(trace(χ^{+-}))')

    for sw in range(nsw):
        plt.plot(iq,jij_sw_data[3+sw], 'g--', markersize=4, linewidth=2)

    plt.ylim([0,omega[-1]]) # only positive omega
    plt.xlabel('q index')
    plt.ylabel('ω (eV)')
    plt.show()
    #plt.savefig('chi-rsa.pdf')

MnO spin susceptibility (LDA, nk=4)

Exercises

How do the spectra change when you use QS $GW$ instead of the LDA? (Nickel oxide has the same structure as MnO, but with lattice constant a=7.88 Bohr; how does its spin wave spectrum look in the LDA, in QS $GW$ ?)
chipm-rsa.h5 also contains the bare susceptibility. Invoke lmf –x0j~band~emax=8.0~de=0.01~ev ctrl.fe (which plots over a large energy range with a coarser mesh) and adapt the plot routine to plot the bare susceptibility instead. What does it show? (Use h5ls to see the contents of the chipm-rsa.h5 file.)
To use a static $U$ , calculated using the sum rule at $\omega=0$ only, invoke lmf adding the option ~umod=2:

lmf --x0j~band~emax=1.0~ev~umod=2 ctrl.fe

Is the result very different for iron? For MnO?

What features are present for $\omega < 0$ in iron, in MnO?

Footnotes and References

¹ T. Kotani & M. van Schilfgaarde, J. Phys.: Condens. Matter 20 (2008) 295214

² S. Lounis, A. T. Costa, R. B. Muniz, and D. L. Mills, Phys. Rev. B 83, 035109