Semi-infinite chain
This tutorial shows how to setup a minimal impurity model, that of a single fermionic state, and solve it using triqs_soehyb. As the first example we pick the analytically solvable case of a semi-infinte chain. (By having the exact solution we will be able to see how the bold hybridization expansion converges as a function of expansion order.)
In the semi-infinite chain with nearest neighbour hopping \(t = 1\) the site with only one neighbour has a single particle Green’s function with semi-circular density of states (DOS). Thus the chain can be reformulated as an impurity problem of a single site (with zero energy) coupled to a hybridization function with semi-cicular DOS.
Initialization
First we spawn a solver instance and setting up a few general system parameters.
[2]:
from triqs_soehyb import Solver
S = Solver(
beta=5.0, # Inverse temperature
gf_struct=[('0', 1)], # Green's function structure 1st index: name, 2nd index: dimension of subspace
eps=1e-8, # Accuracy of Discrete Lehmann Representation (DLR) used for imaginary time response functions
w_max=2.0, # DLR freqiency cut-off (the spectrum of the model must be in the range [-w_max, +w_max]
)
Warning: could not identify MPI environment!
___ ___ _ ___ _____
/ __| ___| __|__| || \ \ / / _ )
\__ \/ _ \ _|___| __ |\ V /| _ \
|___/\___/___| |_||_| |_| |___/ [github.com/TRIQS/soehyb]
beta = 5.0, lamb = 1.00E+01, eps = 1.00E-08, N_DLR = 13
fundamental_operators = [1*c('0',0)]
H_mat.shape = (2, 2)
H_loc = 0
Starting serial run at: 2025-08-17 22:19:09.421767
Hybridization function
In the solver instance the hybridization function is located at S.Delta_tau and we set it to have the semi-circular shape of the semi-infinite chain.
[3]:
from triqs.gf import make_gf_dlr_imtime, make_gf_dlr_imfreq, SemiCircular
for bidx, delta_tau in S.Delta_tau:
delta_w = make_gf_dlr_imfreq(delta_tau)
delta_w << SemiCircular(2.0)
delta_tau[:] = make_gf_dlr_imtime(delta_w)
Local many-body Hamiltonian
The local Hamiltonian of the single-site is in this case trivial with a single state at zero energy.
[4]:
from triqs.operators import c, c_dag
h_int = 0.0 * c_dag('0',0) * c('0',0)
Run solver
Finally we run the solver for a few different expansion orders.
[5]:
from h5 import HDFArchive
max_order = 4
for order in range(1, max_order+1):
S.solve(h_int=h_int, order=order, maxiter=20)
with HDFArchive(f'data_order_{order}.h5', 'w') as A:
A['S'] = S
AdaPol: Hybridization fit accuracy 3.39E-07, using 5 poles.
(Order, N_Diags) = [(1, 1)]
max_order = 1
iter | conv | Z-1
------+----------+-----------
1 | 3.28E-01 | +4.44E-16
2 | 8.18E-02 | +1.89E-14
3 | 1.96E-02 | +1.91E-14
4 | 4.39E-03 | -1.22E-14
5 | 9.45E-04 | +3.06E-14
6 | 1.98E-04 | -4.44E-15
7 | 4.06E-05 | -1.98E-14
8 | 8.23E-06 | -4.55E-14
9 | 1.65E-06 | -1.28E-14
10 | 3.31E-07 | +1.78E-14
11 | 6.58E-08 | -1.50E-14
12 | 1.31E-08 | -8.88E-16
13 | 2.58E-09 | -1.95E-14
14 | 5.11E-10 | +3.11E-15
Timing: incl. excl.
------------------------------------------------------------
Dyson equation: 0.194 0.194 44.5% |-----------------|
Hybridization compression (AAA): 0.079 0.079 18.0% |------|
Pseudo-particle self-energy: 0.006 0.006 1.3% ||
Single particle Green's function: 0.000 0.000 0.0% |
Other: 0.158 0.158 36.2% |-------------|
------------------------------------------------------------
Total: 0.437 100.0%
AdaPol: Hybridization fit accuracy 3.39E-07, using 5 poles.
(Order, N_Diags) = [(1, 1), (2, 10)]
max_order = 2
iter | conv | Z-1
------+----------+-----------
1 | 3.28E-01 | +4.44E-16
2 | 8.18E-02 | +9.10E-15
3 | 1.96E-02 | +2.22E-14
4 | 4.39E-03 | -6.88E-15
5 | 9.45E-04 | +2.89E-15
6 | 1.98E-04 | +1.93E-14
7 | 4.06E-05 | -9.88E-15
8 | 8.23E-06 | -5.55E-15
9 | 1.65E-06 | +1.11E-15
10 | 3.31E-07 | -1.04E-14
11 | 6.58E-08 | +7.55E-15
12 | 1.31E-08 | +1.62E-14
13 | 2.58E-09 | -4.22E-15
14 | 5.11E-10 | -3.67E-14
Timing: incl. excl.
------------------------------------------------------------
Dyson equation: 0.380 0.380 42.1% |----------------|
Hybridization compression (AAA): 0.093 0.093 10.3% |---|
Pseudo-particle self-energy: 0.170 0.170 18.9% |-------|
Single particle Green's function: 0.001 0.001 0.1% |
Other: 0.258 0.258 28.6% |----------|
------------------------------------------------------------
Total: 0.902 100.0%
AdaPol: Hybridization fit accuracy 3.39E-07, using 5 poles.
(Order, N_Diags) = [(1, 1), (2, 10), (3, 100)]
max_order = 3
iter | conv | Z-1
------+----------+-----------
1 | 2.95E-01 | +0.00E+00
2 | 5.22E-02 | -1.52E-14
3 | 6.72E-03 | +2.60E-14
4 | 1.07E-03 | +2.62E-14
5 | 2.58E-04 | -1.31E-14
6 | 6.04E-05 | +9.10E-15
7 | 1.14E-05 | -6.88E-15
8 | 1.57E-06 | +2.13E-14
9 | 2.18E-07 | +1.38E-14
10 | 5.39E-08 | +8.88E-16
11 | 1.26E-08 | +1.78E-15
12 | 2.47E-09 | -2.22E-15
13 | 3.63E-10 | +5.55E-15
Timing: incl. excl.
------------------------------------------------------------
Dyson equation: 0.793 0.793 18.8% |-------|
Hybridization compression (AAA): 0.117 0.117 2.8% ||
Pseudo-particle self-energy: 2.640 2.640 62.6% |------------------------|
Single particle Green's function: 0.193 0.193 4.6% |-|
Other: 0.472 0.472 11.2% |---|
------------------------------------------------------------
Total: 4.215 100.0%
AdaPol: Hybridization fit accuracy 3.39E-07, using 5 poles.
(Order, N_Diags) = [(1, 1), (2, 10), (3, 100), (4, 1000)]
max_order = 4
iter | conv | Z-1
------+----------+-----------
1 | 3.03E-01 | +0.00E+00
2 | 5.92E-02 | -5.77E-15
3 | 9.80E-03 | +5.33E-15
4 | 1.14E-03 | +8.88E-16
5 | 1.69E-04 | -7.44E-15
6 | 3.97E-05 | +6.00E-15
7 | 9.47E-06 | -2.75E-14
8 | 1.90E-06 | -2.31E-14
9 | 3.12E-07 | -8.10E-15
10 | 3.57E-08 | -1.13E-14
11 | 5.58E-09 | +2.00E-14
12 | 1.31E-09 | -4.00E-15
13 | 3.15E-10 | -1.11E-15
Timing: incl. excl.
------------------------------------------------------------
Dyson equation: 1.172 1.172 1.2% |
Hybridization compression (AAA): 0.162 0.162 0.2% |
Pseudo-particle self-energy: 88.449 88.449 93.8% |-------------------------------------|
Single particle Green's function: 3.736 3.736 4.0% |-|
Other: 0.788 0.788 0.8% |
------------------------------------------------------------
Total: 94.308 100.0%
Visualization
To see the convergence with expansion order we plot the resulting single-particle Green’s function.
[6]:
from triqs.plot.mpl_interface import oplot, oplotr, oploti, plt
def plot_dlr_imtime(g_tau, label, n_tau=400, marker='x', linestyle='-', color=None):
from triqs.gf import make_gf_imtime
g_tau_fine = make_gf_imtime(g_tau, n_tau=n_tau)
if color is None:
color = plt.plot([], [], linestyle+marker, label=label)[0].get_color()
label = None
oplotr(g_tau, marker=marker, label=None, color=color)
oplotr(g_tau_fine, label=label, color=color, linestyle=linestyle)
[7]:
import glob
filenames = np.sort(glob.glob('data_order_*.h5'))
results = []
for filename in filenames:
print(f'--> Loading: {filename}')
with HDFArchive(filename, 'r') as A:
results.append(A['S'])
--> Loading: data_order_1.h5
--> Loading: data_order_2.h5
--> Loading: data_order_3.h5
--> Loading: data_order_4.h5
[8]:
for r in results:
plot_dlr_imtime(r.G_tau['0'], label=f'Order {r.order}', linestyle='-' if r.order != 2 else ':')
delta_tau = results[0].Delta_tau['0']
plot_dlr_imtime(delta_tau, label='Exact', marker='', linestyle='-.', color='black')
plt.ylabel(r'$G(\tau)$'); plt.xlabel(r'$\tau$')
plt.legend(loc='best'); plt.grid(True); plt.ylim(top=0);
[9]:
for r in results:
diff_tau = delta_tau - r.G_tau['0']
diff_tau.data[:] = np.abs(diff_tau.data)
plot_dlr_imtime(diff_tau, label=f'Order {r.order}',
linestyle='-' if r.order != 2 else ':')
plt.semilogy([], [])
plt.legend(loc='best')
plt.ylabel(r'max$|G(\tau) - G_{exact}(\tau)|$')
plt.xlabel(r'$\tau$')
plt.grid(True)
The convergence rate with expansion order is fast and the error at order 4 is about \(10^{-3}\).