diff --git a/src/scm/plams/interfaces/adfsuite/crs.py b/src/scm/plams/interfaces/adfsuite/crs.py index 5cd6351ad..a5134da3c 100644 --- a/src/scm/plams/interfaces/adfsuite/crs.py +++ b/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") diff --git a/src/scm/plams/interfaces/molecule/packmol.py b/src/scm/plams/interfaces/molecule/packmol.py index d76a4b1da..0247c478e 100644 --- a/src/scm/plams/interfaces/molecule/packmol.py +++ b/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, " diff --git a/unit_tests/test_packmol.py b/unit_tests/test_packmol.py index ddcebf0f4..61ba65a08 100644 --- a/unit_tests/test_packmol.py +++ b/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: