Local embedding and effective downfolding in the auxiliary-field quantum Monte Carlo method

TL;DR

This paper introduces a local embedding and downfolding scheme in AFQMC that reduces computational costs by focusing on a local cluster and effectively downfolding virtual orbitals, enabling larger system simulations.

Contribution

The methodology extends local embedding in AFQMC to include virtual orbital downfolding, significantly reducing computational effort while maintaining accuracy.

Findings

AFQMC energies converge rapidly with respect to active region sizes.

The approach allows treatment of larger systems with reduced computational resources.

Systematic dependence on active region sizes is characterized.

Abstract

A local embedding and effective downfolding scheme has been developed and implemented in the auxiliary-field quantum Monte Carlo (AFQMC) method. A local cluster in which electrons are fully correlated is defined and the frozen orbital method is used on the remainder of the system to construct an effective Hamiltonian which operates within the local cluster. Local embedding, which involves only the occupied sector, has previously been employed in the context of Co/graphene. Here, the methodology is extended in order to allow for effective downfolding of the virtual sector thus allowing for significant reduction in the computational effort required for AFQMC calculations. The system size which can be feasibly treated with AFQMC is therefore greatly extended as only a single local cluster is explicitly correlated at the AFQMC level of theory. The approximation is controlled by the separate choice of the spatial size of the active occupied region, $R_o$, and of the active virtual region, $R_v$. The systematic dependence of the AFQMC energy on $R_o$ and $R_v$ is investigated and it is found that relative AFQMC energies of physical and chemical interest converge rapidly to the full AFQMC treatment (i.e. using no embedding or downfolding).

Equations

$$\hat{H}=\hat{K}+\hat{V}=\sum _{\mu \nu }{K}_{\mu \nu }{c}_{\mu }^{†}{c}_{\nu }+\frac{1}{2}\sum _{\mu \nu \lambda \rho }{V}_{\mu \nu \lambda \rho }{c}_{\mu }^{†}{c}_{\nu }^{†}{c}_{\lambda }{c}_{\rho },$$

$$\begin{array}{r}\underset{\beta \to \infty }{\lim }{e}^{-\beta \hat{H}}|{\Psi }_T⟩\approx e^{-\tau \hat{H}}{e}^{-\tau \hat{H}}\cdots e^{-\tau \hat{H}}|{\Psi }_T⟩\to |{\Psi }_0⟩,\end{array}$$

$$\begin{array}{r}{E}_0=\frac{⟨{\Psi }_T|\hat{H}|{\Psi }_0⟩}{⟨{\Psi }_T|{\Psi }_0⟩}=\underset{\beta \to \infty }{\lim }\frac{⟨{\Psi }_T|\hat{H}{e}^{-\beta \hat{H}}|{\Psi }_T⟩}{⟨{\Psi }_T|e^{-\beta \hat{H}}|{\Psi }_T⟩}\end{array}$$

$$\begin{array}{r}{e}^{-\tau \hat{H}}\approx e^{-\tau \hat{K}/2}{e}^{-\tau \hat{V}}{e}^{-\tau \hat{K}/2}+\mathcal{O}\left({\tau }^3\right).\end{array}$$

$$\begin{array}{r}{e}^{-\tau \hat{H}}=\int d\mathbf{\sigma }P(\mathbf{\sigma })e^{-\tau \hat{K}/2}{e}^{\sqrt{\tau }\mathbf{\sigma }\cdot \hat{\mathbf{v}}}{e}^{-\tau \hat{K}/2},\end{array}$$

$$\begin{array}{r}{e}^{-\tau \hat{K}/2}{e}^{\sqrt{\tau }\mathbf{\sigma }\cdot \hat{\mathbf{v}}}{e}^{-\tau \hat{K}/2}|\varphi ⟩\equiv e^{-\tau \hat{h}(\mathbf{\sigma })}|\varphi ⟩\to |{\varphi }^{\prime }⟩.\end{array}$$

$$\begin{array}{r}{e}^{-\tau \hat{K}/2}{e}^{\sqrt{\tau }\left(\mathbf{\sigma }-\bar{\mathbf{\sigma }}[\varphi ]\right)\cdot \hat{\mathbf{v}}}{e}^{-\tau \hat{K}/2}|\varphi ⟩\to |{\varphi }^{\prime }⟩,\end{array}$$

$$\begin{array}{r}\bar{\mathbf{\sigma }}[\varphi ]\equiv -\sqrt{\tau }\frac{⟨{\Psi }_T|\hat{\mathbf{v}}|\varphi ⟩}{⟨{\Psi }_T|\varphi ⟩}.\end{array}$$

$$\begin{array}{r}{E}_0\approx \frac{{\sum }_{\varphi }{w}_{\varphi }{E}_{\mathrm{L}}[\varphi ]}{{\sum }_{\varphi }{w}_{\varphi }},\end{array}$$

$$\begin{array}{r}{E}_{\mathrm{L}}[\varphi ]\equiv \frac{⟨{\Psi }_T|\hat{H}|\varphi ⟩}{⟨{\Psi }_T|\varphi ⟩}\end{array}$$

$$\Psi \approx \mathcal{A}({\Psi }_{\mathbb{I}}{\Psi }_{\mathbb{A}}),$$

$$E=⟨\Psi |\hat{H}|\Psi ⟩=⟨{\Psi }_{\mathbb{A}}|{\hat{H}}_{\mathbb{A}}|{\Psi }_{\mathbb{A}}⟩,$$

$$\begin{gathered} {\hat{H}}_{\mathbb{A}}=\sum _{ij\in \mathbb{A}}{K}_{ij}{c}_i^{†}{c}_j+\frac{1}{2}\sum _{ijkl\in \mathbb{A}}{V}_{ijkl}{c}_i^{†}{c}_j^{†}{c}_k{c}_l \\ +\sum _{ij\in \mathbb{A}}{V}_{ij}^{\mathbb{I}-\mathbb{A}}{c}_i^{†}{c}_j+E^{\mathbb{I}}. \end{gathered}$$

