Skip to content

Add support for high/low symmetry variants in tetragonal and trigonal systems #82

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Draft
Copilot wants to merge 4 commits into
base: master
Choose a base branch
from
Draft

Add support for high/low symmetry variants in tetragonal and trigonal systems #82

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
4 commits
Select commit Hold shift + click to select a range
df757ca
Initial plan
Copilot Jan 4, 2026
4835ac2
Add high/low symmetry variants for trigonal and tetragonal systems
Copilot Jan 4, 2026
906fd60
Add test suite and export new symmetry functions
Copilot Jan 4, 2026
d32d7d2
Fix tetragonal_low matrix coefficient for consistency
Copilot Jan 4, 2026
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
3 changes: 3 additions & 0 deletions elastic/__init__.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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__
221 changes: 183 additions & 38 deletions elastic/elastic.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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::
Expand All @@ -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*uyz, 0, 0 ],
[0, 0, 0, 0, 2*uxz, 0, 0 ],
[2*uxy, 0, -2*uxy, 0, 0, 2*uxy, uxx-uyy]])


def orthorombic(u):
'''
Equation matrix generation for the orthorombic lattice.
Expand Down Expand Up @@ -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
Expand All @@ -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 ],
Expand All @@ -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
Expand Down Expand Up @@ -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):
Expand Down Expand Up @@ -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.

Expand Down Expand Up @@ -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:
Expand Down Expand Up @@ -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 = []
Expand All @@ -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])
Expand Down
Loading