THE JOURNAL OF CHEMICAL PHYSICS 142, 044111 (2015)

Self-consistent embedding of density-matrix renormalization group wavefunctions in a density functional environment Thomas Dresselhaus,1 Johannes Neugebauer,1,a) Stefan Knecht,2,b) Sebastian Keller,2 Yingjin Ma,2 and Markus Reiher2,c) 1

Westfälische Wilhelms-Universität Münster, Theoretische Organische Chemie, Organisch-Chemisches Institut and Center for Multiscale Theory and Computation, Corrensstraße 40, 48149 Münster, Germany 2 ETH Zürich, Laboratorium für Physikalische Chemie, Vladimir-Prelog-Weg 2, 8093 Zürich, Switzerland

(Received 5 September 2014; accepted 5 January 2015; published online 27 January 2015) We present the first implementation of a density matrix renormalization group algorithm embedded in an environment described by density functional theory. The frozen density embedding scheme is used with a freeze-and-thaw strategy for a self-consistent polarization of the orbital-optimized wavefunction and the environmental densities with respect to each other. C 2015 AIP Publishing LLC. [http://dx.doi.org/10.1063/1.4906152]

I. INTRODUCTION

Chemical reactions are local phenomena, in which usually only one or two bonds are formed or broken at a time. For this reason, the cluster approach was established in computational chemistry, in which the atomistic rearrangements under study are embedded in a structural model that is as large as necessary to comprise all electronic effects on the reactive subsystem (see, for instance, the extensive work on such cluster models for active sites of metalloproteins by Siegbahn and Himo1). In general, long-range electrostatic interactions may not be neglected and schemes have been devised to incorporate them (most prominent are quantummechanical/molecular-mechanical embedding2,3 and polarizable continuum models4). For the appropriate description of extensive electronic coupling of a subsystem to its environment (mediated, for example, through π-conjugated molecular substructures), more involved embedding schemes are needed. Clearly, the embedding of a quantum system of interest into a more or less electronically noninnocent environment is the realm of open-system quantum theory.5 The central quantity of this well-known theory is the reduced density matrix of a subsystem rather than the wavefunction of a stationary state. This formalism has, recently, been adapted for an embedding framework.6,7 However, it is not straightforward to employ (for instance, the environment had to be approximated in a mean-field fashion). Practical suggestions for the applicability to large chemical systems have been made.8 Considering the success of density functional theory (DFT) for routine applications in computational chemistry, a feasible approach that is able to describe chemically complex and structured environments should also be based on DFT. In the frozen density embedding9 (FDE) scheme, which is a)Electronic mail: [email protected] b)Electronic mail: [email protected] c)Electronic mail: [email protected]

0021-9606/2015/142(4)/044111/7/$30.00

based on subsystem DFT,10,11 a subsystem is calculated in a Kohn–Sham (KS) approach by adding an effective potential to the KS potential of the isolated subsystem. It takes into account all interactions with the environment, assuming a constant (frozen) environmental electron density. In FDE, not only the exchange–correlation energy contribution but also the non-interacting kinetic energy term is calculated from an approximate explicit density functional. No knowledge of orbitals of the supersystem is needed for a FDE/subsystem DFT calculation, nor do the orbitals of the active subsystem need to be orthogonal to those of the environmental subsystems. In this way, large environments can be considered. Due to the limited accuracy of kinetic energy functionals for chemical systems, FDE cannot be directly applied to subsystems which are connected via bonds with covalent character.12 Extensions to resolve this limitation are available. Often, potential-reconstruction methods have been proposed in this context.13–17 A related embedding scheme by Manby et al.18 resolves this issue by imposing orthogonality constraints between orbitals of the active subsystem and the environment. While this strategy avoids the need for explicit density functionals for the evaluation of the nonadditive kinetic energy, it requires that an orbital optimization for the full system (active subsystem plus environment) is feasible, which can be a severe constraint. However, a straightforward basis set truncation scheme for the actual embedding calculation is readily available.19 In FDE, all electron-correlation effects between different subsystems are captured by the density functional. In principle, all correlation effects, including the polarization of the wavefunction of the active subsystem, can be accurately described for the (admittedly unknown) exact functional. Considering the fact that approximate density functionals may yield rather inaccurate results, it is desirable to describe a (active) subsystem by systematically improvable wavefunction theory (WFT). Although based on DFT, the FDE framework also allows for treating one subsystem with a correlated wavefunction method and the rest of the system

142, 044111-1

© 2015 AIP Publishing LLC

044111-2

Dresselhaus et al.

J. Chem. Phys. 142, 044111 (2015)

