Extend MDANSE to support complex neutron scattering lengths
Description of work Users reported problems trying to analyse MD simulations of systems containing 6Li. This PR is meant to enable calculations on atoms with a non-negligible imaginary component of the scattering length.
Fixes
- Added
complexto supported variable types in atom database. - Replaced missing scattering lengths with values from
periodictablewhere possible, and from the NIST table in other cases (typically for the incoherent scattering length). - Removed the
.realpart from the correlation calculations. - Changed the unit of scattering length to fm.
- Added unit information to the atom database.
- Modified the code of the
get_weightsfunction to use the complex conjugate when multiplying b_i*b_j. - Removed b_incoherent2 from the database.
- Updated the TrajectoryWriter to write the new scattering lengths correctly.
- Extended the plotter, which now separates a dataset into two parts (real and imaginary) if it detects a non-zero imaginary component in it.
- Renamed the atomic mass unit from "uma" to "Da".
- Switched from
rffttofftin the jobs that calculate spectra (PPS, CCF, DOS, IR). - Modified the
compare_hdf5function in the tests to compare only the real part of the results. - Performance improvements for reading atom properties from trajectory (by caching).
Note: while FFT is calculated instead of RFFT, the negative side of the results is still discarded in the code.
To test All tests must pass. Convert a trajectory using the protos branch. Convert the same trajectory using this branch. Run DCSF, DISF and DOS analysis for both trajectories. The results should be plottable and nearly identical.
We will need to update the atom database editor as you can now enter any string for their scattering lengths.
@ChiCheng45 I added a ComplexValidator to the atom property table, and replaced all the None values with nan. If #780 gets merged, we could try to add at least a warning in the GUI if the atoms in the trajectory have unknown scattering lengths.
The output data are np.float64 again. Also, assign_weights now takes a complex conjugate of every other pair of atom types, which I thought would be an implementation of the equations you presented above, giving us a $2Re[(\beta_A^{*}\beta_B)]$ prefactor for $F_{AB}$ but please have another look and check if you agree with it.
The output data are
np.float64again. Also,assign_weightsnow takes a complex conjugate of every other pair of atom types, which I thought would be an implementation of the equations you presented above, giving us a 2 R e [ ( β A ∗ β B ) ] prefactor for F A B but please have another look and check if you agree with it.
Not exactly, the equations were
F_{AB}(q,t) = \frac{2}{N\sqrt{c_{A}c_{B}}} \mathrm{Re} \left[\beta_{A}^{*}\beta_{B} C_{AB}(q,t) \right]
I'm not sure that the imaginary part of $C_{AB}(q,t)$ is always small. $C_{AB}(q,t) + C_{BA}(q,t)$ on the other hand should be real.
The output data are
np.float64again. Also,assign_weightsnow takes a complex conjugate of every other pair of atom types, which I thought would be an implementation of the equations you presented above, giving us a 2 R e [ ( β A ∗ β B ) ] prefactor for F A B but please have another look and check if you agree with it.
I had another think about this, and I think this should be fine as long as specific q-vector generators are used. For example, if we used the spherical q-vectors with sufficient sampling, then in our averaging, we should have something like this,
C_{AB}(\vert q \vert, t) = \frac{1}{N_{q}} [\cdots + C_{AB}(-q, t) + C_{AB}(q, t) + \cdots ]
since
C_{AB}(-q, t) = C_{AB}^{*}(q, t) = C_{BA}(q, t)
$C_{AB}(\vert q \vert, t)$ is real and
F_{AB}(\vert q \vert,t) = \frac{2}{N\sqrt{c_{A}c_{B}}} \mathrm{Re} \left[\beta_{A}^{*}\beta_{B} C_{AB}(\vert q \vert,t) \right] = \frac{2}{N\sqrt{c_{A}c_{B}}} \mathrm{Re} \left[\beta_{A}^{*}\beta_{B} \right] C_{AB}(\vert q \vert,t)
The only problem is that the above equation might not be true for other q-vector generation schemes since the sampling may not have generated $q$ and $-q$.
I did some more thinking. This is related to the above comment. Since,
2 \mathrm{Re} [C_{AB}(q, t)] = C_{AB}(q, t) + C_{AB}(-q, t)
so in this branch, we do the following,
\frac{2}{N\sqrt{c_{A}c_{B}}} \mathrm{Re} \left\{\beta_{A}^{*}\beta_{B} \right\} \mathrm{Re} \left\{ C_{AB}(q,t) \right\}= \frac{1}{N\sqrt{c_{A}c_{B}}} \mathrm{Re} \left\{\beta_{A}^{*}\beta_{B} [C_{AB}(q,t) + C_{AB}(-q,t) ] \right\} = \frac{1}{2}[F_{AB}(q,t) + F_{AB}(-q,t)]
in other words when we use any q-vector $q$ we actually average over the values of $q$ and $-q$. I guess there are a few things we could do, update the spherical lattice q-vectors so it always includes $q$ and $-q$, update it so it only includes $q$ or just leave it and rely on the sampling to generate enough $q$ and $-q$.
This has also been solved by #963 instead.