# MASKE: A kinetic simulator of coupled chemo-mechanical processes driving
microstructural evolution

Enrico Masoero

###### Abstract

The microstructure of materials evolves through chemical reactions and
mechanical stress, often strongly coupled in phenomena such as pressure
solution or crystallization pressure. This article presents MASKE: a simulator
to address the challenge of modelling coupled chemo-mechanical processes in
microstructures. MASKE represents solid phases as agglomerations of particles
whose off-lattice displacements generate mechanical stress through interaction
potentials. Particle precipitation and dissolution are sampled using Kinetic
Monte Carlo, with original reaction rate equations derived from Transition
State Theory and featuring contributions from mechanical interactions.
Molecules in solution around the solid are modelled implicitly, through
concentrations that change during microstructural evolution and define the
saturation indexes for user-defined chemical reactions. The structure and
implementation of the software are explained first. Then, two examples on a
nanocrystal of calcium hydroxide address its chemical equilibrium and its
mechanical response under a range of imposed strain rates, involving stress-
driven dissolution and recrystallization. These examples highlight MASKE’s
distinctive ability to simulate strongly coupled chemo-mechanical processes.
MASKE is available, open-source, on GitHub.

###### keywords:

A. chemo-mechanical processes , C. numerical algorithms , A. microstructures ,
Kinetic Monte Carlo.

††journal: …

organization=Department of Civil and Environmental Engineering, Politecnico di
Milano,addressline=Piazza Leonardo da Vinci 32, city=Milan, postcode=20133,
state=MI, country=Italy

MASKE is presented: off-lattice Kinetic Monte Carlo on discretized
microstructures

Original rate equations hard-coupling chemical reactions with mechanical
stress

Dissolution/growth of a calcium hydroxide nanocrystal is simulated over long
time

The crystal is put under realistic strain rates to compute stress-strain
curves

A complex balance between solution chemistry, pressure solution, and plastic
flow

##  1 Introduction

A strong coupling between mechanical stress and chemical reactions often
underlies the microstructural evolution of materials. Examples of stress-
driven nucleation and growth of materials are commonly found in nature, e.g.
rock metamorphism or the adaptive growth of bones (Hart et al., 2017) and
plants (Kouhen et al., 2023), but such processes are also exploited in nano-
manufacturing, e.g. the heterogeneous precipitation of islands on top of
mismatching crystals, which is a challenge in semiconductor production
(Dubrovskii et al., 2010). Further examples concern the degradation of
materials, such as crystallization pressure (Scherer, 2004), pressure solution
(viz. pressure-induced dissolution, Rutter (1983)), creep from the relaxation
of eigenstress originated during the formation of concrete (Bažant et al.,
1997), and stress-corrosion cracking (Raja and Shoji, 2011).

In Solid Mechanics, the modelling of chemo-mechanical processes is quite
established at the macroscopic continuum level, with fully coupled
formulations of reactive transport and linear momentum conservation, e.g. in
Gawin et al. (2009). At the constitutive level, these formulations sometimes
consider the stress or strain stemming from phase transformations, e.g. local
expansion during alkali-silica reaction and sulphate attack in concrete
(Pesavento et al., 2012; Cefis and Comi, 2017), whereas the rates of chemical
transformation are rarely stress-dependent (Bellomo et al., 2012). Chemo-
mechanical modelling is less advanced at the smaller microstructural level,
where the focus in on the evolving morphology of single, or multiple,
crystalline and amorphous phases. The length and time scales of such
microstructural processes are too large for direct atomistic simulations,
hence appropriate coarse-graining approaches are needed. However, most current
simulations of microstructural evolution either describe the geometry of
growing/dissolving solids without addressing their stress field, or they focus
on the mechanical interactions prompting structural aggregation and collapse,
but considering chemical reactions only qualitatively. Among the former, with
focus on chemistry, there are microstructural development codes such as in Van
Breugel (1995); Bishnoi and Scrivener (2009); Bullard et al. (2011), as well
as others simulators based on cellular automata (Bentz et al., 1994; Bullard
et al., 2010) or on Kinetic Monte Carlo (KMC) (Martin et al., 2020). Examples
of simulators with a focus on mechanics are instead those in Ioannidou et al.
(2014, 2016, 2022), based on interacting particles. Exceptions exist, but with
limitations. For example, (Li et al., 2015) and (Feng et al., 2017) staggered
finite element simulations of mechanical response with a reactive module to
simulate solid precipitation, respectively of stress-free cement hydrates
during creep and of expanding ettringite during sulphate attack in a cement
paste. In both cases, however, chemo-mechanical coupling was not intrinsic to
the model, but imposed by alternating chemical and mechanical modules.
Differently, Meakin and Rosso (2008) employed stress-dependent reaction rates
in KMC simulations of crystal dissolution near stress-inducing dislocations.
However the solid in these simulations was bound to a fixed lattice geometry
(a cubic Kossel crystal), hence the stress field near the dislocation was pre-
imposed and non-relaxing during dissolution. Another noteworthy method is the
Phase Field, which does offer a framework to fully couple chemistry and
mechanics; however, its high computational cost is preventing application to
complex 3D structures, with multiple phases, liquid and solid, undergoing
nonlinear mechanical processes (Tourret et al., 2022).

This work presents MASKE: a simulator of fluid-solid microstructures
discretized via mechanically interacting particles whose dissolution and
precipitation are sampled using off-lattice Kinetic Monte Carlo. MASKE allows
the user to define multiple chemical reactions, along with their usual
thermodynamic and kinetic parameters such as equilibrium constants and
activation energies. Chemo-mechanical coupling is introduced through rate
equations for the chemical reactions that feature the change in mechanical
interaction energy accompanying particle dissolution or precipitation.
Precursors of these rate equations where used in previous works. Shvab et al.
(2017) captured the experimental hydration rate curve of a cement paste by
simulating the nucleation and stress-driven agglomeration of calcium-silicate-
hydrate nanoparticles. Coopamootoo and Masoero (2020) simulated the
experimentally measured critical undersaturation of a tricalcim silicate (C3S)
paste, below which stress-driven stepwave dissolution is initiated at screw
dislocations; recently, Coopamootoo K (2023) have extended the result to
finite-sized crystal, capturing the full experimental relationship between
dissolution rate and saturation index of C3S. Alex et al. (2023) applied the
coupled rate equations to the carbonation of a model cement paste, featuring
multiple solid phases, ions in solution, and chemical reactions. These
microscale simulations provided apparent rate constants that, employed in
macroscale simulations of carbonation, matched a range of experimental results
(Freeman et al., 2019). All these applications showed that the coupled
reaction rates in MASKE can capture chemo-mechanical processes at
experimentally realistic length and time scales. However, a detailed
presentation and discussion of the architecture, capabilities, and limitations
of MASKE are still due.

The algorithm and software architecture of MASKE are described in detail,
accompanying its recent release, open-source, on GitHub (Masoero, 2023a). The
main article focuses on the most important parts, concept,s and steps in the
algorithm; much more details are provided in several appendices, which guide
the interested reader to a deeper understanding of the method, while also
describing the input, output, and some data structures of the software. In the
main article, particular attention is paid to: (i) the stress-dependent rates
of chemical reaction implemented in MASKE, and their relationship with
Transition State Theory (TST) and its predictions; (ii) the KMC algorithm to
sample particle dissolution and off-lattice precipitation; (iii) the
possibility to integrate continuous kinetic processes in time along with
discrete KMC events. These key features are demonstrated via two sample
applications, both considering a spheroidal nano-crystal of calcium hydroxide
Ca(OH)2 in an aqueous solution of its ions. The chemical kinetics of Ca(OH)2
under mechanical stress is of practical importance, as the instability and
dissolution of this phase is responsible for pH drop and steel corrosion in
reinforced concrete, and contribute to carbonation creep in cement pastes
(Parrott, 1975). Here however the focus is not on experimental validation,
which has been already addressed in part by the applications summarised in the
previous paragraph. Instead, the aim here is to showcase the three key
features of MASKE listed above. Hence, the first application considers the
chemical equilibrium, dissolution, and growth of the Ca(OH)2 crystal in
stress-free conditions; the results are compared with expectations from
classical TST, and clarify the effect of some parameters in the rate equations
and in the off-lattice sampling of precipitation. The second example puts the
crystal under compression for a wide range of imposed strain rates; the
results show that the interplay between mechanical deformations and stress-
driven dissolution and recrystallization trigger rich and sometimes counter-
intuitive mechanical responses, which are nevertheless explained by analysing
the evolution of local stresses during the tests.

##  2 Methodology

This section presents for the first time the main structure and salient
features of MASKE: a C++, object-oriented, parallel software that combines
continuous processes, which are explicitly integrated in time, with discrete
events, such as dissolution/precipitation reactions, sampled using off-lattice
Kinetic Monte Carlo (KMC). In the Appendices, the interested reader can find
full details of the MASKE algorithm, its implementation, the commands in its
input script, and the parameters in its chemical database input file. In this
section, attention is also paid to mechanical interactions between particles
that are later used in two sample applications on Ca(OH)2 nanocrystals. The
methodological details of these applications are provided here too.

###  2.1 Model description

A typical simulation features spherical particles in a prismatic simulation
box with volume $V_{box}$ and periodic or fixed boundary conditions: see Fig.
1. Box and particles are generated and tracked using the LAMMPS C++ library
interface (Thompson et al., 2022). The particles may represent any explicit
phase, hereafter assumed to be solid ones. The particles interact mechanically
with each other through user-defined interaction potentials, which are also
defined and managed in LAMMPS: e.g. the $U_{ij}(r_{ij})$ interaction in Fig.
1. The current version of MASKE can only manage pairwise interactions, but
extensions to any of the potentials in LAMMPS would be quite straightforward:
see A for more discussion on this. Particle displacements caused by the
mechanical interactions are also computed in LAMMPS, via usual methods such as
molecular dynamics or energy minimization (respectively using the $run$ and
$minimize$ commands in LAMMPS). This interface with LAMMPS gives the user
flexibility to study complex systems with various mechanical interaction
potentials and loading conditions.

Figure 1: A typical simulation box in MASKE, with an implicit solution and
solid particles interacing via effective potentials (in this example, pairwise
energy $U_{ij}$ as a function of interparticle distance $r_{ij}$).

The portion of a simulation box that is not occupied by solid particles
contains a volume $voidV$ of empty space, and a volume $V_{sol}$ of implicit
solution made up of various, user-defined molecular species. The current
version of MASKE assumes that the molar concentrations of the solvated species
are uniform in $V_{sol}$, with no diffusion taking place. There is also an
additional volume $\Delta V_{box}$ attached to the simulation box, which can
receive some of the molecules released during dissolution reactions, or
provide some of the molecules consumed during precipitation reactions. This
$\Delta V_{box}$ can be used to keep $V_{box}$ limited, which is where time-
consuming operations take place; for example, Fig. 2 depicts a scenario where
this approach allows focussing on dissolution and precipitation near the
surface of a large grain, while preserving realistic concentrations in
solution and surface-to-volume ratio of suspended grains at a larger scale.
The initial molar concentrations of all the solvated molecules in $V_{box}$
and in $\Delta V_{box}$ are provided by the user in the input script, along
with fixing the temperature of the solution and the $A$ and $B$ solvent
parameters to be used in Debye-Hückel calculations of activity coefficients.

Figure 2: Rationale of adding a solution volume $\Delta V_{box}$, connected to
the simulation box of volume $V_{box}$.

MASKE simulates discrete events, which happen once every $\Delta t_{k}$ units
of time (where $\Delta t_{k}$ varies during a simulation depending on
mechanical interactions and solution composition), but also continuous
processes that are numerically integrated over a user-specified time interval
$\Delta t_{i}$ with time step $dt_{i}$ (for the $i^{th}$ continuous process).
In the current version of MASKE, the possible discrete events are full
dissolution (called $deletion$) of any existing particle, or full
precipitation (called $nucleation$) of a new particle at any possible location
in $V_{box}$. The only continuous process implemented so far is called
$CmstoreLMP$: it repeatedly executes a set of stored LAMMPS commands,
specified by the user in the MASKE input script (see C for more details on the
input/output structure of the code). A $CmstoreLMP$ process will be used in an
example here, to impose a constant strain rate. The occurrence of discrete
events and continuous processes is staggered onto a unique time line, $t^{*}$,
which is advanced by whichever event or process has the shortest time of
future occurrence, or $FOT$. The details of how this unique time line is
managed can be found in B. The occurrence of discrete events is sampled using
Kinetic Monte Carlo (KMC), thus the list of all possible discrete events is
associated with a unique $FOT_{k}$, which is the future occurrence time of any
one out of all the possible discrete events. If $FOT_{k}$ is smaller than
$FOT_{i}$ of all the other continuous processes, i.e. if a discrete event is
next in line to be executed, then one discrete event is picked out from the
list of all the possible ones, with a probability that is proportional to its
individual rate. This is the usual approach in KMC simulations; specifically,
MASKE uses a rejection-free KMC algorithm that is time-dependent (Prados et
al., 1997), to account for changes in solution composition or stress state
that continuous processes may induce between two successive discrete events.
The equations of time-dependent KMC are also presented and discussed in B.

In the KMC method, $FOT_{k}$ is obtained from the cumulative rate:

$R_{k}=\sum_{j}^{N+M}r_{k,j}$ (1)

which is the sum of the individual rates $r_{k,i}$ of all $N$ possible
particle deletions plus all $M$ particle nucleations at all possible locations
in the simulation box. The expressions for $r_{k,j}$ are presented in the next
section. $N$ is finite, whereas the number of possible locations for particle
precipitation is infinite and, in principle, one should sample them all to use
KMC. To treat this problem, MASKE asks the user to select a region within
$V_{box}$ and to create a 3D lattice of $M$ trial particles of the desired
type, viz. with desired chemical composition and interaction potentials. Each
of these $M$ trial particles is then driven towards its closest local minimum
of interaction energy with the other already-existing, so-called $real$
particles (no interactions between trial particles are considered, and the
positions of the real particles are kept fixed while moving the trial ones).
This energy minimization is performed using a bespoke algorithm, called
$quickmaske$, which the $quickmin$ algorithm in LAMMPS except that each trial
particles is independently relaxed. The fineness of the user-specified lattice
determines the volume $\Delta V$ of the generic lattice cell. The energy
minimum obtained for a trial particle is assumed to characterize all the
particles that may nucleate anywhere else in $\Delta V$. In this way, the
cumulative rate of all the $M$ sampled nucleation events is effectively a
numerical integral in the sampled region, with resolution $\Delta V$. The most
accurate result is obtained when $\Delta V\rightarrow 0$, but in practice a
small but finite $\Delta V$ is already sufficient to estimate the cumulative
nucleation rate accurately. How small $\Delta V$ should be, it depends on the
employed interaction potentials and on the morphology of the simulated
microstructure. An assumption underling this approach is that the energy
minima alone fully control the precipitation rate cumulated over the whole
$V_{box}$. In reality, thermal fluctuations would add contributions also from
the locations surrounding the minima; future versions of MASKE may sample such
fluctuations by relaxing the trial particles using molecular dynamics instead
of energy minimization. F offers more discussion on interaction energy
landscapes and their sampling for dissolution and precipitation events.

