# Towards Globally Consistent Stochastic Human Motion Prediction via Motion
Diffusion

Jiarui Sun, Girish Chowdhary
Dept. of Electrical and Computer Engineering
University of Illinois at Urbana, Champaign, USA
{jsun57<EMAIL_ADDRESS>

###### Abstract

Stochastic human motion prediction aims to predict multiple possible upcoming
pose sequences based on past human motion trajectories. Prior works focused
heavily on generating diverse motion samples, leading to inconsistent,
abnormal predictions from the immediate past observations. To address this
issue, in this work, we propose DiffMotion, a diffusion-based stochastic human
motion prediction framework that considers both the kinematic structure of the
human body and the globally temporally consistent nature of motion.
Specifically, DiffMotion consists of two modules: 1) a transformer-based
network for generating an initial motion reconstruction from corrupted motion,
and 2) a multi-stage graph convolutional network to iteratively refine the
generated motion based on past observations. Facilitated by the proposed
direct target prediction objective and the variance scheduler, our method is
capable of predicting accurate, realistic and consistent motion with an
appropriate level of diversity. Our results on benchmark datasets demonstrate
that DiffMotion outperforms previous methods by large margins in terms of
accuracy and fidelity while demonstrating superior robustness.

##  1 Introduction

Human Motion Prediction (HMP) aims to predict future human movements from
observed motion trajectories of the past. This task has a wide range of
applications in autonomous driving [1], robotics [2], animation creation [3]
and healthcare [4]. Many prior works address the HMP problem from the
deterministic point of view, trying to forecast one most likely future pose
sequence [5, 6, 7]. However, for complex scenarios where uncertainty must be
considered, deterministic HMP models may not be sufficient. For instance, in a
scenario where a robot needs to navigate a crowded street, a deterministic
model is unlikely to accurately predict the future movement of pedestrians and
cause the task to fail.

Recent works have focused on the stochastic HMP problem, addressing the
limitations of deterministic models by learning the conditional motion
distribution through generative models such as Generative Adversarial Networks
(GANs), Variational Autoencoders (VAEs), and Diffusion Models (DMs). Although
these methods can generate motion future stochastically, their effectiveness
is hampered by several issues. First, most stochastic HMP methods are hard to
train. Due to the complex nature of adversarial learning, GAN-based models
often require extensive tuning to stabilize training. While VAE-based methods
are relatively easy to optimize, most of them require multiple training stages
for better motion mode coverage and motion validity. Recent DM-based methods
also require additional training stages for either 1) output postprocessing or
2) motion encoding and decoding. Second, most recent stochastic HMP works are
diversity-centric. To generate diverse motion samples, these methods either
incorporate explicit diversity losses or construct additional sampling spaces
to avoid posterior collapse. While helping HMP models achieve better
performance on diversity-related metrics, such aggressive regularization
techniques often lead to predictions that are inconsistent and anomalous with
the immediate past observations. Although this issue is recently considered
and mitigated through modeling motion in the behavioral space, the problem of
predicting consistent and natural motion with respect to the provided history
remains.

To address these issues, we propose DiffMotion, a DM-based stochastic HMP
framework as shown in Figs. 1 and 2. DiffMotion can generate motion samples
with high quality and appropriate diversity without explicitly promoting
diversity through complex techniques such as dataset-level filtering. More
importantly, the motion generated by DiffMotion is realistic and consistent
with the motion history in both short-term (consistent with the recent past)
and long-term (consistent in terms of behavior). Furthermore, DiffMotion does
not require multiple training stages and is easy to optimize.

Our DiffMotion is novel in several aspects. First, instead of using the common
U-net [8] or a purely transformer-based backbone to model the reverse
diffusion process, our predictor consists of 1) a transformer-based module and
2) a multi-stage Graph Convolutional Network (GCN)-based module to learn in
the Discrete Cosine Transform (DCT) space. The transformer-based and GCN-based
modules together explicitly model the structural-temporal dependencies of
motion data, which is missing in prior DM-based works. Second, instead of
predicting the common noise target used in computer vision domain [9],
DiffMotion tries to directly predict motion as in [10]. This enables
DiffMotion to utilize losses that consider skeleton structure for supervision.
In addition, we propose to adjust the de facto cosine variance scheduler to
further improve sample diversity and accuracy.

Our main contributions are summarized as follows: 1) we propose a novel
diffusion-based framework for stochastic HMP that is able to generate more
coherent and realistic future motion compared to previous methods while
achieving an appropriate level of diversity; 2) we design a motion generator
that models structural-temporal dynamics of human motion explicitly at both
local and global levels. Our experiments on benchmark datasets demonstrate
that DiffMotion outperforms prior methods by significant margins on multiple
evaluation metrics, achieving state-of-the-art performance.

##  2 Related Work

Figure 1: The proposed DiffMotion predictor architecture. The motion predictor
first utilizes the transformer-based motion reconstruction module $F(\cdot)$
to produce an initial reconstruction $\tilde{y}_{t}$ from the $t^{th}$-level
corrupted motion $y_{t}$, then uses the GCN-based motion refinement module
$R(\cdot)$ to generate the entire motion sequence based on motion history $x$.
PE stands for Positional Encoding.

#### Human Motion Prediction

Pioneer data-driven methods on HMP relied on Markov models [11], binary latent
variables [12] or Gaussian processes [13] to capture the temporal dynamics of
human motion. These approaches can model the dynamics of simple motion such as
walking or golf swinging. However, they are inadequate in modeling complex
activities presented in modern large-scale datasets, such as boxing or playing
basketball. After the emergence of deep learning, modern HMP methods have been
relying on neural networks due to their superior effectiveness. The
deterministic HMP problem can be formulated as a regression task, since the
goal is to map from one sequence of historical poses to one predicted motion
trajectory. Due to the sequential nature of motion data, many works [14, 15,
16] employ Recurrent Neural Networks (RNNs) for temporal modeling. However,
these methods often suffer from problems such as error accumulation and
prediction discontinuity. To address these issues, feedforward methods are
proposed. Among these approaches, some employ Convolutional Neural Networks
(CNNs) [17, 18] or Transformers [19, 5] by treating motion as regular grid
data. Some other works leverage GCNs [5, 20] to model motion explicitly as a
spatio-temporal graph, achieving great results. Starting from the seminal work
[21], many feedforward methods have achieved great performance by transforming
motion trajectories into DCT coefficients to model human motion in the
frequency space. However, all deterministic HMP methods cannot model
multimodal data distributions and are thus not suitable for solving the
stochastic HMP problem.

To this end, deep generative methods are proposed to model the motion
distribution, and some early ones are GAN-based [22, 23]. However, the
adversarial nature of GANs makes them hard to train. In contrast, the
likelihood-based objective of VAE-based models is more interpretable and
easier to optimize, making these models [24, 25] popular. However, the recent
trend of focusing on motion diversity makes these methods either design
complex losses that often require dataset-level filtering, or involve multiple
training stages. More importantly, the diversity-centric objective causes
these methods to generate wild and inconsistent motion that are not aligned
with the observed history. Recently, few DM-based approaches are proposed for
stochastic HMP as diffusion-based methods have demonstrated great sample
quality and diversity in the computer vision domain. Wei _et al_. [26]
proposed MotionDiff, a two-stage learning framework that consists of a
transformer-based diffusion network for motion generation, and a GCN-based
network to enhance diversity as a postprocessing step. To address the problem
of history-future inconsistency, Barquero _et al_. [27] proposed to model the
reverse diffusion process in the behavioral space rather than the motion
coordinate space. To achieve this, an auxiliary VAE is designed for motion
encoding and decoding. Although these DM-based methods achieve better motion
diversity and prediction accuracy, they still require multiple training
stages. In addition, they always use RNN-based modules to encode the motion
history into fixed-size latent representations. We argue that such encoding
would inevitably damage the temporal consistency between the history motion
and future motion, making the predictions less accurate. We address this by
modeling the reverse diffusion process through a one-stage transformer- and
GCN-based motion generator that explicitly models structural and temporal
dependencies of motion.

#### Diffusion Model

Denoising Diffusion Probabilistic Models (DDPMs) [9] have recently received
significant attention and applied in many fields such as computer vision [28,
29] due to their superior sample quality and diversity. DDPMs define a
diffusion process in which random noise is gradually added to the data using a
Markov chain, and then learn how to reverse the process to generate desired
data samples. For the context of stochastic HMP that requires modeling
conditionals, prior works first use separate encoders to project the motion
history to intermediate representations and then employ either simple
concatenation [26] or the cross-attention mechanism [27] to incorporate
conditionals into noised latents. However, since these methods treat motion
history separately from the future and the motion history is always compressed
without structure, they fail to exploit the structure-temporal coherence of
motion. We address the problem by using a transformer-based architecture to
generate an initial reconstruction of the corrupted latent which is well-
aligned with the clean motion history.

One downside of DDPMs is their low efficiency. They typically require a large
number of denoising steps - often in the thousands - to generate a single
realistic data sample with high resolution, and this process cannot be
parallelized. Many efforts have been devoted to solving this problem,
including designing techniques for fast sampling [30, 31] and reducing data
resolution [32]. Inspired by [32], Barquero _et al_. proposed BeLFusion [27]
that generate motion samples through behavioral space modeling. Benefiting
from 1) learning motion instead of noise and 2) a motion predictor that
reduces the learning difficulty through frequency-space transformation and
direct motion conditioning, DiffMotion achieves state-of-the-art performance
on both prediction accuracy and fidelity with only a few diffusion steps,
while maintaining an appropriate level of sample diversity.

##  3 Our Approach

Figure 2: The sampling process of DiffMotion. At each step $t=T$ to $1$, the
motion predictor $G(\cdot)$ directly generate motion trajectory
$\hat{y}_{0}^{1:F}$, which is then diffused $t-1$ steps to repeat the process.

In this section, we present our framework DiffMotion, which is shown in Figs.
1 and 2. Following prior works [25], we are given motion history
$x^{1:H}=\\{x^{i}\\}_{i=1}^{H}$ of length $H$, and we aim to predict $F$
following poses $x^{H+1:H+F}=\\{x^{i}\\}_{i=H+1}^{H+F}$. In the stochastic HMP
setting, multiple future pose sequences are predicted given on historical
trajectory, and we use $y^{1:F}\coloneqq x^{H+1:H+F}$ to denote one such
predicted sample. The particular human pose at time step $i$ is denoted by
$x^{i}\in{\mathbb{R}}^{J\times 3}$ where $J$ is the number of body joints, as
we are interested in predicting 3D joint positions.

###  3.1 Conditional Motion Diffusion

Let $\\{x_{t}\\}_{t=0}^{T}$ denote a general Markov noising process where
$x_{0}$ represents the true data samples. In the unconditional case, the
forward diffusion transitions are denoted as:

$q(x_{t}|x_{t-1})=\mathcal{N}(x_{t};\sqrt{\alpha_{t}}x_{t-1},(1-\alpha_{t})\textrm{I}),$
(1)

where $\\{\alpha_{t}\\}_{t=1}^{T}\in(0,1)$ control the noise level, and can
either be fixed [9] or learned [33]. In order to estimate the true data
distribution, the parameterized reverse diffusion process is constructed to
progressively denoise the corrupted data samples $x_{T}$ from $t=T$ to $t=0$:

$p_{\theta}(x_{t-1}|x_{t})=\mathcal{N}(x_{t-1};\mu_{\theta}(x_{t},t),\sigma_{\theta}^{2}(x_{t},t)\textrm{I}).$
(2)

