On-the-Fly ab Initio Semiclassical Calculation of Glycine Vibrational Spectrum.

We present an on-the-fly ab initio semiclassical study of vibrational energy levels of glycine, calculated by Fourier transform of the wavepacket correlation function. It is based on a multiple coherent states approach integrated with monodromy matrix regularization for chaotic dynamics. All four lowest-energy glycine conformers are investigated by means of single-trajectory semiclassical spectra obtained upon classical evolution of on-the-fly trajectories with harmonic zero-point energy. For the most stable conformer I, direct dynamics trajectories are also run for each vibrational mode with energy equal to the first harmonic excitation. An analysis of trajectories evolved up to 50 000 atomic time units demonstrates that, in this time span, conformers II and III can be considered as isolated species, while conformers I and IV show a pretty facile interconversion. Therefore, previous perturbative studies based on the assumption of isolated conformers are often reliable but might be not completely appropriate in the case of conformer IV and conformer I for which interconversion occurs promptly.

Introduction

Glycine, the simplest among amino acids, in addition to its evident importance in biology and medical sciences, has long played a prominent role in both experimental and theoretical chemistry. On one hand, it is the prototypical structural unit of proteins, and it has been detected in the interstellar medium[1,2] with important implications on theories about the origin of life on Earth; on the other, the presence of multiple shallow minima in the glycine potential energy surface (PES) represents a challenge for any theoretical application. One aspect shared by theory and experiment when investigating this simple but quaint amino acid is the elusiveness of its conformers. In fact, experimental microwave spectra have not always come up with consistent conclusions about the relative stability of glycine conformers,[3,4] while theoretical studies have pointed out the presence of several minima (see Figure ) with their relative stability being dependent on the level of electronic theory employed.[5−8] Eventually, at least for some conformers, equilibrium structural parameters and fundamental transition (0 → 1) frequencies have been determined experimentally and confirmed by theoretical geometry optimization and anharmonic frequency estimates.[9] Anyway, the interest in this molecule is still very vivid and it keeps on stimulating new investigations that sometimes may even involve more elusive conformers.[10]

Figure 1

Schematic representation of the glycine PES (cm–1) with its four main conformers. The energy of conformer I (Conf I) is set to 0.

A well-established protocol to identify and characterize a molecular species (especially in research fields such as astrochemistry or esobiology) lies in the analysis of its vibrational spectrum, which features high-frequency and fingerprint regions. However, the determination of vibrational frequencies of big or even medium-sized molecules is an open challenge in modern chemistry. On one hand, experimental techniques (such as IR or Raman spectroscopies) are often not conclusive and necessitate the assistance of quantum theories to assign peaks or discriminate between different hypotheses. On the other, exact theoretical and computational approaches are difficult (if not impossible) to perform when the size of the molecule exceeds a handful of atoms, and it becomes unavoidable to rely on approximations to ease the computational overhead. Furthermore, accuracy of computational results may also be adversely affected by the need to treat the underlying electronic problem by means of affordable low- or medium-level methods.

One basic theoretical procedure to estimate molecular vibrational frequencies consists in the harmonic approximation. Once the equilibrium geometry of the molecule has been optimized, and the corresponding Hessian calculated, the square roots of the nonzero Hessian eigenvalues yield the harmonic vibrational frequencies. The clear drawback of this approach is that any mode anharmonicity is neglected, and possible resonances are missed. However, in their pioneering study (based on simulations at the DFT/B3LYP and MP2 levels of theory with several types of Dunning’s correlation-consistent double-ζ basis sets[11]), Adamowicz et al.[12] have successfully employed the harmonic approximation to validate the assignments of glycine IR spectra obtained for the three most stable conformers isolated in low-temperature (13 K) argon matrixes. Anharmonic effects have been later included in other studies. For instance, Gerber et al.[13] performed an investigation of glycine conformer I based on the vibrational self-consistent field (VSCF) method. They employed a semiempirical electronic structure method with a coordinate scaling procedure to match harmonic frequencies, and further approximated the potential as a sum of single-mode and coupled two-mode contributions to ease computational costs. They got good agreement with experimental values, detecting frequencies with an average error of about 40 cm–1, and demonstrated the portability of their semiempirical potential to alanine and proline. In another work, Fernandez-Clavero and co-workers[14] focused their attention on the far-infrared (low frequency) modes of glycine conformer I. They were able to include anharmonicity by adopting a variational method based on reduced-dimensionality model Hamiltonians[15] with the remaining degrees of freedom allowed to relax. They concluded that glycine large amplitude motions require at least a four-dimensional model to be properly described, but underestimated the frequencies of other modes (treated with lower dimensional models) by as much as 150 cm–1. Hobza et al.[16] provided calculations of all 24 quantum anharmonic frequencies of conformer I. They employed a perturbative second-order vibrational perturbation theory (VPT2) approach[17,18] with electronic calculations performed at the MP2 level of theory with aVDZ basis set. Furthermore, they also got accurate estimates for six high-frequency modes (i.e., O–H, N–H2, C–H2, and C=O stretching modes) of the other three main conformers. They did not reduce the dimensionality of the problem and did not scale frequencies. Remarkably, the calculated frequencies are within 20–30 cm–1 of the experiments performed at 13 K.[12] Finally, a comprehensive set of studies concerning the four most stable glycine conformers has been recently presented by Barone and co-workers, ranging from purely spectroscopical studies[9,19] to reaction dynamics.[20] In their spectra simulations, Barone et al. adopted the VPT2 approach with variational treatment of resonances (the so-called deperturbed second-order vibrational perturbation theory (DVPT2)[21,22]) associated with a fourth-order representation of the PES computed at the DFT/B3LYP level of theory with medium sized basis sets. The calculated frequencies are in excellent agreement with experimental data from ref (12).

