NCrystal : a library for thermal neutron transport

NCrystal : a library for thermal neutron transport

X.-X. Cai T. Kittelmann European Spallation Source ERIC, Sweden DTU Nutech, Technical University of Denmark, Denmark
Abstract

An open source software package for modelling thermal neutron transport is presented. The code facilitates Monte Carlo-based transport simulations and focuses in the initial release on interactions in both mosaic single crystals as well as polycrystalline materials and powders. Both coherent elastic (Bragg diffraction) and incoherent or inelastic (phonon) scattering are modelled, using basic parameters of the crystal unit cell as input.

Included is a data library of validated crystal definitions, standalone tools and interfaces for C++, C and Python programming languages. Interfaces for two popular simulation packages, Geant4 and McStas, are provided, enabling highly realistic simulations of typical components at neutron scattering instruments, including beam filters, monochromators, analysers, samples and detectors. All interfaces are presented in detail, along with the end-user configuration procedure which is deliberately kept user-friendly and consistent across all interfaces.

An overview of the relevant neutron scattering theory is provided, and the physics modelling capabilities of the software are discussed. Particular attention is given here to the ability to load crystal structures and form factors from various sources of input, and the results are benchmarked and validated against experimental data and existing crystallographic software. Good agreements are observed.

journal: Computer Physics Communications\pdfsuppresswarningpagegroup

=1

PROGRAM SUMMARY

Manuscript Title: NCrystal : a library for thermal neutron transport
Authors: X.-X. Cai and T. Kittelmann
Program Title: NCrystal
Journal Reference:
Catalogue identifier:
Licensing provisions: Apache License, Version 2.0 (for core NCrystal).
Programming language: C++, C and Python
Operating system: Linux, OSX, Windows
Keywords:
Thermal neutron scattering, Simulations, Monte Carlo, Crystals, Bragg diffraction
Classification: 4, 7.6, 8, 17
External routines/libraries: Geant4, McStas
Nature of problem:
Thermal neutron transport in structured materials is inadequately supported in popular Monte Carlo transport applications, preventing simulations of a range of otherwise interesting setups.
Solution method:
Provide models for thermal neutron transport in flexible open source library, to be used standalone or as backend in existing Monte Carlo packages. Facilitate validation and work sharing by making it possible to share material configurations between supported applications.

1 Introduction

The modelling of neutron transport through matter dates back to the efforts aimed at tackling neutron diffusion problems in the middle of the twentieth century, closely tied to the introduction of general purpose computers and the inception of the method of Monte Carlo simulations Metropolis (1987). Since then, the needs for transport simulations involving neutrons have expanded, with a wide range of applications in radiation protection, reactor physics, radiotherapy, and scattering instruments at spallation or reactor sources. To facilitate such simulations a plethora of Monte Carlo simulation applications exists today, which for the purposes here can be divided into two categories: general purpose applications Waters et al. (2007); X-5 Monte Carlo Team (2003); Goorley et al. (2013); Agostinelli et al. (2003); Allison et al. (2006, 2016); Böhlen et al. (2014); Sato et al. (2018); Leppänen et al. (2015); Romano et al. (2015), capable of modelling a variety of particle types in flexible geometrical layouts, and specialised applications Lefmann and Nielsen (1999); Willendrup et al. (2004); Zsigmond et al. (2002); Lieutenant et al. (2004); Šaroun and Kulda (1997); Seeger et al. (2002) dealing with the design of neutron scattering instruments. Applications in the latter group is focused on aspects of individual scatterings of thermal neutrons, but generally lack capabilities for dealing with arbitrarily complex geometries and the inclusion of physics of particles other than thermalised neutrons.

On the other hand, the general purpose Monte Carlo applications are typically oriented towards use-cases in fields such as reactor physics, high energy physics, and radiation protection, and as a rule neglect all material structure and effects of inter-atomic bindings – resorting instead to the approximation of treating all materials as a free gas of unbound atoms with no mutual interactions. This approximation is suitable at higher incident energies, but for neutrons it breaks down at the thermal scale where their wavelengths and kinetic energies are comparable to the typical distances and excitational energies resulting from inter-atomic bindings. The only presently available option for improved modelling of thermal neutron interactions involves the utilisation of specially prepared data files of differential cross section data. These “scattering kernels” are only readily available in neutron evaluation libraries, e.g. Brown et al. (2018), for the dozen or so materials which have traditionally played an important role for shielding or moderation purposes at nuclear facilities. For other materials, scattering kernels must be carefully crafted using non-trivial procedures (like the application of delicate and computationally expensive molecular dynamics models Damián et al. (2014)) and in practice this is rarely done. Furthermore, some aspects of thermal neutron scatterings are in practice ignored when modelling is exclusively based on scattering kernels: since nuclear reactor physics has historically mostly been focused on inelastic and multiple-scattering phenomena, it is usually not possible to include a precise model of the a priory highly significant process of Bragg diffraction into the setup for crystalline materials. At thermal energies, Bragg diffraction is the dominant process for many relevant materials, and although it is an elastic scattering process which does not change the neutron energy, it sends neutrons to characteristic solid angles and can therefore critically influence the geometrical reach of neutrons through a given geometry.

It is thus not currently straight-forward to carry out simulations which at the same time incorporate consistent and realistic physics models for neutrons at both high and low energy scales in general materials, while simultaneously supporting flexible geometrical layouts and the treatment of particles other than neutrons. Such simulations would nonetheless be remarkably useful, in particular when considering aspects of neutron scattering instruments Kanaki et al. (2018); Cherkashyna et al. (2014). Here, precise modelling of individual neutron scatterings is crucial in all beam-line components, while detailed understanding of background levels, the effects of shielding, or the workings of neutron detectors, require incorporation of detailed geometrical layouts and additional physics such as modelling of gammas, fast neutrons, and energy depositions resulting from secondary particles released upon neutron capture. Other examples include the design of advanced moderators or reflectors for novel reactors or neutron spallation sources, and proton or neutron-based radiotherapy, in which adverse irradiation from neutrons are abundant and an increasing concern Xu et al. (2008); Newhauser and Durante (2011).

The NCrystal toolkit presented here is aimed at remedying the situation. Rather than introducing an entirely new application, it is designed as a flexible backend, able to plug the gaps where existing Monte Carlo applications are lacking in capabilities for treating thermal neutrons in structured materials. It was originally developed to facilitate the design and optimisation of neutron detectors for the European Spallation Source, presently under construction in Lund, Sweden Garoby et al. (2018); Peggs et al. (2012, 2013); Kittelmann et al. (2014); Kanaki et al. (2018), but is intended to be as widely useful as possible – not only as a service to the community, but also to ultimately ensure a higher level of quality for the toolkit itself due to feedback and validation from a larger potential userbase. Although dedicated utilities for calculating various properties of thermal neutrons exists Kropff and Granada (1977); MacFarlane et al. (2012); Boin (2012); Kittelmann and Boin (2015), the scope of NCrystal is different, delivering at the same time an ambitious set of physics models, a flexible object oriented design, multiple interfaces and bindings, an open approach to development, and not the least a user-friendly method of configuration. The core functionality is implemented as an open source and highly portable C++ library with no third-party dependencies, and comes with language bindings for all versions of C++, C or Python in widespread usage today. The current code is released with version number 1.0.0 and under a highly liberal open source license The Apache Software Foundation (2004).111Optional components adding support for file formats discussed in Sections 4.2 and 4.3 rely on third-party code, available under different open source licenses. Included in the release is additionally tools, examples, data files, and a configuration file for CMake Martin and Hoffman (2015) with which everything can be built and installed. A website Cai and Kittelmann (2018) has been set up for the NCrystal project, on which users will be able to interact with developers, locate future updates, and access relevant documentation.

