DMFT for Hubbard and Hubbard-like models
We now present a Python script for performing DMFT calculations with the NRG
as the impurity solver. The script is fairly complete and robust, and as such it is suitable for production runs.
It supports a variety of models with different symmetry types (including cases
with block structure and with matrix-valued Green’s functions), controls
the chemical potential to achieve the selected band occupancy, performs linear or Broyden mixing
to accelerate the convergence, stores intermediate results (checkpoints) to enable
restarts, and supports arbitrary density of states in the band (e.g. tabulated in a file).
Here is the listing:
# DMFT(NRG) for Hubbard and Hubbard-like models (including Hubbard-Holstein model)
# RZ, Feb 2020
from triqs.gf import *
from triqs.operators import *
from h5 import *
from triqs.utility import mpi
from nrgljubljana_interface import Solver, MeshReFreqPts, hilbert_transform_refreq
import math, os, warnings
import numpy as np
from scipy import interpolate, integrate, special, optimize
from collections import OrderedDict
# Parameters
model = "SIAM" # impurity model (from the template library)
dos = "Bethe" # "Bethe" for a Bethe lattice (semicircular) DOS or a file name for reading tabulated DOS data from a file;
# format of the DOS file: lines with "eps rho(eps)" pairs; arbitrary bandwidth, but should be within the NRG mesh
Bethe_unit = 1.0 # half-bandwidth in the Bethe lattice DOS for dos="Bethe"; 1.0 for D=1 units, 2.0 for t=1 units
U = 2.0 # Hubbard electron-electron repulsion
occupancy_goal = 0.8 # band occupancy (summed over spin)
T = 1e-4 # temperature
B = None # magnetic field for spin-polarized calculation with U(1) symmetry;
# None = non-spin-polarized calculation with SU(2) symmetry
omega = None # Holstein phonon frequency
g1 = None # electron-phonon coupling strength
n1 = None # shift n1 in the Holstein model definition
verbose = True # show info messages during the iteration
verbose_nrg = False # show detailed output from the NRG solver
store_steps = True # store intermediate results to files (one subdirectory per iteration)
min_iter = 5 # minimum number of iterations (prevents premature exit from the DMFT loop)
max_iter = 50 # maximum number of iterations (signals convergence failure)
eps_prev = 1e-5 # convergence criterium: integrated diff between two consecutive local spectral functions
eps_loc_imp = 1e-5 # convergence criterium: integrated diff between local and impurity spectral function
eps_occupancy = 1e-4 # convergence criterium: difference from the target occupancy (occupancy_goal)
mixing_method = "broyden" # "linear" or "broyden" mixing; the quantity being mixed is the hybridisation function
occup_method = "adjust" # "adjust" or None; adjust=shift mu so that the current GF meets the occupancy goal
alpha = 0.5 # mixing parameter from both linear and Broyden mixing
max_mu_adjust = 10 # number of cycles for adjusting the value of the chemical potential
mesh_max = 10.0 # logarithmic mesh: maximum frequency (should be large enough to contain the full spectrum)
mesh_min = 1e-5 # logarithmic mesh: minimum frequency (should be smaller than the temperature T)
mesh_ratio = 1.01 # logarithmic mesh: common ratio of the mesh (1.01 is usually a good choice)
Delta_min = 1e-5 # minimal value for the hybridisation function; if too large, it produces incorrect spectral distribution,
# if too small, it leads to discretization problems (1e-5 is usually a good and rather safe choice)
normalize_to_one = True # scaling of the spectral function (relevant in the case of block structure and/or matrix structure of GF)
solution_filename = "solution.h5"
checkpoint_filename = "checkpoint.h5"
# Additional quantities of interest
observables = ["n_d", "n_d^2"]
if B is not None:
observables.extend(["SZd"])
if "Holstein" in model:
observables.extend(["nph", "displ", "displ^2"])
# Run-time global variables
itern = 0 # iteration counter
verbose = verbose and mpi.is_master_node() # output is only produced by the master node
store_steps = store_steps and mpi.is_master_node() # output is only produced by the master node
symtype = ("QS" if B is None else "QSZ")
# Set up the Solver
S = Solver(model=model, symtype=symtype, mesh_max=mesh_max, mesh_min=mesh_min, mesh_ratio=mesh_ratio)
S.set_verbosity(verbose_nrg)
newG = lambda : S.G_w.copy() # Creates a BlockGf object of appropriate structure
nr_blocks = lambda bgf : len([bl for bl in bgf.indices]) # Returns the number of blocks in a BlockGf object
block_size = lambda bl : len(S.G_w[bl].indices[0]) # Matrix size of Green's functions in block 'bl'
identity = lambda bl : np.identity(block_size(bl)) # Returns the identity matrix in block 'bl'
# Solve Parameters
sp = { "T": T, "Lambda": 2.0, "Nz": 4, "Tmin": 1e-5, "keep": 10000, "keepenergy": 10.0 }
# Model Parameters
mp = { "U1": U, "B1": B, "omega": omega, "g1": g1, "n1": n1 }
mp = { k: v for k,v in mp.items() if v is not None }
sp["model_parameters"] = mp
# Update the chemical potential
def set_mu(new_mu):
global mu
mu = new_mu
sp["model_parameters"]["eps1"] = -mu
# Low-level NRG Parameters (optional tweaks)
nrgp = {}
#nrgp["bandrescale"] = 5.0
S.set_nrg_params(**nrgp)
# Initialize a function ht0 for calculating the Hilbert transform of the DOS
if (dos == "Bethe"):
ht1 = lambda z: 2*(z-1j*np.sign(z.imag)*np.sqrt(1-z**2)) # Analytical expression for Hilbert transform of Bethe lattice DOS
global ht0
ht0 = lambda z: ht1(z/Bethe_unit)
else:
table = np.loadtxt(dos)
global dosA
dosA = Gf(mesh=MeshReFreqPts(table[:,0]), target_shape=[])
for i, w in enumerate(dosA.mesh):
dosA[w] = np.array([[ table[i,1] ]])
ht0 = lambda z: hilbert_transform_refreq(dosA, z)
# Wrapper for ht0 that ensures that the imaginary part of Sigma is negative
EPS = 1e-20 # Cut-off parameter for fixing causality violation
ht = lambda z: ht0(z.real+1j*(z.imag if z.imag>0.0 else EPS)) # Fix problems due to causality violation
# Calculate a GF from hybridisation and self-energy
def calc_G(Delta, Sigma, mu):
G = newG()
for bl in G.indices:
for w in G.mesh:
G[bl][w] = np.linalg.inv( (w + mu)*identity(bl) - Delta[bl][w] - Sigma[bl][w] ) # !!!
return G
# Index range of a GF
index_range = lambda G : range(len(G.indices[0]))
# Return an interpolation-object representation of a spectral function for GF G
def interp_A(G):
lx = np.array(list(G.mesh.values()))
ly = sum( sum( -(1.0/math.pi)*np.array(G[bl].data[:,i,i].imag) for i in index_range(G[bl]) ) # matrix trace
for bl in G.indices ) # sum over blocks
if normalize_to_one:
nr = sum( sum( 1 for i in index_range(G[bl]) ) for bl in G.indices ) # number of contributions
ly = ly/nr # rescale
return interpolate.interp1d(lx, ly, kind='cubic', bounds_error=False, fill_value=0)
# Calculate occupancy for given hybridisation, self-energy and chemical potential
def calc_occupancy(Delta, Sigma, mu):
Gtrial = calc_G(Delta, Sigma, mu)
f = interp_A(Gtrial)
with warnings.catch_warnings():
warnings.simplefilter("ignore")
n = integrate.quad(lambda x : 2*f(x)*special.expit(-x/T), -mesh_max, mesh_max)
return n[0]
# Update mu towards reaching the occupancy goal
def update_mu(Delta, Sigma):
sol = optimize.root_scalar(lambda x : calc_occupancy(Delta, Sigma, x)-occupancy_goal, x0=mu, x1=mu-0.1)
set_mu(sol.root)
# Calculate the local lattice GF and the hybridisation function for the effective impurity model
# for a given self-energy function
def self_consistency(Sigma):
Gloc = newG()
for bl in Gloc.indices:
for w in Gloc.mesh:
for i in range(block_size(bl)):
for j in range(block_size(bl)): # assuming square matrix
if i == j:
Gloc[bl][w][i,i] = ht(w + mu - Sigma[bl][w][i,i]) # Hilbert-transform
else:
assert abs(Sigma[bl][w][i,j])<1e-10, "This implementation only supports diagonal self-energy"
Gloc[bl][w][i,j] = 0.0
Glocinv = Gloc.inverse()
Delta = newG()
for bl in Delta.indices:
for w in Delta.mesh:
Delta[bl][w] = (w+mu)*identity(bl) - Sigma[bl][w] - Glocinv[bl][w] # !!!
return Gloc, Delta
# Iteratively adjust mu, taking into account the self-consistency equation.
# Returns an improved estimate for the hybridisation function.
def adjust_mu(Delta_in, Sigma):
old_mu = mu
Delta = Delta_in.copy()
for _ in range(max_mu_adjust):
update_mu(Delta, Sigma)
Gloc, Delta = self_consistency(Sigma)
new_mu = mu
if verbose: print("Adjusted mu from %f to %f\n" % (old_mu,new_mu))
return Gloc, Delta
# Difference between two Green's functions evaluated as the integrated squared difference between the
# corresponding spectral functions.
def gf_diff(a, b):
f_a = interp_A(a)
f_b = interp_A(b)
with warnings.catch_warnings():
warnings.simplefilter("ignore")
diff = integrate.quad(lambda x : (f_a(x)-f_b(x))**2, -mesh_max, mesh_max)
return diff[0]
# Save a Green's function to a tabulated ASCII file
def save_Gf(fn, gf):
f = open(fn, "w")
for w in gf.mesh:
z = gf[w]
f.write("%f %f %f\n" % (w, z.real, z.imag))
# Save all blocks for a block GF to tabulated ASCII files
def save_BlockGf(fn, bgf):
for bl in bgf.indices:
save_Gf(fn + "_" + bl + ".dat", bgf[bl])
# Save a spectral function (-1/Pi Im GF) to a tabulated ASCII file
def save_A(fn, gf):
f = open(fn, "w")
for w in gf.mesh:
z = gf[w]
f.write("%f %f\n" % (w, -1.0/math.pi*z.imag))
# Save spectral functions for all blocks of the block GF
def save_BlockA(fn, bgf):
for bl in bgf.indices:
save_A(fn + "_" + bl + ".dat", bgf[bl])
# Store the complete set of results
def store_result(fn, S):
with HDFArchive(fn, 'w') as arch:
arch["S"] = S
# Global variables
arch["Gself"] = Gself
arch["Gloc"] = Gloc
arch["mu"] = mu
# Load the minimal set of stored results for restarting the calculation
def load_Sigma_mu(fn):
with HDFArchive(fn, 'r') as arch:
Sigma = arch["S"].Sigma_w
mu = arch["mu"]
return Sigma, mu
# Initial Sigma and mu from a file
def restart_calculation(fn):
if verbose: print("Starting from stored results in file %s\n" % fn)
return load_Sigma_mu(fn) # Load data from an HDF5 archive
# Initial Sigma and mu when starting from scratch
def new_calculation():
if verbose: print("Starting from scratch\n")
Sigma = newG()
for bl in Sigma.indices:
for w in Sigma.mesh:
Sigma[bl][w] = U*occupancy_goal/2.0 # Initialize self-energy with the Hartree shift
mu = U/2.0 # initial approximaiton for the chemical potential
return Sigma, mu
# Provides the initial Sigma and mu
def initial_Sigma_mu():
if os.path.isfile(solution_filename): # continue DMFT iteration after failed convergence
return restart_calculation(solution_filename)
elif os.path.isfile(checkpoint_filename): # continue DMFT iteration after interruption
return restart_calculation(checkpoint_filename)
else: # start from scratch
return new_calculation()
# Exception to raise when convergence is reached
class Converged(Exception):
def __init__(self, message):
self.message = message
# Exception to raise when convergence is not reached (e.g. maximum nr of iterations exceeded)
class FailedToConverge(Exception):
def __init__(self, message):
self.message = message
# Formatting of the header in stats.dat
def fmt_str_header(nr_val):
str = "{:>5}" # itern
for _ in range(nr_val-1): str += " {:>15}"
return str + "\n"
# Formatting of the results in stats.dat
def fmt_str(nr_val):
str = "{:>5}" # itern
for _ in range(nr_val-1): str += " {:>15.8g}"
return str + "\n"
# Adjust Im(Delta) so that the hybridisation strength is not too small for the NRG discretization
# scheme to break down.
def fix_hyb_function(Delta, Delta_min):
Delta_fixed = Delta.copy()
for bl in Delta.indices:
for w in Delta.mesh:
for n in range(block_size(bl)): # only diagonal parts
r = Delta[bl][w][n,n].real
i = Delta[bl][w][n,n].imag
Delta_fixed[bl][w][n,n] = r + 1j*(i if i<-Delta_min else -Delta_min)
# Possible improvement: re-adjust the real part so that the Kramers-Kronig relation is maintained
return Delta_fixed
# Perform a DMFT step. Input is the hybridization function for solving the effective impurity model,
# output is a new hybridization function resulting from the application of the DMFT self-consistency equation.
def dmft_step(Delta_in):
global itern
itern += 1
if verbose:
print("Iteration %i min_iter=%i max_iter=%i\n" % (itern, min_iter, max_iter))
Delta_in_fixed = fix_hyb_function(Delta_in, Delta_min)
S.Delta_w << Delta_in_fixed
S.solve(**sp) # Solve the impurity model
global Gself, Gloc, Gloc_prev
Gself = calc_G(Delta_in_fixed, S.Sigma_w, mu) # impurity GF ("self-energy-trick" improved)
Gloc, Delta = self_consistency(S.Sigma_w) # apply the DMFT self-consistency equation
diff_loc_imp = gf_diff(Gself, Gloc) # difference between impurity and local lattice GF
diff_prev = gf_diff(Gloc, Gloc_prev) # difference between two consecutively computed local latice GFs
Gloc_prev = Gloc.copy()
occupancy = calc_occupancy(Delta, S.Sigma_w, mu)
diff_occupancy = abs(occupancy-occupancy_goal) # this difference is used as the measure of deviation
stats = OrderedDict([("itern", itern), ("mu", mu), ("diff_loc_imp", diff_loc_imp), ("diff_prev", diff_prev),
("diff_occupancy", diff_occupancy), ("occupancy", occupancy)])
for i in observables:
stats[i] = S.expv[i]
header_string = fmt_str_header(len(stats)).format(*[i for i in stats.keys()])
stats_string = fmt_str(len(stats)).format(*[i for i in stats.values()])
if mpi.is_master_node():
if itern == 1: stats_file.write(header_string)
stats_file.write(stats_string)
if verbose: print("stats: %sstats: %s" % (header_string, stats_string))
if store_steps:
os.mkdir(str(itern)) # one subdirectory per iteration
save_BlockGf(str(itern)+"/Delta", Delta_in_fixed)
save_BlockGf(str(itern)+"/Sigma", S.Sigma_w) # self-energy
save_BlockGf(str(itern)+"/G", Gloc) # local lattice Green's function
save_BlockA(str(itern)+"/A", Gloc) # spectral function of local lattice GF
store_result(str(itern)+"/"+solution_filename, S)
if mpi.is_master_node():
store_result(checkpoint_filename, S) # for checkpoint/restart functionality
# Check for convergence. The only way to exit the DMFT loop is by generating exceptions.
if (diff_loc_imp < eps_loc_imp and
diff_prev < eps_prev and
diff_occupancy < eps_occupancy and
itern >= min_iter):
raise Converged(stats_string)
if (itern == max_iter):
raise FailedToConverge(stats_string)
if occup_method == "adjust":
Gloc, Delta = adjust_mu(Delta, S.Sigma_w) # here we update mu to get closer to target occupancy
return Delta
# DMFT driver routine with linear mixing
def solve_with_linear_mixing(Delta, alpha):
Delta_in = Delta.copy()
while True:
Delta_out = dmft_step(Delta_in)
newDelta = alpha*Delta_out + (1-alpha)*Delta_in
Delta_in << newDelta
# Convert GF object to a linear nparray
def gf_to_nparray(gf):
return gf.data.flatten()
# Stack elements from all blocks in a block BF
def bgf_to_nparray(bgf):
return np.hstack((bgf[bl].data.flatten() for bl in bgf.indices))
# Convert a linear numpy array to a single GF object
def nparray_to_gf(a, gf):
b = a.reshape(gf.data.shape)
gf.data[:,:,:] = b[:,:,:]
# Extract blocks from linear numpy array to a block GF
def nparray_to_bgf(a):
G = newG()
split = np.split(a, nr_blocks(G)) # Here we assume all blocks are equally large
for i, bl in enumerate(G.indices):
nparray_to_gf(split[i], G[bl])
return G
# DMFT driver routine with Broyden mixing: the goal is to find a root of the function F(Delta)=dmft_step(Delta)-Delta.
F = lambda Delta : dmft_step(Delta)-Delta
npF = lambda x : bgf_to_nparray(F(nparray_to_bgf(x)))
def solve_with_broyden_mixing(Delta, alpha):
xin = bgf_to_nparray(Delta)
optimize.broyden1(npF, xin, alpha=alpha, reduction_method="svd", max_rank=10, verbose=verbose, f_tol=1e-99) # Loop forever (small f_tol!)
def solve(Delta):
if (mixing_method=="linear"):
solve_with_linear_mixing(Delta, alpha)
if (mixing_method=="broyden"):
solve_with_broyden_mixing(Delta, alpha)
# Initialization
Sigma0, mu0 = initial_Sigma_mu() # Get initial self-energy and chemical potential
set_mu(mu0) # Set global variable mu and update the impurity level (epsilon_d = -mu)
Gloc, Delta = self_consistency(Sigma0) # Initialize local GF Gloc and hybridization function Delta
Gloc, Delta = adjust_mu(Delta, Sigma0) # Adjust mu (and return updated Glocal and Delta as well)
Gloc_prev = Gloc.copy() # Store copy for checking convergence
if mpi.is_master_node():
global stats_file
stats_file = open("stats.dat", "w", buffering=1) # line buffered
try:
solve(Delta)
except Converged as c:
if verbose: print("Converged: %s" % c.message)
if mpi.is_master_node():
store_result(solution_filename, S) # full converged results as an HDF5 file
os.remove(checkpoint_filename) # checkpoint file is no longer needed
save_BlockA("A", Gloc) # converged spectral function for quick plotting
except FailedToConverge as c:
if verbose: print("Failed to converge: %s" % c.message) # ... but restart is possible from the checkpoint file
mpi.barrier() # Synchronized exit
Running this script takes between an hour or several hours, depending on the
degree of parallelization that is enabled and the available computing resources.
During the run, some basic information about the DMFT iteration is stored in the file stats.dat:
itern mu diff_loc_imp diff_prev diff_occupancy occupancy n_d n_d^2
1 0.64216213 0.008410854 0.1639845 0.065470795 0.86547079 0.88853941 1.0170551
2 0.57416456 0.0030978054 0.0033493961 0.048255692 0.84825569 0.86718719 0.97416922
3 0.52350142 0.00032006683 0.0045037239 0.044544862 0.84454486 0.85230857 0.94144344
4 0.47628807 0.00017268295 0.00071735455 0.027302305 0.82730231 0.8329584 0.9156406
5 0.44771388 5.3787042e-05 0.0006119636 0.015770286 0.81577029 0.81874584 0.89894272
6 0.4314115 3.43912e-05 9.5803891e-05 0.0092612995 0.8092613 0.81131439 0.88949259
7 0.42188467 3.6268579e-05 7.0325948e-05 0.0066369482 0.80663695 0.80579098 0.88292743
8 0.41508851 4.5411105e-05 1.8109524e-05 0.0033722505 0.80337225 0.80338313 0.87950144
9 0.41163793 5.8842502e-05 3.9934759e-06 0.0018214352 0.80182144 0.80201381 0.87764421
10 0.4097753 6.3801883e-05 1.2231238e-06 0.0010891672 0.80108917 0.80114589 0.8765355
11 0.40866208 6.6673362e-05 4.6656617e-07 0.00065926119 0.80065926 0.80061656 0.87586398
12 0.40798848 6.8445133e-05 1.7684942e-07 0.00039989342 0.80039989 0.80029503 0.87545665
13 0.40757997 6.9588067e-05 6.2929404e-08 0.00024353334 0.80024353 0.8000997 0.87520934
14 0.40733122 6.6379215e-05 1.9177454e-05 0.0011155475 0.80111555 0.79886188 0.87412851
15 0.40619387 4.4607833e-05 3.3092449e-05 0.0015344904 0.79846551 0.79964663 0.87487525
16 0.40776042 9.3749908e-07 2.8603095e-05 0.00053041152 0.79946959 0.79906779 0.87497248
17 0.40830088 2.4812637e-07 1.0216056e-07 0.00028239426 0.79971761 0.79922373 0.87523217
18 0.40858862 1.3398751e-07 1.5633394e-07 0.00031811496 0.79968189 0.79939042 0.87547037
19 0.40891278 4.8884804e-09 6.2918295e-08 0.00011984463 0.79988016 0.79949007 0.87561332
20 0.4090349 1.3450417e-09 1.0571106e-08 6.0784735e-05 0.79993922 0.79955157 0.87568519
The first column is the iteration number, the second the current estimate of the chemical potential, the following three columns provide information about the convergence (difference between local lattice and impurity Green’s functions, difference between two consecutive local lattice Green’s function, and the difference between target and current band occupancy). This is followed by the occupancy calculed from the spectral function and by the expectation values of a range of local operators (the set depends on the model considered and on the symmetry type).
Note: NRG impurity solver produces big temporary files during the calculation. By default, these are stored in the working directory that is created in the current directory. To override this behavior, set the environment variable NRG_WORKDIR. A suitable choice is a local scratch hard drive or SSD, or a global scratch mount on a high-performance filesystem.
When convergence is reached, the final results are stored in an HDF5 archive file called solution.h5.
This file contains all quantities of interest. Let us plot the spectral function:
import numpy as np
import matplotlib as mpl
import matplotlib.pyplot as plt
import math
from triqs.gf import *
from h5 import *
from nrgljubljana_interface import MeshReFreqPts
def GF_to_nparrays(A):
lx = np.array(list(A.mesh.values()))
ly = np.array(-1.0/math.pi*A.data[:,0,0].imag)
return lx, ly
with HDFArchive('solution.h5','r') as ar:
# Expectation values
print("<n>=",ar['expv']['n_d'])
print("<n^2>=",ar['expv']['n_d^2'])
# Spectral function (of local lattice GF)
A_w = ar['Gloc']['imp']
lx, ly = GF_to_nparrays(A_w)
plt.plot(lx, ly)
plt.show()
Let us now go through the script in some more detail. We first import a number of useful modules:
from triqs.gf import *
from triqs.operators import *
from h5 import *
from triqs.utility import mpi
from nrgljubljana_interface import Solver, MeshReFreqPts, hilbert_transform_refreq
import math, os, warnings
import numpy as np
from scipy import interpolate, integrate, special, optimize
from collections import OrderedDict
We make use of scipy interpolation and integration routines to obtain highly accurate results based on
the real-frequency spectra, and we make use of the Broyden solver from scipy.optimize.
The model name is a full path to the model definition in the template library, except the symmetry part (e.g. “SIAM”, “Holstein/Nph=10”, etc.).
model = "SIAM" # impurity model (from the template library)
The code comes with a built-in support for the Bethe lattice (with half-bandwidth set by the parameter Bethe_unit),
or one can use tabulated DOS (example script for generating the DOS file).
dos = "Bethe" # "Bethe" for a Bethe lattice (semicircular) DOS or a file name for reading tabulated DOS data from a file;
# format of the DOS file: lines with "eps rho(eps)" pairs; arbitrary bandwidth, but should be within the NRG mesh
Bethe_unit = 1.0 # half-bandwidth in the Bethe lattice DOS for dos="Bethe"; 1.0 for D=1 units, 2.0 for t=1 units
We then set key parameters, such as Hubbard \(U\), the occupancy, the temperature, and optionally the magnetic field.
The parameters that are not used should be set to None.
U = 2.0 # Hubbard electron-electron repulsion
occupancy_goal = 0.8 # band occupancy (summed over spin)
T = 1e-4 # temperature
B = None # magnetic field for spin-polarized calculation with U(1) symmetry;
# None = non-spin-polarized calculation with SU(2) symmetry
One can control the level of verbosity and the amount of intermediate results that are stored to disk:
verbose = True # show info messages during the iteration
verbose_nrg = False # show detailed output from the NRG solver
store_steps = True # store intermediate results to files (one subdirectory per iteration)
The DMFT iteration parameters are set next. Note that there are several convergence criteria: they must all be satisfied.
min_iter = 5 # minimum number of iterations (prevents premature exit from the DMFT loop)
max_iter = 50 # maximum number of iterations (signals convergence failure)
eps_prev = 1e-5 # convergence criterium: integrated diff between two consecutive local spectral functions
eps_loc_imp = 1e-5 # convergence criterium: integrated diff between local and impurity spectral function
eps_occupancy = 1e-4 # convergence criterium: difference from the target occupancy (occupancy_goal)
To accelerate the convergence to the converged result, one can make use of two different mixing schemes: linear mixing (with parameter \(\alpha\)) or Broyden mixing (with initial Jacobian \(J=-1/\alpha\)). The direct DMFT iteration corresponds to linear mixing with \(\alpha=1\) which actually corresponds to no mixing at all.
mixing_method = "broyden" # "linear" or "broyden" mixing; the quantity being mixed is the hybridisation function
occup_method = "adjust" # "adjust" or None; adjust=shift mu so that the current GF meets the occupancy goal
alpha = 0.5 # mixing parameter from both linear and Broyden mixing
For adjusting the chemical potential, it is useful to go through several DMFT self-consistency cycles (using the same self-energy obtained in the last NRG run) to obtain an improved estimate. In practice order 10 cycles are quite sufficient in all cases.
max_mu_adjust = 10 # number of cycles for adjusting the value of the chemical potential
The logarithmic discretization mesh should span the full frequency range where the spectral function has non-negligible value. The common ratio should be sufficiently small to produce a fine discretization mesh.
mesh_max = 10.0 # logarithmic mesh: maximum frequency (should be large enough to contain the full spectrum)
mesh_min = 1e-5 # logarithmic mesh: minimum frequency (should be smaller than the temperature T)
mesh_ratio = 1.01 # logarithmic mesh: common ratio of the mesh (1.01 is usually a good choice)
The NRG has issues with regions of very low density of states in the band. To control those we define a minimal value of the hybridisation strenth \(\Gamma=-\mathrm{Im}(\Delta)\):
Delta_min = 1e-5 # minimal value for the hybridisation function; if too large, it produces incorrect spectral distribution,
# if too small, it leads to discretization problems (1e-5 is usually a good and rather safe choice)
The Python dictionary observables contains the list of all operators whose expectation values we want to evaluate.
Some standard choices are the impurity occupancy and impurity occupancy squared, in the presence of magnetic field
also the impurity magnetization, and in the presence of phonon modes we might be interested in the phonon number
and oscillator displacement operators.
observables = ["n_d", "n_d^2"]
if B is not None:
observables.extend(["SZd"])
if "Holstein" in model:
observables.extend(["nph", "displ", "displ^2"])
We then set some further global variables that are used in the DMFT iteration:
itern = 0 # iteration counter
verbose = verbose and mpi.is_master_node() # output is only produced by the master node
store_steps = store_steps and mpi.is_master_node() # output is only produced by the master node
symtype = ("QS" if B is None else "QSZ")
One may notice that in the presence of the magnetic field, the code automatically switches to the appropriate
symmetry type (QSZ, meaning conserved total charge and the z-component of total spin).
We are now ready to construct the impurity solver object:
S = Solver(model=model, symtype=symtype, mesh_max=mesh_max, mesh_min=mesh_min, mesh_ratio=mesh_ratio)
S.set_verbosity(verbose_nrg)
We usually want to suppress the detailed output from the NRG by setting verbose_nrg to False. (In fact,
what this does is to redirect the output to a file named log in the temporary directory, so that
the calculation can still be monitored.)
We define some useful functions for working with BlockGf objects:
newG = lambda : S.G_w.copy() # Creates a BlockGf object of appropriate structure
nr_blocks = lambda bgf : len([bl for bl in bgf.indices]) # Returns the number of blocks in a BlockGf object
block_size = lambda bl : len(S.G_w[bl].indices[0]) # Matrix size of Green's functions in block 'bl'
identity = lambda bl : np.identity(block_size(bl)) # Returns the identity matrix in block 'bl'
The key NRG solver parameters that control the discretization and the truncation are set here:
# Solve Parameters
sp = { "T": T, "Lambda": 2.0, "Nz": 4, "Tmin": 1e-5, "keep": 10000, "keepenergy": 10.0 }
The default value of the discretization parameter \(\Lambda=2\) is an excellent choice for most cases, as long as the calculation remains
feasible in reasonable time with the given computing resources; if not, one should increase \(\Lambda\).
We average over Nz=4 interleaves discretization grids; this is a good number for \(\Lambda=2\), but should be
increased for coarses discretization to reduce the artifacts. Furthermore, it needs to be increased for high-resolution calculations
where the broadening parameters are lower.
The parameter Tmin controls the length of the Wilson chain and it has to be smaller than the physical temperature
\(T\). The truncation at rescaled energy 10 is sufficient to achieve good convergence of the spectral functions,
while beyond 12 there are hardly any improvements. keep (the maximum number of multiplets kept) should
be some large number and it is basically constrained by the available memory.
Model parameters are contained as a separate Python dictionary. We set it here:
mp = { "U1": U, "B1": B, "omega": omega, "g1": g1, "n1": n1 }
mp = { k: v for k,v in mp.items() if v is not None }
sp["model_parameters"] = mp
We remove parameters that are not relevant (those set to None). Furthermore, the interface checks if
all required parameters for a given template are in fact defined, to prevent user errors that are
difficult to detect (e.g. typos).
We provide a function for updating the chemical potential. It sets the global variable mu, as well
as updates the impurity level to \(\epsilon_d=-\mu\).
global mu
mu = new_mu
sp["model_parameters"]["eps1"] = -mu
If some details of the NRG calculation need to be tweaked, this can be achieve by setting the appropriate
parameters in nrg_params.
#nrgp["bandrescale"] = 5.0
S.set_nrg_params(**nrgp)
Now we define the function ht0 which evaluates the Hilbert transform of the density of states for
a given complex value \(z\). The analytical expression for the Bethe lattice is hard-coded, while the
generic case is handled by reading the DOS from a file and doing the Hilbert transform numerically.
ht1 = lambda z: 2*(z-1j*np.sign(z.imag)*np.sqrt(1-z**2)) # Analytical expression for Hilbert transform of Bethe lattice DOS
global ht0
ht0 = lambda z: ht1(z/Bethe_unit)
else:
table = np.loadtxt(dos)
global dosA
dosA = Gf(mesh=MeshReFreqPts(table[:,0]), target_shape=[])
for i, w in enumerate(dosA.mesh):
dosA[w] = np.array([[ table[i,1] ]])
ht0 = lambda z: hilbert_transform_refreq(dosA, z)
We wrap this in another function ht which ensures that in cases of causality violation (imaginary
part of the self-energy being positive) things do no go terribly wrong. Tiny causality violations in NRG
commonly occur at very low temperatures close to the Fermi level, where \(\mathrm{Im}\Sigma\)
overshoots zero. This problem is well understood and, while annoying, it is usually not an issue, unless
one is explicitly interested in the low-temperature asymptotics.
ht = lambda z: ht0(z.real+1j*(z.imag if z.imag>0.0 else EPS)) # Fix problems due to causality violation
In addition, another potential problem in NRG calculation is the artifacts arising from regions of very low DOS.
We alleviate this problem by providing a minimal cutoff value Delta_min for the hybridisation function:
Delta_fixed = Delta.copy()
for bl in Delta.indices:
for w in Delta.mesh:
for n in range(block_size(bl)): # only diagonal parts
r = Delta[bl][w][n,n].real
i = Delta[bl][w][n,n].imag
Delta_fixed[bl][w][n,n] = r + 1j*(i if i<-Delta_min else -Delta_min)
# Possible improvement: re-adjust the real part so that the Kramers-Kronig relation is maintained
return Delta_fixed
In the following we will make use of a function which calculates a Green’s function given hybridisation, self-energy and chemical potential:
G = newG()
for bl in G.indices:
for w in G.mesh:
G[bl][w] = np.linalg.inv( (w + mu)*identity(bl) - Delta[bl][w] - Sigma[bl][w] ) # !!!
return G
For accurate integration over spectral functions, we make use of cubic interpolation. Note that for frequencies outside the mesh range, the spectral function is assummed to be strictly zero.
lx = np.array(list(G.mesh.values()))
ly = sum( sum( -(1.0/math.pi)*np.array(G[bl].data[:,i,i].imag) for i in index_range(G[bl]) ) # matrix trace
for bl in G.indices ) # sum over blocks
if normalize_to_one:
nr = sum( sum( 1 for i in index_range(G[bl]) ) for bl in G.indices ) # number of contributions
ly = ly/nr # rescale
return interpolate.interp1d(lx, ly, kind='cubic', bounds_error=False, fill_value=0)
This interpolation object is then used e.g. to compute the occupancy given hybridisation, self-energy and chemical potential:
Gtrial = calc_G(Delta, Sigma, mu)
f = interp_A(Gtrial)
with warnings.catch_warnings():
warnings.simplefilter("ignore")
n = integrate.quad(lambda x : 2*f(x)*special.expit(-x/T), -mesh_max, mesh_max)
return n[0]
And this, in turn, is used for occupancy control. The following function determines the optimal value of the
chemical potential for given hybridisation and self-energy. It makes use of root_scalar, which implements
the secant method with initial values shifted by 0.1. This seems to be robust and works quite well in practice.
sol = optimize.root_scalar(lambda x : calc_occupancy(Delta, Sigma, x)-occupancy_goal, x0=mu, x1=mu-0.1)
set_mu(sol.root)
Now comes the key function which calculates the lattice Green’s function, and from this a new
hybridisation function for the effective impurity model. The input is the self-energy (and implicitly the
chemical potential mu):
Gloc = newG()
for bl in Gloc.indices:
for w in Gloc.mesh:
for i in range(block_size(bl)):
for j in range(block_size(bl)): # assuming square matrix
if i == j:
Gloc[bl][w][i,i] = ht(w + mu - Sigma[bl][w][i,i]) # Hilbert-transform
else:
assert abs(Sigma[bl][w][i,j])<1e-10, "This implementation only supports diagonal self-energy"
Gloc[bl][w][i,j] = 0.0
Glocinv = Gloc.inverse()
Delta = newG()
for bl in Delta.indices:
for w in Delta.mesh:
Delta[bl][w] = (w+mu)*identity(bl) - Sigma[bl][w] - Glocinv[bl][w] # !!!
return Gloc, Delta
Note that this implementation is limited to problems with diagonal self-energy (and we test for this explicitly, as a precaution). This means that this code is only suitable for channel/orbital conserving models.
Now we are in position to provide an improved method for adjusting the chemical potential, which
repeatedly recalculates the hybridisation function and readjust the chemical potential. We do this
for max_mu_adjust cycles, with 10 being a good choice.
old_mu = mu
Delta = Delta_in.copy()
for _ in range(max_mu_adjust):
update_mu(Delta, Sigma)
Gloc, Delta = self_consistency(Sigma)
new_mu = mu
if verbose: print("Adjusted mu from %f to %f\n" % (old_mu,new_mu))
return Gloc, Delta
We now define some miscelaneous functions. First we need a way of estimating the difference between two Green’s functions. This will be used in the following for convergence testing. Here we again make use of interpolation, and we define the difference as the integral over all frequencies (within the discretization mesh window) of the difference squared. We ignore possible warnings about loss of accuracy in integration.
f_a = interp_A(a)
f_b = interp_A(b)
with warnings.catch_warnings():
warnings.simplefilter("ignore")
diff = integrate.quad(lambda x : (f_a(x)-f_b(x))**2, -mesh_max, mesh_max)
return diff[0]
We now define some functions for saving Green’s function object to tabulated ASCII files. The filename is just
the prefix. To this we append the block name and the suffix dat.
def save_Gf(fn, gf):
f = open(fn, "w")
for w in gf.mesh:
z = gf[w]
f.write("%f %f %f\n" % (w, z.real, z.imag))
# Save all blocks for a block GF to tabulated ASCII files
def save_BlockGf(fn, bgf):
for bl in bgf.indices:
save_Gf(fn + "_" + bl + ".dat", bgf[bl])
# Save a spectral function (-1/Pi Im GF) to a tabulated ASCII file
def save_A(fn, gf):
f = open(fn, "w")
for w in gf.mesh:
z = gf[w]
f.write("%f %f\n" % (w, -1.0/math.pi*z.imag))
# Save spectral functions for all blocks of the block GF
def save_BlockA(fn, bgf):
for bl in bgf.indices:
save_A(fn + "_" + bl + ".dat", bgf[bl])
In addition, the full set of results is stores as an HDF5 file. We store the relevant items from the solver object, as well as some additional quantities (self-energy-improved impurity Green’s function, local lattice Green’s function, chemical potential).
with HDFArchive(fn, 'w') as arch:
arch["S"] = S
# Global variables
arch["Gself"] = Gself
arch["Gloc"] = Gloc
arch["mu"] = mu
When restarting from stored results, we only need to read in the self-energy and the chemical potential, since this is sufficient to initialize the DMFT calculation.
with HDFArchive(fn, 'r') as arch:
Sigma = arch["S"].Sigma_w
mu = arch["mu"]
return Sigma, mu
if verbose: print("Starting from stored results in file %s\n" % fn)
return load_Sigma_mu(fn) # Load data from an HDF5 archive
If not starting from a file, the initial self-energy includes the Hartree term, and the chemical potential is set to \(U/2\):
if verbose: print("Starting from scratch\n")
Sigma = newG()
for bl in Sigma.indices:
for w in Sigma.mesh:
Sigma[bl][w] = U*occupancy_goal/2.0 # Initialize self-energy with the Hartree shift
mu = U/2.0 # initial approximaiton for the chemical potential
return Sigma, mu
The DMFT loop is terminated by raising exceptions. We define them here:
class Converged(Exception):
def __init__(self, message):
self.message = message
# Exception to raise when convergence is not reached (e.g. maximum nr of iterations exceeded)
class FailedToConverge(Exception):
def __init__(self, message):
self.message = message
Now comes the main part, the function which performs a DMFT step. Let us disect it. First we update the iteration counter and initialize the hybridisation function in the impurity solver object after fixing it to a minimal value.
global itern
itern += 1
if verbose:
print("Iteration %i min_iter=%i max_iter=%i\n" % (itern, min_iter, max_iter))
Delta_in_fixed = fix_hyb_function(Delta_in, Delta_min)
S.Delta_w << Delta_in_fixed
We now solve the impurity model, calculate an improved estimate of the impurity GF using the “self-energy trick”, and apply the self-consistency equation:
global Gself, Gloc, Gloc_prev
Gself = calc_G(Delta_in_fixed, S.Sigma_w, mu) # impurity GF ("self-energy-trick" improved)
Gloc, Delta = self_consistency(S.Sigma_w) # apply the DMFT self-consistency equation
We calculate some quantities for estimating the degree of convergence:
diff_prev = gf_diff(Gloc, Gloc_prev) # difference between two consecutively computed local latice GFs
Gloc_prev = Gloc.copy()
occupancy = calc_occupancy(Delta, S.Sigma_w, mu)
diff_occupancy = abs(occupancy-occupancy_goal) # this difference is used as the measure of deviation
("diff_occupancy", diff_occupancy), ("occupancy", occupancy)])
for i in observables:
stats[i] = S.expv[i]
header_string = fmt_str_header(len(stats)).format(*[i for i in stats.keys()])
stats_string = fmt_str(len(stats)).format(*[i for i in stats.values()])
if mpi.is_master_node():
if itern == 1: stats_file.write(header_string)
stats_file.write(stats_string)
if verbose: print("stats: %sstats: %s" % (header_string, stats_string))
Furthermore, if store_steps=True, we save the information about each DMFT step in a separate subdirectory:
os.mkdir(str(itern)) # one subdirectory per iteration
save_BlockGf(str(itern)+"/Delta", Delta_in_fixed)
save_BlockGf(str(itern)+"/Sigma", S.Sigma_w) # self-energy
save_BlockGf(str(itern)+"/G", Gloc) # local lattice Green's function
save_BlockA(str(itern)+"/A", Gloc) # spectral function of local lattice GF
store_result(str(itern)+"/"+solution_filename, S)
At each step, we also store the HDF5 to permit the checkpoint/restart functionality:
store_result(checkpoint_filename, S) # for checkpoint/restart functionality
We now test for convergence (and for failure to reach the convergence within the maximal number of iterations):
diff_prev < eps_prev and
diff_occupancy < eps_occupancy and
itern >= min_iter):
raise Converged(stats_string)
if (itern == max_iter):
raise FailedToConverge(stats_string)
We are nearly done now. We adjust the chemical potential and return the optimized hybridisation function:
Gloc, Delta = adjust_mu(Delta, S.Sigma_w) # here we update mu to get closer to target occupancy
return Delta
There are two types of driving routines for the full DMFT iteration. The first is based on linear mixing:
Delta_in = Delta.copy()
while True:
Delta_out = dmft_step(Delta_in)
newDelta = alpha*Delta_out + (1-alpha)*Delta_in
Delta_in << newDelta
The second one is based on Broyden mixing. The idea is that the DMFT iteration can be thought of as a
procedure for finding the fixed point of the mutli-dimensional function \(F(\Delta_\mathrm{in}) = \Delta_\mathrm{out}[\Delta_\mathrm{in}] - \Delta_\mathrm{in}\),
where \(\Delta_\mathrm{in}\) is the input while \(\Delta_\mathrm{out}\) the output of one DMFT step.
For ths purpose we make use of optimize.broyden1 from scipy, but in order to do so, we
need to wrap/unwrap the Green’s function objects to/from numpy arrays, which adds a degree of complexity:
def gf_to_nparray(gf):
return gf.data.flatten()
# Stack elements from all blocks in a block BF
def bgf_to_nparray(bgf):
return np.hstack((bgf[bl].data.flatten() for bl in bgf.indices))
# Convert a linear numpy array to a single GF object
def nparray_to_gf(a, gf):
b = a.reshape(gf.data.shape)
gf.data[:,:,:] = b[:,:,:]
# Extract blocks from linear numpy array to a block GF
def nparray_to_bgf(a):
G = newG()
split = np.split(a, nr_blocks(G)) # Here we assume all blocks are equally large
for i, bl in enumerate(G.indices):
nparray_to_gf(split[i], G[bl])
return G
# DMFT driver routine with Broyden mixing: the goal is to find a root of the function F(Delta)=dmft_step(Delta)-Delta.
F = lambda Delta : dmft_step(Delta)-Delta
npF = lambda x : bgf_to_nparray(F(nparray_to_bgf(x)))
def solve_with_broyden_mixing(Delta, alpha):
xin = bgf_to_nparray(Delta)
optimize.broyden1(npF, xin, alpha=alpha, reduction_method="svd", max_rank=10, verbose=verbose, f_tol=1e-99) # Loop forever (small f_tol!)
Here we pick the method:
if (mixing_method=="linear"):
solve_with_linear_mixing(Delta, alpha)
if (mixing_method=="broyden"):
solve_with_broyden_mixing(Delta, alpha)
The problem is initialized here:
set_mu(mu0) # Set global variable mu and update the impurity level (epsilon_d = -mu)
Gloc, Delta = self_consistency(Sigma0) # Initialize local GF Gloc and hybridization function Delta
Gloc, Delta = adjust_mu(Delta, Sigma0) # Adjust mu (and return updated Glocal and Delta as well)
Gloc_prev = Gloc.copy() # Store copy for checking convergence
And the solver is started here:
solve(Delta)
except Converged as c:
if verbose: print("Converged: %s" % c.message)
if mpi.is_master_node():
store_result(solution_filename, S) # full converged results as an HDF5 file
os.remove(checkpoint_filename) # checkpoint file is no longer needed
save_BlockA("A", Gloc) # converged spectral function for quick plotting
except FailedToConverge as c:
if verbose: print("Failed to converge: %s" % c.message) # ... but restart is possible from the checkpoint file
This completes the exposition. Have fun running DMFT(NRG) calculations!
Note: this script has not yet been tested for the case of matrix-valued Green’s functions.