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

We present a ﬁ rst 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 (cid:1) 2416 cm (cid:3) 1 . We observe our calculated average O – D vibrational stretch frequency in aqueous glycine, (cid:1) 2416 cm (cid:3) 1 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 (cid:3) 1 , respectively. Contrary to the frequency range in hydrophilic solvation, the average frequency around hydrophobic methylene is (cid:1) 2550 cm (cid:3) 1 , and the hydrogen bond and dangling lifetimes are 0.42 and 0.64 ps, respectively. These data con ﬁ rm 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 di ﬀ erences in local environments surrounding the hydrophilic species. Our calculations of the frequency distribution, two-dimensional infrared (2D-IR) spectrum, the decay of frequency ﬂ uctuations, the lifetime of hydrogen bonds, and dangling modes unveil the distinct behavior of water molecules in solvation regions.


Introduction
In biological systems, molecular-level understanding of solutesolvent interactions is considered to be impossible without explicit aqueous solvent molecules. [1][2][3][4] The basic building blocks of proteins, the amino acids are organic compounds with the general formulae NH 2 -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, 7 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; 12 at least two water molecules 13 were required to stabilize the zwitterionic conformer. However, a minimum of four water molecules was required for stabilizing the zwitterionic state of glycine. 14 Recently, eight or nine water molecules were theoretically predicted to stabilize glycine solvation. [15][16][17] Like many biomolecules such as proteins contain non-polar hydrophobic and polar hydrophilic regions, 18 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][20][21] However, the mechanism of stabilization of the globular proteins by osmolytes [22][23][24] is still debatable 25 within the experimental and modeling approaches. [26][27][28][29][30] One hypothesis [31][32][33][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. 35 Glycine serves as a protein precursor, 36 buffering agent 37 in antacids, analgesics, and antiperspirants. It is used as a metal complexing agent, 38 and inhibits postsynaptic potential acting as a neurotransmitter in the central nervous system. 39 Glycine exhibits nonlinear optical properties 40 and is widely applicable in the pharmaceutical industry. 41 Several studies addressed the solvation of glycine in water 17,42 and its relevance in the eld of biochemistry to investigate the process of protein-water interactions, and protein folding. 18 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][44][45] Raman, 46 dielectric relaxation spectroscopy, 47 neutron diffraction, 48 and calorimetry 49 as well as the theoretical studies like ab initio simulations using different basis sets. 10,13 Furthermore, Monte Carlo 50 and molecular dynamics (MD) 51 studies served crucial for providing the evidence for the structural changes brought about by the solvation and hydrogenbonding networks. Hydrogen bonds play a vital role in explaining the interactions in organic-water molecules, DFT PBE0 level rst 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. 54 In 1990, Basch and Stevens rst published the ab initio calculations on the different possible structures of glycine$H 2 O. 55 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 shis induced by solvent water molecules around glycine were estimated by the vibrational frequency analysis (VFA) with the free energy Hessian (FE Hessian). 56 Experimentally, it was a challenging task to probe the water dynamics in the respective solvation shells in the vicinity of the solute. 57 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.

