From c231cbd3d5147ee920a37b6ee9dd236b376bcf5a Mon Sep 17 00:00:00 2001
From: Matthew Horton <mkhorton@users.noreply.github.com>
Date: Sun, 11 Feb 2024 18:03:53 -0800
Subject: [PATCH] Update `JonesFaithfulTransformation` to use sympy

Signed-off-by: Matthew Horton <mkhorton@users.noreply.github.com>
---
 pymatgen/symmetry/settings.py | 26 ++++++++++++++++++--------
 1 file changed, 18 insertions(+), 8 deletions(-)

diff --git a/pymatgen/symmetry/settings.py b/pymatgen/symmetry/settings.py
index a0d3050984f..87d4fde6f35 100644
--- a/pymatgen/symmetry/settings.py
+++ b/pymatgen/symmetry/settings.py
@@ -6,6 +6,8 @@
 from fractions import Fraction
 
 import numpy as np
+from sympy import Matrix
+from sympy.parsing.sympy_parser import parse_expr
 
 from pymatgen.core.lattice import Lattice
 from pymatgen.core.operations import MagSymmOp, SymmOp
@@ -99,21 +101,29 @@ def parse_transformation_string(
             b_change, o_shift = transformation_string.split(";")
             basis_change = b_change.split(",")
             origin_shift = o_shift.split(",")
+
             # add implicit multiplication symbols
             basis_change = [
                 re.sub(r"(?<=\w|\))(?=\() | (?<=\))(?=\w) | (?<=(\d|a|b|c))(?=([abc]))", r"*", string, flags=re.X)
                 for string in basis_change
             ]
-            # should be fine to use eval here but be mindful for security
-            # reasons
-            # see http://lybniz2.sourceforge.net/safeeval.html
-            # could replace with regex? or sympy expression?
-            P = np.array([eval(x, {"__builtins__": None}, {"a": a, "b": b, "c": c}) for x in basis_change])
-            P = P.transpose()  # by convention
+
+            # basic input sanitation
+            allowed_chars = "0123456789+-*/.abc()"
+            basis_change = ["".join([c for c in string if c in allowed_chars]) for string in basis_change]
+
+            # requires round-trip to sympy to evaluate
+            # (alternatively, `numexpr` looks like a nice solution but requires an additional dependency)
+            basis_change = [
+                parse_expr(string).subs({"a": Matrix(a), "b": Matrix(b), "c": Matrix(c)}) for string in basis_change
+            ]
+            # convert back to numpy, perform transpose by convention
+            P = np.array(basis_change, dtype=float).T[0]
+
             p = [float(Fraction(x)) for x in origin_shift]
             return P, p
-        except Exception:
-            raise ValueError("Failed to parse transformation string.")
+        except Exception as exc:
+            raise ValueError(f"Failed to parse transformation string: {exc}")
 
     @property
     def P(self) -> list[list[float]]: