Parametrization of DFTB3/3OB for magnesium and zinc for chemical and biological applications.

We report the parametrization of the approximate density functional theory, DFTB3, for magnesium and zinc for chemical and biological applications. The parametrization strategy follows that established in previous work that parametrized several key main group elements (O, N, C, H, P, and S). This 3OB set of parameters can thus be used to study many chemical and biochemical systems. The parameters are benchmarked using both gas-phase and condensed-phase systems. The gas-phase results are compared to DFT (mostly B3LYP), ab initio (MP2 and G3B3), and PM6, as well as to a previous DFTB parametrization (MIO). The results indicate that DFTB3/3OB is particularly successful at predicting structures, including rather complex dinuclear metalloenzyme active sites, while being semiquantitative (with a typical mean absolute deviation (MAD) of ∼3-5 kcal/mol) for energetics. Single-point calculations with high-level quantum mechanics (QM) methods generally lead to very satisfying (a typical MAD of ∼1 kcal/mol) energetic properties. DFTB3/MM simulations for solution and two enzyme systems also lead to encouraging structural and energetic properties in comparison to available experimental data. The remaining limitations of DFTB3, such as the treatment of interaction between metal ions and highly charged/polarizable ligands, are also discussed.

Introduction

class="Chemical">Magnesium and zinc are two of the most class="Chemical">pan class="Chemical">common metal cofactors of numerous metalloenzymes and other biomolecules.[1−3] The binding of Mg2+ to DNA and RNA helps to maintain their structures by stabilizing highly negatively charged phosphates. ATP, the main source of energy in cells, is bound to Mg2+ in its biologically active form.[4] In addition, a large number of enzymes involved in the biochemistry of nucleic acids bind Mg2+ for catalytic activity; for example, Mg2+ is essential for DNA replication and phosphoryl transfers.[5] Zinc is the second most abundant cation in the human body, so Zn2+ has important biological functions in numerous biomolecular systems[6,7] such as zinc-finger class of transcription factors, signaling proteins, transport/storage proteins, as well as enzymes. The zinc ion contributes to catalysis in more than 300 enzymes and serves a structural role in many proteins and enzymes.[8−10]

To describe the structural properties of class="Chemical">Mg2+ and class="Chemical">pan class="Chemical">Zn2+ in solution and biomolecules, molecular mechanical (MM) force field parameters have been developed.[11−15] Some of the advanced force field models such as SIBFA,[16] AMOEBA,[17,18] and Drude oscillator models[19] treat electronic polarization and charge transfer[20] terms explicitly and have been demonstrated to describe potential energy surfaces in close agreement with quantum mechanics (QM) calculations.[21,22] For chemical reactions that involve these metal ions, however, MM models are no longer appropriate. Even for structural properties, when the binding mode of ligands is flexible,[23−25] it is challenging to describe the metal site with a MM model.

An appropriate modeling approach in tclass="Chemical">his class="Chemical">pan class="Chemical">context is the QM/MM framework, where the metal site of interest is treated with QM while the rest of the enzyme/solution is described by an empirical force field. For relatively rigid metal sites, ab initio QM and QM/MM methods based on either reaction path or free energy simulations have been demonstrated to be powerful for mechanistic analysis.[26−30] The cost of these calculations for flexible systems, however, remains very high. For instance, one key feature of zinc coordination is its flexibility:[23−25] a zinc ion can adopt multiple binding modes with the coordination number ranging from 4 to 6, underlining the importance of developing efficient QM(/MM) methods that complement ab initio based methods. The dynamical nature of the metal binding mode is particularly important to enzymes that feature solvent accessible active sites.[31−33] For these applications, semiempirical QM-based QM/MM simulations can potentially strike the balance between computational accuracy and conformational sampling.

The traditional semiempirical QM methods based on the Neglect-of-Diatomic-Differential-Overlap (class="Chemical">NDDO) aclass="Chemical">pclass="Chemical">proximation, such as MNDO,[34,35] MNDO/d,[36] class="Chemical">pan class="Species">AM1,[37] and PM3,[38] however, give in general rather poor geometries for Mg/Zn-containing molecules. For example, metal–ligand lengths computed from AM1,[39,40] PM3, MNDO, and MNDO/d[40] deviate by roughly 0.1 Å on average in comparison to B3LYP/6-311++G** for small Mg complexes.[41] Recently developed AM1/d parameters[42] for magnesium provide a clear improvement in accuracy compared to AM1 by incorporating d orbitals. According to a recent test,[43] the mean absolute deviations (MADs) of metal–ligand distances for AM1,[44] PM3, PM3(ZnB),[45] MNDO/d,[46] and PM6[47] are more than 0.05 Å compared to CCSD(T) values for small zinc-containing compounds. The importance of a reliable treatment for zinc sites to the description of enzyme catalysis is illustrated by our recent study of alkaline phosphatase (AP),[48,59] where we showed that previous AM1/d-PhoT/MM calculations[50,51] likely led to incorrect transition states due to distortions of the bimetallic zinc motif.

The Self-class="Chemical">Consistent-Chclass="Chemical">pan class="Chemical">arge Density-Functional Tight-Binding (SCC-DFTB) method provides a promising alternative to NDDO approaches.[52−54] It is derived by expanding the DFT total energy up to second order around a reference density, which is usually taken to be the sum of atomic charge densities. The computational speed of SCC-DFTB is comparable to NDDO, and it has been rather extensively benchmarked for reaction energies, geometries, and vibrational frequencies for a large number of small molecules.[55−61] The fact that SCC-DFTB is formulated in a DFT framework suggests that it is potentially more reliable for metal ions than NDDO-based methods, which are based on Hartree–Fock. Specifically for Mg[41] and Zn,[62,63] SCC-DFTB(/MM)[64−66] has been successfully applied to a rather diverse set of biological and chemical applications that include these ions.[48,49,67−74] More recently, DFTB has been extended to include third-order contributions, termed as DFTB3.[75] Compared to the second-order formulation, DFTB3 considers the charge dependence of the atomic Hubbard parameter and therefore provides an improved description of chemical properties such as proton affinity,[53,75,76] which is important in many biological applications. With the most recent parametrization referred to as 3OB,[77,78] the electronic properties and therefore atomization energies have been substantially improved; as a result, the DFTB3 approach is shown to have comparable performance as DFT with double-ζ-plus-polarization basis sets for many systems of biological and chemical interest, although certain limitations remain.[78,79]

pan class="Chemical">Considering the recent success of class="Chemical">pan class="Chemical">DFTB3/3OB, here we extend the 3OB parametrization of DFTB3 to magnesium and zinc, again focusing on chemical and biological applications. In the next section, we briefly summarize the DFTB3 methodology, the procedures for parametrization, and simulation details for benchmark calculations. We then present data from multiple levels of test sets in the gas phase, which are followed by additional benchmark calculations in both solution and enzymes. Finally, we summarize this work with several conclusions.

Computational Methods

Theory

The pan class="Chemical">DFTB3 class="Chemical">parametrization of class="Chemical">pan class="Chemical">magnesium and zinc follows the protocol outlined in the previous publications[77,80] and extends the 3OB parameter set that so far includes elements C, H, N, O,[77] S, and P.[78] We first briefly review the formulation of DFTB3.[54,75]

The total energy of pan class="Chemical">DFTB3 is given by

Here H0 is the effective Kohn–Sham Hamiltonian evaluated at the atomic reference density ρ0, which is determined by solving the Kohn–Sham equations for the atom in the presence of an additional harmonic class="Chemical">confining class="Chemical">potential with a class="Chemical">pan class="Chemical">confining radius rwf for the atomic basis set and a different radius rdens for the density. Off-diagonal H0 integrals are computed in a two-center approximation together with the PBE[81] exchange-correlation functional and an atomic basis set consisting of Slater functions. The Slater functions are defined by lmax as the highest orbital angular momentum included, nmax, and α0, α1, ..., α4, the exponents of the Slater functions. In principle, basis functions for each orbital angular momentum can be compressed differently; in this work, however, s and p orbitals are confined with the same compression radius for Mg, and s, p, and d orbitals (compression radius for d-orbitals is defined by the parameter rwfd) are confined with the same compression radius for Zn. The diagonal elements H0 are equal to the calculated atomic eigenvalues ε (x = s, p, or d orbital). No confinement of the orbitals is applied here in order to yield the correct dissociation limit. We optimize εp as discussed in the following section.

The atomic Hubbard parameter Ua in the seclass="Chemical">cond-order term of the energy Eγ describes the electron–electron interaction within one atom and enters the γ-function, which interclass="Chemical">polates the on-site and the long-range electronic interaction of atom class="Chemical">pairs. The Hubbard class="Chemical">parameters are calculated from DFT as the first derivative of the highest occuclass="Chemical">pied orbital with resclass="Chemical">pect to its occuclass="Chemical">pation number. For the third-order term EΓ, an additional class="Chemical">parameter class="Chemical">per element Ud, the derivative of the Hubbard class="Chemical">parameter with resclass="Chemical">pect to chclass="Chemical">pan class="Chemical">arge has to be determined. For CHNO it was calculated as the second derivative of the highest occupied orbital with respect to its occupation number;[77] for Mg and Zn we find a significant advantage to optimize it.

Reference Systems and Parameter Fitting

As emphasized earlier,[54] most parameters in the class="Chemical">DFTB3 aclass="Chemical">pclass="Chemical">proach are directly class="Chemical">pan class="Chemical">computed based on atoms or atom pairs at the DFT (PBE) level, while some parameters are optimized based on specific reference systems. In this way, accuracy appropriate for chemical applications can be obtained without significantly reducing the transferability of the DFTB3 model. The parametrization for Mg/Zn used herein follows that specified previously for the parametrization of CHNOPS;[77,78] two distinct groups of parameters have to be optimized, the electronic parameters appearing in the first three terms of the DFTB3 total energy (eq 1) and the repulsive potential.

The electronic parameters are summarized in Table 1. The following parameters are optimized to improve the overall performance:

Table 1

Overview of Electronic Parameters (in Atomic Units if Not Unitless)a

parameter	Mg	Zn
lmax	1	2
nmax	2	2
α0	0.50	0.50
α1	1.11	1.39
α2	2.45	3.87
α3	5.42	10.78
α4	12.00	30.00
rwf	5.5	4.2
rwfd	0.0	4.2
rdens	14.0	9.0
εs	–0.17267741	–0.21408495
εp	–0.02	0.02
εd	–	–0.38831893
Espin	0.0	0.0
U	0.2246	0.2663
Ud	–0.02	–0.03

As described in the text, U, Espin, εs,d, nmax, and α are consistent with the standard choices of DFTB parametrization and not subject to optimizations. By contrast, rwf, rwfd, rdens, εp, and Ud are adjusted based on properties of molecules in the fitting set.

•rwf and rdens: The general procedure was to start from the electronic parameters used in the class="Chemical">MIO set[75,77] (which are rclass="Chemical">pan class="Chemical">Mgwf = 5.0, rMgdens = 12.0, rZnwf = 4.9, rZndens = 10.0). rwf and rdens are fitted in a “brute force” manner to achieve good performance on bond (metal–ligand), angle, vibrational frequencies, as well as reaction energies.

•εp: The eigenvalues are in principle class="Chemical">comclass="Chemical">puted directly from atomic DFT calculations and enter as diagonal elements into the Hamiltonian matrix without the additional class="Chemical">pan class="Chemical">confining potential. Exceptions are made here to magnesium and zinc since the p orbitals are often not heavily involved in the coordination of Mg2+ and Zn2+. εp is calculated to be −0.04914 au for Mg and −0.04261 au for Zn. εp are fitted in the way that p orbitals are less involved in molecular bonding, so as to improve sequential ligation energies and proton affinities.

•Ud: The Hubbard derivative is generally class="Chemical">comclass="Chemical">puted from atomic DFT calculations. In the case of class="Chemical">pan class="Chemical">Mg and Zn, they are −0.0496 au and −0.0544 au, respectively. We found a significant improvement of proton affinities after optimizing Ud for Mg and Zn. It is important to note that fitting of the Hubbard derivatives does not significantly alter most properties of neutral molecules like equilibrium geometries, so with all the other parameters held fixed, Ud is tuned as the last step in the parametrization procedure.

Although the repulsive potential, which class="Chemical">contains nuclear–nuclear reclass="Chemical">pulsion and double-class="Chemical">pan class="Chemical">counting terms that depend on the reference density, can be directly computed, practical application to molecular systems with reasonable accuracy requires fitting the repulsive potential. The repulsive potential between atom type A and type B (Vabrep in eq 1) is described by a spline function in the covalent binding region. A cutoff is introduced at which the function and first three derivatives reach zero. The remaining degrees of freedom for the spline function are fitted to high-level ab initio values of geometries, reaction energies, and vibrational stretching frequencies. The overall quality of the repulsive potential is sensitive to the chosen division points. Therefore, establishing reliable repulsive potential parameters remains a time-consuming step despite the progress of partially automatized fitting procedures.[80,82]

As described in the text, U, Espin, εs,d, nmax, and α are class="Chemical">consistent with the standard choices of class="Chemical">pan class="Chemical">DFTB parametrization and not subject to optimizations. By contrast, rwf, rwfd, rdens, εp, and Ud are adjusted based on properties of molecules in the fitting set.

The repulsive potentials are fitted to several properties as summarized in Table 2. Reference values are calculated from ab initio methods that are known to produce results close to experimental data. For class="Chemical">magnesium, oclass="Chemical">ptimized geometries of small closed-shell molecules are calculated from B3LYP[83−85]/aug-cc-class="Chemical">pVTZ, and reaction energies are calculated from the G3B3[86−88] method. For zinc, both oclass="Chemical">ptimized geometries and reaction energies are calculated from B3LYP/aug-cc-class="Chemical">pVTZ. Additional equations for the fitting class="Chemical">process are class="Chemical">preclass="Chemical">pared using a few vibrational stretching frequencies determined from BLYP/aug-cc-class="Chemical">pVTZ (unscaled) using the harmonic aclass="Chemical">pclass="Chemical">proximation. B3LYP with good quality basis set has been shown to describe structures of class="Chemical">pan class="Chemical">Zn complexes in good agreement with X-ray[89] and CCSD(T)[43,90] data. However, heats of formation (HOF) predicted from B3LYP/aug-cc-pVTZ deviate ∼10 kcal/mol on average from experimental data of eight small complexes even though CCSD(T)/aug-cc-pVTZ does not do much better with a mean unsigned error of ∼7 kcal/mol.[91] B97-1,[92] B97-2,[93] and B98[94,95] have also been tested because they have been reported to be reliable for reproducing experimental data for small transition metal complexes;[96] geometries from the B97 family functionals and B3LYP are very similar, and reaction energies differ by ∼2 kcal/mol (data shown in the Supporting Information). Despite shortcomings in HOF, we have chosen to use the B3LYP functional because it is computationally less intensive so that large sets of consistent reference data are accessible, and it systematically gives reliable results for binding energies and proton affinities, which are of major importance to our work. Since the systems studied here are relatively small in size, dispersion interactions are not expected to impact the structural and energetic properties of interest. For a fairly large active site model for a zinc enzyme, we do discuss briefly the impact of dispersion. For even larger cases, we note that empirical dispersion models have been developed and tested for both SCC-DFTB[97,98] and DFTB3.[99]

Table 2

Parameters Defining the Repulsive Potentialsa

reference geometries for force equations
MgH2, Mg(CH3)2, MgO, [Mg(H2O)2]2+, [Mg(H2O)6]2+, [Mg(NH3)4]2+, MgS, [Mg(SH2)3(SH)]+, Mg(PH2)2, [Mg(PH3)4]2+, ZnH2, [Zn(CH3)]+, [Zn(CO)]2+, [Zn(H2O)3]2+, Zn(OH)2, Zn(H)(NH2), [Zn(NH3)4]2+, [Zn(SH)]+, [Zn(SH2)3]2+, [Zn(SH2)3(SH)]+, [Zn(PH2)]+, [Zn(PH3)4]2+, [Zn(PH3)3(PH2)]+
reference potential energy differences in kcal/mol
[Mg(H2O)2]2+	→	[Mg(H2O)]2+ + H2O	73.3
[Mg(H2O)3]2+	→	[Mg(H2O)2]2+ + H2O	59.2
[Mg(CO)4]2+	→	[Mg(CO)3]2+ + CO	30.9
[Mg(NH3)]2+	→	Mg2+ + NH3	97.7
[Mg(NH3)2]2+	→	[Mg(NH3)]2+ + NH3	83.9
[Mg(NH3)3]2+	→	[Mg(NH3)2]2+ + NH3	63.1
[Mg(NH3)4]2+	→	[Mg(NH3)3]2+ + NH3	50.1
[Mg(NH3)5]2+	→	[Mg(NH3)4]2+ + NH3	28.6
[Mg(NH3)6]2+	→	[Mg(NH3)5]2+ + NH3	26.3
[Mg(SH2)2]2+	→	[Mg(SH2)]2+ + SH2	62.4
[Mg(SH2)3]2+	→	[Mg(SH2)2]2+ + SH2	43.9
[Mg(SH2)4]2+	→	[Mg(SH2)3]2+ + SH2	34.5
[Mg(SH2)5]2+	→	[Mg(SH2)4]2+ + SH2	20.0
[Mg(PH3)2]2+	→	[Mg(PH3)]2+ + PH3	66.9
[Mg(PH3)4]2+	→	[Mg(PH3)3]2+ + PH3	34.8
[Zn(H2O)]2+	→	Zn2+ + H2O	104.3
[Zn(H2O)4]2+	→	[Zn(H2O)3]2+ + H2O	43.0
[Zn(CO)]2+	→	Zn2+ + CO	81.6
[Zn(CO)2]2+	→	[Zn(CO)]2+ + 2CO	149.8
[Zn(CO)3]2+	→	[Zn(CO)]2+ + 3CO	190.1
[Zn(SH2)]2+	→	Zn2+ + SH2	126.1
[Zn(SH2)2]2+	→	Zn2+ + 2SH2	210.5
[Zn(PH3)4]2+	→	[Zn(PH3)3]2+ + PH3	30.5

See text for the methods used as references for different properties.

Condensed-Phase Benchmark Calculations

To supplement the gas-phase calculations, simulations are carried out with class="Chemical">DFTB3/MM to investigate if the class="Chemical">parametrized model works in an exclass="Chemical">plicit class="Chemical">pan class="Chemical">condensed-phase environment. They are important tests because the ultimate goal is to use DFTB3/3OB in enzyme simulations.

Mg2+ and Zn2+ Ion Solvation and pKa Calculations

Both MM and QM/MM simulations are performed using CHARMM.[100] The generalized solvent boundary potential (class="Chemical">GSBP)[101,102] is used to treat long-range electrostatic interactions in MD simulations. The system is set uclass="Chemical">p with an inner region of 20 Å. All molecules are located in the inner region, and the outer region is described with a class="Chemical">pan class="Chemical">constant dielectric continuum with εw = 80. To be consistent with the GSBP protocol, electrostatic interactions among inner region atoms are treated with the extended electrostatic model[103] in which interactions beyond 12 Å are modeled with multipolar expansions, including the dipolar and quadrupolar terms. The reaction field matrix M is evaluated using 400 spherical harmonics for GSBP simulations. The ion is fixed at the center of the sphere. Water molecules are described with the modified version of TIP3P.[104] In the classical MM simulations, standard CHARMM force field[105] parameters are used for Mg2+ and Zn2+. The QM/MM calculations employ the same water sphere as in MM simulations. The QM region includes the metal ion together with the nearest 6 water molecules (referred to as small QM region) or the nearest 25 water molecules (referred to as large QM region) which include the first and the second solvation shells. All the remaining water molecules are represented by the TIP3P water model. In the large QM region simulations, to ensure the QM water molecules are always closest to the metal ion, Flexible Inner Region Ensemble Separator (FIRES[106]) restraining potential is imposed to any outer sphere MM water molecules that are closer to the ion than the most distant inner sphere QM water molecule. The FIRES potential takes the form EFIRES = (1/2)∑kFIRES(r – Rinner)2 where Rinner is defined as the distance between the metal ion and the most distant inner QM water molecule from it. kFIRES is set to be 100 kcal/(mol·Å2). FIRES restraints are not necessary for the small QM region simulations because water molecules in the first solvation shell are strongly bound to the metal ion. All bonds involving hydrogen are constrained with the SHAKE algorithm,[107] allowing a 1 fs time step for MD propagation. Simulations are equilibrated for 500 ps followed by a production run of 500 ps.

To further evaluate the performance of class="Chemical">DFTB3/3OB for class="Chemical">pan class="Chemical">Mg and Zn from the energetics perspective, we perform calculations of relative pKa of Mg2+/Zn2+-bound water in a water sphere using the thermodynamic integration (TI)[108,109] approach within the dual topology single coordinate (DTSC) scheme; absolute pKa results are discussed in the Supporting Information.[66,110] In this approach, the dominant contribution to the total free energy of deprotonation is from the electrostatic free energy change (ΔFM) associated with converting the acidic proton to a dummy atom (D), i.e., transforming M2+·(H2O)6 to M2+·(H2O)5(HOD)−, where M represents the metal ion, Mg2+ or Zn2+. The corresponding free energy derivative is given bywhich represents the QM/MM energy difference averaged for a specific coupling parameter λ using the same set of coordinates for both protonation states. The total electrostatic energy contribution (ΔFM) is determined by integrating the free energy derivatives (∂F/∂λ) over λ from 0 to 1. Estimation of the absolute pKa is in principle possible with the DTSC-TI protocol, but the solvation free energy of a proton has to be taken into account. To avoid the difficulty of getting accurate experimental/calculated values for this term, we compute the relative pKa between Mg2+·(H2O)6 and Zn2+·(H2O)6 in water so that other contributions are canceled out to a great extent, such as the zero-point energy difference between the protonated and deprotonated states and van der Waals interactions involving the acidic proton.[66,110] Since ΔFM is electrostatic in nature, we expect that the free energy derivative depends largely linearly on the coupling parameter λ in this simple system; this is supported by actual calculations (see below). We carry out pKa calculations with both small QM region and large QM region as defined above. To avoid dissociation of hydroxide from the metal ion in the fully deprotonated state (λ = 1 window M2+·(H2O)5(HOD)−), a NOE potential is added to the distance between M2+ and OH–. The NOE potential takes the formin which rmin and rmax set the interval between which the restraining potential is zero; they are taken to be 1.8 and 2.5 Å, respectively. kmax is set to be 1000 kcal/(mol·Å2). To prevent the dummy atom D from being too close to the metal site in the fully deprotonated state, a weak angular constraint is applied to the M–O–D angle to keep it no smaller than 120°.

Myosin-II Enzyme Setup: Test for a Mg Site

