Empirical Self-Consistent Correction for the Description of Hydrogen

Article pubs.acs.org/JCTC

Empirical Self-Consistent Correction for the Description of Hydrogen Bonds in DFTB3 Jan Ř ezác*̌ Institute of Organic Chemistry and Biochemistry, Czech Academy of Sciences, 166 10 Prague, Czech Republic S Supporting Information *

ABSTRACT: The description of hydrogen bonds in the density-functional tightbinding (DFTB) method continues to be a challenging task because the approximations that make the method computationally efficient compromise already the first-order electrostatic contribution to the interaction. So far, the best results have been achieved with fully empirical corrections such as the recently reparametrized DFTB3-D3H4 method. This approach has, however, important limitations that arise from its independence of the actual electronic structure. Here, we present a novel correction denoted as D3H5, which is integrated deeper in the DFTB method, correcting the problem at the place of its origin. It is applied within the self-consistent evaluation of electrostatic interactions, where it empirically models the missing contributions of atomic multipoles and polarization. Despite being very simple and using fewer parameters than D3H4, it is both more accurate and more robust. In data sets of small model systems, it yields errors below 1 kcal/mol, and it performs comparably well in larger systems. Unlike D3H4, it can describe cooperativity in H-bond networks, which makes it more transferable to more complex systems.

1. INTRODUCTION Recent improvements in the description of noncovalent interactions in semiempirical quantum-mechanical methods (SQM, also comprising the density-functional tight-binding method, DFTB) enabled their application to large, complex systems including biomolecules (for a recent review, see ref 1). While this field is still dominated by molecular mechanics, SQM methods provide several advantages. Not only can they describe a wide range of effects of quantum-mechanical origin that are difficult to capture in classical force fields, but they are also more practical because they do not require any systemspecific parametrization. SQM calculations are, of course, more demanding than molecular mechanics, which limits their use to cases where extensive sampling of configurational space is not necessary. One of the successful applications of modern SQM methods is computer-aided drug design, where they may be used for large-scale calculations of protein−ligand interactions.2−4 The approximations involved in all the SQM methods lead to a rather poor description of noncovalent interactions unless specific corrections are introduced. Like in all mean-field QM methods, London dispersion is missing, but it can be added via an empirical correction.5−11 The second important issue is the description of hydrogen bonds, for which various corrections have been introduced.12−17 We have recently reviewed the existing approaches to hydrogen bonding in the self-consistentcharge density-functional tight-binding (SCC-DFTB, abbreviated here as DFTB) method.18 Using an interaction energy decomposition, we have shown that the most severe errors are already introduced at the level of first-order electrostatics, where the monopole approximation inherent to the DFTB © XXXX American Chemical Society

methodology leads to a substantial underestimation of this component of interaction energy. The proper solution would be to extend the DFTB formalism to atomic multipoles. This approach is already under development, but it will take some time before it is ready for practical applications. Even then, it may not be the final solution to this issue, because e.g. all the MNDO-based SQM methods use atomic multipoles, but they still underestimate H-bonds badly. This has motivated our search for a new method that would combine the accuracy and simplicity of fully empirical corrections with the robustness and transferability of those that are more physics-based. The new correction proposed here outperforms the best method available so far, the recent reparametrization18 of the D3H4 correction,15 while being simpler and less empirical. Since it enters the self-consistent treatment of the electrostatics in DFTB, it properly captures many-body effects, which leads to a significant improvement in the description of the dense networks of hydrogen bonds such as in liquid water. Hydrogen-bond correction has been parametrized for the DFTB3 version19 of DFTB with the 3OB parameter set20−22 and D3 dispersion;15,18,23 the resulting method has been denoted DFTB3-D3H5. This paper presents its parametrization for the hydrogen bonds of O, N, and S. The performance of the method has been thoroughly tested on multiple data sets. While the parametrization is, in principle, easily extensible, various issues in the DFTB parameters for heavier elements make it more complicated, as demonstrated on the example of sulfur. These issues become more difficult to Received: June 19, 2017

A

DOI: 10.1021/acs.jctc.7b00629 J. Chem. Theory Comput. XXXX, XXX, XXX−XXX

Article

Journal of Chemical Theory and Computation solve in halogens, where the description of hydrogen bonds is in conflict with the repulsive correction introduced to correct the description of halogen bonds.22 An accurate description of halogen and hydrogen bonding at the same time requires a modification of the DFTB parameters; this extension, which will bring an equally accurate description of all interactions of halogens, is under development.

γab =

γabH =

(1)

1 2

∑ ΔqaΔqbγab ab

(2)

where Δqa and Δqb are self-consistent atomic charges relative to the neutral atom reference. The atomic electron densities are modeled with Slater-type spherical functions, and the function γab is thus a Coulomb integral over these densities. At longrange, the atoms can be approximated well by point charges, and γab converges to 1/rab. In the other limit, when rab approaches zero, the gamma function converges to the onsite self-repulsion energy, which enters the calculation as the Hubbard parameter Ua (which equals twice the chemical hardness). The latter limit imposes the relation between the exponent of the Slater charge density τa (which is closely related to the covalent radius of the atom) and Ua:

τa =

16 Ua 5

⎡ ⎛ U + U ⎞ζ ⎤ 1 b⎟ 2 r ⎥ − S(rab , Ua , Ub) × exp⎢ −⎜ a ⎢⎣ ⎝ 2 ⎠ ab⎥⎦ rab

(5)

This modification is called XH-damping or γ-damping. The exponent ζ in the damping function is an empirical parameter. As it has been found that this modification improves interaction energies in hydrogen bonds, the damping of the γ function has become widely adopted as a correction for this purpose. The value of ζ has been adjusted in order to reproduce the interaction energy in the water dimer, and in the DFTB3 method, the value of ζ = 4.0 is used. While it brings a significant improvement, the energies of hydrogen bonds are still underestimated in some cases, especially in H-bonds involving nitrogen. Moreover, when this correction is made stronger, the distance dependence of the interaction energy is changed, and the resulting method predicts too short hydrogen bonds (what is already observed in some systems with the standard parametrization18). This modification of the γ function corrects only one of the problems encountered in hydrogen bonds − the radius of the hydrogen (although fitting the ζ allows some compensation for the other issues empirically). The most important issue is the monopole approximation; it is the use of spherical electron densities for all atoms. In hydrogen bonds, the electron donor interacts via a lone pair which could only be described properly using an atomic dipole and possibly higher-order multipoles. It is possible to extend DFTB in this direction,26 but no generalpurpose method using this formalism is available yet. Another issue is the limited atomic polarizability caused by the use of a minimal basis set. Polarization correction has been introduced recently,17,27,28 but our tests indicate that it cannot be parametrized reliably because the limited polarizability causes only a fraction of the overall error.18 We have been successful with introducing fully empirical corrections for hydrogen bonds, and the recent reparametrization18 of the H4 correction15 outperforms all the other variants of DFTB in the description of hydrogen bonds. Nevertheless, this approach has its limitations too; since the correction is independent of the electronic structure, it does not transfer well between e.g. neutral and charged species and underestimates the many-body effects necessary for a proper description of larger systems. 2.2. D3H5 Correction. Our experience with the application of DFTB to noncovalent interactions and careful evaluation of existing approaches18 has led us to the proposal of a novel method which is a synthesis of the two most successful ones. Since it would be very complex to correct all the issues encountered in hydrogen bonds by removing the approx-

The first term, EH0, is the tight-binding electronic energy calculated using Hamiltonian and overlap matrix elements precalculated at the DFT level and tabulated. The last term, Erep, accounts for core−core repulsion and corrects doublecounting in other terms. It is obtained empirically by fitting the overall energy to DFT reference results. The derivation of these terms and details on their calculations are covered well in the literature.24 What is most important for us is the second term, Eγ, which represents the correction arising from the electrostatic interaction of the atoms in the second-order expansion of DFT energy. It reads Eγ =

(4)

where the function S is evaluated analytically (see Supporting Information of ref 19 for formulas). The use of a single parameter Ua for both on-site terms, where it has the meaning of chemical hardness, and in the gamma function, where its relation to the size of the atom (eq 3) is used, is clearly a rough approximation. These two quantities correlate well only in a single period of the periodic table, and especially hydrogen behaves differently from the other elements (see Figure 2 in ref 19). This has been corrected empirically by modifying the gamma functions for all atom pairs involving at least one hydrogen.25 In such a case, the function S is multiplied by damping functions so that the gamma function reads

2. METHODS 2.1. Electrostatics in DFTB. In order to proceed, it is necessary to introduce the relevant part of the DFTB methodology.19,24 The method can be derived from DFT formalism as a Taylor expansion of the DFT energy around a reference state with electron density constructed from neutral atom densities. When the expansion is truncated in the second order, the DFTB energy can be expressed as E = E H0 + E γ + Erep

1 − S(rab , Ua , Ub) rab

(3)

The value of Ua for each element is derived from DFT calculations. To account for changes in the chemical hardness in different charge states, a third-order DFTB expansion can be introduced. An additional parameter is needed to describe an atom, namely the derivative ∂Ua/∂qa, which is also obtained from DFT calculations. It is used in the calculation of the derivatives of the function γab with respect to the atomic charges that enter the expression for the third-order DFTB energy in an additional term, EΓ (see ref 19 for details). This formalism mainly improves the description of charged species. The gamma function can be written in terms of a shortranged correction S to the interaction of point charges as B


Article

Journal of Chemical Theory and Computation

