Polarizable Multipole-Based Force Field for Dimethyl and Trimethyl Phosphate.

Phosphate groups are commonly observed in biomolecules such as nucleic acids and lipids. Due to their highly charged and polarizable nature, modeling these compounds with classical force fields is challenging. Using quantum mechanical studies and liquid-phase simulations, the AMOEBA force field for dimethyl phosphate (DMP) ion and trimethyl phosphate (TMP) has been developed. On the basis of ab initio calculations, it was found that ion binding and the solution environment significantly impact both the molecular geometry and the energy differences between conformations. Atomic multipole moments are derived from MP2/cc-pVQZ calculations of methyl phosphates at several conformations with their chemical environments taken into account. Many-body polarization is handled via a Thole-style induction model using distributed atomic polarizabilities. van der Waals parameters of phosphate and oxygen atoms are determined by fitting to the quantum mechanical interaction energy curves for water with DMP or TMP. Additional stretch-torsion and angle-torsion coupling terms were introduced in order to capture asymmetry in P-O bond lengths and angles due to the generalized anomeric effect. The resulting force field for DMP and TMP is able to accurately describe both the molecular structure and conformational energy surface, including bond and angle variations with conformation, as well as interaction of both species with water and metal ions. The force field was further validated for TMP in the condensed phase by computing hydration free energy, liquid density, and heat of vaporization. The polarization behavior between liquid TMP and TMP in water is drastically different.

Introduction

The backbone building blocks of the genetic materials DNA and RNA contain ionic phosphate groups. Due, in part, to their negative charge, nucleic acid molecules can be retained within a pan class="Chemical">lipid membrane and their phosphodiester bonds are very stable against hydrolysis.[1] The tertiary structure and flexibility of DNA and RNA, which are central to their functions, also stems from the rotation of the phosphodiester bonds along the backbone. To study the structural and energetic properties of the backbone of DNA/RNA, the DMP (dimethyl phosphate) anion containing the same phosphodiester linkage as that in DNA/RNA has often been employed as a simple model compound.[2] DMP has also been a popular anion for ionic liquids.[3] TMP (trimethyl phosphate), which has three phosphoester bonds, is a neutral molecule and a liquid at room temperature. TMP finds use as a solvent[4] and as a mild methylating agent.[5]

There are three dominant conformations for both negatively charged DMP and neutral pan class="Chemical">TMP, as depicted in Figure . It is difficult to accurately describe the electrostatic potential around all conformations of these molecules with a single set of atomic partial charges. AMOEBA[6,7] utilizes atomic permanent electrostatic multipole moments through the quadrupole, which we have shown can accurately model the electrostatic potential around various peptide conformations.[8] In addition, many-body polarization effects are explicitly treated with atomic dipole induction. Phosphorus, located in period 3 of the periodic table, is larger and softer than the elements from period 2 and is even more polarizable. In the AMOEBA force field, molecular polarizability is modeled via a Thole-style[9] damped interactive induction model based upon distributed atomic polarizabilities.

Figure 1

Minimum energy conformations of DMP and TMP monomers.

Different secondary or tertiary conformations of DNA/RNA are formed by the rotation of phosphodiester linkages of the backbone, and incorrect nucleic acid torsional parameters may result in significant structural distortion.[10] The three conformations of DMP and TMP (Figure ) have been elucidated and investigated in experimental and quantum mechanical studies.[11,12] Theoretical studies also show that solvent and metal ions may affect the geometry and transition dynamics among different conformations of DMP.[13−15] There are large periodic variations in bond lengths and bond angles around the phosphate O–P–O linkage as a function of phosphodiester torsional rotation. These structural changes in DMP and TMP are exactly analogous to the well-known anomeric effect seen in carbohydrates. Pinto et al. have provided ab initio calculations and a perturbational molecular orbital framework that extends the anomeric effect for bond lengths to account for angle changes.[16]

The AMOEBA polarizable force field for water,[6,17] organic molecules,[18] peptides, and proteins[8] has been developed previously. In this work, as a first step toward generating a polarizable nucleic acid force field for biomolecular simulations, we report the development of AMOEBA models for pan class="Chemical">DMP and TMP based on comprehensive quantum mechanical studies of the molecular properties of DMP and TMP. We present the complete conformational energy surface map for both DMP and TMP, including stable conformations and their interconversion pathways. High-level quantum mechanical (QM) methods were used to further study the geometry and potential energy of DMP and TMP in different environments (e.g., in solution or bound to metal ions). The molecular electrostatic potential, vibration frequencies, polarizability, and interaction energy with water of DMP and TMP were also examined with QM calculations. On the basis of these results, the AMOEBA force field for DMP and TMP was developed. To validate it, the predicted liquid properties of TMP, including density and heat of vaporization, as well as the TMP hydration free energy were compared with reported experimental data.

Computational Methods

All ab initio QM calculations were performed with Gaussian 09,[19] using basis sets as specified in the relevant sections below. The polarizable continuum model (PCM)[20] was applied to introduce solvent effects into the QM calculations. For computation of pan class="Chemical">TMP and DMP conformational energies (Figures , 8, S2, and S3), structures were first optimized using MP2/cc-pVTZ with PCM, followed by a single-point energy calculation at the MP2/aug-cc-pVTZ level, with or without PCM. Atomic multipole moments for DMP or TMP were initially assigned from QM electron density calculated at the MP2/6-311G** level and using Stone’s distributed multipole analysis (GDMA v2.2).[21] The Switch 0 and Radius H 0.65 options were used with GDMA to access the original DMA procedure.

Figure 2

DMP and TMP QM conformational energy maps in implicit solvent. (A, C) For DMP, the two Os–P–Os–C torsion angles (named χ1 and χ2) were sampled every 10°, and energy maps on a 37 × 37 grid were computed. (B, D) For TMP, the three O–P–Os–C torsion angles (named ψ1, ψ2, and ψ3) were sampled every 45° from 0° to 360°, and a 9 × 9 × 9 grid was constructed. Three slices of the 3D map at ψ3 = 45°, 180°, and 315° and the contour surface at 1.3 kcal/mol level are presented. The potential wells of the two maps (C, D) are labeled in white. The possible transition paths for DMP are marked by dashed lines.

