Skip to content

Commit

Permalink
GW embedding (#80)
Browse files Browse the repository at this point in the history 
This PR makes a couple changes on the GW embedding routine. The main changes include:
* Use dlr_wmax and dlr_eps in convert_gw_output()
* extracting H0_loc and Delta_time from g_weiss numerically on a DLR mesh.
* customized dlr_wmax and dlr_eps for crm_dyson_solver to avoid over-fitting

---------

Co-authored-by: Alexander Hampel <[email protected]>
  • Loading branch information
@cnyeh @the-hampel
cnyeh and the-hampel authored Jun 12, 2024
1 parent 45c8de8 commit d488ebe
Show file tree
Hide file tree
Showing 9 changed files with 177 additions and 61 deletions.
34 changes: 29 additions & 5 deletions python/solid_dmft/dmft_tools/solver.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -469,6 +469,14 @@ def solve(self, **kwargs):
self.triqs_solver.solve(h_int=self.h_int, **self.triqs_solver_params)
# *************************************

# dump Delta_tau constructed internally from cthyb when delta_interface = False
if self.general_params['store_solver'] and mpi.is_master_node():
with HDFArchive(self.general_params['jobname'] + '/' + self.general_params['seedname'] + '.h5',
'a') as archive:
if not self.solver_params['delta_interface']:
archive['DMFT_input/solver/it_-1'][f'Delta_time_{self.icrsh}'] = self.triqs_solver.Delta_tau
mpi.barrier()

# call postprocessing
self._cthyb_postprocessing()

Expand Down Expand Up @@ -1111,8 +1119,16 @@ def set_Gs_from_G_l():
mpi.report('\nCRM Dyson solver to extract Σ impurity\n')
# fit QMC G_tau to DLR
if mpi.is_master_node():
G_dlr = fit_gf_dlr(self.triqs_solver.G_tau,
w_max=self.general_params['dlr_wmax'], eps=self.general_params['dlr_eps'])
if self.solver_params['crm_dlr_wmax'] is not None:
dlr_wmax = self.solver_params['crm_dlr_wmax']
else:
dlr_wmax = self.general_params['dlr_wmax']
if self.solver_params['crm_dlr_eps'] is not None:
dlr_eps = self.solver_params['crm_dlr_eps']
else:
dlr_eps = self.general_params['dlr_eps']
mpi.report(f"crm_dyson_solver with (wmax, eps) = ({dlr_wmax}, {dlr_eps}). ")
G_dlr = fit_gf_dlr(self.triqs_solver.G_tau, w_max=dlr_wmax, eps=dlr_eps)
self.G_time_dlr = make_gf_dlr_imtime(G_dlr)

# assume little error on G0_iw and use to get G0_dlr
Expand Down Expand Up @@ -1438,8 +1454,16 @@ def set_Gs_from_G_l():
self.Sigma_Hartree = {}
# fit QMC G_tau to DLR
if mpi.is_master_node():
G_dlr = fit_gf_dlr(self.triqs_solver.results.G_tau,
w_max=self.general_params['dlr_wmax'], eps=self.general_params['dlr_eps'])
if self.solver_params['crm_dlr_wmax'] is not None:
dlr_wmax = self.solver_params['crm_dlr_wmax']
else:
dlr_wmax = self.general_params['dlr_wmax']
if self.solver_params['crm_dlr_eps'] is not None:
dlr_eps = self.solver_params['crm_dlr_eps']
else:
dlr_eps = self.general_params['dlr_eps']
mpi.report(f"crm_dyson_solver with (wmax, eps) = ({dlr_wmax}, {dlr_eps}). ")
G_dlr = fit_gf_dlr(self.triqs_solver.results.G_tau, w_max=dlr_wmax, eps=dlr_eps)
self.G_time_dlr = make_gf_dlr_imtime(G_dlr)
self.G_freq = make_gf_imfreq(G_dlr, n_iw=self.general_params['n_iw'])

Expand Down Expand Up @@ -1475,7 +1499,7 @@ def set_Gs_from_G_l():
gf, _, _ = minimize_dyson(G0_dlr=G0_dlr_iw[block],
G_dlr=G_dlr[block],
Sigma_moments=tail[block][0:1]
)
)
else:
for deg_shell in self.sum_k.deg_shells[self.icrsh]:
for i, block in enumerate(deg_shell):
Expand Down
47 changes: 34 additions & 13 deletions python/solid_dmft/gw_embedding/bdft_converter.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -90,6 +90,18 @@ def _get_dlr_from_IR(Gf_ir, ir_kernel, mesh_dlr_iw, dim=2):
return Gf_dlr


