Investigations of stacked DNA base-pair steps: highly-accurate

Subscriber access provided by University of Winnipeg Library

Quantum Electronic Structure

Investigations of stacked DNA base-pair steps: highly-accurate stacking interaction energies, energy decomposition and many-body stacking effects Holger Kruse, Pavel Banáš, and Jiri Sponer J. Chem. Theory Comput., Just Accepted Manuscript • DOI: 10.1021/acs.jctc.8b00643 • Publication Date (Web): 29 Nov 2018 Downloaded from http://pubs.acs.org on November 30, 2018

Just Accepted “Just Accepted” manuscripts have been peer-reviewed and accepted for publication. They are posted online prior to technical editing, formatting for publication and author proofing. The American Chemical Society provides “Just Accepted” as a service to the research community to expedite the dissemination of scientific material as soon as possible after acceptance. “Just Accepted” manuscripts appear in full in PDF format accompanied by an HTML abstract. “Just Accepted” manuscripts have been fully peer reviewed, but should not be considered the official version of record. They are citable by the Digital Object Identifier (DOI®). “Just Accepted” is an optional service offered to authors. Therefore, the “Just Accepted” Web site may not include all articles that will be published in the journal. After a manuscript is technically edited and formatted, it will be removed from the “Just Accepted” Web site and published as an ASAP article. Note that technical editing may introduce minor changes to the manuscript text and/or graphics which could affect content, and all legal disclaimers and ethical guidelines that apply to the journal pertain. ACS cannot be held responsible for errors or consequences arising from the use of information contained in these “Just Accepted” manuscripts.

is published by the American Chemical Society. 1155 Sixteenth Street N.W., Washington, DC 20036 Published by American Chemical Society. Copyright © American Chemical Society. However, no copyright claim is made to original U.S. Government works, or works produced by employees of any Commonwealth realm Crown government in the course of their duties.

Page 1 of 43 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Journal of Chemical Theory and Computation

Investigations of stacked DNA base-pair steps: highly-accurate stacking interaction energies, energy decomposition and many-body stacking effects Holger Kruse,1,* Pavel Banas,1,2 Jiri Sponer1,2

1Institute

of Biophysics of the Czech Academy of Sciences, Královopolská 135, 612 65 Brno, Czech Republic

2Regional

Centre of Advanced Technologies and Materials, Department of Physical Chemistry, Faculty of Science, Palacky University, 17. listopadu 12, 77146 Olomouc, Czech Republic

*corresponding author: [email protected]

Abstract The stacking energies of ten unique B-DNA base-pair steps were calculated with highly-accurate quantum chemistry and used as reference values in a thorough benchmark of (dispersion-corrected) DFT, wavefunction methods, tight-binding methods, and different force fields, including charge-variants thereof. The reference values were computed using a focal-point energy function based on extrapolated explicitly-correlated MP2-F12 and conventional CCSD(T) data at the triple-ζ level. A collection of 29 different density functionals, sometimes with multiple dispersion-corrections (D3(BJ), D3M(BJ) and VV10) were evaluated, including recent functionals like B97M-V, ωB97M-V and SCAN-D3(BJ), which perform excellently. The double-hybrid DSD-BLYP-NL was found to be the best DFT method. Common wavefunction methods (MP2, SCS-MP2 and MP2.5) and the SNS-MP2 protocol were tested as well, where only the latter and DLPNO-CCSD(T)/CBS were competitive with DFT. The tight-binding methods DFTB3-D3 and GFN-xTB revealed a comparatively low-accuracy. The AMBER force field outperformed CHARMM and GROMOS but still showed systematic gas-phase overbinding, which could be traced ACS Paragon Plus Environment

Journal of Chemical Theory and Computation 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

back to the electrostatic term, as revealed by comparison of different sets of point-charges. High-order SAPT, e.g., SAPT2+3δ(MP2) was not only benchmarked, but also used to study the nature of the stacking interactions to high accuracy. The δ(MP2) term turned out to be crucially important to reach high-accuracy. Finally, we investigated four-body stacking effects with DLPNO-CCSD(T) and DFT, which were found to be significant and strongest for the CpC base-pair step where they reached almost 30% of the total stacking energy.

Introduction Non-covalent interactions play a crucial role in biomolecular systems1–4 and their accurate description is thus essential for molecular modelling attempts with predictive power. In nucleic acids, e.g., one of the dominant non-covalent interactions is base stacking. The interaction of stacked aromatic moieties in particular is governed by strong London-type dispersion interactions.5–8 For nucleic acids systems where hetero-aromatic nucleobases are fundamental buildings blocks one often finds that charge-cloud Coulomb interactions and chargepenetration effects play a strong role as well.9–13 Numerous quantum-chemical (QM) studies have investigated the interaction strengths and types in various dimer complexes, analyzing the nature and magnitude of the direct (intrinsic) interactions between the aromatic rings.14–23 The ultimate contribution of base stacking to thermodynamic stabilities of nucleic acids then arises from an intricate interplay between the intrinsic stacking interactions and many other (free) energy contributions.4,24 The structural and stabilization roles of base stacking in nucleic acids thus cannot be straightforwardly predicted from knowledge of the intrinsic interactions alone, i.e., from evaluation of the potential energy surface (PES) of base stacking. Still, an accurate description of the intrinsic stacking interactions is crucially important for the understanding of the basic principles of molecular interactions in nucleic acids and any type of quantitative modelling of nucleic acids. In the past, accurate QM calculations of base stacking were usually done for nucleobase dimers5,7,32–35,24–31 and can nowadays also be found in countless studies for quantum chemistry benchmarking of non-covalent interactions (for benchmarking see refs. 14,15). The first electron-correlation calculations on base stacking were reported almost 25 years ago and were executed at the MP2 level with small diffuse-polarized basis sets.5,26,27,36 These calculations quite unexpectedly demonstrated that in the first approximation base stacking can be reasonably well modelled by an appropriate Lennard-Jones van der Waals potential combined with a Coulombic term using atom-centered electrostatic-potential-fitted (ESP) atomic charges.5 This has been a major justification for the emerging field of explicit-solvent molecular dynamics (MD) simulations of nucleic acids. The surprising ACS Paragon Plus Environment

Page 2 of 43

Page 3 of 43 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


success of a simple molecular mechanics (MM) force-field description of base stacking can be attributed to the fact that, when decomposing the base stacking interaction into physical terms, their variations on the PES compensate each other to a large degree.12,28 This does not suggests that the MM description of base stacking is flawless. We now know that the early QM calculations were far from being converged, which has motivated many subsequent studies refining the QM picture of base stacking.9,10,22,23,37–40 While (small) dimer complexes are ideal models for accurate analyses of the magnitude of base stacking, benchmarking of QM methods and fine-tuning of empirical potentials (MM force fields), the correct description of nucleic acids requires accurate calculations of interactions between multiple neighboring moieties at varying distances. Here, many-body effects, i.e., nonadditivity of stacking, could play an important role.41–43 A number of studies have appeared over the recent years that investigated not just stacked nucleobases, but stacked base-pairs, i.e., systems with four nucleobases.29,37,44 In a B-DNA double helix, two consecutive stacked base pairs are known as base-pair steps. There are ten unique base-pair step sequences in B-DNA. The base-pair step stacking interactions are assumed to play a significant role in inducing a sequence-dependent variability of the B-DNA double helix conformation,45 a phenomenon which is, despite years of intense investigations, far from being satisfactorily understood.9 Standardly, QM benchmarking for non-covalent interactions is done for structures at equilibrium distances, partly because it is the chemically most important conformation, partly because of habits. However, from a methodological point of view it is equally important to look at the whole interaction energy potential, including geometries with unoptimized molecular separation, to capture the asymptotic decay of non-covalent interactions, which is vital in the treatment of large molecules. Long-range correlation effects such as London dispersion interactions play an important role for essentially all flexible molecules larger than a dozen atoms.46 For DFT this requires the correct asymptotic decay of either the xc-kernel, which is typically not achievable as the density decays exponentially, or the use of an additional dispersion correction.47 Benchmarking for the S66x8 benchmark set of biomolecular non-covalent interactions with multiple intermolecular distances clearly showed that an incorrect asymptotic description of dispersion leads to large errors.20,48 Even for carefully-tuned modern functionals like M06-2X that can describe dispersion bound complexes at the equilibrium distance well enough it is important to include a long-range dispersion correction.20,48 An often-overlooked feature of spin-scaled theories (SCS-MP2 and variants) is the resulting incorrect asymptotic behavior for intermolecular dissociation due to the different scaling of singlet and triplet electron pairs.49 Moreover, spin-scaled theories are often fitted

ACS Paragon Plus Environment


to a certain electron correlation length, e.g. at equilibrium distances of small- to medium sized molecular complexes, adding another empirical component affecting the description of the potential-energy curves. These aspects will play a role if scans for the base-pair step stacking geometries are performed to elucidate the dependencies of stacking energies on helical parameters of the double helix.29,37 In other words, while it is relatively straightforward today to perform reliable semi-quantitative QM calculations on base stacking, achieving a quantitative accuracy over the whole PES is still not trivial. In this work, we present a thorough quantum chemical benchmark study for natural base-pair step stacking geometries as found in B-DNA. We first describe the base-pair step geometries selected for this study. Then, we discuss our choice of benchmark reference method for the accurate computation of pair-wise stacking energies. Following, we benchmark a large variety of dispersion-corrected density functional approximations (DFAs), variants of MP2 and MP2-F12, few high-order SAPT methods, and selected low-cost and semi-empirical methods. We focus particularly on the benchmarking of the inter- and intra-strand stacking energies. The stacking energies of the full base-pair steps are investigated as well and compared to previous studies. Alongside, many-body effects in the stacking interactions are evaluated. Finally, the nature of the inter- and intra-strand stacking energies is revealed to high accuracy with modern SAPT methods. Based on the results, among the general conclusions, we provide recommendations for the selection of methods for various levels of accuracy for nucleobase stacking calculations.

Computational Details

Figure 1. Nomenclature and highlighting of stacking interaction pairwise terms using the CpG base-pair step (stack) with canonical base pairs as an example. Note that we use the standard abbreviation NpM (N,M is any nucleobase; the base-pair step sequence is determined either as 3p1 or 2p4) to mark the base-pair steps despite excluding the sugar-phosphate backbone from the computations.


Page 4 of 43

Page 5 of 43 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Definition of interaction energies. Figure 1 shows the base stacking between two base pairs. The two-body stacking energies for the base-pair steps ΔEstack are calculated by summing the four pair-stacking energies as depicted in Figure 1 and can be written as follows: ΔEstack = ΔE13 + ΔE24 + ΔE14 + ΔE23 where ΔEXY denotes the interaction energy between the nucleobases X and Y. 𝑋,𝑌

𝛥𝐸𝑋𝑌 = 𝐸𝑋𝑌 ―

∑𝐸

𝑖

𝑖

The four-body contribution ΔEABCD is the difference between the four-body and the two-body stacking energy (using the pair-approximation for many-body terms) 𝛥𝐸4𝑠𝑡𝑎𝑐𝑘 = 𝐸1234 ― 𝐸12 ― 𝐸34 𝛥𝐸𝐴𝐵𝐶𝐷 = 𝛥𝐸4𝑠𝑡𝑎𝑐𝑘 ― 𝛥𝐸𝑠𝑡𝑎𝑐𝑘 Here ΔE4stack indicates the inclusion of the four-body interaction term. For the nucleobases adenine (A), guanine (G), cytosine (C) and thymine (T) there are 10 symmetry-unique base-pair step stacks involving canonical CG and AT base pairs.28 To use familiar biomolecular language we speak of, e.g., CpG steps even though we only consider the nucleobases C and G and not the backbone. Elsewhere the notation CG:CG has been introduced to designate a GC base pair stacked upon a CG base-pairs.37 Comparing these nomenclatures leads to: ApA = AA:TT; ApC = AC:GT; ApG = AG:CT; ApT = AT:AT; CpC = CC:GG; CpG = CG:CG; GpC = GC:GC; TpA = TA:TA; TpC = TC:GA; TpG = TG:CA.27–29,36

Preparation of structures. The structures of the base-pair steps were prepared with an in-house software introduced in our recent study.9 This software is able to generate base-pair step geometries from given helical parameters and reference-frame geometries of AT and GC base pairs in a similar fashion to, e.g., the 3DNA program.50 Following the procedure in our earlier work,9 the reference geometries of the AT and GC base pairs used for this purpose were obtained by geometry optimization in the gas phase at the MP2/cc-pVTZ level of theory constrained to CS symmetry. The symmetry constraint is needed for idealized, planar base pairs, preventing a potential non-planarity due to the amino-group pyramidalization.36 The optimal orientation of the thymine methyl group was obtained by a potential energy surface scan at the same level of theory. For creating model B-DNA base-pair step geometries for the benchmarking we chose the 'natural' arrangement using a helical twist of 36° and a Propeller of 0°. The vertical separation (Rise) was optimized for each base-pair step by a



potential energy scan using cc-pV(D,T)Z extrapolated and counterpoise corrected SCS(MI)-MP2 as outlined in previous work.9 All geometries and the results from the vertical separation scans can be found in the Supporting Information. The 2006 set of geometries presented in ref. 29 by one of us can be considered similar to the presently used geometries, but not identical. They have also been constructed to create idealized B-DNA base-pair steps with a helical twist of 36°, using geometrical centers of the C6-C8 base pair axes as the x, y = 0, 0 coordinate points. However, some of the base-pair steps considered non-zero Propeller (optimized at the MM level) and also the vertical separation between the base pairs has been adjusted at the MM level, as in detail explained in ref.29 (cf. Table S10 in the Supporting Information). Further, some unavoidable differences in the stacking geometries between the present and the preceding study arise from the fact that the reference base-pair geometries in ref. 29

were obtained using HF/6-31G** gradient optimizations. Another set of geometries was presented by Parker and Sherrill10 in 2015 that use AT and GC base pairs

optimized at the B3LYP/aug-cc-pVDZ level of theory and a different convention for the reference frame [CEHS (Cambridge University Engineering Department Helix computation Scheme) versus 3DNA in our work]. Moreover, they use helical and base-pair parameters of an average B-DNA geometry, rather than the idealized ones as used in the presented work. In particular, this leads to a Propeller of -11° and potentially unoptimized vertical distance between the base pairs due to the fixed Rise value of 3.32 Å; it is generally recommended to relax the vertical distance between stacked base pairs for each combination of the remaining geometrical parameters.9,28,29 All together this makes an interaction energy comparison with our structures difficult. Numerical values are collected in the Supporting Information Table S9 for the interested reader. However, all three sets of structures (those from the current work, the 2006 and the 2015 ones) are valid geometries to study the energetics of basepair stacking, as they all correspond to the low-energy conformational region sampled by thermal motions in a B-DNA double helix.


Page 6 of 43

Page 7 of 43 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Figure 2. Top-down view of four selected base-pair step stacks showcasing the geometrical overlap (stacking) of the nucleobases. General technical notes. Only valence electrons are correlated in the wavefunction-based methods (frozen core approximation). Auxiliary basis sets for RI-J and RI-JK (DF-JK in PSI4 language) are used as provided by the respective programs. SCF energies were converged to at least 1e-7 Eh for DFT calculations. Additionally, density changes were converged to at least 1e-6 au for correlated methods. All DFT calculations, and wavefunction methods with CP in the name, use the respective supermolecular basis to remove basis set superposition error (BSSE) via the counterpoise correction (CP).51 Molecular visualization was done with Chimera.52 Basis sets. The abbreviation VXZ denotes the correlation-consistent cc-pVXZ basis sets. A lower case a indicates a full set of diffuses functions, e.g. aVTZ = aug-cc-pVTZ, while jun indicates Truhlar’s calendar basis notation for a reduced set of diffuse functions.53–55 Similarly, the abbreviation VTZ-F12 indicates the special-purpose F12 basis set cc-pVTZ-F12.56 Extrapolation using the Ahlrichs’ basis sets57 was done with the double-polarized sets def2TZVPP and def2-QZVPP (abbreviated to def2T and def2Q) using the exponents proposed by Neese.58 The following programs were employed in this work: PSI4. Except when explicitly noting other programs, all DFT and correlated calculations were performed with PSI4.59 An unpruned DFT quadrature grid with 590 angular and 75 radial points was employed. Density-fitting was applied throughout. The D3 correction was calculated through an interface to Grimme’s dftd3 program. No D3(ABC) correction was included for the pair-energy benchmark data. An ongoing development version of the PSI4 program and the official v.1.2.1 release was used. PSI4 v.1.2.1 interfaces the LibXC library60 version 4.0 for the implementation of all density functional equations. The DFT-NL scheme61 has been implemented in PSI4 by one of us (H.K.). SAPT2+3 and SAPT2+(3) plus higher-order coupling through supermolecular δ(MP2) and δ(HF) corrections (denoted e.g. as SAPT2+3δ(MP2) in ref.62) were performed as implemented in PSI4 using v1.2.1 defaults using the jun-cc-pVTZ basis set. The δ(HF) term accounts for higher-order induction effects and is beneficial for polar ACS Paragon Plus Environment


