Vapor Pressures and Heats of Sublimation of Crystalline β

Article pubs.acs.org/JPCB

Vapor Pressures and Heats of Sublimation of Crystalline β‑Cellobiose from Classical Molecular Dynamics Simulations with Quantum Mechanical Corrections Jakob Wohlert* Wallenberg Wood Science Center, KTH Royal Institute of Technology, SE-10044 Stockholm, Sweden ABSTRACT: In this paper, we report the calculation of the enthalpy of sublimation, ΔsubH, as a function of temperature of crystalline β-cellobiose from molecular dynamics (MD) simulations, using two popular carbohydrate force fields. Together with the entropy difference between the solid and the vapor, ΔS, evaluated at atmospheric pressure, ΔsubH gives the vapor pressure of cellobiose over the solid phase as a function of T. It is found that when quantum mechanical corrections to the enthalpy calculated from the distribution of normal modes is applied both force fields give ΔsubH close to experiments. The entropy change, ΔS, which is calculated within a harmonic approximation becomes too small, leading to vapor pressures that are too low. These findings are relevant to MD simulations of crystalline carbohydrates in general, e.g., native cellulose.

■

INTRODUCTION In recent years, there has been a renewed interest in using lingnocellulosic biomass, both as a renewable feedstock for energy production, especially in the form of bioethanol production from enzymatic degradation and subsequent fermentation,1 and, not the least, due to the great potential that natural fibers and biopolymers have for replacing plastics in a large variety of applications.2 This has led to an increased activity in the field of molecular dynamics (MD) simulations of carbohydrates in the solid crystalline phase, due to that naturally occurring polysaccharides, such as cellulose and chitin, which both are linear polymers of glucose, are synthesized in nature as crystalline microfibrils. As far as MD simulations are concerned, most carbohydrate force fields can be expected to work well for mono- and disaccharides in the gas phase or in aqueous solution, simply because these systems are frequently used in the parametrization process. Undoubtedly, a carefully developed force field should be sufficiently general to work well also for more complex carbohydrate systems, such as carbohydrate polymers, or in combination with other molecular systems, e.g., glycoproteins. However, the ability of these parameter sets to faithfully reproduce properties of the solid crystalline state may be questioned. For the specific case of predicting the crystal structure of cellulose, these doubts have already been shown to be well founded,3,4 but when it comes to their ability to reproduce physical and thermodynamic properties of the crystals, the question is still open. Experimental thermodynamical data for crystal-phase carbohydrates are scarce, but Oja and Suuberg5 have measured vapor pressures over the solid phase for a few mono- and © 2014 American Chemical Society

disaccharides, including cellobiose, which is a disaccharide of glucose (see Figure 1) with the same type of linkage (β-1,4) as is present in both cellulose and chitin. The vapor pressure is given by the difference in chemical potential of the two phases at a given temperature. This means that from a computational point of view the vapor pressure is also intimately connected to carbohydrate solubility since dissolution is commonly modeled as a sublimation step followed by a solvation step (see, e.g., Thompson et al.6). The recalcitrance to dissolution of crystalline cellulose is a fundamental issue in many practical applications involving biomass and has been the subject of several theoretical studies, employing classical MD simulations.7−9 It is not always straightforward to compare thermodynamic quantities from MD simulations to experiments due to that the reality to a large extent is quantum mechanical, while the computations are based on classical mechanics. One prominent example of this difficulty is water, for which standard classical models fail to reproduce both heat capacities and internal energy. However, Berens et al.11 have shown that it is possible to achieve good agreement between MD and experiments for liquid water by applying quantum corrections to the simulated classical values, which are based on a harmonic approximation of the normal modes. Although their method is general and applicable to almost any molecular or biomolecular system, it has been very sparsely recognized as such. Quantum corrections have since their original publication been applied Received: February 21, 2014 Revised: April 25, 2014 Published: May 2, 2014 5365

dx.doi.org/10.1021/jp501839k | J. Phys. Chem. B 2014, 118, 5365−5373

The Journal of Physical Chemistry B

Article

∂ln Pvap ∂(1/T )

=−

ΔsubH kB

(3)

The integrated form of the Clapeyron equation (eq 2), sometimes called the Clausius−Clapeyron equation, is frequently used in chemical engineering since it provides a convenient way of calculating latent heats from vapor pressure measurements. Assuming a constant heat of sublimation, ΔsubH is obtained by linear regression in a plot of ln Pvap vs 1/T. Alternatively, on the coexistence line, the chemical potentials of the two phases are identical, which means that the difference in Gibbs’ free energy between the phases is zero, ΔsubG = ΔsubH − TΔsubS = 0, where S is the entropy. Thus, the enthalpy of sublimation is given by ΔsubH = TΔsubS. We can now write ΔsubS = Strans + Srot + ΔSint, where Strans and Srot are the diffusive translational and rotational parts of the entropy, which are zero in the solid, and ΔSint represents the change in entropy associated with internal degrees of freedom: vibrations of bonds and angles and rotations around bonds. Thus, we have that ΔsubH − T (ΔSint + Srot + Strans) = 0

Figure 1. Projection of the β-cellobiose crystal structure (from ref 10) in the x/y plane. The black box represents the crystal unit cell, and green lines show the hydrogen bond network.

(4)

If one then assumes that the translational and rotational degrees of freedom of the vapor phase are accurately described by an ideal gas, Strans is given by the Sackur−Tetrode equation ⎡ ⎛k T ⎞ 5⎤ Strans = kB⎢ln⎜ B 3 ⎟ + ⎥ ⎣ ⎝ PΛ ⎠ 2⎦