def check_iaft_accuracy(Aw, ir_kernel, stats,
beta, dlr_wmax, dlr_prec, data_name):
mpi.report('============== DLR mesh check ==============\n')
mpi.report(f'Estimating accuracy of the user-defined (wmax, eps) = '
f'({dlr_wmax}, {dlr_prec}) for the DLR mesh\n')
ir_imp_kernel = IAFT(beta=beta, lmbda=beta * dlr_wmax, prec=dlr_prec)
Aw_imp = ir_kernel.w_interpolate(Aw, ir_imp_kernel.wn_mesh('f'), 'f')

ir_imp_kernel.check_leakage(Aw_imp, stats, data_name, w_input=True)
mpi.report('=================== done ===================\n')


def calc_Sigma_DC_gw(Wloc_dlr, Gloc_dlr, Vloc, verbose=False):
r"""
Calculate the double counting part of the self-energy from the screened Coulomb interaction
Expand Down Expand Up @@ -288,7 +300,8 @@ def calc_W_from_Gloc(Gloc_dlr, U):
return W_dlr


def convert_gw_output(job_h5, gw_h5, it_1e=0, it_2e=0, ha_ev_conv = False):
def convert_gw_output(job_h5, gw_h5, dlr_wmax=None, dlr_eps=None,
it_1e=0, it_2e=0, ha_ev_conv = False):
"""
read bdft output and convert to triqs Gf DLR objects
Expand All @@ -298,6 +311,10 @@ def convert_gw_output(job_h5, gw_h5, it_1e=0, it_2e=0, ha_ev_conv = False):
path to solid_dmft job file
gw_h5: string
path to GW checkpoint file for AIMBES code
dlr_wmax: float
wmax for dlr mesh, defaults to the wmax from the IR basis
dlr_eps: float
precision for dlr mesh, defaults to the precision from the IR basis
it_1e: int, optional
iteration to read from gw_h5 calculation for 1e downfolding, defaults to last iteration
it_2e: int, optional
Expand Down Expand Up @@ -338,20 +355,21 @@ def convert_gw_output(job_h5, gw_h5, it_1e=0, it_2e=0, ha_ev_conv = False):
gw_data['beta'] = ar['imaginary_fourier_transform']['beta']
gw_data['lam'] = ar['imaginary_fourier_transform']['lambda']
gw_data['gw_wmax'] = gw_data['lam'] / gw_data['beta']
gw_data['gw_dlr_wmax'] = gw_data['gw_wmax'] if dlr_wmax is None else dlr_wmax
gw_data['number_of_spins'] = ar['system/number_of_spins']
assert gw_data['number_of_spins'] == 1, 'spin calculations not yet supported in converter'

prec = ar['imaginary_fourier_transform']['prec']
if prec == 'high':
# set to highest DLR precision possible
gw_data['gw_ir_prec'] = 1e-15
gw_data['gw_dlr_prec'] = 1e-13
gw_data['gw_dlr_prec'] = 1e-13 if dlr_eps is None else dlr_eps
elif prec == 'mid':
gw_data['gw_ir_prec'] = 1e-10
gw_data['gw_dlr_prec'] = 1e-10
gw_data['gw_dlr_prec'] = 1e-10 if dlr_eps is None else dlr_eps
elif prec == 'low':
gw_data['gw_ir_prec'] = 1e-6
gw_data['gw_dlr_prec'] = 1e-6
gw_data['gw_dlr_prec'] = 1e-6 if dlr_eps is None else dlr_eps

# 1 particle properties
g_weiss_wsIab = ar[f'downfold_1e/iter{it_1e}']['g_weiss_wsIab']
Expand Down Expand Up @@ -393,19 +411,26 @@ def convert_gw_output(job_h5, gw_h5, it_1e=0, it_2e=0, ha_ev_conv = False):
# get IR object
mpi.report('create IR kernel and convert to DLR')
# create IR kernel
mpi.report("")
ir_kernel = IAFT(beta=gw_data['beta'], lmbda=gw_data['lam'], prec=gw_data['gw_ir_prec'])

if dlr_wmax is not None or dlr_eps is not None:
# check user-defined dlr_wmax and dlr_eps for the impurity
check_iaft_accuracy(g_weiss_wsIab, ir_kernel, 'f',
gw_data['beta'], gw_data['gw_dlr_wmax'], gw_data['gw_dlr_prec'],
"fermionic Weiss field g")

