Om $#\mathrm{2}$

# Bayesian Estimation of the Hydroxyl Radical Diffusion Coefficient at Low
Temperature and High Pressure from Atomistic Molecular Dynamics

Carter T. Butts<EMAIL_ADDRESS><EMAIL_ADDRESS>Departments of Sociology,
Statistics, Computer Science, and EECS, University of California, Irvine
Rachel W. Martin  Departments of Chemistry and Molecular Biology and
Biochemistry; University of California, Irvine

(7/28/21)

###### Abstract

The hydroxyl radical is the primary reactive oxygen species produced by the
radiolysis of water, and is a significant source of radiation damage to living
organisms. Mobility of the hydroxyl radical at low temperatures and/or high
pressures is hence a potentially important factor in determining the
challenges facing psychrophilic and/or barophilic organisms in high-radiation
environments (e.g., ice-interface or undersea environments in which radiative
heating is a potential heat and energy source). Here, we estimate the
diffusion coefficient for the hydroxyl radical in aqueous solution, using a
hierarchical Bayesian model based on atomistic molecular dynamics trajectories
in TIP4P/2005 water over a range of temperatures and pressures.

††preprint: AIP/123-QED

Ionizing radiation is a feature of both terrestrial and extraterrestrial
environments, presenting challenges as well as opportunities for living
organisms. On one hand, radiation can damage biological molecules either
directly or via the production of reactive chemical species that modify or
degrade them (Blanco  _et al._ , 2018). On the other, ionizing radiation can
act as a non-photosynthetic energy source for microbial communities, in some
cases producing chemical species suitable for chemotrophy (Onstott  _et al._ ,
2003) or in other cases maintaining a habitable environment by temperature
elevation and/or maintenance of liquid water pockets via ice melting (Tarnas
_et al._ , 2018; Ojha  _et al._ , 2021; Tarnas  _et al._ , 2021). While some
terrestrial organisms are known to tolerate high levels of radiation either
transiently or on an ambient basis (White  _et al._ , 1999; Cavicchioli, 2002;
Munteanu, Uivarosi, and Andries, 2015), the role of ionizing radiation in
determining habitability in a broader biological context remains largely open.

For conventional microbial organisms in aqueous environments, the predominant
threat posed by ionizing radiation is the formation of reactive oxygen species
due to the radiolysis of water (LaVerne, 2000). Many species are produced,
including hydrogen peroxide, the superoxide radical, and the hydroperoxyl
radical; (Matheson, 1964) however, the hydroxyl radical ($\cdot\mathrm{OH}$)
is the dominant source of prompt radiation damage for solvated biomolecules
due to its high production rate, reactivity, and unsuitability for enzymatic
processing (Ghosal  _et al._ , 2005). In terrestrial mesophilic organisms
under typical cellular conditions, $\cdot\mathrm{OH}$ survives on sufficiently
long time scales to diffuse to and damage biological macromolecules (Roots and
Okada, 1975), prominently including proteins and DNA. Although the problems
associated with DNA damage are well-appreciated, proteins are the major
cellular targets of $\cdot\mathrm{OH}$-mediated radiation damage (Du and
Gebicki, 2004). Post-translational modification of proteins due to radical
interactions can lead to a wide array of potentially lethal consequences,
including formation of reactive peroxide species (Davies, Fu, and Dean, 1995),
formation of insoluble aggregates (Barnham, Masters, and Bush, 2004), loss of
enzymatic function, and destabilization of functional complexes; at minimum,
such damage increases the rate of protein expression and controlled
degradation required for homeostasis, thereby raising the metabolic cost of
cellular survival. Although $\cdot\mathrm{OH}$ production is a threat to
irradiated organisms in any environment, some environmental conditions may
partially or substantially ameliorate it. In particular, environments that
favor the scavenging of $\cdot\mathrm{OH}$ by other chemical species reduce
the level of chemical stress to which organisms are subject, and may thus
allow for greater radiation tolerance. Because many reactions involving
$\cdot\mathrm{OH}$ are (or are near) diffusion-limited kinetics, the
$\cdot\mathrm{OH}$ diffusion coefficient is of particular relevance to
extremophile biochemistry in irradiated environments.

Among the environments of particular interest for novel biochemistry are those
involving low temperatures and/or high pressures. Such environments occur in
deep ocean and ice/rock interfaces on Earth, and in subsurface oceans in the
outer solar system. Because the diffusion coefficient for $\cdot\mathrm{OH}$
at low temperature and high pressure has not been measured to date, we here
estimate it using atomistic molecular dynamics (MD) simulations, employing a
novel hierarchical Bayesian inference scheme to infer the diffusion
coefficient while correcting for finite size effects. In the process, we also
parameterize a CHARMM-compatible model for $\cdot\mathrm{OH}$ for use with the
TIP4P/2005 water model (chosen for its ease of implementation with standard MD
platforms and its performance in reproducing properties of bulk water over a
wide temperature and pressure range). We summarize our posterior inference in
the form of a simple log-log polynomial model that can be used to reproduce
our simulation-based estimates of $D_{\cdot\mathrm{OH}}$ over a range of
temperatures and pressures.

The remainder of the paper is structured as follows: Sec. I describes our
procedures, in particular including the parameterization of the
$\cdot\mathrm{OH}$ model (Sec. I.1) and inference for the diffusion
coefficient (Sec. I.3). Our results are summarized in Sec. II, and Sec. III
concludes the paper.

##  I Methods

Our interest is inferring the diffusion coefficient of the hydroxyl radical in
aqueous solution, $D_{\cdot\mathrm{OH}}$, as a function of temperature and
pressure. We begin with parameterization of a model for $\cdot\mathrm{OH}$ in
TIP4P/2005 water, followed by our simulation design. We then describe our
approach for inferring the diffusion coefficient, $D_{\cdot\mathrm{OH}}$, from
simulated water and radical trajectories. Results are shown in Sec. II.

###  I.1 Parameterization of the $\cdot\mathrm{OH}$ Model

To perform atomistic simulations of $\cdot\mathrm{OH}$ in solution, we
parameterize a CHARMM-based (Best  _et al._ , 2012) model for
$\cdot\mathrm{OH}$ in TIP4P/2005 water (Absacal and Vega, 2005). We employ
TIP4P/2005 because of its strong performance in reproducing the diffusion
constant of water over a wide range of temperatures and pressures
(Tsimpanogiannis  _et al._ , 2019). Partial charges, bond length, mass, and
force constant for $\cdot\mathrm{OH}$ are taken from Pabis, Szala-Bilnik, and
Swiatla-Wojcik (2011), who performed combined DFT and MD studies of
$\cdot\mathrm{OH}$ in BJH water at physiological temperature and pressure
(Bopp, Jancsó, and Heinzinger, 1983); unfortunately, the non-bonded
interactions employed cannot be directly adapted to the CHARMM forcefield, and
hence it is necessary to parameterize them directly. The non-bonded
interactions in question are defined by a Lennard-Jones potential of the form

$\begin{split}\epsilon_{LJ}(R_{AB}|\epsilon_{A},\epsilon_{B},r_{\min
A}{/\scriptstyle 2},r_{\min B}{/\scriptstyle 2})=\\\
\sqrt{\epsilon_{A}\epsilon_{B}}\left[\left(\frac{r_{\min A}{/\scriptstyle
2}+r_{\min B}{/\scriptstyle 2}}{R_{AB}}\right)^{12}+2\left(\frac{r_{\min
A}{/\scriptstyle 2}+r_{\min B}{/\scriptstyle
2}}{R_{AB}}\right)^{6}\right],\end{split}$

where $R_{AB}$ is the distance between atoms $A$ and $B$, $\epsilon_{A}$ and
$\epsilon_{B}$ are species-specific well-depth parameters, and $r_{\min
A}{/\scriptstyle 2}$ and $r_{\min B}{/\scriptstyle 2}$ are “half-radii” that
determine the zero-point of the interatomic force. Here, we must determine
these parameters for the two respective atoms of $\cdot\mathrm{OH}$, given the
TIP4P/2005 parameters (which we take as fixed). As our interest is in
$D_{\cdot\mathrm{OH}}$, we optimize the $\cdot\mathrm{OH}$ non-bonded
interaction parameters ($\epsilon_{\mathrm{O}}$, $\epsilon_{\mathrm{H}}$,
$r_{\min\mathrm{O}}{/\scriptstyle 2}$, and $r_{\min\mathrm{H}}{/\scriptstyle
2}$) so as to reproduce the measured value of $D_{\cdot\mathrm{OH}}$
=0.23Å2/ps for $\cdot\mathrm{OH}$ in water at 298K and atmospheric pressure
(Dorfman and Adams, 1973), holding all other factors constant.

Our protocol proceeded as follows. We began with a quasi-random search of the
parameter space, drawing 250 points from the intervals
$\epsilon_{\mathrm{O}},\epsilon_{\mathrm{H}}\in(-0.5,0)$ kcal/mol,
$r_{\min\mathrm{O}}{/\scriptstyle 2}\in(0,3)$ Å, and
$r_{\min\mathrm{H}}{/\scriptstyle 2}\in(0,2)$ Å using a four-dimensional
Halton sequence (bases 2, 3, 5, and 7). For each parameter vector, a 1 ns
atomistic simulation of one $\cdot\mathrm{OH}$ in TIP4P/2005 water under
periodic boundary conditions was performed, with frames sampled every 0.5 ps
(integrator step size 2 fs). Simulations were initialized with a cubic box of
$\approx 20$ Å side length at 1 atm and 298K; two adjustment phases of 100 ps
each were performed (with box sizes adjusted for PME calculations after each
phase) prior to the production run, with both adjustment and production phases
performed with the $NpT$ ensemble. Langevin dynamics with an interval of 1/ps
were employed for temperature control, and a Langevin-Nosé-Hoover piston with
a period of 100 ps was used to maintain constant pressure (Martyna, Tobias,
and Klein, 1994; Feller  _et al._ , 1995). Rigid bonds were maintained for all
waters, with the O-H bond of $\cdot\mathrm{OH}$ left flexible. All simulations
were performed using NAMD (Phillips  _et al._ , 2005), with initial conditions
created using VMD (Humphrey, Dalke, and Schulten, 1996), psfgen (Ribeiro  _et
al._ , 2020), and Packmol (Martínez  _et al._ , 2009). Each simulated
trajectory was then unwrapped using the protocol of von Bülow, Bullerjahn, and
Hummel (2020) to account for changing box sizes, and each frame was centered
at its centroid to correct for net drift. Molecular positions were extracted
via the resulting oxygen atom coordinates.

To obtain initial estimates for $D_{\cdot\mathrm{OH}}$ at each parameter
value, the covariance-based estimator of Bullerjahn, von Bülow, and Hummel
(2020) was applied to each processed $\cdot\mathrm{OH}$ trajectory; this
estimator is computationally efficient, and was found to work well in pilot
runs using both $\cdot\mathrm{OH}$ and $\mathrm{H}_{2}\mathrm{O}$ under these
simulation conditions. These raw estimates were then corrected for finite
sample sizes using the analytical correction factor of Yeh and Hummer (2004),

$D_{\infty}=D_{\mathrm{MD}}+\frac{\zeta k_{B}T}{2\pi\eta L},$

where $D_{\infty}$ is the diffusion constant at infinite size,
$D_{\mathrm{MD}}$ is the diffusion constant under PBC obtained from MD
simulation, $\zeta\approx 2.8373$ is a numerical constant, $\eta$ is the shear
viscosity of the solvent, and $L$ is the box length. For $\eta$, the
TIP4P/2005 viscosity of 8.55$\times 10^{-4}$ Js/m3 from González and Abascal
(2010) was employed, and $L$ was taken to be the cube root of the mean box
volume over the simulation. The size-corrected estimates of
$D_{\cdot\mathrm{OH}}$ were then retained for further analysis.

To obtain an initial estimate of the non-bonded parameters, the ARE in the
$D_{\cdot\mathrm{OH}}$ estimate (i.e.
$|0.23-\widehat{D_{\cdot\mathrm{OH}}}|/0.23$) was computed for each parameter
vector. Inspection of the resulting estimates revealed two regions of high
performance (Fig. S1). To distinguish among the competing regions, energies
were calculated for the interaction of a single $\cdot\mathrm{OH}$ and
$\mathrm{H}_{2}\mathrm{O}$ in both conventional hydrogen bonding and “flipped”
(i.e., unfavorable H-H or O-O configurations) over a range of 1-5 Å
(Figs.S2-S4). The single region for which properly oriented hydrogen bonds
were more favorable and the ARE was low was selected for further analysis.
(The most favorable point in this stage was
$\epsilon_{\mathrm{O}}\approx-0.224$, $\epsilon_{\mathrm{H}}\approx-0.291$,
$r_{\min\mathrm{O}}{/\scriptstyle 2}\approx 0.802$,
$r_{\min\mathrm{H}}{/\scriptstyle 2}\approx 0.752$.)

Following this first stage of calibration, we performed a subsequent stage of
refinement by repeating the procedure with search using larger simulated
systems over a smaller parameter range; this secondarily confirmed that
estimates were robust to system size. As before, a four-dimensional Halton
sequence was used to select parameters, with 25 draws taken over the range
$\epsilon_{\mathrm{O}}\in(-0.26,-0.18)$ kcal/mol,
$\epsilon_{\mathrm{H}}\in(-0.31,-0.27)$ kcal/mol,
$r_{\min\mathrm{O}}{/\scriptstyle 2}\in(0.6,1.1)$ Å, and
$r_{\min\mathrm{H}}{/\scriptstyle 2}\in(0.72,0.78)$ Å. For each parameter
vector, the above simulation and analysis protocol was followed, with the
exception that initial box size was increased to 50 Å. Estimated size-
corrected AREs were obtained for each parameter value, and the parameter
vector yielding the minimum error was selected for final use. (As above,
energy calculations were employed to verify that the model favored the correct
$\cdot\mathrm{OH}$-$\mathrm{H}_{2}\mathrm{O}$ orientation.)

The final parameters for the $\cdot\mathrm{OH}$ model (including both
predetermined and calibrated parameters) are $r_{\mathrm{OH}}=0.9751$ Å,
$E_{\mathrm{OH}}=475.6$ kcal/mol/Å2, $m_{\mathrm{O}}=15.994$ Da,
$m_{\mathrm{O}}=1.008$ Da, $\epsilon_{\mathrm{O}}=-0.235$ kcal/mol,
$\epsilon_{\mathrm{H}}=-0.3068$ kcal/mol, $r_{\min\mathrm{O}}{/\scriptstyle
2}=0.78185$ Å, and $r_{\min\mathrm{H}}{/\scriptstyle 2}=0.746939$ Å. The
relatively small value of $r_{\min\mathrm{O}}{/\scriptstyle 2}$ seems to be
necessary for proper interaction with the TIP4P/2005 water model, as larger
values lead either to inappropriate donor/acceptor orientations for inaccurate
values of $D_{\cdot\mathrm{OH}}$ (Figs. S3,S4); the selected value leads to
reasonable water/radical interactions, as shown in Fig. 1. The size-corrected
diffusion constant estimate for the final model in the last selection round
was $\widehat{D_{\cdot\mathrm{OH}}}\approx 0.222$ Å2/ps, for an ARE of
approximately 4% versus the target value of 0.23 Å2/ps.

Figure 1:  Non-bonded interaction energy for $\mathrm{H}_{2}\mathrm{O}$ vs.
$\cdot\mathrm{OH}$ as a function of configuration and distance between
proximate atoms. The final model correctly favors standard donor/acceptor
orientations over “flipped” orientations, with the $\mathrm{H}_{2}\mathrm{O}$
donor interaction slightly favored over the $\cdot\mathrm{OH}$ donor
interaction.

###  I.2 Diffusion Simulation

Given the parameterized model, we simulate the diffusion of $\cdot\mathrm{OH}$
in TIP4P/2005 water at multiple temperatures and pressures. The basic
simulation strategy is equivalent to that described in Section I.1, with the
following modifications. First, in order to ensure a high level of precision
in our trajectory calculations, we sample trajectories every 0.25 ps instead
of every 0.5 ps; additionally, we add an extra 100 ps to the pre-production
run for each trajectory, which is not included in the 1 ns employed for
estimation. Rather than employing the covariance estimator for
$D_{\cdot\mathrm{OH}}$, we use a Bayesian approach as described in Section
I.3.2. And, finally, because the analytical correction for system size depends
on the viscosity (which is not known for TIP4P/2005 over the range of
temperatures and pressures studied here), we employ statistical corrections
involving variable system sizes as explained below.

We perform simulations at 263, 273, 283, and 298K, and at 1, 10, 100, 1000,
and 10000 atm pressure (full factorial design). To allow statistical
correction for finite size effects (and to reduce simulation-related error) we
perform 30 replicate simulations for each condition, with initial box sizes
evenly spaced from 20 to 50 Å. The unfolded trajectories from each size
replicate are then used to infer the diffusion coefficient, as described
below.

###  I.3 Inference for the Diffusion Coefficient

To infer $D_{\cdot\mathrm{OH}}$ from simulation, we must account for both
transients that affect the observed within-trajectory diffusion rate, and
finite size effects that lead to systematic variation across trajectories for
systems of differing size. Here, we use a two-stage Bayesian inference
strategy, first obtaining local posterior estimates of $D_{\mathrm{MD}}$ using
a modified Brownian motion process model, and then integrating these local
estimates via a hierarchical model that combines estimates of
$D_{\mathrm{MD}}^{\cdot\mathrm{OH}}$ and
$D_{\mathrm{MD}}^{\mathrm{H}_{2}\mathrm{O}}$ across systems of varying size to
obtain final estimates of $D_{\cdot\mathrm{OH}}$.

####  I.3.1 Local Estimation of $D_{\mathrm{MD}}$

Although the covariance estimator of Bullerjahn, von Bülow, and Hummel (2020)
is both computationally and statistically efficient at high temperatures, it
performs less well when the diffusion coefficient becomes small relative to
the the background noise for which it controls (see discussion in Bullerjahn,
von Bülow, and Hummel (2020)). Here, we thus use a strategy of Bayesian
estimation for the local (size uncorrected) diffusion constant
$D_{\mathrm{MD}}$, which both makes more complete use of data and provides
regularization of the resulting estimator. The model employed here is based on
the model (somewhat tacitly) underlying the generalized least squares
estimators of Bullerjahn, von Bülow, and Hummel (2020), namely a latent
Brownian motion process with a Gaussian observation mechanism. Given regularly
spaced observations $Y=(Y_{1},Y_{2},\ldots)$ at times $1,2,\ldots$, the
process may be defined in one dimension by

$\displaystyle X_{t+1}=X_{t}+Z_{t}$ (1) $\displaystyle Y_{t+1}=X_{t+1}+W_{t},$
(2)

where $Z_{t}$ is iid $N(0,\sigma^{2})$, and $W_{t}$ is iid $N(0,a^{2})$.
Physically, $X$ here represents a “true” or idealized Brownian motion with
independent perturbations given by $Z$, while $W$ reflects idiosyncratic noise
factors arising from non-Brownian transients. The diffusion constant
corresponds to $D=\sigma^{2}/2$ (in the squared distance units of $Y$ divided
by the time between steps). We may observe that this is a (discretely
measured) Gaussian process with covariance function
$K(Y_{i},Y_{j})=\min(i,j)\sigma^{2}+a^{2}I(i=j)$, and hence the likelihood is
given by (conditioning on and centering the first observation)

$p(Y=y|\sigma^{2},a^{2})=MVN(y|\mathbf{0},K(y)),$ (3)

where $\mathbf{0}$ is the 0-vector, and $K(y)$ is the Gram matrix of the
observed sequence. This is straightforward to work with, although
computationally expensive when the number of time points becomes large, and
can be pooled across dimensions in the isotropic case (as is done here).

To define priors on the variance parameters, we first observe that on our
physical scale of interest (Å2/ps) and over the range of conditions considered
here, it is _a priori_ unlikely that $\sigma^{2}$ will exceed 0.5; likewise
(per Bullerjahn, von Bülow, and Hummel, 2020) it is reasonable to expect
$a^{2}$ to be comparable to (or possibly smaller than) $\sigma^{2}$. We thus
use independent half-Gaussian priors for $\sigma^{2}$ and $a^{2}$, with a
scale of 0.5, which is relatively flat over the region of interest while
discouraging strongly unphysical values. (Note that this is equivalent to
$L_{2}$ regularization of the variance parameters.)

For parameter estimation, it is natural here due to both the Gaussian
structure of the problem and the large data size to employ MAP estimation,
invoking the Laplace approximation (Gelman  _et al._ , 2003) to obtain
posterior standard deviations. We perform estimation using a custom R (R Core
Team, 2021) script, with direct optimization of the log posterior surface
using BFGS (Nash, 1990); the mclust package (Scrucca  _et al._ , 2016) was
used for efficient calculation of the multivariate Gaussian log-likelihood.
The posterior standard deviation of $\hat{D}$ was obtained via the Hessian of
the negative log posterior about the posterior mode (exploiting the linear
relationship between $D$ and $\sigma^{2}$). Due to the cost of computing the
Gram matrix, the initial unfolded and drift-corrected trajectories were
downsampled from 0.25 ps to 0.5 ps resolution, and split into two segments of
500 ps length (i.e., 1000 observations); these were pooled in the likelihood
calculation. In the case of $\mathrm{H}_{2}\mathrm{O}$ trajectories,
trajectories for all water molecules were pooled and jointly analyzed. This
process led to estimates of the posterior mode (assumed equal to the mean,
under the Laplace approximation) and standard deviation for $D_{\mathrm{MD}}$
for both $\cdot\mathrm{OH}$ and $\mathrm{H}_{2}\mathrm{O}$ under each
condition and at each box size. These estimates were then integrated to
estimate $D_{\mathrm{H}_{2}\mathrm{O}}$ and $D_{\cdot\mathrm{OH}}$ in each
condition as described below.

####  I.3.2 Estimation of $D_{\cdot\mathrm{OH}}$

Estimation of small-molecule diffusion constants is challenging due both to
the need to correct for finite-size effects and a high level of idiosyncratic
variation between trajectories that is difficult to account for; moreover, a
single molecule trajectory provides relatively little information per
simulation run (as opposed to the large number of solvent molecule
trajectories obtained on each run). Here, we address both issues via a
hierarchical Bayesian model that pools information between
$\mathrm{H}_{2}\mathrm{O}$ and $\cdot\mathrm{OH}$ trajectories, and that
incorporates multiple sources of variation. The model (whose structure is
described pictorially in the plate diagram of Figure 2) is defined as follows.

$L_{i}$$\alpha$$\widehat{D_{\mathrm{MD}}^{\mathrm{H}_{2}\mathrm{O}}}_{i},\widehat{S_{\mathrm{MD}}^{\mathrm{H}_{2}\mathrm{O}}}_{i}$$\widehat{D_{\mathrm{MD}}^{\cdot\mathrm{OH}}}_{i},\widehat{S_{\mathrm{MD}}^{\cdot\mathrm{OH}}}_{i}$$D_{\mathrm{H}_{2}\mathrm{O}}$$D_{\cdot\mathrm{OH}}$${\tau_{\mathrm{H}_{2}\mathrm{O}}}_{i}$${\tau_{\cdot\mathrm{OH}}}_{i}$$\gamma_{\mathrm{H}_{2}\mathrm{O}}$$\gamma_{\cdot\mathrm{OH}}$$\mu_{\mathrm{H}_{2}\mathrm{O}}$$\mu_{\cdot\mathrm{OH}}$$i\in
1,\ldots,N$ Figure 2:  Structure of the diffusion constant model. White
circles indicate latent quantities, while observed random quantities are
indicated by shaded circles; fixed values are uncircled. Elements within the
central plate are replicated $N$ times, while those outside are pooled.

We begin with the observation that, if $D_{\cdot\mathrm{OH}}$ and
$D_{\mathrm{H}_{2}\mathrm{O}}$ are the respective bulk diffusion constants for
$\cdot\mathrm{OH}$ and $\mathrm{H}_{2}\mathrm{O}$, then

$\displaystyle
D_{\cdot\mathrm{OH}}={D_{\mathrm{MD}}^{\cdot\mathrm{OH}}}_{i}-\alpha/L_{i}$
(4) $\displaystyle
D_{\mathrm{H}_{2}\mathrm{O}}={D_{\mathrm{MD}}^{\mathrm{H}_{2}\mathrm{O}}}_{i}-\alpha/L_{i},$
(5)

where ${D_{\mathrm{MD}}^{\cdot\mathrm{OH}}}_{i}$ and
${D_{\mathrm{MD}}^{\mathrm{H}_{2}\mathrm{O}}}_{i}$ are the diffusion
coefficients for a PBC system with length scale $L_{i}$, and $\alpha$ is the
system size scaling coefficient. Of these, only $L$ is observed. We do,
however, have local estimates of ${D_{\mathrm{MD}}^{\cdot\mathrm{OH}}}_{i}$
and ${D_{\mathrm{MD}}^{\mathrm{H}_{2}\mathrm{O}}}_{i}$, which we model as

$\displaystyle\widehat{{D_{\mathrm{MD}}^{\cdot\mathrm{OH}}}_{i}}={D_{\mathrm{MD}}^{\cdot\mathrm{OH}}}_{i}+E_{\cdot\mathrm{OH}}$
(6)
$\displaystyle\widehat{{D_{\mathrm{MD}}^{\mathrm{H}_{2}\mathrm{O}}}_{i}}={D_{\mathrm{MD}}^{\mathrm{H}_{2}\mathrm{O}}}_{i}+E_{\mathrm{H}_{2}\mathrm{O}},$
(7)

where $E_{\cdot\mathrm{OH}}\sim
N(0,\widehat{S_{\mathrm{MD}}^{\cdot\mathrm{OH}}}_{i}+{\tau_{\cdot\mathrm{OH}}}^{2}_{i})$
and $E_{\mathrm{H}_{2}\mathrm{O}}\sim
N(0,\widehat{S_{\mathrm{MD}}^{\mathrm{H}_{2}\mathrm{O}}}_{i}+{\tau_{\mathrm{H}_{2}\mathrm{O}}}^{2}_{i})$
represent deviations from the idealized local diffusion coefficients. The
error variance is modeled via two components: the posterior variance from the
local model of Section I.3.1
($\widehat{S_{\mathrm{MD}}^{\cdot\mathrm{OH}}}_{i},\widehat{S_{\mathrm{MD}}^{\mathrm{H}_{2}\mathrm{O}}}_{i}$);
and the excess variances ${\tau_{\cdot\mathrm{OH}}}^{2}_{i}$ and
${\tau_{\mathrm{H}_{2}\mathrm{O}}}^{2}_{i}$ representing trajectory-specific
idiosyncratic deviations not reflected by the within-trajectory estimates. We
take the square roots of the excess variances to be generated by 0-truncated
normal distributions, i.e. $\tau_{\cdot\mathrm{OH}}\sim
TN_{0}(\mu_{\cdot\mathrm{OH}},\gamma_{\cdot\mathrm{OH}}^{2}),\tau_{\mathrm{H}_{2}\mathrm{O}}\sim
TN_{0}(\mu_{\mathrm{H}_{2}\mathrm{O}},\gamma_{\mathrm{H}_{2}\mathrm{O}}^{2})$,
with weakly informative standard half-Cauchy priors on the $\mu_{*}$ and
$\gamma_{*}$ parameters. We observe that this prior structure can be seen as
flexibly generalizing several standard regression-like models: in the limit as
$E_{*}\to 0$, we recover a model akin to weighted least squares, with weights
based on the locally estimated variances; when $E_{*}\gg\hat{S_{*}}_{i}$ but
$\gamma_{*}\to 0$, recover a model akin to a standard homoskedastic
regression; and when $\gamma_{*}\gg 0$ we obtain a robust regression with a
heavy-tailed error distribution. Finally, we take $\alpha$ to be _a priori_
uniform on (0,0.75) (as the 0.75 is expected to be strictly larger than the
value of $\alpha$ for the conditions studied here).

Given the above, we perform posterior simulation using the No-U-Turn
Hamiltonian Monte Carlo algorithm (Homan and Gelman, 2014) from the Stan
library (Stan Development Team, 2020a, b). 4 chains were employed for each
condition, with $10^{5}$ burn-in iterations per chain followed by $10^{5}$
additional iterations from which 1000 were retained (i.e., a thinning interval
of 100) for a final sample size of 4000 draws per condition. Convergence was
assessed with $\hat{R}$ (Gelman and Rubin, 1992). Posterior means and 95%
posterior intervals were obtained for $D_{\cdot\mathrm{OH}}$ and
$D_{\mathrm{H}_{2}\mathrm{O}}$ for each condition for subsequent analysis, as
discussed below.

##  II Results and Discussion

Fig. 3 shows the estimated bulk diffusion coefficients for $\cdot\mathrm{OH}$
and $\mathrm{H}_{2}\mathrm{O}$ at the simulated temperatures and pressures.
Although $D_{\cdot\mathrm{OH}}$ and $D_{\mathrm{H}_{2}\mathrm{O}}$
approximately coincide under ambient conditions, we observe some differences
in their response to temperature and pressure. While TIP4P/2005 reproduces the
experimentally observed increase in $\mathrm{H}_{2}\mathrm{O}$ diffusion rate
under moderate pressure, $\cdot\mathrm{OH}$ diffusion slows with pressure at
moderate temperatures (eventually gaining very modest acceleration in the low
temperature regime). At pressures approaching 10000 atm, both water and the
hydroxyl radical slow considerably (as expected from the well-known increase
in water viscosity in this regime), with $D_{\cdot\mathrm{OH}}$ diffusing more
slowly than water despite its smaller size. This difference in the high-
pressure behavior of $\cdot\mathrm{OH}$ versus its behavior under ambient
conditions may be related to the observation of an MD study in BJH water at
relatively high temperature (310K) and 1 atm that $\cdot\mathrm{OH}$ tends to
occupy cavities in the $\mathrm{H}_{2}\mathrm{O}$ hydrogen bonding network
(Pabis, Szala-Bilnik, and Swiatla-Wojcik, 2011). At high pressure, occupancy
of such cavities will be highly favorable, and may retard diffusion.

Changes in $\cdot\mathrm{OH}$ diffusion may have implications for
radiochemistry in high pressure/low temperature environments. The direct
products of $\gamma$-irradiation of water are formed in “spurs,” localized
regions of high concentration for these highly reactive species (most
directly, $\cdot\mathrm{OH}$, $\cdot\mathrm{H}$, and hydrated electrons)
(Schwarz, 1969; Parajon  _et al._ , 2008). When a $\gamma$-ray interacts with
a water molecule, it generates a highly reactive excited state,

$\mathrm{H}_{2}\mathrm{O}+\gamma\rightarrow\mathrm{H}_{2}\mathrm{O}*$ (8)

which can either directly decay to form a hydrogen atom and a hydroxyl
radical,

$\mathrm{H}_{2}\mathrm{O}*\rightarrow\cdot\mathrm{H}+\cdot\mathrm{OH}$ (9)

or it can lose an electron to the surrounding solution, resulting in a radical
cation that in turn decomposes to yield a hydroxyl radical and a proton:

$\mathrm{H}_{2}\mathrm{O}*\rightarrow\mathrm{H}_{2}\mathrm{O}^{+}+e^{-}_{\mathrm{aq}}$
(10)

$\mathrm{H}_{2}\mathrm{O}^{+}\rightarrow\mathrm{H}^{+}+\cdot\mathrm{OH}$ (11)