$${\mu }_p^i=⟨{\varphi }_i|(\hat{r}-⟨{\varphi }_i|\hat{r}|{\varphi }_i⟩{)}^p|{\varphi }_i⟩;\left({\sigma }_p^i=\sqrt[p]{{\mu }_p^i}\right),$$

$${\xi }_p^m=\sum _i{\left(u_p^i\right)}^m.$$

Full text

I Introduction

The solution of the quantum many-electron problem is of fundamental importance to the prediction of the properties of both molecular and condensed matter systems. It is one of the outstanding problems and often a bottleneck in the simulation of materials. Due to electron-electron interactions, the computational cost of exact calculations scales exponentially with system size; exact calculations are possible for only the smallest systems. Approximate many-particle methods, such as coupled cluster (CC) methods in quantum chemistry, scale as high-order polynomials in the number of electrons $N$ ( $N^{7}-N^{8}$), which limits their scope. Unlike wave-function-based CC and other quantum chemistry approaches, Quantum Monte Carlo (QMC) techniques use stochastic methods to directly compute observables. The scaling with system size is reduced to $N^{3}-N^{4}$, as in density functional theory or Hartree-Fock methods, making applications to larger systems feasible. Moreover, QMC has demonstrated high accuracy in systems with strong correlated-electron effects.

Many systems of interest include a spatially localized cluster in which electron correlation effects are expected to be more important than in the rest of the system. It is possible to treat only the local cluster with a highly accurate, but expensive, computational method while treating the rest of the system with a less expensive, but less accurate, method. The effective system size at the many-body level of theory is therefore smaller than that of the original system allowing for an overall larger system to be treated. Local embedding methods were recently reviewed in great detail in the context of a variety of computational methods and embedding strategies Gordon (2017). The auxiliary-field QMC (AFQMC) method Zhang and Krakauer (2003); Al-Saidi et al. (2006); Motta and Zhang (2017), which scales favorably with system size, as compared with methods capable of a similar degree of accuracy, is an ideal tool for embedded calculations. AFQMC can use any set of one-particle orbitals. This can be exploited to use unitarily localized orbitals which can be spatially assigned to either the embedding or the host region.

It was previously shown that local embedding could be employed to reduce the size of the occupied sector in Co/coronene (C24H12) Virgus et al. (2014), leading to highly accurate results directly comparable to experiment. However, with local embedding alone, the entire virtual space is included at the AFQMC level of theory. The virtual space becomes extremely large with high-quality basis sets, especially as the complete basis set limit is approached, which still limits the size of systems which can be treated using local embedding alone. It is therefore desirable to employ an effective local downfolding scheme in order to reduce the size of the virtual sector as well.

In this paper, we show that local embedding, in the occupied space, and local effective downfolding, in the virtual space, can be employed in combination in order to treat strongly correlated local clusters at the AFQMC level of theory at reduced computational cost. Furthermore, the accuracy of embedded and downfolded results can be controlled by the choice of embedding and downfolding parameters. Significant computational savings are achieved in a number of applications while maintaining a high degree of accuracy.

The remainder of the paper is organized as follows. Sec. II reviews the AFQMC method and presents a procedure for achieving local embedding combined with local effective downfolding, and concludes with computational details. A simple application which illustrates the local embedding and downfolding approximation is provided in Sec. III. Additionally, the dependence of the systematic errors on the spatial size of the embedding region is studied. In Sec. IV, local embedding and effective downfolding is applied to Ti-capped carbon chain systems within a graphitic environment. First, the systematic convergence of the AFQMC energy in the size of the embedding region is studied. Next, basis-set and finite-size effects in the Ti-capped carbon chain system are considered. Finally, it is demonstrated that similar systematics are observed for Ti-capped carbon chain systems which interact with a graphitic environment. Additional issues and practical limitations of local embedding and downfolding are discussed in Sec. V. We then conclude with some general remarks in Sec. VI.

II Methodology

In this section, we first review key points of the AFQMC method. Zhang and Krakauer (2003); Al-Saidi et al. (2006); Motta and Zhang (2017) We then describe our embedding/downfolding approach, followed by computational details used in the calculations.

II.1 AFQMC Overview

The Hamiltonian in second-quantized form is given by,

[TABLE]

where ${\hat{K}}$ and ${\hat{V}}$ include all one-body and two-body terms, respectively. The indices run over $M$ orthonormal single-particle basis functions, with $c^{\dagger}_{\mu}$ and $c_{\mu}$ the respective creation and annihilation operators. This form of the Hamiltonian applies to all approaches for interacting electron systems, from realistic orbital-based approaches for materials and molecular systems to lattice-based models.

The ground-state projection AFQMC formalism uses a trial wave function $|{{\Psi_{T}}}\rangle$ from a less-expensive method, such as HF or DFT. The resulting canonical orbitals could be used as the orthonormal basis functions in Eq. (1), but for embedding/downfolding, the canonical orbitals are first unitarily transformed to localized basis functions, as described in the next subsection. Iterative imaginary-time projection is used to obtain the ground state $|{{\Psi_{0}}}\rangle$ from a trial wave function $|{{\Psi_{T}}}\rangle$:

[TABLE]

where $\beta$ is the total imaginary projection time, and $\tau$ is a small imaginary time step. The projection is cast as a branching random walk with Slater determinants, $|{{\phi}}\rangle$. The ground state energy, $E_{0}$, is computed using the mixed estimator,

[TABLE]

The small imaginary time step allows a Trotter-Suzuki decomposition of the imaginary time electron propagator,

[TABLE]

The exponential of a one-body operator acting on a Slater determinant simply produces another Slater determinant. $e^{-\tau{\hat{V}}}$ can be cast as a high-dimensional intergral over many auxiliary fields, $\bm{\sigma}$, by the application of a Hubbard-Stratonovich transformation Zhang et al. (1997); Zhang and Krakauer (2003) as given by

[TABLE]

where $P(\bm{\sigma})$ is a normal distribution function and $\hat{\bm{v}}$ is a one-body operator. The operation of the integrand on a Slater determinant $|{{\phi}}\rangle$ is given by,

[TABLE]

where $|{{\phi^{\prime}}}\rangle$ is another Slater determinant. An initial population of walkers is prepared (usually set equal to $|{{\Psi_{T}}}\rangle$) and the mixed estimator, given in equation 3, is stochasically sampled. In general, the one-body operator $\hat{h}(\bm{\sigma})$ is complex and as the projection proceeds a complex phase is accumulated for each walker $|{{\phi}}\rangle$. This causes statistical fluctuations to increase exponentially with projection time. The phaseless approximation was introduced Zhang and Krakauer (2003) in order to control this problem. First, an importance function, given by $\langle{{\Psi_{T}}}|{{\phi}}\rangle$, is introduced and an importance sampling transformation is applied. Eq. 6 then becomes,

[TABLE]

where the force bias, $\bar{\bm{\sigma}}[\phi]$ is given by

[TABLE]

The mixed estimator at each time step is then given by a weighted sum over random walkers

[TABLE]

expressed in terms of the local energy,

[TABLE]

of each walker and the weight, $w_{\phi}$, is accumulated over the random walk. The procedure is described in more detail in Refs. [2; 3; 7; 8; 4].

II.2 Embedding and Downfolding Approach

Many systems of interest include a spatially localized cluster where electron-correlation effects are expected to be more important than in the rest of the system. In this section, we describe the partitioning of the Hilbert space into active $\mathbb{A}$ and inactive $\mathbb{I}$ sectors, which are used to define an effective Hamiltonian, ${\hat{H}}_{\mathbb{A}}$, to be used in high-accuracy AFQMC calculations. The construction of ${\hat{H}}_{\mathbb{A}}$ is schematically represented in Fig. 1. Since ${\hat{H}}_{\mathbb{A}}$ includes only one-body embedded contributions from the $\mathbb{I}$ sector, as described below, this greatly extends the system size which can be accurately treated by AFQMC.

The embedding/downfolding effective Hamiltonian, ${\hat{H}}_{\mathbb{A}}$, operates on a reduced Hilbert space, spanned by the unitarily localized occupied and virtual $\mathbb{A}$ orbitals, depicted in the third column of Fig. 1. ${\hat{H}}_{\mathbb{A}}$ is obtained using a separability approximation for the many-body wave function: Kahn et al. (1976); Huzinaga (1991)

[TABLE]

where $\Psi_{\mathbb{I}}$ and $\Psi_{\mathbb{A}}$ depend, respectively, only on the $\mathbb{I}$ and $\mathbb{A}$ orbitals in Fig. 1. Both $\Psi_{\mathbb{I}}$ and $\Psi_{\mathbb{A}}$ are separately antisymmetric and normalized, and mutually orthogonal; the antisymmetrizer $\mathcal{A}$ permutes electrons between $\Psi_{\mathbb{I}}$ and $\Psi_{\mathbb{A}}$. This approximation allows the energy of the total system to be mapped onto an equivalent system which explicitly includes only the $\mathbb{A}$ orbitals:

[TABLE]

where ${\hat{H}}_{\mathbb{A}}$ is given by Purwanto et al. (2013):

[TABLE]

The first two terms in Eq. 13 are one and two-body terms that depend only on the $\mathbb{A}$ orbitals. The third term is a one-body embedding potential that describes Coulomb and exchange interactions between the $\mathbb{A}$ orbitals and the environment, which is spanned by the $\mathbb{I}$ orbitals. Formally, the third term is a nonlocal pseudopotential. The constant-energy fourth term $E^{\mathbb{I}}$ includes the kinetic and internal interactions of orbitals in $\mathbb{I}$. The effective Hamiltonian ${\hat{H}}_{\mathbb{A}}$ has the same form as Eq. (1), so AFQMC calculations proceed exactly as outlined in Sec. II.1.

A unitary transformation of the occupied and virtual orbitals to localized orbitals is accomplished by minimizing a cost function that depends on the localized orbitals. The most used methods are Foster-Boys (FB) Boys (1960, 1966), Pipek-Mezey (PM) Pipek and Mezey (1989) and Edmiston-Ruedenberg (ER) Edmiston and Ruedenberg (1963). Some methods, such as FB, use the $p^{th}$ moments $\mu^{i}_{p}$, (and its spread $\sigma^{i}_{p}$) Høyvik et al. (2012):

[TABLE]

with the cost function specified by:

[TABLE]

The FB method, for example, uses $\xi^{1}_{2}$. Compared to the occupied manifold, FB, PM and ER have more difficulty optimizing the virtual orbitals. In recent work, the use of $\xi^{1}_{4}$ based on the fourth central moments (FM) has been shown to be effective for localizing virtual orbitals Høyvik et al. (2012). FM localization can be further improved using $\xi^{2}_{4}$ (FM2), with a penalty-exponent $m=2$ (as previously done for the FB cost function Jansík et al. (2011)). FM2 was found to satisfactorily achieve a uniformly localized set of virtual orbitals. Høyvik et al. (2012). FM and FM2 place more weight on the orbital tail and therefore tend to suppress orbitals with long tails.

A more stringent criterion for a set of localized orbitals that has been considered Høyvik et al. (2012) is the spread of the least-local orbital. This can be useful in divide-and-conquer local-correlation methods that divide the system into local clusters.Knizia and Chan (2013); Yoshikawa and Nakai (2017) Here we focus on systems with a single local cluster, so the efficiency of the method is affected only by the localization of virtual orbitals in the cluster neighborhood and not by those far from the cluster.

II.3 Computational details