to liquid methanol12 and solid argon,13 and recently Waheed and Edholm14 showed that it was possible to reproduce both the heat capacity and the heat of melting of ice using a rigid three-site water model when quantum corrections were included. It is reasonable to suspect that classical models for more complex molecules too need similar treatment to properly reproduce the thermodynamics of the corresponding experimental system. In the present work, we report MD simulations of cellobiose in both the crystal and in the vapor phase, using two different, widely used carbohydrate force fields, namely, GLYCAM06,15 and the latest version of the CHARMM carbohydrate force field.16,17 These two force fields (together with their respective predecessors) cover a large portion of computer simulations of carbohydrates during the past decade. From the simulations the quantum-corrected free energy difference between the crystal and the vapor is calculated at a range of temperatures. These differences are converted into a temperature-dependent vapor pressure, which is compared to experimental data.5

(5)

Here Λ is the thermal de Broglie wavelength, Λ = (h2/ 2πmkBT)1/2 where h is Planck’s constant, and m is the mass per particle. In the present case, P is of course the equilibrium vapor pressure, P = Pvap. For a polyatomic ideal gas, one can add the rotational part of the entropy as the entropy of a rigid rotor (see, e.g., Zhou and Gilson18) Srot = kB ln[(π 7/2e 3/2)(IAIBIC)1/2 (8kBT /h2)3/2 ]

(6)

where IX is the moment of inertia for the principal direction X. Combining eqs 4, 5, and 6 now gives ln Pvap = −

ΔsubH ΔSint S kT 5 + + rot + ln B3 + kBT kB kB 2 Λ

(7)

The derivative with respect to 1/T becomes

■

∂ln Pvap

THEORY In a P−T diagram, the Clapeyron equation relates the slope of the solid/vapor coexistence line to the enthalpy of sublimation, ΔsubH dPvap Δ H = sub (1) dT T ΔV where Pvap is the equilibrium vapor pressure, T absolute temperature, and ΔV the volume change of the phase transition. Assuming an ideal gas, ΔV ≈ Vgas and ΔsubH independent of T, this can be integrated to give Δ H ln Pvap = − sub + C kBT (2)

∂(1/T )

=−

∂ΔSint 1⎛ 1 ∂ΔsubH − ⎜ΔsubH + kB ⎝ T ∂(1/T ) ∂(1/T )

⎞ + 4kBT ⎟ ⎠

(8)

It is easy to verify that by using ΔSint = 0 and ΔsubH = 4kBT for a polyatomic ideal gas eq 8 becomes equal to −ΔsubH/kBT, i.e., the same result as for the Clapeyron equation. For a more complex molecule, however, such as a carbohydrate, there might also be contributions to the slope of ln Pvap vs 1/T from internal degrees of freedom as well, both from a temperaturedependent entropy and/or enthalpy. Finally, we note that assuming that the free energy of the solid phase is independent of pressure, the only P dependence in eq 4 comes from eq 5. Thus, for a general free energy difference between the solid and the vapor we can write

where kB is Boltzmann’s constant, and C is a constant of the integration. This has a derivative with respect to 1/T 5366

dx.doi.org/10.1021/jp501839k | J. Phys. Chem. B 2014, 118, 5365−5373

The Journal of Physical Chemistry B Pvap Δ G ΔG′ = sub − ln kBT kBT P′

Article

values for both energy and heat capacity are readily available from a MD simulation, the entropy is not. Therefore, in the present work, the entropy was estimated directly from the quantum mechanical harmonic oscillators, thus neglecting any anharmonic contributions, by integrating the expression for the entropy over D(ν)

(9)

where the prime means that the free energy is evaluated at a reference pressure P′ which may be different from the vapor pressure. This gives, since the free energy of sublimation (at the coexistence line) is zero, ΔsubG = 0 Pvap

S = kB

ΔG′ =− ln P′ kBT

(10)

∫0

⎞ βhν + − 1⎟D(ν)dν ⎜ exp(βhν) − 1 ⎠ ⎝ 2

∞ ⎛ βhν

where β = 1/kBT. Here the first two terms are the energies of the ground state and the excited states, respectively, of the harmonic oscillators, and the last term is the classical value which becomes NkBT, where N is the number of degrees of freedom, provided of course that D(ν) is normalized such that ∫ D(ν)dν = N. Contributions to this correction grow with increasing frequency. This means that δUqm can be expected to be smaller for models that apply hard constraints to atomic bonds than for fully flexible systems. In a similar manner, a correction to the classical heat capacity can be estimated from the normal mode distribution by integrating the quantum mechanical expression for the heat capacity over D(ν)

∫0

⎞ (βhν)2 exp(βhν) ⎜ − 1⎟D(ν)dν 2 ⎠ ⎝ (exp(βhν) − 1)

⎞ βhν ⎜ βhν − ln(1 − e−βhν)⎟D(ν)dν ⎝e − 1 ⎠

(13)

COMPUTATIONAL DETAILS Molecular dynamics simulations of cellobiose in both its crystal and its vapor phase were run with GROMACS 4.6.5,19 using two different carbohydrate force fields, namely, GLYCAM0615 and the latest version of the CHARMM carbohydrate force field.16,17 Simulations were run at constant NVT or NPT at eight different temperatures, 298, 328, 363, 383, 403, 428, 453, and 488 K. Lennard-Jones interactions were made to go smoothly to zero between 0.8 and 1.2 nm using a shift function, and electrostatic interactions were calculated using PME20,21 with a cutoff of 1 nm for the real space part. In the vapor-phase simulations, however, all interactions were truncated at 1.4 nm, which is enough to include all pairs in a single cellobiose molecule. The crystal phase was constructed from the X-ray structure of Chu and Jeffrey10 (Figure 1), which was replicated and translated in all three principal directions to create a computational box of roughly 3.3 × 4.0 × 3.2 nm, containing 108 cellobiose molecules. The equations of motion were integrated using a leapfrog algorithm with time step 1 fs. Temperature was maintained using stochastic velocity rescaling with a time constant of 1.0 ps.22 After 100 ps of MD at constant NPT at each temperature using anisotropic pressure scaling with a Berendsen barostat23 with 1 atm reference pressure to reach equilibrium, the systems were further equilibrated at constant NVT for 100 ps each, after which a production run of 100 ps was performed to use for calculating normal mode distributions or continued for 20 ns using a Parrinello−Rahman barostat24 instead, to use for calculating heat capacities from fluctuations. The vapor-phase simulations were performed on a single cellobiose molecule for 1 ns using stochastic dynamics to improve sampling.25 However, the simulations for calculating heat capacities were extended to 100 ns and simulated using Replica Exchange MD26 to improve conformational sampling. Here, exchanges between neighboring temperatures were attempted every 10 steps. In the vapor-phase simulations, center-of-mass motions (translations and rotations) were removed at every integration step. Simulations were performed with flexible bonds, but at 298 K simulations were also performed with all bonds constrained to their equilibrium lengths,27 for comparison. In the simulations used for normal mode distributions, velocities were stored every 3 fs, which is necessary to cover the high-frequency region of the normal mode spectrum. The normal mode distributions were calculated from the Fourier transform of the mass-weighted velocity autocorrelation functions of all atoms, as implemented in the GROMACS analysis tool g−dos.28