The enzyme model is built from a high-resolution X-ray structure for the class="Chemical">Mg·class="Chemical">pan class="Chemical">ADP·vanadata bound myosin[111] (PDB code 1VOM), which corresponds to the hydrolyzing state of the D. discoideum myosin-II motor domain.[112,113] As described in our earlier work,[68,114,115] starting from the PDB structure, the ADP·vanadata is replaced by ATP or ADP·Pi. The QM region includes the triphosphate and part of the ribose group of ATP, three water molecules (include the lytic water) in the active site, Ser 181, Ser 236, the Mg2+ ion, and all its ligands (Thr 186 and Ser 237 and two water molecules). Only side chains of protein residues are included in the QM region, and link atoms are introduced to saturate the valence of the QM boundary atoms. ATP is assumed to be fully deprotonated, and the titratable groups are kept in their normal protonation states (i.e., Lys, Arg are protonated, Asp, Glu are deprotonated) which are consistent with a simple estimate of pKa using the Poisson–Boltzmann (PB) approach.[116] GSBP is used as the simulation protocol. The system is divided into a 24 Å spherical inner region centered at the Mg ion. Newtonian equations-of-motion are solved for the MD region (within 20 Å), and Langevin equations-of-motion are solved for the buffer region (20–24 Å) with a temperature bath of 300 K. Protein atoms in the buffer region are harmonically constrained with force constants determined from the crystallographic B-factors.[117] Protein atoms in the MM region are described by the all-atom CHARMM22 force field,[105] and water molecules are described with the TIP3P model. All bonds involving hydrogen are constrained using SHAKE, and the time step is set to 1 fs. All the water molecules in the inner region are subject to a weak GEO type of restraining potential to keep them inside the inner sphere with the MMFP module in CHARMM. The static field due to outer-region atoms, ϕs(io), is evaluated with the linear PB equations. The reaction field matrix M is evaluated with 625 spherical harmonics. In the PB calculations, the protein dielectric constant is set as 1; the water dielectric constant is set as 80; and 0.0 M salt concentration is used. Our previous[118] and recent[119] analyses indicated that the GSBP protocol is reliable for a site well shielded from the bulk solvent.

To study the structural flexibility in the active site of myosin before and after class="Chemical">ATP hydroclass="Chemical">pan class="Chemical">lysis, potential of mean force (PMF) simulations have been carried out. The reaction coordinate is defined as the distance between Mg2+ and atom O in Ser 237 (see Figure 5). The umbrella sampling approach[120] is used to restrain the system along the reaction coordinate with a force constant of 150 kcal/mol·Å–2. Thirty-two windows are used for PMF, and a 500 ps simulation is performed for each window. Convergence of PMF is monitored by examining the overlap of reaction coordinate distributions sampled in different windows. Weighted histogram analysis (WHAM)[121] is used to postanalyze the probability distribution to obtain the PMF along the reaction coordinate. Errors are calculated from block average.

Figure 5

Structural properties of the Myosin active site during equilibrium QM/MM MD simulations with different nucleotide chemical states (ATP or ADP·Pi). (a) Overlay of average DFTB3/MM structure and crystal structure (gray) of Myosin-II motor domain active site in the ATP state. The VO4 moiety in the crystal structure (PDB code 1VOM) is replaced by a PO3 phosphate group and a presumptive lytic water molecule. Instantaneous distances between Mg2+ and oxygen atoms in its four nonwater ligands at (b) ATP state and (c) ADP·Pi state. Ser237 refers to OSer237 and Thr186 refers to OThr186. (d) PMF calculations for the ATP/ADP·Pi states. The reaction coordinate is the distance between Mg2+ and OSer237.

Alkaline Phosphatase Enzyme Setup: Test for a Bimetallic Zinc Site

The enzyme model is class="Chemical">constructed based on the X-ray structures for the class="Chemical">pan class="Species">E. coli Alkaline Phosphatase (AP) mutant R166S with bound inorganic phosphate at 2.05 Å resolution (PDB code 3CMR(122)). The substrate methyl p-nitrophenyl phosphate (MpNPP–) is first “mutated” to the orientation with the −OMe group oriented toward the magnesium ion (denoted as α orientation) starting from the PDB structure. All basic and acidic amino acids are kept in their physiological protonation states except for Ser102 in AP, which is accepted to be the nucleophile and deprotonated in the reactive complex. Water molecules are added following the standard protocol of superimposing the system with a water droplet of 25 Å radius centered at Zn12+ (see Figure 6 for atomic labels) and removing water molecules within 2.8 Å from any heavy atoms resolved in the crystal structure. GSBP is used to treat long-range electrostatic interactions in MD simulations. In the QM/MM simulations, as described in details in our previous work,[48,49] the QM region includes the two zinc ions and their six ligands (Asp51, Asp369, His370, Asp327, His412, His331), Ser102, and the substrate MpNPP–. Only side chains of protein residues are included in the QM region, and link atoms are added between Cα and Cβ atoms. A NOE potential is added to the C–O bonds in Asp51 in AP to prevent artificial polarization.[48] The integration time step is 1 fs, and all bonds involving hydrogens are constrained using SHAKE. The DFTB3/MM results are also compared to MM simulations using a conventional nonbonded zinc model[14] (referred to as a Coulomb scheme) or short–long effective functions (SLEF1)[25] model developed by Zhang and co-workers. Protein atoms in the MM region are described by the all-atom CHARMM22 force field, and water molecules are described with the TIP3P model.

Figure 6

Structural properties of the AP active site during equilibrium QM/MM MD Simulations. A snapshot for the reactant state with (a) Coulomb MM force field, (b) SLEF1 zinc MM force field, (c) DFTB3/3OB, and (d) overlay of crystal and DFTB3/MM structure (blue), with average key distances in angstroms labeled. Asp369, His370, and His412 are omitted for clarity. Substrate inorganic phosphate in the crystal structure (3CMR) is not shown.

To further benchmark the applicability of class="Chemical">DFTB3/3OB to the reaction of interest, we also study an active site model that includes all QM atoms in the QM/MM enzyme model. The Cβ class="Chemical">pan class="Chemical">carbons are fixed at their positions in the crystal structure during geometry optimization; the positions of the link atoms used to saturate the valence of the Cβ atoms in the active site model are also fixed. The reactant (Michaelis) complex and transition state are located for MpNPP– (methyl p-nitrophenyl phosphate), MmNPP– (methyl m-nitrophenyl phosphate), and MPP– (methyl phenyl phosphate) using DFTB3 and B3LYP with the 6-31+G(d,p) basis set. The minimum energy path (MEP) calculations are carried out by one-dimensional adiabatic mapping at the DFTB3 level; the reaction coordinate is the antisymmetric stretch involving the breaking and forming P–O bonds. Subsequently, the saddle point is further refined by conjugated peak refinement (CPR[123]) to obtain precise transition state structure. Single-point energy calculations at DFTB3 and B3LYP geometries are then carried out using B3LYP, M06,[124] PBE, and MP2 methods using a larger basis set of 6-311++G(d,p). The D3[125] dispersion corrections are added for B3LYP, M06, and PBE methods.

Results and Discussions

In the following subsections, we present benchmark calculations for class="Chemical">DFTB3/3OB regarding geometries, class="Chemical">pan class="Chemical">metal–ligand dissociation energies, and proton affinities for Mg/Zn containing molecules and compare the results with commonly used ab initio methods in the gas phase. We also report DFTB3 results when the old set of MIO parameters is used:[75] DFTB3/MIO use parameters as defined in ref (75), Ud = −0.23, −0.16, −0.13, −0.19, −0.14, and −0.09 au for C, H, N, O, P, and S and ζ = 4.2; Ud for Mg/Zn is fitted to −0.1 to achieve the best sequential ligand dissociation energies and proton affinities. Unless noted otherwise, energies are given for geometries that are optimized at the respective level of theory. For PM6, we note that the choice of heat of formation for protons is essential to proton affinities; we take −54.0 kcal/mol for protons as suggested by the MOPAC program. To compare the potential energy surfaces between different computation levels, zero-point energies (ZPEs) are not included. For the calculations of noncovalent interactions, dispersion interactions are not included because most test molecules discussed here are small, although we are aware of the importance of dispersion in general. We explore explicitly the effects of dispersion to a nontrivial active site model of AP that involves zinc. For DFTB3 calculations, we have used our in-house DFTB code, and all PM6, DFT, and MP2 calculations are performed using the Gaussian09 software package.[126] Following the gas-phase benchmarks, we show results on Mg2+/Zn2+-water structural properties and relative pKa in solution. Finally we apply DFTB3/3OB to enzyme systems to further test the new parameters.

Magnesium

Geometries

The “small” class="Chemical">magnesium test set is class="Chemical">pan class="Chemical">composed of 55 molecules, which are simplified biologically relevant ligated magnesium species of the form Mg2+[X] with X = H2O, NH3, SH2, PH3, CO for n = 1–4/6 and X = OH–, SH–, NH2–, and PH2– for n = 1–2. Table 3 summarizes the results, and further details are included in the Supporting Information. Good agreement is seen between MP2 and B3LYP with a small or large basis set. Of the semiempirical methods, geometries are described very well by DFTB3 models and particularly well by DFTB3/3OB. The MADs in metal–ligand distances are 0.013, 0.017, and 0.050 Å for DFTB3/3OB, DFTB/MIO, and PM6, respectively. 3OB improves over MIO, although the magnitude of improvement is modest. The largest metal–ligand distance error of 0.072 Å from DFTB3/3OB is found in Mg–S distance in [Mg(SH2)]2+, which is likely due to the minimal basis nature of DFTB3 for a polarizable ligand. Other deviations in metal–ligand distances are substantially smaller in DFTB3/3OB. Molecules containing S or P are described poorly by PM6 within our test case;[78] PM6 systematically underestimates Mg–S distance by 0.2 Å and Mg–P distance by 0.1 Å. For angles, DFTB3/3OB differs from the reference data by a MAD of 2.9°, which is of the same order with other semiempirical methods. The largest angle derivation of 45.3° is for angle Mg–O–H in [Mg(OH)]+, for which DFTB3/3OB predicts a linear rather than bent structure; interestingly, MP2 also predicts a structure almost linear, while PM6 is close to the B3LYP result (see Supporting Information). Future development of DFTB3 to include multipoles may improve the geometry in this case. It is worth pointing out that the Mg–S–H angle is generally overestimated by ∼13° in DFTB3/3OB. More statistics including MAD for different metal–ligand types are given in the Supporting Information.

Table 3

Mean and Maximum Absolute Deviation of Metal–Ligand Lengths and Angles in Comparison to B3LYP/aug-cc-pVTZ for the Magnesium Test Set

propertya	Nb	MP2c	B3LYPd	PM6	DFTB3/MIO	DFTB3/3OB
r (Å)	120	0.009	0.009	0.050	0.017	0.013
rmax (Å)		0.051	0.024	0.284	0.141	0.072
a (deg)	117	0.9	0.5	3.4	3.6	2.9
amax (deg)		23.7	10.2	38.7	41.3	45.3

Max stands for maximum absolute deviation.

Number of comparisons.

Basis set aug-cc-pVTZ.

Basis set 6-31+G(d,p).

Number of pan class="Chemical">comclass="Chemical">parisons.

Sequential Metal–Ligand Dissociation Energies

We define sequential ligand dissociation energies (pan class="Chemical">sBDEs) class="Chemical">pan class="Chemical">corresponding to the reactionwhere M is Mg2+/Zn2+ and L is a neutral or charged ligand.

Table 4 shows the sequential ligand dissociation energies for class="Chemical">comclass="Chemical">pounds class="Chemical">pan class="Chemical">containing Mg2+. Values obtained for DFTB3/3OB and other methods are compared to those obtained using G3B3. The differences between B3LYP, MP2, and G3B3 are relatively small; however, complexes with charged ligands give several deviations greater than 5 kcal/mol. DFTB3/3OB shows a MAD of 2.4 kcal/mol for complexes with neutral ligands and 13.8 kcal/mol for complexes with charged ligands. Compared to DFTB3/MIO, DFTB3/3OB improves substantially for complexes involving P, S elements and charged ligands. For the other semiempirical method PM6, the deviation is 9.9 kcal/mol compared to G3B3. The large error is anticipated because geometries containing S and P are often poor. Note that we compare potential energy differences for all the methods except PM6 which intrinsically considers heats of formation. We have also carried out single-point G3B3 calculations at DFTB3/3OB geometries, and the MAD for Mg2+ compounds with neutral and charged ligands is less than 1 kcal/mol. These results highlight the similarities between DFTB3 and B3LYP optimized structures, suggesting that DFTB3/3OB is of particular value for the study of biologically relevant magnesium-containing molecules.

Table 4

Sequential Ligand Dissociation Energies of Magnesium -Containing Molecules Compared to G3B3a