Ho _et al_. [9] proposed to 1) fix $\sigma_{\theta}^{2}(x_{t},t)$, 2) predict
the perturbation $\epsilon$ through reparameterization, and 3) optimize the
model with a simple objective
$L(\theta)=\mathbb{E}_{t}[\|\epsilon-\epsilon_{\theta}(x_{t},t)\|^{2}]$. These
techniques are quickly adopted by later DM-based works in various fields due
to the great performance they offer.

In our HMP context, we need to extend the above formulation to the conditional
case. Specifically, the reverse diffusion transition becomes:

$p_{\theta}(y_{t-1}|y_{t},x)=\mathcal{N}(y_{t-1};\mu_{\theta}(y_{t},x,t),\sigma_{\theta}^{2}(y_{t},x,t)\textrm{I}),$
(3)

where $x$ represents the motion history and $y_{t}$ represents $t$-step noised
future motion.

###  3.2 Motion Diffusion Pipeline

With the above diffusion formulation under the HMP setting, the popular setup
from Ho _et al_. has several disadvantages. First, predicting the noise
prevents models from exploiting different geometric losses that have been
extensively studied in HMP research, since the noise does not contain
geometry-related information. Second, different from works from other fields
such as [32] where the conditional signal and data are always multimodal, the
conditional signal $x$ comes from the same coordinate space as $y_{0}$ that
they directly form a complete motion trajectory seamlessly. This demonstrates
the potential of using a unified framework to perform stochastic HMP, where
previous approaches have failed. Lastly, only treating $x$ as conditionals
ignores the temporal consistency between history and future: the entire motion
trajectory should be globally consistent in terms of behavior in most cases;
it also needs to be locally smooth and avoid strong discontinuities close to
the prediction horizon.

#### Prediction target

To address the issues, we first propose to predict the future motion directly,
which is the same objective used in VAE-based works [24, 25]. Specifically,
instead of predicting $\epsilon$, we choose another reparameterization such
that $\hat{y}_{0}\leftarrow G_{\theta}(y_{t},x,t)$, where $\hat{y}_{0}$ is the
learned approximation of target future motion $y_{0}$. After we obtain
$\hat{y}_{0}$, such predicted motion future is diffused $t-1$ steps, and used
as input to generate the next prediction on the sampling chain. The entire
sampling process is illustrated in Fig. 2. The forward diffusion process can
be simplified as:

$q(y_{t-1}|y_{0})=\mathcal{N}(y_{t-1};\sqrt{\bar{\alpha}_{t-1}}y_{0},(1-\bar{\alpha}_{t-1})\textrm{I}),$
(4)

where $\bar{\alpha}_{t-1}=\prod_{s=1}^{t-1}\alpha_{s}$. Importantly,
predicting $y_{0}$ directly allow us to use structure-aware losses.

#### Generator architecture

With the $y_{0}$-prediction objective, our motion generator can directly
benefit from works that address the deterministic HMP problem. Specifically,
$G_{\theta}(\cdot)$ can leverage the spatio-temporal graph structure of motion
data. Our designed network, as illustrated in Fig. 1, consists of two parts: a
time-dependent transformer-based module $F(\cdot)$ for denoising the noised
future motion, and a multi-stage GCN-based refinement module $R(\cdot)$ that
considers the conditional motion history.

Inspired by [10], we employ a transformer-based architecture $F(y_{t},t)$ for
generating an “initial reconstruction” $\tilde{y}_{0}$ from the corrupted
future motion $y_{t}$. Based on the transformer encoder, $F(\cdot)$ can
explicit model temporal correlations across noised frames. The time step $t$
is first mapped to an embedding and then projected to the same latent
dimension as the motion frames $\\{y_{t}^{i}\\}_{i=1}^{F}$ using a feedforward
network to produce the time token $z_{t}$. The motion frames
$\\{y_{t}^{i}\\}_{i=1}^{F}$ are projected to the transformer latent dimension
and then summed with positional encoding [34] to retain positional
information. The transformed motion frames are then prepended with the time
token $z_{t}$ and fed to the transformer encoder. To get $\tilde{y}_{0}$ from
$F(\cdot)$, we first exclude the first output token corresponding to $t$, and
then project the learned representations back to the pose dimension $J\times
3$.

For the multi-stage GCN-based refinement module $R(x,\tilde{y}_{0})$, we
simply adopt the state-of-the-art architecture [35] from the deterministic HMP
field due to its superior performance. $R(\cdot)$ first takes as input the
concatenation of $x$ and $\tilde{y}_{0}$, and then iteratively refines the
entire motion trajectory by converting representations between the pose space
and the frequency space. Specifically, the motion trajectory is first
transformed into DCT coefficients at the start of each learning stage. After
being fed into multiple GCN layers to capture structural dependencies between
joints, the latent motion representation is projected back to the pose space,
serving as input for the next refinement round. The learned future motion
$\hat{y}_{0}$ is finally obtained by excluding the portion corresponding to
the motion history. That is:

$\hat{y}_{0}\leftarrow[\hat{x};\hat{y}_{0}]=R(x,F(y_{t},t))\coloneqq
G_{\theta}(y_{t},x,t).$ (5)

It is important to note that DiffMotion only employs one neural network
$G_{\theta}(\cdot)$, different from prior DM-based stochastic HMP methods [26,
27]. In addition, instead of only generating the predicted future motion
$\hat{y}_{0}$, our motion generator also produces the history estimation
$\hat{x}$. This is particularly beneficial for generating globally consistent
motion trajectories as previous methods show [21, 25].

#### Variance scheduler

The variance scheduler $\\{1-\alpha_{t}\\}_{t=1}^{T}$ is essential to the
performance of DMs. The linear scheduler [9] and the cosine scheduler [36] are
commonly used due to their simplicity and performance they provide. However,
we argue whether these schedulers are optimal for our case. As mentioned, our
generator $G_{\theta}(y_{t},x,t)$ tries to 1) directly predicts motion, and 2)
directly takes $x^{1:H}$ into account, which is different from the multimodal
conditionals used in DM-based works from other fields with its clean temporal
behavioral guidance.

To this end, we study how the standard schedulers are designed. Eq. 4 is
equivalent to:

$y_{t}=\sqrt{\bar{\alpha}_{t}}y_{0}+\sqrt{1-\bar{\alpha}_{t}}\epsilon,$ (6)

where $\epsilon\sim{\mathcal{N}}(0,\textrm{I})$. In the theoretical setup,
$\bar{\alpha}_{0}=1$ and $\bar{\alpha}_{T}=0$ such that the clean data is
completely corrupted to standard Gaussian noise. Governed by this principle,
both linear scheduler and cosine scheduler are empirically designed to have
$\bar{\alpha}_{0}\approx 1$.

However, we find that the results obtained by using standard schedulers often
indicate a lower sample diversity and accuracy comparing to VAE-based methods.
We hypothesize that having $\bar{\alpha}_{0}$ close to $1$ prevents
$G_{\theta}(y_{t},x,t)$ from producing diverse and accurate samples.
Specifically, since we feed the clean motion history directly into
$G_{\theta}(\cdot)$ and aim to predict $y_{0}$ by modeling $x^{1:H}$
explicitly, the strong guidance from $x^{1:H}$ and the strong regularization
of $G_{\theta}(\cdot)$ makes the prediction task overly simple as the reverse
process progresses. Due to the poor mode coverage of the current HMP datasets
with their limited size, the predictions may not be multimodal, resulting in
poor sample diversity and accuracy.

To address the issue, instead of designing complex techniques to promote
sample diversity, we simply adjust the original cosine scheduler such that it
is suitable for our $y_{0}$-prediction HMP setting. Specifically, we have:

$\bar{\alpha}_{t}=\frac{f(t)}{f(0)},\quad
f(t)=\cos\left(\frac{t/T+1}{2}\cdot\frac{\pi}{2}\right)^{2}.$ (7)

Essentially, we use an offset of $1$ for the cosine scheduler instead of the
default offset $0.008$. Doing so, we have $\bar{\alpha}_{0}=0.5$, which makes
strongly guided $y_{0}$-prediction not a trivial task when $t$ is close to
$0$. From another perspective, we purposefully make sure that the motion
history is noised to a certain level across all diffusion steps such that
$G_{\theta}(\cdot)$ always need to learn the structural-temporal dependencies
from the data, instead of relying heavily on motion history guidance.

###  3.3 Learning Algorithm

As previously mentioned, predicting future motion directly (or aiming to
reconstruct the entire motion trajectory) enables us to utilize geometric
losses to supervise the model $G_{\theta}(\cdot)$ for better motion
consistency. As such, we utilize a structure-aware loss [35] to take into
account the details of human motion structure. The structure-aware loss
weights all joints differently to reflect varying degrees of importance in the
motion context. That is:

${\mathcal{L}}_{\theta}(G,y_{0},x)=\mathbb{E}_{\begin{subarray}{c}y_{0}\sim
q(y_{0}\mid x),\\\
t\sim[1,T]\end{subarray}}({\mathcal{L}}_{\textrm{rec}}(G_{\theta}(y_{t},x,t),y_{0},x)),$
(8)

that

${\mathcal{L}}_{\textrm{rec}}=\frac{1}{J}\sum_{j=1}^{J}(\|(x^{j}-\hat{x}^{j})\cdot\lambda^{j}\|_{1}+\|(y_{0}^{j}-\hat{y}_{0}^{j})\cdot\lambda^{j}\|_{1}),$
(9)

where the superscript $j$ denotes joint index, and $\lambda^{j}$ represents
the specific assigned weight. The weights are computed based on the specifics
of the kinematic structure of the human body, and they do not require
learning. In particular, from the structural point of view, higher weights are
assigned to joints that can exhibit vigorous motions (external joints).

Additionally, we use the relaxation technique [25, 27] to further promote the
prediction diversity. Specifically, we sample $k$ future motions for each
historical motion trajectory and only optimize ${\mathcal{L}}(\theta)$ towards
the most accurate prediction. That is:

${\mathcal{L}}_{\textrm{final}}=\min_{k}{\mathcal{L}}_{\theta}(G^{k},y_{0},x),$
(10)

where $G^{k}$ is the $k^{th}$ reconstructed motion trajectory based on the
original pose sequence $[x;y_{0}]$.

##  4 Experiments

Table 1: Quantitative results for Human3.6M and AMASS datasets. The best
results are highlighted in bold, and the second best results are underlined.
FID is only supported by Human3.6M as AMASS cannot be separated by classes.