(11)

δC Vqm = kB

∞⎛

■

which is equivalent to eq 7 From the perspective of molecular dynamics computer simulations, the enthalpy of sublimation can be calculated from two separate simulations: one of the solid phase and one of the vapor phase. The heat of sublimation is then given by the gas solid solid enthalpy difference, ΔsubH = Ugas pot + Ukin + kBT − Upot − Ukin + PVsolid, where pot and kin refer to potential and kinetic energy, respectively, and Vsolid is the volume per molecule in the solid state. However, the PV term only contributes around 0.02 kJ mol−1 per cellobiose and is thus completely negligible. Moreover, since the kinetic energy contributes 1/2kBT per degree of freedom to the energy in a classical system, the kinetic energy terms cancel, and the enthalpy of sublimation solid reduces to ΔsubH = Ugas pot + kBT − Upot . However, as mentioned in the Introduction, a large part of the energy comes from motions that are, at least to some extent, of a quantum mechanical nature. Berens et al.11 showed that it is possible to account for quantum mechanical effects by using a correction to the classical energy based on the distribution of normal modes, D(ν), where ν is the normalmode frequency. If each normal mode is treated as a quantum mechanical harmonic oscillator, the correction to the energy reads δUqm = kBT

∫0

∞⎛

(12)

In the classical limit (high temperatures or low frequencies) the first term of the integrand goes to one, and the correction thus goes to zero. Contrary to the energy correction above, δCqm V is dominated by low-frequency contributions, meaning that it should be essentially the same for models both with and without bond constraints. Berens et al.11 also give an expression for the quantum correction to the classical entropy. However, while classical 5367

dx.doi.org/10.1021/jp501839k | J. Phys. Chem. B 2014, 118, 5365−5373

The Journal of Physical Chemistry B

■

Article

RESULTS AND DISCUSSION Molecular dynamics simulations were performed of cellobiose in its crystal and its vapor phase, at eight different temperatures using two different force fields. The cellobiose crystal structure was stable in all simulations. Crystal lattice parameters at room temperature are close to the experimental ones, for both force fields (see Table 1). This is in line with what has been Table 1. Crystal Unit Cell Parameters and Density of the Cellobiose Crystala a (nm) b (nm) c (nm) β (deg) ρ (g/cm3)

exptl

GLYCAM

CHARMM

1.097 (0.0004) 1.305 (0.0005) 0.509 (0.0003) 90.83 (0.005) 1.560 (0.0)

1.093 (0.014) 1.332 (0.015) 0.498 (0.016) 89.9 (1.1) 1.544 (0.0)

1.080 (0.015) 1.319 (0.013) 0.516 (0.014) 90.9 (0.7) 1.535 (0.1)

Figure 3. Normal mode distributions for a molecule in the gas phase (red) and in the crystal (black) at 298 K. The distributions and their respective cumulative integrals (dashed lines) are similar for both phases.

a

Absolute errors are given in parentheses. Experimental data from ref 10.

The heat capacity at constant pressure can be obtained from fluctuations in enthalpy, Cclass = ⟨δH2⟩/kBT2, from a MD P simulation at constant pressure. This value is usually higher than expected from the corresponding harmonic model since it includes anharmonic contributions to the heat capacity but also since it neglects quantum mechanical effects, which can be corrected for by adding a correction term (eq 12). The heat capacities of cellobiose in both its solid and its vapor phase, calculated from MD simulations employing the GLYCAM06 force field, are shown in Figure 4, Table 2, and

prevoiusly reported for the case of crystalline methyl α-Dglucopyranose,15 as well as a range of different monosaccharides,16 using the present force fields. Normal mode distributions were calculated as described in the previous section. The distributions of the cellobiose crystal at 298 K from the two different force fields are shown in Figure 2 and for GLYCAM06 both in its crystal and its vapor phase in Figure 3.

Figure 2. Normal mode distributions for the GLYCAM06 and CHARMM cellobiose crystals at 298 K. The density at low wavenumbers corresponds to translations and rotations of whole molecules. As one proceeds to the right in the figure, contributions from internal degrees of freedom, bond rotations, angle vibrations, and bond vibrationsin that orderstart to dominate the spectrum. The modes farthest to the right, around 3000 and 3800 cm−1, come from vibrations of bonded protons.

Figure 4. Heat capacities at constant pressure for cellobiose in its solid phase (black) and vapor phase (red), calculated from energy fluctuations with a quantum mechanical correction obtained from the normal mode distributions.

Table 3. An additional correction of 4kB, which is CP of a polyatomic ideal gas, was added to the vapor-phase value to account for the fact that all center-of-mass motions were removed in the gas-phase simulations. From Figure 4 it is clear that the heat capacities of the two phases are similar in this temperature interval. This can be understood by realizing that the largest contribution to CP comes from internal degrees of freedom, which all contribute kB in the classical limit, assuming a purely harmonic model. Quantum mechanical effects lead to a decrease in CP, which becomes approximately the same in both the crystal and the vapor. This is due to the fact that the normal mode distributions of the two phases are quantitatively very similar (see Figure 3). Deviations from a purely harmonic

