Computing the frequency fluctuation dynamics of highly coupled vibrational transitions using neural networks

The extensive use of lithium ion batteries in consumer electronics and electric vehicles has drawn significant attention to the study of the inner workings of this energy storage technology at the molecular level.^1–6 In particular, the structure and dynamics of the lithium ion electrolyte have drawn considerable experimental and theoretical interest^7–9 since many safety issues are related to this critical component, but its molecular level understanding has started to emerge very recently.^6,10–12 The solvation structure and motions of lithium ions in organic-based electrolytes have been investigated using x-ray scattering¹³ and nuclear magnetic resonance,^14–18 along with Raman and IR spectroscopies.^18–23 The spectral intensities in the Raman and IR bands have been used to determine the ratios of free to bound organic carbonates or anions and therefore to estimate the ion–solvent and ion–ion interactions.^{18,19,21,22,24,25} However, these conventional methodologies have not been able to probe the local structure surrounding the lithium ion (solvation shell) because they lack the time resolution needed to investigate the fast motions in the electrolyte solution, which are usually on the timescale of femtoseconds to picoseconds.

Recently, time resolved infrared spectroscopies have been used to close the knowledge gap in the structure and the ultrafast dynamics of the ionic and molecular components in the lithium ion electrolyte. In particular, studies have used femtosecond pump–probe and two-dimensional infrared (2DIR) spectroscopies to gain molecular information of the dynamical processes occurring in these systems.^26–31 Moreover, these experiments have focused on either the solvent or the anion to investigate the solvation structure of the lithium ion since the cation itself does not have a signature in the infrared region. Using different infrared signatures, previous studies shed light on the ultrafast dynamics of the lithium ion solvation shell deformation, solvent exchange process, preferential solvation, and ion pair formation.^{26–28,31–34} In addition, the studies also revealed the existence of the excitonic vibrational states formed by the carbonyl stretches of the solvent molecules solvating the lithium ion.^27,31,35 While the excitonic nature of some of the vibrational transitions has allowed us to observe processes that cannot be seen in solvents lacking strong vibrational coupling, it has also made the interpretation of the results difficult in terms of individual (localized) vibrational states.^31,36,37

Highly coupled vibrational transitions are not limited to the solvent molecules in the solvation shell of a cation. Excitonic vibrational transitions have been observed in various systems, including, but not limited to, ions,^38–40 molecules,^40–44 solvents,^45–52 peptides,^53–57 proteins,^58,59 and polymers.^60–62 Two-dimensional IR spectroscopy has been proven to be an excellent tool for characterizing the structure and tracing the dynamics in a system with coupled vibrational transitions.^27,54,59 For example, using the 2DIR anisotropy and quantum mechanical calculations, the angle among the amide I modes in trialanine was obtained.^53,55 Note that a similar strategy was used to derive the structure of other biomolecular systems.⁵⁷ However, the existence of excitonic transitions leads to difficulties in determining the different components of the vibrational Hamiltonian. In this case, the 2DIR spectra do not provide a direct picture of the underlying molecular processes, and hence, it is necessary to compute the vibrational frequency fluctuations and coupling constants for connecting the spectra with the details of the molecular structure and dynamics. A particular approach commonly used to investigate the local site when participating in excitonic states relies on shifting the vibrational transition of the modes of interest. The most widely used approaches to investigate local sites include shifting the vibrational mode by isotopic labeling,^{55,60,63–75} replacing the molecular system with a “similar” system,^76–78 and using vibrational transitions with small transition dipoles.⁷⁹ Aside from all major drawbacks associated with these approaches (such as cost and detection), the most problematic issue is that they all eliminate the possibility of obtaining structural information from the system since this information is encoded in the coupling constant.⁸⁰ Thus, if the objective is to obtain the geometrical and/or dynamical structural information of the molecular system, one must work with coupled vibrational transitions.

The vibrational Hamiltonian of the highly coupled transition is extremely complicated. For example, the Hamiltonian describing the vibrationally coupled carbonyl stretches of the organic carbonates in the first solvation shell of the lithium ion can be represented in the site basis as

where β^ij is the coupling constant between the i-th and j-th sites and δωⁱ is the frequency fluctuation of the i-th site.³¹ This vibrational Hamiltonian can be simplified using strict approximations, such as coupling constants being equal (β^ij = β) and frequency fluctuations being smaller than the coupling constant (β ≫ δω).³¹ The instantaneous picture of the vibrational Hamiltonian without simplifications can be obtained by its direct diagonalization in the site basis representation. However, this approach requires the knowledge of the instantaneous site frequency fluctuations (δωⁱ) and coupling constants (β^ij). To this end, the approach relies on molecular dynamics (MD) simulations and frequency maps to compute the evolution of vibrational frequencies and coupling constants as a function of time.^81,82 This approach has been successfully used to compute systems such as peptides and pure solvents, where the coupling between sites can be modeled using dipole–dipole coupling,⁸³ transition dipole coupling,^84–87 transition charge coupling,^86,88 and vibrational coupling maps.^{82,85,86,88–91} However, in systems with strong quantum effects, such as charge delocalization, it is likely that the proposed coupling theories will not work and a different approach involving the explicit computations of the vibrational potential is needed.^92,93 Under the limit of computational capacity and efficiency, it is not practical to apply the latter approach for many different systems.

In this paper, a new approach is presented to investigate the dynamics of highly coupled vibrational transitions. The new approach involves using neural networks to compute vibrational frequencies. Typically, vibrational transitions are computed in ab initio quantum chemical computations by computing the second derivatives of the potential energy surface (Hessian) at the minimum of energy.⁹⁴ The requirement of being in the minimum of the potential energy surface hinders the calculation of the vibrational modes in snapshots of the MD simulation since the system is clearly not in the minimum of the potential energy surface at any time of the MD simulation. Neural networks can overcome this problem since they can be trained to compute vibrational frequencies from geometrical coordinates of the system in similarity to the method used for computing potential energies at the density functional theory (DFT)-level.^95,96 Moreover, it has been previously shown that neural networks are capable of predicting the vibrational frequencies of a single transition without any physical model and mathematical formula.^50,97

In this work, neural networks are used to obtain the distribution and dynamics of the instantaneous frequencies corresponding to the highly coupled carbonyl stretches of organic carbonates in the first solvation shell of the lithium ion. In contrast to the previous methodologies, the study focuses on directly computing the vibrational frequencies of the excitonic states. This methodology avoids the inversion problem and allows a one to one comparison with the experiments. In addition, the flexibility of the neural network methodology permits us to evaluate the use of geometrical factors of the molecular system to compute the instantaneous frequencies.

The training datasets of hydrated N-methylacetamide (NMA) were obtained from the classical MD simulation. The MD simulation was performed using the SANDER module of the AMBER 14 program package.⁹⁸ NMA was modeled with the ff14SB protein force field⁹⁹ and water with TIP3P.¹⁰⁰ The simulation system consisted of a 15 Å box containing one NMA molecule and 1522 water molecules. Periodic boundary conditions were imposed in the simulation, while particle mesh Ewald methodology was used to describe the long range electrostatic interaction with the cutoff of 8 Å. A Langevin thermostat was used to maintain the temperature of the system. The SHAKE algorithm was also used to constrain the bonds involving hydrogen. A production run was performed in the NVT ensemble at 300 K for 1 ns with a time step of 2 fs after energy minimization and an isothermal isobaric ensemble run.

The radial distribution function [g(r)] between the oxygen atom of NMA and the oxygen atom of water showed the first solvation shell with a coordination number of 2.7 at the cutoff distance of 3.4 Å (Fig. S1 of the supplementary material). Note that this selection is based on a previous work that showed that the frequency of NMA is correctly represented by three water molecules and a polarizable continuum.¹⁰¹ From the NVT simulation, 100 frames were selected for producing the neural network training set. The frequency of the amide I mode for each cluster (NMA + 3 water solvation shells) was calculated using density functional theory (DFT) at the B3LYP functional level with a 6-31+G^* basis set.¹⁰² All DFT computations were performed with the Gaussian 09 package.¹⁰³ 

The neural network was trained using the Cartesian coordinates of the NMA first solvation shell as the input and the frequency of the amide I mode of NMA as the output. The training set, consisting of 100 NMA configurations and their first solvation shells, was randomly selected from 1000 snapshots extracted every 10 ps from the production run. To minimize the effects of translational and rotational motions on the coordinates, the different NMA configurations and their water clusters were aligned by minimizing the root-mean-square deviation (RMSD) of atomic positions with respect to the gas phase geometry of NMA. The training was performed by randomly picking 70 configurations as the training set, 15 configurations as the validation set, and 15 configurations as the test set from the original 100 configurations. The neural network with the best performance was selected from various tests on neural networks with different algorithms and different architectures using the MATLAB network toolbox.¹⁰⁴ 