possible to compensate for the missing terms such as on-atom polarization. This correction can be extended straightforwardly to third-order DFTB; unlike the γ-damping (eq 5), the additional term does not depend on U, so that no additional derivatives are needed to calculate the third-order terms. To calculate the gradient of the DFTB energy, a derivative ∂γH5 ab / ∂rab is needed, but it is straightforward to calculate it using the chain rule from the variables already available and a derivative of the correction term. The resulting function (with values of the parameters as listed below) is plotted in Figure 1 in comparison with the

imations that cause them, we keep the correction empirical. On the other hand, we apply it at the level of the γ function so that it becomes an integral part of the SCC procedure. As the next step in the series of corrections developed by us, we have named it H5. It is clear that the interaction of a hydrogen with the electron donor in a hydrogen bond is underestimated in DFTB. The main reasons, discussed above, are the monopole approximation, the limited polarizability of atoms, and an incorrect radius of hydrogen. To correct for these empirically within the present formalism, we should increase the electrostatic interaction in these pairs of atoms. This can be easily achieved by scaling the whole γ function. In order to correct only hydrogen bonds, the correction should only be applied at the distances where H-bonding occurs. To avoid altering the length of the H-bonds (which happens in the common γ-damping approach), the correction should be strongest around the equilibrium distance and decay again at shorter distances. This approach has a physical meaning too − at longer distances, the monopole approximation works well, with the multipolar components and polarization becoming important only in the short-range. The short-ranged scaling of the Coulomb term can thus be viewed as an effective means to account for both of these effects, which decay faster than the electrostatic interaction of monopoles. Using such a correction also for the missing polarization can be supported by further evidence. The contribution of polarization is stronger in systems that are highly polar, namely where already the first-order electrostatics is strong. There may be exceptions to this, e.g. interactions of polar molecules with nonpolar but highly polarizable ones, but this is not the case of hydrogen bonds. To verify this hypothesis, we have used the DFT-SAPT interaction energy decomposition that we have recently reported for the S66x8 data set,29 taking only the 66 equilibrium geometries from it. When the induction energy (calculated as the sum of E2ind and E2exch‑ind) is expressed as a linear function of the first-order Coulomb term E1pol, very good correlation with R2 = 0.96 is obtained (the correlation is even stronger but slightly nonlinear, R2 = 0.99 is obtained with a quadratic fit). This suggests that a correction proportional to the Coulomb interaction provides a very good model also for polarization if the different distance dependence is accounted for. The simplest correction that builds upon these principles is the addition of a Gaussian function, which amplifies the interaction, centered around the equilibrium distance of the hydrogen bond:

Figure 1. Plot of the γ function for an oxygen−hydrogen pair as a function of the interatomic distance. The H5 correction is compared to the γ-damping approach and an unmodified γ function. The vertical line denotes the equilibrium O−H distance of the hydrogen bond in the water dimer.

original gamma function (eq 4) and the one using the standard γ-damping (eq 5). It is clear that the γ-damping gets stronger as the distance decreases which leads to underestimated hydrogen bond lengths. The H5 correction acts only around the equilibrium distance of hydrogen bonds, and the function approaches the unmodified γ again at shorter distances. As a result, the H5 correction could be made stonger without compromising the geometry of the system. The H5 correction has been developed for third-order DFTB with the 3OB parameter set. It is applied in combination with the D3 dispersion correction15 that we have recently reparametrized for this setup.18 Since the correction is only applied to hydrogen bonds, it should be examined how the γdamping affects other hydrogen interactions. The most pronounced issue is the hydrogen−hydrogen repulsion in aliphatic hydrocarbons, where DFTB (with dispersion correction) predicts too short HH contacts and too attractive an interaction. This issue is, however, of not only electrostatic origin. In the neopentane dimer from the S66 data set, where this error is very large, the difference between interaction energies calculated using DFTB3 with and without γ-damping is only 0.1 kcal/mol, while about 0.5 kcal/mol is needed to approach the CCSD(T) benchmark. This difference becomes much larger at distances below the equilibrium. We have addressed this lack of repulsion empirically in the dispersion correction by adding a repulsive term to all H−H pairs,15 and we use this version of the D3 correction here as well. 2.3. Benchmark Data Sets Used. The D3H5 correction was parametrized on the S66x8 data set, a data set of 66

⎛ (r − r )2 ⎞⎞ ⎛1 ⎞ ⎛ γabH 5 = ⎜ − S(rab , Ua , Ub)⎟ × ⎜⎜1 + k × exp⎜ − ab 2 0 ⎟⎟⎟ 2w ⎝ rab ⎠ ⎝ ⎝ ⎠⎠

(6)

Three parameters are introduced: the scaling factor k determines the magnitude of the correction, the distance r0 sets the position of its maximum, and the parameter w affects the width of the Gaussian. Since we can observe large differences between the hydrogen bonds of O and N, all the parameters have been set separately for each element. We tried more complex functions with different shapes as well as asymmetric ones, but no significant improvement that would warrant the added complexity was achieved. Note that the H5 formalism, in contrast to the common γ-damping, allows making the interaction stronger than 1/rab, which makes it C


Article


ones; here, we use the recent update of the data set where the cluster decomposition energies were obtained using a state-ofthe-art composite CSCD(T)/CBS setup.41 Since the only available data set of noncovalent complexes of sulfur had been used for the parametrization, it became necessary to build a new validation set for testing the correction. The training set contains complexes of sulfane, thiol, and thioether; to make the validation independent, we have chosen a thioketone and a sulfur heterocycle as our model systems. More specifically, we have built complexes of thioacetone and thiophene with the sulfur forming hydrogen bonds with water, methanol, pyrrole, ammonia, methylamine, sulfane, and methanethiol. As a crosscheck, we have also added two more complexes where these sulfur compounds interact with hydrogen without forming a hydrogen bond, namely complexes with molecular hydrogen and methane. Analogously to our other data sets, the complexes were optimized at the MP2/cc-pVTZ level with counterpoise correction.42 We have designed the complexes to feature a desired interaction motif. The geometries are not guaranteed to be global minima. The interaction energies in this set were calculated using exactly the same protocol as in the training sets, i.e. a composite CCSD(T)/CBS scheme43,44 with counterpoise correction.42 The composite CCSD(T)/CBS result is constructed from three contributions: Hartree−Fock energy (EHF), MP2 correlation energy (EMP2), and the correction term ΔCCSD(T) = ECCSD(T) − EMP2. The Hartree−Fock energy is calculated in the aug-cc-pVQZ basis. The MP2 correlation energy is extrapolated45 to the CBS limit from aug-cc-pVTZ and augcc-pVQZ basis sets. Finally, the ΔCCSD(T) correction is calculated in the aug-cc-pVDZ basis. The interaction energies are provided below in Table 6, and the geometries are available in the Supporting Information. The same methodology was applied to the angular scans performed in the water dimer. To test the effect of the proposed correction on thermochemistry, we use the recently published W4-17 data set of very accurate atomization energies.46 From this data set, we have selected only closed-shell organic molecules (69 items), and we use the nonrelativistic reference energies as the benchmark for DFTB. The errors against the benchmark interaction energies are evaluated as root-mean-square error (RMSE), which we consider to be the most robust error measure.33 To make the errors comparable between data sets featuring interactions of different strength, we also evaluate the relative RMSE (the coefficient of variance of the RMSE), which is calculated as the RMSE divided by the average magnitude of the interaction energy in the set and expressed in percent. 2.4. Model Systems for Proton Tranfer. To test the effects of the H5 correction on proton transfer reaction barriers, we have built four model systems covering both cations and anions: H3O+···H2O, NH4+···NH3, OH−···H2O, and HCOO−··· HCCOH complexes. Because unconstrained optimization of such dimers in gas phase would lead to a symmetric structure with the proton being shared, we use an artificial geometry where the distance between the heavy atoms is set to 2.5 Å (2.7 Å in the case of ammonia). This is a reasonable estimate of the sum of the lengths of a covalent bond and a ionic hydrogen bond in a system where the proton is localized on one of the molecules. The proton is then moved along the axis of the complex, while the remaining coordinates are optimized at the MP2/cc-pVTZ level. This approach is similar to the one used previously for testing of DFTB3.20 The benchmark reaction

dissociation curves of diverse noncovalent complexes calculated using a composite CCSD(T)/CBS scheme.30 We used an earlier parametrization of the dispersion correction,18 which had been parametrized using the dispersion-bound and mixedcharacter systems from the data set. The H5 correction was parametrized on the subset of 23 hydrogen-bonded systems from S66x8. The H5 correction for sulfur was parametrized on a complementary data set of sulfur compounds constructed in exactly the same way as S66x8.31 The reference interaction energies in both of these data sets are constructed from MP2 energy extrapolated to CBS from aug-cc-pVTZ and aug-ccpVQZ basis sets, and the CCSD(T) correction is calculated in the aug-cc-pVDZ basis. This setup has been shown to yield an average error of 2% with respect to the CBS limit.32,33 The method is tested on multiple established data sets of benchmark interaction energies which cover various aspects of the studied problem. The HB104 data set12 is the primary validation set here as it covers a large number of diverse hydrogen bonds. We employ the recently updated benchmark calculated using a CCSD(T)/CBS setup very similar to the one used in the S66x8 training set.18 To test the performance of the methods in the hydrogen bonds of charged organic species, we use the data set of 15 ionic hydrogen bonds (labeled as IHB15).15 The original data set covers dissociation curves calculated at the same level as S66x8; here we use only the equilibrium geometries (15 complexes). To sample also the angular degrees of freedom of hydrogen bonds, we use the set of hydrogen bonded systems from the S66a8 data set;34 it is 184 nonequilibrium geometries with CCSD(T)/CBS interaction energies evaluated analogously to the S88x8 data set. In addition to pairwise interactions, we also test their cooperativity by evaluating three-body energies in trimers. For this, we use the 3B69 data set that we have recently developed.35 To keep the focus on hydrogen bonds, we have selected only the subset of the most polar systems (labeled as ”low dispersion” in the original work). The three-body energy ΔE3 in a trimer ABC is defined as ΔE3(ABC) = E(ABC) − E(A) − E(B) − E(C) − ΔE2(AB) − ΔE2(BC) − ΔE2(CA)

(7)

where the ΔE2 terms are interaction energies in the respective dimers. Interaction energies in larger systems are tested in three more data sets. The L7 36 and S12L 37 sets comprise noncovalent complexes with up to about 120 atoms. From the latter set, we use only the neutral systems. The interaction energies for the S12L data set were originally derived from experimental binding free energies in solution.37 Subsequently, they were recalculated using multiple QM methods.38,39 All of these references are less accurate than the CCSD(T)/CBS benchmark used for smaller systems. Moreover, their error is hard to estimate because there is no reliable benchmark to compare to. We have tested multiple sets of reference values, and we discuss the results obtained with them in more detail. The last data set consists of water clusters up to dodecamers. Analogously to prior tests of DFTB, we have taken only the 14 neutral systems from ref 40, and we label the set WATER14. Here, the calculated quantity is cluster decomposition energy, which is the difference between the energy of the cluster and of a corresponding number of isolated water molecules. In the original work, the benchmark results were obtained using CCSD(T) in the smaller clusters and with MP2 in the large D


Article


artifacts in other systems where such a pair of atoms may not form a hydrogen bond at all (for example, this occurs in the S66 set in the stacked uracil dimer). Our past experience with the H4 correction suggests that the solution is to make the correction more short-ranged. At the cost of a small loss of accuracy in stretched H-bonds, the method becomes much more robust in more complex systems and yields better geometries. For this reason, we have tried to reduce the value of sw, and we have found that the value of sw = 0.25 leads to the most balanced results. Finally, the remaining parameters, the global scaling factor sr and the element parameters kXH, were reoptimized on the same data sets. The value of sr was optimized to 0.714. This approach represents an important reduction of the number of parameters as compared to the H4 correction, which uses pairwise parameters. The effect of the correction on the potential energy surface of not only hydrogen bonds but also of other complexes can be assessed by analyzing the dissociation curves in the S66x8 data set. The plots of the dissociation curves obtained with DFTB3D3γ, DFTB3-D3H5, and CCSD(T)/CBS benchmark are provided in the Supporting Information. The H5 correction systematically improves H-bonds, and the additional H−H repulsion used in the dispersion corrrection brings a better description of hydrocarbon dimers. Only in a few complexes featuring interaction of many nonpolar hydrogen atoms (as in aliphatic hydrocarbons) with O or N electron donor, the H5 correction leads to somewhat overestimated interaction energies at short distances when compared to the benchmark. This behavior is similar to the DFTB3-D3γ approach because both these corrections are not specific to H-bonds but act also in any other close pairs of these atoms. Nevertheless, these correction perform better than a fully empirical approach where a correction of constant size (independent of electronic structure) is added because they only amplify the electrostatic interaction already present in DFTB (which is weaker for nonpolar hydrogens). Overall, the D3H5 corrections bring a very balanced description of the whole data set. Only one exception to the transferability of the parameters has been found in the testing − the hydrogen bonds of azide nitrogen are overestimated already in uncorrected DFTB3 and the H5 only exacerbates the problem. This is a manifestation of the well-known difficulties in describing all the different chemical states of nitrogen with a single set of DFTB parameters. To avoid this issue (among the systems studied here, it occurs in the 3B69 data set), we have defined a separate atom type for the terminal azide nitrogen (easily recognized in the code as monovalent) that has the same DFTB parameters but to which the H5 correction is not applied. Fitting the D3H5 correction for sulfur has led to a substantial decrease of the overall error, but a closer inspection has revealed that in some systems, the intermolecular potential does not have a correct shape around the equilibrium because it lacks part of the repulsion. Despite that, the optimized scaling factor of 0.280 is larger than for oxygen or nitrogen because in other systems, the interaction strength is severely underestimated. Using the correction with this parametrization would result in very short intermolecular contacts and overestimated interaction energies in the cases where the repulsion is underestimated. It is known that DFTB has problems with the description of noncovalent interactions of sulfur that originate from both the electronic part and the empirical repulsive potential in DFTB.61 While this issue has improved with the 3OB parameters, some errors are still present. In order

profiles in these model systems were calculated at the CCSD(T)/CBS level using the setup described above. 2.5. Computational Details. The parametrization and testing of the D3H5 correction were carried out using an experimental implementation of DFTB in the Cuby framework.47,48 The remaining DFTB calculations were performed with the DFTB+ program,49,50 using Cuby to add dispersion and H4 corrections. The MP2 and CCSD(T) calculations in the validation set of sulfur complexes, as well as all DFT calculations reported in this paper, were performed in Turbomole 6.6.51 The newly developed method was also compared to semiempirical PM652 and PM753 calculations performed in MOPAC54 and OM355 calulations carried out using the MNDO99 program. For PM6, we use the D3H4 corrections for dispersion and hydrogen bonding,15 and OM3 was coupled with the D3 dispersion.56 The GFN-xTB57 calculations were performed using the software provided by the authors of the method. The HF-3c method58 was used as implemented in the Orca package.59

3. PARAMETRIZATION OF THE CORRECTION We use the D3 dispersion that we have already parametrized for DFTB3 without any modification. The method is based on Grimme’s D3 dispersion for DFT-D, from which we have adopted the functional form as well as the atomic parameters (C6 coefficient and atomic radii). While also other parametrizations of D3 for DFTB exist,11 we have extended it with an H−H repulsion term needed for the correction of too short intermolecular contacts in some hydrocarbons.15 The dispersion correction was parametrized on a subset of the S66x8 data set with excluded H-bonds. The parameters for DFTB3 with the 3OB parameter set have been reported in our recent paper,18 and a complete implementation of our version of the correction is available in the Cuby framework.48 The H5 correction was parametrized in several steps. First, a preliminary optimization with all three parameters freely adjustable for each element was carried out. For H-bonds of nitrogen and oxygen, the dissociation curves of H-bonded complexes from the S66x8 data set were used. For sulfur, a complementary data set of dissociation curves containing multiple H-bonded complexes was utilized. We have minimized the root-mean-square error to the CCSD(T)/CBS benchmark in each data set by means of gradient optimization based on numerically calculated derivatives of the RMSE with respect to the parameters. After the first round of optimization, we examined the dependence of the geometric parameters (position r0 and width w of the Gaussian) on the sum of the van der Waals radii60 of the interacting atoms rvdW(X) and rvdW(H). It was found that these two parameters can be well approximated as r0 = sr (rvdW(X ) + rvdW(H ))

(8)

w = sw(rvdW(X ) + rvdW(H ))

(9)

using global scaling factors sr = 0.75 and sw = 0.5. These relations hold not only for O and N but also for S, and it can be expected they are transferable to other elements because they reflect the general trends in the length of hydrogen bonds. However, the optimal width of the Gaussian, determined by the parameter sw, makes the correction act over too long a range. While such a parametrization describes well also the hydrogen bonds at longer than equilibrium distances, it would lead to E


Article


is compared to the recently reparametrized D3H4 correction18 and to Grimme’s parametrization of D3 for standard DFTB3 with γ-damping (ζ = 4) labeled here as DFTB3-D3γ.11 These methods have been chosen among multiple possible variants of H-bond corrections for DFTB because the latter represents the default treatment of H-bonds in DFTB3 and the former yields the best results at the cost of being fully empirical. The results, expressed as RMSE in each of the test sets, are summarized in Table 3; the corresponding relative errors are listed in Table 4

to balance the correction better, we have reoptimized the parameter kSH so that it yields the lowest error after geometry optimization, which leads to a lower value of the parameter. The final values of the parameters are listed in Table 1. To Table 1. Values of the Parameters − Global Parameters sr and sw and the Scaling Factors kXH Applied to the Electrostatics Interaction of Element X and Hydrogen (All Dimensionless) parameter

value

sr sw kOH kNH kSH

0.714 0.25 0.06 0.18 0.21

Table 3. Errors of the Studied Methods Including the Newly Developed DFTB3-D3H5 in Multiple Data Sets of Benchmark Interaction Energies (RMSE, in kcal/mol) data set S66x8 S66

provide data for the validation of future implementations of the method, we have provided a detailed list of interaction energies both without and with the corrections in the S66 data set as Supporting Information.

HB104 IHB15 3B69/low disp WATER14 L7 S12L/neutral sulfur

4. RESULTS AND DISCUSSION 4.1. Effect of the Corrections on Thermochemistry. Before we pass to noncovalent interactions, it is necessary to verify that the corrections do not break the description of covalent bonds. The DFTB3 method was parametrized on atomization energies, and we use the same quantity for our tests. For that, we use the recent W4-17 data set of atomization energies obtained using state-of-the-art correlated QM methods.46 The errors of DFTB3 with different corrections are listed in Table 2 (the DFTB3-D3 and DFTB3-D3γ listed there use the same dispersion correction as DFTB3-D3H5).

a

RMSE

DFTB3 DFTB3-D3 DFTB3-γ DFTB3-D3γ DFTB3-D3H5 DFTB3-H5

12.19 12.29 12.09 12.15 12.01 11.93

DFTB3D3H4

DFTB3D3H5

training set training set, equilibrium diverse ON Hbonds ionic H-bonds three-body energies

