Quantum Monte Carlo with variable spins: fixed-phase and fixed-node approximations

TL;DR

This paper explores the fixed-phase spin-orbit diffusion Monte Carlo method, demonstrating its relation to the fixed-node approach, its ability to handle complex wave functions, and its potential for improved accuracy in electronic structure calculations.

Contribution

It introduces a fixed-phase approach that generalizes the fixed-node method, allowing for complex wave functions and potentially more accurate electronic structure calculations.

Findings

Fixed-phase method recovers fixed-node results with similar biases.

Complex wave functions can be constructed without assigning spin labels.

Fixed-phase approach smooths out divergences at wave function nodes.

Abstract

We study several aspects of the recently introduced fixed-phase spin-orbit diffusion Monte Carlo (FPSODMC) method, in particular, its relation to the fixed-node method and its potential use as a general approach for electronic structure calculations. We illustrate constructions of spinor-based wave functions with the full space-spin symmetry without assigning up or down spin labels to particular electrons, effectively "complexifying" even ordinary real-valued wave functions. Interestingly, with proper choice of the simulation parameters and spin variables, such fixed-phase calculations enable one to reach also the fixed-node limit. The fixed-phase solution provides a straightforward interpretation as the lowest bosonic state in a given effective potential generated by the many-body approximate phase. In addition, the divergences present at real wave function nodes are smoothed out to lower dimensionality, decreasing thus the variation of sampled quantities and making the sampling also more straightforward. We illustrate some of these properties on calculations of selected first-row systems that recover the fixed-node results with quantitatively similar levels of the corresponding biases. At the same time, the fixed-phase approach opens new possibilities for more general trial wave functions with further opportunities for increasing accuracy in practical calculations.

Tables

Table 1. Table 1: Total energies in Ha for the first-row elements using FN(FP) DMC with HF nodes (phases). FN calculations are extrapolated to zero time step. FP calculations take a spatial time step of 0.001 and decrease the spin time step until the energy is unchanged. NRL is the estimated nonrelativistic exact energy.

Atom	HF	NRLDavidson et al. (1991)	FN	FP
Li	$-7.43272$	$-7.47806$	$-7.47794\left(2\right)$	$-7.47804\left(7\right)$
Be	$-14.57302$	$-14.66736$	$-14.65720\left(6\right)$	$-14.6574\left(2\right)$
B	$-24.52906$	$-24.65393$	$-24.64030\left(9\right)$	$-24.64016\left(8\right)$
C	$-37.68861$	$-37.84500$	$-37.8300\left(3\right)$	$-37.8291\left(3\right)$
N	$-54.40093$	$-54.58930$	$-54.5750\left(6\right)$	$-54.5754\left(1\right)$
O	$-74.80939$	$-75.06700$	$-75.049\left(1\right)$	$-75.0513\left(1\right)$
F	$-99.40934$	$-99.73400$	$-99.7164\left(6\right)$	$-99.7175\left(1\right)$

Table 2. Table 2: Equilibrium bond lengths ( r e subscript 𝑟 𝑒 r_{e} ), dissociation energies ( D e subscript 𝐷 𝑒 D_{e} ) and vibrational frequencies ( ν 0 subscript 𝜈 0 \nu_{0} ) for the various approximations compared to experiment using a HF nodal surface/phase. Parameters and uncertainties are obtained from a fit to the Morse potential.

Method	$r_e$ (Å)	$D_e$ (eV)	${\nu }_0$ (cm-1)
FN	1.0895(8)	9.616(5)	2402(28)
FP	1.0879(7)	9.555(6)	2396(23)
Expt.nis	1.098	9.758(6)	2358.57(9)

Equations

$$-\frac{\partial \rho (\mathbf{R},\tau )}{\partial \tau }$$

$$-\frac{\partial \Phi (\mathbf{R},\tau )}{\partial \tau }$$

$$\Phi (\mathbf{R},\tau )\stackrel{!}{=}{\Phi }_T(\mathbf{R}).$$

$$\Gamma =\left\{\mathbf{R}\in {\mathbb{R}}^{dN}|{\Psi }_T(\mathbf{R})=0\right\}.$$

$${\tilde{\Psi }}_T={\Psi }_T+i\epsilon {\Psi }_B$$

$$V_{ph}(\mathbf{R})=V_{\infty }\delta (\mathbf{R}-{\mathbf{R}}_{\Gamma })$$

$$\chi (\mathbf{r},s)=\alpha {\phi }^{\uparrow }(\mathbf{r}){\chi }^{\uparrow }(s)+\beta {\phi }^{\downarrow }(\mathbf{r}){\chi }^{\downarrow }(s)$$

$$\begin{array}[]{l r}\langle s_{j}|\chi^{\uparrow}\rangle=e^{is_{j}},&\langle s_{j}|\chi^{\downarrow}\rangle=e^{-is_{j}}\end{array}$$

$$g({\mathbf{R}}^{\prime },t+\tau )=\int \text{d}\mathbf{R}\frac{{\rho }_T({\mathbf{R}}^{\prime })}{{\rho }_T(\mathbf{R})}G({\mathbf{R}}^{\prime }\leftarrow \mathbf{R},\tau )g(\mathbf{R},t)$$

$$T_i^s=-\frac{1}{2{\mu }_s}\left[\frac{{\partial }^2}{\partial s_i^2}+1\right].$$

$$\tilde{G}({\mathbf{X}}^{\prime }\leftarrow \mathbf{X};\tau )\simeq T_{{\mathbf{X}}^{\prime },\mathbf{X}}{e}^{-\tau [E_L(\mathbf{X})+E_L({\mathbf{X}}^{\prime })-2E_T]/2}$$

$$T_{{\mathbf{X}}^{\prime },\mathbf{X}}$$

$${\Psi }_T(\mathbf{X})=e^{U(\mathbf{R})}\sum _{\alpha }{c}_{\alpha }{\text{det}}_{\alpha }\left[\dots ,{\chi }_i({\mathbf{r}}_k,s_k),\dots \right].$$

$${\chi }^{+}$$

$${\chi }^{-}$$

$$\det \left[\begin{array}{ccccc}{\chi }_1^{+}(1)& {\chi }_1^{+}(2)& {\chi }_1^{+}(3)& {\chi }_1^{+}(4)& ...\\ {\chi }_1^{-}(1)& {\chi }_1^{-}(2)& {\chi }_1^{-}(3)& {\chi }_1^{-}(4)& ...\\ {\chi }_2^{+}(1)& {\chi }_2^{+}(2)& {\chi }_2^{+}(3)& {\chi }_2^{+}(4)& ...\\ {\chi }_2^{-}(1)& {\chi }_2^{-}(2)& {\chi }_2^{-}(3)& {\chi }_2^{-}(4)& ...\\ & ...& & & \end{array}\right].$$

$${\chi }^{+}$$

$${\chi }^{-}$$

$$\det \left[\begin{array}{ccccc}{\phi }_1(1)e^{i{s}_1}& {\phi }_1(2)e^{i{s}_2}& {\phi }_1(3)e^{i{s}_3}& {\phi }_1(4)e^{i{s}_4}& ...\\ {\phi }_1(1)e^{-i{s}_1}& {\phi }_1(2)e^{-i{s}_2}& {\phi }_1(3)e^{-i{s}_3}& {\phi }_1(4)e^{-i{s}_4}& ...\\ {\phi }_2(1)e^{i{s}_1}& {\phi }_2(2)e^{i{s}_2}& {\phi }_2(3)e^{i{s}_3}& {\phi }_2(4)e^{i{s}_4}& ...\\ {\phi }_2(1)e^{-i{s}_1}& {\phi }_2(2)e^{-i{s}_2}& {\phi }_2(3)e^{-i{s}_3}& {\phi }_2(4)e^{-i{s}_4}& ...\\ & ...& & & \end{array}\right]$$

$$\det \left[\begin{array}{ccccc}{\phi }_1(1)e^{is}& {\phi }_1(2)e^{i{s}^{\prime }}& {\phi }_1(3)e^{is}& {\phi }_1(4)e^{i{s}^{\prime }}& ...\\ {\phi }_1(1)e^{-is}& {\phi }_1(2)e^{-i{s}^{\prime }}& {\phi }_1(3)e^{-is}& {\phi }_1(4)e^{-i{s}^{\prime }}& ...\\ {\phi }_2(1)e^{is}& {\phi }_2(2)e^{i{s}^{\prime }}& {\phi }_2(3)e^{is}& {\phi }_2(4)e^{i{s}^{\prime }}& ...\\ {\phi }_2(1)e^{-is}& {\phi }_2(2)e^{-i{s}^{\prime }}& {\phi }_2(3)e^{-is}& {\phi }_2(4)e^{-i{s}^{\prime }}& ...\\ & ...& & & \end{array}\right]$$

$$\det \left[\begin{array}{ccccc}0& c_0{\phi }_1(2)& 0& c_0{\phi }_1(4)& ...\\ {\phi }_1(1)e^{-is}& {\phi }_1(2)e^{-i{s}^{\prime }}& {\phi }_1(3)e^{-is}& {\phi }_1(4)e^{-i{s}^{\prime }}& ...\\ 0& c_0{\phi }_2(2)& 0& c_0{\phi }_2(4)& ...\\ {\phi }_2(1)e^{-is}& {\phi }_2(2)e^{-i{s}^{\prime }}& {\phi }_2(3)e^{-is}& {\phi }_2(4)e^{-i{s}^{\prime }}& ...\\ & ...& & & \end{array}\right]$$

$$c_0=[e^{i{s}^{\prime }}-e^{i(2s-s^{\prime })}]=e^{is}[e^{i(s^{\prime }-s)}-e^{-i(s^{\prime }-s)}]=2i{e}^{is}\sin (s^{\prime }-s).$$

$$\propto [\sin (s^{\prime }-s){]}^{N/2}\det \left[\begin{array}{ccccc}{\phi }_1(1)& {\phi }_1(3)& 0& 0& ...\\ {\phi }_2(1)& {\phi }_2(3)& 0& 0& ...\\ 0& 0& {\phi }_1(2)& {\phi }_1(4)& ...\\ 0& 0& {\phi }_2(2)& {\phi }_2(4)& ...\\ & ...& & & \end{array}\right].$$

$$\Psi (1,2,3)$$

$${\Psi }_T(1,2,3)={\text{det}}_{1,2}^{\uparrow }[1s,2s]{\text{det}}_3^{\downarrow }[1s]$$

$$\Psi$$

$$\Psi (1,2)=r_1{r}_{2}g(r_1)g(r_2)[Y_{11}(1)-Y_{11}(2)]$$

$$V_{ph}=\frac{1}{2[(x_1-x_2{)}^2+(y_1-y_2{)}^2]}$$

$$V_{ph}=(1/2)[\nabla {\Phi }_T{]}^2$$

$$(T+V)\Psi =E\Psi \to (T+V+V_{ph})\rho =E\rho$$

$$|{\Psi }_T⟩=c_0|1s^{2}2s^2⟩+c_1\sum _{i\in \{x,y,z\}}|1s^{2}2p_i^2⟩$$

$$V(r)=D_e\left[e^{-2a(r-r_e)}-2e^{-a(r-r_e)}\right]$$

Full text

I Introduction

Recently, we introduced a projector quantum Monte Carlo (QMC) method for calculating quantum systems with both spatial and spin degrees of freedom Melton et al. (2016). The approach is based on an overcomplete representation for spin variables such that the sampling is similar to the spatial variables. Given our choice of spin representation, the method involves the fixed-phase Ortiz, Ceperley, and Martin (1993) approximation, hence its acronym fixed-phase spin-orbit/spinor diffusion Monte Carlo (FPSODMC). This approach enabled us to carry out QMC calculations of atoms and molecules with spin-orbit interactions in the spinor formalism including cases where high accuracy was needed for both spin-orbit and electron correlation effects. In a subsequent work we explored simple cases of fixed-phase vs. fixed-node Foulkes et al. (2001) approximations in order to compare the corresponding biases in these two related possibilities Melton, Bennett, and Mitas (2016); Melton and Mitas . We constructed simple cases where both fixed-phase and fixed-node conditions were equivalent or very similar and we found comparable biases in the total energies using the two approximations.

In this work we explore this direction further by investigating a clear unification and smooth transition between these approaches. It has been known for some time that the fixed-node approximation is a special case of the fixed-phase approximation. Our method makes this relationship explicit through the construction of trial wave functions that in a particular limit recover the fixed-node trial function. We use this property for QMC fixed-phase calculations of several systems and we directly compare the fixed-phase biases to the corresponding fixed-node biases. In addition, we explicitly show how one can obtain the fixed-node result as a limit of the fixed-phase calculation. We elaborate on the advantages and disadvantages of the fixed-phase approach as far as further QMC developments are concerned.

II Fixed-Phase Spinor Diffusion Monte Carlo

Let us briefly outline the key notions of the FPSODMC approach: fixed-phase approximation, continuous spin representation and the corresponding importance sampling approach.

II.1 Fixed-Phase Approximation

For complex wave functions we present a brief sketch of the fixed-phase method (FPDMC) Ortiz, Ceperley, and Martin (1993) and its relation to the fixed-node flavor of DMC.

Let us consider the Born-Oppenheimer Hamiltonian $H=-(1/2)\nabla^{2}+V(\mathbf{R})$, where $\nabla=(\nabla_{1},\nabla_{2},\ldots,\nabla_{N})$ and $V$ denotes the electron-ion and electron-electron Coulomb interactions. We denote spatial configurations as $\mathbf{R}=(\mathbf{r}_{1},\mathbf{r}_{2},\ldots,\mathbf{r}_{N})\in\mathbb{R}^{dN}$, where $N$ is the number of particles and $d$ is dimensionality (here we assume $d=3$). We assume a complex wave function $\Psi(\mathbf{R},\tau)=\rho(\mathbf{R},\tau)e^{i\Phi(\mathbf{R},\tau)}$ and substitute it into the imaginary-time Schrödinger equation. For the amplitude $\rho(\mathbf{R},\tau)$ and phase $\Phi(\mathbf{R},\tau)$ we obtain

