Merge pull request #272 from OpenBioSim/feature_amber_fep

Add support for AMBER as a free-energy perturbation engine.
[ci skip]

doc/source/tutorials/hydration_freenrg.rst

@@ -167,6 +167,20 @@ Let's examine the directory for the :math:`{\lambda=0}` window of the free leg:
    gromacs.err  gromacs.mdp  gromacs.out.mdp  gromacs.tpr
    gromacs.gro  gromacs.out  gromacs.top
 
+Similarly, we can also set up the same simulations using AMBER:
+
+.. code-block:: python
+
+   free_amber = BSS.FreeEnergy.Relative(solvated, protocol, engine="amber", work_dir="freenrg_amber/free")
+   vac_amber = BSS.FreeEnergy.Relative(merged.toSystem(), protocol, engine="amber", work_dir="freenrg_amber/vacuum")
+
+Let's examine the directory for the :math:`{\lambda=0}` window of the free leg:
+
+.. code-block:: bash
+
+   $ ls freenrg_amber/free/lambda_0.0000
+   amber.cfg  amber.err  amber.out  amber.prm7  amber_ref.rst7  amber.rst7
+
 There you go! This tutorial has shown you how BioSimSpace can be used to easily set
 up everything that is needed for complex alchemical free energy simulations. Please
 visit the :data:`API documentation <BioSimSpace.FreeEnergy>` for further information.

python/BioSimSpace/Align/_merge.py

@@ -27,6 +27,7 @@
 __all__ = ["merge"]
 
 from sire.legacy import Base as _SireBase
+from sire.legacy import IO as _SireIO
 from sire.legacy import MM as _SireMM
 from sire.legacy import Mol as _SireMol
 from sire.legacy import Units as _SireUnits
@@ -1389,3 +1390,86 @@ def _is_on_ring(idx, conn):
 
     # If we get this far, then the atom is not adjacent to a ring.
     return False
+
+
+def _removeDummies(molecule, is_lambda1):
+    """
+    Internal function which removes the dummy atoms from one of the endstates
+    of a merged molecule.
+
+    Parameters
+    ----------
+
+    molecule : BioSimSpace._SireWrappers.Molecule
+        The molecule.
+
+    is_lambda1 : bool
+       Whether to use the molecule at lambda = 1.
+    """
+    if not molecule._is_perturbable:
+        raise _IncompatibleError("'molecule' is not a perturbable molecule")
+
+    # Always use the coordinates at lambda = 0.
+    coordinates = molecule._sire_object.property("coordinates0")
+
+    # Generate a molecule with all dummies present.
+    molecule = molecule.copy()._toRegularMolecule(
+        is_lambda1=is_lambda1, generate_intrascale=True
+    )
+
+    # Remove the parameters property, if it exists.
+    if "parameters" in molecule._sire_object.propertyKeys():
+        molecule._sire_object = (
+            molecule._sire_object.edit().removeProperty("parameters").commit()
+        )
+
+    # Set the coordinates to those at lambda = 0
+    molecule._sire_object = (
+        molecule._sire_object.edit().setProperty("coordinates", coordinates).commit()
+    )
+
+    # Extract all the nondummy indices
+    nondummy_indices = [
+        i
+        for i, atom in enumerate(molecule.getAtoms())
+        if "du" not in atom._sire_object.property("ambertype")
+    ]
+
+    # Create an AtomSelection.
+    selection = molecule._sire_object.selection()
+
+    # Unselect all of the atoms.
+    selection.selectNone()
+
+    # Now add all of the nondummy atoms.
+    for idx in nondummy_indices:
+        selection.select(_SireMol.AtomIdx(idx))
+
+    # Create a partial molecule and extract the atoms.
+    partial_molecule = (
+        _SireMol.PartialMolecule(molecule._sire_object, selection).extract().molecule()
+    )
+
+    # Remove the incorrect intrascale property.
+    partial_molecule = (
+        partial_molecule.edit().removeProperty("intrascale").molecule().commit()
+    )
+
+    # Recreate a BioSimSpace molecule object.
+    molecule = _Molecule(partial_molecule)
+
+    # Parse the molecule as a GROMACS topology, which will recover the intrascale
+    # matrix.
+    gro_top = _SireIO.GroTop(molecule.toSystem()._sire_object)
+
+    # Convert back to a Sire system.
+    gro_sys = gro_top.toSystem()
+
+    # Add the intrascale property back into the merged molecule.
+    edit_mol = molecule._sire_object.edit()
+    edit_mol = edit_mol.setProperty(
+        "intrascale", gro_sys[_SireMol.MolIdx(0)].property("intrascale")
+    )
+    molecule = _Molecule(edit_mol.commit())
+
+    return molecule