Questaal Home

Calculating isotropic Heisenberg exchange parameters via the enhanced susceptibility


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):


 XCFUN=0 001 007 ELIND=-1.0 NSPIN=2 TOL=1e-10


 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
 PLAT=sqrt(1/3) 0 {cbya/3}
      -sqrt(1/12)  1/2 {cbya/3}
      -sqrt(1/12) -1/2 {cbya/3}

 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

 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:
 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:
 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 “” 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

(3) correcting the phase

The map-results-irr-to-fbz executable (this is not built automatically in the master branch, please switch to the lm branch – this executable will then be built automatically and can be found in lm’s main build directory) reads in the susceptibility data and applies a phase to account for the geometrical positions of the atoms in the cell.

map-results-irr-to-fbz -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:

  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 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.
methodFe moment () (eV)

BFO Jij LDA and qsGW