Infrared Spectroscopy and BornâOppenheimer Molecular Dynamics

Subscriber access provided by UNIV OF CALIFORNIA SAN DIEGO LIBRARIES

Article

Infrared Spectroscopy and Born-Oppenheimer Molecular Dynamics Simulation Study on Deuterium Substitution in the Crystalline Benzoic Acid Maciej G#ug, Mateusz Zbigniew Brela, Marek Boczar, Andrzej M. Turek, #ukasz Boda, Marek Janusz Wojcik, Takahito Nakajima, and Yukihiro Ozaki J. Phys. Chem. B, Just Accepted Manuscript • DOI: 10.1021/acs.jpcb.6b10617 • Publication Date (Web): 22 Dec 2016 Downloaded from http://pubs.acs.org on December 29, 2016

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

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 48

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

Infrared Spectroscopy and Born-Oppenheimer Molecular Dynamics Simulation Study on Deuterium Substitution in the Crystalline Benzoic Acid. Maciej Gług1*, Mateusz Z. Brela1,3*, Marek Boczar1, Andrzej M. Turek1, Łukasz Boda1, Marek J. Wójcik1, Takahito Nakajima2, Yukihiro Ozaki3 1

Faculty of Chemistry, Jagiellonian University, 30-060 Kraków, Ingardena 3, Poland

2

RIKEN, Advanced Institute for Computational Science, 7-1-26, Minatojima-minami-machi,

Chuo-ku, Kobe, Hyogo 650-0047, Japan 3

Department of Chemistry, School of Science and Technology, Kwansei Gakuin University,

Sanda, Hyogo 669-1337, Japan

*To whom all correspondence should be addressed. email: M.G. [email protected], M.Z.B. [email protected]

1 Environment ACS Paragon Plus


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 In this study we present complementary computational and experimental studies of hydrogen bond interaction in crystalline benzoic acid and its deuterated and partially deuterated derivatives. The experimental part of the presented work includes preparation of partially deuterated samples and measurement of attenuated total reflection (ATR)-FTIR spectra. Analysis of the geometrical parameters and time course of dipole moment of crystalline benzoic acid and its deuterated and partially deuterated derivatives by Born-Oppenheimer molecular dynamics (BOMD) enabled us to deeply analyze the IR spectra. Presented simulations based on BOMD gave us opportunity to investigate individual motion and its contribution to the IR spectra. The band contours calculated using Fourier transform of autocorrelation function are in quantitative agreement with the experimental spectra. Characterization of single bands was carried out by ‘normal coordinate analysis’. The salient point of our study is a comparison of the spectra of the deuterated and partially-deuterated crystalline benzoic acid with that of the non-deuterated one. Furthermore, we have applied the principal component analysis for analysis of the number of components in partially deuterated systems. In this study, we reveal that the arrangements of hydrogen and deuterium atoms in partially deuterated samples are random.


Page 2 of 48

Page 3 of 48

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


I. Introduction

The crystalline structures are very often formed by the dimeric structures with intrinsic hydrogen bond (HB) interactions. Further, the HB, arguably the most important non-covalent interaction in chemistry, plays a pivotal role in a broad variety of systems and phenomena in material science and engineering, biochemistry,1-3

and

supramolecular

chemistry4-6,

among

other

fields.

Additionally, new materials mainly held by HB have recently been developed as hydrogen bonded polymers7-8 and macromolecules. Recently, a notable interest in studies of the hydrogen-bonded carboxylic acid dimers is observed.9-12 One of the important examples is the crystal structure of benzoic acid, that is the simple aromatic carboxylic acid.13-21 This system has been studied over and over again, e.g. Ernst et al. have reported the results of the measurements of the temperature dependent NMR relaxation time22. Based on the experimental data and the considered models, these authors have assumed that thermal activation of the hydrogen pair transfer process occurs mainly as a result of promotion of tunneling by heavy atom rocking. Vibrational spectroscopy is one of the most useful methods employed for studying systems forming hydrogen bonds.23-25 Infrared spectra (IR) provide valuable information related to the dynamics of hydrogen bonds. 26-29 The most conspicuous effects observed in vibrational spectra of such systems are concerned with the bands originating from the X-H modes (νs). For these bands,



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 following effects are usually observed: decrease of their frequencies, increase of their intensity, and appearance of a complex fine structure. Further, the proton transfer is a foundation of the idea of the promoting modes developed by Kamerlin et al.31-32 The fundamental transitions are formally related to tunneling between two potential minima. Such approach is widely used in studies on enzyme catalysis.33-36 The protein modes which are involved in tunneling are linked to the protein promoting vibrations.32 As a complementary to the above, an analysis of the proton potentials in benzoic acid by post-MD quantization of motion will be our future work. The studies on the isotopically substituted cyclic dimers were triggered by Durlak et al.36 The authors performed CPMD calculation for acetic acid cyclic dimers assuming several deuterium substitution rates. They compared the obtained results with the collected experimental data and results of static calculations. This work has revealed a strong superiority of the molecular dynamics simulations over the static calculations especially in terms of reproducing isotope effects. After that Latajka et al.37-39 discussed many disadvantages of the applied approach including the one unit cell limitations or a dependence of the stability of the results on the time step and the energy cutoff parameter. Theoretical methods have proved highly effective in the interpretation of the experimental spectrum.16-20,27-29 They allow to understand the origins of


Page 4 of 48

Page 5 of 48

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


particular spectroscopic effects: picks positions, band shifts and couplings. Over the last decades, several models describing the hydrogen bond systems were developed.40-52 The extension of these models are molecular dynamics (MD) simulations. Particularly common choices are ab initio versions of MD such as Born-Oppenheimer Molecular Dynamics (BOMD)53.54 and Car-Parrinello Molecular Dynamics (CPMD)55. CPMD explicitly includes the electrons in the Nuclei Hamiltonian by fictitious dynamical variables –a fictitious electron mass. Such approximation excludes minimization of the electron wavefunction at every step in the trajectory. CPMD uses fictitious dynamics to keep the electrons close to the ground state, preventing the need for a costly self-consistent iterative minimization at each time step. In order to stay in the BornOppenheimer surface, the time step should be relatively small. The BOMD does not have that kind of restriction and provides information about electronic properties in every step. It should be noted that BOMD is a start point for other methods like Path Integral Molecular Dynamics (PIMD)56.57 or hybrid methods (QM/MM)58 that are also commonly used for studying systems with hydrogen bonds. Based on the trajectory obtained from ab initio molecular dynamics, it is possible to determine the dipole moment as a function of time. Upon Fourier transformation of the dipole moment the autocorrelation function can be obtained and then, the IR transition energies and relative intensities of bands can



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