| Human3.6M | AMASS
---|---|---
| APD $\uparrow$ | APDE $\downarrow$ | ADE $\downarrow$ | FDE $\downarrow$ | MMADE $\downarrow$ | MMFDE $\downarrow$ | CMD $\downarrow$ | FID $\downarrow$ | APD $\uparrow$ | APDE $\downarrow$ | ADE $\downarrow$ | FDE $\downarrow$ | MMADE $\downarrow$ | MMFDE $\downarrow$ | CMD $\downarrow$
Zero-Velocity | 0.000 | 8.079 | 0.597 | 0.884 | 0.683 | 0.909 | 22.812 | 0.606 | 0.000 | 9.292 | 0.755 | 0.992 | 0.814 | 1.015 | 39.262
HP-GAN [22] | 7.214 | - | 0.858 | 0.867 | 0.847 | 0.858 | - | - | - | - | - | - | - | - | -
DeLiGAN [37] | 6.509 | - | 0.483 | 0.534 | 0.520 | 0.545 | - | - | - | - | - | - | - | - | -
GMVAE [38] | 6.769 | - | 0.461 | 0.555 | 0.524 | 0.566 | - | - | - | - | - | - | - | - | -
TPK [39] | 6.723 | 1.906 | 0.461 | 0.560 | 0.522 | 0.569 | 6.326 | 0.538 | 9.283 | 2.265 | 0.656 | 0.675 | 0.658 | 0.674 | 17.127
MT-VAE [40] | 0.403 | - | 0.457 | 0.595 | 0.716 | 0.883 | - | - | - | - | - | - | - | - | -
BoM [41] | 6.265 | - | 0.448 | 0.533 | 0.514 | 0.544 | - | - | - | - | - | - | - | - | -
Motron [42] | 7.168 | 2.583 | 0.375 | 0.488 | 0.509 | 0.539 | 40.796 | 13.743 | - | - | - | - | - | - | -
DSF [43] | 9.330 | - | 0.493 | 0.592 | 0.550 | 0.599 | - | - | - | - | - | - | - | - | -
DLow [24] | 11.741 | 3.781 | 0.425 | 0.518 | 0.495 | 0.531 | 4.927 | 1.255 | 13.170 | 4.243 | 0.590 | 0.612 | 0.618 | 0.617 | 15.185
MultiObj [44] | 14.240 | - | 0.414 | 0.516 | - | - | - | - | - | - | - | - | - | - | -
GSPS [25] | 14.757 | 6.749 | 0.389 | 0.496 | 0.476 | 0.525 | 10.758 | 2.103 | 12.465 | 4.678 | 0.563 | 0.613 | 0.609 | 0.633 | 18.404
DivSamp [45] | 15.310 | 7.479 | 0.370 | 0.485 | 0.475 | 0.516 | 11.692 | 2.083 | 24.724 | 15.837 | 0.564 | 0.647 | 0.623 | 0.667 | 50.239
BeLFusion [27] | 7.602 | 1.662 | 0.372 | 0.474 | 0.473 | 0.507 | 5.988 | 0.209 | 9.376 | 1.977 | 0.513 | 0.560 | 0.569 | 0.585 | 16.995
Ours | 7.753 | 1.570 | 0.351 | 0.458 | 0.495 | 0.507 | 3.279 | 0.106 | 10.320 | 2.117 | 0.490 | 0.543 | 0.467 | 0.465 | 10.502

In this section, we evaluate the effectiveness of DiffMotion. Instead of
merely focusing on the diversity aspect, we aim to study how DiffMotion
performs in terms of generating diverse and realistic motion samples that are
well aligned with the provided history. We also seek to examine DiffMotion’s
generalization capabilities.

For this purpose, we first present the datasets we used in Section 4.1. Then,
we introduce our experimental setup in Section 4.2, including evaluation
metrics and the implementation details of our model. Next, we assess
DiffMotion both quantitatively and qualitatively through comparing it with
state-of-the-art methods in Section 4.3. Lastly, a detailed ablation analysis
is carried out and presented in Section 4.4 to evaluate the effectiveness of
the proposed model components.

###  4.1 Datasets

Human3.6M [46], the most widely used dataset for stochastic HMP, contains
motion clips of 7 subjects performing 15 different actions recorded at 50 Hz.
For a fair comparison with most prior works [25, 27], we adopt the evaluation
protocol of [24], where a 16-joint skeleton is used for human structure
modeling. Specifically, we aim to predict 2s (100 frames) into the future
based on 0.5s (25 frames) of observation.

AMASS [47] unifies multiple Mocap datasets, such as the popular HumanEva-I
dataset [48] using shared SMPL [49] parameterization for human skeleton
modeling. Compared with HumanEva-I which is often used in prior works, AMASS
is much larger, and can be used to perform cross-dataset evaluation [27] to
examine a model’s generalization ability. To this end, we aim to predict 2s
(30 frames) into the future based on 0.5s (120 frames) of observation after
downsampling. Other evaluation settings such as frame-rate,
train/validate/test partition are all set to be the same as previous work [27]
for fair comparisons.

###  4.2 Experimental Setup

Evaluation metrics. To evaluate stochastic HMP models quantitatively, most
prior works adopt the metrics proposed in [24] to measure sample diversity and
accuracy. The Average Pairwise Distance (APD) computes the averaged $\ell_{2}$
distance between all generated sample pairs given one motion history,
measuring the sample diversity. The Average and the Final Displacement Error
(ADE and FDE) calculate the averaged all-time and the last-frame $\ell_{2}$
distances respectively between the groundtruth and the closest prediction,
measuring the sample accuracy. In addition, the multimodal versions of ADE and
FDE, MMADE, and MMFDE assess a method’s ability to produce multimodal
predictions. To do so, MMADE and MMFDE first obtain multimodal groundtruth
future motion through grouping similar last observed poses $x^{H}$ in terms of
$\ell_{2}$ distance, and then compare all predictions with the multimodal
groundtruth.

While these metrics can measure a method’s ability to generate diverse and
accurate samples to some extent, several issues are evident. First, while the
APD metric captures the diversity level of the predicted motion, it fails to
consider consistency between the input motion history and future motion
prediction. As previously mentioned, prior methods that focus heavily on the
diversity aspect (_i.e_., aim to obtain high APD) often generate motion with
abnormal sudden changes, which cannot be reflected by these metrics. Second,
the MMADE and MMFDE metrics are computed based on the multimodal groundtruth,
which is collected based solely on $x^{H}$. However, motion of very distinct
actions or behaviors can contain similar intermediate poses. Hence, these
multimodal metrics can not accurately reflect a method’s ability to generate
consistent motion.

To address these issues and better quantify the realism of motion, the Fréchet
Inception Distance (FID) metric is used by [45, 50, 51, 52] to assess the
similarity between the distributions of generated and real motion based on the
statistics of feature representations extracted by a pre-trained network.
Recently, [27] proposed several additional metrics to better quantify a
model’s capacity of producing globally behavior-consistent motion. The area of
the Cumulative Motion Distribution (CMD) metric measures the difference
between the areas under the cumulative true motion and predicted motion
distributions, measuring the plausibility of predicted motion at a global
level. Additionally, aiming to analyze to which extent the diversity is
properly modeled, the Average Pairwise Distance Error (APDE) metric is defined
as the absolute error between the APD of the multimodal groundtruth and the
APD of the predictions.

Implementation details. We implement the transformer-based module $F(y_{t},t)$
using 8 transformer encoder layers with latent dimension 512 and dropout ratio
0.1. Both the time step $t$ embedding and the positional encoding are
sinusoidal. We use a 2-layer MLP to project the time step embedding to the
transformer dimension. The GCN-based refinement module $R(x,\tilde{y}_{0})$
consists of 3 learning stages, each of which contains 2 GCN-based learning
blocks. The latent dimension and the dropout ratio of the refinement module
are 256 and 0.5.

Adam [53] is used as the optimizer for all experiments, with 0.0001 as the
initial learning rate. For both datasets, the model is trained for 500 epochs
with a batch size of 64, and the learning rate is multiplied by $0.99$ after
the 200th epoch. We use NVIDIA V100 GPUs for all experiments.

###  4.3 Results

In this section, we compare DiffMotion with current state-of-the-art
stochastic HMP methods on Human3.6M and AMASS. Besides the zero-velocity
baseline, we characterize these methods into three categories: 1) general, 2)
diversity-centric and 3) behavior-driven. Different categories are separated
by lines in Table 1.

Figure 3: Qualitative results of DiffMotion compared with state-of-the-art
methods. DiffMotion can consistently generate context-aware motion while
maintaining a good level of diversity.

We first study DiffMotion’s modeling ability quantitatively. As shown in Table
1, DiffMotion significantly outperforms previous methods on accuracy metrics
(ADE and FDE) by large margins on both datasets. More importantly, our method
is capable of generating much more globally consistent and realistic future
motion compared to other methods as it excels at realism metrics (CMD and
FID). On the contrary, methods such as GSPS and DivSamp exhibit high CMD and
FID values as these methods explicitly promote diversity. As a result, they
tend to generate diverse but behaviourally inconsistent motion. Besides, our
method also exhibits a low APDE value, indicating that DiffMotion can properly
model the stochasticity level of future motion given the past. In terms of
diversity (APD), we can observe that our model performs better than the other
DM-based model BeLFusion. While falling short in MMADE on Human3.6M compared
with the most recent approaches, we argue that 1) the multimodal metrics
cannot accurately reflect the ability of the model to generate consistent
behavioral samples as the collected multimodal groudtruth can be wildly
divergent, and 2) the dataset is limited in size thus can be not comprehensive
compared with AMASS.

We also perform a qualitative analysis to further examine our model’s
generative power. Specifically, in Fig. 3, we compare DiffMotion against
multiple state-of-the-art methods on both Human3.6M and AMASS. Two actions
(Directions and Eating) are selected from Human3.6M, and two sub-datasets
(DanceDB and GRAB) are used from AMASS. From the visualization, we can first
confirm that DiffMotion is able to generate accurate predictions given the
motion history. In particular, the first few poses predicted by DiffMotion are
much more aligned with the groundtruth, while others often exhibit sudden
motion discontinuities. Second, we can observe that DiffMotion can produce
diverse predictions that are adapted to the context. Specifically, the motion
predicted for the eating action is characterized by diverse arm movements,
while the legs remain stationary. The qualitative results further confirm the
advantage of DiffMotion over prior methods.

###  4.4 Ablation study

In this section, We conduct an ablation analysis on the Human3.6M dataset to
investigate how the proposed components of DiffMotion may affect its motion
modeling ability.

General model components. We first generally show how each design choice
improves the effectiveness of DiffMotion, including 1) the transformer-based
motion reconstruction module $F(\cdot)$, 2) the multi-stage GCN-based
refinement module $R(\cdot)$ and 3) the proposed variance scheduler.
Specifically, We first establish a plain baseline model, which utilizes the
popular cosine variance scheduler instead of the proposed scheduler.
Additionally, the baseline model does not utilize the reconstruction module
and only tries to learn $y_{0}$ through a single-stage GCN-based module,
similar to [21]. Then, we improve upon the baseline by adding the motion
reconstruction module $F(\cdot)$. Next, we incorporate multi-stage refinement
by using the full refinement module $R(\cdot)$ with 3 refinement stages.
Lastly, we examine the full model with the proposed variance scheduler instead
of the standard cosine scheduler.

All ablation results are summarized in Table 2. First, we can observe that the
motion reconstruction module is indispensable for DiffMotion. The initial
reconstruction provides crucial guidance for the model to produce accurate and
consistent predictions that closely adhere to motion history. Without
$F(\cdot)$, the model exhibits a remarkably high CMD value, indicating that
its predictions usually contain large displacements between consecutive poses.
Second, the multi-stage refinement module boosts the prediction accuracy and
the consistency between motion history and predicted motion greatly. Lastly,
we can see that the proposed variance scheduler brings the ADE, FDE and FID
values to another level, which maintaining a good level of sample diversity.
In general, all components of DiffMotion contribute to it performance.

$F(y_{t},t)$ | $R(x,\tilde{y}_{0})$ | sched. | APD | APDE | ADE | FDE | CMD | FID
---|---|---|---|---|---|---|---|---
✗ | ✗ | ✗ | 13.571 | 5.522 | 1.000 | 1.033 | 1044.016 | 1.070
✓ | ✗ | ✗ | 7.612 | 1.402 | 0.407 | 0.500 | 2.871 | 0.401
✓ | ✓ | ✗ | 7.530 | 1.473 | 0.380 | 0.486 | 3.041 | 0.297
✓ | ✓ | ✓ | 7.753 | 1.570 | 0.351 | 0.458 | 3.279 | 0.106
Table 2: Ablation on DiffMotion’s general architecture. $k$ | APD | ADE | FDE | CMD | FID
---|---|---|---|---|---
1 | 4.104 | 0.345 | 0.474 | 5.930 | 0.159
2 (ours) | 7.753 | 0.351 | 0.458 | 3.279 | 0.106
3 | 9.336 | 0.374 | 0.482 | 4.971 | 0.218
Table 3: Ablation on implicit diversity relaxation.