A potential drawback of previously outlined approaches is that the single conformers are treated as independent species in isolated energy wells. However, the shape of the glycine PES (with its already mentioned shallow wells) certainly warrants a deeper analysis of the influence of conformer interconversion on vibrational frequencies.

For this purpose, a possible approach able to explore different molecular conformations is ab initio molecular dynamics (AIMD), which yields vibrational frequencies through a Fourier transform of the dipole (or velocity) autocorrelation function. The method has become very popular in recent years with a large variety of applications such as, to cite a few relevant examples, those to methonium,[23,24] aqueous solutions,[25,26] or biomolecules.[27−30] Although AIMD is often able to provide results in agreement with experimental data and to account for classical vibrational anharmonicity, and in some approaches even quantum thermalization, it misses real-time dynamical quantum effects which are relevant for certain systems or when dealing with low temperature experiments. Differently from AIMD, we note that semiclassical methods are by construction capable of yielding quantum features starting from classical trajectories evolved on a global, full-dimensional PES.[31−50] They are based on a stationary-phase approximation to the exact Feynman quantum propagator. Over time, developments of semiclassical methods have made them more and more user-friendly, as well as a more powerful tool for molecular dynamics investigations. For instance, the original, difficult root search problem with double boundary condition can be effectively avoided by employing the initial value representation (IVR).[51] Unitarity of the semiclassical propagator is better preserved with the Herman–Kluk pre-exponential factor.[32,52,53] The chaotic nature of classical trajectories, which leads to instability, has been addressed by means of filtering techniques[101−91] and very recently addressed via rigorous approximations to the pre-exponential factor or regularization of the monodromy stability matrix.[54] Monte Carlo convergence has been facilitated by adoption of the coherent state representation,[55,56] and by introduction of techniques based on the time averaging of the oscillatory integrand.[57,58] In particular, of relevance for the present work is the multiple-coherent time-averaged semiclassical initial value representation (MC-TA-SCIVR), which permits accurate reproduction of molecular vibrational spectra by employing a tailored reference state and a few or even just a single classical trajectory.[59−64]

In addition to the possibility of exploring multiple wells through classical dynamics, semiclassical approaches can quite reliably reproduce quantum effects such as anharmonic zero-point energies (zpe), overtones, and tunneling splittings, as well as Fermi and Darling–Dennison resonances.[57,58,62,65,66] Furthermore, they can be associated with ab initio direct dynamics with excellent results,[60,61,67−70] a key feature when dealing with molecules the size of glycine or larger. In fact, even if a large variety of very precise PES fitting techniques is available (see, for instance, refs (71−83).) the computational burden to generate the hundreds of thousands of ab initio energies needed to accurately fit the PESs of such molecules is generally not affordable. For all these reasons, a semiclassical approach appears to be a straightforward choice for an investigation of “multiwell” effects, as indeed recently demonstrated by two of the authors in studying the ammonia vibrational spectrum.[66]

In this work, we employ on-the-fly MC-TA-SCIVR to estimate the vibrational frequencies of the four principal glycine conformers, and to investigate the role of conformer interconversion. The research also demonstrates that semiclassical methods are suitable to treat larger molecular systems and that they are not limited to model systems or small molecules. The paper is structured as follows: section is dedicated to the detailed description of the theoretical and computational methods employed, while in section we report results and discuss them. The paper ends with a summary and some conclusions (section ).

Theoretical and Computational Details

All electronic energy calculations including on-the-fly dynamics have been performed by means of the free NWChem suite of codes[84] at the DFT/B3LYP level of theory with aVDZ basis set. This choice was encouraged by several past studies, which have demonstrated that density functional theory with hybrid functionals and medium-size basis sets is capable of accounting very well for anharmonic effects.[9,85−87]Figure reports a sketchy representation of the glycine PES with the four lowest-energy conformers obtained upon optimization. The corresponding energy data for both the conformers and the transition states (TS) between the wells are presented in Table together with results from previous investigations. Following the footprints of Barone et al., the energy difference between conformers corrected for the calculated anharmonic zero-point energies has also been reported. We note that all calculations employing augmented Dunning’s basis sets yield very similar energy outcomes. As for geometric parameters, Table reports the main values (bond lengths, planar angles, and dihedral angles). Conformers have planar, C symmetry (conformer I (Conf I)), or nonplanar, C1 symmetry (Conf II, Conf III, Conf IV). While excellent agreement between the different studies is found for Conf I and Conf II, Conf III has sometimes been described as planar.[19] Finally, harmonic frequencies have been obtained upon Hessian diagonalization at the four equilibrium geometries, and are reported in Table . Vibrational modes are labeled in decreasing order of frequency, according to Csaszar’s notation.[6]

Table 1

