Prediction of the aqueous redox properties of functionalized quinones using a new QM/MM variational formulation

TL;DR

This work develops a variational grand-potential formulation for a quantum/classical QM/MM system in a semi-grand canonical solvent, enabling explicit solvent effects on redox-active solutes. By combining Born-Oppenheimer and zero electronic entropy approximations with a mean-field QM/MM coupling and a harmonic treatment of nuclear motion, the authors reduce the problem to coupled one-body densities and a tractable nuclear-density functional. They apply the framework to aqueous reduction of functionalized quinones, showing good agreement with electronic DFT/C-PCM predictions and experimental data, while providing rich solvent-structure details through MDFT. The approach highlights the value of explicit solvation for redox chemistry and outlines clear paths for improving QM/MM coupling and solvent polarization in future work.

Abstract

We recently proposed a method coupling quantum mechanics (QM) methods and molecular density functional theory (MDFT) to describe mixed quantum-classical systems [J. Chem. Phys. 161, 014113 (2024)]. This approach is particularly appropriate to account for solvent effect into QM calculations. We introduce a new variational formulation for the grand potential of a mixed quantum-classical system. Within the Born-Oppenheimer approximation and neglecting electronic entropy, the quantum solute is described by a product of electronic and nuclear density matrices, both depending parametrically on coordinates of the classical solvent. It can then be shown that a functional of the total density matrix satisfies a variational principle for the grand potential. Using a mean-field approximation, we express the grand potential of the mixed quantum-classical system as a variational problem which depends only on the nuclear density matrix, which experiences an external field generated by the electronic and classical one-particle densities. In practice, the grand potential is computed by a series of coupled classical and quantum DFT calculations, together with geometry optimization.

Figures (10)

• Figure 1: Quinones studied in this paper. The reference Benzoquinone, labeled 0 is circled in yellow. Quinones from 1 to 11, circled in blue, are compared to experimental data and to DFT/C-PCM calculations in figure \ref{['fig:red_pot_exp']} while quinones from 12 to 23, circled in purple are only compared to DFT/C-PCM predictions in figure \ref{['fig:red_pot_comp']}, due to the lack of experimental data.
• Figure 2: Schematic representation of the overall QM/MM scheme. The cycle start by a eDFT calculation in vacuum, which produces the electrostatic potential, $v_{\text{ext}}^{\text{MM}}$, of Eq. \ref{['eq:mf_qmmm_2_mm']}. This potential is read by the MDFT program to optimize the MM density. Then, the electrostatic potential generated by the MM region, $v_{\text{ext}}^{\text{QM}}$, of Eq. \ref{['eq:mf_qmmm_2_qm']} is read by the QM software, and the cycle is iterated until convergence of the QM/MM energy. The nuclear gradient is then computed and a geometry optimization step is performed. Once the optimal nuclear geometry is obtained, the thermodynamical corrections are computed, to obtain the final estimate of the grand potential of the system.
• Figure 3: Comparison of 2e$^-$/2H$^+$ reduction potentials (vs.Q,2H$^+$/H$_2$Q) for various quinone couples in water. The results obtained with eDFT/MDFT simulations are reported in red, while the results obtained with eDFT/C-PCM are reported in green and experimental data in black HuyAnsCavStaHam-JACS-16.
• Figure 4: Comparison of 2e$^-$/2H$^+$ reduction potentials (vs.Q,2H$^+$/H$_2$Q) for various quinone couples in water. Our results obtained with eDFT/MDFT simulations in red and the references obtained with eDFT/C-PCM in blackHuyAnsCavStaHam-JACS-16.
• Figure 5: Electronic densities (a.u.)$^{-3}$ computed with eDFT/MDFT around Q (top) and H$_2$Q (bottom) are depicted in green. The associated solvent charge density (atomic unit) are displayed in blue and red for the positive and negative regions respectively.