The selected neural network was then used to compute the instantaneous vibrational frequency of the amide I mode from the NVT production run.

The evolution of the lithium-carbonate clusters was derived from an ab initio molecular dynamics (AIMD) simulation carried out with the CP2K/Quickstep package (version 3.0).¹⁰⁵ The details of the computational simulation were described previously in Ref. 81. In short, the AIMD was computed using the Perdew-Burke- Ernzerhof (PBE) functional with the D2¹⁰⁶ dispersion scheme and the TZ2VP basis set with Goedecker–Teter–Hutter (GTH) pseudopotentials.^107–109 The simulated box consisted of 1 lithium ion (Li⁺), 1 hexafluorophosphate ion (PF₆⁻), and either 19 butylene carbonate (BC) or 23 dimethyl carbonate (DMC) molecules. Periodic boundary conditions were applied in both systems. These two systems were equilibrated and simulated in the NVT ensemble for 200 ps with a constant temperature of 300 K using a Nosé–Hoover thermostat. The coordinates of the lithium ion solvation shell were extracted from these trajectories and aligned with the optimized geometry of the isolated clusters.

The training datasets used in the neural network consisted of clusters with a lithium ion and four carbonates as derived from the ab initio molecular dynamics simulations.⁸¹ To this end, two kinds of clusters were investigated. The first was composed of a lithium ion and four butylene carbonate (BC) molecules, while the other comprised four dimethyl carbonate (DMC) molecules surrounding the cation. The structures of these clusters were derived by energy minimization using DFT at the B3LYP functional level with a 6-311++G^** basis set.¹¹⁰ All DFT computations were performed with the Gaussian 09 package.¹⁰³ 

A diverse training set was obtained by preparing different solvation structures with specific perturbations in the molecular arrangement of the solvent molecules surrounding the lithium ion due to the lack of a sufficiently long ab initio MD simulation.⁸¹ In particular, the geometrical factors chosen for altering the solvation shell structure were the Li–O distance, the Li–O–C angle, and the O–Li–O angle (Scheme 1) since they were determined to be the dominant contributions to the frequencies of the carbonyl stretch modes.^27,31,81 The magnitude of the perturbation of the geometrical factor was determined from the thermal distribution of these geometrical elements obtained from their ab initio molecular dynamics simulation.⁸¹ Thus, 100 different geometrical configurations of the lithium ion solvation shells were produced by the combination of four Li–O distances, four Li–O–C angles, six O–Li–O angles, and the different magnitudes of perturbation.

These various solvation structures were first energy-minimized in constraint optimization, where only the selected geometrical factors varied (perturbation) were constrained, followed by a normal mode calculation using the same level of theory as previously described. Note that the constraint optimization was necessary to assure that the complex will stay with the changed geometrical factors during the optimization because the perturbations involved motions with small barriers (thermal displacements). The geometrical factors, such as the atomic position, were obtained by minimizing the RMSD with respect to the original cluster to discard the effects of translational and rotational motions on the coordinates. The coordinates of the atoms, as well as other geometrical factors of the cluster, and the corresponding frequencies were used as the neural network training set.

The neural networks were trained using either Cartesian coordinates of some of the atoms or geometrical factors as the input and the four frequencies of the coupled carbonyl stretch transitions (located around 1750–1850 cm⁻¹) as the output. In the case of Cartesian coordinates, only the atomic coordinates of the four carbonyl groups (i.e., carbon and oxygen atoms) and the lithium ion were used as input data for obtaining the neural network model. Note that the number of atomic coordinates used as input was selected from the neural network model that produced the highest correlation using the lowest number of atomic coordinates. The training was performed by randomly picking 70 configurations as the training set, 15 configurations as the validation set, and 15 configurations as the test set from the original 100 configurations. The neural network with the best performance was selected from various tests on neural networks with different algorithms and different architectures using the MATLAB neural network toolbox.¹⁰⁴ For each neural network algorithm and architecture, the neural network was trained at least 20 times. The selection of the best neural network was based exclusively on the correlation coefficient for the whole training set.

The neural network was trained using the Levenberg–Marquardt (LM) and Bayesian Regularization (BR) algorithms. The BR algorithm implemented here was found to be a good alternative for the creation of neural networks for small datasets.^111,112 However, the high efficiency of the BR algorithm is accompanied by the problem of overfitting because this algorithm does not employ a validation set. Other commonly used algorithms, such as the LM algorithm, prevent the problem of overfitting by employing validation sets.¹¹³ 

