Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

updates on gw embedding and ctseg solver interface #85

Merged
merged 8 commits into from
Sep 11, 2024
Prev Previous commit
1. fix a bug when both improved_estimator and tail_fit are activated;…
… 2. fix missing conflicts with the latest unstable branch; 3. check adapol installation
  • Loading branch information
cnyeh committed Aug 30, 2024
commit 9570cab35f8bda35505addf45c82d88bca44b9f2
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
@@ -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,
@@ -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')
@@ -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
@@ -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):
20 changes: 9 additions & 11 deletions python/solid_dmft/gw_embedding/gw_flow.py
Original file line number Diff line number Diff line change
@@ -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']:
@@ -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']:
@@ -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
@@ -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