The flowchart in Fig. 3 depicts the algorithm for the kinetic simulations in
MASKE. A brief description is given here, while more details can be found in
E. The left side of the flowchart summarises the sampling of $FOT$ for all the
discrete and continuous processes in a simulation. The shortest $FOT$
determines which type of process or event is to be executed next. If a
continuous process is next, it is directly integrated numerically over its
time period $\Delta t_{i}$ (user-specified in the input script). After that,
the individual rates $r_{k,j}$ of the sampled discrete events are multiplied
times $\Delta t_{i}$ and these products are stored in a per-particle array
called $ratesT$, of size $N+M$. The $ratesT$ array is used in the next kinetic
step for computing $FOT_{k}$. If the next event to execute is instead a
discrete one, the $ratesT$ values are immediately updated, adding to each
entry the product $r_{k,j}(FOT_{k}-t^{*})$. The $ratesT$ values are then used
to select which discrete event to carry out specifically, out of all the
possible ones. After executing a discrete event, the concentrations of the
solvated molecules are updated accordingly, as discussed in the next section,
and the $ratesT$ values are reset to zero. Subsequently, irrespective of
whether the executed process is continuous or discrete, MASKE restores the
initial position of all trial particles and executes a third type of possible
processes, called $Every$. Unlike continuous or discrete processes, whose
occurrence is associated to a time line, processes of type $Every$ are simply
executed once every $n$ kinetic steps, where a step is one full loop in Fig.
3. $n$ is user-specified as an input, and presently MASKE implements only one
type of $Every$ process, called $EmstoreLMP$, which executes a set of LAMMPS
command that the user can specify and store in the input script. A typical use
of such an $Every$ process is to execute a full energy minimization of the
system with $n=1$, i.e. after each execution of a continuous or discrete
process.

Figure 3: Flow chart of the main kinetic loop in MASKE.

The algorithm in Fig. 3 is well suited to parallelisation, since each discrete
event (i.e. each of the $N$ particle deletions and $M$ particle nucleations)
can be sampled independently. Furthermore, spatial domain decomposition can be
exploited during energy minimization. To exploit this, MASKE features a two-
level parallelisation. On the first level, the user can specify multiple
groups of processors, called $subcommunicators$, and attribute the sampling of
one or more processes to each subcommunicator (e.g. one subcommunicator
samples particle deletion, another one samples nucleation). On the second
level, each subcommunicator uses all its processors to run any invoked LAMMPS
command, e.g. energy calculations and minimisation using domain decomposition.
This article does not report on the scalability and efficiency of the
parallelisation, which is amenable to optimisation in the current version of
MASKE. However, previously published works have already exploited both levels
of parallelisation to simulate hundreds of thousands particles (Coopamootoo
and Masoero, 2020) and multiple types of solid particles, all sampled for
deletion and nucleation (Alex et al., 2023).

###  2.2 Chemical reactions: rates and mechanisms

MASKE computes the rate of any discrete event, such as particle deletion or
nucleation, based on two user-defined entities: a chemical reaction and a
reaction mechanism, both specified in a chemical database text file, called
$ChemDB$, whose syntax is presented in D. The chemical reaction is assumed to
be one-step and is defined primarily through its formula and stoichiometric
coefficients, e.g. :

$\displaystyle Ca(OH)_{2}\rightarrow Ca^{2+}+2OH^{-}$ (2)

for the dissolution of a calcium hydroxide molecule. The stoichiometric
coefficients are specified in the $ChemDB$ file. The molecules appearing in
the reaction formulas are also user-defined in the $ChemDB$ file, along with a
set of useful parameters. For the $i^{th}$ solvated molecule, these parameters
are its apparent volume $v_{app,i}$, hydrated radius $a^{o}_{i}$, ionic charge
$z_{i}$, and the $b_{i}$ parameter for Debye-Hückel calculations of activity
coefficients, discussed later. For solid molecular species, only linear sizes
are required, which are used to compute their molecular volumes. $ChemDB$
gathers also the thermodynamic and kinetic parameters associated to each
chemical reaction: standard free energy barriers $\Delta G^{o{\ddagger}}$,
standard concentrations of activated complexes $c^{o{\ddagger}}$ (in the sense
of Transition State Theory, which will be used below to express rate
equations), equilibrium constants $K_{eq}$, and interfacial energies $\gamma$.
See D for more details on this quantities.

The $ChemDB$ file also contains the definition of reaction mechanisms, which
are the sets of assumptions that MASKE makes when assembling single chemical
reactions to delete or nucleate a full particle. Indeed, a particle may be a
coarse-grained entity whose deletion or creation may entail more than just a
single reaction, e.g. $n_{r,j,diss}$ or $n_{r,j,nuc}$ reactions respectively
to dissolve or nucleate a particle. Classical Nucleation Theory is an example
of multiple individual reactions being assembled to express a reaction rate at
the particle scale. In the current version of MASKE, the only implemented
mechanism is called $allpar$ and it assumes that, irrespective of how big a
particle is, its full deletion or nucleation only takes the time for one
single-step chemical reaction to occur. This is exact if the volume $V_{p,j}$
of the generic particle $j$ equals the sum of the volumes of solid molecules
dissolving or precipitating in one reaction. For bigger particles instead, the
$allpar$ mechanism effectively assumes that the same chemical reaction takes
places everywhere in $V_{p,j}$ at the same time. This is clearly unrealistic
for particles beyond the nanometre or so, therefore in this article all the
examples will feature particles whose $V_{p,j}$ entail
$n_{r,j,diss}=n_{r,j,nuc}=1$. Hereafter this resolution is called
’unimolecular’, although in principle a single chemical reaction might also
involve multiple solid molecules at once.

Whenever a particle is deleted or nucleated, a chemical reaction is carried
out as part of an $allpar$ mechanisms. MASKE ensures mass conservation by
consistently managing molecule exchanges between solid phases and solution,
following the user-specified reaction formulas in $ChemDB$. Chemical reactions
may entail differences in volume between reactants and products; these are
accommodated in the $voidV$ volume introduced in the previous section, which
may thus be positive or negative. The current version of MASKE does not
simulate chemical reactions taking place in the implicit solution, hence it
does not minimize the free energy of the solution itself. An interface between
MASKE and the thermodynamic simulator PHREEQC (Parkhurst, 1995) is currently
being developed, to address this limitation.

The rate equations for the $allpar$ mechanis in MASKE are rooted in Transition
State Theory (TST), but with an extension to account for excess enthalpy
coming from the mechanical interactions between particles. This excess
enthalpy term, originally proposed in Shvab et al. (2017), is essential for
capturing the effect of local morphology and stress state on reaction rates.
The rate equations for one-step dissolution and precipitation reactions are:

$\displaystyle r_{1,d}$
$\displaystyle=\kappa\frac{k_{B}T}{h}\frac{c^{o{\ddagger}}_{d}}{\gamma^{\ddagger}_{d}}\exp\left(\frac{-\Delta
G^{o{\ddagger}}_{d}}{k_{B}T}\right)\exp\left[-(1-\chi)\frac{\Delta
U_{d}-\gamma\Omega}{k_{B}T}\cdot\frac{1}{n_{r,j,diss}}\right]Q_{r,d}V_{t}^{\alpha/3}$
(3) $\displaystyle r_{1,pr}$
$\displaystyle=\kappa\frac{k_{B}T}{h}\frac{c^{o{\ddagger}}_{pr}}{\gamma^{\ddagger}_{pr}}\exp\left(\frac{-\Delta
G^{o{\ddagger}}_{pr}}{k_{B}T}\right)\exp\left[-\chi\frac{\Delta
U_{pr}+\gamma\Omega}{k_{B}T}\cdot\frac{1}{n_{r,j,nuc}}\right]Q_{r,pr}\Delta
V\cdot V_{t}^{\alpha/3-1}$ (4)

Eqs. 3 and 4 will be hereafter called ’straight’ rate equations, as opposed to
’net’ rates defined later in this section. $\kappa$ is the probability of an
activated complex evolving into reaction products; MASKE makes the common
assumption that $\kappa=1$ (Lasaga, 2014). $k_{B}$ and $h$ are the Boltzmann
and Planck constants, and $T$ is the temperature in degrees Kelvin.
$c^{o{\ddagger}}_{d}$ and $c^{o{\ddagger}}_{pr}$ are the standard
concentrations of activated complexes for the dissolution and precipitation
reactions, $\gamma^{{\ddagger}}_{d}$ and $\gamma^{o{\ddagger}}_{pr}$ are their
activity coefficients, and $\Delta G^{o{\ddagger}}_{d}$ and $\Delta
G^{o{\ddagger}}_{pr}$ are their standard free energies of activation.
$c^{o{\ddagger}}$ and $\Delta G^{o{\ddagger}}$ can be estimated
experimentally, noting that the whole prefactors in Eqs. 3 and 4 defines the
rate constant $k$:

$k=\kappa\frac{k_{B}T}{h}\frac{c^{o{\ddagger}}}{\gamma^{\ddagger}_{d}}\exp\left(\frac{-\Delta
G^{o{\ddagger}}}{k_{B}T}\right)$ (5)

Eq. 5 implies that the arbitrary choice of $c^{o{\ddagger}}$ must be
compensated by the value of $\Delta G^{o{\ddagger}}$, so that a physical,
univocally defined rate constant is preserved. The details of this
compensation are discussed in F.

The $V_{t}$ term in Eqs. 3 and 4 is the tributary volume of the particle being
sampled for deletion or nucleation, viz. the volume of the interaction energy
basin, whose local energy minimum is where the sampled particle sits. The user
of MASKE must specify $V_{t}$ in the $ChemDB$ file; for solid phases governed
by short-range interactions, $V_{t}$ is typically in the range of one particle
volume (see F for more discussion and examples of $V_{t}$). $\alpha$ is the
spatial dimensionality of $c^{o{\ddagger}}_{d}$ or $c^{o{\ddagger}}_{pr}$;
namely, $\alpha=3$ if the reactions are per unit volume or, as typical for
dissolution/precipitation, $\alpha=2$ indicating reactions per unit surface.
$\Delta V$ is the cell volume of the user-defined lattice of trial particles
sampling candidate nucleation sites. When the KMC algorithm sums all the
nucleartion rates for trial particles pertaining to the same interaction
energy basin, the sum of all $\Delta V$’s equals $V_{t}$ (within the limits
posed by the finite, and not infinitesimal, value of $\Delta V$). As a result,
$V_{t}^{\alpha/3}$ in Eqs. 3 equals the sum of all $\Delta V\cdot
V_{t}^{\alpha/3-1}$ pertaining to the same interaction energy basin. This
aspect is discussed in more detail in G.

The square brackets in Eqs. 3 and 4 quantify the excess enthalpy of a
particle, resulting from two contributions: $\Delta U_{d}$, which is the
change in total interaction energy in the system caused by the particle
deletion ($\Delta U_{pr}$ is the analogous term for a nucleation even),
and$\gamma$, which is the solid-solution interfacial energy attributed to the
sampled particle, and user-defined in the $ChemDB$ file. The $\Omega$ term,
which multiplies $\gamma$ in the excess enthalpy, is the surface area of the
sampled particle. $n_{r,j,diss}$ and $n_{r,j,nuc}$ appear in the excess
enthalpy because the $allpar$ mechanism may be used for particles whose
deletion or nucleation entails multiple chemical reactions; for the
unimolecular resolution in this article, $n_{r,j,diss}=n_{r,j,nuc}=1$. The
user-defined parameter $\chi$, with value between 0 and 0.5, specifies which
fraction of excess enthalpy $U+\gamma\Omega$ is still present in the system
when the reaction is in its activated complex state. $\chi=0$ is the usual
assumption in TST, but results in this article will show that $\chi>0$ may
sometimes improve the efficiency of a simulation.

The last terms to be defined in the rate equations are $Q_{r,d}$ and
$Q_{r,pr}$, which are the activity products of the reactants in the
dissolution and precipitation reactions. If $i$ is the generic solvated
reactant:

$Q_{r}=\prod_{i}(\gamma_{i}c_{i})^{\nu_{i}}$ (6)

where $c_{i}$ is the concentration in solution of the molecular species $i$,
$\nu_{i}$ is its stoichiometric coefficient in the reaction, and $\gamma_{i}$
is its activity coefficient. MASKE computes $\gamma_{i}$ using the WATEQ
Debye-Hückel formula (Truesdell and Jones, 1974):

$log_{10}(\gamma_{i})=\frac{-z_{i}^{2}A\sqrt{I}}{1+Ba_{i}^{o}\sqrt{I}}+b_{i}I$
(7)

$z_{i}$ is the charge of the $i^{th}$ species. $A$ and $B$ are solvent-
specific constants (0.51 and 3.29 nm-1 for water). $a_{i}^{o}$ is the hydrated
radius of the molecule in solution, and $b_{i}$ is a molecule-specific
dimensionless constant; both these parameters are provided by the user in the
$ChemDB$ file. $I$ is the ionic strength of the fluid, featuring $n_{s}$
solvated species:

$I=\sum_{j=1}^{n_{s}}c_{j}z_{j}^{2}$ (8)

When $a_{i}^{o}$ or $b_{i}$ are unknown, MASKE computes the activity
coefficients using established simplified methods (Langmuir, 1997), detailed
in F.

The straight rate expressions in Eqs. 3 and 4 sample all the fluctuations
between particle deletion and nucleation that, on average, guide the overall
morphology evolution of the system. However, many such fluctuations may cancel
out, especially near equilibrium, leading to inefficient simulations where the
overall morphology makes very little progress on average. In such scenarios it
may be convenient to remove the fluctuations by considering net rates instead,
defined as the difference between straight rates:

$\displaystyle r_{1,d}^{net}$ $\displaystyle=\max\langle
0,r_{1,d}-r_{1,pr}\rangle=$ $\displaystyle=\max\left\langle
0,\kappa\frac{k_{B}T}{h}\frac{c^{o{\ddagger}}_{d}}{\gamma^{\ddagger}_{d}}\exp\left(\frac{-\Delta
G^{o{\ddagger}}_{d}}{k_{B}T}\right)V_{t}^{\alpha/3}\cdot\right.$
$\displaystyle\hskip 42.67912pt\left.\left\\{\exp\left[-(1-\chi)\frac{\Delta
U_{d}-\gamma\Omega}{k_{B}T}\cdot\frac{1}{n_{r,j,diss}}\right]Q_{r,d}-\exp\left[-\chi\frac{-\Delta
U_{d}+\gamma\Omega}{k_{B}T}\cdot\frac{1}{n_{r,j,diss}}\right]\frac{Q_{p,d}}{K_{eq,d}}\right\\}\right\rangle$
(9) $\displaystyle r_{1,pr}^{net}$ $\displaystyle=\max\langle
0,r_{1,pr}-r_{1,d}\rangle=$ $\displaystyle=\max\left\langle
0,\kappa\frac{k_{B}T}{h}\frac{c^{o{\ddagger}}_{pr}}{\gamma^{\ddagger}_{pr}}\exp\left(\frac{-\Delta
G^{o{\ddagger}}_{pr}}{k_{B}T}\right)\Delta V\cdot
V_{t}^{\alpha/3-1}\cdot\right.$ $\displaystyle\hskip
28.45274pt\left.\left\\{\exp\left[-\chi\frac{\Delta
U_{pr}+\gamma\Omega}{k_{B}T}\cdot\frac{1}{n_{r,j,nuc}}\right]Q_{r,pr}-\exp\left[-(1-\chi)\frac{-\Delta
U_{pr}-\gamma\Omega}{k_{B}T}\cdot\frac{1}{n_{r,j,nuc}}\right]\frac{Q_{p,pr}}{K_{eq,pr}}\right\\}\right\rangle$
(10)