Implicit diversity relaxation. Then, we perform a small ablation study on the
implicit diversity relaxation technique. Specifically, as shown in Table 3, we
train DiffMotion with $k=1,2,3$, which is the relaxation hyperparameter
described in Eq. 10. Unsurprisingly, as $k$ increases, we observe a rising APD
value. In addition, for most of the metrics, $k=2$ provides us with the best
results. This indicates that our DiffMotion framework can indeed regularize
its predictions to some extent, as our method does not utilize additional
techniques like pose priors to supervise out-of-distribution predictions.
However, a strong supervision is still necessary.

##  5 Conclusion

In this work, we present DiffMotion, a novel framework for stochastic HMP.
Instead of explicitly focusing on generating diversified motions like prior
works, we aim to produce natural, globally behaviourally-consistent, and
properly diverse human motions. Specifically, we propose 1) a motion
reconstruction module that generates initial reconstructed motion from
corrupted pose sequences, and 2) a multi-stage refinement module that
generates future motion in a progressive manner using the observed motion
history and denoised future motion. Importantly, We aim to directly learn
future motion under the diffusion framework which is better suited for the
stochastic HMP problem, and we propose to adjust the cosine variant scheduler
for better performance. The results obtained from extensive experiments and
detailed analysis show that DiffMotion has significant performance gains over
state-of-the-art baselines on benchmark datasets, demonstrating the
effectiveness of our method.

## References

  * [1] B. Paden, M. Cáp, S. Z. Yong, D. S. Yershov, and E. Frazzoli, “A survey of motion planning and control techniques for self-driving urban vehicles,” IEEE Trans. Intell. Veh., vol. 1, no. 1, pp. 33–55, 2016.
  * [2] L. Gui, K. Zhang, Y. Wang, X. Liang, J. M. F. Moura, and M. Veloso, “Teaching robots to predict human motion,” in 2018 IEEE/RSJ International Conference on Intelligent Robots and Systems, IROS 2018, Madrid, Spain, October 1-5, 2018, pp. 562–567, 2018.
  * [3] H. van Welbergen, B. J. H. van Basten, A. Egges, Z. Ruttkay, and M. H. Overmars, “Real time animation of virtual humans: A trade-off between naturalness and control,” Comput. Graph. Forum, vol. 29, no. 8, pp. 2530–2554, 2010.
  * [4] W. Taylor, S. A. Shah, K. Dashtipour, A. Zahid, Q. H. Abbasi, and M. A. Imran, “An intelligent non-invasive real-time human activity recognition system for next-generation healthcare,” Sensors, vol. 20, no. 9, p. 2653, 2020.
  * [5] W. Mao, M. Liu, and M. Salzmann, “History repeats itself: Human motion prediction via motion attention,” in Computer Vision - ECCV 2020 - 16th European Conference, Glasgow, UK, August 23-28, 2020, Proceedings, Part XIV, pp. 474–489, 2020.
  * [6] T. Ma, Y. Nie, C. Long, Q. Zhang, and G. Li, “Progressively generating better initial guesses towards next stages for high-quality human motion prediction,” in IEEE/CVF Conference on Computer Vision and Pattern Recognition, CVPR 2022, New Orleans, LA, USA, June 18-24, 2022, pp. 6427–6436, 2022.
  * [7] L. Dang, Y. Nie, C. Long, Q. Zhang, and G. Li, “MSR-GCN: multi-scale residual graph convolution networks for human motion prediction,” in 2021 IEEE/CVF International Conference on Computer Vision, ICCV 2021, Montreal, QC, Canada, October 10-17, 2021, pp. 11447–11456, 2021.
  * [8] O. Ronneberger, P. Fischer, and T. Brox, “U-net: Convolutional networks for biomedical image segmentation,” in Medical Image Computing and Computer-Assisted Intervention - MICCAI 2015 - 18th International Conference Munich, Germany, October 5 - 9, 2015, Proceedings, Part III, pp. 234–241, 2015.
  * [9] J. Ho, A. Jain, and P. Abbeel, “Denoising diffusion probabilistic models,” in Advances in Neural Information Processing Systems 33: Annual Conference on Neural Information Processing Systems 2020, NeurIPS 2020, December 6-12, 2020, virtual, 2020.
  * [10] G. Tevet, S. Raab, B. Gordon, Y. Shafir, D. Cohen-Or, and A. H. Bermano, “Human motion diffusion model,” CoRR, vol. abs/2209.14916, 2022.
  * [11] A. M. Lehrmann, P. V. Gehler, and S. Nowozin, “Efficient nonlinear markov models for human motion,” in 2014 IEEE Conference on Computer Vision and Pattern Recognition, CVPR 2014, Columbus, OH, USA, June 23-28, 2014, pp. 1314–1321, IEEE Computer Society, 2014.
  * [12] G. W. Taylor, G. E. Hinton, and S. T. Roweis, “Modeling human motion using binary latent variables,” in Advances in Neural Information Processing Systems 19, Proceedings of the Twentieth Annual Conference on Neural Information Processing Systems, Vancouver, British Columbia, Canada, December 4-7, 2006 (B. Schölkopf, J. C. Platt, and T. Hofmann, eds.), pp. 1345–1352, MIT Press, 2006.
  * [13] J. M. Wang, D. J. Fleet, and A. Hertzmann, “Gaussian process dynamical models for human motion,” IEEE Trans. Pattern Anal. Mach. Intell., vol. 30, no. 2, pp. 283–298, 2008.
  * [14] J. Martinez, M. J. Black, and J. Romero, “On human motion prediction using recurrent neural networks,” in 2017 IEEE Conference on Computer Vision and Pattern Recognition, CVPR 2017, Honolulu, HI, USA, July 21-26, 2017, pp. 4674–4683, 2017.
  * [15] K. Fragkiadaki, S. Levine, P. Felsen, and J. Malik, “Recurrent network models for human dynamics,” in 2015 IEEE International Conference on Computer Vision, ICCV 2015, Santiago, Chile, December 7-13, 2015, pp. 4346–4354, 2015.
  * [16] A. Jain, A. R. Zamir, S. Savarese, and A. Saxena, “Structural-rnn: Deep learning on spatio-temporal graphs,” in 2016 IEEE Conference on Computer Vision and Pattern Recognition, CVPR 2016, Las Vegas, NV, USA, June 27-30, 2016, pp. 5308–5317, 2016.
  * [17] C. Li, Z. Zhang, W. S. Lee, and G. H. Lee, “Convolutional sequence to sequence model for human dynamics,” in 2018 IEEE Conference on Computer Vision and Pattern Recognition, CVPR 2018, Salt Lake City, UT, USA, June 18-22, 2018, pp. 5226–5234, 2018.
  * [18] J. Bütepage, M. J. Black, D. Kragic, and H. Kjellström, “Deep representation learning for human motion prediction and classification,” in 2017 IEEE Conference on Computer Vision and Pattern Recognition, CVPR 2017, Honolulu, HI, USA, July 21-26, 2017, pp. 1591–1599, 2017.
  * [19] E. Aksan, M. Kaufmann, P. Cao, and O. Hilliges, “A spatio-temporal transformer for 3d human motion prediction,” in International Conference on 3D Vision, 3DV 2021, London, United Kingdom, December 1-3, 2021, pp. 565–574, 2021\.
  * [20] W. Mao, M. Liu, M. Salzmann, and H. Li, “Multi-level motion attention for human motion prediction,” Int. J. Comput. Vis., vol. 129, no. 9, pp. 2513–2535, 2021.
  * [21] W. Mao, M. Liu, M. Salzmann, and H. Li, “Learning trajectory dependencies for human motion prediction,” in 2019 IEEE/CVF International Conference on Computer Vision, ICCV 2019, Seoul, Korea (South), October 27 - November 2, 2019, pp. 9488–9496, 2019.
  * [22] E. Barsoum, J. Kender, and Z. Liu, “HP-GAN: probabilistic 3d human motion prediction via GAN,” in 2018 IEEE Conference on Computer Vision and Pattern Recognition Workshops, CVPR Workshops 2018, Salt Lake City, UT, USA, June 18-22, 2018, pp. 1418–1427, 2018.
  * [23] J. N. Kundu, M. Gor, and R. V. Babu, “Bihmp-gan: Bidirectional 3d human motion prediction GAN,” in The Thirty-Third AAAI Conference on Artificial Intelligence, AAAI 2019, The Thirty-First Innovative Applications of Artificial Intelligence Conference, IAAI 2019, The Ninth AAAI Symposium on Educational Advances in Artificial Intelligence, EAAI 2019, Honolulu, Hawaii, USA, January 27 - February 1, 2019, pp. 8553–8560, 2019.
  * [24] Y. Yuan and K. Kitani, “Dlow: Diversifying latent flows for diverse human motion prediction,” in Computer Vision - ECCV 2020 - 16th European Conference, Glasgow, UK, August 23-28, 2020, Proceedings, Part IX, pp. 346–364, 2020.
  * [25] W. Mao, M. Liu, and M. Salzmann, “Generating smooth pose sequences for diverse human motion prediction,” in 2021 IEEE/CVF International Conference on Computer Vision, ICCV 2021, Montreal, QC, Canada, October 10-17, 2021, pp. 13289–13298, 2021.
  * [26] D. Wei, H. Sun, B. Li, J. Lu, W. Li, X. Sun, and S. Hu, “Human joint kinematics diffusion-refinement for stochastic motion prediction,” CoRR, vol. abs/2210.05976, 2022.
  * [27] G. Barquero, S. Escalera, and C. Palmero, “Belfusion: Latent diffusion for behavior-driven human motion prediction,” CoRR, vol. abs/2211.14304, 2022\.
  * [28] P. Dhariwal and A. Q. Nichol, “Diffusion models beat gans on image synthesis,” in Advances in Neural Information Processing Systems 34: Annual Conference on Neural Information Processing Systems 2021, NeurIPS 2021, December 6-14, 2021, virtual, pp. 8780–8794, 2021.
  * [29] F. Croitoru, V. Hondru, R. T. Ionescu, and M. Shah, “Diffusion models in vision: A survey,” CoRR, vol. abs/2209.04747, 2022.
  * [30] J. Song, C. Meng, and S. Ermon, “Denoising diffusion implicit models,” in 9th International Conference on Learning Representations, ICLR 2021, Virtual Event, Austria, May 3-7, 2021, 2021.
  * [31] Q. Zhang and Y. Chen, “Fast sampling of diffusion models with exponential integrator,” CoRR, vol. abs/2204.13902, 2022.
  * [32] R. Rombach, A. Blattmann, D. Lorenz, P. Esser, and B. Ommer, “High-resolution image synthesis with latent diffusion models,” in IEEE/CVF Conference on Computer Vision and Pattern Recognition, CVPR 2022, New Orleans, LA, USA, June 18-24, 2022, pp. 10674–10685, 2022.
  * [33] D. P. Kingma, T. Salimans, B. Poole, and J. Ho, “Variational diffusion models,” CoRR, vol. abs/2107.00630, 2021.
  * [34] A. Vaswani, N. Shazeer, N. Parmar, J. Uszkoreit, L. Jones, A. N. Gomez, L. Kaiser, and I. Polosukhin, “Attention is all you need,” in Advances in Neural Information Processing Systems 30: Annual Conference on Neural Information Processing Systems 2017, December 4-9, 2017, Long Beach, CA, USA, pp. 5998–6008, 2017.
  * [35] J. Sun and G. Chowdhary, “Towards accurate human motion prediction via iterative refinement,” 2023.
  * [36] A. Q. Nichol and P. Dhariwal, “Improved denoising diffusion probabilistic models,” in Proceedings of the 38th International Conference on Machine Learning, ICML 2021, 18-24 July 2021, Virtual Event, pp. 8162–8171, 2021.
  * [37] S. Gurumurthy, R. K. Sarvadevabhatla, and R. V. Babu, “Deligan: Generative adversarial networks for diverse and limited data,” in 2017 IEEE Conference on Computer Vision and Pattern Recognition, CVPR 2017, Honolulu, HI, USA, July 21-26, 2017, pp. 4941–4949, 2017.
  * [38] N. Dilokthanakul, P. A. M. Mediano, M. Garnelo, M. C. H. Lee, H. Salimbeni, K. Arulkumaran, and M. Shanahan, “Deep unsupervised clustering with gaussian mixture variational autoencoders,” CoRR, vol. abs/1611.02648, 2016.
  * [39] J. Walker, K. Marino, A. Gupta, and M. Hebert, “The pose knows: Video forecasting by generating pose futures,” in IEEE International Conference on Computer Vision, ICCV 2017, Venice, Italy, October 22-29, 2017, pp. 3352–3361, 2017.
  * [40] X. Yan, A. Rastogi, R. Villegas, K. Sunkavalli, E. Shechtman, S. Hadap, E. Yumer, and H. Lee, “MT-VAE: learning motion transformations to generate multimodal human dynamics,” in Computer Vision - ECCV 2018 - 15th European Conference, Munich, Germany, September 8-14, 2018, Proceedings, Part V, pp. 276–293, 2018.
  * [41] A. Bhattacharyya, B. Schiele, and M. Fritz, “Accurate and diverse sampling of sequences based on a “best of many" sample objective,” in 2018 IEEE Conference on Computer Vision and Pattern Recognition, CVPR 2018, Salt Lake City, UT, USA, June 18-22, 2018, pp. 8485–8493, 2018.
  * [42] T. Salzmann, M. Pavone, and M. Ryll, “Motron: Multimodal probabilistic human motion forecasting,” in IEEE/CVF Conference on Computer Vision and Pattern Recognition, CVPR 2022, New Orleans, LA, USA, June 18-24, 2022, pp. 6447–6456, 2022.
  * [43] Y. Yuan and K. M. Kitani, “Diverse trajectory forecasting with determinantal point processes,” in 8th International Conference on Learning Representations, ICLR 2020, Addis Ababa, Ethiopia, April 26-30, 2020, 2020\.
  * [44] H. Ma, J. Li, R. Hosseini, M. Tomizuka, and C. Choi, “Multi-objective diverse human motion prediction with knowledge distillation,” in IEEE/CVF Conference on Computer Vision and Pattern Recognition, CVPR 2022, New Orleans, LA, USA, June 18-24, 2022, pp. 8151–8161, 2022.
  * [45] L. Dang, Y. Nie, C. Long, Q. Zhang, and G. Li, “Diverse human motion prediction via gumbel-softmax sampling from an auxiliary space,” in MM ’22: The 30th ACM International Conference on Multimedia, Lisboa, Portugal, October 10 - 14, 2022, pp. 5162–5171, 2022.
  * [46] C. Ionescu, D. Papava, V. Olaru, and C. Sminchisescu, “Human3.6m: Large scale datasets and predictive methods for 3d human sensing in natural environments,” IEEE Trans. Pattern Anal. Mach. Intell., vol. 36, no. 7, pp. 1325–1339, 2014.
  * [47] N. Mahmood, N. Ghorbani, N. F. Troje, G. Pons-Moll, and M. J. Black, “AMASS: archive of motion capture as surface shapes,” in 2019 IEEE/CVF International Conference on Computer Vision, ICCV 2019, Seoul, Korea (South), October 27 - November 2, 2019, pp. 5441–5450, 2019.
  * [48] L. Sigal, A. O. Balan, and M. J. Black, “Humaneva: Synchronized video and motion capture dataset and baseline algorithm for evaluation of articulated human motion,” Int. J. Comput. Vis., vol. 87, no. 1-2, pp. 4–27, 2010\.
  * [49] M. Loper, N. Mahmood, J. Romero, G. Pons-Moll, and M. J. Black, “SMPL: a skinned multi-person linear model,” ACM Trans. Graph., vol. 34, no. 6, pp. 248:1–248:16, 2015.
  * [50] X. Bie, W. Guo, S. Leglaive, L. Girin, F. Moreno-Noguer, and X. Alameda-Pineda, “Hit-dvae: Human motion generation via hierarchical transformer dynamical VAE,” CoRR, vol. abs/2204.01565, 2022.
  * [51] C. Guo, X. Zuo, S. Wang, S. Zou, Q. Sun, A. Deng, M. Gong, and L. Cheng, “Action2motion: Conditioned generation of 3d human motions,” in MM ’20: The 28th ACM International Conference on Multimedia, Virtual Event / Seattle, WA, USA, October 12-16, 2020, pp. 2021–2029, 2020.
  * [52] M. Petrovich, M. J. Black, and G. Varol, “Action-conditioned 3d human motion synthesis with transformer VAE,” in 2021 IEEE/CVF International Conference on Computer Vision, ICCV 2021, Montreal, QC, Canada, October 10-17, 2021, pp. 10965–10975, 2021.
  * [53] D. P. Kingma and J. Ba, “Adam: A method for stochastic optimization,” in ICLR, 2015.


