Charge density distribution and theoretical analysis of low and high energy phosphate esters

Abstract

Charge density analysis has been performed for crystals of ammonium dihydrogen phosphoenolpyruvate (I) and diammonium hydrogen phosphoglycolate (II). The former compound is a high-energy phosphate ester that stores and provides energy for various biochemical processes, whereas the latter is a low energy one. The main difference between both molecules is mainly manifested by topological parameters of the P–O and C–O bonds within the ester C–O–P bridge. The charge density at the bond critical point in I is larger for C–O and smaller for P–O than in II. Theoretical DFT and NBO calculations show that the ester O atom has an approximate sp² hybridization with two lone pairs: the sp^x (x = 1.2–1.6) and the p one. The latter, which is perpendicular to the C–O–P plane, interacts with the π*(C–C) orbital of the enolpyruvate fragment in I. Extended DFT calculations, performed for 15 model organic phosphates showed that the interaction of the p lone pair (called tlp hereinafter) with the organic part determines predominantly the geometry of the C–O–P fragment as well as the energy of hydrolysis. The interaction of tlp with a neighboring π* or a low-lying σ* orbital from the organic part creates an additional bonding effect along C–O, while the interaction of tlp with the σ* orbitals of the terminal P–O bonds strengthens the ester P–O bond. These two interactions are negatively correlated and the former is determinant. The P–O and C–O distances are additionally modified by interactions of the ester σ(P–O) and the sp² pair with the organic part. On the other hand the hydrolysis energy is proportional to the interaction of tlp with the organic part. Regression analysis of the DFT results shows that the hydrolysis energy may be interpreted as being released mainly from the sp^x pair and the σ(P–O) orbital. It has been also shown that this energy may be in principle deduced from the topological parameters retrieved from the experimental charge density distribution.

---

Introduction

Organic phosphates are important metabolic intermediates. A particular role is played by the so called “high energy phosphates” i.e. phosphate esters which, upon hydrolysis, fuel numerous, otherwise energetically unfavourable, reactions and processes of living cells.¹ These compounds are characterised by relatively high absolute values of standard Gibbs free energy or free enthalpy of hydrolysis ΔG^o_hydr.

One of the most widely recognized examples of energy storing phosphates is adenosine triphosphate or ATP. Its hydrolysis to adenosine mono- or diphosphate (AMP or ADP) provides the energy for numerous biological processes e.g. muscle tension, cellular motion, proton gradient generation and the synthesis of various biological molecules and components. Regeneration of ATP requires, however, energy that must be provided by hydrolysis of another high-energy molecule, e.g. phosphoenolpyruvic acid (H₂O₃POC(CH₂)COOH, abbreviated as H₃PEP). Table 1 shows some of the metabolic organic phosphates with their corresponding standard ΔG^o_hydr.^1,2 Conventionally, ATP delimits the high-energy phosphates from the low-energy ones.

Table 1 Comparison of various organic phosphate esters and their corresponding standard free enthalpy of hydrolysis ΔG^o_hydr

Compound	ΔG^ohydr
	kJ mol^−1	kcal mol^−1	Ref.
Glycerol-3-phosphate	−9.2	−2.2	2
Methyl monophosphate	−9.2	−2.2	2
Glucose-3-phosphate	−9.2	−2.2	1
Glucose-6-phosphate	−13.8	−3.3	1 and 2
Glucose-1-phosphate	−20.9	−5.0	1 and 2
Adenosine triphosphate	−30.5	−7.3	1
Acetyl phosphate	−43.1	−10.3	2
1,3-Bisphosphoglycerate	−49.4	−11.8	1 and 2
Carbamoyl phosphate	−49.4	−11.8	2
Phosphoenolpyruvate	−61.9	−14.8	1

---

Because of the role the phosphate esters play in the metabolism, the question about the reasons behind the variation of their ΔG^o_hydr has attracted the interest of many researchers (e.g. ref. 3–10 and many others).

For instance, one of the widely used textbook explanations in regard to ATP in its physiological form ATP⁴⁻ emphasizes the role of resonance stabilization of the products, electrostatic repulsion between the neighboring phosphate groups and stabilization due to solvation effects;¹ nevertheless this explanation is rather of qualitative character and cannot offer quantitative predictions of the ΔG^o_hydr for other phosphate compounds. One of the proposed explanations, known as the ‘competition for the lone pair’, stresses the role of the competition for the electron density of lone pairs of the bridging ester oxygen atom between the phosphate moiety and the organic rest.¹¹ Following this rationalization, greater ability to attract the electron density by the organic substituent would result in diminished strength of interaction between the bridging O atom and the phosphate moiety, thus leading to the creation of less stable, more exothermic molecules.

Several reports have focused on the role of negative hyperconjugation, referred to as the generalized anomeric effect, in phosphate reactivity.^12–14 It has been shown by Kirby et al.^15–17 on grounds of crystallographic and kinetic studies that a structure–reactivity correlation can be drawn in reference to the lengths of scissile ester P–O bonds. First observations made for another class of compounds, heterocyclic acetal esters, indicated that a relationship between the effective electronegativity of the leaving group, as measured by the pK_a of the appropriate alcohol and the length of the scissile C–O bonds exists. It has been proposed that both acetal oxygen atoms (endocyclic and exocyclic one) take part in n(O)–σ*(C–O) donor–acceptor interactions, weakening thus the scissile bonds and reinforcing the neighboring ones.

Kirby et al. hypothesized that the mechanism behind the observed correlation of structure and reactivity would rest in energy changes of the bridging σ*(C–O) orbital. It was argued that for labile, high-energy compounds lowering of σ*(C–O) energy caused by the increased electronegativity of the substituents would promote additional donation from the oxygen lone pairs. This would also result in the elongation and weakening of the scissile bridging C–O bond on one hand, and in shortening and strengthening of the endocyclic acetal one on the other. An interplay scheme was thus proposed where the two bonds were either reinforced or weakened in an opposing manner by donor–acceptor interactions involving oxygen lone pairs. These observations were then generalized to phosphate esters, where the scissile ester P–O bonds play in a sense a role similar to the exocyclic C–O bonds.

The influence of oxygen electron lone pairs on P–O bonds has been demonstrated in a number of studies. For instance, in early Gorenstein and coworkers' papers^18–20 on phosphate and phosphorane compounds the ester P–O bonds were shown to change depending on the orientation of the lone pairs of the ester O atoms.

Extensive evidence of the influence of the anomeric or hyperconjugation effects on the P–O bond lengths and ΔG^o_hydr was given by Ruben et al.² Based on the NBO analysis^21,22 of a variety of phosphate compounds, these authors reached the conclusion that there is a positive correlation between the extent of the generalized n(O) → σ*(P–O) anomeric effect and the values of ΔG^o_hydr. This was evidenced by the analysis of the σ*(P–O) orbital populations and the energies of orbital interactions based on second order perturbation correction (E⁽²⁾_ij). According to these authors, similar correlations also hold with regard to the P–O bond lengths, the σ*(P–O) orbital energies and the polarity of ester P–O bonds. In other words, their results indicated that the more polar the studied σ*(P–O)_ester orbitals and the lower their energies are, the higher their populations will be. This corresponds to longer ester P–O bonds and greater ΔG^o_hydr values for the studied set of model phosphates. In this way they obtained results which were in line with the mechanism proposed by Kirby et al. for scissile C–O and P–O bonds.

Studies on the role of coordination and solvation of the phosphate group^23–25 demonstrated the role of stabilisation of phosphate terminal oxygen lone pairs in various media.

These studies have, nevertheless, focused mainly on the local interactions that could largely be attributed to the phosphate group. Factors concerning the organic parts of studied molecules were accounted for in a very limited degree only. Our aim is therefore to investigate the role of these fragments. Our interest has been sparked by the fact that the observed high values of ΔG^o_hydr are often associated with molecules that have a double bond in the nearest vicinity of the ester P–O–C bridge.¹¹ This is the case of phosphoenolpyruvate, acetyl phosphate and carbamoyl phosphate, as it may be seen in Table 1.

The scope of this paper includes a study of experimental charge density and theoretical calculations of two molecules of some biological significance: dihydrogen phosphoenolpyruvate, H₂PEP⁻, in a crystal of ammonium dihydrogen phosphoenolpyruvate and hydrogen phosphoglycolate, HPG²⁻ (HO₃POCH₂COO²⁻), in a crystal of diammonium hydrogen phosphoglycolate. Structurally these two compounds differ from each other by the presence of an additional methylidene substituent at the ester C atom in the former. Phosphoenolpyruvic acid has the highest value of ΔG^o_hydr of the organic phosphates involved in glycolysis.¹ On the other hand phosphoglycolic acid is a low-energy phosphate that takes part in the plant photorespiration process;²⁶ it has also been found in the tissues of humans.²⁷ This choice was dictated by comparable ionization states of the phosphate group and the presence of the same countercation i.e. NH₄⁺. In order to get a wider insight into the interactions between the phosphate group and the organic rest we have performed DFT and NBO calculations for 15 model molecules. The results have allowed to rationalize the experimental results and to propose a connection between the molecular structure and such features as the geometry of the ester C–O–P bridge or the hydrolysis energy.

Experimental part

The crystals of ammonium dihydrogen phosphoenolpyruvate (I) and diammonium hydrogen phosphoglycolate (II) had been obtained as described previously.^28,29

Data collection and multipole refinement

Both crystals were subjected to preliminary check of quality. The diffraction data were collected on Xcalibur diffractometers equipped with a Mo anode, a graphite monochromator, an Oxford Cryosystem cooling device and a Ruby (I) or a Sapphire (II) CCD area detector. The data were corrected for Lorentz, polarization and absorption factors, and then scaled and merged with SORTAV³⁰ with the following results: I: R_int was 0.0233 for all 116 955 reflections, 0.0159 for 32 128 reflections below sin(θ)/λ = 0.7 Å⁻¹, 0.0454 for 61 909 reflections in the range 0.7 ≤ sin(θ)/λ ≤ 1.15 Å⁻¹; II: R_int was 0.0243 for all 134 992 reflections, 0.0164 for 27 559 reflections below sin(θ)/λ = 0.7 Å⁻¹, 0.0386 for 81 535 reflections in the range 0.7 ≤ sin(θ)/λ ≤ 1.15 Å⁻¹. For both crystals data up to sin(θ)/λ = 1.15 Å⁻¹ were used for the multipole refinement. Initial atomic data were taken from²⁸ and,²⁹ the formalism adopted was that of Hansen and Coppens³¹ and the refinements were performed with XD2006.³²

The approach adopted for multipole refinements was as follows. The C–H and O–H bond lengths were reset after each cycle to values obtained from the DFT cluster calculations; the N–H distances were taken from Allen et al.³³

The imposed constraints were: 1.09 Å for C–H bonds in I, 1.103 Å in II, 1.03 Å for N–H bonds; 1.03 Å for the phosphate hydroxyl O–H in I and 1.06 Å in II, as well as 1.05 Å for the carboxylic O–H in II. For the H atoms bonded to the main molecular skeleton of the anions, the harmonic anisotropic displacement parameters were estimated by the SHADE2.1 server³⁴ and kept fixed during the refinement. The ammonium H atoms were refined with isotropic atomic displacement parameters. Multipolar functions of phosphorus, oxygen, nitrogen and carbon atoms were expanded up to l = 3 and hydrogen atoms up to l = 2. However, for C-bonded H atoms only the quadrupole terms parallel to the C–H bond were refined, whereas for the O- and N-bonded H atoms all quadrupolar terms were refined in order to reproduce the deformations caused by strong hydrogen bonds. The oxygen atoms in both compounds were divided into classes based on their vicinity, with separate expansion/contraction parameters κ and κ′; see ESI† for the details. For the terminal phosphate O atoms, only the spherical κ was refined, while the multipolar κ′ values were fixed at 0.75 to prevent uncontrolled expansion of the radial functions. In the case of I an individual electric charge constraint (Q = +1) was applied to the NH₄⁺ cation The radial deformation functions R_l(r) were in the form of single Slater functions [ζⁿ⁺³/(n + 2)!]rⁿ exp(−ζr). For phosphorus n was 6, 7, 8 for dipoles, quadrupoles and octupoles, respectively, and the initial ζ value was 5.19 bohr⁻¹.

Correction for isotropic extinction (mosaic distribution dominated, Gaussian distribution) was included for both compounds, following the scheme given by Becker & Coppens.³⁵ The Hirshfeld³⁶ rigid bond test results indicated good quality of the model, with values for bonds between non-H atoms not exceeding 5 × 10⁻⁴ Ų. The crystal data and refinement details for both crystals are presented in Table 2.

Table 2 Data collection and refinement details

Crystal data	I	II
Chemical formula	C3H8NO6P	C2H11N2O6P
Mr	185.07	190.10
Crystal system, space group	Orthorhombic, Pbca	Monoclinic, P21/c
Temperature (K)	85	85
a, b, c (Å)	14.6460(7), 13.0224(9), 7.4486(6)	13.189(4), 8.079(3), 13.804(4)
β (°)		91.27(3)
V (Å^3)	1420.64(17)	1470.5(8)
Z	8	8
Radiation type	Mo Kα	Mo Kα
μ (mm^−1)	0.37	0.37
Crystal size (mm)	0.54 × 0.44 × 0.18	0.37 × 0.32 × 0.28
Data collection
Tmin, Tmax	0.845, 0.941	0.907, 0.928
No. of measured, independent and observed [>3σ(I)] reflections	116 955, 15 165, 13 369	134 992, 29 301, 24 953
Rint	0.023	0.024
Refinement
R[F^2 > 3σ(F^2)], wR(F^2), S	0.014, 0.014, 1.59	0.018, 0.017, 1.47
No. of reflections	8332	17 039
No. of parameters	356	760
No. of restraints	0	0
Δρmax, Δρmin (e Å^−3)	0.23, −0.21	0.29, −0.28

---

Theoretical calculations

The DFT calculations based on Slater-type functions were performed with the ADF suite of programs.³⁷ The functional PW91 ³⁸ and the basis – TZ2P, composed of double ζ functions for core electrons, triple ζ for valence ones and two polarization functions, with unfrozen cores, were used throughout. For the calculations the following systems were taken: I: (a) free dihydrogen phosphoenolpyruvate anion in vacuo; (b) a cluster of 5 dihydrogen phosphoenolpyruvate anions and 7 ammonium cations, so that the central anion might be provided its full environment. The cluster was placed in a cavity of dielectric medium simulating acetic acid to mimic the crystal environment (the COSMO model³⁹); II: (c) and (d) both symmetry independent (see p. 3.1) hydrogen phosphoglycolate anions in vacuo; (e) and (f) clusters of 3 hydrogen phosphoglycolate anions and 6 ammonium cations in either case, devised and set as in (b). Atomic coordinates were taken from the crystal geometry, but before proceeding to the final calculations the O–H and C–H bond lengths were optimized. The N–H distances were set at 1.03 Å, as given by Allen et al. The natural bond orbitals (NBO) analysis was performed with NBO 6.0.⁴⁰ The topological analysis of the theoretically calculated electron density was performed with Xaim.⁴¹

In order to analyze the relation between the energy of hydrolysis and the molecular structure of organic phosphates DFT calculations for a series of model molecules have been performed. The considered reaction was:

R–OPO₃H⁻ + H₂O → R–OH + H₂PO₄⁻

and molecules were as follows: dihydrogen phosphoenolpyruvate (PEPH₂⁻ or 1), hydrogen acetyl phosphate (2), hydrogen glucose-6-phosphate (3), hydrogen but-1-en-2-yl phosphate (4), hydrogen but-1-en-3-yl phosphate (5), hydrogen thioacetyl phosphate (6), hydrogen 1-glyceryl phosphate (7), hydrogen phosphoglycolate (8), hydrogen methyl phosphate (9), hydrogen fluoromethyl phosphate (10), hydrogen chloromethyl phosphate (11), hydrogen bromomethyl phosphate (12), hydrogen trifluoromethyl phosphate (13), hydrogen trichloromethyl phosphate (14) and hydrogen tribromomethyl phosphate (15). The formulae of these compounds are recapitulated in Scheme 1. The compounds 1, 2, 3, 7 and 8 are important biochemical intermediates, 1, 2, 4 and 6 have a double bond adjacent to the phosphate group and the molecules with halogen substituents (10–15) have been chosen to study the possible influence of increased positive charge of the ester C atom. The pair 4 and 5 represents the isomers differing in location of the double bond. As previously, the functional used was PW91 and the basis – TZ2P, with the exception of Br, where the basis was augmented (ATZ2P). Detailed tables containing the results and the quantities on which the discussion will be based may be found in the ESI.† We performed also parallel calculations with augmented bases AUG-ATZ2P; their results are briefly commented in the ESI.†

Scheme 1 Formulae of the model molecules used for the DFT calculations.

Results and discussion

Structures

Both structures have been described previously.^28,29 For the reader's convenience the dihydrogen phosphoenolpyruvate anion, together with its surrounding, is shown in Fig. 1, and both crystallographically independent hydrogen phosphoglycolate anions – in Fig. 2.

Fig. 1 Dihydrogen phosphoenolpyruvate (PEPH₂) molecule together with its environment. The labels are shown only for symmetry independent atoms, the dashed lines represent hydrogen bonds and only fragments of other H₂PEP⁻ anions linked by the hydrogen bonds are shown.

Fig. 2 Two symmetry independent hydrogen phosphoglycolate, (HPG²⁻) anions from II together with their surroundings. The labels are shown only for symmetry independent atoms, the dashed lines represent hydrogen bonds and only fragments of other HPG²⁻ anions linked by the hydrogen bonds are shown.

In II there are two symmetry independent hydrogen phosphoglycolate anions somewhat differing in conformation. For instance the torsion angle O2–P–O1–H1 in molecule A (i.e. the angle O2A–P1A–O1A–H1A) is −147°, while for molecule B – −175°. The bonds and angles in the ester parts of the molecules are listed in Table 3.

Table 3 Distances (Å) and angles (°) within the ester parts of phosphoenolpyruvate (I) and phosphoglycolate (II) anions

I
P–O1	1.5639(2)	O1–P–O2	107.328(9)
P–O2	1.4918(1)	O1–P–O3	112.724(9)
P–O3	1.4995(1)	O1–P–O4	105.758(10)
P–O4	1.6268(2)	O2–P–O3	118.960(9)
O4–C2	1.3721(2)	O2–P–O4	103.054(9)
		O3–P–O4	107.845(9)
		P–O4–C2	125.176(14)

---

II
P1A–O1A	1.5647(2)	O1A–P1A–O2A	99.932(10)
P1A–O2A	1.6037(2)	O1A–P1A–O3A	111.479(10)
P1A–O3A	1.5093(2)	O1A–P1A–O4A	111.783(11)
P1A–O4A	1.5069(2)	O2A–P1A–O3A	109.915(11)
P2B–O1B	1.5566(2)	O2A–P1A–O4A	108.794(11)
P2B–O2B	1.6016(2)	O3A–P1A–O4A	113.993(11)
P2B–O3B	1.5181(2)	P1A–O2A–C2A	118.147(15)
P2B–O4B	1.5025(2)	O1B–P2B–O2B	98.904(10)
O2A–C2A	1.4253(3)	O1B–P2B–O3B	112.200(11)
O2B–C2B	1.4212(3)	O1B–P2B–O4B	113.186(11)
		O2B–P2B–O3B	109.044(11)
		O2B–P2B–O4B	110.609(11)
		O3B–P2B–O4B	112.096(11)
		P2B–O2B–C2B	118.642(15)

