Modeling catalytic reaction on ligand-protected metal nanoclusters

Abstract

Monolayer protected metal nanoclusters (MPC) show great potential for catalysis, highlighting the need for unraveling their mechanistic intricacies to enhance catalyst performance. However, due to their inherent complexity, an in-depth understanding of how nanoclusters facilitate chemical transformation has remained elusive. In this work, we employed a combination of classical molecular dynamics (MD) and quantum mechanics/molecular mechanics (QM/MM) MD simulations, augmented with the on-the-fly probability-based enhanced sampling (OPES) method, to elucidate substrate binding to the MPC and a chemical reaction on the MPC surface. We investigated the oxidation of amino alcohol to oxazolidine, catalyzed by a peptide ligand functionalized Au₂₅ gold nanocluster, as a prototypical example. Classical MD simulations unveiled the crucial role of solvent and peptide-based ligands in substrate binding, while QM/MM simulations elucidated the mechanistic pathway and provided insight into the free energy landscape of the chemical reaction. This work demonstrates a computational workflow that can be applied to study similar MPC-catalyzed chemical reactions.

---

Introduction

Monolayer-protected metal nanoclusters (MPCs) have garnered considerable attention in recent years due to their unique structural and electronic properties that lie at the interface between molecules and nanoparticles.¹ These nanoclusters consist of a precisely defined metal core, typically composed of gold, silver, copper, or other transition metals, which is stabilized by a monolayer of organic ligands such as thiols, phosphines, or peptides. Owing to strong quantum confinement effects and discrete energy levels, MPCs exhibit molecule-like behavior with atomic-level precision in both structure and reactivity.² This combination of atomic precision, tunable ligand environments, and size-dependent electronic properties endows MPCs with distinctive features that set them apart from conventional metal nanoparticles. These properties have enabled MPCs to find versatile applications in fields such as catalysis, energy conversion, sensing, and bioimaging.^3–5 Among these, the application of MPCs as catalysts has attracted significant attention due to their atomic precision important for catalytic efficiency, product selectivity, and ability to fine-tune electronic properties.^6,7 MPCs have been employed across a wide range of catalytic domains, including electro-catalysis,^8–10 heterogeneous catalysis,^11,12 photo-catalysis,^13,14 and enzyme-mimic catalysis.¹⁵ Various metal (Au, Ag, Pt, Cu, etc.) nanoclusters have been found to mimic natural enzymes such as peroxidase, catalase, and superoxide dismutase.^16–22

An important advantage of MPCs as catalysts is their stability under a wide range of conditions. They retain structural integrity under high temperatures and across varying pH levels, which makes them attractive for practical catalytic applications.²³ Moreover, their modular design allows fine-tuning of catalytic behavior by altering core composition, ligand identity, and solvent interactions. This synthetic tunability provides a powerful platform for the rational design of catalysts targeting specific chemical transformations. Despite their widespread utilization as a catalyst, a comprehensive understanding of how MPCs catalyze chemical reactions has remained elusive.^24,25 An in-depth understanding of MPC-catalyzed chemical reaction for designing and synthesizing new MPCs with tuneable ligand properties for specific chemical transformations producing valuable chemicals and pharmaceutically relevant molecules.^5,7

Computer simulation is a convenient tool to decipher the atomistic and dynamic details of a chemical reaction. However, the inherent complexity of MPCs makes it challenging to study chemical reactions on MPC surfaces.^26,27 In our work, we investigated a rather complex reaction – photocatalytic oxidative cyclization of an amino alcohol, catalyzed by a peptide ligand-protected Au₂₅(S–DOP)₁₈ reported by Isozaki et al.²⁸ (DOP signifies dendritic ornithine peptide). Here the MPC plays a dual role: first, as a photocatalyst, where the gold core (Au₂₅) absorbs visible light, exciting oxygen from its ground state (³O₂) to the first excited state (¹O₂),²⁹ initiating the oxidation of amino alcohol. Second, the peptide ligands create an enzyme-like active binding site on the MPC, facilitating efficient catalysis. Taking this as a stereotypical example, we employed classical MD and QM/MM MD with enhanced sampling simulations to investigate the chemical reaction on the MPC.

