Revisiting DebyePDFCalculator vs PDFCalculator

[Niklas Thompson]

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

[Mikkel Juelsholt]

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

[Simon Billinge]

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

--
You received this message because you are subscribed to the Google Groups "diffpy-users" group.
To unsubscribe from this group and stop receiving emails from it, send an email to diffpy-users...@googlegroups.com.
To view this discussion on the web visit https://groups.google.com/d/msgid/diffpy-users/d0e986f9-af6f-436c-b4e9-e099f9dad144n%40googlegroups.com.

--

Simon Billinge

Professor, Columbia University

Physicist, Brookhaven National Laboratory

[Mikkel Juelsholt]

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 

[Niklas Thompson]

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

[Simon Billinge]

It may be a bug then.   We should maybe take this off this thread and move to a discussion on the issue at diffpy.structure on GitHub, both so we don't spam everyone and also for teh sake of good documentation.   I didn't end up opening an issue.  Nik, please could you post your bug report as an issue there and (and also cut/paste the rest of the discussion from here) and we can continue over there (https://github.com/diffpy/diffpy.structure/issues)?

Thanks so much.

S

To view this discussion on the web visit https://groups.google.com/d/msgid/diffpy-users/6deda54d-3f80-4e34-ba31-322f3c1253adn%40googlegroups.com.