be reasonably reconstructed. Moreover, the trajectory might be used for a postMD analysis like quantization of nuclei motion that reflects the distribution of analyzed frequencies. It should be pointed out that the NVT ensemble (canonical ensemble) is the statistical ensemble that represents the possible states of a considered system at some ‘fixed’ temperature. The canonical ensemble provides information about motion at a selected temperature which is a big advantage in vibrational analysis. Our previous works on hydrogen bonds in the crystal phases of vitamin C (L-ascorbic acid)28 and aspirin29 crystal forms showed a good agreement with their IR spectra after the inclusion of the thermal motion of the nuclei. The ab initio molecular dynamics approaches used by us revealed that the key parameter of a simulation is the trajectory length that provides the statistics for the snapshot methodology. However, the comparison between the CPMD and BOMD performed by Kuo et al.59 elucidated that the fictitious electronic energy increases even for a relatively small fictitious electronic mass (µ = 100 au.), which reflects a decrease in ionic temperature. Therefore, we applied the BOMD with the NVT ensemble for analyzing the hydrogen bonding interaction. In the present work, we have carried out some complementary computational and experimental studies of hydrogen bond interaction in crystalline benzoic acid and its deuterated and partially deuterated derivatives. ATR FTIR spectra were measured for benzoic acid and its deuterated and


Page 6 of 48

Page 7 of 48

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


partially deuterated samples. Afterwards, the BOMD simulations were performed for rationalization of the experimental IR spectroscopic data (Fig. 1). The Fourier transform of the autocorrelation functions calculated by TRAVIS60, describing atom positions (power spectrum) from CP2K61 simulations allowed us to analyze individual motion and its contribution to the IR spectra. Finally, as a complementary mrthod, the principal component analysis was performed. This technics gave us a possibility to determine the number of significant factors responsible for the majority of the variance in the collected spectral data. These factors represent the mixed signals (contributions) from different forms of benzoic acid dimer (fully deuterated, partially deuterated or not deuterated). Altogether, the applied methods provide a lot of useful information on motion coupling in the hydrogen-bonded cyclic benzoic acid dimers. This article is further divided in four parts as: computational details, experimental, results and discussion, and conclusions presented in the final section.



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

II. Computational Details The structure used as a starting point in our calculations was crystalline structure of benzoic acid taken from Ref. 62. It is characterized by P 21/c symmetry space group. Cell parameters and positions of atoms are listed in Table 1, and the arrangement of benzoic acid crystal cell is shown in Fig. 1. The Born-Oppenheimer molecular dynamics (BOMD) simulations for periodic systems were performed by CP2K, open source package61. The electronic structure was analyzed in terms of density functional theory (DFT). The BLYP functional63 was combined with Goedecker pseudopotential64,65. To increase accuracy we assumed a large value of the cutoff parameter – 450 Ry and tight SCF convergence criteria (10-7). The canonical ensemble (NVT) was examined by “Canonical Sampling Through Velocity” algorithm (CSVR).66 The 300 K was chosen as a target temperature for all calculations. The analysis of partially deuterated systems implicates symmetry reduction to C1, only translational symmetry elements were retained. The time step was set to 1 fs. Total time of each simulation was 40 ps (40 000 steps). It should be pointed out that tenth of ps is sufficient for an analysis of the vibrational motions. The IR spectra were calculated by the autocorrelation function of the dipole moment. The total dipole moment was determined in terms of maximally-localized Wannier functions.67


Page 8 of 48

Page 9 of 48

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 geometry optimization in periodic boundary condition was carried out as a start point. The optimized structure was fixed as the initial configuration for BOMD simulations of all systems, listed in Table 2: hydrogenated HHHH, partially deuterated HHDD and HDHD and fully deuterated DDDD. The partial deuteration (50%) was analyzed by using two forms. The first form, HHDD represents the deuteration in one benzoic acid dimer. Such distribution was analyzed by Flakus et al. in terms of the theory ofself-organization.68,69 The second form, HDHD stems from the half deuteration of the cyclic dimer. The difference lies in the way of arrangement of the deuterium atoms. It should be emphasized that in the HDHD system, the carboxylic group in the dimer contains one hydrogen atom and one deuterium atom. The principal component analysis was performed by means of authors’ own codes operating in the MATLAB environment70 (version 7.13). The standard procedures e.g. orthogonalization scheme, were implemented from software libraries.



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

III. Experimental The seven kinds of samples of crystalline benzoic acid with various degrees [15%, 30%, 50%, 60%, 85%, 100%] of deuterium substitution were prepared. The preparation method is described in Supplementary Information. It should be noted, however, that the conditions for isotopic exchange occurred only in the carboxylic group. The ATR attachment on a Thermo Scientific Nicolet IR200 spectrometer was used for measuring the IR spectra. The IR spectra were collected at a spectral resolution of 2 cm-1 with 512 scans between 550 and 4000 cm-1.

IV Results and discussion A. Geometry optimization Table 3 summarizes geometric parameters of crystalline benzoic acid after geometry optimization. The data are compared with the X-ray crystallographic measurements.62 The main differences are observed for hydrogen positions which is typical for X-ray crystal diffraction. The explanation for this finding is limitation of X-ray crystallography in detecting of the positions of hydrogen atoms. The differences may also result from thermal expansion – the X-ray measurement was performed at 293 K whereas the calculationsin the optimization procedure do not at all take temperature into account.


Page 10 of 48

Page 11 of 48

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 data in Table 3 concern only one molecule because corresponding values (lengths and angles) for four molecules in the crystal unit cell are almost the same. Anyway, during the geometry optimization the symmetry was disabled.

B. Molecular dynamics simulations First of all, we analyzed the non-deuterated system (HHHH). Table 4 presents the mean values of the geometric parameters obtained in our molecular dynamics simulations. It should be pointed out that the values of the equivalent structural parametersof all the molecules are very similar. Furthermore, these quantities compare quite well with those derived from the geometry optimization. However, the average values of the O···O lengths are greater than those for the optimized structure. That is expected in terms of the thermal motions included in the molecular dynamics. Another difference is apparent in the case of the O-H lengths for the atoms forming hydrogen bonds. It is connected with the proton transfer process (thermal hopping). In case of the cyclic dimer two protons can change their localization in synchronous (at the same time) or asynchronous (one is delayed relative to the other) process. There are several investigations concerning the proton transfer process in benzoic acid. Both mechanisms were reported in the literature.17,18 The analysis of our



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 48

simulations leads to the conclusion that this process occurs twice within 40 ps and is definitely synchronous at least for the used time resolution – 1 fs.

C. Infrared spectra of non-deuterated system The IR spectra were calculated theoretically by the formula describing the absorption coefficient α(ω)71:  2βπω 2  ∞  ⋅ ∫ exp(−iωt ) M (0)M (t ) dt , α (ω ) =   3Vcn(ω )  − ∞

(1)

where T is temperature, V is volume of cell, β is 1/(kT), c is speed of light in vacuum, n(ω) is the refractive index of the medium and M is the total dipole moment. The brackets denote the autocorrelation function. A temporal function M(t) is the sum of the electronic and ionic dipole moments computed in every step of the simulations. Function α(ω) can be treated as approximation of the measured IR spectrum. These calculations were performed using ‘Fourier’ program72. Figure 2 compares the calculated IR spectrum of non-deuterated benzoic acid (HHHH) with its ATR FT-IR spectrum. It can be noted that not only frequencies but also intensities of the IR bands were calculated with good accordance. However, some features in the theoretical spectrum are missing – several bands existing in the experimental spectrum are absent in the theoretical one. Especially in the low frequency region where some experimental bands do


Page 13 of 48

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


