Dispersion-Free Component of Non-Covalent Interaction via Mutual Polarization of Fragments ’ Densities

Fragments’ Densities Marcin Modrzejewski,1, a) Łukasz Rajchel,2, b) Małgorzata M. Szczęśniak,3 and Grzegorz Chałasiński1 1)Faculty of Chemistry, University of Warsaw, 02-093 Warsaw, Pasteura 1, Poland 2)Interdisciplinary Centre for Mathematical and Computational Modelling, University of Warsaw, 02-093 Warsaw, Pawińskiego 5a, Poland 3)Department of Chemistry, Oakland University, Rochester, Michigan 48309-4477, USA


I. INTRODUCTION
The density functional theory (DFT) has emerged as a powerful and efficient methodology for studying chemical systems since Hohenberg and Kohn 1 laid out the milestone founding theorems of the theory and came up with the ingenious idea on how to put them to work, back in the 1960s. 2 The plethora of approximations to the exchange-correlation (xc) functional, the crucial ingredient of the DFT, has been proposed and thoroughly benchmarked (as a recent example, see Ref. 3

and the Refs. therein).
On the one hand, the density functional approximations perform satisfactorily for the determination of geometries, energetics, and static properties of the chemical systems, as well as the interaction energies of hydrogen bonds and electrostatically-bound complexes. 4,5 the other hand, the existing xc functionals fail spectacularly to describe the dispersion interactions. 6,7The local (in case of local spin-density approximations), or semi-local (for generalized gradient approximations) character of current xc functionals makes them unable to describe the long-range intermonomer correlation effects manifesting themselves as the dispersion interaction. 8,91][12] Other approaches for solving the dispersion problem within DFT are generally attempts to model the dispersion energy separately.The most successful examples are dispersion-corrected atom-centered potentials, 13 DFT+dispersion method, 14,15 exchange hole dipole model, 16 and DFT-based symmetry adapted perturbation theory (SAPT). 17,18See Ref. 19 for a recent comprehensive review of the DFT-based dispersion approaches.
For DFT-based dispersion models to be successful, the proper description of the dispersionfree interaction energy is necessary.In the above-mentioned approaches (with the exception of SAPT) dispersion-free energy is obtained in a supermolecular fashion.Such an energy, however, already includes some fraction of correlation effects which can become doubly counted upon combining with dispersion correction.Moreover, the exchange contribution must be compatible with dispersion correction in the regions important for the interaction energy. 20These issues have motivated a few groups to design the DFT method that would rigorously exclude the dispersion interaction.The examples are: • The dispersionless functional of Pernal et al. 21The xc functional is constructed through the fitting to the set of interaction energies from which the dispersion energy is subtracted.
• The Pauli blockade (PB) method described in Refs.22,23.The aims of the method are to include the intramonomer correlation, to rigorously exclude the intermonomer correlation.
The former goal is accomplished by the KS description of the interacting monomers while the latter is achieved by the use of the HF exchange between the monomers.
Consequently, the PB energy is dispersion-free.
• A long-range correction scheme for the exchange functional [24][25][26][27][28] which was proven to correct the inaccuracies of local exchange potentials in regions of large reduced density gradients, important for intermolecular interactions. 20Such calculations lead to energies which can be supplied with a proper dispersion contribution, as was first shown by Kamiya, Tsuneda, and Hirao 20 (see also Ref. 29).
Below we will focus entirely on the PB approach.Our previous works on the subject did not include any calculations performed within a systematic set of noncovalently bound systems.For this reason, no definite conclusions could be made to date on the general performance of the PB method.In this work, we perform calculations for the well-established test suite of noncovalently bound systems, 30,31 which will help identify systems for which the PB method performs satisfactorily, and to find the systems for which it should not be applied.Additionally, an efficient algorithm for converging main equation of the PB method is presented.The new algorithm corrects divergent and oscillatory behavior of the previously used penalty function method. 32
All matrices are written in the OAO basis, except in cases where the atomic orbital (AO) label is present.
Consider a system of two weakly interacting molecules, A and B, referred to as monomers later on.Density matrices describing the electronic density localized at the monomers will be denoted as D ζ , ζ = A, B. These can be used to define the monomer F ζ matrix, which differs from the isolated monomer Kohn-Sham matrix by the presence of nuclear attraction operator of the other interacting molecule: where T is the matrix of the kinetic energy operator, U ne A and U ne B are matrices of the nucleielectron attraction operator of monomer A and B, respectively, J is the matrix of the Hartree potential, and finally, V xc represents the matrix of the exchange-correlation potential.Using Eq. ( 1) we set up a system of equations for interacting monomers: 22,23    where C i is a i-th column of the MO coefficients matrix of monomer A, C A , and similarly for C j .Eq. ( 2) is equivalent to a single Kohn-Sham equation for the dimer, provided that the occupied orbitals are orthogonal: In the A monomer part of Eq. ( 2), matrix J B describes direct Coulomb interaction with electrons localized at monomer B. Its matrix element in OAO basis is computed as where (τ υ|ρσ) denotes a two-electron integral in Mulliken notation and S −1/2 µτ is an element of inverse square root of the AO overlap matrix.∆V xc A is the matrix of the non-additive exchange-correlation potential: with the density matrix of dimer AB being By C occ we denote the matrix of MO coefficients expressed in OAO basis.The orthogonality of orbitals that belong to different monomers implies idempotency of supermolecule density matrix: as should be for the solution of the Kohn-Sham equation.We emphasize that Eq. ( 2) is equivalent to a single Kohn-Sham equation for the AB supermolecule.However, the separation into monomer equations is a good starting point for introducing approximations in the context of intermolecular interactions.This work relies on dispersion-free approximation, 22 in which interaction of monomers is treated at the HF exchange-only level, as the exchangecorrelation non-additivity is approximated by the HF exchange matrix: where The dispersion-free energy of the AB system that corresponds to the approximation of Eq. ( 8) can be expressed as where monomers' density matrices, D A and D B , are obtained by converging Eq. ( 2).E HF and E KS are Hartree-Fock and Kohn-Sham energies.The KS energy is evaluated as where U nn A denotes the nuclear repulsion energy and xc is the exchange-correlation energy density.Substituting Hartree-Fock exchange for the KS exchange-correlation term in the above expression yields E HF A .The solution of Eq. ( 2) is performed in a freeze-and-thaw manner, 33 that is, at a given step a density of one monomer is optimized with the density of the other monomer frozen.
From the numerical point of view, the difference between our approach and the approach of Wesołowski and Weber 33 is that we enforce orthogonality of the monomers' occupied orbitals, and the potential of non-additive kinetic energy does not appear in our model.In order to efficiently converge Eq. ( 2) in a self-consistent field (SCF) manner, and to minimize the number of intermediate matrices, the iterative solution of Eq. ( 2) is carried out with full non-additivity of the xc potential, i.e. without the approximation of Eq. ( 8).Thus, in Eq. ( 2) the KS matrix of the whole AB system is present in each of two monomers' equations and Eq. ( 2) becomes a mere rearrangement of the KS equation for the dimer.The rationale behind such a procedure is as follows.If we start the freeze-and-thaw iterations from isolated monomers' densities slightly changed due to initial symmetric orthogonalization of the occupied orbitals, and then continuously deform each monomer's density in the field of its partner (keeping the orthogonality constraint), we converge to the orbitals which are well localized on the respective monomer.They can be interpreted as orbitals distorted by the intermolecular interaction.The above scheme is essentially identical to the Hartree-Fock Pauli blockade approach of Gutowski and Piela, 32 who used equation analogous to Eq. ( 2) (with the dimer Fock matrix present in each of the monomer equations) to meaningfully decompose the interaction energy.In our scheme, upon convergence, the localized orbitals are utilized in Eq. ( 10) to yield dispersion-free energy of the dimer.As the dispersion-free energy of Eq. ( 10) is not invariant with respect to the unitary transformation mixing A and B orbitals, the localization of orbitals on monomers is indispensable in our approach.In principle, post-SCF localization methods 34,35 could be used, however, the conceptual simplicity and computational scaling favors our approach of converging orbitals while keeping them localized on the monomers.Finally, the total interaction energy is obtained by combining E PB AB with the dispersion contribution: where E KS A and E KS B are KS energies of isolated monomers.

