# %%
################################################################################
#
# solid_dmft - A versatile python wrapper to perform DFT+DMFT calculations
# utilizing the TRIQS software library
#
# Copyright (C) 2018-2020, ETH Zurich
# Copyright (C) 2021, The Simons Foundation
# authors: A. Carta, A. Hampel, M. Merkel, and S. Beck
#
# solid_dmft is free software: you can redistribute it and/or modify it under the
# terms of the GNU General Public License as published by the Free Software
# Foundation, either version 3 of the License, or (at your option) any later
# version.
#
# solid_dmft is distributed in the hope that it will be useful, but WITHOUT ANY
# WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A
# PARTICULAR PURPOSE. See the GNU General Public License for more details.
# You should have received a copy of the GNU General Public License along with
# solid_dmft (in the file COPYING.txt in this directory). If not, see
# <http://www.gnu.org/licenses/>.
#
################################################################################
# pyright: reportUnusedExpression=false
'''
abstract DMFT solver class for solid_dmft
'''
import numpy as np
from triqs.gfs import MeshImTime, MeshReTime, MeshReFreq, MeshLegendre, Gf, BlockGf
from triqs.operators import Operator
# import abc
from abc import ABC, abstractmethod
# circular import of all the solver subclasses
[docs]
class AbstractDMFTSolver(ABC):
"""
Abstract base class for DMFT solvers
This class defines the template for solvers for solving the DMFT impurity problem for the icrsh impurity.
Parameters
----------
general_paramuters: dict
general parameters as dict
solver_params: dict
solver-specific parameters as dict
sum_k: triqs.dft_tools.sumk object
SumkDFT instance
icrsh: int
correlated shell index
h_int: triqs.operator object
interaction Hamiltonian of correlated shell
iteration_offset: int
number of iterations this run is based on
"""
# ********************************************************************
# General initialization common to all solvers
# ********************************************************************
[docs]
def __init__(
self, general_params, solver_params, sum_k, icrsh, h_int, iteration_offset, deg_orbs_ftps, gw_params=None, advanced_params=None
):
self.general_params = general_params
self.solver_params = solver_params
self.gw_params = gw_params
self.advanced_params = advanced_params
self.sum_k = sum_k
self.icrsh = icrsh
self.h_int = h_int
self.iteration_offset = iteration_offset
self.deg_orbs_ftps = deg_orbs_ftps
if self.solver_params.get('crm_dyson_solver'):
self.G_time_dlr = None
self.Sigma_dlr = None
if self.solver_params.get('measure_density_matrix'):
self.density_matrix = None
self.h_loc_diagonalization = None
if self.solver_params.get('measure_chi') is not None:
self.O_time = None
if self.solver_params.get('delta_interface'):
self.Hloc_0 = Operator()
# ********************************************************************
# initialize Freq and Time objects
# ********************************************************************
[docs]
def _init_ImFreq_objects(self):
r"""
Initialize all ImFreq objects
"""
# create all ImFreq instances
self.n_iw = self.general_params['n_iw']
self.G_freq = self.sum_k.block_structure.create_gf(ish=self.icrsh, gf_function=Gf, space='solver', mesh=self.sum_k.mesh)
# copy
self.Sigma_freq = self.G_freq.copy()
self.G0_freq = self.G_freq.copy()
self.G_freq_unsym = self.G_freq.copy()
self.Delta_freq = self.G_freq.copy()
# create all ImTime instances
self.n_tau = self.general_params['n_tau']
self.G_time = self.sum_k.block_structure.create_gf(
ish=self.icrsh, gf_function=Gf, space='solver', mesh=MeshImTime(beta=self.sum_k.mesh.beta, statistic='Fermion', n_tau=self.n_tau)
)
# copy
self.Delta_time = self.G_time.copy()
# create all Legendre instances
if self.solver_params.get('measure_G_l') or self.solver_params.get('legendre_fit'):
self.n_l = self.solver_params['n_l']
self.G_l = self.sum_k.block_structure.create_gf(
ish=self.icrsh,
gf_function=Gf,
space='solver',
mesh=MeshLegendre(beta=self.general_params['beta'], max_n=self.n_l, statistic='Fermion'),
)
# move original G_freq to G_freq_orig
self.G_time_orig = self.G_time.copy()
[docs]
def _init_ReFreq_objects(self):
r"""
Initialize all ReFreq objects
"""
# create all ReFreq instances
self.n_w = self.general_params['n_w']
self.G_freq = self.sum_k.block_structure.create_gf(ish=self.icrsh, gf_function=Gf, space='solver', mesh=self.sum_k.mesh)
# copy
self.Sigma_freq = self.G_freq.copy()
self.G0_freq = self.G_freq.copy()
self.Delta_freq = self.G_freq.copy()
self.G_freq_unsym = self.G_freq.copy()
# create another Delta_freq for the solver, which uses different spin indices
n_orb = self.sum_k.corr_shells[self.icrsh]['dim']
n_orb = n_orb // 2 if self.sum_k.corr_shells[self.icrsh]['SO'] else n_orb
gf = Gf(target_shape=(n_orb, n_orb), mesh=MeshReFreq(n_w=self.n_w, window=self.general_params['w_range']))
self.Delta_freq_solver = BlockGf(name_list=tuple([block[0] for block in self.gf_struct]), block_list=(gf, gf), make_copies=True)
# create all ReTime instances
# FIXME: dummy G_time, since time_steps will be recalculated during run
# time_steps = int(2 * self.solver_params['time_steps'] * self.solver_params['refine_factor']) if self.solver_params['n_bath'] != 0 else int(2 * self.solver_params['time_steps'])
time_steps = int(2 * 1 * self.solver_params['refine_factor']) if self.solver_params['n_bath'] != 0 else int(2 * 1)
self.G_time = self.sum_k.block_structure.create_gf(
ish=self.icrsh,
gf_function=Gf,
space='solver',
mesh=MeshReTime(n_t=time_steps + 1, window=[0, time_steps * self.solver_params['dt']]),
)
# define abstract solve and postprocess methods
[docs]
@abstractmethod
def solve(self, **kwargs):
r"""
Solve the DMFT impurity problem
"""
pass
[docs]
@abstractmethod
def postprocess(self):
r"""
Postprocess the DMFT impurity problem
"""
pass
[docs]
@staticmethod
def _crm_dyson_solver(G0_dlr_iw, G_dlr, Sigma_moments, deg_shells):
"""
Self-energy of every block from the CRM Dyson solver
Each degenerate shell is minimized once, for its first block, and the
result is copied to the other blocks of the shell. Blocks that are in
no degenerate shell are minimized individually.
Parameters
----------
G0_dlr_iw : BlockGf on MeshDLRImFreq
non-interacting impurity Green's function
G_dlr : BlockGf on MeshDLR
interacting impurity Green's function
Sigma_moments : dict of np.ndarray
high-frequency moments of Sigma per block, ``[Sigma_0, (Sigma_1)]``
deg_shells : list
degenerate shells of block names, i.e. ``sum_k.deg_shells[ish]``
Returns
-------
Sigma_dlr : BlockGf on MeshDLRImFreq
self-energy without the constant ``Sigma_moments[block][0]``
"""
from triqs.gfs.dlr_crm_dyson_solver import minimize_dyson
source = {}
for shell in deg_shells:
first = next(iter(shell))
for block in shell:
source[block] = first
Sigma_dlr = G0_dlr_iw.copy()
Sigma_dlr.zero()
for block, _ in Sigma_dlr:
if source.get(block, block) == block:
np.random.seed(85281)
print('Minimizing Dyson via CRM for Σ[block]:', block)
Sigma_dlr[block], _, _ = minimize_dyson(
G0_dlr=G0_dlr_iw[block], G_dlr=G_dlr[block], Sigma_moments=Sigma_moments[block]
)
for block, _ in Sigma_dlr:
if source.get(block, block) != block:
print(f'Copying result from first block of deg list to {block}')
Sigma_dlr[block] << Sigma_dlr[source[block]]
return Sigma_dlr
[docs]
@staticmethod
def _gf_fit_tail_fraction(Gf, fraction=0.4, replace=None, known_moments=[]):
"""
fits the tail of Gf object by making a polynomial
fit of the Gf on the given fraction of the Gf mesh
and replacing that part of the Gf by the fit
0.4 fits the last 40% of the Gf and replaces the
part with the tail
Parameters
----------
Gf : BlockGf (Green's function) object
fraction: float, optional default 0.4
fraction of the Gf to fit
replace: float, optional default fraction
fraction of the Gf to replace
known_moments: np.array
known moments as numpy array
Returns
-------
Gf_fit : BlockGf (Green's function) object
fitted Gf
"""
Gf_fit = Gf.copy()
# if no replace factor is given use the same fraction
if not replace:
replace = fraction
for i, bl in enumerate(Gf_fit.indices):
Gf_fit[bl].mesh.set_tail_fit_parameters(tail_fraction=fraction)
if known_moments == []:
tail = Gf_fit[bl].fit_hermitian_tail()
else:
tail = Gf_fit[bl].fit_hermitian_tail(known_moments[i])
nmax_frac = int(len(Gf_fit[bl].mesh) / 2 * (1 - replace))
Gf_fit[bl].replace_by_tail(tail[0], n_min=nmax_frac)
return Gf_fit
[docs]
@staticmethod
def _fit_tail_window(
Sigma_iw, fit_min_n=None, fit_max_n=None, fit_min_w=None, fit_max_w=None, fit_max_moment=None, fit_known_moments=None
):
"""
Fit a high frequency 1/(iw)^n expansion of Sigma_iw
and replace the high frequency part with the fitted high frequency expansion.
Either give frequency window to fit on in terms of matsubara frequencies index
(fit_min_n/fit_max_n) or value (fit_min_w/fit_max_w).
Parameters
----------
Sigma_iw : Gf
Self-energy.
fit_min_n : int, optional, default=int(0.8*len(Sigma_iw.mesh))
Matsubara frequency index from which tail fitting should start.
fit_max_n : int, optional, default=int(len(Sigma_iw.mesh))
Matsubara frequency index at which tail fitting should end.
fit_min_w : float, optional
Matsubara frequency from which tail fitting should start.
fit_max_w : float, optional
Matsubara frequency at which tail fitting should end.
fit_max_moment : int, optional
Highest moment to fit in the tail of Sigma_iw.
fit_known_moments : ``ndarray.shape[order, Sigma_iw[0].target_shape]``, optional, default = None
Known moments of Sigma_iw, given as an numpy ndarray
Returns
-------
tail_barr : dict of arr
fitted tail of Sigma_iw
"""
from triqs.gfs.gf_fnt import fit_hermitian_tail_on_window, replace_by_tail
# Define default tail quantities
if fit_min_w is not None:
fit_min_n = int(0.5 * (fit_min_w * Sigma_iw.mesh.beta / np.pi - 1.0))
if fit_max_w is not None:
fit_max_n = int(0.5 * (fit_max_w * Sigma_iw.mesh.beta / np.pi - 1.0))
if fit_min_n is None:
fit_min_n = int(0.8 * len(Sigma_iw.mesh) / 2)
if fit_max_n is None:
fit_max_n = int(len(Sigma_iw.mesh) / 2)
if fit_max_moment is None:
fit_max_moment = 3
if fit_known_moments is None:
fit_known_moments = {}
for name, sig in Sigma_iw:
shape = [0] + list(sig.target_shape)
fit_known_moments[name] = np.zeros(shape, dtype=complex) # no known moments
# Now fit the tails of Sigma_iw and replace the high frequency part with the tail expansion
tail_barr = {}
Sigma_fit = Sigma_iw.copy()
for name, sig in Sigma_fit:
tail, err = fit_hermitian_tail_on_window(
sig,
n_min=fit_min_n,
n_max=fit_max_n,
known_moments=fit_known_moments[name],
# set max number of pts used in fit larger than mesh size, to use all data in fit
n_tail_max=10 * len(sig.mesh),
expansion_order=fit_max_moment,
)
tail_barr[name] = tail
replace_by_tail(sig, tail, n_min=fit_min_n)
return Sigma_fit, tail_barr
def _make_spin_equal(self, Sigma):
# if not SOC than average up and down
if not self.general_params['magnetic'] and not self.sum_k.SO == 1:
Sigma['up_0'] = 0.5 * (Sigma['up_0'] + Sigma['down_0'])
Sigma['down_0'] = Sigma['up_0']
return Sigma