Efficient Computation of Entropy and Other Thermodynamic Properties

Subscriber access provided by TULANE UNIVERSITY

B: Liquids, Chemical and Dynamical Processes in Solution, Spectroscopy in Solution

Efficient Computation of Entropy and Other Thermodynamic Properties for Two-Dimensional Systems Using Two-Phase Thermodynamic Model Sindhana Selvi Pannir Sivajothi, Shiang-Tai Lin, and Prabal K. Maiti J. Phys. Chem. B, Just Accepted Manuscript • DOI: 10.1021/acs.jpcb.8b07147 • Publication Date (Web): 10 Dec 2018 Downloaded from http://pubs.acs.org on December 18, 2018

Just Accepted “Just Accepted” manuscripts have been peer-reviewed and accepted for publication. They are posted online prior to technical editing, formatting for publication and author proofing. The American Chemical Society provides “Just Accepted” as a service to the research community to expedite the dissemination of scientific material as soon as possible after acceptance. “Just Accepted” manuscripts appear in full in PDF format accompanied by an HTML abstract. “Just Accepted” manuscripts have been fully peer reviewed, but should not be considered the official version of record. They are citable by the Digital Object Identifier (DOI®). “Just Accepted” is an optional service offered to authors. Therefore, the “Just Accepted” Web site may not include all articles that will be published in the journal. After a manuscript is technically edited and formatted, it will be removed from the “Just Accepted” Web site and published as an ASAP article. Note that technical editing may introduce minor changes to the manuscript text and/or graphics which could affect content, and all legal disclaimers and ethical guidelines that apply to the journal pertain. ACS cannot be held responsible for errors or consequences arising from the use of information contained in these “Just Accepted” manuscripts.

is published by the American Chemical Society. 1155 Sixteenth Street N.W., Washington, DC 20036 Published by American Chemical Society. Copyright © American Chemical Society. However, no copyright claim is made to original U.S. Government works, or works produced by employees of any Commonwealth realm Crown government in the course of their duties.

Page 1 of 44 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

Efficient Computation of Entropy and Other Thermodynamic Properties for Two-Dimensional Systems Using Two-Phase Thermodynamic Model Sindhana Selvi PS,† Shiang-Tai Lin,‡ and Prabal K. Maiti*† †Centre for Condensed Matter Theory, Department of Physics, Indian Institute of Science, Bangalore 560012, India ‡ Department of Chemical Engineering, National Taiwan University, Taipei 10617, Taiwan

1

ACS Paragon Plus Environment

The Journal of Physical Chemistry 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Abstract We present a method based on the two-phase thermodynamic model (2PT) to calculate entropy and free energy of various molecular systems in two dimensions (2D) using Molecular Dynamics (MD) simulations. The 2PT method has been used widely to calculate absolute entropy in a variety of molecular systems in three dimensions. When applying the idea to 2D systems, we found that the fluidicity that determines the decomposition of the vibrational density of states (DoS) into a solid-like and a gas-like component needs to be revised. The solid part is treated using quantum statistics and the gas part is treated as a hard-disk fluid. We validate this method by computing thermodynamic properties of a two-dimensional Lennard-Jones fluid over a range of densities and temperatures and find excellent agreement with these quantities computed from the equation of state. More importantly, this method allows for the calculation of entropy and free energy of 2D systems efficiently from a single MD trajectory of less than 50 ps; therefore, it can be a new, powerful way of assessing the thermodynamic properties in 2D problems.

1

Introduction

Our ability to understand and predict complex chemical phenomena from molecular simulations improves as our methods for free energy computation become more accurate and efficient. Thermodynamic integration (TI), 1 free energy perturbation (FP), 2 widom particle insertion 3 and umbrella sampling 4 are well-established methods for free energy computation. These methods allow direct computation of free energy difference between two states. Most of these methods are computationally demanding. Another serious drawback is the requirement of a well-defined reversible path connecting the two states which is sometimes difficult to get. 5 Calculation of absolute free energy requires access to the both enthalpy and entropy estimates. If sampling issues are taken care of by a variety of means, estimate of absolute enthalpy for a given phase is easily available from simulations. To compute absolute entropy 2


Page 2 of 44

Page 3 of 44 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


from simulations, few methods are available in literature. For example, to compute the configurational entropy of macromolecules the method based on covariance matrix 6 is often used. Another interesting approach to compute entropy is based on permutation reduction (PR). 7,8 Methods based on cell theory 9,10 and the two phase thermodynamic model (2PT) 11 are also used to compute entropy and are computationally very efficient. Two-dimensional (2D) systems have received world wide attention in recent years after the discovery of several materials such as graphene, 12 MoS2 , hexagonal boron-nitride sheet (hBN) , transition metal dichalcogenides. 13 Dimensional reduction gives rise to several unique properties which makes possible potential application of these materials in diverse fields. For example, the melting transition in 2D is very different than their 3D counterpart and believed to happen through an intermediate hexatic phase. 14 Fluids confined in 2D system of graphene, 15 16 lipid bilayer, MoS2 and hBN offers another class of systems which exhibit many interesting phenomena. For example, Srivastava et al. 17 showed that Lennard-Jones fluid confined inside slit-pore geometry of graphene show hexagonal close packed 2D solid phase with a very high tangential pressure near the wall. Similar phenomena was observed in a recent experiment by Vasu et al. 18 where they reported that trapped molecules inside 2D graphene confinement experience vdW pressure as high as 1 GPa. To understand microscopic origin of many such phenomena in 2D systems we need a very robust and efficient method to compute various thermodynamic properties such as entropy, free energy and pressure. The two-phase thermodynamic model (2PT model) is a fast and accurate method used to estimate the entropy and free energy of a system from short molecular dynamics (MD) runs. The 2PT method has been applied to compute the entropies and heat capacities of organic liquids 19,20 and the results are in good agreement with experiment. It has also been applied to compute the entropy of water with great success. 21–23 In this paper, we present a method based on 2PT for calculating entropy in two-dimensional systems. This new method will be referred to as 2PT(2D). The 2PT(2D) method has been tested by applying it to a two-dimensional Lennard-Jones (2D LJ) fluid and the results have been compared with the

3



Page 4 of 44

modified Benedict-Webb-Rubin equation of state (MBWR EOS) 24 of a 2D LJ fluid. The 2PT(2D) method uses a fluidicity parameter f to compute the entropy of a system. It defines the fluidicity parameter (the fraction of the system treated as a fluid) as the ratio of the diffusivity of the system to the diffusivity of a hard-disk fluid at the same temperature and density (in the zero pressure limit). However, the diffusivity of a twodimensional system depends on system size. 25 As a result, we might expect the entropy calculated from the 2PT(2D) method to depend on system size. A discussion on how the thermodynamic properties calculated by the 2PT(2D) method depend on system size is presented in Section 4.6. The paper is organized in the following way. In Section 2, we give the theoretical background of the 2PT(2D) method. In Section 3, we describe the computational details of the MD simulations for the two-dimensional LJ fluid. In Section 4, we compare the results from the 2PT(2D) method for a 2D LJ fluid with the equation of state. Finally, in Section 5, we conclude with a discussion on major benefits of the method and possible application.

2

Theory

Any given thermodynamic property P can be expressed in terms of the vibrational density of states (DoS) S(ν) and the weighting function for that thermodynamic property WPsystem (ν) as shown in eq(1). Thus any thermodynamic property of a system can be computed if the vibrational density of states and the weighting function for the system are known. Z P =

∞

S(ν)WPsystem (ν)dν

(1)

0

The details on what the DoS and weighting functions are, and how they can be obtained, are explained in Section 2.2 and Section 2.3. In the 2PT(2D) method, the DoS of the system is computed from MD simulations. However, analytical expressions are used for weighting functions. The 2PT(2D) method utilizes the weighting functions of known model systems 4


Page 5 of 44 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


(harmonic oscillator and hard-disk fluid) to compute thermodynamic properties of any given system.

2.1

Velocity autocorrelation function

The total velocity autocorrelation function (VACF) can be computed from MD simulations and serves as the input for the 2PT(2D) model. The total velocity autocorrelation function, C(t), is the mass weighted sum of the particle velocity autocorrelation functions, ckj (t).

C(t) =

N X 2 X

mj ckj (t)

(2)

j=1 k=1

ckj (t)

1 = lim τ →∞ 2τ

Z

τ

vjk (t0 )vjk (t0 + t)dt0

(3)

−τ

ckj (t) is the velocity autocorrelation function of the j th atom in the k th direction and N is the total number of particles in the system. The particle velocity autocorrelation function is related to the diffusion coefficient D by the Green-Kubo relations as shown below. 26 We obtain the equations below for a system containing only one type of particle with mass m. 1 D= 4

Z

∞

2 X

ck (t)dt

−∞ k=1

Z

N

∞

2

1 XX k mcj (t)dt −∞ N m j=1 k=1 Z ∞ 1 C(t)dt = 4mN −∞ 1 = 4

2.2

(4)

Vibrational density of states (DoS)

As defined by Lin et al. 11 “The density of states function S(ν) is defined as the distribution of vibrational normal modes of a system”. Properties of the density of states function are derived and explained in the three-dimensional 2PT paper. 11 In this section, we will just state the properties of S(ν). 5



Page 6 of 44

1. S(ν)dν is a measure of the number of normal modes with frequency between ν and ν + dν. Thus when we integrate S(ν) over all frequencies we will get the total degrees of freedom of the system. Z

∞

S(ν)dν = 2N

(5)

0

2. The vibrational density of states, S(ν), can be obtained from the Fourier transform of the total velocity autocorrelation function C(t). 2 S(ν) = kB T

Z

∞

C(t)e−i2πνt dt

(6)

−∞

3. The zero frequency value of the density of states S(0) is related to the diffusion coefficient D as shown below. From here on we will refer to S(0) as s0 for compactness of notation. Z ∞ 2 C(t)dt S(0) = kB T −∞ 2 = 4mN D kB T 8mN D = kB T

(7)

Hence, the diffusion coefficient can be expressed in terms of s0 as shown in eq(8). Here, ρ refers to the number of particles per unit area and T refers to the temperature.

D(T, ρ) =

s0 kB T 8mN

(8)

Since the velocity autocorrelation function C(t) and the density of states function S(ν) are related by eq(6) and the VACF can be obtained from MD simulations, the density of states function S(ν) can also be obtained from simulations.

6


Page 7 of 44 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


2.3

Weighting functions

Weighting functions are a measure of the per particle thermodynamic properties (at a particular frequency ν) of a system. Therefore, when one multiplies them with S(ν)dν, and integrates over all ν, one gets the thermodynamic property of interest. A point to note would be that we make our weighting functions dimensionless. Weighting functions of harmonic oscillator and hard-disk fluid are described in this section.

2.3.1

Harmonic oscillator

If each normal mode in the vibrational DoS of the system is treated as an independent harmonic oscillator, then the partition function Q of the whole system can be calculated using the harmonic oscillator partition function q(ν) given by eq(9). Z ln Q =

∞

S(ν) ln q(ν)dν

(9)

0

Either the quantum harmonic oscillator partition function q Q (ν) or the classical harmonic oscillator partition function q C (ν) can be used to obtain the total partition function of the system. Depending on whether the quantum harmonic oscillator expressions or the classical harmonic oscillator expressions are used, the 2PT(2D) method is referred to as 2PT(Q) or 2PT(C). exp(−βhν/2) 1 − exp(−βhν) 1 q C (ν) = βhν

q Q (ν) =

(10a) (10b)

The expressions for thermodynamic properties as obtained in the three-dimensional 2PT

7



Page 8 of 44

paper 11 are given below.

−1

Z

∞

E = V0 + β S(ν)WE (ν)dν 0 Z ∞ S(ν)WS (ν)dν S=k 0 Z ∞ −1 S(ν)WA (ν)dν A = V0 + β

(11a) (11b) (11c)

0

V0 = E M D − β −1 2N

(11d)

In the set of equations above, E is the energy, S is the entropy, A is the Helmholtz free energy and V0 is the reference energy. The reference energy V0 is computed by equating the energy of a system of 2N classical harmonic oscillators to the energy E M D obtained from molecular dynamics simulation. The weighting functions for the quantum harmonic oscillator are given by the following set of equations. βhν βhν + 2 exp(βhν) − 1 βhν − ln[1 − exp(−βhν)] WSQ (ν) = exp(βhν) − 1 1 − exp(−βhν) WAQ (ν) = ln exp(−βhν/2)

WEQ (ν) =

(12a) (12b) (12c)

The weighting functions for the classical harmonic oscillator are given below.

2.3.2

WEC (ν) = 1

(13a)

WSC (ν) = 1 − ln(βhν)

(13b)

WAC (ν) = ln(βhν)

(13c)

Hard-disk fluid

Since the frequency of normal modes have no relevance to the thermodynamic properties of a hard-disk fluid, the weighting functions of a hard-disk fluid system are independent of 8


Page 9 of 44 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


frequency ν. Hence, to determine the weighting functions of a hard-disk fluid system, we require expressions for the thermodynamic properties of a hard-disk fluid. The relationship between weighting functions and thermodynamic properties are given below, for the special case where the weighting functions are independent of frequency ν.

E

HD

=β

−1

WEHD

Z

∞

S(ν)dν β 1 = 2N kB T 2N 2 0

β HD E 2N 1 = 2

WEHD =

(14)

The equipartition theorem is used to obtain the energy E HD of the hard-disk fluid. The entropy weighting function is computed using the entropy per particle sHD (T, ρ) for a system of hard-disks.

S

HD

=

kB WSHD

Z

∞

S(ν)dν 0

WSHD

=

WSHD =

S HD 1 N 2kB

sHD 2kB

(15)

The weighting function for the Helmholtz free energy is computed from the difference between the weighting function for energy and entropy.

WAHD = WEHD − WSHD =

1 sHD − 2 2kB

(16)

The total entropy of a hard-disk fluid is the sum of excess entropy sex (T, ρ) and ideal gas

9



Page 10 of 44

entropy sIG (T, ρ).

sHD (T, ρ) = sIG (T, ρ) + sex (T, ρ)

(17)

The entropy of an ideal gas in two dimensions is given in eq(18). sIG (T, ρ) = 2 − ln kB

ρh2 2πmkB T

(18)

The expression for excess entropy for a hard-disk fluid 27 is given by eq(19). The expression for excess entropy sex (T, ρ) in eq(19) does not contain the spurious ln(Z) term that was pointed out in the three-dimensional 2PT paper 11 by Sun et al. 28 sex (T, ρ) =− kB

Z

y

(Z(y 0 ) − 1)d(ln y 0 )

(19)

0

2

y = (π/4)ρσ HD is the packing fraction of the hard-disk fluid. For a two-dimensional harddisk fluid, the compressibility 29 is given by the following expression.

Z(y) =

1 + y 2 /8 P = ρkB T (1 − y)2

(20)

Using eq(19) and eq(20), we arrive at the expression for the excess entropy of a hard-disk fluid. sex (T, ρ) 9 y 7 =− + ln(1 − y) kB 8 (1 − y) 8

10


(21)

Page 11 of 44 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


The weighting functions of a hard-disk fluid are given by the following set of equations.

(22a)

WSHD

(22b)

WAHD

2.4

1 2 1 ρh2 9 y 7 = 2 − ln − + ln(1 − y) 2 2πmkB T 8 (1 − y) 8 1 1 ρh2 9 y 7 = − 2 − ln − + ln(1 − y) 2 2 2πmkB T 8 (1 − y) 8

WEHD =

(22c)

Decomposition of DoS

Computing entropy and other thermodynamic properties of liquids is a difficult task. The 2PT(2D) model addresses this problem by decomposing the DoS of a liquid into two parts, the gas component, S g (ν), and the solid component, S s (ν) as shown in eq(23). By doing this, the method effectively treats some fraction of the degrees of freedom as gas-like and the remaining as solid-like. The solid component can be modeled using harmonic oscillators and the gas component using hard-disks. A liquid system is modeled as a combination of the two. S(ν) = S g (ν) + S s (ν)

(23)

Once the procedure to decompose the density of states is determined, thermodynamic properties of any system can be calculated using eq(24) because the weighting functions of a harmonic oscillator and a hard-disk fluid are already known to us. WPHO (ν) is the weighting function for harmonic oscillator and WPHD (ν) is the weighting function for hard-disk fluid, for a given thermodynamic property P . Z P =

∞

S

s

(ν)WPHO (ν)dν

0

Z +

∞

S g (ν)WPHD (ν)dν

(24)

0

The same approach of decomposing the density of states into a gas and solid component is used as in the three-dimensional 2PT method. The major difference between 2PT and 2PT(2D) is in the way gas-like degrees of freedom are modeled. The 2PT model uses hard-

11



Page 12 of 44

spheres to model the gas-like degrees of freedom, whereas the 2PT(2D) model uses harddisks. This difference arises because the 2PT(2D) model is used to compute thermodynamic properties of two-dimensional systems. As explained in Section 2.2, integrating S(ν) from 0 to ∞ gives us the number of degrees of freedom 2N . Similarly, integrating S g (ν) over the same range would give us 2Ng . Here, Ng is the number of particles that are effectively treated as a gas. Using eq(5), we arrive at the conclusion that the number of particles effectively treated as solid would be N − Ng . ∞

Z

S g (ν) = 2Ng

Z0 ∞

S s (ν) = 2(N − Ng )

(25)

0

The availability of analytical expressions for various thermodynamic properties of a harddisk fluid was the motivating factor for choosing it as a model system for the gas phase in 2PT(2D). The DoS of a hard-disk fluid can be calculated analytically. The per particle veHD locity autocorrelation function cHD (t) jk (t) and the total velocity autocorrelation function C

for hard-disks in two dimensions are given by the equations below. kB T exp(−αt) m Ng X 2 X HD mcHD C (t) = jk (t) cHD jk (t) =

(26a) (26b)

j=1 k=1

In the above equation T is the temperature, m is the mass of one particle, kB is the Boltzmann constant and α is the Enskog friction constant. cHD jk (t) is the velocity autocorrelation function of the j th hard-disk in the k th direction. The DoS of a hard-disk fluid is used for the gas part of the DoS and is derived from the

12


Page 13 of 44 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


cosine Fourier transform of the total velocity autocorrelation for hard-disks. 4 S (ν) = kB T

Z

4 = kB T

Z

g

∞

C HD (t)cos(2πνt)dt

0 Ng X 2 ∞X

0

Z

mcHD jk (t)cos(2πνt)dt

j=1 k=1 ∞

4 2N g exp(−αt)kB T cos(2πνt)dt kB T 0 8N g α = 2 α + 4π 2 ν 2 =

(27)

The zero frequency value of the density of states S(0) is related to the diffusion coefficient as indicated in eq(8). Thus S(0) = 0 when the diffusion coefficient is zero. Since solids don’t diffuse, S s (0) = 0 for a solid. Therefore, the zero frequency value of the DoS of the system s0 will be equal to S g (0). This gives us the condition, S g (0) = s0 . Using this condition, we can eliminate α from eq(27). S g (ν) = 1+

s0 πs0 ν 2

(28)

4N g

We still need to find the number of particles effectively treated as gas N g , if we wish to calculate S g (ν). Once we know S g (ν), we can calculate S s (ν) using the eq(23), as we can obtain S(ν) from simulation. Fluidicity factor f is introduced in Section 2.5, in order to determine the number of particles effectively treated as gas N g .

2.5

Fluidicity factor

The fluidicity factor f is the fraction of particles in the system that are treated as a gas. Ng f= N

13


(29)


Page 14 of 44

The hard-disk density of states expressed in terms of fluidicity f is given by the equation below. S g (ν) = 1+

s0 πs ν 2

(30)

0

4f N

As we can see from eq(30), if the fluidicity factor is known, then all thermodynamic properties of the system can be computed. In the three-dimensional 2PT paper 11 fluidicity factor f is defined in a way such that it satisfies the following conditions. 1. f = 1 in the low density or high temperature limit, i.e. the system is treated completely as a hard-disk fluid. 2. f = 0 in the high density limit, i.e. the system is treated completely as a solid. We follow the same principle in two dimensions and define fluidicity f as the ratio of the diffusion coefficient of the system to the diffusion coefficient of a hard-disk fluid at the same temperature and density in the zero pressure limit.

f=

D(T, ρ) D0HD (T, ρ; σ HD )

(31)

The diffusion coefficient is a physically relevant parameter that differentiates solid from gas and hence a suitable quantity to define the fluidicity factor f . The above given definition for f works extremely well in three dimensions. The dependence of the diffusion coefficient on system size in two dimensions forces the fluidicity factor f also to become system size dependent. The total vibrational density of states, s0 and f (calculated from s0 and other system parameters) once specified, completely determine all thermodynamic properties of the system in this model. In the rest of this section we will discuss about how to obtain f from other known quantities. In eq(31) D is the diffusion constant of the system and D0HD is the diffusion constant of a hard-disk fluid in the zero pressure limit. The analytical expression for the diffusion coefficient 30 in the zero pressure limit in two

14


Page 15 of 44 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


dimensions from Chapman-Enskog theory is given by eq(32).

D0HD (T, ρ; σ HD )

1/2 1 kB T = 2ρσ HD πm

(32)

The diffusion coefficient 30 of a hard-disk fluid DHD in two dimensions is given by eq(33) , 2

where y = (π/4)ρσ HD is the packing fraction of the system.

DHD (T, f ρ) = D0HD (T, f ρ, σ HD )

(1 − f y)2 7 1 − 16 fy

(33)

The diffusion coefficient for the hard-disk fluid DHD can also be obtained by integrating the VACF as shown in eq(34a). These equations can together then be used to eliminate σ HD from our final expressions. s0 kB T 8mf N s0 kB T D(T, ρ) = 8mN

DHD (T, f ρ) =

(34a) (34b)

From eq(34a) and eq(34b), we get DHD (T, f ρ) = D(T, ρ)/f . From eq(32), we get D0HD (T, f ρ) = D0HD (T, ρ)/f . Combining eq (31) with these two results we get equation below. (y 2 +

7 y)f 2 + (−2y − 1)f + 1 = 0 16

(35)

Using eq (31), eq(32) and eq(34b) we get the following set of equations. The parameter ∆ (normalized diffusivity) depends only on the material properties of the system.

f = ∆(T, ρ, m, s0 )y 1/2

(36)

1/2 s0 ρ1/2 kB T ∆(T, ρ, m, s0 ) = 2N m

(37)

15



Expressing eq(35) in terms of ∆ and f we obtain the following equation for f .

16f 6 ∆−4 + 7f 4 ∆−2 − 32f 3 ∆−2 − 16f + 16 = 0

(38)

The above equation can be solved for the fluidicity factor f once we know ∆. Normalized diffusivity ∆ can be calculated from information available to us from simulations. The reference energy V0 that was discussed in Section 2.3, is calculated for a system in the 2PT(2D) model by equating the energy obtained from MD simulations E M D with the sum of the energy of the classical harmonic oscillator for the solid degrees of freedom and energy of the hard-disk fluid for the gas degrees of freedom.

V0 = E

MD

−β

−1

Z

∞

h

i S s (ν)WEHO (ν) + S g (ν)WEHD (ν) dν

0

= E M D − β −1 2N (1 − 0.5f )

(39)

Normalized diffusivity ∆ quickly converges to the finite bulk value as system size is increased 1

0.8

Fluidicity f

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 16 of 44

0.6

0.4

0.2

0 1e-06

1e-04

1e-02

1e+00

1e+02

Normalized Diffusivity ∆

Figure 1: Dependence of fluidicity f on normalized diffusivity ∆ for 2PT(2D) in three dimensions. The system size dependence of normalized diffusivity is more complex in two dimensions. The convergence of normalized diffusivity in two dimensions is much 16


Page 17 of 44 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


slower. So fluidicity f is weakly system size dependent. Thus we would expect that even intensive thermodynamic quantities we calculate might be system size dependent. The fluidicity f in Fig 1 is obtained by solving eq(38) for different values of ∆. When solving eq(38) for fluidicity f , we find that there are two solutions for each ∆. One solution satisfies the boundary conditions that f →1 as ∆ → ∞ and the other one diverges with increasing ∆. We plotted the solution that satisfies the boundary conditions in Fig 1 and we use this solution for all calculations. In Table 2, 1PT(Q) refers to thermodynamic quantities computed by treating the entire system as a solid (S(ν) = S s (ν)) and using the quantum weighting functions of the harmonic oscillator (eq(12a), eq(12b) and eq(12c)). 2PT(Q) refers to thermodynamic quantities computed by splitting the density of states into solid and gas part (S(ν) = S s (ν)+S g (ν)). Quantum weighting functions of the harmonic oscillator are used for the solid part and hard-disk weighting functions are used for the gas part. 2PT(C) refers to thermodynamic quantities computed by splitting the density of states into solid and gas parts (S(ν) = S s (ν) + S g (ν)). Classical weighting functions of the harmonic oscillator are used for the solid part and harddisk weighting functions are used for the gas part.

3

Computational details

The 2PT(2D) model was validated by comparing the results produced by it for a two dimensional Lennard-Jones fluid with the equation of state predictions. Particles in this system interact with a pair-wise 12-6 Lennard-Jones potential given by the equation below. σ 12 σ 6 − V (r) = 4 r r

(40)

In the above equation, r is the distance between the particles. is a parameter associated with the depth of the potential well and σ is diameter of the particles. The parameters of Argon ( = 0.238kcal/mol, σ = 3.405˚ A and mass m=39.948 g/mol) were used for MD 17



1.1

2D Lennard-Jones Phase Diagram

1

m

0.9 Temperature T*

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 18 of 44

0.8 0.7

s

0.6 0.5 g

0.4

ll

u

s

0.3 0

0.2

0.4

0.6

0.8

1

ρ*

Figure 2: The phase diagram of a two-dimensional Lennard-Jones system (the phase boundaries have been obtained by fitting the data taken from Barker et al., 31 Smit and Frenkel, 32 Singh et al. 33 ). All the thermodynamic quantities presented in the paper have been calculated for the state points marked with open circles. These state points have been marked as g-gas; s-solid; l-liquid; m-metastable; u-unstable. The state points in the super critical region have not been marked for clarity. The solid lines represent the phase boundaries. simulations. But the final results were reported in reduced units. In reduced units: distance r∗ = r/σ, number density ρ∗ = ρσ 2 , mass m∗ = M/m, temperature T ∗ = kB T /, pressure P ∗ = P σ 2 /, energy E ∗ = E/, entropy S ∗ = S/kB , Helmholtz free energy A∗ = A/, Gibbs free energy G∗ = G/, Diffusivity D∗ = D(m/)1/2 /σ, Planck’s constant h∗ = h/((m)1/2 σ). We consider a range of 5 densities and 3 temperatures (Fig 2) including:

1) 2 stable solid phases; 2) 2 stable liquid phases; 3) 1 stable gas phase; 4) 8 super critical fluid phases; 5) 1 metastable solid phase; 6) 1 unstable fluid phase; The MD simulations were performed using LAMMPS. 34 They were performed in the

18


Page 19 of 44 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Table 1: Normalized Diffusivity ∆ and Fluidicity Factor f Calculated for the 15 State Points Discussed in This Study

ρ∗ 0.01 0.30 0.70 0.77 0.90

ρ∗ 0.01 0.30 0.70 0.77 0.90 a

Normalized Diffusivity ∆ T* 0.45 0.70 1.00 7.452 8.179 9.040 1.055a 1.810 1.866 0.356 0.613 0.620 0.209 0.387 0.438 0.001 0.001 0.005b Fluidicity f T* 0.45 0.70 0.974 0.978 0.644a 0.780 0.383 0.506 0.286 0.401 0.000 0.000

1.00 0.982 0.788 0.509 0.427 0.000b

Unstable state; b Metastable state.

constant number of particles, volume and temperature (NVT) ensemble. The Nose-Hoover thermostat was used to maintain the temperature (time constant of 0.05 ps). The simulations were performed for a system of 2500 particles. The lj/charmmf sw/coul/charmmf sh pair style was used to perform the simulations. 35 36 The parameters of Argon ( = 0.238kcal/mol, σ = 3.405˚ A and mass m=39.948 g/mol) were used. Inner cut-off, rin = 14.5˚ A and outer cut-off, rout = 15˚ A were used in the force switching scheme. The trajectory length for the 2PT(2D) results in Table 2 and associated figures is 80,000 MD steps (∆t = 8f s) or 640 ps. The normalized diffusivity ∆ is calculated from s0 and other state parameters like temperature. From ∆, the fluidicity factor f is obtained by solving eq(38). The normalized diffusivity ∆ and fluidicity factor f for all the 15 systems have been listed in Table 1.

19



4

Results and Discussions

4.1

Pressure, energy and quantum effects 6

m

a) Pressure

MD

5

Pressure P*

4

T*=1.0

s

3 s

2 1 0

g

u

0

0.2

l T*=0.45

l 0.4

0.6

0.8

1

ρ* 1.5 b) Total Energy

MD 2PT(Q)

1 0.5 Total Energy E*

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 20 of 44

g

0 T*=1.00

-0.5 -1 -1.5

m

u T*=0.45

-2

l

l

s s

-2.5 -3 0

0.2

0.4

0.6

0.8

1

ρ*

Figure 3: The pressure and energy(per particle) for a two-dimensional Lennard-Jones system. The curves are based on values obtained from the equation of state. (dotted, T*=0.45; dashed, T*=0.7; solid, T*=1.0). The open circles are from MD simulations and the solid circles are quantum corrected energies. The pressure and energy of two-dimensional Lennard-Jones systems were calculated using MD simulations for the 15 state points. These are compared with the values from the equation of state for 2D Lennard-Jones. The modified Benedict-Webb-Rubin equation of 20


Page 21 of 44 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Table 2: Comparison of Properties of Lennard-Jones Systems Calculated from Different Methods Pressure P* ρ∗ 0.01 0.30 0.70 0.77 0.90

Methoda MD MBWR EOS MD MBWR EOS MD MBWR EOS MD MBWR EOS MD MBWR EOS

0.45 0.004±0.000 0.004 0.008±0.001 -0.008 -0.025±0.002 -0.054 0.215±0.005 0.205 1.992±0.002 0.087

T* 0.70 0.007±0.000 0.007 0.119±0.001 0.117 0.748±0.002 0.727 1.480±0.003 1.456 3.806±0.006 3.869

1.00 0.010±0.000 0.010 0.268±0.001 0.263 1.696±0.002 1.677 2.810±0.003 2.800 5.866±0.011 7.168

T* 0.70 0.654 0.653±0.002 0.657 -0.318 -0.331±0.003 -0.332 -1.308 -1.338±0.000 -1.345 -1.516 -1.554±0.000 -1.567 -1.951 -2.016±0.001 -2.116

1.00 0.966 0.966±0.001 0.967 0.153 0.147±0.001 0.132 -0.889 -0.914±0.001 -0.927 -1.063 -1.096±0.001 -1.113 -1.417 -1.479±0.001 -1.355

Energy E* ρ∗ 0.01

0.30

0.70

0.77

0.90

a

Method 2PT(Q) MD MBWR EOS 2PT(Q) MD MBWR EOS 2PT(Q) MD MBWR EOS 2PT(Q) MD MBWR EOS 2PT(Q) MD MBWR EOS

0.45 0.360 0.365±0.005 0.382 -1.523 -1.559±0.007 -1.046 -1.717 -1.754±0.001 -1.771 -1.916 -1.963±0.001 -1.977 -2.383 -2.469±0.000 -2.437

21



Page 22 of 44

Table 2: (Continued) Entropy S* ρ∗ 0.01

0.30

0.70

0.77

0.90

a

Method 1PT(Q) 2PT(Q) 2PT(C) MBWR EOS 1PT(Q) 2PT(Q) 2PT(C) MBWR EOS 1PT(Q) 2PT(Q) 2PT(C) MBWR EOS 1PT(Q) 2PT(Q) 2PT(C) MBWR EOS 1PT(Q) 2PT(Q) 2PT(C) MBWR EOS

0.45 13.495±0.004 10.628±0.007 10.626±0.007 10.929 4.508±0.007 4.961±0.006 4.926±0.006 5.917 3.967±0.005 4.301±0.043 4.265±0.043 4.433 3.558±0.001 3.797±0.025 3.753±0.025 3.962 2.528±0.002 2.528±0.002 2.443±0.002 2.742

T* 0.70 14.496±0.009 11.316±0.026 11.316±0.026 11.418 7.139±0.008 6.841±0.004 6.835±0.004 7.220 4.817±0.002 5.160±0.019 5.143±0.019 5.198 4.319±0.003 4.673±0.030 4.651±0.030 4.689 3.251±0.001 3.251±0.001 3.208±0.001 3.212

1.00 15.141±0.016 11.721±0.001 11.721±0.001 11.787 7.984±0.004 7.514±0.003 7.512±0.003 7.778 5.393±0.006 5.688±0.044 5.679±0.044 5.696 4.896±0.001 5.228±0.003 5.214±0.003 5.228 3.839±0.002 3.839±0.002 3.813±0.002 4.122

Helmholtz free energy A* ρ∗ 0.01

0.30

Methoda 1PT(Q) 2PT(Q) 2PT(C) MBWR EOS 1PT(Q) 2PT(Q) 2PT(C) MBWR EOS

0.45 -5.291±0.002 -4.417±0.003 -4.418±0.003 -4.536 -3.382±0.003 -3.760±0.003 -3.776±0.003 -3.708

22

T* 0.70 -8.817±0.005 -7.268±0.018 -7.268±0.018 -7.335 -4.860±0.006 -5.111±0.003 -5.115±0.003 -5.386


1.00 -13.199±0.017 -10.755±0.001 -10.755±0.001 -10.820 -7.123±0.004 -7.363±0.004 -7.366±0.004 -7.646

Page 23 of 44 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Table 2: (Continued) 0.70

0.77

0.90

1PT(Q) 2PT(Q) 2PT(C) MBWR EOS 1PT(Q) 2PT(Q) 2PT(C) MBWR EOS 1PT(Q) 2PT(Q) 2PT(C) MBWR EOS

-3.337±0.003 -3.657±0.019 -3.673±0.019 -3.766 -3.397±0.002 -3.632±0.011 -3.653±0.011 -3.760 -3.528±0.001 -3.528±0.001 -3.568±0.001 -3.671

-4.386±0.002 -4.927±0.014 -4.939±0.014 -4.983 -4.299±0.002 -4.794±0.021 -4.810±0.021 -4.849 -4.231±0.001 -4.231±0.001 -4.262±0.001 -4.364

-5.839±0.002 -6.583±0.044 -6.593±0.044 -6.622 -5.585±0.002 -6.297±0.003 -6.310±0.003 -6.341 -5.267±0.003 -5.267±0.003 -5.293±0.003 -5.477

Gibbs free energy G* ρ∗ 0.01

0.30

0.70

0.77

0.90

a

Method 1PT(Q) 2PT(Q) 2PT(C) MBWR EOS 1PT(Q) 2PT(Q) 2PT(C) MBWR EOS 1PT(Q) 2PT(Q) 2PT(C) MBWR EOS 1PT(Q) 2PT(Q) 2PT(C) MBWR EOS 1PT(Q) 2PT(Q) 2PT(C) MBWR EOS

0.45 -4.875±0.004 -4.001±0.006 -4.002±0.006 -4.118 -3.356±0.006 -3.734±0.006 -3.751±0.006 -3.734 -3.373±0.006 -3.693±0.023 -3.710±0.023 -3.842 -3.118±0.008 -3.353±0.017 -3.374±0.017 -3.494 -1.315±0.003 -1.315±0.003 -1.355±0.003 -3.575 23

T* 0.70 -8.137±0.006 -6.588±0.020 -6.588±0.020 -6.654 -4.465±0.008 -4.716±0.005 -4.720±0.005 -4.996 -3.317±0.004 -3.858±0.016 -3.870±0.016 -3.944 -2.376±0.005 -2.872±0.024 -2.888±0.024 -2.957 -0.002±0.008 -0.002±0.008 -0.033±0.008 -0.065


1.00 -12.210±0.017 -9.766±0.002 -9.766±0.002 -9.830 -6.230±0.006 -6.470±0.005 -6.473±0.005 -6.770 -3.416±0.005 -4.160±0.047 -4.169±0.047 -4.227 -1.936±0.005 -2.648±0.006 -2.661±0.006 -2.706 1.251±0.015 1.251±0.015 1.225±0.015 2.488


Page 24 of 44

Table 2: (Continued) Self-diffusion coefficient D* ρ∗ 0.01 0.30 0.70 0.77 0.90

Methoda S(0) MSD S(0) MSD S(0) MSD S(0) MSD S(0) MSD

0.45 12.498 13.933 0.323 0.302 0.071 0.085 0.040 0.038 0.000

a

T* 0.70 17.109 18.765 0.691 0.705 0.153 0.162 0.092 0.086 0.000

1.00 22.601 24.452 0.852 0.881 0.185 0.215 0.125 0.122 0.000

MD, molecular dynamics results; MBWR EOS, the modified Benedict-Webb-Rubin equation of state proposed by Reddy and O’Shea 24 for a two-dimensional Lennard-Jones fluid; 2PT(Q), thermodynamic quantities calculated by decomposing the DoS as discussed in this paper and using quantum harmonic oscillator statistics for the solid part; 2PT(C), same as previous but using classical harmonic oscillator statistics for the solid part; 1PT(Q), thermodynamic properties calculated using harmonic approximation to DoS; S(0), diffusion coefficient calculated using eq(34b); MSD, diffusion coefficient calculated from the mean square displacement.

24


Page 25 of 44

state proposed by Reddy and O’Shea 24 is used for comparison in the fluid phase. As solid state points were not included while parametrizing this equation of state, we cannot compare the results for the density ρ∗ = 0.9, which lies in the solid phase. The total energy and pressure obtained from MD simulations seem to be in good agreement with the equation of state except for one state point. The total energy for this state point (T*=0.45, ρ∗ = 0.3) deviates significantly from the equation of state. As this state point lies in the unstable region of the 2D Lennard-Jones phase diagram, the deviation from the EOS was expected. One might notice from Fig 3 that the classical energies obtained from MD simulations fit more closely to the EOS than the quantum corrected energies obtained from 2PT(2D). This is because the EOS was parametrized by fitting to classical MD and MC simulation data for 2D Lennard-Jones fluid. The quantum correction to energy is about 1.5% in the gas phase, 2.5% in the liquid phase and 3.5% in the solid phase.

4.2

Velocity autocorrelation and density of state distribution 250 Velocity autocorrelation function (*1000)

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


T*=0.45

200

ρ*=0.01 (gas)

150 100 50

ρ*=0.7 (liquid) ρ*=0.3 (unstable)

0 ρ*=0.77 (liquid)

-50

ρ*=0.9 (solid)

-100 0

0.5

1

1.5

2

2.5

3

Time (ps)

Figure 4: The mass weighted velocity autocorrelation function (VACF) for a two-dimensional Lennard-Jones system at T*=0.45 and densities ρ∗=0.01,0.3,0.7,0.77,0.9. The problem of long-time tails in velocity autocorrelation functions of two-dimensional fluids is well known. It generally decays as t−1 at long times 37 38 and its integral does not 25



Page 26 of 44

converge. Thus the diffusion constant in eq(4) would be ill-defined. However, analysis of VACF of the system studied in this paper (2500 particles interacting with a Lennard-Jones potential) showed that at long times, it decays faster than t−1 and the integral in eq(4) converges. The reason for faster than t−1 decay could be because the interaction potential used was different from Alder and Wainwright 37 (hard-disks). The convergence of the integral in eq(4) can also be attributed in part to finite-size effects, since a small system size (N=2500) was used and at long times the resulting velocity autocorrelation function decays to zero within fluctuations. The plots analyzing long-time behavior of velocity autocorrelation function are included in the supporting information. Fig 4 shows the short-time behavior of velocity autocorrelation function. Table 6: Excess Degrees of Freedom Calculated for the 15 State Points Discussed in This Study. Excess Degrees of Freedom T* 0.45 0.70 1 (0.0%) 2 (0.0%) 30 (0.6%)a 33 (0.7%) 52 (1.0%) 58 (1.2%) 47 (0.9%) 56 (1.1%) 0 (0.0%) 0 (0.0%)

ρ∗ 0.01 0.30 0.70 0.77 0.90 a

1.00 2 (0.0%) 38 (0.8%) 60 (1.2%) 59 (1.2%) 0 (0.0%)b

Unstable state; b Metastable state.

Often DoS of the gas component (hard-disk fluid) exceeds total DoS at high frequencies and the difference between them is called excess DoS. The results reported in this paper are computed by ignoring excess DoS of the gas component. The excess degrees of freedom in Table 6 were calculated by integrating total DoS (only over high frequencies > 80 cm−1 ) and subtracting it from the integral of gas-component of DoS (only over high frequencies > 80 cm−1 ). The integral was performed only over frequencies above 80 cm−1 so as to avoid issues of convergence of DoS at low frequencies because of long-time tail of VACF. The excess degrees of freedom are less than 2.0% of the total degrees of freedom as shown in Table 6. The memory function approach of Desjarlais 39 can be used to correct for error due to excess 26


Page 27 of 44 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


DoS. When the system is in gas state the VACF decreases monotonically. At such low densities the collision frequency is very low and thus the velocity of a particle remains correlated for a longer period of time and VACF decays slowly. The density of states (DoS) in the gas state has been plotted in Fig 5. From the DoS of the gas state we can see that the fluidicity factor f is close to 1. This says that a Lennard-Jones system at high temperature and low density has thermodynamic properties very close to that of a hard-disk fluid. As explained in the three-dimensional 2PT paper, 11 at higher temperatures the repulsive part of the potential plays a more dominant role. The state point in the solid phase (T*=0.45 ρ∗ = 0.9) has a VACF that quickly falls and then oscillates about zero. This is characteristic of crystal ordering. As seen from Fig 5, the DoS of the solid state points is zero at zero frequency. A solid crystal does not have diffusive modes and thus one would expect S(0) to be zero. The fluidicity factor f is nearly zero in the solid state. The DoS of the liquid state point (T*=0.45 ρ∗ = 0.77) has characteristics of both solid and gas. As explained in the three-dimensional 2PT paper, 11 the similarity of the DoS of the liquid state to both the solid and the gas was the motivation behind the two-phase thermodynamic model. Fig 5 shows the decomposition of S(ν) into the fraction that is treated as a hard-disk fluid S g (ν) and the fraction that is treated as a solid S s (ν).

4.3

Entropy, free energies and quantum effects

The entropy calculated from the 1PT(Q), 2PT(C) and 2PT(Q) methods are compared with the equation of state in Fig 6. Apart from the unstable point, the entropy is calculated fairly accurately (about 5-6% error) using the 2PT(Q) method. The 2PT(Q) method overestimates the entropy in the liquid phase. The entropy is accurately described in the gas phase. We cannot compare the solid phase entropy with the equation of state as solid state points were not included while parameterizing the equation. The entropy calculated from the 2PT(Q) 27



400 a) Gas (T*=0.45 ρ*=0.01)

350 300

S(ν) (m)

250 200 150

total gas solid

100 50 0 0

0.2

0.4

0.6

0.8

1

1.2

1.4

-1

ν (cm ) 140 b) Liquid (T*=0.45 ρ*=0.77)

120 total

S(ν) (cm)

100 solid

80 60 40 20

gas

0 0

20

40

60

80

100

ν (cm-1) 160 c) Solid (T*=0.45 ρ*=0.9)

140 120 S(ν) (cm)

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 28 of 44

100

total solid gas

80 60 40 20 0 0

20

40

60

80

100

ν (cm-1)

Figure 5: The density of states of a 2D LJ system and its decomposition for (a) gas (T*=0.45 ρ∗ = 0.01) ,(b) liquid (T*=0.45 ρ∗ = 0.77 and (c) solid (T*=0.45 ρ∗ = 0.9). The total density of states S(ν) is given by the solid line, the gas component S g (ν) is given by the dashed line and the solid component S s (ν) is given by the dotted line. 28 ACS Paragon Plus Environment

Page 29 of 44

16

a) Entropy - 1PT

1PT(Q)

14

g

Entropy S*

12 10 8 6 u

4

l

l

2 0

0.2

0.4

0.6

m s s

0.8

1

ρ*

16

b) Entropy - 2PT

2PT(Q) 2PT(C)

14 12 Entropy S*

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


g 10 8 6 u

l

4

l

2 0

0.2

0.4

0.6

0.8

m s s 1

ρ*

Figure 6: Entropy(per particle) determined using a) 1PT(Q), one-phase quantum model b) 2PT(C), two-phase classical model and 2PT(Q), two-phase quantum model. The curves are based on the equation of state predictions for a 2D LJ system. (dotted, T*=0.45; dashed, T*=0.7; solid, T*=1.0).

29



-2

a) Helmholtz Free Energy - 1PT

Helmholtz Free Energy A*

l

u

-4

l

g -6

s s m

-8 -10 -12 -14 1PT(Q) -16 0

0.2

0.4

0.6

0.8

1

ρ*

-2

b) Helmholtz Free Energy - 2PT

-4 Helmholtz Free Energy A*

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 30 of 44

u

g

l

l

-6

s s m

-8 -10 -12 2PT(Q) 2PT(C)

-14 -16 0

0.2

0.4

0.6

0.8

1

ρ*

Figure 7: Helmholtz free energy(per particle) determined using a) 1PT(Q), one-phase quantum model b) 2PT(C), two-phase classical model and 2PT(Q), two-phase quantum model. The curves are based on the equation of state predictions for a 2D LJ system. (dotted, T*=0.45; dashed, T*=0.7; solid, T*=1.0).

30


Page 31 of 44

2 m

a) Gibbs Free Energy - 1PT

Gibbs Free Energy G*

0

s s

-2 u

-4

l

l

g -6 -8 -10 -12

1PT(Q)

-14 0

0.2

0.4

0.6

0.8

1

ρ*

2 m

b) Gibbs Free Energy - 2PT

0 Gibbs Free Energy G*

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


s s

-2 -4

u

g

l

l

-6 -8 -10 2PT(Q) 2PT(C)

-12 -14 0

0.2

0.4

0.6

0.8

1

ρ*

Figure 8: Gibbs free energy(per particle) determined using a) 1PT(Q), one-phase quantum model b) 2PT(C), two-phase classical model and 2PT(Q), two-phase quantum model. The curves are based on the equation of state predictions for a 2D LJ system. (dotted, T*=0.45; dashed, T*=0.7; solid, T*=1.0).

31



and 2PT(C) are almost equal, so we can conclude that the quantum corrections for a LJ fluid are quite small. Quantum corrections are expected to be larger for larger molecules with covalent bonds. Just like in the 3D case, we note that the 1PT(Q) method overestimates the entropy in the gas phase and underestimates the entropy in the liquid phase. The weighting function of the harmonic oscillator is much larger than that of a hard-disk fluid at low frequencies and this leads to overestimation of entropy in the gas phase. The Helmholtz free energy calculated from all the three methods mentioned above are plotted in Fig 7. The Helmholtz free energy is overestimated in the gas phase (by about 2%) using the 2PT(Q) method. It is underestimated by less than 1% in the liquid phase by the 2PT(Q) method. The 1PT(Q) method overestimates the Helmholtz free energy in the liquid phase and underestimates it in the gas phase. The 2PT(2D) method is far more accurate in estimating the Helmholtz free energy than the one-phase method. Fig 8 compares the Gibbs free energy calculated by all the three methods. The 2PT(Q) data is in good agreement with the equation of state. The value of Gibbs free energy is quite accurately described in the liquid phase (within 0.5% error) using 2PT(2D) method. The Gibbs free energy is overestimated in the gas phase by about 2.5% using 2PT(2D) method. The results from the 1PT model are quite inaccurate when compared with the 2PT(2D) method. In the solid phase, both the 1PT(Q) and 2PT(Q) give almost exactly the same results because the fluidicity f is nearly zero and both the methods compute thermodynamic properties using only the quantum harmonic oscillator weighting functions. The entropy, Helmholtz free energy and Gibbs free energy computed for all the 15 state points are given in Table 2.

4.4

Self-diffusion coefficient

The diffusion coefficient is calculated from the zero frequency value of the density of states S(0) using eq(34b) in the 2PT(2D) method for a 640 ps run. The diffusion coefficient is also calculated from the mean square displacement obtained from a MD run. The mean square 32


Page 32 of 44

Page 33 of 44

100 2PT MSD g

10 Self-Diffusivity D*

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


1 u 0.1

l l

0.01 -0.2

0

0.2

0.4 ρ*

0.6

0.8

1

Figure 9: Reduced diffusivity calculated from the zero frequency value of the DoS (2PT(2D)) (640 ps trajectory) is compared with the diffusivity obtained from mean-square displacement calculations (MSD) (by linearly fitting MSD vs time between 540 ps and 740 ps). displacement vs time plot is fit linearly between 540 ps and 740 ps and the slope of this linear fit is used to obtain the diffusion coefficient. The diffusion coefficient obtained by these two methods are compared in Fig 9. The diffusion coefficient data is presented in Table 2.

4.5

Convergence efficiency

The 2PT method in 3D shows convergence to within 0.2% in 2500 MD steps (20 ps) in the gas phase and 1.5% in the liquid phase. The 2PT(2D) method requires longer simulations for convergence. In this method, the entropy converges to within 2% by 5000 MD steps (40ps) in the gas phase and to within 5% by 5000 MD steps (40ps) in the liquid phase. Fig 10 plots the entropy calculated by the 2PT(2D) method for different run lengths. The data for convergence of entropy is provided in Table 3.

4.6

System size dependence

The self-diffusion coefficient in two-dimensional systems are dependent on system size 25 and diverge as the system size increases as can be seen in Fig 11. Since the 2PT(2D) method 33



Table 3: Convergence of Entropy(Per Particle) from 2PT(2D) Method for Liquid (ρ*=0.77, T*=0.45) and Gas (ρ*=0.01, T*=0.45) MD steps time (ps) 2PT(Q) deva 2PT(C) deva MBWR EOSb

500 4 7.564 0.010 7.543 0.010

1000 2500 5000 10000 20000 40000 80000 8 20 40 80 160 320 640 ∗ ρ =0.01, T*=0.45 (gas) 7.656 9.822 10.412 10.617 10.653 10.634 10.628 0.004 0.012 0.011 0.016 0.014 0.013 0.007 7.647 9.820 10.411 10.616 10.652 10.632 10.626 0.004 0.012 0.011 0.016 0.014 0.013 0.007 10.929

ρ∗ =0.77, T*=0.45 (liquid) 2PT(Q) 3.700 3.720 3.738 3.819 3.767 3.836 a dev 0.018 0.020 0.014 0.028 0.014 0.039 2PT(C) 3.652 3.674 3.693 3.774 3.723 3.792 deva 0.018 0.020 0.015 0.028 0.014 0.039 MBWR EOSb 3.962 a Standard deviation. Three samples were used for each run length. b Modified Benedict-Webb-Rubin equation of Reddy and O’Shea 24

3.768 0.015 3.724 0.015

3.797 0.025 3.753 0.025

12 Convergence of entropy: 2PT

11 10

2PT(Q) 2PT(C) MBWR EOS

9 Entropy S*

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 34 of 44

8 gas (ρ*=0.01 T*=0.45) 7 6

liquid (ρ*=0.77 T*=0.45)

5 4 3 10

100

1000

Time (ps)

Figure 10: Convergence of entropy(per particle) with length of MD run for 2PT(2D) method. Integration time step of 8fs was used for MD simulations. The solid lines represent the MBWR EOS results for entropy for the liquid (ρ*=0.77, T*=0.45) and gas (ρ*=0.01, T*=0.45) phases.

34


Page 35 of 44

depends upon the diffusion coefficient for calculating entropy, we analyze the system size dependence of the entropy calculated by this method. We do see a gradual increase in entropy with system size from Fig 12. The results in Fig 11 and Fig 12 were obtained from 640 ps MD runs in the liquid phase (ρ*=0.77 T*=0.45). It should be noted that the size dependence in entropy observed here is a result of using the Chapman-Enskog theory eq(32) for the hard-disk diffusivity. This theory does not address the fact that diffusivity in two-dimensional systems diverges with system size. However, we are not aware of other analytical models that can address this issue properly. While the method proposed here suffers from system size dependence, our results show that the 2PT(2D) does provide good temperature and density dependence of entropy in two dimensions for a fixed system size. 0.055 System size dependence

0.05

0.045 Diffusivity D*

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


0.04

0.035

0.03

0.025 0

20

40

60

80

100

120

140

Box length L*

Figure 11: Diffusion coefficient calculated using 2PT(2D) method for different system sizes for 2D LJ fluid in liquid (ρ*=0.77 T*=0.45) phase.

4.7

Chemical potential

The Gibbs free energy G*(N,V,T) was plotted as a function of number of particles N in Figures 13,14 and the chemical potential µ∗tot was obtained from the slope of the plots. ∗ The part of the chemical potential due to the momentum, T ∗ ln λ∗2 T , is subtracted from µtot

to obtain the reduced chemical potential µ∗ as given in eq(41). 40 The chemical potential is 35



4 System size dependence

3.95

3.9 Entropy S*

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 36 of 44

3.85

3.8

3.75

3.7 0

20

40

60 80 100 Box length L*

120

140

Figure 12: Entropy(per particle) calculated using 2PT(2D) method for different system sizes for 2D LJ fluid in liquid (ρ*=0.77 T*=0.45) phase. The points are results from 2PT(2D) and the dashed line is the MBWR EOS value of entropy. The percentage error in the 2PT(2D) entropy decreases from 20% (N=225 and L∗ = 17.09) to 1% (N=10,000 and L∗ = 113.96) with increase in system size. compared with the value of chemical potential calculated using the Gibbs ensemble technique by Smit and Frenkel, 32 and also with the MBWR equation of state 24 in Table 4.

µ∗ = µ∗tot − T ∗ ln λ∗2 T λ∗T = √

h∗ 2πm∗ T ∗

(41) (42)

Table 4: Chemical Potential from 2PT(2D), MBWR Equation of State, Smit and Frenkel 32 ρ∗ T∗ µ∗ (2PT(2D)) µ∗ (MBWR)a µ∗ (Smit et al )b 0.030 0.45 -1.90 -1.76 -1.77 0.722 0.45 -1.66 -1.79 -1.82 a Modified Benedict-Webb-Rubin equation of Reddy and O’Shea 24 b Chemical potential taken from Smit and Frenkel 32

36


Page 37 of 44

Chemical potential at T*=0.45 rho*=0.03 Linear fit

T*=0.45 ρ*=0.03

Gibbs free energy G*

0

-10000

-20000

-30000

slope=-3.8841 +/- 0.0030

-40000 0

2000

4000

6000

8000

10000 12000

Number of particles N

Figure 13: Gibbs free energy calculated using 2PT(2D) method for different system sizes for 2D LJ fluid in gas (ρ*=0.03 T*=0.45) phase.

Chemical potential at T*=0.45 rho*=0.722 Linear fit

T*=0.45 ρ*=0.722

0 Gibbs free energy G*

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


-10000

-20000

-30000

slope=-3.6432 +/- 0.0034

-40000 0

2000

4000

6000

8000

10000 12000

Number of particles N

Figure 14: Gibbs free energy calculated using 2PT(2D) method for different system sizes for 2D LJ fluid in liquid (ρ*=0.722 T*=0.45) phase.

37



5

Conclusions and Outlook

A new method, 2PT(2D), has been developed for computing the entropy and free energy of two-dimensional systems. This method has been tested by computing the thermodynamic properties of a two-dimensional Lennard-Jones fluid over several densities and temperatures and comparing with the MBWR equation of state. The computed entropy and free energy are in excellent agreement with those obtained from the EOS. The chemical potential has also been computed using the 2PT(2D) method for two state points and compares well with the existing values in the literature. The 2PT(2D) method is computationally very efficient and gives converged entropy values using 30-50 ps long MD trajectories. The 2PT(2D) method can also be extended to molecular systems in a way similar to that proposed for the 2PT method by Lin et al. 21 We recognize that the diffusion coefficient of two-dimensional systems are system size dependent and that this might affect the thermodynamic properties predicted by this method. The variation of thermodynamic quantities calculated using the 2PT(2D) method with system size has also been analyzed. We believe this method will help address wide range of problem in two dimensions including confined fluid in slit-pore geometry of various 2D materials, fluid on lipid surface and on 2D DNA lattice. Future works will involve demonstrating entropy and free energy calculation in some of these systems.

Acknowledgement We would like to thank Prof. Chandan Dasgupta, Prof. Biman Bagchi and Prof. Shankar Prasad Das for useful discussions and comments. We thank DAE, India for financial support.

38


Page 38 of 44

Page 39 of 44 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Supporting Information Available Figures showing the convergence of entropy with simulation length for different temperatures and long-time behavior of velocity autocorrelation function.

References (1) Kirkwood, J. G. Statistical Mechanics of Fluid Mixtures. The Journal of Chemical Physics 1935, 3, 300–313. (2) Zwanzig, R. W. High-Temperature Equation of State by a Perturbation Method. I. Nonpolar Gases. The Journal of Chemical Physics 1954, 22, 1420–1426. (3) Widom, B. Some Topics in the Theory of Fluids. The Journal of Chemical Physics 1963, 39, 2808–2812. (4) Torrie, G. M.; Valleau, J. P. Nonphysical Sampling Distributions in Monte Carlo FreeEnergy Estimation: Umbrella Sampling. Journal of Computational Physics 1977, 23, 187–199. (5) Frenkel, D.; Smit, B. Understanding Molecular Simulation: From Algorithms to Applications, 2nd ed.; Academic Press: San Diego, 2001; pp 167–196. (6) Karplus, M.; Kushick, J. N. Method for Estimating the Configurational Entropy of Macromolecules. Macromolecules 1981, 14, 325–332. (7) Reinhard, F.; Grubm¨ uller, H. Estimation of Absolute Solvent and Solvation Shell Entropies via Permutation Reduction. The Journal of chemical physics 2007, 126, 014102. (8) Sasikala, W. D.; Mukherjee, A. Single Water Entropy: Hydrophobic Crossover and Application to Drug Binding. The Journal of Physical Chemistry B 2014, 118, 10553– 10564. 39



(9) Henchman, R. H. Free Energy of Liquid Water from a Computer Simulation via Cell Theory. The Journal of chemical physics 2007, 126, 064504. (10) Gerogiokas, G.; Calabro, G.; Henchman, R. H.; Southey, M. W.; Law, R. J.; Michel, J. Prediction of Small Molecule Hydration Thermodynamics with Grid Cell Theory. Journal of chemical theory and computation 2013, 10, 35–48. (11) Lin, S.-T.; Blanco, M.; Goddard III, W. A. The Two-Phase Model for Calculating Thermodynamic Properties of Liquids from Molecular Dynamics: Validation for the Phase Diagram of Lennard-Jones Fluids. The Journal of chemical physics 2003, 119, 11792–11805. (12) Geim, A. K.; Novoselov, K. S. Nanoscience and Technology: A Collection of Reviews from Nature Journals; World Scientific, 2010; pp 11–19. (13) Geim, A. K.; Grigorieva, I. V. Van Der Waals Heterostructures. Nature 2013, 499, 419. (14) Chen, K.; Kaplan, T.; Mostoller, M. Melting in Two-Dimensional Lennard-Jones Systems: Observation of a Metastable Hexatic Phase. Physical review letters 1995, 74, 4019. (15) Raghav, N.; Chakraborty, S.; Maiti, P. K. Molecular Mechanism of Water Permeation in a Helium Impermeable Graphene and Graphene Oxide Membrane. Physical Chemistry Chemical Physics 2015, 17, 20557–20562. (16) Kumar, H.; Dasgupta, C.; Maiti, P. K. Phase Transition in Monolayer Water Confined in Janus Nanopore. Langmuir 2018, 34, 12199–12205, PMID: 30216072. (17) Srivastava, D.; Santiso, E. E.; Gubbins, K. E. Pressure Enhancement in Confined Fluids: Effect of Molecular Shape and Fluid–Wall Interactions. Langmuir 2017, 33, 11231– 11245.

40


Page 40 of 44

Page 41 of 44 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


(18) Vasu, K. S.; Prestat, E.; Abraham, J.; Dix, J.; Kashtiban, R. J.; Beheshtian, J.; Sloan, J.; Carbone, P.; Neek-Amal, M.; Haigh, S. J. e. a. Van Der Waals Pressure and Its Effect on Trapped Interlayer Molecules. Nature communications 2016, 7, ncomms12168. (19) Pascal, T. A.; Lin, S.-T.; Goddard III, W. A. Thermodynamics of Liquids: Standard Molar Entropies and Heat Capacities of Common Solvents from 2PT Molecular Dynamics. Physical chemistry chemical physics 2011, 13, 169–181. (20) 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. Journal of chemical theory and computation 2011, 8, 61–74. (21) Lin, S.-T.; Maiti, P. K.; Goddard III, W. A. Two-Phase Thermodynamic Model for Efficient and Accurate Absolute Entropy of Water from Molecular Dynamics Simulations. The Journal of Physical Chemistry B 2010, 114, 8191–8198. (22) Kumar, H.; Mukherjee, B.; Lin, S.-T.; Dasgupta, C.; Sood, A.; Maiti, P. K. Thermodynamics of Water Entry in Hydrophobic Channels of Carbon Nanotubes. The Journal of chemical physics 2011, 134, 124105. (23) Lin, S.-T.; Maiti, P. K.; Goddard, W. A. Dynamics and Thermodynamics of Water in PAMAM Dendrimers at Subnanosecond Time Scales. The Journal of Physical Chemistry B 2005, 109, 8663–8672. (24) Reddy, M. R.; O’Shea, S. F. The Equation of State of the Two-Dimensional LennardJones Fluid. Canadian Journal of Physics 1986, 64, 677–684. (25) Vogele, M.; Hummer, G. Divergent Diffusion Coefficients in Simulations of Fluids and Lipid Membranes. The Journal of Physical Chemistry B 2016, 120, 8722–8732. 41



(26) McQuarrie, D. Statistical mechanics; Harper’s chemistry series; Harper & Row, 1975. (27) Carnahan, N. F.; Starling, K. E. Thermodynamic Properties of a Rigid-Sphere Fluid. The Journal of Chemical Physics 1970, 53, 600–603. (28) Sun, T.; Xian, J.; Zhang, H.; Zhang, Z.; Zhang, Y. Two-phase thermodynamic model for computing entropies of liquids reanalyzed. The Journal of Chemical Physics 2017, 147, 194505. (29) Henderson, D. A Simple Equation of State for Hard Discs. Molecular Physics 1975, 30, 971–972. (30) Garc´ıa-Rojo, R.; Luding, S.; Brey, J. J. Transport Coefficients for Dense Hard-Disk Systems. Physical Review E 2006, 74, 061305. (31) Barker, J.; Henderson, D.; Abraham, F. F. Phase Diagram of the Two-Dimensional Lennard-Jones System; Evidence for First-Order Transitions. Physica A: Statistical Mechanics and its Applications 1981, 106, 226–238. (32) Smit, B.; Frenkel, D. Vapor–Liquid Equilibria of the Two-Dimensional Lennard-Jones Fluid(s). The Journal of chemical physics 1991, 94, 5663–5668. (33) Singh, R. R.; Pitzer, K. S.; de Pablo, J. J.; Prausnitz, J. M. Monte Carlo Simulation of Phase Equilibria for the Two-Dimensional Lennard-Jones Fluid in the Gibbs Ensemble. The Journal of chemical physics 1990, 92, 5463–5466. (34) Plimpton, S. Fast Parallel Algorithms for Short-Range Molecular Dynamics. Journal of computational physics 1995, 117, 1–19. (35) Steinbach, P. J.; Brooks, B. R. New Spherical-Cutoff Methods for Long-Range Forces in Macromolecular Simulation. Journal of Computational Chemistry 1994, 15, 667–683.

42


Page 42 of 44

Page 43 of 44 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


(36) Brooks, B.; Brooks III, C.; Mackerell Jr, A.; Nilsson, L.; Petrella, R.; Roux, B.; Won, Y.; Archontis, G.; Bartels, C.; Boresch, S. et al. HL Woodcock. J. Chem 2009, 30, 1545– 1614. (37) Alder, B. J.; Wainwright, T. E. Decay of the Velocity Autocorrelation Function. Phys. Rev. A 1970, 1, 18–21. (38) Ferrario, M.; Fionino, A.; Ciccotti, G. Long-time tails in two-dimensional fluids by molecular dynamics. Physica A: Statistical Mechanics and its Applications 1997, 240, 268 – 276, Proceedings of the Euroconference on the microscopic approach to complexity in non-equilibrium molecular simulations. (39) Desjarlais, M. P. First-principles calculation of entropy for liquid metals. Phys. Rev. E 2013, 88, 062145. (40) Smit, B.; Frenkel, D. Calculation of the chemical potential in the Gibbs ensemble. Molecular Physics 1989, 68, 951–958.

43



6

Table of Contents Image

Figure 15: For Table of Contents Use Only

44


Page 44 of 44

Efficient Computation of Entropy and Other Thermodynamic Properties

Recommend Documents