FABLE-3DXRD / ImageD11

ImageD11 is a python code for identifying individual grains in spotty area detector X-ray diffraction images.
https://imaged11.readthedocs.io/
GNU General Public License v2.0
15 stars 25 forks source link

Stress/Strain elastic constants testing and validation #194

Closed jonwright closed 6 days ago

jonwright commented 11 months ago

Moving the discussion in merged pull #192 into a new issue. The problem is which numerical values to input for the numbers in the calculations of stress from strain, and whether the numbers can be checked.

Then for trigonal quartz as a low symmetry example. We had an indexing ambiguity, the (100) direction is not distinguished from the (110) direction until you look at the peak intensities. For high-temperature beta quartz, the structure is hexagonal. For stresses in alpha quartz, it seems to come down to the sign of C14, which goes to zero when the structure goes from trigonal to hexagonal (Ohno, I., Harada, K. & Yoshitomi, C. Temperature variation of elastic constants of quartz across the α - β transition. Phys Chem Minerals 33, 1–9 (2006). https://doi.org/10.1007/s00269-005-0008-3). For crystals that are twinned, this might not matter.

There is a description of the problem of conventions here: https://www.comsol.com/blogs/piezoelectric-materials-understanding-standards/

In the case of quartz, the C14 parameter appears to change sign when going from the IRE 1949 convention to the one from IEEE 1978. Some intensity ratios are discussed here: https://royalsocietypublishing.org/doi/pdf/10.1098/rspa.1926.0025 ... and their relation to macroscopic crystal axes is here: http://www.minsocam.org/ammin/AM30/AM30_326.pdf ... and here, BSTJ 22: 3. October 1943: Use of X-Rays for Determining the Orientation of Quartz Crystals. Bond, W.L.; Armstrong, E.J.: https://archive.org/details/bstj22-3-293

Somehow we need to summarise all that literature into something that matches a set of elastic constants to a crystal structure that gives the diffracted intensities. Maybe taking the numbers from the materials project is unambiguous, as they are providing both together:

https://next-gen.materialsproject.org/materials/mp-7000?crystal_system=Trigonal&formula=SiO2

jadball commented 2 weeks ago

While this is fresh in my mind after a chat with @AxelHenningsson and @jonwright today:

https://github.com/materialsproject/pymatgen/blob/77195153effda5dc97e28a47507ab64a52246ffb/src/pymatgen/core/tensors.py#L406

# Take the B0 of the cell
C_c = np.linalg.inv(B0).T
E  = np.column_stack((C_c[:, 0], np.cross(C_c[:, 2], C_c[:, 0]), C_c[:, 2]))
E /= np.linalg.norm(E, axis=0)

strain_crystal_tensor_IEEE =  E.T @ ( strain_crystal_tensor_imaged11 ) @ E
jonwright commented 2 weeks ago

The materials project api is giving me back two elastic tensors, one "raw" and one "ieee_format". Also a set of lattice vectors in the structure field (with a surprising selection of axes).

The code above looks like it does not work for monoclinic or triclinic? C_c look like they are meant to be the a,b,c real space vectors, then it is forming x||a, z||c and making a right-handed set with axc || b* along y. The matrix E is not a rotation unless a and c are at right angles. I suspect we need something like:

z || c y || axc x || (axc)xc

so:

a,b,c = np.linalg.inv(B0)
bstar = np.cross( c, a )
E = np.transpose( ( np.cross( bstar, c ), bstar, c ) )
E /= np.linalg.norm(W, axis=0)
xfab.checks._check_rotation_matrix( E )

... but it all depends on what the IEEE people have defined in their top-secret document.

I have a vague feeling we would be better off just defining the IEEE reference lattice when computing strains. For example:

dzero_cell = ImageD11.grain.grain( ubi = [ ieee_a_vector, ieee_b_vector, ieee_c_vector ] )
strain_crystal_tensor_IEEE = grain.eps_grain_matrix(self, dzero_cell, m=0.5)

Like this, the strains ought to be with respect to the IEEE reference state (or we have to go back to debugging that finite strain code again). It was meant to be set up to let you input the reference as a vector lattice.

jadball commented 2 weeks ago

It seems the IEEE rules are defined here: https://github.com/materialsproject/pymatgen/blob/77195153effda5dc97e28a47507ab64a52246ffb/src/pymatgen/core/tensors.py#L406

To use the IEEE reference lattice when computing strains, wouldn't we have to rotate the B matrix of grain to match too? If we can get that worked out, I can see it working. @AxelHenningsson what do you think?

jonwright commented 1 week ago

It seems that elastic constants need to come with a reference "lattice" setting. I wonder if it makes more sense to add a property to unitcell that gives back unitcell.IEEE for a reference lattice in the IEEE convention. This leaves it open to add some other convention when it comes up, and lets you raise an error for unit cells that need a transformation (primitive to rhombohedral, or a,b,c not in the right order, etc). Several important unit cell settings that are allowed by ImageD11 are not permitted by IEEE, so we would have to transform the UBI to the new setting as well.

For the case of quartz, the problem seems to be the sign of C14. Materials project has both mp-6390 and mp-7000, and they have opposite signs for C14. To get this right, we need to compute the structure factors using the two crystal structures that are in materials project and cross-check with their atom co-ordinates. The quartz lattice has a 6-fold axis along c, but the structure only has a 3-fold. The IEEE rotation they are returning looks like a 6 fold, so they are flipping the sign of their C14. The scatter on the numbers (raw vs IEEE and mp-6930 vs mp-7000) suggests there is some other difference or problem. But I don't think you can compute stresses for quartz unless you have checked the intensity of (101) vs (011), etc, and picked the orientation with an extra indexing step that uses intensities.

Demo code to get data from materials project is below. Something I don't understand is that the unit cell is completely wrong (5.1 vs 4.9 for a). Anyway, the "lattice" matrix they return looks like the same one as ImageD11, but they are choosing a [110] direction. So their "raw" setting looks like the one from Busing and Levy, but with a 3-fold axis applied.

import requests, pprint, numpy as np, ImageD11.unitcell,  pymatgen.core,  pymatgen.core.tensors

headers = {'accept': 'application/json', 
           'X-API-KEY': 'YourKeyHere' }

for mpid in 6930, 7000:
    url = f'https://api.materialsproject.org/materials/elasticity/?material_ids=mp-{mpid}&_per_page=100&_skip=0&_limit=100&_fields=structure%2Celastic_tensor&_all_fields=false'
    req = requests.get(url, headers=headers)
    # print( req.content )
    ans = eval( req.content.decode().replace('true','True') )[ 'data' ] # for demo only
    assert len(ans) == 1
    ans = ans[0]
    print("Lattice (for raw elastic constants)")
    pprint.pprint( ans['structure']['lattice'] )
    mpstr = pymatgen.core.Structure( lattice=ans['structure']['lattice']['matrix'],
                         species=[atom['species'][0]['element'] for atom in ans['structure']['sites']],
                         coords=[atom['abc'] for atom in ans['structure']['sites']] )
    print(mpstr)
    print( 'Raw elastic constants: \n%s'%(str(np.round(np.array( ans['elastic_tensor']['raw'] )))))
    print( 'IEEE symmetrized: \n%s'%(str(np.array( ans['elastic_tensor']['ieee_format'] ))))
    print( 'pymatgen.core.tensors.Tensor.get_ieee_rotation:\n%s'%(str(
          pymatgen.core.tensors.Tensor.get_ieee_rotation( s ))))
    uc = ImageD11.unitcell.unitcell( [ ans['structure']['lattice'][axis] for axis in 
                                     ['a','b','c','alpha','beta','gamma' ] ] )
    UB0 = uc.B
    ubi_0 = np.linalg.inv( UB0 ) # rows
    print("Busing and Levy reference lattice: \n%s"%(str(ubi_0)))
    print("a+b from ImageD11:", ubi_0[0] + ubi_0[1] )

and output:

Lattice (for raw elastic constants)
{'a': 5.024354086551485,
 'alpha': 90.0,
 'b': 5.024354086551485,
 'beta': 90.0,
 'c': 5.51356489,
 'gamma': 120.00000061297857,
 'matrix': [[2.51217702, -4.35121829, -0.0],
            [2.51217702, 4.35121829, 0.0],
            [0.0, -0.0, 5.51356489]],
 'pbc': [True, True, True],
 'volume': 120.53789302383238}
Full Formula (Si3 O6)
Reduced Formula: SiO2
abc   :   5.024354   5.024354   5.513565
angles:  90.000000  90.000000 120.000001
pbc   :       True       True       True
Sites (9)
  #  SP           a         b          c
---  ----  --------  --------  ---------
  0  Si    0.522818  0.522818  -0
  1  Si    0.477182  0          0.666667
  2  Si    0         0.477182   0.333333
  3  O     0.585256  0.839586   0.870772
  4  O     0.160414  0.745671   0.537438
  5  O     0.254329  0.414744   0.204105
  6  O     0.839586  0.585256   0.129228
  7  O     0.745671  0.160414   0.462562
  8  O     0.414744  0.254329   0.795895
Raw elastic constants: 
[[ 93.  -3.   8. -12.  -0.  -0.]
 [ -3.  93.  10.   8.  -0.  -0.]
 [  8.  10.  86.  -1.  -0.  -0.]
 [-12.   8.  -1.  62.  -0.   0.]
 [ -0.  -0.  -0.  -0.  50. -17.]
 [ -0.  -0.  -0.   0. -17.  45.]]
IEEE symmetrized: 
[[ 91.  -1.   9.  13.  -0.  -0.]
 [ -1.  91.   9. -13.   0.  -0.]
 [  9.   9.  86.  -0.  -0.  -0.]
 [ 13. -13.  -0.  56.  -0.   0.]
 [ -0.   0.  -0.  -0.  56.  13.]
 [ -0.  -0.  -0.   0.  13.  46.]]
pymatgen.core.tensors.Tensor.get_ieee_rotation:
[[ 0.49999994 -0.86602544  0.        ]
 [ 0.86602544  0.49999994 -0.        ]
 [ 0.          0.          1.        ]]
Busing and Levy reference lattice: 
[[ 4.35121825e+00 -2.51217709e+00  1.51206732e-15]
 [ 0.00000000e+00  5.02435409e+00  3.07652957e-16]
 [ 0.00000000e+00  0.00000000e+00  5.51356489e+00]]
a+b from ImageD11: [4.35121825e+00 2.51217700e+00 1.81972028e-15]
Lattice (for raw elastic constants)
{'a': 5.019645263159945,
 'alpha': 90.0,
 'b': 5.019645263159945,
 'beta': 90.0,
 'c': 5.50894152,
 'gamma': 120.00000794971314,
 'matrix': [[2.50982233, -4.34714049, 0.0],
            [2.50982233, 4.34714049, -0.0],
            [0.0, -0.0, 5.50894152]],
 'pbc': [True, True, True],
 'volume': 120.21116681490265}
Full Formula (Si3 O6)
Reduced Formula: SiO2
abc   :   5.019645   5.019645   5.508942
angles:  90.000000  90.000000 120.000008
pbc   :       True       True       True
Sites (9)
  #  SP            a          b         c
---  ----  ---------  ---------  --------
  0  Si     0.523687   0.523687  0
  1  Si    -0          0.476313  0.666667
  2  Si     0.476313  -0         0.333333
  3  O      0.256054   0.4148    0.794569
  4  O      0.5852     0.841253  0.127902
  5  O      0.158747   0.743946  0.461236
  6  O      0.4148     0.256054  0.205431
  7  O      0.743946   0.158747  0.538764
  8  O      0.841253   0.5852    0.872098
Raw elastic constants: 
[[ 87. -10.   9.  12.   0.   0.]
 [-10.  83.   6. -13.   0.   0.]
 [  9.   6.  94.  -0.   0.   0.]
 [ 12. -13.  -0.  57.   0.   0.]
 [  0.   0.   0.   0.  56.  15.]
 [  0.   0.   0.   0.  15.  44.]]
IEEE symmetrized: 
[[ 83.  -9.   7. -14.   0.  -0.]
 [ -9.  83.   7.  14.  -0.  -0.]
 [  7.   7.  94.   0.   0.  -0.]
 [-14.  14.   0.  57.  -0.  -0.]
 [  0.  -0.   0.  -0.  57. -14.]
 [ -0.  -0.  -0.  -0. -14.  46.]]
pymatgen.core.tensors.Tensor.get_ieee_rotation:
[[ 0.49999994 -0.86602544  0.        ]
 [ 0.86602544  0.49999994 -0.        ]
 [ 0.          0.          1.        ]]
Busing and Levy reference lattice: 
[[ 4.34713997e+00 -2.50982323e+00  1.51065005e-15]
 [ 0.00000000e+00  5.01964526e+00  3.07364625e-16]
 [ 0.00000000e+00  0.00000000e+00  5.50894152e+00]]
a+b from ImageD11: [4.34713997e+00 2.50982203e+00 1.81801468e-15]