moleculeb	G3B3	MP2c	B3LYPc	PM6	DFTB3/MIO	DFTB3/3OB	G3B3//DFTB3d
[Mg(H2O)2]2+	73.3	–3.5	–0.1	–10.9	–3.4	+0.5	+0.6
[Mg(H2O)3]2+	59.2	–2.0	–1.1	–11.7	–1.5	+1.4	–0.3
[Mg(H2O)4]2+	49.1	–2.3	–3.0	–11.0	–3.3	–1.9	+1.2
[Mg(H2O)5]2+	34.7	–2.2	–5.0	–7.4	–6.7	–6.4	–0.2
[Mg(H2O)6]2+	32.9	–2.5	–5.9	–11.0	–6.8	–6.4	+0.4
[Mg(CO)2]2+	46.3	–0.9	–0.3	+10.0	–2.9	–3.8	–0.6
[Mg(CO)3]2+	36.1	+0.5	–0.9	+1.3	–0.3	–2.3	–0.7
[Mg(CO)4]2+	30.9	–0.1	–2.8	–1.5	–0.1	–2.6	–0.9
[Mg(NH3)2]2+	83.9	–3.4	–0.3	–4.4	–13.8	–6.4	+0.5
[Mg(NH3)3]2+	63.1	–1.8	–2.0	–9.2	–7.5	–1.4	+0.2
[Mg(NH3)4]2+	50.1	–1.8	–3.7	–8.8	–3.9	+1.2	+0.2
[Mg(NH3)5]2+	28.6	–1.0	–5.3	–2.9	–5.8	–3.5	–0.2
[Mg(NH3)6]2+	26.3	–1.5	–6.2	–5.5	–5.6	–2.5	–0.2
[Mg(SH2)2]2+	62.4	–3.0	–0.3	–7.0	–7.9	–2.9	–1.2
[Mg(SH2)3]2+	43.9	–1.1	–2.6	–12.1	–4.5	+0.1	–0.6
[Mg(SH2)4]2+	34.5	–0.8	–4.4	–13.3	–2.8	+1.2	+1.7
[Mg(SH2)5]2+	20.0	+0.2	–6.1	–12.6	+0.0	+1.1	–0.8
[Mg(SH2)6]2+	20.7	–0.5	–7.6	–14.0	–1.8	–0.1	+0.5
[Mg(PH3)2]2+	66.9	–2.3	–0.5	–2.4	–8.6	–0.1	+0.7
[Mg(PH3)3]2+	44.5	–0.8	–2.9	–18.6	–5.4	+2.3	+0.7
[Mg(PH3)4]2+	34.8	–0.6	–4.5	–17.9	–3.8	+3.2	+0.6
[Mg(OH)2]	243.1	–7.2	–4.6	–7.1	–21.0	–14.0	+0.1
[Mg(NH2)2]	237.9	–7.3	–6.6	–4.6	–9.4	–12.0	+0.9
[Mg(SH)2]	201.2	–3.1	–5.0	–30.3	–16.1	–15.8	+0.1
[Mg(PH2)2]	188.8	–1.8	–5.2	+11.0	–20.9	–13.3	+0.6
MAD		2.1	3.5	9.9	6.6	4.3	0.6
MAX		7.3	7.6	30.3	21.0	15.8	1.7

Energies except PM6 are calculated at 0 K excluding zero-point energy and thermal corrections. All numbers are given in kcal/mol.

The nondissociated is listed.

Basis set aug-cc-pVTZ.

G3B3 single-point calculations on top of DFTB3/3OB geometries.

Proton Affinities

Table 5 lists the proton affinities for 22 biological relevant molecules class="Chemical">Mg2+(XH) where X = O, N, S, and P obtained from G3B3, class="Chemical">pan class="Chemical">DFTB3/3OB, and other methods. MADs for MP2 and B3LYP are generally small (∼1 kcal/mol). The discrepancy between DFTB3/3OB and G3B3 is 5 kcal/mol, whereas DFTB3/MIO and PM6 yield slightly larger MADs of 7.8 and 6.6 kcal/mol, respectively. We note that G3B3 single points at DFTB3/3OB optimized structures lead to small MAD of 1.5 kcal/mol and MAX of 3.9 kcal/mol (due to the unsatisfactory deprotonated [Mg(SH2)]2+ structures) compared to the original G3B3, again highlighting the good quality of DFTB3/3OB structures. Future improvements should focus on the nitrogen-containing species, an issue we noted in previous studies.[75,77]

Table 5

Gas-Phase Proton Affinities of Magnesium-Containing Molecules Compared to G3B3a

moleculeb	G3B3	MP2c	B3LYPc	PM6	DFTB3/MIO	DFTB/3OB	G3B3//DFTB3d
[Mg(H2O)1]2+	95.6	+0.6	–3.5	+0.1	–5.5	–13.9	–0.8
[Mg(H2O)2]2+	115.6	+0.0	–1.4	–5.4	–2.2	–9.0	–0.9
[Mg(H2O)3]2+	133.1	+0.6	+0.2	–7.1	+3.0	–1.2	–0.4
[Mg(H2O)4]2+	149.2	–0.6	–0.1	–11.4	+4.8	+2.3	+0.8
[Mg(H2O)5]2+	160.3	–2.2	–1.9	–15.7	+3.2	+2.0	+1.1
[Mg(H2O)6]2+	170.5	–3.0	–2.2	–24.3	+0.2	–0.2	+0.8
[Mg(NH3)1]2+	120.2	+1.5	–3.5	+2.6	–10.5	–15.7	+0.0
[Mg(NH3)2]2+	144.1	+0.8	–0.4	+1.5	–10.2	–14.1	+0.2
[Mg(NH3)3]2+	163.9	+0.5	+0.3	+1.9	–10.6	–13.6	+0.5
[Mg(NH3)4]2+	181.3	–0.1	+0.6	+1.0	–6.4	–7.5	+1.0
[Mg(NH3)5]2+	190.5	–0.4	+0.9	+4.9	–5.2	–6.1	+1.4
[Mg(NH3)6]2+	200.7	+1.8	–0.5	+5.0	–6.2	–5.4	+2.3
[Mg(SH2)1]2+	83.4	+2.1	–1.8	+2.4	–6.7	–0.7	–3.6
[Mg(SH2)2]2+	107.0	+0.8	+0.5	+4.6	–6.2	+0.8	–2.9
[Mg(SH2)3]2+	123.1	+0.4	+1.5	+3.3	–6.2	+2.1	–3.5
[Mg(SH2)4]2+	135.4	+0.0	+2.0	+1.0	–5.7	+4.2	–1.7
[Mg(SH2)5]2+	143.1	–0.2	+2.0	–5.7	–6.9	+3.9	–3.9
[Mg(SH2)6]2+	151.1	–0.6	+1.8	–12.4	–8.7	+3.5	–3.0
[Mg(PH3)1]2+	96.8	+4.1	–1.8	+12.6	–12.8	–1.1	+0.6
[Mg(PH3)2]2+	129.4	+1.7	+0.4	+10.2	–15.3	–0.4	+1.4
[Mg(PH3)3]2+	150.1	+1.0	+1.2	+8.2	–17.6	–1.7	+1.4
[Mg(PH3)4]2+	164.2	+0.7	+1.7	+3.0	–18.1	–0.7	+1.2
MAD		1.1	1.4	6.6	7.8	5.0	1.5
MAX		4.1	3.5	24.3	18.1	15.7	3.9

Energies except PM6 are calculated at 0 K excluding zero-point energy and thermal corrections. All numbers are given in kcal/mol.

The protonated form is listed.

Basis set aug-cc-pVTZ.

G3B3 single-point calculations on top of DFTB3/3OB geometries.

Energies except pan class="Chemical">PM6 are calculated at 0 K excluding zero-class="Chemical">point energy and thermal class="Chemical">pan class="Chemical">corrections. All numbers are given in kcal/mol.

G3B3 single-point calculations on top of pan class="Chemical">DFTB3/3OB geometries.

Energies except pan class="Chemical">PM6 are calculated at 0 K excluding zero-class="Chemical">point energy and thermal class="Chemical">pan class="Chemical">corrections. All numbers are given in kcal/mol.

G3B3 single-point calculations on top of pan class="Chemical">DFTB3/3OB geometries.

Large Test Set

To further test class="Chemical">DFTB3/3OB class="Chemical">parameters we have class="Chemical">pan class="Chemical">compiled a “large test set” that models biologically relevant ligands of magnesium with larger molecules than those discussed above. The set mainly consists of model magnesium compounds that are simplified active sites of crystal structures collected in a survey.[127] Although Mg2+ is generally 6-coordinated, complexes with coordination number of 4 and 5 are also investigated. All complexes are constructed in the way that all ligands are coordinated to the central magnesium ion rather than via hydrogen bonds. CH3COO–/HCOO– is used to model the Asp/Glu side chain, CH3OH to represent the Ser/Thr side chain, and OCHNH2 to represent carbonyls. In total, the set contains 58 magnesium compounds (see Supporting Information).

pan class="Chemical">Metal–ligand distances calculated with class="Chemical">pan class="Chemical">DFTB3/3OB deviate only by 0.018 Å on average and 0.125 Å at most (Table 6). The case with the largest error is the Mg–S distance in [Mg(SH2)4(SH)]−. Besides this outlier, all other metal–ligand length errors are substantially smaller within 0.1 Å. DFTB3/3OB outperforms DFTB3/MIO in the large test set. B3LYP with the small basis set 6-31+G(d,p) deviates on average by 0.011 Å, which is comparable with DFTB3/3OB, and the maximum deviation is 0.033 Å, which is much better than 3OB. PM6 strongly underestimates Mg–O distance by 0.04 Å. More statistics including mean absolute deviations for different metal–ligand types are given in the Supporting Information. For the angle, the MADs are 3.1, 4.5, 0.7, and 11.2°, respectively, for DFTB3/3OB, DFTB3/MIO, B3LYP/6-31+G(d,p), and PM6. In the O–Mg–O angles calculated with DFTB3/3OB, about 6% out of 505 angles deviate more than 10°; the large errors stem mainly from erroneous orientation of hydroxide.

Table 6

Mean and Maximum Absolute Deviation of Metal–Ligand Lengths and Angles in Comparison to B3LYP/aug-cc-pVTZ for the Magnesium Large Test Set

propertya	Nb	B3LYPc	PM6	DFTB3/MIO	DFTB3/3OB
r (Å)	329	0.011	0.063	0.041	0.018
rmax (Å)		0.033	1.414	0.481	0.125
a (deg)	658	0.7	11.2	4.5	3.1
amax (deg)		17.4	64.0	37.5	30.8

Max stands for maximum absolute deviation.

Number of comparisons.

Basis set 6-31+G(d,p).

Due to the relatively high affinity of class="Chemical">ATP/class="Chemical">pan class="Chemical">ADP to Mg2+ and the importance of this metal chelate in the ATP-binding enzyme, we have constructed model molecules involving Mg and phosphate/model AT(D)P as an additional test (Figure 1). The precise binding mode of Mg2+ to ATP remains unclear despite extensive spectroscopic studies. NMR studies proposed several models for the Mg·ATP complex: β-monodentate,[128] β,γ-bidentate,[129] α,β,γ-tridentate,[130,131] and a mixture of α,β-, β,γ-, and α,γ-bidentate.[132] Raman and infrared spectra suggested a mixture of β,γ-bidentate and α,β,γ-tridentate.[133] Here two stable binding modes are tested, one with Mg coordinating β- and γ-phosphates and the other with Mg coordinating α-, β-, and γ-phosphates.[134] Geometries are described very well by DFTB3/3OB. The Mg–O distance agrees on average within 0.04 Å with B3LYP/aug-cc-pVTZ calculations. The O–Mg–O, O–P–O, and P–O–P angles differ on average 5, 1, and 7°, respectively.

Figure 1

