Car–Parrinello and path integral molecular dynamics study of the intramolecular hydrogen bonds in the crystals of benzoylacetone and dideuterobenzoylacetone

The dynamics of the intramolecular short hydrogen bond in the molecular crystal of benzoylacetone and its deuterated analogue are investigated using ab initio molecular dynamics simulations. A study on intramolecular hydrogen bonding in 1-phenyl-1,3-butadione (I) and 1-deuteroxy-2-deutero-1-phenylbut-1-en-3-one (II) crystals has been carried out at 160 K and 300 K on the CPMD method level and at 300 K on the PIMD method level. The analysis of the two-dimensional free-energy landscape of reaction coordinate d -parameter and R O (cid:2)(cid:2)(cid:2) O distances shows that the hydrogen (deuter) between the two oxygen atoms adopts a slightly asymmetrical position in the single potential well. When the nuclear quantum eﬀects are taken into account, very large delocalization of the bridging proton is observed. These studies indicate that hydrogen bonds in the crystal of benzoylacetone have characteristic properties for the type of bonding model resonance-assisted hydrogen bonds (RAHB) without existing the equilibrium of the two tautomers. The infrared spectrum has been calculated, and a comparative vibrational analysis has been performed. The CPMD vibrational results appear to qualitatively agree with the experimental ones.


Introduction
Hydrogen bond is the most important of all directional intermolecular interactions. 1We can distinguish several groups of hydrogen bonds due to their properties.Most often, we share hydrogen bonds of the three main classes: strong, moderate and weak.Among the listed classes, the most important parameters are the length and angle of the bond, the energy of the bond, the position of the proton in the hydrogen bridge and the IR frequency shift of X-H stretching vibration.The most common classification of H bonds is the Jeffrey's classification. 2Proton transfer and hydrogen bonds are fundamental for the function, stability, structure and dynamics of chemically and biologically relevant systems.Very strong hydrogen bonds are among the most interesting features of hydrogen bond categories, which have essential roles in biochemical reactions and enzyme catalysis as transition state. 3,4In some specific molecular systems, for example, amides, 5,6 ketones, 7 pyrroles, 8 phenols, 9 proteins and nucleic acids, 10 the p-bond cooperativity phenomena were observed.These effects are explained by Gilli and Gilli [11][12][13][14][15][16] and coworkers who described and named the new class of hydrogen bonds in the series of b-diketones.This name, resonance-assisted hydrogen bond (RAHB), operates in the scientific environment, and is used when we talk about strong or very strong hydrogen bonds with a neutral donor and acceptor, which are connected by p-conjugated double bonds.][19][20][21][22] Special attention has been paid to the keto-enol tautomerism of the b-diketones, the structural properties of both keto and enol forms and the nature of the strong O-HÁ Á ÁO hydrogen bond in the enol form. 23urrently, in the literature debate continues regarding the shape of the potential well for the b-diketones, especially in the context of the occurrence of the enol forms of b-diketones and the related intramolecular enol-enol forms.In this regard, the scientific community is divided into two groups.Some researchers report a well with a single minimum [24][25][26] and others claim that it is a well with a double minimum function. 27,2825]29 Previous theoretical studies mainly on the DFT level on the benzoylacetone molecule in vacuum found the potential for proton transfer to be a double minimum with a barrier height of ca. 2 kcal mol À1 . 23,30After the addition the correction of the zero-point vibrational energy to the total potential energy, the internal barrier vanished and effective potential for the proton motion is a single one. 30n the present article, we have described the structures, proton motion, two-dimensional free-energy landscapes and the shapes of the O-HÁ Á ÁO potential function in benzoylacetone (BA) (I) and dideuterobenzoylacetone (II) (D 2 BA) in the solid state.Neutron diffraction data shows that in the solid state the enol structure of benzoylacetone with an almost symmetrical enol ring has been found. 24,25We have also discussed the influence of isotope effect on the shape of the free-energy function by analyzing the results of deuterium substitution obtained from MD simulations.To perform these tasks we used the PIMD method, 31,32 which takes into account the quantum nature of the proton and enables the thorough description of the movement of proton tunneling over the potential barrier.

