Thermal neutron scattering law calculations using ab initio molecular dynamics

In recent years, methods for the calculation of the thermal scattering law (i.e. S(α, β), where α and β are dimensionless momentum and energy transfer variables, respectively) were developed based on ab initio lattice dynamics (AILD) and/or classical molecular dynamics (CMD). While these methods are now mature and efficient, further advancement in the application of such atomistic techniques is possible using ab initio molecular dynamics (AIMD) methods. In this case, temperature effects are inherently included in the calculation, e.g. phonon density of states (DOS), while using ab initio force fields that eliminate the need for parameterized semi-empirical force fields. In this work, AIMD simulations were performed to predict the phonon spectra as a function of temperature for beryllium and graphite, which are representative nuclear reactor moderator and reflector materials. Subsequently, the calculated phonon spectra were utilized to predict S(α, β) using the LEAPR module of the NJOY code. The AIMD models of beryllium and graphite were 5 × 5 × 5 crystal unit cells (250 atoms and 500 atoms respectively). Electronic structure calculations for the prediction of Hellman-Feynman forces were performed using density functional theory with a GGA exchange correlation functional and corresponding core electron pseudopotentials. AIMD simulations of 1000– 10,000 time-steps were performed with the canonical ensemble (NVT thermostat) for several temperatures between 300 K and 900 K. The phonon DOS were calculated as the power spectrum of the AIMD predicted velocity autocorrelation functions. The resulting AIMD phonon DOS and corresponding inelastic thermal neutron scattering cross sections at 300 K, where anharmonic effects are expected to be small, were found to be in reasonable agreement with the results generated using traditional AILD. This illustrated the validity of the AIMD approach. However, since the impact of the temperature on the phonon DOS (e.g. broadening of spectral peaks) was observed in AIMD analysis, this technique may be envisioned as the approach for deriving the needed atomistic data for thermal scattering law calculations under realistic temperature and structural conditions for a given material.


Introduction and theory
Thermal neutrons are characterized as low energy electrons with energies on the order of the excitations (e.g. vibrations, rotations) of the medium in which they interact. Moreover these neutrons have De Broglie wavelengths on the order of the distance between atoms and are therefore highly influenced by the atomic structure and dynamics. Consequently, the structure and dynamics of atomic systems may be examined through the use of thermal neutron scattering. The thermal neutron scattering cross-section is typically defined within the first Born-Approximation [1]. In this case, the double-differential scattering cross-section is dependent on the thermal scattering law (TSL), S(α,β) (i.e. dynamic structure factor), as, (1) The double differential cross-section describes the probability of a thermal neutron of incident energy E scattering to energy E through solid angle . Momentum and energy transfer are expressed in terms of the a e-mail: ayman.hawari@ncsu.edu dimensionless variables α and β respectively [2]. The probability of scattering with an atomic nucleus is defined by the coherent and incoherent nuclear scattering crosssections (σ coh and σ inc ), whereas the outcome is defined by the TSL. The scattering law has contributions from noninterference effects (self-part, S s ) and inter-atomic interference effects (distinct-part, S d ), Although the self-part and distinct part both describe elastic and inelastic scattering, the contribution the distinct part to the total inelastic scattering cross-section is typically small. As a consequence of this disparity the TSL are frequently computed using the incoherent approximate, where-by atoms are assumed to move (e.g. vibrate in solid) relatively independent of each other.
If the atomic motion is assumed to be well described by harmonic vibrations, then the self-part of the TSL may be decomposed into an expansion describing increasing orders of vibrational excitation. In crystalline materials this decomposition is the phonon expansion,

ND2016
and requires the phonon density of states (DOS), ρ(β), as a fundamental input. The super-script, P, indicates the number of phonons emitted or absorbed during a single scattering event [1][2][3][4]. Each order of the expansion is a temperature dependent functional of the DOS, and may be computed as a convolution of the DOS with the previous member of phonon order P-1.

