h-CMD: An efficient hybrid fast centroid and quasi-centroid molecular dynamics method for the simulation of vibrational spectra Open Access

Liquid–air, liquid–liquid, and liquid–solid interfaces play critical roles in various chemical, biological, materials, and environmental processes. Vibrational infrared (IR) spectroscopy provides a wealth of information about interactions of molecules in condensed phases and at different interfaces, especially in aqueous environments.^1–4 In particular, IR spectra can reveal key fundamental insights into the complex hydrogen bond (HB) network of water and how it is influenced by temperature, pH, reducing or oxidizing potentials, and the presence of different ions, among others. Reliable atomistic molecular dynamics (MD) simulation of vibrational spectra requires an efficient way of incorporating nuclear quantum effects (NQEs). However, while the role of NQEs in chemical and material sciences is a widely acknowledged fact, the development of efficient computational methods and software platforms capable of incorporating them into MD simulations of complex systems has been a challenging task.⁵ Path-integral approaches, including path-integral MD (PIMD),⁶ centroid MD (CMD),^7–10 and ring polymer MD (RPMD),^11–13 as well as a myriad of more sophisticated methods branched from them,¹⁴ have emerged as the leading platform for the incorporation of NQEs in large-scale condensed-phase simulations. In these methods, NQEs are captured using Feynman’s imaginary-time path-integral formalism^15,16 by taking advantage of the so-called isomorphism between the quantum-mechanical partition function and the partition function of a classical ring polymer (RP). Every nuclear degree of freedom (DOF) is represented by an RP comprised of n copies of the system, known as “beads,” connected by harmonic springs. The extended phase space of this RP captures NQEs, such as zero-point energy and nuclear tunneling.¹⁷ Hence, accurate yet computationally efficient path integral (PI) techniques provide powerful alternatives to prohibitively expensive full quantum-mechanical dynamical simulations.

However, when it comes to calculating the IR spectra of complex heterogeneous condensed phases and interfaces, the quantum dynamic community has struggled to propose a computationally efficient PI method that can accurately capture NQEs,^18,19 which are more pronounced in spectra than other structural and dynamical properties. This is particularly notable in combination and overtone bands,^20,21 which are beyond the scope of this work. One of the most well-known formalisms to approximate the real-time dynamics while incorporating NQEs in the PI framework is CMD. In CMD, the dynamics of the particles are performed under the effective mean-field potential of the imaginary-time path-integral whose centroid is constrained to be at the position of the particle. For low dimensional systems, this effective potential can be calculated on a grid prior to dynamical simulations, while for large systems, it must be calculated “on-the-fly,” which is referred to as adiabatic CMD (ACMD).⁹ This is achieved by introducing an adiabatic separation between the centroid and internal modes by using the physical mass of the particle as the centroid mass, scaling down the mass of the internal modes of the RP, and attaching them to a thermostat so that they sample the equilibrium distribution while being constrained to the position of the slower moving centroid.⁹ Full adiabatic separation in CMD requires a very small mass for the internal modes and, thus, a very small simulation time step for the dynamics to be accurate. The computational cost associated with such choices led to the development of the partially adiabatic CMD (PA-CMD) method, where the mass scaling is not as extreme, allowing for a larger time step while still providing accurate dynamics.²² Further work in reducing the computational cost of CMD comes in the form of fast CMD (f-CMD), where the effective potential is learned beforehand from PIMD simulations using force matching to create a new classical analytical force field²³ or neural network potential (NNP).²⁴ 