Static calculations
To localize the key stationary points on the potential energy surface (PES) of the 1-phenyl-1,3-butadione (benzoylacetone) (I) in the solid state, a series of full geometry and cell optimizations with the London-type empirical correction for dispersion interactions, as proposed by Grimme, 33 together with the vibrational harmonic frequency calculations were undertaken.For the crystal, the structural data was obtained from the neutron diffraction study by Madsen and co-workers. 24,25Calculations were performed using the CRYSTAL09 program 34,35 utilizing the DFT method with PBE functional 36 with the two shrinking factors (7,7) to generate a commensurate grid of k-points in reciprocal space, following the Monkhorst-Pack method. 37Calculations were carried out with the consistent Gaussian basis sets of triple-zeta valence with polarization quality for solid-state calculations (pob_TZVP_2012), as proposed by Peintinger, Vilela Oliveira, and Bredow. 38For both the crystals, 1-phenyl-1,3butadione (benzoylacetone) (I) (BA) and 1-deuteroxy-2-deutero-1phenylbut-1-en-3-one (dideuterobenzoylacetone) (II) (D 2 BA), calculation of vibrational frequencies in CRYSTAL09 was performed at the G-point. 34,35

Dynamic and quantum dynamic simulations
Molecular dynamics calculations were carried out using the CPMD program, version 3.15.3. 39The initial molecular configuration for the 1-phenyl-1,3-butadione and 1-deuteroxy-2-deutero-1-phenylbut-1-en-3-one crystals was optimized by the preconditioned conjugate gradient (PCG) method with the London-type empirical correction for dispersion interactions, as proposed by Grimme 33 along with periodic boundary conditions (PBCs), which was employed for the solid state calculations.In these cases, real space Ewald summation of electrostatic interactions was carried out taking into account 8 cell replicas in each direction.We also generated the Monkhorst-Pack mesh (3,3,3) for calculated k-points in reciprocal space in each direction.The crystals data from the neutron diffraction study by Madsen et al. 24,25 were selected as the starting points.The crystals are monoclinic, benzoylacetone (P2  (0.072566 fs) coupled to a Nose ´-Hoover chains thermostat 40 at a frequency of 3200 cm À1 .The experimental values of unit cell parameters were used in the CPMD and PIMD simulations.An electronic mass parameter of 600 a.u. was employed.Electronic exchange and correlation have been modeled using the gradient-corrected functional of Perdew, Burke and Ernzerhof (PBE). 36Core electrons were treated using the norm-conserving atomic pseudopotentials (PP) of Troullier and Martins, 41 while valence electrons were represented in a This journal is © the Owner Societies 2014 plane-wave basis set truncated at an extended energy cut-off of 80 Ry.Following the initial equilibration period (about 8 ps), data was accrued for further 55 ps for crystals (I) and (II) Car-Parrinello at two different levels of temperature and for both crystals imaginary-time (55 ps) for the path integral dynamics simulation, carried out only at 300 K, respectively for eight (PI8) Trotter replicas (polymer-beads) using the normal mode variable transformation.The data was visualized using the VMD 42 and Gnuplot 43 programs with the path integral data first processed using a script by Kohlmeyer to calculate the centroid position of each set of polymer-beads. 44he vibrational spectrum was also calculated using the Fourier transformation of the dipole autocorrelation function obtained from dipole trajectories generated by the CPMD simulation facilitated by the scripts of Forbert and Kohlmeyer. 45It is important to point out that this approach includes anharmonic effects and all the vibrational modes for molecules in crystals.The lattice vibrations were not received because fixed experimental values of unit cell were used.

Structure analysis
The molecular geometry (in the unit cell) and atoms labeling of 1-phenyl-1,3-butadione (benzoylacetone, BA) in the solid state and its deuterated analogue are shown in Fig. 1.Red labeling was used for the second structure when two hydrogen atoms (H1, H5) in a single molecule has been substituted by deuterium.In Table 1, the calculated selected interatomic distances and angles after optimization compared with the average geometrical parameters from CPMD (standard deviation in brackets) and with the existing experimental data for both the crystals are presented.As can be inferred from Table 1, distances between heavy atoms in the hydrogen bond (O1-O2) in the studied crystals are very short (ca.2.5 Å) and have a similar length regardless of the method of analysis.We can also observe that the average lengths of the deuterated bridge in the analogue (II) are slightly longer than the average lengths of the H-bridge of the structure (I).][26]    geometrical parameters with the neutron diffraction experiment. 24,25owever, the geometry optimization on the static methods level It can be concluded that the equivalence of bond lengths may be attributed to delocalization over a conjugated p-bonded system in the keto-enol group. 15The quoted values of the bonds lengths clearly indicate different chemical bonds character.O2-C4, C2-C3 bonds are single and O1QC2, C4QC3 bonds have a double bond character.In addition, we have a relatively short hydrogen bridge and a non-ideal symmetrical position of the proton.All of these features fit well in the assumptions of the RAHB model proposed by Gilli. 15The large non-linearity of the hydrogen bond is noteworthy.Each of the methods of calculation indicates that the angle values of the hydrogen bridge are between 150 and 157 degrees, which are close to the experimental data. 24,25Moreover, the computational results are in good agreement with the experimental data of the dihedral angles, as shown in Table 1.In our opinion this confirms the correctness of the implementation of the simulation because  these parameters are particularly sensitive to structural errors. 46oreover, the optimized parameters of the unit cell are in reasonable agreement with the experimental values. 24,25The largest error (ca.10%) is noted for a parameter.For the remaining lattice parameters, errors between optimized and experimental values are smaller than 5%.This appears to confirm the need to use dispersion corrections to the full energy for calculations on the DFT method level. 47