Optimized structures of Mg·AT(D)P·water model molecules at different QM levels. Distances are given in angstroms. The numbers without parentheses were obtained at the B3LYP aug-cc-pVTZ level; those in the parentheses are DFTB3 values. The B3LYP optimized structures are depicted in colors, whereas DFTB3 optimized structures are colored in blue.

Number of pan class="Chemical">comclass="Chemical">parisons.

Optimized structures of pan class="Chemical">Mg·AT(D)P·class="Chemical">pan class="Chemical">water model molecules at different QM levels. Distances are given in angstroms. The numbers without parentheses were obtained at the B3LYP aug-cc-pVTZ level; those in the parentheses are DFTB3 values. The B3LYP optimized structures are depicted in colors, whereas DFTB3 optimized structures are colored in blue.

Ligand dissociation energies are class="Chemical">comclass="Chemical">pared in Table 7 (for detailed data, see Suclass="Chemical">pclass="Chemical">porting Information). One neutral class="Chemical">pan class="Chemical">oxygen ligand is removed from the magnesium complexes selected from the large test set. Because the reactions are chosen somewhat arbitrarily, the statistics given in Table 7 have to be considered with care. MADs are reasonably small for all semiempirical methods; the MAD for DFTB3/3OB is 4.6 kcal/mol, even smaller than B3LYP/aug-cc-pVTZ with 5.3 kcal/mol. Interestingly, DFTB3/3OB and B3LYP turn out to systematically underestimate the dissociation energy of the oxygen ligand; the errors seem to be inherited from the sBDE reactions where one water molecule is taken away from [Mg(H2O)5]2+ or [Mg(H2O)6]2+. Further fitting Ud or εp cannot improve these reaction energies significantly. To demonstrate the good performance of DFTB3 for structural properties, we have carried out single-point G3B3 energies at DFTB3 geometries, and the MAD is only 0.8 kca/mol.

Table 7

Error Statistics for the Ligand Dissociation Energies and Proton Affinities of Magnesium Containing Large Molecules Compared to G3B3a

	B3LYPb	PM6	DFTB3/MIO	DFTB3/3OB	G3B3//DFTB3c
Ligand Dissociation	Energies
MAD	5.3	5.1	3.5	4.6	0.8
MAX	7.9	13.0	6.6	7.7	2.2
Proton Affinities
MAD	1.9	17.2	2.7	2.6	0.9
MAX	4.8	25.0	8.0	6.1	2.3

For detailed data, see Supporting Information.

Basis set aug-cc-pVTZ.

G3B3 single-point calculations on top of DFTB3/3OB geometries.

Proton affinities are also class="Chemical">comclass="Chemical">piled in Table 7. The MAD for class="Chemical">pan class="Chemical">DFTB3/3OB is 2.6 kcal/mol, similar to B3LYP/aug-cc-pVTZ with 1.9 kcal/mol. Despite the structural difference for the deprotonated structures mentioned above, the G3B3 single-point calculations at DFTB3/3OB geometries show excellent statistics for proton affinities, demonstrating the strength of DFTB3/3OB in reproducing B3LYP/aug-cc-pVTZ geometries or at least reproducing relative energies for multiple points on the potential energy surface. Similar to findings from our previous studies,[78] PM6 appears to be less reliable for proton affinities with large errors.

G3B3 single-point calculations on top of pan class="Chemical">DFTB3/3OB geometries.

Zinc

Typical biological ligands for class="Chemical">Zn2+ include class="Chemical">pan class="Chemical">water, cysteine, histidine, and glutamic/aspartic acids. Therefore, we performed calculations for a series of complexes that involve the zinc ion and small molecule models for these ligands. Similar to what has been done for magnesium, we have compiled a test set in the form of Zn2+[X] with X = H2O, NH3, SH2, PH3, CO, NH=CH2, and SH–CH3 for n = 1–4/6 and X = OH–, SH–, NH2–, PH2–, NCH2–, and SCH3– for n = 1–2. Here two extra ligands NH=CH2 and SH–CH3 (and their deprotonated form) are added to model histidine and cysteine side chains. The structural comparisons are given in Table 8.

Table 8

Mean and Maximum Absolute Deviation of Metal–Ligand Lengths and Angles in Comparison to B3LYP/aug-cc-pVTZ for the Zinc Test Set

propertya	Nb	MP2c	B3LYPd	PM6	DFTB3/MIOe	DFTB3/3OB
r (Å)	163	0.017	0.004	0.065	0.027	0.012
rmax (Å)		0.077	0.011	0.553	0.139	0.091
a (deg)	180	0.6	0.2	4.2	3.0	2.3
amax (deg)		4.3	4.2	64.0	16.2	14.6

Max stands for maximum absolute deviation.

Number of comparisons.

Basis set aug-cc-pVTZ.

Basis set 6-31+G(d,p).

Molecules containing PH2– or H– do not converge and are excluded from statistics.

The class="Chemical">DFTB3/3OB results are tested against B3LYP/aug-cc-class="Chemical">pVTZ for geometries and energetics. Overall good geometries are obtained from class="Chemical">pan class="Chemical">DFTB3/3OB with metal–ligand distance MAD of 0.012 Å and angle MAD of 2.3°, which are notable improvements over DFTB3/MIO. The largest deviation 0.091 Å of metal–ligand length in DFTB3/3OB is for [Zn(N=CH2)]+. The distance between Zn and N is underestimated by ∼0.1 Å if r(ZnN) is shorter than 1.95 Å due to the imperfect Zn–N repulsive potential at the short-range. Similar to cases for magnesium, angles Zn–S–H are uniformly shifted up by ∼10° in complexes containing ligand SH2 or SH–CH3 in DFTB3/3OB optimized structures. B3LYP/6-31+G(d,p) is structurally very similar to B3LYP/aug-cc-pVTZ with a distance MAD of merely 0.004 Å. Nevertheless, MP2/aug-cc-pVTZ is rather different from B3LYP/aug-cc-pVTZ with a MAD of 0.017 Å; metal–ligand distances between zinc and C, N, O, P, S are uniformly underestimated by 0.02–0.03 Å. Molecular geometries are described rather poorly by PM6 in our test set; for example, the metal–ligand lengths between zinc and phosphorus deviate on average 0.3 Å. More statistics for different metal–ligand types are given in the Supporting Information.

Number of pan class="Chemical">comclass="Chemical">parisons.

Molecules pan class="Chemical">containing PH2– or H– do not class="Chemical">pan class="Chemical">converge and are excluded from statistics.

Sequential Ligand Dissociation Energies

Sequential ligand dissociation energies for zinc species with neutral and negatively chclass="Chemical">arged ligands are summarized in Table 9. class="Chemical">pan class="Chemical">DFTB3/3OB shows a MAD of 3.6 kcal/mol compared to B3LYP/aug-cc-pVTZ. Again, deviations close to or greater than 10 kcal/mol are found for charged ligands. The difference between the small and large basis set in B3LYP is rather small, and even charged ligands produce errors no larger than 2.7 kcal/mol, indicating B3LYP with small basis set may be useful in ab initio QM/MM calculations. There is a slightly large discrepancy between MP2 and B3LYP; it may be attributed to the differences in structures between these two methods as discussed above. PM6 shows large deviations for sBDE with a MAD of 16.8 kcal/mol. With the unphysical Zn–P geometries removed from statistics, the MAD for PM6 drops to 9.1 kcal/mol, still substantially larger than that for DFTB3. We have also performed single-point calculations of B3LYP/aug-cc-pVTZ on top of DFTB3 geometries. The results are remarkably good with a MAD of only 0.5 kcal/mol.

Table 9

Sequential Ligand Dissociation Energies of Zinc Containing Molecules Compared to B3LYP/aug-cc-pVTZa

moleculeb	B3LYPc	MP2c	B3LYPd	PM6	DFTB3/MIO	DFTB3/3OB	B3LYP//DFTB3e
[Zn(H2O)2]2+	89.7	–2.9	+2.0	–15.7	–5.2	+4.3	+0.2
[Zn(H2O)3]2+	56.4	+1.2	+1.5	+1.4	–1.1	–2.2	–0.1
[Zn(H2O)4]2+	43.0	+2.2	+1.9	+3.2	–2.0	–0.9	+0.8
[Zn(H2O)5]2+	25.3	+3.2	+2.6	+1.8	–3.8	–1.7	–0.4
[Zn(H2O)6]2+	23.4	+3.7	+2.3	–3.2	–3.9	–1.5	+0.3
[Zn(CO)2]2+	68.2	+0.4	–0.4	+28.0	–5.8	–3.7	+0.4
[Zn(CO)3]2+	40.3	+3.6	+0.1	+15.3	+1.2	–7.3	+0.3
[Zn(CO)4]2+	30.6	+4.5	+0.3	+16.9	+3.9	–4.6	–0.2
[Zn(NH3)2]2+	108.3	–0.6	+1.9	–8.1	–7.8	–4.3	+0.4
[Zn(NH3)3]2+	59.7	+2.9	+2.4	+8.0	+0.2	–2.6	+0.3
[Zn(NH3)4]2+	43.8	+3.6	+2.0	+9.4	+4.1	+2.8	+0.1
[Zn(NH3)5]2+	15.3	+3.5	+2.1	+1.7	–2.7	+3.2	–0.2
[Zn(NH3)6]2+	13.1	+4.3	+1.8	+1.4	–2.3	+4.1	–0.4
[Zn(SH2)2]2+	84.4	+0.0	+0.5	–5.8	–2.0	+1.0	–1.4
[Zn(SH2)3]2+	42.9	+3.9	+0.5	–1.6	+3.5	+1.8	–0.1
[Zn(SH2)4]2+	29.1	+5.2	+1.2	+0.4	+4.8	+4.9	+2.0
[Zn(PH3)2]2+	89.4	+1.9	+1.1	+55.7	–7.1	+1.2	+0.7
[Zn(PH3)3]2+	43.0	+4.9	+0.8	+51.9	+0.6	+3.7	+0.6
[Zn(PH3)4]2+	30.5	+7.0	+0.9	+65.9	+3.0	+5.5	+0.5
[Zn(NH=CH2)2]2+	111.9	+0.5	+0.9	–2.2	–1.7	+0.1	+0.1
[Zn(NH=CH2)3]2+	60.1	+4.4	+1.2	+10.5	+1.3	–1.5	–0.4
[Zn(NH=CH2)4]2+	43.5	+5.9	+1.1	+9.8	+3.7	+2.8	+0.4
[Zn(SH–CH3)2]2+	92.2	+2.1	+0.1	–8.2	–5.8	–2.8	–0.3
[Zn(SH–CH3)3]2+	46.0	+6.4	+0.5	–1.2	+0.5	–0.4	+0.4
[Zn(SH–CH3)4]2+	31.1	+8.1	+1.1	–0.1	+2.3	+2.9	+0.3
[Zn(OH)2]	253.8	+1.9	+2.7	+12.0	–15.1	–3.4	+0.5
[Zn(NH2)2]	248.6	+4.9	+2.6	+15.9	+11.5	–4.5	+3.1
[Zn(SH)2]	213.2	+6.1	+0.2	–22.1	–6.7	–9.1	–0.2
[Zn(PH2)2]	199.3	+7.1	+0.8	+108.8	-f	–11.7	+0.9
[Zn(OCH3)2]	235.4	+7.6	+2.2	+6.9	–10.0	–0.5	+0.2
[Zn(NCH2)2]	220.6	+10.9	+0.6	+25.7	+13.8	+8.0	+1.1
[Zn(SCH3)2]	207.4	+8.2	+0.5	–18.7	–7.8	–7.5	+0.0
MAD		4.2	1.3	16.8	4.7	3.6	0.5
MAX		10.9	2.7	108.8	15.1	11.7	3.1

Energies except PM6 are calculated at 0 K excluding zero-point energy and thermal corrections. All numbers are given in kcal/mol.

The nondissociated is listed.

Basis set aug-cc-pVTZ.

Basis set 6-31+G(d,p).

B3LYP/aug-cc-pVTZ single-point calculations on top of DFTB3/3OB geometries.

