# Generalizing the self-healing diffusion Monte Carlo approach to finite
temperature: a path for the optimization of low-energy many-body bases

Fernando A. Reboredo  Jeongnim Kim  Materials Science and Technology Division,
Oak Ridge National Laboratory, Oak Ridge, TN 37831, USA

###### Abstract

A statistical method is derived for the calculation of thermodynamic
properties of many-body systems at low temperatures. This method is based on
the self-healing diffusion Monte Carlo method for complex functions [F. A.
Reboredo J. Chem. Phys. 136, 204101 (2012)] and some ideas of the correlation
function Monte Carlo approach [D. M. Ceperley and B. Bernu, J. Chem. Phys. 89,
6316 (1988)]. In order to allow the evolution in imaginary time to describe
the density matrix, we remove the fixed-node restriction using complex
antisymmetric guiding wave functions. In the process we obtain a parallel
algorithm that optimizes a small subspace of the many-body Hilbert space to
have maximum overlap with the subspace spanned by the lowest-energy
eigenstates of a many-body Hamiltonian. We show in a model system that the
partition function is progressively maximized within this subspace. We show
that the subspace spanned by the small basis systematically converges towards
the subspace spanned by the lowest energy eigenstates. Possible applications
of this method to calculate the thermodynamic properties of many-body systems
near the ground state are discussed. The resulting basis can be also used to
accelerate the calculation of the ground or excited states with Quantum Monte
Carlo.

###### pacs:

02.70.Ss,02.70.Tt

##  I Introduction

There is a significant interest in thermodynamical properties observed as
$T\rightarrow 0$. Many physical phenomena that cover superconductivity,
magnetic and structural transitions, chemical reactions etc. require an
adequate treatment of thermal effects. These effects are crucial in systems
where there is a large number of low-energy excitations within an energy
window $1/\beta=k_{B}T$ above the ground state. Electronic thermal effects are
expected to be larger in metals and magnets than in insulators. QTS  In metals
there is a significant number of excitations with vanishing energy. The
magnetic excitations energies frequently go to zero in the long wave limit. A
significant fraction of spectroscopic techniques probe the electronic or
magnetic excitations near the ground state. The development of first-
principles techniques to obtain excitations has historically received a
significant theoretical attention. runge84 ; onida02 ; filippi09  Monte Carlo
methods used to calculate excitation energiesfilippi09  will be accelerated
with basis that retain the physics at the relevant energies.

A first-principles finite-temperature description of many-body systems is also
relevant to describe chemical reactions.mazzola12  Ionic dynamics are usually
calculated within the Born-Oppenheimer approximation. This decouples the wave
function of the “quantum” electrons from the wave function of the ions. Within
this standard approximation, electrons are at zero temperature while the ions
can move with kinetic energies that often exceed the electronic excitations.
Even within the Born-Oppenheimer-ground-state approximation, the standard
approach based on density functional theory (DFT) shows significant
differences with diffusion Monte Carlo (DMC) grossman97 ; saccani13  or
Quantum Chemistry benchmarks. At the transition saddle points, when some
chemical bonds are broken and new ones are formed, the spacing of the
corresponding electronic eigenenergies is minimum, or even zero at the conical
intersections. levine07  Electronic thermal effects are seldom included in
many-body calculations. morales10 ; mazzola12  In order to routinely include
thermal effects, significant improvements in the theory beyond the standard
approach are required.

Most ab-initio calculations in the literature of condensed matter electronic
structure are based in the ground state quantum Monte Carlo calculations of
the homogeneous electron gasceperley80  which made possible the first
approximations of DFT. hohenberg ; kohn  DFT has been extended to finite
temperature long ego. mermin65 ; kohn83 ; weinert92  Fermi occupations of
Kohn-Sham eigenstates and the addition of an entropy term have been
shownweinert92  to provide a variational density functional. However, even
nowadays, the zero temperature approximation for the exchange-correlation
potential is widely used. This approach has been long known to severely
underestimate the critical Curie temperatures of magnetic systems. statunton84
; gyorffy85 ; staunton85  Including temperature for magnetic systems is
possible for cases where the magnetic excitations can be treated classically
wang95 ; landau04  and the electrons can be assumed to be in the ground state
for constrained configurations of the spins. eisenbach09  But an adequate
description of the electronic entropy in the subspace that preserves the spin
is still lacking. yin12  Finite temperatures benchmarks of a quality
comparable to Ref. ceperley80,  are the key ingredients required to
parametrize a finite temperature density functional. Without a reliable
approximation, most work done under a DFT framework still uses a zero
temperature approximation for the exchange correlation functional.

Accurate many-body calculations at high temperatures can be performed within
the path integral Monte Carlo approach (PIMC). ceperley95  Since the cost of
PIMC diverges as $T\rightarrow 0$, it has been mainly used in the hot and
dense regime, magro96 ; militzer00 ; driver12  with a temperature comparable
to the interaction potential. An alternative approach that could start from
the zero temperature limit would be desirable.

The most accurate techniques to describe a large number of electrons
($N_{e}>1000$) at zero temperature are based in projection approaches.
ceperley80 ; purwanto09 ; booth12  One could potentially extend these methods
to finite temperature, limiting the projector $e^{-\beta\mathcal{\hat{H}}}$ to
finite $\beta$. Thermodynamical averages can be later obtained from
derivatives of the Helmholtz free energy
$F(\beta)=-1/\beta\ln[tr(e^{-\beta\mathcal{\hat{H}}})]$, where
$\mathcal{\hat{H}}$ is the many-body Hamiltonian operator and $tr(X)$ the
trace of $X$ over the complete many-body Hilbert space.

The standard diffusion Monte Carlo Method with importance sampling (DMC)
ceperley80 ; mfoulkesrmp2001 ; austin12  constrains the sign or the phase of
the wave function by imposing the nodes or the phase ortiz93  of a guiding
wave function $\Psi_{T}({\bf R})$, where the many-body coordinate ${\bf
R}=\\{{\bf r}_{1},{\bf r}_{2},\cdots,{\bf r}_{N_{e}}\\}$ is the set of
coordinates of $N_{e}$ electrons. These constraints while enforcing an anti-
symmetric fermionic wave function introduce a variational error. The quality
of the wave function and its nodes can be improved with several methods within
a variational Monte Carlo (VMC) context, mfoulkesrmp2001 ; ortiz95 ; jones97 ;
guclu05 ; umrigar07 ; toulouse08 ; petruzielo12 ; rios06  or at the DMC level.
ceperley88 ; jones97 ; keystone ; rollingstones ; stonehenge  In DMC, the
energy of the ground state is exact if the exact nodes or phase are provided.
anderson79 ; ortiz93  Improving the nodes is computationally intensive.
Avoiding this cost is key for finite temperature calculations.

In standard DMC calculations, $\mathcal{\hat{H}}$ is replaced by the fixed-
node Hamiltonian $\mathcal{\hat{H}}_{FN}$ or the fixed-phase Hamiltonian
$\mathcal{\hat{H}}_{FP}$. The use of the fixed-node or fixed phase
approximation can have undesired effects on the calculation of thermal
effects. It has been found that many fermionic systems have a ground state
with two nodal pockets. bressanini12  That is, if the ground state wave
function is real the nodal surface separates the Hilbert space in only two
pockets for positive and negative values respectively. It has been conjectured
ceperley91  that this is a general property of fermionic ground states. In the
fixed-node case, the excitations of $\mathcal{\hat{H}}_{FN}$ are forced to
share the nodes of the ground state. To be orthogonal to the fixed-node ground
state, the fixed-node excited states have to have at least an additional node.
It is known, however, that in many systems there are several fermionic excited
states near the ground state with also two nodal pockets.rasch12  Accordingly,
$tr(e^{-\beta{\mathcal{\hat{H}}}_{FN}})$ does not describe the low temperature
physics. It is easy to see that the same happens in the fixed-phase case.
Therefore, if one wishes to use a DMC-like algorithm to obtain thermodynamical
properties, one must go beyond the usual fixed-node or fixed-phase
approximations. For practical reasons, a parallel approach that can handle a
large number of excitations near the ground state would also be beneficial.

In this paper, we restart the debate on how to calculate low temperature
properties within a many-body ab-initio context taking into account recent
theoretical developments. keystone ; rockandroll ; stonehenge  A method is
derived that introduces temperature within an importance sampling procedure
that shares most of the computational tools developed for projection MC
approaches. The errors in the evolution operator $e^{-\beta\mathcal{\hat{H}}}$
resulting from the fixed-node restriction are eliminated by using complex
linear combinations of eigenstates, which do not have nodes except at the
electronic coincidental points (see Fig 1). Instead of optimizing a single
many-body wave-function so it better describes the ground state of the system,
a basis of several wave functions is optimized to maximize the overlap with
the small subspace spanned by the lowest energy eigenstates of the many-body
Hamiltonian. We show in a model system that the overlap of the optimized
subspace with the lowest energy subspace calculated with a configuration
interaction (CI) approach, increases systematically as the iterations increase
and that the partition function is maximized.

The rest of the paper is organized as follows: In Section II we describe the
general formalism; some of the formulae developed in Ref. stonehenge,  for
complex wave function is repeated here for completeness. Section III outlines
the basic algorithm. In Section IV we describe the results for a model
calculation; and finally, in Section V we discuss the possible applications
and summarize. This paper also has three appendices: A describes how to go
beyond the locality and local-time approximations; B describes how to take
advantage of the eigenstates when they are complex; Finally, C describes how
to work with eigenstate pairs to minimize the variance of the weights of the
walkers while keeping the wave function complex.

##  II A low-energy expansion of the partition function

This section extends the DMC approach ceperley80  for the calculation of the
partition function of a many-body system. We first provide background material
required to understand the rest of the paper. We generalize the upper bound
property of the energy in DMC to an upper bound property of the free energy.
We next give the general outline of our approach and describe how to avoid the
fixed-node approximation in DMC. Finally, we describe the details: basic
formulae and numerical approach.

###  II.1 The upper bound property of the truncated Helmholtz free energy

Thermal effects can be obtained by calculating all excitations within a
thermal energy window above the ground state larger than $\Delta E\sim
1/\beta=k_{B}T$ and then evaluating the density matrix reichl  as:

$\displaystyle\hat{\rho}(\beta)=$ $\displaystyle e^{-\beta\mathcal{\hat{H}}}$
(1) $\displaystyle=$ $\displaystyle\sum_{n}\left|\Psi_{n}\right\rangle
e^{-\beta E_{n}}\left\langle\Psi_{n}\right|,$

where$E_{n}<E_{0}+\Delta E$ is the eigenvalue with eigenvector
$|\Psi_{n}\rangle$ of $\mathcal{\hat{H}}$. In general $\mathcal{\hat{H}}$ is
given by

$\hat{\mathcal{H}}=\sum_{j}^{N_{e}}\frac{(\nabla_{j}+{\bf
A}_{j})^{2}}{2}+\hat{V}({\bf R})$ (2)

where ${\bf A}_{j}={\bf A}({\bf r}_{j})$ is a vector potential at point ${\bf
r}_{j}$ with magnetic field ${\bf B}({\bf r}_{j})=\nabla_{j}\times{\bf
A}_{j}$, and $\hat{V}({\bf R})$ includes the electron-electron interaction,
the interactions of the electrons spins with the magnetic field and any
external potential, local or non-local.

In a closed system that can exchange energy with a bath or reservoir
(canonical ensemble) all thermodynamical averages can be obtained using the
density matrix. The trace of the density matrix
$Z(\beta)=tr[\hat{\rho(\beta)}]$ is the partition function, whereas
$F(\beta)=-1/\beta\ln[Z(\beta)]$ is the Helmholtz free energy.

In general, $\mathcal{\hat{H}}$ has an infinite number of eigenvectors
$|\Psi_{n}\rangle$ that can be ordered with increasing eigenenergy $E_{n}$. If
$\beta(E_{n}-E_{0})\gg 1$, the contribution to $Z(\beta)$ of the eigenstate
$|\Psi_{n}\rangle$ becomes negligible. Therefore, a usual approximation is to
truncate the trace to a finite matrix with a finite number of eigenstates
$M_{S}$.

