Efficient Computation of Difference Vibrational Spectra in Isothermal

Subscriber access provided by CORNELL UNIVERSITY LIBRARY

Article

Efficient Computation of Difference Vibrational Spectra in Isothermal-Isobaric Ensemble Tatsuya Joutsuka, and Akihiro Morita J. Phys. Chem. B, Just Accepted Manuscript • DOI: 10.1021/acs.jpcb.6b07121 • Publication Date (Web): 10 Oct 2016 Downloaded from http://pubs.acs.org on October 16, 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.

The Journal of Physical Chemistry B 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 31

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 Difference Vibrational Spectra in Isothermal-Isobaric Ensemble Tatsuya Joutsuka† and Akihiro Morita∗,†,‡ †Department of Chemistry, Graduate School of Science, Tohoku University, Sendai 980-8578, Japan ‡Elements Strategy Initiative for Catalysts and Batteries (ESICB), Kyoto University, Kyoto 615-8520, Japan E-mail: [email protected]

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

Abstract Difference spectroscopy between two close systems is widely used to augment its selectivity to the different part of the observed system, though the molecular dynamics calculation of tiny difference spectra would be computationally extraordinary demanding by subtraction of two spectra. Therefore, we have proposed an efficient computational algorithm of difference spectra without resorting to the subtraction. The present paper reports our extension of the theoretical method in the isothermal-isobaric (NPT ) ensemble. The present theory expands our applications of analysis including pressure dependence of the spectra. We verified that the present theory yields accurate difference spectra in the NPT condition as well, with remarkable computational efficiency over the straightforward subtraction by several orders of magnitude. This method is further applied to vibrational spectra of liquid water with varying pressure, and succeeded in reproducing tiny difference spectra by pressure change. The anomalous pressure dependence is elucidated in relation to other properties of liquid water.

1 Introduction Difference spectroscopy is a very useful and versatile means to augment the selectivity of the spectral measurement. 1–8 By imposing a change on the system under investigation and taking the difference of spectra, we can selectively probe the part of the observed system associated with the change. There are many applications of the difference spectroscopy using various kinds of change. Just for illustrations, these changes include isotope exchange to assign complicated vibrational spectra, 7 solute species embedded in solutions, 1,2 voltage of electrodes in electrochemical systems to perturb the electric double layer, 3,4 partial structure of protein molecules, 6 etc. On some occasions, thermodynamic conditions such as temperature or pressure are perturbed to observe the difference spectra. 5,8–10 In contrast to these wide applications of the difference spectral measurements, computational analysis of these difference spectra by molecular dynamics (MD) simulation has not been developed so far. Although the MD simulation allows for calculating vibrational spectra and providing a molecular-level insight into the structure and dynamics of the 2 ACS Paragon Plus Environment

Page 2 of 31

Page 3 of 31

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


system of interest, 11–33 the straightforward calculation of the tiny difference spectra would require prohibitively large computational cost of yielding sufficient statistical convergence compared to that of the original spectra. This difficulty becomes particularly serious when the perturbation is quite small so that the relative difference spectra to the original ones are the order of ∼1/100 or less. In order to overcome the difficulty, we have proposed efficient computational algorithms using MD simulation to calculate difference spectra, rather than to calculate two spectra and to take the difference. The high accuracy and efficiency of this theory were demonstrated in the cases of liquid argon and water. 34–36 The new method opens an avenue to analyze various difference spectra with modest cost of computation by MD simulation. In our previous papers 34–36 the theory of difference spectroscopy was given in the microcanonical (NV E) and canonical (NV T ) ensembles. The present paper extends the computational algorithm to the isothermal-isobaric (NPT ) ensemble, and examine the efficiency, accuracy and applicability. Use of the NPT ensemble extends our applications in two ways. First, it allows us to deal with pressure change. The difference spectra with pressure change have been investigated in various phenomena, such as the stability of proteins, 5 effect on the hydrogen-bond network of water, 37–41 etc. Second, use of the isobaric ensemble allows us to discuss various difference spectra under a constant pressure rather than under a constant volume. In many experimental situations a condition of constant pressure is more relevant than that of constant volume. The difference in the two ensembles is emphasized in a limited size of MD system. For example, if we consider a situation that solute molecule is added to a solvent, the solute molecule should increase the pressure in a cell of a constant volume, or expand the size of the cell under a constant pressure. The NV T and NPT ensembles would lead to different results of the calculated spectra. Since we deal with tiny difference in the spectra, the difference in the ensembles would have substantial influence on the calculated results. Furthermore, we apply the present theory to liquid water to demonstrate its validity. Liquid water is known to exhibit a number of anomalous properties, 38,39,42–47 and some anomalous properties are observed in their pressure dependence. Typical examples are seen in the infrared (IR)



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

and Raman spectra of liquid water, which exhibit a turn-over behavior in response to increasing pressure. 38,39 We elucidate the anomalous behavior of vibrational spectra with the present theory, as detailed in Section 4. We note that relative changes of these vibrational spectra are in the order of ∼ 0.01 in the range of pressure up to 10000 atm, and hence the present theory is quite useful to analyze the mechanism of the spectral perturbation. The remainder of this paper is constructed as follows. Section 2 presents the theory and methodology of difference spectra in the NPT ensemble. The accuracy and efficiency of the theory in the NPT condition are examined and validated through MD calculation of liquid argon in Section 3. Then the present method is applied to elucidate the pressure dependence of liquid water in Section 4. Finally, the conclusion is given in Section 5.

2 Formulation in N P T Ensemble The essential idea of the theory of difference spectra has been given in our previous papers, 34–36 and thus we focus on the extension to the NPT ensemble in this section. First the equations of motion are presented in Section 2.1, and the extended theory of difference spectra is given in Section 2.2.

2.1 Equation of Motion We suppose N-particle system in a d-dimensional space, which is described with Cartesian coordinates r j and their conjugate momenta p j (= m j r˙ j with the mass m j ) of the jth particle ( j = 1 ∼ N). In the NPT ensemble, this system is attached to the Nosé-Hoover thermostat 48,49 and a barostat, 49–53 which apply the external temperature T and pressure P to the system, respectively. The thermostat controls the temperature by introducing the additional coordinate ξ . The barostat allows the volume of a simulation cell V to fluctuate, and the volume fluctuation is described with the piston coordinate ε = lnV . The equations of motion for the NPT ensemble are given by the Martyna-Tobias-Klein algo-


Page 4 of 31

Page 5 of 31

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


rithm, 51–53 ( ) pξ d pε p˙ j = F j − 1 + p j − p j, g W Q

p j pε r˙ j = + r j, mj W

p2j

N

Qξ˙ = pξ ,

p˙ξ =

∑ mj

+

j=1

W ε˙ = d pε ,

p2ε − (g + 1)kB T, W

p˙ε = dV [Pint (t) − P] +

d g

N