These net rate equations are also implemented in MASKE, as a possible option
for the $allpar$ mechanism. Their derivation is detailed in F. $Q_{p,d}$ and
$Q_{p,pr}$ are the activity products of the products of the dissolution and
precipitation reactions, with $K_{eq,d}$ and $K_{eq,pr}$ their equilibrium
constants. The reactants of a dissolution reaction and the products of a
precipitation reaction are solid phases, conventionally in standard state when
stress-free, thus $Q_{r,d}=Q_{p,pr}=1$. G shows how the net rates in Eqs. 9
and 10 effectively scale linearly with the saturation index $\beta$ of the
solution for the associated chemical reaction; namely,
$r_{1,d}^{net}\sim(1-\beta)$ and $r_{1,pr}^{net}\sim(\beta-1)$, as usual in
Transition State Theory. The definition of $\beta$ is:

$\beta=\frac{Q_{r,pr}}{K_{eq,d}}=\frac{Q_{p,d}}{K_{eq,d}}$ (11)

noting that both the reactants of precipitation and the products of
dissolution are simply the solvated species participating in the reaction. In
a scenario where there is no excess enthalpy from changes in mechanical stress
and interfacial energies, $\beta=1$ would imply equilibrium, i.e. no
dissolution nor precipitation taking place, whereas $\beta>1$ would favour
precipitation and $\beta<1$ dissolution.

###  2.3 Interaction potentials between identical particles

The energy $U$ from mechanical interactions plays an important role in the
reaction rates in MASKE, via the excess enthalpy terms $\Delta
U_{d}-\gamma\Omega$ and $\Delta U_{pr}+\gamma\Omega$ in Eqs. 3 and 4, as
implemented in the $allpar$ mechanism. All the examples in this article will
employ a truncated harmonic potential of interaction between pairs of
identical, spherical particles, $i$ and $j$:

$\displaystyle U_{ij}$
$\displaystyle=\frac{1}{2}k_{0}\left(r_{ij}-r_{0}\right)^{2}-\varepsilon_{0}\;\;\;\mathrm{if}\;r_{ij}<r_{c}$
$\displaystyle U_{ij}$ $\displaystyle=0\;\;\;\mathrm{if}\;r_{ij}\geq r_{c}$
(12)

Figure 4: Truncated harmonic potential between two identical particles $i$ and
$j$, with interparticle distance $r_{ij}$, equilibrium distance $r_{0}$ with
minimum $U_{ij}=-\varepsilon_{0}$, and interaction cutoff $r_{c}$.

plotted in Fig. 4. $r_{ij}$ is the interparticle distance. When the particles
are at their equilibrium distance $r_{0}$, the interaction energy $U_{ij}$ is
minimum and equal to the bond energy $-\varepsilon_{0}$. The stiffness of the
interaction is dictated by the $k_{0}$ parameter. The cutoff distance for the
interaction, $r_{c}$, is set here as the distance at which $U_{ij}$ returns to
zero when $r_{ij}>r_{0}$; because of this definition, $r_{c}$ is not an
independent parameter, but it is fixed once $r_{0}$, $k_{0}$, and
$\varepsilon_{0}$ are decided:

$r_{c}=r_{0}+\sqrt{\frac{2\varepsilon_{0}}{k_{0}}}$ (13)

$k_{0}$ is obtained from the Young modulus $E$ of the material of the
particles:

$k_{0}=\pi\frac{E(D_{0}/2)^{2}}{r_{0}}$ (14)

The particle diameter $D_{0}$ is close but not equal to $r_{0}$, and their
relationship depends on the modelling choice of the ‘limit crystal lattice’ of
the solid phase. Depending on the type and range of the employed interaction
potential, there are only several ordered arrangements, viz. crystal lattices,
that are mechanically stable and where all the interacting particles sit at
equilibrium distances. For instance, limit crystal lattices for the short-
ranged, first-neighbour-only, radial interactions in Fig. 4, would be face
centred cubic (FCC) or hexagonal close packing (HCP), both with $n_{b}=12$
first neighbours in the bulk. The mismatch between $r_{0}$ and $D_{0}$ emerges
because any limit crystal lattice of spherical particles features a porosity
and fills the space only up to a volume fraction $\eta<1$. For example, both
FCC and HPC lattices have packing fraction $\eta=0.74$, hence a 26% porosity.
However, if these limit lattices are used to discretise a non-porous continuum
in a simulation, the maximum packing $\eta$ should represent a solid with same
density as a non-porous phase, viz. with the same number of molecules, or
mass, per unit volume. To exactly compensate for the artificial porosity
$1-\eta$, Coopamootoo and Masoero (2020) have proposed to use an $r_{0}$
smaller than $D_{0}$ by:

$r_{0}=D\sqrt[3]{\eta}$ (15)

The expressions of excess enthalpy in Eqs. 3 and 4 ($\Delta U\pm\gamma\Omega$)
impose a dependence between interaction potential $\varepsilon_{0}$ and solid-
solution interfacial energy $\gamma$, which must be respected in order to
produce physically meaningful simulations. Consider a particle in kink
position in the crystal lattice, which by definition has $n_{b}/2$ first
neighbours. Precipitating or dissolving a stress-free particle in a kink
position should not change the excess enthalpy of the system, because it
creates the same amount of solid-solution interfacial area as it removes. Such
a dissolution (and similarly, in reverse, for precipitation) can be modelled,
first, as a detachment of the kink particle from the solid, which causes a
loss of $n_{b}/2$ inter-particle bonds each with energy $\varepsilon_{0}$
(accounted for in the $\Delta U$ term), and second, as a dissolution of the
detached particle which removes $\gamma\Omega$ interfacial energy. This leads
to:

$\varepsilon_{0}\frac{n_{b}}{2}=\gamma\Omega\;\;\rightarrow\;\;\varepsilon_{0}=\frac{2\gamma\Omega}{n_{b}}$
(16)

###  2.4 Interaction potentials between different particles

One example in this article will feature mechanical interactions between
particles of two different types, with different sizes and representing
different materials. For such scenarios, Alex et al. (2023) have proposed an
extension of the truncated harmonic potential in Fig. 4:

$\displaystyle U_{ij}$
$\displaystyle=\frac{1}{2}k_{12}\left(r_{ij}-r_{12}\right)^{2}-\varepsilon_{12}\;\;\;\mathrm{if}\;r_{ij}<r_{c,12}$
$\displaystyle U_{ij}$ $\displaystyle=0\;\;\;\mathrm{if}\;r_{ij}\geq r_{c,12}$
(17)

The equilibrium distance $r_{12}=(r_{0,1}+r_{0,2})/2$ is the average of the
equilibrium distances between identical particles for the two different types,
1 and 2. The stiffness term is expressed as
$k_{12}=2\frac{k_{1}k_{2}}{k_{1}+k_{2}}$, i.e. idealizing the interacting
particles as two springs in series, with individual stiffness as per Eq. 14.
The bond energy is
$\varepsilon_{12}=\frac{\gamma_{12}(\Omega_{1}+\Omega_{2})}{n_{b}}$, with
$\gamma_{12}$ the solid-solid interfacial energy between the particles and
$\Omega_{1}$ and $\Omega_{2}$ their individual surface areas. Eq. 13 is still
valid for the cutoff distance, using $r_{12}$, $\varepsilon_{12}$, and
$k_{12}$ in it.

###  2.5 Example 1: dissolution, growth, and equilibrium of a stress-free
crystal

A nanocrystal of calcium hydroxide, Ca(OH)2, is considered in this example.
The crystal is immersed in an aqueous solution of Ca2+ and OH- ions at room
temperature $T=298$ K. The crystal is made of 4,055 spherical particles, each
with the volume of one Ca(OH)2 molecule,
$V_{p,CH}=\frac{V_{M,CH}}{N_{a}}=0.0557$ nm3, where $N_{a}$ is the Avogadro
number and $V_{M,CH}=33.53$ cm3/mol is the molar volume of calcium hydroxide
(indicated as $CH$ in the subscript). This brings to a particle diameter
$D_{CH}=0.4737$ nm. The particles are initially arranged on an FCC lattice
with cell size $a=\sqrt{2}\;r_{0,CH}$, which implies that $r_{0,CH}$ is the
distance between first-neighbour particles in the lattice. Given the packing
density $\eta=0.74$ of the FCC lattice, Eq. 15 provides
$r_{0,CH}=\sqrt[3]{\eta}\;D_{CH}=0.4285$ nm.

Figure 5: (a) Initial configuration of the Ca(OH)2 nano-crystal, obtain from
pre-dissolving a cubic crystal of Ca(OH)2 particles in FCC arrangement; (b) A
section of the same nano-crystal surrounded by trial particles, used to sample
precipitation of new particles. The trial particles are shown after
minimization to their closest local energy minima, hence those in the inner
part of the spherical shell region leave the perfect cubic arrangement of the
initial lattice, driven by mechanical interactions with the other, already
existing Ca(OH)2 particles (in red). Instead, the trial particles in the outer
part of the spherical shell do not interact with any existing particle, hence
they maintain the ordered cubic arrangement of the trial lattice.

The interaction potential between particles is the truncated harmonic one from
Eq. 12. From the Young modulus $E_{CH}=48$ GPa (Wittmann, 1986) of calcium
hydroxide, Eq. 14 provides $k_{0,CH}=19,745$ MPa$\cdot$nm. The bond energy
$\varepsilon_{0,CH}$ comes from Eq. 16, having considered an interfacial
energy between calcium hydroxide and surrounding solution $\gamma=68.4$ mJ/m2
(Estrela et al., 2005) and $n_{b}=12$ first neighbours for a bulk particle in
the FCC lattice. The interaction cutoff distance thus becomes $r_{c,CH}=0.457$
nm, from Eq. 13; this entails a failure strain
$\frac{r_{c,CH}-r_{0,CH}}{r_{0,CH}}\approx 6.5\%$, which is reasonable at the
molecular scale. With these interactions, an initially cubic grain of FCC-
arranged calcium hydroxide particles is pre-dissolved as in Masoero (2023b),
obtaining the rounded grain with diameter of ca. 7.3 nm in Fig. 5, which is
the configuration for this example.

A very large volume of solution is placed in contact with the solid grain, by
setting $\Delta V_{box}=10^{20}$ nm3. In this way, particle dissolution and
precipitation negligibly change the concentration of solvated ions, which
retain their initially assigned values 111The apparent volumes of solvated
Ca2+ and OH- ions, set to -0.025129603 nm3 and 0.001030404 nm3 respectively
(Lothenbach et al., 2019) in the $ChemDB$ database, have a negligible effect
in this example due to the large $\Delta V$ attached to the simulation box.
The chemical reactions underlying dissolution and precipitation here are:

$\displaystyle Ca(OH)_{2}\rightarrow Ca^{2+}+2OH^{-}\;\;\;$
$\displaystyle\mathrm{with}\;K_{eq,d}=6.30866\cdot 10^{-6}$ $\displaystyle
Ca^{2+}+2OH^{-}\rightarrow Ca(OH)_{2}\;\;\;$
$\displaystyle\mathrm{with}\;K_{eq,pr}=\frac{1}{K_{eq,d}}=1.58512\cdot 10^{5}$
(18)

The value of $K_{eq,d}$ is taken from Lothenbach et al. (2019). The solvated
ions are assigned the following parameters to calculate their activities as
well as the ionic strength of the solution: $z_{Ca}=+2$, $a_{Ca}^{o}=0.486$
nm, $b_{Ca}=0.15$, and $z_{OH}=-1$, $a_{OH}^{o}=1.065$ nm, $b_{OH}=0.064$
(Lothenbach et al., 2019). In this way, the activity product of the solvated
ions for Eq. 18 can be computed for any initially imposed set of
concentrations. The explored concentrations range from 0 for both Ca2+ and
$OH^{-}$ (saturation index $\beta=0$) to $c_{Ca}=39.1$ mmol/L, $c_{OH}=78.2$
mmol/L (i.e. pH $=12.89$ and $\beta=10$). Chemical equilibrium in stress-free
conditions is expected at $\beta=1$, which implies $c_{Ca}=16.43$ mmol/L and
$c_{OH}=32.86$ mmol/L (pH $=12.51$) for the fixed $c_{OH}/c_{Ca}=2$ ratio
assumed here to ensure charge neutrality of the solution.

For the dissolution reaction in Eq. 18, the activation energy $\Delta
G^{o{\ddagger}}_{d}$ and standard state concentration of the activated complex
$c^{{\ddagger}}_{d}$ are both computed from the rate constant of Ca(OH)2
dissolution in Bullard (2007): $k=0.66\;\mu$mol m-2 s-1. The conversion from
$1\;\mu$mol m-2 to number of reactions per $nm^{2}$, which are the simulation
units, yields $c^{{\ddagger}}_{d}=10^{-6}N_{a}/10^{18}=0.6022$ nm-2. Eq. 5
then provides $\Delta G^{o{\ddagger}}_{d}=74$ kJ/mol $=122.9$ ag nm2 ns-2 in
simulation units. For the precipitation reaction in Eq. 18,
$c^{{\ddagger}}_{pr}=c^{{\ddagger}}_{d}$ is assumed, thus $\Delta
G^{o{\ddagger}}_{pr}=\Delta G^{o{\ddagger}}_{d}+k_{B}T\ln(K_{eq,d})=45$ kJ/mol
$=74$ ag nm2 ns-2. The activity coefficients of the activated complex are
taken as $\gamma^{\ddagger}_{d}=\gamma^{\ddagger}_{pr}=1$, as usual for
reactions involving solids (Domínguez et al., 1998).

Deletion of existing particles and nucleation of trial particles are sampled
using two different subcommunicators: $subA$ with 1 processor, and $subB$ with
2 processors. For particle nucleation, a cubic lattice of 50,656 trial
particles is employed, with lattice spacing $r_{0,CH}/2$ in all three
directions. The trial lattice is only created in a spherical shell region with
inner and outer radii of 1.8 and 5 nm respectively, as shown in Fig. 5; this
avoids the computational cost of sampling unlikely precipitation events far
from the surface of the initial crystal, while leaving room for enough growth
or dissolution to reach a steady state regime with constant average reaction
rate. For the trial particles to move into a local minimum of interaction
energy, 600 steps of energy minimisation are performed using the
aforementioned $quickmaske$ algorithm, with time step of 0.45 ps and a maximum
displacement per step of $r_{0,CH}/200$. Such a short minimization is
appropriate because, for the interactions used here, a particle only needs to
move by $r_{0,CH}$ or less, to find its closest local energy minimum. The ions
added/removed to/from the solution after particle dissolution/precipitation
are distributed between the volume of solution $V_{sol}$ in the simulation box
and the additional $\Delta V_{box}$ in proportion to their volumes. However,
since $\Delta V_{box}\gg V_{sol}$, changes in ion concentrations and
saturation index saturation $\beta$ within the simulation box are negligible
here. After each accepted dissolution or precipitation event, a process of
type $Every$ is invoked, which minimizes the total interaction energy by
displacing all the real particles in the system; this is done via 10,000 steps
of the LAMMPS $quickmin$ algorithm, with timestep of 1 ps and a maximum
displacement per step of 0.2 nm.

