Improving intermolecular interactions in DFTB3 using extended polarization from chemical-potential equalization

Multipole Expansion of Atomic Electron Density Fluctuation Interactions in the Density-Functional Tight-Binding Method

The accuracy of the density-functional tight-binding (DFTB) method in describing noncovalent interactions is limited due to its reliance on monopole-based spherical charge densities. In this study, we present a multipole-extended second-order DFTB (mDFTB2) method that takes into account atomic dipole and quadrupole interactions. Furthermore, we combine the multipole expansion with the monopole-based third-order contribution, resulting in the mDFTB3 method. To assess the accuracy of mDFTB2 and mDFTB3, we evaluate their performance in describing noncovalent interactions, proton transfer barriers, and dipole moments. Our benchmark results show promising improvements even when using the existing electronic parameters optimized for the original DFTB3 model. Both mDFTB2 and mDFTB3 outperform their monopole-based counterparts, DFTB2 and DFTB3, in terms of accuracy. While mDFTB2 and mDFTB3 perform comparably for neutral and positively charged systems, mDFTB3 exhibits superior performance over mDFTB2 when dealing with negatively charged systems and proton transfers. Overall, the incorporation of the multipole expansion significantly enhances the accuracy of the DFTB method in describing noncovalent interactions and proton transfers.

Reparameterization of the chemical-potential equalization model with DFTB3: A practical balance between accuracy and transferability

To improve the performance of the third-order density-functional tight-binding method (DFTB3) for non-covalent interactions involving organic and biological molecules, a chemical-potential equalization (CPE) approach was introduced [J. Phys. Chem. A, 116, 9131 (2012)] and parameterized for the H, C, N, O, and S chemical elements [J. Chem. Phys., 143, 084123 (2015)]. Based largely on equilibrium structures, the parameterized DFTB3/CPE models were shown to exhibit improvements in molecular polarizabilities and intermolecular interactions. With more extensive analyses, however, we observe here that the available DFTB3/CPE models have two critical limitations: (1) they lead to sharply varying potential energy surfaces, thus causing numerical instability in molecular dynamics (MD) simulations, and (2) they lead to spurious interactions at short distances for some dimer complexes. These shortcomings are attributed to the employed screening functions and the overfitting of CPE parameters. In this work, we introduce a new strategy to simplify the parameterization procedure and significantly reduce free parameters down to four global (i.e., independent of element type) ones. With this strategy, two new models, DFTB3/CPE(r) and DFTB3/CPE(r†) are parameterized. The new models lead to smooth potential energy surfaces, stable MD simulations, and alleviate the spurious interactions at short distances, thus representing consistent improvements for both neutral and ionic hydrogen bonds.

Modeling pyranose ring pucker in carbohydrates using machine learning and semi-empirical quantum chemical methods

Pyranose ring pucker is a key coordinate governing the structure, interactions and reactivity of carbohydrates. We assess the ability of the machine learning potentials, ANI-1ccx and ANI-2x, and the GFN2-xTB semiempirical quantum chemical method, to model ring pucker conformers of five monosaccharides and oxane in the gas phase. Relative to coupled-cluster quantum mechanical calculations, we find that ANI-1ccx most accurately reproduces the ring pucker energy landscape for these molecules, with a correlation coefficient r² of 0.83. This correlation in relative energies lowers to values of 0.70 for ANI-2x and 0.60 for GFN2-xTB. The ANI-1ccx also provides the most accurate estimate of the energetics of the ⁴ C₁ -to-¹ C₄ minimum energy pathway for the six molecules. All three models reproduce chair more accurately than non-chair geometries. Analysis of small model molecules suggests that the ANI-1ccx model favors puckers with equatorial hydrogen bonding substituents; that ANI-2x and GFN2-xTB models overstabilize conformers with axially oriented groups; and that the endo-anomeric effect is overestimated by the machine learning models and underestimated via the GFN2-xTB method. While the pucker conformers considered in this study correspond to a gas phase environment, the accuracy and computational efficiency of the ANI-1ccx approach in modeling ring pucker in vacuo provides a promising basis for future evaluation and application to condensed phase environments.

A quantum chemical molecular dynamics repository of solvated ions

