# A noise-robust Multivariate Multiscale Permutation Entropy for two-phase
flow characterisation

###### Abstract

Using a graph-based approach, we propose a multiscale permutation entropy to
explore the complexity of multivariate time series over multiple time scales.
This multivariate multiscale permutation entropy ($\operatorname{MMPE_{G}}$)
incorporates the interaction between channels by constructing an underlying
graph for each coarse-grained time series and then applying the recent
permutation entropy for graph signals. Given the challenge posed by noise in
real-world data analysis, we investigate the robustness to noise of
$\operatorname{MMPE_{G}}$ using synthetic time series and demonstrating better
performance than similar multivariate entropy metrics.

We also apply $\operatorname{MMPE_{G}}$ to study two-phase flow data, an
important industrial process characterised by complex, dynamic behaviour. To
this end, we process multivariate Electrical Resistance Tomography (ERT) data
and extract multivariate multiscale permutation entropy values.
$\operatorname{MMPE_{G}}$ characterises the flow behaviour transition of two-
phase flow by incorporating information from different scales. The
experimental results show that $\operatorname{MMPE_{G}}$ is sensitive to the
dynamic of flow patterns, allowing us to distinguish between different flow
patterns. We show that our method is noise-robust, which is suitable for
analysis of the complexity of multivariate time series and characterising two-
phase flow recordings.

Index Terms—  permutation entropy, graph signal, entropy metrics, complexity,
noise, two-phase flow.

##  1 Introduction

Two-phase flow is important in many industries, including chemical processes,
petroleum exploitation, nuclear engineering, and transportation [1].
Theoretical and experimental studies have been sought to characterise two-
phase flow, including mathematical approaches [2, 3], fluid dynamics analysis
[4], using high-speed camera [5], and others. These methods have not fully
resolved the complexity and dynamic behaviours of the flow patterns,
especially in the interaction between different channels in the multivariate
signals.

The analysis of a complex system like two-phase flow can benefit from
nonlinear analysis metrics. Univariate permutation entropy [6] has been used
as a nonlinear measure of complexity; it is a computationally fast algorithm,
its performance under noise conditions has been investigated [7], has been
used for discriminating two-phase flow dynamics [8] and characterise
autoregressive processes [9]. Most of the physical systems’ signals are
multivariate. Therefore, several univariate entropy metrics have been
generalised to a multivariate setting, including multivariate sample entropy
[10], multivariate dispersion entropy [11], and multivariate permutation
entropy [12], among others. Some of them have been recently used to analyse
phase flow and are applied to characterise the behaviour of the two-phase flow
system [13, 14, 15]. A prior implementation of multivariate permutation
entropy exists, but it has a major limitation (i.e., it does not consider
cross-channel information).

The effect of noise on the multivariate permutation entropy can lead to
inaccurate values. Some improvements have been proposed to deal with this
problem, such as a multivariate weighted permutation entropy [16]. However,
these multivariate methods analyse each time series separately; hence, the
cross-channel information will be lost and increase the number of parameters
used for the entropy computation. To this end, we present a multiscale
multivariate permutation entropy based on constructing a Cartesian graph
product.

_Contributions:_ This paper introduces a multiscale algorithm to analyse
multivariate time series based on the permutation entropy for graph signals.
We apply the algorithm to a set of synthetic data and two-phase flow data. We
show that it improves the performance of univariate permutation entropy and
classical multivariate permutation entropy because $\operatorname{MMPE_{G}}$
consider the interaction between the different data channels, it is robust to
noise, and it is useful to detect the complexity at different scales of the
phase flow patterns.

_Structure of the manuscript:_ The outline of the manuscript is as follows:
Sec. 2 introduces the permutation entropy for graph signals and Cartesian
graph product. Sec. 3 presents the multiscale multivariate permutation entropy
$\operatorname{MMPE_{G}}$. In Sec. 4, we use synthetic signals to show
$\operatorname{MMPE_{G}}$ is robust to noise and it is applied to analyse
phase flow data (Sec. 5). The conclusions and future work are presented in
Sec. 6.

##  2 Graphs and permutation entropy

This section presents the Cartesian graph product and the recently permutation
entropy for analyse graph signals.

Cartesian graph product. The _Cartesian product_ of two graphs
${G}=(\mathcal{V},\mathcal{E})$ and
${G}^{\prime}=(\mathcal{V}^{\prime},\mathcal{E}^{\prime})$, denoted
${G}\square{G}^{\prime}$, is the graph defined by the vertex set:
$\mathcal{V}({G}\square{G}^{\prime})=\mathcal{V}\times\mathcal{V}^{\prime}=\\{\,(v,v^{\prime})\,|\,v\in\mathcal{V}\text{
and }v^{\prime}\in\mathcal{V}^{\prime}\,\\}\>$. Two vertices $(v,v^{\prime})$
and $(u,u^{\prime})$ are adjacent in ${G}\square{G}^{\prime}$ if and only if
either $v=u$ and $v^{\prime}$ is adjacent to $u^{\prime}$ in ${G}^{\prime}$,
or $v^{\prime}=u^{\prime}$ and $v$ is adjacent to $u$ in ${G}$. The graph
product is a useful structure to model multidomain signals [17] and it is the
perturbation of a periodic graph [18].

Permutation entropy for graph signals ($\operatorname{PE_{G}}$). Let
$G=(\mathcal{V},\mathcal{E})$ be a graph, $\mathbf{A}$ its adjacency matrix
and $\textbf{X}=\left\\{x_{i}\right\\}_{i=1}^{n}$ be a signal on the graph,
$\operatorname{PE_{G}}$ is defined in [19, 20] as follows:

  1. 1.

For $2\leq m\in\mathbb{N}$ the _embedding dimension_ , $L\in\mathbb{N}$ the
_delay time_ and for all $i=1,2,\dots,n$, we define
$y_{i}^{kL}=\frac{1}{\lvert\mathcal{N}_{kL}(i)\rvert}\sum_{j\in\mathcal{N}_{kL}(i)}x_{j}=\frac{1}{\lvert\mathcal{N}_{kL}(i)\rvert}(\mathbf{A}^{kL}\mathbf{X})_{i}\;$,
where $\mathcal{N}_{k}(i)=\\{\,j\in\mathcal{V}\,|\,\scalebox{0.9}[1.0]{it
exists a walk on $k$ edges joining $i$ and $j$}\,\\}\;$. Then, we construct
the embedding vector $\textbf{y}_{i}^{m,L}\in\mathbb{R}^{m}$ given by

$\textbf{y}_{i}^{m,L}=\left(y_{i}^{kL}\right)_{k=0}^{m-1}=\left(y_{i}^{0},y_{i}^{L},\dots
y_{i}^{(m-1)L}\right)\;.$

  2. 2.

The embedding vector $\textbf{y}_{i}^{m,L}$ is arranged in increasing order
and is assigned to one of $k=m!$ permutation (or patterns)
$\pi_{1},\pi_{2},\dots,\pi_{k}$.

  3. 3.

