Combining Static and Dynamical Approaches for Infrared Spectra

Subscriber access provided by EAST TENNESSEE STATE UNIV

Article

Combining static and dynamical approaches for infrared spectra calculations of gas phase molecules and clusters Daria Ruth Galimberti, Alberto Milani, Matteo Tommasini, Chiara Castiglioni, and Marie-Pierre Gaigeot J. Chem. Theory Comput., Just Accepted Manuscript • DOI: 10.1021/acs.jctc.7b00471 • Publication Date (Web): 27 Jun 2017 Downloaded from http://pubs.acs.org on June 28, 2017

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 39

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

Combining static and dynamical approaches for infrared spectra calculations of gas phase molecules and clusters Daria R. Galimberti,∗,†,¶ Alberto Milani,† Matteo Tommasini,† Chiara Castiglioni,∗,† and Marie-Pierre Gaigeot∗,‡ †Dip. Chimica, Materiali, Ing. Chimica "G. Natta", Politecnico di Milano, Milan, Italy ‡LAMBE CNRS UMR8587, Université d’Evry val d’Essonne, Evry, France & Université Paris-Saclay, France ¶Current address: LAMBE CNRS UMR8587, Université d’Evry val d’Essonne, Evry, France & Université Paris-Saclay, France E-mail: [email protected]; [email protected]; [email protected] Phone: +33 (0)169470140

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 Four models for the calculation of the IR spectrum of gas phase molecules and clusters from molecular dynamics simulations are presented with the aim to reduce the computational cost of the usual Fourier transform (FT) of the time correlation function of the dipole moment. These models are based on the VDOS, FT of time correlation function of velocities, and Atomic Polar Tensors (APT). The models differ from each other by the number of APTs inserted into the velocities correlation function. Excellent accuracy is achieved by the model adopting a weighted linear combination of a few selected APTs adapted for the rotation of the molecule (model D). The achieved accuracy relates to band-positions, band-shapes and band-intensities. Depending on the degree of actual dynamics of the molecule, rotational motion, conformational isomerisation and large amplitude motions that can be seen during the finite temperature trajectory, one could also apply one of the other models (models A, B, C ), but with caution. Model D is therefore found simple and accurate, with appealing computational cost and should be systematically applied. Its generalization to condensed phase systems should be straightforward.

Introduction Vibrational spectroscopy (e.g. IR, Raman, SFG...) is one of the most powerful technique applied for the characterization of materials, it is versatile and can be used on a large variety of systems, this being one of the reasons for its widespread application not only in fundamental research but also in the industry. IR spectroscopy has been (and is still) one crucial interrogative tool in chemistry, physics, biology, astrochemistry, astrophysics, and material sciences, where structural characterization is of prime importance, not only for basic research but also for industrial applications. It has been widely used to investigate gas phase molecules, liquids, crystals, semi-crystals, amorphous materials, as well as solid/liquid, solid/solid, solid/gas interfaces, characterizing the architecture of the material and providing 2


Page 2 of 39

Page 3 of 39

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


relationships between (microscopic) structure(s) and (macroscopic) properties that rule its behavior, of interest for technological applications. 1–5 In particular, IR vibrational spectroscopy is certainly one of the most routinely used method in analytical chemistry for the structural characterization of molecules and assemblies, in the gas phase (with action spectroscopies 6–9 ), in liquids 10,11 and at surfaces and interfaces 12–16 (even though vibrational Sum Frequency Generation combining IR and Raman selection rules is definitely more powerful, 17–20 however more complex to interpret). In many applications of vibrational spectroscopy not only qualitative but quantitative measurements are required, for instance to determine the concentration of molecular species. To this aim, it is not only important to understand and interpret the positions of the bands in the IR spectra in terms of molecular species but it is also important to be able to interpret band-intensities and their variations from one sample to another. The interpretation of spectroscopic data based on experimental measurements alone can lead to many uncertainties and ambiguities due to the complexity of numerous phenomena involved. Theoretical calculations of IR spectra have always been used in support to experiments in order to help rationalize positions, intensities and shapes of the spectral bands, and therefore provide the necessary information for a final detailed structural characterization and conformational assignment. In this respect, one can fairly state that it is the predicted position of the bands that is first and foremost used into the spectral assignment process, but predicted intensities are of importance when separating isomeric forms of molecules with subtle structural differences, and are certainly of prime importance whenever a quantitative analysis of the spectra is needed. In this latter instance, a good prediction of band-shapes is also definitely helpful in making final spectral assignments. Because of the remarkable effect of the inter-molecular environment (e.g. specific intermolecular interactions in condensed phases or solvent effects), the accurate prediction of band shapes allows to understand several supramolecular phenomena. There are two families of theoretical methods applied for the calculation of IR vibrational

3



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

spectra (the same discussion applies for Raman): static methods (harmonic and anharmonically corrected) and dynamical methods (which are intrinsically anharmonic). We would like to refer the readers to recent reviews 21–28 for a detailed description of these methods, their respective advantages and disadvantages, this being not the purpose of the present work. The issue we want to raise here is the importance of being able to predict IR vibrational spectra in terms of reliable band-positions, band-intensities and band-shapes, with a theoretical method that is applicable as easily and as successfully for gas phase molecules and clusters, for liquids and solutes immersed in liquids, including biomolecules, and for complex inhomogeneous interfacial systems as solid/liquid interfaces typically probed in electrochemistry. A dynamical approach applied to vibrational spectroscopy is certainly the method of choice responding to these criteria, and simultaneously assuring that temperature and environmental effects are directly included in the spectral calculation without a posteriori models applied. We will briefly present and discuss this methodology in this paper, and we refer the readers to a limited selection of papers 21–23,29–31 for illustrations of dynamical anharmonic IR spectroscopy in the gas phase, liquid phase and in material sciences. The point we want discuss here is the importance of being able to predict accurate IR spectra (and therefore be able to interpret experimental data for conformational assignments) at low enough computational cost. This is the goal of the present paper. The first thing coming into mind would be to run molecular dynamics simulations using classical force fields, cheap to evaluate and therefore providing long enough trajectories at relatively low computational cost. Beside the issue of transferability of classical force fields, it is also well-recognized that these force fields are unable to predict vibrational spectra, unless one uses sophisticated force fields including interaction terms that couple internal motions (e.g. stretch-bending, stretch-torsion, ...) but also including electrostatic terms that not only are going beyond the usual charge-charge interactions (e.g. polarization terms in their various implementations 32–35 ), and furthermore including electrostatic terms that explicitely take care of charge and dipole fluxes. 36 Such fluxes have been shown to be the most important

4


Page 4 of 39

Page 5 of 39

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