[TABLE]

where $T_{kin}=-(1/2)\nabla^{2}$. The fixed-phase approximation is given by imposing the phase to be equal to the phase of trial or variational wave function $\Psi_{T}({\bf R})=\rho_{T}({\bf R})e^{i\Phi_{T}({\bf R})}$ that is independent of $\tau$

[TABLE]

so that the second equation is not considered any further. On the other hand, the stationary trial phase enables us to solve the equation for the non-negative amplitude $\rho$ and the corresponding energy eigenvalue. Clearly, both are now dependent on the trial phase through the additional potential $V_{ph}=1/2|\nabla\Phi_{T}(\mathbf{R})|^{2}$.

II.2 Fixed-phase upper bound property

The fixed-phase approximation is variational since the repulsive potential $V_{ph}$ can only raise the energy for an approximate phase Ortiz, Ceperley, and Martin (1993). This is easy to see from the energy expectation with $\rho\exp(i\Phi_{T})$ that must be an upper bound to the exact energy for an arbitrary symmetric $\rho\geq 0$. The accuracy of this method clearly depends on the accuracy of the trial phase and the convergence towards the exact eigenvalue scales with the square of the difference between the exact and approximate trial function.

II.3 Fixed-phase as a special case of the fixed-node in general

The fixed-phase approximation is a generalization of the more familiar fixed-node approximation, what can be demonstrated in several ways. Let us present perhaps the simplest such construction Melton, Bennett, and Mitas (2016), where we add a complex amplitude to a real-valued $\Psi_{T}(\mathbf{R})$ as follows. We denote the nodes of $\Psi_{T}$ as the set of configurations

[TABLE]

Now we add to $\Psi_{T}$ another function (for simplicity, a nonnegative bosonic ground state of $H$)

[TABLE]

Taking the limit $\varepsilon\to 0$ leads toMelton, Bennett, and Mitas (2016)

[TABLE]

where $\mathbf{R}_{\Gamma}\in\Gamma$ and $V_{\infty}$ diverges as $\propto 1/\varepsilon^{2}$, therefore $V_{ph}$ enforces any wave function to vanish at the node $\Gamma$, i.e., it is equivalent to the fixed-node boundary condition.

II.4 Spin Representation

Let us denote one-particle spinors as

[TABLE]

where $s$ is the coordinate of the spin projection along the $z-$axis. In its minimal representation the spin variables have discrete values $s=\pm 1/2$ so that for $S_{z}$ eigenstates $\chi^{\uparrow}(1/2)=\chi^{\downarrow}(-1/2)=1$, $\chi^{\downarrow}(1/2)=\chi^{\uparrow}(-1/2)=0$. Clearly, the spin configuration space is non-compact and imposes potentially large variations of important quantities during the stochastic updates. Besides the fluctuations of various quantities of interest (local energy, drifts, values of the wave function, etc.) the method looses its efficiency in the many-particle limit.

One possibility to address this obstacle is to make the spin configuration space compact and continuous, which allows for continuous evolution as well as importance sampling Foulkes et al. (2001); Melton, Bennett, and Mitas (2016). We choose an overcomplete spin representation through the utilization of a 1D ring (i.e. a $U(1)$ representation) with the lowest pair of degenerate, orthogonal eigenstates as follows:

[TABLE]

where the spin variable $s_{j}\in[0,2\pi)$. Clearly, the paths in this space are continuous and resemble paths for spatial coordinates.

II.5 Importance sampling

Rewriting the Schrödinger equation in an integral form with importance sampling by $\rho_{T}$ leads to the following equation for the mixed distribution $g=\rho\rho_{T}$

[TABLE]

which is well-known from the fixed-node QMC Foulkes et al. (2001); Melton, Bennett, and Mitas (2016).

Spin variables are sampled by introducing a spin “kinetic” energy with a corresponding energy offset such that for all $s_{i}$, $i\in\{1,2,...,N\}$ we write

[TABLE]

where $\mu_{s}$ is an effective mass. The full Hamiltonian then becomes $H^{\prime}=H+\sum_{i=1}^{N}T_{i}^{s}$. Clearly, $T_{i}^{s}\psi(\mathbf{r}_{i},s_{i})=0$ due to the introduced offset so that there is no energy contribution from the spin Laplacian. The inclusion of the spin kinetic energy leads to the following importance sampled Green’s function

[TABLE]

with

[TABLE]