Energetics (kcal/mol) of the Four Lowest-Energy Conformers of Glycine Calculated at Different Levels of Theory and Basis Setsa

	B3LYP/aVDZb	B3LYP/aug-N07Dc	MP2/aVDZd	MP4/cc-pVTZe	MP2/DZPf
Conf I	0.00	0.00	0.00	0.00	0.00
Conf II	0.49 (0.58)	0.48 (0.72)	0.54	1.05	1.02
Conf III	1.63 (1.80)	1.62 (1.60)	1.59	1.71
Conf IV	1.31 (1.08)		1.25	1.53
TS1–4	1.44
TS1–3	2.21
TS2–3	13.62

Values in parentheses refer to energies corrected for anharmonic zero point energies.

This work.

Reference (19).

Reference (16).

Reference (14).

Reference (13).

Table 2

Geometrical Parameters for the Equilibrium Configuration of the Four Conformers Calculated at the DFT/B3LYP Level of Theory with aVDZ Basis Seta

	Conf I	Conf II	Conf III	Conf IV
C1–C2	1.52	1.54	1.53	1.51
C2–O2	1.36	1.34	1.36	1.36
C2–O1	1.21	1.21	1.21	1.21
C1–N	1.45	1.47	1.45	1.46
H5–O2	0.97	0.99	0.97	0.97
C1–H3	1.10	1.10	1.10	1.10
C1–H4	1.10	1.10	1.10	1.10
N–H1	1.02	1.01	1.02	1.01
N–H2	1.02	1.01	1.02	1.01
C1–C2–O1	125.85	122.75	124.06	125.32
C1–C2–O2	111.35	113.93	113.41	111.67
O2–C2–O1	122.81	123.33	122.53	122.97
N–C1–C2	115.91	111.58	119.42	110.58
C2–O2–H5	107.18	104.89	106.50	106.99
C2–C1–H3	107.51	107.00	105.90	108.47
N–C1–H3	109.84	111.85	109.56	110.68
C2–C1–H4	107.51	107.00	105.93	105.44
N–C1–H4	109.84	111.85	109.60	114.93
H4–C1–H5	105.70	106.70	105.52	106.40
C1–N–H1	110.32	112.56	111.09	111.15
C1–N–H2	110.32	112.56	111.09	110.11
H1–N–H2	105.73	107.86	106.55	108.39
O2–C2–C1–N	180.00	1.00	1.06	162.77
O1–C2–C1–H4	123.30	58.24	56.89	105.55
O2–C2–C1–H4	–56.70	–121.65	–123.01	–72.41
O1–C2–C1–N	0.00	–179.11	–179.04	–19.27
O1–C2–C1–H3	–123.30	–55.86	–54.87	–140.81
O2–C2–C1–H3	56.70	124.24	125.23	41.23
H5–O2–C2–C1	180.00	–0.24	179.87	176.78
H5–O2–C2–O1	0.00	179.87	–0.03	–1.24
H1–N–C1–C2	58.22	117.21	59.35	38.43
H1–N–C1–H3	179.71	119.09	178.54	158.66
H2–N–C1–H3	–63.86	–122.96	–63.03	–81.22
H2–N–C1–C2	–58.22	–120.65	–59.08	158.55
H2–N–C1–H4	63.86	–3.05	63.20	39.33
H1–N–C1–H4	–179.71	–0.83	–178.37	–80.80

Bond distances are in angstroms. Angles are in degrees.

Table 3

Harmonic Frequencies and Zero-Point Energies (cm–1) for the Four Glycine Conformers (DFT/B3LYP Level of Theory, aVDZ Basis Set)a

mode	Conf I	Conf II	Conf III	Conf IV
ν1	3735	3612	3735	3740
ν2	3568	3528	3583	3594
ν3	3495	3448	3504	3501
ν4	3089	3112	3091	3070
ν5	3051	3061	3054	2957
ν6	1804	1830	1796	1808
ν7	1656	1646	1653	1618
ν8	1438	1449	1437	1474
ν9	1384	1416	1370	1435
ν10	1371	1333	1344	1318
ν11	1294	1328	1338	1252
ν12	1175	1212	1180	1209
ν13	1158	1159	1166	1137
ν14	1120	1068	1133	1105
ν15	911	915	904	1014
ν16	908	886	876	840
ν17	816	869	791	813
ν18	647	809	678	658
ν19	629	639	590	620
ν20	510	547	514	519
ν21	458	508	494	462
ν22	249	312	255	276
ν23	208	238	243	168
ν24	56	27	14	95
harmonic zpe	17364	17478	17372	17341

Mode labels are after Csaszar.[6]

On-the-fly multiple-coherent time-averaged semiclassical dynamics has been employed to calculate glycine quantum vibrational frequencies. The method has its roots in the basic Herman–Kluk (HK) version of the semiclassical propagator[52]where the integration is over the points (p0, q0) of the 2F-dimensional phase space, with F being the number of degrees of freedom of the system. C is the Herman–Kluk pre-exponential factor, and S is the instantaneous classical action calculated along the trajectory evolved from (p0, q0). Equation is completed by coherent states (|p, q⟩), which form an overcomplete basis set and are characterized by the nice property to have a Gaussian representation in both coordinate and momentum space. The explicit representation of a coherent state in coordinate space iswhere Γ is the jth element of the diagonal width matrix. In vibrational spectroscopy calculations Γ is chosen equal to the harmonic frequency of the corresponding jth vibrational mode. Hessian calculations are necessary for the evolution of the monodromy matrix elements which define the pre-exponential factor C.[57,88]