Results and discussion

Ligand dynamics

Before delving into the mechanistic details of the reaction, it is crucial to understand the dynamics of the MPCs in solution. To this end, we carried out classical MD simulations of a single Au₂₅(S–DOP)₁₈ molecule in explicit chloroform (CHCl₃) and tetrahydrofuran (THF) solvents for 1 μs each. The MPC ligands exhibited contrasting dynamical behaviors in these two solvents. In the less polar solvent (CHCl₃), ligand aggregation was observed, resulting in an exposed gold surface and a non-uniform distribution of ligand density around the gold core (Fig. 1c). Conversely, in the more polar solvent (THF), the ligands interacted with the solvent molecules, leading to a much more uniform distribution around the gold core (Fig. 1d, see Video S1† for comparative dynamics in both solvents). The interactions are mostly mediated by hydrogen bonds (H-bonds); thus, more H-bonds indicate stronger ligand–ligand interactions. H-bond analysis in DMF revealed that on average ∼15–20 H-bonds were observed between the ligands, while in CHCl₃, this number increased to ∼30–35, resulting in a non-uniform ligand distribution (Fig. 1e). Additionally, in DMF, ∼15–20 H-bonds were also observed to form between the MPC ligands and the solvent molecules (Fig. 1f), which further weaken the MPC's ligand–ligand interactions. Further, root mean square fluctuations (RMSF) values for the MPC ligands are higher in the case of DMF compared to CHCl₃ (Fig. 1g), indicating that the ligands stabilize each other in CHCl₃ while they are more dynamic in DMF. Moreover, the radial distribution function (RDF) calculated for solvents around the MPC center reveals a stronger correlation in CHCl₃ than in DMF, indicating a more defined active site in CHCl₃ (Fig. 1h). This differential ligand dynamics in the two solvents impacts substrate binding, which we shall discuss in the latter section.

Fig. 1 (a) 2D representation of MPC functionalized with DOP–SH ligand, (b) structure of DOP–SH ligand, (c) representative snapshot showing structure of MPC in CHCl₃, (d) representative snapshot showing the structure of MPC in DMF (color code: gold-yellow, sulfur-green, oxygen-red, nitrogen-blue, carbon-grey, hydrogen-white), (e) Kernel density estimation (KDE) of the distribution of H-bond counts between ligands and (f) KDE of the distribution of H-bond counts between the ligands and DMF solvent, (g) root mean square fluctuation (RMSF) of MPC's ligands in presence of CHCl₃ and DMF, (h) radial distribution function (RDF) analysis. RMSF, H-bond, and RDF analysis were performed for the last 50 ns of the simulation.

Substrate binding

To understand the effect of ligands and the solvents on substrate binding, we chose three systems for studying substrate binding: (a) Au–DOP–CHCl₃, (b) Au–DOP–DMF, and (c) Au–PET–CHCl₃. The choice of DOP–SH and PET (2-phenyl ethanediol) ligands in our study is inspired by the experimental observation of the relative efficiency of DOP and PET as discussed in ref. 28. Two factors chiefly impact the overall catalytic efficiency – substrate binding and the chemical reaction, especially its cyclization step.

The substrate-bound state was obtained from an OPES simulation in which the substrate molecule was drawn from the solution to the MPC surface (Fig. S2 and Video S2, see ESI Section 1.2.2† for details). The active site for substrate binding was characterized by the presence of six Au atoms: three inner gold atoms forming a triangle (orange) and three outer gold atoms forming a second triangle (pink) positioned above the inner triangle (Fig. S3†). Subsequently, a set of 25 independent, unbiased MD simulations was carried out starting from the substrate-bound state for a maximum simulation time of 100 ns, and the trajectories were analyzed to calculate the residence time of the substrate at the binding site (Table S1 and Fig. S5†). The average residence time of the substrate in the Au–DOP–CHCl₃ system was found to be approximately 85 ns (18 out of 25 simulations rendered 100% residence time), which is significantly higher than the other two systems with resident times ranging from 3 to 12 ns (Fig. 2b). One of the simulations having 100% residence time was extended up to 1 μs and was analyzed for COM(Au)–COM(substrate) distance, which came out to be ∼7 Å (Fig. 2c, blue curve). The higher residence time in the Au–DOP–CHCl₃ system is attributed to the facile H-bonding interactions between the substrate and DOP ligands (Fig. 2c, orange curve). In contrast, the Au–PET–CHCl₃ system exhibited no H-bonding and substrate confinement due to the small size and hydrophobic nature of the PET ligand. For the Au–DOP–DMF system, the extensive H-bonding interactions between the ligands and the solvent, as well as between the substrate and the solvent, led to a lower residence time (Fig. 2b). These findings are consistent with experimental observations in ref. 28. Our study demonstrates that the differential dynamics of ligands in various solvents significantly affect the substrate binding, which is crucial for the subsequent chemical reaction, especially the cyclization step.

