# Calculating isotropic Heisenberg exchange parameters via the enhanced susceptibility

### Overview

Main reference: Kotani and MvS J. Phys: Cond. Mat. 2008

The calculation is essentially straightforward, consisting of two main steps: the calculation of the spin-dependent susceptibility, $\chi_0^{+-}$ and then the calculation of the response function (“enhanced susceptibility”) and the corresponding Heisenberg exchange parameters. $\chi_0^{+-}$ is obtained using a special mode of the GW susceptibility maker; the evaluation of the response via the TD-DFT-like approach with local kernel $U$ (see main reference) is done by the “chi editor” mode of lmf.

Unfortunately the procedure for obtaining and tabulating $J_{ij}$ needs a little cleaning up and a little care is needed in performing the calculation correctly. Here the example of BFO is presented step-by-step as a typical prototype.

### LDA calculation

Ferroelectric perovskite BFO in the RT hexagonal R3c structure (ICSD Collection Code 15299):

VERS LM=7 FP=7
SYMGRP find

HAM
AUTOBAS[MTO=14 LMTO=5 EH=-0.5]
XCFUN=0 001 007 ELIND=-1.0 NSPIN=2 TOL=1e-10
GMAX=12.1

BZ
METAL=5 NKABC=8

ITER
MIX=B6,b=0.5 NIT=50 CONV=1e-7

# a=b=5.5876,c=13.867AA
% const a=10.5591
% const cbya=2.481
STRUC
NBAS=10 NSPEC=4
ALAT={a}
PLAT=sqrt(1/3) 0 {cbya/3}
-sqrt(1/12)  1/2 {cbya/3}
-sqrt(1/12) -1/2 {cbya/3}

SPEC
ATOM=Bi  Z=83 R=2.6099 LMX=3 LMXA=4 KMXA=5 KMXV=5 A=0.01
ATOM=Fe  Z=26 R=2.0504 LMX=3 LMXA=4 KMXA=5 KMXV=5 A=0.01 MMOM=0 0 4 PZ=0 0 4.3
ATOM=FeX Z=26 R=2.0504 LMX=3 LMXA=4 KMXA=5 KMXV=5 A=0.01 MMOM=0 0 -4 PZ=0 0 4.3
ATOM=O   Z=8  R=1.6155 LMX=2 LMXA=4 KMXA=5 KMXV=5 A=0.01

SITE
ATOM=Bi  XPOS=0.0000000 0.0000000 0.0000000
ATOM=Bi  XPOS=0.5000000 0.5000000 0.5000000
ATOM=Fe  XPOS=0.2212000 0.2212000 0.2212000
ATOM=FeX XPOS=0.7212000 0.7212000 0.7212000
ATOM=O   XPOS=0.3973000 0.5233000 0.9423000
ATOM=O   XPOS=0.0233000 0.8973000 0.4423000
ATOM=O   XPOS=0.4423000 0.0233000 0.8973000
ATOM=O   XPOS=0.8973000 0.4423000 0.0233000
ATOM=O   XPOS=0.5233000 0.9423000 0.3973000
ATOM=O   XPOS=0.9423000 0.3973000 0.5233000


The magnetic ordering is collinear antiferromagnetic.

The details of the basis are not important here, set it up with:

lmfa --usebasp bfo |tee atom_log


Check the symmetry reported by lmf (invoke: lmf –quit=ham bfo):

[scarf587@cn755 LDA]\$ lmf --quit=ham bfo
---------------  START LMF - ref: b99ceea ----------------

rdctrl: reset global max nl from 4 to 5

rdctrl: reading basis parameters from file basp
ioorbp: read species Bi        RSMH,EH RSMH2,EH2
ioorbp: read species Fe        RSMH,EH RSMH2,EH2
ioorbp: read species FeX       RSMH,EH RSMH2,EH2
ioorbp: read species O         RSMH,EH RSMH2,EH2
reset nkaph from 2 to 3

LMF:      nbas = 10  nspec = 4  verb 35
pot:      spin-pol
libxc     LDA X (Slater exchange) + LDA C (Vosko, Wilk & Nusair (VWN5))
float:    float P LDA-style
bz:       metal(5), tetra, invit