where ${\bf X}=({\bf r}_{1},{\bf r}_{2},...,{\bf r}_{N},s_{1},s_{2},...,s_{N})=({\bf R},{\bf S})$. Here we have introduced a spin time step $\tau_{s}=\tau/\mu_{s}$ as well as ${\bf v}_{D}^{\bf R}=\nabla_{\mathbf{R}}\ln\rho_{T}(\mathbf{X})$ and ${\bf v}_{D}^{\bf S}=\nabla_{\mathbf{S}}\ln\rho_{T}(\mathbf{X})$ which correspond to the spatial and spin drifts. The local energy is given by $E_{L}=\textrm{Re}[(H^{\prime}\Psi_{T})/\Psi_{T}^{*}]$ Foulkes et al. (2001); Melton et al. (2016).

Note that there are two possible limiting cases with regard to $\tau,\tau_{s}$, namely, $\tau_{s}>>\tau$ and vice versa. If $\tau_{s}$ is much larger than the spatial step the spin degrees of freedom are evolving much faster so that effectively the spins are integrated out for each spatial step. That guarantees to provide the fixed-phase limit and it is expected that it will lead to the largest bias. Indeed, this is what we have observed Melton, Bennett, and Mitas (2016).

The opposite limit corresponds to very slow spin evolution so that the spin configuration appears as an almost static external field for (relatively) much faster spatial evolution. In the results section we come to this point again and we show that in this mode the simulations will enable us to recover the fixed-node solutions.

III Trial Wave Functions

The FPSODMC trial functions are built from spinors $\chi(\mathbf{r},s)=\alpha\varphi^{\uparrow}({\bf r})\chi^{\uparrow}(s)+\beta\varphi^{\downarrow}({\bf r})\chi^{\downarrow}(s)$ where orbitals $\varphi^{\uparrow},\varphi_{\downarrow}$ are calculated in spinor-based DFT, HF/DF or correlated methods. The full configuration space for particles is $\mathbf{X}=\{(\mathbf{r}_{1},s_{1}),\ldots,(\mathbf{r}_{N},s_{N})\}\in\mathbb{R}^{3N}\times[0,2\pi)^{N}$ and we write the trial wave function as

[TABLE]

with $i,k=1,...,N$. The particle correlations are explicitly approximated by the Jastrow factor $U(\mathbf{R})$ that captures two-particle and, possibly, higher order correlations, as customary in QMC calculations Foulkes et al. (2001); Bajdich and Mitas (2009); Kolorenč and Mitas (2011); Melton, Bennett, and Mitas (2016).

III.1 From fixed-phase to fixed-nodes

In this work we are particularly focused on the limit of vanishing spin-orbit and how the single-reference spinor determinant simplifies to the product of spin-up and spin-down determinants, i.e., to the usual fixed-node form. Let us show that this is indeed what happens for our spin representation as briefly sketched earlier Melton, Bennett, and Mitas (2016). Note that our previous exposition of this aspect was not formulated precisely Melton, Bennett, and Mitas (2016), so that we clarify it in detail here. For the sake of consistency with the previous paper we consider $N$ occupied spinors that can be grouped as $N/2$ Kramer’s pairs (for simplicity assuming $N$ to be even). We can write the Kramer’s pair as

[TABLE]

where the $\Delta\varphi$ is the spin-orbit induced splitting of the spatial orbital $\varphi$. The block of the first four rows/columns from the corresponding Slater determinant reads as follows

[TABLE]

Now we assume that the spin-orbit splitting $\Delta\varphi\to 0$ and then

[TABLE]

by elementary rearrangements (adding, subtracting rows with the same $\varphi$). Explicitly, this gives

[TABLE]

This effectively complexifies the usual real wave function with the additional difference that the Slater matrix is of size $N\times N$. Clearly, for arbitrary spins variables this wave function is different from the usual spin-up and spin-down product although in what follows we will demonstrate how to recover such the fixed-node form using an appropriate choice for the spin coordinates. Let us now assume that $s_{i}=s_{1},s_{3},...$ will become the spin-up channel while $s_{i}=s_{2},s_{4},...$ will end up being the spin-down channel. In order to reach this spin-up/down partitioning explicitly we restrict $s_{1},s_{3},s_{5},...,=s$, and $s_{2},s_{4},...=s^{\prime}$ where $s,s^{\prime}$ are distinct. Then we can write the determinant

[TABLE]

and eliminating elements in each odd row

[TABLE]

where

[TABLE]

Furthermore, by reshuffling the first two rows and columns and factorizing out the spins from the determinant we get

[TABLE]

After reshuffling the rest of rows and columns, the single determinant of spinors factorizes into the product of two determinants of spin-up and spin-down block matrices. Generalization to odd $N$ with unpaired spinor(s) is straightforward. Therefore this decomposition strictly depends on the fact that all the spins have to acquire one of the two distinct values as expected when going from continuous to the fixed-label representation.

III.2 Wave functions with full space-spin symmetries