The importance of ion-solvent interactions in predicting specific ion effects in contexts ranging from viral activity through to electrolyte viscosity cannot be underestimated. Moreover, investigations of specific ion effects in nonaqueous systems, highly relevant to battery technologies, biochemical systems and colloid science, are severely limited by data deficiency. Here, we report IonSolvR – a collection of more than 3,000 distinct nanosecond-scale ab initio molecular dynamics simulations of ions in aqueous and non-aqueous solvent environments at varying effective concentrations. Density functional tight binding (DFTB) is used to detail the solvation structure of up to 55 solutes in 28 different protic and aprotic solvents. DFTB is a fast quantum chemical method, and as such enables us to bridge the gap between efficient computational scaling and maintaining accuracy, while using an internally-consistent simulation technique. We validate the database against experimental data and provide guidance for accessing individual IonSolvR records.

Measurement(s)	solvation structure
Technology Type(s)	quantum chemistry computational method • Molecular Dynamics

BENCHMARKS OF THE DENSITY FUNCTIONAL TIGHT-BINDING METHOD FOR REDOX, PROTONATION AND ELECTRONIC PROPERTIES OF QUINONES

Organic materials with controllable molecular design and sustainable resources are promising electrode materials. Crystalline quinones have been investigated in a variety of rechargeable battery chemistries due to their ubiquitous nature,...

Efficient calculation of protein-ligand binding free energy with GFN methods: the power of cluster model

The protein-ligand interactions are crucial in many biochemical processes and biomedical applications, yet it still remains challenging to accurately calculating the binding free energy of their interactions. In this work,...

Doubly Polarized QM/MM with Machine Learning Chaperone Polarizability

A major shortcoming of semiempirical (SE) molecular orbital methods is their severe underestimation of molecular polarizability compared with experimental and ab initio (AI) benchmark data. In a combined quantum mechanical and molecular mechanical (QM/MM) treatment of solution-phase reactions, solute described by SE methods therefore tends to generate inadequate electronic polarization response to solvent electric fields, which often leads to large errors in free energy profiles. To address this problem, here we present a hybrid framework that improves the response property of SE/MM methods through high-level molecular-polarizability fitting. Specifically, we place on QM atoms a set of corrective polarizabilities (referred to as chaperone polarizabilities), whose magnitudes are determined from machine learning (ML) to reproduce the condensed-phase AI molecular polarizability along the minimum free energy path. These chaperone polarizabilities are then used in a machinery similar to a polarizable force field calculation to compensate for the missing polarization energy in the conventional SE/MM simulations. Because QM atoms in this treatment host SE wave functions as well as classical polarizabilities, both polarized by MM electric fields, we name this method doubly polarized QM/MM (dp-QM/MM). We demonstrate the new method on the free energy simulations of the Menshutkin reaction in water. Using AM1/MM as a base method, we show that ML chaperones greatly reduce the error in the solute molecular polarizability from 6.78 to 0.03 ų with respect to the density functional theory benchmark. The chaperone correction leads to ∼10 kcal/mol of additional polarization energy in the product region, bringing the simulated free energy profiles to closer agreement with the experimental results. Furthermore, the solute-solvent radial distribution functions show that the chaperone polarizabilities modify the free energy profiles through enhanced solvation corrections when the system evolves from the charge-neutral reactant state to the charge-separated transition and product states. These results suggest that the dp-QM/MM method, enabled by ML chaperone polarizabilities, provides a very physical remedy for the underpolarization problem in SE/MM-based free energy simulations.

Interactions between large molecules pose a puzzle for reference quantum mechanical methods

Quantum-mechanical methods are used for understanding molecular interactions throughout the natural sciences. Quantum diffusion Monte Carlo (DMC) and coupled cluster with single, double, and perturbative triple excitations [CCSD(T)] are state-of-the-art trusted wavefunction methods that have been shown to yield accurate interaction energies for small organic molecules. These methods provide valuable reference information for widely-used semi-empirical and machine learning potentials, especially where experimental information is scarce. However, agreement for systems beyond small molecules is a crucial remaining milestone for cementing the benchmark accuracy of these methods. We show that CCSD(T) and DMC interaction energies are not consistent for a set of polarizable supramolecules. Whilst there is agreement for some of the complexes, in a few key systems disagreements of up to 8 kcal mol^-1 remain. These findings thus indicate that more caution is required when aiming at reproducible non-covalent interactions between extended molecules.

Quantum chemical benchmark databases of gold-standard dimer interaction energies

