Effects of Global and Local Anharmonicities on the ... - ACS Publications

Subscriber access provided by UNIV NEW ORLEANS

Article

Effects of global and local anharmonicities on the thermodynamic properties of sulfuric acid monohydrate Lauri Johannes Partanen, Vesa Hänninen, and Lauri Halonen J. Chem. Theory Comput., Just Accepted Manuscript • DOI: 10.1021/acs.jctc.6b00683 • Publication Date (Web): 23 Sep 2016 Downloaded from http://pubs.acs.org on September 29, 2016

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 free 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 accessible to all readers and 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.

Journal of Chemical Theory and Computation 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 47

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

Journal of Chemical Theory and Computation

Effects of global and local anharmonicities on the thermodynamic properties of sulfuric acid monohydrate Lauri Partanen,∗ Vesa Hänninen, and Lauri Halonen Laboratory of Physical Chemistry, Department of Chemistry, P.O. Box 55 (A.I. Virtasen aukio 1), FIN-00014 University of Helsinki, Finland E-mail: [email protected]

Abstract We use state-of-the-art electronic structure calculation methods and large basis sets to obtain reliable values for the thermodynamic properties of sulfuric acid monohydrate and study the effects of vibrational anharmonicity on these properties. We distinguish between two forms of vibrational anharmonicity: local anharmonicity which refers to the anharmonicity of the vibrational modes of a given cluster conformer and global anharmonicity which originates from accounting for the presence of different conformers in the first place. In our most accurate approach, we solve the nuclear Schrödinger equation variationally for the intermolecular large-amplitude motions, thus quantum mechanically accounting for the presence of higher energy conformers for both reactants and products, while using standard vibrational perturbational approach for the other vibrational modes. This results in a value of -11.0 kJ/mol for the reaction Gibbs free energy at 298.15 K. When standard vibrational perturbational approaches are employed the effects of local anharmonicity depend heavily on the choice of the electronic structure calculation basis set. In fact, better results can often be achieved by

1 ACS Paragon Plus Environment


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

combining a simple harmonic treatment for the vibrational partition function with a statistical mechanical accounting of global anharmonicity. Thus, we recommend that future studies that intend to include anharmonicity start by accounting for the presence of higher energy conformers and only then consider if local anharmonicity calculations are feasible and necessary.

1

Introduction

The detailed mechanisms of cloud formation continue to elude the scientific community. Even though clouds fundamentally affect a variety of atmospheric phenomena, the predictive models of cloud formation have generally met with a lack of success. 1 One of the established facts in cloud formation is the central importance of sulfuric acid 2–4 as much of new particle formation in the atmosphere begins with the complexation of sulfuric acid and water with trace species such as ammonia, 5–9 amines, 10–15 ions, 16–24 and various organic species 25–34 in the gas phase. A significant portion of the sulfuric acid molecules in the atmosphere are contained in hydrates 35 with the exact amount changing with the altitude. 36 Because these hydrates stabilize the vapor, it is crucial that the thermodynamic properties of their formation are known accurately if reliable nucleation rates are desired. 37,38 Both ab initio and density functional theory calculations have been previously used to obtain structural and thermodynamical properties for small sulfuric acid hydrates. 36,39–43 The proper accounting of both local and global anharmonicities has been problematic in these and other ab initio studies on thermodynamic properties of atmospheric clusters, but according to Kathmann et al. 38 it is vital for obtaining accurate results at least in the case of aqueous ionic clusters. In this article, local anharmonicity refers to the anharmonicity of the vibrational modes of a given cluster conformer whereas global anharmonicity originates from accounting for the presence of different conformers in the first place. Because of the challenges in treating anharmonicity, the variance in the thermodynamic properties between even the two most recent studies on H2 O·H2 SO4 has been about 3.0 kJ/mol. 16,44 Most importantly, to the best of 2 ACS Paragon Plus Environment

Page 2 of 47

Page 3 of 47

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


our knowledge it still remains an open question how large of an impact global anharmonicity, i.e., the presence of multiple conformers, has on the thermodynamic properties of H2 O·H2 SO4 and larger clusters. This problem has been partly addressed in a recent publication 45 which explored the effect of the increasingly popular 16,24,44,46–61 but incorrect Boltzmann averaging formulae compared to the correct statistical mechanical method 62 for the incorporation of higher energy conformers. In this study, one of our aims is to compare the accuracy of the statistical mechanical method with a quantum mechanical accounting of multiple conformers for one of the most important atmospheric systems - the sulfuric acid monohydrate. With advances made in electronic stucture methods and computation of anharmonic vibrational frequencies, calculation of thermodynamic properties can nowadays be done accurately by utilizing standard quantum chemistry software. As the equilibrium constant is exponentially dependent on the reaction energy, which consists of the stoichiometrically weighted differences between the electronic energies and the zero point energies (ZPEs) of the reactants and products, 63 it is crucial that both of these properties are calculated with high precision. Accurate computation of interaction energies can be performed even for mediumsized systems, due to major developments in explicitly correlated methods that take into account the singularities in the potential energy at points where two electrons coincide. 64–68 Specifically, the explicitly correlated coupled-cluster singles and doubles method with perturbative triples (CCSD(T)-F12) has made it possible to obtain accurate electronic energies and equilibrium structures with reasonably small basis sets. This has substantially reduced computational time when compared with the standard CCSD(T) method. 68,69 Furthermore, when the CCSD(T)-F12 approach is combined with relatively large basis sets and complete basis set extrapolation techniques for the correlation energies, 70 it is possible to calculate interaction energies with high accuracies. 71 There are several popular methods that can be used to go beyond the harmonic approximation for ZPE and vibrational partition function calculations. Many of these approaches including vibrational second-order Møller-Plesset pertubation theory (VMP2), configuration



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

interaction (VCI), and coupled-cluster (VCC) start with the vibrational self-consistent field (VSCF) calculation where, analogously to the SCF methods in electronic structure calculations, the coupling between modes is treated in an average way. 72–74 While these methods yield good ZPEs and anharmonic frequencies for the fundamental modes, their nonlinear scaling with the number of normal modes and the number of modes being correlated can make them expensive to utilize. 75 A second problem is that VSCF often gives unreasonable anharmonic frequencies for large amplitude low-frequency modes. 76–78 Another popular family of methods is based on the vibrational second-order perturbation theory (VPT2), where the anharmonic corrections are obtained from higher order derivatives of the potential energy surface (PES) along the normal mode coordinates. 79,80 The advantage of VPT2 is that it scales linearly with the number of normal modes, and is thus typically more than an order of magnitude cheaper than the VSCF-based methods. Furthermore, there has been numerous refinements to VPT2 over the years such as the generalized VPT2 (GVPT2) method which is able to account for resonance effects arising from the singularities caused when one frequency is either close the sum of two others or twice the value another. 79,81–84 A more recent development was the extension of the GVPT2 into a hybrid degeneracy-corrected VPT2 (HDCPT2) by Bloino, 85 which unlike its predecessors performs well even when the couplings between the high and low frequency modes are large. As these methods yield both the fundamental anharmonic frequencies and ZPEs, they are typically combined with the simple perturbation theory (SPT) approach to obtain the vibrational partition function. 86–90 While popular, SPT is not the only possible way to calculate the vibrational partition functions and other approaches have been employed even in the case of H2 SO4 ·H2 O. 16 To the best of our knowledge, the different approaches discussed have not been compared with one another. It is nowadays commonly known that even approaches like VMP2 and GVPT2 often perform poorly for the floppy intermolecular large-amplitude motions present in many atmospheric complexes, which due to their low-frequency degrees of freedom tend to possess


Page 4 of 47