1.10 1.23

0.65 0.66

0.75 0.58

2.19

1.03

0.92

5.14 0.29

4.95/0.93a 0.37

4.41 0.25

H-bond networks large complexes large complexes validation set

6.94 2.04 3.43 0.67

22.67 2.76 3.29 0.98

8.79 2.73 3.05 0.45

Without and with specific correction for ionic H-bonds.

Table 4. Errors of the Studied Methods in Benchmark Data Sets Expressed as the Relative RMSE in Percent data set S66x8 S66

Table 2. Root Mean Square Error of Atomization Energies (in kcal/mol) in Closed-Shell Organic Molecules from the W4-17 Data Set method

DFTB3D3γ

description

HB104 IHB15 3B69/low disp WATER14 L7 S12L/ neutral sulfur averageb

DFTB3D3γ

DFTB3-D3H4

DFTB3D3H5

training set training set, equilibrium diverse ON Hbonds ionic H-bonds three-body energies H-bond networks large complexes large complexes

27.3% 22.3%

16.2% 12.1%

18.6% 10.5%

34.1%

16.0%

14.2%

26.3% 34.3%

25.3%/4.7%a 43.4%

22.6% 29.4%

7.5% 11.2% 10.3%

24.5% 15.2% 9.9%

9.5% 15.0% 9.2%

validation set

31.5% 22.2%

46.1% 25.8%/22.8%a

21.1% 17.3%

description

a

Without and with specific correction for ionic H-bonds. bExcluding the training set S66x8 and S66.

All the errors are very similar, and the effect of the corrections on the thermochemistry can be neglected. The small diffences observed can be explained in the following way: Adding a dispersion correction (which was not included in the parameterziation of DFTB3) leads to a minor increase of the error. On the other hand, the γ-damping was used in the parametrization, and it thus improves the results compared to plain DFTB3. Surprisingly, the H5 correction improves the results even more, and the DFTB3-D3H5 method is thus slightly more accurate than the original DFTB3 with γdamping. 4.2. Interaction Energies. We have tested the performance of the D3H5 correction first on interaction energies in fixed geometries. Besides the training set (S66x8) and the corresponding set of equilibrium geometries (S66), we report the results for multiple independent validation sets that highlight specific areas of the problem. The D3H5 correction

and plotted in Figure 2. The relative errors, which are comparable across the different data sets, represent the best measure for an overall assessment of the tested methods. It is clear from the results that D3H5 not only yields the best results when averaged over the data sets (the data sets used in the parametrization, S66x8 and S66, have been excluded from the average) but also performs the most consistently. The most important check is the test of the transferability of the parameters to different hydrogen-bond motifs not covered by the training set. For this purpose, we use the HB104 data set,12 a set of 104 diverse hydrogen bonds of O and N. As a benchmark, we use the CCSD(T)/CBS interaction energies that we have recalculated recently.18 Here, the D3H5 corrections with the RMSE of 0.92 kcal/mol outperform the F


Article


Figure 2. Errors of the studied methods in benchmark data sets expressed as the relative RMSE in percent. The dotted line represents the average over the validation sets. For DFTB3-D3H4, results are plotted without and with specific correction for ionic H-bonds.

Figure 3. Distribution of the DFTB3-D3H5 errors for OH-O, NH-O, OH-N, and NH-N interaction motifs in the HB104 data set. Center lines show the medians; box limits indicate the 25th and 75th percentiles; whiskers extend 1.5 times the interquartile range from the 25th and 75th percentiles; outliers are represented by dots.

best DFTB-based method available so far, DFTB3-D3H4 (a RMSE of 1.03 kcal/mol). Both of these methods are significantly better than the default γ-damping in DFTB3-D3 (RMSE 2.19 kcal/mol). It is important to note that D3H5 uses only element-wise parameters (in this case two) while D3H4 uses (four) pairwise parameters, yet it performs comparably well. This can be explained by the fact that unlike H4, the H5 correction takes into account the partial atomic charges of the atoms involved in the hydrogen bond. The HB104 data set is a useful tool for analyzing the difference between the hydrogen bonds of oxygen and nitrogen. The set can be divided into four groups, featuring OH-O, NHO, OH-N, and NH-N interaction motifs; the errors of the studied methods for these groups are reported in Table 5. It is

however, not the main limitation of the method, and the possible slight improvement in accuracy does not justify the additional complexity and loss of generality. The same feature suggests that the H5 correction should work better for the hydrogen bonds of ionic species. This has been tested using a data set of hydrogen bonds covering interactions of ionic amino-acid side chains with small molecules.15 DFTB3-D3γ works the worst here with a RMSE of 5.14 kcal/mol. In the H4 correction, a special scaling term is applied to each ionic group covered by this data set (carboxylate, ammonium, guanidinium and imidazolium), which is parametrized on this set. With this scaling, it is possible to achieve excellent results in these systems (RMSE of 0.93 kcal/mol) or similar binary complexes. The transferability of this correction is, however, rather poor, and it does not work well e.g. in the cases where a single ionic group forms multiple H-bonds. Without this scaling, the H4 correction yields a RMSE of 4.95 kcal/mol. The H5 correction yields a RMSE of 4.41 kcal/mol without any specific parametrization, underestimating the strength of the interaction in 14 cases out of 15. The interaction energies are, however, substantially larger in the ionic systems, so that in relative terms, the error is not much larger than in neutral systems. The relative RMSE in the IHB15 data set is about 23%, while in the HB104 set it is 14%. In the ionic H-bonds, the polarization term becomes more important, and it cannot be fully recovered by a correction to the firstorder electrostatics. Here, a polarization-based correction based on the chemical potential equilibration17 (CPE) has worked very well although it fails in other cases.18 The inclusion of the correction in the SCC cycle means that it affects the electron density in the interacting molecules, which should improve the description of many-body effects. We have tested this hypothesis on the data set of three-body energies 3B69,35 from which we took the subset of 27 most polar systems (the Low-dispersion group as defined in ref 35). In this test, the three-body term in the D3 dispersion23 was switched on, but its effect in these polar complexes is negligible. Here, both DFTB3-D3γ and D3H5 approaches work comparably well with a RMSE of 0.29 and 0.25 kcal/mol, respectively, while the pairwise D3H4 approach exposes the underestimation of the many-body effects in plain DFTB3 with an error of 0.37 kcal/mol. The errors below 0.3 kcal/mol are

Table 5. Errors in the HB104 Data Set Separated for All the Combinations of the Donor and Acceptor Elements (RMSE, in kcal/mol) interaction motif

DFTB3-D3γ

DFTB3-D3H4

DFTB3-D3H5

OH-O NH-O OH-N NH-N

0.75 1.25 3.84 2.70

0.71 0.95 1.52 1.03

0.58 1.17 0.80 1.10

obvious that there is only a little difference, caused by the element at the side of the electron acceptor, because it is only the hydrogen atom that directly takes part in the H-bond. A significant difference can be, however, observed at the side of the electron donor. The errors in the hydrogen bonds of oxygen are substantially smaller than in those of nitrogen. This is most pronounced in the DFTB3-D3γ method, where the correction does not distinguish between the elements. In both D3H4 and D3H5 approaches, the H-bonds of oxygen are described only slightly better, but the error in the H-bonds of nitrogen drops to less than a half. DFTB3-D3H5 provides the most balanced description of all the combinations of donor and acceptor element despite not using pairwise parameters. The distributions of the errors for OH-O, NH-O, OH-N, and NH-N hydrogen bonds in the HB104 data set are plotted in Figure 3. The hydrogen bonds of nitrogen (especially the OH-N bonds, which tend to be stronger) are still the most difficult cases to describe, which is in line with the general difficulties in describing nitrogen in DFTB. Some additional improvement could be achieved if atom types were introduced, distinguishing between hydrogens bound to different elements. This is, G


Article

Journal of Chemical Theory and Computation very good in the realm of SQM methods − for comparison, GGA DFT functionals yield an error of about 0.13 kcal/mol in the same set.35 The many-body effects studied separately in trimers become important when we pass from small models to larger clusters or condensed phase, which contain complex networks of hydrogen bonds. We have evaluated these effects in a set of 14 water clusters consisting of up to 20 water molecules,40 for which accurate CCSD(T)/CBS reference became available recently.41 This set was already used in the testing of DFTB methods under the name WATER14.17,18 Again, the worst results are obtained with the pairwise H4 correction, which does not describe the cooperativity of the H-bonds and thus yields a large error of 22.7 kcal/mol, corresponding to the relative error of 25%. The H5 correction brings a substantial improvement with a RMSE of 8.8 kcal/mol, where the relative error of 9.5% is better than the one in the more diverse set of small model complexes. The DFTB3-D3γ method performs slightly better here with a RMSE of 7.0 kcal/mol (the relative error of 7.5%), which can be attributed to the fact that the exponent ζ in the γdamping is fitted to a single system, the water dimer. These results highlight the importance of the inclusion of many-body effects in the correction. While both D3H4 and D3H5 describe small complexes including the water dimer with similar accuracy, D3H4 fails in the water clusters because it is pairwise additive, while D3H5 does cover the many-body effects. Other tests of the transferability of the method to larger systems are the L7 and S12L data sets. From the latter, we have used only the nine neutral complexes. These data sets are, however, dominated mainly by dispersion (mostly π−π interactions): the L7 data set contains only one system in which H-bonds are mixed with dispersion, and the subset of the S12L used here features only three H-bonded systems. In the L7, all the tested methods perform comparably; DFTB3-D3γ has the smallest error because it describes the π−π interactions best, whereas the D3 correction in D3H4 and D3H5 methods overestimates it slightly. On the other hand, the most serious error of DFTB3-D3γ is the overestimated dispersion in the single aliphatic hydrocarbon dimer in the set. These errors (in the range of 2−3 kcal/mol) translate to relative errors comparable to those obtained in small model systems, which indicates that the dispersion correction is well transferable to larger systems. The situation is more complex in the S12L data set, where the benchmark interaction energies are less reliable. Previously, we have obtained very large errors (close to 10 kcal/mol) for all the DFTB3-based methods, attributing this to the composition of the data set, which contains mostly strong π−π complexes. Nevertheless, it was suspicious that only slightly smaller systems of similar nature present in the L7 set were described much better. These results were obtained using Grimme’s original benchmark interaction energies derived from experimental dissociation energies in solution. When we pass to the best gas-phase computational benchmark available, the DFTSAPT interaction energies calculated by Hesselmann and Korona, the results change significantly, and the errors drop to a range comparable to the L7 data set. In this work, we have thus used only the DFT-SAPT interaction energies as a benchmark. In DFT-D3 calculations of these large systems, it is beneficial to include the three-body dispersion energy in the D3 correction23 which improves the transferability of the correction parametrized in small systems to large ones. We