not have their counterparts in the calculated spectrum. For example, three bands at 800, 1180, and 1400 cm-1 are not reproduced. This might be an intrinsic feature of our model - one unit cell limitation. In Figure S1 the system chosen for our calculations is depicted. It contains two benzoic acid cyclic dimers; the rest of the crystal environment is approximated by the periodic boundary conditions. All the moves in the crystal lattice are merely estimated by moves in our model – one unit cell model that may be responsible for some dropouts in structural motions (reflected in the low frequency region). The assignments of IR bands were performed using ‘normal coordinate analysis’ (NCA) proposed by Thomas et al.73. The calculations were done in Travis package.60 The main idea of ‘normal coordinate analysis’ is based on conversion of Cartesian coordinates to generalized normal coordinates related to a selected reference structure. The methods based on coordinate transformation are not new. However, the NCA method implemented in the Travis package60 is an extension of the transformation idea to the autocorrelation functions obtained in molecular dynamics simulations. In this work the reference structure was chosen as the optimized structure of benzoic acid dimer. The power spectra calculated for specific new coordinates usually should give narrow bands. The maximum of a band can be interpreted as a frequency of a normal mode. The results of the normal coordinate analysis are compared with the IR theoretical and experimental spectra in Fig. 3. The interpreted bands are listed 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

Table 4. Additionally, we present assignments proposed by Yaremko et. all.17 Except w,x which represent ν(O-H) modes all bands obtained by ‘Normal coordinate analysis’ are fairly narrow. In cases of w,x bands, a considerable band broadening is observed. It is obviously connected with the existence of hydrogen bonds. In the 0 - 1700 cm-1 region the calculated bands are characterized by slightly lower wavenumbers than experimental ones. On the other hand, the high frequency bands (ν(O-H), ν(C-H) ) in the calculated spectrum have higher wavenumbers than their counterparts observed in the experimental IR spectrum.

D. Infrared spectra of deuterated compounds In the previous section, we have demonstrated that combining the BOMD simulations with the normal coordinate analysis is very useful in the analysis and interpretation of experimental spectra. Therefore, we have used this approach for interpretation of the IR spectra of partially deuterated systems. Note that such systems can contain three following types of dimers: the dimer with two hydrogen atoms in hydrogen bonds (HH), that with one hydrogen atom and one deuterium atom in hydrogen bonds (HD) and that with two deuterium atoms in hydrogen bonds (DD). By a ‘ratio of substitution’ (α) we understand the fraction of deuterium atoms to total amount of hydrogen and deuterium atoms within the carboxylic groups. The important issue is how many HH, HD


Page 14 of 48

Page 15 of 48

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 DD dimers are formed for aspecific ‘ratio of substitution’ when random distribution of the deuterium atoms is assumed. Figure 4 displays the mole fractions of the considered kinds of dimers as a function of the ratio ofsubstitutionobtained for a random distribution generated by using the Monte Carlo method. As shown in the figure there are 0.25 HH dimers, 0.25 DD dimers and 0.50 HD dimers for a 50 % deuterated system. Figure 4 depicts concentrations of three kinds of dimers: HH, HD and DD in the samples under the study as a function of the substation ratio (α). The sample characterized by α close to 0 has xHH ≈ 1-α, xHD ≈ α and xDD ≈ 0 whereas that characterized by α close to 1 has xHH ≈ 0, xDD ≈ α and xHD ≈ 1 - α. In other words, the samples which are softly substituted by deuterium have only two kinds of dimers – HH and HD (DD is practically absent), while almost completely deuterated systems consist of two types of dimers – DD and HD (HH is practically absent). Another important observation is that if α = 1/3 then xHH = xHD and if α = 2/3 then xDD = xHD. Figure 5 depicts the ATR-FTIR spectra of crystalline benzoic acid with various ratios of H→D isotopic substitutions. Obviously, there are significant differences among these spectra. When the non-deuterated compound is compared with the fully deuterated one, the various forms of the spectra are evidently connected with the differences between vibrational modes of HH and DD dimers. However, in the context of the self-organization theory68.69, more



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

interesting are discrepancies originating from differences between the vibrational modes related to HH and HD or DD and HD dimers. In order to describe credibly the experimental spectra of the examined systems, we performed ‘normal coordinates analysis’ for trajectories of the systems HDHD and HHDD. In the case of the system HHDD we selected two reference structures – dimer HH and DD while for the system HDHD the reference structure was the HD dimer. We obtained normal modes of these three kinds of dimers. Based on these outcomes it is possible to assess differences in the modes involving D atoms between the HD and DD dimers or differences in the modes involving H atoms between the HD and HH dimers. In here, we will discuss more deeply few interesting spectral regions. The first interesting region is the region between 1040 and 1070 cm-1 (Fig. 5, part b). The observed bands are due to in-plane O-D bending modes ( δ(OD) ).17 From the normal coordinate analysis we obtained following frequencies of these modes: 1027 cm-1 (antisymmetric) and 1062cm-1 (symmetric) in case of the DD dimer and 1046 cm-1 in case of the HD dimer. In the IR spectra the band assigned to the symmetric mode (1062cm-1) is not visible, and thus it is expected that the band related to the DD dimer should be characterized by a lower frequency than the band originating from the fact of the presence of the HD dimer. Indeed, there are two bands in the experimental spectra – one at 1047 cm1

and another at 1065 cm-1. It is worth pointing out that the difference between


Page 16 of 48

Page 17 of 48

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 frequency difference (1065-1047=18) and calculated frequency difference (1046-1027=19) is comparable. Another part of the spectra which should be analyzed in more detail is that containing the bands related to the δ(OH) modes, see Figure 5, part c. In this case, the deuteration process causes disappearance of the O-H bending band. The normal coordinate analysis leads to following frequencies for bending modes: 1430 cm-1 (antisymmetric) and 1434cm-1 (symmetric) in the case of the HH dimer and 1434 cm-1 in the case of the HD dimer. Hence, in the IR spectra the bands due to the HH dimer should have lower frequencies than the bands associated with the HD dimer. In the experimental spectra there are bands characterized by a frequency around 1425 cm-1. The decrease in the intensities with the increase in the ratio of substitution (α) is evidently observed. For small α, according to Fig. 5, it is expected that there are more HH dimers than HD dimers. However, for large value of α, this relation becomes reversed. Hence, with the increase in α, apart from the reduction of the intensity, the band maximum should be shifted to a higher frequency. In fact, such effects are observed. This phenomenon is not so pronounced as in the case of the δ(OD) bands, which derives from significantly smaller frequency difference between the δ(OH) modes of the symmetric and non-symmetric dimers. The above discussion has shown that the partially deuterated systems contain significant amounts of the asymmetric dimers (HD). The spectral region



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 48

of the δ(OD) bands was selected for quantitative analysis. First, we created the chart which includes the ratio of xDD to xHD according to random distribution. It is presented in Fig. 6. The integral intensities (Ix) of the bands at 1065 and 1047 cm-1 can be estimated as follows I1047 = k DD xDD

(2)

I1065 = k HD x HD

(3)

where kDD and kHD are unknown constants. Further, the experimental ratio of xDD to xHD can be written as I x DD = k rel 1047 x HD I1065

(4)

where krel = kHD/kDD.

(5)