Molecular dynamics simulations
The CP dynamic simulations were performed at two different temperatures because the experimental data for crystals are reported in this manner, 26 while PIMD simulations were conducted only at 300 K due to a very high computational cost.
The time evolution behavior of bond lengths directly involved in the selected intramolecular hydrogen bonds O2-H1(D1)Á Á ÁO1 from crystals (I) and (II) at 160 K and 300 K is presented in Fig. 2. The careful analysis of the strong O2-H1(D1)Á Á ÁO1 intramolecular hydrogen bond shows very large fluctuation of the bridging proton/deuter between the two oxygen atoms in crystals (I) and (II) at 160 K (see Fig. 2a and c) and significantly larger mobility of proton/deuter in hydrogen bridge in the simulations at 300 K (see Fig. 2b and d).As expected, higher temperature of simulation causes more frequent proton/deuter hops (transfer).At 160 K, proton/deuter jumps about 80 times while at 300 K this occurs almost three times more often.According to our estimates, the full process of proton hopping in the crystal of benzoylacetone takes about 12-16 fs, which is close to the experimental and theoretical data. 5,6,48,49When we analyze the results of the PIMD simulation taking into account the quantum effects, it can be seen  that the proton/deuter is completely delocalized between the two carbonyl oxygens (see Fig. 3).From the CPMD simulations at 160 K, the O2-H1(D1) bonds length changes in the ranges 0.953-1.895Å for the crystal (I) and 0.896-2.052Å for the crystal (II), while at 300 K the O2-H1(D1) bonds length changes in the ranges 0.945-1.887Å for the crystal (I) and 0.917-2.108Å for the crystal (II).Averaged O2-H1(D1) distances from the molecular dynamics method (see Table 1) are strongly dependent on temperature and their values are 1.164/1.161Å and 1.225/1.205Å for the crystal (I) and (II) at 160 K and 300 K, respectively.Following the distance of the hydrogen bridge, we observe that the O-H bonds are extremely sensitive to changes in temperatures.This phenomenon has been described by Wilson and co-workers 50 on the basis of single neutron diffraction studies for the benzoic acid crystal.The analysis of the distribution function of proton/deuter in the hydrogen bridge shows the importance of quantum effects.Their inclusion in the calculation completely changes the picture of physical processes.Fig. 4 and 5 present the distribution functions (d-parameter) of proton/deuter position along the reaction coordinate on the basis of our results from the CPMD and PIMD method.To show the influence of quantum effects in a better manner, we have prepared a series of graphs of free energy.Fig. 6 and 7 show the free energy profiles for proton motion obtained from the CPMD and PIMD results.Free energy profiles were calculated following the equation: where k is the Boltzmann constant, T is the simulation temperature, and P(d) is the proton distribution as a function of d (the reaction coordinate), which is defined as the difference between the distances r O-H(D) À r H(D)-O and is a measure of the proton/deuter transfer degree in the hydrogen bond.The value of zero indicated the midpoint of the hydrogen/deuter position in H/D-bond.The profiles of free energy in Fig. 6 demonstrate the absence of an effective barrier for proton transfer even without an inclusion of quantum effects for the CPMD simulation.Calculations show that the quantity of the barrier for proton transfer in the crystal of BA and its deuterated analogue are about 0.25 kcal mol À1 and decreases to about 0.15 kcal mol À1 as temperature increases.This would suggest that during proton transfer we observed the potentials wells with a double minimum; however, when the quantum effects are taken into account, the effective potential shape drastically changes, as shown in Fig. 7a  and b.These strong HBs have also been named low-barrier hydrogen bonds (LBHB). 51It is interesting to note that the results of CPMD simulations (the double well potential) are very similar to the static DFT calculations of BA molecule in vacuum, 23,30 whereas the CPMD simulations are in very good agreement with the experimental results showing almost symmetrical H-bond. 24,25or a deeper analysis of proton mobility in the studied intramolecular H-bonds, we constructed two-dimensional free-energy landscapes from two-dimensional distribution of d-parameter and OÁ Á ÁO distances, as illustrated in Fig. 8, Fig. 1S-3S (in ESI †), for every hydrogen bridge in the crystal unit, for further analysis for the proton/deuter transfer.As shown in the series of figures from 8 to 3S (ESI †), the panels Let us now briefly discuss the abovementioned landscapes.As shown in sections a, b, d and e of figures, the free energy surfaces obtained from CP simulations have two minima, while sections (c) and (f) obtained from PIMD simulations have one minimum in each case.As shown in Fig. 7a and b, the nuclear quantum effects remove the proton transfer barrier, allowing barrier-less proton transfer.Considering the two-dimensional free-energy landscapes, we can also see a substantial mobility and relocation of proton position at higher temperatures of simulations and negligible impact of isotope effect on the shape and height of the energy barrier of proton/deuter transfer.Moreover, despite the fact that the potential well has one minimum, the position of the proton/deuter is not perfectly symmetrical.It is a very subtle effect, which is often accompanied by typical resonance-assisted hydrogen bond.
According to Gilli and Gilli, 15,16 for the description resonanceassisted hydrogen bonded systems, the very important so-called p-delocalization index is defined as: where the . In our studies, the delocalization parameters (Q) for BA crystal and its deuterated analogue were described by relevant molecular distances Q = (O2-C4 -O1-C2) + (C2-C3 À C3-C4).According to Gilli and co-workers, 15 Q = 0 corresponds to the fully p-delocalized structure and Q = À0.32 or +0.32 Å corresponds to the completely localized bonding.The L value of 0.5 indicates a fully delocalized keto-enolic system, whereas L value closer to 0 or 1 indicates localized bonding.Two-dimensional plots of delocalization parameter Q versus R O1Á Á ÁO2 distance for crystals (I) and (II) are shown in Fig. 9. Values of L and Q are calculated on the basis of average bond lengths from CPMD simulation.In both the cases, the L-parameter is equal to 0.445.After a precise analysis of Fig. 9, we can conclude that the intramolecular hydrogen bond in the crystal of BA can be classified as the group of hydrogen bond with RAHB character.In addition, Fig. 9 allows the observation of the correlation between the Q and L parameters and the length of the hydrogen/deuter bridge throughout the course of the simulation.Therefore, according to the classification proposed by Gilli and co-workers, 15 it can be concluded that changes in the hydrogen and other bonds coupled with the molecular system are characterized by the type of strong p-delocalized structure (enol-enol forms).However, it should be noted that in this disordered molecular system there are two forms with slightly more contribution from the p-localized enol-keto (EK) form than p-localized keto-enol (KE) form.We expect that these mainly result from a little asymmetrical position of the proton/ deuter in the bridge.The calculated and experimental value of the O1Á Á ÁO2 distance close to 2.5 Å might suggest strong O-HÁ Á ÁO intramolecular bond.Estimation of the strength of intramolecular H-bonds for systems in a crystal is not trivial.One of the simplest methods for determining H-bond strength in vacuum is the calculation of the difference in energy between closed (X-HÁ Á ÁY) and open configuration.However, this approach cannot be applied for a system in a crystal.Therefore, for estimating this intramolecular H-bond, we used the very simple relation proposed by Musin and Marian: 52 E HB = (5.554Â 10 5 ) exp(À4.12R), ( where R is the OÁ Á ÁO distance.Using the value of 2.5 Å, which is close to the experimental and average values from CPMD simulations, the estimated H-bond strength is equal to À18.7 kcal mol À1 .The strength of H-bond in a crystal estimated by us is in good agreement with the value of À16.3 kcal mol À1 obtained in vacuum on the basis of DFT calculations by proposed by Schiøtt and coworkers. 30It implies that this H-bond is strong in the group of intramolecular hydrogen bonds.

IR spectra
For crystal (I) and (II), selected characteristics of vibrational frequencies from CPMD, together with harmonic frequencies from static calculations and available experimental data, are presented in Table 2. Fig. 10 presents the CPMD infrared spectrum for both the examined crystals BA and D 2 BA with experimental spectra for BA crystal. 53,54For further comparative vibrational analysis, we also compared the simulated and experimental spectra for BA crystal and its deuterated analogue, as illustrated in Fig. 11.The assignments of the calculated frequencies are based on the experimental IR 53 spectra for benzoylacetone and its deuterated analogue.As we know, infrared spectroscopy is a very useful experimental method for studying the strength of the hydrogen bond and it is also very important because it confirms that the results of molecular dynamics simulation are correct.In the context of this paper, the analysis of the vibration frequency of -OH groups appears to be the most important.6][57][58] The intensity and broadness of this band are dependent on the strength of the intramolecular hydrogen bond. 53We have observed that the stretching vibrations -OH/-OD are shifted in comparison to the -OH/-OD frequencies of bonds not involved in H-bond to lower frequencies of about 2217/1681 cm À1 and 2514/1922 cm À1 for crystals (I) and (II), respectively, in both the calculation methods.The calculated stretching vibrations -OH/-OD are in good agreement with experimental results (2650/1960 cm À1 ) but  it appears that the CP molecular dynamics method reproduces the frequencies of vibrations better than static methods in harmonic approximation in every case.As compared with the experiment, most of the frequencies from CPMD methods are slightly shifted toward a lower frequency; however, no clear tendency can be observed because some frequencies are a slightly higher.As expected, we observed very large red shift effects of -OH stretching vibrations in these molecular systems.This phenomenon is typical for the incidence of enol-enol resonance structure with the strong intramolecular hydrogen bond.
A very important parameter related with the IR spectra of H-bonded systems is the isotope effect defined by the isotopic ratio (ISR) nOH/nOD.Because of the anharmonicity in a majority of H-bonded systems, the ISR value is different from O2.It has been shown that in the case of O-HÁ Á ÁOH-bonds, the maximum value of ISR does not exceed 1.5. 59The calculated ISR values equal 1.32, 1.31, 1.35 from the static calculations, CPMD simulations and experiment, respectively.These values also indicate a very strong hydrogen bond and high mobility of proton in the hydrogen bridge.We observed a good agreement of the experimental data (3110/2313 cm À1 ) with the calculated (3186/2361 cm À1 and 3107/2336 cm À1 ) values for -CH/-CD stretching vibration frequency.In this case, we can also see a very large isotopic effect for these vibrations.The OÁ Á ÁO stretching mode for BA in the experiment appears at about 396/384 cm À1 .According to our static and CP calculations, these modes are located at about 410/395 cm À1 and 345/279 cm À1 .The nOÁ Á ÁO mode is coupled with dC-CH 3 and the position of this band relates to other b-diketones, which also suggests the existence of very strong HB.

