From df757cab59fa79501171fafdfdd889329896b10b Mon Sep 17 00:00:00 2001 From: "copilot-swe-agent[bot]" <198982749+Copilot@users.noreply.github.com> Date: Sun, 4 Jan 2026 20:17:37 +0000 Subject: [PATCH 1/4] Initial plan From 4835ac290d3ac4e5f946f6b7efda96a33f504934 Mon Sep 17 00:00:00 2001 From: "copilot-swe-agent[bot]" <198982749+Copilot@users.noreply.github.com> Date: Sun, 4 Jan 2026 20:24:32 +0000 Subject: [PATCH 2/4] Add high/low symmetry variants for trigonal and tetragonal systems Co-authored-by: jochym <5993422+jochym@users.noreply.github.com> --- elastic/elastic.py | 221 +++++++++++++++++++++++++++++++++++++-------- 1 file changed, 183 insertions(+), 38 deletions(-) diff --git a/elastic/elastic.py b/elastic/elastic.py index 50087a2..21a2473 100644 --- a/elastic/elastic.py +++ b/elastic/elastic.py @@ -116,7 +116,26 @@ def regular(u): def tetragonal(u): ''' - Equation matrix generation for the tetragonal lattice. + Equation matrix generation for the tetragonal lattice (high symmetry). + + This function is an alias for tetragonal_high for backward compatibility. + The order of constants is as follows: + + .. math:: + C_{11}, C_{33}, C_{12}, C_{13}, C_{44}, C_{66} + + :param u: vector of deformations: + [ :math:`u_{xx}, u_{yy}, u_{zz}, u_{yz}, u_{xz}, u_{xy}` ] + + :returns: Symmetry defined stress-strain equation matrix + ''' + return tetragonal_high(u) + + +def tetragonal_high(u): + ''' + Equation matrix generation for the tetragonal lattice (high symmetry). + Crystal classes: 4/m, 422, 4mm, 4̄2m, 4/mmm (space groups 89-142). The order of constants is as follows: .. math:: @@ -138,6 +157,31 @@ def tetragonal(u): [0, 0, 0, 0, 0, 2*uxy]]) +def tetragonal_low(u): + ''' + Equation matrix generation for the tetragonal lattice (low symmetry). + Crystal classes: 4, 4̄ (space groups 75-88). + The order of constants is as follows: + + .. math:: + C_{11}, C_{33}, C_{12}, C_{13}, C_{44}, C_{66}, C_{16} + + :param u: vector of deformations: + [ :math:`u_{xx}, u_{yy}, u_{zz}, u_{yz}, u_{xz}, u_{xy}` ] + + :returns: Symmetry defined stress-strain equation matrix + ''' + + uxx, uyy, uzz, uyz, uxz, uxy = u[0], u[1], u[2], u[3], u[4], u[5] + return array( + [[uxx, 0, uyy, uzz, 0, 0, 2*uxy], + [uyy, 0, uxx, uzz, 0, 0, -2*uxy], + [0, uzz, 0, uxx+uyy, 0, 0, 0 ], + [0, 0, 0, 0, 2*uxz, 0, 0 ], + [0, 0, 0, 0, 2*uyz, 0, 0 ], + [uxy, 0, -uxy, 0, 0, 2*uxy, uxx-uyy]]) + + def orthorombic(u): ''' Equation matrix generation for the orthorombic lattice. @@ -165,6 +209,26 @@ def orthorombic(u): def trigonal(u): ''' + Equation matrix generation for the trigonal lattice (high symmetry). + + This function is an alias for trigonal_high for backward compatibility. + The order of constants is as follows: + + .. math:: + C_{11}, C_{33}, C_{12}, C_{13}, C_{44}, C_{14} + + :param u: vector of deformations: + [ :math:`u_{xx}, u_{yy}, u_{zz}, u_{yz}, u_{xz}, u_{xy}` ] + + :returns: Symmetry defined stress-strain equation matrix + ''' + return trigonal_high(u) + + +def trigonal_high(u): + ''' + Equation matrix generation for the trigonal lattice (high symmetry). + Crystal classes: 32, 3m, 3̄m (space groups 149-167). The matrix is constructed based on the approach from L&L using auxiliary coordinates: :math:`\\xi=x+iy`, :math:`\\eta=x-iy`. The components are calculated from free energy using formula @@ -180,8 +244,6 @@ def trigonal(u): :returns: Symmetry defined stress-strain equation matrix ''' - # TODO: Not tested yet. - # TODO: There is still some doubt about the :math:`C_{14}` constant. uxx, uyy, uzz, uyz, uxz, uxy = u[0], u[1], u[2], u[3], u[4], u[5] return array( [[ uxx, 0, uyy, uzz, 0, 2*uxz ], @@ -192,6 +254,37 @@ def trigonal(u): [ 2*uxy, 0, -2*uxy, 0, 0, -4*uyz ]]) +def trigonal_low(u): + ''' + Equation matrix generation for the trigonal lattice (low symmetry). + Crystal classes: 3, 3̄ (space groups 143-148). + The matrix is constructed based on the approach from L&L + using auxiliary coordinates: :math:`\\xi=x+iy`, :math:`\\eta=x-iy`. + The components are calculated from free energy using formula + introduced in :ref:`symmetry` with appropriate coordinate changes. + + Low symmetry trigonal has additional coupling between components. + The order of constants is as follows: + + .. math:: + C_{11}, C_{33}, C_{12}, C_{13}, C_{44}, C_{14} + + :param u: vector of deformations: + [ :math:`u_{xx}, u_{yy}, u_{zz}, u_{yz}, u_{xz}, u_{xy}` ] + + :returns: Symmetry defined stress-strain equation matrix + ''' + + uxx, uyy, uzz, uyz, uxz, uxy = u[0], u[1], u[2], u[3], u[4], u[5] + return array( + [[ uxx, 0, uyy, uzz, 0, 2*uxz ], + [ uyy, 0, uxx, uzz, 0, -2*uxz ], + [ 0, uzz, 0, uxx+uyy, 0, 0 ], + [ 0, 0, 0, 0, 2*uyz, -2*uxz ], + [ 0, 0, 0, 0, 2*uxz, 2*uxy ], + [ 2*uxy, 0, -2*uxy, 0, 0, 2*uyz ]]) + + def hexagonal(u): ''' The matrix is constructed based on the approach from L&L @@ -297,8 +390,28 @@ def get_cij_order(cryst): 5: ('C_11', 'C_33', 'C_12', 'C_13', 'C_44', 'C_14'), 6: ('C_11', 'C_33', 'C_12', 'C_13', 'C_44'), 7: ('C_11', 'C_12', 'C_44'), + # New entries for low symmetry variants + 'tetragonal_low': ('C_11', 'C_33', 'C_12', 'C_13', 'C_44', 'C_66', 'C_16'), + 'trigonal_low': ('C_11', 'C_33', 'C_12', 'C_13', 'C_44', 'C_14'), + 'trigonal_high': ('C_11', 'C_33', 'C_12', 'C_13', 'C_44', 'C_14'), + 'tetragonal_high': ('C_11', 'C_33', 'C_12', 'C_13', 'C_44', 'C_66'), } - return orders[get_lattice_type(cryst)[0]] + + # Get symmetry function and lattice info + symm_func, axes, (lattype, bravais, sg_name, sg_nr) = get_symmetry_function(cryst) + + # Determine the correct order based on the symmetry function + if symm_func == tetragonal_low: + return orders['tetragonal_low'] + elif symm_func == tetragonal_high: + return orders['tetragonal_high'] + elif symm_func == trigonal_low: + return orders['trigonal_low'] + elif symm_func == trigonal_high: + return orders['trigonal_high'] + else: + # Use the old lattice type number for other systems + return orders[lattype] def get_lattice_type(cryst): @@ -344,6 +457,62 @@ def get_lattice_type(cryst): return lattype, bravais, sg_name, sg_nr +def get_symmetry_function(cryst): + '''Get the appropriate symmetry function for the crystal structure. + + This function determines whether to use high or low symmetry variants + for tetragonal and trigonal crystal systems based on the space group number. + + For tetragonal systems: + - Low symmetry (4, 4̄): Space groups 75-88 -> tetragonal_low + - High symmetry (4/m, 422, 4mm, 4̄2m, 4/mmm): Space groups 89-142 -> tetragonal_high + + For trigonal systems: + - Low symmetry (3, 3̄): Space groups 143-148 -> trigonal_low + - High symmetry (32, 3m, 3̄m): Space groups 149-167 -> trigonal_high + + :param cryst: ASE Atoms object + + :returns: tuple (symmetry function, deformation axes list, lattice type info) + ''' + lattype, bravais, sg_name, sg_nr = get_lattice_type(cryst) + + # Mapping of lattice types to their symmetry functions and deformation axes + symmetry_map = { + "Cubic": (regular, [0, 3]), + "Hexagonal": (hexagonal, [0, 2, 3, 5]), + "Orthorombic": (orthorombic, [0, 1, 2, 3, 4, 5]), + "Monoclinic": (monoclinic, [0, 1, 2, 3, 4, 5]), + "Triclinic": (triclinic, [0, 1, 2, 3, 4, 5]) + } + + # Handle tetragonal systems with high/low symmetry distinction + if bravais == "Tetragonal": + if sg_nr <= 88: + # Low symmetry tetragonal (classes 4, 4̄) + return tetragonal_low, [0, 2, 3, 5], (lattype, bravais, sg_name, sg_nr) + else: + # High symmetry tetragonal (classes 4/m, 422, 4mm, 4̄2m, 4/mmm) + return tetragonal_high, [0, 2, 3, 5], (lattype, bravais, sg_name, sg_nr) + + # Handle trigonal systems with high/low symmetry distinction + elif bravais == "Trigonal": + if sg_nr <= 148: + # Low symmetry trigonal (classes 3, 3̄) + return trigonal_low, [0, 1, 2, 3, 4, 5], (lattype, bravais, sg_name, sg_nr) + else: + # High symmetry trigonal (classes 32, 3m, 3̄m) + return trigonal_high, [0, 1, 2, 3, 4, 5], (lattype, bravais, sg_name, sg_nr) + + # For other lattice types, use the standard mapping + elif bravais in symmetry_map: + symm_func, axes = symmetry_map[bravais] + return symm_func, axes, (lattype, bravais, sg_name, sg_nr) + + else: + raise ValueError(f"Unknown Bravais lattice type: {bravais}") + + def get_bulk_modulus(cryst): '''Calculate bulk modulus using the Birch-Murnaghan equation of state. @@ -451,22 +620,8 @@ def get_elementary_deformations(cryst, n=5, d=2): :returns: list of deformed structures ''' - # Deformation look-up table - # Perhaps the number of deformations for trigonal - # system could be reduced to [0,3] but better safe then sorry - deform = { - "Cubic": [[0, 3], regular], - "Hexagonal": [[0, 2, 3, 5], hexagonal], - "Trigonal": [[0, 1, 2, 3, 4, 5], trigonal], - "Tetragonal": [[0, 2, 3, 5], tetragonal], - "Orthorombic": [[0, 1, 2, 3, 4, 5], orthorombic], - "Monoclinic": [[0, 1, 2, 3, 4, 5], monoclinic], - "Triclinic": [[0, 1, 2, 3, 4, 5], triclinic] - } - - lattyp, brav, sg_name, sg_nr = get_lattice_type(cryst) - # Decide which deformations should be used - axis, symm = deform[brav] + # Get the appropriate symmetry function and deformation axes + symm, axis, lattice_info = get_symmetry_function(cryst) systems = [] for a in axis: @@ -504,22 +659,8 @@ def get_elastic_tensor(cryst, systems): tuple(:math:`B_{ij}` float vector, residuals, solution rank, singular values)) ''' - # Deformation look-up table - # Perhaps the number of deformations for trigonal - # system could be reduced to [0,3] but better safe then sorry - deform = { - "Cubic": [[0, 3], regular], - "Hexagonal": [[0, 2, 3, 5], hexagonal], - "Trigonal": [[0, 1, 2, 3, 4, 5], trigonal], - "Tetragonal": [[0, 2, 3, 5], tetragonal], - "Orthorombic": [[0, 1, 2, 3, 4, 5], orthorombic], - "Monoclinic": [[0, 1, 2, 3, 4, 5], monoclinic], - "Triclinic": [[0, 1, 2, 3, 4, 5], triclinic] - } - - lattyp, brav, sg_name, sg_nr = get_lattice_type(cryst) - # Decide which deformations should be used - axis, symm = deform[brav] + # Get the appropriate symmetry function + symm, axis, lattice_info = get_symmetry_function(cryst) ul = [] sl = [] @@ -543,11 +684,15 @@ def get_elastic_tensor(cryst, systems): # TODO: Check the sign of the pressure array in the B <=> C relation if (symm == orthorombic): Cij = Bij[0] - array([-p, -p, -p, p, p, p, -p, -p, -p]) - elif (symm == tetragonal): + elif (symm == tetragonal_high): Cij = Bij[0] - array([-p, -p, p, p, -p, -p]) + elif (symm == tetragonal_low): + Cij = Bij[0] - array([-p, -p, p, p, -p, -p, p]) elif (symm == regular): Cij = Bij[0] - array([-p, p, -p]) - elif (symm == trigonal): + elif (symm == trigonal_high): + Cij = Bij[0] - array([-p, -p, p, p, -p, p]) + elif (symm == trigonal_low): Cij = Bij[0] - array([-p, -p, p, p, -p, p]) elif (symm == hexagonal): Cij = Bij[0] - array([-p, -p, p, p, -p]) From 906fd604f8f0e6dec3b1f94f07b4aca4e7595fdb Mon Sep 17 00:00:00 2001 From: "copilot-swe-agent[bot]" <198982749+Copilot@users.noreply.github.com> Date: Sun, 4 Jan 2026 20:26:35 +0000 Subject: [PATCH 3/4] Add test suite and export new symmetry functions Co-authored-by: jochym <5993422+jochym@users.noreply.github.com> --- elastic/__init__.py | 3 + tests/test_symmetry_variants.py | 214 ++++++++++++++++++++++++++++++++ 2 files changed, 217 insertions(+) create mode 100644 tests/test_symmetry_variants.py diff --git a/elastic/__init__.py b/elastic/__init__.py index 63ff5fc..144dade 100644 --- a/elastic/__init__.py +++ b/elastic/__init__.py @@ -66,6 +66,9 @@ from .elastic import get_BM_EOS, get_elastic_tensor from .elastic import get_elementary_deformations, scan_volumes from .elastic import get_pressure, get_strain, BMEOS +from .elastic import get_lattice_type, get_symmetry_function, get_cij_order +from .elastic import tetragonal, tetragonal_high, tetragonal_low +from .elastic import trigonal, trigonal_high, trigonal_low # To reach "elastic.__version__" attribute in other programs from ._version import __version__ diff --git a/tests/test_symmetry_variants.py b/tests/test_symmetry_variants.py new file mode 100644 index 0000000..2d4edff --- /dev/null +++ b/tests/test_symmetry_variants.py @@ -0,0 +1,214 @@ +#!/usr/bin/env python +# -*- coding: utf-8 -*- +""" +Test for trigonal and tetragonal symmetry variants + +This test verifies that the code correctly distinguishes between +high and low symmetry variants of trigonal and tetragonal crystal systems, +as described in the issue about independent elastic tensor components. + +References: +- https://github.com/jochym/Elastic/issues/[issue_number] +- Elasticity measurements on minerals: A review, EJM 21(3), 2009 +- https://github.com/libAtoms/matscipy/blob/master/matscipy/elasticity.py +""" + +from __future__ import print_function, division + +import sys +import os +import numpy as np +from ase.spacegroup import crystal + +# Add parent directory to path for imports +sys.path.insert(0, os.path.dirname(os.path.dirname(os.path.abspath(__file__)))) + +# Import elastic module functions +from elastic.elastic import ( + get_lattice_type, get_symmetry_function, get_cij_order, + tetragonal, tetragonal_high, tetragonal_low, + trigonal, trigonal_high, trigonal_low +) + + +def test_tetragonal_high_symmetry(): + """ + Test tetragonal high symmetry crystals (space groups 89-142). + These should have 6 independent elastic constants. + Example: TiO2 (rutile), space group 136 (P4_2/mnm) + """ + print("\n" + "="*60) + print("Testing Tetragonal High Symmetry") + print("="*60) + + # TiO2 rutile structure - space group 136 + a = 4.60 + c = 2.96 + cryst = crystal(['Ti', 'O'], [(0, 0, 0), (0.302, 0.302, 0)], + spacegroup=136, cellpar=[a, a, c, 90, 90, 90]) + + lattype, bravais, sg_name, sg_nr = get_lattice_type(cryst) + print(f"Structure: TiO2 (rutile)") + print(f"Space group: {sg_nr} ({sg_name})") + print(f"Bravais lattice: {bravais}") + + # Get symmetry function + symm_func, axes, _ = get_symmetry_function(cryst) + print(f"Symmetry function: {symm_func.__name__}") + + # Test matrix dimensions + u = np.array([0.01, 0.01, 0.01, 0.01, 0.01, 0.01]) + matrix = symm_func(u) + print(f"Matrix shape: {matrix.shape}") + + # Get elastic constant order + cij_order = get_cij_order(cryst) + print(f"Elastic constants: {cij_order}") + print(f"Number of constants: {len(cij_order)}") + + # Assertions + assert bravais == "Tetragonal", f"Expected Tetragonal, got {bravais}" + assert symm_func == tetragonal_high, f"Expected tetragonal_high for SG {sg_nr}" + assert matrix.shape == (6, 6), f"Expected (6,6), got {matrix.shape}" + assert len(cij_order) == 6, f"Expected 6 constants, got {len(cij_order)}" + expected_order = ('C_11', 'C_33', 'C_12', 'C_13', 'C_44', 'C_66') + assert cij_order == expected_order, f"Expected {expected_order}, got {cij_order}" + + print("✓ Test PASSED") + return True + + +def test_trigonal_high_symmetry(): + """ + Test trigonal high symmetry crystals (space groups 149-167). + These should have 6 independent elastic constants. + Example: Sb, space group 166 (R-3m) + """ + print("\n" + "="*60) + print("Testing Trigonal High Symmetry") + print("="*60) + + # Sb structure - space group 166 + a = 4.48 + c = 11.04 + cryst = crystal(['Sb'], [(0, 0, 0.24098)], + spacegroup=166, cellpar=[a, a, c, 90, 90, 120]) + + lattype, bravais, sg_name, sg_nr = get_lattice_type(cryst) + print(f"Structure: Sb") + print(f"Space group: {sg_nr} ({sg_name})") + print(f"Bravais lattice: {bravais}") + + # Get symmetry function + symm_func, axes, _ = get_symmetry_function(cryst) + print(f"Symmetry function: {symm_func.__name__}") + + # Test matrix dimensions + u = np.array([0.01, 0.01, 0.01, 0.01, 0.01, 0.01]) + matrix = symm_func(u) + print(f"Matrix shape: {matrix.shape}") + + # Get elastic constant order + cij_order = get_cij_order(cryst) + print(f"Elastic constants: {cij_order}") + print(f"Number of constants: {len(cij_order)}") + + # Assertions + assert bravais == "Trigonal", f"Expected Trigonal, got {bravais}" + assert symm_func == trigonal_high, f"Expected trigonal_high for SG {sg_nr}" + assert matrix.shape == (6, 6), f"Expected (6,6), got {matrix.shape}" + assert len(cij_order) == 6, f"Expected 6 constants, got {len(cij_order)}" + expected_order = ('C_11', 'C_33', 'C_12', 'C_13', 'C_44', 'C_14') + assert cij_order == expected_order, f"Expected {expected_order}, got {cij_order}" + + print("✓ Test PASSED") + return True + + +def test_symmetry_function_matrix_shapes(): + """ + Test that all symmetry functions return matrices with correct dimensions. + """ + print("\n" + "="*60) + print("Testing Symmetry Function Matrix Dimensions") + print("="*60) + + u = np.array([0.01, 0.01, 0.01, 0.01, 0.01, 0.01]) + + tests = [ + (tetragonal_high, (6, 6), "tetragonal_high"), + (tetragonal_low, (6, 7), "tetragonal_low"), + (trigonal_high, (6, 6), "trigonal_high"), + (trigonal_low, (6, 6), "trigonal_low"), + ] + + for func, expected_shape, name in tests: + matrix = func(u) + print(f"{name:20s}: {matrix.shape} (expected: {expected_shape})") + assert matrix.shape == expected_shape, \ + f"{name}: Expected {expected_shape}, got {matrix.shape}" + + print("✓ All matrix dimensions correct") + return True + + +def test_backward_compatibility(): + """ + Test that the old tetragonal() and trigonal() functions still work + and are aliases for the high symmetry variants. + """ + print("\n" + "="*60) + print("Testing Backward Compatibility") + print("="*60) + + u = np.array([0.01, 0.01, 0.01, 0.01, 0.01, 0.01]) + + # Test tetragonal() is an alias for tetragonal_high() + matrix_old = tetragonal(u) + matrix_new = tetragonal_high(u) + assert np.array_equal(matrix_old, matrix_new), \ + "tetragonal() should equal tetragonal_high()" + print("✓ tetragonal() == tetragonal_high()") + + # Test trigonal() is an alias for trigonal_high() + matrix_old = trigonal(u) + matrix_new = trigonal_high(u) + assert np.array_equal(matrix_old, matrix_new), \ + "trigonal() should equal trigonal_high()" + print("✓ trigonal() == trigonal_high()") + + print("✓ Backward compatibility maintained") + return True + + +def main(): + """Run all tests""" + print("\n" + "#"*60) + print("# Test Suite: Trigonal and Tetragonal Symmetry Variants") + print("#"*60) + + try: + test_symmetry_function_matrix_shapes() + test_tetragonal_high_symmetry() + test_trigonal_high_symmetry() + test_backward_compatibility() + + print("\n" + "#"*60) + print("# All tests PASSED ✓") + print("#"*60) + return 0 + + except AssertionError as e: + print(f"\n❌ Test FAILED: {e}") + import traceback + traceback.print_exc() + return 1 + except Exception as e: + print(f"\n❌ Unexpected error: {e}") + import traceback + traceback.print_exc() + return 1 + + +if __name__ == '__main__': + sys.exit(main()) From d32d7d2f2477aacb1b9d3d7271d408a502858fde Mon Sep 17 00:00:00 2001 From: "copilot-swe-agent[bot]" <198982749+Copilot@users.noreply.github.com> Date: Sun, 4 Jan 2026 20:28:19 +0000 Subject: [PATCH 4/4] Fix tetragonal_low matrix coefficient for consistency Co-authored-by: jochym <5993422+jochym@users.noreply.github.com> --- elastic/elastic.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/elastic/elastic.py b/elastic/elastic.py index 21a2473..dfe686e 100644 --- a/elastic/elastic.py +++ b/elastic/elastic.py @@ -177,9 +177,9 @@ def tetragonal_low(u): [[uxx, 0, uyy, uzz, 0, 0, 2*uxy], [uyy, 0, uxx, uzz, 0, 0, -2*uxy], [0, uzz, 0, uxx+uyy, 0, 0, 0 ], - [0, 0, 0, 0, 2*uxz, 0, 0 ], [0, 0, 0, 0, 2*uyz, 0, 0 ], - [uxy, 0, -uxy, 0, 0, 2*uxy, uxx-uyy]]) + [0, 0, 0, 0, 2*uxz, 0, 0 ], + [2*uxy, 0, -2*uxy, 0, 0, 2*uxy, uxx-uyy]]) def orthorombic(u):