In our recent paper we have probed into the behavior of such wave functions for simple cases Melton and Mitas . It is useful to use an example such as the Li atom wave function to illustrate various wave function forms we consider here (assuming usual nucleus-electrons Hamiltonian without spin terms). The full symmetry exact wave function for the Li atom doublet is given byWhite and Stillinger (1970)

[TABLE]

where the function $F$ depends only on the spatial coordinates. The function $F$ is the exact, irreducible, spatial variables-only eigenstate for the three electrons in the doublet state. Indeed, it corresponds to the exact fixed-node solution sought after, say, in the FNDMC method.

The single-configuration trial wave function in the fixed-node framework would look like

[TABLE]

where the electrons 1 and 2 are assigned as spin-up while the electron 3 is spin-down. Clearly, this is just a projection onto the spin state $|\uparrow\rangle_{1}|\uparrow\rangle_{2}|\downarrow\rangle_{3}$ with the single-reference term approximating the spatial part $F(1,2,3)$.

Our wave function with variable spins is given by

[TABLE]

It therefore results in determinantal approximations to the function $F$ with the phase factors from varying spins as coefficients. If one chooses, $s_{1}=s_{2}=s$ and $s_{3}=s^{\prime}$, the wave function collapses to a single determinant with a spin variable dependent coefficient as described above.

Sampling both spin and position spaces enables one to evolve between the spatial wave functions with permuted coordinates, i.e., eventually sampling all such equivalent possibilities. Note that the overall structure of the exact wave function and our variable spin formulation are analogous. While in this example the variable spins and corresponding phase factors appear superfluous, the form becomes fully meaningful whenever spin-dependent terms in the Hamiltonian are switched on.

IV Fixed-phase variable spins QMC as a general method

In fixed-node QMC calculations with real wave functions the node improvement is often very challenging since any general method proposed so far appears to have very unfavorable scaling. In several papers we have made some progress in understanding the relations between electron density, multiplicity of bonds and node curvatures that appear to be related with increased fixed-node bias Rasch and Mitas (2012). In addition, we found relationships between nodes and eigenvalues that show the nodes carry information about the spectrum as presented elsewhere Melton, Bennett, and Mitas .

In this respect, the fixed-phase approximation opens new perspectives both in a better understanding of related issues with regard to antisymmetry and the corresponding fermion sign problem as well as possibilities for new constructions of more efficient approximations.

One important property of the fixed-phase approximation is that the sampled distribution $\rho$ is non-negative everywhere and, as we mentioned, generically its zero locus is a subset of configurations with codimension 2, i.e., two dimensions lower than the full configuration space. In that case the sampling is of the configuration space is ergodic. One then solves for the bosonic ground state in a given, state-dependent potential. A simple toy example is an atomic two-particle ${}^{3}P$ state with the wave function

[TABLE]

where $g$ are positive radial functions. Its phase-generated potential is given by

[TABLE]

while the corresponding non-negative amplitude $\rho(r_{1},r_{2})$ vanishes only at $x_{1}=x_{2},y_{1}=y_{2}$ Melton and Mitas . This has also other consequences that make it favorable in comparison with the fixed-node approach, namely, the divergences of the local energy and drift are significantly diminished making them much smoother. For example, the drift for the importance sampled distribution given by $\nabla\ln\rho$ is smooth except at the point of vanishing $\rho$. This removes complications around nodes of real functions such as large local energy fluctuations, non-zero probability of crossing/re-crossing the node within a given time step, possible occurrences of stuck walkers and others, due to the fact that $\ln\Psi_{T}$ is non-analytic at the node. All these difficulties can be brought under the control by decreasing the time step in the fixed-node formalism. However, here these complications are simply absent in the fixed-phase formulation by being smoothed out into the lower dimension.

We note that in low-dimensional systems or for particular symmetry constraints one can end up with special or non-generic cases having zeros of $\rho$ with codimension 1, i.e., one dimension higher than the generic codimension 2 mentioned above. A simple example is the lowest two-particle triplet in a periodic box with the wave function ${\rm det}[1,e^{ikx}]$. This leads to $\rho(x_{1},x_{2})=2|\sin[(x_{1}-x_{2})/2]|$ that has a $(2d-1)$-dimensional zero locus regardless of $d$, i.e., the dimensionality of the box. The reason is that this particular state effectively behaves as having 1D nodal structure that is non-generic. Interestingly enough, for $d>1$ this node volume is smaller than in the corresponding fixed-node wave function given by the real (or imaginary) part, ${\rm Re}\{{\rm det}[1,e^{ikx}]\}$. This aspect is more thoroughly investigated in our subsequent work Melton, Bennett, and Mitas that explore the corresponding properties of nodes in such cases and further generalizations.

Perhaps the most appealing and yet unexplored property is that the approximation has a form of an additive effective many-body potential