# Programming patchy particles to form three-dimensional dodecagonal
quasicrystals

Daniel F. Tracey  Physical and Theoretical Chemistry Laboratory, Department of
Chemistry, University of Oxford, Oxford, OX1 3QZ, United Kingdom  Eva G. Noya
Instituto de Química Física Rocasolano, Consejo Superior de Investigaciones
Científicas, CSIC, Calle Serrano 119, 28006 Madrid, Spain  Jonathan P. K. Doye
<EMAIL_ADDRESS>Physical and Theoretical Chemistry Laboratory,
Department of Chemistry, University of Oxford, Oxford, OX1 3QZ, United Kingdom

###### Abstract

Model patchy particles have been shown to be able to form a wide variety of
structures, including symmetric clusters, complex crystals and even two-
dimensional quasicrystals. Here, we investigate whether we can design patchy
particles that form three-dimensional quasicrystals, in particular targeting a
quasicrystal with dodecagonal symmetry that is made up of stacks of two-
dimensional quasicrystalline layers. We obtain two designs that are able to
form such a dodecagonal quasicrystal in annealing simulations. The first is a
one-component system of 7-patch particles but with wide patches that allow
them to adopt both 7- and 8-coordinated environments. The second is a ternary
system that contains a mixture of 7- and 8-patch particles, and is likely to
be more realizable in experiments, for example, using DNA origami. One
interesting feature of the first system is that the resulting quasicrystals
very often contain a screw dislocation.

##  I Introduction

Quasicrystals are characterized by long-range order (exemplified by sharp
Bragg peaks in their diffraction patterns) in the absence of translational
periodicity, and often exhibit symmetries not feasible in periodic crystals.
The initial examples were found in metallic alloy systems with the most common
symmetries being icosahedral[1] and decagonal,[2] but with metastable
octagonal[3] and dodecagonal[4] quasicrystals also being observed. More
recently, an increasing number of soft-matter quasicrystals have been
discovered,[5] nearly all of which are dodecagonal; these include polymeric
micelles,[6, 7, 8] and nanoparticle superlattices.[9, 10]

Why do quasicrystals, rather than periodic crystals, sometimes form, and what
are the key features of the interparticle interactions that cause this? Theory
and simulations can play an important role in answering these questions.
Quasicrystals have been most commonly observed in simulations for isotropic
potentials with multiple features in the potential whose positions are tuned
to favour quasicrystalline order.[11, 12, 13, 14, 15, 16, 17, 18] For
ultrasoft potentials (i.e. particles without a hard repulsive core), it can be
shown more rigorously that quasicrystals with a particular symmetry can be
favoured when the ratio of two inverse length scales in the Fourier transform
of the potential has a specific value.[19, 20, 21]

Simulation examples that are exceptions to this “two-length-scale”
stabilization mechanism include quasicrystals formed from hard particles
(tetrahedra[22, 23] and triangular bipyramids[24] can form dodecagonal
quasicrystals) and from “patchy” particles with directional bonding;[25, 26,
27] the latter are the focus of this paper. In particular, two-dimensional
patchy particles with five and seven regularly arranged patches have been
found to form dodecagonal quasicrystals that are square-triangle tilings where
the bonds are equally likely to be oriented along twelve equivalent
directions.[25]

The 5-patch system has become a model system to better understand quasicrystal
behaviour. For example, phase diagrams have been computed, showing that the
dodecagonal quasicrystal is the thermodynamically stable phase for an
intermediate range of temperature and patch width, and its stabilization is
due to its greater entropy compared to crystalline approximants.[28] The
growth dynamics of this quasicrystal have also been studied.[29]

Figure 1:  (a) A dodecagonal motif that is a key motif in two-dimensional
dodecagonal quasicrystalline square-triangle tilings. (b) Three possible local
motifs in square-triangle tilings. The labels mirror those used for Frank-
Kasper phases.[30, 31]

A key structural motif in these two-dimensional dodecagonal quasicrystals is
the dodecagon shown in Fig. 1(a). The 5- and 6-coordinate environments in the
dodecagon are both feasible when the particles’ patches are sufficiently wide.
To realize such a quasicrystal, for particles with instead a fixed maximum
number of interaction partners, requires a mixture of 5- and 6-patch
particles,[26] this being achieved experimentally in systems of multi-arm DNA
tiles on a surface.[32]

Here, we extend this approach to obtain 3D dodecagonal quasicrystals,
exploring the additional complexities that the added dimension brings. It is
noteworthy that many of the dodecagonal quasicrystals previously observed in
simulations have been for 2D[12, 25, 13, 14, 26] or quasi-2D[33, 17] systems.

##  II Methods

###  II.1 Potential

The interaction between particles is described using the patchy-particle model
introduced in Refs. 34, 35. In this pair potential, the interaction is
described by a Lennard-Jones repulsive core and an attractive tail modulated
by angular and torsional dependent functions:

$V_{ij}(\mathbf{r}_{ij},\mathbf{\Omega}_{i},\mathbf{\Omega}_{j})=\begin{cases}V_{\mathrm{LJ}}^{\prime}(r_{ij})&:r_{ij}<\sigma_{\mathrm{LJ}}^{\prime}\\\
V_{\mathrm{LJ}}^{\prime}(r_{ij})\underset{{\mathrm{patch~{}pairs~{}}\alpha,\beta}}{\max}\left[\varepsilon_{\alpha\beta}V_{\mathrm{ang}}(\mathbf{\hat{r}}_{ij},\mathbf{\Omega}_{i},\mathbf{\Omega}_{j})V_{\mathrm{tor}}(\mathbf{\hat{r}}_{ij},\mathbf{\Omega}_{i},\mathbf{\Omega}_{j})\right]&:r_{ij}\geqslant\sigma_{\mathrm{LJ}}^{\prime}\end{cases},$
(1)

