Merge pull request #178 from SCM-NV/DavidOrmrodMorley/molecule_add_ha…

…toms_fix

Add delete_atoms method and use bonding info for add_hatoms SO107

CHANGELOG.md

@@ -17,6 +17,7 @@ This changelog is effective from the 2025 releases.
 * 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.
+* `Molecule.delete_atoms` method to delete multiple atoms with partial success 
 * Example on `MoleculeFormats`
 * Script `generate_example.sh` to generate documentation pages from notebook examples
 * GitHub workflows for CI and publishing to PyPI
@@ -44,6 +45,7 @@ This changelog is effective from the 2025 releases.
 * `SingleJob.load` returns the correctly loaded job
 * `AMSJob.check` handles a `NoneType` status, returning `False`
 * `MultiJob.run` locking resolved when errors raised within `prerun` and `postrun` methods
+* `Molecule.add_hatoms` to use bonding information if available when adding new hydrogen atoms
 
 ### Deprecated
 * `plams` launch script is deprecated in favour of simply running with `amspython`

interfaces/adfsuite/ams.py

@@ -2355,7 +2355,11 @@ def copy_mol(mol):
         super().__init__(molecule, *args, **kwargs)
 
     def run(
-        self, jobrunner: "JobRunner" = None, jobmanager: "JobManager" = None, watch: bool = False, **kwargs
+        self,
+        jobrunner: Optional["JobRunner"] = None,
+        jobmanager: Optional["JobManager"] = None,
+        watch: bool = False,
+        **kwargs,
     ) -> AMSResults:
         """Run the job using *jobmanager* and *jobrunner* (or defaults, if ``None``).
 

mol/molecule.py

@@ -336,6 +336,26 @@ def delete_atom(self, atom: Atom) -> None:
             self.delete_bond(b)
         atom.mol = None
 
+    def delete_atoms(self, atoms: Iterable[Atom]) -> None:
+        """Delete multiple *atom* from the molecule.
+
+        *atom* should be an iterable of |Atom| instances which belong to the molecule. All bonds containing these atoms will be removed too.
+
+        Note that this method employs partial success, such that if deleting any atom results in an error, the remaining atoms will
+        still be deleted. An aggregate error will then be raised at the end of the operation if any errors were encountered.
+        """
+
+        errors = []
+        materialised_atoms = [at for at in atoms]  # make sure to materialise the atoms before starting deletion
+        for atom in materialised_atoms:
+            try:
+                self.delete_atom(atom)
+            except MoleculeError as err:
+                errors.append(f"{err}")
+        if any(errors):
+            error_details = str.join("\n", errors)
+            raise MoleculeError(f"Encountered one or more errors when deleting atoms:\n{error_details}")
+
     def add_bond(self, arg1: Union[Bond, Atom], arg2: Optional[Atom] = None, order: float = 1) -> None:
         """Add a new bond to the molecule.
 
@@ -3341,20 +3361,31 @@ def add_hatoms(self) -> "Molecule":
         from subprocess import DEVNULL, Popen
         from tempfile import NamedTemporaryFile
 
-        with NamedTemporaryFile(mode="w+", suffix=".xyz", delete=False) as f_in:
-            self.writexyz(f_in)
+        # Pass an input file to amsprep which contains current geometry and bonding information
+        with NamedTemporaryFile(mode="w+", suffix=".in", delete=False) as f_in:
+            self.writein(f_in)
             f_in.close()
+
+            # Get .xyz file from amsprep containing the geometry to the same precision (.mol file causes rounding)
+            # And then load the bonding information from the output
             with NamedTemporaryFile(mode="w+", suffix=".xyz", delete=False) as f_out:
-                f_out.close()
-                amsprep = os.path.join(os.environ["AMSBIN"], "amsprep")
-                p = Popen(
-                    f"sh {amsprep} -t SP -m {f_in.name} -addhatoms -exportcoordinates {f_out.name}",
-                    shell=True,
-                    stdout=DEVNULL,
-                )
-                p.communicate()
-                retmol = self.__class__(f_out.name)
-                os.remove(f_out.name)
+                with NamedTemporaryFile(mode="w+", suffix=".out", delete=False) as f_out_bonds:
+                    f_out.close()
+                    f_out_bonds.close()
+                    amsprep = os.path.join(os.environ["AMSBIN"], "amsprep")
+                    p = Popen(
+                        f"sh {amsprep} -t SP -m {f_in.name} -addhatoms -exportcoordinates {f_out.name} -bondsonly > {f_out_bonds.name}",
+                        shell=True,
+                        stdout=DEVNULL,
+                    )
+                    p.communicate()
+                    retmol = self.__class__(f_out.name)
+                    with open(f_out_bonds.name) as bonds_file:
+                        for line in bonds_file:
+                            _, i, j, bo = line.split()
+                            retmol.add_bond(retmol[int(i)], retmol[int(j)], float(bo))
+                    os.remove(f_out.name)
+                    os.remove(f_out_bonds.name)
             os.remove(f_in.name)
         return retmol
 

