# Nuclear Quantum Effects in the Acetylene:Ammonia Plastic Co-crystal

Atul C. Thakur  Richard C. Remsing<EMAIL_ADDRESS>Department of
Chemistry and Chemical Biology, Rutgers University, Piscataway, NJ 08854

###### Abstract

Organic molecular solids can exhibit rich phase diagrams. In addition to
structurally unique phases, translational and rotational degrees of freedom
can melt at different state points, giving rise to partially disordered solid
phases. The structural and dynamic disorder in these materials can have a
significant impact on the physical properties of the organic solid,
necessitating a thorough understanding of disorder at the atomic scale. When
these disordered phases form at low temperatures, especially in crystals with
light nuclei, the prediction of materials properties can be complicated by the
importance of nuclear quantum effects. As an example, we investigate nuclear
quantum effects on the structure and dynamics of the orientationally-
disordered, translationally-ordered plastic phase of the acetylene:ammonia
(1:1) co-crystal that is expected to exist on the surface of Saturn’s moon
Titan. Titan’s low surface temperature ($\sim 90$ K) suggests that the quantum
mechanical behavior of nuclei may be important in this and other molecular
solids in these environments. By using neural network potentials combined with
ring polymer molecular dynamics simulations, we show that nuclear quantum
effects increase orientational disorder and rotational dynamics within the
acetylene:ammonia (1:1) co-crystal by weakening hydrogen bonds. Our results
suggest that nuclear quantum effects are important to accurately model
molecular solids and their physical properties in low temperature
environments.

##  I Introduction

Molecular crystals play important roles in a wide variety of fields, including
pharmaceuticals [1, 2], semiconductor technologies [3, 4, 5], energy storage
and transport [6, 7, 8], pesticides [9], and atmospheric and planetary
sciences [10, 11, 12, 13, 14]. Like the more familiar class of atomic
crystals, molecular crystals are crystalline solids made of repeating
molecular subunits. When these solids are made of two or more types of
molecules, they are typically called molecular co-crystals [15, 2, 12].
Molecular co-crystals often exhibit thermodynamic and mechanical properties
that are different from the corresponding single component solids, offering
unique opportunities for tunable materials [16, 17, 18].

Molecular crystals are usually held together by weak interactions like
hydrogen bonding and dispersion interactions. As a result, molecular crystals
often display rich phase behavior. In addition to the usual translationally
and orientationally ordered crystalline solid phase, translationally ordered
and orientationally disordered (or vice versa) phases can also appear. Solids
that exhibit translational order but lack orientational order are known as
plastic crystals [19, 20, 21, 22]. The plasticity, as suggested by the name,
refers to the additional mechanical flexibility endowed by the presence of
orientational disorder. Plastic phases can also exhibit different chemical and
physical properties than the corresponding orientationally ordered phases. For
instance, some plastic phases exhibit different thermodynamic and electrical
properties than the ordered phase [23]. Importantly, these differences in
physical properties can lead to significant macroscopic consequences.

Many molecular co-crystals are predicted to form on the cold ($\sim$90 K)
surface of Saturn’s moon Titan [12, 13, 14, 11, 24, 25, 26, 27, 14, 28]. As a
result, molecular crystals and co-crystals on Titan play an analogous role to
minerals on Earth and have been termed cryominerals. On Titan and similar
planetary bodies, the changes in physical properties arising from plastic
crystal phases could have implications for geochemical processes and even
surface geology [10, 12, 13, 14]. To begin to understand plastic phases of
these cryominerals, we use molecular dynamics (MD) simulations to investigate
the acetylene:ammonia (1:1) plastic co-crystal [10, 12, 25]. In particular,
the quantum mechanical nature of atomic nuclei is expected to be important at
the cold temperatures on Titan’s surface, and we focus on quantifying the
impact of nuclear quantum effects on the structure and dynamics of the
acetylene:ammonia (1:1) plastic co-crystal.

Nuclear quantum effects (NQEs) on structure and dynamics are significant for
light nuclei like hydrogen even at room temperature (300 K) [29, 30, 31]. As a
result, the common approximation of using classical nuclei in molecular
simulations may be inadequate for many molecular crystals. Indeed, NQEs can
play important roles in determining thermodynamic stability [32, 33], thermal
conductivity [32], expansivity, and mechanical properties [34]. At the
cryogenic conditions on Titan, NQEs are anticipated to play an even greater
role in determining the properties of molecular co-crystals.

Here, we quantify the role of NQEs in determining the structure and dynamics
of the acetylene:ammonia (1:1) co-crystal using neural network potential-based
thermostatted ring polymer MD simulations. We first discuss simulation details
and the theory of approximate real time dynamics within the ring polymer
molecular dynamics framework. Then, we discuss the translational and
orientational structure of the acetylene:ammonia (1:1) co-crystal as obtained
from the classical and quantum simulations, where classical and quantum refer
to the treatment of nuclear degrees of freedom, and show that NQEs increase
rotational disorder within the co-crystal. We then quantify NQEs on dynamics
within the co-crystal through examination of vibrational densities of states
and rotational time correlation functions, and we connect the observed
increase in rotational diffusion by NQEs to a similar increase in hydrogen
bond dynamics and a lowering of rotational free energy barriers. Our results
suggest that NQEs can be significant in plastic crystals, especially at low
temperatures, and should be taken into account when modeling these disordered
solids.

##  II Simulation Details

We performed neural network (NN)-based MD simulations of the acetylene:ammonia
(1:1) co-crystal using the LAMMPS software package [35]. We trained a neural
equivariant interatomic potential (NequIP) [36] on ab initio density
functional theory (DFT) data generated using CP2K [37, 38, 39]. We also
trained a deep potential (DeePMD) model [40, 41, 42, 43] on a larger data set,
but this potential resulted in inaccuracies in dynamical properties, as
discussed further in the Appendix. The training data was obtained from the
DFT-based ab initio molecular dynamics (AIMD) simulations [44] at 90 K and 600
K. The ab initio energy and forces were computed using the Perdew-Burke-
Ernzerhof (PBE) treatment of exchange-correlation (XC) interactions [45] along
with Grimme’s D3 correction for long range dispersion interactions [46, 47].
We used molecularly optimized (MOLOPT) Gödecker-Teter-Hutter (GTH)
triple-$\zeta$ single polarization (TZVP-MOLOPT-GTH) basis sets along with
Gödecker-Teter-Hutter (GTH) relativistic, norm-conserving, separable, dual-
space pseudopotentials [48, 49, 50]. More details on the DFT simulations can
be found in our earlier work [10].

The NN potential was constructed using the NequIP code to optimize the force
and energy loss terms according to [36]

$\displaystyle\mathcal{L}=\lambda_{E}\left|\left|\hat{E}-E\right|\right|^{2}+\frac{\lambda_{F}}{3N}\sum_{i=1}^{N}\sum_{\alpha=1}^{3}\left|\left|-\frac{\partial\hat{E}}{\partial
r_{i,\alpha}}-F_{i,\alpha}\right|\right|,$ (1)

where $N$ refers to the total number of atoms in the system, $\alpha$
indicates each independent Cartesian direction, $\hat{E}$ is the predicted
potential energy and $\lambda$ refers to the weights for the energies
($\lambda_{E}$) and forces ($\lambda_{F}$). We fixed the maximum rotation
order for the network’s features at $l=1$ and allowed the features with odd
parity to be included with the network parity turned on. The locality for the
features was imposed with a radial cutoff of 6 Å. We uniformly sampled the DFT
dataset to select 2000 configurations, half from simulations at 90 K and half
from simulations at 600 K. From the total data set, a validation set of 400
snapshots was picked randomly. The energy data was internally shifted for each
particle with per atom total energy mean while forces were rescaled with force
root mean square deviation. Within the network, we used four interaction
blocks with an initial learning rate of 0.005 and a batch size of 5. The
learning rate gets rescaled by a factor of 0.5 if the validation loss of
atomic forces has not improved over 100 consecutive epochs. We terminated the
training process after around 4000 epochs because the errors in the forces and
energies were sufficiently small to give good agreement with the AIMD
simulations. The training process continuously updates the model with the best
force validation loss and we froze the model with the overall best validation
loss to be used for MD simulations.