molecules. The δ(MP2) term accounts for higher-order coupling between induction and dispersion and is for PSI4 included in the induction part. The detailed description of the included terms and used notation herein can be found in works by Sherrill.62,63 The third-order dispersion terms have a CCSD-like scaling and thus a significant computational demand. However, various approximations such as density-fitting and virtual orbital cutoffs based on a natural orbital expansion64,65 make the computations feasible for routine application to stacked nucleobases. The employed SCS(MI)-MP2 variant uses fitted parameters for the basis set extrapolation level CBS(VDZ,VTZ).66 SAPT0 calculations were performed with the slightly augmented double-ζ basis (jun-cc-pVDZ) as suggested to exploit error compensation effects of its individual components.65 The SNS-MP2 protocol for interaction energies was computed with the provided plugin for PSI4.67 FNO in FNO-CCSD(T) stands for frozen natural orbitals that reduce the virtual orbital space with negligible loss of accuracy.17 Turbomole. The SCAN meta-GGA approximation and the low-cost composite methods PBEh-3c, HF-3c and B97-3c were computed with Turbomole v.7.2.1.68,69 The grid m4 was used and for SCAN the radial grid points were vastly increased using the keyword radsize 50 as recommended in the manual leading, e.g., to over 285 radial grid points for the elements K-Kr. For comparison, the default is only about 50 radial grid points for these elements. Density fitting was applied for the SCAN and B97-3c calculations. The composite methods employ the three-body D3(ABC) correction by default. Turbomole was also used to calculate the canonical CCSD(T) results for the CBS(T)-F12-CP method. ORCA 4. The (2015 version of) DLPNO-CCSD(T)70 and DLPNO-MP2-F12/D71 calculations were done with ORCA v4.0.1.72 The PNO threshold was increased with the keyword TightPNO.73 Extrapolations follow the ORCA internal procedures with optimized extrapolation coefficients.58 The SCF threshold was increased with the keyword TightSCF. Also, all canonical MP2-F1274 calculations were performed with ORCA using density fitting with optimized basis sets.75 Others. GFN-xTB76 was used with Grimme's xtb program using default settings. DFTB3-D3(BJ)77,78 was used with the dftb+ program79 and dftd3. PM6-D3H+80 was used with mopac16.81 Note that the H+ correction was not included (zero contribution to the stacking energy). Force fields. The Cornell et al. (AMBER) force field82 stacking energies were calculated using four different sets of charges (see Table S14 in the Supporting Information). The first set consisted of the charges taken directly from the force-field library, including the charges of the R(N9)/Y(N1) nitrogens that are part of the glycosidic bonds in the context of a DNA molecule, while their capping hydrogens R(H9)/Y(H1) were set to neutralize the


Page 8 of 43

Page 9 of 43 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


total charge of each nucleobase (N and Y stand for purine and pyrimidine, respectively). Although this approach yields charges that are imbalanced for the R(N9-H9)/Y(N1-H1) bonds, it resembles most closely the charges used in MD simulations. To indicate the origin of the charges we call this force-field version AMBER(library). The other three set of charges were obtained by the RESP fitting procedure with aid of Antechamber83 and Gaussian 09.84 The HF denoted set of charges [AMBER(HF)] was calculated by RESP fitting at the HF/6-31G(d) level of theory for monomers optimized at the same level, with a low density of grid points (four layers with one grid point per Bohr2). This approach corresponds to the settings used for the derivation of Cornell et al. force field and thus mimics the AMBER(library) charges while simultaneously providing balanced charges of the R(N9-H9)/Y(N1-H1) bonds due to charge equilibration via the RESP fitting procedure. Finally, the last two sets of charges were calculated by RESP fitting at the PBE0/def2-QZVP level of theory so that the MM Coulombic term was as close as possible to the QM electrostatic interaction. The first set of PBE0 charges, denoted PBE0stacks calculates the wavefunction and RESP charges for all individual pair-wise stacking geometries, i.e., we have carried out 40 separate charge-fitting calculations. This captures polarization effects for base dimers individually, but not for the whole base-pair steps. These charge sets are, at least formally, the most physically-accurate for a given stacking geometry. Since the charges are computed for each stacked nucleobase dimer, they are not constant. The second set of PBE0 charges is denoted PBE0AT,GC and uses fixed charges derived for isolated AT and GC Hbonded base pairs. This captures the H-bond polarization effects in the base pairs that are particularly pronounced in the CG base pair and are missing in the default AMBER charges. To compare the AMBER force-field stacking energies with those calculated by some other MM force fields, we performed additional calculations using the CHARMM3685 and Gromos 53a686 force fields. As in the case of AMBER(library) stacking energies, we calculated the CHARMM(library) and GROMOS(library) energies using the appropriate van der Waals parameters, combination rules and charges of the particular force field. As in the case of AMBER(library), the charges of capping R(H9)/Y(H1) hydrogens were set to neutralize the nucleobases; see the previous paragraph. CHARMM(HF) and GROMOS(HF) indicate calculations carried out using the same RESP HF/6-31G(d) charges as applied for the AMBER(HF) computations. Additional DFT details. The full list of functionals with respective citation reads: B2PLYP,87 B2GPPLYP,88 DSDBLYP,89 DSD-PBEB86,90 PTPSS,91 DSD-PBEB95,92 B3LYP,93,94 M06-2X,95 M05-2X,96 M11,97 SOGGA11X,98 MN15,99 PBE0,100 revPBE0,101 LC-VV10,102 ωB97M-V,103 PW6B95,104 LC-wPBE,105 CAM-B3LYP,106 BLYP,107,108 revTPSS,109 TPSS,110 B97-D,111 VV10,102 PBE,112 SCAN,113 BP86,107,114 MN15-L115 and B97M.116



An overview of the applied D3 and NL parameters is given in the Supporting Information Tables S5-S8. The NL or -V suffixes indicate dispersion correlation from the VV10 kernel102 developed for van der Waals density functionals like VV10 and LC-VV10. The conventional xc-part in VV10 is rPW86PBE, i.e., VV10 equals rPW86PBENL using the DFT-NL notation. In the literature VV10 is often used for both the density-dependent kernel calculating the dispersion contribution and for the functional name itself. We adhere to the names used in the LibXC where the functional name VV10 stands for rPW86PBE including the VV10 correlation kernel. The -NL correction is an a posteriori addition to the existing functional while for the -V functionals, like ωB97M-V, it was included in the functional fit from the beginning.103 We applied the VV10 kernel in a fully self-consistent manner. A recent systematic study, however, suggests that this is not necessary for accurate results.117 New VV10 parameter for the revTPSS functional. As standard for the -NL procedure the C parameter is fixed at 0.0093 and only the b parameter is varied. Using a line search procedure over a fit set consisting of the (revised) S66x848 set of biomolecular non-covalent interactions and the RG18 rare-gas set taken from the GMTKN55,15 we obtained for revTPSS-NL a b parameter of 5.465.

RESULTS & DISCUSSION Choice of the reference method. Inspired by the recent benchmarking efforts of Martin and coworkers19,48 we arrived at the following focal-point based energy function labeled CBS(T)-F12-CP to estimate the CCSD(T) basis set limit : 𝐸[𝐶𝐵𝑆(𝑇)–𝐹12–𝐶𝑃] = 𝐸[𝑀𝑃2–𝐹12/𝐶𝐵𝑆(𝑉𝑇𝑍–𝐹12,𝑉𝑄𝑍–𝐹12)–𝐶𝑃] + 𝛥𝐻𝐿[𝐶𝐶𝑆𝐷(𝑇)/𝑗𝑢𝑛–𝑉𝑇𝑍] where ΔHL stands for the high-level (HL) correction that accounts for high-order correlation effects by taking the differences between the indicated high-level in brackets and MP2. 𝛥𝐻𝐿 = 𝐸[𝐶𝐶𝑆𝐷(𝑇)] ― 𝐸[𝑀𝑃2] |𝑗𝑢𝑛–𝑉𝑇𝑍 The MP2-F12/CBS extrapolation employs the generalized formula by Martin: 𝐸(𝐶𝐵𝑆) = 𝐸(𝐿) +

𝐸(𝐿) ― 𝐸(𝐿 ― 1) 𝐿 𝛼 ( ) ―1 𝐿―1

with the exponents 7.6070 for the SCF and 4.3548 for the correlation energy (taken from ref 48) and L denoting the basis set ζ-level. Over the last years the consensus has been reached that the smallest basis set for the highlevel correction should be aug-cc-pVDZ, and even then larger outliers are to be expected.19,118 Thus, the general


Page 10 of 43

Page 11 of 43 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


recommendation is to use the largest affordable basis set, which for us was the augmented triple-ζ basis jun-ccpVTZ, which is very similar to a heavy-augmented version often used in the past and commonly showing robust results.14 Even using the memory-efficient implementation in Turbomole, the calculations for all 40 pairwise base stacks required considerable computational efforts. Lastly, the “CP” in CBS(T)-F12-CP indicates that the MP2-F12 calculations were performed in the supermolecular basis set to remove spurious basis set superposition error (BSSE) effects.119 The conventional MP2 and SCS-MP2 extrapolation using PSI4 employed the Halkier two-point extrapolation120 with α=1.63 and β=3.0, and using the highest basis set for the SCF value, following the automatized procedure within PSI4. A note on the expected accuracy. The intra-strand stacking interactions are expected to be more sensitive to the one-particle basis set error and level of theory than the inter-strand stacks since they have a larger density overlap between each other and are generally stronger. Methods to describe them need to accurately capture the strong London-type dispersion interactions. It is also well-known that stacking energies are prone to the BSSE, especially for correlated methods.119 The MP2 BSSE in particular, and basis set convergence in general, have been studied for benzene dimers, stacked nucleobases and other biomolecular dimers to great detail.48,119,121,122 Due to the interplay between BSSE and basis set incompleteness error (BSIE) one often finds that applying the CP correction for small to medium-sized basis sets is not accurate enough, especially for MP2, as BSSE and BSIE are opposing effects.121,122 Then a scaling factor is applied to the CP correction to empirically compensate for the basis set incompleteness (usually by a factor of 0.5). Extrapolations to the complete basis set limit, either by a fixed exponent of 3 or by an empirical fitted exponent for a specific basis set combination, reduce this problem, but were designed to fix the basis set incompleteness error (BSIE) rather than the BSSE. Note that the separation of the general basis set error into BSIE and BSSE is a conceptual construct and thus often arbitrary. Combining a fitted CBS extrapolation with a scaled CP correction, however, leads to a loss of the systematic improvement of the one-particle basis expansion and less predictable results due to the introduction of additional empirical elements. A recent trend to overcome this problem are F12-based explicitly-correlated methods that inherently reduce the BSIE.123 A positive side effect of applying F12 methods is the use of the specialized F12 basis sets ccpVXZ-F12 that are larger compared to their cc-pVXZ companions due to an increase in the amount of valence basis functions. The cc-pVDZ-F12 basis set, e.g., essentially possesses a quadruple-ζ level for the valence orbitals, but only a double-ζ polarization set. Both the larger basis and the F12 correction tremendously reduce the BSIE



and the CBS extrapolation quickly leads to the "true" basis set limit. After extensive testing Martin and co-workers advocate to apply the full CP correction on MP2-F12/CBS results in their study of the S66x8 benchmark set.48 The general recommendation for a ΔHL correction is to use the largest basis set possible, as mentioned above. A CP correction and CBS extrapolation are generally less important here as the difference between MP2 and CCSD(T) data is less sensitive to basis set errors.118 F12 methods are also recommended for the ΔHL correction,19 but they introduce complications with the perturbative triples correction as the standard (T) contributions contain no direct F12 correction creating an imbalance between CCSD-F12 and (T) that is typically solved by scaling the (T) correction by some scheme.48,124,125 Herein we decided to apply canonical CCSD(T)/jun-cc-pVTZ calculation, which is around the upper limit of calculations that can be routinely performed for stacked nucleobases on typical academic research clusters. We assume that the main source of inaccuracy in our selected reference method lies in remaining basis set error of MP2-F12/CBS, where a small remaining BSIE for the intra-strand stacking energies still interferes. The ΔHL correction carries an unknown, but likely even smaller amount of error. Conservatively, we expect an error of 1-2% for the intra-strand stacking energies. Given an average interaction energy of about 6 kcal/mol this would amount to 0.06 - 0.12 kcal/mol. For comparison we also computed full-extrapolated DLPNO-CCSD(T)/CBS(def2T,def2Q) results without CP correction. These calculations are in good agreement with our reference and give a percentage deviation for the more difficult intra-strand stacking of 3.5% (cf. Table 4 in the wavefunction benchmark section). This result is encouraging from two perspectives: first, the DLPNO-CCSD(T)/CBS calculations are significantly less computationally demanding than the canonical CCSD(T)/jun-cc-pVTZ calculations employed in the CBS(T)-F12-CP scheme, making them an excellent alternative reference method for even larger systems, and second, the treatment of the higher-order correlation effects in the extrapolated DLPNO-CCSD(T)/CBS(3,4) scheme is more rigorous than the ΔHL correction at the jun-cc-pVTZ level. The DLPNO-CCSD(T)/CBS calculations, while capturing less correlation energies than the conventional treatment by construction with the same basis set, have the advantage of being able to use larger basis sets. New developments using a PNO-based full (T) correction over the currently in ORCA available semi-canonical (T0) version promise to capture even higher amounts of canonical correlation energy.126–128 In our experience, the updated DLPNO implementation70 of 2015 also seems superior in handling delocalized systems such as nucleobases than the older 2013 DLPNO version,129 which can also be ascribed to an overall improved and more robust localization procedure. Moreover, the DLPNO-CCSD(T)/CBS method has also been used to provide reference values for various subsets of the GMTKN55 benchmark study.15


Page 12 of 43

Page 13 of 43 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


From the above points we are confident that our DLPNO-CCSD(T)/CBS results can reliably serve as a benchmark method for computing the four-body interaction term needed to study non-additive effects in the full base-pair step stacking complexes.

Benchmarking of pair-stacking energies. In the following we will discuss the new reference pair-stacking energies ΔEstack for the ten B-DNA base-pair step stacks and compare them with some earlier results. We wish to point out that our comparison of MM and QM methods is slightly biased against the MM data for gas-phase calculations, when the electrostatic term of the force field is used without any modification. Most biomolecular MM parametrizations aim to reproduce structures in explicit-solvent simulations and thus often use scaled point-charges, or deliberately over-polarized charges like the HF/6-31G(d) ESP-fitted charges for the AMBER force field; the reason is to mimic typical polarization of the solute electronic structure by chargescreening effects from polar bulk solvation and also to pair with the pair-additive explicit-solvent water models.2,82,130 Typical MM potentials should be understood as 'effective' potentials whose errors for a particular type of interaction, or the complete absence of energy terms like e.g. polarization, are balanced implicitly by the parametrization of the available force-field terms.12 Even though the atom-pair wise potentials are strictly additive, the overall parametrization of a biomolecular force field tries to capture many-body effects implicitly. The utilization of an unmodified biomolecular force field for small isolated gas-phase dimers during benchmarking can thus be considered a 'worst case' scenario as it is, regarding the coulombic term, furthest from the target parametrization.5 However, even such computations are important, as too large or systematic deviations between QM and MM would inevitably transfer to errors in the description of the true target systems in the explicit-solvent simulations. Nevertheless, it has been suggested that, in order to best illustrate the capability of the pair-additive force field to capture base stacking, it is advisable to reparametrize the point charges at a computational level which is as close as possible to the QM benchmark method.5 In Table 1 we report our new reference interaction energies and compare them with several other calculations. The full set of reference values for all individual pair-stacking energies can be found in the Supporting Information. Our new benchmark is listed in the first data column while the second column shows our alternative benchmark DLPNO-CCSD(T)/CBS which achieves very consistent results with the CBS(T)-F12-CP reference. Moreover, we