Page 5 of 47

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 largest impacts on the vibrational partition function. This implies that if accurate thermodynamic properties are required, an alternative treatment is necessary for these modes. Our aim in this study is to obtain reliable thermodynamic properties for the sulfuric acid monohydrate by accurately accounting for both the local and global anharmonicities using the so called anharmonic domain approximation. 91,92 In this approach, anharmonic effects for the most important vibrational motions, like the large amplitude motions coupling different conformers to one another, are treated by limiting the dimensionality of the potential energy surface (PES) to small one- to three-dimensional subspaces that contain a few strongly coupled vibrational degrees of freedom. In these subspaces, the vibrational problem is solved separately from the rest of the vibrational degrees of freedom using the variational method. We also study how different approaches to the vibrational partition function affect the thermodynamic properties. We compare our anharmonic domain calculations with results from more approximate treatments of global anharmonicity, such as Boltzmann averaging, to obtain a better estimate of how accurate these methods are. In the next section, we introduce the vibrational model used for the anharmonic domain calculations, and give a brief overview of the statistical thermodynamical methods used to account for global and local anharmonicities in this study. This is followed by computational details of the calculations in Section 3. The results of our quantum chemical and variational vibrational calculations are given in Section 4. The article ends with conclusions and implications for further research in Section 5.

2 2.1

Theory and Background Vibrational Model for Large Amplitude Motions

When the potential energy function consists of a single deep minimum the harmonic approximation can yield satisfactory results for the lowest vibrational states. Problems arise when the vibrational quantum number gets larger because of the increasing density of vibrational 5 ACS Paragon Plus Environment


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

states. The problem becomes more demanding for systems with multiple low-lying potential energy minima because the tunneling splittings of states can cause the harmonic approximation to fail even for the lowest vibrational states. A further complication arises as the dimensionality of the PES increases by three with each added atom. Because a complete variational anharmonic treatment requires knowledge of vast areas of the PES, the rapid escalation in the complexity of the PES already makes it difficult to treat molecular systems such as H2 SO4 and H2 SO4 ·H2 O fully variationally. To deal with the dimensionality problem, we divide the vibrational problems of sulfuric acid and its monohydrate into smaller isolated sub-systems that can be treated with the variational approach. In this article as in our previous contributions, we have called this local mode -type approach 93 the anharmonic domain (AD) approximation. Our aim is to include in the anharmonic domains the vibrational modes which result in a change of configuration, are strongly coupled to one another, and for which many excited states need to be calculated to obtain accurate thermodynamic properties. This includes, for example, the intermolecular modes which in general show substantially more anharmonicity than the intramolecular stretching and bending modes. 94 It has been demonstrated that with a careful choice of reduced-dimensional model, reasonably accurate energy levels can be obtained. 95,96 Because the high-frequency vibrations such as the O-H -stretching modes mainly affect the thermodynamic properties through their zero-point energies, the vibrational partition functions for these degrees of freedom can be accurately described by truncating the sum over states after the fundamental frequency. 63 In these cases, it is sufficient to employ a simple perturbative anharmonic method such as HDCPT2, which can be used to obtain the zero point energy and the lowest energy levels accurately. In this study, the anharmonic domain approximation was used for the large amplitude motions of both sulfuric acid and sulfuric acid monohydrate, while the rest of the vibrational motions were treated using HDCPT2. In contrast to the normal coordinates employed in the HDCPT2 calculations, curvilinear internal coordinates were employed in the anharmonic


Page 6 of 47

Page 7 of 47

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: Definitions of the coordinates included within the anharmonic domains AD H2 SO4 {φ1 , φ2 }

coordinate

definition

φ1 φ2

H1, O1, S1, and O2 H2, O2, S1, and O1

d4 d5 d6

H2, O2, S1 and O1 O5, O2, S1, and O1 H3, O5, O2, and S1

{a6 , a7 , d7 }

a6 a7 d7

O5, O2 and S1 H3, O5, and O2 H4, O5, H3, and O2

{d1 }

d1

H1, O1, S1, and O2

{r7 }

r7

O5 and O2

H2 SO4 ·H2 O {d4 , d5 , d6 }

domain calculation, as these coordinates are better suited for the description of large amplitude motions and usually decrease the couplings between different modes. 97 The drawback of curvilinear internal coordinates compared to rectilinear normal ones is that the kinetic energy operator becomes complicated as discussed in Section 2.2. The curvilinear internal coordinates employed in this study for H2 SO4 and H2 O·H2 SO4 are defined in terms of the atom designations shown in Figure 1 and are summarized in Table 1. For sulfuric acid, the active coordinates that formed the anharmonic domain were the two OH-dihedral angles φ1 and φ2 , where φ1 is the dihedral angle between atoms H1, O1, S1, and O2, while φ2 is the dihedral angle between atoms H2, O2, S1, and O1. These OH-rotations connect the two conformers of sulfuric acid. For sulfuric acid monohydrate, eight coordinates were considered for the anharmonic treatment. Five of these were dihedral angles: The d1 coordinate was the angle between H1, O1, S1, and O2; the d4 coordinate was the angle between H2, O2, S1 and O1; the d5 coordinate was the angle between O5, O2, S1, and O1; the d6 coordinate was the angle between H3, O5, O2, and S1, while the d7 coordinate was the angle between H4,



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 8 of 47

O5, H3, and O2. Two coordinates in the domains were bond angles: the a6 coordinate was the angle between O5, O2, and S1, while the a7 coordinate was the angle between H3, O5, and O2. Finally, the r7 coordinate was the distance between O5 and O2. Based on a series of constrained geometry optimizations where one of the coordinates was fixed to a given value and the changes in the other coordinates were monitored, we were able to divide these coordinates into two three-dimensional and two one-dimensional domains. The first of the three-dimensional domains consisted of coordinates a6 , a7 , and d7 , where d7 was the coordinate connecting the two lowest energy conformers of sulfuric acid monohydrate. The second three-dimensional domain consisted of the coordinates d4 , d5 , and d6 . The one-dimensional domains consisted of the coordinates r7 and d1 , where d1 is the coordinate connecting the the lowest energy conformer with the third-lowest energy conformer.

O3

O4

O4

O3

S1 H2

S1

O2

O2

H3

O1

O1

H1

H2

H1

H4

O5

Figure 1: The atomic designations used for the definition of the anharmonic domains.

2.2

Kinetic Energy Operator

In the AD approximation, the 3N -dimensional coordinate space is partitioned into a group of active coordinates A = {q1 , q2 , . . . qA } that constitute the anharmonic domain and a group of P = 3N − A passive coordinates P = {qA+1 , qA+2 , . . . , qA+P }. The presence of rigid 8 ACS Paragon Plus Environment

Page 9 of 47

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


constraints affects the kinetic energy operator of the anharmonic domain. In particular, the rigid constraints change the contravariant metric tensors of the AD space which are defined as g (i,j) =

N X

gα(i) · gα(j) ,

(1)

α=1 (i)

where the sum goes over all nuclei in the molecule N and i, j ∈ A. The quantities gα are called measuring vectors. In the unconstrained case, the measuring vectors for the vibra(i)

tional coordinates hα are obtained by taking the mass-weighted gradient of the vibrational (i)

coordinate hα =

√ 1 ∇ α qi . mα

However, in the rigidly constrained case, the measuring vectors

are obtained by adding the froze-mode correction term derived in ref. 98 to the unconstrained (i)

measuring vector hα : gα(i)

=

hα(i)

−

P X

([I])

λ(i[I])

kα , k ([I][I])

(i)

(j)

(3)

(i)

(j)

(4)

I=1

(2)

where k

(ij)

=

N X

kβ · kβ ,

β=1

λ(ij) =

N X

hβ · kβ ,

β=1

and the bracketed index counter is defined by the equation [I] = A + P + 1 − I. The vectors ([1])

([2])

([P ])

kα , kα , . . . , kα ([1])

are obtained from the unconstrained measuring vectors of the passive

([2])

([P ])

coordinates hα , hα , . . . , hα

k([I+1]) α

with the recursive formula

=

h([I+1]) α

−

I X

([J])

([I+1][J])

λ

J=1

([1])

with kα

([1])

= hα

kα , k ([J][J])

(5)

for all nuclei α.

With the measuring vectors calculated from equation (2), the kinetic energy operator can be written for the zero angular momentum states as 99



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

A−3 A−3 2 X −1/2 X ∂ 1 1 ∂g h ¯ ∂ † (ij) Tˆ = + −1/2 g (ij) , pi g pj = − 2 ij 2 ij ∂qi g ∂qi ∂qj

Page 10 of 47

(6)