The krel characterizes each system. This parameter can be estimated as equal 0.5. It is connected with the fact that the mode δ(OD) for the DD dimer involves two deuterium atoms while the HD dimer has only one deuterium atom. Moreover, the change of the dipole moment is approximately two times greater in the case ofthe DD dimer than the HH dimer. This assumption was further validated by static quantum mechanics calculations. We performed the frequency analysis for both DD and HD dimers. The calculated frequencies and intensities of the analyzed spectra have been included in Supplementary 18 Environment ACS Paragon Plus

Page 19 of 48

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


Information, Table S2. The value of krel calculated in this way is about 0.6. In Table S3 we present the values of I1047/I1065 calculated for different substitution ratios, α. Figure 6 shows the calculated ratio of (xDD / xHD) obtained for krel = 0.5 and a random distribution curve. The points lie quite close to the random curve, suggesting that the distribution of dimers in the considered samples with different concentrations of the deuterium atoms is indeed a completely random. In that point, it should be stressed that our approach is based on the assumption that dynamics of a specific kind of dimers is independent of the kind of neighbouring dimers. To check that hypothesis we compared the power spectra of the HH dimer calculated from the simulations of the HHHH and HHDD systems with those of the DD dimer obtained from the simulations of the DDDD and HHDD systems. This procedure can be helpful in the evaluation of the influence of existing DD dimers on nearby HH dimers and vice versa during molecular dynamics. Figure 7 presents the HH dimer spectra connected with the HHHH and HHDD systems and the DD dimer spectra associated with the HHDD and DDDD systems. The analysis of the band positions leads to the conclusion that this effect is negligible – the corresponding power spectra are very much alike. Naturally, there are some discrepancies but this may partially be due to randomness of the simulations that are intrinsic to the thermostatic method. Finally, we have compared these results with those obtained from the



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

normal coordinate analysis related to dimers in specific systems and found that the difference was a few cm-1. E. Principal Component Analysis Principal Component Analysis (PCA) is frequently used for various purposes in analytical chemistry74. As a complementary method for diverse theoretical studies it offers a broad field of new applications. Recently, for instance, it has been combined with the regression methods to evaluate the correlation between the constituents of the OH stretching band and the corresponding geometry parameters of the snapshot structures of the crystalline sodium hydrogen bis(sulfate)75. In this paper, according to our best knowledge, the PCA is employed, probably for the first time, to shed a light on the spectral aspects of deuteration of benzoic acid. The first step in our analysis was preprocessing of the collected data. From the recorded IR spectra the blank spectra were subtracted (see Fig. 5) . The resulting spectra were then interpolated to provide the spectra with evenly spaced intervals (0.5 cm-1). These were collected in a matrix containing in columns, from left to right, the spectra measured for increasing degree of deuteration. The criteria of the PCA, designed by Malinowski74 to determine the proper number of significant eigenvalues retaining the majority of the data variance, reveal that the matrix is a three component system, however, when divided into two submatrices, the “low wavenumber range” matrix (< 1750


Page 20 of 48

Page 21 of 48

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


cm-1) appears to be a three component while that including the “high wavenumber range” matrix (> 1750 cm-1) shows up as a two component system (see below). The employed criteria are the most convincing in the case of the Malinowski’s indicator function (that should assume a minimum for the correct number of significant factors) and the statistical hypothesis based on one of Fisher’s tests that assigns a high significance level to the eigenvalues associated with noise. Assuming that some additional preprocessing of the data could possibly improve the Malinowski’s criteria for the overall data matrix the spectra were subjected to a pairwise shifting (i.e. the second spectrum was shifted with respect to the first, and then the third relative to that already shifted from the second spectrum and so on). In this approach the correlation optimizing shifting code devised by F. van der Berg and loosely based on V. G. van Mispelaar et. al.76 was used. The optimal shifts were typically found at 4 or 5 data points. In addition, some of the obtained spectra were also vertically shifted (the applied offsets were actually very tiny – not exceeding 0.012 absorbance units). The resulting spectra were normalized to a unit area and put, from left to right, as columns of a data matrix X , preserving as before the proper sequence of deuteration. The matrix X , sized m x n, when subjected to Malinowski’s PCA criteria appeared to be a two component spectral data matrix (see Figure 8). It can be factorized by means of Singular Value Decomposition (SVD) as follows



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 22 of 48

X = USV = USVT + E = UAT + E ,

(6)

where U (m x n) is a column orthonormal matrix that contains eigenvectors (abstract spectra) of the matrix XXT , V (n x n) is a column and row orthonormal matrix containing eigenvectors of the matrix XT X , S (n x n) is a diagonal matrix with singular values of X (squared singular values equal to eigenvalues

of

the

covariance

matrices

XXT

and

XT X

),

U (m x f ), V (n x f ), and S ( f x f ) are truncated matrices of U , V , and S ,

associated with the first f significant singular values, and E is a matrix of residual instrumental noise, typically called the error matrix. Consequently, for a two component system (f = 2) each recorded spectrum, xk , (k = 1,…,n), as well as the spectra of pure components, xA , and x B , can be represented as a linear combination of the first two abstract spectra, uα , and u β embedded in matrix U, x k = α k uα + β k uγβ ,

(7)

where α k and βk are equivalent to the entries ak,1 and ak,2 in the matrix A, (see Eq. (6)). If all elements of the vector x k in Eq. (7) are summed up then, due to the applied unit area normalization, the following expression is obtained 











1 = α k  ∑ u j ,α  + β k  ∑ u j ,β  . 

j =1

j =1



(8)

The above equation describes a so called normalization line, β k = β k (α k ) , in which all the points representing the measured three component and resolved pure component spectra should be located. To find specific values of α k and β k 22 Environment ACS Paragon Plus

Page 23 of 48

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 a two component spectrum, it suffices to make use of orthogonality of the vectors contained in the matrix U (premultiplication of both sides of Eq. (7) by uαT gives α k = uαT xk and by uTβ provides β k = u Tβ x k ). The points assigned to the

unit area normalized spectra of pure constituents define the turning points in a straight normalization line encompassing in-between all points representing the 1-norm spectra of the studied two component mixtures. This concept was originally introduced by Lawton and Sylvestre.77 The α–β normalization line obtained for the spectral set of gradually deuterated benzoic acid is shown in Figure 9. Moving along this line through the points in the direction opposite to the location of the point representing the spectrum of the completely hydrogenated species one can enter the region where a point assigned to the completely deuterated species should be located. This point called the LS point (named such after Lawton and Sylvestre) corresponds to a spectrum with the smallest possible number of negativities and thus it constitutes the best estimation of a true completely deuterated spectrum (see point D in Figure 9). The hypothetical completely deuterated spectrum and the recorded spectrum of benzoic acid prior to the deuteration process are displayed in Figure 10. Both spectra are normalized to unit area. The whole collection of unit area normalized experimental spectra and the unit area normalized LS estimated spectrum of completely deuterated benzoic acid are presented in Supplementary Information Figure S2.



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

Eventually, it is worth recalling that PCA is a useful method for analyzing the IR spectra of deuterated compounds. The number of principal components suggests the presence or absence of coupling between modes.

