# Building an Efficient Cluster Cosmology Software Package for Modeling
Cluster Counts and Lensing

M. Aguena1 O. Alves2 J. Annis3 D. Bacon4 S. Bocquet5 D. Brooks6 A. Carnero
Rosell1,7,8 C. Chang9,10 M. Costanzi11,12,13 C. Coviello3 L. N. da Costa1 T.
M. Davis14 J. De Vicente15 H. T. Diehl3 P. Doel6 J. Esteves2 S. Everett16 I.
Ferrero17 A. Ferté18 D. Friedel19 J. Frieman9,3 M. Gatti20 G. Giannini21 D.
Gruen5 R. A. Gruendl22,19 G. Gutierrez3 K. Herner3 S. R. Hinton14 D. L.
Hollowood23 K. Honscheid24,25 D. J. James26 T. Jeltema23 M. Kirby27 K.
Kuehn28,29 O. Lahav6 P. Li30 J. L. Marshall31 T. McClintock27 D. Mellor32 J.
Mena-Fernández15 R. Miquel33,21 J. O’Donnell23,9,3 A. Palmese30 M. Paterno3 M.
E. S. Pereira34 A. Pieres1,35 A. A. Plazas Malagón36,18 M. Rodriguez-Monroy15
A. K. Romer37 A. Roodman18,36 E. Sanchez15 M. Schubnell2 I. Sevilla-Noarbe15
T. Shin38 M. Smith39 E. Suchyta40 M. E. C. Swanson41 G. Tarle2 J. Weller42,43
P. Wiseman39 H.-Y. Wu44 Y. Zhang45 C. Zhou23,46
(DES Collaboration)

###### Abstract

We introduce a software suite developed for galaxy cluster cosmological
analysis with the Dark Energy Survey Data. Cosmological analyses based on
galaxy cluster number counts and weak-lensing measurements need efficient
software infrastructure to explore an increasingly large parameter space, and
account for various cosmological and astrophysical effects. Our software
package is designed to model the cluster observables in a wide-field optical
survey, including galaxy cluster counts, their averaged weak-lensing masses,
or the cluster’s averaged weak-lensing radial signals. To ensure maximum
efficiency, this software package is developed in C++ in the CosmoSIS software
framework, making use of the CUBA integration library. We also implement a
testing and validation scheme to ensure the quality of the package. We
demonstrate the effectiveness of this development by applying the software to
the Dark Energy Survey Year 1 galaxy cluster cosmological data sets, and
acquired cosmological constraints that are consistent with the fiducial Dark
Energy Survey analysis.

##  1 Introduction

The abundance of galaxy clusters, the most massive gravitationally-bound
structures in the Universe, are important cosmology probes utilized by wide-
field cosmic surveys like the Dark Energy Survey (DES) (Abbott et al., 2020;
To et al., 2021a). Together with large scale structure correlation functions
(Frieman et al., 2008; Weinberg et al., 2013), galaxy cluster cosmology
studies are projected to yield significantly improved constraints on dark
energy in coming years, especially with the Rubin Observatory Legacy Survey of
Space and Time (LSST) (The LSST Dark Energy Science Collaboration et al.,
2018).

To derive cosmological constraints from galaxy cluster abundance observables
(see a review in Allen et al., 2011), studies rely on a few theoretical
ingredients to model the observations. First, gravity theories predict that
galaxy clusters form through linear gravitational accretion and nonlinear
collapse (see a review in Kravtsov & Borgani, 2012). Different cosmology
models predict different abundances of galaxy clusters, and those abundances
can be modeled with a dark matter halo mass function (e.g., Press & Schechter,
1974; Tinker et al., 2008; McClintock et al., 2019a). This mass function
depends on cosmological parameters such as the matter fluctuation amplitude
and its growth, $\sigma_{8}$, and the total matter density, $\Omega_{m}$, and
is subject to effects such as the number of neutrino species and the neutrino
matter densities (Costanzi et al., 2013), as well as baryonic effects at the
small scales (Bocquet et al., 2016). Further, the cosmic geometry predicted by
cosmological models allows us to translate a dark matter halo mass function
into dark matter halo counts (Weinberg et al., 2013).

To compare the abundance of dark matter halos to galaxy cluster observations,
we also need to forward model how imaging surveys identify and measure galaxy
clusters using a set of observational criteria (e.g., Rykoff et al., 2014,
2016; Palmese et al., 2020; Aguena et al., 2021). The key to this task is a
cluster observable-mass relation which models the relation between a cluster
mass “proxy” and the dark matter halo mass. This relation takes many forms in
literature – the Log-Normal relation is typical, but new forms of the relation
have been suggested in recent years (e.g., Mao et al., 2015; Costanzi et al.,
2019a; Zhang & Annis, 2022).

Parameter inference in cluster cosmology is often done by adopting a Bayesian
framework comparing theoretical predictions to observations. The posterior
distributions of the cosmological parameters and nuisance parameters are often
sampled using Markov Chain Monte Carlo (MCMC) methods, with a typical analysis
requiring hundreds to thousands of CPU hours (Lemos et al., 2022). Additional
systematic effects, such as incorrectly choosing a cluster’s central galaxy
known as “mis-centering” (Zhang et al., 2019), the projection of correlated
and uncorrelated structures around the clusters (Costanzi et al., 2019a), the
weak lensing cluster member contamination (Varga et al., 2019), and the shape
and orientation of real non-spherical clusters also need to be considered as
corrections to the theoretical expectation from spherically symmetric cluster
density distributions. Given all of those elements, forward-modeling the
galaxy cluster abundance and their shear measurements can be approached as an
integration problem because of the need to marginalize over various
theoretical and observational effects. Solving this problem poses a
significant computational challenge when the integrations have large
dimensions and need to be computed at every step of a Monte Carlo parameter
sampling analysis.

In this paper, we introduce our efforts to develop a software package to
compute models for cluster abundance and weak lensing measurements. The latter
introduces a few more dimensions of integrations than the fiducial analysis
adopted in DES Y1 (Abbott et al., 2020). As the MCMC sampling methods used for
cosmological parameter inference often need to compute those models with
varying cosmological and nuisance parameters for many thousands of times, the
efficiency of those computations becomes critically important. The change
necessitates the development of a new pipeline that is focused on computing
efficiency and the application on high-performance computing resources. In the
rest of this paper, we describe the design background of this new pipeline
including the theoretical models, the structure of our software package and
our testing mechanisms. We also use this new package to reproduce the fiducial
DES Y1 analysis in Section 6 and acquire consistent cosmological results. This
analysis also serve as an independent check of the fiducial DES Y1 analysis.

The package described in this paper will be made openly available.

##  2 Cluster Cosmology Theoretical Framework

###  2.1 Number Counts and Average Masses

In galaxy cluster cosmology studies performed by optical surveys like Sloan
Digital Sky Survey (SDSS) and DES, the clusters are discovered through their
richness, denoted as $\lambda$ (Rykoff et al., 2012). To-date, the cluster
cosmology analyses performed with large cluster samples from SDSS and DES
(Rozo et al., 2010; Costanzi et al., 2019b; Abbott et al., 2020) use summary
statistics of the selected samples, i.e. the number of clusters in a few
richness and redshift ranges, and their average masses in these
richness/redshift ranges as the observable of cosmology. In this section, we
describe the analytical formalism for modeling those cluster number and mass
observables.

We first present the theoretical predictions for the total number and the
average mass of the clusters in one redshift and richness bin. Those
quantities can be written as

$\begin{split}\left<N(\Delta\lambda_{i},\Delta z_{j})\right>=&\int_{\Delta
z_{j}}dz^{\text{ob}}\int_{\Delta\lambda_{i}}d\lambda^{\text{ob}}\int_{0}^{\infty}dM\int_{0}^{\infty}dz^{\text{true}}\\\
&P(z^{\text{ob}}|z^{\text{true}},\Delta\lambda_{i})P(\lambda^{\text{ob}}|M,z^{\text{true}})\\\
&n(M,z^{\text{true}})\
\Omega(z^{\text{true}})\frac{dV}{d\Omega\,dz^{\text{true}}}(z^{\text{true}})\
.\\\ \left<NM(\Delta\lambda_{i},\Delta z_{j})\right>=&\int_{\Delta
z_{j}}dz^{\text{ob}}\int_{\Delta\lambda_{i}}d\lambda^{\text{ob}}\int_{0}^{\infty}dM\int_{0}^{\infty}dz^{\text{true}}\\\
&P(z^{\text{ob}}|z^{\text{true}},\Delta\lambda_{i})P(\lambda^{\text{ob}}|M,z^{\text{true}})\\\
&M\cdot n(M,z^{\text{true}})\
\Omega(z^{\text{true}})\frac{dV}{d\Omega\,dz^{\text{true}}}(z^{\text{true}})\
.\\\ \left<M(\Delta\lambda_{i},\Delta
z_{j})\right>=&\left<NM(\Delta\lambda_{i},\Delta
z_{j})\right>/\left<N(\Delta\lambda_{i},\Delta z_{j})\right>.\end{split}$ (1)

In this set of equations, $\left<N(\Delta\lambda_{i},\Delta z_{j})\right>$ is
the predicted number of galaxy clusters in a redshift range of $\Delta z_{j}$
and a richness range of $\Delta\lambda_{i}$, while
$\left<NM(\Delta\lambda_{i},\Delta z_{j})\right>$ is the predicted total mass
of those clusters and $\left<M(\Delta\lambda_{i},\Delta z_{j})\right>$ is
their average mass.

The predictions for those three quantities consider several theoretical
elements. First, the volume density of galaxy clusters in terms of their
theoretical quantities, mass $M$ and redshift $z^{\text{true}}$, is known as
the halo mass functions $n(M,z^{\text{true}})$ . Analytical approximations for
the halo mass function can be calculated from the matter power spectrum and
the cosmological model, and depends most strongly on the values of
$\Omega_{m}$ and $\sigma_{8}$ in a $\mathrm{\Lambda CDM}$ cosmology (our
applications are based upon Tinker et al., 2008). Second,
$\Omega(z^{\text{true}})$ is the total solid angle area of the cosmic survey,
which may have a redshift dependence as the usable area for cluster searching
changes with redshift. Further,
$\frac{dV}{d\Omega\,dz^{\text{true}}}(z^{\text{true}})$ is the cosmological
volume element in solid angle $d\Omega$ and a redshift slice
$dz^{\text{true}}$, which also depends on the cosmological parameters. On top
of those cosmology and geometry related quantities, we also consider a few
relations linking cluster observable to their theoretical quantities. Among
those, $P(z^{\text{ob}}|z^{\text{true}},\Delta\lambda_{i})$ considers the
redshift uncertainty of the clusters, as its measured value of $z^{\text{ob}}$
may deviate from its true redshift characterized by a (sometimes Gaussian)
probabilistic distribution. $P(\lambda|M,z^{\text{true}})$ relates the
cluster’s observational selection quantity $\lambda$ to the cluster’s
theoretical mass quantity, $M$, and true redshift $z^{\text{true}}$.
$P(\lambda|M,z^{\text{true}})$ can be further expanded to marginalize over
latent parameters to consider astrophysical effects like projection and
cluster triaxial shapes. A more detailed description can also be found in
Abbott et al. (2020).

###  2.2 Number Counts and Average Lensing Signals

In observations, to measure the averaged cluster masses, the averaged cluster
lensing signals are analyzed for the clusters in the richness ranges
(McClintock et al., 2019b). In this software development, our end goal is to
forward model this lensing observable together with the cluster number counts,
instead of relying on the cluster’s average mass measurements and number
counts. Realizing this goal requires us to further consider the relations
between cluster mass, concentration ($c$), and the cluster’s projected radial
mass density from the center $\Sigma_{cen}(r|M,z^{\text{true}},c)$.

$\displaystyle\left<N\Sigma_{\mathrm{cen}}(r,\Delta\lambda_{i},\Delta
z_{j},c)\right>=$ (2) $\displaystyle\int_{\Delta
z_{j}}dz^{\text{ob}}\int_{\Delta\lambda_{i}}d\lambda\int_{0}^{\infty}dM\int_{0}^{\infty}dz^{\text{true}}$
$\displaystyle
P(z^{\text{ob}}|z^{\text{true}},\Delta\lambda_{i})P(\lambda|M,z^{\text{true}})$
$\displaystyle\Sigma_{cen}(r|M,z^{\text{true}},c)\cdot n(M,z^{\text{true}})\
\Omega(z^{\text{true}})\frac{dV}{d\Omega dz^{\text{true}}}(z^{\text{true}})$

Further, those projected mass densities are susceptible to an additional
physical effect – mis-identification of galaxy cluster centers, that the
cluster’s real centers can be offset from the identified centers by a
distance, $R_{\text{cen}}$. This mis-centering offset can be described by a
probability distribution $P(R_{\text{cen}})$:

$\begin{split}&P(R_{\text{cen}})=f_{\text{cen}}\delta(R_{\text{cen}})+(1-f_{\text{cen}})P_{mis}(R_{\text{cen}}).\\\
\end{split}$ (3)

In this equation, we recognize that a significant fraction ($f_{\text{cen}}$)
of the clusters are actually well-centered. For those well-centered clusters,
$R_{\text{cen}}$ can be described as a $\delta$ function, as often adopted in
cluster lensing modeling (e.g., McClintock et al., 2019b). For the clusters
that are indeed mis-centered, their mis-centering distance can be described by
a continuous probabilistic function $P_{mis}(R_{\text{cen}})$, constrained by
multi-wavelength studies (e.g., Zhang et al., 2019). Further, mis-
identification of cluster centers have also been shown to further affect the
cluster richness observable by introducing dependence on the cluster’s mis-
centering offset, which can be modeled by an additional latent parameter
$\lambda^{\text{cen}}$ describing a cluster’s would-be richness without the
mis-centering effect (Johnston et al., 2007; Zhang et al., 2019).

$\begin{split}&P(\lambda^{\text{ob}}|M,z^{\text{true}},R_{\text{cen}})=\\\
&\int_{0}^{\infty}d\lambda^{\text{cen}}\,P(\lambda^{\text{ob}}|\lambda^{\text{cen}},R_{\text{cen}})P(\lambda^{\text{cen}}|M,z^{\text{true}}).\end{split}$
(4)

This equation marginalizes over the latent $\lambda^{\text{cen}}$ parameter,
assuming that the cluster’s observed richness $\lambda^{\text{ob}}$ depends on
$\lambda^{\text{cen}}$ and the cluster’s mis-centering offset $R_{\text{cen}}$
with a probability distribution
$P(\lambda^{\text{ob}}|\lambda^{\text{cen}},R_{\text{cen}})$.

Thus, the prediction for the clusters’ projected surface mass density as a
function of radius can be written as,

$\displaystyle\left<\Sigma(r|\Delta\lambda_{i},\Delta z_{j},c)\right>=$ (5)
$\displaystyle\frac{f_{\text{cen}}}{\left<N(\Delta\lambda_{i},\Delta
z_{j})\right>}\left<N\Sigma_{\mathrm{cen}}(r,\Delta\lambda_{i},\Delta
z_{j},c)\right>\ +$
$\displaystyle\frac{1-f_{\text{cen}}}{\left<N(\Delta\lambda_{i},\Delta
z_{j})\right>}\left<N\Sigma_{\mathrm{mis}}(r,\Delta\lambda_{i},\Delta
z_{j},c)\right>\ .$

While $\left<N\Sigma_{\mathrm{cen}}(r,\Delta\lambda_{i},\Delta
z_{j},c)\right>$ has been defined in Equation2,
$\left<N\Sigma_{\mathrm{mis}}(r,\Delta\lambda_{i},\Delta z_{j},c)\right>$
needs to azimuthally average over the mis-centering distance distribution,
which can be written as:

$\displaystyle\left<N\Sigma_{\mathrm{mis}}(r,\Delta\lambda_{i},\Delta
z_{j},c)\right>=$ (6) $\displaystyle\int_{\Delta
z_{j}}dz^{\text{ob}}\int_{\Delta\lambda_{i}}d\lambda^{\text{ob}}\int_{0}^{\infty}dM\int_{0}^{\infty}dz^{\text{true}}$
$\displaystyle\int_{0}^{\infty}d\lambda^{\text{cen}}\int_{0}^{\infty}dR_{\text{cen}}\int_{0}^{2\pi}\frac{d\alpha}{2\pi}$
$\displaystyle P(z^{\text{ob}}|z^{\text{true}},\Delta\lambda_{i})\
P(\lambda^{\text{ob}}|\lambda^{\text{cen}},R_{\text{cen}})\ $ $\displaystyle
P(\lambda^{\text{cen}}|M,z^{\text{true}})\ P_{mis}(R_{\text{cen}})$ (7)
$\displaystyle n(M,z^{\text{true}})\ \Omega(z^{\text{true}})\frac{dV}{d\Omega
dz^{\text{true}}}(z^{\text{true}})$
$\displaystyle\Sigma_{cen}(\sqrt{r^{2}+R_{\text{cen}}^{2}-2rR_{\text{cen}}.\mathrm{cos}\alpha}|M,z^{\text{true}},c).$

Given that mis-centering affects the cluster richness-mass relation, we shall
also update the writing of the cluster number counts as

$\displaystyle\left<N(\Delta\lambda_{i},\Delta z_{j})\right>=\ $
$\displaystyle
f_{\text{cen}}\cdot\left<N_{\mathrm{cen}}(\Delta\lambda_{i},\Delta
z_{j})\right>$ (8)
$\displaystyle+(1-f_{\text{cen}})\left<N_{\mathrm{mis}}(\Delta\lambda_{i},\Delta
z_{j})\right>,$ $\displaystyle\left<N_{\mathrm{cen}}(\Delta\lambda_{i},\Delta
z_{j})\right>$ $\displaystyle=$ $\displaystyle\int_{\Delta
z_{j}}dz^{\text{ob}}\int_{\Delta\lambda_{i}}d\lambda^{\text{ob}}\int_{0}^{\infty}dM\int_{0}^{\infty}dz^{\text{true}}$
$\displaystyle P(z^{\text{ob}}|z^{\text{true}},\Delta\lambda_{i})\
P(\lambda^{\text{ob}}|M,z^{\text{true}})$ $\displaystyle n(M,z^{\text{true}})\
\Omega(z^{\text{true}})\frac{dV}{d\Omega\,dz^{\text{true}}}(z^{\text{true}}),$
$\displaystyle\left<N_{\mathrm{mis}}(\Delta\lambda_{i},\Delta z_{j})\right>$
$\displaystyle=$ $\displaystyle\int_{\Delta
z_{j}}dz^{\text{ob}}\int_{\Delta\lambda_{i}}d\lambda^{\text{ob}}\int_{0}^{\infty}dM\int_{0}^{\infty}dz^{\text{true}}$
$\displaystyle\int_{0}^{\infty}d\lambda^{\text{cen}}\int_{0}^{\infty}dR_{\text{cen}}$
$\displaystyle P(z^{\text{ob}}|z^{\text{true}},\Delta\lambda_{i})\
P(\lambda^{\text{ob}}|\lambda^{\text{cen}},R_{\text{cen}})\ $ $\displaystyle
P(\lambda^{\text{cen}}|M,z^{\text{true}})\ P_{mis}(R_{\text{cen}})$
$\displaystyle n(M,z^{\text{true}})\
\Omega(z^{\text{true}})\frac{dV}{d\Omega\,dz^{\text{true}}}(z^{\text{true}})\
.$

Finally, the results of the cluster’s averaged projected densities can be
further converted into the so-called excess surface mass density
$\Delta\Sigma$ as with:

$\left<\Delta\Sigma(R|\Delta\lambda_{i},\Delta z_{j},c)\right>=\\\
\frac{2}{R^{2}}\int_{0}^{R}r\left<\Sigma(r|\Delta\lambda_{i},\Delta
z_{j},c)\right>\mathrm{d}r-\left<\Sigma(R|\Delta\lambda_{i},\Delta
z_{j},c)\right>.$ (9)

or into reduced tangential shear signals, $\gamma_{t}$, which are observables
from cluster weak lensing analyses (Umetsu, 2020).

##  3 Parameter Estimation and CosmoSIS

Figure 1: Components of CosmoSIS pipeline using a Markov Chain Monte Carlo
sampling method to derive the posterior cosmological constraints.

Our development is guided by the Bayes’ theorem to derive a posterior
probability distribution of a set of model parameters $\theta$, upon a set of
observations $x$.

$P(\theta\,|\,x,I)\propto P(x\,|\,\theta,I)\,P(\theta\,|\,I),$ (10)

In this equation, our knowledge of the parameters before consideration of the
data $x$, but based on our prior information $I$, is described by the _joint
prior probability_ $P(\theta\,|\,I)$. On the other hand, $P(x\,|\,\theta,I)$,
often called the _likelihood_ , is the probability of observing the data $x$,
given the parameters $\theta$ and the same prior information $I$.

In cluster cosmology studies that rely on summary statistics, such as the
cluster number counts and their lensing masses or their average lensing
profiles described in the previous section, the likelihood is often assumed to
be a multi-dimensional Gaussian as studied and tested in literature (Gruen et
al., 2015; Wu et al., 2019; Fumagalli et al., 2021). The mean of the Gaussian
is the expectation values of the predictions for the data in a given cosmology
($\mu(\theta)$) while the spread of the distribution is described by a
covariance matrix $C$ that characterizes the uncertainties in both observation
and theory. Thus the likelihood function can be written as:

$P(x\,|\,\theta,I)=\mathcal{N}(x-\mu(\theta),C(\theta))$ (11)

To derive the posterior cosmological constraints given the observational
measurements, theoretical predictions, and the covariance matrix, we rely on
MCMC methods, as well as other sampling techniques. These samplings and the
development of the prediction code are performed in the CosmoSIS software
framework 111https://bitbucket.org/joezuntz/cosmosis. Figure 1 demonstrates a
typical workflow to derive posterior cosmological constraints in CosmoSIS.

CosmoSIS Zuntz et al. (2015) is a modular system for cosmological parameter
estimation. It solves integrals like that in Equation 10 using MCMC. MCMC
solves an $n$-dimensional integral by generating a collection of points,
called _samples_ , in an $n$-dimensional space; the points are distributed
according to the functional form of the solution. A _sampler_ drives the
exploration of the parameter space. A sampler is the software component
responsible for determining the location in the parameter space for the next
step of the MCMC chain based upon the evaluation of the user-supplied
likelihood function, and the location of the current sample. The user provides
the likelihood function by supplying one or more modules (_Forward Model_
modules and/or _Gaussian Likelihood_ modules in Figure 1) which, given a point
in the parameter space, evaluate the likelihood of the data given those
parameters.

In CosmoSIS, a _module_ is user-written code responsible for calculating
physical quantities of interest (e.g., predicted cluster number counts) based
on the parameters of the current sample. CosmoSIS modules can be written in
several programming languages, including Python, C, Fortran, C++ and Julia.
Calculations from those modules are stored in a structure known as a DataBlock
object, which is passed through the list of modules used in a CosmoSIS
pipeline setup, and thus functions to pass results between modules.

##  4 Pipeline structure

Our pipeline delivers prediction code that can be implemented in a CosmoSIS
_module_ for the likelihood calculation (Equation 11), $\mu(\theta)$. In
particular, we calculate the predictions for the cluster number counts,
cluster masses and/or shear profiles using models and the integrals they
involve. The efficiency of the numerical integration library we use is of
primary importance, as the calculation of the integrals of interest can
involve $10^{5}$—$10^{6}$ function evaluations. In the rest of the section, we
describe the structures of this pipeline as well as how those structures
function in an analysis implementation.

###  4.1 “Models”

The smallest software units in our development of the pipeline are the _model_
s. Examples include the survey area function $\Omega(z^{\text{true}})$, the
halo mass function $n(M,z^{\text{true}})$, and the cluster observable-mass
relation $P(\lambda|M,z^{\text{true}})$ from Equation 1.

In each case, the essence of one _model_ is a function which we need to
evaluate as part of our integrand, and which can be tested for correctness in
isolation from the rest of the integrand. Furthermore, there may be more than
one sensible implementation of the function—for example, in the case of the
survey area function, the implementation would be different for different
surveys. Finally, the evaluation of the function may require the use of some
additional data, other than the values of the function arguments. A very
common example of this is the use of interpolation tables, which can be output
arrays from other CosmoSIS modules (such as the CosmoSIS CAMB module), or
derivations from supporting analyses not as part of this efforts.

Figure 2 shows most of the implementation of one of our _model_ s, the survey
area for the SDSS redMaPPer cluster sample, which is redshift dependent. The
name of the class is OMEGA_Z_SDSS. The state of each object of this type
consists only of a polynomial object named fit_, which itself is a callable
object representing a polynomial interpolation function. The _default
constructor_ of the class, which initializes an object of a class, establishes
the state of the object, in this case a polynomial used to evaluate the
function. There is a second constructor, which takes a CosmoSIS DataBlock
object; this use of the constructor will be described later. In this case, the
constructor does nothing differently than the default constructor—it is
acceptable for a constructor to ignore its argument(s). The function call
operator shows that this class represents a function of one argument, zt, and
the implementation is using the interpolation table to calculate the function.
Importantly, the function call operator is implemented _inline_. This means
the compiler is free to put the code of the function directly into the place
where the function is called. This is important because it provides the
compiler the opportunity to see the implementation of the function at any
place in the code where this function is called, which allows an optimizing
compiler to make use of this information to do more optimizations than would
be possible without this information. This member function is const-qualified,
which means calling the function does not alter the state of the object.
Marking the function helps assure the correctness of our code; if an error in
implementation is made that perhaps inadvertently modifies the state of the
object (for example, by changing the values in an interpolation table), a
compilation failure will occur. Marking the function const can also provide a
compiler additional opportunities for optimization of the generated machine
code.

The example class shown in Figure 2 computes the survey area for SDSS, but
additional classes/_model_ s that represent the survey area of different data
sets, e.g., DES year 1 data sets, or year 6 data sets can be constructed. Each
class that represents a survey area has a function call operator with the same
signature: one argument, of type double representing the redshift at which the
area element is to be calculated, and returning a double representing the area
of the survey at that redshift.

In addition, every _model_ (not just those representing cosmological area
calculations) has a constructor that takes a CosmoSIS DataBlock object. This
constructor is used when the _model_ is invoked in a CosmoSIS pipeline, and is
used to provide any necessary initialization values to the _model_. In this
case, no such initialization parameters are needed, and this constructor does
not need to consult the DataBlock which is given. However, the constructor
must still be present, so that the generic code that creates a _model_ (and
which provides the DataBlock to be used) is well-formed. Our library of models
includes two volume elements, one for the DES and one for the SDSS.

⬇

private:

polynomial<12> fit_{...};

public:

OMEGA_Z_SDSS() = default;

OMEGA_Z_SDSS(cosmosis::DataBlock&) {}

double

operator()(double zt) const

{

return fit_(zt - 0.2);

}

};

Figure 2: Abridged implementation of the class OMEGA_Z_SDSS, which represents
the cosmological survey area for SDSS. Details of namespaces and the
implementation of the initialization of fit_ have been suppressed for brevity.

We also require that those _model_ s, as they are the smallest development
unit, must have a corresponding test whenever possible. This is described in
more detail in Section 5.1.

###  4.2 CosmoSIS modules in general

The models are embedded into CosmoSIS modules. The module interface provided
by CosmoSIS has a natural implementation when translated into C++. Each module
is represented as an instance of a _module class_. We note that no inheritance
is necessary, because there is nothing in CosmoSIS that would ever interact
with the base class of an inheritance hierarchy.

A module object is created at program start-up time, using the constructor of
the class. Any configurable parameters required by the class are passed to the
constructor through a DataBlock object, which contains the configuration
information read by the CosmoSIS program. At each step of the MCMC process,
CosmoSIS calls the module’s _execute_ function; this is naturally implemented
in C++ as a member function execute, which is passed an instance of the class
_DataBlock_ carrying the cosmological parameters and derived physical
quantities for the current MCMC sample. For our prediction modules, these
quantities are used as inputs for calculating the predictions of observables,
given the cosmology of the current sample. The prediction modules in turn
record these predictions in the DataBlock for later processing by our
likelihood module. Our likelihood module retrieves the predictions from the
DataBlock and calculates the data likelihood for the sample cosmology,
recording it in the DataBlock. Figure 3 shows a flowchart of the operations
between integrands, modules and the DataBlock described in this section.
Finally, the module objects are destroyed at program shutdown. In a well-
designed class, there is usually no user-written code necessary for the
destructor; the compiler-generated destructor will assure that all allocated
resources are properly released.

###  4.3 Our physics modules

The prediction modules calculate their predictions based on integrands formed
from the models described in Section 4.1. These models are typically created
for a given cosmology, and thus model instances are created anew for each MCMC
sample. (In some cases, as described below, models have an even shorter
lifetime.) Equation 2 shows one such integral, which calculates the centered
shear profile. This integral is done over five variables of integration. Two
of the variables of integration ($z^{\text{ob}}$, $\lambda^{\text{ob}}$) are
each integrated over several distinct, possibly overlapping, ranges. Thus
there are distinct 5-dimensional volumes of integration to be calculated for
each sample. In addition, the integral is evaluated for a distinct set of
values of $r$ and $c$, which form a grid of locations in $(r,c)$ space at
which each of the integrals are calculated. Typically, several hundred
integrals are calculated for a single integrand—the product of the number of
volumes of integration and the number of grid points at which the integral is
evaluated.

Figure 3: Flowchart of the CosmoSIS cluster pipeline summarizing the
operations between integrands, modules and DataBlock.

###  4.4 The physicist supplied integration problem class

To simplify the task of writing a prediction module, we have created the class
template CosmoSISScalarIntegrand that takes a user provided template parameter
and supplies the integrand. The template provides several facilities. It
configures the model upon construction, when possible. It re-configures the
model for each sample, if needed. It re-configures the model for the grid
points at which the integrands are integrated. It configures the volumes of
the integrand, and calculates all the integration results. It then organizes
the integration results into simple data structures and puts them into the
DataBlock. The type used for this template argument must be a class (or
struct).

We have developed a class template (CosmoSISScalarIntegrationModule) that
employs these tools to make the writing of a physics module that performs
integrals like those described above simpler. The physicist who wishes to
write an integration module can do so with the help of this template, and then
needs to write a simpler class (named X, for example) that provides the
following features; the type CosmoSISScalarIntegrationModule<X> is then a
fully-functional CosmoSIS module.

The required features of the physicist-provided class include not only the
implementation of the function to be integrated (which is implemented through
the function call operator of the class), but also functions that provide the
information needed to make the module configuration code function for the
specific integrand.

The physicist inventing a new module must write a class with several features
that are relied upon by the class template CosmoSISScalarIntegrationModule.
These features include some required nested types (type definitions within the
class), static and non-static member functions, and static member data. In the
remainder of this subsection we describe these features.

The class must have two nested type declarations, grid_t and volume_t. The
type grid_t is to be defined using the helper template y3_cluster::grid_t<M>,
using the appropriate value of M for the dimensionality of the desired grid
points (as described in Section 4.3, grid points are set of variable values at
which the integrals are evaluated). The type volume_t is to be defined using
the helper template cubacpp::IntegrationVolume<N>, and using the appropriate
value of N to define the dimensionality of the volume of integration
(otherwise also known as the boundaries of the parameters to be integrated).

The class must have a constructor that accepts a single DataBlock object. This
object carries the configuration parameters provided by the CosmoSIS user’s
configuration file. Parameters that do not change with MCMC sample or with the
grid point at which the integral is being evaluated should be set in this
constructor.

The class must provide a function call operator that defines the actual
integrand. This must define a (non-static) member function with a return type
of double and which accepts any fixed number of double arguments. At
compilation time, this member function will be probed to determine the
dimensionality of the resulting integrand, and to verify the dimensionality of
the integration volume specified in the type volume_t.

The class must provide a member function set_sample that accepts a DataBlock
object. This function is called with every MCMC sample; the argument carries
all the data (most importantly, the cosmology) corresponding to the current
sample. The integrand object can read whatever data are needed to adjust model
parameters to reflect the current sample.

The class must provide a static member function make_integration_volumes,
which accepts a DataBlock. This function is to be implemented using one of two
provided helper function templates, make_integration_volumes_wall_of_numbers
and make_integration_volumes_cartesian_product. The choice of helper used by
the integrand class determines how the integration volumes used by resulting
CosmoSIS module are configured. Both helper templates are called with a number
of C-style strings specifying the names of the parameters that are to be read
from the CosmoSIS module configuration; the number of strings specified must
match the declared dimensionality of the integrand’s grid. The compiler will
reject as invalid code any class for which these two do not agree. The values
of these parameters read at the startup of a CosmoSIS run determines the set
of points at which the integrals are to be evaluated. Note that the author of
the class determines the dimensionality of the integration volume(s) and the
name(s) of the configuration parameter(s) to be read from the CosmoSIS
initialization file, but the user of CosmoSIS determines the number of
integration volumes to be used in the run of CosmoSIS, and their sizes.

The class must provide a non-static member function set_grid_point that is
called to allow the setting of any model parameters that define the grid of
points at which the integral is calculated for each sample. This function is
passed an object of type grid_point_t, carrying the data that specifies the
location of the grid point at which the integrand will next be evaluated. The
CosmoSISScalarIntegrationModuleclass template implements the loop over the
grid points specified in the CosmoSIS user’s configuration file.

###  4.5 Using the module

We rely upon several of the integration routines provided by the C-language
library CUBA222http://www.feynarts.de/cuba/(Hahn, 2007). Primarily, we use the
(deterministic) algorithm CUHRE, described in (Berntsen et al., 1991a); we
also make some use of the Monte Carlo algorithm VEGAS (Lepage, 1978). The CUBA
library requires an inconvenient form for the functions it integrates, and
also integrates only over the unit hypercube. It is left to the CUBA user to
provide the necessary transformation of the integrand to handle the
integration volume of interest to the user. To solve these issues, we use the
cubacpp333https://bitbucket.org/mpaterno/cubacpp C++ library. This allows the
user-supplied integrand to be any “callable object” that takes any number of
double arguments, and returns double. It also allows the user to supply an
arbitrary volume of integration, and assures at compilation time that the
volume of integration is of the same dimensionality as the integrand. The user
of the algorithm can specify both _relative_ and _absolute_ tolerances, and
the integration algorithm stops when either condition is met (or when a
specified maximum number of integrand evaluations has been done, providing a
coarse means of limiting the time allowed to calculate an integral). The value
returned from an integration problem includes the estimate of the integral,
and an approximate error in that estimate, as well as some diagnostic
information, most importantly whether the integration algorithm has converged
or not.

A physics module created using the class template described above,
instantiated using a physicist-supplied integrand class, has several
configuration parameters that are to be set by the user of CosmoSIS. The
author of the integrand class has specified the mathematical models to be
used, the dimensionality of the integrand in question, and the parameters that
are used to define the locations in parameter space where the integrals are to
be calculated. The user of CosmoSIS, not the author of the integrand class,
chooses how many grid locations are to be used, and where they are. The user
of CosmoSIS also defines the set of volumes over which the integral should be
evaluated. The user also chooses which of the available integration routines
(CUHRE or VEGAS) should be used, and configures the relative and absolute
error tolerances for the integration. The user also specifies a maximum number
of function evaluations to be used for each integration task, as a means of
limiting the running time of the module for cases when the algorithm chosen is
not converging rapidly enough.

##  5 Tests and Physics Validation

To ensure a robust pipeline, we utilize both unit tests and physics-based
validations. In this set-up, a unit test is designed to test the models or
infrastructural functions (e.g., 1-d and 2-d interpolators), while in cases
that a crucial calculation makes use of multiple models and/or infrastructural
functions, we rely on a physics validation document to check the quality of
the calculations.

###  5.1 Unit Tests

⬇