Advances in computational chemistry create an ongoing need for larger and higher-quality datasets that characterize noncovalent molecular interactions. We present three benchmark collections of quantum mechanical data, covering approximately 3,700 distinct types of interacting molecule pairs. The first collection, which we refer to as DES370K, contains interaction energies for more than 370,000 dimer geometries. These were computed using the coupled-cluster method with single, double, and perturbative triple excitations [CCSD(T)], which is widely regarded as the gold-standard method in electronic structure theory. Our second benchmark collection, a core representative subset of DES370K called DES15K, is intended for more computationally demanding applications of the data. Finally, DES5M, our third collection, comprises interaction energies for nearly 5,000,000 dimer geometries; these were calculated using SNS-MP2, a machine learning approach that provides results with accuracy comparable to that of our coupled-cluster training data. These datasets may prove useful in the development of density functionals, empirically corrected wavefunction-based approaches, semi-empirical methods, force fields, and models trained using machine learning methods.

Multi-level free energy simulation with a staged transformation approach

Combining multiple levels of theory in free energy simulations to balance computational accuracy and efficiency is a promising approach for studying processes in the condensed phase. While the basic idea has been proposed and explored for quite some time, it remains challenging to achieve convergence for such multi-level free energy simulations as it requires a favorable distribution overlap between different levels of theory. Previous efforts focused on improving the distribution overlap by either altering the low-level of theory for the specific system of interest or ignoring certain degrees of freedom. Here, we propose an alternative strategy that first identifies the degrees of freedom that lead to gaps in the distributions of different levels of theory and then treats them separately with either constraints or restraints or by introducing an intermediate model that better connects the low and high levels of theory. As a result, the conversion from the low level to the high level model is done in a staged fashion that ensures a favorable distribution overlap along the way. Free energy components associated with different steps are mostly evaluated explicitly, and thus, the final result can be meaningfully compared to the rigorous free energy difference between the two levels of theory with limited and well-defined approximations. The additional free energy component calculations involve simulations at the low level of theory and therefore do not incur high computational costs. The approach is illustrated with two simple but non-trivial solution examples, and factors that dictate the reliability of the result are discussed.

Density-functional tight-binding: basic concepts and applications to molecules and clusters

The scope of this article is to present an overview of the Density Functional based Tight Binding (DFTB) method and its applications. The paper introduces the basics of DFTB and its standard formulation up to second order. It also addresses methodological developments such as third order expansion, inclusion of non-covalent interactions, schemes to solve the self-interaction error, implementation of long-range short-range separation, treatment of excited states via the time-dependent DFTB scheme, inclusion of DFTB in hybrid high-level/low level schemes (DFT/DFTB or DFTB/MM), fragment decomposition of large systems, large scale potential energy landscape exploration with molecular dynamics in ground or excited states, non-adiabatic dynamics. A number of applications are reviewed, focusing on -(i)- the variety of systems that have been studied such as small molecules, large molecules and biomolecules, bare orfunctionalized clusters, supported or embedded systems, and -(ii)- properties and processes, such as vibrational spectroscopy, collisions, fragmentation, thermodynamics or non-adiabatic dynamics. Finally outlines and perspectives are given.

What Next for Quantum Mechanics in Structure-Based Drug Discovery?

There is significant potential for electronic structure methods to improve the quality of the predictions furnished by the tools of computer-aided drug design, which typically rely on empirically derived functions. In this perspective, we consider some recent examples of how quantum mechanics has been applied in predicting protein-ligand geometries, protein-ligand binding affinities and ligand strain on binding. We then outline several significant developments in quantum mechanics methodology likely to influence these approaches: in particular, we note the advent of more computationally expedient ab initio quantum mechanical methods that can provide chemical accuracy for larger molecular systems than hitherto possible. We highlight the emergence of increasingly accurate semiempirical quantum mechanical methods and the associated role of machine learning and molecular databases in their development. Indeed, the convergence of improved algorithms for solving and analyzing electronic structure, modern machine learning methods, and increasingly comprehensive benchmark data sets of molecular geometries and energies provides a context in which the potential of quantum mechanics will be increasingly realized in driving future developments and applications in structure-based drug discovery.

Analysis of Density Functional Tight Binding with Natural Bonding Orbitals