V Conclusions We carried out a research study on the crystalline benzoic acid including deuterium substituents. In the article we have presented the results of the simulations of the IR spectra of the crystal of benzoic acid. The IR spectra calculated by the BOMD method reproduce reasonably well the experimental spectra. The identification of bands in the calculated IR spectra was carried out by Fourier transformation of the autocorrelation function of the dipole moment. After that, the normal coordinate analysis (NCA) proposed by Thomas et. al. enabled us to characterize some selected single bands. We investigated the effects of the isotopic substitution in benzoic acid crystals on the IR spectra. We prepared partially deuterated samples and measured the ATR-FTIR spectra. The BOMD simulations allowed us to interpret the spectra of the partially deuterated systems. In addition, the analysis of the systems with various degrees of deuterium substitution gave us a possibility to discuss the coupling of the modes within those partially deuterated systems.


Page 24 of 48

Page 25 of 48

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 PCA analysis reveals two components in the IR spectra of the partially deuterated systems. Furthermore, the undertaken thorough examination of the theoretical IR spectra of different kinds of dimers unveils that the coupling between the cyclic dimers is negligible. Altogether, the PCA analysis and MD studies of the benzoic acid systems confirm the independence of the OH and O-D motions. In particular, we have proven that the arrangement of thehydrogen and deuterium atoms in the partially deuterated samples is random. To the best of our knowledge, the presented work is the first approach that combines the IR spectroscopy method, the principle component analysis and the BOMD simulations to solve a complicated problem of hydrogen-deuterium exchange in the dimeric systems of benzoic acid.

Supporting Information The Supporting Information is available free of charge on the ACS Publications website at DOI: .The experimental procedure for the deuterium substitution, geometric parameters analysis of the crystalline benzoic acid and the selected band analysis for different values of the deuterium substitution ratio.

Acknowledgements The results of the calculations presented in this paper were obtained using PL-Grid Infrastructure and resources provided by ACC Cyfronet AGH (Academic Computer Centre Cyfronet, University of Science and Technology).



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 26 of 48

This work was partially supported by the National Science Center. We would like to thank dr. Jarosław Wilamowski, Faculty of Chemistry, Jagiellonian University, Kraków, for helpful comments and suggestions experimental part of this study.


related to the

Page 27 of 48

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


References: 1.

Elias, M.; Wellner, A.; Goldin-Azulay, K.; Chabriere, E.; Vorholt, J.; Tobiasand, J.; Tawfik, D. S. The Molecular Basis of Phosphate Discrimination in Arsenate-Rich Environments Nature 2012, 491, 134–137.

2.

Olsson, M. H. M.; Siegbahn, P.E.M.; Warshel, A. Simulating Large Nuclear Quantum Mechanical Corrections in Hydrogen Atom Transfer Reactions in Metalloenzymes J Biol Inorg Chem 2004, 9, 96.

3.

Binder, W. Hydrogen Bonded Polymers 2007, Springer Berlin Heidelberg

4.

Kitano, H.; Imai, M.; Sudo, K.; Ide, M. Hydrogen-Bonded Network Structure of Water in Aqueous Solution of Sulfobetaine Polymers J. Phys. Chem. B 2002, 106, 11391-11396.

5.

Such, G. K.; Johnston A. P. R.; Caruso, F. Engineered Hydrogen-Bonded Polymer Multilayers: from Assembly to Biomedical Applications, Chem. Soc. Rev. 2011, 40, 1929.

6.

ten Brinke, G.; Ruokolainen, J.; Ikkala, O. Supramolecular Materials Based On Hydrogen-Bonded Polymers in Hydrogen Bonded Polymers edited by Binder, Wolfgang, 2007, Springer Berlin Heidelberg.

7.

Yamada, M.; Hamada, F. Supramolecular Assemblies with Water-Soluble pSulfonatothiacalix[6]arene Lanthanide Complexes: One-Dimensional ‘ladder-type’ Coordination Polymers and One-Dimensional Hydrogen-Bonded Polymers CrystEngComm 2013, 15, 5703-5712.

8.

Brocorens, P.; Linares, M.; Guyard-Duhayon, C.; Guillot, R.; Andrioletti, B.; Suhr, D.; Isare, B.; Lazzaroni, R.; Bouteiller, L. Conformational Plasticity of Hydrogen Bonded Bis-urea Supramolecular Polymers J. Phys. Chem. B 2013, 117, 5379-5386.

9.

Schuster, P. The hydrogen bond: recent developments in theory and experiments. 1. Theory. North-Holland, 1976.

10. Lundberg, H.; Tinnis, F.; Selander, N.; Adolfsson, H. Catalytic Amide Formation from Non-Activated Carboxylic Acids and Amines Chem. Soc. Rev. 2014, 43, 2714-2742. 11. Shi, H.; Wanga, Y.; Huaa, R. Acid-Catalyzed Carboxylic Acid Esterification and Ester Hydrolysis Mechanism: Acylium Ion as a Sharing Active Intermediate via a Spontaneous Trimolecular Reaction Based on Density Functional Theory Calculation and Supported by Electrospray Ionization-Mass Spectrometry Phys. Chem. Chem. Phys. 2015, 17, 30279-30291. 12. Bhattacharya, A.; Shin, J.-W.; Clawson, K. J.; Bernstein, E. R. Conformation Specific and Charge Directed Reactivity of Radical Cation Intermediates of α-Substituted (amino, hydroxy, and keto) Bioactive Carboxylic Acids Phys. Chem. Chem. Phys. 2010, 12, 9700-9712.



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

13. Marechal, Y.; Witkowski, A. Infrared Spectra of H-Bonded Systems J. Chem. Phys. 1968, 48, 3697–3705. 14. Blaise, P.; Wojcik, M. J.; Henri-Rousseau, O. Theoretical Interpretation of the Line Shape of the Gaseous Acetic Acid Cyclic Dimer J. Chem. Phys. 2005, 122, 064306– 064306. 15. Brela, M. Z.; Wójcik, M. J.; Boczar, M.; Hashim, R. Car–Parrinello Simulation of the Vibrational Spectra of Strong Hydrogen Bonds with Isotopic Substitution Effects: Application to Oxalic Acid Dehydrate Chem. Phys. Lett. 2013, 558, 88−92. 16. Boczar, M.; Boda, Ł.; Wójcik, M. J. Theoretical Modeling of the O–H Stretching IR Bands of Hydrogen-Bonded Dimers of Benzoic Acid in s0 and s1 Electronic States J. Chem. Phys. 2007, 127, 084307. 17. Yaremko, A.; Ratajczak, H.; Barnes, A.; Baran, J.; Durlak, P.; Latajka, Z. Fermi Resonance and Strong Anharmonic Effects in the Absorption Spectra of the ν-oh (ν-od) Vibration of Solid h-and d-Benzoic Acid, Chem. Phys. 2009, 364, 51–63. 18. Antoniou, D.; Schwartz, S. D. Proton Transfer in Benzoic Acid Crystals: Another Look using Quantum Operator Theory J. Chem. Phys. 1998, 109, 2287–2293 19. Flakus, H. T.; Chemecki, M. Infrared spectra of the Hydrogen Bond in Benzoic Acid Crystals: Temperature and Polarization Effects Spectrochim. Acta, Part A 2002, 58, 179– 196. 20. Boczar, M. Szczeponek, K. Wójcik, M. J. Paluszkiewicz, C. Theoretical modeling of infrared spectra of benzoic acid and its deuterated derivative J. Mol. Struct. 2004, 700, 39–48. 21. Schuster, P. The hydrogen bond: recent developments in theory and experiments. 1. Theory. North-Holland, 1976. 22. Stöckli, A.; Meier, B. H.; Kreis, R.; Meyer, R.; R. R. Ernst, R. R. Hydrogen Bond Dynamics in Isotopically Substituted Benzoic Acid Dimers J. Chem. Phys. 1990, 93, 1502. 23. Siesler, H. W.; Ozaki, Y.; Kawata, S.; Heise, H. M. Eds. Near-Infrared Spectroscopy; Wiley-VCH: Weinheim, Germany, 2002. 24. Ozaki, Y.; McClure, W. F.; Christy, A. A. Eds. Near-Infrared Spectroscopy in Food Science and Technology, Wiley-Interscience: Hoboken, NJ, 2007. 25. Beć, K. B.; Futami, Y.; Wójcik, M. J.; Nakajima, T.; Ozaki, Y. Spectroscopic and Computational Study of Acetic Acid and Its Cyclic Dimer in the Near-Infrared Region J. Phys. Chem. A 2016, 120, 6170–6183.