with DFT.20 Several research groups implemented such a WFT-in-DFT embedding; for instance, to facilitate completeactive-space self-consistent-field (CASSCF) in periodic21,22 and non-periodic23 DFT, Møller–Plesset perturbation theory in DFT,24 coupled cluster in DFT,25 CASSCF with secondorder perturbation theory (CASPT2) in DFT,23,26 and most recently quantum Monte Carlo in DFT27 calculations (for overviews see Refs. 28–30). Within a comparatively short time, the density matrix renormalization group (DMRG) algorithm31,32 has become a standard WFT approach in quantum chemistry (for reviews see Refs. 33–38), which is under continuous development (see Refs. 39 and 40 for most recent work), providing interesting new perspectives for the description of electronic structures dominated by static electron correlation.41–49 DMRG is able to iteratively converge to the exact solution of the electronic Schrödinger equation in a given active orbital space with polynomial rather than factorial cost by which the full configuration (and thus the CASSCF) approach is plagued. However, the benefits of this iterative protocol come with a caveat that has long prohibited the extension of its applicability to compact and strongly correlated molecules like transition metal complexes and clusters. This caveat is rooted in the matrix-product-state (MPS) structure of the optimized DMRG wavefunction. As a consequence, DMRG should only be applied to pseudo-one-dimensional molecular structures. It could be shown50 that its benefits outweigh these formal objections; still, care must be taken to control the parameters that determine the accuracy of an iterative DMRG calculation.51 It has now been fully recognized that accuracy- and performance-wise, DMRG can be a substitute for the CASSCF (and CASPT2) methodology pioneered by the Lund group52 for cases which require much larger active orbital spaces than those accessible by this standard approach. With DMRG, active orbital spaces are in reach that are about five times larger than those accessible to CASSCF. This feature makes DMRG an ideal target for a multireference WFT-in-DFT embedding method. Here, we report the first implementation of a DMRG-inDFT embedding (including orbital relaxation in the DMRG step; DMRG-SCF), which follows the FDE framework in a freeze-and-thaw53 manner, to establish a fully self-consistent embedding. Often (but not always27) is the environment polarized by employing a DFT-in-DFT approach, then kept fixed in the WFT calculations, and a fixed embedding potential is applied.23,25 Alternatively, Carter and co-workers imposed constraints on the total54 or environmental22 density.

optimizing the energy of the total system with respect to the density ρact(r). Through an iterative exchange of the role of active and environmental densities known as “freeze-andthaw” (f&t), FDE becomes equivalent to subsystem DFT.28 The method can be extended in such a way that the active subsystem is described in terms of a wavefunction Ψact, while the environment is described in terms of density functional theory. This idea was introduced in a hybrid-energy fashion,24 but can also be more formally regarded as a non-standard DFT variant in which the density of the active system is represented in terms of a correlated many-electron wavefunction, while a Kohn–Sham surrogate system of noninteracting fermions is taken to represent the density of the environment.20 The total electronic energy Etot can then be calculated in the following way: WFT KS–DFT OF–DFT Etot = Eact + Eenv + Eint .

(2)

WFT Here, Eact is the energy of the active subsystem, which, for most wavefunction-based quantum chemical methods, is calculated as the expectation value, WFT Eact = ⟨Ψact| Hˆ act|Ψact⟩,

(3)

with Hˆ act being the isolated-system Hamiltonian for the WFT here as active part. To be precise, we can identify Eact the energy of the isolated active subsystem, but having assumed a wavefunction that is fully relaxed with respect to the environment. KS–DFT Eenv is the Kohn–Sham energy evaluated for the OF–DFT environmental density, and Eint is the interaction-energy, which is given as an explicit density functional in FDE (and thus calculated as an orbital free (OF-)DFT contribution), OF–DFT Eint [ρact, ρenv] = Tsnad[ρact, ρenv] nad int + Exc [ρact, ρenv] + Eelstat [ρact, ρenv],

(4) with the non-additive non-interacting kinetic energy Tsnad[ρact, ρenv] = Ts[ρact + ρenv] −Ts[ρact] −Ts[ρenv],

(5)

the non-additive exchange-correlation energy nad Exc [ρact, ρenv] = Exc[ρact + ρenv] − Exc[ρact] − Exc[ρenv],

(6)

and the electrostatic interaction terms   ρact(r1)ρenv(r2) int dr1dr2 Eelstat [ρact, ρenv] = |r1 − r2|  +

env ρact(r1)vext (r1)dr1

 +

II. WFT-IN-DFT EMBEDDING THEORY

We briefly review the theoretical framework that is the basis of our implementation. In FDE, the total electron density is represented by a sum of subsystem densities, here, denoted as an active (ρact) and environmental (ρenv) densities, ρtot(r) = ρact(r) + ρenv(r).

(1)

A goal of FDE is determining ρact(r), given a certain ρenv(r). The active subsystem density is found by variationally

act ρenv(r1)vext (r1)dr1,

(7)

I where vext is the potential due to the nuclei of subsystem I (I = env, act). WFT Eact in Eq. (3) is calculated by diagonalizing a matrix representation of Hˆ act in a matrix-product-state basis,55,56 which is iteratively optimized by the DMRG algorithm. By contrast to almost all previous quantumchemical DMRG implementations, our implementation51,57 calculates the matrix representation of Hˆ act directly from its matrix-product-operator (MPO) form (it is thus a second

044111-3

Dresselhaus et al.

J. Chem. Phys. 142, 044111 (2015)

generation quantum-chemical DMRG implementation). The MPS structure of the electronic state can be transferred to the expression for the Hamiltonian  Hσσ ′|σ⟩⟨σ ′|, (8) Hˆ act = σ,σ ′