The description of chemical bonding by the density functional tight binding (DFTB) model is analyzed using natural bonding orbitals (NBOs) and compared to results from density functional theory (B3LYP/aug-cc-pVTZ) calculations. Several molecular systems have been chosen to represent fairly diverse bonding scenarios that include standard covalent bonds, hypervalent interactions, multicenter bonds, metal-ligand interactions (with and without the pseudo-Jahn-Teller effect), and through-space donor-acceptor interactions. Overall, the results suggest that DFTB3/3OB provides physically sound descriptions for the different bonding scenarios analyzed here, as reflected by the general agreement between DFTB3 and B3LYP NBO properties, such as the nature of the NBOs, the magnitudes of natural charges and bond orders, and the dominant donor-acceptor interactions. The degree of ligand-to-metal charge transfer and the ionic nature of pentavalent phosphate are overestimated, likely reflecting the minimal-basis nature of DFTB3/3OB. Moreover, certain orbital interactions, such as geminal interactions, are observed to be grossly overestimated by DFTB3 for hypervalent phosphate and several transition metal compounds that involve copper and nickel. The study indicates that results from NBO analysis can be instructive for identifying electronic structure descriptions at the approximate quantum-mechanical level that require improvement and thus for guiding the systematic improvement of these methods.

GFN2-xTB-An Accurate and Broadly Parametrized Self-Consistent Tight-Binding Quantum Chemical Method with Multipole Electrostatics and Density-Dependent Dispersion Contributions

An extended semiempirical tight-binding model is presented, which is primarily designed for the fast calculation of structures and noncovalent interaction energies for molecular systems with roughly 1000 atoms. The essential novelty in this so-called GFN2-xTB method is the inclusion of anisotropic second order density fluctuation effects via short-range damped interactions of cumulative atomic multipole moments. Without noticeable increase in the computational demands, this results in a less empirical and overall more physically sound method, which does not require any classical halogen or hydrogen bonding corrections and which relies solely on global and element-specific parameters (available up to radon, Z = 86). Moreover, the atomic partial charge dependent D4 London dispersion model is incorporated self-consistently, which can be naturally obtained in a tight-binding picture from second order density fluctuations. Fully analytical and numerically precise gradients (nuclear forces) are implemented. The accuracy of the method is benchmarked for a wide variety of systems and compared with other semiempirical methods. Along with excellent performance for the "target" properties, we also find lower errors for "off-target" properties such as barrier heights and molecular dipole moments. High computational efficiency along with the improved physics compared to its precursor GFN-xTB makes this method well-suited to explore the conformational space of molecular systems. Significant improvements are furthermore observed for various benchmark sets, which are prototypical for biomolecular systems in aqueous solution.

Exploring the applicability of density functional tight binding to transition metal ions. Parameterization for nickel with the spin-polarized DFTB3 model

In this work, we explore the applicability and limitations of the current third order density functional tight binding (DFTB3) formalism for treating transition metal ions using nickel as an example. To be consistent with recent parameterization of DFTB3 for copper, the parametrization for nickel is conducted in a spin-polarized formulation and with orbital-resolved Hubbard parameters and their charge derivatives. The performance of the current parameter set is evaluated based on structural and energetic properties of a set of nickel-containing compounds that involve biologically relevant ligands. Qualitatively similar to findings in previous studies of copper complexes, the DFTB3 results are more reliable for nickel complexes with neutral ligands than for charged ligands; nevertheless, encouraging agreement is noted in comparison to the reference method, B3LYP/aug-cc-pVTZ, especially for structural properties, including cases that exhibit Jahn-Teller distortions; the structures also compare favorably to available X-ray data in the Cambridge Crystallographic Database for a number of nickel-containing compounds. As to limitations, we find it is necessary to use different d shell Hubbard charge derivatives for Ni(I) and Ni(II), due to the distinct electronic configurations for the nickel ion in the respective complexes, and substantial errors are observed for ligand binding energies, especially for charged ligands, d orbital splitting energies and splitting between singlet and triplet spin states for Ni(II) compounds. These observations highlight that future improvement in intra-d correlation and ligand polarization is required to enable the application of the DFTB3 model to complex transition metal ions. © 2018 Wiley Periodicals, Inc.

Alchemical and structural distribution based representation for universal quantum machine learning