p2j

pξ

∑ m j − Q pε .

(1) (2) (3)

j=1

Eqs 1, 2 and 3 are the equations of motion for the system (r, p), thermostat (ξ , pξ ) and barostat (ε , pε ), respectively. F j is the force applied to the jth particle and defined as the derivative of the potential energy U by F j = −∂ U(r,V )/∂ r j . g denotes the number of the degrees of freedom for the original system, and g = dN with no constraint. Q and W are the effective masses for the thermostat and barostat, respectively, and kB is the Boltzmann constant. In eq 3, the internal pressure Pint (t) is defined as 1 Pint (t) = dV

[

N

∑

j=1

p2j mj

] +rj ·Fj −

∂ U(r,V ) . ∂V

(4)

In the following we consider isotropic fluctuations of volume for simplicity, though the extension to an anisotropic case is straightforward. The above Martyna-Tobias-Klein algorithm is proven to generate the correct isothermal-isobaric distribution at equilibrium, 54

ρ (r, p, pξ , pε ,V ) = exp[−β H]/Z

(5)

where β = 1/(kB T ) and p2ξ

p2j

p2ε +U(r,V ) + H = ∑ + + PV 2Q 2W j 2m j ∫ ∞

Z =

∫

dV 0

∫

N

d r

∫

N

d p

∫

d pξ

d pε exp[−β H] ≡


(6) ∫

dΓ exp[−β H(Γ)]

(7)


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 6 of 31

Here we define the set of extended variables Γ = (r, p, pξ , pε ,V ). The NPT ensemble is generated also when no external force is involved (∑ j F j = 0). The Martyna-Tobias-Klein algorithm is proved to give the correct NPT ensemble, unlike others such as Berendsen 55,56 or the original Hoover 49 thermostat and barostat. We found that correct ensemble is crucial to the theory of difference spectra, arguably because the difference spectra are sensitive to tiny deviation of the ensemble by perturbation. This is consistent with the tendency that fluctuation is generally more sensitive to the statistical ensemble than average. 57

2.2 Difference Spectra Difference spectra are essentially defined between a pair of slightly different states. Here we suppose that both states obey the NPT ensemble. The difference between the two states is attributed to some different parameter(s) involved in those equations. In our previous papers 34–36 we assumed that the difference of the systems is represented with their different Hamiltonians and constructed the theory of difference spectra. Here we present a slightly different formulation in the NPT ensemble since the present equations of motion are not necessarily based on the Hamiltonian. To begin with, the line shape function of an arbitrary spectrum I(ω ) is defined by the Fourier transform of the corresponding time correlation function C(t) = ⟨B(0)B(t)⟩, 1 I(ω ) = 2π

∫

dte−iω t C(t),

(8)

where B denotes the relevant physical quantity to the spectrum in question, such as velocity, dipole moment or polarizability. ω is the angular frequency, and the angle bracket ⟨· · · ⟩ means the ensemble average. C(t) in eq 8 is represented as C(t) = ⟨B(0)B(t)⟩ =

∫

dΓρ (Γ)B(Γ) exp(iLt)B(Γ)

(9)

where the equilibrium distribution function ρ (Γ) is given in eq 5. exp(iLt) is the time evolution


Page 7 of 31

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


