Chargetransfer from Regularized SymmetryAdapted Perturbation Theory
Abstract
The chargetransfer (CT) together with the polarization energy appears at second and higher orders in symmetryadapted perturbation theory (SAPT). At present there is no theoretically compelling way of isolating the chargetransfer energy that is simultaneously basisset independent and applicable for arbitrary intermolecular separation. We argue that the chargetransfer can be interpreted as a tunneling phenomenon, and can therefore be defined via regularized SAPT. This leads to a physically convincing, basisindependent definition of the chargetransfer energy that captures subtleties of the process, such as the asymmetry in the forward and backward charge transfer, as well as secondary transfer effects. With this definition of the chargetransfer the damping needed for polarization models can be determined with a level of confidence hitherto not possible.
pacs:
34.20.Gj 31.10.+z 31.15.pI Introduction
Symmetryadapted perturbation theory (SAPT) remains one of the most accurate and versatile methods for calculating intermolecular interaction energies. Within SAPT the interaction energy is decomposed into physical components such as the electrostatic, exchangerepulsion, induction and dispersion energies. This decomposition is useful both as an interpretative tool and as a means for developing models for the interaction energy. The latter is made possible as, with the exception of the shortranged exponentially decaying exchangerepulsion energy, each of these components is associated with a welldefined multipole expansion, the coefficients of which can be calculated from monomer properties alone.
There is, however, an important exception to the list of physical components of the interaction energy as described by SAPT: the chargetransfer energy is not defined as a separate component, but is included, together with the polarization energy, in the induction energy at second and higher orders in the interaction operator. While these components are not usually separated in a SAPT calculation of the interaction energy, there are approximate methods, such as that introduced by Stone Stone (1993) and used by Stone & Misquitta Stone and Misquitta (2009) in the context of a variant of SAPT based on a densityfunctional description of the monomers (termed SAPT(DFT)) Stone and Misquitta (2009), which enable us to make a reasonable partitioning of the induction energy—at least at secondorder—into its chargetransfer and polarization constituents. Though these methods, and their supermolecular counterparts Khaliullin et al. (2008), provide sensible chargetransfer energies, particularly for complexes at their equilibrium configurations, they exhibit serious drawbacks, all related to their basisset dependence for small intermolecular separations. Furthermore, as we shall demonstrate in this paper, the Stone & Misquitta procedure leads to chargetransfer energies that do not seem to reflect the physical nature of the process.
There are several reasons why we require a rigorous definition of the chargetransfer energy.

Polarization models: The polarization energy can be unambiguously defined through molecular polarizabilities and multipole moments, but only for large molecular separations where chargedensity overlap is negligible. At shorter separations the expansion must be damped, and to determine the damping we require a knowledge of the polarization energy, and hence the chargetransfer energy.

Validation of density functionals: One of the remaining issues with the exchange functionals used in density functional theory (DFT) is their inability to describe the charge transfer energy accurately Cohen et al. (2008). A considerable effort is being made to attempt to remedy this problem, but for this to work accurate benchmark chargetransfer energies are needed.

Energy decomposition analysis (EDA) techniques are a powerful means of developing an understanding of the nature of molecular interactions. However, with the exception of perturbative techniques like SAPT, these methods contain an element of ambiguity: energy components, particularly the chargetransfer and polarization energies, obtained from different EDAs can be significantly different. To benchmark EDA methods we need a physical and basisindependent definition of the interaction energy components. SAPT already provides such a decomposition for the electrostatic, exchangerepulsion and dispersion energies, but not the chargetransfer and polarization. This disadvantage needs to be overcome.
The second and third points are linked: EDA techniques are used to determine the physical content of DFT interaction energies, but to do so reliably, the EDA methods themselves need to be calibrated against more rigorous decomposition methods, such as SAPT. In this paper we will focus on the first and second: the role the chargetransfer plays in determining accurate polarization models, and the provision of benchmark chargetransfer energies.
Ii What is the chargetransfer energy?
Intermolecular chargetransfer is a groundstate phenomenon. When two molecules are in proximity, their electronic charge density can be delocalised, that is, shared between them, leading to a stabilisation by an amount that is the chargetransfer energy. However, from the point of view of the interacting monomers, the delocalisation can be viewed in terms of excitations partly localised on the partner monomer: these are the socalled intermolecular chargetransfer excitations. This is why the chargetransfer energy first appears at second order in intermolecular perturbation theory. Common to both views is the idea of charge delocalisation leading to stability. One may even suggest that the term ‘chargetransfer’ is misleading, and rather, ‘chargedelocalisation’ may be a more appropriate term for this phenomenon. We will return to this issue of nomenclature later in the paper.
This simple physical picture of the charge transfer process lies at the heart of current attempts to define the chargetransfer energy: if we can calculate interaction energies both by allowing and suppressing the delocalisation process, the chargetransfer energy can be defined as an energy difference.
ii.1 Basisspace definitions
Currently, definitions of the chargetransfer energy rely on basisspace localisation methods. There are many flavours of this kind of approach, all of which are based on the idea: that with a suitably localised basis set we suppress the chargetransfertype excitations, and consequently isolate the chargetransfer from the polarization energy. In the Stone & Misquitta Stone and Misquitta (2009) technique the chargetransfer energy is defined as the difference in the secondorder induction energy, , calculated in a dimer centered (DC) and monomer centered (MC) basis set. That is,
(1) 
Here, unless otherwise specified, the induction energy will be assumed to be the sum of the polarization and exchange components as Misquitta and Stone (2008a): . In the MC basis the molecule is described using basis functions located on its own nuclei only, therefore suppressing any excitation that could give rise to an increased electronic density on the partner monomer. That is, chargetransfer would be suppressed. By contrast, in the DC type of basis, the monomers are additionally described using basis functions located on their partners, therefore allowing chargetransfer type excitations. Hence we are led to the above definition of the CT energy which we will henceforth term CTSM09.
Note that in practical calculations, instead of the DC type of basis we often use an equivalent, but smaller, MC+ type of basis which includes midbond and the socalled farbond functions Williams et al. (1995). The farbond functions are a subset of the basis of the partner monomer; usually the s and pfunctions only. Extensive tests have shown that these two types of basis set result in nearly identical induction energies Williams et al. (1995). Consequently we will use the DC and MC+ types interchangeably in eq. (1) and will consider induction energies calculated in either basis type to include all CT effects.
There are, however, a few objections that one could raise with the CTSM09 definition:

Equation (1) leads to a basisdependent definition of the CT energy. As the monomer basis gets larger and more complete, the distinction between the MC and DC basis types grows smaller. This would lead to smaller apparent CT energies. This issue has been acknowledged by Stone & Misquitta who have argued that this effect is indeed present, but is small enough that for practical basis sets this is not and issue.