ingredients in vibrational spectroscopy analyses. 37,38 Such force fields are extremely rare, as they are complicated to develop and parameterize (already for one given class of molecules, let us not mention transferability issues), and furthermore become computationally rather costly. 39,40 Evidently, all these items point towards using ab initio based molecular dynamics simulations methods that include these interactions and fluxes, by construction. The downside of these methods is their higher computational cost (even for semi-empirical methods). Several groups have shown within the past decade that DFT-based molecular dynamics simulations are powerful for reliable predictions of IR spectra for diverse applications in the gas phase, liquids, materials, nowadays applicable to systems of several hundreds (even thousands 41,42 ) of atoms. Semi-empirical and tight-binding electronic based MD are also found more and more applied to vibrational spectroscopy, simulating even much bigger systems at low computational cost. 43–46 As will be discussed in this paper, dynamical IR spectra are based on the Fourier transformation of dipole time correlation functions, and the strength of the electronic-based MD methods is to calculate these dipole moments without pre-established/parameterized models. This however remains computationally expensive, usually based on the Berry phase method or on the localization of the (delocalized) wavefunction into localized functions as is done with the Wannier center localization. 29,47–49 Alternatives have also been recently developed. 50,51 We present here a method that circumvent any such procedure, without including any a priori parameterization model for the calculation of molecular dipoles. Our strategy is fully based on the ab initio MD trajectories and makes use of previously developed concepts on Atomic Polar Tensor (APT) (both developed for static and dynamical spectra calculations) 52–62 that elegantly include charge and dipole fluxes, and on the well-known property of fast convergence of correlation functions of velocities. The method is described in the first part of this paper, tests and validations are presented afterwards. We restrict this paper to IR spectroscopy of gas phase molecules and clusters as a

5



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 39

first demonstration, including the important topic of spectra calculation when isomerisation events occur along the trajectory. Applications of the same strategy to more complex liquids and materials should be straightforward but is not the object of the present paper.

Background on dynamical IR spectroscopy calculation The usual method for the prediction of the IR spectrum of a system from a molecular dynamics trajectory is based on the Fourier transformation of the dipole autocorrelation function: 63,64 2πβω 2 I(ω) = 3cV

Z

+∞

dt eiωt hδµ(t) · δµ(0)i

(1)

−∞

where β = 1/kT , ω is the frequency of the absorbed light, c is the speed of light in vacuum, V the volume of the system, δµ(t) = µ(t) − hµi is the instantaneous fluctuation of the dipole moment. This has been successfully applied in the gas phase and in the condensed phase in the literature, especially using First-Principles Molecular Dynamics (FPMD) simulations. 21,22,30,65,66 This procedure relies on an accurate determination of the time evolution of the dipole moment, which can be rather computationally involved. One of the issue in this respect is that FPMD programs, such as CP2K 42,67 and CPMD 68–70 codes, work in periodic boundary conditions and the correct definition of the dipole moment of an infinite periodic system is not trivial and unambiguous. Different solutions have been proposed in the literature, the Berry phase approach and the localization of the wave function by means of Wannier functions are most frequently used. The Berry phase method 71,72 allows to predict the dipole moment of the system up to a phase factor, in a relatively cheap computational way. However, only the total dipole of the system can be evaluated, and it is not possible to split the partial contributions, for example arising from a solute from its solvent environment in the condensed phase. Such spliting can be achieved using the localized Wannier function(s) formalism. 29,47–49,73 Within this procedure, the electronic wave function is localized at each step of the dynamics, according to an appropriate transformation. The 6


Page 7 of 39

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


expectation value of the vector position ri of the i-th localized Wannier orbital, is interpreted as the site where an electron pair is located, in this way a set of points with formal electronic charge qi = −2|e| is provided (or qi = −1|e| in the singly occupied orbitals of an open shell system), where |e| is the elementary charge of an electron. After the localization of the electronic wave function, the total dipole moment of the system can be evaluated P P P P as µ = j runs i runs over all the Wannier centers and j Zj Rj where i qi ri + |e|

over all the nuclei (Rj and Zj are respectively the j-th nuclear position and the associated

atomic number). Knowing the details in the positions of the atoms and Wannier centers at each time step of the dynamics, one can easily decompose the total dipole moment of the system µ into individual contributions arising from each molecule of the system µmol . The Wannier localization is however a more expensive method compared to the Berry phase one, and alternative method(s) to eq 1 to calculate the IR spectrum, at lower computational cost but with the same accuracy, is necessary. This is where VDOS spectra, Velocity Density Of States calculated as the Fourier transform of time correlation functions of the velocities, come into play. Velocities are intrinsic components of the trajectories, there is therefore no supplementary effort into their evalutation along time. Furthermore, it is well known that VDOS spectra converge over shorter time-scale trajectories than the ones required for spectra calculated from dipole moment correlations. The downfall is that VDOS spectra do not possess any selection rules, one cannot know which vibrational bands are IR/Raman active/inactive. This information is however encoded in Atomic Polar Tensors (APT) which components are the first derivatives of the dipole with respect to positions in cartesian coordinates. We therefore develop an alternative method for the calculation of IR absorption spectra in which the correlation function of the dipole moment is not used anymore, but the spectrum is predicted by the Fourier transformation of the correlation of atomic velocities (VDOS) convoluted by the correlation of atomic polar tensors P . The starting point is similar to refs 52 and goes beyond these works, proposing a general method.

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

This choice is motivated by the following facts: 1. In the case of static calculations in the double harmonic approximation, the absolute absorption IR intensities are computed from the Atomic Polar Tensor (APT). The APTs are easily available since they are printed in the output of standard calculation packages of harmonic vibrational spectra. 2. The APTs can be easily transferred among different molecules sharing the same chemical groups. 74 Thus it is possible to calculate APTs for small model systems (at a high level of theory), and then transfer the APTs to larger and more complex systems. 3. The APTs have a partially local character: in fact each APT is associated to displacements of a given atom in the molecule. Therefore, it is relatively easy to separately evaluate the contributions due, for instance, to a solute molecule or to the solvent molecules. 4. Using our model, it is possible to decouple the calculation of the velocities correlation function (and hence the vibrational frequencies) from the calculation of the APTs (intensities of bands). On the one hand, VDOS spectra are well described by trajectories much shorter than those required for the prediction of accurate spectra from the correlation function of the dipole moment. On the other hand, accurate band intensities are dependent on the quality of the dipole moments which in turn depend on the quality of the wavefunction of the system, i.e. its calculation with tight enough convergence criteria. A method in which one evaluates the electronic response separately from the atomic velocity response, would allow to use less stringent convergence criteria for the FPMD and would enable predicting IR spectra through shorter trajectories, without loosing accuracy. In the next sections, we illustrate the analytical derivation and some applications of this method.

8


Page 8 of 39

Page 9 of 39

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


Models for IR spectra prediction Within the time-correlation function formalism, 63,64 the IR absorption spectrum of a given system is equal to: 2πβω 2 I(ω) = 3cV

Z

+∞

dt e −∞

iωt

2πβ hδµ(t) · δµ(0)i = 3cV

Z

+∞

dt eiωt hδ( −∞

dµ(t) dµ(0) ) · δ( )i dt dt

(2)

where β = 1/kT , ω is the frequency of the absorbed light, c is the speed of light in vacuum, V the volume of the system, µ(t) is the instantaneous dipole moment of the system at time t, δµ(t) = µ(t) − hµi is the fluctuation with respect to the mean value and h. . . i refers to the equilibrium time correlation function. In equation (2) a quantum correction factor β~/(1 − exp(−β~ω)) has been applied to correct the classical line shape. 29 Equivalently, we write: I(ω)/k1 = where k1 =

2πβ . 3cV

Z

+∞

