rearrangement; thus, they cannot be captured by any standard pairwise interaction. A
simple solution for modelling bond swaps is to embed Monte Carlo hops into hybrid
molecular dynamics or Monte Carlo (MD,MC) simulations [77–80], but this breaks
the continuous time of MD, and it requires fine tweaking in order to follow the real
dynamics of the system. We propose instead a fully MD method based on the
implementation of a three-body potential [81] that already has been shown to
provide meaningful results in the context of vitrimers [66, 67, 82] as well as a
trick to speed up the equilibration of strong network formers [83]. This additional
three-body potential is built upon the specific pairwise interaction of the system and
enriches it with a mechanism to rearrange bonds formed by this pairwise attraction.
At the same time, it also provides single bond per site condition, preventing the
clusterization of the whole system around a highly connected region. To show how
the three-body mechanism works, we will start picking an optimal pairwise interaction in the form of a generalized Lennard-Jones:
U
2
ð Þ
GLJ r ij
À Á ¼ 4ε
σ
r ij
2n
À
σ
r ij
n
"
#
r < r cut ,
ð43Þ
where as usual r ij ¼ j r
!
i À r
!
j j and r cut is some cut-off distance. Its parameters
have to mimic the physics of the system, so if the bonds are covalent, we will pick
the bond energy ε ) k B T and n ! 6 in order the make short-range bonds with an
equilibrium distance of r min ¼ 2
1/n
σ. The choice of a large n is particularly relevant in
order to achieve computational efficiency: since the main cost associated with the
evaluation of the three-body potential comes from the number of interacting triplets
to account, it is desirable to reduce it by selecting the smallest cut-off r c , consistent
with the speed at which the interaction vanishes.
After the pairwise interaction is set up, the three-body term quantifies how much
the force between monomer i and j is affected by the presence of an additional
monomer k within range of interaction. The effect of this interacting triplet is
captured by the following three-body exchange (E3) potential
U
3
ð Þ
E3 r
!
i , r
!
j , r
!
k ; λ
¼ λε b
U
2
ð Þ
E2 r ij
À Á Â b
U
2
ð Þ
E2 r ik
ð Þ,
ð44Þ
which is the product of two two-body exchange potentials related to the generalized
Lennard-Jones potential as
b
U
2
ð Þ
E2 r ij
À Á ¼
1
r r min
À 1=ε
ð ÞU
2
ð Þ
GLJ r ij
À Á
r > r min :
&
ð45Þ
In Eq. (44) we introduced the three-body energy parameter λ that directly controls
the swap rate by tuning the energy required for a swap event to happen, thus
mimicking the process that controls swap, e.g. catalyst concentration. Its role is
made clear in the following example where Eq. (44) is rewritten for a monomer a at
Rheology, Rupture, Reinforcement and Reversibility: Computational Approaches. . .
87
simple solution for modelling bond swaps is to embed Monte Carlo hops into hybrid
molecular dynamics or Monte Carlo (MD,MC) simulations [77–80], but this breaks
the continuous time of MD, and it requires fine tweaking in order to follow the real
dynamics of the system. We propose instead a fully MD method based on the
implementation of a three-body potential [81] that already has been shown to
provide meaningful results in the context of vitrimers [66, 67, 82] as well as a
trick to speed up the equilibration of strong network formers [83]. This additional
three-body potential is built upon the specific pairwise interaction of the system and
enriches it with a mechanism to rearrange bonds formed by this pairwise attraction.
At the same time, it also provides single bond per site condition, preventing the
clusterization of the whole system around a highly connected region. To show how
the three-body mechanism works, we will start picking an optimal pairwise interaction in the form of a generalized Lennard-Jones:
U
2
ð Þ
GLJ r ij
À Á ¼ 4ε
σ
r ij
2n
À
σ
r ij
n
"
#
r < r cut ,
ð43Þ
where as usual r ij ¼ j r
!
i À r
!
j j and r cut is some cut-off distance. Its parameters
have to mimic the physics of the system, so if the bonds are covalent, we will pick
the bond energy ε ) k B T and n ! 6 in order the make short-range bonds with an
equilibrium distance of r min ¼ 2
1/n
σ. The choice of a large n is particularly relevant in
order to achieve computational efficiency: since the main cost associated with the
evaluation of the three-body potential comes from the number of interacting triplets
to account, it is desirable to reduce it by selecting the smallest cut-off r c , consistent
with the speed at which the interaction vanishes.
After the pairwise interaction is set up, the three-body term quantifies how much
the force between monomer i and j is affected by the presence of an additional
monomer k within range of interaction. The effect of this interacting triplet is
captured by the following three-body exchange (E3) potential
U
3
ð Þ
E3 r
!
i , r
!
j , r
!
k ; λ
¼ λε b
U
2
ð Þ
E2 r ij
À Á Â b
U
2
ð Þ
E2 r ik
ð Þ,
ð44Þ
which is the product of two two-body exchange potentials related to the generalized
Lennard-Jones potential as
b
U
2
ð Þ
E2 r ij
À Á ¼
1
r r min
À 1=ε
ð ÞU
2
ð Þ
GLJ r ij
À Á
r > r min :
&
ð45Þ
In Eq. (44) we introduced the three-body energy parameter λ that directly controls
the swap rate by tuning the energy required for a swap event to happen, thus
mimicking the process that controls swap, e.g. catalyst concentration. Its role is
made clear in the following example where Eq. (44) is rewritten for a monomer a at
Rheology, Rupture, Reinforcement and Reversibility: Computational Approaches. . .
87