In this paper, the canonical occupied and virtual orbitals are obtained from restricted Hartree-Fock (RHF) calculations using NWChem Valiev et al. (2010). The RHF wave functions are also used as the trial wave function in all AFQMC calculations. Standard gaussian-type orbital (GTO) cc-pVXZ basis setsSchuchardt et al. (2007) are used, as described in the sections below. AFQMC calculations were done using the generic second-quantized Hamiltonian in Eq. (1) for the full system and Eq. (13) for the active cluster, as described by Al-Saidi et al. Al-Saidi et al. (2006).

The ERKALELehtola et al. (2012) computer code is used to apply FB localization to the occupied orbitals and FM localization to the virtual orbitals, unless otherwise stated. Although ERKALE has not implemented, to date, trust-region minimization Høyvik et al. (2012); Lehtola and Jónsson (2013), which is more robust than gradient-based minimization algorithms (particularly for FM and FM2 cost functions), its localization engine is capable of producing orbitals which are sufficiently local for the purpose of the present work, as indicated in the sections below.

III Illustrative Application

The methodology of the local embedding and effective downfolding method is illustrated in this section for a simple model system consisting of a single oxygen atom and a linear chain of 20 hydrogen atoms as depicted in Fig. 2. The atomic coordinates are given in Table LABEL:sup-tab:O+H20-geom. The AFQMC total energy is calculated as a function of $h$, the position of the O atom above the H-H bridge site at one end of the linear chain. RHF calculations were performed using a cc-pVDZ basis for both O and H atoms. The Hamiltonian is represented in the basis of localized RHF orbitals, using FB and FM localization on the occupied and virtual orbitals respectively. FM2 localization failed during the line search procedure as implemented in ERKALE Lehtola and Jónsson (2013) for the virtual orbitals of O/H20 with default settings.

We regard the active region as roughly centered on the O atom and including some neighboring H atoms and define the localization radii $R_{o}$ and $R_{v}$ from the position of the O atom, projected onto the chain, to the centroid of the orbital. Fig. 3 shows the numbers of occupied and virtual orbitals assigned to $\mathbb{A}$, $N_{o}$ and $N_{v}$, as a function of $R_{o}$ and $R_{v}$, respectively, for two values of $h$. The plotted distribution is based on the centroids of the localized orbitals. Beyond about 1-2 Bohr, the distributions for the two heights are nearly identical. As shown below, this leads to relative insensitivity of total energy differences when $R_{o}$ and $R_{v}$ are varied. The unit-step feature for the occupied localized orbitals reflects that they are centered on the H-H bond centers away from the O atom.

Figure 4 displays the calculated AFQMC energy of the O/H20 chain as a function of $R_{v}$ for four values of $R_{o}$ and two values of $h$. All the curves show similar qualitative convergence behavior. For $R_{v}<R_{o}+C$, where $C\approx 3$ Bohr, the total energy decreases approximately linearly with increasing $R_{v}$. For $R_{v}>R_{o}+C$ the total energy is converged very nearly to the full virtual space treatment (i.e. all virtual orbitals are included in $\mathbb{A}$) for the chosen $R_{o}$, which is indicated at the $R_{v}$ value corresponding to the full length of the H chain in the figure. In general, $C$ is some system-dependent constant as will be shown in the next Section. A key point is that, for the same $R_{o}$ but different $h$, the total energy curves are nearly parallel across almost the entire range of $R_{v}$, even before the plateau region; as $R_{o}$ is increased, the energy difference between pairs of curves with different $h$ quickly converges. Thus, relative energies such as binding energies or potential energy curves converge rapidly with increasing localization radii, as we demonstrate next.

The computed potential energy curves (PECs) for the O/H20 chain are shown as a function of the O-atom height, $h$, in Fig. 5. Embedding/downfolding PECs are shown for three choices of the cutoff radii and compared to exact AFQMC results, (ı.e., AFQMC for the entire O/H20 chain without downfolding/embedding approximations). The choice of the ($R_{o}$, $R_{v}$) cutoffs was guided by the cutoff-radii convergence properties shown in Fig. 4. Quantitatively good agreement is found across all values of $h$ for these cutoffs. The smallest cutoffs shown, $(R_{o},R_{v})=(3.7,6.2)$ Bohr, correspond to an active region consisting of the O atom and the 4 nearest H atoms. The corresponding effective Hamiltonian operates on a greatly reduced Hilbert space, yielding significant computational savings. For ($R_{o}$, $R_{v}$)$=$($7.2$, $9.8$) Bohr, the dimension of the single particle Hilbert space, $M$, is reduced from 113 for the full system (with atomic core frozen), to only 43 and the number of active electrons, $N$, is reduced from 13 to only 6 per spin sector. Based on the theoretical scaling of AFQMC ($\propto M^{2}N^{2}$), this reduction in the size of the active Hilbert space corresponds to a reduction in computational cost by a factor, $f\approx 27$. For ($R_{o}$, $R_{v}$)$=$($3.8$, $8.3$) Bohr, the effective $M$ = $42$ and $N$ = $5$ giving $f\approx 50$. ($R_{o}$, $R_{v}$)$=$($3.8$, $6.2$) Bohr, has $M$ = $37$ and $N$ = $5$ giving $f\approx 63$. The reduction in computational cost for O/H20, even for the largest choice of cutoff radii, is significant in its own right but is modest when compared to the results of the next section.

IV Applications

In this section our embedding method is applied to carbon chain systems in a graphitic environment. Our choice is motivated by the recent successful fabrication of long linear carbon chains (LLCCs) confined within double walled carbon nanotubes (DWCNTs),Shi et al. (2016) and the possible technological and scientific applications of TM-functionalized LLCCs.

IV.1 1-Dimensional correlated system : Ti capped linear carbon chain