We introduce a representation of any atom in any chemical environment for the automatized generation of universal kernel ridge regression-based quantum machine learning (QML) models of electronic properties, trained throughout chemical compound space. The representation is based on Gaussian distribution functions, scaled by power laws and explicitly accounting for structural as well as elemental degrees of freedom. The elemental components help us to lower the QML model's learning curve, and, through interpolation across the periodic table, even enable "alchemical extrapolation" to covalent bonding between elements not part of training. This point is demonstrated for the prediction of covalent binding in single, double, and triple bonds among main-group elements as well as for atomization energies in organic molecules. We present numerical evidence that resulting QML energy models, after training on a few thousand random training instances, reach chemical accuracy for out-of-sample compounds. Compound datasets studied include thousands of structurally and compositionally diverse organic molecules, non-covalently bonded protein side-chains, (H₂O)₄₀-clusters, and crystalline solids. Learning curves for QML models also indicate competitive predictive power for various other electronic ground state properties of organic molecules, calculated with hybrid density functional theory, including polarizability, heat-capacity, HOMO-LUMO eigenvalues and gap, zero point vibrational energy, dipole moment, and highest vibrational fundamental frequency.

Intermolecular interactions in the condensed phase: Evaluation of semi-empirical quantum mechanical methods

To facilitate further development of approximate quantum mechanical methods for condensed phase applications, we present a new benchmark dataset of intermolecular interaction energies in the solution phase for a set of 15 dimers, each containing one charged monomer. The reference interaction energy in solution is computed via a thermodynamic cycle that integrates dimer binding energy in the gas phase at the coupled cluster level and solute-solvent interaction with density functional theory; the estimated uncertainty of such calculated interaction energy is ±1.5 kcal/mol. The dataset is used to benchmark the performance of a set of semi-empirical quantum mechanical (SQM) methods that include DFTB3-D3, DFTB3/CPE-D3, OM2-D3, PM6-D3, PM6-D3H+, and PM7 as well as the HF-3c method. We find that while all tested SQM methods tend to underestimate binding energies in the gas phase with a root-mean-squared error (RMSE) of 2-5 kcal/mol, they overestimate binding energies in the solution phase with an RMSE of 3-4 kcal/mol, with the exception of DFTB3/CPE-D3 and OM2-D3, for which the systematic deviation is less pronounced. In addition, we find that HF-3c systematically overestimates binding energies in both gas and solution phases. As most approximate QM methods are parametrized and evaluated using data measured or calculated in the gas phase, the dataset represents an important first step toward calibrating QM based methods for application in the condensed phase where polarization and exchange repulsion need to be treated in a balanced fashion.

Parametrization and Benchmark of Long-Range Corrected DFTB2 for Organic Molecules

We present the parametrization and benchmark of long-range corrected second-order density functional tight binding (DFTB), LC-DFTB2, for organic and biological molecules. The LC-DFTB2 model not only improves fundamental orbital energy gaps but also ameliorates the DFT self-interaction error and overpolarization problem, and further improves charge-transfer excited states significantly. Electronic parameters for the construction of the DFTB2 Hamiltonian as well as repulsive potentials were optimized for molecules containing C, H, N, and O chemical elements. We use a semiautomatic parametrization scheme based on a genetic algorithm. With the new parameters, LC-DFTB2 describes geometries and vibrational frequencies of organic molecules similarly well as third-order DFTB3/3OB, the de facto standard parametrization based on a GGA functional. LC-DFTB2 performs well also for atomization and reaction energies, however, slightly less satisfactorily than DFTB3/3OB.

Quantum mechanical force fields for condensed phase molecular simulations

Molecular simulations are powerful tools for providing atomic-level details into complex chemical and physical processes that occur in the condensed phase. For strongly interacting systems where quantum many-body effects are known to play an important role, density-functional methods are often used to provide the model with the potential energy used to drive dynamics. These methods, however, suffer from two major drawbacks. First, they are often too computationally intensive to practically apply to large systems over long time scales, limiting their scope of application. Second, there remain challenges for these models to obtain the necessary level of accuracy for weak non-bonded interactions to obtain quantitative accuracy for a wide range of condensed phase properties. Quantum mechanical force fields (QMFFs) provide a potential solution to both of these limitations. In this review, we address recent advances in the development of QMFFs for condensed phase simulations. In particular, we examine the development of QMFF models using both approximate and ab initio density-functional models, the treatment of short-ranged non-bonded and long-ranged electrostatic interactions, and stability issues in molecular dynamics calculations. Example calculations are provided for crystalline systems, liquid water, and ionic liquids. We conclude with a perspective for emerging challenges and future research directions.