In what follows we defined $tr()$ as the trace of a truncated square matrix
with size $M_{S}$. We also relate $Z(\beta)$ and $F(\beta)$ to that truncated
trace. Since $e^{-\beta\mathcal{\hat{H}}}$ is positive definite, for a given
basis, $Z(\beta)$ increases and $F(\beta)$ decreases as $M_{S}$ increases.

The trace of any linear operator is invariant for linear transformations of
the form ${\hat{B}}e^{-\beta\mathcal{\hat{H}}}{\hat{B}}^{-1}$ with
${\hat{B}}{\hat{B}}^{-1}=1$. Thus, in principle, one does not need to obtain
the eigenstates of $e^{-\beta\mathcal{\hat{H}}}$ or equivalently
$\mathcal{\hat{H}}$ to calculate the free energy. Any linearly independent
basis that spans the same subspace can be used to obtain $Z(\beta)$.
Thermodynamical properties only require us to evaluate $Z(\beta)$ in a
linearly independent basis $\\{|\chi^{S}_{m}\rangle\\}$. However, if
statistical methods are used, then each element contributing to the truncated
trace also increases the statistical error bar. Therefore, it is
computationally more efficient to use the most compact basis, with minimum
$M_{S}$, that retains the low-energy properties.

Any eigenstate $|\Psi_{n}\rangle$ can be written in a complete basis
${|\chi^{S}_{n}}\rangle$ as

$\displaystyle|\Psi_{n}\rangle=e^{-\hat{J}}\left[\sum_{m=0}^{M_{S}-1}\lambda_{n}^{m}|\chi^{S}_{m}\rangle+\sum_{m=M_{S}}^{N_{B}\rightarrow\infty}\lambda_{n}^{m}|\chi^{S}_{m}\rangle\right],$
(3)

where $\langle{\bf R}|e^{-\hat{J}}|{\bf R}\rangle=e^{-J({\bf R})}$ is a
Jastrow factor that introduces adequate cusp conditions kato1957  and are
$\langle{\bf R}|\chi^{S}_{m}\rangle=\chi_{m}^{S}({\bf R})$ linear combinations
of an infinite orthogonal set $\\{\langle{\bf R}|n\rangle\\}$ [e.g. Slater
determinants, or Pfaffians bajdich08 , or symmetry constrained functions
(SCF), etc]. In practice we restrict the Hilbert space to a finite number
$N_{B}$. We denote the subspace spanned by $N_{B}\gg M_{S}$ functions
$|m\rangle$ as the big subspace. The big subspace has to be large enough to
describe the low temperature physics of the complete Hilbert space, which is
in general infinite. We define the small subspace as the subspace spanned by
the first $M_{S}$ basis functions $|\chi^{S}_{n}\rangle$.

Within the small basis $\\{|\chi^{S}_{m}\rangle\\}$, the free energy will be
minimum if all $|\Psi_{n<M_{S}}\rangle$ can be spanned in the small basis,
namely, $\lambda_{n\geq M_{S}}^{m}=\lambda_{n}^{m\geq M_{S}}=0$. Errors in the
small basis will result in higher values of the free energy. Thus the free
energy in the truncated basis is an upper bound to the true Helmholtz free
energy. Optimization of the Helmholtz free energy in the small basis is
analogous to the variational principle of the ground state. Likewise, the
partition function in the truncated basis is a lower bound of the exact
partition function. Improved bounds may be obtained with a basis that better
describes the lower energy eigenstates. $M_{S}$ has to be large enough to
include all the relevant physics for a given temperature.

Most of the optimization methods in the QMC literature focus, on optimizing
the eigenstates. Several methods have been proposed to obtain low-energy
excited states within the linear method Monte Carlofilippi09 ; toulouse12
(LMMC) or diffusion Monte Carlo. ceperley88 ; rockandroll  However, since the
eigenstates are sometimes difficult calculate, and we only need an average, we
argue that one might save computational time by optimizing the many-body basis
first as in the correlation function diffusion Monte Carlo (CFDMC) ceperley88
method. Optimizing the basis directly could be more practical than obtaining
accurate eigenstates energies, if the number of excitations near the ground
state is large (e. g. typically the case in metallic or magnetic systems).

###  II.2 Guiding ideas of the finite temperature SHDMC method and definitions

Instead of performing the usual projection for infinite imaginary time of a
single trial wave-function, we run DMC for multiple guiding wave-functions
(forming a linearly independent basis) for finite imaginary time, which is
equivalent to finite temperature. Instead of using a single real guiding
function with nodes, we use a set of complex antisymmetric guiding functions
without nodes. Therefore, the Hamiltonian $\mathcal{\hat{H}}$ is not altered
at the nodes as in the standard importance sampling DMC approachceperley80
with the fixed-node approximation.anderson79  As explained in the
introduction, extending DMC to finite temperatures requires to go beyond those
fixed schemes. Complex-valued antisymmetric wave functions, that do not have
nodal pockets, can be constructed as a linear combination of two real wave
function with different nodes (see Fig. 1). We go beyond the standard fixed-
phase approximation and the local-time approximation. stonehenge  As in the
SHDMC method for complex wave functions, stonehenge  the infamous sign problem
is avoided with complex antisymmetric guiding functions. The result is acurate
as long as enough statistical information is collected.

SHDMCkeystone ; rockandroll ; stonehenge  systematically improves a trial wave
function by maximizing the overlap with the wave function propagated in
imaginary time in DMC. Here instead of maximizing the overlap of a single wave
function we will maximize the overlap with the basis.

Following the ideas of the SHDMC method, we use a recursive approach. In every
iteration $\ell$, importance sampling DMC ceperley80  is performed and
statistical data of the evolution in $\beta$ of a set of $M_{S}$ guiding wave
functions is projected on the many-body bases
$\\{|\chi^{S,\ell}_{n}\rangle\\}$ and $\\{|n\rangle\\}$ . The statistical data
is used to improve the small basis $\\{|\chi^{S,\ell+1}_{n}\rangle\\}$ and the
guiding functions for the next iteration In what follows, we will omit the
iteration index $\ell$ in the notation for clarity, when the basis is not
changed. The small basis $\\{|\chi^{S}_{n}\rangle\\}$ is orthonormal:
$\langle\chi^{S}_{n}|\chi^{S}_{n}\rangle=\delta_{n,m}$.

For numerical efficiency, depending on the problem, we choose guiding
functions related to the eigenstates of
$\mathcal{\hat{U}}=e^{\hat{J}}e^{\mathcal{\hat{H}}}e^{-\hat{J}}$. (i)
$\\{|\chi^{U}_{n}\rangle\\}$ is formed by the Slater expansions of the
eigenstates in the small basis. (ii) $\\{|\chi^{V}_{n}\rangle\\}$ is formed by
the Slater expansion of linear combinations of eigenstates pairs.

We construct wave functions of the form

$\displaystyle\Psi^{T}_{n}({\bf R})=$ $\displaystyle e^{-J({\bf
R})}\chi^{X}_{n}({\bf R});$ (4)

where the super index $X$ refers to either $S$, $U$, or $V$, depending on the
case. To simplify the notation, we omit $X$ in $\Psi^{T}_{n}({\bf R})$.

We assume the Jastrow factor operator $e^{-\hat{J}}$ to be diagonal in the
many-body configuration space ${\bf R}$, and positive, which implies that it
must have an inverse. fn:back_flow  The Jastrow factor is fixed in SHDMC but
can be optimized variationaly so that the free energy of the system is
minimized.

In contrast with the CFDMC ceperley88  and the released phase jones97
methods, we use anti-symmetric guiding functions, which are improved
recursively with a maximum overlap criterion. Since the exponential growth of
the bosonic ground state is prevented by the guiding functions, the free
energy obtained is an upper bound. This approach is different to the
correlated linear method filippi09  because the wave function is optimized at
the DMC level and we use anti-symmetic guiding functions. In variance with the
original SHDMC approach for excited states rockandroll  multiple wave
functions are propagated in parallel. A serial orthogonalization step in the
original SHDMC method for excited states rockandroll ; stonehenge  is
postponed in this new approach until DMC has been run for the all basis
functions. fn:complications

###  II.3 Working with complex guiding wave functions to avoid the fixed-node
approximation

Figure 1:  (Color online) A schematic representation of the nodes of the real
and imaginary parts of a complex antisymmetric function. This type of complex
functions must be zero only at the coincidental points. Though real
antisymmetric functions must have nodal pockets a complex linear combination
must not. Note also that two different antisymmetric real wave functions
approaching different eigenstates will have different nodes and produce
different fixed-node Hamiltonians.

While complex guiding wave functions allow us to avoid the fixed-node
approximation, they introduce additional complications stonehenge  that are
discussed here. Once these complications are dealt with, the sign problem is
avoided as in the fixed-node. As in standard SHDMC the result is accurate as
long as enough statistics is obtained.

Following Refs. ortiz93,  and stonehenge, , $\Psi^{T}_{n}({\bf R})=\langle{\bf
R}|\Psi_{n}^{T}\rangle$ can be written as an explicit product of a complex
phase and an amplitude $\Psi^{T}_{n}({\bf R})=\Phi_{n}^{T}({\bf R})e^{{\bf
i}\phi_{n}({\bf R})}$.

The expressions

$\displaystyle\Phi_{n}^{T}({\bf R})$ $\displaystyle=\sqrt{\Psi_{n}^{T}({\bf
R})\Psi^{T*}_{n}({\bf R})},\text{ and }$ (5) $\displaystyle\phi_{n}({\bf R})$
$\displaystyle=\ln[\Psi_{n}^{T}({\bf R})/\Psi_{n}^{T*}({\bf R})]/(2i)+\pi m$
(6)

allow the computation of all the gradients and Laplacians in terms of those of
an arbitrary complex function $\chi_{n}^{X}({\bf R})$ and $e^{-J({\bf R})}$.

In Eq. (6) $m$ is an arbitrary integer that changes the Riemann branch of the
natural logarithm $\ln$ of a complex number. $m$ only contributes to the
gradient or Laplacian at the Reimann cuts. fn:reimann  Since the position of
the Reimann cuts is an arbitrary mathematical convention, their contribution
to gradients and Laplacians is unphysical and ignored.

The dependence in $\beta$ of $e^{-\beta\mathcal{\hat{H}}}\Psi_{n}^{T}({\bf
R})$ is given by

$\displaystyle\Psi_{n}^{T}({\bf R},\beta)=$
$\displaystyle\;e^{-\beta\hat{\mathcal{H}}}\Psi_{n}^{T}({\bf R})$ (7)
$\displaystyle=$
$\displaystyle\;e^{-\beta\hat{\mathcal{H}}}\left[\Phi_{n}^{T}({\bf R})e^{{\bf
i}\phi_{n}({\bf R})}\right]$ (8) $\displaystyle=$
$\displaystyle\;\Phi_{n}^{T}({\bf R},\beta)e^{{\bf i}\phi_{n}({\bf R})}.$ (9)

Equation (9) includes by definition all the temperature dependence in
$\Phi_{n}^{T}({\bf R},\beta)~{}=~{}\Psi_{n}^{T}({\bf R},\beta)/e^{{\bf
i}\phi_{n}({\bf R})}$, denoted as free-amplitude, since it can be complex,
stonehenge  with $\Phi_{n}^{T}({\bf R},0)=\Phi_{n}^{T}({\bf R})$. The phase
$\phi_{n}({\bf R})$ remains fixed in the interval $[0,\beta)$ as in Ref.
ortiz93,

Following Ref. stonehenge,  we define the quantity

$\displaystyle f_{n}({\bf R},\beta)=\Psi_{n}^{T*}({\bf R})\Psi^{T}_{n}({\bf
R,\beta})e^{\beta E^{X}_{n}},$ (10)

where $E^{X}_{n}$ is a reference energy adjusted numerically to satisfy the
condition
$\langle\chi^{X}_{n}|e^{\hat{J}}|\Psi^{T}_{n}({\beta})\rangle~{}=~{}1$. This
reference energy is different from the one commonly used to obtain the ground
state. In practice, $E^{X}_{n}$ depends on the Slater determinant expansion
$\chi^{X}_{n}({\bf R})$ used to construct the guiding wave function and
contains the relevant information required to calculate thermodynamical
averages.