For the distinct permutation, the relative frequency is denoted by
$p(\pi_{1}),p(\pi_{2}),\dots,p(\pi_{k})$. The permutation entropy
$\operatorname{PE_{G}}$ for the graph signal X is computed as the normalised
Shannon entropy

$\operatorname{PE_{G}}=-\dfrac{1}{\ln(m!)}\sum_{i=1}^{m!}p(\pi_{i})\ln
p(\pi_{i})\;.$

##  3 Multiscale Multivariate Permutation Entropy ($\operatorname{MMPE_{G}}$)

Time series from phase flow data contain multiple temporal scale structures.
Consider only a single scale have limited capability and they can only assess
the system’s irregularity at a single temporal scale. Multiple scales need to
be analysed to understand the dynamics in the signals and to describe the
properties of its model.

Here, we propose a nonlinear multivariate multiscale methodology based on
graph signals to analyse such data. Let
$\mathbf{X}=\\{x_{t,s}\\}_{t=1,2,\dots,n}^{s=1,2,\dots,p}$ be a multivariate
time series with $p-$channels of length $n$ and with $I_{p}$ the graph of
interactions between channels. The $\operatorname{MMPE_{G}}$ algorithm relies
on a three-step procedure:

  1. 1.

Coarse-grained procedure. From the original univariate signal $\mathbf{X}$, we
derive multiple successive coarse-grained versions by averaging the time data
points within non-overlapping time segments of increasing length, $\epsilon$ ,
referred to as the scale factor. Each element of the coarse-grained time
series, $\mathbf{U}^{\epsilon}=\\{u_{i,j}^{\epsilon}\\}_{i=1,2,\dots,\lfloor
n/\epsilon\rfloor}^{j=1,2,\dots,p}$ , is calculated as:

$u_{i,j}^{\epsilon}=\dfrac{1}{\epsilon}\sum_{t=(i-1)\epsilon+1}^{i\epsilon}x_{t,j}\;,\quad\text{for
}1\leq j\leq p\text{ and }1\leq i\leq n/\epsilon\;.$

The length of each coarse-grained time series is $\epsilon$ times shorter than
the original one. For $\epsilon=1$, we get the original series, i.e,
$\mathbf{U}^{1}=\mathbf{X}$.

We use this coarse-grained procedure for simplicity, however other approaches
to construct the coarse-graining in multiscale entropy exist [21, 22].

  2. 2.

Graph construction associated to a multivariate signal. For each coarse-
grained multivariate $\mathbf{U}^{\epsilon}$, we construct the graph
$\mathcal{G}_{\mathbf{U}^{\epsilon}}$ given by:

$\mathcal{G}_{\mathbf{U}^{\epsilon}}:=\overrightarrow{P}_{\lfloor
n/\epsilon\rfloor}\square I_{p}\;.$

Usually, we will consider two basic graph structures, $I_{p}=K_{p}$ the
complete graph with $p$ vertices and $I^{\prime}_{p}=\emptyset_{p}$ the empty
graph on $p$ vertices i.e., without channels-interaction (see Fig. 1 for an
example with $p=3$).

Fig. 1: Graph structures considering as interaction graphs.

The directed graph $\overrightarrow{P}_{n}$ has $n$ vertices with adjacency
matrix ${\mathbf{A}}_{\overrightarrow{P}_{n}}$ (size $n\times n$), and the
graph $I_{p}$ has $p$ vertices with adjacency matrix ${\mathbf{A}}_{I_{p}}$
(size $p\times p$), then the adjacency matrix
${\mathbf{A}}_{{\overrightarrow{P}_{n}}\square I_{p}}$ (size $np\times np$) of
the Cartesian product of both graphs ${\overrightarrow{P}_{n}}\square I_{p}$
is given by

${\mathbf{A}}_{{\overrightarrow{P}_{n}}\square
I_{p}}=\mathbf{A}_{{\overrightarrow{P}_{n}}}\otimes\mathbb{I}_{p}+\mathbb{I}_{n}\otimes\mathbf{A}_{I_{p}}\;;$

where $\otimes$ denotes the Kronecker product of matrices and $\mathbb{I}_{n}$
denotes the $n\times n$ identity matrix.

  3. 3.

$\mathbf{PE}$ for graph signals. We consider $\mathbf{U}^{\epsilon}$ as a
signal defined on the graph $\mathcal{G}_{\mathbf{U}^{\epsilon}}$, i.e.,
$\mathbf{U}^{\epsilon}\colon\mathcal{V}(\mathcal{G}_{\mathbf{U}^{\epsilon}})\longrightarrow\mathbb{R}\;.$
The _multivariate permutation entropy ( $\operatorname{MMPE_{G}}$)_ is defined
as the permutation entropy $\operatorname{PE_{G}}$ for the graph signal
$\mathbf{U}^{\epsilon}$ defined on the graph
$\mathcal{G}_{\mathbf{U}^{\epsilon}}$, i.e.,

$\displaystyle\operatorname{MMPE_{G}}=\operatorname{PE_{G}}(\mathbf{U}^{\epsilon})\;.$

To be noted, if we consider $I_{p}=\emptyset$, our method leads to the
multiscale multivariate permutation entropy previously presented in [12].
However, $\operatorname{MMPE_{G}}$ uses information on between channel
interaction (represented in the graph $I_{p}$) and leads to a more robustness
method in the presence of noise, as we will see in the next section.

##  4 Robustness to noises

To demonstrate the effectiveness of the method, we will analyse the effect of
additive Gaussian noise on the performance of $\operatorname{MMPE_{G}}$ on the
Lorenz system, this system has important applications in mechanics, biology,
and circuit theory [23]. The Lorenz system is is given by the system of
ordinary differential equations:

$\displaystyle{x^{\prime}}$ $\displaystyle=\sigma(y-x),$
$\displaystyle{y^{\prime}}$ $\displaystyle=x(\rho-z)-y,$
$\displaystyle{z^{\prime}}$ $\displaystyle=xy-\beta z.$

A simulation for the values $\rho=50,\sigma=19$ and $\beta=6$, with initial
state $x=-1,y=0$ and $z=1$ is depicted in Fig. 2.

Fig. 2: Time series or sequences of Lorenz system.

We calculate the multivariate permutation entropy, considering the two
previous graph structures, $I_{3}$ and $I^{\prime}_{3}$, and we plot the
entropy values as a function of the scale factor $\epsilon$ (see Fig. 3). In
both graph structures, the entropy values increase with the scale factor.
However, our $\operatorname{MMPE_{G}}$ shows more stability to the change of
scale factor, and the complexity of the multivariate system is more constant,
while $\operatorname{MMPE}$ [12, 16] shows more variability in its values,
changing from no complexity at lower scales to almost random behaviour at
higher scales for the same Lorenz system. This robustness would be important
for noise.

Fig. 3: Entropy values as a function of the temporal scale $\epsilon$.