Perspective: Quantum mechanical methods in biochemistry and biophysics

In this perspective article, I discuss several research topics relevant to quantum mechanical (QM) methods in biophysical and biochemical applications. Due to the immense complexity of biological problems, the key is to develop methods that are able to strike the proper balance of computational efficiency and accuracy for the problem of interest. Therefore, in addition to the development of novel ab initio and density functional theory based QM methods for the study of reactive events that involve complex motifs such as transition metal clusters in metalloenzymes, it is equally important to develop inexpensive QM methods and advanced classical or quantal force fields to describe different physicochemical properties of biomolecules and their behaviors in complex environments. Maintaining a solid connection of these more approximate methods with rigorous QM methods is essential to their transferability and robustness. Comparison to diverse experimental observables helps validate computational models and mechanistic hypotheses as well as driving further development of computational methodologies.

Metal Ion Modeling Using Classical Mechanics

Metal ions play significant roles in numerous fields including chemistry, geochemistry, biochemistry, and materials science. With computational tools increasingly becoming important in chemical research, methods have emerged to effectively face the challenge of modeling metal ions in the gas, aqueous, and solid phases. Herein, we review both quantum and classical modeling strategies for metal ion-containing systems that have been developed over the past few decades. This Review focuses on classical metal ion modeling based on unpolarized models (including the nonbonded, bonded, cationic dummy atom, and combined models), polarizable models (e.g., the fluctuating charge, Drude oscillator, and the induced dipole models), the angular overlap model, and valence bond-based models. Quantum mechanical studies of metal ion-containing systems at the semiempirical, ab initio, and density functional levels of theory are reviewed as well with a particular focus on how these methods inform classical modeling efforts. Finally, conclusions and future prospects and directions are offered that will further enhance the classical modeling of metal ion-containing systems.

Challenging Dogmas: Hydrogen Bond Revisited

Hydrogen bond directionality in the water dimer is explained on the basis of symmetry-adapted intermolecular perturbation theory which directly separates the intermolecular interaction energy into four physically interpretable components: electrostatics, exchange-repulsion, dispersion, and induction. Analysis of these four main contributions to the binding energy allows a deeper understanding of the dominant factors ruling the mutual arrangement of the two monomers. A preference for the linear configuration is shown to be due to a subtle interplay of all four energy components. While the first-order terms, electrostatic and exchange-repulsion, almost perfectly cancel each other near the equilibrium geometry of the dimer, the importance of the second- and higher-order terms, induction and dispersion, becomes evident.

Semiempirical Quantum Mechanical Methods for Noncovalent Interactions for Chemical and Biochemical Applications

Semiempirical (SE) methods can be derived from either Hartree-Fock or density functional theory by applying systematic approximations, leading to efficient computational schemes that are several orders of magnitude faster than ab initio calculations. Such numerical efficiency, in combination with modern computational facilities and linear scaling algorithms, allows application of SE methods to very large molecular systems with extensive conformational sampling. To reliably model the structure, dynamics, and reactivity of biological and other soft matter systems, however, good accuracy for the description of noncovalent interactions is required. In this review, we analyze popular SE approaches in terms of their ability to model noncovalent interactions, especially in the context of describing biomolecules, water solution, and organic materials. We discuss the most significant errors and proposed correction schemes, and we review their performance using standard test sets of molecular systems for quantum chemical methods and several recent applications. The general goal is to highlight both the value and limitations of SE methods and stimulate further developments that allow them to effectively complement ab initio methods in the analysis of complex molecular systems.

Towards a barrier height benchmark set for biologically relevant systems

