# 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, and then the calculation of the response function (“enhanced susceptibility”) and the corresponding Heisenberg exchange parameters. 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 (see main reference) is done by the “chi editor” mode of **lmf**.

Unfortunately the procedure for obtaining and tabulating 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
autoread: mto basis(4)
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 -mesh for the GW -mesh (used here for ).

```
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 -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 -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 for each -point with an 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 .)

For precise calculations, the 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 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 -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 -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 -points for ) - 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 -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 (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 -space) tabulation
- -magat=3,4 just relabels the output indexes to be consistent with the site numbering used by
**lmf**

#### (4) printing using **lmgf**

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

- add
*ASA=7*to section*VERS* - make the section
*OPTIONS*and add*ASA[]* - make the section
*GF*and add*MODE=11* - 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 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 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: , i.e. positive means ferromagnetic, 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.

method | Fe moment () | (eV) |
---|---|---|

LDA | 3.58 | 0.70 |

qsGW | 4.15 | 3.58 |