Using Eqs (2) and (7), one can easily obtain

$\displaystyle\frac{\partial}{\partial\beta}f_{n}({\bf R},\beta)=$
$\displaystyle\Psi_{n}^{T*}({\bf
R})\frac{\partial}{\partial\beta}\left[\Psi^{T}_{n}({\bf R,\beta})e^{\beta
E^{X}_{n}}\right]$ $\displaystyle=$ $\displaystyle-\Psi_{n}^{T*}({\bf
R})\left[(\hat{\mathcal{H}}-E^{X}_{n})\Psi^{T}_{n}({\bf R,\beta})e^{\beta
E^{X}_{n}}\right]$ $\displaystyle=$ $\displaystyle-\left[E_{n}^{T}({\bf
R},\beta)-E^{X}_{n}\right]f_{n}({\bf R},\beta)$ (11)

with

$\displaystyle E_{n}^{T}({\bf R},\beta)$
$\displaystyle=\frac{\mathcal{\hat{H}}\Psi^{T}_{n}({\bf
R},\beta)}{\Psi^{T}_{n}({\bf R},\beta)}$
$\displaystyle=-\frac{1}{2}\\!\\!\sum_{j}^{N_{e}}\frac{\nabla_{j}^{2}\Phi_{n}^{T}({\bf
R},\beta)}{\Phi_{n}^{T}({\bf R},\beta)}$ (12)
$\displaystyle+\frac{1}{2}\sum_{j}^{N_{e}}\left|{\bf
A}_{j}+{\bf\nabla}_{j}\phi_{n}({\bf R})\right|^{2}+V({\bf R},\beta)$
$\displaystyle-{\bf
i}\sum_{j}^{N_{e}}\left\\{\frac{{\bf\nabla}_{j}\Phi_{n}^{T}({\bf
R},\beta)}{\Phi_{n}^{T}({\bf R},\beta)}\cdot\left[{\bf
A}_{j}\\!+\\!{\bf\nabla}_{j}\phi_{n}({\bf R})\right]\right.$
$\displaystyle\;\;\;\;\;\;\;\;\;\;\;\;\;+\left.\frac{{\bf\nabla}_{j}\cdot\left[{\bf
A}_{j}+{\bf\nabla_{j}}\phi_{n}({\bf R})\right]}{2}\right\\}.$

In order to perform an importance sampling using $\Phi_{n}^{T}({\bf R})$ as a
guiding wave function (as in Refs. ceperley80,  and ortiz93, ), it is
convenient to express the kinetic energy in terms of $\Phi_{n}^{T}({\bf R})$.
The term including ${\nabla_{j}^{2}\Phi_{n}^{T}({\bf
R},\beta)}/{\Phi_{n}^{T}({\bf R},\beta)}$ can be rewritten using the
identityceperley80 ; HLRbook ; mfoulkesrmp2001

$\displaystyle\frac{\nabla^{2}_{j}\Phi^{T}_{n}({\bf
R},\beta)}{\Phi_{n}^{T}({\bf R},\beta)}=$
$\displaystyle\frac{\nabla^{2}_{j}f_{n}({\bf R,\beta})}{f_{n}({\bf
R,\beta})}+\frac{\nabla^{2}_{j}\Phi^{T}_{n}({\bf R})}{\Phi^{T}_{n}({\bf R})}$
$\displaystyle-\frac{{\bf\nabla}_{j}\cdot\left[f_{n}({\bf R},\beta){\bf
F}_{Q}^{j}\right]}{f_{n}({\bf R,\beta})}\;,$ (13)

where

$\displaystyle{\bf F}_{Q}^{j}={\bf\nabla}_{j}ln\left|\Phi_{n}^{T}({\bf
R})\right|^{2}\;.$ (14)

Note that Eqs. (II.3) and (14) are valid as long as $\Phi^{T}_{n}({\bf R})\neq
0$ and $f_{n}({\bf R},\beta)\neq 0$. In practice, any divergence of ${\bf
F}_{Q}^{j}$ at the nodes enforces $f_{n}({{\bf R},\beta})$ to be zero. Figure
1 shows that complex antisymmetric wave functions can be constructed so that
they have nodes only at points ${\bf R}$ with ${\bf r}_{i}={\bf r}_{j}$. In
this case the nodal error is avoided but errors in the phase introduce a phase
shiftortiz93 ; stonehenge  and a complex contribution to $E^{L}_{n}({\bf
R},\beta)$. However, for complex wave functions without zeros, Eq. (II.3) is
always valid, except at the coincidental points (if cusp conditions are not
satisfied). To satisfy Eq. (II.3) at the coincidental points, a Jastrow factor
is introduced in Eq. (4). While using complex wave functions involves some
complications, the advantage is that the evolution in imaginary time $\beta$
describes the thermodynamical properties with $\beta^{-1}=k_{B}T$. However,
going beyond the fixed-phase approximation ortiz93 ; meyer13  is required to
obtain the thermodynamics. In this work the phase is not “released” in the
same sense of Ref. jones97 , it is only free within the small subspace.

Replacing Eq. (II.3) into Eq (II.3) and then into Eq. (II.3) one obtains

$\displaystyle\frac{\partial f({\bf R},\beta)}{\partial\beta}=$
$\displaystyle\frac{1}{2}\sum_{j}^{N_{e}}\left\\{\nabla_{j}^{2}f_{n}({\bf
R},\beta)-{\bf\nabla}_{j}\cdot\left[f_{n}({\bf R},\beta){\bf
F}_{Q}^{j}\right]\right\\}$ $\displaystyle-\left[E_{n}^{L}({\bf
R},\beta)-E^{X}_{n}\right]f_{n}({\bf R},\beta),$ (15)

with $E_{n}^{L}({\bf R},\beta)$ where

$\displaystyle E_{n}^{L}({\bf R},\beta)=$
$\displaystyle-\frac{1}{2}\\!\\!\sum_{j}^{N_{e}}\frac{\nabla_{j}^{2}\Phi_{n}^{T}({\bf
R})}{\Phi_{n}^{T}({\bf R})}$ (16)
$\displaystyle+\frac{1}{2}\sum_{j}^{N_{e}}\left|{\bf
A}_{j}+{\bf\nabla}_{j}\phi_{n}({\bf R})\right|^{2}+V({\bf R},\beta)$
$\displaystyle-{\bf
i}\sum_{j}^{N_{e}}\left\\{\frac{{\bf\nabla}_{j}\Phi_{n}^{T}({\bf
R},\beta)}{\Phi_{n}^{T}({\bf R},\beta)}\cdot\left[{\bf
A}_{j}\\!+\\!{\bf\nabla}_{j}\phi_{n}({\bf R})\right]\right.$
$\displaystyle\;\;\;\;\;\;\;\;\;\;\;\;\;+\left.\frac{{\bf\nabla}_{j}\cdot\left[{\bf
A}_{j}+{\bf\nabla_{j}}\phi_{n}({\bf R})\right]}{2}\right\\}.$

That is, the local energy now depends on two guiding functions (i)
$\Phi_{n}^{T}({\bf R})$, and (ii) $\Phi_{n}^{T}({\bf R},\beta)$ which is an
approximation that must be obtained and improved for $\beta\neq 0$.

The use of complex valued guiding functions originates the gradients
${\bf\nabla}_{j}\phi_{n}({\bf R})$ that appear in the local energy in Eq.
(16). Their contribution prevents the result from reaching the bosonic
solution and enforces an upper bound on the fermionic ground state. ortiz93
In addition, the contribution of ${\bf\nabla}_{j}\Phi_{n}^{T}({\bf R},\beta)$
must be taken into account in the presence of a magnetic field (see term
between the $\\{\\}$, when ${\bf A}_{j}\neq 0$), even when using a real-
bosonic guiding wave function with ${\bf\nabla}_{j}\phi_{n}({\bf R})=0$.

A locality approximation locality  has been used in the past when a non-local
pseudo potential is included in $\mathcal{\hat{H}}$ in the potential term
$\hat{V}({\bf R})$. It consists in replacing $V({\bf
R},\beta)=\frac{\hat{V}({\bf R})\Psi^{T}_{n}({\bf R},\beta)}{\Psi^{T}_{n}({\bf
R},\beta)}$ by $V({\bf R},0)$ in Eq.(II.3). Since $V({\bf R},\beta)\rightarrow
V({\bf R},0)$ for $\beta\rightarrow 0$, we will use the locality approximation
in the first iteration. However, we will improve it in subsequent iterations
(see Appendix A).

A local-time approximation, analogous to the locality approximation, locality
, was introduced in Ref. stonehenge,  to estimate the ratio

$\displaystyle\frac{{\bf\nabla}_{j}\Phi_{n}^{T}({\bf
R},\beta)}{\Phi_{n}^{T}({\bf R},\beta)}=$
$\displaystyle\frac{{\bf\nabla}_{j}\Phi_{n}^{T}({\bf R})}{\Phi_{n}^{T}({\bf
R})}+{\bf\nabla}_{j}\ln\left[\frac{\Phi_{n}^{T}({\bf
R},\beta)}{\Phi_{n}^{T}({\bf R})}\right].$ (17)

One can neglect the last term in Eq. (17) for $\beta\rightarrow 0$ where
${\Phi_{n}^{T}({\bf R},\beta)}/{\Phi_{n}^{T}({\bf R})}\rightarrow 1$. In the
present work we will use the local-time approximation only in the first
iteration. For subsequent iterations we improve the evaluation of Eq. (17) by
using the sampling of the dependence in $\beta$ of $\Phi_{n}^{T}({\bf
R},\beta)$ obtained in the previous iteration.

The locality and local-time approximations have little impact in optimization
methods that focus on eigenstates because the dependence on $\beta$ of
$\Psi^{T}_{n}({\bf R,\beta})$ is minimized when the optimization progresses as
$\Psi^{T}_{n}({\bf R})\rightarrow\Psi_{n}({\bf R})$. Going beyond these
approximations is required, however, to circumvent the nodes of the
eigenstates with complex wave functions. Fortunately, optimization in the
small subspace allows an easy sampling of the $\beta$ dependence. The approach
is exact if the big basis is large enough and if enough statistical data is
collected fn:fixedphase  as $\ell\rightarrow\infty$.

Circumventing the nodes with complex wave functions is necessary in this case
because, in standard DMC calculations using real-valued wave functions with
nodes, if any walker crosses a node, it is either killedceperley80  or the
move is rejected. umrigar93  This introduces an artificial divergent potential
at the nodal surface, which adds a kink at the node (a step for the rejection
case). Since there is a one-to-one correspondence between energy of one
eigenstate and its nodes, rosetta  eigenstates with different energies must
have different nodes. As a consequence, two real wave wave functions that
approach different eigenstates introduce different nodal potentials. Since the
fixed-node Hamiltonian is different for different eigenstates, and affect the
dynamics at the node in the evolution in imaginary time, the $\beta$
dependence obtained using the fixed-node approximation will not describe the
thermodynamics even if the exact nodes of the ground state are provided.

###  II.4 Differences with other DMC-like projection methods

The implementation of this method follows essentially the same approach
developed for DMC or SHDMC, with some key numerical changes.

Equation (15) is identical to Eq. (1) in Ref. ceperley80,  except for the
local energy, which now has an explicit dependence in $\beta$. Unlike Eq. (13)
in Ref. stonehenge, , Eq. (15) is now valid for $\beta\geq 0$. As in Ref
ceperley80, , Eq (15) describes the evolution of an ensemble. Each member of
the ensemble of configurations ${\bf R}_{i}$ (walker) undergoes (i) a random
diffusion and (ii) drifting by the quantum force $\sum_{j}{\bf F}_{Q}^{j}$
(which depends only on $\Phi_{T}({\bf R})$ and not on the phase). Following
Ref. stonehenge, , (iii) each walker carries a complex phase factor. In a
nonbranching algorithm, the complex weight of the walkers is multiplied by
$\exp\\{-\left[E^{L}_{n}({\bf R},\beta)-E^{X}_{n}\right]\delta\beta\\}$ at
every diffusion step. In contrast with Ref. jones97, , the phase factor of the
walkers starts in 1 and evolves towards the difference between the guiding
phase and the phase of $\Psi({\bf R},\beta)$, while in the release phase
approach the initial phase of each walker depends on the initial positions of
the walkers ${\bf R}_{i}$ but remains constant in $\beta$.

If the $\chi_{n}^{X}({\bf R})$ are linear combinations of antisymmetric
functions, with arbitrary complex coefficients, the $\Phi_{n}^{T}({\bf R})$ in
Eq. (16) do not have nodal surfaces (see Fig. 1). Therefore, ${\bf F}_{q}$ is
not divergent, but at the coincidental points.

All walkers must add the same inverse temperature $\beta$ after $k$ steps. The
standard time-step correction to minimize time step errors at the nodes [Eq.
(33) in Ref. umrigar93, ] is not strictly necessary since the divergences in
${\bf F}_{q}$ are removed. If one introduces it, one must readjust
$\delta\beta$ during the time evolution so that all the walkers add up to the
same $\beta$. For the same reason, the standard accept/reject scheme that
enforces detailed balanceHLRbook  is modified: other moves are retried after
rejection until a move is accepted.

###  II.5 Calculation of the partition function

This section shows that the partition function can be obtained as
$Z(\beta)\simeq\sum_{n=0}^{M_{S}-1}e^{-\beta E^{S}_{n}}$ where the $E^{S}_{n}$
are reference energies instead of eigen energies.

Numerically, it is convenient to start the calculation with a distribution of
walkers proportional to $f_{n}({\bf R},0)=|\Phi^{T}_{n}({\bf R})|^{2}$. In the
importance sampling approach ceperley80  setting the second line of Eq. (15)
equal to zero, provides an equilibrium distribution proportional to
$|\Phi^{T}_{n}({\bf R})|^{2}$.

As in the DMC and SHDMC methods, the evolution in inverse temperature $\beta$
is discretized into $k$ finite steps $\delta\beta=\beta/k$. Following the
SHDMC approachkeystone ; rockandroll ; rollingstones ; stonehenge  the
weighted distribution of the walkers can be written as

$\displaystyle f_{n}({\bf R},k\delta\beta)$
$\displaystyle=\lim_{N_{c}\rightarrow\infty}\frac{1}{N_{c}}\sum_{i=1}^{N_{c}}W_{i}\delta\left({\bf
R-R}_{i}\right).$ (18)

In Eq. (18), ${\bf R}_{i}$ corresponds to the position of the walker $i$ , and
$N_{c}$ is the number of equilibrated configurations. The complex weights
$W_{i}$ are given by

$\displaystyle W_{i}=e^{-\left[E_{i}^{k}-E^{X}_{n}\right]k\delta\beta}$ (19)

with

$\displaystyle E_{i}^{k}=\frac{1}{k}\\!\sum_{j=0}^{k-1}E^{L}_{n}({\bf
R}_{i}^{-j}),$ (20)

Where $k$ is a number of steps and $E_{L}({\bf R}_{i}^{-j})$ is the previous
value of the local energy obtained $j$ steps $\delta\beta$ earlier for the
walker $i$.

The evolution in inverse temperature $\beta$ of the guiding wave function
$\Psi_{n}^{T}({\bf R},\beta)$ can be written, without loss of generality, as

$\displaystyle\Psi_{n}^{T}({\bf R},\beta)=e^{-\beta E^{X}_{n}}e^{-J({\bf
R})}\left[\chi^{X}_{n}({\bf R})+\delta\chi^{X}_{n}({\bf R},\beta)\right].$
(21)

That is, the product of an average decay factor $e^{-\beta E^{X}_{n}}$ times
the Slater determinant part. The Slater part is given by the the one at
$\beta=0$ plus an orthogonal displacement $\delta\chi^{X}_{n}({\bf R},\beta)$.
The $X$ in term $\delta\chi^{X}_{n}({\bf R},\beta)$ denotes the explicit
dependence on $\chi^{X}_{n}({\bf R})$. At least one overlap
$\langle\chi^{X}_{m}|\delta\chi^{X}_{n},\beta\rangle$ must be different from
zero for $n\neq m$ , if $\Psi^{T}_{n}({\bf R})$ is not an eigenstate.

Using Eq. (21) we can correct equation (16) beyond the locality and local-time
approximations. The displacement $\delta\chi^{X}_{n}({\bf R},\beta)$ can be
sampled from the DMC run as follows: From Eqs. (7), (10) and (18), one can
obtain

$\displaystyle e^{-J({\bf R})}\delta\chi^{X}_{n}({\bf R},\beta)=$
$\displaystyle\langle{\bf R}|\left[e^{\beta
E^{X}_{n}}|\Psi^{T}_{n}(\beta)\rangle-|\Psi^{T}_{n}\rangle\right]$
$\displaystyle=$ $\displaystyle\;\langle{\bf
R}|\left[e^{-\beta(\hat{\mathcal{H}}-E^{X}_{n})}-1\right]e^{-\hat{J}}|\chi^{X}_{n}\rangle\;$
(22) $\displaystyle=$ $\displaystyle\;\frac{\left[f_{n}({\bf
R},\beta)-f_{n}({\bf R},0)\right]}{\Psi^{T*}_{n}({\bf R})}$ (23)
$\displaystyle\simeq$ $\displaystyle\frac{1}{N_{c}}\sum_{i=1}^{N_{c}}e^{J({\bf
R})}\frac{[W_{i}-1]}{\chi^{X*}_{n}({\bf R})}\delta({\bf R}-{\bf R}_{i}).$ (24)

Within the subspace spanned by the basis $\\{\chi^{S}_{n}({\bf R})\\}$, the
identity operator $\hat{E}$ is given by

$\displaystyle\langle{\bf R}^{\prime}|\hat{E}|{\bf R}\rangle$
$\displaystyle=\sum_{m=0}^{M_{S}-1}e^{-J({\bf R^{\prime}})}\chi^{S}_{m}({\bf
R^{\prime}})\chi^{S*}_{m}({\bf R})e^{J({\bf R})}.$ (25)

Applying Eq. (25) to both sides of Eq. (24), and integrating over ${\bf R}$,
one can easily obtain an expression of the diffusion displacement within the
basis $\\{|\chi^{S}_{n}\rangle\\}$ as

$\displaystyle\delta\chi^{X}_{n}({\bf R},\beta)$
$\displaystyle=\sum_{m=0}^{M_{S}-1}\lambda_{n}^{m}(\beta)\chi^{S}_{m}({\bf
R})$ (26)

with

$\displaystyle\lambda_{n}^{m}(\beta)=\frac{1}{N_{c}}\sum_{i=1}^{N_{c}}e^{2J({\bf
R}_{i})}\frac{\chi^{S*}_{m}({\bf R}_{i})}{\chi^{X*}_{n}({\bf
R}_{i})}[W_{i}-1]$

where $N_{c}=\sum_{i=1}^{N_{c}}W_{i}$.

Using $X=S$ in Eqs. (18) -(21) one can prove that

$\displaystyle\langle\chi^{S}_{n}|\mathcal{\hat{U}}|\chi^{S}_{m}\rangle=$
$\displaystyle\int{\bf dR}e^{J({\bf R})}\chi^{S*}_{n}({\bf
R})\Psi_{m}^{T}({\bf R},\beta)$ $\displaystyle=$ $\displaystyle e^{-\beta
E^{S}_{n}}\left[\delta_{n,m}+\int{\bf dR}\chi^{S*}_{n}({\bf
R})\delta\chi^{S}_{m}({\bf R},\beta)\right]$ $\displaystyle=$ $\displaystyle
e^{-\beta E^{S}_{n}}(\delta_{n,m}+\lambda_{n}^{m})$ (28)

with $\mathcal{\hat{U}}$ having the structure of the transcorrelated method
ten-no02

$\displaystyle\mathcal{\hat{U}}=$ $\displaystyle
e^{\hat{J}}e^{-\beta\mathcal{\hat{H}}}e^{-\hat{J}}.$ (29)

We use condition $\lambda_{n}^{n}=0$ [See Eq (LABEL:eq:deltalambda)] to
determine the value of $E^{S}_{n}$. In practice, we adjust the reference
energy of the guiding functions every iteration as $e^{-\beta
E^{X,\ell+1}_{n}}=e^{-\beta E^{X,\ell}_{n}}(\delta_{n,m}+\lambda_{n}^{m})$.

Since $tr(e^{-\beta\mathcal{\hat{H}}})=tr(\mathcal{\hat{U}})$, the
contribution to the Helmholtz free energy of the small subspace is given by

$\displaystyle
F(\beta)=-\frac{1}{\beta}\ln{\left[\sum_{n=0}^{M_{S}-1}e^{-\beta
E^{S}_{n}}\right]},$ (30)

where the expression inside the brackets is the partition function $Z(\beta)$.

In general for an arbitrary guiding function $\Psi_{m}^{T}({\bf R})$, the
variance will grow with $M_{S}$. An energy span larger than $\beta^{-1}$ must
be retained in the basis to calculate thermodynamical properties. Arbitrary
trial wave functions spanned by this space might have significant variance in
the walkers weights. To reduce the variance we use guiding functions that are
approximately linear combinations of a pair of neighboring eigenstates.

When using guiding functions that are different from the small basis
functions, the contributions to the trace of the density matrix in the small
basis can be obtained with

$\displaystyle e^{-\beta
E^{S}_{n}}=\sum_{i}\left[|\langle\chi^{X}_{i}|\chi^{S}_{n}\rangle|^{2}+\langle\chi^{X}_{i}|\chi^{S}_{n}\rangle\langle\chi^{S}_{n}|\delta\chi^{X}_{i}\rangle\right]e^{-\beta
E^{X}_{i}}\;.$ (31)

The details of the derivation are in Appendix B.

In a recent paper, Mazzola, Zen and Sorellamazzola12  proved that

$\displaystyle\langle\Psi^{T}_{n}|e^{-\beta\hat{\mathcal{H}}}|\Psi^{T}_{n}\rangle\geq
e^{-\beta\langle\Psi^{T}_{n}|\hat{\mathcal{H}}|\Psi^{T}_{n},\rangle}.$ (32)

Ref. mazzola12,  used the righthand side of Eq. (32) to approximate the free
energy obtaining a lower bound for $F(\beta)$. Reference mazzola12,  can be
considered a VMC approach to the evaluation of the free energy. That
approximation becomes exact if all the $|\Psi^{T}_{n}\rangle$ are eigenstates
of $\hat{\mathcal{H}}$. However, that method is very poor for an arbitrary
random guiding function. In the present approach, we go beyond Ref. mazzola12,
by evaluating the lefthand side of Eq. (32) directly using DMC.

In many situations, the excitations of a mean field method based on
approximations DFT might be good enough to obtain the low energy
thermodynamical properties using Eq. (30). If that were the case, at least two
DMC runs for each function of the basis are required. One to obtain the
$\beta$ dependence and a second to evaluate the reference energies beyond the
local-time approximation. However, in the so-called highly correlated
materials, usual approximations of DFT fail to describe the low energy
physics. In those cases a method that could optimize the basis is more
important. That method is described in the following subsections.

###  II.6 The first iteration: Construction of the first small basis
$\\{|\chi^{S,1}_{n}\rangle\\}$

While the present approach will optimize the basis from any starting basis
set, the calculation will be more efficient starting from a good basis. A
procedure to generate a good starting set is described here.

The only restriction for the small basis $\\{|\chi^{S}_{n}\rangle\\}$ is to
avoid the nodes associated with real wave functions. In this work we choose
the initial basis with a Lanczos-like procedure combined with the SHDMC
approach.

The big subspace basis set ${|m\rangle}$ is constructed by symmetry
constrained functions (linear combinations of Slater determinants with the
same symmetry of the ground state) ordered with increasing mean field energy.

We choose the first basis function of the small subspace to be

