# Analytic Gradients and Periodic Boundary Conditions for Direct Reaction Field Polarizable QM/MM with Electrostatic Potential Fitting

Thomas Fay, Miquel Huix-Rotllant, Nicolas Ferré

► To cite this version:

Thomas Fay, Miquel Huix-Rotllant, Nicolas Ferré. Analytic Gradients and Periodic Boundary Conditions for Direct Reaction Field Polarizable QM/MM with Electrostatic Potential Fitting. Journal of Chemical Theory and Computation, 2025, 21 (18), pp.8773-8786. 10.1021/acs.jctc.5c00863. hal-05297211

HAL Id: hal-05297211

https://cnrs.hal.science/hal-05297211v1

Submitted on 4 Oct 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>

Distributed under a Creative Commons Attribution 4.0 International License
# Analytic gradients and periodic boundary conditions for direct reaction field polarizable QM/MM with electrostatic potential fitting

Thomas P. Fay,¹,², a) Miquel Huix-Rotllant,¹ and Nicolas Ferré¹

¹) Aix Marseille Univ, CNRS, ICR, 13397 Marseille, France

²) Department of Chemistry and Biochemistry, University of California, Los Angeles, California 90095, United States

Our recently developed Direct Reaction field with ESPF Embedding Model (DREEM) method offers an efficient and physically rigorous framework for incorporating polarizable molecular mechanics (MM) environments into quantum mechanics/molecular mechanics (QM/MM) simulations. By coupling the QM and MM regions through the instantaneous MM electrostatic polarization response to QM charge density fluctuations, DREEM enables consistent treatment of ground and excited electronic states, capturing electronic state-specific polarization and dispersion effects absent in conventional mean-field or linear response approaches. The use of electrostatic potential fitting (ESPF) approximation method to describe charge density fluctuations significantly improves the computational efficiency compared the integral-exact direct reaction field. In this work, we present two methodological advancements to extend the applicability of DREEM to realistic condensed-phase simulations. Firstly the development of efficient analytic energy gradients, enabling geometry optimization, transition state searches, and molecular dynamics, and secondly a formulation of periodic boundary conditions (PBC) compatible with the DREEM framework. These capabilities are implemented in the open-source OpenESPF code, interfacing PySCF and OpenMM for high-performance QM and MM calculations. We demonstrate that the resulting implementation enables practical simulations of excited-state optical properties in periodic polarizable environments, where we calculate the fluorescence spectrum of acetone in water including quantum vibronic and non-Condon effects. This paves the way for predictive modeling of photochemical reactivity and spectroscopy in complex systems where environment polarization is important.

## I. INTRODUCTION

Accurate computational models for chemical reactivity, optical properties, and photochemical processes in complex condensed phase are important in a understanding a huge range of chemical, biological and materials systems. Quantum Mechanics/Molecular Mechanics (QM/MM) models enable one to describe the molecular fragments of interest using accurate QM methods, where electrons and nuclei are both treated explicitly, whilst describing the condensed phase environment with in atomic detail but with a far simpler and more computationally efficient MM force field.¹,² These methods circumvent the large computational cost that would come with treating electrons in the whole system explicitly, and therefore molecular dynamics and accurate thermal configuration sampling is routinely possible with QM/MM methods.²⁻⁵

Most MM force fields use a fixed charge description of electrostatic interaction between atoms, which cannot account for fluctuating electronic polarisation. In order to account for this, polarizable force fields have been developed, such as the AMOEBA force field,⁶,⁷ which use a coarse-grained description of electronic polarization through the use of induced dipole or charge equilibration models. Combining polarizable MM force fields with QM/MM therefore in principle accounts for the response of electronic polarization in condensed phase environments to changes in electron density in the QM region, which could be important for describing the interaction between the QM and MM regions.⁸,⁹ Several methods exist for

describing the embedding of the QM region in a polarizable MM environment.¹⁰⁻¹⁶ The most commonly used approach is to use a mean-field approach, which couples the QM and MM regions through the average charge density in the QM region $\langle\rho(r)\rangle$. The apparent simplicity of this scheme is appealing, however one runs into conceptual and computational difficulty when considering electronic excited states. If one only uses the ground state charge density to couple to the polarizable MM region, state-specific electronic polarization effects are completely missed, but coupling to the excited state charge density leads to broken orthogonality between ground and excited states, and also adds considerable computational overhead because the excited states and MM polarization need to be solved for self-consistently. Linear response schemes¹³,¹⁷ have been proposed but these cannot account for the full excited state specific polarization.¹⁷

In order to avoid some of the aforementioned issues, there has recently been a renewed interest in the alternative Direct Reaction Field (DRF) scheme.¹⁴,¹⁵,¹⁸⁻²¹ In this approach the MM polarization couples to the instantaneous QM region charge density, as encoded in the charge density operator $\hat{\rho}(r)$, and therefore the QM/MM embedding is treated entirely at the QM Hamiltonian operator level. This means ground and excited states are all treated on the same footing, and state specific polarization is automatically included, together with state-specific QM-MM dispersion interactions. The cost of the DRF method however is the need to calculate the MM polarization response to any possible charge density fluctuation in the QM region, rather than just the average. This introduces new screened two-electron interactions and a sizeable computational overhead.¹⁵ This means that for self-consistent field (SCF) electronic structure methods like Hartree-Fock (HF) or density functional theory (DFT), the computational cost of

<footnote>a) Corresponding author: thomaspfay@ucla.edu; Present address: Department of Chemistry and Biochemistry, University of California, Los Angeles, California 90095, United States</footnote>
the integral exact formulation of DRF (IEDRF)¹⁴,¹⁵ scales as $O(N_{AO}^3 N_{MM})$ when the size of atomic-orbital basis set ($N_{AO}$) is large, or at least $O(N_{AO}^2 N_{MM}^2)$ when solving the induction equations for the polarization response for the $N_{MM}$ MM atoms dominates the cost of the calculation.²¹ This motivated us to use electrostatic potential fitted (ESPF) multipole operators²² to approximate the DRF embedding hamiltonian, and in a previous paper we introduced the approximate Direct Reaction with the ESPF Embedding Model (DREEM) method (which we previously gave the much more cumbersome abbreviation ESPF-DRF).²¹ We describe the DREEM method in more detail below, but in short the QM charge density operator is approximated with an atom-centered multipole expansion $\hat{\rho}(\mathbf{r}) \approx \sum_a \hat{q}_a \rho_a(\mathbf{r})$, which yields a compact, low-rank approximation to the DRF hamiltonian. For SCF methods this reduces the computational scaling to $O(N_{AO}^3 N_{QM})$ when the orbital basis is large, or at least $O(N_{QM} N_{MM}^2)$ when solution of the MM induction equations is the computationally limiting step. For large MM systems and large basis sets this represents a huge potential computational saving.

Furthermore many QM/MM schemes simplify long-range electrostatics in condensed phase systems by setting interactions beyond a certain cut-off distance to zero. However electrostatic interactions decay very slowly with a $1/r$ distance dependence, so the validity of models with finite ranged electrostatics can be questionable,²³ particularly for polarizable environments.²⁴ As a result of this, finite sized models can suffer from spurious effects arising from the unphysically large surface area to volume ratio. To overcome this MM studies typically use periodic boundary conditions (PBC), where the system is replicated infinitely in all directions, to eliminate phase boundaries from the simulated system. The use of PBC can also speed-up calculations by reducing the size of the system needed to represent the condensed phase environment and enabling efficient treatment of electrostatics through Ewald based methods, such as the particle-mesh Ewald (PME) method,²⁵ which reduces the computational cost of evaluating electrostatic energies from $O(N^2)$ to $O(N \log N)$.²⁶ Whilst simulations with PBC have many advantages in condensed phase systems,³,²³ it can be complex to combine QM/MM with periodic boundary conditions. There have been several recently developed approaches to combine QM/MM using fixed charge force fields with PBC,³,²³,²⁴,²⁶⁻²⁸ but the IEDRF and DREEM methods have not yet been extended to periodic systems.

In order to extend the DREEM method, and enable more realistic and complex simulations to be performed, a PBC implementation is needed, as well as analytic gradients of the DREEM energy. This would make many chemically useful applications of the DREEM method possible like geometry optimizations, transition state searches, and condensed phase molecular dynamics. In this paper we derive both efficient analytical gradients and a formulation of periodic boundary conditions for the DREEM method. These developments are implemented in a new code, OpenESPF,²⁹ which uses PySCF and OpenMM as back-ends for efficiently performing QM and MM calculations respectively. The DREEM method as implemented in OpenESPF is efficient enough to perform excited state molecular dynamics simulations with PBC, as we shall

demonstrate below. We start with a brief description of the DRF and DREEM methods, and move on to describe how we combine the DREEM method with periodic boundary conditions. We then describe how to efficiently calculate analytic energy gradients and how these are implemented in the ESPF code. We then demonstrate out implementation through a set of example applications comparing PBC and droplet models, and demonstrating that stable NVE (constant energy) molecular dynamics can be performed with our method and implementation. In our final example we apply the DREEM method to run excited state molecular dynamics for acetone in water and from this calculate quantum vibronic, non-Condon and solvent effects on acetone's fluorescence spectrum.

## II. THEORY

### A. Direct Reaction Field QM/MM

In a QM/MM model the system is partitioned into a set of QM atoms, which are described using an explicit treatment of the electron density, and a set of MM atoms described with a far simpler molecular mechanics force-field which implicitly accounts for the electron density and energy for these atoms. The total energy can therefore be expressed as a sum of parts for the isolated QM and MM regions, and an interaction term

$$E = E_{QM} + E_{MM} + E_{int}. \quad (1)$$

When the MM system is described using a polarizable force field, there are several possible ways to describe the interaction between the QM and MM region, which vary in their complexity. Here we focus on the direct reaction field approach, detailed below.