operator to propagate B for time t. The Liouville operator L is defined as ) ( ∂ ∂ ∂ ∂ ∂ ∂ + p˙ j + ξ˙ + p˙ε iL = ∑ r˙ j + p˙ξ + V˙ ∂rj ∂pj ∂ξ ∂ pξ ∂V ∂ pε j

(10)

(Note that ε˙ (∂ /∂ ε ) = V˙ (∂ /∂ V ).) We note that C(t) in eq 8 could be replaced with a cross-correlation function (C(t) = ⟨B(0)B′ (t)⟩ with B′ ̸= B) in some cases, e.g. the sum frequency generation spectroscopy. 17,58 The extension to the cross-correlation function is straightforward, and we employ the autocorrelation function in the following discussion. Then, we suppose two states, namely the reference state and the perturbed state. The difference in the states could be generally defined with any different parameter except for the Γ variables themselves, such as the pressure P, temperature T , and potential U(r,V ). The difference spectrum is accordingly given with the difference in the time correlation functions ∆C(t) as 1 ∆I(ω ) = 2π

∫

dte−iω t ∆C(t),

(11)

where ∆C(t) = ⟨B(0)B(t)⟩ − ⟨B(0)B(t)⟩0 ∫

=

dΓρ (Γ)B(Γ) exp(iLt)B(Γ) −

∫

dΓρ0 (Γ)B(Γ) exp(iL0t)B(Γ).

(12)

ρ and L are the distribution function and Liouville operator, respectively, in the perturbed state, and ρ0 and L0 are those in the reference state. To distinguish the two states, the quantities defined at the reference state are represented with subscript 0 in the following. Eq 12 indicates that the difference in C(t) is expressed through that in the distribution functions (ρ and ρ0 ) or that in the Liouville operators (L and L0 ). Accordingly we can decompose ∆C(t) in



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 31

eq 12 into two terms, ∆C(t) = ∆Cc (t) + ∆Ct (t),

(13)

where ∫

∆Cc (t) =

∫

∆Ct (t) =

dΓ(ρ (Γ) − ρ0 (Γ))B(Γ) exp(iL0t)B(Γ)

(14)

dΓρ (Γ)B [exp(iLt) − exp(iL0t)] B(Γ),

(15)

The first term ∆Cc (t) including (ρ − ρ0 ) is called the configuration term, and the second term ∆Ct (t) including [exp(iLt) − exp(iL0t)] is called the time evolution term. The decomposition in eq 13 is the essence of the present theory of difference spectroscopy. The details of the calculation procedure of the two terms were described in our previous paper, 36 and the procedure is straightforwardly applied to the present case of NPT ensemble. The configuration term is rewritten as ∆Cc (t) =

⟨exp[−β ∆H(0)]B(0)B(t)⟩0 − ⟨B(0)B(t)⟩0 ⟨exp[−β ∆H(0)]⟩0

(16)

where ∆H = H − H0 is the difference of H in eq 6 between the two states. The time evolution term is calculated by the non-equilibrium MD simulation. We employ the exact, non-perturbative calculation procedure of the two terms without resorting to the first-order perturbation expansion, 34 since the non-perturbative procedure is numerically more accurate and stable than the first-order expansion with almost comparable cost of computation. 36 In practical MD calculations of the difference spectra, the calculations are performed in both the forward and backward directions. The forward calculation treats the change from the reference state to the perturbed state, and the backward vice versa. The mutual treatment effectively doubles the statistical sampling without a significant increase in the computational cost, and also allows us to check the coincidence of the results in the two directions. The coincidence provides an evidence


Page 9 of 31

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


for reliability of the calculation. We also note that the present calculation of difference spectra provides the vibrational spectra at the reference and perturbed states with no additional cost of computation.

3 Validation with Liquid Argon The accuracy and efficiency of the present theory in the NPT ensemble are examined by calculating difference vibrational spectra of liquid argon for some typical perturbations, i.e. pressure, potential, temperature and mass. The validation of the present theory (analytical method) is performed in comparison to the straightforward subtraction of two spectra (numerical method). The latter is regarded as a rigorous method to yield the difference spectra, and thus used as a reference to validate the accuracy and efficiency of the former. In this section, a small model system of liquid argon is employed to allow for the numerical method, which is very time consuming as we mentioned in Section 1. After we confirm the reliability of the present analytical method in comparison to the numerical one, we apply the present analytical method to larger systems in Section 4.

3.1 MD conditions The liquid argon is modeled with N = 30 molecules enclosed in a cubic cell with the threedimensional periodic boundary conditions. The NPT ensemble is employed with the temperature T = 86 K and the pressure P = 2500 atm unless otherwise noted. The intermolecular interactions are represented with the pairwise Lennard-Jones (LJ) function, U(r) = 4ε

[( ) σ 12 r

−

( σ )6 ] r

,


(17)


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 10 of 31

˚ 34,59 and the mass of each molecule is 40 amu. The vibrational where ε = 119.8 K, σ = 3.405 A, spectra are monitored with the spectrum of velocity autocorrelation function, I(ω ) =

1 2π

∫ ∞ −∞

dte−iω t

1 N

N

∑ ⟨v j (t)v j (0)⟩

(18)

j=1

The difference spectra in the NPT condition are considered for four types of perturbations, (1) pressure, (2) potential, (3) temperature and (4) mass, as summarized in Table 1. In the case (2), the parameter of one molecule is changed while those of the other molecules unchanged. In order to compare the analytical method with the numerical one, the amount of each parameter change is chosen at the convenience of both methods. A small perturbation could be readily treated by the analytical method, though it becomes extraordinary demanding to the numerical method to be converged. Thus the perturbation should be large enough at the convenience of the numerical method. On the other hand, it should be small enough to make the analytical method reliable. If the perturbation were too large, the coincidence between the forward and backward calculations would be worse in the analytical method. Consequently, the amplitude of the difference spectrum for the purpose of present validation is typically set in the order of 1/100 of the original spectrum. Table 1: Four types of perturbations (1)-(4), and change in each corresponding parameter. perturbation (1) pressure (2) potential(∗,∗∗) (3) temperature(∗∗) (4) mass(∗∗) P /atm ε /K T /K m /amu reference 500 119.8 86.0 40 perturbed 490 131.78 86.172 40.8 (∗)

The parameter of one molecule is changed in (2).

(∗∗)

P = 2500 atm in (2)-(4).

˚ in case In the MD calculations, the cutoff distance for the LJ interaction is set to rcut = 5.0 A ˚ in cases (2)-(4) (P = 2500 atm) so that twice the cutoff (1) (P = 490 ∼ 500 atm) or rcut = 4.5 A distance 2rcut does not exceed the varying box length in each case. The calculation of the internal pressure takes account of the long-range correction beyond the cutoff distance. 60 We also adopted


Page 11 of 31

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 tapering function f (x) to the LJ potential near the cutoff distance,    1 (x ≤ 0)    f (x) = 1 − x3 {10 − x(15 − 6x)} (0 < x < 1)      0 (1 ≤ x)

(19)

where x = (r − 0.95rcut )/(0.05rcut ). f (x) in eq 19 is multiplied to the LJ potential when the distance r is in the range of 0.95 rcut < r < rcut (0 < x < 1) to avoid discontinuous force. In the case (2), one molecule has a slightly different LJ parameter ε than the others. Thus the LJ interaction in eq 17 is calculated using the Lorentz-Berthelot combining rule (arithmetic mean for

ε ). The time propagation was carried out by the velocity Verlet algorithm 60 with a time step ∆t of 5 fs. The effective masses for the thermostat and barostat are taken to be Q = gkT /ω p2 and W = (g + d)kT /ωb2 , respectively. 51 The thermostat and barostat parameters are ω p = 0.64 and

ωb = 1.6 × 10−3 ps−1 in case (1), and ω p = 1.6 and ωb = 8.0 × 10−4 ps−1 in cases (2)-(4). For each case of (1)-(4), the difference spectrum is calculated by both the numerical and analytical methods. Regarding the numerical method, we prepare 320-1600 different initial configurations, and generate independent MD trajectories for 2.5 ns of equilibration per each configuration with both the reference and sample states (time developments). Then the production runs at the reference and sample states are carried out to calculate the difference in the vibrational spectra. Regarding the analytical method, 160-320 independent equilibration runs are generated for 0.5 ns, and the subsequent production run is conducted. The statistical sampling by either method is taken sufficiently in relation to the convergence time discussed in Section 3.3.

3.2 Accuracy Pressure change:

We begin by discussing the accuracy of calculated difference spectrum in case

(1) pressure change. Figure 1(a) is the reference vibrational spectrum I0 (ω ) for liquid argon at 500 atm, and Figure 1(b) is the difference spectrum with the pressure change from 500 to 490 atm, ∆I(ω ) = I(ω ) − I0 (ω ). The original lineshape I0 (ω ) has a maximum at ω = 0, and shows no vi11 ACS Paragon Plus Environment


brational resonance at this pressure. The difference spectrum ∆I(ω ) has a relative intensity on the order of 1/100 of the original at ω = 0. The lineshape of ∆I(ω ) indicates that the high-frequency component decreases and the low-frequency increases at lower pressure. This change with decreasing pressure is intuitively obvious because the density decreases and the intermolecular interactions weakens at lower pressure. Figure 1(b) shows the analytical and numerical results in black and red lines, respectively, which almost coincide with each other. The analytical spectrum is the average of the forward and backward calculations, and they are separately shown in the inset of Figure 1(b). The forward and backward results overlap with each other, indicating reliability of the present calculations. The analytical result is decomposed into the configuration term (blue dashed line) and the time evolution term (green dashed line) in the main panel of Figure 1(b). This decomposition shows that the configuration term dominates the difference spectrum. The time evolution term has little contribution in the pressure change, arguably because the pressure change has indirect influence on the force acting on the molecules and thus on their time evolution. (a)

(b)

1

Norm. intensity

Norm. intensity

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 31

0.015

0.01

0.005

0

0

0

50

100

150

200

0

Frequency / cm -1

50

100

150

200

Frequency / cm -1

Figure 1: (a) Vibrational spectrum I0 (ω ) for liquid argon at P = 500 atm. (b) Difference spectrum with pressure change from P = 500 to 490 atm, ∆I(ω ) = I(ω ) − I0 (ω ). The analytical and numerical results are plotted by black and red lines, respectively. The analytical result is decomposed into the configuration term (blue dashed line) and the time evolution term (green dashed line). The inset in (b) shows the forward (black) and backward (red) results. The ordinates in (a) and (b) are commonly normalized with the maximum intensity of I0 (ω ).

Other parameters:

We next discuss the accuracy in the other cases (2)-(4) at constant pressure.

Figure 2(a) is the reference vibrational spectrum I0 (ω ) for liquid argon at P = 2500 atm. Here 12 ACS Paragon Plus Environment

Page 13 of 31

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 intensity maximum moved to ω = 33 cm−1 as a consequence of dense packing. We set the reference state at a higher pressure to let the difference spectrum sensitive to the parameter change. Figures 2(b)-(d) display the difference spectra for the perturbation on (2) potential, (3) temperature and (4) mass in Table 1, respectively. The notations of the lines and colors are the same as in Figure 1. In each case, the analytical result agrees with the numerical one within the error bar of the latter. The statistical uncertainties of the numerical method are evaluated as the standard deviations of the sampling in Section 3.3. In fact, the error bar of the numerical method arises from very slow convergence of statistical sampling, as we argue in Sec. 3.3. The inset of each panel shows that the forward and backward results of the analytical method reasonably agree with each other though the deviation is discernible. Regarding the band shapes of the difference spectra, panel (b) indicates a blue shift while panel (c) and (d) red shifts. In panel (b) the increasing

ε of one molecule strengthens the intermolecular interaction with other molecules, while in panel (c) increasing T randomizes the structure and in panel (d) increasing mass slows the vibrational motion. Regarding the relative contributions of configuration and time evolution terms, the configuration term has major contribution in (b) and (c) though the time evolution term is not negligibly small. In case (d), the configuration and time evolution terms have the opposite signs. Comparing the present results with previous ones in the NV T ensemble, 34,36 we notice in case (c) that the previous NV T calculation showed a negligible time evolution term, whereas the present NPT calculation gives a minor but non-negligible contribution of the time evolution term.

3.3 Efficiency The efficiency of calculation is evaluated by the amount of statistical sampling to get converged results of the difference spectra. We monitor the Student confidence interval 36,61 δ (ω ,t) of 95 % confidence for the calculated difference spectra ∆I(ω ) as a function of the sampling time t, and confirmed that δ (ω ,t) decreases with increasing sampling time t. The convergence time is thereby defined with the sampling time t to realize δ (ω ,t)/|∆I(ω )| = 0.1 (10 % relative uncertainty) at 13 ACS Paragon Plus Environment


(a)

1

(b)

Norm. intensity

Norm. intensity

0.005

0

0

−0.005

−0.01 0

50

100

150

200

0

Frequency / cm -1

50

100

150

200

−1

(c)

(d)

0.02

Norm. intensity

Frequency / cm

0

0.01 Norm. intensity

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 31

0.005

-0.02

0 -0.04

0

50

100

150

0

50

200

−1

100

150

200

Frequency / cm -1

Frequency / cm

Figure 2: (a) Calculated vibrational spectrum for liquid argon at 2500 atm, and difference spectra for (b) potential (c) temperature (d) mass changes at the pressure. The notations of the lines and colors are the same as in Figure 1. These difference spectra are normalized with the maximum intensity in panel (a). Each panel of (b)-(d) shows the error bar of the numerical method at the frequency of maximum absolute amplitude.


Page 15 of 31

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 frequency ω of maximum absolute amplitude |∆I(ω )|. This procedure of error estimation is detailed in our previous publications. 36 The results of convergence times in four cases (1)-(4) by the numerical and analytical methods are summarized in Table 2. Table 2: Convergence times of the difference spectra in four cases (1)-(4) using the numerical and analytical methods. The ratios of the two methods are also presented. The values in the parentheses are the convergence times in the N V T ensemble for comparison. perturbation (1) pressure (2) potential (3) temperature (4) mass numerical / µ s 129 1470 (375) 1050 47.4 analytical / ns 34.6 47.8 (33.2) 31.3 (1.12) 0.831 (0.187) ratio 3730 30800 (11300) 33500 57000

It is clear in Table 2 that the analytical method achieves the convergence much faster than the numerical method by factors ranging from 3730 to 57000. These results indicate superior efficiency of the analytical method to the numerical one by several orders of magnitude in the NPT ensemble. Table 2 also shows the convergence times of difference spectra calculated in the NV T ensemble, and thus allows us to compare the computational efficiency in the NPT ensemble to NV T . We find that the NPT ensemble generally requires several times as much sampling time as NV T either by numerical or analytical method, arguably because that the volume fluctuation in the NPT ensemble requires additional computational effort to sample the phase space. It is also noted that the original vibrational spectrum itself I0 (ω ) converges in 85.1 and 110 ps in the NV T and NPT ensembles, respectively, by the same criterion (10 % relative uncertainty). In either ensemble, the analytical method of calculating difference spectra exhibits remarkably improved efficiency over the numerical method.

4 Pressure Dependence of Liquid Water In this section we apply the present theory to difference spectra of liquid water by pressure change, and elucidate the pressure dependence of experimental vibrational spectra and related properties. Some of the anomalous properties of liquid water are observed in their pressure dependence. 15 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

For example, the viscosity 42,43 shows a minimum and the self-diffusion coefficient measured by NMR 44–46 shows a maximum around 1000–2000 atm at room temperature. X-ray diffraction experiment 47 reported an anomaly in intermolecular structure around 2000 atm. This anomaly has been understood by the following two competing effects of pressurization: [1] As pressure rises, the tetrahedral hydrogen-bond structure characteristic of water tends to be destroyed. Accordingly, the intermolecular hydrogen bond interactions become weaker. [2] Pressurization tends to decrease intermolecular distances, which is common to normal liquids, and accordingly the intermolecular interactions become stronger. 62 With increasing pressure the former effect initially prevails around ambient pressure, yet at about 2000 atm the latter effect takes over, leading to the anomaly for pressurization in liquid water. This picture was supported by infrared (IR) 38 and Raman 39 measurements of bulk water. However, some other experiments, such as neutron scattering 63 and Raman spectroscopy, 64,65 did not observe such behavior, presumably because the small spectral change is buried in the experimental uncertainty. In addition, a Monte Carlo simulation 66 has failed to confirm this picture. Therefore, we apply our analytical method of difference vibrational spectra for liquid water to resolve this issue in comparison with experimental IR 38,40,41,67,68 and Raman 37,39,69,70 spectroscopies.

4.1 MD conditions The MD simulation of liquid water was conducted with N = 512 water molecules in a cubic cell with the three-dimensional periodic boundary conditions. The model water is flexible and polarizable, developed previously by us with the charge response kernel. 71 The cutoff distance of LJ ˚ and the electrostatic interactions were calculated by the Ewald method 60 with interactions was 9 A ˚ −1 . The simulation time of the equilibration and production the Ewald coefficient set to 0.349 A runs for each perturbation was 16 ns. 160 independent simulations are initiated with an equilibration run for 100 ps. The temperature T was set to 300 K in the NPT ensemble, and the pressure was set to various values from P = 1 atm to 10000 atm. Otherwise, the computational procedure is the same as in Section 3.1. 16 ACS Paragon Plus Environment

Page 16 of 31

Page 17 of 31

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 calculations of IR and Raman spectra were performed with the time correlation function of dipole and polarizability. 72 The detailed formulas for these spectra are described in ref 36, where the self terms of correlation were used in the following discussion. Besides the vibrational spectra, we also calculated the site-site radial distribution functions gOO (r), gOH (r) and gHH , and the self diffusion coefficient D of water with varying pressure. The diffusion coefficient D is calculated by the following formula, 60

⟨ 1 D = lim t→∞ 6t

⟩ 1 N |r j (t) − r j (0)|2 . N∑ j

(20)

where r j is the Cartesian coordinate of the center of mass of the jth water molecule.

4.2 Vibrational spectra Figure 3 displays the difference vibrational spectra of liquid water in the OH-stretch region (3000– 4000 cm−1 ) for various pressure changes, P = 1 atm → 1000 atm, 1000 → 2000, 2000 → 3000, 3000 → 4000, 4000 → 5000 and 9000 → 10000. The four panels show (a) velocity correlation, (b) infrared, (c) VV Raman and (d) VH Raman difference spectra. Figure 3(a) manifests an intriguing trend in the difference spectra of velocity correlation with increasing pressure. These lineshapes indicate that the vibrational spectra show blue shift in relatively low pressures (1 → 1000 and 1000 → 2000), whereas red shift in high pressures (4000 → 5000 and 9000 → 10000). This turn-over behavior supports the picture of pressure dependence of water properties mentioned above. With increasing pressure, the hydrogen-bond structure is disrupted, resulting in reduced red shift (blue shift) of the OH stretch frequency, and then the dense packing at high pressure strengthen the intermolecular interaction to augment the red shift. The turning point is seen at around 2000 atm, which is consistent with the experimental spectroscopic measurements, 38,39 and this turning point has been observed by MD simulation for the first time. The present method can readily calculate the small difference spectra in Figure 3 on the order of 0.01 compared to the original spectrum. The IR difference spectra in Figure 3(b) show the analogous transition, consistent with the



(a)

(b) 1→1000 atm 1000→2000 atm 2000→3000 atm 3000→4000 atm 4000→5000 atm 9000→10000 atm

0.01

0.01 Normalized intensity

Normalized intensity

0.02

0 -0.01

0

-0.01

-0.02

-0.02 3000

3200

3400

3600

3800

4000

3000

3200

Frequency / cm -1

3400

3600

3800

4000

3800

4000

Frequency / cm -1

(c)

(d) 0

0



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 31

-0.02

-0.04 3000

3200

3400

3600

3800

4000

-0.02

-0.04

-0.06 3000

3200

Frequency / cm -1

3400

3600

Frequency / cm -1

Figure 3: Difference vibrational spectra in the OH-stretch region (ω = 3000 ∼ 4000cm−1 ) for the pressure change from P = 1 to 10000 atm in liquid water. (a) velocity correlation, (b) infrared, (c) VV polarized and (d) VH polarized Raman difference spectra. The black, red, green, blue, magenta and cyan lines denote the spectra for the changes 1 → 1000, 1000 → 2000, 2000 → 3000, 3000 → 4000, 4000 → 5000 and 9000 → 10000 atm, respectively. The thick arrows in (a) illustrate the direction of pressurization. The intensity of each difference spectrum is normalized with the maximum of the original spectrum at 1 atm.


Page 19 of 31

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


experimental IR measurements. 38,67 Comparing the line shapes of velocity correlation in (a) to those of IR spectra in (b), we find that the amplitudes in low frequency region below 3400 cm−1 are generally more emphasized in the latter. This feature is understood from the fact that the transition moment of OH stretch is augmented by strong hydrogen bonds in low-frequency region, 9,14,73 the trend of which has been commonly observed for various hydrogen-bonded systems. 74 In the Raman difference spectra in Figure 3(c) and Figure 3(d), the change in the low-frequency region is further emphasized so that the initial blue shift in low pressures (1 → 1000, 1000 → 2000) becomes hardly discernible. These calculation results may explain the previous experiment of Raman spectroscopy 64 that such transition was not observed clearly within the experimental uncertainty.

4.3 Other properties Related anomaly in other properties of structure and dynamics of liquid water is examined by the present MD simulation with varying pressure. Figure 4 shows the MD results of site-site radial distribution functions under the pressures ranging from P = 1 to 10000 atm. The most striking ˚ change with increasing pressure is the disappearance of the second peak in gOO (r) at r ≃ 4.5 A, which is indicative of the characteristic tetrahedral structure of water. The second peak around ˚ which exists at 1 atm becomes flat up to 3000 atm, and eventually a valley at 10000 atm. 4.5 A This change implies that the tetrahedral hydrogen bond structure is disrupted by pressurization. This overall trend with increasing pressure is in accord with previous Monte Carlo simulation and experiment. 66 Furthermore, the diffusion coefficient is known to be sensitve to the strength of intermolecular interaction and hydrogen bond. 75,76 The calculated self diffusion coefficient D of water molecules with varying pressure is displayed in Figure 5. The diffusion coefficient D at 1 atm and 300 K is calculated to 1.2 × 10−9 m2 /s, which underestimates the experimental value of 2.4 × 10−9 m2 /s. 77 Such underestimation is also seen in other polarizable water models. 78,79 Nevertheless, Figure 5 shows that the present MD simulation reproduces the turn-over behavior of the diffusion coeffi19 ACS Paragon Plus Environment


(a)

(b)

4

1 1000 2000 3000 5000 10000

3.5 3

atm atm atm atm atm atm

1.8 1.6 1.4 1.2

gOH

2.5 gOO

2 1.5

1 0.8 0.6

1

0.4

0.5

0.2

0

0

1

2

3

4

5

6

7

0

8

1

2

3

4

r/Å (c)

5

6

7

8

r/Å 1.8 1.6 1.4

gHH

1.2 1 0.8 0.6 0.4 0.2 0

1

2

3

4

5

6

7

8

r/Å

Figure 4: Site-site radial distribution functions (rdf) of liquid water under P = 1 ∼ 10000 atm. (a) OO (gOO ), (b) OH (gOH ), and (c) HH (gHH ) rdf. The black, red, green, blue, magenta and cyan lines denote the rdfs under P = 1, 1000, 2000, 3000, 4000, 5000, 10000 atm, respectively.

1.6 Diffusion coefficient / 10 -9 m 2/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 20 of 31

1.4

1.2

1

0

2000

4000

6000

8000

10000

Pressure / atm

Figure 5: Diffusion coefficient D under various pressures from 1 to 10000 atm.


Page 21 of 31

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


cient. 44,45 The pressure dependence shows a turning point around 2000 atm, while the experiment shows the turning point somewhat in low pressure side around 1000 atm. 44,45 Such turn-over behavior of the diffusion coefficient is quite consistent with the pressure dependence of the vibrational spectra discussed in Section 4.2. With increasing pressure, the hydrogen bonds are disrupted and accordingly the self diffusion of water molecules becomes faster. Conversely, for a higher pressure the dense packing leads to a reduced diffusion coefficient. All these results support the qualitative trend of the turn over calculated in the difference vibrational spectra with varying pressure.

5 Conclusion We extended the theory and efficient computational algorithm of difference vibrational spectra in the isothermal-isobaric (NPT ) ensemble. This extension also allows for treating general difference spectra by pressure change. The present theory retains its accuracy and remarkable efficiency in the NPT ensemble as well. We have verified that it reproduces rigorous difference spectra obtained by the numerical subtraction, while it reduce the computational cost by several orders of magnitude. We further applied this theory to the difference spectra of liquid water, and succeeded in reproducing and elucidating the turn over behavior by the pressure change. As pressure increases, the OH stretching band first shows blue shift and then turn to red shift. The present theory is demonstrated to be quite powerful to obtain such small difference spectra in a reliable manner. Such anomaly in the liquid water is further elucidated in relation to the radial distribution functions and diffusion coefficient. The present extension expands the applicability of molecular dynamics calculation of difference spectra in various realistic conditions. Now the present theory has grown to be a robust and practically useful tool to analyze a variety of experimental difference spectra.



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

Acknowledgement We thank Prof. Dor Ben-Amotz for useful discussion. This work was supported by the Grantsin-Aid for Scientific Research (Nos. JP25104003, JP26288003) by the Japan Society for the Promotion of Science (JSPS), Japan. The computations were carried out using the supercomputers in Research Center for Computational Science, Okazaki, Japan.

References (1) Davis, J. G.; Gierszal, K. P.; Wang, P.; Ben-Amotz, D. Water Structural Transformation at Molecular Hydrophobic Interfaces. Nature 2012, 491, 582–585. (2) Davis, J. G.; Rankin, B. M.; Gierszal, K. P.; Ben-Amotz, D. On the Cooperative Formation of Non-Hydrogen-Bonded Water at Molecular Hydrophobic Interfaces. Nat. Chem. 2013, 5, 796–802. (3) Yamakata, A.; Soeta, E.; Ishiyama, T.; Osawa, M.; Morita, A. Real-Time Observation of the Destruction of Hydration Shells under Electrochemical Force. J. Am. Chem. Soc. 2013, 135, 15033–15039. (4) Yamakata, A.; Shimizu, H.; Oiki, S. Surface-Enhanced IR Absorption Spectroscopy of the KcsA Potassium Channel upon Application of an Electric Field. Phys. Chem. Chem. Phys. 2015, 17, 21104–21111. (5) Taniguchi, Y.; Takeda, N. In Biological Systems Under Extreme Conditions; Taniguchi, Y., Stanley, H. E., Ludwig, H., Eds.; Biological and Medical Physics Series; Springer: Berlin, 2002; pp 101–120. (6) Furutani, Y.; Shimizu, H.; Asai, Y.; Fukuda, T.; Oiki, S.; Kandori, H. ATR-FTIR Spectroscopy Revealing the Different Vibrational Modes of the Selectivity Filter Interacting with


Page 22 of 31

Page 23 of 31

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


K+ and Na+ in the Open and Collapsed Conformations of the KcsA Potassium Channel. J. Phys. Chem. Lett. 2012, 3, 3806–3810. (7) Wang, R.; Sivakumar, V.; Johnson, T. W.; Hastings, G. FTIR Difference Spectroscopy in Combination with Isotope Labeling for Identification of the Carbonyl Modes of P700 and P700+ in Photosystem I. Biophys. J. 2004, 86, 1061–1073. (8) Walrafen, G. E.; Fisher, M. R.; Hokmabadi, M. S.; Yang, W.-H. Temperature Dependence of the Low- and High-Frequency Raman Scattering from Liquid Water. J. Chem. Phys. 1986, 85, 6970–6982. (9) Corcelli, S. A.; Skinner, J. L. Infrared and Raman Line Shapes of Dilute HOD in Liquid H2 O and D2 O from 10 to 90 ◦ C. J. Phys. Chem. A 2005, 109, 6154–6165. (10) Nagata, Y.; Hasegawa, T.; Backus, E. H. G.; Usui, K.; Yoshimune, S.; Ohto, T.; Bonn, M. The Surface Roughness, but not the Water Molecular Orientation Varies with Temperature at the Water-Air Interface. Phys. Chem. Chem. Phys. 2015, 17, 23559–23564. (11) Torii, H. Time-Domain Calculations of the Infrared and Polarized Raman Spectra of Tetraalanine in Aqueous Solution. J. Phys. Chem. B 2007, 111, 5434–5444. (12) Buch, V.; Tarbuck, T.; Richmond, G. L.; Groenzin, H.; Li, I.; Shultz, M. J. Sum Frequency Generation Surface Spectra of Ice, Water, and Acid Solution Investigated by an Exciton Model. J. Chem. Phys. 2007, 127, 204710. (13) Laage, D.; Stirnemann, G.; Sterpone, F.; Rey, R.; Hynes, J. T. Reorientation and Allied Dynamics in Water and Aqueous Solutions. Annu. Rev. Phys. Chem. 2011, 62, 395–416. (14) Skinner, J. L.; Auer, B. M.; Lin, Y.-S. Vibrational Line Shapes, Spectral Diffusion, and Hydrogen Bonding in Liquid Water. Adv. Chem. Phys. 2009, 142, 59–103. (15) Laage, D.; Stirnemann, G.; Hynes, J. T. Water Jump Reorientation and Ultrafast Vibrational Spectroscopy. J. Photochem. Photobiol. A Chem. 2012, 234, 75–82. 23 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

(16) Skinner, J. L.; Pieniazek, P. A.; Gruenbaum, S. M. Vibrational Spectroscopy of Water at Interfaces. Acc. Chem. Res. 2012, 45, 93–100. (17) Morita, A.; Ishiyama, T. Recent Progress in Theoretical Analysis of Vibrational Sum Frequency Generation Spectroscopy. Phys. Chem. Chem. Phys. 2008, 10, 5801–5816. (18) Ishiyama, T.; Imamura, T.; Morita, A. Theoretical Studies of Structures and Vibrational Sum Frequency Generation Spectra at Aqueous Interfaces. Chem. Rev. 2014, 114, 8447–8470. (19) Perry, A.; Neipert, C.; Space, B.; Moore, P. B. Theoretical Modeling of Interface Specific Vibrational Spectroscopy: Methods and Applications to Aqueous Interfaces. Chem. Rev. 2006, 106, 1234–1258. (20) Rey, R.; Møller, K. B.; Hynes, J. T. Ultrafast Vibrational Population Dynamics of Water and Related Systems: A Theoretical Perspective. Chem. Rev. 2004, 104, 1915–1928. (21) Gaigeot, M.-P.; Martinez, M.; Vuilleumier, R. Infrared Spectroscopy in the Gas and Liquid Phase from First Principle Molecular Dynamics Simulations: Application to Small Peptides. Mol. Phys. 2007, 105, 2857–2878. (22) Mennucci, B. Modeling Environment Effects on Spectroscopies Through QM/Classical Models. Phys. Chem. Chem. Phys. 2013, 15, 6583–6594. (23) Zhuang, W.; Hayashi, T.; Mukamel, S. Coherent Multidimensional Vibrational Spectroscopy of Biomolecules: Concepts, Simulations, and Challenges. Angew. Chem., Int. Ed. 2009, 48, 3750–3781. (24) l. C. Jansen, T.; Knoester, J. Waiting Time Dynamics in Two-Dimensional Infrared Spectroscopy. Acc. Chem. Res. 2009, 42, 1405–1411. (25) McRobbie, P. L.; Hanna, G.; Shi, Q.; Geva, E. Signatures of Nonequilibrium Solvation Dynamics on Multidimensional Spectra. Acc. Chem. Res. 2009, 42, 1299–1309.


Page 24 of 31

Page 25 of 31

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


(26) Rhee, H.; Choi, J.-H.; Cho, M. Infrared Optical Activity: Electric Field Approaches in Time Domain. Acc. Chem. Res. 2010, 43, 1527–1536. (27) Kim, H.; Cho, M. Infrared Probes for Studying the Structure and Dynamics of Biomolecules. Chem. Rev. 2013, 113, 5817–5847. (28) Yagasaki, T.; Saito, S. Molecular Dynamics Simulation of Nonlinear Spectroscopies of Intermolecular Motions in Liquid Water. Acc. Chem. Res. 2009, 42, 1250–1258. (29) Yagasaki, T.; Saito, S. Fluctuations and Relaxation Dynamics of Liquid Water Revealed by Linear and Nonlinear Spectroscopy. Ann. Rev. Phys. Chem. 2013, 64, 55–75. (30) Reppert, M.; Tokmakoff, A. Computational Amide I 2D IR Spectroscopy as a Probe of Protein Structure and Dynamics. Ann. Rev. Phys. Chem. 2016, 67, 359–386. (31) Cisneros, G. A.; Wikfeldt, K. T.; Ojamäe, L.; Lu, J.; Xu, Y.; Torabifard, H.; Bartók, A. P.; Csányi, G.; Molinero, V.; Paesani, F. Modeling Molecular Interactions in Water: From Pairwise to Many-Body Potential Energy Functions. Chem. Rev. 2016, 116, 7501–7528. (32) Paesani, F. Getting the Right Answers for the Right Reasons: Toward Predictive Molecular Simulations of Water with Many-Body Potential Energy Functions. Acc. Chem. Res. 2016, in press. (33) Nagata, Y.; Ohto, T.; Backus, E. H. G.; Bonn, M. Molecular Modeling of Water Interfaces: From Molecular Spectroscopy to Thermodynamics. J. Phys. Chem. B 2016, 120, 3785–3796. (34) Sakaguchi, S.; Ishiyama, T.; Morita, A. Theory and Efficient Computation of Differential Vibrational Spectra. J. Chem. Phys. 2014, 140, 144109. (35) Sakaguchi, S.; Ishiyama, T.; Morita, A. Erratum: “Theory and Efficient Computation of Differential Vibrational Spectra” [J. Chem. Phys.140, 144109 (2014)]. J. Chem. Phys. 2014, 141, 149901.



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) Joutsuka, T.; Morita, A. Improved Theory of Difference Vibrational Spectroscopy and Application to Water. J. Chem. Theory Comput. 2016, in press. (37) Walrafen, G. E. Raman Spectral Studies of the Effects of Solutes and Pressure on Water Structure. J. Chem. Phys. 1971, 55, 768–792. (38) Gorbaty, Y. E.; Bondarenko, G. V.; Kalinichev, A. G.; Okhulkov, A. V. The Effect of Pressure on Hydrogen Bonding in Water: IR Study of νOD HDO at Pressures of up to 1500 bar. Mol. Phys. 1999, 96, 1659–1665. (39) Cavaille, D.; Combes, D.; Zwick, A. Effect of High Hydrostatic Pressure and Additives on the Dynamics of Water: a Raman Spectroscopy Study. J. Raman Spectrosc. 1996, 27, 853–857. (40) Inoue, A.; Kojima, K.; Taniguchi, Y.; Suzuki, K. Near-Infrared Spectra of Water and Aqueous Electrolyte Solutions at High Pressures. J. Solution Chem. 1984, 13, 811–23. (41) Khoshtariya, D. E.; Zahl, A.; Dolidze, T. D.; Neubrand, A.; van Eldik, R. Local Dense Structural Heterogeneities in Liquid Water from Ambient to 300 MPa Pressure: Evidence for Multiple Liquid-Liquid Transitions. ChemPhysChem 2004, 5, 1398–1404. (42) Stanley, E. M.; Batten, R. C. Viscosity of Water at High Pressures and Moderate Temperatures. J. Phys. Chem. 1969, 73, 1187–1191. (43) Först, P.; Werner, F.; Delgado, A. The Viscosity of Water at High Pressures – Especially at Subzero Degrees Centigrade. Rheol. Acta 2000, 39, 566–573. (44) Wilbur, D. J.; DeFries, T.; Jonas, J. Self-Diffusion in Compressed Liquid Heavy Water. J. Chem. Phys. 1976, 65, 1783–1786. (45) Harris, K. R.; Newitt, P. J. Self-Diffusion of Water at Low Temperatures and High Pressure. J. Chem. Eng. Data 1997, 42, 346–348. (46) Jonas, J.; Lee, Y. T. NMR and Laser Raman Scattering Studies of Fluids at High Pressure. J. Phys. Condens. Matter 1992, 4, 305. 26 ACS Paragon Plus Environment