$\displaystyle|\Psi_{0}^{S}\rangle=\frac{1}{\sqrt{2}}(|0\rangle+{\bf
i}|1\rangle)$ (33)

being $|0\rangle$ and $|1\rangle$the ground and first excited states of a non-
interacting solution of the system.

Using Eq. (33) as guiding function in Eqs. (II.5)-(24) and replacing
$\chi^{S}_{m}({\bf R})$ by $\langle{\bf R}|m\rangle$ in Eq. (25), but not on
$\Psi^{T}_{n}({\bf R})$, we obtain an expression similar to Eq. (26)

$\displaystyle\delta\tilde{\chi}_{n}^{S}({\bf R})=\langle{\bf
R}|\delta\tilde{\chi}^{S}_{n}\rangle=\tilde{\sum}_{m=0}c_{n}^{m}\langle{\bf
R}|m\rangle.$ (34)

The tilde in $\delta\tilde{\chi}_{n}^{S}({\bf R})$ means that the expansion is
in the big basis $\\{|m\rangle\\}$ with

$\displaystyle c_{n}^{m}(\beta)=\frac{1}{N_{c}}\sum_{i=1}^{N_{c}}e^{2J({\bf
R}_{i})}\frac{\langle m|{\bf R}_{i}\rangle}{\chi^{S*}_{n}({\bf
R}_{i})}[W_{i}-1].$

The symbol $\tilde{\sum}$ in Eq. (34) means that the sum is restricted to the
coefficients $c_{n}^{m}$ with an error bar smaller than $25\%$ of the absolute
value (this is the standard recipe of the SHDMC algorithm keystone ).

We define the next basis function $|\chi^{S}_{n+1}\rangle$ recursively as

$\displaystyle|\chi^{S}_{n+1}\rangle=$
$\displaystyle\frac{1}{\mathcal{N}_{n+1}}\hat{P}_{n+1}|\delta\tilde{\chi}^{S}_{n}\rangle\text{
with }$ (36) $\displaystyle\hat{P}_{n}=$ $\displaystyle
1-\sum_{m=0}^{n-1}|\chi^{S}_{m}\rangle\langle\chi^{S}_{m}|$ (37)

where $\mathcal{N}_{n+1}$ is a normalization constant. Equations (37) and (36)
mean that $|\chi^{S}_{n}\rangle$ is the projection of the displacement
$|\delta\tilde{\chi}^{S}_{n+1}\rangle$ orthogonal to the subspace spanned by
the $n$ basis functions found previously.

One repeats this procedure until a basis of $M_{S}$ functions
$\\{|\chi^{S}_{n}\rangle\\}$ is constructed. This Lanczos-like procedure
grants that the initial small subspace basis times the Jastrow factor has a
large projection onto the lowest-energy eigenstates of $\mathcal{\hat{H}}$ or
the largest eigenstates of $e^{-\beta\mathcal{\hat{H}}}$.

Since the evolution in inverse temperature $\beta$ is not known during the
initialization step, we use the local-time approximation discussed in the
previous section. However, once a basis is generated, we can go beyond the
local-time approximation in successive iterations. Note that by construction
any $\chi^{S}_{n}({\bf R},\beta)$ can be approximated as a linear combination
of the basis function $\chi^{S}_{m}({\bf R})$ with $m<n+1$, since the finite
temperature projection of one wave function of the basis into the other was
used to construct the small basis. The details on how to approximate the
evolution in $\beta$ are in Appendix A.

###  II.7 Systematic improvement of the small basis
$\\{|\chi^{S}_{n}\rangle\\}$

One of the goals of this work is to obtain a much smaller basis of $M_{S}$
functions $\\{|\chi^{S}_{n}\rangle\\}$ than $\\{|m\rangle\\}$, the big set of
$N_{B}$ basis functions. The small basis $\\{|\chi^{S}_{n}\rangle\\}$ should
retain the lowest energy physics of $\mathcal{\hat{H}}$ and
$e^{-\beta\mathcal{\hat{H}}}$. While for some purposes (e.g. the ground state
calculation), the initial basis set described in the previous section might be
enough, in this subsection we describe how to further optimize the small basis
so it better describes the low-energy excitations of $\mathcal{\hat{H}}$ and
thermodynamical properties.

Note that $\chi^{S}_{n}({\bf
R},\beta)=e^{\hat{J}}e^{-\beta\mathcal{\hat{H}}}e^{-\hat{J}}\chi^{S}_{n}({\bf
R})$ will converge to the antisymmetric part of the ground state wave-function
as $\beta\rightarrow\infty$. In order to avoid every state in the basis
collapsing to the same function we (i) remove the projection into the other
states of the basis, (ii) add the diffusion displacement orthogonal to the
small subspace, and (iii) perform a Gram Schmidt orthogonalization as follows:

$\displaystyle|\chi^{S,\ell+1}_{0}\rangle=$
$\displaystyle|\chi^{S,\ell}_{0}\rangle+|\delta\tilde{\chi}_{0}^{S}\rangle-|\delta\chi^{S}_{0}(\beta)\rangle$
(38) $\displaystyle|\chi^{S,\ell+1}_{n}\rangle=$
$\displaystyle\hat{P}^{\ell+1}_{n}\left[|\chi^{S,\ell}_{n}\rangle+|\delta\tilde{\chi}_{n}^{S}\rangle-|\delta\chi^{S}_{1}(\beta)\rangle\right]\;\;n>0,$

with $\hat{P}^{\ell+1}_{n}$ given by Eq. (37) replacing $\chi^{S}_{m}$ by
$\chi^{S,\ell+1}_{m}$. Note in Eq. (38) that
$|\delta\tilde{\chi}_{n}^{S}\rangle$ given by Eq. (34) is a direct projection
of the diffusion displacement into the big basis of SCFs $\\{|m\rangle\\}$,
whereas $|\delta\chi^{S}_{n}(\beta)\rangle$ given by Eq. (26) is an indirect
projection (since the $|\delta\chi^{S}_{n}(\beta)\rangle$ are projected into
the small basis $\\{|\chi^{S}_{n}\rangle\\}$ which in turn are linear
expansions of functions that belong to $\\{|m\rangle\\}$). The difference
$|\Delta\tilde{\chi}^{S}_{n}\rangle=|\delta\tilde{\chi}^{S}_{n}\rangle-|\delta\chi^{S}_{n}(\beta)\rangle$
is by construction orthogonal to the small subspace. Accordingly, it describes
the decay of the small subspace basis into the eigenstates with the lowests
energies.

##  III Algorithm

The goal of this algorithm is to optimize a minimal basis
$\\{|\chi^{S}_{n}\rangle\\}$ to span the lowest energy excitations of
$\mathcal{\hat{H}}$ (equivalently the eigenstates of
$e^{-\beta\mathcal{\hat{H}}}$ with the largest eigenvalues). That basis can be
used to calculate finite temperature expectation values of thermodynamical
properties and accelerate the calculation of the ground and lower excited
states. In this section we summarize how the theory described in detail
earlier can be implemented.

Initialization: the small basis $\\{|\chi^{S}_{n}\rangle\\}$, a set of
orthogonal linear combinations of many-body functions $|m\rangle$ is
constructed using the procedure described in II.6. Once this procedure is
concluded, a set of linearly independent guiding functions
$e^{-\hat{J}}|\chi^{V}_{n}\rangle$ is constructed as linear combinations of
pairs of approximated eigenstates of $\mathcal{\hat{U}}$.

Subsequently, we can use for the evaluation of the local energy
$E^{L}_{n}({\bf R},\beta)$ an approximate dependence of the guiding functions
in $\beta$ given by Eq. (47).

Basis update iteration: each iteration $\ell$ is composed by (1) a parallel
diffusion of intermediate functions $|\chi^{V}_{n}\rangle$. (2) A linear
transformation to obtain $|\delta\tilde{\chi}^{S}_{n}\rangle$ and
$|\delta\chi^{S}_{n}\rangle$. (3) Recalculation of $\mathcal{\hat{U}}$ in the
small basis. (4) Update of the small basis for the next iteration
$|\tilde{\chi}^{S,\ell+1}_{n}\rangle$. (5) Update of the intermediate
functions $|\chi^{V,\ell+1}_{n}\rangle$. (6) Finally, the algorithm decides to
increase the number of samples in the next iteration or not. These six steps
are repeated recursively.

The individual steps of the iteration are described below in more detail:

(1) Parallel diffusion: Each displacement

$\displaystyle|\delta\chi^{V}_{n}\rangle=e^{\hat{J}}\left[e^{-\beta\mathcal{\hat{H}}}-1\right]e^{-\hat{J}}|\chi^{V}_{n}\rangle$
(39)

is projected into the small basis $\\{|\chi^{S}_{n}\rangle\\}$ and the big
basis $\\{|m\rangle\\}$ using Eqs. (LABEL:eq:deltalambda) and Eq.
(LABEL:eq:cnm) respectively replacing $\chi^{X}_{n}({\bf R})$ by
$\chi^{V}_{n}({\bf R})$.

Each diffusion contains $M_{b}$ sampling subblocks. For each sampling
subblock, uncorrelated walker positions are generated from the previous one
with a VMC algorithm. Next, DMC is run for $k$ steps with a shorter time step
$\delta\beta$. The coefficients of $|\delta\tilde{\chi}^{V}_{n}\rangle$ and
$|\delta\chi^{V}_{n}\rangle$ are sampled at the end of each subblock using
Eqs. (LABEL:eq:deltalambda) and (LABEL:eq:cnm). Statistical data is collected
for $M_{b}$ subblocks for each parallel diffusion before an update of the
small basis.

(2) Linear transformation: The $|\delta\tilde{\chi}^{S}_{n}\rangle$ and
$|\delta\chi^{S}_{n}\rangle$ can be constructed in terms of
$|\delta\tilde{\chi}^{V}_{n}\rangle$ and $|\delta\chi^{V}_{n}\rangle$ using
Eq. (B) replacing the super index $X$ by $V$.

(3) Calculation of $\mathcal{\hat{U}}$, $\\{|\Psi^{U,\ell+1}_{n}\rangle\\}$,
and $\\{|\Psi^{V,\ell+1}_{n}\rangle\\}$ : A matrix representation of
$\mathcal{\hat{U}}$ is obtained in the small subspace
$\\{|\Psi^{U,\ell}_{n}\rangle\\}$, using Eqs. (LABEL:eq:deltalambda), (II.5)
and (B). The left and right eigenvectors are obtained by diagonalizing
$\mathcal{\hat{U}}$.

(4) Update of the small Basis: Equation (38) is used to improve the small
basis.

(5) We perform the correspondence
$\chi^{S,\ell}_{n}\rightarrow\chi^{S,\ell+1}_{n}$ (see Appendix A ) and
construct the basis $\\{|\chi^{U,\ell+1}_{n}\rangle\\}$, and
$\\{|\chi^{V,\ell+1}_{n}\rangle\\}$ with the coefficients of the eigenvectors
of $\mathcal{\hat{U}}$ in the iteration $\ell$.

(6) Updating $M_{b}$: At first, the number of sampling subblocks $M_{b}$ is
set to a small number (e.g. $M_{b}=3$). When the noise is dominant
$\sum_{n}\langle\Delta\tilde{\chi}^{S,\ell+1}_{n}|\Delta\tilde{\chi}^{S,\ell}_{n}\rangle\leq
0$, we increase $M_{b}$ by a factor larger than 1. As a result, the total
number of configurations $N_{c}$ sampled increases as the iteration $\ell$
increases and the statistical error is reduced. Hence, as the statistics are
improved the number of basis functions retained in the expansion Eq. (25)
increases over time.

(7) We use Eq. (30) to calculate thermodynamical properties.

##  IV Results in model calculations

This section describes the results obtained for a model system with an applied
magnetic field. The results are compared with configuration interaction (CI)
calculations in the same model used in Refs. rosetta, ; keystone, ;
rockandroll,  and stonehenge, . The lowest-energy eigenstates were found for
two polarized electrons ($J=1$) moving in a two-dimensional square with a side
length $1$ and a repulsive interaction potential of the form $\hat{V}({\bf
r},{\bf
r^{\prime}})=8\pi^{2}\gamma\cos{[\alpha\pi(x-x^{\prime})]}\cos{[\alpha\pi(y-y^{\prime})]}$
with $\alpha=1/\pi$ and $\gamma=4$.

The main advantage of the model is that fully converged CI calculations can be
performed which are nearly analytical. In order to perform the CI calculations
the many-body wave function of the small basis $\\{\chi^{S}_{n}\rangle\\}$ are
spanned in a basis of functions $\\{|m\rangle\\}$ that are eigenstates of the
noninteracting system. They are linear combinations of functions of the form
$\prod_{\nu}\sin(m_{\nu}\pi x_{\nu})$ with $m_{\nu}\leq 7$. Converged CI
calculations were performed to obtain a nearly exact expression of the lowest
energy states of the system $|\Psi_{n}\rangle=\sum_{m}a_{m}^{n}|m\rangle$. The
matrix elements involving the magnetic vector potential ${\bf A}$ (in the
symmetric gauge) were calculated analytically. The result of the CI
calculations were used to evaluate the partition functions and to quantify the
convergence of the basis.

The same basis used to construct the CI Hamiltonian is used as the big basis
to test our finite temperature version of SHDMC. All the calculations reported
are with $J({\bf R})=0$, which increases the statistically noise, makes the
test more difficult and facilitates the comparison with the CI results. The
results presented here are a proof of principle on the validity of the
algorithm, which is necessary before requesting and using the massive amount
of computing time required for realistic finite temperature calculations in
solids. While clearly a demonstration in a realistic system is required in the
future, a comparison with an exact model is the first essential step to
validate the scheme. This includes not only the value obtained for the
partition function but also a detailed analysis of the convergence of the
basis.

In the absence of magnetic fields there are two degenerate solutions: one that
transforms line $x$, and the other that transform like $y$. This degeneracy is
broken with a magnetic field. The eigenstates transform like $x+{\bf i}y$ and
$x-{\bf i}y$. Figure 2 shows the evolution of
$Z_{+}(\beta)=tr_{+}(\mathcal{\hat{U}})$ of the model system with $M_{S}=20$.
The subindex “$+$” in $Z_{+}$ and $tr_{+}$ means that the results of Fig. (2)
were obtained considering only the subspace of the Hamiltonian that transforms
like $x+{\bf i}y$. The calculation of thermodynamical properties requires,
however, the inclusion of all possible symmetries of the wave function, which
implies that $Z_{-}$, the trace in a small basis that transforms as $x-{\bf
i}y$, should also be added. In order to calculate $Z_{+}(\beta)$ we have
defined the zero of energy to be the ground state of the CI. The calculations
were run using $\delta\beta=0.00002$ and $\beta=0.004$. We have used a
magnetic field of $B=0.6283$.

Figure 2 shows the value of $Z_{+}(\beta)$, calculated with different methods,
relative to the full CI value obtained with the $M_{S}=20$ lowest eigenvalues.
The blue cycles were obtained with SHDMC using Eq. (30). The red rhombi were
obtained by evaluating
$Z_{+}(\beta)=\sum_{n,m}|\langle\Psi^{T}_{m}|\Psi_{n}\rangle|^{2}e^{-\beta
E_{n}}$, being the $|\Psi_{n}\rangle$ and $E_{n}$ the eigenvectors and
eigenvalues of the full CI. The empty squares were obtained as
$Z_{+}(\beta)\simeq\sum_{n}e^{-\beta\langle|\Psi^{T}_{n}|\hat{\mathcal{H}}|\Psi^{T}_{n}\rangle}$
being the $|\Psi^{T}_{n}\rangle$ linear combinations of pairs of approximated
eigenstates obtained with SHDMC. Therefore, the squares correspond to the
result that one would had obtained using the approximation of Ref. mazzola12,
in Eq. (32) for a very good set of functions. Finally, the up triangles mark
the result obtained with
$Z_{+}(\beta)\simeq\sum_{n}e^{-\beta\langle|\Psi^{T}_{n}|\hat{R^{-1}}\hat{\mathcal{H}}\hat{R}|\Psi^{T}_{n}\rangle}$,
being $\hat{R}$ a random rotation defined in the small subspace.

Figure 2 shows that as the iteration $\ell$ increases, the $Z_{+}(\beta)$
obtained with all methods increases. Note, that the scale of the $y$ axis does
not start from zero. The initialization scheme already produces a basis that
retains 90 % of the exact value of the truncated partition function. Similar
to previous SHDMC methods, the present generalization optimizes the small
basis overlap, not their average energy. The trace of $\mathcal{\hat{U}}$
increases indirectly as the small basis approaches to the subspace of the
eigenvectors with lowest energy. Note that the SHDMC results are within $\sim
1.5$% below the values obtained analytically by projection into the CI data
(red rhombi). This difference is due to the remaining errors in the evolution
phase which neglects the projection orthogonal to the small subspace. The
method used in Ref. mazzola12,  applied to approximated pairs of eigenstates
(empty squares) gives results only slightly below the SHDMC values, because
the energy difference between eigenstates is much smaller that $\beta^{-1}$.
However, a random rotation of the basis that spans the same subspace (up
triangles) would had produced a significantly worse result. This shows that
the approach described in Ref. 12 is significantly worse if each element of
the basis does not have a large projection on each eigenstate.

Note in the inset of Fig. 2 that the number of sampling blocks remains very
low ($M_{b}=3$) for the first 25 iterations and starts increasing when the
noise becomes dominant around $\ell=32$. If one considers only the first 32
iterations a significant improvement of the partition function is achieved
with little computational cost as compared with the calculation of individual
eigenstates.

Figure 2:  (Color online) Evolution of the trace of
$\mathcal{\hat{U}}=e^{\hat{J}}e^{-\beta\mathcal{\hat{H}}}e^{-\hat{J}}$ in the
subspace defined by the small basis $\\{|\chi^{S,\ell}_{n}\rangle\\}$ as a
function of the iteration index $\ell$ for $M_{S}=20$ relative to the
truncated trace obtained with full CI eigenvalues. Note that the scale starts
from 86% of the full CI value. The blue circles are results obtained with
SHDMC using Eq. (30). The red rhombi denote the exact evaluation of
$\mathcal{\hat{U}}=e^{\hat{J}}e^{-\beta\mathcal{\hat{H}}}e^{-\hat{J}}$ in the
small basis, using full CI data. The empty squares (up triangles) were
obtained by evaluating
$\langle\Psi^{T}_{n}|\hat{\mathcal{H}}|\Psi^{T}_{n},\rangle$ for linear
combination of eigenstate pairs (random rotation within the small subspace) in
the full CI basis while $Z(\beta)$ was obtained using an lower bound
approximation. mazzola12  The inset shows the computational cost, the number
of sampling subblocks $M_{b}$, as a function of the iteration.

Figure 3 quantifies the convergence of the small basis
$\\{|\chi^{S,\ell}_{n}\rangle\\}$ towards the basis defined by the eigenstates
of the CI $\\{|\Psi_{n}\rangle\\}$. For that purpose we define the logarithm
of the residual subspace projection as

$\displaystyle
LR_{sp}=\ln\left\\{1-|Det\left[\langle\Psi_{n}|\chi^{S,\ell}_{m}\rangle\right]|^{1/M_{S}}\right\\}.$
(40)

In Eq. (40) $Det\left[\langle\Psi_{n}|\chi^{S,\ell}_{m}\rangle\right]$ is the
Determinant of a square matrix of size $M_{S}$ formed by the overlap
$\langle\Psi_{n}|\chi^{S,\ell}_{m}\rangle$. The determinant of the matrix is a
complex number of modulus $1$ in the limit when any eigenstate
$|\Psi_{n<M_{S}}\rangle$ can be written as a linear combination of
$|\chi^{S,\ell}_{m}\rangle$. Any error in any member of the small basis
$\\{|\chi^{S,\ell}_{m}\rangle\\}$ reduces the modulus of the determinant by a
factor. The exponent $1/M_{S}$ in Eq. (40) is a standard geometric average. A
large negative value in Eq. (40) indicates a very good small basis with a
determinant that is approaching $1$.

Figure 3 shows the evolution of $LR_{sp}$ given by Eq. (40) as a function of
the iteration index $\ell$ for the same system described in Fig. 2. Note that
$LR_{sp}$ becomes increasingly negative as a function of $\ell$, which implies
a global improvement of the basis approaching to the one described by the
eigenstates of the full CI.

Figure 3:  Logarithm of the residual subspace projection $LR_{sp}$ [see Eq.
40] of the small basis $\\{|\chi^{S,\ell}_{n}\rangle\\}$ into the basis
$\\{|\Psi_{n}\rangle\\}$ defined by the eigenstates of the full CI as a
function of the iteration index $\ell$ for the same system described in Fig 2.
A larger negative number correspond to in improved small basis.

We next need to characterize how well an individual eigenstate can be
described by the small basis. To measure this we define the logarithm of the
residual projection $LR_{n}$ as

$\displaystyle
LR_{n}=\ln\left[1-\sum_{m=0}^{M_{S}-1}|\langle\Psi_{n}|\chi^{S,\ell}_{m}\rangle|^{2}\right].$
(41)

Note in Eq. (41) that, if the normalized eigenstate can be written as a linear
combination of the small basis $\\{|\chi^{S,\ell}_{m}\rangle\\}$, the
expression in the brackets should be zero. A large negative number in $LR_{n}$
implies that the eigenstate $|\Psi_{n}\rangle$ is very well described in the
small basis.

Figure 4 describes the evolution of $LR_{n}$ for different eigenstates of the
CI as a function of the iteration index $\ell$.

The blue (red) contribution to the color decreases (increases) as the index
$n$ increases. The continuous line follows the ground state. One can clearly
observe that the ground state of the CI is already very well described at the
initialization stage within the Lanczos-like setup. As the iteration $\ell$
increases, the small basis describes the lowest-energy excitations better
while higher excitations require more iterations. Note that for 25 iterations
15 eigenvectors are very well described within a basis of 20 states. The total
cost at this point is 300 000 individual DMC steps. The calculations of 15
eigenstates with the original SHDMC algorithm for excited statesrockandroll
would had cost at least twice as much. The current algorithm becomes
competitive, in addition, if one considers that it is parallel, which allows
to distribute this cost in multiple tasks ($M_{S}$) reducing the time to
solution to 2% as compared with the original SHDMC algorithm for excited
states.

Figure 4:  (Color online) Evolution of the logarithm of residual projection
$LR_{n}$ as a function of the iteration index $\ell$. Blue (red) color denotes
lower (higher) values of $n$. The continuous line follows the ground state
eigenstate.

Finally, for infinite statistics one could in principle obtain the
eigenenergies of $\mathcal{\hat{H}}$ from the eigenvalues of $u_{n}$ as
$E_{n}\approx-\ln(u_{n})/\beta$. This procedure is known to be inefficient to
obtain the eigenenergies which are better described by sampling
$\mathcal{\hat{H}}$ as in the CFDMC approach. The off diagonal noise in the
matrix elements of $\mathcal{\hat{U}}$ has a perverse effect on the magnitude
of the eigenvalues. Therefore, while this method is an efficient one to
optimize the basis, it should be combined with other methods to obtain the
eigenvalue spectra.

##  V Summary and discussions

In this paper we have presented a general framework aimed to calculate
thermodynamical properties of many-body system in an importance sampling DMC
context. ceperley80

We showed that a many-body basis can be optimized to describe a small subspace
maximizing the overlap with the subspace described by the lowest eigenstates
of the Hamiltonian. The Helmholtz free energy obtained within this truncated
basis is an upper bound of the exact free energy of a system. The
corresponding partition function is a lower bound of the exact partition
function.

This generalization of the SHDMC method for finite temperature takes advantage
of complex wave functions that do not have nodal pockets. Accordingly, we
avoid the appearance of the artificial potentials in the standard fixed-node
approximation when the amplitude of the importance sampling guiding function
is zero. The antisymmetric properties of the wave-function are enforced by a
complex phase. This introduces a complex contribution to the local energy. The
complex local energy is handled in the complex weight of the walkers. Going
beyond the local-time approximation used in Ref. stonehenge, , the evolution
of the complex phase factor for $\beta>0$ is now approximately taken into
account in the evaluation local energy. The evolution in $\beta$ becomes exact
as the statistical error is reduced as the sampling increases.