III. ALGORITHMIC DETAILS
In this Section we derive and discuss an algorithm that solves Eq. ( 2).The new scheme alleviates the convergence problems of previously reported penalty function approach. 22,23,32 our experience the original penalty function approach to PB method suffers from a divergent and oscillatory behavior of the SCF iterations in many cases, e.g., for most of the dipole-interaction and charge transfer complexes reported in this work.The new algorithm presented below has convergence properties similar to those of standard single point Kohn-Sham calculations for the AB dimer.
In order to optimize occupied orbitals of one monomer in the presence of its counterpart, we adopt an exponential ansatz of orbital rotation, Anti-symmetric matrix X A is defined as where i and a denote occupied and virtual orbitals, respectively.It follows from Eq. ( 14) that X A ia values determine the change of the occupied orbitals of monomer A due to interaction with monomer B. In the following, we will focus on how to obtain X ζ ia parameters variationally.As occupied orbitals localized at monomer B are excluded from summation in Eq. ( 15), the orthogonality between the occupied orbitals of A and B is maintained through the orbital rotation after it has been established in the zeroth iteration of Eq. ( 2).An iterative solution of Eq. ( 2) may proceed as follows: 1. Generate a set of orthogonal occupied orbitals by symmetric orthogonalization of converged occupied MO vectors of isolated monomers.Compute virtual orbitals as a complement to the occupied space.

