Rapid determination of entropy for flexible molecules in condensed phase from the two-phase thermodynamic model

Abstract

The conformation of a flexible molecule varies mainly via changes in the angle of dihedral torsions. Here we show that the entropy and free energy associated with such torsional rotations can be efficiently and accurately determined from separating the torsional density of states (DoS) into a solid-like (harmonic vibration) and a gas-like (free rotation) component, and applying suitable quantum statistics to the corresponding components of the DoS. We examined the calculation of internal rotation entropy of ethane and methanol molecules in liquid and vapor phases over a wide range of temperatures. Furthermore, the torsional entropy of larger, flexible molecules, such as butane and hexane, are studied at the standard condition. It is found that the proposed two-phase thermodynamic (2PT) model can provide converged internal rotation entropy from a short, about 20 ps, trajectory of molecular dynamic simulation. The 2PT determined entropy is found to fall between those based on the theoretical models by Pitzer and Gwinn and by Truhlar. Our results suggest that this method can be a powerful method for understanding thermodynamic driving forces for conformation changes of flexible molecular structures.

---

1. Introduction

Understanding the entropy and free energy changes accompanied with the conformational change of large flexible molecules, such as polymers and proteins, is of great importance as it may elucidate the thermodynamic driving forces for conformation transitions. While existing methods¹ (e.g., thermodynamic integration,² thermodynamic perturbation theory,^3,4 Widom particle insertion,⁵ overlapping distribution method,⁶ umbrella sampling,⁷ etc.) can be used for extracting thermodynamic information from molecular simulations, most of them require a rather long simulation time or additional simulations for a hypothetical path of integration. Therefore, specific techniques have been developed for the treatment of conformational entropy of flexible molecules. One class of methods^8–21 determines the properties based on the equilibrium molecular geometry and vibrational frequencies of isolated molecules. For condensed phase systems, Benjamin J. Killian et al.,^22,23 showed that the entropy associated with dihedral torsions can be determined based on the torsional angle probability distribution extracted from Monte Carlo or molecular dynamics simulations. The challenges here are the need for higher order corrections due to the coupling among torsions and slow convergence of the probability distribution for large molecules. The quasiharmonic approximation²⁴ is a popular method for more rapid convergence for configurational entropy of flexible molecules in condensed phases. However, such a method is known to result in significant errors for low frequency tortional motions.^25,26

A new method, the two-phase thermodynamic (2PT) model,^27–29 has recently emerged and shown to provide converged thermodynamic properties based on a short, about 20 ps, trajectory of regular MD simulations, for a variety of systems.^30,31 It has also been used to understand many interesting phenomena, such as the penetration of water in carbon nanotube,^32,33 water dynamics in biological systems,^34,35 amorphous polyethylene glycol–water mixtures,³⁶ sugar at the membrane,³⁷ cellodextrin chain.³⁸ Despites its success, the theory for describing the thermodynamic properties of conformationally flexible molecules has not be developed.

The 2PT method determines the thermodynamic properties of a system from its vibrational density of states (DoS or normal modes), obtained from the Fourier transform of the velocity autocorrelation function. For solids, all the vibrational modes are nearly harmonic and the thermodynamic properties can be calculated by treating each mode as an independent harmonic oscillator (e.g., the Einstein–Debye model³⁹). For fluids, the low frequency modes, corresponding to diffusion and libration, are highly anharmonic. In 2PT, the DoS of fluids are considered to be superposition of that of a solid-like component (containing all the harmonic modes) and that of a gas-like component, which include all the low frequency anharmonic modes. Such a separation of the normal modes allows for applying suitable quantum statistics to different modes, and therefore accurate thermodynamic properties can be obtained.

The idea of decomposing the DoS into gas and solid-like components has been very successful in considering the anharmonic effects in translational and rotational motions, as evidenced by the various applications noted previously. However, for large flexible molecules, such polymers and proteins, the conformational changes as a result changes in dihedral torsions could also be highly anharmonic and needs to be properly taken into account.⁴⁰

In this work, we develop the 2PT theories for internal rotations. At low temperatures, internal rotations are hindered by the relatively high torsional barrier and the motions are nearly harmonic. As temperature increases and the kinetic energy overcome the rotational barrier, the dihedral torsion becomes a free internal rotor. Therefore, the DoS of internal rotations can be regarded as superposition of that of a free rotor and a harmonic oscillator with different weighting at different temperatures. Our results show that the 2PT can be as accurate as other theoretical methods developed by Pitzer and Gwinn (1942),⁴¹ Truhlar (1990),⁴² and McClurg et al. (1997).⁴³ Combined with its efficiency, the 2PT model is a powerful tool for study thermodynamic properties associated with conformation changes in flexible molecules.

2. Method and theory

2.1 The vibrational density of states

The vibrational density of states (DoS) of component i, is defined as the mass weighted sum of velocity spectral density from all atoms in the system,⁴⁴where m_j is the mass of atom j. The velocity spectral density s_j^k(ν) of atom j in the k-th coordinate (k = x, y, and z in the Cartesian coordinate) is determined from the square of the Fourier transform of the velocities as

The DoS can be also calculated from the Fourier transform of the velocity autocorrelation function (VAC).⁴⁴

2.2 Thermodynamic properties from two-phase thermodynamic (2PT) model

For polyatomic species, the total atomic velocity can be decomposed into translational, rotational and vibrational velocities (v_tot = v_trans + v_rot + v_vib).²⁸ Consequently, the total density of states is the sum of contributions from molecular translation (S_trn(ν)), rotation (S_rot(ν)), and vibration (S_vib(ν))