In what follows we will describe the DRF method in the specific case of a simple induced dipole model for MM region polarizability, and leave a description of the technicalities that arise when using more complex MM force fields like the AMOEBA model⁶ until later. For a polarizable force field, where the electronic polarization is represented with a set of atom-centred induced dipoles $\mu_A$, the total electrostatic energy as a functional of the set of induced dipoles $\mu$ and the fixed charge density $\rho(\mathbf{r})$ is¹⁸,¹⁹

$$\mathcal{F}[\mu, \rho] = E_0[\rho] - \mathcal{E}[\rho] \cdot \mu + \frac{1}{2} \mu \cdot (\alpha^{-1} + \mathbf{T}) \cdot \mu \quad (2)$$

where $\alpha$ is a diagonal matrix of polarizabilities of each induced dipole, $\mathbf{T}$ is a tensor describing the interaction between dipoles. $E_0[\rho]$ describes the energy of the fixed charge density $\rho$ on its own, given by

$$E_0[\rho] = \frac{1}{2} \int d\mathbf{r} \int d\mathbf{r}' \rho(\mathbf{r}) \rho(\mathbf{r}') \phi(\mathbf{r} - \mathbf{r}'), \quad (3)$$

where $\phi(\mathbf{r}) = 1/|\mathbf{r}|$ is the Coulomb potential. $\mathcal{E}[\rho]$ is a vector of the electric fields $\mathcal{E}_A[\rho]$ at each induced dipole, $A$, generated by the charge density $\rho$. These electric fields are given by

$$\mathcal{E}_A[\rho] = - \int d\mathbf{r} \rho(\mathbf{r}) \nabla_{\mathbf{r}} \phi(\mathbf{R}_A^{\text{MM}} - \mathbf{r}). \quad (4)$$
It is important to note that the electric fields are a linear functional of $\rho(\mathbf{r})$. In order to obtain the polarization energy, it is assumed that the induced dipoles equilibrate instantaneously to the fixed charge density $\rho$, and therefore the total electrostatic energy, including polarization is

$$E_{\text{pol}}[\rho] = \min_{\mu} \mathcal{F}[\mu, \rho] = E_0[\rho] - \frac{1}{2} \mathcal{E}[\rho] \cdot \mu[\rho] \quad (5)$$

$$\mu[\rho] = (\alpha^{-1} + \mathbf{T})^{-1} \mathcal{E}[\rho]. \quad (6)$$

We again note here that the induced dipoles $\mu[\rho]$ are a linear functional of the fixed charge density, so overall the polarization energy is a quadratic functional of the fixed charge density.

Within the DRF framework, we assume that the MM dipoles respond instantaneously to the QM charge density, so the interaction Hamiltonian between QM and MM regions can be obtained by replacing the total charge density $\rho_{\text{tot}}(\mathbf{r})$ of the system with a sum of the fixed MM charge density $\rho_{\text{MM}}(\mathbf{r})$ and the QM charge density operator $\hat{\rho}(\mathbf{r})$

$$\hat{\rho}_{\text{tot}}(\mathbf{r}) = \rho_{\text{MM}}(\mathbf{r}) + \hat{\rho}(\mathbf{r}), \quad (7)$$

where the QM charge density operator $\hat{\rho}(\mathbf{r})$ (in atomic units) is

$$\hat{\rho}(\mathbf{r}) = \sum_{A=1}^{N_{\text{QM}}} Z_A \delta(\mathbf{r} - \mathbf{R}_A) - \sum_{i=1}^{N_c} \delta(\mathbf{r} - \hat{\mathbf{r}}_i), \quad (8)$$

in which $\mathbf{R}_A$ are positions of the QM nuclei and $\hat{\mathbf{r}}_i$ are QM electron position operators. After subtracting the electrostatic interactions involving only QM and only MM charge densities, since these are already accounted for in the $E_{\text{QM}}$ and $E_{\text{MM}}$ terms in the QM/MM energy, Eq. (1), the remaining interaction Hamiltonian is

$$\hat{V}_{\text{int}} = \int d\mathbf{r} \int d\mathbf{r}' \rho_{\text{MM}}(\mathbf{r}') \phi(\mathbf{r} - \mathbf{r}') \hat{\rho}(\mathbf{r}) \\ - \mu[\rho_{\text{MM}}] \cdot \mathcal{E}[\hat{\rho}] - \frac{1}{2} \mu[\hat{\rho}] \cdot \mathcal{E}[\hat{\rho}] \quad (9)$$

Overall this is a quadratic functional of the QM charge density, and therefore the QM/MM interaction Hamiltonian involves both one and two-electron terms, where the two electron terms describe the screening of the electron-electron repulsion by MM induced dipoles.¹⁴,¹⁵

In the integral exact formulation of the DRF method (IEDRF),¹⁴,¹⁵ zero-, one- and two-electron terms in the interaction Hamiltonian, $V^{(0)}$, $\hat{V}^{(1)}$ and $\hat{V}^{(2)}$, are evaluated exactly (after appropriate damping functions are introduced to account for the finite size of the MM charge densities and make integrals finite).¹⁴,¹⁵ For example, additional two-electron integrals in an atomic orbital $\chi_\mu$ basis have the form

$$(\mu\nu|\hat{V}^{(2)}|\mu'\nu') = -\mu[\rho_{\mu\nu}] \cdot \mathcal{E}[\rho_{\mu'\nu'}], \quad (10)$$

where $\rho_{\mu\nu}(\mathbf{r}) = -\chi_\mu^*(\mathbf{r})\chi_\nu(\mathbf{r})$. So for IEDRF evaluation of modified exchange matrices for self-consistent field electronic structure methods initially requires $O(N_{\text{AO}}^2)$ solutions of the MM induction equations, which typically scales at least as $O(N_{\text{MM}}^2)$, followed by $O(N_{\text{AO}}^3 N_{\text{MM}})$ operations to finally construct the exchange matrix.¹⁴,¹⁵

## B. Approximate DRF with DREEM

In order to obtain a more efficient but approximate direct-reaction field method we take a step back and first assume the charge density of the QM region (treated classically for now) can be approximated with an atom-centred multipole expansion of the form²²

$$\rho(\mathbf{r}) \approx \sum_{a=1}^{N_Q} q_a \rho_a(\mathbf{r}) \quad (11)$$

where $q_a$ is multipole operator centred at a particular point, and $\rho_a(\mathbf{r})$ is a point multipole charge density, for example for a point charge at $\mathbf{R}_A$, $\rho_a(\mathbf{r}) = \delta(\mathbf{r} - \mathbf{R}_A)$. We choose these multipoles to be centred on atoms in the QM region, but they could also be placed on additional virtual sites. If atom centred charges are used for every atom $N_Q = N_{\text{QM}}$ and if charges and dipoles are used then $N_Q = 4N_{\text{QM}}$. If we insert the multipole expansion for the QM charge density into the equation for the polarization energy, with the multipoles collected into a vector $\mathbf{q}$, we obtain the following expansion for the polarization energy (excluding direct interactions between the multipoles),

$$E_{\text{pol}}(\mathbf{q}) = u^{(\text{MM})} + \sum_{a} u_{a}^{(\text{QM/MM})} q_{a} + \frac{1}{2} \sum_{ab} u_{ab}^{(\text{QM})} q_{a} q_{b} \quad (12)$$

where the terms in the expansion are given by

$$\begin{align}
u^{(\text{MM})} &= E_{\text{pol}}[\rho_{\text{MM}}] \tag{13} \\
u_{a}^{(\text{QM/MM})} &= \int d\mathbf{r} \int d\mathbf{r}' \rho_{\text{MM}}(\mathbf{r}') \phi(\mathbf{r} - \mathbf{r}') \rho_{a}(\mathbf{r}) - \mu[\rho_{\text{MM}}] \cdot \mathcal{E}[\rho_{a}] \tag{14}
\end{align}$$

$$u_{ab}^{(\text{QM})} = -\mu[\rho_a] \cdot \mathcal{E}[\rho_b]. \quad (15)$$

So we see that in order to calculate the energy for any set of multipoles $\mathbf{q}$, we only need to solve the induction equations for $N_Q$ point multipole densities, in order to obtain the quadratic expansion above.

Returning now to the DRF method, we can apply the atom-centred multipole expansion to the charge density arising from every $\rho_{\mu\nu}$ density to construct the ESPF approximation to the electron charge density²¹

$$\rho_{\mu\nu}(\mathbf{r}) \approx \sum_{a=1}^{N_Q} Q_{a\mu\nu} \rho_a(\mathbf{r}). \quad (16)$$

The ESPF multipole operator matrix elements, $Q_{a\mu\nu}$, are fitted to reproduce the electrostatic potential generated by the $\rho_{\mu\nu}(\mathbf{r})$ charge density on a grid of points around the QM region, ensuring that the ESPF operators accurately approximate electrostatic interactions.²² The atom centred multipole operators are constrained to correctly reproduce the total charge for the $\rho_{\mu\nu}(\mathbf{r})$ charge distribution, and dipole moment if dipoles are included, which ensures their accuracy at long-range (with $O(R^{-2})$ error for charges, and $O(R^{-3})$ error when dipoles are additionally included).²¹ Details of the ESPF operators are
given in the supporting information. Combining these with the nuclear point charge densities, the QM charge density operator can be approximated as

$$ \hat{\rho}(\mathbf{r}) \approx \sum_a \left( Z_a + \sum_{i=1}^{N_c} \hat{Q}_{a,i} \right) \rho_a(\mathbf{r}) = \sum_a \hat{q}_a \rho_a(\mathbf{r}). \quad (17) $$

Here $Z_a$ is the nuclear charge contribution to point multipole $a$, which is $Z_A$ for a point charge at nucleus A and zero otherwise, and $\hat{Q}_{a,i}$ is the multipole operator $a$ acting on electron $i$ in the QM system. Using this to approximate the QM charge density, and combining this with the polarisation energy for an arbitrary set of multipoles, we obtain the DREEM interaction Hamiltonian,²¹

$$ \hat{V}_{\text{int}} = \sum_{a} u_{a}^{(\text{QM/MM})} \hat{q}_{a} + \frac{1}{2} \sum_{ab} u_{ab}^{(\text{QM})} \hat{q}_{a} \hat{q}_{b}. \quad (18) $$

This has zero, one and two-electron contributions, denoted $V^{(0)}$, $\hat{V}^{(1)}$, and $\hat{V}^{(2)}$ respectively. The zero-electron is given explicitly by

$$ V^{(0)} = \sum_a Z_a u_a^{(\text{QM/MM})} + \frac{1}{2} \sum_{ab} u_{ab}^{(\text{QM})} Z_a Z_b. \quad (19) $$

The one-electron term involves the interaction between electrons and polarization of MM induced dipoles in the environment by nuclei and the MM fixed charges, and a self-interaction term from direct interaction with the polarizable environment. The matrix elements of this one-electron Hamiltonian term in a given orbital basis are

$$ V_{\mu\nu}^{(1)} = \sum_{a} \left( u_{a}^{(\text{QM/MM})} + \sum_{b} u_{ab}^{(\text{QM})} Z_{b} \right) Q_{a\mu\nu} + \frac{1}{2} \sum_{ab} u_{ab}^{(\text{QM})} [Q_{a}Q_{b}]_{\mu\nu}, \quad (20) $$

where the matrix elements of $\hat{Q}_a \hat{Q}_b$ are

$$ [Q_a Q_b]_{\mu\nu} = \sum_{\mu'\nu'} Q_{a\mu\mu'} [S^{-1}]_{\mu'\nu'} Q_{b\nu'\nu}, \quad (21) $$

in which $[S]_{\mu\nu} = S_{\mu\nu} = \langle \chi_{\mu} | \chi_{\nu} \rangle$ is the overlap matrix for the orbital basis. Lastly, the two-electron term $\hat{V}^{(2)}$ gives rise to additional contributions to the two-electron integrals which are given by

$$ (\mu\nu|\hat{V}^{(2)}|\mu'\nu') = \sum_{ab} u_{ab}^{(\text{QM})} (\mu\nu|\hat{Q}_{a,1} \hat{Q}_{a,2}|\mu'\nu') \quad (22) $$

$$ = \sum_{ab} u_{ab}^{(\text{QM})} Q_{a\mu\nu} Q_{b\mu'\nu'}. \quad (23) $$

Eqs. (18)–(23) define the DREEM method. This is essentially a low-rank approximation to the full IEDRF two-electron integral. Exchange matrices involving this term can now be evaluated with a computational cost of $O(N_{AO}^3 N_Q)$, and also calculating $u_{ab}^{(\text{QM})}$ only requires $N_Q$ solutions to the MM induction equations, as opposed to $O(N_{AO}^2)$ for IEDRF.

Because the DRF method introduces corrections to two-electron interactions, care must be taken when combining DRF with density functional theory.¹⁵ Since the two-electron interaction is changed in a non-local way, formally the entire density functional changes (and would also now be non-local), however to use existing density functional approximations without reparametrization we simply include a 100% contribution of exact DRF exchange, together with the DRF contribution coulomb term.²¹ The DRF two-electron contribution to the total electron energy is typically small compared to other terms, which justifies the use of this scheme, although it may be possible to develop alternative approaches to combining DFT and DRF. For example it has been proposed that optimal tuning of range-separated hybrids can be used to approximately account for changes in the exchange-correlation functional.¹⁵

### C. Periodic boundary conditions

In order to obtain a method that treats periodic boundary conditions rigorously, we simply replace the standard Coulomb interactions in the above expressions $\phi(\mathbf{r}) = 1/|\mathbf{r}|$ with the Ewald potential for a periodic box³⁰

$$ \phi(\mathbf{r}) = \sum_{\mathbf{n}=0} \frac{1}{|\mathbf{r} + \mathbf{n}\mathbf{L}|} \text{erfc}(\beta|\mathbf{r} + \mathbf{n}\mathbf{L}|) + \frac{4\pi}{L^3} \sum_{\mathbf{k} \neq 0} \frac{e^{-|\mathbf{k}|^2/4\beta^2}}{|\mathbf{k}|^2} e^{i\mathbf{k} \cdot \mathbf{r}} - \frac{\pi}{L^3 \beta^2}. \quad (24) $$

On using this interaction in the previous expressions for the electrostatic energy [Eq. (3), Eq. (4), Eq. (9)], we obtain the PBC polarisation energy and the PBC version of the DREEM interaction Hamiltonian. This energy expression has essentially the same form as the non-PBC DREEM energy, only now there is an additional interaction between the charge density and its periodic images. If we further apply the multipole expansion for the QM charge density, we obtain essentially the same expression as Eq. (12) for the polarisation energy of a set of atom-centred multipoles $q$ with periodic boundary conditions, but where the interactions are now given by the periodic Ewald potential. By using ESPF multipole operators constrained to reproduce the total charge (and dipole moment) of the $\rho_{\mu\nu}(\mathbf{r})$ charge densities, we rigorously ensure the accuracy of the method for long-range interactions up to $O(R^{-1})$, where $l$ is the order of the multipole expansion. In the SI we analyse the long range error in the ESP arising from use of ESPF operators including up to dipoles ($l = 1$), and we conclude that the error introduced is likely within an order of magnitude of typical truncation error in using the PME method for long-range electrostatics.²⁵ Modification of the ESPF operators by constraining them to also reproduce higher moments of the electron density may offer a strategy for further reducing this error in future implementations.²⁸

The only additional terms that arise in the DREEM interaction Hamiltonian are additional contributions to quadratic interaction term $u_{ab}^{(\text{QM})}$ which appear due to interactions of the ESPF multipoles with their periodic replicas. These additional
terms are given by,

$$ u_{ab,\text{PBC}}^{(\text{QM})} = \int dr \int dr' \rho_a(r) \times \left( \phi(r-r') - \frac{1}{|r-r'|} \right) \rho_b(r'). \quad (25) $$