Temperature Dependence of ΔsubH. As mentioned above, the Clapeyron equation rests on the assumption of a constant heat of sublimation, ΔsubH, which of course is an approximation. Since CP = (∂H/∂T)P, one way to investigate the severity of this approximation is to integrate the difference in heat capacities between the two phases ΔsubH =

∫T

T2

1

(CPvapor − CPsolid)dT + ΔsubH(T1)

(14) 5368

dx.doi.org/10.1021/jp501839k | J. Phys. Chem. B 2014, 118, 5365−5373

The Journal of Physical Chemistry B

Article

Table 2. Heat Capacity Per Molecule in J mol−1 K−1 as a Function of Absolute Temperature, For Crystalline Cellobiose Using GLYCAM06a T [K]

Cqm V

δCqm V

298 328 363 383 403 428 453 488

390.4 423.1 463.1 484.9 506.0 530.4 553.3 584.2

−731.5 −698.8 −658.8 −637.0 −615.9 −591.5 −568.6 −537.7

Cclass P 1139.2 1193.1 1186.2 1261.7 1191.6 1219.1 1204.9 1206.3

CP

(39.0) (46.2) (27.9) (27.6) (39.8) (42.8) (34.8) (28.5)

407.7 494.3 527.4 624.7 575.7 627.6 636.3 668.6

a

The heat capacity is calculated as CP = Cclass + Cqm P V . The statistical error in the classical value is estimated from block averaging. Figure 5. Experimental vapor pressures Pvap vs 1/T with linear fits from Oja and Suuberg.5 Dashed lines indicate the reported uncertainty of the fitted line. The red curve is the linear fit which has been corrected for a temperature-dependent ΔsubH, estimated from the difference in heat capacity between the two phases (see text).

Table 3. Heat Capacity Per Molecule in J mol−1 K−1 as a Function of Absolute Temperature for Cellobiose in the Gas Phase Using GLYCAM06a

a

T [K]

Cqm V

δCqm V

298 328 363 383 403 428 453 488

361.7 394.1 429.4 453.0 470.8 496.1 512.8 540.5

−735.2 −702.8 −667.6 −643.9 −626.1 −600.8 −584.1 −556.4

Cclass V 1170.4 1172.0 1191.4 1207.6 1222.4 1235.7 1241.5 1237.4

(1.3) (3.6) (3.8) (4.0) (4.1) (4.2) (3.4) (2.9)

Cclass V

CV

CP

435.2 469.2 523.8 563.7 596.3 634.9 657.4 681.0

443.5 477.5 532.1 572.0 604.6 643.2 665.7 689.3

Table 4. Classical Energies and Quantum Mechanical Correction Terms (Equation 11) in kJ mol−1 at P = 1 atm, As a Function of Absolute Temperaturea CHARMM T [K]

Usolid class

δUsolid qm

298 328 363 383 403 428 453 488

807.4 824.8 845.1 856.5 868.6 884.0 900.0 923.4

657.2 637.2 611.6 594.1 584.2 564.3 549.8 529.7

T [K]

Usolid class

δUsolid qm

Ugas class

298 328 363 383 403 428 453 488

854.7 872.0 891.6 902.5 915.9 929.9 945.7 966.4

689.7 666.7 638.5 627.3 612.2 595.1 584.0 565.3

1125.3 1144.3 1163.9 1177.1 1188.8 1203.1 1197.7 1214.3

δCqm V

The heat capacity is calculated as CP = + + 4kB. The statistical error in the classical value is estimated from block averaging.

model tend to increase CP. Perhaps a bit surprisingly, it can be concluded that also the anharmonic contributions are more or less the same in both phases, even though one might have expected the larger conformational flexibility in the vapor phase to have an impact. However, the error bars of the crystal-phase values are fairly large, and longer simulations are evidently needed for the difference in heat capacity to be analyzed in greater detail. When the difference Cvapor − Csolid is integrated numerically, P P it shows that ΔsubH is increasing slightly as the temperature goes up, but the total difference between 298 and 488 K amounts to only 1.4 kJ mol−1. The integration constant can be determined from the experimental value of the enthalpy,5 ΔsubH(T1) = ΔsubHexptl = −301.5 kJ mol−1, which is determined at T ≈ 488 K. In Figure 5 the experimental vapor pressures are plotted together with both a linear fit with slope ΔsubHexptl and the same curve corrected for the temperature dependence. As can be seen, the difference at room temperature between the curves is negligible. It thus seems fair to assume a constant heat of sublimation for this system, and the integrated form of the Clapeyron equation with a linear regression all the way down to room temperature evidently seems a good approximation. Vapor Pressure. Next, the potential energy and quantum mechanical correction terms, as well as the vibrational, rotational, and translational entropy of cellobiose, were calculated as described above, in both the solid and the vapor phase, at eight different temperatures and atmospheric pressure. The results are presented in Table 4 and Table 5. From the calculated energies and entropies, the free energy differences, ΔG′ = ΔH − TΔS per cellobiose molecule, and subsequently the vapor pressure as a function of temperature were calculated

Ugas class

δUgas qm

ΔsubHclass

ΔsubHcorr

210.8 211.2 211.6 212.4 212.3 212.0 211.5 209.5

303.6 299.8 303.3 306.7 305.1 307.4 308.0 306.1

δUgas qm

ΔsubHclass

ΔsubHcorr

724.7 705.8 681.5 667.0 654.1 636.9 628.1 612.5

273.1 275.0 275.3 277.8 276.3 276.8 255.8 252.0

308.1 314.1 318.3 317.5 318.2 318.6 299.9 299.2

