Skip to content

Commit 0b49571

Browse files
HugoStrandHugoStrand
committed
[bench] dimer example
1 parent d2262c7 commit 0b49571

2 files changed

Lines changed: 227 additions & 0 deletions

File tree

‎benchmarks/dimer/calc_ppgf.py‎

Lines changed: 160 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,160 @@
1+
2+
3+
import numpy as np
4+
5+
from triqs.gfs import BlockGf, Gf, MeshDLRImTime, MeshImTime, make_gf_imtime, make_gf_dlr
6+
from triqs.operators import c, c_dag, Operator
7+
8+
from triqs.utility import mpi
9+
from triqs_xca.block_sparse_solver import BlockSparseSolver
10+
11+
12+
beta = 1.0
13+
14+
# -- Full dimer calculation
15+
16+
t = -np.sqrt(2.0)
17+
18+
H = t * ( c_dag('0', 0) * c('0', 1) + c_dag('0', 1) * c('0', 0) )
19+
20+
fundamental_operators = [c('0', 0), c('0', 1)]
21+
22+
Sd = BlockSparseSolver(
23+
H_loc=H, beta=beta, w_max=4.0, eps=1e-10,
24+
gf_struct=[['0', 2]], conserved_operators=[])
25+
Sd.solve(max_order=0, spgf_max_order=1, maxiter=0, verbose=True)
26+
27+
rho_d = Sd.many_body_density_matrix(Sd.G0)[0][1]
28+
print(f'rho_d =\n{rho_d.real}')
29+
30+
G_tau_ref = Sd.G_tau['0'][0, 0]
31+
32+
33+
# -- Trace out one site
34+
35+
# The reduced ppgf is given by the diagonal trace over the bath
36+
# up to an additional normalization Tr[G(beta)] = -1
37+
38+
G_00 = Sd.G0['0'][0, 0]
39+
G_01 = Sd.G0['0'][1, 1]
40+
41+
#G_ref = S.G.copy()
42+
G_ref = BlockGf(name_list=['0'], block_list=[Gf(mesh=Sd.G.mesh, target_shape=[2, 2])])
43+
G_ref['0'][0, 0] = 0.5 * (G_00 + G_01)
44+
G_ref['0'][1, 1] = 0.5 * (G_00 + G_01)
45+
46+
# Normalize G_test
47+
rho_ref = -make_gf_dlr(G_ref['0'])(beta)
48+
print(f'rho_ref =\n{rho_ref}')
49+
50+
eta = np.log(np.trace(rho_ref)) / beta
51+
tau = np.array([float(t) for t in G_ref.mesh])
52+
53+
G_ref_norm = G_ref.copy()
54+
G_ref_norm['0'].data[:] *= np.exp(- eta * tau)[:, None, None]
55+
rho_ref_norm = -make_gf_dlr(G_ref_norm['0'])(beta)
56+
57+
print(f'rho_ref_norm =\n{rho_ref_norm}')
58+
#exit()
59+
60+
G_ref = G_ref_norm
61+
62+
63+
# -- Perturbative calculation -- Expanding in one site
64+
65+
H0 = 0. * c_dag('0', 0) * c('0', 0)
66+
S = BlockSparseSolver(H_loc=H0, beta=beta, w_max=4.0, eps=1e-10,
67+
gf_struct=[['0', 1]], conserved_operators=[])
68+
S.Delta_tau['0'].data[:] = -0.5 * t**2
69+
S.solve(max_order=1, verbose=True)
70+
rho = S.many_body_density_matrix(S.G)[0][1]
71+
print(f'rho =\n{rho.real}')
72+
73+
74+
from triqs.plot.mpl_interface import oplot, plt
75+
76+
plt.figure(figsize=(3.25*2, 8))
77+
78+
subp = [2, 2, 1]
79+
80+
plt.subplot(*subp); subp[-1] += 1
81+
oplot(make_gf_imtime(G_tau_ref, n_tau=100).real, label='exact')
82+
oplot(make_gf_imtime(S.G_tau, n_tau=100).real, ':', label='o1')
83+
plt.ylim(top=0)
84+
85+
plt.subplot(*subp); subp[-1] += 1
86+
oplot(make_gf_imtime(G_tau_ref - S.G_tau['0'][0, 0], n_tau=100).real, label='exact')
87+
88+
#plt.subplot(*subp); subp[-1] += 1
89+
#oplot(make_gf_imtime(Sd.G0, n_tau=100).real, label=None)
90+
91+
plt.subplot(*subp); subp[-1] += 1
92+
#oplot(make_gf_imtime(S.G0, n_tau=100).real, label='G0')
93+
oplot(make_gf_imtime(S.G['0'][0, 0], n_tau=100).real, label='G')
94+
oplot(make_gf_imtime(G_ref['0'][0, 0], n_tau=20).real, '.', label='G (ref)')
95+
96+
plt.subplot(*subp); subp[-1] += 1
97+
oplot(make_gf_imtime((S.G - G_ref)['0'][0, 0], n_tau=100).real, label=None)
98+
plt.ylabel('Error ppgf')
99+
100+
if False:
101+
subp = [3, 1, 1]
102+
for bidx in ['0', '1', '2']:
103+
plt.subplot(*subp); subp[-1] += 1
104+
oplot(make_gf_imtime(Sd.G0[bidx], n_tau=100).real)
105+
plt.title(bidx)
106+
plt.tight_layout()
107+
plt.show()
108+
109+
exit()
110+
111+
112+
113+
ed = TriqsExactDiagonalization(H, fundamental_operators, beta)
114+
115+
#m = MeshDLRImTime(beta=beta, statistic='Fermion', eps=1e-12, w_max=2.0)
116+
m = MeshImTime(beta=beta, statistic='Fermion', n_tau=400)
117+
118+
G_tau = Gf(mesh=m, target_shape=[])
119+
ed.set_g2_tau(G_tau, c('0', 0), c_dag('0', 0))
120+
121+
# -- XCA approximation
122+
123+
H0 = 0. * c_dag('0', 0) * c('0', 0)
124+
125+
Ss = []
126+
127+
orders = [1]
128+
for order in orders:
129+
S = BlockSparseSolver(H_loc=H0, beta=beta, w_max=4.0, eps=1e-10, gf_struct=[['0', 1]])
130+
S.Delta_tau['0'].data[:] = -0.5 * t**2
131+
132+
S.solve_bare(max_order=order, use_dyson=False, verbose=True)
133+
S.G_tau_bare = S.G_tau.copy()
134+
S.eta = 0
135+
136+
S.solve_bare(max_order=order, use_dyson=True, verbose=True)
137+
S.G_tau_bare_dyson = S.G_tau.copy()
138+
139+
S.G = S.G0.copy() # reset
140+
S.eta = 0.0 # reset
141+
S.solve(max_order=order, tol=1e-8, maxiter=20, verbose=True, mix=1.)
142+
Ss.append(S)
143+
144+
from triqs.plot.mpl_interface import oplot, plt, oplotr, oploti
145+
146+
if mpi.is_master_node():
147+
for S in Ss:
148+
if S.max_order % 2 == 1:
149+
s = '-'
150+
else:
151+
s = ':'
152+
153+
n_tau = 40
154+
oplot(make_gf_imtime(S.G_tau_bare, n_tau=n_tau).real, '.'+s, label=f'O{S.max_order} bare')
155+
oplot(make_gf_imtime(S.G_tau_bare_dyson, n_tau=n_tau).real, s+'x', label=f'O{S.max_order} bare Dys')
156+
oplot(make_gf_imtime(S.G_tau, n_tau=n_tau).real, s+'+', label=f'O{S.max_order} sc')
157+
158+
oplot(G_tau.real, '--', lw=4., alpha=0.5, label='ED')
159+
#plt.ylim(top=0.)
160+
plt.show()