In the case of the BR neural network, the optimum neural network model contained 3 hidden layers, and the Pearson correlation coefficient (R value) of the training set was found to be 1.00, while the R value for the whole test set was 0.80 (see the supplementary material). The decrease in the R value indicates the existence of overfitting in the training set. Frequencies of the amide I mode of NMA predicted from the different snapshots of the MD trajectories presented a Gaussian distribution with its central frequency located at 1800 cm⁻¹, which is 8 cm⁻¹ higher than the gas phase frequency. The Gaussian-like distribution of frequencies computed by the neural network is in agreement with previous computations of the NMA line shape using frequency maps (Fig. 1). The frequency–frequency correlation function (FFCF) of the amide I vibrational transition of NMA shows a decay that is well modeled with a bi-exponential decay (Fig. 1). Fast dynamics with the decay time of ∼50 fs and relatively slower dynamics on the timescale of 1.1 ps are deduced from the fitting. The observed dynamics are in agreement with the previous reports of the FFCF computed from time resolved infrared spectroscopy and the ab initio molecular dynamics simulation.^114–116 However, there is a large discrepancy in the short timescale, which could be attributed to the incorrect sampling of the water librations when building the neural network model due to the small displacements of these motions.

The neural network trained with the LM algorithm was also computed (see the supplementary material). In this case, the R value after the training was found to be 0.88, which was lower than the one from the BR algorithm. The difference between the neural network models might arise from the overfitting problem of the BR algorithm.¹¹⁷ The LM neural network also produced a Gaussian distribution of frequencies (Fig. 1). Moreover, the FFCF of instantaneous frequencies computed from the LM neural network also showed a bi-exponential decay. The amplitudes of the frequency fluctuations obtained from the fitting (Table I) present a similar ratio for the two neural network algorithms, but the BR neural network predicts amplitudes of the fluctuation that are ∼7 times larger than the neural network using the LM algorithm. Although the LM neural network presented a narrower distribution in frequency, it also predicted that the frequency fluctuations have a decorrelation time of 1.0 ps on the same timescale as the BR algorithm. It is important to note that the LM algorithm is susceptible to find the local minimum rather than the global minimum.¹¹⁸ Hence, it is expected that the LM algorithm will be less accurate as compared to the BR algorithm.

Compared to the 2DIR experimental results, the decorrelation dynamics of the frequency fluctuations are predicted very accurately by either neural network model, but the amplitudes of the fluctuations computed from the neural network models are smaller than the experimental results.¹¹⁹ In addition, there is a difference in the central frequency of the distribution for both models, which arises from the selected functional and the gas phase conditions used to compute the frequencies (see Sec. II). While the difference could be attributed to either the method or the approximations, the computation of the frequency fluctuation using the same coordinates (NMA + 3 water molecules) and the maps developed by Skinner (SM) and Tokmakoff (TM) showed the distribution of frequencies with widths similar to those of the neural network, providing support to our new methodology (Fig. 2).^120,121 Moreover, the FFCF dynamics derived from both frequency maps appear (SM and TM) to be faster for the long time decay than those predicted using the neural network (Table I). In contrast, the short time FFCF dynamics predicted by the maps are very similar to those of the neural network models. Thus, it appears that the neural network methodology is correct, but the DFT frequencies used to build the neural network model are not sufficiently accurate to appropriately describe the NMA–water system, given that the neural network models predict substantially smaller amplitudes of the frequency fluctuations, while the FFCF correlation decays are similar to those predicted by the maps (Table I). The shortcomings of the neural network model are likely caused by the assumptions made here. First, the amide I frequency in the NMA molecule is only considered to be under the influence of the water molecules in the first solvation shell of the carbonyl group. This is clearly not true in the NMA system since the water outside the first solvation shell of the carbonyl group is shown to be important for the prediction of the NMA-amide I frequency.^101,122 Second, the DFT B3LYP functional and the 6-31+G^* basis set appear to be critical in building the neural network model as when computing the frequency fluctuations from DFT computations for frequency maps.¹²³ Thus, a more accurate model should take into consideration the water molecules interacting with the N–H group of the amide. In addition, a more accurate DFT functional/basis set pair and larger sampling for building the neural network model should increase the overall prediction performance of the neural network model. However, it is not expected that this problem will hinder the application of the method to the case of the solvated lithium ion because the carbonyl modes in the Li(solvent)₄⁺ cluster are more isolated. Their isolation arises from the positioning of these solvent molecules such that their carbonyl groups point toward Li⁺, while the solvent molecules in the second solvation shell have their carbonyl stretch pointing away from Li⁺. Hence, the arrangement of the first and second solvation shells minimizes their vibrational interaction.⁸¹ 