These reactions serve as the starting point for the formation of a complex
mixture of reactive oxygen species Ershov and Gordeev (2008) and/or reactions
of these primary products with biomolecules (Omar, Hasnaoui, and de la Lande,
2021). The recent advent of experimental techniques enabling direct detection
of the long-hypothesized $\mathrm{H}_{2}\mathrm{O}^{+}$ intermediate (Loh  _et
al._ , 2020) and observation of attosecond dynamics in liquid water (Jordan
_et al._ , 2020), as well as the discovery of new reaction pathways (Thürmer
_et al._ , 2013; Ren  _et al._ , 2018) have led to renewed interest in the
fundamental radiochemistry of aqueous solutions. Analysis of the relevant
reactions requires reasonable estimates of the diffusion coefficients of key
chemical species, as small solutes can undergo anomalous diffusion in water
Kirchner, Stubbs, and Marx (2002); Roberts  _et al._ (2009); Marx, Chandra,
and Tuckerman (2010), particularly under extreme conditions.
$\cdot\mathrm{OH}$, the focal species of this work, has been the subject of
previous studies; however, these have generally focused on high temperature,
e.g. (Svishchev and Plugatyr, 2005). The different patterns of
$\cdot\mathrm{OH}$ diffusion seen in low-temperature, high-pressure regime
suggest considerable value in further experimental studies of these
environments.

Figure 3: Posterior mean estimates and 95% posterior intervals for
$D_{\mathrm{H}_{2}\mathrm{O}}$ and $D_{\cdot\mathrm{OH}}$. (In some cases,
interval width is narrower than point sizes.) Both behave similarly at low and
very pressures, but $D_{\cdot\mathrm{OH}}$ lacks the accelerated diffusion
seen in $\mathrm{H}_{2}\mathrm{O}$ at moderately high pressure.

To better summarize the posterior surface of $D_{\cdot\mathrm{OH}}$ in
response to temperature and pressure, we fit a least-squares log-log
polynomial approximation to the posterior means. The resulting summary model
($R^{2}=0.995$, standard error 0.031 on 13 residual degrees of freedom (log
scale), RMSE $0.0035\unit{\square\per\pico}$ on the phenomenological scale) is
given by

$\begin{split}\widehat{D_{\cdot\mathrm{OH}}}=\exp\left[$-4.410\,203\text{\times}{10}^{02}$+$1.476\,812\text{\times}{10}^{02}$T^{\prime}\right.\\\
\left.+$8.225\,372\text{\times}{10}^{-01}$P^{\prime}+$-1.238\,130\text{\times}{10}^{01}$T^{\prime
2}\right.\\\ \left.+$2.814\,801\text{\times}{10}^{-02}$P^{\prime
2}+$-1.582\,325\text{\times}{10}^{-01}$T^{\prime}P^{\prime}\right.\\\
\left.+$-2.622\,113\text{\times}{10}^{-03}$P^{\prime
3}\right]\unit{\square\per\pico}\end{split}$ (12)

where $T^{\prime}$ and $P^{\prime}$ are the dimensionless quantities
$T^{\prime}=\log\left(T\;\unit{\per}\right)$ and
$P^{\prime}=\log\left(P\;\unit{\per}\right)$. The resulting surface is shown
in Fig. 4. As expected from the known coincidence of $D_{\cdot\mathrm{OH}}$
and $D_{\mathrm{H}_{2}\mathrm{O}}$ under ambient conditions,
$D_{\cdot\mathrm{OH}}$ has broadly similar behavior to
$D_{\mathrm{H}_{2}\mathrm{O}}$. However, $\cdot\mathrm{OH}$ shows less
accelerated diffusion at moderate pressures, and slows considerably at very
high pressures (even in excess of $\mathrm{H}_{2}\mathrm{O}$).

Figure 4: $\widehat{D_{\cdot\mathrm{OH}}}$ as a function of temperature and
pressure, based on posterior inference from simulations in TIP4P/2005 water.
$\cdot\mathrm{OH}$ diffusion is slightly faster at moderate pressures when
temperatures are low (though less so than $\mathrm{H}_{2}\mathrm{O}$), but
slows considerably at high pressures.

##  III Conclusion

Diffusion rates for the hydroxyl radical differ substantially at low
temperatures and high pressures. For the most part, the pattern of variation
in $D_{\cdot\mathrm{OH}}$ qualitatively follows that of
$D_{\mathrm{H}_{2}\mathrm{O}}$ (as would be expected for a small molecule in
aqueous solution), but we do not observe the same extent of enhanced diffusion
at moderate pressures (particularly at higher temperatures), and diffusion
slows more substantially as pressures approach 10000 atm. This may have an
impact on the efficiency of radical-scavenging mechanisms in high-pressure
terrestrial environments, and hypothetically in the even higher pressure
conditions predicted to obtain within some outer solar system oceans.

To estimate the $\cdot\mathrm{OH}$ diffusion coefficient, we combined
atomistic molecular dynamics with hierarchical Bayesian inference that allows
us to easily pool information between solvent and radical, and across multiple
trajectories of varying size. The modular structure of Bayesian models lends
itself naturally to this application, as does the ease with which one can e.g.
leverage _a priori_ information on model parameters or account for variance
across trajectories. This strategy (which builds on other recent work on
statistical inference for the diffusion coefficient (Bullerjahn, von Bülow,
and Hummel, 2020)) is readily adapted to the estimation of diffusion
coefficients for other systems, or to other, related physical quantities.

###### Acknowledgements.

This work was supported by NASA award 80NSSC20K0620 to R.W.M. and C.T.B. We
thank George Miller for helpful discussions about the chemistry of water
radiolysis.

## Data Availability Statement

Data sharing not applicable – no new data generated. CHARMM-compatible
parameters for the $\cdot\mathrm{OH}$ model are included in the main text.

