Article pubs.acs.org/JPCB
Effects of Temperature on the Properties of Glycerol: A Computer Simulation Study of Five Different Force Fields David A. Jahn, Frederick O. Akinkunmi, and Nicolas Giovambattista* Department of Physics, Brooklyn College of the City University of New York, 2900 Bedford Avenue, Brooklyn, New York 11210, United States S Supporting Information *
ABSTRACT: We perform molecular dynamics simulations of glycerol (propane-1,2,3-triol) at normal pressure and a wide range of temperatures (300−460 K) and study the sensitivity of simulation results to the force field (FF) considered. We employ five commonly used FFs: (i) AMBER, (ii) CHARMM22, and (iii) three versions of the OPLS-AA FF (OPLS1, OPLS2, and OPLS3). We study thermodynamic (density ρ(T), thermal expansion coefficient αP(T), isobaric specific heat cP(T)), dynamic (diffusion coefficient D(T)), as well as structural properties (molecular conformations and hydrogen-bond statistics). In comparison with experiments, FFs i and iii provide reasonable estimations of ρ(T) with deviations of ≤4.5%; for FF ii, deviations in density are more pronounced, ≤9%. Values of αP(T) vary considerably among the FFs; e.g., deviations are ≤9% for OPLS1-FF and ≤60% for FF ii. For all models studied, values of cP(T) are approximately twice the corresponding experimental values. Diffusion coefficients are very sensitive to the FFs considered. Specifically, for FFs i and ii and OPLS3, the values of D(T) are remarkably close to the experimental values over the whole range of temperatures studied. Instead, in the cases of OPLS1 and OPLS2-FFs, D(T) is underestimated by approximately 2 orders of magnitude. Interestingly, in all cases, D(T) can be well described by a Vogel−Tamman−Fulcher equation, as observed in experiments. We present a detailed characterization of glycerol backbone conformation based on the traditional classification introduced by Bastiansen, defined in terms of glycerol’s OCCC dihedral angles. All FFs indicate that the conformer population varies smoothly with temperature. However, the FFs provide very different conformer distributions. This implies that, from the microscopic point of view, these glycerol models may provide very different liquid environments for, for example, guest biomolecules and hence may play a relevant role in interpreting simulation results involving glycerol-based solutions. We also discuss the statistics of inter- and intramolecular hydrogen bonds (HBs). The FFs are qualitatively comparable regarding HB statistics; however, quantitative differences remain. For example, molecules form a total of 5.5−7 HBs at T = 350 K, depending on the FF considered, including at least one intramolecular HB. Glycerol is also a unique cryoprotectant3 and antifreeze, chosen by many organisms as a tool to survive subzero temperatures.4 In addition, glycerol is a sustainable solvent for green chemistry, with useful properties for catalysis, organic synthesis, separations, and materials chemistry.5 We note that glycerol is a byproduct of the biodiesel industry, and hence, an increasing demand of glycerol is expected to benefit this industry, reducing the price of biodiesels.1 The molecular structure of glycerol is relatively simple, and yet glycerol exhibits intriguing thermodynamic properties, such as high surface tension γ = 63.4 dyn/cm 6 at T = 293 K, close to the surface tension of water (γ ≈ 72.5 dyn/cm) and much higher than the surface tension of most organic molecules.7,8 Glycerol’s melting temperature is relatively low, Tm = 291 K, and its boiling temperature is relatively high, Tm = 573 K.
1. INTRODUCTION Glycerol (propane-1,2,3-triol, Figure 1) is a sugar alcohol that plays a very important role in many technological and scientific applications. It is used intensively in the pharmaceutical, cosmetics, and food industries, as well as in the manufacturing of inks, adhesives, synthetic plastics, and even explosives;1 a list of more that 1500 applications of glycerol can be found in ref 2.
Received: June 13, 2014 Revised: August 22, 2014 Published: September 4, 2014
Figure 1. Glycerol molecule with atom labels used in this work. © 2014 American Chemical Society
11284
dx.doi.org/10.1021/jp5059098 | J. Phys. Chem. B 2014, 118, 11284−11294
The Journal of Physical Chemistry B
Article
Interestingly, glycerol does not crystallize easily;9,10 it can be supercooled by as much as 100 K and has a glass transition temperature of Tg = 190 K.11 The glass forming ability of glycerol makes it an ideal substance for the study of the glass and supercooled state, as well as the phenomenology associated with the glass transition (see refs12−16.). At the molecular level, glycerol exhibits a rather complex behavior. Having three OH groups, glycerol can form intra- as well as intermolecular hydrogen bonds.17 This allows molecules in the liquid and glass states to form a complex hydrogen bonded network,18 with molecules adopting a wide range of conformations. Glycerol has a total of 126 conformations which are defined by the relative location of the O, C, and H atoms.19 If only the heavy atoms (O and C) are considered, glycerol structure is characterized by six different backbone conformations. Not surprisingly, determining the structure of glycerol has been rather challenging both experimentally and in computer simulations. Disagreement in the conformation distribution of noncrystalline glycerol exists between experiments and computer simulations (see refs 20 and 21), among independent simulations (see refs17, 22, and 23), as well as among experimental works (see refs 24−27). Interestingly, computational studies of glycerol are not abundant. This is in contrast with the large number of experimental studies, which include thermodynamic techniques,28 infrared spectroscopy,17 NMR,29 Raman spectroscopy,29,30 and neutron scattering experiments.18,25,31,32 The first computer simulation study of glycerol was performed by Root and Stillinger approximately 25 years ago33 in an attempt to characterize the structure of the liquid and amorphous solid state. Because of the limited computer simulation capabilities at the time, this work employed a “united atoms” approximation wherein each of the three C atoms of glycerol and their directly attached H atoms are treated as single interacting sites. United atom models of glycerol were used in many subsequent studies (see refs 34 and 35), and the first computer simulation study of glycerol using a full-atomistic model was performed in 1999 by Chelli and collaborators.23 In a series of publications, Chelli et al. studied different properties of glycerol in the crystal, liquid, and glass states, including structural (e.g., radial distribution functions and conformation distributions), thermodynamic (e.g., density and specific heat), and dynamical (e.g., mean square displacements) properties.20,23,36 In the past few years, a few full-atom computer simulations of glycerol and its aqueous solutions have been performed (see refs 21, 22, and 37−41) at limited thermodynamic conditions (pressure, temperature, and concentration) and focused on different properties. In addition, most of these works employ different force fields (FFs). It follows that it is difficult, if not impossible, to compare results obtained from different studies, and it is unclear how successful many of the available FFs are in reproducing available experimental data. On account of its molecular flexibility and ability to form hydrogen bonds (HBs), it is very important to compare the conformations adopted by glycerol using different FFs. Different glycerol conformations may result in different microscopic environments for guest molecules, such as biomolecules in glycerol-based solutions. If so, computational studies of glycerol and its solutions may be dependent, even at the qualitative level, on the FF employed. In this work, we perform molecular dynamics simulations of glycerol at P = 1 atm and temperatures in the range 300−460 K, using five different FFs. Our primary goal is to compare the performance of the FFs with available experimental data. Our
study is focused on thermodynamic (density, thermal expansion coefficient, specific heat), dynamical (diffusion coefficient), and structural (conformation distribution and HB statistics) properties. This work is organized as follows. In section 2, we discuss the computational details and describe the FFs studied. The thermodynamic and dynamical properties of glycerol are presented in sections 3 and 4, respectively. Section 5 includes a discussion of the structural properties of glycerol, including glycerol backbone conformations and HB statistics. A summary is presented in section 6.
2. COMPUTER SIMULATION DETAILS We perform molecular dynamics simulations (MD) of glycerol at P = 0.1 MPa and T = 300, 320, 350, 380, 400, 420, 440, and 460 K. Five different glycerol FFs are considered: CHARMM22, AMBER, and three versions of the OPLS-AA force field (described below). In all cases, the system is composed of N = 700 molecules placed in a cubic box with periodic boundary conditions. The box side length L varies with temperature and the FF considered, but it is found to be in the range of ∼40−55 Å. All simulations are performed using the LAMMPS software package.42 The temperature and pressure are controlled using a Nosé−Hoover style thermostat and barostat employing the equations of motion described by Shinoda et al.43 The thermostat uses a time constant for temperature coupling of 0.1 ps,44 and the barostat uses a pressure relaxation time constant of 5 ps. Long-range electrostatic interactions are calculated using a particle−particle particle-mesh solver45 with an accuracy value of 10−4. Lennard-Jones and Coulombic interactions are brought smoothly to zero from r1 = 10 Å to r2 = 11 Å using a polynomial switching function. The starting glycerol configuration for a given FF is prepared as follows. First, an independent MD simulation of a single glycerol molecule in a water bath is performed at T = 300 K, using the SPC/E model46 for water, from which 700 glycerol configurations (atom positions/velocities) are obtained. These configurations are used to locate randomly 700 molecules in a cubic box of L = 80 Å, ensuring that no atom of a given molecule is within a distance 3.5 Å from the atom of any other molecule in the system. The box is then compressed rapidly at 300 K until the corresponding experimental liquid density at normal pressure (∼1.26 g/cm3 47) is reached. The resulting glycerol configuration is then used as the starting configuration for simulation runs at constant temperature T and pressure P = 0.1 MPa. Simulations are run for up to Δtsim = 200 ns; see Table 1. For all temperatures studied, Δtsim is chosen so that Table 1. Total Simulation Times for the Different FFs and Temperatures Studieda
a
11285
T (K)
R-FF/BC-FF/OPLS3-FF
OPLS1-FF/OPLS2-FF
460 440 420 400 380 350 320 300
5 (3) 5 (3) 5 (3) 5 (3) 5 (3) 15 (8) 40 (15) 200 (80)
5 (3) 5 (3) 5 (3) 15 (8) 35 (20) 120 (60) 200 (80)
Equilibration times are given in parentheses. Times are given in ns. dx.doi.org/10.1021/jp5059098 | J. Phys. Chem. B 2014, 118, 11284−11294
The Journal of Physical Chemistry B
Article
molecules diffuse for at least 2−3 times the size of the glycerol molecule (approximately 10 Å), allowing the system to relax to equilibrium. This allows us to study temperatures as low as T = 300 or 320 K, depending on the FF considered. In all cases, the simulation time step is 1.0 fs. We confirm that our results are independent of the starting configuration by replicating the described procedure for every FF at T = 400 K, varying the initial location of the molecules of the system. Additionally, for two different models (Blieck− Chelli and Reiling FF) we perform simulations using glycerol molecules extracted from vapor phase simulations, all having identical conformations. In all cases, after thermalization we obtain the same results (density, diffusivity, conformations, etc.), indicating that the system has reached equilibrium and our results are independent of the starting configuration. Glycerol Force Fields. Five all-atom (AA) glycerol force fields are considered: the model developed by Reiling et al. (RFF)48 based on the CHARMM22 FF,49,50 the model introduced by Chelli et al.23 and modified by Blieck et al.37,51 (BC-FF) based on the AMBER FF,52 and three models based on the OPLS-AA FF.53 In the last case, we consider the original glycerol parametrization of Jorgensen et al.53,54 (OPLS1-FF), a second glycerol model (OPLS2-FF) based on OPLS1-FF where the values assigned to the CCCO dihedral interaction parameters are modified,55 and a third model used in ref 38 (OPLS3-FF) where the CCCO, OCCO, and HOCH dihedral interaction parameters are modified from the original parametrization.56 The parameters that define these glycerol models are provided in the Supporting Information. We choose to study three different versions of OPLS-AA FF because we found it to be somewhat ambiguous to identify glycerol parameters from the OPLS-AA parameter tables. In addition, independent computational studies of glycerol, all using the OPLS-AA FF, may differ in the parameters employed (see refs 38, 41, and 56−58). We stress that there are models of glycerol available other than those studied in this work (see refs 22, 38, and 41). The glycerol models considered in this work are based on FFs commonly used in the literature to study molecular systems. The target FFs all employ equivalent functional forms for the intramolecular potential energies associated with bond stretching, angular, and torsional interactions, as well as employ Lennard-Jones and Coulombic interactions. They differ in the potential parameters used as well as the partial charges assigned to each glycerol atom. The parameters and partial charges used in this study are given in the Supporting Information. Additionally, these FFs include different parameters to weight the so-called 1−4 pair interactions (see ref 59). Another important difference among the FFs is the rule used to combine the Lennard-Jones parameters between different atoms. Specifically, in the cases of the R-FF and BC-FF, the Lennard-Jones parameters are obtained using the Lorentz− Berthelot combination rule, ϵi,j = (ϵiϵj)1/2 and σi,j = 0.5(σi + σj), where i and j refer to the interacting atoms. In the OPLS-AA FFs, these parameters are defined as ϵi,j = (ϵiϵj)1/2 and σi,j = (σiσj)1/2.59
Figure 2. (a) Density and (b) thermal expansion coefficient of glycerol as a function of temperature at P = 0.1 MPa for the different FFs studied. Experimental data in (a) is from ref 60. Error bars in (a) are smaller than the corresponding symbol size and result in relative errors in (b) being smaller than 10%.
the experimental data, with differences of approximately δρ < 0.05 g/cm3 (i.e.,