as the coefficients Hσσ ′, defined with respect to subsystem occupation number vectors |σ⟩ and |σ ′⟩ of the total state defined on the CAS, can be encoded in matrix-product form  σl σ′ σ 1σ ′ σ Lσ′ (9) Hσσ ′ = H1b 1 ··· Hb bl ··· Hb 1L , 1

b1, ...,b L−1

l−1 l

L−1

where the indices bi label the local states on orbitals from the CAS and L denotes the total number of orbitals in the CAS. Combining Eqs. (8) and (9) yields   σl σ′ σ 1σ ′ σ Lσ′ H1b 1 ··· Hb bl ··· Hb 1L |σ⟩⟨σ ′|. (10) Hˆ act = σ,σ ′ b1, ...,b L−1

1

l−1 l

L−1

Contracting over the orbital-state indices σi in σ by defining  σ σ′ l Hbl l −1b l = Hb bl |σl ⟩⟨σl′|, (11) σ l ,σ ′l

l −1 l

finally yields the product formulation of the Hamiltonian  1 Hˆ act = H1b ··· Hbl l −1b l ··· HbLL−11. (12) 1 b1, ...,b L−1

In which, the Hamiltonian is constructed from elementary creation and annihilation operators defined on the orbitals. OF–DFT The functional derivative of Eint with respect to a subsystem density yields an effective embedding potential for this subsystem, which includes all interactions with the other subsystems. For the active subsystem, this reads as act vemb [ρact, ρenv](r) = vkin[ρact + ρenv](r) − vkin[ρact](r)

+ vxc[ρact + ρenv](r) − vxc[ρact](r)  ρenv(r2) env + dr2 + vext (r), (13) |r2 − r| where vkin results from the functional derivative of the noninteracting kinetic energy and vxc is the exchange–correlation potential. Finding the optimum energy of the total system thus requires finding a density that is the ground-state density for the Hamiltonian n act  ′ act Hˆ act = Hˆ act + vemb [ρact, ρenv](ri ), (14) i

with nact being the number of electrons in the active system. A typical solution for the active system will proceed by variationally minimizing the pseudo-energy WFT ′ Eact = ⟨Ψact| Hˆ act |Ψact⟩,

(15)

and then inserting the resulting wavefunction into Eq. (3). This pseudo-energy can be identified15 as the hypothetical energy of the active subsystem with an additional external potential act corresponding to vemb [ρact, ρenv](r). For post-Hartree–Fock methods like DMRG, this means in practice that one has to add this extra potential to the one-electron part of the (generalized) Fock operator. Its eigenstates are the molecular orbitals from which then one- and two-electron integrals are calculated

that ultimately enter the MPO representation of the secondquantized Hamiltonian in Eq. (12). Similarly, the environment env will experience the additional potential vemb [ρenv, ρact](r) in the Kohn–Sham equations determining the orbitals for the environment. In the following, we will discuss several aspects relevant for practical calculations of this type and describe how they are handled in our implementation. act A. Self-consistent evaluation of vemb

act Since vemb depends on ρact and on ρenv, a fully consistent act evaluation requires that vemb is updated whenever one of the two densities changes in an (iterative) calculation. The env same holds for vemb . The active density ρact changes in every self-consistency cycle of the wavefunction calculation (called “inner self-consistency” in the following), and ρenv changes in every f&t iteration (“outer self-consistency”). Clearly, there are several possible strategies to arrive at a fully self-consistent solution.58 In the past, many WFT-in-DFT approaches have either introduced constraints on the total54 or environmental22 density, or assumed that the wavefunction-derived density, upon convergence of the freeze-and-thaw iterations, can be approximated and thus substituted by a KS–DFT density for the active part calculated in a subsystem-DFT manner.25 In this latter strategy, one replaces the expensive outer WFT-in-DFT self-consistency steps by a much simpler outer DFT-in-DFT self-consistency. Moreover, the inner self-consistency steps are completely avoided. But this scheme may fail if the DFT density for the active part differs substantially from the wavefunction-based density. Particularly, for spin-polarized systems, where DFT tends to over-delocalize the spin density, this can introduce uncontrollable errors. Hence, we have implemented a fully self-consistent variant of WFT-in-DFT embedding. Since evaluating the density in WFT can be cumbersome or computationally expensive, here, we chose a route to fully self-consistent act solutions that avoids building vemb during the self-consistent iterations of the WFT part. I.e., we avoid the inner selfact consistency cycles for vemb . Instead, we construct the embedding potential in the i-th f&t iteration, denoted as act,(i) vemb , from the fixed density ρ(i−1) act obtained in the (i − 1)-th f&t iteration. Upon convergence of the f&t iterations, the densities of two subsequent iterations will be the same, and also the inner self-consistency cycles will not change ρact anymore (within the numerical threshold chosen). We should stress that also the environmental embedding potential at the end of this procedure is constructed with the fully relaxed wavefunction-derived density for the active part.