where $\mathbf{r}_{ij}$ is the interparticle vector, $\alpha$ and $\beta$ are
patches on particles $i$ and $j$, respectively, $\mathbf{\Omega}_{i}$ is the
orientation of particle $i$, $V_{\mathrm{LJ}}^{\prime}(r)$ is a cut-and-
shifted Lennard-Jones (LJ) potential and $\sigma_{\mathrm{LJ}}^{\prime}$
corresponds to the distance at which $V_{\mathrm{LJ}}^{\prime}$ passes through
zero. We set the cutoff distance $r_{\mathrm{cut}}=2.5\,\sigma_{\mathrm{LJ}}$.
$\varepsilon_{\alpha\beta}$ is a measure of the relative strength of the
interactions between patches $\alpha$ and $\beta$; $\varepsilon_{\alpha\beta}$
is set to zero for patches that do not interact.

The angular modulation term $V_{\mathrm{ang}}$ is a measure of how directly
the patches $\alpha$ and $\beta$ point at each other, and is given by

$V_{\mathrm{ang}}(\mathbf{\hat{r}}_{ij},\mathbf{\Omega}_{i},\mathbf{\Omega}_{j})=\exp\left(-\frac{\theta_{\alpha
ij}^{2}}{2\sigma_{\mathrm{ang}}^{2}}\right)\exp\left(-\frac{\theta_{\beta
ji}^{2}}{2\sigma_{\mathrm{ang}}^{2}}\right),$ (2)

$\theta_{\alpha ij}$ is the angle between the patch vector
$\mathbf{\hat{P}}_{i}^{\alpha}$, representing the patch $\alpha$, and
$\mathbf{\hat{r}}_{ij}$. $\sigma_{\mathrm{ang}}$ is a measure of the angular
width of the patch.

The torsional modulation term $V_{\mathrm{tor}}$ describes the variation in
the potential as either of the particles is rotated about the interparticle
vector $\mathbf{r}_{ij}$ and is given by

$V_{\mathrm{tor}}(\mathbf{\hat{r}}_{ij},\mathbf{\Omega}_{i},\mathbf{\Omega}_{j})=\exp\left(-\frac{1}{2\sigma_{\mathrm{tor}}^{2}}\left[\min\limits_{\phi^{\mathrm{offset}}_{\alpha\beta}}\left(\phi_{\alpha\beta}-\phi^{\mathrm{offset}}_{\alpha\beta}\right)\right]^{2}\right).$
(3)

where $\phi^{\mathrm{offset}}$ is the preferred value of the torsional angle
$\phi$. To define the torsional angle $\phi_{\alpha\beta}$, a unique reference
vector (here, this is always one of the other patch vectors) is associated
with each patch. In order to capture the symmetry of an environment, more than
one equivalent offset angle can be defined; in these cases we find the minimum
value of $\phi-\phi_{\mathrm{offset}}$ across the set of equivalent offset
angles. To turn off the torsional component of a particular interaction we set
$V_{\mathrm{tor}}=1$. We use $\sigma_{\mathrm{tor}}=2\sigma_{\mathrm{ang}}$
throughout.

We use $\sigma_{\mathrm{LJ}}$ as our unit of length, and the LJ well depth
$\varepsilon_{\mathrm{LJ}}$ as our unit of energy. Values of
$\sigma_{\mathrm{ang}}$ and $\sigma_{\mathrm{tor}}$ are given in radians.
Temperatures are given in reduced form,
$T^{*}=k_{B}T/\varepsilon_{\mathrm{LJ}}$. Full details of each patchy particle
model that we study are tabulated in the Supporting Information.

###  II.2 Simulations

We use Metropolis Monte Carlo to simulate the patchy-particle systems in the
canonical ensemble. The diffusive dynamics this algorithm generates is
appropriate for modelling the dynamics of colloidal particles in solution. To
simulate the larger systems we use a GPU-enabled code that uses the
parallelization strategy introduced in Ref. 36. Like in Ref. 35, assembly is
achieved by slow annealing simulations in which the temperature is linearly
decreased with time. The system is initiated from a disordered, low-density
($\rho=0.1\sigma^{-3}_{\mathrm{LJ}}$) fluid state at a temperature just above
that for which nucleation of a condensed-state cluster can occur on the
simulation time scales.

###  II.3 Structural characterization

We analyse the structure of the clusters that grow in a number of ways. To
identify the different local particle environments we use a common-neighbour
analysis,[37, 26] focussing on those environments that are relevant to the
target dodecagonal quasicrystals. In particular, we identify the three-
dimensional equivalents of the local environments in Fig. 1(b), where the in-
plane bonding is the same as in the figure, but with out-of-plane bonding to
particles directly above and below the particle of interest. Environments not
matching these three are labelled as U (unidentified); this category will
include surface particles with incomplete bonding and particles at defect
sites.

To identify quasicrystallinity, we compute diffraction patterns for our
assembled structures. Dodecagonal quasicrystallinity is indicated by a 12-fold
symmetric pattern when viewed from a direction perpendicular to the
quasicrystalline layers. Specifically, we compute the real part of the
structure factor $S(\mathbf{q})$:

$S(\mathbf{q})=\frac{1}{N}\sum_{j=1}^{N}\sum_{k=1}^{N}e^{-i\mathbf{q}\cdotp\mathbf{r}_{kj}}.$
(4)

We define the normal to the layers as the average direction of the axial
patches on the particles in a cluster, and, after orienting this direction
along $z$, we evaluate the structure factor for $\mathbf{q}=(q_{x},q_{y},0)$.
We evaluate the diffraction patterns averaged over a set of configurations at
the end of our annealing runs.

##  III Results

The basic design approach that we use is to target structures that are stacked
crystalline or quasicrystalline layers, where atoms in adjacent layers lie
directly above each other. The patch designs that mediate interactions within
the layers will be similar to the previous two-dimensional examples.[25, 26]
In addition, axial patches perpendicular to the above mediate interactions
between the layers. The torsional component of the interactions facilitates
the layer formation. Specifically, for the reference vector of the equatorial
patches we use one of the axial patches with preferred offset angles of 0 or
180∘.

Figure 2: Stacked $\sigma$-crystal. (a) The patchy particle design and its
relationship to the $\sigma$ local environment, and the interaction matrix
between the different types of patches. (Note there is a second axial patch of
type 1 (silver) on the opposite side of the particle shown). (b) Two views of
an 864-particle crystalline cluster grown in an annealing simulation. Only the
bond network is visualized. The stacked layers and the in-plane $\sigma$-like
order are both apparent. $\sigma_{\mathrm{ang}}=0.3$.

###  III.1 Stacked-$\sigma$ crystal design

We first illustrate our approach for a crystal made of stacked
“$\sigma$-layers”, where each particle has an in-plane environment identical
to the $\sigma$ environment in Fig. 1(b). The crystal is tetragonal with space
group $P4/mbm$ (127) with the $4g$ Wyckoff sites occupied. As each environment
in the crystal is equivalent, only one particle type is needed. Three types of
equatorial patches are used for the three unique types of in-plane bonds, with
the inter-patch angles exactly matching those of the $\sigma$ environment,
i.e. they are 60∘ or 90∘ if they are part of in-plane triangles or squares,
respectively (Fig. 2(a)). This design fully specifies the target crystal [35]
and, as expected, the crystal forms straightforwardly in annealing simulations
for moderately narrow patches (we use $\sigma_{\mathrm{ang}}=0.3$) where
7-fold coordination is clearly favoured. An example crystallite is illustrated
in Fig. 2(b); note the anisotropic crystal shape reflects the greater
stabilization provided by the larger number of in-plane bonds.

###  III.2 One-component quasicrystal design

We next consider a potential one-component quasicrystal-forming design. The
particle design is very similar to the above, but all five equatorial patches
are now equivalent and regularly-spaced in the equatorial plane, i.e. the in-
plane inter-patch angles are 72∘ (Fig. 3(a)).

Figure 3:  One-component dodecagonal quasicrystal assembly. (a) The patchy
particle design. (b) An example final configuration from an annealing run at
$\sigma_{\mathrm{ang}}=0.52$, viewed from the side and above. The graduated
coloring is to add visual clarity. In the top view each dodecagon is
highlighted by a black circle, and yellow lines connect the centres of edge-
sharing or overlapping dodecagons. (c) Example diffraction patterns for the
above configuration (left) and for a cluster from another annealing runs at
the same $\sigma_{\mathrm{ang}}$ (right) viewed perpendicular to the stacked
layers. Each diffraction pattern is averaged over a set of configurations at
the end of the annealing run. (d) The fraction of particles in the largest
cluster with a given local environment as a function of annealing temperature
in a system of 20 000 particles at $\sigma_{\mathrm{ang}}=0.49$ and 0.52. The
annealing rate for the temperature was 1.2 and $3.87\times 10^{-4}$ per
$10^{5}$ MC cycles, respectively. Further properties for these two annealing
runs are shown in Supplementary Fig. S2.

We first consider assembly at $\sigma_{\mathrm{ang}}=0.3$ where we expect a
strong preference for 7-fold coordination (5-fold coordination in the layers);
$\sigma$ environments are more likely than H environments because the $\sigma$
environment has a slightly better match to the inter-patch angles.[25] The
crystallites resulting from our annealing simulations are very similar to
those in Section III.1, but with some particles in H environments (see
Supplementary Fig. S1).

The more interesting behaviour is anticipated for larger
$\sigma_{\mathrm{ang}}$, because, as the potential becomes less anisotropic,
an increasing propensity for 6-fold in-plane coordination is expected. For the
previously studied 2D systems, this leads to a small region of stability for
the dodecagonal quasicrystal in the ($\sigma_{\mathrm{ang}},T$) phase diagram,
sandwiched between regions where the fluid, $\sigma$ crystal or hexagonal
crystal is most stable.[28] On annealing, two types of assembly pathways were
seen for this 2D system depending on the value of $\sigma_{\mathrm{ang}}$. At
lower $\sigma_{\mathrm{ang}}$, the quasicrystal assembled direct from the
fluid, whereas at slightly larger $\sigma_{\mathrm{ang}}$, the system first
formed a hexagonal crystal before transforming into a quasicrystal at lower
temperature.

We find somewhat similar behaviour in the annealing simulations for our 3D
systems, observing both direct and indirect pathways for quasicrystal
assembly. For example, at $\sigma_{\mathrm{ang}}=0.49$ a quasicrystalline
cluster grows directly from the fluid with the order increasing as the
temperature decreases (Fig. 3(d)). By contrast, at
$\sigma_{\mathrm{ang}}=0.52$ the cluster that first assembles from the gas
phase has a much higher proportion of hexagonal in-plane environments and
unidentified environments and a much lower fraction of $\sigma$ environments.
The cluster in this temperature regime exhibits significant disorder both
within the layers and in terms of a reduced tendency to form planar layers.
However, at low temperature the cluster then transforms into a quasicrystal
with the final cluster having a very similar distribution of environments to
the quasicrystals that grow direct from the gas phase.

Diffraction patterns for the final clusters from two of the annealing runs are
shown in Fig. 3(c); they have clear 12-fold symmetry confirming the
quasicrystalline character. From the top view of this configuration (Fig.
3(b)) the planar layers can be seen to be tilings of squares and triangles,
with the different local environments and the dodecagonal motifs of Fig. 1
evident. Although there is no translational order, orientational order is
maintained through the whole system with the inter-plane bonds equally likely
in twelve equivalent directions. The dodecagonal motifs are often seen to
share edges, forming larger motifs such as square or triangular arrangements
of dodecagons; also present in the structure are examples where the dodecagons
overlap. The network of interconnected dodecagons propagates through the whole
structure.

