NJOY+NCrystal: an open-source tool for creating thermal neutron scattering libraries

In this work we present NJOY+NCrystal, a tool to generate thermal neutron scattering libraries with support for coherent and incoherent elastic components for crystalline solid materials. This tool, which is a customized version of NJOY, was created by modifying the nuclear data processing program NJOY to call the thermal scattering software library NCrystal, and includes a proposed change in the ENDF-6 format to store both the coherent and incoherent elastic components. Necessary changes to enable this format in NJOY, as well as to sample it in the OpenMC Monte Carlo code, are detailed here. Examples of materials that are coherent-dominant, incoherent-dominant, and mixed elastic scatterers are presented, as well as the creation of novel libraries for MgH$_2$ and MgD$_2$, that are under consideration as advanced neutron reflectors in the HighNESS project at the European Spallation Source. NJOY+NCrystal simplifies greatly the process to generate thermal scattering libraries (TSL) and this is exemplified with 213 new and updated TSL evaluations.

data that is later used by NJOY to produce thermal scattering libraries for poly-crystalline (or sometimes referred to as crystals treated in the powder approximation) materials. Additionally, we include modifications to support the proposed mixed elastic format for tsl-ENDF files and produce thermal scattering .ACE files, in the here proposed format with both coherent and incoherent elastic scattering components. With NJOY+NCrystal, the user isn't anymore limited to materials with hard-coded options, but any solid poly-crystalline material can be calculated with support for coherent and incoherent elastic components in the current ENDF format, as well as in the proposed mixed elastic format. The Monte Carlo code OpenMC [4] was modified as well to support this new format of .ACE files, which makes this a complete implementation, starting from microscopic calculations and ending with Monte Carlo sampling.

Theory
As presented in [9], the scattering of neutrons can be represented in terms of the double differential scattering cross section as: where is the cross section, dΩ ′ is the scattering angle, dE' is the scattered energy of neutron, and ′ are wavenumbers for the incoming and scattered neutron, ℏ ⃗ is the momentum transfer, equal to ℏ ⃗ = ℏ ⃗ ′ −ℏ ⃗ . ( ⃗ , ) is the scattering function, which for a crystalline system is divided into an equivalent subsystem and can be represented as: where represents the subsystem size, ′ are the average scattering lengths over the subsystem, and ⟨ , ′ ⟩ is the expectation value at thermal equilibrium which correlates the position of the nucleus at time with the position of the nucleus ' at time 0 and it has form ⟨ − ⃗ ⋅ ⃗ ′ (0) − ⃗ ⋅ ⃗ ( ) ⟩. The scattering function of such a system can then be represented as a sum of coherent and incoherent parts: where ⟨ , ′ ⟩ is dependent on positions of the nuclei in a crystal. This correlation can be further simplified by applying the harmonic approximation which consists in separating the position operator in displacements ⃗ from equilibrium positions ⃗ . The resulting correlation between displacements can be expanded in a Taylor series: The first term in this expansion corresponds to elastic scattering. The correlation in this term is known as the Debye-Waller function: For isotropic systems (i.e. polycrystalline samples), integration over all orientations results in a Q-dependent Debye-Waller function: where ⟨ 2 ⟩ is the mean-squared displacement (MSD).

Incoherent elastic scattering
For elastic scattering, with = ′ , the incoherent elastic double differential scattering cross section is equal to: where the bound incoherent cross section, , is equal to = 4 ( 2 − ( ) 2 ). In NJOY, the incoherent elastic cross section is calculated in the THERMR module as: with: where E is the incident energy of the neutron, is the cosine of the scattering angle, and ENDF is the Debye-Waller integral in eV −1 and is computed from the phonon spectrum in the LEAPR module of NJOY: ENDF and are stored in the incoherent elastic section of the ENDF file by LEAPR if the incoherent elastic cross section is present. The integrated incoherent elastic cross section is equal to:

Coherent elastic scattering
The scattering function for the coherent elastic component can be written as: where is the volume of the unit cell,| (⃗ ℎ )| is the form factor of the unit cell, and ⃗ ℎ is a point in the reciprocal lattice. In the powder approximation, where the crystal grains appear with uniformly randomised orientations, the integrated coherent elastic scattering cross section can be written as: where ℎ = ℏ 2 2 ℎ 8 is the energy threshold of the Bragg edge caused by the family of planes with Miller indices (ℎ ), and ℎ is the d-spacing of the planes. In the ENDF-6 format this information is stored as pairs ℎ , ℎ ℎ el,coh .

