Bugfix for get_points_in_sphere (Superceded by #60)

[pc494]

---

name: Bugfix for get_points_in_sphere
about: Full details in #56

---

Release Notes

major
bugfix
Summary: Internal rotation corections (bundle with #46)

Comments
@EirikOpheim if you could run this on your test case before we merge that would be great.

[EirikOpheim]

Hello Phillip,

I tried to run the program with the suggested changes, but the outcome is the same as before the changes on my computer. I was wondering if converting:

 a, b, c = reciprocal_lattice.a, reciprocal_lattice.b, reciprocal_lattice.c  ----->
 standard_base = diffpy.structure.lattice.Lattice(base=reciprocal_lattice.stdbase)
 a, b, c = standard_base.a, standard_base.b, standard_base.c

actually changes anything? I tried with the rotation matrices for (-60,-60,0) and (30,-60,0) to see if there were any difference in the values of a, b and c between the original and changed implementations. I found them to be the same for both implementations and they are also equal for both rotation matrices.

I also tested the distance function for both reciprocal_lattice and standard_base. It turns out they give the same distance between points of fractional coordinates [0.5,0.5,0.5] and [0.8,0.8,0.8] (Arbitrarily chosen relative coordinates. ). So I believe the bugfix does not actually change how the code runs.

I have to admit that I am not really sure how the process should run. Do we want the standard base to be the direct lattice? I believe standard_base behaves like reciprocal_lattice in every instance where it is called in the function. They are not equal, as reciprocal_lattice.base is not equal to standard_base.base, but their function-specific behavior seems to be equal.

Below is the code I have run to test the functions when orientations (-60,-60,0) and (30,-60,0) still provided the same diffraction pattern.

Regards,
Eirik

Code I have run:

from transforms3d.euler import euler2mat
import numpy as np
import diffpy.structure
from diffsims.generators.diffraction_generator import DiffractionGenerator

orientation1 = (30,-60,0)
orientation2 = (-60,-60,0)

matrix1 = euler2mat(*np.deg2rad(orientation1), 'rzxz')
matrix2 = euler2mat(*np.deg2rad(orientation2), 'rzxz')

latt = diffpy.structure.lattice.Lattice(4.08, 4.08, 4.08, 90, 90, 90)
atom = diffpy.structure.atom.Atom(atype='Au', xyz=[0,0,0], lattice=latt)
atom2 = diffpy.structure.atom.Atom(atype='Au', xyz=[0,1/2,1/2], lattice=latt)
atom3 = diffpy.structure.atom.Atom(atype='Au', xyz=[1/2,0,1/2], lattice=latt)
atom4 = diffpy.structure.atom.Atom(atype='Au', xyz=[1/2,1/2,0], lattice=latt)
structure = diffpy.structure.Structure(atoms=[atom,atom2,atom3,atom4], lattice=latt)

#Lattice 1
stdbase = structure.lattice.stdbase
stdbase_inverse = np.linalg.inv(stdbase)
rotation_matrix_diffpy = stdbase_inverse @ matrix1 @ stdbase
    
lattice_rotated = diffpy.structure.lattice.Lattice(
    *structure.lattice.abcABG(),
    baserot=rotation_matrix_diffpy)
structure_rotated = diffpy.structure.Structure(structure)
structure_rotated.placeInLattice(lattice_rotated)
latt1 = structure_rotated.lattice

#Lattice 2
stdbase = structure.lattice.stdbase
stdbase_inverse = np.linalg.inv(stdbase)
rotation_matrix_diffpy = stdbase_inverse @ matrix2 @ stdbase
    
lattice_rotated = diffpy.structure.lattice.Lattice(
    *structure.lattice.abcABG(),
    baserot=rotation_matrix_diffpy)
structure_rotated = diffpy.structure.Structure(structure)
structure_rotated.placeInLattice(lattice_rotated)
latt2 = structure_rotated.lattice

print(latt1)
print(latt2) #Not equal to latt1

#Lattice 1
recip_latt = latt1.reciprocal()
standard_base = diffpy.structure.lattice.Lattice(base=recip_latt.stdbase)
a, b, c = standard_base.a, standard_base.b, standard_base.c #New Method
d, e, f = recip_latt.a, recip_latt.b, recip_latt.c  #Old Method 
print(a,b,c)
print(d,e,f) # Equal to print(a,b,c)
print(standard_base.dist([0.5,0.5,0.5],[0.8,0.8,0.8]))
print(recip_latt.dist([0.5,0.5,0.5],[0.8,0.8,0.8]))  #Equal to standard_base.dist(..)

#Lattice 2 
recip_latt = latt2.reciprocal()
standard_base = diffpy.structure.lattice.Lattice(base=recip_latt.stdbase)
a, b, c = standard_base.a, standard_base.b, standard_base.c
d, e, f = recip_latt.a, recip_latt.b, recip_latt.c
print(a,b,c)
print(d,e,f) #Equal to print(a,b,c) and equal to print(a,b,c) from Lattice 1. 

[pc494]

I tried to run the program with the suggested changes, but the outcome is the same as before the changes on my computer.

Just to confirm, how far did you run the changes. Did you install this branches code?

I was wondering if converting:

 a, b, c = reciprocal_lattice.a, reciprocal_lattice.b, reciprocal_lattice.c  ----->
 standard_base = diffpy.structure.lattice.Lattice(base=reciprocal_lattice.stdbase)
 a, b, c = standard_base.a, standard_base.b, standard_base.c

actually changes anything? I

It looks like it doesn't, you're right. But there is more than that in the PR :)

I also tested the distance function for both reciprocal_lattice and standard_base. It turns out they give the same distance between points of fractional coordinates [0.5,0.5,0.5] and [0.8,0.8,0.8] (Arbitrarily chosen relative coordinates. ). So I believe the bugfix does not actually change how the code runs.

Here I have to protest! The docstrings for the .distance()method say "in fractional coordinates" and so the two cases should be the same distance apart. You can rotate as much as you like the midpoint and the corner of a cell will be the same distance apart.

I have to admit that I am not really sure how the process should run. Do we want the standard base to be the direct lattice? I believe standard_base behaves like reciprocal_lattice in every instance where it is called in the function. They are not equal, as reciprocal_lattice.base is not equal to standard_base.base, but their function-specific behavior seems to be equal.

standard_base is the reciprocal lattice with (this is the cubic, but the generalisation should be handled by diffpy) a along x, b along y and c along z. This means that using lists of integers (for h,k and l) will behave as we would expect. In another basis (say a 45 degree rotation about z) the 100 spot will lie at cartesians of [1/sqrt(2),1/sqrt(2),0].

The reasons this code is better (I hope) is that now you find that the 100 is present (say) you can use reciprocal_lattice_cartesian() to get the correct location in cartesians.

Here is an MWE

# a = 4, b = 2, c = 1
recip_latt = diffpy.structure.lattice.Lattice(0.25,0.5,1, 90, 90, 90)
standard_base = diffpy.structure.lattice.Lattice(base=recip_latt.stdbase)
a, b, c = standard_base.a, standard_base.b, standard_base.c #New Method
potential_points = np.asarray([(1,0,0)])
in_sphere = np.abs(recip_latt.dist(potential_points, [0, 0, 0])) < 10
spot_indicies = potential_points[in_sphere]
spot_coords = recip_latt.cartesian(spot_indicies)

base_rotated = euler2mat(0,0,np.pi/4,axes='rzxz') @ recip_latt.base
recip_latt_2 = diffpy.structure.lattice.Lattice(base=base_rotated)
a, b, c = standard_base.a, standard_base.b, standard_base.c #New Method
potential_points = np.asarray([(1,0,0)])
in_sphere = np.abs(recip_latt_2.dist(potential_points, [0, 0, 0])) < 10
spot_indicies = potential_points[in_sphere]
spot_coords = recip_latt_2.cartesian(spot_indicies)

which sadly doesn't work (although I think it would for the cubic case).

My proposal is that we instead work entirely in the standard basis, and then at the final step rotate all our cartesians by the rotation matrix. For this to work we need a R that is applied on the left of the cartesian coordinates. I've not figured out what the best way to extract R from diffpy objects is yet.

[EirikOpheim]

Just to confirm, how far did you run the changes. Did you install this branches code?

I am still quite new to git, so I have just copied the changes seen in "files changed" in this PR and added it manually to the file in the diffsims library. If there is an easy way to switch to a specific branch I would be very happy to hear about it! :)