1015.7 750.0 1033.3 725.8 1053.7 703.3 1065.7 688.4 1077.5 677.0 1092.4 659.7 1107.7 646.3 1129.0 626.1 GLYCAM06

Statistical errors are estimated to approximately 1 kJ mol−1 for Uclass and 3 kJ mol−1 for the QM correction terms. This gives a statistical error of 5 kJ mol−1 in ΔsubHcorr. a

from eq 10. The enthalpy difference can be directly identified as the enthalpy of sublimation, ΔH = ΔsubH. The results without quantum mechanical corrections are shown in Figure 6 together with experimental data from ref 5. Here it can be seen that the CHARMM vapor pressure agrees well with experimental results at high temperatures but deviates from the experimental curve at low T. The GLYCAM06 vapor pressure on the other hand is too low for all temperatures. For both force fields, ln Pvap vs 1/T is linear over the whole temperature range. The slope of that line should, according to the Clapeyron equation (eq 2) correspond to the heat of sublimation, ΔsubH. 5369

dx.doi.org/10.1021/jp501839k | J. Phys. Chem. B 2014, 118, 5365−5373

The Journal of Physical Chemistry B

Article

Table 5. Entropies in kJ mol−1 K−1 at P = 1 atm, As a Function of Absolute Temperaturea

Table 6. Average Enthalpy of Sublimations and the Slopes −d ln P/d(1/T) × kB, from Figure 6 and Figure 7, in kJ mol−1a

CHARMM Sgas trans

ΔS

0.362 0.146 0.401 0.147 0.444 0.149 0.470 0.149 0.492 0.150 0.525 0.150 0.554 0.151 0.595 0.152 GLYCAM06

0.182 0.184 0.186 0.187 0.188 0.189 0.190 0.192

0.301 0.304 0.303 0.303 0.300 0.297 0.298 0.289

Ssolid harm

Sgas harm

Sgas rot

Sgas trans

ΔS

0.379 0.418 0.469 0.498 0.523 0.556 0.580 0.633

0.358 0.398 0.443 0.475 0.494 0.527 0.531 0.562

0.147 0.147 0.149 0.149 0.150 0.151 0.152 0.152

0.182 0.184 0.186 0.187 0.188 0.189 0.190 0.192

0.308 0.311 0.309 0.313 0.309 0.311 0.293 0.273

T [K]

Ssolid harm

Sgas harm

298 328 363 383 403 428 453 488

0.389 0.428 0.476 0.503 0.530 0.567 0.597 0.650

T [K] 298 328 363 383 403 428 453 488

Sgas rot

without QM corr. CHARMM GLYCAM06 exptl

with QM corr.

slope

⟨ΔsubHclass⟩

slope

⟨ΔsubHcorr⟩

203.1 279.9 301.5 ± 43.9

211.4 (0.9) 270.3 (10.2)

284.6 300.4

305.0 (2.7) 311.7 (8.3)

a

Standard deviations are given in parentheses. Ideally, the slope should be equal to −ΔsubH, but this is not the case since neither ΔsubH nor ΔS is really constant. Experimental value (calculated in the range 474− 488 K) from ref 5.

the thermal energy, kBT. At the temperatures investigated here, this means that quantum mechanical effects start to become important for wavenumbers around 100 cm−1 and larger. Thus, the normal-mode frequencies present in this system reside to a large extent in the quantum mechanical regime. For that reason, it makes sense to investigate quantum mechanical corrections to the energy. Calculated correction terms are presented in Table 4, and the vapor pressure as a function of temperature including these corrections is shown in Figure 7.

a

Harmonic contributions are calculated from eq 13. The translational entropy of the gas phase comes from the Sackur−Tetrode equation (eq 5) and the rotational part from a rigid rotor approximation (eq 6). The statistical error in ΔS is estimated from block averaging to approximately 0.003 kJ mol−1 K−1, which corresponds to an error of between 0.9 and 1.5 kJ mol−1 in ΔG″.

Figure 7. Calculated and experimental vapor pressures Pvap vs 1/T with linear fits. Experimental data with fitted linear function from Oja and Suuberg.5 Dashed lines indicate the reported uncertainty of the fitted line.

The effect of using quantum mechanical corrections can be seen in Figure 7. Now, the resulting slopes of ln Pvap vs 1/T from the two different force fields are very similar, and furthermore, since ΔS too is approximately the same in both cases (see Table 5), also the calculated vapor pressures, in absolute terms, become the same with both force fields. Thus, despite giving totally different results initially, they give the same answer when correction terms are applied. However, both force fields now give vapor pressures that are much too low compared to the experimental data. Evidently there is some fairly large term, ∼20kBT, that is missing in ΔG′. By comparing eq 2 to eq 10, one realizes that the constant in the integrated form of Clapeyron’s equation can be identified as ΔS. This is discussed in more detail below. As can be seen in Table 4, the enthalpy of sublimation is approximately constant over the range of temperatures investigated also after the inclusion of a quantum mechanical correction. The values have risen, and the average values are now −305.0 and −311.7 kJ mol−1 for CHARMM and

Figure 6. Experimental and calculated vapor pressures, Pvap vs 1/T, without QM corrections. The experimental data are fitted with the linear function from Oja and Suuberg.5 Dashed lines indicate the reported uncertainty of the fit.

Indeed, the slopes of the fitted lines in Figure 6 correspond well to the averages of the calculated ΔsubHclass, which are approximately constant over the whole temperature range (see Table 6). From the slope of the experimental data,5 the experimental ΔsubH has been estimated to 301.5 ± 43.9 kJ mol−1. This means that the slopes from simulations, which are −203.1 and −279.9 kJ mol−1 for CHARMM and GLYCAM06, respectively, are significantly smaller than the experimental value, even though the result using GLYCAM06 falls within the given experimental error. Quantum Mechanical Corrections to ΔsubH. Quantum mechanical effects are expected to be important when the quantum mechanical energy hν is comparable to or larger than 5370