We first apply our method to a H-Ti-C29-H linear chain. The initial geometry was obtained as follows. First, an H-C30-H chain was fully relaxed with DFT PW91 calculations using NWChem. One H-C end-unit was replaced by an H-Ti unit and relaxed, keeping fixed the geometry of the remaining C29-H unit. The Ti-C bond is a triple bond. The coordinates of atomic centers are listed in Table LABEL:sup-tab:H-Ti-C29-H-geom. We used cc-pVDZ basis sets for Ti and C, and the 6-31g basis set for the terminating H atoms. Due to near linear dependence of the basis functions in some of the carbon chain systems, we slightly modified the C cc-pVDZ basis by changing some of the p-function exponents on the C atoms (given in Fig. LABEL:sup-bas:custom-cc-pVDZ-basis). The basis sets are given in Fig. LABEL:sup-bas:custom-cc-pVDZ-basis (used for most Ti-C bond lengths) and in Fig. LABEL:sup-bas:custom-cc-pVDZ-basis-1.6 (used for the compressed bond length of 3.0 Bohr). Although accurate treatment of transition metal atoms in many-body calculations usually requires larger basis sets and multi-determinant trial wave functions Purwanto et al. (2016), we focus here on the systematics of the downfolding and embedding method, and restrict ourselves to the modest Ti cc-pVDZ basis and single-determinant trial wave functions. (We do, however, examine basis size-dependence below for a shorter H-Ti-C7-H linear chain.)

Figure 6 shows bar plots of the fourth central moment orbital spreads, $\sigma_{4}$, of the localized virtual (upper panel) and localized occupied (lower panel) orbitals. Since FM2 localization of the virtuals did not converge, FB localization was used for the occupied and FM localization was used for the virtual orbitals. Although many virtual states in Fig. 6 have rather large $\sigma_{4}$ values, most of these orbitals are far from the Ti atom. At the same time, occupied and virtual orbitals near Ti are seen to have much smaller $\sigma_{4}$ than most other orbitals. It is therefore expected that convergence with respect to the localization radii, $R_{o}$ and $R_{v}$, will only weakly depend on $\sigma_{4}$ of the least local orbital. This is confirmed below.

The H-Ti-C29-H PECs are shown in Fig. 7 for three sets of ($R_{o}$, $R_{v}$). For each $R_{o}$ in Fig. 7, Fig. 8 shows the convergence for fixed $R_{o}$ as a function of $R_{v}$, for two Ti-C bond lengths ($\delta_{\textrm{Ti-C}}=3.40$ and 4.35 Bohr). The energy in Fig. 8 is seen to quickly plateau as $R_{v}$ becomes greater than $R_{o}$. Thus, the $R_{o}$ cutoff effectively sets the size of the number of one-particle states that must be correlated in the effective Hamiltonian. The rapid convergence with increasing numbers of virtual orbitals is similar for curves with the same $R_{o}$ but different $\delta_{\textrm{Ti-C}}$. The curves for $\delta_{\textrm{Ti-C}}=3.40$ and 4.35 Bohr are nearly parallel across almost the entire range of $R_{v}$, reflecting good cancellation of errors for the Ti-C chain, similar to O/H20 (Fig. 4). This is shown by the inset of Fig. 8. The PECs in Fig. 7 are converged in $R_{o}$ at about $R_{o}$ = 9.7 Bohr for all bond lengths except $\delta_{\textrm{Ti-C}}=3.0$ Bohr.

IV.2 Finite-Size and Basis-Set Effects in the Linear Chain

In the next subsection, we present embedding/downfolding calculations for a H-Ti-Cn-H chain on a graphitic substrate. These were done for a shorter H-Ti-C7-H chain. Calculations for shorter chains facilitate direct comparisons with full AFQMC calculations, allowing us to focus on the embedding/downfolding approximations, while reducing the additional computational cost due to the graphitic substrates. Starting from the equilibrium structure H-Ti-C29-H chain, 22 interior C atoms were removed, keeping the remaining bond lengths the same. Atomic positions are given in Table LABEL:sup-tab:H-Ti-C7-H-geom. Unless otherwise stated, the custom C cc-pVDZ basis used to remove near linear dependence in H-Ti-C29-H (Fig. LABEL:sup-bas:custom-cc-pVDZ-basis) is used for H-Ti-C7-H as well. FB localization is used for both the occupied and virtual orbitals as FM did not converge for H-Ti-C7-H on the graphitic substrate studied in the next section.

We first calculated energies for H-Ti-C7-H as a function of the Ti-C bondlength $\delta_{\textrm{Ti-C}}$. Figure 9 compares the H-Ti-C7-H embedding PEC to that of a fully correlated AFQMC calculation without embedding or downfolding approximations; the inset of the figure shows that deviations from full AFQMC are less than about 0.1 eV, except at the shortest $\delta_{\textrm{Ti-C}}$, where it is $\simeq 0.24$ eV, and at an intermediate bondlength, $\delta_{\textrm{Ti-C}}=3.78$ Bohr, where it is $\simeq 0.15$ eV. The error bars in the inset reflect the stochastic error of the individual energy measurements (as opposed to the joint stochastic error). Fig. 10 shows the H-Ti-C7-H energy vs. $R_{v}$ for two fixed values of $R_{o}$ and a few values of $\delta_{\textrm{Ti-C}}$. The convergence in $R_{v}$ for fixed $R_{o}$ is similar to the longer H-Ti-C29-H chain (compare to Fig. 8) for $\delta_{\textrm{Ti-C}}=3.40$ and $4.35$ Bohr. However, for $\delta_{\textrm{Ti-C}}=3.0$ Bohr, using $R_{o}=9.7$ Bohr, the convergence in $R_{v}$ differs from the other two $\delta_{\textrm{Ti-C}}$ plotted. This is reflected in the inset of Fig. 9 in which it is observed that the greatest deviation from the fully correlated result occurs at this particular value for $\delta_{\textrm{Ti-C}}$ when using ($R_{o}$, $R_{v}$) = (9.7, 13.7) Bohr. It seems that for a slightly larger $R_{v}$ of 16.1 Bohr that the relative energy between $\delta_{\textrm{Ti-C}}=3.0$ Bohr and other $\delta_{\textrm{Ti-C}}$ may be more converged.

