Skip to content

Update for frictional contacts #2184

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
wants to merge 3 commits into
base: trunk
Choose a base branch
from
Open

Update for frictional contacts #2184

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
3 commits
Select commit Hold shift + click to select a range
9f5240d
Changed the moment of inertia x/y to z in the calculation of nu and t…
Nov 20, 2025
7fa9f12
Updated the Documentation
Nov 21, 2025
5d0c92d
Updated Documentation
Nov 24, 2025
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
8 changes: 4 additions & 4 deletions hoomd/md/FrictionPair.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -353,8 +353,8 @@ void FrictionPair<friction_evaluator>::computeForces(uint64_t timestep)
quat<Scalar> angmom_i(h_angmom.data[i]);
quat<Scalar> qxP_i = conj(quat_i) * angmom_i;
vec3<Scalar> bf_vel_i = Scalar(0.5)
* vec3(qxP_i.v.x / h_moment_inertia.data[i].x,
qxP_i.v.y / h_moment_inertia.data[i].y,
* vec3(qxP_i.v.x / h_moment_inertia.data[i].z,
qxP_i.v.y / h_moment_inertia.data[i].z,
qxP_i.v.z / h_moment_inertia.data[i].z);

// Rotate angular velocity into global frame
Expand Down Expand Up @@ -412,8 +412,8 @@ void FrictionPair<friction_evaluator>::computeForces(uint64_t timestep)
quat<Scalar> angmom_j(h_angmom.data[j]);
quat<Scalar> qxP_j = conj(quat_j) * angmom_j;
vec3<Scalar> bf_vel_j = Scalar(0.5)
* vec3(qxP_j.v.x / h_moment_inertia.data[j].x,
qxP_j.v.y / h_moment_inertia.data[j].y,
* vec3(qxP_j.v.x / h_moment_inertia.data[j].z,
qxP_j.v.y / h_moment_inertia.data[j].z,
qxP_j.v.z / h_moment_inertia.data[j].z);

// Rotate angular velocity into global frame
Expand Down
12 changes: 6 additions & 6 deletions hoomd/md/FrictionPairGPU.cuh
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -270,8 +270,8 @@ gpu_compute_pair_friction_forces_kernel(Scalar4* d_force,
// calculate angular momentum of i-th particle in the body frame
quat<Scalar> qxP_i = conj(quati) * angmomi;
vec3<Scalar> bf_vel_i = Scalar(0.5)
* vec3<Scalar>(qxP_i.v.x / moment_inertia_i.x,
qxP_i.v.y / moment_inertia_i.y,
* vec3<Scalar>(qxP_i.v.x / moment_inertia_i.z,
qxP_i.v.y / moment_inertia_i.z,
qxP_i.v.z / moment_inertia_i.z);

// Rotate angular velocity into global frame
Expand Down Expand Up @@ -320,8 +320,8 @@ gpu_compute_pair_friction_forces_kernel(Scalar4* d_force,
// calculate angular momentum of i-th particle in the body frame
quat<Scalar> qxP_j = conj(quatj) * angmomj;
vec3<Scalar> bf_vel_j = Scalar(0.5)
* vec3<Scalar>(qxP_j.v.x / moment_inertia_j.x,
qxP_j.v.y / moment_inertia_j.y,
* vec3<Scalar>(qxP_j.v.x / moment_inertia_j.z,
qxP_j.v.y / moment_inertia_j.z,
qxP_j.v.z / moment_inertia_j.z);

// Rotate angular velocity into global frame
Expand Down Expand Up @@ -372,8 +372,8 @@ gpu_compute_pair_friction_forces_kernel(Scalar4* d_force,
{
// Calculate nu for the Ito formalism
Scalar nu_ito = ((Scalar(1.0) / massi) + (Scalar(1.0) / massj))
+ (((diaj * diaj / Scalar(4.0)) / moment_inertia_j.x)
+ ((diai * diai / Scalar(4.0)) / moment_inertia_i.x));
+ (((diaj * diaj / Scalar(4.0)) / moment_inertia_j.z)
+ ((diai * diai / Scalar(4.0)) / moment_inertia_i.z));
eval.setNu(nu_ito);
}

Expand Down
182 changes: 111 additions & 71 deletions hoomd/md/pair/friction.py
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -1,60 +1,100 @@
# Copyright (c) 2009-2025 The Regents of the University of Michigan.
# Part of HOOMD-blue, released under the BSD 3-Clause License.

r"""Frictional pair force classes apply a force, and torque on every
r"""Frictional pair force classes apply a force and torque on every
particle in the simulation state. The following general expression
for Markovian tangential friction forces is implemented for
interactions between two spherical particles and is discussed in detail in
`Hofmann et. al. 2025`_. For two particles :math:`i` and
:math:`j` with radii :math:`R_{i,j}`, center positions :math:`\mathbf{r}_{i,j}`,
angular velocities :math:`\mathbf{\omega}_{i,j}`, and translational velocities
:math:`\mathbf{v}_{i,j}`, their surface velocities at the contact point are given by
`Hofmann et al. 2025`_.

.. _Hofmann et. al. 2025: https://doi.org/10.48550/arXiv.2507.16388
For two particles :math:`i` and :math:`j` with radii :math:`R_{i,j}`,
center positions :math:`\mathbf{r}_{i,j}`, angular velocities
:math:`\mathbf{\omega}_{i,j}`, and translational velocities :math:`\mathbf{v}_{i,j}`,
their surface velocities at the contact point are given by

.. _Hofmann et al. 2025: https://doi.org/10.48550/arXiv.2507.16388

.. math::

\begin{align*}
\mathbf{u}_i &= \mathbf{v}_i+\mathbf{\omega}_i \times \mathbf{\hat{r}}_{ij}R_i \\
\mathbf{u}_i &= \mathbf{v}_i+\mathbf{\omega}_i \times \mathbf{\hat{r}}_{ij}R_i, \\
\mathbf{u}_j &= \mathbf{v}_j-\mathbf{\omega}_j \times \mathbf{\hat{r}}_{ij}R_j \, ,
\end{align*}

where :math:`\mathbf{\hat{r}}_{ij}=\mathbf{r}_{ij}/r_{i,j}`. With these expressions,
where :math:`\mathbf{\hat{r}}_{ij}=\mathbf{r}_{ij}/r_{ij}`. With these expressions,
we calculate the relative tangential velocity :math:`\mathbf{u}^\perp_{i,j}` at the
contact point

.. math::

\mathbf{u}^\perp_{i,j} = \mathbf{P}(\mathbf{\hat{r}}_{ij})(\mathbf{v}_j
-\mathbf{v}_i) -(\mathbf{\omega}_iR_i+\mathbf{\omega}_jR_j)
\times \mathbf{\hat{r}}_{ij}\, ,
\times \mathbf{\hat{r}}_{ij},

with the projection operator

.. math::

\mathbf{P}(\mathbf{\hat{r}}_{ij})=1-\mathbf{\hat{r}}_{ij}\mathbf{\hat{r}}_{ij}.

wit the projection operator :math:`\mathbf{P}(\mathbf{\hat{r}}_{ij})=1
-\mathbf{\hat{r}}_{ij}\mathbf{\hat{r}}_{ij}`. We model the tangential
friction force at the contact point very general as
We model the tangential friction force at the contact point very generally as

.. math::

\mathbf{F}^\mathrm{f,contact}_i = -\mathbf{F}^\mathrm{f,contact}_j = f(u^\perp_{i,j}
,r_{i,j})\mathbf{\hat{u}}^\perp_{i,j}

with :math:`\mathbf{\hat{u}}^\perp_{i,j}=\mathbf{u}^\perp_{i,j}/u^\perp_{i,j}`, where
:math:`f(u^\perp_{i,j},r_{i,j})` is an arbitrary function. The surface force
:math:`\mathbf{F}^\mathrm{f,contact}_i` generates a center-of-mass force and torque
acting on particle :math:`i`,
where :math:`\mathbf{\hat{u}}^\perp_{i,j}=\mathbf{u}^\perp_{i,j}/u^\perp_{i,j}`, and
:math:`f(u^\perp_{i,j},r_{i,j})` is an arbitrary scalar function. The functional form of
:math:`f(u^\perp_{i,j},r_{i,j})` specifies the frictional model.

In addition, a stochastic force satisfying the fluctuation-dissipation relation can
be included. It has the form

.. math::

\mathbf{F}^\mathrm{R,contact}_{i} = -\mathbf{F}^\mathrm{R,contact}_j =
\sqrt{D(u^\perp_{ij},r_{ij})}\Big[\mathbf{P}
(\mathbf{\hat{r}}_{ij})\mathbf{\xi}_{ij}
- \mathbf{\hat{r}}_{ij} \times \mathbf{N}
_{ij}\Big]

where :math:`\mathbf{\xi}_{ij}` and :math:`\mathbf{N}_{ij}` are Gaussian white noise
vectors with correlations

.. math::

\begin{align*}
\mathbf{F}^\mathrm{f}_{ij} &= \mathbf{F}^\mathrm{f,contact}_i \\
\mathbf{\tau}^\mathrm{f}_{ij} &= R_i\mathbf{\hat{r}}_{ij} \times
\mathbf{F}^\mathrm{f,contact}_i\, ,
\end{align*}
\langle \mathbf{\xi}_{ij} \mathbf{\xi}_{kl}
\rangle &= \mathbf{1}kT
(\delta_{ik}\delta_{jl}-\delta_{il}
\delta_{jk})/\delta t \\
\langle \mathbf{N}_{ij} \mathbf{N}_{kl} \rangle
&= \mathbf{1}kT(\delta_{ik}
\delta_{jl}+\delta_{il}\delta_{jk})/
\delta t \, .
\end{align*}

The function :math:`D(u,r)` is calculated as

which is a pair friction force and torque resulting from the friction with the particle
:math:`j`.
.. math::

D(u,r) = \frac{1}{kT\nu}\int_u^\infty \mathrm{d}u'f(u',r)\mathrm{exp}(-\frac{u'^2
-u^2}{2kT\nu})

with :math:`\nu=(1/m_i+1/m_j)+(R^2_i/I_i+R_j^2/I_j)`.

The suface force :math:`\mathbf{F}_i^\mathrm{contact}=\mathbf{F}^\mathrm{f,contact}_i
+\mathbf{F}^\mathrm{R,contact}_i` generates a center-of-mass force and a torque acting
on particle :math:`i`,

.. math::

\mathbf{F}_{ij} = \mathbf{F}_i^\mathrm{contact},\quad \mathbf{\tau}_{ij}=R_i\hat{
\mathbf{r}}_{ij}\times\mathbf{F}^\mathrm{contact}_i,

The functional form of :math:`f(u^\perp_{i,j},r_{i,j})` specifies the frictional model.
which is the pair frictional contact force and torque resulting from the friction with
particle :math:`j`.

`FrictionalPair` does not support any shifting modes.

Expand Down Expand Up @@ -101,9 +141,9 @@ class FrictionLJCoulomb(FrictionalPair):

`FrictionLJCoulomb` computes the frictional interaction
between pairs of particles with a Coulomb friction model
as described in `Hofmann et. al. 2025`_.
as described in `Hofmann et al. 2025`_.

.. _Hofmann et. al. 2025: https://doi.org/10.48550/arXiv.2507.16388
.. _Hofmann et al. 2025: https://doi.org/10.48550/arXiv.2507.16388

The Coulomb friction model is defined by the function

Expand All @@ -129,38 +169,38 @@ class FrictionLJCoulomb(FrictionalPair):
.. math::

\begin{align*}
D_\mathrm{C}(u,r) &= \frac{w(r)\kappa_\mathrm{f}\sqrt{\pi}}{\sqrt{2k_\mathrm{B}
T\nu}}\mathrm{e}^{\frac{u^2}{2k_\mathrm{B}T\nu}}
D_\mathrm{C}(u,r) &= \frac{w(r)\kappa_\mathrm{f}\sqrt{\pi}}{\sqrt{2kT\nu}}
\mathrm{e}^{\frac{u^2}{2kT\nu}}
\mathrm{Erfc}\big(\frac{u}
{\sqrt{2k_\mathrm{B}T\nu}} \big) \\
{\sqrt{2kT\nu}} \big) \\
\mathbf{F}^\mathrm{R}_{ij} = -\mathbf{F}^\mathrm{R}_{ji} &= \sqrt{D_\mathrm{C}
(u^\perp_{ij},r_{ij})}\Big[\mathbf{P}
(\mathbf{\hat{r}}_{ij})\mathbf{\xi}^\mathrm{f}_{ij}
- \mathbf{\hat{r}}_{ij} \times \mathbf{N}^\mathrm{f}
(\mathbf{\hat{r}}_{ij})\mathbf{\xi}_{ij}
- \mathbf{\hat{r}}_{ij} \times \mathbf{N}
_{ij}\Big] \\
\frac{\mathbf{\tau}^\mathrm{R}_{ij}}{R_i} =
\frac{\mathbf{\tau}^\mathrm{R}_{ji}}{R_j} &=
\sqrt{D_\mathrm{C}(u^\perp_{ij},r_{ij})}\Big[
\mathbf{\hat{r}}_{ij} \times
\mathbf{\xi}^\mathrm{f}_{ij}+\mathbf{P}
(\mathbf{e}_{ij})\mathbf{N}^\mathrm{f}_{ij}\Big]\, ,
\mathbf{\xi}_{ij}+\mathbf{P}
(\mathbf{e}_{ij})\mathbf{N}_{ij}\Big]\, ,
\end{align*}

where :math:`\nu=(1/m_i+1/m_j)+(R_i^2/\mathbf{I}_i+R_j^2/\mathbf{I}_j))`,
and :math:`\mathbf{\xi}^\mathrm{f}_{ij}` and :math:`\mathbf{N}^\mathrm{f}_{ij}`
where :math:`\nu=(1/m_i+1/m_j)+(R_i^2/I_i+R_j^2/I_j))`,
and :math:`\mathbf{\xi}_{ij}` and :math:`\mathbf{N}_{ij}`
are three dimensional Gaussian white noise vectors with correlations

.. math::

\begin{align*}
\langle \mathbf{\xi}^\mathrm{f}_{ij}(t) \mathbf{\xi}^\mathrm{f}_{kl}(t')
\rangle &= \mathbf{1}k_\mathrm{B}T
\langle \mathbf{\xi}_{ij} \mathbf{\xi}_{kl}
\rangle &= \mathbf{1}kT
(\delta_{ik}\delta_{jl}-\delta_{il}
\delta_{jk})\delta(t-t') \\
\langle \mathbf{N}^\mathrm{f}_{ij}(t) \mathbf{N}^\mathrm{f}_{kl}(t') \rangle &=
\mathbf{1}k_\mathrm{B}T(\delta_{ik}
\delta_{jk})/\delta t \\
\langle \mathbf{N}_{ij} \mathbf{N}_{kl} \rangle &=
\mathbf{1}kT(\delta_{ik}
\delta_{jl}+\delta_{il}\delta_{jk})
\delta(t-t')\, .
/\delta t\, .
\end{align*}

.. rubric:: Example:
Expand Down Expand Up @@ -204,9 +244,9 @@ class FrictionLJCoulombNewton(FrictionalPair):

`FrictionLJCoulombNewton` computes the frictional interaction
between pairs of particles with a Coulomb-Newton friction model as described in
`Hofmann et. al. 2025`_.
`Hofmann et al. 2025`_.

.. _Hofmann et. al. 2025: https://doi.org/10.48550/arXiv.2507.16388
.. _Hofmann et al. 2025: https://doi.org/10.48550/arXiv.2507.16388

The Coulomb-Newton friction model is defined by the function

Expand Down Expand Up @@ -234,48 +274,48 @@ class FrictionLJCoulombNewton(FrictionalPair):
\begin{align*}
D_\mathrm{CN}(u,r) &= \begin{cases}
\frac{w(r)\kappa_\mathrm{f}\sqrt{\pi}}
{\sqrt{2k_\mathrm{B}T\nu}}\mathrm{e}^{\frac{u^2}
{2k_\mathrm{B}T\nu}}\mathrm{Erfc}\big(\frac{u}
{\sqrt{2k_\mathrm{B}T\nu}} \big) & ,\, u\ge
{\sqrt{2kT\nu}}\mathrm{e}^{\frac{u^2}
{2kT\nu}}\mathrm{Erfc}\big(\frac{u}
{\sqrt{2kT\nu}} \big) & ,\, u\ge
\frac{w(r)\kappa_\mathrm{f}}{\gamma_\mathrm{f}} \\
\gamma_\mathrm{f}\left(1-\mathrm{e}^{\left(u^2
-(\frac{w(r)\kappa_\mathrm{f}}{\gamma_\mathrm{f}})^2
\right)/2k_\mathrm{B}T\nu} \right) + \frac{w(r)
\kappa_\mathrm{f}\sqrt{\pi}}{\sqrt{2k_\mathrm{B}T\nu}}
\mathrm{e}^\frac{u^2}{2k_\mathrm{B}T\nu}\mathrm{Erfc}
\right)/2kT\nu} \right) + \frac{w(r)
\kappa_\mathrm{f}\sqrt{\pi}}{\sqrt{2kT\nu}}
\mathrm{e}^\frac{u^2}{2kT\nu}\mathrm{Erfc}
\left(\frac{w(r)\kappa_\mathrm{f}/\gamma_\mathrm{f}}
{\sqrt{2k_\mathrm{B}T\nu}}\right) &,\, u < \frac{w(r)
{\sqrt{2kT\nu}}\right) &,\, u < \frac{w(r)
\kappa_\mathrm{f}}{\gamma_\mathrm{f}}
\end{cases} \\
\mathbf{F}^\mathrm{R}_{ij} = -\mathbf{F}^\mathrm{R}_{ji} &= \sqrt{D_\mathrm{CN}
(u^\perp_{ij},r_{ij})}\Big[\mathbf{P}
(\mathbf{\hat{r}}_{ij})
\mathbf{\xi}^\mathrm{f}_{ij} -
\mathbf{\xi}_{ij} -
\mathbf{\hat{r}}_{ij} \times
\mathbf{N}^\mathrm{f}_{ij}\Big] \\
\mathbf{N}_{ij}\Big] \\
\frac{\mathbf{\tau}^\mathrm{R}_{ij}}{R_i} = \frac{\mathbf{\tau}^\mathrm{R}_{ji}}
{R_j} &= \sqrt{D_\mathrm{CN}
(u^\perp_{ij},r_{ij})}\Big[\mathbf{\hat{r}}_{ij}
\times \mathbf{\xi}^\mathrm{f}_{ij}+\mathbf{P}
(\mathbf{e}_{ij})\mathbf{N}^\mathrm{f}_{ij}
\times \mathbf{\xi}_{ij}+\mathbf{P}
(\mathbf{e}_{ij})\mathbf{N}_{ij}
\Big]\, ,
\end{align*}

where :math:`\nu=(1/m_i+1/m_j)+(R_i^2/\mathbf{I}_i+R_j^2/\mathbf{I}_j))`, and
:math:`\mathbf{\xi}^\mathrm{f}_{ij}` and :math:`\mathbf{N}^\mathrm{f}_{ij}` are
where :math:`\nu=(1/m_i+1/m_j)+(R_i^2/I_i+R_j^2/I_j))`, and
:math:`\mathbf{\xi}_{ij}` and :math:`\mathbf{N}_{ij}` are
three dimensional Gaussian white noise vectors with correlations

.. math::

\begin{align*}
\langle \mathbf{\xi}^\mathrm{f}_{ij}(t) \mathbf{\xi}^\mathrm{f}_{kl}(t')
\rangle &= \mathbf{1}k_\mathrm{B}T
\langle \mathbf{\xi}_{ij} \mathbf{\xi}_{kl}
\rangle &= \mathbf{1}kT
(\delta_{ik}\delta_{jl}-\delta_{il}
\delta_{jk})\delta(t-t') \\
\langle \mathbf{N}^\mathrm{f}_{ij}(t) \mathbf{N}^\mathrm{f}_{kl}(t') \rangle &=
\mathbf{1}k_\mathrm{B}T
\delta_{jk})/\delta t \\
\langle \mathbf{N}_{ij} \mathbf{N}_{kl} \rangle &=
\mathbf{1}kT
(\delta_{ik}\delta_{jl}+\delta_{il}
\delta_{jk})\delta(t-t')\, .
\delta_{jk})/\delta t\, .
\end{align*}

.. rubric:: Example:
Expand Down Expand Up @@ -326,9 +366,9 @@ class FrictionLJLinear(FrictionalPair):

`FrictionLJLinear` computes the frictional interaction
between pairs of particles with a linear friction model as described in
`Hofmann et. al. 2025`_.
`Hofmann et al. 2025`_.

.. _Hofmann et. al. 2025: https://doi.org/10.48550/arXiv.2507.16388
.. _Hofmann et al. 2025: https://doi.org/10.48550/arXiv.2507.16388

The linear friction model is defined by the function

Expand Down Expand Up @@ -368,21 +408,21 @@ class FrictionLJLinear(FrictionalPair):
\Big]\, ,
\end{align*}

where :math:`\nu=(1/m_i+1/m_j)+(R_i^2/\mathbf{I}_i+R_j^2/\mathbf{I}_j))`, and
:math:`\mathbf{\xi}^\mathrm{f}_{ij}` and :math:`\mathbf{N}^\mathrm{f}_{ij}` are
where :math:`\nu=(1/m_i+1/m_j)+(R_i^2/I_i+R_j^2/I_j))`, and
:math:`\mathbf{\xi}_{ij}` and :math:`\mathbf{N}_{ij}` are
three dimensional Gaussian white noise vectors with correlations

.. math::

\begin{align*}
\langle \mathbf{\xi}^\mathrm{f}_{ij}(t) \mathbf{\xi}^\mathrm{f}_{kl}(t')
\rangle &= \mathbf{1}k_\mathrm{B}T
\langle \mathbf{\xi}_{ij} \mathbf{\xi}_{kl}
\rangle &= \mathbf{1}kT
(\delta_{ik}\delta_{jl}-\delta_{il}
\delta_{jk})\delta(t-t') \\
\langle \mathbf{N}^\mathrm{f}_{ij}(t) \mathbf{N}^\mathrm{f}_{kl}(t') \rangle &=
\mathbf{1}k_\mathrm{B}T(\delta_{ik}
\delta_{jk})/\delta t \\
\langle \mathbf{N}_{ij} \mathbf{N}_{kl} \rangle &=
\mathbf{1}kT(\delta_{ik}
\delta_{jl}+\delta_{il}\delta_{jk})
\delta(t-t')\, .
/\delta t\, .
\end{align*}

.. rubric:: Example:
Expand Down
Loading