Ultrasonic study and molecular simulation of propylene glycol at pressure up to 1.4 GPa

We report an ulsrasonic measurements of density and bulk modulus of propylene glycol at room temperature and at the temperature of liquid nitrogen combined with molecular dynamics simulations with two different force fields. We find that experimental density of propylene glycol at room temperature is well described within COMPASS force fields simulations, while the bulk modulus from simulation deviates from the experimental one. Number of hydrogen bonds in propylene glycol is also evaluated.


INTRODUCTION
Propylene glycol (PG) is an important liquid for many technological applications, including food industry, resin production, pharmacy, etc. It is also widely used as antifreeze. Al of this makes understanding of the properties of PG to be of a great importance. And indeed there is a plenty of publications devoted to different properties of PG. At the same time there are still a lot of properties of PG which are still badly understood.
Propylene glycol has a chemical formula C 3 H 8 O 2 . It is a molecular liquid, which can be easily transformed into glass (the glass transition temperature at ambient pressure is T g = 165 K). A molecule of PG contains two hydroxyl groups, which makes it possible to form hydrogen bonds (H-bonds) in PG. Although H-bonds are weaker than other types of bonding they can strongly affect the properties of a liquid [1][2][3][4]. H-bonds can affect the orientational correlations between molecules of a liquid and glasses. The relaxational dynamics of H-bonds strongly affects the dynamical properties of liquids. Hbonds also can be responsible for the properties of a liquid as a solvent, as it is assumed in the case of water [5]. Apparently, a lot of attention of researchers was paid to investigation of H-bonding in different liquids, both experimentally and theoretically. For instance, in Refs. [6][7][8][9][10][11][12][13][14] experimental study of H-bonds in liquid alcohols was reported. These works were complemented by molecular dynamics simulation [15][16][17]. However, the problem of the behavior of H-bonds still remains unclear, since different works give different, sometimes opposite results (see, for instance, Refs. [18][19][20][21][22] for hydrogen bonds in methanol).
In the present work we perform ulstrasonic measurements and molecular simulations of PG at room temperature and at the temperature of liquid nitrogen up to the * Corresponding author: fomin314@mail.ru pressure of 1.4 GPa. We measure the equation of state, compressibility and its baric derivative and compare the results to the ones from computer simulation. Molecular simulation also allows us to monitor microscopic structure of PG and evaluate the number of H-bonds in a wide range of pressures.

Experimental methods
The experimental part of the work was carried out using a high-pressure ultrasonic piezometer. It was created on the basis of a press and a piston-cylinder high-pressure chamber, placed in a thermostat for low-temperature measurements [23,24]. This setup makes it possible to carry out experiments up to pressures of 2 GPa and in the temperature range of 77-310 K both in the isothermal compression and isobaric heating regimes. In this work, we performed isothermal compression up to 1 GPa at constant temperatures of 77 K (liquid nitrogen temperature) and 295 K (room temperature). The pressure measurement error was no more than 0.05 GPa, that is, 5%, and the temperature was controlled with an accuracy of 1 K. To measure the temperature, we used 4 copper-constantan thermocouples: 2 were glued in the immediate vicinity of the sample and 2 more were glued to punches that pressed on sample, to control the temperature gradient. Since propylene glycol is a liquid under normal conditions, we used capsules to contain the substance. The capsule is a thin-walled Teflon cylinder with an outer diameter of 18 mm and an inner diameter of 16 mm, which is closed on both sides with copper caps 1 mm thick and a rubber ring as a seal. Propylene glycol 99.5% pure was purchased from Sigma-Aldrich. The capsule with the substance was placed in a high-pressure chamber and compressed from both sides by punches, on the ends of which piezoplates of lithium niobate LiN bO 3 were glued. They were both generators and receivers of the ultrasonic signal. For longitudinal and transverse ultrasonic waves, different plates were used, with different thicknesses and different resonant frequencies: 5 MHz for transverse and 10 MHz for longitudinal ultrasonic waves. The transit time of ultrasound through the sample was measured with an accuracy of 1 ns. The length of the sample was also measured -for this we used 2 micrometers with a measurement accuracy of 5 µm. From the length and transit time, the ultrasound velocity was obtained, then, in the approximation of a homogeneous isotropic medium [25,26] the adiabatic bulk modulus B s was calculated: where ρ is the density of a substance, v l is longitudinal ultrasonic wave velocity and v t is the transverse ultrasonic wave velocity.
In Supplementary materials we describe some details of the experimental derivation of density of PG.