where g = det g (ij) , pj = −i¯h ∂q∂ j . Because the total angular momentum quantum number is set to zero, the three rotational coordinates present in the active coordinate space make no contributions and the sum is only over the active vibrational degrees of freedom. The Jacobian of the coordinate transform is obtained from J = g −1/2 . With equation (6), the volume element of integration is of the form dτ = Jdq1 dq2 . . . dqA . If one wishes to use a weight function of one for integrations, i.e., dτ = dq1 dq2 . . . dqA , the kinetic energy becomes 100 Tˆ1 = g −1/4 Tˆg 1/4 .

(7)

After simplification, this transforms into 97,101,102 A−3

1 X (ij) Tˆ1 = pi g (q)pj . 2 ij

(8)

Because its effect on the energies is usually small, 97,103 in equation (8), we have neglected the pseudopotential term Vˆ 0 (q) which arises from setting the weight function equal to one. By comparing the relative magnitudes of the effects of the g-matrix-elements, this kinetic energy operator can be further simplified by neglecting those elements g (ij) which do not significantly affect the energy levels. In a previous work, it was found that for the calculation of thermodynamic properties it is sufficient to include only the frozen mode corrected diagonal g-matrix elements. 92 Furthermore, these can be treated as constants by taking the average over the g (ij) obtained with different values for the active coordinates. In this case, the kinetic energy operator is A−3

1 X (ii) 2 Tˆ1 = g p, 2 i a i

(9)

where the subindex a indicates that the g-matrix element is a constant, and not a function of the active coordinates A. 10 ACS Paragon Plus Environment

Page 11 of 47

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

Potential Energy Operator

For sulfuric acid, the potential energy surface had an identical form to the one we used in our previous article, 92 i.e.,

Vˆ (φ1 , φ2 ) =

3 X

A(1) n cos(nφ1 ) cos(nφ2 ) +

n=0

+

3 X

A(3) n

+

4 X

A(2) n [cos(nφ1 ) + cos(nφ2 )]

n=1

sin(nφ1 ) sin(nφ2 ) +

n=1

5 X

4 X

A(4) n [cos(nφ1 ) cos(φ2 ) + cos(φ1 ) cos(nφ2 )]

n=2

A(5) n [cos(nφ1 ) cos(2φ2 )+cos(2φ1 ) cos(nφ2 )]+

n=3

4 X

A(6) n [sin(nφ1 ) sin(φ2 )+sin(φ1 ) sin(nφ2 )].

n=2

+

4 X

A(7) n [sin(nφ1 ) sin(2φ2 ) + sin(2φ1 ) sin(nφ2 )]

(10)

n=3

The sines and cosines were grouped so that each term in equation (10) is totally symmetric. A one-dimensional free rotor basis was used for the two active coordinates φ1 and φ2 , which makes the calculation of Hamiltonian matrix elements as easy as possible with equations (9) and (10). For the four anharmonic domains of sulfuric acid monohydrate, several different forms of potential energy surfaces and basis set types were employed. For the r7 stretch, the potential energy operator was expressed in terms of the Morse variable z7 which is defined as z7 = 1 − e−aR R7 , where R7 = r7 − r7,e , aR is the Morse parameter, and r7,e is the equilibrium bond length for the coordinate r7 . The potential energy surface for the r7 coordinate had the form V (z7 ) =

10 X

Bk z7k ,

(11)

k=2

and harmonic oscillator eigenfunctions were used as the basis for this calculation. For the



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 12 of 47

one-dimensional d1 domain, the potential energy surface had the form Vˆ (d1) =

6 X

Cn(1) cos(nd1 ) +

n=0

4 X

Cn(2) sin(nd1 ),

(12)

n=1

with the one-dimensional free rotor basis employed in the calculation. Finally, for both of the two three-dimensional domains, the potential energy surface had the form:

Vˆ (q1 , q2 , q3 ) = Vˆ1 (q1 ) + Vˆ2 (q2 ) + Vˆ3 (q3 ) + Vˆ12 (q1 , q2 ) + Vˆ23 (q2 , q3 ) + Vˆ13 (q1 , q3 ),

(13)

where qi = Qi − Qi,e is the bond angle displacement coordinate, Qi is an instantaneous bond angle, and Qi,e is its equilibrium value. The terms on the right of equation (13) are defined by Vî (qi ) =

