Comparison of approximate intermolecular potentials for ab initio fragment calculations on medium sized N‐heterocycles

Abstract The ground state intermolecular potential of bimolecular complexes of N‐heterocycles is analyzed for the impact of individual terms in the interaction energy as provided by various, conceptually different theories. Novel combinations with several formulations of the electrostatic, Pauli repulsion, and dispersion contributions are tested at both short‐ and long‐distance sides of the potential energy surface, for various alignments of the pyrrole dimer as well as the cytosine–uracil complex. The integration of a DFT/CCSD density embedding scheme, with dispersion terms from the effective fragment potential (EFP) method is found to provide good agreement with a reference CCSD(T) potential overall; simultaneously, a quantum mechanics/molecular mechanics approach using CHELPG atomic point charges for the electrostatic interaction, augmented by EFP dispersion and Pauli repulsion, comes also close to the reference result. Both schemes have the advantage of not relying on predefined force fields; rather, the interaction parameters can be determined for the system under study, thus being excellent candidates for ab initio modeling.


| INTRODUCTION
Ab initio quantum chemistry is now an indispensable part of chemical research. Although a large arsenal of methods is available, the cost of calculations grows rapidly with system size; and, the more advanced the method, the faster the increase in cost. A possible way to overcome this problem is the fragmentation of the system into smaller components and applying an appropriate formulation to account for their interaction. The number of methods of this kind is tremendous, yet there are still many challenges. One of the biggest issues is how to cut chemical bonds during fragmentation.
A specific situation is presented by molecular complexes, whose fragments are molecules connected by non-covalent interactions. In this case the Hamiltonian of a two-component system is trivially par- preliminary calculations showed that the most often used, purely electrostatic approximations fail to provide the correct potential energy surface, because further contributions are significant. 10 This is particularly true when the short-distance part of the potential is of interest, for example for stacked interactions in DNA.
To find the best-suited approximation for these terms, the present paper investigates the ground states of non-covalently interacting bimolecular complexes and compares various methods for different contributions to the interaction energy. Keeping our focus on nitrogen-containing heterocycles, our test systems include pyrrole dimers (denoted hereafter as Pyr-Pyr) in different orientations, as well as stacked cytosine-uracil base pairs (Cyt-Ura). These complexes have low-lying Frenkel-type pairs of excited states for which the proper handling of non-covalent interactions is needed to obtain accurate splitting.
The paper is organized as follows. Section 2 gives a short summary of the available methodologies, with emphasis on the differences in theoretical formulations and the possibility of incorporating them into ab initio calculations. Based on this analysis, in Section 3 we select a set of methods to be used in the tests. Section 4 describes the computational details, and Section 5 presents the results.

| METHODOLOGIES OF MODELING INTERMOLECULAR INTERACTIONS
Already at the dawn of quantum chemistry, it became clear that quantum effects give rise to important short-range forces. 11,12 By now, the treatment of intermolecular interactions has a broad theoretical background. Systematic overviews of the topic have been presented in comprehensive texts. [13][14][15][16][17] The contributions to the interaction energy are usually referred to as 1,15 : 1. Electrostatic (eventually including polarization effects);
One way to distinguish these terms is offered by the SAPT theory 1,2 that, in different orders of perturbation theory, provides approximate definitions for the individual terms of interaction. Of course, one should always keep in mind that distinguishing these contributions is physically not strict since an exact (or accurate) calculation on the complete system would give just one single, interaction energy. Nevertheless, from the perspective of chemical and physical intuition, and also for developing reasonable approximations, this distinction is well justified.
In SAPT, the electrostatic energy enters in first order in the interaction operator, with sizable contributions from the correlation terms. 1,2 Thus, a careful modeling of this contribution is essential.
The electrostatic interaction is traditionally incorporated into the model via charges placed on the interacting molecules. 18  The ways of defining the partial charges have been reviewed recently by Janeček et al. 22 Here we just mention that the population analysis concept used in the past is nowadays replaced by fitting directly to ab initio results on interaction energies, dipole moments or electrostatic potentials. For the present study CHELPG 23 is of special importance since it represents a readily available algorithm to calculate atomic charges for the particular monomer of interest. A recent advance in this field is the W-RESP potential of Janeček et al. 22 A more flexible parametrization of the electrostatic interaction can be achieved using atomic multipoles, first suggested by Stone. 14,24 Pioneering work in applying this technique was done within the Effective Fragment Potential (EFP) framework. 3 39 or as a standalone extension 40 ; it has also been implemented in the Gaussian program package 41 for QM/MM type electrostatic modeling.
Two additional energy components decay fast with the distance of the fragments but proved to be important in SAPT theory. 1 One is the dispersion energy, a result of correlation of the respective electrons in the two fragments, which appears first in the second order perturbation terms, thus correlation is clearly essential in it. Second is the Pauli (exchange) repulsion, which arises when antisymmetrizing the wave function of the complete system. 42 In molecular force fields, dispersion and Pauli repulsion are often considered together as one van der Waals potential, traced back to Lennard-Jones. 43 A comprehensive historical summary of the subject was published recently. 44 Lennard-Jones (LJ)   and EFP. [47][48][49] For dispersion, EFP2 calculates the so-called C 6 coefficient of the London expansion as the instantaneous dipole-induced-dipole interaction, obtained from the frequency dependent polarizabilities; it also includes a C 8 term estimated as one-third of C 6 . 50 The different versions of D corrections (D, D2, D3) on density functional theory (DFT), to include dispersion, were developed by Grimme et al. [51][52][53] for molecules, but these can also be used to obtain intermolecular dispersion energy. Similarly to the respective EFP2 term, it is based on a London expansion up to C 8 , with C 6 calculated from hydrides of the given atoms. It also includes a theoretically justifiable damping factor and three-body terms.
Beyond these a posteriori calculations, a more sophisticated way to treat intermolecular interactions is to explicitly include them into the Hamiltonian of the fragments. These models belong to the general family of embedding methods. 54 In these schemes the subsystem(s) is To overcome these problems, an embedding version of EFP2, termed QM/EFP2 or ab initio /EFP2 has been introduced recently. 25 The embedding scheme also facilitates the treatment of all subsystems at QM level, with the advantage that the interaction with the electron density of the environment can directly be calculated, without discretization of the charge distribution or parametrization of other contributions. Usually, the density of the surroundings is calculated at a lower level of theory. The Kohn-Sham (KS) DFT offers an excellent framework for this QM based embedding as it scales favorably with the system size while providing a relatively accurate description of the environment's density. In a typical DFT embedding scheme the energy of the system is written as where This technique employs a subtractive embedding scheme where a DFT calculation is performed on the entire system, which is then divided into an active and an environment part. The energy of an embedding calculation of this type can be given as where E WF e A is calculated in the presence of the environment's density, including thereby the electrostatic interaction of the subsystems.
Note that the exchange interaction between the two subsystems is retained from the original super system calculation and included in the final energy. However, since the underlying DFT calculation on the super system does not describe dispersion correctly, this contribution needs to be added afterwards.
In this work we test the two main categories of methods for describing intermolecular interactions. First, a posteriori corrections sourced from representative force fields including EFP2, the latter considered as a system-specific FF. Second, embedding techniques, which, to our knowledge, were not yet studied systematically for intermolecular interactions beyond QM/MM. As comparison, SAPT interaction terms are also evaluated for the two stacked complexes.
Note, however, that SAPT can only be used as a ground state benchmark since it is not yet set to be generally used for excited states, although one should acknowledge a recent promising effort by Hapka et al. 70

| INTERACTION POTENTIALS
As outlined above, we want to test individual contributions to intermolecular non-covalent energy, as treated in various methods. The aim is to find a combination of contributions, which is accurate enough for ground states, can be combined with various QM methods, and eventually can also be extended to excited states. Our approach is pragmatic in the sense that instead of seeking general conclusions, we intend to identify the best potentials/procedures to describe complexes of N-heterocycles at the vicinity of their equilibrium structures. Table 1 summarizes the methods we are going to explore.

| Ab initio calculations
The reference interaction energies were calculated at the CCSD(T)/ aug-cc-pVDZ level 71

| Molecular systems
The test systems under study are illustrated in Figure 1. In the discussion below, the distance of monomers refers to that between the centers of masses.
For the pyrrole dimers six conformers were studied. Three of them are co-planar (side to side orientations), denoted Pyr-Pyr (ILn)). In Pyr-Pyr

| Dispersion and Pauli repulsion
In Figure 2  contribution, which results in a slight overestimation of the magnitude of this term (see Figure S1 of the Supplementary material.) case this means that the choice between D3 and EFP2 is less relevant for the co-planar structures, while for the T-shaped and sandwich configurations the difference can go up to 0.002 a.u. or even larger (0.01 a.u.), as in the case of Cyt-Ura (S).
The dispersion plus Pauli repulsion potential of EFP2 (EFP2 Rep +-EFP2 Disp) corresponds, as also seen in Table 1, to the effects covered by the GAFF Lennard-Jones potential. The behavior of the GAFF LJ potential (blue curves on Figure 2) with respect to the corresponding EFP2 curves (red curve) is not systematic. While EFP2 is generally more repulsive, for Pyr-Pyr (IL2) the opposite is true. Although the general shape of the curves appears to be similar, substantial energy differences can occur even around the equilibrium distances of some structures. For to that of other structures (À0.002 a.u.). This is because for this system the dispersion contribution in GAFF LJ is small, similar to D3, at least near equilibrium (see the corresponding curves in Figure 2). This finding warrants the favorization of EFP2 potentials over Lennard-Jones type ones for the dispersion contribution in similar systems. One has to keep in mind, however, that other force fields can possibly show a better performance than GAFF LJ parameters. which support many of the conclusions above. Not only appears EFP2 dispersion as preferred over D3 but also when combined with the Pauli repulsion, the EFP2 follows the SAPT curve better than GAFF LJ does. Nevertheless, the repulsive nature kicks in faster with SAPT than with EFP2.

| Electrostatic interaction
In Figure