S(ν) = S_trn(ν) + S_rot(ν) + S_vib(ν) (3)

where S_trn(ν) is determined from the center of mass velocities of all the molecules (v_trans) and S_vib(ν) is determined from intramolecular vibration velocities (v_vib). The rotational density of state is determined from the angular velocity (in analogous to eqn (1) and (2))where I_l^k is the k-th principle moment of inertia of molecule l; ω_l^k is the angular velocity along the k principal axis; and M is the total number of molecules in the system. The density of state of translation and rotation can be further decomposed into a gas and a solid-like component

S_i,trn(ν) = S^g_i,trn(ν) + S^S_i,trn(ν)(5)

S_i,rot(ν) = S^g_i,rot(ν) + S^S_i,rot(ν) (6)

where the gas-like diffusive component S^g_i,trn(ν) corresponds to 3N^g_i,trn = 3f_i,trnN_i degrees of freedom with f_i,trn being the gas fraction of component i and the remainder, S^s_i,trn(ν), describes a solid-like part (non-diffusive) of component i for translation. Therefore, there are 3N^s_i,trn = 3N_i,trn − 3N^g_i,trn = 3N_i,trn(1 − f_i,trn) solid-like degrees of freedom for component i. The degree of freedom for rotation is defined in the similar way. The thermodynamic properties P_i of the system are determined from the individual DoS components with proper weighting functions.²⁷where and are the weighting functions of a harmonic oscillator and the corresponding the gas part of component i, respectively.

The decomposition of DoS is achieved by considering the gas-component as a hard sphere fluid, whose density of states is known

The solid component S^s_i(ν) is determined from subtracting S^g_i(ν) from the total density of state S_i(ν). The DoS for the gas component is completely determined with two parameters: s_0,i and f_i (for translation and rotation, respectively). In order to include all the diffusive modes to the gas component, s_0,i is set to be zero frequency DoS value for component i, S_i(0). This ensures that the solid component has no contribution to the diffusivity. The “fluidicity” factor f_i that determines the conceptual partition of each component between solid and gas parts can be calculated from the equation below (readers are referred to ref. 27 for derivation details).

where Δ_i is some normalized diffusivity, whose value depends on the temperature, volume, the particle mass and s⁰_i.

In the case of pure fluids, V_i is the same as the system total volume V. For a mixture, the system volume is shared by all components, and the partial molar volume (V̄_i) should be used, i.e.,²⁹

V_i = N_iV̄_i (12)

In addition, the original 2PT model treats all the vibration motions as harmonic oscillations. Therefore, S_i,vib(ν) = S^S_i,vib(ν). In the next section, the internal rotation, or the gas-portion of vibration will be considered.

To complete the 2PT model, we need to specify the weighting functions. The solid-like portion is calculated from the quantum partition function, which gives the harmonic oscillator weighting functions for energy, entropy, and Helmholtz free energy as follows⁴⁴

where β = 1/kT and h is Planck constant. For the gas-component, the weighting functions are derived from the corresponding properties of hard sphere gas for translation²⁷and free rigid rotor for rotation,where S^HS_i and S^R_i are the entropy of hard sphere and rigid rotors, respectively.²⁸

Once the decomposition of the DoS is done, the partial molar thermodynamic properties of each component can be obtained

where m = trn (center of mass translation), rot (rotation), or vib (intramolecular vibration). E_0,i is the reference energy and takes the form of

Ē_0,i = Ē^MD_i − β⁻¹3N_i (1 − 0.5f_i,trn − 0.5f_i,rot) (17)

The thermodynamic properties are determined as a sum of contributions from the corresponding translation, rotation and vibration.

Ē_i = Ē_0,i + Ē_i,trn + Ē_i,rot + Ē_i,vib (18a)

S̄_i = S̄_i,trn + S̄_i,rot + S̄_i,vib (18b)

Ā_i = Ē_0,i + Ā_i,trn + Ā_i,rot + Ā_i,vib (18c)

2.3 Density of states (DoS) for internal rotations

In order to identify the normal modes associated with internal rotations, it is necessary to use the internal coordinates. In other words, we would like to express the intramolecular vibrations (v_vib) in terms of contributions from bond stretching, angle bending, dihedral torsion and improper torsion, etc. Assuming infinitesimal amplitudes of vibration, the atomic displacement in Cartesian coordinates (x̲) depends on the in internal coordinates (q̲)⁴⁵where is a matrix whose elements are the change in internal coordinates corresponding to an infinitesimal perturbation in Cartesian coordinates, i.e.,

Note that there are 3N_atom degrees of freedom in Cartesian coordinates (x̲) and 3N_atom − 6 (or 3N_atom − 5 for linear molecules) degrees of freedom in internal coordinates (q̲). The dimension of q̲ are made equal to 3N_atom by adding six dummy internal coordinates T_x, T_y, T_z, R_x, R_y and R_z.^45–47 The additional rows in matrix has the following relationships with x̲.

and