dx.doi.org/10.1021/jp501839k | J. Phys. Chem. B 2014, 118, 5365−5373

The Journal of Physical Chemistry B

Article

rotations of hydroxyl and hydroxymethyl groups, as well as rotations around the glycosidic bond, which are far from harmonic. Second, the entropy of the vapor phase is based on fairly short (1 ns) sampling times. Some motions within a complex molecule naturally have longer correlation times than others. Again, for the case of a disaccharide, rotation around the glycosidic bond is an example of this. It has been shown that to reach convergence of the configurational entropy for a disaccharide in water even on a 50 ns time scale some method to enhance the sampling of configurational space is needed,29 largely due to slow rotations around the glycosidic linkage. Unlike the case for the heat capacities above, Replica Exchange MD cannot be used here to calculate the entropy since it does not produce continuous trajectories. There are other methods to calculate the entropy that potentially could resolve these issues, for instance, the quasiharmonic approximation by Karplus and Kushick,30 which can be applied on much longer time scales than the method used here and also does not require a continuous trajectory. On the other hand, it scales very badly with system size and would therefore be computationally inefficient for the crystal-phase calculations. One thus has to bear in mind that there is an unknown contribution to ΔS that comes from both a (potentially) insufficient sampling of the conformational space, which is positive, and deviations from a perfectly harmonic model, which is negative.31 However, the fairly large error in ΔS of more than 50% of the calculated values makes it seem quite unlikely that this could be remedied by more extensive sampling and/or anharmonic corrections alone. Outlook. There is no doubt that there is great potential for improvement of these force fields. Considering that they were all designed for mono- and disaccharides in aqueous solution, it is almost surprising that they work so well also for carbohydrates in the solid phase. It has been noted before that the force fields used in the present study both have problems reproducing the experimental crystal structure of cellulose,3,4,32 which has been determined by X-ray and neutron diffraction.33 However, Hadden et al.4 showed that assigning new partial charges to their GLYCAM06 cellulose model, based on QM calculations on larger cellulose fragments (cellotriose) in an environment that mimicked a cellulose crystal, resulted in a structure that was considerably more stable. The stability arose from higher partial charges especially on the OH groups, leading to stronger hydrogen bonds. The same effect was also achieved by using a force field that explicitly models oxygen lone pairs.4 It is quite obvious that a reparametrized charge distribution would also affect the thermodynamics of the crystal. It is however difficult to speculate about what the effect would be in terms of equilibrium vapor pressures and the analysis presented in this work. Stronger hydrogen bonds would obviously lower the energy of the crystal phase but at the same time also lower the entropy due to that fluctuations would be further restricted. Furthermore, the effect on the thermodynamics of the vapor phase, or indeed the aqueous phase, is also something that needs to be considered. It may very well be that a different parameter set is warranted for each different phase. An alternative could perhaps be to use a polarizable force field, which is sensitive to changes in the electrostatic environment. On a final note, it is important to remember that the discussion presented here is based on comparison with only

GLYCAM06, respectively, which in both cases are in much better agreement with the experimental result. The slopes of the linear fits are fairly close to the average enthalpy, but there is a difference of about 7% in the case of CHARMM. The reason why there is not perfect agreement between the calculated enthalpy and the slopes of ln Pvap vs 1/T is that variations of both enthalpy and entropy with temperature will give contributions to the total derivative (see eq 8), but the discrepancies are sufficiently small for the assumption of a constant ΔsubH to remain valid. The analysis shows that the largest contribution to δUqm in all cases comes from the ground state energy (the first term of the integrand in eq 11), and the dominating contribution to this term is the high-frequency part of the normal mode distributions, which mainly comes from vibrations of covalent bonds. This is also where the largest difference between the solid and the gas phase appears. This also has the effect that the correction terms are quite insensitive to simulation length since the fast motions converge rapidly. In contrast, simulations that are performed with all bonds constrained to their respective equilibrium value, which is standard procedure in MD simulations, lack the high-frequency modes. In this case, the quantum correction to ΔsubH is around −10 kJ mol−1, which is fairly small, and also negative, which means that the correction only makes the situation worse compared to the experimental value. As one example, for GLYCAM06 at 298 K with bond constraints, the total enthalpy of sublimation becomes ΔsubHcorr = 248.1 kJ mol−1. It is important to realize that the level of agreement between simulations and experiements presented here is not a coincidence but rather a result of successful force field parametrizations. Both GLYCAM06 and CHARMM are to some extent validated against experimental vibrational frequencies of single bonded terms,15,16 and in the case of GLYCAM06, this has been shown to lead to good agreement also for the more complex normal modes present in a carbohydrate crystal,15 which is a fact that gives additional credability to the present calculations. Entropy of Sublimation. The change in entropy between the solid and the vapor phase at atmospheric pressure is calculated to around 0.3 kJ mol−1 K−1 per molecule (see Table 5). It is approximately constant in temperature and between force fields. It is notable that the largest part of the difference comes from the translational and rotational entropy of the vapor phase. The vibrational entropy change is small, and moreover, it is negative. It might seem strange that the vibrational energy decreases when the molecule goes from the solid to the vapor, but one must remember that in the solid also vibrations around its center of mass are included in the vibrational entropy, whereas in the vapor, these degrees of freedom are accounted for explicitly. In that respect, these numbers are not directly comparable. As mentioned above, the calculated vapor pressures end up being too low for both force fields, with a difference to the experimental data corresponding to an error in ΔG′ of around 20kBT, indepentent of T. This corresponds to an error in ΔS of 0.17 kJ mol−1 K−1 per molecule, over 50% of the calculated values. Therefore, it is reasonable to ask whether the entropy change of the phase transition can be accurately calculated. The entropy from eq 13 suffers from two potential drawbacks. First, it relies on a harmonic approximation, although the normal modes are certainly anharmonic to some extent. Specifically, for a disaccharide, there are modes that are associated with 5371