Inelastic scattering
The ≥ 1 terms in the expansion of eq. 6 correspond to inelastic scattering. In NJOY+NCrystal this is handled by the LEAPR module, in which the double differential scattering cross section is defined in the incoherent approximation as: where is the cosine of the scattering angle, and LEAPR ( , ) is the scattering function. This is different from the previous definition in eq. 1 because in the incoherent approximation the bound atom scattering cross section is factored out. The variables and are related respectively to the momentum transfer and energy transfer: where A is the ratio of the mass of the scattering atom to the neutron mass. In the incoherent and Gaussian approximation the ( , ) is defined as: where: with: where ( ) is the phonon spectrum, and̂ is the time measured in units of ℏ∕ .
The total inelastic cross section is obtained by integrating eq. 16 in outgoing energy and scattering angle:

Total cross section
The total interaction probability for an incident energy will be given by the sum of the absorption cross section (representing the probability of all events that are not scattering), the inelastic cross section (eq. 21, and the coherent (eq. 15 and incoherent (eq. 13) elastic cross sections, if present:

Implementation
NJOY+NCrystal was motivated mainly by a need to implement the mixed elastic format proposal in an open source code, such as NJOY. Additionally, the motivation was to provide an implementation of both a proper and user-friendly handling of the coherent and incoherent elastic cross section calculations, in both the current and proposed format, in NJOY by utilizing NCrystal. Therefore, NJOY+NCrystal is a customized version of NJOY that relies on NCrystal to provide the necessary information for the calculations of the coherent and incoherent elastic components, whereas the inelastic component will be calculated in LEAPR module of NJOY+NCrystal. Additionally, both NJOY and NCrystal utilize CMake and GitHub Repository's for the installation, hence making the integration smoother.
The changes to NJOY can be split into two parts. The first part is to provide an interface between LEAPR and

LEAPR modifications
In the current LEAPR module of NJOY, which is used to prepare the scattering law in ENDF-6 format, there are a few hard coded options for the coherent elastic component, given by the iel flag in the input deck: iel = 0 for no coherent elastic contribution, and options iel = 1-6 for hard coded materials (Table 1). For the implementation of the mixed-elastic proposal, a new option (iel = 100, hereinafter referred to as MEF option) was added, which calculates and stores the parameters for the coherent and incoherent elastic components, as described in the proposed format. Since this change has not yet been widely adopted, an additional (iel = 99) option is provided to control the output of elastic data in the tsl-ENDF file using the current ENDF format (CEF), and hereinafter iel = 99 option will be referred to as CEF option.  Beryllium (hard coded in LEAPR). 3 Beryllium oxide (hard coded in LEAPR). 4 Aluminum (hard coded in LEAPR). 5 Lead (hard coded in LEAPR). 6 Iron (hard coded in LEAPR). … Reserved.

99
Compute TSL in the current ENDF format (CEF). 100 Compute TSL in the proposed mixed elastic format (MEF).

Mixed elastic format, MEF
Coherent scattering is a property of the system (Equation 4), whereas incoherent scattering can be assigned individually to atoms (Equation 5). In the ENDF-6 format, the scattering has to be decomposed into the contribution from the different elements. With this in mind, in the MEF option, the thermal scattering library for each element that composes a material contains two parts: coherent scattering divided by the total number of atoms in the system, and the incoherent scattering that corresponds to the atom (as can be seen in Figure 1): NCrystal LEAPR ENDF-6 TSL File (mixed elastic)

Current ENDF format, CEF
In the CEF option, only the coherent or incoherent scattering model can be used, but not both. For systems with only one type of atom (e.g. Al, Fe, Si), the major component is stored but it is scaled to the total bound scattering cross section (as can be seen in Figure 2) to ensure the proper asymptotic behavior in the epithermal range: NCrystal LEAPR ENDF-6 TSL File (coherent elastic) ENDF-6 TSL File (incoherent elastic) or Figure 2: Simplified flow of NJOY+NCrystal when the CEF option is used for a monatomic scatterer.
For systems with more than one atom, the contributions are sorted and the atom with the minimum incoherent contribution (called designated coherent, DC) is stored in the coherent approximation. Its incoherent contribution, if it exists, is distributed among the rest of the elements which are stored in the incoherent approximation (as can be seen in Figure 3): The first term in this equation is stored in the elastic component of the DC atom in the coherent approximation, and each term of the sum in the second term is stored in the rest of the atoms in the incoherent approximation.
Since incoherent scattering of the designated coherent element is being redistributed between the other atoms in the molecule/compound, the approximation made in the last step impĺies using the Debye-Waller factor of the rest of the components to account for the missing incoherent scattering of the designated coherent element. If the mixed-elastic proposal cannot be used, the calculation of the CEF option using NJOY+NCrystal provides a reasonable approximation within the limitations of the format. It is important to note that the CEF option is currently supported in standard versions of MCNP and PHITS, whereas the MEF option requires modifications. Examples of this are provided in Section 4. It is also important to note that both the CEF and the MEF options will only provide the right answer if they are used in the correct stoichiometry.

Wrapper function
For the CEF and MEF options, LEAPR calls a wrapper function (since NJOY2016 is written in Fortran while NCrystal is in C++ programing language) that for the coherent elastic component obtains a list of Bragg edges, given as pairs of energy and cross section data for each Bragg edge from NCrystal. The wrapper function also passes to LEAPR the values for the bound coherent and incoherent elastic cross sections, as well as the free atom cross section for the nuclide that is being processed in LEAPR, which is used for consistency between the inelastic and elastic components. These cross sections are obtained from a builtin library inside NCrystal, which includes data for all the major elements and isotopes. Since the elastic component is handled by NCrystal and the inelastic by LEAPR module of NJOY+NCrystal, an additional consistency check is performed by comparing the mean squared displacement (MSD) values from the two codes for the nuclide that is being processed in LEAPR. The wrapper function passes to LEAPR the NCrystal MSD value and both NCystal and NJOY+NCrystal values for MSD are stored in the output file of the LEAPR module. Since both the inelastic and elastic component exhibit dependence on the phonon spectrum (for the elastic component the dependence is through Debye-Waller integral as seen in Equations 8,11,12, and for the inelastic dependence can be seen in Equations 18,19,20), the check of the MSD values is also a check on whether the same phonon spectrum is used in both the NCrystal and NJOY+NCrystal input files. The LEAPR value of MSD can be calculated from ENDF , as defined in Equation 12, as follows: where is the mass of neutron in units of eV ⋅ ps 2 ⋅ Å −2 , and ℏ 2 is Planck's constant in units of eV ⋅ ps.

THERMR/ACER modifications
The modules THERMR and ACER of NJOY were also modified to implement processing of the tsl-ENDF proposed mixed elastic format, as well as to produce .ACE files in a format that is being proposed here, since at the time of the development of the code no official format for mixed elastic implementation in .ACE files has been developed or announced. THERMR reconstructs the double differential cross section and produces pointwise ENDF files (PENDF tapes) from tsl-ENDF files, and ACER processes PENDF files into .ACE files.
THERMR reads the tsl-ENDF file and handles the elastic section according to the lthr flag: lthr=1 for coherent  elastic, lthr=2 for incoherent elastic or here proposed lthr=3 for mixed elastic. There are no changes to the THERMR input deck and the lthr flag is set during the writing process of the tsl-ENDF file in LEAPR.
In the ACER input deck, an additional option for the ielas flag was added. Option ielas=0 processes the PENDF tape as coherent elastic, ielas=1 processes the file as incoherent elastic, and the here proposed ielas=2 as mixed elastic. The .ACE format [11] requires changes to store the additional data for the mixed elastic contributions, and the changes we propose in this work can be seen in Figure 4. In the proposed format, the flag IDPNC=5 is reserved for mixed elastic, and four additional integers are stored in the NXS and JXS control arrays. NCL2 determines the dimension of the second elastic block, ITCE2 gives the position of the energy table for the second elastic block, ITCX2 gives the position of the cross section table for the second elastic block, and ITCA2 gives the position of the angular distribution. In this implementation of the mixed elastic format of .ACE files, the first elastic block is coherent and the second elastic block is an incoherent elastic component.

OpenMC modifications
As a proof of concept implementation we modified OpenMC to read and use .ACE files in the proposed mixed elastic format. The Python API was modified to read mixed elastic .ACE files and as well to read and write HDF5 files in mixed elastic format. The C++ calculation engine was modified to add the second elastic component to the total cross section if it is present and to sample both elastic reactions. The modifications to OpenMC can be found at [12].

Example evaluations
Using NJOY+NCrystal, 213 tsl-ENDF evaluations were created for 112 new or updated materials. A summary of the NJOY+NCrystal library is given in Tables 3-7. In this section, a few materials will be highlighted, which include tin as an example of a coherent scatterer, vanadium as an incoherent scatterer, nickel as mixed elastic scatterer, and  [13,14,15,16]. In summary, the files from the database were used with Phonopy [17] to calculate phonon eigenvalues and eigenfrequencies, which were then utilized with oClimax [18] to extract partial phonon spectra. A detailed explanation of how the phonon spectrum and crystal structure was obtained for each material is provided in the comments section of the tsl-ENDF files.
The experimental data for the validation of the tsl-ENDF files is scarce. Wherever the experimental data were available, tsl-ENDF files were validated against total cross section measurements and diffraction data, while all materials were validated against specific heat capacity curves as a minimum standard for acceptance. The last column of Tables 3-7 contains information on the valid temperature range for each material, references for VDOS curves, and a list of validations applied. For the VDOS curves a reference is provided, where for the ones calculated with Phonopy and oClimax the reference to Kyoto University phonon database is given. The link to the exact files used for each material at Kyoto University database can be found inside the tsl-ENDF files in the Highness Github repository. For validation, the correct crystal structure was a starting point (a reference for each crystal structure can be found inside tsl-ENDF files as well), followed by comparison with the experimental specific heat capacity and if available, total cross section measurements. For and SiC, diffraction data was used as well as a means of validation of the libraries.

Tin
With inc = 0.022 b and coh = 4.870 b, tin is a mainly coherent scatterer. Using NJOY+ NCrystal we have created tsl-ENDF and .ACE files for tin, and compared the calculated total cross section with measurements by Mayer [19,20], retrieved from EXFOR [21] (Figure 5). In this comparison, both scattering and absorption are included, using 2200 = 0.626 b. For completeness, the total cross section from the .ACE files, prepared by NJOY+ NCrystal, was plotted using OpenMC and the results for CEF and MEF options are shown, as well as each scattering component from the MEF option. In Figure 5 it can be observed that the agreement between the experimental total cross sections and the calculated cross sections for NCrystal, and using both CEF and MEF options in NJOY+NCrystal is very good (all three lines are overlaid on top of each other). Therefore, for predominantly coherent materials like tin, the CEF option performs nearly as well as the MEF option, and it can be used even without the mixed-elastic format.

Vanadium
Vanadium is a common example of a mostly incoherent material, with inc = 5.08 b and coh = 0.018 b. The calculated cross sections were obtained in the same manner as for tin.
As can be seen from Figure 6, the calculated total cross sections, including 2200 = 5.08 b, are in excellent agreement with the experimental total cross section measurements from Vertebnyj [22] and Palevsky [23]. The coherent elastic component is negligible, while the incoherent elastic component is dominant. Again, the CEF option provides an adequate representation of the total cross section for Vanadium.

Nickel
Nickel, with inc = 5.2 b and coh = 13.3 b, has a significant contribution from both coherent and incoherent scattering. This, in the current ENDF-6 format, can only be approximated correctly to a certain degree. In Figure 7 we compare the calculated total cross sections, including 2200 = 4.49 b, with experimental data from Allen [24,25] and Palevsky [26,27]. In this plot it can be seen that the agreement of the NCrystal and the total cross section curve calculated with the MEF option with the experimental data is excellent, while the agreement for the CEF option curve is good in the region above the first Bragg edge. Since the CEF option scales the major elastic component, in this case the coherent elastic part (Equation 24), by the total bound incoherent plus coherent scattering cross section, the total cross section below the first Bragg edge only contains the contribution from the inelastic component and hence differs significantly. There is also some discrepancy above the first Bragg edge as well for CEF option, which is due to the fact that coherent elastic scaling provides the correct epithermal asymptotic behavior but only approximates the low energy region. Thus this demonstrates the importance for the ENDF-6 format to support the capability to store both the coherent and incoherent elastic scattering components.
As mentioned in Section 3.3, changes to OpenMC were implemented so that the new format can be used to sample thermal scattering events. Figure 8 shows the results of a simple transmission Monte Carlo calculation, which includes a mono-energetic neutron beam on a 5 mm thick slab of nickel. The total cross section curve was obtained from the ratio of the incident and the transmitted neutron spectrum and we can see that the agreement between the calculation and the experimental data, as well as the total cross section calculated in the MEF option, is excellent. The discrepancy at the lower energies is caused by differences in the absorption cross section in ENDF/B-VIII.0 library ( 2200 = 4.09 b) and NCrystal ( 2200 = 4.49 b).

MgD 2 /MgH 2
As an additional example of the capabilities of NJOY+ NCrystal, related to the HighNESS project, we computed the total scattering neutron cross sections of polycrystalline MgH 2 and its deuterated variant MgD 2 . MgH 2 recently emerged as a good candidate material for cold neutron reflectors [28], while MgD 2 , because of an even lower absorption cross section, is also investigated here.
Phonon spectrum calculations were performed using Density Functional Perturbation Theory (DFPT) [29] with the Perdew-Burke-Ernzerhof (PBE) [30] approximation to the exchange and correlation energy functional. Ultrasoft pseudopotentials [31,32] and a plane wave expansion of Kohn-Sham orbitals up to a kinetic cutoff of 60 Ry was employed, as implemented in the Quantum-ESPRESSO package [33]. In the electronic structure calculations, the Brillouin zone (BZ) integration was performed over a uniform Γ-centered 10×10×14 k-point mesh [34].
MgH 2 crystallizes in the rutile structure (P4 2 /mmm space group) with six atoms per unit cell [35]. Geometry optimization yields the theoretical equilibrium lattice parameters a=4.512 Å and c=3.010 Å which are very close to the experimental values of a=4.5025 Å and c=3.0123 Å as measured from X-ray scattering [35]. The dynamical matrix is computed within DFPT on a 4×4×6 q-mesh. Short range force constants were then obtained by Fourier transforming the dynamical matrices over the q-mesh after subtracting the contribution from point charges as described in [29]. The force constants then allow computing the phonon spectrum on a fine q-mesh [29].
The phonon dispersion relations and the DOS of MgH 2 and MgD 2 are compared in Figures 9-10. The results for MgH 2 are very similar to previous ab-initio works [36,37]. In particular, in [36] the overall validity of the isotropic approximation of the mean square displacement is assessed at different temperatures. The high-frequency phonons, due to the motion of hydrogen, are almost exactly scaled by a factor √ 2 by isotopic substitution while the acoustic part of the spectrum is essentially unchanged.
The neutron-weighted phonon density of states, computed by multiplying the projected DOS with the corresponding total neutron scattering cross section [38], is compared in Figure 11 with the experimental spectrum from inelastic neutron scattering data [39]. The agreement between theory and experiment is overall good, although a sizable error is present in the position of the peak at about 600 cm −1 . A good agreement between theory and experiments is also found for the frequency of the Raman active modes (Γ point) as shown in Table 2.  From the calculated phonon spectra, tsl-ENDF files were created for both MgH 2 and MgD 2 , using the NJOY + NCrystal tool. In Figure 12 it can be seen that the agreement between the measured total cross section data from Granada [28] and the calculated curves (which include an absorption cross section ( H 2200 = 0.3326 b, Mg 2200 = 0.063 b) is excellent.
Since MgH 2 is an incoherent scatterer, due to the large incoherent cross section of hydrogen, iel=99 is a good option 0 200 400 600 800 1000 1200 1400 Frequency [cm 1 ] a.u. Theory Exp. Figure 11: Neutron-weighted theoretical phonon DOS compared with an experimental spectrum from inelastic neutron scattering [39]. The theoretical density of states has been computed using the tetrahedron method using a 20×20×30 q-grid for interpolation and then convoluted with a Gaussian function with standard deviation of 5 cm −1 to better match the experimental resolution.
as well.  Figure 13 shows a comparison between the NCrystal and the NJOY+NCrystal calculated total cross sections for MgD 2 , as well as different scattering components as calculated in NJOY+NCrystal and NCrystal. The total cross section for MgD 2 is orders of magnitude lower than for MgH 2 . This is due to incoherent hydrogen, with a cross section of 82.02 barns, being replaced with coherent deuterium with the cross section of 7.64 barns. From Figure 13 we can see that both the coherent and incoherent elastic components are not negligible in the thermal region and that the agreement between NCrystal, CEF and MEF calculated total cross sections is excellent.

Conclusions
In Although NJOY+NCrystal represents significant improvement upon currently available methods for treatment of thermal neutron scattering in Monte Carlo particle codes, it is still limited by what can be stored in the ENDF-6 format. For a complete utilization of the capabilities of NCrystal, such as modelling single crystals, small angle neutron scattering, coherent one-phonon approximation, multi-phase and defects support, direct calls to NCrystal by Monte Carlo codes would need to be made. This will be the focus of future work.
All the source code presented in this work, as well as the thermal scattering libraries in tsl-ENDF and .ACE format are available in the Github repository of the HighNESS project: https://github.com/highness-eu.

Supplemental information
Supplemental information is available in the online version of this article, including plots of the total cross and heat capacity for all materials, and a text file containing the comments for each model. Numbers in the cross section plots are EXFOR entries. Experimental cross section data in the plots has not been renormalized.

Acknowledgements
This work was funded by the HighNESS project at the European Spallation Source. HighNESS is funded by the European Framework for Research and Innovation Horizon 2020, under grant agreement 951782.