The left panel of Fig. 11 shows the basis set dependence of the fully correlated AFQMC PECs. The cc-pVTZ basis used for C atoms was modified similarly to the C cc-pVDZ basis (Fig. LABEL:sup-bas:custom-cc-pVTZ-basis). As seen in the figure, all basis sets give results within about 0.05 eV near equilibrium. Only small basis-set effects are found for the fully correlated treatment. In the right panel of Fig. 11, the deviation of embedded/downfolded results from the fully correlated results are plotted (similar to the inset of the left panel). Again, only small basis-set effects are observed and the embedding and downfolding effects are seen to be largely independent of basis set zeta-quality.

Results for the H-Ti-C7-H and H-Ti-C29-H chains are compared in Fig. 12. The custom cc-pVDZ basis set given in Fig. LABEL:sup-bas:custom-cc-pVDZ-basis-1.6 is used for H-Ti-C7-H (only at $\delta_{\textrm{Ti-C}}=3.0$ Bohr) to facilitate comparison to H-Ti-C29-H at that particular bond length. Cutoff radii for both systems were chosen as $R_{o}=9.7$ and $R_{v}=13.7$ Bohr. As indicated by Figs. 8 and 10, this choice of cutoffs yields well converged results with the possible exception of the compressed bond length, $\delta_{\textrm{Ti-C}}=3.0$ Bohr. The PECs of both systems are in good qualitative agreement. We note that the potential well is slightly deeper for H-Ti-C7-H than for H-Ti-C29-H despite the fact that the active subsystems, at the AFQMC level of theory, are nearly identical in both systems. We discuss this point further in Section V.

IV.3 Interaction of the Linear Chain with a Graphitic Substrate

There has been recent interest in studying the role of environmental interactions on carbon chain systems Wanko et al. (2016). Recent experiments have observed long linear carbon chains (LLCCs) confined within double walled carbon nanotubes (DWCNTs). Shi et al. (2016) Here we consider a simpler model of carbon chains near a graphitic planar substrate. Given the similar properties of the Ti-capped C29 and C7 chains in the previous section, we study H-Ti-C7-H chains on a finite-size graphitic ribbon, H22C56, as depicted in Fig. 13, where the dangling bonds of the boundary C atoms are saturated by H atoms. The C-C bond lengths are fixed to the experimentally observed value in graphene $\delta_{\textrm{C-C}}=2.68$ Bohr; the C-H bond lengths are set to 2.06 Bohr. The lateral dimensions of the ribbon are $14.0\times 29.9$ Bohr, excluding H terminations. The chain is placed parallel to the ribbon at a vdW height of 5.99 Bohr, taken from a carbon chain in Ref. 29, and within the reflection symmetry plane of the ribbon. For the equilibrium $\delta_{\textrm{Ti-C}}$ found above for H-Ti-C7-H in vacuum, we found only small residual atomic forces $\simeq 0.2$ eV/Bohr using vdW-including DFT/PBE-D2 NWCHEM calculations. This indicates that the system is close to its equilibrium geometry for the vacuum-relaxed H-Ti-C7-H chain geometry. The positions of all atoms, including both the H-Ti-C7-H chain atoms and the H22C56 graphitic substrate atoms, are given in Table LABEL:sup-tab:H-Ti-C7-HatH22-C56-geom.

RHF calculations for H-Ti-C7-H/H22C56 were done as follows. The atoms in the H-Ti-C7-H chain used the same basis sets as in vacuum (Fig. LABEL:sup-bas:custom-cc-pVDZ-basis). For the H22C56 substrate, however, we excluded d-functions for the C atoms to reduce the computational cost (Fig. LABEL:sup-bas:custom-cc-pVDZ-basis-substrate). FB localization was used for the RHF occupied sector as in the other calculations. FM localization for the virtual orbitals did not converge in ERKALE, however, so FB localization of the virtuals was used instead. Although the substrate FB virtuals have somewhat large $\sigma_{4}$ values, they are well separated from the chain FB virtuals. Moreover, in H-Ti-C7-H/H22C56, the chain localized orbitals closely resemble those of the vacuum H-Ti-C7-H chain, both in their centroid positions and in their $\sigma_{4}$ values. It was therefore straightforward to select $R_{o}$ and $R_{v}$ cutoffs within the chain (similar to those used in vacuum) and separate cutoffs to select the active region within the substrate.

The H-Ti-C7-H chain components of the active space are defined using similar in-chain radial cutoffs, $R_{o}$ and $R_{v}$, as in the previous section. If we required a similar criterion for the ribbon, we would define the in-plane active region using substrate radii, $R_{o}^{sub}=\sqrt{R_{o}^{2}-h^{2}}$ and $R_{v}^{sub}=\sqrt{R_{v}^{2}-h^{2}}$, where $h=5.99$ Bohr is the chain height above the ribbon, and where $R_{o}^{sub}$ and $R_{v}^{sub}$ are measured from the projected position of the Ti atom. In this case, occupied (virtual) orbitals in the ribbon would only be correlated if $R_{o}>h\simeq 6$ Bohr ($R_{v}>h\simeq 6$ Bohr). Thus, for $R_{o}=4.9$ Bohr in Fig. 10, all the ribbon occupied states would be frozen. For the choice $(R_{o},R_{v})=(9.7,13.7)$ Bohr in Fig. 9, this would yield $R_{o}^{sub}=7.6$ and $R_{v}^{sub}=12.3$ Bohr. The latter choice would result in larger calculations than we wanted to undertake in this paper. Instead, calculations were done setting $R_{o}^{sub}=R_{v}^{sub}=R^{sub}$ , where $R^{sub}$ is a common radius measured from the 6-fold hollow site. Calculations were done for three values of $R^{sub}=0$, 2.83 and 4.72 Bohr, where the latter two are illustrated by the circles in the inset of Fig. 14 (left panel). The $R^{sub-nu}=0$ case corresponds to the embedding approximation for all of the occupied ribbon states; all the virtual ribbon orbitals are excluded from the active space. Thus for $R^{sub}=0$, chain-ribbon interactions are present only in the 1-body contributions to the effective Hamiltonian ($\hat{V}^{\mathbb{I}-\mathbb{A}}$ in Eq. 13).