have tested the same approach in DFTB3-D3H5, but the results in the larger data sets become slightly worse (RMSE 3.7 kcal/mol in both L7 and S12L). This is because the dispersion was parametrized so that interactions at longer distances (which are more important in large systems) are slightly underestimated in order to conserve physically correct global scaling factor in the dispersion correction.18 Reducing the dispersion contribution further by means of the three-body term thus leads to worse results. Not using the three-body dispersion is also beneficial for computational efficiency because this calculation, and especially the gradient of the three-body dispersion energy, becomes nontrivial in larger systems. Finally, we have tested the description of the noncovalent interactions of sulfur in the new validation set constructed for this purpose. It should be mentioned that it contains sulfur atoms of different atom types (thioketone and aromatic heterocycle) from those used in the training set (thioether, thiol, sulfane). This validation set should hence provide information on not only the accuracy but also transferability of the correction. In DFTB3-D3H4, only the dispersion correction is applied to sulfur, but no correction is used for its hydrogen bonds. It yields a RMSE of 0.98 kcal/mol, which is a large error in a data set where the average interaction energy amounts to −2.1 kcal/mol. Both the γ-damping and D3H5 apply H-bonding correction to sulfur, and the errors drop to 0.67 and 0.45 kcal/mol, respectively. Even the performance of D3H5 is not completely satisfactory as it translates to the relative error of 21%, but it cannot be improved further by this correction as it is caused by errors inherent to the DFTB parameters. The individual results along with the benchmark interaction energies are listed in Table 6. The errors range from positive to negative values with an average of 0.24 kcal/mol. Such behavior is expected as the correction was intentionally made weaker in order to avoid extremely close contacts in some of the systems in the training set. The errors are very similar in the complexes of thioacetone and thiophene although the Table 6. DFTB3-D3H5 Interaction Energies in a Validation Set of Sulfur-Containing Complexes along with the CCSD(T)/CBS Reference and Errors with Respect to It (All Values in kcal/mol)

H

complex

CCSD(T)/ CBS

DFTB3-D3H5

error

thioacetone···H2 thioacetone···ammonia thioacetone···methane thioacetone···methanethiol thioacetone···methanol thioacetone···methylamine thioacetone···pyrrole thioacetone···sulfane thioacetone···water thiophene···H2 thiophene···ammonia thiophene···methane thiophene···methanethiol thiophene···methanol thiophene···methylamine thiophene···pyrrole thiophene···sulfane thiophene···water

−0.46 −3.35 −0.63 −2.55 −5.43 −3.94 −5.58 −2.61 −5.02 −0.28 −0.65 −0.45 −1.1 −1.48 −0.88 −1.52 −1.04 −1.26

−0.35 −2.58 −0.57 −2.85 −5.22 −3.14 −4.55 −3.07 −5.06 −0.32 −0.54 −0.48 −0.92 −0.94 −0.78 −0.98 −0.81 −0.75

0.11 0.77 0.05 −0.30 0.21 0.80 1.03 −0.46 −0.03 −0.04 0.12 −0.03 0.17 0.53 0.10 0.54 0.23 0.51


Article


geometry is a result of their interplay. Our model was the water dimer, in which we performed scans of two angles shown in Figure 4. In both cases, all the DFTB-based methods yielded similar profiles of interaction energy, although the curve was shifted because each of the methods yielded a slightly different interaction energy in the equilibrium. In the H−O···H angle α, the DFTB potential is more flat than the CCSD(T)/CBS benchmark. Since the H4 correction includes an angular term defined in terms of this angle, it describes better the rise of the energy at values of α far from the equilibrium. On the other hand, the angular function used there is rather flat around the equilibrium, and this translates to a potential in this region that is flatter than the one from D3γ or D3H5 corrections. The predicted geometry is, however, rather good; at the benchmark level, the equilibrium value of α is 173°, while all the DFTB methods yield the minimum around 180°. The angle β represents the rotation around the oxygen participating in the H-bond, and the angular profile thus describes the role of the lone pair in the interaction. Here, none of the corrections make use of this angle, and the electrostatics in DFTB are based on monopole approximation, so that the lone pair cannot be accounted for. The only means of description for this angular dependence are the dipole moment of the whole molecule (which would yield the most favorable interaction with the hydrogen at β = 180°) and the terms coming from the zeroth-order Hamiltonian, which is resolved at the level of orbitals. Since DFTB describes the geometry rather well, it is obvious that the latter plays an important role here, and the form of the correction does not play any significant role. None of the DFTB methods, however, correctly describes the flattening of the potential at the lowest values of β. In summary, however, DFTB describes the angular dependence of the H-bond rather well, despite the severe approximations that are involved. To test the description of angular intermolecular degrees of freedom in a more diverse set of systems, we use the hydrogenbonded systems from the S66a8 data set.34 This data set comprises nonequilibrium geometries prepared by rotating the monomers around the atom involved in the interaction by 30 degrees. For each system from the S66 data set, eight nonequilibrium geometries were prepared (each monomer is rotated in both directions along the principal axes perpendicular to the intermolecular vector). To assess the description of the potential energy surface out of equilibrium, the error calculated on these displaced geometries is compared to the error obtained in the equilibrium. The results are listed in Table 9. Because the H-bonds become weaker in the nonequilibrium geometries, we will discuss only the relative errors (calculated as RMSE divided by the average interaction energy in the set). Both DFTB3-D3 and DFTB3-D3γ yield rather large errors that are almost the same in equilibrium and displaced geometries. DFTB3-D3H4 and DFTB3-D3H5 perform better, but the relative error grows from 7% to 12−13% when we pass to the nonequilibrium geometries. This can be expected because we have already shown on the water dimer that the angular potential in these methods is more flat that it should be (what is also confirmed by the fact that the strength of the interactions in the S66a8 set is overestimated). Nevertheless, the error is still acceptably small. Comparing the H4 and H5 corrections, the latter performs remarkably well when we take into account that H4 contains an angular term while H5 does not.

former is substantially more polar and forms stronger complexes. 4.3. Geometries. Other important characteristics of the method are related to the geometries obtained with it. First, the method should reproduce the geometry itself; second, the interaction energy in the minimum obtained with the tested method should be close to the benchmark. We have tested these properties in the S66 data set, analyzing separately the three groups of systems of which it consists, hydrogen bonds, dispersion-dominated and mixed-character complexes. Each system in the set was optimized using the tested methods, and we evaluated the root-mean-square difference (RMSD) of the optimized geometry to the benchmark one (optimized using MP2/cc-pVTZ, with the intermolecular distance corrected at the CCSD(T)/CBS level). The geometry differences in the individual groups, as well as in the whole set, are listed in Table 7. Here, D3H4 and D3H5 corrections Table 7. Geometry Difference from the Benchmark (RMSD, in Å) Averaged over the Three Groups of Complexes in the S66 Data Set H-bonds dispersion other all

DFTB3-D3γ

DFTB3-D3H4

DFTB3-D3H5

0.14 0.15 0.17 0.15

0.19 0.18 0.18 0.18

0.19 0.18 0.21 0.19

perform similarly, slightly worse than the D3γ approach. This is because in the latter, the correction is weaker, and thus it introduces fewer artifacts in the geometries. If only the intermolecular distance is evaluated, D3H5 provides the best description of H-bonds with an average (signed) error of 0.00 Å and D3γ with an error of −0.01 Å; D3H4 overestimates the length of the H-bonds by 0.06 Å on average. The error of the geometries is acceptably small in all the tested methods. For practical use, it is even more important to evaluate also the error of interaction energies in optimized geometries (Table 8). Here, D3γ correction yields the largest Table 8. RMSE of Interaction Energies (in kcal/mol) in the S66 Data Set after Geometry Optimization (by Group and Total) H-bonds dispersion other all

DFTB3-D3γ

DFTB3-D3H4

DFTB3-D3H5

1.60 0.80 0.60 1.18

0.92 0.77 0.89 0.86

0.94 0.90 0.63 0.84

error, while D3H4 and D3H5 work better, especially for hydrogen bonds. These results also show one weakness of the H5 correction, slight deterioration of the results in dispersiondominated complexes. This is because the correction is applied to all pairs of hydrogen with oxygen or nitrogen, so that it affects also interactions of CH groups with these heteroatoms. The increase of the error from 0.8 to 0.9 kcal/mol is, however, not important enough to justify the introduction of the atom types that would make it possible to correct only true H-bonds. Finally, we should also examine the angular dependence of the interaction energy in a hydrogen bond. This is an important characteristics that may affect the description of more complex systems where other interactions are present and the final I