dt eiωt hδ( −∞

dµ(t) dµ(0) ) · δ( )i dt dt

(3)

For the sake of simplicity in the rest of the text, we will systematically write

µ(t) instead of δµ(t) in the equations, but the fluctuations are indeed taken into account in the calculations. Similarly, fluctuations on the derivatives δ( dµ(t) )= dt

dµ(t) dt

− h dµ i are not dt

written anymore in the equations. If we introduce µu , the u = x, y, z components of µ in the fixed cartesian referential, equation (2) becomes: I(ω)/k1 =

X Z

dt eiωt h

u=x,y,z

dµu (t) dµu (0) · i dt dt

(4)

We now introduce vector ξ that collects the 3N cartesian coordinates of the N atoms of the system: ξ = [x1 , y1 , z1 , x2 , . . . , zN ]T . In the absence of an external field,

dµu (t) dt

can be

expanded as1 3N dµu (t) X ∂µu (t) ∂ξm (t) = · dt ∂ξm ∂t m=1 1

When an external field is present, additional terms must be considered in the expansion.

9


(5)


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

Recognizing that

∂ξm (t) ∂t

= vm (t) is the mth cartesian velocity while

∂µu (t) ∂ξm

Page 10 of 39

is the ’um’ element

of the Atomic Polar Tensor P (t) ( ∂µ∂ξum(t) = Pum (t)), equation (5) becomes: 3N dµu (t) X = Pum (t) · vm (t) dt m=1

(6)

Substituting (6) in (4) provides:

I(ω)/k1 =

X Z

u=x,y,z

=

+∞

dt e

iωt

−∞

3N X 3N Z X X

u=x,y,z m=1 l=1

h

3N X

Pum (t) · vm (t)

m=1

3N X

Pul (0) · vl (0)i

l=1

+∞

dt eiωt h Pum (t)vm (t)Pul (0)vl (0) i

(7)

−∞

While the velocities of all the atoms are easily available at each time step of the trajectory (by construction), it would be computationally too expensive to evaluate also the elements of P (t) at each time step. If the "electrical harmonic" approximation holds, P (t) can be assumed to be time independent, provided a "molecular fixed" reference system is used. This is supported by the fact that the velocity correlation function has fluctuations on short time scales, while the dipole moment derivatives correlation function is almost constant over the same short time scale. This behavior can be exploited to reduce the computational cost for the evaluation of equation (7), by decorrelating the two time-scales. This is the property exploited in the models described below (models A-D). While in principle the same methodology introduced below can be applied to any kind of systems (isolated molecules, condensed phase, ...), for simplicity we will report the equations referring only to one single isolated molecule, and we will demonstrate their validity and limitation on gas phase trajectories only. We present below four models, named model A-D, which aim at reducing the computational cost of equation (7), and we test their validity and domain of applications in the next

10


Page 11 of 39

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


sections. • Model A: one single APT If the trajectory is "not too long" one can assume that the atomic polar tensor does not change over time: Pum (t) ∼ Pum (0)

(8)

This approximation is severe and can hold only at very low temperature or in the solid phase, where neither rigid rotations of the molecule nor conformational isomerisations take place, as indeed it will be shown in our applications in the next sections. • Model B: one single APT adapted for molecular rotation In the gas phase and liquid phase, molecules can rotate along time. Accordingly, the internal reference system of the molecule rotates with respect to the fixed external reference. In order to follow the rotation of the molecule during the trajectory - translations are not considered here since P (t) is translational invariant and the center of mass translation of the molecule is systematically removed from the atomic velocities - one can apply a rotation matrix to Pum (0) (eq. 8) calculated at the initial time t=0, in order to get Pum (t) at any forward time t: P (t) = R(0 → t)P (0)RT (0 → t)

(9)

If we call X(t) the (3N) vector collecting the cartesian coordinates of the N atoms of the molecule in the center of mass reference, R(0 → t) is the rotation matrix2 describing the best overlap between R(0 → t)X(0) and X(t) (RT is the transpose). To evaluate R, we used a quaternion fit code 75,76 that minimizes the sum of the squared distances between the mass weighted coordinates of corresponding atoms. It has been demonstrated that such a rotation satisfies the Eckart axis conditions in case of small displacements. 77 In the case of flexible molecules and/or to avoid the noise generated by fast vibrating atoms (like the H atoms of 2

R(0 → t) is a (3N,3N) diagonal block matrix: each block is the (3,3) rotational matrix for a single atom.

11



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 39

a CH3 rotor), the use of an appropriate subset of atoms for the definition of the positions matrix X, can be a smart choice for the quaternion fit (for instance by not selecting the H atoms of the CH3 rotor into the subset). Another possible strategy to determine R(0 → t) is to consider P (0) in the principal axis reference of the molecule. Hence, R(0 → t) can then be estimated as the rotation matrix which relates the principal axes of the molecule at time t to the reference principal axes of the molecule. It is important to recall that both models A and B do not consider the effects of isomerisation of the molecules. However it is reasonable to assume that, for small displacements (deformations due to vibrational motions only), this error could be negligible for the quantity P (t). We will see in the applications whether this hypothesis holds. • Model C: time dependent APTs possibly adapted for molecular rotation While the method proposed in Model B can account for the rotational motion of the molecule as a rigid body, it cannot describe its conformational changes. Considering as example the Ala − Ala − H + molecule (discussed in more details in section Results), this peptide has two stable conformers at room temperature separated by a low energy barrier. During room temperature trajectories both conformations are explored, as shown in ref., 78 and the APT of the two isomers are inherently different. To take into consideration the conformational dynamics during the trajectory, one can define P (t) within time intervals ∆τi . To this end, the trajectory is divided into time intervals ∆τi = (tb − ta ) (∆τi >> δt), and for each ∆τi interval, one assumes: Pum (t) = Pum (ta )

if

ta < t < tb

(10)

Hence, for each ∆τi interval one needs to compute the APT once (i.e., Pum (ta )). The maximum length of the ∆τi intervals required to take into account the variations of P (t) is strongly dependent on the dynamics of the system under study. Slow evolving systems, like rigid molecules at low temperature (see the cis-NMA example in the section Results), will require ∆τi of the order of the picosecond, while other systems with a fast dynamics might 12


Page 13 of 39

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


demand the APTs to be evaluated each 100 fs. Furthermore, obviously the transformation described in model B can also be applied here: during each ‘macro’ interval ∆τi , one can rotate P (ta ) at each time step in order to follow the “rigid” rotation of the molecule.

P (t) = R(ta → t)P (ta )RT (ta → t)

if

ta < t < tb

(11)

In many cases we have seen that this allows to greatly increase the time period ∆τi (i.e. less APTs to compute) without loss of accuracy (see for example the trans − N M A and the Ala − Ala − H + IR spectra in section Results). Adapting the APT to the molecule rotation will be the standard strategy for model C as shown in section Results. We will explicity note when we will not apply it, with the label ’NoRot’. • Model D: weighted linear combination of few selected APTs adapted for molecular rotation With model C, using the proper time interval ∆τi between the evaluation of two subsequent APTs, one should be able to get a reliable description of changes in the molecular conformation or large amplitude vibrations. However, there are cases (for example in presence of fast rotating CH3 groups) in which the required time interval would be too small and the calculation of too many APTs would be required, as will be shown in the section Results. To avoid this issue, we propose a further improvement in the evaluation of P (t). Considering again a single molecule that can change its conformation several times over the trajectory (e.g. due to isomerisation, CH3 rotation, proton transfer, large amplitude motions, very flexible conformations, · · · ), a reasonable hypothesis is that P (t) at instant t can be properly described by a linear combination of the atomic polar tensors Pj ref of a set of appropriate

