Apparent inconsistency between TRIQS/solid_dmft and Amulet

[Vahid999]

Bulk Ca2N has been studied with DMFT using the Amulet code as reported in J. Phys. Chem. C 2021, 125, 15724. According to their Fig.4 , at \Gamma, DMFT is higher than DFT but along \Gamma-S0, DFT is higher and the two have different curvature. At T, DMFT predicts lower energies for the band crossing the Femi level.

In trying to reproduce the above results with solid_DMFT, I use 4 wannier functions: three non-interacting p-type for N2 and one single interacting s-like orbital crossing the Fermi level and accommodating a single electron. My .inp file is

` 0 12 12 12

1.0

1

0 0 0 1 0 0

7.2309`

and my .toml is

`[general]
seedname = "Ca2N"
jobname = "triqs"
csc = true

eta = 0.5
n_iw = 1025
n_tau = 10001

n_iter_dmft_first = 10
n_iter_dmft_per = 2
n_iter_dmft = 40

block_threshold = 1e-05

h_int_type = "simple_intra"
U = 2.5
J = 0.0
beta = 60
prec_mu = 0.0001

sigma_mix = 1.0
g0_mix = 1.0
dc_type = 0
dc = true
dc_dmft = true
calc_energies = true

h5_save_freq = 1

store_solver = false
enforce_off_diag = false
h_int_basis = "qe"

[solver]
type = "hubbardI"
n_l = 15
measure_G_l = false
measure_density_matrix = true

[dft]
dft_code = "qe"
n_cores = 64
mpi_env = "openmpi"
projector_type = "w90"
dft_exec = "pw.x"
w90_tolerance = 1e-2`

The attached band structure shows a comparison between my DFT and DMFT+DFT results. Comparing with the Amulet results, the single interacting band behaves differently between the two codes. Should we expect consistency between the two codes or am I the source of this inconsistency? I like to think the latter is true but am not sure what needs to change. Any input will be appreciated.

Thanks,
Vahid
Ca2N_elband.pdf

[the-hampel]

Hi @Vahid999 ,

I do not have time to study the paper you mentioned in detail, but I can assure you that you will be able to get qualitatively the same results as the AMULET code if all parameters are chosen similar and MLWFs match. The first thing that I noticed in your settings is that you are using the hubbardI solver. This will definitely give drastically different results compared to the AMULET code which uses some form of hybridization expansion solver. Please switch to the triqs cthyb solver type = "cthyb" as a first step. Then, are you sure that in the reference paper the p-orbitals are treated as non-interacting? Let me know if you need help with that. The rest of the parameters (dc, interaction etc) looks okay.

Best,
Alex

[Vahid999]

