Add CB atom in sidechain reconstruction if necessary

modelling/pymod/_reconstruct_sidechains.py

@@ -396,6 +396,26 @@ def _GetDisulfidBridges(frame_residues, keep_sidechains, cystein_indices,
     return disulfid_indices, disulfid_rotamers
 
 
+def _AddCB(prot):
+    '''Checks if every residue has a CB atom. Constructs them if necessary.'''
+    edi = prot.handle.EditXCS(mol.BUFFERED_EDIT)
+    for r in prot.residues:
+        olc = r.one_letter_code
+        if olc not in ['G', 'g']:
+            cb = r.handle.FindAtom('CB')
+            if not cb.IsValid():
+                n = r.handle.FindAtom('N')
+                ca = r.handle.FindAtom('CA')
+                c = r.handle.FindAtom('C')
+                if n.IsValid() and ca.IsValid() and c.IsValid():
+                    cb_pos = core.ConstructCBetaPos(n.GetPos(), ca.GetPos(), 
+                                                    c.GetPos())
+                    cb = edi.InsertAtom(r.handle, "CB", cb_pos, "C");
+                    # prot is a view on the underlying handle, so we also have
+                    # to add it to the view
+                    r.AddAtom(cb)
+
+
 ###############################################################################
 
 def ReconstructSidechains(ent, keep_sidechains=False, build_disulfids=True,
@@ -484,6 +504,9 @@ def ReconstructSidechains(ent, keep_sidechains=False, build_disulfids=True,
 
     # take out ligand chain and any non-peptides
     prot = ent.Select("peptide=true and cname!='_'")
+
+    # make sure that we have all CB atoms
+    _AddCB(prot)
     
     # parse residues (all lists of length len(prot.residues))
     rotamer_ids = _GetRotamerIDs(prot.residues)