Page 26 of 31

Page 27 of 31

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


(47) Okhulkov, A. V.; Demianets, Y. N.; Gorbaty, Y. E. X-ray Scattering in Liquid Water at Pressures of up to 7.7 kbar: Test of a Fluctuation Model. J. Chem. Phys. 1994, 100, 1578–1588. (48) Nosé, S. A Unified Formulation of the Constant Temperature Molecular Dynamics Methods. J. Chem. Phys. 1984, 81, 511–519. (49) Hoover, W. G. Canonical Dynamics: Equilibrium Phase-Space Distributions. Phys. Rev. A 1985, 31, 1695–1697. (50) Melchionna, S.; Ciccotti, G.; Holian, B. L. Hoover NPT Dynamics for Systems Varying in Shape and Size. Mol. Phys. 1993, 78, 533–544. (51) Martyna, G. J.; Tobias, D. J.; Klein, M. L. Constant Pressure Molecular Dynamics Algorithms. J. Chem. Phys. 1994, 101, 4177–4189. (52) Martyna, G. J.; Tuckerman, M. E.; Tobias, D. J.; Klein, M. L. Explicit Reversible Integrators for Extended Systems Dynamics. Mol. Phys. 1996, 87, 1117–1157. (53) Lippert, R. A.; Predescu, C.; Ierardi, D. J.; Mackenzie, K. M.; Eastwood, M. P.; Dror, R. O.; Shaw, D. E. Accurate and Efficient Integration for Molecular Dynamics Simulations at Constant Temperature and Pressure. J. Chem. Phys. 2013, 139, 164106. (54) Tuckerman, M. E.; Liu, Y.; Ciccotti, G.; Martyna, G. J. Non-Hamiltonian Molecular Dynamics: Generalizing Hamiltonian Phase Space Principles to Non-Hamiltonian Systems. J. Chem. Phys. 2001, 115, 1678–1702. (55) Berendsen, H. J. C.; Postma, J. P. M.; van Gunsteren, W. F.; DiNola, A.; Haak, J. R. Molecular Dynamics with Coupling to an External Bath. J. Chem. Phys. 1984, 81, 3684–3690. (56) Morishita, T. Fluctuation Formulas in Molecular-Dynamics Simulations with the Weak Coupling Heat Bath. J. Chem. Phys. 2000, 113, 2976–2982. (57) Paci, E.; Marchi, M. Constant-Pressure Molecular Dynamics Techniques Applied to Complex Molecular Systems and Solvated Proteins. J. Phys. Chem. 1996, 100, 4314–4322. 27 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