act B. Numerical calculation of vemb

act As is apparent from Eq. (13), vemb can be easily evaluated for a set of real-space grid points if explicitly densitydependent approximations are used for vkin and for vxc. In fact, even reconstructed potentials17 can be used in this way. Contributions to (generalized) Fock matrix elements due to the embedding potential can then be obtained similarly to the corresponding exchange–correlation contributions, namely, by numerical integration on a grid. The matrix representation

044111-4

Dresselhaus et al.

J. Chem. Phys. 142, 044111 (2015)

of the operator associated with this potential is, thus, obtained by numerical integration as act act (vemb ) µ,ν = ⟨ χ µ |vemb | χν ⟩  act ≈ χ µ (rk )vemb (rk ) χν (rk )wk ,

(16)

k

for a basis set { χ µ } and an integration grid with weights {wk } at grid points {rk }. The choice of the integration grid can have a strong influence on the calculated properties. To ensure comparability between energy terms of the full system and the subsystems, we use the same supermolecular integration grid for all numerically evaluated contributions to the embedding potential and energies. This grid is created from atom-centered grids for all atoms in the total system. The potentials vkin and vxc will be easy to calculate for local density approximation (LDA) and generalized gradient approximation (GGA) if the density, the density gradient, and its Hessian are available on a grid. This can be problematic in practice if the density of the WFT method is not readily accessible, which can be the case for correlated coupledcluster or multi-reference perturbation theory approaches. For the terms involving the sum of active and environment densities, the corresponding ingredients can just be added if the same grid is used for all quantities. In addition, the electrostatic potential for the corresponding subsystem needs to be available on the same grid. This is straightforward if the reduced one-particle density matrix and the corresponding one-electron integrals are available. C. Calculation of the energy

Once the f&t procedure is considered converged, the total energy can be calculated from the wavefunctions and densities obtained. As noted above and discussed, e.g., in Ref. 23, it is for computational convenience tempting to exploit the quantities calculated for the pseudo-energy in Eq. (15) when optimizing the fully relaxed wavefunction of the active part. WFT Eact is related to the energy of the embedded subsystem WFT Eact by  act WFT WFT + (vemb ) µ,ν Pact Eact = Eact µ,ν µ,ν

 WFT + = Eact

act vemb (r)ρact(r)dr  WFT act ≈ Eact + vemb (ri )ρact(ri )wi ,

(17)

i

where Pact is the reduced one-particle density matrix of the active subsystem. As discussed in Ref. 23, the integral over the embedding potential contains several terms OF–DFT also appearing in Eint , either in exact (electrostatic) or approximate, linearized (exchange–correlation, kinetic energy) form. Nevertheless, the correct total energy requires to subtract the corresponding integral from the pseudo-energy,  tot WFT act E = Eact − vemb (r)ρact(r)dr KS–DFT OF–DFT + Eenv + Eint .

(18)

In our implementation, the matrix representation of the embedding potential is calculated alongside the core Hamiltonian and added to the corresponding matrix. In this way, no modification of the SCF algorithm is necessary act to include the effects of vemb on the wavefunction. The WFT corresponding pseudo-energy Eact does not carry a direct physical meaning. Still, it provides a simple tool for analyzing the strength of the influence of the environment. The WFT subsystem energy Eact can easily be obtained by subtracting act the contribution due to vemb . III. COMPUTATIONAL METHODOLOGY

Our implementation provides an interface between a development version of MOLCAS,59 the quantum-chemical version of the MAQUIS DMRG program,51,57 and the ADF program package.58,60,61 Our quantum chemical version of MAQUIS51,57 explicitly implements an MPO-based formalism for the full second-quantized electronic Hamiltonian and exploits all technical features of the underlying MAQUIS program.62 We extended the scripting framework PyADF63 accordingly to allow for easy scripting of the calculations. In our implementation, an active subsystem described by MPSs as optimized with DMRG is embedded into an environment described with DFT. The following steps are carried out: (0) 1. Calculate initial DFT densities ρ(0) env and ρact (as well as DFT density gradients, DFT Hessians, and electrostatic potentials) for the isolated subsystems; set iteration counter i = 1. 2. Calculate an embedding potential for the environmental env,(i) (i−1) subsystem, vemb [ρenv , ρ(i−1) act ](r). 3. Perform an embedded DFT calculation on the environenv,(i) mental subsystem with vemb and then calculate ρ(i) env (as well as its gradient, Hessian, and the electrostatic potential). 4. Calculate an embedding potential for the active subsystem, act,(i) (i−1) vemb [ρact , ρ(i) env](r). 5. Carry out the 4-index transformation to obtain the twoelectron integrals in the molecular orbital basis of the active subsystem; the basis-functions to molecular-orbitals transformation of the one-electron integrals yields the WFT part that carries the embedding information. 6. Perform an embedded DMRG calculation on the active act,(i) subsystem using vemb and then calculate ρ(i) act (as well as its gradient, Hessian, and the electrostatic potential) from the DMRG wavefunction. 7. Increase iteration counter i by 1. Repeat steps 2–6 until convergence is achieved.

