From 6815b891fb9485c8a8bd1cdb6d3883c059c61cfc Mon Sep 17 00:00:00 2001 From: Lester Hedges Date: Thu, 20 Nov 2025 13:47:20 +0000 Subject: [PATCH] Backport fix from PR #471. [ci skip] --- python/BioSimSpace/FreeEnergy/_atm.py | 79 ++++++++++++++++++++++++ python/BioSimSpace/Process/_atm_utils.py | 30 +++++++-- 2 files changed, 103 insertions(+), 6 deletions(-) diff --git a/python/BioSimSpace/FreeEnergy/_atm.py b/python/BioSimSpace/FreeEnergy/_atm.py index dc7a997d..bf245bc9 100644 --- a/python/BioSimSpace/FreeEnergy/_atm.py +++ b/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( diff --git a/python/BioSimSpace/Process/_atm_utils.py b/python/BioSimSpace/Process/_atm_utils.py index 6af19b9a..08311fa5 100644 --- a/python/BioSimSpace/Process/_atm_utils.py +++ b/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"