packmol_around docs, clean up workdir SO103

CHANGELOG.md

@@ -16,6 +16,7 @@ This changelog is effective from the 2025 releases.
 * `AMSJob` can accept an AMS `ChemicalSystem` instead of a PLAMS `Molecule` as an input system
 * Specific `ConfigSettings` and related settings classes with explicitly defined fields
 * Support for work functions: `AMSResults.get_work_function_results` and `plot_work_function`
+* New `packmol_around` function for packing in non-orthorhombic boxes.
 * Example on `MoleculeFormats`
 * Script `generate_example.sh` to generate documentation pages from notebook examples
 * GitHub workflows for CI and publishing to PyPI

doc/source/components/mol_packmol.rst

@@ -13,6 +13,7 @@ Packmol (`Packmol website <https://m3g.github.io/packmol/download.shtml>`__) is
 The following functions eixst:
 
 * ``packmol`` (for fluids with 1 or more components)
+* ``packmol_around`` (for fluids with 1 or more components, used to pack around an existing system in AMS2025+)
 * ``packmol_on_slab`` (for solid/liquid or solid/gas interfaces with 1 or more components in the fluid)
 * ``packmol_in_void`` (for packing molecules inside crystal voids)
 * ``packmol_microsolvation`` (for microsolvation of a solute with a solvent)
@@ -27,6 +28,8 @@ the packmol program included with the Amsterdam Modeling Suite will be used.
 
 .. autofunction:: packmol
 
+.. autofunction:: packmol_around
+
 .. autofunction:: packmol_on_slab
 
 .. autofunction:: packmol_in_void

doc/source/examples/PackMolExample/PackMol.ipynb

@@ -255,11 +255,7 @@
     "print(\"water-5.xyz: pure liquid in non-orthorhombic box (requires AMS2025 or later)\")\n",
     "# Non-orthorhombic boxes use UFF MD simulations behind the scenes\n",
     "# You can pack inside any lattice using the packmol_around function\n",
-    "from scm.plams import init, Settings\n",
     "\n",
-    "s = Settings()\n",
-    "s.log.stdout = 0\n",
-    "init(config_settings=s)\n",
     "box = Molecule()\n",
     "box.lattice = [[10.0, 2.0, -1.0], [-5.0, 8.0, 0.0], [0.0, -2.0, 11.0]]\n",
     "out = packmol_around(box, molecules=[water], n_molecules=[32])\n",

doc/source/examples/PackMolExample/PackMol.rst.include

@@ -123,11 +123,7 @@ First, create the gasphase molecule:
     print("water-5.xyz: pure liquid in non-orthorhombic box (requires AMS2025 or later)")
     # Non-orthorhombic boxes use UFF MD simulations behind the scenes
     # You can pack inside any lattice using the packmol_around function
-    from scm.plams import init, Settings
 
-    s = Settings()
-    s.log.stdout = 0
-    init(config_settings=s)
     box = Molecule()
     box.lattice = [[10.0, 2.0, -1.0], [-5.0, 8.0, 0.0], [0.0, -2.0, 11.0]]
     out = packmol_around(box, molecules=[water], n_molecules=[32])

doc/source/general.rst

@@ -112,6 +112,7 @@ Added
 
 * Support for AMS ``ChemicalSystem`` within |AMSJob| and |AMSResults|. |AMSJob| can accept a ``ChemicalSystem`` as an input system, and the methods :meth:`~scm.plams.interfaces.adfsuite.ams.AMSResults.get_system`, :meth:`~scm.plams.interfaces.adfsuite.ams.AMSResults.get_input_system` and :meth:`~scm.plams.interfaces.adfsuite.ams.AMSResults.get_main_system` on |AMSResults| return a ``ChemicalSystem``. These provide the option to use a ``ChemicalSystem`` in place of a PLAMS ``Molecule``.
 * Support for work functions through :meth:`~scm.plams.interfaces.adfsuite.ams.AMSResults.get_work_function_results` and :func:`~scm.plams.tools.plot.plot_work_function`
+* Added :meth:`~scm.plams.interfaces.molecule.packmol.packmol_around` method to pack around an existing molecule or pack molecules into a non-orthorhombic box.
 
 Changed
 ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

examples/PackMol.py

@@ -12,7 +12,6 @@
 
 # ## Helper functions
 
-
 def printsummary(mol, details=None):
     if details:
         density = details["density"]