To demonstrate the robustness of the method, we will add white Gaussian noise
to the multivariate signals defined by the Lorenz system. Fig. 4 shows the
effects of WGN on multivariate multiscale $\operatorname{MMPE}$ and
$\operatorname{MMPE_{G}}$ for the Lorenz system. The signal-to-noise ratio of
WGN is set for all the channels as 10dB, 20dB, 30dB, 40dB and 50dB,
respectively. The noise has low influence with both methods for scales higher
than ten. The influence of the noise is very important for
$\operatorname{MMPE}$ in lower scales, while the impact is smaller for our
algorithm, showing the robustness of $\operatorname{MMPE_{G}}$. Similar
results are obtained when the noise is added to one or two channels.

Fig. 4: Entropy values when WGN is added.

The analysis of dynamics with $\operatorname{MMPE_{G}}$ includes temporal
information and channel interaction. Experiments where a noise channel is
added to the Lorenz system confirmed that $\operatorname{MMPE_{G}}$ has a good
performance. The entropy profile did not change noticeably with the addition
of the noise channel (results not shown due to space constraints). Preserving
all channels (without deleting noise) is important to compare entropy values
for different systems but with the same number of channels.

##  5 Flow analysis with $\operatorname{MMPE_{G}}$

The two-phase flow experiment was carried out at Tianjin University [13],
using Electrical Resistance Tomography (ERT). Based on the principle that the
conductivity of medium differs, ERT collects boundary voltages between
electrodes placed around the pipe by applying electric currents to obtain
conductivity distribution of two-phase flow inside the pipe. A constant
electrical current of $50$ kHz is adopted as the exciting signal, and the data
acquisition rate is 120 frames/s. A $16-$electrode ERT obtains $16\times
13=208$ voltage data. Given the redundancy in the measurements and similar to
[24], we extract feature vectors $V_{Ri}$ from each electrode to reduce the
dimension as follows:
$V_{Ri}=\frac{1}{13}\sum_{j=1}^{13}(V_{ij}-V_{i{j_{0}}})/V_{i{j_{0}}}$ where
$V_{ij}$ is the measure voltage value, $V_{i{j_{0}}}$ is the $V_{ij}$ when the
pipe is full of water and $1\leq R\leq 16$. Finally, to reduce the
computational time, the 16 features vectors $V_{Ri}$ are compressed into 4
time series, by average four $V_{Ri}$ electrode belonging to the same set
group [13]. Then, we applied multivariate multiscale permutation entropy to
characterise the two-phase flow.

Fig. 5: a) Signals of $4$ set of electrodes under annular flow and b) is the
graph modelling cross-channel interaction.

Two-phase flow series require preprocessing to eliminate the noise, but our
method is noise-robust as shown in Sec. 4 and has better results than
classical $\operatorname{MMPE}$. Hence we do not need filtering to obtain a
good characterisation of the dynamic in phase flow. The noise causes only
small variations of entropy values in the lower scales; hence, we will work
with the original data without additional preprocessing.

For two-phase flow, gas and water were mixed. Water velocity ranged from $0.4$
m/s to $2.9$ m/s and gas velocities from $0.06$ m/s to $5.64$ m/s. The flow
pattern observed with this experiments and analysed are patterns known as
Bubble flow, Slug flow, Churn flow, and Annular flow. For the analysis, $89$
experiments carried out under different conditions of flow rate of gas and
water and the typical length of each recording is $\sim$1400 time samples. We
perform the $\operatorname{MMPE_{G}}$ on the signals of the four flow
patterns, and plot the corresponding $\operatorname{MMPE_{G}}$ versus scale
factor in Fig. 6. The coarse-grain process reduces the time series length;
consequently, the results show more variability in high-scale factors.

Fig. 6: Mean and standard deviation from $\operatorname{MMPE_{G}}$ values of
signals for different flow patterns.

Entropy of flow patterns. The bubble flow shows the highest
$\operatorname{MMPE_{G}}$ value in almost all scales compared with the other
three patterns. The presence of small bubbles [25] on the regime produces
complex time series on all four sets of nodes. Because of gravity, the bubbles
affect more the top electrodes than the bottom. Then, the high complexity of
the time series and different response of the sensors leads to the highest
$\operatorname{MMPE_{G}}$ values.

The signals acquired in slug flow show a periodic fluctuation induced by the
repeated occurrence of big gas bubbles. Moreover, the length of the gas bubble
could be visually identified in the signal because the four electrode sets are
all affected by the big gas bubbles flowing over the measured cross-section
and show high voltages resulting from the bubble. Hence, the periodicity on
temporal dimension and similar effects in all the electrodes produce the
lowest $\operatorname{MMPE_{G}}$ values for all scales.

The symmetrical distribution of the liquid film around the pipe perimeter in
the annular flow leads to the electrodes in the annular flow showing similar
fluctuations [26], but with different amplitudes. The liquid film at the
bottom of the pipe is thicker than the top due to the gas velocity; see Fig.
5. The complexity depends more on the temporal dimension than the structural
dimension. Hence, the values of $\operatorname{MMPE_{G}}$ overlap in low
scales for slug and annular flow, making them indistinguishably but resulting
in less complexity than bubble flow. For higher scale,
$\operatorname{MMPE_{G}}$ is able to distinguish between the slug, annular,
and bubble flow.

Churn flow is the flow with more dynamic changes along the scale values
because it is a highly unstable flow [15]. In lower scales, churn flow is
similar to the slug and annular flow. The presence of bubbles in the churn
flow and the interaction and coalesce with each other produces the highest
entropy values than slug and annular flow for scales between 3-14 and more
complexity than bubble for scales between 7-11 because of the presence of no
complete slugs and small waves. Periodic waves are relevant for higher scales,
decreasing the values of $\operatorname{MMPE_{G}}$ and making them similar to
the annular flow.

##  6 Conclusions and future work

This paper proposes a multiscale nonlinear methodology to analyse multivariate
time series using the concept of permutation entropy. Contrary to the previous
state of the art, our method allows for the consideration of cross-channel
interactions, thanks to the exploitation of graph products and our recent
formulation of permutation entropy for graph signals.
$\operatorname{MMPE_{G}}$ is robust with respect to additive noise, making it
suitable for analysis of the complexity of multivariate time series and
characterisation of two-phase flow recordings.

Some future lines of research are: additional statistical analysis of
$\operatorname{MMPE_{G}}$ values (for example, analysis of the slopes in the
scale factor), to explore other underlying graph construction (tensor product,
strong product), and to consider other types of data as well, including
biomedical signals.

