# Elastic Collision Based Dynamic Partitioning Scheme for Hybrid Simulations

Björn Kirchhoff  [  Elvar Örn Jónsson  [  Asmus Ougaard Dohn  [  Timo Jacob  [
Hannes Jónsson<EMAIL_ADDRESS>[

###### Abstract

The scattering-adapted flexible inner region ensemble separator (SAFIRES) is a
partitioning scheme designed to divide a simulation cell into two regions to
be treated with different computational methodologies. SAFIRES prevents
particles from crossing between regions and resolves boundary events through
elastic collisions of the particles mediated by the boundary, conserving
energy and momenta. A multiple-time-step propagation algorithm is introduced
where the time step is scaled automatically to identify the moment a collision
occurs. If the length of the time step is kept constant, the new propagator
reduces to a regular algorithm for Langevin dynamics, and to the velocity
Verlet algorithm for classical dynamics if the friction coefficient is set to
zero. SAFIRES constitutes the exact limit of the premise behind boundary-based
methods such as FIRES, BEST, and BCC which take advantage of the
indistinguishability of molecules on opposite sides of the separator. It gives
correct average ensemble statistics despite the introduction of an ensemble
separator. SAFIRES is tested in simulations where the molecules on the two
sides are treated in the same way, for a Lennard-Jones (LJ) liquid and a LJ
liquid in contact with a surface, as well as for liquid modelling simulations
using the TIP4P force field. Simulations using SAFIRES are shown to reproduce
the unconstrained reference simulations without significant deviations.

###### keywords:

hybrid simulations, partitioning scheme, theoretical electrochemistry, atomic
simulation environment

HI] Science Institute and Faculty of Physical Sciences, University of Iceland,
VR-III, 107 Reykjavík (Iceland). HI] Science Institute and Faculty of Physical
Sciences, University of Iceland, VR-III, 107 Reykjavík (Iceland). HI] Science
Institute and Faculty of Physical Sciences, University of Iceland, VR-III, 107
Reykjavík (Iceland). [DTU] Technical University of Denmark, Lyngby, Denmark.
UUlm] Institute of Electrochemistry, Ulm University, Albert-Einstein-Allee 47,
89081 Ulm (Germany). [HIU] Helmholtz-Institute Ulm (HIU) Electrochemical
Energy Storage, Helmholtz-Straße 16, 89081 Ulm (Germany). [KIT] Karlsruhe
Institute of Technology (KIT), P.O. Box 3640, 76021 Karlsruhe (Germany). HI]
Science Institute and Faculty of Physical Sciences, University of Iceland, VR-
III, 107 Reykjavík (Iceland).

##  1 Introduction

With the advent of and widespread access to high-performance computing
resources over the last three decades, computational methods have become
increasingly important in predicting material properties and understanding
chemical reaction mechanisms. However, the desire to correctly describe
processes at interfaces has been straining the limits of standard
computational methods based on density functional theory (DFT). The solid-
liquid interface is of particular interest and has proven especially
troublesome since solvents — and especially water — often take part in
reactions rather than acting as neutral bystanders1. The influence of
solvation can alter reaction behavior significantly2, 3, 4, 5 and can even
cause structural changes in the catalyst through solvent-induced surface
rearrangement6, 7. Correct description of solvent interactions is oftentimes
crucial when trying to understand complex natural phenomena such as the water
splitting reaction in photosystem II8 or substrate binding in protein active
sites9.

Unfortunately, only few explicit solvent molecules can be included in DFT
calculations due to quickly inflating computational effort. Furthermore,
static DFT calculations, even when human bias is limited via the use of global
optimization methods10, 11, 12, may only serve as an approximation for the
first rigid ice-like layers of H2O on transition metal surfaces or for ordered
water clusters. However, this approach is less useful for less strongly
interacting molecules and for the description of bulk-like properties of the
solvent where less ordered structures are expected.

When studying the electronic interactions of, for example, a solid-liquid
interface system, the interactions that are most important for the correct
description of the reactive properties of a model are those between the
surface, the adsorbate, and the first few layers of solvent. Relevant
interactions include charge transfer, hybridization of electronic states, and
proton transfer reactions. The widely used implicit solvation models, which
represent the solvent as a potential based on key bulk properties such as the
dielectric constant and dipole moment, in combination with DFT calculations
fail to capture these important effects and are often reported to not improve
results over calculations without implicit solvation13, 14, 15, 16.

One possible solution to this problem is to divide the simulation into two
coupled regions which are treated using different computational methodologies.
A prominent example of this approach are QM/MM17, 18, 19 methods. Here, the
surface, the active species, and few solvent molecules are calculated using
accurate and expensive QM methods while the majority of the (entirely
explicit) solvent molecules is computed using efficient force field methods.
This significant computational speedup allows users to run molecular dynamics
(MD) or Monte Carlo (MC) methods to sample the phase space and to extract the
solvation energy as a thermodynamic average.

A crucial facet in any coupling scheme is a sound description of the boundary
between inner and outer region. Several different approaches exist to
implement this necessary but inherently nonphysical concept while causing the
least amount of disturbance to the system. In adaptive methods, particles are
allowed to exchange between inner and outer regions20, 21, 22, 23, 24, 25, 26,
27, 28. When different potentials are used in the description of the inner and
outer region, this approach can lead to discontinuities where the potentials
come into contact28. In order to retain smooth dynamics and correct ensemble
properties when particles cross between regions, a large number of different
configurations of particles in a buffer region around the boundary have to be
calculated if consistency of the total energy and forces is to be achieved.
This corrective process increases computational costs by orders of magnitude.

Another boundary approach is presented by Rowley and Roux in form of the
Flexible Inner Region Ensemble Separator (FIRES)29. With FIRES, particles
cannot transfer between regions but the boundary can expand and contract based
on the position of the outermost particle in the inner region. This implies
that a typical FIRES simulation contains three distinct parts: a solute that
acts as the central point of the simulation, an inner region of particles
around the solute, and an outer region of particles; the solute and the inner
region are treated using the same computational methodology. The boundary
radius is defined as the distance between the solute and the outermost
particle in the inner region and enforced using a Hookean force constraint
acting on the inner and outer particles.

While this method sacrifices microscopic trajectories of individual particles
in the system, it can be mathematically shown to deliver correct average
thermodynamic values30, 29 if the same chemical species is present in both
regions. FIRES therefore establishes the main premise that all boundary-based
separation algorithms rely on: it is assumed that a situation where a particle
from the outer region travels to the inside et vice verse is statistically
equivalent to a situation where the particles are redirected back where they
came from at the point of equal distance to the solute. Thus, correct average
thermodynamic properties can be retained this way at the cost of sacrificing
continuity of individual particle trajectories. FIRES was used to good effect
in benchmark calculations for DFTB3 paramterization of Zn and Mg31 and to
calculate the solvation of a transition state of a nucleophilic carbonyl
reaction32.

Notably, the finite force constraint used in the Rowley and Roux
implementation is actually an approximation to the FIRES premise, which in its
mathematical formalism implies instantaneous redirection at the boundary. It
is conceivable that this approximation could lead to what Bulo et al. reported
as an accumulation of molecules at the boundary when using FIRES, analogous to
what is observed for coupled simulations with a static, hard sphere
boundary33.

Aside from FIRES, other boundary-based methods have been developed, in
particular the Boundary based on Exchange Symmetry (BEST) method 34, 35 and
the Boundary Constraint with Correction (BCC) method36. BEST aims to find a
generalized approach to the separating potential. The BEST separating
potential is constructed by assigning particles to the inner or outer region
and calculating a penalty function which is unity for a separated particle
pair and approaches zero as the particle configuration becomes more
undesirable, i.e. as a particle pair travels further into the respective other
region34, 35. The bias function is then obtained as the product of penalty
functions of all permutations of particle pairs; hence, the reference to
exchange symmetry. However, since particle pairs that are unproblematic (i.e.
in different phases) do not contribute significantly to the bias function,
Shiga and Masia find that construction of the bias function can be simplified
by only taking into account single or double exchanges of the particle pairs
with the penalty functions closest to zero.

The BCC method utilizes Fermi-derived bias potential inspired by BEST to
separate particles in the regions36. As the authors remark, the basic premise
behind boundary-based methods — instantaneous redirection at the boundary —
cannot be achieved using a constraining potential. The BCC method therefore
attempts to rectify this statistical error, which is related to the
accumulation artifacts at the boundary discussed for FIRES, by performing an
additional set of calculations to remove contributions of the bias potential
from thermodynamic properties of the erroneous configurations. This approach
therefore adds additional computational overhead to simulations.

In this study, a new separation algorithm, the scattering-adapted flexible
inner region ensemble separator (SAFIRES), is introduced. SAFIRES divides a
computational model into two regions to treat them using different
methodologies. Analogous to other boundary-based methods, the boundary between
the inner and outer region is defined flexibly between the solute (which can
be molecular or periodic) and the outermost particle in the inner region. To
enforce the boundary, SAFIRES performs elastic collisions between particles
from the inner and outer region when they are at the same distance from the
solute. SAFIRES therefore fully conserves energy and momenta. This
instantaneous approach to particle redirection however necessitates that time
steps be flexible during the simulation for collisions to be exact. Therefore,
fractional time steps are calculated on-the-fly as conflicts arise and a new
multiple-time-step propagation algorithm is introduced. The new propagator
reduces to the Vanden-Eijnden / Ciccotti implementation of the Langevin
algorithm37 for constant time steps and further reduces to the Velocity-Verlet
algorithm for constant time steps and a friction coefficient of zero. As a
final requirement, forces are only updated at full intervals of the initial
time step in order to retain the temporal consistency of the calculation and
to satisfy time-reversibility of the Taylor expressions that the propagation
algorithm is derived from. Figure 1 illustrates the SAFIRES approach for
molecular and surface model systems in a schematic way.

Figure 1: Illustration of the SAFIRES approach for a molecular (a) and an
infinite-surface based model system (b) with water. The flexible boundary is
shown as a dotted black circle or line. Molecules in the inner region are
shown in red, those in the outer region are shown in orange. The surface is
shown in grey. $r$: vector connecting the solute (a) or a plane parallel to
the surface (b) to the particles in the inner (in) and outer (out) region. The
greyed out and dotted $\vec{v}$ vectors represent the initial velocity of the
molecule in the inner (in) and outer (out) region before the SAFIRES algorithm
redirects them. The vectors $\vec{v}^{\prime}$ represent the velocity vectors
of the molecule in the inner (in, green) and outer (out, blue) region after
the SAFIRES algorithm redirected them via an elastic collision mediated by the
boundary.

This study is structured as follows: section 2 details the implementation of
SAFIRES. Therein, subsection 2.1 gives an overview of the SAFIRES algorithm,
subsection 2.2 details propagation and how fractional time steps are
calculated to find the exact point of collision, subsection 2.3 examines the
collision handling, and subsection 2.4 delves into advanced details of the
implementation. Subsequently, sections 3 and 4 describe the computational
details and test calculations of SAFIRES, respectively. The SAFIRES algorithm
is first illustrated based on a Lennard-Jones (LJ) liquid model system with a
molecular solute (subsection 4.1) before expanding the approach to a LJ
infinite-surface model as the solute for applications at the solid-liquid
interface (subsection 4.2). Inert noble gas surfaces have been used in the
past, for example in a QM/MM study by Daru et al. in order to study the
rearrangement of ion-doped water clusters on an inert surface38. Partial RDFs
can be used to study the structure of a solvent near a surface39. Finally,
SAFIRES is tested in an MM*/MM simulation of a water-in-water model using the
TIP4P force field (subsection 4.3). Radial distribution functions (RDFs) are
used in all cases to evaluate the ability of SAFIRES to reproduce the correct
average structural features of unconstrained simulations.

##  2 The SAFIRES Method

###  2.1 Algorithm overview

The key steps in SAFIRES for detecting and resolving a crossing of particles
from the outer to the inner region are illustrated with a spherically-
symmetric model system, as depicted schematically in Figure 2.

Figure 2: Schematic overview of the steps involved in the SAFIRES scheme. a
The initial, conflict-free state of the system. $\vec{r}_{\text{b}}$ indicates
the boundary radius, which is defined by the distance between the solute and
the outermost inner region (in) particle. $\vec{v}_{\text{in}}$ and
$\vec{v}_{\text{out}}$ indicate the velocity vectors of an inner and outer
region particle, respectively. b After the system is propagated using the
default time step $dt$, the SAFIRES algorithm checks if one or more outer
particles are closer to the solute than the outermost inner particle. If a
conflict is detected, the SAFIRES algorithm proceeds to resolve the issue. c
The first step to resolve the issue is to reset the configuration to the last
conflict-free state. d A fractional time step $dt^{\prime}$ is solved for and
executed which exactly propagates the conflicting in/out particle pair to the
boundary. $\vec{v}_{\perp}$ and $\vec{v}_{\parallel}$ indicate the normal and
tangential velocity components, respectively. e An elastic collision between
the conflicting in/out pair is performed, exchanging propagating components
normal to the boundary surface. f After the collision has changed the
trajectories of the conflicting in/out pair, the system is propagated by the
remaining time step, $dt^{\prime}-dt$, required to fulfill a complete regular
time step. Afterwards, SAFIRES checks if another outer particle has propagated
beyond the outermost inner particle. If this is the case, SAFIRES repeats
steps c-f until there no more conflicts detected and a full time step $dt$ is
achieved.

For simplicity it is assumed that all particles are monoatomic and that the
solute particle is fixed; particles with more than one atom are discussed in
section 2.4. The SAFIRES routine for a 2D-periodic system is largely the same
with some simplifications as outlined in Section 2.4.2.

At the start of an MD simulation, three sets of indices are assigned. The
first, $s$, identifies the solute particle, whose center of mass acts as the
origin for the boundary radius. The second set, $S_{\text{in}}$, identifies
the particles within this radius, and lastly, $S_{\text{out}}$ identifies the
particles outside of this boundary.

Given the 3$n$-dimensional position vector field of all $n$ particles of the
system, $\mathbf{R}^{i}$, at iteration $i=0$, the index assignment has to
satisfy

$\\{\mid\mid\mathbf{R}^{i}(\alpha)-\mathbf{R}^{i}(s)\mid\mid\ \forall\alpha\in
S_{\text{in}}\\}<\ \\{\mid\mid\mathbf{R}^{i}(\beta)-\mathbf{R}^{i}(s)\mid\mid\
\forall\beta\in S_{\text{out}}\\}\,,$ (1)

i.e. all of the outer particles should be further away from the solute
compared to all of the inner particles. $||\mathbf{R}||$ denotes the Euclidian
norm. Note that the condition simply involves a change of origin of the input
coordinate space, and any fixed point of origin can in principle be used – it
does not have to be associated with a particle.

The SAFIRES algorithm is called after an integrated time step $dt>0$ of the
superordinate MD simulation is performed which propagates the particle
positions, leading to

$\mathbf{R}^{i}\xrightarrow{dt}\mathbf{R}^{i+1}\,.$ (2)

First, the radius of the boundary sphere, $r_{b}$, which is shown as a dotted
circle in Figure 2, is updated. The radius is based on the farthest distance
between the fixed solute (s) and inner (in) particles

$\alpha_{b}\leftarrow
r_{b}=\text{max}\left\\{\mid\mid\mathbf{R}^{i+1}(\alpha)-\mathbf{R}^{i+1}(s)\mid\mid\forall\alpha\in
S_{\text{in}}\right\\}$ (3)

Here $\alpha_{b}$ indexes the outermost inner particle. SAFIRES then checks if
any of the outer (out) particles have crossed the boundary

$S^{\prime}_{\text{out}}=\left\\{\mid\mid\mathbf{R}^{i+1}(\beta)-\mathbf{R}^{i+1}(s)\mid\mid<r_{b}\
\forall\beta\in S_{\text{out}}\right\\}$ (4)

where $S^{\prime}_{\text{out}}$ is a subset of $S_{\text{out}}$ and identifies
the conflicting outer particles (it can have more than one element). If no
conflicts are detected (i.e. $S^{\prime}_{\text{out}}=\O$), the systems state
is stored and the next iteration of the MD integrator performed.

If a conflict is detected (Figure 2b), the simulation is reverted to the last
conflict-free system state (Figure 2c). From this state, a fractional time
step $dt^{\prime}\leq dt$ is determined – for all conflicting particle pairs
as indexed by $\alpha_{b}$ and $S^{\prime}_{\text{out}}$ – and the smallest
$dt^{\prime}$ executed in order to propagate the conflicting in/out particle
pair to the same distance from the solute

$\mathbf{R}^{i}\xrightarrow{dt^{\prime}}\mathbf{R}^{i+\delta}\,,$ (5)

where $\delta$ indicates that only a partial iteration is performed (which can
be defined as $\delta=dt^{\prime}/dt$ for completeness). The smallest value of
$dt^{\prime}$ is resolved first since it corresponds to the first boundary
crossing in the time step interval $dt$. The time step extrapolation that
yields $dt^{\prime}$ is discussed in section 2.2.

The partial propagation, leading to the state shown in Figure 2d, is handled
by SAFIRES. In this state, an elastic collision is performed between the
conflicting in/out particle pair resulting in an exchange of propagating
vectors normal to the boundary surface, $\mathbf{v}_{\perp}^{b}$, resulting in
Figure 2e. This step is detailed in section 2.3. After the collision has been
performed successfully, the resulting system state is saved in case the
ensemble needs to be reverted again.