We performed classical MD simulations of the acetylene:ammonia (1:1) co-
crystal using the above-described NequIP NN potential. All simulations were
performed using LAMMPS with the nequip pair style provided by the pair-nequip
library [36]. The crystal structure for the acetylene:ammonia (1:1) co-crystal
was taken from the experimental crystal structure reported by Boese  _et al._
(Cambridge Structural Database ID: FOZHOS) [51]. We used a $4\times 4\times 4$
supercell consisting of 1024 atoms; 128 acetylene and 128 ammonia molecules.
The system was equilibrated at 90 K for at least 500 ps in the canonical (NVT)
ensemble using a Nosé-Hoover thermostat with three thermostat chains [52, 53,
54, 55]. We also performed simulations with classical nuclei at 30 K, where
the crystal is not plastic and is orientationally ordered, to contrast with
the plastic phase and highlight the effects of orientational disorder. For
simulations at 30 K, the temperature was scaled down gradually starting from
an equilibrated configuration at 90 K using Nosé-Hoover chains. To gather
accurate dynamics, we sampled the system in the microcanoncal (NVE) ensemble
for at least 1.5 ns at 90 K and for 500 ps at 30 K. The equations of motion
were integrated using the velocity Verlet algorithm as implemented in LAMMPS
with a timestep of 1 fs.

To quantify NQEs, we performed NN potential-based thermostatted ring polymer
molecular dynamics (TRPMD) simulations using LAMMPS interfaced with i-PI [56].
The i-PI package exploits a client-server interface in which the integration
of the equations of motion for the ring polymer is handled by i-PI while the
computationally intensive energy and force calculations are offloaded to
client codes. We used the LAMMPS client interface along with the NequIP neural
network potential for evaluation of energies and forces. Convergence of the
ring polymer discretization with respect to Trotter bead number was evaluated
using NN-based simulations of the acetylene:ammonia (1:1) co-crystal and
achieved at 128 beads at 90 K, in agreement with previous simulations of
hydrocarbons at similar temperatures [57]. To avoid numerical instabilities
and enable the use of a larger timestep of 0.5 fs in the TRPMD simulations, we
used the Cayley integrator with the BCOCB factorization scheme [58, 59, 60].

Our TRPMD [61, 62] simulations employed the same $4\times 4\times 4$
simulation supercell as the classical simulations. The starting configuration
for the TRPMD simulations was chosen randomly from the classical simulation
and was further equilibrated at 90 K in the canonical (NVT) ensemble with a
global stochastic velocity rescaling (SVR) thermostat [63] for over 2 ps. To
ensure ergodicity in the TRPMD simulations, we shoot two trajectories, each
starting from a randomly selected state point along the NVT RPMD trajectory.
Both TRPMD trajectories were sampled for 25 ps. The internal ring polymer
modes were thermostatted during sampling using the white noise PILE-G
thermostat, while all of the physical modes of the system remained
unthermostatted. This procedure is expected to be accurate for extracting
approximate quantum dynamics while still mitigating ergodicity issues arising
in completely unthermostatted RPMD simulations [61, 62]. The results presented
here reflect the average over both TRPMD simulations.

##  III Real time Dynamics with TRPMD

We begin with a brief review of the definitions of various quantum time
correlations functions and their corresponding RPMD approximations [62]. We
use the Kubo-transformed definition of time correlation functions as our exact
quantum correlations which we approximate by time correlations from RPMD
simulations. However, we note that Kubo transforms are not the only definition
of exact quantum time correlations and symmetrized quantum correlations have
also been proposed that are related to the Kubo-transformed definition [64,
62].

Translational dynamics can be probed using the classical velocity
autocorrelation function

$C_{v}(t)=\frac{\left<\mathbf{v}(0)\cdot\mathbf{v}(t)\right>}{\left<\mathbf{v}(0)\cdot\mathbf{v}(0)\right>},$
(2)

where $\mathbf{v}(t)$ indicates the instantaneous velocity of an atom at time
$t$ and $\left<\cdots\right>$ indicates an ensemble average. The vibrational
power spectrum, $C_{v}(\omega)$, is the Fourier transform of $C_{v}(t)$. The
quantum mechanical Kubo-transformed velocity autocorrelation function is [65]

${C}_{v}(t)=\frac{1}{\beta}\int_{0}^{\beta}d\lambda\left<\mathbf{v}(-i\lambda\hbar)\cdot\mathbf{v}(t)\right>,$
(3)

where $\beta=1/k_{\rm B}T$ is the inverse Boltzmann temperature and
$\mathbf{v}(t)$ refers to the Heisenberg evolved velocity operator of an atom
at time $t$. The Kubo-transformed correlation function is the most classical
version of an exact quantum mechanical correlation function that one can
construct, and it obeys the same symmetries as the classical time correlation
functions [64, 62]. To compute quantum time correlation functions, we use the
RPMD approximation to the Kubo-transformed time correlation function. The RPMD
approximation to the Kubo-transformed velocity autocorrelation function is
given by [65]

$\displaystyle{C}_{v}(t)$
$\displaystyle\simeq\frac{1}{(2\pi\hbar)^{9NP}\mathcal{Z}_{P}}\int\int\prod_{k=1}^{3N}\prod_{\alpha=1}^{P}d\mathbf{p}_{k}^{(\alpha)}d\mathbf{q}_{k}^{(\alpha)}$
(4) $\displaystyle\times
e^{-\beta_{P}\mathcal{H}_{P}(\\{\mathbf{p}_{k}^{(\alpha)}\\},\\{\mathbf{q}_{k}^{(\alpha)}\\})}\mathbf{\overline{v}}(0)\cdot\mathbf{\overline{v}}(t),$

where $\beta_{P}=\beta/P$, $\mathcal{H}_{P}$ is the ring polymer Hamiltonian,
$\mathbf{q}^{(\alpha)}$ and $\mathbf{p}^{(\alpha)}$ refer to the position and
momentum coordinates of bead index $\alpha$, and $\mathcal{Z}_{P}$ is the
canonical partition function of the path integral ring polymer system with $P$
beads [65]. The operator $\mathbf{\overline{v}}(t)$ is the instantaneous bead-
averaged velocity operator of an atom at time $t$ given by

$\mathbf{\overline{v}}(t)=\frac{1}{P}\sum\limits_{\alpha=1}^{P}\mathbf{v}^{(\alpha)}(t).$
(5)

Note that the RPMD approximation to the Kubo-transformed quantum correlation
function becomes exact in the short-time limit ($t\rightarrow 0$) and
coincides with the classical correlation function in the limit that the ring-
polymer collapses onto a single bead ($P\rightarrow 1$).

The rotational dynamics of molecules can be quantified using rotational (or
orientational) time correlation functions [66]. For classical nuclei, the
rotational time correlation functions we examine are of the form

$C_{l}(t)=\left<P_{l}(\mathbf{\hat{n}}(0)\cdot\mathbf{\hat{n}}(t))\right>,$
(6)

where $\mathbf{\hat{n}}(t)$ is the instantaneous unit vector pointing along
the N-H bond for ammonia molecules or the C-C bond for acetylene molecules,
and $P_{l}(x)$ is the $l$th Legendre polynomial. We quantify rotational motion
with $C_{1}(t)$ and $C_{2}(t)$, which correspond to the first and second order
Legendre polynomials. These correlation functions can be experimentally
measured through IR and Raman spectroscopies, respectively [67, 68, 69]. The
classical expression for $C_{1}(t)$ is

$C_{1}(t)=\left<\mathbf{\hat{n}}(0)\cdot\mathbf{\hat{n}}(t)\right>,$ (7)

and $C_{2}(t)$ is given by

$C_{2}(t)=\left<\frac{3}{2}(\mathbf{\hat{n}}(0)\cdot\mathbf{\hat{n}}(t))^{2}-\frac{1}{2}\right>.$
(8)

The corresponding generalized Kubo-transformed expression for the $l$th order
quantum rotational time correlation function is [65]

${C}_{l}(t)=\int_{0}^{\beta}d\lambda\left<P_{l}(\mathbf{\hat{n}}(-i\lambda\hbar)\cdot\mathbf{\hat{n}}(t))\right>.$
(9)

Therefore, the $l=1$ Kubo-transformed rotational time correlation function is

${C}_{1}(t)=\int_{0}^{\beta}d\lambda\left<\mathbf{\hat{n}}(-i\lambda\hbar)\cdot\mathbf{\hat{n}}(t)\right>,$
(10)

and the corresponding $l=2$ expression is

${C}_{2}(t)=\int_{0}^{\beta}d\lambda\left<\frac{3}{2}(\mathbf{\hat{n}}(-i\lambda\hbar)\cdot\mathbf{\hat{n}}(t))^{2}-\frac{1}{2}\right>.$
(11)

The RPMD approximation to the Kubo-transformed $l=1$ rotational correlation
can then be written as [65, 70]

$\displaystyle{C}_{1}(t)$
$\displaystyle\simeq\frac{1}{(2\pi\hbar)^{9NP}\mathcal{Z}_{P}}\int\int\prod_{k=1}^{3N}\prod_{\alpha=1}^{P}d\mathbf{p}_{k}^{(\alpha)}d\mathbf{q}_{k}^{(\alpha)}$
(12) $\displaystyle\times
e^{-\beta_{P}\mathcal{H}_{P}(\\{\mathbf{p}_{k}^{(\alpha)}\\},\\{\mathbf{q}_{k}^{(\alpha)}\\})}\mathbf{\overline{n}}(0)\cdot\mathbf{\overline{n}}(t),$

where $\mathbf{\overline{n}}$ is a bead-averaged operator,

$\mathbf{\overline{n}}(t)=\frac{1}{P}\sum\limits_{\alpha=1}^{P}\mathbf{\hat{n}}^{(\alpha)}(t).$
(13)

The RPMD approximation to the $l=2$ rotational time correlation function is
given by [65, 70]

$\displaystyle{C}_{2}(t)$
$\displaystyle\simeq\frac{3}{2(2\pi\hbar)^{9NP}\mathcal{Z}_{P}}\int\int\prod_{k=1}^{3N}\prod_{\alpha=1}^{P}d\mathbf{p}_{k}^{(\alpha)}d\mathbf{q}_{k}^{(\alpha)}$
(14) $\displaystyle\times
e^{-\beta_{P}\mathcal{H}_{P}(\\{\mathbf{p}_{k}^{(\alpha)}\\},\\{\mathbf{q}_{k}^{(\alpha)}\\})}\operatorname{Tr}[\mathbf{\overline{M}}(0)\mathbf{\overline{M}}(t)]-\frac{1}{2},$

where $\operatorname{Tr}$ indicates the trace of a matrix and
$\mathbf{\overline{M}}(t)$ is a bead-averaged third rank tensor,

$\mathbf{\overline{M}}(t)=\frac{1}{P}\sum\limits_{\alpha=1}^{P}\mathbf{\hat{n}}^{(\alpha)}(t)\mathbf{\hat{n}}^{(\alpha)}(t)^{T}.$
(15)

The operator $\mathbf{\hat{n}}^{(\alpha)}(t)$ is a $3\times 1$ dimensional
vector at real time $t$ and at imaginary time (or bead index) $\alpha$.

Figure 1:  Snapshots showing the zig-zag structure of the acetylene:ammonia
(1:1) co-crystal from (a) classical NN based simulations and (b) quantum TRPMD
simulations, where every fifth bead in the ring polymer is shown. Carbon atoms
are colored gold, nitrogen atoms are colored black, and hydrogen atoms are
colored white. Both C-H$\cdots$N and N-H$\cdots\pi$ hydrogen bonds are shown
and colored according to the respective donor atoms. Both panels depict views
from the crystal $z$-axis.

We also quantify hydrogen bond dynamics with the time correlation function
[71, 72, 73]

$\tilde{C}_{\rm{HB}}(t)=\left<h(0)h(t)\right>,$ (16)

where the hydrogen bond indicator function is given by
$h(t)=\Theta[R_{c}-R(t)]\Theta[\phi_{c}-\phi(t)]$, $\Theta(x)$ is the
Heaviside step function, $R(t)$ is the distance between the two heavy atoms
involved in the H-bond at time $t$, and $\phi(t)$ is the donor-acceptor-H
angle at time $t$. The indicator function, $h(t)$, equals 1 if the H-bond
distance and angle are less than the respective cutoffs $R_{c}$ and
$\phi_{c}$. We use $R_{c}=4$ Å and $\phi_{c}=30^{\circ}$ for both types of
H-bonds. The Kubo-transformed hydrogen bond time correlation function is

$\tilde{C}_{\rm
HB}(t)=\frac{1}{\beta}\int_{0}^{\beta}d\lambda\left<h(-i\lambda\hbar){h}(t)\right>.$
(17)

Then, the RPMD approximation to the quantum hydrogen bond correlation function
is given by

$\displaystyle\tilde{C}_{\text{HB}}(t)$
$\displaystyle\simeq\frac{1}{(2\pi\hbar)^{9NP}\mathcal{Z}_{P}}\int\int\prod_{k=1}^{3N}\prod_{\alpha=1}^{P}d\mathbf{p}_{k}^{(\alpha)}d\mathbf{q}_{k}^{(\alpha)}$
(18) $\displaystyle\times
e^{-\beta_{P}\mathcal{H}_{P}(\\{\mathbf{p}_{k}^{(\alpha)}\\},\\{\mathbf{q}_{k}^{(\alpha)}\\})}\overline{h}(0)\overline{h}(t),$

where $\overline{h}(t)$ is a bead-averaged indicator function,

$\overline{h}(t)=\frac{1}{P}\sum\limits_{\alpha=1}^{P}\overline{h}^{(\alpha)}(t).$
(19)

We normalize the hydrogen bond time correlation functions according to

${C}_{\text{HB}}(t)=\frac{\tilde{C}_{\text{HB}}(t)}{\tilde{C}_{\text{HB}}(0)}.$
(20)

In general, $\tilde{C}_{\rm HB}(0)$ is smaller in a quantum system than a
classical one. The reduction of $\tilde{C}_{\rm HB}(t)$ at $t=0$ results from
the diffuse nature of the quantum particle, such that some beads are H-bonded
while others are not. Normalization of the H-bonding correlation function
allows us to better compare H-bond dynamics in the quantum and classical
systems.

##  IV Results and Discussion

###  IV.1 Nuclear Quantum Effects on Translational and Orientational Structure

The acetylene:ammonia (1:1) co-crystal consists of antiparallel planes of
acetylene and ammonia molecules arranged in a zig-zag fashion [10, 25, 74, 51,
12], as shown in Fig. 1. The zig-zag arrangement of the molecules within the
co-crystal is primarily held together by two types of interactions: non-
directional van der Waals interactions and directional hydrogen bonding
interactions. The acetylene:ammonia (1:1) co-crystal exhibits C-H$\cdots$N and
N-H$\cdots\pi$ hydrogen bonds [10, 25, 74, 51]. C-H$\cdots$N hydrogen bonds
result from electrostatic attractions between acetylene hydrogens and lone
pairs of nitrogen atoms, and N-H$\cdots\pi$ hydrogen bonds result from
electrostatic attractions between ammonia hydrogens and $\pi$ systems of
acetylene molecules. N-H$\cdots\pi$ hydrogen bonds are weaker than
C-H$\cdots$N hydrogen bonds. Thermal fluctuations at 90 K can disrupt the
N-H$\cdots\pi$ hydrogen bonds, causing the ammonia molecules to be transiently
free [10]. Additionally, the ammonia molecules within the co-crystal are
situated at positions of higher symmetry than the $C_{3}$ symmetry inherent to
the ammonia molecules. An ammonia molecule cannot hydrogen bond to all six
neighbors at once, and the presence of symmetry equivalent orientations can
lead to orientational disorder.

Figure 2:  Site-site radial distribution functions (RDFs), $g(r)$, computed
for (a) acetylene carbons and (b) ammonia nitrogens within the classical
(NNPMD) acetylene:ammonia (1:1) co-crystal at $T=30$ and $90$ K. The results
with the quantum nuclei at $T=90$ K are also shown (TRPMD).

We quantify the translational structure of the crystal by computing site-site
radial distribution functions, $g_{\rm XY}(r)$, where ${\rm X}$ and ${\rm Y}$
indicate atomic or electronic sites. The classical C-C and N-N RDFs evaluated
at 30 K are much sharper than the corresponding classical structures at 90 K,
Fig. 2. This is expected due to thermal broadening at higher temperatures. At
90 K, NQEs have little impact on the C-C and N-N RDFs aside from a slight
broadening of the peaks, Fig. 2. This suggests that NQEs do not significantly
impact the structure of heavy atoms in the acetylene:ammonia (1:1) co-crystal
at 90 K.

Figure 3:  Site-site radial distribution functions (RDFs), $g(r)$, relevant
for sites forming (a, c) C-H$\cdots$N hydrogen bonds and (b, d) N-H$\cdots\pi$
hydrogen bonds within the classical (NNPMD) acetylene:ammonia (1:1) co-crystal
at $T=30$ and $90$ K. The RDFs computed with quantum nuclei at $90$ K at are
also shown (TRPMD).