## References

  * [1] S. Levy,  Two-phase flow in complex systems,  John Wiley & Sons, 1999.

  * [2] D.A. Drew,  “Mathematical modeling of two-phase flow,”  Annual review of fluid mechanics, vol. 15, no. 1, pp. 261–291, 1983\.
  * [3] H. B. Stewart and B. Wendroff,  “Two-phase flow: models and methods,”  Journal of Computational Physics, vol. 56, no. 3, pp. 363–409, 1984\.
  * [4] M. Ishii and T. Hibiki,  Thermo-fluid dynamics of two-phase flow,  Springer Science & Business Media, 2010.
  * [5] H. Li, Y. Zhou, B. Sun, and Y. Yang,  “Multi-scale chaotic analysis of the characteristics of gas-liquid two-phase flow patterns,”  Chinese Journal of Chemical Engineering, vol. 18, no. 5, pp. 880–888, 2010.
  * [6] C. Bandt and B. Pompe,  “Permutation entropy: a natural complexity measure for time series,”  Physical review letters, vol. 88, no. 17, pp. 174102, 2002.
  * [7] I. Veisi, N. Pariz, and A. Karimpour,  “Fast and robust detection of epilepsy in noisy EEG signals using permutation entropy,”  in 2007 IEEE 7th international symposium on bioinformatics and bioengineering. IEEE, 2007, pp. 200–203.
  * [8] C.L. Fan, N.D. Jin, X.T. Chen, and Z.K. Gao,  “Multi-scale permutation entropy: A complexity measure for discriminating two-phase flow dynamics,”  Chinese Physics Letters, vol. 30, no. 9, pp. 090501, 2013.
  * [9] A. Dàvalos, M. Jabloun, P. Ravier, and O. Buttelli,  “Multiscale permutation entropy: Statistical characterization on autoregressive and moving average processes,”  in 2019 27th European Signal Processing Conference (EUSIPCO), 2019, pp. 1–5.
  * [10] M.U. Ahmed and D. P. Mandic,  “Multivariate multiscale entropy: A tool for complexity analysis of multichannel data,”  Physical Review E, vol. 84, no. 6, pp. 061918, 2011.
  * [11] H. Azami, A. Fernández, and J. Escudero,  “Multivariate multiscale dispersion entropy of biomedical times series,”  Entropy, vol. 21, no. 9, pp. 913, 2019.
  * [12] F.C. Morabito, D. Labate, F. Foresta, A. Bramanti, G. Morabito, and I. Palamara,  “Multivariate multi-scale permutation entropy for complexity analysis of Alzheimer’s disease EEG,”  Entropy, vol. 14, no. 7, pp. 1186–1202, 2012.
  * [13] C. Tan, J. Zhao, and F. Dong,  “Gas–water two-phase flow characterization with electrical resistance tomography and multivariate multiscale entropy analysis,”  ISA transactions, vol. 55, pp. 241–249, 2015.
  * [14] Z.K. Gao, M.S. Ding, H. Geng, and N.D. Jin,  “Multivariate multiscale entropy analysis of horizontal oil–water two-phase flow,”  Physica A: Statistical Mechanics and Its Applications, vol. 417, pp. 7–17, 2015.
  * [15] Z.K. Gao, Y.X. Yang, L.S. Zhai, M.S. Ding, and N.D. Jin,  “Characterizing slug to churn flow transition by using multivariate pseudo wigner distribution and multivariate multiscale entropy,”  Chemical Engineering Journal, vol. 291, pp. 74–81, 2016.
  * [16] Y. Yin and P. Shang,  “Multivariate weighted multiscale permutation entropy for complex time series,”  Nonlinear Dynamics, vol. 88, no. 3, pp. 1707–1722, 2017.
  * [17] S.K. Kadambari and S.P. Chepuri,  “Learning product graphs from multidomain signals,”  in ICASSP 2020-2020 IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP). IEEE, 2020, pp. 5665–5669.
  * [18] J.F. Fabila-Carrasco,  The discrete magnetic laplacian: geometric and spectral preorders with applications,  Ph.D. thesis, Universidad Carlos III de Madrid, 2020.
  * [19] J.S. Fabila-Carrasco, C. Tan, and J. Escudero,  “Permutation entropy for graph signals,”  IEEE Transactions on Signal and Information Processing over Networks, vol. 8, pp. 288–300, 2022.
  * [20] J.S. Fabila-Carrasco, C. Tan, and J. Escudero,  “Multivariate permutation entropy, a Cartesian graph product approach,”  in 30th European Signal Processing Conference, 2022.
  * [21] H. Azami and J. Escudero,  “Coarse-graining approaches in univariate multiscale sample and dispersion entropy,”  Entropy, vol. 20, no. 2, pp. 138, 2018.
  * [22] J.F. Valencia, A. Porta, M. Vallverdu, F. Claria, R. Baranowski, E. Orlowska-Baranowska, and P. Caminal,  “Refined multiscale entropy: Application to 24-h holter recordings of heart period variability in healthy and aortic stenosis subjects,”  IEEE Transactions on Biomedical Engineering, vol. 56, no. 9, pp. 2202–2213, 2009.
  * [23] M.W. Hirsch, S. Smale, and R.L. Devaney,  Differential equations, dynamical systems, and an introduction to chaos,  Academic press, 2012.
  * [24] C. Tan, Y. Shen, K. Smith, F. Dong, and J. Escudero,  “Gas–liquid flow pattern analysis based on graph connectivity and graph-variate dynamic connectivity of ERT,”  IEEE Transactions on Instrumentation and Measurement, vol. 68, no. 5, pp. 1590–1601, 2018.
  * [25] X. Dong, C. Tan, and F. Dong,  “Gas–liquid two-phase flow velocity measurement with continuous wave ultrasonic doppler and conductance sensor,”  IEEE Transactions on Instrumentation and Measurement, vol. 66, no. 11, pp. 3064–3076, 2017.
  * [26] C. Tan, W. Dai, H. Yeung, and F. Dong,  “A Kalman estimation based oil–water two-phase flow measurement with CRCC,”  International Journal of Multiphase Flow, vol. 72, pp. 306–317, 2015.


# Testing the quasicentroid molecular dynamics method on gas-phase ammonia

Christopher Haggard  Yusuf Hamied Department of Chemistry, University of
Cambridge, Lensfield Road, Cambridge, CB2 1EW, UK.  Vijay Ganesh Sadhasivam
Yusuf Hamied Department of Chemistry, University of Cambridge, Lensfield Road,
Cambridge, CB2 1EW, UK.  George Trenins  Laboratory of Physical Chemistry, ETH
Zürich, 8093 Zürich, Switzerland.
Stuart C. Althorpe  Yusuf Hamied Department of Chemistry, University of
Cambridge, Lensfield Road, Cambridge, CB2 1EW, UK.

(November 2021)

###### Abstract

Quasicentroid molecular dynamics (QCMD) is a path-integral method for
approximating nuclear quantum effects in dynamics simulations, which has given
promising results for gas- and condensed-phase water. Here, by simulating the
infrared spectrum of gas-phase ammonia, we test the feasibility of extending
QCMD beyond water. Overall, QCMD works as well for ammonia as for water,
reducing or eliminating blue shifts from the classical spectrum without
introducing the artificial red-shifts or broadening associated with other
imaginary-time path-integral methods. However, QCMD gives only a modest
improvement over the classical spectrum for the position of the symmetric bend
mode, which is highly anharmonic (since it correlates with the inversion
pathway). We expect QCMD to have similar problems with large-amplitude degrees
of freedom in other molecules, but otherwise to work as well as for water.

##  I Introduction

