accept scm.base.ChemicalSystem as input for packmol function SO--

src/scm/plams/interfaces/adfsuite/crs.py

@@ -353,7 +353,7 @@ def get_x_axis(array: np.ndarray, x_axis: Optional[Union[str, np.ndarray]]) -> n
             import matplotlib
 
             if plot_fig:
-                if terminal == "jupyter":
+                if terminal == "jupyter" and ipython is not None:
                     ipython.run_line_magic("matplotlib", "inline")
                 else:
                     matplotlib.use("TkAgg")

src/scm/plams/interfaces/molecule/packmol.py

@@ -23,6 +23,8 @@
         pass
 
 T = TypeVar("T")
+MoleculeLike = Union[Molecule, "ChemicalSystem"]
+MoleculesInput = Union[List[MoleculeLike], MoleculeLike]
 
 __all__ = [
     "packmol",
@@ -38,6 +40,25 @@ def tolist(x: Union[T, List[T]]) -> List[T]:
     return x if isinstance(x, list) else [x]
 
 
+def _to_plams_molecule(molecule: MoleculeLike) -> Molecule:
+    if isinstance(molecule, Molecule):
+        return molecule
+
+    try:
+        from scm.base import ChemicalSystem
+    except ImportError:
+        pass
+    else:
+        if isinstance(molecule, ChemicalSystem):
+            parsed_molecules = cast(Dict[str, Molecule], AMSJob.from_input(str(molecule)).molecule)
+            return parsed_molecules[""]
+
+    raise TypeError(
+        "molecules must contain only PLAMS Molecule or ChemicalSystem objects, "
+        f"but received {type(molecule).__name__}"
+    )
+
+
 class PackMolError(MoleculeError):
     pass
 
@@ -482,7 +503,7 @@ def _guess_molecular_volume(molecule: Union[Molecule, "ChemicalSystem"]) -> floa
 
 @overload
 def packmol(
-    molecules: Union[List[Molecule], Molecule],
+    molecules: MoleculesInput,
     mole_fractions: Optional[List[float]] = ...,
     density: Optional[float] = ...,
     n_atoms: Optional[int] = ...,
@@ -501,7 +522,7 @@ def packmol(
 ) -> Molecule: ...
 @overload
 def packmol(
-    molecules: Union[List[Molecule], Molecule],
+    molecules: MoleculesInput,
     mole_fractions: Optional[List[float]] = ...,
     density: Optional[float] = ...,
     n_atoms: Optional[int] = ...,
@@ -520,7 +541,7 @@ def packmol(
 ) -> Tuple[Molecule, Dict[str, Any]]: ...
 @overload
 def packmol(
-    molecules: Union[List[Molecule], Molecule],
+    molecules: MoleculesInput,
     mole_fractions: Optional[List[float]] = ...,
     density: Optional[float] = ...,
     n_atoms: Optional[int] = ...,
@@ -539,7 +560,7 @@ def packmol(
 ) -> Tuple[None, Dict[str, Any]]: ...
 @overload
 def packmol(
-    molecules: Union[List[Molecule], Molecule],
+    molecules: MoleculesInput,
     mole_fractions: Optional[List[float]] = ...,
     density: Optional[float] = ...,
     n_atoms: Optional[int] = ...,
@@ -558,7 +579,7 @@ def packmol(
 ) -> Molecule: ...
 @overload
 def packmol(
-    molecules: Union[List[Molecule], Molecule],
+    molecules: MoleculesInput,
     mole_fractions: Optional[List[float]] = ...,
     density: Optional[float] = ...,
     n_atoms: Optional[int] = ...,
@@ -576,7 +597,7 @@ def packmol(
     _return_only_details: bool = False,
 ) -> Tuple[Molecule, Dict[str, Any]]: ...
 def packmol(
-    molecules: Union[List[Molecule], Molecule],
+    molecules: MoleculesInput,
     mole_fractions: Optional[List[float]] = None,
     density: Optional[float] = None,
     n_atoms: Optional[int] = None,
@@ -600,8 +621,8 @@ def packmol(
     It is *strongly recommended* to specify ``density`` and/or ``box_bounds``. Otherwise you will
     get a (very inaccurate) guessed density in a cubic box (experimental feature).
 
-    molecules : |Molecule| or list of Molecule
-        The molecules to pack
+    molecules : |Molecule|, ChemicalSystem, or list
+        The molecules to pack. A list may contain any mixture of PLAMS ``Molecule`` and ``ChemicalSystem`` objects.
 
     mole_fractions : list of float
         The mole fractions (in the same order as ``molecules``). Cannot be combined with ``n_molecules``. If not given, an equal (molar) mixture of all components will be created.
@@ -738,6 +759,11 @@ def packmol(
         if region_names is not None and isinstance(region_names, list):
             raise ValueError("Illegal combination of arguments: region_names is a list, when molecules is not")
 
+    if isinstance(molecules, list):
+        molecules = [_to_plams_molecule(molecule) for molecule in molecules]
+    else:
+        molecules = _to_plams_molecule(molecules)
+
     if n_atoms is not None and n_molecules is not None and not one_n_molecules_missing:
         raise ValueError(
             "Illegal combination of arguments: n_atoms and n_molecules are mutually exclusive, "

unit_tests/test_packmol.py

@@ -863,6 +863,54 @@ def test_pack_mol_happy(self, test_case):
             for i, (m, _) in enumerate(mols.values()):
                 assert all(a.properties.region == {f"mol{i}"} for a in m.atoms)
 
+    @pytest.mark.parametrize("use_mixed_inputs", [False, True], ids=["chemical-systems", "mixed"])
+    def test_pack_mol_accepts_chemical_systems(self, use_mixed_inputs: bool) -> None:
+        skip_if_no_scm_base()
+
+        from scm.utils.conversions import plams_molecule_to_chemsys
+
+        water_ucs = plams_molecule_to_chemsys(self.water)
+        acetonitrile_ucs = plams_molecule_to_chemsys(self.acetonitrile)
+        molecules = [water_ucs, self.acetonitrile if use_mixed_inputs else acetonitrile_ucs]
+
+        _, details = packmol(
+            molecules=molecules,
+            n_molecules=[2, 3],
+            density=1.0,
+            executable=".",
+            _return_only_details=True,
+        )
+
+        assert details["n_atoms"] == 2 * len(self.water) + 3 * len(self.acetonitrile)
+        assert details["n_molecules"] == [2, 3]
+        assert details["mole_fractions"] == [0.4, 0.6]
+
+    def test_pack_mol_accepts_single_chemical_system(self) -> None:
+        skip_if_no_scm_base()
+
+        from scm.utils.conversions import plams_molecule_to_chemsys
+
+        water_ucs = plams_molecule_to_chemsys(self.water)
+        _, details = packmol(
+            molecules=water_ucs,
+            n_molecules=2,
+            density=1.0,
+            executable=".",
+            _return_only_details=True,
+        )
+
+        assert details["n_atoms"] == 2 * len(self.water)
+        assert details["n_molecules"] == [2]
+
+    def test_pack_mol_rejects_unsupported_molecule_type(self) -> None:
+        with pytest.raises(TypeError, match="only PLAMS Molecule or ChemicalSystem objects.*str"):
+            packmol(
+                molecules=[self.water, "not a molecule"],
+                n_molecules=[1, 1],
+                density=1.0,
+                executable=".",
+            )
+
 
 class TestGuessDensity: