Distinctive behavior and two-dimensional vibrational dynamics of water molecules inside glycine solvation shell

Abstract

We present a first principles molecular dynamics study of a deuterated aqueous solution of a single glycine moiety to explore the structure, dynamics, and two-dimensional infrared spectra of water molecules found in the solvation shell of glycine. Structural aspects of the solvation shell are analyzed by calculating radial, spatial, and combined distribution functions. Wavelet analysis of trajectories computes the average vibrational stretch frequency of all the OD modes at ∼2416 cm⁻¹. We observe our calculated average O–D vibrational stretch frequency in aqueous glycine, ∼2416 cm⁻¹ matches with the earlier reported values. Computational tools enable the calculation of frequencies inside the solvation shell as well as the bulk OD modes. In the normalized frequency distribution, the average frequency inside the hydrophilic solvation shell of aminic nitrogen and carboxyl oxygen is found to be 2425 and 2423 cm⁻¹, respectively. Contrary to the frequency range in hydrophilic solvation, the average frequency around hydrophobic methylene is ∼2550 cm⁻¹, and the hydrogen bond and dangling lifetimes are 0.42 and 0.64 ps, respectively. These data confirm the presence of more non-hydrogen bonded (free OD) or “dangling” water molecules in the methylene solvation shell. Our results also show that the carboxylate–water hydrogen bonds are stronger than water molecules interacting with protons of the amine group. Two-dimensional infrared (2D-IR) spectroscopy probes femtosecond dynamics in chemical processes. The 2D-IR peak shapes in the respective ammonium and carboxylate hydrophilic solvation shells are found to be almost similar. However, a closer view reveals a slightly broader component of the diagonal line width for the OD stretch in the carboxylate hydration shell, which exhibits the differences in local environments surrounding the hydrophilic species. Our calculations of the frequency distribution, two-dimensional infrared (2D-IR) spectrum, the decay of frequency fluctuations, the lifetime of hydrogen bonds, and dangling modes unveil the distinct behavior of water molecules in solvation regions.

---

1. Introduction

In biological systems, molecular-level understanding of solute–solvent interactions is considered to be impossible without explicit aqueous solvent molecules.^1–4 The basic building blocks of proteins, the amino acids are organic compounds with the general formulae NH₂–CH(R)–COOH. When the alkyl substituent R is replaced by hydrogen (R = H), the simplest amino acid is glycine. In neutral solution (pH = 7), glycine exists in zwitterionic form,^5,6 which is amphoteric in nature,⁷ and acts as an acid as well as a base. Calculations on the stability of the zwitterionic form of glycine indicated the process of intramolecular proton transfer from the carboxylic oxygen to the amine nitrogen.^8,9 Computational results suggested water-mediated proton transfer.^8,10,11 One water molecule was not enough to stabilize glycine;¹² at least two water molecules¹³ were required to stabilize the zwitterionic conformer. However, a minimum of four water molecules was required for stabilizing the zwitterionic state of glycine.¹⁴ Recently, eight or nine water molecules were theoretically predicted to stabilize glycine solvation.^15–17 Like many biomolecules such as proteins contain non-polar hydrophobic and polar hydrophilic regions,¹⁸ glycine provides a variety of chemical environments around itself because of its hydrophilic and hydrophobic functional groups. Small organic compounds or osmolytes such as the N-methylated glycine derivatives are produced in the living cells under unfavorable environmental conditions.^19–21 However, the mechanism of stabilization of the globular proteins by osmolytes^22–24 is still debatable²⁵ within the experimental and modeling approaches.^26–30 One hypothesis^31–34 demonstrates stabilization in terms of preferential exclusion of osmolytes from the protein surface, while the other considers solution structure and solvent properties to understand protein stability.³⁵ Glycine serves as a protein precursor,³⁶ buffering agent³⁷ in antacids, analgesics, and antiperspirants. It is used as a metal complexing agent,³⁸ and inhibits postsynaptic potential acting as a neurotransmitter in the central nervous system.³⁹ Glycine exhibits nonlinear optical properties⁴⁰ and is widely applicable in the pharmaceutical industry.⁴¹

Several studies addressed the solvation of glycine in water^17,42 and its relevance in the field of biochemistry to investigate the process of protein–water interactions, and protein folding.¹⁸ To have a complete understanding of the behavior of glycine in aqueous solution, there is a need to throw light into the dynamical aspect together with the structural features. The structural arrangement of water molecules in aqueous glycine and its N-methyl derivatives were described by experimental IR,^14,43–45 Raman,⁴⁶ dielectric relaxation spectroscopy,⁴⁷ neutron diffraction,⁴⁸ and calorimetry⁴⁹ as well as the theoretical studies like ab initio simulations using different basis sets.^10,13 Furthermore, Monte Carlo⁵⁰ and molecular dynamics (MD)⁵¹ studies served crucial for providing the evidence for the structural changes brought about by the solvation and hydrogen-bonding networks. Hydrogen bonds play a vital role in explaining the interactions in organic–water molecules, DFT PBE0 level first principle theory^52,53 on the electronic and vibrational properties of hydrogen bonds in glycine water revealed covalent-like properties of H bonds but there wasn't any report on the timescales of the breaking and reforming dynamics of these hydrogen bonds.⁵⁴ In 1990, Basch and Stevens first published the ab initio calculations on the different possible structures of glycine·H₂O.⁵⁵

However, the heterogeneity in the solvation nature of glycine can be divided into two regions. Broadly the hydrophilic parts in zwitterionic glycine molecule consist of the solvation shell around aminic nitrogen, and carboxylic oxygen and the region around the methylene group represent the hydrophobic zone. This unique diversity withdraws interest to explore the interactions of glycine with its local environment. The vibrational frequency shifts induced by solvent water molecules around glycine were estimated by the vibrational frequency analysis (VFA) with the free energy Hessian (FE Hessian).⁵⁶ Experimentally, it was a challenging task to probe the water dynamics in the respective solvation shells in the vicinity of the solute.⁵⁷ Spectroscopic techniques such as dielectric relaxation spectroscopy,^58,59 X-ray diffraction, and neutron scattering displays translational, rotational dynamics of water molecules, and provide information about the structural patterns.^59,60 Environmental alterations cause the vibrational frequencies of the vibrational chromophore to change, and the spectral response to its surroundings need to be monitored in the most fundamental amino acid for explaining the complex chemical and biological process involving aqueous glycine.

2. Methodology

First principles molecular dynamics (FPMD) study was performed using CPMD simulation package⁶¹ employing the Car–Parrinello method.⁶² A deuterated system containing a single glycine and 30 water molecules was simulated in a cubic box, having a box length of 10.13 Å. Periodic boundary conditions were applied in three dimensions to eliminate the surface effects, and the electronic structure of this system was represented by the self-consistent Kohn–Sham equations⁶³ of density functional theory^64,65 within a plane wave basis. The plane wave cut off of 80 Ry was used, and we used Troullier–Martins atomic pseudopotentials⁶⁶ to treat core electrons. BLYP^67,68 functional with dispersion correction⁶⁹ (BLYP-D2) was used. Previous reports on water simulations exhibited a dramatic rise in density with the addition of dispersion corrections to the density functionals, but the choice of BLYP-D2 density functional reported a density value close to the liquid density of water.⁷⁰ The calculation of oxygen–oxygen radial distribution function (RDF) using the same functional reflected perfect tetrahedral hydrogen bond associations.⁷¹ The formamide–water system simulated using BLYP-D2 density functional adopting the same protocol provided an good agreement between the theoretical calculations and the experimental results.⁷² Hence, we were motivated to account for the dispersion forces in our simulations with BLYP-D2 density functional. A time step of 5 a.u. was used for the simulations with a fictitious mass of 700 a.u. to describe the electronic degrees of freedom at 300 K. We carried out wavelet analysis^73,74 of trajectories obtained from the first principles molecular dynamics simulations to generate the time-dependent instantaneous frequencies of the probe OD groups. The probe encounters different environments, and owing to varied interactions, we see the change in OD vibrational stretch frequency with time, called the phenomena of spectral diffusion. Fourier analysis generates frequencies from time-dependent signals. However, the frequency representations do not establish time-frequency relationships or impart spectral information relative to time. The time-frequency dependence introduced the concept of instantaneous frequencies represented by the following functional form f(t) = a(t)e^iØt and the derivative of the phase function Ø defines frequency. Time-series wavelet method^73,74 determines the instantaneous frequency at each time step incorporating a time localization parameter. Wavelet analysis is an appropriate method of calculating the frequencies fluctuating with time and it provides a better time localization and spectral information. A complex function of mother wavelet f(t), is expressed as a combination of time-dependent real and imaginary terms

