diff --git a/burnman/classes/anisotropicmineral.py b/burnman/classes/anisotropicmineral.py index 99bfa947..fa7805ec 100644 --- a/burnman/classes/anisotropicmineral.py +++ b/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): diff --git a/burnman/utils/unitcell.py b/burnman/utils/unitcell.py index 34a8e2fc..112e4107 100644 --- a/burnman/utils/unitcell.py +++ b/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]) diff --git a/tests/test_anisotropicmineral.py b/tests/test_anisotropicmineral.py index 9048a975..65a39705 100644 --- a/tests/test_anisotropicmineral.py +++ b/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 diff --git a/tests/test_anisotropicsolution.py b/tests/test_anisotropicsolution.py index 922ebbf4..13e8514c 100644 --- a/tests/test_anisotropicsolution.py +++ b/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 diff --git a/tests/test_anisotropy.py b/tests/test_anisotropy.py index 38e8860e..8bf5e408 100644 --- a/tests/test_anisotropy.py +++ b/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()