Source code for adapol.hybfit

import numpy as np
from .fit_utils import pole_fitting, eval_with_pole

[docs] def hybfit( Delta, iwn_vec, tol=None, Np=None, svdtol=1e-7, solver="lstsq", maxiter=500, mmin=4, mmax=50, verbose=False, ): """ The function for hybridization fitting. Examples: ---------- - Fitting with :math:`N_p` poles: :code:`hybfit(Delta, iwn_vec, Np = Np)` - Fitting with fixed error tolerance tol : :code:`hybfit(Delta, iwn_vec, tol = tol)` Parameters: ------------ :code:`Delta`: np.array, :math:`(N_w, N_\mathrm{orb}, N_\mathrm{orb})` The input hybridization function in Matsubara frequency. :code:`iwn_vec`: np.array, :math:`(N_w)` The Matsubara frequency vector, complex-valued :code:`svdtol`: float, optional Truncation threshold for bath orbitals while doing SVD of weight matrices in hybridization fitting default:1e-7 :code:`tol`, :code:`Np`, :code:`solver`, :code:`maxiter`, :code:`mmin`, :code:`mmax`, :code:`verbose`: see below in anacont Returns: --------- :code:`bathenergy` :math:`E`: np.array, :math:`(N_b)` Bath energy :code:`bathhyb` :math:`V`: np.array, :math:`(N_b,N_{\mathrm{orb}})` Bath hybridization :code:`final_error`: float final fitting error :code:`func`: function Hybridization function evaluator :math:`f(z) = \sum_n V_{ni}V_{nj}^*/(z-E_n).` """ # Check dimensions assert len(iwn_vec.shape) == 1 or len(iwn_vec.shape) == 2 if len(iwn_vec.shape) == 2: assert iwn_vec.shape[1] == 1 iwn_vec = iwn_vec.flatten() assert len(Delta.shape) == 3 or len(Delta.shape) == 1 if len(Delta.shape) == 1: assert Delta.shape[0] == iwn_vec.shape[0] Delta = Delta[:, None, None] if len(Delta.shape) == 3: assert Delta.shape[0] == iwn_vec.shape[0] assert Delta.shape[1] == Delta.shape[2] solver = solver.lower() assert solver == "lstsq" or solver == "sdp" # Check input tol or Np if tol is None and Np is None: raise ValueError("Please specify either tol or Np") if tol is not None and Np is not None: raise ValueError( "Please specify either tol or Np. One can not specify both of them." ) wn_vec = np.imag(iwn_vec) if Np is not None: pol, weight, fitting_error = pole_fitting( Delta, wn_vec, Ns=Np + 1, maxiter=maxiter, solver=solver, disp=verbose ) elif tol is not None: pol, weight, fitting_error = pole_fitting( Delta, wn_vec, tol=tol, mmin=mmin, mmax=mmax, maxiter=maxiter, solver=solver, disp=verbose, ) bathenergy, bathhyb, bath_mat = obtain_orbitals(pol, weight, svdtol=svdtol) def func(Z): return eval_with_pole(bathenergy, Z, bath_mat) Delta_reconstruct = func(iwn_vec) final_error = np.max(np.abs(Delta - Delta_reconstruct)) return bathenergy, bathhyb, final_error, func
[docs] def hybfit_triqs( Delta_triqs, tol=None, Np=None, svdtol=1e-7, solver="lstsq", maxiter=500, mmin=4, mmax=50, verbose=False, debug=False ): """ The triqs interface for hybridization fitting. The function requires triqs package in python. Examples: ---------- - Fitting with :math:`N_p` poles: :code:`hybfit_triqs(delta_triqs, Np = Np)` - Fitting with fixed error tolerance tol : :code:`hybfit_triqs(delta_triqs, tol = tol)` Parameters: ------------ :code:`Delta_triqs`: triqs Green's function container The input hybridization function in Matsubara frequency :code:`debug`: bool return additional outputs for debugging. Default: False :code:`svdtol`: float, optional Truncation threshold for bath orbitals while doing SVD of weight matrices in hybridization fitting default: 1e-7 :code:`tol`, :code:`Np`, :code:`solver`, :code:`maxiter`, :code:`mmin`, :code:`mmax`, :code:`verbose`: same as in hybfit Returns: --------- :code:`bathhyb`: np.array :math:`(N_b, N_\mathrm{orb})` Bath hybridization :code:`bathenergy`: np.array :math:`(N_b,)` Bath energy :code:`Delta_fit`: triqs Gf or BlockGf Discretized hybridization function The input hybridization function in Matsubara frequency if debug is True: :code:`final_error`: float final fitting error :code:`weight`: np.array :math:`(N_p, N_\mathrm{orb}, N_\mathrm{orb})` weights obtained from fitting """ try: from triqs.gf import Gf, BlockGf, MeshImFreq, MeshDLRImFreq except ImportError: raise ImportError("Failed to import the triqs package (https://triqs.github.io/triqs/latest/). " "Please ensure it is installed.") if isinstance(Delta_triqs, Gf) and isinstance(Delta_triqs.mesh, (MeshImFreq, MeshDLRImFreq)): iwn_vec = np.array([iw.value for iw in Delta_triqs.mesh.values()]) eps_opt, V_opt, final_error, func = hybfit(Delta_triqs.data, iwn_vec, tol, Np, svdtol, solver, maxiter, mmin, mmax, verbose) print('optimization finished with fitting error {:.3e}'.format(final_error)) delta_fit = Gf(mesh=Delta_triqs.mesh, target_shape=Delta_triqs.target_shape) delta_fit.data[:] = func(iwn_vec) if debug: return V_opt.T.conj(), eps_opt, delta_fit, final_error else: return V_opt.T.conj(), eps_opt, delta_fit elif isinstance(Delta_triqs, BlockGf) and isinstance(Delta_triqs.mesh, (MeshImFreq, MeshDLRImFreq)): V_list, eps_list, delta_list, error_list = [], [], [], [] for j, (block, delta_blk) in enumerate(Delta_triqs): res = hybfit_triqs(delta_blk, tol, Np, svdtol, solver, maxiter, mmin, mmax, verbose, debug) V_list.append(res[0]) eps_list.append(res[1]) delta_list.append(res[2]) if debug: error_list.append(res[3]) if debug: return V_list, eps_list, BlockGf(name_list=list(Delta_triqs.indices), block_list=delta_list), error_list else: return V_list, eps_list, BlockGf(name_list=list(Delta_triqs.indices), block_list=delta_list) else: raise RuntimeError("Error: Delta_triqs.mesh must be an instance of MeshImFreq or MeshDLRImFreq.")
def obtain_orbitals(pol, weight, svdtol=1e-7): """ obtaining bath orbitals through svd """ polelist = [] veclist = [] matlist = [] for i in range(weight.shape[0]): eigval, eigvec = np.linalg.eig(weight[i]) for j in range(eigval.shape[0]): if eigval[j] > svdtol: polelist.append(pol[i]) veclist.append(eigvec[:, j] * np.sqrt(eigval[j])) matlist.append( (eigvec[:, j, None] * np.conjugate(eigvec[:, j].T)) * (eigval[j]) ) return np.array(polelist), np.array(veclist), np.array(matlist)