Eqn (21) and (22) correspond to three translational and three rotational coordinates, respectively. These relationships are equivalent to the well-known Eckart conditions with M_tot^−1/2 and I^−1/2 being introduced as normalized factors. Because internal motions do not change molecule's overall translation and rotation, it is worth noting that T_x = T_y = T_z = R_x = R_y = R_z = 0. Therefore, matrix is a (3N_atom × 3N_atom) square matrix. The 3N_atom − 6 internal degrees of freedom are chosen to be N_atom − 1 bond stretching, N_atom − 2 angle bending and N_atom − 3 dihedral torsions, expressed in terms of the Z-matrix.⁴⁸ The elements in are obtained from numerical differentiation in terms of the change in internal coordinates in Z-matrix with respect to perturbation in Cartesian coordinates. Once the matrix is obtained, the velocity in the internal coordinates () is determined from⁴⁹

where the velocity in Cartesian coordinates is equivalent to v_vib obtained from velocity decomposition of the MD trajectory.

In additional to the velocities in the internal coordinates (), the corresponding “mass” for each internal mode is also needed for the calculation of the DoS. This is obtained from the fact that the kinetic energy (T) is independent of coordinate system, i.e.,

where is the mass matrix (M_ij = m_iδ_ij with m_i being the mass of atom associated with degree of freedom i), is the well-known Wilson G-matrix,⁵⁰ and its inverse can be regarded as effective mass in internal coordinates

The mass weighted internal velocity is calculated from

The DoS associated with the internal vibrations are

where is the mass weighted internal velocity of molecule m in n th internal coordinates. The total internal vibrational DoS is the sum of all the internal degrees of freedom and molecules.

This is the vibration density of state calculated from internal coordinates.

2.4 Two phase thermodynamic model for hindered rotors

Following the idea of 2PT, the DoS from the internal rotations is decomposed into a gas-like (free rotor) and solid-like (harmonic vibration) portion. Eqn (27) is used to represent the DoS of the free rotator, with the normalized diffusivity and the fluidicity factor also calculated from eqn (10) and (11).

2.4.1 Weighting functions of free rotor. The weighting functions of free internal rotors can be derived from the partition function

E^free = 1/2RT (29d)

C^free_P = 1/2R (29e)

where β is the reduced temperature (1/kT), h is Planck constant, σ is symmetric number, I_r is the reduced moment of inertia. The weighting functions for free rotation are

The reference energy in eqn (17) is modified as

Ē_0,i = Ē^MD_i − β⁻¹3N_i(1 − 0.5f_i,trn − 0.5f_i,rot − 0.5f_i,free) (31)

It is useful to note that many efforts have been made in developing a single partition function for hindered rotors. The representatives ones are the Pitzer and Gwinn model developed in 1942,⁴¹

the Truhlar model in 1990,⁴²and more recently the model by McClurg et al. in 1997 (ref. 43)where is the partition function of a quantum harmonic oscillator, is that for classical harmonic oscillator, u = βhν, and I₀ is the modified Bessel function. V₀ is the hindered rotor barrier height, which can be obtained from V₀ = 8πI_rν²/σ². The thermodynamic properties determined from the three partition functions are provided in the Appendix.

2.5 Reduced moment of inertia

The reduced moment of inertia for internal rotations can be rigorously determined from the equilibrium structure (minimum energy) of a molecule according to the theories established by Pitzer and co-workers.^41,51–53 For non-equilibrium structures (as the case for most molecules in MD simulations), East and Radom⁵⁴ proposed different levels of approximations to obtain reduced moment of inertia I^(m,n). The superscript (m, n) is used for bookkeeping different approximations. The n indicates the level of approximation for a rotor attached to a fixed frame before reduction, while the m indicates the level of approximation of the coupling reduction. Readers are referred to the original paper⁵⁴ for details. In this work, we choose n = 1, the moment of inertia of the rotating group is computed about the axis containing the twisted bond, and m = 2, the reduced moment due to coupling with overall molecular rotation is approximated by

One can arbitrarily choose either end as “left (L)” or “right (R)” of the twisted bond to be the rotating group. Note that when m = 3 and n = 4, the reduced moment of inertia is equivalent to that from Pitzer and Gwinn.⁵¹ For multi-torsional molecules, the reduced moment of inertia of each torsion can be defined in a similar manner.

3. Computational details

The vibrational entropy of methanol and ethane are used to examine the accuracy of the proposed method. These molecules are chosen for their different symmetries (symmetric tops for ethane and asymmetric tops for methanol). The interactions between different particles are described by the OPLS-AA force field.^55,56 Open software GROMACS⁵⁷ is used for molecular dynamic simulations. For each simulation, 256 molecules are randomly inserted into a 3-dimensional periodic box of desired density (0.79 g cm⁻³ for methanol and 0.54 g cm⁻³ for ethane). The system is first stabilized with energy minimization and then equilibrated under constant temperature, volume and number of particles (NVT) for 1 ns with a timestep of 1 fs. Additional 100 ps was subsequently performed and the trajectory file was recorded every 4 fs for the 2PT analysis. Nose–Hoover thermostat with time constant 0.2 ps was used to control the temperature. The NVT simulations were performed at different temperatures ranging from 200 K to 800 K. For methanol, the box length is 25.82 Å, while for ethane, the box length is 28.63 Å. The nonbond and electrostatic cutoff are both 9 Å. Particle-Mesh Ewald (PME) is used to calculate long range interactions. The Fourier spacing is 1.2 Å, and the PME order is 4. Additional NPT simulations are performed at standard condition (298.15 K and 1 atm) in order to compare results with experiments. Once the simulation is finished, the atomic velocity in each recorded time step is further decomposed into translational, rotational and vibrational velocities. The vibrational velocity in Cartesian coordinate is then converted to that in internal coordinate (eqn (23)). Once internal vibrational velocity and effective mass (eqn (26)) are obtained, the internal DoS are determined by eqn (27) and (28). The reduced moment of inertia needed to be used in partition function is calculated by eqn (35). For Pitzer and Gwinn, Truhlar, and McClurg's partition function approximations, the internal DoS is used directly in the integration. For the 2PT method, the internal DoS need to be separated into gas-portion (free rotor) and solid-portion (harmonic oscillator).