(58) Morita, A.; Hynes, J. T. A Theoretical Analysis of the Sum Frequency Generation Spectrum of the Water Surface. II. Time-Dependent Approach. J. Phys. Chem. B 2002, 106, 673–685. (59) Hirschfelder, J. O.; Curtiss, C. F.; Bird, R. B. Molecular Theory of Gases and Liquids; Wiley: New York, 1954. (60) Allen, M. P.; Tildesley, D. J. Computer Simulation of Liquids; Oxford University Press: New York, 1989. (61) Ingram, J. Introductory Statistics; Cummings Publishing Company: Menlo Park, CA, 1974. (62) Radnai, T.; Ohtaki, H. X-ray Diffraction Studies on the Structure of Water at High Temperatures and Pressures. Mol. Phys. 1996, 87, 103–121. (63) Bellissent-Funel, M.-C.; Bosio, L. A Neutron Scattering Study of Liquid D2 O under Pressure and at Various Temperatures. J. Chem. Phys. 1995, 102, 3727–3735. (64) Walrafen, G. Raman Spectra from Partially Deuterated Water and Ice VI to 10.1 kbar at 28◦ C. J. Solution Chem. 1973, 2, 159–171. (65) Sun, Q.; Zheng, H.; Xu, J.; Hines, E. Raman Spectroscopic Studies of the Stretching Band from Water up to 6 kbar at 290 K. Chem. Phys. Lett. 2003, 379, 427–431. (66) Kalinichev, A.; Gorbaty, Y.; Okhulkov, A. Structure and Hydrogen Bonding of Liquid Water at High Hydrostatic Pressures: Monte Carlo NPT -Ensemble Simulations up to 10 kbar. J. Mol. Liq. 1999, 82, 57–72. (67) Scott, J. N.; Vanderkooi, J. M. Evidence of a Structural Defect in Ice VII and the Side Chain Dependent Response of Small Model Peptides to Increased Pressure. Appl. Spectrosc. 2011, 65, 756–764. (68) Tassaing, T.; Danten, Y.; Besnard, M. Infrared Spectroscopic Study of Hydrogen-Bonding in Water at High Temperature and Pressure. J. Mol. Liq. 2002, 101, 149–158. 28 ACS Paragon Plus Environment