| Total interaction energy
In Figure Table 2 where the different contributions to the F I G U R E 3 Distance dependence of various potentials describing electrostatic and Pauli-repulsion contribution of the interaction energy. d e marks the reference equilibrium distance total interaction energy are listed at the respective equilibrium distance of the complex under study.
Pyr-Pyr (IL1) structure: Despite the fact that the dipole moments of the two monomers are aligned in this system, there is a very weak interaction between the monomers. Contrary to e.g. the sandwich systems discussed below, here all three contributions to the interaction energy are small (see Table 2). Both the EFP2 Tot and the (CHE- Pyr-Pyr (IL2) structure: The total interaction energy is even smaller in this structure than in Pyr-Pyr (IL1). This is partly because the electrostatic part is repulsive, thus compensating dispersion. Here the best agreement with the reference is provided by (Embed + EFP2 Disp), although the minimum is at a somewhat larger distance. The note that in this case the dispersion parts of the interactions were also found to agree above.
Pyr-Pyr (T1) structure: The total interaction energy in this system is the second largest among the investigated structures, with a strong attractive electrostatic component (see Table 2 Pyrrole-pyrrole sandwich(Pyr-Pyr (S)): Describing the interaction in this system is a special challenge for theoretical methods. 86 Near the equilibrium, dispersion is the strongest contribution to the interaction energy (À5 mE h ), but it is largely canceled by the Pauli repulsion and the (also repulsive) electrostatics, resulting in a very small total interaction energy. It is thus an example where good result can only be achieved with a method where a delicate balance between the different contributions is present. Indeed, as Figure 4 shows, there is a substantial discrepancy between the two embed- where GAFF LJ fails badly. Typically, the bonding provided by this combination is a bit weaker than in the reference, representing a quality comparable to EFP2 Tot and somewhat less accurate than embedding methods.
The dispersion corrected embedding curves (Embed + EFP2 Disp) and (Embed + D3 Disp) are in most cases closer to the reference than other methods investigated here. From these two, the former with the EFP2 dispersion clearly gives the best agreement, outperforming the latter, which uses D3 dispersion.

| CONCLUSIONS
Various approximate methods for the description of non-covalent intermolecular interactions were investigated on certain conformers of pyrrole-pyrrole and cytosine-uracil complexes as test systems. Specifically, the significance of various terms in the models and combinations thereof were analyzed. Despite the approximate, in some sense arbitrary definition of some contributions, we confirmed their general importance. In fact, it was found that all contributions of a particular model need to be included even for qualitatively correct potential energy curves.
With future studies in mind, we paid special attention to models that can be combined with any electronic structure methods. In this respect, QM/MM type formulations with point charges are ideal since they can easily be implemented in the one-electron Hamiltonian, without affecting the execution of any level of calculation. However, despite its widespread, successful applications to, e.g., solvent effects, QM/MM is not suited per se for weakly bonded dimers. In our test systems, other terms, such as dispersion or Pauli repulsion give nonnegligible, often dominating contributions. Thus, the electrostatic part needs to be augmented with the latter. Most force fields contain a Lennard-Jones type potential for this purpose; our calculations indicate, however, that more sophisticated methods, e.g. the Effective Fragment Potential technique, may be needed.
At the same time, one has to keep in mind that different models address the division of interaction into individual terms by different strategies. Therefore, it is not surprising that combinations of different approaches tend to give inconsistent or controversial results.
Still, we found certain setups of promising performance, approaching the reference CCSD(T) and, where available, also the SAPT results. In particular, a DFT/CCSD embedding model covering electrostatic and Pauli repulsion effects, if augmented with dispersion via the respective EFP2 term, proved to be a reasonable working strategy. It has shown remarkable accuracy specifically for the stacked interactions. A QM/MM style electrostatics using CHE-LPG charges, with the point charges determined at the partner molecule's nuclei is also a viable alternative to the widely adopted force field charges, also liberating the description from the dependence on precalculated parameters and atom types. Combination with either the EFP2 dispersion plus Pauli repulsion or with a Lennard-Jones potential of the GAFF force field gave similar results, the latter being slightly better in some examples; but the latter finding is probably specific to the choice of this particular FF and may not be achievable for systems with different atom types. Also, the LJ contribution is disturbingly random: depending on the specific arrangement, it may severely overestimate, in other cases underestimate the interaction. Thus, its application would be difficult to support and we feel that the idea of a general, system-adopted description provided by EFP2 is comforting enough to be favored over the LJ approach of any force field.
The D3 dispersion correction, a strategy also founded on transferable atomic parameters, was clearly outperformed by EFP2 in our results (This finding is, of course, not related to the original applicability of D-type corrections in DFT).
In summary, among the many alternative formulations arising from the synthesis of different theories we have found combinations, which not only offer good performance, but also retain the generality and theoretical desirability of ab initio modeling. For future studies, especially our planned work on excited states, the sophisticated description offered by DFT/CCSD embedding and the easily implementable QM/MM approach with CHELPG charges, both augmented with missing terms from EFP, are the recommended tools.