Fig. 2 (a) Pictorial representation of MPC–substrate complex (inset: zoomed-in picture of substrate bound to MPC), (b) average resident time of Au–DOP–CHCl₃, Au–DOP–DMF, and Au–PET–CHCl₃ systems (*18 simulations out of 25 had 100% retention time), (c) KDE analysis of important distances from classical MD simulation.

QM/MM simulation

Now we shift our attention from classical MD to the hybrid QM/MM MD simulation to model the chemical reaction occurring inside the binding pocket on the MPC. The reaction occurs in two steps. In the first step, oxygen is activated from its triplet (³O₂) ground state to the singlet (¹O₂) excited state via direct photosensitization by the MPC. It is well known that the optical band gap of Au₂₅(SR)₁₈ is larger (∼1.3 eV) than the energy of ¹O₂ (0.97 eV) which allows efficient energy transfer to ³O₂.²⁹ In the second step, the amino alcohol substrate cyclizes in the presence of a singlet oxygen to give the final product. Here we have focused on modeling the cyclization step; therefore, we directly take singlet oxygen along with amino-alcohol substrate inside the MPC's binding pocket (Fig. 3a). For QM/MM MD simulations, the substrate and singlet oxygen were included in the QM region (Fig. S6a†) and treated at the density functional theory (DFT) level, while the MPC and surrounding solvent molecules were modeled using molecular mechanics (MM). To accurately describe the chemical reaction, we employed the BLYP functional,^30,31 augmented with Grimme's D3 dispersion correction³² and Becke–Johnson damping.^33,34 A double-ζ valence polarized (DZVP) basis set was used for the QM atoms. The plane-wave expansion was carried out with a density cut-off of 600 Ry, and core electrons were represented using Goedecker–Teter–Hutter (GTH) pseudopotentials.³⁵

Fig. 3 (a) A representative snapshot from the equilibrium QM/MM MD simulation showing the interactions among the substrate, oxygen molecule, and the MPC ligand, (b) time evolution of root mean square deviation (RMSD) calculated using heavy atoms from the MPC's core and ligand interacting with the substrate, the substrate itself, and the oxygen molecule, (c) time evolution of the distance between the center of mass (COM) of the substrate and the COM of the gold core, and (d) time evolution of the distance between the nitrogen atom of the substrate and one of the oxygen atoms of the oxygen molecule.

To begin, we carried out a 50 ps unbiased QM/MM MD simulation to ensure proper binding of the substrate within the MPC's active site and to allow the system to reach equilibrium (Fig. 3a). In this context, a 50 ps QM/MM simulation is generally considered to be significantly long and often adequate for equilibrating a system.^36–39 The root mean square deviation (RMSD) analysis confirmed structural equilibration (Fig. 3b). The resulting trajectory was subsequently analyzed to characterize key interactions among the substrate, MPC, and the oxygen molecule (Fig. 3b–d, S7a and b†). Throughout the simulation, both the substrate and the oxygen molecule remained stably confined within the binding pocket of the MPC. This was supported by monitoring the time evolution of the distance between the center of mass (COM) of the substrate and the COM of the gold atoms in the MPC (Fig. 3c), which indicates stable binding. These results indicate that the system is well-equilibrated, and thus, we did not extend the simulation further. In any case, here our main focus is not on the unbiased simulation; rather, we are interested in sampling the chemical reaction, for which we apply enhanced sampling simulations that allow us to sample much longer timescales within affordable computational resources.