Plat                                  Qlat
0.577350   0.000000   0.827000        1.154701   0.000000   0.403063
-0.288675   0.500000   0.827000       -0.577350   1.000000   0.403063
-0.288675  -0.500000   0.827000       -0.577350  -1.000000   0.403063
alat = 10.5591  Cell vol = 843.173308

LATTC:  as= 2.000   tol=1.00e-8   alat=10.55910   awald= 0.212   lmax=6
r1=  2.481   nkd= 85       q1=  3.661   nkg= 157

SGROUP: 1 symmetry operations from 0 generators
SYMLAT: Bravais system is rhombohedral with 12 symmetry operations.
SYMCRY: crystal invariant under 3 symmetry operations for tol=1e-5
GROUPG: the following are sufficient to generate the space group:
r3z
r3z
MKSYM:  found 3 space group operations; adding inversion generated 6 ops

SPLCLS:  4 species split into 6 classes
Species  Class      Sites...
Bi       1:Bi       1
5:Bi2      2
O        4:O        5 9 10
6:O2       6 7 8
BZMESH: 90 irreducible QP from 512 ( 8 8 8 )  shift= F F F


Converge the LDA:

mpirun -np 24 lmf bfo |tee scf_log


### Bare susceptibility calculation

#### (1) run lmfgwd

Prepare the system for invoking hx0fp0 (part of the GW code). If GW_NKABC is not defined, lmfgwd will default to using the DFT $k$-mesh for the GW $q$-mesh (used here for $\chi_0^{+-}(q)$ ).

lmfgwd --job=-1 bfo


note that, this time, lmfgwd reports:

 GROUPG: the following are sufficient to generate the space group:
r3z
r3z
MKSYM:  found 3 space group operations
no attempt to add inversion symmetry


and the number of $q$-points in the IBZ for the susceptibility calculation is 176 (and not 90 as in the LDA case).

#### (2) edit GWinput

GWinput requires modification before running the calculation, we need to specify:

• which sites are to be considered magnetic: we choose to restrict our attention to the Fe sites only
• we need to request a particular compatibility mode
• we need to request the susceptibility to be calculated at all of the $q$-points listed
sed -i "1s/.*/MagAtom 3 4/" GWinput               # specify magnetic sites
sed -i "s/GWversion   12/GWversion 0/" GWinput    # only NLF version supported
sed -i "s/QforEPSIBZ off/QforEPSIBZ on/" GWinput  # for evaluating response at all QIBZ


#### (3) run susceptibility maker

lmgw --mpi=24,24 --chixNLF bfo


This may take some time (probably an hour or two on a decent modern xeon) and it will produce a large number of files. The main log file for the calculation is lx0, and the main output consists of the files ChiPM0001.nlfc.mat, of which there should be 176. They contain $\chi_0^{+-}(\omega)$ for each $q$-point with an $\omega$ mesh as given by the usual GWinput specification. (The files are ASCII and the format is straightforward: number of sites, magnetic moments, positions, q-vector, omega, site indexed $\chi$.)

For precise calculations, the $\omega$ mesh for the susceptibility should be made (much) finer than is typical for qsGW calculations. Pay attention to the results of the pole search printed in the next step. The $\omega$ mesh is controlled by dw and omg_c in GWinput.

### Response calculation

#### (1) invoking the “chi editor”

lmf '--chimedit~new 176~read tkrs~export qx0q~exchange~a' bfo |tee chimedit_log


The quoted string contains options for the “chi editor” built into lmf. This editor can also be run interactively. All of the quantities discussed in the main reference paper are printed to the log file. We read all the 176 chi files prepared earlier which have a specific format (Takao Kotani Real Space) – the chi editor acts on each $q$-point individually and produces tabulations of the site-resolved susceptibility, response and exchange integral. Study the output! The remaining steps are only necessary in order to tabulate the data more usefully.

#### (2) mapping of IBZ to FBZ