Note that for better statistics, internal modes of the same types are combined, especially dihedral rotations. Take methanol for example. The methyl group has three dihedral rotations of the same type. Thus, the s_0,i in eqn (11) is the summation of the three, and T is the average of the three. The rest of the calculation is similar to that for translation and rotation. Moreover, the gas–solid decomposing is needed for the dihedral torsions. Bond stretching and angle bending are treated as harmonic oscillations.

4. Results and discussion

4.1 The density of state of internal rotation

Fig. 1 illustrates the DoS calculated from vibrational velocities (v_vib). In other words, the contributions from molecular center of mass translation and rotations are removed in Cartesian (green line) and internal coordinates (red and blue lines) for ethane at 800 K. The peaks around 3000 cm⁻¹ correspond to CH₃ stretching, the one around 1000 cm⁻¹ is CC stretching, the one around 250 cm⁻¹ is the torsion, and others range from 700 to 1800 cm⁻¹ are associated with the deform and rock of CH3 (angle bending and coupled torsion and angle bending). In Cartesian coordinates, the hindered rotation is featured by a broad peak around 250 cm⁻¹ and a finite zero frequency intensity (see the green line in inset of Fig. 1). The use of internal coordinates allows for the deconvolution of the internal rotation modes from others and treating it based on its angular velocity. Therefore, the hindered internal rotation is featured by a rapid decay from zero-frequency intensity followed by a peak near 300 cm⁻¹ in internal coordinates (red line in inset of Fig. 1). This behavior is similar to molecular rotation observed in previous works.²⁸ The zero intensities of other vibrational modes (blue line) at low frequencies (<400 cm⁻¹) indicate that the torsional mode is properly removed. The torsional peaks at higher frequencies (red line) are coupled torsional modes with other vibrations. These coupled modes are treated as harmonic motions both in 2PT and other theories. Also shown in the inset of Fig. 1 are the gas (blue dashed line) and solid (purple dashed line) components of the internal rotation. As can be seen that the 2PT method nicely attributes the low-frequency modes to the gas-like component and the peak around 300 cm⁻¹ to solid-like component.

Fig. 1 Vibrational DoS for ethane at 800 K from Cartesian coordinate (green) and vibrational DoS from internal coordinate (red: internal rotation; blue: other vibrations; light blue: gas-component of internal rotation; purple: solid-component of internal rotation).

4.2 The entropy of internal rotation

Fig. 2(a) and 3(a) show the internal rotation entropy for methanol and ethane over a range of temperatures (200 K to 800 K). For convenience, the entropy of harmonic approximation (i.e., treating the internal rotation as harmonic oscillation) is subtracted in these figures. The numerical values are listed in Tables 1 and 2. It is expected that at low temperatures, the internal rotation behaves more like a harmonic oscillator, so that the entropy (referenced to harmonic approximation) should be small. On the other hand, at high temperatures, the internal rotation approaches to a free rotor. Since the harmonic approximation overestimates entropy of hindered rotation at high temperatures (see below), the entropy difference in Fig. 2(a) and 3(a) decreases with increasing temperature.

Fig. 2 (a) The torsional entropy of ethane from methods of Pitzer and Gwinn (circle), Truhlar (square), 2PT (diamond) using internal DoS in reference to harmonic approximation value. Also shown are the weighting function of entropy (W_s) for free rotor (black), harmonic oscillator (blue), Pitzer and Gwinn (green), Truhlar (pink), torsional DoS (red), gas component (light blue), and solid component (purple) for ethane at (b) 200 K and (c) 800 K.

Fig. 3 (a) The torsional entropy of ethane from methods of Pitzer and Gwinn (circle), Truhlar (square), 2PT (diamond) using internal DoS in reference to harmonic approximation value. Also shown are the weighting function of entropy (W_s) for free rotor (black), harmonic oscillator (blue), Pitzer and Gwinn (green), Truhlar (pink), torsional DoS (red), gas component (light blue), and solid component (purple) for methanol at (b) 200 K and (c) 800 K.

Table 1 Torsional energy and entropy of ethane from 200 K to 800 K and constant density of 0.54 g cm⁻³ determined based on Pitzer and Gwinn (PG),⁴¹ Truhlar,⁴² McClurg, Flagan, and Goddard (MFG),⁴³ 2PT, harmonic approximation (HO), and free rotor (free)

