© J. Sweezy et al., Published by EDP Sciences, 2025

1. Introduction

The Monte Carlo Application Toolkit (MCATK) is a C++ component-based Monte Carlo particle transport capability developed by Los Alamos National Laboratory. The goals of MCATK are multifold, including (but not limited to): (1) provide components for general-purpose Monte Carlo codes, such as the MCNP^® [3] software, (2) provide components for developing Monte Carlo tools for specific bespoke applications, (3) provide well-tested and reliable components, and (4) provide a modern code basis for use on advanced computing architectures. This paper provides an update on the recent additions to MCATK’s capabilities, especially those that have been added since the prior update papers in 2015 and 2017 [1, 2].

The development of MCATK began in 2008 is a part of the response to change the computing landscape. In 2008 the Roadrunner [4] supercomputer was installed at LANL, breaking the petaflop barrier. Roadrunner used AMD Opteron CPUs accelerated with the IBM Cell Broadband Engine (also used in the Sony PlayStation (3) [5]. The architecture and its challenges foretold the current state of high performance scientific computing that is now prevalent in the form of GPU accelerated HPC platforms. The architecture presented multiple programming challenges to software developers, including that no Fortran compiler was initially available. While Fortran continues to be well supported for CPU support, programming GPUs through Fortran remains a challenge.

Also at the same time, there came a need to develop bespoke Monte Carlo applications that provided significant performance, reduced time to delivery, with targeted verification and validation. It became clear that meeting all these demands with the MCNP^® software, while also supporting the wide user base, would be difficult.

Fig. 1. MCATK geometry plotter, showing two views of the Tinkertoy 2 critical assembly.

1.1. Development methodology

The early development of MCATK (2008–2013) [1] was performed using Extreme Programming (XP) [6] and Agile [7] software development methodologies. Hallmarks of this development included, paired programming, co-location of developers in the same room, test-driven development, unit testing, frequent commits to the head of the repository, continuous building and testing, a goal of testing and building in less than 10 minutes, and frequent delivery of products to users.

The introduction of the Git version control system [8], a growing development team, and multiple years with the team working from home due to the coronavirus pandemic – created an environment that changed the software development methods used by the MCATK team. Currently, the team uses a more long-term branching strategy, thanks to easier branching management with Git. The use of paired programming has decreased with an increase in the reliance on code review that has been made easier with tools like GitLab.

1.2. Capabilities

The capabilities in the MCATK code base have been steadily growing to make MCATK a more capable Monte Carlo software library. Currently MCATK’s major capabilities include:

• mesh-based and solid-body geometries

• multiple solution algorithms that include: static k_eff and α eigenvalue, time-dependent, fixed-source, radiography, and stochastic systems solution algorithms

• continuous-energy neutron, photon, and photonuclear transport physics

• track length and next-event estimators for tallies

• HDF5 output for tallies

• GPU acceleration support of ray-tracing tallies via the MonteRay library

• MPI parallelism, including: on-node shared memory and domain decomposition

• Python, C++, C and Fortran bindings

2. Features

The current major features of the MCATK library are listed in the following sections.

2.1. Geometry

MCATK supports both mesh geometry and solid-body geometry [9], allowing users flexibility in their geometry modeling. The mesh geometries available in MCATK include 3-D Cartesian, 2-D cylindrical, and 1-D spherical meshes. The solid-body geometry is based on a hierarchical scene graph structure, where each solid in the scene is assigned to a node in a tree structure. Each node is a child of, or contained by, the node above it, similar to the structure used in CERN’s (Conseil Européen pour la Recherche Nucléaire) ROOT geometry [10]. Nodes in this geometry hierarchy have the properties:

• each node is assigned a geometric primitive.

• Each node can have other nodes registered as children.

• A node can have an optional transformation that will be inherited by its children.

• Children are allowed to overlap, but MCATK applies an order of precedence.

• Parent containment and child precedence can be used to replicate combinatorial operators.

• Additionally, a node can be assigned a mesh object as the geometric primitive, but these nodes cannot have any children.

MCATK can import its own geometry mesh format and the LNK3DNT mesh geometry format. The LNK3DNT mesh geometry is a mesh format supported by the PARTISN [11] deterministic transport software at Los Alamos, historically from the THREEDANT code [12]. The LNK3DNT format is a common format for transport codes at Los Alamos. It is a structured mesh format that supports both 3-D Cartesian, 2-D cylindrical, and 1-D spherical mesh representations. Both MCATK and MCNP can import LNK3DNT mesh geometry. MCNP can also create LNK3DNT mesh files, providing a mechanism to convert MCNP geometry to a format MCATK can read.

MCATK provides both a 2-D plotting tool and a 3-D rendering tool. The 2-D plotting tool is shown in Figure 1. A rendering from the 3-D rendering tool is shown in Figure 2.

Fig. 2. Rendering of the Tinkertoy 2 critical assembly, created with the MCATK rendering tool.

2.2. Solution algorithms

MCATK has three main solution algorithms, or modes: a static k_eff eigenvalue mode, a static α eigenvalue mode, and a time-dependent mode [13]. MCATK contains a number of other solution types, which are subsets of these 3 solution types, including: fixed source, radiography, and stochastic system solution algorithms.

The k_eff eigenvalue mode is analogous to the kcode mode in the MCNP software [14]. k_eff is the effective multiplication factor, which is generally considered to be the number of neutrons in a given generation divided by the number of neutrons in the previous generation. Or in terms of the Boltzmann transport equation, k_eff can be thought of as a scaling factor on the multiplication term:

$$\begin{array}{c}\hfill \mathbf{L}\psi +{\mathrm{\Sigma }}_t\psi =\mathbf{S}\psi +\frac{1}{k_{}\text{eff}}\mathbf{F}\psi \end{array}$$(1)

where Lψ is the leakage term, Σ_tψ is the absorption term, Sψ is the scattering term, and Fψ is the fission term. MCATK has three estimators for k_eff, including a track length estimator, a collision estimator, and an absorption estimator. MCATK uses a radius of gyration estimate as a measure of stationarity [15].

MCATK’s k_eff eigenvalue mode has been validated using 1-group analytic solutions developed by Sood [16]. And the k_eff eigenvalue mode has been compared with MCNP software for criticality benchmarks from the International Criticality Safety Benchmark Evaluation Project (ICSBEP) Handbook [17].

MCATK’s k_eff algorithm mode has been validated against MCNP software, version 6.3.0, for criticality benchmarks from MCNP’s original 31 benchmark model criticality validation suite [18, 19], which model experiments from the ICSBEP. These 31 problems have been modeled with solid-body geometry in both MCNP software and MCATK. Some of the problems have also been modeled in both MCNP software and MCATK with the LNK3DNT mesh representation. And some spherical models have also been modeled directly in MCATK’s spherical mesh representation. Thus the additional mesh representations increased the original 31 problems from the MCNP benchmark model criticality validation suite to 55 problems.

The k_eff and CPU timing results for 37 of the 55 problems are provided in Table 1 for the solid-body and LNK3DNT mesh geometries. The spherical mesh results are not provided as they are nearly identical to the results of the solid-body geometry – due to statistical variation. The timings for single CPU performance are provided in the last column relative to MCNP. The results show that MCATK is generally comparable to the performance of MCNP software for most problem types. MCATK is approximately 2–4 times faster for very simple geometries like Godiva. This indicates that MCATK’s cross-section lookup may be faster than the MCNP code, however, the MCNP code is also tallying more information than MCATK which may explain the difference. For future work the performance difference will be investigated. But, MCATK is 7 times slower than MCNP software for the test problem SB5RN3, which is a closely packed hexagonal lattice. As MCATK does not have a lattice geometry capability it cannot reach the performance of MCNP software for closely packed lattices. But MCATK can approach the performance of MCNP software for lattices that can be arranged in non-overlapping rows, such as BAWXI2, ICT2C3, and PNL33, by nesting the cells in a lattice row into row nodes. This approach is not possible for closely packed lattices. MCATK is approximately 2–4 times faster than MCNP software for LNK3DNT mesh geometries. It is unclear why MCATK is faster than MCNP in some cases but is likely due to other simulation physics rather than geometry tracking algorithm differences.

The results of all 55 problems were analyzed statistically and found to be in agreement with MCNP software within expected normal statistics. To perform this analysis, we compared the number of tests agreeing within 1, 2, and 3 standard deviations to the expected values based on a normal distribution. The metric is the k_eff difference between the simulations divided by the combined uncertainty:

$$\begin{array}{c}\hfill \frac{\mathrm{\Delta }{k}_{}\text{eff}}{{\sigma }_{\mathrm{\Delta }}}=\frac{k_{}\text{eff,MCATK}-k_{}\text{eff,MCNP}}{\sqrt{{\sigma }_{}{\text{MCATK}}^2+{\sigma }_{}{\text{MCNP}}^2}}\end{array}$$(2)

The results from this statistical analysis are presented in Table 2.

While the use of k_eff calculations is common for criticality safety applications, in the experimental world k_eff is inferred. What is actually measured is reactor period (λ) or subcritical multiplication. The system time constant, usually denoted as α where α = 1/λ, describes the rate of change of neutrons in a static system. Using the expression for α in the Boltzmann transport equation from [14] results in the transport equation with α providing a loss term:

$$\begin{array}{c}\hfill \mathbf{L}\psi +({\mathrm{\Sigma }}_t+\frac{\alpha }{\nu })\psi =\mathbf{S}\psi +\frac{1}{k^{\prime }}\mathbf{F}\psi ,\end{array}$$(3)

where Lψ is the leakage term, Sψ is the scattering term, (1/k′)Fψ is the fission term, Σ_tψ is the absorption term. (α/ν)ψ is often referred to as the “time loss” term. MCATK calculates α through an iterative procedure until k′=1.

2.3. Stochastic systems solution algorithms

MCATK contains several solution modes that are useful for describing stochastic systems, including multiplication, probability of initiation (POI), probability of extinction (POE) [20, 21], probability of survival (POS) [22], and quasi-static [23] modes. These solution modes provide dynamic analysis for quantities of interest that cannot be obtained with static analysis. Such analyses include time-integrated or time-dependent quantities on a chain-by-chain basis such as total dose or instantaneous dose rate, burst wait time of pulsed reactor experiments, and the probability of survival (probability that a fission chain persists until some designated time). The POE and POS as a function of time for a simulated criticality excursion are shown in Figure 3. These solution modes are useful for modeling criticality accident scenarios, criticality experiments, transient reactor analysis, pulsed reactors, and subcritical systems.

Fig. 3. Probability of initiation (POI) and probability of survival (POS) for a criticality excursion, from Trahan et al [22].

The C5G7-TD benchmark has been used to verify MCATK’s quasi-static solution modes. C5G7-TD is a light water reactor with 16 UO2 and MOX fuel assemblies [24]. The details of MCATK’s modeling of this benchmark are provided in reference [23].

2.4. Sources

MCATK now supports a wide array of source options. MCATK allows sources to be specified in a variety of formats. The source can be provided on the geometry mesh or defined explicitly by the user. The available source definitions are listed in Table 3. MCATK also supports a surface source read/write capability in order to import or export sources to MCNP software.

2.5. Nuclear data and the Avalanche library

MCATK uses continuous-energy neutron and photon data from data processed through the NJOY nuclear data processing code. Traditionally MCATK directly read ACE (A Compact ENDF) [25] formatted Evaluated Nuclear Data Files (ENDF) [26]. Recently, the MCATK toolkit has been updated to utilize a new nuclear data interface library called Avalanche. Avalanche accesses nuclear data through the use of the ACEtk [27] library which reads ACE-formatted nuclear data files. The Avalanche C++ library provides a standardized interface for accessing continuous-energy photoatomic and nuclear data across various Monte Carlo particle transport codes. Avalanche offers several key benefits, including:

• a uniform interface for accessing both photoatomic and nuclear data

• support for both C₊₊ and Python interfaces

• data format independence, ensuring that the API remains the same regardless of the underlying data storage format.

While Avalanche is agnostic to the data storage format, it still requires populating its data structures from a specific format. Currently, only ACE files are supported. Notably, Avalanche enables access to its data on Graphics Processing Units (GPUs). To optimize performance, Avalanche stores its data in GPU-friendly containers with minimal pointer indirections, leveraging modern C++17 features like variants to achieve polymorphism without relying on virtual base classes and inheritance patterns. Additionally, users can opt to store data in single precision, reducing memory requirements and potentially improving performance on GPUs. Additionally, Avalanche allows for linearizing log-log interpolated photoatomic data that is provided by ACE files. Finally, Avalanche provides all necessary data structures and accessors used during transport in header files, allowing for maximum inlining and optimized performance on GPUs.

2.6. Hash-based cross-section lookup

MCATK uses a hash-based cross-section lookup to reduce the computational cost of cross-section access. The cross-section hash data table is based on the concept described by Brown [28]. By exploiting the binary representation of floating-point numbers according to the IEEE 754 standard, MCATK replaces the computationally expensive mathematical logarithm function with a bit manipulation operation. This approach significantly reduces the computational cost of cross-section access. The bit shifting hash function used in MCATK is listed in pseudocode in Figure 4. This capability has recently been moved into to the Avalanche nuclear reaction library.

Fig. 4. Pseudocode for the bit-shift-based hash function used to convert energy to a logarithmic integer index used in place of the mathematical log function.

2.7. Tallies

MCATK supports various types of tallies, including collision event tallies, surface tallies, track length tallies, next-event estimators (also known as point detectors) [29], and volumetric ray casting estimators [30].

Tally binning is facilitated through the use of tally filters. A filter is a function that maps a particle from anywhere in phase space to a specific tally bin, if applicable. Users can employ filters to control binning based on properties such as energy, time, particle type, direction, material, and cell location.

MCATK offers several score modifiers to customize tallies. Available options include:

• flux

• heating

• neutron production (fission-induced and net)

• photon production (neutron-induced and total)

• absorption

• cross-section multiplier

• energy multiplier

• inverse speed multiplier

• histogram and linear response functions.

MCATK’s tallies have been verified against MCNP software for most tally types. The photon next-event estimator has been compared to MCNP software using a series of cylinder tests [31]. For these cylinder tests, shown in Figure 5, photons are incident on a cylinder, 400 cm long and 40 cm in diameter. A point detector is located 400 cm off axis and 175 cm from the face of the cylinder. Cylinders of differing material and density comprise multiple tests. The MCATK results compared with MCNP software are shown in Table 4. The MCATK results are within two standard deviations of the MCNP code results, except the 7.5 × 10⁻⁵ g/cc Tungsten case which has a large difference in the number of standard deviations but the result is calculated to the limits of numeric precision and the difference is on the order of 10⁻¹⁴.

Fig. 5. Geometry of photon next-event estimator cylinder tests [31].

2.8. GPU acceleration with MonteRay

MCATK can use the MonteRay ray tracing engine to accelerate the calculation of ray tracing tallies, such as next-event estimators [30]. This is done by using a GPU to perform the ray tracing calculations while the random walk continues to be performed on the CPU. MonteRay can provide significant performance gains for ray tracing tallies on GPU hardware. A simulated radiograph of the Zubal head phantom created with MCATK utilizing MonteRay on GPU hardware is shown in Figure 6 [32, 33]. For this simulation the collided component of the radiograph converges much slower than the source component, thus the performance of the calculation of the collided component dominates. The performance improvement of the collided component with MonteRay is 15×.

Fig. 6. Radiograph of a head phantom generated with MonteRay on GPU hardware providing a 15× performance increase for the collided component. The source component is more than 100× more intense than the collided component, but the calculation time of the collided component is 10.3× more than the source component.

2.9. Parallel computing capabilities

MCATK has historically used MPI for parallel computing capabilities for message passing between distributed memory multiprocessors using the Boost C++ MPI interface [34, 35]. MCATK uses MPI for distributing particle histories to all MPI processes. MCATK also has multiple options for distributing the geometry information to the processes. This includes replicating the entire geometry among the computing nodes, known as domain replication (see Fig. 7b). For mesh based geometries only, MCATK can also distribute the mesh among computing nodes in order to support modeling of very large mesh geometries, known as domain decomposition (see Fig. 7c). For domain decomposition, MCATK passes the particles among the nodes owning the local mesh geometry data. MCATK can also support mixing both domain decomposition and domain replication, referred to as domain hybrid, allowing some domains to be replicated as specified by the user (see Fig. 7d).

Fig. 7. A notional example of geometry mesh parallel decomposition options available in MCATK. (a) Original 3 × 4 rectangular mesh geometry. The original mesh is shown with no decomposition applied. (b) Domain Replication: The mesh is replicated on four processes, resulting in each process holding an identical copy of the entire geometry. (c) Domain Decomposition: The mesh is decomposed into four segments, each containing a unique portion of the geometry. Each process holds one segment. (d) Hybrid Decomposition: A combination of domain decomposition and replication, where the mesh is divided into two segments, with each process holding two replicas of each segment.

Since 2017, MCATK began to address the issue of threading on multi-processor nodes with shared memory. In contrast to the MCNP software, which uses a combined MPI+OpenMP approach for its hybrid programming model that requires developers to identify private data for each thread, MCATK has taken a different approach by utilizing shared memory to designate which data is shared among all processors on a shared memory compute node. This method allows developers to concentrate on sharing large amounts of static data, like geometry and material information, that do not change during particle transport simulations. This approach to the hybrid “MPI+X” parallel programming model [36], which combines “distributed memory” (MPI) and “shared memory” (X) parallelism, is termed “MPI+MPI” [37]. By focusing on static data, this approach also reduces the risk of race conditions, making it easier to debug. Initially, MCATK employed the Boost interprocess library [35] to share mesh geometry data using shared memory; however, it is now transitioning to use MPI-3 shared memory for sharing both nuclear data and geometry data. In the future, this shared memory method will also be applied to optionally share read/write tally data, offering additional memory savings for large tallies. A detailed study of the performance impacts of storing nuclear data and geometry information in shared memory has not yet been conducted – the performance impacts will be investigated in the future.

3. Future work

MCATK continues to evolve. The near term features to be added to MCATK include:

• measurement of performance efficiency of the shared memory threading implementation.

• Supporting full particle transport on GPU hardware.

• Thick target bremsstrahlung photon creation.

• Electron transport.

• Tracking in boundary representation meshes, such as STL and STEP geometry files.

MCATK and Avalanche are not yet available for users outside of Los Alamos National Laboratory. It is expected that MCATK will be distributed through the Radiation Shielding Information Computational Center (RSICC) in the future, when MCATK’s feature set and user documentation can support widespread use.

Acknowledgments

The authors would like to acknowledge the following people for their contribution to the development of MCATK: Rob Aulwes, George P. Burdell, Simon Greene, H. Grady Hughes, JP Morgan, Lori Pritchett-Sheats, and Chris Werner. AI acknowledgements: OpenAI’s ChatGPT Enterprise (a large language model trained by OpenAI, based on the GPT-4 architecture) [38] and Meta’s open source Llama 3.1 [39] were used to create BibTeX entries, format tables for LaTeX, and provide grammar corrections throughout the document.

Funding

This work is supported by the U.S. Department of Energy Advanced Scientific Computing Program. This work is supported by the U.S. Department of Energy through Los Alamos National Laboratory (LANL) operated by Triad National Security, LLC, for the National Nuclear Security Administration (NNSA) under Contract No. 89233218CNA000001.

Conflicts of interest

The authors declare that they have no competing interests to report.

Data availability statement

This article has no associated data generated and/or analyzed.

Author contribution statement

Writing – original draft: J.S.; Writing – review and editing: T.B., T.A., J.G., S.B., A.T.; Visualization: T.A., J.G., S.N., J.S., T.T.; Formal analysis: T.B., T.A., A.T., T.T., J.S.; Validation: T.A., S.N., T.T., A.T., J.S.; Software: T.B., T.A., S.B., J.G., A.L., S.N., C.S., A.T., T.Z., A.M., T.T., J.S.; Project administration: T.T., A.L., T.A., S.N.

References

All Tables

All Figures

	Fig. 1. MCATK geometry plotter, showing two views of the Tinkertoy 2 critical assembly.
In the text

	Fig. 2. Rendering of the Tinkertoy 2 critical assembly, created with the MCATK rendering tool.
In the text

	Fig. 3. Probability of initiation (POI) and probability of survival (POS) for a criticality excursion, from Trahan et al [22].
In the text

	Fig. 4. Pseudocode for the bit-shift-based hash function used to convert energy to a logarithmic integer index used in place of the mathematical log function.
In the text

	Fig. 5. Geometry of photon next-event estimator cylinder tests [31].
In the text

	Fig. 6. Radiograph of a head phantom generated with MonteRay on GPU hardware providing a 15× performance increase for the collided component. The source component is more than 100× more intense than the collided component, but the calculation time of the collided component is 10.3× more than the source component.
In the text

Fig. 7. A notional example of geometry mesh parallel decomposition options available in MCATK. (a) Original 3 × 4 rectangular mesh geometry. The original mesh is shown with no decomposition applied. (b) Domain Replication: The mesh is replicated on four processes, resulting in each process holding an identical copy of the entire geometry. (c) Domain Decomposition: The mesh is decomposed into four segments, each containing a unique portion of the geometry. Each process holds one segment. (d) Hybrid Decomposition: A combination of domain decomposition and replication, where the mesh is divided into two segments, with each process holding two replicas of each segment.

In the text