Hubbard model Green's function
Author Martin Ulaga
Uses the Lanczos algorithm[1] to calculate the Green's function of the Hubbard model. See also the documentation page for the spin structure factor. The Green's function is given by
\[
G({\bf k}, \omega)=-i\int dt e^{-i\omega t}\langle \lbrace c_{\bf k}(t),c^{\dagger}_{\bf k}\rbrace\rangle.
\]
To achieve a mesh of momentum space, we can use either use the generated momenta inside the Wigner-Seitz cell allowed by the finite cluster or twisted boundary conditions[2], where the lattice \({\bf k}\)-point is shifted to \({\bf k + \boldsymbol{\theta}}\) by introducing a flux via Peierls substitution in the Hamiltonian, \(t_{ij} \rightarrow t_{ij}\exp(i{\bf \boldsymbol{\theta}\cdot r}_{ij})\). For this example, we're calculating the spectral function along the cut \(\Gamma\)-\(M\), i.e., the diagonal of the reciprocal lattice unit cell.
Example code
using XDiag
using HDF5
function main()
say_hello()
# IO
latticeInput = "../../misc/data/square.8.hubbard.ttprime.toml"
lfile = FileToml(latticeInput)
filename = "../../misc/data/examples_output/hubbard_greens_f.h5"
outfile = h5open(filename, "w")
# Define the Hubbard model
N = 8
nup = 3
ndn = 3
ops = read_opsum(lfile, "Interactions")
t = 1.0
tp = 0.0
ops["Ty"] = t
ops["Tx"] = t
ops["Tprime+"] = tp
ops["Tprime-"] = tp
ops += "U" * Op("HubbardU")
ops["U"] = -6.0
@show(ops)
opsTBC = ops
irrep = read_representation(lfile, "Gamma.C1.A")
# compute groundstate (known to be at k=0)
println("Computing ground state ...")
block = Electron(N, nup, ndn, irrep)
e0, gs = eig0(ops, block)
make_complex!(gs)
println("done.")
println("Ground state energy: $e0")
outfile["e0"] = e0
function tbc_mesh(kx, ky)
# compare to twisted boundary condition calculation,
# cf. Zemljic & Prelovsek, PRB 75, 104514 (2007)
# Tohyama, PRB 70, 174517 (2004).
println("\n considering TBC momentum [$kx,$ky]")
opsTBC["Tx"] = t * exp(1im * kx)
opsTBC["Ty"] = t * exp(1im * ky)
opsTBC["Tprime+"] =
tp * exp(1im * (kx + ky)) # r_ij = (1,1)
opsTBC["Tprime-"] =
tp * exp(1im * (kx - ky)) # r_ij = (1,-1)
e0, gs_tbc = eig0(opsTBC, block)
# @show e0
c_q_tbc = symmetrize(sqrt(N) * Op("Cup", 1), irrep)
Av = apply(c_q_tbc, gs_tbc)
@show nrm = norm(Av)
Av /= nrm
res_tbc = eigvals_lanczos_inplace(opsTBC, Av)
@show res_tbc.eigenvalues[1]
return nrm, res_tbc
end
# loop through momenta
irreps = Dict{String,Vector{Float64}}(
"Gamma.C1.A" => [0.0, 0.0],
"M.C1.A" => [3.1415926535897931, 3.1415926535897931],
"Sigma0.C1.A" => [1.5707963267948966, 1.5707963267948966],
"Sigma1.C1.A" => [1.5707963267948966, -1.5707963267948966],
"Sigma2.C1.A" => [-1.5707963267948966, 1.5707963267948966],
"Sigma3.C1.A" => [-1.5707963267948966, -1.5707963267948966],
"X0.C1.A" => [3.1415926535897931, 0.0000000000000000],
"X1.C1.A" => [0.0000000000000000, 3.1415926535897931])
for (name, momentum) in irreps
println("\n considering irrep $name, momentum $momentum")
aq_irrep = read_representation(lfile, name)
c_q = symmetrize(sqrt(N) * Op("Cup", 1), aq_irrep)
Av = apply(c_q, gs)
@show nrm = norm(Av)
Av /= nrm
res = eigvals_lanczos_inplace(ops, Av)
@show res.eigenvalues[1]
outfile["$name/norm"] = nrm
outfile["$name/alphas"] = res.alphas
outfile["$name/betas"] = res.betas
outfile["$name/eigs"] = res.eigenvalues
tbc_mesh(momentum[1], momentum[2])
end
for (k, kx) in enumerate(0:pi/101:pi)
nrm, res_tbc = tbc_mesh(kx, kx)
ik = "$(k-1)" # mimics c++ output
# group = create_group(outfile, ik)
# subgroup = create_group(group, ik)
outfile["$ik/$ik/norm"] = nrm
outfile["$ik/$ik/alphas"] = res_tbc.alphas
outfile["$ik/$ik/betas"] = res_tbc.betas
outfile["$ik/$ik/eigs"] = res_tbc.eigenvalues
end
end
main()
#include <xdiag/all.hpp>
using namespace xdiag;
using namespace arma;
using fmt::format;
using namespace std::complex_literals;
int main() try {
say_hello();
// IO
std::string latticeInput =
XDIAG_DIRECTORY "/misc/data/square.8.hubbard.ttprime.toml";
auto lfile = FileToml(latticeInput);
std::string filename =
XDIAG_DIRECTORY "/misc/data/examples_output/hubbard_greens_f.h5";
auto outfile = FileH5(filename, "w!");
// Define the Hubabrd model
int N = 8;
int nup = 3;
int ndn = 3;
auto ops = read_opsum(lfile, "Interactions");
auto t = -1.;
auto tp = 0.;
ops["Ty"] = t;
ops["Tx"] = t;
ops["Tprime+"] = tp;
ops["Tprime-"] = tp;
ops += "U" * Op("HubbardU");
ops["U"] = -6.;
auto opsTBC = ops;
auto irrep = read_representation(lfile, "Gamma.C1.A");
// compute groundstate (known to be at k=0)
Log("Computing ground state ...");
auto block = Electron(N, nup, ndn, irrep);
auto [e0, gs] = eig0(ops, block);
gs.make_complex();
Log("Ground state energy: {:.12f}", e0);
outfile["e0"] = e0;
auto tbc_mesh =
[&opsTBC, &t, &tp, &block, &N,
&irrep](double const &kx,
double const &ky) -> std::pair<double, EigvalsLanczosResult> {
// lambda to perform TBC calc.
// compare to twisted boundary condition calculation,
// cf. Zemljic & Prelovsek, PRB 75, 104514 (2007);
// Tohyama, PRB 70, 174517 (2004).
Log("TBC calculation, momentum [{}, {}]", kx, ky);
opsTBC["Tx"] = t * exp(1i * kx);
opsTBC["Ty"] = t * exp(1i * ky);
opsTBC["Tprime+"] = tp * exp(1i * (kx + ky)); // r_ij = (1,1)
opsTBC["Tprime-"] = tp * exp(1i * (kx - ky)); // r_ij = (1,-1)
// XDIAG_SHOW(opsTBC);
auto [e0_tbc, gs_tbc] = eig0(opsTBC, block, 1e-12, 100);
// XDIAG_SHOW(e0_tbc);
auto c_q_tbc = symmetrize(sqrt((double)N) * Op("Cup", 0), irrep);
auto Av = apply(c_q_tbc, gs_tbc);
auto nrm = norm(Av);
XDIAG_SHOW(nrm);
Av /= nrm;
auto res_tbc = eigvals_lanczos_inplace(opsTBC, Av, 1, 1e-12, 100);
XDIAG_SHOW(res_tbc.eigenvalues(0));
Log("n iterations {}", res_tbc.niterations);
return {nrm, res_tbc};
};
// loop through momenta of the cluster reciprocal lattice
std::vector<std::pair<std::string, arma::vec>> irreps = {
{"Gamma.C1.A", {0., 0.}},
{"M.C1.A", {3.1415926535897931, 3.1415926535897931}},
{"Sigma0.C1.A", {1.5707963267948966, 1.5707963267948966}},
{"Sigma1.C1.A", {1.5707963267948966, -1.5707963267948966}},
{"Sigma2.C1.A", {-1.5707963267948966, 1.5707963267948966}},
{"Sigma3.C1.A", {-1.5707963267948966, -1.5707963267948966}},
{"X0.C1.A", {3.1415926535897931, 0.0000000000000000}},
{"X1.C1.A", {0.0000000000000000, 3.1415926535897931}}};
for (const auto &[name, momentum] : irreps) {
Log("considering irrep {}, momentum [{}, {}]", name, momentum(0),
momentum(1));
auto aq_irrep = read_representation(lfile, name);
auto c_q = symmetrize(sqrt((double)N) * Op("Cup", 0),
aq_irrep); // normalize to 1/sqrt(N)
auto Av = apply(c_q, gs);
auto nrm = norm(Av);
Av /= nrm;
XDIAG_SHOW(nrm);
auto res = eigvals_lanczos_inplace(ops, Av);
XDIAG_SHOW(res.eigenvalues(0));
outfile[format("{}/norm", name)] = nrm;
outfile[format("{}/alphas", name)] = res.alphas;
outfile[format("{}/betas", name)] = res.betas;
outfile[format("{}/eigs", name)] = res.eigenvalues;
// compare to tbc calc
auto [nrm_tbc, res_tbc] = tbc_mesh(momentum(0), momentum(1));
XDIAG_SHOW(nrm_tbc);
XDIAG_SHOW(res_tbc.eigenvalues(0));
}
// TBC mesh, 10x10 grid
int nkx = 101;
int nky = nkx;
auto momentum = [&nkx, &nky](int const &ikx,
int const &iky) -> std::pair<double, double> {
return {(double)ikx * 3.1415926535897931 / (double)(nkx - 1),
(double)iky * 3.1415926535897931 / (double)(nky - 1)};
};
for (int ikx = 0; ikx < nkx; ++ikx) {
auto iky = ikx;
// for (int iky = 0; iky < nky; ++iky) {
auto [kx, ky] = momentum(ikx, iky);
auto [nrm, res] = tbc_mesh(kx, ky);
XDIAG_SHOW(nrm);
XDIAG_SHOW(res.eigenvalues(0));
outfile[format("{}/{}/norm", ikx, iky)] = nrm;
outfile[format("{}/{}/alphas", ikx, iky)] = res.alphas;
outfile[format("{}/{}/betas", ikx, iky)] = res.betas;
outfile[format("{}/{}/eigs", ikx, iky)] = res.eigenvalues;
outfile[format("{}/{}/kx", ikx, iky)] = kx;
outfile[format("{}/{}/ky", ikx, iky)] = ky;
// }
}
} catch (Error e) {
error_trace(e);
}
References
[1] Prelovšek, P., & Bonča, J. (2013). Ground state and finite temperature Lanczos methods. Strongly Correlated Systems: Numerical Methods, 1-30.
[2] Tohyama, T. (2004). Asymmetry of the electronic states in hole-and electron-doped cuprates: Exact diagonalization study of the \(t-t′-t ″-J\) model. Phys. Rev. B, 70(17), 174517