{

// read in comparison data file

std::ifstream infile =

std::vector<double> zs;

std::vector<double> ys;

while (infile) {

double z, y;

infile >> z >> y;

zs.push_back(z);

ys.push_back(y);

}

OMEGA_Z_SDSS omega;

//Perform the unit test comparison

for (std::size_t i = 0,

sz = zs.size(); i != sz; ++i) {

double const fz = omega(zs[i]);

double constexpr epsrel = 1.0e-6;

CHECK(fz ==

Approx(ys[i]).epsilon(epsrel));

}

}

Figure 4: Abridged Unit test code for OMEGA_Z_SDSS, which computes the
cosmological survey area for the redMaPPer cluster catalog from the Sloan
Digital Sky Survey (SDSS).

We use unit tests to ensure the accuracies of the models in our pipeline.
These unit tests are implemented and ran through the CTest tool
444https://cmake.org/cmake/help/latest/manual/ctest.1.html. An example of the
unit test is shown in Figure 4.

The structures of the unit tests are organized as the follows. A comparison
data set is read in and the relevant function is called to execute at the
comparison data points. The output from the function is compared to the values
in the comparison data sets. If the two values are comparable within an
acceptance limit for every comparison point, the test passes. Otherwise, the
test fails.

To set up the unit tests for the physical models, we acquire expected values
of those model functions at a given set of input values from previously
implemented code in other analysis, e.g., Fortran90 code implemented for SDSS
cluster cosmology analysis (Costanzi et al., 2019b), or Python code
implemented for quantifying cluster mis-centering effects (Zhang et al.,
2019). We require the model output and the comparison value sets to agree
either within a relative or an absolute (in the case of values being zeros)
accuracy limit. The relative accuracy limit, meaning the allowed relative
differences between the values, are by default set at $10^{-6}$ while the
absolute accuracy limit, which is the absolute differences between the values,
is set at $10^{-12}$ by default. In the cases that the comparison values are
calculated through approximations or interpolations, we relax the accuracy
requirements on a case-by-case basis.

These unit tests are designed to be re-ran with a ctest command when the code
re-compiles because of updates or revision and give a pass or fail indication.
Output from some of the unit tests, especially those implementing physics
models, are saved and can be used for debugging or visualization later.

###  5.2 Physics Validation

Figure 5: Excerpt of the physics validation output implemented in R. The
example shows the predictions from our pipeline for the DES Y1 cluster number
counts, compared to the model output in (Abbott et al., 2020) using the same
set of values for cosmological and nuisance parameters. The two pipelines give
highly consistent predictions. The plots and values of relative and absolute
differences are automatically compiled from CosmoSIS output.

We also compile a physics validation document, which visualizes the unit test
results and in addition, visualizes a selected set of more complicated
calculations from the physics modules. Specifically, we compare the module
outputs that calculate the theoretical predictions for a set of SDSS/DES
cluster observables, to the theoretical outputs calculated by previous
analyses (Costanzi et al., 2019b; Abbott et al., 2020).

In designing the physics validation document, our original goal was to
implement the validation in a system that enables version control to easily
track the comparisons. We have considered a few different systems including
static pages such as word or PDF. Jupyter notebooks are also considered, but
ruled out because they are saved in XML which makes it difficult to interpret
the version differences. In the end, the validation is implemented as R
markdown notebooks, which makes it easier to track versions and present the
differences.

The physics validation R notebooks can be compiled at a single command. The
generated document will automatically use the outputs from unit tests and
module outputs from the latest runs, producing a file that visualizes
additional checks of the system. Thus, the validation document provide a more
thorough test of the pipeline. We require agreement with a relative accuracy
of $10^{-6}$ and absolute accuracy of $10^{-12}$ by default in the validation
document, which maybe relaxed if necessary. In Fig. 5, we show an excerpt of
the document comparing predictions for cluster observable in the first year’s
of DES observation data (DES Y1), as an example of what to expect in this
document. The output from our pipeline and the fiducial DES Y1 analysis
pipeline, gives highly consistent predictions for the same set of cosmological
and nuisance parameters.

###  5.3 Instances of unit tests and physics validation

We stress that both unit tests and physics validation are important to ensure
the successful development and configuration of the pipeline. For example, the
unit test of the 2-d interpolator caught an array index switch; finding this
in a physics validation plot would be difficult. On the other hand, physics
validation plots have caught errors that we were not able to uncover in unit
tests. We illustrate this aspect with our experience of catching a mis-match
in the installation of the pipeline.

  1. 1.

Our analyses need to account for the effect of massive neutrinos in the galaxy
cluster mass function, thus we developed modifications to the power spectrum
implementation in the CAMB module in CosmoSIS which only uses the dark and
baryonic matter components – excluding the neutrino component – of the power
spectrum to compute the halo mass function. This is suggested in Costanzi et
al. (2013) as massive neutrinos do not cluster at small scales. Because the
modifications are carried out by a CosmoSIS module that is not part of the
package, we attempted to establish a branch in the CosmoSIS repository that
contains the changes to be used in our development.

  2. 2.

The development of our package continued in the CosmoSIS branch. A few months
later, we onboarded a new member to the team to use and participate in the
development. The new member installed the software and ran the initial tests
and validations.

  3. 3.

All unit tests pass with no issues.

  4. 4.

Upon running the physics-validation files in R, the predictions for the DES Y1
data vectors do not match comparison values. During the investigation, it was
also discovered that the halo mass function output from the CosmoSIS pipeline
do not match previous values.

  5. 5.

We realized that our earlier modifications to the CAMB (Lewis et al., 2000;
Howlett et al., 2012) module 44footnotetext: https://camb.info44footnotetext:
https://bitbucket.org/joezuntz/cosmosis/wiki/default_modules/camb_Jan15 a few
months earlier have not been committed to the CosmoSIS repository. We promptly
committed the branch and requested the CosmoSIS development team for approval,
which since then has been merged into CosmoSIS.

In this case, the success of our pipeline depends on successfully capturing
the dependences which are also undergoing constant developments and updates.
Given that our modification of the CAMB module in CosmoSIS is not part of this
package, those differences were not caught by unit tests. The separate
validation testing which captures individual physics tests makes it possible
to trace and debug this dependence.

##  6 Analysis Cases

###  6.1 DES Y1 cluster Cosmology “Mock” Data

⬇

y1_analysis_mor::operator()(double lo, double zt, double lnM) const

{

std::vector<double> results(2 * zo_low_.size());

double mass = std::exp(lnM);

double common_term = (*mor)(lo, lnM, zt) *

(*dv_do_dz)(zt) * (*hmf)(lnM, zt) * (*omega_z)(zt);

// Number counts first in the returned results,

// followed by the masses

for (std::size_t i = 0; i != zo_low_.size(); i++) {

results[i] = (*int_zo_zt)(zo_low_[i], zo_high_[i], zt) * common_term;

results[i + zo_low_.size()] = mass * results[i];

}

return results;

}

Figure 6: Snippet of the module operator function implementing the integrand
in Equation 1. This operator function is part of the integration class,
described in Section 4, implementing the DES Year 1 analysis. It makes use of
the models described in Section 4 as well: for example ( omega_z)() is a call
to the operator that implements the footprint function, similar to the one
shown in Figure 2.  Figure 7: Comparison between the theoretical predictions
from this pipeline, and the reference pipeline used in DESY1 fiducial
analyses. The upper panels show the predictions (left figure) of cluster
number counts in three redshift ranges (see legend) and four richness ranges
(x-axis), and the relative difference between the two pipelines (right
figure). The lower panels show the predictions of clusters’ average masses in
the redshift/richness ranges (left figure) and the relative differences
between the two pipelines (right figure). Note that the difference comparisons
are highly consistent in the right figures, that the lines representing
different redshift bins overlap. Figure 8: $\Omega_{\mathrm{m}}-\sigma_{8}$
posteriors derived from DES Y1 “Mock” data. The “Mock”data are generated with
$\Omega_{\mathrm{m}}=0.3406$, and $\sigma_{8}=0.8049$. Our pipeline
successfully recovers the truth parameters with the 1$\sigma$ level.

In this section, we demonstrate the efficacy of this pipeline by implementing
a DES Year 1 cluster abundance cosmology analysis Abbott et al. (2020) and
compute prediction data vectors from the pipeline. In this demonstration, we
model the number of clusters, and their average masses in richness/redshift
subsamples.

Specifically, we use the new pipeline to derive theoretical predictions of the
cluster counts and average masses described by Equation 1. As discussed in
Section 4.3 and 4.4, we use the class template to implement the integrals of
Equation 1 in CosmoSIS. Figure 6 shows a code snippet implementing the
integrands in those equations, which is a member function of the physics
module class that will be used to calculate the integrals. Notably, we did not
adopt a projection model (Costanzi et al., 2019a) in the richness–mass
relation in this analysis. The projection model implementation does not
significantly increase the computing time for a DES Y1-like analysis (Equation
1) as discussed in the next subsection, which is dominated by the running time
of other cosmological modules outside of this development such as CAMB, that
can take several seconds to run. Conversely, our modules dominate the running
time when calculating cluster shear profiles in Equation 5, which are much
slower to calculate and can takes more than 100s of seconds using CPU
machines; in this case, skipping the projection model greatly improves the
running speed.

We compare the predictions of number counts and masses computed from our
pipeline and the DES fiducial pipeline using one same set of cosmology and
richness-mass relation parameters. Those comparisons are done five times with
five sets of parameters, and the predictions are all highly consistent. Figure
7 shows the prediction data points for one set of the parameters (with
$\Omega_{\mathrm{m}}=0.3406$, $\sigma_{8}=0.8049$): the number count
predictions from the two pipelines agree within 3% in the lowest richness bin,
or within 1% in the highest richness bin, and the mass predictions agree
within 1%. The prediction differences between the two pipelines are negligible
compared to the measurement uncertainties.

To test the capability of our pipeline to recover unbiased parameter
posteriors, we further use the prediction data vector as a “Mock”
observational data vector (from the DES fiducial Year 1 pipeline), and derive
the cosmological constraints for this “Mock” data set. Because the “Mock”
observations are derived with the same theoretical modelling code in this
pipeline, we expect to recover the cosmological and richness-mass relation
parameters that was used to derive those “Mock” data. To derive the posterior
cosmological constraints, we set up the likelihood calculation: the integral
output from the modules, which are the predictions for the cluster counts and
average masses, are passed to a separate CosmoSIS likelihood module, as
described by equation 11 in Section 3. The likelihood is assumed to be a
multi-dimensional Gaussian distribution with a covariance matrix (as described
in Abbott et al. (2020)).

The posterior parameter distributions are sampled with the “polychord” sampler
implemented in CosmoSIS, as recommended in (Lemos et al., 2022, private
communications with the DES TCP Working Group) which is more accurate compared
to the EMCEE sampler (Foreman-Mackey et al., 2013) used in Abbott et al.
(2020) at sampling a large parameter space with the same number of samples. We
vary and sample a total of 12 parameters (same with Abbott et al. (2020) but
setting both $\Omega_{b}$ and $\Omega_{\nu}$ to be 0), including cosmological
parameters, cluster scaling relation parameters and nuisance parameters.
Finally, Figure 8 shows the posterior constraints on the $\Omega_{m}$ and
$\sigma_{8}$ parameters. For all of the 12 cosmological and richness-masss
relation parameters, our pipeline has successfully recovered the values used
to create the “Mock” observations, within the $1\sigma$ level.

###  6.2 DES Y1 cluster cosmology constraints

Figure 9: $\Omega_{\mathrm{m}}-\sigma_{8}$ posteriors from the DES Y1 cluster
abundance cosmology analysis reanalyzed with a polychord sampler (red dotted
line) and the result from the pipeline described in this paper (black dashed
line). From DES Y1: $\Omega_{\mathrm{m}}$’s constraint is $0.191\pm 0.043$,
while $\sigma_{8}$’s constraint is $0.834\pm 0.059$. From this work:
$\Omega_{\mathrm{m}}$’s constraint is $0.210\pm 0.046$, while $S_{8}$’s
constraint is $0.813\pm 0.057$. The results from the two pipelines are
consistent within $0.5\sigma_{8}$ for both $\Omega_{\mathrm{m}}$ and
$\sigma_{8}$.

