loos/DQMC
1
0
Fork 1
Code Issues Pull Requests Projects Releases Wiki Activity

revision Michel

Browse Source
This commit is contained in:
Pierre-Francois Loos 2022-12-16 13:47:59 +01:00
parent fbf6bf5607
commit 0d964aaae6
2 changed files with 165 additions and 149 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

129
g.bib
View File

@ -1,75 +1,73 @@
%% This BibTeX bibliography file was created using BibDesk.
%% http://bibdesk.sourceforge.net/
%% Created for Pierre-Francois Loos at 2022-09-30 16:13:18 +0200
%% Created for Pierre-Francois Loos at 2022-12-16 13:28:28 +0100
%% Saved with string encoding Unicode (UTF-8) *
%% Saved with string encoding Unicode (UTF-8)
@article{Ceperley_1983,
title = {The simulation of quantum systems with random walks: A new algorithm for charged systems},
journal = {Journal of Computational Physics},
volume = {51},
number = {3},
pages = {404-422},
year = {1983},
issn = {0021-9991},
doi = {https://doi.org/10.1016/0021-9991(83)90161-4},
url = {https://www.sciencedirect.com/science/article/pii/0021999183901614},
author = {D Ceperley},
abstract = {Random walks with branching have been used to calculate exact properties of the ground state of quantum many-body systems. In this paper, a more general Green's function identity is derived which relates the potential energy, a trial wavefunction, and a trial density matrix to the rules of a branched random walk. It is shown that an efficient algorithm requires a good trial wavefunction, a good trial density matrix, and a good sampling of this density matrix. An accurate density matrix is constructed for Coulomb systems using the path integral formula. The random walks from this new algorithm diffuse through phase space an order of magnitude faster than the previous Green's Function Monte Carlo method. In contrast to the simple diffusion Monte Carlo algorithm, it is an exact method. Representative results are presented for several molecules.}
}
abstract = {Random walks with branching have been used to calculate exact properties of the ground state of quantum many-body systems. In this paper, a more general Green's function identity is derived which relates the potential energy, a trial wavefunction, and a trial density matrix to the rules of a branched random walk. It is shown that an efficient algorithm requires a good trial wavefunction, a good trial density matrix, and a good sampling of this density matrix. An accurate density matrix is constructed for Coulomb systems using the path integral formula. The random walks from this new algorithm diffuse through phase space an order of magnitude faster than the previous Green's Function Monte Carlo method. In contrast to the simple diffusion Monte Carlo algorithm, it is an exact method. Representative results are presented for several molecules.},
author = {D Ceperley},
date-modified = {2022-12-16 13:14:50 +0100},
doi = {https://doi.org/10.1016/0021-9991(83)90161-4},
issn = {0021-9991},
journal = {J. Comput. Phys.},
number = {3},
pages = {404-422},
title = {The simulation of quantum systems with random walks: A new algorithm for charged systems},
url = {https://www.sciencedirect.com/science/article/pii/0021999183901614},
volume = {51},
year = {1983},
bdsk-url-1 = {https://www.sciencedirect.com/science/article/pii/0021999183901614},
bdsk-url-2 = {https://doi.org/10.1016/0021-9991(83)90161-4}}
@article{Kalos_1962,
title = {Monte Carlo Calculations of the Ground State of Three- and Four-Body Nuclei},
author = {Kalos, M. H.},
journal = {Phys. Rev.},
volume = {128},
issue = {4},
pages = {1791--1795},
numpages = {0},
year = {1962},
month = {Nov},
publisher = {American Physical Society},
doi = {10.1103/PhysRev.128.1791},
url = {https://link.aps.org/doi/10.1103/PhysRev.128.1791}
}
author = {Kalos, M. H.},
doi = {10.1103/PhysRev.128.1791},
issue = {4},
journal = {Phys. Rev.},
month = {Nov},
numpages = {0},
pages = {1791--1795},
publisher = {American Physical Society},
title = {Monte Carlo Calculations of the Ground State of Three- and Four-Body Nuclei},
url = {https://link.aps.org/doi/10.1103/PhysRev.128.1791},
volume = {128},
year = {1962},
bdsk-url-1 = {https://link.aps.org/doi/10.1103/PhysRev.128.1791},
bdsk-url-2 = {https://doi.org/10.1103/PhysRev.128.1791}}
@article{Kalos_1970,
title = {Energy of a Boson Fluid with Lennard-Jones Potentials},
author = {Kalos, M. H.},
journal = {Phys. Rev. A},
volume = {2},
issue = {1},
pages = {250--255},
numpages = {0},
year = {1970},
month = {Jul},
publisher = {American Physical Society},
doi = {10.1103/PhysRevA.2.250},
url = {https://link.aps.org/doi/10.1103/PhysRevA.2.250}
}
author = {Kalos, M. H.},
doi = {10.1103/PhysRevA.2.250},
issue = {1},
journal = {Phys. Rev. A},
month = {Jul},
numpages = {0},
pages = {250--255},
publisher = {American Physical Society},
title = {Energy of a Boson Fluid with Lennard-Jones Potentials},
url = {https://link.aps.org/doi/10.1103/PhysRevA.2.250},
volume = {2},
year = {1970},
bdsk-url-1 = {https://link.aps.org/doi/10.1103/PhysRevA.2.250},
bdsk-url-2 = {https://doi.org/10.1103/PhysRevA.2.250}}
@article{Moskowitz_1986,
author = {Moskowitz,Jules W. and Schmidt,K. E. },
title = {The domain Green’s function method},
journal = {The Journal of Chemical Physics},
volume = {85},
number = {5},
pages = {2868-2874},
year = {1986},
doi = {10.1063/1.451046},
URL = {
https://doi.org/10.1063/1.451046
},
eprint = {
https://doi.org/10.1063/1.451046
}
}
author = {Moskowitz,Jules W. and Schmidt,K. E.},
date-modified = {2022-12-16 13:14:34 +0100},
doi = {10.1063/1.451046},
journal = {J. Chem. Phys.},
number = {5},
pages = {2868-2874},
title = {The domain Green's function method},
volume = {85},
year = {1986},
bdsk-url-1 = {https://doi.org/10.1063/1.451046}}
@article{Willow_2012,
author = {Willow,Soohaeng Yoo and Kim,Kwang S. and Hirata,So},
date-modified = {2022-09-30 16:10:46 +0200},
@ -339,13 +337,12 @@ eprint = {
year = {2012}}
@book{Ceperley_1979,
author = {D.M. Ceperley and M.H Kalos},
editor={K.Binder},
chapter={4},
publisher = {Springer, Berlin},
title = {Monte Carlo Methods in Statistical Physics},
year = {1979}}
author = {D.M. Ceperley and M.H Kalos},
chapter = {4},
editor = {K.Binder},
publisher = {Springer, Berlin},
title = {Monte Carlo Methods in Statistical Physics},
year = {1979}}
@article{Carlson_2007,
author = {J. Carlson},
@ -1698,10 +1695,10 @@ eprint = {
@article{Assaraf_2000,
author = {Assaraf, Roland and Caffarel, Michel and Khelif, Anatole},
date-modified = {2022-09-15 10:17:18 +0200},
date-modified = {2022-12-16 13:13:55 +0100},
doi = {10.1103/physreve.61.4566},
issn = {1095-3787},
journal = {Phsys. Rev. E},
journal = {Phys. Rev. E},
month = {Apr},
number = {4},
pages = {4566--4575},

185
g.tex
View File

@ -127,7 +127,7 @@ breaklinks=true
\author{Anthony Scemama}
\affiliation{\LCPQ}
\author{Michel Caffarel}
\email{Corresponding author: caffarel@irsamc.ups-tlse.fr}
\email[Corresponding author: ]{caffarel@irsamc.ups-tlse.fr}
\affiliation{\LCPQ}
\begin{abstract}
@ -139,7 +139,7 @@ Here, we extend this idea to the general case of a walker trapped within domains
The equations of the resulting effective stochastic dynamics are derived.
The larger the average (trapping) time spent by the walker within the domains, the greater the reduction in statistical fluctuations.
A numerical application to the Hubbard model is presented.
Although this work presents the method for finite linear spaces, it can be generalized without fundamental difficulties to continuous configuration spaces.
Although this work presents the method for (discrete) finite linear spaces, it can be generalized without fundamental difficulties to continuous configuration spaces.
\end{abstract}
\maketitle
@ -154,7 +154,7 @@ DMC can be used either for systems defined in a continuous configuration space (
Here, we shall consider only the discrete case, that is, the general problem of calculating the lowest eigenvalue and/or eigenstate of a (very large) matrix.
The generalization to continuous configuration spaces presents no fundamental difficulty.
In essence, DMC is based on \textit{stochastic} power methods, a family of well established numerical approaches able to extract the largest or smallest eigenvalues of a matrix (see, \eg, Ref.~\onlinecite{Golub_2012}).
In essence, DMC is based on \textit{stochastic} power methods, a family of well-established numerical approaches able to extract the largest or smallest eigenvalues of a matrix (see, \eg, Ref.~\onlinecite{Golub_2012}).
These approaches are particularly simple as they merely consist in applying a given matrix (or some simple function of it) as many times as required on some arbitrary vector belonging to the linear space.
Thus, the basic step of the corresponding algorithm essentially reduces to successive matrix-vector multiplications.
In practice, power methods are employed under more sophisticated implementations, such as, \eg, the Lancz\`os algorithm (based on Krylov subspaces) \cite{Golub_2012} or Davidson's method where a diagonal preconditioning is performed. \cite{Davidson_1975}
@ -176,18 +176,17 @@ It is shown how to define an effective stochastic dynamics describing walkers mo
The equations of the effective dynamics are derived and a numerical application for a model (one-dimensional) problem is presented.
In particular, it shows that the statistical convergence of the energy can be greatly enhanced when domains associated with large average trapping times are considered.
It should be noted that the use of domains in quantum Monte Carlo is not new. Domains have been introduced within the context of the
Green's function Monte Carlo (GFMC) method pioneered by Kalos\cite{Kalos_1962,Kalos_1970}
and later developped and applied by Kalos and others.\cite{Kalos_1974,Ceperley_1979,Ceperley_1983,Moskowitz_1986}
In GFMC, an approximate Green's function that can be sampled is needed to propagate stochastically the wavefunction.
In the so-called Domain's GFMC version of GFMC introduced in Ref.\onlinecite{Kalos_1970} and \onlinecite{Kalos_1974}
It should be noted that the use of domains in quantum Monte Carlo is not new. Domains have been introduced within the context of
Green's function Monte Carlo (GFMC) pioneered by Kalos \cite{Kalos_1962,Kalos_1970}
and later developed and applied by Kalos and others. \cite{Kalos_1974,Ceperley_1979,Ceperley_1983,Moskowitz_1986}
In GFMC, an approximate Green's function that can be sampled is required for the stochastic propagation of the wave function.
In the so-called domain GFMC version of GFMC introduced in Ref.~\onlinecite{Kalos_1970} and \onlinecite{Kalos_1974}
the sampling is realized by using the restriction of the Green's function
to a small domain consisting of the Cartesian product of small spheres around each particle, the potential being considered as constant within the domain.
Fundamentally, the method presented in this work is closely related to the Domain's GFMC, although the way we present the formalism in terms of walkers
to a small domain consisting of the cartesian product of small spheres around each particle, the potential being considered constant within the domain. Fundamentally, the method presented in this work is closely related to the domain GFMC, although the way we present the formalism in terms of walkers
trapped within domains
and derive the equations may appear as different.
However, a key difference here is that we show how to use
domains of arbitrary size, a point which can greatly enhance the efficiency of the simulations when domains are suitably chosen, as illustrated in our numerical
and derive the equations that may appear different.
However, we show here how to use
domains of arbitrary size, a new feature that greatly enhances the efficiency of the simulations when domains are suitably chosen, as illustrated in our numerical
application.
@ -339,7 +338,7 @@ as it should.
The rewriting of $\bar{T}_{ij}$ as a product of a stochastic matrix times some general real weight
does not introduce any constraint on the choice of the stochastic matrix, so that, in theory, any stochastic matrix could be used. However,
in practice, it is highly desirable that the magnitude of the fluctuations of the weight during the Monte Carlo simulation be as small as
possible. A natural solution is to choose a stochastic matrix equal to a good approximation of $\bar{T}_{ij}$. This is done as follows.
possible. A natural solution is to choose a stochastic matrix as close as possible to $\bar{T}_{ij}$. This is done as follows.
Let us introduce the following operator
\be
@ -404,14 +403,14 @@ Note that one can eschew this condition via a simple generalization of the trans
\ee
This new transition probability matrix with positive entries reduces to Eq.~\eqref{eq:pij} when $T^+_{ij}$ is positive as $\sum_j \PsiG_{j} T^+_{ij} = 1$.
Now, we need to make the connection between the transition probability matrix, $p_{i \to j}$, defined from the approximate Hamiltonian $H^{+}$ via $T^+$ and the operator $T$ associated with the exact Hamiltonian $H$.
Now, we need to make the connection between the transition probability matrix, $p_{i \to j}$, defined from the Hamiltonian $H^{+}$ via $T^+$ and the operator $T$ associated with the exact Hamiltonian $H$.
This is done thanks to Eq.~\eqref{eq:defTij} that connects $p_{i \to j}$ and $T_{ij}$ through the weight
\be
w_{ij}=\frac{T_{ij}}{T^+_{ij}},
\ee
derived from Eqs.~\eqref{eq:defT} and \eqref{eq:pij}.
To calculate the probabilistic averages, an artificial (mathematical) ``particle'' called walker (or psi-particle) is introduced.
To calculate the probabilistic averages, an artificial (mathematical) ``particle'' called a walker (or psi-particle) is introduced.
During the Monte Carlo simulation, the walker moves in configuration space by drawing new states with
probability $p_{i_k \to i_{k+1}}$, thus realizing the path of probability $\text{Prob}_{i_0}$.
Note that, instead of using a single walker, it is common to introduce a population of independent walkers and to calculate the averages over this population.
@ -439,8 +438,8 @@ We shall not insist here on these practical details that are discussed, for exam
\label{sec:single_domains}
%=======================================%
During the simulation, walkers move from state to state with the possibility of being trapped a certain number of times on the same state before
exiting to a different state. This fact can be exploited in order to integrate out some part of the dynamics, thus leading to a reduction of the statistical
During the simulation, walkers move from state to state with the possibility of being trapped a certain number of times in the same state before
exiting to a different state. This fact can be exploited in order to integrate out some parts of the dynamics, thus leading to a reduction of the statistical
fluctuations. This idea was proposed some time ago and applied to the SU($N$) one-dimensional Hubbard model.\cite{Assaraf_1999A,Assaraf_1999B,Caffarel_2000}
Considering a given state $\ket{i}$, the probability that a walker remains exactly $n$ times in $\ket{i}$ (with $ n \ge 1 $) and then exits to a different state $j$ (with $j \neq i$) is
@ -541,7 +540,11 @@ Since the sums are restricted to states belonging to the domain, it is convenien
\label{eq:pi}
P_I = \sum_{\ket{i} \in \cD_I} \dyad{i}{i},
\ee
as well as the projection of $T^+$ over $\cD_I$, \ie, $T^+_I = P_I T^+ P_I$, which governs the dynamics of the walkers trapped in this domain.
as well as the projection of $T^+$ over $\cD_I$
\be
T^+_I = P_I T^+ P_I,
\ee
which governs the dynamics of the walkers trapped in this domain.
%see Eq.~\eqref{eq:pij} where $T^+$ is now restricted to the domain.
Using Eqs.~\eqref{eq:pij} and \eqref{eq:eq1C}, the probability can be rewritten as
\be
@ -559,7 +562,7 @@ where the operator $F^+_I = P_I T^+ (1-P_I)$, corresponding to the last move co
\ee
Physically, $F$ may be seen as a flux operator through the boundary of $\cD_{I}$.
Knowing the probability of remaining $n$ times in the domain and, then, to exit to some state, it is possible to obtain
Knowing the probability of remaining $n$ times in the domain and, then, exiting to some state, it is possible to obtain
the probability of being trapped $n$ times in $\cD_{I}$, just by summing over all possible exit states:
\be
\label{eq:PiN}
@ -583,7 +586,7 @@ The average trapping time defined as ${\bar t}_{I}={\bar n}_{I} \tau$ where $ {\
{\bar t}_{I}=\frac{1}{\PsiG_I} \mel{I}{ { \qty[ P_I \qty( H^+ - \EL^+ \Id ) P_I ] }^{-1} }{ \PsiG }.
\ee
In practice, the various quantities restricted to the domain will be computed by inverting the matrix $(H^+-\EL^+ \Id)$ in $\cD_{I}$.
Note that it is possible only if the dimension of the domains is not too large (say, less than a few thousands).
Note that it is possible only if the dimension of the domains is not too large (say, less than a few thousand).
%=======================================%
\subsection{Time-dependent Green's matrix using domains}
@ -594,23 +597,40 @@ Note that it is possible only if the dimension of the domains is not too large (
%\subsubsection{Time-dependent Green's matrix}
%\label{sec:time}
%--------------------------------------------%
In this section we generalize the path-integral expression of the Green's matrix, Eq.~\eqref{eq:G}, to the case where domains are used.
To do so, we introduce the Green's matrix associated with each domain as follows:
In this section, we generalize the path-integral expression of the Green's matrix, Eq.~\eqref{eq:G}, to the case where domains are used.
To do so, we need to introduce the Green's matrix restricted to each domain as follows:
\be
G^{(N),\cD}_{ij}= \mel{ i }{ T_i^N }{ j }.
\ee
%Starting from Eq.~\eqref{eq:cn}
%\be
%G^{(N)}_{i_0 i_N}= \sum_{i_1,...,i_{N-1}} \prod_{k=0}^{N-1} \langle i_k| T |i_{k+1} \rangle.
%\ee
Starting from Eq.~\eqref{eq:cn} and using the representation of the paths in terms of exit states and trapping times [see Eq.~\eqref{eq:eff_series}], we get
where $T_i$ is the projection of the operator $T$ over the domain $\cD_i$ of $\ket{i}$
\be
T_i = P_i T^+ P_i,
\ee
and $P_i$ is the projector, $P_i=\sum_{\ket{k} \in \cD_i} \dyad{k}{k}$, Eq.~\eqref{eq:pi}.
Using the representation of the paths in terms of exit states and trapping times [see Eq.~\eqref{eq:eff_series}],
the set of paths can be partitioned according to the number $p$ of exit states.
As already noted above, the number of exit states ranges from $p=0$ (no exit of the initial domain) to $p=N$ (exit of the current domain at each time step).
We shall denote $\cC_p$ the set of paths with $p$ exit states.
The time-dependent Green's matrix, Eq.~\eqref{eq:cn}, is then rewritten as
\be
G^{(N)}_{i_0 i_N} = \sum_{p=0}^N
\sum_{\cC_p} \sum_{(i_1,...,i_{N-1}) \in \cC_p}
\prod_{k=0}^{N-1} \mel{ i_k }{ T }{ i_{k+1} } ,
\prod_{k=0}^{N-1} \mel{ i_k }{ T }{ i_{k+1} } .
\ee
where $\cC_p$ is the set of paths with $p$ exit states, $\ket{I_k}$, and trapping times $n_k$ with the constraints that $\ket{I_k} \notin \cD_{k-1}$ (with $1 \le n_k \le N+1$ and $\sum_{k=0}^p n_k= N+1$).
It follows that
$\cC_p$ can be partitioned further by considering the subset of paths, denoted by $\cC_p^{I,n}$,
having $(\ket{I_k}, n_k)$ for $k=0$ to $p$ as exit states and trapping times.
We recall that, by definition of the exit states, $\ket{I_k} \notin \cD_{I_{k-1}}$. The total time must be conserved, so the relation $\sum_{k=0}^p n_k= N+1$
must be fulfilled.
Now, the contribution of $\cC_p^{I,n}$ to the path integral is obtained by summing over all elementary paths $(i_1,...,i_{N-1})$
of this set, thus giving
\begin{multline}
\sum_{(i_1,...,i_{N-1}) \in \cC_p^{I,n}} \prod_{k=0}^{N-1} \mel{ i_k }{ T }{ i_{k+1} }
\\
= \qty[ \prod_{k=0}^{p-1} \mel{ I_k }{ \qty(T_{I_k})^{n_k-1} F_{I_k} }{ I_{k+1} } ]
\mel{ I_k }{ \qty(T_{I_p})^{n_p} }{ i_N }
\end{multline}
Finally, collecting all contributions the Green's matrix is written as
\begin{multline}
\label{eq:Gt}
G^{(N)}_{i_0 i_N}= G^{(N),\cD}_{i_0 i_N} +
@ -680,7 +700,7 @@ where $\expval{\cdots}_p$ is the probabilistic average defined over the set of p
\\
\cS_{n_p,I_p} & = \frac{1}{\PsiG_{I_p}} \sum_{i_N} {G}^{(n_p-1), \cD}_{I_p i_N} (\Psi_T)_{i_N},
\end{align}
two quantities taking into account the contribution of the trial wave function at the end of the path.
two quantities taking into account the contribution of the trial wave function at the end of the path.
In practice, the present Monte Carlo algorithm is a simple generalization of the standard diffusion Monte Carlo algorithm.
Stochastic paths starting at $\ket{I_0}$ are generated using the probability $\cP_{I_k \to I_{k+1}}(n_k)$ and are stopped when $\sum_k n_k$ is greater than some target value $N$.
@ -706,7 +726,7 @@ First, the integer $n_k$ is drawn using the probability $P_{I_k}(n)$ [see Eq.~\e
The aim of this section is to show that it is possible to go further by integrating out the trapping times, $n_k$, of the preceding expressions, thus defining a new effective stochastic dynamics involving now only the exit states.
Physically, it means that we are going to compute exactly the time evolution of all stochastic paths trapped in each domain.
We shall present two different ways to derive the new dynamics and renormalized probabilistic averages.
The first one, called the pedestrian way, consists in starting from the preceding time expression for $G$ and make the explicit integration over the $n_k$'s.
The first one, called the pedestrian way, consists in starting from the preceding time expression for $G$ and making the explicit integration over the $n_k$'s.
The second, more direct and elegant, is based on the Dyson equation.
@ -719,8 +739,8 @@ Let us define the energy-dependent Green's matrix
\ee
The denomination ``energy-dependent'' is chosen here since
this quantity is the discrete version of the Laplace transform of the time-dependent Green's function in a continuous space,
usually known under this name.\footnote{As $\tau \rightarrow 0$ and $N \rightarrow \infty$ with $N\tau=t$, the operator $T^N$ converges to $e^{-t(H-E \Id)}$. We then have $G^E_{ij} \rightarrow \int_0^{\infty} dt \mel{i}{e^{-t(H-E \Id)}}{j}$, which is the Laplace transform of the time-dependent Green's function $\mel{i}{e^{-t(H-E \Id)}}{j}$.}
The remarkable property is that, thanks to the summation over $N$ up to infinity, the constrained multiple sums appearing in Eq.~\eqref{eq:Gt} can be factorized in terms of a product of unconstrained sums, as follows
usually known under this name.\footnote{As $\tau \to 0$ and $N \to \infty$ with $N\tau=t$, the operator $T^N$ converges to $e^{-t(H-E \Id)}$. We then have $G^E_{ij} \to \int_0^{\infty} dt \mel{i}{e^{-t(H-E \Id)}}{j}$, which is the Laplace transform of the time-dependent Green's function $\mel{i}{e^{-t(H-E \Id)}}{j}$.}
The remarkable property is that thanks to the summation over $N$ up to infinity, the constrained multiple sums appearing in Eq.~\eqref{eq:Gt} can be factorized in terms of a product of unconstrained sums, as follows
\begin{multline}
\sum_{N=1}^\infty \sum_{p=1}^N \sum_{n_0 \ge 1} \cdots \sum_{n_p \ge 1} \delta_{\sum_{k=0}^p n_k,N+1} F(n_0,\ldots,n_N)
\\
@ -823,8 +843,10 @@ In practical Monte Carlo calculations, the DMC energy is obtained by computing a
For $ p\ge 1$, Eq.~\eqref{eq:final_E} gives
\begin{align}
H^\text{DMC}_p & = \PsiG_{i_0}\expval{ \qty(\prod_{k=0}^{p-1} W^E_{I_k I_{k+1}}) {\cal H}_{I_p} },
\label{eq:defHp}
\\
S^\text{DMC}_p & = \PsiG_{i_0} \expval{ \qty(\prod_{k=0}^{p-1} W^E_{I_k I_{k+1}}) {\cal S}_{I_p} },
\label{eq:defSp}
\end{align}
where
\begin{align}
@ -841,7 +863,7 @@ For $p=0$, the two components are exactly evaluated as
\end{align}
Note that the evaluation of $(H_0,S_0)$ is possible as long as the size of the domain is small enough to allow the calculation of the inverse matrix.
Avoiding the stochastic calculation of quantities, such as $H_0$ or $S_0$, that can be evaluated analytically is computationally very appealing as the statistical error associated with these quantities is completely removed.
We thus suggest to extend the exact calculation of $H_p$ and $S_p$ to higher $p$ values, up to the point where the exponential increase of the number of intermediate states involved in the explicit sums makes the calculation unfeasible.
We thus suggest extending the exact calculation of $H_p$ and $S_p$ to higher $p$ values, up to the point where the exponential increase of the number of intermediate states involved in the explicit sums makes the calculation unfeasible.
%\begin{multline}
%H_p= \sum_{I_1 \notin \cD_0, \hdots , I_p \notin \cD_{p-1}}
% \qty[ \prod_{k=0}^{p-1} \mel{ I_k }{ {\qty[ P_k \qty( H-E \Id ) P_k ] }^{-1} (-H)(\Id-P_k) }{ I_{k+1} } ]
@ -854,6 +876,7 @@ Finally, $\cE^\text{DMC}(E,p_\text{ex},p_\text{max})$ is evaluated in practice a
\be
\cE^\text{DMC}(E,p_\text{ex},p_\text{max})= \frac{ \sum_{p=0}^{p_\text{ex}-1} H_p +\sum_{p=p_\text{ex}}^{p_\text{max}} H^\text{DMC}_p }
{ \sum_{p=0}^{p_\text{ex}-1} S_p +\sum_{p=p_\text{ex}}^{p_\text{max}} S^\text{DMC}_p },
\label{eq:E_pex}
\ee
where $p_\text{ex}$ is the number of components computed exactly.
Let us emphasize that, since the magnitude of $H_p$ and $S_p$ decreases as a function of $p$, most of the statistical error is removed by computing the dominant terms analytically.
@ -863,7 +886,7 @@ It is easy to check that, in the vicinity of the exact energy $E_0$, $\cE(E)$ is
Therefore, in practice, we compute the value of $\cE(E)$ for several values of $E$, and fit these data using a linear, quadratic or a more complicated function of $E$ in order to obtain, via extrapolation, an estimate of $E_0$.
In order to have a precise extrapolation of the energy, it is best to compute $\cE(E)$ for values of $E$ as close as possible to the exact energy.
However, as $E \to E_0$, both the numerators and denominators of Eq.~\eqref{eq:calE} diverge.
This is reflected by the fact that one needs to compute more and more $p$-components with an important increase of statistical fluctuations.
This is reflected by the fact that one needs to compute more and more $p$-components with an important increase in statistical fluctuations.
Thus, from a practical point of view, a trade-off has to be found between the quality of the extrapolation and the amount of statistical fluctuations.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
@ -880,7 +903,7 @@ Let us consider the one-dimensional Hubbard Hamiltonian for a chain of $N$ sites
H= -t \sum_{\expval{ i j } \sigma} \hat{a}^+_{i\sigma} \hat{a}_{j\sigma}
+ U \sum_i \hat{n}_{i\uparrow} \hat{n}_{i\downarrow},
\ee
where $\expval{ i j }$ denotes the summation over two neighboring sites, $\hat{a}^+_{i\sigma}$ ($\hat{a}_{i\sigma}$) is the fermionic creation (annihilation) operator of an spin-$\sigma$ electron (with $\sigma$ = $\uparrow$ or $\downarrow$) on site $i$, $\hat{n}_{i\sigma} = \hat{a}^+_{i\sigma} \hat{a}_{i\sigma}$ the number operator, $t$ the hopping amplitude, and $U$ the on-site Coulomb repulsion.
where $\expval{ i j }$ denotes the summation over two neighboring sites, $\hat{a}^+_{i\sigma}$ ($\hat{a}_{i\sigma}$) is the fermionic creation (annihilation) operator of a spin-$\sigma$ electron (with $\sigma$ = $\uparrow$ or $\downarrow$) on site $i$, $\hat{n}_{i\sigma} = \hat{a}^+_{i\sigma} \hat{a}_{i\sigma}$ the number operator, $t$ the hopping amplitude, and $U$ the on-site Coulomb repulsion.
We consider a chain with an even number of sites and open boundary conditions at half-filling, that is, $N_{\uparrow}=N_{\downarrow}=N/2$.
In the site representation, a general vector of the Hilbert space can be written as
\be
@ -890,7 +913,7 @@ where $n_{i \sigma}=0$ or $1$ is the number of electrons of spin $\sigma$ on sit
For the one-dimensional Hubbard model with open boundary conditions, the components of the ground-state vector have the same sign (say, $c_i > 0$).
It is then possible to equate the guiding and trial vectors, \ie, $\ket{c^+} = \ket{c_\text{T}}$.
As trial wave function, we shall employ a generalization of the Gutzwiller wave function \cite{Gutzwiller_1963}
As a trial wave function, we shall employ a generalization of the Gutzwiller wave function \cite{Gutzwiller_1963}
\be
\braket{ n }{ c_\text{T} }= e^{-\alpha n_D(n)-\beta n_A(n)},
\ee
@ -900,7 +923,7 @@ where $n_D(n)$ is the number of doubly occupied sites for the configuration $\ke
\ee
The parameters $\alpha$ and $\beta$ are optimized by minimizing the variational energy of the trial wave function, \ie,
\be
E_\text{T}(\alpha,\beta) = \frac{\mel{ c_\text{T} }{H }{c_\text{T}}} {\braket{ c_\text{T} }{ c_\text{T} }}.
E_\text{v}(\alpha,\beta) = \frac{\mel{ c_\text{T} }{H }{c_\text{T}}} {\braket{ c_\text{T} }{ c_\text{T} }}.
\ee
%======================================
@ -982,20 +1005,21 @@ In what follows, we shall restrict ourselves to the case of the Green's function
Let us begin with a small chain of 4 sites with $U=12$.
From now on, we shall take $t=1$.
The size of the linear space is ${\binom{4}{2}}^2 = 36$ and the ground-state energy obtained by exact diagonalization is $E_0=-0.768068\ldots$.
The two variational parameters of the trial vector have been optimized and fixed at the values of $\alpha=1.292$ and $\beta=0.552$ with a variational energy $E_\text{T}=-0.495361\ldots$.
The two variational parameters of the trial vector have been optimized and fixed at the values of $\alpha=1.292$ and $\beta=0.552$ with a
variational energy $E_\text{v}=-0.495361\ldots$.
In what follows, $\ket{I_0}$ is systematically chosen as one of the two N\'eel states, \eg, $\ket{I_0} = \ket{\uparrow,\downarrow, \uparrow,\ldots}$.
Figure \ref{fig3} shows the convergence of $H_p$ as a function of $p$ for different values of the reference energy $E$.
Figure \ref{fig3} shows the convergence of $H^{DMC}_p$, Eq.~\eqref{eq:defHp}, as a function of $p$ for different values of the reference energy $E$.
We consider the simplest case where a single-state domain is associated to each state.
Five different values of $E$ have been chosen, namely $E=-1.6$, $-1.2$, $-1.0$, $-0.9$, and $-0.8$.
Only $H_0$ is computed analytically ($p_\text{ex}=0$).
At the scale of the figure, error bars are too small to be seen.
At the scale of the figure, the error bars are too small to be seen.
When $E$ is far from the exact value, the convergence is very rapid and only a few terms of the $p$-expansion are necessary.
In contrast, when $E$ approaches the exact energy, a slower convergence is observed, as expected from the divergence of the matrix elements of the Green's matrix at $E=E_0$ where the expansion does not converge at all.
The oscillations of the curves as a function of $p$ are due to a parity effect specific to this system.
In practice, it is not too much of a problem since a smoothly convergent behavior is nevertheless observed for the even- and odd-parity curves.
The ratio, $\cE^\text{DMC}(E,p_\text{ex}=1,p_\text{max})$ as a function of $E$ is presented in Fig.~\ref{fig4}.
The ratio, $\cE^\text{DMC}(E,p_\text{ex}=1,p_\text{max})$, Eq.~\eqref{eq:E_pex}, as a function of $E$ is presented in Fig.~\ref{fig4}.
Here, $p_\text{max}$ is taken sufficiently large so that the convergence at large $p$ is reached.
The values of $E$ are $-0.780$, $-0.790$, $-0.785$, $-0.780$, and $-0.775$.
For small $E$, the curve is extrapolated using the so-called two-component expression:
@ -1008,7 +1032,7 @@ which considers only the first two terms of the exact expression of $\cE(E)$ [se
%\be
%\cE(E) = \frac{ \frac{E_0 c_0}{E_0-E} + \frac{E_1 c_1}{E_1-E}}{\frac{c_0}{E_0-E} + \frac{c_1}{E_1-E} }
%\ee
Here, the fitting parameters that needs to be determined are $c_0$, $\epsilon_0$, $c_1$, and $\epsilon_1$.
Here, the fitting parameters that need to be determined are $c_0$, $\epsilon_0$, $c_1$, and $\epsilon_1$.
The estimate of the energy obtained from $\cE(E)=E$ is $-0.76807(5)$ in full agreement with the exact value of $-0.768068\ldots$.
%%% FIG 1 %%%
@ -1036,9 +1060,10 @@ Error bars are smaller than the symbol size.}
Table \ref{tab1} illustrates the dependence of the Monte Carlo results upon the choice of the domain.
The reference energy is $E=-1$.
The first column indicates the various domains consisting of the union of some elementary domains as explained above.
The first line of the table gives the results when one uses a single-state domain for all states, and the last one for the maximal domain containing the full linear space.
The first line of the table gives the results when one uses a single-state domain for all states and the last one for the maximal domain containing the full linear space.
The size of the various domains is given in the second column, the average trapping time for the state $\ket{I_0}$ in the third column, and an estimate of the speed of convergence of the $p$-expansion for the energy in the fourth column.
To quantify the rate of convergence, we report the quantity, $p_\text{conv}$, defined as the smallest value of $p$ for which the energy is converged to six decimal places.
To quantify the rate of convergence, we report the quantity, $p_\text{conv}$, defined as the smallest value of $p$ for which the energy is roughly stabilized
to five decimal places.
The smaller $p_\text{conv}$, the better the convergence is.
Although this is a rough estimate, it is sufficient here for our purpose.
As clearly seen, the speed of convergence is directly related to the magnitude of $\bar{t}_{I_0}$.
@ -1060,12 +1085,13 @@ Using this domain leads to a reduction in the statistical error as large as abou
\caption{One-dimensional Hubbard model with $N=4$, $U=12$, $E=-1$, $\alpha=1.292$, $\beta=0.552$, and $p_\text{ex}=4$.
The simulation is performed with 20 independent blocks and $10^5$ stochastic paths starting from the N\'eel state.
$\bar{t}_{I_0}$ is the average trapping time for the N\'eel state.
$p_\text{conv}$ is a measure of the convergence of $\cE^\text{DMC}(p)$ as a function of $p$.
$p_\text{conv}$ is a measure of the convergence of $\cE^\text{DMC}$ as a function of $p$
(smallest value of $p$ for which the energy is roughly stabilized to five decimal places)
See text for more details.}
\label{tab1}
\begin{ruledtabular}
\begin{tabular}{lrrrl}
Domain & Size & $\bar{t}_{I_0}$ & $p_\text{conv}$ & \mcc{$\cE^\text{DMC}$} \\
Domain & Size & $\bar{t}_{I_0}$ & $p_\text{conv}$ & \mcc{$\cE^\text{DMC}(E=-1)$} \\
\hline
Single & 1 & 0.026 & 88 &$-0.75276(3)$\\
$\cD(0,3)$ & 2 & 2.1 & 110 &$-0.75276(3)$\\
@ -1087,11 +1113,12 @@ $\cD(0,1)\cup\cD(1,0)\cup\cD$(2,0) & 36 & $\infty$&1&$-0.75272390$\\
\end{ruledtabular}
\end{table}
As explained above, it is very advantageous to calculate exactly as many $(H_p,S_p)$ as possible in order to avoid the statistical error on the largest components.
As explained above, it is very advantageous to calculate exactly as many $(H_p,S_p)$ as possible in order to avoid statistical error on the largest components.
Ideally, this should be done up to the value of $p$ for which the calculation of these quantities whose cost increases exponentially is possible.
Table \ref{tab2} shows the results both for the case of a single-state domain and for the domain having the largest average trapping time, namely $\cD(0,1) \cup \cD(1,1)$ (see Table \ref{tab1}).
Table \ref{tab2} reports the statistical fluctuations of the energy for the simulation of Table \ref{tab1}.
Results show that it is indeed worth computing exactly as many components as possible.
For the single-state domain, the statistical error is reduced by a factor two when passing from no analytical computation, $p_\text{ex}=0$, to the case where eight components for $H_p$ and $S_p$ are computed exactly.
For the single-state domain, the statistical error is reduced by a factor of two when passing from no analytical computation, $p_\text{ex}=0$, to the case where eight components for $H_p$ and $S_p$ are computed exactly.
For the best domain, the impact is much more important with a huge reduction of about three orders of magnitude in the statistical error.
Table \ref{tab3} reports the energies convergence as a function of $p$ alongside their statistical error on the last digit for $E=
@ -1099,9 +1126,10 @@ Table \ref{tab3} reports the energies convergence as a function of $p$ alongside
The values are displayed in Fig.~\ref{fig4}.
As seen, the behavior of $\cE$ as a function of $E$ is almost perfectly linear.
The extrapolated values obtained from the five values of the energy with three different fitting functions are reported.
Using the linear fitting function leads to an energy of $-0.7680282(5)$ (to be compared with the exact value of $-0.768068\ldots$).
Fitting the data using a simple linear function leads to an energy of $-0.7680282(5)$ (to be compared with the exact value of $-0.768068\ldots$).
A small bias of about $4 \times 10^{-5}$ is observed.
This bias vanishes within the statistical error when one relies on the quadratic and two-component fitting functions.
This bias vanishes within the statistical error when resorting to more flexible fitting functions, such as a quadratic function of $E$ or the
two-component representation given by Eq.(\ref{eq:2comp}).
Our final value is in full agreement with the exact value with about six decimal places.
Table \ref{tab4} shows the evolution of the average trapping times and extrapolated energies as a function of $U$ when using $\cD(0,1) \cup \cD(1,0)$ as the main domain.
@ -1115,10 +1143,10 @@ No careful construction of domains maximizing the average trapping time has been
As seen, as the number of sites increases, the average trapping time for the chosen domains decreases.
This point is, of course, undesirable since the efficiency of the approach may gradually deteriorate when considering large systems.
A more elaborate way of constructing domains is clearly desirable in this case.
The exact DMC energies extrapolated using the two-component function are also reported.
The exact energies extrapolated using the two-component function are also reported.
Similarly to what has been done for $N=4$, the extrapolation is performed using about five values of the reference energy.
The extrapolated DMC energies are in full agreement with the exact value within error bars.
However, an increase of the statistical error is observed when the system size increases.
The extrapolated DMC energies are in full agreement with the exact value within the error bar.
However, an increase in statistical error is observed when the system size increases.
To get lower error bars, more accurate trial wave functions may be considered, better domains, and also larger simulation times.
Of course, it will also be particularly interesting to take advantage of the fully parallelizable character of the algorithm to get much lower error bars.
All these aspects will be considered in a forthcoming work.
@ -1126,7 +1154,8 @@ All these aspects will be considered in a forthcoming work.
%%% TABLE II %%%
\begin{table}
\caption{One-dimensional Hubbard model with $N=4$, $U=12$, and $E=-1$.
Dependence of the statistical error on the energy with the number of $p$-components calculated analytically.
Dependence of the statistical error on the energy with the number of $p$-components $p_{ex}$ calculated analytically in the expression
of the energy, $\cE^\text{DMC}(E,p_\text{ex},p_\text{max})$, Eq.~\eqref{eq:E_pex}.
The simulation is performed the same way as in Table \ref{tab1}.
Results are presented when a single-state domain is used for all states and when $\cD(0,1) \cup \cD(1,0)$ is chosen as the main domain.}
\label{tab2}
@ -1152,11 +1181,11 @@ $8$ & $2.2 \times10^{-5}$ &$ 0.05 \times 10^{-8}$\\
\caption{One-dimensional Hubbard model with $N=4$, $U=12$, $\alpha=1.292$, $\beta=0.552$, and $p_\text{ex}=4$.
The main domain is $\cD(0,1) \cup \cD(1,0)$.
The simulation is performed with 20 independent blocks and $10^6$ stochastic paths starting from the N\'eel state.
The fits are performed with the five values of $E$.}
The fits are performed with the five values of $E$ reported in this table.}
\label{tab3}
\begin{ruledtabular}
\begin{tabular}{lc}
$E$ & $E^\text{DMC}$ \\
$E$ & $\cE^\text{DMC}(E)$ \\
\hline
$-0.8$ &$-0.7654686(2)$\\
$-0.795$&$-0.7658622(2)$\\
@ -1174,12 +1203,14 @@ $E_0$ exact & $-0.768068\ldots$\\
%%% TABLE IV %%%
\begin{table}
\caption{One-dimensional Hubbard model with $N=4$. Average trapping time as a
function of $U$.
function of $U$. $(\alpha,\beta)$ and $E_\text{v}$ are the parameters and variational energy
of the trial wave function, respectively. $E_\text{ex}$ is the exact energy obtained by diagonalization of the
Hamiltonian matrix.
The main domain is $\cD(0,1) \cup \cD(1,0)$.}
\label{tab4}
\begin{ruledtabular}
\begin{tabular}{rcccr}
$U$ & $(\alpha,\beta)$ & $E_\text{T}$ & $E_\text{ex}$ & $\bar{t}_{I_0}$ \\
$U$ & $(\alpha,\beta)$ & $E_\text{v}$ & $E_\text{ex}$ & $\bar{t}_{I_0}$ \\
\hline
8 & (0.908,\,0.520) & $-0.770342\ldots$ &$-1.117172\ldots$ & 33.5\\
10 & (1.116,\,0.539) & $-0.604162\ldots$ &$-0.911497\ldots$ & 63.3\\
@ -1194,18 +1225,20 @@ $U$ & $(\alpha,\beta)$ & $E_\text{T}$ & $E_\text{ex}$ & $\bar{t}_{I_0}$ \\
%%% TABLE V %%%
\begin{table*}
\caption{One-dimensional Hubbard model with $U=12$.
The extrapolated DMC energies are obtained using the two-component fitting function.}
\caption{Ground-state energy of the one-dimensional Hubbard model for different sizes and $U=12$.
The parameters $(\alpha,\beta)$ of the trial wave function and the corresponding variational energy $E_\text{v}$ are reported.
$\bar{t}_{I_0}$ is the average trapping time for the N\'eel state.
The extrapolated DMC energies, $E_\text{DMC}$, are obtained using the two-component fitting function, Eq.~\eqref{eq:2comp}. $E_\text{DMC}$ values are in full agreement with the exact values $E_\text{ex}$ computed by diagonalization of the Hamiltonian matrix.}
\label{tab5}
\begin{ruledtabular}
\begin{tabular}{crcrccccr}
$N$ & Hilbert space size & Domain & Domain size & $(\alpha,\beta)$ &$\bar{t}_{I_0}$ & $E_\text{T}$ & $E_\text{ex}$ & $E^\text{DMC}$\\
$N$ & Hilbert space size & Domain & Domain size & $(\alpha,\beta)$ &$\bar{t}_{I_0}$ & $E_\text{v}$ & $E_\text{DMC}$ & $ E_\text{ex}$\\
\hline
4 & 36 & $\cD(0,1) \cup \cD(1,0)$ & 30 &(1.292,\,0.552)& 108.7 & $-0.495361$ & $-0.768068$ & $-0.7680676(5)$\\
6 & 400 & $\cD(0,1) \cup \cD(1,0)$ &200 &(1.124,\,0.689)&57.8 & $-0.633297$ & $-1.215395$& $-1.215389(9)$\\
8 & 4 900 & $\cD(0,1) \cup \cD(1,0)$ & 1 190 &(0.984,\,0.788)& 42.8 & $-0.750995 $& $-1.66395$& $-1.6637(2)$\\
10 & 63 504 & $\cD(0,5) \cup \cD(1,4)$ & 2 682 &(0.856,\,0.869)& 31.0 & $-0.855958$ & $-2.113089$& $-2.1120(7)$\\
12 & 853 776 & $\cD(0,8) \cup \cD(1,7)$ & 1 674 &(0.739,\,0.938)& 16.7 & $-0.952127$ & $-2.562529$& $-2.560(6)$\\
4 & 36 & $\cD(0,1) \cup \cD(1,0)$ & 30 &(1.292,\,0.552)& 108.7 & $-0.495361$ & $-0.7680676(5)$ & $-0.768068$ \\
6 & 400 & $\cD(0,1) \cup \cD(1,0)$ &200 &(1.124,\,0.689)&57.8 & $-0.633297$ & $-1.215389(9)$ & $-1.215395$ \\
8 & 4 900 & $\cD(0,1) \cup \cD(1,0)$ & 1 190 &(0.984,\,0.788)& 42.8 & $-0.750995 $ & $-1.6637(2)$ & $-1.66395$ \\
10 & 63 504 & $\cD(0,5) \cup \cD(1,4)$ & 2 682 &(0.856,\,0.869)& 31.0 & $-0.855958$ & $-2.1120(7)$ & $-2.113089$ \\
12 & 853 776 & $\cD(0,8) \cup \cD(1,7)$ & 1 674 &(0.739,\,0.938)& 16.7 & $-0.952127$ & $-2.560(6)$ & $-2.562529$ \\
\end{tabular}
\end{ruledtabular}
\end{table*}
@ -1294,16 +1327,6 @@ and
\end{align}
Without loss of generality, let us choose, for example, $\ket{I_0} = \ket{1}$.
Then, it is easy to show that
%\begin{align}
% \cI_1 & = A_1 C_2,
% &
% \cI_2 & = A_1 A_2 C_1,
% &
% \cI_3 & = A_1 A_2 A_1 C_2,
% &
% \qq{etc.}
%\end{align}
%that is
\begin{align}
\cI_{2k+1} & = \frac{C_2}{A_2} (A_1 A_2)^{2k+1},
&
@ -1340,10 +1363,6 @@ Using Eqs.~\eqref{eq:defA1} and \eqref{eq:defA2}, one gets, for $i = 1$ or $2$,
which finally yields
\be
\mel{ 1 }{ \qty(H-E \Id)^{-1} }{ \Psi}
% = \qty(H_{11}-E)^{-1} \Psi_1
% - H_{12} \frac{ \qty( \Psi_2 - \frac{ H_{12}\Psi_1}{H_{11}-E} )}
% { (H_{11}-E)(H_{22}-E) - H^2_{12}}
% \\
= \frac{ H_{22}-E}{\Delta} \Psi_1 - \frac{H_{12}}{\Delta} \Psi_2,
\label{eq:final}
\ee