Figure 8

Comparison of conformational energy surfaces, including implicit solvation, calculated by QM and AMOEBA for DMP and TMP. (A) DMP χ1–χ2 2D potential energy map. (B) A slice with ψ3 = 45° from the TMP 3D conformational energy map.

All force field-based calculations were performed using TINKER 6.[22] The VALENCE program in TINKER was used to derive initial valence force field parameters based on QM-optimized structure and frequencies. The TINKER POLARIZE program was used to compute molecular polarizabilities based on atomic polarizability parameters. The TINKER MINIMIZE program was used for structure optimization, with the convergence criteria set to 0.001 kcal/mol/Å. The TINKER POTENTIAL program was used to obtain electrostatic potentials (ESP) around a molecule from Gaussian 09 cube files or based on AMOEBA multipole parameters. POTENTIAL was also used to optimize the electrostatic parameters against QM electrostatic potential grids (MP2/aug-cc-pVQZ). Only atomic dipole and quadrupole moments were allowed to deviate from the DMA values during the ESP optimization, which stops when a gradient of 0.1 kcal/mol/electron[2] was achieved for all conformations included for either TMP or DMP. The generalized Kirkwood (GK)[23] implicit solvent model was used to account for solvent effects as part of AMOEBA structure optimization and conformational energy calculations. AMOEBA force field parameters for water[6] and metal ions[24,25] from previous studies were used here. The TINKER ANALYZE program was used to calculate the total potential energy and the energy components as well as the multipole moments.

The nonlinear optimization method lsqnonlin in Matlab[26] was used to fit the torsion parameters to the conformational potential energy surface and to optimize the torsion, stretch-torsion, and angle-torsion parameters in the final refinement step against conformational energies of the dominant configurations.

Molecular dynamics (MD) simulations were carried out using TINKER 6.[22] The Bussi–Parrinello thermostat[27] and RESPA integration method were used in all simulations. Spherical energy cutoffs for van der Waals and Ewald direct-sum were 12.0 and 8.0 Å, respectively. Default PME cutoff and grid sizes were used in the reciprocal space.

The TMP liquid simulations were performed at four temperatures, 227 K (the melting temperature), 298.15 K, 373 K, and 470.35 K (the boiling temperature). For each temperature, a cubic box containing 300 pan class="Chemical">TMP molecules (initial conformations were all C3) was built. The initial box sizes were set to match the initial density to the corresponding experimental value. All systems were relaxed via energy minimization and then simulated for 2 ns in the NVT ensemble. NPT simulations, utilizing a Berendsen barostat,[28] were performed for 4 ns at each of the four temperatures and at a target pressure of 1 atm to evaluate the density and heat of vaporization of TMP. After discarding 1 ns of the simulation for equilibration, the final 3 ns of each trajectory was used for the calculation of the liquid properties of TMP.

Molecular dynamics (MD) simulations of gas-phase TMP at 470.35 K (boiling point) were carried out for heat of vaporization calculation. NVT simulations of 10 TMP molecules in a box of 62.5 × 62.5 × 62.5 Å3 were performed for 2 ns.

The Bennett acceptance ratio (BAR) method was used for TMP hydration free energy calculation, using a protocol similar to that in previous studies.[29,30] A pan class="Chemical">TMP molecule was inserted into a 40 × 40 × 40 Å3 cubic box containing 2138 water molecules. After system equilibration using NPT simulations at 293.15 K, the electrostatic interaction between TMP and water molecules was turned off at a constant interval of λ = 0.1, where λ is the interaction strength. A total of 11 NVT simulations (293.15 K) were carried out for each λ value for 2 ns. Similarly, the soft-core van der Waals (vdW) interactions were weakened in a series of steps: 1.0, 0.9, 0.8, 0.75, 0.7, 0.65, 0.6, 0.55, 0.5, 0.4, 0.3, 0.2, 0.1, 0.0. In the gas-phase recharging process, the λ interval was set to 0.1.

Potential Energy Model and Parameterization Strategy

The AMOEBA potential energy function (eq ) comprises valence and nonbonded interaction terms. The former includes van der Waals, permanent, and induced electrostatic interactions. Bond-stretching, angle-bending, and torsional energy are the traditional bonded terms. New stretch-torsion and angle-torsion coupling terms have been added to accurately model the structure and energy of pan class="Chemical">DMP and TMP in different conformations.

Both DMP and TMP are represented by five types of atoms: phosphorus (P), double-bonded oxygen (O), single-bonded oxygen (Os), methyl carbon (C), and methyl hydrogen (H).

As with previous force field development for small organic molecules,[31] the parameters for each potential term were determined step-by-step (illustrated in Scheme S1, Supporting Information) based mainly on QM calculations of various properties using Gaussian 09.[19] First, the structure for each molecular conformation was obtained by high-level QM optimization in gas and solution (implicit solvent) environments. On the basis of these structures, the molecular polarizability, electrostatic potential field, and vibration frequencies were calculated by QM. The electrostatic parameters were derived from distributed multipole analysis (DMA)[21] and then optimized to the QM electrostatic potential surface (step 2 in Scheme S1), a procedure described in detail previously.[8] Similar to the work of Kramer et al.,[32−34] we fit multiple conformers to the ESP simultaneously to ensure conformational transferability. Atomic polarizabilities were initially transferred from previous Thole-type models,[9] with values for atoms P and O fine-tuned to match the QM-derived molecular polarizability tensor. Bond and angle parameters were determined by fitting to QM structures and vibrational frequencies (step 3). The QM interaction energy curves were computed for DMP–water and TMP–water and used to optimize the van der Waals parameters (step 4). Next, torsional potential parameters were developed to reproduce the QM conformational energy surface (step 5), and stretch-torsion and angle-torsion coupling terms were parametrized to capture the variations of bond length and angle with torsions (step 6). Subsequently, all valence potential parameters (bond, angle, torsion, and the coupling terms) were refined via fitting to both the QM energies and structures (step 7). Lastly, as validation of the model, the final force field parameters were applied to compute the hydration free energy, liquid density, and heat of vaporization of TMP, as well as DMP–ion interaction structures and energetics.

Important Conformers of DMP and TMP

Minimum-energy geometries for the important stable conformations of DMP and pan class="Chemical">TMP were obtained from QM and compared with crystal structures.[35,36] The whole conformational energy surface for DMP and TMP was explored to locate minimum-energy conformations and their transition routes in solution. The two dimensions used for the DMP 2D conformational energy map were the two Os–P–Os–C torsions (named χ1 and χ2) (Figure A). The three dimensions for TMP 3D conformational map were the three O–P–Os–C torsions (named ψ1, ψ2, and ψ3) (Figure B). The three most stable conformations of DMP were gg, gt, and tt (as shown in Figure C). The transition state configurations between gg and gt (gg → gt) were near (70°,125°) or (−70°,−125°), and the transition states between gt and tt (gt → tt) were about (180°/125°), or (180°/−125°), with the two torsion angles interchangeable due to symmetry. The TMP 3D conformational map (Figure D) shows that the three most stable conformations are C3, C1, and Cs. The energy gaps between these three conformations are small and comparable to the thermal energy at room temperature (0.593 kcal/mol). Transition state configurations between C1 and C3 (C1 → C3) are approximately at (45°,120°,45°) or (−45°,−120°,−45°), and transition configurations between C1 and Cs (C1 → Cs) are near (0°,180°,45°).

Next, we investigate the effect of the chemical environment on the molecular geometry and energy of these important conformations using QM methods. Counterions, usually metal ions, are typically present in solutions of pan class="Chemical">DMP anion and nucleic acids. Using MP2/cc-pVQZ level QM optimization, the bond and angle values of DMP were calculated for different environments, including gas phase, with one water (Figure A), with three waters (Figure B), in solution (PCM[20]), and bound to different metal ions mediated by a water molecule (Figure C). The optimized structures show that more polar or ionic environments make the P–Os bond shorter, the P–O and Os–C bonds longer, and the Os–P–Os angle larger (Table ). The difference among gas phase, solution phase, and ionic environments is rather significant. Compared with the crystal structures (the last four lines in Table ), the optimized structures with metal ions comprise the best representation for DMP in the condensed phase. We note that the environment affects both the structures and the potential energies. Table lists the relative potential energy of conformations in the gas phase, solution (PCM), and bound to sodium and magnesium (Figure C). The solution and ionic environments tend to smooth the conformational energy surface, lowering the energy gaps and transition barriers between different conformations. For TMP, the geometry and potential energy differences between the gas and solution phases are less significant than those of DMP, but the overall trends are similar (Table S1). On the basis of the above analysis, solution-phase DMP conformations with bound sodium/water and solution-phase TMP conformations were used in force field parametrization.

Figure 3

Three different environments used in DMP geometry studies.

Table 1

Calculated Bond Length and Angle Values of the DMP gg Conformation in Different Environments and the Corresponding Values Sampled from Crystal Structuresa

	bonds/angles
environment	P–O (Å)	P–Os (Å)	Os–C (Å)	Os–P–Os (deg)	P–Os–C (deg)
gas	1.490	1.653	1.413	99.6	115.5
1 water molecule (without PCM)	1.492	1.642	1.414	100.3	115.7
3 water molecules (without PCM)	1.495	1.633	1.423	101.4	116.1
PCM	1.497	1.628	1.427	102.1	117.4
Na+/water (PCM)	1.500	1.619	1.428	102.8	117.4
Mg2+/water (PCM)	1.505	1.610	1.431	103.4	117.6
Ni2+/water (PCM)	1.508	1.603	1.434	104.0	117.8
[Ni (H2O)6]2+	1.495	1.589	1.443	99.4	118.9
[Ni(C2H8N2)3]2+	1.487	1.604	1.436	104.1	118.9
AMPI+	1.472	1.602	1.430	103.3	118.3
Z-DNA (PDB: 3P4J)	1.496/1.490	1.600/1.589	1.438	103.7	120.8

Structures of DMP and its complexes (illustrated in Figure ) were optimized using MP2/cc-pVQZ with or without PCM. The crystal structural values for DMP complexed with the cations hexa-aqua-nickel(II) ([Ni(H2O)6]2+), tris(ethylenediamine-N,N′)-nickel(II) ([Ni(C2H8N2)3]2+), and 1-allyl-1-methylpyrrolidinium (AMPI+) are listed after the QM values. The last row shows values from a Z-DNA crystal structure at 0.55 Å resolution (PDB ID 3P4J), where the phosphates of the DNA adopt a gt conformation and alternative bond values were measured and averaged.

Table 2

DMP Conformational Energy in Different Environmentsa

	conformation
environment	gg	gg → gt	gt	gt → tt	tt
gas	0	2.428	1.430	3.499	3.295
PCM	0	2.688	1.408	3.499	2.975
Na+/water (PCM)	0	2.636	1.127	3.063	2.374
Mg2+/water (PCM)	0	2.261	0.571	2.446	1.885

Energy is in units of kcal/mol. The three most stable conformations and the two transition state configurations are included. All configurations were optimized using MP2/cc-pVQZ with or without PCM, and the energies were calculated using the same method as that during optimization. For the two transition states, one Os–P–Os–C torsion angle was restrained to be 125° during the structural optimization.

Structures of DMP and its complexes (illustrated in Figure ) were optimized using pan class="Gene">MP2/cc-pVQZ with or without PCM. The crystal structural values for DMP complexed with the cations hexa-aqua-nickel(II) ([Ni(H2O)6]2+), tris(ethylenediamine-N,N′)-nickel(II) ([Ni(C2H8N2)3]2+), and 1-allyl-1-methylpyrrolidinium (AMPI+) are listed after the QM values. The last row shows values from a Z-DNA crystal structure at 0.55 Å resolution (PDB ID 3P4J), where the phosphates of the DNA adopt a gt conformation and alternative bond values were measured and averaged.

