Merge pull request #13 from kousuke-nakano/devel

Updated SHRY. Fixed the bug reported in the discussion #12

.github/workflows/shry-pytest.yml

@@ -6,6 +6,8 @@ name: SHRY pytest
 on:
   push:
     branches: [ "master", "devel" ]
+  pull_request:
+    branches: [ "main", "devel" ]
 
 jobs:
   build:
@@ -25,7 +27,7 @@ jobs:
       run: |
         python -m pip install --upgrade pip
         python -m pip install flake8 pytest
-        python -m pip install .    
+        python -m pip install .
     - name: Lint with flake8
       run: |
         # stop the build if there are Python syntax errors or undefined names
@@ -34,6 +36,6 @@ jobs:
         flake8 . --count --exit-zero --max-complexity=10 --max-line-length=127 --statistics
     - name: Test with pytest
       run: |
-        pytest
+        pytest -v
 
 

shry/core.py

@@ -31,6 +31,7 @@
 from scipy.special import comb
 from sympy.utilities.iterables import multiset_permutations
 from tabulate import tabulate
+from pymatgen.core.periodic_table import get_el_sp
 
 # shry modules
 from . import const
@@ -47,7 +48,9 @@
 
 
 def get_integer_formula_and_factor(
-    self, max_denominator: int = 10000, iupac_ordering: bool = False
+    self,
+    max_denominator: int = int(1 / const.DEFAULT_ATOL),
+    iupac_ordering: bool = False,
 ) -> Tuple[str, float]:
     """
     The default Composition groups together different ox states which is not ideal...
@@ -69,14 +72,19 @@ def to_int_dict(self):
     Returns:
         Dict with element symbol and integer amount
     """
-    _, factor = self.get_integer_formula_and_factor()
+    _, factor = self.get_integer_formula_and_factor(
+        max_denominator=int(1 / const.DEFAULT_ATOL)
+    )
     int_dict = {e: int(a) for e, a in (self / factor).as_dict().items()}
+
     # be safe: Composition groups together different ox states which is not ideal...
-    if not all(
-        np.isclose(x * factor, y)
-        for x, y in zip(int_dict.values(), self.as_dict().values())
-    ):
-        raise ValueError("Composition is not rational!")
+    for x, y in zip(int_dict.values(), self.as_dict().values()):
+        if not np.isclose(x * factor, y, atol=const.DEFAULT_ATOL):
+            raise ValueError(
+                "Composition (Occupancy) is not rational! Please try to increase significant digits "
+                "e.g., 1/3 = 0.3333 -> 1/3 = 0.3333333333333."
+            )
+
     return int_dict
 
 
@@ -85,23 +93,65 @@ def inted_composition(self):
     """
     Return Composition instance with integer formula
     """
-    _, factor = self.get_integer_formula_and_factor()
+    _, factor = self.get_integer_formula_and_factor(
+        max_denominator=int(1 / const.DEFAULT_ATOL)
+    )
     int_comp = self / factor
 
     # be safe
     int_dict = {e: int(a) for e, a in int_comp.as_dict().items()}
     if not all(
-        np.isclose(x * factor, y)
+        np.isclose(x * factor, y, atol=const.DEFAULT_ATOL)
         for x, y in zip(int_dict.values(), self.as_dict().values())
     ):
-        raise ValueError("Composition is not rational!")
+        raise ValueError(
+            "Composition (Occupancy) is not rational! Please try to increase significant digits "
+            "e.g., 1/3 = 0.3333 -> 1/3 = 0.3333333333333."
+        )
 
     return int_comp
 
 