Molecular simulation
In the computational part of the work we simulated a system of 1000 molecules of PG in a cubic box with periodic boundary conditions. Two different force fields (FF) were used: COMPASS [27] and Charmm [28]. Figure 1 shows a molecule of PG with notation of all atoms.
In both FFs the system was equilibrated for 1 · 10 7 steps with the time step dt = 0.1 fs. Then the system was simulated for 2 ·10 7 step for calculation of its properties. The bonds involving hydrogens were constrained via Shake algorithm [29]. Firstly the system was simulated at temperature T = 295 K at several densities in canonical ensemble (constant number of particles N, volume V and temperature T). Equation of state was calculated and compared to the experimental data. This equation of state was approximated by polynomial functions to calculate the isothermal bulk modulus B T = ρ dP dρ T and the isothermal baric derivative dB dP T . The details of approximation are given in Supplementary materials.
When the simulation at T = 295 K was finished the system was quenched to a low temperature T = 77 K. The cooling was made with within 5 · 10 7 steps, i.e. with the cooling rate 4.36 · 10 9 K/s. Equation of state (EoS), bulk modulus and baric derivative were calculated at the low temperature too. Figure 2 shows a definition of hydrogen bonds in the system. We assume that a hydrogen bond is formed if two conditions are fulfilled: i) the distance r < 3.5Å and (ii) the angle α < 30 o .
The structure of the system was characterized by center of mass radial distribution functions (RDFs) and par- Hydrogen bonding in PG is discussed. The presence of H-bonds was analyzed based on geometrical criterion. Several definitions of H-bonds in water were discussed in Ref. [30]. In the present paper we adopted the r − α definition from this paper, i.e. we assume that an Hbond is formed if the distance r is below some cut-off value and the angle α is below some cut-off angle (see Fig. 2). Following Ref. [31] we choose r cut = 3.5Å and α cut = 30 o .
All simulations were performed using the LAMMPS simulation package [32]. Figure 3 shows a comparison of experimental data with molecular simulation with two different force fields -COMPASS and Charmm. One can see that at room temperature COMPASS force field perfectly coincide with experimental data up to the pressure of 0.5 GPa and stays within the error bars of experiment up to about 0.8 GPa. At the same time one can see that starting from P ≈ 0.5 GPa the experimental density is higher than the one from simulation. Several reasons can induce this disagreement. In particular, the ultrasonic measurements of PG were carried out both on compression and decompression modes. The intrinsic friction of the installation leads to overestimation of the pressure in the compression cycle and underestimation in the decompression mode. As a result we observe a hysteresis loop which is eliminated by averaging of the results. However, this method works accurately in the pressure interval of 0.2 − 0.7 GPa only. Errors can appear due to the uniaxial compression in low pressure regime and due to very strong friction in the high pressure one. Thus, the discrepancy between the theoretical and experimental density at high pressures can be explained by the increased experimental error in pressure measurements. Basing on this reasoning we may suppose that experimental densities are slightly overestimated at high pressures, and the ones from COMPASS FF modelling are closer to the real ones. Unfortunately, we are not aware of any measurements of EoS of PG at elevated pressures by another experimental method and therefore we cannot verify our assumption by direct comparison of the data.