Temperature	200 K	300 K	400 K	500 K	600 K	700 K	800 K
E^PG (kJ mol^−1)	186.44	187.74	190.27	194.43	199.60	206.18	213.65
E^Trular (kJ mol^−1)	186.24	187.54	190.08	194.19	199.24	205.64	212.86
E^MFG (kJ mol^−1)	185.48	186.78	189.31	193.47	198.65	205.23	212.71
E^2PT (kJ mol^−1)	186.19	187.49	190.07	194.22	199.31	205.75	213.00
E^HO (kJ mol^−1)	186.29	187.75	190.56	195.01	200.47	207.26	214.91
E^free (kJ mol^−1)	4.15	6.22	8.27	10.33	12.40	14.46	16.58
S^PG (J mol^−1 K^−1)	6.05	11.53	19.06	27.99	37.40	47.24	57.18
S^Truhlar (J mol^−1 K^−1)	3.65	9.16	16.81	25.75	35.04	44.67	54.33
S^MFG (J mol^−1 K^−1)	6.05	11.53	19.06	27.99	37.40	47.24	57.18
S^2PT (J mol^−1 K^−1)	4.07	10.15	18.35	27.58	37.07	46.85	56.59
S^OH (J mol^−1 K^−1)	4.07	10.60	19.91	30.58	41.67	52.62	63.62
S^free (J mol^−1 K^−1)	30.78	39.15	45.04	49.64	53.41	56.59	59.42
Torsional fluidicity	0.01	0.04	0.07	0.11	0.14	0.16	0.19

---

Table 2 Torsional energy and entropy of methanol from 200 K to 800 K and constant density of 0.79 g cm⁻³ determined based on Pitzer and Gwinn (PG),⁴¹ Truhlar,⁴² McClurg, Flagan, and Goddard (MFG),⁴³ 2PT, harmonic approximation (HO), and free rotor (free)

Temperature	200 K	300 K	400 K	500 K	600 K	700 K	800 K
E^PG (kJ mol^−1)	128.74	129.42	131.18	133.71	137.16	141.53	146.04
E^Trular (kJ mol^−1)	128.69	129.27	130.91	133.37	136.81	141.16	145.64
E^MFG (kJ mol^−1)	128.25	128.93	130.69	133.22	136.68	141.05	145.56
E^2PT (kJ mol^−1)	128.54	129.04	130.64	133.12	136.59	141.03	145.56
E^HO (kJ mol^−1)	128.75	129.41	131.17	133.83	137.53	142.21	147.06
E^free (kJ mol^−1)	2.50	3.70	4.96	6.18	7.40	8.71	9.92
S^PG (J mol^−1 K^−1)	2.44	6.54	12.86	19.70	26.54	33.64	40.11
S^Truhlar (J mol^−1 K^−1)	1.78	5.42	11.25	17.84	24.62	31.71	38.18
S^MFG (J mol^−1 K^−1)	2.44	6.54	12.86	19.70	26.54	33.64	40.11
S^2PT (J mol^−1 K^−1)	2.21	5.90	11.86	18.72	25.77	33.24	39.90
S^OH (J mol^−1 K^−1)	2.56	6.89	13.37	20.78	28.44	36.53	44.02
S^free (J mol^−1 K^−1)	20.69	25.51	29.08	31.75	33.96	35.93	37.54
Torsional fluidicity	0.039	0.06	0.08	0.09	0.11	0.12	0.14

---

In all cases, the internal entropies determined based on Pitzer and Gwinn are the largest, those based on Truhlar are the lowest, and those based on 2PT fall in between. The different behaviors from the three methods can be understood from the weighting function of entropy.

Fig. 2(b) and (c), 3(b) and (c) show the entropy weighting function of entropy and DoS of internal rotation. The weighting function of free rotor (black line, eqn (30b)) is a constant and that of a harmonic oscillator (blue line, eqn (13b)) decreases monotonically from infinity at zero frequency to zero at infinite frequency. Therefore the weighting function of free rotor always intersects with that of harmonic oscillator at some intermediate frequency. The entropy weighting function based on Truhlar (purple line) is always lower than those from free rotor and harmonic oscillator. Therefore, the Truhlar value of entropy should be the lower bound for internal rotation. The weighting function of Pitzer and Gwinn (green line) always lies above that of Truhlar, and intersects with that of the harmonic oscillator. Therefore, the entropy of internal rotation from Pitzer and Gwinn is always higher than that from Truhlar and may exceed that of the harmonic approximation for hindered rotations. Unlike these methods, the entropy of internal rotation in 2PT is determined from the different weightings of gas (free rotor) and solid (harmonic oscillator) component, determined by the fluidicity factor. It is the difference in the weighting functions and the distribution of DoS from internal rotation that leads to the different results seen in Fig. 2(a) and 3(a).

For example, the DoS of internal rotation of ethane has a libration peak at around 310 cm⁻¹ and a small rotation diffusion intensity of 0.1 cm at 200 K (Fig. 2(b)). As the temperature increases to 800 K (Fig. 2(c)), the libration peak broadens with intensity reduced and the intensity of rotational diffusivity increases to 1.4 cm. The coexistence of rotational diffusion and libration indicates hindered rotation at 800 K. The entropy is calculated as the integration of product of the DoS and the weighting function (eqn (31)). At 200 K, the libration peak locates in the region where the weighting function of W^PG_S > W^T_S ≅ W^HO_S and the calculated entropy follows the same order. In this case, the fluidicity of internal rotation is nearly zero (0.01), and therefore the 2PT entropy is close to that of Truhlar. This also explains the positive deviation of Pitzer and Gwinn value from that of harmonic oscillator at low temperatures (Fig. 2(a)). At 800 K, the libration peak red shifts to about 290 cm⁻¹, close to the cross point of W^PG_S and W^HO_S. Interestingly, the 2PT decomposition (with fluidicity increased to 0.19) provides a weighting between free rotor and harmonic oscillator, resulting in an internal rotation energy very close to that of Pitzer and Gwinn.