In order to obtain the DREEM interaction Hamiltonian in periodic boundary conditions, we simply swap out the non-PBC interaction expansion coefficients, $u_a^{(\text{QM/MM})}$ and $u_{ab}^{(\text{QM})}$, for their PBC counterparts in Eq. (18). With this approach to introducing periodic boundary conditions to the DREEM method, the interaction with both the polarisable environment and periodic images is treated consistently at the level of an operator on the QM system involving up to two-electron interactions. Furthermore because the form of the interaction Hamiltonian is identical for PBC and non-PBC calculations, the formal computational scaling for the modified QM calculation is unchanged by introducing periodic boundary conditions. The computational cost of evaluating electrostatic energies and solving the induction equations with the Ewald potential can be reduced to $O(N_{\text{MM}} \log N_{\text{MM}})$ using the particle mesh Ewald method to evaluate these terms, and therefore this PBC implementation can be very efficient for large systems.²⁵

### D. Analytic energy gradients

In order perform calculations like geometry optimizations and molecular dynamics analytic gradients of the QM/MM energy are needed. Because the DRF interaction energy takes the form of an operator on the QM system, this is relatively easy to evaluate. As a concrete case we consider deriving analytic energy gradients for Hartree-Fock (HF) theory,³¹ but the generalisation to DFT is straightforward. We have derived and implemented analytic gradients for linear-response methods like Configuration Interaction Singles (CIS) and TDDFT similarly, although the implementation is more complicated.³²,³³ The total energy, using HF to describe the QM system can be written in terms of the one-electron reduced density matrix $D_{\mu\nu}$, where $\mu$ and $\nu$ label spin-orbitals, as

$$ E = \sum_{\mu\nu} H_{\mu\nu} D_{\mu\nu} + \frac{1}{2} \sum_{\mu\nu} D_{\mu\nu} [(\mu\nu|\mu'\nu') - (\mu\mu'|\nu\nu')] D_{\mu'\nu'} + \sum_{\mu\nu} V_{\mu\nu}^{(1)} D_{\mu\nu} + \frac{1}{2} \sum_{\mu\nu} D_{\mu\nu} [(\mu\nu|\hat{V}^{(2)}|\mu'\nu')] - (\mu\mu'|\hat{V}^{(2)}|\nu\nu')] D_{\mu'\nu'} + V_{\text{QM,nuc}} + E_{\text{MM}} + V^{(0)}. \quad (26) $$

Its derivative with respect to a parameter $\xi$, either a QM or MM atom coordinate, involves two types of terms: those that depend on the derivative of the density matrix $D_{\mu\nu}^{\xi}$, and those that involve derivatives of zero, one and two electron matrix elements.³¹ The terms depending on $D_{\mu\nu}^{\xi}$ can be handled with the standard energy-weighted density matrix approach,³¹

$$ E_W^{\xi} = - \sum_{\mu\nu} W_{\mu\nu} S_{\mu\nu}^{\xi} = - \sum_{\mu\nu} \sum_{i \in \text{occ}} \epsilon_i C_{\mu i} C_{\nu i} S_{\mu\nu}^{\xi}. \quad (27) $$

Derivatives of the isolated QM and MM system energies with $D_{\mu\nu}$ held constant have also been derived previously.³¹ The remaining interaction term derivative with $D_{\mu\nu}$ held constant, which we denote $E_{\text{int}}^{(\xi)}$, can be written in terms of the one-electron density matrix and two-electron density matrix, $\Gamma_{\mu\nu\mu'\nu'} = D_{\mu\nu}D_{\mu'\nu'} - D_{\mu\mu'}D_{\nu\nu'}$, as³¹,³³

$$ E_{\text{int}}^{(\xi)} = \frac{\partial V^{(0)}}{\partial \xi} + \sum_{\mu\nu} D_{\mu\nu} \frac{\partial V_{\mu\nu}^{(1)}}{\partial \xi} + \frac{1}{2} \sum_{\mu\nu\mu'\nu'} \frac{\partial (\mu\nu|\hat{V}^{(2)}|\mu'\nu')}{\partial \xi}. \quad (28) $$