###  2.6 Example 2: strain-rate effect on a crystal under stress-driven
dissolution

The nanocystal of Ca(OH)2 from the previous example is now brought into
initial contact with two platens, via a procedure detailed in Masoero (2023b):
see Fig. 6.a. The platens are arbitrarily discretized with particles of
diameter $D_{S}=0.6$ in an FCC lattice, hence with equilibrium distance
$r_{0,S}=\sqrt[3]{\eta}D_{S}=0.5427$. The mechanical interactions between
particles in the platens are of the form in Eq. 12, with elastic modulus
typical of steel, $E_{S}=200$ GPa, which yields a high interaction stiffness,
$k_{0,S}=208,397$ MPa$\cdot$nm A high interfacial energy is then assumed
between platen particles and the solution, $\gamma_{S}=200$ mJ m${}^{-}2$,
which leads to $\varepsilon_{0,S}=37.70$ MPa$\cdot$nm3. The resulting cutoff
for the platen particles is $r_{c,S}=0.5617$ nm, viz. a strain at failure of
ca. 3.5%. The mechanical interactions between platen and Ca(OH)2 particles are
of the form in Eq. 17, with equilibrium distance $r_{0,CH-S}=0.4856$ nm and
interaction stiffness $k_{CH-S}=33,200$ MPa$\cdot$nm, obtained from the
already computed, individual $r_{0}$ and stiffness of the two particle types.
The platen-Ca(OH)2 bond energy is set to $\varepsilon_{CH-S}=0$, making the
interaction purely repulsive, i.e. $r_{c,CH-S}=r_{0,CH-S}$ and no bond under
tension.

Figure 6: (a) Crystal of Ca(OH)2 in initial contact with two platens; (b)
Section of the crystal showing the distribution of per-particle virial stress
$\sigma_{z,i}$ after a small downward displacement of the top platen.

The volume of the simulation box is set to $V_{box}=8,000$ nm3, with an
additional volume in contact with it $\Delta V_{box}=10,000$ nm3. The initial
volume of the crystal is given by the initial number of Ca(OH)2 particles
times their individual volume, $V_{0,CH}=N_{0,CH}\cdot V_{p,CH}=4055\cdot
0.0557=226$ nm3, hence the initial volume of solution (in $V_{box}$ plus in
$\Delta V_{box}$) is $V_{sol,tot}\approx 18,000$ nm3. The change in
concentration of Ca2+ and OH- ions in solution caused by the dissolution of
one Ca(OH)2 can be approximately estimated by neglecting the changes in
solution volume caused by the reaction (MASKE accounts for them, through the
volume of the Ca(OH)2 particle and the apparent volumes of solvated Ca2+ and
OH- ions given in Section 2.5). For Ca2+, a dissolution event releases one
molecule, hence changing $c_{Ca}$ approximately by $1/V_{sol,tot}\approx 92$
$\mu$mol/L; for OH- ions the change is $2/V_{sol,tot}\approx 184$ $\mu$mol/L.
Therefore, in this example the dissolution of Ca(OH)2 particles during a
simulation may significantly change the concentrations of ions in solution,
and thus the saturation index $\beta$.

Only the Ca(OH)2 particles are allowed to dissolve and to nucleate, via the
same chemical reactions, thermodynamic and kinetic parameters, and trial
precipitation lattice as in Example 1. The initial ion concentrations are set
to $c_{Ca}=0.01643$ mol/L and $c_{OH}=0.03286$ mol/L, implying $\beta=1$. This
means that no sustained dissolution nor precipitation is expected in stress-
free conditions. All the particles in both platens are constrained not to move
under the action of forces; this is done using the $setforce$ fix in LAMMPS.
However, the particles in the top platen are rigidly displaced towards the
bottom platen by 0.01 nm (via the $displace\\_atoms$ command in LAMMPS) every
$dt$ nanoseconds of simulated time in MASKE. This is done by using a
continuous process of $mstoreLMP$ type in the MASKE input script. The chosen
value of $dt$ controls the downward velocity of the top platen:
$v_{z}=0.01/dt$ nm/ns. The initial distance between the centres of mass of the
two platens is 9.12 nm, hence the imposed strain rate is estimated in
$\dot{\epsilon_{z}}\approx 0.01/dt/9.12=1.1\cdot 10^{-3}/dt$ ns-1.

The platens put the Ca(OH)2 crystal under compression: see Fig. 6.b. The
average compressive stress in the crystal, $\sigma_{z}$, is estimated as the
per-particle virial stress in $z$ direction for each Ca(OH)2 particles,
$\sigma_{z,i}$ (computed in LAMMPS via a $compute$ of type $stress/atom$),
averaged over all the $N_{CH}$ particles of Ca(OH)2 existing at any time in
the simulation:

$\sigma_{z,i}=\frac{1}{V_{p,i}}\sum_{j=1}^{N_{n,i}}\left(z_{i}\cdot
F_{z,j}+z_{j}\cdot F_{z,i}\right)$ (19)

$V_{p,i}$ is the volume of the ith particle, j is the index of a neighbouring
particle (of any type), $N_{n,i}$ is the number of neighbours mechanically
interacting with particle $i$, $z$ is the Z coordinate of the two interacting
particles, and $F_{z}$ is the Z component of the force on each particle
stemming from the pairwise interaction. The presence of per-particle stresses
entails a nonzero excess enthalpy in the rate equations for dissolution and
precipitation. As a result, the compression is expected to trigger dissolution
of Ca(OH)2 even if the solution is initially at $\beta=1$. As the compressive
strain increases, dissolution will help relaxing the increasing $\sigma_{z}$,
but the concurring increase in ion concentrations in solution, and thus of
$\beta$, will try to slow down dissolution and even cause precipitation of new
Ca(OH)2 particles in stress-free regions of the crystal surface. All this will
result in a complex kinetic balance between solution chemistry and mechanical
stress, which will depend on the imposed strain rate and generate a non-
trivial evolution of the crystal morphology through dissolution/precipitation
events and plastic deformations.

##  3 Results

###  3.1 Example 1: dissolution, growth, and equilibrium of a stress-free
crystal

This section considers the Ca(OH)2 nanocrystal described in Section 2.5, with
the aim of assessing: (1) whether the simulations correctly predict the
transition from crystal dissolution to growth that, for stress-free crystals,
is expected as the solution goes from undersaturated ($\beta<1$) to
supersaturated ($\beta>1$); (2) the variability in the results stemming from
the random numbers involved in the Kinetic Monte Carlo algorithm (see B to
appreciate the role of random numbers); (3) how the fineness of the lattice of
trial particles may impact the results; (4) the pros and cons of using
straight or net rates (viz. Eqs. 3 and 4 rather than Eqs. 9 and 10).

Fig. 7 shows precipitation and dissolution rates obtained from simulations
using net rates. The rates are expressed per unit area of crystal surface $S$
and per unit time. To estimate $S$, knowing the volume $V_{p}$ of a Ca(OH)2
particle and the number $N_{CH}$ of particles at any give time in a
simulation, one can compute the total volume of the crystal $V_{p}N_{CH}$, and
set $S=\pi(6V_{p}N_{CH}/\pi)^{2/3}$, having approximated the whole crystal as
a sphere. Fig. 7.a shows a typical curve of precipitation rate (averaged over
the 100 KMC events before the current one) vs. simulated time. The time scale
extends to various seconds, despite the nanometre scale of the system: this
ability to sample long time scales irrespective of the length scale being
considered, is a strength of the KMC approach. Fig. 7.a also shows a steady
state of constant rate in the initial part of the simulation; this happens for
all $\beta\leq 1$ and $\beta\geq 8$ values, but not when $1<\beta<8$. The
steady state values are plotted against $\beta$ in Fig. 7.b; for the
$1<\beta<8$ range, the plotted rate is the one averaged over the first 100 KMC
steps of a simulation. Fig. 7.b shows that: (i) dissolution (viz. a negative
precipitation rate) is correctly predicted when $\beta<1$, as well as
precipitation when $\beta<1$ and and chemical equilibrium (viz. zero rate)
when $\beta=1$. The figure also shows that the sequence of pseudo-random
numbers used in the KMC algorithm has a negligible effect on the average rates
(‘seed 1’ and ‘seed 2’ indicate two difference sequences).

Figure 7: (a) Precipitation rate as a function of time for a Ca(OH)2 crystal
in a supersaturated solution with $\beta=10$. The rate is normalized per unit
area of crystal surface and per unit time; (b) Initial rates of dissolution
and precipitation at various saturation levels, showing the expected
transition across chemical equilibrium at $\beta=1$ for two different
sequences of pseudo-random numbers in the KMC algorithm (seed 1 and seed 2);
(c) Final crystal configuration after precipitation at $\beta=10$.

Transition State Theory (TST) predicts a linear relationship between rates per
unit surface and $\beta$; by contrast the simulations, despite using rates
that are based on TST, produce a nonlinearity in the $1<\beta<8$ range: see
Fig. 7.b. This effect is due to the crystal morphology. An overall linear
relationship, as per TST, requires that kink sites are not exhausted during
precipitation; this is possible if there is a very large number of such sites
in the initial crystal, which is not the case for the nanocrystal simulated
here. Alternatively, new kink sites must be created rapidly, which is possible
if $\beta$ is large compared to the energy barriers of precipitating particles
with fewer neighbours $n$ than a kink site, $\Delta
G_{n}=\varepsilon_{0}(n_{kink}-n)$, or if rare, energetically unfavourable
fluctuations favour precipitation in under-coordinated sites, such as adatoms,
which provide starting points for the growth of new layers. The latter
scenario is impossible when net rates are used, as they average out all the
energetically unfavourable fluctuations. As a result in Fig. 7.b, when
$1<\beta<8$, precipitation is short-lived as it only consumes the initial kink
sites. The continuous exhaustion of kink sites causes a continuous decrease in
rate, which explains why a steady state is not established in said range. When
$\beta\geq 8$, instead, new kink sites can be formed on certain crystal facets
and precipitation can proceed further and faster. However, when the layers
growing from facets on which new kinks can form at $\beta=10$ are exhausted,
precipitation stops and the crystal ends up with the clear-cut morphology in
Fig. 7.c. Even higher $\beta$ values would enable formation of new kinks on
other facets too, thus sustaining further precipitation; a systematic
exploration of these $\beta$ thresholds is not of primary interest here, but
an indication of this process at play will be given later in Example 2. A
close inspection of Fig. 7.b reveals a nonlinear increase in dissolution rate
also as $\beta$ approaches 0. This happens because a very low $\beta$ enables
the dissolution of sites that are more coordinated than kinks, creating
vacancies (or ‘pits’) on the crystal surfaces where new dissolution fronts can
emanate, accelerating the overall process.

The effect of energetically unfavourable fluctuations is shown in Fig. 8.a,b,
with simulation results from using straight rates. The $\chi$ parameter in the
rate equations controls how prevalent the effects of fluctuations are. The
explored values are $\chi=0$ (as usual in Transition State Theory), which
attributes all the excess enthalpy from mechanical interactions to the
dissolution reaction only and maximizes the fluctuations, to $\chi=0.5$, which
equally distributes the excess enthalpy between dissolution and precipitation,
minimizing fluctuations but not removing them (as opposed to net rates). Fig.
8.a shows that all the results from straight rates predict a linear
relationship between overall rate and $\beta$, also in the $1\leq\beta\leq 8$
range where net rates were unable to nucleate new crystal layers. Straight
rates instead always enable the formation and stabilization of adatoms (e.g.
those circled in Fig. 8.b) which allow for indefinite growth. A close
inspection of Fig. 8.a reveals a slight shift of equilibrium (i.e. rate = 0)
towards $\beta$ values that are slightly greater than 1; this is the effect of
relying on an estimated value of $V_{t}$ when explicitly sampling both
dissolution and precipitation events (here $V_{t}$ was set to
$V_{p,CH}\frac{r_{c,CH}}{r_{0,CH}}$).

Figure 8: (a) Precipitation rate vs. $\beta$ curves from using straight rates;
(b) a typical final morphology from simulations with straight rates and
$\chi=0.5$ (this one is for $\beta=10$). A few adatoms are circled in black,
which result from energetically unfavourable fluctuations and provide
nucleation spots for the growth of new layers; (c) Effect of the fineness of
the cubic trial lattice for precipitation, whose spacing is scaled here by a
factor $\zeta$, such that the spacing is $\zeta\cdot r_{0,CH}/2$. The kinetics
in (c) are for simulations with net rates, at $\beta=10$.

All the simulations with straight rates, $\beta>1$, and $\chi=0.5$ produce the
same morphology as in Fig. 8.b, which is reached when the entire nucleation
lattice is covered by real particles; further growth is expected if bigger
lattices were used. The same morphology is expected to emerge when $\chi=0$
and $\chi=0.05$, but in this case the simulations were stopped after a few
thousands atoms precipitated, because of the high computational cost of
sampling many unlikely fluctuations. Indeed, the precipitation of ca. 500
particles at $\beta=10$ required ca. 10,000 KMC steps in simulations with
$\chi=0$ and 8,000 steps when $\chi=0.05$, as opposed to 800 steps with
$\chi=0.5$, and 500 steps for net rates. On the other hand, Fig. 8.a shows
that the precipitation rates compute with $\chi=0$ are significantly
underestimated in all the other simulations with $\chi>0$, as well as when net
rates are used. This confirms that rare fluctuations are important for
quantitatively capturing the nucleation of new layers. Therefore, when using
MASKE one should judge on a case-by-case basis whether to use straight or net
rates and, in case, which value of $\chi$, in order to balance the
requirements of running rapid calculations, predicting qualitatively correct
mechanisms, and obtaining quantitatively correct rates.

A last parametric study here concerns the impact of the fineness of the
nucleation lattice on the resulting average rates. For simulations using net
rates, Fig. 8.c shows the temporal increase in $N_{CH}$ at $\beta=10$ when
cubic nucleation lattices of different spacing $\zeta\cdot r_{0,CH}/2$ are
employed. All simulations with lattice spacing less than $r_{0,CH}/2$, viz.
with $\zeta\leq 1$, return very similar results 222The simulation with
$\zeta=0.33$ was stopped earlier because it entailed the sampling of ca. 1.4
million trial particles, as opposed to ca. 50,000 when $\zeta=1$ and ca. 6,000
when $\zeta=2$. Coarser lattices however, with $\zeta>1$, underestimate
precipitation because they miss energetically convenient minima in the
interaction energy landscape. This indicates that a nucleation lattice spacing
of $r_{0,CH}/2$ or less is appropriate for the simulations and the mechanical
interaction potentials considered here.

###  3.2 Example 2: strain-rate effect on a crystal under stress-driven
dissolution

In this example, the Ca(OH)2 crystal from the previous section is surrounded
by a relatively small volume of solution, which is initially at chemical
equilibrium with the crystal itself, viz. $\beta=1$. The crystal is then
subjected to a compressive strain with fixed rate, which should favour the
dissolution of particles under high mechanical stress. The dissolving
particles increase the concentration of Ca2+ and OH- in solution, hence its
saturation index $\beta$ too; this may cause precipitation of new Ca(OH)2
particles in regions of the crystal surface that are under low mechanical
stress. At high imposed strain rates, plastic deformations may occur too.