Based on the success with the DES Y1 “Mock” observations, we further use the
pipeline to reanalyze the DES Year 1 data vector. We use the same prediction
modules described in the previous subsection, but this time, using the real
DES Year 1 data vector in the likelihood instead of the “Mock” observations.
We sample a total f 14 cosmological and richness-mass parameters as presented
in Abbott et al. (2020), using the same prior distribution for those
parameters.

Figure 9 shows the posterior constraints on the $\Omega_{m}$ and $\sigma_{8}$
parameters from our re-analysis. For comparison, we also show the $\Omega_{m}$
and $\sigma_{8}$ constraints from the fiducial DES analysis pipeline of the
same data vectors but using the polychord sampler. The results are consistent
within 0.5$\sigma$ of the fidual DES analysis despite the differences in the
software pipelines. We note that we have tried different prior setups for the
parameters, and found that in one case, the differences can be further reduced
to within 0.3$\sigma$, but the posteriors presented here are based on the
closest implementation of the DES fiducial analysis (excluding projection
effect in the richness-mass relation). Because our pipeline is an independent
implementation from the DES fiducial analysis, using different programming
languages (C++ VS Fortran), integration algorithms, and pipeline design, the
consistency in the results between the two pipelines also highlight the
robustness of the cosmological constraints from the DES fiducial analysis.

The CosmoSIS module we implemented for this re-analysis, which makes
theoretical predictions of the cluster counts and average masses in a variety
of integration volumes using Equation 1 is fast. Using a single processor on a
Macbook pro (2GHz intel core i5), the calculation for all of the 12 richness
and redshift bins (corresponding integration volumes in our pipeline) listed
in the DES Y1 analysis takes less than 0.1 second for one set of cosmological
and nuisance parameters. Including the projection effect model increases the
computing time to about 1.9 seconds. Both impelmentations are significantly
faster than the other CosmoSIS modules needed for the DES Y1 analysis, such as
the CosmoSIS CAMB module (which can take $\sim$ 4 seconds). Notably, the speed
of those computations change with the values of cosmological and nuisance
parameters, as the integrations in Equation 1 may be slow to converge if the
parameters are far away from a likely range (but necessary to sample in the
MCMC process).

###  6.3 DES Y1 shear profile modeling

Figure 10: Predictions from our pipeline for cluster number counts (upper
panel) and their weak lensing measurement signals (lower panel) in multiple
redshift (each subplot) and richness ranges (each data point/line) for one set
of cosmological and nuisance parameters (not at “best-fit”). In cosmological
analysis, those predictions will be compared to constrain cosmological
parameters. We defer that analysis to a future study.

The cluster masses that we used as an observable in DES Y1 cluster abundance
analysis are further derived from the cluster-weak lensing shear measurements
McClintock et al. (2019b). Given an efficient method, we can directly model
cosmology from cluster number counts and cluster weak-lensing measurements,
bypassing the step to derive masses. Here, we further demonstrate that the
pipeline can be used to compute the predictions of the DES Y1 weak-lensing
observables based on Equations 5 to 9. Using a set of cosmological and
nuisance parameters from the posterior distribution in the previous section
(but not at “best-fit”), and assuming a concentration parameter, we show those
predictions in Figure 10.

Ideally, we will also constrain cosmological parameters through modeling the
measurements, and sampling using MCMC. However, a fully accurate pipeline will
require a good understanding of the selection effects (Abbott et al., 2020) as
well as a robust procedure of testing on simulations (DeRose et al., 2021; To
et al., 2021b) and a blinding scheme (Muir et al., 2020). We defer this
analysis to an ongoing study. In addition, the analysis pipeline developed for
this purpose is also being implemented for GPU computing clusters to speed up
the analysis. We will discuss those implementations in other publications.

##  7 Discussion

In this paper, we describe a software development effort to deliver
theoretical modeling of galaxy cluster cosmology observables implemented in
CosmoSIS for cosmological likelihood calculations. Given the set of
astrophysical and cosmological factors to be taken into account and thus the
challenging computing demands, this sofware package is developed to optimize
speed and accuracy. Specifically, we address several challenges in the
development efforts.

First, the prediction of the cluster cosmological observables depend on a few
models for different physics elements. Each of the models, for example the
cluster mass-observable relation, halo mass function etc., is still under
active development, and subject to active updates. In addition, it is often
desirable to explore several different forms of those models to assess the
robustness of the analysis due to different physical assumptions. Two examples
are the models for the mass-observable relations and the halo mass function.
To make it easy to change a model implementation, we make the models modular,
that one model implementation can be replaced by another that shares use of
the same set of physical parameters. Each model is represented as a C++ class,
which is _callable_ by implementing the function call operator.

Second, the predictions can be implemented as numerical integrals as shown in
Section 2. The dimensions of the integrals are often larger than two, and in
our case, reaches as many as eight. Common integration tools such as those
implemented in Python, or GSL, have only implemented algorithms to handle 1-D
or 2-D integrals. The dimension of those integrals can become even larger for
future analyses if marginalizing over parameters representing additional
systematic effects. We make use of the algorithms described in Berntsen et al.
(1991b), and implemented in the CUBA library Hahn (2005) in our development.

Thirdly, the integrals mentioned above need to be evaluated for each sample of
cosmological and nuisance parameters. If using MCMC sampling, those integrals
often need to be iterated $\sim 10^{6}$ times to evaluate the posterior
distribution of cosmological parameters. In order to finish the sampling in a
reasonable timescale (ideally in a matter of days), the integrals need to be
performed very quickly. We implement all our models and our integration
routines in C++: it provides the best combination of expressiveness and
performance available. We also implement procedures to ensure the
reproducibility of the code, achieved through testing.

In this paper, we demonstrate the effectiveness of this software package by
reproducing the cluster cosmology analysis based on the public DES Year 1 data
sets of cluster number counts and weak-lensing masses. This analysis has also
served as an independent check of the DES Y1 fiducial analysis, verifying the
robustness of the previous analysis. In the future, we expect to further adapt
the software for GPU resources (Sakiotis et al., 2021). Using this package, we
will further extend the DES Year 1 analysis by forward modeling cluster number
counts and weak lensing radial signals. The analysis is currently under-way
and will be described in a future publication.

##  8 Acknowledgement

MESP is funded by the Deutsche Forschungsgemeinschaft (DFG, German Research
Foundation) under Germany’s Excellence Strategy – EXC 2121 ”Quantum Universe”
– 390833306.

Funding for the DES Projects has been provided by the U.S. Department of
Energy, the U.S. National Science Foundation, the Ministry of Science and
Education of Spain, the Science and Technology Facilities Council of the
United Kingdom, the Higher Education Funding Council for England, the National
Center for Supercomputing Applications at the University of Illinois at
Urbana-Champaign, the Kavli Institute of Cosmological Physics at the
University of Chicago, the Center for Cosmology and Astro-Particle Physics at
the Ohio State University, the Mitchell Institute for Fundamental Physics and
Astronomy at Texas A&M University, Financiadora de Estudos e Projetos,
Fundação Carlos Chagas Filho de Amparo à Pesquisa do Estado do Rio de Janeiro,
Conselho Nacional de Desenvolvimento Científico e Tecnológico and the
Ministério da Ciência, Tecnologia e Inovação, the Deutsche
Forschungsgemeinschaft and the Collaborating Institutions in the Dark Energy
Survey.

The Collaborating Institutions are Argonne National Laboratory, the University
of California at Santa Cruz, the University of Cambridge, Centro de
Investigaciones Energéticas, Medioambientales y Tecnológicas-Madrid, the
University of Chicago, University College London, the DES-Brazil Consortium,
the University of Edinburgh, the Eidgenössische Technische Hochschule (ETH)
Zürich, Fermi National Accelerator Laboratory, the University of Illinois at
Urbana-Champaign, the Institut de Ciències de l’Espai (IEEC/CSIC), the
Institut de Física d’Altes Energies, Lawrence Berkeley National Laboratory,
the Ludwig-Maximilians Universität München and the associated Excellence
Cluster Universe, the University of Michigan, NSF’s NOIRLab, the University of
Nottingham, The Ohio State University, the University of Pennsylvania, the
University of Portsmouth, SLAC National Accelerator Laboratory, Stanford
University, the University of Sussex, Texas A&M University, and the OzDES
Membership Consortium.

Based in part on observations at Cerro Tololo Inter-American Observatory at
NSF’s NOIRLab (NOIRLab Prop. ID 2012B-0001; PI: J. Frieman), which is managed
by the Association of Universities for Research in Astronomy (AURA) under a
cooperative agreement with the National Science Foundation.

The DES data management system is supported by the National Science Foundation
under Grant Numbers AST-1138766 and AST-1536171. The DES participants from
Spanish institutions are partially supported by MICINN under grants
ESP2017-89838, PGC2018-094773, PGC2018-102021, SEV-2016-0588, SEV-2016-0597,
and MDM-2015-0509, some of which include ERDF funds from the European Union.
IFAE is partially funded by the CERCA program of the Generalitat de Catalunya.
Research leading to these results has received funding from the European
Research Council under the European Union’s Seventh Framework Program
(FP7/2007-2013) including ERC grant agreements 240672, 291329, and 306478. We
acknowledge support from the Brazilian Instituto Nacional de Ciência e
Tecnologia (INCT) do e-Universo (CNPq grant 465376/2014-2).

This manuscript has been authored by Fermi Research Alliance, LLC under
Contract No. DE-AC02-07CH11359 with the U.S. Department of Energy, Office of
Science, Office of High Energy Physics.