Methods of generating the thermal scattering law
The generation of TSL requires the phonon DOS as the fundamental input, which may be calculated from atomistic models. In particular ab initio lattice dynamics (AILD) and classical molecular dynamics (CMD) represent predictive approaches that may be utilized to calculate the phonon DOS for generation of the scattering law in crystalline materials [3]. These methods are mature and efficient and have been utilized in the generation of ENDF TSL libraries [5,6]; however, these methods may be further advanced with the use of ab initio molecular dynamics (AIMD). Modern computational power and techniques have recently made AIMD accessible for calculations of the phonon DOS. This method is an emerging technique for the calculation of phonon DOS for use in the generation of the TSL. In this work the phonon DOS was generated as a function of temperature for beryllium and graphite using AIMD methods, and was compared to established ab initio lattice dynamics methods (AILD) [3,4]. Subsequently, the TSL and inelastic scattering cross-section were evaluated for both methods.

Ab initio molecular dynamics
Ab initio molecular dynamics represents a combination of the AILD and CMD techniques. In AIMD Hellmann-Feynman forces derived from electronic structure calculations are utilized to evolve the atoms. As in CMD, this AIMD method uses classical mechanics to evolve an atomic system; however, the forces are not semiempirical but are in fact predictive. Equilibrium dynamics of an atomic system may be evolved using Car-Parrinello method, in which the Born-Oppenheimer approximation is invoked such that the dynamics of the simulation are defined on the time scale of atomic motion [7]. In these simulations, the atomic motion is dominantly harmonic such that the phonon DOS may be calculated as a Fourier transform of the velocity auto-correlation function (VACF), where ω is the phonon frequency (i.e. energy, E = ω) and v j (t) are the atom velocities at a time t [8].

Calculation of phonon density of states
Beryllium metal and graphite were selected as an exemplar moderator and reflector materials in this work to test the capability of AIMD to predict the phonon DOS. Beryllium is a simple HCP metal (space group 192), whereas graphite is a semi-metal (space group 186), see Fig. 1. Due to the small number of core and valence electrons (beryllium 1s 2 and 2s 2 , and graphite 1s 2 and 2s 2 2p 2 ) these materials are well suited to modeling using modern density functional theory methods with atomic-core pseudopotentials. To explore the phonon behavior of these materials AIMD and AILD simulations were performed using the VASP density functional theory (DFT) code [9,10] with PAW-GGA pseudo-potential libraries [11][12][13]. AIMD simulations of phonons in beryllium and graphite were performed using 5 × 5 × 5 supercells (250 atoms and 500 atoms respectively), as illustrated in Fig. 1. Before performing AIMD simulations DFT simulations on the unit-cell were performed to predict the 0 K lattice parameters, summarized in Tables 1  and 2 for beryllium and graphite respectively. Subsequently, the lattice parameters of the supercell at each temperature were determined from these 0 K parameters using the experimental thermal expansion coefficients [14,15]. Initially the supercell was equilibrated using AIMD at the simulation temperature for 1000 time-steps of 1fs using a canonical thermodynamic ensemble (NVT). To reduce the limitations that periodic boundary conditions impose on the permitted phonon wavelengths. a Langevin thermostat was added to all atoms adjacent to the boundary of the supercell [16]. Following equilibration, each temperature was simulated for an additional 10000 time-steps for beryllium and 1000 time-steps for graphite to generate the evolution of the atom velocities needed for the calculation of the phonon density of states. These simulation lengths were found capture the oscillations in the VACF necessary for the calculation of the phonon density of states.