Conclusion
We have presented the results of theoretical studies on the intramolecular hydrogen bonds in 1-phenyl-1,3-butadione and its deuterated analogue 1-deuteroxy-2-deutero-1-phenylbut-1-en-3-one in the solid state using static as well as ab initio molecular dynamics methods.In our study of both the crystals, we found evidence for the existence of strong p-delocalized structure (enol-enol forms) with the location of the enol hydrogen in a very flat slightly asymmetric (almost symmetric) single minimum potential, which is in accordance with neutron diffraction studies performed by Madsen 24,25 and co-workers.The analysis of average distances in the geometric structure of four bonds in the keto-enol parts of the molecule indicated that the equivalence of the bond lengths may be attributed to delocalization over a conjugated p-bonded system in the keto-enol group.We also calculated the p-delocalization index (which is equal to 0.445) for the both crystals taking as the basis for calculating the average bond lengths from CPMD simulation.In addition, the calculated correlation between the Q and the L parameters and the length of the hydrogen/deuter bridge throughout the simulation shows that the system is characterized by the type of strong p-delocalized structure (enol-enol forms) with slightly more contribution from the p-localized enol-keto (EK) form than p-localized keto-enol (KE) form.In our opinion, all the arguments presented above (a relatively short hydrogen bridge, non-ideal symmetrical position of the proton, the shape of the potential function of proton, the correlation between the Q and L parameters and very big red shift effects of -OH stretching vibrations) indicate that the hydrogen bonds in the benzoylacetone crystals have properties characteristic of the type of bonding model resonance-assisted hydrogen bonds without the existing equilibrium of the two tautomers.Fig. 11 Comparison of the experimental and simulated IR spectra for benzoylacetone in the solid state.The numbers of spectra were assigned as follows: (1) benzoylacetone and (2) dideuterobenzoylacetone in the solid state from the DFT (PBE) harmonic calculation with London-type empirical correction for dispersion interactions, (3) benzoylacetone and (4) dideuterobenzoylacetone in the solid state from CP molecular dynamics simulation at 300 K and number (5) experimental IR absorbance spectrum 53,54 for the benzoylacetone in the solid.
gratefully acknowledge the Academic Computer Centre in Gdansk (CI TASK) for the use of the Galera-ACTION Cluster and the Wroclaw Centre for Networking and Supercomputing (WCSS) for the use of the Supernova Cluster.Dr Przemysław Dopieralski is acknowledged for helpful discussion as well as help with preparing the AWK's scripts.

