  The nature of correlated states in twisted bilayer graphene (TBG) at the
magic angle has received intense attention in recent years. We present a
numerical study of an interacting Bistritzer-MacDonald (IBM) model of TBG using
a suite of methods in quantum chemistry, including Hartree-Fock, coupled
cluster singles, doubles (CCSD), and perturbative triples (CCSD(T)), as well as
a quantum chemistry formulation of the density matrix renormalization group
method (DMRG). Our treatment of TBG is agnostic to gauge choices, and hence we
present a new gauge-invariant formulation to detect the spontaneous symmetry
breaking in interacting models. To benchmark our approach, we focus on a
simplified spinless, valleyless IBM model. At integer filling ($\nu=0$), all
numerical methods agree in terms of energy and $C_{2z} \mathcal{T}$ symmetry
breaking. Additionally, as part of our benchmarking, we explore the impact of
different schemes for removing ``double-counting'' in the IBM model. Our
results at integer filling suggest that cross-validation of different IBM
models may be needed for future studies of the TBG system. After benchmarking
our approach at integer filling, we perform the first systematic study of the
IBM model near integer filling (for $|\nu|< 0.2$). In this regime, we find that
the ground state can be in a metallic and $C_{2z} \mathcal{T}$ symmetry
breaking phase. The ground state appears to have low entropy, and therefore can
be relatively well approximated by a single Slater determinant. Furthermore, we
observe many low entropy states with energies very close to the ground state
energy in the near integer filling regime.
