Return receptor and ligand interfaces from cpydock scoring function

Author: mvidalpinol

lightdock/scoring/cpydock/driver.py

@@ -119,13 +119,13 @@ def __call__(self, receptor, receptor_coordinates, ligand, ligand_coordinates):
         """Computes the pyDock scoring energy using receptor and ligand which are
         instances of DockinModel.
         """
-        elec, vdw, solv_rec, solv_lig = cpydock.calculate_energy(receptor_coordinates, ligand_coordinates,
-                                                                 receptor.charges, ligand.charges,
-                                                                 receptor.vdw_energy, ligand.vdw_energy,
-                                                                 receptor.vdw_radii, ligand.vdw_radii,
-                                                                 receptor.hydrogens, ligand.hydrogens,
-                                                                 receptor.sasa, ligand.sasa,
-                                                                 receptor.des_energy, ligand.des_energy)
+        elec, vdw, solv_rec, solv_lig, interface_receptor, interface_ligand = cpydock.calculate_energy(receptor_coordinates, ligand_coordinates,
+                                                                              receptor.charges, ligand.charges,
+                                                                              receptor.vdw_energy, ligand.vdw_energy,
+                                                                              receptor.vdw_radii, ligand.vdw_radii,
+                                                                              receptor.hydrogens, ligand.hydrogens,
+                                                                              receptor.sasa, ligand.sasa,
+                                                                              receptor.des_energy, ligand.des_energy)
         solv = -1*(solv_rec + solv_lig)
         return (elec + parameters.scoring_vdw_weight * vdw + solv)*-1.
 

lightdock/scoring/cpydock/energy/c/cpydock.c