We probe the hydrogen bonding structure within the co-crystal by computing the
RDFs for the sites involved in hydrogen bonding interactions, Fig. 3. For
C-H$\cdots$N hydrogen bonds, we compute the C-N and (C)H-N RDFs, where the
notation (C)H-N indicates that the RDFs were computed between hydrogens of the
acetylene molecules and the nitrogen atoms on the ammonia. Analogously, for
N-H$\cdots\pi$ hydrogen bonds, we compute the N-$\pi$ and (N)H-$\pi$
structures, where the notation (N)H indicates the hydrogen atoms of the
ammonia. We define the $\pi$ site as the center of the C-C bond [10]. The 30 K
RDFs are sharper than those at 90 K, suggesting more structured hydrogen bonds
within the crystal at low temperatures. Upon increasing the temperature from
30 to 90 K, the first three peaks in the (N)H-$\pi$ RDF broaden, hinting at
the plastic crystal phase transition of the co-crystal. The C-N structure
within the co-crystal is minimally perturbed with the inclusion of NQEs.
However, the (C)H-N RDF slightly broadens upon inclusion of NQEs. Similarly,
the N-$\pi$ RDF remains relatively unaffected by NQEs, while the (N)H-$\pi$
RDF is slightly broadened. The peak broadening upon inclusion of NQEs suggests
that the quantum hydrogen bond structure is more disordered than the
corresponding classical network. This effect is primarily due to the hydrogen
nuclei, which experience the largest NQEs due to their light mass and large
thermal radius.

Figure 4:  Probability distributions of the ammonia H atom positions projected
onto the crystal $xz$ plane at (a) $T=30$ K and (b,c) $T=90$ K. Results for
both (b) classical and (c) quantum nuclei are shown at 90 K. The contour line
is drawn at a probability of 0.2. Peak heights are larger than one but the
color scales are zoomed in to highlight the connectivity of the distributions
at 90 K.

RDFs provide little information on the orientational order of the molecules
within the co-crystal, although some information on the orientational order
can be gleaned from the (N)H-$\pi$ correlations. We can better understand the
orientational order of the molecules by projecting their positions on the
crystallographic $xz$ plane, Fig. 4. The orientational distributions at 30 K
show the presence of discrete disconnected regions of finite probability,
indicating orientational order. In contrast, the continuous orientational
probability distribution of the classical system at 90 K shows that ammonia
molecules can adopt all angular positions in the $xz$ plane, albeit some
orientations are occupied more than others. These preferred regions correspond
to hydrogen bonded configurations with six distinct regions that arise from
two different orientations of ammonia that hydrogen bond to neighboring
acetylenes; these two orientations are discussed in more detail below. When
NQEs are included, the probability distribution further broadens along the
$xz$ plane and the ammonia molecules preferential orientation in specific
directions is decreased. The broadening of the probability distribution
suggests that rotational disorder of ammonia molecules is even more pronounced
in the quantum system. The continuous nature of the orientational
distributions at 90 K suggests that ammonia molecules dynamically transition
between different hydrogen bonded configurations and that NQEs increase these
transitions. We quantify the classical and quantum dynamics associated with
these transitions in the next section.

###  IV.2 Nuclear Quantum Effects on Translational and Rotational Dynamics

Figure 5:  Total all-atom vibrational power spectrum evaluated for
acetylene:ammonia (1:1) co-crystal at (a) $T=30$ K and (b) $T=90$ K with
(TRPMD) and without (NNPMD) nuclear quantum effects.

To further understand the increased rotational disorder from NQEs, we
investigate the dynamics within the crystal. We begin with a discussion of
NQEs on vibrations within the co-crystal as quantified by vibrational power
spectra, Fig. 5. The vibrational frequencies predicted are in good agreement
with experimental results, as we showed previously [10]. When heating the
system from 30 to 90 K, the peaks slightly broaden and increase in intensity,
reflecting increased fluctuations.

The inclusion of NQEs significantly impacts the vibrational densities of
states, especially the high frequency modes. The highest frequency peaks,
resulting from C-H and N-H stretching modes, are broadened into a single wide
peak. Broadening of high frequency vibrational bands can be attributed to
anharmonicity of the intermolecular vibrations, which becomes more important
when zero point contributions are included [75]. In general, all vibrational
peaks are red shifted and broadened by NQEs. The red shift with respect to the
classical spectrum brings the vibrational frequencies in better agreement with
experimental results [10, 25, 74]. Red shifting of vibrational bands by NQEs
also occurs in liquid water [76] and has been rationalized as softening of the
potential due to contributions from zero point fluctuations [77].

Figure 6:  Rotational time correlation functions, (a, b) $C_{1}(t)$ and (c, d)
$C_{2}(t)$, evaluated with (TRPMD) and without (NNPMD) NQEs for (a, c)
acetylene molecules and for (b, d) ammonia molecules within the
acetylene:ammonia (1:1) co-crystal. Results at both $T=30$ and 90 K are shown.
The dotted lines show the exponential fits to the correlation functions at
long times. The insets in panels (b, d) shows the short time decay of ammonia
rotational time correlation functions.

We now turn to investigating orientational dynamics within the
acetylene:ammonia (1:1) co-crystal. For the acetylene molecules within the co-
crystal, we compute $C_{1}(t)$ and $C_{2}(t)$ orientational time correlation
functions for the C-C bond vector, Fig. 6a,c. For the ammonia molecules, the
N-H bond vectors are used to define both orientational time correlation
functions, Fig. 6b,d.

For the acetylene molecules, both correlation functions show no appreciable
long time decay. In the short time limit, $t<500$ fs, the classical
correlations show a slight decay due to librational motion. This minor decay
is more apparent in the quantum correlation functions due to increased
librational motion from zero-point fluctuations [70]. Finally, in the long
time limit, both quantum and classical correlation functions saturate to a
constant value. This overall lack of decay indicates that acetylene molecules
do not undergo rotational diffusion and remain both translationally and
orientationally ordered.

We apply a similar analysis to understand the orientational dynamics of
ammonia molecules within the acetylene:ammonia (1:1) co-crystal. Unlike the
acetylene molecules, the $C_{1}(t)$ and $C_{2}(t)$ orientational correlations
for the ammonia molecules exhibit an exponential decay at long times,
indicating rotational diffusion. Following our previous work [10], we estimate
the correlation times $\tau_{1}$ and $\tau_{2}$ from $C_{1}(t)$ and $C_{2}(t)$
by numerically integrating the time correlation function from 0 to 500 fs and
integrating an exponential fit to the correlation function for longer times.
Inclusion of NQEs increases the rate of decay of these correlations and
shortens the correlation times; $\tau_{1}^{\rm c}=1.5$ ps while $\tau_{1}^{\rm
q}=0.79$ ps and $\tau_{2}^{\rm c}=1.3$ ps while $\tau_{2}^{\rm q}=0.52$ ps,
where c and q indicate classical and quantum correlation times, respectively.

At short times ($t<500$ fs), $C_{1}(t)$ and $C_{2}(t)$ both exhibit rapid
decay followed immediately by the onset of a characteristic glitch in the
correlation functions (inset Fig. 6b,d). A similar glitch is also observed in
the quantum time correlation functions and is suggestive of a librational
rebound in cage-like environments [65]. In the diffusive regime ($t>500$ fs),
both classical and quantum correlations decay exponentially as the molecules
eventually break through the cage overcoming the free energy barrier
associated with rotations. Overall, both classical and quantum correlation
functions suggest that the ammonia molecules within the crystal exhibit a
plastic phase, _i.e._ they remain translationally ordered but are
orientationally disordered.

Comparing the classical and quantum correlations for ammonia molecules at all
times indicates that quantum correlations decay faster than the classical
correlations. The rapid decay of quantum correlations is suggestive of
increased orientational disorder of ammonia molecules within RPMD simulations
and is a manifestation of NQEs. Another interesting NQE is the deviation of
the value of both $C_{1}(t)$ and $C_{2}(t)$ at $t=0$ from unity. Such
deviations from unity arise from quantum dispersions in the orientational
coordinates [65]. We note that the RPMD approximation coincides with the exact
Kubo-transformed quantum correlation function in the limit $t\rightarrow 0$,
implying that quantum dispersions within the orientational coordinates is a
purely quantum mechanical effect and not a consequence of the RPMD
approximation.

###  IV.3 Nuclear Quantum Effects on Hydrogen Bond Dynamics

