2.3. QM/MM Molecular Dynamics

In order to choose a representative configuration to carry out the QM/MM MD study, we performed a clustering analysis [34] on the frames extracted from the above-mentioned classical MD simulation. Non-overlapping clusters were assigned using a local root mean square deviation (RMSD) (computed taking into account all protein residues within a 5 Å cutoff from the chromophore) as the measure of the distance between configurations. The central structure from the most populated cluster was then used to initiate the QM/MM MD simulation.

The QM/MM implementation employed in this study combines the use of the QM program QUICKSTEP [35] and the MM driver FIST, both forming part of the CP2K package [36]. In this code, the general QM/MM scheme is based on a real space multigrid technique to compute the electrostatic coupling between the QM and MM regions [37,38]. A quantum region (comprising the pCA chromophore and the side chains of the residues either involved in short HBs with pCA phenolic oxygen (Y42 and E46) or directly attached to it (C69)) was treated at the DFT level. The remaining part of the system, including water molecules and counterions, was modeled at the classical level using the AMBER force field. The valence of the terminal QM atoms was saturated by the addition of capping hydrogen atoms. A dual basis set, Gaussian and plane-wave (GPW) formalism, was employed to compute the interaction energy within the atoms belonging to the QM region. In particular, a double-ζ valence basis set augmented with a set of polarization functions (DZVP) [39] was used in order to obtain an accurate description of the wave function of the QM subsystem, while an auxiliary plane-wave basis set expanded up to a density cutoff of 320 Ry was utilized to converge the electron density in conjunction with Goedecker–Teter–Hutter (GTH) pseudopotentials [40,41] to describe the core electrons. Exchange and correlation energies were computed within the Generalized Gradient Approximation (GGA) using the BLYP functional [42,43]. The QM/MM MD simulations were performed with Born−Oppenheimer dynamics in the canonical (NVT) ensemble and using an integration time step of 0.5 fs. A self-consistent field (SCF) convergence threshold of 10⁻⁶ au was used to optimize the wave function at every step of the dynamics. Starting from the snapshot extracted from the forcefield-based dynamics, the system was equilibrated at the QM/MM level for about 2 ps. Then, the simulation was extended by 10 ps using stochastic velocity rescaling thermostats [44], independently for the QM and MM regions, to maintain the temperature of both at 303.15 K. [32].