include two methods that are among the best performers in our benchmark study (vide infra). The highly-reliable methods SNS-MP2 and DSD-BLYP-NL in columns five and six provide results very close to the benchmark. Their respective statistical performance is discussed in detail in later parts of this work. The third column of Table 1 shows that the AMBER(library) force field on average over-estimates the strength of gas phase base stacking. Similar results were found by Parker et al. for both AMBER and CHARMM using Sherrill's 'silver standard' method based on dispersion-weighted CCSD(T)-F12.10 However, the AMBER(HF) charge re-fit (column four) has seemingly no such systematic error and can lead to both over- and underestimation of the stacking energies. The right side of the Table shows the decade-old benchmark energies29 together with a recalculation of the energies using the 2006 geometries and our best DFA method DSD-BLYP-NL. This was done to enable an easier comparison between the two different sets of geometries. Comparing both DSD-BLYP-NL columns one recognizes non-negligible differences. They can be caused by the fact that the 2006 geometries have considered a non-zero Propeller twist in some base-pair steps, they have differently adjusted the vertical separation and the base-pair geometries were obtained only at the HF-quality (cf. Table S10, Supporting Information). This, however, has no effect on using both geometry sets for a benchmark to compare different methods. Thus, the 2006 CBS(T) data shows more negative (stabilizing) stacking compared to the DSD-BLYP-NL method which indicates that the CBS(T) method is somewhat over-stabilizing the base stacking (considering the close correspondence between the DSDBLYP-NL and CBS(T)-F12-CP methods, cf. columns six and one). For the sake of completeness, the interaction energies for the H-bonded AT and GC dimers are −16.87 and −32.12 kcal/mol at the CBS(T)-F12-CP level of theory.

Table 1. Explicit sum pair-stacking energies (i.e., base-pair step stacking) ΔE stack including the new CBS(T)-F12-CP reference values. base pair stack

CBS(T)-F12-CP

DLPNO-CCSD(T)a

AMBER(library)

AMBER(HF) SNS-MP2 DSD-BLYP-NL CBS(T)b DSD-BLYP-NL 2006 geometriesc ApA AA:TT −13.37 −13.16 −14.72 −13.52 −13.37 −13.59 −13.1d −12.65 ApC AC:GT −13.33 −13.24 −14.35 −13.60 −13.33 −13.50 −14.2 −13.39 ApG AG:CT −13.96 −13.67 −14.76 −12.95 −14.11 −14.05 −14.4 −13.83 ApT AT:AT −11.66 −11.62 −13.36 −12.52 −12.00 −11.97 −13.3 −12.89 CpC CC:GG −13.81 −13.50 −13.65 −12.04 −13.03 −13.84 −13.7 −13.18 CpG CG:CG −18.44 −18.11 −19.50 −16.28 −18.57 −18.34 −18.4 −17.55 GpC GC:GC −16.14 −15.77 −16.51 −15.57 −16.14 −16.07 −16.6 −15.61 TpA TA:TA −13.74 −13.58 −16.13 −14.43 −14.17 −13.85 −13.0 −12.51 TpC TC:GA −13.22 −12.90 −14.27 −13.25 −13.06 −13.29 −13.6 −12.82 TpG TG:CA −15.95 −15.73 −18.25 −15.94 −15.92 −16.02 −16.0 −15.31 a CBS(def2T,def2Q); b MP2/CBS(aVDZ,aVTZ) + ΔHL[CCSD(T)/6-31G*(0.25)] data taken from Table 4 in ref. 29 no four-body term; c calculated on the 2006 structures from ref. 29; d geometry with 0 Propeller has been used.


Page 14 of 43

Page 15 of 43 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Separate Evaluation of Intra- and Inter-strand Base Stacking Pairs. As already noted by Parker et al.,37 benchmarking the sum stacking energies ΔEstack is susceptive to error cancellation effects between the individual stacking energies. An equally accurate treatment of the intra- and inter-strand stacks is important for a balanced description and perhaps the most difficult aspect for this benchmark. The interaction energies and interaction nature between the intra- and inter-strand bases are substantially different. The average interaction energies are −6.0 kcal/mol for the intra-strand stacks and −1.1 kcal/mol for the inter-strand stacks. The intra-strand stacks have a stronger London-dispersion interaction component than the inter-strand stacks (vide infra), especially considering the ΔE23 term as the helical twist moves the corresponding bases away from each other. Because of possible error compensation and enhancement effects while summing up all four stacking energies for each base-pair step stack we think that statistical benchmarking for the method evaluation is more justified for the individual dimer base stacking energies. This further allows to evaluate the inter- and intra-strand stacks separately. The statistical analysis of the pair stacking energies is listed in Table 2. For the statistics over all base stacks we use the mean absolute deviation (MAD) as it is an established measure in the QM community. The additional listing of the mean deviation (MD) serves the purpose to reveal any systematic over- or underbinding of the methods. The MD is calculated as ΔEmethod – ΔEref, so that a negative MD indicates overbinding while a positive MD underbinding. Moreover, we report the maximum absolute error (|MAX|). For the analysis of the inter- and intra-strand stacking energies we use the mean absolute percentage error (MAPD, cf. ref 20) that takes the magnitude of interaction into account. The respective intra-strand and inter-strand MAD and MDs are listed in the Supporting Information Tables S11-S13. The percentage-based MAPD evaluation enables us to compare the errors for the inter and intra-strand stacks with their different average stacking energies on a comparable basis. Table 2 is ordered according to Jacob's ladder131 starting with double-hybrid density functional approximations (DHDFA) as the highest rung. These DHDFAs incorporate a perturbative treatment of the virtual Kohn-Sham energies to obtain correlation energy, analogous to MP2 theory (PT2) but without singles excitations.87 The correlation part of the exchange-correlation kernel (xc) is accordingly scaled down as originally proposed. A recent study by Hait and Head-Gordon denoted this approach tDH (truncated double-hybrid) and showed that it is the more robust approach compared with the alternative xDH, e.g. used in XYG3.132 Only double-hybrids of the tDH class are tested.