It has been documented that vibrational spectra calculated using CMD methods exhibit redshifting of peaks, primarily of stretch bands, which is worsened with the lowering of temperatures and/or for lighter particles.²⁵ This phenomenon, dubbed the “curvature problem,” is a result of the RP spreading out along the angular coordinates of a system when approaching the inside of the stretching coordinates. This stretching reduces the force along the stretching coordinate, lowering the oscillation frequency of that coordinate and redshifting the spectra. One approach to overcome the curvature problem is the quasi-centroid MD (QCMD) method introduced by Althorpe and co-workers.^26–28 QCMD works by defining a quasi-centroid (QC) for the molecule using curvilinear coordinates and determining the effective mean-field potential using an imaginary time RP that has its QC constrained to the particle position as opposed to its Cartesian centroid.²⁶ The resulting potential does not exhibit the same flattening on the inside of stretching curves as in CMD and thus does not have the redshifting of vibrational spectra. The original adiabatic implementation of QCMD does not lend itself to scaling to complex chemical systems due to computational complexity and cost.^26,28 Inspired by the f-CMD approach, Manolopoulos and co-workers developed the f-QCMD method,²⁹ which has been used to study vibrational spectra of gas phase molecules,²⁹ condensed phase systems,^29,30 and cavity confined molecules.³¹ In f-QCMD, a corrective potential to the classical force field is learned to create QCMD-like dynamics in classical simulations, but the corrections are found using iterative Boltzmann inversion (IBI)^32,33 on a set of distribution functions instead of force matching. While f-QCMD has shown reduced computational cost compared to QCMD,^29,30 determining the distribution functions for large molecules/frameworks poses a challenge. In a similar vein, Kapil and co-workers developed the temperature elevation PI coarse-graining simulations (T_e PIGS) method.³⁴ Similar to some f-CMD implementations, a classical NNP is learned through force matching a PIMD simulation to incorporate NQEs. However, the PIMD simulations are performed at a high temperature where the RP distribution is more compact.³⁴ The resulting NNP is then used for calculating vibrational spectra at lower temperatures without the curvature problem of CMD^34,35 and could provide an efficient method for use in systems with complex molecules. An additional method for incorporating NQEs at a cost very close to classical simulations is the adaptive quantum thermal bath (adQTB) method,^21,36 where vibrational spectra of bulk water do not show the curvature problem of CMD at room temperature.

Here, to circumvent the curvature problem and create an accurate yet efficient scheme for calculating vibrational spectra of complex condensed phases and interfaces, we have devised a new method coined hybrid CMD (h-CMD) that combines f-QCMD and f-CMD methods. This choice is motivated by the demonstrated applicability of the QCMD method to a number of chemical systems and allows us to simulate interfacial systems involving complex frameworks with high accuracy for the most important degrees of freedom and reduced complexity for the others. To bring f-QCMD and f-CMD methods on equal computational footing, the latter f-CMD method is implemented with the regularized IBI method used in f-QCMD. We showcase the accuracy of our new hybrid h-CMD method on simulating IR spectra of D₂O confined in the zeolitic-imidazolate framework-90 (ZIF-90) compared to the experiment as reference. By comparing the spectra using our hybrid h-CMD method with PA-CMD and T-RPMD, as well as f-CMD and classical MD at different temperatures, we provide detailed insights into the applicability of this new scheme to complex interfaces.

This work is structured as follows: in Sec. II, we outline our h-CMD method, followed by the methods used for calculating vibrational spectra. Simulation details are provided in Sec. III. In Sec. IV, we provide our detailed analyses of the IR spectra of interfacial D₂O confined in ZIF-90 at different temperatures, followed by conclusions and future studies.

The CMD framework, particularly PA-CMD, stands as a popular method for incorporating NQEs in calculating vibrational spectra despite its known issue of the curvature problem. QCMD has been shown to remedy the curvature problem for some simple chemical systems,^26–28 but its computational complexity hinders its applicability for general use. While f-QCMD does a lot of work to simplify the QCMD method,^29,30 it has its own limitations for its use for complex systems. Here, we introduce a new method that we refer to as hybrid CMD (h-CMD). It combines the f-QCMD and f-CMD methods to create a scheme that is applicable to heterogeneous systems containing large complex molecules and/or extended frameworks. In this section, we give an overview of the QCMD and f-QCMD methods and discuss their limitations for complex molecules and extended frameworks. We then give our rationale for the h-CMD method and discuss its formulation.

For gas phase systems, the forces on the QCs can be found directly through their relationship with the curvilinear coordinates.²⁶ The original QCMD simulations used an adiabatic implementation (AQCMD) where two systems were evolved simultaneously, one for the particles of interest and another for the corresponding RPs, which had their QCs constrained to the positions of the particles.^26–28 In condensed phase systems, multiple molecules interact with each other, so the positions of the QCs of each molecule relative to each other must be known. While the set of curvilinear coordinates is enough to define the positions of the atoms in a molecule with respect to each other, it cannot uniquely define the molecule’s location and orientation in 3D space. QCMD uses an Eckart-like frame^37,38 to orient the molecules by rotating the QCs into the frame of the Cartesian centroid.²⁶ These Eckart-like conditions are met through the use of SHAKE and RATTLE algorithms,^39–41 while the forces between molecules are introduced through rotational forces.²⁶ Approximations to the torque on the QCs are needed,²⁶ with work to improve the estimation of the torque being done.²⁸ 

The AQCMD method suffers from a high computational expense associated with its need to use a small mass and, therefore, a small time step to achieve the adiabatic separation for properly sampling the effective mean-field potential. In an effort to make QCMD more feasible for large-scale simulations, the f-QCMD method has been developed.^29,30 f-QCMD is inspired by the f-CMD method, where the effective potential is learned beforehand, but the potential is meant to mimic that of QCMD. This allows for classical-like simulations with the same dynamics as those of AQCMD at a much lower cost.

Using this implementation of f-CMD for complex molecules reduces the complexity of calculating the reference distributions, as no rotation matrices will be needed for those molecules. As mentioned above, for heterogeneous systems, treating small molecules where the accuracy of the incorporation of NQEs is important to the spectra with f-QCMD and larger molecules and frameworks with f-CMD is the basis of our h-CMD method. It should be emphasized that the difference in the treatment of the molecules is in how the reference distributions are calculated; all the molecules are treated the same in the actual PIMD simulations, the IBI process, and the final classical-like dynamics under the correction potentials.

In order to test the accuracy and applicability of our new h-CMD scheme, we use it to simulate IR spectra of interfacial D₂O molecules confined in the ZIF-90 framework. We specifically selected this system as experimental IR spectra were available for comparison. The natural division of this system in the h-CMD framework is to treat the D₂O molecules with f-QCMD and ZIF-90 with f-CMD, as schematically shown in Fig. 1.

We would like to emphasize that this method allows for improved accuracy in terms of the removal of the curvature problem for the relevant molecules while being simple to implement for large, more complex, and heterogeneous condensed phases and interfaces. The h-CMD method will suffer in accuracy for systems where spectral features of interest are part of the large molecule/framework that is being treated with f-CMD, as that feature will still exhibit the curvature problem. We still believe that the method has wide use for complex condensed phase and interfacial systems where the spectra of small molecule adsorbates are the focus.

All classical and quantum PI simulations are performed using a development version of our software package DL_POLY Quantum 2.0⁴⁷ using the experimentally obtained 1 × 1 × 1 crystal structure of ZIF-90,⁴⁸ utilizing periodic boundary conditions. DL_POLY Quantum 2.0 is a highly parallelized, modular, and user-friendly computational platform for classical and advanced PI simulations in condensed phases. The flexible anharmonic q-TIP4P/F⁴³ quantum water potential was used for modeling water, which is shown to be very successful in reproducing different properties of water, including its phase diagram,⁴⁹ melting point,^50–52 and temperature of maximum density^43,51 as well as the magnitude of the isotope fractionations⁵³ and vibrational infrared spectra.⁴³ The equations of motion were propagated using the velocity–Verlet algorithm with a time step of 0.2 fs, with the temperature being kept constant at 300 K. The path integral Langevin equation (PILE) thermostat⁵⁴ with a relaxation time of 0.05 ps and barostat with a relaxation time of 2.0 ps were employed as implemented in DL_POLY Quantum 2.0. All bonded and non-bonded parameters of the ZIF-90 framework were obtained from Ref. 55 (see the supplementary material, Table S1). Water–framework interactions were included using the non-bonded electrostatic and Lennard-Jones potentials. Lorentz–Berthelot mixing rules were used to drive cross-interaction terms with parameters for water taken from the q-TIP4P/F water model, which is the same as the original TIP4P/2005.⁵⁶ To account for properly capturing the interactions between the water molecules and ZIF-90 framework, we modify the non-bonded Lennard-Jones (LJ) parameters between the carbonyl oxygen (O) of the ZIF-90 framework and the oxygen atoms of D₂O (O_W) adsorbate molecules by increasing $εO−OW$ by 50% and reducing $σO−OW$ by 10%. This adjustment decreases the O–O_W repulsion, allowing D₂O molecules to approach the framework more closely and form stronger hydrogen bonds with the oxygen of the carbonyl groups, which is essential for reproducing the correct IR line shapes observed in the experimental spectra. A similar modification to this LJ cross term was made in Ref. 55 using the MB-pol water model.⁵⁷ However, in contrast to Ref. 55, where the $εO−OW$ parameter was decreased by 20%, we increased $εO−OW$ by 50% which was based on fine-tuning the interaction to match the spectral line shapes and relative intensities more closely with experimental observations as mentioned above (see the supplementary material, Sec. S1 and Fig. S1, for more details).

