Initial import of required modules and outsourced parts#
[2]:
import os
import numpy as np
from typing import Dict, cast
import scce.loop
from scce.cluster_bath_solver.DMRG import ClusterBathSolver_DMRG_sfo_matrix # only for mypy
from scce import save_result, LoopSCCE
from lattice_builder.kpoints import bands
import matplotlib.pyplot as plt
import matplotlib as matplotlib
from matplotlib.colors import ListedColormap
# Outsourced reusable parts
from base_parameters import spr, cpr
from compute_bounds import get_lower_bound, get_upper_bound
from compute_chebyshev import compute_chebyshev
from self_energy import SigmaCluster_U
plt.style.use('HQStyle')
matplotlib.rcParams['figure.dpi'] = 200
matplotlib.rcParams['figure.figsize'] = [3.15, 2.362]
matplotlib.rcParams.update({'font.size': 10})
matplotlib.rcParams.update({'legend.fontsize': 10})
mutedmagma = ListedColormap(np.loadtxt('/home/bschoenauer/.config/matplotlib/mutedmagma.csv',
delimiter=','))
np.set_printoptions(linewidth=np.inf)
Update of calculation parameters#
[ ]:
skipDMRG = False
Plot = True
spr.lattice = "1D_CHAIN"
cpr.ChebyOrder = 200
cpr.Use_EOM = True
cpr.Use_EOM2 = False
spr.UU = 5.0
spr.UV = 0.0
spr.NumCluster = 64
spr.Mx = 16
spr.phs_symmetric = False
spr.update()
DMRG_dirResults = "/dmrg_results/" + spr.lattice + "_U" + str(spr.UU) + "_" + str(spr.UV)
os.makedirs(DMRG_dirResults, exist_ok=True)
# Number of omega points
Omegas = np.linspace(-5.0, 5.0, 1001)
Nu = (spr.N + spr.Sz) // 2
Nd = (spr.N - spr.Sz) // 2
eta_scale = float(os.environ.get('eta_scale', 1.25))
System-specific parameters and system description#
[ ]:
unitcell = {'atoms': [['0', 'A', [0, 0, 0], -spr.UV, spr.UU, spr.Delta_S], ],
'lattice_vectors': [[1, 0, 0],
[0, 0, 0],
[0, 0, 0]],
'bonds': [['0', '0', [1, 0, 0], [-spr.tu, -spr.td, spr.tud, spr.tdu],
[spr.UV, spr.Jz, spr.Jp],
[0.0, 0.0, spr.Delta_Sud, spr.Delta_Sdu],
'Hx']]
}
# specification of the system (cluster) for the SCCE-DMRG run
spr.N = spr.MC * len(unitcell.get('atoms'))
cluster = {'cluster_size': [spr.Mx, spr.My, 1],
'system_size': [1, 1, 1],
'system_boundary_condition': 'hard-wall',
'site_type': 'spinful_fermions',
'N': spr.N,
'Sz': spr.Sz,
'mod_N': spr.mod_N,
'mod_Sz': spr.mod_Sz
}
spr.MC = spr.MC * len(unitcell.get('atoms'))
dict_optional: Dict[str, int] = {}
Initialization of the ClusterBathSolver#
[ ]:
dict_run = {'states': 1}
indict = {'unitcell': unitcell,
'system': cluster,
'run': dict_run,
'optional': dict_optional}
loop_dir = os.path.dirname(scce.loop.__file__)
binary_path = str(os.path.join(loop_dir, '../dmrg/SuWaves-mt'))
script_path = str(os.path.join(loop_dir, '../dmrg/doDMRG_SCCE'))
loop = LoopSCCE(indict,
MF='iHF',
run=dict_run,
optional=dict_optional,
DMRG_dirResults=DMRG_dirResults,
DMRG_dirScratch="/dmrg_scratch/SCCE",
DMRG_script=script_path,
DMRG_binary=binary_path)
if loop.h0.ndim == 2:
loop.h0[0, 0] += spr.UV
loop.h0[spr.MC - 1, spr.MC - 1] += spr.UV
else:
loop.h0[0, 0, 0] += spr.UV
loop.h0[0, spr.MC - 1, spr.MC - 1] += spr.UV
loop.h0[1, 0, 0] += spr.UV
loop.h0[1, spr.MC - 1, spr.MC - 1] += spr.UV
Prepare Chebyshev expansion#
[ ]:
loop.dbg_plot = False
loop(skipCalc=skipDMRG,
DMRG_Scale_AB=20,
DMRG_NumThreads=6,
NumStates=6
)
if loop.solver is None:
raise Exception("DMRG solver not available")
if not isinstance(loop.solver, ClusterBathSolver_DMRG_sfo_matrix):
raise Exception("Wrong DMRG solver selected.")
E0 = loop.En[0]
Dn = loop.En[:] - E0
Zn = np.exp(-1.0 * loop.beta_mix * Dn)
Z0 = cast(float, np.sum(Zn))
cpr.En = np.array([])
cpr.Zn = np.array([])
for en, zn in zip(loop.En, Zn):
if abs(en - E0) < 1. / loop.beta_mix:
cpr.En = np.append(cpr.En, en)
cpr.Zn = np.append(cpr.Zn, zn)
save_result(loop, filedir='./', filename="Res_EL_" + spr.lattice, overwrite=True)
Calculate extended spectrum for scaling to (-1, 1)#
[ ]:
cpr.ELp, cpr.ELm = get_lower_bound(loop, spr, cpr, skipDMRG)
cpr.EHp, cpr.EHm = get_upper_bound(loop, spr, cpr, skipDMRG)
loop.N = spr.N
loop.Sz = spr.Sz
Chebyshev expansion#
[ ]:
calc_indices = [n for n in range((spr.MC + 1) // 2)]
symmetries = [lambda x: x, lambda x: (spr.MC - 1 - x)]
cpr.get_scaling()
if not os.path.exists('./Res'):
os.makedirs('Res/')
P_chebyshev_d, P_chebyshev_Fd = compute_chebyshev(loop,
calc_indices,
symmetries,
spr,
cpr,
skipDMRG)
Determine k-Point path#
[ ]:
lat_vectors = loop.lb.config['unitcell']['lattice_vectors']
if len(lat_vectors) == 1 or (abs(lat_vectors[1]) == 0 and abs(lat_vectors[2] == 0)):
r1 = np.array(lat_vectors[0])
kDict, kPath, primitiveVectors, *_ = bands.returnLatticeAndPath1D(r1)
kPath = 'X1-' + kPath
totalSamplePoints = spr.M
elif len(lat_vectors) == 2 or abs(lat_vectors[2] == 0):
r1 = np.array(lat_vectors[0])
r2 = np.array(lat_vectors[1])
kDict, kPath, primitiveVectors, *_ = bands.returnLatticeAndPath2D(r1, r2)
totalSamplePoints = spr.M
elif len(lat_vectors) == 3:
r1 = np.array(lat_vectors[0])
r2 = np.array(lat_vectors[1])
r3 = np.array(lat_vectors[2])
kDict, kPath, primitiveVectors, *_ = bands.returnLatticeAndPath3D(r1, r2, r3)
totalSamplePoints = spr.M
else:
print("Dimensionality of the system unclear!")
(kPathPoints,
plotHelper,
locationSpecialPoints,
indexSpecialPoints) = bands.constructPath(kDict,
kPath,
primitiveVectors,
totalSamplePoints=totalSamplePoints)
kLabels = kPath.split('-')
Perform CPT (DCA) and plot the Spectral function#
[ ]:
PlotFile = "Spectral_Function_" + spr.lattice + "_M" + str(spr.MC) + "_" + str(spr.NumCluster)
PlotFile += "_U" + str(spr.UU)
if cpr.Use_EOM:
PlotFile += "_EOM1"
if cpr.Use_EOM2:
PlotFile += "_EOM2"
PlotFile += "_P" + str(cpr.ChebyOrder) + ".png"
# Calculate the cluster self-energy, which is needed as the input for "spectral"
SigmaC = np.zeros((Omegas.size, spr.MC, spr.MC), dtype=complex)
for n in range(Omegas.size):
SigmaC[n, :, :] = SigmaCluster_U(Omegas[n],
spr.UU,
eta_scale,
loop,
spr,
cpr,
P_chebyshev_d,
P_chebyshev_Fd)
A, _ = loop.lb.spectral(energyWindow=Omegas, kPoints=kPathPoints, Sigma_Array=SigmaC, eta=1e-1)
kLabels = kPath.split('-')
kLabels[1] = 'Γ'
fig = plt.figure(dpi=200)
ax = fig.add_subplot(111)
plt.imshow(A, aspect='auto', origin='lower', interpolation='lanczos', # interpolation='lanczos'
extent=(0, 1, Omegas[0], Omegas[-1]), cmap=mutedmagma)
for ii, _ in enumerate(kLabels):
plt.axvline(x=indexSpecialPoints[ii] / indexSpecialPoints[-1], ymin=0, ymax=1,
ls='--', c='#B4B4B4', alpha=0.6)
plt.xticks(indexSpecialPoints / indexSpecialPoints[-1], kLabels)
plt.ylabel('ω / t', fontsize=9)
plt.xlabel('k-point path', fontsize=9)
plt.tight_layout()
plt.clim(0, )
plt.colorbar()
plt.savefig(PlotFile)
plt.show()