A problem that arises when applying the Herman–Kluk propagator to multidimensional systems is the difficulty to converge the integration. To help alleviate the issue, Miller and co-workers have exploited Liouville’s theorem to work out an effective and very accurate approach based on the time averaging (TA) of the integrand, making it less oscillatory and achieving convergence with orders of magnitude fewer trajectories.[57] Furthermore, the integrand of this time-averaged approach (TA-SCIVR) can be made positive-definite by implementing the so-called “separable approximation”, which consists in approximating C to its phase.[57,58]

In general, a molecular power spectrum (from which vibrational frequencies are obtained) can be calculated as the Fourier transform (FT) of the survival amplitude, or, in terms of TA-SCIVR with separable approximation, aswhere the notation is wholly similar to eq with the addition of ϕ(p0,q0), which is the phase of the Herman–Kluk prefactor; |χ⟩, which is a generic reference state; and T that indicates the total evolution time employed in performing the time average. Despite the valuable advance introduced by the TA technique, the number of trajectories needed for convergence is still too demanding for on-the-fly ab initio applications. To simulate power spectra of complex molecules starting from just a handful of trajectories or even a single trajectory, it is necessary to resort to MC-TA-SCIVR.[59,62] The basic idea of MC-TA-SCIVR rests on two pillars. First, among all trajectories that might be generated starting from random points in the 2F-dimensional phase space, those that are likely to give the major contribution to the spectrum are the trajectories with energy equal or close to the quantum vibrational one. Second, the reference state |χ⟩ can be tailored in a way to maximize the survival amplitude and get an enchanced signal in the spectrum. A general expression for the MC-TA-SCIVR reference state iswhere ε(j) are coefficients chosen to select peaks in the spectrum according to symmetry, parity, or overtone order.[61] In the MC-TA-SCIVR approach, an on-the-fly trajectory is generated for each state starting from the initial phase-space point (peq(k), qeq(k)).

We have performed single-trajectory MC-TA-SCIVR calculations with the reference state chosen (upon coherent state duplication[61]) asqeq, was derived from the conformer equilibrium geometry, and peq, was determined through the harmonic approximation peq,2 = 2ℏω(n + 1/2), where ω is the harmonic frequency of the jth mode and n are the quanta of excitation (in this work n = 0 or n = 1). ε(j) was either set equal to +1 for all degrees of freedom when the goal was to detect the zero-point-energy signal, or ε(j) = −1 for a specific mode j and ε(j′) = +1 for all others when the goal was to evaluate the frequency of the fundamental transition of mode j. Single-trajectory MC-TA-SCIVR is rather accurate in detecting the frequencies of fundamental transitions, even if the trajectory energy is somewhat shifted from their values.[59] Conversely, a proper estimate of mode overtones able to take into account anharmonicity requires that the trajectory energy be centered on the harmonic estimate of that specific mode.[62] A solid basis to the effectiveness of these single-trajectory simulations is given by the relevant work of De Leon and Heller, who demonstrated that quantum eigenvalues and eigenfunctions can be semiclassically obtained from a trajectory that not necessarily has to reside on the quantizing torus.[89]

Full-dimensional on-the-fly trajectories have been generated starting from initial conditions (i.e., atomic equilibrium positions and atomic velocities) selected according to the specific conformer and the chosen internal energy. The molecule was given no rotation with all energy concentrated on vibrational modes. Direct dynamics evolution has been performed by means of the NWChem software, which employs a velocity-Verlet integrator. Trajectories have been evolved with time steps of 10 atomic time units each, for a total of 50 000 au (approximately 1.2 ps). At each step, any instantaneous rotational and translational energy was subtracted. Energy is typically conserved with an accuracy of about 50 cm–1, while the total energy drift is about 1 kcal/mol mainly due to the rotational and translational energies subtracted. The drift value is about 2% of the total energy, but, as anticipated, semiclassical methods do not need to rely on classical trajectories run at the exact quantum energy level and the effect of the drift is limited to an enlargement of the peaks and consequently to an increase of the uncertainty associated with the results. As demonstrated in the relevant literature,[57,58,62,66] an evolution time of 25 000 au is usually long enough for a reasonable spectral resolution and for convergence of semiclassical results. As expected, we found out that for Conf I our estimate of a large majority of mode frequencies was already converged after 25 000 au. For four of the modes, though, convergence was achieved only after 50 000 au. For this reason, in the following, we will present results for all conformers based on a dynamics 50 000 au long. For every conformer, a single ab initio trajectory was run with harmonic zero-point energy, i.e. with no quanta of excitation in any vibrational modes. For conformer I, the investigation has been also performed by generating a new trajectory for each of the first 23 modes. Each of these trajectories was started with one quantum of harmonic excitation in the specific mode under examination. The lowest-frequency mode, mode ν24, is an internal torsion that in past studies has often been neglected (see, for instance, refs (9 and 19)). Following those investigations, and reckoning that when ν24 is singularly excited the total energy is anyway very close to the zero-point one, we have decided not to consider mode ν24 in the procedure based on excited trajectories and report only its estimate as obtained from zero-point-energy simulations.

