ICE-MAN the Integrated Computational Environment for Modeling and Analysis for Neutrons at ORNL

. ICE-MAN is a modeling and analysis workbench for multi-modal studies, designed with neutron science in mind. It streamlines the workflow between di ff erent experimental techniques, computer modeling, and databases and reduces the time and learning curve needed to access them thus making a holistic approach to data interpretation more amenable and e ffi cient.


Introduction
One of the most significant strengths of neutron scattering is the direct correlation between theoretical models and experiment; it is possible to calculate the neutron scattering cross-section rigorously from models [1].However, most tools used to analyze data from atomistic modeling are qualitative rather than quantitative.Many standard fitting tools do not necessarily include atomistic detail, particularly those extracted from electronic structure calculations.This paper describes a platform for analyzing neutron scattering data that can be broadly applied across instrument suites at SNS and HFIR, building a suite of individual data analysis/processing engines that streamline workflows in data analysis for neutron scattering, the Integrated Computational Environment for Modeling and Analysis for Neutrons (ICEMAN).

The OCLIMAX module
The OCLIMAX software calculates the neutron scattering synthetic spectra, S(Q, ω), for solids.When doing lattice dynamics calculations, and the harmonic approximation the results are rigurous, and it is necessary to provide the eigenvectors (displacement vectors) and eigenvalues (frequencies) for the calculation.OCLIMAX continues the naming tradition first established by Kearley and Tomkinson [38] and later by A.J. Ramirez-Cuesta [1,39].The initial version used force constants to fit INS spectra of simple molecules, the second was able to produce INS spectrum from the output of DFT modeling, for details see [2].The conventional approach using density functional theory and lattice dynamics often falls short when the materials of interest are complex (e.g., defective, disordered, heterogeneous, amorphous, large-scale), for which molecular dynamics driven by interatomic force fields are a more common approach.OCLIMAX can directly convert molecular dynamics trajectories into simulated INS spectra, including not only fundamental but also higher order excitations [3].

Main features
OCLIMAX can calculate the spectra for extended solids, molecules in solids (using the isolated molecule approach).It can introduce solid state effects in the isolated molecule approach by using the shape of the phonon wings to approximately calculate overtones and combinations.For the calculation of the INS spectra of molecules  it uses the incoherent approximation, that essentially consists of using the total cross section of each atom.
For extended solids with periodic boundary conditions it is possible to calculate the incoherent and coherent contributions to S(Q, ω).The contribution from elastic scattering can also added.The sum of all contributions is the total scattering function.The program can calculate single crystal and powder averaged spectra of S(Q, ω), the whole map and along a given trajectory, like VISION.It can calculate direct and indirect geometry instrumentation.Instrument resolution can be entered as an empirical polynomial, or as in the case of direct geometry spectrometers, using instrumental parameters.It is also possible to calculate generic S(Q, ω) maps covering given areas of Q and ω with arbitrary resolution.

Examples
A series of examples are shown in figure 1.On the left column is the experimental and calculated S(Q, ω) for graphite as seen is SEQUOIA.At the bottom of the left column is the spectra of graphite as measured in VISION, with the calculated spectra using the incoherent approximation (blue trace) and the full coherent calculation (red trace).It can be easily seen that the coherent effects have to be taken into consideration to fully reproduce the experimental data.
In the mid column are the S(Q, ω) from DFT calculation top and middle, at the bottom is the experimental data of a single crystal of α − T eO 2 measured at HySPEC.Upon doing the same integration and range as measured data, i.e. including the correct momentum resolution, the calculation effectively reproduces the experimental data.
On the right column, the top two panels show the neutron energy gain for graphite and diamond at 300 K.The bottom panel shows the calculated total neutron cross section for both materials.The origin of the different transmission at low incident energy is clearly illustrated.Graphite has a low-lying phonon branch below 15 meV, which is the LA/LO mode and is highly populated at room temperature.The ω, Q trajectory at low incident energy cuts through the branch in between 5 and 15 meV, contributing to its inelastic scattering cross-section.In con-trast, diamond does not have a major phonon branching lying in the same ω, Q trajectory.It only crosses a weak band between 15 and 20 meV at higher Q, resulting in the lower cross-section.It is worth noting that such insight can only be obtained through a full calculation of this kind including the coherent effects; the conventional approach using the phonon density of states will be far less revealing [4].