dx.doi.org/10.1021/jp501839k | J. Phys. Chem. B 2014, 118, 5365−5373

The Journal of Physical Chemistry B

■

one experimental study, which to our knowledge is the only one available. Obviously more experimental thermodynamical data for carbohydrates in the solid phase would be useful to further improve the performance of classical force fields for simulations of carbohydrates.

Article

AUTHOR INFORMATION

Corresponding Author

*E-mail: [email protected]. Notes

■

The authors declare no competing financial interest.

■

CONCLUSIONS In the present paper the equilibrium vapor pressure over the solid phase of crystalline cellobiose is calculated from molecular dynamics simulations using two different classical carbohydrate force fields, GLYCAM06 and CHARMM, at eight different temperatures from 298 to 488 K. The vapor pressures are calculated from the difference in chemical potential between the two phases, which in turn is calculated from the differences in enthalpy and entropy separately. The calculated vapor pressures are compared to experimental data5 acquired at high temperatures, which is extrapolated to lower temperatures using Clapeyron’s equation. By integrating the difference in heat capacity obtained from simulations, it is shown that the temperature dependence of the heat of sublimation, ΔsubH, is small, only 1.4 kJ mol−1, i.e., approximately 0.5%, between 488 K and room temperature. This suggests that it is indeed a reasonable approximation to assume constant ΔsubH and use linear regression of ln Pvap vs 1/ T to extrapolate the vapor pressure to lower temperatures. This is also corroborated by the potential energies which are obtained from the simulations. The difference in enthalpy between the two phases that comes directly from the simulations is approximately constant over the whole temperature range, but the value differs substantially from the experimental ΔsubH for both force fields. However, when these values are supplemented with quantum mechanical corrections obtained from the distribution of normal modes, they are significantly improved, ending up close to the experimental value of 301.5 kJ mol−1. The quantum mechanical correction is dominated by the zero-point energy of the highfrequency motions, i.e., vibrations of covalent bonds. Entropies too are estimated from the normal mode densities. The difference between the solid and the vapor phase is approximately constant, around 0.3 kJ mol−1 K−1, and is almost entirely given by the translational and rotational entropy of the (ideal) vapor. This value seems to be too low, however, since the calculated vapor pressures are smaller than the experimental value, at all temperatures. Insufficient sampling of the conformational space of the disaccharides in the vapor phase, in combination with the shortcomings of using a purely harmonic model, are identified as possible reasons for this. Poor transferability of the static partial charges between different phases is another. The results presented here show that it is not a trivial task to use classical molecular dynamics simulations to study the phase behavior of carbohydrates, or any complex molecule for that matter, since quantum mechanical effects have to be accounted for to faithfully reproduce the latent heat of the solid/vapor phase transition. When corrections are applied, classical carbohydrate force fields provide good agreement with experimental data. This is important to remember when classical simulations are used to study processes in which carbohydrates are transferred between different phases, as in, e.g., dissolution. To make accurate predictions of the solubility or the stability of the solid phase, the free energy may have to be quantum corrected.

ACKNOWLEDGMENTS The author wishes to express his gratitude to Olle Edholm for helpful discussions and comments. Comments from Malin Wohlert and Lars Berglund are also gratefully acknowledged. Computational resources were provided by the PDC Center for High Performance Computing. Molecular graphics were produced using VMD.34

■

REFERENCES

(1) Himmel, M. E.; Ding, S.-Y.; Johnson, D. K.; Adney, W. S.; Nimlos, M. R.; Brady, J. W.; Foust, T. D. Biomass Recalcitrance: Engineering Plants and Enzymes for Biofuels Production. Science 2007, 315, 804−807. (2) Moon, R. J.; Martini, A.; Nairn, J.; Simonsen, J.; Youngblood, J. Cellulose Nanomaterials Review: Structure, Properties and Nanocomposites. Chem. Soc. Rev. 2011, 40, 3941−3994. (3) Matthews, J. F.; Beckham, G. T.; Bergenstråhle-Wohlert, M.; Brady, J. W.; Himmel, M. E.; Crowley, M. F. Comparison of Cellulose I Beta Simulations With Three Carbohydrate Force Fields. J. Chem. Theory Comput. 2012, 8, 735−748. (4) Hadden, J. A.; French, A. D.; Woods, R. J. Unraveling Cellulose Microfibrils: A Twisted Tale. Biopolymers 2013, 99, 746−756. (5) Oja, V.; Suuberg, E. M. Vapor Pressures and Enthalpies of Sublimation of D-glucose, D-xylose, Cellobiose, and Levoglucosan. J. Chem. Eng. Data 1999, 44, 26−29. (6) Thompson, J. D.; Cramer, C. J.; Truhlar, D. G. Predicting Aqueous Solubilities from Aqueous Free Energies of Solvation and Experimental or Calculated Vapor Pressures of Pure Substances. J. Chem. Phys. 2003, 119, 1661−1670. (7) Bergenstråhle, M.; Wohlert, J.; Himmel, M. E.; Brady, J. W. Simulation Studies of the Insolubility of Cellulose. Carbohydr. Res. 2010, 345, 2060−2066. (8) Payne, C. M.; Himmel, M. E.; Crowley, M. F.; Beckham, G. T. Decrystallization of Oligosaccharides from the Cellulose I Beta Surface With Molecular Simulation. J. Phys. Chem. Lett. 2011, 2, 1546−1550. (9) Gross, A. S.; Bell, A. T.; Chu, J.-W. Thermodynamics of Cellulose Solvation in Water and the Ionic Liquid 1-butyl-3-methylimidazolium Chloride. J. Phys. Chem. B 2011, 115, 13433−13440. (10) Chu, S. S. C.; Jeffrey, G. A. The refinement of the Crystal Structures of β-D-Glucose and Cellobiose. Acta Crystallogr. 1968, B24, 830−838. (11) Berens, P. H.; Mackay, D. H. J.; White, G. M.; Wilson, K. R. Thermodynamics and Quantum Corrections from Molecular Dynamics for Liquid Water. J. Chem. Phys. 1983, 79, 2375−2389. (12) Hawlicka, E.; Pálinkás, G.; Heizinger, K. A Molecular Dynamics Study of Liquid Methanol with a Flexible Six-site Model. Chem. Phys. Lett. 1989, 154, 255−259. (13) Hardy, R. J.; Lacks, D. J.; Shukla, R. C. Quantum Corrections to the Simulated Properties of Solids. Phys. Rev. B 1998, 57, 833−838. (14) Waheed, Q.; Edholm, O. Quantum Corrections to Classical Molecular Dynamics Simulations of Water and Ice. J. Chem. Theory Comput. 2011, 7, 2903−2909. (15) Kirschner, K. N.; Yongye, A. B.; Tschampel, S. M.; GonzálezOuteiriño, J.; Daniels, C. R.; Foley, B. L.; Woods, R. J. GLYCAM06: A Generalizable Biomolecular Force Field. Carbohydrates. J. Comput. Chem. 2007, 29, 622−655. (16) Guvench, O.; Greene, S. N.; Kamath, G.; Brady, J. W.; Venable, R. M.; Pastor, R. W.; MacKerell, A. D., Jr. Additive Empirical Force Field for Hexopyranose Monosaccharides. J. Comput. Chem. 2008, 29, 2543−2564. 5372

