Fix xtb opt with cartesian coordinates (duartegroup#247)

* Fix xtb opt with cartesian coordinates

* Update changelog [skip actions]

autode/wrappers/XTB.py

@@ -50,29 +50,31 @@ def print_cartesian_constraints(self, inp_file, molecule):
         if molecule.constraints.cartesian is None:
             return None
 
-        constrained_atom_idxs = [i + 1 for i in molecule.constraints.cartesian]
-        list_of_ranges, used_atoms = [], []
-
-        for i in constrained_atom_idxs:
-            atom_range = []
-            if i not in used_atoms:
-                while i in constrained_atom_idxs:
-                    used_atoms.append(i)
-                    atom_range.append(i)
-                    i += 1
-                if len(atom_range) in (1, 2):
-                    list_of_ranges += str(atom_range)
-                else:
-                    list_of_ranges.append(f"{atom_range[0]}-{atom_range[-1]}")
+        atom_idxs = list(
+            sorted(int(i) + 1 for i in molecule.constraints.cartesian)
+        )
+        list_of_ranges = []
 
+        for atom_idx in atom_idxs:
+            last_range = (
+                list_of_ranges[-1] if len(list_of_ranges) > 0 else None
+            )
+            if last_range is not None and atom_idx - 1 == last_range[-1]:
+                last_range.append(atom_idx)
+            else:
+                list_of_ranges.append([atom_idx])
+
+        list_of_ranges_str = [
+            f"{idxs[0]}-{idxs[-1]}" if len(idxs) > 1 else str(idxs[0])
+            for idxs in list_of_ranges
+        ]
         print(
             f"$constrain\n"
             f"force constant={self.force_constant}\n"
-            f'atoms: {",".join(list_of_ranges)}\n'
+            f'atoms: {",".join(list_of_ranges_str)}\n'
             f"$",
             file=inp_file,
         )
-
         return None
 
     @staticmethod

doc/changelog.rst

@@ -20,6 +20,7 @@ Bug Fixes
 *********
 - Fixes :code:`ERROR` logging level being ignored from environment variable :code:`AUTODE_LOG_LEVEL`
 - Fixes :code:`autode.values.Value` instances generating items with units on division, and throw a warning if multiplying
+- Fixes the printing of cartesian constraints in XTB input files, meaning they are no longer ignored
 
 
 1.3.4

tests/test_wrappers/test_xtb.py

@@ -302,6 +302,11 @@ def test_xtb_opt_non_contiguous_range_cart_constraints():
     )
     calc.run()
 
+    assert len(calc.input.additional_filenames) > 0
+    xcontrol_lines = open(calc.input.additional_filenames[-1], "r").readlines()
+    expected_range = "1-3,6"
+    assert sum(expected_range in line for line in xcontrol_lines) == 1
+
     assert calc.optimiser.converged
     assert mol.energy is not None
 
@@ -362,3 +367,23 @@ def run_calc(_mol):
     run_calc(mol)
 
     assert mol.energy != unconstrained_energy
+
+
+@testutils.requires_with_working_xtb_install
+@work_in_tmp_dir()
+def test_xtb_cartesian_constrained_opt():
+
+    init_r = 0.9
+    h2 = Molecule(atoms=[Atom("H"), Atom("H", x=init_r)])
+
+    h2_unconstrained = h2.new_species(name="unconstrained_h2")
+    h2_unconstrained.optimise(method=XTB())
+    # expected minimum for H2 is ~0.77 Å
+    assert abs(h2_unconstrained.distance(0, 1) - init_r) > 0.1
+
+    h2.constraints.cartesian = [0, 1]
+    h2.optimise(method=XTB())
+
+    # if the coordinates are constrained then the distance should be
+    # close to the initial
+    assert abs(h2.distance(0, 1) - init_r) < 0.1