Lastly, the system is propagated by the remaining fractional time step $dt-
dt^{\prime}$ that is required to complete a regular time step $dt$, see Figure
2f. The propagation is performed by the SAFIRES propagator and uses the post-
collision velocity vectors

$\mathbf{R}^{i+\delta}\xrightarrow{(dt-
dt^{\prime}),\mathbf{v}_{\perp}^{b}}\mathbf{R}^{i+1}_{\text{c}}$ (6)

After this second fractional propagation step, the boundary is updated based
on the new positions and another check for conflicts is performed using
$\mathbf{R}^{i+1}_{c}$ in equations (3)-(4). This updates both the index
$\alpha_{b}$ and index set $S^{\prime}_{\text{out}}$ since having resolved the
first boundary crossing any other inner and outer particle pair can become a
new conflicting pair. If a new conflict arises the steps depicted in Figure
2c-f are repeated. Multiple conflicts are detailed in section 2.4.

After a conflict free system state is reached, the SAFIRES propagator
completes the iteration. Control is then returned to the superordinate MD
propagator and a new iteration is started.

###  2.2 Time step extrapolation

The relationship for extrapolation of the time step required to propagate the
conflicting in/out particle pair to the same distance from the solute is
obtained from the Langevin propagator37

$\displaystyle\mathbf{V}^{i+\frac{1}{2}}$ $\displaystyle=$
$\displaystyle\mathbf{V}^{i}+\frac{1}{2}dt\mathbf{F}(\mathbf{R}^{i})-\frac{1}{2}dt\gamma\mathbf{V}^{i}+\frac{1}{2}\sqrt{dt}\mathbf{\Sigma}^{i}\bm{\sigma}$
(7) $\displaystyle-$
$\displaystyle\frac{1}{8}dt^{2}\gamma\left(\mathbf{F}(\mathbf{R}^{i})-\gamma\mathbf{V}^{i}\right)-\frac{1}{4}dt^{\frac{3}{2}}\gamma\left(\frac{1}{2}\mathbf{\Sigma}^{i}+\frac{1}{\sqrt{3}}\mathbf{\Lambda}^{i}\right)\bm{\sigma}$
$\displaystyle\mathbf{R}^{i+1}$ $\displaystyle=$
$\displaystyle\mathbf{R}^{i}+dt\mathbf{V}^{i+\frac{1}{2}}+dt^{\frac{3}{2}}\frac{1}{2\sqrt{3}}\mathbf{\Lambda}^{i}\bm{\sigma}$
(8) $\displaystyle\mathbf{V}^{i+1}$ $\displaystyle=$
$\displaystyle\mathbf{V}^{i+\frac{1}{2}}+\frac{1}{2}dt\mathbf{F}(\mathbf{R}^{i+1})-\frac{1}{2}dt\gamma\mathbf{V}^{i+\frac{1}{2}}+\frac{1}{2}\sqrt{dt}\mathbf{\Sigma}^{i}\bm{\sigma}$
(9) $\displaystyle-$
$\displaystyle\frac{1}{8}dt^{2}\gamma\left(\mathbf{F}(\mathbf{R}^{i+1})-\gamma\mathbf{V}^{i+\frac{1}{2}}\right)-\frac{1}{4}dt^{\frac{3}{2}}\gamma\left(\frac{1}{2}\mathbf{\Sigma}^{i}+\frac{1}{\sqrt{3}}\mathbf{\Lambda}^{i}\right)\bm{\sigma}$
$\displaystyle\bm{\sigma}$ $\displaystyle=$
$\displaystyle\sqrt{\frac{2T\gamma}{\mathbf{M}}}$ (10)

Here, $\mathbf{V}$ and $\mathbf{F}$ are the 3$n$-dimensional velocity and
force vector fields, respectively, for all $n$ particles. Similarly and for
the following sections, lowercase boldfaced letters – for example $\mathbf{v}$
and $\mathbf{r}$ – refer to their 3-dimensional single particle analogs.
$\gamma$ is the thermostat friction coefficient, and $\mathbf{\Sigma}$ and
$\mathbf{\Lambda}$ are 3$n$-dimensional covariant Gaussian distributions that
are randomized at the start of each iteration $i$. $T$ is the target
temperature for the thermostat and $\mathbf{M}$ is a $n$-dimensional column
vector of the particle masses. Setting $\gamma=0$ reduces this set of
equations to the Velocity-Verlet propagator.

The propagation during the SAFIRES routine is split into two parts — first to
propagate towards the boundary, and then to propagate away from it again.
Forces are only evaluated again after a full time step has been completed and
the coordinate space is without additional conflicts. Therefore, equation (9)
is performed after all conflicts are resolved. As a result, only equations
(7)-(8) need to be considered for the time step extrapolation of conflicting
in/out particle pairs.

In order to satisfy the time reversibility of the Taylor expressions that the
Velocity Verlet and Langevin propagators are derived from, all components in
equation (7) and the third term in equation (8) are added based on the full
time step, $dt$. Equation (8) is then reduced to

$\displaystyle\mathbf{R}^{i+1}$ $\displaystyle=$
$\displaystyle\mathbf{R}^{i}+dt\mathbf{A}^{i}$ (11)
$\displaystyle\mathbf{A}^{i}$ $\displaystyle=$
$\displaystyle\mathbf{V}^{i+\frac{1}{2}}+\mathbf{B}^{i}$ (12)
$\displaystyle\mathbf{B}^{i}$ $\displaystyle=$ $\displaystyle
dt^{\frac{1}{2}}\frac{1}{2\sqrt{3}}\mathbf{\Lambda}^{i}\bm{\sigma}$ (13)

and the propagation of the coordinate space in equation (11) is now linear in
terms of $dt$. Given $\mathbf{R}^{i+1}$ as defined above, the outermost inner
particle is identified via equation (3), and the condition of equation (4) is
checked. For $\alpha_{b}$ and elements $\beta$ in $S^{\prime}_{\text{out}}$
all conflicting in/out particle pairs are considered. For the time step
extrapolation, the key quantities are, however, not the particle positions
$\mathbf{R}^{i}(\alpha_{b})$ and $\mathbf{R}^{i}(S^{\prime}_{\text{out}})$ but
the distances of the in/out particle pair from the solute, which are more
naturally traced with a change in origin
$\mathbf{R}^{i}_{s}=\mathbf{R}^{i}-\mathbf{R}^{i}(s)$. Hence, in this frame
the norm of the positions of the conflicting in/out pairs must be equal at the
point of collision:

$\displaystyle||\mathbf{R}^{i+\delta}_{s}(\alpha_{b})||$
$\displaystyle\overset{!}{=}$
$\displaystyle||\mathbf{R}^{i+\delta}_{s}(\beta)||\ \ \forall\beta\in
S^{\prime}_{\text{out}}.$ (14)

Using equation (14) in combination with equation (11) leads to

$||\mathbf{R}_{s}^{i}(\alpha_{b})+\delta
t^{\prime}(\beta)\mathbf{A}^{i}(\alpha_{b})||=||\mathbf{R}_{s}^{i}(\beta)+\delta
t^{\prime}(\beta)\mathbf{A}^{i}(\beta)||\ \ \forall\beta\in
S^{\prime}_{\text{out}}$ (15)

and is solved analytically for $\delta t^{\prime}(\beta)$, which are the
partial time steps required to propagate each conflicting in/out particle
pairs to the same distance from the solute. Using

$\beta_{b}\leftarrow dt^{\prime}=\text{min}\\{\delta
t^{\prime}(\beta)\mid\forall\delta t^{\prime}(\beta)\in\mathbb{R}_{+}\\}$ (16)

the temporally first boundary conflict and corresponding outer particle index,
$\beta_{b}$, is identified, i.e. the smallest positive real valued solution.
The system’s coordinate space is evolved accordingly

$\mathbf{R}^{i+\delta}=\mathbf{R}^{i}+dt^{\prime}\mathbf{A}^{i}$ (17)

leading to the state depicted in Figure 2d. Note that a positive real valued
solution where $dt^{\prime}\leq dt$ is always guaranteed since $dt>0$ and due
to the condition set by equations (1) and (4).

###  2.3 Collision handling

Once a conflict has been detected and the system propagated such that the
inner ($\alpha_{b}$) and outer ($\beta_{b}$) particle are both at the
boundary, the particles need to be redirected in order to resolve the
conflict. This is achieved with an exchange of momentum mediated by adjusting
the propagating vectors normal to the boundary surface. In the case of
Langevin dynamics the vectors $\mathbf{a}^{i}_{\text{in}}$ =
$\mathbf{A}^{i}(\alpha_{b})$ and
$\mathbf{a}^{i}_{\text{out}}=\mathbf{A}^{i}(\beta_{b})$ are considered. These
components reduce to $\mathbf{v}^{i+\frac{1}{2}}_{\text{in}}$ and
$\mathbf{v}^{i+\frac{1}{2}}_{\text{out}}$ in Velocity-Verlet. In the following
the MD step index $i$ is omitted for clarity.

In case of a solute with a spherical boundary, a transformation of the
coordinates of one of the involved particles is required to bring them into
the same frame of reference. Geometrically speaking, the propagating vector of
one particle needs to be rotated on the surface of the boundary sphere to the
position of the other particle. To this end, consider the two vectors
$\mathbf{r}_{\text{in,s}}=\mathbf{R}(\alpha_{b})-\mathbf{R}(s)$ and
$\mathbf{r}_{\text{out,s}}=\mathbf{R}(\beta_{b})-\mathbf{R}(s)$ connecting the
solute and the in/out particle pair. The angle $\phi$ between these vectors is
given by Vincenty’s formula40 as

$\phi=\text{arctan2}\left(||\mathbf{r}_{\text{in,s}}\times\mathbf{r}_{\text{out,s}}||,(\mathbf{r}_{\text{in,s}}\cdot\mathbf{r}_{\text{out,s}})\right).$
(18)

Assuming rotation of the outer particle the rotated propagating vector is
given by

$\mathbf{a}^{\text{rot}}_{\text{out}}=X(\mathbf{c},\phi)\cdot\mathbf{a}_{\text{out}},$
(19)

with a rotational matrix, $X(\mathbf{c},\phi)$, of the Euler-Rodrigues form
and the rotational axis

$\mathbf{c}=\frac{\mathbf{r}_{\text{in,s}}\times\mathbf{r}_{\text{out,s}}}{||\mathbf{r}_{\text{in,s}}\times\mathbf{r}_{\text{out,s}}||}.$
(20)

After $\mathbf{a}_{\text{out}}$ is transformed into the same reference frame
as $\mathbf{a}_{\text{in}}$, an elastic collision is performed between the two
particles according to

$\displaystyle\mathbf{a}_{\perp}^{b}$ $\displaystyle=$
$\displaystyle\frac{((\mathbf{a}^{\text{rot}}_{\text{out}}-\mathbf{a}_{\text{in}})\cdot\mathbf{r}_{\text{in,s}})\
\mathbf{r}_{\text{in,s}}}{||\mathbf{r}_{\text{in,s}}||^{2}}$ (21)
$\displaystyle\mathbf{a}_{\text{out}}^{\prime\text{\ rot}}$ $\displaystyle=$
$\displaystyle\mathbf{a}_{\text{out}}^{\text{rot}}-\frac{2\
m_{\text{in}}}{m_{\text{in}}+m_{\text{out}}}\mathbf{a}_{\perp}^{b}$ (22)
$\displaystyle\mathbf{a}_{\text{in}}^{\prime}$ $\displaystyle=$
$\displaystyle\mathbf{a}_{\text{in}}+\frac{2\
m_{\text{out}}}{m_{\text{in}}+m_{\text{out}}}\mathbf{a}_{\perp}^{b}$ (23)

where $\mathbf{a}_{\perp}^{b}$ is normal to the boundary surface, since it is
is parallel to $\mathbf{r}_{\text{in,s}}$ which defines the boundary radius,
and connects the solute to the particles $\alpha_{\text{b}}$ and
$\beta_{\text{b}}$ which coincide after the rotation. In the case of Velocity
Verlet this component is reduced to $\mathbf{v}_{\perp}^{b}$.

After the collision, $\mathbf{a}_{\text{out}}^{\prime\text{\ rot}}$ is rotated
back to the outer particle reference frame

$\mathbf{a}_{\text{out}}^{\prime}=X(\mathbf{c},2\pi-\phi)\cdot\mathbf{a}_{\text{out}}^{\prime\text{\
rot}}$ (24)

and the system’s propagating vector field is updated such that

$\\{\mathbf{A}^{i}(\alpha_{b})=\mathbf{a}^{\prime}_{\text{in}},\mathbf{A}^{i}(\beta_{b})=\mathbf{a}_{\text{out}}^{\prime}\\}\to\mathbf{A}^{i}_{c}$
(25)

which implies

$\\{\mathbf{V}^{i+\frac{1}{2}}(\alpha_{b})=\mathbf{v}^{\prime}_{\text{in}},\mathbf{V}^{i+\frac{1}{2}}(\beta_{b})=\mathbf{v}_{\text{out}}^{\prime}\\}\to\mathbf{V}^{i+\frac{1}{2}}_{c}$
(26)

since the components of equation (12) are additive and both propagate the
coordinate space linearly in terms of $dt$, such that the exchange as outlined
above can be applied separately on each component.

The system coordinate space is evolved by the remainder of the time step

$\mathbf{R}^{i+1}_{c}=\mathbf{R}^{i+\delta}+(dt-
dt^{\prime})\mathbf{A}^{i}_{c}$ (27)

At this stage the updated coordinate space is checked for additional
conflicts, first by updating the boundary radius, equation (3), resulting in
an update of $\alpha_{b}$, followed by a check if any outer particles have
crossed the new boundary, equation (4), resulting in an update of
$S^{\prime}_{\text{out}}$. Multiple conflicts are discussed in section 2.4.1.

If no additional conflicts are detected in the interval $(dt-dt^{\prime})$,
the SAFIRES propagator triggers a force calculation for the new configuration
$\mathbf{F}(\mathbf{R}^{i+1})$ and the velocity vectors are updated to
$\mathbf{V}^{i+1}$ according to equation 9, applied to
$\mathbf{V}_{c}^{i+\frac{1}{2}}$. This completes the update of all vector
fields.

###  2.4 Special cases

####  2.4.1 Multiple conflicts during the same time step

Multiple conflicts do not add any additional complexity, rather, the same
procedures are applied as outlined in the preceding section, using the
$\mathbf{R}^{i+\delta}$ and $\mathbf{A}^{i}_{c}$ vector fields as a starting
point. A new partial time step, $dt^{\prime\prime}$, is solved for following
equations (14)–(16), and is $dt^{\prime\prime}\leq(dt-dt^{\prime})$, where the
norm of the vector field
$\mathbf{R}^{i+\delta}_{s}=\mathbf{R}^{i+\delta}-\mathbf{R}^{i+\delta}(s)$ is
now considered. Namely the condition

$\displaystyle||\mathbf{R}^{i+\delta+\epsilon}_{s}(\alpha_{b})||$
$\displaystyle\overset{!}{=}$
$\displaystyle||\mathbf{R}^{i+\delta+\epsilon}_{s}(\beta)||\ \ \forall\
\beta\in S^{\prime\prime}_{\text{out}}\;,$ (28)

where $\epsilon$ is defined as $dt^{\prime\prime}/(dt-dt^{\prime})$, is
imposed and solved for using

$||\mathbf{R}_{s}^{i+\delta}(\alpha_{b})+\delta
t^{\prime\prime}(\beta)\mathbf{A}^{i}_{c}(\alpha_{b})||=||\mathbf{R}^{i+\delta}_{s}(\beta)+\delta
t^{\prime\prime}(\beta)\mathbf{A}^{i}_{c}(\beta)||\ \ \forall\ \beta\in
S_{\text{out}}$ (29)

and

$\beta_{b}\leftarrow dt^{\prime\prime}=\text{min}\\{\delta
t^{\prime\prime}(\beta)\mid\forall\ \delta
t^{\prime\prime}(\beta)\in\mathbb{R}_{+}\\}$ (30)

where $dt^{\prime\prime}$ now brings the temporally first conflict,
corresponding to particles $\alpha_{b}$ and $\beta_{b}$, to the boundary
radius as defined by the propagation of the vector field

$\mathbf{R}^{i+\delta+\epsilon}=\mathbf{R}^{i+\delta}+dt^{\prime\prime}\mathbf{A}_{\text{c}}^{i}$
(31)

An elastic collision is again performed using the vectors
$\mathbf{a}^{i}_{\text{c,in}}$ and $\mathbf{a}^{i}_{\text{c,out}}$ (i.e.
$\mathbf{A}_{\text{c}}^{i}(\alpha_{b})$ and
$\mathbf{A}_{\text{c}}^{i}(\beta_{b})$ respectively).

This process of resolving multiple conflicts and solving for partial time
steps in smaller and smaller remaining time intervals can in principle be
continued indefinitely. However, in the practical examples presented in this
work (see Section X) we find it is a rare event for conflicts to exceed one,
even if multiple possible outer particles are identified through equation (4).
This occurs since the resolution of the temporally first conflict often
results in the resolution of other possible conflicts within the time step
interval due to the change in the trajectory of $\alpha_{\text{b}}$.

