Quantum Dynamical Effects in 2-(2-Furyl)-3-Hydroxychromone Using Path Integral Molecular Dynamics ()
1. Introduction
Understanding the mechanism of fluorescent compounds in solution requires accurate modeling approaches that capture both structural and dynamical features. Classical molecular dynamics (MD) simulations have long been used for this purpose due to their computational efficiency and scalability. However, such methods neglect nuclear quantum effects (NQEs), which can significantly alter the behavior of light atoms, particularly hydrogen. Quantum phenomena such as zero-point energy, delocalization, and tunneling play a crucial role in modulating molecular interactions, including hydrogen bonding, vibrational modes, and diffusion [1].
Path Integral Molecular Dynamics (PIMD) provides a robust framework for incorporating NQEs into molecular simulations. In this approach, each nucleus is represented as a ring polymer consisting of several fictitious beads linked by harmonic springs, enabling the sampling of quantum fluctuations at finite temperatures by mapping the quantum system onto an extended classical-like phase space [2]. PIMD has been instrumental in uncovering reaction mechanisms inaccessible to classical transition state theory. In contrast, classical MD, which treats nuclei as point masses, often overestimates free energy barriers in processes such as proton transfer and fails to reproduce kinetic isotope effects, rendering it unsuitable in cases where NQEs are dominant [3].
Another strategy is the hybrid Quantum Mechanics/Molecular Mechanics (QM/MM) framework, which combines quantum mechanical precision in the reactive region with the computational efficiency of classical force fields for the environment. However, standard QM/MM implementations can suffer from abrupt energy changes and unphysical force discontinuities at the region boundary, which can severely distort the potential energy landscape for sensitive processes like proton transfer [4]. For instance, a recent QM/MM study on FHC itself reported significant artifacts at the boundary, underscoring the challenge for modeling its ESIPT process accurately. When paired with quantum dynamics techniques such as PIMD for the QM region, this method can more effectively capture NQEs in complex condensed-phase processes, including Excited State Intramolecular Proton Transfer (ESIPT) [5]. Nevertheless, the computational expense of such ab initio PIMD simulations remains prohibitive for extensive sampling of the solvent environment.
Recent studies highlight the limitations of standard QM/MM methods, where abrupt changes in energy and forces occur at the QM/MM boundary, distorting the potential energy landscape for sensitive processes like proton transfer. While hybrid QM/MM paired with PIMD can mitigate some issues, the computational cost and challenges in modeling solvent quantum effects remain significant. This study focuses on a full PIMD approach with a carefully parameterized force field to capture NQEs uniformly across the entire system, avoiding such artifacts and providing a consistent quantum-mechanical treatment of the solvent microenvironment. FHC’s ESIPT process is highly sensitive to its environment, with experimental observations showing slowed proton transfer kinetics in protic solvents due to enhanced hydrogen bonding. This slowdown has been attributed to a stronger, more structured hydrogen-bonding network that alters the ESIPT potential energy barrier. This work aims to quantify how NQEs in the ground state contribute to this sensitivity by shaping this precise microenvironment.
In this study, we use both classical MD and PIMD to explore the influence of NQEs on 2-(2-furyl)-3-hydroxychromone (FHC), a fluorescent molecule with ESIPT character, solvated in methanol. FHC belongs to the family of 3-hydroxychromone derivatives and exhibits ESIPT an ultrafast photophysical process where a proton is transferred from the hydroxyl group to an adjacent carbonyl oxygen upon photoexcitation, generating a tautomeric form with distinct emission characteristics [6]. Figure 1 shows the molecular structure of FHC in its normal form. Its pronounced sensitivity to the surrounding medium, such as the experimentally observed slowdown of ESIPT kinetics in protic environments, makes it an effective probe for supramolecular assemblies and biological environments [6] [7]. This sensitivity suggests that a quantum treatment of nuclear motion, even in the ground state, is crucial for a realistic model of its solvation structure and dynamics.
Through a comparative analysis of structural distributions, diffusion characteristics, and infrared spectra obtained from classical MD and PIMD simulations, we seek to assess the impact of NQEs on the ground-state dynamical and spectroscopic behavior of FHC in the condensed phase. This work provides an essential foundation for understanding how quantum effects modulate the microenvironment that governs FHC’s photophysical properties and will inform future more expensive QM/MM-PIMD studies of the excited-state reaction itself.
Figure 1. FHC in the normal form.
2. Methodology
We studied a single molecule of 2-(2-furyl)-3-hydroxychromone (FHC) solvated in methanol by combining classical molecular dynamics (MD) and path integral molecular dynamics (PIMD) simulations. All simulations were carried out at 300 K and 1 atm to closely reproduce ambient experimental conditions [8].
2.1. System Preparation and Force Field Parameterization
The studied system consisted of one FHC molecule solvated in a cubic box of methanol with an initial side length of 20 Å, containing approximately 313 solvent molecules [9]. The final box size was determined during the NPT equilibration phase to achieve the correct experimental density of methanol (0.791 g/cm3 at 300 K) [10].
Molecular mechanics parameters (bonds, angles, dihedrals) for FHC were assigned using the GAFF2 (Generalized Amber Force Field 2) within the AMBER force field family [1]. The GAFF2 parameters were validated against ab initio calculations at the PBE0/cc-pVTZ level to ensure accuracy for hydrogen bonding interactions [2]. Dispersion corrections were incorporated using the DFT-D3 method to account for van der Waals interactions accurately [3].
The atomic partial charges were derived quantum mechanically to ensure accuracy for electrostatic interactions. The electrostatic potential (ESP) was computed at the PBE0/cc-pVTZ level of theory on the ground-state (S0) optimized geometry [4]. Atomic partial charges were then obtained by fitting the ESP using the restrained electrostatic potential (RESP) procedure [5] [6].
The topology (prmtop) and initial coordinate (inpcrd) files were generated using the tleap module of AmberTools22 [7]. The solvent box was packed around the solute using Packmol, ensuring uniform solvation. Prior to dynamics, the system underwent energy minimization using the steepest descent algorithm to relax steric clashes and optimize the initial geometry, which is essential for numerical stability. Periodic boundary conditions were applied in all three spatial directions to mimic a bulk liquid environment [8].
2.2. Classical Molecular Dynamics Protocol
After energy minimization, the system was equilibrated in two phases, NPT Ensemble with 1 ns equilibration at 300 K and 1 atm to allow the density to relax to the correct experimental value and NVT Ensemble with subsequent 1 ns equilibration at 300 K using the average volume from the NPT phase [9].
Production dynamics were then carried out in the NVT ensemble for >20 ns to ensure adequate sampling for dynamical properties.
All simulations were performed using the PMEMD.CUDA engine of the AMBER20 software suite [6]. Long-range electrostatic interactions were treated using the Particle Mesh Ewald (PME) method [7] with a 9.0 Å direct space cutoff. A Langevin thermostat with a collision frequency (gamma_ln) of 2.0 ps−1 was used for temperature control, which provides robust thermal regulation suitable for dynamical properties [8]. The equations of motion were integrated with a 2 fs timestep, with all bonds involving hydrogen atoms constrained using the SHAKE algorithm [9].
2.3. Path Integral Molecular Dynamics (PIMD) Protocol
PIMD was employed to incorporate nuclear quantum effects (NQEs), which are critical for accurately modeling hydrogen bonding and proton delocalization in this system [2] [10]. The methodology is based on Feynman’s path integral formulation, where each quantum nucleus is represented by a ring polymer of P classical beads connected by harmonic springs, effectively mapping the quantum system onto an extended classical system [11] [12].
In this study, we used P = 32 beads per atom, which provides a satisfactory compromise between quantum accuracy and computational cost for properties at 300 K [13]. All other simulation conditions (box size, number of solvent molecules, force field) were kept identical to the classical MD simulations to ensure a direct and consistent comparison.
A critical aspect of PIMD is efficient sampling of the quantum canonical ensemble. We applied the Path Integral Langevin Equation (PILE) thermostat to the normal modes of the ring polymers [14] [15]. This approach optimally thermostats the internal modes without perturbing the physical dynamics of the centroid, which is essential for accurate dynamical properties [16].
The equations of motion were integrated using the multi-timestep RESPA algorithm [17]. The timestep (dt) was selected based on the fastest vibrational mode of the ring polymer. The fundamental period Tp is given by Tp = h/(kB * T * P). To ensure numerical stability and proper resolution of quantum fluctuations, the timestep was set to dt ≈ Tp/20 ≈ 0.5 fs [18].
2.3.1. PIMD Equilibration
Before production, the PIMD system underwent careful equilibration under NPT conditions (300 K, 1 atm) using the same protocol as classical MD [19]. The system was considered equilibrated once potential energy and density had stabilized [20].
2.3.2. PIMD Production
The production simulation was conducted in the NVT ensemble using the average volume from PIMD-NPT equilibration [21]. The PILE thermostat and 0.5 fs timestep were maintained. The total production time was >4 ns (representing >128 ns cumulative bead dynamics), allowing for statistically robust estimates of key observables [20].
The resulting trajectories from both classical MD and PIMD were analyzed to quantify the impact of NQEs on structural, dynamical, and spectroscopic behavior of FHC in methanol [22].
3. Results and Discussions
3.1. Structural Properties and System Equilibration
The structural stability of 2-(2-furyl)-3-hydroxychromone (FHC) in methanol was evaluated by monitoring key thermodynamic parameters. Figure 2 shows the evolution of temperature and density during the equilibration and production phases, which stabilized at 300 K and 0.79 g/cm3, respectively [23]. The potential and kinetic energy components remained stable throughout the production phase. Statistical analysis of the final 10 ns of the trajectory showed no systematic drift and Gaussian distributions for these quantities (Figure 3), confirming the system had reached robust thermal and structural equilibrium [24].
Figure 2. Evolution of temperature and density during the equilibration and production phases of the simulations.
Figure 3. Temperature and density distributions from classical molecular dynamics (MD) simulations.
A direct comparison of the potential energy evolution from classical MD and PIMD simulations reveals a crucial finding, the average potential energy is consistently lower in the PIMD simulation (Figure 4). This is a direct manifestation of the inclusion of nuclear quantum effects (NQEs), primarily zero-point energy (ZPE). The ZPE effectively lifts the system to a higher vibrational baseline, allowing it to sample lower regions of the classical potential energy well that are inaccessible to classical nuclei, which are frozen at their minimum at 0 K. This result underscores the necessity of accounting for quantum fluctuations when modeling systems dominated by light atoms like hydrogen [24] [25].
Figure 4. Evolution of potential energies of MD and PIMD simulations.
Further evidence of NQEs is found in the atomic coordinate and velocity distributions. As illustrated in Figure 5, the distributions derived from PIMD are broader than those from classical MD. This is the expected signature of quantum mechanical delocalization, where nuclei are not point particles but are instead represented by smeared probability distributions [26]. The PIMD framework, by representing each nucleus as a ring polymer, explicitly captures this delocalization, leading to a wider exploration of phase space. This fundamental difference is critical for accurately modeling properties that depend on nuclear position, such as hydrogen bond lengths and vibrational spectra [27].
3.2. Diffusion Coefficient from Velocity Autocorrelation
Function (VACF)
The molecular diffusion coefficient (D) was calculated from the velocity autocorrelation function (VACF) of the center of mass of FHC using the Green-Kubo relation D = (1/3) ∫⟨v(0)*v(t)⟩ dt [28].
The VACF decay (Figure 6) shows a clear difference between the two simulation methods. The VACF from the PIMD simulation decays more slowly than its
Figure 5. Distribution of atomic coordinates and velocities from MD and PIMD simulations.
Figure 6. Velocity Autocorrelation Function (VACF) of the center of mass of FHC.
classical counterpart. A slower decay of the VACF indicates that the molecule’s velocity remains correlated for a longer time, which is directly linked to a reduction in the computed diffusion coefficient [29].
The calculated diffusion coefficients are D_MD = 1.8 × 10−5 cm2/s and D_PIMD = 1.4 × 10−5 cm2/s. A summary of these and other key properties is provided in Table 1. This ~22% reduction in diffusion under PIMD is a significant quantum effect [30]. It can be attributed to the interplay of zero-point energy and delocalization. Stronger Effective Hydrogen Bonding, Quantum delocalization allows protons to sample shorter O-H-O distances, effectively strengthening and tightening the hydrogen bond network between FHC and the methanol solvent [31]. The more structured solvation shell around FHC, reinforced by quantum effects, creates a more restrictive cage, hindering the translational motion of the solute molecule [32].
This result aligns with the experimental observation that FHC’s ESIPT kinetics are slowed in protic, hydrogen-bonding environments [33]. The PIMD simulation suggests that part of this environmental modulation may stem from quantum-mechanically enhanced solvent-solute interactions that restrict mobility even in the ground state.
Table 1. Comparison of structural and dynamical properties from classical MD and PIMD simulations.
Properties |
Classical MD |
PIMD |
Temperature (K) |
300 |
300 |
Density (g/cm3) |
0.78 |
0.78 |
Potential Energy (kcal/mol) |
Higher |
Lower |
Velocity Autocorrelation |
Faster |
Slower |
Diffusion Coefficient (cm2/s) |
1.8 × 10−5 |
1.4 × 10−5 |
3.3. Infrared Absorption Spectra
The infrared (IR) absorption spectrum was computed from the Fourier transform of the autocorrelation function of the total dipole moment derivative of the system [21]. As shown in Figure 7, the PIMD-derived spectrum exhibits broader absorption bands and a general red-shift compared to the classical MD spectrum.
Figure 7. IR spectra from MD and PIMD.
These spectral changes are classic signatures of nuclear quantum effects, arising from the sampling of a wider distribution of molecular geometries due to quantum nuclear delocalization, effectively incorporating anharmonicity [34]. Primarily caused by zero-point energy, which lowers the effective barrier for vibrational modes, reducing their fundamental frequencies.
These effects are most pronounced for vibrations involving hydrogen atoms due to their large quantum fluctuations. The improved agreement of the PIMD spectrum with expected experimental line shapes highlights its superiority over classical MD for predicting spectroscopic properties. It provides a more realistic representation of vibrational modes, anharmonicity, and the influence of the quantum mechanical hydrogen bond network [35].
3.4. Impact of Nuclear Quantum Effects on FHC’s Microenvironment
The incorporation of nuclear quantum effects through PIMD paints a consistent picture of FHC’s behavior in methanol, NQEs promote a tighter, more structured, and more strongly interacting solvation environment. Energetically, the system benefits from zero-point energy, sampling lower potential energies. Structurally, nuclei are delocalized, leading to smoother, broader distributions of atomic positions, particularly for hydrogen bonds. Dynamically, this translates to stronger solute-solvent interactions that restrict translational diffusion. Spectroscopically, the vibrations are softened and broadened, reflecting the anharmonic quantum nature of the potential energy surface.
This comprehensive description of the ground-state microenvironment is crucial for understanding FHC’s photophysics. The quantum-mechanically enhanced hydrogen bonding captured by PIMD likely creates a solvation shell that pre-organizes the molecule and modifies the potential energy landscape for the ESIPT reaction. This provides a foundational, atomistic explanation for the profound sensitivity of FHC’s ESIPT kinetics to its solvent environment, as observed experimentally [36].
4. Conclusions
This study demonstrates that NQEs fundamentally alter the solvation environment of FHC, leading to stronger hydrogen bonding, reduced diffusion, and modified spectroscopic properties. These findings provide an atomistic explanation for experimental observations of solvent-dependent ESIPT kinetics. Future work will involve a QM/MM-PIMD approach to directly model the excited-state dynamics, combining the accuracy of electronic structure methods with the quantum treatment of nuclei. This multi-scale strategy will pave the way for a comprehensive understanding of photochemical processes in complex environments.
Our comparative analysis reveals that NQEs induce consistent and significant changes, energetically, the system samples lower potential energies due to the inclusion of zero-point energy. Structurally, atomic position and velocity distributions are broader, reflecting quantum mechanical smearing, particularly for hydrogen atoms. Dynamically, the self-diffusion coefficient of FHC is reduced by approximately 22% under PIMD. This is a direct consequence of quantum effects strengthening the hydrogen-bonding network and enhancing solvent caging. Spectroscopically, the simulated infrared spectrum exhibits broader and red-shifted absorption bands, aligning more closely with expected experimental results by incorporating vibrational anharmonicity and quantum fluctuations.
These findings demonstrate that NQEs are not a minor perturbation but a fundamental factor defining the solvation microenvironment of FHC. The PIMD approach provides a more realistic depiction of intermolecular interactions in environments dominated by hydrogen bonding and light atoms.
While the ESIPT reaction itself is an excited-state process involving intrinsic quantum dynamics like proton tunneling, our results provide crucial ground-state context. The quantum-mechanically enhanced solvation shell observed here likely pre-organizes the FHC molecule and modifies the potential energy landscape it experiences upon photoexcitation. This offers a foundational, atomistic explanation for the experimentally observed sensitivity of FHC’s ESIPT kinetics to its solvent environment, where protic solvents slow down the reaction rate.
Looking forward, this work establishes a robust classical PIMD framework for understanding the solvent’s role. The logical and essential next step is to combine these insights with a quantum mechanical treatment of the electrons. A multi-scale QM/MM-PIMD approach, where the FHC core is treated with an electronic structure method and the solvent is modeled with the quantum-accurate force field used here, represents the optimal strategy to directly simulate the quantum dynamics of the ESIPT reaction in a realistic, condensed-phase environment. This would unite the accurate solvent representation achieved in this study with a quantum mechanical description of the reaction itself, paving the way for a complete quantum-mechanical understanding of photochemical processes in complex systems.
Data Availability Statement
The molecular dynamics simulation data, including trajectory files, input parameters, and analysis scripts generated during this study, are available from the corresponding author upon reasonable request. Due to the size of the datasets, they are not publicly archived but can be provided to facilitate verification and further research.
Ethics Statement
This study did not involve any research on human participants, human data, or animals. Therefore, ethical approval was not required.