Allow to optimize subrotamers when FRM rotamer model is used in Reconstruct

sidechain/pymod/_reconstruct_sidechains.py

@@ -348,14 +348,43 @@ def _GetDisulfidBridges(frame_residues, cystein_indices, res_list, rotamer_libra
         disulfid_rotamers.append(cystein_rot_groups[a[1]][b[1]])
 
     return disulfid_indices, disulfid_rotamers
-        
+
+
+def RefineFRMRotamerGroups(solution, frm_rotamer_groups,
+                           graph_max_complexity=100000000 , 
+                           graph_initial_epsilon=0.02):
+    
+    rrm_rotamer_groups = list()
+    for i, rg in enumerate(frm_rotamer_groups):
+        frm_rotamer = rg[solution[i]]
+        rrm_rotamers = list()
+        original_active_subrotamer = frm_rotamer.GetActiveSubrotamer()
+        for j in range(frm_rotamer.GetNumSubrotamers()):
+            frm_rotamer.SetActiveSubrotamer(j)
+            new_rrm_rotamer = frm_rotamer.ToRRMRotamer()
+            if(j != original_active_subrotamer):
+                new_rrm_rotamer.SetInternalEnergy(0.0)
+            else:
+                new_rrm_rotamer.SetInternalEnergy(-0.5)
+            rrm_rotamers.append(new_rrm_rotamer)
+        new_rotamer_group = sidechain.RRMRotamerGroup(rrm_rotamers, 0)
+        rrm_rotamer_groups.append(new_rotamer_group)
+    graph = sidechain.RotamerGraph.CreateFromRRMList(rrm_rotamer_groups)
+
+    solution = graph.TreeSolve(max_complexity=graph_max_complexity,
+                              initial_epsilon=graph_initial_epsilon)[0]
+
+    return (solution, rrm_rotamer_groups)
+
+
+
 
 ###############################################################################
 
 def Reconstruct(ent, keep_sidechains=False, build_disulfids=True,
                 rotamer_model="frm", consider_ligands=True, 
-                rotamer_library=None, graph_max_complexity=100000000,
-                graph_intial_epsilon=0.02):
+                rotamer_library=None, optimize_subrotamers = False,
+                graph_max_complexity=100000000, graph_initial_epsilon=0.02):
     '''Reconstruct sidechains for the given structure.
 
     :param ent: Structure for sidechain reconstruction. Note, that the sidechain
@@ -385,6 +414,20 @@ def Reconstruct(ent, keep_sidechains=False, build_disulfids=True,
                             library.
     :type rotamer_library:  :class:`BBDepRotamerLib` / :class:`RotamerLib`
 
+    :param optimize_subrotamers: Only considered when **rotamer_model** 
+                                 is "frm".
+                                 If set to True, the FRM solution undergoes 
+                                 some postprocessing. The final FRM rotamers get 
+                                 turned into :class:`RRMRotamerGroup` objects 
+                                 and fed into a second run of graph solving to 
+                                 find the optimal subrotamers from the FRM 
+                                 model. This mainly improves the reconstruction 
+                                 performance of bulky sidechains such as 
+                                 PHE/TYR/TRP. Internal energies of the 
+                                 :class:`RRMRotamer` objects are set to 0.0,
+                                 -0.5 if they represent the active subrotamer 
+                                 in the :class:`FRMRotamer`.
+
     :param graph_max_complexity: Max. complexity for
                                  :meth:`RotamerGraph.TreeSolve`.
     :type graph_max_complexity:  :class:`int`
@@ -474,7 +517,13 @@ def Reconstruct(ent, keep_sidechains=False, build_disulfids=True,
         graph = sidechain.RotamerGraph.CreateFromRRMList(rotamer_groups)
 
     solution = graph.TreeSolve(max_complexity=graph_max_complexity,
-                               initial_epsilon=graph_intial_epsilon)[0]
+                               initial_epsilon=graph_initial_epsilon)[0]
+
+    if use_frm and optimize_subrotamers:
+        solution, rotamer_groups = RefineFRMRotamerGroups(solution, 
+                                                          rotamer_groups,
+                                                          graph_max_complexity,
+                                                          graph_initial_epsilon)
 
     # update structure
     for i, rot_group, sol in zip(residues_with_rotamer_group, rotamer_groups,