Figure 4:  Ternary dodecagonal quasicrystal assembly. (a) The three types of
patchy particles, their relationship to the basic dodecagonal motif, and the
interaction matrix. (b) An example final configuration from an annealing run
at $\sigma_{\mathrm{ang}}=0.3$ for a system of $N=20,007$ particles. The
relative proportions of particle types 1, 2 and 3 were 6/13, 6/13 and 1/13,
respectively. The annealing rate for the temperature was $5\times 10^{-5}$ per
$10^{5}$ MC cycles. The colors match the particle types in (a). (c) The
diffraction pattern for this configuration viewed perpendicular to the stacked
layers. (d) The number of particles in the largest cluster broken down by
particle type. Snapshots of the growth of this quasicrystalline cluster are
given in Supplementary Figure S3.

###  III.3 Ternary quasicrystal design

One disadvantage of the above one-component system is that the quasicrystal
only forms for a relatively small range of $\sigma_{\mathrm{ang}}$ in which
multiple coordination environments compete as the potential becomes less
anisotropic. Instead, it is probably experimentally easier to realize patchy
particles with well-defined valence. Therefore, we also explored quasicrystal
formation for multi-component designs where the environments with different
coordination numbers are achieved by particles with different numbers of
patches. These systems are somewhat analogous to the 2D examples studied in
Ref. 26. Such systems have the advantage that quasicrystal formation is likely
to be relatively insensitive to the precise value of the patch width; we
choose to use $\sigma_{\mathrm{ang}}=0.3$; this is in the regime where
directional bonding is dominant. However, multi-component designs also have
potential disadvantages. Having multiple components is likely to have some
adverse affects on the kinetics. Assembly is likely to be somewhat slower as
particles must search for their correct partners, and rearrangements of the
particles when in the condensed state may be significantly slower due to the
greater selectivity in the bonding. Also, there is the potential that the
particles might prefer to phase-separate rather than form an ordered phase
that involves all particle types.

After exploratory testing of different models, the most promising design that
we then studied in more detail is given in Fig. 4(a). The system has three
particle types. Particle type 3 has six in-plane patches and is designed to be
at the centre of the dodecagonal motif. The other two particle types have five
in-plane patches. One (type 2) is designed to bind around the eight-patch
particle to form the hexagonal motif at the centre of the dodecagon. The other
(type 1) is geometrically equivalent to the particle studied in the above one-
component example, with all equatorial patches equivalent, and will typically
be located on the outside of the dodecagonal motifs. Two patch types can bind
to multiple partners. This flexibility allows the particles to organize in
multiple ways, reducing the tendency for crystalline structures to form. The
magenta patches on particle type 2, as well as bonding to the equatorial
patches of particle 1 can also self-interact. The latter allows overlapping
dodecagonal motifs to form. Similarly, the equatorial patches on particle type
1, as well as binding to the magenta patches on particle type 2, can also
self-interact.

Thus, particle type 1 can potentially form a fully-bonded $\sigma$-crystal on
its own, and particle types 2 and 3 could together form a crystal of
overlapping dodecagons without particle type 1. Therefore, this system has the
potential to demix rather than form an ordered ternary phase. In the previous
sections, all attractive patch-patch interactions were of equal strength (i.e.
$\varepsilon_{\alpha\beta}=1$). Here, we used the values of
$\varepsilon_{\alpha\beta}$ to tune the behaviour of the system to prevent
phase separation. In particular, we reduced the strength of the self-
interactions of patches 2 and 6 by 10% (i.e.
$\varepsilon_{22}=\varepsilon_{66}=0.9$ and $\varepsilon_{26}=1$), to slightly
favour bonding between particle types 1 and 2. Without this change, the
initial clusters tended to predominantly involve particles 1, with particles 2
and 3 only being added to the outside of the clusters once the concentration
of particles 1 in the fluid phase had significantly diminished.

In this model we also reduced the strength of the axial interactions on all
particles by 5% (i.e.
$\varepsilon_{11}=\varepsilon_{33}=\varepsilon_{77}=0.95$). This design choice
was simply to change the shape of the assembled clusters, so that they had a
greater lateral extent, thus allowing easier identification of the potential
quasicrystalline order. As in Ref. 26, the composition was chosen to match
that for a triangular crystal of edge-sharing dodecagonal units.

Fig. 4(b) shows one of the clusters resulting from our annealing simulations,
and its diffraction pattern. Again, the diffraction pattern has clear 12-fold
order confirming its quasicrystalline character (Fig. 4(c)). Examination of
the top view of the configuration shows that the planes are square-triangle
tilings, with dodecagons being a common motif although somewhat less prevalent
than for the one-component system. While the system clearly naturally forms a
ternary cluster — the hexagonal motifs are distributed throughout the cluster
— the composition seems to be slightly higher in particle type 1 towards the
centre of the planes. This tendency is also evident in Fig. 4(b); the number
of particles of type 2 in the growing cluster somewhat lags that of type 1,
even though there are equal numbers in the simulation. The initial nucleation
event always occurred around a small hexagonal stack (Supplementary Fig. S3).

Figure 5:  Configurations exhibiting screw dislocations for one-component
quasicrystals formed at (a) $\sigma_{\mathrm{ang}}$=0.43, and (b)
$\sigma_{\mathrm{ang}}$=0.46, and (c) for the stacked $\sigma$-crystal at
$\sigma_{\mathrm{ang}}$=0.3. On the left are the complete clusters with the
exit points of the screw dislocations highlighted by orange circles. On the
right are sections through the clusters that allow the path of the dislocation
through the structures to be seen (highlighted by the orange ellipses). These
configurations were prepared by taking the final configuration from an
annealing simulation, and then quenching it to $T^{*}=0.01$ to remove thermal
noise. The particles are colored according to their potential energy per
patch.

###  III.4 Screw dislocations

An especially interesting feature of the assembled one-component quasicrystals
was the ubiquitous presence of screw dislocations. They were present in the
final configurations for all clusters assembled at
$\sigma_{\mathrm{ang}}=0.43$ and 0.46, and many at 0.49 and 0.52. Their lower
likelihood of occurrence at larger $\sigma_{\mathrm{ang}}$ suggests their
formation is somewhat suppressed when the quasicrystal does not form direct
from the fluid, but rather via an intermediate condensed state.

Two examples are shown in Fig. 5. The screw dislocation runs through the
clusters roughly perpendicular to the quasicrystalline planes and leads to a
ledge on the top surface of the quasicrystal at the exit point of the
dislocation. During growth the ledge spirals around the defect core, leading
to continuous growth without the need for the secondary nucleation of an
island on top of a flat layer. For this reason, screw dislocations are
associated with enhanced crystal growth rates.

In a crystal, the Burgers vector is parallel to the screw dislocation and
perpendicular to the plane of the Burgers circuit. By analogy, the Burgers
vectors in our quasicrystals would be in the axial direction perpendicular to
the quasicrystalline layers. However, defining a Burgers circuit in a
quasicrystal is less straightforward, and the formal theory of dislocations in
quasicrystals makes use of higher-dimensional crystallography.[38]

To aid the visualization of the dislocations, the final configurations from
the annealing simulations were quenched to reduce thermal noise (by simulating
the systems at $T^{*}=0.01$ for $5\times 10^{5}$ MC cycles). The particles
were then colored according to their potential energy. This leads to expected
differences in color between particles in bulk, surface or edge sites, but
also allows the path of the screw dislocation to be easily identified in cuts
through the structures. This color scheme also highlights other defects that
are present in the bulk of the structures, but these are not our focus here.

The first example in Fig. 5(a) has two screw dislocations emerging on the top
surface. Consequently, there is a noticeable tilt to the layers. The cut
clearly shows the helical pattern of bonds around the dislocation core. The
stress associated with the dislocations is quite localized, with clear
deviations from the bulk energy only noticeable close to the dislocation core.
The screw dislocations are generally very straight but with occasional jogs or
kinks. Thus, the screw dislocation in the cut in Fig. 5(a) is only in the
plane of the cut for the top third of the cluster, whereas for the example in
Fig. 5(b) it is visible passing through the entire structure.

Why do screw dislocations occur for the one-component quasicrystal but not the
ternary example? It is probably for two reasons. Firstly, the wide patches
that we use in the one-component case (to facilitate multiple environments)
give rise to smaller energetic penalties for the deformed bonds around the
dislocation core. When we performed simulations for this system with
$\sigma_{\mathrm{ang}}=0.3$, we never observed screw dislocations. By
contrast, when we simulated larger systems of the stacked $\sigma$ crystal at
$\sigma_{\mathrm{ang}}=0.3$, screw dislocations sometimes formed (Fig. 5(c)).
In this case, the more anisotropic shape of the crystals may be relevant as
the total energetic cost of the dislocation will be proportional to the
crystal thickness, thus partially offsetting the greater local cost of the
screw dislocation. Secondly, it may be harder to create a pattern of bonding
at a dislocation core that is compatible with the greater specificity of the
interactions in the ternary case.

Another important question is why and at what stage do the screw dislocations
form. Clearly, for systems of this size, it is very unlikely that screw
dislocation formation is thermodynamically favourable, because the entropy
gain from their presence is unlikely to outweigh their unfavourable energy.
For the same reason, they are unlikely to be present in pre-critical nuclei
which are in quasi-equilibrium with the fluid. This is especially so given the
critical nucleus size is likely to be small at the temperatures at which we
see spontaneous nucleation in our simulations. Therefore, it is likely that
they become incorporated into growing post-critical clusters. The question
then is what happens when, during the growth, a local structural perturbation
forms on the surface of the cluster that could potentially lead to a fully-
formed screw dislocation: does it get annealed out because of its less
favourable energy, or does it get incorporated and locked in due to the boost
it provides to the growth kinetics? Under some conditions the latter seems to
happen sufficiently frequently that the resulting quasicrystals invariably
include screw dislocations.

Figure 6:  Configurations from the early stages of cluster growth that show
the presence of screw dislocations. (a) $N\approx 540$ at
$\sigma_{\mathrm{ang}}$=0.46. (b) $N\approx 1990$ at
$\sigma_{\mathrm{ang}}$=0.43. As in Fig. 5, the cluster configuration was
quenched to $T^{*}=0.01$ to remove thermal noise, and the particles colored
according to their potential energy per patch.

Analysis of the growing clusters indicates that the screw dislocations
generally form in the early stages of cluster growth. Fig. 6(a) shows a
cluster with just over 500 particles that has a screw dislocation near to the
edge of the cluster. The larger cluster in Fig. 6(b) (around 2 000 particles)
already has two screw dislocations passing through the complete cluster. It is
likely that as the clusters grow larger it becomes less likely that new screw
dislocations are incorporated, because of the increased formation costs due to
the greater thickness of the clusters.

Even once fully formed, the screw dislocations can potentially be removed if
they diffuse to the edge of the clusters and are then expelled. However,
dislocation motion was slow on our computational time scales, so we rarely saw
this occur.

##  IV Conclusions

We have shown that patchy particles can be programmed to form 3D dodecagonal
quasicrystals. We have introduced two designs that achieve this. The first is
a one-component system with wide patches that enable the particles to adapt to
multiple environments, in particular ones with different coordination numbers.
The second is a ternary mixture, where particles with seven and eight patches
explicitly match the different coordination environments in the quasicrystal;
the well-defined valence of these particles is likely to make them easier to
realize experimentally. In this system, the patch-patch interaction strengths
had to be tuned in order to prevent phase separation of the ternary mixture.
The quasicrystallinity of both systems was confirmed by the 12-fold symmetry
of their diffraction patterns.