unit_tests/test_molecule.py

@@ -14,6 +14,7 @@
 from scm.plams.mol.molecule import Molecule, MoleculeError
 from scm.plams.core.functions import read_all_molecules_in_xyz_file
 from scm.plams.interfaces.molecule.rdkit import from_smiles
+from scm.plams.unit_tests.test_helpers import skip_if_no_ams_installation
 
 
 class MoleculeTestBase(ABC):
@@ -131,8 +132,28 @@ def test_add_delete_atoms_and_bonds_happy(self, hydroxide, water):
         assert oh2.mol == peroxide4
         assert o1.mol == peroxide4
 
+        # 1) Add water to water
+        # 2) Delete extra H atoms (with partial success)
+        # 3) Add O-O bond and O-H bonds
+        peroxide5 = water.copy()
+        peroxide5.add_molecule(water, copy=True, margin=1.0)
+
+        def atoms_to_delete():
+            yield peroxide5[3]
+            yield peroxide5[6]
+            yield peroxide5[6]
+
+        with pytest.raises(MoleculeError):
+            peroxide5.delete_atoms(atoms_to_delete())
+
+        peroxide5.add_bond(peroxide5[1], peroxide5[3])
+        peroxide5.add_bond(peroxide5[1], peroxide5[2])
+        peroxide5.add_bond(peroxide5[3], peroxide5[4])
+
         # Assert the same
-        assert peroxide1.label(3) == peroxide2.label(3) == peroxide3.label(3) == peroxide4.label(3)
+        assert (
+            peroxide1.label(3) == peroxide2.label(3) == peroxide3.label(3) == peroxide4.label(3) == peroxide5.label(3)
+        )
 
     def test_add_delete_atoms_and_bonds_unhappy(self, water):
         water2 = water.copy()
@@ -159,6 +180,10 @@ def test_add_delete_atoms_and_bonds_unhappy(self, water):
         with pytest.raises(MoleculeError):
             water.delete_atom(h2)
 
+        # Cannot delete multiple atoms which do not form part of the molecule
+        with pytest.raises(MoleculeError):
+            water.delete_atoms([water2[1], water2[2]])
+
         # Cannot add bond which has invalid arguments
         water2.guess_bonds()
         water = water2.copy()
@@ -216,6 +241,12 @@ def test_system_and_atomic_charge(self, mol):
         with pytest.raises(MoleculeError):
             assert mol.guess_atomic_charges() == [1, 1, 0, 0]
 
+    def test_add_hatoms(self, mol, hydroxide):
+        skip_if_no_ams_installation()
+
+        water = hydroxide.add_hatoms()
+        assert mol.label(4) == water.label(4)
+
     def test_get_complete_molecules_within_threshold(self, mol):
         m0 = mol.get_complete_molecules_within_threshold([2], 0)
         m1 = mol.get_complete_molecules_within_threshold([2], 1)
@@ -589,6 +620,17 @@ def test_locate_rings(self, mol):
         assert_rings_equal(mol.locate_rings_networkx(False), expected)
         assert_rings_equal(mol.locate_rings_networkx(True), expected)
 
+    def test_add_hatoms(self, mol):
+        skip_if_no_ams_installation()
+
+        mol.guess_bonds()
+        mol2 = mol.copy()
+
+        mol2.delete_atoms([at for at in mol2 if at.symbol == "H"])
+        mol2 = mol2.add_hatoms()
+        assert len(mol2.bonds) == len(mol.bonds) == 12
+        assert mol.label(4) == mol2.label(4)
+
 
 def assert_rings_equal(actual, expected):
     assert len(actual) == len(expected)
@@ -919,6 +961,70 @@ def test_locate_rings(self, mol):
         assert_rings_equal(mol.locate_rings_networkx(True), expected)
 
 
