Automated workflow for accurate scattering function calculation of light water using quantum corrections

1. Introduction

Neutron scattering is widely used across various fields, such as materials science, chemistry, energy, and condensed matter physics. Compared to synchrotron radiation, neutrons have a unique advantage in detecting light elements such as hydrogen. However, the strong incoherent scattering signal from hydrogen dominates the total signal, obscuring the coherent scattering that carries structural information and complicating data analysis. As a result, materials containing hydrogen often require deuterium substitution in neutron experiments, which limits the range of samples that can be studied. Light water represents a particularly important and challenging case due to its ubiquity and the complexity of its hydrogen dynamics. To overcome these limitations and expand the application of neutron scattering, accurate theoretical calculations of the scattering function for light water materials are essential.

Theoretical approaches provide a foundation for developing robust data analysis and correction methods, enabling more precise interpretation of experimental results. Methods for calculating the neutron scattering function of solids have been developed for decades, primarily based on the harmonic model.^1–3 Currently, NCrystal library ⁴ is widely used to support a range of physical calculations, including coherent, incoherent, elastic, and inelastic scattering in solid materials. However, the development of methods for calculating neutron scattering functions in liquids has lagged behind. Unlike solids, there is no universal approach for liquids, and different models are adopted depending on specific applications.^5–7 Due to the limited accuracy and uncertain applicability of these models, their selection often relies on empirical judgement.

Van Hove’s theory ⁸ provides the fundamental framework for calculating the neutron scattering function through space–time correlation functions. For the computed scattering function to be physically meaningful, it must satisfy three fundamental conditions ⁹ : (a) the sum rule, (b) detailed balance, and (c) numerical stability over a large dynamic range to avoid Fourier-transform artifacts. Advances in molecular dynamics (MD) simulations have made it possible to perform calculation using ab initio methods.^10–14 In 2014, the Granada team proposed the so-called Centro Atómico Bariloche (CAB) model, ¹⁵ named after their institution, which utilizes the density of states as an input and employs a Gaussian approximation (GA) to describe the incoherent scattering function. This model addresses a series of numerical issues in Fourier transforms through three physical approximations and was adopted by the U.S. Nuclear Data Library in 2018.

However, the CAB model requires a significant number of experimental parameters, limiting its applicability to calculations for liquid materials other than light water. Additionally, the CAB model relies on classical MD simulations results as input, which are only valid at scales much larger than the quantum regime.

The need for quantum corrections in neutron scattering calculations has long been recognized. While several approaches have been proposed over the years,^16–18 this paper focuses on the recent Gaussian approximation-assisted quantum correction (GAAQC) method introduced by Abe and Tasaki, ¹⁹ which approximates the width function Γ ( t ) as Gaussian in time t , offers a clear and efficient framework. However, Abe and Tasaki’s method exhibits significant limitations. First, its effective analysis window is restricted to only 10 ps, which is insufficient to capture the full temporal evolution of the diffusion process typically requiring windows of at least 50 ps. More critically, as pointed out by the Granada group, ²⁰ the direct calculation using two Gaussian models within Abe and Tasaki’s framework becomes numerically unstable in the high momentum transfer (large- Q ) region. As a result, in Abe and Tasaki’s original work, only partial results within a limited ( Q , ω ) range were presented, and the full scattering function S ( Q , ω ) could not be reliably reconstructed.

In this work, we build upon a recently developed high-precision Fourier transform technique to overcome these numerical limitations for light water. We achieve a complete and robust open-source tool that applies the GAAQC quantum correction to classical neutron scattering data. The tool takes two inputs—both obtainable from standard MD trajectories: (a) the vibrational density of states (VDOS) and (b) the classical scattering function S cl ( Q , ω ) . And it outputs the quantum-corrected S ( Q , ω ) over a wide range of ( Q , ω ) . This work thus transform Abe and Tasaki’s theoretical formalism into a practical computational resource, providing a reliable foundation for future studies of hydrogenous materials.

2. Theory of quantum correction

The scattering function can be divided into coherent and incoherent components. In hydrogen-rich liquid materials such as light water, the contribution of coherent scattering between different atoms is less than 5% of the self-scattering. ⁹ Therefore, the incoherent approximation is adopted, which neglects cross-correlations between distinct atoms. The scattering function is defined as

S ( Q , ω ) = 1 2 π ∫ − ∞ ∞ d t F ( Q , t ) e − i ω t

(1)

A classical approach to computing the scattering function involves utilizing data derived from MD simulations. Specifically, this function can be evaluated based on atomic trajectory information obtained from such simulations:

G c ( r , t ) = ⟨ δ ( r − R i ( t + τ ) + R i ( τ ) ) ⟩ c

(2)

C ( Q , ω ) = 1 2 π ∫ e − i ω t d t ∫ e − i k r G c ( r , t ) d r

(3)

C ( Q , ω ) = lim ℏ → 0 S ( Q , ω )

(4)

S ( Q , ω ) = exp ( ℏ ω k B T ) S ( Q , − ω )

(5)

For non-quantum liquids, the relationship between the quantum-mechanical and classical neutron scattering functions can be expressed as: ¹⁸

S ( Q , ω ) = ∫ − ∞ ∞ R ( Q , ω ′ ) C ( Q , ω − ω ′ ) d ω ′

(6)

R ( Q , ω ) = 1 2 π ∫ F q ( Q , t ) F c ( Q , t ) e − i ω t d t

(7)

Using the Gaussian approximation (GA) framework, a complete computational pathway from molecular dynamics trajectories to the dynamic structure factor S ( Q , ω ) can be achieved. The physical justification for the GA is that the self-intermediate scattering function is exactly Gaussian in two important limiting cases: for recoil scattering from free atoms at large scattering vectors, and for quasi-elastic neutron scattering (QENS) from diffusing atoms at small scattering vectors. ¹⁹ The F ( Q , t ) can be defined as:

F ( Q , t ) = exp [ − Q 2 2 Γ ( t ) ]

(8)

Γ q ( t ) = ℏ m ∫ 0 ∞ d ω ρ ( ω ) ω [ coth ( ℏ ω 2 k B T ) ( 1 − cos ω t ) − i sin ω t ]

(9)

Γ c ( t ) = ℏ 2 k B T m ∫ 0 ∞ d ω ρ ( ω ) ω 2 ( 1 − cos ω t )

(10)

ρ ( ω ) = 2 m 3 π k B T ∫ 0 ∞ d t ⟨ v ( 0 ) ⋅ v ( t ) ⟩ c cos ( ω t )

(11)

According these equations above, the entire computational workflow of the method is shown in Figure 1. The following section details the key steps within this workflow, along with the specific challenges encountered and the solutions developed to address them.

OPEN IN VIEWER

Figure 1.

The workflow of quantum correction for light water.

The GAAQC tool requires two inputs obtainable from classical MD: the VDOS and the classical scattering function S cl ( Q , ω ) . To generate the inputs used for validation in this work, we performed MD simulations of liquid water using LAMMPS ²² with the TIP4P/2005 model. A system of 512 water molecules was simulated at 293 K in a cubic periodic box. After NPT equilibration, a 100 ps NVT production run was conducted with a 0.2 fs timestep. The VDOS was computed from the VAF, and S cl ( Q , ω ) was obtained via Fourier transform of the intermediate scattering function.

The computational workflow encounters two major challenges in evaluating the width function Γ ( t ) . To address the issues of high oscillations and the removable singularity at ω = 0 , a hybrid integration strategy is employed: (a) for the singular region near ω = 0 , the contribution is computed via a “Taylor limit” treatment, where the integrand at ω 0 = 0 is evaluated using its analytic small-argument limit, and the integral over [ 0 , ω 1 ] is approximated by the trapezoidal rule; (b) the remaining oscillatory integral ( ω ≥ ω 1 ) is evaluated using a custom implementation of Filon’s quadrature, ²³ following the formulas detailed in Appendix A of Du and Cai. ²⁴ The complete C++ source code for both the low-frequency handler (cal_limit) and the Filon routine (cal_integral) is provided in the accompanying repository under the GPL License. The calculated width function of light water for both quantum and classical regimes is shown in Figure 2.

OPEN IN VIEWER

Figure 2.

The width function of light water in its classic and quantum state, calculated using equations (9) and (10). The insert highlights the differences in the real part, imaginary part, and amplitude over a short-time interval.

The second major challenge arises in the application of the Fourier transform to obtain the scattering function S ( Q , ω ) . Note that in this work, the energy transfer Δ E is defined as E f − E i . Therefore, positive Δ E on the horizontal axis corresponds to neutron energy gain. At high momentum transfers, a Fast Fourier Transform (FFT) of F ( Q , t ) becomes numerically unstable. As Q increases, the exponential factor causes F ( Q , t ) to decay below the noise level for all but the shortest times (see in Figure 3), leading to severe artifacts in the resulting S ( Q , ω ) . These artifacts manifest as spectral aliasing and unphysical oscillations, rendering the standard FFT approach unreliable for large Q .

OPEN IN VIEWER

Figure 3.

The intermediate scattering function F ( Q , t ) for various momentum transfers Q , calculated using equation (8). At high Q values (e.g., Q = 40  Å⁻¹), F ( Q , t ) decays within a very short time.

To overcome this limitation, we exploit a key property of the Gaussian approximation: the convolution theorem allows us to express the Fourier transform of exp [ − Q 2 Γ ( t ) / 2 ] as a convolution of results computed at smaller Q values. Specifically,

S ( Q , ω ) = 1 2 π F { exp [ − Q 2 2 Γ ( t ) ] } = 1 2 π F { exp [ − ( Q / 2 ) 2 2 Γ ( t ) ] 2 } = S ( Q 2 , ω ) ⊗ S ( Q 2 , ω )

(12)