Table 1
Calculated selected geometrical parameters after optimization compared with the average geometrical parameters from CPMD (standard deviation in brackets) and the existing experimental data for benzoylacetone (I) and dideuterobenzoylacetone (II) crystals (bond lengths and unit cell parameters are in Å, angles and dihedrals are in degrees)
slightly underestimates the length of the hydrogen bridge of about 0.051 Å for the structure (I).For analyzing the geometric structure of the molecule in terms of the prevalence of characteristics model RAHB, we should carefully consider the length of the four bonds in the keto-enol parts of the molecule.Taking into the account our simulations carried out on the static and especially on dynamic levels, we can conclude that the bond lengths (O1QC2/O2-C4 and C4QC3/C2-C3) coincide with typical lengths for the bonds determined experimentally.For example, let us analyze the data from both simulations for the crystal (I), static optimization on the DFT method level and molecular dynamics in terms of Car-Parrinello.The resulting values are presented in order O1QC2/O2-C4 (1.286/1.336Å) and C4QC3/C2-C3 (1.388/1.418Å); O1QC2/O2-C4 (1.290/1.316Å) and C4QC3/C2-C3 (1.406/1.415Å) for DFT and CPMD methods, respectively.

Fig. 4
Fig. 4 Comparison of the distribution functions (d-parameter) from CPMD simulations for the intramolecular H-bonds for (a) benzoylacetone at 160 K, (b) benzoylacetone at 300 K for (c) dideuterobenzoylacetone at 160 K, and for (d) dideuterobenzoylacetone at 300 K.