---

Charges

The experimental topological atomic charges are fairly consistent with their theoretical analogues (see Table 4 for the COPO₃ fragments and ESI† for full data). The atomic charges of both compounds display no significant changes from other related compounds studied so far.^42–46 The influence of the organic moiety on the phosphate part is fairly small, as the atomic charges of phosphorus and oxygen atoms are similar in both molecules (H₂PEP⁻ and HPG²⁻). The P atoms have positive Bader charges of approx. 3.7e and the Hirshfeld ones of 0.3–0.4e and as such are consistent with previously studied phosphorus compounds, where the atomic charges of P atoms locate in the range of 2.92–3.78e (Bader) and 0.19–0.49 (Hirshfeld). Also the charges of O atoms are very stable, and almost independent of their chemical surroundings. In almost all compounds studied thus far the Bader charges of terminal O atoms are in the range of −1.2 to −1.8e, a rough estimate of the average being −1.5e, and the Hirshfeld charges locate within the range of −0.2 to −0.4e. However the experimental charges of P atoms in both anions of HPG²⁻ are somewhat smaller than those predicted by theoretical calculations. This is most likely an artefact of the multipole refinement, as a previous study of H₂PG⁻ (ref. 45) shows experimental and theoretical charges closer to current theoretical HPG²⁻ ones.

Table 4 Hirshfeld (Q_H) and Bader (Q_B) atomic charges for the COPO₃ fragments. Subscripts isol and cl indicate the charges obtained for isolated molecules and clusters, respectively

Atom	Multipole refinements		DFT calculations
	QH	QB	QB.isol	QH.isol	QB.cl	QH.cl
I (NH4H2PEP)
P	0.41	3.68	3.65	0.33	3.69	0.46
O1	−0.27	−1.48	−1.33	−0.25	−1.44	−0.23
O2	−0.31	−1.49	−1.49	−0.47	−1.49	−0.30
O3	−0.34	−1.51	−1.51	−0.48	−1.47	−0.28
O4	−0.24	−1.39	−1.23	−0.16	−1.30	−0.15
C2	0.06	0.46	0.52	0.05	0.50	0.06
II [(NH4)2HPG]
P1A	0.13	3.42	3.65	0.30	3.68	0.42
P2B	0.11	3.38	3.68	0.30	3.71	0.42
O1A	−0.36	−1.49	−1.33	−0.27	−1.43	−0.27
O2A	−0.33	−1.46	−1.30	−0.20	−1.30	−0.16
O3A	−0.39	−1.45	−1.55	−0.54	−1.49	−0.35
O4A	−0.38	−1.46	−1.54	−0.52	−1.49	−0.34
O1B	−0.37	−1.53	−1.33	−0.26	−1.45	−0.27
O2B	−0.30	−1.45	−1.29	−0.20	−1.30	−0.16
O3B	−0.40	−1.42	−1.55	−0.54	−1.50	−0.37
O4B	−0.39	−1.47	−1.55	−0.53	−1.51	−0.34
C2A	−0.10	0.15	0.47	−0.04	0.47	−0.01
C2B	−0.11	0.12	0.45	−0.04	0.43	−0.01

---

Also the ester O atomic charges of the organic moieties show no significant difference between the studied molecules, although in this study the H₂PEP⁻ ester O atom is more positive than its counterparts in HPG²⁻. It is perhaps also worth to notice that the experimental charges of the phosphate ester O atoms are systematically lower than those predicted from the theoretical calculations; such a difference was observed previously in KH₂PG⁴⁵ and sodium bis(p-nitrophenyl)phosphate.⁴⁶ Both Hirshfeld and Bader charges show that the ester C atom in H₂PEP⁻ experiences more electron withdrawal and is thus more positive than that in HPG²⁻, what corresponds with shorter C–O_e (O_e is the ester O atom in the C–O–P fragment) bonds in the former molecule. The methylene C atom is weakly negative, whereas the carboxylic C atoms, as in other carboxylate compounds,⁴⁷ are strongly positive.

Calculations performed in vacuo and in the cluster environment convey an additional opportunity to look into the shifts of the electron density induced by intermolecular interactions, allowing thus a crude evaluation of the effects of coordination, hydrogen bonding etc. to be made. In particular, comparison of the results for the isolated molecules and the clusters shows that inclusion of the hydrogen bonds brings about a slight depletion of negative charge of terminal O atoms.

Topology

The main topological parameters, namely the charge densities at the bond critical points (bcps), ρ_c, and their Laplacians, ∇²ρ_c, for selected bonds in I are collected in Table 5 and those in II – in Table 6 (full tables may be found in ESI†).

Table 5 Values of ρ_c (e Å⁻³) and ∇²ρ_c (e Å⁻⁵) for the for the –COPO₃ and –COO fragments in ammonium dihydrogen phosphoenolpyruvate (I)

Bond	Experiment		Theoretical, cluster		Theoretical, free anion
	ρc	∇^2ρc	ρc	∇^2ρc	ρc	∇^2ρc
P–O1	1.352(17)	16.03(6)	1.375	15.81	1.333	16.65
P–O2	1.686(21)	20.12(8)	1.622	27.55	1.638	27.55
P–O3	1.609(21)	19.90(8)	1.591	25.69	1.604	25.81
P–O4	1.138(16)	10.70(6)	1.179	11.69	1.145	12.47
O4–C2	1.940(13)	−15.29(6)	1.898	−12.79	1.958	−15.41
C1–O5	2.359(14)	−26.82(8)	2.310	−17.97	2.243	−14.19
C1–O6	2.882(19)	−30.30(11)	2.763	−12.83	2.800	−15.28

---

Table 6 Values of ρ_c (e Å⁻³) and ∇²ρ_c (e Å⁻⁵) for the –COPO₃ and –COO fragments in diammonium hydrogen phosphoglycolate (II)

Bond	Experiment		Theoretical, clusters		Theoretical, free anions
	ρc	∇^2ρc	ρc	∇^2ρc	ρc	∇^2ρc
P1A–O1A	1.38(3)	14.99(8)	1.40	16.41	1.34	17.74
P1A–O2A	1.30(3)	10.40(8)	1.26	13.03	1.28	12.64
P1A–O3A	1.65(3)	15.89(9)	1.57	24.44	1.57	24.63
P1A–O4A	1.70(3)	17.04(9)	1.57	24.34	1.58	24.51
P2B–O1B	1.42(3)	14.95(8)	1.42	17.35	1.36	18.61
P2B–O2B	1.30(3)	11.90(8)	1.27	12.80	1.28	12.49
P2B–O3B	1.62(3)	15.73(9)	1.54	22.92	1.55	23.11
P2B–O4B	1.67(3)	18.23(10)	1.58	25.45	1.59	25.64
O2A–C2A	1.700(19)	−8.88(6)	1.68	−11.65	1.67	−11.48
O2B–C2B	1.676(19)	−8.21(6)	1.68	−11.37	1.68	−11.52
O5A–C1A	2.74(3)	−31.88(14)	2.63	−15.03	2.65	−16.19
O6A–C1A	2.60(3)	−31.43(13)	2.55	−19.18	2.52	−17.90
O5B–C1B	2.78(3)	−32.86(14)	2.63	−15.89	2.65	−17.05
O6B–C1B	2.58(3)	−30.96(13)	2.54	−19.06	2.52	−17.97

---

Considering the ρ_c values, the general agreement between the experiment and theory is good. Certain discrepancies may be observed for the terminal P–O (P1A–O3A, –O4A, P2B–O3B and –O4B) and carboxylate C–O bonds, namely the experimental ρ_c values are slightly higher than the theoretical ones, but the differences are only slightly greater than the commonly accepted criterion of 3σ. In any case, as far as the three types of P–O bonds (i.e. the terminal ones, e.g. P–O2 or P–O3 in I; the hydroxyl P–OH like P–O1 in I; the ester P–OC like P–O4 in I) are concerned, the quantities obtained for the present compounds resemble the relevant values reported previously.^{43–45,48–52} The ester P–O bond in phosphoenolpyruvate (P–O4) is less populated (0.16 in e Å⁻³ terms of ρ_c) than in phosphoglycolate molecules (P1A–O2A and P2B–O2B). Contrary to that the ester C–O bonds are more populated in I, than in II, the difference being approximately 0.25 e Å⁻³. Out of other intramolecular bonds the carboxylate C–O ones merit a short comment. In I the carboxylic group is unionized –COOH, while in II these groups have lost their protons. Consequently, in I there is a difference of approximately 0.52 e Å⁻³ between the charge densities at the bcps of the formally double and single bonds (C1–O6 and C1–O5, respectively). This value locates between analogous differences for dihydrogen phosphoglycolate⁴⁵ (0.41 e Å⁻³), trans-cinnamic acid⁴⁷ (0.31 e Å⁻³), picolinic acid N-oxide⁵³ (0.55 e Å⁻³) or coumarin-3-carboxylic acid⁴⁷ (0.75 e Å⁻³). In II the C–O bonds are formally equal, with the bond order of 1.5. The different values of ρ_c result here from the influence of hydrogen bonds, namely atoms O6A and O6B are involved in strong H bonds, what brings about certain decrease of relevant charge density in comparison with ρ_c of C1A–O5A or C1B–O5B. All these quantities are within the range of the literature values.^54,55

In the case of Laplacian the question of agreement between the experimental and theoretical values is less straightforward. First of all it may be said that all the experimental quantities found for the present compounds are in reasonable agreement with the Laplacians reported so far for related bonds. Smaller or greater discrepancies appear if one compares the experimental data with the theoretical ones. It has been noticed that experimental and theoretical Laplacians have a tendency to differ, particularly in the case of polar bonds⁵⁶ and this is what may be observed in the present structures. In particular there is a good concurrence of both Laplacians for the C–C, C=C and ester C–O (i.e. C2–O4 in I, C2A–O2A and C2B–O2B in II) bonds. Also for the P–O_e, P–OH and O–H ones the differences are not serious. The discrepancy for C–H and N–H is larger, but it should be taken into account that these bonds are characterized by rather large negative Laplacians, that means that a small modification of the charge distribution near the bcp may give rise to a significant change of ∇²ρ_c. Anyway, the experimental and theoretical quantities found for the present systems are consistent with the previously reported.⁵⁶ The terminal P–O and carboxylic C–O bonds are therefore the most problematic. As it was previously commented in the literature (e.g. in the quoted paper) the difference results almost exclusively from the difference of Hessian eigenvalues associated with the eigenvectors parallel to the respective bonds (typically denoted as λ₃); such is also our case (see ESI† for the details). To get a deeper insight into the nature of these discrepancies we have calculated profiles of experimental and theoretical (the cluster) electron density, ρ, and its Laplacian, ∇²ρ, along the P1A–O3A and C1A–O5A bonds; they are presented in Fig. 3 and 4.

Fig. 3 Profiles of experimental (blue) and theoretical (dark red) Laplacian (upper diagram) and charge density (lower diagram) along the P1A–O3A bond. The x values show the distance from the phosphorus atom. The green line represents the difference between the experimental and theoretical charge density functions.

Fig. 4 Profiles of experimental (blue) and theoretical (dark red) Laplacian (upper diagram) and charge density (lower diagram) along the C1A–O5A bond. The x values show the distance from the carbon atom. The green line represents the difference between the experimental and theoretical charge density functions.

It may be noticed that the shape of the experimental charge density is almost perfectly reproduced by the theoretical function for the P–O bond; in the case of C–O the agreement is somewhat worse. Approximately 0.2–0.4 Å from the oxygen atom there seems to be a local increase of ρ_exp − ρ_theor; this is more evident in the P–O diagram. Its reason is not clear, but this feature illustrates slightly flatter experimental ρ near the minimum and steeper experimental ∇²ρ. The theoretical Laplacians, apart from being flatter, show similar shapes as the experimental ones. Large differences between both these functions at the bcps to a great extent reflect the fact that the critical points are located in the range of steep decrease of ∇²ρ.

Similar differences between the experimental and theoretical shapes of ρ and ∇²ρ have been recently reported for the N–N bond in pigment Yellow 110 molecule.⁵⁷

The discussed discrepancies may in part ensue from simplified (single Slater) form of the deformation radial functions used for fitting the experimental charge densities.⁵⁸

Theoretical study of PEP and PG

The NBO orbitals in the phosphate part of both molecules are presented in Table 7. Because the atom labels differ between I and II, in each line the chemically matching atoms and bonds are assembled (as far as it is possible; e.g. the ester oxygen atoms, O4 from I and O2A from II, or the hydroxyl phosphate oxygen atoms, i.e. O1 from I and O1A from II, etc.). For the sake of brevity only the data for molecule A from II are shown.

Table 7 The NBO orbitals involved in bonding in the phosphate group^a

H2PEP^−			HPG^2−(A)
Orbital(s)*	Occ.	Composition	Orbital(s)	Occ.	Composition
a lp – lone pair, dlp – diagonal lone pair, tlp – transversal lone pair; see text for the explanation.
σ(P–O1)	1.985	0.454P(sp^3.64) + 0.891O1(sp^1.86)	σ(P1A–O1A)	1.984	0.452P1A(sp^3.75) + 0.892O1A(sp^1.77)
σ(P–O2)	1.979	0.510P(sp^2.19) + 0.861O2(sp^2.19)	σ(P1A–O3A)	1.985	0.493P1A(sp^2.46) + 0.868O3A(sp^2.08)
σ(P–O3)	1.981	0.504P(sp^2.35) + 0.864O3(sp^2.18)	σ(P1A–O4A)	1.984	0.498P1A(sp^2.43) + 0.867O4A(sp^2.10)
σ(P–O4)	1.954	0.441P(sp^4.16) + 0.898O4(sp^3.05)	σ(P1A–O2A)	1.973	0.467P1A(sp^3.36) + 0.884O2A(sp^2.81)
lp(O1)	1.964	sp^1.34	lp(O1A)	1.968	sp^1.29
lp(O1)	1.924	p	lp(O1A)	1.930	p
lp(O2)	1.975	sp^0.46	lp(O3A)	1.975	sp^0.48
lp(O2)	1.818	p	lp(O3A)	1.858	p
lp(O2)	1.816	p	lp(O3A)	1.836	p
lp(O3)	1.974	sp^0.46	lp(O4A)	1.974	sp^0.48
lp(O3)	1.837	p	lp(O4A)	1.846	p
lp(O3)	1.820	p	lp(O4A)	1.835	p
lp(O4) (dlp)	1.933	sp^1.61	lp(O2A) (dlp)	1.944	sp^1.25
lp(O4) (tlp)	1.832	p	lp(O2A) (tlp)	1.891	p
σ*(P–O1)	0.220	0.891P(sp^3.64) − 0.454O1(sp^1.86)	σ*(P1A–O1A)	0.208	0.892P1A(sp^3.75) − 0.452O1A(sp^1.77)
σ*(P–O2)	0.119	0.861P(sp^2.19) − 0.510O2(sp^2.19)	σ*(P1A–O3A)	0.148	0.868P1A(sp^2.46) − 0.493O3A(sp^2.08)
σ*(P–O3)	0.143	0.864P(sp^2.35) − 0.504O3(sp^2.18)	σ*(P1A–O4A)	0.144	0.867P1A(sp^2.43) − 0.498O4A(sp^2.10)
σ*(P–O4)	0.256	0.898P(sp^4.16) − 0.441O4(sp^3.05)	σ*(P1A–O2A)	0.214	0.884P1A(sp^3.36) − 0.467O2A(sp^2.81)

---

Inspection of Table 7 leads to the following observations:

1. The populations of the purely p lone pairs in the two terminal and the hydroxyl oxygen atoms are depleted in comparison to the sp^x pairs, as expected, because the former ones are oriented perpendicularly to the P–O bonds and therefore their interactions with the σ*(P–O) orbitals are more favored;

2. The populations of the p lone pairs in H₂PEP⁻ are systematically lower than populations of their analogues in HPG²⁻, indicating thus stronger donor–acceptor interactions in the former molecule;

3. The population of the ester σ*(P–O4) orbital in H₂PEP⁻ is distinctly higher than that of σ*(P1A–O2A) in HPG²⁻ (in accord with previous findings for high-energy phosphates²).

4. The ester oxygen atom (O4 and O2A, respectively) has two lone pairs. The first one is oriented more or less along the diagonal of the C–O–P angle and for this reason will be hereinafter labelled as diagonal or dlp; its population is only slightly smaller than 2.0. The other pair is perpendicular to the COP plane and therefore will be identified as transversal or tlp; its population is reduced in a degree comparable to the other p orbitals of the terminal phosphate O atoms – more for H₂PEP⁻ than for HPG²⁻. A view of dlp and tlp for molecule A from II is shown in Fig. 5.

Fig. 5 (a) View of the diagonal lone pair of O2A. The view direction is perpendicular to the (P1A, O2A, C2A) plane; (b) a view of the transversal lone pair of O2A. The view direction is along the diagonal of the P1A–O2A–C2A angle.

Energy of hydrolysis and related questions

In this section we would like to discuss a few questions related with the hydrolysis energy of phosphate esters. Before we proceed to more specific issues a few assumptions will be made. First of all it will be assumed that the difference between the bonding energies of substrates and products, ΔE_h, is a good approximation of −ΔG^o_hydr. It has already been argued that the hydrolysis of organic phosphates is an isodesmic reaction⁵⁹ (i.e. a reaction in which the number and types of bonds do not change) and for this reason both these quantities do not differ significantly. Further on, it has been known that total ΔG^o_hydr of PEP (−14.8 kcal mol⁻¹) results from two consecutive processes: hydrolysis of phosphoenolpyruvate to enolpyruvate, and a subsequent tautomerization of the product to pyruvate. ΔG of the latter process is −6.9 kcal mol⁻¹;⁶⁰ by subtracting this quantity from −14.8 kcal mol⁻¹ the energy of the very hydrolysis step is −7.9 kcal mol⁻¹. Apart from that the following symbols will be used: O_t – phosphate terminal oxygen atom and R – a relevant correlation coefficient. Interactions of tlp or dlp with the organic fragments of the esters will be denoted as tlp-C or dlp-C, respectively. The following discussion will be kept within the limits of the NBO scheme.

The ester P–O and C–O bond lengths – factors affecting them and their relation with the energy of hydrolysis

There is a strong correlation between these two bond lengths (−0.950), what indicates that factors affecting these two quantities are essentially the same. A similar strong correlation was observed for experimentally determined P–O_e and C–O_e distances.⁶¹ For this reason we discuss these parameters jointly.

It has been suggested that there is a positive correlation between the length of this bond and the free energy of hydrolysis. We have tried to correlate both these quantities over the set of the 15 molecules specified above and the result is presented in Fig. 6.

Fig. 6 Correlation between the calculated ester P–O bond lengths and the energy of hydrolysis. The red points represent the molecules with a double bond vicinal to the ester O atom, the green ones – the trihalogenomethyl phosphates, the blue ones – monohalogenomethyl phosphates and the black ones – the remaining compounds. The black line represents the regression and R is the correlation coefficient.

The ester P–O bond length is somewhat less correlated with ΔE_h for the present data than it was reported previously by Ruben et al.² for a set of 10 molecules, for which the correlation coefficient between d(P–O) and −ΔG of hydrolysis was 0.88 (N.B. throughout this paper we take the hydrolysis energy as the difference between the bonding energies of substrates minus the energies of products, i.e. as the quantity related to the phosphorylation potential,⁶² therefore we have changed the sign of R accordingly). Inspection of the diagram presented in Fig. 6 reveals that there are two types of molecules departing from the general trend. Those in which the C atom bonded to the phosphate group is involved in a double bond (anions 1, 2, 4, 6) are located above the regression line, whereas the trihalogenomethyl phosphates (13, 14, 15) are underneath. Generally saying, four separate subgroups of compounds i.e. monohalides, trihalides, compounds with alkyl ester C atom and those with a double bond in the vicinity of the ester C atom can be discerned in the diagram.

The moderate correlation, illustrated in Fig. 6, indicates that the P–O_e bond length and the energy of hydrolysis are influenced by a somewhat different set of factors. On the other hand there is a very strong correlation between P–O_e bond length and the population or the energy of the σ*(P–O_e), as evidenced by the respective R factors of 0.995 and −0.999, what is in turn in good accordance with the results reported by Ruben et al. (0.98 and −0.99, respectively).

Accordingly, these authors argued that the donation from the terminal phosphate oxygen lone pairs to the σ*(P–O_e) orbitals – the so called generalized anomeric effect – should play a prevalent role in the scale of P–O bond lengthening or weakening. However, in the course of our investigation we have found a similarly strong correlation of P–O bond lengths with the orbital energy of the ester O atom transversal lone pair, tlp (R = −0.955). Assuming that the observed and calculated differences must originate from the chemical character and constitution of the residue bonded to the phosphate group (as within the latter there seems to rest no apparent reason to bring about variation of these bonds), we have decided to focus on the ester O atom and its interactions.

Indeed, a chain of correlations leading from the interactions of tlp with the ester C atom and its neighbourhood to the P–O_e distance can be drawn. A diagram of such correlations, showing a proposed logical sequence of factors, is presented in Chart 1. Interactions of lone pairs with empty antibonding orbitals are quantified here as the second order perturbation energy^21,22 calculated as

where q_i is the population of the occupied lone pair, F_ij – the element of the Fock matrix between the i-th occupied orbital and the j-th empty antibonding one, E_j and E_i – energies of the respective orbitals. The interactions of tlp with C (rectangle 1) are understood, accordingly, as the sum of E⁽²⁾ between tlp and the empty antibonding orbitals located in the vicinity of the ester C atom [among others π*(C–C) in H₂PEP⁻, π*(C–S) in thioacetyl phosphate, 6, σ*(C–Br) in 15 and so on].

Chart 1 A proposed scheme of relations that connect the interaction of the transversal lone pair of the ester oxygen atom, tlp, with the organic moiety and the lengths of the P–O_e and C–O_e bonds. The numbers represent the correlation coefficients between the respective parameters. Interactions are quantified in terms of relevant second order perturbation energies, E⁽²⁾, or their sums; see the text for more explanations.

The following interactions should be commented, when analyzing the presented chart:

(a) Interaction of the occupied orbital with an adjacent empty one (e.g. step 1 → 2 in the diagram). This leads to lowering of the energy of the occupied orbital. This should appear, in principle, as a result of the perturbations described above, but it should be kept in mind that the NBO orbitals are retrieved from the canonical DFT ones that have already included this effect.

(b) There may be a logically inverse relation: lowering of the energy of a considered orbital leads to weakening of its interactions with neighboring empty antibonding ones (step 2 → 3).

(c) Interaction of a lone pair located on A with σ*(B–C) leads to strengthening of the bond A–B (the case of 3 → 4 or 6 → 7 → 8) and to

(d) weakening of B–C (7 → 4).

To illustrate the donor–acceptor interactions referred to in p. (c) an exemplary overlapping of tlp with the π*(C=O) and one of the σ*(P–O_t) orbitals in the molecule of hydrogen acetyl phosphate (1) is shown in Fig. 7.

Fig. 7 Overlapping of tlp (blue), π*(C=O) (orange) and one of the σ*(P–O_t) (green) orbitals in hydrogen acetyl phosphate molecule. The upper drawing shows a projection approximately perpendicular to the C–O–P plane; the lower one represents a view roughly along the diagonal of the P–O–C angle.

It results from Chart 1 that interactions of the transversal lone pair of O_e with the ester carbon atom and its neighborhood may be regarded as the main factor shaping the P–O_e bond. As these interactions manifestly grow when low lying (in terms of energy) antibonding orbitals, such as π* or σ*(C–Br, Cl, F) are available (see Table 8), it may be noticed that the molecules with a double bond by the ester C atom, as well as the trihalogenomethyl phosphates, have the longest P–O_e bonds (the detailed values may be found in the ESI†). It may be also noticed that a few steps are needed to connect tlp-C interactions [expressed as the sum of second order perturbation energies, ∑E⁽²⁾(tlp-C)] with the length of the P–O_e bond and for this reason the correlation between these two parameters is somewhat blurred: the 1–4 correlation coefficient is 0.849; see Fig. 8 for the ∑E⁽²⁾(tlp-C)–d(P–O_e) diagram. The analogical correlation concerning the C–O_e bond is stronger – −0.913, what results from more direct connection between the C–O_e distance and ∑E⁽²⁾(tlp-C); see Chart 1. Loops 4 → 3 → 4 and 4 → 5 → 6 → 7 → 4 provide additional feedback that strengthens the bond elongation.

Table 8 Phosphate esters ordered according to interaction energies between tlp and the organic fragments, R-

Number from Scheme 1	R–	∑ΔE^(2)(tlp-C) (kcal mol^−1)	ΔEh (kcal mol^−1)
6	CH3–C(=S)–	46.08	7.5
2	CH3–C(=O)–	43.75	7.7
15	CBr3–	40.38	6.4
14	CCl3–	35.77	5.0
13	CF3–	33.30	4.1
1	HOOC–C(=CH2)–	31.69	4.0
4	CH3–CH2–C(=CH2)–	28.70	5.2
12	CBrH2–	25.00	3.8
10	CFH2–	21.56	2.9
11	CClH2–	21.31	3.5
8	HOOC–CH2–	12.09	1.5
5	CH2=CH–CH(–CH3)–	10.46	1.9
9	CH3–	10.37	0.9
7	CH2OH–CHOH–CH2–	10.21	1.6
3	Glucose-6-O-yl-	9.75	2.0

---

Fig. 8 Correlation between the sum of the second order perturbation energies, ∑E⁽²⁾(tlp-C), and the P–O_e bond lengths. The colors show the subgroups as in Fig. 6. For the purpose of interpretation the ∑E⁽²⁾ may be treated as the reason, and P–O_e length as a consequence. Likewise, in following figures the quantities shown in the x axes may be treated as causes, as only they may be reasonably interpreted as such.

In the case of the other lone pair of O_e it is not possible to draft a similar diagram. For instance, the correlation coefficient between the second order interactions of dlp with the organic moiety and the energy of this orbital (analogously to the step 1 → 2) is −0.296; the next step, equivalent to 2 → 3 gives R = −0.203 and so on. On the other hand, these interactions are strongly correlated with the C–O_e length (R = −0.908); this means that the parameters of dlp merely reflect shortening or lengthening of this bond and the orbital does not take essential part in the interplay of the organic moiety with the phosphate group.

Also the charge of the ester carbon atom does not seem to influence the P–O_e distance; the correlation coefficient between P–O_e and the Bader charge of C, Q_B, is 0.528. The multiple correlation coefficient, obtained from linear regression of P–O_e as a function of both ∑E⁽²⁾(tlp-C) and Q_B, is 0.876 – this means that the carbon charge does not affect the P–O_e bond length effectively.

Another factor, namely the steric interactions, merits analysis. In the NBO approach there are two ways of description of these interactions. The first consists in calculating the difference between the energies of orthogonalized and pre-orthogonalized natural localized molecular orbitals (NLMO) and the other – in calculating repulsions between pairs of disjoint orbitals.²² The correlation obtained by fitting the P–O_e distance as a function of ∑E⁽²⁾(tlp-C) and the steric repulsion is strong – 0.971 when using the first approach. The other method allows to separate the repulsions between various fragments of the molecule. In this way it may be found (see details in the ESI†) that the correlation is improved the most by the tlp-organic part repulsions (R = 0.944). Taking this into account and noting that the coefficient by the steric repulsion energies in the regression equation is positive, we are inclined to think that these interactions are not actually an independent factor that might influence the P–O_e distance (in such a case the coefficients should be negative reflecting thus a kind of braking force that would act oppositely to the tlp-C attraction), but they merely reflect lengthening or shortening of the C–O_e bond.