dx.doi.org/10.1021/jp501839k | J. Phys. Chem. B 2014, 118, 5365−5373

The Journal of Physical Chemistry B

Article

(17) Guvench, O.; Hatcher, E.; Venable, R. M.; Pastor, R. W.; MacKerell, A. D., Jr. CHARMM Aditive All-atom Force Field for Glycosidic Linkages Between Hexopyranoses. J. Chem. Theory Comput. 2009, 5, 2353−2370. (18) Zhou, H.-X.; Gilson, M. K. Theory of Free Energy and Entropy in Noncovalent Binding. Chem. Rev. 2009, 109, 4092−4107. (19) Hess, B.; Kutzner, C.; van der Spoel, D.; Lindahl, E. GROMACS 4: Algorithms for Highly Efficient, Load-balanced, and Scalable Molecular Simulation. J. Chem. Theory Comput. 2008, 4, 435−447. (20) Darden, T.; York, D.; Pedersen, L. Particle Mesh Ewald: An Nlog(N) Method for Ewald Sums in Large Systems. J. Chem. Phys. 1993, 98, 10089−10092. (21) Essmann, U.; Perera, L.; Berkowitz, M. L.; Darden, T.; Lee, H.; Pedersen, L. G. A Smooth Particle Mesh Ewald Method. J. Chem. Phys. 1995, 103, 8577−8592. (22) G. Bussi, D. D.; Parrinello, M. Canonical Sampling Through Velocity Rescaling. J. Chem. Phys. 2007, 126, 014101. (23) Berendsen, H. J. C.; adn A. DiNola, J. P. M. P.; Haak, J. R. Molecular Dynamics with Coupling to an External Bath. J. Chem. Phys. 1984, 81, 3684−3690. (24) Parrinello, M.; Rahman, A. Polymorphic Transitions in Single Crystals: a New Molecular Dynamics Method. J. Appl. Phys. 1981, 52, 7182−7190. (25) van Gunsteren, W. F.; Berendsen, H. J. C. A Leap-frog Algorithm for Stochastic Dynamics. Mol. Simul. 1988, 1, 173−185. (26) Sugita, Y.; Okamoto, Y. Replica-exchange Molecular Dynamics Method for Protein Folding. Chem. Phys. Lett. 1999, 314, 141−151. (27) Hess, B. P-LINCS: A Parallel Linear Constraint Solver for Molecular Simulation. J. Chem. Theory Comput. 2008, 4, 116−122. (28) Caleman, C.; van Maaren, P. J.; Hong, M.; Hub, J. S.; Costa, L. T.; van der Spoel, D. Force Field Benchmark of Organic Liquids: Density, Enthalpy of Vaporization, Heat Capacities, Surface Tension, Isothermal Compressibility, Volumetric Expansion Coefficient, and Dielectric Constant. J. Chem. Theory Comput. 2012, 8, 61−74. (29) Perić-Hassler, L.; Hansen, H. S.; Baron, R.; Hünenberger, P. H. Conformational Properties of Glucose-based Disaccharides Investigated Using Molecular Dynamics Simulations with Local Elevation Umbrella Sampling. Carbohydr. Res. 2010, 345, 1781−1801. (30) Karplus, M.; Kushick, J. N. Method for Estimating the Configurational Entropy of Macromolecules. Macromolecules 1981, 14, 325−332. (31) Baron, R.; Hünenberger, P. H.; McCammon, J. A. Absolute Single-molecule Entropies from Quasi-harmonic Analysis of Microsecond Molecular Dynamics: Correction Terms and Convergence Properties. J. Chem. Theory Comput. 2009, 5, 3150−3160. (32) Paavilainen, S.; Rog, T.; Vattulainen, I. Analysis of Twisting of Cellulose Nanofibrils in Atomistic Molecular Dynamics Simulations. J. Phys. Chem. B 2011, 115, 3747−3755. (33) Nishiyama, Y.; Langan, P.; Chanzy, H. Crystal Structure and Hydrogen-Bonding System in Cellulose Iβ from Synchrotron X-ray and Neutron Fiber Diffraction. J. Am. Chem. Soc. 2002, 124, 9074− 9082. (34) Humphrey, W.; Dalke, A.; Schulten, K. VMD - Visual Molecular Dynamics. J. Mol. Graphics 1996, 14, 33−38.

5373

dx.doi.org/10.1021/jp501839k | J. Phys. Chem. B 2014, 118, 5365−5373