Interactions¶
Interactions are used to calculate forces on individual particles due to their neighbours. Pairwise shortrange interactions are currently supported, and membrane forces.
Summary¶
Interaction () 
Base interaction class 
MembraneForces () 
Abstract class for membrane interactions. 
ObjBinding () 
Forces attaching a ParticleVector to another via harmonic potentials between the particles of specific pairs. 
ObjRodBinding () 
Forces attaching a RodVector to a RigidObjectVector . 
Pairwise () 
Generic pairwise interaction class. 
RodForces () 
Forces acting on an elastic rod. 
Details¶

class
Interaction
¶ Bases:
object
Base interaction class

__init__
()¶ Initialize self. See help(type(self)) for accurate signature.


class
MembraneForces
¶ Bases:
mmirheo.Interactions.Interaction
Abstract class for membrane interactions. Meshbased forces acting on a membrane according to the model in [Fedosov2010]
 The membrane interactions are composed of forces comming from:
 bending of the membrane, potential \(U_b\)
 shear elasticity of the membrane, potential \(U_s\)
 constraint: area conservation of the membrane (local and global), potential \(U_A\)
 constraint: volume of the cell (assuming incompressible fluid), potential \(U_V\)
 membrane viscosity, pairwise force \(\mathbf{F}^v\)
 membrane fluctuations, pairwise force \(\mathbf{F}^R\)
The form of the constraint potentials are given by (see [Fedosov2010] for more explanations):
\[\begin{split}U_A = \frac{k_a (A_{tot}  A^0_{tot})^2}{2 A^0_{tot}} + \sum_{j \in {1 ... N_t}} \frac{k_d (A_jA_0)^2}{2A_0}, \\ U_V = \frac{k_v (VV^0_{tot})^2}{2 V^0_{tot}}.\end{split}\]The viscous and dissipation forces are central forces and are the same as DPD interactions with \(w(r) = 1\) (no cutoff radius, applied to each bond).
Several bending models are implemented. First, the Kantor enrgy reads (see [kantor1987]):
\[U_b = \sum_{j \in {1 ... N_s}} k_b \left[ 1\cos(\theta_j  \theta_0) \right].\]The Juelicher energy is (see [Juelicher1996]):
\[\begin{split}U_b = 2 k_b \sum_{\alpha = 1}^{N_v} \frac {\left( M_{\alpha}  C_0\right)^2}{A_\alpha}, \\ M_{\alpha} = \frac 1 4 \sum_{<i,j>}^{(\alpha)} l_{ij} \theta_{ij}.\end{split}\]It is improved with the areadifference model (see [Bian2020]), which is a discretized version of:
\[U_{AD} = \frac{k_{AD} \pi}{2 D_0^2 A_0} \left(\Delta A  \Delta A_0 \right)^2.\]Currently, the stretching and shear energy models are:
WLC model:
\[U_s = \sum_{j \in {1 ... N_s}} \left[ \frac {k_s l_m \left( 3x_j^2  2x_j^3 \right)}{4(1x_j)} + \frac{k_p}{l_0} \right].\]Lim model: an extension of the Skalak shear energy (see [Lim2008]).
\[\begin{split}U_{Lim} =& \sum_{i=1}^{N_{t}}\left(A_{0}\right)_{i}\left(\frac{k_a}{2}\left(\alpha_{i}^{2}+a_{3} \alpha_{i}^{3}+a_{4} \alpha_{i}^{4}\right)\right.\\ & +\mu\left(\beta_{i}+b_{1} \alpha_{i} \beta_{i}+b_{2} \beta_{i}^{2}\right) ),\end{split}\]where \(\alpha\) and \(\beta\) are the invariants of the strains.
[Fedosov2010] (1, 2) Fedosov, D. A.; Caswell, B. & Karniadakis, G. E. A multiscale red blood cell model with accurate mechanics, rheology, and dynamics Biophysical journal, Elsevier, 2010, 98, 22152225 [kantor1987] Kantor, Y. & Nelson, D. R. Phase transitions in flexible polymeric surfaces Physical Review A, APS, 1987, 36, 4020 [Juelicher1996] Juelicher, Frank, and Reinhard Lipowsky. Shape transformations of vesicles with intramembrane domains. Physical Review E 53.3 (1996): 2670. [Bian2020] Bian, Xin, Sergey Litvinov, and Petros Koumoutsakos. Bending models of lipid bilayer membranes: Spontaneous curvature and areadifference elasticity. Computer Methods in Applied Mechanics and Engineering 359 (2020): 112758. [Lim2008] Lim HW, Gerald, Michael Wortis, and Ranjan Mukhopadhyay. Red blood cell shapes and shape transformations: newtonian mechanics of a composite membrane: sections 2.1–2.4. Soft Matter: Lipid Bilayers and Red Blood Cells 4 (2008): 83139. 
__init__
(name: str, shear_desc: str, bending_desc: str, filter_desc: str='keep_all', stress_free: bool=False, **kwargs) → None¶ Parameters:  name – name of the interaction
 shear_desc – a string describing what shear force is used
 bending_desc – a string describing what bending force is used
 filter_desc – a string describing which membranes are concerned
 stress_free – if True, stress Free shape is used for the shear parameters