f(t) = r(t) + ip(t) (1)

where r(t) is the real part representing instantaneous fluctuations of the OD bond, and ip(t) representing momentum in the above equation. The complex equation is a function of two variables a and b, which expresses translations and dilations in terms of the mother wavelet:ψ_a,b(t) has compact support or decays to 0 rapidly for at t → ±∞. The wavelet transform of f(t), gives the coefficients a and b. The coefficients of the expansion are given by the wavelet transform of the function:

In the above equation, a > 0, b is a real term and [small psi, Greek, macron] represents the Fourier transform of ψ. Wavelet transform considers a time-frequency representation f; the parameter b shifts the wavelets so as to store the local information of f in L_ψf(a,b) at time t = b. L_ψf(a,b) depicts the frequency of f over the time interval t = b. Wavelet theorem expresses itself in terms of two functions frequency variable a (also known as scale factor) and time-variable b. The frequency variations are regulated by the dilation parameter a. The wavelet ψ_a,b shrinks for |a| ≪ 1, mostly in the higher frequency range and for |a| ≫ 1, ψ_a,b dilates with the frequency content in the low frequency range. After wavelet analysis at each time window t = b, we get the frequency of the OD mode from the inverse of scale factor a. The scale value a governs the width of the time window. At high frequency, the small value of a automatically narrows the time window, and the time window widens at a low frequency due to the large a value. Morlet–Grossman wavelet⁷⁵ is used for analysis

Expressing eqn (4) in terms of

Substituting in eqn (2),

ψ must be real and should satisfy the admissible condition, , and λ = 1 and σ = 2 are the values for constants. These parameters are adjusted to get better resolution in the calculated frequencies. The snapshots of trajectory frames undergo Fourier transformation⁷⁶ of time-domain data to generate frequencies. Wavelet code runs through the entire trajectory and computes the frequency^77,78 of all OD modes.

The linear spectroscopy gives a time-averaged broad vibrational absorption spectrum owing to the variety of surrounding local environments. Two-dimensional infrared (2D-IR) vibrational echo spectroscopy^79,80 is a useful tool to track the molecular interactions occurring at sub-picosecond to picosecond timescales in the fastest chemical exchange reactions, solute/solvent interactions and solution dynamics,^81–83 water dynamics,^84,85 and hydrogen network evolution.⁸⁶ For accessing the dynamics in the aqueous solvation shells, a vibrational chromophore known as the probe entity is chosen. Here, OD vibrational stretch is the probe of interest. Probes respond to the changes in the local environment by interacting with them, resulting in its vibrational stretch frequency changes referred to as the phenomena of spectral diffusion.⁸⁷ Spectroscopic techniques like IR, NMR can track the chemical processes if timescales are slow. However, ultrafast vibrational spectroscopy can interrogate the molecular structure and ultrafast intermolecular motions. The time-dependent frequencies of the OD modes obtained from the time-series wavelet analysis^73,74 were used for the calculation of non-linear vibrational spectral features of the water molecules in the solvation shell of ammonium and carboxylate groups of glycine presented in Fig. 5. We employed second-order cumulant approximation of the response function^88–90 in our calculation of the 2D-IR spectrum. In this approximation, it is assumed that the fluctuating dynamics describe a Gaussian process and is valid for observing the peak shift. However, in strongly hydrogen-bonded systems like liquid water, non-Condon line shape is more realistic as it takes into account non-Gaussian frequency fluctuations and environmental impacts on the transition dipole (non-Condon effects).⁸⁸

Theoretical calculations of the 2D-IR spectroscopic observables involve response function formulation.^88,90 Three ultrashort IR pulses with wave vectors k₁, k₂, and k₃ interacts with three input electric fields E₁, E₂, and E₃ to produce an echo signal E_sig. Previously reported literature has the details of the 2D-IR experimental setup^79,85 and derivations of the mathematical formulae.^88–90 Ignoring, the excited state and rotational lifetime effects, the equations for response functions in second-order cumulant approximation is given by

R₁(t₃,t₂,t₁) = R₂(t₃,t₂,t₁) = μ₁₀⁴ exp(−i〈ω₁₀〉t₃ + i〈ω₁₀〉t₁) × exp(−g₁(t₁) + g₁(t₂) − g₁(t₃) − g₁(t₂ + t₁) − g₁(t₂ + t₃) + g₁(t₁ + t₂ + t₃)), (7)

R₃(t₃,t₂,t₁) = −μ₁₀⁴κ² exp(iΔt₃) × exp(−i〈ω₂₁〉t₃ + i〈ω₁₀〉t₁) × exp(−g₁(t₁) + g₂(t₂) − g₃(t₃) − g₂(t₂ + t₁) − g₂(t₂ + t₃) + g₂(t₁ + t₂ + t₃)), (8)

R₄(t₃,t₂,t₁) = R₅(t₃,t₂,t₁) = μ₁₀⁴ exp(−i〈ω₁₀〉t₃ − i〈ω₁₀〉t₁) × exp(−g₁(t₁) − g₁(t₂) − g₁(t₃) + g₁(t₂ + t₁) + g₁(t₂ + t₃) − g₁(t₁ + t₂ + t₃)), (9)

R₆(t₃,t₂,t₁) = −μ₁₀⁴κ² exp(iΔt₃) × exp(−i〈ω₂₁〉t₃ − i〈ω₁₀〉t₁) × exp(−g₁(t₁) − g₂(t₂) − g₃(t₃) + g₂(t₂ + t₁ ) + g₂(t₂ + t₃) − g₂(t₁ + t₂ + t₃)). (10)

The time delay between the first and second pulses is referred to as t₁, between second and third pulses as t₂, t₂ = T_w(waiting time), and the time after the third pulse as t₃. g₁(t), g₂(t) and g₃(t) are the time correlations of the fluctuating transition frequencies given by the following expressions

where Δ = 〈ω₂₁〉 − 〈ω₁₀〉 and κ is equal to μ₂₁/μ₁₀. μ₁₀ and μ₂₁ represent the transition dipole moment corresponding to 0–1 and 1–2 transitions. ω₁₀ and ω₂₁ correspond to the vibrational frequencies of the mode in the ground and excited state. Δ value gives the measure of anharmonicity. 2D IR spectrum is obtained by Fourier transformation of the response functions along the two frequency axes in terms of waiting time (t₂ = T_w) expressed as follows^88–90where,for i = 1 to 3 andfor i = 4 to 6, respectively.

At initial T_w, the final frequencies are found to be similar to the initial frequencies. The frequency correlation initiates the phenomena of spectral diffusion, and the 2D-IR line shape is seen elongated along the diagonal axis. At higher T_w, the final frequencies are different from the initial value, and the frequencies decorrelate with time. Thus, with time more spectral diffusion occurs, resulting in frequency decorrelation quantified by a symmetrical 2D-IR line shape.⁹¹ Fig. 3 depicts the diagrammatic scheme of frequency calculations based on wavelet theory and the mathematical formulation of the 2D IR spectrum from the calculated time-dependent fluctuating frequencies.

3. Results and discussion

3.1. Structural aspects of the solvation regime in aqueous glycine

To understand the structure in aqueous glycine solution, we calculated the pair correlation function or radial distribution function (RDF). The RDF, g(r), determines the atomic density at a certain distance away from the reference molecule in a system of molecules^92,93 given bywhere r_i and r_j depict the magnitude of the position vectors of the i^th and j^th particle. RDF was calculated by constructing histograms,⁹³ inserting atom pair distances in each trajectory timeframe, and it provides a detailed description of the condensed phase liquid structure.