RESULTS AND DISCUSSION
The results of molecular simulations with Charmm FF strongly underestimate the density at T = 295 K. The disagreement is about 10%, which we consider as unsatisfactory result.
At the same time one can see that EoS obtained by COMPASS FF at T = 295 K has larger derivative than the one from experiment, while the slope of Eos obtained by Charmm FF is close to the experimental one. Figure  4 (a) shows a comparison of bulk modulus B T obtained withing the framework of these two force fields with the experimental data for the modulus B S . Although we calculate isothermal modulus in case of molecular simulation and adiabatic in case of experimental work, as it was shown in Supplementary materials of Ref. [24], these moduli are almost indistingluishable at T = 77 K and deviate from each other within 10% interval at room temperature. For this reason such a comparison is possible. One can see almost perfect agreement of the data from Charmm FF with experimental data for the pressures above 0.2 GPa, while COMPASS FF overestimates B T comparing to the experimental data.
Experimental data for bulk modulus suggest that B T is almost constant withing the studied interval of pressure and takes value close to 8.0. This value of baric derivative corresponds to LJ system. Basing on this reasoning one can assume that the interaction between PG molecules can be roughly approximated by the LJ potential.
The baric derivatives from molecular simulation demonstrate more complex behavior (see Fig. 4 (b)). Although baric derivative from simulation is not constant for the case of COMPASS FF at T = 295 K, but demonstrates assymptotic tendency to the value of 8.0. From this reasoning one can expect that the interaction potential of PG molecules becomes more LJ-like with pressure. In order to check the assumption above we calculate the RDFs of centers of mass of PG molecules along T = 295 K isotherms with COMPASS FF. From this figure one can see that RDF of centers of mass of PG molecules cannot be described by LJ model by two reasons: (i) there is a pre-peak next to the first peak of g(r) and (ii) the first peak of g(r) demonstates a shoulder at the right-hand side of the peak. At the same time one can see that both effects become lower as the pressure increases. Basing on this one can expect that further increase of the pressure may lead the center of mass RDF of PG to the shape character that the RDF of LJ system.
Finally we are going to discuss the hydrogen bonding in PG. Formation of a hydrogen bond requires a polar hydrogen. In PG molecule it corresponds to the H 12 and H 13 atoms, which are bonded to oxygens within hydroxyl groups (see Fig. 1). The hydrogens belonging to CH 3 groups do not form hydrogen bonds due to their small polarity [31].
Our calculations show that at T = 295 K the number of hydrogen bonds does not depend on pressure: at all studied pressures the number of hydrogen bonds is with-  ing 410-460 in the system of 1000 molecules, i.e. less than 0.5 bonds on a molecule. Such small number of hydrogen bonds does not allow to form a network. A common problem os H-bonds definition from geometrical criteria is their strong dependence on the choice of cut-off parameters. A usual approach to selection of r cut is based on the first minimum of the O − H partial RDF. However, there is no such a simple criterion for the angle cut-off. In our simulations we observe that the angle cut-of strongly influence the results. For instance, if we take α cut = 40 o the number of H-bonds doubles and becomes about 800-860 bonds per 1000 molecules of PG. However, independently on the choosen cut-offs we observe that the number of H-bonds demonstrates very tiny dependence on pressure.

SUPPLEMENTARY MATERIALS EXPERIMENTAL METHODS
Here we provide more detailed information on the experimental protocol employed in the work.
There were problems with determining the density experimentally through a change in the length of the sample, because at low pressures (up to 0.2 GPa) it is difficult to provide all-round uniform compression, it is uniaxial in nature. Therefore, we calculated the density by integrating the equation of state. Using the definition of the isothermal bulk modulus, where V is volume, we obtain that ∂ρ ∂P T = ρ/B T and, subsequently, for pressure dependence of density Here one can take the atmospheric value of the initial density ρ 0 in the case of our pressure range. Adiabatic bulk modulus Bs can be calculated as: where ρ is the density of a substance, v l islongitudinal ultrasonic wave velocity and v t transverse ultrasonic wave velocity. The adiabatic (isentropic) bulk modulus B S related to B T by the equation [33] B s B T = 1 + α 2 P T B s ρc P = 1 + ξ where α P is the volume thermal expansion coefficient, c P is the specific heat at constant pressure, we can calculate density from experimental ultrasonic transverse v t and longitudinal v l velocities as functions of pressure. From Eqs. 2 -4 one can obtain direct relation between density and the experimentally measured ultrasonic velocities: For ordinary (non-viscous) liquids Eq. 5 simplifies to In the first approximation, we can take ξ ≈ ξ 0 = const , where ξ 0 can be easily calculated from the atmospheric tabulated parameters.

COMPARISON WITH LITERATURE DATA
Most of experimental works deal either with mixtures of PG with other liquids (see, for instance, Refs. [34][35][36][37][38]), or with PG oligomers (see, for instance, Ref. [39]). To the best of our knowledge there is just a single work [40] where equation of state of pure PG is measured up to 101 MPa (0.101 GPa). Figure 6 shows a comparison of our data with the ones from Ref. [40]. One can see that the data are in perfect agreement: the maximum deviation of the data is about 0.3%.

COMPUTATIONAL METHODS
In molecular dynamics simulation we obtaine pressure as a function of density, but one of objectives of the present work is to calculate the baric derivative of pressure, i.e. the second derivative of the data from simulation. As a result, even small fluctuations of the computational data may lead to great errors in the baric derivative. In order to avoid these problems we approximate obtained equation of state (EoS) by 4-th order polynomial functions: where the density ρ is experssed in g/cm 3 and pressure P in GP a. The computational part of this work was carried out using computing resources of the federal collective usage center "Complex for simulation and data processing for mega-science facilities" at NRC "Kurchatov Institute", http://ckp.nrcki.ru, and supercomputers at Joint Supercomputer Center of the Russian Academy of Sciences (JSCC RAS). The experimental part of the work was supported by Russian Science Foundation (Grant 22-22-00530) and the computational part was supported by the Council of the President of the Russian Federation for State Support of Young Scientist (Grant MD-6103.2021.1.2).