The number of available DFAs in the literature is vast and some selection has to be done for practical matters. We chose functionals that i) are historically performing well for non-covalent interactions in benchmarks, vii) are commonly encountered in application work for nucleic acids and biomolecular systems or iii) have been only recently introduced and more benchmarking experience is necessary. Of course, the functional selection is also technically limited to what is implemented in the quantum chemistry programs available within our group. Double-Hybrid DFA. Best-performing DHDFAs are DSD-BLYP-NL and DSD-PBEP86-NL, with MADs of only 0.08 and 0.09 kcal/mol, respectively, and the only DFAs tested herein with an MAD below 0.10 kcal/mol. While they are also best-performing for the inter-strand stacking, the perhaps more important intra-stand stacking is best described by B2PLYP-D3M(BJ) with an MAPD of 1.7 %, although DSD-PBEP86-D3(BJ) follows very closely with an MAPD of 1.9 % (note that all MAPD values discussed in the text are either for intra- or inter-strand stacks). All tested double-hybrids except PTPSS-D3(BJ) and B2GPPLYP-D3(BJ) show excellent MADs below 0.3 kcal/mol and can generally be recommended for the calculation of nucleic acids base stacking. The latter two functionals show a systematic underbinding. Here one can imagine that a re-fit of the D3 correction might improve the results, similarly as D3M(BJ) improves the results for B2PLYP. (meta-)hybrid DFA. Global or range-separated hybrid functionals include a percentage of exact Fock-exchange. While this (partially) solves many problems of GGA functionals like the delocalization or many-electron selfinteraction error,133 it comes with an increased computational cost making their application to large biomolecular fragments expensive with large basis sets. The hybrid DFAs tested herein show a larger variance in their performance than the DHDFAs, partially because some DFAs like M06-2X are tested also without additional dispersion correction, as they were initially designed to incorporate dispersion in their functional form.95 However, also the missing PT2 correlation causes a generally reduced robustness of the DFAs. Best performing DFAs are B3LYP-D3(BJ), B3LYP-NL, revPBE0-NL and ωB97M-V, all of which stay below an MAD of 0.2 kcal/mol and are thus comparable in performance to double-hybrid DFAs. A D3(BJ) variant of ωB97M-V (ωB97M-D3(BJ) was recently proposed, in which the D3(BJ) replaces the -V correction without re-fitting the functional parameters.117 However, for this benchmark it shows an underbinding tendency and yields only an MAD of 0.35 kcal/mol. Other well-performing functionals, but slightly above an MAD of 0.2 kcal/mol, are LC-ωPBE-D3(BJ) and LC-VV10, both long-range corrected hybrids. All the best-performing functionals have in common that they perform around or below an MAPD of 5% for the intra-strand stacking energies. It has been recognized that many of Minnesota-type functionals that were designed to include dispersion in their functional forms are only able to


Page 16 of 43

Page 17 of 43 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


do so for regions of overlapping densities as found in equilibrium distances of small dimers.20,48 Long-range dispersion requires an add-on correction that D3 can readily provide.134 Then the M05-2X, M06-2X, M11 and MN15 functionals are significantly improved from MADs sometimes above 1 kcal/mol with heavy underbinding to a reasonable performance with MADs from 0.30 to 0.68 kcal/mol. The improvement is especially large for the intra-strand stacking. For SOGGA11-X a dispersion correction is mandatory, as indicated by the uncorrected MAD of over 4 kcal/mol. SOGGA11-X-D3(BJ) results are largely improved; an MAD of 0.49 kcal/mol is obtained over all base stacks, while the intra-strand base stacks are well described with an MAPD of 5.4%. However, the interstrand performance of only 4.5% indicates some problems that worsen the overall MAD, otherwise the functional would be rated among the top performing hybrid DFAs. The MN15 functional, a new type of functional based on NGA (non-separable gradient approximation)99 instead of GGA formalism, is a global-hybrid meta-NGA that performs merely acceptably without dispersion correction, giving an MAD of 0.63 kcal/mol and an MAPD of 11% for intra-strand stacks. MN15-D15-D3(BJ) significantly improves the performance to a reasonable MAD of 0.30 kcal/mol. A more in-depth study of some recent Minnesota functionals, and why for the newer functionals the BJ-damping form of the D3 correction is recommended can be found in the work by Goerigk.134 (meta-)GGA/NGA. The performance of several (meta-)GGAs and the NGA functional MN15-L are tested next. The latter is the local, i.e. no nonlocal Fock exchange, version of the hybrid MN15. However, unlike for its hybrid counterpart the performance is not acceptable with an MAD over all base stacks of 1.43 kcal/mol. The cause is a heavy overbinding (MD < 0, Table S11) of the intra-strand stacking energies (MAPD 49%). Conversely, the interstrand pairs are slightly underbound. MN15-L shows together with the SCAN functional the highest MAPD for the inter-strand pairs among all (meta-)GGA/NGA functionals (value over 9%). This indicates some sort of imbalanced description. The D3(0) correction is not able to correct such strong overbinding by construction and the actual dispersion contributions to the energy are essentially zero. Thus, MN15-L cannot be recommended for evaluation of nucleobase stacking either with or without dispersion correction. For BLYP-D3(BJ) the excellent performance of LYP combined with D3(BJ)20 carries over from the hybrid B3LYP and the GGA version shows an MAD of just 0.20 kcal/mol. Five GGA DFAs yield an MAPD for the intra-strand stacking of around or below 5% [B97-D3(BJ), TPSS-NL, VV10, B97M-V and SCAN-D3(BJ)], all of which also show an excellent MAD close or below 0.2 kcal/mol. Among those B97M-V performs best with an MAD of only 0.10 kcal/mol, closely followed by our new revTPSS-NL and SCAN-D3(BJ) (MADs of 0.13 and 0.14, respectively). Again, this constitutes essentially the same performance as achieved for the double-hybrid methods. The computationally cheaper B97M-D3(BJ)



Page 18 of 43

variant still offers good performance with an MAD of 0.27. Note that the dispersion-uncorrected SCAN shows the worst performance among all (meta-)GGA/NGAs and cannot be recommended for the calculation of nucleobase stacking. Both the intra- and inter-strand stacking is strongly underbound and a dispersion correction is mandatory. A serious drawback of the SCAN functional is a very strong dependency on the radial quadrature grid size;135 see also the note about the strongly increased radial grid in the Computational Details. Commonly, meta-GGAs often show an increased, but less extreme dependency for the angular (spherical) grid.136 The increased demand on the radial grid severely increases the computation time of energies and gradients. As local functional's main case of application is geometry optimization, a strong grid dependency is a notable drawback for application work on large molecules.

Table 2. Statistical evaluation of individual base stacking energies. a

Method

MAD

all base-base stacks MD |MAX|

inter-strand MAPD [%]

intra-strand MAPD [%]

double-hybrids B2PLYP-D3(BJ) B2PLYP-D3M(BJ) B2PLYP-NL B2GPPLYP-D3(BJ) B2GPPLYP-NL DSD-BLYP-D3(BJ) DSD-BLYP-NL DSD-PBEP86-D3(BJ) DSD-PBEP86-NL PTPSS-D3(BJ) DSD-PBEB95-D3(BJ)

0.18 0.12 0.23 0.35 0.17 0.19 0.08 0.11 0.09 0.43 0.18

0.18 0.05 −0.22 0.35 −0.14 −0.10 −0.02 0.07 −0.08 0.43 −0.04

0.39 0.36 0.69 0.61 0.55 0.61 0.38 0.32 0.37 0.85 0.63

1.0 1.1 1.6 1.8 1.2 1.8 0.8 0.9 0.8 2.1 1.8

4.2 1.7 6.9 9.5 4.6 5.1 2.6 1.9 2.4 11.0 4.4

0.16 0.15 1.02 0.34 1.67 0.39 1.81 0.68 4.36 0.49 0.63 0.30 0.40 0.41 0.46 0.15 0.22 0.18 0.35 0.38 0.21 0.28 0.29

0.00 0.01 1.02 0.34 1.67 0.33 1.81 0.68 4.36 0.22 0.63 0.03 0.39 0.37 0.44 −0.10 0.16 −0.08 0.35 0.36 0.14 0.24 0.26

0.53 0.50 1.89 1.14 3.68 1.38 3.30 1.74 8.46 2.02 1.46 0.90 1.27 1.31 1.26 0.46 0.87 0.56 1.05 0.91 0.64 0.86 1.03

1.0 1.0 5.3 1.5 9.8 2.1 9.7 3.1 25.4 4.5 3.1 2.2 2.3 2.3 2.5 1.0 1.3 1.4 1.6 1.5 0.7 1.3 1.9

3.4 2.9 24.5 7.8 49.5 12.8 47.4 15.6 128.8 5.4 11.1 5.2 13.3 13.9 15.8 3.1 5.2 3.8 6.2 11.2 5.4 7.4 8.9

0.20

−0.08

0.51

1.6

6.0

(meta-) hybrids B3LYP-D3(BJ) B3LYP-NL M06-2X M06-2X-D3(0) M05-2X M05-2X-D3(0) M11 M11-D3(BJ) SOGGA11X SOGGA11X-D3(BJ) MN15 MN15-D3(BJ) PBE0-D3M(BJ) PBE0-D3(BJ) PBE0-NL revPBE0-NL LC-VV10 ωB97M-V ωB97M-D3(BJ) PW6B95-D3(BJ) LC-ωPBE-D3(BJ) LC-ωPBE-D3M(BJ) CAM-B3LYP-D3(BJ) (meta)-GGAs/NGA BLYP-D3(BJ)


Page 19 of 43 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


revTPSS-D3(BJ) 0.37 0.33 1.14 1.8 11.8 revTPSS-NL 0.13 −0.08 0.46 1.1 2.9 TPSS-D3(BJ) 0.44 0.34 1.18 2.2 13.7 TPSS-NL 0.16 0.04 0.39 0.4 3.7 B97-D3(BJ) 0.21 −0.10 0.63 0.6 5.4 B97-D3M(BJ) 0.52 0.51 1.39 3.0 16.7 B97M-D3(BJ) 0.27 0.27 0.79 1.6 6.2 B97M-V 0.10 0.04 0.55 0.8 1.4 VV10 0.18 −0.18 0.40 0.9 5.6 PBE-D3(BJ) 0.46 0.37 1.23 2.6 14.5 PBE-D3M(BJ) 0.45 0.40 1.20 2.6 14.5 SCAN 1.65 1.65 3.49 9.9 48.3 SCAN-D3(BJ) 0.14 0.03 0.51 0.8 3.0 BP86-D3(BJ) 0.37 −0.15 1.01 3.3 10.8 MN15-L [-D3(0)]b 1.43 −1.09 3.71 9.6 48.9 a Mean (absolute) deviation (MD, MAD) and maximum absolute error (|MAX|) in kcal/mol and mean absolute percentage deviation (MAPD) in percentage; all values CP-corrected. b D3(0) dispersion contribution is negligible; basis set: def2-QZVP.

Comparison of dispersion correction variants. For some functionals multiple types of dispersion correction were applied, which allows a comparison between the -D3 and –NL approaches, as shown in Figure 3. For B2PLYP and PBE0 we chose the D3M(BJ) variant as it gave better results. However, D3M(BJ) is not consistently better as seen for B97-D3 and LC-wPBE-D3. The largest accuracy gain is found for PBE0-D3, as already reported in the D3M paper.137 More interesting is the D3 versus NL comparison. Few earlier comparisons for benchmarks can be found in the literature.48,138,139 In general we find that the VV10 correction does a better job at capturing the dispersion for our stacking benchmark than the D3 correction. Largest differences are found for the set of meta-GGA functionals and ωB97M-V. Our newly proposed VV10 parameters for revTPSS work also very well and revTPSSNL is even the second best (meta)GGA tested, slightly behind the excellent B97M-V. A major drawback of these comparisons is that the fitting of each respective correction is not always done on equal footing, i.e., either different sizes and types of training sets have been used or a different quality of level of theory for the reference data was applied (details can be found in the references of the dispersion parameters given in the Supporting Information). In other words, some fitting procedures were perhaps more general than other more specialized fits. Proper assessment of the performance of different dispersion corrections could potentially require a re-fit of some of the parametrizations. At least for stacking interactions, the B2PLYP-NL parametrization seems not optimal.



Figure 3. Comparison of the MADs for D3 (black) versus VV10 (blue) dispersion corrections. Low-cost QM methods, semi-empirical methods and MM force fields. We group the low-cost and semiempirical methods together with the force fields for convenience as only a limited set of methods has been tested. The low-cost DFT methods HF-3c, PBEh-3c and B97-3c are composite methods defining DFA, basis sets and correction potentials. The three correction "3c" suffix indicates some form of a treatment for dispersion and basis set effects. Unfortunately, the "3c" corrections differ among the methods. They all have in common the D3 dispersion correction and correction potentials for basis sets effects. HF-3c includes both the geometrical counterpoise correction140 (gCP) and a short-range basis set correction (SRB), while PBEh-3c only includes gCP (in a modified variant) and B97-3c only includes the SRB correction. The basis sets increase from a quasi-minimal one for HF-3c, to a polarized double-ζ for PBEh-3c and to a polarized triple-ζ for B97-3c. A recent review over these composite or simplified DFT methods was given by Caldeweyher and Brandenburg.141 Although a substantial amount of empiricism is included in the construction of these simplified DFT methods, they are far from being conventional semi-empirical methods, which parametrize Hamiltonian matrix elements directly. Because of the way they are constructed no Boys & Bernardi counterpoise correction was applied. The gCP correction takes care of BSSE for HF-3c and PBEh-3c, while B97-3c employs a special triple-ζ basis and absorbs part of the BSSE in the D3 correction and functional re-parametrization. For biomolecular systems only few semiempirical methods are accurate enough to consider.142 Table 3 lists three which we believe are in general applicable for biomolecules, building upon previous experience in our lab, namely PM6-D3H+, DFTB3-D3 and GFN-xTB.143–145


Page 20 of 43

Page 21 of 43 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Best performing low-cost method is B97-3c with an MAD of 0.31 kcal/mol over all base stacks and excellent intra-strand stacking accuracy with an MAPD of 6.9%. HF-3c and PBEh-3c show a similar performance, with HF3c being somewhat better at the intra-strand stacking. Here HF-3c has the advantage of avoiding any doublecounting of medium-range correlation effects due to the absence of a correlation functional. However, for electrostatics, which plays a non-negligible role for the stacking of nucleobases, HF is less suitable compared to DFT due to its tendency to overpolarize the nucleobases. The error sources of HF-3c and PBEh-3c are fundamentally different and their similar results are thus mostly coincidental. Table 3. Statistical evaluation of individual pair-stacking energies for low-cost QM methods and selected force-field versions. a Method

MAD


inter-strand intra-strand MAPD MAPD

low-cost methods PM6-D3H+ DFTB3-D3 HF-3c PBEh-3c B97-3c GFN-xTB

0.38 0.53 0.42 0.40 0.31 0.77

−0.20 0.18 0.25 0.18 0.09 0.72

1.80 2.76 1.33 1.35 0.96 2.05

2.8 2.6 1.6 1.7 1.0 3.0

7.5 11.7 8.5 10.2 6.9 24.8

force-field

Details.

AMBER(library) b 0.64 AMBER(HF) b 0.34 AMBER(PBE0stacks)b 0.49 AMBER(PBE0AT,GC)b 0.64 CHARMM(library) 1.18 CHARMM(HF) 0.44 GROMOS(library) 1.26 GROMOS(HF) 1.40 avalues for MAPD in percentage, otherwise in

−0.30 1.78 4.6 0.09 1.11 1.2 −0.34 1.24 3.4 −0.14 1.81 3.4 −0.57 4.20 8.6 −0.07 1.16 1.9 0.86 5.18 5.9 1.39 3.67 8.9 kcal/mol; bbrackets indicate the set of

16.2 10.5 12.4 13.6 29.2 12.1 48.9 46.5 point-charges – see Computational

PM6-D3H+ performs well for the calculation of stacking energies and the description of the inter- and intramolecular stacks seems balanced. Compared to the other low-cost methods PM6-D3H+ shows over- instead of underbinding. The dispersion-corrected tight-binding method DFTB3-D3(BJ) yields better energies than the conceptually similar GFN-xTB method. GFN-xTB shows the strongest underbinding of all low-cost methods, which is especially visible for the intra-strand stacking where the MAPD reaches almost 25%. Note that other parametrization and correction potentials for PM6 and DFTB have been proposed in the literature.146,147 Among empirical force-field methods, AMBER is the most successful force field for nuclei acids. Its performance is of major interest for the classical molecular dynamics simulation community, as it is the most widely used force field for extended explicit-solvent simulations of nucleic acids.148 In addition, AMBER has a straightforward procedure to derive the point-charge distributions.82 Thus, four different sets of charges for the AMBER force field are used in our work, as described in the Computational Details section. For the sake of



completeness, we also tested the CHARMM36 and GROMOS 53a6 force fields. Obviously, the presented benchmark data allow the readers to perform a straightforward testing of any other force field or method. The AMBER force field with library charges [AMBER(library)], corresponding most closely to those actually used in MD simulations (except of the charge on the capping hydrogen atom as mentioned in the Computational Details), performs at a similar level of accuracy as mediocre DFAs with an overall MAD of 0.64 kcal/mol. A slight tendency for overbinding for intra-strand stacking is observed. This is somewhat counterbalanced by an underbinding tendency of the intra-strand stacks. Note that in the library-charge parametrization of the AMBER force field, the charges are derived by fitting to reproduce the molecular ESP derived at the HF level around the nucleosides (and a dimethylphosphate), so that the R(N9)/Y(N1) nitrogens of glycosidic bonds are bound to ribose C1' carbons.130 Refitting the charges at the same HF/6-31G(d) level of theory just for the nucleobases [AMBER(HF)] thus follows the AMBER charge-fitting philosophy. It mimics the AMBER-like charges while it simultaneously allows the charge to redistribute, considering the hydrogen capping of the R(N9-H9)/Y(N1-H1) groups. This refitting leads to a significant improvement of gas phase stacking energies and better agreement with the QM reference data. The MAD drops to 0.34 kcal/mol and the MD shows only marginal underbinding (0.09 kcal/mol). In other words, the AMBER(HF) results show no significant systematic under- or overbinding of the stacking energies. That electrostatics (including polarization effects) plays a substantial role becomes also visible when performing the force-field calculations with the AMBER(PBE0stacks) charges, which are expected to i) yield better gas phase MM electrostatics allowing a fairer comparison against QM data and ii) include pair-wise polarization effects implicitly as the RESP fit was repeated for all respective inter- and intra-strand stacks. Such charge distributions should shift the force-field electrostatics as close as possible to the QM calculations. Indeed, the overbinding of the intra-strand stacks is reduced to about half of its value compared to AMBER(library) (MD is shifted from −0.77 to −0.44 kcal/mol). However, now also the inter-strand stacks show overbinding and the error compensation effects seen with the library charges are no longer present. Consequently, the MD evaluated over all base stacks is essentially the same for both sets of charges. The MAD is slightly reduced to 0.49 kcal/mol. The second set of PBE0 charges, AMBER(PBE0AT,GC) that uses charges determined on the indicated H-bonded basepairs (see Computational Details), leads to similar MAPDs for both the inter- and intra-strand stacking as PBE0stacks, but has little impact on the overall MAD. However, the systematic tendency to overbind is reduced as visible from the smaller MDs, particularly for the intra-strand stacking (Table S12, Supporting Information). Nevertheless, as the |MAX| value does not improve compared to the “library” charges, and the MAD is similar, also this set of


Page 22 of 43

Page 23 of 43 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


charges is not fully satisfying. Since the non-bonded potential in AMBER is heavily parametrized and an 'effective' potential that relies on error compensation between its components, any imbalance can have strong effects on its accuracy. Here, it seems that the AMBER(library) version does not work well for base stacking calculations on truncated model systems using just the nucleobases. The AMBER(HF) charge set is definitely more balanced for base stacking calculations. Overall, our results support the view that while the AMBER force field provides an impressively good description of the base stacking, taking into account its simplicity, it is definitely far from being quantitatively accurate.149 Even the remarkably well-performing AMBER(HF) version yields a quite large |MAX| of over 1 kcal/mol. The simplicity of the MM terms may preclude achieving further significant improvements, especially when taking into consideration that the force field would have to be then still balanced with all the other energy terms. The CHARMM(library) force field shows a similar overbinding tendency as the AMBER force field but with an overall increased error leading to an MAD of 1.18 kcal/mol. GROMOS(library) parameter set, on the other hand, shows underbinding behavior and a slightly worse MAD of 1.26 kcal/mol than CHARMM. For the present benchmark set, neither CHARMM nor GROMOS seem competitive to the AMBER force field. A comparison of the stacking energies obtained with AMBER(HF), CHARMM(HF), and GROMOS(HF) parameters suggests that the main source of error in the CHARMM(library) stacking calculations is the electrostatics, while the GROMOS FF does not seem to improve with the HF charges at all, tentatively indicating a larger imbalance of the non-bonded terms for nucleic acids. The present results confirm that the good performance of the AMBER force field to describe stacking of conformationally rigid nucleobases comes from a suitable charge-fitting model aiming to correctly reproduce the ESP surrounding the molecule.149 However, it is important to mention that various variants of CHARMM and GROMOS force fields are widely used for simulations of protein dynamics2 and that the CHARMM force field provides also a very good overall performance of B-DNA dynamics.85,150 In addition, the present results cannot be directly extrapolated to predict the performance for explicit-solvent simulations of nucleic acids. Table 4. Statistical evaluation of individual pair-stacking energies for wavefunction methods.a Method

MAD


inter-strand intra-strand MAPD MAPD

supermolecular methods DLPNO-CCSD(T)/CBS(def2T,def2Q) DLPNO-MP2-F12D/VQZ-F12-CP DLPNO-SCS-MP2-F12D/VQZ-F12-CP MP2/CBS(aVTZ,aVQZ)

0.11 0.98 0.53 1.03

0.06 −0.98 0.53 −1.03


0.25 2.85 1.16 2.94

0.5 6.7 2.5 7.0

3.5 33.6 15.5 34.9


SCS-MP2/CBS(aVTZ,aVQZ) SCS(MI)-MP2/CBS(VDZ,VTZ)-CP MP2-F12/CBS(VTZ-F12,VQZ-F12) MP2-F12/CBS(VTZ-F12,VQZ-F12)-CP MP2.5/VTZ-CP MP2.5/VTZ MP2/CBS(aVTZ,aVQZ)+ΔHL(MP2.5/VTZ) MP2/CBS(aVTZ,aVQZ)+ ΔHL(FNO-CCSD(T)/aVDZ)

0.51 0.30 1.12 0.97 1.09 0.75 0.24 0.39

0.51 0.27 −1.12 −0.97 1.09 −0.75 −0.03 0.39

1.15 0.68 3.12 2.84 2.11 1.79 0.90 0.78

Page 24 of 43

2.3 1.7 7.5 6.6 6.3 4.6 0.7 2.2

15.0 7.9 37.8 33.1 33.8 22.7 9.0 12.1

0.08 0.00 0.42 0.4 0.66 −0.65 1.62 4.1 0.30 −0.28 0.93 2.2 0.48 −0.48 1.18 3.1 0.24 −0.21 0.79 1.8 0.61 0.46 2.36 1.8 avalues for MAPD in percentage, otherwise in kcal/mol; b junVTZ c junVDZ.

3.0 20.7 8.1 14.1 6.3 17.4

intermolecular methods SNS-MP2 SAPT2+3b SAPT2+3δ(MP2)b SAPT2+(3)b SAPT2+(3)δ(MP2)b SAPT0c

Wavefunction methods. The statistics for a variety of wavefunction-based methods are listed in Table 4. Due to the size of the benchmark set only a limited number of wavefunction-based methods was tested. We have included the recently proposed spin-network-scaled MP2 (SNS-MP2) protocol.67 Note that the protocol is a significantly more elaborate approach than just a simple rescaling of the spin-components of MP2. SNS-MP2 is inherently CP corrected and also employs SAPT-like terms. As such it is meant to compute interaction energies directly and not through absolute ground state energies. A major drawback is a missing gradient for geometry optimizations. The performance of SNS-MP2 is exceptionally good and essentially on par with the DLPNOCCSD(T)/CBS results. A particular strength of SNS-MP2 is that all MDs are close to zero, indicating no systematic over- or underbinding. SNS-MP2 turns out to be a highly accurate method viable for reference calculations for nucleobase stacking. The |MAX| error of DLPNO-CCSD(T)/CBS, however, is still smaller than for SNS-MP2 (0.25 versus 0.42 kcal/mol, respectively). Such robustness is an important aspect for a potential reference method and the more rigorous DLPNO-CCSD(T) has a principle advantage here. MP2 and MP2-F12 reveal their well-known tendency of overestimating stacking interactions.49,151 The MP2 basis set limit can be deduced from the comparison of MP2/CBS, MP2-F12/CBS and MP2-F12/CBS-CP data. The difference between the MP2-F12/CBS and MP2-F12/CBS-CP data is minimal. It is below 0.05 kcal/mol for the MAD calculated for all base stacks and around 0.2 kcal/mol for the intra-strand interaction energies (cf. Supporting Information Table S13). Anticipating a spurious BSIE even at our CBS extrapolation that usually causes underbinding,121,122 the "true" MP2 basis set limit should lie somewhere between those values. This results into a rather small error range for the MP2-F12 values used in constructing our CBS(T)-F12-CP reference method. The conventional MP2/CBS results are somewhat closer to the MP2-F12 results, indicating that the commonly employed aVTZ to aVQZ MP2 extrapolation is already near convergence for our nucleobase stacks. The basis set ACS Paragon Plus Environment

Page 25 of 43 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


convergence directly translates to the SCS-MP2 results, and the conventional SCS-MP2/CBS results are decent, but not competitive with the best density functional methods. Especially the MAPD of the intra-strand stacking of 15% is rather poor. The parametrized SCS variant SCS(MI)-MP2 surpasses standard SCS-MP2 for the stacking energies and yields the overall third best results for the wavefunction methods. We note a slight systematic underbinding as already reported previously.9 Nonetheless those results are good, and the energy scan for the vertical separation based on SCS(MI)-MP2 is deemed reliable for our purposes. It does also have the advantage of a gradient implementation. The DLPNO-MP2-F12/D results (no extrapolation, but CP correction) compare very well with the respective canonical MP2-F12 results, which is a promising outlook to reach the basis set limit for large molecules. It can be easily envisioned to construct a focal-point energy function of DLPNO-MP2/CBS with a ΔHL(DLPNO-CCSD(T)) correction. However, potential outliers from problematic systems for local methods, e.g. strongly delocalized or low band-gap systems, must be carefully studied. As discussed above, the DLPNO-CCSD(T)/CBS results show an excellent agreement with the canonical CBS(T)-F12-CP method and usage of the former method as a reference method as extensively done in the GMTKN55 study15 appears justified. Not suitable as a reference method are the energy functions relying on a lower level of theory ΔHL correction using MP2.5/VTZ152 or FNO-CCSD(T)/aVDZ with a standard MP2/CBS(aVTZ,aVQZ) extrapolation. The MP2.5 based variant yields a decent MAD of 0.24 kcal/mol but ultimately revealing the limitations of using MP2.5 in the ΔHL correction. Using the FNO-CCSD(T)/aVDZ ΔHL correction leads to an even worse MAD of 0.39 kcal/mol. This shows once again the insufficient quality of a mere double-ζ CCSD(T) high-level correction even with added diffuse functions. Both methods are not competitive with the best dispersion-corrected DFAs which shows that for benchmarking modern DFAs care must be taken that the reference method is accurate enough. We also include high-order SAPT calculations as a non-supermolecular approach that is inherently free of BSSE effects. The statistics in Table 4 clearly demonstrate the importance of the supermolecular δ(MP2) correction. It significantly reduces the MADs of both the SAPT2+3/junVTZ and SAPT+(3)/junVTZ. The δ(MP2) term introduces additional coupling between induction and dispersion. One finds an almost equal improvement for both the intra- and inter-strand MAPDs. Due to the significant improvements the inclusion of the supermolecular δ(MP2) correction seems mandatory for base stacking in a triple-ζ basis set. SAPT in general can benefit from error compensation effects of its individual components and it is important to include interaction terms in a balanced way. We see such error compensation effects comparing SAPT2+(3) and SAPT2+3 data, where the latter as the



higher-order SAPT treatment shows a reduced accuracy. Favorable error compensation effects lead to the costeffective SAPT0 treatment, which works best in a slightly augmented double-ζ basis (junVDZ).65 The SAPT0 method shows overall similar results for the total interaction energy as the SAPT2+3 level at only a fraction of its cost. However, SAPT0 is systematically underbinding the stacked nucleobase dimers while the higher-order SAPT2 methods are overbinding. To the best of our knowledge an appropriate δ(MP2) correction has not been proposed but based on a comparison with the higher-order SAPT results it may be beneficial.

Table 5. Statistical evaluation of the sum of pair-stacking energies ΔEstack for selected methods. a method MAD DLPNO-CCSD(T)/CBS 0.23 DSD-BLYP-D3(BJ) 0.41 DSD-BLYP-NL 0.12 B2PLYP-D3M(BJ) 0.21 B3LYP-D3(BJ) 0.15 B3LYP-NL 0.13 ωB97M-V 0.32 M06-2X 4.07 M06-2X-D3(0) 1.35 revTPSS-NL 0.34 B97M-V 0.17 PBEh-3c 0.74 B97-3c 0.38 SAPT2+3δ(MP2) 1.10 SAPT0 1.83 SNS-MP2 0.21 AMBER(library) 1.22 AMBER(HF) 0.75 AMBER(PBE0stacks) 1.66 AMBER(PBE0AT,GC) 0.86 CHARMM(library) 2.57 areference: CBS(T)-F12-CP; values in kcal/mol;

MD 0.23 −0.41 −0.09 0.19 0.02 −0.03 −0.32 4.07 1.35 −0.31 0.17 0.74 0.35 −1.10 1.83 −0.01 −1.19 0.35 −1.38 −0.54 -2.29

A short summary of the results from the previous subsections is collected in Figure 4, presenting the MAD/MD over all bases for a selected subset of methods containing the best methods of Tables 2-4 and key results, such as the AMBER(HF) re-fit and the newly fitted DFA revTPSS-NL.


Page 26 of 43

Page 27 of 43 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Figure 4. Graphical summary of few key results from Tables 2-4. Benchmarking of ΔEstack for the full base-pair step stacks. Table 5 reports the statistics for the sum of all four pair-stacking energies (see Computational Details) for each base-pair step stack, i.e., the ΔEstack data. As noted above, the summation can lead to both error cancellation and error enhancement effects. However, for the sake of completeness we present also the statistics for the ΔEstack values for a subset of methods. The trends for the DFT and wavefunction methods largely agree with the statistics over the individual base stacks and will not be discussed in greater detail. Still, there are some deviations, e.g., SNS-MP2 had the lowest MAD of all methods – together with DSD-BLYP-NL – for the individual stacks benchmarking, but not for the MADs in Table 5, where other methods like B3LYP-NL are performing better. One can expect that methods with a MD close to zero in the individual stack benchmarking show less error compensation or enhancement. As such, one can especially recommend DSD-BYLYP-NL, B3LYP-NL and M97M-V from the DFT section, SNS-MP2 and DLPNO-CCSD(T) from the wavefunction section, and B97-3c from the low-cost method section for the calculations of ΔEstack as they show very low MADs and balanced MDs in the individual stacks benchmarking. The AMBER force-field calculations reveal an interesting behavior. While the individual stack analysis shows that the PBE0stacks charge re-fit improves the overall statistics, for the sum ΔEstack the opposite is found. The library charges perform better than the PBE0stacks ones. This is a good example of error compensation effects. The AMBER(library) charge set leads to an overbinding of the intra-strand stacks, and an underbinding of the interstrand stacks (see Supporting Information Table S12). For the PBE0stacks re-fit both types of stacked pairs show overbinding. Thus, the PBE0stacks charge set does not benefit from the error compensation. However, the MDs of



both charge sets show overbinding by about 1.3 kcal/mol on average for ΔEstack, revealing a systematic error. The opposite is observed for the PBE0AT,GC set of charges. In the pair-wise benchmark it showed a similar MAD over all base stacks as the library AMBER charges. For ΔEstack, however, AMBER(PBE0AT,GC) error compensation effects are favorable, so that the MAD is significantly reduced from 1.22 for AMBER(library) to 0.86 kcal/mol, and the overbinding tendency is reduced by more than half of its value to −0.54 kcal/mol. However, the best results are again obtained with the nucleobase-optimized AMBER(HF) charges, for ΔEstack showing even a slight underbinding. Note, that this does not necessarily mean molecular dynamics simulations would benefit from these charges, which anyway would have to be merged with the backbone charge set. It just means the HF re-fit suits better for calculations of isolated stacks while the AMBER(library) charges derived for nucleosides are likely better suited for nucleic acids simulations. However, when stacking geometries are subtracted from simulation trajectories to monitor base stacking energies, one might consider using the HF charges rather than the library charges. Nevertheless, based on all the data it seems that the effect of charge distributions on MM stacking calculations would deserve a broader analysis for a much larger set of geometries which we plan to do in future. Four-body effects in stacked base-pair steps. London-type dispersion interactions are known to be nonadditive153 and dispersion-corrections for DFT should ideally capture many-body dispersion. The herein applied D3 correction is pair-additive but can optionally include a three-body term [D3(ABC)].154 The nonlocal dispersion correction (DFT-NL), which is based in the VV10 correlation kernel is a density-dependent dispersion correction and captures many-body dispersion if applied self-consistently (as done herein). For wavefunction-based electron correlation methods non-additive effects arise earliest in the hierarchy of excitation classes where doubles excitations are coupled with other excitations. This takes place in perturbation theory earliest at the MP3 and for coupled-cluster at the CCD level. MP2, on the other hand, shows similar non-additivity as HF as the pair-correlation energies are simple sums over (uncoupled) double excitations and do not directly contribute non-additive correlation energy. Both VV10 (if applied in a post-SCF fashion) and MP2, however, contain indirect many-body effects through SCF variations in the HF orbitals or DFA densities. Previous investigations of the four-body contribution to the stacking interaction of stacked base pairs by one of us found the largest values for the CpC step (2.2 kcal/mol) and values of around 1 kcal/mol also for CpG and GpC steps (cf. Table 6).27,29 The large repulsive (anti-cooperative) non-additivity of the CpC step is related to the presence of two almost parallel intra-strand GG and CC homo-stacks, as both G and C bases are very polar.


Page 28 of 43

Page 29 of 43 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Results from the DLPNO-CCSD(T)/CBS calculations show generally larger four-body term values than the older data computed at the MP2/aVDZ level.29 While the increased correlation energy non-additivity revealed in the present work is substantial, the earlier trends were confirmed, since also at the DLPNO-CCSD(T)/CBS level the strongest four-body contribution is found for the CpC (CC:GG) base-pair step stack. Note that with our new benchmark calculation the four-body term reaches in absolute value almost 30% of the total stacking energy of the CpC step, which is very significant. The benchmark of the DFT-D3(BJ) methods in Table 2 has been performed without the three-body dispersion term D3(ABC). The main motivation behind this is that most D3(BJ) parameters have been initially parametrized without the D3(ABC). The D3(ABC) contribution is destabilizing for all nucleobase stacks, which is in line with previous benchmark studies on large systems.155 While Risthaus and Grimme recommend the inclusion of D3(ABC), they do note that all their DFT-D3 results are slightly overbinding.155 They further use a "half CP" correction as their Pauling point. As CP works in the same direction as D3(ABC) the choice of a scaled CP correction makes sense, but the interplay between spurious BSSE in quadruple-ζ DFT calculations and D3(ABC) contributions is ambiguous. As we employ a full CP correction the exclusion of D3(ABC) for pairwise stacking energies is a sensible choice. For our systems all this indicates that D3(ABC) is only beneficial for those few D3(BJ) functionals that show overbinding behavior like, e.g., B97-D3(BJ) and BP86-D3(BJ). The D3(ABC) contribution to the intra- and inter-strand stacking is small, on average 0.14 kcal/mol, ranging from practically zero to 0.3 kcal/mol. This brings the D3(ABC) contribution to about a similar magnitude as the CP-correction in the DFTD3(BJ) calculations (average BSSE for B3LYP is 0.17 kcal/mol), and in the same energy direction. For the ΔEABCD term D3(ABC) contributes 0.3 to 0.4 kcal/mol. While our data certainly proves a weak, but notable influence of the D3(ABC) energy, we are reluctant to formulate a final statement about the importance of the inclusions of D3(ABC) contributions in DFT-D3 for base-pair step stacking calculations. We find an overall qualitative agreement of the four-body terms of DLPNO-CCSD(T) and the tested DFAs. For the most part, the relative magnitudes seem to be reliable with DFT as the MADs are between 0.70 and 0.93 kcal/mol. A somewhat systematic shift towards smaller four-body effects is seen for the DFAs. If the four-body term is small, as it is for the ApA, ApT and TpA steps, this error can lead to opposite signs of the non-additivity. The origin of this systematically smaller non-additivity for the tested DFAs is not well-understood and further studies are required.



Finally a short note on many-body BSSE156 effects should be made. This point has recently been systematically investigated by Richard et al.157 We employed in the calculation the conventional CP correction, i.e. the base-pair electronic energies E12 and E34 are evaluated in the supermolecular basis sets of the whole base-pair step. For ΔEstack calculations, however, the supermolecular basis sets include only the nucleobase dimers, not the full tetramer. This means we might miss many-body BSSE effects in the calculation of the four-body term ΔEABCD. To estimate an error bound for our values we computed the BSSE for the total interaction energy of the base-pair stacks as a 'worst case' scenario. We compared the conventional CP approach, which is an approximation for many-bodies called site-site functional counterpoise (SSFC) correction,158 with the correct BSSE expansion based on the Valiron and Mayer functional correction (VMFC).159 The VMFC expansion requires 65 single-point calculations for our four-body system and constitutes a sizable computational effort. We tested the VMFC expansion for the total interaction energy of the GpC base-pair step stack at the BLYP-D3(BJ)/def2-QZVP level of theory. The VMFC interaction (not stacking!) energy of the four bases is -79.35 kcal/mol. SSFC differs by only 0.05 kcal/mol to it while the total amount of (VMFC) BSSE is 1.07 kcal/mol. As we are limited to stacking energies instead of interaction energies 0.05 kcal/mol is an upper bound for the many-body BSSE error of the four-body term. It shows that the SSFC approach is sufficiently accurate and our values for the four-body terms thus should be reliable. The BSSE is typically very similar across GGAs and hybrid DFAs140 so that this result should be reasonably transferrable to other DFAs. Table 6. Four-body non-additivity ΔE ABCD of different methods. DLPNOCCSD(T)

B3LYPB3LYPB3LYP- DLPNO- ωB97M- B3LYPB3LYPB3LYPa D3(ABC)a MP2b NL D3(BJ) D3(ABC,BJ) CCSD(T) V NL D3(ABC,BJ) ΔE4stack ΔEABCD ApA −12.95 −14.04 −13.81 −13.72 −13.37 0.21 −0.32 −0.57 −0.29 0.35 −0.35 0.0 ApC −11.88 -13.00 −12.85 −12.86 −12.51 1.36 0.67 0.43 0.66 0.35 0.50 0.8 ApG −12.54 −13.68 −13.71 −13.47 −13.07 1.14 0.53 0.28 0.56 0.39 0.47 0.8 ApT −11.00 -12.05 −12.02 −12.05 −11.73 0.62 -0.05 −0.28 −0.03 0.32 −0.14 0.0 CpC −10.53 -11.76 −11.68 −11.38 −10.98 2.97 2.30 2.07 2.33 0.41 2.14 2.2 CpG −16.23 -17.80 −17.85 −17.02 −16.59 1.88 0.99 0.78 1.06 0.44 0.94 1.1 GpC −13.94 −15.38 −14.96 −14.97 −14.58 1.83 1.11 0.92 1.14 0.39 1.00 1.2 TpA −12.92 −14.08 −14.07 −13.70 −13.33 0.66 -0.08 −0.30 −0.06 0.36 −0.14 0.2 TpC −11.58 −12.96 −12.65 −12.51 −12.15 1.32 0.59 0.30 0.56 0.36 0.39 0.7 TpG −14.17 −15.50 −15.54 −14.99 −14.59 1.56 0.80 0.60 0.85 0.40 0.75 0.9 MAD 0.0 1.25 1.14 0.89 0.52 0.0 0.70 0.93 0.68 0.52 0.56 Values in kcal/mol; DFT calculations are CP corrected. a Contribution of the three-body D3 term.b MP2/aVDZ taken from ref. 29 ωB97M-V

Interaction energy decomposition of base-pair stacking interaction energies. The interaction energy of stacked nucleobases has been studied in detail over the past decades and it is generally accepted that the stabilization is dominated by London dispersion while the Twist dependence of the stacking energy is to a high degree affected by electrostatics.5,24 Nucleobase pair stacks show an increased component of electrostatics over, e.g., benzene or anthracene stacks due to their heteroaromatic nature, leading to a dipolar shape of the ACS Paragon Plus Environment

Page 30 of 43

Page 31 of 43 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


molecular electrostatic potential. For the stacking of base pairs an energy decomposition was recently done by Parker et al. who investigated the dependence of the interaction energy components on the base pair stacking parameters Rise and Twist, and partially also Slide.37 They found that variations in the Twist parameter cause most changes in the net dispersion (exchange + dispersion) for low Rise values (vertical compression of the basepair step), while for high values of Rise changes in electrostatics dominate. Low Rise values are associated with attractive electrostatics and significant charge-penetration effects. The charge-penetration effects are essentially absent in commonly applied (point-charge) biomolecular force fields and their modelling is still an active area of research.13,160 Other works, however, emphasized that when making energy scans for varying helical or basepair parameters, the vertical separation (Rise) between the base pairs should be optimized for each value of the parameters, due to high energy gradients associated with any vertical compression or stretching of the extended stacked systems.9,24 Further, we reiterate that one has to be careful and not over-interpret the impact of QM energy decompositions for parametrizations and evaluations of force fields, as the individual terms in the decompositions show opposite signs with exponential behaviour upon changing the inter-monomer separations. Many terms, despite looking at first sight dramatically large and variable, may mutually compensate each other. Then it is possible to include their interplay reasonably well into even simple MM energy functions.9,12,28 For computational efficiency Parker et al. applied the SAPT0 level of theory, which, while quite accurate for the total interaction energy due to error cancellation effects, is less accurate for the individual decomposition terms. Herein we try to get accurate SAPT-based decomposition values by applying SAPT2+3δ(MP2). As can be seen from Table 4 the inclusion of the supermolecular δ(MP2) correction is greatly beneficial for the accuracy of the total interaction energy.



Figure 5. Percentage-wise attraction energy decomposition of the ΔEstack. Figure 5 displays the percentage-wise decomposition of the attractive parts of the stacking interaction energy (ΔEstack) for all ten base-pair stacks (full numerical data is in Table S2, Supporting Information). While the dispersion energy is heavily dominating by 70-85%, the electrostatic component is significant and in the range of 10-20%. The ApT step is dominated by the strongest dispersion, while the CpG step shows strongest electrostatics and induction. However, the induction and electrostatic terms do not correlate. Both the ApT and CpC step show similarly weak electrostatics, but the induction in the CpC step is about twice as strong as in ApT (3.6% vs 7.8%, respectively). For this particular case it is known27 that the C and G homodimer intra-strand stacks are electrostatically unfavorable (see also the intra-strand EL terms in Figure 6), reducing the overall electrostatic stabilization. This is also related to the large repulsive (anti-cooperative) four-body term in this step discussed above. In general, because of the different strengths and signs of the inter- and intra-strand stacking energy components (vide infra), the simple data of the sum of attractive components as done in Figure 5 should not be over-interpreted. To understand the nature of ΔEstack in further detail we also looked at the nucleobase pair-wise SAPT decomposition displayed in Figure 6 (numerical data is in Table S3 and S4, Supporting Information). One immediately sees the large impact the DISP component has on the stacking of the intra-strand stacks. However, it is accompanied by a similarly large EX component indicating overlapping charge densities. For the inter-strand stacks DISP and EX are both dramatically smaller in magnitude. The smaller magnitude of DISP causes missing dispersion in DFA to be somewhat less impactful. For the GpC, TpG and TpA steps, the intra-strand DISP component is somewhat smaller than for the other base-pair steps. However, this is accompanied by an increase


Page 32 of 43

Page 33 of 43 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


in the respective inter-strand DISP. The same is found for the EX component that follows the same trend as DISP. As the EL component is not always attractive there is no unambiguous way to make a similar percentage-wise decomposition as for ΔEstack. However, from a visual inspection it seems that IND and EL have a somewhat higher contribution to the interaction energy of the inter-strand stacks, at least for a subset of the ten base-pair stacks like GpC. There is no simple correlation between the IND and EL energy components.

Figure 6. Pair-wise SAPT2+3δ(MP2) decomposition into electrostatics (EL), induction (IND), exchange (EX) and dispersion (DISP) components.

Conclusions Extensive benchmarking of ten idealized B-DNA base-pair step stacks with DFT, correlation methods, semiempirics and force fields was presented. The inter- and intra-stand stacking energies were separately statistically evaluated and many-body stacking effects were studied. High-level SAPT-based energy decomposition provided insights into the governing intermolecular interaction components of base-pair stacking. We obtained a new set of highly-accurate reference base-stacking energies based on a focal-point energy function comprising triple- to quadruple-ζ extrapolated and counterpoise-corrected MP2-F12 data improved by a CCSD(T)/jun-cc-pVTZ high-level correction, which is at the edge of currently feasible calculations with conventional wavefunction methods. Conservatively, we estimate an error of 1-2% of the interaction energies.



The majority of high-quality DFAs show a very good performance for base stacking. Thus, we have a wide spectrum of DFAs available for reliable QM calculations of base stacking. Most Minnesota functionals like M062X, M11 or SOGGA11X require a dispersion-correction to yield acceptable results, with MN15 being the exception. However, also for MN15 the D3(BJ) correction is beneficial. The MN15 performance is not transferred to its local counterpart MN15-L that is heavily overbinding and thus cannot be improved by a dispersion correction. Based on a SAPT2+3δ(MP2) decomposition of the stacking interaction energy we reaffirm the general notion that base stacking is dominated by London-type dispersion interaction (70-85%) and a significant electrostatic component (10-20%) and provide accurate energy decomposition data. The ApT step shows the strongest dispersion interactions, while the CpG step has the strongest electrostatics and induction contributions. Strong four-body effects of over 1 and even 2 kcal/mol were found in many base-pair steps. For the CpC step they reach almost 30% of the total stacking energy. Calculations of the four-body effects using three different DFAs show an overall good qualitative agreement with the DLPNO-CCSD(T) reference, albeit with a systematic shift towards smaller values by 0.5 to 1 kcal/mol, which, for weak four-body effects results in a change of sign. While the inclusion of the D3(ABC) correction, at least for B3LYP, slightly improves the overall statistics for the full ΔE4stack base-pair step stacking energies, it is not necessarily improving the four-body term in all the individual base-pair steps. Thus, due to the generally destabilizing nature of D3(ABC), its inclusion is beneficial only in stacking calculations using slightly overbinding DFAs. A strong charge-dependence of the AMBER force field results was found. The charges taken directly from the library with neutralization of the capping hydrogens are less suitable for stacking calculation than RESP HF/631G(d) charges (HF charges in short) derived for isolated nucleobases. The more physically accurate PBE0 charges perform better than the library charges, but are inferior to the HF charges, at least for the present set of geometries. The CHARMM FF was also improved using the same HF charges as for AMBER(HF). For calculations of isolated nucleobase stacking the three tested force fields rank as AMBER > CHARMM > GROMOS in accuracy.

Below is a short summary of the best performing methods for each computational level: double-hybrids: Outstanding performance of DSD-BLYP-NL as the overall best DFA and the best performing method overall (together with SNS-MP2). Other excellent choices are B2PLYP-D3M(BJ), DSD-BYLP-D3(BJ) and DSD-PBEP86-NL.


Page 34 of 43

Page 35 of 43 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


global-hybrids: B3LYP with either -NL or -D3(BJ) correction and revPBE0-NL perform best. B3LYP-NL has a slight advantage. range-separated hybrids: Clear winner is ωB97M-V with both LC-VV10 and LC-ωPBE-D3(BJ) being the second choice. meta-GGAs: B97M-V performs best while revTPSS-NL and SCAN-D3(BJ) follow closely. GGAs/NGAs: The three functionals VV10, BLYP-D3(BJ) and B97-D3(BJ) can all be recommended. wavefunction methods: The best performing supermolecular approach is DLPNO-CCSD(T)/CBS, for which we reaffirm that it is a viable reference method on its own. SNS-MP2 is an excellent choice and tied as best method overall with DSD-BLYP-NL, albeit being limited to interaction energies only. From the SAPT hierarchy the SAPT2+(3)δ(MP2) level is recommended, where the δ(MP2) correction is found to be crucial for accurate results. DLPNO-MP2-F12/D is established as an accurate alternative to canonical MP2-F12. low-cost methods: The most robust simplified DFT method is B97-3c and for semi-empirical methods DFTB3D3(BJ) performs the best. In summary, for accurate single point energies of nucleobase stacking we recommend the double-hybrid DFA DSD-BLYP-NL with quadruple-ζ basis set and CP correction, being more cost-effective than the DLPNOCCSD(T)/CBS and SNS-MP2 methods. The functionals B3LYP-D3/NL, revPBE0-NL, B97M-V, revTPSS-NL and B97D3(BJ) are all excellent choices if a trade off in accuracy for computational efficiency is acceptable.

Acknowledgements This work was supported by the project SYMBIT reg. number: CZ.02.1.01/0.0/0.0/15_003/0000477 financed by the ERDF; grant 16-13721S from the Czech Science Foundation and by project LO1305 of the Ministry of Education, Youth and Sports of the Czech Republic.

Supporting Information Complementary numerical data for the SAPT figures. MM point-charge values. Additional statistics for the inter- and intra-strand stacking. Details for the vertical separation interpolation. Values and origin of all D3 and NL parameters. Information about the former (2006) reference structure. (pdf).



Cartesian coordinates of all 10 base-pair steps and numerical reference data for all 40 base-stacking energies. (zip).

References

(1) (2) (3)

(4) (5) (6) (7) (8) (9)

(10) (11) (12) (13) (14) (15) (16) (17) (18) (19) (20) (21)

Neidle, S. Principles of Nucleic Acid Structure; Elsevier, 2008. Nerenberg, P. S.; Head-Gordon, T. New Developments in Force Fields for Biomolecular Simulations. Curr. Opin. Struct. Biol. 2018, 49, 129–138. Yildirim, I.; Turner, D. H. RNA Challenges for Computational Chemists †. Biochemistry 2005, 44, 13225– 13234. Šponer, J.; Šponer, J. E.; Petrov, A. I.; Leontis, N. B. Quantum Chemical Studies of Nucleic Acids: Can We Construct a Bridge to the RNA Structural Biology and Bioinformatics Communities? J. Phys. Chem. B 2010, 114, 15723–15741. Sponer, J.; Leszczynski, J.; Hobza, P. Nature of Nucleic Acid-Base Stacking: Nonempirical Ab Initio and Empirical Potential Characterization of 10 Stacked Base Dimers. Comparison of Stacked and H-Bonded Base Pairs. J. Phys. Chem. 1996, 100, 5590–5596. Grimme, S. Do Special Noncovalent π-π Stacking Interactions Really Exist? Angew. Chemie - Int. Ed. 2008, 47, 3430–3434. Lipparini, F.; Scalmani, G.; Mennucci, B. Non Covalent Interactions in RNA and DNA Base Pairs: A Quantum-Mechanical Study of the Coupling between Solvent and Electronic Density. Phys. Chem. Chem. Phys. 2009, 11, 11617. Hobza, P. Stacking Interactions. Phys. Chem. Chem. Phys. 2008, 10, 2581. Banáš, P.; Mládek, A.; Otyepka, M.; Zgarbová, M.; Jurečka, P.; Svozil, D.; Lankaš, F.; Šponer, J. Can We Accurately Describe the Structure of Adenine Tracts in B-DNA? Reference Quantum-Chemical Computations Reveal Overstabilization of Stacking by Molecular Mechanics. J. Chem. Theory Comput. 2012, 8, 2448–2460. Parker, T. M.; Sherrill, C. D. Assessment of Empirical Models versus High-Accuracy Ab Initio Methods for Nucleobase Stacking: Evaluating the Importance of Charge Penetration. J. Chem. Theory Comput. 2015, 11, 4197–4204. Morgado, C. A.; Jurecka, P.; Svozil, D.; Hobza, P.; Sponer, J. Balance of Attraction and Repulsion in NucleicAcid Base Stacking: CCSD(T)/Complete-Basis-Set-Limit Calculations on Uracil Dimer and a Comparison with the Force-Field Description. J. Chem. Theory Comput. 2009, 5, 1524–1544. Zgarbová, M.; Otyepka, M.; Šponer, J.; Hobza, P.; Jurečka, P. Large-Scale Compensation of Errors in Pairwise-Additive Empirical Force Fields: Comparison of AMBER Intermolecular Terms with Rigorous DFTSAPT Calculations. Phys. Chem. Chem. Phys. 2010, 12, 10476. Wang, Q.; Rackers, J. A.; He, C.; Qi, R.; Narth, C.; Lagardere, L.; Gresh, N.; Ponder, J. W.; Piquemal, J.-P.; Ren, P. General Model for Treating Short-Range Electrostatic Penetration in a Molecular Mechanics Force Field. J. Chem. Theory Comput. 2015, 11, 2609–2618. Řezáč, J.; Hobza, P. Benchmark Calculations of Interaction Energies in Noncovalent Complexes and Their Applications. Chem. Rev. 2016, 116, 5038–5071. Goerigk, L.; Hansen, A.; Bauer, C.; Ehrlich, S.; Najibi, A.; Grimme, S. A Look at the Density Functional Theory Zoo with the Advanced GMTKN55 Database for General Main Group Thermochemistry, Kinetics and Noncovalent Interactions. Phys. Chem. Chem. Phys. 2017, 19, 32184–32215. Zhao, Y.; Truhlar, D. G. How Well Can New-Generation Density Functional Methods Describe Stacking Interactions in Biological Systems? Phys. Chem. Chem. Phys. 2005, 7, 2701. Deprince, A. E.; Sherrill, C. D. Accurate Noncovalent Interaction Energies Using Truncated Basis Sets Based on Frozen Natural Orbitals. J. Chem. Theory Comput. 2013, 9, 293–299. Rezac, 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. Kesharwani, M. K.; Karton, A.; Sylvetsky, N.; Martin, J. M. L. The S66 Non-Covalent Interactions Benchmark Reconsidered Using Explicitly Correlated Methods Near the Basis Set Limit. Aust. J. Chem. 2018, 71, 238–248. Goerigk, L.; Kruse, H.; Grimme, S. Benchmarking Density Functional Methods against the S66 and S66x8 Datasets for Non-Covalent Interactions. ChemPhysChem 2011, 12, 3421–3433. 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. ACS Paragon Plus Environment

Page 36 of 43

Page 37 of 43 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


(22) (23) (24) (25) (26) (27) (28) (29) (30) (31) (32) (33) (34) (35) (36) (37) (38) (39) (40) (41) (42) (43) (44) (45)

Poater, J.; Swart, M.; Bickelhaupt, F. M.; Fonseca Guerra, C. B-DNA Structure and Stability: The Role of Hydrogen Bonding, π-π Stacking Interactions, Twist-Angle, and Solvation. Org. Biomol. Chem. 2014, 12, 4691–4700. Hamlin, T. A.; Poater, J.; Fonseca Guerra, C.; Bickelhaupt, F. M. B-DNA Model Systems in Non-Terran BioSolvents: Implications for Structure, Stability and Replication. Phys. Chem. Chem. Phys. 2017, 19, 16969– 16978. Šponer, J.; Šponer, J. E.; Mládek, A.; Jurečka, P.; Banáš, P.; Otyepka, M. Nature and Magnitude of Aromatic Base Stacking in DNA and RNA: Quantum Chemistry, Molecular Mechanics, and Experiment. Biopolymers 2013, 99, 978–988. Morgado, C. A.; Jurečka, P.; Svozil, D.; Hobza, P.; Sponer, J. Balance of Attraction and Repulsion in NucleicAcid Base Stacking: CCSD(T)/Complete-Basis-Set-Limit Calculations on Uracil Dimer and a Comparison with the Force-Field Description. J. Chem. Theory Comput. 2009, 5, 1524–1544. Hobza, P.; Sponer, J.; Polasek, M. H-Bonded and Stacked DNA-Base Pairs - Cytosine Dimer - an Ab-Initio 2nd-Order Moller-Plesset Study. J. Am. Chem. Soc. 1995, 117, 792–798. Sponer, J.; Gabb, H. A.; Leszczynski, J.; Hobza, P. Base-Base and Deoxyribose-Base Stacking Interactions in B-DNA and Z-DNA: A Quantum-Chemical Study. Biophys. J. 1997, 73, 76–87. Šponer, J.; Riley, K. E.; Hobza, P. Nature and Magnitude of Aromatic Stacking of Nucleic Acid Bases. Phys. Chem. Chem. Phys. 2008, 10, 2595. Šponer, J.; Jurečka, P.; Marchan, I.; Luque, F. J.; Orozco, M.; Hobza, P. Nature of Base Stacking: Reference Quantum-Chemical Stacking Energies in Ten Unique B-DNA Base-Pair Steps. Chem. - A Eur. J. 2006, 12, 2854–2865. Leininger, M. L.; Nielsen, I. M. B.; Colvin, M. E.; Janssen, C. L. Accurate Structures and Binding Energies for Stacked Uracil Dimers. J. Phys. Chem. A 2002, 106, 3850–3854. Waller, M. P.; Robertazzi, A.; Platts, J. A.; Hibbs, D. E.; Williams, P. A. Hybrid Density Functional Theory for π-Stacking Interactions: Application to Benzenes, Pyridines, and DNA Bases. J. Comput. Chem. 2006, 27, 491–504. Cybulski, H.; Sadlej, J. Symmetry-Adapted Perturbation-Theory Interaction-Energy Decomposition for Hydrogen-Bonded and Stacking Structures. J. Chem. Theory Comput. 2008, 4, 892–897. McDonald, A. R.; Denning, E. J.; Mackerell, A. D. Impact of Geometry Optimization on Base-Base Stacking Interaction Energies in the Canonical A- and B-Forms of DNA. J. Phys. Chem. A 2013, 117, 1560–1568. Frey, J. A.; Müller, A.; Losada, M.; Leutwyler, S. Isomers of the Uracil Dimer: An Ab Initio Benchmark Study. J. Phys. Chem. B 2007, 111, 3534–3542. Hunter, R. S.; Van Mourik, T. DNA Base Stacking: The Stacked Uracil/Uracil and Thymine/Thymine Minima. J. Comput. Chem. 2012, 33, 2161–2172. Hobza, P.; Sponer, J. Structure, Energetics, and Dynamics of the Nucleic Acid Base Pairs: Nonempirical Ab Initio Calculations. Chem. Rev. 1999, 99, 3247–3276. Parker, T. M.; Hohenstein, E. G.; Parrish, R. M.; Hud, N. V.; Sherrill, C. D. Quantum-Mechanical Analysis of the Energetic Contributions to ?? Stacking in Nucleic Acids versus Rise, Twist, and Slide. J. Am. Chem. Soc. 2013, 135, 1306–1316. Hill, J. G.; Platts, J. A. Spin-Component Scaling Methods for Weak and Stacking Interactions. J. Chem. Theory Comput. 2007, 3, 80–85. Mukherjee, S.; Kailasam, S.; Bansal, M.; Bhattacharyya, D. Stacking Interactions in RNA and DNA: RollSlide Energy Hyperspace for Ten Unique Dinucleotide Steps. Biopolymers 2014, 103, 1–30. Barone, G.; Fonseca Guerra, C.; Bickelhaupt, F. M. B-DNA Structure and Stability as Function of Nucleic Acid Composition: Dispersion-Corrected DFT Study of Dinucleoside Monophosphate Single and Double Strands. ChemistryOpen 2013, 2, 186–193. von Lilienfeld, O. A.; Tkatchenko, A. Two- and Three-Body Interatomic Dispersion Energy Contributions to Binding in Molecules and Solids. J. Chem. Phys. 2010, 132. Řezáč, 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. Mahadevi, A. S.; Sastry, G. N. Cooperativity in Noncovalent Interactions. Chemical Reviews. American Chemical Society March 9, 2016, pp 2775–2825. Hill, J. G.; Platts, J. A. Calculating Stacking Interactions in Nucleic Acid Base-Pair Steps Using SpinComponent Scaling and Local Second Order Møller–Plesset Perturbation Theory. Phys. Chem. Chem. Phys. 2008, 10, 2785. Dickerson, R. E. Base Sequence and Helix Structure Variation in B and A DNA. J. Mol. Biol. 1983, 166, 419– 441. ACS Paragon Plus Environment


(46) (47) (48) (49) (50) (51) (52) (53) (54) (55) (56) (57) (58) (59)

(60) (61) (62) (63) (64) (65) (66) (67) (68)

Grimme, S. Dispersion Interaction and Chemical Bonding. In The Chemical Bond: Chemical Bonding Across the Periodic Table; Franking, G., Shaik, S., Eds.; Wiley-VCH Verlag GmbH & Co. KGaA: Weinheim, Germany, 2014; pp 477–500. Grimme, S. Density Functional Theory with London Dispersion Corrections. Wiley Interdiscip. Rev. Mol. Sci. 2011, 1, 211–228. Brauer, B.; Kesharwani, M. K.; Kozuch, S.; Martin, J. M. L. The S66x8 Benchmark for Noncovalent Interactions Revisited: Explicitly Correlated: Ab Initio Methods and Density Functional Theory. Phys. Chem. Chem. Phys. 2016, 18, 20905–20925. Grimme, S.; Goerigk, L.; Fink, R. F. Spin-Component-Scaled Electron Correlation Methods. Wiley Interdiscip. Rev. Comput. Mol. Sci. 2012, 2, 886–906. Lu, X. J.; Olson, W. K. 3DNA: A Software Package for the Analysis, Rebuilding and Visualization of ThreeDimensional Nucleic Acid Structures. Nucleic Acids Res. 2003, 31, 5108–5121. Boys, S. F.; Bernardi, F. The Calculation of Small Molecular Interactions by the Differences of Separate Total Energies. Some Procedures with Reduced Errors. Mol. Phys. 1970, 19, 553–566. Pettersen, E. F.; Goddard, T. D.; Huang, C. C.; Couch, G. S.; Greenblatt, D. M.; Meng, E. C.; Ferrin, T. E. UCSF Chimera - A Visualization System for Exploratory Research and Analysis. J. Comput. Chem. 2004, 25, 1605–1612. Papajak, E.; Truhlar, D. G. Convergent Partially Augmented Basis Sets for Post-Hartree−Fock Calculations of Molecular Properties and Reaction Barrier Heights. J. Chem. Theory Comput. 2011, 7, 10–18. Dunning, T. H. Gaussian Basis Sets for Use in Correlated Molecular Calculations. I. The Atoms Boron through Neon and Hydrogen. J. Chem. Phys. 1989, 90, 1007–1023. Kendall, R. A.; Dunning, T. H.; Harrison, R. J. Electron Affinities of the First-row Atoms Revisited. Systematic Basis Sets and Wave Functions. J. Chem. Phys. 1992, 96, 6796–6806. Peterson, K. A.; Adler, T. B.; Werner, H. J. Systematically Convergent Basis Sets for Explicitly Correlated Wavefunctions: The Atoms H, He, B-Ne, and Al-Ar. J. Chem. Phys. 2008, 128. Weigend, F.; Ahlrichs, R. Balanced Basis Sets of Split Valence, Triple Zeta Valence and Quadruple Zeta Valence Quality for H to Rn: Design and Assessment of Accuracy. Phys. Chem. Chem. Phys. 2005, 7, 3297– 3305. Neese, F.; Valeev, E. F. Revisiting the Atomic Natural Orbital Approach for Basis Sets: Robust Systematic Basis Sets for Explicitly Correlated and Conventional Correlated Ab Initio Methods? J. Chem. Theory Comput. 2011, 7, 33–43. Parrish, R. M.; Burns, L. A.; Smith, D. G. A.; Simmonett, A. C.; DePrince, A. E.; Hohenstein, E. G.; Bozkaya, U.; Sokolov, A. Y.; Di Remigio, R.; Richard, R. M.; et al. Psi4 1.1: An Open-Source Electronic Structure Program Emphasizing Automation, Advanced Libraries, and Interoperability. J. Chem. Theory Comput. 2017, 13, 3185–3197. Marques, M. A. L.; Oliveira, M. J. T.; Burnus, T. Libxc: A Library of Exchange and Correlation Functionals for Density Functional Theory. Comput. Phys. Commun. 2012, 183, 2227–2281. Hujo, W.; Grimme, S. Performance of the van Der Waals Density Functional VV10 and (Hybrid)GGA Variants for Thermochemistry and Noncovalent Interactions. J. Chem. Theory Comput. 2011, 7, 3866– 3871. Parker, T. M.; Burns, L. A.; Parrish, R. M.; Ryno, A. G.; Sherrill, C. D. Levels of Symmetry Adapted Perturbation Theory (SAPT). I. Efficiency and Performance for Interaction Energies. J. Chem. Phys. 2014, 140, 094106. Hohenstein, E. G.; Sherrill, C. D. Wavefunction Methods for Noncovalent Interactions. Wiley Interdiscip. Rev. Comput. Mol. Sci. 2012, 2, 304–326. Hohenstein, E. G.; Sherrill, C. D. Efficient Evaluation of Triple Excitations in Symmetry-Adapted Perturbation Theory via Second-Order Møller–Plesset Perturbation Theory Natural Orbitals. J. Chem. Phys. 2010, 133, 104107. Hohenstein, E. G.; Parrish, R. M.; Sherrill, C. D.; Turney, J. M.; Schaefer, H. F. Large-Scale SymmetryAdapted Perturbation Theory Computations via Density Fitting and Laplace Transformation Techniques: Investigating the Fundamental Forces of DNA-Intercalator Interactions. J. Chem. Phys. 2011, 135, 174107. Distasio, R. A.; Head-Gordon, M. Optimized Spin-Component Scaled Second-Order Mller-Plesset Perturbation Theory for Intermolecular Interaction Energies. Mol. Phys. 2007, 105, 1073–1083. McGibbon, R. T.; Taube, A. G.; Donchev, A. G.; Siva, K.; Hernández, F.; Hargus, C.; Law, K.-H.; Klepeis, J. L.; Shaw, D. E. Improving the Accuracy of Møller-Plesset Perturbation Theory with Neural Networks. J. Chem. Phys. 2017, 147, 161725. TURBOMOLE V7.2.1. a development of University of Karlsruhe and Forschungszentrum Karlsruhe GmbH, 1989-2007, TURBOMOLE GmbH, since 2007; available from http://www.turbomole.com. ACS Paragon Plus Environment

Page 38 of 43

Page 39 of 43 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


(69) (70) (71) (72) (73) (74) (75) (76) (77) (78) (79) (80) (81) (82) (83) (84) (85) (86) (87) (88)

(89) (90) (91) (92) (93)

Ahlrichs, R.; Bar, M.; Marco, H.; Horn, H.; Kolmel, C. ELECTRONIC STRUCTURE CALCULATIONS ON WORKSTATION COMPUTERS: THE PROGRAM SYSTEM TURBOMOLE Reinhart AHLRICHS, Michael BAR, Marco H;ISER, Hans HORN and Christoph KtjLMEL. Chem. Phys. Lett. 1989, 162, 165. Riplinger, C.; Pinski, P.; Becker, U.; Valeev, E. F.; Neese, F. Sparse Maps - A Systematic Infrastructure for Reduced-Scaling Electronic Structure Methods. II. Linear Scaling Domain Based Pair Natural Orbital Coupled Cluster Theory. J. Chem. Phys. 2016, 144. Pavošević, F.; Pinski, P.; Riplinger, C.; Neese, F.; Valeev, E. F. SparseMaps - A Systematic Infrastructure for Reduced-Scaling Electronic Structure Methods. IV. Linear-Scaling Second-Order Explicitly Correlated Energy with Pair Natural Orbitals. J. Chem. Phys. 2016, 144. Neese, F. Software Update: The ORCA Program System, Version 4.0. Wiley Interdiscip. Rev. Comput. Mol. Sci. 2017, e1327. Liakos, D. G.; Neese, F. Is It Possible to Obtain Coupled Cluster Quality Energies at near Density Functional Theory Cost? Domain-Based Local Pair Natural Orbital Coupled Cluster vs Modern Density Functional Theory. J. Chem. Theory Comput. 2015, 11, 4054–4063. Liakos, D. G.; Izsák, R.; Valeev, E. F.; Neese, F. What Is the Most Efficient Way to Reach the Canonical MP2 Basis Set Limit? Mol. Phys. 2013, 111, 2653–2662. Yousaf, K. E.; Peterson, K. A. Optimized Auxiliary Basis Sets for Explicitly Correlated Methods. J. Chem. Phys. 2008, 129. Grimme, S.; Bannwarth, C.; Caldeweyher, E.; Pisarek, J.; Hansen, A. A General Intermolecular Force Field Based on Tight-Binding Quantum Chemical Calculations. J. Chem. Phys. 2017, 147, 161708. Brandenburg, J. G.; Grimme, S. Accurate Modeling of Organic Molecular Crystals by Dispersion-Corrected Density Functional Tight Binding (DFTB). J. Phys. Chem. Lett. 2014, 5, 1785–1789. Gaus, M.; Cui, Q.; Elstner, M. DFTB3: Extension of the Self-Consistent-Charge Density-Functional TightBinding Method (SCC-DFTB). J. Chem. Theory Comput. 2011, 7, 931–948. Aradi, B.; Hourahine, B.; Frauenheim, T. DFTB+, a Sparse Matrix-Based Implementation of the DFTB Method †. J. Phys. Chem. A 2007, 111, 5678–5684. Kromann, J. C.; Christensen, A. S.; Steinmann, C.; Korth, M.; Jensen, J. H. A Third-Generation Dispersion and Third-Generation Hydrogen Bonding Corrected PM6 Method: PM6-D3H+. PeerJ 2014, 2, e449. Stewart, J. J. P. MOPAC2016. MOPAC2016: Colorado Springs, CO, USA 2016. Cornell, W. D.; Cieplak, P.; Bayly, C. I.; Gould, I. R.; Merz, K. M.; Ferguson, D. M.; Spellmeyer, D. C.; Fox, T.; Caldwell, J. W.; Kollman, P. A. A Second Generation Force Field for the Simulation of Proteins, Nucleic Acids, and Organic Molecules. J. Am. Chem. Soc. 1995, 117, 5179–5197. Wang, J.; Wang, W.; Kollman, P. A.; Case, D. A. Automatic Atom Type and Bond Type Perception in Molecular Mechanical Calculations. J. Mol. Graph. Model. 2006, 25, 247–260. Frisch, M. J.; Trucks, G. W.; Schlegel, H. B.; Scuseria, G. E.; Robb, M. A.; Cheeseman, J. R.; Scalmani, G.; Barone, V.; Mennucci, B.; Petersson, G. A.; et al. Gaussian 09. Gaussian, Inc.: Wallingford, CT, USA 2009. Hart, K.; Foloppe, N.; Baker, C. M.; Denning, E. J.; Nilsson, L.; Mackerell Jr., A. D. Optimization of the CHARMM Additive Force Field for DNA: Improved Treatment of the BI/BII Conformational Equilibrium. J Chem Theory Comput 2012, 8, 348–362. Oostenbrink, C.; Villa, A.; Mark, A. E.; Van Gunsteren, W. F. A Biomolecular Force Field Based on the Free Enthalpy of Hydration and Solvation: The GROMOS Force-Field Parameter Sets 53A5 and 53A6. J. Comput. Chem. 2004, 25, 1656–1676. Grimme, S. Semiempirical Hybrid Density Functional with Perturbative Second-Order Correlation. J. Chem. Phys. 2006, 124. Karton, A.; Tarnopolsky, A.; Lamère, J. F.; Schatz, G. C.; Martin, J. M. L. Highly Accurate First-Principles Benchmark Data Sets for the Parametrization and Validation of Density Functional and Other Approximate Methods. Derivation of a Robust, Generally Applicable, Double-Hybrid Functional for Thermochemistry and Thermochemical . J. Phys. Chem. A 2008, 112, 12868–12886. Kozuch, S.; Gruzman, D.; Martin, J. M. L. DSD-BLYP: A General Purpose Double Hybrid Density Functional Including Spin Component Scaling and Dispersion Correction †. J. Phys. Chem. C 2010, 114, 20801–20808. Kozuch, S.; Martin, J. M. L. DSD-PBEP86: In Search of the Best Double-Hybrid DFT with Spin-Component Scaled MP2 and Dispersion Corrections. Phys. Chem. Chem. Phys. 2011, 13, 20104. Goerigk, L.; Grimme, S. Efficient and Accurate Double-Hybrid-Meta-GGA Density Functionals-Evaluation with the Extended GMTKN30 Database for General Main Group Thermochemistry, Kinetics, and Noncovalent Interactions. J. Chem. Theory Comput. 2011, 7, 291–309. Kozuch, S.; Martin, J. M. L. Spin-Component-Scaled Double Hybrids: An Extensive Search for the Best Fifth-Rung Functionals Blending DFT and Perturbation Theory. J. Comput. Chem. 2013, 34, 2327–2344. Becke, A. D. Density-Functional Thermochemistry. III. The Role of Exact Exchange. J. Chem. Phys. 1993, ACS Paragon Plus Environment


(94) (95)

(96) (97) (98) (99) (100) (101) (102) (103) (104) (105) (106) (107) (108) (109) (110) (111) (112) (113) (114) (115) (116) (117)

98, 5648. Stephens, P. J.; Devlin, F. J.; Chabalowski, C. F.; Frisch, M. J. Ab Initio Calculation of Vibrational Absorption and Circular Dichroism Spectra Using Density Functional Force Fields. J. Phys. Chem. 1994, 98, 11623– 11627. Zhao, Y.; Truhlar, D. G. The M06 Suite of Density Functionals for Main Group Thermochemistry, Thermochemical Kinetics, Noncovalent Interactions, Excited States, and Transition Elements: Two New Functionals and Systematic Testing of Four M06-Class Functionals and 12 Other Function. Theor. Chem. Acc. 2008, 120, 215–241. Zhao, Y.; Schultz, N. E.; Truhlar, D. G. Design of Density Functionals by Combining the Method of Constraint Satisfaction with Parametrization for Thermochemistry, Thermochemical Kinetics, and Noncovalent Interactions. J. Chem. Theory Comput. 2006, 2, 364–382. Peverati, R.; Truhlar, D. G. Improving the Accuracy of Hybrid Meta-GGA Density Functionals by Range Separation. J. Phys. Chem. Lett. 2011, 2, 2810–2817. Peverati, R.; Truhlar, D. G. Communication: A Global Hybrid Generalized Gradient Approximation to the Exchange-Correlation Functional That Satisfies the Second-Order Density-Gradient Constraint and Has Broad Applicability in Chemistry. J. Chem. Phys. 2011, 135. Yu, H. S.; He, X.; Li, S. L.; Truhlar, D. G. MN15: A Kohn–Sham Global-Hybrid Exchange–correlation Density Functional with Broad Accuracy for Multi-Reference and Single-Reference Systems and Noncovalent Interactions. Chem. Sci. 2016, 7, 5032–5051. Adamo, C.; Barone, V. Toward Reliable Density Functional Methods without Adjustable Parameters: The PBE0 Model. J. Chem. Phys. 1999, 110, 6158. Zhang, Y.; Yang, W. Comment on “Generalized Gradient Approximation Made Simple.” Physical Review Letters. 1998, p 890. Vydrov, O. A.; Van Voorhis, T. Nonlocal van Der Waals Density Functional: The Simpler the Better. J. Chem. Phys. 2010, 133. Mardirossian, N.; Head-Gordon, M. ω B97M-V: A Combinatorially Optimized, Range-Separated Hybrid, Meta-GGA Density Functional with VV10 Nonlocal Correlation. J. Chem. Phys. 2016, 144. Zhao, Y.; Truhlar, D. G. Design of Density Functionals That Are Broadly Accurate for Thermochemistry, Thermochemical Kinetics, and Nonbonded Interactions. J. Phys. Chem. A 2005, 109, 5656–5667. Weintraub, E.; Henderson, T. M.; Scuseria, G. E. Long-Range-Corrected Hybrids Based on a New Model Exchange Hole. J. Chem. Theory Comput. 2009, 5, 754–762. Yanai, T.; Tew, D. P.; Handy, N. C. A New Hybrid Exchange-Correlation Functional Using the CoulombAttenuating Method (CAM-B3LYP). Chem. Phys. Lett. 2004, 393, 51–57. Becke, A. D. DENSITY-FUNCTIONAL EXCHANGE-ENERGY APPROXIMATION WITH CORRECT ASYMPTOTICBEHAVIOR. Phys. Rev. A 1988, 38, 3098–3100. Lee, C. T.; Yang, W. T.; Parr, R. G. DEVELOPMENT OF THE COLLE-SALVETTI CORRELATION-ENERGY FORMULA INTO A FUNCTIONAL OF THE ELECTRON-DENSITY. Phys. Rev. B 1988, 37, 785–789. Sun, J.; Marsman, M.; Csonka, G. I.; Ruzsinszky, A.; Hao, P.; Kim, Y. S.; Kresse, G.; Perdew, J. P. SelfConsistent Meta-Generalized Gradient Approximation within the Projector-Augmented-Wave Method. Phys. Rev. B - Condens. Matter Mater. Phys. 2011, 84, 035117. Tao, J. M.; Perdew, J. P.; Staroverov, V. N.; Scuseria, G. E. Climbing the Density Functional Ladder: Nonempirical Meta-Generalized Gradient Approximation Designed for Molecules and Solids. Phys. Rev. Lett. 2003, 91, 4. Grimme, S. Semiempirical GGA-Type Density Functional Constructed with a Long-Range Dispersion Correction. J. Comput. Chem. 2006, 27, 1787–1799. Perdew, J. P.; Burke, K.; Ernzerhof, M. Generalized Gradient Approximation Made Simple. Phys. Rev. Lett. 1996, 77, 3865–3868. Sun, J.; Ruzsinszky, A.; Perdew, J. Strongly Constrained and Appropriately Normed Semilocal Density Functional. Phys. Rev. Lett. 2015, 115, 1–6. Perdew, J. P. Density-Functional Approximation for the Correlation Energy of the Inhomogeneous Electron Gas. Phys. Rev. B 1986, 33, 8822–8824. Yu, H. S.; He, X.; Truhlar, D. G. MN15-L: A New Local Exchange-Correlation Functional for Kohn-Sham Density Functional Theory with Broad Accuracy for Atoms, Molecules, and Solids. J. Chem. Theory Comput. 2016, 12, 1280–1293. Mardirossian, N.; Head-Gordon, M. Mapping the Genome of Meta-Generalized Gradient Approximation Density Functionals: The Search for B97M-V. J. Chem. Phys. 2015, 142, 074111. Najibi, A.; Goerigk, L. The Nonlocal Kernel in van Der Waals Density Functionals as an Additive Correction: An Extensive Analysis with Special Emphasis on the B97M-V and ΩB97M-V Approaches. J. Chem. Theory ACS Paragon Plus Environment

Page 40 of 43

Page 41 of 43 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


(118) (119) (120) (121) (122) (123) (124) (125) (126) (127) (128) (129) (130) (131) (132) (133) (134) (135) (136) (137) (138) (139) (140) (141) (142)

Comput. 2018, 14, 5725–5738. Marshall, M. S.; Burns, L. A.; Sherrill, C. D. Basis Set Convergence of the Coupled-Cluster Correction, ΔMP2CCSD(T): Best Practices for Benchmarking Non-Covalent Interactions and the Attendant Revision of the S22, NBC10, HBC6, and HSG Databases. J. Chem. Phys. 2011, 135, 194102. van Duijneveldt, F. B.; van de Rijdt, J. G. C. M. va. D.; van Lenthe, J. H. State of the Art in Counterpoise Theory. Chem. Rev. 1994, 94, 1873–1885. Halkier, A.; Helgaker, T.; Jorgensen, P.; Klopper, W.; Koch, H.; Olsen, J.; Wilson, A. K. Basis-Set Convergence in Correlated Calculations on Ne, N-2, and H2O. Chem. Phys. Lett. 1998, 286, 243–252. Burns, L. A.; Marshall, M. S.; Sherrill, C. D. Comparing Counterpoise-Corrected, Uncorrected, and Averaged Binding Energies for Benchmarking Noncovalent Interactions. J. Chem. Theory Comput. 2014, 10, 49–57. Brauer, B.; Kesharwani, M. K.; Martin, J. M. L. Some Observations on Counterpoise Corrections for Explicitly Correlated Calculations on Noncovalent Interactions. J. Chem. Theory Comput. 2014, 10, 3791– 3799. Kong, L.; Bischoff, F. A.; Valeev, E. F. Explicitly Correlated R12/F12 Methods for Electronic Structure. Chem. Rev. 2011, 112, 75–107. Marchetti, O.; Werner, H.-J. Accurate Calculations of Intermolecular Interaction Energies Using Explicitly Correlated Wave Functions. Phys. Chem. Chem. Phys. 2008, 10, 3400. Köhn, A. Explicitly Correlated Connected Triple Excitations in Coupled-Cluster Theory. J. Chem. Phys. 2009, 130, 131101. Schmitz, G.; Hättig, C. Perturbative Triples Correction for Local Pair Natural Orbital Based Explicitly Correlated CCSD(F12*) Using Laplace Transformation Techniques. J. Chem. Phys. 2016, 145, 234107. Ma, Q.; Werner, H.-J. Scalable Electron Correlation Methods. 5. Parallel Perturbative Triples Correction for Explicitly Correlated Local Coupled Cluster with Pair Natural Orbitals. J. Chem. Theory Comput. 2018, 14, 198–215. Guo, Y.; Riplinger, C.; Becker, U.; Liakos, D. G.; Minenkov, Y.; Cavallo, L.; Neese, F. Communication: An Improved Linear Scaling Perturbative Triples Correction for the Domain Based Local Pair-Natural Orbital Based Singles and Doubles Coupled Cluster Method [DLPNO-CCSD(T)]. J. Chem. Phys. 2018, 148, 011101. Riplinger, C.; Sandhoefer, B.; Hansen, A.; Neese, F. Natural Triple Excitations in Local Coupled Cluster Calculations with Pair Natural Orbitals. J. Chem. Phys. 2013, 139, 134101. Cieplak, P.; Cornell, W. D.; Bayly, C.; Kollman, P. A. Application of the Multimolecule and Multiconformational RESP Methodology to Biopolymers: Charge Derivation for DNA, RNA, and Proteins. J. Comput. Chem. 1995, 16, 1357–1377. Perdew, J. P.; Schmidt, K. Jacob’s Ladder of Density Functional Approximations for the ExchangeCorrelation Energy. AIP Conf. Proc. 2001, 577, 1–20. Hait, D.; Head-Gordon, M. XDH Double Hybrid Functionals Can Be Qualitatively Incorrect for NonEquilibrium Geometries: Dipole Moment Inversion and Barriers to Radical-Radical Association Using XYG3 and XYGJ-OS. J. Chem. Phys. 2018, 148, 171102. Cohen, A. J.; Mori-Sánchez, P.; Yang, W. Challenges for Density Functional Theory. Chem. Rev. 2012, 112, 289–320. Goerigk, L. Treating London-Dispersion Effects with the Latest Minnesota Density Functionals: Problems and Possible Solutions. J. Phys. Chem. Lett. 2015, 6, 3891–3896. Brandenburg, J. G.; Bates, J. E.; Sun, J.; Perdew, J. P. Benchmark Tests of a Strongly Constrained Semilocal Functional with a Long-Range Dispersion Correction. Phys. Rev. B 2016, 94, 115144. Johnson, E. R.; Becke, A. D.; Sherrill, C. D.; DiLabio, G. A. Oscillations in Meta-Generalized-Gradient Approximation Potential Energy Surfaces for Dispersion-Bound Complexes. J. Chem. Phys. 2009, 131. Smith, D. G. A.; Burns, L. A.; Patkowski, K.; Sherrill, C. D. Revised Damping Parameters for the D3 Dispersion Correction to Density Functional Theory. J. Phys. Chem. Lett. 2016, acs.jpclett.6b00780. Hujo, W.; Grimme, S. Comparison of the Performance of Dispersion-Corrected Density Functional Theory for Weak Hydrogen Bonds. Phys. Chem. Chem. Phys. 2011, 13, 13942. Goerigk, L. How Do DFT-DCP, DFT-NL, and DFT-D3 Compare for the Description of London-Dispersion Effects in Conformers and General Thermochemistry? J. Chem. Theory Comput. 2014, 140206112919009. Kruse, H.; Grimme, S. A Geometrical Correction for the Inter- and Intra-Molecular Basis Set Superposition Error in Hartree-Fock and Density Functional Theory Calculations for Large Systems. J. Chem. Phys. 2012, 136, 154101. Caldeweyher, E.; Brandenburg, J. G. Simplified DFT Methods for Consistent Structures and Energies of Large Systems. J. Phys. Condens. Matter 2018, 30, 213001. Thiel, W. Semiempirical Quantum-Chemical Methods. Wiley Interdiscip. Rev. Comput. Mol. Sci. 2014, 4, ACS Paragon Plus Environment


(143) (144) (145) (146) (147) (148) (149)

(150) (151) (152) (153) (154) (155) (156) (157) (158) (159) (160)

145–157. Kruse, H.; Havrila, M.; Šponer, J. QM Computations on Complete Nucleic Acids Building Blocks: Analysis of the Sarcin–Ricin RNA Motif Using DFT-D3, HF-3c, PM6-D3H, and MM Approaches. J. Chem. Theory Comput. 2014, 10, 2615–2629. Kruse, H.; Mladek, A.; Gkionis, K.; Hansen, A.; Grimme, S.; Sponer, J. Quantum Chemical Benchmark Study on 46 RNA Backbone Families Using a Dinucleotide Unit. J. Chem. Theory Comput. 2015, 11, 4972–4991. Szabla, R.; Havrila, M.; Kruse, H.; Šponer, J. Comparative Assessment of Different RNA Tetranucleotides from the DFT-D3 and Force Field Perspective. J. Phys. Chem. B 2016, 120, 10635–10648. Miriyala, V. M.; Řezáč, J. Description of Non-Covalent Interactions in SCC-DFTB Methods. J. Comput. Chem. 2017, 38, 688–697. Miriyala, V. M.; Řezáč, J. Testing Semiempirical Quantum Mechanical Methods on a Data Set of Interaction Energies Mapping Repulsive Contacts in Organic Molecules. J. Phys. Chem. A 2018, 122, 2801– 2808. Sponer, J.; Bussi, G.; Krepl, M.; Banas, P.; Bottaro, S.; Cunha, R. A.; Gil-Ley, A.; Pinamonti, G.; Poblete, S.; Jurečka, P.; et al. RNA Structural Dynamics as Captured by Molecular Simulations: A Comprehensive Overview. Chem. Rev. 2018, 118, 4177–4338. Hobza, P.; Kabelac, M.; Sponer, J.; Mejzlik, P.; Vondrasek, J. Performance of Empirical Potentials (AMBER, CFF95, CVFF, CHARMM, OPLS, POLTEV), Semiempirical Quantum Chemical Methods (AM1, MNDO/M, PM3), and Ab Initio Hartree-Fock Method for Interaction of DNA Bases: Comparison with Nonempirical beyond Hartree-Fock Res. J. Comput. Chem. 1997, 18, 1136–1150. Šponer, J.; Banáš, P.; Jurečka, P.; Zgarbová, M.; Kührová, P.; Havrila, M.; Krepl, M.; Stadlbauer, P.; Otyepka, M. Molecular Dynamics Simulations of Nucleic Acids. from Tetranucleotides to the Ribosome. J. Phys. Chem. Lett. 2014, 5, 1771–1782. Šponer, J.; Hobza, P. MP2 and CCSD(T) Study on Hydrogen Bonding, Aromatic Stacking and Nonaromatic Stacking. Chem. Phys. Lett. 1997, 267, 263–270. Pitonak, M.; Neogrady, P.; Cerny, J.; Grimme, S.; Hobza, P. Scaled MP3 Non-Covalent Interaction Energies Agree Closely with Accurate CCSD(T) Benchmark Data. Chemphyschem 2009, 10, 282–289. Robert Jr., A. D.; Vivekanand, V. G.; Alexandre, T. Many-Body van Der Waals Interactions in Molecules and Condensed Matter. J. Phys. Condens. Matter 2014, 26, 213202. 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. Risthaus, T.; Grimme, S. Benchmarking of London Dispersion-Accounting Density Functional Theory Methods on Very Large Molecular Complexes. J. Chem. Theory Comput. 2013, 9, 1580–1591. Ouyang, J. F.; Bettens, R. P. A. Many-Body Basis Set Superposition Effect. J. Chem. Theory Comput. 2015, 11, 5132–5143. Richard, R. M.; Bakr, B. W.; Sherrill, C. D. Understanding the Many-Body Basis Set Superposition Error: Beyond Boys and Bernardi. J. Chem. Theory Comput. 2018, 14, 2386–2400. Wells, B.; Wilson, S. Van Der Waals Interaction Potentials Many-Body Effects. Mol. Phys. 1985, 55, 199– 210. Valiron, P.; Mayer, I. Hierarchy of Counterpoise Corrections for N-Body Clusters: Generalization of the Boys-Bernardi Scheme. Chem. Phys. Lett. 1997, 275, 46–55. Rackers, J. A.; Wang, Q.; Liu, C.; Piquemal, J.-P.; Ren, P.; Ponder, J. W. An Optimized Charge Penetration Model for Use with the AMOEBA Force Field. Phys. Chem. Chem. Phys. 2017, 19, 276–291.


Page 42 of 43

Page 43 of 43 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


for Table of Contents use only

Investigations of stacked DNA base-pair steps: highly-accurate stacking interaction energies, energy decomposition and many-body stacking effects

Holger Kruse*, Pavel Banas, Jiri Sponer


Investigations of stacked DNA base-pair steps: highly-accurate

Recommend Documents