In constrast to the CHARMM[37] and Amber ff03[38] force fields, which use only one set of bonded parameters for phosphates in different charge states, we adopt two different sets of valence parameters for DMP and TMP. We found that the geometry of methyl phosphate ion (MP2–), dimethyl phosphate ion (DMP–), trimethyl phosphate (TMP), and dimethyl phosphoric acid (DMPH) is significantly different across the various charge states (Table S2). Phosphates with a higher charge have longer P–O and P–Os bonds but shorter Os–C bonds. In DNA and RNA, the diester bonds O3′–C3′ and O5′–C5′, similar to those of the DMP ion, are shorter and stronger than those of the neutral phosphate molecules and exhibit greater resistance to hydrolysis.[39]

Electrostatic Parameters

The electrostatic parameters for pan class="Chemical">TMP and DMP were determined from data for multiple minimum-energy conformations. AMOEBA electrostatic interactions comprise permanent and induced components. The permanent electrostatic energy is expressed as the sum of multipole–multipole interactions between atom pairs, excluding or scaling nearby through-bond pairs. For atomic dipole and quadrupole moments, a local frame is defined at each atom based upon the neighboring atoms. Three conventions, Z-then-X, Z-bisector, and Z-only, were used for the local frame definitions for DMP/TMP, as illustrated in Figure .[31] The phosphate in DMP has 2-fold symmetry; thus, the Z-bisector convention was used for the central phosphorus atom. The bisector of the Os–P–Os angle defines the Z axis, and the X axis falls on the Os–P–Os plane. The phosphate in TMP and the methyl carbons have 3-fold symmetry, and the Z-only convention is used, where the X and Y axes are on the plane perpendicular to the Z axis. Other atoms follow the Z-then-X convention. When possible a non-hydrogen-bonded atom is selected to define the Z axis, another neighboring non-hydrogen atom defines the ZX plane and the positive X direction.

Figure 4

Local frame definitions for (left) DMP phosphorus, (middle) single-bonded oxygen, and (right) TMP phosphorus.

Local frame definitions for (left) DMP phosphorus, (middle) single-bonded pan class="Chemical">oxygen, and (right) TMP phosphorus.

Electron density matrices for the DMP gt conformation and the pan class="Chemical">TMP C1 conformation were computed at the MP2/6-311G** level, from which the initial DMP and TMP atomic multipole moments were obtained using Stone’s distributed multipole analysis (see Computational Methods).[21] To ensure that the DMP force field is generally applicable, DMP conformations gg, gt, and tt in three environments (solution only, with Na+/water, and with Mg2+/water), a total of 9 conformations, were included in the DMP electrostatic potential (ESP) optimization. Stable conformations C3, C1, and Cs in a solution environment were employed for the TMP electrostatic potential ESP optimization. The electrostatic potential on a grid of points around each conformer is constructed from MP2/aug-cc-pVQZ calculations. Grid points of four shells (0.35 Å apart) were generated around the molecule with an offset of 1 Å from the van der Waals surface. Then, a single set of multipole parameters for DMP or TMP was determined by fitting to the electrostatic potential grids of multiple conformers simultaneously using the TINKER POTENTIAL program. The initial GDMA partial charge values were fixed during the ESP optimization.

In the AMOEBA model, induced electrostatic interactions are implemented via Thole’s damped induction approach.[9] The atomic polarizabilities of single-bonded pan class="Chemical">oxygen Os and methyl C and H in DMP and TMP were assigned typical values taken from other molecules in the AMOEBA force field. However, the polarizabilities of phosphorus and double-bonded oxygen were reoptimized to better reproduce the QM molecular polarizability of DMP and TMP. We found that double-bonded oxygen polarizabilities for the DMP ion (1.724 Å3) needed to be much larger than those for TMP (0.892 Å3), whereas the phosphorus polarizabilities of the two molecules could be the same (1.788 Å3). Table lists the fitted and theoretical molecular polarizabilities. The Thole damping coefficients for all atom types of TMP and DMP were set to 0.390, as in previous work.

Table 3

Comparison of AMOEBA and QM Molecular Polarizabilities of DMP and TMP in Various Conformationsa

(A) DMP molecular polarizability tensor eigenvalues (Å3)
gg QM	gg AMOEBA	gt QM	gt AMOEBA	tt QM	tt AMOEBA
11.09	11.14	11.50	11.45	12.24	12.36
10.68	10.80	10.24	10.20	9.65	9.70
9.26	9.04	9.42	9.37	9.31	9.09

The QM values were calculated using the wB97xD/aug-cc-pVTZ method. With fitted atomic polarizability parameters (P and O), the AMOEBA molecular polarizabilities were calculated using the TINKER POLARIZE program.

Initial Bond Stretch and Angle Bend Parameters

The initial parameters for equilibrium bond lengths and bond angles were assigned as the average QM value in the stable conformations (gg, gt, and tt for DMP; C3, C1, and pan class="Chemical">Cs for TMP). The initial force constant parameters for bond stretches and angle bends were assigned to best match AMOEBA and QM vibrational frequencies. Except for the force constants of the O–P–Os and Os–P–Os angles, other force constant parameter were assumed to be the same between DMP and TMP. Theoretical vibrational frequencies were computed using the MP2/cc-pVTZ method. The TINKER VALENCE program was used to calculate the modeled molecular frequency. The frequencies were most strongly related to the force constants of the C–H and Os–C bonds and the H–C–H and Os–C–H angles. Thus, these four parameters were first optimized starting from typical values from prior AMOEBA force fields. These parameters were then fixed while the other force constants parameters were optimized. Figure shows the vibrational frequencies of TMP with final AMOEBA force field parameters and the corresponding QM vibrational frequencies as well as the values from infrared spectroscopy (IR) for gas-phase TMP.[40]

Figure 5

Comparison of TMP vibrational frequencies calculated using AMOEBA and QM as well as the values from infrared spectroscopy (IR) for gas-phase TMP. The QM values were calculated using the MP2/cc-pVTZ method. The AMOEBA vibrational frequencies were calculated with the C3 conformation using the TINKER VALENCE program.

van der Waals Parameters

Six vdW parameters were optimized by fitting to the DMP and TMP interaction curves with water, shown in Figure : the van der Waals radius and well depth of single-bonded oxygen (Os), double-bonded oxygen (O), and phosphorus (P). In our present model, DMP and TMP share the same van der Waals (vdW) parameters for phosphate atoms. The DMP–water or TMP–water dimer structures (Figure A,C) with one or two hydrogen bonds between water and phosphate oxygen atoms were used. The optimal vdW radius and well depth of phosphorus are 2.225 Å and 0.300 kcal/mol, respectively. The double-bonded oxygen is smaller than single-bonded oxygen, but the potential well depth is larger for both DMP and TMP. The oxygen atoms of TMP are slightly smaller than those of the DMP anion. With these parameters, the AMOEBA interaction curves are in good agreement with the QM results (Figure B,D).

Figure 6

Ab initio and AMOEBA interaction energy curves for DMP–water (A, B) and TMP–water (C, D). (A) For DMP, waters form two hydrogen bonds with atoms Os/Os of the tt conformer, O/Os of the gg conformer, and O/O of the gt conformer. (C) For TMP, waters form one or two hydrogen bonds with atoms Os(g)–Os(−g) of the Cs conformer, O of the C3 conformer, and Os(t)–O of the C1 conformer. The two hydrogen-bond distances in configurations with double H-bonds were made to be equal when sampling was performed around the equilibrium values. Single-point energies of monomers and dimers were calculated using MP2/cc-pVQZ methods, and the interaction energy is given by the difference between the dimer energy and the sum of the energy of the monomers (with basis set superposition error corrections). The force field parameters of water were taken from the 2003 model,[6] distributed as part of the TINKER 6 software package.

Ab initio and AMOEBA interaction energy curves for DMP–pan class="Chemical">water (A, B) and TMP–water (C, D). (A) For DMP, waters form two hydrogen bonds with atoms Os/Os of the tt conformer, O/Os of the gg conformer, and O/O of the gt conformer. (C) For TMP, waters form one or two hydrogen bonds with atoms Os(g)–Os(−g) of the Cs conformer, O of the C3 conformer, and Os(t)–O of the C1 conformer. The two hydrogen-bond distances in configurations with double H-bonds were made to be equal when sampling was performed around the equilibrium values. Single-point energies of monomers and dimers were calculated using MP2/cc-pVQZ methods, and the interaction energy is given by the difference between the dimer energy and the sum of the energy of the monomers (with basis set superposition error corrections). The force field parameters of water were taken from the 2003 model,[6] distributed as part of the TINKER 6 software package.

Torsional Parameters

The rotation around P–Os ester bonds defines the conformational energy surface of DMP and TMP (Figure ). The conformational energy in a classical force field is a sum of intermolecular (1–4 and beyond) vdW and electrostatic interaction energies plus the intrinsic torsional energy that is typically represented by a Fourier series. There are three types of torsion angles present in the DMP and TMP systems: O–P–Os–C, Os–P–Os–C, and P–Os–C–H.

The AMOEBA conformational energy surface was first constructed with the relevant torsional parameters set to zero. Then, these parameters were obtained by fitting the torsional energy function to the energy difference between the QM and AMOEBA conformational energy surface. For the P–Os–C–H torsion, the three C–H bonds have 3-fold symmetry, and only the n = 3 term is needed (kpan class="Chemical">P-Os-C-H,3 = 0.12 kcal/mol in both DMP and TMP). Because the O–P–Os–C and Os–P–Os–C torsion angles are coupled, their parameters were optimized together, and different parameters were used for DMP and TMP. As illustrated in Figure , the conformational energy surface was sampled as a 37 × 37 grid for DMP and a 9 × 9 × 9 grid for TMP. In this initial parametrization step, the QM conformational energy surface was calculated at MP2/cc-pVTZ without PCM solvation. The nonlinear optimization method lsqnonlin in Matlab[26] was used to optimize the torsion parameters. The objective function of the optimization was the sum of all of the squared differences between the torsional potential energy (O–P–Os–C + Os–P–Os–C calculated) and corresponding QM target values. Configurations with energies higher than 8.0 kcal/mol when compared to the most stable configuration were excluded during the fitting process.

Stretch-Torsion and Angle-Torsion Coupling

Using only the above traditional force field valence terms, it is difficult to reproduce the geometry and conformational energy accurately due to anomeric effects associated with the phosphate group[41,42,2] (for example, see Figure B), resulting in significant coupling among bonds, angles, and torsions. Figure shows the general geometry variation maps generated by pan class="Gene">MP2/cc-pVTZ with PCM for DMP and TMP. To capture such effects, new stretch-torsion and angle-torsion coupling terms have been introduced. The parameters for these new coupling terms were assigned by fitting to the corresponding bond or angle variations.

Figure 7

Geometry variation with torsion in DMP and TMP generated by the anomeric effect. The anomeric effect is explained in (B), where the blue arrow indicates the effective donation of the lone pair of Os1 to the antibonding orbital in the gt conformation. DMP geometries: (A) P–Os1 both length and (C) P–Os1–C1 and (D) O1–P–Os1 angles are shown as functions of torsion angles χ1 (Os2–P–Os1–C1) and χ2 (Os1–P–Os2–C2). The atom labels and optimized structures are the same as those in Figure A. (F) Coupling relations between Os–P–Os–C torsion and O1/O2–P–Os–C torsion in DMP. (E) TMP O–P–Os and P–Os–C angles that vary with O–P–Os–C torsion ψ3 have 180° and 120° periods. The TMP O–P–Os–C torsions ψ1 and ψ2 were both fixed to 45°. The atomic labels and optimized structures are the same as those in Figure B.