This term can now further be split into terms that depend on derivatives of the multipole operators, $\partial Q_{a\mu\nu}/\partial\xi$, and classical terms that depend on derivatives of the polarisation energy expansion coefficients, $\partial u_a^{(\text{QM/MM})}/\partial\xi$ etc.,

$$ E_{\text{int}}^{(\xi)} = E_{\text{int}}^{(\xi, \text{op})} + E_{\text{int}}^{(\xi, \text{cl})}. \quad (29) $$

The derivatives of the multipole operators are only non-zero for QM atom positions, because the positions of the multipoles and the grid of points used in the ESPF procedure depend on QM atom coordinates. Expressions for the ESPF operator derivatives have been reported previously and are given in the supporting information.³,²² The one-electron contribution to $E_{\text{int}}^{(\xi, \text{op})}$ is

$$ E_{\text{int},1}^{(\xi, \text{op})} = \sum_a \left( u_a^{(\text{QM/MM})} + \sum_b u_{ab}^{(\text{QM})} Z_b \right) \left\langle \frac{\partial Q_a}{\partial \xi} \right\rangle^{(1)} + \frac{1}{2} \sum_{ab} u_{ab}^{(\text{QM})} \left\langle \frac{\partial}{\partial \xi} Q_a Q_b \right\rangle^{(1)} \quad (30) $$

and the two-electron contribution is

$$ E_{\text{int},2}^{(\xi, \text{op})} = \frac{1}{2} \sum_{ab} u_{ab}^{(\text{QM})} \left\langle \frac{\partial}{\partial \xi} Q_{a,1} Q_{b,2} \right\rangle^{(2)}. \quad (31) $$

Here we have introduced the shorthand $\langle A \rangle^{(1)} = \sum_{\mu\nu} D_{\mu\nu} A_{\mu\nu}$ to denote the expectation value of a one-electron operator and $\langle A \rangle^{(2)} = \sum_{\mu\nu\mu'\nu'} \Gamma_{\mu\nu\mu'\nu'} (\mu\nu|\hat{A}|\mu'\nu')$ to denote the expectation value of a two-electron operator. The derivatives of the products of multipole operators can be obtained with a simple application of the chain rule. Note also that the two-electron density matrix $\Gamma_{\mu\nu\mu'\nu'}$ never has to be constructed explicitly, because the two-electron expectation values can be decomposed into contractions of lower rank tensors.

The final classical derivative term (which includes the $\partial V^{(0)} / \partial \xi$ term above) can be written as

$$ E_{\text{int}}^{(\xi, \text{cl})} = \sum_a \frac{\partial u_a^{(\text{QM/MM})}}{\partial \xi} \langle q_a \rangle + \frac{1}{2} \sum_{ab} \frac{\partial u_{ab}^{(\text{QM})}}{\partial \xi} \langle q_a q_b \rangle \quad (32) $$

where we have introduce the expectation value of a general operator involving up to two-electron terms, $\hat{A} = A^{(0)} + \hat{A}^{(1)} + \hat{A}^{(2)}$, as

$$ \langle A \rangle = A^{(0)} + \langle A^{(1)} \rangle^{(1)} + \langle A^{(2)} \rangle^{(2)}, \quad (33) $$
so $\langle q_a \rangle$ and $\langle q_a q_b \rangle$ are

$$ \langle q_a \rangle = Z_a + \langle Q_a \rangle^{(1)} \quad (34) $$

$$ \langle q_a q_b \rangle = Z_a Z_b + Z_a \langle Q_b \rangle^{(1)} + Z_b \langle Q_a \rangle^{(1)} + \langle Q_{a,1} Q_{a,2} \rangle^{(2)}. \quad (35) $$

From this expression it would seem that we need to evaluate all of the $\partial_{\xi} u_{ab}^{(QM)}$ derivatives, which would require $O(N_{QM}^2)$ solutions of the induction equations. An alternative is to explicitly differentiate the expression for $u_{ab}^{(QM)}$, but this can be complicated to implement in existing MM codes. Instead we take an approach that only requires routine energy and force calculations for the polarizable MM system.

We start by noting that the derivatives of $u_a^{(QM/MM)}$ can be calculated by evaluating the energy derivatives for the polarizable MM system with a test multipole $\pm q_{a,0}$ at QM site a and all others are zero. $q_{a,0}$ is a constant and its value can be chosen arbitrarily. Denoting this vector of multipoles as $\pm \mathbf{q}_0^a$, the energy expansion term can be evaluated as

$$ u_a^{(QM/MM)} = \frac{1}{2q_{a,0}} (E_{\text{pol}}(+\mathbf{q}_0^a) - E_{\text{pol}}(-\mathbf{q}_0^a)) \quad (36) $$

and the derivative can therefore be evaluated as

$$ \frac{\partial u_a^{(QM/MM)}}{\partial \xi} = \frac{1}{2q_{a,0}} \left( \frac{\partial}{\partial \xi} E_{\text{pol}}(+\mathbf{q}_0^a) - \frac{\partial}{\partial \xi} E_{\text{pol}}(-\mathbf{q}_0^a) \right). \quad (37) $$

This requires $2N_Q$ solutions to the induction equations and force evaluations.

The derivatives of the quadratic terms can also be evaluated with only $N_Q$ additional force calculations. To see this we note that the quadratic term in the interaction energy can be written in terms of the symmetric multipole correlation matrix,

$$ [\mathbf{Q}_2]_{ab} = \frac{1}{2} (\langle q_a q_b \rangle + \langle q_b q_a \rangle). \quad (38) $$

This symmetric matrix can be decomposed into its eigenvalues and eigenvectors

$$ \mathbf{Q}_2 = \sum_n \lambda_n \mathbf{v}_n \mathbf{v}_n^\top, \quad (39) $$

and then the interaction energy can be written as

$$ E_{\text{int}} = \sum_a \langle q_a \rangle u_a^{(QM/MM)} + \frac{1}{2} \sum_{nab} \lambda_n u_{ab}^{(QM)} v_{na} v_{nb} \quad (40) $$

$$ = \sum_a u_a^{(QM/MM)} \langle q_a \rangle + \sum_n \lambda_n E_{\text{pol}}(\mathbf{v}_n) \quad (41) $$

$$ - \sum_n \lambda_n (u^{(MM)} + \sum_a u_a^{(QM/MM)} v_{na}) $$

Differentiating this (while holding the electronic state and multipole moment operators constant) yields the classical energy derivative term,

$$ E_{\text{int}}^{(\xi, \text{cl})} = \sum_a \frac{\partial u_a^{(QM/MM)}}{\partial \xi} \langle q_a \rangle + \sum_n \lambda_n \frac{\partial}{\partial \xi} E_{\text{pol}}(\mathbf{v}_n) - \sum_n \lambda_n \left( \frac{\partial u^{(MM)}}{\partial \xi} + \sum_a \frac{\partial u_a^{(QM/MM)}}{\partial \xi} v_{na} \right). \quad (42) $$

<figure><img src="image_2.png" /><figcaption>FIG. 1. A schematic of the computational workflow implemented in the OpenESPF code for calculating DREEM energies and forces.</figcaption></figure>

This only requires an additional $N_Q$ force calculations, and thus overall the classical contribution to the analytic energy derivative only requires $O(N_Q)$ solutions to the MM polarizable force calculations. This means both energy and force calculations only require a linear scaling number of MM force evaluations with the QM system size. Furthermore this approach to evaluating the classical energy derivatives is identical for both PBC and non-PBC calculations.

Combining all of these terms the full energy gradient is

$$ E^{\xi} = E_{W}^{\xi} + E_{QM}^{(\xi)} + E_{MM}^{(\xi)} + E_{\text{int}}^{(\xi, \text{op})} + E_{\text{int}}^{(\xi, \text{cl})}. \quad (43) $$

The analytic energy derivative expression is identical for DFT methods, where there is simply an additional exchange-correlation functional contribution in $E_{QM}^{(\xi)}$. For TDDFT methods the analytic gradient can also be written in this general form, only now the $E_W^{\xi}$ term involves additional terms from solving the coupled-perturbed equations for the Z vector, which describes orbital relaxation, and the expressions for excited state relaxed one- and two-electron density matrices are more complicated, involving the ground-state density matrix, excitation amplitudes and the Z vector. Also as in SCF methods the full rank-four two-electron density matrices never need to be constructed in full for TDDFT because expectation values can be evaluated more efficiently through contractions of lower-rank tensors.

E. Implementation Details

We have implemented the DREEM method in the OpenESPF code.²⁹ This provides a Python interface between the PySCF code³⁴⁻³⁶ for electronic structure quantum calculations, and the OpenMM code³⁷,³⁸ for performing MM calculations for the polarization energy and other MM energy terms. The OpenESPF code uses as its core class a QM/MM system class, which calls methods in separate QM and MM system objects which in turn call PySCF and OpenMM respectively for performing the core calculations such as solving MM induction equations and solving the QM self-consistent field equations. The QM/MM system class passes the minimum information
needed between QM and MM systems for DREEM calculations. The general workflow for a DRF energy and force calculation in OpenESPF is depicted in Fig. 1. Firstly QM and MM coordinates are passed to the MM system, where OpenMM is called to calculate the polarization energy expansion terms, $u^{(MM)}$, $u_a^{(QM/MM)}$, and $u_{ab}^{(QM)}$. The polarization energy expansion is then passed to the QM system, which constructs the ESPF operators (handled by a separate Multipole class) and modifies the PySCF one-electron Hamiltonian and routines for calculating Coulomb and exchange matrices, and calls PySCF routines for calculating electronic energies. The resulting density matrix is then used to calculate QM electronic energy derivative terms, also by calling PySCF routines, and the ESPF operator derivative contributions to the gradient, as well as expectation values of the atom-centred multipole moments. These are then passed to the MM system where the classical interaction energy derivatives are calculated, and final forces are then computed by combining all these terms. Because the OpenESPF code implements a framework where minimal information is needed to be passed between QM and MM codes, in principle other QM and MM codes can be interfaced with relative ease.

