Source code for pytriqs.gf.tools

################################################################################
#
# TRIQS: a Toolbox for Research in Interacting Quantum Systems
#
# Copyright (C) 2011 by M. Ferrero, O. Parcollet
# Copyright (C) 2019 The Simons Foundation
# Authors: H. U.R. Strand
#
# TRIQS 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.
#
# TRIQS 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
# TRIQS. If not, see <http://www.gnu.org/licenses/>.
#
################################################################################
import lazy_expressions, descriptors, gf_fnt
from meshes import MeshImFreq, MeshReFreq, MeshImTime, MeshReTime, MeshLegendre
from block_gf import BlockGf
from gf import Gf
import numpy as np
from itertools import product
from backwd_compat.gf_refreq import GfReFreq 
from map_block import map_block

[docs]def inverse(x): """ Return the inverse of a Green's function """ if descriptors.is_lazy(x): return lazy_expressions.lazy_function("inverse", inverse) (x) assert hasattr(x,'inverse') return x.inverse()
[docs]def conjugate(x): """ Return the conjugate of a Green's function """ if descriptors.is_lazy(x): return lazy_expressions.lazy_function("conjugate", conjugate) (x) assert hasattr(x,'conjugate') return x.conjugate()
[docs]def transpose(x): """ Return the transpose of a Green's function """ if descriptors.is_lazy(x): return lazy_expressions.lazy_function("transpose", transpose) (x) assert hasattr(x,'transpose') return x.transpose()
[docs]def delta(g): """ Compute Delta_iw from G0_iw. CAUTION: This function assumes the following properties of g * The diagonal components of g should decay as 1/iOmega * g should fullfill the property g[iw][i,j] = conj(g[-iw][j,i]) Parameters ---------- g : BlockGf (of GfImFreq) or GfImFreq Non-interacting Green's function. Returns ------- delta_iw : BlockGf (of GfImFreq) or GfImFreq Hybridization function. """ if isinstance(g, BlockGf): return BlockGf(name_block_generator = [ (n, delta(g0)) for n,g0 in g], make_copies=False) elif isinstance(g.mesh, MeshImFreq): assert len(g.target_shape) in [0,2], "delta(g) requires a matrix or scalar_valued Green function" assert gf_fnt.is_gf_hermitian(g), "delta(g) requires a Green function with the property g[iw][i,j] = conj(g[-iw][j,i])" delta_iw = g.copy() delta_iw << descriptors.iOmega_n - inverse(g) tail, err = gf_fnt.fit_hermitian_tail(delta_iw) delta_iw << delta_iw - tail[0] if err > 1e-5: print "WARNING: delta extraction encountered a sizeable tail-fit error: ", err return delta_iw else: raise TypeError, "No function delta for g0 object of type %s"%type(g)
# Determine one of G0_iw, G_iw and Sigma_iw from other two using Dyson's equation
[docs]def dyson(**kwargs): """ Solve Dyson's equation for given two of G0_iw, G_iw and Sigma_iw to yield the third. Parameters ---------- G0_iw : Gf, optional Non-interacting Green's function. G_iw : Gf, optional Interacting Green's function. Sigma_iw : Gf, optional Self-energy. """ if not (len(kwargs)==2 and set(kwargs.keys())<set(['G0_iw','G_iw', 'Sigma_iw'])): raise ValueError, 'dyson: Two (and only two) of G0_iw, G_iw and Sigma_iw must be provided to determine the third.' if 'G0_iw' not in kwargs: G0_iw = inverse(kwargs['Sigma_iw'] + inverse(kwargs['G_iw'])) return G0_iw elif 'G_iw' not in kwargs: G_iw = inverse(inverse(kwargs['G0_iw']) - kwargs['Sigma_iw']) return G_iw elif 'Sigma_iw' not in kwargs: Sigma_iw = inverse(kwargs['G0_iw']) - inverse(kwargs['G_iw']) return Sigma_iw
[docs]def read_gf_from_txt(block_txtfiles, block_name): """ Read a GfReFreq from text files with the format (w, Re(G), Im(G)) for a single block. Notes ----- A BlockGf must be constructed from multiple GfReFreq objects if desired. The mesh must be the same for all files read in. Non-uniform meshes are not supported. Parameters ---------- block_txtfiles: Rank 2 square np.array(str) or list[list[str]] The text files containing the GF data that need to read for the block. e.g. [['up_eg1.dat']] for a one-dimensional block and [['up_eg1_1.dat','up_eg2_1.dat'], ['up_eg1_2.dat','up_eg2_2.dat']] for a 2x2 block. block_name: str Name of the block. Returns ------- g: GfReFreq The real frequency Green's function read in. """ block_txtfiles = np.array(block_txtfiles) # Must be an array to use certain functions N1,N2 = block_txtfiles.shape mesh = np.genfromtxt(block_txtfiles[0,0],usecols=[0]) # Mesh needs to be the same for all blocks g = GfReFreq(indices=range(N1),window=(np.min(mesh),np.max(mesh)),n_points=len(mesh),name=block_name) for i,j in product(range(N1),range(N2)): data = np.genfromtxt(block_txtfiles[i,j],usecols=[1,2]) g.data[:,i,j] = data[:,0] + 1j*data[:,1] return g
[docs]def write_gf_to_txt(g): """ Write a GfReFreq or GfImFreq to in the format (w/iw, Re(G), Im(G)) for a single block. Parameters ---------- g: GfReFreq or GfImFreq The real/imaginary frequency Green's function to be written out. """ if isinstance(g.mesh, MeshReFreq): mesh = np.array([w.real for w in g.mesh]).reshape(-1,1) elif isinstance(g.mesh, MeshImFreq): mesh = np.array([w.imag for w in g.mesh]).reshape(-1,1) else: raise ValueError, 'write_gf_to_txt: Only GfReFreq and GfImFreq quantities are supported.' for i,j in product(range(g.target_shape[0]),range(g.target_shape[1])): txtfile = '%s_%s_%s.dat'%(g.name,i,j) redata = g.data[:,i,j].real.reshape((g.data.shape[0],-1)) imdata = g.data[:,i,j].imag.reshape((g.data.shape[0],-1)) mesh_and_data = np.hstack((mesh,redata,imdata)) np.savetxt(txtfile,mesh_and_data)
[docs]def make_zero_tail(g, n_moments=10): """ Return a container for the high-frequency coefficients of a given Green function initialized to zero. Parameters ---------- g: GfImFreq or GfReFreq or GfImTime or GfReTime The real/imaginary frequency/time Green's function that we create the tail-array for. n_moments [default=10]: The number of high-frequency moments in the tail (starting from order 0). """ if isinstance(g, Gf) and isinstance(g.mesh, (MeshImFreq, MeshReFreq, MeshImTime, MeshReTime)): n_moments = max(1, n_moments) return np.zeros((n_moments,) + g.target_shape, dtype = np.complex128) elif isinstance(g, BlockGf): return map_block(lambda g_bl: make_zero_tail(g_bl, n_moments), g) else: raise RuntimeError, "Error: make_zero_tail has to be called on a frequency or time Green function object"
[docs]def fit_legendre(g_t, order=10): """ General fit of a noisy imaginary time Green's function to a low order Legendre expansion in imaginary time. Only Hermiticity is imposed on the fit, so discontinuities has to be fixed separately (see the method enforce_discontinuity) Author: Hugo U.R. Strand Parameters ---------- g_t : TRIQS imaginary time Green's function (matrix valued) Imaginary time Green's function to fit (possibly noisy binned data) order : int Maximal order of the fitted Legendre expansion Returns ------- g_l : TRIQS Legendre polynomial Green's function (matrix valued) Fitted Legendre Green's function with order `order` """ import numpy.polynomial.legendre as leg if isinstance(g_t, BlockGf): return map_block(lambda g_bl: fit_legendre(g_bl, order), g_t) assert isinstance(g_t, Gf) and isinstance(g_t.mesh, MeshImTime), "fit_legendre expects imaginary-time Green function objects" assert len(g_t.target_shape) == 2, "fit_legendre currently only implemented for matrix_valued Green functions" # -- flatten the data to 2D N_tau x (N_orb * N_orb) shape = g_t.data.shape fshape = [shape[0], np.prod(shape[1:])] # -- extend data accounting for hermiticity mesh = g_t.mesh tau = np.array([ t.value for t in mesh ]) # Rescale to the interval (-1,1) x = 2. * tau / mesh.beta - 1. data = g_t.data.reshape(fshape) data_herm = np.transpose(g_t.data, axes=(0, 2, 1)).conjugate().reshape(fshape) # -- Separated real valued linear system, with twice the number of RHS terms data_re = 0.5 * (data + data_herm).real data_im = 0.5 * (data + data_herm).imag data_ext = np.hstack((data_re, data_im)) c_l_ext = leg.legfit(x, data_ext, order - 1) c_l_re, c_l_im = np.split(c_l_ext, 2, axis=-1) c_l = c_l_re + 1.j * c_l_im # -- make Legendre Green's function of the fitted coeffs lmesh = MeshLegendre(mesh.beta, mesh.statistic, order) # Nb! We have to scale the actual Legendre coeffs to the Triqs "scaled" Legendre coeffs # see Boehnke et al. PRB (2011) l = np.arange(len(lmesh)) scale = np.sqrt(2.*l + 1) / mesh.beta scale = scale.reshape([len(lmesh)] + [1]*len(g_t.target_shape)) g_l = Gf(mesh=lmesh, target_shape=g_t.target_shape) g_l.data[:] = c_l.reshape(g_l.data.shape) / scale return g_l