Notably, we observed a persistent interaction between the nitrogen atom of the substrate and one of the oxygen atoms (denoted as O2) of the oxygen molecule, maintaining an average distance of approximately 2.0 Å (Fig. 3d, blue curve). This close contact likely results from the donation of a lone pair of electrons from the nitrogen atom to the empty π* orbital of the singlet oxygen, consistent with previous reports on singlet oxygen reactivity.⁴⁰ In addition to this interaction, O₂ was also found to form a H-bond with a hydrogen atom (H1) from one of the peptide ligands, contributing to the stabilization of the oxygen molecule within the binding site (Fig. S7a†). Moreover, the other oxygen atom in the O₂ molecule (denoted as O1) engaged in a H-bond with a hydrogen atom (H) covalently bonded to the hydroxyl group of the substrate (Fig. S7b†). These non-covalent interactions collectively help stabilize the reactive configuration of the substrate and oxygen molecule within the MPC binding pocket, setting the stage for the ensuing chemical transformation.

To model the chemical reaction on the MPC surface and compute the associated free energy profile, we employed the on-the-fly probability enhanced sampling (OPES) method.⁴¹ In our simulations, we have used a surrogate approach of the path collective variable (σ) implementation by Leines and Ensing^42,43 to define the bias for the OPES simulation and monitor the reaction progress (Fig. 4a). A curve, s(σ) represents the average transition path that connects the reactant state R and product state P. The parameter σ(d), where d = d_i a set of descriptors indicates the progress along the path from R (reactant) to P (product) such that s(0) ∈ R and s(1) ∈ P. In our case, we used seven distance descriptors as shown in Fig. 4b. The initial guess path (s_g(σ_g)) was obtained using a one-dimensional Harmonic Linear Discriminant Analysis (HLDA)⁴⁴ CV (Fig. 4). This guess path in the d_i space consists of M nodes (23 in our case) represented as s_g(σ_g) → {s_j}, with j = 1, 2, …, M. The progress along the average transition path is given by

where v₁, v₂ and v₃ are defined as, v₁ = s_m − d, v₂ = d − s_m−1, and v₃ = s_m+1 − s_m; s_m is the closest path node, and s_m+1 and s_m−1 are its neighboring nodes.

Fig. 4 (a) Chemical reaction showing transformation of reactant (R) to product (P) via intermediate (I), (b) depiction of distances used as descriptors for path CV, (c) graph showing the evolution of distance descriptors along the number of nodes (M).

Path-metadynamics⁴² was performed using the OPES metadynamics⁴¹ with ‘σ’ as the path collective variable (pCV). The evolution of the distance descriptors with respect to the number of nodes is given in Fig. 4c. Additionally, we applied a restrain potential on the distance from the path |s(σ(d_k)) − d_k|, where k is the MD step, to efficiently sample the transition from the reactant to product and avoid sampling of less probable high energy states.

Multiple back-and-forth transitions between the reactant and product states were observed in the OPES simulation (Fig. S8 and Video S3†). To ensure better free energy convergence and accuracy, we conducted three independent OPES metadynamics simulations and plotted the average free energy profile and the error associated with it in Fig. 5a. The free energy profile converged smoothly, as shown in the free energy difference (ΔG) vs. simulation time plot (Fig. S11†). Furthermore, to assess the impact of a larger QM region, we performed another QM/MM MD simulation that incorporates a small portion of the ligand, particularly the region containing the peptide –NH– group and a few gold and sulfur atoms in the QM region (Fig. S6b†). The integrated Molecular Orbital Molecular Mechanics (IMOMM)⁴⁵ method was used to link atoms in the QM region with atoms in the MM region. However, this inclusion did not significantly affect the free energy profile (Fig. 5a, red curve). Further, we selected 200 structures along the reaction pathway of the larger QM region. To refine the analysis, we expanded the larger QM region to include all gold and sulfur atoms (S–CH₂) in the QM region and recalculated the QM/MM energies. The correlation plot showed reasonable consistency between the two QM region energy evaluations (Fig. S12†). This manifests that the QM model used in our study is sufficient to obtain the reaction profiles.