13



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 39

reference structures {j}:

P (t) =

X

wj (t)Rj (ref → t)Pjref (ref )Rj T (ref → t)

(12)

j

where wj (t) is the weight of the j-th reference structure in the whole combination at time t. The rotation matrix R(ref → t) must be introduced in the equation to ensure the consistency of the internal reference system of the {j} structures with the one of the molecule at time t along the trajectory (as in models B and C). The reference structures can be chosen following different strategies. For example it is possible to use the equilibrium conformations of the molecule, if available, or these stable conformations with the addition of some transition points along the transition pathway (if known), see for instance the Cl− · · · methanol example treated in section Results. Another possibility is to randomly choose some structures explored along the trajectory and to use them as reference. To evaluate the weights wj in the combination, we followed a strategy similar to the one proposed by Mathias et. al 79 in a different context. We assumed a Gaussian distribution around reference structures: d

exp(− 2σjc ) wj (t) = P dj j exp(− 2σc ) where dj is a properly defined ’distance’ measuring the difference between the geometry of the molecule at time t and the reference structure j, and σc is the width of the Gaussian distribution. To define dj , one possibility is to directly use cartesian coordinates, i.e., dj = ||R(ref → t)X ref − X(t)||. However, this choice can generate numerical noise, because it j includes changes over all the internal coordinates, such as vibrational displacements along bonds and valence angles, which should be rather independent on conformational changes. To avoid this problem, we introduce a distance based on selected internal coordinates. In particular, we have chosen to define the distance based only on torsional angles (this can be anyway easily generalized to any combination of internal coordinates, as proposed by

14


Page 15 of 39

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


Mathias et. al . 79 ): v u Nt uX dj (t) = t (θk (t) − θkj ) k

where θk (t) is the k-th torsional angle of the molecule at step t of the trajectory, and θkj is the corresponding torsional angle value for the reference structure j. In most of the cases a small number Nt of torsional angles (even only one) allows to univocally characterize the different reference structures.

Results Below we benchmark models A-D for the prediction of IR spectra of gas phase molecules and clusters, taking as a reference the spectra obtained from the correlation function of the dipole moment (equation 2). To this aim, we used the same FPMD simulations to predict the spectra both from the dipole correlation function and from our simplified models. The systems tested were selected among some of our previously published applications of DFTMD for IR spectroscopy in the gas phase, for which the IR spectrum from the "standard" FPMD method (equation 2), nicely reproduced the experimental spectra. All the FPMD have been generated with the CP2K 42,67 code. The computational details of the trajectories for each case are reported in each section below. The evaluation of the APTs has been done by means of the Gaussian code. 80 For consistency with the FPMD simulations, the functional BLYP 81,82 has been adopted, together with the 6-311++G** basis set, chosen as a compromise between the required accuracy 83 and the need to limit the computational cost (but the level of theory could be augmented with reasonable additional effort). Note that dispersion corrections were not included as they were shown in our previous papers for the systems of interest in the present work, not to be of importance for the spectra. The four examples chosen illustrate the different accuracy/relevance of models A to D, increasing the complexity of the examined system: 1) low temperature dynamics, 2) rotation of the CH3

15



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 39

group in the molecule, 3) isomeric conformational dynamics, 4) large amplitude torsional vibrational motion. Note that each time a discussion on ’computational cost’ will be reported, it systematically refers to the cost of calculating APT tensors only, since the MD trajectory is common to all models tested. Eq. 7 generates the whole absorption spectrum, which in principle includes both the contribution of vibrational and rotational transitions. Often the rotational subspace is not well sampled along trajectories, which can lead to artifacts in the spectrum. For this reason, it could be advisable to remove rotational contributions from the simulated spectra. However, in all the cases analyzed here, the molecules have a sufficiently high moment of inertia, which ensures that rotational artifacts do not appear in the simulated IR absorption spectrum. For this reason, for all the cases presented, it has not been necessary to introduce any correction aimed at the removal of the rotational contribution.

1- IR spectrum of the cis isomer of N-Methylacetamide (NMA) at 20 K At 20 K, cis-NMA only vibrates around its equilibrium geometry and conformational changes Table 1: Information on the FPMD trajectory of the cis isomer of the NMethylacetamide (NMA) molecule. Trajectory information: Time step: 0.4 fs Trajectory length: 10 ps T=20K BLYP/TZV2PX-MOLOPT-GHT + PW energy cutoff 340Ry REF. 21,84

are inhibited. 21,84 Therefore, as shown in Figure 1, our model for the IR spectrum in its simplest version (model A) already predicts a spectrum practically superimposable to the one evaluated with the dipole correlation function. Model A is therefore an excellent choice 16


Page 17 of 39

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 the IR spectrum calculation of such a ’non-dynamical’ molecule.

Figure 1: Computed infrared spectrum of the N-Methylacetamide (NMA) cis-isomer at 20K. Bottom: reference spectrum from equation 2. Top: spectrum from model A (equation 8.)

2- IR spectrum of the trans isomer of N-Methylacetamide (NMA) at 300 K Table 2: Information on the FPMD trajectory of the trans isomer of the NMethylacetamide (NMA) molecule. Trajectory information: Time step: 0.4 fs Trajectory length: 10ps T=300K BLYP/TZV2PX-MOLOPT-GHT+ PW energy cutoff 340Ry REF. 21,84

Along the 300 K trajectory, the NMA molecule does not isomerize but it rotates as a rigid 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 ACS Paragon Plus Environment

Page 18 of 39

Page 19 of 39

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


modes of the NMA molecule is very sensitive to the effects of the rotation of the whole molecule and less sensitive to the torsion of the CH3 groups. Therefore, as demonstrated by the IR spectrum calculated from model B, this band is better described by using one single APT rotated to follow the trajectory (model B), rather than by using multiple non-rotated APTs obtained every ∆τi = 100 fs along the trajectory (i.e., model C with no rotation applied). Indeed, considering the temperature, the molecule rotates in a non-negligible way within a time interval of 100 fs. A similar behavior is found for band (b) at 1650 cm−1 (amide I band). Its intensity is underestimated in all the cases where the APT rotation (i.e., B type correction) is not used. Let us now consider the group of bands between 1600 and 1300 cm−1 (c), related to the CH3 bending motion. Model B is necessary for a proper description of the (c) bands in the spectrum. However, the relative intensities of this group of bands is better predicted by using multiple APTs to calculate the spectrum, as in model C. This is due to the large amplitude torsional motions of the CH3 groups, which are coupled with the bending modes. Similar arguments also apply to the convolution of bands at 1200 cm−1 (d), which is associated to the NH-bending, mechanically coupled to the torsional motions of the CH3 groups. However, while (c) is not very sensitive to the details in the torsional motion of the CH3 (just one calculated APT every ∆τi = 1 ps is enough), a correct description of the shape of (d) requires at least the evaluation of one APT every ∆τi = 100 fs, although an acceptable prediction is already accomplished with model B alone. Let us finally consider the high frequencies region, i.e. above 2800 cm−1 . The NH stretching band at 3500 cm−1 (e) is decoupled from the torsional motions of the CH3 groups. Therefore it is already well described by using just one APT, provided that the APT rotation correction is applied (model B). On the other hand, the CH-stretching modes (bands between 3100 - 2800 cm−1 ) are highly coupled to the CH3 torsional motions and only model C with an APT calculated every ∆τi = 100 fs can predict the correct shape and details of these bands.