For this study, we focused on the highest relative humidity (RH) of 60%, which is also considered to be the saturation limit of this material at room temperature. The initial packing of D₂O molecules within the ZIF-90 framework was performed using PACKMOL,⁵⁸ ensuring random distribution inside the cavities based on geometrical constraints. To achieve a fully randomized initial configuration of the D₂O molecules within the ZIF-90 framework, simulated annealing was performed using three short consecutive canonical NVT simulations at varying temperatures. Initially, the system was simulated at 1000 K for 50 ps to facilitate rapid diffusion and avoid local clustering. Subsequently, the temperature was lowered to 500 K for 30 ps, and finally, the system was equilibrated at the target temperature of 300 K for 20 ps. This step-wise reduction in temperature was designed to ensure uniform distribution of the confined D₂O molecules. Atom–atom cutoff distances were set to 8.5 Å to truncate the short-range interactions, whereas the long-range electrostatics were calculated using the Ewald summation technique.⁵⁹ Following simulated annealing, the system was equilibrated in the isobaric–isothermal (NPT) ensemble at 300 K and 1 atm for 500 ps, allowing the volume to relax (see the supplementary material, Table S2, for the MD and PIMD equilibrated cell lengths compared to the experiment). The system was further equilibrated using 32-bead PIMD simulations in the NPT ensemble for 200 ps. The NPT equilibrated box length of 17.352 Å was used for all subsequent NVT and NVE simulations. PIMD equilibrations were then performed in the NVT ensemble for 50 ps using a PILE thermostat, followed by an additional 50 ps production simulation for generating sampling initial configurations for subsequent NVE simulations.

For h-CMD simulations, we partition the system so that the adsorbed water molecules are treated with the f-QCMD method while the ZIF-90 framework atoms are treated with f-CMD as schematically shown in Fig. 1. In addition, we further approximate our treatment of the framework atoms by setting all correction potentials for interactions involving just the framework atoms to zero. This reduces our intra-molecular corrections to be just for the water O–D bond and D-O-D angle and our inter-molecular corrections to be the water–water and water–ZIF interactions. We justify this treatment on several bases. First, we are focused on the O–D stretch region in our spectra calculation, so it is more important to correct for the curvature problem in the water molecules as opposed to the framework. Second, the framework is composed primarily of heavy atoms in which NQEs are negligible at room temperature, meaning the correction potentials within the framework should be small enough to be ignored. Finally, while there are hydrogen atoms within the framework that do exhibit NQEs, they do not strongly interact with the water molecules as the water molecules donate hydrogen bonds to the framework rather than accepting them.

This treatment results in a total of 19 pair interactions, for which we calculate the correction potentials. To further test the importance of the different types of pair interactions, we also consider the case where we only correct the water-water pairs, leaving us with only 3 RDFs that are part of the inter-molecular corrections in the IBI process. To generate the distribution functions for the IBI process, we averaged over 50 independent PIMD simulations with 32 beads and an integration time step of 0.2 fs. The initial configurations of these runs were obtained from snapshots from the NVT production simulation described at the end of Sec. III A. The PILE thermostat with a time constant of 0.05 ps is used to achieve the NVT ensemble at 300 K. Each trajectory is first equilibrated for 50 ps, followed by 150 ps of production with the distribution functions calculated every 100 time steps using histogram binning.

During the IBI process, the distribution functions are generated from classical-like simulations under the correction potentials. The regularization parameter, ɛ, is set to 10, and a total of 30 IBI iterations are performed to converge the correction potentials; see the supplementary material Fig. S8, for the RMSE convergences with respect to the number of IBI iterations. For each iteration, the distribution functions were averaged over 50 independent NVT simulations. Each trajectory is evolved with a time step of 0.2 fs and thermostatted with a massive Nosé–Hoover chain (m-NHC) thermostat^60,61 of length 3 and a time constant of 0.05 ps. Trajectories were equilibrated for 50 ps and an additional 50 ps for data collection, where the distribution functions were calculated every 100 time steps, with the RDFs calculated using force sampling. The correction potential updates for each iteration were calculated using a modified version of the IBI code provided by the authors of Ref. 30.

