Support inelastic vdos expansions also at very low temperatures

[dddijulio]

We see a strange behavior in the cross-section as it goes to the free atom region when calculating at low temperatures. For example,

ncrystal_inspectfile MgH2_sg136_MagnesiumHydride.ncmat -c temp=1.5K, appears OK, while ncrystal_inspectfile MgH2_sg136_MagnesiumHydride.ncmat -c temp=1.0K, produces the behavior we are seeing.

It happens also in Al for instance, but the temperature needs to be lower. ncrystal_inspectfile Al_sg225.ncmat -c temp=0.005K

Douglas

[tkittel]

I confirm that there is something fishy going on, possibly related to the SAB grid. For instance one can play with the following command and change temp from, say, 1 to 0.005 and see what happens and how it depends on the grid (temp=0.001 triggers an assert in debug builds - which should have been a CalcError exception instead):

ncrystal_inspectfile 'Al_sg225.ncmat;vdoslux=5'  'Al_sg225.ncmat;vdoslux=4' 'Al_sg225.ncmat;vdoslux=3' 'Al_sg225.ncmat;vdoslux=2' 'Al_sg225.ncmat;vdoslux=1' 'Al_sg225.ncmat;vdoslux=0' -c comp=inelas -c temp=0.005

There are likely some assumptions, since I didn't imagine crazy people going to around 1K when I was developing a lot of the code. But certainly it is worthwhile to try to fix this, or possibly at least find a workaround.

[tkittel]

The warnings from the above inspectfile command are likely a hint:

NCrystal WARNING: Discarding 1326 edges of provided kernel data due to missing S values.
NCrystal WARNING: Discarding 618 edges of provided kernel data due to missing S values.
NCrystal WARNING: Discarding 275 edges of provided kernel data due to missing S values.
NCrystal WARNING: Discarding 117 edges of provided kernel data due to missing S values.
NCrystal WARNING: Discarding 39 edges of provided kernel data due to missing S values.
NCrystal WARNING: Discarding 19 edges of provided kernel data due to missing S values.

They shouldn't really occur when using vdos expanded sabs, so it tells me that we might be looking at a numerical issue (e.g. exp(-largenum) becomes strictly 0.

I also quickly verified that disabling this trimming does not affect the final cross section curves, so at least the trimming itself is safe and does not introduce artifacts.

[XuShuqi7]

Hi Thomas,

I tried with the additional command "vdoslux", but it seems that the problems of the inelastic behaviors at very low temperatures still exist. Here we can see two examples: Beryllium and Magnetism hydride: We can see clearly the cutoff energies of the vdoslux options (0.5, 1, 3, 5 eV,...). Beyond these cutoff energies, the neutron scattering cross sections converge to the free cross sections since the Free Gas Model (FGM) is applied. In these two examples, interesting things are that the neutron cross sections obtained by integrating the S(alpha, beta) tables are not coincide with the FGM, if we continue to increase the grid cutoff. For possible reasons, I guess so, that the problem is perhaps due to some numerical issues, since we have 1/k_BT in many exponential functions.

Best regards, Shuqi

[tkittel]

Indeed, that was my point. However there might be more too it, since there is obviously a vdoslux dependency...

[tkittel]

Possibly related to this, the sampled beta distribution in Al at T=5mK and Eneutron=10kT shows clear issues with inadequate grid granularity (and this distribution does not appear to depend much on vdoslux):

[tkittel]

Looking at the generated SAB for Al@1mK I see that the "free gas like" slope towards higher energies is not correct (the data does not like centered on the dashed white line):

Increasing the temperature slowly also slowly aligns the two curves, and at 1K they are certainly on top of each other:

[tkittel]

I suspect we are witnessing two different problems here. One problem (seen in the cross-section curves) is related to the wrong slopes in the SAB plane, which is caused by too large numbers in the argument of an exponential function (I suspect that the reason that it works at all for beta values in the 1000s is because of some tricks I already implemented around https://github.com/mctools/ncrystal/blob/master/ncrystal_core/src/NCVDOSToScatKnl.cc#L146).

The other problem affecting the precision in the beta sampling at Eneutron ~= kT is likely due to too few points added around here: https://github.com/mctools/ncrystal/blob/master/ncrystal_core/src/NCVDOSToScatKnl.cc#L665 . Btw., this issue would also go away if I ever found time to look at my most-favourite-thing-that-I-havent-time-for, a.k.a. #40.

[tkittel]

So I just looked at the formula, and if I did them right, the β value where the dashed while line (a.k.a. β=-α/A) intersects the kinematic parabola (a.k.a. α+(β)=2e+β+2sqrt(e(e+β)) with e=E/kT) is β=-4eA/(A+1)^2. So if we want to keep beta higher than (for instance) -1000, we have to ensure that E<1000((A+1)^2/4A)kT. I will try to proceed and see if I can adjust the algs to take that into account.

[tkittel]

So here is an issue. Currently, we never allow the beta-grid to be contracted so much that the first order phonon contribution will not be fully encompassed. So for Al the VDOS goes to 0.04eV, meaning that at T=5mK the beta grid can not be contracted to less than approximately 0.04eV/kT=93000. We can of course override this, but doing so would mean that we force a transition between a calculated VDOS-based SAB and into the free-gas regime, after only a few phonons. This does not seem very proper unfortunately.

[tkittel]

So another random thought of mine was that perhaps this issue is related to how we still don't actually use the effective temperature for anything. If I simply trust the current code to calculate it correctly, it should be approx 149.5K, i.e. Teff/T~=30000. However, the observed slope in the SAB place is only wrong by a factor of approximately 50, so I am not really sure how this all relates...

[XuShuqi7]

Hi Thomas, is it possible that the contribution from phonon is not very significant at very low temperature? So maybe it is safe to contract the beta grid to a less important value? Shuqi

[tkittel]

@XuShuqi7 hmm... not really, our model is that the entire inelastic cross section comes from phonons. The transition to a free-gas model is just a more efficient (and numerically stable) way to represent maaaaaany phonons at once.

[tkittel]

So one thing that is definitely messed up here is the G1 function at 5mK:

Compared with how it looks at 1K:

Perhaps I simply need to make the G1 function more numerically stable and it will all be fine! (too optimistic perhaps...)

[XuShuqi7]

@tkittel thank you for your comments. I understand that the scattering kernels are computed from phonon expansion. My point is that the free gas model assumptions may be fulfilled at temperatures near to zero for the inelastic part. The G1 plots seem confirm my guess. The G1 function at 5mK looks like a Dirac function, which is the same as the phonon density of states for the free gas model.

[tkittel]

I don't think it is a dirac function exactly, but I am trying to fix it to be more precise and then we will se. Perhaps you are right, it is not really clear to me.

However, the problem with your implied solution is that then we need to come up with an algorithm for deciding when to transfer and not to a pure free gas algorithm, which implies a lot of room for getting such an algorithm wrong in general...

[XuShuqi7]

Yes, I agree on that how to decide if we meet the free gas model criteria is cumbersome compare to the direct use of phonon expansion algorithm.

[tkittel]

I think/hope the issue will be solved with a suitable rebinning of G1. However, rather than just doing that and looking at the final curve, I am currently trying to verify that such a rebinning actually is enough to give a correct G1. With the help of high-res floating points (python's mpmath module).

[tkittel]

FYI I am working all day on this today. I believe I have a fix for the G1 issue, just hunting down a sign error or something like that somewhere atm.

Btw., in debug builds it is hard to hunt down this since NCrystal emits:

NCrystal.NCLogicError: Assertion failure: m_temperature.get()>=1.0&&m_temperature.get()<1e5

Which kinds of moves this from "bug" to "feature request" territory :grin:

[tkittel]

Thinking further, I think we do need to consider this a new feature. Because who knows how many places have to be double-checked to actually support such low temperatures. However, it is also a bug since the cfg parameter should have indicated the error up front, rather than emitting an assert alter.

[tkittel]

Yeah, trying to lower the temperature limit from 1K to 0.001K I indeed trigger some assert, in this case related to how many orders of a Taylor expansion is sufficient, which is again dependent on the lowest possible T value.

[tkittel]

Ok, I rewrote (still needs cleanup) all those Taylor expansions now, in a way which will make it more robust and able to work with lower T values. I also made sure that all temperature range checks in the code are now consistently applying 0.001K - 1e6K as the allowed range.

I still have a lot of cleanup + unit tests to fix, and I am also currently investigating a fix for the inadequate beta grid for sampling lowE and lowT (i.e. to get rid of those "triangles" in the beta sample plot above).

[tkittel]

Since I already have a fix for the G1-evaluation issue, I opened a separate issue (#112) for solving artifacts in the beta distribution. A proper fix for the latter is more long-term, while I can already now prepare a release fixing the former.

[tkittel]

So I will try to release NCrystal v3.5.1 later today with a fix for the G1 numerical precision issue, which is indeed the underlying cause of the cross section curves that do not tend to the free scattering cross section value.

For reference, here are a few performance plots with the new release. First the G1 function for Al@5mK which looked like a "delta function" in the v3.5.0 plots pasted in an earlier comment:

Now it looks sane again, and indeed is indentical to the high res reference curve that I have now implemented. The upscattering tail is also OK on a log scale:

(for the record, the G1 function is NOT supposed to go to a delta function as T->0, but Gn on the other hand will become Gaussian as n->infinity).

For more reference, here is also the Al@5mk vdoslux dependency plot redone with v3.5.1, this plot showed a lot of weird inconsistencies above, but is now perfectly sane:

The Be and MgH2 plots from @XuShuqi7 are also OK now (although the MgH2 shows that it benefits well from vdoslux not being too low - which is related to #112 ):

[tkittel]

Release v3.5.1 is out now, with these fixes.

[XuShuqi77]

Thanks! @tkittel With the new version, we can produce correctly the XS at very low temperatures.