The stretch-torsion coupling potential energy for an A–B–C–D torsion χ is expressed as shown in eq where b1, b2, and b3 are the observed bond lengths A–B, B–C, and C–D, respectively, and b10, b20, and b30 are the corresponding equilibrium lengths, taken to be equal to the bond stretch energy equilibrium values. If the coupling coefficient k is positive, then a smaller bond length b is more favorable; otherwise, a larger bond length b is more favorable. Figure A shows how the P–Os1 bond length changes with χ1 and χ2, resulting in map with periodic variation. The P–Os1 bond is the first bond (b1) along the χ2 torsion and reaches a minimum length at χ2 = 180°. The P–Os1 bond is also the second bond (b2) along the χ1 torsion and reaches a maximum length at χ1 = 180° or 0° (especially when χ2 is near ±70°). On the basis of these observations, the terms m = 1 and n = 1 with a negative coupling coefficient and m = 2 and n = 2 with a positive coupling coefficient were chosen. The P–Os bond lengths of high-level QM-optimized structures (Tables B and S3B) were used to determine the optimal stretch-torsion parameters. Note that it is only necessary to couple the Os–P–Os–C torsion with its first and second bonds. With k = −4.25 kcal/mol/Å and k = 1.80 kcal/mol/Å, the AMOEBA-optimized structures of both DMP and TMP are in optimal agreement with the ab initio Os–P bond variation. In addition, the equilibrium bond length for b was shifted to a 0.010 Å shorter length than the initial value following addition of the stretch-torsion coupling term.

Table 4

DMP Geometry Comparison between QM and AMOEBA with Implicit Solvent Models (PCM or GK Solvation)

(A) structural root-mean-square deviation (RMSD)
conformation	gg	gt	tt	gg → gt	gt → tt
RMSD Å (heavy atoms)	0.030	0.039	0.019	0.036	0.106
RMSD Å (all atoms)	0.066	0.071	0.029	0.061	0.118

An angle-torsion coupling potential term was also introduced to capture the variations in the bond angle with the torsional angle in DMP and pan class="Chemical">TMPwhere a1 and a2 are bond angles A–B–C and B–C–D, respectively, contained in torsion A–B–C–D, and a10 and a20 are the equilibrium angles, which take the same value as in the angle bend potential term. By analyzing the variation of the three coupled phosphate angles O–P–O, Os–P–Os, and O–P–Os (Figures D and S1), we found periodicities analogous to the stretch-torsion case. The two O–P–Os1–C torsional angles are tightly coupled with the Os2–P–Os1–C torsion since they rotate about the same P–Os1 bond (Figure F). Thus, it is convenient to have all angle-torsion coupling terms associated with only the O–P–Os–C torsional angle. Also note that there are four O–P–Os–C torsions and just two Os–P–Os–C torsions in DMP. From Figure D,F, it can be seen that when O1–P–Os1–C1 equals 180° then Os2–P–Os1–C1 (χ1) is ∼65° (Figure F) and the O1–P–Os1 angle reaches a minimum value (Figure D). When O1–P–Os1–C1 equals 90°, Os2–P–Os1–C1 (χ1) is ∼160°, and the O1–P1–Os1 angle reaches a maximum. Thus, the period of O1–P–Os1 with respect to O1–P–Os–C is approximately 180°. On the basis of this observation, the terms m = 1 and n = 2 with a negative coupling coefficient were selected to model coupling between O–P–Os and O–P–Os–C.

As shown in Figure C, the P–Os1–C1 angle varies most with rotation of the χ1 torsion. The period of P–Os1–C1 with respect to Os2–P–Os1–C (or O–P–Os–C) is about 120°, exhibiting minima at χ1 = 60°, 180°, and 300°. Similarly, the variation of the P–Os–C angle with the O–P–Os–C torsion in TMP also showed a 120° period (Figure E). Thus, the P–Os–C angle and O–P–Os–C torsion coupling terms m = 2 and n = 3 were assigned with a negative coupling coefficient. High-level ab initio structures (Table C and S3B) were used to determine the two coefficients (kangtor,12 = −0.115 kcal/mol/degree and kangtor,23 = −0.011 kcal/mol/degree) for coupling the O–P–Os and P–Os–C angles with the O–P–Os–C torsion, respectively. The same parameters are used for both DMP and TMP. The equilibrium values for angles O–P–Os (a10), P–Os–C (a20), and O–P–O were optimized together with these two coupling parameters. The final equilibrium value for the O–P–O angle was 2.6° greater than the initial value obtained before the introduction of angle-torsion coupling.

Bonded Parameter Refinement

The five types of bonded potential parameters (bond stretch, angle bend, torsional angle, stretch-torsion, and angle-torsion coupling) were optimized separately and by different methods detailed in the previous sections. Initial bond and angle parameters were fit to the molecular frequency, torsional parameters were fit to the conformational energy surface, and coupling term parameters were fit to the structures. The final set of parameters was then fine-tuned together against various structural and energetic properties of both DMP and pan class="Chemical">TMP, including the structures of stable conformations and the activation energy and energy gap for conformational transitions. The details for this final refinement process, which involved an iterative parameter optimization protocol, are included in the Supporting Information.

With the final parameters in hand, the differences in bond lengths between AMOEBA- and QM-optimized structures are less than 0.004 Å, and the differences in bond angles were below 2.0°. Tables and S3 provide detailed structure comparisons between QM and AMOEBA for DMP and pan class="Chemical">TMP, respectively. The maximum root-mean-squared deviation (RMSD) between QM and AMOEBA structures of the three DMP conformations is 0.039 Å, and the maximum RMSD of the three TMP stable conformations is 0.078 Å (hydrogen atoms not included). It is noted that the difference between the two P–Os bond lengths in the DMP gt conformation and the P–Os bond variations between conformations, caused by the anomeric effect, are correctly captured by AMOEBA. Tables and S4 show the relative conformational energies for DMP and TMP using the final AMOEBA force field parameter sets. The root-mean-squared error (RMSE) between AMOEBA and QM, with or without implicit solvent, for the seven DMP configurations list in Table is 0.129 and 0.193 kcal/mol, respectively. The AMOEBA conformational energy surfaces of DMP and TMP were compared with the corresponding QM energy surfaces (Figures , S2, and S3). The QM and AMOEBA energy maps are in excellent agreement, showing the same local energy minimum positions and interconversion pathways between minima. Except for the very high-energy area on the potential map (i.e., the four corners in the 2D DMP conformational energy map in Figures A and the center of the 3D TMP map near 180°/180°/180°), the energy differences are well within 1.0 kcal/mol.

Table 5

Relative Conformational Energies of DMPa

	methods
