First tutorial on QSGW+DMFT
Setting up the DMFT loop
The spin-polarized QSGW starting point
This tutorial assumes you have terminated a spin-polarized QSGW calculation to be corrected with DMFT. Your QSGW calculation is supposed to be spin-polarized even for non-magnetic materials. For the purpose of this tutorial, we will refer to a QSGW calculation on ferromagnetic Nickel.
From this link you can download the files you will need to continue the tutorial on Nickel, if instead you want to produce them by your own, you can follow the commands below.
The file init.ni to start from scratch is:
LATTICE
SPCGRP=225
A=3.524 UNITS=A
SITE
ATOM=Ni X=0 0 0
SPEC
ATOM=Ni MMOM=0.0 0.0 0.6
To run a full QSGW calculation follow the commands below:
# the flag --mag for spin-polarized calculations
blm ni --gw --wsitex --mag --nit=20 --nk=12 --nkgw=8 --gmax=8.7
mv actrl.ni ctrl.ni
lmfa ni # Starting guess is the atomic density
mv basp0.ni basp.ni
lmf ni # DFT calculation. At convergence mmom = 0.59
lmfgwd ni --jobgw=-1 # GWinput
lmgwsc --wt --sym -maxit=20 --metal --getsigp --tol=2e-5 ni
The value of the parameters are a pretty low on purpose to run a QSGW loop in a reasonable time. We recommend to run the last step on a parallel machine (use the --openmp or the --mpi flag).
Note: Of course you can also do LDA+DMFT. The procedure is basically the same, but you can ignore all reference to any sigm file.
Input folders, files and programs
Once you have a converged spin-polarized QSGW calculation you still need some additional file to run lmfdmft and ctqmc. Input files can be download at this link.
In this tutorial, we use CTQMC solver developed by Kristjan Haule at Rutgers University which can be downloaded here. Once, the package is installed, make sure to add the executable directory in your $PATH. Wien2k package is not needed to run this tutorial. WIEN_DMFT_ROOT global variable has to be set :
export WIEN_DMFT_ROOT=/where/the/package/is/installed
export PATH=/where/the/package/is/installed:$PATH
Let qsgw the folder with the QSGW calculation and dmft-input the one where you extracted the content of the .tar.gz file linked above, then you dispatch relevant input files into two folders:
mkdir lmfinput qmcinput # input folders
cp qsgw/{ctrl,basp,site,rst}.ni lmfinput # copy relevant QSGW output files
cp qsgw/sigm.ni lmfinput/sigm_old
cp dmft-input/indmfl_input lmfinput/indmfl.ni
# copy files and programs for CTQMC runs
cp dmft-input/{atom_d.py,broad_sig.f90,Trans.dat,PARAMS} qmcinput/
Edit the ctrl file
You need to add some tokens to ctrl.ni.
# add some tokens at the end of ctrl.ni
cd lmfinput
echo 'DMFT PROJ=2 NLOHI=1,8 BETA=50 NOMEGA=2000 KNORM=0' >> ctrl.ni
The token DMFT_NLOHI defines the projection window in band index, DMFT_BETA is the inverse temperature in eV and DMFT_NOMEGA is the number of Matsubara frequencies in the mesh. Some details of the projection procedure are controlled by DMFT_PROJ and DMFT_KNORM, but you are not meant to change their value.
Moreover we recommend to add % const bxc0=0 and BXC0={bxc0} in the HAM section of ctrl.ni file. Setting HAM_BXC0 to TRUE, tells lmf to construct from the spin-averaged charge density. The reason for this will be clarified in the fourth tutorial.
At the end, you can see how your ctrl.ni should look like the following.
# Autogenerated from init.ni using:
# blm ni --gw --wsitex --mag --nit=20 --nk=12 --nkgw=8 --gmax=8.7
# Variables entering into expressions parsed by input
% const nit=20
% const met=5
% const nsp=2 so=0
% const lxcf=2 lxcf1=0 lxcf2=0 # for PBE use: lxcf=0 lxcf1=101 lxcf2=130
% const pwmode=0 pwemax=3 # Use pwmode=1 or 11 to add APWs
% const sig=12 gwemax=2 gcutb=3.3 gcutx=2.7 # GW-specific
% const nkabc=12 nkgw=8 gmax=8.7
% const bxc0=0
VERS LM:7 FP:7 # ASA:7
IO SHOW=f HELP=f IACTIV=f VERBOS=35,35 OUTPUT=*
EXPRESS
# Lattice vectors and site positions
file= site
# Basis set
gmax= {gmax} # PW cutoff for charge density
autobas[pnu=1 loc=1 lmto=5 mto=4 gw=1 pfloat=2]
# Self-consistency
nit= {nit} # Maximum number of iterations
mix= B2,b=.3,k=7 # Charge density mixing parameters
conv= 1e-5 # Convergence tolerance (energy)
convc= 3e-5 # tolerance in RMS (output-input) density
# Brillouin zone
nkabc= {nkabc} # 1 to 3 values. Use n1<0 => |n1| ~ total number
metal= {met} # Management of k-point integration weights in metals
# Potential
nspin= {nsp} # 2 for spin polarized calculations
so= {so} # 1 turns on spin-orbit coupling
xcfun= {lxcf},{lxcf1},{lxcf2} # set lxcf=0 for libxc functionals
#SYMGRP i r4x r3d
HAM
PWMODE={pwmode} PWEMIN=0 PWEMAX={pwemax} # For APW addition to basis
FORCES={so==0} ELIND=-0.7
RDSIG={sig} SIGP[EMAX={gwemax}] # Add self-energy to LDA
BXC0={bxc0}
GW NKABC={nkgw} GCUTB={gcutb} GCUTX={gcutx} DELRE=.01 .1
GSMEAR=0.003 PBTOL=1e-3
SPEC
ATOM=Ni Z= 28 R= 2.354453 LMX=3 LMXA=4 MMOM=0 0 0.6
DMFT PROJ=2 NLOHI=1,8 BETA=50 NOMEGA=2000 KNORM=0
Prepare spin-averaged self-energy
Although you have done a spin-polarized calculation, the starting point of the DMFT loop has to be non-magnetic. To do that you have to produce a spin-averaged sigm.ni.
cp sigm_old sigm.ni
# read sigm, make spin-average, write it down, and quit
lmf --rsig~spinav --wsig -vbxc0=1 ni > log
# rename sigm2: you will work with this spin-averaged sigm
mv sigm2.ni sigm.ni
cd ..
Check that at among the last lines of the log you find
replace sigma with spin average ...
and
Exit 0 done writing sigma, file sigm2
Compile the broadening program
The statistical noise of Quantum Monte Carlo calculations can be source of instabilities. Because of this, you need to broad the output of the ctqmc software at the end of each iteration.
In the file dmft-input.tar.gz you should have downloaded, you will find broad_sig.f90 which has precisely this purpose. However you can use whatever method you prefer (but be careful in not spoiling the low- and the high-frequency limits).
cd qmcinput
gfortran -o broad_sig.x broad_sig.f90
cd ..
The tutorial will continue assuming you are using broad_sig.x to broaden the impurity self-energy.
Prepare a vanishing impurity self-energy
You start the loop from scratch by creating an empty impurity self-energy:
mkdir siginp0
cd siginp0
cp ../lmfinput/* .
lmfdmft ni --ldadc=71.85 -job=1 -vbxc0=1 > log
cd ..
You can check that the line
Missing sigma file : create it and exit ...
is written at the end of the log.
The calculation has stopped just after reading the indmfl.ni. A text file called sig.inp has been created. It is formatted with the first column being the Matsubara frequencies (in eV) and then 0.0 repeated for a number of columns equal to twice the number of m channels (e.g. ten columns for d-type impurity grouped in real and imaginary parts).
Of course, if you want you can start from non-vanishing sig.inp files (e.g. from a previously converged DMFT loop).
…Ready to go!
The list of relevant files in the two input directories is
$ ls -l lmfinput
basp.ni # basis set used in the QSGW calculation
ctrl.ni # amended ctrl.ni file with DMFT category and HAM_BXC0={bxc0} token
indmfl.ni # instructions for the correlated subsystem
rst.ni # electronic density of the spin-polarized QSGW loop
sigm.ni # spin-averaged sigm.ni from converged QSGW loop
site.ni
$ ls -l qmcinput
atom_d.py # initialise the atomic problem in a d-electron system
broad_sig.x # broadens Sig.out at the end of each ctqmc run
PARAMS # main parameters for the ctqmc calculation
Trans.dat # translation table needed by atom_d.py
You are now ready to start the DMFT loop, following the link to the next tutorial.