[ENHANCEMENT] Look into using `periodictable`
Is your feature request related to a problem? Please describe.
Current MDANSE ships with its own custom JSON tables of atomic properties (c.f. MDANSE/Src/MDANSE/Chemistry/atoms.json), periodictable is a library which packages many of these and is kept up to date with standards.
Describe the solution you'd like
It may be worth investigating this to determine whether it would be possible to replace atoms.json with periodictable entirely.
Describe alternatives you've considered N/A
Additional context N/A
I am looking at the atom properties now, in connection to #713 and #718
For neutron properties, I would normally refer to the NIST tables: https://www.ncnr.nist.gov/resources/n-lengths/list.html According to them, 3He, 6Li and 10B isotopes have complex scattering lengths, both coherent and incoherent.
It seems to me that periodictable provides only the coherent scattering length, and for everything else it only provides the scattering cross sections. (I got this information from https://periodictable.readthedocs.io/en/latest/api/nsf.html ).
In any case, the NIST tables show some complicated values, like for 151Eu where b_incoherent is (+/-)4.5-2.14i, which could not be stored in the current MDANSE atom database either.
Jacob Wilkins
Hi Duc, I've added Franz as the interested party w.r.t. MDMC, and Adam, because I'm not sure whether/how they've dealt with this in Euphonic and he might have some useful input.
So the issue in question here is that some elements (e.g. He-3, Li-6, Gd-155/157, etc.) have imaginary coherent and incoherent scattering factors and we don't know quite how this comes into play or how to use them; in MDMC this has obviously been dealt with incorrectly due to a slight misunderstanding of the periodictable library, in MDANSE we're not 100% sure if we're dealing with them correctly, and 70% sure we're not.
When we were chatting at lunch yesterday a couple of things were mentioned:
- The materials with these imaginary factors are highly absorbing, so are unlikely to be used in experiment
- The factors end up being $b^2$ so you end up with real numbers anyway.
However, based on Keen (https://doi.org/10.1107/S0021889800019993) for e.g. the $G(r)$ the form just uses $b_{i}b_{j}$, which doesn't tally with what we talked about at lunch. The paper also mentions that the $\bar{b}$ is the average over all isotopes and spin states, how do the complex scattering factors get included here? (edited)
Journal of Applied CrystallographyJournal of Applied Crystallography A comparison of various commonly used correlation functions for describing total scattering The formalism that describes total-scattering correlation functions is developed using a consistent notation and is compared with other, frequently encountered, formalisms. It is hoped that this will lead to increased transparency for newcomers to the field of total scattering.
Franz Lang
Just to put some additional context into this. Both MDANSE and MDMC calculate neutron scattering signatures from molecular dynamics trajectories. There are well-established algorithms for these calculations that we implement, for example for the dynamics structure factor from https://doi.org/10.1016/0010-4655(95)00048-K:
[EQUATIONS]
Or for the radial distribution function, as Jacob already mentioned, from https://doi.org/10.1107/S0021889800019993:
[EQUATIONS]
In all of these cases you get these products of scattering lengths. Both MDANSE and MDMC look these up using the periodictable module (https://github.com/python-periodictable/periodictable), which in turn implements the "values collected by the Atomic Institute of the Austrian Universities", which are published at https://www.ati.ac.at/~neutropt/scattering/table.html and contain tables for the real and imaginary ports of the coherent scattering lengths. These look similar (if not identical) to the ones published by NIST at https://www.ncnr.nist.gov/resources/n-lengths/.
This leaves us with the following problem: we have publications for the algorithms to calculate the neutron signatures and we have a Python module that we import to get the values of the scattering lengths that appear in those algorithms. However, we do not know how we should treat complex scattering lengths in these algorithms. Because they should in general lead to complex results, which makes no sense physically.
For some algorithms, like the one for S(Q,w) above, I can see that if you assume that one of the scattering lengths appearing should really be a complex conjugate, then the overall result would end up being a real number due to the symmetry of swapping the labels alpha and beta. But I don't know if that is indeed the correct approach or not.
Interestingly, the ILL publishes the scattering lengths in their "Neutron Data Booklet" (https://www.ill.eu/fileadmin/user_upload/ILL/1_About_ILL/Documentation/ILL_brochure/NeutronDataBooklet.pdf) as the real values only when compared to say the NIST data, but their cross-sections are the same (i.e. the cross sections must include the fact that there is an imaginary part to the scattering length. Which I find confusing.
Duc Le
So I think that you can use the modulus |b| of the complex b scattering lengths in your code. The cross-sections σ should just be 4π b ² and this makes sense if you use either the full complex values or the modulus. The ILL data booklet just lists the real part and I think this is incorrect. I'll try to dig up some literature...
Adam Jackson
If I recall correctly, Abins uses cross sections (from Mantid) directly. Euphonic uses scattering length in an expression where it gets squared
I should check if it is actually brought in as a complex number though
Duc Le [actual answer]
So Eq. 4.21 in Boothroyd is the same as eq (10) in the NMOLDYN paper but with b_alpha complex-conjugated - so I think the NMOLDYN equation is not quite right...
Discussing it with Franz we think that this would give you a real structure factor when the phase factors exp(-iq*R) and the sums alpha->beta and beta->alpha are included... but this won't necessarily work for the radial distribution function...
Right, Boothroyd Chapter 2.4 talks about the pair distribution function and eq 2.43 is the same as the J. Appl. Cryst eq (10)-(11), but again with a complex conjugate on the i term. The $\sin(Qr_{ij}) / Qr_{ij}$ term is real however, which means that if you include both i->j and j->i terms in the sum the imaginary parts cancel and you only have the real parts, effectively $\sum_{ij} \mathcal{R}(b_i)*\mathcal{R}(b_j) \frac{sin(Qr_{ij})}{Qr_{ij}}$ (because rij = rji)
So, I think you should probably use conj(b_i) b_j and complex values for b in all parts of the code and it should give you real outputs... (!)
Franz Lang
I agree. I just need to check now that the various things we calculate all indeed have the same symmetry upon swapping the i,j subscripts. Thank you very much Duc!
Duc Le [actual answer]
Just spotted that you also asked about the averages: The paper also mentions that the $\bar{b}$ is the average over all isotopes and spin states, how do the complex scattering factors get included here?
The $\bar{b}$ values are tabulated as the non-isotope value in Sears' table (reproduced on the NIST site). E.g. for Gd, it is 6.5-13.82i and is just the concentration-weighted average of the scattering lengths of all the isotopes.
The spin states refers to two possible final spin states of the nucleus after interacting with the neutron which is I+1/2 or I-1/2 where I is the nuclear spin quantum number of the specific isotope. The probabilities of the I+1/2 and I-1/2 states are not equal (50/50) but is determined by p+=(I+1)/(2I+1) and p-=I/(2I+1) as described in Boothroyd chapter 1.8.2 (eq. 1.48). The spin probability weighted scattering lengths for each isotope is what is tabulated by NIST/Sears. The individual b+ and b- values are in the ILL data booklet.