For comparison, we also computed the vibrational spectra for the entire system using f-CMD. As with h-CMD, we only correct the water intra-molecular potential and the water–water and water-framework inter-molecular potentials. The simulation details are the same as the h-CMD discussed above but with the Cartesian centroid used in all distribution function calculations.

The initial configurations for the real-time dynamics trajectories for the spectra calculations were taken from the final set of PIMD trajectories discussed in Sec. III A, with configurations saved every 1 ps, providing a total of 50 independent initial configurations for subsequent dynamics simulation using PA-CMD and T-RPMD methods. Spectra calculations for all methods follow the same basic procedure. Configurations for dynamics trajectories are generated from NVT simulations at 300 K that involve 50 ps of equilibration and then 50 ps of sampling, with configurations taken every 1 ps. Sampling for h-CMD and f-CMD are all done using the m-NHC thermostat with a chain length of 3 and a time constant of 0.05 ps. Sampling for other PI methods is done with PIMD simulations using 32 beads and the PILE thermostat with a time constant of 0.05 ps. Finally, sampling for the classical MD simulation is done similarly to the PI methods, except that a single-bead RP is used, effectively treating the system classically.

Real-time dynamics for all methods are performed with 50 independent NVE trajectories for 50 ps. PA-CMD simulations use 32 beads, and the internal RP modes are thermostatted at 300 K with the m-NHC thermostats of length 3. T-RPMD simulations are done with 32 beads, and the internal modes are thermostatted with the PILE thermostat. IR spectra for D₂O are calculated using the total dipole derivative correlation function of the water molecules in the system. The first 5 ps of each trajectory is dropped, with the second 45 ps used, to ensure the system is well equilibrated, especially for methods that include thermostatting, and to avoid any artifacts associated with the initial transient dynamics. A Hann window with τ = 1 ps is applied to the correlation function.

We also consider a number of temperatures other than 300 K. Starting configurations for sampling at other temperatures come from the NPT equilibrated structure at 300 K, as described in Sec. III A. For h-CMD and f-CMD simulations, the IBI process of Sec. III B is performed at each temperature to obtain the necessary correction potentials. All other sampling and spectra calculations are the same as described above, with the exception of 64 beads being used in all PI simulations at 250 K.

As mentioned, IR spectroscopy is an excellent tool for probing the complex HB network of water in different condensed phases and interface environments. Experimentally, the deep inelastic neutron scattering (DINS) technique can provide direct information on proton momentum distribution and has been applied to both liquid bulk water⁶² and water confined in hydrophobic Nafion⁶³ with average pore diameters of ∼20 Å, silica nanopores⁶⁴ with average pore diameters of ∼24 and ∼82 Å as well as to carbon nanotubes with ∼16 Å diameters.⁶⁵ The data from DINS on water confined in Nafion showed that the nanoconfined water is in a distinctly different quantum state than that of bulk water.⁶³ From the theory standpoint, despite recent theoretical studies attempting to provide molecular-level insights into water structures at different interfaces, the impact of NQEs on the vibrational spectral signatures of the interfacial confined water still remains an open question. Here, we use the complex system of interfacial D₂O molecules confined in the ZIF-90 framework to test the accuracy of our new h-CMD method for infrared spectral simulations of large complex systems compared to a myriad of classical and quantum PI simulation methods, as implemented in our in-house developed software DL_POLY Quantum 2.0. We then discuss the ability of each approach to predict the infrared spectra at different temperatures of 250–600 K.