Page 28 of 48

Page 29 of 48

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. Pirc, G.; J. Mavri, J.; Stare, J. Program Package for Numerical Solving Time-Independent Schrödinger Equation for Vibrational Problems: Inclusion of Coordinate Dependent Reduced Masses Vib. Spectrosc. 2012, 58, 153-162. 27. Brela, M.; Stare, J.; Pirc, G.; Sollner Dolenc, M.; Boczar, M.; Wójcik, M. J.; Mavri, J. Car-Parrinello Simulation of the Vibrational Spectrum of a Medium Strong Hydrogen Bond by Two-Dimensional Quantization of the Nuclear Motion: Application to 2hydroxy-5-nitrobenzamide J. Phys. Chem. B 2012, 116, 4510- 4518. 28. Brela, M. Z.; Wójcik, M. J.; Boczar, M.; Witek, Ł.; Yasuda, M.;Ozaki, Y. Car–Parrinello Molecular Dynamics Simulations of Infrared Spectra of Crystalline Vitamin C with Analysis of Double Minimum Proton Potentials for Medium-Strong Hydrogen Bonds J. Phys. Chem. B 2015, 119, 7922–7930. 29. Brela, M. Z.; Wójcik, M. J.; Witek, Ł. J.; Boczar, M. ; Wrona, E.; Hashim, R.; Ozaki, Y. Born–Oppenheimer Molecular Dynamics Study on Proton Dynamics of Strong Hydrogen Bonds in Aspirin Crystals, with Emphasis on Differences between Two Crystal Forms J. Phys. Chem. B 2016, 120, 3854–3862. 30. Meyer, M. P.; Klinman, J. P. Modeling Temperature Dependent Kinetic Isotope Rffects for Hydrogen Transfer in a Series of Soybean Lpoxygenase Mutants: The Effect of Anharmonicity upon Transfer Distance. Chem. Phys. 2005, 319, 283-296. 31. Kamerlin, S.C.L, Mavri J., Warshel AExamining the Case for the Effect of Barrier Compression on Tunneling, Vibrationally Enhanced Catalysis, Catalytic Entropy and Related Issues FEBS Letters 2010, 584, 2759-2766. 32. Saen-Oon, S.; Ghanem, M.; Schramm, V. L.; Schwartz, S. D. Remote Mutations and Active Site Dynamics Correlate with Catalytic Properties of Purine Nucleoside Phosphorylase Biophysical Journal 2008, 94, 4078-4088. 33. Sudhir C. Sharma, S. C.; and Judith P. Klinman, J. P. Experimental Evidence for Hydrogen Tunneling when the Isotopic Arrhenius Prefactor (AH/AD) is Unity J Am. Chem. Soc. 2008, 130, 17632-3. 34. Knapp, M.J.; Rickert, K.; Klinman, J. P. Temperature-Dependent Isotope Effects in Soybean Lipoxygenase-1: Correlating Hydrogen Tunneling with Protein Dynamics J. Am. Chem. Soc. 2005, 124, 3865-3874. 35. Kržan, M.; Vianello, R.; Maršavelski, A.; Repič, M.; Zakšek, M.; Kotnik, K.; Fijan, E.; Mavri, J. The Quantum Nature of Drug-Receptor Interactions: Deuteration Changes Binding Affinities for Histamine Receptor Ligands. PLoS ONE 2016, 11, e0154002. 36. Durlak, P.; Latajka, Z. Car-Parrinello Molecular Dynamics and Density Functional Theory Simulations of Infrared Spectra for Acetic Acid Monomers and Cyclic Dimers Chem. Phys. Lett. 2009, 477, 249-254.



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. Dopieralski, P.; Latajka, Z.; Olovsson, I. Proton Transfer Dynamics in Crystalline Maleic Acid from Molecular Dynamics Calculations J. Chem. Theory Comput. 2010, 6, 14551461. 38. Dopieralski, P.; Latajka, Z.; Olovsson, I. Proton Transfer Dynamics in the (HCO3-)2 Dimer of KHCO3 from Car-Parrinello and Path Integrals Molecular Dynamics Calculations Acta Cryst. B 2010, 66, 222-228. 39. Durlak, P.; Latajka, Z.; Berski, S. A Car-Parrinello and Path Integral Molecular Dynamics Study of the Intramolecular Lithium Bond in the Lithium 2-pyridyl-N-oxide Acetate J. Chem. Phys. 2009, 131, 024308-024316. 40. Wójcik, M. J. Theory of the Infrared Spectra of the Hydrogen Bond in Molecular Crystals International Journal of Quantum Chemistry 1976, 10, 747–760. 41. Witkowski, A.;Wójcik, M. Infrared Spectra of Hydrogen Bond a General Theoretical Model, Chem. Phys. 1973, 1, 9–16. 42. Wójcik, M. J. Fermi Resonance in Dimers: a Model Study Molecular Physics 1978, 36, 1757–1767. 43. Henri-Rousseau, O. Chamma, D. IR Spectral Density of Weak H-Bonded Complexes Involving Damped Fermi Resonances. i. Quantum Theory Chem. Phys. 1998, 229, 37– 50. 44. Yaremko, A.; Ratajczak, H.; Baran, J.; Barnes, A.; Mozdor, E.; Silvi, B. Theory of Profiles of Hydrogen Bond Stretching Vibrations: Fermi–Davydov Resonances in Hydrogen-Bonded Crystals Chem. Phys. 2004, 306, 57–70. 45. Rösch, N.; Ratner, M. A. Model for the Effects of a Condensed Phase on the Infrared Spectra of Hydrogen-Bonded Systems J. Chem. Phys. 1974, 61, 3344–3351. 46. Bratos, S. Profiles of Hydrogen Stretching IR Bands of Molecules with Hydrogen bonds: a Stochastic Theory. i. Weak and Medium Strength Hydrogen Bonds J. Chem. Phys., 1975, 63, 3499–3509. 47. Robertson, G.; Yarwood, J. Vibrational Relaxation of Hydrogen-Bonded Species in Solution. i. Theory Chem. Phys. 1978, 32, 267–282. 48. Henri-Rousseau O. ; Blaise, P. The Infrared Spectral Density of Weak Hydrogen Bonds within the Linear Response Theory Advances in Chemical Physics, 1998, 103, 1–186. 49. Boczar, M.; Boda, Ł.; Wójcik, M. J. Theoretical Model for a Tetrad of Hydrogen Bonds and Its Application to Interpretation of Infrared Spectra of Salicylic Acid J. Chem. Phys. 2006, 124, 084306.