In our implementation, all embedding potentials are calculated with ADF. A supermolecular integration grid is used throughout, as described in Sec. II B. For the embedded DMRG calculations, the embedding potential is read in by MOLCAS, transformed to a matrix representation, and added to the core Hamiltonian in matrix representation as a potential which is unchanged in the DMRG calculation. Using the reduced one-particle density matrix of the converged DMRG calculation, we obtain the DMRG density, its gradient and

044111-5

Dresselhaus et al.

J. Chem. Phys. 142, 044111 (2015) TABLE I. DMRG-in-DFT calculations DMRG(8,8)[1000]-SCF and the PW91/PW91k functional on methane embedded in an ammonium-ion enviWFT denotes ronment at different distances. Energies are given in Hartree. E iso the DMRG(8,8)[1000]-SCF energy of the isolated active subsystem. The energy difference in the last column is a measure for the polarization of the DMRG wavefunction.

FIG. 1. A methane molecule (left) in the presence of an ammonium ion (right) at the equilibrium distance of 5.68 Å.

Hessian, and the electrostatic potential. These quantities are exported on the integration grid. These data are not only needed for in the calculation of the embedding potential for the environment, but also for the embedding potential in the DMRG calculation of the next f&t cycle. The embedded DMRG calculations are performed in two steps. First, orbitals are produced in an embedded Hartree–Fock calculation. These are then taken as a starting guess in the DMRG-SCF calculation (or in a CASSCF calculation, which we choose as a reference in this work). For the DMRG calculations, we adopt the short-hand notation DMRG(e,o)[M]-SCF to denote the numbers of active electrons e and active orbitals o in the active space as well as the number of renormalized block states M, which determines the accuracy of the DMRG results. Hence, WFT methods accessible in our implementation include Hartree–Fock, CASSCF, and DMRG (with and without orbital optimization: DMRG-SCF and DMRG-CI, respectively). Although the embedding potential for the WFT subsystem depends on the density calculated from the wavefunction, we keep it fixed during the wavefunction calculation, which simplifies the implementation. The error introduced by this approximation vanishes as soon as the wavefunction converges with respect to the f&t procedure; thus, no additional error is produced in the final results. All DFT calculations with ADF are performed with the PW9164 functional in a TZP65 basis. Also, the embedding potentials were calculated with ADF with the PW91k66 functional for the nonadditive kinetic energy. A supermolecular integration grid was used throughout. All DMRG-SCF calculations were performed with MOLCAS using a cc-pVTZ67 basis.

IV. RESULTS A. Embedding in a weakly interacting environment

As a first benchmark, for which we can still obtain CASSCF reference data, we considered the polarization of a methane molecule embedded in the environment of an ammonium ion as depicted in Figure 1. The structure of the methane–ammonia system was optimized with the ORCA program68 employing the BP8669,70-D3BJ71 functional plus dispersion correction and a Def2-TZVP basis set72 in

r(C–N) (Å)

WFT − E WFT Eact iso

WFT − E WFT E act iso

10 5.68 (eq.)

−5.31 × 10−1 −9.34 × 10−1

7.44 × 10−5 5.97 × 10−4

TABLE II. Comparison of dipole moments in debye of methane embedded in an ammonium-ion environment obtained from DMRG-in-DFT with DMRG(8,8)[1000]-SCF and the PW91/PW91k functional and analogous DFT-in-DFT calculations, respectively. r(C–N) (Å) 10 5.68 (eq.)

DMRG-in-DFT

DFT-in-DFT

0.11 0.33

0.12 0.36

combination with density fitting.73 This model system is a prototypical case to which an FDE-based embedding strategy with presently available approximate kinetic energy density functionals can be reliably applied. All DMRG-SCF results for an active space of eight electrons (the 1s electrons of C are kept frozen) in eight orbitals, CAS(8,8), are compiled in Table I. At equilibrium distance (“eq.”), the effect on the subsystem energy is by an order of magnitude more pronounced than at a distance of 10 Å between the molecule centers. As expected, the polarization of the wavefunction is more pronounced for the calculation on a closer distance. The wavefunction-polarization effect is more obvious when investigating the dipole moment of the methane molecule, as shown in Table II. Analogous CASSCF-in-DFT calculations (not shown) yield the same results up to the presented accuracy. Furthermore, the DMRG-SCF-in-DFT data compare very well with the corresponding DFT-in-DFT results (right column of Table II). B. The dipole moment of a chain of HCN molecules

To demonstrate that our new implementation reaches beyond the capabilities of a CASSCF-in-DFT embedding approach, we studied the dipole moment of a HCN molecule induced by the presence of further HCN molecules (see Figure 2 for the full tetrameric super-structure and the internuclear distances chosen). We selected a HCN tetramer as a benchmark system as it is known that HCN chains show a strong cooperative behaviour in this regard.28 Table III comprises the calculated dipole moments obtained using

FIG. 2. Structure of the (HCN)4 benchmark system. Internuclear distances, which specify this linear HCN chain, are given in Å. The shaded regions denote the two HCN molecules that serve as an environment for the dimer (HCN)2 (i.e., A + B) in the DMRG-in-DFT embedding calculation.

044111-6

Dresselhaus et al.

J. Chem. Phys. 142, 044111 (2015)

TABLE III. Dipole moment µ in debye per HCN molecule in the presence and absence of other HCN molecules as depicted in Figure 2. Isolated HCN (central) dimer and tetramer calculations are compared to the corresponding ones for the embedded monomer and dimer, respectively, in which the environment “Env.” is represented by a DFT density.

Method

DMRG(10,9)[4096]-SCF DFT DMRG(10,9)[4096]-SCF DFT DMRG(10,9)[4096]-SCF DFT DMRG(20,18)[4096]-SCF DMRG(20,18)[4096]-SCF DFT DMRG(40,36)[2048]-SCF

Active fragment(s) Monomer A A A A B B Dimer A+B A+B A+B Tetramer All

Env.

µ (per HCN)

... ... DFT DFT DFT DFT

3.09 2.96 3.45 3.42 3.29 3.32

... DFT DFT

3.44 3.77 3.93

...

3.81

either a supermolecular or the DMRG/DFT-in-DFT approach. The reference DMRG-SCF calculation on the (central) supermolecular dimer comprised an active space of 20 electrons in 18 orbitals (twice the valence complete active space for one HCN molecule, which is a CAS(10,9)). We find only small differences in the dipole moments between the DFT and the DMRG calculations. A qualitatively correct polarization of the embedded HCN molecule compared to its isolated counterpart can be observed. An advantage of the subsystem approach is the possibility to investigate the properties of only one subsystem in the presence of other subsystems. To demonstrate this feature, we collected the different dipole moments for the two different HCN molecules A and B in Table III. As our benchmark system can be systematically enlarged, we next consider the HCN tetramer, in which we embed the central HCN dimer into two flanking HCN molecules (see Figure 2). The reference DMRG-SCF calculation on

the full tetramer structure comprised an active space of 40 electrons in 36 orbitals. The data in Table III clearly illustrate the cooperative effect on the dipole moment in HCN chains. The dipole moment per HCN unit increases by approximately 10% from a single HCN molecule, to the isolated dimer. It further increases by a similar percentage from the isolated dimer to the tetramer. We further observe a qualitatively correct polarization of the embedded dimer (HCN)2 compared to its isolated counterpart. We also note that the DFT-in-DFT data deviate only slightly from the DMRG-in-DFT results. The importance of taking into account a fully selfconsistent polarization of the environment by means of the f&t procedure is indicated in Figure 3. The f&t procedure shows a rapid convergence. The error decreases by remarkable 98.6% in each f&t cycle compared to the converged result. V. CONCLUSIONS

Our new DMRG-in-DFT embedding approach combines the advantage of DMRG to be able to handle large active orbital spaces with the reliable quantum mechanical FDE scheme in a f&t approach for the polarization of the embedding density. This approach facilitates accurate calculations on systems with strong static correlation embedded in environments whose effects are important beyond a classical (purely electrostatic) description. The pilot applications reported in this work are mainly proofof-principle calculations. However, we propose the HCNchain model system as a suitable benchmark system for WFT-in-DFT calculations in general as (i) its size can be easily extended to create a computationally more demanding benchmark, (ii) the not strongly interacting HCN molecules do not lead to errors that are due to the currently available approximation for the kinetic energy functional in a FDE framework, and (iii) they even pose an interesting physical problem as the dipole moment is significantly affected by neighboring HCN molecules. Our algorithm offers various options for future developments. For instance, ground- and excited-state calculations on molecules that demand the consideration of large active orbital spaces embedded into a highly structured environment are accessible. An example is an extended chromophore embedded in a bio-macromolecule (such as retinal in rhodopsin or the special pair in photosystem II). Clearly, dynamic-correlation effects also need to be considered for quantitative results—either on top of the converged DMRGin-DFT results by a multi-reference perturbation calculation or before the DMRG-in-DFT calculation through a dressed Hamiltonian—and such developments are currently underway in our laboratories. ACKNOWLEDGMENTS

FIG. 3. Error of the pseudo-energy in Hartree with respect to a converged f&t procedure for a DMRG(10,9)-SCF-in-DFT calculation on a HCN molecule embedded into the environment of another HCN molecule as a function of the f&t cycle number.

T.D. and J.N. acknowledge funding by the Deutsche Forschungsgemeinschaft through SFB 858. This work was supported by ETH Zürich (Research Grant No. ETH-34 12-2) and by the Schweizerischer Nationalfonds (Project No. 200020_156598). Support by COST Action CODECS is gratefully acknowledged.

044111-7 1P.

Dresselhaus et al.