Page 28 of 31

Page 29 of 31

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


(69) Hédoux, A.; Guinet, Y.; Paccou, L. Analysis of the Mechanism of Lysozyme Pressure Denaturation from Raman Spectroscopy Investigations, and Comparison with Thermal Denaturation. J. Phys. Chem. B 2011, 115, 6740–6748. (70) Okada, T.; Komatsu, K.; Kawamoto, T.; Yamanaka, T.; Kagi, H. Pressure Response of Raman Spectra of Water and its Implication to the Change in Hydrogen Bond Interaction. Spectrochim. Acta, Part A 2005, 61, 2423–2427. (71) Ishiyama, T.; Morita, A. Analysis of Anisotropic Local Field in Sum Frequency Generation Spectroscopy with the Charge Response Kernel Water Model. J. Chem. Phys. 2009, 131, 244714. (72) McQuarrie, D. Statistical Mechanics; University Science Books: Sausalito, CA, 2000. (73) Schmidt, J. R.; Corcelli, S. A.; Skinner, J. L. Pronounced Non-Condon Effects in the Ultrafast Infrared Spectroscopy of Water. J. Chem. Phys. 2005, 123, 044513. (74) Pimentel, G.; McClellan, A. The Hydrogen Bond; A Series of chemistry books; W.H. Freeman, 1960. (75) Xu, H.; Stern, H. A.; Berne, B. J. Can Water Polarizability Be Ignored in Hydrogen Bond Kinetics? J. Phys. Chem. B 2002, 106, 2054–2060. (76) Morita, A.; Kato, S. Molecular Dynamics Simulation with the Charge Response Kernel: Diffusion Dynamics of Pyrazine and Pyrazinyl Radical in Methanol. J. Chem. Phys. 1998, 108, 6809–6818. (77) Prielmeier, F. X.; Lang, E. W.; Speedy, R. J.; Lüdemann, H.-D. The Pressure Dependence of Self Diffusion in Supercooled Light and Heavy Water. Ber. Bunsen-Ges. Phys. Chem. 1988, 92, 1111–1117. (78) Stern, H. A.; Rittner, F.; Berne, B. J.; Friesner, R. A. Combined Fluctuating Charge and



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

Polarizable Dipole Models: Application to a Five-Site Water Potential Function. J. Chem. Phys. 2001, 115, 2237–2251. (79) Ren, P.; Ponder, J. W. Polarizable Atomic Multipole Water Model for Molecular Mechanics Simulation. J. Phys. Chem. B 2003, 107, 5933–5947.


Page 30 of 31

Page 31 of 31

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 of contents graphic.


Efficient Computation of Difference Vibrational Spectra in Isothermal

Recommend Documents