Overall, it has been demonstrated that the neural network has a reliable performance in predicting frequencies at a given DFT level and tracing the dynamics of the frequency fluctuation of a single vibrational mode, without any additional information or assumptions on the functional form of the electric field acting on the selected site of the solute molecules. The results presented here agree with the work by Cho’s group, in which they reproduced the frequency fluctuation of the NMA amide I mode using neural networks with different algorithms and different architectures.⁹⁷ 

Previous work established that the carbonyl stretches of the organic carbonates solvating Li⁺ are highly coupled vibrational transitions.⁸¹ To compute the frequencies of these highly coupled vibrational transitions, neural networks were trained using different algorithms and different numbers of layers. As in the case of NMA, the performance of the neural network was assessed using the correlation coefficient (R). All the algorithms produced an optimum model with 8–11 layers for both linear and cyclic carbonates (see the supplementary material).

It is observed that the model produces an R value higher than 0.99 in all the neural network architectures irrespective of algorithms, demonstrating the good performance of the method as compared to common DFT frequency maps.⁵⁸ In addition, the BR-based neural network has the best performance, which agrees with this algorithm being designed for noisy, small, and difficult datasets.^112,124,125 However, the BR algorithm often produces over-fitted models, which do not match the results of other modeling methodologies, such as frequency maps, as shown in Sec. III A. Therefore, the neural networks with the LM algorithm (the second best performer according to the R value) were adopted for comparison. Since all the different neural network models perform exceedingly well with 11 layers (see the supplementary material), the instantaneous frequencies were derived using these number of layers and different algorithms.

The frequency distributions of four coupled carbonyl stretching modes in both [Li(BC)₄]⁺ and [Li(DMC)₄]⁺ for the neural networks using the BR and LM algorithms are presented in Fig. 3. The results show the presence of four coupled vibrational transitions, in which three of the vibrational transitions are downshifted from the fourth one for both BC and DMC irrespective of the neural network algorithms. In addition, all the instantaneous frequencies present Gaussian-like distributions and show substantial overlap among the three frequency distributions with the lower frequencies. As in the case of the NMA frequencies, the LM algorithm produces narrower frequency distributions for both the DMC and BC systems (Fig. 3), but the effect is more noticeable for the BC sample. The difference between the two algorithms might be associated with the drawbacks of the LM algorithm as previously described. Overall, the results of the neural network modeling agree with the results of IR spectroscopy and quantum calculations, in which only a pair of peaks was observed due to the spectral overlap of three of the transitions.²⁷ The low frequency transitions correspond to the infrared-allowed asymmetric stretches of the four coupled transitions, while the high frequency transition is infrared forbidden and overlaps with the band of the free solvent.²⁷ 

The average coupling constant was also calculated from the difference between the low frequencies and the high frequency. Using a vibrational Hamiltonian in which the coupling constants among sites are the same and the fluctuations of the site frequencies are small,²⁷ the excitonic frequencies can be approximated to

and

where ω₁₀ is the site frequency, β is the coupling constant, and δω_i is the frequency fluctuation of the i-th site. Hence, within the framework of this model, one can compute the average coupling constant from

because by definition, $δωi=0$⁠. The computed coupling constant is found to be 13.2 and 9.4 cm⁻¹ for BC and DMC systems, respectively. These coupling constants are in the same order of magnitude as the experimental values of 6.0 and 7.9 cm⁻¹ for BC and DMC systems, respectively.²⁷ The observed error of the computed coupling constant is likely to arise from the selected DFT functional and basis set, which might not be the most accurate for describing [Li(solvent)₄]⁺

