fix: correct usage of nonbond energy for close contacts

deeprankcore/features/contact.py

@@ -1,4 +1,4 @@
-from typing import List
+from typing import List, Tuple
 import logging
 import warnings
 import numpy as np
@@ -8,48 +8,62 @@
 from deeprankcore.molstruct.pair import ResidueContact, AtomicContact
 from deeprankcore.domain import edgestorage as Efeat
 from deeprankcore.utils.parsing import atomic_forcefield
+import numpy.typing as npt
 
 _log = logging.getLogger(__name__)
 
-MAX_COVALENT_DISTANCE = 2.1
+# cutoff distances for 1-3 and 1-4 pairing. See issue: https://github.com/DeepRank/deeprank-core/issues/357#issuecomment-1461813723
+covalent_cutoff = 2.1
+cutoff_13 = 3.6
+cutoff_14 = 4.2
 
-def _get_coulomb_potentials(atoms: List[Atom], distances: np.ndarray) -> np.ndarray:
-    """Calculate Coulomb potentials between between all Atoms in atom.
+
+def _get_nonbonded_energy(atoms: List[Atom], distances: npt.NDArray[np.float64]) -> Tuple [npt.NDArray[np.float64], npt.NDArray[np.float64]]:
+    """Calculates all pairwise electrostatic (Coulomb) and Van der Waals (Lennard Jones) potential energies between all atoms in the structure.
 
     Warning: there's no distance cutoff here. The radius of influence is assumed to infinite.
     However, the potential tends to 0 at large distance.
+
+    Args:
+        atoms (List[Atom]): list of all atoms in the structure 
+        distances (npt.NDArray[np.float64]): matrix of pairwise distances between all atoms in the structure 
+            in the format that is the output of scipy.spatial's distance_matrix (i.e. a diagonally symmetric matrix)
+
+    Returns:
+        Tuple [npt.NDArray[np.float64], npt.NDArray[np.float64]]: matrices in same format as `distances` containing
+            all pairwise electrostatic potential energies and all pairwise Van der Waals potential energies
     """
 
+    # ELECTROSTATIC POTENTIAL
     EPSILON0 = 1.0
     COULOMB_CONSTANT = 332.0636
     charges = [atomic_forcefield.get_charge(atom) for atom in atoms]
-    coulomb_potentials = np.expand_dims(charges, axis=1) * np.expand_dims(charges, axis=0) * COULOMB_CONSTANT / (EPSILON0 * distances)
-    return coulomb_potentials
+    electrostatic_energy = np.expand_dims(charges, axis=1) * np.expand_dims(charges, axis=0) * COULOMB_CONSTANT / (EPSILON0 * distances)
+    # remove for close contacts
+    electrostatic_energy[distances < cutoff_13] = 0
 
 
-def _get_lennard_jones_potentials(atoms: List[Atom], distances: np.ndarray) -> np.ndarray:
-    """Calculate Lennard-Jones potentials between all Atoms in atom.
-
-    Warning: there's no distance cutoff here. The radius of influence is assumed to infinite.
-    However, the potential tends to 0 at large distance.
-    """
-
-    # calculate intra potentials
-    sigmas = [atomic_forcefield.get_vanderwaals_parameters(atom).intra_sigma for atom in atoms]
-    epsilons = [atomic_forcefield.get_vanderwaals_parameters(atom).intra_epsilon for atom in atoms]
+    # VAN DER WAALS POTENTIAL
+    # calculate main vdw energies
+    sigmas = [atomic_forcefield.get_vanderwaals_parameters(atom).sigma_main for atom in atoms]
+    epsilons = [atomic_forcefield.get_vanderwaals_parameters(atom).epsilon_main for atom in atoms]
     mean_sigmas = 0.5 * np.add.outer(sigmas,sigmas)
     geomean_eps = np.sqrt(np.multiply.outer(epsilons,epsilons))     # sqrt(eps1*eps2)
-    intra_potentials = 4.0 * geomean_eps * ((mean_sigmas / distances) ** 12 - (mean_sigmas / distances) ** 6)
-    
-    # calculate inter potentials
-    sigmas = [atomic_forcefield.get_vanderwaals_parameters(atom).inter_sigma for atom in atoms]
-    epsilons = [atomic_forcefield.get_vanderwaals_parameters(atom).inter_epsilon for atom in atoms]
+    vdw_energy = 4.0 * geomean_eps * ((mean_sigmas / distances) ** 12 - (mean_sigmas / distances) ** 6)
+
+    # calculate energies for 1-4 pairs
+    sigmas = [atomic_forcefield.get_vanderwaals_parameters(atom).sigma_14 for atom in atoms]
+    epsilons = [atomic_forcefield.get_vanderwaals_parameters(atom).epsilon_14 for atom in atoms]
     mean_sigmas = 0.5 * np.add.outer(sigmas,sigmas)
     geomean_eps = np.sqrt(np.multiply.outer(epsilons,epsilons))     # sqrt(eps1*eps2)
-    inter_potentials = 4.0 * geomean_eps * ((mean_sigmas / distances) ** 12 - (mean_sigmas / distances) ** 6)
+    energy_14pairs = 4.0 * geomean_eps * ((mean_sigmas / distances) ** 12 - (mean_sigmas / distances) ** 6)
 
-    lennard_jones_potentials = {'intra': intra_potentials, 'inter': inter_potentials}
-    return lennard_jones_potentials
+    # adjust vdw energy for close contacts
+    vdw_energy[distances < cutoff_14] = energy_14pairs[distances < cutoff_14]
+    vdw_energy[distances < cutoff_13] = 0
+
+
+    return electrostatic_energy, vdw_energy
 
 
 def add_features(pdb_path: str, graph: Graph, *args, **kwargs): # pylint: disable=too-many-locals, unused-argument
@@ -76,8 +90,7 @@ def add_features(pdb_path: str, graph: Graph, *args, **kwargs): # pylint: disabl
         warnings.simplefilter("ignore")
         positions = [atom.position for atom in all_atoms]
         interatomic_distances = distance_matrix(positions, positions)
-        interatomic_electrostatic_potentials = _get_coulomb_potentials(all_atoms, interatomic_distances)
-        interatomic_vanderwaals_potentials = _get_lennard_jones_potentials(all_atoms, interatomic_distances)
+        interatomic_electrostatic_energy, interatomic_vanderwaals_energy = _get_nonbonded_energy(all_atoms, interatomic_distances)
 
     # assign features
     if isinstance(graph.edges[0].id, AtomicContact):
@@ -90,12 +103,9 @@ def add_features(pdb_path: str, graph: Graph, *args, **kwargs): # pylint: disabl
             edge.features[Efeat.SAMERES] = float( contact.atom1.residue == contact.atom2.residue)  # 1.0 for True; 0.0 for False
             edge.features[Efeat.SAMECHAIN] = float( contact.atom1.residue.chain == contact.atom1.residue.chain )  # 1.0 for True; 0.0 for False
             edge.features[Efeat.DISTANCE] = interatomic_distances[atom1_index, atom2_index]
-            edge.features[Efeat.COVALENT] = float( edge.features[Efeat.DISTANCE] < MAX_COVALENT_DISTANCE )  # 1.0 for True; 0.0 for False
-            edge.features[Efeat.ELECTROSTATIC] = interatomic_electrostatic_potentials[atom1_index, atom2_index]
-            if edge.features[Efeat.SAMERES]:
-                edge.features[Efeat.VANDERWAALS] = interatomic_vanderwaals_potentials['intra'][atom1_index, atom2_index]
-            else:
-                edge.features[Efeat.VANDERWAALS] = interatomic_vanderwaals_potentials['inter'][atom1_index, atom2_index]
+            edge.features[Efeat.COVALENT] = float( edge.features[Efeat.DISTANCE] < covalent_cutoff )  # 1.0 for True; 0.0 for False
+            edge.features[Efeat.ELECTROSTATIC] = interatomic_electrostatic_energy[atom1_index, atom2_index]
+            edge.features[Efeat.VANDERWAALS] = interatomic_vanderwaals_energy[atom1_index, atom2_index]
 
     elif isinstance(contact, ResidueContact):
         for edge in graph.edges:        
@@ -106,6 +116,6 @@ def add_features(pdb_path: str, graph: Graph, *args, **kwargs): # pylint: disabl
             ## set features
             edge.features[Efeat.SAMECHAIN] = float( contact.residue1.chain == contact.residue2.chain )  # 1.0 for True; 0.0 for False
             edge.features[Efeat.DISTANCE] = np.min([[interatomic_distances[a1, a2] for a1 in atom1_indices] for a2 in atom2_indices])
-            edge.features[Efeat.COVALENT] = float( edge.features[Efeat.DISTANCE] < MAX_COVALENT_DISTANCE )  # 1.0 for True; 0.0 for False
-            edge.features[Efeat.ELECTROSTATIC] = np.sum([[interatomic_electrostatic_potentials[a1, a2] for a1 in atom1_indices] for a2 in atom2_indices])
-            edge.features[Efeat.VANDERWAALS] = np.sum([[interatomic_vanderwaals_potentials['inter'][a1, a2] for a1 in atom1_indices] for a2 in atom2_indices])
+            edge.features[Efeat.COVALENT] = float( edge.features[Efeat.DISTANCE] < covalent_cutoff )  # 1.0 for True; 0.0 for False
+            edge.features[Efeat.ELECTROSTATIC] = np.sum([[interatomic_electrostatic_energy[a1, a2] for a1 in atom1_indices] for a2 in atom2_indices])
+            edge.features[Efeat.VANDERWAALS] = np.sum([[interatomic_vanderwaals_energy[a1, a2] for a1 in atom1_indices] for a2 in atom2_indices])

deeprankcore/utils/parsing/vdwparam.py

@@ -1,17 +1,17 @@
 class VanderwaalsParam:
     def __init__(
             self,
-            inter_epsilon: float,
-            inter_sigma: float,
-            intra_epsilon: float,
-            intra_sigma: float):
-        self.inter_epsilon = inter_epsilon
-        self.inter_sigma = inter_sigma
-        self.intra_epsilon = intra_epsilon
-        self.intra_sigma = intra_sigma
+            epsilon_main: float,
+            sigma_main: float,
+            epsilon_14: float,
+            sigma_14: float):
+        self.epsilon_main = epsilon_main
+        self.sigma_main = sigma_main
+        self.epsilon_14 = epsilon_14
+        self.sigma_14 = sigma_14
 
     def __str__(self) -> str:
-        return f"{self.inter_epsilon}, {self.inter_sigma}, {self.intra_epsilon}, {self.intra_sigma}"
+        return f"{self.epsilon_main}, {self.sigma_main}, {self.epsilon_14}, {self.sigma_14}"
 
 
 class ParamParser:
@@ -26,17 +26,17 @@ def parse(file_):
                 (
                     _,
                     type_,
-                    inter_epsilon,
-                    inter_sigma,
-                    intra_epsilon,
-                    intra_sigma,
+                    epsilon_main,
+                    sigma_main,
+                    epsilon_14,
+                    sigma_14,
                 ) = line.split()
 
                 result[type_] = VanderwaalsParam(
-                    float(inter_epsilon),
-                    float(inter_sigma),
-                    float(intra_epsilon),
-                    float(intra_sigma),
+                    float(epsilon_main),
+                    float(sigma_main),
+                    float(epsilon_14),
+                    float(sigma_14),
                 )
             elif len(line.strip()) == 0:
                 continue

tests/features/test_contact.py

@@ -32,12 +32,12 @@ def _wrap_in_graph(edge: Edge):
     return g
 
 
-def _get_contact(residue_num1: int, atom_name1: str, residue_num2: int, atom_name2: str, residue_level: bool = False):
-    pdb_path = "tests/data/pdb/101M/101M.pdb"
+def _get_contact(pdb_id: str, residue_num1: int, atom_name1: str, residue_num2: int, atom_name2: str, residue_level: bool = False) -> Edge:
+    pdb_path = f"tests/data/pdb/101M/{pdb_id}.pdb"
 
     pdb = pdb2sql(pdb_path)
     try:
-        structure = get_structure(pdb, "101m")
+        structure = get_structure(pdb, pdb_id)
     finally:
         pdb._close() # pylint: disable=protected-access
 
@@ -68,38 +68,61 @@ def _get_contact(residue_num1: int, atom_name1: str, residue_num2: int, atom_nam
     return edge_obj
 
 
-def test_vanderwaals_positive():
-    """MET 0: N - CA, very close. Should have positive vanderwaals energy.
+def test_covalent_pair():
+    """MET 0: N - CA, covalent pair (at 1.49 A distance). Should have 0 vanderwaals and electrostatic energies.
     """
 
