Behavior of Structure.Bisoequiv

[niklasbt]

@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.

[Juelsholt]

Hi Nik

I have a possible answer to your 2nd issue, but it is just a guess and I cannot find out if it is real. But it might be that PDFcalculator strictly only takes periodic structures, aka crystal structures. So the long-range correlations you see are because the calculator repeats your input file over and over again. However, I cannot see what type of file types the PDFcalculator takes if there are any limits at all. I know this could be an issue if we were in the fitting part of Diffpy. But I am not sure about the calculators.

[niklasbt]

Hi Mikkel,

The thing is, nominally this structure is periodic. It was taken from an MD simulation in a periodic box, hence the Structure.placeInLattice command in my original example.

In fact, you can see the 'real' periodic image effect if you set rmax to some value larger than the nominal cell size. One of the reasons I have been comparing the reciprocal vs. real space summation approaches is because I was curious what the difference would be between treating these periodic calculations as discrete versus extended structures.

Moreover, the spurious oscillations are only apparent when you calculate the partial PDFs. Somehow these oscillations 'cancel out' when you sum the partial PDFs to the total, i.e., the oscillations in g(H-H) almost cancel the oscillations in g(O-H) to create a 'normal' looking PDF for water.

[sbillinge]

I saw that anticorrelation of the oscillatrions but I need to think what is an explanation. The other thing going on is what we have been discussing...the Uiso's seemed to be not behaving the same way in each calculator. This may be the bug, but needs a bit more looking.

[Juelsholt]

Correct me if I am wrong but in MD simulations periodic boundary conditions "moves" atoms at an edge to the opposite edge. However in Diffpy periodic boundaries creates another cell adjacent to the original cell. Since we have a static situation there will be some correlation over the border between neighbouring unit cells. See crude and fast illustration.

But I have no idea if I am correct, just guessing.

[niklasbt]

I think that is a correct description of periodic boundary conditions in MD calculations. However, the forces acting on the particles do take into account an 'infinite' tiling of the simulation cell, in that atoms at the boundary experience forces due to periodic images from their nearest neighbors. Thus, (I believe) the arrangement of atoms in the cell should be consistent with how Diffpy calculates the pair correlations, if that makes sense.

I can calculate the O-O RDF using, e.g. VMD, and get the following with and without periodic conditions enforced in the calculation.

Without PBC:

With PBC:

These were calculated from exactly the same frame as in my Diffpy example. As you can see, there are no undamped oscillations in the baseline when using VMD, which is why I think the behavior of the PDFCalculator is spurious.

[niklasbt]

Here's a side-by-side comparison of the three calculation methods for the partial O-O PDF:

[pavoljuhas]

structure.Bisoequiv = 0.04

I did not have time to follow the discussion in detail, but note that B=0.04 is a very small value, it is usually about 10 times larger. The relation between B and U is just a scaling factor of 8*pi^2. The mean square displacement for a 2 atom separation at B=0.04 is sqrt(2*U) = 0.03, which is just ~3 datapoints on PDFCalculator calculation grid spaced at 0.01;
not enough to represent peak function.

Perhaps that could be the cause of inaccuracies, try to adjust the B value or, if correct, use a smaller rstep in PDFCalculator.

[sbillinge]

Thanks Pavol, good catch.

…

On Wed, Aug 31, 2022 at 4:18 PM Pavol Juhas ***@***.***> wrote: structure.Bisoequiv = 0.04 I did not have time to follow the discussion in detail, but note that B=0.04 is a very small value, it is usually about 10 times larger. The relation between B and U is just a scaling factor of 8*pi^2. The mean square displacement for a 2 atom separation at B=0.04 is sqrt(2*U) = 0.03, which is just ~3 datapoints on PDFCalculator calculation grid spaced at 0.01; not enough to represent peak function. Perhaps that could be the cause of inaccuracies, try to adjust the B value or, if correct, use a smaller rstep in PDFCalculator. — Reply to this email directly, view it on GitHub <https://urldefense.proofpoint.com/v2/url?u=https-3A__github.com_diffpy_diffpy.structure_issues_51-23issuecomment-2D1233379919&d=DwMCaQ&c=009klHSCxuh5AI1vNQzSO0KGjl4nbi2Q0M1QLJX9BeE&r=6Pin9WIW2wJZzVMqXib07fokbkHzCKPr15sA4fQTrA8&m=mziszj3YF6HIci7Wzwn88SECLubJiY6cNzfeRdje0GGH5MGe5C7jD4IwQpAxpegu&s=6bzd_tX06QYQBHE1uCyzrnsSTcdZfVz9V9hXOlKpXlI&e=>, or unsubscribe <https://urldefense.proofpoint.com/v2/url?u=https-3A__github.com_notifications_unsubscribe-2Dauth_ABAOWUP6KVQ35QB6YAJY5N3V364ZBANCNFSM576V6LLA&d=DwMCaQ&c=009klHSCxuh5AI1vNQzSO0KGjl4nbi2Q0M1QLJX9BeE&r=6Pin9WIW2wJZzVMqXib07fokbkHzCKPr15sA4fQTrA8&m=mziszj3YF6HIci7Wzwn88SECLubJiY6cNzfeRdje0GGH5MGe5C7jD4IwQpAxpegu&s=r9KejEo-C8QfHKF6Nvj8FPhWhVKEwnJtjSjp4sk-P7o&e=> . You are receiving this because you were mentioned.Message ID: ***@***.***>

-- Simon Billinge Professor, Columbia University Physicist, Brookhaven National Laboratory

[niklasbt]

Hi Pavol,

I don't think the thermal parameters being too small is creating the oscillations in the partial PDFs. E.g., here is a series of calculations on the same structure:

For a given rmax, the oscillations seem invariant to setting the U value for the oxygens over three orders of magnitude. Rather, rmax seems to control the 'frequency' of the oscillation in real space.