A variety of imaginary-time path-integralChandler and Wolynes (1981);
Parrinello and Rahman (1984) methods, including (thermostatted) ring-polymer
molecular dynamics ([T]RPMD) Craig and Manolopoulos (2004); Habershon  _et
al._ (2013); Rossi, Ceriotti, and Manolopoulos (2014); Rossi, Kapil, and
Ceriotti (2018); Welsch  _et al._ (2016) and centroid molecular dynamics (CMD)
Cao and Voth (1994a, b); Hone, Rossky, and Voth (2006); Medders and Paesani
(2015); Suleimanov, Aoiz, and Guo (2016), have been developed for
approximating nuclear quantum effects in simulations of dynamical properties.
These methods resemble classical MD in the extended phase-space of the
imaginary-time Feynman paths or ‘ring-polymers’. They were introduced
heuristically, but can also be thought of as approximations to ‘Matsubara
dynamics’ Hele  _et al._ (2015a); Trenins and Althorpe (2018); Jung, Videla,
and Batista (2019); Althorpe (2021); Chowdhury and Huo (2021), which is the
classical, quantum-Boltzmann-conserving, dynamics that emerges when jagged
imaginary-time paths are smoothed to remove real-time quantum coherence;
[T]RPMD is a short-time approximation to Matsubara dynamics; CMD is a mean-
field approximation, equivalent to averaging out the quantum fluctuations
about the ring-polymer centroids Hele  _et al._ (2015b); Hele (2017); Althorpe
(2021).

The recently developed quasicentroid molecular dynamics (QCMD) Trenins,
Willatt, and Althorpe (2019) method has been applied so far only to water, and
is similar to CMD, except that the fluctuations are mean-field averaged around
a curvilinear ‘quasicentroid’ rather than a cartesian centroid. The resulting
QCMD infrared spectra (of gas-phase water and of the q-TIP4P/F Habershon,
Markland, and Manolopoulos (2009) model of liquid water and ice) show none of
the artificial red-shifting and broadening that affects the CMD stretch band
in water (below 400 K in the gas-phase and 300 K in the liquid), with the
positions of the bands lining up almost perfectly with the exact quantum
spectrum (with a small temperature-independent blue shift), and the
intensities of the fundamental bands approximating well the intensities of the
quantum bands (in gas-phase water) Trenins, Willatt, and Althorpe (2019);
Benson, Trenins, and Althorpe (2019). These results are very promising,
although QCMD is currently roughly 100 times more expensive than TRPMD, which
gives band positions for water which are as good as QCMD but artificially
broadens the lineshapes.

The reason QCMD works so well for water is that each point on the
quasicentroid potential of mean force corresponds to an ensemble of ring
polymers with the same average bond length. At low temperatures, this
condition prevents polymers that stretch around rotational (or librational)
curves from spuriously lowering the free energy near inner turning points by
forming artificial instantons, as happens in CMD Trenins and Althorpe (2018);
Witt  _et al._ (2009); Ivanov  _et al._ (2010). As a result, QCMD gives more
compact ring-polymer distributions than CMD, and hence a better approximation
to Matsubara dynamics.

Here, we investigate whether these advantages can be generalised to molecules
other than water by extending QCMD to treat gas-phase ammonia. A (non-
dissociating) water molecule is especially simple to treat because it has no
internal degrees of freedom capable of highly anharmonic large-amplitude
motion. We know from recent studies of overtones, combination bands and Fermi
resonances, that coupling between the centroid and the fluctuation modes can
have a major effect when anharmonicity is important Plé  _et al._ (2021);
Benson and Althorpe (2021). Such coupling is neglected by QCMD (and also by
CMD, TRPMD and classical MD). We therefore pay particular attention in what
follows to how well QCMD performs for the symmetric stretch band in ammonia,
which correlates with the inversion pathway.

##  II Methodology

It is straightforward to extend the QCMD treatment of water in ref. Trenins,
Willatt, and Althorpe, 2019 to ammonia. We define curvilinear centroids

$\displaystyle
R_{i}=\frac{1}{N}\sum_{k=1}^{N}r_{i}^{(k)}\;\;\;\;\;\;\Theta_{i}=\frac{1}{N}\sum_{k=1}^{N}\theta_{i}^{(k)}$
(1)

where $i=1,2,3$ denotes the bond angle or bond length, and $k=1,\dots,N$
denotes the imaginary-time replica or ‘bead’. For each replica, the bond
lengths and angles are functions of the cartesian bead coordinates of the four
atoms, which we will refer to collectively as ${\bf
q}\equiv\\{q^{(k)}_{l}\\}$, with $k=1,\dots,N$ and $l=1,\dots,12$. These
coordinates roughly correspond to the stretch and bend vibrational normal mode
coordinates. In addition to these 6 internal coordinates (DOFs), 6 external
coordinates are needed to uniquely define the quasicentroid configuration.
These are obtained from Eckart constraints Eckart (1935); Jørgensen (1978);
Kudin and Dymarsky (2005) similar to those in ref. Trenins, Willatt, and
Althorpe, 2019, which orient the ring-polymer with respect to the
quasicentroid unit. We will denote the full set of quasicentroid coordinates
($\\{R_{i},\Theta_{i}\\}$ plus the Eckart constraints) by $\boldsymbol{\xi}$.

The QCMD potential of mean force ${\cal F}(\boldsymbol{\xi})$ is defined to be
the free energy obtained by averaging the ring-polymer fluctuations around
$\boldsymbol{\xi}$, namely

$e^{-\beta{\cal
F}(\boldsymbol{\xi})}=Z_{qc}(\boldsymbol{\xi})=\int\\!d\textbf{p}^{\prime}\int\\!d\textbf{q}^{\prime}\,e^{-\beta_{N}W_{N}(\textbf{q}^{\prime},\textbf{p}^{\prime})}\;\delta(\boldsymbol{\xi}^{\prime}-\boldsymbol{\xi})$
(2)

where $W_{N}(\textbf{q},\textbf{p})$ is the ring-polymer Hamiltonian, and
${\bf p}^{\prime}$ are the cartesian momenta conjugate to the cartesian bead
coordinates ${\bf q}^{\prime}$. The mean forces and torques on the
quasicentroid are thus

$\displaystyle f_{R_{i}}(\boldsymbol{\xi})$
$\displaystyle=-\frac{\partial{\cal F}(\boldsymbol{\xi})}{\partial
R_{i}}\approx-\frac{1}{N}\left\langle\sum_{k=1}^{N}\frac{\partial U_{N}({\bf
q})}{\partial r_{i}^{(k)}}\right\rangle_{\boldsymbol{\xi}}$ (3) $\displaystyle
f_{\Theta_{i}}(\boldsymbol{\xi})$ $\displaystyle=-\frac{\partial{\cal
F}(\boldsymbol{\xi})}{\partial\Theta_{i}}\approx-\frac{1}{N}\left\langle\sum_{k=1}^{N}\frac{\partial
U_{N}({\bf q})}{\partial\theta_{i}^{(k)}}\right\rangle_{\boldsymbol{\xi}}$ (4)

where $U_{N}({\bf q})$ is the ring-polymer potential energy minus the spring
term, and