MC-TA-SCIVR requires the evaluation of the Hessian matrix at each step along the trajectory to evolve the Herman–Kluk prefactor and its phase,[88] a computationally costly procedure that represents a bottleneck for the entire simulation. In the performed simulations, each Hessian calculation is about 25 times more expensive than a dynamics step. However, previous research has shown that the Hessian can be approximated by means of a compact finite-difference scheme.[88,90] In practice, the Hessian is calculated ab initio only once every few steps, and evaluated in between by means of a computationally cheaper gradient-based approximation. We have tested and employed this approach for glycine and found out that to keep accuracy and avoid artificially shifted peaks in the spectrum, for the low-frequency excitations the Hessian can be calculated ab initio every three steps, while higher-frequency ones require exact Hessian evaluations every two steps.

Finally, we have faced the monodromy matrix instability issue typical of semiclassical methods. In classical dynamics, the monodromy matrix () is defined as a four-element matrix, with each element being itself a matrix made of the partial derivatives of instantaneous positions or momenta with respect to initial positions or momenta. Eigenvalues of yield an estimate of the trajectory stability, with large real eigenvalues associated with chaotic motion. In semiclassical dynamics, the elements of the monodromy matrix are used to evaluate the Herman–Kluk prefactor, and their numerical stability is then strictly related to the reliability and feasibility of the whole semiclassical calculation. In particular, the determinant of should be unitary at all times during the simulation, but finite-precision numerical integration is not accurate enough to keep it constant in the instance of a highly chaotic trajectory. A rigorous way to keep track of the monodromy matrix stability is based on the evaluation of the determinant of the positive-definite matrix .[91] In a basic MC-TA-SCIVR simulation, we would need to stop the trajectory (and thus the simulation) as soon as the determinant of differed from unity more than 10–2. This generally happened between 15 000 and 30 000 au depending on the internal excitation (usually the higher the excitation, the faster the rejection), type of conformer (nonplanar conformers have been found to be more dynamically unstable), and the Hessian approximation adopted (which hinders preservation of determinant unitarity). To overcome the problem and be able to perform longer-time simulations, we have employed a monodromy matrix regularization technique based on the discard of the largest eigenvalue of whenever it becomes greater than a chosen threshold value (generally of the order of 103–104). Recently, this procedure has been tested on a set of model and molecular systems, demonstrating that it is able to preserve the accuracy of semiclassical results.[54] The technique prevents the semiclassical pre-exponential factor from becoming numerically unstable, thus allowing semiclassical spectra to be based on longer-time dynamics.

Results and Discussion

Our first investigation concerns the calculation of zero-point energy and fundamental frequencies for the global minimum conformer. Figure shows the peaks obtained by means of a single trajectory with internal energy equal to the harmonic zpe. The figure is divided into three parts, respectively a low-frequency part (zpe peak plus modes ν24–ν16), a medium-frequency section (modes ν15–ν7), and a high-frequency part (modes ν6–ν1). The intensity of all peaks has been normalized to that of the zpe peak, and harmonic frequencies (previously listed in Table ) are reported as dashed, vertical lines to better appreciate the anharmonicity of the quantum frequencies. Additional calculations for conformer I are performed by employing a different trajectory with specific harmonic excitation (one quantum) for each mode in the attempt to improve frequency estimates by means of trajectories with an energy content which is closer than zpe to the actual vibrational eigenvalue. Figure reports the outcome of this second approach, and Table shows a detailed numerical comparison between MC-TA-SCIVR, DVPT2, and experimental results.

Figure 2

Composition of glycine semiclassical power spectra calculated via MC-TA-SCIVR and a single classical trajectory with harmonic zero-point energy. Upper panel: zero-point energy and modes ν24–ν16. Middle panel: modes ν15–ν7. Bottom panel: modes ν6–ν1. The trajectory was evolved on-the-fly for 50 000 au. The anharmonic zero-point energy has been set to 0. Vertical dashed lines indicate the harmonic frequencies. Peak intensities have been individually normalized to that of the main peak.

Figure 3

Composition of refined MC-TA-SCIVR power spectra calculated from single trajectories with one quantum of harmonic excitation. Upper panel: zero-point energy and modes ν23–ν16. Middle panel: modes ν15–ν7. Bottom panel: modes ν6–ν1. Trajectories were evolved on-the-fly for 50 000 au. The semiclassical zero-point energy has been set to 0. Vertical dashed lines indicate the harmonic frequencies. Peak intensities have been individually normalized to that of the zpe peak.

Table 4

Comparison between Calculated Anharmonic Vibrational Frequencies for Conformer I and Corresponding Experimental Valuesa

Conf I	harmonic	MC-TA-SCIVR (zpe)	MC-TA-SCIVR (1 exc)	DVPT2b	exptc
ν1	3735	3650	3565	3575	3585
ν2	3568	3390	3395	3418	3410
ν3	3495	3395	3405	3367	3359
ν4	3089	2920	2885	2961	2969
ν5	3051	2920	2885	2947	2943
ν6	1804	1785	1765	1774	1779
ν7	1656	1675	1625	1612	1608
ν8	1438	1410	1380	1435	1429
ν9	1384	1375	1330	1387	1405
ν10	1371	1345	1330	1353	1340
ν11	1294	1300	1250	1286	1297
ν12	1175	1165	1150	1164	1166
ν13	1158	1120	1105	1144	1136
ν14	1120	1100	1090	1103	1101
ν15	911	905	900	863	883
ν16	908	875	845	907	907
ν17	816	795	785	802	801
ν18	647	645	600	603	619
ν19	629	625	625	633	615
ν20	510	490	485	494	500
ν21	458	460	470	461	458
ν22	249	260	275	255	250
ν23	208	220	180	203	204
ν24	56	85	–	–	–
zpe	17364	17160	17160	–	–