The present paper will present the NCrystal framework itself, including interfaces, configuration, and data input options. Although a brief overview of physics capabilities will be provided, detailed discussions of the implementation of actual thermal neutron scattering models are beyond the scope of the present paper, and will be the subject of dedicated publications at a later date. After a review of the relevant theoretical concepts in Section 2, the design and implementation of the core NCrystal code is presented in Section 3. Section 4 presents the capabilities for initialising relevant parameters of modelled crystals from various data sources, introduces the library of data files included with NCrystal and discusses how the loaded results have been validated. Section 5 presents the uniform method for material configuration intended for most end-users, and an overview of the existing language bindings and application interfaces is provided in Section 6. Finally, a discussion of future directions and planned improvements is carried out in Section 7.

2 Theoretical background

Rather than intending to be an exhaustive treatment, this section will provide a brief overview of crystals and neutron scattering, introducing relevant concepts, models and terminology needed for the purposes of the present discussions of NCrystal. Naturally, relevant references to more extensive background literature are provided for readers seeking more information.

2.1 Crystals

Although support for other types of materials is likely to be added eventually, the present scope of NCrystal is restricted to crystals. Consequently, a few basic concepts will be introduced in the following. Further details can be found in Sands (1993); Kittel (2004); Waseda et al. (2011); Hahn and of Crystallography (2002).

Informally, a crystal can be described as an arrangement of bound atoms which is built up by the periodic repetition of a basic element in all three spatial directions. This repeated element is known as the unit cell and is typically chosen such that it is the smallest such cell that reflects the symmetry of the structure. In an idealised setting, the crystal structure would be infinite, but in practice it is sufficient that the unit cell is much smaller than the entire structure. The choice of unit cell for a given crystal is not always unique, but it is always possible to select one which is a parallelepiped spanned by three linearly independent basis vectors , and .222Vectors are here and throughout the text donated with arrows (), while the absence of arrows indicate the corresponding scalar magnitudes (). Additionally, unit vectors are donated with hats () and quantum mechanical operators are shown in bold (, ). By convention they are defined by their respective lengths, , , and , and the angles between them: , , and . The crystal structure is completed by the additional specification of the contents of the unit cell, often provided in coordinates relative to the crystal axes. The position of the th constituent in the unit cell is thus given by:

(1)

Crystals exhibit a discrete symmetry under spatial translation with the vectors:

(2)

Apart from a trivial global offset, coordinates of all constituents in the entire crystal are given by for all combinations of and . The translational symmetry means that any function representing features of the crystal system (such as particle densities) will obey:

(3)

Additional symmetries under rotations, reflections or inversions are possible, depending on the detailed shape and contents of the unit cell. The symmetry properties of a crystal are described by the concept of space groups, and it is possible to describe any spatial crystal symmetry by one of only 230 such groups. These groups can be divided into 7 general classes of crystal systems. Loosely ordered from lowest to highest symmetry they are: triclinic, monoclinic, orthorhombic, tetragonal, trigonal, hexagonal, and cubic. The space groups are defined and described exhaustively in Hahn and of Crystallography (2002).

The space group of a given crystal structure naturally depends on the unit cell positions of constituents, . Conversely, it is possible to construct the full list of these positions by applying the symmetry operators of a given space group to a smaller number of so-called Wyckoff positions Wyckoff (1922). It is indeed very common to find crystal structures defined solely in terms of their space group, unit cell shape and dimensions, and a list of elements and associated Wyckoff positions.

The set of points at positions for all integers , , and , constitutes the crystal lattice. In each such lattice, there exist an infinite number of families of evenly spaced parallel planes passing through the lattice points. Each such family of lattice planes is described by a common plane normal and an interplanar distance (also known as the -spacing), with planes passing through the points for all . Particle diffraction in a crystal occurs in the form of reflections from such families of lattice planes. The task of finding and classifying these families is often done in the so-called reciprocal lattice in momentum-space, which is the Fourier transform of the direct lattice. The basis vectors of the reciprocal lattice are given by:

(4)

Where is the unit cell volume. The points in the reciprocal lattice are thus given by:

(5)

As will be motivated in Section 2.3, it can be shown that each point in the reciprocal lattice corresponds to a family of equidistant lattice planes in the direct lattice. The family of planes corresponding to will be indexed by the value (known as its Miller index), and has interplanar spacing and normal .

In most real systems, regions in which the crystal lattice is continuous and the defining translational invariance of Eq. 3 is unbroken, exists only at the microscopic scale.333Exceptions to this exist, like in some synthetic silicon crystals. However, the scattering theory discussed in Section 2.2 is in any case not strictly applicable to such systems. Macroscopic crystals are then built up from these independently oriented grains of perfect crystals (also known as “crystallites”), and the actual distribution of orientations will not only be related to various macroscopic material properties, but will influence interaction probabilities when the system is probed with impinging particles. As long as the grain sizes are small compared to the regions of material being probed, these effects can be accounted for in a statistical manner – for instance by integrating microscopic interaction cross sections over the distribution of crystallite orientations in order to arrive at effective macroscopic cross sections. Such integrations are particularly trivial to perform in the extreme case where the orientation of individual crystallites are completely independent and uniformly distributed over all solid angles. This model, known as the powder approximation, is not only suitable for modelling actual powdered crystalline samples, but can also be used to approximate interactions in polycrystalline materials like metals or ceramics – especially when the level of correlation in crystallite orientation (“texture”) is low or when the setup involves sufficiently large geometries or spread in incoming particles that effects due to local correlations are washed out. An example of this would be the case where polycrystalline metal support structures are placed in various places throughout a neutron instrument. On the other hand, interactions in a highly textured polycrystalline sample placed in a tightly focused beam of probe particles might not be very well described by the powder approximation.

In another extreme, all crystallites are almost co-aligned in so-called mosaic crystals, with the exact distribution around the common axis of alignment given by the mosaicity distribution. In the simple isotropic case, this distribution is a Gaussian function of the angular displacement, with the corresponding width (referred to as the mosaicity of the crystal) having typical values anywhere from a few arc seconds to a few degrees. In NCrystal, such Gaussian mosaic crystals are referred to as “single crystals”, due to the large degree of crystallite alignment throughout the material.

Some materials are better described by other mosaicity distributions, including anisotropic ones where the distribution is not symmetric with respect to all axes. A particular anisotropic distribution supported by NCrystal is one in which crystallites are aligned around one particular axis, but appearing with random rotations around the same axis. Such distributions occur in stratified or layered crystals, like pyrolytic graphite, in which certain planes of atoms tend to have strongly aligned normals in all crystallites, but no strong alignment for the rotation of the same planes around their normals.

2.2 Scattering theory