$\left\langle\dots\right\rangle_{{\boldsymbol{\xi}}}=\frac{1}{Z_{qc}(\boldsymbol{\xi})}\int\\!d\textbf{p}^{\prime}\int\\!d\textbf{q}^{\prime}\,e^{-\beta_{N}W_{N}(\textbf{q}^{\prime},\textbf{p}^{\prime})}\dots\;\delta(\boldsymbol{\xi}^{\prime}-\boldsymbol{\xi})$
(5)

The approximations in Eqs. (3) and (4) follow from assuming that the ring-
polymer distribution is sufficiently compact that the polymer-spring
contribution to the force is negligible. As discussed in ref. Trenins,
Willatt, and Althorpe, 2019, this assumption of compactness also allows one to
propagate the quasicentroid dynamics in cartesian coordinates (so that one
does not need to work with curvilinear momenta and mass matrices), and (most
crucially) to use the dynamics of the quasicentroid as a proxy for the
dynamics of the centroid. Each of these approximations contributes an error to
the static quantum Boltzmann distribution (sampled by the quasicentroid
dynamics), and it is thus necessary to check that the QCMD static properties
agree well with those of standard PIMD.

##  III Numerical details

To propagate the QCMD dynamics of (gas-phase) ammonia, we converted the
curvilinear forces of Eqs. (3) and (4) to cartesians (using the relations
given in the Supplementary Material), and propagated the dynamics using the
adiabatic QCMD (AQCMD) algorithm of Trenins et al. Trenins, Willatt, and
Althorpe (2019) at 300 K (using $N=32$ replicas) and 150 K ($N=64$), using the
ammonia potential energy and dipole-moment surfaces of Yurchenko et al.
Yurchenko  _et al._ (2009) At each temperature, a total of 32 initial
geometries were equilibrated for 15 ps using TRPMD with a global path-integral
Langevin thermostat (PILE-G) Ceriotti  _et al._ (2010) ($\lambda=0.5,\
\tau_{0}=100$ fs), after which they followed QCMD dynamics for 2 ps (to re-
equilibrate to the QCMD distribution), then QCMD for a further 5 ps
(production runs). The AQCMD algorithm used a timestep of $0.1/\gamma$ fs with
$\gamma=32$ (300 K) and 64 (150 K), with a local path-integral Langevin
thermostat (PILE-L) attached to the ring-polymers and a global Langevin
thermostat Bussi and Parrinello (2008) to the quasicentroids, as in ref.
Trenins, Willatt, and Althorpe, 2019.

During QCMD simulations, we observed (rare) numerical instabilities occurring
in the vicinity of inversion events.111The inversion process is a thermally
activated rare event along the symmetric bend coordinate, promoted primarily
by ring polymer fluctuations. In QCMD, these fluctuations do not sufficiently
carry over to the quasicentroid due to strong anharmonicity, so, beyond the
harmonic limit the quasicentroid is no longer close enough to the centroid and
the QCMD approximation breaks down. A small fraction of the QCMD trajectories
($\sim 1.6\%$) were hence discarded when the ammonia attempted to invert,
since the Eckart constraint cannot handle large amplitude motion (LAM) Szalay
(1988); Sayvetz (1939) 222Here we are referring to the LAM of the centroid
from the quasicentroid, not the inversion pathway (which is itself a LAM)..

Figure 1:  Infrared absorption spectra for gas-phase ammonia at 300 K,
computed by extending QCMD as described in the text (black line), compared
with the damped quantum spectrum adapted from ref. Yurchenko  _et al._ , 2009
(grey shading). Also shown are the results of standard CMD, TRPMD and
classical MD simulations. The dashed lines were obtained by multiplying the
$\nu_{1}/\nu_{3}+\nu_{2}$ and $\nu_{1}/\nu_{3}+\nu_{4}$ combination bands
(except for the TRPMD bands, which have very wide tails) by the post-
processing correction factor of ref. Plé  _et al._ , 2021; Benson and
Althorpe, 2021. The absorption intensities in the two panels are scaled in the
ratio 1:13. Figure 2: Same as Fig. 1 for 150 K. The absorption intensities in
the two panels are scaled in the ratio 1:19.

The spectra were obtained from the dipole derivative time-autocorrelation
function $C_{\mathbf{\dot{\mu}}\cdot\mathbf{\dot{\mu}}}(t)$. Since volume is
not defined for the gas-phase simulations, the power spectrum

$\tilde{I}_{\mathbf{\dot{\mu}}}\left(\omega\right)=\beta\int_{-\infty}^{\infty}\\!dt\
e^{-i\omega t}C_{\mathbf{\dot{\mu}}\cdot\mathbf{\dot{\mu}}}(t)f(t)$ (6)

was calculated in lieu of the IR spectrum Habershon, Fanourgakis, and
Manolopoulos (2008), where $f(t)$ is a Hann window of width $\tau=750$ fs.

##  IV Results

The resulting QCMD spectra are shown in Figs. 1 (300 K) and 2 (150 K), where
they are compared with the results of CMD and TRPMD (computed using standard
PIMD methods), classical MD, and quantum dynamics (obtained by Boltzmann-
weighting the line list found in ref. Yurchenko  _et al._ , 2009 and
convolving with $f(t)$). The combination-band regions are also shown
multiplied by the perturbative correction factor of ref. Yao and Overend,
1976, which was found in refs. Benson and Althorpe, 2021 and Plé  _et al._ ,
2021 to account for most of the increase in intensity of the non-fundamental
bands that results from the neglected coupling between the Matsubara dynamics
of the centroid and the fluctuation modes.

With the exception of the symmetric bend ($\nu_{2}$), QCMD performs almost as
well for ammonia as for gas-phase water. The QCMD asymmetric bend ($\nu_{4}$)
and stretch bands ($\nu_{1}/\nu_{3}$) line up with the TRPMD bands, but are
not artificially broadened. Both QCMD and TRPMD give a similar (temperature-
independent) blue shift for the stretch band for ammonia (35 cm-1) to gas-
phase water (60 cm-1), because they are both affected by the neglect of
centroid-fluctuation coupling and real-time coherence. The QCMD combination
bands are also close in position to the quantum bands and the intensities are
in good agreement when multiplied by the perturbative correction factor
(dashed lines in Figs. 1 and 2). This suggests that the ring-polymer
distributions around the quasi-centroids are compact, such that QCMD gives a
good approximation to Matsubara dynamics, and also samples a good
approximation to the exact quantum Boltzmann statistics. This is borne out by
the QCMD average bond lengths and bond angles, which are within better than
$0.2\%$ of the values computed using standard PIMD.

Figure 3: Expanded plot of the symmetric-bend bands of Figs. 1 and 2. Figure
4: Schematic plot of the potential energy (green) and the potentials of mean
force for QCMD (black) and CMD (blue) along the inversion coordinate $\rho$.
The QCMD curvilinear centroid constraint artificially prevents the formation
of the tunnelling instanton at the barrier.