-    edge_close=_get_contact(0,"N",0,"CA")
-    assert edge_close.features[Efeat.VANDERWAALS] > 0.0
+    edge_covalent = _get_contact('101M', 0, "N", 0, "CA")
+    assert edge_covalent.features[Efeat.VANDERWAALS] == 0.0, 'nonzero vdw energy for covalent pair'
+    assert edge_covalent.features[Efeat.ELECTROSTATIC] == 0.0, 'nonzero electrostatic energy for covalent pair'
+    assert edge_covalent.features[Efeat.COVALENT] == 1.0, 'covalent pair not recognized as covalent'
+
+
+def test_13_pair():
+    """MET 0: N - CB, 1-3 pair (at 2.47 A distance). Should have 0 vanderwaals and electrostatic energies.
+    """
+
+    edge_13 = _get_contact('101M', 0, "N", 0, "CB")
+    assert edge_13.features[Efeat.VANDERWAALS] == 0.0, 'nonzero vdw energy for 1-3 pair'
+    assert edge_13.features[Efeat.ELECTROSTATIC] == 0.0, 'nonzero electrostatic energy for 1-3 pair'
+    assert edge_13.features[Efeat.COVALENT] == 0.0, '1-3 pair recognized as covalent'
+
+
+def test_14_pair():
+    """MET 0: N - CG, 1-4 pair (at 4.12 A distance). Should have non-zero electrostatic energy and small non-zero vdw energy.
+    """
+
+    edge_14 = _get_contact('101M', 0, "CA", 0, "SD")
+    assert edge_14.features[Efeat.VANDERWAALS] != 0.0, '1-4 pair with 0 vdw energy'
+    assert abs(edge_14.features[Efeat.VANDERWAALS]) < 0.1, '1-4 pair with large vdw energy'
+    assert edge_14.features[Efeat.ELECTROSTATIC] != 0.0, '1-4 pair with 0 electrostatic'
+    assert edge_14.features[Efeat.COVALENT] == 0.0, '1-4 pair recognized as covalent'
 
 
 def test_vanderwaals_negative():
-    """MET 0 N - ASP 27 CB, very far. Should have negative vanderwaals energy.
+    """MET 0 N - ASP 27 CB, very far (29.54 A). Should have negative vanderwaals energy.
     """
 
-    edge_far = _get_contact(0,"N",27,"CB")
+    edge_far = _get_contact('101M', 0, "N", 27, "CB")
     assert edge_far.features[Efeat.VANDERWAALS] < 0.0
 
 
 def test_vanderwaals_morenegative():
-    """MET 0 N - PHE 138 CG, intermediate distance. Should have more negative vanderwaals energy than the far interaction.
+    """MET 0 N - PHE 138 CG, intermediate distance (12.69 A). Should have more negative vanderwaals energy than the far interaction.
     """
 
-    edge_intermediate=_get_contact(0,"N",138,"CG")
-    edge_far = _get_contact(0,"N",27,"CB")
+    edge_intermediate = _get_contact('101M', 0, "N", 138, "CG")
+    edge_far = _get_contact('101M', 0, "N", 27, "CB")
     assert edge_intermediate.features[Efeat.VANDERWAALS] < edge_far.features[Efeat.VANDERWAALS]
 
 
 def test_edge_distance():
     """Check the edge distances.
     """
 
-    edge_close=_get_contact(0,"N",0,"CA")
-    edge_intermediate=_get_contact(0,"N",138,"CG")
-    edge_far = _get_contact(0,"N",27,"CB")
+    edge_close = _get_contact('101M', 0, "N", 0, "CA")
+    edge_intermediate = _get_contact('101M', 0, "N", 138, "CG")
+    edge_far = _get_contact('101M', 0, "N", 27, "CB")
 
     assert (
         edge_close.features[Efeat.DISTANCE]
@@ -112,19 +135,19 @@ def test_edge_distance():
 
 
 def test_attractive_electrostatic_close():
-    """ARG 139 CZ - GLU 136 OE2, very close. Should have attractive electrostatic energy.
+    """ARG 139 CZ - GLU 136 OE2, very close (5.60 A). Should have attractive electrostatic energy.
     """
 
-    close_attracting_edge = _get_contact(139,"CZ",136,"OE2")
+    close_attracting_edge = _get_contact('101M', 139, "CZ", 136, "OE2")
     assert close_attracting_edge.features[Efeat.ELECTROSTATIC] < 0.0
 
 
 def test_attractive_electrostatic_far():
-    """ARG 139 CZ - ASP 20 OD2, far. Should have attractive more electrostatic energy than above.
+    """ARG 139 CZ - ASP 20 OD2, far (24.26 A). Should have attractive more electrostatic energy than above.
     """
 
-    far_attracting_edge=_get_contact(139,"CZ",20,"OD2")
-    close_attracting_edge = _get_contact(139,"CZ",136,"OE2")
+    far_attracting_edge = _get_contact('101M', 139, "CZ", 20, "OD2")
+    close_attracting_edge = _get_contact('101M', 139, "CZ", 136, "OE2")
     assert (
         far_attracting_edge.features[Efeat.ELECTROSTATIC] < 0.0
     ), 'far electrostatic > 0'
@@ -135,19 +158,20 @@ def test_attractive_electrostatic_far():
 
 
 def test_repulsive_electrostatic():
-    """GLU 109 OE2 - GLU 105 OE1. Should have repulsive electrostatic energy.
+    """GLU 109 OE2 - GLU 105 OE1 (9.64 A). Should have repulsive electrostatic energy.
     """
 
-    opposing_edge=_get_contact(109,"OE2",105,"OE1")
+    opposing_edge = _get_contact('101M', 109, "OE2", 105, "OE1")
     assert opposing_edge.features[Efeat.ELECTROSTATIC] > 0.0
 
 
 def test_residue_contact():
     """Check that we can calculate features for residue contacts.
     """
 
-    res_edge = _get_contact(0, '', 1, '', residue_level = True)
+    res_edge = _get_contact('101M', 0, '', 1, '', residue_level = True)
     assert res_edge.features[Efeat.DISTANCE] > 0.0, 'distance <= 0'
     assert res_edge.features[Efeat.DISTANCE] < 1e5, 'distance > 1e5'
     assert res_edge.features[Efeat.ELECTROSTATIC] != 0.0, 'electrostatic == 0'
     assert res_edge.features[Efeat.VANDERWAALS] != 0.0, 'vanderwaals == 0'
+    assert res_edge.features[Efeat.COVALENT] == 1.0, 'neighboring residues not seen as covalent'