Probe2 fixes (#1032)

* probe2 v4.3.0 adds False for stop_for_unknown and flip_symmetric_amino_acids flags

* Probe2 v4.4.0 handles element name capitalization consistently, supporting two-letter element names properly

* More verbose reporting of failure to find atom information

* probe2 v4.3.0 adds False for stop_for_unknown and flip_symmetric_amino_acids flags

* Probe2 v4.4.0 handles element name capitalization consistently, supporting two-letter element names properly

* More verbose reporting of failure to find atom information

* Version 4.5.0 properly right-adjusts the residue names in a 3-character field to handle DNA and RNA properly

* Version 4.6.0 properly reports atom classes rather than elements in raw output

* Version 4.7.0 of Probe2 computes the spike lengths and gaps correctly in raw output

mmtbx/probe/Helpers.py

@@ -465,14 +465,14 @@ def getExtraAtomInfo(model, bondedNeighborLists, useNeutronDistances = False, pr
               else:
                 fullName = (chain.id + ' ' + a.parent().resname.strip() + ' ' +
                   str(a.parent().parent().resseq_as_int()) + ' ' + a.name.strip())
-                raise Sorry("Could not find atom info for "+fullName+
-                  "(perhaps interpretation was not run on the model?)\n")
+                raise Sorry("Unrecognized specific H bond type for "+fullName+", got "+hb_type+"\n")
 
             except Exception as e:
               fullName = (chain.id + ' ' + a.parent().resname.strip() + ' ' +
                 str(a.parent().parent().resseq_as_int()) + ' ' + a.name.strip())
               raise Sorry("Could not find atom info for "+fullName+
-                " (perhaps interpretation was not run on the model?)\n")
+                " (perhaps interpretation was not run on the model?)\n"+
+                "  Exception: "+str(e))
 
   return getExtraAtomInfoReturn(extras, warnings)
 

mmtbx/programs/probe2.py

@@ -29,7 +29,7 @@
 # @todo See if we can remove the shift and box once reduce_hydrogen is complete
 from cctbx.maptbx.box import shift_and_box_model
 
-version = "4.2.0"
+version = "4.7.0"
 
 master_phil_str = '''
 profile = False
@@ -712,10 +712,14 @@ def _atom_class_for(self, a):
       :param a: Atom whose class is to be specified
       :return: If our parameters have been set to color and sort by NA base,
       then it returns the appropriate base name.  Otherwise, it returns the
-      element of the atom.
+      element of the atom with any second letter in the element name lower-case
+      to match the values in the _allAtomClasses list.
     '''
     if not self.params.output.color_by_na_base:
-      return a.element
+      val = a.element
+      if len(val) > 1:
+        val = val[0] + val[1:].lower()
+      return val
     else:
       resName = a.parent().resname
       cl = common_residue_names_get_class(name = resName)
@@ -1134,7 +1138,7 @@ def _writeRawOutput(self, groupName, masterName):
           chainID = a.parent().parent().parent().id
           iCode = a.parent().parent().icode
           alt = a.parent().altloc
-          ret += "{:>2s}{:>4s}{}{} {}{:1s}:".format(chainID, resID, iCode, resName, a.name, alt)
+          ret += "{:>2s}{:>4s}{}{:>3s} {}{:1s}:".format(chainID, resID, iCode, resName, a.name, alt)
 
           # Describe the target atom, if it exists
           t = node.target
@@ -1146,12 +1150,12 @@ def _writeRawOutput(self, groupName, masterName):
             chainID = t.parent().parent().parent().id
             iCode = t.parent().parent().icode
             alt = t.parent().altloc
-            ret += "{:>2s}{:>4s}{}{} {:<3s}{:1s}:".format(chainID, resID, iCode, resName, t.name, alt)
+            ret += "{:>2s}{:>4s}{}{:>3s} {:<3s}{:1s}:".format(chainID, resID, iCode, resName, t.name, alt)
 
             r1 = self._extraAtomInfo.getMappingFor(a).vdwRadius
             r2 = self._extraAtomInfo.getMappingFor(t).vdwRadius
-            sl = (Helpers.rvec3(a.xyz)-Helpers.rvec3(t.xyz)).length()
-            gap = sl - (r1 + r2)
+            sl = (node.loc-node.spike).length()
+            gap = (Helpers.rvec3(a.xyz)-Helpers.rvec3(t.xyz)).length() - (r1 + r2)
             dtgp = node.gap
             score = 0.0
 
@@ -1174,12 +1178,12 @@ def _writeRawOutput(self, groupName, masterName):
               node.spike[0], node.spike[1], node.spike[2], sl, score/density)
 
           try:
-            tName = t.element
+            tName = self._atomClasses[t]
             tBVal = "{:.2f}".format(t.b)
           except Exception:
             tName = ""
             tBVal = ""
-          ret += ":{}:{}:{:.3f}:{:.3f}:{:.3f}".format(a.element, tName,
+          ret += ":{}:{}:{:.3f}:{:.3f}:{:.3f}".format(self._atomClasses[a], tName,
             node.loc[0], node.loc[1], node.loc[2])
 
           ret += ":{:.2f}:{}\n".format(a.b, tBVal)
@@ -1314,7 +1318,7 @@ def _writeOutput(self, groupName, masterName):
           else:
             ptmast = " '{}' ".format(node.ptmaster)
 
-          pointid = "{}{:1s}{} {:>3d} {:1s}{}".format(a.name, a.parent().altloc, a.parent().resname,
+          pointid = "{}{:1s}{:>3s} {:>3d} {:1s}{}".format(a.name, a.parent().altloc, a.parent().resname,
             a.parent().parent().resseq_as_int(), a.parent().parent().icode,
             a.parent().parent().parent().id)
           if pointid != lastPointID:
@@ -1845,11 +1849,17 @@ def run(self):
     # Get the bonding information we'll need to exclude our bonded neighbors.
     allAtoms = self.model.get_atoms()
     make_sub_header('Compute neighbor lists', out=self.logger)
+
+    self.model.set_stop_for_unknowns(False)
     p = mmtbx.model.manager.get_default_pdb_interpretation_params()
     p.pdb_interpretation.use_neutron_distances = self.params.use_neutron_distances
     p.pdb_interpretation.allow_polymer_cross_special_position=True
     p.pdb_interpretation.clash_guard.nonbonded_distance_threshold=None
     p.pdb_interpretation.proceed_with_excessive_length_bonds=True
+    # We need to turn this on because without it the interpretation is
+    # renaming atoms to be more correct.  Unfortunately, this causes the
+    # dot names to no longer match the input file.
+    p.pdb_interpretation.flip_symmetric_amino_acids=False
     try:
       self.model.process(make_restraints=True, pdb_interpretation_params=p, logger=self.logger) # make restraints
       geometry = self.model.get_restraints_manager().geometry