Skip to content

Commit

Permalink
Prototype for w2dyn_cthyb
Browse files Browse the repository at this point in the history
  • Loading branch information
@hmenke
hmenke committed Nov 10, 2023
1 parent 22e3125 commit 4c8c193
Show file tree
Hide file tree
Showing 3 changed files with 200 additions and 11 deletions.
3 changes: 3 additions & 0 deletions python/solid_dmft/dmft_tools/results_to_archive.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -98,6 +98,9 @@ def _compile_information(sum_k, general_params, solver_params, solvers,
write_to_h5['F_freq_{}'.format(icrsh)] = solvers[icrsh].F_freq
write_to_h5['F_time_{}'.format(icrsh)] = solvers[icrsh].F_time

if general_params['solver_type'] == 'w2dyn_cthyb':
pass # TODO

return write_to_h5

def write(archive, sum_k, general_params, solver_params, solvers, it, is_sampling,
Expand Down
146 changes: 145 additions & 1 deletion python/solid_dmft/dmft_tools/solver.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -231,6 +231,16 @@ def __init__(self, general_params, solver_params, advanced_params, sum_k, icrsh,
self.git_hash = triqs_ctseg_hash
self.version = version

elif self.general_params['solver_type'] == 'w2dyn_cthyb':
from w2dyn_cthyb.version import w2dyn_cthyb_hash, version

# sets up necessary GF objects on ImFreq
self._init_ImFreq_objects()
# sets up solver
self.triqs_solver = self._create_w2dyn_cthyb_solver()
self.git_hash = w2dyn_cthyb_hash
self.version = version

# ********************************************************************
# initialize Freq and Time objects
# ********************************************************************
Expand Down Expand Up @@ -264,7 +274,8 @@ def _init_ImFreq_objects(self):
or self.general_params['solver_type'] == 'cthyb' and self.general_params['legendre_fit']
or self.general_params['solver_type'] == 'ctseg' and self.solver_params['measure_gl']
or self.general_params['solver_type'] == 'ctseg' and self.general_params['legendre_fit']
or self.general_params['solver_type'] == 'hubbardI' and self.solver_params['measure_G_l']):
or self.general_params['solver_type'] == 'hubbardI' and self.solver_params['measure_G_l']
or self.general_params['solver_type'] == 'w2dyn_cthyb' and self.solver_params['measure_G_l']):

self.n_l = self.general_params['n_l']
self.G_l = self.sum_k.block_structure.create_gf(ish=self.icrsh, gf_function=Gf, space='solver',
Expand Down Expand Up @@ -610,6 +621,50 @@ def make_positive_definite(G):
# call postprocessing
self._ctseg_postprocessing()

if self.general_params['solver_type'] == 'w2dyn_cthyb':

# TODO: factor this block out of cthyb and w2dyn_cthyb
if self.general_params['cthyb_delta_interface']:
mpi.report('\n Using the delta interface for cthyb passing Delta(tau) and Hloc0 directly.')
# prepare solver input
sumk_eal = self.sum_k.eff_atomic_levels()[self.icrsh]
solver_eal = self.sum_k.block_structure.convert_matrix(sumk_eal, space_from='sumk', ish_from=self.sum_k.inequiv_to_corr[self.icrsh])
# fill Delta_time from Delta_freq sum_k to solver
for name, g0 in self.G0_freq:
self.Delta_freq[name] << iOmega_n - inverse(g0) - solver_eal[name]
known_moments = make_zero_tail(self.Delta_freq[name], 1)
tail, err = fit_hermitian_tail(self.Delta_freq[name], known_moments)
# without SOC delta_tau needs to be real
if not self.sum_k.SO == 1:
self.triqs_solver.Delta_tau[name] << make_gf_from_fourier(self.Delta_freq[name], self.triqs_solver.Delta_tau.mesh, tail).real
else:
self.triqs_solver.Delta_tau[name] << make_gf_from_fourier(self.Delta_freq[name], self.triqs_solver.Delta_tau.mesh, tail)


# Make non-interacting operator for Hloc0
Hloc_0 = Operator()
for spin, spin_block in solver_eal.items():
for o1 in range(spin_block.shape[0]):
for o2 in range(spin_block.shape[1]):
# check if off-diag element is larger than threshold
if o1 != o2 and abs(spin_block[o1,o2]) < self.solver_params['off_diag_threshold']:
continue
else:
# TODO: adapt for SOC calculations, which should keep the imag part
Hloc_0 += spin_block[o1,o2].real/2 * (c_dag(spin,o1) * c(spin,o2) + c_dag(spin,o2) * c(spin,o1))
self.solver_params['h_loc0'] = Hloc_0
else:
# fill G0_freq from sum_k to solver
self.triqs_solver.G0_iw << self.G0_freq

# Solve the impurity problem for icrsh shell
# *************************************
self.triqs_solver.solve(h_int=self.h_int, **{ **self.solver_params, **random_seed })
# *************************************

# call postprocessing
self._w2dyn_cthyb_postprocessing()

return

# ********************************************************************
Expand Down Expand Up @@ -798,6 +853,24 @@ def _create_ctseg_solver(self):

return triqs_solver

def _create_w2dyn_cthyb_solver(self):
r'''
Initialize w2dyn_cthyb solver instance
'''
from triqs_w2dyn_cthyb import Solver as w2dyn_cthyb_solver

gf_struct = self.sum_k.gf_struct_solver_list[self.icrsh]
# Construct the triqs_solver instances
if self.solver_params['measure_G_l']:
triqs_solver = w2dyn_cthyb_solver(beta=self.general_params['beta'], gf_struct=gf_struct,
n_iw=self.general_params['n_iw'], n_tau=self.general_params['n_tau'],
n_l=self.general_params['n_l'])
else:
triqs_solver = w2dyn_cthyb_solver(beta=self.general_params['beta'], gf_struct=gf_struct,
n_iw=self.general_params['n_iw'], n_tau=self.general_params['n_tau'])

return triqs_solver

def _create_ftps_solver(self):
r'''
Initialize ftps solver instance
Expand Down Expand Up @@ -1173,3 +1246,74 @@ def set_Gs_from_G_l():
self.perturbation_order = self.triqs_solver.histogram

return

def _w2dyn_cthyb_postprocessing(self):
r'''
Organize G_freq, G_time, Sigma_freq and G_l from w2dyn_cthyb solver
'''

def set_Gs_from_G_l():

# create new G_freq and G_time
for i, g in self.G_l:
g.enforce_discontinuity(np.identity(g.target_shape[0]))
# set G_freq from Legendre and Fouriertransform to get G_time
self.G_freq[i].set_from_legendre(g)
self.G_time[i].set_from_legendre(g)

# Symmetrize
self.G_freq << make_hermitian(self.G_freq)
self.G_freq_unsym << self.G_freq
self.sum_k.symm_deg_gf(self.G_freq, ish=self.icrsh)
self.sum_k.symm_deg_gf(self.G_time, ish=self.icrsh)
# Dyson equation to get Sigma_freq
self.Sigma_freq << inverse(self.G0_freq) - inverse(self.G_freq)

return

# get Delta_time from solver
self.Delta_time << self.triqs_solver.Delta_tau

# if measured in Legendre basis, get G_l from solver too
if self.solver_params['measure_G_l']:
# store original G_time into G_time_orig
self.G_time_orig << self.triqs_solver.G_tau
self.G_l << self.triqs_solver.G_l
# get G_time, G_freq, Sigma_freq from G_l
set_Gs_from_G_l()

else:
self.G_freq << make_hermitian(self.triqs_solver.G_iw)
self.G_freq_unsym << self.G_freq
self.sum_k.symm_deg_gf(self.G_freq, ish=self.icrsh)
# set G_time
self.G_time << self.triqs_solver.G_tau
self.sum_k.symm_deg_gf(self.G_time, ish=self.icrsh)

if self.general_params['legendre_fit']:
self.G_time_orig << self.triqs_solver.G_tau
# run the filter
self.G_l << legendre_filter.apply(self.G_time, self.general_params['n_l'])
# get G_time, G_freq, Sigma_freq from G_l
set_Gs_from_G_l()
elif self.solver_params['perform_tail_fit'] and not self.general_params['legendre_fit']:
# if tailfit has been used replace Sigma with the tail fitted Sigma from cthyb
self.Sigma_freq << self.triqs_solver.Sigma_iw
self.sum_k.symm_deg_gf(self.Sigma_freq, ish=self.icrsh)
else:
# obtain Sigma via dyson from symmetrized G_freq
self.Sigma_freq << inverse(self.G0_freq) - inverse(self.G_freq)

# if density matrix is measured, get this too
if self.solver_params['measure_density_matrix']:
self.density_matrix = self.triqs_solver.density_matrix
self.h_loc_diagonalization = self.triqs_solver.h_loc_diagonalization

if self.solver_params['measure_pert_order']:
self.perturbation_order = self.triqs_solver.perturbation_order
self.perturbation_order_total = self.triqs_solver.perturbation_order_total

if self.general_params['measure_G2_iw_ph']:
pass # TODO

return
62 changes: 52 additions & 10 deletions python/solid_dmft/read_config.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -109,6 +109,7 @@
* 'hubbardI'
* 'hartree'
* 'ctseg'
* 'w2dyn_cthyb'
n_iw : int, optional, default=1025
number of Matsubara frequencies
Expand Down Expand Up @@ -345,6 +346,29 @@
Sigma_iw will automatically be calculated via
http://dx.doi.org/10.1103/PhysRevB.85.205106
w2dyn_cthyb parameters
======================
measure_G2_iw_ph : bool, optional, default=False
Measure two-particle Green's function in particle-hole
frequency convention
measure_G2_n_bosonic : int, optional, default=30
Number of bosonic Matsubara frequencies for two-particle
Green's function measurement
measure_G2_n_fermionic : int, optional, default=30
Number of fermionic Matsubara frequencies for
two-particle Green's function measurement
worm_components : iterable of int, optional, default=None
Overrides worm components to measure with compound indices
(see ``w2dyn.auxiliaries.compound_index``) in iterable if not
None (Advanced setting: do not change if unsure)
move_global_prob : float, optional, default=0.005
Overall probability of the global moves
flavourchange_moves : int, optional, default=0
Use flavourchange-moves instead of 4-operator insertions;
they seemed to be inferior to 4 operator moves...
statesampling : int, optional, default=0
Activate state-sampling algorithm
hartree parameters
================
with_fock : bool, optional, default=False
Expand Down Expand Up @@ -499,7 +523,7 @@


'n_l': {'converter': int, 'valid for': lambda x, _: x > 0,
'used': lambda params: params['general']['solver_type'] in ['cthyb', 'inchworm', 'hubbardI', 'ctseg']
'used': lambda params: params['general']['solver_type'] in ['cthyb', 'inchworm', 'hubbardI', 'ctseg', 'w2dyn_cthyb']
and (params['solver']['measure_G_l'] or params['solver']['legendre_fit'])},

'n_iw': {'converter': int, 'valid for': lambda x, _: x > 0,
Expand Down Expand Up @@ -692,19 +716,19 @@
# cthyb parameters
#
'length_cycle': {'converter': int, 'valid for': lambda x, _: x > 0,
'used': lambda params: params['general']['solver_type'] in ['cthyb', 'ctint', 'ctseg']},
'used': lambda params: params['general']['solver_type'] in ['cthyb', 'ctint', 'ctseg', 'w2dyn_cthyb']},

'n_warmup_cycles': {'converter': lambda s: int(float(s)),
'valid for': lambda x, _: x > 0,
'used': lambda params: params['general']['solver_type'] in ['cthyb', 'ctint', 'ctseg']},
'used': lambda params: params['general']['solver_type'] in ['cthyb', 'ctint', 'ctseg', 'w2dyn_cthyb']},

'n_cycles_tot': {'converter': lambda s: int(float(s)),
'valid for': lambda x, _: x >= 0,
'used': lambda params: params['general']['solver_type'] in ['cthyb', 'ctint', 'ctseg']},
'used': lambda params: params['general']['solver_type'] in ['cthyb', 'ctint', 'ctseg', 'w2dyn_cthyb']},

'max_time': {'converter': int, 'valid for': lambda x, _: x >= 0,
'default': None,
'used': lambda params: params['general']['solver_type'] in ['cthyb', 'ctint', 'ctseg']},
'used': lambda params: params['general']['solver_type'] in ['cthyb', 'ctint', 'ctseg', 'w2dyn_cthyb']},

'imag_threshold': {'converter': float, 'default': None,
'used': lambda params: params['general']['solver_type'] in ['cthyb']},
Expand All @@ -716,13 +740,13 @@
'used': lambda params: params['general']['solver_type'] in ['cthyb']},

'measure_G_tau': {'converter': BOOL_PARSER, 'default': True,
'used': lambda params: params['general']['solver_type'] in ['hubbardI', 'ctseg']},
'used': lambda params: params['general']['solver_type'] in ['hubbardI', 'ctseg', 'w2dyn_cthyb']},

'measure_G_iw': {'converter': BOOL_PARSER, 'default': False,
'used': lambda params: params['general']['solver_type'] in ['ctseg']},

'measure_G_l': {'converter': BOOL_PARSER, 'default': False,
'used': lambda params: params['general']['solver_type'] in ['cthyb', 'hubbardI', 'ctseg']},
'used': lambda params: params['general']['solver_type'] in ['cthyb', 'hubbardI', 'ctseg', 'w2dyn_cthyb']},

'measure_density_matrix': {'converter': BOOL_PARSER, 'default': False,
'used': lambda params: params['general']['solver_type'] in ['cthyb', 'hubbardI']},
Expand All @@ -731,13 +755,13 @@
'used': lambda params: params['general']['solver_type'] in ['cthyb', 'ctint']},

'measure_pert_order': {'converter': BOOL_PARSER, 'default': False,
'used': lambda params: params['general']['solver_type'] in ['cthyb', 'ctint', 'ctseg']},
'used': lambda params: params['general']['solver_type'] in ['cthyb', 'ctint', 'ctseg', 'w2dyn_cthyb']},

'move_shift': {'converter': BOOL_PARSER, 'default': False,
'used': lambda params: params['general']['solver_type'] in ['cthyb']},

'random_seed': {'converter': str, 'default': None,
'used': lambda params: params['general']['solver_type'] in ['cthyb', 'ctint', 'ctseg']},
'used': lambda params: params['general']['solver_type'] in ['cthyb', 'ctint', 'ctseg', 'w2dyn_cthyb']},

'perform_tail_fit': {'converter': BOOL_PARSER,
'used': lambda params: params['general']['solver_type'] in ['cthyb']
Expand Down Expand Up @@ -863,6 +887,24 @@
'used': lambda params: params['general']['solver_type'] in ['ftps']},
'dmrg_tw': {'converter': float, 'valid for': lambda x, _: x > 0, 'default': 1e-9,
'used': lambda params: params['general']['solver_type'] in ['ftps']},
#
# w2dyn_cthyb parameters
#
'measure_G2_iw_ph': {'converter': BOOL_PARSER, 'default': False,
'used': lambda params: params['general']['solver_type'] in ['w2dyn_cthyb']},
'measure_G2_n_bosonic': {'converter': int, 'valid for': lambda x, _: x > 0, 'default': 30,
'used': lambda params: params['general']['solver_type'] in ['w2dyn_cthyb']},
'measure_G2_n_fermionic': {'converter': int, 'valid for': lambda x, _: x > 0, 'default': 30,
'used': lambda params: params['general']['solver_type'] in ['w2dyn_cthyb']},
'worm_components': {'converter': lambda s: list(map(int, s.split(','))), 'default': None,
'used': lambda params: params['general']['solver_type'] in ['w2dyn_cthyb']},
'move_global_prob': {'converter': float, 'default': 0.005,
'used': lambda params: params['general']['solver_type'] in ['w2dyn_cthyb']},
'flavourchange_moves': {'converter': int, 'default': 0,
'used': lambda params: params['general']['solver_type'] in ['w2dyn_cthyb']},
'statesampling': {'converter': int, 'default': 0,
'used': lambda params: params['general']['solver_type'] in ['w2dyn_cthyb']},

},
'advanced': {'dc_factor': {'converter': float, 'used': True, 'default': 'none'},

Expand Down Expand Up @@ -1235,7 +1277,7 @@ def read_config(config_file):
if isinstance(parameters['advanced']['map_solver_struct'], dict):
parameters['advanced']['map_solver_struct'] = [parameters['advanced']['map_solver_struct']]

if parameters['general']['solver_type'] in ['cthyb', 'ctint', 'ctseg']:
if parameters['general']['solver_type'] in ['cthyb', 'ctint', 'ctseg', 'w2dyn_cthyb']:
parameters['solver']['n_cycles'] = parameters['solver']['n_cycles_tot'] // mpi.size
del parameters['solver']['n_cycles_tot']

Expand Down

0 comments on commit 4c8c193

Please sign in to comment.