Add seed to packmol and rdkit (#12)

* run packmol with seed

* [pre-commit.ci] auto fixes from pre-commit.com hooks

for more information, see https://pre-commit.ci

* add and test rdkit seeds

* [pre-commit.ci] auto fixes from pre-commit.com hooks

for more information, see https://pre-commit.ci

* test and allow custom path

* [pre-commit.ci] auto fixes from pre-commit.com hooks

for more information, see https://pre-commit.ci

---------

Co-authored-by: pre-commit-ci[bot] <66853113+pre-commit-ci[bot]@users.noreply.github.com>

README.md

@@ -23,7 +23,9 @@ print(atoms)
 ```
 
 ## Packmol Interface
-If you have [packmol](https://github.com/m3g/packmol) (at least `v20.15.0`) you can use the rdkit2ase interface.
+
+If you have [packmol](https://github.com/m3g/packmol) (at least `v20.15.0`) you
+can use the rdkit2ase interface.
 
 ```py
 from rdkit2ase import pack, smiles2conformers

rdkit2ase/rdkit2ase.py

@@ -6,10 +6,10 @@
 import rdkit.Geometry
 
 
-def rdkit2ase(mol) -> ase.Atoms:
+def rdkit2ase(mol, seed: int = 42) -> ase.Atoms:
     """Convert an RDKit molecule to an ASE atoms object."""
     mol = rdkit.Chem.AddHs(mol)
-    rdkit.Chem.AllChem.EmbedMolecule(mol)
+    rdkit.Chem.AllChem.EmbedMolecule(mol, randomSeed=seed)
     rdkit.Chem.AllChem.UFFOptimizeMolecule(mol)
 
     return ase.Atoms(

rdkit2ase/utils/smiles.py

@@ -5,7 +5,7 @@
 from ..rdkit2ase import rdkit2ase
 
 
-def smiles2atoms(smiles: str) -> ase.Atoms:
+def smiles2atoms(smiles: str, seed: int = 42) -> ase.Atoms:
     """
     Convert a SMILES string to an ASE Atoms object.
 
@@ -16,7 +16,7 @@ def smiles2atoms(smiles: str) -> ase.Atoms:
         atoms (ase.Atoms): The Atoms object.
     """
     mol = Chem.MolFromSmiles(smiles)
-    return rdkit2ase(mol)
+    return rdkit2ase(mol, seed=seed)
 
 
 def smiles2conformers(

rdkit2ase/utils/solvate.py

@@ -39,6 +39,7 @@ def pack(
     seed: int = 42,
     tolerance: float = 2,
     logging: bool = False,
+    packmol: str = "packmol",
 ) -> ase.Atoms:
     """
     Pack the given molecules into a box with the specified density.
@@ -57,6 +58,8 @@ def pack(
         The tolerance for the packing algorithm, by default 2.
     logging : bool, optional
         If True, enables logging of the packing process, by default False.
+    packmol : str, optional
+        The path to the packmol executable, by default "packmol".
 
     Returns
     -------
@@ -93,6 +96,7 @@ def pack(
 tolerance {tolerance}
 filetype xyz
 output mixture.xyz
+seed {seed}
 pbc 0 0 0 {" ".join([f"{x:.6f}" for x in cell])}
     """
     for category, indices in enumerate(selected_idx):
@@ -114,7 +118,7 @@ def pack(
                 )
         (tmpdir / "pack.inp").write_text(file)
         subprocess.run(
-            "packmol < pack.inp",
+            f"{packmol} < pack.inp",
             cwd=tmpdir,
             shell=True,
             check=True,

tests/test_rdkit2ase.py

@@ -1,4 +1,6 @@
 import ase
+import numpy as np
+import numpy.testing as npt
 import pytest
 import rdkit
 
@@ -65,3 +67,25 @@ def test_smiles_to_atoms(smiles, formula):
     atoms = smiles2atoms(smiles)
     assert isinstance(atoms, ase.Atoms)
     assert atoms.get_chemical_formula() == formula
+
+
+def test_seeded_smiles_to_atoms():
+    atoms = smiles2atoms("CCO", seed=42)
+    npt.assert_almost_equal(
+        atoms.get_positions(),
+        [
+            [-0.92537086, 0.07420817, 0.03283984],
+            [0.51231816, -0.4191819, -0.07431482],
+            [1.3778212, 0.44937969, 0.60442827],
+            [-1.02252523, 1.07307188, -0.44285695],
+            [-1.60443695, -0.63678765, -0.48322232],
+            [-1.22359694, 0.14724363, 1.10020568],
+            [0.80577952, -0.5060158, -1.14510451],
+            [0.58521288, -1.42739618, 0.3853352],
+            [1.49479822, 1.24547818, 0.0226896],
+        ],
+    )
+    # now test with different seeds
+    assert not np.allclose(
+        atoms.get_positions(), smiles2atoms("CCO", seed=43).get_positions()
+    )

tests/test_solvate.py

@@ -17,6 +17,18 @@ def test_pack_pbc():
     assert len(atoms) == np.sum(atoms_dist < min_mol_dist * 0.99)
 
 
+def test_pack_seeded():
+    water = smiles2conformers("O", 1)
+
+    atoms1 = pack([water], [1], 1000, seed=42)
+    atoms2 = pack([water], [1], 1000, seed=42)
+
+    assert np.all(atoms1.get_positions() == atoms2.get_positions())
+
+    atoms3 = pack([water], [1], 1000, seed=43)
+    assert not np.all(atoms1.get_positions() == atoms3.get_positions())
+
+
 def test_pack_density():
     ethanol = smiles2conformers("CCO", 1)
     water = smiles2conformers("O", 2)