The theory of thermal neutron scatterings in condensed matter is thoroughly treated in the literature, see for instance Squires (2012); Schober (2014); Marshall and Lovesey (1971). The discussion in the present and following sections is to a large degree inspired by Squires (2012); Schober (2014).

Scattering in materials by thermal neutrons, taken here to loosely mean energies below or at the scale or wavelengths above or at the scale, can be described via point-like neutron-nuclei interactions or magnetic dipole interactions between the neutron and any unpaired electrons in the target material. For unpolarised neutrons and samples, the latter can be largely ignored Sears (1986) and the current discussion will not concern itself with such magnetic interactions, nor are they currently modelled in NCrystal – although there is nothing fundamental preventing their inclusion in the future.

Thermal neutrons carry such low energy that in practice no emission of secondary particles occurs during pure scattering interactions. As a consequence, scattering of unpolarised neutrons can be completely described in terms of the differential cross section for scattering the incident neutron with wavevector into the wavevector . Based on the Born approximation or Fermi’s golden rule, this differential cross section is related to quantum mechanical transition amplitudes by the so-called master equation for thermal neutron scattering:

(6)

Here is the neutron mass, indices and denotes initial and final states respectively, represents a complete set of eigenstates of the target system, represents the initial occupation levels of target states, and is the interaction potential. Furthermore, and represents the total energy change of the neutron and target system respectively. At this point it might be useful to recall that in addition to wavenumber, , the energy state of a non-relativistic neutron might equally be characterised in terms of energy , momentum , velocity , or wavelength , related by , , and .

Before continuing with the evaluation of Eq. 6, it is important to understand its validity and applicability. Based on the Born approximation, the perturbation experienced by the initial neutron wave function during the scattering is assumed to be small, so it is possible to neglect secondary scattering of waves created as a result of the primary scattering. The sum in Eq. 6 implies that the resulting cross sections will grow with the size of the target system, and therefore the validity of the equation will eventually break down as the size of the target system grows and the probability for secondary scattering becomes non-negligible. In reality this issue is circumvented in the context of Monte Carlo transport simulations in which neutron transport is modelled in a series of steps between independent scatterings, with each simulated step length depending on the actual mean free path length predicted by the cross section in Eq. 6.444Specifically, if is the mean free path length (depending on cross sections and material density) and is a pseudo random number in the unit interval, then the step length until next interaction can be modelled as . High cross sections will automatically lead to small step lengths, and therefore only a small part of the sample will be traversed for each application of Eq. 6. Simulations employing such stepping will therefore to some degree account for multiple scattering phenomena like extinction effects in crystal diffraction. This effective subdivision of a target system into smaller decoherent systems is, however, only possible if the linear scale over which any symmetries in the material structure exists is small compared to the mean free path lengths involved, so wavefunctions originating from scattering in separate subsystems will be mutually incoherent. As described in Section 2.1 most real macroscopic crystals consists of independently aligned microscopic crystal grains, and in practice this ensures the necessary decoherence between different parts of the macroscopic system. The exception is the case of perfect macroscopic crystals (or mosaic crystals with vanishing mosaic spread) like some synthetic silicon crystals. For an excellent and detailed discussion about these issues, refer to (Schober, 2014, Sec. 11.5).

In order to proceed with the evaluation of Eq. 6, one must provide suitable representations of both interaction potential and distributions of target states, . When dealing with thermal neutrons and their relatively long wavelength, the short-range neutron-nuclei interactions can effectively be described as being point-like. An effective potential which describes them as such was introduced by Fermi Fermi (1936). For a collection of nuclei located at positions it looks like:

(7)

Here is the scattering length of the th nucleus, an effective parameter capturing the details of the neutron-nucleus interaction which must be determined through measurements. The term scattering length stems from the fact that the corresponding cross section for scattering from a single fixed nucleus is , equal to the classical cross section of a hard sphere of radius . In the presence of nuclear resonances, the scattering length becomes strongly dependent on the neutron energy. However, at the energy scale of thermal neutrons, resonances only exist for a few rare isotopes – and even then mostly at the high end of the thermal scale van Sluijs et al. (2015). For the present purposes, the scattering length thus attain constant values, specific to each type of nuclei. It is additionally possible to use non-zero imaginary components of scattering lengths to represent absorption physics, however in many contexts – including the present – the two types of interactions are dealt with separately and the scattering lengths will therefore be treated as real and constant numbers. Such separation between scattering and absorption processes is valid at the level for thermal neutrons Sears (1984), which is considered acceptable. In fact, as the current scope of NCrystal does not cover physics of nuclear resonances, absorption processes will be dealt with using the simple – but in most cases accurate van Sluijs et al. (2015) – model of absorption cross sections being inversely proportional to the neutron velocity. This scaling can be intuitively interpreted as having the probability of absorption by a given nucleus proportional to the time spent by the neutron near the nucleus, but also follows from more careful reasoning Sears (1984); Breit and Wigner (1936); Bethe (1935).

It is possible to bring Eq. 6 into a form in which the sum over final states of the target system is removed and the potential is replaced with its Fourier transform. This is particularly convenient for potentials involving -functions like Eq. 7, as these transform into constants. The resulting form is:

(8)

where is the momentum transfer, , and is the scattering function defined by:

(9)

The notation is here used to donate operator expectation values in the target system, , with the target state weights, , usually defined by a requirement for the target system to be in thermal equilibrium. The expectation value under the integral in Eq. 9 correlates the position of the nucleus at time with the position of the nucleus at time , and will here be abbreviated as:

(10)

Leading to the shorter expression for the scattering function:

(11)

This definition of the scattering function contains a sum over the target system under consideration (which as discussed must be of a linear size compatible with the applied Born approximation). Most such systems of interests can be divided into a number of statistically equivalent subsystems. That this is possible for crystals is obvious, since each unit cell forms such a subsystem, but subdivisions are generally possible in almost all systems, be they liquid, polymeric or gaseous in nature. The scattering function for a subsystem can be expressed as:

(12)

where now refers to the subsystem size and might for instance represent the number of nuclei in a crystal unit cell, and represents an average performed over an ensemble of equivalent subsystems. Now, neutron-nuclei scattering has the peculiarity that the scattering lengths are isotope and spin-state dependent, whereas the positions of scatterers, the nuclei, are determined by chemical properties, depending on the nuclear charge but otherwise independent of isotope or spin state.555Of course, this statement becomes invalid for material temperatures near absolute zero where the energy difference between different nuclear spin states becomes comparable to thermal energies. This independence means that the ensemble average of products of scattering lengths found at positions and obeys:

(13)

With this, Eq. 12 can be written as the sum of two distinct contributions:

(14)
(15)
(16)

The terms in the coherent scattering function directly depend on material structure as they involve pair correlations between nuclei at all positions in the material, whereas the incoherent scattering function solely contains self-correlations terms. The effective scattering length used for the nuclei at position in is equal to the statistical average of the scattering lengths at the position in the entire ensemble, and is referred to as the coherent scattering length, . The corresponding effective scattering length in is instead given as the statistical spread of the scattering lengths at the position in the entire ensemble, and is referred to as the incoherent scattering length, . Using , one might also talk about the related coherent and incoherent cross sections. Coherent and incoherent scattering lengths or cross sections are defined for each natural element (or any other well defined composition of isotopes), and can be directly applied to target systems in which a given position is always occupied by a nuclei of the same element in all subsystem replicas. It can be shown that for thermal neutrons, the sum of the incoherent and coherent scattering cross sections is to a good approximation equal to the unbound scattering cross section, , for which the total scattering cross section will converge at shorter wavelengths Sears (1984).666At very short wavelengths this conversion is ultimately disrupted by nuclear resonance physics or P-wave contributions.