In variance CFDMC, ceperley88  the solution remains an upper bound of the
fermionic free energy. While the CFDMC approach uses a bosonic-trial wave
function without nodes, SHDMC uses a complex linear combination of anti-
symmetric functions without nodal pockets. The walker distribution is
prevented to fall into the bosonic ground state solution by the phase factor
of the guiding function, which remains antisymmetric, and introduces an
effective potential in the local energy. ortiz93 ; stonehenge

The present approach shares many aspects of the SHDMC method for complex wave
functions, but it also incorporates a key advantage of the CFDMC ceperley88
approach: several wave functions are optimized at the same time. In systems
where many excitations can be approximated by a single amplitude $\Phi({\bf
R})$ and a different phase factor $e^{{\bf i}\phi_{n}({\bf R})}$, a correlated
sampling approach that reweighs the walkers in Eq. (18) as $W_{i}\rightarrow
W_{i}\Phi^{T}_{n}({\bf R})/\Phi({\bf R})$ and changes the phase contribution
to the local energy in Eq. (16), would save significant time. If that
approximation were used, this generalization of the SHDMC method would look
very similar in spirit with CFDMC, the main difference being the use of a
complex guiding functions that prevents the exponential growth of the bosonic
ground state.

The displacement of each wave function in the small basis during the DMC
process is decomposed into a displacement within the subspace already
described by the other elements of the small basis plus a contribution
orthogonal to the small subspace. The displacement included within the basis
is used to improve upon the locality locality  and local-timestonehenge
approximations. The displacement orthogonal to the small subspace is used to
correct the small basis used in the next iteration.

The serial orthogonalization step required in the original SHDMC algorithm for
excited states rockandroll ; stonehenge  is avoided with a method that allows
the calculation of multiple wave functions in parallel. In addition, the
complications of inequivalent nodal pockets of excited states rockandroll  is
avoided using complex trial wave functions without nodes.

The scaling of the cost of an individual iteration of this method is
proportional to $M_{S}\times N_{c}$; $M_{S}$ the size of the small basis and
$N_{c}$ the number of configurations of the DMC run. The cost of and
individual DMC step is dependent of the number of basis functions $N_{b}$ and
electrons.

It is well known that as the size of the system increases, the number of basis
functions $N_{b}$ required to maintain a fixed error bar for a given
eigenstate must increase factorially. But in practice, the error required for
evaluation of thermodynamical properties is determined by $\beta^{-1}$: errors
must be much smaller than the temperature of the system. Therefore, as the
temperature increases and averages of multiple eigenstates are obtained, the
detail required by calculations of the ground state energies with chemical
accuracy is no longer necessary. The present approach can take advantage of
the acceleration of the algorithms used to evaluate large numbers of
determinants. nukala09 ; clark11  For a very large $N_{b}$, the cost of these
algorithms scales as $N_{e}\times N_{b}$.

The total cost is dependent on the physical system and the goal of the
calculation. If the goal is to converge the entire basis or to optimize the
Free energy, the bottleneck for convergence is the energy gap
$E_{M_{S}+1}-E_{M_{S}}$ which determines the convergence of the basis towards
the highest eigenstates considered. Accordingly, in this case, the ideal
situation for this method would be a system with ($M_{S}$) nearly degenerate
eigenstates well-separated from the rest of the spectra within an energy scale
of $k_{b}T$. If the goal, instead, is to converge the small basis so as only
the lowest $M_{L}$ eigenstates are well described, the convergence of the
algorithm is much faster and it is limited by the number of statistical
samples and the exponential decay $e^{-\beta(E_{M_{S}+1}-E_{M_{L}})}$. The
cost is reduced as compared with the calculation of eigenstates if one accepts
an error in the higher excitations. If one wishes to retain the physics of
higher eigenstates in the basis, it is computationally more efficient to
increase $M_{S}$ (which increases the cost linearly), instead of improving the
basis for the higher excitations which increases the cost exponentially.

Comparisons of the method with full CI calculations show that SHDMC can be
used to optimize many-body basis sets to maximize the overlap with the lowest
energy excitations of the Hamiltonian. Each eigenenergy obtained with this
method has lower quality than those obtained with alternative approaches such
as LMMC or the standard SHDMC for excited states. However, this method could
be a useful tool to optimize the basis, minimizing the size of the matrices
used in LMMC and thus reducing the effects of numerical noise in LMMC.

###### Acknowledgements.

The authors would like to thank J. Krogel and P. R. C. Kent for a critical
reading of the manuscript and discussions. This work has been supported by the
grant ERKCS92 Materials Science and Engineering division of Basic Energy
Sciences, Department of Energy.

##  Appendix A Going beyond the locality and local-time approximations

Approximate coefficients for the Slater determinant expansion of the
eigenstates of $\mathcal{\hat{H}}$ and $e^{-\beta\mathcal{\hat{H}}}$ can be
obtained from the eigenstates of
$\mathcal{\hat{U}}=e^{\hat{J}}e^{-\beta\mathcal{\hat{H}}}e^{-\hat{J}}$ in the
$\\{\chi^{S}_{n}\\}$ basis [see Eq. (II.5)]. But, $\mathcal{\hat{U}}$ is not
hermitian, since
$\mathcal{\hat{U}}^{\dagger}=e^{2\hat{J}}\mathcal{\hat{U}}e^{-2\hat{J}}$.
Nevertheless, as long as the Jastrow factor operator $e^{-\hat{J}}$ has an
inverse, $\mathcal{\hat{U}}$ has a set of right eigenvectors
$|\chi^{U}_{i}\rangle=e^{\hat{J}}|\Psi_{i}\rangle$ and left eigenvectors
$\langle\chi^{U}_{i}|=\langle\Psi_{i}|e^{-\hat{J}}$. Since $\mathcal{\hat{H}}$
is Hermitian, its eigenstates $|\Psi_{i}\rangle$ are orthogonal, which implies
that $\langle\chi^{U}_{i}|\chi^{U}_{j}\rangle=\delta_{i,j}$. Within
statistical error bars, in the small subspace, the matrix elements of
$\mathcal{\hat{U}}$ obtained with the basis $|\chi^{S,\ell}_{n}\rangle$ of the
iteration $\ell$ are given by Eq. (II.5). In the first iteration the matrix
elements of $\mathcal{\hat{U}}$ can be obtained directly from the Lanczos-like
procedure.

Within the small subspace, $\mathcal{\hat{U}}$ can be written as

$\displaystyle\mathcal{\hat{U}}=\sum_{i}u_{i}|\chi^{U}_{i}\rangle\langle\chi^{U}_{i}|.$
(42)

Since the $u_{i}$ are also the eigenvalues of $e^{-\beta\mathcal{\hat{H}}}$
their dependence with $\beta$ is exponential. Thus for an arbitrary
$\beta^{\prime}$ the eigenvalue will be $u_{i}^{\beta^{\prime}/\beta}$.

Provided that the difference
$|\Delta\tilde{\chi}^{S,\ell+1}_{n}\rangle=|\chi^{S,\ell+1}_{n}\rangle-|\chi^{S,\ell}_{n}\rangle$
is small [which is always valid for $\beta\rightarrow 0$ see Eq. (38)], the
dependence in $\beta$ of $|\chi^{S,\ell+1}_{n}(\beta)\rangle$ can be
approximated as follows: Let’s first define the operator
$\hat{R}=\sum_{n=0}^{M_{S}-1}|\chi^{S,\ell}_{n}\rangle\langle\chi^{S,\ell+1}_{n}|$
and its inverse within the small subspace
$\hat{R}^{\dagger}=\sum_{n=0}^{M_{S}-1}||\chi^{S,\ell+1}_{n}\rangle\langle\chi^{S,\ell}_{n}|$.

Accordingly, the dependence in $\beta^{\prime}$ of the new basis can be
approximated as

$\displaystyle\left|\chi^{S,\ell+1}_{n}\left(\beta^{\prime}\right)\right\rangle=$
$\displaystyle\hat{R}^{\dagger}\mathcal{\hat{U}}^{{\beta^{\prime}}/{\beta}}\hat{R}\left|\chi^{S,\ell+1}_{n}\right\rangle,$
(43)

with $\mathcal{\hat{U}}^{{\beta^{\prime}}/{\beta}}$ given by Eq. (42)
replacing $u_{i}$ by $u_{i}^{\beta^{\prime}/\beta}$.

##  Appendix B Working with eigenstates of $\mathcal{\hat{U}}$

While in some systems eigenstates of $\mathcal{\hat{H}}$ are always real (e.g.
confined systems with time reversal symmetry), in many cases the wave function
of the eigenstates is known to be complex, without nodal pockets. In those
cases Eq. (II.3) is valid and no walker needs to be killed or rejected because
the eigenstate wave function does not have a nodal surface. stonehenge  In
these cases it might be advantageous to propagate single eigenstates of
$\hat{R}^{\dagger}\mathcal{\hat{U}}^{\frac{j}{k}}\hat{R}$, since the variance
of the weights is minimized and lower (larger) $T$ ($\beta$) can be reached
with less statistical data. An additional advantage of working with functions
that are closer to the eigenstates is that the locality and local-time
approximations can be used.

Since any
$|\chi^{S}_{n}\rangle=\sum_{i=0}^{M_{S}-1}\langle\chi^{X}_{i}|\chi^{S}_{n}\rangle|\chi^{X}_{i}\rangle$
can be written as a linear combination of $|\chi^{X}_{i}\rangle$ and vise
versa, when eigenstates are complex one can use as trial wave function
$e^{\hat{J}}|\chi^{U}_{i}\rangle$ in the SHDMC propagation and sampling. Then
the propagation of $|\chi^{S}_{n}\rangle$ can be written as a linear
combination of the propagation of the eigenstates of $\mathcal{\hat{U}}$. as:

$\displaystyle|\delta\chi^{S}_{n}\rangle$
$\displaystyle=\sum_{i=0}^{M_{S}}\left\\{\langle\chi^{X}_{i}|\chi^{S}_{n}\rangle
e^{-\beta(E^{X}_{i}-E^{S}_{n})}\left[|\chi^{X}_{i}\rangle+|\delta\chi^{X}_{i}\rangle\right]\right\\}-|\chi^{S}_{n}\rangle\text{
and }$ $\displaystyle|\delta\tilde{\chi}^{S}_{n}\rangle$
$\displaystyle=\sum_{i=0}^{M_{S}}\left\\{\langle\chi^{X}_{i}|\chi^{S}_{n}\rangle
e^{-\beta(E^{X}_{i}-E^{S}_{n})}\left[|\chi^{X}_{i}\rangle+|\delta\tilde{\chi}^{X}_{i}\rangle\right]\right\\}-|\chi^{S}_{n}\rangle$
(44)

$E^{S}_{n}$ results from the condition
$\langle\chi^{S}_{n}|\delta\chi^{S}_{n}\rangle=0$ and it is given by

$\displaystyle e^{-\beta
E^{S}_{n}}=\sum_{i}\left[|\langle\chi^{X}_{i}|\chi^{S}_{n}\rangle|^{2}+\langle\chi^{X}_{i}|\chi^{S}_{n}\rangle\langle\chi^{S}_{n}|\delta\chi^{X}_{i}\rangle\right]e^{-\beta
E^{X}_{i}}\;,$ (45)

replacing $X$ by $U$ in Eqs. (B) and (45). The brackets
$\langle\chi^{U}_{i}|\chi^{S}_{n}\rangle$ are obtained by diagonalizing
$\mathcal{\hat{U}}$: the coefficient $n$ of the left eigenvectors $i$ in the
small basis $\\{\langle\chi^{S}_{n}|\\}$.

The disadvantages is that complex eigenstates appear only in certain
Hamiltonians or boundary conditions. Albeit without nodes, they might have
large variation in probability density, in particular for small magnetic
fields or twist boundary conditions close to high symmetry points. Large
variations in the probability density hinder correlating sampling.