####  2.4.2  Semi-infinite Surfaces

In the case of a 2D-periodic system the SAFIRES ’check for conflicts’ and
’elastic collision’ steps are simplified. Given an orthorhombic left-handed
Cartesian coordinate system with origin $O=(0,0,0)$ and periodicity along the
$x$\- and $y$-axis, an $xy$-plane boundary is defined as the numerically
largest $z$ position of the inner particles. The condition imposed on the
indexing of the inner and outer particles in the initial position vector space
is then as follows

$\\{\mathbf{R}_{z}^{i}(\alpha)\ \forall\alpha\in S_{\text{in}}\\}<\
\\{\mathbf{R}_{z}^{i}(\beta)\ \forall\beta\in S_{\text{out}}\\}$ (32)

such that all outer particles are above the $xy$-plane boundary. Here
$\mathbf{R}_{z}^{i}(\alpha)$ is shorthand for $\mathbf{R}^{i}(\alpha)(0,0,z)$,
i.e. the z-coordinate of particle $\alpha$. After an integrated time step of
the MD simulation the resulting particle positions are used to identify the
outermost inner particle

$\alpha_{b}\leftarrow r_{z}=\text{max}\left\\{\mathbf{R}_{z}^{i+1}(\alpha)\
\forall\alpha\in S_{\text{in}}\right\\}$ (33)

and the check for conflicts is then simply

$S^{\prime}_{\text{out}}=\left\\{\mathbf{R}_{z}^{i+1}(\beta)<r_{z}\
\forall\beta\in S_{\text{out}}\right\\}$ (34)

The partial time step required to propagate the conflicting particles to the
boundary is then solved for according to

$\displaystyle\mathbf{R}^{i+\delta}_{z}(\alpha_{b})$
$\displaystyle\overset{!}{=}$ $\displaystyle\mathbf{R}^{i+\delta}_{z}(\beta)\
\ \forall\beta\in S^{\prime}_{\text{out}}.$ (35)

resulting in

$\mathbf{R}_{z}^{i}(\alpha_{b})+\delta
t^{\prime}(\beta)\mathbf{A}_{z}^{i}(\alpha_{b})=\mathbf{R}_{z}^{i}(\beta)+\delta
t^{\prime}(\beta)\mathbf{A}_{z}^{i}(\beta)\ \ \forall\beta\in
S^{\prime}_{\text{out}}$ (36)

At the boundary the collision between the temporally first in/out particle
pair results in an exchange of the $z$-components of the propagating vectors.
The exchange is

$\displaystyle\mathbf{a}_{\perp}^{b}=$ $\displaystyle=$
$\displaystyle(\mathbf{a}_{\text{out}}-\mathbf{a}_{\text{in}})_{z}$ (37)
$\displaystyle\mathbf{a}_{\text{out},z}^{\prime}$ $\displaystyle=$
$\displaystyle\mathbf{a}_{\text{out},z}-\frac{2\
m_{\text{in}}}{m_{\text{in}}+m_{\text{out}}}\mathbf{a}_{z}$ (38)
$\displaystyle\mathbf{a}_{\text{in},z}^{\prime}$ $\displaystyle=$
$\displaystyle\mathbf{a}_{\text{in},z}+\frac{2\
m_{\text{out}}}{m_{\text{in}}+m_{\text{out}}}\mathbf{a}_{z}$ (39)

followed by an update of the propagating vector fields and the coordinate
space in the same manner as described in the preceding sections.

No rotation operation is required in the case of an asymmetric surface model –
where solvent is in contact with the surface one side – such as the one
studied in this work, see Section 4.2. In the case of a symmetric surface
model with solvent particles in contact with the surface on both sides,
$\mathbf{a}_{\text{out}}$ needs to be inverted if the conflicting in/out pair
$\alpha_{\text{b}}$ and $\beta_{\text{b}}$ are located on opposite sides of
the surface model. $\mathbf{a}_{\text{out}}$ is then inverted once more after
an exchange of the propagating vectors.

####  2.4.3 Multiatomic particles

Applications require multiatomic particles to be handled as well. An example
are the rigid TIPnP family of water force fields.41, 42, 43 In such cases
SAFIRES works with reduced vector field spaces, corresponding to the position
and momentum at the center of mass for each molecular species. Effectively,
each molecule is considered as a pseudo-particle and conflicts are monitored
relative to the center of mass for each pseudo-particle, not individual atoms.
In the case of a system of $n_{\text{mol}}$ rigid water molecules the vector
fields transform from $3n$\- to $3n_{\text{mol}}$-dimensional fields, which
are for the position and velocity

$\displaystyle\mathbf{R}^{i}_{\text{cm}}=\left\\{\sum_{a}^{S_{\text{mol}}}m_{a}\mathbf{R}^{i}(a)/M\
\forall S_{\text{mol}}\in S\right\\}$ (40)
$\displaystyle\mathbf{V}^{i}_{\text{cm}}=\left\\{\sum_{a}^{S_{\text{mol}}}m_{a}\mathbf{V}^{i}(a)/M\
\forall S_{\text{mol}}\in S\right\\}$ (41)

where the subset $S_{\text{mol}}=\\{\ce{O},\text{H1},\text{H2}\\}$ indexes the
atoms of the each molecule in the set $S$ of all atom indices. The force
vector field is reduced as

$\mathbf{F}_{\text{cm}}(\mathbf{R}_{\mathrm{cm}}^{i})=\left\\{\sum_{a}^{S_{\text{mol}}}\mathbf{F}(\mathbf{R}^{i}(a))\
\forall S_{\text{mol}}\in S\right\\}$ (42)

$m_{a}$ and $M$ is the mass of the atom and molecule, respectively.
Additionally, for the Langevin integrator the random variable vector fields
transform as

$\mathbf{\Sigma}^{i}_{\text{cm}}=\left\\{\sum_{a}^{S_{\text{mol}}}\sqrt{m_{a}}\mathbf{\Sigma}^{i}(a)/M\
\forall S_{\text{mol}}\in S\right\\}$ (43)

and similarly for $\mathbf{\Lambda}^{i}_{\text{cm}}$. The square root
dependence of the per-atom mass-weighing follows from the multiplicative
factor $\bm{\sigma}$ as defined in equation (10). Using the center of mass
transformed fields above, the corresponding propagating field,
$\mathbf{A}^{i}_{\text{cm}}$, is evaluated with equation (12).

The SAFIRES processes operate on the reduced dimensional vector fields and
corresponding molecular index set in such a way that the ’checking for
conflicts’, equations (3)-(4), followed by ’time step extrapolation’,
equations (15)-(16) result in identifying conflicting in/out pseudo-particle
pair indices $\alpha_{b}$ and $\beta_{b}$, as well as index sets
$S_{\text{mol}}^{\alpha_{b}}=\\{\alpha^{\text{O}}_{b},\alpha^{\text{H1}}_{b},\alpha^{\text{H2}}_{b}\\}$
and
$S_{\text{mol}}^{\beta_{b}}=\\{\beta^{\text{O}}_{b},\beta^{\text{H1}}_{b},\beta^{\text{H2}}_{b}\\}$.
That is, for a conflicting in/out pseudo-particle pair the molecular indices,
and corresponding atomic indices are identified.

With the molecular index identified the SAFIRES ’elastic collision’, as
described in Section 2.3, results in an exchange of center of mass normal
components. For example, given conflicting psuedo-particle pairs and
corresponding propagating fields $\mathbf{a}^{i}_{\text{in,cm}}$ =
$\mathbf{A}_{\text{cm}}^{i}(\alpha_{b})$ and
$\mathbf{a}^{i}_{\text{out,cm}}=\mathbf{A}_{\text{cm}}^{i}(\beta_{b})$ the
resulting center of mass components become

$\displaystyle\mathbf{a}_{\text{out,cm}}^{\prime\text{\ rot}}$
$\displaystyle=$
$\displaystyle\mathbf{a}_{\text{out,cm}}^{\text{rot}}-\frac{2\
M_{\text{in}}}{M_{\text{in}}+M_{\text{out}}}\mathbf{a}_{\perp,\text{cm}}^{b}$
(44) $\displaystyle\mathbf{a}_{\text{in,cm}}^{\prime}$ $\displaystyle=$
$\displaystyle\mathbf{a}_{\text{in,cm}}+\frac{2\
M_{\text{out}}}{M_{\text{in}}+M_{\text{out}}}\mathbf{a}_{\perp,\text{cm}}^{b}$
(45)

