Stress/Strain elastic constants testing and validation

[jonwright]

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.

• Some unit tests could be added.

• we do not have the old FitAllB strain convention in ImageD11 yet. This could be added for checking against fitallb.

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]

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

• Stiffness tensors have a coordinate system, which relates to how the grain Cartesian reference frame is aligned with the lattice vectors (x along a*, etc etc). As far as I can tell, Materials Project (via pymatgen) always provides the Stiffness Tensor in the IEEE 1978 reference:

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

• The IEEE 1978 system is different from the one used in ImageD11
  Axel has worked this out and determined an additional tensor rotation to apply to ensure agreement:

# 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

• The stiffness tensor will change depending on what vector notation system you use for the strain tensor. Therefore, there are "Voigt" stiffness tensors (overwhelmingly common, used on Materials Project) and "Mandel" stiffness tensors (much less common). This means I have some changes to make in ImageD11/sinograms/tensor_map.py, and that we should clearly document which format we expect the stiffness tensor to be in

• Eventually I would like to merge/combine my Numba-accelerated kernels with stress.py, rather than specifying that we are working on flat Python lists of UBI objects which is currently the case

[jonwright]

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]

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]

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]