E. Siegbahn and F. Himo, WIREs Comput. Mol. Sci. 1, 323 (2011). M. Senn and W. Thiel, Top. Curr. Chem. 268, 173 (2007). 3H. M. Senn and W. Thiel, Angew. Chem., Int. Ed. 48, 1198 (2009). 4B. Mennucci, WIREs Comput. Mol. Sci. 2, 386 (2012). 5A. Amann and U. Muller-Herold, Offene Quantensysteme: Die Primas Lectures (Springer-Verlag, Berlin, Heidelberg, 2011). 6G. Knizia and G. K.-L. Chan, Phys. Rev. Lett. 109, 186404 (2012). 7G. Knizia and G. K.-L. Chan, J. Chem. Theory Comput. 9, 1428 (2013). 8Q. Sun and G. K.-L. Chan, J. Chem. Theory Comput. 10, 3784 (2014). 9T. A. Wesolowski and A. Warshel, J. Phys. Chem. 97, 8050 (1993). 10G. Senatore and K. R. Subbaswamy, Phys. Rev. B 34, 5754 (1986). 11P. Cortona, Phys. Rev. B 44, 8454 (1991). 12S. Fux, K. Kiewisch, C. R. Jacob, J. Neugebauer, and M. Reiher, Chem. Phys. Lett. 461, 353 (2008). 13O. Roncero, M. P. de Lara-Castells, P. Villarreal, F. Flores, J. Ortega, M. Paniagua, and A. Aguado, J. Chem. Phys. 129, 184104 (2008). 14J. D. Goodpaster, N. Ananth, F. R. Manby, and T. F. Miller, J. Chem. Phys. 133, 084103 (2010). 15C. Huang and E. A. Carter, J. Chem. Phys. 135, 194104 (2011). 16C. R. Jacob and L. Visscher, J. Chem. Phys. 128, 155102 (2008). 17S. Fux, C. R. Jacob, J. Neugebauer, L. Visscher, and M. Reiher, J. Chem. Phys. 132, 164101 (2010). 18F. R. Manby, M. Stella, J. D. Goodpaster, and T. F. Miller, J. Chem. Theory Comput. 8, 2564 (2012). 19T. A. Barnes, J. D. Goodpaster, F. R. Manby, and T. F. Miller, J. Chem. Phys. 139, 024103 (2013). 20T. A. Wesolowski, Phys. Rev. A 77, 012504 (2008). 21T. Klüner, N. Govind, Y. A. Wang, and E. A. Carter, Phys. Rev. Lett. 86, 5954 (2001). 22P. Huang and E. A. Carter, J. Chem. Phys. 125, 084102 (2006). 23C. Daday, C. König, O. Valsson, J. Neugebauer, and C. Filippi, J. Chem. Theory Comput. 9, 2355 (2013). 24N. Govind, Y. A. Wang, A. J. R. da Silva, and E. A. Carter, Chem. Phys. Lett. 295, 129 (1998). 25A. S. P. Gomes, C. R. Jacob, and L. Visscher, Phys. Chem. Chem. Phys. 10, 5353 (2008). 26D. K. Kanan, S. Sharifzadeh, and E. A. Carter, Chem. Phys. Lett. 519-520, 18 (2012). 27C. Daday, C. König, J. Neugebauer, and C. Filippi, ChemPhysChem 15, 3205 (2014). 28C. R. Jacob and J. Neugebauer, WIREs Comput. Mol. Sci. 4, 325 (2014). 29F. Libisch, C. Huang, and E. A. Carter, Acc. Chem. Res. 47, 2768 (2014). 30A. S. P. Gomes and C. R. Jacob, Annu. Rep. Prog. Chem., Sect. C: Phys. Chem. 108, 222 (2012). 31S. R. White, Phys. Rev. Lett. 69, 2863 (1992). 32S. R. White, Phys. Rev. B 48, 10345 (1993). 33Ö. Legeza, R. M. Noack, J. Sólyom, and L. Tincani, Lect. Notes Phys. 739, 653 (2008). 34K. H. Marti and M. Reiher, Z. Phys. Chem. 224, 583 (2010). 35S. Sharma and G. K.-L. Chan, J. Chem. Phys. 136, 144105 (2012). 36G. K.-L. Chan, WIREs: Comput. Mol. Sci. 2, 907 (2012). 37Y. Kurashige, Mol. Phys. 112, 1485 (2014). 38S. Wouters and D. van Neck, Eur. Phys. J. D 68, 272 (2014). 39S. Wouters, W. Poelmans, P. W. Ayers, and D. V. Neck, Comput. Phys. Commun. 185, 1501 (2014). 40S. Knecht, O. Legeza, and M. Reiher, J. Chem. Phys. 140, 041101 (2014). 41K. Boguslawski, K. H. Marti, Ö. Legeza, and M. Reiher, J. Chem. Theory Comput. 8, 1970 (2012). 2H.

J. Chem. Phys. 142, 044111 (2015) 42K.