We have collected computed barrier heights and reaction energies (and associated model structures) for five enzymes from studies published by Himo and co-workers. Using this data, obtained at the B3LYP/6- 311+G(2d,2p)[LANL2DZ]//B3LYP/6-31G(d,p) level of theory, we then benchmark PM6, PM7, PM7-TS, and DFTB3 and discuss the influence of system size, bulk solvation, and geometry re-optimization on the error. The mean absolute differences (MADs) observed for these five enzyme model systems are similar to those observed for PM6 and PM7 for smaller systems (10-15 kcal/mol), while DFTB results in a MAD that is significantly lower (6 kcal/mol). The MADs for PMx and DFTB3 are each dominated by large errors for a single system and if the system is disregarded the MADs fall to 4-5 kcal/mol. Overall, results for the condensed phase are neither more or less accurate relative to B3LYP than those in the gas phase. With the exception of PM7-TS, the MAD for small and large structural models are very similar, with a maximum deviation of 3 kcal/mol for PM6. Geometry optimization with PM6 shows that for one system this method predicts a different mechanism compared to B3LYP/6-31G(d,p). For the remaining systems, geometry optimization of the large structural model increases the MAD relative to single points, by 2.5 and 1.8 kcal/mol for barriers and reaction energies. For the small structural model, the corresponding MADs decrease by 0.4 and 1.2 kcal/mol, respectively. However, despite these small changes, significant changes in the structures are observed for some systems, such as proton transfer and hydrogen bonding rearrangements. The paper represents the first step in the process of creating a benchmark set of barriers computed for systems that are relatively large and representative of enzymatic reactions, a considerable challenge for any one research group but possible through a concerted effort by the community. We end by outlining steps needed to expand and improve the data set and how other researchers can contribute to the process.

Copper Oxidation/Reduction in Water and Protein: Studies with DFTB3/MM and VALBOND Molecular Dynamics Simulations

We apply two recently developed computational methods, DFTB3 and VALBOND, to study copper oxidation/reduction processes in solution and protein. The properties of interest include the coordination structure of copper in different oxidation states in water or in a protein (plastocyanin) active site, the reduction potential of the copper ion in different environments, and the environmental response to copper oxidation. The DFTB3/MM and VALBOND simulation results are compared to DFT/MM simulations and experimental results whenever possible. For a copper ion in aqueous solution, DFTB3/MM results are generally close to B3LYP/MM with a medium basis, including both solvation structure and reduction potential for Cu(II); for Cu(I), however, DFTB3/MM finds a two-water coordination, similar to previous Born-Oppenheimer molecular dynamics simulations using BLYP and HSE, whereas B3LYP/MM leads to a tetrahedron coordination. For a tetraammonia copper complex in aqueous solution, VALBOND and DFTB3/MM are consistent in terms of both structural and dynamical properties of solvent near copper for both oxidation states. For copper reduction in plastocyanin, DFTB3/MM simulations capture the key properties of the active site, and the computed reduction potential and reorganization energy are in fair agreement with experiment, especially when the periodic boundary condition is used. Overall, the study supports the value of VALBOND and DFTB3(/MM) for the analysis of fundamental copper redox chemistry in water and protein, and the results also help highlight areas where further improvements in these methods are desirable.

DFTB3 Parametrization for Copper: The Importance of Orbital Angular Momentum Dependence of Hubbard Parameters

We report the parametrization of a density functional tight binding method (DFTB3) for copper in a spin-polarized formulation. The parametrization is consistent with the framework of 3OB for main group elements (ONCHPS) and can be readily used for biological applications that involve copper proteins/peptides. The key to our parametrization is to introduce orbital angular momentum dependence of the Hubbard parameter and its charge derivative, thus allowing the 3d and 4s orbitals to adopt different sizes and responses to the change of charge state. The parametrization has been tested by applying to a fairly broad set of molecules of biological relevance, and the properties of interest include optimized geometries, ligand binding energies, and ligand proton affinities. Compared to the reference QM level (B3LYP/aug-cc-pVTZ, which is shown here to be similar to the B97-1 and CCSD(T) results, in terms of many properties of interest for a set of small copper containing molecules), our parametrization generally gives reliable structural properties for both Cu(I) and Cu(II) compounds, although several exceptions are also noted. For energetics, the results are more accurate for neutral ligands than for charged ligands, likely reflecting the minimal basis limitation of DFTB3; the results generally outperform NDDO based methods such as PM6 and even PBE with the 6-31+G(d,p) basis. For all ligand types, single-point B3LYP calculations at DFTB3 geometries give results very close (∼1-2 kcal/mol) to the reference B3LYP values, highlighting the consistency between DFTB3 and B3LYP structures. Possible further developments of the DFTB3 model for a better treatment of transition-metal ions are also discussed. In the current form, our first generation of DFTB3 copper model is expected to be particularly valuable as a method that drives sampling in systems that feature a dynamical copper binding site.