‎benchmarks/dimer/calc_spgf.py‎

Lines changed: 67 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,67 @@
1+
2+
3+
import numpy as np
4+
5+
from triqs.gfs import Gf, MeshDLRImTime, MeshImTime, make_gf_imtime
6+
from triqs.operators import c, c_dag, Operator
7+
8+
from pyed.TriqsExactDiagonalization import TriqsExactDiagonalization
9+
10+
from triqs.utility import mpi
11+
from triqs_xca.block_sparse_solver import BlockSparseSolver
12+
13+
t = -np.sqrt(2.0)
14+
H = t * ( c_dag('0', 0) * c('0', 1) + c_dag('0', 1) * c('0', 0) )
15+
16+
fundamental_operators = [c('0', 0), c('0', 1)]
17+
18+
beta = 4.0
19+
20+
ed = TriqsExactDiagonalization(H, fundamental_operators, beta)
21+
22+
#m = MeshDLRImTime(beta=beta, statistic='Fermion', eps=1e-12, w_max=2.0)
23+
m = MeshImTime(beta=beta, statistic='Fermion', n_tau=400)
24+
25+
G_tau = Gf(mesh=m, target_shape=[])
26+
ed.set_g2_tau(G_tau, c('0', 0), c_dag('0', 0))
27+
28+
# -- XCA approximation
29+
30+
H0 = 0. * c_dag('0', 0) * c('0', 0)
31+
32+
Ss = []
33+
34+
orders = [1]
35+
for order in orders:
36+
S = BlockSparseSolver(H_loc=H0, beta=beta, w_max=4.0, eps=1e-10, gf_struct=[['0', 1]])
37+
S.Delta_tau['0'].data[:] = -0.5 * t**2
38+
39+
S.solve_bare(max_order=order, use_dyson=False, verbose=True)
40+
S.G_tau_bare = S.G_tau.copy()
41+
S.eta = 0
42+
43+
S.solve_bare(max_order=order, use_dyson=True, verbose=True)
44+
S.G_tau_bare_dyson = S.G_tau.copy()
45+
46+
S.G = S.G0.copy() # reset
47+
S.eta = 0.0 # reset
48+
S.solve(max_order=order, tol=1e-8, maxiter=20, verbose=True, mix=1.)
49+
Ss.append(S)
50+
51+
from triqs.plot.mpl_interface import oplot, plt, oplotr, oploti
52+
53+
if mpi.is_master_node():
54+
for S in Ss:
55+
if S.max_order % 2 == 1:
56+
s = '-'
57+
else:
58+
s = ':'
59+
60+
n_tau = 40
61+
oplot(make_gf_imtime(S.G_tau_bare, n_tau=n_tau).real, '.'+s, label=f'O{S.max_order} bare')
62+
oplot(make_gf_imtime(S.G_tau_bare_dyson, n_tau=n_tau).real, s+'x', label=f'O{S.max_order} bare Dys')
63+
oplot(make_gf_imtime(S.G_tau, n_tau=n_tau).real, s+'+', label=f'O{S.max_order} sc')
64+
65+
oplot(G_tau.real, '--', lw=4., alpha=0.5, label='ED')
66+
#plt.ylim(top=0.)
67+
plt.show()

0 commit comments

Comments
 (0)