To understand how NQEs increase the rotational disorder of ammonia molecules,
we first quantify their impact on hydrogen bond dynamics within the
acetylene:ammonia (1:1) co-crystal, Fig. 7. In hydrogen bonded systems, NQEs
can promote or obstruct the diffusional motion of molecules by strengthening
or weakening the hydrogen bond network [30, 75]. For instance, quantum
delocalization of the nuclei along the hydrogen bonding direction results in
stronger hydrogen bonds and slower dynamics while delocalization along the
perpendicular directions results in weaker hydrogen bonds and faster dynamics
[30, 75].

Figure 7:  Hydrogen bonding time correlation functions evaluated at $T=30$ and
90 K for (a) C-H$\cdots$N and (b) N-H$\cdots\pi$ hydrogen bonds within the
classical acetylene:ammonia (1:1) co-crystal (NNPMD). The hydrogen bond
dynamics with quantum nuclei at $T=90$ K are also shown (TRPMD). Dotted lines
correspond to exponential fits to the correlation functions at long times.

We probe NQEs on both types of hydrogen bonds in the acetylene:ammonia (1:1)
co-crystal. Classical hydrogen bonding time correlation functions at 30 and 90
K exhibit similar behavior for C-H$\cdots$N hydrogen bonds. The C-H$\cdots$N
hydrogen bond correlations undergo a slight decay at short times followed by
constant correlations. The quantum correlation functions show a similar decay,
although the magnitude of the short-time decay is slightly larger due to
increased positional disorder of the quantum nuclei. The overall lack of decay
in the classical and quantum hydrogen bonding correlations suggest that
C-H$\cdots$N hydrogen bonds do not undergo significant breakage and
reformation. This lack of hydrogen bond dynamics is consistent with the lack
of acetylene rotations.

For the N-H$\cdots\pi$ hydrogen bonds, both classical and quantum correlations
decay on a finite timescale and resemble the decay of the N-H bond
orientational correlations, suggesting that N-H$\cdots\pi$ hydrogen bonds
break and reform as ammonia molecules rotate. With quantum nuclei, the
correlations of the N-H$\cdots\pi$ hydrogen bonds decay faster than the
classical correlations. Fitting to an exponential decay yields a classical
H-bond correlation time of 0.75 ps and a quantum correlation time of 0.28 ps.
The faster hydrogen bond dynamics in the quantum system can be tied to the
faster reorientational motion of the ammonia molecules as a consequence of
NQEs.

###  IV.4 Nuclear Quantum Effects on Rotational Thermodynamics

Figure 8:  Total free energy surface averaged over all ammonia molecules. Both
classical (NNPMD) and quantum (TRPMD) free energy surfaces at $T=90$ K are
shown. The dotted line indicates the thermal energy, $k_{\rm B}T$.

We can rationalize the faster reorientational dynamics of ammonia molecules in
the quantum system by quantifying NQEs on thermodynamics. NQEs can affect
dynamics by decreasing relevant free energy barriers or nuclear tunneling
through free energy barriers. We define a one dimensional free energy, $\Delta
G(\phi_{z})$, describing the rotations of ammonia as

$\Delta G(\phi_{z})=-k_{\rm B}T\ln P(\phi_{z}),$ (21)

where $P(\phi_{z})$ is the probability distribution of the azimuthal angle
made by the N-H bond vector with the crystal $z$ axes, shown in Fig. 8 The
free energy exhibits six minima and barriers when one might anticipate only
three due to the symmetry of the ammonia molecule.

The six minima in the free energy profile arise from the ability of the
ammonia molecules to adopt two hydrogen bonding orientations related to each
other by a $60^{\circ}$ or $180^{\circ}$ rotation around the principle axis
and a small rotation around a perpendicular axis, Fig. 9c,d. We label these
orientations as Type I and Type II, although apart from the rotation angle
these two conformations are completely equivalent. At 30 K, when ammonia
molecules do not rotate, the free energy profile of the two conformations can
be separated, Fig. 9. The resulting free energy profiles for Type I and Type
II ammonias are essentially the same but offset by 90∘. At 90 K, however, the
ammonia molecules can transition between the two orientations and a
distinction between the two types cannot be made. As a result, $\Delta
G(\phi_{z})$ shows features for both types. The six-fold symmetry in the total
free energies in Fig. 8 is therefore a result of two free energy surfaces with
three-fold symmetry, one for Type I and one for Type II.

Figure 9:  Classical free energy surface as computed for the two types of
orientations, (a) Type I and (b) Type II, of ammonia molecules within the
acetylene:ammonia (1:1) co-crystal at the two temperatures studied here.
Snapshots of (c) Type I and (d) Type II orientations are also shown, with
N-H$\cdots\pi$ hydrogen bonds shown as black dashed lines.

Now that we understand the origins of the six-fold symmetry in
$\Delta~{}G(\phi_{z})$, we can compare the classical and quantum free
energies, Fig. 8. In the classical systems, the free energy barriers to
rotation are slightly larger than the thermal energy. The ammonia molecules
transition between hydrogen bonding configurations through a jump diffusion
process, mainly with 120∘ jumps, but also less probable 60∘ jumps that
transition the ammonia molecules between Type I and Type II. Upon inclusion of
NQEs, the barriers to rotation are significantly lowered, making the barriers
easier to cross. As a result, this lowering of free energy barriers by NQEs
increases the dynamic rotational disorder in the quantum system.

We note that the transition between Type I and Type II orientations requires
rotation around an axis perpendicular to $z$, and the ammonia rotations in
general are likely coupled to positional fluctuations of neighboring acetylene
molecules [23]. As a result, $\phi_{z}$ is not the only order parameter needed
to fully describe the free energy landscape governing ammonia rotations, and
the free energies presented here may not be sufficient to estimate kinetics.
However, $\Delta G(\phi_{z})$ adequately illustrates NQEs on free energy
barriers to rotation; barriers are lowered by NQEs. NQEs will similarly lower
free energy barriers on higher dimensional free energy landscapes involving
additional order parameters.

##  V Conclusion

In this work, we used classical and quantum molecular dynamics simulations
with neural network potentials to quantify NQEs on the structure and dynamics
of the acetylene:ammonia (1:1) plastic co-crystal at Titan surface conditions.
We find that NQEs increase structural and dynamic disorder within the co-
crystal by weakening the hydrogen bond network, especially those hydrogen
bonds formed between ammonia molecules and the $\pi$-system of the acetylene
molecules. This weakening of the hydrogen bond network is supported by shorter
hydrogen bond lifetimes and lower ammonia rotational free energy barriers in
the system with NQEs included. Because dynamic rotational disorder often
softens the mechanical properties of plastic crystals, e.g. by reducing
relevant elastic moduli, we expect that NQEs will alter the mechanical
response of plastic crystals on Titan’s surface. In addition, the weakening of
hydrogen bonds holding the crystal together by NQEs may also promote processes
like dissolution at liquid-solid interfaces. As a result, we expect that NQEs
may need to be taken into account when predicting the properties of Titan
cryominerals, especially those that may exist in a plastic crystal phase.

###### Acknowledgements.

This work is supported by the National Aeronautics and Space Administration
under grant number 80NSSC20K0609, issued through the NASA Exobiology Program.
We acknowledge the Office of Advanced Research Computing (OARC) at Rutgers,
The State University of New Jersey for providing access to the Amarel and
Caliburn clusters and associated research computing resources that have
contributed to the results reported here. This work used the Advanced
Cyberinfrastructure Coordination Ecosystem: Services & Support (ACCESS),
formerly Extreme Science and Engineering Discovery Environment (XSEDE) [78],
which is supported by National Science Foundation grant number ACI-1548562.
Specifically, this work used Stampede2 and Ranch at the Texas Advanced
Computing Center through allocation TG-CHE210081. We thank Harender
Dhattarwal, Josiah Hall, Amanda Liyanaarachchi, and Olivia Martin for feedback
on the manuscript.

##  VI Appendix A: Convergence of Training Neural Network Potentials and
Comparison with Ab Initio Simulations

Figure 10:  Learning curves for the NequIP and DeePMD models showing the
convergence of validation errors in (a) forces and (b) energies per atom on a
log-log plot. NequIP errors are showing the mean absolute errors while the
DeePMD curves are plotted using root mean square errors.