It looks like it doesn't, you're right. But there is more than that in the PR :)

I am a bit confused here. The changes I have seen are five lines added and four lines removed. I wanted to see how these lines would change the logic of which points are selected. I have assumed that moving the local import to global import would not change anything for this function. Furthermore, a,b and c are set to the same numbers. I can only see that the change from reciprocal_base to standard_base is used at one more instance, and that is in calculating the boolean "in-sphere". As the distance function seems to give the same distance (allthough the interpretation is that it is in relative coordinates), I can not see how it will change the behavior of the function. I tried the distance function for relative coordinates against the origin. That is, standard_base.dist([0.3,0.5,0.7],[0,0,0]) and recip_latt.dist([0.3,0.5,0.7],[0,0,0]), as it is used in the function. As far as I can see it still gives the same distance, even when the original lattice is not cubic.

The reasons this code is better (I hope) is that now you find that the 100 is present (say) you can use reciprocal_lattice_cartesian() to get the correct location in cartesians.

I think I understand. It is a utility that we did not have in the old code, but it is not yet in use in the new code?

My proposal is that we instead work entirely in the standard basis, and then at the final step rotate all our cartesians by the rotation matrix. For this to work we need a R that is applied on the left of the cartesian coordinates. I've not figured out what the best way to extract R from diffpy objects is yet.

I think this is a good idea! I will let you know if I should come across a way to extract R from diffpy objects.

[pc494]

In terms of the physics (and code solution) I think #59 is the answer. I'll address the other bits of comments shortly.

[pc494]

Superceded by #60.