In the RDF graph shown in Fig. 1, the location of the first peak minimum is taken as the distance cut-off for defining the solvation boundary of the corresponding intermolecular atomic pair in the glycine solution. We computed intermolecular RDFs between the atoms of glycine molecule and the oxygen (O_W) or deuterium (D_W) atoms of the neighboring water molecules, i.e., the radial distribution of N–O_W, D_N–O_W, D_C–O_W, O_C–O_W, and O_C–D_W pairs to create a solvation shell demarcation⁹⁴ from the bulk around amine nitrogen, methylene groups and carboxylic oxygen in aqueous glycine. D_N, D_C, and O_C represent the deuterium atoms attached to nitrogen, deuterium atoms of the methylene group, and two oxygen atoms adjacent to the carboxylic carbon, respectively. The corresponding number integrals (NIs) are also presented together with the RDFs in Fig. 1 by the dashed lines. The number integral NI(r) is expressed by the equation

Fig. 1 In the glycine–water system, the solid lines blue, maroon, green, magenta, and red represent the radial distribution functions (RDFs) for N–O_W, D_N–O_W, D_C–O_W, O_C–D_W, and O_C–O_W pairs, respectively and the corresponding number integrals (NI) are shown with dashed lines. The SDF of oxygen and deuterium atoms are shown with magenta and yellow colored isosurfaces (isovalue = 78), respectively, in the upper part (inset figure).

The value of the number integral at the first peak minimum of RDF estimates the coordination number or hydrogen bond number depending on the types of atoms involved. Analysis of the RDF provides the structural details in solvated zwitterionic glycine. The first peak maximum, minimum, and corresponding values of the number integral of the N–O_W RDF pair are 2.81, 3.48 Å, and 3.5 in accordance with the N–O_W RDF data given by the ab initio molecular dynamics (AIMD) study of one glycine and 30 water molecules.^8,95 Neutron diffraction experiments⁹⁶ and Car–Parrinello⁸ method reported 3 to 4 solvation water molecules around the ammonium group in glycine. The latter AIMD simulation⁸ slightly overestimates the solvation number to four, but the number is reduced to three when the H_N–O_W interactions are eliminated. A slight hump in noticed in g(N–O_W) at the minimum due to the interactions from the extra solvent water (NI = 3.5) at its vicinity. The first peak minimum of D_N–O_W RDF at 2.44 Å serves as the hydrogen bond cut-off criteria inside the aminic nitrogen solvation shell in the glycine system. Three water molecules, found in the first solvation sphere around the ammonium group, match with the previous CPMD study of glycine intramolecular proton transfer in water.⁸ This seems to be quiet interesting as we can say that each hydrogen bond donor site in the amine group is hydrogen-bonded to one of its nearest neighboring oxygen atom of water displaying strong hydrogen bond interactions. A second peak is seen in the N–O_W radial distribution function plot extending to a minimum at 4.3 Å, and the second solvation shell is not well defined because of its shallowness compared to the first peak minimum. Thus, the second hydration layer around aminic nitrogen lies undefined as in preceding reports.⁹⁵ The RDF of D_C–O_W pair around methylene functionality reflects very weak interaction with the nearby water as there is no sharp peak maximum. The solvation boundary is not well-defined with minimum extending up to 4.3 Å due to the water-repelling nature of CD₂ moiety that inhibits the approach of water molecules. Moving the attention towards the carboxylates, we calculated O_C–O_W RDF, and the first peak maximum appears at 2.85 Å matching with the previously reported values,^97,98 and the first solvation shell is at a distance of 3.30 Å, (at the first peak minimum) around carboxylic oxygen of glycine. To count the number of water molecules hydrogen-bonded to O_C in the COO solvation, we calculated the pair correlation function of O_C–D_W, the first peak minimum is observed at 2.52 Å. Within this distance, the water hydrogens are hydrogen-bonded to the carboxylic oxygen. The number integral at the first minimum is found to be 4.3 yielding one water molecule with two deuterium atoms near each O_C. On average, a total of 4.7 number of hydrogen bonds were reported within water and carboxylate oxygen atoms.⁸ However, both the values were found to be smaller than the hydration numbers 5.3 and 7 predicted by simulations using empirical force field⁹⁹ and OPLS force fields¹⁰⁰ respectively. Classical molecular dynamics study¹⁰¹ with AMBER ff03 force field showed the first peaks in RDFs of O_C–O_W, H_N–O_W, and N–O_W pairs at 2.71, 2.05, and 2.95 Å, respectively, were in good agreement with various organic crystals.¹⁰² The hydration numbers reported by the study of glycine/water mixtures¹⁰¹ in crystal/solution environments around amino and carboxylic groups were 2.97 and 5.38, respectively. Neutron diffraction measurements⁴⁸ with isotopic substitution showed 3 water molecules around the amino group on an average and 4 to 5.5 water molecules representing more hydrophilicity around the carboxyl group.

We visualized the probability of finding the oxygen and deuterium atoms around the glycine molecule in a three-dimensional space by calculating the spatial distribution function (SDF). An isosurface plot shown at the inset of Fig. 1, having isovalue of 78, was used to represent the oxygen and deuterium atoms of water molecules around glycine. The concept of isosurface¹⁰³ comes from contour lines passing through the areas with the same probability in topographic maps. SDF data of water molecules around ND₃, CD₂, and COO groups in aqueous zwitterionic glycine reveal that the atomic oxygen distribution is denser than that of hydrogen atoms around all three sites. The 3D figure representing the spatial distribution of water molecules around two hydrophilic zones, the amine nitrogen, and carboxylic oxygen enables us to determine water's structural arrangement around hydrophilic groups and distinguish it from the distribution pattern of the water molecules around hydrophobic methylene. In the SDF figure, we find a small fraction of water molecules around the apolar methylene group; this justifies the water-repelling nature of hydrophobes. However, there is a dense distribution of water molecules around water-loving amine nitrogen and carboxylic oxygen sites in glycine. A close view of the atom cloud density at the COO group in solvated glycine reveals the non-uniform distribution of water molecules around the two carboxylic oxygen atoms (O_C), which behave little differently. Moreover, splitting of the hydrated HDO infrared band⁷⁶ and ab initio molecular dynamics simulations⁹⁵ revealed carboxylate anion underwent nonequivalent interaction with hydrating water via two oxygen atoms.

In the aqueous solution, we find glycine in the zwitterion state possesses two hydrophilic hydrogen bonding sites around aminic nitrogen and carboxylic oxygen. In order to gain deep insight into the hydrogen bonding patterns, we calculated the combined distribution function (CDF), as shown in Fig. 2. The figure presents the combined plot of distance [R in pm] along the X-axis versus the angular distribution of the hydrogen bond angle formed by two vectors represented by angle [θ in degrees] along the Y-axis. The intercalated models in the top and bottom panel of Fig. 2 provide the information regarding the distance and H-bond angle considered here. SDF and CDF were computed using the TRAVIS software package.¹⁰³ In the contour plot shown in the upper panel, we find the CDF in aminic nitrogen solvation shell, considering N–O_W distance and the angle formed by the bond vectors N–D_N and D_N–O_W. The most probable angle is found to be around 15° at a distance of 300 pm, and the dark-colored angular distribution ranges from 0° to 23° within the N–O_W distance of 280 to 320 pm. Light green colored contour is visible along 90° to 135° in the same N–O_W distance range, but the probability of finding it is very less. In Fig. 2(b), the hydrophilic carboxylic group in glycine, the gnuplot exhibits CDF of O_C–O_W distance, and the angle formed by the hydrogen bond between D_W of water and acceptor site O_C, the O_W–D_W⋯O_C bond angle. In this case, the darker color depicts the most probable O_C–O_W distance at 300 pm, and the angular distribution lies along 165°. Along distance from 280 to 320 pm, the angular distribution stretches, and the O_C–O_W distance 320 pm marks the first minimum (3.30 Å) of the O_C–O_W RDF peak. We also find a light green patch along 45° to 75°, but it is a region that is less probable and not so well defined.

