Prediction of the Aqueous Redox Properties of Functionalized Quinones Using a New QM/MM Variational Formulation

We recently proposed to couple quantum mechanics (QM) methods with 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 effects in QM calculations. Motivated by the growing interest in quinones as potential electrolytes for aqueous redox-flow batteries, we apply the QM/MDFT framework to compute the two-electron redox potentials of a series of benzoquinone/hydroquinone couples in aqueous solution. However, since these molecules are made of several dozens of atoms, their geometries are not trivial. This motivates the development of a geometry optimization procedure within the QM/MDFT framework. To this end, 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. The nuclei experience an external field generated by the electronic and classical one-particle densities. In practice, the computation of the grand potential is reduced to a sequence of density optimizations. First, the classical solvent density and the solute electronic density are optimized for a fixed solute nuclear geometry using the previously reported mixed QM/classical procedure. Subsequently, the solute geometry is optimized for a fixed solvent density. Finally, the redox potentials of a selection of benzoquinone/hydroquinone couples are computed after geometry optimizations. The predictions are in good agreement with the QM calculation using a continuum solvent model and with experimental data.