We assess the accuracy of the NequIP and DeePMD NN potentials by computing the
learning curves and training errors, as well as comparing the models against
ab initio molecular dynamics (AIMD) simulations. Note that the NequIP model
was trained on 2000 configurations sampled from AIMD simulations at 90 K and
600 K, while the DeePMD model was trained on a larger data set consisting of
9000 configurations obtained from simulations at 30 K, 90 K, 200 K, and 600 K.
The errors in the energies and forces, quantified as the mean absolute errors
for NequIP and root mean square errors for DeePMD, are shown in Fig. 10. The
energy errors are converged to a value of $9.38\times 10^{-5}~{}{\rm eV/atom}$
for NequIP and $4.37\times 10^{-5}~{}{\rm eV/atom}$ for DeePMD, while the
force errors respectively are $8.24\times 10^{-3}~{}{\rm eV/\AA}$ and
$2.52\times 10^{-2}~{}{\rm eV/\AA}$, which are expected to yield reasonably
accurate results [36, 79, 80].

Figure 11:  Comparison of the site-site radial distribution functions (RDFs),
$g(r)$, for sites within the acetylene:ammonia (1:1) co-crystal, as obtained
from ab initio simulations and from simulations using the NequIP and DeePMD NN
models.  Figure 12:  Comparison of the $C_{2}(t)$ rotational correlations for
(a) acetylene and (b) ammonia molecules as obtained from simulations using ab
initio and NequIP and DeePMD NN potential at $T=90$ K. Hydrogen bond time
correlation functions for both (c) C-H$\cdots$N and (d) N-H$\cdots\pi$ type
hydrogen bonds are also shown.

We also compare the structure produced by the NN potentials to that from the
reference AIMD simulations. The partial structure for the carbon and nitrogen
sites within the acetylene:ammonia (1:1) computed using both NN potentials
compares well with the predicted structure from AIMD simulations, suggesting
that both NN potentials adequately describe the potential energy surface of
the acetylene:ammonia (1:1) co-crystal, Fig. 11a,b. The RDFs relevant to
C-H$\cdots$N and N-H$\cdots\pi$ hydrogen bonding are also reproduced well by
both NN potentials, Fig. 11c,d. However, DeePMD produces a slightly higher
first minimum in the (N)H-$\pi$ RDF. The agreement between results from AIMD
and NN potentials suggests that the NN models can describe the energetics of
hydrogen bonding in the co-crystal.

As an even stricter test reflecting the accuracy of the NN potential, we
compare the rotational and hydrogen bond dynamics of the ammonia molecules
predicted by the NN potentials to AIMD results. The rotational and hydrogen
bond time correlation functions involving acetylene and C-H$\cdots$N hydrogen
bonds are reproduced well by both NN potentials, Fig. 12a,c. The NequIP
potential also accurately reproduces the rotational and hydrogen bond time
correlation functions, Fig. 12b,d, further suggesting that NequIP can
accurately reproduce the energetics of hydrogen bonding within the co-crystal.
In contrast, the DeePMD potential produces rotational and hydrogen bond
dynamics that are too fast. This suggests that the DeePMD potential produces
N-H$\cdots\pi$ hydrogen bonds that are too weak, resulting in classical
dynamics that are similar to those produced by RPMD dynamics, despite
producing a more structured (N)H-$\pi$ than the quantum system. The reason for
the disagreement between the DeePMD and AIMD simulations is not immediately
apparent, but recent work on water suggests that deficiencies in describing
many body effects within DeePMD can lead to such disagreements [81].

