Skip to content

Commit

Permalink
1. fix a bug when both improved_estimator and tail_fit are activated;…
Browse files Browse the repository at this point in the history
… 2. fix missing conflicts with the latest unstable branch; 3. check adapol installation
  • Loading branch information
cnyeh authored and the-hampel committed Sep 11, 2024
1 parent 012fdb5 commit 70a3117
Show file tree
Hide file tree
Showing 3 changed files with 38 additions and 44 deletions.
55 changes: 23 additions & 32 deletions python/solid_dmft/dmft_tools/solvers/ctseg_interface.py
Original file line number Diff line number Diff line change
Expand Up @@ -278,8 +278,24 @@ def set_Gs_from_G_l():
# get G_time, G_freq, Sigma_freq from G_l
set_Gs_from_G_l()
elif self.solver_params['perform_tail_fit']:
mpi.report('Self-energy post-processing algorithm: tail fitting with analytic static impurity self-energy')
self.Sigma_freq = inverse(self.G0_freq) - inverse(self.G_freq)
if not self.solver_params['improved_estimator']:
mpi.report('Self-energy post-processing algorithm: tail fitting with analytic static impurity self-energy')
self.Sigma_freq = inverse(self.G0_freq) - inverse(self.G_freq)
else:
mpi.report('Self-energy post-processing algorithm: '
'improved estimator + tail fitting with analytic static imppurity self-energy')
self.F_freq = self.G_freq.copy()
self.F_freq << 0.0
self.F_time = self.G_time.copy()
self.F_time << self.triqs_solver.results.F_tau
F_known_moments = make_zero_tail(self.F_freq, n_moments=1)
for i, bl in enumerate(self.F_freq.indices):
self.F_freq[bl] << Fourier(self.triqs_solver.results.F_tau[bl], F_known_moments[i])

for block, fw in self.F_freq:
for iw in fw.mesh:
self.Sigma_freq[block][iw] = self.F_freq[block][iw] / self.G_freq[block][iw]