Perform several microiterations to optimize occupied orbitals of monomer A interacting
with monomer B. Rotate virtual orbitals of the whole system along with the occupied orbitals: 3. Analogously to the above, optimize occupied orbitals of monomer B interacting with monomer A: 4. Go to step 2 if convergence is not achieved.Exit the loop if the density matrix error is acceptable.
Variational parameters contained in the X ζ matrix are determined iteratively through the Kohn-Sham total energy (E KS ) minimization.We cast the augmented Roothaan-Hall master equation 36,37 for the monomer A, in a form that is suitable for restricting the variational space: where ( 22) For clarity, we drop the A index at F, D, and X.The notation introduced in Ref. 37 is used in Eq. (20).Here, Fk and D k are Kohn-Sham and the density matrices, respectively, from the k-th iteration.Eq. ( 21) is written in the basis in which occupied-occupied and virtualvirtual blocks of the F matrix are diagonal.Such vectors can be computed by solving a pair of eigenvalue problems for occupied-occupied and virtual-virtual blocks of the KS matrix, and transforming eigenvectors to the OAO basis: During the microiterations optimizing the A monomer orbitals we solve Eq. ( 21) with the condition X ia = 0 for i ∈ B, that is, we substitute the matrix of Eq. ( 15).This is necessary to prevent the orbitals of monomer A from mixing with the occupied orbitals of monomer B.
To facilitate solving Eq. ( 21), we define the intermediates: With the new definitions, Eq. ( 21) can be reformulated as We can determine σ l (X) on the RHS of Eq. ( 32) if both sides of Eq. ( 32) are multiplied by −2 (D m ) ia and summed over i, a: Eq. ( 33) can be rearranged into the form of a (n − 1) × (n − 1) linear equation: In our implementation the maximum value of n equals 8. Matrix exponential in Eq. ( 14) can be conveniently computed by squaring and scaling method: Optimal (j, k) pairs for a given norm of matrix A were tabulated by Moler and Van Loan. 38

IV. RESULTS AND DISCUSSION
In order to assess the applicability of the dispersion-free approximation within the PB framework and to compare its results with the currently available benchmarks, we employed the database of noncovalent dimers of Zhao and Truhlar 30,31 .Their database gathers interacting molecules in subsets according to the dominant character of the interaction.For the cases in which the relative discrepancy between the PB+dispersion and CCSD(T) results is noticeable, we additionally performed the symmetry-adapted perturbation theory (SAPT) calculations.All SAPT calculations reported in this work are accurate through the second order in the intermolecular interaction operator [see Eq. ( 39)].They were performed with the Molpro package. 39Subroutines to carry out the PB calculations can be obtained directly from the authors.

