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

almot ok with manuscript

Browse Source
This commit is contained in:
Pierre-Francois Loos 2022-10-03 22:11:01 +02:00
parent afac9706af
commit 3e4a3bec96
1 changed files with 67 additions and 104 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

171
g.tex
View File

@ -686,7 +686,7 @@ which can be rewritten probabilistically as
\end{align} \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.
\titou{In practice, the Monte Carlo algorithm is a simple generalization of that used in standard diffusion Monte Carlo calculations.} 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$. 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$.
Averages of the weight products $ \prod_{k=0}^{p-1} W_{I_k I_{k+1}} $ times the end-point contributions ${(\cH/\cS)}_{n_p,I_p}$ are computed for each $p$. Averages of the weight products $ \prod_{k=0}^{p-1} W_{I_k I_{k+1}} $ times the end-point contributions ${(\cH/\cS)}_{n_p,I_p}$ are computed for each $p$.
The generation of the paths can be performed using a two-step process. The generation of the paths can be performed using a two-step process.
@ -760,12 +760,11 @@ where $H_0$ is some arbitrary reference Hamiltonian, we have the Dyson equation
with $G^E_{0,ij} = \mel{i}{ \qty( H_0-E \Id )^{-1} }{j}$. with $G^E_{0,ij} = \mel{i}{ \qty( H_0-E \Id )^{-1} }{j}$.
Let us choose $H_0$ such that $\mel{ i }{ H_0 }{ j } = \mel{ i }{ P_i H P_i }{ j }$ for all $i$ and $j$. Let us choose $H_0$ such that $\mel{ i }{ H_0 }{ j } = \mel{ i }{ P_i H P_i }{ j }$ for all $i$ and $j$.
Then, the Dyson equation \eqref{eq:GE} becomes Then, the Dyson equation \eqref{eq:GE} becomes
\begin{multline} \be
{G}^{E}_{ij} {G}^{E}_{ij}
= {G}^{E,\cD}_{ij} = {G}^{E,\cD}_{ij}
\\
+ \sum_k \mel{ i }{ {\qty[P_i \qty(H-E \Id)P_i ]}^{-1} (H_0-H) }{ k } {G}^{E}_{kj}. + \sum_k \mel{ i }{ {\qty[P_i \qty(H-E \Id)P_i ]}^{-1} (H_0-H) }{ k } {G}^{E}_{kj}.
\end{multline} \ee
Using the following identity Using the following identity
\be \be
\begin{split} \begin{split}
@ -776,42 +775,40 @@ Using the following identity
\end{split} \end{split}
\ee \ee
the Dyson equation may be written under the form the Dyson equation may be written under the form
\begin{multline} \be
G^E_{ij} = {G}^{E,\cD}_{ij} G^E_{ij} = {G}^{E,\cD}_{ij}
\\ + \sum_{k \notin \cD_i} \mel{ i }{ {\qty[P_i \qty(H-E \Id) P_i]}^{-1} (-H)( \Id -P_i) }{ k } {G}^{E}_{kj},
+ \sum_{k \notin \cD_i} \mel{ i }{ {\qty[P_i \qty(H-E \Id) P_i]}^{-1} (-H)( \Id -P_i) }{ k } {G}^{E}_{kj} \ee
\end{multline}
which is identical to Eq.~\eqref{eq:eqfond} when $G^E_{ij}$ is expanded iteratively. which is identical to Eq.~\eqref{eq:eqfond} when $G^E_{ij}$ is expanded iteratively.
Let us use as effective transition probability density Let us use as effective transition probability density
\be \be
P_{I \to J} = \frac{1} {\PsiG(I)} \mel{ I }{ {\qty[ P_I \qty(H^+ - \EL^+ \Id) P_I]}^{-1} (-H^+) (\Id -P_I) }{ J } \PsiG(J) P_{I \to J} = \frac{1} {\PsiG(I)} \mel{ I }{ {\qty[ P_I \qty(H^+ - \EL^+ \Id) P_I]}^{-1} (-H^+) (\Id -P_I) }{ J } \PsiG(J),
\ee \ee
and the weight and the weight
\be \be
W^E_{IJ} = W^E_{IJ} =
\frac{ \mel{ I }{ {\qty[ P_I \qty(H-E \Id) P_I]}^{-1} (-H)( \Id -P_I) }{ J} } \frac{ \mel{ I }{ {\qty[ P_I \qty(H-E \Id) P_I]}^{-1} (-H)( \Id -P_I) }{ J} }
{\mel{ I }{ {\qty[ P_I \qty(H^+ - \EL^+ \Id) P_I]}^{-1} (-H^+)( \Id -P_I) }{ J} } {\mel{ I }{ {\qty[ P_I \qty(H^+ - \EL^+ \Id) P_I]}^{-1} (-H^+)( \Id -P_I) }{ J} }.
\ee \ee
Using Eqs.~\eqref{eq:eq1C}, \eqref{eq:eq3C} and \eqref{eq:relation}, one can easily verify that $P_{I \to J} \ge 0$ and $\sum_J P_{I \to J}=1$. Using Eqs.~\eqref{eq:eq1C}, \eqref{eq:eq3C} and \eqref{eq:relation}, one can easily verify that $P_{I \to J} \ge 0$ and $\sum_J P_{I \to J}=1$.
Finally, the probabilistic expression reads Finally, the probabilistic expression reads
\begin{multline} \be
\label{eq:final_E} \label{eq:final_E}
G^E_{i_0 i_N}= {G}^{E,\cD}_{i_0 i_N} G^E_{i_0 i_N}= {G}^{E,\cD}_{i_0 i_N}
+ \sum_{p=1}^{\infty} \PsiG_{i_0} \expval{ \qty( \prod_{k=0}^{p-1} W^E_{I_k I_{k+1}} ) + \sum_{p=1}^{\infty} \PsiG_{i_0} \expval{ \qty( \prod_{k=0}^{p-1} W^E_{I_k I_{k+1}} )
\frac{{G}^{E,\cD}_{I_p i_N}} { \PsiG_{I_p}} } \frac{{G}^{E,\cD}_{I_p i_N}} { \PsiG_{I_p}} }.
\end{multline} \ee
% HERE
%----------------------------% %----------------------------%
\subsubsection{Energy estimator} \subsubsection{Energy estimator}
%----------------------------% %----------------------------%
To calculate the energy, we introduce the following estimator To calculate the energy, we introduce the following estimator
\be \be
\cE(E) = \frac{ \mel{ I_0 }{ \qty(H-E \Id)^{-1} }{ H\PsiT } } {\mel{ I_0 }{ \qty(H-E \Id)^{-1} }{ \PsiT } }. \cE(E) = \frac{ \mel{ I_0 }{ \qty(H-E \Id)^{-1} }{ H\PsiT } } {\mel{ I_0 }{ \qty(H-E \Id)^{-1} }{ \PsiT } },
\ee \ee
which is the counterpart of the quantity ${\cal E}_N =\frac{ \mel{ I_0 }{ T^N }{ H\PsiT } } {\mel{ I_0 }{T^N}{ \PsiT } }$ used which is the counterpart of the quantity ${\cal E}_N =\frac{ \mel{ I_0 }{ T^N }{ H\PsiT } } {\mel{ I_0 }{T^N}{ \PsiT } }$ used
in the time-dependent formalism. In this case, the energy was easily obtained by taking the large $N$-limit of ${\cal E}_N$, see Eq.(\ref{eq:E0}). in the time-dependent formalism. In this case, the energy was easily obtained by taking the large $N$-limit of ${\cal E}_N$, see Eq.(\ref{eq:E0}).
Here, the situation is not as simple and we need a way to extract the energy from $\cE(E)$. Here, the situation is not as simple and we must find a way to extract the energy from $\cE(E)$.
Using the spectral decomposition of $H$, we have Using the spectral decomposition of $H$, we have
\be \be
@ -822,12 +819,12 @@ where $c_i = \braket{ I_0 }{ \Phi_i } \braket{ \Phi_i}{ \PsiT }$ and $\Phi_i$ a
The important observation is that for all eigenstates we have $\cE(E_i)= E_i$. Thus, to get the ground-state energy we propose to search for the solution The important observation is that for all eigenstates we have $\cE(E_i)= E_i$. Thus, to get the ground-state energy we propose to search for the solution
of the non-linear equation $\cE(E)= E$ in the vicinity of $E_0$. of the non-linear equation $\cE(E)= E$ in the vicinity of $E_0$.
In practical Monte Carlo calculations the DMC energy will be obtained by computing a finite number of components $H_p$ and $S_p$ defined as follows In practical Monte Carlo calculations, the DMC energy is obtained by computing a finite number of components $H_p$ and $S_p$ defined as follows
\be \be
\cE^\text{DMC}(E,p_{max})= \frac{ H_0+ \sum_{p=1}^{p_{max}} H^\text{DMC}_p } \cE^\text{DMC}(E,p_{max})= \frac{ H_0+ \sum_{p=1}^{p_\text{max}} H^\text{DMC}_p }
{S_{\titou{0}} +\sum_{p=1}^{p_{max}} S^\text{DMC}_p }. {S_{\titou{0}} +\sum_{p=1}^{p_\text{max}} S^\text{DMC}_p }.
\ee \ee
For $ p\ge 1$, Eq~\eqref{eq:final_E} gives For $ p\ge 1$, Eq.~\eqref{eq:final_E} gives
\begin{align} \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} }, 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} },
\\ \\
@ -847,11 +844,8 @@ For $p=0$, the two components are exactly evaluated as
S_0 & = \mel{ I_0 }{ {\qty[ P_{I_0} \qty(H-E \Id) P_{I_0} ]}^{-1} }{ \PsiT }. S_0 & = \mel{ I_0 }{ {\qty[ P_{I_0} \qty(H-E \Id) P_{I_0} ]}^{-1} }{ \PsiT }.
\end{align} \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. 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.
To avoid the calculation by Monte Carlo of a quantity, such as $H_0$ or $S_0$, which can be evaluated analytically is computationally 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.
very appealing (the statistical error is removed). We thus propose to extend 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.
the exact calculation of $H_p$ and $S_p$ to higher values of $p$,
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} %\begin{multline}
%H_p= \sum_{I_1 \notin \cD_0, \hdots , I_p \notin \cD_{p-1}} %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} } ] % \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} } ]
@ -860,16 +854,13 @@ of the number of intermediate states involved in the explicit sums makes the cal
%\end{multline} %\end{multline}
%with a similar formula for $S_p$. In practice, $\cE^\text{DMC}(E,p_{ex},p_{max})$ %with a similar formula for $S_p$. In practice, $\cE^\text{DMC}(E,p_{ex},p_{max})$
Finally, $\cE^\text{DMC}(E,p_{ex},p_{max})$ Finally, $\cE^\text{DMC}(E,p_\text{ex},p_\text{max})$ is evaluated in practice as follows
is evaluated in practice as follows
\be \be
\cE^\text{DMC}(E,p_{ex},p_{max})= \frac{ \sum_{p=0}^{p_{ex}-1} H_p +\sum_{p=p_{ex}}^{p_{max}} H^\text{DMC}_p } \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_{ex}-1} S_p +\sum_{p=p_{ex}}^{p_{max}} S^\text{DMC}_p } { \sum_{p=0}^{p_\text{ex}-1} S_p +\sum_{p=p_\text{ex}}^{p_\text{max}} S^\text{DMC}_p },
\ee \ee
where $p_{ex}$ is the number of components exactly computed. 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$, to remove the statistical error \titou{Let us emphasize that, since the magnitude of $H_p$ and $S_p$ decreases as a function of $p$, to remove the statistical error for the first most important contributions can lead to important gains, as illustrated in the numerical application to follow.}
for the first most important contributions can lead
to important gains, as illustrated in the numerical application to follow.
It is easy to check that, in the vicinity of the exact energy $E_0$, $\cE(E)$ is a linear function of $E - E_0$. It is easy to check that, in the vicinity of the exact energy $E_0$, $\cE(E)$ is a linear function of $E - E_0$.
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$. 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$.
@ -878,36 +869,6 @@ However, as $E \to E_0$, both the numerators and denominators of Eq.~\eqref{eq:c
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 of 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. 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.
%Let us describe the 3 functions used here for the fit.
%i) Linear fit\\
%We write
%\be
%\cE(E)= a_0 + a_1 E
%\ee
%and search the best value of $(a_0,a_1)$ by fitting the data.
%Then, $\cE(E)=E$ leads to
%$$
%E_0= \frac{a_0}{1-a_1}
%$$
%ii) Quadratic fit\\
%At the quadratic level we write
%$$
%\cE(E)= a_0 + a_1 E + a_2 E^2
%$$
%leading to
%$$
%E_0 = \frac{1 - a_1 \pm \sqrt{ (a_1 -1)^2 - 4 a_0 a_2}}{2 a_2}
%$$
%iii) Two-component fit\\
%We take the advantage that the exact expression of $\cE(E)$ is known as an infinite series, Eq.(\ref{calE}).
%If we limit ourselves to the first two-component, we write
%\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 variational parameters used for the fit are ($c_0,E_0,c_1,E_1)$.\\
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Numerical application to the Hubbard model} \section{Numerical application to the Hubbard model}
\label{sec:numerical} \label{sec:numerical}
@ -969,20 +930,21 @@ For the other states, we choose a single-state domain as
\ee \ee
To define $\cD$, let us introduce the following set of states To define $\cD$, let us introduce the following set of states
\be \be
{\cS_{ij} = \qty{ \ket{n} : n_D(n)=i \land n_A(n)=j }}. {\cS_{ij} = \qty{ \ket{n} : n_D(n)=i \land n_A(n)=j }},
\ee \ee
which means that $\cD$ contains the set of states having up to $n_D^\text{max}$ double occupations and, for each state with a number of double occupations equal to $m$, a number of nearest-neighbor antiparallel pairs between $n_A^\text{min}(m)$ and $n_A^\text{max}(m)$. which means that $\cD$ contains the set of states having up to $n_D^\text{max}$ double occupations and, for each state with a number of double occupations equal to $m$, a number of nearest-neighbor antiparallel pairs between $n_A^\text{min}(m)$ and $n_A^\text{max}(m)$.
Here, $n_A^\text{max}(m)$ will not be varied and taken to be the maximum possible for a given $m$, $n_A^\text{max}(m)= \max_{N}(N-1-2m,0)$. Here, $n_A^\text{max}(m)$ is fixed and set at its maximum value for a given $m$, \ie, $n_A^\text{max}(m)= \max_{N}(N-1-2m,0)$.
Using these definitions, the main domain is taken as the union of some elementary domains Using these definitions, the main domain is taken as the union of some elementary domains
\be \be
\cD = \bigcup_{n_D=0}^{n_D^\text{max}}\cD(n_D,n_A^\text{min}(n_D)), \cD = \bigcup_{n_D=0}^{n_D^\text{max}}\cD(n_D,n_A^\text{min}(n_D)),
\ee \ee
where the elementary domain is defined as where the elementary domains are defined as
\be \be
\cD(n_D,n_A^\text{min}(n_D))=\bigcup_{ j = n_A^\text{min}(n_D)}^{n_A^\text{max}(n_D)} \cS_{n_D j}. \cD(n_D,n_A^\text{min}(n_D))=\bigcup_{ j = n_A^\text{min}(n_D)}^{n_A^\text{max}(n_D)} \cS_{n_D j}.
\ee \ee
The two quantities defining the main domain are thus $n_D^\text{max}$ and $n_A^\text{min}(m)$. The two quantities defining the main domain are thus $n_D^\text{max}$ and $n_A^\text{min}(m)$.
To give an illustrative example, let us consider the 4-site case. There are 6 possible elementary domains To give an illustrative example, let us consider the 4-site case.
There are 6 possible elementary domains:
\begin{align*} \begin{align*}
\cD(0,3) & = \cS_{03}, \cD(0,3) & = \cS_{03},
& &
@ -1011,7 +973,7 @@ where
\\ \\
\cS_{20} & = \qty{ \ket{\uparrow \downarrow, \uparrow \downarrow, 0 ,0 } + \ldots }. \cS_{20} & = \qty{ \ket{\uparrow \downarrow, \uparrow \downarrow, 0 ,0 } + \ldots }.
\end{align*} \end{align*}
For the three last cases, the dots indicate the remaining states obtained by permuting the position of the pairs. For the three last cases, the dots indicate that one must also consider the remaining states obtained by permuting the position of the pairs.
%====================================== %======================================
\subsection{DMC simulations} \subsection{DMC simulations}
@ -1023,20 +985,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$. Let us begin with a small chain of 4 sites with $U=12$.
From now on, we shall take $t=1$. 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 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 of $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{T}=-0.495361\ldots$.
In what follows $\ket{I_0}$ will be systematically chosen as one of the two N\'eel states, \eg, $\ket{I_0} = \ket{\uparrow,\downarrow, \uparrow,\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_p$ 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. 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$. 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_{ex}=0$). 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, error bars are too small to be seen.
When $E$ is far from the exact value of $E_0=-0.768068\ldots$, the convergence is very rapid and only a few terms of the $p$-expansion are necessary.
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. 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.
Note the oscillations of the curves as a function of $p$ due to a parity effect specific to this system. 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. 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_{ex}=1,p_{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})$ as a function of $E$ is presented in Fig.~\ref{fig4}.
Here, $p_{max}$ is taken sufficiently large so that the convergence at large $p$ is reached. 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$. 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: For small $E$, the curve is extrapolated using the so-called two-component expression:
\be \be
@ -1076,28 +1039,28 @@ 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. Table \ref{tab1} illustrates the dependence of the Monte Carlo results upon the choice of the domain.
The reference energy is $E=-1$. 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 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 using a minimal 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. 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_{conv}$, defined as the smallest value of $p$ for which the energy is converged with 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 converged to six decimal places.
The smaller $p_{conv}$, the better the convergence is. The smaller $p_\text{conv}$, the better the convergence is.
Although this is a rough estimate, it is sufficient here for our purpose. 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}$. As clearly seen, the speed of convergence is directly related to the magnitude of $\bar{t}_{I_0}$.
The longer the stochastic trajectories remain trapped within the domain, the better the convergence. The longer the stochastic trajectories remain trapped within the domain, the better the convergence.
Of course, when the domain is chosen to be the full space, the average trapping time becomes infinite. Of course, when the domain is chosen to be the full space, the average trapping time becomes infinite.
Let us emphasize that the rate of convergence has no reason to be related to the size of the domain. Let us emphasize that the rate of convergence has no reason to be related to the size of the domain.
For example, the domain $\cD(0,3) \cup \cD(1,0)$ has a trapping time for the N\'eel state of $6.2$, while the domain $\cD(0,3) \cup \cD(1,1)$ having almost the same number of states (28 states), has an average trapping time about 6 times longer. For example, the domain $\cD(0,3) \cup \cD(1,0)$ has a trapping time for the N\'eel state of $6.2$, while the domain $\cD(0,3) \cup \cD(1,1)$, having almost the same number of states (28 states), has an average trapping time about 6 times longer.
Finally, the last column gives the energy obtained for $E=-1$. Finally, the last column gives the energy obtained for $E=-1$.
The energy is expected to be independent of the domain and to converge to a common value, which is indeed the case here. The energy is expected to be independent of the domain and to converge to a common value, which is indeed the case here.
The exact value, $\cE(E=-1)=-0.75272390\ldots$, can be found at the last row of Table \ref{tab1} for the case of a domain corresponding to the full space. The exact value, $\cE(E=-1)=-0.75272390\ldots$, can be found in the last row of Table \ref{tab1} for the case of a domain corresponding to the full space.
In sharp contrast, the statistical error depends strongly on the type of domains used. In sharp contrast, the statistical error depends strongly on the type of domains used.
As expected, the largest error of $3 \times 10^{-5}$ is obtained in the case of a single-state domain for all states. As expected, the largest error of $3 \times 10^{-5}$ is obtained in the case of a single-state domain for all states.
The smallest statistical error is obtained for the ``best'' domain having the largest average trapping time. The smallest statistical error is obtained for the ``best'' domain having the largest average trapping time.
Using this domain leads to a reduction in the statistical error as large as about three orders of magnitude, nicely illustrating the critical importance of the domains employed. Using this domain leads to a reduction in the statistical error as large as about three orders of magnitude, which nicely illustrates the critical importance of the domains employed.
%%% TABLE I %%% %%% TABLE I %%%
\begin{table} \begin{table}
\centering \centering
\caption{One-dimensional Hubbard model with $N=4$, $U=12$, $E=-1$, $\alpha=1.292$, $\beta=0.552$, and $p_{ex}=4$. \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. 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. $\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}(p)$ as a function of $p$.
@ -1128,20 +1091,20 @@ $\cD(0,1)\cup\cD(1,0)\cup\cD$(2,0) & 36 & $\infty$&1&$-0.75272390$\\
\end{table} \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 the statistical error on the largest components.
Table \ref{tab2} shows the results both for the case of a single-state main 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} 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}. Table \ref{tab2} reports the statistical fluctuations of the energy for the simulation of Table \ref{tab1}.
Results show that it is indeed interesting to compute exactly as many components as possible. Results show that it is indeed worth computing exactly as many components as possible.
For the single-state domain, a factor 2 of reduction of the statistical error is obtained when passing from no analytical computation, $p_{ex}=0$, to the case where eight components for $H_p$ and $S_p$ are exactly computed. 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 best domain, the impact is much more important with a huge reduction of about three orders of magnitude in the statistical error. 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 converged as a function of $p$ with their statistical error on the last digit for $E=
Table \ref{tab3} reports the energies convergence as a function of $p$ alongside their statistical error on the last digit for $E=
-0.8$, $-0.795$, $-0.79$, $-0.785$, and $-0.78$. -0.8$, $-0.795$, $-0.79$, $-0.785$, and $-0.78$.
The values are displayed in Fig.~\ref{fig4}. The values are displayed in Fig.~\ref{fig4}.
As seen on the figure the behavior of $\cE$ as a function of $E$ is very close to the linearity. As seen, the behavior of $\cE$ as a function of $E$ is very linear.
The extrapolated values obtained from the five values of the energy with the three fitting functions are reported. 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 compare with the exact value of $-0.768068\ldots$. Using the linear fitting 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. A small bias of about $4 \times 10^{-5}$ is observed.
This bias vanishes within the statistical error with the quadratic fit. This bias vanishes within the statistical error when one relies on the quadratic and two-component fitting functions.
It is also true when using the two-component function.
Our final value is in full agreement with the exact value with about six decimal places. 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. 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.
@ -1149,17 +1112,17 @@ We also report the variational and exact energies together with the values of th
As $U$ increases the configurations with zero or one double occupation become more and more predominant and the average trapping time increases. As $U$ increases the configurations with zero or one double occupation become more and more predominant and the average trapping time increases.
For very large values of $U$ (say, $U \ge 12$) the increase of $\bar{t}_{I_0}$ becomes particularly steep. For very large values of $U$ (say, $U \ge 12$) the increase of $\bar{t}_{I_0}$ becomes particularly steep.
Finally, in Table \ref{tab4} we report the results obtained for larger systems at $U=12$ for a size of the chain ranging from $N=4$ (36 states) Finally, in Table \ref{tab4}, we report the results obtained for larger systems at $U=12$ for a chain size ranging from $N=4$ (36 states)
to $N=12$ ($\sim 10^6$ states). to $N=12$ ($\sim 10^6$ states).
No careful construction of domains maximizing the average trapping time has been performed, we have merely chosen domains of reasonable size (not more than 2682) by taking not too large number of double occupations (only, $n_D=0,1$) and not too small number of nearest-neighbor antiparallel pairs. No careful construction of domains maximizing the average trapping time has been performed, we have merely chosen domains of reasonable size (no more than 2682) by taking not too large number of double occupations (only, $n_D=0,1$) and not too small number of nearest-neighbor antiparallel pairs.
As seen, as the number of sites increases, the average trapping time for the domains chosen decreases. 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. 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 called for. 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 DMC 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 for the reference energy. 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. 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 size increases. However, an increase of the statistical error is observed when the system size increases.
To get lower error bars we need to use better trial wave functions, better domains, and also larger simulation times. To get lower error bars, a 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. 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. All these aspects will be considered in a forthcoming work.
@ -1167,7 +1130,7 @@ All these aspects will be considered in a forthcoming work.
\begin{table} \begin{table}
\caption{One-dimensional Hubbard model with $N=4$, $U=12$, and $E=-1$. \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 calculated analytically.
The simulation is performed the same way as Table \ref{tab1}. 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.} 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} \label{tab2}
\begin{ruledtabular} \begin{ruledtabular}
@ -1189,7 +1152,7 @@ $8$ & $2.2 \times10^{-5}$ &$ 0.05 \times 10^{-8}$\\
%%% TABLE III %%% %%% TABLE III %%%
\begin{table} \begin{table}
\caption{One-dimensional Hubbard model with $N=4$, $U=12$, $\alpha=1.292$, $\beta=0.552$, and $p_{ex}=4$. \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 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 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$.}
@ -1218,7 +1181,7 @@ The main domain is $\cD(0,1) \cup \cD(1,0)$.}
\label{tab4} \label{tab4}
\begin{ruledtabular} \begin{ruledtabular}
\begin{tabular}{rcccr} \begin{tabular}{rcccr}
$U$ & $(\alpha,\beta)$ & $E_\text{T}$ & $E_{ex}$ & $\bar{t}_{I_0}$ \\ $U$ & $(\alpha,\beta)$ & $E_\text{T}$ & $E_\text{ex}$ & $\bar{t}_{I_0}$ \\
\hline \hline
8 & (0.908,\,0.520) & $-0.770342\ldots$ &$-1.117172\ldots$ & 33.5\\ 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\\ 10 & (1.116,\,0.539) & $-0.604162\ldots$ &$-0.911497\ldots$ & 63.3\\
@ -1234,11 +1197,11 @@ $U$ & $(\alpha,\beta)$ & $E_\text{T}$ & $E_{ex}$ & $\bar{t}_{I_0}$ \\
%%% TABLE V %%% %%% TABLE V %%%
\begin{table*} \begin{table*}
\caption{One-dimensional Hubbard model with $U=12$. \caption{One-dimensional Hubbard model with $U=12$.
The fits to extrapolate the DMC energies are done using the two-component function} The extrapolated DMC energies are obtained using the two-component fitting function.}
\label{tab5} \label{tab5}
\begin{ruledtabular} \begin{ruledtabular}
\begin{tabular}{crcrccccr} \begin{tabular}{crcrccccr}
$N$ & Hilbert space size & Domain & Domain size & $(\alpha,\beta)$ &$\bar{t}_{I_0}$ & $E_\text{T}$ & $E_{ex}$ & $E^\text{DMC}$\\ $N$ & Hilbert space size & Domain & Domain size & $(\alpha,\beta)$ &$\bar{t}_{I_0}$ & $E_\text{T}$ & $E_\text{ex}$ & $E^\text{DMC}$\\
\hline \hline
4 & 36 & $\cD(0,1) \cup \cD(1,0)$ & 30 &(1.292,\,0.552)& 108.7 & $-0.495361$ & $-0.768068$ & $-0.7680676(5)$\\ 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)$\\ 6 & 400 & $\cD(0,1) \cup \cD(1,0)$ &200 &(1.124,\,0.689)&57.8 & $-0.633297$ & $-1.215395$& $-1.215389(9)$\\