After rotation of the outer pseudo-particle propagating vector – to give
$\mathbf{a}_{\text{out,cm}}^{\prime}$ – the components need to be
redistributed to the atoms. The tangential components of the atoms in each
pseudo-particle are not affected in the collision, hence the redistribution is
done as follows

$\displaystyle\mathbf{A}^{i}_{c}(\alpha)=\mathbf{a}_{\parallel,\text{in}}(\alpha)+\frac{m_{\alpha}}{M}(\mathbf{a}^{\prime}_{\text{in,cm}})_{\perp}\
\forall\ \alpha\in S_{\text{mol}}^{\alpha_{b}}$ (46)
$\displaystyle\mathbf{A}^{i}_{c}(\beta)=\mathbf{a}_{\parallel,\text{out}}(\beta)+\frac{m_{\beta}}{M}(\mathbf{a}^{\prime}_{\text{out,cm}})_{\perp}\
\forall\ \beta\in S_{\text{mol}}^{\beta_{b}}$ (47)

where $\mathbf{a}_{\parallel}$ are the unperturbed tangential components. This
way any tangential modes of the molecule are not quenched in the process.
While this procedure does conserve momentum, some energy in the rotational
degrees of freedom is transformed to energy in the translational degrees of
freedom.

##  3 Computational Details

SAFIRES and FIRES are implemented within the Python-based framework of the
Atomic Simulation Environment (ASE) available under the GNU Lesser General
Public License44. A force constant of $k_{\text{FIRES}}=500$ kcal A-2 is used
for calculations using FIRES.

Test calculations are performed using the LJ potential available in ASE with
simulation parameters for argon ($\sigma=3.4$ Å, $\epsilon=120\ k_{\text{B}}$,
$\rho=1.374$ g cm-3, $T=94.4$ K)45. For MD simulations, the Velocity-Verlet
and Langevin dynamics implemented in ASE are utilized. The $g(r)$ of LJ
systems are sampled over 1 ns, using an NVE ensemble and a time step of 1 fs
with the Velocity Verlet propagator. Starting configurations for LJ liquid
simulations are pre-equilibrated over 1 ns of NVT Langevin dynamics using the
same time step.

In case of the LJ liquid simulation, 25 of the total 512 LJ particles are
included in the inner region, including the fixed central particle (”solute”)
for simulations using FIRES or SAFIRES. For the LJ surface model, three layers
á 16 particles are cut in 111 direction from the most stable fcc Ar crystal
configuration46. The Ar surface model is frozen during simulations. The inner
and outer regions contain 32 and 96 LJ particles in the liquid state,
respectively.

The TIP4P force field is used for water simulations41. The $g(r)$ for water-
in-water calculations are sampled every 1 ps for a total of 20 ns, using an
NVT ensemble and the Langevin propagator with a time step of 0.5 fs. A
friction coefficient of 0.05 is used. Starting configurations for water-in-
water simulations are equilibrated over 250 ps in an NVT ensemble, using the
Langevin propagator with a time step of 2 fs and a friction coefficient of
0.01. The first 20 ps of each run are discarded before sampling. RATTLE
constraints are used to ensure rigid bond lengths and bond angles of the water
molecules47. For calculations with FIRES and SAFIRES, 14 molecules are
included in the inner region, including the fixed central molecule (”solute”).
All pair distributions are sampled using the VMD program48.

##  4 Test calculations

###  4.1 Lennard-Jones liquid

To test SAFIRES, an LJ liquid with argon parameters is used as the model
system. This simple model allows efficient benchmarking of the technical
implementation. Figure 3 depicts a cross section of the liquid with different
colors for the solute (pink), particles in the inner region (blue), and
particles in the outer region (grey).

Figure 3: Illustration of a cross section of the LJ liquid computational model
used in this work. Pink: solute, blue: particles in the inner region, grey:
particles in the outer region.

Both inner and outer region are calculated using the same LJ potential.

First, energy conservation of the SAFIRES algorithm is tested for NVE
Velocity-Verlet dynamics and this model system using the approach presented by
Allen and Tildesley49. It is found that energy conservation with SAFIRES is
significantly improved over FIRES, see Figure S1.

Next, the influence of the ensemble separator on the simulation is quantified.
To this end, the $g(r)$ are calculated from 1 ns of NVE Velocity-Verlet
dynamics each. Velocity-Verlet is used over Langevin to keep the first test as
simple as possible and to monitor energy conservation. Only pairs involving
the solute are considered so that $r$ correlates to the location of the
boundary with respect to the solute. The $g(r)$ for an unconstrained reference
simulation as well as for simulations with FIRES and SAFIRES are shown in
Figure 4.

Figure 4: $g(r)$ of a LJ liquid using argon parameters. Pairs are sampled
between a fixed central LJ particle (”solute”) and the surrounding LJ
particles. Black dotted line: reference calculation without any ensemble
separation; red: FIRES; blue: SAFIRES; dashed light grey: normalized
probability distribution of the SAFIRES boundary location; dark grey:
normalized probability distribution of the FIRES boundary location.

Also shown in Figure 4 are normalized probability distributions of the
locations of the FIRES and SAFIRES boundaries.

SAFIRES reproduces the unconstrained $g(r)$ without significant deviations.
Notably, FIRES introduces accumulation artifacts around the boundary location.
The boundary location distributions as well as the artifact around the
boundary observed with FIRES broaden for a LJ liquid of lower density, see
Figure S2. This means that the artifact can potentially be obfuscated under
certain simulation conditions, as speculated by Bulo and co-workers33. The
observed improvement in case of SAFIRES is likely related to its instantaneous
resolution of conflicts. When using FIRES, particles experience a spring force
when passing the boundary, accelerating them either away or towards the
solute. This means, they will spend several iterations in the other respective
region, being first decelerated and then accelerated again in the opposite
direction. It is suggested that the spring force approach will lead to an
accumulation of density in the $g(r)$ around the border, even for a simplistic
model system like the LJ fluid. Analogous results are obtained with an NVT
ensemble and Langevin dynamics (see Figure S3).

###  4.2 Lennard-Jones surface

SAFIRES is built with interface simulations in mind. To illustrate this
application, a LJ liquid using argon parameters is placed in a simulation cell
with a 111-indexed surface made up of three layers of frozen LJ argon
particles. Figure 5 illustrates the model system.

Figure 5: Illustration of the periodic surface computational model used in
this work. Pink: frozen surface layer, blue: particles in the inner region,
grey: particles in the outer region. Dotted line indicates the periodically
repeating unit. Note that vacuum separates periodic images in z direction.

The vacuum region above this asymmetric surface slab used to avoid interaction
with the next periodic image in z direction is not shown.

The $g(z)$ is calculated from 1 ns of NVE Velocity-Verlet dynamics using this
model. Since the solute in this case is a periodic surface, pairs are
calculated between a xy plane located in the top layer of the surface and the
particles that constitute the liquid. This way, the origin of the abscissa in
Figure 6 coincides with the top layer of the surface model. Figure 6 shows the
resulting $g(z)$ of an unconstrained simulation and of two simulation using
SAFIRES, one simulation with 16 LJ particles in the inner region (”one layer”)
and another simulation with 32 particles in the inner region (”two layers”).

Figure 6: $g(z)$ of a LJ liquid using argon parameters on top of a three-layer
(111)-indexed surface of frozen LJ argon particles. ”One layer”: 16 LJ liquid
particles in the inner region; ”two layers”: 32 LJ particles in the inner
region. Pairs are sampled between an xy plane located in the top layer of the
surface model and the LJ particles that constitute the liquid. Black dotted
line: unconstrained reference calculation; blue: SAFIRES with 16 particles in
the inner region; grey dashed line: normalized probability distribution of the
SAFIRES boundary location corresponding to 16 particles in the inner region;
grey dashed-dotted line: normalized probability distribution of the SAFIRES
boundary location corresponding to 32 particles in the inner region.

Simulations using SAFIRES reproduce the reference RDF without significant
deviations, both in case of 16 or 32 particles within the inner region. This
highlights the robustness of the method against the composition of inner and
outer region.

In case of the simulation with 32 atoms in the inner region (”two layers”),
the normalized probability distribution of the SAFIRES boundary is slightly
broadened and asymmetric compared to the liquid presented in Figure 4. This
difference in shape is a result of the generally lower density of the liquid
in this case as it can expand into the vacuum above the surface.

For the simulation with 16 atoms in the inner region, the normalized
probability distribution of the SAFIRES boundary is significantly broadened as
a result of interface effects. Analyzing the trajectories reveals that the
first layer of LJ liquid particles ontop the surface behaves significantly
more orderly than the remaining liquid, similar to the ice-like water layers
on a Pt(111) surface in contact with liquid water. This leads to on average
larger particle distances within this first layer and between the first and
second layers and thus results in a broadening of the boundary location
distribution in the z direction perpendicular to the surface.