+class TestFragments(MoleculeTestBase):
+
+    @pytest.fixture
+    def mol(self, xyz_folder):
+        mol = Molecule(xyz_folder / "C7H8N2_fragments.xyz")
+        return mol
+
+    @property
+    def expected_atoms(self):
+        return [
+            ("H", 37.97589166, 6.50362513, 6.10947932, {}),
+            ("H", 22.94566879, 5.21944993, 6.11208918, {}),
+            ("H", 24.24811668, 3.58937684, 3.1101599, {}),
+            ("H", 26.46029574, 4.22221813, 1.9970365, {}),
+            ("H", 28.11383033, 5.49145028, 3.3383532, {}),
+            ("H", 33.53518443, 7.87366456, 5.79895021, {}),
+            ("C", 36.81091549, 5.79838054, 2.99659108, {}),
+            ("N", 35.79722124, 6.77372035, 4.90114805, {}),
+            ("H", 27.40840017, 6.85618188, 6.76571824, {}),
+            ("C", 35.70493885, 6.05732335, 2.17566355, {}),
+            ("C", 34.68676663, 6.93124889, 4.1845686, {}),
+            ("C", 34.61783813, 6.60177631, 2.80991515, {}),
+            ("N", 25.60906993, 5.44707201, 5.57733345, {}),
+            ("N", 32.40678698, 7.45865484, 4.11253188, {}),
+            ("N", 28.93739263, 6.57981707, 5.45707573, {}),
+            ("H", 37.67050745, 5.3014948, 2.59000904, {}),
+            ("H", 35.78463338, 5.85496894, 1.08423825, {}),
+            ("H", 33.70046648, 6.8213344, 2.30167641, {}),
+            ("C", 37.95605417, 5.90776431, 5.20779435, {}),
+            ("C", 23.41069028, 4.45252478, 5.49358156, {}),
+            ("C", 24.74504443, 4.71605898, 4.8543161, {}),
+            ("C", 25.04631067, 4.17412898, 3.6262066, {}),
+            ("C", 26.2268722, 4.52245672, 2.98318138, {}),
+            ("C", 26.77975995, 5.66497554, 5.02272776, {}),
+            ("C", 27.10838378, 5.26861749, 3.73156573, {}),
+            ("C", 33.47999143, 7.48497558, 4.78073603, {}),
+            ("C", 36.82761129, 6.16730707, 4.30869146, {}),
+            ("C", 27.73379633, 6.4261697, 5.80579523, {}),
+            ("H", 29.31659198, 7.50946073, 5.34720982, {}),
+            ("H", 32.44075585, 7.41892399, 3.10388547, {}),
+        ]
+
+    def test_add_hatoms(self, mol):
+        skip_if_no_ams_installation()
+
+        # Given two fragments with no bond information
+        # When add H atoms
+        mol2 = mol.add_hatoms()
+
+        # Then adds hydrogens according to guessed bonds
+        assert mol2.get_formula() == "C14H15N4"
+
+        # But given bonding information
+        mol3 = mol.copy()
+        mol3.guess_bonds()
+        mol3[19, 27].order = 1
+
+        # When add H atoms
+        mol4 = mol3.add_hatoms()
+
+        # Then adds hydrogens which respect modified bond orders
+        assert mol4.get_formula() == "C14H16N4"
+
+
 def test_read_multiple_molecules_from_xyz(xyz_folder):
     """Test for read_all_molecules_in_xyz_file"""
 

unit_tests/xyz/C7H8N2_fragments.xyz

@@ -0,0 +1,32 @@
+30
+
+H      37.97589166       6.50362513       6.10947932
+H      22.94566879       5.21944993       6.11208918
+H      24.24811668       3.58937684       3.11015990
+H      26.46029574       4.22221813       1.99703650
+H      28.11383033       5.49145028       3.33835320
+H      33.53518443       7.87366456       5.79895021
+C      36.81091549       5.79838054       2.99659108
+N      35.79722124       6.77372035       4.90114805
+H      27.40840017       6.85618188       6.76571824
+C      35.70493885       6.05732335       2.17566355
+C      34.68676663       6.93124889       4.18456860
+C      34.61783813       6.60177631       2.80991515
+N      25.60906993       5.44707201       5.57733345
+N      32.40678698       7.45865484       4.11253188
+N      28.93739263       6.57981707       5.45707573
+H      37.67050745       5.30149480       2.59000904
+H      35.78463338       5.85496894       1.08423825
+H      33.70046648       6.82133440       2.30167641
+C      37.95605417       5.90776431       5.20779435
+C      23.41069028       4.45252478       5.49358156
+C      24.74504443       4.71605898       4.85431610
+C      25.04631067       4.17412898       3.62620660
+C      26.22687220       4.52245672       2.98318138
+C      26.77975995       5.66497554       5.02272776
+C      27.10838378       5.26861749       3.73156573
+C      33.47999143       7.48497558       4.78073603
+C      36.82761129       6.16730707       4.30869146
+C      27.73379633       6.42616970       5.80579523
+H      29.31659198       7.50946073       5.34720982
+H      32.44075585       7.41892399       3.10388547