Data (cm–1) are reported for MC-TA-SCIVR (single trajectory with zero-point energy (zpe), or with one quantum of harmonic excitation (1 exc)), deperturbed second-order vibrational perturbation theory (from ref (9)), and experiment at low temperature (12–15 K) (from refs (12 and 92)).

From ref (9).

From refs (12 and 92).

Peaks in Figures and 3 are well resolved, and frequencies are generally in qualitative good agreement with previous DVPT2[9,19] and experimental works.[12,92] Semiclassical methods in applications to much smaller molecules such as ammonia[66] or H2O[54] have shown small mean absolute errors (38 and 20 cm–1 respectively) with respect to exact quantum calculations. However, exact quantum mechanical results are not available for glycine and a straightforward quantitative comparison between the reported data is not possible, since DVPT2 studies employed a different basis set from ours, and experimental spectra always come with some uncertainties in their interpretation. For these reasons comparisons are intended to be mainly qualitative. To assign error bars to our glycine semiclassical calculations, we calculated the full width at half-maximum (fwhm) of the peaks presented in Figures and 3. For the zpe simulation the average fwhm is 50 cm–1, while it increases to 99 cm–1 for estimates based on mode-specific excited trajectories. This means that the semiclassical estimates reported in Table have average uncertainties of ±25 and ±49 cm–1 respectively. The higher uncertainty for estimates based on excited trajectories is an expected drawback since coupling and energy flow between different vibrational modes are more likely when the energy is higher with the result of broadening the spectroscopic signal. One relevant exception concerns the O–H stretch which has an uncertainty of ±18 cm–1 in the refined simulation. Complete data for each peak can be found in the Supporting Information. From a further inspection of Figures and 3, we note that higher anharmonicity is usually found for the high-frequency modes involving the O–H (mode ν1), N–H2 (modes ν2 and ν3), and C–H2 (modes ν4 and ν5) stretches, in agreement with previous studies.[16]

We then move to a comparison of vibrational frequencies regarding conformers II, III, and IV. Experimental data are not available for the last, and they are also scarce for the former two. As reported in Table , and with the previously illustrated caveat about data comparisons, MC-TA-SCIVR results (shorthanded as “MC” in Table ) are in good agreement with DVPT2 and the experiment, especially for the mid-frequency modes. Average uncertainties of our semiclassical frequencies for conformers II, III, and IV are ±23, ±26, and ±45 cm–1, respectively. Conformer IV with its very shallow well confirms to be difficult to investigate (see DVPT2 data in Table ), while detailed fwhm values can be found in the Supporting Information. We note that the two symmetry groups to which the four glycine conformers belong (C or C1) have only monodimensional representations, but in spite of this, a few modes are degenerate when calculated semiclassically. This is actually an accidental outcome of our calculations that involves modes which are very similar regarding both frequency and type of motion. A couple of comments support this claim. On one hand, also other approaches are not totally immune from this accidental degeneracy drawback. For instance, the two C–H2 stretches (mode ν4 and ν5) are nearly degenerate in DVPT2 simulations (see data for conformers II and III in Table ). On the other hand, semiclassical approaches in general, and specifically MC-TA-SCIVR ones, are fully able to properly and correctly account for nuclear symmetry, as demonstrated in several previous works.[57,58,61,66,67]

Table 5

Comparison between Harmonic (harm), MC-TA-SCIVR (MC), DVPT2, and Experimental Estimates of Fundamental Frequencies for Glycine Conformers II, III, and IVa

	Conf II				Conf III				Conf IV
	harm	MC	DVPT2b	exptc	harm	MC	DVPT2b	exptc	harm	MC	DVPT2d
ν1	3612	3439	3440	3410	3735	3636	3576	3560	3740	3649	3579
ν2	3528	3409	3373		3583	3492	3425	3410	3594	3484	3441
ν3	3448	3325	3235	∼3275	3504	3414	3371		3501	3415	3361
ν4	3112	2935	2955		3091	2832	2933		3070	2987	2956
ν5	3061	2941	2953	2958	3054	2856	2932	2958	2957	2854	2866
ν6	1830	1807	1821	1790	1796	1776	1779	1767	1808	1789
ν7	1646	1609	1618	1622	1653	1632	1641	1630	1618	1609
ν8	1449	1417	1431	1429	1437	1428	1417	1429	1474	1429
ν9	1416	1357	1377	1390	1370	1350	1339		1435	1381
ν10	1333	1309	1322		1344	1368	1318	1339	1318	1279
ν11	1328	1309	1294		1338	1308	1339		1252	1219
ν12	1212	1165	1169	1210	1180	1164	1153		1209	1183
ν13	1159	1129	1145		1166	1179	1128	1147	1137	1112
ν14	1068	1051	1044		1133	1128	1105	1101	1105	1104
ν15	915	901	899	911	904	900	895		1014	979
ν16	886	847	851	867	876	816	827	852	840	837
ν17	869	829	840		791	810	770	777	813	763
ν18	809	757	805	786	678	666	656	644	658	637
ν19	639	607	633		590	588	604		620	589
ν20	547	541	543		514	531	499		519	523	502
ν21	508	493	502		494	510	488		462	467	467
ν22	312	319	299	303	255	258	260		276	252	278
ν23	238	235	229		243	216	224		168	196	156
zpe	17478	17221			17372	17216			17341	17033