Hi Alex, I used solver type=cthyb. Here are the .inp and .toml files. The attached band structure shows very little change compared to DFT. The expected changes is much larger. 0 12 12 12 1.0 4 0 0 0 1 0 0 1 0 1 1 0 0 2 0 1 1 0 0 3 0 1 1 0 0 7.2309 [general] seedname = "Ca2N" jobname = "triqs" csc = true eta = 0.5 n_iw = 1025 n_tau = 10001 n_iter_dmft_first = 10 n_iter_dmft_per = 2 n_iter_dmft = 40 block_threshold = 1e-05 h_int_type = "simple_intra" U = 2.5 J = 0.0 beta = 60 prec_mu = 0.0001 sigma_mix = 1.0 g0_mix = 1.0 dc_type = 0 dc = true dc_dmft = true calc_energies = true h5_save_freq = 1 store_solver = false enforce_off_diag = false h_int_basis = "qe" [solver] type = "cthyb" length_cycle = 2000 n_warmup_cycles = 5e+3 n_cycles_tot = 1e+7 imag_threshold = 1e-10 perform_tail_fit = true fit_max_moment = 6 fit_min_w = 2 fit_max_w = 16 measure_density_matrix = true [dft] dft_code = "qe" n_cores = 64 mpi_env = "openmpi" projector_type = "w90" dft_exec = "pw.x" w90_tolerance = 1e-1 Best, Vahid On Jan 5, 2025, at 3:15 PM, Alexander Hampel ***@***.***> wrote: CAUTION: The Sender of this email is not from within Dalhousie. Hi @Vahid999<https://github.com/Vahid999> , I do not have time to study the paper you mentioned in detail, but I can assure you that you will be able to get qualitatively the same results as the AMULET code if all parameters are chosen similar and MLWFs match. The first thing that I noticed in your settings is that you are using the hubbardI solver. This will definitely give drastically different results compared to the AMULET code which uses some form of hybridization expansion solver. Please switch to the triqs cthyb solver type = "cthyb" as a first step. Then, are you sure that in the reference paper the p-orbitals are treated as non-interacting? Let me know if you need help with that. The rest of the parameters (dc, interaction etc) looks okay. Best, Alex — Reply to this email directly, view it on GitHub<#96 (comment)>, or unsubscribe<https://github.com/notifications/unsubscribe-auth/ARHYUBEESCDH2WCONK77QJD2JGAF3AVCNFSM6AAAAABUUEQYZGVHI2DSMVQWIX3LMV43OSLTON2WKQ3PNVWWK3TUHMZDKNZRG4ZDIMBZHE>. You are receiving this because you were mentioned.Message ID: ***@***.***>

[Vahid999]

Here is the band structure.

elband.pdf

[the-hampel]

Hi @Vahid999 ,

I think the input file looks reasonable to me now. But I am not sure I understand your plot. From cthyb you should have a frequency dependent self-energy on the Matsubara Axis. How did you get a single-particle band plot from that? I would have expected you to show a spectral function (colored plot) compared to the band structure. Similar to Fig 4 of your referenced paper. It might be also instructive to compare first the Matsubara self-energy similar to Fig 2 of the paper. Could you prepare such plots of your final self-energy?

Best,
Alex

[Vahid999]

Hi Alex,

Attached please see the self-energy versus frequency (Figure 2 of the paper). For plotting Figure 4 of the paper, I have not found a script yet on solid_dmft website. What I plotted was the Ca2N_band.dat_it39 on top of the DFT band structure. I assumed that Ca2N_band.dat_it39 would be the final band structure with DMFT.

Thanks,
Vahid
Figure2.pdf

[the-hampel]

HI Vahid,

okay the self energy looks very decent. Nicely converged. To get more out of the data I suggest to run the calculation with:

fit_min_w = 6
fit_max_w = 12

This will keep QMC data as is up to w=6 and not only 2 eV. Anyway this is just a minor point. Otherwise the results looks well converged.

However, you cannot use Ca2N_band.dat_it39 to assess any DMFT result. This is merely the input for DMFT, i.e. the change in the Wannier functions entering DMFT. The only thing you can read from this is that charge self-consistent steps did not have a significant influence on the result. To actually get access to the spectral function you have to analytically continue the self-energy, then calculate (on the real frequency axis):

$$\hat G(k,\omega) = \left[ \omega + \mu -\hat{\epsilon}_{DFT}(k)-\left( \hat\Sigma(\omega)^{imp}-\hat\Sigma^{dc} \right) \right]^{-1} $$

and then plot the spectral function (what is shown in Fig 4 of the reference paper):

$$ A(k,\omega) = - \frac{1}{\pi} \text{Im} \ G(k, \omega) $$

please take a look at our abinitio dmft tutorials in TRIQS: https://github.com/TRIQS/tutorials/tree/3.3.x/AbinitioDMFT . You can start with:
https://github.com/TRIQS/tutorials/blob/3.3.x/AbinitioDMFT/01-solid_dmft.ipynb
This will also guide you through continuing the self-energy.

However, doing this will require including charge self-consistency will require some extra steps. I would suggest you start with one shot calculations as done in the tutorial, i.e. when running the w90 converter set bloch_mode=False .

Best,
Alex