The main difficulty in the evaluation of the scattering functions in Eqs. 15 and 16 is the integral over states in implicit in Eq. 10, whose evaluation in principle relates to the potentially complicated time-dependent distribution of nuclei in the target system. In a crystal, the nuclear positions can be decomposed in terms of displacements from equilibrium positions :

(17)

And thus:

(18)

Under the assumption that displacements are small compared to inter-atomic distances, motion of nuclei can be described with potentials which are quadratic functions of the displacements. In this so-called harmonic approximation, nuclei are essentially described as being pulled towards their equilibrium positions by linear spring-like forces, with resulting harmonic vibrations. Working in this approximation, it is possible to show that:

(19)

where the Debye-Waller function gives a measure of the time-independent average squared displacement of nuclei along :

(20)

The time dependence in Eq. 19 enters exclusively in the last exponential factor, which thus involves nuclear motion and can be expanded in a Taylor series as:

(21)

The expectation value correlates linear displacements along of two nuclei at different times. It is possible to show that this is directly related to an interaction in which the neutron exchanges energy and momentum with a phonon state, and the expansion in Eq. 21 thus lends itself to physics interpretation in the phonon picture, with the th term corresponding to -phonon interactions. At lower displacements or values the expansion converges more rapidly, and accordingly multi-phonon physics will be less important at lower neutron energies or material temperatures.

2.3 Elastic scattering

The first term in the expansion of Eq. 21 gives rise to elastic () scattering when inserted into Eq. 15 or Eq. 16, since:

(22)

With this, the partial differential cross section for incoherent elastic scattering can be immediately found:

(23)

The only directional dependency enters here in the Debye-Waller factor, , which approaches unity when the neutron wavelength is much larger than the atomic displacements, implying isotropic scattering. At higher neutron energies or temperatures this approximation breaks down, and the scattering will be increasingly focused in the forward direction at low values of the scatter angle , since for elastic scattering .

The coherent scattering terms involve correlations between pairs of nuclei and their evaluation will therefore be considerably more complex. In a crystal the sum over all nuclei can be rewritten as a sum over unit cells first and unit cell contents second: . Here the index assumes all values from Eq. 2 and the index runs over all nuclei in the unit cell. Thus, the equilibrium positions formerly denoted now becomes , with the positions defined in Eq. 1. Reordering sums appropriately and proceeding as for the incoherent case, the expression for coherent elastic scattering in crystals becomes:

(24)

Where the form factor of the unit cell have been introduced:

(25)

The parenthesised factor in Eq. 24 accounts for interference between different unit cells. Given that all for suitable choice of integers , and , the translational invariance of Eq. 3 implies:

(26)

Where is the number of unit cells in the crystal, which can be removed by adopting the convention of providing cross sections normalised per unit cell. The remaining factor can be investigated further by expressing in terms of the basis vectors of the reciprocal lattice from Eq. 4, :

(27)

Using:

(28)

Eq. 27 becomes:

(29)

Thus, coherent elastic scattering occurs only when the momentum transfer, , is exactly identical to one of the points in the reciprocal lattice, , corresponding to interaction with the associated family of lattice planes. The fact that further supports the interpretation of reflection by lattice planes with normal along , since elastic specular reflection by a mirror will always have the normal proportional to the transferred momentum. As constitutes an orthogonal but not an orthonormal base:

(30)

Where the last equality was found by inserting the definitions from Eq. 4 and evaluating the resulting cross product of two cross products with the rule. With this result, the coherent elastic scattering function (normalised to the unit cell) can be written more compactly:

(31)

The Debye-Waller factors will tend to suppress the form factors at high momenta, equivalent to large absolute values of , , and . The summation range of these indices can therefore in practice be limited to a region around 0. This will be explored further in Section 4.1 where the impact on total coherent elastic cross sections due to a lower bound, , on (corresponding to upper bounds on ) is investigated systematically for a large number of crystals. For the values selected for consideration it is then possible to pre-calculate the form factors, at the corresponding indices. As will be discussed in sections Sections 3 and 4, the initialisation of crystal structures in NCrystal indeed involves the preparation of lists of indices with corresponding -spacings and pre-calculated form factors. The Debye-Waller functions in the form factors are evaluated using the Debye model presented in Section 2.5.

The scattering described by Eq. 31 is usually referred to as Bragg diffraction. Reflections from the family of lattice planes indexed by requires . Since the scattering is elastic and where is the scattering angle. The maximal value of is found when and is . Thus, reflection is impossible unless (the Bragg condition):

(32)

And the scattering angle will satisfy the Bragg equation:

(33)

Often the above equation will be stated using the Bragg angle, defined as . It also often contains on the left hand side an integer factor, , denoting the so-called scattering order. However, this is merely a convenient manner in which to include multiple co-oriented lattice plane families in a single equation, utilising the fact that and .

For a crystal powder, in which crystal grains appear with uniformly randomised orientations, the total coherent elastic cross section for scattering on a given family of lattice planes satisfying Eq. 32, can be found as an isotropic average over the orientation of with respect to . Considering just the -functions, denoting the cosine of the angle between and with , and dropping the indices for brevity, this isotropic average gives:

(34)

Where it was used that . Inserting the omitted factors from Eq. 31 in Eq. 34, it is evident that the complete coherent elastic cross section of a crystal powder can be written as the following sum over all lattice plane families satisfying Eq. 32:

(35)

In case of a scattering event, the relative probability for it to happen on a particular plane will depend on its contribution to the sum in Eq. 35, and the scattering angle will subsequently be determined by Eq. 33. The azimuthal scattering angle around the direction of is, however, not constrained and all such angles contribute equally to the cross section. Thus, will be distributed uniformly in a cone around with opening angle . Such cones are known as Debye-Scherrer cones, and are a prominent feature at any powder diffractometer.

2.4 Inelastic scattering

The terms in the expansion of Eq. 21 can be interpreted as inelastic scatterings in which the incoming neutron exchanges energy and momentum with phonon states. The evaluation of these terms can in principle be very complicated depending on the material in question and the required precision. For the purposes of the present publication, inelastic processes play only a minor role, and the details of the implementations of such processes in NCrystal is reserved for a future dedicated publication, with just a brief overview provided here.

For non-oriented materials like liquids and crystal powders, rotational symmetry implies that the scatter functions become independent of the direction of around , and thus become two-dimensional: . As inelastic scattering terms do not include the factors of found in elastic terms, will generally be sufficiently smooth that it is possible to interpolate it from its values on a discrete grid in -space.777While phonon terms contain -functions ensuring momentum and energy conservation, the convolution with phonon state densities and integration over grain orientations in the powder approximation eliminates them from the final scatter functions. Such grids of tabulated values, usually referred to as scattering kernels in the context of Monte Carlo applications like MCNP, can in principle provide a high degree of detail, but their construction requires significant efforts. One appealing option is to directly measure in a neutron scattering experiment, which has obvious advantages but also disadvantages depending on experimental precision and coverage of space, as well as access to an appropriate neutron instrument in the first place. The other option is theoretical calculations, which usually must employ some approximations due to the complexity involved in an ab initio approach. One such powerful if computationally expensive numerical approach is molecular dynamics simulations in which semi-classical modelling of atomic trajectories are used to provide expectation values for atomic displacements and correlations Damián et al. (2014).