+def formula_double_format_tol(
+    afloat, ignore_ones=True, tol: float = const.DEFAULT_ATOL * 10
+):
+    """
+    This function is used to make pretty formulas by formatting the amounts.
+    Instead of Li1.0 Fe1.0 P1.0 O4.0, you get LiFePO4.
+
+    Args:
+        afloat (float): a float
+        ignore_ones (bool): if true, floats of 1 are ignored.
+        tol (float): Tolerance to round to nearest int. i.e. 2.0000000001 -> 2
+
+    Returns:
+        A string representation of the float for formulas.
+    """
+    if ignore_ones and afloat == 1:
+        return ""
+    if abs(afloat - round(afloat)) < tol:
+        return round(afloat)
+    return round(afloat, 8)
+
+
+@property
+def formula(self) -> str:
+    """
+    Returns a formula string, with elements sorted by electronegativity,
+    e.g., Li4 Fe4 P4 O16.
+    """
+    sym_amt = self.get_el_amt_dict()
+    syms = sorted(sym_amt, key=lambda sym: get_el_sp(sym).X)
+    formula = [
+        f"{s}{formula_double_format_tol(sym_amt[s], False)}" for s in syms
+    ]
+    return " ".join(formula)
+
+
 Composition.to_int_dict = to_int_dict
 Composition.get_integer_formula_and_factor = get_integer_formula_and_factor
 Composition.inted_composition = inted_composition
+Composition.formula = formula
 
 
 class PatchedSymmetrizedStructure(SymmetrizedStructure):
@@ -519,9 +569,20 @@ def structure(self, structure):
                 disorder_sites.append(site)
                 # Ad hoc fix: if occupancy is less than 1, stop.
                 # TODO: Automatic vacancy handling
-                if not np.isclose(site.species.num_atoms, 1):
+                if not np.isclose(
+                    site.species.num_atoms, 1.0, atol=self._atol
+                ):
+                    logging.warning(
+                        f"The occupancy of the site {site.species} is {site.species.num_atoms}."
+                    )
+                    logging.warning(
+                        f"This should be 1 within the torelance, atol={self._atol}."
+                    )
+                    logging.warning(
+                        "If you want to consider vacancy sites, please add pseudo atoms."
+                    )
                     raise RuntimeError(
-                        "Please fill vacancy sites with pseudo atoms"
+                        "The sum of number of occupancies is not 1."
                     )
         if not disorder_sites:
             logging.warning("No disorder sites found within the Structure.")
@@ -548,7 +609,10 @@ def structure(self, structure):
                     (
                         f"Can't fit {integer_formula} "
                         f"within {len(sites)} sites "
-                        f"(enlarge by {formula_unit_sum/len(sites):.4f}x)."
+                        f"(enlarge by {formula_unit_sum/len(sites):.4f}x). "
+                        f"If the integer composition ({integer_formula}) is not what you expected, "
+                        "please try to increase the precision of the occupancy "
+                        "e.g., 1/3 = 0.3333 -> 1/3 = 0.3333333333333."
                     )
                 )
             self.disorder_groups[orbit] = sites

shry/script.py

@@ -196,7 +196,7 @@ def main():  # pylint: disable=missing-function-docstring
     group.add_argument(
         "--atol",
         type=float,
-        default=const.DEFAULT_SYMPREC,
+        default=const.DEFAULT_ATOL,
         help="Discretization absolute tolerance (angstrom).",
     )
     group.add_argument(

tests/test_core.py

@@ -23,11 +23,13 @@
     TooBigError,
 )
 from shry.main import LabeledStructure, ScriptHelper
+from shry import const
 from helper import chdir
 
 # Tolerances
-SHRY_TOLERANCE = 0.01  # angstrom
-SHRY_ANGLE_TOLERANCE = 5.0  # degree
+SHRY_TOLERANCE = const.DEFAULT_SYMPREC
+SHRY_ANGLE_TOLERANCE = const.DEFAULT_ANGLE_TOLERANCE
+SHRY_ATOL = const.DEFAULT_ATOL
 
 # PatternMaker basic functions.
 
@@ -510,6 +512,7 @@ def test_benchmark():
             structure,
             symprec=SHRY_TOLERANCE,
             angle_tolerance=SHRY_ANGLE_TOLERANCE,
+            atol=SHRY_ATOL
         )
         count_obtained = s.count()
         count_ref = equivalent