## References

  * Abbott et al. (2020) Abbott, T. M. C., Aguena, M., Alarcon, A., et al. 2020, Phys. Rev. D, 102, 023509, doi: 10.1103/PhysRevD.102.023509
  * Aguena et al. (2021) Aguena, M., Benoist, C., da Costa, L. N., et al. 2021, MNRAS, 502, 4435, doi: 10.1093/mnras/stab264
  * Allen et al. (2011) Allen, S. W., Evrard, A. E., & Mantz, A. B. 2011, Annual Review of Astronomy and Astrophysics, 49, 409, doi: 10.1146/annurev-astro-081710-102514
  * Berntsen et al. (1991a) Berntsen, J., Espelid, T. O., & Genz, A. 1991a, ACM Trans. Math. Softw., 17, 437–451, doi: 10.1145/210232.210233
  * Berntsen et al. (1991b) —. 1991b, ACM Trans. Math. Softw., 17, 437
  * Bocquet et al. (2016) Bocquet, S., Saro, A., Dolag, K., & Mohr, J. J. 2016, MNRAS, 456, 2361, doi: 10.1093/mnras/stv2657
  * Costanzi et al. (2013) Costanzi, M., Villaescusa-Navarro, F., Viel, M., et al. 2013, J. Cosmology Astropart. Phys, 2013, 012, doi: 10.1088/1475-7516/2013/12/012
  * Costanzi et al. (2019a) Costanzi, M., Rozo, E., Rykoff, E. S., et al. 2019a, MNRAS, 482, 490, doi: 10.1093/mnras/sty2665
  * Costanzi et al. (2019b) Costanzi, M., Rozo, E., Simet, M., et al. 2019b, MNRAS, 488, 4779, doi: 10.1093/mnras/stz1949
  * DeRose et al. (2021) DeRose, J., Wechsler, R. H., Becker, M. R., et al. 2021, arXiv e-prints, arXiv:2105.13547.  https://arxiv.org/abs/2105.13547
  * Foreman-Mackey et al. (2013) Foreman-Mackey, D., Conley, A., Meierjurgen Farr, W., et al. 2013, emcee: The MCMC Hammer.  http://ascl.net/1303.002
  * Frieman et al. (2008) Frieman, J. A., Turner, M. S., & Huterer, D. 2008, ARA&A, 46, 385, doi: 10.1146/annurev.astro.46.060407.145243
  * Fumagalli et al. (2021) Fumagalli, A., Saro, A., Borgani, S., et al. 2021, A&A, 652, A21, doi: 10.1051/0004-6361/202140592
  * Gruen et al. (2015) Gruen, D., Seitz, S., Becker, M. R., Friedrich, O., & Mana, A. 2015, MNRAS, 449, 4264, doi: 10.1093/mnras/stv532
  * Hahn (2005) Hahn, T. 2005, Computer Physics Communications, 168, 78, doi: https://doi.org/10.1016/j.cpc.2005.01.010
  * Hahn (2007) Hahn, T. 2007, Computer Physics Communications, 176, 712, doi: 10.1016/j.cpc.2007.03.006
  * Howlett et al. (2012) Howlett, C., Lewis, A., Hall, A., & Challinor, A. 2012, J. Cosmology Astropart. Phys, 2012, 027, doi: 10.1088/1475-7516/2012/04/027
  * Johnston et al. (2007) Johnston, D. E., Sheldon, E. S., Tasitsiomi, A., et al. 2007, ApJ, 656, 27, doi: 10.1086/510060
  * Kravtsov & Borgani (2012) Kravtsov, A. V., & Borgani, S. 2012, ARA&A, 50, 353, doi: 10.1146/annurev-astro-081811-125502
  * Lemos et al. (2022) Lemos, P., Weaverdyck, N., Rollins, R. P., et al. 2022, arXiv e-prints, arXiv:2202.08233.  https://arxiv.org/abs/2202.08233
  * Lepage (1978) Lepage, G. P. 1978, Journal of Computational Physics, 27, 192, doi: 10.1016/0021-9991(78)90004-9
  * Lewis et al. (2000) Lewis, A., Challinor, A., & Lasenby, A. 2000, The Astrophysical Journal, 538, 473, doi: 10.1086/309179
  * Mao et al. (2015) Mao, Y.-Y., Williamson, M., & Wechsler, R. H. 2015, ApJ, 810, 21, doi: 10.1088/0004-637X/810/1/21
  * McClintock et al. (2019a) McClintock, T., Rozo, E., Becker, M. R., et al. 2019a, ApJ, 872, 53, doi: 10.3847/1538-4357/aaf568
  * McClintock et al. (2019b) McClintock, T., Varga, T. N., Gruen, D., et al. 2019b, MNRAS, 482, 1352, doi: 10.1093/mnras/sty2711
  * Muir et al. (2020) Muir, J., Bernstein, G. M., Huterer, D., et al. 2020, MNRAS, 494, 4454, doi: 10.1093/mnras/staa965
  * Palmese et al. (2020) Palmese, A., Annis, J., Burgad, J., et al. 2020, MNRAS, 493, 4591, doi: 10.1093/mnras/staa526
  * Press & Schechter (1974) Press, W. H., & Schechter, P. 1974, ApJ, 187, 425, doi: 10.1086/152650
  * Rozo et al. (2010) Rozo, E., Wechsler, R. H., Rykoff, E. S., et al. 2010, ApJ, 708, 645, doi: 10.1088/0004-637X/708/1/645
  * Rykoff et al. (2012) Rykoff, E. S., Koester, B. P., Rozo, E., et al. 2012, ApJ, 746, 178, doi: 10.1088/0004-637X/746/2/178
  * Rykoff et al. (2014) Rykoff, E. S., Rozo, E., Busha, M. T., et al. 2014, ApJ, 785, 104, doi: 10.1088/0004-637X/785/2/104
  * Rykoff et al. (2016) Rykoff, E. S., Rozo, E., Hollowood, D., et al. 2016, ApJS, 224, 1, doi: 10.3847/0067-0049/224/1/1
  * Sakiotis et al. (2021) Sakiotis, I., Arumugam, K., Paterno, M. F., et al. 2021, Proceedings of the International Conference for High Performance Computing, Networking, Storage and Analysis
  * The LSST Dark Energy Science Collaboration et al. (2018) The LSST Dark Energy Science Collaboration, Mandelbaum, R., Eifler, T., et al. 2018, arXiv e-prints, arXiv:1809.01669.  https://arxiv.org/abs/1809.01669
  * Tinker et al. (2008) Tinker, J., Kravtsov, A. V., Klypin, A., et al. 2008, ApJ, 688, 709, doi: 10.1086/591439
  * To et al. (2021a) To, C., Krause, E., Rozo, E., et al. 2021a, Phys. Rev. Lett., 126, 141301, doi: 10.1103/PhysRevLett.126.141301
  * To et al. (2021b) To, C.-H., Krause, E., Rozo, E., et al. 2021b, MNRAS, 502, 4093, doi: 10.1093/mnras/stab239
  * Umetsu (2020) Umetsu, K. 2020, A&A Rev., 28, 7, doi: 10.1007/s00159-020-00129-w
  * Varga et al. (2019) Varga, T. N., DeRose, J., Gruen, D., et al. 2019, MNRAS, 489, 2511, doi: 10.1093/mnras/stz2185
  * Weinberg et al. (2013) Weinberg, D. H., Mortonson, M. J., Eisenstein, D. J., et al. 2013, Phys. Rep., 530, 87, doi: 10.1016/j.physrep.2013.05.001
  * Wu et al. (2019) Wu, H.-Y., Weinberg, D. H., Salcedo, A. N., Wibking, B. D., & Zu, Y. 2019, MNRAS, 490, 2606, doi: 10.1093/mnras/stz2617
  * Zhang & Annis (2022) Zhang, Y., & Annis, J. 2022, arXiv e-prints, arXiv:2201.02167.  https://arxiv.org/abs/2201.02167
  * Zhang et al. (2019) Zhang, Y., Jeltema, T., Hollowood, D. L., et al. 2019, MNRAS, 487, 2578, doi: 10.1093/mnras/stz1361
  * Zuntz et al. (2015) Zuntz, J., Paterno, M., Jennings, E., et al. 2015, Astronomy and Computing, 12, 45, doi: 10.1016/j.ascom.2015.05.005

##  Appendix A Institutions

1 Laboratório Interinstitucional de e-Astronomia - LIneA, Rua Gal. José
Cristino 77, Rio de Janeiro, RJ - 20921-400, Brazil
2 Department of Physics, University of Michigan, Ann Arbor, MI 48109, USA
3 Fermi National Accelerator Laboratory, P. O. Box 500, Batavia, IL 60510, USA
4 Institute of Cosmology and Gravitation, University of Portsmouth,
Portsmouth, PO1 3FX, UK
5 University Observatory, Faculty of Physics, Ludwig-Maximilians-Universität,
Scheinerstr. 1, 81679 Munich, Germany
6 Department of Physics & Astronomy, University College London, Gower Street,
London, WC1E 6BT, UK
7 Universidad de La Laguna, Dpto. Astrofisica, E-38206 La Laguna, Tenerife,
Spain
8 Instituto de Astrofísica de Canarias, E-38205 La Laguna, Tenerife, Spain
9 Kavli Institute for Cosmological Physics, University of Chicago, Chicago, IL
60637, USA
10 Department of Astronomy and Astrophysics, University of Chicago, Chicago,
IL 60637, USA
11 Astronomy Unit, Department of Physics, University of Trieste, via Tiepolo
11, I-34131 Trieste, Italy
12 INAF-Osservatorio Astronomico di Trieste, via G. B. Tiepolo 11, I-34143
Trieste, Italy
13 Institute for Fundamental Physics of the Universe, Via Beirut 2, 34014
Trieste, Italy
14 School of Mathematics and Physics, University of Queensland, Brisbane, QLD
4072, Australia
15 Centro de Investigaciones Energéticas, Medioambientales y Tecnológicas
(CIEMAT), Madrid, Spain
16 Jet Propulsion Laboratory, California Institute of Technology, 4800 Oak
Grove Dr., Pasadena, CA 91109, USA
17 Institute of Theoretical Astrophysics, University of Oslo. P.O. Box 1029
Blindern, NO-0315 Oslo, Norway
18 SLAC National Accelerator Laboratory, Menlo Park, CA 94025, USA
19 Center for Astrophysical Surveys, National Center for Supercomputing
Applications, 1205 West Clark St., Urbana, IL 61801, USA
20 Department of Physics and Astronomy, University of Pennsylvania,
Philadelphia, PA 19104, USA
21 Institut de Física d’Altes Energies (IFAE), The Barcelona Institute of
Science and Technology, Campus UAB, 08193 Bellaterra (Barcelona) Spain
22 Department of Astronomy, University of Illinois at Urbana-Champaign, 1002
W. Green Street, Urbana, IL 61801, USA
23 Santa Cruz Institute for Particle Physics, Santa Cruz, CA 95064, USA
24 Center for Cosmology and Astro-Particle Physics, The Ohio State University,
Columbus, OH 43210, USA
25 Department of Physics, The Ohio State University, Columbus, OH 43210, USA
26 Center for Astrophysics $|$ Harvard & Smithsonian, 60 Garden Street,
Cambridge, MA 02138, USA
27 Department of Physics, University of Arizona, Tucson, AZ 85721, USA
28 Lowell Observatory, 1400 Mars Hill Rd, Flagstaff, AZ 86001, USA
29 Australian Astronomical Optics, Macquarie University, North Ryde, NSW 2113,
Australia
30 Department of Physics, Carnegie Mellon University, Pittsburgh, Pennsylvania
15312, USA
31 George P. and Cynthia Woods Mitchell Institute for Fundamental Physics and
Astronomy, and Department of Physics and Astronomy, Texas A&M University,
College Station, TX 77843, USA
32 Jodrell Bank Center for Astrophysics, School of Physics and Astronomy,
University of Manchester, Oxford Road, Manchester, M13 9PL, UK
33 Institució Catalana de Recerca i Estudis Avançats, E-08010 Barcelona, Spain
34 Hamburger Sternwarte, Universität Hamburg, Gojenbergsweg 112, 21029
Hamburg, Germany
35 Observatório Nacional, Rua Gal. José Cristino 77, Rio de Janeiro, RJ -
20921-400, Brazil
36 Kavli Institute for Particle Astrophysics & Cosmology, P. O. Box 2450,
Stanford University, Stanford, CA 94305, USA
37 Department of Physics and Astronomy, Pevensey Building, University of
Sussex, Brighton, BN1 9QH, UK
38 Department of Physics and Astronomy, Stony Brook University, Stony Brook,
NY 11794, USA
39 School of Physics and Astronomy, University of Southampton, Southampton,
SO17 1BJ, UK
40 Computer Science and Mathematics Division, Oak Ridge National Laboratory,
Oak Ridge, TN 37831
41
42 Universitäts-Sternwarte, Fakultät für Physik, Ludwig-Maximilians
Universität München, Scheinerstr. 1, 81679 München, Germany
43 Max Planck Institute for Extraterrestrial Physics, Giessenbachstrasse,
85748 Garching, Germany
44 Department of Physics, Boise State University, Boise, ID 83725, USA
45 NSF’s National Optical-Infrared Astronomy Research Laboratory, 950 N.
Cherry Ave., Tucson, AZ, 85719, USA
46 Department of Physics, Duke University Durham, NC 27708, USA


# A portable and flexible implementation of the Wang–Landau algorithm in order
to determine the Density of States

Felipe Moreno<EMAIL_ADDRESS>Sergio Davis<EMAIL_ADDRESS>Joaquín Peralta<EMAIL_ADDRESS>Departamento de Física, Facultad de
Ciencias Exactas, Universidad Andrés Bello, Sazié 2212, Santiago, Chile.
Comisión Chilena de Energía Nuclear, Casilla 188-D, Santiago 7600713, Chile.

###### Abstract

In this work we develop an implementation of the Wang-Landau algorithm [Phys.
Rev. Lett. 86, 2050-2053 (2001)]. This algorithm allows us to find the density
of states (DOS), a function that, for a given system, describes the proportion
of states that have a certain energy. The implementation uses the Python
language for the algorithm itself, and it can take advantage of any library,
such as the powerful LAMMPS library, for the computation of energy. Therefore,
the resulting implementation is simple and flexible without sacrificing
efficiency. This implementation also considers recent developments in the
parallelization of the code for faster computation. We establish the soundness
and effectiveness of our implementation by studying well-known systems such as
the Ising model, the Lennard-Jones and EAM solids. We have found that our
implementation can find the DOS with very good precision in a reasonable
amount of time. Therefore, we are equipped with a very powerful and flexible
implementation that can be easily used in order to study more realistic models
of matter.

††journal: Computer Physics Communications

##  1 Introduction