Single trajectories with zero-point energy and 50 000 au long have been employed for the MC-TA-SCIVR calculations. Data are in cm–1. Experimental data at 12–15 K.

From ref (19).

From refs (12 and 92).

From ref (9).

A peculiar feature of MC-TA-SCIVR is the capability to detect quantum overtones at no additional cost. In fact, the reference state in eq depends on the choice of the ε(j) parameters. As described in section , when a trajectory is given an initial single harmonic excitation in the jth mode, then the reference state is built with a value ε(j) set equal to −1. In this way, not only the peak relative to the fundamental transition of the jth mode is selected, but it is also possible to identify spectral signals associated with an odd number of excitations in that mode. As a result, in Figure we show four representative examples of second order overtone peaks (i.e., at frequencies corresponding to a triple excitation of the mode) .

Figure 4

Fundamental and second overtone frequencies obtained via MC-TA-SCIVR for conformer I from trajectories with one quantum of harmonic excitation. Peak frequencies are (cm–1): 201 = 485, 203 = 1400; 191 = 625, 193 = 1865; 151 = 900, 153 = 2705; 141 = 1090, 143 = 3265.

To complete this section, we focus on the possibility to have a substantial “multiwell” effect on vibrational frequencies. For a molecule like glycine, the conformers should not be considered as isolated ones due to their relatively low interconversion barriers. Modes related to large-amplitude motions along the interconversion path, in fact, could in principle be largely influenced by the presence of multiple wells. For this reason, precise perturbative treatments should rely on a global surface, which, however, is computationally prohibitive for a molecule the dimension of glycine. One advantage of MC-TA-SCIVR over perturbative methods is that it naturally includes the “multiwell” effect if the classical trajectories visit several wells. To assess the importance of the effect in the present study, we developed a procedure based first on trajectory investigation, and then on a comparison between frequencies calculated from trajectories of different lengths. Figure shows the equilibrium normal-mode coordinates for most vibrations. For some of the modes, the equilibrium coordinate does not change moving from a conformer to another. We interpret these modes as those not linked to any interconversion paths. On the contrary, there are modes that have substantially different equilibrium coordinates for different conformers. These are likely related to the interconversion path and, when excited, have a chance to lead to a change in glycine conformation. Following this idea, we have run and analyzed trajectories with single harmonic excitation starting from all conformers up to the 50 000 atomic time units used in the semiclassical calculations, and then even longer to 200 000 au. In the 50 000 au range, we found no interconversion between conformer II and conformer III, nor between conformer III and conformer I. The interconversion between conformer I and conformer IV was pretty facile, instead, sometimes taking place even within an evolution time of just 25 000 au (as in the case of trajectories starting from conformer IV, and of modes ν2, ν4, ν5, ν6, ν8, ν9, ν17, ν18, and ν23 when excited from conformer I). For some other modes (ν1, ν3, ν7, ν13, ν14, ν19, ν21, and ν22 starting from conformer I), the I–IV interconversion happens after an evolution time larger than 25 000 but shorter than 50 000 au. Consequently, our semiclassical frequency estimates at t = 50 000 au are based either on isolated trajectories within the starting conformer well, or on trajectories that experience (I ↔ IV) conformer interconversion, thus differentiating our approach from VPT2 calculations based on the assumption of isolated wells.

Figure 5

Mass-scaled normal mode equilibrium coordinates for the four conformers. For clarity of the figure, some modes are not reported since they are flat around Qeq = 0.

The presence of a substantial “multiwell” effect should reveal itself in two main ways. First, the mode frequencies should be independent from the geometry of the starting conformer; second, a substantial modification of the frequencies is expected when multiple wells are visited compared to isolated well estimates. Actually, we found different semiclassical frequency estimates by starting our calculations from different conformers as seen in Tables and 5. To check out the second clue for the presence of a “multiwell” effect, we investigated longer time classical dynamics (up to 200 000 au) and observed that eventually interconversion among all the conformers does take place, as expected, but it is a rare event. At such times a semiclassical calculation is not affordable for glycine with ab initio on-the-fly dynamics. For this reason we have moved to classical dipole–dipole simulations. These are based on short-time (up to about 5 ps) and not thermalized classical trajectories, but they may nevertheless be helpful to check out if the interconversion events (even if rare) are able to produce any frequency shifts due to a substantial “multiwell” effect. Table reports the outcome of the dipole–dipole simulations, calculated at different evolution times, for some relevant modes that undergo conformer interconversion. Trajectories have been started with single harmonic excitation in the interested mode and evolved up to 200 000 au. At short times (up to 50 000 au) conformer interconversion has still to take place. By comparing these transient frequency values at different times, we can see that the value of the frequency peaks is not drastically changed even when trajectories are longer. Changes are within a few wavenumbers. Furthermore, by comparing the same classical mode frequency calculated starting from different conformers that interconvert during the longer dynamics, we cannot see any merging of the frequency, in agreement with the semiclassical estimates. Also experiments performed at low temperature (13 K)[12] and reported in Tables and 5 yield different frequency values for the different conformers. A regime in which a single frequency for each mode is detected independently from the starting conformer is expected to be reached at much higher energies than those required by a semiclassical investigation, since conformer interconversion then will not be a rare but rather an instantaneous and frequent event.

