Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Update WK_scattering_factors.py #700

Merged

Conversation

ZacharyVarley
Copy link
Contributor

@ZacharyVarley ZacharyVarley commented Dec 8, 2024

Relevant WK computation in others.f90 in Prof. Marc De Graef's EMsoft reads:

DO J=1,4
  FPHON = FPHON + A1(J)*A1(J)*(DEWA * RI1(B1(J),B1(J),G) - RI2(B1(J),B1(J),G,UL))
  DO I=1,J-1
   FPHON = FPHON + 2.*A1(J)*A1(I)*(DEWA * RI1(B1(I),B1(J),G)-RI2(B1(I),B1(J),G,UL))

When J=1:
The inner loop condition is DO I=1,0 (since J-1 = 0). In Fortran, when the upper bound is less than the lower bound, the loop is not executed at all. So for J=1, only the first FPHON calculation happens, the inner loop is skipped. In the current py4DSTEM code, the inner loop is executed for ii = jj = 0 making the diagonal computed a second time with a factor of 2.0, for an overall prefactor of 3.0 instead of 1.0 when ii equals jj. If I am missing something and this python code has already been verified against others.f90, then completely ignore this.

Relevant WK computation in others.f90 in Prof. Marc De Graef's EMsoft reads: 

DO J=1,4
  FPHON = FPHON + A1(J)*A1(J)*(DEWA * RI1(B1(J),B1(J),G) - RI2(B1(J),B1(J),G,UL))
  DO I=1,J-1
   FPHON = FPHON + 2.*A1(J)*A1(I)*(DEWA * RI1(B1(I),B1(J),G)-RI2(B1(I),B1(J),G,UL))

When J=1:
The inner loop condition is DO I=1,0 (since J-1 = 0)
In Fortran, when the upper bound is less than the lower bound, the loop is not executed at all
So for J=1, only the first FPHON calculation happens, the inner loop is skipped. In the current py4DSTEM code, the inner loop is executed for ii = jj = 0 making the diagonal be computed a second time with a factor of 2.0 and take an overall prefactor of 3.0 instead of 1.0. If I am missing something and this has been verified against others.f90, then completely ignore this.
@sezelt
Copy link
Member

sezelt commented Dec 10, 2024

Good catch on the discrepancy. It seems we do differ from the relevant snippet in EMSoftLib:
https://github.com/EMsoft-org/EMsoft/blob/2c071d6975265a423b1e83cf5436a8995293c1d3/Source/EMsoftLib/others.f90#L210-L236

I haven't yet had a chance to test this update, but I think we can also fix this more cleanly by just fixing the range limits:

for jj in range(4):
    for ii in range(jj):
        [...]

Curiously, during initial implementation I used to have yet another different bounds (both of which probably resulted from misunderstanding the Fortran behavior): a74d2fb#diff-652c64a84248c41f4cc70b79f26dd9b2f72223f4079c9a52668e886ebc2bf820R190

@ZacharyVarley
Copy link
Contributor Author

I have written a completely vectorized version (over accelerating voltage, Z & sigma pairs per atom type, and g) in PyTorch and I will try to port it to numpy / scipy for you all at some point. It might be worth putting some sort of disclaimer for those who have used these WK computations thus far as the delta is significant.

@sezelt sezelt merged commit 966a41c into py4dstem:dev Dec 18, 2024
6 checks passed
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

Successfully merging this pull request may close these issues.

2 participants