Step-by-step custom constraints implementation

Users can customize various types of constraints, such as distance constraints, angle constraints, dihedral angle constraints, hydrogen bond constraints, etc

Some of the common scoring functions could be:

distance constraints
angle constraints
hydrogen bond constraints

The current framework implements distance constraints. In order to guide users to implement custom constraints, taking the example of adding angle constraints, in the following part,we will demonstrate in detail how to implement angle constraints

1. Define angle constraint class–AngleConstraintSF

The added custom constraints can be understood as a new scoring function, so first go to the scorer folder and open constraints.py

vim opendock/scorer/constraints.py

In the constrains.py file, it can be seen that distance constraints have been implemented, such as DistanceConstraintSF and DistanceMatrixConstraintSF. To add angle constraints, the first step is to implement an angle calculation function located in the base class ConstraintSF. The input is 3 coordinate points, and for ease of calculation, the output is the radian value of the angle.

def _angle(self, x, y, z):
    # Calculate the two vectors
    vec1 = x - y
    vec2 = z - y

    # Compute the dot product of the two vectors
    dot_product = torch.sum(vec1 * vec2)

    # Compute the magnitudes of the vectors
    mag1 = torch.sqrt(torch.sum(torch.pow(vec1, 2)))
    mag2 = torch.sqrt(torch.sum(torch.pow(vec2, 2)))

    # Calculate the cosine of the angle using the dot product formula
    cos_theta = dot_product / (mag1 * mag2)

    # Ensure cos_theta is within the valid range [-1, 1] to avoid numerical errors
    cos_theta = torch.clamp(cos_theta, -1.0, 1.0)

    # Calculate the angle (in radians) using the arccosine function
    angle = torch.acos(cos_theta)

    return angle

Then we can define AngleConstraintSF, which inherits from the base class ConstraintSF

class AngleConstraintSF(ConstraintSF):

  def __init__(self,
             receptor = None,
             ligand = None,
             **kwargs):
    super(AngleConstraintSF, self)\
    .__init__(receptor=receptor, ligand=ligand)

    self.grpA_mol_ = kwargs.pop('groupA_mol', "receptor")

    self.grpB_mol_ = kwargs.pop('groupB_mol', "receptor") # receptor or ligand

    self.grpC_mol_ = kwargs.pop('groupC_mol', "ligand")

    self.grpA_idx_ = kwargs['grpA_ha_indices']
    self.grpB_idx_ = kwargs['grpB_ha_indices']
    self.grpC_idx_ = kwargs['grpC_ha_indices']

    self.constraint_type_ = kwargs.pop('constraint', 'harmonic')
    self.force_constant_ = kwargs.pop('force', 1.0)
    #self.constraint_reference_ = kwargs.pop('reference', None)
    # angle boundary, unit is angstrom
    self.bounds_ = kwargs.pop('bounds', [0, np.pi])


    assert (len(self.grpA_idx_) > 0 and len(self.grpB_idx_) > 0 and len(self.grpC_idx_) > 0)

  def scoring(self):

    if self.grpA_mol_.lower() in ['receptor', 'protein']:
        _grpA_xyz = self.receptor.rec_heavy_atoms_xyz
    elif self.grpA_mol_.lower() in ['ligand', 'molecule']:
        _grpA_xyz = self.ligand.pose_heavy_atoms_coords[0]


    if self.grpB_mol_.lower() in ['receptor', 'protein']:
        _grpB_xyz = self.receptor.rec_heavy_atoms_xyz
    elif self.grpB_mol_.lower() in ['ligand', 'molecule']:
        _grpB_xyz = self.ligand.pose_heavy_atoms_coords[0]

    if self.grpC_mol_.lower() in ['receptor', 'protein']:
        _grpC_xyz = self.receptor.rec_heavy_atoms_xyz
    elif self.grpC_mol_.lower() in ['ligand', 'molecule']:
        _grpC_xyz = self.ligand.pose_heavy_atoms_coords[0]

    pairs = list(zip(self.grpA_idx_, self.grpB_idx_,self.grpC_idx_))

    self._angle_paired_ = []
    for i, (atm1, atm2,atm3) in enumerate(pairs):
        _a = self._angle(_grpA_xyz[atm1], _grpB_xyz[atm2], _grpC_xyz[atm3])
        self._angle_paired_.append(_a)

    self.angle_paired = torch.stack(self._angle_paired_)

    score = self._apply_constraint(torch.mean(self.angle_paired))

    return score.reshape((1, -1))

Thus, we successfully defined the angle constraint through the class AngleConstrainSF. Next, we will demonstrate how to use the newly defined angle constraint

2. Use newly defined angle constraints

Consistent with previous constraint cases, taking pdb: 3gzj as an example, we chose ligand atom CAF, protein side chain atoms ser87-OG, and GLY19-O as cases

Note

Please note that the residue index (residx) is generally 1-based as indicated in the PDB file. The above atomic names have some differences between protein PDB files and PDBQT files, but the atoms are the same.

vim opendock/example/3gzj/atom_angle_constraint_example.py

indices_r1 = asl.select_atom(atomnames=['O', ], chains=['A'], residx=['10'], resnames=['GLY'])
print("indices_r1",indices_r1)
print(indices_r1, receptor.dataframe_ha_.head())

indices_r2 = asl.select_atom(atomnames=['OG', ], chains=['A'], residx=['78'], resnames=['SER'])
asl = AtomSelection(molecule=ligand)

indices_l = asl.select_atom(atomnames=['CAF', ])

cnstr = AngleConstraintSF(receptor, ligand,
                             grpA_ha_indices=indices_r1,
                             grpB_ha_indices=indices_r2,
                             grpC_ha_indices=indices_l,
                             constraint='wall',
                             bounds=[1.0698, 1.0698]  #1.0698 indicate 61.3°, choose your custom angle upper and lower limits
                             )
print(cnstr.scoring())

# vina scoring function
sf1 = VinaSF(receptor, ligand)
vs = sf1.scoring()
print("Vina Score ", vs)

# combined scoring function
sf = HybridSF(receptor, ligand, scorers=[sf1, cnstr], weights=[0.5, 0.5])
vs = sf.scoring()
print("HybridSF Score ", vs)

So far, we have successfully defined the angle constraint and used it. Let’s develop your own custom constraints!