The density of states (DOS), denoted by $\Omega(E)$ and defined as
$\Omega(E)=\int d\mathbf{r}\,\delta(E-H(\mathbf{r}))$, is a fundamental
function of a system such that $\Omega(E)dE$ is the number of states that have
an energy between $E$ and $E+dE$. This function can be used to compute all of
the thermodynamic properties of a system. Examples of such properties are the
caloric curve and the heat capacity, both being very important in the study of
any system. Therefore, computing the DOS is of fundamental importance, and
consequently several methods for the calculation of this function for a given
system have been developed [1, 2, 3], some of them based in generalized
ensembles and reweighting techniques.

Among the generalized-ensemble methods, Wang-Landau (WL) is a well known Monte
Carlo technique for computing the DOS of a system [4]. This method was
originally developed for discrete systems and single processors in mind.
However, since then the method has been expanded to off-lattice systems and
parallel processing[5]. This has been done in order to take full advantage of
the powerful multicore processing units that are the standard today. Despite
the widespread use of the algorithm and a large number of variants and
improvements [6, 7, 8], there are still practical difficulties, ranging from
the fine-tuning needed for an optimal operation, to the lack of a universal
implementation able to be “plugged into” traditional simulation codes.

In this work, we present a state-of-the-art implementation of WL using the
flexibility and portability of the Python programming language [9]. This
implementation takes advantage of two well known libraries: The powerful
LAMMPS [10] library, which will be used to compute the energy of a complex
configuration of atoms, and the easy-to-use Python multiprocessing library for
parallel computing [11]. Furthermore, this implementation supports both
discrete and continuous systems. Despite previous efforts to improve the
algorithm such as the optimization of flatness criteria [12], parallelization
scheme by replica-exchange [5], and more efficient Monte Carlo trial moves
[13] among others, to the best of our knowledge there is no standard library
that incorporates a large set of improvements of the algorithm. Here we build
up a well-modularized system that allows the scientific community to make an
efficient use and incorporation of different aspects of the WL algorithm.

The applicability of the WL algorithm is quite broad to different research
topics, such as integration [14], protein folding [15], polymers [16], free
energy profile [17], and spin-glass models [18, 19, 20], among others. Here we
focus the technique mainly on the calculation of the caloric-curve for
different cases, to expand its use towards the condensed matter community. We
tested this implementation using a variety of well-established benchmark
systems such as the Ising model, the Potts model, classical and quantum
harmonic oscillator and the Lennard-Jones system. The goal of these tests was
to prove the soundness of our implementation through the comparison with
previously known results.

This paper is organized as follows. Following this introduction, in Section 2
we describe the inner workings of the Wang–Landau algorithm. In Section 3 we
offer a brief description of this implementation. In Section 4 we provide
details of the benchmark simulations performed. A discussion of our findings
is given in Section 5.

##  2 A brief review of the Wang-Landau algorithm

The WL algorithm is used to find the DOS of a many-particle system with $N$
degrees of freedom. This system can be described by a configurational state
$X=(X_{1},X_{2},\ldots,X_{N})$ and a function of $X$ that describes the system
energy $E=\mathcal{H}(X)$. This function is called the Hamiltonian.

In order to produce the DOS of a system, the Wang-Landau algorithm takes
advantage of the following property. Suppose we sample a large number of
possible states according to a distribution that depends only on the
Hamiltonian of each state $P(X)=\rho(\mathcal{H}(X))$. If this distribution is
chosen to be

$\rho(\mathcal{H}(X))\propto\frac{1}{\Omega(\mathcal{H}(X))},$ (1)

then the energies will be distributed according to

$P(E)=\rho(E)\Omega(E)=\text{constant}.$ (2)

Consequently, we would build up a flat histogram. Thereby, we can use this
property to test how good a candidate DOS is, and even better, we can
determine how to make a small change to it so that it is closer to the true
DOS. We keep on making those small changes to the candidate DOS until the
histogram is flat enough. The family of algorithms that takes advantage of
this property are known as flat histogram methods [21].

The specific steps to reproduce the above procedure for discrete systems is
detailed now. At the start of our simulation we define a candidate DOS of the
system, $\Omega(E)=1$ for all $E$. We define a _walker_ as an entity that
wanders over the configuration space standing on an unique configuration at
every simulation step. The energy of a walker is the energy of the
configuration where this walker is currently standing on. We start out with a
walker on a random state $X$, and then we produce a trial move from $X$ to
$X^{\prime}$. The transition acceptance that is in concordance with Eq. 1 and
the detailed balance condition is

$W(X\to
X^{\prime})=\min\left(1,\frac{\Omega(\mathcal{H}(X))}{\Omega(\mathcal{H}(X^{\prime}))}\right).$
(3)

If the trial move is accepted, the walker moves to the new state $X^{\prime}$,
otherwise, it remains on its previous state $X$. Whether the change was
accepted or not, we update the value of the DOS at the current energy $E$
according to $\Omega(E)\leftarrow\Omega(E)\cdot f$, where $f>1$. Additionally,
we update the histogram that counts up the number of times a certain energy
has been visited. This process continues until the histogram is sufficiently
flat. The flatness criterion is that the minimum value of the histogram must
be at least 80% of its mean value. Once this criterion has been met, the
process is interrupted, the value of $f$ is decreased according to
$f\leftarrow\sqrt{f}$, the histogram is reset, and the same process is started
over again. Finally, the algorithm is stopped when $f$ reaches a value close
enough to unity, such as $f\approx 1+10^{-8}$, and the final $\Omega(E)$ is
considered to be the true DOS.

For continuous systems, the energies are discretized into small bins of width
$\Delta E$. Each bin has its own $\Omega$ and histogram, and they are updated
when the energy of the walker lies between $E-\Delta E/2$ and $E+\Delta E/2$.
The transition probability is computed after previously performing a cubic
interpolation on the discrete $\Omega$, thus obtaining $\Omega(E)$ from the
interpolated function. This last step provides far better results than merely
choosing the discrete $\Omega$ of the bin that contains $E$. This is because
the walker will prefer a transition to a less probable energy, even if this
new energy and the current one belong to the same bin. Furthermore, instead of
the discrete case where the way in which we produce a change is obvious, we
must select an appropriate step size, or alternatively, we can use a variable
step size through the definition of a target acceptance. A target acceptance
of 30% has been used in the continuous systems.

Figure 1: The entire energy range is split up into overlapping windows.

The algorithm is parallelized in the following way. Firstly, the entire energy
range is split up into overlapping windows of energy (Fig. 1). An overlap of
75% has been reported to provide good results [5]. Secondly, every window is
given its own walker and then a so-called “sweep” is executed—the performing
of a certain number of simulation steps inside of every window. These steps
run in parallel and they are completely independent. Thirdly, a replica
exchange procedure [5] is tried. This procedure consists in the exchange of
walkers between two windows. For this exchange to occur, the windows must be
adjacent and their walkers must have an energy that lies over the overlapping
interval of their respective windows. If the previous conditions are met for
windows $i$ and $j$, the exchange is performed according to the probability

$P(w_{i}\leftrightarrow
w_{j})=\min\left(1,\frac{\Omega_{i}(E_{i})\Omega_{j}(E_{j})}{\Omega_{i}(E_{j})\Omega_{j}(E_{i})}\right),$
(4)

where $w$, $E$ and $\Omega$ are the walker, the walker’s current energy and
the current DOS of their respective windows. This procedure allow walkers to
travel through the entire energy landscape thus avoiding the presence of
energy traps. Fourthly, we repeat the last two steps until every window has
ended its own simulation. The windows that have already ended their
simulations are required to continue running them without updating their DOS.
This is required because we need that those windows keep on performing replica
exchanges with the active ones. Finally, we reconstruct the full DOS from the
pieces of DOS obtained at every window.

In order to reconstruct $\Omega$, we must consider that we want to produce the
smoothest curve possible, in particular, we want the microcanonical caloric
curve $\beta(E)=\partial(\log\Omega(E))/\partial E$ to be smooth. Having this
in mind, we select the pieces of $\Omega$ that belong to the first
($\Omega_{1}$) and second ($\Omega_{2}$) windows, and then we paste them
through the following procedure. Firstly, we compute
$\beta_{1}(E)=\partial(\log\Omega_{1}(E))/\partial E$ and
$\beta_{2}(E)=\partial(\log\Omega_{2}(E))/\partial E$ and we select the energy
$E^{*}$ where the best match between these two functions is. Secondly, we
define $\bar{\Omega}_{2}=[\Omega_{1}(E^{*})/\Omega_{2}(E^{*})]\Omega_{2}$ so
that $\Omega_{1}$ and $\bar{\Omega}_{2}$ match their values at $E^{*}$.
Thirdly, we paste $\Omega_{1}$ and $\bar{\Omega}_{2}$ at $E^{*}$, cutting out
the leftover “tails” of each function. Finally, the above steps are repeated
between the previous pasted $\Omega$ and the next piece of DOS until there are
no more pieces left.

##  3 Implementation

Figure 2: Diagram of the main process splitted in different subprocess with a

This implementation is based on several modules. Firstly, we have a module
that loads the simulation data such as the energy range, step size, target
acceptance and so on from an input file. Secondly, there is the module that
runs the actual simulations. This module creates the subprocesses that will
run the WL simulation on every window, perform the replica exchanges, merge
the DOS pieces and write the result into a file (Fig. 2). Finally, there are
the walker modules, one for each system. These walkers keep track of the state
of the system, perform changes, and compute the energy. Furthermore, those
walkers are totally independent of the algorithm itself, so they can be used
for any other algorithm and they can be coded in another language. For
example, a LAMMPS walker has been developed which contains a wrapper for the
C++ LAMMPS library. Any other library can be used through a wrapper class.
Because of this, this implementation does not suffer of slowness due to the
fact of being coded in Python, since for a typical problem the time consuming
function is the walker method that computes the energy of the system, which
can be implemented in a faster language like C++ or Julia. The algorithm
itself adds negligible overhead.

##  4 Simulations and computational procedures

A set of different simulations were performed in this work in order to
determine the accuracy and practical usefulness of our implementation. For
every system, instead of computing the full DOS, we only computed the
configurational density of states (CDOS). This means that instead of using the
full Hamiltonian $\mathcal{H}=K+\phi$, we only used the configurational part
$\phi$. This is done because the contribution of the kinetic part $K$ is
always the same and it is easy to compute analytically.

The caloric curve for both microcanonical and canonical ensembles has been
computed using the CDOS ($\Omega(\phi)$) obtained in the simulation. If $d$ is
the number of degrees of freedom in the system, the microcanonical curve has
been computed according to

$\expectationvalue{\beta}_{E}=\partialderivative{E}\log\eta(E),$ (5)

where

$\eta(E)\equiv\int\Omega(\phi)\left[E-\phi\right]^{\frac{d}{2}-1}d\phi$ (6)

is a function that is proportional to the full DOS [22], and the canonical
curve has been computed according to the well known equation

$\expectationvalue{E}_{\beta}=\frac{d}{2\beta}-\partialderivative{\beta}\log
Z(\beta),$ (7)

where

$Z(\beta)\equiv\int\Omega(\phi)e^{-\beta\phi}d\phi$ (8)

is the configurational part of the partition function, and $d/2\beta$ is the
kinetic energy term.

The results are presented below. The discrete models do not use the LAMMPS
library because they are not atomic models. They are presented just as a proof
of the soundness of our implementation. The fifth one, the Lennard–Jones
solid, is the most important result because it shows how our code works
flawlessly with the LAMMPS library which enables us to study a great quantity
of systems.

###  4.1 Discrete Models

Figure 3: CDOS (left) and caloric curve (right) for the 2-dimensional Ising
model on a $10\times 10$ lattice. Note that the caloric curve unusually
contains the contribution of the kinetic energy of the system, only for
purposes of comparison with atomistic models.

The Ising model is the first discrete model studied, it is a well known model
of interacting spins. The spins can take two values—up and down. This is one
of the most simple model that exhibits a phase transition. The interaction
between the spins is described by

$\phi=-J\sum_{\langle i,j\rangle}\bm{\sigma_{i}}\cdot\bm{\sigma_{j}},$ (9)