## References

  * Beran [2016] G. J. Beran, “Modeling polymorphic molecular crystals with electronic structure theory,” Chem. Rev. 116, 5567–5613 (2016).
  * Kumar and Nanda [2018] S. Kumar and A. Nanda, “Approaches to design of pharmaceutical cocrystals: A review,” Molecular Crystals and Liquid Crystals 667, 54–77 (2018).
  * Coropceanu  _et al._ [2007] V. Coropceanu, J. Cornil, D. A. da Silva Filho, Y. Olivier, R. Silbey, and J.-L. Brédas, “Charge transport in organic semiconductors,” Chem. Rev. 107, 926–952 (2007).
  * Zhang  _et al._ [2018a] Z. Zhang, R. C. Remsing, H. Chakraborty, W. Gao, G. Yuan, M. L. Klein, and S. Ren, “Light-induced dilation in nanosheets of charge-transfer complexes,” PNAS 115, 3776–3781 (2018a).
  * Bronstein  _et al._ [2020] H. Bronstein, C. B. Nielsen, B. C. Schroeder, and I. McCulloch, “The role of chemical design in the performance of organic semiconductors,” Nat. Rev. Chem. 4, 66–77 (2020).
  * Hirshberg, Gerber, and Krylov [2014] B. Hirshberg, R. B. Gerber, and A. I. Krylov, “Calculations predict a stable molecular crystal of N8,” Nat. Chem. 6, 52–56 (2014).
  * Zhang  _et al._ [2020] H. Zhang, Y. Fang, F. Yang, X. Liu, and X. Lu, “Aromatic organic molecular crystal with enhanced $\pi$–$\pi$ stacking interaction for ultrafast Zn-ion storage,” Energy. Environ. Sci. 13, 2515–2523 (2020).
  * Liu  _et al._ [2020] S. Liu, L. Zhao, M. Yao, M. Miao, and B. Liu, “Novel All-nitrogen Molecular Crystals of Aromatic N10,” Adv. Sci. 7, 1902320 (2020).
  * Zhu  _et al._ [2021] X. Zhu, C. T. Hu, B. Erriah, L. Vogt-Maranto, J. Yang, Y. Yang, M. Qiu, N. Fellah, M. E. Tuckerman, M. D. Ward, _et al._ , “Imidacloprid crystal polymorphs for disease vector control and pollinator protection,” J. Am. Chem. Soc. 143, 17144–17152 (2021).
  * Thakur and Remsing [2023] A. C. Thakur and R. C. Remsing, “Molecular Structure, Dynamics, and Vibrational Spectroscopy of the Acetylene: Ammonia (1: 1) Plastic Co-crystal at Titan Conditions,” ACS Earth Space Chem. 7, 479–489 (2023).
  * Czaplinski  _et al._ [2023] E. C. Czaplinski, T. H. Vu, M. L. Cable, M. Choukroun, M. J. Malaska, and R. Hodyss, “Experimental Characterization of the Pyridine: Acetylene Co-crystal and Implications for Titan’s Surface,” ACS Earth Space Chem. 7, 597–608 (2023).
  * Cable  _et al._ [2021] M. L. Cable, T. Runčevski, H. E. Maynard-Casely, T. H. Vu, and R. Hodyss, “Titan in a Test Tube: Organic Co-crystals and Implications for Titan Mineralogy,” Acc. Chem. Res. 54, 3050–3059 (2021).
  * Maynard-Casely  _et al._ [2018] H. E. Maynard-Casely, M. L. Cable, M. J. Malaska, T. H. Vu, M. Choukroun, and R. Hodyss, “Prospects for Mineralogy on Titan,” Am. Mineral. 103, 343–349 (2018).
  * Maynard-Casely  _et al._ [2016] H. E. Maynard-Casely, R. Hodyss, M. L. Cable, T. H. Vu, and M. Rahm, “A Co-crystal between Benzene and Ethane: A Potential Evaporite Material for Saturn’s Moon Titan,” IUCrJ 3, 192–199 (2016).
  * Kirchner, Bläser, and Boese [2010] M. T. Kirchner, D. Bläser, and R. Boese, “Co-crystals with Acetylene: Small Is Not Simple!” Chem.-Eur. J. 16, 2131–2146 (2010).
  * Sun  _et al._ [2019] L. Sun, Y. Wang, F. Yang, X. Zhang, and W. Hu, “Cocrystal engineering: a collaborative strategy toward functional materials,” Adv Mater. 31, 1902328 (2019).
  * Karki  _et al._ [2009] S. Karki, T. Friščić, L. Fábián, P. R. Laity, G. M. Day, and W. Jones, “Improving mechanical properties of crystalline solids by cocrystal formation: New compressible forms of paracetamol,” Adv Mater. 21, 3905–3909 (2009).
  * Dai, Chen, and Lu [2018] X.-L. Dai, J.-M. Chen, and T.-B. Lu, “Pharmaceutical cocrystallization: An effective approach to modulate the physicochemical properties of solid-state drugs,” CrystEngComm 20, 5292–5316 (2018).
  * Klein and Lewis [1990] M. L. Klein and L. J. Lewis, “Simulation of Dynamical Processes in Molecular Solids,” Chem. Rev. 90, 459–479 (1990).
  * Klein [1985] M. L. Klein, “Computer simulation studies of solids,” Annu. Rev. Phys. Chem. 36, 525–548 (1985).
  * Bounds, Klein, and Patey [1980] D. G. Bounds, M. L. Klein, and G. Patey, “Molecular Dynamics Simulation of the Plastic Phase of Solid Methane,” J. Chem. Phys. 72, 5348–5356 (1980).
  * Nosé and Klein [1983] S. Nosé and M. L. Klein, “Structural Transformations in Solid Nitrogen at High Pressure,” Phys. Rev. Lett. 50, 1207 (1983).
  * Lynden-Bell and Michel [1994] R. Lynden-Bell and K. Michel, “Translation-rotation Coupling, Phase Transitions, and Elastic Phenomena in Orientationally Disordered Crystals,” Rev. Modern Phys. 66, 721 (1994).
  * Cable  _et al._ [2019] M. L. Cable, T. H. Vu, M. J. Malaska, H. E. Maynard-Casely, M. Choukroun, and R. Hodyss, “A Co-crystal between Acetylene and Butane: A Potentially Ubiquitous Molecular Mineral on Titan,” ACS Earth Space Chem. 3, 2808–2815 (2019).
  * Cable  _et al._ [2018] M. L. Cable, T. H. Vu, H. E. Maynard-Casely, M. Choukroun, and R. Hodyss, “The Acetylene-Ammonia Co-crystal on Titan,” ACS Earth Space Chem. 2, 366–375 (2018).
  * Cable  _et al._ [2020] M. L. Cable, T. H. Vu, M. J. Malaska, H. E. Maynard-Casely, M. Choukroun, and R. Hodyss, “Properties and Behavior of the Acetonitrile-Acetylene Co-crystal under Titan Surface Conditions,” ACS Earth Space Chem. 4, 1375–1385 (2020).
  * Vu  _et al._ [2014] T. H. Vu, M. L. Cable, M. Choukroun, R. Hodyss, and P. Beauchamp, “Formation of a New Benzene-ethane Co-crystalline Structure under Cryogenic Conditions,” J. Phys. Chem. A 118, 4087–4094 (2014).
  * Vu  _et al._ [2022] T. H. Vu, H. E. Maynard-Casely, M. L. Cable, M. Choukroun, M. J. Malaska, and R. Hodyss, “1, 3-Butadiene on Titan: Crystal Structure, Thermal Expansivity, and Raman Signatures,” ACS Earth Space Chem. (2022).
  * Markland and Ceriotti [2018] T. E. Markland and M. Ceriotti, “Nuclear quantum effects enter the mainstream,” Nat. Rev. Chem. 2, 0109 (2018).
  * Ceriotti  _et al._ [2016] M. Ceriotti, W. Fang, P. G. Kusalik, R. H. McKenzie, A. Michaelides, M. A. Morales, and T. E. Markland, “Nuclear quantum effects in water and aqueous systems: Experiment, theory, and current challenges,” Chem. Rev. 116, 7529–7550 (2016).
  * Remsing [2023] R. C. Remsing, “Modeling Nuclear Quantum Effects on Long Range Electrostatics in Nonuniform Fluids,” (2023), arXiv:2309.04613 [physics.chem-ph] .
  * Rossi, Gasparotto, and Ceriotti [2016] M. Rossi, P. Gasparotto, and M. Ceriotti, “Anharmonic and Quantum Fluctuations in Molecular Crystals: A First-principles Study of the Stability of Paracetamol,” Phys. Rev. Lett. 117, 115702 (2016).
  * Kapil and Engel [2022] V. Kapil and E. A. Engel, “A complete description of thermodynamic stabilities of molecular crystals,” PNAS 119, e2111769119 (2022).
  * Ko  _et al._ [2018] H.-Y. Ko, R. A. DiStasio Jr, B. Santra, and R. Car, “Thermal Expansion in Dispersion-bound Molecular Crystals,” Phys. Rev. Mater. 2, 055603 (2018).
  * Thompson  _et al._ [2022] A. P. Thompson, H. M. Aktulga, R. Berger, D. S. Bolintineanu, W. M. Brown, P. S. Crozier, P. J. in’t Veld, A. Kohlmeyer, S. G. Moore, T. D. Nguyen, _et al._ , “LAMMPS- A Flexible Simulation Tool for Particle-based Materials Modeling at the Atomic, Meso, and Continuum Scales,” Comput. Phys. Commun. 271, 108171 (2022).
  * Batzner  _et al._ [2022] S. Batzner, A. Musaelian, L. Sun, M. Geiger, J. P. Mailoa, M. Kornbluth, N. Molinari, T. E. Smidt, and B. Kozinsky, “E (3)-equivariant graph neural networks for data-efficient and accurate interatomic potentials,” Nat. Commun. 13, 2453 (2022).
  * Kühne  _et al._ [2020] T. D. Kühne, M. Iannuzzi, M. Del Ben, V. V. Rybkin, P. Seewald, F. Stein, T. Laino, R. Z. Khaliullin, O. Schütt, F. Schiffmann, D. Golze, J. Wilhelm, S. Chulkov, M. H. Bani-Hashemian, V. Weber, U. Borštnik, M. Taillefumier, A. S. Jakobovits, A. Lazzaro, H. Pabst, T. Müller, R. Schade, M. Guidon, S. Andermatt, N. Holmberg, G. K. Schenter, A. Hehn, A. Bussy, F. Belleflamme, G. Tabacchi, A. Glöß, M. Lass, I. Bethune, C. J. Mundy, C. Plessl, M. Watkins, J. VandeVondele, M. Krack, and J. Hutter, “CP2K: An Electronic Structure and Molecular Dynamics Software Package-Quickstep: Efficient and Accurate Electronic Structure Calculations,” J. Chem. Phys. 152, 194103 (2020).
  * Hutter  _et al._ [2014] J. Hutter, M. Iannuzzi, F. Schiffmann, and J. VandeVondele, “CP2K: Atomistic Simulations of Condensed Matter Systems,” WIREs Comput. Mol. Sci. 4, 15–25 (2014).
  * VandeVondele and Hutter [2003] J. VandeVondele and J. Hutter, “An Efficient Orbital Transformation Method for Electronic Structure Calculations,” J. Chem. Phys. 118, 4365–4369 (2003).
  * Zhang  _et al._ [2018b] L. Zhang, J. Han, H. Wang, W. Saidi, R. Car, and W. E, “End-to-end symmetry preserving inter-atomic potential energy model for finite and extended systems,” in  _Advances in Neural Information Processing Systems_, Vol. 31, edited by S. Bengio, H. Wallach, H. Larochelle, K. Grauman, N. Cesa-Bianchi, and R. Garnett (Curran Associates, Inc., 2018).
  * Lu  _et al._ [2021] D. Lu, H. Wang, M. Chen, L. Lin, R. Car, E. Weinan, W. Jia, and L. Zhang, “86 PFLOPS Deep Potential Molecular Dynamics Simulation of 100 Million Atoms with Ab Initio Accuracy,” Comput. Phys. Commun. 259, 107624 (2021).
  * Zeng  _et al._ [2023] J. Zeng, D. Zhang, D. Lu, P. Mo, Z. Li, Y. Chen, M. Rynik, L. Huang, Z. Li, S. Shi, _et al._ , “DeePMD-kit v2: A Software Package for Deep Potential Models,” J. Chem. Phys. 159 (2023).
  * Wang  _et al._ [2018] H. Wang, L. Zhang, J. Han, and E. Weinan, “DeePMD-kit: A Deep Learning Package for Many-body Potential Energy Representation and Molecular Dynamics,” Comput. Phys. Commun. 228, 178–184 (2018).
  * VandeVondele  _et al._ [2005] J. VandeVondele, M. Krack, F. Mohamed, M. Parrinello, T. Chassaing, and J. Hutter, “QUICKSTEP: Fast and Accurate Density Functional Calculations Using a Mixed Gaussian and Plane Waves Approach,” Comput. Phys. Commun. 167, 103–128 (2005).
  * Perdew, Burke, and Ernzerhof [1996] J. P. Perdew, K. Burke, and M. Ernzerhof, “Generalized Gradient Approximation Made Simple,” Phys. Rev. Lett. 77, 3865 (1996).
  * Grimme  _et al._ [2010] S. Grimme, J. Antony, S. Ehrlich, and H. Krieg, “A Consistent and Accurate Ab Initio Parametrization of Density Functional Dispersion Correction (DFT-D) for the 94 Elements H-Pu,” J. Chem. Phys. 132, 154104 (2010).
  * Grimme, Ehrlich, and Goerigk [2011] S. Grimme, S. Ehrlich, and L. Goerigk, “Effect of the Damping Function in Dispersion Corrected Density Functional Theory,” J. Comput. Chem. 32, 1456–1465 (2011).
  * Krack [2005] M. Krack, “Pseudopotentials for H to Kr Optimized for Gradient-corrected Exchange-correlation Functionals,” Theor. Chem. Acc. 114, 145–152 (2005).
  * Hartwigsen, Gœdecker, and Hutter [1998] C. Hartwigsen, S. Gœdecker, and J. Hutter, “Relativistic Separable Dual-space Gaussian Pseudopotentials from H to Rn,” Phys. Rev. B 58, 3641 (1998).
  * Goedecker, Teter, and Hutter [1996] S. Goedecker, M. Teter, and J. Hutter, “Separable Dual-space Gaussian Pseudopotentials,” Phys. Rev. B 54, 1703 (1996).
  * Boese, Bläser, and Jansen [2009] R. Boese, D. Bläser, and G. Jansen, “Synthesis and Theoretical Characterization of an Acetylene-Ammonia Cocrystal,” J. Am. Chem. Soc. 131, 2104–2106 (2009).
  * Nosé [1984a] S. Nosé, “A Unified Formulation of the Constant Temperature Molecular Dynamics Methods,” J. Chem. Phys. 81, 511–519 (1984a).
  * Nosé [1984b] S. Nosé, “A molecular dynamics method for simulations in the canonical ensemble,” Mol. Phys. 52, 255–268 (1984b).
  * Hoover [1985] W. G. Hoover, “Canonical Dynamics: Equilibrium Phase-space Distributions,” Phys. Rev. A 31, 1695 (1985).
  * Martyna, Klein, and Tuckerman [1992] G. J. Martyna, M. L. Klein, and M. Tuckerman, “Nosé–hoover Chains: The Canonical Ensemble Via Continuous Dynamics,” J. Chem. Phys. 97, 2635–2643 (1992).
  * Kapil  _et al._ [2019] V. Kapil, M. Rossi, O. Marsalek, R. Petraglia, Y. Litman, T. Spura, B. Cheng, A. Cuzzocrea, R. H. Meißner, D. M. Wilkins, _et al._ , “i-PI 2.0: A Universal Force Engine for Advanced Molecular Simulations,” Comput. Phys. Commun. 236, 214–223 (2019).
  * Uhl, Marx, and Ceriotti [2016] F. Uhl, D. Marx, and M. Ceriotti, “Accelerated Path Integral Methods for Atomistic Simulations at Ultra-low Temperatures,” J. Chem. Phys. 145 (2016).
  * Korol, Bou-Rabee, and Miller [2019] R. Korol, N. Bou-Rabee, and T. F. Miller, “Cayley Modification for Strongly Stable Path-integral and Ring-polymer Molecular Dynamics,” J. Chem. Phys. 151 (2019).
  * Rosa-Raíces  _et al._ [2021] J. L. Rosa-Raíces, J. Sun, N. Bou-Rabee, and T. F. Miller, “A Generalized Class of Strongly Stable and Dimension-free T-RPMD Integrators,” J. Chem. Phys. 154 (2021).
  * Korol  _et al._ [2020] R. Korol, J. L. Rosa-Raíces, N. Bou-Rabee, and T. F. Miller, “Dimension-free Path-integral Molecular Dynamics without Preconditioning,” J. Chem. Phys. 152 (2020).
  * Rossi, Ceriotti, and Manolopoulos [2014] M. Rossi, M. Ceriotti, and D. E. Manolopoulos, “How to Remove the Spurious Resonances from Ring Polymer Molecular Dynamics,” J. Chem. Phys. 140 (2014).
  * Craig and Manolopoulos [2004] I. R. Craig and D. E. Manolopoulos, “Quantum Statistics and Classical Mechanics: Real Time Correlation Functions from Ring Polymer Molecular Dynamics,” J. Chem. Phys. 121, 3368–3373 (2004).
  * Bussi, Zykova-Timan, and Parrinello [2009] G. Bussi, T. Zykova-Timan, and M. Parrinello, “Isothermal-isobaric Molecular Dynamics Using Stochastic Velocity Rescaling,” J. Chem. Phys. 130 (2009).
  * Habershon  _et al._ [2013] S. Habershon, D. E. Manolopoulos, T. E. Markland, and T. F. Miller III, “Ring-polymer Molecular Dynamics: Quantum Effects in Chemical Dynamics from Classical Trajectories in an Extended Phase Space,” Annu. Rev. Phys. Chem. 64, 387–413 (2013).
  * Miller and Manolopoulos [2005] T. F. Miller and D. E. Manolopoulos, “Quantum Diffusion in Liquid Water from Ring Polymer Molecular Dynamics,” J. Chem. Phys. 123 (2005).
  * Hansen and McDonald [2006] J. P. Hansen and I. R. McDonald,  _Theory of Simple Liquids_ (Elsevier Ltd., 2006).
  * Thakur and Remsing [2021] A. C. Thakur and R. C. Remsing, “Distributed charge models of liquid methane and ethane for dielectric effects and solvation,” Mol. Phys. 119, e1933228 (2021).
  * Gordon [1965a] R. Gordon, “Molecular Motion in Infrared and Raman Spectra,” J. Chem. Phys. 43, 1307–1312 (1965a).
  * Gordon [1965b] R. Gordon, “Relations between Raman Spectroscopy and Nuclear Spin Relaxation,” J. Chem. Phys. 42, 3658–3665 (1965b).
  * Wilkins  _et al._ [2017] D. M. Wilkins, D. E. Manolopoulos, S. Pipolo, D. Laage, and J. T. Hynes, “Nuclear Quantum Effects in Water Reorientation and Hydrogen-bond Dynamics,” J. Phys. Chem. Lett. 8, 2602–2607 (2017).
  * Luzar and Chandler [1996a] A. Luzar and D. Chandler, “Effect of Environment on Hydrogen Bond Dynamics in Liquid Water,” Phys. Rev. Lett. 76, 928 (1996a).
  * Luzar and Chandler [1996b] A. Luzar and D. Chandler, “Hydrogen-bond Kinetics in Liquid Water,” Nature 379, 55–57 (1996b).
  * Luzar [2000] A. Luzar, “Resolving the Hydrogen Bond Dynamics Conundrum,” J. Chem. Phys. 113, 10663–10675 (2000).
  * Preston and Signorell [2012] T. C. Preston and R. Signorell, “The Formation and Stability of Co-crystalline NH3·C2H2 Aerosol Particles,” Mol. Phys. 110, 2807–2815 (2012).
  * Habershon, Markland, and Manolopoulos [2009] S. Habershon, T. E. Markland, and D. E. Manolopoulos, “Competing Quantum Effects in the Dynamics of a Flexible Water Model,” J. Chem. Phys. 131 (2009).
  * Paesani  _et al._ [2006] F. Paesani, W. Zhang, D. A. Case, T. E. Cheatham, and G. A. Voth, “An Accurate and Simple Quantum Model for Liquid Water,” J. Chem. Phys. 125 (2006).
  * Paesani, Iuchi, and Voth [2007] F. Paesani, S. Iuchi, and G. A. Voth, “Quantum Effects in Liquid Water from an Ab Initio-based Polarizable Force Field,” J. Chem. Phys. 127 (2007).
  * Towns  _et al._ [2014] J. Towns, T. Cockerill, M. Dahan, I. Foster, K. Gaither, A. Grimshaw, V. Hazlewood, S. Lathrop, D. Lifka, G. D. Peterson, R. Roskies, J. R. Scott, and N. Wilkins-Diehr, “XSEDE: Accelerating Scientific Discovery,” Comput. Sci. Eng. 16, 62–74 (2014).
  * Pagotto, Zhang, and Duignan [2022] J. Pagotto, J. Zhang, and T. Duignan, “Predicting the properties of salt water using neural network potentials and continuum solvent theory,” (2022).
  * Yue  _et al._ [2021] S. Yue, M. C. Muniz, M. F. Calegari Andrade, L. Zhang, R. Car, and A. Z. Panagiotopoulos, “When Do Short-range Atomistic Machine-learning Models Fall short?” J. Chem. Phys. 154 (2021).
  * Zhai  _et al._ [2023] Y. Zhai, A. Caruso, S. L. Bore, Z. Luo, and F. Paesani, “A “short Blanket” Dilemma for a State-of-the-art Neural Network Potential for Water: Reproducing Experimental Properties or the Physics of the Underlying Many-body interactions?” J. Chem. Phys. 158 (2023).