add spacegroup determination to geometry optimisation (#153)

* add spacegroup determination to geometry optimisation

---------

Co-authored-by: ElliottKasoar <45317199+ElliottKasoar@users.noreply.github.com>

janus_core/calculations/geom_opt.py

@@ -17,13 +17,15 @@
 
 from janus_core.helpers.janus_types import ASEOptArgs, ASEWriteArgs
 from janus_core.helpers.log import config_logger
-from janus_core.helpers.utils import none_to_dict
+from janus_core.helpers.utils import none_to_dict, spacegroup
 
 
 def optimize(  # pylint: disable=too-many-arguments,too-many-locals,too-many-branches
     struct: Atoms,
     fmax: float = 0.1,
     steps: int = 1000,
+    symmetry_tolerance: float = 0.001,
+    angle_tolerance: float = -1.0,
     filter_func: Optional[Callable] = DefaultFilter,
     filter_kwargs: Optional[dict[str, Any]] = None,
     optimizer: Callable = LBFGS,
@@ -45,6 +47,12 @@ def optimize(  # pylint: disable=too-many-arguments,too-many-locals,too-many-bra
         Default is 0.1.
     steps : int
         Set maximum number of optimization steps to run. Default is 1000.
+    symmetry_tolerance : float
+        Atom displacement tolerance for spglib symmetry determination, in Å.
+        Default is 0.001.
+    angle_tolerance : float
+        Angle precision for spglib symmetry determination, in degrees. Default is -1.0,
+        which means an internally optimized routine is used to judge symmetry.
     filter_func : Optional[callable]
         Apply constraints to atoms through ASE filter function.
         Default is `FrechetCellFilter` if available otherwise `ExpCellFilter`.
@@ -93,6 +101,13 @@ def optimize(  # pylint: disable=too-many-arguments,too-many-locals,too-many-bra
     log_kwargs.setdefault("name", __name__)
     logger = config_logger(**log_kwargs)
 
+    s_grp = spacegroup(struct, symmetry_tolerance, angle_tolerance)
+    message = f"Before optimisation spacegroup {s_grp}"
+    struct.info["initial_spacegroup"] = s_grp
+
+    if logger:
+        logger.info(message)
+
     if filter_func is not None:
         filtered_struct = filter_func(struct, **filter_kwargs)
         dyn = optimizer(filtered_struct, **opt_kwargs)
@@ -114,13 +129,18 @@ def optimize(  # pylint: disable=too-many-arguments,too-many-locals,too-many-bra
 
     converged = dyn.run(fmax=fmax, steps=steps)
 
+    s_grp = spacegroup(struct, symmetry_tolerance, angle_tolerance)
+    message = f"After optimisation spacegroup {s_grp}"
+    struct.info["final_spacegroup"] = s_grp
+
     # Calculate current maximum force
     if filter_func is not None:
         max_force = linalg.norm(filtered_struct.get_forces(), axis=1).max()
     else:
         max_force = linalg.norm(struct.get_forces(), axis=1).max()
 
     if logger:
+        logger.info(message)
         logger.info("Max force: %.6f", max_force)
 
     if not converged:

janus_core/helpers/utils.py

@@ -3,6 +3,42 @@
 from pathlib import Path
 from typing import Optional
 
+from ase import Atoms
+from spglib import get_spacegroup
+
+
+def spacegroup(
+    struct: Atoms, sym_tolerance: float = 0.001, angle_tolerance: float = -1.0
+) -> str:
+    """
+    Determine the spacegroup for a structure.
+
+    Parameters
+    ----------
+    struct : Atoms
+        Structure as an ase Atoms object.
+    sym_tolerance : float
+        Atom displacement tolerance for spglib symmetry determination, in Å.
+        Default is 0.001.
+    angle_tolerance : float
+        Angle precision for spglib symmetry determination, in degrees. Default is -1.0,
+        which means an internally optimized routine is used to judge symmetry.
+
+    Returns
+    -------
+    str
+        Spacegroup name.
+    """
+    return get_spacegroup(
+        cell=(
+            struct.get_cell(),
+            struct.get_scaled_positions(),
+            struct.get_atomic_numbers(),
+        ),
+        symprec=sym_tolerance,
+        angle_tolerance=angle_tolerance,
+    )
+
 
 def none_to_dict(dictionaries: list[Optional[dict]]) -> list[dict]:
     """

pyproject.toml

@@ -34,6 +34,7 @@ typer = "^0.9.0"
 typer-config = "^1.4.0"
 phonopy = "^2.23.1"
 seekpath = "^2.1.0"
+spglib = "^2.3.0"
 
 [tool.poetry.group.dev.dependencies]
 coverage = {extras = ["toml"], version = "^7.4.1"}

tests/data/NaCl-sg.cif

@@ -0,0 +1,33 @@
+data_image0
+_chemical_formula_structural       NaClNaClNaClNaCl
+_chemical_formula_sum              "Na4 Cl4"
+_cell_length_a       5.65
+_cell_length_b       5.64
+_cell_length_c       5.64
+_cell_angle_alpha    90.0
+_cell_angle_beta     90.0
+_cell_angle_gamma    90.0
+
+_space_group_name_H-M_alt    "P 1"
+_space_group_IT_number       1
+
+loop_
+  _space_group_symop_operation_xyz
+  'x, y, z'
+
+loop_
+  _atom_site_type_symbol
+  _atom_site_label
+  _atom_site_symmetry_multiplicity
+  _atom_site_fract_x
+  _atom_site_fract_y
+  _atom_site_fract_z
+  _atom_site_occupancy
+  Na  Na1       1.0  0.0  0.0  0.0  1.0000
+  Cl  Cl1       1.0  0.5  0.0  0.0  1.0000
+  Na  Na2       1.0  0.0  0.5  0.5  1.0000
+  Cl  Cl2       1.0  0.5  0.5  0.5  1.0000
+  Na  Na3       1.0  0.5  0.0  0.5  1.0000
+  Cl  Cl3       1.0  0.0  0.0  0.5  1.0000
+  Na  Na4       1.0  0.5  0.5  0.0  1.0000
+  Cl  Cl4       1.0  0.0  0.5  0.0  1.0000

tests/test_geom_opt.py

@@ -187,3 +187,21 @@ def test_restart(tmp_path):
     )
     final_energy = single_point.run("energy")["energy"]
     assert final_energy < intermediate_energy
+
+
+def test_space_group():
+    """Test spacegroup of the structure."""
+
+    single_point = SinglePoint(
+        struct_path=DATA_PATH / "NaCl-sg.cif",
+        architecture="mace_mp",
+        calc_kwargs={"model": MODEL_PATH},
+    )
+
+    optimize(
+        single_point.struct,
+        fmax=0.001,
+    )
+
+    assert single_point.struct.info["initial_spacegroup"] == "I4/mmm (139)"
+    assert single_point.struct.info["final_spacegroup"] == "Fm-3m (225)"