On the other end of the accuracy spectrum is the usage of various empirical closed form expressions Freund (1983); Cassels (1950); Binder (1970) describing the total inelastic cross section as a function of neutron energy or wavelength. Not only do such formulas not provide details on scatter angles or energy transfers, they usually require additional tuned parameters and tend to be valid either at very long or very short neutron wavelengths. Although it is possible to patch two such formulas together Kittelmann and Boin (2015); Adib and Fathalla (2007) in order to cover both ends of the wavelength spectrum, the resulting behaviour at intermediate wavelengths tends in general to be highly inaccurate.

Returning to Eq. 21, it is possible to employ various approximations in order to extract results. This approach is in particular useful at longer neutron wavelengths, where the single-phonon () term dominates, and indeed experimental techniques like neutron spectroscopy tend to focus on using the single-phonon term to access information about the dynamics of investigated samples, with contributions from multi-phonon terms () usually seen as unwanted background to be estimated and subtracted.

It is planned for NCrystal to support the modelling of inelastic scattering with enhanced data like scattering kernels or phonon density histograms when available, and prototype code with novel capabilities for precision sampling already exist Cai et al. (2018). It is nonetheless desirable that reasonably accurate modelling of inelastic scattering components should be provided by NCrystal even for materials where no more data than that required for Bragg diffraction is provided. Thus, an iterative procedure suggested in Sjölander (1958) has been adopted for NCrystal, in which the contribution to the scatter function involving phonons are estimated from the contribution involving phonons, starting from single-phonon contributions determined by the Debye model (cf. Section 2.5). Additionally working under the so-called incoherent approximation Placzek (1952), in which off-diagonal elements of are assumed to cancel each other out in coherent inelastic scattering, allows the prescribed method to estimate both coherent and incoherent components of inelastic scattering.

2.5 The Debye model

It is clear from the discussion so far that any actual evaluation of scattering functions or cross sections involves evaluation of the Debye-Waller function defined in Eq. 20, which is not surprising since this function captures the dynamics of the system due to thermal fluctuations. In principle the evaluation requires highly non-trivial material-specific information, either based on theory, numerical work, measurements, or a combination of those. In order to be able to provide meaningful results for user defined materials which might lack such specialised information, NCrystal code is currently evaluating the Debye-Waller factors using a simplified model, in which expected displacements are isotropic around the equilibrium positions and whose dynamics are otherwise governed by a model introduced by Debye Debye (1912). The Debye model assumes that phonons in crystals always propagate at a fixed velocity (neglecting certain effects like anisotropy and polarisation) and only exist below some frequency threshold, , with the associated Debye temperature, , corresponding to the temperature of the most energetic phonon. Loosely speaking, a high Debye temperature indicates a strongly bound material structure and vice versa. The model is additionally based on derivations by Glauber Glauber (1955), who applied the Debye model in order to estimate non-isotropic displacements.

With denoting a Cartesian coordinate and averaging over phonon polarisations, (Glauber, 1955, Eqs. 4-5) in the notation of the present paper implies:

(36)

Where the index runs over all phonon states, is the number of phonon states, is the atomic mass and is the material temperature. The sum over phonon states can be replaced by an integral with the state density, :

(37)

Now, as phonons in the Debye model propagate at constant velocity, the frequency of a phonon will be proportional to its momentum. Assuming that phonon states carry momenta distributed uniformly in momentum space, this implies . Fixing the normalisation yields . Inserting into Eq. 37 and performing a variable change gives:

(38)

Where:

(39)

A result which is identical to (Glauber, 1955, Eq. 11). It is interesting to note that the displacements predicted by Eq. 38 do not vanish at , which is a quantum mechanical effect related to the wave functions in the ground states. For a given atomic mass, temperature and Debye temperature, it is straight-forward to evaluate numerically and thus approximate the Debye-Waller function as:

(40)

where is the (isotropic) mean-squared displacement provided by Eq. 38. In principle the Debye temperature could be allowed to vary for each site in the unit cell, but in practice it is usual to allow it to be specified for each type of element found in the crystal and NCrystal accordingly supports the specification of Debye temperatures either globally or per type of element. It is often the case that Debye temperatures for a given crystal can be found in the literature (e.g. Peng et al. (1996); Wang et al. (2013)). If not, one might instead be able to obtain values for mean-squared displacements, which can then be used with Eq. 38 to estimate the Debye temperature, or one might alternatively return to the original purpose of the Debye model and determine the global Debye temperature for the material from its heat capacity.

Figure 1 shows the displacements predicted by the Debye model as a function of temperature for a number of crystals. Indeed, as required by the harmonic approximation, it is generally true that the predicted displacements are much smaller than typical inter-atomic distances, although the validity is strongest at lower temperatures.

Figure 1: Root-mean-squared atomic displacements in a number of crystals predicted by the Debye model as discussed in the text. The applied Debye temperatures are indicated in parentheses, and curves are terminated at the melting points where relevant.

3 Core framework overview

At its core, all capabilities of the NCrystal toolkit are implemented in an object oriented manner in a C++ library, providing both clearly defined interfaces for clients and internal separation between code implementing physics models, code loading data, and code providing infrastructure needed for integration into final applications. As will be discussed in further detail in Sections 5 and 6, it should be noted that while it is certainly possible to use the C++ classes discussed in the present section directly, typical users are expected to use other more suitable interfaces for their work – employing also a simpler and generic approach to material configuration.

In Figure 2 is shown the most important classes available in the release of NCrystal presented here, along with their internal inheritance relationships. At the root of the tree sits a few infrastructure classes not directly related to actual physics modelling, starting with the universal base class RCBase, which provides reference counted memory management for all derived classes.888Although C++11 provides alternatives for such reference counting in the form of modern smart pointers, it is for the time being the aim of NCrystal to support also C++98, in which support for such is incomplete. One level below this sits the CalcBase class, from which all classes actually implementing physics modelling code must inherit. The main feature of CalcBase is currently that it provides derived classes with access to pseudo-random numbers, in a manner which can be configured either globally or separately for each CalcBase instance. In order to do so, each CalcBase instance keeps a reference to an instance of a class deriving from RandomBase – a class which provides a generic interface for pseudo-random number generation. The primary purpose of this setup is to let NCrystal code use external sources of random numbers when embedded into existing frameworks managing their own random number generators, like those discussed in Sections 6.3 and 6.4. If no particular random generator is otherwise enabled, the system will fall back to using a xoroshiro128+ Blackman and Vigna (2018) generator implemented in the RandXRSR class.

Figure 2: NCrystal class hierarchy. Abstract classes are indicated with rounded edges while concrete classes are shown in octagons. The most important interfaces providing end-user results are indicated in red and final classes providing physics modelling are shown in blue. Internal utility and factory classes are not shown.