One interesting feature of the one-component systems was that the majority of
the quasicrystals that we grew had screw dislocations passing through the
structures in the axial direction. The features of this system that
facilitated their formation were the wide patches and the layer-like
structure. The screw dislocations were incorporated into the quasicrystals in
the early stages of growth. In these instances, when the initial structural
perturbation at the cluster surface that could potentially lead to a screw
dislocation occurred, the tendency for the screw dislocation to facilitate
growth must have won out over the tendency for the defect to be annealed out
because of its less favourable potential energy. It is perhaps noteworthy
that, although there have been many simulations of screw dislocations,[39] we
are unaware of any simulation studies where screw dislocations have
spontaneously appeared during crystal growth.

The quasicrystals that have been observed for soft matter systems have almost
invariably been of dodecagonal symmetry.[5] Thus, it would be particularly
interesting if patchy particle designs could be developed that form three-
dimensional quasicrystals with other symmetries. We have now very recently
achieved this for icosahedral quasicrystals, where the target quasicrystalline
structures were generated by projection from six-dimensional cubic
lattices.[27] Generating appropriate target structures for other axial
quasicrystals, however, may be challenging. In this study, the task was made
easier firstly because we chose the simplest way one could conceive of making
a 3D axial quasicrystal, namely by a simple stacking of two-dimensional
quasicrystalline layers. Secondly, the square-triangle character of the
dodecagonal quasicrystalline layers is straightforwardly compatible with a
particle model because the edges of the tiling are all of the same length.
Thirdly, we also already knew that two-dimensional patchy particles could form
such quasicrystalline layers.[25, 26]

Experimental examples for which the structure of the quasicrystal has been
solved may provide one source of target structures.[40] However, the high
coordination numbers typically present in metallic alloys may present
additional challenges for a patchy-particle approach. Another potential source
of target structures is from simulation studies where 3D axial quasicrystals
have formed in particle-based models with complex isotropic potentials. These
include dodecagonal quasicrystals that have a Frank-Kasper-like particle
decoration of square and triangular prismatic units,[11, 18] two decagonal
examples[16, 17] and one octagonal example.[17]

How might such patchy particles be realized experimentally? Although there has
been much progress in synthesizing colloidal analogues of patchy
particles,[41, 42, 43] the degree of control over the geometry of the patches
is not yet sufficient to realize particles like those envisaged here. An
alternative might be to use the techniques of DNA nanotechnology.[44] For
example, DNA origami polyhedra[45] whose vertices have been decorated with
short DNA single strands that facilitate inter-vertex binding have been used
to form a variety of simple crystal structures.[46] For this system, the
inter-vertex binding is quite flexible, but we have also proposed an approach
in which a six-helix bundle extends from each vertex to provide more
directional and torsionally specific interactions between the DNA origami.[27]
The sequences of the single strands that allow inter-origami binding can be
tuned to provide the requisite specificity, and it would be relatively
straightforward to choose the type of polyhedron and tune their edge lengths
to match the geometry of the patchy particle. As the origami particles would
have a well-defined maximum number of binding partners, they would be more
appropriate for realizing the ternary quasicrystal. For example, in this
scheme, a hexagonal bipyramidal DNA origami would be used to create an
analogue of particle type 3 with its six in-plane and two out-of-plane
patches, and pentagonal bipyramidal origamis would be used for particle types
1 and 2.

###### Acknowledgements.

D.F.T. is grateful for funding via the ESPRC Centre for Doctoral Training in
Theory and Modelling in Chemical Sciences, under grant EP/L015722/1. E.G.N.
acknowledges funding from the Agencia Estatal de Investigación (AEI) and the
Fondo Europeo de Desarrollo Regional (FEDER) under grant numbers
FIS2015-72946-EXP(AEI) and FIS2017-89361-C3-2-P(AEI/FEDER,UE), and the
European Union’s Horizon 2020 research and innovation programme under the
Marie Skłodowska-Curie grant agreement No. 734276. We acknowledge the use of
the University of Oxford Advanced Research Computing (ARC) facility and the
Cambridge Service for Data Driven Discovery (CSD3).

## Data Availability

The data that support the findings of this study are available from the
corresponding author upon reasonable request.

## References

  * Shechtman et al. [1984] D. Shechtman, I. Blech, D. Gratias, and J. W. Cahn, Phys. Rev. Lett. 53, 1951 (1984).
  * Bendersky [1985] L. Bendersky, Phys. Rev. Lett. 55, 1461 (1985).
  * Wang et al. [1987] N. Wang, H. Chen, and K. H. Kuo, Phys. Rev. Lett. 59, 1010 (1987).
  * Ishimasa et al. [1985] T. Ishimasa, H.-U. Nissen, and Y. Fukano, Phys. Rev. Lett. 55, 511 (1985).
  * Dotera [2011] T. Dotera, Isr. J. Chem. 51, 1197 (2011).
  * Zeng et al. [2004] X. Zeng, G. Ungar, Y. Liu, V. Percec, A. E. Dulcey, and J. K. Hobbs, Nature 428, 157 (2004).
  * Fischer et al. [2011] S. Fischer, A. Exner, K. Zielske, J. Perlich, S. Deloudi, W. Steurer, P. Lindner, and S. Förster, Proc. Natl. Acad. Sci. USA 108, 1810 (2011).
  * Zhang and Bates [2012] J. Zhang and F. S. Bates, J. Am. Chem. Soc. 134, 7636 (2012).
  * Talapin et al. [2009] D. V. Talapin, E. V. Shevchenko, M. I. Bodnarchuk, X. Ye, J. Chen, and C. B. Murray, Nature 461, 964 (2009).
  * Ye et al. [2017] X. Ye, J. Chen, M. E. Irrgang, M. Engel, A. Dong, S. C. Glotzer, and C. B. Murray, Nat. Mater. 16, 214 (2017).
  * Dzugutov [1993] M. Dzugutov, Phys. Rev. Lett. 70, 2924 (1993).
  * Engel and Trebin [2007] M. Engel and H.-R. Trebin, Phys. Rev. Lett. 98, 225505 (2007).
  * Dotera et al. [2014] T. Dotera, T. Oshiro, and P. Ziheri, Nature 506, 208 (2014).
  * Pattabhiraman et al. [2015] H. Pattabhiraman, A. P. Gantapara, and M. Dijkstra, J. Chem. Phys. 143, 164905 (2015).
  * Engel et al. [2015] M. Engel, P. F. Damasceno, C. L. Phillips, and S. C. Glotzer, Nat. Mater. 14, 109 (2015).
  * Ryltsev et al. [2015] R. Ryltsev, B. Klumov, and N. Chtchelkatchev, Soft Matter 11, 6991 (2015).
  * Damasceno et al. [2017] P. F. Damasceno, S. C. Glotzer, and M. Engel, J. Phys: Condens. Matter 29, 234005 (2017).
  * Ryltsev and Chtchelkatchev [2017] R. Ryltsev and N. Chtchelkatchev, Soft Matter 13, 5076 (2017).
  * Subramanian et al. [2016] P. Subramanian, A. H. Archer, E. Knobloch, and A. M. Rucklidge, Phys. Rev. Lett. 117, 075501 (2016).
  * Savitz et al. [2018] S. Savitz, M. Babadi, and R. Lifshitz, IUCrJ 5, 247 (2018).
  * Ratliff et al. [2019] D. J. Ratliff, A. H. Archer, P. Subramanian, and A. M. Rucklidge, Phys. Rev. Lett. 123, 148004 (2019).
  * Haji-Akbari et al. [2009] A. Haji-Akbari, M. Engel, A. S. Keys, X. Zheng, R. G. Petschek, P. Palffy-Muhoray, and S. C. Glotzer, Nature 462, 773 (2009).
  * Je et al. [2021] K. Je, S. Lee, E. G. Teich, M. Engel, and S. C. Glotzer, Proc. Natl. Acad. Sci. USA 118, e2011799118 (2021).
  * Haji-Akbari et al. [2011] A. Haji-Akbari, M. Engel, and S. C. Glotzer, Phys. Rev. Lett. 107, 215702 (2011).
  * van der Linden et al. [2012] M. N. van der Linden, J. P. K. Doye, and A. A. Louis, J. Chem. Phys. 136, 054904 (2012).
  * Reinhardt et al. [2017] A. Reinhardt, J. S. Schreck, and J. P. K. Doye, J. Phys.: Condens. Matter 29, 014006 (2017).
  * Noya et al. [2020] E. G. Noya, C. K. Wong, P. Llombart, and J. P. K. Doye, submitted (2020).
  * Reinhardt et al. [2013] A. Reinhardt, F. Romano, and J. P. K. Doye, Phys. Rev. Lett. 110, 255503 (2013).
  * Gemeinhardt et al. [2018] A. Gemeinhardt, M. Martinsons, and M. Schmiedeberg, Eur. Phys. J. E 41, 126 (2018).
  * Frank and Kasper [1958] F. C. Frank and J. S. Kasper, Acta Crystallogr. 11, 184 (1958).
  * Frank and Kasper [1959] F. C. Frank and J. S. Kasper, Acta Crystallogr. 12, 483 (1959).
  * Liu et al. [2019] L. Liu, Z. Li, Y. Li, and C. Mao, J. Am. Chem. Soc. 141, 4248 (2019).
  * Metere et al. [2016] A. Metere, P. Oleynikov, M. Dzugutov, and S. Lidin, Soft Matter 12, 8869 (2016).
  * Wilber et al. [2009] A. W. Wilber, J. P. K. Doye, A. A. Louis, and A. C. F. Lewis, J. Chem. Phys. 131, 175102 (2009).
  * Tracey et al. [2019] D. F. Tracey, E. G. Noya, and J. P. K. Doye, J. Chem. Phys. 151, 224506 (2019).
  * Anderson et al. [2013] J. A. Anderson, E. Jankowski, T. L. Grubb, M. Engel, and S. C. Glotzer, J. Comput. Phys. 254, 27 (2013).
  * Clarke and Jonsson [1993] A. S. Clarke and H. Jonsson, Phys. Rev. E 47, 3975 (1993).
  * Sandbrink and Schmeideberg [2014] M. Sandbrink and M. Schmeideberg, Phys. Rev. B 90, 064108 (2014).
  * Maresca et al. [2018] F. Maresca, D. Dragoni, G. Csányi, N. Marzari, and W. A. Curtin, NPJ Comput. Mater. 4, 69 (2018).
  * Steurer and Deloudi [2009] W. Steurer and S. Deloudi, _Crystallography of quasicrystals: Concepts, methods and structures_ (Springer, 2009).
  * Glotzer and Solomon [2007] S. C. Glotzer and M. Solomon, Nat. Mater. 6, 557 (2007).
  * Li et al. [2020] W. Li, H. Palis, R. Mérindol, J. Majimel, S. Ravaine, and E. Duguet, Chem. Soc. Rev. 49, 1955 (2020).
  * Liu et al. [2020] M. Liu, X. Zheg, V. Grebe, D. J. Pine, and M. Weck, Nat. Mater. 19, 1354 (2020).
  * Ma et al. [2020] N. Ma, B. Minevich, J. Liu, M. Ji, Y. Tian, and O. Gang, Top. Curr. Chem. 379, 36 (2020).
  * Veneziano et al. [2016] R. Veneziano, S. Ratanalert, F. Zhang, H. Yan, W. Chiu, and M. Bathe, Science 352, 1534 (2016).
  * Tian et al. [2020] Y. Tian, J. R. Lhermitte, L. Bai, T. Vo, H. L. Xin, H. Li, R. Li, M. Fukuto, K. G. Yager, J. S. Kahn, et al., Nat. Mater. 19, 789 (2020).