19



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

As a conclusion, we can state that model B taking care of the rotational motion of the APT tensors along time is mandatory to apply, while, depending on the vibrational mode and its coupling to the CH3 torsional motions, either one unique APT (model B) or multiple APTs (model C) could be used. In this latter case, the time interval for the evaluation of the APTs can range from 1 ps to 100 fs or less, depending on the vibrational transition of interest.

3- IR spectrum of the protonated Alanine Dipeptide Ala-Ala-H+ , trans-conformers at 300 K Table 3: Information on the FPMD trajectory of trans conformers of the protonated Alanine Dipeptide (Ala-Ala-H+ ) Trajectory information: Time step: 0.5 fs Trajectory length: 15 ps T=300K BLYP/DZVP-GHT + PW energy cutoff 320Ry REF. 21,78,85

As shown in ref 78, at room temperature this molecule shows a continuous isomerization between two conformers, labelled as TransA1 and TransA2. These two conformers are characterized by a different value of the torsional angle τ (-C-N-C-COOH), namely τ = 201.5o for TransA1 and τ = 279.5o for TransA2 (from Gaussian BLYP/6-311++G** geometry optimization). This continuous isomerization is due to the fact that the (harmonic) potential energy surface has a small barrier between TransA1 and TransA2 (see Figure 3) and that the free energy profile extracted from the dynamics (including anharmonic effects) is barrierless. 78 Moreover, while the energy difference between the two conformers at 0 K is about 1.5 kcal/mol (Gaussian calculations BLYP/6-311++G**), entropic contributions at 300 K lower the free energy difference to 0.3 kcal/mol (see ref. 78 for details). Before going into 20


Page 20 of 39

Page 21 of 39





Page 22 of 39

Page 23 of 39

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


Let us focus now on the group of bands between 1600 and 1700 cm−1 . The shoulder at 1611 cm−1 (b), related to the bending of the CH3 group bonded to the carbon C* (see scheme in table 3), is slightly sensitive to the rotation of the molecule and it is not sensitive to isomerization; for this reason it is predicted with almost the same shape and intensity by all the approximations. On the other hand, the two bands (c) at 1714 cm−1 (C = O stretching of the acid group), and 1663 cm−1 (backbone amide I band), are predicted with the right shape in all cases, but their predicted intensity is generally weaker than in the reference spectrum. In particular, model C and a calculation of the APTs every ∆τi = 250 fs or every 1 ps provides the correct intensity ratio with respect to the group of bands in the region below 1500 cm−1 . This is due to the fact that these bands are sensitive to the changes in internal conformation of the peptide and their consequences on the APT contributions. The 1600-1500 cm−1 region (d), related to the NH bending of both amide groups and the N H 3+ terminal, gives a further demonstration of the importance of the B model. Indeed, this region is well-predicted just by using model B and one APT, while it is too much structured when employing the C model, also by using APTs computed every ∆τi = 250 fs. The same behavior is also shown by the two bands (e) between 1200 cm−1 and 1100 cm−1 . Since these bands are due to NH and CH bending of the N H 3+ and CH3 groups bonded to carbon C*, they are not much affected by isomerization. Another example of this behavior is the doublet at 700 cm−1 (f), related to skeleton vibrations (in particular torsions around the amide group), which is too intense if the B-type model is not adopted. Considering the structured region between 1500 cm−1 and 1300 cm−1 , related to CH and NH bendings, we found that this group of bands is very sensitive to both the rotation of the molecule and conformational changes. Therefore, the spectral pattern and intensity ratios with respect to the other regions of the spectrum can be predicted reliably only by using the C model and APTs computed each ∆τi = 250 fs. Similar considerations can be applied to the band at 400 cm−1 (h), related to skeleton vibrations.

23




Page 24 of 39

Page 25 of 39

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


it, at a cost which is 1/25 (i.e. it uses 2 APTs instead of 50 !) This finding indicates that model D, when an appropriate number of reference geometries is used (only two here), can correctly take into account the evolution in the conformation of the molecule over time and its influence on the spectral patterns, at a much lower computational cost than model C.

4- Cl− · · · methanol ionic complex at 100K Table 4: Information on the FPMD trajectory of Cl− · · · methanol ionic complex Trajectory information: Time step: 0.4 fs Trajectory length: 4.5 ps T=100K BLYP/aug-TZV2P-GHT + PW energy cutoff 360Ry REF. 86

The Cl− · · · methanol cluster allows to assess further the ability of model D to predict the spectrum of a molecule in the presence of a large torsional motion, as the CH3 group rotates practically freely around the CO bond at 100 K. 86 Also the low frequency mode related to the Cl− · · · H intermolecular stretching of the H-bond can be mechanically coupled with other vibrations. We tested two possible choices for the D model: i) Using only the APTs corresponding to the three minima positions on the torsional potential energy surface of the CH3 group, labeled hereafter model D(min). ii) Including also the contribution of transition geometries between two minima (maximum of the torsional potential energy surface), labeled hereafter model D (min + max) As reported in Figure 6, even by employing model C with a frequency of change in APTs every ∆τi = 100 fs, the calculated spectrum shows spurious bands compared to the reference spectrum obtained from the dipole correlation. The doublet around 750 cm−1 (a) in the spectrum given by model C is described as one single band in the reference spectrum. This band, related to the CH3 torsional motion, is well described by model D(min), which includes 25




Page 26 of 39

Page 27 of 39

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


molecule), is very sensitive to the description of the CH3 torsional contributions due to mechanical coupling. Therefore this region can be described satisfactorily only in model D by using both the minimum and the maximum points. However, even in this case, the model does not perfectly reproduce the reference spectrum, and presumably a few more selected conformations into the weighted combination are necessary. The results on the Cl− · · · methanol cluster together with those on Ala-Ala-H+ indicate that the D model allows to predict an accurate spectrum, both in flexible molecules characterized by more than one stable conformation explored over time, and/or where large amplitude vibrational motions occur. The key point is considering enough geometries to compute the different APTs required for an accurate sampling of the dipole derivatives along the whole trajectory.

Conclusions We have presented simplified and low computational cost models to compute IR spectra from the Fourier transform (FT) of the correlation function of velocities (characterized by fast convergence, i.e. allowing short MD trajectories) modulated by (suitable) APT (Atomic Polar Tensor) contributions. This approach is alternative to the FT of the correlation function of the total dipole moment of the system. Different strategies for selecting APTs along time have been proposed to make the calculation both computationally feasible (low computational cost) and accurate. Hence four models have been presented and assessed on different gas phase and cluster molecular systems. In particular, model D (as denoted in the text), which is based on the weighted linear combination of few selected APT tensors adapted for the rotational motion of the molecule along time, allows to obtain a spectrum which is in excellent agreement (in terms of band positions, shapes and intensities) with the reference spectrum obtained through the correlation function of the dipole moment. For instance, compared to the method proposed in ref., 52 which required the evaluation of a new set of APTs

27



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

at every ∆ti = 1 fs, our model remarkably allows predicting a very accurate IR spectrum of the highly flexible protonated Ala-Ala-H+ trans-peptide by using just two APTs. Similarly, the IR spectrum of Cl− · · · methanol cluster, and especially the bands at low frequency related to large amplitude intermolecular motions, have been obtained with high accuracy by model D, including only six selected APT tensors. Interestingly, the Cl− · · · methanol cluster shows that also the transition states between two stable conformations may give important contributions and model D naturally includes these. The cost of model D is very limited, as we have shown that it is the cheaper and the more accurate method (between the model tested). The methods developed here have been applied to gas phase IR spectroscopy. We are currently extending these concepts to other vibrational spectroscopies for more complex and inhomogeneous molecular systems, such as solids, liquids, or solid/liquid interfaces. Such generalization should be straightforward although not so trivial either, as the APTs should also include the effects of the molecular environment in the condensed phase. The concepts introduced here (i.e. a prefactor quantity obtained from high level quantum chemistry calculation, weighted in time, and time-correlated with velocities or positions) should also be easily transferred to the calculation of other observables from molecular dynamics, reducing their computational cost. Other theoretical considerations should also be assessed, such as the length of the trajectory needed to reliably calculate vibrational spectra. Also, can we relax the accuracy of the electronic representation in the first-principles MD trajectories (e.g., Gaussian basis set size, plane-wave basis set size, convergence cut-off values for the wave function optimization, etc.) and still obtain reliable velocities, therefore still providing accurate band-positions ? Considering a very precise force field developed for vibrational spectroscopy, could we even use classical trajectories for the velocities and (high) level quantum APTs and get the same accuracy for spectra ? These are typical issues opened by the present work.

28


Page 28 of 39

Page 29 of 39

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 This work was performed using HPC resources from GENCI-France (CINES/IDRIS Grant 072484). MPG also acknowledges that part of this work was developed during her period at Institut Universitaire de France (IUF).

References (1) Deglmann, P.; Schafer, A.; Lennartz, C. Application of quantum calculations in the chemical industry. An overview. Int. J. Quant. Chem. 2015, 115, 107–136. (2) Castiglioni, C. Handbook of Vibrational Spectroscopy; John Wiley & Sons, Ltd, 2006. (3) Zerbi, G. Molecular Vibrations of High Polymers. Applied Spect. Rev. 1969, 2, 193–261. (4) Painter, P. C.; Coleman, M. M.; Koenig, J. L. Theory of vibrational spectroscopy and its application to polymeric materials; Wiley, 1982. (5) Coleman, M. M.; Painter, P. C.; Graf, J. F. Specific interactions and the miscibility of polymer blends; CRC Press, 1995. (6) Schermann, J. P. Spectroscopy and Modeling of Biomolecular Building Blocks; Elsevier Science, Amsterdam, The Netherlands, 2007. (7) Rizzo, T. R.; Stearns, J. A.; Boyarkin, O. V. Spectroscopic studies of cold, gas-phase biomolecular ions. Int. Rev. Phys. Chem. 2009, 28, 481–515. (8) Rijs, A. M.; Oomens, J. Gas-Phase IR Spectroscopy and Structure of Biological Molecules. Topics Curr. Chem. 2015, 364, 1–42. (9) Heine, N.; Asmis, K. Cryogenic ion trap vibrational spectroscopy of hydrogen-bonded clusters relevant to atmospheric chemistry. Int. Rev. Phys. Chem. 2015, 34, 1–34.

29



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

(10) Reppert, M.; Tokmakoff, A. Computational Amide I 2D IR Spectroscopy as a Probe of Protein Structure and Dynamics. Annu. Rev. Phys. Chem. 2016, 67, 359–386. (11) Kim, H.; Cho, M. Infrared Probes for Studying the Structure and Dynamics of biomolecules. Chem. Rev. 2013, 113, 5817–5847. (12) Costa, D.; Pradier, C.; Tielens, F.; Savio, L. Adsorption and self-assembly of bioorganic molecules at model surfaces: A route towards increased complexity. Surf. Sci. Rep. 2015, 70, 449–553. (13) Hamamoto, M.; Katsura, M.; Nishiyama, N.; Tononue, R.; Nakashima, S. Transmission IR Micro-Spectroscopy of Interfacial Water between Colloidal Silica particles. Surf. Sci. Nanotech. 2015, 13, 301–306. (14) Davantes, A.; Lefevre, G. Molecular orientation of molybdate ions adsorbed on goethite nanoparticles revealed by polarized in situ ATR-IR spectroscopy. Surf. Sci. 2016, 653, 88–91. (15) Tang, M.; Cziczi, D.; Grassian, V. Interactions of Water with Mineral Dust Aerosol: Water Adsorption, Hygroscopicity, Cloud Condensation, and Ice Nucleation. Chem. Rev. 2016, 116, 4205–4259. (16) Griffith, E.; Guizado, T.; Pimentel, A.; Tyndall, G.; Vaida, V. Oxidized Aromaticaliphatic mixed films at the air-aqueous solution interface. J. Phys. Chem. C. 2013, 117, 22341–22350. (17) Shen, Y. R.; Ostroverkhov, V. Sum-Frequency Vibrational Spectroscopy on Water Interfaces:Polar Orientation of Water Molecules at Interfaces. Chem. Rev. 2006, 106, 1140–1154. (18) Tian, C.; Shen, Y. R. Recent progress on sum-frequency spectroscopy. Surf. Sci. Rep. 2014, 69, 105–131. 30


Page 30 of 39

Page 31 of 39

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


(19) Jubb, A.; Hua, W.; Allen, H. Environmental Chemistry at Vapor/Water Interfaces: Insights from Vibrational Sum Frequency Generation Spectroscopy. Annu. Rev. Phys. Chem. 2012, 63, 107–130. (20) Covert, P.; Hore, D. Geochemical Insight from Nonlinear Optical Studies of mineralwater interfaces. Annu. Rev. Phys. Chem. 2016, 67, 233–257. (21) Gaigeot, M.-P. Theoretical spectroscopy of floppy peptides at room temperature. A DFTMD perspective: gas and aqueous phase. Phys. Chem. Chem. Phys. 2010, 12, 3336–3359. (22) Thomas, M.; Brehm, M.; Fligg, R.; Vohringer, P.; Kirchner, B. Computing vibrational spectra from ab initio molecular dynamics. Phys. Chem. Chem. Phys. 2013, 15, 6608– 6622. (23) Baldauf, C.; Rossi, M. Going clean: structure and dynamics of peptides in the gas phase and paths to solvation. J. Phys.: Condens. Matt 2015, 27, 493002. (24) Gaigeot, M. P.; Spezia, R. Theoretical Methods for Vibrational Spectroscopy and Collision Induced Dissociation in the Gas Phase. Topics Curr. Chem. 2015, 364, 99–151. (25) Roy, T. K.; Gerber, R. B. Vibrational self-consistent field calculations for spectroscopy of biological molecules: new algorithmic developments and applications. Phys. Chem. Chem. Phys. 2013, 15, 9468–9492. (26) Chaban, G. M.; Gerber, R. B. Anharmonic vibrational spectroscopy calculations with electronic structure potentials: comparison of MP2 and DFT for organic molecules. Theor. Chem. Account. 2008, 120, 273–279. (27) Qu, C.; Prosmiti, R.; Bowman, J. MULTIMODE calculations of the infrared spectra of H7+ and D7+ using ab initio potential energy and dipole moment surfaces. Theor. Chem. Acc. 2013, 132, 1413. 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