Molecule does not converge and is excluded from the statistics.

Table 10 shows the protonation energies calculated for 24 biologically relevant ligand model molecules obtained using B3LYP, class="Species">MP2, and semiemclass="Chemical">pirical methods. As exclass="Chemical">pected, B3LYP/6-31+G(d,class="Chemical">p) class="Chemical">performs very well class="Chemical">pan class="Chemical">compared to B3LYP with a large basis set. The MAD for DFTB3 is 4.2 kcal/mol, which is substantially better than DFTB3/MIO and PM6; the errors are small enough to render this approach useful. High-level single-point correction at DFTB3/3OB geometries once again delivers excellent proton affinities except the ones involving the SH– ligand whose geometries are less well described as discussed above.

Table 10

Gas-Phase Proton Affinities of Zinc Containing Molecules Compared to B3LYP/aug-cc-pVTZa

moleculeb	B3LYPc	MP2c	B3LYPd	PM6	DFTB3/MIO	DFTB/3OB	B3LYP//DFTB3e
[Zn(H2O)1]2+	61.7	+5.4	–0.6	+7.8	–1.2	–7.7	+0.8
[Zn(H2O)2]2+	94.2	+0.6	–0.3	–2.5	–1.1	–2.8	+0.5
[Zn(H2O)3]2+	121.3	–0.7	–0.5	–9.6	+1.9	–1.2	+0.1
[Zn(H2O)4]2+	141.5	–1.6	–0.7	–16.9	+2.3	–1.3	+0.7
[Zn(NH3)1]2+	86.9	+4.1	+0.1	+6.5	–13.6	–10.8	+1.8
[Zn(NH3)2]2+	127.3	–0.8	+0.7	–2.9	–18.5	–10.7	+0.6
[Zn(NH3)3]2+	155.0	–1.7	+1.0	–9.7	–18.0	–10.6	+0.2
[Zn(NH3)4]2+	178.1	–2.2	+1.0	–4.8	–10.3	–8.4	+0.7
[Zn(SH2)1]2+	56.6	+1.6	–1.7	+16.1	–13.0	–5.4	–4.0
[Zn(SH2)2]2+	93.2	–2.0	–1.5	+12.5	–10.6	–0.2	–4.0
[Zn(SH2)3]2+	115.0	–2.6	–2.0	+7.7	–9.8	–0.5	–4.0
[Zn(SH2)4]2+	131.1	–3.0	–2.2	+6.0	–9.1	–1.0	–2.1
[Zn(PH3)1]2+	66.9	+3.2	–0.7	+21.8	-f	–1.5	+0.3
[Zn(PH3)2]2+	115.6	–0.2	–0.3	–5.7	-f	–0.9	+0.8
[Zn(PH3)3]2+	140.8	–1.0	–0.1	–30.0	-f	–2.7	+0.7
[Zn(PH3)4]2+	158.0	–1.6	+0.1	+14.3	–18.1	–3.2	+1.0
[Zn(NH=CH2)1]2+	89.1	+6.1	–0.9	+2.3	–16.2	–2.8	+1.1
[Zn(NH=CH2)2]2+	140.9	–0.4	–0.2	–12.6	–14.1	–2.1	+0.4
[Zn(NH=CH2)3]2+	168.6	–1.5	+0.1	–18.1	–14.2	–2.2	+0.0
[Zn(NH=CH2)4]2+	188.1	–2.1	–0.1	–14.3	–14.8	–1.6	+2.4
[Zn(SH–CH3)1]2+	74.2	+0.3	–2.5	+12.1	–4.6	–8.3	–2.5
[Zn(SH–CH3)2]2+	114.9	–3.0	–2.3	+5.6	–10.0	–4.2	–2.1
[Zn(SH–CH3)3]2+	136.7	–4.1	–2.4	+1.1	–10.5	–4.8	–1.7
[Zn(SH–CH3)4]2+	152.6	–6.0	–2.2	+0.3	–2.6	–5.9	–1.4
MAD		2.3	1.0	10.1	10.2	4.2	1.4
MAX		6.1	2.5	30.0	18.5	10.8	4.0

Energies except PM6 are calculated at 0 K excluding zero-point energy and thermal corrections. All numbers are given in kcal/mol.

The protonated form is listed.

Basis set aug-cc-pVTZ.

Basis set 6-31+G(d,p).

B3LYP/aug-cc-pVTZ single-point calculations on top of DFTB3/3OB geometries.

Molecule does not converge and is excluded from the statistics.

Energies except pan class="Chemical">PM6 are calculated at 0 K excluding zero-class="Chemical">point energy and thermal class="Chemical">pan class="Chemical">corrections. All numbers are given in kcal/mol.

B3LYP/aug-cc-pVTZ single-point calculations on top of pan class="Chemical">DFTB3/3OB geometries.

Molecule does not pan class="Chemical">converge and is excluded from the statistics.

Energies except pan class="Chemical">PM6 are calculated at 0 K excluding zero-class="Chemical">point energy and thermal class="Chemical">pan class="Chemical">corrections. All numbers are given in kcal/mol.

B3LYP/aug-cc-pVTZ single-point calculations on top of pan class="Chemical">DFTB3/3OB geometries.

Molecule does not pan class="Chemical">converge and is excluded from the statistics.

Mixed Ligands

One goal to develop class="Chemical">DFTB3/3OB is to describe enzyme active sites. class="Chemical">pan class="Chemical">Zn-containing structures are mainly tetrahedrally bound in metalloenzymes. The most abundant primary ones are CCCC, zinc ion with four coordinating cysteines, followed by CCCH, CCHH, CHHH, HHHH, HHHO, HHOO, HOOO, HHHD, and HHDD, where O stands for water and the remaining are 1-letter amino acid codes.[12] In this subsection, we compare geometries and energetics for 43 zinc complexes with mixed types of simple ligands (including the deprotonated form) to model representative active sites in proteins. Here, H2O is used to represent the Asp/Glu side chain, NH3 to represent the His side chain, and H2S to represent the Cys side chain. More complicated/realistic ligands are discussed in the next subsection on “Large Test Set”. Table 11 summarizes the geometrical properties. The MADs in comparison to B3LYP/aug-cc-pVTZ for 161 metal–ligand lengths are 0.005, 0.019, 0.053, and 0.118 Å for B3LYP/6-31+G(d,p), DFTB3/3OB, DFTB3/MIO, and PM6, respectively. Metal–ligand angles on average deviate by 0.9, 2.6, 3.0, and 6.7°. Again, the 3OB set is a major improvement over the MIO set in metal–ligand lengths. PM6 is not able to give reliable geometries due to a large number of unreasonable metal–ligand distances between zinc and other main group elements.

Table 11

Mean and Maximum Absolute Deviation of Metal-Ligand Lengths and Angles in Comparison to B3LYP/aug-cc-pVTZ for Zinc Mixed Ligand Test Set

propertya	Nb	B3LYPc	PM6	DFTB3/MIOd	DFTB3/3OB
r (Å)	161	0.005	0.118	0.053	0.019
rmax (Å)		0.027	0.514	0.179	0.070
a (deg)	231	0.9	6.7	3.0	2.6
amax (deg)		22.2	57.2	21.6	24.7

max stands for maximum absolute deviation.

number of comparisons.

basis set 6-31+G(d,p).

Molecules containing PH2– do not converge and are excluded from statistics.

number of pan class="Chemical">comclass="Chemical">parisons.

Molecules pan class="Chemical">containing PH2– do not class="Chemical">pan class="Chemical">converge and are excluded from statistics.

class="Chemical">DFTB3/3OB energetics are class="Chemical">predicted with similar accuracy as for class="Chemical">pan class="Chemical">complexes with the same type of ligand as shown in Table 12 (for detailed data, see Supporting Information). Once again, the reactions are collected somewhat arbitrarily so the error analysis needs to be treated with care. B3LYP/6-31+G(d,p) gives small error compared to B3LYP with large basis set, demonstrating basis set superposition error is not severe here. DFTB3 gives rather good results in comparison with the reference (MAD of ∼3 kcal/mol), especially after single-point energies using B3LYP/aug-cc-pVTZ (MAD of 0.9 kcal/mol), demonstrating that the 3OB set is transferable from one type ligand to complicated coordinating environment with multiple types of ligands.

Table 12

Error Statistics for the Ligand Dissociation Energies and Proton Affinities of Zinc Containing Mixed Ligand Molecules Compared to B3LYP/aug-cc-pVTZa

	B3LYPb	PM6	DFTB3/MIO	DFTB3/3OB	B3LYP//DFTB3c
Ligand Dissociation	Energies
MAD	1.9	8.2	3.1	3.2	0.9
MAX	2.8	31.6	6.8	6.5	2.6
Proton Affinities
MAD	0.9	10.4	5.5	3.8	0.9
MAX	1.9	21.7	11.3	7.1	2.4

For detailed data, see Supporting Information.

Basis set 6-31+G(d,p).

B3LYP/aug-cc-PVTZ single-point calculations on top of DFTB3/3OB geometries.

B3LYP/aug-cc-PVTZ single-point calculations on top of pan class="Chemical">DFTB3/3OB geometries.

To further test the 3OB parameter set, we have chosen lclass="Chemical">arger molecules to model the biologically relevant zinc class="Chemical">pan class="Chemical">coordinating environment. In this “large test set”, CH3COO–/HCOO– is used to model the Asp/Glu acid side chain, imidazole/NCH2CH3 to represent His, and SHCH3 to model the Cys side chain. Cys can adopt both protonated and deprotonated form in the active site so that a random mixture of SHCH3 and SCH3– is tested in modeling the Cys containing active site. In total we have compiled 51 zinc molecules. Geometries are compared to B3LYP/aug-cc-pVTZ optimized structures. B3LYP with the small basis set deviates on average by only 0.005 Å, and the maximal deviation is only 0.063 Å (Table 13). The MAD for metal–ligand length from DFTB3/3OB optimized structures is 0.024 Å. MAD of metal–ligand lengths Zn–O is 0.033 Å and Zn–S is 0.027 Å; they appear to be worse than the small test set but are good enough for modeling enzyme active sites. The largest deviations for metal–ligand lengths are found in Zn–O distances in [Zn(Im)(H2O)(SCH3)2] and [Zn(H2O)(SHCH3)(SCH3)2], which are underestimated by 0.141 and 0.103 Å, respectively. All others errors are smaller than 0.1 Å. The largest angle deviation 25.8° occurs in angle Zn–N–C in [Zn(C3N2H3)]+ (i.e., Zn2+ bound to a single deprotonated imidazole), which is not encountered often in biological systems.

Table 13

Mean and Maximum Absolute Deviation of Metal–Ligand Lengths and Angles in Comparison to B3LYP/aug-cc-pVTZ for the Zinc Large Test Set

propertya	Nb	B3LYPc	PM6	DFTB3/MIO	DFTB3/3OB
r (Å)	173	0.005	0.128	0.053	0.024
rmax (Å)		0.063	0.955	0.177	0.141
a (deg)	236	0.6	5.3	4.0	3.3
amax (deg)		5.1	63.3	25.8	25.8

Max stands for maximum absolute deviation.

Number of comparisons.

Basis set 6-31+G(d,p).

Number of pan class="Chemical">comclass="Chemical">parisons.

Ligand dissociation energies are class="Chemical">comclass="Chemical">piled in Table 14 (for detailed data, see Suclass="Chemical">pclass="Chemical">porting Information). class="Chemical">pan class="Chemical">DFTB3 shows that removal of neutral imidazole may cause an error of more than 10 kcal/mol, indicating interactions between Zn–sp2 hybridized N are problematic. Nevertheless, when calculating single-point energies using B3LYP/aug-cc-pVTZ the errors turn out to be small with the MAD of merely 0.5 kcal/mol. Proton affinities are also compared in Table 14. DFTB3/3OB shows excellent results with a MAD of 3.3 kcal/mol, which is fairly remarkable in this set of rather complicated molecules; single-point B3LYP calculations lead to an even smaller MAD of 0.6 kcal/mol.

