MP2 and CCSD(T) complete basis set (CBS) limit interaction energies and geometries for more than 100 DNA base pairs, amino acid pairs and model complexes are for the first time presented together. Extrapolation to the CBS limit is done by using two-point extrapolation methods and different basis sets (aug-cc-pVDZ – aug-cc-pVTZ, aug-cc-pVTZ – aug-cc-pVQZ, cc-pVTZ – cc-pVQZ) are utilized. The CCSD(T) correction term, determined as a difference between CCSD(T) and MP2 interaction energies, is evaluated with smaller basis sets (6-31G** and cc-pVDZ). Two sets of complex geometries were used, optimized or experimental ones. The JSCH-2005 benchmark set, which is now available to the chemical community, can be used for testing lower-level computational methods. For the first screening the smaller training set (S22) containing 22 model complexes can be recommended. In this case larger basis sets were used for extrapolation to the CBS limit and also CCSD(T) and counterpoise-corrected MP2 optimized geometries were sometimes adopted.

We report a reparameterization of the glycosidic torsion χ of the Cornell et al. AMBER force field for RNA, χOL. The parameters remove destabilization of the anti region found in the ff99 force field and thus prevent formation of spurious ladder-like structural distortions in RNA simulations. They also improve the description of the syn region and the syn–anti balance as well as enhance MD simulations of various RNA structures. Although χOL can be combined with both ff99 and ff99bsc0, we recommend the latter. We do not recommend using χOL for B-DNA because it does not improve upon ff99bsc0 for canonical structures. However, it might be useful in simulations of DNA molecules containing syn nucleotides. Our parametrization is based on high-level QM calculations and differs from conventional parametrization approaches in that it incorporates some previously neglected solvation-related effects (which appear to be essential for obtaining correct anti/high-anti balance). Our χOL force field is compared with several previous glycosidic torsion parametrizations.

Abstract Density functional theory (DFT) methods, including nonlocal density gradient terms in the exchange and correlation energy functionals, were applied to various types of molecular clusters: H‐bonded, ionic, electrostatic, and London. Reliable results on the structure and stabilization energy were obtained for the first two types of cluster as long as Becke3LYP and Becke3P86 functionals and basis sets of at least DZ + P quality were used. DFT methods with currently available functionals failed completely, however, for London‐type clusters, for which no minimum was found on the potential energy surfaces. DFT interaction energy exhibits the same basis set extension dependence as the Hartree‐Fock (HF) interaction energy. Therefore, the Boys‐Bernardi function counterpoise procedure should be employed for elimination of the DFT basis set superposition error. © 1995 John Wiley & Sons, Inc.

Z-DNA duplexes are a particularly complicated test case for current force fields. We performed a set of explicit solvent molecular dynamics (MD) simulations with various AMBER force field parametrizations including our recent refinements of the ε/ζ and glycosidic torsions. None of these force fields described the ZI/ZII and other backbone substates correctly, and all of them underpredicted the population of the important ZI substate. We show that this underprediction can be attributed to an inaccurate potential for the sugar–phosphate backbone torsion angle β. We suggest a refinement of this potential, βOL1, which was derived using our recently introduced methodology that includes conformation-dependent solvation effects. The new potential significantly increases the stability of the dominant ZI backbone substate and improves the overall description of the Z-DNA backbone. It also has a positive (albeit small) impact on another important DNA form, the antiparallel guanine quadruplex (G-DNA), and improves the description of the canonical B-DNA backbone by increasing the population of BII backbone substates, providing a better agreement with experiment. We recommend using βOL1 in combination with our previously introduced corrections, εζOL1 and χOL4, (the combination being named OL15) as a possible alternative to the current β torsion potential for more accurate modeling of nucleic acids.

Hydrogen-bonded nucleic acids base pairs substantially contribute to the structure and stability of nucleic acids. The study presents reference ab initio structures and interaction energies of selected base pairs with binding energies ranging from −5 to −47 kcal/mol. The molecular structures are obtained using the RI-MP2 (resolution of identity MP2) method with extended cc-pVTZ basis set of atomic orbitals. The RI-MP2 method provides results essentially identical with the standard MP2 method. The interaction energies are calculated using the Complete Basis Set (CBS) extrapolation at the RI-MP2 level. For some base pairs, Coupled-Cluster corrections with inclusion of noniterative triple contributions (CCSD(T)) are given. The calculations are compared with selected medium quality methods. The PW91 DFT functional with the 6-31G** basis set matches well the RI-MP2/CBS absolute interaction energies and reproduces the relative values of base pairing energies with a maximum relative error of 2.6 kcal/mol when applied with Becke3LYP-optimized geometries. The Becke3LYP DFT functional underestimates the interaction energies by few kcal/mol with relative error of 2.2 kcal/mol. Very good performance of nonpolarizable Cornell et al. force field is confirmed and this indirectly supports the view that H-bonded base pairs are primarily stabilized by electrostatic interactions.