12 X

(i)

Dk qik ,

(14)

k=2

and Vîj (qi , qj ) =

11 X 11 X

(i,j)

Dk,l qik qjl .

(15)

k=1 l=1

2.4

Approaches to the vibrational partition function–local anharmonicity

If the full vibrational problem can be solved with sufficient accuracy, for example, by a variational calculation, then the need for separate treatment of local and global anharmonities becomes superfluous. In this case under the ideal gas assumption, the vibrational partition function qvib as the sum over molecular vibrational energy levels already contains all appropriate information about the different conformers. When this kind of solution of the vibrational problem is unfeasible, one can combine the anharmonic domain approximation with a less accurate anharmonic treatment for those passive vibrational degrees of freedom excluded from the AD. Typically, these degrees of freedom are treated as Morse or Rosen-Morse (Pöschl-Teller) oscillators. While in the diatomic case the exact solution for 12 ACS Paragon Plus Environment

Page 13 of 47

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 partition function can be obtained by using the Jacobi theta function and the Poisson summation formula, 104 in the case of polyatomic molecules it is necessary to make simplifications. All approaches adopted in this study are based on the anharmonic expression for the energy of the molecule which in non-degenerate cases is given by

E({ni }) = X0 +

P X

1 hνi (ni + ) + 2 i=A+1

P X j=A+1, A+1≤i≤j

1 1 Xij (ni + )(nj + ), 2 2

(16)

where the sum goes over the vibrational degrees of freedom excluded from the anharmonic domains, except in cases where the AD approximation was not employed or when couplings between the AD and the passive degrees of freedom were included in the calculations. The quantity ni is the quantum number of mode i, {ni } denotes the group of quantum numbers for all the different modes, X0 is the anharmonicity term for the ZPE, νi is the harmonic frequency of mode i, and Xij is the anharmonic coefficient between modes i and j. Upon extraction of the ZPE from the sums, this expression becomes

E({ni }) = EZPE +

P X i=A+1

P X

hνi ni +

j=A+1, A+1≤i≤j

1 Xij (ni nj + (ni + nj )). 2

(17)

An intuitive approach for calculating the partition function is to sum only over those energy levels of equation (17) with energies below a certain threshold value. While this is feasible for the full dimensional vibrational problem for medium sized systems like sulfuric acid with threshold values of about 3000 cm−1 the presence of several large amplitude motions can make it intractable, as is the case for sulfuric acid monohydrate. An additional problem with Morse type-energy expressions presented in equations (16) and (17) is that if the Xij coefficients are negative, then the corresponding terms in the sum may cause the partition function to diverge. 104 To overcome the difficulties of the direct summing scheme, perhaps the most common



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 14 of 47

approach is to use SPT 84,86–90 where the vibrational partition function is written as

qvib = ΠPi=A+1

e−EZPE /kB T . 1 − e−hνi /kB T

(18)

The quantity EZPE is the zero point energy of the modes P and νi is the fundamental anharmonic frequency. In equation (18), the equal spacing of the energy levels which is present in the harmonic approximation is retained, but the distance between neighbouring energy levels corresponds to the anharmonic fundamental frequency, typically obtained from VSCF, VPT2 or HDCPT2 calculations, for example. A second important change compared to the harmonic model is the more accurate treatment of the zero point energy in SPT because the equilibrium constant is exponentially dependent on it. A final alternative to the calculation of qvib is to insert equation (17) to the definition of the partition function and use a Taylor expansion up to the first order for the exponential term corresponding to the anharmonic energy correction. This approach was utilized, for example, by K´ urten et al. 16 and results in a vibrational partition function of the form

qvib = qvib,0 e−EZPE /kB T

P X e−ξi e−ξi (1 + e−ξi ) χ − χii 1− ii (1 − e−ξi )2 1 − e−ξi i=A+1 i=A+1 P X

 −EZPE /kB T

−qvib,0 e

−ξi −ξj P  X e e  χij +  −ξ i 1−e 1 − e−ξj j=A+1, A+1

Effects of Global and Local Anharmonicities on the ... - ACS Publications

Recommend Documents