ONIOM capping with chemical groups (#141)

* Implementation for other chemical groups capping in ONIOM. Tests, docs.

tangelo/molecule_library.py

@@ -126,3 +126,16 @@
     ("H", (-1.1260, -0.2065,  2.4095)),
     ("H", ( 4.1118, -0.2131, -1.6830))
 ]
+
+
+# Ethane (H3CCH3).
+xyz_ethane = [
+    ("C", ( 0.7166,  0.8980,  0.6425)),
+    ("H", ( 0.5397,  1.7666, -0.0025)),
+    ("H", ( 0.4899,  0.0005,  0.0551)),
+    ("H", (-0.0078,  0.9452,  1.4640)),
+    ("C", ( 2.1541,  0.8745,  1.1659)),
+    ("H", ( 2.3313,  0.0053,  1.8100)),
+    ("H", ( 2.8785,  0.8284,  0.3444)),
+    ("H", ( 2.3805,  1.7715,  1.7542))
+]

tangelo/problem_decomposition/oniom/_helpers/capping_groups.py

@@ -0,0 +1,58 @@
+# Copyright 2021 Good Chemistry Company.
+#
+# Licensed under the Apache License, Version 2.0 (the "License");
+# you may not use this file except in compliance with the License.
+# You may obtain a copy of the License at
+#
+#   http://www.apache.org/licenses/LICENSE-2.0
+#
+# Unless required by applicable law or agreed to in writing, software
+# distributed under the License is distributed on an "AS IS" BASIS,
+# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+# See the License for the specific language governing permissions and
+# limitations under the License.
+
+"""Module to regroup the list of possible atoms and built-in chemical groups for
+broken link capping. The geometries for the capping groups are taken from the
+NIST Chemistry WebBook database (https://webbook.nist.gov/chemistry/).
+"""
+
+
+elements = [
+    "X",
+    "H", "He", "Li", "Be",  "B",  "C",  "N",  "O",  "F", "Ne",
+    "Na", "Mg", "Al", "Si",  "P",  "S", "Cl", "Ar",  "K", "Ca",
+    "Sc", "Ti",  "V", "Cr", "Mn", "Fe", "Co", "Ni", "Cu", "Zn",
+    "Ga", "Ge", "As", "Se", "Br", "Kr", "Rb", "Sr",  "Y", "Zr",
+    "Nb", "Mo", "Tc", "Ru", "Rh", "Pd", "Ag", "Cd", "In", "Sn",
+    "Sb", "Te",  "I", "Xe", "Cs", "Ba", "La", "Ce", "Pr", "Nd",
+    "Pm", "Sm", "Eu", "Gd", "Tb", "Dy", "Ho", "Er", "Tm", "Yb",
+    "Lu", "Hf", "Ta",  "W", "Re", "Os", "Ir", "Pt", "Au", "Hg",
+    "Tl", "Pb", "Bi", "Po", "At", "Rn", "Fr", "Ra", "Ac", "Th",
+    "Pa",  "U", "Np", "Pu", "Am", "Cm", "Bk", "Cf", "Es", "Fm",
+    "Md", "No", "Lr", "Rf", "Db", "Sg", "Bh", "Hs", "Mt", "Ds",
+    "Rg", "Cn", "Nh", "Fl", "Mc", "Lv", "Ts", "Og"
+]
+
+chemical_groups = {
+    "CH3": [  # From ethane.
+        ["X", [  2.15410,  0.87450,  1.16590]],
+        ["C", [  0.71660,  0.89800,  0.64250]],
+        ["H", [  0.53970,  1.76660, -0.00250]],
+        ["H", [  0.48990,  0.00050,  0.05510]],
+        ["H", [ -0.00780,  0.94520,  1.46400]]
+    ],
+    "CF3": [  # From hexafluoroethane.
+        ["X", [  2.03950,  1.01440,  0.00010]],
+        ["C", [  0.49550,  1.01460, -0.00010]],
+        ["F", [  0.04100, -0.24410,  0.00000]],
+        ["F", [  0.04120,  1.64380, -1.09040]],
+        ["F", [  0.04080,  1.64420,  1.08980]]
+    ],
+    "NH2": [  # From ammonia.
+        ["X", [ -0.01300, -0.00270,  0.00580]],
+        ["N", [  0.87140,  0.50430, -0.01000]],
+        ["H", [  1.59990, -0.20890,  0.00580]],
+        ["H", [  0.92660,  0.93550, -0.93240]]
+    ]
+}

tangelo/problem_decomposition/oniom/_helpers/helper_classes.py