In addition to the DREEM interaction, we also include the exchange model of Pauli repulsion to account explicitly for how QM electrons are repelled by electrons in the MM region (which are only treated implicitly). This model is based on the empirical observation that the total repulsion energy between two molecules is approximately $-1/2$ of the exchange energy, so Pauli repulsion can be approximated by an additional repulsive one-electron term in the Hamiltonian given by$^{39–41}$

$$V_{\mu\nu}^{\text{rep}} = \frac{1}{2} \sum_{\mu'\nu'} D_{\mu'\nu'}^{\text{MM}} (\chi_{\mu}\varphi_{\mu'}^{\text{MM}}|\varphi_{\nu'}^{\text{MM}}\chi_{\nu}) \quad (44)$$

where $D_{\mu'\nu'}^{\text{MM}}$ is an effective model density matrix for the MM atoms, and $\varphi_{\mu'}^{\text{MM}}$ is a model orbital on MM atoms. There are various ways to parametrize the MM orbitals and density matrices. We use a simple scheme where the orbitals are parametrized as a STO-3G approximation to a 1s orbital with the width adjusted to agree with density derived van der Waals' radii for the atoms, and the density matrix elements are treated as diagonal with a population given by the valence shell population adjusted to be consistent with the total charge on the atom.$^{21}$ This repulsion is very short ranged so can be ignored for MM atoms far from the QM region. Analytic derivatives of this repulsion energy are easy to derive because they only involve standard two-electron derivative integrals, e.g. $((\nabla\chi_{\mu})\varphi_{\mu'}^{\text{MM}}|\varphi_{\nu'}^{\text{MM}}\chi_{\nu})$ and $(\chi_{\mu}(\nabla\varphi_{\mu'}^{\text{MM}})|\varphi_{\nu'}^{\text{MM}}\chi_{\nu})$. The analytic derivatives of this repulsion term are also implemented in OpenESPF.

In all calculations of the polarization energy between multipoles in the QM and MM region the Thole damping scheme is used, where interactions with induced dipoles are damped to avoid divergence in the polarization energy at short range.$^{6,42}$ We also damp the interactions between permanent multipoles in the MM region and multipoles in the QM region using the same scheme. Thole damping requires assigning an effective

atomic polarizability $\alpha_A$ and thole parameters $a_A$ to define the pairwise damping radius $r_{AB}^{-1} = \min(a_A, a_B)/(\alpha_A\alpha_B)^{1/6}$ to each QM atom, which we assign by selecting values from the AMOEBA force field for chemically similar atoms.$^6$ By damping the QM charge- MM charge interactions, we miss some charge-penetration effects,$^{43}$ which we account for by adding an extra damped repulsion term of the form,

$$V_{\text{rep,cp}} = \sum_{A \in \text{QM}} \sum_{B \in \text{MM}} \frac{Z_A Z_{B,\text{eff}}}{\|\mathbf{R}_A - \mathbf{R}_B^{\text{MM}}\|} f_{\text{cp}}^{\text{AB}} (\|\mathbf{R}_A - \mathbf{R}_B^{\text{MM}}\|), \quad (45)$$

where $Z_A$ is the nuclear charge associated with atom A, $Z_{B,\text{eff}} = Z_B - N_{\text{core},B}$ is the nuclear charge on MM atom B with the number of core electrons removed, and $f_{\text{cp}}^{\text{AB}}(R)$ is a damping factor that approaches 1 at small R and 0 very rapidly at large R. The form of this damping is simply $1 - f_{\text{damp}}^{\text{AB}}(R)$ for the Thole damped charge-charge interactions, so no additional parameters are needed for this charge-penetration repulsion term. Further justification of this term is given in the supporting information.

The OpenESPF code currently supports efficient energy and gradient calculations for spin-restricted, unrestricted, and restricted-open SCF methods, including pure, hybrid and range-separated hybrid density functionals. The ESPF operators are currently limited to atom centered charges and dipoles, but future extensions to include higher order multipoles and virtual charge sites is planned. Furthermore efficient single-point energies and gradients are implemented for spin-restricted TDDFT methods, both with and without the Tamm-Dancoff approximation. Use of density-fitting for the QM calculation together with DREEM is also supported. Single point energy calculations with almost any post-SCF method available in PySCF, is also supported through in-core calculation and storage of the modified two-electron interaction integrals, including CASSCF, MP2, and CCSD(T) methods. GPU acceleration of MM energy calculations is possible through the use of OpenMM. Several examples are included including extensive validation of the analytic energy gradients.

### III. RESULTS

Having described how to combine the DREEM method with PBC, and how to obtain analytic derivatives of the total QM/MM energy, we now illustrate the method and its implementation in the OpenESPF code with applications to small organic molecules in aqueous solution.

#### A. Excitation energies in PBC

As a first example of the DREEM method in periodic boundary conditions as implemented in the OpenESPF code, we have calculated the excitation energies for acrolein in water. The coordinates for acrolein in a PBC water box were taken from our previous study$^{21}$ where they were sampled using a bespoke
<figure><img src="image_3.png" /><figcaption>FIG. 2. (A) Illustration of the droplet and PBC models for acrolein in water. B and C) excitation energies as a function of cut-off distance for the droplet model (red squares) and real-space cut-off for the PME PBC model (blue circles) for the S₁ and S₂ states of acrolein. C) Ground state acrolein dipole moment as a function of cut-off distance.</figcaption></figure>

parameterized force-field for acrolein based on the OPLS/AA force field⁴⁴,⁴⁵ together with the TIP3P water model. Other computational details for sampling can be found in Ref. 21. As a demonstration of the utility of the PBC approach, we compare the first two excitation energies and the ground-state dipole moment of acrolein as a function of both the real-space cut-off distance in a PME-based PBC calculation, and as a function of cut-off distance for a finite, non-PBC droplet model. The QM region included only the acrolein molecule, whilst water molecules were treated at the MM level with the polarizable model from Ref. 46 with Thole damping. This is illustrated in Fig. 2A. The PBC water box contained 5000 water molecules, with a total of 15000 polarizable sites (three per water) and a periodic box dimension of 52.937Å was used. In all calculations the first three excited state energies of acrolein are calculated with TDA-TDDFT⁴⁷ using the ωB97X-D3 functional⁴⁸ and pc-1 (double zeta) basis set. Atom-centred charges and dipoles are used in the QM region ESPF operators. The Ewald α parameter and k-space cut-off was determined automatically determined by OpenMM for the PBC calculations. The ESPF representation of the QM region used atom-centred charges and dipoles for each QM atom.

In Fig. 2B and Fig. 2C we show the calculated excitation energies for the S₁ and S₂ state for the PBC and finite-sized droplet calculations as a function of cut-off distance, as well as the vacuum excitation energy (where the same geometry for the acrolein molecule is used as in the QM/MM calculations). First we see that the PBC excitation energies are effectively independent of the real-space cut-off. This is because the full periodic system is still treated for all values of the real-space cut-off, so as long as an appropriate Ewald α parameter and

k-space cut-off are used results should be independent of this cut-off. We see that the droplet QM/MM excitation energies are relatively slow to converge with respect to droplet size, with errors of about 0.1 eV for a 10 Å cut-off and even with a 22 Å cut-off excitation energies do not exactly agree with the PBC calculation. The same behaviour can be seen for the ground state dipole moment as a function of cut-off shown in Fig. 2D. Here we only show results for a single set of QM and MM atom coordinates, but we have observed very similar behaviour for other snapshots from the same MD trajectory.

In Fig. 3 we show the computational timings for the different components of the DREEM calculations as a function of cut-off distance. All calculations were run on an Apple MacBook Pro with M2 processor and 8GB of memory, and the OpenCL platform was used for the OpenMM calculations. We see that in both PBC and droplet calculations the QM calculation time is approximately independent of the cut-off distance, and the additional terms in the DREEM calculation only add about computational 10% overhead compared to the vacuum calculation. The independence of cut-off is expected because the embedding terms only enter into the QM calculation through the polarization response terms $u_a^{(QM/MM)}$ and $u_{ab}^{(QM)}$ which scale with the size of the QM system. The evaluation of these response terms however does depend on the cut-off distance. The droplet MM energy calculations are considerably faster than the PBC calculations, primarily because the droplet calculations involve fewer MM atoms, 5745 in the largest droplet calculation compared to 15000 in all the PBC calculations, plus the PBC calculations require additional computational overhead associated with the PME calculation. For an appro-
priately chosen cut-off distance (in this case 10Å) this overhead for PBC calculations can be minimised, and the fastest PBC calculation time is only about 2.4 times slower than the largest droplet calculation, for which calculated properties have still not fully converged. The non-monotonic dependence of the MM PBC calcualtion time is due to different reciprocal space cut-offs assigned automatically by OpenMM for the differ- ent real-space cut-off distances; for small real-space cut-offs a large reciprocal space cut-off is automatically selected, but for large cut-offs there are more real-space interactions to be evaluated. The computation times for the QM and MM com- ponents of the DREEM calculation will of course be system and computing platform specific, but these results serve to il- lustrate that the DREEM components do not add considerable overhead to the QM calculation, despite modifying both one and two-electron terms, and that PBC calculations are only moderately more expensive than droplet calculations, whilst offering a more realistic description of condensed phase envi- ronments with long-range electrostatics.

B. Analytic gradient validation

Having demonstrated the potential utility of our PBC
DREEM implementation, we move on to validate the analytic
gradient implementation in OpenESPF. As a test system we
consider a system consisting of a single QM water molecule
embedded in a periodic box of 299 AMOEBA 2013 water
molecules. The configuration was taken from an equilibrated
NPT molecular dynamics run using the AMOEBA 2013 force
field with 300 water molecules at 298K and 100 kPa. In order
to validate the analytic gradients we consider the x force as
a function of x displacement for the QM water O atom and
the nearest MM water H atom, which lies 1.9 Å from the QM