For the symmetric bend, QCMD gives a $\sim$50 cm-1 temperature-dependent blue
shift with respect to the quantum band, making this the only QCMD band that
does not line up with TRPMD at 150 K (see Fig. 3). Discarding the inverting
QCMD trajectories (see above) 333The Eckart conditions fail to invert the
quasicentroid unit when the ring-polymer inverts, which necessitated the
discarding of inverting trajectories. Since the timescale of this inversion
process is higher than any of the vibration timescales, this procedure ought
to have a minimal effect on the spectrum.  may have biased the sampling of
this band and thus skewed its shape, but such an error is likely to be minor,
since discarding the inverting trajectories in TRPMD does not affect its
spectrum. There may also be an error in this band from the inability of QCMD
and TRPMD to describe coherent tunnelling and thus to reproduce the 35 cm-1
bend-tunnelling splitting Yurchenko  _et al._ (2005), although such an error
would be unlikely to shift the overall position of the band, and presumably
has a similar effect on both QCMD and TRPMD. Nevertheless, the (planar)
inversion barrier does give a clue as to what causes the QCMD blue shift. In
CMD (and TRPMD), the ring-polymer lowers its potential energy at the barrier
by forming a (non-artificial) instanton Richardson and Althorpe (2009) along
the inversion coordinate $\rho$ Yurchenko  _et al._ (2009),

$\displaystyle\rho=\arccos\frac{2}{\sqrt{3}}\sin{\left(\frac{\theta_{1}+\theta_{2}+\theta_{3}}{6}\right)}$
(7)

as illustrated in Fig. 4 444Note that $\rho$ defined here and the inversion
coordinate $\overline{\rho}$ in ref.30 are related as
$\rho=\frac{\pi}{2}-\overline{\rho}$. However, in QCMD, the curvilinear
centroid constraint at the barrier ($\Theta_{1}+\Theta_{2}+\Theta_{3}=2\pi$)
forces each of the $N$ ammonia replicas to be planar, which prevents the
instanton from forming.555Strictly, the quasicentroid constraint does not
prevent the instanton geometry from forming, but assigns it to a quasicentroid
ensemble that corresponds to a bent geometry, where its contribution to the
free energy is negligible.  Although the barrier itself is not important, we
can expect that a similar constriction in the range of geometries sampled also
increases the QCMD potential of mean force on the barrier-side of the
potential wells.

In other words, the QCMD curvilinear centroid constraints, which annihilate
the artificial instantons, also annihilate the genuine instanton at the
barrier top, thus giving an artificial blue shift in the symmetric bend. The
CMD cartesian centroid constraints permit both the artificial and the genuine
instantons to form, which seems to have resulted in an almost perfect
cancellation of errors in the frequency of the symmetric bend at 150 K (but
not at other temperatures). Of the methods tested, the best estimate of the
symmetric bend frequency is given by TRPMD, which does not interfere with the
instanton, and thus exerts the correct force on the centroid (although not on
the fluctuation modes). Note that a variety of other semiclassical methods,
including path-integral Liouville dynamicsLiu and Zhang (2016) and semi-
classical initial value representation (SC-IVR) methodsMiller (2001); Conte,
Aspuru-Guzik, and Ceotto (2013), also give a good estimate of this frequency.

It is tempting to try to fix this problem by designing a centroid coordinate
that annihilates artificial instantons but preserves real ones. However,
besides being messy and semi-empirical, such a proposed solution neglects the
real problem here, which is that the anharmonic Matsubara dynamics of the
symmetric bend involves so much centroid-fluctuation coupling that is probably
necessary to remove at least some of the fluctuation modes from the mean-
field, and to explicitly couple their dynamics to that of the quasicentroid.
Such an approach would re-introduce the Matsubara phase, and is thus unlikely
to be practical.

However, we should keep in mind that the ‘bad’ symmetric-bend frequency is no
worse than that predicted by classical MD, and that overall the agreement
between QCMD and the quantum spectrum (Figs. 1 and 2) is a big improvement on
that provided by TRPMD and CMD.

##  V Conclusions

These tests on gas-phase ammonia suggest that the advantages of QCMD can be
generalised to many systems other than water, with the important caveat that
degrees of freedom that are both (Boltzmann-statistically) quantum and highly
anharmonic (such as the ammonia inversion mode) are unlikely to be treated
better by QCMD than by other path-integral methods.

We expect these findings to transfer straightforwardly to the condensed phase,
since intermolecular degrees of freedom (such as librations and centre-of-mass
vibrations) tend to be Boltzmann-statistically classical. QCMD has in fact
already been applied to the q-TIP4P/F model of liquid water and iceTrenins,
Willatt, and Althorpe (2019). Similar calculations would be straightforward
for liquid ammonia (since the only modification needed to Sec. II would be the
addition of the torque estimator $\boldsymbol{\tau}$, as discussed in ref.
Trenins, Willatt, and Althorpe, 2019). More generally, applications of QCMD to
mixtures of water and organic molecules, e.g. clathrates or water-organic
solvent interfaces should be possible. However, QCMD is unlikely to be able to
treat proton-transfer processes (in solution or in the gas-phase), since
proton-transfer coordinates are typically (Boltzmann-statistically) quantum
and highly anharmonicYu and Bowman (2019).

The ‘elephant in the room’ is that QCMD is typically 10-100 times more
expensive than competing path-integral dynamics methods (such as TRPMD and
CMD), on account of the adiabatic algorithm used to generate the quasicentroid
potential of mean force. More efficient algorithms will need to be developed
if QCMD is to become widely applicable.

## Supplementary Material

See the supplementary material for further details of the cartesian to
curvilinear coordinate transformations.

###### Acknowledgements.

It is a pleasure to thank Venkat Kapil for insightful discussions about path-
integral simulation. C.H. acknowledges support from the EPSRC Centre for
Doctoral Training in Computational Methods for Materials Science (Grant No.
EP/L015552/1); V.G.S. from St. John’s College, University of Cambridge,
through a Dr Manmohan Singh Scholarship; G.T. from the Cambridge Philosophical
Society and St Catharine’s College, University of Cambridge.

## Conflicts of Interest

The authors have no conflicts to disclose.

## Author Contributions

C.H. and V.G.S. contributed equally to this work.

## Data Availability

The data that supports the findings of this study are available within the
article.

## Credit Line

This article may be downloaded for personal use only. Any other use requires
prior permission of the author and AIP Publishing. This article appeared in C.
Haggard, V. Sadhasivam, G. Trenins and S. Althorpe, J. Chem. Phys., 155,
174120 (2021). and may be found at https://doi.org/10.1063/5.0068250.