Fig. 2 In the top panel (a) the x-axis represents N–O_W distance (in pm), and the y-axis represents the hydrogen bond cut-off angle between two bond vectors N–D_N and D_N–O_W. The bottom panel (b) represents the CDF of O_C–O_W distance and the hydrogen bond angle that originates from O_W–D_W and O_C–D_W bond vectors.

3.2. Distribution of O–D stretch frequency inside hydrophilic and hydrophobic hydration shell of glycine

The analysis of RDF, SDF, and CDF depicts the structural arrangement and behavior of the different types of water molecules around glycine. The distribution of water molecules in glycine solvation seems to be heterogeneous, but this topic has not been explored in great detail so far. Therefore, we classify the solvation shell of aqueous zwitterionic glycine into three sections to understand the non-uniform solvation of water molecules. We categorize these as two hydrophilic regions, one around amine nitrogen and the other carboxylic oxygen, and hydrophobic methylene based on the hydrogen bonding motifs. The RDF data already shows the presence of strong hydrogen bonding environments around water-loving –ND₃ (D_N⋯O_W) and –COO (O_C⋯D_W) groups and the weakly interacting water molecules around hydrophobic methylene.

A division is made in between the solvation shell and bulk¹⁰⁴ based on RDF calculations, and the solvation boundary was defined around aminic nitrogen taking N–O_W RDF peak minimum at 3.48 Å, and the criteria for hydrogen-bonded water molecules in this nitrogen solvation was taken from D_N–O_W minimum at 2.44 Å. In a similar way, peak minima of radial distribution, O_C–O_W at 3.30 Å, and O_C–D_W at 2.52 Å were taken as the solvation shell condition and cut-off distance⁹⁴ within which water acts as the proton donor to participate in hydrogen bonding with the C-terminal oxygen in –COO. A sharp peak is absent in D_C–O_W RDF, and the solvation shell around methylene is not pronounced, but still, the distance 4.3 Å marks the limit of the solvation zone of water molecules around –CD₂. Bulk water refers to those solvent molecules that are outside the amine, carboxylic, and methylene hydration sphere.

The time-dependent vibrational stretch frequencies of all the OD modes were calculated using the wavelet method of Arevalo and Wiggins⁷⁴ and was found to be 2416 cm⁻¹. Raman spectrum recorded a bifurcated D₂O stretching region with peaks at 2190 and 2390 cm⁻¹ and revealed the former band had more intensity than the later with the increase in the concentration of dissolved molecules.¹⁰⁵ IR spectra of deuterated glycine derivatives in low-temperature argon matrix exhibited spectral features in three most stable glycine conformers for the first time and assigned the frequencies to the stretching and bending vibrations.¹⁰⁶ These matrix-isolated IR spectra of the deuterated glycine portrayed a spectral signature at 2399 cm⁻¹ for the OD stretch, and the peak characteristics were dependent on the conformation of glycine. DFT calculations, carried out using the B3LYP density functional and aug-cc-pVDZ basis set, reported the O–D stretch band at 2411 cm⁻¹ and helped to identify the IR spectral characteristics of the glycine conformers in the experiment. Our calculated value of the average frequency of all the OD bonds in deuterated glycine ∼2416 cm⁻¹ is close to the computed OD stretch frequency 2411 cm⁻¹. Solute-affected HDO spectra and DFT calculations using the polarizable continuum model showed that splitting of HDO band was due to the asymmetric distribution of water molecules around ammonium and carboxylate groups of glycine.¹⁰⁷ It is difficult to determine the frequencies inside each of the hydration regions and outside it by applying the solvation shell conditions experimentally,¹⁰⁸ but the computational tools enable solvation zone wise calculation of the normalized frequency distribution, which is displayed in Fig. 4. Inside the N–O_W solvation shell, the vibrational stretch frequency is 2425 cm⁻¹. Inside aminic nitrogen solvation shell, OD modes that form D_N⋯O_W hydrogen bond has an average frequency of 2418 cm⁻¹, exhibiting a bathochromic shift of 7 cm⁻¹. We also observe a higher value of the average frequency 2550 cm⁻¹ around methylene of the amino acid; therefore, we predict the dominant contribution of non-hydrogen bonded water molecules or the free OD modes due to weak interaction of the methylene hydrogen and nearest neighboring oxygen leading to hydrophobic hydration in glycine. A characteristic sharp peak at ∼2580 cm⁻¹ in the vibrational power spectrum marked the presence of dangling OD modes together with the hydrogen-bonded OD modes in water–carbon tetrachloride (liquid–liquid) interface.¹⁰⁹ Stretch frequency distribution in aqueous ionic solutions like carbonate,¹¹⁰ and phosphate ion¹¹¹ showed a broad spectral feature lying in the range 2500–2700 cm⁻¹ representing dangling OD groups. Inside the solvation regime of carboxylic oxygen, the average frequency is 2423 cm⁻¹, and the solvation shell O–D stretch frequency of aminic nitrogen is 2425 cm⁻¹. Similarly, in the aqueous solution of N-methyl acetamide, the average stretch frequency of the OD modes using BLYP and BLYP-D3 density functionals were reported to be 2459 and 2465 cm⁻¹, respectively, inside amine nitrogen solvation shell; moreover, 2341 and 2346 cm⁻¹, respectively, inside carbonyl oxygen hydration shell.¹¹² Till now, the observed frequency shifts are within ±7 cm⁻¹ for hydrophilic solvation and ∼25 cm⁻¹ around hydrophobic –CD₂ due to the existence of dangling OD modes of water near it. Carboxylate–water hydrogen bonds with average frequency 2397 cm⁻¹ are found to be stronger than amine–water hydrogen bonds having a value of 2418 cm⁻¹. The previous study showed that stronger hydrogen bonds were formed by the water molecules, which interact with the carboxylate group than with the hydrogens of the amine group of glycine and its N-methylated derivatives.⁷⁶ Bulk OD modes with frequency 2410 cm⁻¹ are in agreement with the computed average O–D stretch frequency value (∼2380 cm⁻¹) of heavy water,⁷⁸ and the experiments reported vibrational stretch of the OD band centered at ∼2500 cm⁻¹.¹¹³ The discrepancies are a result of the difference in electronic structure calculations with different density functionals, electronic fictitious mass, and finite basis set cut-offs.^114,115 Comparing the bulk average O–D stretch frequency (∼2410 cm⁻¹) with the average frequency of OD modes hydrogen-bonded to carboxylate (∼2397 cm⁻¹), and ammonium (∼2418 cm⁻¹); we observe that the OD groups hydrogen-bonded to carboxylate are red-shifted with respect to the water molecules hydrogen-bonded to other water at bulk, and the hydrogen-bonded OD in amine solvation is slightly blue-shifted relative to bulk. First principles molecular dynamics simulations of aqueous formamide solvation revealed carbonyl–water hydrogen bonds are stronger than hydrogen-bonded water molecules at bulk, and amine–water hydrogen bonds are weaker than bulk based on the frequency shifts.¹¹⁶ The trend observed in the glycine solution indicates the correlation between frequency and hydrogen bond strength matching the formamide–water system. In aqueous methylamine, a blue shift of 30 cm⁻¹ was reported ongoing from bulk towards amine hydration shell.¹⁰⁴ The OD modes within the solvation shell cut-off in amine hydration (2425 cm⁻¹) show a similar 15 cm⁻¹ frequency shift corresponding to the bulk OD stretch (2410 cm⁻¹) in glycine. In glycine–water system, bulk water molecules exhibit a red shift in the average vibrational stretch frequency as compared to the carboxylate hydration shell and this shift is comparable to the 34 cm⁻¹ shift observed ongoing from the carbonyl oxygen solvation towards bulk in aqueous acetone.¹¹⁷