<figure><img src="image_4.png" /><figcaption>FIG. 3. Timings for the QM and MM components of the DREEM calculation for acrolein in a PBC box of 5000 water molecules as a function of real-space cut-off in the evaluation of MM energy terms (blue lines) and acrolein in finite sized water droplets as a function of droplet radius (red lines).</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>FIG. 4. Panels A) and C) show interaction forces for QM and MM atoms for a QM/MM simulation of a periodic water box containing 300 water atoms, where one is treated as the QM system. Panels B) and D) show force errors for the interaction and QM or MM forces calculated with different OpenMM platforms which work in either double-precision (Reference) or GPU-accelerated single-precision (OpenCL).</figcaption></figure>

water O atom. The QM atom is described at the HF/6-31G* level of theory with an SCF energy convergence of 10<sup>-10</sup> Ha, and atom centred ESPF multipoles up to dipoles are included to describe the QM/MM interaction. The forces are calculated with a fourth order central finite difference approximation with a step size of 0.015 Å ≈ 0.028 Bohr. In each case the Particle Mesh Ewald (PME) method is used to describe the periodic electrostatic interactions with fourth order spline interpolation with β = 3.29 nm<sup>-1</sup> and k<sub>max</sub> = 49 and a real-space cut-off of 1.0 nm (as automatically determined by OpenMM to achieve a desired force error of 10<sup>-5</sup> kJ mol<sup>-1</sup> nm<sup>-1</sup>), and the dipole convergence was set to 10<sup>-6</sup> e nm.

In Fig. 4A we show the numerical and analytic inter-
action force on the MM water H atom, where we define
F<sub>int</sub> = F<sub>tot</sub> − F<sub>MM</sub> where F is the total force and F<sub>MM</sub> is
the force in the absence of the QM system, where we see ex-
cellent agreement between numerical and analytical results.
These are calculated using OpenMM’s double-precision Ref-
erence platform. The error is further analysed in Fig. 4B,
where we compare the error in the interaction force (solid
lines) to the numerical error in the MM (dashed lines), us-
ing both OpenMM’s Reference (CPU) platform, which runs
in double precision, and the OpenCL (GPU) platform which
runs in single precision. For the Reference double-precision
platform, we see very similar error for the MM and interaction
forces at about 10<sup>-8</sup> Ha/Bohr, whereas for the single-precision
OpenCL platform the interaction force error is typically around
10<sup>-5</sup> Ha/Bohr, about an order of magnitude higher than the
MM force error calculated with the same platform. This is
<figure><img src="image_6.png" /><figcaption>FIG. 5. Energy change during a DREEM simulation of one QM water molecule in a periodic box of 99 AMOEBA 2013 water molecules, with different time-steps and using either A) OpenMM's Reference double precision platform or B) the OpenCL single precision platform.</figcaption></figure>

likely because the interaction force involves combining $N_Q$ MM force calculations, so the error is compounded for the interaction force.

The numerical and analytic interaction force calculation for the QM water O atom is shown in Fig. 4C, again using the Reference platform, where excellent agreement between analytic and numerical forces can be seen (in this case interaction force is defined as the total force with the QM contribution in the absence of the MM system subtracted $F_{int} = F_{tot} - F_{QM}$). Fig. 4C we plot the errors for the QM energy calculation (using PySCF's double-precision CPU implementation) and interaction energy with the Reference double-precision platform, where we see very similar error of around $10^{-7}$ Ha/Bohr. As for the MM interaction force, using the single-precision OpenCL platform deteriorates the quality of the force, increasing the error to around $10^{-5}$ Ha/Bohr.

As further validation of analytic forces, we check the energy conservation for an NVE molecular dynamics simulation of a QM water molecule in a periodic box of 99 AMOEBA 2013 water molecules using the DREEM method. The initial configuration is taken from an NPT MM simulation at 298K and 100 kPa pressure, and velocities are sampled from a Boltzmann distribution at 298K, and Newton's equations of motion for the system are integrated using the velocity Verlet algorithm. The PBC box dimension was 15.9 Å and again the PME method is used for electrostatics with a real-space cut-off of 0.75 nm, with $\beta = 4.39 \text{ nm}^{-1}$ and $k_{max} = 47$. The induced dipole convergence threshold was set to $10^{-6}$ e nm and the QM SCF energy convergence threshold was set to $10^{-10}$ Ha and atom-centred charges are used for the ESPF operators. Overall 20 ps of NVE dynamics were run with two different time-steps: 0.5 fs and 1.0 fs.

In Fig. 5A we show the change in total energy in meV/atom for the simulation using the Reference double-precision platform for the MM calculations. Overall we see good energy conservation, with energy fluctuations reducing by about a factor of 4 on reducing the time-step by a factor of 2, as expected for the velocity Verlet algorithm with a global error of $O(\Delta t^2)$. The overall energy drift from a linear fit is small at -0.001 meV/atom/ps for $\Delta t = 0.5$ fs and -0.006 meV/atom/ps for $\Delta t = 1.0$ fs. We have also repeated these simulations with the same initial positions and velocities using the single-precision OpenCL platform. These calculations are significantly faster, but the use of single-precision reduces the quality of the forces, as shown above. Reassuringly we still see reasonable energy conservation, with short time fluctuations similar to the use of double-precision, but with larger, longer time-scale fluctuations on top of these fluctuations. Additionally the energy drifts are larger, at +0.006 meV/atom/ps for $\Delta t = 0.5$ fs and +0.019 meV/atom/ps for $\Delta t = 1.0$ fs. We note in passing that when we run these simulations just using the MM force field and the same time-steps and platforms, we see near identical energy fluctuations and drifts. The choice of between using the slower, more accurate Reference platform and a lower accuracy GPU-accelerated single-precision platform will depend on the desired accuracy for the problem at hand.

### C. Excited state MD and fluorescence spectra

As a final demonstration of the DREEM method, and its implementation in the OpenESPF code, we have used DREEM to assess the impact of non-Condon and dynamical effects on the fluorescence spectrum of acetone in aqueous solution, and we compare this to the gas phase fluorescence spectrum. In this example we also use the Amoeba 2013 water model for the MM water molecules and PBE0/def2-SV(P) for the QM acetone molecule, with density fitting using the def2/JK auxiliary basis set. Atom-centred ESPF charge operators are used to describe the QM/MM interaction. The $S_1$ state energies and forces are obtained using TDA-TDDFT. Starting from the same acetone MM equilibrated geometry used in the previous example, we first ran 50 steps of energy minimization for the $S_1$ excited state of acetone with the DREEM method using the L-BFGS-B method. From this geometry we ran 12 ps of excited state molecular dynamics using the fast-forward Langevin algorithm$^{49}$ in order to obtain an equilibrated starting geometry at 298K. From this geometry velocities were resampled from a Boltzmann distribution and fast-forward Langevin dynamics were run for a further 10 ps for production. In all fast-forward Langevin dynamics simulations a time-step of 1 fs and damping constant of $10 \text{ ps}^{-1}$ were used. At every time-step the $S_0-S_1$ excitation energy and transition dipole moment were calculated. Using the procedure described in Ref. 50 the fluorescence spectra including non-Condon effects were calculated. In short, the energy gap and transition dipole moment time-correlation functions were calculated and Fourier transformed to give a set of spectral densities and from this the fluorescence spectrum with quantum vibronic effects were calculated both with and without non-Condon effects,
<figure><img src="image_7.png" /><figcaption>FIG. 6. Joint probability distributions for the square transition dipole moment and S₁-S₀ energy gap for A) gas phase acetone and B) acetone in water calculated with PBC DREEM for the aqueous model. Both distributions come from sampling the S₁ excited state potential energy surface.</figcaption></figure>

essentially by mapping onto a displaced quantum harmonic oscillator model. For further details we refer the reader to Ref. 50. The transition moments were all rotated into a molecular frame defined by minimising mean square displacement of the non-hydrogen atom positions in the QM region for this procedure.

In Fig. 6 we show the joint probability distributions of the square transition dipole moment $|\mu_{01}|^2$ for thr S₁ → S₀ transition for acetone in the gas phase (A) and aqueous solution (B), calculated with DREEM for the aqueous phase. These show a clear negative correlation between the transition dipole moment and the energy gap in both gas and solution phases, so configurations with smaller energy gaps have larger transition probabilities. The maximum in the joint distribution shifts slightly from the gas to solution phase from a region with higher $|\mu_{01}|^2$ and lower energy gap, to a region with a larger energy gap and lower transition moment. These differences are also reflected in the mean energy gap, and the mean transition-moment weighted energy gap for the two distributions, as show in Table I. On accounting for the $|\mu_{01}|^2$ weighting the mean transition energy shifts by about 0.07 eV in the gas phase and 0.10 eV in the solution phase.

To further account for both dynamical and non-Condon effects in the fluorescence spectra, we have calculated the fluorescence spectra using our excited state molecular dynamics simulations together with the Gaussian-Condon (GC) and Gaussian-Non-Condon (GNC) theories in Ref. 50. In Fig. 7A we show the calculated spectral densities, $J(\omega)$, which describe the vibrational frequency dependence of the harmonic model reorganisation energy. The peak between 1300 and 1500 cm⁻¹ in both the gas and solution phases arises from the C–O bond stretching in the electronic excited state, where it is partially broken, hence it appears at a lower frequency than the ground state C=O stretch in acetone at about 1700 cm⁻¹.