Boguslawski, P. Tecmer, O. Legeza, and M. Reiher, J. Phys. Chem. Lett. 3, 3129 (2012). 43K. Boguslawski, P. Tecmer, G. Barcza, Ö. Legeza, and M. Reiher, J. Chem. Theory Comput. 9, 2959 (2013). 44Y. Kurashige, G. K.-L. Chan, and T. Yanai, Nat. Chem. 5, 660 (2013). 45T. V. Harris, Y. Kurashige, T. Yanai, and K. Morokuma, J. Chem. Phys. 140, 054303 (2014). 46T. N. Lan, Y. Kurashige, and T. Yanai, J. Chem. Theory Comput. 10, 1953 (2014). 47S. Sharma, T. Yanai, G. H. Booth, C. J. Umrigar, and G. K.-L. Chan, J. Chem. Phys. 140, 104112 (2014). 48P. Tecmer, K. Boguslawski, Ö. Legeza, and M. Reiher, Phys. Chem. Chem. Phys. 16, 719 (2014). 49S. Wouters, T. Bogaerts, P. Van Der Voort, V. Van Speybroeck, and D. Van Neck, J. Chem. Phys. 140, 241103 (2014). 50K. H. Marti, I. M. Ondik, G. Moritz, and M. Reiher, J. Chem. Phys. 128, 014104 (2008). 51S. Keller and M. Reiher, Chimia 68, 200 (2014). 52B. O. Roos, “Multiconfigurational quantum chemistry for ground and excited states,” Radiation Induced Molecular Phenomena in Nucleic Acids, Challenges and Advances in Computational Chemistry and Physics, edited by M. K. Shukla and J. Leszczynski (Springer Science+Business Media B.V, Berlin, Heidelberg, 2008), Vol. 5, pp. 125–156. 53T. A. Wesolowski, H. Chermette, and J. Weber, J. Chem. Phys. 105, 9182 (1996). 54N. Govind, Y. A. Wang, and E. A. Carter, J. Chem. Phys. 110, 7677 (1999). 55S. Östlund and S. Rommer, Phys. Rev. Lett. 75, 3537 (1995). 56S. Rommer and S. Östlund, Phys. Rev. B 55, 2164 (1997). 57S. Keller, M. Dolfi, M. Reiher, and M. Troyer, “QC-MAQUIS: A massively parallel, matrix-product-operator based quantum-chemical DMRG program,” (unpublished). 58C. R. Jacob, J. Neugebauer, and L. Visscher, J. Comput. Chem. 29, 1011 (2008). 59F. Aquilante, L. D. Vico, N. Ferré, G. Ghigo, P.-Å. Malmqvist, P. Neogrády, T. B. Pedersen, M. Pitonak, M. Reiher, B. Roos, L. Serrano-Andrés, M. Urban, V. Veryazov, and R. Lindh, J. Comput. Chem. 31, 224 (2010). 60G. te Velde, F. M. Bickelhaupt, E. J. Baerends, S. J. A. van Gisbergen, C. Fonseca Guerra, J. G. Snijders, and T. Ziegler, J. Comput. Chem. 22, 931 (2001). 61J. Neugebauer, C. R. Jacob, T. A. Wesolowski, and E. J. Baerends, J. Phys. Chem. A 109, 7805 (2005). 62M. Dolfi, B. Bauer, S. Keller, A. Kosenkov, T. Ewart, A. Kantian, T. Giamarchi, and M. Troyer, Comput. Phys. Commun. 185, 3430 (2014). 63C. R. Jacob, S. M. Beyhan, R. E. Bulo, A. S. P. Gomes, A. W. Götz, K. Kiewisch, J. Sikkema, and L. Visscher, J. Comput. Chem. 32, 2328 (2011). 64J. P. Perdew, J. A. Chevary, S. H. Vosko, K. A. Jackson, M. R. Pederson, D. J. Singh, and C. Fiolhais, Phys. Rev. B 46, 6671 (1992). 65E. Van Lenthe and E. J. Baerends, J. Comput. Chem. 24, 1142 (2003). 66A. Lembarki and H. Chermette, Phys. Rev. A 50, 5328 (1994). 67T. H. Dunning, Jr., J. Chem. Phys. 90, 1007 (1989). 68F. Neese, WIREs Comput. Mol. Sci. 2, 73 (2012). 69A. D. Becke, Phys. Rev. A 38, 3098 (1988). 70J. P. Perdew, Phys. Rev. B 33, 8822 (1986). 71S. Grimme, J. Antony, S. Ehrlich, and H. Krieg, J. Chem. Phys. 132, 154104 (2010). 72F. Weigend and R. Ahlrichs, Phys. Chem. Chem. Phys. 7, 3297 (2005). 73K. Eichkorn, O. Treutler, H. Öhm, M. Häser, and R. Ahlrichs, Chem. Phys. Lett. 240, 283 (1995).

Journal of Chemical Physics is copyrighted by AIP Publishing LLC (AIP). Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. For more information, see http://publishing.aip.org/authors/rights-and-permissions.

Self-consistent embedding of density-matrix renormalization group wavefunctions in a density functional environment.

We present the first implementation of a density matrix renormalization group algorithm embedded in an environment described by density functional the...
631KB Sizes 1 Downloads 5 Views