# SPDX-FileCopyrightText: 2023-2025 OpenMC contributors and Paul Romano
# SPDX-License-Identifier: MIT
"""This module is for reading ACE-format cross sections. ACE stands for "A
Compact ENDF" format and originated from work on MCNP_. It is used in a number
of other Monte Carlo particle transport codes.
ACE-format cross sections are typically generated from ENDF_ files through a
cross section processing program like NJOY_. The ENDF data consists of tabulated
thermal data, ENDF/B resonance parameters, distribution parameters in the
unresolved resonance region, and tabulated data in the fast region. After the
ENDF data has been reconstructed and Doppler-broadened, the ACER module
generates ACE-format cross sections.
.. _MCNP: https://laws.lanl.gov/vhosts/mcnp.lanl.gov/
.. _NJOY: http://t2.lanl.gov/codes.shtml
.. _ENDF: http://www.nndc.bnl.gov/endf
"""
from __future__ import annotations
from collections import OrderedDict
import enum
from pathlib import Path
import struct
from typing import Tuple, List, Union, Optional, Iterable, TextIO, Any
import numpy as np
from .data import ATOMIC_SYMBOL, gnds_name, EV_PER_MEV, K_BOLTZMANN
from .fileutils import PathLike
from .records import ENDF_FLOAT_RE
import endf
def get_metadata(zaid: int, metastable_scheme: str = 'nndc') -> Tuple[str, str, int, int, int]:
"""Return basic identifying data for a nuclide with a given ZAID.
Parameters
----------
zaid
ZAID (1000*Z + A) obtained from a library
metastable_scheme : {'nndc', 'mcnp'}
Determine how ZAID identifiers are to be interpreted in the case of
a metastable nuclide. Because the normal ZAID (=1000*Z + A) does not
encode metastable information, different conventions are used among
different libraries. In MCNP libraries, the convention is to add 400
for a metastable nuclide except for Am242m, for which 95242 is
metastable and 95642 (or 1095242 in newer libraries) is the ground
state. For NNDC libraries, ZAID is given as 1000*Z + A + 100*m.
Returns
-------
name : str
Name of the table
element : str
The atomic symbol of the isotope in the table; e.g., Zr.
Z : int
Number of protons in the nucleus
mass_number : int
Number of nucleons in the nucleus
metastable : int
Metastable state of the nucleus. A value of zero indicates ground state.
"""
Z = zaid // 1000
mass_number = zaid % 1000
if metastable_scheme == 'mcnp':
if zaid > 1000000:
# New SZA format
Z = Z % 1000
if zaid == 1095242:
metastable = 0
else:
metastable = zaid // 1000000
else:
if zaid == 95242:
metastable = 1
elif zaid == 95642:
metastable = 0
else:
metastable = 1 if mass_number > 300 else 0
elif metastable_scheme == 'nndc':
metastable = 1 if mass_number > 300 else 0
while mass_number > 3 * Z:
mass_number -= 100
# Determine name
element = ATOMIC_SYMBOL[Z]
name = gnds_name(Z, mass_number, metastable)
return (name, element, Z, mass_number, metastable)
def ascii_to_binary(ascii_file: PathLike, binary_file: PathLike):
"""Convert an ACE file in ASCII format (type 1) to binary format (type 2).
Parameters
----------
ascii_file
Filename of ASCII ACE file
binary_file
Filename of binary ACE file to be written
"""
# Read data from ASCII file
with open(str(ascii_file), 'r') as ascii_file:
lines = ascii_file.readlines()
# Set default record length
record_length = 4096
# Open binary file
with open(str(binary_file), 'wb') as binary_file:
idx = 0
while idx < len(lines):
# check if it's a > 2.0.0 version header
if lines[idx].split()[0][1] == '.':
if lines[idx + 1].split()[3] == '3':
idx = idx + 3
else:
raise NotImplementedError('Only backwards compatible ACE'
'headers currently supported')
# Read/write header block
hz = lines[idx][:10].encode()
aw0 = float(lines[idx][10:22])
tz = float(lines[idx][22:34])
hd = lines[idx][35:45].encode()
hk = lines[idx + 1][:70].encode()
hm = lines[idx + 1][70:80].encode()
binary_file.write(struct.pack(str('=10sdd10s70s10s'),
hz, aw0, tz, hd, hk, hm))
# Read/write IZ/AW pairs
data = ' '.join(lines[idx + 2:idx + 6]).split()
iz = np.array(data[::2], dtype=int)
aw = np.array(data[1::2], dtype=float)
izaw = [item for sublist in zip(iz, aw) for item in sublist]
binary_file.write(struct.pack(str('=' + 16*'id'), *izaw))
# Read/write NXS and JXS arrays. Null bytes are added at the end so
# that XSS will start at the second record
nxs = [int(x) for x in ' '.join(lines[idx + 6:idx + 8]).split()]
jxs = [int(x) for x in ' '.join(lines[idx + 8:idx + 12]).split()]
binary_file.write(struct.pack(str('=16i32i{}x'.format(record_length - 500)),
*(nxs + jxs)))
# Read/write XSS array. Null bytes are added to form a complete record
# at the end of the file
n_lines = (nxs[0] + 3)//4
start = idx + _ACE_HEADER_SIZE
xss = np.fromstring(' '.join(lines[start:start + n_lines]), sep=' ')
extra_bytes = record_length - ((len(xss)*8 - 1) % record_length + 1)
binary_file.write(struct.pack(str('={}d{}x'.format(
nxs[0], extra_bytes)), *xss))
# Advance to next table in file
idx += _ACE_HEADER_SIZE + n_lines
[docs]
def get_table(filename: PathLike, name: str = None):
"""Read a single table from an ACE file
Parameters
----------
filename : str
Path of the ACE library to load table from
name : str, optional
Name of table to load, e.g. '92235.71c'
Returns
-------
endf.ace.Table
ACE table with specified name. If no name is specified, the first table
in the file is returned.
"""
tables = get_tables(filename, name)
if name is not None and not tables:
raise ValueError(f'Could not find ACE table with name: {name}')
return tables[0]
# The beginning of an ASCII ACE file consists of 12 lines that include the name,
# atomic weight ratio, iz/aw pairs, and the NXS and JXS arrays
_ACE_HEADER_SIZE = 12
[docs]
def get_tables(
filename: PathLike,
table_names: Optional[Union[str, Iterable[str]]] = None,
verbose: bool = False
):
"""Get all tables from an ACE-formatted file.
Parameters
----------
filename : str
Path of the ACE library file to load.
table_names : None, str, or iterable, optional
Tables from the file to read in. If None, reads in all of the
tables. If str, reads in only the single table of a matching name.
verbose : bool, optional
Determines whether output is printed to the stdout when reading a
Library
Attributes
----------
tables : list
List of :class:`Table` instances
"""
if isinstance(table_names, str):
table_names = [table_names]
if table_names is not None:
table_names = set(table_names)
tables = []
# Determine whether file is ASCII or binary
filename = str(filename)
try:
fh = open(filename, 'rb')
# Grab 10 lines of the library
sb = b''.join([fh.readline() for i in range(10)])
# Try to decode it with ascii
sb.decode('ascii')
# No exception so proceed with ASCII - reopen in non-binary
fh.close()
with open(filename, 'r') as fh:
return _read_ascii(fh, table_names, verbose)
except UnicodeDecodeError:
fh.close()
with open(filename, 'rb') as fh:
return _read_binary(fh, table_names, verbose)
def _read_binary(ace_file, table_names, verbose=False, recl_length=4096, entries=512):
"""Read a binary (Type 2) ACE table.
Parameters
----------
ace_file : file
Open ACE file
table_names : None, str, or iterable
Tables from the file to read in. If None, reads in all of the
tables. If str, reads in only the single table of a matching name.
verbose : str, optional
Whether to display what tables are being read. Defaults to False.
recl_length : int, optional
Fortran record length in binary file. Default value is 4096 bytes.
entries : int, optional
Number of entries per record. The default is 512 corresponding to a
record length of 4096 bytes with double precision data.
"""
tables = []
while True:
start_position = ace_file.tell()
# Check for end-of-file
if len(ace_file.read(1)) == 0:
return tables
ace_file.seek(start_position)
# Read name, atomic mass ratio, temperature, date, comment, and
# material
name, atomic_weight_ratio, kT, *_ = \
struct.unpack('=10sdd10s70s10s', ace_file.read(116))
name = name.decode().strip()
# Read ZAID/awr combinations
data = struct.unpack('=' + 16*'id', ace_file.read(192))
pairs = list(zip(data[::2], data[1::2]))
# Read NXS
nxs = list(struct.unpack('=16i', ace_file.read(64)))
# Determine length of XSS and number of records
length = nxs[0]
n_records = (length + entries - 1)//entries
# verify that we are supposed to read this table in
if (table_names is not None) and (name not in table_names):
ace_file.seek(start_position + recl_length*(n_records + 1))
continue
if verbose:
kelvin = round(kT * EV_PER_MEV / K_BOLTZMANN)
print(f"Loading nuclide {name} at {kelvin} K")
# Read JXS
jxs = list(struct.unpack('=32i', ace_file.read(128)))
# Read XSS
ace_file.seek(start_position + recl_length)
xss = list(struct.unpack(f'={length}d', ace_file.read(length*8)))
# Insert zeros at beginning of NXS, JXS, and XSS arrays so that the
# indexing will be the same as Fortran. This makes it easier to
# follow the ACE format specification.
nxs.insert(0, 0)
nxs = np.array(nxs, dtype=int)
jxs.insert(0, 0)
jxs = np.array(jxs, dtype=int)
xss.insert(0, 0.0)
xss = np.array(xss)
# Create ACE table with data read in
table = Table(name, atomic_weight_ratio, kT, pairs, nxs, jxs, xss)
tables.append(table)
# Advance to next record
ace_file.seek(start_position + recl_length*(n_records + 1))
def _read_ascii(
ace_file: TextIO,
table_names: Optional[Union[str, Iterable[str]]] = None,
verbose: bool = False
):
"""Read an ASCII (Type 1) ACE table.
Parameters
----------
ace_file : file
Open ACE file
table_names : None, str, or iterable
Tables from the file to read in. If None, reads in all of the
tables. If str, reads in only the single table of a matching name.
verbose : str, optional
Whether to display what tables are being read. Defaults to False.
"""
tables = []
tables_seen = set()
lines = [ace_file.readline() for i in range(_ACE_HEADER_SIZE + 1)]
while len(lines) != 0 and lines[0].strip() != '':
# Read name of table, atomic mass ratio, and temperature. If first
# line is empty, we are at end of file
# check if it's a 2.0 style header
if lines[0].split()[0][1] == '.':
words = lines[0].split()
name = words[1]
words = lines[1].split()
atomic_weight_ratio = float(words[0])
kT = float(words[1])
commentlines = int(words[3])
for _ in range(commentlines):
lines.pop(0)
lines.append(ace_file.readline())
else:
words = lines[0].split()
name = words[0]
atomic_weight_ratio = float(words[1])
kT = float(words[2])
datastr = ' '.join(lines[2:6]).split()
pairs = list(zip(map(int, datastr[::2]),
map(float, datastr[1::2])))
datastr = '0 ' + ' '.join(lines[6:8])
nxs = np.fromstring(datastr, sep=' ', dtype=int)
# Detemrine number of lines in the XSS array; each line consists of
# four values
n_lines = (nxs[1] + 3)//4
# Ensure that we have more tables to read in
if (table_names is not None) and (table_names <= tables_seen):
break
tables_seen.add(name)
# verify that we are supposed to read this table in
if (table_names is not None) and (name not in table_names):
for _ in range(n_lines - 1):
ace_file.readline()
lines = [ace_file.readline() for i in range(_ACE_HEADER_SIZE + 1)]
continue
# Read lines corresponding to this table
lines += [ace_file.readline() for i in range(n_lines - 1)]
if verbose:
kelvin = round(kT * EV_PER_MEV / K_BOLTZMANN)
print(f"Loading nuclide {name} at {kelvin} K")
# Insert zeros at beginning of NXS, JXS, and XSS arrays so that the
# indexing will be the same as Fortran. This makes it easier to
# follow the ACE format specification.
datastr = '0 ' + ' '.join(lines[8:_ACE_HEADER_SIZE])
jxs = np.fromstring(datastr, dtype=int, sep=' ')
datastr = '0.0 ' + ''.join(lines[_ACE_HEADER_SIZE:_ACE_HEADER_SIZE + n_lines])
xss = np.fromstring(datastr, sep=' ')
# When NJOY writes an ACE file, any values less than 1e-100 actually
# get written without the 'e'. Thus, what we do here is check
# whether the xss array is of the right size (if a number like
# 1.0-120 is encountered, np.fromstring won't capture any numbers
# after it). If it's too short, then we apply the ENDF float regular
# expression. We don't do this by default because it's expensive!
if xss.size != nxs[1] + 1:
datastr = ENDF_FLOAT_RE.sub(r'\1e\2\3', datastr)
xss = np.fromstring(datastr, sep=' ')
assert xss.size == nxs[1] + 1
table = Table(name, atomic_weight_ratio, kT, pairs, nxs, jxs, xss)
tables.append(table)
# Read all data blocks
lines = [ace_file.readline() for i in range(_ACE_HEADER_SIZE + 1)]
return tables
[docs]
class TableType(enum.Enum):
"""Type of ACE data table."""
NEUTRON_CONTINUOUS = 'c'
NEUTRON_DISCRETE = 'd'
THERMAL_SCATTERING = 't'
DOSIMETRY = 'y'
PHOTOATOMIC = 'p'
PHOTONUCLEAR = 'u'
PROTON = 'h'
DEUTERON = 'o'
TRITON = 'r'
HELIUM3 = 's'
ALPHA = 'a'
[docs]
@classmethod
def from_suffix(cls, suffix: str) -> TableType:
"""Determine ACE table type from a suffix.
Parameters
----------
suffix : str
Single letter ACE table designator, e.g., 'c'
Returns
-------
TableType
ACE table type
"""
for member in cls:
if suffix.endswith(member.value):
return member
raise ValueError(f"Suffix '{suffix}' has no corresponding ACE table type.")
[docs]
class Table:
"""ACE cross section table
Parameters
----------
name : str
Full ACE table identifier, e.g., '92235.70c'.
atomic_weight_ratio : float
Atomic mass ratio of the target nuclide.
kT : float
Temperature of the target nuclide in [MeV]
pairs : list of tuple
16 pairs of ZAIDs and atomic weight ratios. Used for thermal scattering
tables to indicate what isotopes scattering is applied to.
nxs : numpy.ndarray
Array that defines various lengths with in the table
jxs : numpy.ndarray
Array that gives locations in the ``xss`` array for various blocks of
data
xss : numpy.ndarray
Raw data for the ACE table
Attributes
----------
data_type : TableType
Type of the ACE data
temperature : float
Temperature of the target nuclide in [K]
zaid : int
ZAID identifier of the table, e.g., 92235
nxs : numpy.ndarray
Array that defines various lengths with in the table
jxs : numpy.ndarray
Array that gives locations in the ``xss`` array for various blocks of
data
xss : numpy.ndarray
Raw data for the ACE table
"""
def __init__(self, name: str, atomic_weight_ratio: float, kT: float,
pairs: List[Tuple[int, float]],
nxs: np.ndarray, jxs: np.ndarray, xss: np.ndarray):
self.name = name
self.atomic_weight_ratio = atomic_weight_ratio
self.kT = kT
self.pairs = pairs
self.nxs = nxs
self.jxs = jxs
self.xss = xss
@property
def zaid(self) -> int:
return int(self.name.split('.')[0])
@property
def data_type(self) -> TableType:
xs = self.name.split('.')[1]
return TableType.from_suffix(xs[-1])
@property
def temperature(self) -> float:
return self.kT * EV_PER_MEV / K_BOLTZMANN
def __repr__(self) -> str:
return f"<ACE Table: {self.name} at {self.temperature:.1f} K>"
[docs]
def interpret(self, **kwargs) -> Any:
"""Get high-level interface class for the ACE table
Parameters
----------
**kwargs
Keyword-arguments passed to the high-level interface class
Returns
-------
Instance of a high-level interface class, e.g.,
:class:`endf.IncidentNeutron`.
"""
if self.data_type == TableType.NEUTRON_CONTINUOUS:
return endf.IncidentNeutron.from_ace(self, **kwargs)
else:
raise NotImplementedError(f"No class implemented for {self.data_type}")
[docs]
def get_libraries_from_xsdir(path: PathLike) -> List[Path]:
"""Determine paths to ACE files from an MCNP xsdir file.
Parameters
----------
path
Path to xsdir file
Returns
-------
List of paths to ACE libraries
"""
xsdir = Path(path)
# Find 'directory' section
with open(path, 'r') as fh:
lines = fh.readlines()
for index, line in enumerate(lines):
if line.strip().lower() == 'directory':
break
else:
raise RuntimeError("Could not find 'directory' section in MCNP xsdir file")
# Handle continuation lines indicated by '+' at end of line
lines = lines[index + 1:]
continue_lines = [i for i, line in enumerate(lines)
if line.strip().endswith('+')]
for i in reversed(continue_lines):
lines[i] = lines[i].strip()[:-1] + lines.pop(i + 1)
# Create list of ACE libraries
libraries = {}
for line in lines:
words = line.split()
if len(words) < 3:
continue
lib = (xsdir.parent / words[2]).resolve()
if lib not in libraries:
# Value in dictionary is not used, so we just assign None. Below a
# list is created from the keys alone
libraries[lib] = None
return list(libraries.keys())
def get_libraries_from_xsdata(path: PathLike) -> List[Path]:
"""Determine paths to ACE files from a Serpent xsdata file.
Parameters
----------
path
Path to xsdata file
Returns
-------
List of paths to ACE libraries
"""
xsdata = Path(path)
with open(xsdata, 'r') as xsdata_file:
# As in get_libraries_from_xsdir, we use a dict for O(1) membership
# check while retaining insertion order
libraries = OrderedDict()
for line in xsdata_file:
words = line.split()
if len(words) >= 9:
lib = (xsdata.parent / words[8]).resolve()
if lib not in libraries:
libraries[lib] = None
return list(libraries.keys())