Merge pull request #585 from bobmyhill/conventions

added required frame_convention argument to cell parameter conversions

burnman/classes/anisotropicmineral.py

@@ -152,6 +152,11 @@ class AnisotropicMineral(Mineral, AnisotropicMaterial):
     4D array, please refer to the code
     or to the original paper.
 
+    For non-orthotropic materials, the argument frame_convention
+    should be set to define the orientation of the reference frame
+    relative to the crystallographic axes
+    (see :func:`~burnman.utils.unitcell.cell_parameters_to_vectors`).
+
     If the user chooses to define their parameters as a dictionary,
     they must also provide a function to the psi_function argument
     that describes how to compute the tensors Psi, dPsidf and dPsidPth
@@ -188,6 +193,7 @@ def __init__(
         isotropic_mineral,
         cell_parameters,
         anisotropic_parameters,
+        frame_convention=None,
         psi_function=None,
         orthotropic=None,
     ):
@@ -206,6 +212,14 @@ def __init__(
             self.anisotropic_params = anisotropic_parameters
             self.psi_function = psi_function
 
+        if self.orthotropic:
+            self.frame_convention = [0, 1, 2]
+        else:
+            assert (
+                len(frame_convention) == 3
+            ), "frame_convention must be set for this non-orthotropic mineral"
+            self.frame_convention = frame_convention
+
         Psi_Voigt0 = self.psi_function(0.0, 0.0, self.anisotropic_params)[0]
         if not np.all(Psi_Voigt0 == 0.0):
             raise ValueError(
@@ -215,7 +229,9 @@ def __init__(
             )
 
         # cell_vectors is the transpose of the cell tensor M
-        self.cell_vectors_0 = cell_parameters_to_vectors(cell_parameters)
+        self.cell_vectors_0 = cell_parameters_to_vectors(
+            cell_parameters, self.frame_convention
+        )
 
         if (
             np.abs(np.linalg.det(self.cell_vectors_0) - isotropic_mineral.params["V_0"])
@@ -405,7 +421,7 @@ def cell_parameters(self):
             spatial coordinate axes.
         :rtype: numpy.array (1D)
         """
-        return cell_vectors_to_parameters(self.cell_vectors)
+        return cell_vectors_to_parameters(self.cell_vectors, self.frame_convention)
 
     @material_property
     def shear_modulus(self):

burnman/utils/unitcell.py

@@ -25,7 +25,7 @@ def molar_volume_from_unit_cell_volume(unit_cell_v, z):
     return V
 
 
-def cell_parameters_to_vectors(cell_parameters):
+def cell_parameters_to_vectors(cell_parameters, frame_convention):
     """
     Converts cell parameters to unit cell vectors.
 
@@ -37,38 +37,58 @@ def cell_parameters_to_vectors(cell_parameters):
         (:math:`\\gamma`) to the angle between the first and second vectors.
     :type cell_parameters: numpy.array (1D)
 
+    :param frame_convention: A list (c) defining the reference frame
+        convention.  This function dictates that the c[0]th cell vector
+        is colinear with the c[0]th axis, the c[1]th cell vector is
+        perpendicular to the c[2]th axis,
+        and the c[2]th cell vector is defined to give a right-handed
+        coordinate system.
+        In common crystallographic shorthand, x[c[0]] // a[c[0]],
+        x[c[2]] // a[c[2]]^* (i.e. perpendicular to a[c[0]] and a[c[1]]).
+
+    :type frame_convention: list of three integers
+
     :returns: The three vectors defining the parallelopiped cell [m].
-        This function assumes that the first cell vector is colinear with the
-        x-axis, and the second is perpendicular to the z-axis, and the third is
-        defined in a right-handed sense.
+
     :rtype: numpy.array (2D)
     """
-    a, b, c, alpha_deg, beta_deg, gamma_deg = cell_parameters
-    alpha = np.radians(alpha_deg)
-    beta = np.radians(beta_deg)
-    gamma = np.radians(gamma_deg)
-
-    n2 = (np.cos(alpha) - np.cos(gamma) * np.cos(beta)) / np.sin(gamma)
-    M = np.array(
-        [
-            [a, 0, 0],
-            [b * np.cos(gamma), b * np.sin(gamma), 0],
-            [c * np.cos(beta), c * n2, c * np.sqrt(np.sin(beta) ** 2 - n2**2)],
-        ]
+    lengths = cell_parameters[:3]
+    angles_deg = cell_parameters[3:]
+    angles = np.radians(angles_deg)
+
+    c = frame_convention
+
+    n2 = (np.cos(angles[c[0]]) - np.cos(angles[c[2]]) * np.cos(angles[c[1]])) / np.sin(
+        angles[c[2]]
     )
-    return M
+    MT = np.zeros((3, 3))
+    MT[c[0], c[0]] = lengths[c[0]]
+    MT[c[1], c[0]] = lengths[c[1]] * np.cos(angles[c[2]])
+    MT[c[1], c[1]] = lengths[c[1]] * np.sin(angles[c[2]])
+    MT[c[2], c[0]] = lengths[c[2]] * np.cos(angles[c[1]])
+    MT[c[2], c[1]] = lengths[c[2]] * n2
+    MT[c[2], c[2]] = lengths[c[2]] * np.sqrt(np.sin(angles[c[1]]) ** 2 - n2**2)
+
+    return MT
 
 
-def cell_vectors_to_parameters(M):
+def cell_vectors_to_parameters(vectors, frame_convention):
     """
     Converts unit cell vectors to cell parameters.
 
-    :param M: The three vectors defining the parallelopiped cell [m].
-        This function assumes that the first cell vector is colinear with the
-        x-axis, the second is perpendicular to the z-axis, and the third is
-        defined in a right-handed sense.
-    :type M: numpy.array (2D)
+    :param vectors: The three vectors defining the parallelopiped cell [m].
+    :type vectors: numpy.array (2D)
 
+    :param frame_convention: A list (c) defining the reference frame
+        convention. This function dictates that the c[0]th cell vector
+        is colinear with the c[0]th axis, the c[1]th cell vector is
+        perpendicular to the c[2]th axis,
+        and the c[2]th cell vector is defined to give a right-handed
+        coordinate system.
+        In common crystallographic shorthand, x[c[0]] // a[c[0]],
+        x[c[2]] // a[c[2]]^* (i.e. perpendicular to a[c[0]] and a[c[1]]).
+
+    :type frame_convention: list of three integers
 
     :returns: An array containing the three lengths of the unit cell vectors [m],
         and the three angles [degrees].
@@ -79,22 +99,30 @@ def cell_vectors_to_parameters(M):
     :rtype: numpy.array (1D)
     """
 
-    assert M[0, 1] == 0
-    assert M[0, 2] == 0
-    assert M[1, 2] == 0
+    c = frame_convention
+    assert vectors[c[0], c[1]] == 0
+    assert vectors[c[0], c[2]] == 0
+    assert vectors[c[1], c[2]] == 0
 
-    a = M[0, 0]
-    b = np.sqrt(np.power(M[1, 0], 2.0) + np.power(M[1, 1], 2.0))
-    c = np.sqrt(
-        np.power(M[2, 0], 2.0) + np.power(M[2, 1], 2.0) + np.power(M[2, 2], 2.0)
-    )
+    lengths = np.empty(3)
+    angles = np.empty(3)
 
-    gamma = np.arccos(M[1, 0] / b)
-    beta = np.arccos(M[2, 0] / c)
-    alpha = np.arccos(M[2, 1] / c * np.sin(gamma) + np.cos(gamma) * np.cos(beta))
+    lengths[c[0]] = vectors[c[0], c[0]]
+    lengths[c[1]] = np.sqrt(
+        np.power(vectors[c[1], c[0]], 2.0) + np.power(vectors[c[1], c[1]], 2.0)
+    )
+    lengths[c[2]] = np.sqrt(
+        np.power(vectors[c[2], c[0]], 2.0)
+        + np.power(vectors[c[2], c[1]], 2.0)
+        + np.power(vectors[c[2], c[2]], 2.0)
+    )
 
-    gamma_deg = np.degrees(gamma)
-    beta_deg = np.degrees(beta)
-    alpha_deg = np.degrees(alpha)
+    angles[c[2]] = np.arccos(vectors[c[1], c[0]] / lengths[c[1]])
+    angles[c[1]] = np.arccos(vectors[c[2], c[0]] / lengths[c[2]])
+    angles[c[0]] = np.arccos(
+        vectors[c[2], c[1]] / lengths[c[2]] * np.sin(angles[c[2]])
+        + np.cos(angles[c[2]]) * np.cos(angles[c[1]])
+    )
 
-    return np.array([a, b, c, alpha_deg, beta_deg, gamma_deg])
+    angles_deg = np.degrees(angles)
+    return np.array([*lengths, *angles_deg])

tests/test_anisotropicmineral.py

@@ -45,7 +45,7 @@ def make_forsterite(orthotropic=True):
         ]
     )
 
-    m = AnisotropicMineral(fo, cell_parameters, constants)
+    m = AnisotropicMineral(fo, cell_parameters, constants, [0, 1, 2])
     return m
 
 

tests/test_anisotropicsolution.py

@@ -19,7 +19,10 @@ def make_nonorthotropic_mineral(a, b, c, alpha, beta, gamma, d, e, f):
     cell_parameters = np.array(
         [cell_lengths[0], cell_lengths[1], cell_lengths[2], alpha, beta, gamma]
     )
-    fo.params["V_0"] = np.linalg.det(cell_parameters_to_vectors(cell_parameters))
+    frame_convention = [0, 1, 2]
+    fo.params["V_0"] = np.linalg.det(
+        cell_parameters_to_vectors(cell_parameters, frame_convention)
+    )
     constants = np.zeros((6, 6, 3, 1))
     constants[:, :, 1, 0] = np.array(
         [
@@ -42,7 +45,7 @@ def make_nonorthotropic_mineral(a, b, c, alpha, beta, gamma, d, e, f):
         ]
     )
 
-    m = AnisotropicMineral(fo, cell_parameters, constants)
+    m = AnisotropicMineral(fo, cell_parameters, constants, frame_convention)
     return m
 
 

tests/test_anisotropy.py

@@ -4,6 +4,10 @@
 import numpy as np
 
 from burnman.classes import anisotropy
+from burnman.utils.unitcell import (
+    cell_parameters_to_vectors,
+    cell_vectors_to_parameters,
+)
 
 
 class test_anisotropy(BurnManTest):
@@ -135,6 +139,34 @@ def test_rhombohedral_ii(self):
             ],
         )
 
+    def test_cell_parameters_conversions(self):
+        cell_parameters = np.array([1.0, 4.2, 2.2, 80.0, 85.0, 88.0])
+        lengths_orig = cell_parameters[:3]
+        cos_a = np.cos(np.deg2rad(cell_parameters[3:]))
+        for convention in [
+            [0, 1, 2],
+            [0, 2, 1],
+            [1, 0, 2],
+            [1, 2, 0],
+            [2, 1, 0],
+            [2, 0, 1],
+        ]:
+            v = cell_parameters_to_vectors(cell_parameters, convention)
+            p = cell_vectors_to_parameters(v, convention)
+            self.assertArraysAlmostEqual(cell_parameters, p)
+
+            # check angles
+            lengths = np.linalg.norm(v, axis=1)
+            cos_a2 = np.array(
+                [
+                    np.dot(v[1], v[2]) / lengths[1] / lengths[2],
+                    np.dot(v[0], v[2]) / lengths[0] / lengths[2],
+                    np.dot(v[0], v[1]) / lengths[0] / lengths[1],
+                ]
+            )
+            self.assertArraysAlmostEqual(cos_a, cos_a2)
+            self.assertArraysAlmostEqual(lengths_orig, lengths)
+
 
 if __name__ == "__main__":
     unittest.main()