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 the QSGW stating point of FeSe which could be found in the directory:

cp /work3/training/questaal/data/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 follows:

dmft    proj=2 nlohi=3,20 betaev=55 nomb=2 knorm=1 nomega=999
    block:
        l=2
        sites=1
        sidxd=1 2 3 4 5
    block:
        l=2
        sites=2
        sidxd=6 7 8 9 10

Here we specify beta in 1/eV.

Set k-mesh

The k-mesh of the GW calculation and the k-mesh of used in DMFT (lmfdmft, which uses the lmf mesh) need to be aligned. Here we use a 4x4x4 GW mesh: we change the lmf mesh to match:

% const nkabc=4 nkgw=4

You could eventually change nomb which is the number of matsubara point where W is computed and betaev which is the inverse temperature in 1/eV.

Finally, use lmgw to compute $W_{loc}$

lmgw.sh --crpa --mpirun-n='mpirun -n 12' ctrl.fese

It will dump 2 files containing $W_{loc}$ on the real and matsubara axis.

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

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.3684708281800133
U_xy_xy_z2_z2 (w = 0) = 1.8265819812892294
U_xy(site1)_xy(site1)_xy(site2)_xy(site2) (w = 0) = 0.46928827233936765

(These values are slightly smaller than the converged result with a satisfactory k-mesh, such as 6x6x6 or 8x8x8.)

Invoking lmgw.sh with –wproj (instead of –crpa) calculates $W^{RPA}$ instead of $W^{cRPA}$ to allow comparison, writing instead to file “Wp_real.h5.”