Considering the lower frequency HB peak from the experiment, which appears at 2535 cm⁻¹, our classical MD simulated spectrum shows a blue shift of 56 cm⁻¹. T-RPMD reduces this blue shift, positioning the first peak at the experimental value, but fails to accurately resolve the doublet feature due to the inherent artificial broadening of spectra caused by excessive friction in the dynamics. On the other hand, the PA-CMD calculated spectrum exhibits a redshift of 12 cm⁻¹ for the first peak compared to the experiment stemming from its well-known curvature problem. In contrast, our newly introduced h-CMD method successfully reproduces the doublet feature demonstrating this method’s ability to overcome the curvature problem with its O–D peak positioning in line with T-RPMD. Meanwhile, f-CMD behaves rather similarly to PA-CMD with slightly more redshifted O–D peaks, a behavior that was also observed in bulk water spectra for the force-matching f-CMD method.⁶⁷ 

The second higher-frequency experimental peak, which corresponds to the interfacial water molecules interacting with the ZIF-90 framework, appears at 2655 cm⁻¹. Classical MD once again demonstrates a notable blue shift of 49 cm⁻¹, reflective of its limitations in capturing the subtleties of HBs in confined water due to not including NQEs. As mentioned, with smeared peak structures, the fine details of the spectrum are lost in T-RPMD simulated spectra while PA-CMD redshifts the second peak by 36 cm⁻¹. On the other hand, h-CMD accurately captures both peak positions, with minimal shifts from the experimental data, and successfully resolves the doublet feature. The close agreement in both the positions and the relative intensities of the peaks between h-CMD and the experiment highlights its accuracy for calculating vibrational spectra while incorporating NQEs.

In our simulations, we have employed the q-TIP4P/F water model, while a similar study reported in Ref. 55 utilized the MB-pol model. MB-pol is a highly accurate water model constructed using a many-body expansion approach that is fitted to accurate quantum mechanical calculations.⁵⁷ MB-pol excels at describing both water clusters and bulk water, making it ideal for high-fidelity simulations. Unlike MB-pol, q-TIP4P/F is parameterized based on CMD simulations of bulk water properties, which limits its accuracy compared to the MB-pol water model. However, q-TIP4P/F has demonstrated its ability to reproduce key features of bulk D₂O spectra and to capture the line shapes of confined D₂O in ZIF-90, as shown in the simulated spectra presented in this work. Detailed comparisons between the simulated spectra obtained from the two water models for the liquid bulk and confined water systems are provided in the supplementary material, Fig. S3. Overall, a good agreement exists between the two water models, with only minor differences in the simulated IR spectra of both the liquid bulk and confined water systems.

As mentioned in Sec. III B, in h-CMD, intra-molecular corrections are considered only for the water O–D and D–O–D bonds and angles, with inter-molecular corrections included for all the water–water and water–framework interactions. Here, we consider the case of h-CMD with further simplification to the inter-molecular correction potentials. In this approximate h-CMD method, which we have labeled as h-CMD-D₂O, we only include correction potentials for the 3 water–water pairs. The resulting IR spectra for h-CMD-D₂O is shown in Fig. 3, where it is compared to the case with all water–water and water–framework interactions corrected as well as the experimental spectrum as reference.⁵⁵ 

As expected, the full h-CMD method, which considers all inter-molecular interactions, results in an IR spectrum that closely resembles that of the experiment. Interestingly, the h-CMD-D₂O method does provide a spectrum that is very close to the full h-CMD method, with only minor redshifting of the higher frequency peak. This is as expected since this peak corresponds to the water molecules donating HBs to the carbonyl oxygen of the ZIF-90 framework, and this interaction has not been corrected. However, this result does indicate that for a given system, the number of correction potentials needed can be optimized to reduce the computational cost without sacrificing too much accuracy.

Currently, experimental IR spectra for the temperatures investigated in this study, apart from 300 K, are unavailable. In the absence of experimental data for validation, our results serve as predictions for the temperature dependence of IR spectra beyond 300 K. As discussed above, comparison with spectra at 300 K gives us confidence that our methodology is able to provide qualitative predictions for the spectra at other temperatures. However, since our fine-tuning for the force field parameters is based upon reproducing the spectrum at 300 K, this modified potential may not accurately capture the non-bonded interactions at other temperatures, meaning that we cannot determine if the shifts in the spectra are quantitatively accurate. Nonetheless, we believe this approach is useful for observing the general trends in spectra for varying temperatures.

