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

Maxime Labat, Guillaume Jeanmairet, Emmanuel Giner

► To cite this version:

Maxime Labat, Guillaume Jeanmairet, Emmanuel Giner. Prediction of the aqueous redox properties of functionalized quinones using a new QM/MM variational formulation. Journal of Chemical Theory and Computation, 2025, 10.1021/acs.jctc.5c00976. [https://doi.org/10.1021/acs.jctc.5c00976](https://doi.org/10.1021/acs.jctc.5c00976)

HAL Id: hal-05374802

https://hal.science/hal-05374802v1

Submitted on 20 Nov 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>
# Prediction of the aqueous redox properties of functionalized quinones using a new QM/MM variational formulation.

Maxime Labat,¹ Guillaume Jeanmairet*,¹,², a) and Emmanuel Giner*³, b)

¹) Sorbonne Université, CNRS, Physico-Chimie des Electrolytes et Nanosystèmes Interfaciaux, PHENIX, F-75005 Paris, France

²) Réseau sur le Stockage Électrochimique de l'Énergie (RS2E), FR CNRS3459, 80039 Amiens Cedex, France

³) Sorbonne Université, CNRS, Laboratoire de Chimie Théorique, Sorbonne Université, F-75005 Paris, France

## Abstract

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 effect into 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-electrons 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 quantum mechanics/classical procedure. Subsequently, the solute geometry is optimized for a fixed solvent configuration. Finally, the redox potentials of a selection of Benzoquinone/Hydroquinone couples are computed after geometry optimizations. The predictions are in good agreement with QM calculation using a continuum solvent model and with experimental data.

## I. INTRODUCTION

Recent advancements in redox flow batteries (RFBs) highlight their potential for large-scale energy storage, particularly for intermittent renewable sources like solar and wind¹⁻³. These batteries utilize external storage tanks for liquid-phase electrolytes. This design facilitates maintenance and allows for scalable energy accumulation². Traditional RFBs primarily employ metal-based active materials such as vanadium, iron, and zinc⁴. However, concerns over their cost and toxicity have hindered their widespread utilization⁵,⁶. Organic redox-active materials, particularly quinones, offer a promising alternative due to their abundance and highly tunable electrochemical windows¹,²,⁷⁻⁹.

Quinones are particularly noted for their rapid redox kinetics, making them suitable for advanced RFB technologies⁷,¹⁰. During their reduction process, each quinone can accept up to two electrons due to their characteristic two carbonyl groups. The precise mechanism

of this reduction process depends strongly on the nature of the solvent. In organic aprotic solvents, where quinones are highly soluble, the reduction is a sequence of two single-electron transfers¹¹⁻¹⁴. Alternatively, it has been highlighted that quinone derivatives undergo a simultaneous double-electron reduction and double proton transfer in acidic aqueous solutions⁷,¹⁵,¹⁶. These proton-coupled electron transfers (PCET) have encouraged the development of aqueous RFBs (ARFBs). However, practical applications of quinones in ARFBs are often limited by their low solubility and suboptimal half-cell potentials¹⁷,¹⁸.

To increase the redox potential window of quinones, a common strategy is to create new quinone derivatives by functionalizing them with electron-withdrawing or donating groups. Due to gigantic chemical space of possible quinones, experimental screening would be too time-consuming and costly. To circumvent this problem, the idea of high-throughput computational screening has emerged¹,²,¹⁶,¹⁹⁻²³. In most studies, the chemical properties of each compound are computed using electronic Density Functional Theory (DFT), while the solvent is represented with a continuum solvent model, among which the polarizable continuum model (PCMs) is very popular.

<footnote>a) Electronic mail: guillaume.jeanmairet@sorbonne-universite.fr</footnote>

<footnote>b) Electronic mail: emmanuel.giner@lct.jussieu.fr</footnote>
Continuum solvent models allow accounting for solva-
tion effects at a very moderate computational cost but
suffer from several limitations. The major one is that
they replace a series of interactions between molecules
with a macroscopic electrostatic picture. Classical den-
sity functional theory (cDFT) is a promising alternative
approach that can account for solvation effects at a com-
petitive cost while retaining a molecular description.

In short, cDFT allows computing the grand-potential
by minimizing a functional of the liquid density²⁴⁻²⁶.
This minimum is reached for the equilibrium density.
Similarly to its electronic counterpart, cDFT is theoret-
ically well grounded but approximations in the choice of
the functional describing the interactions between par-
ticles, here solvent molecules, are required to make its
utilisation practical. Several attempts to use cDFT as a
solver for QM/MM calculations have been proposed²⁷⁻³⁰.

We recently proposed to use molecular density func-
tional theory (MDFT), a flavour of cDFT that is designed
to study solvation of classical solutes in rigid, classical
solvents to account for solvent effects in QM calculations.
Formally, one can show that there is a variational prin-
ciple for the free energy in terms of the full QM/MM
density matrix. The computation of the free energy can
be reformulated using constrained-search methods à la
Levy-Lieb into a variational problem involving only the
quantum and classical one-body densities³¹. In practice,
the QM/MDFT approach can be used with any QM tech-
niques such as electronic DFT³² or selected CI³³, and has
been used to predict ground- and excited-states proper-
ties of organic molecules, and chemical reactivity.

In this paper, electronic DFT and MDFT are used
to compute the aqueous redox potentials of a series of
quinones. In particular, we focus on the simultane-
ous double proton charge electron transfer reduction.
While our previous applications of the QM/MDFT ap-
proach were limited to small molecules such as water or
formaldehyde, the considered quinones are composed of
several dozens of atoms, and some functional groups are
quite flexible. This makes the choice of the geometry of
solute to be used in the calculation less straightforward.
Thus, we propose, in Sec. II, a theoretical framework to
carry on the geometry optimization of the solute within
the QM/MDFT formalism. This is done by formulat-
ing a variational principle for the grand potential as de-
pending only on the nuclear variables, *i.e.* the nuclear
density matrix. The nuclei, however, experience an effec-
tive potential that is due to the electrons and the solvent
molecules. We then propose a series of practical approx-
imations that allow performing the computation of the
grand potential as a series of quantum and classical den-
sities optimizations followed by geometry optimization.
The harmonic approximation is used to account for the
vibrational motion of the nuclei and requires the com-
putation of the Hessian matrix of the nuclei-nuclei inter-
action energy. In Sec. III, we describe some theoretical
and computational details before focusing on the redox
properties of quinones in Sec. IV. The quality of our re-

sults is assessed by a thorough comparison with previous
predictions obtained by Huynh et al.<sup>16</sup> using eDFT in a
conductor-like PCM (C-PCM) and to experimental data
reported by the same group.

II. THEORY: VARIATIONAL FRAMEWORK FOR QM/MM CALCULATIONS

This section is dedicated to the description of the
grand potential of a QM/MM system where both the
electrons and nuclei of the solute are treated quantum
mechanically, while the solvent is treated at the MM
level. Starting from a very general QM/MM framework
for a nuclei-electron/classical solvent system and thanks
to well-defined approximations, we obtain a formalism
similar to that used in "thermodynamical corrections"
in standard QM calculations using a continuous solvent
model. We therefore derive an approximation for the
grand potential of an electron-nuclei QM system solvated
by classical molecules, based on geometry optimization
and frequency calculations.

We first present in Sec. II A the Hamiltonian of the
QM/MM system and establish the exact equilibrium den-
sity matrix (EDM) in the semi-grand canonical ensemble,
together with the corresponding grand potential. Then,
in Sec. II B we show how the exact EDM can be simpli-
fied within the Born Oppenheimer (BO) and zero elec-
tronic entropy (ZEE) approximations, which correspond
to an adiabatic electronic ground state approximation
while keeping the entropic term for both the QM nu-
clei and classical solvent. Then, inspired by a previous
work³¹, in Sec. II C we propose a variational formulation
of the grand potential in terms of high-dimensional ob-
jects. In Sec. IIC & IID, by using a mean mean-field
approximation, which ignore the explicit dependency of
the quantum objects on the classical degrees of free-
dom, it becomes possible to reuse the recently proposed
eDFT/cDFT mean-field coupling. Eventually, by ap-
proximating the nuclear effective potential as a harmonic
potential, we obtain a practical computational scheme to
approximate the Helmholtz free energy of the QM/MM
system, which is presented in Sec. IIE.

A. Exact QM/MM Equilibrium Density Matrix in the semi-grand canonical ensemble and grand potential

We consider a solute/solvent QM/MM system at
finite temperature *T* in a fixed volume *V*. The
solute is described at the QM level and is con-
stituted of *N*<sub>e</sub> electrons, whose positions are de-
noted by **r**<sup>*N*<sub>e</sub></sup> = {*r*<sub>1</sub>, ..., *r*<sub>*N*<sub>e</sub></sub>} and *N*<sub>N</sub> nuclei identi-
fied by their positions **R**<sup>*N*<sub>N</sub></sup> = {*R*<sub>1</sub>, ..., *R*<sub>*N*<sub>N</sub></sub>}, charges
**Z**<sup>*N*<sub>N</sub></sup> = {*Z*<sub>1</sub>, ..., *Z*<sub>*N*<sub>N</sub></sub>}, and masses {*M*<sub>1</sub>, ..., *M*<sub>*N*<sub>N</sub></sub>}.
On the other hand, the solvent environment is described by
classical rigid molecules, whose number *N*<sub>c</sub> is allowed to
fluctuate due to exchange with a reservoir at fixed chem-
ical potential $\mu$. The appropriate statistical ensemble to describe the current system is thus the semi-grand canonical ensemble (N$_e$, N$_N$, $\mu$, V, T).

For a given number N$_c$ of solvent molecules, a solvent configuration is described by the molecules center of mass $\mathbf{Q}^{N_c} = \{\mathbf{Q}_1, \dots, \mathbf{Q}_{N_c}\}$, and their orientations $\mathbf{\Omega}^{N_c} = \{\mathbf{\Omega}_1, \dots, \mathbf{\Omega}_{N_c}\}$, where $\mathbf{\Omega}_i$ is a shorthand for the three Euler angles ($\theta_i, \phi_i, \Psi_i$) of the i$^{th}$ molecule. For conciseness, these positional and orientational variables are collectively denoted as $\tilde{\mathbf{Q}}^{N_c} = \{\mathbf{Q}^{N_c}, \mathbf{\Omega}^{N_c}\}$. The superscript N$_c$ indicates the number of particles, which, as a reminder, is allowed to fluctuate. Additionally, the solvent molecules possess linear momenta $\mathbf{P}^{N_c} = \{\mathbf{P}_1, \dots, \mathbf{P}_{N_c}\}$. From this point on, atomic units are used.

The total Hamiltonian of this QM/MM system with N$_c$ classical particles is defined as:

$$ \hat{\mathbf{H}}_{\text{tot}}(\tilde{\mathbf{Q}}^{N_c}) = \hat{\mathbf{H}}_{\text{QM}}(\tilde{\mathbf{Q}}^{N_c}) + T_{N_c}(\mathbf{P}^{N_c}) + W_{N_c}(\tilde{\mathbf{Q}}^{N_c}), \quad (1) $$

where $T_{N_c}$ is the kinetic energy of N$_c$ solvent molecules and $W_{N_c}$ is the potential energy of interaction. Note that, since the solvent is constituted of rigid molecules, their angular momenta also contribute to Hamiltonian through a rotational kinetic energy term. However, similarly to the linear momenta, these angular momenta are not coupled to the other degrees of freedom and can therefore be integrated analytically in the classical (grand-)partition function. To avoid unnecessarily cumbersome notation, the rotational kinetic energy term is ignored in Eq. (1) and in the following derivation, without any loss of generality.

In Eq. (1), $\hat{\mathbf{H}}_{\text{QM}}(\tilde{\mathbf{Q}}^{N_c})$ is the quantum Hamiltonian describing the electrons and nuclei. Since it also includes the QM/MM interaction, it depends parametrically on the coordinates $\tilde{\mathbf{Q}}^{N_c}$ of the N$_c$ classical particles. The QM Hamiltonian appearing in Eq. (1) is defined as the sum of nuclear and electronic contributions:

$$ \hat{\mathbf{H}}_{\text{QM}}(\tilde{\mathbf{Q}}^{N_c}) = \hat{T}_N + W_{N}^{\text{MM}}(\tilde{\mathbf{Q}}^{N_c}) + \hat{\mathbf{H}}_{\text{BO}}(\tilde{\mathbf{Q}}^{N_c}, \mathbf{R}^{N_N}), \quad (2) $$

where $\hat{T}_N$ is the kinetic energy operator of the nuclei and $W_{N}^{\text{MM}}$ is their potential energy of interaction with N$_c$ molecules of solvent at configuration $\tilde{\mathbf{Q}}^{N_c}$. The electronic Hamiltonian $\hat{\mathbf{H}}_{\text{BO}}(\tilde{\mathbf{Q}}^{N_c})$ is the sum of the usual electronic Born-Oppenheimer (BO) Hamiltonian and the interaction between the electrons and N$_c$ classical solvent molecules at configuration $\tilde{\mathbf{Q}}^{N_c}$, labelled $W_{e}^{\text{MM}}(\tilde{\mathbf{Q}}^{N_c})$, which therefore reads:

$$ \hat{\mathbf{H}}_{\text{BO}}(\tilde{\mathbf{Q}}^{N_c}, \mathbf{R}^{N_N}) = \hat{T}_e + \hat{W}_{ee} + \hat{W}_{Ne}(\mathbf{R}^{N_N}) + \hat{W}_{NN}(\mathbf{R}^{N_N}) \\ + W_{e}^{\text{MM}}(\tilde{\mathbf{Q}}^{N_c}). \quad (3) $$

In Eq. (3), the first four terms are the usual contributions to the electronic BO Hamiltonian, i.e. the electronic kinetic operator $\hat{T}_e$, mutual Coulomb repulsion $\hat{W}_{ee}$, Coulomb attraction with the nuclei $\hat{W}_{Ne}$, and nuclear-nuclear Coulomb repulsion $\hat{W}_{NN}$, respectively.

The QM/MM terms $W_{N}^{\text{MM}}$ and $W_{e}^{\text{MM}}$ appearing in Eqs. (2) and (3) lack a unique, practical, definition and several propositions for these terms have been made in the literature. Therefore, we will carry the general derivations without explicitly defining $W_{N}^{\text{MM}}$ and $W_{e}^{\text{MM}}$, and simply assume the following generic forms,

$$ W_{N}^{\text{MM}} = \sum_{i=1}^{N_c} \sum_{j=1}^{N_N} w_{N}^{\text{MM}}(\mathbf{R}_j, \tilde{\mathbf{Q}}_i), \quad (4) $$

$$ W_{e}^{\text{MM}} = \sum_{i=1}^{N_c} \sum_{j=1}^{N_e} w_{e}^{\text{MM}}(\mathbf{r}_j, \tilde{\mathbf{Q}}_i), \quad (5) $$

where $w_{N}^{\text{MM}}(\mathbf{R}, \tilde{\mathbf{Q}})$ and $w_{e}^{\text{MM}}(\mathbf{r}, \tilde{\mathbf{Q}})$ are functions describing the pairwise nuclear-solvent and electron-solvent interactions, respectively.

The only assumption is that the QM/MM interaction is pairwise additive, both for electrons and nuclei. We further emphasize that, while we chose to include these coupling terms into the Quantum Hamiltonian, they also act on the classical particles.

Within these definitions, the equilibrium density matrix (EDM) for the QM/MM system with N$_c$ classical solvent molecules in the canonical ensemble can be expressed as a quantum density matrix parametrized by the classical degrees of freedom ($\tilde{\mathbf{Q}}^{N_c}, \mathbf{P}^{N_c}$), namely:

$$ \hat{P}_{\text{eq}}^{N_c}(\tilde{\mathbf{Q}}^{N_c}, \mathbf{P}^{N_c}) = \frac{e^{-\beta(\hat{\mathbf{H}}_{\text{QM}}(\tilde{\mathbf{Q}}^{N_c}) + W_{N_c}(\tilde{\mathbf{Q}}^{N_c}))} e^{-\beta T_{N_c}}(\mathbf{P}^{N_c})}{Z_{N_c}}, \quad (6) $$

where the superscript "N$_c$" emphasizes that it is the canonical EDM when there are "N$_c$" classical particles, and $Z_{N_c}$ is the partition function of the QM/MM system with N$_c$ classical solvent molecules,

$$ Z_{N_c} = \iint \text{Tr}_{\text{QM}}[e^{-\beta \hat{\mathbf{H}}_{\text{QM}}(\tilde{\mathbf{Q}}^{N_c})}] e^{-\beta (W_{N_c}(\tilde{\mathbf{Q}}^{N_c}) + T_{N_c})} d\mathbf{P}^{N_c} d\tilde{\mathbf{Q}}^{N_c}. \quad (7) $$

The "Tr$_{QM}$" operator in Eq. (7) denotes the usual trace over the quantum degrees of freedom. The form of the EDM as in Eq. (6) is intuitive, as it follows the standard $\exp(-\beta H)$ form of statistical physics. It is worth emphasizing that this expression can be derived from first principles using Wigner transformation techniques and semi-classical expansions (see Refs. 31 and 34).

Once the EDM and the partition function have been defined in the canonical ensemble it is rather natural to extend the formalism to the semi-grand canonical ensemble where the number of classical particles, N$_c$, is allowed to fluctuate. Therefore, the corresponding EDM can be expressed as the direct sum over the canonical density matrices for all possible particle number N$_c$, weighted by the Boltzman factor associated to the chemical potential $\mu$. This construction of the semi-grand canonical ensemble is conceptually similar to the definition of the Fock space in quantum mechanics, where it is formulated
as the direct sum of all N-particle Hilbert spaces. The
grand canonical equilibrium density matrix associated to
the number of classical particle N<sub>c</sub> is defined as:

$$
\hat{P}_{\text{eq}}^{\text{GC}}(\tilde{\mathbf{Q}}^{\text{N}_c}, \mathbf{P}^{\text{N}_c}) = \frac{e^{-\beta(\hat{H}_{\text{QM}}(\tilde{\mathbf{Q}}^{\text{N}_c}) + W_{\text{N}_c}(\tilde{\mathbf{Q}}^{\text{N}_c}) - \mu_{\text{N}_c})} e^{-\beta T_{\text{N}_c}(\mathbf{P}^{\text{N}_c})}}{\Xi} \quad (8)
$$

where the grand-partition function Ξ is defined as

$$
\Xi = \sum_{N_c=0}^{\infty} \frac{1}{N_c!} Z_{N_c} e^{\beta \mu_{N_c}}. \quad (9)
$$

In the following, as we will constantly be dealing with
quantities both in the canonical and semi-grand canoni-
cal ensembles, it is useful to introduce some conventions.
When considering a quantity, say *O*, without explicit
mention to the number of particles *N*<sub>c</sub>, it will then be
the semi-grand canonical extension of the operator *O*<sup>*N*<sub>c</sub></sup>,
the latter being defined for each number of classical parti-
cles *N*<sub>c</sub>. As usually done in quantum mechanics with the
notion of Fock space, the semi-grand canonical extension
of an operator *O*<sup>*N*<sub>c</sub></sup> can be formally defined in terms of
a direct sum, over all particle number sectors *N*<sub>c</sub>, of the
canonical operators *O*<sup>*N*<sub>c</sub></sup>. Within the present context,
it is also useful to introduce the semi-grand canonical
classical trace operator, Tr<sub>GC</sub>, of a semi-grand canonical
operator *O* as follows

$$
\mathrm{Tr}_{\mathrm{GC}}[O] = \sum_{N_c=0}^{\infty} \frac{1}{N_c!} \mathrm{Tr}[O^{N_c}], \quad (10)
$$

where Tr[$O^{N_c}$] is the canonical trace over the quantum
and N_c classical particles degrees of freedom. Within
these definitions, the partition function of Eq. (9) be-
comes

$$
\Xi = \mathrm{Tr}_{\mathrm{GC}}[e^{-\beta(\hat{H}_{\mathrm{QM}} + W - \mu N)} e^{-\beta T}], \quad (11)
$$

where $\hat{H}_{QM}$, $W$, $T$ and $N$ are the semi-grand canonical extensions of the quantum Hamiltonian, classical potential energy, classical kinetic energy, and classical number of particles, respectively.

One can then naturally define the grand potential, Ω,
that is the thermodynamic state function of this system

$$
\Omega = -\frac{1}{\beta} \log(\Xi). \qquad (12)
$$

One of the objectives of this paper is to propose a vari-
ational formulation of the grand potential, which ex-
tends our previous work in two significant ways. First,
the classical solvent particles are now described within
a grand-canonical framework, allowing their number to
fluctuate. Second, the nuclei of the solute are no longer
treated as clamped, but instead their degrees of free-
dom are explicitly accounted for in the expression of
the grand-canonical functional. We begin by introduc-
ing two well-established approximations that simplify the
expression of the equilibrium density matrix: the Born-
Oppenheimer approximation and the zero electronic en-
tropy (ZEE) approximation.

B. Born Oppenheimer and zero electronic entropy approximations

For a given number of classical particle N<sub>c</sub>, the grand-canonical equilibrium density matrix P̂<sub>eq</sub><sup>GC</sup> can be expressed as the product of a purely classical probability distribution p<sub>eq</sub><sup>GC</sup>(Q̃<sup>N<sub>c</sub></sup>, P<sup>N<sub>c</sub></sup>) with a purely quantum electronic-nuclear density matrix, p̂<sub>eq</sub>(Q̃<sup>N<sub>c</sub></sup>), parametrically depending on the classical variable:

$$
\hat{P}_{\text{eq}}^{\text{GC}}(\tilde{\mathbf{Q}}^{\text{N}_c}, \mathbf{P}^{\text{N}_c}) = p_{\text{eq}}^{\text{GC}}(\tilde{\mathbf{Q}}^{\text{N}_c}, \mathbf{P}^{\text{N}_c}) \hat{\rho}_{\text{eq}}(\tilde{\mathbf{Q}}^{\text{N}_c}), \quad (13)
$$

with

$$
p_{\text{eq}}^{\text{GC}}(\tilde{\mathbf{Q}}^{\text{N}_c}, \mathbf{P}^{\text{N}_c}) = Z_{\text{QM}}(\tilde{\mathbf{Q}}^{\text{N}_c}) \frac{e^{-\beta(W_{\text{N}_c}(\tilde{\mathbf{Q}}^{\text{N}_c}) + T_{\text{N}_c}(\mathbf{P}^{\text{N}_c}) - \mu_{\text{N}_c})}}{\Xi}, \quad (14)
$$

$$
\hat{\rho}_{\text{eq}}(\tilde{\mathbf{Q}}^{\text{N}_c}) = \frac{e^{-\beta \hat{H}_{\text{QM}}(\tilde{\mathbf{Q}}^{\text{N}_c})}}{Z_{\text{QM}}(\tilde{\mathbf{Q}}^{\text{N}_c})} \quad \text{and} \quad (15)
$$

$$
Z_{\text{QM}}(\tilde{\mathbf{Q}}^{\text{N}_c}) = \text{Tr}_{\text{QM}} \{ e^{-\beta \hat{H}_{\text{QM}}(\tilde{\mathbf{Q}}^{\text{N}_c})} \}. \quad (16)
$$

It is important to note that, although ˆP̂<sup>GC</sup><sub>eq</sub> is written as
a product of a QM density matrix and a classical prob-
ability distribution, it explicitly incorporates QM/MM
correlation effects, as the quantum density matrix de-
pends parametrically on the classical solvent coordinates
Q̃<sup>Nc</sup>. Furthermore, the classical probability distribution
is influenced by the quantum subsystem through the nor-
malization condition (i.e. through Z<sub>QM</sub>(Q̃<sup>Nc</sup>)), thus pre-
serving the coupling between the two subsystems.

To simplify the quantum density matrix ˆρ<sub>eq</sub>, several
standard approximations can be employed. First, mak-
ing the BO approximation allows for rewriting ˆρ<sub>eq</sub> as a
weighted sum over electronic-nuclear adiabatic states, la-
belled by indices (α, n), where α denotes the electronic
state and n the nuclear vibrational state, i.e.

$$
\hat{\rho}_{\text{eq}}(\tilde{\mathbf{Q}}^{\text{N}_c}) = \sum_{\alpha, n} \frac{e^{-\beta \mathcal{E}_{n}^{\alpha}(\tilde{\mathbf{Q}}^{\text{N}_c})}}{Z_{\text{QM}}(\tilde{\mathbf{Q}}^{\text{N}_c})} |\Phi_{n}^{\alpha}(\tilde{\mathbf{Q}}^{\text{N}_c})\rangle \langle \Phi_{n}^{\alpha}(\tilde{\mathbf{Q}}^{\text{N}_c})| \quad (17)
$$

In the above expression, $\mathcal{E}_n^\alpha(\tilde{\mathbf{Q}}^{\text{N}_c})$ and $|\Phi_n^\alpha(\tilde{\mathbf{Q}}^{\text{N}_c})\rangle$ are the eigenvalues and eigenvectors of the total electronic-nuclear adiabatic quantum Hamiltonian (i.e. neglecting the non-adiabatic couplings between electrons and nuclei) for a given solvent configuration $\tilde{\mathbf{Q}}^{\text{N}_c}$. This expression can be further simplified in the case of a weakly correlated system at room temperature, where one can neglect the contribution from electronically excited adiabatic states. Under this assumption, only the ground electronic state $\alpha = 0$ is retained in the summation of Eq. (17). This approximation, hereafter referred to as zero electronic entropy (ZEE), is justified by the fact that
the relative statistical weights of excited adiabatic electronic states decay exponentially with increasing energy gaps. Consequently, their contributions become negligible at room temperature.

Under the ZEE approximation, and exploiting the product structure of the electron-nuclei wave function provided by the BO approximation, the quantum electronic-nuclear density matrix can be written as the product of an electronic density matrix and a nuclear density matrix

$$ \hat{\rho}_{\text{eq}}(\tilde{\mathbf{Q}}^{\text{N}_c}) = \hat{\rho}_{\text{N}}^{0}(\tilde{\mathbf{Q}}^{\text{N}_c}) \hat{\rho}_{\text{el}}^{0}(\tilde{\mathbf{Q}}^{\text{N}_c}, \mathbf{R}^{\text{N}_N}), \quad (18) $$

where $\hat{\rho}_{\text{el}}^{0}(\tilde{\mathbf{Q}}^{\text{N}_c})$ is the ground state electronic adiabatic density matrix depending parametrically on both the solvent and nuclear coordinates

$$ \hat{\rho}_{\text{el}}^{0}(\tilde{\mathbf{Q}}^{\text{N}_c}, \mathbf{R}^{\text{N}_N}) = \left| \Psi^{0}(\tilde{\mathbf{Q}}^{\text{N}_c}, \mathbf{R}^{\text{N}_N}) \right\rangle \left\langle \Psi^{0}(\tilde{\mathbf{Q}}^{\text{N}_c}, \mathbf{R}^{\text{N}_N}) \right|, \quad (19) $$

and $\hat{\rho}_{\text{N}}^{0}(\tilde{\mathbf{Q}}^{\text{N}_c})$ is the nuclear density matrix which contains contributions from all excited nuclear state associated to the adiabatic electronic ground state

$$ \hat{\rho}_{\text{N}}^{0}(\tilde{\mathbf{Q}}^{\text{N}_c}) = \sum_{n} \frac{e^{-\beta \varepsilon_n^0(\tilde{\mathbf{Q}}^{\text{N}_c})}}{Z_{\text{QM}}(\tilde{\mathbf{Q}}^{\text{N}_c})} \left| \chi_{\text{N}}^{0,n}(\mathbf{R}^{\text{N}_n}; \tilde{\mathbf{Q}}^{\text{N}_c}) \right\rangle \left\langle \chi_{\text{N}}^{0,n}(\mathbf{R}^{\text{N}_n}; \tilde{\mathbf{Q}}^{\text{N}_c}) \right|, \quad (20) $$

and where the partition function defined in Eq. (16) simplifies into a purely nuclear contribution

$$ Z_{\text{QM}}(\tilde{\mathbf{Q}}^{\text{N}_c}) = \sum_n e^{-\beta \varepsilon_n^0(\tilde{\mathbf{Q}}^{\text{N}_c})}. \quad (21) $$

### C. QM/MM Variational Formulation and mean-field approximation

In the previous sections, we derived an explicit expression for the equilibrium density matrix of a QM/MM system made of $N_e$ electrons and $N_N$ nuclei, and a fluctuating number of classical particles. This derivation was performed within the semi-grand canonical ensemble, employing both the BO and ZEE approximations. We now propose a variational formulation of the grand potential for this system, extending the approach introduced in our previous work within the canonical ensemble$^{31}$.