@@ -28,28 +28,27 @@ static PyObject * cpydock_calculate_energy(PyObject *self, PyObject *args) {
     PyArrayObject *rec_charges, *lig_charges, *rec_vdw, *lig_vdw, *rec_vdw_radii, *lig_vdw_radii = NULL;
     PyArrayObject *rec_hydrogens, *lig_hydrogens, *rec_asa, *lig_asa, *rec_des_energy, *lig_des_energy = NULL;
     double atom_elec, total_elec, total_vdw, total_solvation_rec, total_solvation_lig, vdw_energy, vdw_radius, p6, k, solv_rec, solv_lig;
-    unsigned int rec_len, lig_len, i, j;
-    double **rec_array, **lig_array, x, y, z, distance2;
+    unsigned int rec_len, lig_len, i, j, interface_len, *interface_receptor, *interface_ligand;
+    double **rec_array, **lig_array, x, y, z, distance2, interface_cutoff;
     npy_intp dims[2];
-    PyArray_Descr *descr;
     double *rec_c_charges, *lig_c_charges, *rec_c_vdw, *lig_c_vdw, *rec_c_vdw_radii, *lig_c_vdw_radii = NULL;
     double *rec_c_asa, *lig_c_asa, *rec_c_des_energy, *lig_c_des_energy = NULL;
     unsigned int *rec_c_hydrogens, *lig_c_hydrogens = NULL;
     double *min_rec_distance, *min_lig_distance = NULL;
-    PyObject *result = PyTuple_New(4);
+    PyObject *result = PyTuple_New(6);
 
     total_elec = 0.0;
     atom_elec = 0.0;
     total_vdw = 0.0;
     total_solvation_rec = 0.0;
     total_solvation_lig = 0.0;
+    interface_cutoff = 3.9;
+    interface_len = 0;
 
-    if (PyArg_ParseTuple(args, "OOOOOOOOOOOOOO",
+    if (PyArg_ParseTuple(args, "OOOOOOOOOOOOOO|d",
             &receptor_coordinates, &ligand_coordinates, &rec_charges, &lig_charges,
             &rec_vdw, &lig_vdw, &rec_vdw_radii, &lig_vdw_radii, &rec_hydrogens, &lig_hydrogens,
-            &rec_asa, &lig_asa, &rec_des_energy, &lig_des_energy)) {
-
-        descr = PyArray_DescrFromType(NPY_DOUBLE);
+            &rec_asa, &lig_asa, &rec_des_energy, &lig_des_energy, interface_cutoff)) {
 
         tmp0 = PyObject_GetAttrString(receptor_coordinates, "coordinates");
         tmp1 = PyObject_GetAttrString(ligand_coordinates, "coordinates");
@@ -59,10 +58,10 @@ static PyObject * cpydock_calculate_energy(PyObject *self, PyObject *args) {
 
         dims[1] = 3;
         dims[0] = rec_len;
-        PyArray_AsCArray((PyObject **)&tmp0, (void **)&rec_array, dims, 2, descr);
+        PyArray_AsCArray((PyObject **)&tmp0, (void **)&rec_array, dims, 2, PyArray_DescrFromType(NPY_DOUBLE));
 
         dims[0] = lig_len;
-        PyArray_AsCArray((PyObject **)&tmp1, (void **)&lig_array, dims, 2, descr);
+        PyArray_AsCArray((PyObject **)&tmp1, (void **)&lig_array, dims, 2, PyArray_DescrFromType(NPY_DOUBLE));
 
         // Get pointers to the Python array structures
         rec_c_charges = PyArray_GETPTR1(rec_charges, 0);
@@ -81,6 +80,8 @@ static PyObject * cpydock_calculate_energy(PyObject *self, PyObject *args) {
         // Structures to store the atom at minimal distance of a given atom
         min_rec_distance = malloc(rec_len*sizeof(double));
         min_lig_distance = malloc(lig_len*sizeof(double));
+        interface_receptor = malloc(lig_len*sizeof(unsigned int));
+        interface_ligand  = malloc(lig_len*sizeof(unsigned int));
         for (i = 0; i < rec_len; i++) min_rec_distance[i] = HUGE_DISTANCE;
         for (j = 0; j < lig_len; j++) min_lig_distance[j] = HUGE_DISTANCE;
 
@@ -109,14 +110,25 @@ static PyObject * cpydock_calculate_energy(PyObject *self, PyObject *args) {
                 }
 
                 // Van der Waals energy
-                if (distance2 <= VDW_DIST_CUTOFF2){
+                if (distance2 <= VDW_DIST_CUTOFF2) {
                     vdw_energy = sqrt(rec_c_vdw[i] * lig_c_vdw[j]);
                     vdw_radius = rec_c_vdw_radii[i] + lig_c_vdw_radii[j];
                     p6 = pow(vdw_radius, 6) / pow(distance2, 3);
                     k = vdw_energy * (p6*p6 - 2.0 * p6);
                     if (k > VDW_CUTOFF) k = VDW_CUTOFF;
                     total_vdw += k;
                 }
+
+                if (sqrt(distance2) <= interface_cutoff) {
+                   interface_receptor[interface_len] = i;
+                   interface_ligand[interface_len++] = j;
+                }
+
+            }
+
+            if (!(interface_len%lig_len)) {
+                interface_receptor = realloc(interface_receptor, (interface_len/lig_len + 1)*lig_len*sizeof(unsigned int));
+                interface_ligand = realloc(interface_ligand, (interface_len/lig_len + 1)*lig_len*sizeof(unsigned int));
             }
         }
         // Convert total electrostatics to Kcal/mol:
@@ -146,18 +158,23 @@ static PyObject * cpydock_calculate_energy(PyObject *self, PyObject *args) {
         // Free structures
         PyArray_Free(tmp0, rec_array);
         PyArray_Free(tmp1, lig_array);
-        Py_DECREF(descr);
         Py_DECREF(tmp0);
         Py_DECREF(tmp1);
         free(min_rec_distance);
         free(min_lig_distance);
     }
 
+    interface_receptor = realloc(interface_receptor, interface_len*sizeof(unsigned int));
+    interface_ligand = realloc(interface_ligand, interface_len*sizeof(unsigned int));
+    dims[0] = interface_len;
+
     // Return a tuple with the following values for calculated energies:
     PyTuple_SetItem(result, 0, PyFloat_FromDouble(total_elec));
     PyTuple_SetItem(result, 1, PyFloat_FromDouble(total_vdw));
     PyTuple_SetItem(result, 2, PyFloat_FromDouble(total_solvation_rec));
     PyTuple_SetItem(result, 3, PyFloat_FromDouble(total_solvation_lig));
+    PyTuple_SetItem(result, 4, PyArray_SimpleNewFromData(1, dims, NPY_UINT, interface_receptor));
+    PyTuple_SetItem(result, 5, PyArray_SimpleNewFromData(1, dims, NPY_UINT, interface_ligand));
     return result;
 }