@@ -31,57 +30,53 @@ def printsummary(mol, details=None):
 # First, create the gasphase molecule:
 
 water = from_smiles("O")
-plot_molecule(water)
+plot_molecule(water);
 
 
 print("pure liquid from approximate number of atoms and exact density (in g/cm^3), cubic box with auto-determined size")
 out = packmol(water, n_atoms=194, density=1.0)
 printsummary(out)
 out.write("water-1.xyz")
-plot_molecule(out)
+plot_molecule(out);
 
 
 print("pure liquid from approximate density (in g/cm^3) and an orthorhombic box")
 out = packmol(water, density=1.0, box_bounds=[0.0, 0.0, 0.0, 8.0, 12.0, 14.0])
 printsummary(out)
 out.write("water-2.xyz")
-plot_molecule(out)
+plot_molecule(out);
 
 
 print("pure liquid with explicit number of molecules and exact density")
 out = packmol(water, n_molecules=64, density=1.0)
 printsummary(out)
 out.write("water-3.xyz")
-plot_molecule(out)
+plot_molecule(out);
 
 
 print("pure liquid with explicit number of molecules and box")
 out = packmol(water, n_molecules=64, box_bounds=[0.0, 0.0, 0.0, 12.0, 13.0, 14.0])
 printsummary(out)
 out.write("water-4.xyz")
-plot_molecule(out)
+plot_molecule(out);
 
 
 print("water-5.xyz: pure liquid in non-orthorhombic box (requires AMS2025 or later)")
 # Non-orthorhombic boxes use UFF MD simulations behind the scenes
 # You can pack inside any lattice using the packmol_around function
-from scm.plams import init, Settings
 
-s = Settings()
-s.log.stdout = 0
-init(config_settings=s)
 box = Molecule()
 box.lattice = [[10.0, 2.0, -1.0], [-5.0, 8.0, 0.0], [0.0, -2.0, 11.0]]
 out = packmol_around(box, molecules=[water], n_molecules=[32])
 out.write("water-5.xyz")
-plot_molecule(out)
+plot_molecule(out);
 
 
 print("Experimental feature (AMS2025): guess density for pure liquid")
 print("Note: This density is meant to be equilibrated with NPT MD. It can be very inaccurate!")
 out = packmol(water, n_atoms=100)
 print(f"Guessed density: {out.get_density():.2f} kg/m^3")
-plot_molecule(out)
+plot_molecule(out);
 
 
 # ## Water-acetonitrile mixture (fluid with 2 or more components)
@@ -120,7 +115,7 @@ def printsummary(mol, details=None):
 )
 printsummary(out, details)
 out.write("water-acetonitrile-1.xyz")
-plot_molecule(out)
+plot_molecule(out);
 
 
 # The ``details`` is a dictionary as follows:
@@ -139,7 +134,7 @@ def printsummary(mol, details=None):
 )
 printsummary(out, details)
 out.write("water-acetonitrile-2.xyz")
-plot_molecule(out)
+plot_molecule(out);
 
 
 print("2-1 water-acetonitrile from explicit number of molecules and density, cubic box with auto-determined size")
@@ -151,7 +146,7 @@ def printsummary(mol, details=None):
 )
 printsummary(out, details)
 out.write("water-acetonitrile-3.xyz")
-plot_molecule(out)
+plot_molecule(out);
 
 
 print("2-1 water-acetonitrile from explicit number of molecules and box")
@@ -162,18 +157,18 @@ def printsummary(mol, details=None):
 )
 printsummary(out)
 out.write("water-acetonitrile-4.xyz")
-plot_molecule(out)
+plot_molecule(out);
 
 
 print("Experimental feature (AMS2025): guess density for mixture")
 print("Note: This density is meant to be equilibrated with NPT MD. It can be very inaccurate!")
 out = packmol([water, acetonitrile], mole_fractions=[x_water, x_acetonitrile], n_atoms=100)
 print(f"Guessed density: {out.get_density():.2f} kg/m^3")
-plot_molecule(out)
+plot_molecule(out);
 
 
 # ## Pack inside sphere
-#
+# 
 # Set ``sphere=True`` to pack in a sphere (non-periodic) instead of in a periodic box. The sphere will be centered near the origin.
 
 print("water in a sphere from exact density and number of molecules")
