Jumps in T-dependence of inelastic cross section

[marquezj]

With @dddijulio we came across an issue when plotting XS as a function of temperature. There are fluctuations in the scattering cross section that are bigger than the variation with temperature. This could be important for quantitative wavelength-dependent transmission analysis or other applications where you need the cross section as a function of temperature.

This is an example of what we are seeing:

E_thermal = 0.0253
cfg = 'Polyethylene_CH2.ncmat'
temperatures = np.linspace(250, 350, 101)
plt.figure()
for lux in (3,4):
    xs_thermal = []
    for T in temperatures:
        mat = NC.load(cfg+f';temp={T};vdoslux={lux}')
        xs_thermal.append(mat.scatter.xsect(E_thermal))
    plt.plot(temperatures, xs_thermal, label=f'vdoslux={lux}')
plt.legend()
plt.xlabel('Temperature (K)')
plt.ylabel(f'Cross section @ E={E_thermal}')
plt.show()

A more intense effect is seen with the Debye Model, and perhaps it could be useful for debugging:

E = np.geomspace(1e-4,10,1000)
for T in  (319., 319.001):
    mat = NC.load(f'solid::C/1gcm3/TDebye400K;vdoslux=5;comp=inelas;temp={T}')
    xs = mat.xsect(E) 
    plt.semilogx(E, xs, '-')
plt.xlabel('Energy [eV]')
plt.ylabel('Scattering XS')

temperatures = np.linspace(300,320,101)
E0 = 1e-3
xs25 = [NC.load(f'solid::C/1gcm3/TDebye400K;vdoslux=5;elas=0;temp={t}').xsect(E0) for t in temperatures]
plt.plot(temperatures, xs25, '.-')
plt.xlabel('Temperature [K]')
plt.ylabel('Scattering XS at 1 meV')

[tkittel]

So my gut feeling is that this effect is 100% related to the numerical artifacts induced by the way we force everything on first a VDOS grid and then later a SAB grid. VDOSDebye is getting an incredibly crude grid only, so the fix there is most likely to make it less crude. I was being rather brutal, since I thought VDOSDebye is anyway a super crude model. But if smoothness of the particular curves you show is important, I guess that could be revisited.

For the other example you show I don't actually think the effect is very big, considering the gridding artifacts, and most likely it is much smaller than the actual measurement errors.

A solution would most likely involve changing the procedure in a fundamental way. If you are mainly worried about effects at the large wavelength scale (where indeed the narrow phase space and the coarse grids are really at odds), then the solution might be to move forward with the planned direct sampling of the G functions (at least for G1, possibly G2). Then we could leave G1+G2 out of the SAB entirely which would certainly be nice.

Interesting how such direct sampling would intersect with the liquids of course, and if it would still be meaningful.

[tkittel]

I am removing the bug label, since I consider this a "known deficiency" of the approach we are using. But I am adding the enhancement label since of course there is room for improvement!

[dddijulio]

Just one observation I had when looking at the poly example, if I set, egrid 1e-5 5.0 300, or egrid 1e-5 5.0 900, the curves are smoother, but there is a systematic difference between the two.

[tkittel]

I don't have many cycles available right now, but things that would be helpful for the investigation:

• Check if improving the egrid would help (but be aware of the linear-vs-log interpolation-vs-plot issue that might make it hard to draw conclusions if you only twiddle it a bit - add a huge amount of more points.
• Check if it happens also in non-H materials. Preferably a mono-atomic material (you can always mock one up if needed). Also make sure to run comparisons with comp=inelas to avoid noise from DebyeWaller factors.
• Check if just the G1 contribution would do it (needs a @CUSTOM_... hack in the .ncmat, I'll come back with a recipe.
• For VDOSDebye specifically one can increase the number of VDOS points from 20 to e.g. 400 in line in the NC::createVDOSDebye function at the bottom of the file ncrystal_core/src/dyninfoutils /NCDynInfoUtils.cc (remember you can just pip install ./ncrystal_core ./ncrystal_python)
• Check if abusing the ucnmode parameter to be more careful near the kinematic endpoint of the phase space helps: ucnmode=refine:0.01eV.

[tkittel]

For how to get just the G1 function, you can for instance use the NCMATComposers hidden function: mycomposer._unofficial_vdos2sab_ignore( 2 ) will ignore G2 and above.

Or just add in the NCMAT data:

@CUSTOM_UNOFFICIALHACKS
vdos2sab_ignorecontrib 2

Full syntax is:

@CUSTOM_UNOFFICIALHACKS
vdos2sab_ignorecontrib low [high] [coherent|incoherent]

(this is a feature added to support the BzScope plugin by @xxcai1 et. al.)

[marquezj]

I will leave this here for further analysis: a couple of cases for polyethylene where it is particularly noticeable.

plt.figure()
temps = (306.265, 306.270)
E = np.geomspace(1e-4,10,1000)
cfg = 'Polyethylene_CH2.ncmat'
xs = {}
for t in temps:
    mat = NC.load(cfg+f';vdoslux=3;ucnmode=refine:0.01eV;comp=inelas;temp={t}')
    xs[t] = mat.xsect(E)
plt.semilogx(E, xs[temps[1]]/xs[temps[0]]-1)
plt.xlabel('Energy [eV]')
plt.ylabel('XS (rel. diff)')
plt.show()

plt.figure()
temps = (294.48, 294.49)
E = np.geomspace(1e-4,10,1000)
cfg = 'Polyethylene_CH2.ncmat'
xs = {}
for t in temps:
    mat = NC.load(cfg+f';vdoslux=4;ucnmode=refine:0.01eV;comp=inelas;temp={t}')
    xs[t] = mat.xsect(E)
plt.semilogx(E, xs[temps[1]]/xs[temps[0]]-1)
plt.xlabel('Energy [eV]')
plt.ylabel('XS (rel. diff)')
plt.show()

[marquezj]

The discontinuity is actually in energy, but the temperature variation makes it more obvious:

temperatures = [306.26, 306.28]
E = np.linspace(7e-3,20e-3,1000)
xs = {}
for t in temperatures:
    mat = NC.load(cfg+f';vdoslux=3;comp=inelas;temp={t}')
    xs[t] = mat.xsect(E)
for t, val in xs.items():
    plt.plot(E, val, label=f'{t} K')
plt.legend()
plt.xlabel('Energy [eV]')
plt.ylabel('Scattering XS')

[tkittel]

So a discontinuity in energy is rather interesting!

[tkittel]

For what it is worth, the jumps are truly discontinous:

Here I was forcing the egrid to use 3e4 rather than 3e2 points. So the fact that they appear non-discontinuous is just an interpolation artifact.

[tkittel]

I think I have a pretty good idea where this is coming from. Basically from my suspicion all along. But it might be a bit tricky to fix (although hopefully not impossible).