Table 6

Classical Vibrational Frequencies (cm–1) for Some Relevant Modes Undergoing Conformer Interconversion

	ν6	ν16	ν17
Conf I (50 000 au)	1772	847
Conf I (100 000 au)	1772	844
Conf I (200 000 au)	1773	842
Conf II (50 000 au)			754
Conf II (100 000 au)			756
Conf II (200 000 au)			757
Conf III (50 000 au)	1756		773
Conf III (100 000 au)	1758		769
Conf III (200 000 au)	1760		765
Conf IV (50 000 au)		826
Conf IV (100 000 au)		824
Conf IV (200 000 au)		821

Summary and Conclusions

We have presented a semiclassical study of glycine vibrational frequencies. The molecule has 24 vibrational modes which make it a particularly challenging system for an application of semiclassical dynamics. MC-TA-SCIVR has been employed to produce spectra for the four lowest-energy conformers. As anticipated in the Introduction, the conformer equilibrium geometries are difficult to characterize, with optimized structures significantly dependent on the level of electronic theory and basis set employed. This aspect is confirmed, for instance, in the case of conformer III. We found a nonplanar equilibrium geometry, in agreement with ref (16), where MP2 with aVDZ basis set was employed, but in disagreement with ref (19) based on DFT/B3LYP and N07D basis set, which returned a perfectly planar geometry.

A remarkable feature of our study is that spectra have been obtained from single on-the-fly ab initio trajectories, i.e. by direct dynamics. Furthermore, results based on zero-point-energy trajectories have shown good accuracy in the whole range of calculated frequencies, and allowed us to examine conformers II, III, and IV with much reduced computational overhead. The reason for such flexibility lies in the Gaussian shape of the reference state, which permits having sufficient survival amplitude overlap within a broad shell centered on the energy of the trajectory. For this same reason energy drifts in classical trajectories do not prevent obtaining resolved semiclassical spectra, but have an influence that is limited to the uncertainty of frequency estimates. The latter may be evaluated from the fwhm values of peaks in the spectra. The average fwhm values we found confirm that conformer IV is the most elusive one as predictable by looking at its very shallow well.

A novelty we have brought with this work is that our glycine investigation, differently from previous studies on the same molecule, is based on a time-dependent approach with the vibrational spectrum obtained via Fourier transform of the semiclassical time-dependent wavepacket correlation function. On one hand, this introduces the issue of zero-point-energy leak. It is a difficult problem to solve, and past studies have come up with contrasting conclusions on the most efficient way to address it. Among several techniques developed for this purpose we mention the “quantum kicks” model,[93,94] adiabatic switching,[95,96] and dynamical normal-mode analysis.[97] In the Supporting Information we present an investigation of zero-point-energy leak based on a normal-mode approximation. On the other hand, the time-dependent MC-TA-SCIVR approach permits accounting for the influence of the presence of many accessible potential wells in the molecular PES. In our study, when conformer interconversion is a rare event taking place at long times, there is no appreciable “multiwell” effect on the semiclassical frequency values. Conversely, there are cases (conformer IV and some of conformer I modes) for which the I ↔ IV interconversion happens almost instantaneously (even before 25 000 au) and might be a factor in the frequency estimate. For these reasons, we conclude that previous approaches based on the assumption of isolated wells are often reliable and justified, but the isolated-well approximation they employ could be not very appropriate in the case of conformer IV and partially conformer I.

Our semiclassical calculations were based on 50 000 au long trajectories. However, in a standard application of MC-TA-SCIVR to glycine, instability of the pre-exponential factor would arise well before 50 000 au, hampering the whole simulation. To overcome this issue, we employed a monodromy matrix regularization procedure that, after discarding the influence of chaoticity on the monodromy matrix, allowed us to perform long-time simulations. A further advantage over perturbative approaches brought by MC-TA-SCIVR is the possibility to detect overtones efficiently and without additional simulations. We have reported a few examples of this. On the contrary, overtones are difficult to determine with perturbative methods because of the potential failure of these approaches in the presence of even accidental (quasi-) degeneracies. To overcome the issue, one should integrate perturbative methods with variational ones (see, for instance, ref (19)) at the cost of complicating and hybridizing the approach.

Finally, we want to remark that the present work will have a natural development in the investigation of vibrational frequencies of glycine in its protonated and zwitterionic forms. Results for neutral glycine are encouraging, because they have been obtained with very good accuracy through a computationally cheap approach. A study of vibrational motion of protonated glycine will be highly valuable. The species is found in solution,[98−100] and it is important in biological processes where it acts as an intermediate in many biochemical reactions. MC-TA-SCIVR studies of protonated and zwitterionic glycine are foreseen, and they might also provide an additional, technical “byproduct”. In fact, trajectories for protonated glycine can be either generated with an on-the-fly ab initio approach or evolved through one of the many commercial force fields available for biological systems. Hence, results based on the ab initio calculations will also serve as a benchmark for assessing the force field accuracy.