## References

  * Blanco  _et al._ (2018) Y. Blanco, G. de Diego-Castilla, D. V. dez Moreiras, E. Cavalcante-Silva, J. A. Rodriguez-Manfredi, A. F. Davila, C. P. McKay, and V. Parro, “Effects of gamma and electron radiation on the structural integrity of organic molecules and macromolecular biomarkers measured by microarray immunoassays and their astrobiological implications,” Astrobiology 18, DOI: 10.1089/ast.2016.1645 (2018).
  * Onstott  _et al._ (2003) T. C. Onstott, D. P. Moser, S. M. Pfiffner, J. K. Fredrickson, F. J. Brockman, T. J. Phelps, D. C. White, A. Peacock, D. Balkwill, R. Hoover, L. R. Krumholz, M. Borscik, T. L. Kieft, and R. Wilson, “Indigenous and contaminant microbes in ultradeep mines,” Environmental Microbiology 5, 1168–1191 (2003).
  * Tarnas  _et al._ (2018) J. Tarnas, J. Mustard, B. S. Lollar, M. Bramble, K. Cannon, A. Palumbo, and A.-C. Plesa, “Radiolytic H2 production on Noachian Mars: Implications for habitability and atmospheric warming,” Earth and Planetary Science Letters 502, 133–145 (2018).
  * Ojha  _et al._ (2021) L. Ojha, S. Karunatillake, S. Karimi, and J. Buffo, “Amagmatic hydrothermal systems on mars from radiogenic heat,” Nature Communications 12, https://doi.org/10.1038/s41467–021–21762–8 (2021).
  * Tarnas  _et al._ (2021) J. Tarnas, J. Mustard, B. S. Lollar, V. Stamenković, K. Cannon, J.-P. Lorand, T. Onstott, J. Michalski, O. Warr, A. Palumbo, and A.-C. Plesa, “Earth-like habitable environments in the subsurface of mars,” Astrobiology 21, https://doi.org/10.1089/ast.2020.2386 (2021).
  * White  _et al._ (1999) O. White, J. A. Eisen, J. F. Heidelberg, E. K. Hickey, J. D. Peterson, R. J. Dodson, D. H. Haft, M. L. Gwinn, W. C. Nelson, D. L. Richardson, K. S. Moffat, H. Qin, L. Jiang, W. Pamphile, M. Crosby, M. Shen, J. J. Vamathevan, P. Lam, L. McDonald, T. Utterback, C. Zalewski, K. S. Makarova, L. Aravind, M. J. Daly, K. W. Minton, R. D. Fleischmann, K. A. Ketchum, K. E. Nelson, S. Salzberg, H. O. Smith, V. J. Craig, and C. M. Fraser, “Genome sequence of the radioresistant bacterium Deinococcus radiodurans R1,” Science 286, 1571–1577 (1999).
  * Cavicchioli (2002) R. Cavicchioli, “Extremophiles and the search for extraterrestrial life,” Astrobiology 2, 281–292 (2002).
  * Munteanu, Uivarosi, and Andries (2015) A.-C. Munteanu, V. Uivarosi, and A. Andries, “Recent progress in understanding the molecular mechanisms of radioresistance in Deinococcus bacteria,” Extremophiles 9, 707–719 (2015).
  * LaVerne (2000) J. A. LaVerne, “OH radicals and oxidizing products in the gamma radiolysis of water,” Radiation Research 153, 196–200 (2000).
  * Matheson (1964) M. S. Matheson, “The formation and detection of intermediates in water radiolysis,” Radiation Research Supplement 4, 1–23 (1964).
  * Ghosal  _et al._ (2005) D. Ghosal, M. Omelchenko, E. Gaidamakova, M. V.Y., A. Vasilenko, A. Venkateswaran, M. Zhai, H. Kostandarithes, H. Brim, K. Makarova, L. Wackett, J. Fredrickson, and M. Daly, “How radiation kills cells: Survival of Deinococcus radiodurans and Shewanella oneidensis under oxidative stress,” FEMS Microbiology Reviews 29, 361–375 (2005).
  * Roots and Okada (1975) R. Roots and S. Okada, “Estimation of life times and diffusion distances of radicals involved in X-ray-induced DNA strand breaks or killing of mammalian cells,” Radiation Research 64, 306–320 (1975).
  * Du and Gebicki (2004) J. Du and J. M. Gebicki, “Proteins are major initial cell targets of hydroxyl free radicals,” The International Journal of Biochemistry and Cell Biology 36, 2334–2343 (2004).
  * Davies, Fu, and Dean (1995) M. Davies, S. Fu, and R. T. Dean, “Protein hydroperoxides can give rise to reactive free radicals,” Biochemical Journal 305, 643–649 (1995).
  * Barnham, Masters, and Bush (2004) K. J. Barnham, C. L. Masters, and A. I. Bush, “Neurodegenerative diseases and oxidative stress,” Nature Reviews Drug Discovery 3, 205–214 (2004).
  * Best  _et al._ (2012) R. B. Best, X. Zhu, J. Shim, P. E. M. Lopes, J. Mittal, M. Feig, and J. A. D. Mackerell, “Optimization of the additive CHARMM all-atom protein force field targeting improved sampling of the backbone $\phi$, $\psi$ and side-chain $\chi$(1) and $\chi$(2) dihedral angles,” Journal of Chemical Theory and Computation 8, 3257–3273 (2012).
  * Absacal and Vega (2005) J. L. F. Absacal and C. Vega, “A general purpose model for the condensed phases of water: TIP4P/2005,” Journal of Chemical Physics 123, 234505 (2005).
  * Tsimpanogiannis  _et al._ (2019) I. N. Tsimpanogiannis, O. A. Moultos, L. F. Franco, M. B. M. Spera, M. Erdös, and I. G. Economou, “Self-diffusion coefficient of bulk and confined water: a critical review of classical molecular simulation studies,” Molecular Simulation 45, 425–453 (2019).
  * Pabis, Szala-Bilnik, and Swiatla-Wojcik (2011) A. Pabis, J. Szala-Bilnik, and D. Swiatla-Wojcik, “Molecular dynamics study of the hydration of the hydroxyl radical at body temperature,” Physical Chemistry Chemical Physics 13, 9458–9468 (2011).
  * Bopp, Jancsó, and Heinzinger (1983) P. Bopp, G. Jancsó, and K. Heinzinger, “An improved potential for non-rigid water molecules in the liquid phase,” Chemical Physics Letters 98, 129–133 (1983).
  * Dorfman and Adams (1973) L. M. Dorfman and G. E. Adams,  _Reactivity of the Hydroxyl Radical in Aqueous Solutions_ (U.S. Department of Commerce, National Bureau of Standards, Washington, D.C., 1973).
  * Martyna, Tobias, and Klein (1994) G. J. Martyna, D. J. Tobias, and M. L. Klein, “Constant pressure molecular dynamics algorithms,” Journal of Chemical Physics 101, 4177–4189 (1994).
  * Feller  _et al._ (1995) S. E. Feller, Y. Zhang, R. W. Pastor, and B. R. Brooks, “Constant pressure molecular dynamics simulation: The Langevin piston method,” Journal of Chemical Physics 103, 4613–4621 (1995).
  * Phillips  _et al._ (2005) J. C. Phillips, R. Braun, W. Wang, J. Gumbart, E. Tajkhorshid, E. Villa, C. Chipot, R. D. Skeel, L. Kalé, and K. Schulten, “Scalable molecular dynamics with NAMD,” Journal of Computational Chemistry 26, 1781–1802 (2005).
  * Humphrey, Dalke, and Schulten (1996) W. Humphrey, A. Dalke, and K. Schulten, “VMD: Visual molecular dynamics,” Journal of Molecular Graphics 14, 33–38, 27–28 (1996).
  * Ribeiro  _et al._ (2020) J. a. V. Ribeiro, B. Radak, J. Stone, J. Gullingsrud, J. Saam, and J. Phillips, “psfgen plugin for VMD,” Software File (2020).
  * Martínez  _et al._ (2009) L. Martínez, R. Andrade, E. Birgin, and J. M. Martínez, “Packmol: A package for building initial configurations for molecular dynamics simulations,” Journal of Computational Chemistry 30, 2157–2164 (2009).
  * von Bülow, Bullerjahn, and Hummel (2020) S. von Bülow, J. T. Bullerjahn, and G. Hummel, “Systematic errors in diffusion coefficients from long-time molecular dynamics dimulations at constant pressure,” Journal of Chemical Physics 153, 021101 (2020).
  * Bullerjahn, von Bülow, and Hummel (2020) J. T. Bullerjahn, S. von Bülow, and G. Hummel, “Optimal estimates of self-diffusion coefficients from molecular dynamics simulations,” Journal of Chemical Physics 153, 024116 (2020).
  * Yeh and Hummer (2004) I.-C. Yeh and G. Hummer, “System-size dependence of diffusion coefficients and viscosities from molecular dynamics simulations with periodic boundary conditions,” Journal of Physical Chemistry, B 108 (2004).
  * González and Abascal (2010) M. A. González and J. L. F. Abascal, “The shear viscosity of rigid water models,” Journal of Chemical Physics 132, 096101 (2010).
  * Gelman  _et al._ (2003) A. Gelman, J. B. Carlin, H. S. Stern, and D. B. Rubin,  _Bayesian Data Analysis_ , 2nd ed. (Chapman and Hall, London, 2003).
  * R Core Team (2021) R Core Team,  _R: A Language and Environment for Statistical Computing_, R Foundation for Statistical Computing, Vienna, Austria (2021).
  * Nash (1990) J. C. Nash,  _Compact Numerical Methods for Computers. Linear Algebra and Function Minimization_ (Adam Hilger, 1990).
  * Scrucca  _et al._ (2016) L. Scrucca, M. Fop, T. B. Murphy, and A. E. Raftery, “mclust 5: clustering, classification and density estimation using Gaussian finite mixture models,” The R Journal 8, 289–317 (2016).
  * Homan and Gelman (2014) M. D. Homan and A. Gelman, “The No-U-Turn sampler: Adaptively setting path lengths in Hamiltonian Monte Carlo,” Journal of Machine Learning Research 15, 1593–1623 (2014).
  * Stan Development Team (2020a) Stan Development Team, “RStan: the R interface to Stan,”  (2020a), r package version 2.21.2.
  * Stan Development Team (2020b) Stan Development Team, “Stan: A C++ library for probability and sampling,”  (2020b), software Library.
  * Gelman and Rubin (1992) A. Gelman and D. B. Rubin, “Inference from iterative simulation using multiple sequences (with discussion),” Statistical Science 7, 457–511 (1992).
  * Schwarz (1969) H. A. Schwarz, “Applications of the spur diffusion model to the radiation chemistry of aqueous solutions,” Journal of Physical Chemistry 73, 1928–1937 (1969).
  * Parajon  _et al._ (2008) M. H. Parajon, P. Rajesh, T. Mu, S. M. Pimblott, and J. A. LaVerne, “H atom yields in the radiolysis of water,” Radiation Physics and Chemistry 77, 1203–1207 (2008).
  * Ershov and Gordeev (2008) B. Ershov and A. Gordeev, “A model for radiolysis of water and aqueous solutions of H2 H2O2 and O2,” Radiation Physics and Chemistry 77, 928–935 (2008).
  * Omar, Hasnaoui, and de la Lande (2021) K. A. Omar, K. Hasnaoui, and A. de la Lande, “First-principles simulations of biological molecules subjected to ionizing radiation,” Annual Review of Physical Chemistry 72, 445–465 (2021).
  * Loh  _et al._ (2020) Z.-H. Loh, G. Doumy, C. Arnold, L. Kjellsson, S. H. Southworth, A. A. Haddad, Y. Kumagai, M.-F. Tu, P. J. Ho2, A. M. March, R. D. Schaller, M. S. B. M. Yusof, T. Debnath, M. Simon, R. Welsch, L. Inhester, K. Khalili, K. Nanda, A. I. Krylov, S. Moeller, G. Coslovich, J. Koralek, M. P. Minitti, W. F. Schlotter, J.-E. Rubensson, R. Santra, and L. Young, “Observation of the fastest chemical processes in the radiolysis of water,” Science 367, 179–182 (2020).
  * Jordan  _et al._ (2020) I. Jordan, M. Huppert, D. Rattenbacher, M. Peper, D. Jelovina, C. Perry, A. von Conta, A. Schild, and H. J. Wörner, “Attosecond spectroscopy of liquid water,” Science 369, 974–979 (2020).
  * Thürmer  _et al._ (2013) S. Thürmer, M. Ončák, N. Ottosson, R. Seidel, U. Hergenhahn, S. E. Bradforth, P. Slavíček, and B. Winter, “On the nature and origin of dicationic, charge-separated species formed in liquid water on x-ray irradiation,” Nature Chemistry 5, 590–596 (2013).
  * Ren  _et al._ (2018) X. Ren, E. Wang, A. D. Skitnevskaya, A. B. Trofimov, K. Gokhberg, and A. Dorn, “Experimental evidence for ultrafast intermolecular relaxation processes in hydrated biomolecules,” Nature Physics 14, 1062–1066 (2018).
  * Kirchner, Stubbs, and Marx (2002) B. Kirchner, J. Stubbs, and D. Marx, “Fast anomalous diffusion of small hydrophobic species in water,” Physical Review Letters 89, 215901 (2002).
  * Roberts  _et al._ (2009) S. T. Roberts, P. B. Petersen, K. Ramasesha, A. Tokmakoff, I. S. Ufimtsev, and T. J. Martinez, “Observation of a Zundel-like transition state during proton transfer in aqueous hydroxide solutions,” Proceedings of the National Academy of Sciences of the United States of America 106, 15154–15159 (2009).
  * Marx, Chandra, and Tuckerman (2010) D. Marx, A. Chandra, and M. E. Tuckerman, “Aqueous basic solutions: Hydroxide solvation, structural diffusion, and comparison to the hydrated proton,” Chemical Reviews 110, 2174–2216 (2010).
  * Svishchev and Plugatyr (2005) I. M. Svishchev and A. Y. Plugatyr, “Hydroxyl radical in aqueous solution: Computer simulation,” J. Phys. Chem. B 109, 4123–4128 (2005).

Supplemental Materials for _Bayesian Estimation of the Hydroxyl Radical
Diffusion Coefficient at Low Temperature and High Pressure from Atomistic
Molecular Dynamics_
Carter T. Butts and Rachel W. Martin
University of California, Irvine
7/28/21

Figure S1: AREs for estimated $D_{\cdot\mathrm{OH}}$, 20 Åbox size
simulations, by parameter value. Parameters are indicated as follows:
epsO=$\epsilon_{\mathrm{O}}$, hRminO=$r_{\min\mathrm{O}}{/\scriptstyle 2}$,
epsH=$\epsilon_{\mathrm{H}}$, hRminH=$r_{\min\mathrm{H}}{/\scriptstyle 2}$.
$r_{\min\mathrm{O}}{/\scriptstyle 2}$ is a critical parameter, with values
between apx 1.3-1.8 Åleading to especially poor diffusion performance. Two
broad regions with high-performing parameters are visible, motivating
selection by examination of other aspects of model behavior.  Figure S2:
Lennard-Jones energy for interaction of water oxygen with $\cdot\mathrm{OH}$
atoms (as TIP4P/2005 includes no LJ effects for water protons, these are not
included, and there are no orientation effects). Diffusion-accurate parameter
values are characterized by an exclusion of attractive O-O interactions within
a particular distance band, but have no specific relationship to the O-H
distance.  Figure S3: Winnowing of parameter vectors by hydrogen bonding
behavior with $\mathrm{H}_{2}\mathrm{O}$ as acceptor; curves indicate energies
for specific parameter vectors, with the vector selected in the initial round
of optimization indicated in black. (A) total $\mathrm{H}_{2}\mathrm{O}$,
$\cdot\mathrm{OH}$ interaction energy in an O-O facing configuration by
distance. (B) equivalent energies for the O-H configuration
($\mathrm{H}_{2}\mathrm{O}$ donor). (C) difference between O-O and O-H
energies by trajectory (positive values indicate that the O-O orientation is
correctly unfavorable at the specified distance); note that, while some low-
ARE parameters select the incorrect orientation, the chosen parameter vector
favors the correct state. (D) AREs by parameter (epsO=$\epsilon_{\mathrm{O}}$,
hRminO=$r_{\min\mathrm{O}}{/\scriptstyle 2}$, epsH=$\epsilon_{\mathrm{H}}$,
hRminH=$r_{\min\mathrm{H}}{/\scriptstyle 2}$), with points scaled by extent to
which the corresponding parameter factors the correct versus incorrect
orientation. Restricting consideration to parameters with correct hydrogen
bonding behavior greatly narrows the range of viable parameters.  Figure S4:
Winnowing of parameter vectors by hydrogen bonding behavior with
$\mathrm{H}_{2}\mathrm{O}$ as donor; curves indicate energies for specific
parameter vectors, with the vector selected in the initial round of
optimization indicated in black. (A) total $\mathrm{H}_{2}\mathrm{O}$,
$\cdot\mathrm{OH}$ interaction energy in an H-O facing configuration by
distance ($\mathrm{H}_{2}\mathrm{O}$ donor). (B) equivalent energies for the
H-H configuration). (C) difference between H-H and O-H energies by trajectory
(positive values indicate that the H-H orientation is correctly unfavorable at
the specified distance); the chosen parameter vector favors the correct state
except at unphysically short radii (at which neither configuration is
favorable). (D) AREs by parameter (epsO=$\epsilon_{\mathrm{O}}$,
hRminO=$r_{\min\mathrm{O}}{/\scriptstyle 2}$, epsH=$\epsilon_{\mathrm{H}}$,
hRminH=$r_{\min\mathrm{H}}{/\scriptstyle 2}$), with points scaled by extent to
which the corresponding parameter factors the correct versus incorrect
orientation. Restricting consideration to parameters with correct hydrogen
bonding behavior narrows the range of viable parameters, but less so than the
$\mathrm{H}_{2}\mathrm{O}$ acceptor case.


# Spade: An Expression-Based HDL With Pipelines

Frans Skarman and Oscar Gustafsson  Department of Electrical Engineering,
Linköping University
SE-581 83 Linköping, Sweden
Email: {frans.skarman<EMAIL_ADDRESS>
###### Abstract

Spade is a new open source hardware description language (HDL) designed to
increase developer productivity without sacrificing the low-level control
offered by HDLs. It is a standalone language which takes inspiration from
modern software languages, and adds useful abstractions for common hardware
constructs. It also comes with a convenient set of tooling, such as a helpful
compiler, a build system with dependency management, tools for debugging, and
editor integration.

###### Index Terms:

Hardware description languages, languages and compilers, Design automation

##  I Introduction

Developing digital hardware is traditionally done using Verilog or VHDL, both
languages originating in the 1980s. While they have evolved since their
inception, they are still lacking many subsequent advancements in language
design.

Spade111https://gitlab.com/spade-lang/spade/ is a new open source Hardware
Description Language (HDL) which aims to reduce the development effort of
digital hardware by taking inspiration from software languages, and adding
language-level mechanisms for common hardware structures. This is while still
retaining the low-level control provided by HDLs.

Being inspired by Rust and functional programming languages, Spade is
expression-based and has a rich type system. It supports product-types like
structs and tuples, and sum-types in the form of enums. The language also
supports linear type-checking which can be used to ensure that hardware
resources such as memory ports are used exactly once. †† This work was
presented at the _3rd Workshop on Open-Source Design Automation (OSDA) 2023_ ,
co-hosted with _Design, Automation, and Test in Europe Conference_ _(DATE)_
_2023_ in Antwerp, Belgium, on April 17, 2023. For further information and
material, please visit the paper’s web page on the workshop website at
https://osda.ws/r/2QkTi (just scan the QR code on the right). The copyright is
retained by the authors. This work is licensed under the terms of the Creative
Commons Attribution license (CC BY), version 4.0, see
https://creativecommons.org/licenses/by/4.0. _Proceedings of the 3rd Workshop
on Open-Source Design Automation (OSDA)_ , Pages 7–12, 2023.

Spade has built in constructs and abstractions for common hardware structures
such as pipelines, memories, and registers. Pipelining allows the user to
specify a computation to be performed, with explicit statements for separating
stages of the pipeline, but without the need of separate variables for each
pipeline stage. Retiming such a pipeline does not require changing any
variables, only moving the staging statement. Additionally, the delays of
pipelines are explicit in the language and checked by the compiler to ensure
that changes to a pipeline do not affect the computation results. In order to
more accurately reflect the hardware being described, all logic is
combinatorial by default with the only sequential elements being registers and
memories which are instantiated explicitly.

The rest of the paper is structured as follows. Related works are discussed in
Section II, including highlighting where Spade differs. The basic semantics
are introduced in Section III, while in Section IV the use of linear types to
model input and output ports is described. The software provided for Spade is
discussed in Section V, with concluding remarks in Section VI.

##  II Related Work

In recent years, several new HDLs have been developed. A common approach is to
embed the language as a Domain Specific Language (DSL) inside a conventional
software programming language such as Scala, Python, or Ruby. There, the HDL
consists of a library of hardware constructs which the user instantiates in
the host language in order to describe their hardware. This allows the user to
take advantage of the power of the host language in their hardware
description, for example, by using software control flow structures to
generate parameterized hardware and using object-oriented or functional
programming approaches in the hardware description. Additionally, this
approach exposes all the tooling available for the host language to the
hardware designer, such as dependency managers, build tools, and IDEs.
Examples of this include Chisel [1], SpinalHDL [2] and DFiant [3] embedded in
Scala, Amaranth [4] embedded in Python, ROHD [5] embedded in Dart, and RubyRTL
[6] embedded in Ruby.

Some languages take subsets of conventional programming languages and compile
them to hardware. An example of this is Clash [7] which compiles a large
subset of Haskell to hardware by taking advantage of the natural mapping
between a pure functional language and hardware. Another example is PipelineC
[8] which is a C-like HDL with automatic pipelining. While not regular C, it
is close enough to it to allow many PipelineC programs to be compiled and
“simulated” by a standard C-compiler.

TL-Verilog [9] is a SystemVerilog extension which supports timing abstract
modeling, where behavior and timing are separated. It provides language
features for common hardware constructs such as pipelines and FIFOs.

There are also several languages which are completely independent of existing
languages. One such example is Bluespec [10] in which hardware is described by
guarded atomic actions. Another example is Silice [11] which contains
abstractions for common hardware constructs without losing control over the
generated hardware. In addition, it provides a higher-level description style
where the design is expressed as sequences of operations with software-like
control flow.

Finally, a common alternative for describing hardware is High-Level Synthesis
(HLS), in which higher-level languages, typically designed for software, are
compiled to hardware. This design methodology is quite different from the HDLs
described previously. Specifically, HLS tools generally provide limited
control over the hardware being generated. In particular, the exact timing of
any circuit is generally abstracted away, and synchronization between HLS
generated modules is done at runtime via synchronization signals.

Now, a natural question to ask is whether there is a need for another HDL, and
what Spade offers that the existing work does not? First, it is not an
embedded DSL in another language, which separates it from the likes of Chisel
and Amaranth. This gives more freedom in the design of the language, as it is
not restricted to the expressiveness of a DSL. For example, most embedded HDLs
are forced to invent new “keywords” in order to not clash with the host
language, a typical example being Chisel using when, and otherwise instead of
the more familiar if, and else. This problem runs deeper than surface level
keywords, however. As an example, Scala has pattern matching, but in Chisel
and SpinalHDL, that feature can only be used at “build time” on Scala values,
not in hardware on “runtime” values. Finally, a potential user is not required
to learn both the host language and DSL when initially picking up the
language.

Unlike Clash and PipelineC, Spade is not restricted to conforming to the
execution model of the host language which means that the language can be
designed directly for hardware description without any restrictions.

Spade describes the behavior of circuits in a cycle-to-cycle manner, which
gives more control than Bluespec, where the fundamental abstraction is atomic
rules, and PipelineC where the compiler automatically performs pipelining.

Finally, Spade is similar in spirit to TL-Verilog. However, by designing a new
language instead of building on SystemVerilog, there is more freedom to
explore new language constructs.

##  III Basic Semantics

In order to describe the programming model and semantics of Spade, an example
is used. Listing 1 shows Spade code which blinks an LED at a configurable
interval. Lines 1–2 define an entity called blink which takes a clock signal
of type clock, a reset signal of type bool, and a 20-bit integer, int<20>,
specifying the maximum value for the counter. The return value is a single
bool value which drives the led. The entity keyword specifies that this unit
can contain sequential logic. One can also write a function, fn, which only
allows combinatorial logic, and pipeline which defines a pipeline. The details
of pipelines are discussed later. When instantiating a unit, the syntax
differs between functions, entities and pipelines, giving a reader of the code
an indication of what can happen behind an instantiation.

Listing 1: Spade code which blinks an LED

⬇

1entity blink(clk: clock, rst: bool, max: int<20>)

2 -> bool

3{

4 reg(clk) counter reset (rst: 0) =

5 if counter == max {

6 0

7 } else {

8 trunc(counter + 1)

9 };

10 counter > (max >> 1)

11}

Line 4 defines a register called counter which is clocked by the clk signal
and reset to 0 by the rst signal. Registers in Spade are explicit constructs
rather than being inferred from the use of something like rising_edge.
Registers are, together with writes to memories, the only sequential
constructs in the language, everything else describes combinatorial circuits.

In Spade, the behavior of a circuit is described in a cycle-to-cycle manner.
The new value of a register in the design is given as an expression of the
register values in the current clock cycle. All variables in Spade are
immutable: they can only be assigned once. Single assignment is possible
because Spade is an expression-oriented language, meaning that most statements
are expressions and produce values[12]. For example, the counter register is
set to the value of the expression following the equals sign on line 4, in
this case an if-expression spanning lines 5–10. If the counter has reached the
max value, the value returned by the if expression is 0, otherwise it is the
next value of the counter. Compared to traditional imperative HDLs, the use of
immutable variables and expression-based control flow is closer to the
resulting hardware where variables correspond to wires and “control flow” is
implemented with multiplexers. The language also does not have an explicit
return keyword, the last expression in a unit body is the output of the
entity, here on line 11.

The call to the trunc function on line 9 truncates the value of (counter + 1)
from 21 bits down to 20 bits. This is needed as Spade prevents accidental
overflow by using the largest possible values for most arithmetic operations.
More details on the type system is given in Section III-B

Spade has similar scoping rules to most modern software languages: variables
are only visible below their definition which makes it more difficult to
accidentally create combinatorial loops. Additionally, these scoping rules
make it easier for a developer to find the definition of a variable. The
definition will always be above where it is used, and generally be grouped
with its assignment. Sometimes it is still useful to create a loop of
dependencies between variables, however, this is done explicitly using the
decl keyword.

###  III-A Pipelines

Pipelining is an important construct in most hardware designs, it allows
designs to maintain a high clock frequency and throughput at the cost of
latency. However, despite their importance, most HDLs require the user to
manually build their pipelines, a process that is tedious and error-prone as
one must make sure that computations are performed on values corresponding to
the correct time step. In some cases, designers use patterns on their variable
names, and ad-hoc static checking tools to verify that the pipelining is
correct [13].

⬇

, autogobble=true

, caption={Example of the pipelining construct in Spade.}

, label=lst:pipeline

, numbers=left

, numberstyle=\color{gray}

, xleftmargin=2em

]