The utility of molecular dynamics (MD) simulations to model biomolecular structure, dynamics, and interactions has witnessed enormous advances in recent years due to the availability of optimized MD software and access to significant computational power, including GPU multicore computing engines and other specialized hardware. This has led researchers to routinely extend conformational sampling times to the microsecond level and beyond. The extended sampling time has allowed the community not only to converge conformational ensembles through complete sampling but also to discover deficiencies and overcome problems with the force fields. Accuracy of the force fields is a key component, along with sampling, toward being able to generate accurate and stable structures of biopolymers. The Amber force field for nucleic acids has been used extensively since the 1990s, and multiple artifacts have been discovered, corrected, and reassessed by different research groups. We present a direct comparison of two of the most recent and state-of-the-art Amber force field modifications, bsc1 and OL15, that focus on accurate modeling of double-stranded DNA. After extensive MD simulations with five test cases and two different water models, we conclude that both modifications are a remarkable improvement over the previous bsc0 force field. Both force field modifications show better agreement when compared to experimental structures. To ensure convergence, the Drew–Dickerson dodecamer (DDD) system was simulated using 100 independent MD simulations, each extended to at least 10 μs, and the independent MD simulations were concatenated into a single 1 ms long trajectory for each combination of force field and water model. This is significantly beyond the time scale needed to converge the conformational ensemble of the internal portions of a DNA helix absent internal base pair opening. Considering all of the simulations discussed in the current work, the MD simulations performed to assess and validate the current force fields and water models aggregate over 14 ms of simulation time. The results suggest that both the bsc1 and OL15 force fields render average structures that deviate significantly less than 1 Å from the average experimental structures. This can be compared to similar but less exhaustive simulations with the CHARMM 36 force field that aggregate to the ∼90 μs time scale and also perform well but do not produce structures as close to the DDD NMR average structures (with root-mean-square deviations of 1.3 Å) as the newer Amber force fields. On the basis of these analyses, any future research involving double-stranded DNA simulations using the Amber force fields should employ the bsc1 or OL15 modification.

Ab initio (MP2/6-31G*(0.25)) interaction energies were calculated for almost 240 geometries of 10 stacked nucleic acid−base pairs: A···A, C···C, G···G, U···U, A···C, G···A, A···U, G···C, C···U, and G···U; in some cases uracil was replaced by thymine. The most stable stacked pair is the G···G dimer (−11.3 kcal/mol), and the least stable is the uracil dimer (−6.5 kcal/mol). The interaction energies of H-bonded base pairs range from −25.8 kcal/mol (G···C) to −10.6 kcal/mol (T···T). The stability of stacked pairs originates in the electron correlation, while all the H-bonded pairs are dominated by the HF energy. The mutual orientation of the stacked bases is, however, primarily determined by the HF interaction energy. The ab initio base stacking energies are well reproduced by the empirical potential calculations, provided the atomic charges are derived by the same method as used in the ab initio calculations. Some contributions previously postulated to significantly influence base stacking (induction interactions, π−π interactions) were not found. Base stacking was also investigated in six B-DNA and two Z-DNA base pair steps; their geometries were taken from the oligonucleotide crystal data. The many-body correction was estimated at the HF/MINI-1 level. The sequence-dependent variations of the total base pair step stacking energies range from −9.9 to −14.7 kcal/mol. The range of the calculated many-body corrections to the stacking energy is 2 kcal/mol. The ab initio calculations exclude the consideration that the unusual conformational properties of the CpA(TpG) steps might be associated with attractive induction interactions of the exocyclic groups of DNA bases and the aromatic rings of bases.

Stacking energies in low-energy geometries of pyrimidine, uracil, cytosine, and guanine homodimers were determined by the MP2 and CCSD(T) calculations utilizing a wide range of split-valence, correlation-consistent, and bond-functions basis sets. Complete basis set MP2 (CBS MP2) stacking energies extrapolated using aug-cc-pVXZ (X = D, T, and for pyrimidine dimer Q) basis sets equal to −5.3, −12.3, and −11.2 kcal/mol for the first three dimers, respectively. Higher-order correlation corrections estimated as the difference between MP2 and CCSD(T) stacking energies amount to 2.0, 0.7, and 0.9 kcal/mol and lead to final estimates of the genuine stacking energies for the three dimers of −3.4, −11.6, and −10.4 kcal/mol. The CBS MP2 stacking-energy estimate for guanine dimer (−14.8 kcal/mol) was based on the 6-31G*(0.25) and aug-cc-pVDZ calculations. This simplified extrapolation can be routinely used with a meaningful accuracy around 1 kcal/mol for large aromatic stacking clusters. The final estimate of the guanine stacking energy after the CCSD(T) correction amounts to −12.9 kcal/mol. The MP2/6-31G*(0.25) method previously used as the standard level to calculate aromatic stacking in hundreds of geometries of nucleobase dimers systematically underestimates the base stacking by ca. 1.0−2.5 kcal/mol per stacked dimer, covering 75−90% of the intermolecular correlation stabilization. We suggest that this correction is to be considered in calibration of force fields and other cheaper computational methods. The quality of the MP2/6-31G*(0.25) predictions is nevertheless considerably better than suggested on the basis of monomer polarizability calculations. Fast and very accurate estimates of the MP2 aromatic stacking energies can be achieved using the RI-MP2 method. The CBS MP2 calculations and the CCSD(T) correction, when taken together, bring only marginal changes to the relative stability of H-bonded and stacked base pairs, with a slight shift of ca. 1 kcal/mol in favor of H-bonding. We suggest that the present values are very close to ultimate predictions of the strength of aromatic base stacking of DNA and RNA bases.