The left panel of Fig. 14 shows calculated AFQMC energies of H-Ti-C7-H / H22C56 for three values of $R^{sub}$, as a function of $R_{v}$ and for two Ti-C bond lengths. All the curves are for $R_{o}=9.7$ Bohr (the same $R_{o}$ that was used for the H-Ti-C7-H chain in vacuum in Figs. 11, and 12). For comparison, the analogous H-Ti-C7-H vacuum results from Fig. 10 are also shown. All the curves, for the same $\delta_{\textrm{Ti-C}}$, are seen to be within $\simeq 0.40$ eV of each other. The $R^{sub}=0$ and vacuum H-Ti-C7-H curves are nearly identical, which shows that the 1-body embedding contributions for $R^{sub}=0$ are very small. This is not too surprising given the large chain-ribbon separation of 5.99 Bohr. The $\delta_{\textrm{Ti-C}}=3.40$ and 4.35 Bohr curves are nearly parallel to each other and show similar convergence with increasing $R_{v}$, leading to rapid convergence of energy differences. The similar convergence of relative energies between correlated substrate treatments and the H-Ti-C7-H chain in vacuum indicates that converged in-chain $R_{o}$ and $R_{v}$ can be chosen by considering only the chain system in vacuum. This would be considerably less expensive than performing many trial calculations in the presence of a larger environmental system.

The H-Ti-C7-H / H22C56 PECs corresponding to the left panel of Fig. 14 are shown in the right panel of the same figure. For stretched $\delta_{\textrm{Ti-C}}$, the $R^{sub}=2.83$ and 4.72 Bohr PECs (with active states in the substrate) are about $\simeq 0.2-0.4$ eV larger than the $R^{sub}=0$ (frozen substrate) and vacuum-chain H-Ti-C7-H PECs. The correlated-substrate well depths are therefore deeper than those for the vacuum or frozen substrate case. The frozen substrate and the vacuum chain PECs are in very good agreement except at 3.78 Bohr. At that specific bond length, the RHF solution landscape is particularly complicated and so it is possible that the RHF trial wavefunction used for H-Ti-C7-H in vacuum is slightly different than the analogous RHF trial wavfunction used for H-Ti-C7-H / H22C56 in terms of the single-particle states corresponding to the chain. This conclusion is further supported by the very strong agreement between the frozen substrate and correlated substrate treatments at that particular $\delta_{\textrm{Ti-C}}$. We note that the PECs for the $R^{sub}=2.83$ and 4.72 Bohr treatments are in excellent agreement for all $\delta_{\textrm{Ti-C}}$ computed, suggesting that convergence in $R^{sub}$ has been achieved. Further calculations should be done, however, to check the convergence with increasing $R^{sub}$ to confirm this. Additionally, the effect of separate substrate occupied and virtual localization radii, $R_{o}^{sub}$ and $R_{v}^{sub}$, should be considered. Ribbon finite-size effects should also be examined using larger (especially longer) ribbons.

V Discussion

As shown in Sections III and IV, embeddding/downfolding energies converge rapidly with increasing cutoff radius $R_{v}$ for fixed $R_{o}$ (Figs. 4, 8, 10, and the left panel of Fig. 14). Thus, $R_{v}$ can be chosen naturally as $R_{v}=R_{o}+C$, where $C\simeq 2-6$ Bohr is a system-dependent constant. The convergence rate, in terms of $C$, is observed to be nearly the same for closely related systems. For example, all of the H-Ti-Cn-H systems studied in Sec. IV have $C\approx 4$ Bohr, including H-Ti-C7-H / C56H22 (independent of the in-plane cutoff, $R^{sub}$) as can be seen in Figs. 8, 10, 14 (left panel). Convergence is even more rapid for relative energies (inset of Fig. 8), due to cancellation of errors. In the case of O/H20, the most converged PEC in Fig. 5, with ($R_{o}$, $R_{v}$) = (7.8, 9.8) Bohr, agrees with the full AFQMC treatement to within 0.1 eV for almost all O-atom heights. For H-Ti-C29-H, the inset of Fig. 7 shows that the embedding results are nearly converged; the ($R_{o}$, $R_{v}$) = (9.7, 13.7) Bohr PEC deviates from the ($R_{o}$, $R_{v}$) = (14.7, 18.6) Bohr PEC by no more than 0.1 eV for all bond lengths except for the most compressed bond length where the deviation is about 0.2 eV. Similar behavior is observed for H-Ti-C7-H in Fig. 9 with a maximum deviation of about 0.24 eV, with most bond lengths deviating by 0.1 eV or less. Thus, local embedding/downfolding can be applied with little loss of quantitative accuracy.

Using the embedding/downfolding effective Hamiltonian significantly reduces the computational cost compared to the AFQMC treatment of the full Hamiltonian. For the H-Ti-C29-H chain system studied in Sec. IV.1 using ($R_{o}$, $R_{v}$) $=$ ($9,7$, $13.7$) Bohr, which was found to be converged for most $\delta_{\textrm{Ti-C}}$, the number of electrons per spin, N, is reduced from 65 (core states frozen) to 13, and the number of basis set functions, M, is reduced from 419 to only 100, reducing the computational cost by a factor of about 440, with AFQMC scaling as $\propto M^{2}N^{2}$. For H-Ti-C7-H /H22C56 studied in Section IV.3, the full system has $M$=625 and $N$ (per spin) = 144. Using localization radii ($R_{o}$, $R_{v}$) $=$ ($9,7$, $13.7$) Bohr and allowing the 6-fold hollow site nearest to the Ti atom to be correlated ($R^{sub}=2.83$ Bohr), $M$ becomes 130 and $N$ becomes 22, reducing the cost by nearly 3 orders of magnitude (a factor of 990). With such reductions, we were able to treat, with AFQMC, the H-Ti-C29-H and the H-Ti-C7-H /H22C56 systems using only modest computing resources. Thus, the effective system size which can be treated is greatly extended by the local embedding and effective downfolding approximation. We note that these computations were performed using the modest cc-pVDZ basis. For larger basis sets, the cost may be expected to be reduced by a similar factor, since the fraction of basis functions which are assigned to $\mathbb{A}$ will be roughly the same as for smaller basis sets, assuming a similar distribution of localized single-particle orbitals in both basis sets.