The DoS of internal rotation of methanol behaves quite differently from that of ethane in the range of temperature from 200 K to 800 K. The rotational diffusion (zero frequency intensity) is not negligible (compared to the libration peak) at 200 K. The libration peak (about 140 cm⁻¹) is very close to the cross point of W^PG_S and W^HO_S. The 2PT entropy is close to the Pitzer and Gwinn entropy in this case (similar to the situation of ethane at 800 K). At 800 K, the libration peak broadens and is rather flat up to 500 cm⁻¹. In this case, the entropy from different methods reflects the value of weighting function throughout the spectrum.

It is noteworthy that hydrogen bonding results in the distinct mode dependence on temperature in ethane and methanol. The fluidicity of methane (0.01) is lower than that of methanol (0.04) at 200 K but becomes higher (0.19 vs. 0.14) at 800 K. In other words, rotational “melting” (from hindered rotor to free rotor) is more difficult in hydrogen bonding fluids.

From these analyses, we observe that the 2PT entropies for internal rotation are bounded by that of free rotor at high temperatures (fluidicity goes to unity) and by that of harmonic oscillator at low temperatures (fluidicity goes to zero). At intermediate temperatures, the internal rotation entropy is attenuated by the value of fluidicity, which is determined by the rotational diffusivity (or S(0)), and the density of the fluid (eqn (10) and (11)). The entropy from Truhlar and Pitzer and Gwinn also converges to that of free rotor at high temperatures, in which case the DoS concentrates on low frequencies. However, at low temperatures, the Truhlar entropy is always lower than that of harmonic oscillator. The Pitzer and Gwinn entropy may either be higher or lower than that of harmonic oscillator depending on the distribution of DoS at low temperatures.

4.3 Two-phase thermodynamic model for molecules containing multiple torsions

Liquid butane and hexane are used to illustrate the proposed method for larger, flexible molecules. Fig. 4 shows the density of states of the three torsions (indicated on Fig. 4(a)) for liquid butane at 298.15 K and 0.599 g cm⁻³. The DoS of the first and third torsions (Fig. 4(a) and c) are similar due to conformational symmetry. However, there is a significant red shift of the first peak of the second torsion (110 cm⁻¹) compared to that of the first and third torsions (235 cm⁻¹). The simulation results agree well with experiment, where the center torsion has a frequency around 100 (cm⁻¹),⁵⁸ and the other two torsion are at around 197–225 (cm⁻¹). The results show that the proposed method correctly differentiates different torsions in butane. Furthermore, the fluidicity and entropy of the three torsions are f₁ = 0.036, S₁ = 6.3 (J mol⁻¹ K⁻¹); f₂ = 0.019, S₂ = 8.1 (J mol⁻¹ K⁻¹); f₃ = 0.03, S₃ = 6.4 (J mol⁻¹ K⁻¹), respectively. The lower fluidicity indicates the motion of second rotation is significantly hindered by the two heavy methyl groups, leading to a larger entropy contribution.

Fig. 4 The torsional DoS of the first (a) second (b), and third (c) dihedral angle for liquid butane at 298.15 K (red: total dihedral DoS, light-blue: gas-component, green: solid-component).

Fig. 5 shows the density of state of all 5 torsions for liquid hexane at 298.15 K and 0.655 g cm⁻³. Generally, the characteristic is very similar to that from butane. However, it is observed that a bimodal peak exists for the second and the fourth hindered rotors, which is considered a coupling between the adjacent torsions. For longer chains, it is believed that the coupling between adjacent torsions is more significant, which renders it difficult to analyze. Under this circumstance, the 2PT internal method is an efficient tool to use. The fluidicity and entropy of the 5 torsions are f₁ = 0.028, S₁ = 9.61 (J mol⁻¹ K⁻¹); f₂ = 0.029, S₂ = 14.05 (J mol⁻¹ K⁻¹); f₃ = 0.032, S₃ = 14.17 (J mol⁻¹ K⁻¹); f₄ = 0.03, S₄ = 14.07 (J mol⁻¹ K⁻¹); f₅ = 0.029, S₅ = 9.62 (J mol⁻¹ K⁻¹), respectively. Contrary to butane, the fluidicity of the end and the middle torsions are approximately the same, indicating that the methylene groups can rotate more freely as the chain length increases. Furthermore, the middle torsions has a more significant contribution to the entropy (about 14 J mol⁻¹ K⁻¹) compared to the torsions near the two ends (about 9.6 J mol⁻¹ K⁻¹).

Fig. 5 The torsional DoS of the first to fifth dihedral angles (a to e, respectively) for liquid hexane at 298.15 K (red: total dihedral DoS, light-blue: gas-component, green: solid-component).

Fig. 6 The 2PT torsional entropy determined from different lengths of a trajectory for methanol at 200 K (circle) and 800 K (square).

Table 3 shows the entropy components (translation, rotation, torsion, non-tortional vibrations) of methanol, ethane, butane and hexane determined at standard condition. In general, the 2PT determined entropies are in good agreement with the experimental results.⁵⁹ The vibrational entropy is small compared to translational and rotational contributions, especially for small molecules. The contribution of the internal rotation to entropy becomes more significant for larger, flexible molecules, such as hexane. As noted previously, the harmonic approximation overestimates the torsional entropy for nearly freely rotors but underestimates it for hindered rotors. The less significant contribution of torsional entropy and the counterbalance of harmonic approximations results in similar vibrational entropy value for small sized molecules. However, for longer chain molecules, such as hexane, the hindered rotation of torsion dominates the vibrational entropy and the results from 2PT shows a significant improvement over that form harmonic approximation. The discrepancy in entropy of hexane from 2PT and experiment might be a result of the performance of force field used. The OPLS-AA force field was recently shown to be less accurate to reproduce experimental data for long hydrocarbons.⁶⁰