# without any degenerate shells we run the minimization for all blocks
self.Sigma_freq, tail = self._fit_tail_window(
self.Sigma_freq,
Expand All @@ -295,42 +311,20 @@ def set_Gs_from_G_l():
self.G_freq = inverse(inverse(self.G0_freq) - self.Sigma_freq)

# if improved estimators are turned on calc Sigma from F_tau, otherwise:
elif self.solver_params['improved_estimator']:
if self.solver_params['perform_tail_fit']:
mpi.report('Self-energy post-processing algorithm: '
'improved estimator + tail fitting with analytic static imppurity self-energy')
else:
mpi.report('Self-energy post-processing algorithm: improved estimator')
elif self.solver_params['improved_estimator'] and not self.solver_params['perform_tail_fit']:
mpi.report('Self-energy post-processing algorithm: improved estimator')
self.F_freq = self.G_freq.copy()
self.F_freq << 0.0
self.F_time = self.G_time.copy()
self.F_time << self.triqs_solver.results.F_tau
F_known_moments = make_zero_tail(self.F_freq, n_moments=1)
if mpi.is_master_node():
for i, bl in enumerate(self.F_freq.indices):
self.F_freq[bl] << Fourier(self.triqs_solver.results.F_tau[bl], F_known_moments[i])
for i, bl in enumerate(self.F_freq.indices):
self.F_freq[bl] << Fourier(self.triqs_solver.results.F_tau[bl], F_known_moments[i])

self.F_freq << mpi.bcast(self.F_freq)
self.G_freq << mpi.bcast(self.G_freq)
for block, fw in self.F_freq:
for iw in fw.mesh:
self.Sigma_freq[block][iw] = self.F_freq[block][iw] / self.G_freq[block][iw]

# Further tail fitting on top of Sigma_freq from the improved estimator
if self.solver_params['perform_tail_fit']:
# without any degenerate shells we run the minimization for all blocks
self.Sigma_freq, tail = self._fit_tail_window(
self.Sigma_freq,
fit_min_n=self.solver_params['fit_min_n'],
fit_max_n=self.solver_params['fit_max_n'],
fit_min_w=self.solver_params['fit_min_w'],
fit_max_w=self.solver_params['fit_max_w'],
fit_max_moment=self.solver_params['fit_max_moment'],
fit_known_moments=self.Sigma_moments,
)
# recompute G_freq from Sigma_freq
self.G_freq = inverse(inverse(self.G0_freq) - self.Sigma_freq)

# should this be moved to abstract class?
elif self.solver_params['crm_dyson_solver']:
mpi.report('Self-energy post-processing algorithm: crm dyson solver')
Expand Down Expand Up @@ -408,9 +402,6 @@ def set_Gs_from_G_l():

if self.solver_params['measure_pert_order']:
self.perturbation_order_histo = self.triqs_solver.results.pert_order_Delta
bin_vec = np.arange(0, self.perturbation_order_histo.data.shape[0])
self.avg_pert_order = np.sum(bin_vec * self.perturbation_order_histo.data[:])
if mpi.is_master_node():
print(f'Average perturbation order: {self.avg_pert_order}')
self.avg_pert_order = self.triqs_solver.results.average_order_Delta

return
7 changes: 6 additions & 1 deletion python/solid_dmft/gw_embedding/bdft_converter.py
Original file line number Diff line number Diff line change
Expand Up @@ -181,7 +181,12 @@ def tail_fit_g_weiss(g_weiss_wsIab, ir_kernel, gw_data, wmax_imp=None, eps_imp=N

def bath_fitting(delta_wsIab, iw_mesh, Np=5):
mpi.report(f"Bath fitting for hybridization with nbath/impurity orbital = {Np}")
from adapol import hybfit
try:
from adapol import hybfit
except ImportError:
raise ImportError("bath fitting requires the adapol package (https://github.com/flatironinstitute/adapol). \n"
"Please ensure that it is installed. ")

nspin, nImp = delta_wsIab.shape[1:3]
error = -1
for s in np.arange(nspin):
Expand Down
20 changes: 9 additions & 11 deletions python/solid_dmft/gw_embedding/gw_flow.py
Original file line number Diff line number Diff line change
Expand Up @@ -233,8 +233,9 @@ def embedding_driver(general_params, solver_params, gw_params, advanced_params):
assert gw_params['code'] == 'aimbes', 'Only AIMBES is currently supported as gw code'

# prepare output h5 archive
archive = general_params['jobname']+'/'+general_params['seedname']+'.h5'
if mpi.is_master_node():
with HDFArchive(general_params['jobname'] + '/' + general_params['seedname'] + '.h5', 'a') as ar:
with HDFArchive(archive, 'a') as ar:
if 'DMFT_results' not in ar:
ar.create_group('DMFT_results')
if 'last_iter' not in ar['DMFT_results']:
Expand Down Expand Up @@ -336,7 +337,7 @@ def embedding_driver(general_params, solver_params, gw_params, advanced_params):

# write block structure to h5 archive
if gw_params['it_1e'] == 1 and mpi.is_master_node():
with HDFArchive(general_params['jobname'] + '/' + general_params['seedname'] + '.h5', 'a') as ar:
with HDFArchive(archive, 'a') as ar:
if 'block_structure' not in ar['DMFT_input']:
ar['DMFT_input']['block_structure'] = sumk.block_structure
if 'deg_shells' not in ar['DMFT_input']:
Expand Down Expand Up @@ -558,16 +559,13 @@ def embedding_driver(general_params, solver_params, gw_params, advanced_params):
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
mpi.report('Writing iter {} results to h5 archives {}/{}.h5 and {}.'.format(
iteration, general_params['jobname'], general_params['seedname'], gw_params['h5_file']))
mpi.report('Writing iter {} results to h5 archives {} and {}.'.format(iteration, archive, gw_params['h5_file']))
if mpi.is_master_node():
with HDFArchive(general_params['jobname'] + '/' + general_params['seedname'] + '.h5', 'a') as ar:
results_to_archive.write(ar, sumk, general_params, solver_params, solvers,
map_imp_solver, solver_type_per_imp, iteration,
False, gw_params['mu_emb'], density_mat_pre, density_mat)

results_to_archive.write(archive, sumk, general_params, solver_params, solvers,
map_imp_solver, solver_type_per_imp, iteration,
False, gw_params['mu_emb'], density_mat_pre, density_mat)
with HDFArchive(archive, 'a') as ar:
# store also IR / DLR Sigma
ar['DMFT_results/it_{}'.format(iteration)]['ir_mesh'] = ir_mesh
ar['DMFT_results/it_{}'.format(iteration)]['Sigma_imp_wsIab'] = Sigma_ir
Expand All @@ -576,7 +574,7 @@ def embedding_driver(general_params, solver_params, gw_params, advanced_params):
ar['DMFT_results/it_{}'.format(iteration)][f'Sigma_dlr_{ish}'] = Sigma_dlr[ish]

# write results to GW h5_file
with HDFArchive(gw_params['h5_file'],'a') as ar:
with HDFArchive(gw_params['h5_file'], 'a') as ar:
ar[f'downfold_1e/iter{iteration}']['Sigma_imp_wsIab'] = Sigma_ir
ar[f'downfold_1e/iter{iteration}']['Vhf_imp_sIab'] = Vhf_imp_sIab

Expand Down

0 comments on commit 70a3117

Please sign in to comment.