Article


Figure 4. Angular scans in the water dimer (interaction energy corresponding to geometries where only the selected angle is varied).

Table 9. Errors of the DFTB3-Based Method in Angularly Displaced and Equilibrium Geometries of Hydrogen Bonds (from S66a8 and S66 Data Sets)a method DFTB3-D3 DFTB3-D3γ DFTB3-D3H4 DFTB3-D3H5

displaced 2.08 1.33 0.81 0.91

equilibrium

(30.4%) (19.4%) (11.8%) (13.3%)

2.76 1.72 0.59 0.64

(31.1%) (19.3%) (6.6%) (7.2%)

a

RMSE in kcal/mol and RMSE relative to the average interaction energy in the respective set.

4.4. Proton Transfer. Semiempirical QM methods often have difficulties with describing proton transfer (PT) and yielding correct barrier heights of PT reactions. In the development of DFTB3, special attention was paid to this issue, and proton transfer barriers were used in the fitting of the repulsive potential used in the 3OB parameter set.20 Despite that, the barrier heights are still underestimated. Adding a correction for hydrogen bonding affects the PT energetics, and it should be tested how it changes the results. We have prepared four model systems where the average proton transfer barrier height calculated at the CCSD(T)/CBS level is 1.05 kcal/mol. The benchmark barrier heights and the results obtained with DFTB3-D3, DFTB3-D3γ, and DFTB3D3H5 are listed in Table 10. The complete proton transfer potential in the H3O+···H2O complex is plotted in Figure 5 for illustration.

Figure 5. Proton transfer potential in a H3O+···H2O complex with O− O distance fixed at 2.5 Å.

more than the minima which leads to a barrier lower than in DFTB3 without γ-damping. 4.5. Comparison with Other Methods. So far, we have discussed the performance of the DFTB3-D3H5 method in the context of other approaches based on DFTB. Now, we will briefly compare it to other methods of similar complexity as well as to more accurate but computationally much more demanding methods. At the semiempirical level, we compare DFTB3-D3H5 to the following methods: First, PM6-D3H4 represents a state-of-the-art MNDO-based SQM method with advanced corrections for dispersion and hydrogen bonding.15,52 Second, PM7 is a successor of PM6, and it already includes the corrections for noncovalent interactions.53 Despite that, it describes them with lower accuracy than PM6-D3H4.62 We also include the OM3 method,55 the most recent method from a series of SQM methods that use additional orthogonalization corrections.63 It is coupled with a corresponding parametrization of the D3 dispersion correction.56 Third, we test the new self-consistent tight-binding GFN-xTB method (which is an empirical tight-binding approach, constructed differently than DFTB).57 Subsequently, to demonstrate the strengths and limitations of the SQM methods, we compare them to full SCF calculations corrected for noncovalent interactions. The HF-3c method58 represents the most efficient approach applicable to noncovalent interactions; it combines a Hartree−Fock calculation in minimal basis with three additional corrections: an empirical dispersion and corrections for basis set superposition error and bond lengths.64 Finally, an accurate DFT-D method, namely B-LYP/def2-QZVP calulation combined with

Table 10. Proton Transfer Barriers in Model Complexes (in kcal/mol) Calculated Using CCSD(T)/CBS and Three Variants of DFTB3 system H3O+···H2O NH4+···NH3 OH−···H2O HCOO−···HCCOH average error

CCSD(T)/ CBS

DFTB3D3

DFTB3D3γ

DFTB3D3H5

0.72 2.03 0.66 0.80

0.35 0.23 0.33 0.23 0.77

0.15 0.27 0.25 0.19 0.84

0.63 1.36 0.88 0.46 0.33

The H5 correction improves the barrier heights significantly, and the error in comparison with the CCSD(T)/CBS benchmark is only 0.33 kcal/mol (30% of the average barrier height). This can be explained by the fact that the H5 correction stabilizes the minima where there is a well-defined H-bond more than the transition state. On the other hand, the γ-damping is less specific, and it stabilizes the transition state J


Article

Journal of Chemical Theory and Computation Table 11. Errors of DFTB3-D3H5 Compared to Other Semiempirical Methods and DFT-D3 in Multiple Data Sets of Benchmark Interaction Energies (RMSE, in kcal/mol)

a

data set

description

DFTB3-D3H5

PM6-D3H4

PM7

OM3-D3

GFN-xTB

HF-3c

DFT-D3

S66x8 S66 HB104 IHB15 3B69/low disp WATER14 L7 S12L/neutral sulfur

training set training set, equibrium diverse ON H-bonds ionic H-bonds three-body energies H-bond networks large complexes large complexes validation set

0.75 0.58 0.92 4.41 0.25 8.79 2.73 3.05 0.45

0.66 0.67 0.96 2.57/0.99a 0.51 7.23 3.90 3.27 0.80

0.95 0.98 1.64 1.85 0.53 5.29 5.66 13.15 1.54

1.17 1.27 1.39 2.15 0.43 1.86 3.21

1.05 1.27 1.35 2.57 0.29 13.92 4.59 5.69 0.20

0.53 0.46 0.67 4.38 0.29 10.68 2.07 1.37 0.94

0.22 0.24 0.24 0.63 0.12 4.15 1.63 4.76 0.36

Without and with specific correction for ionic H-bonds.

Table 12. Relative Errors (RMSE in Percent of Average Magnitude of Interaction Energy in Each Data Set) of DFTB3-D3H5 Compared to Other Semiempirical Methods and DFT-D3

a

data set

description

DFTB3-D3H5

PM6-D3H4

PM7

GFN-xTB

HF-3c

DFT-D3

S66x8 S66 HB104 IHB15 3B69/low disp WATER14 L7 S12L/neutral sulfur averageb

training set training set, equibrium diverse ON H-bonds ionic H-bonds three-body energies H-bond networks large complexes large complexes validation set

18.6% 10.5% 14.2% 22.6% 29.4% 9.5% 15.0% 9.2% 21.1% 17.3%

16.4% 12.2% 15.0% 13.2%/5.0%a 60.4% 7.8% 21.4% 9.8% 37.5% 23.6%/22.4%a

23.6% 17.8% 25.5% 9.5% 62.5% 5.7% 31.1% 39.5% 72.5% 35.2%

26.2% 23.2% 21.0% 13.2% 33.9% 15.0% 25.2% 17.1% 9.6% 19.3%

13.1% 8.4% 10.5% 22.4% 34.0% 11.5% 11.4% 4.1% 44.5% 19.8%

5.4% 4.4% 3.8% 3.2% 13.8% 4.5% 9.0% 14.3% 16.8% 9.3%

Without and with specific correction for ionic H-bonds. bExcluding the S66x8 and S66 data sets.

the D3(BJ) dispersion,65 represents one of the most accurate mean-field approaches still applicable to rather large molecular systems. Significantly higher accuracy can be achieved only with correlated quantum-mechanical methods which are substantially more expensive. The results are summarized in Table 11, and corresponding relative errors are listed in Table 12. PM6-D3H4 performs similarly to DFTB-D3H5; the two notable differences are the poor description of the three-body interactions and worse results for the sulfur data set. Here, the main advantage of DFTB over PM6 is not visible − it is the insufficient robustness of PM6 in systems with other elements for which the parametrization is not as good as for organic molecules.4 PM7 works only slightly worse for neutral hydrogen bonds (which is caused by a simpler form of the correction built into the method) but yields better results for charged ones. This method has, however, a serious problem when we pass to larger systems (represented by the L7 and S12L data sets) where the interaction energies are significantly overestimated. This is caused by overestimated dispersion at longer range. We have investigated this issue previously.62 Despite not having any specific corrections for hydrogen bonds, the OM3-D3 method describes them rather well. The extremely small error in the data set of water clusters is probably only an error compensation not occurring in other systems, which is indicated both by the larger errors in more diverse set of H-bonds and by large error in the three-body energies. It would be a promising method for applications, but its parametrization is limited only to organic elements. The GFN-xTB method exhibits a very balanced performance. Overall, it yields errors only slightly worse than DFTB3-D3H5, and it describes most of the data sets reasonably well. The

advantage of this method is a wider parametrization compared to DFTB, and the very good results in the sulfur data set indicate that it remains accurate even for other elements. The HF-3c method yields errors comparable to the best semiempirical methods and DFTB, but it is much more computationally demanding. The use of such a method could be justified if it was more robust than the more approximate methods, but it is not so. The larger error for the sulfur test sets (as well as tests on other elements not reported here) suggests that the parametrization of the corrections is not transferable well to other elements. This is due to the fact that the method is based on empirical corrections of very large errors caused by the small basis set used. As expected, DFT-D3 in large basis sets performs the best in this comparison. The largest errors are observed in the threebody energies (we have shown that DFT has difficulties describing these effects35) and in the L7 and S12L data sets. The latter can be improved further by the use of the three-body term in the D3 correction, but it could also increase the error in the smaller systems on which the correction was parametrized.23 These results show that the best corrected semiempirical methods perform surprisingly well, with errors only about twice as large as the much more expensive DFT-D3 approach. Among these semiempirical methods, the newly developed DFTB3-D3H5 approach yields the smallest error.

5. CONCLUSIONS We have developed a novel correction for hydrogen bonding in DFTB. It is combined with D3 dispersion correction, which we parametrized earlier; we call the resulting method DFTB3D3H5. Unlike our previous empirical correction H4, which was K


Article


however, used in the future as a part of a more physical model where it would allow further fine-tuning by empirical means.