[TABLE]

so that the original Schrödinger equation changes to

[TABLE]

This effective potential formulation offers a clear conceptual understanding of the transformed problem that reminds us of effective potential/field methods used in other areas of quantum and high energy physics. It has a number of desired properties when thinking about the solution of the many-body problem, such as that the solution is non-negative everywhere, the state-dependent potential $V_{ph}$ is purely repulsive (it only raises the energy) and it is explicitly and directly given by the approximate phase. Consequently, it provides a constructive path for improvements with the perspective that the solution really exists, i.e., in the case of the exact phase one obtains the exact solution/eigenstate similarly to the fixed-node approximation (that is its special case). Interestingly, not much is known about the phases of stationary states. It is possible that more thorough analysis of the corresponding effective potentials will lead to a better understanding as well as to better approximations for practical calculations of realistic systems. What follows provides the first attempts to probe some aspects of this formulation.

V Results

We calculate the total energies for the first-row atoms using both the FN and FP approximations. For the FN calculations, we build our trial wave function from HF orbitals generated from Gamess-US Schmidt et al. (1993). For the FP calculations, we build our trial wave function from the one-particle spinors generated from Dirac14 Saue et al. (2014). For the FN calculations, we perform a linear time step extrapolation to zero time step. In all cases, a spatial time step of 0.001 Ha is in agreement with the zero time step limit. Motivated by that, for all FP calculations we hold the spatial time step fixed at $\tau=0.001$ Ha, rather than performing a spatial time step extrapolation.

We have previously studied some of the aspects of spin time steps Melton et al. (2016); Melton, Bennett, and Mitas (2016). The analysis of total energies as function of spin time step leads to the following conclusions:

a) At large spin time step the spins are basically fully integrated out for each spatial step that is assumed to be much smaller. Then one sees higher fixed-phase bias since the repulsive potential acts in the full configuration space unlike the fixed-node condition that applies only on the configuration subspace.

b) At very small spin steps and for small number of electrons, the propagation eventually finds the region(s) close to the pure fixed-node wave function. Apart from small spin fluctuations the energy therefore reaches very closely to the fixed-node solution.

In both limits and also for intermediate time step regimes the energy is an upper bound, since the energy basically limited from below by the fixed-node limit. The complexified wave function and the fixed-phase only increases the energy since it acts in full space instead of fixed-node codimension 1 hypersurface and expands the configuration space in an ad hoc manner through the continuous spin as we argued in previous parts. This has also further implications that single reference wave functions will be, in general, less accurate, as we have actually observed in calculations of several systems Melton, Bennett, and Mitas (2016); Melton and Mitas . Here we are actually focused on the short spin time-step limit that enables us to recover the fixed-node results although the calculations are carried out in FPSODMC setting.

As described in §III, we initialize the spin-configurations to facilitate the decomposition into two independent determinants, as must be the case in a spin independent Hamiltonian. As an illustration of why this is necessary, consider the N atom using HF spatial orbitals and no Jastrow factor. The VMC energy should agree with the HF energy, within the statistical errors. If one randomly generates the spin variables and performs a VMC calculation, the obtained energy is $-54.3341(8)$ Ha which clearly disagrees with the HF value of -54.40093 Ha. However, if we initialize the spin variables such that $s_{1}=s_{3}=s_{5}=s_{6}=s_{7}=s$ and $s_{2}=s_{4}=s^{\prime}$ while using a very small spin time step $\tau_{s}$, we obtain $-54.4003(7)$ Ha, which agrees with the energy obtained via HF. By initializing the spin variables such that the wave function properly decomposes into a product of determinants, we can use a small $\tau_{s}$ such that the spin variables stay close to the original configurations. We allow the spin variables to continue to drift in order to sample the spin configuration space.

When performing FPSODMC, we vary the spin time step until the energy is saturated for a fixed spatial time step. An example of the spin time step extrapolation is shown in Figure 1. For $\tau_{s}$ between $10^{-12}$ and $10^{-9}$, the DMC energies all agree to within the error bars. Performing the same procedure for all atoms, we list the total energies in Table 1. By comparing the FN and FP total energies to the “exact” energies in the non-relativistic limit (NRL)Davidson et al. (1991), we calculate the fixed-node/phase error as the percentage of the total as it is plotted in Figure 2.

Regardless of the approximation, the associated error decreases with atomic number subject to the choice of HF nodes/phases. Additionally, the FN and FP approximations yield essentially identical errors.

