Source code for triqs.sumk.sumk_discrete

# Copyright (c) 2013-2017 Commissariat à l'énergie atomique et aux énergies alternatives (CEA)
# Copyright (c) 2013-2017 Centre national de la recherche scientifique (CNRS)
# Copyright (c) 2017 Hugo Strand
# Copyright (c) 2020 Simons Foundation
#
# This program 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.
#
# This program 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 may obtain a copy of the License at
#     https:#www.gnu.org/licenses/gpl-3.0.txt
#
# Authors: Michel Ferrero, Olivier Parcollet, Hugo U. R. Strand, Nils Wentzell


from triqs.gf import *
import triqs.utility.mpi as mpi
from itertools import *
import inspect
import copy,numpy

[docs] class SumkDiscrete: """ INTERNAL USE The function to compute \[ G \leftarrow \sum_k (\omega + \mu - eps_k - Sigma(k,\omega))^{-1} \] for GF functions with blocks of the size of the matrix eps_k with a discrete sum. The class contains the discretized hoppings and points in the arrays hopping, bz_points,bz_weights,mu_pattern,overlap (IF non orthogonal) It can also generate a grid (ComputeGrid) for a regular grid or a Gauss-Legendre sum. """
[docs] def __init__ (self, dim, gf_struct, orthogonal_basis = True ): """ Just constructs the arrays, but without initializing them - dim is the dimension - gf_struct: Indices of the Green function - orthogonal_basis: True by default """ self.__GFBLOC_Structure = copy.deepcopy(gf_struct) self.orthogonal_basis,self.dim = orthogonal_basis,dim
#-------------------------------------------------------------
[docs] def resize_arrays (self, nk): """ Just constructs the arrays, but without initializing them - nk: total number of k points """ # constructs the arrays. no = len(self.__GFBLOC_Structure) self.hopping = numpy.zeros([nk,no,no],numpy.complex_) # t(k_index,a,b) self.bz_points = numpy.zeros([nk,self.dim],numpy.float_) # k(k_index,:) self.bz_weights = numpy.ones([nk],numpy.float_)/ float(nk) # w(k_kindex) , default normalisation self.mu_pattern = numpy.identity(no,numpy.complex_) if self.orthogonal_basis else numpy.zeros([no,no,nk],numpy.complex_) self.overlap = numpy.array(self.mu_pattern, copy=True)
#------------------------------------------------------------- def __get_GFBloc_Structure(self): """Returns the ONLY block indices accepted for the G and Sigma argument of the SumK function""" return self.__GFBLOC_Structure GFBlocIndices = property(__get_GFBloc_Structure) #------------------------------------------------------------- def __call__ (self, Sigma, mu=0, field=None, epsilon_hat=None, result=None, selected_blocks=()): """ - Computes: result <- \[ \sum_k (\omega + \mu - field - t(k) - Sigma(k,\omega)) \] if result is None, it returns a new GF with the results. otherwise, result must be a GF, in which the calculation is done, and which is then returned. (this allows chain calculation: SK(mu = mu,Sigma = Sigma, result = G).total_density() which computes the sumK into G, and returns the density of G. - Sigma can be a X, or a function k-> X or a function k,eps ->X where: - k is expected to be a 1d-numpy array of size self.dim of float, containing the k vector in the basis of the RBZ (i.e. -0.5< k_i <0.5) - eps is t(k) - X is anything such that X[BlockName] can be added/subtracted to a GFBloc for BlockName in selected_blocks. e.g. X can be a BlockGf(with at least the selected_blocks), or a dictionnary Blockname -> array if the array has the same dimension as the GF blocks (for example to add a static Sigma). Each block of X has to have the same shape as self.hopping or epsilon_hat(self.hopping[i]). - field: Any k independant object to be added to the GF - epsilon_hat: a function of eps_k returning a matrix with the same matrix-dimensions as each block in Sigma - selected_blocks: The calculation is done with the SAME t(k) for all blocks. If this list is not None only the blocks in this list are calculated. e.g. G and Sigma have block indices 'up' and 'down'. if selected_blocks ==None: 'up' and 'down' are calculated if selected_blocks == ['up']: only 'up' is calculated. 'down' is 0. """ assert selected_blocks == (), "selected_blocks not supported for now" #S = Sigma.view_selected_blocks(selected_blocks) if selected_blocks else Sigma #Gres = result if result else Sigma.copy() #G = Gres.view_selected_blocks(selected_blocks) if selected_blocks else Gres # check Sigma # case 1) Sigma is a BlockGf if isinstance(Sigma, BlockGf): model = Sigma Sigma_fnt = False # case 2) Sigma is a function returning a BlockGf else: assert callable(Sigma), "If Sigma is not a BlockGf it must be a function" Sigma_Nargs = len(inspect.getargspec(Sigma)[0]) assert Sigma_Nargs <= 2, "Sigma must be a function of k or of k and epsilon" if Sigma_Nargs == 1: model = Sigma(self.bz_points[0]) elif Sigma_Nargs == 2: model = Sigma(self.bz_points[0], self.hopping[0]) Sigma_fnt = True G = result if result else model.copy() assert isinstance(G,BlockGf), "G must be a BlockGf" assert isinstance(G.mesh, MeshImFreq), "G.mesh must be MeshImFreq but is {}".format(type(G.mesh)) # check input assert self.orthogonal_basis, "Local_G: must be orthogonal. non ortho cases not checked." # check that each block has the same size assert len(list(set([g.target_shape[0] for i,g in G]))) == 1 assert self.bz_weights.shape[0] == self.n_kpts(), "Internal Error" no = list(set([g.target_shape[0] for i,g in G]))[0] # check that the target shape of each block matches self.hopping eps_hat = epsilon_hat(self.hopping[0]) if epsilon_hat else self.hopping[0] assert (no,no) == eps_hat.shape, (f"Target shape of each block in Sigma: {(no,no)} does not to match orbital dimension of the hopping matrix: {eps_hat.shape}.") # Initialize G.zero() tmp,tmp2 = G.copy(),G.copy() mupat = mu * numpy.identity(no, numpy.complex_) tmp << iOmega_n if field != None: tmp -= field if not Sigma_fnt: tmp -= Sigma # substract Sigma once for all # Loop on k points... for w, k, eps_k in zip(*[mpi.slice_array(A) for A in [self.bz_weights, self.bz_points, self.hopping]]): eps_hat = epsilon_hat(eps_k) if epsilon_hat else eps_k tmp2 << tmp tmp2 -= tmp2.n_blocks * [eps_hat - mupat] if Sigma_fnt: if Sigma_Nargs == 1: tmp2 -= Sigma(k) elif Sigma_Nargs == 2: tmp2 -= Sigma(k,eps_k) tmp2.invert() tmp2 *= w G += tmp2 G << mpi.all_reduce(mpi.world,G,lambda x,y: x+y) mpi.barrier() return G #-------------------------------------------------------------
[docs] def n_kpts(self): """ Returns the number of k points""" return self.bz_points.shape[0]