Currently, the only class deriving directly from CalcBase is Process. This is an abstract class representing a physics process, specifying a general interface for calculation of cross sections as a function of incident neutron state. The returned cross section values are normalised to the number of atoms in the sample, and the neutron states are specified in terms of kinetic energy () and direction (), which is equivalent to providing but using parameters typically available in general purpose Monte Carlo applications without conversions. Additionally, for reasons of convenience and computational efficiency, it is possible for processes to indicate if they are only relevant for a given energy range (domain), and processes are divided into those that are oriented in the sense that they actually depend on and those (isotropic) ones that do not. For the latter, one can access the cross section without specifying a direction. For reference, the precise methods available via the Process interface can be seen in Table 1.

Methods available via the Process interface:
bool isOriented() const;
double crossSection(double ekin, const double (&indir)[3] ) const;
double crossSectionNonOriented( double ekin ) const;
void domain(double& ekin_low, double& ekin_high) const;
Additional methods available via the Scatter interface:
void generateScattering( double ekin, const double (&indir)[3],
\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ double (&outdir)[3], double& delta_ekin ) const;
void generateScatteringNonOriented( double ekin, double& angle,
\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ double& delta_ekin ) const;
Table 1: C++ methods provided on the Process and Scatter interface classes.

Two abstract classes further specialise the Process interface, with roles indicated by their names: Absorption and Scatter. The Absorption class does not currently provide any functionality over the Process class, since the current scope of NCrystal does not include any particular description of absorption reactions apart from their cross sections. In principle it would of course be possibly to extend this class in the future, should the need ever arise for NCrystal to be able to model the actual outcome of such events in terms of secondary particles produced.

The Scatter class does on the other hand extend the Process interface, adding methods for random sampling of final states, as can also be seen in Table 1. Specifically, this sampling generates both energy transfers, , and final direction of the neutron, . In the case of isotropic processes, only the energy transfer and the polar scattering angle, , from to will depend on the actual physics implemented, and will itself be independent of . The azimuthal scattering angle will be independent and uniformly distributed in . For isotropic processes, it is thus again possible if desired to avoid the methods with full directional vectors and instead sample energy transfer and polar scatter angle given just the kinetic energy. The specialised base class ScatterIsotropic is provided for developers in order to simplify implementation of isotropic scatter processes. This somewhat complicated design was chosen in order to make it straight-forward to accommodate both oriented and isotropic scatter processes in the same manner, avoiding unnecessary duplication of code and interfaces while still making it possible to retain full computational efficiency allowed by the symmetries in the isotropic cases. The next class in Figure 2, ScatterXSCurve, exists in order to facilitate scatter models which only provide cross sections, adding simple fall-back models for sampling of final states. Rounding up the scattering infrastructure, ScatterComp is a wrapper class playing a special role, representing the composition of multiple Scatter objects into a single one. It is used whenever multiple independent components representing partial cross sections must be combined in order to fully describe the scattering process of interest – such as when a Bragg diffraction component is combined with a component providing inelastic or incoherent scattering. At the bottom of the class hierarchy in Figure 2 sits classes actually implementing physics models of absorption or scattering, shown in blue. They will be discussed in further detail in Section 3.1.

Finally, the Info class is a data structure containing information about crystal structure at the microscopic scale of crystallites. As illustrated in Figure 3, the Info class serves the role of separating sources of crystal definitions from the physics models using the data as input. The data fields available on the Info class are provided in Table 2. Not all input sources will be able to provide all the data fields shown, nor will a particular physics modelling algorithm require all fields to be available. Thus, to ensure maximal flexibility for both providers and consumers of Info objects, all data fields are considered optional. If a given physics algorithm does not find the information it needs, it should simply indicate a configuration error, by throwing an appropriate C++ exception. When configuring materials via the recommended method for end-users (cf. Section 5), the factory code will generally try to avoid triggering undesired errors when the user intent seems obvious: if for instance the input data does not provide any information out of which it would be possible to provide cross sections for inelastic scattering processes, it is most likely that the user is simply only interested in elastic physics and no instantiation of inelastic processes will be attempted.

Figure 3: Flow of crystal data in NCrystal. Different factories are responsible for turning crystal definitions from a variety of sources into full-fledged instantiations of the Info class, which is then consumed by various physics models.
Field Contents Units
StructureInfo Basic info about unit cell:
   Lattice parameters (, , , , , ) ,
   Volume
   Number of atoms
   Space group number
AtomInfo (list) List of atoms in unit cell, each with:
   Atomic number (Z)
   Number per unit cell
   Per-element Debye temperature
   Mean squared displacement
   List of element positions
HKLInfo (list) List of families, each with:
    value (representative)
   -spacing
   Form factor,
   Multiplicity (family size)
   List of all values in family
   List of all normals in family
Threshold -spacing value for list
Absorption cross section at
Unbound scattering cross section
Cross section curve for inelastic/incoh.
scattering components (function object)
Temperature Material temperature
Debye temperature Global Debye temperature of material
Density Material density
Table 2: Data fields of Info objects. A field can represent either a scalar value or a block of related data, and the availability of all fields is optional. Additionally, some parameters inside block fields might themselves be optional, as indicated with ’s. Cross sections are given as per-atom values, averaged over the atoms of the unit cell. An “ family” is here a group of indices for which form factors and -spacings have identical values.

This scheme of decoupling data loading code from physics models allows NCrystal to easily support data read from a variety of input files, and makes it simple to add support for additional input sources in the future. If desired, it would even be possible to specify the information directly in program code or load it from a database. Currently supported data sources are discussed in Section 4.

In Listing 1 is shown an example of C++ code loading a crystal definition and constructing objects able to model absorption and scattering in a polycrystalline or powdered material. The example illustrates in practice how some of the classes discussed in this section are used, but also indicates the complexity arising out of the need to provide model-specific parameters. In Section 5, a more convenient method for initialisation and configuration will be introduced.

#include "NCrystal/NCrystal.hh"
namespace NC = NCrystal;//alias for convenience
int main() {
  //Create crystal Info based on a datafile and other
  //parameters (here temperature and d-spacing cutoff):
  std::string datafile = "somefile.ncmat";
  double temp_kelvin = 400;
  double dcutoff_angstrom = 0.5;
  const NC::Info * info = NC::loadNCMAT( datafile,
                                         temp_kelvin,
                                         dcutoff_angstrom );
  //Compose Scatter and Absorption objects, passing the Info
  //object and other parameters as needed:
  bool bkgd_thermalise = true;
  NC::ScatterComp * scat = new NC::ScatterComp();
  scat->addComponent(new NC::PCBragg(info));
  scat->addComponent(new NC::BkgdPhonDebye( info,
                                            bkgd_thermalise ));
  NC::Absorption * absn = new NC::AbsOOV(info);
  // ...
  // code here using info, scat and absn pointers
  // ...
  return 0;
}
Listing 1: C++ code loading a crystal definition and creating related NCrystal objects.

3.1 Physics models

At the heart of NCrystal is the actual modelling of specific physics processes, provided via the classes shown with blue outlines in Figure 2. The only modelling of absorption cross sections presently available is provided in the AbsOOV class, which implements the scaling model discussed in Section 2.2, by scaling the value of given at the reference velocity of (cf. Table 2). Despite being implemented with a simple model, absorption cross sections are nonetheless provided in NCrystal for completeness, in order to carry out validation against data from measurements of total cross sections like those presented in Section 4.4 or to facilitate the creation of plugins for Monte Carlo simulations like the one presented in Section 6.4. Additionally, it is important to reiterate that for most nuclear isotopes, the modelling is actually remarkably accurate at thermal neutron energies (cf. Section 2.2).