@@ -17,11 +17,14 @@
 which bonds are broken and how to fix them) as well as the solver(s) to use.
 """
 
+import warnings
+
 import numpy as np
+from scipy.spatial.transform import Rotation as R
 
-# Imports of electronic solvers and data structure
-from tangelo.algorithms import CCSDSolver, FCISolver, VQESolver, MINDO3Solver
 from tangelo import SecondQuantizedMolecule
+from tangelo.algorithms import CCSDSolver, FCISolver, VQESolver, MINDO3Solver
+from tangelo.problem_decomposition.oniom._helpers.capping_groups import elements, chemical_groups
 
 
 class Fragment:
@@ -212,19 +215,27 @@ def __init__(self, staying, leaving, factor=1.0, species="H"):
         methods to generate a new bond, appending the intended species.
 
         Args:
-            index1 (int): Order in the molecular geometry of atom retained in
-                model-unit.
-            leaving (int): Order in mol. Geometry of atom lost.
+            staying (int): Atom id retained.
+            leaving (int): Atom id lost.
             factor (float) optional: Rescale length of bond, from that in the
                 original molecule.
-            species (str) optional: Atomic species of appended atom for new
-                link.
+            species (str) optional: Atomic species or a chemical group
+                identifier for new link. Can be a list (first element = "X" to
+                detect the orientation) for a custom chemical group.
         """
 
         self.staying = staying
         self.leaving = leaving
         self.factor = factor
-        self.species = species
+
+        if isinstance(species, str) and species in elements:
+            self.species = [(species, (0., 0., 0.))]
+        elif isinstance(species, str) and species in chemical_groups:
+            self.species = chemical_groups[species]
+        elif isinstance(species, (list, tuple)) and species[0][0].upper() == "X":
+            self.species = species
+        else:
+            raise ValueError(f"{species} is not supported. It must be a string identifier or a list of atoms (with a ghost atom ('X') as the first element).")
 
     def relink(self, geometry):
         """Create atom at location of mended-bond link.
@@ -234,12 +245,29 @@ def relink(self, geometry):
                 [[str,tuple(float,float,float)],...].
 
         Returns:
-            str: Atomic species.
-            tuple: Position (x, y, z) of replacement atom.
+            list: List of atomic species and position (x, y, z) of replacement
+                atom / chemical group.
         """
 
+        elements = [a[0] for a in self.species if a[0].upper() != "X"]
+        chem_group_xyz = np.array([[a[1][0], a[1][1], a[1][2]] for a in self.species if a[0].upper() != "X"])
+
         staying = np.array(geometry[self.staying][1])
         leaving = np.array(geometry[self.leaving][1])
+
+        # Rotation (if not a single atom).
+        if len(elements) > 1:
+            axis_old = leaving - staying
+            axis_new = chem_group_xyz[0] - np.array(self.species[0][1])
+
+            with warnings.catch_warnings():
+                warnings.simplefilter("ignore", UserWarning)
+                rot, _ = R.align_vectors([axis_old], [axis_new])
+            chem_group_xyz = rot.apply(chem_group_xyz)
+
+        # Move the atom / group to the right position in space.
         replacement = self.factor*(leaving-staying) + staying
+        translation = replacement - chem_group_xyz[0]
+        chem_group_xyz += translation
 
-        return (self.species, (replacement[0], replacement[1], replacement[2]))
+        return [(element, (xyz[0], xyz[1], xyz[2])) for element, xyz in zip(elements, chem_group_xyz)]

tangelo/problem_decomposition/oniom/oniom_problem_decomposition.py

@@ -30,7 +30,6 @@
     Chemical Reviews 2015 115 (12), 5678-5796. DOI: 10.1021/cr5004419.
 """
 # TODO: Supporting many (3+) layers of different accuracy.
-# TODO: Capping with CH3 or other functional groups.
 
 from tangelo.problem_decomposition.problem_decomposition import ProblemDecomposition
 from tangelo.toolboxes.molecular_computation.molecule import atom_string_to_list
@@ -95,9 +94,10 @@ def distribute_atoms(self):
 
             # If there are broken_links (other than an empty list nor None).
             # The whole molecule geometry is needed to compute the position of
-            # the capping atom (or functional group in the future).
+            # the capping atom (or functional group).
             if fragment.broken_links:
-                fragment.geometry += [li.relink(self.geometry) for li in fragment.broken_links]
+                for li in fragment.broken_links:
+                    fragment.geometry += li.relink(self.geometry)
 
     def simulate(self):
         r"""Run the ONIOM core-method. The total energy is defined as

tangelo/problem_decomposition/tests/oniom/test_capping.py