###  4.3 Simulation of water-in-water

Finally, exemplary MM*/MM simulations of water-in-water using the TIP4P force
field and an NVT ensemble are performed using FIRES and SAFIRES to separate
the inner and outer regions. The geometry of water molecules is kept rigid as
the present implementation of TIP4P does not support flexible molecules. The
$g(r)$ of a reference calculation without ensemble separation as well as of
simulations using FIRES and SAFIRES are shown in Figure 7.

Figure 7: $g(r)$ of TIP4P water. Pairs are sampled between a fixed central
water molecule (”solute”) and the surrounding water molecules. Black dotted
line: unconstrained reference calculation; red: FIRES; blue: SAFIRES; dashed
light grey and dark grey: normalized probability distribution of the SAFIRES
and FIRES boundary location, respectively.

SAFIRES reproduces the unconstrained $g(r)$ without significant deviations.
The $g(r)$ obtained with FIRES shows an accumulation artifact around the
boundary region, similar to results obtained for the LJ liquid in Figure 4.

##  5 Conclusion

The scattering-adapted flexible inner region ensemble separator (SAFIRES) is a
separation algorithm for hybrid simulations coupling different methodologies.
It has been designed in particular for interface calculations. Like other
boundary-based methods, SAFIRES is built on the premise that if the same type
of particle is present in the inner and outer region, correct average ensemble
statistics can be obtained despite particles being unable to cross the
boundary and exchange between regions. However, unlike other boundary-based
methods that use repulsive forces or bias potentials to keep particles
separated, SAFIRES instantaneously redirects particles at the boundary via
energy- and momentum-conserving elastic collisions. This approach therefore
constitutes the most rigorous implementation of the original flexible boundary
premise by Rowley and Roux to date29. To ensure exact collisions, SAFIRES
introduces a new, multiple-time-step propagator which reduces to the Vanden-
Eijnden / Ciccotti implementation of a Langevin propagator for constant time
steps and to the Velocity-Verlet algorithm for constant time steps and zero
friction.

Using a LJ liquid and a LJ liquid in contact with a surface as exemplary
systems, SAFIRES reproduces unconstrained reference $g(r)$ or $g(z)$,
respectively, without significant deviations while FIRES introduces
accumulation artifacts around the boundary region. Lastly, an exemplary MM*/MM
study of water-in-water using the TIP4P force field is presented. The $g(r)$
obtained with SAFIRES does not significantly deviate from the unconstrained
reference calculation while accumulation artifacts are present in a simulation
using FIRES. SAFIRES will be available as part of the open-source ASE package.

Future development will focus on the inclusion of virtual counter ions in
SAFIRES to model the electrode potential near a surface within a Poisson-
Boltzmann scheme as well as combining SAFIRES with polarizable QM/MM coupling
models50 for electrochemical applications. SAFIRES is also not restricted to
liquid-liquid or solid-liquid interfaces; it could for example be used to
study diffusion in solids for battery research. Future work will therefore
also focus on exploring various interface phenomena using this approach.

BK thanks the University of Iceland Research Fund for funding through a PhD
fellowship. Computations were performed on resources provided by the Icelandic
High Performance Computing Centre at the University of Iceland. This work was
supported by the Icelandic Research Fund and by the DFG (German Science
Foundation) through the Collaborative Research Center SFB-1316 as well as the
cluster of excellence POLiS (project ID 390874152).

The results of the energy conservation test, a RDF of a Lennard-Jones liquid
at lower density, a RDF of a Lennard-Jones liquid sampled using an NVT
ensemble, and pseudocode of the SAFIRES algorithm are presented in the
electronic Supporting Information.