The remaining physics models all concern scattering processes. Three processes implement the coherent elastic physics of Bragg diffraction discussed in Section 2.3, differing by which distribution of crystallite grains they assume as discussed in Section 2.1. PCBragg implements an ideal powder model, with cross sections given by Eq. 35 and scattering into Debye-Scherrer cones. It can be used to model powders and polycrystalline materials. SCBragg implements a model for single crystals with isotropic Gaussian mosaicity, and LCBragg provides a special anisotropic model of layered crystals similar to pyrolytic graphite. For reasons of scope, the three models of Bragg diffraction will be presented in detail in a future dedicated publication, covering details of both theory, implementation and validation.

Two models currently provide incoherent and inelastic physics: BkgdExtCurve and BkgdPhonDebye.999The term “Bkgd” is here used as a short-hand for inelastic and incoherent processes. It is merely used to lighten the notation, and of course is not meant to imply that e.g. inelastic neutron scattering is always to be considered “background” to some other signal (which would be incorrect). The former is a simple wrapper of externally provided cross section curves (cf. Table 2), and exists mostly for reference, allowing nxslib-provided cross section curves to be exposed via NCrystal when using .nxs files (cf. Section 4.2). Unless selected explicitly, the factories discussed in Section 5 will always prefer the more accurate model provided by BkgdPhonDebye. This improved model relies on the incoherent approximation and the iterative procedure discussed in Section 2.4, and will be presented in detail in a future dedicated publication. In the current release of NCrystal it does not provide any detailed sampling of scatter angles or energy transfers, and hence like BkgdExtCurve it derives from the ScatterXSCurve class. When asked to generate a scattering, ScatterXSCurve will scatter isotropically and either model energy transfers as absent (i.e. elastic), or as “fully thermalising”. In the latter case, the energy of the outgoing neutron will be sampled from a thermal (Maxwell) energy distribution specified solely by the temperature of the material simulated. Both of these choices use well defined distributions, but are clearly lacking in realism. It is the plan to improve this situation in the future, in order to provide a complete and consistent treatment of all components involved in thermal neutron scattering. First and foremost by adopting proper sampling models in BkgdPhonDebye, still based on the currently used iterative procedure for cross section estimation. As already discussed in Section 2.4, it is additionally planned to allow for even higher accuracy and reliability when modelling inelastic scattering in select materials, by introducing data-driven models which are able to utilise pre-tabulated scatter-kernels or phonon state densities, if available.

4 Crystal data initialisation

Currently, NCrystal supports the loading of crystal information from four different file formats: the native and recommended .ncmat format which will be introduced in Section 4.1, and the existing .nxs, .laz, and .lau formats which will be discussed in Sections 4.2 and 4.3. The latter formats are mostly supported as a service to the neutron scattering community already using these, for instance in the context of configuring instrument simulations in McStas Lefmann and Nielsen (1999); Willendrup et al. (2004) or VitESS Zsigmond et al. (2002); Lieutenant et al. (2004). Support for direct loading of CIF (“Crystallographic Information File”) files Hall et al. (1991) was considered as they are readily available in online databases Gražulis et al. (2009); Downs and Hall-Wallace (2003), but ultimately abandoned due to the flexible nature of such files, many of which were found in practice to not contain suitable information for the purposes of NCrystal. Manually selected CIF files were, however, read with Larsen et al. (2017) and used to assemble a library of .ncmat files which is shipped with NCrystal. These files, providing the crystal structures listed in Table 3, first and foremost describe many materials of interest to potential NCrystal users, serving both as a convenient starting point and point of reference. Moreover, the library encompasses six of the seven crystal systems, includes mono- and poly-atomic crystals, and includes materials with large variations in quantities like Debye temperatures, neutron scattering lengths and number of atoms per unit cell. Thus, it provides a convenient ensemble for benchmarking and validating NCrystal code. Additionally, where not prevented for technical reasons like usage of per-element Debye temperatures in poly-atomic systems, .nxs files were automatically generated from the provided .ncmat files. Section 4.4 will describe the validation carried out as concerns both the data files themselves, as well as the code responsible for loading crystal information from them. Naturally, it is fully expected that the list of files in Table 3 will be expanded in the future, in response to requests from the user community.