calculated independently, the H5 correction modifies the native description of the electrostatics in DFTB and enters the SCC cycle. It is an extension of the γ-damping approach used in DFTB, which removes its main limitations, the problematic distance dependence and the neglect of the differences between hydrogen-bond acceptor elements. It can also be viewed as an empirical term reproducing the short-ranged electrostatic and polarization contributions lost due to the monopole approximations inherent to DFTB. It uses only a very small number of free parameters (one per each electron donor element and one global parameter), which is a substantial reduction as compared to the D3H4 approach, using pairwise parameters. It was verified that the correction does not deteriorate the description of covalent bonds. We used our previous parametrization of the D3 correction.18 The H5 correction was parametrized on the CCSD(T)/CBS dissociation curves of H-bonds from the S66x8 data set and an analogous data set of sulfur compounds. It was then validated in multiple established data sets, where it consistently performs well. The first important test was transferability to a different chemical context, that is to the H-bonds of functional groups not used in the parametrization. Here, in the data set of 104 diverse hydrogen bonds of O and N as well as in the newly developed validation set of the H-bonds of sulfur, it yields the lowest error of all the DFTB-based methods available. Second, the accuracy achieved in the small model systems transfers well to larger, more complex systems. This is because the D3H5 correction, due to its sounder physical basis, describes the many-body effects better. This is demonstrated not only in the explicit calculation of three-body energies in trimers but also by a significant improvement in the description of water clusters, where the cooperativity of Hbonds plays a crucial role. When compared to the other methods in this class, it is to other semiempirical QM methods with corrections for noncovalent interactions, it outperforms even the best of them by at least a small margin. Our tests also show that these methods had reached accuracy only about two times worse than the much more expensive DFT-D3 in the large basis set. The new method performs well not only for interaction energies but also for the prediction of the geometries of the noncovalent complexes. The overall quality of the geometries measured by their RMSD with respect to a benchmark is comparable to other similar approaches, but DFTB-D3H5 has two more important features. First, it has a negligible systematic error of hydrogen bond lengths. Second, it yields very good interaction energies on the optimized geometries, which is very important for practical use. Finally, we have shown that all the DFTB-based methods tested here describe the angular dependence of the interaction energy of a H-bond rather well, although the potential around the minimum is flatter than it should be. The H5 correction also significantly improves barriers for proton transfer along a hydrogen bond. Systematic improvement toward the CCSD(T)/CBS benchmark was observed in both anionic and cationic complexes with H-bonds involving both oxygen and nitrogen. In addition to its practical value, the H5 method also demonstrates that even a simple correction applied at the right, physically justified place can be more successful than a more complex but empirical approach. On the other hand, it is still not a perfect solution as it does not have all the properties that a true inclusion of atomic multipoles would have. It could be,

■

ASSOCIATED CONTENT

S Supporting Information *

The Supporting Information is available free of charge on the ACS Publications website at DOI: 10.1021/acs.jctc.7b00629. Table of interaction energies in the S66 data set calculated with and without the corrections and plots of the S66x8 dissociation curves (PDF) Coordinates of the complexes forming the data set for the validation of H-bonds of sulfur (ZIP)

■

AUTHOR INFORMATION

Corresponding Author

*E-mail: [email protected]. ORCID

Jan Ř ezác:̌ 0000-0001-6849-7314 Funding

This work has been supported by the Czech Science Foundation by grant no. P208/16-11321Y and is a part of the Research Project RVO 61388963 of the Institute of Organic Chemistry and Biochemistry of the Czech Academy of Sciences. Notes

The author declares no competing financial interest.

■

REFERENCES

(1) Christensen, A. S.; Kubař, T.; Cui, Q.; Elstner, M. Semiempirical Quantum Mechanical Methods for Noncovalent Interactions for Chemical and Biochemical Applications. Chem. Rev. 2016, 116, 5301− 5337. (2) Lepšík, M.; Ř ezác,̌ J.; Kolár,̌ M.; Pecina, A.; Hobza, P.; Fanfrlík, J. The Semiempirical Quantum Mechanical Scoring Function for In Silico Drug Design. ChemPlusChem 2013, 78, 921−931. (3) Pecina, A.; Meier, R.; Fanfrlík, J.; Lepšík, M.; Ř ezác,̌ J.; Hobza, P.; Baldauf, C. The SQM/COSMO Filter: Reliable Native Pose Identification Based on the Quantum-Mechanical Description of Protein-ligand Interactions and Implicit COSMO Solvation. Chem. Commun. 2016, 52, 3312−3315. (4) Pecina, A.; Haldar, S.; Fanfrlik, J.; Meier, R.; Ř ezác,̌ J.; Lepsik, M.; Hobza, P. The SQM/COSMO Scoring Function at the DFTB3-D3H4 Level: Unique Identification of Native Protein-Ligand Poses. J. Chem. Inf. Model. 2017, 57, 127. (5) Elstner, M.; Hobza, P.; Frauenheim, T.; Suhai, S.; Kaxiras, E. Hydrogen Bonding and Stacking Interactions of Nucleic Acid Base Pairs: A Density-Functional-Theory Based Treatment. J. Chem. Phys. 2001, 114, 5149. (6) Zhechkov, L.; Heine, T.; Patchkovskii, S.; Seifert, G.; Duarte, H. A. An Efficient a Posteriori Treatment for Dispersion Interaction in Density-Functional-Based Tight Binding. J. Chem. Theory Comput. 2005, 1, 841−847. (7) Martin, B.; Clark, T. Dispersion Treatment for NDDO-Based Semiempirical MO Techniques. Int. J. Quantum Chem. 2006, 106, 1208−1216. (8) McNamara, J. P.; Hillier, I. H. Semi-Empirical Molecular Orbital Methods Including Dispersion Corrections for the Accurate Prediction of the Full Range of Intermolecular Interactions in Biomolecules. Phys. Chem. Chem. Phys. 2007, 9, 2362. (9) Morgado, C. A.; McNamara, J. P.; Hillier, I. H.; Burton, N. A.; Vincent, M. A. Density Functional and Semiempirical Molecular Orbital Methods Including Dispersion Corrections for the Accurate Description of Noncovalent Interactions Involving Sulfur-Containing Molecules. J. Chem. Theory Comput. 2007, 3, 1656−1664. L


Article


Symmetry-Adapted Perturbation Theory. J. Chem. Theory Comput. 2017, 13, 1638−1646. (30) Ř ezác,̌ J.; Riley, K. E.; Hobza, P. S66: A Well-Balanced Database of Benchmark Interaction Energies Relevant to Biomolecular Structures. J. Chem. Theory Comput. 2011, 7, 2427−2438. (31) Mintz, B. J.; Parks, J. M. Benchmark Interaction Energies for Biologically Relevant Noncovalent Complexes Containing Divalent Sulfur. J. Phys. Chem. A 2012, 116, 1086−1092. (32) Ř ezác,̌ J.; Dubecký, M.; Jurečka, P.; Hobza, P. Extensions and Applications of the A24 Data Set of Accurate Interaction Energies. Phys. Chem. Chem. Phys. 2015, 17, 19268−19277. (33) Ř ezác,̌ J.; Hobza, P. Benchmark Calculations of Interaction Energies in Noncovalent Complexes and Their Applications. Chem. Rev. 2016, 116, 5038−5071. (34) Ř ezác,̌ J.; Riley, K. E.; Hobza, P. Extensions of the S66 Data Set: More Accurate Interaction Energies and Angular-Displaced Nonequilibrium Geometries. J. Chem. Theory Comput. 2011, 7, 3466−3470. (35) Ř ezác,̌ J.; Huang, Y.; Hobza, P.; Beran, G. J. O. Benchmark Calculations of Three-Body Intermolecular Interactions and the Performance of Low-Cost Electronic Structure Methods. J. Chem. Theory Comput. 2015, 11, 3065−3079. (36) Sedlák, R.; Janowski, T.; Pitoňaḱ , M.; Ř ezác,̌ J.; Pulay, P.; Hobza, P. Accuracy of Quantum Chemical Methods for Large Noncovalent Complexes. J. Chem. Theory Comput. 2013, 9, 3364−3374. (37) Grimme, S. Supramolecular Binding Thermodynamics by Dispersion-Corrected Density Functional Theory. Chem. - Eur. J. 2012, 18, 9955−9964. (38) Ambrosetti, A.; Alfè, D.; DiStasio, R. A.; Tkatchenko, A. Hard Numbers for Large Molecules: Toward Exact Energetics for Supramolecular Systems. J. Phys. Chem. Lett. 2014, 5, 849−855. (39) Heßelmann, A.; Korona, T. Intermolecular Symmetry-Adapted Perturbation Theory Study of Large Organic Complexes. J. Chem. Phys. 2014, 141, 094107. (40) Bryantsev, V. S.; Diallo, M. S.; van Duin, A. C. T.; Goddard, W. A. Evaluation of B3LYP, X3LYP, and M06-Class Density Functionals for Predicting the Binding Energies of Neutral, Protonated, and Deprotonated Water Clusters. J. Chem. Theory Comput. 2009, 5, 1016−1026. (41) Manna, D.; Kesharwani, M. K.; Sylvetsky, N.; Martin, J. M. L. Conventional and Explicitly Correlated Ab Initio Benchmark Study on Water Clusters: Revision of the BEGDB and WATER27 Data Sets. J. Chem. Theory Comput. 2017, 13, 3136. (42) Boys, S.; Bernardi, F. Calculation of Small Molecular Interactions by Differences of Separate Total Energies - Some Procedures with Reduced Errors. Mol. Phys. 1970, 19, 553−566. (43) Tsuzuki, S.; Honda, K.; Uchimaru, T.; Mikami, M.; Tanabe, K. Origin of Attraction and Directionality of the π/π Interaction: Model Chemistry Calculations of Benzene Dimer Interaction. J. Am. Chem. Soc. 2002, 124, 104−112. (44) Jurečka, P.; Hobza, P. On the Convergence of the (ΔECCSD(T)-ΔEMP2) Term for Complexes with Multiple H-Bonds. Chem. Phys. Lett. 2002, 365, 89−94. (45) Halkier, A.; Helgaker, T.; Jørgensen, P.; Klopper, W.; Koch, H.; Olsen, J.; Wilson, A. K. Basis-Set Convergence in Correlated Calculations on Ne, N2, and H2O. Chem. Phys. Lett. 1998, 286, 243−252. (46) Karton, A.; Sylvetsky, N.; Martin, J. M. L. W4-17: A Diverse and High-Confidence Dataset of Atomization Energies for Benchmarking High-Level Electronic Structure Methods. J. Comput. Chem. 2017, 38, 2063−2075. (47) Ř ezác,̌ J. Cuby: An Integrative Framework for Computational Chemistry. J. Comput. Chem. 2016, 37, 1230−1237. (48) Ř ezác,̌ J. Cuby 4, Software Framework for Computational Chemistry, 2015. http://cuby4.molecular.cz/ (accessed Sept 19, 2017). (49) Aradi, B.; Hourahine, B.; Frauenheim, T. DFTB+, a Sparse Matrix-Based Implementation of the DFTB Method. J. Phys. Chem. A 2007, 111, 5678−5684. (50) DFTB+. http://www.dftb-plus.info/ (accessed Sept 19, 2017).