pipeline(4) X(clk: clock, a: int<32>, b: int<32>)

-> int<33>

{

’initial

let x = inst(3) subpipe(clk, a);

let p = a * b;

reg * 3;

let s = x + f(a, p);

reg;

s + stage(initial).a

}

pipeline(3) subpipe(...) -> int<32> {...}

Figure 1: Hardware described in pipeline X by the code in Listing
LABEL:lst:pipeline.

Spade on the other hand has language-level support for describing pipelines
where the user describes which computations are performed in each stage, and
the compiler manages the insertion of registers between stages. Listing
LABEL:lst:pipeline shows an example of the pipelining construct in use to
describe the hardware in Fig. 1. Like the entity discussed previously, the
first two lines describe the external interface of the pipeline: name, inputs,
and outputs. It also specifies depth, in this case 4, which is the number of
registers in the pipeline and therefore the delay between input and output.
While the compiler can infer the depth, it has to be specified manually
because it is part of the public API. This allows a user to see the interface
of the pipeline without reading the body.

Line 4 names the current pipeline stage initial, allowing references to it
later in the code. On line 5, another pipeline called subpipe is instantiated
using the inst keyword. When instantiating pipelines, the depth of the
pipeline has to be specified, and the compiler checks this against the depth
specified in the head of the instantiated pipeline. If the depth of the
instantiated pipeline is changed later, the compiler gives an error, forcing
the developer to consider whether that change has an effect in the
instantiating code. The compiler also checks the availability of variables
from sub-pipelines, giving an error if the result is used before it is ready.
As an example, if x is used in the stage right after its definition, it will
emit the error message shown in Listing 2.

Line 6 defines a new variable called p, which is the product of a and b, and,
from type inference, an int<64>. Following that is a reg-statement, which
behaves differently from the reg statement in an entity that was described
previously. In a pipeline, the reg-statement tells the compiler to insert a
pipeline stage, registering all the values visible above it. After the reg-
statement on line 7, any references to p will refer to a delayed version of
the value. The same is true for the x value computed by subpipe, however,
because it is computed by a pipeline of depth three, the first three registers
are present in the sub-pipe and will not be inserted in the instantiating
pipeline. The times three operator on the register, reg * 3, specifies that
three stages should be inserted without any new computations. Finally, on line
10, s is added to the value of a from the first stage. The stage keyword is
used to bypass some pipeline registers and to refer directly to a signal as it
appears in another stage, in this case the stage named initial. Stage
references can refer to both “future” and “past” stages, and can do so by
name, as is done here, or as relative offsets, for example stage(-2).

Listing 2: Example of pipeline error message

⬇

errorr@: Use of x before it is ready

\--> src/main.spade:10:19

|

10 | let sum = x + f(a, product);

| ^ Is unavailable for another 2 stagesr@

= Requesting x from stage 1

= But it will not be available until stage 3

The pipelining feature decouples the description of the computation from the
description of the pipeline itself. In a pipeline without feedback a developer
can easily add and remove pipeline stages as needed without altering the
output value of the pipeline. If such changes potentially affect the outcome
of other parts of the project, the inclusion of the depth at the call site
ensures that the user is made aware of such potential issues via compiler
error messages. In pipelines with feedback the structured description of the
stages still helps during design iteration, even if some manual care is
required to ensure correct computation.

###  III-B Types and Pattern Matching

Spade is a statically typed language with type inference. This means that
types can be omitted in most code since the compiler will infer the
appropriate types from context, and report errors if types can not be
inferred.

Like most languages, Spade has primitive types such as integers and booleans,
and compound types like arrays, tuples, and structs. In addition, the language
supports enum types inspired by software languages like Rust, Haskell, and ML.
Unlike their C or VHDL namesake, these enums have data associated with them in
addition to being one of a set of variants. A common use case for this
construct is the Option type which is defined in the standard library as shown
in Listing 3. It is generic over a contained type T and takes on one of two
values: Some in which case a value, val, of type T is present, or None in
which case no such value is present. One can view the Option type as having a
valid/invalid signal bundled with the data it validates.

The main way to interact with enums in Spade is the match-expression, which
allows pattern matching on values. As an example, Listing 4 shows the match-
expression in use on a tuple of Option-values: a and b. The resulting hardware
is shown in Fig 2. If a is Some, its inner value is bound to val and returned
from the match expression. If a is None but b is Some, its inner value is
returned. Finally, if both a and b are None, 0 is returned. Enums are encoded
as bits that specify the currently active variant, the discriminant, followed
by bits containing the payload. In this case, since the enum consists of two
variants, the discriminant is a single bit while the remaining bits contain
the val when present and are undefined otherwise.

Listing 3: Definition of the Option type.

⬇

enum Option<T> {

None,

Some{ val: T },

}

Listing 4: Example of pattern matching on a tuple of Option values.

⬇

let a: Option<T> = ...;

let b: Option<T> = ...;

let result = match (a, b) {

(Some(val), _) => val,

(_, Some(val)) => val,

_ => 0

}

Figure 2: Hardware generated by Listing 4.

###  III-C Memories

As discussed earlier, Spade requires register definitions to include an
expression for the new value of the whole register as a function of the values
from the previous clock cycle. However, this abstraction becomes problematic
when working with memory-like structures in which only a small part of the
total state is updated in each cycle. In order to mitigate this, Spade has an
explicit construct for memories, currently implemented as entities defined in
the standard library which the compiler handles separately. Instantiation of a
memory is done using the clocked_memory entity which creates a memory with a
fixed set of write ports. Reading from said memory is done via the similarly
defined read_memory entity. This means that unlike VHDL and Verilog, memories
are explicitly instantiated as memories, rather than inferred from the
structure of the code. An example of using memories is included in the next
Section.

##  IV Ports

The discussion so far has been centered around computations on values. A unit
receives a set of values as inputs, and produces a set of values as output.
Internally, values are used in computation and are pipelined by pipelines.
This becomes inconvenient when working with something like a memory which is
external to the current unit. The unit must return an address, a write-enable
signal, and a value to write to the memory. Then the value read from the
memory must be passed as an input to the unit. The memory in turn, produces an
output value which must be fed back into the controlling unit. This approach
has a few issues: the control signals and output are delayed by the pipeline
mechanism. This pipelining introduces additional delays unless manually
mitigated, for example by stage references. Additionally, there is no clear
link between memory control-signals and the corresponding output.

In order to mitigate this issue, the Spade type system contains the concept of
ports and wires. Wires come in two forms: mutable and immutable denoted by
&mut and & respectively. An immutable wire can be used to pass values via
units without them being delayed in pipelines. Mutable wires are similarly not
pipelined but allow setting the value of the wire in a module which takes the
wire as an input. Finally, ports are collections of wires and other ports
which are bundled together.

Listing 5 contains an example of how the port feature can be used to share a
memory between two pipelines. The first lines define port-types containing
read, write, and address wires. Line 8 defines an entity where the actual
memory is instantiated. It returns a read-port, and a write-port. Line 9 has
been trimmed for space but instantiates the mut-wire r_addr and a WPort. Lines
10–11 instantiate a memory with a single write-port using the clocked_memory
entity, and lines 12–14 asynchronously reads one value from the memory. At the
end of the entity, the read-port is assembled from the mutable address, and
the result wire, and returned along with the write-port.

The pipelines reader and writer are two pipelines which use the read-port and
write-port respectively. On line 20, the reader pipeline sets the address it
wants to read from, and reads the resulting value on line 21. Similarly, the
writer sets the target address and write-value on line 27.

Finally, the three units are instantiated on lines 32–34. The ports returned
by the memory are passed to the reader and writer.

Listing 5: Ports being used to share access to a memory between modules.

⬇

1struct port RPort {

2 addr: &mut int<16>, read_val: &int<32>

3}

4struct port WPort {

5 inner: &mut (int<16>, Option<int<32>>)

6}

7

8entity dp_mem(clk: clock) -> (RPort, WPort) {

9 let (r_addr, w) = ...;

10 let w_ports = [inst read_wire(w.inner)];

11 let mem = inst clocked_memory(clk, w_ports);

12 let r_val = inst read_memory(

13 mem, inst read_wire(r_addr)

14 );

15 (RPort(r_addr, &r_val), w)

16}

17

18pipeline(10) reader(clk: clock, r_port: RPort) {

19 reg;

20 set r_port.addr = address;

21 let mem_out = *r_port.read_val;

22 reg; // ...

23}

24

25pipeline(5) writer(clk: clock, w_port: WPort) {

26 reg;

27 set w_port.inner = (address, Some(value))

28 reg; // ...

29}

30

31entity top(clk: clock) -> ... {

32 let (r, w) = inst dp_mem(clk);

33 let reader_out = inst(10) reader(clk, r);

34 let writer_out = inst(5) writer(clk, w);

35 // ...

36}

There are some pitfalls when working with mutable wires: A unit returning a
mutable wire expects a value to be set for that wire, otherwise the value of
the wire may be undefined, or if it is set conditionally, a latch could be
inferred. Similarly, if a wire is driven by multiple drivers, conflicting
values may cause issues. While it is possible to catch this in simulation, it
is better to let the compiler catch such errors.

The solution to this is inspired by [14] which uses affine types to ensure
correct memory access patterns in an HLS tool. Affine types can guarantee that
a value is used at most once, resolving the multiple driver problem, however,
to ensure that all mutable wires are set exactly once, the stronger notion of
linear types[15] is required, and implemented in Spade.