(28) Barone, V.; Biczysko, M.; Bloino, J. Fully anharmonic IR and Raman spectra of medium size molecular systems: accuracy and interpretation. Phys. Chem. Chem. Phys. 2014, 16, 1759–1787. (29) Gaigeot, M. P.; Sprik, M. Ab Initio Molecular Dynamics Computation of the Infrared Spectrum of Aqueous Uracil. J. Phys. Chem. B. 2003, 107, 10344–10358. (30) Sulpizi, M.; Salanne, M.; Sprik, M.; Gaigeot, M.-P. Vibrational Sum Frequency Generation Spectroscopy of the Water LiquidâĂŞVapor Interface from Density Functional Theory-Based Molecular Dynamics Simulations. J. Chem. Phys. Let. 2013, 4, 83–87. (31) Gaigeot, M. P.; Sulpizi, M. In Molecular modeling of geochemical reactions: an introduction - Chapter 8: Mineral-water interaction; Kubicki, J., Ed.; 2016; Chapter 8, pp 271–310. (32) Ren, P.; Wu, C.; Ponder, J. W. Polarizable Atomic Multipole-Based Molecular Mechanics for Organic Molecules. J. Chem. Theory. Comput. 2011, 7, 3143–3161. (33) Chelli, R.; Procacci, P. A transferable polarizable electrostatic force field for molecular mechanics based on the chemical potential equalization principle. J. Chem. Phys. 2002, 117, 9175–9189. (34) Lopes, P. E. M.; Huang, J.; Shim, J.; Luo, Y.; Li, H.; Roux, B.; MacKerell, A. D. Polarizable Force Field for Peptides and Proteins Based on the Classical Drude Oscillator. J. Chem. Theory. Comput. 2013, 9, 5430–5449. (35) Esser, A.; Belsare, S.; Marx, D.; Head-Gordon, T. Mode specific THz spectra of solvated amino acids using the AMOEBA polarizable force field. Phys. Chem. Chem. Phys. 2017, 19, 5579–5590. (36) Palmo, K.; Krimm, S. Electrostatic model for IR intensities in a spectroscopically determined molecular mechanics force field. J. Comput. Chem. 1998, 19, 754–768. 32


Page 32 of 39

Page 33 of 39

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


(37) Galimberti, D.; Milani, A.; Castiglioni, C. Charge mobility in molecules: Charge fluxes from second derivatives of the molecular dipole. J. Chem. Phys. 2013, 138, 164115. (38) Galimberti, D.; Milani, A.; Castiglioni, C. Infrared intensities and charge mobility in hydrogen bonded complexes. J. Chem. Phys. 2013, 139, 074304. (39) Semrouni, D.; Sharma, A.; Dognon, J.; Ohanessian, G.; Clavaguera, C. Finite Temperature Infrared Spectra from Polarizable Molecular Dynamics Simulations. J. Chem. Theory. Comput. 2014, 10, 3190–3199. (40) Kratz, E.; Walker, A.; Lagardere, L.; Lipparini, F.; Piquemal, J.; Cisneros, G. LICHEM: A QM/MM Program for Simulations with Multipolar and Polarizable Force Fields. J. Comput. Chem. 2016, 37, 1019–1029. (41) VandeVondele, J.; Borstnik, U.; Hutter, J. Linear Scaling Self-Consistent Field Calculations with Millions of Atoms in the Condensed Phase. J. Chem. Theory. Comput. 2012, 8, 3565–3573. (42) Hutter, J.; Iannuzzi, M.; Schiffmann, F.; VandeVondele, J. CP2K: atomistic simulations of condensed matter systems. Wiley Interdisciplinary Reviews: Computational Molecular Science 2014, 4, 15–25. (43) Simon, A.; Iftner, C.; Mascetti, J.; Spiegelman, F. Water Clusters in an Argon Matrix: Infrared Spectra from Molecular Dynamics Simulations with a Self-Consistent Charge Density Functional Based Tight Binding/Force-Field Potential. J. Phys. Chem. A. 2015, 119, 2449–2467. (44) Oliveira, L.; Cuny, J.; Moriniere, M.; Dontot, L.; Simon, A.; Spiegelman, F.; Rapacioli, M. Phase changes of the water hexamer and octamer in the gas phase and adsorbed on polycyclic aromatic hydrocarbons. Phys. Chem. Chem. Phys. 2015, 17, 17079–17089.

33



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

(45) Meng, Y.; Wu, Q.; Chen, L.; Wangmo, S.; Gao, Y.; Wang, Z.; Zhang, R.; Ding, D.; Niehaus, T.; Frauenheim, T. Signatures in vibrational and UV-visible absorption spectra for identifying cyclic hydrocarbons by graphene fragments. Nanoscale 2013, 5, 12178–12184. (46) Feng, C.; Zhang, R.; Dong, S.; Niehaus, T.; Frauenheim, T. Signatures in Vibrational Spectra of Ice Nanotubes Revealed by a Density Functional Tight Binding Method. J. Phys. Chem. C. 2007, 111, 14131–14138. (47) Marzari, N.; Vanderbilt, D. Maximally localized generalized Wannier functions for composite energy bands. Phys. Rev. B 1997, 56, 12847–12865. (48) Silvestrelli, P. L.; Parrinello, M. Water Molecule Dipole in the Gas and in the Liquid Phase. Phys. Rev. Let. 1999, 82, 3308–3311. (49) Silvestrelli, P. L.; Parrinello, M. Structural, electronic, and bonding properties of liquid water from first principles. J. Chem. Phys. 1999, 111, 3572–3580. (50) Thomas, M.; Brehm, M.; Kirchner, B. Voronoi dipole moments for the simulation of bulk phase vibrational spectra. Phys. Chem. Chem. Phys. 2015, 17, 3207–3213. (51) Luber, S. Local electric dipole moments for periodic systems via density functional theory embedding. J. Chem. Phys. 2014, 141, 234110. (52) 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. (53) Person, W. B.; Zerbi, G. Vibrational intensities in infrared and Raman spectroscopy; Elsevier Science Ltd, 1982; Vol. 20. (54) Person, W. B.; Newton, J. H. Dipole moment derivatives and infrared intensities. I. Polar tensors. J. Chem. Phys. 1974, 61, 1040–1049. 34


Page 34 of 39

Page 35 of 39

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