Fig. 5 (a) Free energy profile along the s-component of the path collective variable, (b) depiction of the substrate and an interacting ligand, (c) evolution of important distances along the reaction coordinate, (d) snapshots of important intermediates and transition states involved in the reaction.

The reaction occurs in two steps, the first step having a free energy barrier of ∼5 kcal mol⁻¹. In this step, one oxygen atom of ¹O₂ takes one hydrogen from the carbon adjacent to the nitrogen atom. The transition state (TS_1) is stabilized via H-bonding with the MPC ligand as shown in Fig. 5d. The intermediate (I) reacts via the attack of oxygen (O) on the sp² carbon along with the transfer of proton (H) to –OOH, leading to the formation of a cyclized product, oxazolidine and H₂O₂. The evolution of important distances as a function of the reaction coordinate is depicted in Fig. 5b and c. In the intermediate I, O₂ begins to interact strongly with the substrate's hydrogen (H), and the substrate's oxygen (O) interacts with the ligand's hydrogen (H1). Fig. 5d shows important geometries along the reaction pathway. A slight difference in the distance evolution was observed between the small and large QM regions (Fig. S10†), which is expected due to the stochastic nature of molecular simulations. Furthermore, substrate–ligand interactions are inherently dynamic, fluctuating within a range characteristic of the reactant, intermediate, and product states. These variations reflect the flexible nature of the system.

Conclusion

In conclusion, we have developed an advanced modeling approach to study the chemical reactions catalyzed by monolayer-protected metal nanoclusters (MPCs). Utilizing classical molecular dynamics simulations, we explored the ligand dynamics of MPCs in different solvents. The differential reaction rates observed with MPCs protected by a simple ligand (PET) compared to a peptide-based ligand (DOP–SH) were attributed to variations in the substrate's residence time on the MPC's surface. Further, we used the hybrid QM/MM technique, combined with an enhanced sampling method, OPES with path collective variable (pCV) to model the reaction on MPC's surface. The detailed QM/MM modeling provided crucial insights into how the ligand stabilizes the transition state and intermediates. The use of a surrogate model of the path collective variable allowed for efficient sampling of the overall mechanism and accurate determination of the free energy surface. While the converged free energy profile provided the free energy barrier corresponding to the cyclization step, obtaining the reaction rate from it remains challenging. This is due to the limitation of accurately determining the pre-exponential factor in the rate constant equation. Moreover, even if one succeeds, a direct comparison of the computationally obtained rate with the experimentally reported one is difficult, since the former corresponds to the cyclization step while the latter refers to the overall reaction rate. Nevertheless, this study offers a comprehensive understanding of the role of MPCs as catalysts and the effects of ligands and solvents on their catalytic performance. The workflow we developed not only deepens our understanding of MPC-catalyzed reactions but also provides a robust framework that can be applied to other complex catalytic systems.

Computational details

Classical MD simulation

Classical MD simulations were carried out using GROMACS 2021.4 software⁴⁶ patched with plumed 2.8.0 code.^47,48 NanoModeler^49,50 server was used to generate the topology of MPCs, which uses empirical force-field parameters for bonded and non-bonded interactions of gold–sulfur motifs generated by Pohjolainen et al.,⁵¹ and Heinz et al.,⁵² respectively. The ligands and solvents were modeled using the general AMBER force field (GAFF).⁵³ Further details of the systems, force fields, and simulation methods are provided in the ESI.†

QM/MM MD simulation

All QM/MM MD simulations were performed using CP2K 2022.1 software.⁵⁴ The system was first minimized using the steepest descent algorithm, followed by equilibration at 300 K in the NVT ensemble for 50 ps of NVT equilibration at the DFT level (BLYP-D3(BJ)/DZVP). All production QM/MM MD simulations were carried out at the DFT level in the NVT ensemble, using a time step of 0.5 fs. Temperature was maintained using a Nosé–Hoover thermostat with a coupling constant of 0.1 ps. Further details are provided in the ESI.†

Data availability

All input files required to run the simulations presented in this work can be found in a public GitHub repository (https://github.com/vikast282/mpc_QMMM) and PLUMED-NEST⁵⁵ (ID: 24.015).

Author contributions

Project design and planning: VT and TK, execution of the project and formal analysis: VT, and manuscript writing and editing: both authors.

Conflicts of interest

The authors declare no conflict of interest.

Acknowledgements

VT thanks the Ministry of Education, Govt. of India, for the PMR Fellowship. TK acknowledges the Anusandhan National Research Foundation (ANRF) formally known as Science and Engineering Research Board (SERB), New Delhi, India, for the Start-up Research Grant (File No. SRG/2022/000969). We also acknowledge IIT Delhi for the Seed Grant. We thank the IIT Delhi HPC facility (FIST) for computational resources.

Notes and references

1. I. Chakraborty and T. Pradeep, Chem. Rev., 2017, 117, 8208–8271.
2. R. Jin, C. Zeng, M. Zhou and Y. Chen, Chem. Rev., 2016, 116, 10346–10413.
3. K. Saha, S. S. Agasti, C. Kim, X. Li and V. M. Rotello, Chem. Rev., 2012, 112, 2739–2779.
4. Y. Tao, M. Li, J. Ren and X. Qu, Chem. Soc. Rev., 2015, 44, 8636–8663.
5. X. Du and R. Jin, ACS Nano, 2019, 13, 7383–7387.
6. A. Longo, E. J. J. Boed, N. Mammen, M. Linden, K. Honkala, H. Häkkinen, P. E. Jongh and B. Donoeva, Chem.–Eur. J., 2020, 26, 7051–7058.
7. A. Pecina, D. Rosa-Gastaldo, L. Riccardi, S. Franco-Ulloa, E. Milan, P. Scrimin, F. Mancin and M. De Vivo, ACS Catal., 2021, 11, 8736–8748.
8. T. Kawawaki, T. Okada, D. Hirayama and Y. Negishi, Green Chem., 2024, 26, 122–163.
9. M. H. Naveen, R. Khan and J. H. Bang, Chem. Mater., 2021, 33, 7595–7612.
10. K. Kwak and D. Lee, Acc. Chem. Res., 2018, 52, 12–22.
11. L. Liu and A. Corma, Chem. Rev., 2018, 118, 4981–5079.
12. R. Jin, G. Li, S. Sharma, Y. Li and X. Du, Chem. Rev., 2020, 121, 567–648.
13. Y. Du, C. Li, Y. Dai, H. Yin and M. Zhu, Nanoscale Horiz., 2024, 9, 1262–1278.
14. Y. Liu, Y. Wang and N. Pinna, ACS Mater. Lett., 2024, 6, 2995–3006.
15. X. Pan, Y. Yao, M. Zhang, X. Yuan, Q. Yao and W. Hu, Nanoscale, 2024, 16, 8196–8215.
16. H. Liu, Y. Li, S. Sun, Q. Xin, S. Liu, X. Mu, X. Yuan, K. Chen, H. Wang, K. Varga, W. Mi, J. Yang and X.-D. Zhang, Nat. Commun., 2021, 12, 114.
17. S. Wang, W. Chen, A.-L. Liu, L. Hong, H.-H. Deng and X.-H. Lin, ChemPhysChem, 2012, 13, 1199–1204.
18. Y. Tao, Y. Lin, Z. Huang, J. Ren and X. Qu, Adv. Mater., 2013, 25, 2594–2599.
19. L. Hu, Y. Yuan, L. Zhang, J. Zhao, S. Majeed and G. Xu, Anal. Chim. Acta, 2013, 762, 83–86.
20. C.-P. Liu, T.-H. Wu, Y.-L. Lin, C.-Y. Liu, S. Wang and S.-Y. Lin, Small, 2016, 12, 4127–4135.
21. C. Liu, Y. Cai, J. Wang, X. Liu, H. Ren, L. Yan, Y. Zhang, S. Yang, J. Guo and A. Liu, ACS Appl. Mater. Interfaces, 2020, 12, 42521–42530.
22. L. Chen, L. Klemeyer, M. Ruan, X. Liu, S. Werner, W. Xu, A. Koeppen, R. Bücker, M. G. Gonzalez, D. Koziej, W. J. Parak and I. Chakraborty, Small, 2023, 19, 2206772.
23. V. Sudheeshkumar, K. O. Sulaiman and R. W. J. Scott, Nanoscale Adv., 2020, 2, 55–69.
24. Y. Du, H. Sheng, D. Astruc and M. Zhu, Chem. Rev., 2020, 120, 526–622.
25. S. Li, N.-N. Li, X.-Y. Dong, S.-Q. Zang and T. C. W. Mak, Chem. Rev., 2024, 124, 7262–7378.
26. S. Malola and H. Häkkinen, Nat. Commun., 2021, 12, 2197.
27. V. Rojas-Cervellera, L. Raich, J. Akola and C. Rovira, Nanoscale, 2017, 9, 3121–3127.
28. K. Isozaki, R. Ueno, K. Ishibashi, G. Nakano, H. Yin, K. Iseri, M. Sakamoto, H. Takaya, T. Teranishi and M. Nakamura, ACS Catal., 2021, 11, 13180–13187.
29. H. Kawasaki, S. Kumar, G. Li, C. Zeng, D. R. Kauffman, J. Yoshimoto, Y. Iwasaki and R. Jin, Chem. Mater., 2014, 26, 2777–2788.
30. A. D. Becke, Phys. Rev. A, 1988, 38, 3098–3100.
31. C. Lee, W. Yang and R. G. Parr, Phys. Rev. B:Condens. Matter Mater. Phys., 1988, 37, 785–789.
32. S. Grimme, J. Antony, S. Ehrlich and H. Krieg, J. Chem. Phys., 2010, 132, 154104.
33. S. Grimme, S. Ehrlich and L. Goerigk, J. Comput. Chem., 2011, 32, 1456–1465.
34. E. R. Johnson and A. D. Becke, J. Chem. Phys., 2006, 124, 174104.
35. S. Goedecker, M. Teter and J. Hutter, Phys. Rev. B:Condens. Matter Mater. Phys., 1996, 54, 1703–1710.
36. J. P. Rivas-Fernández, M. Vuillemin, B. Pilgaard, L. J. Klau, F. Fredslund, C. Lund-Hanssen, D. H. Welner, A. S. Meyer, J. P. Morth and F. Meilleur, et al. , Nat. Commun., 2025, 16, 2670.
37. Z. Wang, X. Fu, W. Diao, Y. Wu, C. Rovira and B. Wang, Chem. Sci., 2025, 16, 3173–3186.
38. M. Sagiroglugil, Q. Liao, A. Planas and C. Rovira, ChemCatChem, 2024, 16, e202400769.
39. M. Sagiroglugil, A. Nin-Hill, E. Ficko-Blean and C. Rovira, ACS Catal., 2024, 14, 16897–16904.
40. J. S. Nam, Y. Hong, C. G. Lee, T. I. Kim, C. Lee, D.-H. Roh, I. S. Lee, S. Kweon, G. Ahn, S. K. Min, B.-S. Kim and T.-H. Kwon, JACS Au, 2022, 2, 933–942.
41. M. Invernizzi and M. Parrinello, J. Phys. Chem. Lett., 2020, 11, 2731–2736.
42. G. D. Leines and B. Ensing, Phys. Rev. Lett., 2012, 109, 020601.
43. A. Pérez de Alba Ortíz, A. Tiwari, R. C. Puthenkalathil and B. Ensing, J. Chem. Phys., 2018, 149, 072320.
44. D. Mendels, G. Piccini and M. Parrinello, J. Phys. Chem. Lett., 2018, 9, 2776–2781.
45. A. Genest, A. Woiterski, S. Krüger, A. M. Shor and N. Rösch, J. Chem. Theory Comput., 2006, 2, 47–58.
46. H. Berendsen, D. van der Spoel and R. van Drunen, Comput. Phys. Commun., 1995, 91, 43–56.
47. G. A. Tribello, M. Bonomi, D. Branduardi, C. Camilloni and G. Bussi, Comput. Phys. Commun., 2014, 185, 604–613.
48. M. Bonomi, D. Branduardi, G. Bussi, C. Camilloni, D. Provasi, P. Raiteri, D. Donadio, F. Marinelli, F. Pietrucci, R. A. Broglia and M. Parrinello, Comput. Phys. Commun., 2009, 180, 1961–1972.
49. S. Franco-Ulloa, L. Riccardi, F. Rimembrana, M. Pini and M. D. Vivo, J. Chem. Theory Comput., 2019, 15, 2022–2032.
50. S. Franco-Ulloa, L. Riccardi, F. Rimembrana, E. Grottin, M. Pini and M. D. Vivo, J. Chem. Theory Comput., 2023, 19, 1582–1591.
51. E. Pohjolainen, X. Chen, S. Malola, G. Groenhof and H. Häkkinen, J. Chem. Theory Comput., 2016, 12, 1342–1350.
52. H. Heinz, R. A. Vaia, B. L. Farmer and R. R. Naik, J. Phys. Chem. C, 2008, 112, 17281–17290.
53. J. Wang, R. M. Wolf, J. W. Caldwell, P. A. Kollman and D. A. Case, J. Comput. Chem., 2004, 25, 1157–1174.
54. T. D. Kühne, M. Iannuzzi, M. Del Ben, V. V. Rybkin, P. Seewald, F. Stein, T. Laino, R. Z. Khaliullin, O. Schütt, F. Schiffmann, D. Golze, J. Wilhelm, S. Chulkov, M. H. Bani-Hashemian, V. Weber, U. Borštnik, M. Taillefumier, A. S. Jakobovits, A. Lazzaro, H. Pabst, T. Müller, R. Schade, M. Guidon, S. Andermatt, N. Holmberg, G. K. Schenter, A. Hehn, A. Bussy, F. Belleflamme, G. Tabacchi, A. Glöβ, M. Lass, I. Bethune, C. J. Mundy, C. Plessl, M. Watkins, J. VandeVondele, M. Krack and J. Hutter, J. Chem. Phys., 2020, 152, 194103.
55. M. Bonomi, G. Bussi, C. Camilloni, G. A. Tribello, P. Banáš, A. Barducci, M. Bernetti, P. G. Bolhuis, S. Bottaro, D. Branduardi, R. Capelli, P. Carloni, M. Ceriotti, A. Cesari, H. Chen, W. Chen, F. Colizzi, S. De, M. De La Pierre, D. Donadio, V. Drobot, B. Ensing, A. L. Ferguson, M. Filizola, J. S. Fraser, H. Fu, P. Gasparotto, F. L. Gervasio, F. Giberti, A. Gil-Ley, T. Giorgino, G. T. Heller, G. M. Hocky, M. Iannuzzi, M. Invernizzi, K. E. Jelfs, A. Jussupow, E. Kirilin, A. Laio, V. Limongelli, K. Lindorff-Larsen, T. Löhr, F. Marinelli, L. Martin-Samos, M. Masetti, R. Meyer, A. Michaelides, C. Molteni, T. Morishita, M. Nava, C. Paissoni, E. Papaleo, M. Parrinello, J. Pfaendtner, P. Piaggi, G. Piccini, A. Pietropaolo, F. Pietrucci, S. Pipolo, D. Provasi, D. Quigley, P. Raiteri, S. Raniolo, J. Rydzewski, M. Salvalaglio, G. C. Sosso, V. Spiwok, J. Šponer, D. W. H. Swenson, P. Tiwary, O. Valsson, M. Vendruscolo, G. A. Voth and A. White, The PLUMED consortium, Nat. Methods, 2019, 16, 670–673.

---

Footnote

† Electronic supplementary information (ESI) available: Details of the force field, OPES simulation parameters, collective variable, and QM/MM simulation; Video_S1.mp4: video showing ligand dynamics in CHCl3 and DMF; Video_S2.mp4: video showing substrate binding to the MPC; Video_S3.mp4: video showing QM/MM modeled reaction mechanism. See DOI: https://doi.org/10.1039/d5sc00421g

---