Workflows with OCLIMAX
The OCLIMAX program and its manual (including installation and user guide) can be found on the Web [5].OCLIMAX requires an .aclimax(incoherent approximation simulation only) or .oclimax(required by coherent simulation) file, as well as an optional .paramsfile to run.[2] To convert MD trajectories to INS spectra, OCLIMAX can be used in essentially the same way as reported previously, with the same format for the parameter file.Here instead of having the .aclimaxor .oclimaxfile for phonon information, we need a .tclimaxfile for the trajectory.[3] The .aclimax/.oclimaxfile contains information on the atomic structure, isotope, neutron scattering cross sections, and phonons, and the .paramsfile contains parameters defining the simulation (if different from the default).The .aclimax/.oclimaxfile can be obtained from the output file(s) of a phonon calculation (e.g., using Phonopy), by running the conversion script provided in the program.Currently the script supports output files from CASTEP, VASP (either directly from VASP or through Phonopy), CP2K, Quantum Espresso, Phonopy, ORCA, DMol3, NWChem, etc.
In the .paramsfile, users can specify how the calculation should be done.The default setup is for VI-SION/TOSCA spectrum simulation.One can modify the .paramsfile to include combinations and overtones (up to a desired order), temperature effect, coherent effect, isotope substitution, and phonon wing calculation.The users can select to simulate the full S(Q, ω) map or certain sampling trajectories corresponding to a (real or virtual) direct or indirect geometry instrument.A detailed explanation of the input parameters can be found in the manual.The program is currently distributed as a Docker image [40], allowing it to run on various platforms with minimum environment dependence.For most calculations with optimized phonon sampling density and parameters, it should take no more than a few minutes to finish.Also distributed on the Web site is a script to generate quick plot for visualization of the calculated spectra (pclimax.py).
The execution of OCLIMAX is a the command line call: $ OCLIMAX run yourfile.oclimax.A conversion tool is included in the OCLIMAX Docker distributions.For help with these, run $ OCLIMAX convert on the command line, the output will describe the available options.
The standard output data is in .csvformat.The users can also select to generate a restart file (in binary format), an SPE file (for DAVE [41]), or (neutron weighted) phonon density of states.Overall it is much more powerful than its previous version aCLIMAX [39], as well as the recently developed ABINS algorithm [42] in Mantid [43].