## References

  * Chandler and Wolynes (1981) D. Chandler and P. G. Wolynes, J. Chem. Phys. 74, 4078 (1981).
  * Parrinello and Rahman (1984) M. Parrinello and A. Rahman, J. Chem. Phys. 80, 860 (1984).
  * Craig and Manolopoulos (2004) I. R. Craig and D. E. Manolopoulos, J. Chem. Phys. 121, 3368 (2004).
  * Habershon  _et al._ (2013) S. Habershon, D. E. Manolopoulos, T. E. Markland, and T. F. Miller III, Annu. Rev. Phys. Chem. 64, 387 (2013).
  * Rossi, Ceriotti, and Manolopoulos (2014) M. Rossi, M. Ceriotti, and D. E. Manolopoulos, J. Chem. Phys. 140, 234116 (2014).
  * Rossi, Kapil, and Ceriotti (2018) M. Rossi, V. Kapil, and M. Ceriotti, J. Chem. Phys. 148, 102301 (2018).
  * Welsch  _et al._ (2016) R. Welsch, K. Song, Q. Shi, S. C. Althorpe, and T. F. Miller III, J. Chem. Phys. 145, 204118 (2016).
  * Cao and Voth (1994a) J. Cao and G. A. Voth, J. Chem. Phys. 100, 5093 (1994a).
  * Cao and Voth (1994b) J. Cao and G. A. Voth, J. Chem. Phys. 100, 5106 (1994b).
  * Hone, Rossky, and Voth (2006) T. D. Hone, P. J. Rossky, and G. A. Voth, J. Chem. Phys. 124, 154103 (2006).
  * Medders and Paesani (2015) G. R. Medders and F. Paesani, J. Chem. Theory Comput. 11, 1145 (2015).
  * Suleimanov, Aoiz, and Guo (2016) Y. V. Suleimanov, F. J. Aoiz, and H. Guo, J. Phys. Chem. A 120, 8488 (2016).
  * Hele  _et al._ (2015a) T. J. H. Hele, M. J. Willatt, A. Muolo, and S. C. Althorpe, J. Chem. Phys. 142, 134103 (2015a).
  * Trenins and Althorpe (2018) G. Trenins and S. C. Althorpe, J. Chem. Phys. 149, 14102 (2018).
  * Jung, Videla, and Batista (2019) K. A. Jung, P. E. Videla, and V. S. Batista, J. Chem. Phys. 151, 034108 (2019).
  * Althorpe (2021) S. C. Althorpe, Eur. Phys. J. B 94, 155 (2021).
  * Chowdhury and Huo (2021) S. N. Chowdhury and P. Huo, J. Chem. Phys. 154, 124124 (2021).
  * Hele  _et al._ (2015b) T. J. H. Hele, M. J. Willatt, A. Muolo, and S. C. Althorpe, J. Chem. Phys. 142, 191101 (2015b).
  * Hele (2017) T. J. H. Hele, Mol. Phys. 115, 1435 (2017).
  * Trenins, Willatt, and Althorpe (2019) G. Trenins, M. J. Willatt, and S. C. Althorpe, J. Chem. Phys. 151, 054109 (2019).
  * Habershon, Markland, and Manolopoulos (2009) S. Habershon, T. E. Markland, and D. E. Manolopoulos, J. Chem. Phys. 131, 024501 (2009).
  * Benson, Trenins, and Althorpe (2019) R. L. Benson, G. Trenins, and S. C. Althorpe, Faraday Discuss. 221, 350 (2019).
  * Witt  _et al._ (2009) A. Witt, S. D. Ivanov, M. Shiga, H. Forbert, and D. Marx, J. Chem. Phys. 130, 194510 (2009).
  * Ivanov  _et al._ (2010) S. D. Ivanov, A. Witt, M. Shiga, and D. Marx, J. Chem. Phys. 132, 31101 (2010).
  * Plé  _et al._ (2021) T. Plé, S. Huppert, F. Finocchi, P. Depondt, and S. Bonella, J. Chem. Phys. 155, 104108 (2021).
  * Benson and Althorpe (2021) R. L. Benson and S. C. Althorpe, J. Chem. Phys. 155, 104107 (2021).
  * Eckart (1935) C. Eckart, Phys. Rev. 47, 552 (1935).
  * Jørgensen (1978) F. Jørgensen, Int. J. Quantum Chem. 14, 55 (1978).
  * Kudin and Dymarsky (2005) K. N. Kudin and A. Y. Dymarsky, J. Chem. Phys. 122, 224105 (2005).
  * Yurchenko  _et al._ (2009) S. N. Yurchenko, R. J. Barber, A. Yachmenev, W. Thiel, P. Jensen, and J. Tennyson, J. Phys. Chem. A 113, 11845 (2009).
  * Ceriotti  _et al._ (2010) M. Ceriotti, M. Parrinello, T. E. Markland, and D. E. Manolopoulos, J. Chem. Phys. 133, 124104 (2010).
  * Bussi and Parrinello (2008) G. Bussi and M. Parrinello, Comput. Phys. Commun. 179, 26 (2008).
  * Note (1) The inversion process is a thermally activated rare event along the symmetric bend coordinate, promoted primarily by ring polymer fluctuations. In QCMD, these fluctuations do not sufficiently carry over to the quasicentroid due to strong anharmonicity, so, beyond the harmonic limit the quasicentroid is no longer close enough to the centroid and the QCMD approximation breaks down.
  * Szalay (1988) V. Szalay, J. Mol. Spectrosc. 128, 24 (1988).
  * Sayvetz (1939) A. Sayvetz, J. Chem. Phys. 7, 383 (1939).
  * Note (2) Here we are referring to the LAM of the centroid from the quasicentroid, not the inversion pathway (which is itself a LAM).
  * Habershon, Fanourgakis, and Manolopoulos (2008) S. Habershon, G. S. Fanourgakis, and D. E. Manolopoulos, J. Chem. Phys. 129, 074501 (2008).
  * Yao and Overend (1976) S. J. Yao and J. Overend, Spectrochim. Acta Part A Mol. Spectrosc. 32, 1059 (1976).
  * Note (3) The Eckart conditions fail to invert the quasicentroid unit when the ring-polymer inverts, which necessitated the discarding of inverting trajectories. Since the timescale of this inversion process is higher than any of the vibration timescales, this procedure ought to have a minimal effect on the spectrum.
  * Yurchenko  _et al._ (2005) S. N. Yurchenko, J. Zheng, H. Lin, P. Jensen, and W. Thiel, J. Chem. Phys. 123, 134308 (2005).
  * Richardson and Althorpe (2009) J. O. Richardson and S. C. Althorpe, J. Chem. Phys. 131, 214106 (2009).
  * Note (4) Note that $\rho$ defined here and the inversion coordinate $\overline{\rho}$ in ref.30 are related as $\rho=\frac{\pi}{2}-\overline{\rho}$.
  * Note (5) Strictly, the quasicentroid constraint does not prevent the instanton geometry from forming, but assigns it to a quasicentroid ensemble that corresponds to a bent geometry, where its contribution to the free energy is negligible.
  * Liu and Zhang (2016) J. Liu and Z. Zhang, J. Chem. Phys. 144, 34307 (2016).
  * Miller (2001) W. H. Miller, J. Phys. Chem. A 105, 2942 (2001).
  * Conte, Aspuru-Guzik, and Ceotto (2013) R. Conte, A. Aspuru-Guzik, and M. Ceotto, J. Phys. Chem. Lett. 4, 3407 (2013).
  * Yu and Bowman (2019) Q. Yu and J. M. Bowman, J. Phys. Chem. A 123, 1399 (2019).