A more subtle, but related shortcoming of eq. (1) is that it inevitably leads to a separationdependent definition of the chargetransfer. As the intermolecular separation reduces, the MCtype basis set is better able to describe CTtype excitations, that is, the definition of the CT gets progressively worse.

The CTSM09 definition is really only the secondorder chargetransfer energy. There are contributions to the CT from third and higher orders in perturbation theory that, as we shall see, are comparable to the secondorder CT.

This definition relies on the use of localised atomic basis sets and offers no clear route to extension to planewave basis sets. Consequently, eq. (1) is largely limited to programs that use Gaussiantype orbitals. This is certainly not a significant issue, but it would be advantageous to develop a method that was independent of the type of basis set.
The first two points are illustrated in fig. 1 for the water dimer in its minimumenergy, Hbonded conformation. The variation of the CTSM09 energy with basis set is indeed small at the equilibrium separation, particularly for the two larger basis sets. However, as the intermolecular separation reduces, the differences between all basis sets grows, with the augccpVQZ basis resulting in the least amount of (secondorder) chargetransfer energy at short separations.
These short separations are important for several reasons. The separation of water molecules can decrease significantly in clusters due to a combination of manybody polarization and chargetransfer effects. For example, the oxygen–oxygen separation in the water dimer at its global energy minimum configuration is Å, but it is less than Å in the cage hexamer. This range is indicated in fig. 1. Presumably the intermolecular separation could be even smaller for larger clusters with stronger manybody effects, or for water under pressure.
The differences in CTSM09 calculated using the two larger basis sets may seem small even at the short end of the range, but the variation in CTSM09 with basis at even shorter separations is large enough that the polarization damping model—which is fitted largely to these energies—becomes basisset dependent with damping parameters varying from , to , through a.u. for the aDZ, aTZ and aQZ basis sets. Not only is this variation large, but the damping parameters tend to depend strongly on the range of data used—a possible indication that the energies we are fitting to are not entirely polarization, and might be contaminated with chargetransfer. This variation in the damping parameter has no effect on the twobody interaction energy, but causes the manybody polarization energy for clusters of water molecules to vary considerably. For example, with the above damping parameters, the manybody polarization energy of the cyclic water pentamer varies from , to , to . This variation is large enough to make these polarization models unreliable and almost unusable without explicitly fitting them to reproduce the manybody energies of such clusters Sebetci and Beran (2010). While this can be done, this approach is tedious and unsatisfactory as there is always the possibility that a damping models that works for one set of clusters may not work for another.
The Stone and Misquitta technique is but one of many methods that attempt to define the chargetransfer energy. Before continuing we mention two of these: Schenter and Glendening Schenter and Glendening (1996) have proposed a decomposition based on natural bond orbitals, which has the merit that the chargetransfer energies appear to converge with basis set, but these energies are unreasonably large. More recently Khaliullin et al. Khaliullin et al. (2008) developed a method that operates on a principle similar to that of Stone and Misquitta, only this time absolutely localised molecular orbitals (ALMOs) to suppress chargetransfertype excitations. In a recent analysis by Azar et al. Azar et al. (2013) the ALMO method has been demonstrated to exhibit shortcomings similar to those of the Stone and Misquitta definition: there is no complete basis set limit to the chargetransfer energy defined in this procedure.
Iii ChargeTransfer via regularized SAPT
The Stone–Misquitta and related approaches are all basisspace methods; that is, these methods attempt to isolate the chargetransfer energy by manipulating the atomic or molecular basis sets to try to suppress CTtype excitations. To see how we may find an alternative approach consider the RayleighSchrödinger perturbation expression for the secondorder induction energy of molecule A (exchange effects are included in a separate term):
(2) 
where the and are the exact eigenstates and eigenvalues of monomer X=A,B and is the total electrostatic potential of monomer B, that is, it arises from both electrons and nuclei and can be written as
(3) 
where and are the nuclear charges and position vectors of monomer B and is the unperturbed electronic chargedensity of B. In the last step in eq. (2) we have defined the firstorder induction wavefunction correction :
(4) 
Similar expressions exist for the secondorder induction energy of monomer B. The subscript ‘pol’ in eq. (2) indicates that this is the energy obtained from the socalled polarization expansion Jeziorski et al. (1994), that is, without the inclusion of exchange effects. This notation is unfortunate as it can lead to the erroneous conclusion that (A) is the polarization energy.
The third form of (A) in eq. (2) allows us to make the useful interpretation of the secondorder induction as the secondorder energy response to the static potential of the partner monomer (or, more generally, the potential arising from the entire environment of monomer A). This potential, eq. (3), consists of two parts: an attractive, but singular, term arising from the nuclei of B, and a finite, repulsive term arising from the electrons of B. At longrange, monomer A sees the nett effect of these two terms, but at shortrange, when chargedensities overlap, monomer A additionally responds to the singularity arising from the nuclear potential of B. The electronic density of B screens the singularity of the nuclear potential, but the extent of the screening depends on whether B is electronrich or electrondeficient. For an electronrich atom B (such as the electronegative oxygen atom in water) the nuclear potential is heavily screened, leading to little tunneling of the electronic charge density of A into the nuclear potential of B. However, for an electrondeficient atom (such as the electropositive hydrogen atom in a water molecule) the screening is weak, leading to a significant degree of tunneling.
This interpretation of the chargetransfer offers us to a natural way of suppressing the chargetransfer process altogether: if the potential well developed by the nuclear potential can be suitably eliminated, there would be no possibility of tunneling, and consequently no charge transfer. This can be achieved by regularizing the nuclear potential that appears in eq. (3) to suppress the nuclear well.
iii.1 Regularization
Regularization is a technique originally developed to accelerate the convergence of the class of SAPT theories. Here we will be concerned mainly with the version of SAPT based on symmetrized RayleighSchrödinger (SRS) perturbation theory. The subject of regularization and the role it plays in convergence of SAPT is comprehensively discussed in Refs. Patkowski et al., 2001a and Patkowski et al., 2004. Although a full discussion of convergence issues would be out of place in this article, we will discuss some of the important issues here as these will serve to place the concept of regularization in context.
Possibly the earliest major attempt to understand the convergence properties of intermolecular perturbation theories was made by Claverie Claverie (1971) who argued that the polarization expansion Jeziorski et al. (1994) (essentially manybody RayleighSchrödinger perturbation theory) on which SRS theory is based either diverged, or if it did converge, for systems with more than two electrons, it would converge to a Pauliforbidden state which would be more strongly bound than the physical ground state. These ideas were expanded by Kutzelnigg Kutzelnigg (1980) who showed that the polarization expansion indeed diverged. Simultaneously, Morgan and Simon Morgan III and Simon (1980) proved that if any of the interacting atoms had atomic number greater than two, the physical ground state of the system could be buried in a continuum of Pauliforbidden states.
The cause of the divergence of the polarization expansion, and consequently SRS theory Adams (1999), lies in the presence of the unphysical states. Adams explored this issue in a series of papers (see for example Ref. Adams, 1999) and showed that the concept of regularization introduced by Herring Herring (1962) could be used to destabilize these unphysical tunneling states and hence to ensure the convergence of the theory. Simultaneously, Patkowski and collaborators Patkowski et al. (2001b); Patkowski et al. (2004) applied similar ideas to develop regularized perturbation theories that converged to the physical ground state of the dimer. The idea here is to write the singular electron–nuclear potential as a shortranged, singular part and a longranged part that is wellbehaved. In the notation used by Patkowski et al. this is expressed as
(5) 
where is the singular, shortranged part and the longranged, wellbehaved part of the nuclear potential. Various schemes can be used to achieve this splitting, here we will use the Gaussianbased scheme Patkowski et al. (2001b):
(6) 
Patkowski et al. used SRS theory for the nonsingular part () and a strong symmetryforcing theory (using a symmetrization technique completely suppressing the Pauliforbidden terms) for the singular part Patkowski et al. (2001b); Adams (2002), leading to a unified theory proven to be convergent for the van der Waals and as well as chemical bonding separationsPatkowski et al. (2004).
iii.2 Chargetransfer via regularization
The issues discussed above are fundamental to foundation and development of perturbation theories, but are not germane at low orders in perturbation theory where even weak symmetryforcing theories like SRS theory are known to yield sensible results Jeziorski et al. (1994); Patkowski et al. (2004); Adams (2002). There is a tremendous body of work demonstrating the accuracy of loworder SAPT (truncated at second or third order with higherorder corrections estimated in a nonperturbative manner) for numerous systems (for example see Refs. Jeziorski and Szalewicz, 1998, 2002; Misquitta et al., 2005; Podeszwa et al., 2006a; Bukowski et al., 2006). For such a truncated theory regularization has a very different role to play: following on from our discussion of the chargetransfer, we argue that, if used at second and higher orders in perturbation theory, the regularization procedure can be used to destabilize the chargetransfer states and consequently allow us to define a chargetransfer free interaction energy.
At second order in perturbation theory we calculate the regularized induction energy, , that, with an appropriate amount of regularization that is yet to be defined, is free of chargetransfer contributions (details below). Then, in a manner analogous with eq. (1), we may define the secondorder chargetransfer energy as
(7) 
Here is the regularized secondorder induction energy that may be identified with the true polarization energy. There is no basis restriction on the above definition, except that the basis set used needs to be large enough to converge the total induction energy, which is the usual requirement for any energy calculation.
iii.3 Implementation
We have implemented a version of regularized SRS theory (RSRS) derived by Patkowski et al. in Ref. Patkowski et al. (2012). In RSRS theory the dimer wavefunction corrections are obtained in response to the interaction operator with regularized nuclear potentials, but the interaction energy corrections are calculated using the original, unregularized interaction operator. For the secondorder induction energy this means that instead of using , we calculate firstorder induction wavefunction correction in response to the regularized electrostatic potential :
(8) 
to give
(9) 
The regularised secondorder polarization component of the induction energy is then defined as:
(10) 
To this, as always, we need to add the similarly regularized exchangeinduction energy Jeziorski et al. (1994); Patkowski et al. (2004); Adams (2002). Notice that in this step we have used the full electrostatic potential of monomer B.
The SAPT(DFT) expression for the polarization part of the induction energy of monomer A in response to the field of B is
(11) 
where and label occupied and virtual states, are matrix elements of the unperturbed potential of monomer B given in eq. (3), and the amplitudes are obtained, in the case of SAPT(DFT), by solving the coupled Kohn–Sham equations
(12) 
where is the electric Hessian Casida (1995); Colwell et al. (1995) of Kohn–Sham linearresponse theory that is given in full form for hybrid functionals in Ref. Misquitta and Stone, 2008a. The expression for also involves the amplitudes defined above, though the matrix elements multiplying it are more complex and are given in full form in Refs. Jeziorski et al., 1993; Patkowski et al., 2012. Analogous expressions exist for the induction energy of monomer B due to the field of A.
In regularized SAPT(DFT) (RSAPT(DFT)) the amplitudes — the coefficients of the firstorder induction wavefunction — are calculated using eq. (12), but this time with the matrix elements on the R.H.S. replaced by those calculated with the regularized electrostatic potential given in eq. (8). Otherwise, the expressions for and are identical with the nonregularized forms. These equations have been implemented in the CamCASP program Misquitta and Stone (2012) and are available as part of a regular SAPT(DFT) calculation of the interaction energy. Note that the expression for used in this paper does not involve scaling Misquitta and Stone (2008a). The scaling expression seems to work by cancellation of errors (of the secondorder terms with corresponding terms in the estimate of the higherorder induction effects Jeziorski et al. (1994)) and cannot be relied upon to give reasonable results for the regularized induction energy.
Iv Numerical details
All calculations of SAPT(DFT) energies and molecular electrostatic and polarization models have been performed with the CamCASP Misquitta and Stone (2012) program. The Kohn–Sham orbital orbitals and orbital energies used in by the CamCASP program were obtained using the DALTON 2.0 Helgaker et al. (2005) program with a patch distributed as part of the SAPT2008 Bukowski et al. (2002) suite of codes.
SAPT(DFT) induction energies were calculated using the PBE0 Perdew et al. (1996); Adamo and Barone (1999) exchangecorrelation functional asymptotically corrected (AC) using the Fermi–Amaldi Fermi and Amaldi (1934) longrange form and the Tozer–Handy splicing function Tozer and Handy (1998). For details see Ref. Misquitta et al., 2005. The secondorder induction energies were evaluated using the hybrid form of the linearresponse kernel Misquitta et al. (2005); Misquitta and Stone (2008a). was evaluated without scaling Misquitta and Stone (2008a) using the expressions in Refs. Hesselmann et al., 2005; Podeszwa et al., 2006b. This was found to be necessary for the regularized induction energies.
Distributed multipoles have been calculated using the GDMA2 Stone (2005) module that is part of the CamCASP suite. Unless otherwise stated, these include terms to rank 4 (hexadecapole moments) on all sites, including the hydrogen atoms. Distributed (anisotropic) polarizabilities have been calculated using the Williams–Stone–Misquitta (WSM) method Misquitta and Stone (2008a); Misquitta et al. (2008a) with terms to rank 3 on all sites except for the hydrogen atoms for which these were restricted to rank 1. Molecular multipole and polarizability models have been calculated using a daugccpVTZ basis with the PBE0/AC density and density response functions. Model energy evaluations were performed using the Orient Stone et al. (2013) program.
V Results
v.1 Determining the regularization parameter
Regularization involves the parameter that has the units of . Equivalently, has the dimensions of length. From fig. 2 and the discussion in Sec. III we know that this lengthscale will be associated with the width of the screened nuclear potential, which, in principle, will be dependent on the atom type and bonding environment. We are faced with two choices: either the parameter is only weakly dependent on atom type, in which case a fixed value may be used for all calculations, or this parameter exhibits a strong atomtype dependence, in which case the regularization procedure will be potentially cumbersome to apply in practice. The first and most important question is how are we to determine the appropriate value of ?
In principle, could be obtained by examining the screened nuclear potential and determining the appropriate regularization needed to ‘fill’ it in to as to prevent all tunneling states. But it is at present unclear what would constitute sufficient ‘filling’, consequently we have instead opted for two alternate procedures:

We could determine by observing that at the optimal regularization there will be no chargetransfer but the polarization component of the induction energy will be left unchanged. Consequently, all basis sets large enough to describe the pure polarization effect would result in the same regularized induction energy at all separations.

Alternatively, we could determine by requiring that all the induction energy in rare gas dimers is chargetransfer.
Consider the first proposal: Our premise here is that the true polarization energy can be described by a reasonably large monomer basis, that is, without the need of basis functions located on the partner monomer. Consequently, with the correct regularization, i.e., the value of that completely suppresses all chargetransfer tunneling effects without altering the longrange part of the electrostatic potential responsible for the polarization, all basis sets capable of describing the true polarization energy will yield the same, regularized, induction energy. If on the other hand, the regularization is insufficient, then some degree of chargetransfer will be allowed. This will lead to a spread in energies as the larger basis sets will be better able to describe the residual tunneling. Finally, if the regularization is excessive, the longrange part of the potential will be affected. This will also cause a variation in energies calculated with different basis sets, with the larger, more diffuse basis sets being more affected.
In fig. 3 we display regularized secondorder induction energies for the water dimer calculated with various basis sets and three values of and a.u. For the regularization parameter a.u. all but the smallest basis set results in the same energy over the entire range of intermolecular separations. The one exception is the augccpVDZ basis in the monomercentered type which is unable to adequately describe even the polarization component of the induction energy. For larger and smaller values of the spread in the results from the six basis sets is seen to increase. To ensure that this observation is not specific to the OH contact in the water dimer, in fig. 4 we compare regularized induction energies for 414 water dimers: 400 of these were chosen using the pseudorandom algorithm described in Ref. Misquitta et al., 2008b and the remaining 14 are those already shown in the above figures. It is remarkable that with a.u. the regularized induction energies calculated with the augccpVQZ MC basis and augccpVQZ MC+ basis are nearly perfectly in agreement for all 414 dimers. Varying to 2.0 or 4.0 a.u. results in a significantly less correlation between the two sets of energies, not just for the hydrogenbonded dimer configurations with strong chargetransfer, but also for configurations with OO contacts which exhibit small chargetransfer energies. We point out here that with this procedure we are unable to distinguish between in the range a.u. However, it seems fairly clear that the optimum value is very close to a.u.
In fig. 5 we display chargetransfer energies for the HF dimer in its minimum energy hydrogenbonded configuration Peterson and T. H. Dunning (1995) calculated with various basis sets. Only results for a.u. are presented as the behaviour of this system is largely similar to that of the water dimer: here too, a.u. is a good choice for the regularization parameter, though the regularized induction energies from the larger basis sets are not in as good agreement at short separations at which a larger value of may be more appropriate. This may indicate that should depend on the type of atom. Although this discrepancy shows up at separations not easily accessed, this needs further investigation. In fig. 6 we compare regularized induction energies calculated using a few values of for 120 pseudorandom dimers and 39 dimers at specific orientations. We see a good agreement of regularized induction energies calculated with the augccpVQZ MC and MC+ basis sets for a regularization of a.u.
Now consider the second proposal: that we determine the optimum value of the regularization parameter by requiring that all the induction energy in a rare gas dimer is chargetransfer. There are a few issues with this proposal. The total induction energy of the rare gas dimers is nearly zero at and around the equilibrium separation, consequently we need to use very short separations to obtain an appreciable total induction energy. This leads to conundrum: while we can argue that the multipole expansion should not result in any polarization energy for these dimers, at short separations, the expansion is not valid as the chargedensities penetrate. It is worth bearing this in mind in the following discussion.
The induction energy of the argon dimer, like that of all raregas dimers, is almost zero at the equilibrium separation Patkowski et al. (2005) of 7.13 a.u., with the polarization and exchange components almost completely cancelling. However, as shown in fig. 7, for very small interatomic separations this is no longer the case, and we get an exponentially varying induction energy. However, with regularization parameter a.u. the regularized induction energy is close to zero for all separations. Using this value of we have calculated secondorder chargetransfer energies using the augccpVZ, =D,T,Q basis sets. From fig. 8 we see that is insensitive to basis set, but, as with the water dimer example, shows considerable basis variation, with the chargetransfer energy tending to zero as the basis set increases. Unlike the first proposal, the optimum value of obtained in this manner is dependent on the type of system: it is 3.2 a.u. for the ArNe complex, and 2.9 a.u. for the ArHe complex. Furthermore, for the neon and helium dimers it varies from 4.0 to 5.0 a.u. though, in these cases, the total induction energies are considerably smaller even at very short separations, so there is a associated ambiguity in the choice of .
Even if we set aside the second proposal, the first method provides compelling evidence that the regularization parameter a.u. is appropriate, at least for the lighter atoms. This value corresponds to a regularization lengthscale of 0.577 Bohr. As mentioned above, we see some indication that should vary with the type of atom, but numerical evidence suggests that this variation is likely to be small, and is probably manifest at very short separations only. In the rest of this article will be computed through eq. (7) with the regularization parameter a.u.
v.2 Analysis of the secondorder chargetransfer energy
Having determined the manner in which the regularization is to be performed, we will now analyse the energies in detail. In fig. 9 we display and energies for the water dimer. The energies have been calculated using (7) with both terms on the RHS calculated using the MC+ type of basis and therefore show very little sensitivity to the choice of basis set (as long as augmented double or better in quality). Contrast this with the strong basis sensitivity of the energies. There are a few features of these energies worth highlighting: at longrange, the energies are similar to those of with the augccpVQZ basis set. But at shortrange, is closer to with the much smaller augccpVDZ basis set. This behaviour supports the discussion in Sec. II.1: The definition does result in accurate secondorder chargetransfer energies with the augccpVQZ basis set at longrange, but as the intermolecular separation decreases, this large basis set starts to pick up chargetransfer excitations leading to an underestimation of the chargetransfer defined via eq. (1). At short separations, we should expect a smaller basis set, in this case, the augccpVDZ basis, to yield more accurate energies. This is indeed what our definition demonstrates.
In fig. 10 we present the chargetransfer energies calculated in the augccpVQZ basis set on a semilog scale. The chargetransfer is usually thought of as being decaying exponentially with separation, so on this plot it should appear as a straight line. However, this is not the case for either of the methods. The energy does exhibit an exponential behaviour at large intermolecular separations, but becomes subexponential at short separations. This is a consequence of the problem discussed in the above paragraph: the will always result in too little CT at short range, and this effect is larger for the larger basis sets. By contrast, is superexponential at short separations: it exhibits a clear double (possibly even multiple) exponential behaviour. This should be expected. If we accept that the chargetransfer process is a tunneling phenomenon, then we should also expect to see contributions from tunneling into each of the (screened) nuclear wells. The dominant process is expected to be the chargedensity of the electron donor (oxygen) tunneling into the poorly screened nuclear potential of the electron acceptor (hydrogen). (In this paper, we use the terms ‘donor’ and ‘acceptor’ to refer to electrons and not protons.) There will also be tunneling of the acceptor density into the wellscreened nuclear potential of the donor oxygen. However, we may additionally expect weaker processes such as the donor density tunneling into the acceptor oxygen. Each process will be approximately exponential, with the barrier height and width determining the exponent. Perturbation theory allows us make the donor to acceptor and acceptor to donor decomposition. These results are displayed in fig. 10. The decomposition of is qualitatively different from that of . In the former we see two processes both apparently with the same exponent (except at shortrange). These do not reflect the tunneling processes described above. However, exhibits at least two exponential processes: The acceptor to donor energy has a larger exponent and decays rapidly with separation. This is what we would expect as tunneling into the wellscreened donor oxygen potential must proceed through a large barrier, leading to a large exponent.
The regularization procedure allows us to analyse the chargetransfer process in even more detail as we can selectively regularize the nuclear potential of the acceptor hydrogen atoms and isolate the chargetransfer contribution from the donor to the oxygen of the acceptor. From fig. 10 we see that, as might be expected, this is a much weaker process. What is unusual is the relatively small exponent associated with this process. The reason for this is as yet unclear.
In figs. 11 and 12 we present a similar analysis of the secondorder chargetransfer energy for the HF dimer in its hydrogenbonded orientation. The features we have highlighted for the water dimer are also seen here. As with the water dimer, is largely basisindependent (the augccpVTZ/MC+ and augccpVQZ/MC+ results are nearly identical, but the augccpVDZ/MC+ values differ) and interpolates between the values calculated using the augccpVDZ basis (at shortrange) and the augccpVQZ basis (at longrange). Also, as for the water dimer, the acceptor to donor secondorder chargetransfer energy becomes appreciable only for short dimer separations (less than 5 a.u.). This is quite different from the energy for which we see acceptor to donor contributions at all separations.
/Bohr  Khaliullin et al. Ref.Khaliullin et al.,2008  
Water dimer in hydrogenbonded orientation  
—  —  
—  —  
HF dimer in hydrogenbonded orientation  
—  —  
FHCO (linear)  
—  —  
FHOC (linear)  
—  —  
HBCO : augccPVTZ/DC : B3LYP optimized (linear)  
HBNH : augccPVTZ/DC : B3LYP optimized  
Pyridine dimer: augccPVTZ/MC+ : double Hbonded geometry  
—  — 
Table 1 shows numerical values of the secondorder induction and chargetransfer energies for various dimers. The results for the water and HF dimers quantify what is already displayed in the above figures: at the selected geometries, and agree reasonably well; but the differences get larger at shorter separations. With the exception of the HB complexes, the chargetransfer is mainly from the donor to the acceptor; the reverse process (acceptor to donor) is much weaker. This is particularly true for ; although the donor to acceptor energy is dominant for , the acceptor to donor energies can be significantly larger than those from . For the mixed HF and CO system, chargetransfer from the C to H in FHCO is nearly 3 times as large as that from O to H in FHOC, consistent with the strong donor property of C in CO.
The HB complexes with CO and NH are interesting as both of these have very short separations. These separations are short enough that one may question the use of an intermolecular perturbation theory like SAPT(DFT). Perhaps surprisingly, SAPT(DFT) interaction energies are within 5% of CCSD(T) energies for both complexes, with the agreement between the two best for the HBNH complex. The differences between and are quite large for both complexes. As discussed in sec. II.1, the CTSM09 definition should be expected to recover an ever smaller fraction of the true secondorder CT as the intermolecular separation decreases. This seems to be the case here. While correctly identifies the donor (NH, CO) to acceptor (HB) chargetransfer as the larger quantity, the actual amount of CT from this method is substantially smaller than both and the infiniteorder results of Khaliullin et al. Khaliullin et al. (2008). Interestingly, while is consistent with the Khaliullin results for the HBNH complex, the two sets of results differ qualitatively for the HBCO complex. Here we would expect transfer from the CO to the HB to dominate. This is given by , but the ALMO method shows the opposite result.
Finally, the pyridine dimer is an interesting case as it exhibits a double hyhrogen bond between each of the NH pairs in the symmetry, planar complex. The donor (N) to acceptor (H) bond length is the longest of the complexes discussed above. This leads to relatively small total induction energies and even smaller chargetransfer energies. The latter are just over a tenth of the total induction energy—nearly an order of magnitude smaller than the chargetransfer energies in the other complexes. In this case the electron delocalisation process works symmetrically in both directions, consequently there is no nett charge transfered between the two pyridine molecules.
v.3 Polarization models
Now that we have demonstrated both the numerical stability and the physical nature of , we are in a position to use this definition to determine the damping needed in polarization models. The basic idea is to fit the damping parameters so that the (noniterated) classical polarization model matches the true secondorder polarization energies for a large number of dimers. In the following examples the classical polarization model was constructed using a rank 3 (3 on the heavy atoms and 1 on the hydrogen atoms) WSM distributed polarizability model Misquitta and Stone (2008a); Misquitta et al. (2008a) and a rank 4 (all atoms) distributed multipole model Stone (2005). The models are damped using the Tang–Toennies damping functions Tang and Toennies (1984) with isotropic (possibly sitepairdependent) damping parameters. For numerical details see Sec. IV. Note that the damping coefficients obtained in this section depend on the details of the model, consequently a direct comparison to other attempts to determine the damping coefficients Sebetci and Beran (2010) cannot be made.
Fig. 13 shows three polarization models for the water dimer in its hydrogenbonded orientation. These models differ in their damping parameters only. The first uses a damping parameter derived from the expression
(13) 
which was given by Misquitta and Stone Misquitta and Stone (2008a) and is derived from the molecular vertical ionization energies and . In a later paper on dispersion models Misquitta and Stone (2008b), these authors demonstrated that this simple expression resulted in accurate damping models for a variety of systems, ranging from hydrogenbonded complexes to van der Waals systems in a wide range of orientations. For the water dimer, using a vertical ionization energy of 0.4638 a.u. Lias , using eq. (13) we get a.u. From fig. 13 we see that, as demonstrated in Ref. Misquitta and Stone (2008a), the resulting polarization model agrees well with the total, unregularized energy . The agreement is particularly good at the dimer equilibrium separation, though at shorter separations this model tends to underestimate . Notwithstanding this seemingly good performance, there is considerable evidence from Sebetci and Beran Sebetci and Beran (2010) and also from our own numerical experiments that polarization models derived using eq. (13) significantly overestimate the manybody polarization energy. By fitting to the manybody energies of clusters of water molecules, Sebetci and Beran obtain an optimized parameter of a.u., i.e., the damping needs to be considerably enhanced.
The reason for this is that the polarization model should reproduce the true polarization energy and not . If the model is derived to match , it will implicitly include, via the damping, some amount of the twobody chargetransfer energy. While this is, in itself, not a significant problem for the twobody energy, it can lead to large errors in the manybody polarization energy which will then include spurious twobody chargetransfer effects. This would lead to the polarization overbinding of clusters of molecules seen by Sebetci and Beran and discussed in Secs. I and II of this article.
As discussed above, the solution to this problem is to determine the damping of the polarization model by fitting it to reproduce the true polarization energy only. If we do this using a single damping parameter, that is, the damping parameter is independent of type of sites, we get a.u. From fig. 13 we see that this model results in a very good match with , at least for the dimer in its hydrogenbonded orientation. However, as can be seen in fig. 14, the agreement is not as good for other orientations; in particular, at those with close OO separations the polarization model tends to overestimate . A detailed examination of the water dimer polarization energies at orientations with close HH and OO contacts suggests that the damping parameter is strongly dependent on the types of the sites in the interacting pair. A much better fit to is obtained with a three parameter damping model, , in which , and a.u. Of these, is not very well defined and can be fixed to a relatively wide range of values. From figs. 13 and 14 we see that this model is significantly better than the simpler, oneparameter model with a.u.
Even the threeparameter model exhibits somewhat larger errors for the dimers with close OO contacts (highlighted in fig. 14). It is possible that a proper nonlinear optimization of the damping model will result in a model that removes these discrepancies, but it is also possible that at least some of the residual error is from the lack of angular dependence in the damping model. These issues are currently under investigating.
Although it may appear that the singleparameter model with a.u. matches the Sebetci and Beran value of a.u., matters are not as straightforward. First of all, Sebetci and Beran used a somewhat different multipole and polarizability models: their distributed multipole model had terms on the hydrogen atoms limited to rank 1, and in their WSM polarizability model the maximum rank was 2. Both of these, particularly the former, result in smaller polarization energies. Consequently, the damping need not be as large. Indeed, using a multipole model similar to theirs and by fitting to hydrogenbonded dimer orientations only, we obtain a single parameter damping parameter of 1.6 a.u. (a larger means less damping). But these orientations are not fully representative of those found in water clusters such as those used by Sebetci and Beran. Here the contacts are more diverse, particularly in the more compact water clusters. We conjecture that the Sebetci and Beran value of a.u. is a compromise that is the average of the sitepairdependent parameters in model described above. While this needs to be confirmed, we emphasise that all of these results are consistent: the polarization damping should be much smaller than what we would expect from eq. (13).
The HF dimer exhibits many of these results. The polarization model based on eq. (13) and a vertical ionization potential of 0.5669 a.u.Lias results in a damping parameter of 2.2 a.u. From fig. 15 we see that this model closely matches ; in this case the match is even better than for the water dimer. Our second model is obtained by fitting to the true polarization energy . This yields a single parameter model with a.u. which, as we see from fig. 16, is able to match the energies for other dimer orientations, albeit with a larger scatter than we had for the water dimer. At least some of this scatter is due to the unusual shortrange mismatch in and the polarization model seen in fig. 15. But there is also evidence that, as with the water dimer, an accurate polarization model would require three different damping parameters, with the FF interaction needing a considerably stronger damping. We have yet to explore this possibility fully.
Curiously, Sebetci and Beran Sebetci and Beran (2010) found the original damping model with a.u. to be suitable for cyclic clusters of HF molecules. Exploration of this issue shows that with a multipole model similar to the one they used (with terms limited to rank 1 on the hydrogen atoms) we obtain an optimized damping parameter of 2.0 a.u. by fitting to . This parameter should increase still more when the maximum rank of the polarization model is reduced to 2 from the maximum of 3 which we have used. Therefore, in this case too, we believe that our results are fully consistent with those from Sebetci and Beran and indicate that we are indeed now able to derive manybody polarization models from the dimer alone.
Vi Discussion
We have presented a definition of the chargetransfer energy that is based on interpreting the transfer of charge between molecules through the process of tunneling. In this view, as illustrated in fig. 2, intermolecular charge transfer occurs by electron density tunneling into the screened nuclear potential of the partner monomer. This viewpoint leads to a simple way of defining the chargetransfer energy by regularizing the screened nuclear potential so as to suppress tunneling, while still allowing classical electromagnetic polarization.
The Gaussiantype regularization we have used involves a parameter that has the dimensions , or, equivalently, the regularization introduces a nuclear screening lengthscale . This is the single most important parameter in this procedure. A priori the screening lengthscale is expected to depend on the atom type, but this dependence has been shown to be weak. Using a number of systems and two different procedures, we have demonstrated that the regularization parameter, is very close to 3.0 a.u., that is, the regularization occurs on a lengthscale of 0.58 Bohr. This value of the parameter is suitable for systems involving lighter elements, but further work is needed to understand why this value is appropriate and how it may be expected to change for systems containing heavier elements.
Once the value of the regularization parameter has been fixed, the secondorder chargetransfer energy, , has a welldefined completebasisset limit for all intermolecular separations. This strongly contrasts with the definition put forward by Stone and Misquitta Stone and Misquitta (2009) which exhibits a strong basis dependence, particularly at short separations.
With the proposed definition of the chargetransfer through regularization the chargetransfer energy for the water dimer and HF dimer systems exhibits a clear double exponential character. This has been shown to stem from the strong asymmetry in the electron donor to acceptor and electron acceptor to donor components of the energies: the donor to acceptor process is not only dominant, but decays more slowly with separation; indicative of a tunneling process through a smaller barrier. Additionally, using a partial regularization we have been able to decompose the donor to acceptor chargetransfer energy into a primary process from the electronrich donor O or F atom into the electron deficient acceptor H atom, and a much smaller secondary process from the donor O or F into the acceptor O or F. The secondary process is an order of magnitude smaller than the primary donor to acceptor process.
This new definition of the secondorder chargetransfer energy has been used on a variety of systems, including some with very strong hydrogen bonds (HBCO and HBNH) and one doubly hydrogenbonded system (the pyridine dimer in its planar configuration). In all cases the computed chargetransfer energy makes physical sense. The differences between and other definitions such as the Stone and Misquitta and ALMO method from Khaliullin et al. Khaliullin et al. (2008) are particularly dramatic for the most strongly bound HB complexes for which only results in physically meaningful energies. In the doubly hydrogenbonded pyridine example, while there is a stabilisation due to the tunneling, due to symmetry there is no nett charge transferred between the two molecules: they remain neutral. So, in a sense, the term ‘chargetransfer’ is a misnomer.
Finally, we have used regularization to suppress all chargetransfer contributions from the secondorder induction energy thereby defining a pure secondorder polarization energy against which we have fitted the damping in the WSM polarization models. Comparisons against data from Sebetci and Beran Sebetci and Beran (2010) indicates that the resulting models are able to describe the manybody polarization energies accurately. This is a major step forward in our ability to model the major part of the manybody nonadditive energy in polarizable systems from calculations on the dimer alone.
Based on the arguments put forward in the Introduction and
the many examples provided herein we now propose a new definition of the
chargetransfer energy:
The process of charge transfer can be thought of as a charge
delocalisation through electron tunneling into the screened nuclear
potentials of the partner monomers.
The resulting energy of stabilization is the
intermolecular chargetransfer or delocalisation energy.
We suggest that it may be conceptually advantageous to term this
the delocalisation energy and reserve the term ‘chargetransfer’
for the phenomenon of chargetransfer excitations studied by
Mulliken Mulliken (1952).
Though the interpretation of the CT as a tunneling process may appear different from the usual picture of the CT as an excitation from a donor orbital to an acceptor located on the partner monomer, we believe that these two viewpoints are consistent. The secondorder RayleighSchrödinger perturbation expression for the secondorder induction energy given in eq. (2) contains terms that arise from excited states that have a significant component at the sites of the partner monomer. For these states to make a significant contribution to the induction energy, the matrix elements in the numerator need to be significant (and the energy differences in the denominator small). One way this can happen is when the excited state has significant contributions in the regions of the nuclei of the partners where the screened potential is significant. This is just another way of describing tunneling states.
Chargetransfer as incipient chemical bonding: if we accept the physical picture of charge transfer as a tunnelling of electrons into the screened nuclear potentials of the partner monomer, then we must also accept the view the wellestablished view that chargetransfer is a form of incipient covalent bonding.
There are several advantages to the proposed definition of the chargetransfer energy:

It leads to physically appealing definition of CT.

The resulting (electron) acceptor to donor and (electron) donor to acceptor contributions make physical sense.

It leads to polarization models that are applicable to manybody systems although they are based on dimer calculations alone.

The method is implemented as part of regular SAPT(DFT) calculation.

The results are independent of basis set and true basisset convergence can be achieved.

The method is independent of the type of basis set: it could be used, for example, with planewave basis sets. Therefore it should be possible to use this definition to calculate the chargetransfer, or delocalisation energy, in a variety of electronic structure programs.
Amongst the issues that need resolving is our incomplete understanding of the origin of the regularization lengthscale and how it depends on atom type, and the manner in which the present procedure needs to be extended to calculate chargetransfer effects beyond second order in perturbation theory. Both of these issues are under current investigation.
Vii Acknowledgements
The author would like to acknowledge Prof Anthony J Stone and Prof Kenneth Jordan for many stimulating and fruitful discussions, and in particular, Prof Stone for many useful comments and suggestions and a stimulating collaboration. Additionally, the author thanks the School of Physics at Queen Mary, University of London for support and the Thomas Young Centre for providing a stimulating research platform in the London area.
References
 Stone (1993) A. J. Stone, Chem. Phys. Lett. 211, 101 (1993).
 Stone and Misquitta (2009) A. J. Stone and A. J. Misquitta, Chem. Phys. Lett. 473, 201 (2009).
 Khaliullin et al. (2008) R. Z. Khaliullin, A. T. Bell, and M. HeadGordon, J. Chem. Phys. 128, 184112 (2008).
 Cohen et al. (2008) A. J. Cohen, P. MoriSÃ¡nchez, and W. Yang, Science 321, 792 (2008), eprint http://www.sciencemag.org/content/321/5890/792.full.pdf, URL http://www.sciencemag.org/content/321/5890/792.abstract.
 Misquitta and Stone (2008a) A. J. Misquitta and A. J. Stone, J. Chem. Theory Comput. 4, 7 (2008a).
 Williams et al. (1995) H. L. Williams, E. M. Mas, K. Szalewicz, and B. Jeziorski, J. Chem. Phys. 103, 7374 (1995).
 Sebetci and Beran (2010) A. Sebetci and G. J. O. Beran, J. Chem. Theory Comput. 6, 155 (2010).
 Schenter and Glendening (1996) G. K. Schenter and E. D. Glendening, J. Phys. Chem. 100, 17152 (1996).
 Azar et al. (2013) R. J. Azar, P. R. Horn, E. J. Sundstrom, and M. HeadGordon, J. Chem. Phys. 138, 084102 (2013).
 Jeziorski et al. (1994) B. Jeziorski, R. Moszynski, and K. Szalewicz, Chem. Rev. 94, 1887 (1994).
 Patkowski et al. (2001a) K. Patkowski, T. Korona, and B. Jeziorski, J. Chem. Phys. 115, 1137 (2001a).
 Patkowski et al. (2004) K. Patkowski, B. Jeziorski, and K. Szalewicz, J. Chem. Phys. 120, 6849 (2004).
 Claverie (1971) P. Claverie, Int. J. Quantum Chem. 5, 273 (1971), ISSN 1097461X, URL http://dx.doi.org/10.1002/qua.560050304.
 Kutzelnigg (1980) W. Kutzelnigg, J. Chem. Phys. 73, 343(17) (1980).
 Morgan III and Simon (1980) J. D. Morgan III and B. Simon, Int. J. Quantum Chem. 17, 1143 (1980).
 Adams (1999) W. H. Adams, Int. J. Quantum Chem. 72, 393 (1999).
 Herring (1962) C. Herring, Rev. Mod. Phys. 34, 631 (1962), URL http://link.aps.org/doi/10.1103/RevModPhys.34.631.
 Patkowski et al. (2001b) K. Patkowski, B. Jeziorski, and K. Szalewicz, J. Mol. Struct. (Theochem) 547, 293 (2001b).
 Adams (2002) W. H. Adams, Theor. Chim. Acta 108, 225 (2002), ISSN 1432881X, URL http://dx.doi.org/10.1007/s0021400203773.
 Jeziorski and Szalewicz (1998) B. Jeziorski and K. Szalewicz, Encyclopedia of Computational Chemistry (Wiley, Chichester, UK, 1998), vol. 2, p. 1376.
 Jeziorski and Szalewicz (2002) B. Jeziorski and K. Szalewicz, in Handbook of Molecular Physics and Quantum Chemistry, edited by S. Wilson (Wiley, 2002), vol. 8, chap. 8, pp. 37–83.
 Misquitta et al. (2005) A. J. Misquitta, R. Podeszwa, B. Jeziorski, and K. Szalewicz, J. Chem. Phys. 123, 214103 (2005).
 Podeszwa et al. (2006a) P. Podeszwa, R. Bukowski, and K. Szalewicz, J. Phys. Chem. A 110, 10345 (2006a).
 Bukowski et al. (2006) R. Bukowski, K. Szalewicz, G. Groenenboom, and A. van der Avoird, J. Chem. Phys. 125, 044301 (2006).
 Patkowski et al. (2012) K. Patkowski, K. Szalewicz, and B. Jeziorski (2012), private communication.
 Casida (1995) M. E. Casida, in Recent Advances in DensityFunctional Theory, edited by D. P. Chong (World Scientific, 1995), p. 155.
 Colwell et al. (1995) S. M. Colwell, N. C. Handy, and A. M. Lee, Phys. Rev. A 53, 1316 (1995).
 Jeziorski et al. (1993) B. Jeziorski, R. Moszynski, A. Ratkiewicz, S. Rybak, K. Szalewicz, and H. Williams, SAPT: A program for manybody symmetryadapted perturbation theory calculations of intermolecular interaction energies (STEF, Cagliari, 1993), vol. B, p. 79.
 Misquitta and Stone (2012) A. J. Misquitta and A. J. Stone, CamCASP: a program for studying intermolecular interactions and for the calculation of molecular properties in distributed form, University of Cambridge (2012), http://wwwstone.ch.cam.ac.uk/programs.html#CamCASP.
 Helgaker et al. (2005) T. Helgaker, H. J. A. Jensen, P. Joergensen, J. Olsen, K. Ruud, H. Aagren, A. Auer, K. Bak, V. Bakken, O. Christiansen, et al., Dalton, a molecular electronic structure program, release 2.0 (2005), see http://www.kjemi.uio.no/software/dalton/dalton.html.
 Bukowski et al. (2002) R. Bukowski, W. Cencek, P. Jankowski, B. Jeziorski, M. Jeziorska, V. Lotrich, S. Kucharski, A. J. Misquitta, R. Moszynski, K. Patkowski, et al., SAPT2008: an ab initio program for manybody symmetryadapted perturbation theory calculations of intermolecular interaction energies, University of Delaware and University of Warsaw (2002), URL http://www.physics.udel.edu/~szalewic/.
 Perdew et al. (1996) J. Perdew, K. Burke, and M. Ernzerhof, Phys. Rev. Lett. 77, 3865 (1996), also see erratum, Phys. Rev. Lett. 78, 1396 (1996).
 Adamo and Barone (1999) C. Adamo and V. Barone, J. Chem. Phys. 110, 6158 (1999).
 Fermi and Amaldi (1934) E. Fermi and G. Amaldi, Men. R. Acad. Italia 6, 117 (1934).
 Tozer and Handy (1998) D. J. Tozer and N. C. Handy, J. Chem. Phys. 109, 10180 (1998).
 Hesselmann et al. (2005) A. Hesselmann, G. Jansen, and M. Schütz, J. Chem. Phys. 122, 014103 (2005).
 Podeszwa et al. (2006b) R. Podeszwa, R. Bukowski, and K. Szalewicz, J. Chem. Theory Comput. 2, 400 (2006b).
 Stone (2005) A. J. Stone, J. Chem. Theory Comput. 1, 1128 (2005).
 Misquitta et al. (2008a) A. J. Misquitta, A. J. Stone, and S. L. Price, J. Chem. Theory Comput. 4, 19 (2008a).
 Stone et al. (2013) A. J. Stone, A. Dullweber, O. Engkvist, E. Fraschini, M. P. Hodges, A. W. Meredith, D. R. Nutt, P. L. A. Popelier, and D. J. Wales, Orient: a program for studying interactions between molecules, version 4.8, University of Cambridge (2013), http://wwwstone.ch.cam.ac.uk/programs.html#Orient.
 Misquitta et al. (2008b) A. J. Misquitta, G. W. A. Welch, A. J. Stone, and S. L. Price, Chem. Phys. Lett. 456, 105 (2008b).
 Peterson and T. H. Dunning (1995) K. A. Peterson and J. T. H. Dunning, J. Chem. Phys. 102, 2032 (1995).
 Patkowski et al. (2005) K. Patkowski, G. Murdachaew, C.M. Fou, and K. Szalewicz, Mol. Phys. 103, 2031 (2005), eprint http://www.tandfonline.com/doi/pdf/10.1080/00268970500130241, URL http://www.tandfonline.com/doi/abs/10.1080/00268970500130241.
 Tang and Toennies (1984) K. T. Tang and J. P. Toennies, J. Chem. Phys. 80, 3726 (1984).
 Misquitta and Stone (2008b) A. J. Misquitta and A. J. Stone, Mol. Phys. 106, 1631 (2008b).
 (46) S. G. Lias, Ionization energy evaluation in nist chemistry webbook, nist standard reference database number 69, eds. w. g. mallard and p. j. linstrom, gaithersburg, 2000 (http://webbook.nist.gov).
 Mulliken (1952) R. S. Mulliken, J. Am. Chem. Soc. 74, 811 (1952).