Experimental and Modeling Studies of Local and Nanoscale para-Cresol Behavior: A Comparison of Classical Force Fields

The dynamics of bulk liquid para-cresol from 340–390 K was probed using a tandem quasielastic neutron scattering (QENS) and molecular dynamics (MD) approach, due to its relevance as a simple model component of lignin pyrolysis oil. QENS experiments observed both translational jump diffusion and isotropic rotation, with diffusion coefficients ranging from 10.1 to 28.6 × 10–10 m2s–1 and rotational rates from 5.7 to 9.2 × 1010 s–1. The associated activation energies were 22.7 ± 0.6 and 10.1 ± 1.2 kJmol–1 for the two different dynamics. MD simulations applying two different force field models (OPLS3 and OPLS2005) gave values close to the experimental diffusion coefficients and rotational rates obtained upon calculating the incoherent dynamic structure factor from the simulations over the same time scale probed by the QENS spectrometer. The simulations gave resulting jump diffusion coefficients that were slower by factors of 2.0 and 3.8 and rates of rotation that were slower by factors of 1.2 and 1.6. Comparing the two force field sets, the OPLS3 model showed slower rates of dynamics likely due to a higher molecular polarity, leading to greater quantities and strengths of hydrogen bonding.

(a) (b) (c) Figure S1: A molecule of p-cresol with the atoms labelled with (a) the atom names according to table S1, (b) the atom charges of the OPLS2005 model and (c) the atom charges of the OPLS3 model.

S3
The experimental densities were measured on a DMA 4100 M density meter and showed close agreement with values found in the literature. [1][2][3] The experimental densities were only measured up to 363 K due to limitations of the instrument. These densities were compared to the densities of both simulated systems. The MD system densities were obtained from the average density after running each system for 2 ns in an NPT ensemble following an equilibration procedure. All the density values are listed in table S2. A plot of the variation in densities with temperature is displayed in figure S2.  Figure S2: The densities of liquid p-cresol.
The OPLS2005 system densities give a close match to the experimental densities, especially at higher temperatures. The greater density of the OPLS3 systems highlights potential inaccuracies with the intramolecular force-fields and/or the cresol atomic charges.

Quasielastic neutron scattering
Significant broadening of the QENS spectra is shown for p-cresol at 370 K across four values

Modelling localised and translational motions
For modelling localised motions, only confined motions and isotropic rotations of the molecules were considered in depth. Uniaxial rotations were considered unlikely, as the molecules are likely to rotate in multiple orientations in a liquid and methyl rotations were considered to rotate on a timescale that was too fast for the instrument to measure, considering previous data, 4 and may be partially accounted for by the flat background.
A cresol molecule translating within a confined space such as a dynamical basin can be described by a confined diffusion model, depicted in S4 (a). Confined translation assumes that the potential field within a symmetrical sphere is low compared to the infinite potential outside of it. Confined diffusion is indicated by a fixed Lorentzian half-width half-maxima (HWHM,∆ω(Q)) at low Q 2 , corresponding to a lack of movement over long distances, but follows Fickian diffusion at higher Q 2 . The diffusion coefficient (D s ) for confined diffusion S6 can then be calculated. In our results, the Lorentzian full-width half-maxima (FWHM = HWHM × 2) was plot against Q 2 .
Isotropic rotation, describes a molecule's random re-orientation around any axis, shown in figure S4   For purely diffusive motions, Fickian diffusion shows linear proportionality between the Lorentzian broadening and Q 2 , whereas jump diffusion begins to reach a plateau at higher Q 2 .

During jump diffusion the molecule translates a certain distance (d) between time periods (τ )
where the molecule appears static and it is vibrating around an equilibrium position. Jump diffusion often occurs in liquid systems with local order where the translation of molecules is periodically hindered by sterics or intermolecular interactions. Fickian diffusion is modelled below.
Various jump diffusion models can be used to fit suitable Lorentzian HWHM plots. The Hall and Ross (HR) 5 and Singwi and Sjölander (SS) 6 models represent jump lengths as a S7 continuous distribution developed for the characterisation of liquids. This is in opposition to the Chudley and Elliot (CE) 7 model which assumes discrete jump diffusion distances which was deemed unsuitable for largely disordered liquid diffusion, and so it is not discussed here.
All of the models assume that the jump time is negligible.
The HR model is as follows.
The jump lengths are distributed normally according to the equation below.
Whereas, the SS model is as follows.
Here, the jump lengths are characterised by an exponentially decaying distribution.

Experimental QENS analysis
The FWHM of Lorentzian 1 relating to translational diffusion exhibited linear proportionality at low Q 2 but reached a plateau as Q 2 increased, typical of a jump-diffusion process. Little difference was observed qualitatively in the goodness of fit for either the HR or SS model,

S8
shown in figure S5. However, the HR model was chosen on the basis of it having higher mean coefficients of determination (R 2 ) and lower χ 2 values. The R 2 values were calculated as 0.987 and 0.930 and the χ 2 values as 0.114 and 0.189 for the HR and SS models respectively, averaged across all temperatures.

Mean squared displacement
The D s for translational diffusion obtained from the MSD plots for each system applying the two force-fields were used to calculate their respective activation energies E a via the Arrhenius plot shown in figure S7. Figure S7: Arrhenius plot used to calculate the E a of translation calculated from the MSDs of simulated p-cresol systems applying OPLS2005 and OPLS3 models.

Incoherent dynamic structure factor
For the full details of the background, theory and calculations applied to gain the incoherent dynamic structure factors from the simulated trajectories, we direct the reader to the MDANSE manual version 1.0, through the following link: https://epubs.stfc.ac.uk/work/51935555 8 The calculated incoherent dynamic structure factors for the OPLS2005 and OPLS3 systems at 370 K were then fit using DAVE software, 9 with the fittings shown in figures S8 and S9. S10