Page 30 of 48

Page 31 of 48

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


50. Boczar, M.; Boda, Ł.; Wójcik, M. J. Theoretical Model of Infrared Spectra of Hydrogen Bonds in Molecular Crystals and Its Application to Interpretation of Infrared Spectra of 1-Methylthymine J. Chem. Phys. 2006, 125, 084709. 51. Boczar, M.; Kwiendacz, J.; Wójcik, M. J. Theoretical and Spectroscopic Study of Infrared Spectra of Hydrogen-Bonded 1-Methyluracil Crystal and Its Deuterated Derivative. J. Chem. Phys. 2008, 128, 164506. 52. Wójcik, M. J.; Kwiendacz, J.; Boczar, M.; Boda, Ł.; Ozaki, Y. Theoretical and Spectroscopic Study of Hydrogen Bond Vibrations in Imidazole and Its Deuterated Derivative Chem. Phys. 2010, 372, 72. 53. Born, M.; Oppenheimer, R. Annalen der Physik 1927, 84, 457–484. 54. Marx, D.; Hutter, J. Ab Initio Molecular Dynamics, Cambridge University Press, Cambridge, 2009. 55. Car, R.; Parrinello, M. U Unified Approach for Molecular Dynamics and DensityFunctional Theory Phys. Rev. Lett. 1985, 55, 2471-2474. 56. Marx, D.; Parrinello, M. Ab Initio Path Integral Molecular Dynamics: Basic Ideas. J. Chem. Phys. 1996, 104, 4077. 57. Marx, D.; Parrinello, M Ab Initio Path Integral Molecular Dynamics. Zeitschrift für Physik B Condensed Matter 1994, 95, 143-144. 58. Warshel, A.; Levitt, M. Theoretical Studies of Enzymic Reactions: Dielectric, Electrostatic and Steric Stabilization of the Carbonium Ion in the Reaction of Lysozyme J. Mol. Biol. 1976, 103, 227–49. 59. Kuo, I-F. W.; Mundy, C. J.; McGrath, M. J.; Siepmann, J. I. Time-Dependent Properties of Liquid Water: A Comparison of Car−Parrinello and Born−Oppenheimer Molecular Dynamics Simulations. J. Chem. Theory Comput. 2006, 2, 1274–1281. 60. Brehm, M.; Kirchner, B. Travis-a Free Analyzer and Visualizer for Monte Carlo and Molecular Dynamics Trajectories J. Chem. Inf. Model. 2011, 51, 2007–2023. 61. 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.Comput. Phys. Commun.2005,167, 103−128. 62. Bruno, G.; Randaccio, L. A Refinement of the Benzoic Acid Structure at Room Temperature Acta Crystallogr., Sect. B: Struct. Crystallogr. Cryst. Chem. 1980, 36, 1711–1712. 63. 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.



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

64. Goedecker, S.; Teter, M.; Hutter, J. Separable Dual-Space Gaussian Pseudopotentials Phys. Rev. B 1996, 54, 1703. 65. Hartwigsen, C.; Gœdecker, S.; Hutter, J. Relativistic Separable Dual-Space Gaussian Pseudopotentials from H to Rn Phys. Rev. B 1998, 58, 3641. 66. Bussi, G.; Donadio, D.; Parrinello, M. Canonical Sampling Through Velocity Rescaling J. Chem. Phys. 2007, 126, 014101. 67. Marzari, N. Vanderbilt, D. Maximally Localized Generalized Wannier Functions for Composite Energy Bands Phys. Rev. B 1997, 56, 12847. 68. Flakus, H. T. Vibronic Model for H/D Isotopic Self-Organization Effects in Centrosymmetric Dimers of Hydrogen Bonds J. Mol. Struct. 2003, 646, 15–23. 69. Flakus, H. T.; Michta, A. Strong Coupling Effects in the Polarized IR Spectra of the Chain Hydrogen Bond Systems in Imidazole Crystals: H/D Isotopic Self-Organization Effects in the IR Spectra of Isotopically Diluted Imidazole Single Crystals J. Mol. Struct. 2004, 707, 17–31. 70. MATLAB version 7.13, (2011), (computer software), The MathWorks Inc., Natick, Massachusetts. 71. Ramı́rez, R.; López-Ciudad, T.; Kumar P. P.; Marx, D. Quantum Corrections to Classical Time-Correlation Functions: Hydrogen Bonding and Anharmonic Floppy Modes. J. Chem. Phys. 2004, 121, 3973-3983. 72. Forbert, H.; Kohlmeyer, A. Fourier, v 1.1, 2002-2008. 73. Thomas, M.; Brehm, M.; Hollóczki, O.; Kelemen, Z.; Nyulászi, L.; Pasinszki, T.; Kirchner, B. Simulating the Vibrational Spectra of Ionic Liquid Systems: 1-Ethyl-3Methylimidazolium Acetate and its Mixtures J. Chem. Phys. 2014, 141, 024510. 74. Malinowski, E. R. Factor Analysis in Chemistry, John Wiley & Sons, Inc., New York, Edition 3, 2002. 75. Pirc, G.; Mavri, J.; Novič, M.; Stare, J. Virtually Nonexistent Correlation between the OH Stretching Frequency and the Instantaneous Geometry in the Short Hydrogen Bond of Sodium Hydrogen Bis(sulfate): Advanced Chemometrics Analysis J. Phys. Chem. B, 2012, 116, 7221. 76. van Mispelaar, V. G.; Tas, A. C.; Smilde, A. K.; Schoenmakers, P. J.; van Asten, A. C. Quantitative Analysis of Target Components by Comprehensive Two-Dimensional Gas Chromatography J Chromatogr A. 2003, 1019, 15–29. 77. Lawton, W. H.; Sylvestre, E. A. Self Modeling Curve Resolution Technometrics, 1971, 13, 617–633.


Page 32 of 48

Page 33 of 48

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


Table 1. Atomic positions in the system displayed in Fig. 1. The space group is P 21/c. Cell parameters: a = 5.510(5) Å, b = 5.157(6) Å, c = 21.973(8) Å, α = 90o, β = 97.41(8)o, γ = 90o 62 ATOM O1 O2 C1 C2 C3 C4 C5 C6 C7 H1 H2 H3 H4 H5

FRACTIONAL COORDINATES x y z 0.22130 0.23720 0.01320 -0.0908 0.14100 0.06470 0.10010 0.27130 0.05760 0.18160 0.47400 0.10370 0.38490 0.62650 0.09700 0.45960 0.81820 0.14060 0.33170 0.85340 0.19020 0.12990 0.70240 0.19650 0.05330 0.51330 0.15340 0.47300 0.59900 0.06100 0.60600 0.92400 0.13600 0.38200 0.99100 0.22100 0.04100 0.73000 0.23200 -0.09100 0.40500 0.15800



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