Fig. 3 The schematic represents the method of computing the 2D IR spectrum from the time-dependent frequencies obtained from the wavelet-based frequency calculations.

Fig. 4 Normalized distribution of O–D stretch frequencies in the solvation shell of amine nitrogen is shown by the blue line. The maroon line represents the stretch frequencies of hydrogen-bonded OD modes inside the nitrogen hydration shell. Non-hydrogen bonded OD modes, i.e., the distribution of modes in the dangling state around hydrophobic methylene, are shown by the green line. Inside the solvation shell of carboxylic oxygen in glycine, the distribution of all OD modes and hydrogen-bonded OD modes are represented by red and magenta lines, respectively. The black line represents the bulk OD modes.

3.3. Ultrafast non-linear vibrational spectroscopy

The spectral changes observed in the computed 2D-IR spectra around ammonium and carboxylate groups at different waiting times (T_w) 0.01, 0.1, 1, and 2 ps shown in Fig. 5, quantify vibrational spectral diffusion of OD modes of water molecules occurring at each timescale. The 2D-IR peak shapes are found to be similar at each T_w due to dipolar interactions between positively charged ND₃⁺ and negatively charged COO⁻ in aqueous zwitterionic glycine. The peak on the ω₁–ω₃ diagonal originates from stimulated emission and ground-state bleach from ν = 0 → 1 energy levels, and the off-diagonal peak represents vibrational excited state absorption for ν = 1 → 2 transition, shifted by anharmonicity to a lower wavenumber. In the 2D-IR contour plot around the carboxylate group, the diagonal line width of the lobes is slightly broad as compared to the ammonium; this trend matches with the full width at half-maximum (fwhm) of the O–D stretch frequency distribution inside the solvation shells in glycine.

Fig. 5 2D IR spectroscopy of the OD stretch modes at four different time delays (T_w) to monitor the phenomena of spectral diffusion. Left panel (a) represents the 2D IR spectra of the solvation shell water molecules around the ammonium group, and the spectral features in the carboxylate hydration shell are shown in the right panel (b).

3.4. Frequency fluctuations due to environmental variations and the hydrogen bond and dangling dynamics

Frequency changes owing to the differences in the surrounding environment around the vibrational chromophore of interest are reported here considering O–D mode as the probe. As the probe molecule interacts with the environmental alterations, its frequency changes. Thus, frequency fluctuations serve as a reporter of the system dynamics⁷⁸ and can be computed using the following relationwhere C_ω(t) is the frequency–frequency time correlation function calculated by averaging over all OD bonds of interest, δω(t) is the fluctuation obtained by subtracting the average value from instantaneous frequency. The decay of fluctuations in frequency correlation is shown in Fig. 6. Table 1 displays the timescales and weight parameters obtained by fitting the decay curve with the damped oscillatory fitting function^118,119

f(t) = a₀ cos ω_ste^−t/τ₀ + a₁e^−t/τ₁ + (1 − a₀ − a₁)e^−t/τ₂ (20)

Fig. 6 The blue line shows the frequency–frequency time correlation decay of the amine nitrogen solvation shell OD modes. The frequency correlation in the hydration of the methyl group is shown by the green line, and the solvation shell of carboxylic oxygen is shown by the red line. The hazy lines show raw data, and the solid lines give the fitted data, fit by the damped oscillatory tri-exponential fitting function.

Table 1 The frequency–frequency correlation function parameters are listed below. The timescales and the frequencies are in ps and cm⁻¹ units respectively

Cω(t)	τ0	τ1	τ2	ωs	a0	a1
Solvation shell (–ND3)	0.08	0.98	5.10	95.32	0.52	0.31
Solvation shell (–CD2)	0.09	0.58	4.49	98.26	0.55	0.40
Solvation shell (–COO)	0.09	2.07	5.23	99.84	0.46	0.39

---

A faster component appears in the range of 80–90 fs timescale due to the underdamped motions of intact hydrogen bonds. We also find two other time constants, the intermediate time constant and the slowest timescale extending to a few ps. The intermediate timescales are in the range 0.5 to 2.0 ps, and the longest timespan corresponds to 4.49 ps around hydrophobic methylene and approximately 5.0 ps in hydrophilic solvation shells.

We find different types of hydrogen bonding interactions around glycine. Previously covalent-like hydrogen bonds had been confirmed by penetrating molecular-orbitals of glycine–water clusters.⁵⁴ We also calculated the continuous hydrogen bond time correlation function, S_HB(t).^120–125 The auto-correlation function is defined with the help of the population correlation function approach^115,126

Here, we compute the probability of hydrogen bond to remain continuously hydrogen-bonded from time t = 0 up to t = t. h(t) and H(t) refer to the bond parameter. h(t) = 1, when a particular pair of glycine molecules remain hydrogen-bonded at time t and zero if the hydrogen bond is broken. Similarly, H(t) = 1 when all hydrogen bonds fulfill the geometric criterion and remain in the bonded state at all times. In computer simulations, liquid water is modeled as a system consisting of N water molecules, r(t) denotes the position of all atoms at time t, and h represents the average of h[r(t)].^120,127 In the definition of the time-dependent hydrogen bond fluctuations, the denominator h in eqn (21) is presented as the normalization factor h(0)² in other studies.^87,117 Integrating S_HB(t), we compute the hydrogen bond lifetime (τ_HB).⁷⁸ In the apolar hydrophobic solvation shell of methylene of glycine, we explore the predominant presence of dangling OD. S_DH(t)¹²⁸ represent the probability of a dangling mode to remain non-hydrogen bonded in the dangling state. τ_DH gives the time period, the OD mode of water remains dangling. Time dependence of hydrogen bonds and dangling OD modes of water molecules in each of the solvation shell are shown in Fig. 7. S_HB(t) yields hydrogen bond lifetime (τ_HB) for amine–water and carbonyl–water as 1.09 and 2.17 ps, respectively. It is seen that the lifetime of O_C⋯D_W hydrogen bond is twice that of D_N⋯O_W hydrogen bond. This fact verifies the results of our frequency distribution, where the carboxyl–water hydrogen bond is found to be stronger, reflected by the double lifespan of hydrogen bond than the hydrogen bond in amine–water. Comparing the lifespan of hydrogen bonds in hydrophobic methylene solvation, 0.42 ps, the average dangling lifetime of non-hydrogen bonded OD is found to be 0.64 ps. This proves the predominance of free water molecules around the hydrophobic methylene group. The dangling lifetime of OD modes inside hydrophilic solvation is found to be 0.35 ps for aminic nitrogen and 0.46 ps for the carboxylic oxygen hydration sphere. Therefore, near hydrophilic sites, we find the water molecules are forming D_N⋯O_W and O_C⋯D_W hydrogen bonds.

Fig. 7 The upper panel (a) represents the continuous hydrogen bond correlation function, and the lower panel (b) displays the time-dependent decay of the continuous dangling OD probability function. The blue, green, and red lines represent the dynamics of hydrogen bonds and dangling OD groups inside the solvation shell of amine nitrogen, methylene, and carboxylic oxygen, respectively.

As time pass over the hydrogen bonds are broken, then either of the two things may happen: water molecules with the broken hydrogen bond stay inside the solvation shell but remain non-hydrogen bonded with the glycine molecule before it again reforms the hydrogen bond or the water molecules leave the solvation shell and diffuses away. We calculated the residence time correlation function S_R(t). By fitting the decay curve with a bi-exponential equation, we get the residence time (τ_R). Residence time measures the lifetime of water molecules inside the solvation shell. τ_R of water molecules inside N–O_W solvation shell, around methylene group and carboxylic oxygen solvation shell in glycine, are 5.18, 4.22, and 5.38 ps, respectively. We examine the similarity between the timescales of the frequency–frequency correlation function and the dynamics of hydrogen bonds and dangling OD modes. Along with the hydrophilic sites, the hydrogen bond lifetime, 1.09 ps (for H_N⋯O_W), and 2.17 ps (O_C⋯O_W) are comparable to the faster component of picosecond timescale of the decay of frequency–frequency time correlation function. For hydrophobic solvation, the dangling lifetime 0.64 ps is found similar to the intermediate component of frequency fluctuations (0.58 ps), validating the fact that the dangling water molecules are responsible for the frequency alterations in hydrophobes. The longest time constant of frequency correlation corresponds to the escape dynamics of water molecules from the solvation shell, also known as the residence dynamics.