Fig. 5
Fig. 5 Comparison of the distribution functions (d-parameter) from PIMD simulations for the intramolecular H-bonds for (a) benzoylacetone at 300 K and (b) dideuterobenzoylacetone at 300 K.

Fig. 6
Fig.6Comparison of the single proton transfer free-energy DF profiles from CPMD simulations for the intramolecular H-bonds for (a) benzoylacetone at 160 K, (b) benzoylacetone at 300 for (c) dideuterobenzoylacetone at 160 K, and for (d) dideuterobenzoylacetone at 300 K.
(a) and (b) presented the free energy surfaces obtained for crystal (I) from CPMD simulations and panel (c) from PIMD simulations.Similarly, the sections (d)-(f) of figures were obtained for the second crystal (II).

Fig. 7
Fig. 7 Comparison of the single proton transfers free-energy DF profiles from PIMD simulations for the intramolecular H-bonds for (a) benzoylacetone at 300 and for (b) dideuterobenzoylacetone at 300 K.

Fig. 8
Fig. 8 Two-dimensional free-energy landscape of d-parameter (reaction coordinate) and R O1Á Á ÁO2 distances for crystal (I) for (a) CPMD simulation at 160 K, (b) CPMD simulation at 300 K, (c) PIMD simulation at 300 K, and for crystal (II) for (d) CPMD simulation at 160 K, (e) CPMD simulation at 300 K, (f) PIMD simulation at 300 K.The unit of DF free energy (potential of mean force) is kcal mol À1 .

Fig. 9
Fig. 9 Two-dimensional p-delocalization index (L) landscape of R O1Á Á ÁO2 distances and delocalization parameter (Q) for crystal (I) for (a) CPMD simulation at 160 K, (b) CPMD simulation at 300 K, and for crystal (II) for (c) CPMD simulation at 160 K, (d) CPMD simulation at 300 K.