Table 2. Systems under study. Prot → Deueterim exchange (column 2) is related to Fig. 1. Percent of deuteration is correlated only with the presence of atoms forming hydrogen bonds. LABEL PROT → DEUETERIUM EXCHANGE HHHH DDDD 15,45,30,60 HHDD 15,45 HDHD 15,30

DESCRIPTION 0% deuterated 100% deuterated 50% deuterated 50% deuterated*


Page 34 of 48

Page 35 of 48

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 3. Optimized geometry (CALC) of crystalline benzoic acid compared with experimental data (EXP) taken from Ref. 62. BOND C33-O31 C33-O32 O32-H45 C34-C33 C34-C35 C35-C36 C36-C37 C37-C38 C38-C39 C39-C34 C35-H40 C36-H41 C37-H42 C38-H43 C39-H44 O32···O1 H45-O1

CALC [Å] 1.258 1.325 1.040 1.492 1.421 1.409 1.416 1.413 1.405 1.425 1.099 1.099 1.101 1.100 1.098 2.594 1.557

EXP [Å] 1.264 1.275 0.897 1.484 1.392 1.401 1.384 1.379 1.387 1.390 0.990 0.990 0.995 0.983 0.987 2.633 1.736

ANGLE O31-C33-O32 C33-O32-H45 O32-C33-C34 C33-C34-C35 C34-C35-C36 C35-C36-C37 C36-C37-C38 C37-C38-C39 C38-C39-C34 C39-C34-C35 C34-C35-H40 C35-C36-H41 C36-C37-H42 C37-C38-H43 C38-C39-H44

CALC [o] 122.99 110.77 115.84 118.91 119.96 120.08 119.99 120.42 119.78 119.76 119.85 119.55 119.68 120.27 120.81

EXP [o] 123.21 118.16 118.01 120.17 119.72 119.77 120.27 120.46 119.88 119.90 119.19 119.38 120.25 119.24 120.40



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 4. Frequencies of bands obtained by performing ‘Normal coordinate analysis’ of the HHHH system trajectory (column 3) and their counterparts present in the experimental spectrum (column 2). Interpretation of bands is disclosed in Fig. 3. LABEL a b c d e f g h i j k l m n o p r s t u v a

EXPa [cm -1] 553, (547) 665 (667) 683 707 (707) 809 (804) 935 (934) 945 1028 (1028) 1075 (1072) 1127 (1128) 1175 (1179) 1186 1294 (1292) 1327 (1326) 1424 (1423) 1453 (1454) 1498 (1496) 1584 (1584) 1604 (1601) 1687 (1687)

CALC[cm -1] 531 651 675 697 778 920 926 1016 1075 1095 1156 1175 1274 1302 1394 1424 1438 1468 1541 1572 1609

data in brackets taken from Ref. 17.


Page 36 of 48

Page 37 of 48

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


Figure Captions: Figure 1. Benzoic acid cell. Geometric parameters are listed in Table 1. The structure of two dimers used for calculation is shown in Supplementary Information, Fig. S1. Figure 2. The comparison of ATR-FTIR spectrum (red,b) and theoretical IR spectrum (green,a) of crystalline benzoic acid (HHHH system). Figure 3. ATR-FTIR spectrum of crystalline benzoic acid (bottom panel), theoretical IR spectrum of crystalline benzoic acid (HHHH system, central panel) and calculated (observed in IR spectrum) bands obtained by ‘Normal coordinate analysis’ (top part). Figure 4. The concentrations of HH, DD and HD types of benzoic acid dimer as a function of ratio of substitution assuming a random distribution The mole fractions were obtained using conditional probability (Please see section IV.D). Figure 5. ATR-FTIR spectra of crystalline benzoic acid with various ratios of H→D isotopic substitutions. Figure 6. The proportion of xDD to xHD in function of substitution ratio. It is obtained based on the distribution included in Fig. 4. The squares presented the valued calculated from experimental results. Figure 7. The power spectra of DD dimer obtained from trajectories of DDDD and HHDD (part a) systems and the power spectra of HH dimer obtained from trajectory of HHHH and HHDD systems (part b). Figure 8. Malinowski’s PCA criteria applied to the preprocessed data matrix (the spectra in the range 550 – 4000 cm-1). Figure 9. Section of alpha-beta normalization line closed by two extreme points representing a hydrogenated and an LS deuterated species. Figure 10. Pure component unit area normalized absorption spectra of completely hydrogenated (black) and completely deuterated (red) benzoic acid.



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

Figure 1. Benzoic acid cell. Geometric parameters are listed in Table 1. The structure of two dimers used for calculation is shown in Supplementary Information, Fig. S1. Fig. 1. 450x194mm (96 x 96 DPI)

ACS Paragon Plus Environment

Page 38 of 48

Page 39 of 48

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


Figure 2. The comparison of ATR-FTIR spectrum (red,b) and theoretical IR spectrum (green,a) of crystalline benzoic acid (HHHH system). Fig. 2. 254x190mm (96 x 96 DPI)



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

Figure 3. ATR-FTIR spectrum of crystalline benzoic acid (bottom panel), theoretical IR spectrum of crystalline benzoic acid (HHHH system, central panel) and calculated (observed in IR spectrum) bands obtained by ‘Normal coordinate analysis’ (top part). Fig. 3. 254x190mm (96 x 96 DPI)


Page 40 of 48

Page 41 of 48

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


Figure 4. The concentrations of HH, DD and HD types of benzoic acid dimer as a function of ratio of substitution assuming a random distribution The mole fractions were obtained using conditional probability(Please see section IV.D). Fig. 4. 360x278mm (96 x 96 DPI)



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

Figure 5. ATR-FTIR spectra of crystalline benzoic acid with various ratios of H→D isotopic substitutions. Fig. 5. 790x762mm (96 x 96 DPI)


Page 42 of 48

Page 43 of 48

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


Figure 6. The proportion of xDD to xHD in function of substitution ratio. It is obtained based on the distribution included in Fig. 4. The squares presented the valued calculated from experimental results. Fig. 6. 360x278mm (96 x 96 DPI)



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

Figure 7. The power spectra of DD dimer obtained from trajectories of DDDD and HHDD (part a) systems and the power spectra of HH dimer obtained from trajectory of HHHH and HHDD systems (part b). Fig. 7. 361x269mm (96 x 96 DPI)


Page 44 of 48

Page 45 of 48

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


Figure 8. Malinowski’s PCA criteria applied to the preprocessed data matrix (the spectra in the range 550 – 4000 cm-1). Fig. 8. 191x132mm (96 x 96 DPI)



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

Figure 9. Section of alpha-beta normalization line closed by two extreme points representing a hydrogenated and an LS deuterated species. Fig. 9. 212x138mm (96 x 96 DPI)


Page 46 of 48

Page 47 of 48

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


Figure 10. Pure component unit area normalized absorption spectra of completely hydrogenated (black) and completely deuterated (red) benzoic acid. Fig. 10. 515x386mm (95 x 95 DPI)



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

Graphical abstract GA 79x39mm (96 x 96 DPI)


Page 48 of 48

Infrared Spectroscopy and BornâOppenheimer Molecular Dynamics

Recommend Documents

Infrared Spectroscopy and BornâOppenheimer Molecular Dynamics