TPSC approximation#
Author: Niklas Witt
Theory of TPSC#
The Two-Particle Self-Consistent (TPSC) approximation is a non-perturbative semi-analytical method that was first introduced by Vilk and Tremblay [Vilk and Tremblay, 1997]. TPSC can be used to study magnetic fluctuations, while it also obeys the Mermin-Wagner theorem, i.e., a phase transtition at finite temperatures is prohibited in one or two dimensions. In addition, the TPSC method satisfies several conservation laws, sum rules and the Pauli principle (actually, it is constructed in a way to fulfill these, since they are used to determine model parameters self-consistently). TPSC is applicable in the weak to intermediate coupling regime, but it breaks down in the strong coupling regime and it cannot describe the Mott transition unlike other non-perturbative methods like Dynamical Mean-Field Theory (DMFT).
For a (pedagogical) review, please have a look at [Allen et al., 2004, Tremblay, 2012] for the single-orbital case implemented here and [Zantout et al., 2021] for the more complex multi-orbital theory.
Set of TPSC equations#
We review the set of equations that need to be solved in the TPSC approximation assuming a one-band Hubbard model with interaction \(U\) (it is not so easy to extend TPSC to models with more parameters, since sum rules to determine the additional parameters self-consistently need to be found) in the paramagnetic phase (SU(2) symmetric), i.e., \(n = \langle n\rangle = 2n_{\sigma}\) for the electron filling \(n\). TPSC is constructed in a way to fulfill certain sum rules and the Pauli principle in the form \(\langle n^2\rangle = \langle n\rangle\). The control quantitites are spin and charge correlation function (susceptibilities) which are evaluated in a Random-Phase-Approximation (RPA) like fashion
with the irreducible susceptibility (“bubble diagram”)
where \(\nu^\mathrm{B} = n\pi T\) with even \(n\) [\(\nu^\mathrm{F} = n\pi T\) with odd \(n\)] and \(\boldsymbol{q}\) [\(\boldsymbol{k}\)] are bosonic [fermionic] Matsubara frequencies and momentum at temperature \(T\), \(N_{\boldsymbol{k}}\) denotes the number of \(\boldsymbol{k}\)-points, and \(G_0(\mathrm{i}\nu^\mathrm{F},\boldsymbol{k}) = [\mathrm{i}\nu^\mathrm{F} - (\varepsilon_{\boldsymbol{k}}-\mu)]^{-1}\) is the bare (non-interacting) Green function with with single-particle dispersion \(\varepsilon_{\boldsymbol{k}}\) and chemical potential \(\mu\). Please note that sometimes a factor of 2 is included in the definition of \(\chi_0\) leading to slightly different factors in all equations given here. The convolution sum to calculate \(\chi_0\) can be easily evaluated by Fourier transforming to imaginary-time and real space, resulting in a simple multiplication
In our practical implementation, we will perform this step using the sparse-ir package. A similar calculation is necessary to set the chemical potential \(\mu\) for fixed electron density \(n\), as the equation
(factor 2 from spin degeneracy and \(0^+ = \lim_{\eta\to 0+} \eta\)) needs to be solved by using some root finding algorithm like bisection method or Brent’s method. The Fourier transformation to \(\tau=0^+\) can be easily performed with the sparse-ir package.
The above RPA definition of spin and charge susceptibility violate the Pauli principle. In TPSC, we overcome this problem by introducing two effective, renormalized interactions (“irreducible vertices”) \(U_{\mathrm{sp}}\) and \(U_{\mathrm{ch}}\) that enter spin and charge correlation functions as
These two effetive interactions are determined by the two local sum rules
Both sum rules can be exactly derived from the Pauli principle (\(\langle n^2\rangle = \langle n\rangle\)). In principle, we can now determine \(U_{\mathrm{sp}}\) and \(U_{\mathrm{ch}}\) from local-spin and local-charge sum rule if we knew the double occupancy \(\langle n_{\uparrow}n_{\downarrow}\rangle\). TPSC makes the ansatz
which reproduces Kanamori-Brueckner type screening. The four equations above form a set of self-consistent equations for either \(U_{\mathrm{sp}}\) or equivalently \(\langle n_{\uparrow}n_{\downarrow}\rangle\). In practice, we treat \(U_{\mathrm{sp}}\) as the parameter to be determined self-consistently by inserting the ansatz in the local-spin sum rule. Effectively, we then need to find the root of the function
Afterwards we can calculate the double occupancy \(\langle n_{\uparrow}n_{\downarrow}\rangle = \frac{U_{\mathrm{sp}}}{4U} n^2\) and then perform a similar root finding for \(U_{\mathrm{ch}}\) from the function
In TPSC, a self-energy \(\Sigma\) can be derived that is calculated from the interaction [Moukouri et al., 2000]
The self-energy itself is given by a convolution in \((\mathrm{i}\nu^\mathrm{F}, \boldsymbol{k})\) space
which Fourier transformed to \((\tau,\boldsymbol{r})\) space takes the form
The interacting Green function is determined by the Dyson equation
Notes on practical implementation#
When implementing the TPSC, a few points need to be treated carefully which we adress in the following:
The constant Hartree term \(V_{\mathrm{H}} = U\) in the interaction \(V\) and respective self-energy term \(\Sigma_H = U n_{\sigma} = U\frac{n}{2}\) can be absorbed into the definition of the chemical potential \(\mu\). Otherwise we would have to treat this term separately, since the IR basis cannot model the \(\delta\)-peak in frequency space well/compactly. In this case, evaluate the term analytically and add it after the Fourier transformation step.
An upper bound for the renormalized spin vertex \(U_{\mathrm{sp}}\) exists. Since the denominator spin susceptibility \(\chi_{\mathrm{sp}}\) should not diverge, the upper bound is given by the RPA critical interaction value \(U_{\mathrm{crit}} = 1/\mathrm{max}\{\chi_0\}\). Mathematically, the function \(f(U_{\mathrm{sp}}) = 2\sum \chi_{\mathrm{sp}}(U_{\mathrm{sp}}) - n + \frac{U_{\mathrm{sp}}}{2U}n^2\), from which \(U_{\mathrm{sp}}\) is determined, turns unstable for \(U_{\mathrm{sp}} \geq U_{\mathrm{crit}}\) (try plotting \(f(U_{\mathrm{sp}})\)!). At this point, TPSC is not applicable and, e.g., the temperature \(T\) is too low or the (unrenormalized) interaction \(U\) too large.
An internal accuracy check \(\frac{1}{2}\mathrm{Tr}(\Sigma G) = U \langle n_{\uparrow} n_{\downarrow}\rangle\) can be employed to test the validity of TPSC (not done here).
Code implementation#
We are implementing TPSC for the simple case of a square lattice model with dispersion \(\varepsilon_{\boldsymbol{k}} = -2t\,[\cos(k_x) + \cos(k_y)]\) with nearest-neighbor hopping \(t\) which sets the energy scale of our system (bandwidth \(W = 8t\)). First, we load all necessary basic modules that we are going to need in implementing TPSC and visualizing results:
import numpy as np
import scipy as sc
import scipy.optimize
from warnings import warn
import sparse_ir
%matplotlib inline
import matplotlib.pyplot as plt
Parameter setting#
### System parameters
t = 1 # hopping amplitude
W = 8*t # bandwidth
wmax = 10 # set wmax >= W
T = 0.1 # temperature
beta = 1/T # inverse temperature
n = 0.85 # electron filling per lattice site, summed over both spins (n=1: half filling)
U = 4 # Hubbard interaction
### Numerical parameters
nk1, nk2 = 24, 24 # number of k_points along one repiprocal crystal lattice direction k1 = kx, k2 = ky
nk = nk1*nk2
IR_tol = 1e-10 # desired accuary for l-cutoff of IR basis functions
Generating meshes#
We need to generate a \(\boldsymbol{k}\)-mesh as well as set up the IR basis functions on a sparse \(\tau\) and \(\mathrm{i}\nu\) grid. Then we can calculate the dispersion on this mesh. In addition, we set calculation routines to Fourier transform \(k\leftrightarrow r\) and \(\tau\leftrightarrow \mathrm{i}\nu\) (via IR basis).
#### Initiate fermionic and bosonic IR basis objects
IR_basis_set = sparse_ir.FiniteTempBasisSet(beta, wmax, eps=IR_tol)
class Mesh:
"""
Holding class for k-mesh and sparsely sampled imaginary time 'tau' / Matsubara frequency 'iwn' grids.
Additionally it defines the Fourier transform routines 'r <-> k' and 'tau <-> l <-> wn'.
"""
def __init__(self,IR_basis_set,nk1,nk2):
self.IR_basis_set = IR_basis_set
# generate k-mesh and dispersion
self.nk1, self.nk2, self.nk = nk1, nk2, nk1*nk2
self.k1, self.k2 = np.meshgrid(np.arange(self.nk1)/self.nk1, np.arange(self.nk2)/self.nk2)
self.ek = -2*t*( np.cos(2*np.pi*self.k1) + np.cos(2*np.pi*self.k2) ).reshape(nk)
# lowest Matsubara frequency index
self.iw0_f = np.where(self.IR_basis_set.wn_f == 1)[0][0]
self.iw0_b = np.where(self.IR_basis_set.wn_b == 0)[0][0]
### Generate a frequency-momentum grid for iw_n and ek (in preparation for calculating the Green function)
# frequency mesh (for Green function)
self.iwn_f = 1j * self.IR_basis_set.wn_f * np.pi * T
self.iwn_f_ = np.tensordot(self.iwn_f, np.ones(nk), axes=0)
# ek mesh
self.ek_ = np.tensordot(np.ones(len(self.iwn_f)), self.ek, axes=0)
def smpl_obj(self, statistics):
""" Return sampling object for given statistic """
smpl_tau = {'F': self.IR_basis_set.smpl_tau_f, 'B': self.IR_basis_set.smpl_tau_b}[statistics]
smpl_wn = {'F': self.IR_basis_set.smpl_wn_f, 'B': self.IR_basis_set.smpl_wn_b }[statistics]
return smpl_tau, smpl_wn
def tau_to_wn(self, statistics, obj_tau):
""" Fourier transform from tau to iwn via IR basis """
smpl_tau, smpl_wn = self.smpl_obj(statistics)
obj_tau = obj_tau.reshape((smpl_tau.tau.size, self.nk1, self.nk2))
obj_l = smpl_tau.fit(obj_tau, axis=0)
obj_wn = smpl_wn.evaluate(obj_l, axis=0).reshape((smpl_wn.wn.size, self.nk))
return obj_wn
def wn_to_tau(self, statistics, obj_wn):
""" Fourier transform from iwn to tau via IR basis """
smpl_tau, smpl_wn = self.smpl_obj(statistics)
obj_wn = obj_wn.reshape((smpl_wn.wn.size, self.nk1, self.nk2))
obj_l = smpl_wn.fit(obj_wn, axis=0)
obj_tau = smpl_tau.evaluate(obj_l, axis=0).reshape((smpl_tau.tau.size, self.nk))
return obj_tau
def k_to_r(self,obj_k):
""" Fourier transform from k-space to real space, A(r) = sum_k exp(-ik.r) A(k), without a factor 1/nk """
obj_k = obj_k.reshape(-1, self.nk1, self.nk2)
obj_r = np.fft.fftn(obj_k,axes=(1,2))
obj_r = obj_r.reshape(-1, self.nk)
return obj_r
def r_to_k(self,obj_r):
""" Fourier transform from real space to k-space, A(k) = (1/nk^2) sum_r exp(ik.r) A(r) (1/nk from ifftn and 1/nk here),
so that the product of two k_to_r arrays A(r)B(r) gives the convolution (1/nk) sum_q A(q)B(k-q) """
obj_r = obj_r.reshape(-1, self.nk1, self.nk2)
obj_k = np.fft.ifftn(obj_r,axes=(1,2))/self.nk
obj_k = obj_k.reshape(-1, self.nk)
return obj_k
TPSC solver#
We wrap the calculation steps of TPSC (i.e. determining \(U_{\mathrm{sp}},U_{\mathrm{ch}}\)) in the TPSCSolver class. We use the Mesh class defined above to perform calculation steps.
class TPSCSolver:
def __init__(self, mesh, U, n, U_sfc_tol=1e-12, verbose=True):
"""
Solver class to calculate the TPSC method.
After initializing the Solver by `solver = TPSCSolver(mesh, U, n, **kwargs)` it
can be run by `solver.solve()`.
"""
## set internal parameters for the solver
self.U = U
self.n = n
self.mesh = mesh
self.U_sfc_tol = U_sfc_tol
self.verbose = verbose
## set initial Green function and irreducible susceptibility
self.sigma = 0
self.mu = 0
self.mu_calc()
self.gkio_calc(self.mu)
self.grit_calc()
self.ckio_calc()
# determine critical U_crit = 1/max(chi0) as an upper bound to U_sp
self.U_crit = 1/np.amax(self.ckio.real)
#%%%%%%%%%%% Solving instance
def solve(self):
"""
Determine spin and charge vertex self-consistently from sum rules and calculate self-energy.
"""
# determine spin vertex U_sp
self.spin_vertex_calc()
# set double occupancy from Kanamori-Bruckner screening
self.docc_calc()
# determine charge vertex U_ch
self.charge_vertex_calc()
# set spin and charge susceptibility
self.chi_spin = self.RPA_term_calc( self.U_sp)
self.chi_charge = self.RPA_term_calc(-self.U_ch)
# calculate interaction, self-energy and interacting Green function
self.V_calc()
self.sigma_calc()
self.mu_calc()
self.gkio_calc(self.mu)
#%%%%%%%%%%% Calculation steps for self.energy
def gkio_calc(self, mu):
""" Calculate Green function G(iw,k) """
self.gkio = (self.mesh.iwn_f_ - (self.mesh.ek_ - mu) - self.sigma)**(-1)
def grit_calc(self):
""" Calculate real space Green function G(tau,r) [for calculating chi0 and sigma] """
# Fourier transform
grit = self.mesh.k_to_r(self.gkio)
self.grit = self.mesh.wn_to_tau('F', grit)
def ckio_calc(self):
""" Calculate irreducible susciptibility chi0(iv,q) """
ckio = self.grit * self.grit[::-1, :]
# Fourier transform
ckio = self.mesh.r_to_k(ckio)
self.ckio = self.mesh.tau_to_wn('B', ckio)
def V_calc(self):
""" Calculate interaction V(tau,r) from RPA-like spin and charge susceptibility """
V = self.U/4 * (3*self.U_sp*self.chi_spin + self.U_ch*self.chi_charge)
# Constant Hartree Term V ~ U needs to be treated extra, since it cannot be modeled by the IR basis compactly.
# In the single-band case, the Hartree term can be absorbed into the chemical potential.
# Fourier transform
V = self.mesh.k_to_r(V)
self.V = self.mesh.wn_to_tau('B', V)
def sigma_calc(self):
""" Calculate self-energy Sigma(iwn,k) """
sigma = self.V * self.grit
# Fourier transform
sigma = self.mesh.r_to_k(sigma)
self.sigma = self.mesh.tau_to_wn('F', sigma)
#%%%%%%%%%%% Determining spin and charge vertex
def RPA_term_calc(self, U):
""" Set RPA-like susceptibility """
chi_RPA = self.ckio / (1 - U*self.ckio)
return chi_RPA
def chi_qtrace_calc(self, U):
""" Calculate (iv_m, q) trace of chi_RPA term """
# chi_qtrace = sum_(m,q) chi(iv_m,q)
chi_RPA = self.RPA_term_calc(U)
chi_trace = np.sum(chi_RPA, axis=1)/self.mesh.nk
chi_trace_l = self.mesh.IR_basis_set.smpl_wn_b.fit(chi_trace)
chi_trace = self.mesh.IR_basis_set.basis_b.u(0)@chi_trace_l
return chi_trace.real
def docc_calc(self):
""" Calculate double occupancy from Kanamori-Bruckner type screening """
self.docc = 0.25 * self.U_sp/self.U * self.n**2
def spin_vertex_calc(self):
""" Determine U_sp self-consistently from local sum rule """
# interval [U_a, U_b] for root finding
U_a = 0
U_b = np.floor(self.U_crit*100)/100
chi_trace = self.chi_qtrace_calc
sfc_eq = lambda U_sp : 2*chi_trace(U_sp) - self.n + 0.5*(U_sp/self.U)*self.n**2
if sfc_eq(U_b) > 0:
self.U_sp = sc.optimize.brentq(sfc_eq, U_a, U_b, rtol = self.U_sfc_tol)
print("Setting U_sp = {}".format(self.U_sp))
else:
print("Error")
warn("System underwent phase transition, U^sp > U_crit = {}! U is too large or T too low for given doping.".format(self.U_crit))
def charge_vertex_calc(self):
""" Determine U_ch self-consistently from local sum rule """
# interval [U_a, U_b] for root finding
U_a = 0
U_b = 100
chi_trace = self.chi_qtrace_calc
sfc_eq = lambda U_ch : 2*chi_trace(-U_ch) - self.n - 2*self.docc + self.n**2
self.U_ch = sc.optimize.brentq(sfc_eq, U_a, U_b, rtol = self.U_sfc_tol)
#%%%%%%%%%%% Setting chemical potential mu
def calc_electron_density(self, mu):
""" Calculate chemical potential mu from Green function """
self.gkio_calc(mu)
gio = np.sum(self.gkio,axis=1)/self.mesh.nk
g_l = self.mesh.IR_basis_set.smpl_wn_f.fit(gio)
g_tau0 = self.mesh.IR_basis_set.basis_f.u(0)@g_l
n = 1 + np.real(g_tau0)
n = 2*n #for spin
return n
def mu_calc(self):
""" Find chemical potential for a given filling n0 via brentq root finding algorithm """
n_calc = self.calc_electron_density
n0 = self.n
f = lambda mu : n_calc(mu) - n0
self.mu = sc.optimize.brentq(f, np.amax(self.mesh.ek)*3, np.amin(self.mesh.ek)*3)
Execute TPSC solver#
# initialize calculation
IR_basis_set = sparse_ir.FiniteTempBasisSet(beta, wmax, eps=IR_tol)
mesh = Mesh(IR_basis_set, nk1, nk2)
solver = TPSCSolver(mesh, U, n)
# perform TPSC calculation
solver.solve()
Setting U_sp = 2.1011348343400282
Visualize results#
# plot 2D k-dependence of lowest Matsubara frequency of e.g. green function
plt.pcolormesh(2*mesh.k1.reshape(nk1,nk2), 2*mesh.k2.reshape(nk1,nk2), np.real(solver.gkio[mesh.iw0_f].reshape(mesh.nk1,mesh.nk2)), shading='auto')
ax = plt.gca()
ax.set_xlabel(r'$k_x/\pi$')
ax.set_xlim([0,2])
ax.set_ylabel(r'$k_y/\pi$')
ax.set_ylim([0,2])
ax.set_aspect('equal')
ax.set_title(r'Re $G(k,\mathrm{i}\pi T)$')
plt.colorbar()
plt.show()
# plot 2D k-dependence of lowest Matsubara frequency of e.g. self-energy
plt.pcolormesh(2*mesh.k1.reshape(nk1,nk2), 2*mesh.k2.reshape(nk1,nk2), np.imag(solver.sigma[mesh.iw0_f].reshape(mesh.nk1,mesh.nk2)), shading='auto')
ax = plt.gca()
ax.set_xlabel(r'$k_x/\pi$')
ax.set_xlim([0,2])
ax.set_ylabel(r'$k_y/\pi$')
ax.set_ylim([0,2])
ax.set_aspect('equal')
ax.set_title(r'Im $\Sigma(k,\mathrm{i}\pi T)$')
plt.colorbar()
plt.show()
# plot 2D k-dependence of lowest Matsubara frequency of e.g. chi_spin
plt.pcolormesh(2*mesh.k1.reshape(nk1,nk2), 2*mesh.k2.reshape(nk1,nk2), np.real(solver.chi_spin[mesh.iw0_b].reshape(mesh.nk1,mesh.nk2)), shading='auto')
ax = plt.gca()
ax.set_xlabel(r'$k_x/\pi$')
ax.set_xlim([0,2])
ax.set_ylabel(r'$k_y/\pi$')
ax.set_ylim([0,2])
ax.set_aspect('equal')
ax.set_title(r'$\chi_{\mathrm{sp}}(k,\mathrm{i}\nu^\mathrm{B}=0)$')
plt.colorbar()
plt.show()
Example: Interaction dependent renormalization#
As a simple example demonstration of our sparse-ir TPSC code developed above, we will reproduce Fig. 2 of [Vilk and Tremblay, 1997]. It shows the \(U\) dependence of renormalized/effective spin and charge interactions \(U_{\mathrm{sp}}\) and \(U_{\mathrm{ch}}\) (irreducible vertices) at half filling \(n=1\) and \(T>T_{\mathrm{crit}}\) for all considered \(U\) (i.e. \(U_{\mathrm{sp}}<U_{\mathrm{crit}}\) is ensured).
You can simply execute the following two code blocks which will first perform the calculation and then generate a figure like in the reference above.
#%%%%%%%%%%%%%%% Parameter settings
print('Initialization...')
# system parameters
t = 1 # hopping amplitude
n = 1 # electron filling per lattice site, summed over both spins (n=1: half filling)
T = 0.4 # temperature
beta =1/T
U_array = np.linspace(1e-2,5,51) # Hubbard interaction
W = 8*t # bandwidth
wmax = 10 # set wmax >= W
# numerical parameters
nk1, nk2 = 24, 24 # k-mesh sufficiently dense!
nk = nk1*nk2
IR_tol = 1e-8 # desired accuary for l-cutoff of IR basis functions
# initialize meshes
IR_basis_set = sparse_ir.FiniteTempBasisSet(beta, wmax, eps=IR_tol)
mesh = Mesh(IR_basis_set, nk1, nk2)
# set initial self_energy - will be set to previous calculation step afterwards
sigma_init = 0
# empty arrays for results later
U_sp_array = np.empty((len(U_array)))
U_ch_array = np.empty((len(U_array)))
#%%%%%%%%%%%%%%% Calculations for different U values
print("Start TPSC loop...")
for U_it, U in enumerate(U_array):
#print("Now: U = {:.1f}".format(U))
# TPSC solver
solver = TPSCSolver(mesh, U, n, verbose=False)
solver.solve()
# save data for plotting
U_sp_array[U_it] = solver.U_sp
U_ch_array[U_it] = solver.U_ch
print("Finished. Plotting now.")
#%%%%%%%%%%%%%%%% Plot results
plt.plot(U_array, U_ch_array, '-', label=r'$U_{\mathrm{ch}}$')
plt.plot(U_array, U_sp_array, '-', label=r'$U_{\mathrm{sp}}$')
ax = plt.gca()
ax.set_xlabel(r'$U$', fontsize=12)
ax.set_xlim([0,5])
ax.set_ylabel(r'Interaction', fontsize=12)
ax.set_ylim([0,20])
ax.legend(frameon=False, fontsize=12)
Initialization...
Start TPSC loop...
Setting U_sp = 0.009988853647896
Setting U_sp = 0.10844539767253622
Setting U_sp = 0.20462542385193294
Setting U_sp = 0.2984789555689598
Setting U_sp = 0.38995933428489765
Setting U_sp = 0.4790237958356736
Setting U_sp = 0.565634045351875
Setting U_sp = 0.6497568174617774
Setting U_sp = 0.7313644072750964
Setting U_sp = 0.8104351570125639
Setting U_sp = 0.8869538831613245
Setting U_sp = 0.9609122297728379
Setting U_sp = 1.0323089350516337
Setting U_sp = 1.1011500006196007
Setting U_sp = 1.1674487557881674
Setting U_sp = 1.2312258126058024
Setting U_sp = 1.2925089112219597
Setting U_sp = 1.3513326589721484
Setting U_sp = 1.4077381703028404
Setting U_sp = 1.4617726179865584
Setting U_sp = 1.513488708827769
Setting U_sp = 1.5629440990822872
Setting U_sp = 1.6102007660232784
Setting U_sp = 1.6553243524696442
Setting U_sp = 1.6983835006942025
Setting U_sp = 1.7394491910491903
Setting U_sp = 1.7785940990215285
Setting U_sp = 1.8158919824173199
Setting U_sp = 1.8514171081351636
Setting U_sp = 1.8852437256723953
Setting U_sp = 1.9174455922480051
Setting U_sp = 1.9480955523246621
Setting U_sp = 1.9772651724445822
Setting U_sp = 2.0050244307054497
Setting U_sp = 2.0314414589133416
Setting U_sp = 2.0565823344586733
Setting U_sp = 2.080510918245479
Setting U_sp = 2.103288734561642
Setting U_sp = 2.1249748885017796
Setting U_sp = 2.145626016533246
Setting U_sp = 2.1652962658462207
Setting U_sp = 2.1840372983104337
Setting U_sp = 2.2018983151632447
Setting U_sp = 2.2189260988095914
Setting U_sp = 2.2351650684989766
Setting U_sp = 2.2506573469711557
Setting U_sp = 2.265442835514978
Setting U_sp = 2.279559295215254
Setting U_sp = 2.2930424324769056
Setting U_sp = 2.30592598719869
Setting U_sp = 2.3182418222497003
Finished. Plotting now.
<matplotlib.legend.Legend at 0x7ff960707cb0>