(55) Castiglioni, C.; Gussoni, M.; Zerbi, G. Handbook of Vibrational Spectroscopy, edited by J. Chalmers and P. Griffiths; John Wiley and Sons, Chichester, UK, 2001. (56) Decius, J. C. An effective atomic charge model for infrared intensities. J. Mol. Spect. 1975, 57, 348–362. (57) King, W. T.; Mast, G. B. Infrared intensities, polar tensors, and atomic population densities in molecules. J. Phys. Chem. 1976, 80, 2521–2525. (58) Gussoni, M.; Castiglioni, C.; Ramos, M.; Rui, M.; Zerbi, G. Infrared Intensities - from Intensity Parameters to an Overall Understanding of the Spectrum. J. Mol. Struct. 1990, 224, 445–470. (59) Haiduke, R. L. A.; Bruns, R. E. An atomic charge-charge flux-dipole flux atom-inmolecule decomposition for molecular dipole-moment derivatives and infrared fundamental intensities. J. Phys. Chem. A 2005, 109, 2680–2688. (60) Milani, A.; Castiglioni, C. Modeling of Molecular Charge Distribution on the Basis of Experimental Infrared Intensities and First-Principles Calculations: The Case of CH Bonds. J. Phys. Chem. A 2010, 114, 624–632. (61) Milani, A.; Galimberti, D.; Castiglioni, C.; Zerbi, G. Molecular charge distribution and charge fluxes from Atomic Polar Tensors: The case of OH bonds. J. Mol. Struct. 2010, 976, 342–349. (62) Milani, A.; Tommasini, M.; Castiglioni, C. Atomic charges from IR intensity parameters: theory, implementation and application. Theo. Chem. Acc. 2012, 131, 1139. (63) McQuarrie, D. A. Statistical Mechanics; University Science Books, 2000. (64) Kubo, R.; Toda, M.; Hashitsume, N. Statistical Physics II ; Springer Series in SolidState Sciences; Springer Berlin Heidelberg, 1991; Vol. 31.

35



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

(65) Silvestrelli, P. L.; Bernasconi, M.; Parrinello, M. Ab initio infrared spectrum of liquid water. Chem. Phys. Letters 1997, 277, 478 – 482. (66) Pagliai, M.; Cavazzoni, C.; Cardini, G.; Erbacci, G.; Parrinello, M.; Schettino, V. Anharmonic infrared and Raman spectra in CarâĂŞParrinello molecular dynamics simulations. J. Chem. Phys. 2008, 128, 224514. (67) VandeVondele, J.; Krack, M.; Mohamed, F.; Parrinello, M.; Chassaing, T.; Hutter, J. Quickstep: Fast and accurate density functional calculations using a mixed Gaussian and plane waves approach. Computer Physics Communications 2005, 167, 103–128. (68) Marx, D.; Hutter, J. In Modern Methods and Algorithms of Quantum Chemistry; Grotendorst, J., Ed.; NIC; Forschungszentrum Julich, 2000; Chapter 13, pp 301–449. (69) Andreoni, W.; Curioni, A. New Advances in Chemistry and Material Science with CPMD and Parallel Computing. Parallel Computing 2000, 26, 819–842. (70) Marx, D.; Hutter, J. Ab Initio Molecular Dynamics; Cambridge University Press, 2009. (71) King-Smith, R. D.; Vanderbilt, D. Theory of polarization of crystalline solids. Phys. Rev. B 1993, 47, 1651–1654. (72) Vanderbilt, D.; King-Smith, R. D. Electric polarization as a bulk quantity and its relation to surface charge. Phys. Rev. B 1993, 48, 4442–4455. (73) Kirchner, B.; Hutter, J. Solvent effects on electronic properties from Wannier functions in a dimethyl sulfoxide/water mixture. J. Chem. Phys. 2004, 121, 5133–5142. (74) Person, W. B.; Zerbi, G. Vibrational intensities in infrared and Raman spectroscopy; Elsevier Science Ltd, 1982; Vol. 20. (75) Ponder, J. W.; Richards, F. M. An efficient newton-like method for molecular mechanics energy minimization of large molecules. J. Comp. Chem. 1987, 8, 1016–1024. 36


Page 36 of 39

Page 37 of 39

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


(76) Kearsley, S. K. On the orthogonal transformation used for structural comparisons. Acta Crystallogr., Sect. A: Found. Crystallogr. 1989, 45, 208–210. (77) Kudin, N. K. Eckart axis conditions and the minimization of the root-mean-square deviation: Two closely related problems. J. Chem. Phys. 2005, 224105. (78) Marinica, D. C.; Gregoire, G.; Desfrancois, C.; Schermann, J. P.; Borgis, D.; Gaigeot, M. P. Ab Initio Molecular Dynamics of Protonated Dialanine and Comparison to Infrared Multiphoton Dissociation Experiments. J. Phys. Chem. A 2006, 110, 8802–8810. (79) Mathias, G.; Ivanov, S. D.; Witt, A.; Baer, M. D.; Marx, D. Infrared Spectroscopy of Fluxional Molecules from (ab Initio) Molecular Dynamics: Resolving Large-Amplitude Motion, Multiple Conformations, and Permutational Symmetries. J. Chem. Theory. Comput. 2012, 8, 224–234. (80) Frisch, M. J.; Trucks, G. W.; Cheeseman, J. R.; Scalmani, G.; Caricato, M.; Hratchian, H. P.; Li, X.; Barone, V.; Bloino, J.; Zheng, G. et al. Gaussian 09. version 1, Gaussian Inc. Wallingford CT 2009. (81) Becke, A. D. Density-functional exchange-energy approximation with correct asymptotic behavior. Phys. Rev. A 1988, 38, 3098–3100. (82) Lee, C.; Yang, W.; Parr, R. G. Development of the Colle-Salvetti correlation-energy formula into a functional of the electron density. Phys. Rev. B 1988, 37, 785–789. (83) Milani, A.; Castiglioni, C. Modeling of Molecular Charge Distribution on the Basis of Experimental Infrared Intensities and First-Principles Calculations: The Case of CH Bonds. J. Phys. Chem. A 2010, 114, 624–632. (84) Gaigeot, M. P.; Vuilleumier, R.; Sprik, M.; Borgis, D. Infrared Spectroscopy of N-

37



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

Methylacetamide Revisited by ab Initio Molecular Dynamics Simulations. J. Chem. Theory. Comput. 2005, 1, 772–789. (85) Gregoire, G.; Gaigeot, M. P.; Marinica, D. C.; Lemaire, J.; Schermann, J. P.; Desfrancois, C. Resonant infrared multiphoton dissociation spectroscopy of gas-phase protonated peptides. Experiments and CarâĂŞParrinello dynamics at 300 K. Phys. Chem. Chem. Phys. 2007, 9, 3082–3097. (86) Beck, J. P.; Cimas, A.; Lisy, J. M.; Gaigeot, M.-P. OH anharmonic vibrational motions in Cl− · · · CH3 OH ionic clusters. Combined IRPD experiments and AIMD simulations. Spectrochim. Acta, Part A 2014, 119, 12–17.

38


Page 38 of 39

Page 39 of 39



Combining Static and Dynamical Approaches for Infrared Spectra

Recommend Documents