Figure 5 compares the temperature dependence of the IR spectra for PA-CMD and T-RPMD methods compared to that of h-CMD and h-CMD-D₂O. For the highest studied temperature of 600 K, except for h-CMD-D₂O, the simulated spectra calculated using all methods match very closely. All methods show a shifting of the spectra to lower frequencies with lowering temperatures. This shifting is the most dramatic with PA-CMD, highlighting how the curvature problem worsens with lowering the temperature, even for the heavier isotope of water. The broadening of the features of T-RPMD spectra is also demonstrated, with no doublet feature appearing at every temperature. The last comparison of notes is between h-CMD and h-CMD-D₂O. While the relative peak intensities do vary between the two treatments, the peak positions are very close to each other, except in the 600 K case, further showcasing the tunability of our scheme in reducing the number of correction potentials.

In this work, we have introduced a new hybrid f-CMD/f-QCMD scheme, where the accuracy of the more complex f-QCMD method is reserved for molecules that are believed to suffer from the curvature problem more severely, while the rest of the system is treated with the comparatively simple f-CMD method where determining the quasi-centroid coordinates is difficult or is entirely unfeasible. By simulating IR spectra of D₂O molecules confined in the extended ZIF-90 framework as a test case, the efficiency and accuracy of this scheme are showcased for large-scale PI simulations of complex condensed phase and interfacial systems. Our detailed analyses show that the newly introduced h-CMD method is capable of capturing the twin lower and higher frequency bulk-like and interfacial peaks for the nanoconfined water in ZIF-90 compared to the experiment. This is in sharp contrast to the traditionally employed T-RPMD method, where the artificial broadening of the IR spectra leads to the loss of the doublet feature, especially as the temperature is lowered.

We note that our h-CMD method is not the only avenue for efficiently simulating accurate vibrational spectra. The previously mentioned T_e-PIGS method could have success for these types of systems. Considering the largest difference between the two methods is in determining the effective potential used for the final simulations, it is difficult to directly compare the computational expense of the two methods without performing careful benchmarks. The elevated temperature used for the PIMD simulations in the T_e-PIGS method would allow for cheaper reference simulations, leaving the cost and complexity between the NNP fitting and the IBI procedure as a bigger question.

Future studies will focus on integrating our h-CMD scheme with accurate neural network potentials for predictive simulations of the vibrational spectra of the complex heterogeneous condensed phase and interfacial water.

The supplementary material contains details of fine-tuning force field parameters, average simulated box lengths, deconvoluted IR spectra for various methods, calculated IR spectra for bulk and confined water using the MB-pol and q-TIP4P/F water models, convergence of the intra-molecular bond lengths and angles of ZIF-90 and water over IBI iterations, convergence of the inter-molecular RDFs of all 19 pairs of ZIF-90 and water molecules over IBI iterations, convergence of the RMSE between reference and computed RDFs over IBI iterations, temperature-dependent RDFs obtained from h-CMD, and average simulation temperatures using different methods.

This research was partially supported by the National Science Foundation through Award No. CBET-2302618. Simulations used resources from Bridges2 at Pittsburgh Supercomputing Center through Allocation Nos. PHY230099, PHY240170, and CHE240103 from the Extreme Science and Engineering Discovery Environment (XSEDE),⁷⁰ which was supported by the National Science Foundation under Grant No. 1548562. The authors also acknowledge the HPC center at UMKC for providing the use of computing resources and support.

The authors have no conflicts to disclose.

Dil K. Limbu and Nathan London contributed equally to this work.

Dil K. Limbu: Data curation (equal); Formal analysis (equal); Methodology (supporting); Software (equal); Writing – original draft (supporting); Writing – review & editing (equal). Nathan London: Data curation (equal); Formal analysis (equal); Methodology (equal); Software (equal); Writing – original draft (supporting); Writing – review & editing (equal). Md Omar Faruque: Data curation (supporting); Formal analysis (supporting); Methodology (supporting); Software (supporting); Writing – original draft (supporting); Writing – review & editing (supporting). Mohammad R. Momeni: Conceptualization (equal); Data curation (equal); Formal analysis (equal); Funding acquisition (lead); Methodology (equal); Software (supporting); Supervision (lead); Writing – original draft (lead); Writing – review & editing (equal).

The data that support the findings of this study are available from the corresponding authors upon reasonable request. The DL_POLY Quantum software package is available for download at https://github.com/dlpolyquantum/dlpoly_quantum.