Fig. 9 shows the deformation and stress field of the Ca(OH)2 crystal when
compressed at different strain rates $\dot{\epsilon}_{z}$. The KMC algorithm
in MASKE allows exploring a broad range of strain rates, here from $10^{8}$
s-1 to $10^{-4}$ s-1; this is valuable because the experiments typically
impose rates in the order of the lowest values here, which are difficult to
reach with other simulation techniques (e.g. molecular dynamics; Cao et al.
(2019)). All the simulations in Fig. 9 employed net rates; some considerations
on the effect of straight rates will be made later in this section.

At low strain rates, $\dot{\epsilon}_{z}\sim 10^{-4}$ s-1 in Fig. 9, the
dissolution triggered by the per-particle stress (via the excess enthalpy in
the rate equations) is fast compared to the deformation process. As a result,
the dissipation of local stress $\sigma_{z,i}$ through dissolution is
efficient, therefore the local stress remains relatively low even at high
$\epsilon_{z}$, as shown by the dark violet color in the $\epsilon_{z}=40$%
snapshot in Fig. 9. The same snapshot shows significant recrystallization,
driven by the increase of $\beta$ in solution eventually triggering
precipitation of new Ca(OH)2 in low-stress regions.

Higher strain rates induce higher per-particle stresses $\sigma_{z,i}$: see
$\dot{\epsilon}_{z}=1.1$ and $\sim 10^{8}$ s-1 in Fig. 9. For compatibility of
displacements, the rate at which a layer of particles in contact with the
platens dissolves must equal the velocity at which the top platen moves. A
higher strain rate entails a higher platen velocity, which requires a higher
dissolution rate, achieved by allowing a slightly higher strain before
dissolution, which in turn increases the local stress and thus the excess
enthalpy term in the rate. Fig. 9 also shows that little recrystallization
takes place at high strain rates. This is because the precipitation rate
scales linearly with $\beta$, and therefore with the number of dissolved
particles (assuming that no reprecipitation takes places at all, which is a
good approximation at high strain rates). At a given strain level
$\epsilon_{z}$, the number of dissolved particles is basically the same
irrespective of the imposed strain rate (for sufficiently high
$\dot{\epsilon}_{z}$). By contrast, the strain rates have been increased by
orders of magnitude here, as well as the dissolution rates required to
accommodate the strain; in the resulting short timescales of the
deformation/dissolution processes, precipitation events become rare. In the
simulation with highest strain rate, $\sim 10^{8}$ s-1, the animation of the
full chemo-mechanical process in Fig. 9 shows that plastic shear slips along
crystalline planes occur at high strain levels, $\epsilon_{z}\gtrsim 30\%$.
Such slips are not recorded at lower strain rates, where the dissolution-
induced relaxation of local stresses prevents the activation of plastic
deformations.

Figure 9: Snapshots from simulations of a Ca(OH)2 nanocrystal deforming,
dissolving, and reprecipitating under a compression imposed at a constant
strain rate $\dot{\epsilon}_{z}$. The color map highlights the per-particle
virial stress $\sigma_{z,i}$.

Fig. 9 also shows results for a purely mechanical response, i.e. with both
dissolution and precipitation turned off. The simulation was conducted at a
low strain rate $\dot{\epsilon}_{z}\sim 10^{-4}$ s-1, but the same result
would be obtained at higher strain rates because the particles are brought to
static equilibrium between two successive execution of the platen-moving
continuous process in MASKE (static equilibrium is imposed via a MASKE process
of type $Every$ performing energy minimization in MASKE). Therefore inertial
effects do not contribute to rate-dependence here. The purely mechanical
response in Fig. 9 is somewhat similar to the coupled chemo-mechanical one at
high $\dot{\epsilon}_{z}\sim 10^{8}$ s-1, but with some important differences
too. At low strain levels (see $\epsilon_{z}=10\%$ in Fig. 9), even the
limited dissolution taking place at $\dot{\epsilon}_{z}\sim 10^{8}$ s-1 is
sufficient to reduce the local stresses compared to the purely mechanical
test. At $\epsilon_{z}\sim 20\%$ and $40\%$ the purely mechanical test
displays a wider low-stress region near the surface of the crystal; the
particles in this region do not contribute significantly to the overall
mechanical response to the compression. Partial dissolution at
$\dot{\epsilon}_{z}\sim 10^{8}$ instead moulds the crystal layers in contact
with the platens, removing particles with high $\sigma_{z,i}$ that induce
stress concentration, and thus enabling a more uniform mechanical response.
Another distinctive aspect of the purely mechanical test is that extensive
plastic deformations from shear slips already start at $\epsilon_{z}\approx
17\%$, and a fully developed central shear band controls the entire
deformation response at $\epsilon_{z}\gtrsim 30\%$. In the dissolving crystal
with $\dot{\epsilon}_{z}\sim 10^{8}$, instead, small and isolated, plastic
shear slips only start appearing at $\epsilon_{z}\approx 30\%$.

Fig. 10.a shows the stress strain curves computed at different strain rates
$\dot{\epsilon}_{z}$. A higher $\dot{\epsilon}_{z}$ produces a higher average
compressive stress $\sigma_{z}$; this reflects the already-discussed need for
faster dissolution (whose rate depends on local stresses) to ensure
compatibility of displacements at the platen-crystal interface.

Figure 10: (a) Stress-strain curves for the Ca(OH)2 nanocrystal from
simulations at different strain rates $\dot{\epsilon}_{z}$; (b) Evolution of
saturation index $\beta$ during the simulations; (c) Evolution of the
precipitation rate (negative for dissolution), averaged over 100 KMC steps,
from the simulations with $\dot{\epsilon}_{z}\approx 10^{-4}$ and $1$ s-1.

A counter-intuitive result in Fig. 10.a is that crystals that are allowed to
dissolve, when compressed at high rates, are stronger than crystals that are
not allowed dissolve. This is indicated in Fig. 10.a, by the higher peak
$\sigma_{z}$ of the curve for $\dot{\epsilon}_{z}\approx 10^{8}$ s-1, compared
to the ‘mech. only’ curve. The ‘mech. only’ curve lacks the mechanism of
particle dissolution, which helps relaxing the local stresses on particles
while ensuring compatibility of displacement at the crystal-platen interface.
Such compatibility, in the ‘mech. only’ simulation, relies entirely on
mechanical deformations; these are initially elastic and, when local stresses
become too high, plastic shear slips are triggered and relax the local stress.
One may thus expect that the ‘mech. only’ simulation, being more constrained,
would produce a higher $\sigma_{z}$ than simulations where dissolution is
possible too. This is indeed the case at low strain levels: see
$\epsilon_{z}<0.2$ in Fig. 10.a. However, the first extensive shear slips in
the ‘mech. only’ simulation occur already at $\epsilon_{z}\approx 17\%$. After
that point, subsequent shear slips effectively cap the peak stress
$\sigma_{z}$ to ca. 2’000 MPa. By contrast, in the $\dot{\epsilon}_{z}\approx
10^{8}$ s-1 simulation, dissolution of the most stressed particles limits the
maximum local stress, which gets redistributed as a lower stress to the
particles neighbouring the dissolved one. All this leads to a more uniform
stress distribution in the crystal (see Fig. 9) while also delaying the
triggering of plastic slips until $\epsilon_{z}\approx 30\%$. As a result, the
peak stress $\sigma_{z}$ in the $\dot{\epsilon}_{z}\approx 10^{8}$ s-1 is ca.
2’500 MPa, greater than in the ‘mech. only’ case.

Another finding from Fig. 10.a is the stepwise stress-strain curve at the
lowest strain rate, $\dot{\epsilon}_{z}\approx 10^{-4}$ s-1. Specifically,
$\sigma_{z}$ increases with $\epsilon_{z}$ for a while, then it oscillates
around a constant value while $\epsilon_{z}$ keeps increasing, until a new
regime of increasing $\sigma_{z}$ is entered, and so forth. Three such steps
can be counted in Fig. 10.a. This behaviour is the result of the kinetic
chemo-mechanical balance established during the simulations. At low rates
$\dot{\epsilon}_{z}$, no plastic deformations are triggered and particle
dissolution is the only mechanisms to accommodate the externally imposed
$\dot{\epsilon}_{z}$. When compression starts, the first dissolution events
release ions in solution, raising the saturation index $\beta$ for Ca(OH)2
precipitation: see the increasing $\beta$ at $\epsilon_{z}<0.15$ in the
$\dot{\epsilon}_{z}\approx 10^{-4}$ s-1 curve in Fig. 10.b. With the net rate
equations used here, Example 1 in the previous section has shown that no
significant precipitation is possible when $\beta\lesssim 8$, hence $\beta$
keeps increasing. The increasing $\beta$ slows down dissolution but, since the
overall rate $\dot{\epsilon}_{z}$ is imposed and must be matched, a
compensating increase in dissolution rate is made possible through an increase
in mechanical stress (which increases the excess enthalpy term in the
dissolution rate equation). This explains the regimes of increasing
$\sigma_{z}$ in Fig. 10.a, which indeed correspond to regimes of increasing
$\beta$ in Fig. 10.b. At $\beta\approx 8$, particle precipitation consumes any
further ion that dissolution releases in solution, so that the average rate
(precipitation minus dissolution over 100 KMC steps) oscillates around zero:
see the curve for $\dot{\epsilon}_{z}\approx 10^{-4}$ s-1 in Fig. 10.c when
$0.15<\epsilon_{z}<0.3$. As a result, $\beta$ now remains constant while
$\epsilon_{z}$ increases, and the constant overall rate $\dot{\epsilon}_{z}$
implies that the stress $\sigma_{z}$ must remain constant too. This explains
the matching plateaus in Figs. 10.a and 10.b. The first plateau continues
until exhausting all the surface sites where precipitation at $\beta\approx 8$
is possible: cf. the final morphology in Fig. 7.b. Then a new regime of
increasing $\beta$ and $\sigma_{z}$ is entered, until a new plateau occurs at
$\beta\approx 50$, when precipitation on other exposed crystal surfaces
becomes possible too. The new plateau again corresponds to a net dissolution
rate oscillating around zero, in Fig. 10.c. All this describes a complex,
kinetic equilibrium between composition of the solution, stress state in the
crystal, and microstructural evolution.

The stepwise shape of the stress-strain curve is lost at higher strain rates,
as shown by the $\dot{\epsilon}_{z}\approx 1$ s-1 curves in Fig. 10.
Additional results omitted here show a progressive reduction and disappearance
of the plateaus for intermediate rate values between
$\dot{\epsilon}_{z}\approx 10^{-4}$ s-1 and 1 s-1. This means that, despite
some particles do still precipitate when $\beta\gtrsim 8$, they cannot compete
any more with the faster rate of the dissolving particles. Consistently, when
$\dot{\epsilon}_{z}\approx 1$ s-1 the saturation index $\beta$ increases
monotonically in Fig. 10.b and the average dissolution rate never goes to zero
in Fig. 10.c. The $\beta(\epsilon_{z})$ curves in Fig. 10.b for
$\dot{\epsilon}_{z}\approx 1$ and $10^{+4}$ s-1 are similar, meaning that
dissolution progresses at a similar rate. This is expected because, in a
scenario where precipitation is negligible, the total amount of dissolved ions
depends only on the accommodated strain $\epsilon_{z}$. If $\beta$ is
approximately identical for $\dot{\epsilon}_{z}\approx 1$ and $10^{+4}$ s-1,
then $\sigma_{z}$ is necessarily greater at the higher imposed rate, since
only $\sigma_{z}$ now can speed up the dissolution rate to match
$\dot{\epsilon}_{z}$. Less intuitively, the $\beta(\epsilon_{z})$ curve for
$\dot{\epsilon}_{z}\approx 10^{8}$ s-1 in Fig. 10.b is initially similar to
those for $\dot{\epsilon}_{z}\approx 1$ and $10^{+4}$ s-1, but it
significantly diverges from them at $\epsilon_{z}\approx 0.2-0.3$, leading to
less concentrated solutions at high strain levels $\epsilon_{z}$. The reason
for this lies in the triggering of plastic deformations at high rate,
$\dot{\epsilon}_{z}\approx 10^{8}$ s-1, which cap the stress level (see
$\sigma_{z}$ at $\epsilon_{z}\gtrsim 0.3$ in Fig. 10.a) and provide another
mechanism to accommodate the imposed strain. With a capped $\sigma_{z}$, the
increasing $\beta$ reduces the dissolution rate, which in turn reduces the
slope of the $\beta(\epsilon_{z})$ curve in Fig. 10.b. The resulting
interdependence between plastic deformations and composition of the solution
is another non-trivial result of the chemo-mechanical coupling in MASKE.

All the results in this section are from simulations where the $allpar$
mechanism in MASKE employed net rates of dissolution and precipitation, as per
Eqs. 9 and 10. Additional tests were conducted using straight rates too, as
per Eqs. 3 and 4, with $\chi$ values between 0 and 0.5. At high strain strain
rates, $\dot{\epsilon}_{z}\gtrsim 1$ s-1, straight and net rates produce
almost identical results; this is expected, since the main effect of using
straight straight is to favour precipitation, which plays a negligible role at
high $\dot{\epsilon}_{z}$. Instead, at low $\dot{\epsilon}_{z}\lesssim 1$,
simulations using straight rates produced lower $\beta(\epsilon_{z})$ than
those from net rates at same $\dot{\epsilon}_{z}$. This agrees with the
findings from Example 1 in the previous section, where straight rates were
shown to enable precipitation at lower $\beta$ values. Therefore, simulations
with straight rates are indeed expected to produce plateaus in $\beta$ and
$\sigma_{z}$ at lower $\epsilon_{z}$ than simulations with net rates. These
aspect, however, should be studied further because the efficiency loss from
using straight rates significantly limited the maximum $\epsilon_{z}$ that
could be reached during the simulations. For the purpose of this paper, the
results presented here for net rates have already demonstrated how MASKE can
predict the emergence of a complex kinetic balance between chemical and
mechanical processes over long time scales.

##  4 Conclusion and outlook

The article has introduced MASKE, a kinetic simulator of microstructural
evolution driven by chemical reactions and mechanical stress. Solid domains
are modelled as mechanically interacting particles, while the surrounding
environment, which was an aqueous solution in the examples here, is modelled
implicitly through the concentrations of molecular species in it. The
composition of the solution co-evolves with the solid phases, in respect of
mass conservation, all while dissolution and precipitation reactions transform
molecules from the solid into solvated species and vice versa.

The uniqueness of MASKE stems from several key strengths:

  1. 1.

MASKE implements an off-lattice, rejection-free, Kinetic Monte Carlo (KMC)
framework to simulate the dissolution and precipitation of particles. This
framework gives access to realistic time scales that are relevant for the
experiments;

  2. 2.

MASKE features original chemo-mechanical rate equations for dissolution and
precipitation, which capture the effect of mechanical stress on the reaction
rates. To simulate the deformations that induce mechanical stress, the
particles are free to move off-lattice, but this implies that infinite
positions are to be sampled for possible particle precipitation. To treat this
issue, a discretization method has been proposed which approximates the
integral of the precipitation rate over the whole simulation box by
considering only a finite number of trial particles;

  3. 3.

