Questaal Home
Navigation

cRPA

This tutorial explains how to estimate the Coulomb interaction U and Hund’s coupling J using the constrained RPA for FeSe.

This tutorial starts from qsgw stating point of FeSe which could be found in the lm directory :

cp path_to_lm_directory/dmft/test/fese/{ctrl,site,basp,rst,sigm}.fese .

Set dmft subspace

It is better to use the special mode KNORM=1 (just for the cRPA mode)

Each block must contain exactly one site (i.e equivalent sites need to be in different block). In our case of FeSe, the two Fe sites need to be as follow :

DMFT    PROJ=2 NLOHI=3,20 BETA=100 NOMEGA=999 KNORM=1
    BLOCK:
    L=2
    SITES=1
    SIDXD=1 2 3 4 5
    BLOCK:
    L=2
    SITES=2
    SIDXD=6 7 8 9 10

Set k-mesh

The k-mesh of the GW calculation and the k-mesh of used in DMFT need to be aligned. To do so, NKABC in GW token in ctrl file need to be changed by

GW    NKABC={nkabc}

and a 6x6x6 mesh will be used in this tutorial :

% const nkabc=6

then all the file from GW code are created using :

echo -1 |lmfgwd  fese
echo 'beta_crpa 2' >> GWinput
echo 'nomb_crpa 2' >> GWinput

you could eventually change nomb_crpa which is the number of matsubara point where W is computed and beta_crpa which is the inverse temperature (in eV)

Finally, use lmgw to compute

lmgw --crpa   --ht fese

It will dump 2 files containing on the real and matsubara axis.

h5ls Wp_crpa_real.h5
Wloc                     Dataset {29, 10, 10, 10, 10}
nw                       Dataset {1}
omega                    Dataset {29}

The first axis is for the frequency axis and the rest are the orbital indices : 2 atoms of 5 orbitals.

This file can easily be read by a simple python script :

import h5py
import numpy as np
f = h5py.File('Wp_crpa_real.h5', 'r')
W=np.array(f['Wloc']).real*13.6*2 # hatree to eV
Omega = np.array(f['omega'])*13.6*2 # hatree to eV
nw=W.shape[0]
print('U_xy_xy_xy_xy (w = 0) = {}'.format(W[0,0,0,0,0]))
print('U_xy_xy_z2_z2 (w = 0) = {}'.format(W[0,0,0,2,2]))
print('U_xy(site1)_xy(site1)_xy(site2)_xy(site2) (w = 0) = {}'.format(W[0,0,0,5,5]))

which should print :

U_xy_xy_xy_xy (w = 0) = 3.46373647456
U_xy_xy_z2_z2 (w = 0) = 1.91495383902
U_xy(site1)_xy(site1)_xy(site2)_xy(site2) (w = 0) = 0.427998487894