Table 14

Error Statistics for the Ligand Dissociation Energies and Proton Affinities of Zinc Containing Large Molecules Compared to B3LYP/aug-cc-pVTZa

	B3LYPb	PM6	DFTB3/MIO	DFTB3/3OB	B3LYP//DFTB3c
Ligand Dissociation	Energies
MAD	1.3	5.9	3.6	3.4	0.5
MAX	2.6	20.2	14.4	13.8	1.3
Proton Affinities
MAD	1.1	10.1	4.9	3.3	0.6
MAX	2.1	19.3	12.7	6.9	2.6

For detailed data, see Supporting Information.

Basis set 6-31+G(d,p).

B3LYP/aug-cc-PVTZ single-point calculations on top of DFTB3/3OB geometries.

B3LYP/aug-cc-PVTZ single-point calculations on top of pan class="Chemical">DFTB3/3OB geometries.

Dinuclear Mg/Zn Active Sites

An important catalytic motif in class="Chemical">metalloenzymes is the biclass="Chemical">pan class="Chemical">metallic (or dinuclear) Mg/Zn site.[135,136] In this bimetallic motif, metal ions are commonly coordinated through a monatomic bridging ligand like water/hydroxide or a polyatomic group like carboxylate in a bidentate fashion. The two-metal-ion mechanism[5] was proposed for this type of motif, where one metal ion actively participates in facilitating the dissociation of the leaving group and stabilizing the pentavalent intermediate while the second metal ion holds an important structural role during catalysis. The advantage of bimetallic centers is the charge delocalization so as to lower the barrier of binding large substrates and generating a nucleophile.[137] To assess DFTB3/3OB for these important catalytic motifs, we have constructed four model compounds based on bimetallic Mg/Zn centers in enzymes (Figure 2). In our models, certain ligands are protonated to prevent a high net charge in the gas phase.

Figure 2

Optimized structures of dinuclear Mg2+/Zn2+ complexes of the simplified protein active site (a) PDB ID 1A49; (b) PDB ID 4ECQ; (c) PDB ID 3L4K; (d) PDB ID 1AMP. Distance are given in angstroms. The numbers without parentheses are obtained at the B3LYP 6-31+G(d,p) level; those with parentheses are DFTB3 values.

Figure 2(a) is a simplified active site structure taken from class="Chemical">pyruvate kinase crystal structure (class="Chemical">pan class="Disease">PDB ID 1A49(138)). Pyruvate kinase is an enzyme of the glycolytic pathway and phosphorylates ADP to generate ATP and pyruvate. This enzyme can also utilize other divalent metal ions like Mn2+, Ni2+, and Zn2+. Figure 2(b) is simplified from the polymerase η·DNA·dATP complex crystal structure (PDB ID 4ECQ(139)) which is involved in DNA synthesis. The formation of a phosphodiester bond is recently captured by time-resolved X-ray crystallography.[139] In the (a)(b) active sites, one Mg2+ is tridentated to α-, β-, and γ-oxygens in ATP. Figure 2(c) is modified from the type II topoisomerase (topoII)–DNA cleavage complex crystal structure (PDB ID 3L4K(140)). This enzyme cleaves and ligates DNA using tyrosine as the nucleophilic agent. Following the nucleophilic attack toward DNA backbones, a covalent phosphotyrosyl bond that links topoII and partial DNA is formed. It was long believed that a catalytic reaction could be aided via a single Mg2+. More recently, it was found that the canonical two-metal-ion mechanism may be possible.[140,141] Therefore, two Mg2+ ions are modeled in the reactant state. Figure 2(d) is the active site in Aeromonas Protelytica Aminopeptidase (AAP) (PDB ID 1AMP(142)), which is a prototypical member of the bizinc enzyme family.

As shown in Figure 2, pan class="Chemical">DFTB3 is tested against B3LYP/6-31+G(d,class="Chemical">p) oclass="Chemical">ptimized structures. The overall structures from class="Chemical">pan class="Chemical">DFTB3 reproduce B3LYP structures very well, particularly the metal ion–ligand distances. Mg–Mg/Zn–Zn distances are also excellent with the largest deviation of 0.13 Å. The monodentate or bridging binding modes of carboxylate are correctly described by the 3OB parameter set in comparison to B3LYP. These calculations demonstrate that DFTB3/3OB is useful for modeling bimetallic centers in enzymes with good accuracy.

Optimized structures of dinuclear pan class="Chemical">Mg2+/class="Chemical">pan class="Chemical">Zn2+ complexes of the simplified protein active site (a) PDB ID 1A49; (b) PDB ID 4ECQ; (c) PDB ID 3L4K; (d) PDB ID 1AMP. Distance are given in angstroms. The numbers without parentheses are obtained at the B3LYP 6-31+G(d,p) level; those with parentheses are DFTB3 values.

To further benchmark class="Chemical">DFTB3/3OB in describing the transition state and reaction energy barrier, we study an active site model of AP in the gas class="Chemical">phase. AP catalyzes the hydrolytic reactions of various class="Chemical">pan class="Chemical">phosphates via a two-step mechanism: an oxygen neocleophile first attacks the phosphorus, then a water (hydroxide) replaces the leaving group in the subsequent step that is essentially the reverse of the first step.[143] In our previous work,[48] we studied the first step for the hydrolysis of a well-studied phosphate diester MpNPP– in R166S AP and showed that the transition states are synchronous in nature. To help understand the intrinsic errors in DFTB3/3OB for the reactions of interest, we investigate the active site model in the gas phase with QM calculations. Specifically we study the same reaction for two additional phosphate diesters methyl 3-nitrophenyl phosphate (MmNPP–) and methyl phenyl phosphate (MPP–) and compare the energetics and nature of the transition state to MpNPP–. For simplicity, we only study the α orientation of different substrates.

As shown in Figure 3, there is generally good agreement between class="Chemical">DFTB3/3OB and B3LYP/6-31+G(d,class="Chemical">p) results in terms of geometries; the P–O distances, class="Chemical">pan class="Chemical">Zn2+–Zn2+ distances, and coordination structures are similar at the two QM computational levels for both reactant and transition states. The Zn–O distance error is within 0.03 Å and is a major improvement over MIO parameters which systematically overestimate Zn–O distance by 0.1 Å (MIO data not shown here). The Zn2+–Zn2+ distance is underestimated in DFTB3/3OB by ∼0.1 Å in the reactant and is shrunk ∼0.1 Å more in the transition state in comparison to B3LYP distances. In terms of the nature of the transition state, both DFTB3 and B3LYP predict a synchronous type of structure with the tightness coordinate (i.e., P–Oligand + P–Onucleophile) of 4.0–4.2 Å even though the individual P–O bond is overestimated or underestimated by 0.1 Å in DFTB3, indicating DFTB3 provides a good description for the potential energy surface around the transition state for the three substrates.

Figure 3

AP active site gas-phase model with MpNPP–, MmNPP–, and MPP–. Distances are given in angstroms. The numbers without parentheses are obtained at the B3LYP 6-31+G(d,p) level; those in the parentheses are DFTB3 values. (a) MpNPP– reactant state; (b) MpNPP– transition state; (c) MmNPP– reactant state; (d) MmNPP– transition state; (e) MPP– reactant state; and (f) MPP– transition state.

Regarding the reaction energy barrier, Table 15 shows that the qualitative trend is class="Chemical">consistent among various class="Chemical">pan class="Chemical">computational levels although there are substantial differences in the quantitative values. Previous calculations[48] of proton affinities for the leaving group indicate that MP2 gives results close to the experimental values, so MP2/6-311++G(d,p) single-point energies are used as the reference here. Numbers span a broad range with different DFT methods without dispersion, and the results become much closer after adding the D3 dispersion; an explicit consideration of dispersion is important for a large active site discussed here. Compared to MP2, DFTB3 gives a reasonable reaction energy for MpNPP– but underestimates the barrier of MmNPP– by ∼4 kcal/mol and overestimates the barrier of MPP– by ∼9 kcal/mol. However, MP2 single-point energies on top of DFTB3 and B3LYP/6-31+G(d,p) geometries differ by ∼1 kcal/mol for all the substrates, indicating DFTB3 produces satisfactory geometries and also highlighting the importance of single-point energy correction with a reliable ab initio method at DFTB3 geometries. DFTB3 reaction barriers are close to B3LYP results without D3 correction; this is not surprising because DFTB3 is parametrized based on B3LYP (see Methods) and inherits some errors. M06 behaves the best among three DFT functionals as compared to MP2, especially with the D3 dispersion included. PBE systematically underestimates the energy barrier. Further optimization of the D3 dispersion model developed for DFTB3[99] will further improve DFTB3/3OB energies.

Table 15

Calculated Activation Barrier (in kcal/mol) for Diester Hydrolysis in the Active Site Model of AP in the Gas Phasea

		single points at DFTB3 geometries				single points at B3LYP geometries
substrate	DFTB3	B3LYP	M06	PBE	MP2	B3LYP	M06	PBE	MP2
MpNPP–	8.8	9.2 (4.2b)	5.7 (5.1)	5.7 (3.1)	7.6	12.0 (6.7)	6.8 (6.1)	7.2 (4.3)	8.8
MmNPP–	9.7	9.2 (9.4)	12.1 (12.4)	6.8 (7.5)	13.6	15.3 (12.0)	13.2 (12.8)	11.0 (9.3)	13.2
MPP–	20.0	16.2 (9.1)	12.4 (10.5)	11.6 (7.1)	9.6	20.1 (11.6)	13.0 (10.8)	14.0 (8.6)	10.4

The geometries are optimized at either the DFTB3/3OB or B3LYP/6-31+G(d,p) level. Single-point energy calculations are carried out at the B3LYP, M06, PBE, and MP2 method at two levels of geometries using a large basis set of 6-311++G(d,p).

The numbers with parentheses are obtained with corresponding DFT functional plus D3 dispersion.

AP active site gas-phase model with pan class="Chemical">MpNPP–, class="Chemical">pan class="Chemical">MmNPP–, and MPP–. Distances are given in angstroms. The numbers without parentheses are obtained at the B3LYP 6-31+G(d,p) level; those in the parentheses are DFTB3 values. (a) MpNPP– reactant state; (b) MpNPP– transition state; (c) MmNPP– reactant state; (d) MmNPP– transition state; (e) MPP– reactant state; and (f) MPP– transition state.

The geometries are optimized at either the pan class="Chemical">DFTB3/3OB or B3LYP/6-31+G(d,class="Chemical">p) level. Single-class="Chemical">point energy calculations are carried out at the B3LYP, M06, class="Chemical">pan class="Chemical">PBE, and MP2 method at two levels of geometries using a large basis set of 6-311++G(d,p).

The numbers with parentheses are obtained with pan class="Chemical">corresclass="Chemical">ponding DFT functional class="Chemical">plus D3 disclass="Chemical">persion.

Mg2+ and Zn2+ Ions Solvation Structure and Relative pKa

The class="Chemical">metal ion-O radial distribution function (g(r)) and class="Chemical">pan class="Chemical">coordination numbers are shown in Figure 4 for the two cations in solution. A well-defined first solvation shell is observed in both cases. The first peak in DFTB3/MM hydrated Mg2+ occurs at 2.07 Å, in good agreement with X-ray diffraction data of 2.09 Å ± 0.04 Å[144] and with previous AIMD simulations in the range of 2.08–2.13 Å.[145−149] The MM force field predicts the Mg2+–O distance to be shorter than the DFTB3/MM model with a sharper peak at 2.0 Å. The inner coordination sphere remains similar to the gas-phase structure. For the cluster [Mg(H2O)6]2+ in the gas phase, the Mg2+–O distance is 2.10 Å in both DFTB3 and B3LYP/aug-cc-pVTZ optimized structures. For Zn2+, the first maximum from DFTB3/MM is at 2.11 Å, consistent with X-ray absorption measurements 2.06 ± 0.02 Å[150] and earlier AIMD results 2.07–2.13 Å.[145,146,151] The classical MM model gives the first peak at 2.08 Å, although the peak is higher and narrower than in the DFTB3/MM model. The small difference in the location of the first peak in DFTB3/MM may be traced back to the gas-phase model: for the cluster [Zn(H2O)6]2+ in the gas phase, Zn2+–O is 2.16 Å in DFTB3 versus 2.12 Å in B3LYP/aug-cc-pVTZ.