The FFCF of the four coupled transitions calculated for both BC and DMC samples are shown in Fig. 4. Note that the FFCF of low frequency is averaged over three overlapped transitions. The plot showcases that the FFCFs extracted from different transitions in the same system have similar decay time. It is also evident that the FFCFs of the transitions for [Li(BC)₄]⁺ have slower dynamics than those of [Li(DMC)₄]⁺. The characteristic times of FFCF dynamics, obtained by fitting with exponential decays, are shown in Table II. The low frequency transitions present an average decay time that is three times larger for BC than DMC irrespective of the neural network algorithm. Previous experimental results showed decay times for the FFCF of 12.8 and 1.4 ps at x(Li⁺) = 0.09 and 4.3 and 2.3 ps at X(Li⁺) = 0.15 for BC and DMC systems, respectively.²⁷ Hence, when compared to the experimental results, the neural network modeling presents the same trend, where the FFCF of coupled transitions for the coordinated carbonates has slower dynamics for BC than DMC. The match between the predicted and experimental FFCFs demonstrates that the neural network can correctly compute the frequencies of the system in highly coupled systems, such as organic carbonates solvating Li⁺. The difference between the neural network and the experimental time constants is likely to arise from the insufficient sampling of ab initio modeling of the system using a small box and a short production run due to its high computational cost.

The modeling of the FFCF of highly coupled vibrational transitions shows a good agreement with the experimental result. Hence, it was also studied whether geometrical factors can be used to compute the FFCF of the lithium ion coordination shell. Previously, it was shown that the dynamics of the first solvation shell of the lithium ion in BC and in DMC were well described using the time evolution of the geometry factor q, where q describes how tetrahedral the system is.⁸¹ Consequently, a new neural network model that computes the frequencies of the four coupled transitions for the lithium-carbonate complex as a function of both the geometry factor of coordinated carbonyl carbon atoms (qC) and the geometry factor of coordinated carbonyl oxygen atoms (qO) was also developed. In this case, the neural networks were trained with only two inputs (qC and qO) and four output target frequencies. In other words, the detailed information of the atomistic arrangement of the lithium-carbonate clusters was discarded, and only two metrics describing the arrangement of oxygen and carbon atoms around the lithium ion were kept as input. Similar to the training of the neural network using the atomic position, all the new neural network models produced R values that are higher than 0.9, indicating that the simplified model still shows good performance in predicting the frequencies of coupled carbonyl transitions (see the supplementary material). In particular, it is observed that the optimal neural network of the geometrical factor q has an 11-layer architecture and was trained with the LM algorithm, which shows similar performance to the previous network with 27 coordinates as input.

The FFCFs calculated using the geometrical factor neural network again show the same trend (Fig. 5) as previously described, i.e., faster dynamics for DMC than BC (Table II). The agreement in the FFCF dynamics indicates that this simplified model also does a decent job in predicting the frequency fluctuations in the two Li⁺-carbonate systems. While the trend of the FFCF is correctly predicted, the dynamics are not as accurate as the ones predicted with the atomic coordinates of the carbonyl oxygen and carbon atoms. This result is not surprising, given the extreme simplification of the system from 27 (atomic coordinates) to 2 variables (geometrical factors).

In this paper, the use of artificial neural networks as a new approach to the frequency computation of highly coupled transitions is showcased. First, the methodology was validated by comparing the frequency fluctuation of the amide I mode of NMA using neural networks, frequency maps, and experimental data. Second, the neural network methodology was used to compute the frequency dynamics of the highly coupled carbonyl stretch modes for organic carbonates solvating the lithium ion. The high R values of the neural network model, in addition to the good agreement of the FFCF derived from this methodology with the experimental results, demonstrate the good performance of the neural network in the frequency modeling of highly coupled (excitonic) vibrational modes. In addition, the modeling of the system with different geometrical parameters, such as atomic coordinates and tetrahedral geometrical factors, was explored. In both cases, the modeling of the FFCF predicted the correct trend, with the experiments further supporting the suitability of the methodology. While this new approach requires neural network training for each specific system, it can shed light on the frequency dynamics of systems with highly coupled transitions, which otherwise require a very tedious inversion of their complicated vibrational Hamiltonians. Overall, the presented methodology provided a new alternative for studying highly coupled vibrational transitions, such as those observed in polymers, proteins, and solvents.

See the supplementary material for the radial distribution function between the oxygen atom in NMA and the oxygen in water; the comparison of frequency distributions and FFCFs of the amide I mode using neural networks, frequency maps, and the DFT calculation results; the plots showing the performance of the optimum neural networks trained for the hydrated NMA system; the optimal R values for the different neural networks trained with different architectures and algorithms; and the representatives of the neural networks trained for BC or DMC solvating Li⁺ using the BR or LM algorithm.

The data that support the findings of this study are available from the corresponding author upon reasonable request.

The authors would like to acknowledge financial support from the National Science Foundation (Grant No. CHE-1751735). The authors would also like to acknowledge the computer time provided by the High Performance Computing Center at Louisiana State University and the Louisiana Optical Network Initiative (LONI).