dkriegner / xrayutilities

xrayutilities - a package with useful scripts for X-ray diffraction physicists
http://xrayutilities.sourceforge.io
GNU General Public License v2.0
81 stars 29 forks source link

Miller indices in hexagonal systems #104

Closed Pej-ecp closed 4 years ago

Pej-ecp commented 4 years ago

I am interested in the effect of the rhombohedral angle on the intensity of some super-structure peak. I have therefore defined my crystal (sp.group R3m, #160) in hexagonal setting, and made a loop to define crystals (with xu.materials.SGlattice(160, a_H, c_H), atoms=[...], pos=[...], occ=[...])) with the values of the lattice parameters a_H and c_H calculated from a_R (held constant) and alpha (varying between 90 and 89 degrees). I then generate the PowderDiffraction pattern for each with xu.simpack.Powder(...) and create a list of [alpha, xu.simpack.PowderDiffraction(...)] So far so good.

When I try to plot the evolution of the intensity of a given reflection vs alpha, say (212) I get the error message: KeyError: (2, 1, 2) because they key is absent from some patterns (but not all, it does exist in some as verified with individual prints of each PowderDiffraction element).

On further inspection, it turns out that the chosen Miller indices for (2 1 2) are (3 -1 2) in some patterns. The Miller-Bravais indices corresponding to (2 1 2) being (2 1 -3 2), the equivalent directions should be (-3 2 2) and (1 -3 2), not (3 -1 2). This happens for several reflections, e.g. (0 3 0) becomes (3 0 0), (2 1 8) becomes (3 -1 8), etc. but, strangely some reflections remain unchanged ex: (4 -1 4), (4 3 2), (6 1 4)... as if the reference frame was changing (from (a_H, b_H, c_H) to (-b_H, -a_H, c_H) for some reflections.

I therefore have two questions: 1/ is this a bug or I am missing something ? 2/ is there a way for me to access these reflections other than by hand (e.g. through a selection of q or 2theta range) ?

Thanks a lot for your help (I can pm a python script or minimal working example if needed)

dkriegner commented 4 years ago

I am afraid the behavior you describe is how the code works (at the moment). The names of reflections in PowderDiffraction can change when the lattice parameters are different. I, however, want to point out that always an equivalent set of Miller indices is used. For spacegroup 160, 212 and 3-12 are equivalent. See Crystal.equivalent_hkls((2,1,2)).

To get what you want I see two options: 1) define a series of PowderDiffraction objects as you do and then check every equivalent of the Miller indices you look for. This should work but might be cumbersome. 2) define only one PowderDiffraction object and change the lattice parameters in a loop after the definition and recalculate the powder lines after every change. In this way the names of reflections stay the same since the naming is only done once (during initialization of the object)

To change for example 'a' do

p.mat.a = 5 p.set_sample_parameters() p.update_powder_lines()

'p' must be your PowderDiffraction object.

I would not recommend to try to access the reflections via q-range. This will be confusing if you change lattice parameters.

On Mon, 10 Aug 2020, 19:43 Pej-ecp, notifications@github.com wrote:

I am interested in the effect of the rhombohedral angle on the intensity of some super-structure peak. I have therefore defined my crystal (sp.group R3m, #160) in hexagonal setting, and made a loop to define crystals (with xu.materials.SGlattice(160, a_H, c_H), atoms=[...], pos=[...], occ=[...])) with the values of the lattice parameters a_H and c_H calculated from a_R (held constant) and alpha (varying between 90 and 89 degrees). I then generate the PowderDiffraction pattern for each with xu.simpack.Powder(...) and create a list of [alpha, xu.simpack.PowderDiffraction(...)] So far so good.

When I try to plot the evolution of the intensity of a given reflection vs alpha, say (212) I get the error message: KeyError: (2, 1, 2) because they key is absent from some patterns (but not all, it does exist in some as verified with individual prints of each PowderDiffraction element).

On further inspection, it turns out that the chosen Miller indices for (2 1 2) are (3 -1 2) in some patterns. The Miller-Bravais indices corresponding to (2 1 2) being (2 1 -3 2), the equivalent directions should be (-3 2 2) and (1 -3 2), not (3 -1 2). This happens for several reflections, e.g. (0 3 0) becomes (3 0 0), (2 1 8) becomes (3 -1 8), etc. but, strangely some reflections remain unchanged ex: (4 -1 4), (4 3 2), (6 1 4)... as if the reference frame was changing (from (a_H, b_H, c_H) to (-b_H, -a_H, c_H) for some reflections.

I therefore have two questions: 1/ is this a bug or I am missing something ? 2/ is there a way for me to access these reflections other than by hand (e.g. through a selection of q or 2theta range) ?

Thanks a lot for your help (I can pm a python script or minimal working example if needed)

— You are receiving this because you are subscribed to this thread. Reply to this email directly, view it on GitHub https://github.com/dkriegner/xrayutilities/issues/104, or unsubscribe https://github.com/notifications/unsubscribe-auth/ACKZJFKWNM3PENKV6F27KE3SAAWVBANCNFSM4P2GDZFQ .

Pej-ecp commented 4 years ago

Thanks a lot @dkriegner for this quick answer, I do appreciate.

I have defined the following function

def pd_vs_alpha(pd, a_r, alpha_min, nb_int=10, alpha_max=90):
    """
    Generate a list of [alpha, PowderDiff] for Rhombohedral cells

    Parameters
    ----------
    pd : PowderDiffraction object
        PowderDiffraction object generated from reference structure.
    a_r : float
        Rhombohedral angle [A]
    alpha_min : float
        minimum value of alpha [degrees]
    nb_int : integer, optional 
       nb of intervals (i.e. nb_pts-1) to be calculated 
    alpha_max : float, optional
        maximum value of alpha [degrees]. The default is 90.

    Returns
    -------
    list of [alpha, PowderDiff]
    """
    pdlist = []
    for i in range(0,nb_int+1):
        alpha = alpha_max-i*(alpha_max-alpha_min)/nb_int
        alpha_rad = np.deg2rad(alpha)
        a_h = 2*a_r*np.sin(alpha_rad/2)
        c_h = a_r*np.sqrt(3)*np.sqrt(1+2*np.cos(alpha_rad))
        print(alpha, a_h, c_h)
        pd.mat.a = a_h
        pd.mat.c = c_h
        pd.set_sample_parameters()
        pd.update_powder_lines(tt_cutoff=tt_max)
        pdlist.append([alpha, pd])
    return pdlist

However, the DiffractionPowder elements in my list (the pds) are strictly identical: same lattice parameters, same intensities...

I thought that pd.update_powder_lines(tt_cutoff=tt_max) changed (as in "replaced") the PowderDiffraction object. Is it not the case ? The Docstring only mentions that it calculates the powder intensity and positions up to an angle of tt_cutoff (deg) and updates the values.

Side point about Miller-Bravais indices: when I try Crystal.equivalent_hkls((2,1,2)), I get AttributeError: 'Crystal' object has no attribute 'equivalent_hkls'. I have tried with the xu.materials.Crystal object created without more success.

dkriegner commented 4 years ago

Here you seem to not understand how Python works. In your function you always add a reference of the PowderDiffraction object to the return list. In every iteration it is changed. But the reference still points to the one and only PowderDiffraction object you work with. So it must be the same. If you really need the different objects you will need to copy them. You can try the Python copy module (copy.deepcopy) but I am not sure it will work with PowderDiffraction objects.

Regarding the second issue: Did you try

pd.mat.lattice.equivalent_hkls((2,1,2)) ?

On Tue, 11 Aug 2020, 11:58 Pej-ecp, notifications@github.com wrote:

Thanks a lot @dkriegner https://github.com/dkriegner for this quick answer, I do appreciate.

I have defined the following function

def pd_vs_alpha(pd, a_r, alpha_min, nb_int=10, alpha_max=90): """ Generate a list of [alpha, PowderDiff] for Rhombohedral cells

Parameters
----------
pd : PowderDiffraction object
    PowderDiffraction object generated from reference structure.
a_r : float
    Rhombohedral angle [A]
alpha_min : float
    minimum value of alpha [degrees]
nb_int : integer, optional
   nb of intervals (i.e. nb_pts-1) to be calculated
alpha_max : float, optional
    maximum value of alpha [degrees]. The default is 90.

Returns
-------
list of [alpha, PowderDiff]
"""
pdlist = []
for i in range(0,nb_int+1):
    alpha = alpha_max-i*(alpha_max-alpha_min)/nb_int
    alpha_rad = np.deg2rad(alpha)
    a_h = 2*a_r*np.sin(alpha_rad/2)
    c_h = a_r*np.sqrt(3)*np.sqrt(1+2*np.cos(alpha_rad))
    print(alpha, a_h, c_h)
    pd.mat.a = a_h
    pd.mat.c = c_h
    pd.set_sample_parameters()
    pd.update_powder_lines(tt_cutoff=tt_max)
    pdlist.append([alpha, pd])
return pdlist

However, the DiffractionPowder elements in my list (the pds) are strictly identical: same lattice parameters, same intensities...

I thought that pd.update_powder_lines(tt_cutoff=tt_max) changed (as in "replaced") the PowderDiffraction object. Is it not the case ? The Docstring only mentions that it calculates the powder intensity and positions up to an angle of tt_cutoff (deg) and updates the values.

Side point about Miller-Bravais indices: when I try Crystal.equivalent_hkls((2,1,2)), I get AttributeError: 'Crystal' object has no attribute 'equivalent_hkls'. I have tried with the xu.materials.Crystal object created without more success.

— You are receiving this because you were mentioned. Reply to this email directly, view it on GitHub https://github.com/dkriegner/xrayutilities/issues/104#issuecomment-671851411, or unsubscribe https://github.com/notifications/unsubscribe-auth/ACKZJFJQNQAEBV5ZCXYWA7DSAEI5BANCNFSM4P2GDZFQ .

Pej-ecp commented 4 years ago

I acknowledge having still many things to learn about Python... Thank you for your response and your time. I tried several things before responding but I am still stuck.

I am trying to understand how update_powder_lines works in order to use it properly. Looking at its definition (https://xrayutilities.sourceforge.io/_modules/xrayutilities/simpack/powder.html#PowderDiffraction.update_powder_lines), I had the impression that it changed (replaced) the values in the PowderDiffraction object for existing peaks and created new ones if need be. Have I got that right ?

Therefore I have tried a simple loop, not within a function; to avoid confusions between local and global variables:

p = xu.simpack.Powder(ABCx, 1, crystallite_size_gauss=cryst_size)
pd = xu.simpack.PowderDiffraction(p, tt_cutoff=tt_max)

pdlist2 = []
alpha_min, alpha_max = 60, 90
nb_int=10
a_r=8.16
for i in range(0,11):
    alpha = alpha_max-i*(alpha_max-alpha_min)/nb_int
    alpha_rad = np.deg2rad(alpha)
    a_h = 2*a_r*np.sin(alpha_rad/2)
    c_h = a_r*np.sqrt(3)*np.sqrt(1+2*np.cos(alpha_rad))
    pd.mat.a = a_h
    pd.mat.c = c_h
    pd.set_sample_parameters()
    pd.update_powder_lines(tt_cutoff=tt_max)
    pdlist2.append([alpha, pd])

after having defined ABCx as a xu.materials.Crystal object. This gives identical PowderDiffraction objects for each element of pdlist2. Am I again missing something ?

Regarding the second issue: There is no pd.mat.lattice (AttributeError: 'Powder' object has no attribute 'lattice'), but I have found the pd.mat.material.lattice that only has the pd.mat.material.lattice.isequivalent to check whether two hkl Miller indices are equivalent or not. I could not find on x-ray utilities documentation any references to "equivalent_hkls".

dkriegner commented 4 years ago

That's not an issue of global vs local variables. There is simply only one PowderDiffraction object and in every loop iteration this same object is changed. You just look at the same PowderDiffraction object in your list multiple times!

The equivalent_hkls function is very new and not yet included in version 1.6.0. so its only present in git master. Its also possible that the peak naming issue we discussed before manifests itself different in the git master since this part of the code was rewritten after the last release. But i think as a user one should not rely on the naming to be the same for different lattice parameters. Otherwise one needs to define how to select THE name among the equivalent Miller indices. I do not know if a common rule for this exists...

On Wed, 12 Aug 2020, 16:52 Pej-ecp, notifications@github.com wrote:

I acknowledge having still many things to learn about Python... Thank you for your response and your time. I tried several things before responding but I am still stuck.

I am trying to understand how update_powder_lines works in order to use it properly. Looking at its definition ( https://xrayutilities.sourceforge.io/_modules/xrayutilities/simpack/powder.html#PowderDiffraction.update_powder_lines), I had the impression that it changed (replaced) the values in the PowderDiffraction object for existing peaks and created new ones if need be. Have I got that right ?

Therefore I have tried a simple loop, not within a function; to avoid confusions between local and global variables:

p = xu.simpack.Powder(ABCx, 1, crystallite_size_gauss=cryst_size) pd = xu.simpack.PowderDiffraction(p, tt_cutoff=tt_max)

pdlist2 = [] alpha_min, alpha_max = 60, 90 nb_int=10 a_r=8.16 for i in range(0,11): alpha = alpha_max-i(alpha_max-alpha_min)/nb_int alpha_rad = np.deg2rad(alpha) a_h = 2a_rnp.sin(alpha_rad/2) c_h = a_rnp.sqrt(3)np.sqrt(1+2np.cos(alpha_rad)) pd.mat.a = a_h pd.mat.c = c_h pd.set_sample_parameters() pd.update_powder_lines(tt_cutoff=tt_max) pdlist2.append([alpha, pd])

after having defined ABCx as a xu.materials.Crystal object. This gives identical PowderDiffraction objects for each element of pdlist2. Am I again missing something ?

Regarding the second issue: There is no pd.mat.lattice (AttributeError: 'Powder' object has no attribute 'lattice'), but I have found the pd.mat.material.lattice that only has the pd.mat.material.lattice.isequivalent to check whether two hkl Miller indices are equivalent or not. I could not find on x-ray utilities documentation any references to "equivalent_hkls".

— You are receiving this because you were mentioned. Reply to this email directly, view it on GitHub https://github.com/dkriegner/xrayutilities/issues/104#issuecomment-672920339, or unsubscribe https://github.com/notifications/unsubscribe-auth/ACKZJFPKAE2DW3KUNE6V2GDSAKUDTANCNFSM4P2GDZFQ .

dkriegner commented 4 years ago

I now tried this again once myself and have two comments:

1) my previous statement regarding the Python problems of the code sent by @Pej-ecp remain the same.

2) I originally made a slightly wrong comment regarding xrayutilities PowderDiffraction. I said that always an equivalent set of Miller indices is used to name the peak. This is not always true. For the purpose of the calculation of the powder diffraction signal for isotropic systems its not needed to know the "name" of the peak and therefore peaks at the same q-position are merged to one line. Especially in non-centrosymmetric space groups (like 160) one might want to see the difference between e.g. 003 and 00-3.

Therefore if one wants to track the intensity of a specific set of Miller indices one needs to actually force the behavior of an anisotropic PowderDiffraction class. This can be easiest done like this:

import xrayutilities as xu

class FPani(xu.simpack.FP_profile):
    isotropic = False

m = xu.materials.Crystal('test', xu.materials.SGLattice(160, 4, 5))
print(xu.simpack.PowderDiffraction(m, fpclass=FPani))

The code above then lists peaks at equivalent q separately:

Reflections: 
--------------
      h k l     |    tth    |    |Q|    |Int     |   Int (%)
   ---------------------------------------------------------------
    (-1, 0, -1)    31.3903      2.207         0.01      100.00
      (1, 0, 1)    31.3903      2.207         0.01      100.00
    (0, -1, -2)    44.6647      3.099         0.00       45.95
      (0, 1, 2)    44.6647      3.099         0.00       45.95
    (-1, -1, 0)    45.3059      3.142         0.01       88.99
     (0, 0, -3)    55.0555      3.770         0.00        9.52
      (0, 0, 3)    55.0555      3.770         0.00        9.52
....

likely the code given above needs the latest GIT master version to work as shown.

I would propose to close this issue. If there are technical issues after fixing the Python problems in the reporters code separate issues should be filed

Pej-ecp commented 4 years ago

To complement @dkriegner's comment I have tested the following things:

- to track the intensity of a given reflection as a function of the lattice parameter: generate a series of PowderDiffraction objects and get the list of the Miller indices used through hkl_pd = [list(el) for el in [*pd[0].data.keys()]] with pd[0] any of the generated PowderDiffraction object (as they are identical since generated from the same Powder object). For each PowderDiffraction object, update the lattice parameters with pd[i].mat.a = ..., pd[i].mat.c = ... , and pd[i].set_sample_parameters(), then update PowderDiffraction object with pd[i].update_powder_lines(tt_cutoff=tt_max) as suggested by @dkriegner . This enables to refer to given peak(s) in a consistent manner despite the change of lattice parameter. @dkriegner comment above remains absolutely relevant for peaks with different Miller indices but same q value.

- about the (hkl) Miller indices chosen from the set of (hkil) Bravais-Miller indices (obtained through circular permutation of +/- (hki)), I confirm that I could not find any official recommendation as to which to chose. The solution I have used to determine which Miller indices (hkl) are used in a PowderDiffraction pattern that correspond to the reflection I am interested in in the rhombohedral setting is : 1) to calculate from the Miller indices in R setting the Bravais-Miller indices (hkil) with i=-h-k in hexagonal setting, 2) to generate the 6 possible Miller indices from circular permutations and change of signs ( (hkl) , (-h-kl), (kil), (-k-il) etc. ), and 3) to search for these in the list of hkl values used in the PowderDiffraction object data.keys(). 4) to use the (hkl) values in H setting throughout the rest of the script.

I would like to point out that @dkriegner comment is of great use (as usual) when one wishes to convert the Miller indices in H setting back into the R setting as e.g. (306) and (036) in H setting correspond to the same q value, whereas they correspond e.g. to (112) and (330) in R setting.

Unless something's wrong in this comment, I agree to close that issue.

dkriegner commented 4 years ago

what @Pej-ecp writes seems a valid approach to me.

I would only add two more things: for space group 160 in hexagonal setting: 306 and 036 are indeed not the same. When using git master the equivalent_hkls function also shows that:

import xrayutilities as xu
m = xu.materials.Crystal('test', xu.materials.SGLattice(160, 4,5))
m.lattice.equivalent_hkls((3,0,6))
# output: {(-3, 3, 6), (0, -3, 6), (3, 0, 6)}
m.lattice.equivalent_hkls((0,3,6))
# output: {(-3, 0, 6), (0, 3, 6), (3, -3, 6)}

One thing which also could be done if the hexagonal setting is actually a burden: use the rhombohedral setting directly:

m = xu.materials.Crystal('test', xu.materials.SGLattice('160:R', 4, 65))

But I would assume that the other problems described here remain the same.