Because Spade describes behavior in a cycle-to-cycle manner, the
implementation of a linear type system is easier than in the general case.
Each resource of linear type must be consumed exactly once each clock cycle. A
value is consumed when it is set using the set statement, or passed to another
unit, which delegates the consumption requirement to that unit.

Linear type checking happens after normal type checking, meaning that the
compiler knows which resources are of linear type and must be checked. For
each expression which produces a resource of linear type, a tree is created
where leaf nodes represent primitive linear types, and non-leaves represent
compound linear types such as tuples or structs. Any statement that aliases a
resource, such as a let-binding of an expression, or field access on a struct,
creates a pointer from the alias to the corresponding tree and node.

When expressions or statements which consume resources, such as the set-
statement, or a resource being passed to a unit are encountered, the nodes
corresponding to the consumed object are marked as consumed. If the node or
its child nodes are already marked as consumed, the resource is used more than
once and an error is thrown.

This ensures that nodes are not used more than once, but does not guarantee
that they are used exactly once. Therefore, at the end of the process, a final
pass goes over all the trees to ensure that each node is consumed. If the
traversal finds an unconsumed leaf, it represents a resource of linear type
which was not set, and an error is reported.

To exemplify the linear type checking algorithm, Fig. 3 shows this process in
action for the code in Listing 6. In the figure, $e_{x}$ represent anonymous
names given to sub-expressions before they are bound to variables. A dashed
node represents it not being consumed, a solid node means it is consumed, and
a crossed node indicates double consumption.

##  V Spade Software

A good set of tools for a language, both for working with the language itself,
and for integration with existing tools can be of huge help for driving
language adoption. First and foremost, the Spade compiler and the language is
built from the start to produce useful and easy to read error messages such as
the one shown in Listing 2, and unhelpful diagnostics from the compiler are
considered bugs. The language is also designed to prefer emitting errors over
potentially surprising behavior, for example, by requiring explicit truncation
after potentially overflowing arithmetic.

Figure 3: The linear type checking procedure for the code in Listing 6.
Listing 6: Example code for Fig. 3.

⬇

1 let mem: MemPort = memory();

2 let addr = mem.addr;

3 set addr = value;

4 consume(mem);

###  V-A Compiler Architecture

The Spade compiler is a multi-stage compiler written in Rust, which compiles
the input Spade code to a target language, currently a slim subset of
SystemVerilog. The compilation process starts with lexing and parsing to
generate an Abstract Syntax Tree (AST). The AST is then traversed thrice,
first to collect all types in the program, then to collect all units, and
finally to be lowered into a High-Level Intermediate Representation (HIR). The
AST to HIR lowering process retains the tree structure of the AST but resolves
names and scoping rules, and performs initial semantic analysis. Once the HIR
is generated, type inference is performed, followed by linear type checking as
discussed earlier. The HIR along with the type information is used to generate
a Mid-Level Intermediate Representation (MIR). In this step, more semantic
analysis is performed, and the tree structure is flattened to a list of simple
statements. Finally, SystemVerilog is generated from the MIR. SystemVerilog is
chosen as the target language as it is well-supported by both open source and
proprietary simulation and synthesis tools, but the compiler is written in a
target independent way to enable experimenting with, or changing to a
different backend with ease. Especially interesting backend are the CIRCT[16]
dialects, such as LLHD[17] or Calyx[18], which can offer language independent
optimization as well as code generation of other output languages than
SystemVerilog.

###  V-B Tooling and Ecosystem

Spade comes with the build tool Swim222https://gitlab.com/spade-lang/swim
which manages project files, backend build tools, and dependencies. A Swim
project consists of Spade files and a configuration file written in TOML. This
file contains, among other things: build tool parameters, raw Verilog files to
be included in the project, and external dependencies. Swim then manages
namespacing of project files, downloading and versioning of dependencies as
well as calling the synthesis and simulation tools with a convenient
interface. It also has a plugin system for extending the build flow, for
example by running commands to generate Spade or Verilog code, loading
additional Yosys plugins or bundling the output bitstream into an executable
for a microcontroller which in turn programs a target FPGA.

To facilitate integration of Spade with existing projects, units can be
annotated to prevent name mangling. Spade units can also be marked as
external, in order to allow use of existing IP blocks within Spade projects.

A tree-sitter grammar, and a rudimentary Language Server Protocol (LSP)[19]
server is available, enabling an IDE-like experience in any text editor
supporting LSP and/or tree-sitter.

There is limited effort to generate Verilog similar in structure to the input
Spade code. However, the compiler does attempt to keep names readable, and has
functions for mapping names and expressions back to their source location to
aid debugging. Swim automatically translates values in VCD files into their
high-level Spade values, and the compiler includes a source mapping in the
output Verilog which makes things like timing reports readable without looking
at the output Verilog.

###  V-C Test Benches and Simulation

Listing 7: Example of a test bench for Spade written using cocotb.

⬇

1 # top=peripherals::timer::timer_test_harness

2 @cocotb.test()

3 async def timer_works(dut):

4 s = SpadeExt(dut)

5

6 clk = dut.clk_i

7 await start_clock(clk)

8

9 s.i.mem_range = "(1024, 2048)"

10 s.i.addr = "1024 + 0"

11 s.i.memory_command = "Command::Write(10)"

12 await FallingEdge(clk)

13 s.o.assert_eq("10")

The Spade language itself is designed primarily for hardware description and
synthesis, rather than simulation. However, in order to verify correctness of
the resulting hardware, simulation and test benches are essential. Spade tests
are written using cocotb[20], a Python based co-simulation test bench
environment for verification. The cocotb library is extended with features for
writing Spade values as inputs and outputs to the unit under test, while
Verilog generated by the compiler is simulated with off-the-shelf Verilog
simulators.

As an example, Listing 7 shows part of a test bench for a memory mapped timer
peripheral. The first line specifies the module being tested, and the next two
lines are a standard cocotb test case definition. Line 4 wraps the cocotb
design under test in a Spade class which extends the cocotb interface to add
Spade-specific features. Lines 6–7 start a task to drive the clock of the
Design Under Test (DUT). On lines 9–11, the inputs to the module are set to
values which are compiled and evaluated by the Spade compiler. Finally, line
13 asserts that the output is as expected, again passing a string containing a
Spade expression as the expected value.

In order to achieve this, the Python test bench must be able to compile Spade
code, and in order to allow the use of types and units defined in the project
inside a test bench, the state of the compiler must be made available to the
test bench. For this reason, the state of the compiler after building a
project is serialized and stored on disk. Additionally, parts of the Spade
compiler are exported as a Python module which reads the stored state, then
compiles and evaluates the expressions used in the test bench. The resulting
bit vectors are used to drive the inputs of the DUT, or compared to the
expected output.

##  VI Conclusions

Spade is an HDL which attempts to ease the development of FPGAs or ASICs. To
do so, it makes common hardware constructs like pipelines, registers, and
memories explicit and part of the language. It is heavily inspired by modern
software languages, e.g., by integrating a powerful type system and pattern
matching. Unlike most current alternative HDLs, Spade is its own standalone
language with a custom compiler, and does not abstract away the underlying
hardware. Finally, it comes with useful tooling, such as a compiler designed
to emit detailed error messages, a build system with dependency management,
and an LSP for an IDE-like experience.

## References

  * [1] J. Bachrach _et al._ , “Chisel: Constructing hardware in a Scala embedded language,” in _Proc. Des. Automat. Conf._ , Jun. 2012, pp. 1212–1221.
  * [2] SpinalHDL contributors, “SpinalHDL,” https://github.com/SpinalHDL/SpinalHDL, 2022.
  * [3] O. Port and Y. Etsion, “DFiant: A dataflow hardware description language,” in _Proc. Int. Conf. Field Programmable Logic Appl._ , 2017.
  * [4] Amaranth contributors, “Amaranth HDL,” https://github.com/amaranth-lang/amaranth, 2022.
  * [5] M. Korbel, “Rapid open hardware development framework,” in _Proc. Workshop Open-Source EDA Technol._ , 2022.
  * [6] J.-C. Le Lann, H. Badier, and F. Kermarrec, “Towards a hardware DSL ecosystem: RubyRTL and friends,” in _Proc. Des. Automat. Test Europe Conf. – Open Source Des. Automat. Workshop._ , Mar. 2020.
  * [7] C. Baaij, “Digital circuits in c$\lambda$aSH,” PhD. Thesis, University of Twente, Jan. 2015.
  * [8] J. Kemmerer, “PipelineC: Easy open-source hardware description between RTL and HLS,” in _Proc. Workshop Open-Source EDA Technol._ , 2022.
  * [9] S. F. Hoover, “Timing-abstract circuit design in transaction-level Verilog,” in _Proc. IEEE Int. Conf. Comput. Des._ IEEE, Nov. 2017.
  * [10] R. Nikhil, “Bluespec SystemVerilog: efficient, correct RTL from high-level specifications,” in _Proc ACM/IEEE Int. Conf. Formal Methods and Models for Co-Design_. IEEE, 2004.
  * [11] S. Lefebvre, “Silice,” https://github.com/sylefeb/Silice/tree/5003ec72, Nov. 2022.
  * [12] S. Klabnik and C. Nichols, _The Rust programming language_ , 1st ed. No Starch Press, Jun. 2018, ch. Glossary.
  * [13] D. Luu, “Writing safe Verilog,” https://danluu.com/pl-troll/.
  * [14] R. Nigam _et al._ , “Predictable accelerator design with time-sensitive affine types,” in _Proc. ACM SIGPLAN Conf. Program. Language Des. Implementation_. ACM, Jun. 2020.
  * [15] P. Wadler, “Linear types can change the world!” in _Programming Concepts and Methods_ , 1990.
  * [16] CIRCT Project, “CIRCT,” https://circt.llvm.org/, 2022.
  * [17] F. Schuiki _et al._ , “LLHD: a multi-level intermediate representation for hardware description languages,” in _Proc. ACM SIGPLAN Conf. Program. Language Des. Implementation_ , Jun. 2020.
  * [18] R. Nigam _et al._ , “A compiler infrastructure for accelerator generators,” in _Proc. ACM Int. Conf. Arch. Support for Programming Lang. Operating Syst._ , Apr. 2021.
  * [19] Microsoft, “Language server protocol,” https://microsoft.github.io/language-server-protocol/.
  * [20] cocotb contributors, “cocotb,” https://www.cocotb.org/, Nov. 2022.

  *[HDL]: Hardware Description Language
  *[HDLs]: Hardware Description Language
  *[DSL]: Domain Specific Language
  *[HLS]: High-Level Synthesis
  *[AST]: Abstract Syntax Tree
  *[HIR]: High-Level Intermediate Representation
  *[MIR]: Mid-Level Intermediate Representation
  *[LSP]: Language Server Protocol
  *[DUT]: Design Under Test