@@ -0,0 +1,114 @@
+# Copyright 2021 Good Chemistry Company.
+#
+# Licensed under the Apache License, Version 2.0 (the "License");
+# you may not use this file except in compliance with the License.
+# You may obtain a copy of the License at
+#
+#   http://www.apache.org/licenses/LICENSE-2.0
+#
+# Unless required by applicable law or agreed to in writing, software
+# distributed under the License is distributed on an "AS IS" BASIS,
+# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+# See the License for the specific language governing permissions and
+# limitations under the License.
+
+import unittest
+
+import numpy as np
+
+from tangelo.problem_decomposition.oniom._helpers.helper_classes import Fragment, Link
+from tangelo.molecule_library import xyz_ethane
+from tangelo.toolboxes.molecular_computation.molecule import atom_string_to_list
+
+
+class ONIOMCappingTest(unittest.TestCase):
+
+    def test_unsupported_string_species(self):
+        """Test unsupported built-in chemical group (raise ValueError)."""
+
+        with self.assertRaises(ValueError):
+            Link(0, 4, factor=1., species="UNSUPPORTED")
+
+    def test_unsupported_custom_species(self):
+        """Test unsupported custom chemical group (raise ValueError)."""
+
+        # No ghost atom (X) as first element.
+        not_supported_chem_group = [
+            ["C", [ 0.71660,  0.89800,  0.64250]],
+            ["I", [ 0.53970,  1.76660, -0.00250]],
+            ["Cl", [ 0.48990,  0.00050,  0.05510]],
+            ["F", [-0.00780,  0.94520,  1.46400]]
+        ]
+
+        with self.assertRaises(ValueError):
+            Link(0, 4, factor=1., species=not_supported_chem_group)
+
+    def test_cf3(self):
+        """Test capping with CF3 (built-in)."""
+
+        broken_link = Link(0, 4, factor=1., species="CF3")
+        new_pos = broken_link.relink(xyz_ethane)
+
+        # First CH3 + capping group.
+        xyz = xyz_ethane[:4] + new_pos
+
+        # Those position have been verified with a molecular visualization software.
+        ref_capped = atom_string_to_list("""
+            C     0.71660    0.89800    0.64250
+            H     0.53970    1.76660   -0.00250
+            H     0.48990    0.00050    0.05510
+            H    -0.00780    0.94520    1.46400
+            C     2.15410    0.87450    1.16590
+            F     2.20191    1.41954    2.38719
+            F     2.94837    1.57115    0.34441
+            F     2.59308   -0.38799    1.23257
+        """)
+
+        # Every atom must be the same (same order too).
+        # Position can be almost equals.
+        for i, atom in enumerate(xyz):
+            self.assertEqual(atom[0], ref_capped[i][0])
+
+        # Testing only the position of the central atom in the chemical group.
+        np.testing.assert_almost_equal(xyz[4][1], ref_capped[4][1], decimal=4)
+
+    def test_custom_species(self):
+        """Test capping with a custom chemical group (here -CFICl)."""
+
+        cap = [
+            ["X", [ 2.15410,  0.87450,  1.16590]],
+            ["C", [ 0.71660,  0.89800,  0.64250]],
+            ["I", [ 0.53970,  1.76660, -0.00250]],
+            ["Cl", [ 0.48990,  0.00050,  0.05510]],
+            ["F", [-0.00780,  0.94520,  1.46400]]
+        ]
+
+        broken_link = Link(0, 4, factor=1., species=cap)
+        new_pos = broken_link.relink(xyz_ethane)
+
+        # First CH3 + capping group.
+        xyz = xyz_ethane[:4] + new_pos
+
+        # Those position have been verified with a molecular visualization software.
+        ref_capped = atom_string_to_list("""
+            C     0.71660    0.89800    0.64250
+            H     0.53970    1.76660   -0.00250
+            H     0.48990    0.00050    0.05510
+            H    -0.00780    0.94520    1.46400
+            C     2.15410    0.87450    1.16590
+            I     2.22665    0.32520    2.11184
+            Cl    2.83194    0.39157    0.45228
+            F     2.53182    1.88810    1.34418
+        """)
+
+        # Every atom must be the same (same order too).
+        # Position can be almost equals.
+        for i, atom in enumerate(xyz):
+            self.assertEqual(atom[0], ref_capped[i][0])
+
+        # Testing only the position of the central atom in the chemical group.
+        np.testing.assert_almost_equal(xyz[4][1], ref_capped[4][1], decimal=4)
+
+
+if __name__ == "__main__":
+    unittest.main()

tangelo/problem_decomposition/tests/oniom/test_oniom.py

@@ -80,10 +80,10 @@ def test_capping_broken_link(self):
             self.assertEqual(atom[0], PHE_backbone_capped[i][0])
 
             for dim in range(3):
-                self.assertAlmostEqual(atom[1][dim], PHE_backbone_capped[i][1][dim])
+                self.assertAlmostEqual(atom[1][dim], PHE_backbone_capped[i][1][dim], places=4)
 
     def test_energy(self):
-        """Testing the oniom energy with a low accuraccy method (RHF) and an
+        """Testing the oniom energy with a low accuracy method (RHF) and an
         higher one (CCSD) for PHE molecule. The important fragment is chosen to
         be the backbone. The side chain is computed at the RHF level.
         """