To correctly map the IBZ to the full BZ, we ask lmf to provide the used list of $q$-points; a script is provided to do this that reads the printed output from lmf at a high print level. The script assumes the data we would like to map is in the file “results”.

• set the number of iterations in the lmf ctrl file to zero (NIT=0) to avoid any accidents/wasting time.
• add the flag (see listing below) –noinv to prevent lmf adding inversion symmetry (the list must be consistent with the $q$-points for $\chi$)
• print out the atomic positions (which may have been “shortened”, i.e. shifted closer to the origin using PBCs) to the file “pos.bfo”
• if you set GW_NKABC different to BZ_NKABC (the single-particle $k$-mesh) then change BZ_KNABC to be the same as the GW_NKABC you have used
• send the lmf output to the file “out.tmp”
• the script “preAllpoints.sh” is in the utils directory of your Questaal source folder. It produces files called “points”, “allpoints”, etc, which describe the IBZ-FBZ mapping
lmf --pr80 --wpos=pos --noinv bfo >out.tmp
ln -s Jmat_X0inv_w0eb.allq results
~/lm/lm/utils/preAllpoints.sh


#### (3) correcting the phase

The map-results-irr-to-fbz.x executable reads in the susceptibility data and applies a phase to account for the geometrical positions of the atoms in the cell.

~/lm/lm/utils/map-results-irr-to-fbz.x -takao3,2,3.58,3.58 -pos:0,0,0.54879720,-0.28867513,-0.5,0.13529720 -ft -magat=3,4 |tee phase_log

• “takao3” specifies a format; we have 2 magnetic sites and their absolute moments (look in scf_log) are both 3.5763 $\mu_B$ (even though we have an AF, we don’t use signs here)
• the positions are the Cartesian x1,y1,z1,x2,y2,z2 of the two Fe atoms: these are the 3rd and 4th lines of the file “pos.bfo” written in the step above
• -ft means that we would like the data to be Fourier transformed and the resulting “jr.dat” should be a real-space (and not $q$-space) tabulation
• -magat=3,4 just relabels the output indexes to be consistent with the site numbering used by lmf

#### (4) printing $J$ using lmgf

Add the necessary options to allow the ctrl file to be used by lmgf:

1. add ASA=7 to section VERS
2. make the section OPTIONS and add ASA[]
3. make the section GF and add MODE=11
4. add EMESH=10 1 0 0 0 0 to section BZ (the particular energy mesh will not be used)

GF_MODE=11 is a special mode that prints out the $J_{ij}$ tables and does some simple analysis.

• move the exchange data to a file with the expected extension: lmgf reads files jr.extn (extn=bfo)
• setup the structure constants needed by lmgf by invoking lmstr
• make sure you read the correct positions! (–rpos=pos)
• make sure that the symmetry is consistent (–noinv)
mv jr.dat jr.bfo
lmstr --noinv --rpos=pos bfo
lmgf  --noinv --rpos=pos --wrsj:tol=0 --sites:pair:3,4 --amoms=3.58,-3.58 bfo |tee jijprt_log

• option “wrsj” means to write a simplified tabulation to “rsj.bfo”, tol=0 means that even tiny values are printed
• we have only calculated $J_{ij}$ for sites 3 and 4 (the Fe sites), we restrict the pair printout to these sites (eg, there is no data for Fe-O couplings)
• option “amoms” informs lmgf of the moments, including sign

The convention used for the tabulation of the Heisenberg integrals is: $H = - \sum_{i\neq j} J_{ij} \vec{e}_i . \vec{e}_j$, i.e. positive means ferromagnetic, $\vec{e}_i$ are unit spins and note that there is no “factor of 1/2”. The units of the “rsj.bfo” file are lattice constants for the connecting vector and mRyd for the interactions.

• the LDA calculation can be replaced with a qsGW one: the gap opens and Fe moments become larger and more local. This results in a strong decrese in the Heisenberg exchange interactions, particularly for the nearest neighbour. Similar to the antiferromagnets CoO or La2CuO4 – the QSGW gap here is much too large. The BSE gap is much better.
methodFe moment ($\mu_B$)$E_g$ (eV)
LDA3.580.70
qsGW4.153.58