Beryllium phonon DOS
Phonons were examined for beryllium using AIMD simulations at 300 K, 600 K, and 900 K. In these simulations a 1 × 1 × 1 Monkhorst Pack k-mesh and 350 eV plane wave cutoff were used. The phonon DOS were calculated at each temperature from the VACF, Eq. (4), and compared in Fig. 2. As the temperature increases the change in the lattice parameter results in weaker Table 1. Beryllium lattice parameters obtained through DFT calculations (calc.) at 0 K compared to measurement. The orientation of the a and c parameters may be observed in Fig. 1  interatomic forces and a corresponding shift in the phonon DOS toward lower energies. Moreover, broadening of peaks in the DOS occur due to increasing contributions of anharmonic forces with increasing temperature. The fluctuations in the phonon DOS at energies below 0.04 eV may be attributed to periodic conditions of the MD simulation [16]. Periodicity in the simulation limits the maximum phonon wavelength and thus the minimum phonon frequency. Consequently, the number of permitted wavelengths for low energy phonons is finite and increasingly sparse. Although the Langevin bath in the current AIMD simulations reduces this effect, the effect may be further reduced with larger system sizes. In the present work Savitzky-Golay smoothing was used to filter noise from the DOS [17].
The phonon spectrum of beryllium was also calculated using AILD, which is the current method by which phonon DOS for generation of TSL are calculated. The VASP code was used to generate Hellmann-Feynman forces on a 4 × 4 × 3 supercell with a 4 × 4 × 4 k-mesh and 350 eV plane wave cutoff. These Hellmann-Feynman forces were supplied to the lattice dynamics code PHONON [19]. The phonon dispersion relationship predicted using this method, Fig. 3, demonstrates reasonable agreement with the phonons measured using neutron scattering [20,21]. At 300 K the impact of anharmonic forces in beryllium is expected to be small [22]. The AIMD phonon DOS at this temperature is, therefore, compared to AILD (see Fig. 4). Phonon energies below 0.04 eV are expected to have greatest impact on the inelastic thermal neutron crosssection. In this region the AIMD phonon spectrum is in reasonable agreement to that of AILD. Considering the impact of thermal expansion, good agreement between both methods is observed over-all.

Graphite phonon DOS
AIMD simulations of phonons in graphite were performed at 300 K, 600 K, and 900 K with a 1 × 1 × 1 Monkhorst Pack k-mesh and 500 eV plane wave cutoff. In graphite the weak Van der Waals (VDW) forces which exist between each plane are not captured by traditional DFT. This force was included in the current model by adding the VDW interaction between pairs of graphite atoms using semi-empirical parameters [23]. The VDW radius of carbon atoms was parameterized for graphite in this work to predict the c-lattice parameter within reasonable agreement of experiment using the GGA-PBE pseudopotential, which significantly over-predicts this lattice parameter in the absences of VDW. The VDW radius and the corresponding lattice parameters at 0 K are tabulated in Table 2 and compared to theory and measurements.
The phonon DOS calculated at 300 K, 600 K and 900 K calculated from the VACF, Eq. (4), are compared in Fig. 5. As in beryllium, the change in lattice parameter with increasing temperature in graphite results in a shifting of the phonon DOS toward lower energies. Furthermore, the increasing contributions of anharmonic forces result in a broadening of phonon peaks (e.g. 0.02 eV, 0.6 eV, and 0.017-0.02 eV). The low energy region, below 0.02 eV, corresponds to vibrations along the c-axis. The phonon spectrum of graphite was calculated using AILD. The AILD phonons were calculated using the VASP and PHONON codes in a 4 × 4 × 4 supercell with a 4 × 4 × 2 k-mesh and 650 eV plane wave cutoff. As in the AIMD method, the AILD calculations included the additional VDW interactions. The phonon dispersion relationship predicted using this method, Fig. 6, demonstrates good agreement with the phonons measured using neutron and Raman scattering [26][27][28][29][30][31]. The AIMD phonon DOS at 300 K is compared to AILD in Fig. 7. As demonstrated in the dispersion curve, the phonons at small wave-vectors near the -point have energies less than 0.030 eV. The phonon DOS in this region is expected to be sensitive to the AIMD system size and periodic boundary conditions in the basal plane (i.e. plane composed of a-axis and b-axis) where the forces are approximately harmonic [16]. The contributions from phonons propagating along the c-axis ( -A) result from anharmonic forces between graphite layers. Consequently, reasonable agreement in the predicted phonon DOS is found between AIMD and AILD below 0.01 eV, which may be expected to most significantly affect the inelastic scattering cross-section in graphite. Nevertheless, the phonon DOS above this low energy region have similar characteristics.