The reaction rates are implemented in both straight and net forms. Depending
on the application, the user can opt for the computational efficiency of net
rates, or for more computationally demanding but also more precise
(mechanistically and quantitatively) simulations using straight rates;

  4. 4.

The KMC sampling of discrete dissolution/precipitation events is coupled with
the explicit integration of continuous processes. In the article, this
coupling has been exploited to impose a strain rate. Other continuous
processes are currently being implemented, in particular one that simulates
the metabolic activity of bacteria, explicitly modelled as particles as in Li
et al. (2019), for the purpose of modelling biomineralization in self-healing
concrete Bagga et al. (2022);

  5. 5.

In MASKE, the user can define bespoke chemical reactions and mechanisms in a
chemistry database file. MASKE is also interfaced with the LAMMPS library,
through which the user can specify complex initial microstructures, loading
conditions, and interaction potentials. As a result, a wide variety of
problems with disparate chemistries can be simulated;

  6. 6.

MASKE features a two level parallelisation, ready for use in High Performance
Computing. The simulations in this article used both levels of
parallelisation, but on simple examples that required only 3 processors in
total (despite one simulation already featured 1.4 million trial particles).
Applications of MASKE involving more processors have been presented elsewhere
(Coopamootoo and Masoero, 2020; Alex et al., 2023).

A number of additional features provide promising avenues for further
development. These in particular may involve:

  1. 1.

Coarse-graining the reaction rates to reach microstructural scales. The
$allpar$ mechanism presented here is only accurate for unimolecular reactions,
but Shvab et al. (2017) already proposed coarse-grained rate expression for
nanoparticles with diameter of ca. 10 nm. The version of MASKE that is
currently on GitHub already includes a mechanism for particles with diameter
of $\sim 1$ $\mu$m, which will be presented in a separate article;

  2. 2.

Allowing particles to partially grow or dissolve, whereas now they can only
appear or disappear at once. This would be desirable when using coarser
particles, to avoid unrealistically large, sudden changes in stress (and thus
in stress-dependent reaction rates) when a full particle is deleted or added.
Partial dissolution and growth would thus be important for capturing well
certain phenomena, such as crystallization pressure (Masoero, 2023b). Even
better if the partial growth or dissolution of particles were anisotropic.
Such anisotropy would at least require ellipsoidal particles, which would be
rather straightforward to include in MASKE, since ellipsoidal particles and
appropriate potentials for them already exist in LAMMPS (Berardi et al.,
1998). However, partial growth and dissolution (either as discrete KMC events
or as continuous processes) will require additional implementation and a new,
probably less efficient, way of sampling energy changes $\Delta U$ to compute
the excess enthalpy terms in the reaction rates;

  3. 3.

Improving the model of the solution, allowing for heterogenous local
concentrations of solvated molecules, which may also react chemically within
the implicit fluid. For the reactions in solution, there is ongoing work to
couple MASKE with the PHREEQC library (Parkhurst et al., 1999). Heterogeneous
concentrations would be desirable for simulations at larger length scales, and
the resulting concentration gradients would prompt diffusion that can be
implemented in MASKE as a continuous process to be integrated in time.

Two examples in this article have shown that MASKE, already in its current
form, can capture interesting coupled phenomena that challenge other
simulation methods. Example 1 confirmed that MASKE correctly predicts the
dissolution, growthm and chemical equilibrium of a nanocrystal of Ca(OH)2,
depending on the concentrations of Ca2+ and OH- ions in the surrounding
solution. Straight rate equations have been shown to effectively allow the
nucleation of new crystal layers through energetically unfavourable
fluctuations. Example 2 then addressed the effects of chemo-mechanical
coupling, by imposing a compressive strain rate to the crystal. The results
showed that a complex kinetic balance is established between dissolution
rates, chemical composition of the solution surrounding the crystal, elastic
stress state in the crystal, and plastic deformations too. A family of curves
is obtained, which quantifies the stress-strain response of the crystal as a
function of the imposed strain rate. These curves offer a constitutive models
for simulations al larger scales, e.g. continuum simulations based on the
Finite Element Method. All in all, MASKE offers a new way to reconstruct the
chemo-mechanical kinetics of evolving microstructures. This capability can be
leveraged to simulate the formation and behaviour of new materials in
computer-aided design, and to predict the degradation of materials for both
conventional and accidental scenarios of future exposure.

## References

  * Alex et al. (2023) Alex, A., Freeman, B., Jefferson, A., Masoero, E., 2023\.  Carbonation and self-healing in concrete: Kinetic monte carlo simulations of mineralization.  Cement and Concrete Composites , 105281\.
  * Bagga et al. (2022) Bagga, M., Hamley-Bennett, C., Alex, A., Freeman, B.L., Justo-Reinoso, I., Mihai, I.C., Gebhard, S., Paine, K., Jefferson, A.D., Masoero, E., et al., 2022\.  Advancements in bacteria based self-healing concrete and the promise of modelling.  Construction and Building Materials 358, 129412.
  * Bažant et al. (1997) Bažant, Z.P., Hauggaard, A.B., Baweja, S., Ulm, F.J., 1997\.  Microprestress-solidification theory for concrete creep. I: Aging and drying effects.  Journal of Engineering Mechanics 123, 1188–1194.
  * Bellomo et al. (2012) Bellomo, F.J., Armero, F., Nallim, L.G., Oller, S., 2012\.  A constitutive model for tissue adaptation: necrosis and stress driven growth.  Mechanics Research Communications 42, 51–59.
  * Bentz et al. (1994) Bentz, D.P., Coveney, P.V., Garboczi, E.J., Kleyn, M.F., Stutzman, P.E., 1994.  Cellular automaton simulations of cement hydration and microstructure development.  Modelling and Simulation in Materials Science and Engineering 2, 783.
  * Berardi et al. (1998) Berardi, R., Fava, C., Zannoni, C., 1998.  A gay–berne potential for dissimilar biaxial particles.  Chemical physics letters 297, 8–14.
  * Bishnoi and Scrivener (2009) Bishnoi, S., Scrivener, K.L., 2009\.  $\mu$ic: A new platform for modelling the hydration of cements.  Cement and Concrete Research 39, 266–274.
  * Bullard (2007) Bullard, J.W., 2007.  A three-dimensional microstructural model of reactions and transport in aqueous mineral systems.  Modelling and Simulation in Materials Science and Engineering 15, 711.
  * Bullard et al. (2010) Bullard, J.W., Enjolras, E., George, W.L., Satterfield, S.G., Terrill, J.E., 2010.  A parallel reaction-transport model applied to cement hydration and microstructure development.  Modelling and Simulation in Materials Science and Engineering 18, 025007.
  * Bullard et al. (2011) Bullard, J.W., Lothenbach, B., Stutzman, P.E., Snyder, K.A., 2011\.  Coupling thermodynamics and digital image models to simulate hydration and microstructure development of portland cement pastes.  Journal of Materials Research 26, 609–622.
  * Cao et al. (2019) Cao, P., Short, M.P., Yip, S., 2019.  Potential energy landscape activations governing plastic flows in glass rheology.  Proceedings of the National Academy of Sciences 116, 18790–18797.
  * Cefis and Comi (2017) Cefis, N., Comi, C., 2017.  Chemo-mechanical modelling of the external sulfate attack in concrete.  Cement and Concrete Research 93, 57–70.
  * Coopamootoo and Masoero (2020) Coopamootoo, K., Masoero, E., 2020\.  Simulations of crystal dissolution using interacting particles: prediction of stress evolution and rates at defects and application to tricalcium silicate.  The Journal of Physical Chemistry C 124, 19603–19615.
  * Coopamootoo K (2023) Coopamootoo K, M.E., 2023.  Simulations of tricalcium silicate dissolution at screw dislocations: effects of fnite crystal size and mechanical interaction potentials.  Cement and Concrete Research In press.
  * Domínguez et al. (1998) Domínguez, A., Jimenez, R., López-Cornejo, P., Pérez, P., Sánchez, F., 1998.  On the calculation of transition state activity coefficient and solvent effects on chemical reactions.  Collection of Czechoslovak chemical communications 63, 1969–1976.
  * Dubrovskii et al. (2010) Dubrovskii, V., Sibirev, N., Zhang, X., Suris, R., 2010\.  Stress-driven nucleation of three-dimensional crystal islands: from quantum dots to nanoneedles.  Crystal growth & design 10, 3949–3955.
  * Estrela et al. (2005) Estrela, C., Estrela, C.R.d.A., Guimarães, L.F., Silva, R.S., Pécora, J.D., 2005.  Surface tension of calcium hydroxide associated with different substances.  Journal of Applied Oral Science 13, 152–156.
  * Feng et al. (2017) Feng, P., Bullard, J.W., Garboczi, E.J., Miao, C., 2017\.  A multiscale microstructure model of cement paste sulfate attack by crystallization pressure.  Modelling and simulation in materials science and engineering 25, 065013.
  * Freeman et al. (2019) Freeman, B.L., Cleall, P.J., Jefferson, A.D., 2019.  An indicator-based problem reduction scheme for coupled reactive transport models.  International Journal for Numerical Methods in Engineering 120, 1428–1455.
  * Gawin et al. (2009) Gawin, D., Pesavento, F., Schrefler, B.A., 2009.  Modeling deterioration of cementitious materials exposed to calcium leaching in non-isothermal conditions.  Computer Methods in Applied Mechanics and Engineering 198, 3051–3083.
  * Hart et al. (2017) Hart, N.H., Nimphius, S., Rantalainen, T., Ireland, A., Siafarikas, A., Newton, R., 2017\.  Mechanical basis of bone strength: influence of bone material, bone structure and muscle action.  Journal of musculoskeletal & neuronal interactions 17, 114.
  * Ioannidou et al. (2016) Ioannidou, K., Krakowiak, K.J., Bauchy, M., Hoover, C.G., Masoero, E., Yip, S., Ulm, F.J., Levitz, P., Pellenq, R.J.M., Del Gado, E., 2016\.  Mesoscale texture of cement hydrates.  Proceedings of the National Academy of Sciences 113, 2029–2034.
  * Ioannidou et al. (2022) Ioannidou, K., Labbez, C., Masoero, E., 2022.  A review of coarse grained and mesoscale simulations of c–s–h.  Cement and Concrete Research 159, 106857.
  * Ioannidou et al. (2014) Ioannidou, K., Pellenq, R.J.M., Del Gado, E., 2014.  Controlling local packing and growth in calcium–silicate–hydrate gels.  Soft Matter , 1121–1133.
  * Kouhen et al. (2023) Kouhen, M., Dimitrova, A., Scippa, G.S., Trupiano, D., 2023\.  The course of mechanical stress: Types, perception, and plant response.  Biology 12, 217\.
  * Langmuir (1997) Langmuir, D., 1997.  Aqueous environmental.  Geochemistry Prentice Hall: Upper Saddle River, NJ 600\.
  * Lasaga (2014) Lasaga, A.C., 2014.  Kinetic theory in the earth sciences.  Princeton University Press.
  * Li et al. (2019) Li, B., Taniguchi, D., Gedara, J.P., Gogulancea, V., Gonzalez-Cabaleiro, R., Chen, J., McGough, A.S., Ofiteru, I.D., Curtis, T.P., Zuliani, P., 2019\.  Nufeb: A massively parallel simulator for individual-based modelling of microbial communities.  PLoS Computational Biology 15, e1007125.
  * Li et al. (2015) Li, X., Grasley, Z., Garboczi, E.J., Bullard, J.W., 2015\.  Modeling the apparent and intrinsic viscoelastic relaxation of hydrating cement paste.  Cement and Concrete Composites 55, 322–330.
  * Lothenbach et al. (2019) Lothenbach, B., Kulik, D.A., Matschei, T., Balonis, M., Baquerizo, L., Dilnesa, B., Miron, G.D., Myers, R.J., 2019\.  Cemdata18: A chemical thermodynamic database for hydrated portland cements and alkali-activated materials.  Cement and Concrete Research 115, 472–506.
  * Martin et al. (2020) Martin, P., Gaitero, J.J., Dolado, J.S., Manzano, H., 2020\.  Kimera: a kinetic montecarlo code for mineral dissolution.  Minerals 10, 825\.
  * Masoero (2023a) Masoero, E., 2023a.  MASKE GitHub Repository.
  * Masoero (2023b) Masoero, E., 2023b.  Maske: Particle-based chemo-mechanical simulations of degradation processes, in: International RILEM Conference on Synergising expertise towards sustainability and robustness of CBMs and concrete structures, Springer. pp. 159–170.
  * Meakin and Rosso (2008) Meakin, P., Rosso, K.M., 2008\.  Simple kinetic monte carlo models for dissolution pitting induced by crystal defects.  The Journal of chemical physics 129\.
  * Parkhurst (1995) Parkhurst, D.L., 1995.  User’s guide to PHREEQC: A computer program for speciation, reaction-path, advective-transport, and inverse geochemical calculations.  95-4227, US Department of the Interior, US Geological Survey.
  * Parkhurst et al. (1999) Parkhurst, D.L., Appelo, C., et al., 1999\.  User’s guide to phreeqc (version 2): A computer program for speciation, batch-reaction, one-dimensional transport, and inverse geochemical calculations.  Water-resources investigations report 99, 312.
  * Parrott (1975) Parrott, L., 1975.  Increase in creep of hardened cement paste due to carbonation under load.  Magazine of Concrete Research 27, 179–181.
  * Pesavento et al. (2012) Pesavento, F., Gawin, D., Wyrzykowski, M., Schrefler, B.A., Simoni, L., 2012.  Modeling alkali–silica reaction in non-isothermal, partially saturated cement based materials.  Computer Methods in Applied Mechanics and Engineering 225, 95–115.
  * Prados et al. (1997) Prados, A., Brey, J., Sánchez-Rey, B., 1997.  A dynamical monte carlo algorithm for master equations with time-dependent transition rates.  Journal of statistical physics 89, 709–734.
  * Raja and Shoji (2011) Raja, V., Shoji, T., 2011.  Stress corrosion cracking: theory and practice.  Elsevier.
  * Rutter (1983) Rutter, E., 1983.  Pressure solution in nature, theory and experiment.  Journal of the Geological Society 140, 725–740.
  * Scherer (2004) Scherer, G.W., 2004.  Stress from crystallization of salt.  Cement and concrete research 34, 1613–1624.
  * Shvab et al. (2017) Shvab, I., Brochard, L., Manzano, H., Masoero, E., 2017\.  Precipitation mechanisms of mesoporous nanoparticle aggregates: off-lattice, coarse-grained, kinetic simulations.  Crystal Growth & Design 17, 1316–1327.
  * Thompson et al. (2022) Thompson, A.P., Aktulga, H.M., Berger, R., Bolintineanu, D.S., Brown, W.M., Crozier, P.S., in’t Veld, P.J., Kohlmeyer, A., Moore, S.G., Nguyen, T.D., et al., 2022\.  Lammps-a flexible simulation tool for particle-based materials modeling at the atomic, meso, and continuum scales.  Computer Physics Communications 271, 108171.
  * Tourret et al. (2022) Tourret, D., Liu, H., LLorca, J., 2022.  Phase-field modeling of microstructure evolution: Recent applications, perspectives and challenges.  Progress in Materials Science 123, 100810.
  * Truesdell and Jones (1974) Truesdell, A.H., Jones, B.F., 1974\.  Wateq, a computer program for calculating chemical equilibria of natural waters.  J. Res. US Geol. Surv 2, 233–248.
  * Van Breugel (1995) Van Breugel, K., 1995.  Numerical simulation of hydration and microstructural development in hardening cement-based materials (i) theory.  Cement and concrete research 25, 319–331.
  * Wittmann (1986) Wittmann, F.H., 1986.  Estimation of the modulus of elasticity of calcium hydroxide.  Cement and Concrete Research 16, 971–972.

