-
Notifications
You must be signed in to change notification settings - Fork 23
Waasmaier-Kirfel form factor coefficients lose precision due to missing Fortran kind specifiers #29
Description
Summary
The Waasmaier-Kirfel X-ray form factor coefficients in lib_f90/element_data_mod.f90 are declared as real(kind=PREC_DP) arrays, but the literal constants in the initialisers lack _PREC_DP (or D0) kind suffixes. This causes gfortran to parse them as default real (float32) before promoting to float64, silently discarding ~7 digits of precision. The practical impact on typical calculations is negligible (~1e-7 relative), but the fix is mechanical and affects every element and ion.
Minimal reproducer
program precision_test
implicit none
integer, parameter :: DP = selected_real_kind(15, 307)
real(kind=DP), parameter :: with_kind = 17.958397999999998973_DP
real(kind=DP), parameter :: without = 17.958397999999998973
write(*,'(a,ES25.18)') ' With kind suffix: ', with_kind
write(*,'(a,ES25.18)') ' Without kind suffix: ', without
write(*,'(a,ES25.18)') ' Difference: ', with_kind - without
end programOutput (gfortran 14):
With kind suffix: 1.795839799999999897E+01
Without kind suffix: 1.795839881896972656E+01
Difference: -8.189697275895468920E-07
Affected arrays
All per_wa1 through per_wa5, per_wb1 through per_wb5, and per_wc arrays (lines 325–610) use bare literals. The effect is ~1e-7 relative error per coefficient, propagating to ~1e-7 relative error in the computed form factors.
Suggested fix
Add _PREC_DP suffixes to all literal constants in the WK coefficient arrays. For example, change:
17.958397999999998973, 12.063053999999999277, ...to:
17.958397999999998973_PREC_DP, 12.063053999999999277_PREC_DP, ...