Methodology
First principles molecular dynamics (FPMD) study was performed using CPMD simulation package 61 employing the Car-Parrinello method. 62 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 63 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 66 to treat core electrons. BLYP 67,68 functional with dispersion correction 69 (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. 70 The calculation of oxygen-oxygen radial distribution function (RDF) using the same functional reected perfect tetrahedral hydrogen bond associations. 71 The formamidewater system simulated using BLYP-D2 density functional adopting the same protocol provided an good agreement between the theoretical calculations and the experimental results. 72 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 ctitious 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 rst 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 Ø denes 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 uctuating 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 where r(t) is the real part representing instantaneous uctuations 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: j a,b (t) has compact support or decays to 0 rapidly for at t / AEN. 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 jrepresents the Fourier transform of j. Wavelet transform considers a timefrequency representation f; the parameter b shis the wavelets so as to store the local information of f in L j f(a,b) at time t ¼ b. L j 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 timevariable b. The frequency variations are regulated by the dilation parameter a. The wavelet j a,b shrinks for |a| ( 1, mostly in the higher frequency range and for |a| [ 1, j a,b dilates with the frequency content in the low frequency range. Aer 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 75 is used for analysis Expressing eqn (4) in terms of p e 2pilðtÀb=aÞ e ÀðtÀb=aÞ 2 =2s 2 (5) s ffiffiffiffiffiffi 2p p e 2pilðtÀb=aÞ e ÀðtÀb=aÞ 2 =2s 2 (6) j must be real and should satisfy the admissible condition, 0\c j ¼ 2p Ð þN ÀN ð|j ðuÞ| 2 =|u|Þdu\N, and l ¼ 1 and s ¼ 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 76 of timedomain 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. 86 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. 87 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][89][90] in our calculation of the 2D-IR spectrum. In this approximation, it is assumed that the uctuating dynamics describe a Gaussian process and is valid for observing the peak shi. However, in strongly hydrogen-bonded systems like liquid water, non-Condon line shape is more realistic as it takes into account non-Gaussian frequency uctuations and environmental impacts on the transition dipole (non-Condon effects). 88 Theoretical calculations of the 2D-IR spectroscopic observables involve response function formulation. 88,90 Three ultrashort IR pulses with wave vectors k 1 , k 2 , and k 3 interacts with three input electric elds E 1 , E 2 , and E 3 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][89][90] Ignoring, the excited state and rotational lifetime effects, the equations for response functions in secondorder cumulant approximation is given by 10 4 exp(Àihu 10 it 3 À ihu 10 it 1 ) Â exp(Àg 1 (t 1 ) À g 1 (t 2 ) À g 1 (t 3 ) + g 1 (t 2 + t 1 ) + g 1 (t 2 + t 3 ) À g 1 (t 1 + t 2 + t 3 )), (9) R 6 (t 3 ,t 2 ,t 1 ) ¼ Àm 10 4 k 2 exp(iDt 3 ) Â exp(Àihu 21 it 3 À ihu 10 it 1 ) Â exp(Àg 1 (t 1 ) À g 2 (t 2 ) À g 3 (t 3 ) + g 2 (t 2 + t 1 ) + g 2 (t 2 + t 3 ) À g 2 (t 1 + t 2 + t 3 )). (10) The time delay between the rst and second pulses is referred to as t 1 , between second and third pulses as t 2 , t 2 ¼ T w (waiting time), and the time aer the third pulse as t 3 . g 1 (t), g 2 (t) and g 3 (t) are the time correlations of the uctuating transition frequencies given by the following expressions 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 2 ¼ T w ) expressed as follows [88][89][90] where, for i ¼ 1 to 3 and for i ¼ 4 to 6, respectively. At initial T w , the nal 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 nal frequencies are different from the initial value, and the frequencies decorrelate with time. Thus, with time more spectral diffusion occurs, resulting in frequency decorrelation quantied by a symmetrical 2D-IR line shape. 91 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 uctuating frequencies.

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 by where 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, 93 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 rst peak minimum is taken as the distance cut-off for dening 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 94 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 The value of the number integral at the rst 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 rst 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 96 and Car-Parrinello 8 method reported 3 to 4 solvation water molecules around the ammonium group in glycine. The latter AIMD simulation 8 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 rst 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 rst solvation sphere around the ammonium group, match with the previous CPMD study of glycine intramolecular proton transfer in water. 8 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 dened because of its shallowness compared to the rst peak minimum. Thus, the second hydration layer around aminic nitrogen lies undened as in preceding reports. 95 The RDF of D C -O W pair around methylene functionality reects very weak interaction with the nearby water as there is no sharp peak maximum. The solvation boundary is not well-dened with minimum extending up to 4.3Å due to the water-repelling nature of CD 2 moiety that inhibits the approach of water molecules. Moving the attention towards the carboxylates, we calculated O C -O W RDF, and the rst peak maximum appears at 2.85Å matching with the previously reported values, 97,98 and the rst solvation shell is at a distance of 3.30Å, (at the rst 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 rst 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 rst 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. 8 However, both the values were found to be smaller than the hydration numbers 5.3 and 7 predicted by simulations using empirical force eld 99 and OPLS force elds 100 respectively. Classical molecular dynamics study 101 with AMBER ff03 force eld showed the rst 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. 102 The hydration numbers reported by the study of glycine/water mixtures 101 in crystal/solution environments around amino and carboxylic groups were 2.97 and 5.38, respectively. Neutron diffraction measurements 48 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 nding the oxygen and deuterium atoms around the glycine molecule in a threedimensional 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 103 comes from contour lines passing through the areas with the same probability in topographic maps. SDF data of water molecules around ND 3 , CD 2 , 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 gure 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 gure, we nd a small fraction of water molecules around the apolar methylene group; this justies 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 76 and ab initio molecular dynamics simulations 95 revealed carboxylate anion underwent nonequivalent interaction with hydrating water via two oxygen atoms.
In the aqueous solution, we nd 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 gure 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 [q 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 soware package. 103 In the contour plot shown in the upper panel, we nd 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 nding 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 rst minimum (3.30Å) of the O C -O W RDF peak. We also nd a light green patch along 45 to 75 , but it is a region that is less probable and not so well dened.

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 3 (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 104 based on RDF calculations, and the solvation boundary was dened 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 94 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 2 . 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 74 and was found to be 2416 cm À1 . Raman spectrum recorded a bifurcated D 2 O stretching region with peaks at 2190 and 2390 cm À1 and revealed the former band had more intensity than the later with the increase in the concentration of dissolved molecules. 105 IR spectra of deuterated glycine derivatives in low-temperature argon matrix exhibited spectral features in three most stable glycine conformers for the rst time and assigned the frequencies to the stretching and bending vibrations. 106 These matrix-isolated IR spectra of the deuterated glycine portrayed a spectral signature at 2399 cm À1 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 À1 and helped to identify the IR spectral characteristics of the glycine conformers in the experiment. Our calculated value of the average frequency   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 À1 units respectively of all the OD bonds in deuterated glycine $2416 cm À1 is close to the computed OD stretch frequency 2411 cm À1 . 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. 107 It is difficult to determine the frequencies inside each of the hydration regions and outside it by applying the solvation shell conditions experimentally, 108 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 À1 . Inside aminic nitrogen solvation shell, OD modes that form D N /O W hydrogen bond has an average frequency of 2418 cm À1 , exhibiting a bathochromic shi of 7 cm À1 . We also observe a higher value of the average frequency 2550 cm À1 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 À1 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. 109 Stretch frequency distribution in aqueous ionic solutions like carbonate, 110 and phosphate ion 111 showed a broad spectral feature lying in the range 2500-2700 cm À1 representing dangling OD groups. Inside the solvation regime of carboxylic oxygen, the average frequency is 2423 cm À1 , and the solvation shell O-D stretch frequency of aminic nitrogen is 2425 cm À1 . 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 À1 , respectively, inside amine nitrogen solvation shell; moreover, 2341 and 2346 cm À1 , respectively, inside carbonyl oxygen hydration shell. 112 Till now, the observed frequency shis are within AE7 cm À1 for hydrophilic solvation and $25 cm À1 around hydrophobic -CD 2 due to the existence of dangling OD modes of water near it. Carboxylate-water hydrogen bonds with average frequency 2397 cm À1 are found to be stronger than amine-water hydrogen bonds having a value of 2418 cm À1 . 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. 76 Bulk OD modes with frequency 2410 cm À1 are in agreement with the computed average O-D stretch frequency value ($2380 cm À1 ) of heavy water, 78 and the experiments reported vibrational stretch of the OD band centered at $2500 cm À1 . 113 The discrepancies are a result of the difference in electronic structure calculations with different density functionals, electronic ctitious mass, and nite basis set cut-offs. 114,115 Comparing the bulk average O-D stretch frequency ($2410 cm À1 ) with the average frequency of OD modes hydrogen-bonded to carboxylate ($2397 cm À1 ), and ammonium ($2418 cm À1 ); we observe that the OD groups hydrogen-bonded to carboxylate are red-shied with respect to the water molecules hydrogen-bonded to other water at bulk, and the hydrogen-bonded OD in amine solvation is slightly blue-shied relative to bulk. First principles molecular dynamics simulations of aqueous formamide solvation revealed carbonyl-water hydrogen bonds are stronger than hydrogenbonded water molecules at bulk, and amine-water hydrogen bonds are weaker than bulk based on the frequency shis. 116 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 shi of 30 cm À1 was reported ongoing from bulk towards amine hydration shell. 104 The OD modes within the solvation shell cut-off in amine hydration (2425 cm À1 ) show a similar 15 cm À1 frequency shi corresponding to the bulk OD stretch (2410 cm À1 ) in glycine. In glycine-water system, bulk water molecules exhibit a red shi in the average vibrational stretch frequency as compared to the carboxylate hydration shell and this shi is comparable to the 34 cm À1 shi observed ongoing from the carbonyl oxygen solvation towards bulk in aqueous acetone. 117

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 3 + and negatively charged COO À in aqueous zwitterionic glycine. The peak on the u 1 -u 3 diagonal originates from stimulated emission and ground-state bleach from n ¼ 0 / 1 energy levels, and the off-diagonal peak represents vibrational excited state absorption for n ¼ 1 / 2 transition, shied 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 halfmaximum (fwhm) of the O-D stretch frequency distribution inside the solvation shells in glycine.

Frequency uctuations 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 uctuations serve as a reporter of the system dynamics 78 and can be computed using the following relation where C u (t) is the frequency-frequency time correlation function calculated by averaging over all OD bonds of interest, du(t) is the uctuation obtained by subtracting the average value from instantaneous frequency. The decay of uctuations in frequency correlation is shown in Fig. 6. Table 1 displays the timescales and weight parameters obtained by tting the decay curve with the damped oscillatory tting function 118,119 f(t) ¼ a 0 cos u s te Àt/s0 + a 1 e Àt/s1 + (1 À a 0 À a 1 )e Àt/s2 (20) A We nd different types of hydrogen bonding interactions around glycine. Previously covalent-like hydrogen bonds had been conrmed by penetrating molecular-orbitals of glycinewater clusters. 54 We also calculated the continuous hydrogen bond time correlation function, S HB (t). [120][121][122][123][124][125] The autocorrelation function is dened with the help of the population correlation function approach 115,126 S HB ðtÞ ¼ hhð0ÞHðtÞi h Here, we compute the probability of hydrogen bond to remain continuously hydrogen-bonded from time 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 fulll 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 denition of the timedependent hydrogen bond uctuations, the denominator h in eqn (21) is presented as the normalization factor h(0) 2 in other studies. 87,117 Integrating S HB (t), we compute the hydrogen bond lifetime (s HB ). 78 In the apolar hydrophobic solvation shell of methylene of glycine, we explore the predominant presence of dangling OD. S DH (t) 128 represent the probability of a dangling mode to remain non-hydrogen bonded in the dangling state. s 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 (s 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 veries the results of our frequency distribution, where the carboxyl-water hydrogen bond is found to be stronger, reected by the double lifespan of hydrogen bond than the hydrogen bond in aminewater. 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 nd the water molecules are forming D N /O W and O C /D W hydrogen bonds.
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 tting the decay curve with a bi-exponential equation, we get the residence time (s R ). Residence time measures the lifetime of water molecules inside the solvation shell. s 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 frequencyfrequency correlation function and the dynamics of hydrogen bonds and dangling OD modes. Along with the hydrophilic sites, the hydrogen bond lifetime,

Conclusions
Glycine is one of the simplest biologically relevant neutral amino acids, which is stabilized by the formation of zwitterion 129 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 rst 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-dened 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-veries the RDF data.
The wavelet method 74 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-shied to 2550 cm À1 than the hydrophilic solvation sites where the average frequency lie around 2424 cm À1 . However, applying the geometric criteria of the formation of hydrogen bonds, we nd 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 uctuations 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 shi observable. A slight bathochromic shi 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 shi, 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 reects 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 aminewater. 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 classication 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-specic 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 specic 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 uctuations, and spectral shis 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 reected 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 specic solvation shells due to the preferential interaction of water molecules with the surrounding environment.

Conflicts of interest
There are no conicts to declare.