Following the methodological framework established in the preceding sections, this section presents a comprehensive analysis of the scattering function of light water calculated using our refined CDFT approach. The discussion is structured into two main parts to systematically validate the methodology. First, we evaluate the scattering function S ( Q , ω ) obtained within the GA framework. Subsequently, we apply and validate the GAAQC model, confirming its consistency with established theoretical limits and experimental data.

4.1. Evaluation of the scattering function calculated via refined CDFT (GA)

The scattering function S ( Q , ω ) obtained from the refined CDFT methodology via GA model is shown in Figure 4(a) and (b), which display the scattering function S ( Q , ω ) for small and large Q , respectively. The STC approximation, based on instantaneous collisions, requires Q ≫ 2 π / d where d ≈ 3 Å is the characteristic collision length in water. ²⁵ For Q = 10 Å − 1 , this condition is not fully satisfied, explaining the discrepancy with the CDFT result in Figure 4(b). The convergence of the two methods at higher Q ( ≥ 40 Å − 1 ) confirms that CDFT correctly recovers the STC limit when the short-time-collision assumption becomes valid. The consistency of these results with established physical expectations and our previous work ²⁴ validates the robustness of the CDFT methodology.

OPEN IN VIEWER

Figure 4.

(a) Scattering function calculated from refined CDFT via GA model at Q = 1 Å − 1 , demonstrating physically reasonable line shape and detailed-balance compliance. (b) Comparison of CDFT (solid lines) and STC (dashed lines) results for Q = 10 Å − 1 to 80 Å − 1 . At Q = 10 Å − 1 , the discrepancy reflects non-instantaneous collision effects captured by CDFT but not by the STC approximation. The convergence at higher Q ( ≥ 40 Å − 1 ) validates that CDFT correctly approaches the STC limit when short-time collisions dominate.

4.2. Validation of quantum correction of the scattering function

The quantum-corrected scattering function S ( Q , ω ) is shown in Figure 5, which refers primarily to the agreement of the quasi elastic peak positions with the STC predicted trajectory (dashed “STC max” line). The STC model provides a well established relationship between the peak location in and momentum transfer for systems dominated by short time collisions. The alignment of our computed peak maxima with this theoretical curve across a wide range Q validates that quantum-corrected scattering function correctly captures the short time dynamics underlying the STC limit.

OPEN IN VIEWER

Figure 5.

The GAAQC scattering function is in agreement with the STC model.

Moreover, as shown in Figure 6, the half width at half maximum (HWHM) of quasi-elastic peaks at small momentum transfers is compared. The HWHM derived from GA exhibits a linear relationship with Q 2 , consistent with Fick’s law of continuous diffusion. In contrast, GAAQC HWHM is accurately described by the jump-diffusion model, ²⁶ yielding a diffusion coefficient D = 0.23 Å 2 ps − 1 and a residence time τ = 0.47  ps upon fitting. This agreement underscores the necessity of quantum corrections to capture the true microscopic motion of water. Crucially, the line-width correction achieved by GAAQC cannot be reproduced by a simple multiplicative detailed-balance factor. Such a factor only enforces the intensity asymmetry S ( Q , − ω ) = e − β ℏ ω S ( Q , ω ) without altering the Q -dependence of the line shape.

OPEN IN VIEWER

Figure 6.

Comparison of HWHM derived from GA and GAAQC model.

The double differential scattering cross section calculated for an incident energy of E 0 = 0.154  eV agrees with experimental result from Esch et al. ²⁷ at scattering angles of 25 ∘ and 40 ∘ (Figure 7). Discrepancies, particularly in the quasi-elastic peak widths and the inelastic wings, are observed. These are likely influenced by the fact that the instrumental resolution function is not convolved with the theoretical result. The overall agreement validates the accuracy of our quantum-corrected computational framework for capturing the main scattering features, while a more precise quantitative comparison would require convolution with the experimental resolution.

OPEN IN VIEWER

Figure 7.

Comparison of the double differential scattering cross section between simulation and experiment ²⁷ at E 0 = 0.154  eV. Upper panel: scattering angle 25 ∘ . Lower panel: scattering angle 40 ∘ . The measured data (symbols) is compared with GAAQC results (solid lines).

This paper presents a computational workflow that applies the GAAQC framework to classical neutron scattering function S ( Q , ω ) of light water, overcoming long-standing limitations in quantum correction methods. By integrating high-precision Fourier transform techniques with advanced quantum corrections, we have developed a robust approach that maintains numerical stability across an extensive dynamic range in both momentum and energy transfer.

The key achievement of this work is the transformation of theoretical GAAQC framework into a practical, open-source implementation. The tool takes as input the VDOS and the classical scattering function S cl ( Q , ω ) —both obtainable from standard MD simulations—and outputs the quantum-corrected S ( Q , ω ) . Unlike empirical models, our implementation introduces no additional fitted parameters in the correction step and ensures physical consistency through automatic satisfaction of fundamental constraints including the sum rule and detailed balance condition. The complete source code and documentation are publicly available at https://github.com/Jecrevin/sqw, providing a robust computational foundation for future studies of hydrogenous liquids, beginning with light water.