mantid
mantid copied to clipboard
mantidproject
Reworked Abins DMOL3 parser
Description of work.
Addresses #34341 : the DMOL3 parser currently scrambles the eigenvector matrix when the data takes more than one block of 9 columns to present. (i.e. when more than 9 eigenvectors are printed). This data is used to determine intensities in Abins. In some cases the resulting error may not be obvious, as frequencies will still be correct but peak intensities are not.
The problem was detected when a user file happened to yield a large block of zeros when read incorrectly, causing a divide-by-zero error in a normalisation step that should never see a total of zero.
There are a couple of other odd features of the parser that can be addressed at the same time:
- a normalisation constant that is stored as state on the parser and will keep accumulating if the file is read repeatedly. It is unclear if this does any long-term harm, as another normalisation step follows. But there should not be state changes resulting from repeated reader calls from the same file: that is surprising behaviour.
- the second normalisation uses the default Mantid masses, but these may not be consistent with those used in the eigenvector calculation? It would make more sense to normalise by the reported masses, which reflect any isotopic substitutions.
The root cause of the eigenvector problem is that it was difficult to inspect the output eigenvectors and compare them with the input file, because they were re-ordered and normalised at read-time. The new implementation should separate the reading and post-processing so that it is easier to sanity-check the intermediate results. By making the testing more atomic we need rely less on file-based tests to cover all the edge cases.
Report to: Stewart Parker
To test: The reworked parser should have better unit tests than the original version.
The acid test would be comparison with similar data from another code such as CASTEP, to be posted in this Issue.
Fixes #34341
Reviewer
Please comment on the following (full description):
Code Review
- Is the code of an acceptable quality?
- Does the code conform to the coding standards?
- Are the unit tests small and test the class in isolation?
- If there is GUI work does it follow the GUI standards?
- If there are changes in the release notes then do they describe the changes appropriately?
- Are the release notes saved in a separate file, using Issue or PR number for file name and in the correct location?
Functional Tests
- Do changes function as described? Add comments below that describe the tests performed?
- Do the changes handle unexpected situations, e.g. bad input?
- Has the relevant (user and developer) documentation been added/updated?
Does everything look good? Mark the review as Approve. A member of @mantidproject/gatekeepers will take care of it.
Update: I think the new parser is a lot more testable and understandable. It returns eigenvectors that are different to the existing ones, as expected. I would like to take a bit more time to really sanity-check these results against the other data formats. Unfortunately I will be on annual leave for the next week; I appreciate this reduces the likelihood of getting the fix into release-next. We shall see what happens!
Hi @ajjackson,
Do we think we might be able to get this fix in? If it helps at all I've just removed the need to support RHEL7 (python 3.6) so if we rebase this on release-next that will be one less thing have to worry about.
That's great news! I'll run the rebase in a minute then.
Executive summary: code and unit tests already seem good to go, but validation is lacking
I am currently rather conflicted on this fix: on the one hand, the new parser is clearly better-written, better-tested and removes a demonstrable problem.
But empirically the improvement is less obvious. For the small molecule (methane) case which I had previously believed to be unaffected by the scrambling, the DMOL3 results seem to get worse compared to a CRYSTAL17 calculation of the same system. After a day or two digging through the numbers I am fairly convinced that the small-molecule case was also wrong and should also be fixed by the new parser; yet empirically it appears that somehow it has got worse.
This may be a coincidence related to other differences in the results between calculators; "fortuitous cancellation of errors" being painfully undone. From the end of November I will begin host a graduate project that aims to generate a larger test set so that this kind of thing can be investigated much more systematically.
Or perhaps there is a mistake in the new code.
I want to avoid the situation that we announce a "fix" in release that immediately needs revising in nightly or the next point release. The resulting confusion could be just as -or more- harmful than the bug, as it might encourage more people to use a broken feature...
I will try to think of another validation test that can quickly be used to clarify the situation.
I think I've got to the bottom of it. The sqrt(mass) factor applied to the DMOL3 eigenvectors, which I imitated from the existing parser, should not be there. This did have some kind of fortuitous cancellation of errors with the eigenvector scrambling, but now that both things are fixed we get convincing agreement with a CRYSTAL17 calculation for the same system.
Of course the frequency difference remains, but the intensities now match the CRYSTAL calculation quite closely. By visual inspection the traces of the displacement tensors (i.e. b_traces in Abins) are also similar to the results from CASTEP. (There's no point in including the CASTEP results on this plot though as that calculation is dominated by spurious rotational modes which were removed in the open-boundary codes.)
Executive summary: I'm happy with the maths for now
I have run another test calculation and inspected the data and am now fairly confident that
- the raw data from DMOL3 are normalised eigenvectors (also known as "polarisation vectors")
- it is correct to convert these to atomic displacements by dividing these through by sqrt(mass) as Abins does (indirectly) in powdercalculator.py
This puts them on a similar footing to the data from CASTEP, which also does not apply mass weighting in the parser.
CRYSTAL is a bit different; the raw data are atomic displacements with the amplitude of the zero-point vibrational mode. The parser converts these to the expected normalised polarisation vectors. (The exact means by which it does this remains mysterious to me, but is no longer today's problem.)
I need to update the systemtest data, perhaps that can wait until tomorrow while Anthony checks the code?
Here is the systemtest mismatch. Some difference is expected in intensities, while frequencies should match identically. If anything, I'm surprised there is not more difference.
The reference file will be updated with the new output.
There is actually systemtest data for the same system calculated with CASTEP. Making a three-way comparison... we can maybe argue that the new parser intensities agrees better in the high frequencies, but it's a weak claim. The difference in calculated frequencies between CASTEP and DMOL (not the fault of the parsers...) is bigger than the difference between DMOL parsers in the low-mid frequency range.
I hope to make some more rigorous cross-code comparisons next year, working with an SCD graduate. Until then this is the best we have.
Aha, the intensity improvement is more convincing when we simplify the spectrum by only plotting the order-1 Fluorine atom contributions
Release note added, are there any issues making sure it gets into the formatted docs at this stage?
The release note will be added to the appropriate page manually as part of final cleanup by the release editors once the PR is merged in. I'm not imaging there will be any issues but we'll be in touch if anything crops up.