More specifically, we consider a general QM/MM electron-nuclei-solvent density matrix, $\hat{\mathcal{P}}^{\text{GC}}$, which can differ from the equilibrium one, while conserving the product structure of Eq. (13). Therefore, for a given number of classical particles $N_c$, it has the following structure

$$ \hat{\mathcal{P}}^{\text{GC}}(\tilde{\mathbf{Q}}^{\text{N}_c}, \mathbf{P}^{\text{N}_c}) = p^{\text{GC}}(\tilde{\mathbf{Q}}^{\text{N}_c}, \mathbf{P}^{\text{N}_c}) \hat{\rho}(\tilde{\mathbf{Q}}^{\text{N}_c}). \quad (22) $$

In Eq. (22), $p^{\text{GC}}$ is the probability of having $N_c$ classical particles with positions and momenta $(\tilde{\mathbf{Q}}^{\text{N}_c}, \mathbf{P}^{\text{N}_c})$,

and $\hat{\rho} = \hat{\rho}_{\text{N}}\hat{\rho}_{\text{el}}$ is the electron-nuclei density matrix, which has a similar form to Eq. (18), and therefore depends parametrically on the classical solvent configuration $\tilde{\mathbf{Q}}^{\text{N}_c}$. Therefore, $\hat{\mathcal{P}}^{\text{GC}}$ is normalized to unity, i.e. $\text{Tr}_{\text{GC}}[\hat{\mathcal{P}}^{\text{GC}}] = 1$. We can now introduce the grand potential functional for the QM/MM system in the semi-grand canonical ensemble

$$ \omega[\hat{\mathcal{P}}^{\text{GC}}] = \mathcal{U}_{\text{int}}[\hat{\mathcal{P}}^{\text{GC}}] - \frac{1}{\beta} \mathcal{S}_{\text{int}}[\hat{\mathcal{P}}^{\text{GC}}], \quad (23) $$

where $\mathcal{U}_{\text{int}}[\hat{\mathcal{P}}^{\text{GC}}]$ is the internal energy of the QM/MM system

$$ \mathcal{U}_{\text{int}}[\hat{\mathcal{P}}^{\text{N}_c}] = \text{Tr}_{\text{GC}}\{p^{\text{GC}}\hat{\rho}[\hat{\mathbf{H}}_{\text{QM}} + W + T - \mu N]\}, \quad (24) $$

and $\mathcal{S}_{\text{int}}[\hat{\mathcal{P}}^{\text{GC}}]$ is the dimensionless entropy associated with the grand-canonical density matrix $\hat{\mathcal{P}}^{\text{GC}}$

$$ \mathcal{S}_{\text{int}}[\hat{\mathcal{P}}^{\text{GC}}] = -\text{Tr}_{\text{GC}}\{p^{\text{GC}}\hat{\rho}_{\text{N}} \log(p^{\text{GC}}\hat{\rho}_{\text{N}})\}. \quad (25) $$

Notice that the definition of the entropy in Eq. (25) does not contain the electronic density matrix, which is a consequence of the ZEE approximation.

Following the same steps as in Ref. 31 (i.e. using Jensen's inequality), it can be shown that the functional $\omega$ of Eq. (23) satisfies the variational principle for the grand potential, i.e.

$$ \Omega = \min_{\hat{\mathcal{P}}^{\text{GC}}} \omega[\hat{\mathcal{P}}^{\text{GC}}], \quad \hat{\mathcal{P}}_{\text{eq}}^{\text{GC}} = \operatorname*{argmin}_{\hat{\mathcal{P}}^{\text{GC}}} \omega[\hat{\mathcal{P}}^{\text{GC}}], \quad (26) $$

where the minimization is performed over all admissible density matrices $\hat{\mathcal{P}}_{\text{eq}}^{\text{GC}}$ consistent with the assumed product structure and normalization constraints of Eq. (22). The novelty of the present approach lies in formulating the variational principle within the semi-grand canonical ensemble, and explicitly accounting for the quantum degrees of freedom of the solute nuclei.

Although the variational form of Eq. (26) remains highly complex due to the involvement of high-dimensional objects, it provides a foundation for introducing practical approximations that eventually yield to a one-body density-functional formulation. The first step in this strategy is to neglect the explicit electron-solvent QM/MM correlation terms as proposed in earlier works$^{32,33}$, which corresponds to adopting a mean-field approximation.

Therefore, we propose to extend this mean-field approximation to the electron-nuclei-solvent EDM $\hat{\mathcal{P}}^{\text{GC}}$ in Eq. (22), by ignoring the explicit dependence on $\tilde{\mathbf{Q}}^{\text{N}_c}$ in the QM part. This leads to the following mean-field canonical QM/MM electron-nuclei-solvent density matrix in the semi-grand-canonical ensemble

$$ \hat{\mathcal{P}}_{\text{MF}}^{\text{GC}}(\tilde{\mathbf{Q}}^{\text{N}_c}, \mathbf{P}^{\text{N}_c}) = p^{\text{GC}}(\tilde{\mathbf{Q}}^{\text{N}_c}, \mathbf{P}^{\text{N}_c}) \hat{\rho}_{\text{N}} \hat{\rho}_{\text{el}}(\mathbf{R}^{\text{N}_N}). \quad (27) $$

Here, the density matrices $\hat{\rho}_{\text{N}}$ and $\hat{\rho}_{\text{el}}(\mathbf{R}^{\text{N}_N})$ are quantum objects that do not depend explicitly on the classical degrees of freedom, and therefore neglect explicit QM/MM
correlation effects. Inserting this form of $\hat{\mathcal{P}}^{\text{GC}}$ into the variational principle of Eq. (26) yields the mean-field approximation $\Omega_{\text{MF}}$ which provides an upper-bound to the exact grand potential $\Omega$, i.e.

$$ \Omega_{\text{MF}} = \min_{\hat{\mathcal{P}}_{\text{MF}}^{\text{GC}}} \omega[\hat{\mathcal{P}}_{\text{MF}}^{\text{GC}}], \quad \Omega \le \Omega_{\text{MF}}. \qquad (28) $$

### D. Effective potential for the nuclei and QM/MM density functional formalism

The variational formulation of Eq. (28), expressed in terms of electron-nuclei-solvent objects and the product structure of $\hat{\mathcal{P}}_{\text{MF}}$ allows, after some algebraic manipulations which are explicated in Appendix 1, to approximate $\Omega_{\text{MF}}$ by a variational formulation involving only the nuclear density matrix. Nevertheless, the nuclei now experience an effective potential arising from both the electrons and the solvent. More precisely, $\Omega_{\text{MF}}$ can be rewritten as follows

$$ \Omega_{\text{MF}} = \min_{\hat{\mathcal{P}}_{\text{N}}} (\mathcal{U}_{\text{N}}[\hat{\mathcal{P}}_{\text{N}}] - \frac{1}{\beta} \mathcal{S}_{\text{N}}[\hat{\mathcal{P}}_{\text{N}}]), \qquad (29) $$

where $\mathcal{S}_{\text{N}}[\hat{\mathcal{P}}_{\text{N}}]$ is the purely nuclear quantum entropy

$$ \mathcal{S}_{\text{N}}[\hat{\mathcal{P}}_{\text{N}}] = -\text{Tr}\{\hat{\mathcal{P}}_{\text{N}} \log(\hat{\mathcal{P}}_{\text{N}})\}, \qquad (30) $$

and $\mathcal{U}_{\text{N}}$ is the nuclear quantum energy

$$ \mathcal{U}_{\text{N}}[\hat{\mathcal{P}}_{\text{N}}] = \text{Tr}\{\hat{\mathcal{P}}_{\text{N}}[\hat{T}_{\text{N}} + \hat{W}_{\text{NN}}^{\text{tot}}(\mathbf{R}^{\text{NN}})]\}. \qquad (31) $$

The "Tr" symbol in Eqs. (30)-(31) refers to the trace over the nuclear degrees of freedom. In Eq. (31), $\hat{W}_{\text{NN}}^{\text{tot}}(\mathbf{R}^{\text{NN}})$ is the total effective potential felt by the nuclei defined as

$$ \hat{W}_{\text{NN}}^{\text{tot}}(\mathbf{R}^{\text{NN}}) = \hat{W}_{\text{NN}}(\mathbf{R}^{\text{NN}}) + \mathcal{W}_{\text{N}}^{\text{eff}}(\mathbf{R}^{\text{NN}}), \qquad (32) $$

which is the sum of the nucleus-nucleus Coulomb interaction and of $\mathcal{W}_{\text{N}}^{\text{eff}}$, the effective potential experienced by the nuclei due to the presence of the electrons and solvent molecules. The latter is defined as

$$ \mathcal{W}_{\text{N}}^{\text{eff}}(\mathbf{R}^{\text{NN}}) = \min_{\hat{\mathcal{P}}_{\text{el}}, p^{\text{GC}}} (f_{\text{el}}[\hat{\mathcal{P}}_{\text{el}}] + f_{\text{MM}}^{\text{GC}}[p^{\text{GC}}] + E_{\text{QM}}^{\text{MM}}[\hat{\mathcal{P}}_{\text{el}}, p^{\text{GC}}] + V_{\text{el}}^{\text{N}}[\hat{\mathcal{P}}_{\text{el}}; \mathbf{R}^{\text{NN}}] + V_{\text{MM}}^{\text{N}}[p^{\text{GC}}; \mathbf{R}^{\text{NN}}]), \qquad (33) $$

with $f_{\text{el}}[\hat{\mathcal{P}}_{\text{el}}]$ and $f_{\text{MM}}[p^{\text{GC}}]$ being the intrinsic electronic and classical energy functionals, respectively, defined as

$$ f_{\text{el}}[\hat{\mathcal{P}}_{\text{el}}] = \text{Tr}_{\text{QM}}\{\hat{\mathcal{P}}_{\text{el}}(\hat{T}_{\text{e}} + \hat{W}_{\text{ee}})\}, \qquad (34) $$

$$ \begin{aligned} f_{\text{MM}}^{\text{GC}}[p^{\text{GC}}] &= \text{Tr}_{\text{GC}}[p^{\text{GC}}(W + T - \mu N + \frac{1}{\beta} \log(p^{\text{GC}}))] \\ &= \sum_{N_c=0}^{\infty} \iint p^{\text{GC}}(W_{N_c} + T_{N_c} - \mu N_c + \frac{1}{\beta} \log(p^{\text{GC}})) d\tilde{\mathbf{Q}}^{N_c} d\mathbf{P}^{N_c}. \end{aligned} \qquad (35) $$

On the other hand, $E_{\text{QM}}^{\text{MM}}[\hat{\mathcal{P}}_{\text{el}}, p^{\text{GC}}]$ is the electron-solvent QM/MM interaction

$$ \begin{aligned} E_{\text{QM}}^{\text{MM}}[\hat{\mathcal{P}}_{\text{el}}, p^{\text{GC}}] &= \text{Tr}_{\text{GC}}[p^{\text{GC}} \text{Tr}_{\text{QM}}[\hat{\mathcal{P}}_{\text{el}} \mathcal{W}_{\text{e}}^{\text{MM}}]] \\ &= \sum_{N_c=0}^{\infty} \iint p^{\text{GC}}(\tilde{\mathbf{Q}}^{N_c}, \mathbf{P}^{N_c}) \text{Tr}_{\text{QM}}[\hat{\mathcal{P}}_{\text{el}} \mathcal{W}_{\text{e}}^{\text{MM}}] d\tilde{\mathbf{Q}}^{N_c} d\mathbf{P}^{N_c}, \end{aligned} \qquad (36) $$

$V_{\text{el}}^{\text{N}}[\hat{\mathcal{P}}_{\text{el}}; \mathbf{R}^{\text{NN}}]$ and $V_{\text{MM}}^{\text{N}}[p^{\text{GC}}; \mathbf{R}^{\text{NN}}]$ are the electron and classical interactions with the nuclei, respectively, defined as

$$ V_{\text{el}}^{\text{N}}[\hat{\mathcal{P}}_{\text{el}}; \mathbf{R}^{\text{NN}}] = \text{Tr}_{\text{QM}}[\hat{\mathcal{P}}_{\text{el}} \hat{W}_{\text{Ne}}(\mathbf{R}^{\text{NN}})], \qquad (37) $$

$$ \begin{aligned} V_{\text{MM}}^{\text{N}}[p^{\text{GC}}, \mathbf{R}^{\text{NN}}] &= \text{Tr}_{\text{GC}}[p^{\text{GC}} \mathcal{W}_{\text{N}}^{\text{MM}}] \\ &= \sum_{N_c=0}^{\infty} \iint p^{\text{GC}}(\tilde{\mathbf{Q}}^{N_c}, \mathbf{P}^{N_c}) \mathcal{W}_{\text{N}}^{\text{MM}}(\tilde{\mathbf{Q}}^{N_c}, \mathbf{R}^{\text{NN}}) d\tilde{\mathbf{Q}}^{N_c} d\mathbf{P}^{N_c}. \end{aligned} \qquad (38) $$

An important remark is that both the electron-nuclear and the solvent-nuclear potentials are expressed as a sum of pairwise terms (cf. Eq. (4) for the latter). Therefore, for a given nuclear coordinate $\mathbf{R}^{\text{NN}}$, $V_{\text{el}}^{\text{N}}$ and $V_{\text{MM}}^{\text{N}}$ act on the electrons and solvent molecules as external local potentials expressed as sums of one-body terms. Consequently, $\mathcal{W}_{\text{N}}^{\text{eff}}(\mathbf{R}^{\text{NN}})$ is obtained as a minimization of a functional describing a QM/MM system made of interacting electrons and classical particles, both embedded in local external potentials. We can therefore identify $\mathcal{W}_{\text{N}}^{\text{eff}}(\mathbf{R}^{\text{NN}})$ as the mean-field approximation of the grand potential for the electron-solvent system which was formalized in Ref. 31. As shown in the latter, the Levy-Lieb constrained search formalism$^{35-37}$ can be used to express $\mathcal{W}_{\text{N}}^{\text{eff}}(\mathbf{R}^{\text{NN}})$ as a functional of the one-body electronic and classical densities. More precisely, $\mathcal{W}_{\text{N}}^{\text{eff}}(\mathbf{R}^{\text{NN}})$ can be written as

$$ \mathcal{W}_{\text{N}}^{\text{eff}}(\mathbf{R}^{\text{NN}}) = \min_{\rho, n} w_{\text{N}}^{\text{eff}}[\rho, n, \mathbf{R}^{\text{NN}}], \qquad (39) $$

where, $\rho$ and $n$ are respectively the one-body electronic and classical densities. In (39), $w_{\text{N}}^{\text{eff}}[\rho, n, \mathbf{R}^{\text{NN}}]$ is the one-body density grand potential functional defined as

$$ w_{\text{N}}^{\text{eff}}[\rho, n, \mathbf{R}^{\text{NN}}] = \mathcal{F}_{\text{el}}[\rho] + (v_{\text{Ne}}(\mathbf{R}^{\text{NN}})|\rho) + \mathcal{F}_{\text{MM}}[n] \\ + (v_{\text{N}}^{\text{MM}}(\mathbf{R}^{\text{NN}})|n) + (n|v_{\text{e}}^{\text{MM}}|\rho). \qquad (40) $$

In Eq. (40), $\mathcal{F}_{\text{el}}[\rho]$ is the usual electronic Levy-Lieb$^{35-37}$ universal functional, $(v_{\text{Ne}}(\mathbf{R}^{\text{NN}})|\rho)$ is the one-body nuclear-electron interaction, $\mathcal{F}_{\text{MM}}[n]$ is the classical universal functional, as defined by Evans$^{25}$, $(v_{\text{N}}^{\text{MM}}(\mathbf{R}^{\text{NN}})|n)$ is the one-body nuclear-solvent interaction defined as

$$ (v_{\text{N}}^{\text{MM}}(\mathbf{R}^{\text{NN}})|n) = \int d\tilde{\mathbf{Q}} v_{\text{N}}^{\text{MM}}(\tilde{\mathbf{Q}}; \mathbf{R}^{\text{NN}}) n(\tilde{\mathbf{Q}}), \qquad (41) $$

with $v_{\text{N}}^{\text{MM}}(\tilde{\mathbf{Q}}; \mathbf{R}^{\text{NN}}) = \sum_{i=1}^{N_N} w_{\text{N}}^{\text{MM}}(\mathbf{R}_i, \tilde{\mathbf{Q}})$, and $(n|v_{\text{e}}^{\text{MM}}|\rho)$ is the mean-field electron-solvent interaction

$$ (n|v_{\text{e}}^{\text{MM}}|\rho) = \int d\mathbf{r} d\tilde{\mathbf{Q}} \rho(\mathbf{r}) v_{\text{el}}^{\text{MM}}(\mathbf{r}, \tilde{\mathbf{Q}}) n(\tilde{\mathbf{Q}}). \qquad (42) $$
Therefore, by obtaining $W_N^{\text{eff}}(\mathbf{R}^{\text{NN}})$ through a coupled functional minimization over electrons and solvent one-body densities, we can rigorously approximate the grand potential of a QM/MM system composed of quantum electrons and nuclei solely as a functional of the nuclear coordinates. To the best of our knowledge, this represents the first rigorous derivation of such a grand potential formulation explicitly accounting for quantum nuclei within this QM/MM framework.

### E. Harmonic approximation

Although the formulation of Eq. (29) significantly reduces the dimensionality of the problem, it still requires solving a non factorisable quantum problem involving all nuclear degrees of freedom, which is rapidly cumbersome. To address this, and following the treatment proposed by McQuarrie$^{38}$, we approximate the Gibbs free energy of the nuclei by invoking a local description of the external potential and employing a tractable, integrable model. These two criteria naturally lead to the harmonic oscillator (HO) approximation for the effective potential $\hat{W}_{\text{NN}}^{\text{tot}}(\mathbf{R}^{\text{NN}})$ of Eq. (32), along with the rigid-rotor harmonic oscillator (RRHO) model to describe internal rotations of the solute. Several variants of the HO and RRHO approximations have been developed to estimate the Gibbs free energy$^{39-41}$, particularly aimed at mitigating the divergence of entropic contributions in the low-frequency regime. A detailed discussion about the implication of the latter can be found in a recent paper$^{42}$. However, since the primary objective of the present work is to benchmark our solvent model against polarizable continuum models (PCM), we adopt the simplest HO formulation—without low-frequency cutoffs—to remain consistent with widely used implementations in standard electronic structure packages such as Gaussian$^{43,44}$, PySCF$^{45}$, and Orca$^{46}$. Naturally, this choice does not prevent the use of more sophisticated HO or RRHO treatments, which can be incorporated as post-processing steps based on the computed normal modes and their associated vibrational frequencies.

After a series of mathematical derivation detailed in appendix 2, the grand potential can be approximated as

$$ \Omega = \hat{W}_{\text{NN}}^{\text{tot}}(\mathbf{R}_{\text{eq}}^{\text{NN}}) + \sum_{k} \left( \frac{\omega_k}{2} + \frac{1}{\beta} \log(1 - e^{-\beta \omega_k}) \right), \quad (43) $$

where $\hat{W}_{\text{NN}}^{\text{tot}}$ is the effective potential of Eq (32), $\mathbf{R}_{\text{eq}}^{\text{NN}}$ is the equilibrium geometry and $\omega_k$ is the frequency associated to the normal mode $k$.

Evaluating $\Omega$ according to Eq. (43) requires the knowledge of the equilibrium geometry $\mathbf{R}_{\text{eq}}^{\text{NN}}$. Moreover, it is necessary to compute the Hessian matrix of $\hat{W}_{\text{NN}}^{\text{tot}}$ at the equilibrium geometry to compute the normal modes $k$ and their frequencies $\omega_k$. A discussion about the evaluation of the gradient of the total effective potential $\hat{W}_{\text{NN}}^{\text{tot}}$

with respect to the solute nuclear coordinates can be found in Appendix 3.

### III. METHODOLOGY

#### A. Theoretical models

As a practical illustration of the proposed framework, we focus on the computation of redox properties of a series of quinones in aqueous solution. To evaluate the grand potential according to Eq. (43) we must specify how the various terms contributing to the effective potential felt by the nuclei as defined in Eq. (40) are modeled. In practice, the electronic Levy-Lieb functional $\mathcal{F}_{\text{el}}$ is replaced by an approximated electronic functional within the Kohn-Sham formalism. For the classical universal functional $\mathcal{F}_{\text{MM}}$, we chose a MDFT functional which is well suited for the study of solvated systems. Since detailed descriptions of MDFT can be found in our previous work, we only recall here its essential aspects. As mentioned earlier, the classical solvent is constituted of rigid molecules and is described by the solvent density $n(\mathbf{Q}) = n(\mathbf{r}, \mathbf{\Omega})$. The functional of this solvent density is split into two contributions

$$ \mathcal{F}_{\text{MM}}[n] = \mathcal{F}_{\text{id}}[n] + \mathcal{F}_{\text{exc}}[n]. \quad (44) $$

The ideal term $\mathcal{F}_{\text{id}}[n]$ is known exactly and corresponds to the grand potential of a non interacting fluid of density $n$. While no exact practical expression exists for $\mathcal{F}_{\text{exc}}[n]$, it can be approximated by performing a Taylor expansion around a homogeneous reference fluid. This leads to:

$$ \mathcal{F}_{\text{exc}}[n] = -\frac{1}{2\beta} \iint \Delta n(\mathbf{r}) c(|\mathbf{r} - \mathbf{r}'|, \mathbf{\Omega}, \mathbf{\Omega}'; n_0) \Delta n(\mathbf{r}') d\mathbf{r} d\mathbf{r}' d\mathbf{\Omega} d\mathbf{\Omega}' + \mathcal{F}_{\text{B}}[n]. \quad (45) $$

In Eq. (45) the second order term in density involves the second order direct correlation function of the homogeneous fluid at density $n_0$ and $\Delta n = n - n_0$. All the terms of higher order are gathered into the so-called bridge functional $\mathcal{F}_{\text{B}}$. While the second order term can be computed thanks to ingenious computational tricks detailed elsewhere$^{47}$, the higher-order terms are not tractable in practice and require approximations for $\mathcal{F}_{\text{B}}$. The last parts of Eq. (40) that needs to be specified are the QM/MM coupling terms, $v_e^{\text{MM}}$ and $v_N^{\text{MM}}$. Again there is no unique practical definition for these terms. In this work, we adopt a mean-field Coulombic interaction between the solvent and both the electrons and nuclei, supplemented by a Lennard-Jones potential to account for the Pauli repulsion and dispersion forces. More specifically, the electron-solvent interaction term of Eq. (42) is expressed as an electrostatic term

$$ (n|v_e^{\text{MM}}|\rho) = - \int d\mathbf{r}_1 d\mathbf{r}_2 n_c(\mathbf{r}_1) \frac{1}{|\mathbf{r}_1 - \mathbf{r}_2|} \rho(\mathbf{r}_2), \quad (46) $$
where $n_c(\mathbf{r})$ is the classical charge density associated to the classical density $n$, which is given by

$$n_c(\mathbf{r}) = \iint n(\mathbf{r}', \mathbf{\Omega})\gamma(\mathbf{r} - \mathbf{r}', \mathbf{\Omega})d\mathbf{r}'d\mathbf{\Omega}, \quad (47)$$

where $\gamma(\mathbf{r}, \mathbf{\Omega})$ is the point charge distribution of a water molecule, taken at the origin with orientation $\mathbf{\Omega}$.

The nuclear-solvent interaction term, $v_{\text{N}}^{\text{MM}}$, is here expressed as a sum of an electrostatic and Lennard-Jones (LJ) potential contributions, i.e.

$$(v_{\text{N}}^{\text{MM}}(\mathbf{R}^{\text{N}})|n) = (w_{\text{C}}^{\text{MM}}(\mathbf{R}^{\text{N}})|n_c) + (w_{\text{LJ}}^{\text{MM}}(\mathbf{R}^{\text{N}})|n), \quad (48)$$

with

$$(w_{\text{C}}^{\text{MM}}(\mathbf{R}^{\text{N}})|n_c) = \int d\mathbf{r} n_c(\mathbf{r}) \sum_A \frac{Z_A}{|\mathbf{R}_A - \mathbf{r}|}, \quad (49)$$

$$(w_{\text{LJ}}^{\text{MM}}(\mathbf{R}^{\text{N}})|n) = \int d\tilde{\mathbf{Q}} n(\tilde{\mathbf{Q}}) \sum_A \hat{W}_{A}^{\text{LJ}}(\mathbf{R}_A, \tilde{\mathbf{Q}}), \quad (50)$$

where $\hat{W}_{A}^{\text{LJ}}$ is the usual LJ interaction between a nucleus A and a classical water molecule.

Regarding the practical implementation of our QM/MM scheme, it is useful to introduce the following potentials

$$v_{\text{ext}}^{\text{QM}}(\mathbf{r}) = \int d\mathbf{r}' n_c(\mathbf{r}') \frac{1}{|\mathbf{r}' - \mathbf{r}|}, \quad (51)$$

$$v_{\text{ext}}^{\text{MM}}(\mathbf{r}) = \sum_A \frac{Z_A}{|\mathbf{R}_A - \mathbf{r}|} - \int d\mathbf{r}' \rho(\mathbf{r}') \frac{1}{|\mathbf{r}' - \mathbf{r}|}. \quad (52)$$

The potential $v_{\text{ext}}^{\text{QM}}$ is the electrostatic potential created by the MM region and felt by the QM region, and on the other hand $v_{\text{ext}}^{\text{MM}}$ is the electrostatic potential created by the solute (i.e. electrons and nuclei) and felt by the MM region. With these definitions, the electron-solvent interaction of Eq. (46) can then be interpreted as an external electrostatic potential for the QM system, namely

$$(n|v_{\text{e}}^{\text{MM}}|\rho) = - \int d\mathbf{r} v_{\text{ext}}^{\text{QM}}(\mathbf{r}) \rho(\mathbf{r}), \quad (53)$$

and similarly the total electrostatic interaction between the QM and MM region can be obtained as

$$(w_{\text{C}}^{\text{MM}}(\mathbf{R}^{\text{N}})|n_c) + (n|v_{\text{e}}^{\text{MM}}|\rho) = \int d\mathbf{r} v_{\text{ext}}^{\text{MM}}(\mathbf{r}) n_c(\mathbf{r}). \quad (54)$$

The potentials $v_{\text{ext}}^{\text{QM}}$ and $v_{\text{ext}}^{\text{MM}}$ will be used in practice to couple the QM and MM calculations as explained in Sec. III B.

## B. Computational details

This paper focuses on the aqueous solvation of a selection of quinones. The geometries of the 24 quinones displayed in Fig. 1 were extracted from Huynh et al.¹⁶.

For a given nuclear geometry $\mathbf{R}^{\text{N}}$, the methodology of QM/MDFT calculations has already been reported in our previous work³²,³³. In brief, the total effective nuclear potential $\hat{W}_{\text{NN}}^{\text{tot}}(\mathbf{R}^{\text{N}})$ of Eqs. (32) and (39) is obtained from successive QM and MM functional optimizations of both the classical and quantum densities. Each of the QM and MM calculations are coupled to one-another through an effective potential (either $v_{\text{ext}}^{\text{QM}}$ or $v_{\text{ext}}^{\text{MM}}$ of Eqs. (51) and (52)), and the cycle is stopped when some energetic convergence criteria (detailed below) is reached. Then, once the optimal QM and MM densities are obtained, the nuclear gradient is computed to optimize the nuclear geometry of the solute. The solvent contribution to the nuclear gradient (see Sec. 3 a) is then added to the usual electronic gradient contribution, forming the total nuclear gradient which is then used to optimize the geometry. The latter is summarized in Fig. 2.

In more details, the QM calculations are performed using the version 2.6.2 of the PySCF software⁴⁵, employing the B3LYP functional⁴⁸,⁴⁹ and the 6-31++G** basis set³⁰⁻⁵⁴ to approximate the electronic Levy-Lieb functional. This choice of exchange/correlation functional and basis set is similar to that of Huynh et al.¹⁶, such that we can directly compare with the results reported in that previous work. Nevertheless, the QM calculations are influenced by the classical solvent through an additional external electrostatic potential $v_{\text{ext}}^{\text{QM}}(\mathbf{r})$ coming from the classical charge density (see Eq. (51)). The latter is obtained by a set of point-charges lying at the nodes of a 48³ a.u³ cubic box, with 100 grid nodes in each direction. Once the ground state electronic density is obtained, the corresponding electrostatic potential $v_{\text{ext}}^{\text{MM}}(\mathbf{r})$ (see Eq. (52)) generated by the solute is computed on the same cubic grid. This potential is then used as an external potential in the subsequent MDFT calculation.

The solvent density is optimized and computed using our in-house Fortran-based MDFT program. The spatial grid used for the solvent density matches that of the electrostatic potential, and 196 discrete orientations are used to sample the three Euler angles $(\theta, \phi, \psi)$. Water is modeled using the non-polarizable SPC/E model⁵⁵, which then defines the charge density $\gamma(\mathbf{r}, \mathbf{\Omega})$ used in Eq. (47). The $c^{(2)}$ direct correlation function, required to compute the second order term of the excess functional of Eq. (45), was kindly provided by Luc Belloni. It was computed for a homogeneous water density of $n_0 = 0.0332891 \text{Å}^{-3}$ using a combination of Monte Carlo simulation and integral equation theory⁵⁶.

The bridge functional was neglected, i.e. truncating the excess functional at second order, corresponding to the so-called HNC excess functional⁴⁷.

To model dispersion-repulsion interactions and prevent unphysical overlaps between the classical and quan-
tum densities—which could lead to numerical instabilities—LJ sites were added to the solute nuclei following Eq. (50). The LJ parameters are extracted from the second-generation Generalized Amber Force Field (GAFF2)⁵⁷. All selected LJ parameters are provided in SI.

While the gradients of $w_{C}^{MM}$ and $v_{N}^{MM}$ are directly computed in Pyscf, the gradient of the LJ potential $w_{LJ}^{MM}$ is computed in the MDFT program and then added to Pyscf gradient, before performing geometry optimization (see Sec. 3 a for a detailed derivation). A new geometry is proposed using the geomTRIC solver⁵⁸ of PySCF and a new QM/MDFT cycle is performed using this updated structure. This process is iterated until the solute geometry reaches convergence. The full QM/MM calculation is considered complete when the convergence criterion for the solute free energy is satisfied. The threshold used here is $10^{-4}$ Hartree.

Then, the thermodynamic corrections, including zero-point energy contributions, thermochemical corrections at T = 298.15 K and P = 1 atm are computed and added to the solvation free energy and the solute energy yielding the total free energy of the system, which eventually yields to $\Omega$ of Eq. (43). When the full QM/MM calculation is converged, the QM energy, the QM/MM coupling and the classical energy are added to form $\dot{W}_{NN}^{tot}(\mathbf{R}^{NN})$ of Eq. (32).

For the reference quinone, the eDFT/PCM calculations took roughly 5 min while the full eDFT/MDFT calculation, with geometry optimization took approximately 5 hours on the same computer. The numerical cost of our eDFT/MDFT procedure might appear prohibitive but one should keep in mind that it is "naively" implemented and not optimized at this stage. Moreover, MDFT provides a full molecular description of the solvent structure which is not accessible with continuum methods. In fact, the description of the solvent is theoretically equivalent to simulation based methods such as Molecular Dynamics, which are even more computationally demanding.

### C. Calculations of the $2e^- / 2H^+$ Reduction Potentials

We now turn to the computation of the redox properties of a series of quinone derivatives. More specifically, we focus on two electrons reductions of benzoquinone derivatives (QR) into their hydroquinone ($H_2QR$) form:

$$QR + 2e^- + 2H^+ \rightarrow H_2QR. \quad (55)$$

It is not possible to compute directly the reaction free energy, or equivalently the redox potential, associated with a redox half reaction such as Eq. (55) because it involves the chemical potentials of both the electron and proton, which are ill-defined in solution. Instead, we introduce a reference redox couple—the unfunctionalized

benzoquinone/hydroquinone pair $Q/H_2Q$ and define the redox process through the following isodesmic reaction:

$$QR + H_2Q \rightarrow H_2QR + Q \quad (56)$$

This reaction avoids the explicit appearance of electrons and protons, and only involves neutral species, allowing us to compute the reaction free energy based solely on the free energies of the quinones involved.

It is worth noticing that this approach assumes that the two sub-reactions ($QR/H_2QR$ and $Q/H_2Q$) are chemically equivalent in terms of the redox process. In particular, it relies on the assumption that the electrons and protons involved in Eq. (55) are effectively identical and can be canceled out between the two half-reactions. This assumption is reasonable when comparing structurally similar molecules, such as benzoquinone derivatives, where redox mechanisms are expected to be analogous. However, if the electronic or protonation environments differ significantly between the systems, this simplification can introduce non-negligible errors in the computed redox properties.

The redox potential for the half-cell reaction appearing in Eq. (55), $E_{QR/H_2QR}^{\circ}$, is related to the reaction free energy, $\Delta_r G_{QR/H_2QR}$, through

$$\Delta_r G_{QR/H_2QR} = -2FE_{QR/H_2QR}^{\circ} \quad (57)$$

where $F = 96,485.3321$ C.mol$^{-1}$ is the Faraday constant. The overall reaction free energy of Eq. (56) is the difference of the half-cell reaction free energies of the QR/H$_2$QR and of the reference couple, $Q/H_2Q$. By using Eq. (57), this allows for expressing the redox potential of the considered couple as

$$E_{QR/H_2QR}^{\circ} = -\frac{\Delta_r G^{\circ}}{2F} + E_{Q/H_2Q}^{\circ} \quad (58)$$

The reduction potential of $Q/H_2Q$ is set to 0 V, and the redox potentials of functionalized quinones can then be computed from the free energy difference between the reduced, $H_2QR$, and oxydized form, $Q$.

## IV. APPLICATIONS

### A. Redox Potential calculations

We computed the redox potential for a set of 23 quinone/hydroquinone couples displayed in Fig. 1 using Eq. (58). The required reactions free energies were computed using the geometries optimized with the QM/MDFT scheme, including ZPE and thermochemical corrections. The predicted redox potentials are reported in Fig. 3 and compared to the experimental and theoretical results of Huynh et al. who employed electronic DFT (eDFT) combined with the conductor-like polarizable continuum model (C-PCM)¹⁶. As can be seen from Fig. 3, the redox potential predicted using
<figure><img src="image_2.png" /></figure>

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 3 while quinones from 12 to 23, circled in purple are only compared to DFT/C-PCM predictions in figure 4, due to the lack of experimental data.

the QM/MDFT are consistent with the eDFT/C-PCM predictions. With respect to experimental results, both methods underestimate the redox potential by a few tenth of a volt on average. To further test our approach, we extended the study to 12 additional quinone derivatives for which no experimental redox potentials were reported in Huynh's work. The redox potentials for these 12 new compounds, together with the 11 original ones, are presented in Fig. 4 where they are compared to Huynh's predictions using eDFT/C-PCM. Again, the agreement is satisfactory, with an average deviation of only 0.007 V between our results and theirs.

Overall, the prediction of the redox potential with the QM/MDFT method agrees quite well both with the experimental values and the eDFT/C-PCM method. This is rather satisfactory since we employed the simplest mean-field coupling between eDFT and MDFT. Moreover, the dispersion-repulsion interaction has been modeled with a simple Lennard-Jones interaction whose parameters were not optimized. Several potential improvements, such as including explicit polarization effects or refining the QM/MM coupling, remain to be explored to further enhance the accuracy.

## B. Electronic and Solvent Charge Densities

One clear advantage of the QM/MDFT approach over implicit solvation models is its ability to capture the detailed solvent structure around the solute. Specifically, the spatially and orientationally dependent solvent density $n(r, \Omega)$ is obtained directly from the self-consistent functional minimization. A more intuitive and visually accessible quantity is the inhomogeneous solvent charge

<figure><img src="image_3.png" /></figure>

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_{ext}^{MM}$, of Eq. (52). This potential is read by the MDFT program to optimize the MM density. Then, the electrostatic potential generated by the MM region, $v_{ext}^{QM}$, of Eq. (51) 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><img src="image_4.png" /></figure>

Figure 3. Comparison of $2e^-/2H^+$ reduction potentials (vs.Q,2H$^+/H_2Q$) 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$^{16}$.

density $n_c(r)$ (see Eq. (47)), which is derived from the molecular solvent distribution and allows further analysis of the solvation environment. The charge density directly impacts the electronic structure calculation through the electrostatic external potential $v_{ext}^{QM}$ it creates (see Eq. (51)).

The solvent charges and electronic densities in the plane of the aromatic ring of the Q/H$_2$Q couple are
<figure><img src="image_5.png" /><figcaption>Figure 4. Comparison of 2e⁻/2H⁺ reduction potentials (vs.Q,2H⁺/H₂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 black¹⁶.</figcaption></figure>

reported in Fig. 5 for the optimized solute geometries. Starting with the benzoquinone, the simplest molecule of our benchmark, the electronic density exhibits the expected pattern: it is largely concentrated around the atomic nuclei.

In contrast, the solvent charge density shows a more complex structure. Although not clearly visible in Fig. 5, there is a vanishing region in the vicinity of the solute, due to the exclusion of solvent molecules. There are two banana shape lobes of positive charge close to the oxygen atoms, with two marked maxima. These features reflect the formation of hydrogen bonds between the hydrogen atoms of nearby water molecules and the lone pairs of the carbonyl oxygens. Each of these positively charged lobes is followed by a negatively charged region, indicating that the water molecules are oriented such that the hydrogen atoms point toward the carbonyl groups while the oxygen atoms pointing outwards.

In the vicinity of the benzoquinone's hydrogen atoms, three negatively charged regions are observed: two aligned with the C-H bonds and a third, more intense one located on the perpendicular bisector of the H-H segment. This pattern indicates preferential orientation of water molecules such that their oxygen atoms face the hydrogen atoms of Q.

A diffuse region of low positive charge is also visible closer to the solute. This can be an artefact due to the usage of the SPC/E model, in which only the oxygen atom carries a Lennard-Jones site. The excluded volume, i.e. the region where the water density vanishes, is essentially governed by the Lennard-Jones interaction between the oxygen of the solvent water molecules and the solute. However, due to molecular rotations, hydrogen atoms of water molecules may enter this excluded zone. This is confirmed by examining the map of the oxygen density,

<figure><img src="image_6.png" /><figcaption>Figure 5. Electronic densities (a.u.)⁻³ computed with eDFT/MDFT around Q (top) and H₂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.</figcaption></figure>

where the frontier between the first solvation shell and the excluded volume essentially correspond to the beginning of negative charge surfaces.

The second solvation shell is also well captured, as evidenced by the oscillatory pattern of the charge distribution. Finally, as moving away from the solute, the charge density tends toward zero indicating that all orientations of the solvent become equivalent.

The solvation structure around H₂Q differs notably from that around Q as evidenced by the solvent charge density displayed in the bottom of Fig. 5. The most notable change is the diminution of the positively charged region near the oxygen atom, due to the presence of OH bond. It also modifies the negatively charged region, with the apparition of a marked maximum aligned with the OH bond, and of two lateral "bumps". The second solvation shell is also slightly modified.

The functionalization of quinone molecules naturally impacts their solvation structure. As an illustration, the equilibrium solvent charge density and electronic ground-state density around the tetrafluorinated quinone and hydroquinone are shown in Fig. 6. As expected, the electronic density exhibits an increase of the density in the
<figure><img src="image_7.png" /></figure>

<figure><img src="image_8.png" /><figcaption>Figure 6. Same as figure 5 for the QF₄/H₂QF₄.</figcaption></figure>

<figure><img src="image_9.png" /></figure>

<figure><img src="image_10.png" /><figcaption>Figure 7. Difference of solute electronic density (a.u.) in blue and red between eDFT/MDFT simulation and in vacuum. The solutes, represented in green, are Q (top) and Q4F (bottom)</figcaption></figure>

region where fluorine atoms replaced the hydrogen atoms. Regarding the solvent charge, the most notable change is that the QF₄ molecule is now surrounded by a marked region of positive charge, followed by a region of negative charge. In particular, there is a change of sign in the region in the vicinity of fluorine atoms compared to the benzoquinone, indicating a flip in the orientation of water molecules in these regions. It is noteworthy that the solvent charge is more polarized around the O atoms than around the F atoms.

The modification of the charge density in the reduced form are similar to what is observed for the unsubstituted hydroquinone *i.e.* the positive charge is reduced close to the O-H bond while negative charge increases.

Overall, MDFT effectively captures the non-trivial three-dimensional solvent structure around the quantum solute, as evidenced in Fig. 5 and Fig. 6. In particular, the tetrahedral shape of the charge distribution, resulting from hydrogen bonding, is well rendered. This is a clear advantage of MDFT compared to simpler continuum solvation models, which ignore the molecular nature of the solvent.

To assess the impact of the solvent on the electronic structure of the solute, we also computed the difference between the electronic densities obtained in the presence of solvent and in vacuum. These differences are displayed

in Fig. 7 for the Q and QF₄ couple. Blue regions correspond to an increase in electronic density compared to vacuum, while red regions indicate a depletion. Overall, there is a clear connection between the maps of solvent charge density displayed in Fig. 5 and Fig. 6 and the maps of electronic density difference.

In both Q and QF₄, electronic density increases primarily around the oxygen atoms, indicating strong polarization effects induced by the positively charged solvent environment in those regions. This coincides with a decrease in electronic density around the aromatic ring, particularly on the carbon atoms bonded to oxygen.

However, while the entire cycle undergoes a decrease in electronic density in QF₄, the picture is more subtle in Q where some density depletions are observed along the C-H bonds and along the 4 C-C bonds connected to the carbonyl. In Q, the H atoms are surrounded by negative solvent charges which tend to repel the electrons towards the center of the molecule while the F atoms in QF₄ are surrounded by a positive solvent charge density which tends to increase the polarization of the C-F bond with respect to vacuum. The stronger solvent-induced polarization of the oxygen atoms relative to the fluorine
<figure><img src="image_11.png" /><figcaption>Figure 8. Difference of solute electronic density (a.u.) in blue and red between eDFT/C-PCM simulation and in vacuum. The color code is the same as in Fig. 7</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>Figure 9. Difference of solute electronic density (a.u.) between eDFT/MDFT and eDFT/C-PCM simulation. The color code is the same as in Fig. 7</figcaption></figure>

atoms in QF₄ is consistent with the known differences in their electric dipole polarizabilities: the polarizability of oxygen is approximately 30-45% higher than that of fluorine⁵⁹⁻⁶¹.

In Fig. 8, we display the difference between the electronic densities computed using a C-PCM model and in vacuum for the Q/HQ couple. The maps are almost identical to the one displayed in Fig. 7, indicating that modeling the solvent with MDFT yields an electronic density comparable to that obtained with using a C-PCM model. A closer look reveals that the charges are slightly more intense when the C-PCM model is used.

To further examine this behaviour we report in Fig. 9 the difference between the electronic densities obtained with MDFT and C-PCM. Here, the blue regions indicate an excess of electronic density predicted by the eDFT/MDFT calculation with respect to eDFT/C-PCM one. For both molecules, the electronic density predicted around the O atoms is greater with the C-PCM model than with eDFT/MDFT. In benzoquinone, the electronic density surrounding the hydrogen and carbon atoms is lower with C-PCM than with MDFT.

In the case of QF₄, electrons are more localized around O and F atoms when the C-PCM is used than when

MDFT is. The trend is opposite around the aromatic ring, where MDFT predicts an excess of electrons with respect to C-PCM. The prediction of the electronic structure using both solvent models are in good agreement overall.

In summary, modeling the solvent using MDFT calculations clearly impacts the electronic density of the solute, which differs markedly with respect to in-vacuum predictions. These variations are well rationalized by examining the solvent charge computed by MDFT. The comparison with eDFT/C-PCM shows an overall good agreement even if some differences have been highlighted.

Importantly, the QM/MDFT calculations rely on only one set of *ad hoc* parameters—those from the Lennard-Jones force field used to model dispersion-repulsion interactions between solute and solvent. These parameters were taken directly from the generic GAFF2 force field and were not fine-tuned. The fact that MDFT-based predictions closely match those of the well-established and heavily parametrized C-PCM method is therefore quite encouraging.

In addition, MDFT allows for the prediction of the detailed structure of the solvent around the solute which is not possible with C-PCM but requires the use of more
expansive techniques, such as Molecular Dynamics.

V. CONCLUSIONS

In this paper, we propose a variational principle of the grand potential of a QM/MM system as a functional of the one body electronic and classical densities. This framework is constructed through a series of well-defined approximations, eventually leading to a practical mean-field computational scheme. We demonstrate the applicability of this approach by computing the redox properties of a series of quinones.

This work begins by describing in Sec. II A the generic QM/MM system treated here, which consists of a QM solute made of electrons and nuclei, both interacting with classical rigid solvent molecules, whose number is allowed to vary. This QM/MM system is therefore described in a semi-grand canonical statistical ensemble. We introduce the corresponding equilibrium density matrix, partition function and grand potential which are the three key objects characterizing the system. Under the Born-Oppenheimer and zero electronic entropy approximations, we show in Sec. IIB that the equilibrium density matrix can be written as the product of a nuclear and an electronic density matrices. Nevertheless, these QM density matrices are parametrized by the classical solvent coordinates, thereby introducing QM/MM correlation terms.

Based on this findings, we propose in Sec. IIC a variational formulation of the grand potential in terms of both QM density matrices and classical probability distributions. This expression is then simplified using a mean-field approximation which neglects QM/MM correlation effects, allowing us to write an upper bound to the exact grand potential.

By properly integrating over quantum electrons and classical degrees of freedom, we reformulate in Sec. IID the mean-field grand potential only in terms of nuclear density matrices. This naturally defines an effective potential experienced by the nuclei due to the presence of the electrons and classical molecules. This effective potential is then identified as a special case of QM/MM object established by some of us in a previous work³¹. This allows us to reformulate the effective potential as a coupled minimization problem over one-body QM and MM densities.

We then introduce the standard Harmonic approximation in Sec. IIE to obtain a tractable expression of the grand potential. Eventually, computing the grand potential requires the knowledge of the equilibrium geometry of the nuclei and of the Hessian matrix of the total effective potential experienced by the nuclei at this geometry. Details about the computation of nuclear gradients and Hessian matrix are presented in the Appendix 3.

Therefore, accessing the grand potential reduces to a sequence of density optimizations. First, the solvent density and the solute electronic density are optimized

for a fixed solute nuclear geometry. The self-consistent QM/MM procedure allowing the relaxation of the solvent density around the solute was presented in a previous work³³. Subsequently, the solute geometry is optimized for a fixed solvent configuration. This process is iterated until the overall convergence is reached. A schematic description of this QM/MM scheme is presented in Fig. 2.

To validate this methodology, we applied the proposed QM/MM scheme to the calculation of the redox properties of a series of functionalized quinones in aqueous solution. A selection of 23 quinone couples (see Fig. 1) was studied and predictions were compared to experimental data and to computations using electronic DFT in a C-PCM model extracted from Huynh et al.¹⁶. The quantitative results are presented in Sec. IV A.

The redox potentials predicted using QM/MDFT are in good agreement with the eDFT/C-PCM predictions, but both methods underestimate the experimental values. It is worth emphasizing that the MDFT framework offers a more detailed representation of the solvent structure surrounding the solute than conventional PCM or other continuum solvation models. This level of detail allowed us to thoroughly discuss the interplay between the electronic density and the surrounding solvent, as presented in Sec. IVB, where we plotted the solvent charge density around the solute for two quinone couples. Overall, we found that MDFT accurately captures the complex, three-dimensional solvent structure around the quantum solute.

To further assess the interplay between the solvent environment and the electronic density, we computed the difference between the electronic densities computed in the presence of the solvent moldeled by MDFT and in vacuum. As expected, the electronic density is impacted by the presence of a non trivial solvent charge density. These electronic density variation can be understood by simple electrostatic consideration. We also compared the electronic density variations from eDFT/C-PCM calculations with respect to the vacuum.

Eventually, the proposed eDFT/MDFT scheme is able to obtain comparable results to the well-known C-PCM, whose parameters have been subtly optimized over the years⁶²,⁶³, even if the only ad hoc parameters introduced in our model, namely the Lennard-Jones parameters, were not tuned at all. In addition, MDFT provides a detailed picture of the solvent structure around the solute, a feature that is out of reach for continuum models.

Overall, these results are rather encouraging, particularly considering there are still many routes to explore which might improve the proposed QM/MDFT scheme. Among those, the QM/MM coupling scheme can certainly be improved, using for instance some overlap models for the Pauli repulsion instead of LJ potentials. Another obvious limitation is the lack of polarizability of the classical solvent molecules in the current MDFT functional. We believe these are the two main sources of errors, which should be addressed in priority to improve the precision of the QM/MDFT approach.
## ACKNOWLEDGEMENTS

The authors want to thank Luc Belloni for providing the SPC/E direct correlation functions. M.L. acknowledges the Institute of Materials Science (iMAT) of the Alliance Sorbonne Université for financial support. This work was supported by the French National Research Agency under the France 2030 program (Grant ANR-22-PEBA-0002).

## AUTHOR DECLARATIONS

### Conflict of Interest

The authors have no conflicts to disclose.

## AUTHOR CONTRIBUTION

**Maxime Labat:** Conceptualization (equal); Data curation (equal); Investigation (equal); Methodology (equal); Project administration (equal); Software (equal); Supervision (equal); Validation (equal); Visualization (equal); Writing – original draft (equal); Writing – review & editing (equal). **Emmanuel Giner:** Conceptualization (equal); Investigation (equal); Methodology (equal); Project administration (equal); Software (equal); Supervision (equal); Validation (equal); Writing – original draft (equal); Writing – review & editing (equal). **Guillaume Jeanmairet:** Conceptualization (equal); Data curation (equal); Investigation (equal); Methodology (equal); Project administration (equal); Software (equal); Supervision (equal); Validation (equal); Visualization (equal); Writing – original draft (equal); Writing – review & editing (equal).

## DATA AVAILABILITY

The data that support the findings of this study are available from the corresponding author upon reasonable request.

## 1. Appendix: Derivation of $\Omega_{MF}$

We aim at deriving the variational principle for the grand potential involving only the nuclear density matrix. We start by inserting the mean-field approximation of the electron-nuclei-solvent density matrix of Eq. (27) in the grand potential functional for the QM/MM system

of Eq. (23)

$$ 
\begin{aligned}
\Omega_{MF} &= \min_{\hat{\rho}_{\text{MF}}^{\text{GC}}} (\mathcal{U}_{\text{int}}[\hat{\mathcal{P}}_{\text{MF}}^{\text{GC}}] - \frac{1}{\beta}\mathcal{S}_{\text{int}}[\hat{\mathcal{P}}_{\text{MF}}^{\text{GC}}]) \\
&= \min_{p^{\text{GC}}\hat{\rho}_{\text{N}}\hat{\rho}_{\text{el}}} (\mathcal{U}_{\text{int}}[p^{\text{GC}}\hat{\rho}_{\text{N}}\hat{\rho}_{\text{el}}] - \frac{1}{\beta}\mathcal{S}_{\text{int}}[p^{\text{GC}}\hat{\rho}_{\text{N}}]).
\end{aligned}
\quad (59) $$

As a reminder, $\hat{\rho}_{el}$ is not included in $\mathcal{S}_{int}$ due to the ZEE approximation.

Let us start with the entropic term $\mathcal{S}_{int}$ of Eq. (25). Thanks to the mean field approximation, it can be split into a purely nuclear quantum entropy, $\mathcal{S}_N$, of Eq. (30) and a purely classical term, $\mathcal{S}_{MM}$, as follow

$$ 
\begin{aligned}
\mathcal{S}_{\text{int}}[\hat{\mathcal{P}}_{\text{MF}}^{\text{GC}}] &= -\text{Tr}_{\text{GC}}\{p^{\text{GC}}\hat{\rho}_{\text{N}} \log(p^{\text{GC}}\hat{\rho}_{\text{N}})\} \\
&= -\text{Tr}_{\text{GC}}\{p^{\text{GC}}\hat{\rho}_{\text{N}}[\log(p^{\text{GC}}) + \log(\hat{\rho}_{\text{N}})]\} \\
&= -\text{Tr}\{\hat{\rho}_{\text{N}} \log(\hat{\rho}_{\text{N}})\} \\
&\quad - \sum_{N_c=0}^{\infty} \left[ \iint d\tilde{\mathbf{Q}}^{N_c} d\mathbf{P}^{N_c} p^{\text{GC}} \log(p^{\text{GC}}) \right] \\
&= \mathcal{S}_{\text{MM}}[p^{\text{GC}}] + \mathcal{S}_{\text{N}}[\hat{\rho}_{\text{N}}].
\end{aligned}
\quad (60) $$

We now move to the energetic term, $\mathcal{U}_{int}$, i.e.

$$ 
\begin{aligned}
\mathcal{U}_{\text{int}}[\hat{\mathcal{P}}_{\text{MF}}^{\text{GC}}] &= \text{Tr}_{\text{GC}}\{p^{\text{GC}}\hat{\rho}_{\text{N}}\hat{\rho}_{\text{el}}[\hat{T}_{\text{N}} + \mathcal{W}_{\text{N}}^{\text{MM}} \\
&\quad + \hat{T}_{\text{e}} + \hat{W}_{\text{ee}} + \hat{W}_{\text{Ne}}(\mathbf{R}^{\text{NN}}) + \hat{W}_{\text{NN}}(\mathbf{R}^{\text{NN}}) \\
&\quad + \mathcal{W}_{\text{e}}^{\text{MM}} + W + T - \mu N]\},
\end{aligned}
$$

where we explicit the whole Hamiltonian. Isolating the terms which depends only on $\mathbf{R}^{NN}$ leads to

$$ 
\begin{aligned}
\mathcal{U}_{\text{int}}[\hat{\mathcal{P}}_{\text{MF}}^{\text{GC}}] &= \text{Tr}\{\hat{\rho}_{\text{N}}[\hat{T}_{\text{N}} + \hat{W}_{\text{NN}}(\mathbf{R}^{\text{NN}})]\} + \text{Tr}_{\text{GC}}\{p^{\text{GC}}\hat{\rho}_{\text{N}}\hat{\rho}_{\text{el}}[\mathcal{W}_{\text{N}}^{\text{MM}} + \hat{T}_{\text{e}} + \hat{W}_{\text{ee}} + \hat{W}_{\text{Ne}}(\mathbf{R}^{\text{NN}}) + \mathcal{W}_{\text{e}}^{\text{MM}} \\
&\quad + W + T - \mu N]\}.
\end{aligned}
$$

Inserting the electron and classical particles embedding potentials defined in Eqs. (37)-(38)

$$ 
\begin{aligned}
\mathcal{U}_{\text{int}}[\hat{\mathcal{P}}_{\text{MF}}^{\text{GC}}] &= \\
&\quad \text{Tr}\{\hat{\rho}_{\text{N}}[\hat{T}_{\text{N}} + \hat{W}_{\text{NN}}(\mathbf{R}^{\text{NN}}) + V_{\text{el}}^{\text{N}}[\hat{\rho}_{\text{el}}; \mathbf{R}^{\text{NN}}] + V_{\text{MM}}^{\text{N}}[p^{\text{GC}}; \mathbf{R}^{\text{NN}}]]\} + \\
&\quad \text{Tr}_{\text{GC}}\{p^{\text{GC}}\hat{\rho}_{\text{N}}\hat{\rho}_{\text{el}}[\hat{T}_{\text{e}} + \hat{W}_{\text{ee}} + \mathcal{W}_{\text{e}}^{\text{MM}} + W + T - \mu N]\}.
\end{aligned}
$$

Gathering the electronic terms, we recognize the expression of the intrinsic electronic functional, $f_{el}$, given in Eq. (34). If the entropic contribution due to the classical solvent, $\mathcal{S}_{MM}$, is added we can similarly identify the intrinsic electronic functional, $f_{MM}^{GC}$, in Eq. (35)

$$ 
\begin{aligned}
\mathcal{U}_{\text{int}}[\hat{\mathcal{P}}_{\text{MF}}^{\text{GC}}] &= \frac{1}{\beta}\mathcal{S}_{\text{MM}}[p^{\text{GC}}] = \\
&\quad \text{Tr}\{\hat{\rho}_{\text{N}}[\hat{T}_{\text{N}} + \hat{W}_{\text{NN}}(\mathbf{R}^{\text{NN}})]\} + \\
&\quad \text{Tr}\{\hat{\rho}_{\text{N}}[V_{\text{el}}^{\text{N}}[\hat{\rho}_{\text{el}}; \mathbf{R}^{\text{NN}}] + V_{\text{MM}}^{\text{N}}[p^{\text{GC}}; \mathbf{R}^{\text{NN}}] + f_{\text{el}}[\hat{\rho}_{\text{el}}] + f_{\text{MM}}^{\text{GC}}[p^{\text{GC}}] \\
&\quad + E_{\text{QM}}^{\text{MM}}[\hat{\rho}_{\text{el}}, p^{\text{GC}}]]\}.
\end{aligned}
$$
By adding $S_N[\hat{\rho}_N]$, we can now go back to expression of the grand potential, $\Omega_{MF}$, which requires minimizations over $p^{GC}\hat{\rho}_N \hat{\rho}_{el}$

$$ \Omega_{MF} = \min_{p^{GC} \hat{\rho}_N \hat{\rho}_{el}} \left( \mathrm{Tr}\{\hat{\rho}_N[\hat{T}_N + \hat{W}_{NN}(\mathbf{R}^{NN}) + \beta^{-1} \ln(\hat{\rho}_N)] \right. \\ \left. + \hat{\rho}_N[V_{el}^N[\hat{\rho}_{el}; \mathbf{R}^{NN}] + V_{MM}^N[p^{GC}; \mathbf{R}^{NN}] + f_{el}[\hat{\rho}_{el}] \right. \\ \left. + f_{MM}^{GC}[p^{GC}] + f_{QM}^{MM}[\hat{\rho}_{el}, p^{GC}]] \right). $$

Here, if one perform a first optimization over the electronic and classical density matrix, i.e.

$$ \min_{p^{GC} \hat{\rho}_N \hat{\rho}_{el}} = \min_{\hat{\rho}_N} \{ \min_{p^{GC} \hat{\rho}_{el}} \} \quad (61) $$

one recover the expression of the total effective potential felt by the nuclei, $W_N^{eff}(\mathbf{R}^{NN})$ and finally the expression of the grand potential, $\Omega_{MF}$ of Eq.(29).

## 2. Appendix: The Harmonic Approximation

This Appendix briefly summarizes the basic equations and concept of harmonic approximation allowing us to obtain Eq. (43) as an estimate of the Grand potential. We therefore recapitulate here the main steps proposed originally by McQuarrie$^{38}$, which were further explicit by Ochterski in the documentation of the thermodynamic corrections as implemented in the Gaussian software$^{43,64}$.

We consider a system of quantum nuclei evolving in a generic potential, which in our case is $\hat{W}_{NN}^{tot}(\mathbf{R}^{NN})$. The first step of the harmonic approximation is to approximate up to the second-order the potential $\hat{W}_{NN}^{tot}(\mathbf{R}^{NN})$ around the equilibrium geometry, i.e.

$$ \hat{W}_{NN}^{tot}(\mathbf{R}^{NN}) \approx \hat{W}_{NN}^{tot}(\mathbf{R}_{eq}) + \frac{1}{2} \delta_{\mathbf{R}^{NN}}^{\dagger} \mathbf{h} \delta_{\mathbf{R}^{NN}}, \quad (62) $$

where $\mathbf{R}_{eq}$ is the equilibrium geometry verifying

$$ \nabla_{\mathbf{R}^{NN}} \hat{W}_{NN}^{tot}(\mathbf{R}^{NN})|_{\mathbf{R}^{NN}=\mathbf{R}_{eq}} = 0, \quad (63) $$

$\delta_{\mathbf{R}^{NN}} = \mathbf{R}^{NN} - \mathbf{R}_{eq}$ is the displacement vector with respect to $\mathbf{R}_{eq}$, and $\mathbf{h}$ is the Hessian matrix defined as

$$ \mathbf{h}_{ij} = \frac{\partial^2 \hat{W}_{NN}^{tot}(\mathbf{R}^{NN})}{\partial \mathbf{R}_i \partial \mathbf{R}_j} \bigg|_{\mathbf{R}^{NN}=\mathbf{R}_{eq}}. \quad (64) $$

Written in the basis of the $N_ν$ normal modes (i.e. the eigenvectors of the Hessian matrix), the total nuclear Hamiltonian become a sum of uncoupled harmonic oscillators, such that $U_N$ in Eq. (31) can be written as

$$ \begin{aligned} U_N[\hat{\rho}_N] &= \hat{W}_{NN}^{tot}(\mathbf{R}_{eq}) + \sum_{k=1}^{N_ν} \mathrm{Tr}_{QM}\{\hat{\rho}_N \hat{h}_{\omega_k}\}, \\ \hat{h}_{\omega_k} &= -\frac{1}{2\mu_k} \frac{\partial}{\partial q_k^2} + \frac{1}{2}\mu_k (\omega_k)^2 q_k^2, \end{aligned} \quad (65) $$

where $q_k, \mu_k$ are the displacement and effective mass associated to the normal mode $k$, respectively, and $\omega_k$ is the associated frequency defined as

$$ \omega_k = \sqrt{\frac{\lambda_k}{\mu_k}}, \quad (66) $$

where $\lambda_k$ is the corresponding eigenvalue of the Hessian matrix. Provided that $\mathbf{R}_{eq}$ is not a saddle point, the eigenvalues are all non negative, but nevertheless there are small eigenvalues corresponding to rotational and translational degrees of freedom. As such motions cannot be described by the harmonic oscillators, one simply neglects the lowest eigenvalues to retain only vibrational degrees of freedom. Under these conditions, minimizing the functional as in Eq. (29) with the approximated potential of Eq. (65) is equivalent to find the partition function of a sum of uncoupled oscillators. As a consequence, the total partition function is simply the product of the individual partition functions for each vibration mode

$$ Z_{\hat{\mathbf{Q}}^{N_c}} = e^{-\beta \hat{W}_{NN}^{tot}(\mathbf{R}_{eq})} \prod_{k=1}^{N_ν} Z_{\mathrm{vib}}^{k}, \quad (67) $$

where $Z_{\mathrm{vib}}^k$ is the partition function of the k-th vibration mode,

$$ Z_{\mathrm{vib}}^{k} = \frac{e^{-\frac{\beta \omega_{k}}{2}}}{1 - e^{-\beta \omega_{k}}}. \quad (68) $$

When this partition function is inserted into the definition of the grand potential

$$ \Omega = -\frac{1}{\beta} \log \left( Z_{\hat{\mathbf{Q}}^{N_c}} \right), \quad (69) $$

it yields the approximate expression of the grand potential of Eq. (43).

## 3. Appendix: Evaluating the nuclear gradient and Hessian

As shown in Sec. 2, the computation of $\Omega$ requires the gradient and Hessian of the total effective potential $\hat{W}_{NN}^{tot}(\mathbf{R})^{NN}$ with respect to the solute nuclear coordinates, and the non trivial part is that coming from $W_N^{eff}(\mathbf{R}^{NN})$. We provide here a summary of the main equations involved to compute these quantities.

### a. Computation of the total nuclear gradient of $W_N^{eff}(\mathbf{R}^{NN})$

Let us first define the optimal classical and quantum one-body densities for a given nuclear geometry

$$ (\rho_{\mathrm{opt}}(\mathbf{R}^{\mathrm{NN}}), n_{\mathrm{opt}}(\mathbf{R}^{\mathrm{NN}})) = \underset{\rho, n}{\operatorname{argmin}} w_{\mathrm{N}}^{\mathrm{eff}}[\rho, n, \mathbf{R}^{\mathrm{NN}}]. \quad (70) $$
Then, the gradient of $W_N^{\text{eff}}(\mathbf{R}^{\text{NN}})$ with respect to the nuclear geometry is defined as

$$ 
\begin{aligned}
\nabla_{\mathbf{R}^{\text{NN}}} W_N^{\text{eff}}(\mathbf{R}^{\text{NN}}) &= \nabla_{\mathbf{R}^{\text{NN}}} \min_{\rho, n} w_N^{\text{eff}}[\rho, n, \mathbf{R}^{\text{NN}}] \\
&= \nabla_{\mathbf{R}^{\text{NN}}} w_N^{\text{eff}}[\rho_{\text{opt}}(\mathbf{R}^{\text{NN}}), n_{\text{opt}}(\mathbf{R}^{\text{NN}}), \mathbf{R}^{\text{NN}}].
\end{aligned}
\quad (71) $$

The gradient contains therefore two terms: one coming from the explicit dependency of the functional $w_N^{\text{eff}}$ with respect to the nuclear geometry, and another one coming from the implicit dependence of $(\rho_{\text{opt}}(\mathbf{R}^{\text{NN}}), n_{\text{opt}}(\mathbf{R}^{\text{NN}}))$ with respect to the nuclear geometry. Because $(\rho_{\text{opt}}(\mathbf{R}^{\text{NN}}), n_{\text{opt}}(\mathbf{R}^{\text{NN}}))$ minimizes the functional $w_N^{\text{eff}}$, the equivalent of the Hellman-Feynman theorem can be used such that only the contribution from coming from the explicit dependency of the functional is non zero, which yields to

$$ 
\begin{aligned}
\nabla_{\mathbf{R}^{\text{NN}}} W_N^{\text{eff}}(\mathbf{R}^{\text{NN}}) = & (\nabla_{\mathbf{R}^{\text{NN}}} v_{\text{Ne}}(\mathbf{R}^{\text{NN}}) |\rho_{\text{opt}}(\mathbf{R}^{\text{NN}})) \\
& + (\nabla_{\mathbf{R}^{\text{NN}}} v_N^{\text{MM}}(\mathbf{R}^{\text{NN}}) |n_{\text{opt}}(\mathbf{R}^{\text{NN}})).
\end{aligned}
\quad (72) $$

In Eq. (72), one recognizes the usual electronic gradient term $(\nabla_{\mathbf{R}^{\text{NN}}} v_{\text{Ne}}(\mathbf{R}^{\text{NN}})|\rho_{\text{opt}}(\mathbf{R}^{\text{NN}}))$, and the other one is the contribution to the classical density, which is simply computed as an integral

$$ (\nabla_{\mathbf{R}^{\text{NN}}} v_N^{\text{MM}}(\mathbf{R}^{\text{NN}}) |n_{\text{opt}}(\mathbf{R}^{\text{NN}})) = \int d\tilde{Q} n_{\text{opt}}(\tilde{Q}) \nabla_{\mathbf{R}^{\text{NN}}} v_N^{\text{MM}}(\tilde{Q}, \mathbf{R}^{\text{NN}}). \quad (73) $$

Therefore, for each component of the nuclear coordinates, one just needs to compute the integral of the nuclear gradient of the nuclear-solvent interaction potential (i.e. $\nabla_{\mathbf{R}^{\text{NN}}} v_N^{\text{MM}}(\tilde{Q}, \mathbf{R}^{\text{NN}})$ in Eq. (73)) multiplied by the classical solvent density (i.e. $n_{\text{opt}}(\tilde{Q})$). Each of these integrals is typically approximated by a numerical integral as usually done in MDFT, and can then be added to the usual electronic gradient to obtain the total nuclear gradient for the QM/MM system.

### b. Computation of the exact and approximated Hessian matrix of $W_N^{\text{eff}}(\mathbf{R}^{\text{NN}})$

Regarding now the second-order derivative of $W_N^{\text{eff}}(\mathbf{R}^{\text{NN}})$, applying the linear response formalism to the QM/MM functional of Eqs. (39) and (40) leads to the following expression

$$ \frac{\partial^2 W_N^{\text{eff}}(\mathbf{R}^{\text{NN}})}{\partial \mathbf{R}_i \partial \mathbf{R}_j} = \mathbf{h}_{\text{Ne}}(\mathbf{R}_i, \mathbf{R}_j) + \mathbf{h}_N^{\text{MM}}(\mathbf{R}_i, \mathbf{R}_j) + \mathbf{h}_e^{\text{MM}}(\mathbf{R}_i, \mathbf{R}_j), \quad (74) $$

with

$$ 
\begin{aligned}
\mathbf{h}_{\text{Ne}}(\mathbf{R}_i, \mathbf{R}_j) = & \left( \frac{\partial^2 v_{\text{Ne}}(\mathbf{R}^{\text{NN}})}{\partial \mathbf{R}_i \partial \mathbf{R}_j} | \rho_{\text{opt}} \right) + \left( \nabla_{\mathbf{R}_i} v_{\text{Ne}}(\mathbf{R}^{\text{NN}}) | \nabla_{\mathbf{R}_j} \rho_{\text{opt}} \right) \\
& \quad \text{Therefore, under this approximation, no additional computation is needed to obtain the Hessian matrix.}
\end{aligned}
\quad (75) $$

$$ \mathbf{h}_N^{\text{MM}}(\mathbf{R}_i, \mathbf{R}_j) = \left( \frac{\partial^2 v_N^{\text{MM}}(\mathbf{R}^{\text{NN}})}{\partial \mathbf{R}_i \partial \mathbf{R}_j} \right|_{n_{\text{opt}}} + \left( \nabla_{\mathbf{R}_i} v_N^{\text{MM}}(\mathbf{R}^{\text{NN}}) \right) |_{\nabla_{\mathbf{R}_j} n_{\text{opt}}}, \quad (76) $$

$$ \mathbf{h}_e^{\text{MM}}(\mathbf{R}_i, \mathbf{R}_j) = (\nabla_{\mathbf{R}_i} n_{\text{opt}} | v_e^{\text{MM}} | \nabla_{\mathbf{R}_j} \rho_{\text{opt}}). \quad (77) $$

Here, the response vectors $\nabla_{\mathbf{R}^{\text{NN}}} \rho_{\text{opt}}$ and $\nabla_{\mathbf{R}^{\text{NN}}} n_{\text{opt}}$ can be obtained by solving the coupled linear response equations which, after some derivations, yield

$$ 
\begin{aligned}
\nabla_{\mathbf{R}^{\text{NN}}} \frac{\delta}{\delta \rho(\mathbf{r})} w_N^{\text{eff}}[\rho_{\text{opt}}, n_{\text{opt}}, \mathbf{R}^{\text{NN}}] = & \langle \hat{\chi}^{\text{el}} \cdot \nabla_{\mathbf{R}^{\text{NN}}} \rho_{\text{opt}} \rangle_{\mathbf{r}} \\
& + \langle \hat{\chi}_{\text{el}}^{\text{MM}} \cdot \nabla_{\mathbf{R}^{\text{NN}}} n_{\text{opt}} \rangle_{\mathbf{r}},
\end{aligned}
\quad (78) $$

$$ 
\begin{aligned}
\nabla_{\mathbf{R}^{\text{NN}}} \frac{\delta}{\delta n(\tilde{Q})} w_N^{\text{eff}}[\rho_{\text{opt}}, n_{\text{opt}}, \mathbf{R}^{\text{NN}}] = & \langle \hat{\chi}^{\text{MM}} \cdot \nabla_{\mathbf{R}^{\text{NN}}} n_{\text{opt}} \rangle_{\mathbf{r}} \\
& + \langle \hat{\chi}_{\text{MM}}^{\text{el}} \cdot \nabla_{\mathbf{R}^{\text{NN}}} \rho_{\text{opt}} \rangle_{\mathbf{r}},
\end{aligned}
\quad (79) $$

where the notation $\langle \cdot \rangle_{\mathbf{r}}$ denotes the convolution product

$$ \langle \hat{f} \cdot g \rangle_{\mathbf{r}} = \int d\mathbf{r}' f(\mathbf{r}, \mathbf{r}') g(\mathbf{r}'), \quad (80) $$

and where the response functions are defined as

$$ \chi^{\text{el}}(\mathbf{r}, \mathbf{r}')^{-1} = - \frac{\delta^2}{\delta\rho(\mathbf{r})\delta\rho(\mathbf{r}')} w_N^{\text{eff}}[\rho, n, \mathbf{R}^{\text{NN}}], \quad (81) $$

$$ \chi^{\text{MM}}(\tilde{Q}, \tilde{Q}')^{-1} = - \frac{\delta^2}{\delta n(\tilde{Q}) \delta n(\tilde{Q}')'} w_N^{\text{eff}}[\rho, n, \mathbf{R}^{\text{NN}}], \quad (82) $$

$$ \chi_{\text{el}}^{\text{MM}}(\tilde{Q}, \mathbf{r})^{-1} = - \frac{\delta^2}{\delta n(\tilde{Q}) \delta \rho(\mathbf{r})} w_N^{\text{eff}}[\rho, n, \mathbf{R}^{\text{NN}}], \quad (83) $$

$$ \chi_{\text{MM}}^{\text{el}}(\mathbf{r}, \tilde{Q})^{-1} = - \frac{\delta^2}{\delta \rho(\mathbf{r}) \delta n(\tilde{Q})} w_N^{\text{eff}}[\rho, n, \mathbf{R}^{\text{NN}}]. \quad (84) $$

One sees therefore that the total Hessian matrix contains both electron-electron (i.e. $\chi^{\text{el}}^{\text{MM}}$), electron-solvent (i.e. $\chi_{\text{el}}^{\text{MM}}$ and $\chi_{\text{MM}}^{\text{el}}$) and solvent-solvent ($\chi^{\text{MM}}$) response functions.

In the computations presented here, we neglected the response functions $\chi_{\text{MM}}^{\text{el}}$, $\chi_{\text{el}}^{\text{MM}}$ and $\chi^{\text{MM}}$. This leads then to the standard electronic linear response equations, which are nevertheless evaluated at the optimal electronic density $\rho_{\text{opt}}$ which takes into account the influence of the solvent, i.e.
## SUPPORTING INFORMATION

The supporting information file contains i) the redox potentials computed within our eDFT/MDFT scheme together with the experimental and eDFT/C-PCM data of Ref. 16, ii) the Lennard-Jones parameters used in our simulations, and iii) the Cartesian coordinates obtained with our eDFT/MDFT geometry optimization scheme.

## REFERENCES

¹Y. Ding, C. Zhang, L. Zhang, Y. Zhou, and G. Yu, Chem. Soc. Rev. **47**, 69 (2018).

²V. Singh, S. Kim, J. Kang, and H. R. Byon, Chem. Soc. Rev. **12**, 1 (2019).

³D. G. Kwabi, Y. Ji, and M. J. Aziz, Chemical Reviews **120**, 6467 (2020).

⁴J. Noack, N. Roznyatovskaya, T. Herr, and P. Fischer, Angewandte Chemie International Edition **54**, 9776 (2015).

⁵G. Kear, A. A. Shah, and F. C. Walsh, International Journal of Energy Research **36**, 1105 (2012).

⁶J. Winsberg, T. Hagemann, T. Janoschka, M. D. Hager, and U. S. Schubert, Angewandte Chemie International Edition **56**, 686 (2017).

⁷B. Huskinson, M. P. Marshak, C. Suh, S. Er, M. R. Gerhardt, C. J. Galvin, X. Chen, A. Aspuru-Guzik, R. G. Gordon, and M. J. Aziz, Nature **505** (2014), 10.1063/5.0019056.

⁸Q. Zhao, W. Huang, Z. Luo, L. Liu, Y. Lu, Y. Li, L. Li, J. Hu, H. Ma, and J. Chen, Science Advances **4**, 1761 (2018).

⁹B. Lee, Y. Ko, G. Kwon, S. Lee, K. Ku, J. Kim, and K. Kang, Joule **2**, 61 (2018).

¹⁰M. R. Gerhardt, L. Tong, R. Gómez-Bombarelli, Q. Chen, M. P. Marshak, C. J. Galvin, A. Aspuru-Guzik, R. G. Gordon, and M. J. Aziz, Advanced Energy Materials **7**, 1601488 (2017).

¹¹K. Sasaki, T. Kashimura, M. Ohura, Y. Ohsaki, and N. Ohta, Journal of The Electrochemical Society **137**, 2437 (1990).

¹²M. Shamsipur, A. Salimi, S. M. Golabi, H. Sharghi, and M. F. Mousavi, Journal of Solid State Electrochemistry **5**, 68 (2001).

¹³W. Wang, W. Xu, L. Cosimbescu, D. Choi, L. Li, and Z. Yang, Chem. Commun. **48**, 6669 (2012).

¹⁴S. Roy and P. S. Guin, Journal of The Electrochemical Society **162**, H124 (2014).

¹⁵P. S. Guin, S. Das, and P. C. Mandal, International Journal of Electrochemistry **2011**, 816202 (2011).

¹⁶M. T. Huynh, C. W. Anson, A. C. Cavell, S. S. Stahl, and S. Hammes-Schiffer, J Am Chem Soc **138**, 15903 (2016).

¹⁷B. Yang, L. Hoober-Burkhardt, F. Wang, G. K. S. Prakash, and S. R. Narayanan, J Electrochem Soc **161**, A1371 (2014).

¹⁸K. Wedege, E. Dražević, D. Konya, and A. Bentien, Sci Rep **6**, 1 (2016).

¹⁹A. Khetan, Batteries **9**, 24 (2023), number: 1 Publisher: Multi-disciplinary Digital Publishing Institute.

²⁰S. Er, C. Suh, M. P. Marshaka, and A. Aspuru-Guzik, Chem. Sci. **6**, 885 (2015).

²¹H. Kim, T. I. Goodson, and P. M. Zimmerman, The Journal of Physical Chemistry C **120**, 22235 (2016).

²²Q. Zhang, A. Khetan, and S. Er, Scientific Reports **10**, 22149 (2020).

²³S. Elhajj and S. Gozem, Journal of Chemical Theory and Computation (2024), 10.1021/acs.jctc.4c00602.

²⁴N. D. Mermin, Phys Rev **137**, A1441 (1965).

²⁵R. Evans, Adv Phys (1979).

²⁶W. S. B. Dwandaru and M. Schmidt, Physical Review E **83**, 061133 (2011).

²⁷J. Lischner and T. A. Arias, Physical Review B **84**, 125143 (2011).

²⁸R. Sundararaman, K. Letchworth-Weaver, K. A. Schwarz, D. Gunceler, Y. Ozhabes, and T. Arias, SoftwareX **6**, 278 (2017).

²⁹W. Tang, C. Cai, S. Zhao, and H. Liu, The Journal of Physical Chemistry C **122**, 20745 (2018).

³⁰J.-M. Lude, P. Tremblin, I. Chataigner, D. Borgis, and R. Spezia, "In Press" (2025), doi:10.26434/chemrxiv-2025-rb48g.

³¹G. Jeanmairet, M. Labat, and E. Giner, arXiv (2024), 10.48550/arXiv.2411.11821, 2411.11821.

³²G. Jeanmairet, M. Levesque, and D. Borgis, J Chem Theory Comput **16**, 7123 (2020).

³³M. Labat, E. Giner, and G. Jeanmairet, The Journal of Chemical Physics **161**, 014113 (2024).

³⁴S. Nielsen, R. Kapral, and G. Ciccotti, J. Chem. Phys. **115**, 5805 (2001).

³⁵M. Levy, Proc Natl Acad Sci U.S.A **76**, 6062 (1979), 16592733.

³⁶M. Levy, Phys Rev A **26**, 1200 (1982).

³⁷E. H. Lieb, Int J Quantum Chem **24**, 243 (1983).

³⁸D. A. McQuarrie, *Quantum Chemistry*, second edition ed. (Viva Books, 2003).

³⁹R. F. Ribeiro, A. V. Marenich, C. J. Cramer, and D. G. Truhlar, J Phys Chem B **115**, 14556 (2011).

⁴⁰S. Grimme, Chem Eur J **18**, 9955 (2012).

⁴¹Y.-P. Li, J. Gomes, S. Mallikarjun Sharada, A. T. Bell, and M. Head-Gordon, J Phys Chem C **119**, 1840 (2015).

⁴²J. A. Velmiskina, V. I. Malyshev, I. S. Gerasimov, and M. G. Medvedev, J Chem Phys **162**, 124115 (2025).

⁴³J. W. Ochterski, "Thermochemistry in Gaussian | Gaussian.com," (2000).

⁴⁴M. J. Frisch, G. W. Trucks, H. B. Schlegel, G. E. Scuseria, M. A. Robb, J. R. Cheeseman, G. Scalmani, V. Barone, G. A. Petersson, H. Nakatsuji, X. Li, M. Caricato, A. V. Marenich, J. Bloino, B. G. Janesko, R. Gomperts, B. Mennucci, H. P. Hratchian, J. V. Ortiz, A. F. Izmaylov, J. L. Sonnenberg, D. Williams-Young, F. Ding, F. Lipparini, F. Egidì, J. Goings, B. Peng, A. Petrone, T. Henderson, D. Ranasinghe, V. G. Zakrzewski, J. Gao, N. Rega, G. Zheng, W. Liang, M. Hada, M. Ehara, K. Toyota, R. Fukuda, J. Hasegawa, M. Ishida, T. Nakajima, Y. Honda, O. Kitao, H. Nakai, T. Vreven, K. Throssell, J. A. Montgomery, Jr., J. E. Peralta, F. Ogliaro, M. J. Bearpark, J. J. Heyd, E. N. Brothers, K. N. Kudin, V. N. Staroverov, T. A. Keith, R. Kobayashi, J. Normand, K. Raghavachari, A. P. Rendell, J. C. Burant, S. S. Iyengar, J. Tomasi, M. Cossi, J. M. Millam, M. Klene, C. Adamo, R. Cammi, J. W. Ochterski, R. L. Martin, K. Morokuma, O. Farkas, J. B. Foresman, and D. J. Fox, "Gaussian 16 Revision C.01," (2016), gaussian Inc. Wallingford CT.

⁴⁵Q. Sun, X. Zhang, S. Banerjee, P. Bao, M. Barbry, N. S. Blunt, N. A. Bogdanov, G. H. Booth, J. Chen, Z.-H. Cui, J. J. Eriksen, Y. Gao, S. Guo, J. Hermann, M. R. Hermes, K. Koh, P. Koval, S. Lehtola, Z. Li, J. Liu, N. Mardirossian, J. D. McClain, M. Motta, B. Mussard, H. Q. Pham, A. Pulkin, W. Purwanto, P. J. Robinson, E. Ronca, E. R. Sayfutyarova, M. Scheurer, H. F. Schurkus, J. E. T. Smith, C. Sun, S.-N. Sun, S. Upadhyay, L. K. Wagner, X. Wang, A. White, J. D. Whitfield, M. J. Williamson, S. Wouters, J. Yang, J. M. Yu, T. Zhu, T. C. Berkelbach, S. Sharma, A. Y. Sokolov, and G. K.-L. Chan, The Journal of Chemical Physics **153**, 024109 (2020).

⁴⁶F. Neese, F. Wennmohs, U. Becker, and C. Riplinger, J Chem Phys **152** (2020), 10.1063/5.0004608.

⁴⁷L. Ding, M. Levesque, D. Borgis, and L. Belloni, The Journal of Chemical Physics **147**, 094107 (2017).

⁴⁸A. D. Becke, The Journal of Chemical Physics **98**, 5648 (1993).

⁴⁹C. Lee, W. Yang, and R. G. Parr, Phys. Rev. B **37**, 785 (1988).

⁵⁰W. J. Hehre, R. Ditchfield, and J. A. Pople, J. Chem. Phys. **56**, 2257 (1972).

⁵¹P. C. Hariharan and J. A. Pople, Theor. Chim. Acta **28**, 213 (1973).

⁵²M. M. Francl, W. J. Pietro, W. J. Hehre, J. S. Binkley, M. S. Gordon, D. J. DeFrees, and J. A. Pople, J. Chem. Phys. **77**, 3654 (1982).
<figure><img src="image_13.png" /><figcaption>Figure 10. For Table of Contents Only</figcaption></figure>

<footnote n="53">⁵³T. Clark, J. Chandrasekhar, G. W. Spitznagel, and P. V. R. Schleyer, *J. Comput. Chem.* **4**, 294 (1983).</footnote>

<footnote n="54">⁵⁴G. W. Spitznagel, T. Clark, P. v. R. Schleyer, and W. J. Hehre, *J. Comput. Chem.* **8**, 1109 (1987).</footnote>

<footnote n="55">⁵⁵H. J. C. Berendsen, J. R. Grigera, and T. P. Straatsma, *The Journal of Physical Chemistry* **91**, 6269 (1987).</footnote>

<footnote n="56">⁵⁶L. Belloni, *The Journal of Chemical Physics* **147**, 164121 (2017).</footnote>

<footnote n="57">⁵⁷X. He, V. H. Man, W. Yang, T.-S. Lee, and J. Wang, *The Journal of Chemical Physics* **153**, 114502 (2020).</footnote>

<footnote n="58">⁵⁸L.-P. Wang and C. Song, *The Journal of Chemical Physics* **144**, 214108 (2016).</footnote>

<footnote n="59">⁵⁹R. R. Teachout and R. T. Pack, *Atomic Data and Nuclear Data Tables* **3**, 195 (1971).</footnote>

<footnote n="60">⁶⁰K. Andersson and A. J. Sadlej, *Phys. Rev. A* **46**, 2356 (1992).</footnote>

<footnote n="61">⁶¹P. Schwerdtfeger and J. K. Nagle, *Molecular Physics* **117**, 1200 (2019).</footnote>

<footnote n="62">⁶²V. Barone and M. Cossi, *The J. Phys. Chem. A* **102**, 1995 (1998).</footnote>

<footnote n="63">⁶³M. Cossi, N. Rega, G. Scalmani, and V. Barone, *Journal of Computational Chemistry* **24**, 669 (2003).</footnote>

<footnote n="64">⁶⁴"Vibrational Analysis in Gaussian | Gaussian.com," (2025), [Online; accessed 10. Jun. 2025].</footnote>