4. Conclusions

Glycine is one of the simplest biologically relevant neutral amino acids, which is stabilized by the formation of zwitterion¹²⁹ in aqueous solution. Many structural studies on the most stable glycine conformer exist,^{42,44,130,131} but the dynamical aspects and ultrafast non-linear spectroscopy of deuterated hydroxyl groups of water molecules in the solvation shell of glycine remain unexplored. Here, we present dynamics and 2D-IR spectroscopy of water molecules of the solvation shell of glycine using first principles molecular dynamics simulations by tracing the water molecules by their distinctive behavior towards heterogeneous solvation facilitated primarily by two functional groups of amino acid: hydrophobic and hydrophilic. We performed structural calculations like radial distribution function (RDF) to understand the solvation shell demarcation.^94,104 For aminic nitrogen and carboxylic oxygen, the solvation shell criteria are N–O_W distance less than 3.48 Å and O_C–O_W distance less than 3.30 Å, respectively. For methylene, the D_C–O_W RDF peak minimum extends up to 4.2 Å. It limits the hydration shell as there is no sharp and well-defined peak because of poor interaction with the apolar hydrophobe methylene group. An isosurface plot of atomic cloud density reveals the dense population of oxygen atoms around glycine. A close view of SDF shows the uneven distribution, sparsely populated water molecules around water-hating methylene, and dense distribution along with hydrophilic sites. A contour plot of distance versus hydrogen bond angle re-verifies the RDF data.

The wavelet method⁷⁴ was used to compute the instantaneous frequency of the stretching mode of water molecules in various environments. The distribution of normalized intensity versus wavenumber was plotted for each of the solvation zones, and it was found that the average frequency around hydrophobic methylene is blue-shifted to 2550 cm⁻¹ than the hydrophilic solvation sites where the average frequency lie around 2424 cm⁻¹. However, applying the geometric criteria of the formation of hydrogen bonds, we find carboxylate–water hydrogen bonds are stronger than amine–water ones, and water molecules with free (or dangling) OD groups are found around hydrophobic methylene. To report the dynamics prevalent in glycine, the effect of environmental changes on the OD vibrational stretch was monitored from the frequency fluctuations with time. In this present study, the OD vibrational stretch serves as the local probe for 2D-IR to determine the ultrafast structure and dynamics in aqueous solvation shells in glycine. The second-order cumulant approximation was applied to generate the 2D-IR spectrum and expose the bonding patterns in aqueous solvation as a function of peak shift observable. A slight bathochromic shift is observed in the average OD stretch frequency distribution in the carboxylate solvation shell as compared to the hydration shell of amine nitrogen. Thus, the shape of the 2D-IR spectrum at each reported waiting time is similar around both the hydrophilic sites owing to the negligible peak shift, and the diagonal line width of the lobes are in accordance with the fwhm of the frequency distribution plot in the respective solvation shells. The broader diagonal width in the carboxylate solvation shell distinguishes the structural patterns and dynamical response within the two hydrophilic hydration sites in glycine. In the analysis of frequency-time correlation in each solvation shell, we observe three timescales: a rapid decay in femtosecond time owing to the librational motion of intact hydrogen bonds, the intermediate time constant corresponds to the lifetime of hydrogen bonds for hydrophilic functional groups and dangling lifetime for weakly interacting methylene solvation. Lastly, the slowest picosecond time constant reflects the time period the water molecules stay in the corresponding solvation shell. It is found that the lifetime of a carboxylate–water hydrogen bond is twice that of amine–water. Residence times around the hydrophobic methylene and around hydrophiles are 4.22 and ∼5.0 ps; respectively, due to the presence of free water molecules, they readily escape from the hydration shell in water-hating hydrophobes. Thus, the classification of aqueous glycine hydration shells in three solvation regimes enables us to explore both the structure and solvation shell dynamics occurring at the sub-picosecond to picosecond timescales and thereby, comment on the diversity that exists in the vicinity of aminic nitrogen, carboxylic oxygen, and methylene groups.

In this study, we investigate the response of water molecules towards the differential site-specific interactions of glycine. Structural aspects in solvated glycine had been reviewed in the past, but we emphasize the solvent dynamics existing in glycine solution. Our calculations help in unveiling the heterogeneous dynamics induced by the heterogeneity in structure facilitated by both the hydrophilic and the hydrophobic environments. Limitations of the experimental technique in analyzing the structure and the dynamics of vibrational spectral diffusion through spectroscopic measurements inside the solvation shell demanded support from the molecular modeling approach. Our method of calculation enables us to capture the differential solvent interactions with each of the hydrophilic and the hydrophobic sites inside the hydration shell of glycine. Hence, computation serves as an exquisite tool that can overcome the drawbacks of the experiments performing calculations and analyzing the structural, spectral, and dynamical properties inside the solvation shell. The numerous experimental and theoretical results, discussed here, provide structural evidence and exhibit the vibrational O–D stretch spectra of the water molecules in the vicinity of glycine. Our analysis of the vibrational spectrum agrees with the experimental measurements hinting the level of accuracy involved in our calculations. Structural calculations, such as the radial distribution functions (RDFs), demarcate the solvation boundary specific to a particular functional group. We classify the solvation shell of glycine into three regions based on the RDF data: the aminic nitrogen hydration shell, hydrophobic methylene solvation, and the solvation boundary around carboxylate. In each of the solvation regimes, we monitor the spectrum of the OD vibrational stretch band, time-dependent decay of frequency fluctuations, and spectral shifts due to the hydrogen bonded O–D modes and dangling water molecules inside the solvation shell. Ultrafast 2D IR spectroscopy is a powerful experimental technique that interrogates ultrafast molecular motions. In the current work, we employ the time-dependent frequencies obtained from the wavelet theorem and perform the 2D IR spectral calculations to distinguish between the two hydrophilic solvation shells of aminic nitrogen and carboxylic oxygen depending on the diagonal line width reflected in the 2D IR spectrum. We further quantify the strength of hydrogen bonds and dangling OD modes based on their average lifetime in each of the hydration shells. Here, we present a comprehensive study of vibrational spectral diffusion in aqueous glycine, exposing the dissimilarities prevailing inside specific solvation shells due to the preferential interaction of water molecules with the surrounding environment.

Conflicts of interest

There are no conflicts to declare.

Acknowledgements

The authors acknowledge financial support (EMR/2016/004965) for this work from the Department of Science and Technology, India. Aritri Biswas likes to thank Ministry of Human and Resources Development, India for the PhD fellowship.

References