gw_data['mesh_dlr_iw_b'] = MeshDLRImFreq(
beta=gw_data['beta']/conv_fac,
statistic='Boson',
w_max=gw_data['gw_wmax']*conv_fac,
w_max=gw_data['gw_dlr_wmax']*conv_fac,
eps=gw_data['gw_dlr_prec'],
symmetrize=True
)
gw_data['mesh_dlr_iw_f'] = MeshDLRImFreq(
beta=gw_data['beta']/conv_fac,
statistic='Fermion',
w_max=gw_data['gw_wmax']*conv_fac,
w_max=gw_data['gw_dlr_wmax']*conv_fac,
eps=gw_data['gw_dlr_prec'],
symmetrize=True
)
Expand Down Expand Up @@ -473,17 +498,13 @@ def convert_gw_output(job_h5, gw_h5, it_1e=0, it_2e=0, ha_ev_conv = False):
mpi.report(f'writing results in {job_h5}/DMFT_input')

with HDFArchive(job_h5, 'a') as ar:
if 'DMFT_results' in ar and 'iteration_count' in ar['DMFT_results']:
it = ar['DMFT_results']['iteration_count'] + 1
else:
it = 1
if 'DMFT_input' not in ar:
ar.create_group('DMFT_input')
if f'iter{it}' not in ar['DMFT_input']:
ar['DMFT_input'].create_group(f'iter{it}')
if f'iter{it_1e}' not in ar['DMFT_input']:
ar['DMFT_input'].create_group(f'iter{it_1e}')

for key, value in gw_data.items():
ar[f'DMFT_input/iter{it}'][key] = value
ar[f'DMFT_input/iter{it_1e}'][key] = value