kwargs:
 tot_area: total area of the membrane at equilibrium
 tot_volume: total volume of the membrane at equilibrium
 ka_tot: constraint energy for total area
 kv_tot: constraint energy for total volume
 kBT: fluctuation temperature (set to zero will switch off fluctuation forces)
 gammaC: dissipative forces coefficient
 initial_length_fraction: the size of the membrane increases linearly in time from this fraction of the provided mesh to its full size after grow_until time; the parameters are scaled accordingly with time. If this is set, grow_until must also be provided. Default value: 1.
 grow_until: the size increases linearly in time from a fraction of the provided mesh to its full size after that time; the parameters are scaled accordingly with time. If this is set, initial_length_fraction must also be provided. Default value: 0
Shear Parameters, warm like chain model (set shear_desc = ‘wlc’):
 x0: \(x_0\)
 ks: energy magnitude for bonds
 mpow: \(m\)
 ka: energy magnitude for local area
Shear Parameters, Lim model (set shear_desc = ‘Lim’):
 ka: \(k_a\), magnitude of stretching force
 mu: \(\mu\), magnitude of shear force
 a3: \(a_3\), non linear part for stretching
 a4: \(a_4\), non linear part for stretching
 b1: \(b_1\), non linear part for shear
 b2: \(b_2\), non linear part for shear
Bending Parameters, Kantor model (set bending_desc = ‘Kantor’):
 kb: local bending energy magnitude
 theta: spontaneous angle
Bending Parameters, Juelicher model (set bending_desc = ‘Juelicher’):
 kb: local bending energy magnitude
 C0: spontaneous curvature
 kad: area difference energy magnitude
 DA0: area difference at relaxed state divided by the offset of the leaflet midplanes
filter_desc = “keep_all”:
The interaction will be applied to all membranesfilter_desc = “by_type_id”:
The interaction will be applied membranes with a given type_id (see
MembraneWithTypeId
) type_id: the type id that the interaction applies to

class
ObjBinding
¶ Bases:
mmirheo.Interactions.Interaction
Forces attaching a
ParticleVector
to another via harmonic potentials between the particles of specific pairs.Warning
To deal with MPI, the force is zero if two particles of a pair are apart from more than half the subdomain size. Since this interaction is designed to bind objects to each other, this should not happen under normal conditions.

__init__
(name: str, k_bound: float, pairs: List[int2]) → None¶ Parameters:  name – Name of the interaction.
 k_bound – Spring force coefficient.
 pairs – The global Ids of the particles that will interact through the harmonic potential. For each pair, the first entry is the id of pv1 while the second is that of pv2 (see
setInteraction
).


class
ObjRodBinding
¶ Bases:
mmirheo.Interactions.Interaction
Forces attaching a
RodVector
to aRigidObjectVector
.
__init__
(name: str, torque: float, rel_anchor: real3, k_bound: float) → None¶ Parameters:  name – name of the interaction
 torque – torque magnitude to apply to the rod
 rel_anchor – position of the anchor relative to the rigid object
 k_bound – anchor harmonic potential magnitude


class
Pairwise
¶ Bases:
mmirheo.Interactions.Interaction
Generic pairwise interaction class. Can be applied between any kind of
ParticleVector
classes. The following interactions are currently implemented: DPD:
Pairwise interaction with conservative part and dissipative + random part acting as a thermostat, see [Groot1997]
\[\begin{split}\mathbf{F}_{ij} &= \left(\mathbf{F}^C_{ij} + \mathbf{F}^D_{ij} + \mathbf{F}^R_{ij} \right) \mathbf{\hat r} \\ F^C_{ij} &= \begin{cases} a(1\frac{r}{r_c}), & r < r_c \\ 0, & r \geqslant r_c \end{cases} \\ F^D_{ij} &= \gamma w^2(\tfrac{r}{r_c}) (\mathbf{\hat r} \cdot \mathbf{u}) \\ F^R_{ij} &= \sigma w(\tfrac{r}{r_c}) \, \frac{\theta}{\sqrt{\Delta t}} \,\end{split}\]where bold symbol means a vector, its regular counterpart means vector length: \(x = \left\lVert \mathbf{x} \right\rVert\), hated symbol is the normalized vector: \(\mathbf{\hat x} = \mathbf{x} / \left\lVert \mathbf{x} \right\rVert\). Moreover, \(\theta\) is the random variable with zero mean and unit variance, that is distributed independently of the interacting pair ij, dissipation and random forces are related by the fluctuationdissipation theorem: \(\sigma^2 = 2 \gamma \, k_B T\); and \(w(r)\) is the weight function that we define as follows:
\[\begin{split}w(r) = \begin{cases} (1r)^{p}, & r < 1 \\ 0, & r \geqslant 1 \end{cases}\end{split}\]
 MDPD:
Compute MDPD interaction as described in [Warren2003]. Must be used together with “Density” interaction with kernel “MDPD”.
The interaction forces are the same as described in “DPD” with the modified conservative term
\[F^C_{ij} = a w_c(r_{ij}) + b (\rho_i + \rho_j) w_d(r_{ij}),\]where \(\rho_i\) is computed from “Density” and
\[\begin{split}w_c(r) = \begin{cases} (1\frac{r}{r_c}), & r < r_c \\ 0, & r \geqslant r_c \end{cases} \\ w_d(r) = \begin{cases} (1\frac{r}{r_d}), & r < r_d \\ 0, & r \geqslant r_d \end{cases}\end{split}\]
 SDPD:
Compute SDPD interaction with angular momentum conservation, following [Hu2006] and [Bian2012]. Must be used together with “Density” interaction with the same density kernel.
\[\begin{split}\mathbf{F}_{ij} &= \left(F^C_{ij} + F^D_{ij} + F^R_{ij} \right) \\ F^C_{ij} &=  \left( \frac{p_{i}}{d_{i}^{2}}+\frac{p_{j}}{d_{j}^{2}}\right) \frac{\partial w_\rho}{\partial r_{ij}}, \\ F^D_{ij} &=  \eta \left[ \left(\frac{1}{d_{i}^{2}}+\frac{1}{d_{j}^{2}}\right) \frac{\zeta}{r_{ij}} \frac{\partial w_\rho}{\partial r_{ij}}\right] \left( \mathbf{v}_{i j} \cdot \mathbf{e}_{ij} \right), \\ F^R_{ij} &= \sqrt{2 k_BT \eta} \left[ \left(\frac{1}{d_{i}^{2}}+\frac{1}{d_{j}^{2}}\right) \frac{\zeta}{r_{ij}} \frac{\partial w_\rho}{\partial r_{ij}}\right]^{\frac 1 2} \xi_{ij},\end{split}\]where \(\eta\) is the viscosity, \(w_\rho\) is the density kernel, \(\zeta = 2+d = 5\), \(d_i\) is the density of particle i and \(p_i = p(d_i)\) is the pressure of particle i.. The available density kernels are listed in “Density”. The available equations of state (EOS) are:
Linear equation of state:
\[p(\rho) = c_S^2 \left(\rho  \rho_0 \right)\]where \(c_S\) is the speed of sound and \(\rho_0\) is a parameter.
Quasi incompressible EOS:
\[p(\rho) = p_0 \left[ \left( \frac {\rho}{\rho_r} \right)^\gamma  1 \right],\]where \(p_0\), \(\rho_r\) and \(\gamma = 7\) are parameters to be fitted to the desired fluid.
 LJ:
Pairwise interaction according to the classical LennardJones potential
\[\mathbf{F}_{ij} = 24 \epsilon \left( 2\left( \frac{\sigma}{r_{ij}} \right)^{12}  \left( \frac{\sigma}{r_{ij}} \right)^{6} \right) \frac{\mathbf{r}}{r^2}\]As opposed to
RepulsiveLJ
, the force is not bounded from either sides.
 RepulsiveLJ:
Pairwise interaction according to the classical LennardJones potential, truncated such that it is always repulsive.
\[\mathbf{F}_{ij} = \max \left[ 0.0, 24 \epsilon \left( 2\left( \frac{\sigma}{r_{ij}} \right)^{12}  \left( \frac{\sigma}{r_{ij}} \right)^{6} \right) \frac{\mathbf{r}}{r^2} \right]\]Note that in the implementation, the force is bounded for stability at larger time steps.
 GrowingRepulsiveLJ:
Same as RepulsiveLJ, but the length scale is growing linearly in time until a prespecified time, from a specified fraction to 1. This is useful when growing membranes while avoiding overlaps.
 Morse:
Pairwise interaction according to the classical Morse potential
\[\mathbf{F}_{ij} = 2 D_e \beta \left( e^{2 \beta (r_0r)}  e^{\beta (r_0r)} \right) \frac{\mathbf{r}}{r},\]where \(r\) is the distance between the particles.
 Density:
Compute density of particles with a given kernel.
\[\rho_i = \sum\limits_{j \neq i} w_\rho (r_{ij})\]where the summation goes over the neighbours of particle \(i\) within a cutoff range of \(r_c\). The implemented densities are listed below:
kernel “MDPD”:
see [Warren2003]
\[\begin{split}w_\rho(r) = \begin{cases} \frac{15}{2\pi r_d^3}\left(1\frac{r}{r_d}\right)^2, & r < r_d \\ 0, & r \geqslant r_d \end{cases}\end{split}\]kernel “WendlandC2”:
\[w_\rho(r) = \frac{21}{2\pi r_c^3} \left( 1  \frac{r}{r_c} \right)^4 \left( 1 + 4 \frac{r}{r_c} \right)\]
[Groot1997] Groot, R. D., & Warren, P. B. (1997). Dissipative particle dynamics: Bridging the gap between atomistic and mesoscopic simulations. J. Chem. Phys., 107(11), 44234435. doi <https://doi.org/10.1063/1.474784> [Warren2003] Warren, P. B. “Vaporliquid coexistence in manybody dissipative particle dynamics.” Physical Review E 68.6 (2003): 066702. [Hu2006] Hu, X. Y., and N. A. Adams. “Angularmomentum conservative smoothed particle dynamics for incompressible viscous flows.” Physics of Fluids 18.10 (2006): 101702. [Bian2012] Bian, Xin, et al. “Multiscale modeling of particle in suspension with smoothed dissipative particle dynamics.” Physics of Fluids 24.1 (2012): 012002. 
__init__
(name: str, rc: float, kind: str, **kwargs) → None¶ Parameters:  name – name of the interaction
 rc – interaction cutoff (no forces between particles further than rc apart)
 kind – interaction kind (e.g. DPD). See below for all possibilities.
Create one pairwise interaction handler of kind kind. When applicable, stress computation is activated by passing stress = True. This activates virial stress computation every stress_period time units (also passed in kwars)
kind = “DPD”
 a: \(a\)
 gamma: \(\gamma\)
 kBT: \(k_B T\)
 power: \(p\) in the weight function
kind = “MDPD”
 rd: \(r_d\)
 a: \(a\)
 b: \(b\)
 gamma: \(\gamma\)
 kBT: temperature \(k_B T\)
 power: \(p\) in the weight function
kind = “SDPD”
 viscosity: fluid viscosity
 kBT: temperature \(k_B T\)
 EOS: the desired equation of state (see below)
 density_kernel: the desired density kernel (see below)
kind = “LJ”
 epsilon: \(\varepsilon\)
 sigma: \(\sigma\)
kind = “RepulsiveLJ”
 epsilon: \(\varepsilon\)
 sigma: \(\sigma\)
 max_force: force magnitude will be capped to not exceed max_force
 aware_mode:
 if “None”, all particles interact with each other.
 if “Object”, the particles belonging to the same object in an object vector do not interact with each other. That restriction only applies if both Particle Vectors in the interactions are the same and is actually an Object Vector.
 if “Rod”, the particles interact with all other particles except with the ones which are below a given a distance (in number of segment) of the same rod vector. The distance is specified by the kwargs parameter min_segments_distance.