It is well known that using HF nodes, Be has a significant FN error Rasch and Mitas (2012). The ground state symmetry of Be is ${}^{1}S_{0}$, which is obtained with the electron configuration $1s^{2}2s^{2}$. Using a HF trial wave function, the nodal surface $\partial\Omega=\left\{\mathbf{R}\in\mathbb{R}^{3N}|\Psi_{T}(\mathbf{R})=0\right\}$ separates the configuration space into 4 nodal domains, two of which where the wave function is positive and two in which it is negative. It is well-known that by adding just one more configuration that is related to the near-degenerate state of the same symmetry one finds only two nodal domains as expected for generic fermionic ground state Ceperley (1991); Bressanini, Ceperley, and Reynolds (2002); Mitas (2006). Previous calculations have found that it almost completely eliminates the fixed-node bias Umrigar, Wilson, and Wilkins (1988); Bressanini, Ceperley, and Reynolds (2002). The corresponding two-configuration trial function is given by

[TABLE]

as With this choice of trial wave function and full optimization all variational parameters one can reach FN result with almost zero bias Umrigar, Wilson, and Wilkins (1988). Instead we perform an optimization of this wave function with only the Jastrow parameters and expansion coefficient, keeping the HF orbitals fixed with resulting small increase in the energy compared to the nearly exact value. Total energies are shown in Figure 3. Again, by choosing the FP calculations to preserve the spin assignments of $s_{1}=s_{3}=s$ and $s_{2}=s_{4}=s^{\prime}$ by using a small spin time step, the FN and FP calculations agree to within statistical uncertainty.

Thus far, we have only presented results for all-electron systems. We also consider the FN and FP approximation when nonlocal pseudopotentials are included Melton et al. (2016); Melton, Bennett, and Mitas (2016). We calculate the binding curve for the nitrogen dimer in the ${}^{1}\Sigma_{g}$ state, using a single-reference trial wave function for each approximation built from HF spatial orbitals. We utilize a BFD pseudopotential for N Burkatzki, Filippi, and Dolg (2007). In order to calculate the binding curve, we first calculate the isolated N atom both in FN and FP. Under the locality approximation Mitáš, Shirley, and Ceperley (1991), we perform a time step extrapolation within FNDMC and obtain a total energy of $-9.7912(1)$ Ha. Using a fixed spatial time step, we perform a spin time step extrapolation as described above to facilitate decomposition into a product of two independent determinants using FPSODMC and obtained a total energy of $-9.7917(4)$ Ha.

The dimer curve is shown in Figure 4 and shows the binding obtained from the FN and FP methods. The QMC data is fit to the morse potential

[TABLE]

and the vibrational frequency can be obtained via

[TABLE]

where $\mu$ is the reduced mass of the dimer. The FP solution for the dimer has a slightly larger bias in comparision to the FN solution, on the order of 1 mHa across the entire binding curve. Coupled with the slightly lower energy for the individual atom, the overall binding energy differs from the FN result by roughly $\sim$ 0.1 eV, as shown in Table 2. For completeness, we calculated the dimer using an improved nodal surface/phase given by a trial wave function composed of PBE0 nodes at the equilibirum bond length. At $r_{e}=1.09$ Å, the PBE0 nodal surface is lower in energy by only 0.0010(5) Ha, which slightly improves the binding energy prediction to $9.654(8)$ eV. The PBE0 phase has a more significant improvement over the HF phase, yielding a lower energy by 0.0029(5) Ha. The binding energy prediction becomes $9.64(1)$ eV, which is very close to the FN result. Clearly, the differences between the methods are very small, basically similar to variations in the fixed-node biases for different atoms and molecular systems and choices of orbitals used in single-reference trial functions.

VI Conclusions

In this paper we elaborate in detail on a particularly important aspect of the fixed-phase spin-orbit/spinors DMC (FPSODMC) method that we have introduced recently Melton et al. (2016). We highlight some of the key aspects, in particular, how to obtain the fixed-node limit results from the fixed-phase setting both in theory and in practical calculations. We point out the promising features of the fixed-phase method and also show its behavior in our continuous spin formalism. We illustrate the results on first row atoms calculations. The method enables us to write full space-spin symmetry wave functions for Hamiltonians with or without explicit spin terms and opens thus possibilities for further improvements of trial wave functions. We consider the results very encouraging since in a straightforward manner we were able to obtain the fixed-node results in both all-electron and effective core potential settings as well as confirm essentially the same quality of both single and multi-reference trial wave functions. The method opens interesting new perspectives for many-body electronic structure calculations in complex wave function and spinor formalism that takes into account variable nature of the spin degrees of freedom and provides new possibilities for construction of more general trial wave functions.

Acknowledgments. This research was supported by the U.S. Department of Energy (DOE), Office of Science, Basic Energy Sciences (BES) under Award DE-SC0012314. For calculations we used resources at NERSC, a DOE Office of Science User Facility supported by the Office of Science of the U.S. Department of Energy under Contract No. DE-AC02-05CH11231. Most calculations have been carried out at TACC.