Eventually we have decided to consider the interaction of σ(P–O_e) with the organic part. Such an interaction would lead to creation of an additional weak bonding effect along C–O_e at the cost of the strength of the P–O_e bond. The values of ∑E⁽²⁾[σ(P–O_e)–C] are generally lower than ∑E⁽²⁾(tlp-C) and range between 2.8 and 4.2 kcal mol⁻¹ for compounds with aliphatic ester C, 4.0 and 6.3 kcal mol⁻¹ for monohalogenomethyl phosphates, 5.3 and 9.9 kcal mol⁻¹ for esters with a double bond at the ester C atom and between 11.5 and 16.8 kcal mol⁻¹ for trihalogenomethyl phosphates. These quantities are much lower than ∑E⁽²⁾(tlp-C), which are between 9.8 and 46.1 kcal mol⁻¹ (see ESI† for the details), but nevertheless the values for CX₃PO₃H⁻ (X = F, Cl, Br) may play some role in shaping the P–O_e bond. Interestingly, the coefficient of correlation between d(P–O_e) and ∑E⁽²⁾[σ(P–O_e)–C] is 0.954. Low values, however, of these interactions seem to suggest that they essentially mirror stronger interplay of tlps with the organic parts; only in the case of CX₃PO₃H⁻ they may exert considerable influence. This conclusion may be supported by the values of the R(∑E⁽²⁾(tlp-C), ∑E⁽²⁾[σ(P–O_e)–C]) correlation coefficient: 0.776 for all 15 molecules and 0.919 after exclusion of trihalogenomethyl phosphates. An examplary interaction of the σ(P–O_e) with a σ*(C–Br) in CBr₃PO₃H⁻ is illustrated in Fig. 9.

Fig. 9 Overlapping of the σ(P–O_e) (red and blue) with a σ*(C–Br) (yellow and green) in CBr₃PO₃H⁻.

Taking all these regards into account it must be said that the multiple regression of d(P–O_e) against both ∑E⁽²⁾(tlp-C) and ∑E⁽²⁾[σ(P–O_e)–C] gives a rather strong correlation coefficient, R = 0.970. The regression of d(P–O_e) against ∑E⁽²⁾(tlp-C) and ∑E⁽²⁾[σ(P–O_e)–C] is shown in Fig. 10.

Fig. 10 Relation between the P–O_e bond lengths calculated from regression against ∑E⁽²⁾(tlp-C) and ∑E⁽²⁾[σ(P–O_e)–C] and those obtained from DFT. The regression equation is: d(P–O_e) = 0.0008(4)∑E⁽²⁾(tlp-C) + 0.0069(11)∑E⁽²⁾[σ(P–O_e)–C] + 1.6119(59)‡ and the energies are in kcal mol⁻¹.

Similar reasoning may be drawn for C–O_e; in this case also the interaction of dlp with the organic part should be taken into account. As this interaction is strongly correlated with the interaction of σ(P–O_e)–organic part (the respective R = 0.970) both quantities may be treated jointly and the auxiliary interactions may be thus represented as the sum, ∑E⁽²⁾(aux) = ∑E⁽²⁾[σ(P–O_e)–C] + ∑E⁽²⁾(dlp-C)]. The correlation coefficients of interest are: R[d(C–O_e), ∑E⁽²⁾(tlp-C)] = −0.913, R[d(C–O_e), ∑E⁽²⁾(aux)] = −0.904, R[∑E⁽²⁾(tlp-C), ∑E⁽²⁾(aux)] = 0.805 and the correlation coefficient between the DFT data and the values obtained from the multiple regression against ∑E⁽²⁾(tlp-C) and ∑E⁽²⁾(aux) is 0.956. The regression results are shown in Fig. 11.

Fig. 11 Relation between the C–O_e bond lengths calculated from regression against ∑E⁽²⁾(tlp-C) and ∑E⁽²⁾(aux) and those obtained from DFT. The regression equation is: d(C–O_e) = −0.00196(54)∑E⁽²⁾(tlp-C) − 0.00193(58)∑E⁽²⁾(aux) + 1.4748(91) and the energies are in kcal mol⁻¹.

Energy of hydrolysis

Energy of hydrolysis, such as defined earlier, is strongly correlated with the interactions of tlp with the organic moiety: the correlation coefficient between this energy and ∑E⁽²⁾(tlp-C) is 0.959 and the relevant diagram is presented in Fig. 12.

Fig. 12 Correlation of the second order perturbation energy between tlp and the organic part of the phosphate esters, ∑E⁽²⁾(tlp-C), and the calculated hydrolysis energy, ΔE_h. The regression equation is ΔE_h = 0.154(6)∑E⁽²⁾(tlp-C). The constant term was less than 1.0 e.s.d and was therefore removed from the fitting.

This relation manifestly shows that the energy of hydrolysis is basically predetermined by the chemical constitution of the organic moiety; Table 7 presents the studied esters, ordered in descending order of magnitude of the tlp-organic moiety interactions. Inspection of this table reveals that the presence of low energy π* or σ* antibonding orbitals increases ΔE_h; only the former factor seems to be of interest for the biochemical systems.

Another issue arises here, namely, whether the hydrolysis energy may be somehow localized. This question perhaps should not be understood too literally, but nevertheless it would be appealing to identify orbitals which may contribute to release of this energy. In technical terms it means that a correlation of ΔE_h with the differences of the energies of the respective orbitals before and after the hydrolysis is sought. In the following discussion the quantities before hydrolysis will be labelled as pre and those after – as post. We have analyzed the following energies: energy of tlp [E(tlp)], energy of dlp [E(dlp)], total interaction of tlp or dlp with the organic fragment [∑E⁽²⁾(tlp-C) or ∑E⁽²⁾(dlp-C), respectively], steric repulsion energies of tlp or dlp [E_st(tlp) or E_st(tlp)], energy of σ(C–O_e) [symbolized as E(CO)], energy of σ(P–O_e) [E(PO)] and that of σ*(P–O_e). The post-hydrolytic energies (as well as the populations) of the two last orbitals were taken as averages of the respective quantities of the two P–O(H) bonds of the H₂PO₄⁻ anion. All the considered energy differences were weighted with the populations of the involved orbitals, according to the scheme: ΔE = P_preE_pre − P_postE_post, where P and E denote the applicable population and energy (in the case of Δ∑E⁽²⁾ the populations were halved).

The hydrolysis energy is correlated the most with the difference connected with dlp, Δ∑E⁽²⁾(dlp-C) (R = 0.852), but the chemical sense of such correlation may be difficult to interpret. For this reason we have decided to leave the Δ∑E⁽²⁾ parameters aside. In this way the best regression was obtained with three variables: ΔE(dlp), ΔE(PO) and ΔE(CO). The equation has the form (in kcal mol⁻¹​):

ΔE_h = 0.119(11)ΔE(dlp) + 0.077(11)ΔE(PO) + 0.125(23)ΔE(CO) − 3.75(65) (2)

and the correlation coefficient is 0.965. Curiously enough, the energy of tlp does not seem to contribute to the released energy. Substitution of ΔE(dlp) with ΔE(tlp) in the set of variables reduces the correlation coefficient to 0.513 (regression against ΔE(PO) and ΔE(CO) only gives R = 0.450). It is noteworthy that the ΔE(CO) values are generally negative (these values are particularly low in the case of trihalogenomethyl phosphates), hence the contribution originating from the C–O_e bond is destructive. The small constant term in the above equation may be a consequence of the fact that the total population of the three orbitals before and after hydrolysis is not strictly conserved: the post-hydrolytic population sum grows from 0.041 (compound 5) to 0.100e (13) in comparison to the pre-hydrolytic one. A diagram showing DFT vs. fitted values of ΔE_h is shown in Fig. 13.

Fig. 13 Correlation of the hydrolysis energies (ΔE) obtained from DFT and from linear regression with ΔE(dlp), ΔE(PO) and ΔE(CO) as the explanatory variables.

When analyzing eqn (2) it may be noticed that the relative contribution of dlp to ΔE_h, calculated as 0.119ΔE(dlp) is regularly greater than the share of σ(P–O_e) [computed as 0.077ΔE(PO)]. The ratio 0.119ΔE(dlp)/0.077ΔE(PO) is e.g. 6.0 for 4, 3.5 for H₂PEP (1), 2.8 for acetyl phosphate (2) and so on; only for trihalogenomethyl phosphates drops below 1.0. It may be therefore concluded that for most of the compounds majority of the hydrolysis energy is, in a sense, stored at dlp.

Relations of hydrolysis energy with structural and topological parameters of phosphate esters

It has been shown in the preceding paragraph that the energy of hydrolysis is primarily predetermined by the interactions of the transversal lone pair of the ester oxygen atom with the organic part of the molecule. These interactions are expressed as sums of relevant second order perturbation energies and, as such, cannot be experimentally measured. For this reason a question arises whether it is possible to estimate ΔE_h taking advantage of parameters experimentally derived from X-ray diffraction, such as e.g. bond lengths and/or charge densities at the bond critical points etc. To examine this problem we have taken into consideration the following parameters: bond lengths C–O_e and P–O_e, charge densities ρ_c(C–O_e) and ρ_c(P–O_e), together with their Laplacians, ∇²ρ_c(C–O_e) and ∇²ρ_c(P–O_e), as well as the Bader charge of the ester carbon atoms, Q_B(C). The analogous charge of O_e was disregarded, as it varies insignificantly in comparison with Q_B(C). The linear regression of ΔE_h against all the seven variables gave R = 0.978, but a few parameters could be safely eliminated because they were strongly correlated with others. Finally, it was found that, with a good approximation (R = 0.953), ΔE_h could be expressed as a function of three variables – ρ_c(C–O_e), ∇²ρ_c(C–O_e) and Q_B(C). The regression diagram is presented in Fig. 14 and the regression equation is:

ΔE_h = 189(23)ρ_c(C–O_e) + 40.5(64)∇²ρ_c(P–O_e) + 3.9(10)Q_B(C) − 26(3); (3)

ρ_c(C–O_e) and ∇²ρ_c(C–O_e) are in atomic units. For the density and Laplacian expressed in e Å⁻³ and e Å⁻⁵, respectively, the equation takes form:

ΔE_h = 28.1(34)ρ_c(C–O_e) + 1.68(27)∇²ρ_c(P–O_e) + 3.9(10)Q_B(C) − 26(3). (4)

Fig. 14 Correlation of the hydrolysis energy obtained from DFT and from linear regression against ρ_c(C–O_e), ∇²ρ_c(C–O_e) and Q_B(C).

However, a practical application of this equation may be limited, what can be shown when ΔE_h of the two compounds studied experimentally in this work are calculated. Using the parameters obtained from the multipole refinement one gets 4.4 kcal mol⁻¹ for H₂PEP⁻, what is quite a reasonable result. On the contrary, ΔE_h calculated for molecules A and B of PGH²⁻ is 7.3 and 7.6 kcal mol⁻¹ and it may be easily checked that such a large difference between the fitting and the DFT calculations results mainly from the difference between the experimental and theoretical values of the Laplacian in the critical point of the C–O_e bond. It seems, thus, that the precise determination of the Laplacian is crucial for realistic assessment of ΔE_h. Unfortunately, this is not always possible, as it has been discussed above (in the Topology section). Omitting ∇²ρ_c(C–O_e) yields a relation with R = 0.757. However, exclusion of the molecules with a double bond attached to the ester C atom greatly improves the regression results (R = 0.958; see Fig. 15).

Fig. 15 Correlation of the hydrolysis energy obtained from DFT and from linear regression against ρ(C–O_e) and Q_B(C) after elimination molecules with a double bond at the ester carbon atom.

The regression equation in this case is:

ΔE_h = −1.39(47)Q_B(C) + 50.0(59)ρ(C–O_e) − 9.8(15) (5a)

(atomic units) or

ΔE_h = −1.39(47)Q_B(C) + 7.41(88)ρ(C–O_e) − 9.8(15) (5b)

(e Å⁻³ for ρ(C–O_e))

The energy of hydrolysis, calculated from this relation is 2.6 and 2.5 kcal mol⁻¹ for phosphoglycolate molecules A and B, respectively.

Finally, it should be emphasized that the relations derived in this paragraph concern the mono-ionized phosphate esters. The question of whether and how the eqn (3)–(5) may be modified for phosphates at other ionization degrees needs a more extensive study.

Conclusions

Analysis of the experimental charge density distribution within the C–O–PO₃H fragments in HPG²⁻ and H₂PEP⁻ anions shows that the most essential differences between these two molecules are primarily manifested by topological parameters, namely ρ_c(C–O_e) and ρ_c(P–O_e). The first of these quantities is greater for phosphoenolpyruvate, while the latter – for phosphoglycolate. The differences between the atom charges within the phosphate group are less pronounced. There is a considerable difference between the experimental and theoretically calculated Bader charges of the ester carbon atoms in HPG²⁻ anions, the former being smaller. This is probably an effect of different parametrization of the radial functions; the differences between the experimental and theoretical values of certain Laplacians are another consequence of this reason. The DFT calculations show that the ester O atom has an approximate sp² hybridization, and one of its lone pairs, labelled in this paper as dlp, is located on the C–O–P plane and has the sp^x (x = 1.2–1.6) character, while the other, tlp, perpendicular to this plane, is located on the pure p orbital. The population of tlp in both compounds is depleted as compared to dlp. The depletion is larger in the case of H₂PEP⁻, what reflects larger interactions of this pair with the organic part of the molecule, in the first place with the antibonding π*(C–C) orbital. The calculations performed for 15 various model molecules allowed to get a deeper insight in the mechanisms determining the bond lengths within the C–OPO₃H fragment and to identify factors deciding about the energy of hydrolysis of organic phosphates. It results that the key role is played by tlp and its interactions with the carbon ester atom and its surrounding. This electron pair interacts with the phosphate group, mainly with the antibonding σ*(P–O_t) orbitals, on one side and with the organic part of the ester on the other. These two interactions are mutually competitive and the interactions with the organic part seem to be the dominating factor; its increase leads to lowering of the energy of tlp and in consequence to weakening of the tlp interactions with the phosphate group. As these interactions strengthen the P–O_e bond, the increased tlp-C interaction leads eventually to elongation and weakening of the P–O_e bond. This is in line with the previously proposed explanation, known as the “opposing resonance theory”.⁶³

Interaction of tlp with the organic part leads to depopulation of the pair, partial transfer of its electrons onto adjacent antibonding orbitals and in consequence to creation of additional bonding effect along the C–O_e bond. This effect is particularly strong if the energies of the available acceptor antibonding orbitals are low; this is the case of the π* and low lying σ*(C–X) orbitals (X = F, Cl, Br). This, in turn, results in shortening of the C–O_e bond in the trihalogenomethyl phosphates and in the molecules with the ester C atom involved in a double bond.

Additional donor interactions of σ(P–O_e) with the organic fragment also affect the lengths of the P–O_e and C–O_e bonds in the studied compounds; the latter bond seems to be affected by C-dlp interactions too. On the contrary, the hydrolysis energy depends almost exclusively on the tlp-C interactions; moreover it is proportional to them. For these reasons the P–O_e bond length shows only approximate relation to the hydrolysis energy.

The energy release during the hydrolysis is due to scission of the ester P–O bond and consecutive rearrangement of the molecular and electronic structure of the products. It has been attempted in this work to find whether it is possible to attribute the energy liberation to particular regions of the ester. The regression analysis suggests that there are three orbitals that chiefly contribute to the energy production; they are: dlp, σ(P–O_e) and σ(C–O_e). Two first ones experience stabilisation upon hydrolysis and therefore can be considered as those that release the energy, while the last one undergoes destabilisation at the cost of a part of the energy liberated by the former orbitals. Accordingly, the net ΔE_h is the relevant difference. Analysis of the regression equation suggests that dlp, although not taking a leading part in the interplay of the phosphate group with the organic moiety, may be treated as the major storage place of the hydrolysis energy, especially in the case of high energy phosphates with the double bond at the ester C atom.

This work presents analysis of the role of the donor–acceptor interactions of the ester oxygen lone pairs in determining the electronic and geometric structure as well as the thermodynamics of phosphate esters. Previously evidence of Kirby et al. indicated that the kinetics and scissile bond lengthening are strongly dependent on the effective electronegativity of the organic substituent, expressed as pK_a of the corresponding alcohol. We have come to the conclusion, that this effective electronegativity, or electron withdrawal force can be identified with the sum of the donor–acceptor interaction energies of tlp with the organic fragment of the ester.

Thus, following Kirby et al.^15–17 it may be inferred that in systems where strong delocalizing effects between tlp and organic parts exists, addition of a phosphate moiety is energetically less favourable, due to a high level of involvement of tlp in interactions with the organic part, and phosphate molecules thus formed experience greater destabilization with respect to the hydrolysis products.

We have also shown that, taking some precautions, especially connected with the divergent values of the experimental and theoretical Laplacians at the C–O_e bond critical point in some compounds, it is possible to estimate the hydrolysis energy from the topological parameters retrieved from the experimental charge density map.

Acknowledgements