@@ -182,7 +177,7 @@ def printsummary(mol, details=None):
 print(f"Radius  of sphere: {details['radius']:.3f} ang.")
 print(f"Center of mass xyz (ang): {out.get_center_of_mass()}")
 out.write("water-sphere.xyz")
-plot_molecule(out)
+plot_molecule(out);
 
 
 print(
@@ -199,13 +194,13 @@ def printsummary(mol, details=None):
 )
 printsummary(out, details)
 out.write("water-acetonitrile-sphere.xyz")
-plot_molecule(out)
+plot_molecule(out);
 
 
 # ## Packing ions, total system charge
-#
+# 
 # The total system charge will be sum of the charges of the constituent molecules.
-#
+# 
 # In PLAMS, ``molecule.properties.charge`` specifies the charge:
 
 ammonium = from_smiles("[NH4+]")  # ammonia.properties.charge == +1
@@ -219,7 +214,7 @@ def printsummary(mol, details=None):
 tot_charge = out.properties.get("charge", 0)
 print(f"Total charge of packmol-generated system: {tot_charge}")
 out.write("water-ammonium-chloride.xyz")
-plot_molecule(out)
+plot_molecule(out);
 
 
 # ## Microsolvation
@@ -233,7 +228,7 @@ def printsummary(mol, details=None):
 out.write("acetonitrile-microsolvated.xyz")
 
 figsize = (3, 3)
-plot_molecule(out, figsize=figsize)
+plot_molecule(out, figsize=figsize);
 
 
 # ## Solid-liquid or solid-gas interfaces
@@ -244,21 +239,21 @@ def printsummary(mol, details=None):
 
 rotation = "90x,0y,0z"  # sideview of slab
 slab = fromASE(fcc111("Al", size=(4, 6, 3), vacuum=15.0, orthogonal=True, periodic=True))
-plot_molecule(slab, figsize=figsize, rotation=rotation)
+plot_molecule(slab, figsize=figsize, rotation=rotation);
 
 
 print("water surrounding an Al slab, from an approximate density")
 out = packmol_around(slab, water, density=1.0)
 printsummary(out)
 out.write("al-water-pure.xyz")
-plot_molecule(out, figsize=figsize, rotation=rotation)
+plot_molecule(out, figsize=figsize, rotation=rotation);
 
 
 print("2-1 water-acetonitrile mixture surrounding an Al slab, from mole fractions and an approximate density")
 out = packmol_around(slab, [water, acetonitrile], mole_fractions=[x_water, x_acetonitrile], density=density)
 printsummary(out)
 out.write("al-water-acetonitrile.xyz")
-plot_molecule(out, figsize=figsize, rotation=rotation)
+plot_molecule(out, figsize=figsize, rotation=rotation);
 
 
 from ase.build import surface
@@ -274,15 +269,15 @@ def printsummary(mol, details=None):
 
 
 # ## Pack inside voids in crystals
-#
+# 
 # Use the ``packmol_around`` function. You can decrease ``tolerance`` if you need to pack very tightly. The default value for ``tolerance`` is 2.0.
 
 from scm.plams import fromASE
 from ase.build import bulk
 
 bulk_Al = fromASE(bulk("Al", cubic=True).repeat((3, 3, 3)))
 rotation = "-85x,5y,0z"