conformation	QM without PCM	AMOEBA without GK	QM with PCM	AMOEBA with GK
gg	0	0	0	0
gt	0.8666	1.2775	1.5859	1.5494
tt	1.5546	1.5950	2.9875	2.6218
gg → gt 95	0.4613	0.2759	1.4144	1.4600
gg → gt 125	2.0109	1.9469	2.8929	2.9124
gg → gt 155	1.5619	1.8613	2.3712	2.2915
gt → tt	2.6268	2.7868	3.7015	3.6603

Energy is in units of kcal/mol. The QM energy was calculated at the MP2/cc-pVQZ PCM level with or without PCM implicit solvation, after the structures were optimized using the MP2/cc-pVQZ + PCM method with a sodium/water environment (Figure C). The AMOEBA energy was calculated with or without GK implicit solvation,[23] after the structures were optimized with GK solvation contribution. The transition state conformations gg → gt 95, 125, and 155 are explained in the Supporting Information, section II.

Energy is in units of kcal/mol. The QM energy was calculated at the MP2/cc-pVQZ pan class="Chemical">PCM level with or without PCM implicit solvation, after the structures were optimized using the MP2/cc-pVQZ + PCM method with a sodium/water environment (Figure C). The AMOEBA energy was calculated with or without GK implicit solvation,[23] after the structures were optimized with GK solvation contribution. The transition state conformations gg → gt 95, 125, and 155 are explained in the Supporting Information, section II.

Comparison of conformational energy surfaces, including implicit solvation, calculated by QM and AMOEBA for DMP and pan class="Chemical">TMP. (A) DMP χ1–χ2 2D potential energy map. (B) A slice with ψ3 = 45° from the TMP 3D conformational energy map.

TMP Parameter Validation via Liquid and Hydration Simulations

To validate the TMP AMOEBA force field parameters for liquid-phase applications, the hydration free energy, the density of liquid pan class="Chemical">TMP at a few temperatures, and the heat of vaporization at the boiling point (470.35 K) were calculated from molecular dynamics simulations and compared to the experimental values (in Table 4 of Carl Yaws handbook).[43] The results show that the TMP force field reasonably predicts all of these properties (Figure E,F). The TMP solvation free energy at 298 K calculated using the BAR method was 7.4 ± 0.4 kcal/mol, compared with the experimental value of 8.7 kcal/mol.[44,45] The underestimation of HFE is likely related to the overestimated repulsion in the TMP–water dimer interaction at short distance (Figure C,D). The double-bonded oxygen–water hydrogen bond (O–Hw; the middle structure of Figure C) is the dominant interaction in water (see below for further discussion). The repulsion between water and TMP is overestimated by 1 kcal/mol when the O–Hw distance is at 1.6 Å in comparison to QM and becomes worse when the distance is shorter. The shape of the energy profile cannot be improved by reducing the vdW radius while maintaining the minimum energy distance. We believe that this discrepancy can be addressed only by including the charge penetration effect,[46] which requires significant modification of the whole potential energy function. The errors in density at 227, 298.15, 373, and 480.35 K are 2.1, 1.2, 0.5, and −1.4%, respectively (Table S5, experimental values from Carl Yaws handbook[43]). From 4 ns liquid-state MD simulations at the boiling point of 470.35 K, the average liquid potential energy was found to be U = −31.56 kcal/mol per TMP. The average potential energy for a TMP molecule in the vapor phase was U = −21.60 kcal/mol. Thus, the heat of vaporization at 470.35 K is estimated to be 10.89 kcal/mol, according to ΔH = ΔU + pΔV = U – U + RT. The experimental values from Tables 9 and 10 in Carl Yaws handbook[43] are 10.27 and 11.48 kcal/mol, respectively.

Figure 9

Results from molecular dynamics simulations of liquid TMP and TMP in water. (A) Conformation distribution comparison of liquid TMP at 298 K and TMP in water solution at 298 K. The three letters indicate the three O–P–Os–C torsion angles. The letter c (between −30° and 30°) stands for cis; g (between 30° and 90°) and −g (between −30° and −90°) stand for gauche; a (between 90° and 150°) and −a (between −90° and −150°) stand for anticlinal; and t (between 150° and 210°) stands for trans. (B) Comparison of the induced dipole distributions of liquid TMP and TMP in a dilute water environment. The average induced dipole moments of TMP in neat liquid and in dilute water are 0.631 and 1.550 D, respectively. (C) Radial distribution function and coordination number of TMP double-bonded oxygen–water oxygen (O–Ow) and TMP carbon–water oxygen (C–Ow). (D) Radial distribution function and coordination number of phosphorus–phosphorus (P–P), carbon–carbon (C–C), and double-bonded oxygen–carbon (O–C). (E) Comparison of calculated and experimental solvation free energies at 293 K and heat of vaporization of TMP at its boiling point. (F) Predicted density values of TMP compared with the experimental values. Densities were computed by molecular dynamic simulations as described in the text.

Additional structural insight into TMP in a dilute pan class="Chemical">water solution and liquid TMP was obtained from MD simulations. First, the TMP molecules exist mostly in C1, Cs, or C1 → Cs transition conformations in water (corresponding to ggt, g-gt, and ctg in Figure A), which are strongly polar (the dipole moments in the gas phase are ∼3.6 D). However, TMP molecules in neat liquid show a diverse conformational distribution, and the dominate conformations are C3 -like, which are weakly polar (ccg, cgg, cg-g, ccc, and ggg in Figure A; the dipole moments in the gas phase are ∼1.0 D). As the temperature increased, the population of polar conformations in liquid became larger (Figure S4). Second, the induced dipoles of TMP in liquid and in a water environment are also quite different (see the induce dipole distribution in Figure B). The average induced dipoles of TMP in liquid and in dilute water were 0.631 and 1.550 D, respectively. Third, from the radial distribution function curves (Figures C, black line, 9D, blue line, and S6), only the double-bonded oxygen, not single-bonded oxygens, of TMP form hydrogen bonds with water or the methyl group of other TMP molecules both in water and in neat liquid. In neat liquid, the related angles (Figure S5) also support the existence of such hydrogen-bonding interactions. Finally, we obtained additional structural information from the simulations. In TMP liquid, the phosphorus–phosphorus radial distribution function (Figure D, black line) shows that the first shell around TMP, from 4.8 to 8.0 Å, contains about 11 TMP molecules. In TMP liquid, each CH3– is surrounded by about seven other CH3– within 3.6–5.0 Å (Figure D, red lines). In TMP–water solution, about four water molecules surround each CH3– within 3.6–5.0 Å (Figure C, red lines).