Table 3 Entropy of methanol, ethane, butane, and hexane at standard condition determined from 2PT model based on Cartesian and internal coordinates, respectively (all numbers are in unit of J mol⁻¹ K⁻¹)

Method	State of aggregation^a	Strn	Srot	Svib		S^cal		S^exp^b
		2PT	2PT	2PT^c	HO^d	2PT	HO^d
a G: gas, L: liquid.b NIST reference database number 69.^59c The numbers (A + B) in parenthesis indicates entropy contribution from dihedral torsions (A) and non-dihedral torsions (B).d Entropy determined based harmonic approximations.
Methanol	L	64.8	46.0	5.9(4.0 + 1.9)	6.6	116.7 ± 1.1	117.4 ± 1.2	127.0
Ethane	G	131.8	84.6	10.6(7.8 + 2.8)	11.0	227.0 ± 0.3	227.4 ± 0.3	229.0
Butane	G	156.7	104.9	48.0(32.1 + 15.9)	47.6	309.6 ± 0.3	309.2 ± 0.3	310.2
Hexane	L	83.1	65.7	102.3(61.5 + 40.8)	96.1	251.2 ± 1.1	244.9 ± 1.1	289.5

---

4.4 Convergence of two-phase thermodynamic model

One attractive feature of 2PT model is its use of a very short sampling time. From the previous studies, converged properties are obtained from 20 ps trajectory of a MD simulation.^27–29 Here we examine the convergence of 2PT internal method. Fig. 6 shows the 2PT vibrational entropy for methanol at 200 K and 800 K. In both cases, the entropy converges within 20 ps. However, at low temperatures (200 K), the entropy converges within 5 ps. These results indicate that the thermodynamic properties from internal rotations also converge efficiently within 20 ps.

5. Conclusion

In this work, the 2PT method has been extended to consider internal rotations of flexible molecules. The vibrational DoS obtained from internal coordinates is divided into gas-like (free rotor) and solid-like (harmonic oscillator) components, allowing for an efficient determination of the thermodynamic properties via suitable weighting functions for each component. This method can be used to consider the anharmonic effects associated with conformation transitions of flexible molecules. Combined with its efficiency, the 2PT model can be a powerful tool study thermodynamic properties of complex systems containing flexible, macromolecules such as polymers, proteins, and DNAs.

Appendix: weighting functions for hindered rotors

Pitzer and Gwinn:

Truhlar:

McClurg, Flagan, and Goddard:

Acknowledgements

The authors like to acknowledge the funding support from the National Science Council of Taiwan under grants NSC 100-2221-E-002-175 and 101-2628-E-002-014-MY3. The computation resources from the National Center for High-Performance Computing of Taiwan and the Computing and Information Networking Center of the National Taiwan University are acknowledged.

References