##  Appendix A Details on system description

The main article has presented a typical simulation in MASKE as featuring
explicit particles, typically representing solid units, and an implicit fluid
surrounding them: see Fig. 1. The current version of MASKE assumes prismatic
boxes (orthogonal in LAMMPS syntax) but extensions to other box shapes, e.g.
triclinic ones, are within scope for future implementation. MASKE accepts all
boundary conditions in LAMMPS, but for now only fixed and periodic conditions
have been tested. LAMMPS assigns a type to each particle, which MASKE uses to
associate possible events to sets of particles: for example, MASKE may
associate a certain reaction formula to the dissolution of type 1 particles,
and a different formula to particles of type 2. One would often associate
different types to different solid species, but MASKE does not require that,
as it does not track the chemical composition nor the type of individual
particles; the user is responsible for attributing events and reactions
consistently with particle types. Such implementation is appropriate when a
particle represents a homogeneous, mono-phase solid, such as a nanocrystal;
future versions of MASKE may include tracking of the internal structure of
particles, e.g. to simulate partial growth or dissolution of particles each
representing a polycrystal.

Particle shapes are also defined in LAMMPS, via the atom_style command. For
example, the sphere style includes information on per-particle radius, whereas
ellipsoid includes per-particle shape through radii along three principal
axes. MASKE for now assumes spherical particles, hence only uses per-particle
radii. Future implementation may also include other atom_styles from LAMMPS,
e.g. to consider individual atoms and directly apply MASKE at the atomistic
level, or to employ ellipsoidal particles for simulating richer coarse-grained
morphologies.

The particles interact via effective potentials, which are defined and managed
in LAMMPS: e.g. see Fig. 1. The particle displacements caused by the
interactions are also computed by LAMMPS, with typical methods such as energy
minimization or molecular dynamics. Currently MASKE extracts per-particle
interaction energies from LAMMPS, obtained through a pe/atom compute. MASKE
uses these energies to quantify the change in total energy of the system,
$\Delta U_{i}$, generated by adding or deleting the generic particle $i$. This
$\Delta U_{i}$ informs the reaction rates. For pairwise interactions, $\Delta
U_{i}$ is simply twice the per-particle energy of particle $i$ read from
LAMMPS. Currently MASKE assumes that this relationship holds, hence that
pairwise interactions are used in LAMMPS (these include user-defined,
tabulated, pair potentials). Higher-order interactions, such as three-body or
else, will need future implementation and a different, probably more time-
consuming, approach to compute $\Delta U_{i}$ in MASKE.

The implicit solution surrounding the solid particles is defined and managed
directly in MASKE, which tracks the number of molecules and molar
concentrations for various, user-defined, solvated species. Molarities require
the user to specify the molecular volume of each solvated species, so that the
total volume of the solution can be tracked while chemical reactions alter it.
Future versions of MASKE may adopt molalities instead. The concentrations of
solvated species are assumed to be uniform in the simulation box, which is the
same as assuming infinitely fast diffusion. This hypothesis is realistic at
small scales, e.g. below the micrometre. Diffusion algorithms may be included
in future implementations. Currently, MASKE does not minimize the free energy
of the solution through speciation, neither it considers the kinetics of
chemical reactions between molecules in solution, unless these reactions
produce or dissolve solid particles as a result. There is however no obstacle
to implementing these features and, indeed, current development is underway to
couple MASKE with PHREEQC to include speciation in solution.

During a typical simulation, molecules disappear from solution to form new
solid particles, or vice versa. During these transformations, MASKE preserves
the total number of molecules, as long as the user has consistently associated
particle types to solid species, and to chemical reactions with correct
stoichiometry. The volume of molecules in solution typically differs from that
of the same molecules when combined in a solid; as a result, the total volume
of solid particles plus implicit solution may change as the system evolves.
The MASKE variable $voidV$ tracks the difference between the volume of the
simulation box, $V_{box}$, and the volume of solid plus solution in it. If
$voidV<0$, then the total volume of the system has expanded beyond $V_{box}$.
Presently, $voidV$ is just recorded for information, but in the future it may
be used to change the volume of the box or to estimate the pressure of the
fluid.

A last element of the model is an additional volume of solution, $\Delta V$,
implicitly attached to the simulation box. $\Delta V$ initially features the
same molecules with same concentrations as in the solution in $V_{box}$.
Currently there is no process in MASKE to exchange molecules between $V_{box}$
and $\Delta V$, but dissolution/precipitation events taking place in $V_{box}$
can use molecules from $\Delta V$ if allowed by the user.

##  Appendix B Processes and time line in MASKE

MASKE allows for both discrete and continuous processes to take place
contextually. An example of discrete process may be the precipitation of a new
solid particle, which occurs instantaneously after a time $\Delta t$ from a
previous event, e.g. precipitation or dissolution of another particle. A
continuous process, instead, advances gradually with time; examples may be the
diffusion of molecules in solution, displacements being applied to certain
particles with a given rate, or solid particles growing or dissolving only in
part. Continuous processes are typically governed by differential equations
that MASKE integrates over time, with a user-specified time step $dt_{i}$ for
the $i^{th}$ process. The user also provides a number of integration steps
$n_{i}$, so that whenever the $i^{th}$ continuous process is invoked, MASKE
integrates it over a time period $\Delta t_{i}=dt_{i}\cdot n_{i}$ (future
implementations may include adaptive $dt_{i}$ and criteria to stop integrating
processes that have reached equilibrium).

To stagger discrete and continuous processes, MASKE employs the concept of
future time of occurrence, $FOT$. For each $i^{th}$ continuous process,
$FOT_{i}=t_{i}+\Delta t_{i}$, where $t_{i}$ is the last time that the process
has reached through previous integration: see processes 1 and 2 in Fig. 11.a.

Figure 11: Time line in MASKE with continuous processes 1 and 2, and with the
discrete KMC processes grouped under the letter $k$. The subfigures show three
successive steps in MASKE, where the process to be carried out (in red) is the
one with nearest future time of occurrence, FOT. (a) carries out continuous
process 2, (b) carries out continuous process 1, and (c) carries out a
discrete event. In (a) the current time $t^{*}$ was set by a discrete event
that had just been carried out, hence the cumulative rate of discrete events
$R_{k,j}$ was the first one to be calculated since accepting a discrete event,
i.e. $n_{cp}=0$ in Eq. 21. Subfigure (c) eventually brings back a scenario
where a discrete event is accepted, hence $n_{cp}$ is reset after that and the
following step will be analogous to the one in (a). Note how $\Delta t_{i}$
and $FOT_{i}$ of continuous events remain constant (in the current
implementation of MASKE), whereas $\Delta t_{k}$ and $FOT_{k}$ generally
change whenever $t^{*}$ advances.

The $FOT_{k}$ of discrete events, instead, is determined by a single $\Delta
t_{k}$ that accounts for all discrete events, in the spirit of rejection-free
Kinetic Monte Carlo (KMC). Consider, for example, a simulation box containing
$N$ solid particles and two, user-defined, possible discrete processes:
dissolution of any of said $N$ particles, and nucleation of one new particle
at $M$ possible locations; this means $N+M$ possible discrete events. MASKE
computes a rate for each of these events, and cumulates them into a total KMC
rate $R_{k}$. The time increment associated to one discrete event taking place
is thus:

$\Delta t_{k}=R_{k}^{-1}\ln\left(\frac{1}{u}\right)$ (20)

where $u$ is a random number extracted with uniform probability in (0,1]. Eq.
20 is only correct if the rates of all the discrete events remain constant
between two successive realizations of a discrete events; generally this is
not the case in MASKE, since continuous processes taking place meanwhile might
change the composition of the solution or the stress state, both of which can
impact the rate of individual events. To account for this, a more general
expression for $\Delta t_{k}$ will be given shortly below.

Fig. 11 depicts how $FOT$s are used to advance a simulation in MASKE. Without
losing generality, assume that just before Fig. 11.a, a discrete event had
taken place, determining the value of the current time in MASKE to be
$t^{*}=t_{k}$. MASKE computes the $\Delta t_{i}$ and $FOT_{i}$ of all
continuous processes ($i=1,2$ in the figure, with last occurrence times
$t_{1}$ and $t_{2}$ respectively), and the $\Delta t_{k}$ and $FOT_{k}$ of a
generic discrete event, using Eq. 20. In Fig. 11.a, the continuous process 2
has the smallest $FOT$ and is therefore carried out, setting the new $t^{*}$
and $t_{2}$ in Fig. 11.b to coincide with $FOT_{2}$ from Fig. 11.a. MASKE
computes the new $FOT_{2}$ in Fig. 11.b, but it must also update $FOT_{k}$ to
account for the change in rates that the execution of continuous process 2
might have caused. To this end MASKE employs a time-dependent formulation of
KMC, whereby the remaining $\Delta t_{k}$ follows the relation:

$\sum_{j=1}^{n_{cp}}{\Delta t^{*}_{j}R_{k,j}}+\Delta
t_{k}R_{k,n_{cp}+1}=\ln\left(\frac{1}{u}\right)$ (21)

$n_{cp}$ is the number of time increments ($\Delta t^{*}_{j}$) that continuous
processes have brought about since the last discrete event. For the scenario
in Fig. 11.b, $n_{cp}=1$ due to the continuous process that prevailed in Fig.
11.a, and the cumulative rate of all discrete events changed from $R_{k,1}$ in
Fig. 11.a to $R_{k,2}$ in the current state. Therefore Eq. 21 yields $\Delta
t_{k}=R_{k,2}^{-1}\left[\ln\left(\frac{1}{u}\right)-\Delta
t^{*}_{1}R_{k,1}\right]$ as written in Fig. 11.b. In some cases Eq. 21 may
return a negative $\Delta t_{k}$; when this happens MASKE issues a warning and
approximates $\Delta t_{k}=0$, thus ensuring that a discrete event will be
carried out next.

Going back to Fig. 11.b, $FOT_{1}$ is now the smallest hence process 1 is
carried out leading to a new $t^{*}$ and $t_{1}$ in Fig. 11.c. In Fig. 11.c,
$FOT_{1}$ and $FOT_{k}$ are updated and the latter happens to be the smallest
one, hence a discrete event is carried. The individual discrete event to
execute (e.g. removing or inserting one specific particle) is selected out of
the list of all possible discrete events contributing to $FOT_{k}$ with a
probability that, for event $i$ out of the $M+N$ possible ones in our example,
is:

$P_{i}\propto\sum_{j=1}^{n_{cp}}{\Delta t^{*}_{j}r_{k,i,j}}+\Delta
t_{k}r_{k,n_{cp}+1}$ (22)

$r_{k,i,j}$ is the rate of event $i$ when the $j^{th}$ continuous process was
executed since the last discrete event took place (i.e. since $t^{*}$ in Fig.
11.a). $r_{k,i,n_{cp}+1}$ is the latest rate of the same event, which
contributed to the cumulative $R_{k,n_{cp}+1}$ in Eq. 21 and that now gets
multiplied times the $\Delta t_{k}$ from Eq. 21.

The execution of a discrete event in Fig. 11.c brings back an analogous
scenario as in Fig. 11.a. This, however, does not mean that the order in which
processes are carried out will remain the same: for example, after Fig. 11.c,
process 2 will be executed twice before calling again process 1.

##  Appendix C Code structure and input script

Fig. 12 shows the code structure of the MASKE software. MASKE features two
levels of parallelization using the MPI library. The user decides how many
processors to involve, $n_{p}$ and assigns them to $n_{sc}$ different sub-
communicators. This assignment is carried out by the MPI_Commm_Split command
in the universe class. Each sub-communicator may be responsible for carrying
out a subset of commands and for sampling a subset of all the discrete events
and compute their individual rates $r$. This is important because sampling
discrete events is usually the most time-consuming task in MASKE. Furthermore,
each sub-communicator may perform its tasks and sampling using multiple
processors in parallel. For example, the user may use 12 processors in total
for a simulation, assigning 8 to sub-communicator $subA$ and 4 to $subB$.
$subA$ may sample the deletion of $N$ particles, which entails computing per-
particle interaction energy using LAMMPS; for this operation, LAMMPS can
operate in parallel using the 8 processors in $subA$. Meanwhile $subB$ may run
a separate instance of LAMMPS using 4 processors to sample the creation of new
particles at $M$ possible locations.

Figure 12: Scheme of the classes forming MASKE and their main relationships to
user commands in the input script.

The main class in the code is maske, which creates and eventually deletes all
the other objects and which hosts several pointers to output files, flags, and
parameters that are useful across the code: see Fig. 12. Another fundamental
class is in pointers.h, where MASKE employs the same strategy as in LAMMPS to
allow any object to directly include and access any other object. Several
classes produce support objects that are used throughout the code: error logs
error messages and stops a simulation when an error occurs; lammpsIO creates
LAMMPS instances and calls LAMMPS operations via one-line strings like those
in a LAMMPS input script; output creates MASKE output files (more details
later) and features functions to write individual entries into them; randm
features functions to generate pseudo-random sequences of numbers (based on
user provided seeds) and to extract numbers from such sequences when needed
(e.g. $u$ in Eq. 21).

The flow of user-defined operations is dictated by an input script, i.e. a
text file that the object inputmsk reads and executes line by line. The two
only exceptions to the line by line execution are the loop and the store multi
commands. The loop command makes MASKE repeat multiple times a number of input
lines: its working is not detailed here for brevity. The store multi command
will be explained later.

The subcomm command is to create the desired number of sub-communicators and,
for each of them, specify a name, assign a number of processors, create a
LAMMPS instance if needed, and provide a seed for a sequence of pseudo-random
numbers.

The chemDB command loads a text file with all the information required to
describe the chemistry of the system and its chemical kinetics. This
information includes data on individual species in solution or forming the
solids, thermodynamic and kinetic parameters such as interfacial energies,
activation energies, stoichiometric coefficients and equilibrium constants of
possible chemical reactions, and mechanism according to which chemical
reactions are combined together in the rate equations. A later dedicated
section will provide more details on the structure and content of the file
loaded by chemDB. The information gathered by chemDB is stored in the
chemistry object.

The real_types and trial_types commands are for the user to list the LAMMPS
particle types associated to real and trial particles. The latter are used for
sampling the possible formation of new particles at various locations in the
simulation box. The vectors of real and trial particle types are stored in the
maske object.

The lammps command passes a string with LAMMPS syntax directly to LAMMPS for
direct execution. This is a very important command to set up a simulation box,
to decide the atom_style and units in LAMMPS, to create particles if needed or
to import a configuration from a previous LAMMPS dump file, to set up the
interaction potentials, to run any preliminary energy minimization, to set up
LAMMPS variables, computes and output files, etc. The $lammps$ command allows
specifying which sub-communicator should execute the line; if the keyword all
is used, all the processors with an active LAMMPS instance will execute the
line using LAMMPS. A type of LAMMPS command that are known not to work through
single-line execution are those related to loops, hence a bespoke loop command
has been implemented in MASKE, as already mentioned.