DMP–Sodium and DMP–Magnesium Ion Interactons

In a recent report, a liquid simulation of dimethyl phosphoric acid (pan class="Chemical">DMPH) was used to develop and test the OPLS force field of DMP.[47] A revised AMBER parameter was also reported through fitting the pKa and solvation free energy of DMPH.[48] However, as Table S2 shows, the geometry of DMPH and, of course, its charge distribution are very different from those of DMP itself. In this work, the AMOEBA parameter set for DMP was further tested on the interaction with metal ions. A sodium or magnesium ion was placed at the bisector of the O–P–O angle of the DMP ion, as illustrated in Figure S7. These dimers were optimized by both QM calculations at the MP2/aug-cc-pVTZ level and with the proposed AMOEBA force field for DMP. Both the structures and interaction energies of DMP–Na+ or DMP–Mg2+ show very good agreement between QM and AMOEBA (Table ). Note that AMOEBA uses the same Thole damping coefficient (0.39) for all organic molecules and monovalent ions. For divalent ions, such as Ca2+, Mg2+, and Zn2+, we discovered previously that different damping coefficients are necessary to capture the difference in the electronic structure of these ions when they are involved in polarization interactions.[24,25,49] For DMP–Mg2+ polarization, the optimal damping coefficient was found to be 0.15, which is slightly larger than the 0.09 value that we previously reported for water–Mg2+.[24]

Table 6

Comparison of the DMP–Metal Ion Interaction Geometry and Energy Calculated by QM and AMOEBAa

system (methods)	M–P (Å)	M–O1 (Å)	M–O2 (Å)	interaction energy (kcal/mol)
DMP/Na+ (QM)	2.726	2.277	2.290	–138.2
DMP/Na+ (AMOEBA)	2.828	2.058	2.593	–138.3
DMP/Mg2+ (QM)	2.464	1.925	1.928	–391.4
DMP/Mg2+ (AMOEBA)	2.552	1.878	1.886	–391.3

Sodium and magnesium ions were included in this study. The force field parameters of these two metal ions were taken from TINKER 7.1.[22] The Mg2+–DMP/TMP Thole damping coefficient was set to 0.15, whereas 0.0952 was used for Mg2+–water.

Discussion and Conlusions

The AMOEBA polarizable force field for DMP and TMP has been elaborated. The parameters were largely derived from high-level QM calculations. The structure, energy, electrostatic fields, vibration frequencies, and molecular polarizability of DMP and TMP generated by the final AMOEBA models show good agreement with QM results. As a further validation, we randomly selected monomers and dimers from the liquid MD simulations of TMP and calculated intramolecular energy and interaction energy using both QM and AMOEBA force fields. The AMOEBA-calculated interaction energy has an excellent correlation coefficient with the QM interaction energy (Figure S8B).

We found two unique aspects of DMP that require special consideration during force field development. First, the structure of pan class="Chemical">DMP is significantly different between the gas and solution phases as well as with or without a metal ion bound to the phosphate. As we are mainly targeting applications for DMP and nucleic acids in solution environments where counterions are typically present, we have chosen the corresponding ab initio calculations for use in force field parametrization. Second, we noticed that the bond lengths and angles are strongly coupled with the torsion angles in DMP and TMP due to a generalized anomeric effect, which is much stronger than what is typically observed for first-row organic molecules. Such phenomena could be important for the conformational flexibility of DNA and RNA. To address this effect, new stretch-torsion and angle-torsion coupling terms have been introduced into the force field. The resulting model is able to reasonably capture the bond and angle variations with respect to DMP/TMP conformational changes. The difference in bond lengths between AMOEBA- and QM-optimized structures is about 0.005 Å or less, differences in bond angles are mostly less than 0.8°, and the errors in relative conformational energy for important conformations are within 0.6 kcal/mol. Furthermore, the difference between the two P–Os bond lengths in the DMP gt conformation and the P–Os bond variations between conformations, caused by the anomeric effect, is correctly captured by AMOEBA. With the coupling terms, we calculated the correlation coefficients between the coupling energy and valence/torsion energy using the DMP structures that cover the whole conformational space (see Figure S9). It seems that the relationship between the valence energy and cross-terms is rather complicated and highly nonlinear.

On the basis of our previous experience, a solely QM-derived force field does not necessarily perform well in condensed-phase simulations. As a validation of this, liquid and dilute water simulations of pan class="Chemical">TMP were carried out. We found that the dominant conformations and induced dipole moments are quite different between these two systems (Figure A,B). The latter suggests that it is difficult for a nonpolarizable force field that implicitly includes polarization effect to transfer between neat liquid and hydration environments. Nevertheless, polarizable AMEOBA is able to capture thermodynamic properties in both environments. The errors in predicted TMP liquid density at temperatures ranging from the melting point to the boiling point were all below 2.1%. The predicted heat of vaporization and solvation free energy are also in reasonable agreement with experimental values. The latter is somewhat underestimated by the force field, which we believe is a result of the lack of short-range charge penetration in the current AMOEBA function. In addition, the binding geometry and energy of DMP with Na+ or Mg2+ ions compare favorably with the corresponding QM-calculated results. These comparisons provide further validation of this force field, as these ion parameters have been developed previously to work as components of a wide range of molecular systems.

Overall, the AMOEBA polarizable atomic-multipole-based force field for DMP and pan class="Chemical">TMP provides a reliably accurate description of molecular structures, conformational properties, and interaction energies with water and metal ions in the gas phase, as well as liquid-phase thermoproperties for TMP. It should be noted that the parametrizations and validations in this study, particularly for the DMP anion, heavily relied on QM calculations. Thus, future examination of condensed-phase systems of DMP salts and acids will be important for further improvement and validation. Nonetheless, we expect that this force field will be useful for computational studies of DMP and TMP in various applications and in different environments and that it lays the foundation for future model development for other phosphate-containing molecular species such as lipids and nucleic acids.