A. SAPT Overview
The SAPT method allows one to represent the interaction energy as the sum of physically meaningful terms of the intermolecular perturbation theory: 40 elst + E exch + E ind + E (2) The first-order terms on the RHS of Eq. ( 39) are the electrostatic and exchange energies, and the second-order contributions include induction, exchange-induction, dispersion, and exchange-dispersion energies, respectively.The method basically comes in two variants.
In the Hartree-Fock SAPT [SAPT(HF)] the monomers are described by uncorrelated HF orbitals, i.e., at the zeroth-order with respect to the intramonomer correlation operator.This is denoted customarily by introducing the second superscript at the contributions to interaction energy, e.g.E disp .In the DFT-based SAPT [SAPT(DFT)], 17,41,42 the monomers are described by the KS orbitals which are presumed to include the intramonomer correlation.The second-order terms are obtained from the coupled KS (CKS) theory.The SAPT(DFT) interaction energy of Eq. ( 39) is often augmented by the residual HF term defined as elst + E (10)   exch + E ind + E exch-ind , with E HF int being the supermolecular HF interaction energy.The purpose of the δ HF term is to account for higher-order induction effects and other residual terms necessary to obtain a convergent result.Thus, the full SAPT(DFT) interaction energy is Both PB and the cited reference calculations were performed with with counterpoise correction.The interaction energies do not include contributions from deformation of monomers.Reference CCSD(T) interaction energies and SAPT(DFT) dispersion contributions, are taken from Ref. 21.They were computed using the aug-cc-pVTZ basis set with midbond functions.The dispersion energy in Eq. ( 13) is either taken from SAPT(DFT): or from the numerical fit to the dispersion and exchange-dispersion terms of Ref. 43, E KS int stands for the supermolecular KS interaction energy.The SAPT values in Table VI are computed in aug-cc-pVTZ basis set with the asymptoticallycorrected PBE0 functional.SAPT(DFT)δ interaction energies include δ HF which is also listed separately.

HB6/04
The results for six complexes form the HB6/04 database are presented in Table I.In the assessment of the dispersion-free approach E PB int in combination with two formulations of dispersion energy (see Sec. IV B) is compared with the supermolecular interaction energies E KS int and E

CCSD(T) int
. The E PB int values show practically no dependence on the basis set, aug-cc-pVDZ vs. aug-cc-pVTZ and only a weak dependence on the type of the hybrid functional, PBE0 vs. B3LYP.The difference between E KS int and E PB int arises from the different xc potential acting between the interacting monomers: the same as the intramonomer potential in the former and the exact exchange in the latter.The differences between the two quantities are considerable and roughly similar (but not equal) to the dispersion energy.For the first four complexes PB+d or PB+df are in excellent agreement with the CCSD(T) benchmark.For the more strongly-bound complexes with double hydrogen bonds as (HCOOH) 2 the agreement deteriorates to the 20% overbinding.

CT7/04
The group of charge-transfer complexes consists of donor-acceptor pairs ordered according to increasing strength of interaction.The dispersion-free energies, E PB int , are always above the total reference interaction energies.However, after the addition of dispersion contribution, severe overbinding is present except for the weak complexes involving F 2 .
For the complexes of ClF (particularly NH 3 •••ClF) the discrepancy between PB+d (PB+df) and CCSD(T) reaches 55% (B3LYP) to over 70% (PBE0).The question is whether the reason for this behavior is the failure of the PB method for strongly overlapping monomers, or the problem with xc potentials which fail to properly match the donor-acceptor properties of subsystems.This isse is discussed in more detail below in Sec.IV C 7.

DI6/04
The group of molecular complexes in this database is primarily bound by the dispersion energy as E PB int represents only a fraction of the total interaction energy.The PB+d(df) treatment faithfully represents the bonding trends of this group as compared to CCSD(T) with a slight overbinding of about 10-15% (with one outlier, CH 3 SH•••HCl, with 25% overbinding).

WI7/05
In this group of van der Waals complexes E PB int predicts all interactions to be repulsive in both functionals.This is correct behavior expected from a dispersion-free approximation.In combination with the dispersion energy, the PB+d(df) approach yields the interaction energies in excellent agreement with CCSD(T).It should be pointed out that the supermolecular E KS int results obtained in PBE0 and B3LYP show trends that are opposite to one another.This is in striking contrast to E PB int which remains virtually independent of the functional.We note that similar behavior of the dispersion-free energy with respect to the functional was observed by Kamiya, Tsuneda, and Hirao 20 .

PPS5/05
This database contains representative π-stacking complexes.The dispersion free approximation yields strongly repulsive interaction energies for all but one of the complexes.The repulsive interactions for the three conformers of benzene dimer show consistency between two functionals PBE0 and B3LYP and two different basis sets: two stacking configurations, sandwich and displaced have equal interactions and the T-shaped one is less repulsive.When combined with the dispersion energy, PB+d(df) follows closely the CCSD(T) benchmarks.
The slight overbinding of about 10 % still persists.Unlike for the HB6/04 complexes the E KS int interaction energies are in poor agreement with both PB+d(df) and the benchmarks.The ordering of the three conformers of the benzene dimer is wrong in both functionals.

Comparison with SAPT and DFT-based approaches
In Table VI  agreement is contingent on including the δ HF term in the interaction energy.In these three complexes the δ HF is very large in magnitude and in NH 3 •••ClF it is responsible for the entire binding effect in SAPT(DFT)δ.The large δ HF is in itself and indicator that the convergence of SAPT is problematic. 44,45In Ne 2 and T-shaped benzene dimer the PB+d approach is closer to CCSD(T).The performance of SAPT(HF) appears to be erratic.On the one hand, in (HCOOH) 2 it agrees reasonably well with CCSD(T); on the other hand,   In Tables VII the PB+d/PBE0 method is compared against other DFT approaches specifically tailored for intermolecular interactions.Clearly, PB+d approach achieves high accuracy in WI7/05 and PPS5/05 sets when compared to other approximate methods, notably without additional parametrization.Its performance with respect to dispersion-dominated complexes is comparable to that of dlDF+D and B97-D methods, which probably represent the best of what the contemporary well-parametrized functionals can achieve for these interactions.However, inspecting first three columns of data in Tables VII and VIII reveals that the PB+d approximation worsens the results of the underlying functional (PBE0) and performs poorly compared to other approaches for the hydrogen-bonded and dipole-bound complexes.For the charge-transfer group all the approaches, except for dlDF+D, appear to have problems.

Charge-transfer interactions
Let us examine more closely the case of charge-transfer interactions.with CCSD(T) as well as with the supermolecular E KS int term.The PB interaction energy curve has a well depth of about 10 kcal/mol, which is ca. 5 kcal/mol less than the E KS int curve and about 4 kcal/mol more than the E HF int curve.However, when combined with the dispersion energy from SAPT, PB+d descends rapidly.Clearly, these two terms are incompatible.
There are two reasons for this.First, because of a short intersystem distance, the PB approach already includes some dispersion effect (i.e.no longer serves as a dispersion-free approximation) and combining it with the SAPT-based dispersion leads to double counting.
Second, E KS int /PBE0 calculations (and to a lesser extent with B3LYP, see Table II) lead to a considerable overbinding, as evidenced by a 4 kcal/mol overestimation of the well depth of E KS int compared to CCSD(T) in Fig. 1.
We will demonstrate below, using PBE0 as an example, that this overbinding can be traced, at least partially, to the tendency of these functionals to make the donor molecule too good a donor and the acceptor molecule too good an acceptor of electrons.This problem has long been attributed 51 to the mismatch between the HOMO energy of the donor and the LUMO energy of the acceptor leading to a too narrow fundamental gap in some hybrid functionals.Alternatively, this problem can be related to the delocatization error in the approximate xc functionals. 9To detect the presence of the delocalization error we inspect the behavior of the total energy between the integer values of electron numbers.In the exact functional this dependence should be linear from the cation to the neutral and from the neutral to the anion, showing a discontinuity in the first derivative at the neutral. 52partures from the linear behavior with respect to fractional electron numbers indicate the presence of the delocalization error. 9,53,54For the two constituent molecules ClF and NH 3 we performed calculations of energy as a function of fractional electron numbers with the PBE0 functional.The calculations were performed with NWChem package. 55The longrange corrected PBE0 (LC-PBE0) is included for comparison.The results in Fig. 2 and 3 demonstrate that LC-PBE0 yields a linear relationship of the total energy with fractional electron numbers for both ClF and NH3, i.e., it behaves as the exact functional.For NH 3 a positive slope of neutral-to-anion curve indicates that NH 3 does not form a negative ion, unlike ClF which does.By contrast, PBE0 alone leads to convex dependence of energy with respect to fractional electron numbers for both molecules.This means that in the PBE0 functional ClF becomes more prone to accepting a fractional charge compared to the

V. CONCLUSIONS
In this work we presented comprehensive tests of recently-developed dispersion-free functional approach.
Except for dispersion-dominated complexes for which PB+d yields excellent results, systematic overbinding has been observed.This error is most severe for strongly bound dimers with short intersystem distances where the distinction between the intrasystem and intersystem correlation effects becomes problematic.In these complexes combining the PB We also acknowledge computational resources of the Interdisciplinary Centre for Mathematical and Computational Modelling of the University of Warsaw.

2 TABLE
III.Comparison of interaction energies for the complexes from DI6/04 database.All values are in kcal/mol.System PBE0 / aug-cc-pVTZ PBE0 / aug-cc-pVDZ B3LYP / aug-cc-pVDZ the results of the PB+d approach are compared with two versions of the second-order SAPT: one based on the HF description of monomers [SAPT(HF)] and one based on the DFT monomer description [SAPT(DFT)δ].As described in Sec.IV A, SAPT(DFT)δ incorporates the HF-derived residual term, δ HF , into the interaction energy.For the comparison we chose three problematic cases discussed above (HCOOH) 2 , NH 3 •••ClF, CH 3 SH•••HCl, as well as Ne 2 and T-shaped (C 6 H 6 ) 2 which are notoriously difficult for theory.At the first glance SAPT(DFT)δ appears to provide much better agreement with the CCSD(T) benchmarks than PB+d for (HCOOH) 2 , NH 3 •••ClF, CH 3 SH•••HCl.However, this Comparison of interaction energies for the complexes from PPS5/05 database.All values are in kcal/mol.System PBE0 / aug-cc-pVTZ PBE0 / aug-cc-pVDZ B3LYP / aug-cc-pVDZ

Fig. 1
shows the dependence of E PB int and E PB+d int in NH 3 •••ClF on the intersystem distance and the comparison

FIG. 1 .
FIG. 1.The interaction energy curves for the NH 3 •••ClF complex.R denotes the distance between N and Cl.The PB and KS interaction energies are computed with the PBE0 functional and aug-cc-pVTZ basis set.The HF and CCSD(T) terms are computed with the aug-cc-pVQZ basis set.

FIG. 2 .
FIG. 2. Energy of ClF as a function of fractional electron number ∆n relative to the neutral system (∆n = 0).The LC-PBE0 functional employs the range separation parameter 0.3, a portion of 0.25 HF exchange in the short range, and 1.0 HF exchange in the long range.All calculations are performed in aug-cc-pVDZ basis set.

FIG. 3 .
FIG. 3. Energy of NH 3 as a function of fractional electron number ∆n relative to the neutral system (∆n = 0).The LC-PBE0 functional employs the range separation parameter 0.3, a portion of 0.25 HF exchange in the short range, and 1.0 HF exchange in the long range.All calculations are performed in aug-cc-pVDZ basis set.

TABLE I .
Comparison of interaction energies for the complexes from HB6/04 database.All values are in kcal/mol.

TABLE IV .
Comparison of interaction energies for the complexes from WI7/05 database.All values are in kcal/mol.

TABLE VI .
Comparison of the PB+d approach with SAPT(HF) and SAPT(DFT).The PB+d values are obtained with intramonomer functional PBE0 while SAPT(DFT) terms are obtained with asymptotically corrected variant of PBE0, PBE0AC.δ HF denotes the residual HF terms (see the text).All values are in kcal/mol.

TABLE VII .
Mean unsigned percentage error of DFT-based approaches.All results were obtained with aug-cc-pVTZ basis set.Results other than PB+d or supermolecular PBE0 were taken from