The Qclimax module
Qclimax provides users with a flexible, straightforward environment to enable fitting to complex, user-defined functions.The power of QENS arises from the ability to access both momentum and energy transfer domains at the same time, providing information about the characteristic diffusion processes.It simultaneously allows fits to be performed to data from multiple Q values and provides global fitting parameters.
Typically, the measured intensity of the QENS experiments, I(Q, ω), can be described as follows [44]: where, Q is the momentum transfer and ω is the energy transfer, x(Q) is the elastic fraction of the signal; δ(ω) is a Dirac delta function accounting for the referred elastic scattering (zero energy transfer);S (Q, ω) is the dynamic structure factor, i.e., the QENS model scattering function; B(Q, ω) is a linear background; and R(Q, ω) corresponds to the instrument resolution function, the symbol ⊗ represents the convolution operation with the resolution function .In general, B(Q, ω) accounts for dynamical processes that are too fast for the dynamic range of the instrument (that is, with timescales shorter than a few picoseconds) and, possibly, sample-and temperature-dependent backgrounds.As an example, S (Q, ω) can be described by a Lorentzian function: where Γ (Q) is the half width at the half maximum of the curve.
In the simple case of a single Lorentzian process with a linear background, the equation to fit the experimental data is: The sequential fitting procedure consists in binning the data in a number of values in Q and fitting Equation 3 to the data for each individual value of Q.The parameters of the fits as a function of Q, are, in this case: x(Q), Γ (Q), A(Q), and B(Q).The parameters are then fitted to a given functional form in Q.In this case, for example, the generic equation for Γ (Q) is: where α, β, . . .are the parameters of interest that contain the physical meaning of the phenomena.QClimax allows to include the functional form in Q within the fitting process.The results of the fitting are the parameters that fit Equation 3 given the constraints from Figure 2 shows how the use of global fitting.On the top panel, the sequential fit is not good on the last few points at large Q.The global fitting, on the other hand, introduces in the analysis a "prior knowledge" of the system, i.e., it behaves as a Chudley-Elliott model.This "prior knowledge" is an intrinsic Bayesian behavior that makes the fittings more robust.The bottom panel, T = 750K finds similar results using both methods, this is due to the better quality of the sequential fitting.
Qclimax has been written in python and uses the LM-FIT python library [45]; it can be accessed through the ICEMAN homepage [46] (it is necessary to register for XCAMS access [47]

to use it).
There are some predefined peak shapes built in QClimax.To add extra peaks shapes, please contact the authors.QClimax allows the user to create new constraint functions in Python.These functions once created are saved in the user's profile.All the data resulting from the fitting is saved in the profile.It also saves the initial data files, parameters, fitting functions, and constraint functions.The results can be downloaded as a .zipfile and shared.All the information contained can be used to fully reproduce the fitting results.
The experimental files are in DAVE .grpformat [41].It is recommended to use µeV as units in energy.The initialization file .inicontains the required information for the fittings.

Conclusions
The ICEMAN project is still in development.The ICE-MAN software with its two modules, OCLIMAX and QClimax, is available for use.OCLIMAX is available as a Docker container, QClimax can be accessed through the ICEMAN webpage.For QClimax there are Virtual Machine builds and soon it will be a Docker version available.For more information, suggestions, collaborations, etc., please contact any of the authors of this paper.
This research benefited from the use of the BASIS, HySPEC, SEQUOIA ,and VISION beam-lines at the Spallation Neutron Source a DOE Office of Science User Facility operated by the Oak Ridge National Laboratory.The computing and software resources were made available through the VirtuES and the ICEMAN projects, funded by the Laboratory Directed Researchand Development program (LDRDs 7739, 8237,10447) and Compute and Data Environment for Science (CADES) at the Oak Ridge National Laboratory, which is supported by the Office of Science of the DOE under Contract DE-AC05-00OR22725

Figure 1 .
Figure 1.Examples of OCLIMAX applications.left column top and center, the S(Q, ω) map for graphite powder, bottom the INS spectra calculated for graphite along the VISION trajectory.Middle column, phonon dispersion along the [1, 0, L] transverse cut using the HYSPEC spectrometer [13].DFT simulations (a) without integration and (b) with same integration and range as measured data in (c).(c) Measured INS data integrated at −0.12 < ζ < 0.12 in (1 + ζ, 0, 0) and −0.12 < η < 0.12 in (1, η, 0).Right column, OCLIMAX calculation of the neutron energy gain part of the S(Q, ω) map for (top) graphite and (center) diamond.The red line indicates the trajectory where nearby phonons have major contributions to the scattering cross-section for low-energy (< 1meV) neutrons.Note that to see the weak phonon excitations in diamond, the scale of the color bar for the right panel is only 1/10 of the left panel.

Figure 2 .
Figure 2. Example of fitting data using sequential and global fits in the case of BaH 2 [32], for T = 670K at the top and T = 750K at the bottom panel.The black points are the fitted values of Γ (Q) using sequential fitting (both DAVE and QClimax give the same values), the black trace is the fitting of the data to the Chudley-Elliott model, including all values.The red points are the fitted points using global fitting in QClimax, the red trace is the overall model.