## References

  * Kitanosono et al. 2018 Kitanosono, T.; Masuda, K.; Xu, P.; Kobayashi, S. Catalytic Organic Reactions in Water toward Sustainable Society. _Chem. Rev._ 2018, _118_ , 679–746
  * Gould et al. 2020 Gould, N. S.; Li, S.; Cho, H. J.; Landfield, H.; Caratzoulas, S.; Vlachos, D.; Bai, P.; Xu, B. Understanding solvent effects on adsorption and protonation in porous catalysts. _Nat. Commun._ 2020, _11_ , 1060
  * Mellmer et al. 2018 Mellmer, M. A.; Sanpitakseree, C.; Demir, B.; Bai, P.; Ma, K.; Neurock, M.; Dumesic, J. A. Solvent-enabled control of reactivity for liquid-phase reactions of biomass-derived compounds. _Nat. Catal._ 2018, _1_ , 199–207
  * Román-Leshkov et al. 2006 Román-Leshkov, Y.; Chheda, J. N.; Dumesic, J. A. Phase Modifiers Promote Efficient Production of Hydroxymethylfurfural from Fructose. _Science_ 2006, _312_ , 1933–1937
  * Rossin et al. 2006 Rossin, A.; Kovács, G.; Ujaque, G.; Lledós, A.; Joó, F. The Active Role of the Water Solvent in the Regioselective CO Hydrogenation of Unsaturated Aldehydes by [RuH2(mtppms)x] in Basic Media. _Organometallics_ 2006, _25_ , 5010–5023
  * Warzok et al. 2018 Warzok, U.; Marianski, M.; Hoffmann, W.; Turunen, L.; Rissanen, K.; Pagel, K.; A. Schalley, C. Surprising solvent-induced structural rearrangements in large [NI + N] halogen-bonded supramolecular capsules: an ion mobility-mass spectrometry study. _Chem. Sci._ 2018, _9_ , 8343–8351
  * Yang et al. 2019 Yang, K.; Chen, X.; Zheng, Z.; Wan, J.; Feng, M.; Yu, Y. Solvent-induced surface disorder and doping-induced lattice distortion in anatase TiO2 nanocrystals for enhanced photoreversible color switching. _J. Mater. Chem. A_ 2019, _7_ , 3863–3873
  * Hodel and Luber 2016 Hodel, F. H.; Luber, S. Redox-Inert Cations Enhancing Water Oxidation Activity: The Crucial Role of Flexibility. _ACS Catal._ 2016, _6_ , 6750–6761
  * Abel et al. 2008 Abel, R.; Young, T.; Farid, R.; Berne, B. J.; Friesner, R. A. Role of the Active-Site Solvent in the Thermodynamics of Factor Xa Ligand Binding. _J. Am. Chem. Soc._ 2008, _130_ , 2817–2831
  * Burnham and English 2019 Burnham, C. J.; English, N. J. Crystal Structure Prediction via Basin-Hopping Global Optimization Employing Tiny Periodic Simulation Cells, with Application to Water–Ice. _J. Chem. Theory Comput._ 2019, _15_ , 3889–3900
  * Zhang and Dolg 2016 Zhang, J.; Dolg, M. Global optimization of clusters of rigid molecules using the artificial bee colony algorithm. _Phys. Chem. Chem. Phys._ 2016, _18_ , 3003–3010
  * Reda et al. 2018 Reda, M.; Hansen, H. A.; Vegge, T. DFT study of stabilization effects on N-doped graphene for ORR catalysis. _Catal. Today_ 2018, _312_ , 118–125
  * Heenen et al. 2020 Heenen, H. H.; Gauthier, J. A.; Kristoffersen, H. H.; Ludwig, T.; Chan, K. Solvation at metal/water interfaces: An ab initio molecular dynamics benchmark of common computational approaches. _J. Chem. Phys._ 2020, _152_ , 144703
  * E. Skyner et al. 2015 E. Skyner, R.; L. McDonagh, J.; R. Groom, C.; Mourik, T. v.; O. Mitchell, J. B. A review of methods for the calculation of solution free energies and the modelling of systems in solution. _Phys. Chem. Chem. Phys._ 2015, _17_ , 6174–6191
  * Gray et al. 2017 Gray, C. M.; Saravanan, K.; Wang, G.; Keith, J. A. Quantifying solvation energies at solid/liquid interfaces using continuum solvation methods. _Mol. Simul._ 2017, _43_ , 420–427
  * Zhang et al. 2017 Zhang, J.; Zhang, H.; Wu, T.; Wang, Q.; van der Spoel, D. Comparison of Implicit and Explicit Solvent Models for the Calculation of Solvation Free Energy in Organic Solvents. _J. Chem. Theory Comput._ 2017, _13_ , 1034–1043
  * Warshel and Levitt 1976 Warshel, A.; Levitt, M. Theoretical studies of enzymic reactions: Dielectric, electrostatic and steric stabilization of the carbonium ion in the reaction of lysozyme. _J. Mol. Biol._ 1976, _103_ , 227–249
  * Thole and van Duijnen 1980 Thole, B. T.; van Duijnen, P. T. On the quantum mechanical treatment of solvent effects. _Theor. Chim. Acta_ 1980, _55_ , 307–318
  * Field et al. 1990 Field, M. J.; Bash, P. A.; Karplus, M. A combined quantum mechanical and molecular mechanical potential for molecular dynamics simulations. _J. Comp. Chem._ 1990, _11_ , 700–733
  * Waller et al. 2014 Waller, M. P.; Kumbhar, S.; Yang, J. A Density-Based Adaptive Quantum Mechanical/Molecular Mechanical Method. _ChemPhysChem_ 2014, _15_ , 3218–3225
  * Bernstein et al. 2012 Bernstein, N.; Várnai, C.; Solt, I.; Winfield, S. A.; Payne, M. C.; Simon, I.; Fuxreiter, M.; Csányi, G. QM/MM simulation of liquid water with an adaptive quantum region. _Phys. Chem. Chem. Phys._ 2012, _14_ , 646–656
  * Heyden et al. 2007 Heyden, A.; Lin, H.; Truhlar, D. G. Adaptive Partitioning in Combined Quantum Mechanical and Molecular Mechanical Calculations of Potential Energy Functions for Multiscale Simulations. _J. Phys. Chem. B_ 2007, _111_ , 2231–2241
  * Pezeshki and Lin 2015 Pezeshki, S.; Lin, H. Adaptive-Partitioning QM/MM for Molecular Dynamics Simulations: 4. Proton Hopping in Bulk Water. _J. Chem. Theory Comput._ 2015, _11_ , 2398–2411
  * Bulo et al. 2009 Bulo, R. E.; Ensing, B.; Sikkema, J.; Visscher, L. Toward a Practical Method for Adaptive QM/MM Simulations. _J. Chem. Theory Comput._ 2009, _5_ , 2212–2221
  * Field 2017 Field, M. J. An Algorithm for Adaptive QC/MM Simulations. _J. Chem. Theory Comput._ 2017, _13_ , 2342–2351
  * Watanabe et al. 2014 Watanabe, H. C.; Kubař, T.; Elstner, M. Size-Consistent Multipartitioning QM/MM: A Stable and Efficient Adaptive QM/MM Method. _J. Chem. Theory Comput._ 2014, _10_ , 4242–4252
  * Watanabe 2018 Watanabe, H. C. Improvement of performance, stability and continuity by modified size-consistent multipartitioning quantum mechanical/molecular mechanical method. _Molecules_ 2018, _23_ , 1882
  * Watanabe and Cui 2019 Watanabe, H. C.; Cui, Q. Quantitative Analysis of QM/MM Boundary Artifacts and Correction in Adaptive QM/MM Simulations. _J. Chem. Theory Comput._ 2019, _15_ , 3917–3928
  * Rowley and Roux 2012 Rowley, C. N.; Roux, B. The Solvation Structure of Na ${}^{\textrm{+}}$ and K ${}^{\textrm{+}}$ in Liquid Water Determined from High Level ab Initio Molecular Dynamics Simulations. _J. Chem. Theory Comput._ 2012, _8_ , 3526–3535
  * Beglov and Roux 1994 Beglov, D.; Roux, B. Finite representation of an infinite bulk system: Solvent boundary potential for computer simulations. _J. Chem. Phys._ 1994, _100_ , 9050–9063, Publisher: American Institute of Physics
  * Lu et al. 2015 Lu, X.; Gaus, M.; Elstner, M.; Cui, Q. Parametrization of DFTB3/3OB for Magnesium and Zinc for Chemical and Biological Applications. _J. Phys. Chem. B_ 2015, _119_ , 1062–1082
  * Boereboom et al. 2018 Boereboom, J. M.; Fleurat-Lessard, P.; Bulo, R. E. Explicit Solvation Matters: Performance of QM/MM Solvation Models in Nucleophilic Addition. _J. Chem. Theory Comput._ 2018, _14_ , 1841–1852
  * Bulo et al. 2013 Bulo, R. E.; Michel, C.; Fleurat-Lessard, P.; Sautet, P. Multiscale Modeling of Chemistry in Water: Are We There Yet? _J. Chem. Theory Comput._ 2013, _9_ , 5567–5577
  * Shiga and Masia 2013 Shiga, M.; Masia, M. Boundary based on exchange symmetry theory for multilevel simulations. I. Basic theory. _J. Chem. Phys._ 2013, _139_ , 044120
  * Shiga and Masia 2013 Shiga, M.; Masia, M. Erratum: “Boundary based on exchange symmetry theory for multilevel simulations. I. Basic theory” [J. Chem. Phys. 139, 044120 (2013)]. _J. Chem. Phys._ 2013, _139_ , 119901
  * Takahashi et al. 2018 Takahashi, H.; Kambe, H.; Morita, A. A simple and effective solution to the constrained QM/MM simulations. _J. Chem. Phys._ 2018, _148_ , 134119
  * Vanden-Eijnden and Ciccotti 2006 Vanden-Eijnden, E.; Ciccotti, G. Second-order integrators for Langevin equations with holonomic constraints. _Chem. Phys. Lett._ 2006, _429_ , 310–316
  * Daru et al. 2019 Daru, J.; Gupta, P. K.; Marx, D. Restricting Solvation to Two Dimensions: Soft Landing of Microsolvated Ions on Inert Surfaces. _J. Phys. Chem. Lett._ 2019, _10_ , 831–835
  * Agrafonov et al. 2015 Agrafonov, Y.; Petrushin, I.; Damdinov, B.; Tsydypov, S. Radial distribution function for liquid near the solid surface. _Proc. Mtgs. Acoust._ 2015, _24_ , 045002
  * Vincenty 1975 Vincenty, T. Direct and Inverse Solutions of Geodesics on the Ellipsoid with Application of Nested Equations. _Surv. Rev._ 1975, _23_ , 88–93
  * Abascal and Vega 2005 Abascal, J. L. F.; Vega, C. A general purpose model for the condensed phases of water: TIP4P/2005. _J. Chem. Phys._ 2005, _123_ , 234505
  * Jorgensen et al. 1983 Jorgensen, W. L.; Chandrasekhar, J.; Madura, J. D.; Impey, R. W.; Klein, M. L. Comparison of simple potential functions for simulating liquid water. _J. Chem. Phys._ 1983, _79_ , 926–935
  * Mahoney and Jorgensen 2000 Mahoney, M. W.; Jorgensen, W. L. A five-site model for liquid water and the reproduction of the density anomaly by rigid, nonpolarizable potential functions. _J. Chem. Phys._ 2000, _112_ , 8910–8922
  * Larsen et al. 2017 Larsen, A. H.; Mortensen, J. J.; Blomqvist, J.; Castelli, I. E.; Christensen, R.; Dulak, M.; Friis, J.; Groves, M. N.; Hammer, B.; Hargus, C.; Hermes, E. D.; Jennings, P. C.; Jensen, P. B.; Kermode, J.; Kitchin, J. R.; Kolsbjerg, E. L.; Kubal, J.; Kaasbjerg, K.; Lysgaard, S.; Maronsson, J. B.; Maxson, T.; Olsen, T.; Pastewka, L.; Peterson, A.; Rostgaard, C.; Schiøtz, J.; Schütt, O.; Strange, M.; Thygesen, K. S.; Vegge, T.; Vilhelmsen, L.; Walter, M.; Zeng, Z.; Jacobsen, K. W. The atomic simulation environment—a Python library for working with atoms. _J. Phys.: Condens. Matter_ 2017, _29_ , 273002
  * Rahman 1964 Rahman, A. Correlations in the Motion of Atoms in Liquid Argon. _Phys. Rev._ 1964, _136_ , A405–A411
  * Barrett and Meyer 1965 Barrett, C. S.; Meyer, L. The Crystal Structures of Argon and Its Alloys. Low Temperature Physics LT9. Boston, MA, 1965; pp 1085–1088
  * Andersen 1983 Andersen, H. C. Rattle: A “velocity” version of the shake algorithm for molecular dynamics calculations. _J. Comp. Phys._ 1983, _52_ , 24–34
  * Humphrey et al. 1996 Humphrey, W.; Dalke, A.; Schulten, K. VMD: Visual molecular dynamics. _J. Mol. Graph._ 1996, _14_ , 33–38
  * Allen and Tildesley 2017 Allen, M. P.; Tildesley, D. J. _Computer simulation of liquids_ , second edition ed.; Oxford University Press: Oxford, United Kingdom, 2017
  * Dohn et al. 2019 Dohn, A. O.; Jónsson, E. Ö.; Jónsson, H. Polarizable Embedding with a Transferable H2O Potential Function II: Application to (H2O)n Clusters and Liquid Water. _J. Chem. Theory Comput._ 2019, _15_ , 6578–6587