Figure 4

Radial distribution function of water oxygen around Mg2+/Zn2+ in 20 Å water droplet. The pure MM model is based on the CHARMM force field.[105]

The integral of the first solvation shell for both cations gives the class="Chemical">coordination number of 6. The lclass="Chemical">pan class="Chemical">arge QM region results predict very similar first solvation shell as with the small QM region, though they differ slightly for the second solvation shell; the large QM region leads to a broader peak for the second solvation shell than the small QM region and MM model. It is worth noting that in the large QM region simulation of Zn2+ water exchanges between the first and second solvation shells (coordination number remains the same), whereas no exchange of water is seen in Mg2+ simulations. These observations are in agreement with experimental findings; for Mg2+, experimental results suggest that water molecules could reside in the first shell up to a few microseconds,[152,153] while the water residence time for Zn2+ is only hundreds of picoseconds.[154]

To further benchmark the energetics of class="Chemical">DFTB3/3OB, we use the class="Chemical">pan class="Chemical">DTSC-TI approach to calculate the relative pKa of Mg2+/Zn2+-bound water with the DFTB3/MM model. Altogether six intermediate states are used to obtain the ∂F/∂λ values. In each window, the data collected during the first 250 ps are discarded as part of the equilibration process, and the remaining 500 ps data are then subject to statistical analysis to obtain the proper average of free energy derivative. The statistical errors associated with the free energy are on the order of 0.2–0.5 kcal/mol. In all cases, the linear dependence of the free energy derivatives on λ holds well with a correlation coefficient (R2) higher than 0.99 (see Figure S1 in Supporting Information). Integrating ∂F/∂λ over λ gives a value of 132.9 kcal/mol for Mg2+ with the small QM region, 124.1 kcal/mol for Zn2+ with the small QM region, 130.0 kcal/mol for Mg2+ with the large QM region, and 125.5 kcal/mol for Zn2+ with the large QM region, leading to a ΔpKa of 6.4 pKa units with the small QM region and 3.3 pKa units with the large QM region. The latter is in good agreement with the experimental value of 3.2 pKa units (12.4 for a Mg2+ solution and 9.2 for hydrated Zn2+).[145,155] The better performance of the large QM region results highlights the importance of using a QM description for the second solvation shell, which allows coordination number fluctuation in the Zn2+ solution, especially at the deprotonated state. We have also calculated the absolute pKa of Mg2+/Zn2+-bound water as a further benchmark. Compared to the experimental values, the closest computational estimate from current work is 14.0 for Mg2+ and 10.4 for Zn2+ (see Supporting Information).

Radial distribution function of pan class="Chemical">water class="Chemical">pan class="Chemical">oxygen around Mg2+/Zn2+ in 20 Å water droplet. The pure MM model is based on the CHARMM force field.[105]

Structural Flexibility of the Active Site in Myosin

Myosin-II (thereafter referred to as myosin) is one of the best characterized molecular motors, which class="Chemical">couclass="Chemical">ple lclass="Chemical">pan class="Chemical">arge-scale conformation transitions with ATP binding and hydrolysis. Despite many previous efforts from both experimental and theoretical studies,[156−160] the complex nature of “mechanochemical coupling” still remains to be well understood.

As shown in Figure 5(a), binding modes obtained from the class="Chemical">DFTB3/MM simulation reclass="Chemical">produce the key features of the crystal structure. In the discussion of our earlier work,[115] the weakened interaction between class="Chemical">pan class="Chemical">Mg2+-nucleotide and Ser237 as reflected by the larger distance fluctuation was proposed to indicate the destabilization of Switch I following ATP hydrolysis. As a test of the DFTB3/3OB magnesium parameters in the enzyme, DFTB3/MM equilibrium simulations are carried out with two nucleotide chemical states (ATP or ADP·Pi), and the distances between Mg2+ and four nonwater ligands are monitored. As shown in Figures 5(b) and (c), the distances between Mg2+ and ligands are reasonably stable over 1 ns; the most obvious difference observed between the two nucleotide states is the fluctuation of Mg2+–OSer237 distance, for which the average and fluctuation are 2.16 and 0.08 Å in the ATP state and 2.24 and 0.11 Å in the ADP·Pi state. The magnitude of fluctuation observed in current DFTB3/MM simulations is notably less compared with previous work (MM simulations and SCC-DFTBPR/MM simulation, in which the Mg–OSer237 distance reached ∼3.5 Å), which led to coordination number shifting from 6 to 5 after ATP hydrolysis.[115] To further evaluate the situation, PMFs are calculated for the distance between Mg2+ and OSer237 in the ATP and ADP·Pi states. As seen in Figure 5(d), dissociation of Ser237 is indeed energetically more favorable in the ADP·Pi state than the ATP state by a few kcal/mol. Nevertheless, there is a significant energetic cost to fully dissociate the Ser237 side chain from the Mg2+, even in the ADP·Pi state. Further experimental analysis based on infrared spectroscopy can potentially generate additional insights into this issue.[161]

Structural properties of the Myosin active site during equilibrium QM/MM MD simulations with different nucleotide chemical states (class="Chemical">ATP or class="Chemical">pan class="Chemical">ADP·Pi). (a) Overlay of average DFTB3/MM structure and crystal structure (gray) of Myosin-II motor domain active site in the ATP state. The VO4 moiety in the crystal structure (PDB code 1VOM) is replaced by a PO3 phosphate group and a presumptive lytic water molecule. Instantaneous distances between Mg2+ and oxygen atoms in its four nonwater ligands at (b) ATP state and (c) ADP·Pi state. Ser237 refers to OSer237 and Thr186 refers to OThr186. (d) PMF calculations for the ATP/ADP·Pi states. The reaction coordinate is the distance between Mg2+ and OSer237.

Force Fields Comparison in AP Active Site

We have benchmarked class="Chemical">DFTB3/3OB for the AP active site with a gas-class="Chemical">phase model in earlier discussion. As an additional test, class="Chemical">pan class="Chemical">DFTB3/MM simulations for the R166S AP enzyme are carried out, and the results are compared to two popular MM force fields for zinc. The comparison of equilibrated reactant structures by Coulomb,[14] SLEF1,[25] and DFTB3/MM is shown in Figure 6, with the crystal structure with bound phosphate as a reference. Binding modes obtained from the DFTB3/MM simulation reproduce the key features of the crystal structure and previous SCC-DFTBPR/MM simulations[48] very well (shown in Figure 6(c) and (d)). The conventional Coulomb scheme (Figure 6(a)) leads to a much shortened Zn2+–Zn2+ distance and significant changes in the active site geometry: Ser102 and Asp327 are stabilized by two zinc ions simultaneously; Asp369 (not shown in Figure 6 for clarity) is bidentate with the Zn12+ in the Coulomb model but monodentate in the crystal and DFTB3/MM structures. By comparison, coordination geometry is notably improved by the SLEF1 force field (Figure 6(b)), even though the Zn2+–Zn2+ (3.7 Å) distance is still substantially shorter than that in the crystal structure (4.3 Å). It is perhaps not surprising that SLEF1 does not provide a correct distance between zinc ions as SLEF1 was not parametrized for the binuclear zinc motif and designed to be compatible with the Amber99SB force field rather than with the CHARMM force field. Nevertheless, SLEF1 correctly describes the binding modes of Ser102 and Asp369 to the zinc ions; yet Asp327 still has a strong preference to bridge the two zinc cations due to the short distance between them. Moreover, an extra His372 is pulled closer to bind to Zn22+; this might have been caused by the stronger attractive interaction between the zinc ion and nitrogen in the short-range part of the SLEF1 force field. In short, DFTB3/MM produces encouraging results for modeling a challenging bimetallic zinc active site, demonstrating the greater applicability of DFTB3/MM to metalloenzymes even when only structural properties are considered.

Structural properties of the AP active site during equilibrium QM/MM MD Simulations. A snapshot for the reactant state with (a) class="Chemical">Coulomb MM force field, (b) SLEF1 zinc MM force field, (c) class="Chemical">pan class="Chemical">DFTB3/3OB, and (d) overlay of crystal and DFTB3/MM structure (blue), with average key distances in angstroms labeled. Asp369, His370, and His412 are omitted for clarity. Substrate inorganic phosphate in the crystal structure (3CMR) is not shown.

Concluding Remarks

pan class="Chemical">Magnesium and zinc class="Chemical">play many imclass="Chemical">portant roles in class="Chemical">pan class="Chemical">metalloenzymes and various chemical systems. An efficient and accurate QM description for Mg/Zn is critical for the use of theoretical methods to study the reaction mechanism of these systems. In this work, we extend the parametrization of the approximate density functional tight binding theory, DFTB3, to magnesium and zinc. The parametrization is performed in a framework consistent with the DFTB3/3OB set for CHNOPS; hence this parameter set can be used to study the chemistry of Mg/Zn moieties in solution and biological systems.

Benchmark calculations in the gas phase and in the class="Chemical">condensed class="Chemical">phase with a QM/MM framework demonstrate that the current class="Chemical">parametrization generally leads to reliable structures and semiquantitative energetics (with tyclass="Chemical">pically a MAD of ∼3–5 kcal/mol class="Chemical">pan class="Chemical">compared to high level ab initio calculations). In the gas phase, DFTB3/3OB are tested with a diverse range of molecules that reflect the typical coordination environment of magnesium and zinc in biomolecules, including dinuclear metal sites. We focus on the geometries, ligand dissociation energies, and ligand proton affinities. The results calculated at the DFTB3/3OB level are compared to B3LYP, ab initio (G3B3, MP2), and a popular semiempirical method PM6 as well as the older parametrization, MIO. In general, DFTB3/3OB shows substantial improvement over DFTB3/MIO and outperforms PM6 in many aspects within our test sets, particularly in terms of structural properties and ligand proton affinities. Compared to high-level calculations, DFTB3/3OB is very successful at predicting geometries. Large errors can still be found in certain cases, especially with charged ligands. Even for these cases, single-point energy calculations with high level QM methods generally give very reliable energetics, suggesting that changes of structures during typical processes of interest (e.g., ligation dissociation or ligand deprotonation) are well captured with DFTB3/3OB.

Due to class="Chemical">comclass="Chemical">putational efficiency and general robustness for geometry class="Chemical">predictions, class="Chemical">pan class="Chemical">DFTB3/MM calculations are expected to be effective for studies in condensed-phase systems. For example, our solution benchmark study for the magnesium and zinc ion has reproduced experimental solvation structures and relative pKa of metal bound water, especially when both the first and second solvation shells are treated as DFTB3/3OB. For the enzyme benchmark, DFTB3/3OB is successful at producing active site structures in myosin and AP relative to available crystal structures; for AP, DFTB3/3OB also leads to semiquantitative energetics for an active site model, although our calculations using several DFT functionals and MP2 highlight the importance of careful calibration of energetics for large metal sites in enzymes, including consideration of dispersion interactions. These studies have laid the groundwork for using DFTB3/MM simulations with state-of-the-art sampling techniques such as metadynamics[162] and finite-temperature string methods[163−165] to study the chemical processes in myosin and AP; high-level single-point QM/MM calculations are still required for improving the quantitative nature of the free energy surfaces. Finally, we emphasize that there are still remaining limitations in DFTB3, such as the treatment of interaction between Mg/Zn and highly charged/polarizable ligands (e.g., hydroxide and SH–), thus it is essential to continue the formal development of DFTB3, such as adding multipole terms for the charge fluctuations and formulating a better description of polarization and short-range Pauli repulsions.