Periodic subsystem density-functional theory

Periodic Subsystem Density-Functional Theory Alessandro Genova,1 Davide Ceresoli,1,2 and Michele Pavanello1 arXiv:1406.7803v1 [cond-mat.mtrl-sci] 30 Jun 2014 1 Department of Chemistry, Rutgers University, Newark, NJ 07102, USA 2 CNR-ISTM: Institute of Molecular Sciences and Technologies, Milano, Italy Date: July 1, 2014 Status: Submitted to J. Chem. Phys. 1 Abstract By partitioning the electron density into subsystem contributions, the Frozen Density Embedding (FDE) formulation of subsystem DFT has recently emerged as a powerful tool for reducing the computational scaling of Kohn–Sham DFT. To date, however, FDE has been employed to molecular systems only. Periodic systems, such as metals, semiconductors, and other crystalline solids have been outside the applicability of FDE, mostly because of the lack of a periodic FDE implementation. To fill this gap, in this work we aim at extening FDE to treat subsystems of molecular and periodic character. This goal is achieved by a dual approach. On one side, the development of a theoretical framework for periodic subsystem DFT. On the other, the realization of the method into a parallel computer code. We find that periodic FDE is capable of reproducing total electron densities and (to a lesser extent) also interaction energies of molecular systems weakly interacting with metallic surfaces. In the pilot calculations considered, we find that FDE fails in those cases where there is appreciable density overlap between the subsystems. Conversely, we find FDE to be in semiquantitative agreement (but still within chemical accuracy) with Kohn–Sham DFT when the inter-subsystem density overlap is low. We also conclude that to make FDE a suitable method for describing molecular adsorption at surfaces, kinetic energy density functionals that go beyond the GGA level must be employed. 2 1 Introduction An important theoretical simplification to electronic structure theory was provided by the Hohemberg and Kohn theorems [1], according to which either the electron density, ρ(r), or the electronic wavefunction, can be considered as the central quantity in an electronic structure calculation. In other words, ρ(r) provides all information one needs to evaluate the electronic energy and observables of a molecular system. In practical calculations ρ(r) is computed with the Kohn–Sham DFT (KS-DFT) method [2], which has by now been implemented in a multitude of computer codes worldwide. KS-DFT is nowadays the most popular quantum mechanical method for calculating the electronic structure of molecules and materials. Because of the N 3 scaling with the sistem size, KS-DFT is limited to the type of sistems it can approach. Tipically, when theoretically modeling of a physicochemical process, a careful choice of the model system precedes the computations. A model system approachable by KS-DFT computations generally should not exceed 500 atoms, so that the simulations can be completed in a reasonable time. However, it is now understood that if chemical accuracy in the predictions is sought, the KS-DFT simulations must take into account the complexity of the environment surrounding the model system of interest [3–6]. Thus, the formulation of algorithms that simplify the KS-DFT problem by exploiting the locality of the electronic structure have been a popular avenue of research [7]. A natural way to simplify any problem is by applying the principle of divide and conquer. Such an approach is achieved in KS-DFT by partitioning the very basic quantity of DFT: the electron density. Namely, ρ(r) = NS X ρI (r). (1) I As we explain in more detail in the next section, such a partitioning results into a set of coupled differential equations which solve for each subsystem electron density, ρI . This is not a new idea. For example, the frozen density embedding (FDE) formulation of subsystem DFT introduced by Wesolowski and Warshel [8], has emerged as the leading approach to exploit Eq. (1). Since that first formulation, the FDE method has been implemented in several quantum chemistry packages that use atomic orbital basis sets, such as ADF [9] and TurboMole [10]. FDE has made possible the all-electron description of such large systems as pigment aggregates [11], bulk liquids [12], solvated 3 chromophores [13, 14], molecular crystals [15], as well as biological systems [16–18] to name a few applications. Another appealing feature of Eq. (1), is that it provides us with a molecule-centric description of nature. As many systems naturally occur as a set of weakly interacting molecules, a subsystem formulation of KS-DFT offers an intuitive description of the underlying physics and an easier interpretation of the chemistry as advocated by the earliest attempts [19, 20]. Such a molecule-centric perspective, for example, has motivated many other types of partitioning, e.g. the ones occurring at the wavefunction level [21, 22]. The possibility of applying FDE to periodic systems is of particular interest. Systems such as metals and semiconductors are involved in a plethora of applications important to chemistry and materials science, ranging from photovoltaics to catalysis. A major motivation of this work resides in the success of early implementations of Wavefunction in periodic DFT embedding schemes [23]. The CASTEP implementation [23], was employed in simulations involving two subsystems, of which one was treated at the correlated wavefunction level [24–26]. There is another periodic implementation of FDE in the CP2K code. This implementation, however, features a Brillouin zone sampling of the Γ-point only [12] and, to the best of our knowledge, it was only applied to the dynamics of liquid water. In this work, we build on the mentioned early work, and extend the FDE method to approach a system composed of multiple subsystems (2 or more) in which one or more subsystems can be periodic (e.g. metals or semiconductors). For this purpose, we introduce a novel implementation of FDE in the Q UANTUM ESPRESSO (QE) [27] suite of softwares. QE uses pseudopotentials and expands the so-called pseudowaves and the pseudodensity in a plane waves basis set. This attractive framework allows us to account for the non-vanishing Brillouin zone using a k-point mash so that we can model accurately the electronic structure of metal and semiconducting surface slab subsystems [28]. There are known weaknesses of the FDE method, which are related to the fact that FDE employs orbital–free Kinetic energy density functionals (KEDF). Unlike the exchange– correlation, semilocal KEDF approximants employed in a pure orbital–free framework are unable to reproduce basic aspects of the electronic structure of atoms, such as the shell structure [29]. In FDE, the KEDF are used only at the non-additive level (vide infra). In this role, they are shown to be successful provided that there is weak overlap between electron densities of the subsystems. Weakly interacting subsystems, such as hydrogen-bonded systems [30–35], van der Waals complexes [30,36–38], ionic bonds, and 4 even charge transfer systems [16, 17, 39] are shown to be successfully reproduced by FDE. Due to the lack of an extension of FDE to periodic systems (including k-point sampling), the ability of FDE to reproduce those interactions arising between molecular adsorbates and metallic or semiconducting surfaces has so far not been investigated. In this work, we make significant steps forward to fill this gap, and test a number of non-additive KEDF specifically for molecules at metal surfaces. This paper is organized as follows. The following section is completely devoted to the theoretical background of FDE when periodic systems are considered. Section 3 is devoted to the computational details of our calculations, while Section 4 presents numerical results of our FDE implementation applied to molecular and periodic systems. Section 5 collects conclusions and future directions of this work. 2 Theory of Periodic Frozen Density Embedding In this section, we will discuss the theory behind periodic subsystem DFT. As our goal is to tackle periodic systems, the equations derived for practical calculations (which we call Frozen Density Embedding, or FDE) are generalized here to a non-integer subsystem orbital occupation, and to a non-vanishing size of the Brillouin zone. Although we recognize a possible semantical inconsistency [40], hereafter, we will consider FDE a synonym of subsystem DFT. The foundamental idea behind any subsystem DFT method is that the electron density of the supersystem can be partitioned into the electron densities of the NS subsystems, as already shown in Eq. (1): NS X ρI (r) (2) ρtot (r) ≡ ρ(r) = I The subsystem densities are defined from the orbitals of the corresponding subsystem (for the sake of clarity only the closed-shell case is reported here): 2 ρI (r) = ΩBZ Z dk BZ X nIj,k uIj,k (r) 2 (3) j where ΩBZ is the volume of the BZ, and 0 ≤ nIj,k ≤ 1 are the occupation numbers of orbitals belonging to fragment I, and uIj,k (r) is such that the Bloch function is ψj,k = eik·r uIj,k (r). 5 The occupation numbers are formulated such that X nIj,k = NI ∀ k ∈ BZ, (4) j with NI the number of electrons assigned to subsystem I. The above definition of subsystem density differs from what has been considered in previous formulations of FDE. Here, besides an integral over the BZ, partial orbital occupations are invoked. There are two reasons for this. First a practical one, when tackling metallic systems, due to the continuous density of states at the Fermi energy, the SCF procedure would simply not converge without smearing the occupations across the Fermi energy. Second, when working in a finite temperature regime, the non-pure state resulting form the statistical occupation of excited states can be formally accounted for with a non-integer orbital occupation [41]. Since the orbitals of one subsystem are not required to be orthogonal to those in another subsystem, complications in the computation of the total non-interacting kinetic energy arise. We can still define a subsystem non-interacting (Janak) kinetic energy as: 2 TJ [ρI ] = ΩBZ Z dk BZ X nIj,k j * (∇ + ik)2 I − uj,k 2 uIj,k + (5) We are adopting this notation rather than Ts , because the latter is defined for integer occupations only [41]. As the supersystem total Janak kinetic energy is not simply the sum of the subsystem kinetic energies, it needs to be corrected. Namely, TJ [ρ] = X TJ [ρI ] + TJnad [{ρI }] (6) I The non-additive kinetic energy term is defined as: TJnad [{ρI }] = T̃J [ρ] − X T̃J [ρI ] (7) I where T̃J differs from TJ in the fact that it is a pure functional of the electron density. This approach is particularly useful in practical implementations of the method as it allows us to avoid the diagonalization step of the Fock matrix of the supersystem, effectively making FDE a linear scaling method. 6 We can also write other energy contributions in a similar fashion: F [ρ] = F [{ρI }] = X F [ρI ] + F̃ [ρ] − I X F̃ [ρI ] = I X F [ρI ] + F nad [{ρI }], (8) I where F [ρ] can be any functional of the total electron density, such as the Hartree energy EH or the exchange–correlation energy EXC . This trick allows us to write the total energy of the supersystem as the KS-DFT energies of the single subsystems, plus a contribution arising from the interaction with the other subsystems: EFDE = X TJ [ρI ] + I EeN [ρI ] + EH [ρI ] + EXC [ρI ] + I NS NS X X I J nad nad EeN [ρI ]+TJnad +EH +EXC +VN N J6=I (9) where is the electron–nuclear interaction energy due to the nuclei of subsystem J with the electrons of subsystem I, namely J EeN [ρI ] J EeN [ρI ] X Z Zα ρI (r) dr. =− |r − Rα | α∈J (10) In addition, VN N is the nuclear repulsion energy. We obtain the subsystem orbitals by solving self-consistently the following coupled equations:   1 2 I − (∇ + ik) + Veff (r) uIj,k (r) = ǫIj,k uIj,k (r). (11) 2 In the above equation, as we takle periodic systems the effective Hamiltonian for periodic systems is invoked which includes the k vector belonging to the BZ, solving for the I periodic part of the Bloch wave. The effective potential Veff is given by I I I Veff (r) = VeN (r) + VH [ρI ](r) + VXC [ρI ](r) + Vemb (r), (12) I and Vemb appearing above is the so-called embedding potential, which takes the following form # "Z NS ′ X X Zα ∂ T̃J [ρ] ∂ T̃J [ρI ] ∂ ẼXC [ρ] ∂ ẼXC [ρI ] ρJ (r ) ′ I dr − − + − . + Vemb (r) = ′ |r − r | |r − Rα | ∂ρ ∂ρI ∂ρ ∂ρI α∈J J6=I (13) 7 We refer to other review publications [42–44] regarding the more profound theoretical ramifications of using the partitioning in Eq. (2), and the assumptions related to the non-interacting v-representability of each subsystem density. We will limit ourselves in stating that the more accurate the non-additive KEDF is the closer the FDE results converge to KS-DFT of the supersystem [43]. As the exact functional is unknown, the most important avenue of research in FDE is the systematical improvement of the available KEDF approximants. [Figure 1 about here.] We should remark that the plane wave basis set provides a very efficient mean of calculating the electron–electron repulsion through the solution of the Poisson equation for the corresponding potential in Fourier space. By exploiting the fast Fourier transform (FFT), this problem is solved with an algorithm that scales as O(N log N ) with N being the number of plane waves used to span the charge density. This is in contrast to localized orbital-based codes, where the Hartree potential is computed in real-space. This computation inherently scales non-linearly with the size of the atomic orbital basis set. The so-called “freeze-and-thaw” procedure [40] by-passes this problem by evaluating the Hartree potential due to the frozen density only at the beginning of the calculation. This is made even more computationally amenable by employing density fitting techniques [9]. We have, therefore, re-casted the FDE equations in such a way that the Hartree potential due to the total electron density can be calculated as often as each SCF cycle. This is achieved by imposing that the SCF cycle in Figure 1b is executed only once. In the implementation, this is achieved using Message-Passing Interface (MPI) as indicated in Figure 1a. For those interested in the details of the working equations, we include an appendix at the end of the paper. 3 Computational Details All the calculations (KS-DFT and FDE) have been performed with Q UANTUM ESPRESSO. We have used Ultrasoft pesudopotentials from the original QE library and from GBRV [45] library. The plane wave cutoffs are 40 Ry and 400 Ry, for the wave functions and density, respectively. Unless stated otherwise, the PBE [46] functional has been employed for the exchange-correlation, and either the LC94 [47] or revAPBEK [48] functionals have been employed for the non-additive KEDF. 8 Regarding sampling of the BZ, in the calculations invoving molecular systems, we have only sampled the Γ point, while for periodic systems we have used the Monkhorst– Pack [28] sampling with a 2 × 2 × 1 mesh. In addition, we used a Methfessel–Paxton smearing [49] with 1 × 10−2 Ry of smearing parameter. The smearing was only applied to the metallic subsystems, whereas the molecular subsystems were treated with integer orbital occupations. We have benchmarked FDE calculations against KS-DFT references by comparing physisorption energies and corresponding electron densities obtained with the two methods. The simulations consist of single point calculations with the same geometry for both KS-DFT and FDE (e.g. the equilibrium geometry of the system calculated with KS-DFT with exception of the 2PDI on Au). Moreover, a more insightful comparison is made by calculating the number of electrons misplaced by FDE, h∆ρi, defined as: 1 h∆ρi = 2 Z |∆ρ(r)| dr (14) where ∆ρ(r) = ρFDE (r) − ρKS (r) (15) h∆ρi is an important quantity, as it vanishes only when FDE and KS-DFT electron densities coincide. As h∆ρi is a size-sensitive quantity, we always compare it to the total number of subsystems and the total number of electrons. 4 Results We have chosen a set of molecular systems (water dimer and ammonia borane) as well as molecules on metal surfaces (CH4 on Pt, water on Pt and PDI on Au) as our set of pilot calculations. The goal of the calculations involving the molecular systems is to show that our FDE implementation reproduces the already reported behavior of FDE for these systems. Generally, FDE delivers results close to KS-DFT when the intersubsystem density overlap is low. Thus, water dimer is expected to be well described by FDE, while for the ammonia borane system, due to the partial charge transfer and covalent character of the B–N dative bond, FDE is expected to fail. We will show that when tackling molecules adsorbed on surfaces, similar considerations to the molecular case apply. I.e. the larger the density in the region between subsystems, the less accurate FDE will be. 9 [Table 1 about here.] 4.1 H2 O dimer The first system we used to benchmark our FDE code is the water dimer, with the two water molecules arranged to form a single hydrogen bond, see Figure 2a. It has been shown in previous publications [12, 32] that localized basis set FDE is capable of obtaining very accurate results for this system when a GGA KEDF is employed for the non-additive kinetic energy part. [Figure 2 about here.] We have calculated an attractive interaction energy of 4.317 kcal with regular KS-DFT, while FDE just slightly overestimates this energy by about 0.3 kcal and in this particular case is well within what’s considered chemical accuracy (1 kcal). The very good agreement between between KS-DFT and FDE for the interaction energy is also reflected in the number of misplaced electrons h∆ρi. As reported in Table 1 only 2.85 × 10−2 electrons out of 16 have been displaced by the selfconsistent FDE calculation. In Figure 2a the difference between the FDE and KS-DFT densities, ∆ρ(r), is plotted as defined in Eq. (15). The figure shows that the FDE method localizes too much density on the oxygen lone pairs in turn depleting of electron density the hydrogen bond between the two dimers. It has to be pointed out that the 1 × 10−3 isosurface is a very strict threshold. As a reference, Figure 2b reports the difference between the KS-DFT density of the dimer and the sum of the density of the two isolated fragments. Using the same isosurface we can see that the difference is much larger than that of FDE. Thus, we conclude that for this system FDE performs very well, recovering almost exactly the same electronic structure obtained with a KS calculation. 4.2 Ammonia Borane As shown previously [50], the ammonia borane complex is expected to represent a challenge for FDE. It is a Lewis acid-base complex that exhibits a partial charge transfer from the ammonia lone pair to the borane, yielding a so-called dative bond characterized by a dissociation energy of 34.783 kcal/mol as calculated by us using KS-DFT, in good agreement with coupled cluster calculations [51]. The FDE interaction energy of 10 80.865 kcal/mol and 82.468 kcal/mol significantly overestimates the reference when using LC94 or revAPBEK non-additive KEDF, respectively. More important is the analysis of the electron density: the difference between the FDE density (and superposition of the isolated fragments densities) with the supersystem KS density is reported in Figure 3a (Figure 3b). As the electron density undergoes significant changes upon formation of the complex, FDE only qualitatively reproduces the KS-DFT density — there is a general underestimation of electron density along the B-N bond, while the hydrogens bonded to N (B) are too acidic (basic). Similar results are obtained from localized basis formulations of FDE [50]. [Figure 3 about here.] Quantitatively, the number of electrons misplaced by FDE in our calculations is 1.940 × 10−1 and 1.981 × 10−1 when employing the non-additive KEDF LC94 and revAPBEK, respectively. As expected, this is not satisfactory considering that the supersystem only has 14 electrons. In conclusion, we confirm previous FDE calculations of this complex [50] and we also show that our plane wave implementation maintains quantitatively the FDE electronic structure predicted by atomic orbital based implementations. 4.3 CH4 on Pt(100) Catalitic activation of methane by adsorption on a transition metal surface represents an intresting problem, as it could represent a cheap way to produce hydrogen, or a precursor in the synthesis of more complex molecules. On a Pt(100 ) surface methane can be adsorbet on several sites, so-called top, bridge, and hollow. In addition the molecule can have several configurations with respect to the surface. We have tested the behavior of FDE against the 8 configurations shown in figure 4. As it has been obeserved in Ref. [52] and confirmed in our calculations (Table 1) the CH4 −Pt interaction is generally very weak, and the potential energy surface is relatively uncorrugated (i.e. insensitive to the methane orientation). The only clear potential well is represented by the T1 configurations, which has an interaction energy about 0.5 kcal/mol higher than all the others. From an analysis of the local density of states of this configuration, it has been concluded that in the T1 case we are dealing with a real chemisorption (electron backdonation from the the metal d orbitals to the carbon), rather than a physisorption (polarization of the metal-methan electron densities). 11 We can see that FDE reproduces the supersystem KS calculation in a very accurate way for all configurations but the T1: interaction energies are always within chemical accuracy, and the number of displaced electrons is very small, ranging from 2.1 × 10−2 to 3.2 × 10−2 for the LC94 functional, and from 1.5 × 10−2 to 2.5 × 10−2 for the revAPBEK functional. It has to be remarked that the revAPBEK functional provided the best results also for magnitude of the interaction energy. As expected the T1 configuration is the one for which FDE yielded the worst results, because of the chemical nature of the bond between the hydrocarbon and the surface: we can see that FDE predicts a small repulsive interaction, equal to 1.302 kcal/mol for LC94 and 0.987 kcal/mol for revAPBEK, whereas KS has an attractive energy of 1.823 kcal/mol. The number of electrons misplaced ∆ρ is about 9 × 10−2 , which is which is much higher than what abserved for all the other configurations, but still reasonably small. In Figure 5, an isosurface plot of ∆ρ(r) for the T1 configuration is reported. Again, we see that the selfconsistent FDE density is smaller than the KS-DFT one in the interfragment bonding region. [Figure 4 about here.] [Figure 5 about here.] 4.4 H2 O on Pt(111) The most complex system we have used to benchmark this periodic implementation of FDE is a water bilayer adsorbed on a Pt(111 ) surface. The system consists of 13 fragments, 12 water molecules plus the metal surface. The geometry we have employed to model the water bilayer is the so-called RT3, as described in Refs. [53, 54] On the layer closer to the surface, we can distinguish between two kinds of water molecules: the ones with the dipole moment parallel to the metal surface, and the ones with a non-zero dipole component in the perpendicularly to the surface. The intermolecular interactions playing an important role in keeping the the bilayer and the metal together are three: (1) the hydrogen bonds network, (2) the interaction between the surface and the parallel dipole water molecules, and (3) the interaction between the surface and the water molecules having non-zero perpendicular component of the dipole. We have already seen in Section 4.1 that FDE describes well hydrogen bonds. This system, however, will probe the ability of FDE to describe the other two cases. Thus, we have run two additional simulations: 12 1. A single parallel water molecule adsorbed on the surface; 2. A single non-parallel water molecule adsorbed on the surface; In both cases the position of the molecule is the same as that in the bilayer, without performing any additional optimization. We will describe these two additional calculations and the one pertaining the full bi-layer separately below. 4.4.1 Parallel H2 O For the parallel water molecule, KS predicts an adsorption energy of 3.168 kcal while according to FDE the molecule is unbound, with and adsorption energy very close to zero, see Table 1. This discrepancy is also reflected in the number of misplaced electrons ∆ρ, that for this system is the highest so far and equal to 1.643 × 10−1 . Plots of ∆ρ(r) in Figure 6a show that once again FDE places too much density on the oxygen lone pairs and not enough in the bonding region between the molecule and the surface. 4.4.2 Perpendicular H2 O The vertical water molecules are further from the surface, and interact with it very weakly, 0.523 kcal as predicted by KS, and 0.162 kcal according to FDE. Being the density of the two fragments almost non-overlapping, obviously FDE has no problems recovering the same total density as KS, with only 1.62 × 10−2 electrons displaced. 4.4.3 Water Bilayer As we have seen in the previous sections, we expect that for the bilayer the largest source of error is due to the water molecules in the first layer parallel to the surface and to a lesser extent due to the hydrogen bond network. [Figure 6 about here.] The total intermolecular interactions as predicted by KS-DFT is 135.147 kcal, while FDE underestimates it by just 10% yielding an interaction equal to 122.594 kcal. As shown in Table 1 only 6.487 × 10−1 out of the total 368 are misplaced by FDE in this complex 13 fragments system. An isosurface plot of ∆ρ(r) is presented in Figure 6b. In the figure we see that the expected maximum discrepancy between the KS-DFT and FDE densities is 13 in the bond between the parallel water molecules and the surface, as well as in the many hydrogen bonds keeping the bilayer together. 4.5 Perylene Diimide on Au(111) Perylene diimides (PDIs) are a class of compounds whose electron-accepting character and charge transport properties make them useful materials in the design of new photovoltaic cells based on organic compounds. We have run a single point KS simulation of two N,N’-dihydro perylene diimides stacked on top of a Au surface. The geometry has been optimized at the KS-DFT level constraining the PDI–PDI distance and the PDI–Au distance to be 3 Å. From the KS-DFT calculation of the supersystem, it has been predicted a 8.678 kcal/mol repulsive interaction energy. The positive sign in the interaction energy is not surprising since van der Waals (vdW) forces are the most important intermolecular forces in this case, and it is known that KS-DFT calculations employing semilocal GGA exchange– correlation functionals is generally unable to recover them. Conversely, FDE predicts an attractive interaction energy of 17.078 kcal/mol and 16.064 kcal/mol with LC94 and revAPBEK, respectively. Despite not reproducing the KS-DFT interaction energy, this result is correctly attractive. Such a behavior from FDE has been reported previously [30, 55]. As reported in Table 1, the number of displaced electrons is 3.4 × 10−1 with LC94 and 3.9 × 10−1 with revAPBEK – a negligible value considering that each PDI molecule has 140 electrons. 5 Conclusions In this work, we have extended the Frozen Density Embedding formulation of subsystem DFT to treat inherently periodic subsystems. In FDE, the subsystems are treated as isolated Kohn–Sham systems which interact through an effective embedding potential. This potential is a functional of the electron density of all subsystems. We introduce periodic boundary conditions by applying the Bloch theorem to the Kohn–Sham orbitals of the subsystems. Thus, sampling over the Brillouin zone (the so-called k-point sampling) enters the working Kohn–Sham-like equations needed to recover the electronic structure (electron density) of each subsystem. 14 The ultimate goal of this approach is to reduce the computational complexity of KSDFT for systems composed of molecules interacting with surfaces of metals and semiconductors. Several open questions have been answered here. When semilocal GGA kinetic energy functionals are employed, our implementation of FDE reproduces the results of existing molecular implementations (i.e. implementations employing a localized orbital basis) for interacting molecules. We have established this connection for those systems which are well described by FDE and for those which semilocal FDE is known to fail (such as the ammonia borane complex). We have also run pilot calculations to establish how well FDE tackles weakly interacting molecules at metal surfaces. We have established that FDE yields chemical accuracy when, in the KS-DFT calculation, there is no hybridization of the orbitals of the molecules with the bands of the metal. With such encouraging results, this work sets the stage to applying FDE to a wide array of molecular systems interacting with surfaces with applications to energy-related materials as well as catalysis and photovoltaics. Because of the inapplicability of semilocal FDE to molecule–surface interactions beyond weak physisorption, we also conclude that non-locality in the energy functionals (both at the kinetic and the exchange–correlation level) must be included if FDE aspires to become a true alternative to Kohn–Sham DFT. Efforts in this directions have been taken in our lab regarding the exchange–correlation [55], and ongoing efforts are placed in the search for computationally amenable ways to include non-locality at the kinetic energy level as well. 6 Acknowledgements We thank Dr. Oliviero Andreussi for stimulating discussions. We acknowledge startup funds provided by Rutgers University–Newark, for support of this research. 7 Appendix 7.1 Recast of the Effective Potential expression Due to the fact that in molecular codes employing localized orbital basis sets the computation of two-electron integrals can be expensive (Gaussian basis), or only possible through density fitting (Slater basis), the Hartree energy and potential in the subsystem 15 DFT framework is usually computed as: N S 1X EH [ρ] = 2 I Z VH [ρI ](r)ρI (r) dr + Z VHIemb [ρJ6=I ](r)ρI (r) dr  , (16) where ρI (r′ ) dr′ |r − r′ | NS Z NS X ρJ (r′ ) ′ X Iemb VH [ρJ6=I ](r) = VH [ρJ ](r) ≡ V [ρ](r) − V [ρI ](r), dr ≡ ′| |r − r J6=I J6=I VH [ρI ](r) = Z (17) (18) from which is clear that we are indeed calculating the total Hartree energy in a similar way as the kinetic energy: EH [ρ] = NS X I EH [ρI ] + EH [ρ] − NS X EH [ρI ] ≡ I NS X nad EH [ρI ] + EH . (19) I This approach is useful for the freeze-and-thaw scheme implemented in molecular codes. In this scheme, VHIemb (r) is calculated once and it doesn’t change until the end of the SCF cycle of the active fragment. The plane wave basis set allows us to efficiently calculate the potential due to the total electron density, VH [ρ]. Thus, we have the possibility to do so on the fly at each SCF step (and we almost always do). In this way, there is no need to calculate the embedding potential in Eq. (13), and the effective potential appearing in the KS-like equations, Eq. (11), becomes: I nad nad Veff [ρ, ρI ](r) = VeN (r) + VH [ρ](r) + VXC [ρI ](r) + Vkin [ρ, ρI ](r) + VXC [ρ, ρI ](r) (20) where in VeN (r) we have grouped the interactions with all nuclei in the system. 7.2 Recast of the Total energy expression To better exploit the Q UANTUM ESPRESSO code base, we found useful to express the FDE energy as a sum of band energies (i.e. the sum of eigenvalues of the occupied states), and a sum of the energy of the single fragments. 16 For a single fragment the band energy is: X k 1 I [ρ, ρI ]|ψkI i hψkI | − ∇2 + Veff 2 k Z Z = TJ [ρI ] + VeN (r) ρI (r) dr + VH [ρ](r) ρI (r) dr Z Z + VXC [ρI ](r) ρI (r) dr + ṼXC [ρ](r) ρI (r) dr Z Z − ṼXC [ρI ](r) ρI (r) dr + Ṽkin [ρ](r) ρI (r) dr Z − Ṽkin [ρI ](r) ρI (r) dr ǫIk = X (21) If we define the quantity VK ρK (found in the code as the deband variable) as V K ρK = − Z   VH [ρ] + VXC [ρK ] + ṼXC [ρ] − ṼXC [ρK ] + Ṽkin [ρ] − Ṽkin [ρK ] ρK (r) dr (22) we see from Eq. (21) that: nK X ǫiK + VK ρK = Ts [ρK ] + EeN [ρK ] (23) i And defining the fragment energy as: EK = nK X ǫiK + VK ρK + EXC [ρK ] + EH [ρK ] + VN N i (24) ≡ Ts [ρK ] + EeN [ρK ] + EXC [ρK ] + EH [ρK ] + VN N we finally obtain that the FDE energy is simply: EFDE = N X nad EK + Tsnad + EXC − (N − 1)VN N K 17 (25) 7.3 Kinetic energy functional and potential We have implemented an array of semilocal kinetic energy functional in our code. The general form of these functionals is: Tsnad [ρ, ∇ρ] = CT F Z dr ρ5/3 (r)F (s(r)) (26) where s(r) = Cs |∇ρ(r)| ρ4/3 (r) (27) 3 (3π 2 )2/3 10 1 Cs = 2(2π 2 )1/3 CT F = (28) (29) For a generic functional F [ρ, ∇ρ] = Z dr G(ρ(r), ∇ρ(r)) (30) its functional derivative is given by ∂G δF = δρ(r) ∂ρ − ∇· ρ=ρ(r) ∂G ∂∇ρ (31) ρ=ρ(r) We can therefore calculate the kinetic energy potential for any of these functionals   ∂G 5 2/3 4 ∂F 5/3 ∂F ∂s 2/3 5 = ρ F (s) + ρ =ρ F (s) − s(r) ∂ρ 3 ∂s ∂ρ 3 3 ∂s ∂F ∂s ∇ρ ∂F ∇ρ ∂G = ρ5/3 = Cs ρ1/3 ∂∇ρ ∂s ∂|∇ρ| |∇ρ| ∂s |∇ρ| (32) (33) Therefore ∂Tsnad [ρ, ∇ρ] = CT F ∂ρ  2/3 ρ     5 4 ∂F 1/3 ∂F ∇ρ F (s) − s(r) − ∇ · Cs ρ 3 3 ∂s ∂s |∇ρ| ≡ Ṽkin [ρ](r) 18 (34) 7.4 Handling of non-local pseudopotentials In all Q UANTUM ESPRESSO calculations carried out in this work, either Ultrasoft (US) [56] or Projector Augmented Wave (PAW) [57] pseudopotentials are employed. In the PAW formalism the true all-electron single particle wavefunction ψn can be recovered from a smooth auxiliary wavefunction ψ̃n through the transformation |ψn i = T̂ |ψ̃n i (35) where the short label n identifies a k, band, and spin index. The T̂ operator can be written as: T̂ = 1 + XX a |φai i − |φ̃ai i i  hp̃ai | (36) where a labels the augmentation sphere of the ions (i.e. the nuclei), the so-called partial waves |φai i form a complete basis set, and the auxiliary |φ̃ai i obey |φai i = T̂ |φ̃ai i, and the set of the smooth projectior functions hp̃ai | has the following properties: hp̃ai |φ̃aj i = δij X |φ̃ai i hp̃ai | = 1 (37) (38) i p̃ai (r) ≡ hr|p̃ai i for |r − Ra | > rca = 0, (39) In the frozen core approximation of PAW the core electron are considered as frozen, and the eigen value problem is restricted only to the valence electrons. The full valence wave function obtained from Eq. (35) is guaranteed to be orthogonal with respect to the core states. Since in the FDE formalism the KS-DFT orbitals of one subsystem are not required to be orhogonal to those (valence nor core) in another subsystem (the non-additive kinetic energy potential is responsible of recovering those properties that in the KS method are guaranteed from orthonormality), the operator T̂ describing the core states of one subsystem does not need to include projectors belonging to other subsystems. Thus, a subsystem-specific PAW projector must be used. Namely, I T̂ = 1 + XX a∈I |φai i i 19 − |φ̃ai i  hp̃ai | . (40) References [1] P. Hohenberg and W. Kohn, Phys. Rev. 136, 864 (1964). [2] W. Kohn and L. J. Sham, Phys. Rev. 140, 1133 (1965). [3] D. Hartmann Douma, B. M’Passi-Mabiala, and R. Gebauer, J. Chem. Phys. 137, 154314 (2012). [4] F. De Angelis, S. Fantacci, and R. Gebauer, J. Chem. Phys. Lett. 2, 813 (2011). [5] B. Mennucci, Phys. Chem. Chem. Phys. 15, 6583 (2013). [6] J. Neugebauer, C. Curutchet, A. Munioz-Losa, and B. Mennucci, J. Chem. Theory Comput. 6, 1843 (2010). [7] D. R. Bowler and T. Miyazaki, Rep. Prog. Phys. 75, 036503 (2012). [8] T. A. Wesolowski and A. Warshel, J. Phys. Chem. 97, 8050 (1993). [9] C. R. Jacob, J. Neugebauer, and L. Visscher, J. Comput. Chem. 29, 1011 (2008). [10] S. Laricchia, E. Fabiano, and F. Della Sala, J. Chem. Phys. 133, 164111 (2010). [11] C. König and J. Neugebauer, Phys. Chem. Chem. Phys. 13, 10475 (2011). [12] M. Iannuzzi, B. Kirchner, and J. Hutter, Chem. Phys. Lett. 421, 16 (2006). [13] J. Neugebauer, M. J. Louwerse, P. Belanzoni, T. A. Wesolowski, and E. J. Baerends, J. Chem. Phys. 123, 114101 (2005). [14] J. Neugebauer, M. J. Louwerse, E. J. Baerends, and T. A. Wesolowski, J. Chem. Phys. 122, 094115 (2005). [15] R. Kevorkyants, X. Wang, D. M. Close, and M. Pavanello, J. Phys. Chem. B 117, 13967 (2013). [16] P. Ramos and M. Pavanello, J. Chem. Theory Comput. 10, 2546 (2014). [17] A. Solovyeva, M. Pavanello, and J. Neugebauer, J. Chem. Phys. 136, 194104 (2012). [18] C. R. Jacob and L. Visscher, J. Chem. Phys. 128, 155102 (2008). 20 [19] R. G. Gordon and Y. S. Kim, J. Chem. Phys. 56, 3122 (1972). [20] W. Kolos and E. Radzio, Int. J. Quantum Chem. 13, 627 (1978). [21] B. Jeziorski, R. Moszynski, and K. Szalewicz, Chem. Rev. 94, 1887 (1994). [22] L. Rajchel, P. S. Zuchowski, M. M. Szczceniak, and G. Chałasinski, Phys. Rev. Lett. 104, 163001 (2010). [23] N. Govind, Y. Wang, A. da Silva, and E. Carter, Chem. Phys. Lett. 295, 129 (1998). [24] P. Huang and E. A. Carter, J. Chem. Phys. 125, 084102 (2006). [25] T. Klüner, N. Govind, Y. A. Wang, and E. A. Carter, J. Chem. Phys. 116, 42 (2002). [26] T. Klüner, N. Govind, Y. A. Wang, and E. A. Carter, Phys. Rev. Lett. 86, 5954 (2001). [27] P. Giannozzi et al., J. Phys.: Condens. Matter 21, 395502 (2009). [28] H. J. Monkhorst and J. D. Pack, Phys. Rev. B 13, 5188 (1976). [29] Y. A. Wang and E. A. Carter, Orbital-Free Kinetic-Energy Density Functional Theory, in Theoretical Methods in Condensed Phase Chemistry, edited by S. D. Schwartz, pages 117–184, Kluwer, Dordrecht, 2000. [30] R. Kevorkyants, M. Dulak, and T. A. Wesolowski, J. Chem. Phys. 124, 024104 (2006). [31] T. A. Wesolowski, H. Chermette, and J. Weber, J. Chem. Phys. 105, 9182 (1996). [32] T. A. Wesolowski, J. Chem. Phys. 106, 8516 (1997). [33] K. Kiewisch, G. Eickerling, M. Reiher, and J. Neugebauer, J. Chem. Phys. 128, 044114 (2008). [34] A. W. Götz, S. M. Beyhan, and L. Visscher, J. Chem. Theory Comput. 5, 3161 (2009). [35] T. A. Wesolowski, J. Am. Chem. Soc. 126, 11444 (2004), PMID: 15366883. [36] T. A. Wesolowski, Y. Ellinger, and J. Weber, J. Chem. Phys. 108, 6078 (1998). [37] T. A. Wesolowski, P.-Y. Morgantini, and J. Weber, J. Chem. Phys. 116, 6411 (2002). [38] C. R. Jacob, T. A. Wesolowski, and L. Visscher, J. Chem. Phys. 123, 174104 (2005). 21 [39] A. Solovyeva, M. Pavanello, and J. Neugebauer, J. Chem. Phys. 140, 164103 (2014). [40] T. A. Wesolowski and J. Weber, Chem. Phys. Lett. 248, 71 (1996). [41] R. G. Parr and W. Yang, Density-Functional Theory of Atoms and Molecules, Oxford University Press, Oxford, 1989. [42] T. A. Wesolowski, One-Electron Equations for Embedded Electron Density: Challenge for Theory and Practical Payoffs in Multi-Level Modeling of Complex Polyatomic Systems, in Computational Chemistry: Reviews of Current Trends, edited by J. Leszczynski, volume 10, pages 1–82, World Scientific, Singapore, 2006. [43] C. R. Jacob and J. Neugebauer, WIREs: Comput. Mol. Sci. , Accepted for publication (2013). [44] O. V. Gritsenko, On the Principal Difference Between the Exact and Approximate Frozen-Density Embedding Theory, in Recent Advances in Orbital-Free Density Functional Theory, edited by T. A. Wesolowski and Y. A. Wang, chapter 12, pages 355–365, World Scientific, Singapore, 2013. [45] K. F. Garrity, J. W. Bennett, K. M. Rabe, and D. Vanderbilt, Computational Materials Science 81, 446 (2014). [46] J. P. Perdew, K. Burke, and M. Ernzerhof, Phys. Rev. Lett. 78, 1396 (1997). [47] A. Lembarki and H. Chermette, Phys. Rev. A 50, 5328 (1994). [48] L. A. Constantin, E. Fabiano, S. Laricchia, and F. Della Sala, Phys. Rev. Lett. 106, 186406 (2011). [49] M. Methfessel and A. T. Paxton, Phys. Rev. B 40, 3616 (1989). [50] S. Fux, K. Kiewisch, C. R. Jacob, J. Neugebauer, and M. Reiher, Chem. Phys. Lett. 461, 353 (2008). [51] C. W. Bauschlicher Jr. and A. Ricca, Chem. Phys. Lett. 237, 14 (1995). [52] P. Moussounda, M. Haroun, B. M’Passi-Mabiala, and P. Legare, Surface Science 594, 231 (2005). [53] D. L. Doering and T. E. Madey, Surface Science 123, 305 (1982). 22 [54] S. Meng, E. G. Wang, and S. Gao, Phys. Rev. B 69, 195404 (2004). [55] R. Kevorkyants, H. Eshuis, and M. Pavanello, J. Chem. Phys. (2014), submitted. [56] D. Vanderbilt, Phys. Rev. B 41, 7892 (1990). [57] P. E. Blöchl, Phys. Rev. B 50, 17953 (1994). 23 List of Figures 1 2 3 4 6 Parallel architecture Water dimer . . . . Ammonia borane . CH4 Conformations Water on Pt . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 24 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 25 26 27 28 30 MPI_COMM_WORLD MPI_Comm_split ... Frag. #1 ρ1(r) Frag. #2 ρ2(r) intra_image_comm inter_image_comm Frag. #Ns-1 Frag. #Ns ρNs-1(r) ρNs(r) ... MPI_Allreduce ρ(r) (a) (b) Figure 1: (a): The structure of the MPI routines. The MPI routines are responsible for “summing” all the subsystem densities, {ρI }, I = 1, 2 . . . NS and broadcasting the toal density, ρ, to all processes. Subsystems communicate through the inter-image communicator, while the various processes solving for the single subsystem electronic structure communicate with the intra-image communicator. (b): The subsystem SCF cycle. “FDE equations” stands for Eq. (11). If the SCF is executed to selfconsistency before broadcasting the {ρI }, the freeze-and-thaw preocedure is recovered. In our code, we allow the total density to be updated at every SCF cycle so that when the SCF has converged all subsystem densities are computed selfconsistently to each other. 25 (a) (b) Figure 2: Water dimer. (a): Difference between the FDE density and the KS density; (b): Difference between the sum of the densities of isolated fragments and the KS density. In both cases a 1 × 10−3 isosurface is plotted. 26 (a) (b) Figure 3: Ammonia borane complex. Borane on top. (a): Difference between the FDE density and the KS density; (b): Difference between the sum of the densities of isolated fragments and the KS density. In both cases a 1 × 10−3 isosurface is plotted. 27 (a) T1 (b) T2 (c) T3 (d) B2a (e) B2p (f) H1 (g) H2 Figure 4: Configurations of methane adsorbed on Pt(100) considered in the calculations. 28 (h) H3 Figure 5: Methane on Pt ∆ρ(r) 0.001 isosurface plot. Yellow identifies positive regions (excess of FDE density), blue identifies negative regions (excess of KS density). 29 (a) (b) Figure 6: Water on Pt (a) and water bilayer on Pt (b) ∆ρ(r) 0.001 isosurface plot. Yellow identifies positive regions (excess of FDE density), blue identifies negative regions (excess of KS density). 30 List of Tables 1 Results summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 32 31 Table 1: Results summary Table. For each system studied, we have calculated the interaction energy, and the number of electrons that has been displaced by FDE with respect to the KS reference Int. En. KS Int. En. LC94 h∆ρi LC94 (kcal/mol) (kcal/mol) Int. En. rAPBEK h∆ρi rAPBEK (kcal/mol) H2 O Dimer BH3 −NH3 −4.317 −34.783 −4.645 −80.865 0.0285 0.1940 −4.472 −82.468 0.0291 0.1981 CH4 on Pt T1 CH4 on Pt T2 CH4 on Pt T3 CH4 on Pt B2a CH4 on Pt B2p CH4 on Pt H1 CH4 on Pt H2 CH4 on Pt H3 −1.823 −1.375 −1.368 −1.366 −1.313 −1.271 −1.306 −1.346 +1.302 −0.492 −0.493 −0.433 −0.606 −0.546 −0.540 −0.535 0.0917 0.0291 0.0235 0.0282 0.0280 0.0318 0.0262 0.0213 +0.987 −0.690 −0.687 −0.601 −0.712 −0.726 −0.651 −0.640 0.0903 0.0248 0.0180 0.0238 0.0231 0.0285 0.0215 0.0155 H2 O on Pt (h) H2 O on Pt (v) 12 H2 O on Pt −3.168 −0.523 −135.470 +0.040 −0.162 −122.594 0.1643 0.0162 0.6487 +0.074 +0.046 −117.824 0.1634 0.0135 0.6421 +8.678 −17.078 0.3407 −16.064 0.3877 2 PDI on Au 32