1. P. M. Wiggins, Microbiol. Rev., 1990, 54, 432–449.
2. R. Cooke and I. D. Kuntz, Annu. Rev. Biophys. Bioeng., 1974, 3, 95–126.
3. B. Jayaram and T. Jain, Annu. Rev. Biophys. Biomol. Struct., 2004, 33, 343–361.
4. P. Ball, Chem. Rev., 2008, 108, 74–108.
5. F. R. Tortonda, J.-L. Pascual-Ahuir, E. Silla, I. Tuñón and F. J. Ramírez, J. Chem. Phys., 1998, 109, 592–603.
6. O. M. D. Lutz, C. B. Messner, T. S. Hofer, L. R. Canaval, G. K. Bonn and C. W. Huck, Chem. Phys., 2014, 435, 21–28.
7. M. Hashiguchi, Y. Nishi, T. Kanie, S. Ban and E. Nagaoka, Dent. Mater. J., 2009, 28, 307–314.
8. K. Leung and S. B. Rempe, J. Chem. Phys., 2005, 122, 184506.
9. B. Balta and V. Aviyente, J. Comput. Chem., 2003, 24, 1789–1802.
10. A. Chaudhari, P. K. Sahu and S.-L. Lee, J. Chem. Phys., 2003, 120, 170–174.
11. B. Balta and V. Aviyente, J. Comput. Chem., 2004, 25, 690–703.
12. S. M. Bachrach, J. Phys. Chem. A, 2008, 112, 3722–3730.
13. J. H. Jensen and M. S. Gordon, J. Am. Chem. Soc., 1995, 117, 8159–8170.
14. R. Ramaekers, J. Pajak, B. Lambie and G. Maes, J. Chem. Phys., 2004, 120, 4182–4193.
15. R. Pérez de Tudela and D. Marx, J. Phys. Chem. Lett., 2016, 7, 5137–5142.
16. H. Kayi, R. I. Kaiser and J. D. Head, Phys. Chem. Chem. Phys., 2012, 14, 4942–4958.
17. J.-Y. Kim, D.-S. Ahn, S.-W. Park and S. Lee, RSC Adv., 2014, 4, 16352–16361.
18. J. Li, T. Liu, X. Li, L. Ye, H. Chen, H. Fang, Z. Wu and R. Zhou, J. Phys. Chem. B, 2005, 109, 13639–13648.
19. P. H. Yancey, M. E. Clark, S. C. Hand, R. D. Bowlus and G. N. Somero, Science, 1982, 217, 1214–1222.
20. M. S. da Costa, H. Santos and E. A. Galinski, Adv. Biochem. Eng./Biotechnol., 1998, 61, 117–153.
21. P. H. Yancey, Am. Zool., 2001, 41, 699–709.
22. Q. Zou, B. J. Bennion, V. Daggett and K. P. Murphy, J. Am. Chem. Soc., 2002, 124, 1192–1202.
23. A. T. Russo, J. Rösgen and D. W. Bolen, J. Mol. Biol., 2003, 330, 851–866.
24. M. Ishibashi, K. Sakashita, H. Tokunaga, T. Arakawa and M. Tokunaga, J. Protein Chem., 2003, 22, 345–351.
25. P. Bruździak, B. Adamczak, E. Kaczkowska, J. Czub and J. Stangret, Phys. Chem. Chem. Phys., 2015, 17, 23155–23164.
26. D. R. Canchi and A. E. García, Annu. Rev. Phys. Chem., 2013, 64, 273–293.
27. S. S. Cho, G. Reddy, J. E. Straub and D. Thirumalai, J. Phys. Chem. B, 2011, 115, 13401–13407.
28. M. V. Athawale, J. S. Dordick and S. Garde, Biophys. J., 2005, 89, 858–866.
29. J. Krakowiak, J. Wawer and A. Panuszko, J. Chem. Thermodyn., 2013, 60, 179–190.
30. L. R. Singh and T. A. Dar, Cellular Osmolytes: From Chaperoning Protein Folding to Clinical Perspectives, Springer, 2017.
31. K. Gekko and S. N. Timasheff, Biochemistry, 1981, 20, 4667–4676.
32. S. N. Timasheff, Proc. Natl. Acad. Sci. U. S. A., 2002, 99, 9721–9726.
33. T. V. Chalikian, J. Chem. Phys., 2014, 141, 22D504.
34. A. Mukherjee, M. K. Santra, T. K. Beuria and D. Panda, FEBS J., 2005, 272, 2760–2772.
35. O. Miyawaki, M. Dozen and K. Nomura, Biophys. Chem., 2014, 185, 19–24.
36. J. M. Berg, J. L. Tymoczko and L. Stryer, Biochem., 5th edn, 2002.
37. K. A. Pikal-Cleland, J. L. Cleland, T. J. Anchordoquy and J. F. Carpenter, J. Pharm. Sci., 2002, 91, 1969–1979.
38. K. Brajter and K. Słonawska, Talanta, 1983, 30, 471–474.
39. B. López-Corcuera, A. Geerlings and C. Aragón, Mol. Membr. Biol., 2001, 18, 13–20.
40. B. Narayana Moolya, A. Jayarama, M. R. Sureshkumar and S. M. Dharmaprakash, J. Cryst. Growth, 2005, 280, 581–586.
41. S. H. Baie and K. A. Sheikh, J. Ethnopharmacol., 2000, 71, 93–100.
42. Z. Wei, D. Chen, H. Zhao, Y. Li, J. Zhu and B. Liu, J. Chem. Phys., 2014, 140, 085103.
43. G. Fischer, X. Cao, N. Cox and M. Francis, Chem. Phys., 2005, 313, 39–49.
44. S. Kumar, A. K. Rai, V. B. Singh and S. B. Rai, Spectrochim. Acta, Part A, 2005, 61, 2741–2746.
45. N. Derbel, B. Hernández, F. Pflüger, J. Liquier, F. Geinguenaud, N. Jaïdane, Z. B. Lakhdar and M. Ghomi, J. Phys. Chem. B, 2007, 111, 1470–1477.
46. A. Di Michele, M. Freda, G. Onori, M. Paolantoni, A. Santucci and P. Sassi, J. Phys. Chem. B, 2006, 110, 21077–21085.
47. Y. Hayashi, Y. Katsumoto, I. Oshige, S. Omori and A. Yasuda, J. Phys. Chem. B, 2007, 111, 11858–11863.
48. M. Sasaki, Y. Kameda, M. Yaegashi and T. Usuki, Bull. Chem. Soc. Jpn., 2003, 76, 2293–2299.
49. M. T. Parsons and Y. Koga, J. Chem. Phys., 2005, 123, 234504.
50. C. M. Aikens and M. S. Gordon, J. Am. Chem. Soc., 2006, 128, 12835–12850.
51. M. Sironi, A. Fornili and S. L. Fornili, Phys. Chem. Chem. Phys., 2001, 3, 1081–1085.
52. C. Adamo, G. E. Scuseria and V. Barone, J. Chem. Phys., 1999, 111, 2889–2899.
53. S. Liu and C. K. Schauer, J. Chem. Phys., 2015, 142, 054107.
54. Y. Shi, Z. Zhang, W. Jiang and Z. Wang, Chem. Phys. Lett., 2017, 684, 53–59.
55. L. Kócs, E. E. Najbauer, G. Bazsó, G. Magyarfalvi and G. Tarczay, J. Phys. Chem. A, 2015, 119, 2429–2437.
56. Y. Kitamura, N. Takenaka, Y. Koyano and M. Nagaoka, J. Chem. Theory Comput., 2014, 10, 3369–3379.
57. P. K. Gupta, A. Esser, H. Forbert and D. Marx, Phys. Chem. Chem. Phys., 2019, 21, 4975–4987.
58. W. Wachter, R. Buchner and G. Hefter, J. Phys. Chem. B, 2006, 110, 5147–5154.
59. L. Lupi, L. Comez, M. Paolantoni, D. Fioretto and B. M. Ladanyi, J. Phys. Chem. B, 2012, 116, 7499–7508.
60. H. Ohtaki and T. Radnai, Chem. Rev., 1993, 93, 1157–1204.
61. https://www.cpmd.org.
62. R. Car and M. Parrinello, Phys. Rev. Lett., 1985, 55, 2471–2474.
63. W. Kohn and L. J. Sham, Phys. Rev., 1965, 140, A1133–A1138.
64. J. VandeVondele, M. Krack, F. Mohamed, M. Parrinello, T. Chassaing and J. Hutter, Comput. Phys. Commun., 2005, 167, 103–128.
65. M. Eichinger, P. Tavan, J. Hutter and M. Parrinello, J. Chem. Phys., 1999, 110, 10452–10467.
66. N. Troullier and J. L. Martins, Phys. Rev. B: Condens. Matter Mater. Phys., 1991, 43, 1993–2006.
67. A. D. Becke, Phys. Rev. A: At., Mol., Opt. Phys., 1988, 38, 3098–3100.
68. C. Lee, W. Yang and R. G. Parr, Phys. Rev. B: Condens. Matter Mater. Phys., 1988, 37, 785–789.
69. S. Grimme, J. Comput. Chem., 2006, 27, 1787–1799.
70. S. Biswas, D. Chakraborty and B. S. Mallik, ACS Omega, 2018, 3, 2010–2017.
71. A. Bankura, A. Karmakar, V. Carnevale, A. Chandra and M. L. Klein, J. Phys. Chem. C, 2014, 118, 29401–29411.
72. S. Biswas and B. S. Mallik, J. Mol. Liq., 2015, 212, 941–946.
73. I. Daubechies, IEEE Trans. Inf. Theory, 1990, 36, 961–1005.
74. L. V. Vela-Arevalo and S. Wiggins, Int. J. Bifurcation Chaos Appl. Sci. Eng., 2001, 11, 1359–1380.
75. R. Carmona, W.-L. Hwang and B. Torresani, Practical Time-Frequency Analysis: Gabor and Wavelet Transforms, with an Implementation in S, Academic Press, 1998.
76. A. Panuszko, M. Śmiechowski and J. Stangret, J. Chem. Phys., 2011, 134, 115104.
77. A. Semparithi and S. Keshavamurthy, Phys. Chem. Chem. Phys., 2003, 5, 5051–5062.
78. B. S. Mallik, A. Semparithi and A. Chandra, J. Phys. Chem. A, 2008, 112, 5104–5112.
79. J. Zheng, K. Kwak and M. D. Fayer, Acc. Chem. Res., 2007, 40, 75–83.
80. M. D. Fayer, Ultrafast Infrared Vibrational Spectroscopy, CRC Press, 2013.
81. J. Zheng, K. Kwak, J. Asbury, X. Chen, I. R. Piletic and M. D. Fayer, Science, 2005, 309, 1338–1343.
82. Y. S. Kim and R. M. Hochstrasser, Proc. Natl. Acad. Sci. U. S. A., 2005, 102, 11185–11190.
83. J. Zheng, K. Kwak, X. Chen, J. B. Asbury and M. D. Fayer, J. Am. Chem. Soc., 2006, 128, 2977–2987.
84. C. J. Fecko, J. D. Eaves, J. J. Loparo, A. Tokmakoff and P. L. Geissler, Science, 2003, 301, 1698–1702.
85. J. B. Asbury, T. Steinel, K. Kwak, S. A. Corcelli, C. P. Lawrence, J. L. Skinner and M. D. Fayer, J. Chem. Phys., 2004, 121, 12431–12446.
86. J. B. Asbury, T. Steinel and M. D. Fayer, J. Phys. Chem. B, 2004, 108, 6544–6554.
87. B. S. Mallik, A. Semparithi and A. Chandra, J. Phys. Chem. A, 2008, 112, 5104–5112.
88. J. R. Schmidt, S. A. Corcelli and J. L. Skinner, J. Chem. Phys., 2005, 123, 044513.
89. S. Mukamel, Principles of Nonlinear Optical Spectroscopy, Oxford University Press, Oxford, New York, 1995.
90. P. Hamm and M. T. Zanni, Concepts and methods of 2D infrared spectroscopy, Cambridge University Press, New York, 2011.
91. P. L. Kramer, J. Nishida, C. H. Giammanco, A. Tamimi and M. D. Fayer, J. Chem. Phys., 2015, 142, 184505.
92. J. G. Kirkwood and E. M. Boggs, J. Chem. Phys., 1942, 10, 394–402.
93. M. P. Allen and D. J. Tildesley, Computer Simulation of Liquids, Oxford University Press, Oxford, 2nd edn, 2017.
94. D. Ojha and A. Chandra, J. Chem. Sci., 2017, 129, 1069–1080.
95. J. Sun, D. Bousquet, H. Forbert and D. Marx, J. Chem. Phys., 2010, 133, 114508.
96. Y. Kameda, H. Ebata, T. Usuki, O. Uemura and M. Misawa, Bull. Chem. Soc. Jpn., 1994, 67, 3159–3164.
97. M. G. Campo, J. Chem. Phys., 2006, 125, 114511.
98. S. Banerjee and H. Briesen, J. Chem. Phys., 2009, 131, 184705.
99. G. Alagona, C. Ghio and P. A. Kollman, J. Mol. Struct.: THEOCHEM, 1988, 166, 385–392.
100. J. Gao, J. Am. Chem. Soc., 1995, 117, 8600–8607.
101. S. Gnanasambandam, Z. Hu, J. Jiang and R. Rajagopalan, J. Phys. Chem. B, 2009, 113, 752–758.
102. L. N. Kuleshova and P. M. Zorkii, Acta Crystallogr., Sect. B: Struct. Crystallogr. Cryst. Chem., 1981, 37, 1363–1366.
103. M. Brehm and B. Kirchner, J. Chem. Inf. Model., 2011, 51, 2007–2023.
104. S. Biswas and B. S. Mallik, ChemistrySelect, 2017, 2, 74–83.
105. K. Furić, V. Mohaček, M. Bonifačić and I. Štefanić, J. Mol. Struct., 1992, 267, 39–44.
106. S. G. Stepanian, I. D. Reva, E. D. Radchenko, M. T. S. Rosado, M. L. T. S. Duarte, R. Fausto and L. Adamowicz, J. Phys. Chem. A, 1998, 102, 1041–1054.
107. A. Panuszko, M. Śmiechowski and J. Stangret, J. Chem. Phys., 2011, 134, 115104.
108. K. Yuan, H. Bian, Y. Shen, B. Jiang, J. Li, Y. Zhang, H. Chen and J. Zheng, J. Phys. Chem. B, 2014, 118, 3689–3695.
109. D. Chakraborty, B. S. Mallik and A. Chandra, Curr. Sci., 2014, 106, 1207–1218.
110. S. Yadav and A. Chandra, J. Phys. Chem. B, 2018, 122, 1495–1504.
111. B. Sharma and A. Chandra, J. Phys. Chem. B, 2017, 121, 10519–10529.
112. V. K. Yadav and A. Chandra, J. Phys. Chem. B, 2015, 119, 9858–9867.
113. L. De Marco, W. Carpenter, H. Liu, R. Biswas, J. M. Bowman and A. Tokmakoff, J. Phys. Chem. Lett., 2016, 7, 1769–1774.
114. I.-F. W. Kuo and D. J. Tobias, J. Phys. Chem. A, 2002, 106, 10969–10976.
115. H.-S. Lee and M. E. Tuckerman, J. Chem. Phys., 2007, 126, 164501.
116. S. Biswas and B. S. Mallik, J. Mol. Liq., 2015, 212, 941–946.
117. B. Mallik and A. Chandra, J. Chem. Sci., 2012, 124, 215–221.
118. E. T. J. Nibbering and T. Elsaesser, Chem. Rev., 2004, 104, 1887–1914.
119. H. J. Bakker and J. L. Skinner, Chem. Rev., 2010, 110, 1498–1517.
120. A. Luzar, Faraday Discuss., 1996, 103, 29–40.
121. A. Luzar and D. Chandler, Phys. Rev. Lett., 1996, 76, 928–931.
122. A. Luzar, J. Chem. Phys., 2000, 113, 10663–10675.
123. D. C. Rapaport, Mol. Phys., 1983, 50, 1151–1162.
124. A. Chandra, Phys. Rev. Lett., 2000, 85, 768–771.
125. S. Balasubramanian, S. Pal and B. Bagchi, Phys. Rev. Lett., 2002, 89, 115505.
126. R. Gnanasekaran, Y. Xu and D. M. Leitner, J. Phys. Chem. B, 2010, 114, 16989–16996.
127. A. Luzar and D. Chandler, Nature, 1996, 379, 55–57.
128. V. K. Yadav and A. Chandra, J. Chem. Phys., 2013, 138, 224501.
129. J. S. Alper, H. Dothe and D. F. Coker, Chem. Phys., 1991, 153, 51–62.
130. V. Barone, M. Biczysko, J. Bloino and C. Puzzarini, J. Chem. Theory Comput., 2013, 9, 1533–1547.
131. A. Kuffel and J. Zielkiewicz, J. Phys. Chem. B, 2008, 112, 15503–15512.

---