single crystals dir2

[dddijulio]

When using single crystals, we always need to specify dir1 and dir2 What if I don't care about the secondary orientation?

[tkittel]

Yeah, this is something I have wondered if I should make more obvious. Of course, the physics WILL depend on the secondary axis, even if you profess to not caring. The exception is materials like pyrolytic graphite (with lcbragg set), where the need to specify dir2 is truly just a nuisance.

But I also find myself mostly only caring about dir1, which is why the dirtol parameter exists.

So while you DO have to specify dir2, if you set "dirtol=180deg" you can put it to some nonsense (anything non-parallel) and get a more or less arbitrary direction for dir2.

Check "dirtol" on https://github.com/mctools/ncrystal/wiki/CfgRefDoc - the 180deg (ab)usage is noted there 🙂

Additionally, if you would like to somehow get physics representing an average over all choices of dir2, you can set the lcbragg axis to the crystal axis of dir1. However, that might seriously slow you down, and still be unrealistic for most single crystals.

Finally, I have considered in the past if we should make this simpler by supporting dir2=any or dir2=auto, or even just allowing dir2 to be omitted. I will open an issue to track such a feature.

[tkittel]

If you are putting dirtol=180deg you are saying that NCrystal is free to move the secondary axis around as needed, so if you expect dir2 to have an impact, you should not do the dirty trick of using dirtol=180deg (the price is that you need to put a consistent dir2 value, or you will get an error!).

Even with dirtol=180deg, dir2 can still have an impact though. To be more precise, carefully read the dir2 description on https://github.com/mctools/ncrystal/wiki/CfgRefDoc:

Type: crystal axis orientation
Description: Secondary orientation axis of a single crystal. This is specified using the same syntax as for the dir1 parameter. In general the opening angle between the dir1 and dir2 vectors must be nonzero and identical in the crystal and lab frames, but a discrepancy up to the value of the dirtol parameter is allowed. In any case, the components of the dir2 vectors parallel to the dir1 vectors are ignored. When this parameter is set, the parameters mos and dir1 must also be provided.

The crucial point here is that dirtol only says whether or not you will get an error or not (so it doesn't actually change the result, assuming you do not get an error), and the sentence "In any case, the components of the dir2 vectors parallel to the dir1 vectors are ignored".

So if you modify your dir2 value only in the components parallel to dir1, then indeed you will not see an effect. This is also what your plots show, as far as I can tell (I was a bit puzzled, because from your description, I was expecting to see two identical curves in the plot you provided).

If you want to see how the words above are actually applied in code, you can have a look at the verifyLatticeOrientDefAndConstructCrystalRotation function here: https://github.com/mctools/ncrystal/blob/master/ncrystal_core/src/NCLatticeUtils.cc#L415

Specifically look at the last 4-5 lines.

[XuShuqi7]

Hi @tkittel,

thank you for your responses!

So, if I well understood, NCrystal necessitates identical angles between the provided vectors in the crystal and laboratory frames. Then I calculated the secondary direction in order to remove the option dirtol. Three cases are tested using
Cu = NC.createScatter('Cu_sg225.ncmat;temp={}K;{}'.format(T, cfg))
xs = Cu.crossSection(E, (0, 0, 1))
Cua = NC.createAbsorption('Cu_sg225.ncmat;temp={}K;{}'.format(T, cfg))
xsa = Cua.crossSection(E, (0, 0, 1))

The cfgs are indicated in the following figure. Nevertheless, I don't understand why the case 3 does not overlap with one of the others.

[tkittel]

Why would it overlap?

[tkittel]

@XuShuqi7 I now added a print statement to the code to show you what happens in your cfg-strings above (btw., next time please don't just include them in a picture, but please make it easy for me to copy-and-paste the entire cfg-string from text).

case 1 ("Cu_sg225.ncmat;mos=0.05deg;dir1=@crys_hkl:1,1,1@lab:1,1,0;dir2=@crys_hkl:2,0,0@lab:1,1,2"):
NCrystal: REQUESTED ORIENTATION: dir1=@crys_hkl:1,1,1@lab:1,1,0, dir2=@crys_hkl:2,0,0@lab:1,1,2 dirtol=0.0001
NCrystal: FINAL ORIENTATION:     c1={ 0.57735, 0.57735, 0.57735 }, l1={ 0.707107, 0.707107, 0 } c2={ 0.816497, -0.408248, -0.408248 }, l2={ 6.7987e-17, 6.7987e-17, 1 }

case2 ("Cu_sg225.ncmat;mos=0.05deg;dir2=@crys_hkl:1,1,1@lab:1,2,1;dir1=@crys_hkl:2,0,0@lab:1,0,1"):
NCrystal: REQUESTED ORIENTATION: dir1=@crys_hkl:2,0,0@lab:1,0,1, dir2=@crys_hkl:1,1,1@lab:1,2,1 dirtol=0.0001
NCrystal: FINAL ORIENTATION:     c1={ 1, 0, 0 }, l1={ 0.707107, 0, 0.707107 } c2={ 2.71948e-16, 0.707107, 0.707107 }, l2={ 6.7987e-17, 1, 6.7987e-17 }

case3 ("Cu_sg225.ncmat;mos=0.05deg;dir1=@crys_hkl:1,1,1@lab:1,1,0;dir2=@crys_hkl:2,0,0@lab:1,0,1;dirtol=180deg"):
NCrystal: REQUESTED ORIENTATION: dir1=@crys_hkl:1,1,1@lab:1,1,0, dir2=@crys_hkl:2,0,0@lab:1,0,1 dirtol=3.14159
NCrystal: FINAL ORIENTATION:     c1={ 0.57735, 0.57735, 0.57735 }, l1={ 0.707107, 0.707107, 0 } c2={ 0.816497, -0.408248, -0.408248 }, l2={ 0.408248, -0.408248, 0.816497 }

As you can see, you are asking for three different directions (and that is also what it looked like to me).

Also note that in the second of the three cfg-strings, you have swapped the order of dir1 and dir2, I do not know if that was intentional. If I rename dir1 and dir2 in that one, I get:

case2B ("Cu_sg225.ncmat;mos=0.05deg;dir1=@crys_hkl:1,1,1@lab:1,2,1;dir2=@crys_hkl:2,0,0@lab:1,0,1"):
NCrystal: REQUESTED ORIENTATION: dir1=@crys_hkl:1,1,1@lab:1,2,1, dir2=@crys_hkl:2,0,0@lab:1,0,1 dirtol=0.0001
NCrystal: FINAL ORIENTATION:     c1={ 0.57735, 0.57735, 0.57735 }, l1={ 0.408248, 0.816497, 0.408248 } c2={ 0.816497, -0.408248, -0.408248 }, l2={ 0.57735, -0.57735, 0.57735 }

Which again seems different (and again, I am not sure why you would expect any of the cfgs you provided to end up with the same orientation, when clearly they are different).

Cheers,
Thomas

ps. Be aware that today is a national holiday here, so you might not get more replies.

[XuShuqi7]

Hi @tkittel,

thank you very much for your responses! Before I thought that if dirtol=180deg was set, then NCrystal would use one of the two given directions and ignore the other. Now I realize that my thought was wrong. As you commented, I will avoid using dirtol=180deg but define the directions properly.

Thanks also for printing the vectors of final orientation. While I am trying to calculate the rotation matrix based on these vectors by following the function in NCRotMatrix.cc (lines 25-63), because the rotation matrix is more direct for me. In addition, understanding the conversion from dir1=...;dir2=... to the rotation matrix is quite important, because I would like to stack some single crystal blocks in Monte Carlo simulations.

The Python code is

import numpy as np
import sympy

def calcRotMat(v1, u1, v2, u2):
    u2 -= u1 * np.dot(u2, u1)
    v2 -= v1 * np.dot(v2, v1)
    u2 /= np.linalg.norm(u2)
    v2 /= np.linalg.norm(v2)
    
    v1 /= np.linalg.norm(v1)
    u1 /= np.linalg.norm(u1)
    v2 /= np.linalg.norm(v2)
    u2 /= np.linalg.norm(u2)

    print(v1, u1, v2, u2)
    
    cross1 = np.cross(v1, u1)
    cross2 = np.cross(v2, u2)

    aMat = [ [ v2[0], u2[0], cross2[0], 1, 0, 0 ],
             [ v2[1], u2[1], cross2[1], 0, 1, 0 ],
             [ v2[2], u2[2], cross2[2], 0, 0, 1 ] ]

    aMat = sympy.Matrix(aMat).rref()[0]
    aMat = np.array(aMat).astype(np.float64)

    b_matrix = [ [ aMat[0][3], aMat[0][4], aMat[0][5] ],
                 [ aMat[1][3], aMat[1][4], aMat[1][5] ],
                 [ aMat[2][3], aMat[2][4], aMat[2][5] ] ]

    lab_matrix = [ [ v1[0], u1[0], cross1[0] ],
                   [ v1[1], u1[1], cross1[1] ],
                   [ v1[2], u1[2], cross1[2] ] ]

    M = np.matmul(lab_matrix, b_matrix)

    return M

I tested for a trivial case:

v1 = np.array([ 0., 1, 0 ])
u1 = np.array([ 1., 0, 0 ])
v2 = np.array([ 1., 0, 0 ])
u2 = np.array([ 0., 0, 1 ])
M1 = calcRotMat(v1, u1, v2, u2)

M1 is a rotation matrix representing the correct rotation between (u1, u2) and (v1,v2). Neverlethess, when I tested for

v1 = np.array([ 0., 0, 1 ])
u1 = np.array([ 1., 0, 0 ])
v2 = np.array([ -0.5*np.sqrt(2), 0.5*np.sqrt(2), 0 ])
u2 = np.array([ 0.5*np.sqrt(2), 0.5*np.sqrt(2), 0 ])
M2 = calcRotMat(v1, u1, v2, u2)

M2 is not a rotation matrix, but det(M2)=2 and $M_2^TM_2=[[3., 1., 0.],[1., 1., 0.],[0., 0., 2.]])\neq I_3$.

I cannot jump out of the mistake, probably I misunderstand something when reading NCRotMatrix.cc. Do you know what the Python code misses? Thx!

[tkittel]

Opened #156 related to this.