(10) Rapacioli, M.; Spiegelman, F.; Talbi, D.; Mineva, T.; Goursot, A.; Heine, T.; Seifert, G. Correction for Dispersion and Coulombic Interactions in Molecular Clusters with Density Functional Derived Methods: Application to Polycyclic Aromatic Hydrocarbon Clusters. J. Chem. Phys. 2009, 130, 244304−244304−10. (11) Grimme, S. Towards First Principles Calculation of Electron Impact Mass Spectra of Molecules. Angew. Chem., Int. Ed. 2013, 52, 6306−6312. (12) Ř ezác,̌ J.; Fanfrlík, J.; Salahub, D.; Hobza, P. Semiempirical Quantum Chemical PM6Method Augmented by Dispersion and HBonding Correction Terms Reliably Describes Various Types of Noncovalent Complexes. J. Chem. Theory Comput. 2009, 5, 1749− 1760. (13) Korth, M.; Pitoňaḱ , M.; Ř ezác,̌ J.; Hobza, P. A Transferable HBonding Correction for Semiempirical Quantum-Chemical Methods. J. Chem. Theory Comput. 2010, 6, 344−352. (14) Korth, M. Third-Generation Hydrogen-Bonding Corrections for Semiempirical QM Methods and Force Fields. J. Chem. Theory Comput. 2010, 6, 3808−3816. (15) Ř ezác,̌ J.; Hobza, P. Advanced Corrections of Hydrogen Bonding and Dispersion for Semiempirical Quantum Mechanical Methods. J. Chem. Theory Comput. 2012, 8, 141−151. (16) Domínguez, A.; Niehaus, T. A.; Frauenheim, T. Accurate Hydrogen Bond Energies within the Density Functional Tight Binding Method. J. Phys. Chem. A 2015, 119, 3535−3544. (17) Christensen, A. S.; Elstner, M.; Cui, Q. Improving Intermolecular Interactions in DFTB3 Using Extended Polarization from Chemical-Potential Equalization. J. Chem. Phys. 2015, 143, 084123. (18) Miriyala, V. M.; Ř ezác,̌ J. Description of Non-Covalent Interactions in SCC-DFTB Methods. J. Comput. Chem. 2017, 38, 688−697. (19) Gaus, M.; Cui, Q.; Elstner, M. DFTB3: Extension of the SelfConsistent-Charge Density-Functional Tight-Binding Method (SCCDFTB). J. Chem. Theory Comput. 2011, 7, 931−948. (20) Gaus, M.; Goez, A.; Elstner, M. Parametrization and Benchmark of DFTB3 for Organic Molecules. J. Chem. Theory Comput. 2013, 9, 338−354. (21) Gaus, M.; Lu, X.; Elstner, M.; Cui, Q. Parameterization of DFTB3/3OB for Sulfur and Phosphorus for Chemical and Biological Applications. J. Chem. Theory Comput. 2014, 10, 1518−1537. (22) Kubillus, M.; Kubař, T.; Gaus, M.; Ř ezác,̌ J.; Elstner, M. Parameterization of the DFTB3Method for Br, Ca, Cl, F, I, K, and Na in Organic and Biological Systems. J. Chem. Theory Comput. 2015, 11, 332−342. (23) Grimme, S.; Antony, J.; Ehrlich, S.; Krieg, H. A Consistent and Accurate Ab Initio Parametrization of Density Functional Dispersion Correction (DFT-D) for the 94 Elements H-Pu. J. Chem. Phys. 2010, 132, 154104. (24) Elstner, M.; Porezag, D.; Jungnickel, G.; Elsner, J.; Haugk, M.; Frauenheim, T.; Suhai, S.; Seifert, G. Self-Consistent-Charge DensityFunctional Tight-Binding Method for Simulations of Complex Materials Properties. Phys. Rev. B: Condens. Matter Mater. Phys. 1998, 58, 7260. (25) Yang; Yu, H.; York, D.; Cui, Q.; Elstner, M. Extension of the Self-Consistent-Charge Density-Functional Tight-Binding Method: Third-Order Expansion of the Density Functional Theory Total Energy and Introduction of a Modified Effective Coulomb Interaction. J. Phys. Chem. A 2007, 111, 10861−10873. (26) Finnis, M. W.; Paxton, A. T.; Methfessel, M.; van Schilfgaarde, M. Crystal Structures of Zirconia from First Principles and SelfConsistent Tight Binding. Phys. Rev. Lett. 1998, 81, 5149−5152. (27) Giese, T. J.; York, D. M. Density-Functional Expansion Methods: Grand Challenges. Theor. Chem. Acc. 2012, 131, 1145. (28) Kaminski, S.; Giese, T. J.; Gaus, M.; York, D. M.; Elstner, M. Extended Polarization in Third-Order SCC-DFTB from ChemicalPotential Equalization. J. Phys. Chem. A 2012, 116, 9131−9141. (29) Sedlak, R.; Ř ezác,̌ J. Empirical D3 Dispersion as a Replacement for Ab Initio Dispersion Terms in Density Functional Theory-Based M


Article

Journal of Chemical Theory and Computation (51) TURBOMOLE v6.6; 2015. (52) Stewart, J. J. P. Optimization of Parameters for Semiempirical Methods V: Modification of NDDO Approximations and Application to 70 Elements. J. Mol. Model. 2007, 13, 1173−1213. (53) Stewart, J. J. P. Optimization of Parameters for Semiempirical Methods VI: More Modifications to the NDDO Approximations and Re-Optimization of Parameters. J. Mol. Model. 2013, 19, 1−32. (54) Stewart, J. J. P. MOPAC 2016; 2016. (55) Scholten, M. Semiempirische Verfahren Mit Orthogonalisierungskorrekturen: Die OM3Methode, Ph.D. Thesis, Universität Düsseldorf, Germany, 2003. (56) Risthaus, T.; Grimme, S. Benchmarking of London DispersionAccounting Density Functional Theory Methods on Very Large Molecular Complexes. J. Chem. Theory Comput. 2013, 9, 1580−1591. (57) Grimme, S.; Bannwarth, C.; Shushkov, P. A Robust and Accurate Tight-Binding Quantum Chemical Method for Structures, Vibrational Frequencies, and Noncovalent Interactions of Large Molecular Systems Parametrized for All Spd-Block Elements (Z = 1−86). J. Chem. Theory Comput. 2017, 13, 1989. (58) Sure, R.; Grimme, S. Corrected Small Basis Set Hartree-Fock Method for Large Systems. J. Comput. Chem. 2013, 34, 1672−1685. (59) Neese, F. The ORCA Program System. WIREs Comput. Mol. Sci. 2012, 2, 73−78. (60) Bondi, A. Van Der Waals Volumes and Radii. J. Phys. Chem. 1964, 68, 441−451. (61) Petraglia, R.; Corminboeuf, C. A Caveat on SCC-DFTB and Noncovalent Interactions Involving Sulfur Atoms. J. Chem. Theory Comput. 2013, 9, 3020. (62) Hostaš, J.; Ř ezác,̌ J.; Hobza, P. On the Performance of the Semiempirical Quantum Mechanical PM6 and PM7Methods for Noncovalent Interactions. Chem. Phys. Lett. 2013, 568−569, 161−166. (63) Weber, W.; Thiel, W. Orthogonalization Corrections for Semiempirical Methods. Theor. Chem. Acc. 2000, 103, 495−506. (64) Kruse, H.; Grimme, S. A Geometrical Correction for the Interand Intra-Molecular Basis Set Superposition Error in Hartree-Fock and Density Functional Theory Calculations for Large Systems. J. Chem. Phys. 2012, 136, 154101. (65) Grimme, S.; Ehrlich, S.; Goerigk, L. Effect of the Damping Function in Dispersion Corrected Density Functional Theory. J. Comput. Chem. 2011, 32, 1456−1465.

N


Empirical Self-Consistent Correction for the Description of Hydrogen

Recommend Documents