<table><caption>TABLE I. Mean Condon excitation energies $\frac{1}{N} \sum_{r=1}^{N} \Delta E^{(r)}$, mean non-Condon excitation energies $\frac{\sum_{r=1}^{N} |\mu_{01}^{(r)}|^2 \Delta E^{(r)}}{\sum_{r=1}^{N} |\mu_{01}^{(r)}|^2}$, and fluorescence maxima for GC and GNC theories in the gas and aqueous solution phases calculated in this work. Aqueous solution results come from PBC DREEM simulations. Experimental values from Ref. 55.</caption><thead><tr><th></th><th>Gas</th><th>Aqueous</th><th>Shift</th></tr></thead><tbody><tr><td>Mean (Condon)</td><td>3.289 ± 0.016</td><td>3.336 ± 0.020</td><td>0.049</td></tr><tr><td>Mean (Non-Condon)</td><td>3.220 ± 0.018</td><td>3.237 ± 0.020</td><td>0.017</td></tr><tr><td>Max (Condon)</td><td>3.356 ± 0.016</td><td>3.403 ± 0.020</td><td>0.047</td></tr><tr><td>Max (Non-Condon)</td><td>3.335 ± 0.016</td><td>3.372 ± 0.020</td><td>0.037</td></tr><tr><td>Experiment</td><td>3.00 – 3.01</td><td>3.06</td><td>0.05 – 0.07</td></tr></tbody></table>

This peak is slightly red-shifted in solution due to interactions with the solvent. The larger contributions to the spectral density in the low frequency, between 0 and 100 cm⁻¹, region in the aqueous solution spectral density arises from the solvent reorganisation energy, where the polar water solvent interacts differently with the ground state and the $n\pi^*$ excited state charge distributions.

In Fig. 7B we show the calculated spectra in both gas and aqueous phases, with and without non-Condon effects. In both phases non-Condon effects slightly broaden and red-shift the spectra, shifting the maximum by about 0.02 and 0.03 eV in gas and solution phases respectively, which is a slightly smaller shift than that predicted by the simple transition-moment weighted average transition energy if 0.07-0.10 eV. The aqueous spectra are slightly broadened compared to the gas phase due to additional solvent contributions to the reorganisation energy. Including quantum vibronic effects through the GC and GNC theories blue shifts the spectra compared to spectra predicted with a simple classical average of the transition energy. This is because a significant portion of the reorganisation energy is associated with high frequency vibrational modes, so nuclear quantum effects become important. The calculated gas and solution phase spectra are both quite significantly blue shifted compared to the experimental spectra, in both phases by about 0.33 eV, which is primarily due to the use of the relatively low-level TDA/PBE0/def2-SV(P) electronic structure method. Comparing the TDA/PBE0/def2-SV(P) excitation energy at the minimum energy geometry for the S₁, 3.21 eV, to that predicted with the higher level DLPNO-STEOM-CCSD/aug-cc-pVQZ⁵¹,⁵² (calculated with Orca 6.0.0⁵³) at the same geometry, 2.91 eV, suggests a blue-shift of about 0.30 eV in our calculated spectra is to be expected. Previous studies of acetone S₁ fluorescence using non-polarizable QM/MM have predicted larger gas to aqueous solvatochromic shifts,⁵⁴ whereas other studies accounting for state-specific solvent polarization predicted shifts much closer to our own predictions. Upon transition from the $n\pi^*$ state to the ground state, charge density moves from the carbonyl carbon atom back to the oxygen atom, and non-poarizable models cannot account for the instantaneous solvent response to the change in acetone charge density, so excitation energy shifts will tend to be overestimated.
<figure><img src="image_8.png" /><figcaption>FIG. 7. A) Spectral densities for the gas (grey dashed) and aqueous solution (blue solid line) fluorescence calculated from excited state molecular dynamics. B) Gas and aqueous solution fluorescence spectra calculated with GC and GNC theories for acetone, together with the gas phase spectrum calculated at the Franck-Condon point with a Gaussian broadening of 0.3 eV. Spectra are normalised to have the same integral. Aqueous solution results come from PBC DREEM simulations.</figcaption></figure>

## IV. CONCLUSIONS

In this work we have presented an extension to the DREEM method to include periodic boundary conditions and analytic energy gradients, enabling QM/MM simulations of molecules in condensed phases. The direct reaction field framework that DREEM employs facilitates the inclusion of excited state specific environment polarization and dispersion interactions, and the use of ESPF multipole operators gives a computationally efficient description of QM/MM interactions. By combining DREEM with periodic boundary conditions, we can now perform simulations of condensed phases without artefacts that accompany simple finite droplet models, whilst also using the efficient particle mesh Ewald method to accelerate calculations of electrostatic interactions. The analytic energy gradient implementation we have described has enabled us to perform molecular dynamics simulations with DREEM for both ground and electronic excited states of molecules in polarizable MM environments. These developments are implemented in the open-source OpenESPF code, which interfaces the PySCF quantum chemistry code with the OpenMM code for MM calculations. Previous studies have shown that environment electronic polarization can have a significant impact of excited state energies and properties calculated in QM/MM simulations, and neglecting these effects can in some cases give rise to bigger errors than those originating from the approximate electronic structure methods being used. Overall the DREEM method offers a good balance between efficiency and accuracy when accounting for environment polarization effects in QM/MM simulations, and the development of analytic gradients and a PBC-compatible version of the method presented here paves the way for simulations of more complex systems, particularly for electronic excited states.

## Acknowledgements

This work was supported by "Agence Nationale de la Recherche" through the project MAPPLE (ANR-22-CE29-0014-01).

## Supporting information

Supporting information provides additional technical details around the DREEM method: ESPDF operators, Thole damping and the AMOEBA force field. Examples of numerical vs analytic forces are also provided, and force calculation timings.

## Data availability

The OpenESPF code, provided for free at www.github.com/tomfay/openespf. Examples of for performing typical calculations are also provided in the GitHub repository, together with example scripts for performing calculations like those performed in this paper. Data not provided in the manuscript are available from the corresponding author upon a reasonable request.

## Code availability

The OpenESPF code, provided for free at www.github.com/tomfay/openespf. Examples of for performing typical calculations are also provided in the GitHub repository.

## Conflicts of interest

The authors declare no conflicts of interest.

<footnote n="1">¹H. M. Senn and W. Thiel, "QM/MM Methods for Biomolecular Systems," Angewandte Chemie International Edition **48**, 1198–1229 (2009).</footnote>

<footnote n="2">²C. E. Tzeliou, M. A. Mermigki, and D. Tzeli, "Review on the QM/MM Methodologies and Their Application to Metalloproteins," Molecules **27**, 2660 (2022).</footnote>

<footnote n="3">³S. Bonfrate, N. Ferré, and M. Huix-Rotllant, "Analytic Gradients for the Electrostatic Embedding QM/MM Model in Periodic Boundary Conditions Using Particle-Mesh Ewald Sums and Electrostatic Potential Fitted Charge Operators," Journal of Chemical Theory and Computation **20**, 4338–4349 (2024).</footnote>

<footnote n="4">⁴S. Bonfrate, W. Park, D. Trejo-Zamora, N. Ferré, C. Ho Choi, and M. Huix-Rotllant, "Assessment of Free Energies From Electrostatic Embedding Density Functional Tight Binding-Based/Molecular Mechanics in Periodic Boundary Conditions," Journal of Computational Chemistry **46**, e70107 (2025).</footnote>

<footnote n="5">⁵K. Csizi and M. Reiher, "Universal qm/mm approaches for general nanoscale applications," WIREs Computational Molecular Science **13**, e1656 (2023).</footnote>

<footnote n="6">⁶P. Ren and J. W. Ponder, "Polarizable Atomic Multipole Water Model for Molecular Mechanics Simulation," The Journal of Physical Chemistry B **107**, 5933–5947 (2003).</footnote>

<footnote n="7">⁷Y. Shi, Z. Xia, J. Zhang, R. Best, C. Wu, J. W. Ponder, and P. Ren, "Polarizable Atomic Multipole-Based AMOEBA Force Field for Proteins," Journal of Chemical Theory and Computation **9**, 4046–4063 (2013), publisher: American Chemical Society (ACS).</footnote>
⁸M. Bondanza, M. Nottoli, L. Cupellini, F. Lipparini, and B. Mennucci, "Polarizable embedding QM/MM: the future gold standard for complex (bio)systems?" Physical Chemistry Chemical Physics 22, 14433-14448 (2020).

⁹M. Bondanza, T. Nottoli, M. Nottoli, L. Cupellini, F. Lipparini, and B., "The OpenMMPol library for polarizable QM/MM calculations of properties and dynamics," The Journal of Chemical Physics 160, 134106 (2024).

¹⁰A. Marini, A. Munoz-Losa, A. Biancardi, and B. Mennucci, "What is Solvatochromism?" The Journal of Physical Chemistry B 114, 17128-17135 (2010).

¹¹L. Nicoli, T. Giovannini, and C. Cappelli, "Assessing the quality of QM/MM approaches to describe vacuo-to-water solvatochromic shifts," The Journal of Chemical Physics 157, 214101 (2022).

¹²L. V. Slipchenko, "Solvation of the Excited States of Chromophores in Polarizable Environment: Orbital Relaxation versus Polarization †," The Journal of Physical Chemistry A 114, 8824-8830 (2010).

¹³F. Lipparini, C. Cappelli, and V. Barone, "Linear Response Theory and Electronic Transition Energies for a Fully Polarizable QM/Classical Hamiltonian," Journal of Chemical Theory and Computation 8, 4153-4165 (2012).

¹⁴X. Liu, A. Humeniuk, and W. J. Glover, "Conical Intersections in Solution with Polarizable Embedding: Integral-Exact Direct Reaction Field," Journal of Chemical Theory and Computation 18, 6826-6839 (2022).

¹⁵A. Humeniuk and W. J. Glover, "Multistate, Polarizable QM/MM Embedding Scheme Based on the Direct Reaction Field Method: Solvatochromic Shifts, Analytical Gradients and Optimizations of Conical Intersections in Solution," Journal of Chemical Theory and Computation 20, 2111-2126 (2024).

¹⁶R. Improta, V. Barone, G. Scalmani, and M. J. Frisch, "A state-specific polarizable continuum model time dependent density functional theory method for excited state calculations in solution," The Journal of Chemical Physics 125, 054103 (2006).