##  Appendix C Working with eigenstates pairs

It is well known that in many physical systems the energy spacing between
eigenstates decreases as the size of the system increases. It is also known
that the variance of the local energy, which is related to the statistical
error in the energy, increases as the size of the system increases. nemec10
Accordingly as the size of the system increases, it becomes more difficult to
obtain eigenstates. As the size of the system increases the error in the
variance introduced by a linear combination of eigenstates in the Slater part
of the wave function becomes smaller than the variance introduced by short
range correlations. These short range correlations cannot be accounted by the
Slater part, even with a very large basis $\\{|n\rangle\\}$, or with simple
Jastrow factors. Accordingly, in large systems little is lost by using a
linear combination of eigenstates, since their contribution to the variance is
proportional to the energy separation that decreases as the system become
larger. In contrast, much is gained avoiding the divergences at the nodes by
using complex linear combinations of eigenstates, in particular, to obtain
average of thermodynamical properties. However, to propagate for larger
$\beta$ with a branching algorithm, it will be necessary to minimize the
variance of the local energy.

If the eigenstates are real, one must use a linear combination of eigenstates.
The minimum variance will be reached by constructing the Slater part of the
guiding wave function with linear combinations of eigenstates with consecutive
eigenvalues of $u_{i}$. In this work we use

$\displaystyle|\chi^{V}_{2j-1}\rangle$
$\displaystyle=\frac{1}{\sqrt{2}}\left(|\chi^{U}_{2j-1}\rangle+e^{{\bf
i}\phi}|\chi^{U}_{2j}\rangle\right)$ (46) $\displaystyle|\chi^{V}_{2j}\rangle$
$\displaystyle=\frac{1}{\sqrt{2}}\left(|\chi^{U}_{2j-1}\rangle-e^{{\bf
i}\phi}|\chi^{U}_{2j}\rangle\right),$

where $e^{{\bf i}\phi}$ is a complex phase which can be adjusted to minimize
the variance of the amplitude of the complex wave function $\Phi^{T}_{n}({\bf
R})$. The conjugate vectors $\langle\chi^{V}_{2j-1}|$ are constructed using
complex conjugate coefficients and the left eigenvectors of
$\mathcal{\hat{U}}$ in the small basis $\\{\langle\chi^{S}_{n}|\\}$.

Their evolution in imaginary time is given by

$\displaystyle|\chi^{V}_{2j-1}(\beta^{\prime})\rangle$
$\displaystyle=\frac{1}{\sqrt{2}}\left[(u_{2j-1})^{\frac{\beta^{\prime}}{\beta}}|\chi^{U}_{2j-1}\rangle+e^{{\bf
i}\phi}(u_{2j})^{\frac{\beta^{\prime}}{\beta}}|\chi^{U}_{2j}\rangle\right]$
(47) $\displaystyle|\chi^{V}_{2j}(\beta^{\prime})\rangle$
$\displaystyle=\frac{1}{\sqrt{2}}\left[(u_{2j-1})^{\frac{\beta^{\prime}}{\beta}}|\chi^{U}_{2j-1}\rangle-e^{{\bf
i}\phi}(u_{2j})^{\frac{\beta^{\prime}}{\beta}}|\chi^{U}_{2j}\rangle\right]$

## References

  * (1) C. Kittel Quantum Theory of Solids (John Wiley & Sons, New York, 1987).
  * (2) E. Runge, E. K. U. Gross, Phys. Rev. Lett. 52, 997-1000 (1984).
  * (3) G. Onida, L. Reining, A. Rubio, Rev. Mod. Phys. 74, 601-659 (2002).
  * (4) C. Filippi, M. Zaccheddu and F. Buda. Chem. Theory Comput.5, 2074 2087 (2009).
  * (5) G. Mazzola, A. Zen, and S. Sorella J. Chem. Phys. 137, 134112 (2012).
  * (6) J.C. Grossman and L. Mitas, Phys. Rev. Lett. 79, 4353-4356, (1997).
  * (7) S. Saccani, C. Filippi, S. Moroni http://arxiv.org/abs/1211.5462 (2013).
  * (8) B. G. Levine, T. J. Martinez, in Annual Review of Physical Chemistry 58, 613-634 (2007).
  * (9) M. A. Morales, C. Pierleoni, D. M. Ceperley, Phys. Rev. E. 81, 021202 (2010).
  * (10) D. M. Ceperley and B. J. Alder, Phys. Rev. Lett. 45, 566 (1980).
  * (11) P. Hohenberg and W. Kohn, Phys. Rev. 136, B864 (1964).
  * (12) W. Kohn and L. J. Sham, Phys. Rev. 140, A1133 (1965).
  * (13) N. D. Mermin, Phys. Rev. 137, A1441 (1965).
  * (14) W. Kohn and P. Vashishta. in Theory of the inhomogeneous electron gas (S. Lundqvist and N. H. March eds.) New York: Plenum, pp. 79-147 (1983).
  * (15) M. Weinert and J. W. Davenport Phys. Rev. B 45, 13709 (1992).
  * (16) J. Staunton, B. L. Gyorffy, A. J. Pindor, et al. J. of Mag. and Mag. Mat. 45, 15-22 (1984).
  * (17) B. L. Gyorffy, A. J. Pindor, J Staunton, et al. J. of Phys. F 15, 1337-1386 (1985).
  * (18) J. Staunton, B. L. Gyorffy, A. J. Pindor, et al. J. of Phys. F 15, 1387-1404 (1985).
  * (19) Y. Wang, G. M. Stocks, W. A. Shelton, D. M. C. Nicholson, Z. Szotek and W. M. Temmerman, Phys. Rev. Lett. 75, 2867 (1995).
  * (20) D. P. Landau, S.-H. Tsai, and M. Exler, American Journal of Physics, 72, 1294 1302, (2004).
  * (21) M. Eisenbach, C.-G. Zhou, D. M. Nicholson, G. Brown, J. Larkin, and T. C. Schulthess, SC ’09: Proceedings of the Conference on High Performance Computing Networking, Storage and Analysis, ACM (2009).
  * (22) J. Yin, M. Eisenbach, D. M. Nicholson, A. Rusanu, Phys. Rev. B 86, 214423 (2012).
  * (23) D. M. Ceperley, Rev. of Mod. Phys. 67, 279-355 (1995).
  * (24) W. R. Magro, D. M. Ceperley, C, Pierleoni, and B. Bernu, Phys. Rev. Lett. 76,1240-1243 (1996).
  * (25) B. Militzer, D. M. Ceperley, Phys. Rev. Lett. 85, 1890-1893 (2000).
  * (26) K.P. Driver, B. Militzer, B. , Phys. Rev. Lett. 108, 15502 (2012).
  * (27) W. Purwanto, S. Zhang, and H. Krakauer, J. Chem. Phys. 130, 094107 (2009).
  * (28) G. H. Booth, A. Gr neis, G. Kresse, and A. Alavi Nature 493, 365 370 (2012).
  * (29) W. M. C. Foulkes, L. Mitas, R. J. Needs, and G. Rajagopal, Rev. Mod. Phys. 73, 33 (2001).
  * (30) B. M. Austin, D. Y. Zubarev,W. A. Lester, Chem. Rev. 112, 263-288 (2012).
  * (31) G. Ortiz, D. M. Ceperley, and R. M. Martin, Phys Rev. Lett. 71, 2777 (1993).
  * (32) G. Ortiz, and D. M. Ceperley Phys. Rev. Lett. 75, 4642 (1995).
  * (33) M. D. Jones, G. Ortiz, and D. M. Ceperley, Phys. Rev. E, 55, 6202, (1997).
  * (34) A. D. Güçlü and C. J. Umrigar, Phys. Rev. B, 72, 045309 (2005); A. D. Güçlü, G. S. Jeon, C. J. Umrigar and J. K. Jain, Phys. Rev. B 72, 205327 (2005); G. S. Jeon, A. D. Güçlü, C. J. Umrigar, and J. K. Jain, Phys. Rev. B 72, 245312, (2005).
  * (35) C. J. Umrigar, J. Toulouse, C. Filippi, S. Sorella, and R. G. Hennig, Phys. Rev. Lett. 98, 110201 (2007).
  * (36) J. Toulouse and C. J. Umrigar, J. Chem. Phys. 128, 174101 (2008).
  * (37) F. R. Petruzielo, Julien Toulouse, C. J. Umrigar, J. of Chem. Phys. 136, 124116 (2012).
  * (38) P. Lòpez-Rios, A. Ma, N. D. Drummond, M. D. Towler, and R. J. Needs, Phys. Rev. E 74, 066701 (2006).
  * (39) D. M. Ceperley and B. Bernu, J. Chem. Phys. 89, 6316 (1988); B. Bernu, D. M. Ceperley, and W. A. Lester, Jr., J. Chem. Phys. 93, 552 (1990).
  * (40) F. A. Reboredo, R. Q. Hood, and P. R. C. Kent, Phys. Rev. B 79, 195117 (2009).
  * (41) M. Bajdich, M. L. Tiago, R. Q. Hood, P. R. C. Kent, and F. A. Reboredo, Phys. Rev. Lett. 104, 193001 (2010).
  * (42) F. A. Reboredo J. Chem. Phys. 136, 204101 (2012).
  * (43) J. B. Anderson, Int. J. Quantum Chem. 15, 109 (1979).
  * (44) D. Bressanini Phys. Rev. B 86, 115120 (2012) and references there in.
  * (45) D. M. Ceperley, J. Stat. Phys. 63, 1237 (1991).
  * (46) K. M. Rasch and L. Mitas, Chem. Phys. Lett. 528, 59 (2012).
  * (47) F. A. Reboredo, Phys. Rev. B 80, 125110 (2009).
  * (48) L. E. Reichl A modern Course in Statistical Physics (John Wiley & Sons, INC., New York,1997)
  * (49) T. Kato, Commun. Pure Appl. Math. 10, 151 (1957).
  * (50) M. Bajdich,* L. Mitas, and L. K. Wagner and K. E. Schmidt, Phys. Rev. B 77, 115112 (2008).
  * (51) J. Toulouse, M. Caffarel, P. Reinhardt, P.E. Hoggan, C. J. Umrigar, In Advances in theory of quantum systems in chemistry and physics Book Series: Progress in Theoretical Chemistry and Physics 22, 343-351 (2012).
  * (52) The theory can be extended to a back-flow operator rios06  as long at it has an inverse.
  * (53) Since we used guiding functions without nodal pockets, the complications of inequivalent nodal pockets rockandroll  are also avoided.
  * (54) The complex logarithm is a multivalued function. See Donald Sarason,”Complex function theory”, 2nd ed., Amer. Math. Society, 2007, Section IV.9. Also en.wikipedia.org/wiki/Complex_logarithm.
  * (55) B. L. Hammond, W. A. Lester, Jr., and P. J. Reynolds, Monte Carlo Methods in Ab Initio Quantum Chemistry (World Scientific, Singapore-New Jersey-London-Hong Kong, 1994).
  * (56) D. Meyer, S. Boblest, G. Wunner, Phys. Rev. A 87, 032515 (2013).
  * (57) L. Mitas, E.L. Shirley and D.M. Ceperley J. Chem. Phys. 95, 3467 (1991).
  * (58) The fixed-phase approximation ortiz93 , which corresponds to neglecting the term inside the $\\{\\}$ in Eq. (II.3), is never used in this work.
  * (59) C. J. Umrigar, M. P. Nightingale, and K. J. Runge, J. Chem. Phys. 99, 2865 (1993).
  * (60) F. A. Reboredo and P. R. C. Kent, Phys. Rev. B 77, 245110 (2008).
  * (61) Seiichiro Ten-no, Osamu Hino Int. J. Mol. Sci. 3, 459-474 (2002).
  * (62) P. K. V. V. Nukala and P. R. C. Kent, J. Chem. Phys. 130, 204105 (2009).
  * (63) Clark, B.K., Morales, M.A., McMinis, J., Kim, and G. E. Scuseria, J. of Chem. Phys. 135, 244105 (2011).
  * (64) N. Nemec, Phys. Rev. B 81, 035119 (2010).