mpi.report(f'finished writing results in {job_h5}/DMFT_input')
return gw_data, ir_kernel
Expand Down
77 changes: 49 additions & 28 deletions python/solid_dmft/gw_embedding/gw_flow.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -34,6 +34,8 @@
from triqs.utility import mpi
from triqs.gf.tools import inverse
from triqs.gf import (
iOmega_n,
fit_hermitian_tail,
Gf,
BlockGf,
make_hermitian,
Expand Down Expand Up @@ -253,6 +255,7 @@ def embedding_driver(general_params, solver_params, gw_params, advanced_params):
gw_data, ir_kernel = convert_gw_output(
general_params['jobname'] + '/' + general_params['seedname'] + '.h5',
gw_params['h5_file'],
general_params['dlr_wmax'], general_params['dlr_eps'],
it_1e = gw_params['it_1e'],
it_2e = gw_params['it_2e'],
)
Expand Down Expand Up @@ -358,21 +361,32 @@ def embedding_driver(general_params, solver_params, gw_params, advanced_params):

# prepare solver input
imp_eal = sumk.block_structure.convert_matrix(gw_params['Hloc0'][ish], ish_from=ish, space_from='sumk', space_to='solver')
delta_dlr = sumk.block_structure.convert_gf(gw_params['delta_dlr'][ish], ish_from=ish, space_from='sumk', space_to='solver')
# fill Delta_time from Delta_freq sumk to solver
for name, g0 in delta_dlr:
# make non-interacting impurity Hamiltonian hermitian
imp_eal[name] = (imp_eal[name] + imp_eal[name].T.conj())/2
for name, g0 in G0_dlr:
# Estimate the HF shift
G0_iw = solvers[ish].G0_freq[name]
Delta_iw = Gf(mesh=G0_iw.mesh, target_shape=G0_iw.target_shape)
Delta_iw << iOmega_n - inverse(G0_iw)
tail, err = fit_hermitian_tail(Delta_iw)
# overwrite H0 using estimation from high-frequency tail
imp_eal[name] = tail[0]
mpi.report(f"Tail fitting for extracting the HF shift in g_weiss with error {err}")

if mpi.is_master_node():
print('H_loc0[{:2d}] block: {}'.format(ish, name))
fmt = '{:11.7f}' * imp_eal[name].shape[0]
for block in imp_eal[name]:
print((' '*11 + fmt).format(*block.real))

G0_dlr_iw = make_gf_dlr_imfreq(g0)
Delta_dlr_iw = Gf(mesh=G0_dlr_iw.mesh, target_shape=g0.target_shape)
for iw in G0_dlr_iw.mesh:
Delta_dlr_iw[iw] = iw.value - inverse(G0_dlr_iw[iw]) - imp_eal[name]
Delta_dlr = make_gf_dlr(Delta_dlr_iw)
# create now full delta(tau)
Delta_tau = make_hermitian(make_gf_imtime(Delta_dlr, n_tau=general_params['n_tau']))

# without SOC delta_tau needs to be real
if not sumk.SO == 1:
# create now full delta(tau)
Delta_tau = make_hermitian(make_gf_imtime(delta_dlr[name], n_tau=general_params['n_tau']))
solvers[ish].Delta_time[name] << Delta_tau.real
else:
solvers[ish].Delta_time[name] << Delta_tau
Expand Down Expand Up @@ -444,37 +458,45 @@ def embedding_driver(general_params, solver_params, gw_params, advanced_params):
sumk.symm_deg_gf(Sigma_dlr_iw[ish],ish=ish)
Sigma_dlr[ish] = make_gf_dlr(Sigma_dlr_iw[ish])

for i, (block, gf) in enumerate(Sigma_dlr[ish]):
# print Hartree shift
print('Σ_HF {}'.format(block))
fmt = '{:11.7f}' * solvers[ish].Sigma_Hartree[block].shape[0]
for vhf in solvers[ish].Sigma_Hartree[block]:
print((' '*11 + fmt).format(*vhf.real))

# average Hartree shift if not magnetic
if not general_params['magnetic']:
if general_params['enforce_off_diag']:
solvers[ish].Sigma_Hartree['up_0'] = 0.5*(solvers[ish].Sigma_Hartree['up_0']+
solvers[ish].Sigma_Hartree['down_0'])
solvers[ish].Sigma_Hartree['down_0'] = solvers[ish].Sigma_Hartree['up_0']
else:
for iorb in range(gw_params['n_orb'][ish]):
solvers[ish].Sigma_Hartree[f'up_{iorb}'] = 0.5*(solvers[ish].Sigma_Hartree[f'up_{iorb}']+
solvers[ish].Sigma_Hartree[f'down_{iorb}'])
solvers[ish].Sigma_Hartree[f'down_{iorb}'] = solvers[ish].Sigma_Hartree[f'up_{iorb}']
for i, (block, gf) in enumerate(Sigma_dlr[ish]):
# print Hartree shift
print('Σ_HF {}'.format(block))
fmt = '{:11.7f}' * solvers[ish].Sigma_Hartree[block].shape[0]
for vhf in solvers[ish].Sigma_Hartree[block]:
print((' '*11 + fmt).format(*vhf.real))

# average Hartree shift if not magnetic
if not general_params['magnetic']:
if general_params['enforce_off_diag']:
solvers[ish].Sigma_Hartree['up_0'] = 0.5*(solvers[ish].Sigma_Hartree['up_0']+
solvers[ish].Sigma_Hartree['down_0'])
solvers[ish].Sigma_Hartree['down_0'] = solvers[ish].Sigma_Hartree['up_0']
else:
for iorb in range(gw_params['n_orb'][ish]):
solvers[ish].Sigma_Hartree[f'up_{iorb}'] = 0.5*(solvers[ish].Sigma_Hartree[f'up_{iorb}']+
solvers[ish].Sigma_Hartree[f'down_{iorb}'])
solvers[ish].Sigma_Hartree[f'down_{iorb}'] = solvers[ish].Sigma_Hartree[f'up_{iorb}']

iw_mesh = solvers[ish].Sigma_freq.mesh
# convert Sigma to sumk basis
Sigma_dlr_sumk = sumk.block_structure.convert_gf(Sigma_dlr[ish], ish_from=ish, space_from='solver', space_to='sumk')
Sigma_Hartree_sumk = sumk.block_structure.convert_matrix(solvers[ish].Sigma_Hartree, ish_from=ish, space_from='solver', space_to='sumk')
# store Sigma and V_HF in sumk basis on IR mesh
ir_nw_half = len(ir_mesh_idx)//2
for i, (block, gf) in enumerate(Sigma_dlr_sumk):
Vhf_imp_sIab[i,ish] = Sigma_Hartree_sumk[block]
for iw in range(len(ir_mesh_idx)):
Sigma_ir[iw,i,ish] = gf(iw_mesh(ir_mesh_idx[iw]))
# Make sure sigma_ir[iw].conj() = sigma_ir[-iw]
for n in range(ir_nw_half):
iw_pos = ir_nw_half+n
iw_neg = ir_nw_half-1-n
Sigma_ir[iw_pos,i,ish] = gf(iw_mesh(ir_mesh_idx[iw_pos]))
Sigma_ir[iw_neg,i,ish] = gf(iw_mesh(ir_mesh_idx[iw_pos])).conj()

if not general_params['magnetic']:
break
if mpi.is_master_node():
print("\nChecking impurity self-energy on the IR mesh...")
ir_kernel.check_leakage(Sigma_ir, stats='f', name="impurity self-energy", w_input=True)

# Writes results to h5 archive
if mpi.is_master_node():
Expand All @@ -495,7 +517,6 @@ def embedding_driver(general_params, solver_params, gw_params, advanced_params):
ar[f'downfold_1e/iter{iteration}']['Sigma_imp_wsIab'] = Sigma_ir
ar[f'downfold_1e/iter{iteration}']['Vhf_imp_sIab'] = Vhf_imp_sIab


mpi.report('*** iteration finished ***')
mpi.report('#'*80)
mpi.barrier()
Expand Down
Loading

0 comments on commit d488ebe

Please sign in to comment.