From 557477bee3d4f542f37355b53f68c9b4281bdacc Mon Sep 17 00:00:00 2001
From: Gabriel Studer <gabriel.studer@unibas.ch>
Date: Fri, 16 Oct 2020 15:39:44 +0200
Subject: [PATCH] Add CB atom in sidechain reconstruction if necessary

---
 modelling/pymod/_reconstruct_sidechains.py | 23 ++++++++++++++++++++++
 1 file changed, 23 insertions(+)

diff --git a/modelling/pymod/_reconstruct_sidechains.py b/modelling/pymod/_reconstruct_sidechains.py
index 73563312..47228fb0 100644
--- a/modelling/pymod/_reconstruct_sidechains.py
+++ b/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)
-- 
GitLab