1. D. Frenkel and B. Smit, Understanding Molecular Simulation From Algorithms to Applications, Academic press, New York, 2002.
2. D. A. Kofke, Mol. Phys., 1993, 78, 1331–1336.
3. M. S. Wertheim, J. Stat. Phys., 1986, 42, 477–492.
4. M. S. Wertheim, J. Stat. Phys., 1984, 35, 35–47.
5. B. Widom, J. Chem. Phys., 1963, 39, 2808–2812.
6. C. H. Bennett, J. Comput. Phys., 1976, 22, 245–268.
7. G. M. Torrie and J. P. Valleau, J. Comput. Phys., 1977, 23, 187–199.
8. T. F. Miller and D. C. Clary, J. Chem. Phys., 2003, 119, 68–76.
9. T. F. Miller and D. C. Clary, J. Chem. Phys., 2002, 116, 8262–8269.
10. T. F. Miller and D. C. Clary, J. Phys. Chem. B, 2004, 108, 2484–2488.
11. T. F. Miller and D. C. Clary, Mol. Phys., 2005, 103, 1573–1578.
12. A. Chakraborty and D. G. Truhlar, J. Chem. Phys., 2006, 124, 224105.
13. V. van Speybroeck, R. Gani and R. J. Meier, Chem. Soc. Rev., 2010, 39, 1764–1779.
14. P. Vansteenkiste, D. Van Neck, V. Van Speybroeck and M. Waroquier, J. Chem. Phys., 2006, 125, 049902.
15. P. Vansteenkiste, D. Van Neck, V. Van Speybroeck and M. Waroquier, J. Chem. Phys., 2006, 124, 044314.
16. P. Vansteenkiste, T. Verstraelen, V. Van Speybroeck and M. Waroquier, Chem. Phys., 2006, 328, 251–258.
17. J. J. Zheng, R. Meana-Paneda and D. G. Truhlar, Comput. Phys. Commun., 2013, 184, 2032–2033.
18. J. J. Zheng, S. L. Mielke, K. L. Clarkson and D. G. Truhlar, Comput. Phys. Commun., 2012, 183, 1803–1812.
19. J. J. Zheng and D. G. Truhlar, J. Chem. Theory Comput., 2013, 9, 1356–1367.
20. J. J. Zheng and D. G. Truhlar, Faraday Discuss., 2012, 157, 59–88.
21. J. J. Zheng, T. Yu, E. Papajak, I. M. Alecu, S. L. Mielke and D. G. Truhlar, Phys. Chem. Chem. Phys., 2011, 13, 10885–10907.
22. V. Hnizdo, J. Tan, B. J. Killian and M. K. Gilson, J. Comput. Chem., 2008, 29, 1605–1614.
23. B. J. Killian, J. Y. Kravitz and M. K. Gilson, J. Chem. Phys., 2007, 127, 024107.
24. M. Karplus and J. N. Kushick, Macromolecules, 1981, 14, 325–332.
25. J. Numata and E. W. Knapp, J. Chem. Theory Comput., 2012, 8, 1235–1245.
26. J. Numata, M. Wan and E. W. Knapp, Genome informatics, Int. Conf. Genome Inf., 2007, 18, 192–205.
27. S.-T. Lin, M. Blanco and W. A. Goddard, J. Chem. Phys., 2003, 119, 11792.
28. S.-T. Lin, P. K. Maiti and W. A. Goddard, J. Phys. Chem. B, 2010, 114, 8191–8198.
29. P. K. Lai, C. M. Hsieh and S. T. Lin, Phys. Chem. Chem. Phys., 2012, 14, 15206–15213.
30. S.-N. Huang, T. A. Pascal, W. A. Goddard, P. K. Maiti and S.-T. Lin, J. Chem. Theory Comput., 2011, 7, 1893–1901.
31. T. A. Pascal and W. A. Goddard, J. Phys. Chem. B, 2012, 116, 13905–13912.
32. Y. Zhao, C.-C. Ma, L.-H. Wong, G. Chen, Z. Xu, Q. Zheng, Q. Jiang and A. T. Chwang, J. Comput. Theor. Nanosci., 2006, 3, 852–856.
33. H. Kumar, B. Mukherjee, S.-T. Lin, C. Dasgupta, A. K. Sood and P. K. Maiti, J. Chem. Phys., 2011, 134, 124105.
34. V. Vasumathi and P. K. Maiti, Macromolecules, 2010, 43, 8264–8274.
35. B. Nandy and P. K. Maiti, J. Phys. Chem. B, 2011, 115, 217–230.
36. H. Shin, T. A. Pascal, W. A. Goddard and H. Kim, J. Phys. Chem. B, 2013, 117, 916–927.
37. J. H. Tian, A. Sethi, B. I. Swanson, B. Goldstein and S. Gnanakaran, Biophys. J., 2013, 104, 622–632.
38. P. K. GhattyVenkataKrishna, E. M. Alekozai, G. T. Beckham, R. Schulz, M. F. Crowley, E. C. Uberbacher and X. L. Cheng, Biophys. J., 2013, 104, 904–912.
39. D. A. McQuarrie, Statistical mechanics, University Science Books, 2000.
40. O. M. Becker, A. D. Mackerell Jr, B. Roux and M. Watamabe, Computational Biochemistry and Biophysics, Marcel Dekker, Inc. New York, 2001.
41. K. S. Pitzer and W. D. Gwinn, J. Chem. Phys., 1942, 10, 428–440.
42. D. G. Truhlar, J. Comput. Chem., 1991, 12, 266–270.
43. R. B. McClurg, R. C. Flagan and W. A. Goddard, J. Chem. Phys., 1997, 106, 6675–6680.
44. P. H. Berens, D. H. J. Mackay, G. M. White and K. R. Wilson, J. Chem. Phys., 1983, 79, 2375–2389.
45. E. B. Wilson, J. C. Decius and P. C. Cross, Molecular vibrations: The theory of infrared and Raman vibrational spectra, 1980.
46. P. R. Bunker and P. Jensen, Fundamental of Molecular Symmetry, 2005.
47. P. R. Bunker and P. Jensen, Molecular Symmetry and Spectroscopy, 1998.
48. J. Parsons, J. B. Holmes, J. M. Rojas, J. Tsai and C. E. M. Strauss, J. Comput. Chem., 2005, 26, 1063–1068.
49. M. E. Castro, A. Nino and C. Munoz-Caro, Comput. Phys. Commun., 2010, 181, 1469–1473.
50. E. B. Wilson, J. Chem. Phys., 1941, 9, 76–84.
51. K. S. Pitzer, J. Chem. Phys., 1946, 14, 239–243.
52. J. E. Kilpatrick and K. S. Pitzer, J. Chem. Phys., 1949, 17, 1064–1075.
53. J. C. M. Li and K. S. Pitzer, J. Phys. Chem., 1956, 60, 466–474.
54. A. L. L. East and L. Radom, J. Chem. Phys., 1997, 106, 6655–6674.
55. W. L. Jorgensen, D. S. Maxwell and J. Tirado-Rives, J. Am. Chem. Soc., 1996, 118, 11225–11236.
56. W. L. Jorgensen and J. Tirado-Rives, J. Am. Chem. Soc., 1988, 110, 1657–1666.
57. H. J. C. Berendsen, D. van der Spoel and R. van Drunen, Comput. Phys. Commun., 1995, 91, 43–56.
58. P. H. Berens, D. H. J. Mackay, G. M. White and K. R. Wilson, J. Chem. Phys., 1983, 79, 2375–2389.
59. C. W. NIST, National Institute of Standards and Technology, Reference Database Number 69 edn, 2000.
60. S. W. I. Siu, K. Pluhackova and R. A. Bockmann, J. Chem. Theory Comput., 2012, 8, 1459–1470.

---