One of us (AM) acknowledges the financial support of the Subaction 8.2.2, Action 8.2 of Priority VIII of the Operational Programme ‘Human Resources Development’, co-financed by the European Union, the budget of Poland and the Lubuskie Province, grant number DFS.VI.052.4.62.19.2013. The ADF calculations have been carried out at the Wrocław Centre for Networking and Supercomputing (http://www.wcss.wroc.pl), grant No. 58.

References

1. J. L. Berg, J. Tymoczko and L. Stryer, Biochemistry, W. H. Freeman and Co., Polish edition, PWN, New York and Basingstoke, 5th edn, 2002, pp. 373–382.
2. E. A. Ruben, J. A. Plumley, M. S. Chapman and J. D. Evanseck, J. Am. Chem. Soc., 2008, 130, 3349–3358.
3. F. Lipmann, Adv. Enzymol. Relat. Areas Mol. Biol., 1941, 1, 99–162.
4. J. R. Knowles, Annu. Rev. Biochem., 1980, 49, 877–919.
5. F. H. Westheimer, Science, 1987, 235, 1173–1178.
6. W. W. Cleland and A. C. Hengge, FASEB J., 1995, 9, 1585–1594.
7. N. Diaz and M. J. Field, J. Am. Chem. Soc., 2004, 126, 529–542.
8. M. Dittrich, S. Hayashi and K. Schulten, Biophys. J., 2003, 85, 2253–2266.
9. X. Lopez, M. Schaefer, A. Dejaegere and M. Karplus, J. Am. Chem. Soc., 2002, 124, 5010–5018.
10. P. P. Dzeja and A. Terzic, J. Exp. Biol., 2003, 206, 2039–2047.
11. H. M. Kalckar, Chem. Rev., 1941, 28, 71–178.
12. M. E. Rudbeck, S. O. N. Lill and A. Barth, J. Phys. Chem. B, 2012, 116, 2751–2757.
13. J. C. Summerton, J. D. Evanseck and M. S. Chapman, J. Phys. Chem. A, 2012, 116, 10209–10217.
14. Y. Yang and Q. Cui, J. Phys. Chem. A, 2009, 113, 12439–12446.
15. F. H. Allen and A. J. Kirby, J. Am. Chem. Soc., 1984, 106, 6197–6200.
16. A. J. Briggs, R. Glenn, P. G. Jones, A. J. Kirby and P. Ramaswamy, J. Am. Chem. Soc., 1984, 106, 6200–6206.
17. P. G. Jones and A. J. Kirby, J. Am. Chem. Soc., 1984, 106, 6207–6212.
18. D. G. Gorenstein, B. A. Luxon, J. B. Findlay and R. Momii, J. Am. Chem. Soc., 1977, 99, 4170–4172.
19. D. G. Gorenstein, B. A. Luxon and J. B. Findlay, J. Am. Chem. Soc., 1977, 99, 8048–8049.
20. K. Taira, T. Fanni and D. G. Forenstein, J. Am. Chem. Soc., 1984, 106, 1521–1523.
21. F. Weinhold and C. R. Landis, Discovering Chemistry with Natural Bond Orbitals, Wiley, Hoboken, New Jersey, Canada, 2012, pp. 96–105.
22. J. K. Badenhoop and F. Weinhold, J. Chem. Phys., 1997, 107, 5406–5421.
23. L. de Meis and V. A. Suzano, FEBS Lett., 1988, 232, 73–77.
24. H. Cheng, I. Nikolic-Hughes, H. J. Wang, H. Deng, P. J. O'Brien, L. Wu, Z.-Y. Zhang, D. Herschlag and R. Callender, J. Am. Chem. Soc., 2002, 124, 11295–11306.
25. P. K. Grzyska, P. G. Czyryca, J. Golightly, K. Small, P. Larsen, R. H. Hoff and A. C. Hengge, J. Org. Chem., 2002, 67, 1214–1220.
26. C. Peterhansel, K. Krause, H.-P. Braun, G. S. Espie, A. R. Fernie, D. T. Hanson, O. Keech, V. G. Maurino, M. Mielewczik and R. F. Sage, Plant Biol., 2013, 15, 754–758.
27. Å. Örström, Arch. Biochem. Biophys., 1951, 33, 484–485.
28. A. Weichsel and T. Lis, Pol. J. Chem., 1994, 68, 2079–2096.
29. T. Lis, Acta Crystallogr., Sect. C: Cryst. Struct. Commun., 1996, 52, 2326–2332.
30. R. H. Blessing, Acta Crystallogr., Sect. A: Found. Crystallogr., 1995, 51, 33–38.
31. N. K. Hansen and P. Coppens, Acta Crystallogr., Sect. A: Cryst. Phys., Diffr., Theor. Gen. Crystallogr., 1979, 39, 909–921.
32. A. Volkov, P. Macchi, L. J. Farrugia, C. Gatti, P. Mallinson, T. Richter and T. Koritsánszky, XD2006-a computer program for multipole refinement, topological analysis of charge densities and evaluation of intermolecular energies from experimental or theoretical structure factors, 2006.
33. F. A. Allen, O. Kennard, D. G. Watson, L. Brammer, A. G. Orpen and R. Taylor, J. Chem. Soc., Perkin Trans. 2, 1987, S1–S19.
34. A. Ø. Madsen, J. Appl. Crystallogr., 2006, 39, 757–758.
35. (a) P. J. Becker and P. Coppens, Acta Crystallogr., Sect. A: Cryst. Phys., Diffr., Theor. Gen. Crystallogr., 1974, 30, 129–147; (b) P. J. Becker and P. Coppens, Acta Crystallogr., Sect. A: Cryst. Phys., Diffr., Theor. Gen. Crystallogr., 1974, 30, 148–153.
36. F. L. Hirshfeld, Acta Crystallogr., Sect. A: Cryst. Phys., Diffr., Theor. Gen. Crystallogr., 1976, 32, 239–244.
37. E. J. Baerends, J. Autschbach, A. Bérces, C. Bo, P. M. Boerrigter, L. Cavallo, D. P. Chong, L. Deng, R. M. Dickson, D. E. Ellis, M. van Faassen, L. Fan, T. H. Fischer, C. Fonseca Guerra, S. J. A. van Gisbergen, J. A. Groeneveld, O. V. Gritsenko, M. Grüning, F. E. Harris, P. van den Hoek, H. Jacobsen, L. Jensen, G. van Kessel, F. Kootstra, E. van Lenthe, D. A. McCormack, A. Michalak, V. P. Osinga, S. Patchkovskii, P. H. T. Philipsen, D. Post, C. C. Pye, W. Ravenek, P. Ros, P. R. T. Schipper, G. Schreckenbach, J. G. Snijders, M. Solà, M. Swart, D. Swerhone, G. te Velde, P. Vernooijs, L. Versluis, O. Visser, F. Wang, E. van Wezenbeek, G. Wiesenekker, S. K. Wolff, T. K. Woo, A. L. Yakovlev and T. Ziegler, ADF2013.01, SCM, Theoretical Chemistry, Vrije Universiteit, Amsterdam, The Netherlands, 2008, http://www.scm.com.
38. J. P. Perdew, J. A. Chevary, S. H. Vosko, K. A. Jackson, M. R. Pederson, D. J. Singh and C. Fiolhais, Phys. Rev. B: Condens. Matter Mater. Phys., 1992, 46, 6671–6687.
39. (a) A. Klamt and G. Schüürmann, J. Chem. Soc., Perkin Trans. 2, 1993, 799–805; (b) A. Klamt, J. Phys. Chem., 1995, 99, 2224–2235; (c) A. Klamt and V. Jones, J. Chem. Phys., 1996, 105, 9972–9981.
40. E. D. Glendening, C. R. Landis and F. Weinhold, J. Comput. Chem., 2013, 34, 1429–1437.
41. J. C. Ortiz and C. Bo, Xaim – X atoms in molecule interface. Version 1.0., Universitat Rovira i Virgili, Tarragona, Spain, 1998.
42. V. Kinzhybalo, A. Mermer, T. Lis and P. Starynowicz, Acta Crystallogr., Sect. B: Struct. Sci., Cryst. Eng. Mater., 2013, 69, 344–355.
43. R. Janicki and P. Starynowicz, Acta Crystallogr., Sect. B: Struct. Sci., 2010, 66, 559–567.
44. A. Mermer and P. Starynowicz, Acta Crystallogr., Sect. B: Struct. Sci., 2011, 67, 399–408.
45. A. Mermer and P. Starynowicz, Acta Crystallogr., Sect. B: Struct. Sci., 2012, 68, 625–635.
46. P. Starynowicz and T. Lis, Acta Crystallogr., Sect. B: Struct. Sci., Cryst. Eng. Mater., 2014, 70, 723–731.
47. J. A. K. Howard, M. F. Mahon, P. R. Raithby and H. A. Sparkes, Acta Crystallogr., Sect. B: Struct. Sci., 2009, 65, 230–237.
48. E. Espinosa, C. Lecomte, E. Molins, S. Veinemillas, A. Cousson and W. Paulus, Acta Crystallogr., Sect. B: Struct. Sci., 1996, 52, 519–534.
49. N. Pérès, A. Boukhris, M. Souhassou, G. Gavoille and C. Lecomte, Acta Crystallogr., Sect. A: Found. Crystallogr., 1999, 55, 1038–1048.
50. M. Slouf, A. Holy, V. Petricek and I. Cisarova, Acta Crystallogr., Sect. B: Struct. Sci., 2002, 58, 519–529.
51. E. Aubert, F. Porcher, M. Souhassou and C. Lecomte, Acta Crystallogr., Sect. B: Struct. Sci., 2003, 59, 687–700.
52. M. Małecka, S. Mondal, S. van Smaalen and C. Paulmann, Acta Crystallogr., Sect. B: Struct. Sci., 2013, 69, 621–628.
53. A. V. Shishkina, V. V. Zhurov, A. I. Stash, M. V. Vener, A. A. Pinkerton and V. G. Tsirelson, Cryst. Growth Des., 2013, 13, 816–828.
54. E. Rödel, M. D. B. Messerschmidt and P. Luger, Org. Biomol. Chem., 2006, 4, 475–481.
55. L. Chęcińska, S. Mebs, C. B. Hübschle, D. Förster, W. Morgenroth and P. Luger, Org. Biomol. Chem., 2006, 4, 3242–3251.
56. A. E. Whitten, P. Turner, W. T. Klooster, R. O. Piltz and M. A. Spackman, J. Phys. Chem. A, 2006, 110, 8763–8776.
57. J. J. Du, L. Váradi, J. Tan, Y. Zhao, P. W. Groundwater, J. A. Platts and D. E. Hibbs, Phys. Chem. Chem. Phys., 2015, 17, 4677–4686.
58. J.-M. Gillet and T. Koritsanszky, in Modern Charge-Density Analysis, ed. C. Gatti and P. Macchi, Springer, 2012, pp. 181–211.
59. D. M. Hayes, G. L. Kenyon and P. Kollman, J. Am. Chem. Soc., 1978, 100, 4331–4340.
60. Y. Chiang, A. J. Kresge and P. Pruszynski, J. Am. Chem. Soc., 1992, 114, 3103–3107.
61. T. A. S. Brandão, J. Phys. Org. Chem., 2014, 27, 358–366.
62. A. L. Lehninger, Bioenergetics, the Molecular Basis of Biological Energy Transformations, W. A. Benjamin, Inc., German edition, Georg Thieme Verlag, New York, Amsterdam, Stuttgart, 1965.
63. E. A. Ruben, M. S. Chapman and J. D. Evanseck, J. Am. Chem. Soc., 2005, 127, 17789–17798.

---

Footnotes

† Electronic supplementary information (ESI) available: Parameters of the Slater deformation functions used in the multipole refinements, complete list of atom charges, full lists of topological parameters, results of the DFT calculations and regressions of the P–O_e distances obtained by inclusion of steric repulsion. CCDC 1403706 and 1403707. For ESI and crystallographic data in CIF or other electronic format see DOI: 10.1039/c5ra12605c

‡ In this and the following regression equations the estimated standard errors are given in parentheses and its order is the order of the last digit of the relevant estimated quantity, e.g. 1.6119(59) means that the e.s.d = 0.0059. A short comment on the statistical aspects of the regression analyses presented in this paper is given in the ESI.†

---