Generation of thermal scattering law
The TSL for beryllium and graphite were evaluated at each temperature using the LEAPR module of the NJOY nuclear data processing package [2,33]. This module utilizes the phonon expansion, Eq. (3), to generate the TSL with the phonon DOS as the fundamental input. To examine the impact of the temperature dependent phonon DOS from AIMD, the inelastic TSL was generated at each simulation temperature using the corresponding DOS  and compared to the TSL using the AILD phonon DOS. Specifically, to demonstrate the impact of the temperature dependent AIMD phonon DOS, the inelastic cross-section and secondary neutron energy spectrum were calculated from TSL using the THERMR module of NJOY and compared to AILD. The inelastic scattering cross-section is an indicator of the rate at which neutrons thermalize, whereas the secondary spectrum describes the probability distribution of energy loss during these scattering events.

Beryllium TSL
The inelastic scattering cross-sections and secondary spectrum in beryllium are compared for the AIMD and AILD methods in Figs. 8 and 9, respectively. As temperature increases, the AIMD simulation predicts an increasingly higher contribution to inelastic scattering. This is consistent with the temperature dependent broadening and shifting of the phonon DOS, which result from an increase in the contribution of anharmonicity to atomic motion.
The secondary spectrum is compared at 300 K. In this case, the contribution of temperature effects to the phonon DOS (e.g. anharmonicity) are expected to be  small [22]; therefore, this temperature provides a basis for evaluation of the AIMD method in comparison to the established AILD method. Reasonable agreement is observed between the AIMD and AILD methods at 300 K in both the inelastic scattering cross-section and secondary spectrum. At energies between 0.1 and 5 eV both methods predict approximately equal secondary distributions, and may be expected to result in comparable predictions of neutron thermalization. For neutrons with incident energies near thermal energy (0.0253 eV) the methods differ; however, these deviation occur for secondary energies with probabilities an order of magnitude less than the most probable.

Graphite TSL
The inelastic scattering cross-sections and secondary spectrum in graphite are compared for the AIMD and AILD methods in Figs. 10 and 11 respectively. Due to the sensitivity of the low energy region of the phonon DOS in graphite to periodic boundary conditions present in AIMD, deviations of at most 30% occur in the inelastic scattering cross-section at energies below the thermal energy. However, above this energy range the predicted cross-section and secondary spectrum are in reasonable  agreement between the two methods. Furthermore, it may be observed that the low energy features of the phonon DOS, see Fig. 7, become more prominent with increasing temperature. This effect may partially explain the increasing agreement between AIMD and AILD in graphite with increasing temperature which is in contrast to beryllium.
Despite deviations in the inelastic cross-section, the secondary energy spectrum predicted at 300 K demonstrates reasonable agreement between the methods above the thermal energy. At the thermal energy the behavior predicted by both methods is consistent; however, deviations in the magnitude of the secondary spectrum predicted by AIMD and AILD occur due to the deviations in the low energy phonon DOS.

Conclusion
An ab initio molecular dynamics technique to produce phonon density of states for the generation of inelastic thermal neutron scattering libraries was evaluated using beryllium and graphite as exemplar solid moderators. Phonon DOS predicted by AIMD simulations were within reasonable agreement with the established method of ab initio lattice dynamics which suggests this method may be used to generate reliable phonon DOS [3,4]. The inelastic thermal neutron scattering cross sections generated from the phonon DOS for the AIMD and AILD methods were also in reasonable agreement for beryllium; however, some disagreement was found for graphite. Nevertheless, this disparity is attributed to size constraints of the current model, and may be further refined. Consequently, the AIMD method presents a first principles approach for the prediction of the dynamic structure factor (i.e. scattering law). This approach may be envisioned as a predictive approach for the generation of the thermal scattering law. Moreover, this approach may be of particular importance in the generation of thermal scattering cross-sections in materials (and/or under conditions) that exhibit significant anharmonic motion that prevents the use of AILD.