-plot_molecule(bulk_Al, rotation=rotation, radii=0.4)
+plot_molecule(bulk_Al, rotation=rotation, radii=0.4);
 
 
 out = packmol_around(
@@ -297,9 +292,9 @@ def printsummary(mol, details=None):
 
 
 # ## Bonds, atom properties (force field types, regions, ...)
-#
+# 
 # The ``packmol()`` function accepts the arguments ``keep_bonds`` and ``keep_atom_properties``. These options will keep the bonds defined for the constitutent molecules, as well as any atomic properties.
-#
+# 
 # The bonds and atom properties are easiest to see by printing the System block for an AMS job:
 
 water = from_smiles("O")
@@ -322,7 +317,7 @@ def printsummary(mol, details=None):
 print(AMSJob(molecule=out).get_input())
 
 
-# By default, the ``packmol()`` function assigns regions called ``mol0``, ``mol1``, etc. to the different added molecules. The ``region_names`` option lets you set custom names.
+# By default, the ``packmol()`` function assigns regions called ``mol0``, ``mol1``, etc. to the different added molecules. The ``region_names`` option lets you set custom names. 
 
 out = packmol(
     [water, n2],
@@ -333,7 +328,7 @@ def printsummary(mol, details=None):
 print(AMSJob(molecule=out).get_input())
 
 
-# Below, we also set ``keep_atom_properties=False``, this will remove the previous regions (in this example "oxygen_atom") and mass.
+# Below, we also set ``keep_atom_properties=False``, this will remove the previous regions (in this example "oxygen_atom") and mass. 
 
 out = packmol([water, n2], n_molecules=[2, 1], density=0.5, keep_atom_properties=False)
 print(AMSJob(molecule=out).get_input())
@@ -350,3 +345,4 @@ def printsummary(mol, details=None):
     keep_atom_properties=False,
 )
 print(AMSJob(molecule=out).get_input())
+

interfaces/molecule/packmol.py

@@ -13,6 +13,7 @@
 from scm.plams.interfaces.molecule.rdkit import readpdb, writepdb
 from scm.plams.core.functions import requires_optional_package, log
 from scm.plams.core.settings import Settings
+from scm.plams.core.jobmanager import JobManager
 
 if TYPE_CHECKING:
     try:
@@ -816,24 +817,27 @@ def _run_uff_md(
         s.input.ams.MolecularDynamics.Deformation.TargetLattice._1 = target_lattice_str
 
     previous_config = config.copy()
-    config.job.pickle = False
-    config.log.stdout = 0
-    # TODO: fix this so that it doesn't leave plams_workdir on disk if it is created
-    job = AMSJob(settings=s, molecule=md_ucs, name="shakemd")
-    job.run()
-    job.results.wait()
-    config.job.pickle = previous_config.job.pickle
-    config.log.stdout = previous_config.log.stdout
-    if not job.ok():
-        raise PackMolError(
-            f"Try a lower density or a less skewed cell! Original file in {job.path} . "
-            + str(job.results.get_errormsg())
-        )
-    my_packed = job.results.get_main_system()
-    if not keepjob:
-        delete_job(job)
+    try:
+        config.job.pickle = False
+        config.log.stdout = 0
+
+        with tempfile.TemporaryDirectory() as tmp_dir:
+            job_manager = JobManager(config.jobmanager, path=tmp_dir)
+            job = AMSJob(settings=s, molecule=md_ucs, name="shakemd")
+            job.run(jobmanager=job_manager)
+
+            if not job.ok():
+                raise PackMolError(
+                    f"Try a lower density or a less skewed cell! Original file in {job.path} . "
+                    + str(job.results.get_errormsg())
+                )
+            my_packed = job.results.get_main_system()
+            my_packed.remove_region("PACKMOL_thermostatted")
+
+    finally:
+        config.job.pickle = previous_config.job.pickle
+        config.log.stdout = previous_config.log.stdout
 
-    my_packed.remove_region("PACKMOL_thermostatted")
     return my_packed
 
 
@@ -855,7 +859,7 @@ def packmol_around(
 
     For all other arguments, see the ``packmol`` function.
 
-    In the returned ``Molecule`, the system will be mapped to [0..1]. It has the same lattice has ``current``.
+    In the returned ``Molecule``, the system will be mapped to ``[0..1]``. It has the same lattice has ``current``.
     """
     from scm.libbase import (
         UnifiedChemicalSystem as ChemicalSystem,
@@ -924,16 +928,17 @@ def packmol_around(
         supercell.supercell_trafo(trafo)
         supercell.map_atoms(0)
         system_for_packing = supercell
+        tolerance = kwargs.get("tolerance", 1.5)
     else:
         # now distort the original system to the target lattice
         distorted = original_ucs.copy()
         distorted.lattice.vectors = np.diag(maxcomponents)
         distorted.set_fractional_coordinates(original_frac_coords)
         system_for_packing = distorted
+        # in general we need higher tolerance here since we may be expanding the original system,
+        # and we do not want the added molecules to enter in artificial "voids"
+        tolerance = kwargs.get("tolerance", 1.5) * 1.3  # should depend on distortion_vol_expansion_factor somehow
 
-    # in general we need higher tolerance here since we may be expanding the original system,
-    # and we do not want the added molecules to enter in artificial "voids"
-    tolerance = kwargs.get("tolerance", 1.5) * 1.3  # should depend on distortion_vol_expansion_factor somehow
     log(f"{system_for_packing_type=}", loglevel)
     log(f"{n_molecules=}, {box_bounds=}, {tolerance=}", loglevel)
     my_packed, details = packmol(