Behavior of Structure.Bisoequiv #51

@sbillinge I am porting over a discussion from the diffpy-users google group on the behavior of Structure.Bisoequiv for calculating PDFs via real-space summation with diffpy.srreal.pdfcalculator.PDFCalculator.

The original thread:

Hello all,

I have posted questions about calculating scattering from periodic (MD) models before, but didn't get much traction. Here I'm attaching a small example of the issues I've been encountering in the form of a jupyter notebook in hopes that this will be more useful in addressing my specific concerns.

In this example, I want to calculate the scattering from a single frame of a MD simulation of water. The frame is from a trajectory of a periodic water box with dimension 34.691 Angstroms. First, I calculate the scattering in reciprocal space treating the frame as an aperiodic structure using the DebyePDFCalculator. Despite ignoring periodic boundary conditions, the calculation 'looks good' to me, in that it looks like what I would expect from a total scattering measurement on water. Moreover, all of the individual pair correlations (O-O, O-H, and H-H correlations) make sense in real and reciprocal space.

However, when I try to repeat this calculation using real-space summation with the PDFCalculator, I get huge oscillations in the PDF baseline. Back-transforming the component G(r)'s back to reciprocal space, it's clear that somehow this method fails to properly treat the low-Q region for the O-O and H-H pairs, which produce the oscillations in the baseline (erroneous sharp peaks at low-Q). Moreover, the relative amplitude/sharpness of the O-H pair correlation seems incorrect. Even though the O-H distances are constrained in the structure, so the peak should be a delta-function in real space, I feel like the amplitude is nonetheless 'too sharp' given what should be the average scattering power for this pair.

I've been stumped for quite some time as to why trying to do a periodic calculation in real-space produces these issues. Any insights would be greatly appreciated! I will note that, anecdotally, doing similar real-space summations with models that don't contain H atoms do not seem to suffer from these pathologies.

Thanks,
Nik

Hi Nik

I think you will be a little bit annoyed by the answer ;)
There is a bug in the PDFcalculator. I think it also exists in the Debye version as well.

The structure.Bisoequiv = 0.04 line does not really set the thermal parameters for your model.
Instead you should use:

structure[structure.element == 'Element which you want to set Uiso for'].U = 0.005

If you prefer I think you can use B by doing it like this: structure[structure.element == 'Element which you want >to set Biso for'].B = 0.005

If you add the lines to the top of your calculation:

structure[structure.element == 'H'].U = 0.005

structure[structure.element == 'O'].U = 0.005

the ripples go away.

I would also use this for the Debye calculator as atoms are always moving no matter how you calculate the scattering.

I have attached the curves I get with those two lines in your script.

Cheers Mikkel

Thanks so much Mikkel (and Nik),

I posted this as an issue on the diffpy.structure repo on GH so we can look into it. Please can you check that it is actually a bug though? I wonder if we are giving enough information to diffpy.structure for it to know what to update with structure.Bisoequiv = 0.04. Is it ambiguous which elements' Uiso's should be updated? Of course, this is resolved by your code.

What behavior would you want to happen if the user types structure.Bisoequiv = 0.04? That ALL ADPs are updated to that value? This may be too much "magic" and lead to unintended results for the user if, in their head, they are updating just one atom type, let's say, but the program is doing something else. What do you think?

My usual preference is less magic as it is harder to debug. Clearly, one issue here was perhaps (I didn't check) that the ADPs weren't set but no message was returned to the user that this was the case?

S

Hi Simon

You are right that it might not be a bug, but just a documentation issue. I ran into this issue because I too thought `structure.Bisoequiv could be used with the PDFcalculator. After some trial and error I found this solution. So maybe updating the description is enough to solve the issue.

Cheers Mikkel

Hi Simon and Mikkel,

Thanks for your insights. I am still a little confused about the documentation on 'Structure.Bisoequiv' still. I.e.,

help(Structure.Bisoequiv)
Help on property:
Array of Debye-Waller isotropic thermal displacement or equivalent values. Assignment updates the U attribute of all atoms.

And when I test this:

structure = loadStructure('test.xyz')
structure.U
array([[[0., 0., 0.], [0., 0., 0.], [0., 0., 0.]], ...
structure.Bisoequiv = 0.04
structure.U
array([[[0.00050661, 0. , 0. ], [0. , 0.00050661, 0. ], [0. , 0. , 0.00050661]], ...

It certainly does appear that updating Bisoequiv for a structure object puts thermal parameters on all atoms in the structure object. Nevertheless, Mikkel, your solution addresses, somewhat, the odd behavior for the O-H correlations using the PDFCalculator.

However, I still do not understand, for example, the discrepancy between calculating the O-O correlations with the two methods. The partial PDF from the Debye calculator looks correct, but the real-space summation produces unphysical (and by the looks of it, undamped) oscillations in the PDF. Is there some reason the real-space calculator doesn't behave properly when using setTypeMask?

Cheers,
Nik

To be clear, I think there may be two issues at hand: (1) the expected behavior of Structure.Bisoequiv for setting thermal parameters of atoms in Structure objects, and (2) the expected behavior of using setTypeMask with a PDFCalculator object. the latter might more properly be an issue to file under diffpy.srreal.