The kB and hpl commands specify the values of the Boltzmann and Planck
constants in the units chosen by the user; these units should be consistent
with those set in LAMMPS. The values of the constants are recorded in the
maske object.

The sol_start command initializes the implicit solution in MASKE by listing
the names and corresponding molar concentrations of the solvated species, as
well as the of the solvent species. The names listed here must correspond to
species that have been defined in the chemistry file previously read by the
chemDB command. sol_start also provides the temperature of the solution in
Kelvin degrees and the $A$ and $B$ constants of the solvent as per Debye-
Hückel theory, in consistent units (e.g. $A=0.51$ and $B=3.29$ nm-1 for
water). The command also assigns the initial volume of void space in the
simulation box, the additional volume $\Delta V$ of implicit solution, the
volume of voids in $\Delta V$, and a unit conversion factor to convert the
user-decided units of volume into liters (e.g. 10-24 if the user is working in
nm). Having read these inputs, sol_start calls the solution object to compute
the initial volume occupied by the solution in the simulation box and in
$\Delta V$, respectively $SVol$ and $dVSVol$. The former is obtained as
$SVol=V_{box}-voidV-SolidV$, where the volume of the simulation box $V_{box}$
and the volume occupied by solid particles $SolidV$ are both extracted from
LAMMPS, assuming spherical particles. This implies that, in the current
version of MASKE, the user should not change the simulation box size nor
insert, delete, or alter the radii of particles after using the sol_start
command, or MASKE will return an inconsistent $SVol$. Once $SVol$ is computed,
the solution object employs the user-provided initial concentrations to
compute the number of molecules in solution for each solvated species in the
box and then, knowing the apparent volume of each solvated species including
the solvent, the number of solvent molecules is obtained to match $SVol$
itself. An analogous process is followed to compute $dVSVol$ in $\Delta V$ and
the number of solvated molecules for each species there, except that there are
no solid particles in $\Delta V$.

The store command stores information for later, repeated use. Five styles of
storing are currently implemented: four of them store a one-line string each
and are used by the fix_nucleate object, described later, when sampling the
insertion of new particles. These four styles are: region and lattice, which
stores the LAMMPS commands to create respectively a region and a lattice, DV
which stores the LAMMPS command to create a variable quantifying the volume of
a cell in the lattice, and minimize which stores parameters to perform energy
minimization in LAMMPS using a new and bespoke minimizer that MASKE adds as a
user-defined package to LAMMPS when compiling the code. The package is called
USER-MASKE and can be found in the src folder of the LAMMPS compiled as a
submodule of MASKE. This minimizer, called quickmaske, is analogous to
quickmin in LAMMPS but it treats each particle separately when applying and
damping their displacements. This overcomes a limitation of all the minimizers
in LAMMPS, which perform poorly when some particles occupy positions with very
high energy, e.g. largely overlapping with other particles. This is a frequent
scenario when sampling insertion of many new particles at random places and
running a single energy minimization for all of them, as the fix_nucleate
object does in the current version of the code (more details later). The fifth
and last store style is multi, which stores multiple strings for objects such
as fix_CmstoreLMP and fix_EmstoreLMP, described later.

The fix command sets the various processes determining the evolution of the
system in MASKE. Three types of fixes are currently implemented: KMC-free
which sample and execute discrete events using the rejection-free KMC approach
described above, Cont which integrate continuous processes, and Every which
are invoked after $n_{s}$ processes, either discrete or continuous, have been
carried out ($n_{s}$ is input by the user). Depending on the type of fix,
MASKE expects a different set of parameters, which the user must provide in
the input script and that MASKE records in the fix object. Two parameters that
are common to all fixes are the name assigned to the fix and the name of the
sub-communicator that is responsible for sampling and executing it, when
needed. fix Every is the only type of fix which can, and often should, be run
by all the sub-communicators via the all keyword.

Two types of fix KMC-free are currently implemented in MASKE: delete and
nucleate, which respectively sample and execute the dissolution and
precipitation of solid particles through the corresponding objects fix_delete
and fix_nucleate. These objects are crucial in MASKE, hence they will be
described in detail in later sections. For both of them, the user must
specify: (i) the LAMMPS type of real particles that the fix would delete or
create if an event of theirs is executed, (ii) a chemical mechanism from those
defined in the chemistry file loaded by chemDB, to compute the rates of
individual events; (iii) a sol_out option specifying whether the molecules
released into the solution, or removed from it, after executing an event will
be released or removed only in the simulation box (via the box option) or also
in the additional volume $\Delta V$ in proportion to the $SVol/dVSVol$ ratio
(box+dV option). In addition to these, for the nucleate type, the user must
also indicate the names of a stored region, lattice, and lattice cell size, as
well as the LAMMPS type of trial particles for sampling precipitation, and the
diameter of the particles that may precipitate (under the current assumption
in MASKE of spherical particles).

Only one type of fix Cont and one type of fix Every are currently implemented,
both of which simply run a sequence of LAMMPS commands previously stored using
a store multi command. The LAMMPS commands are executed respectively by the
fix_CmstoreLMP and fix_EmstoreLMP objects when appropriate.

The thermo command creates a text file to print output data every $n$ steps in
MASKE, i.e. one entry is recorded whenever $n$ events (discrete or continuous)
have been carried out ($n$ is user-defined, as is the file name). The thermo
file records $m$ user-chosen quantities, organized in a table with $m+2$
columns. The first two columns always record the MASKE step number and the
MASKE time $t^{*}$ corresponding to the quantities in the current entry, i.e.
in the rest of the row being recorded. The quantities in the other $m$ columns
may be of three types: conc_name prints the concentration in solution of a
molecule defined in the chemDB and called name (for example, name may be H2O
or CaCO3); lmp_v_name prints the content of a scalar variable previously
defined in LAMMPS and called name (for example, the user might have defined a
variable called $press$ in LAMMPS, monitoring the pressure in the box);
lmp_c_name is analogous to lmp_v_name but it prints the scalar resulting from
a compute in LAMMPS instead of a variable.

The dump command is the same as the dump command in LAMMPS, in that it creates
a file recording the configuration of the system (i.e. particle positions and
types) as well as various user-decided per-particle quantities such as radii,
stress per atom, etc. The only difference here is that the user inputs the
number of MASKE steps between two successive records of the configuration, as
opposed to the dump command in LAMMPS which prints the configuration every so
many steps in LAMMPS. Both the thermo and dump commands are carried out by the
output object.

The Krun command invokes the krun object to launch a simulation, performing
MASKE steps until covering a user-specified period of time. This means that
the user does not decide how many steps MASKE will perform, but rather by how
much the current Krun will advance $t^{*}$.

##  Appendix D Chemical database

The previous section mentioned that the chemDB command loads a file containing
data on the chemistry of the system. This file is loaded line by line, each
line starting with a specific keyword.

The molecule keyword creates species that can be used to initiate the solution
(see the sol_start command in the previous section) or to define chemical
reactions later in the chemDB file. The keyword is followed by a user-defined
name of the molecule, some quantities that are useful if the molecule is later
associated to a solid phase (namely quantities that describe the size and
shape of the molecule in a solid) and other quantities that are useful for
molecules that are later associated to the solution: apparent volume in
solution, hydrated radius, ionic charge, and the parameter $b$ used in the
Debye-Hückel theory and discussed in the main body of this article.

The gammax keyword is to record a value for the activity coefficient of an
activated complex in standard state, $\gamma^{\ddagger}$. The DGx keyword is
to record the standard state activation free energy of an activated complex,
$\Delta G^{\ddagger}$, as well as the concentration of said complex
$c^{\ddagger}$, and the dimensionality of the concentration: 2 if the per unit
area, 3 if per unit volume. All quantities recorded via gammax and DGx can
later be associated to certain chemical reactions, and will inform the rate
equations from Transition State Theory described in the main body of the
article.

The reax keyword defines a chemical reaction, assigning a user-defined name to
it. The attribute simple indicates that a chemical reaction is just one-step
(a chain attribute can also be used, to combine a series of simple reactions,
but this will not be discussed here for brevity). Take for example the simple
reaction of calcium hydroxide dissolution:

$\mathrm{Ca(OH)_{2}(s)\rightarrow Ca^{2+}(aq)+2OH^{-}(aq)}$ (23)

The reax keyword imports previously recorded values of $\gamma^{\ddagger}$,
$\Delta G^{\ddagger}$, and $c^{\ddagger}$ to associate to this reaction, as
well as an input value for the equilibrium constant $K_{eq}$. The
stoichiometry in Eq. 23 is then recorded by assigning the molecule Ca(OH)2 to
the solid (called ‘foreground’ in MASKE and assigned via the $fgd$ attribute)
with a stoichiometric coefficient -1, because the reaction deletes one such
molecule, and assigning the molecules Ca2+ and OH- to the solution
(‘background’ in MASKE, via the $bkg$ attribute) with stoichiometric
coefficients +1 and +2 respectively. An additional parameter is also needed: a
real number $\chi$ between 0 and 1 used in the rate equations to distribute
the excess free energy from mechanical stress between forward and backward
reactions.

The sen keyword records values of solid-solution interfacial energies that
will later be associated to reaction mechanisms. Such mechanisms are defined
using the mech keyword; they are the last set of entries needed in the chemDB
file and they inform the fix commands of KMC-free type in the main input
script, as discussed in the previous section. A mech keyword requires a user-
defined name for the mechanism and the name of a previously-defined chemical
reaction to attach to it. The user then specifies the type of mechanism, which
is important to determine how a number of chemical reactions of specified time
would assemble to generate (or dissolve) the target particle. For now, only
one type of mechanism is implemented in MASKE, called allpar. This mechanism
assumes that all the reactions required to delete or create a particle take
place in parallel, hence the time to dissolve or form a new particle equals
the time for a single unimolecular reaction to take place. This is exact only
if the volume of a particle is the same as the volume of the individual
molecules involved in a single reaction, e.g. the volume of one Ca(OH)2 solid
molecule in the example above. The allpar mechanism also assumes that the
reaction rates follow the modified version of Transition State Theory
presented in the main body of this article. The last attribute to associate to
a mechanism is whether, in calculating the rates, it should consider straight
or net reactions, also discussed in the main body of the article.

##  Appendix E Krun

The main loop in MASKE, invoked by the Krun command, advances the timescale
through a combination of discrete KMC events and stepwise integration of
continuous processes. As an example, consider a system made of three particle
types, $A$, $B$, and $C$, whose deletion or nucleation are sampled as
discrete, KMC-free events. Consider also that a continuous process (Cont type)
is taking place, e.g. some particles are displaced with a given rate to apply
strain. Finally, consider that after any process is executed, the interaction
energy of all real particles in the system is minimized through a process of
type Every. Let us then assume that two subcommunicators are employed, $subA$
and $subB$, each running two processors, for a total of 4 processors being
employed in the simulation. $subA$ is tasked to sample events for particles of
types $A$ and $C$, namely the deletion of $N_{A}$ and $N_{C}$ real particles
and the nucleation of new particles of type $A$ only at $M_{A}$ possible
locations; let us call these events $fix\\_delA$, $fix\\_delC$, and
$fix\\_nucA$. Similarly, $subB$ samples $N_{B}$ deletions and $M_{B}$
nucleation sites for particles of type $B$, through events $fix\\_delB$ and
$fix\\_nucB$. Processes of Cont and Every type are also tasked to one
subcommunicator (or to all of them) but now it is not important to specify
which one(s) since MASKE eventually synchronizes all processors after a
process is executed.

Fig. 3, in the main body of the article, shows a flow chart for the $Krun$
loop. The first step is for the $fix\\_nucleate$ object to create trial
particles at $M$ locations for each discrete KMC-free event of type $nucleate$
being considered. In our example, $subA$ and $subB$ respectively create
$M_{A}$ and $M_{B}$ trial particles at various locations. The creation of such
particles is managed by the LAMMPS interface, which enables parallelization
via domain decomposition; this means that, since each subcommunicator features
two processors in our example, each processor will manage a subset of the
trial particles, i.e. processors $p1$ and $p_{2}$ forming $subA$ will
respectively create and manage $M_{A,1}$ and $M_{A,2}$ trial particles with
$M_{A,1}+M_{A,2}=M_{A}$, and analogously for $p_{3}$ and $p_{4}$ in $subB$
with its $M_{B}$ trial particles. The $Krun$ loop then proceeds with sampling
all the possible, discrete, KMC-free events to compute their rates. This is
done by the $fix\\_nucleate$ and $fix\\_delete$ objects, which perform
operations that depend on user-specified mechanisms and involve one particle
per event. For example, to inform the rate of a $delete$ event involving
particle $i$ of type $A$ the processor hosting it in $subA$ may compute the
change in interaction energy associated to removing that particle. In this
way, each processor computes the rates $r_{k,i}$ of the events that are
associated to the particles that LAMMPS hosts in that processor itself; Fig.
13 helps clarifying this, with $p1$ for example computing rates for the
deletion of $N_{A,1}$ and $N_{C,1}$ particles and the nucleation of the
$M_{A,1}$ trial particles in it. Each processors stores these rates into two
$rates\\_each$ arrays, one in $fix\\_delete$ and the other one in
$fix\\_nucleate$, as also depicted in Fig. 13. The $rates\\_each$ arrays in
each subcommunicator are then assembled by the submaster processor in it,
leading to two arrays per subcommunicator containing respectively the $rates$
of all $delete$ and $nucleate$ events sampled there: see Fig. 13. The $rates$
arrays are sorted in increasing order of particle ID, these latter being
defined in LAMMPS and extracted from there.

Figure 13: Rate arrays for four discrete KMC-free processes ($fix\\_delA$,
$fix\\_delB$, $fix\\_nucA$, and $fix\\_nucB$) distributed between two
subcommunicators, $subA$ and $subB$, and four processors, $p1-p4$. Processors
$p1$ and $p3$ are the submasters of $subA$ and $subB$; they gather and sort by
particle ID all the rates coming from processor-specific $rate\\_each$ arrays.

The cumulative of the $rates$ arrays across all subcommunicators is the
$R_{k,n_{cp+1}}$ term in Eq. 21, which for the first step in the loop
coincides with $R_{k}$ in Eq. 20 hence suffices to compute $\Delta t_{k}$ and
thus the $FOT_{k}$ of any one discrete event. The $Krun$ loop then samples the
$\Delta t_{i}$ and thus $FOT_{i}$ of all continuous processes, which simply
consists in multiplying together the user-defined integration time step
$dt_{i}$ and number of integration steps $n_{i}$ of each $Cont$ process $i$.
The overall simulation time is thus advanced to
$t^{*}=\min\limits_{k,i}(FOT)$, i.e. the smallest FOT out of $FOT_{k}$ for a
discrete event and $FOT_{i}$ of all $i$ continuous processes. The process with
smallest FOT is then carried out. If the process is of $Cont$ type, its
execution simply entails the integration over $n_{i}$ steps with operations
specified in an object that depends on the type of continuous process, e.g.
$fix\\_CmstoreLMP$ in Fig. 12. The execution of a continuous process is
carried out either by one or by all subcommunicators; for example, processes
of type $fix\\_CmstoreLMP$ are always executed by all subcommunicators.