where the sum is over adjacent pairs. Usually, there is an additional term
with the contribution of an external magnetic field. We use an 2D Ising Model
of $N=100$ spins in a $10\times 10$ lattice without the magnetic field. The
continuous transition can be observed in the caloric curve computed in the
simulation shown in Fig. 3. The transition temperature is $k_{B}T_{c}/J=2.332$
which is close to the well-known exact transition temperature
$k_{B}T_{c}/J=2.269$ of this model [23].

Figure 4: CDOS (left) and caloric curve (right) for the Potts model with
$q=10$ on a $10\times 10$ lattice. Note that the caloric curve unusually
contains the contribution of the kinetic energy of the system, only for
purposes of comparison with atomistic models.

A most generalized discrete study case is the Potts model [24]. Potts is a
generalized Ising model of interacting classical spins $q$ on a regular
lattice. In this model the spins can take $q$ different “orientations”. This
model is described by

$\phi=-J\sum_{\langle i,j\rangle}\delta(\sigma_{i},\sigma_{j})$ (10)

where the sum is over adjacent pairs and $\delta$ is the usual Kronecker
delta. In the particular case of $q=2$ we recover the 2D-Ising model up to a
multiplicative factor. In what follows, we use a value of $q=10$ on a
$10\times 10$ 2-dimensional lattice. This model exhibits a first-order phase
transition which can be observed looking at the caloric curve obtained in the
simulation (Fig. 4). This model has been widely used in multiple research
areas [25, 26, 27, 19, 20, 21].

We obtained the same result than in our previous work where the phase
transitions and metastable phases of this model have been characterized [28].
The metastable region lies inside the S-shaped loop in the microcanonical
curve.

Figure 5: CDOS (left) and caloric curve (right) for a system of 100 Quantum
Harmonic Oscillators.

The third, and last, discrete model is a system of $N$ identical 1-dimensional
quantum harmonic oscillators described by

$\phi=\frac{1}{2}\hbar\omega\left(\sum_{i=1}^{N}n_{i}+\frac{1}{2}\right),$
(11)

where $n_{i}$ is a non-negative integer that is equal to the quantum number of
the i-th oscillator. We consider $N=100$ oscillators. The plot (Fig. 5) shows
how the caloric curve starts with a very low slope and then this slope
gradually reaches the well-known [29] classical value $E/Nk_{B}T=1$. Despite
the applicability in the particular case of phase transition of discrete
systems, and the wide research during last decades [30, 31, 32, 25, 26, 27,
19, 20, 21], these prevail as validation in a large variety of research.

###  4.2 Continous Systems

Figure 6: CDOS (left) and caloric curve (right) for a system of 100 Classical
Harmonic Oscillators.

As a first example of a continuous system we consider a system of $N$
identical 1-dimensional classical harmonic oscillators described by

$\phi=\frac{1}{2}m\omega^{2}\sum_{i=1}^{N}x_{i}^{2}.$ (12)

where $x_{i}$ describes the position of the i-th oscillator. The energy is
computed in units of $m\omega x_{0}^{2}$ where $x_{0}$ is an arbitrary length.
At every step a random oscillator is displaced a random quantity drawn from an
uniform distribution $[-x_{0},x_{0}]$. As expected, in Fig. 6 we obtain a
straight line as the caloric curve [29], with the correct slope of 1.

Figure 7: CDOS (left) and caloric curve (right) for a Lennard–Jones system
with $N=120$ atoms.

The second continuous system studied is a Lennard–Jones (LJ) solid by using
the energy provided by the classical LJ potential [33]. The LJ potential is
described by

$\phi=4\epsilon\sum_{i.j<i}\left[\left(\frac{\sigma}{r_{ij}}\right)^{12}-\left(\frac{\sigma}{r_{ij}}\right)^{6}\right].$
(13)

where the sum is over every pair of atoms. We utilize a cut-off value of
$2.5\sigma$ which it means that we neglect the contribution of pairs where
$r_{ij}>2.5\sigma$. The structure is made of $N=120$ atoms with masses $m$ in
a cubic cell of side $15\sigma$. This potential and the EAM potential are the
most important model studied since they use the LAMMPS library. For this
model, an initial minimization has been performed in order to find the low
energy solid states, thus determining an appropriate energy range for our
simulation. Looking at the caloric curve (Fig. 7) we observe a transition
temperature of $k_{B}T/\epsilon\approx 0.7$ which agrees with the literature.
This model produces very accurate results for noble gases, for example, for
the Argon solid in which $\epsilon=1.03\times 10^{-2}\text{ eV}$, this model
reports a temperature of melting close to $T=83.59\text{ K}$.

Figure 8: CDOS (left) and caloric curve (right) for an EAM system with $N=432$
atoms.

As a last study case, we use the embedded atom model (EAM) interatomic
potential. The EAM is widely used in metallic systems, and incorporate a local
density on the atomic interaction. A generalized form of this type of
potential es given by

$E_{i}=F_{\alpha}\left(\sum_{i\neq
j}\rho_{\beta}(r_{ij})\right)+\frac{1}{2}\sum_{i\neq
j}\phi_{\alpha\beta}(r_{ij})$ (14)

where $r_{ij}$ is the distance between atoms $i$ and $j$, $\phi_{\alpha\beta}$
is a pair-wise potential function, $\rho_{\beta}$ is the contribution to the
electron charge density from atom $j$ of type $\beta$ at the location of atom
$i$, and $F$ is an embedding function that represents the energy required to
place atom $i$ of type $\alpha$ into the electron cloud. Despite this kind of
potentials have been widely used [34, 35, 36], they are still highly time-
consuming on computational procedures. Testing this results, help us to define
in an accurate way the scalability of our implementation of the WL-algorithm.
The parameters used in the simulations are the ones that fit the energies for
Fe atoms [37]. The results are presented in Fig. 8. To the best of our
knowledge, there are no previous WL simulations on EAM systems.

##  5 Concluding remarks

Our work presents a modular implementation of the original Wang-Landau
algorithm enhanced with state-of-the-art developments such as off-lattice
systems and parallelization. It also integrates seamlessly with the LAMMPS
library for fast computing of the energy of atomic configurations, which means
that we can, in principle, study any system that can be described by a
potential supported by LAMMPS, or, in principle, any other atomistic
simulation package for which a wrapper can be written. The results obtained in
our simulation of well known systems matched what we expected, including the
specially important cases of the Lennard-Jones system and the EAM system
because these are the models that run with the LAMMPS library. Consequently,
our implementation of the WL algorithm can be used from now on to efficiently
study the properties of more realistic models of matter.

## Acknowledgments

SD acknowledges partial support by ANID FONDECYT 1171127 and ANID PIA
ACT172101 grants. FM thanks to Universidad Andrés Bello for providing the
doctorate scholarship. Computational work was supported by the supercomputing
infrastructures of the NLHPC (ECM-02), and FENIX (UNAB).

## References

  * Shimizu [2004] H. Shimizu, Phys. Rev. E 70, 056704 (2004).
  * Habeck [2007] M. Habeck, Phys. Rev. Lett. 98, 200601 (2007).
  * Montecinos et al. [2021] A. Montecinos, C. Loyola, J. Peralta, and S. Davis, Phys. A 562, 125279 (2021).
  * Wang and Landau [2001] F. Wang and D. P. Landau, Phys. Rev. Lett. 86, 2050 (2001).
  * Vogel et al. [2013] T. Vogel, Y. W. Li, T. Wüst, and D. P. Landau, Phys. Rev. Lett. 110, 210603 (2013).
  * Liang [2005] F. Liang, Journal of the American Statistical Association 100, 1311 (2005).
  * Atchadé and Liu [2010] Y. F. Atchadé and J. S. Liu, Statistica Sinica 20, 209 (2010), ISSN 10170405, 19968507\.
  * Bornn et al. [2013] L. Bornn, P. E. Jacob, P. D. Moral, and A. Doucet, Journal of Computational and Graphical Statistics 22, 749 (2013).
  * van Rossum [1995] G. van Rossum, Tech. Rep. CS-R9526, Centrum voor Wiskunde en Informatica (CWI), Amsterdam (1995).
  * Plimpton [1995] S. Plimpton, J. Comput. Phys. 117, 1 (1995).
  * Dalcin et al. [2011] L. D. Dalcin, R. R. Paz, P. A. Kler, and A. Cosimo, Advances in Water Resources 34, 1124 (2011), ISSN 0309-1708, new Computational Methods and Software Tools.
  * Belardinelli and Pereyra [2007] R. E. Belardinelli and V. D. Pereyra, Phys. Rev. E 75, 046701 (2007).
  * Wüst and Landau [2009] T. Wüst and D. P. Landau, Phys. Rev. Lett. 102, 178101 (2009).
  * Tröster and Dellago [2005] A. Tröster and C. Dellago, Phys. Rev. E 71, 066705 (2005).
  * Rathore et al. [2003] N. Rathore, T. A. Knotts, and J. J. de Pablo, The Journal of Chemical Physics 118, 4285 (2003).
  * Taylor et al. [2009] M. P. Taylor, W. Paul, and K. Binder, The Journal of Chemical Physics 131, 114907 (2009).
  * Calvo [2002] F. Calvo, Molecular Physics 100, 3421 (2002).
  * Brown and Schulthess [2005] G. Brown and T. C. Schulthess, Journal of Applied Physics 97, 10E303 (2005).
  * Torbrügge and Schnack [2007] S. Torbrügge and J. Schnack, Phys. Rev. B 75, 054403 (2007).
  * Alder et al. [2004] S. Alder, S. Trebst, A. K. Hartmann, and M. Troyer, Journal of Statistical Mechanics: Theory and Experiment 2004, P07008 (2004).
  * Landau and Binder [2014] D. P. Landau and K. Binder, _A Guide to Monte Carlo Simulations in Statistical Physics_ (Cambridge, 2014).
  * Davis et al. [2020] S. Davis, J. Jain, and B. Bora, Phys. Rev. E 102, 042137 (2020).
  * Baxter [2016] R. J. Baxter, _Exactly solved models in Statistical Mechanics_ (Elsevier, 2016).
  * Wu [1982] F. Y. Wu, Rev. Mod. Phys. 54, 235 (1982).
  * Turner and Sherratt [2002] S. Turner and J. A. Sherratt, J. Theor. Biol. 216, 85 (2002).
  * Davis et al. [2014] S. Davis, Y. Navarrete, and G. Gutiérrez, Eur. Phys. J. B 87, 78 (2014).
  * Farías and Davis [2020] C. Farías and S. Davis, arXiv e-prints arXiv:2008.04368 (2020), 2008.04368.
  * Moreno et al. [2018] F. Moreno, S. Davis, C. Loyola, and J. Peralta, Phys. A 509, 361 (2018).
  * Greiner et al. [2012] W. Greiner, L. Neise, and H. Stöcker, _Thermodynamics and statistical mechanics_ (Springer Science & Business Media, 2012).
  * Carollo et al. [2020] A. Carollo, D. Valenti, and B. Spagnolo, Physics Reports 838, 1 (2020), ISSN 0370-1573, geometry of quantum phase transitions.
  * Dekker [1981] H. Dekker, Physics Reports 80, 1 (1981), ISSN 0370-1573.
  * Um et al. [2002] C.-I. Um, K.-H. Yeon, and T. F. George, Physics Reports 362, 63 (2002), ISSN 0370-1573.
  * Lennard-Jones [1931] J. E. Lennard-Jones, Proceedings of the Physical Society 43, 461 (1931).
  * Daw and Baskes [1984] M. S. Daw and M. I. Baskes, Phys. Rev. B 29, 6443 (1984).
  * Daw et al. [1993] M. S. Daw, S. M. Foiles, and M. I. Baskes, Materials Science Reports 9, 251 (1993), ISSN 0920-2307.
  * Ackland and Bonny [2020] G. J. Ackland and G. Bonny, pp. 544–572 (2020).
  * Malerba et al. [2010] L. Malerba, M. Marinica, N. Anento, C. Björkas, H. Nguyen, C. Domain, F. Djurabekova, P. Olsson, K. Nordlund, A. Serra, et al., Journal of Nuclear Materials 406, 19–38 (2010).