\ssmall
Crystal Space group References Formats Validations
Ag 225 (Cubic) (Downs and Hall-Wallace, 2003, 11135)Wyckoff (1963)Peng et al. (1996) .ncmat, .nxs N
AlO (Corundum) 167 (Trigonal) (Downs and Hall-Wallace, 2003, 09327)Kirfel and Eichhorn (1990) .ncmat GF
Al 225 (Cubic) (Downs and Hall-Wallace, 2003, 11136)Wyckoff (1963)Peng et al. (1996) .ncmat, .nxs NT
Au 225 (Cubic) (Downs and Hall-Wallace, 2003, 11140)Wyckoff (1963)Peng et al. (1996) .ncmat, .nxs N
Ba 229 (Cubic) (Downs and Hall-Wallace, 2003, 11207)Wyckoff (1963)Peng et al. (1996) .ncmat, .nxs NF
BeF (Be. fluoride) 152 (Trigonal) Ghalsasi and Ghalsasi (2011) .ncmat F
Be 194 (Hexagonal) (Downs and Hall-Wallace, 2003, 11165)Wyckoff (1963)Peng et al. (1996) .ncmat, .nxs NT
C (Pyr. graphite) 194 (Hexagonal) (Downs and Hall-Wallace, 2003, 14675)Trucano and Chen (1975) .ncmat NTF
C (Diamond) 227 (Cubic) (Downs and Hall-Wallace, 2003, 11242)Wyckoff (1963)Peng et al. (1996) .ncmat, .nxs NF
CaCO (Aragonite) 62 (Orthorhombic) (Downs and Hall-Wallace, 2003, 06300)Antao and Hassan (2009) .ncmat GF
Ca 225 (Cubic) (Downs and Hall-Wallace, 2003, 11141)Wyckoff (1963)Peng et al. (1996) .ncmat, .nxs N
Ca (-calcium) 229 (Cubic) (Downs and Hall-Wallace, 2003, 11208)Wyckoff (1963)Peng et al. (1996) .ncmat, .nxs NF
Cr 229 (Cubic) (Downs and Hall-Wallace, 2003, 11209)Wyckoff (1963)Peng et al. (1996) .ncmat, .nxs NT
CuO (Cuprite) 224 (Cubic) (Downs and Hall-Wallace, 2003, 09326)Kirfel and Eichhorn (1990) .ncmat F
Cu 225 (Cubic) (Downs and Hall-Wallace, 2003, 11145)Wyckoff (1963)Peng et al. (1996) .ncmat, .nxs NT
Fe (-iron) 229 (Cubic) (Downs and Hall-Wallace, 2003, 11214)Wyckoff (1963)Peng et al. (1996) .ncmat, .nxs N
Fe (-iron) 229 (Cubic) (Downs and Hall-Wallace, 2003, 11215)Wyckoff (1963)Peng et al. (1996) .ncmat, .nxs NF
Ge 227 (Cubic) (Downs and Hall-Wallace, 2003, 11245)Wyckoff (1963)Peng et al. (1996) .ncmat, .nxs NT
MgO (Periclase) 225 (Cubic) (Downs and Hall-Wallace, 2003, 00501)Hazen (1976) .ncmat F
Mg 194 (Hexagonal) (Downs and Hall-Wallace, 2003, 11183)Wyckoff (1963)Peng et al. (1996) .ncmat, .nxs N
Mo 229 (Cubic) (Downs and Hall-Wallace, 2003, 11221)Wyckoff (1963)Peng et al. (1996) .ncmat, .nxs NT
Na 229 (Cubic) (Downs and Hall-Wallace, 2003, 11223)Wyckoff (1963)Peng et al. (1996) .ncmat, .nxs NF
NaSiAlOCl (Sodalite) 218 (Cubic) (Downs and Hall-Wallace, 2003, 06211)Antao et al. (2008) .ncmat F
Nb 229 (Cubic) (Downs and Hall-Wallace, 2003, 11224)Wyckoff (1963)Peng et al. (1996) .ncmat, .nxs NT
Ni 225 (Cubic) (Downs and Hall-Wallace, 2003, 11153)Wyckoff (1963)Peng et al. (1996) .ncmat, .nxs NT
Pb 225 (Cubic) (Downs and Hall-Wallace, 2003, 11154)Wyckoff (1963)Peng et al. (1996) .ncmat, .nxs N
Pd 225 (Cubic) (Downs and Hall-Wallace, 2003, 11155)Wyckoff (1963)Peng et al. (1996) .ncmat, .nxs N
Pt 225 (Cubic) (Downs and Hall-Wallace, 2003, 11157)Wyckoff (1963)Peng et al. (1996) .ncmat, .nxs N
Rb 229 (Cubic) (Downs and Hall-Wallace, 2003, 11228)Wyckoff (1963)Peng et al. (1996) .ncmat, .nxs NF
Sc 194 (Hexagonal) (Downs and Hall-Wallace, 2003, 11192)Wyckoff (1963)Peng et al. (1996) .ncmat, .nxs N
SiLuO 15 (Monoclinic) Gustafsson et al. (2001) .ncmat RF
SiO (Quartz) 154 (Trigonal) (Downs and Hall-Wallace, 2003, 06212)Antao et al. (2008) .ncmat GF
Si 227 (Cubic) (Downs and Hall-Wallace, 2003, 11243)Wyckoff (1963)Peng et al. (1996) .ncmat, .nxs NT
Sn 141 (Tetragonal) (Downs and Hall-Wallace, 2003, 11248)Wyckoff (1963)Peng et al. (1996) .ncmat, .nxs NT
Sr 225 (Cubic) (Downs and Hall-Wallace, 2003, 11161)Wyckoff (1963)Peng et al. (1996) .ncmat, .nxs NF
Ti 194 (Hexagonal) (Downs and Hall-Wallace, 2003, 11195)Wyckoff (1963)Peng et al. (1996) .ncmat, .nxs N
UO (Uraninite) 225 (Cubic) (Downs and Hall-Wallace, 2003, 11728)Wyckoff (1963)Wang et al. (2013) .ncmat, .nxs NT
V 229 (Cubic) (Downs and Hall-Wallace, 2003, 11235)Wyckoff (1963)Peng et al. (1996) .ncmat, .nxs N
W 229 (Cubic) (Downs and Hall-Wallace, 2003, 11236)Wyckoff (1963)Peng et al. (1996) .ncmat, .nxs N
YO (Ytr. oxide) 206 (Cubic) Park et al. (2011) .ncmat GF
Y 194 (Hexagonal) (Downs and Hall-Wallace, 2003, 11199)Wyckoff (1963)Peng et al. (1996) .ncmat, .nxs NT
Zn 194 (Hexagonal) (Downs and Hall-Wallace, 2003, 11200)Wyckoff (1963)Peng et al. (1996) .ncmat, .nxs NT
Zr 194 (Hexagonal) (Downs and Hall-Wallace, 2003, 11201)Wyckoff (1963)Peng et al. (1996) .ncmat, .nxs N
Table 3: Crystal definitions available in the data library currently shipped with NCrystal. Crystal structures and Debye temperature data were obtained from the listed references, with specific database IDs indicated where structures were taken from Downs and Hall-Wallace (2003). The last column indicates the validations performed, as discussed in Section 4.4.

4.1 NCrystal material files (.ncmat)

Keeping in mind both human readers and parsing algorithms, the layout of NCrystal material files is deliberately kept simple, as illustrated with the sample file in Listing 2. These simple text files always begin with the keyword NCMAT and the version of the format. Next comes optionally a number of #-prefixed lines with free-form comments, intended primarily as a place to document the origin, purpose or suitability of the file. The data itself follows after this introduction, and is placed in four clearly denoted sections as shown in the listing. The @CELL and @ATOMPOSITIONS sections define the unit cell layout directly, and the definition must be compatible with the space group number which can optionally be provided in the @SPACEGROUP section. Finally, the @DEBYETEMPERATURE section contains the effective Debye temperature of the crystal (cf. Section 2.5) – either a single global value, or as per-element values, indicated with the name of each element (for mono-atomic crystals, there is of course no difference between specifying a per-element or a global value). Where relevant, values are specified in units of , and for lengths, temperature and angles respectively.

NCMAT v1
#Converted from the CIF file of the entry 0011136 in the AMCSD
#reference: Wyckoff R W G, Crystal Structures, vol. 1, p. 7-83,
#1963. The Debye temperature is derived from the Debye-Waller
#factor at 293K compiled in the supplement of Acta Cryst., A52,
#p. 456-470, 1996.
@CELL
    lengths 4.04958 4.04958 4.04958
    angles 90. 90. 90.
@SPACEGROUP
    225
@ATOMPOSITIONS
    Al 0. 0.5 0.5
    Al 0. 0. 0.
    Al 0.5 0.5 0.
    Al 0.5 0. 0.5
@DEBYETEMPERATURE
    Al 410.35
Listing 2: Sample .ncmat file defining the face-centred cubic structure of a pure aluminium crystal, with lattice lengths of and a Debye temperature of .

The atomic positions are specified using relative lattice coordinates (cf. Eq. 1), and the name of an element is also required at each such position. It is not necessary, nor currently possible, to specify element-specific data such as masses, cross sections or scattering lengths. Instead, NCrystal currently ships with an embedded database of such numbers for all natural elements with , based on Jericha and Hainbuchner (2001); Rauch and Waschkowski (2000); Koester et al. (1991). Although somewhat inflexible, this scheme provides a high level of convenience, consistency and robustness by significantly lowering the amount of parameters required in each .ncmat file. It is envisioned that a future version of NCrystal would support more specialised use-cases by supporting not only the optional specification of element- and isotope-specific data in .ncmat files, but also to allow for specification of features like chemical disorder, impurities, doping or enrichment.

With the information parsed directly from the .ncmat file and the specification of a material temperature, the remaining information shown in Table 2 will be derived numerically at initialisation time (with the exception of