In Section IV.2, basis-convergence effects on the local embedding and downfolding approximation were considered. For fixed ($R_{o}$, $R_{v}$), the deviation of AFQMC results using embedding and downfolding from the full AFQMC treatment seems to be largely independent of basis set quality as seen in the right panel of Fig. 11. This can be exploited in order to study convergence in the localization radii while using a smaller, and therefore less expensive, basis set before using higher quality basis sets to obtain quantitatively accurate results.

Finite-size effects were studied in Section IV.2, comparing the PECs for H-Ti-C7-H and H-Ti-C29-H (Fig. 12). Good qualitative agreement was observed, with slight quantitative differences in the depth of the potential well. The geometry of H-Ti-C7-H is based on the geometry of H-Ti-C29-H with the chain truncated after 7 C atoms and with a H termination applied. As can be seen from the distribution and spatial extent of local orbitals in both systems (Fig. LABEL:sup-fig:H-Ti-Cn-H-centHist), the active regions have spatially similar localized occupied and virtual orbitals. As a result, the active regions of both chains are the same spatial size and have the same number of electrons and have basis set sizes that differ by only one. The inactive host regions differ, of course, containing 4 and 26 carbon atoms in the H-Ti-C7-H and H-Ti-C29-H chains respectively. Therefore the effective Hamiltonians of the two chains differ mostly in their mean-field embedding contributions, i.e., in ${\hat{V}}^{\mathbb{I}-\mathbb{A}}$ from Eq. 13. Unitarily localized Kohn-Sham orbitals, obtained from DFT, or natural orbitals, obtained using some many-body method could lead to different one-body size effects in ${\hat{V}}^{\mathbb{I}-\mathbb{A}}$ than the orbitals used here.

Orbital localization using the FM2 cost function with Hessian-based trust-region minimization is regarded as the method of choice for obtaining well-localized virtual orbitals.Høyvik et al. (2012) As mentioned, we used the ERKALE program, which uses gradient-based minimization, to localize both occupied and virtual orbitals in the present work. The consequence of using gradient-based methods to minimize the FM2 cost function is possible non-convergence, as is the case for the systems studied here. FM did not converge for H-Ti-C7-H /H22C56 and so FB, which is expected to produce virtual orbitals with long tails, was used for all H-Ti-C7-H systems. Despite the use of such virtual orbitals, rapid convergence in $R_{v}$ at fixed $R_{o}$ is observed. This suggests that local embedding/downfolding for a single active cluster is robust against basis sets which have not been optimally localized. Further study is needed to determine the effects of the system dimensionality on the convergence in $R_{v}$. In this regard, the one-dimensional geometry of the chain systems is certainly advantageous. It is also observed in all H-Ti-Cn-H systems studied here that the least local orbitals tend to be localized far from the Ti-C cluster (see Figs. 6 and  LABEL:sup-fig:H-Ti-Cn-H-centHist), providing some additional explanation for the rapid convergence observed in these systems. Using FM2 virtual orbitals is likely to be especially efficient for two- and three-dimensional system geometries.

VI Summary

In summary, a local embedding and effective downfolding approximation has been developed and implemented in which a local cluster is explicitly treated at the AFQMC level of theory. Previously, only local embedding of the occupied sector had been used in AFQMC, which requires the large virtual sector to be explicitly treated in all AFQMC computations. Here, we have developed a local effective downfolding scheme for the virtual sector and combined it with the local embedding approximation. The computational effort of performing AFQMC computations using embedding/downfolding is greatly reduced (by as much as 3 orders of magnitude in the examples studied here) from the full AFQMC treatment. Therefore, the effective system size which can be feasibly treated with AFQMC is significantly increased. We have studied the convergence systematics in the localization radii, $R_{o}$ and $R_{v}$, which define the embedding and downfolding regions. Relative energies computed using local embedding and effective downfolding are observed to converge rapidly to the full AFQMC result. Furthermore, the rate of convergence is seen to be similar for closely related systems such as the H-Ti-Cn-H systems studied in this paper. The convergence is seen to be largely independent of the zeta-quality of the basis set used.

We have used only RHF trial wave functions, $|{{\Psi_{T}}}\rangle$, to study the convergence systematics in the cutoff radii, $R_{o}$ and $R_{v}$. In some cases open-shell or multi-determinant $|{{\Psi_{T}}}\rangle$ may be required to achieve high accuracy AFQMC results. The embedding/downfolding approximation can be generalized to these cases. For example, in multi-determinant $|{{\Psi_{T}}}\rangle$’s, the determinants are typically expressed in terms of a single set of orthonormal canonical orbitals, which can then be localized by the same unitary transformation. Furthermore, it is expected that the systematics in $R_{o}$ and $R_{v}$ will be similar in these cases since the action of the embedding and downfolding transformation merely reduces the effective size of the single-particle Hilbert space which is spanned by a set of unitarily localized orbitals. Further study is necessary to determine the systematics for more general $|{{\Psi_{T}}}\rangle$.

Acknowledgements.

This work is supported by DOE (Grant No. DE-SC0001303), ONR (Contract No. N00014-17-1-2237), and the Simons Foundation. This work was performed using computing facilities at William & Mary which were provided by contributions from the National Science Foundation, the Commonwealth of Virginia Equipment Trust Fund and the Office of Naval Research. The Flatiron institute is a division of the Simons Foundation.