kind = “GrowingRepulsiveLJ”
 epsilon: \(\varepsilon\)
 sigma: \(\sigma\)
 max_force: force magnitude will be capped to not exceed max_force
 aware_mode:
 if “None”, all particles interact with each other.
 if “Object”, the particles belonging to the same object in an object vector do not interact with each other. That restriction only applies if both Particle Vectors in the interactions are the same and is actually an Object Vector.
 if “Rod”, the particles interact with all other particles except with the ones which are below a given a distance (in number of segment) of the same rod vector. The distance is specified by the kwargs parameter min_segments_distance.
 init_length_fraction: tnitial length factor. Must be in [0, 1].
 grow_until: time after which the length quantities are scaled by one.
kind = “Morse”
 De: \(D_e\)
 r0: \(r_0\)
 beta: \(\beta\)
 aware_mode: See “RepulsiveLJ” kernel description.
kind = “Density”
 density_kernel: the desired density kernel (see below)
The available density kernels are “MDPD” and “WendlandC2”. Note that “MDPD” can not be used with SDPD interactions. MDPD interactions can use only “MDPD” density kernel.
For SDPD, the available equation of states are given below:
EOS = “Linear” parameters:
 sound_speed: the speed of sound
 rho_0: background pressure in \(c_S\) units
EOS = “QuasiIncompressible” parameters:
 p0: \(p_0\)
 rho_r: \(\rho_r\)

class
RodForces
¶ Bases:
mmirheo.Interactions.Interaction
Forces acting on an elastic rod.
 The rod interactions are composed of forces comming from:
 bending energy, \(E_{\text{bend}}\)
 twist energy, \(E_{\text{twist}}\)
 bounds energy, \(E_{\text{bound}}\)
The form of the bending energy is given by (for a bisegment):
\[E_{\mathrm{bend}}=\frac{l}{4} \sum_{j=0}^{1}\left(\kappa^{j}\overline{\kappa}\right)^{T} B\left(\kappa^{j}\overline{\kappa}\right),\]where
\[\kappa^{j}=\frac {1} {l} \left((\kappa \mathbf{b}) \cdot \mathbf{m}_{2}^{j},(\kappa \mathbf{b}) \cdot \mathbf{m}_{1}^{j}\right).\]See, e.g. [bergou2008] for more details. The form of the twist energy is given by (for a bisegment):
\[E_{\mathrm{twist}}=\frac{k_{t} l}{2}\left(\frac{\theta^{1}\theta^{0}}{l}\overline{\tau}\right)^{2}.\]The additional bound energy is a simple harmonic potential with a given equilibrium length.
[bergou2008] Bergou, M.; Wardetzky, M.; Robinson, S.; Audoly, B. & Grinspun, E. Discrete elastic rods ACM transactions on graphics (TOG), 2008, 27, 63 
__init__
(name: str, state_update: str='none', save_energies: bool=False, **kwargs) → None¶ Parameters:  name – name of the interaction
 state_update – description of the state update method; only makes sense for multiple states. See below for possible choices.
 save_energies – if True, save the energies of each bisegment
kwargs:
 a0 (real): equilibrium length between 2 opposite cross vertices
 l0 (real): equilibrium length between 2 consecutive vertices on the centerline
 k_s_center (real): elastic force magnitude for centerline
 k_s_frame (real): elastic force magnitude for material frame particles
 k_bending (real3): Bending symmetric tensor \(B\) in the order \(\left(B_{xx}, B_{xy}, B_{zz} \right)\)
 kappa0 (real2): Spontaneous curvatures along the two material frames \(\overline{\kappa}\)
 k_twist (real): Twist energy magnitude \(k_\mathrm{twist}\)
 tau0 (real): Spontaneous twist \(\overline{\tau}\)
 E0 (real): (optional) energy ground state
state update parameters, for state_update = ‘smoothing’:
(not fully implemented yet; for now just takes minimum state but no smoothing term)state update parameters, for state_update = ‘spin’:
 nsteps number of MC step per iteration
 kBT temperature used in the acceptancerejection algorithm
 J neighbouring spin ‘dislike’ energy
The interaction can support multiple polymorphic states if kappa0, tau0 and E0 are lists of equal size. In this case, the E0 parameter is required. Only lists of 1, 2 and 11 states are supported.