Merge pull request #471 from OpenBioSim/fix_470

Fix for ATM centre of mass restraints

Author: lohedges

python/BioSimSpace/FreeEnergy/_atm.py

@@ -525,6 +525,8 @@ def prepare(
             self._setLig2ComAtoms(ligand_free_com_atoms)
 
             self._findAtomIndices()
+            self._setLig1ProteinDisplacement()
+            self._setLig2ProteinDisplacement()
             self._makeData()
             serialisable_disp = [
                 self._displacement.x(),
@@ -555,6 +557,8 @@ def prepare(
             self._setLig1ComAtoms(ligand_bound_com_atoms)
             self._setLig2ComAtoms(ligand_free_com_atoms)
             self._findAtomIndices()
+            self._setLig1ProteinDisplacement()
+            self._setLig2ProteinDisplacement()
             self._makeData()
             serialisable_disp = [
                 self._displacement.x(),
@@ -878,6 +882,79 @@ def _setLig2ComAtoms(self, lig2_com_atoms):
                 for a in ligand_free._sire_object[f"atoms within 11 angstrom of {com}"]
             ]
 
+    @staticmethod
+    def _find_separation(com1, com2):
+        """
+        Finds the separation vector between two sets of coordinates.
+
+        Parameters
+        ----------
+        com1 : Coordinates
+            The coordinates of the first molecule of interest.
+        com2 : Coordinates
+            The coordinates of the second molecule of interest.
+
+        Returns
+        -------
+        list
+            A list of floats representing the separation vector between the two molecules.
+        """
+        separation_vector = com2 - com1
+        separation = [i.value() for i in separation_vector]
+        return separation
+
+    def _getLig1ProteinDisplacement(self):
+        """
+        Get the calculated displacement between ligand 1 and the protein.
+
+        Returns
+        -------
+        list
+            A list of floats representing the displacement between ligand 1 and the protein.
+        """
+        return self._lig1_protein_displacement
+
+    def _setLig1ProteinDisplacement(self):
+        """
+        Set the calculated displacement between ligand 1 and the protein.
+        Must be called after MakeSystemFromThree and setting of indices.
+        """
+        # find the ligand atoms that define its center of mass
+        lig1 = self._system[self._ligand_bound_index]
+        lig1_com_coords = lig1._sire_object.atoms()[self._lig1_com_atoms].coordinates()
+        # protein com coords (assumes that the first index in protein_index is representative of the whole protein)
+        prot = self._system[self.protein_index[0]]
+        prot_com_coords = prot._sire_object.atoms()[self._mol1_com_atoms].coordinates()
+        self._lig1_protein_displacement = self._find_separation(
+            prot_com_coords, lig1_com_coords
+        )
+
+    def _getLig2ProteinDisplacement(self):
+        """
+        Get the calculated displacement between ligand 2 and the protein.
+
+        Returns
+        -------
+        list
+            A list of floats representing the displacement between ligand 2 and the protein.
+        """
+        return self._lig2_protein_displacement
+
+    def _setLig2ProteinDisplacement(self):
+        """
+        Set the calculated displacement between ligand 2 and the protein.
+        Must be called after MakeSystemFromThree and setting of indices.
+        """
+        # find the ligand atoms that define its center of mass
+        lig2 = self._system[self._ligand_free_index]
+        lig2_com_coords = lig2._sire_object.atoms()[self._lig2_com_atoms].coordinates()
+        # protein com coords (assumes that the first index in protein_index is representative of the whole protein)
+        prot = self._system[self.protein_index[0]]
+        prot_com_coords = prot._sire_object.atoms()[self._mol1_com_atoms].coordinates()
+        self._lig2_protein_displacement = self._find_separation(
+            prot_com_coords, lig2_com_coords
+        )
+
     def _makeData(self):
         """
         Make the data dictionary for the ATM system
@@ -901,6 +978,8 @@ def _makeData(self):
         self.data["protein_com_atoms"] = self._mol1_com_atoms
         self.data["ligand_bound_com_atoms"] = self._lig1_com_atoms
         self.data["ligand_free_com_atoms"] = self._lig2_com_atoms
+        self.data["lig1_protein_displacement"] = self._lig1_protein_displacement
+        self.data["lig2_protein_displacement"] = self._lig2_protein_displacement
 
     @staticmethod
     def viewRigidCores(

python/BioSimSpace/Process/_atm_utils.py

@@ -72,6 +72,10 @@ def findAbsoluteCOMAtoms(self):
             self.lig2_first_atomnum, self.data["ligand_free_com_atoms"]
         ).tolist()
 
+    def findLigProtDisplacements(self):
+        self.lig1_protein_displacement = self.data["lig1_protein_displacement"]
+        self.lig2_protein_displacement = self.data["lig2_protein_displacement"]
+
     def getATMForceConstants(self, index=None):
         from .. import Protocol as _Protocol
 
@@ -404,6 +408,7 @@ def createCOMRestraint(self, force_group=None):
         """
 
         self.findAbsoluteCOMAtoms()
+        self.findLigProtDisplacements()
         # Groups contained within the constraint
         protein_com = self.protein_com_atoms
         lig1_com = self.lig1_com_atoms
@@ -416,6 +421,19 @@ def createCOMRestraint(self, force_group=None):
         output += "protein_com = {}\n".format(protein_com)
         output += "lig1_com = {}\n".format(lig1_com)
         output += "lig2_com = {}\n".format(lig2_com)
+
+        output += "# Displacement values for protein-ligand COM restraint\n"
+        # round to 3 decimal places for clarity
+        output += "displacement_bound = {}\n".format(
+            [round(x, 3) for x in self.lig1_protein_displacement]
+        )
+        output += "displacement_free = {}\n".format(
+            [round(x, 3) for x in self.lig2_protein_displacement]
+        )
+        output += "# Convert displacement to nm from angstrom\n"
+        output += "displacement_bound = [x * 0.1 for x in displacement_bound]\n"
+        output += "displacement_free = [x * 0.1 for x in displacement_free]\n"
+
         output += "# Constants for the CM-CM force in their input units\n"
         output += "kfcm = {} * kilocalorie_per_mole / angstrom**2\n".format(kf_cm)
         output += "tolcm = {} * angstrom \n".format(tol_cm)
@@ -434,9 +452,9 @@ def createCOMRestraint(self, force_group=None):
         output += """parameters_bound = (
         kfcm.value_in_unit(kilojoules_per_mole / nanometer**2),
         tolcm.value_in_unit(nanometer),
-        0.0 * nanometer,
-        0.0 * nanometer,
-        0.0 * nanometer,
+        displacement_bound[0] * nanometer,
+        displacement_bound[1] * nanometer,
+        displacement_bound[2] * nanometer,
         )\n"""
         output += "force_CMCM.addBond((1,0), parameters_bound)\n"
         output += "numgroups = force_CMCM.getNumGroups()\n"
@@ -446,9 +464,9 @@ def createCOMRestraint(self, force_group=None):
         output += """parameters_free = (
         kfcm.value_in_unit(kilojoules_per_mole / nanometer**2),
         tolcm.value_in_unit(nanometer),
-        displacement[0] * nanometer,
-        displacement[1] * nanometer,
-        displacement[2] * nanometer,
+        displacement_free[0] * nanometer,
+        displacement_free[1] * nanometer,
+        displacement_free[2] * nanometer,
         )\n"""
 
         output += "force_CMCM.addBond((numgroups+1,numgroups+0), parameters_free)\n"