¹⁷D. Loco, E. Polack, S. Caprasecca, L. Lagardere, F. Lipparini, J.-P. Piquemal, and B. Mennucci, "A QM/MM Approach Using the AMOEBA Polarizable Embedding: From Ground State Energies to Electronic Excitations," Journal of Chemical Theory and Computation 12, 3654-3661 (2016).

¹⁸B. T. Thole and P. T. Duijnen, "On the quantum mechanical treatment of solvent effects," Theoretica Chimica Acta 55, 307-318 (1980).

¹⁹B. Thole and P. Van Duijnen, "The direct reaction field hamiltonian: Analysis of the dispersion term and application to the water dimer," Chemical Physics 71, 211-220 (1982).

²⁰A. Humeniuk and W. J. Glover, "Efficient CPU and GPU implementations of multicenter integrals over long-range operators using Cartesian Gaussian functions," Computer Physics Communications 280, 108467 (2022).

²¹T. P. Fay, N. Ferré, and M. Huix-Rotllant, "Efficient polarizable qm/mm using the direct reaction field hamiltonian with electrostatic potential fitted multipole operators," Journal of Chemical Theory and Computation 21, 183-201 (2025).

²²N. Ferré and J. G. Ángyán, "Approximate electrostatic interaction operator for QM/MM calculations," Chemical Physics Letters 356, 331-339 (2002).

²³J. P. Pederson and J. G. McDaniel, "DFT-based QM/MM with particle-mesh Ewald for direct, long-range electrostatic embedding," The Journal of Chemical Physics 156, 174105 (2022).

²⁴S. Kvedaravičiūtė, D. Carrasco-Busturia, K. B. Møller, and J. M. H. Olsen, "Polarizable Embedding without Artificial Boundary Polarization," Journal of Chemical Theory and Computation 19, 5122-5141 (2023).

²⁵T. Darden, D. York, and L. Pedersen, "Particle mesh Ewald: An N log(N) method for Ewald sums in large systems," The Journal of Chemical Physics 98, 10089-10092 (1993).

²⁶S. Bonfrate, N. Ferré, and M. Huix-Rotllant, "An efficient electrostatic embedding QM/MM method using periodic boundary conditions based on particle-mesh Ewald sums and electrostatic potential fitted charge operators," The Journal of Chemical Physics 158, 021101 (2023).

²⁷Z. C. Holden, B. Rana, and J. M. Herbert, "Analytic gradient for the QM/MM-Ewald method using charges derived from the electrostatic potential: Theory, implementation, and application to ab initio molecular dynamics simulation of the aqueous electron," The Journal of Chemical Physics 150, 144115 (2019).

²⁸C. Li and G. K.-L. Chan, "Accurate QM/MM Molecular Dynamics for Periodic Systems in GPU4P <span style="font-variant:small-caps;">> y</span>/span> SCF with Applications to Enzyme Catalysis," Journal of Chemical Theory and Computation 21, 803-816 (2025).

²⁹www.github.com/tomfay/openespf (2025), OpenESPF python package implementing DREEM, date accessed: 27-05-2025.

³⁰W. Smith, "Point Multipoles in the Ewald Summation," CCP5 Newsletter (1982).

³¹J. A. Pople, R. Krishnan, H. B. Schlegel, and J. S. Binkley, "Derivative studies in hartree-fock and møller-plesset theories," International Journal of Quantum Chemistry 16, 225-241 (2009).

³²F. Furche and R. Ahlrichs, "Adiabatic time-dependent density functional methods for excited state properties," The Journal of Chemical Physics 117, 7433-7447 (2002).

³³Q. Zeng, J. Liu, and W. Liang, "Molecular properties of excited electronic state: Formalism, implementation, and applications of analytical second energy derivatives within the framework of the time-dependent density functional theory/molecular mechanics," The Journal of Chemical Physics 140, 18A506 (2014).

³⁴Q. Sun, "Libcint: An efficient general integral library for Gaussian basis functions," Journal of Computational Chemistry 36, 1664-1671 (2015).

³⁵Q. Sun, T. C. Berkelbach, N. S. Blunt, G. H. Booth, S. Guo, Z. Li, J. Liu, J. D. McClain, E. R. Sayfutyarova, S. Sharma, S. Wouters, and G. K. Chan, "PySCF: the Python-based simulations of chemistry framework," WIREs Computational Molecular Science 8, e1340 (2018).

³⁶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, "Recent developments in the PySCF program package," The Journal of Chemical Physics 153, 024109 (2020).

³⁷P. Eastman, J. Swails, J. D. Chodera, R. T. McGibbon, Y. Zhao, K. A. Beauchamp, L.-P. Wang, A. C. Simmonett, M. P. Harrigan, C. D. Stern, R. P. Wiewiora, B. R. Brooks, and V. S. Pande, "OpenMM 7: Rapid development of high performance algorithms for molecular dynamics," PLOS Computational Biology 13, e1005659 (2017).

³⁸P. Eastman, R. Galvelis, R. P. Peláez, C. R. A. Abreu, S. E. Farr, E. Gallicchio, A. Gorenko, M. M. Henry, F. Hu, J. Huang, A. Krämer, J. Michel, J. A. Mitchell, V. S. Pande, J. P. Rodrigues, J. Rodriguez-Guerra, A. C. Simmonett, S. Singh, J. Swails, P. Turner, Y. Wang, I. Zhang, J. D. Chodera, G. De Fabritiis, and T. E. Markland, "OpenMM 8: Molecular Dynamics Simulation with Machine Learning Potentials," The Journal of Physical Chemistry B 128, 109-116 (2024).

³⁹C. Amovilli and R. McWeeny, "A matrix partitioning approach to the calculation of intermolecular potentials. General theory and some examples," Chemical Physics 140, 343-361 (1990).

⁴⁰C. Amovilli and B. Mennucci, "Self-Consistent-Field Calculation of Pauli Repulsion and Dispersion Contributions to the Solvation Free Energy in the Polarizable Continuum Model," The Journal of Physical Chemistry B 101, 1051-1057 (1997).

⁴¹T. Giovannini, P. Lafiosca, and C. Cappelli, "A General Route to Include Pauli Repulsion and Quantum Dispersion Effects in QM/MM Approaches," Journal of Chemical Theory and Computation 13, 4854-4870 (2017).

⁴²B. Thole, "Molecular polarizabilities calculated with a modified dipole interaction," Chemical Physics 59, 341-350 (1981).

⁴³J.-P. Piquemal, N. Gresh, and C. Giessner-Prettre, "Improved Formulas for the Calculation of the Electrostatic Contribution to the Intermolecular Interaction Energy from Multipolar Expansion of the Electronic Distribution," The Journal of Physical Chemistry A 107, 10353-10359 (2003), publisher: American Chemical Society (ACS).

⁴⁴T. P. Fay and D. T. Limmer, "Unraveling the mechanisms of triplet state formation in a heavy-atom free photosensitizer," Chemical Science 15, 6726-6737 (2024).

⁴⁵W. L. Jorgensen, D. S. Maxwell, and J. Tirado-Rives, "Development and Testing of the OPLS All-Atom Force Field on Conformational Energetics and Properties of Organic Liquids," Journal of the American Chemical Society 118, 11225-11236 (1996).

⁴⁶L. Jensen, P. T. Van Duijnen, and J. G. Snijders, "A discrete solvent reaction field model within density functional theory," The Journal of Chemical Physics 116, 084101 (2002).
Physics **118**, 514–521 (2003).

<footnote n="47"><sup>47</sup>S. Hirata and M. Head-Gordon, “Time-dependent density functional theory within the Tamm–Dancoff approximation,” Chemical Physics Letters **314**, 291–299 (1999).</footnote>

<footnote n="48"><sup>48</sup>Y.-S. Lin, G.-D. Li, S.-P. Mao, and J.-D. Chai, “Long-Range Corrected Hybrid Density Functionals with Improved Dispersion Corrections,” Journal of Chemical Theory and Computation **9**, 263–272 (2013).</footnote>

<footnote n="49"><sup>49</sup>M. Hijazi, D. M. Wilkins, and M. Ceriotti, “Fast-forward Langevin dynamics with momentum flips,” The Journal of Chemical Physics **148**, 184109 (2018).</footnote>

<footnote n="50"><sup>50</sup>Z. R. Wiethorn, K. E. Hunter, T. J. Zuehlsdorff, and A. Montoya-Castillo, “Beyond the Condon limit: Condensed phase optical spectra from atomistic simulations,” The Journal of Chemical Physics **159**, 244114 (2023).</footnote>

<footnote n="51"><sup>51</sup>A. K. Dutta, F. Neese, and R. Izsák, “Towards a pair natural orbital coupled cluster method for excited states,” The Journal of Chemical Physics **145**,</footnote>

<footnote n="034102">034102 (2016).</footnote>

<footnote n="52"><sup>52</sup>A. K. Dutta, M. Nooijen, F. Neese, and R. Izsák, “Exploring the Accuracy of a Low Scaling Similarity Transformed Equation of Motion Method for Vertical Excitation Energies,” Journal of Chemical Theory and Computation **14**, 72–91 (2018).</footnote>

<footnote n="53"><sup>53</sup>F. Neese, “Software Update: The <span style="font-variant:small-caps;">&gt;ORCA</span><sup>&lt;/span&gt;</sup> Program System—Version 6.0,” WIREs Computational Molecular Science **15**, e70019 (2025).</footnote>

<footnote n="54"><sup>54</sup>K. Aidas, K. V. Mikkelsen, B. Mennucci, and J. Kongsted, “Fluorescence and phosphorescence of acetone in neat liquid and aqueous solution studied by QM/MM and PCM approaches,” International Journal of Quantum Chemistry **111**, 1511–1520 (2011), publisher: Wiley.</footnote>

<footnote n="55"><sup>55</sup>Y. Orozco-González, K. Coutinho, and S. Canuto, “Excited state electronic polarization and reappraisal of the n ← π* emission of acetone in water,” Chemical Physics Letters **499**, 108–112 (2010).</footnote>