diff --git a/sidechain/doc/index.rst b/sidechain/doc/index.rst
index bfca8940ed6af8f0c1897b3274a679ea58eee138..98c570f303b9a79bdea34dd6a866864acb93ce05 100644
--- a/sidechain/doc/index.rst
+++ b/sidechain/doc/index.rst
@@ -20,43 +20,7 @@ function:
 
 .. literalinclude:: ../../../tests/doc/scripts/sidechain_reconstruct.py
 
-.. method:: Reconstruct(prot, keep_sidechains=False, build_disulfids, \
-                        rotamer_model="frm", consider_hbonds=True, \
-                        rotamer_library=None)
-
-  The function takes a structure and reconstructs its sidechains given the input 
-  parameters.
-
-  :param prot:          Structure for sidechain reconstruction. Note, that the
-                        sidechain reconstruction gets directly applied on the
-                        structure itself.
-
-  :param keep_sidechains: Flag, whether complete sidechains in *prot* (i.e. 
-                          containing all required atoms) should be kept rigid
-                          and directly be added to the frame.
-
-  :param build_disulfids: Flag, whether possible disulfid bonds should be 
-                          searched. If a disulfid bond is found, the two
-                          participating cysteins are fixed and added to
-                          the frame.
-
-  :param rotamer_model: Rotamer model to be used, can either be "frm" or "rrm"
-
-  :param consider_hbonds: Flag, whether hbonds should be evaluated in the energy
-                          function. If set to False, no hydrogens will be built
-                          when building rotamers and frame and the **add_polar_hydrogens**
-                          flag won't have any consequences.
-
-  :param rotamer_library: A rotamer library to extract the 
-                          rotamers from.
-
-
-  :type prot:           :class:`ost.mol.EntityHandle`
-  :type keep_sidechains: :class:`bool`
-  :type build_disulfids: :class:`bool`
-  :type rotamer_model:  :class:`str`
-  :type consider_hbonds: :class:`bool`
-  :type rotamer_library: :class:`BBDepRotamerLib` / :class:`RotamerLib`
+.. autofunction:: Reconstruct
 
 
 Sidechain Module Functionality
diff --git a/sidechain/pymod/_reconstruct_sidechains.py b/sidechain/pymod/_reconstruct_sidechains.py
index c385aac2f632bb716d769501f334b0338494067f..bb8f5410857019c9a64358643ac01d1e37b603f9 100644
--- a/sidechain/pymod/_reconstruct_sidechains.py
+++ b/sidechain/pymod/_reconstruct_sidechains.py
@@ -2,54 +2,36 @@ from . import _sidechain as sidechain
 from ost import geom
 from ost import mol
 
-def Reconstruct(ent, keep_sidechains = False, build_disulfids = True, 
-                rotamer_model = "frm", consider_hbonds = True,
-                rotamer_library = None):
-
-    if rotamer_model.lower() == "frm":
-        use_frm = True
-    elif rotamer_model.lower() == "rrm":
-        use_frm = False
-    else:
-        raise RuntimeError("Only \"rrm\" and \"frm\" allowed for rotamer_model!")
-
-
-    rotamer_settings = sidechain.RotamerSettings()
-    rotamer_settings.consider_hbonds = consider_hbonds
-    if rotamer_library == None: 
-        rotamer_library = sidechain.LoadDunbrackLib()
-    bbdep = False
-    if type(rotamer_library) is sidechain.BBDepRotamerLib:
-        bbdep = True
-    residues_with_rotamer_group = list()
-    rotamer_groups = list()
-    # ignore ligand chain and any non-peptides
-    prot = ent.Select("peptide=true and cname!='_'")
-    incomplete_sidechains = list()
-
-    #extract rotamer ids
-    rotamer_ids = [sidechain.ALA] * len(prot.residues)
-    for i,r in enumerate(prot.residues):
+###############################################################################
+# helper functions
+def _GetRotamerIDs(res_list):
+    '''Return list (length = len(res_list)) of rotamer IDs for all residues.'''
+    rotamer_ids = [sidechain.ALA] * len(res_list)
+    for i,r in enumerate(res_list):
         rot_id = sidechain.TLCToRotID(r.GetName())
         if rot_id == sidechain.XXX:
             continue # no idea what it is, so we stick with ALA 
                      # => don't model sidechain
         rotamer_ids[i] = rot_id
-
-    #extract dihedral angles
-    phi_angles = [0.0] * len(prot.residues)
-    psi_angles = [0.0] * len(prot.residues)
-  
-    for i,r in enumerate(prot.residues):
-
+    return rotamer_ids
+
+def _GetDihedrals(res_list):
+    '''Extract dihedral angles for all residues.
+    Returns phi and psi angles as 2 lists with length = len(res_list).
+    '''
+    phi_angles = [0.0] * len(res_list)
+    psi_angles = [0.0] * len(res_list)
+    for i,r in enumerate(res_list):
+        # def. fallback = helix
         phi = -1.0472 
         psi = -0.7854
-  
+
+        # try to get phi from torsion angles
         tor = r.GetPhiTorsion()
         if tor.IsValid():
             phi = tor.GetAngle()
         elif i > 0:
-            r_prev = prot.residues[i-1]
+            r_prev = res_list[i-1]
             if mol.InSequence(r_prev.handle,r.handle):
                 c_prev = r_prev.FindAtom("C")
                 n = r.FindAtom("N")
@@ -58,12 +40,13 @@ def Reconstruct(ent, keep_sidechains = False, build_disulfids = True,
                 if c_prev.IsValid() and n.IsValid() and ca.IsValid() and c.IsValid():
                     phi = geom.DihedralAngle(c_prev.GetPos(),n.GetPos(),
                                              ca.GetPos(),c.GetPos())
-      
+        
+        # try to get psi from torsion angles
         tor = r.GetPsiTorsion()
         if tor.IsValid():
             psi = tor.GetAngle()
-        elif i < (len(prot.residues) - 1):
-            r_next = prot.residues[i+1]
+        elif i < (len(res_list) - 1):
+            r_next = res_list[i+1]
             if mol.InSequence(r.handle,r_next.handle):
                 n = r.FindAtom("N")
                 ca = r.FindAtom("CA")
@@ -72,222 +55,154 @@ def Reconstruct(ent, keep_sidechains = False, build_disulfids = True,
                 if n.IsValid() and ca.IsValid() and c.IsValid() and n_next.IsValid():
                     psi = geom.DihedralAngle(n.GetPos(), ca.GetPos(),
                                              c.GetPos(), n_next.GetPos())
-
-
+        # store
         phi_angles[i] = phi
         psi_angles[i] = psi
+    return phi_angles, psi_angles
 
-
-    #set nter and cter
-    for c in prot.chains:
-        c.residues[0].SetIntProp("n_ter",1)
-        c.residues[-1].SetIntProp("c_ter",1)
-
-
-    #build up frame and cysteins
-    frame_residues = list()
-
-    #build up backbone frame
-    for i,r in enumerate(prot.residues):
+def _AddBackboneFrameResidues(frame_residues, res_list, rotamer_ids,
+                              rotamer_settings, phi_angles):
+    '''Update frame_residues (list) with BackboneFrameResidues for res_list.'''
+    for i,r in enumerate(res_list):
         try:
-            frame_residue = sidechain.ConstructBackboneFrameResidue(r.handle,
-                                                                    rotamer_ids[i],
-                                                                    i,rotamer_settings,
-                                                                    phi_angles[i],
-                                                                    r.HasProp("n_ter"),
-                                                                    r.HasProp("c_ter"))
+            frame_residue = sidechain.ConstructBackboneFrameResidue(\
+                                r.handle, rotamer_ids[i], i, rotamer_settings,
+                                phi_angles[i], r.HasProp("n_ter"), r.HasProp("c_ter"))
             frame_residues.append(frame_residue)
         except:
             continue
 
+def _AddLigandFrameResidues(frame_residues, ent_lig, rotamer_settings, offset):
+    '''Update frame_residues (list) with FrameResidues for res. in ent_lig.
+    Set offset >= number of non-ligand residues (used for residue_index).
+    '''
+    # TODO? Split peptides and non-peptides?
+    # for now: add as static carbons
+    for i,r in enumerate(ent_lig.residues):
+        try:
+            frame_residue = sidechain.ConstructFrameResidue(r.handle, offset+i)
+            frame_residues.append(frame_residue)
+        except:
+            continue
 
-    frame = None
-
-    if build_disulfids:
-
-        cystein_indices = list()
-
-        if keep_sidechains:
-            #let's add all complete sidechains except the cysteins to the frame
-            for i,r in enumerate(prot.residues):
-
-                if rotamer_ids[i] == sidechain.CYS:
-                    cystein_indices.append(i)
-                    continue
-
-                if rotamer_ids[i] == sidechain.ALA or rotamer_ids[i] == sidechain.GLY:
-                    continue #no sidechain to model
-
-                try:
-                    frame_residue = sidechain.ConstructSidechainFrameResidue(r.handle,
-                                                                             rotamer_ids[i],
-                                                                             i,rotamer_settings)
-                    frame_residues.append(frame_residue)
-                except:
-                    incomplete_sidechains.append(i)
-        else:
-            #let's add everything except cysteins to the incomplete sidechains
-            for i,r in enumerate(prot.residues):
-
-                if rotamer_ids[i] == sidechain.CYS:
-                    cystein_indices.append(i)
-                    continue
+def _AddSidechainFrameResidues(frame_residues, incomplete_sidechains,
+                               keep_sidechains, res_list, rotamer_ids,
+                               rotamer_settings, cystein_indices=None):
+    '''Update frame_residues (list) with SidechainFrameResidues for res_list,
+    incomplete_sidechains (list of indices) with sidechains to be constructed,
+    and (if given) cystein_indices (list of indices) with all CYS (appended).
+    Each residue can only end up in one of the 3 lists.
+    '''
+    if keep_sidechains:
+        # try to generate frame residues for all existing side chains
+        # skip non-existing sidechains and CYS (if cystein_indices) and update
+        # incomplete_sidechains and cystein_indices
+        for i,r in enumerate(res_list):
+
+            if cystein_indices is not None and rotamer_ids[i] == sidechain.CYS:
+                cystein_indices.append(i)
+                continue
 
-                if rotamer_ids[i] == sidechain.ALA or rotamer_ids[i] == sidechain.GLY:
-                    continue #no sidechain to model
+            if rotamer_ids[i] == sidechain.ALA or rotamer_ids[i] == sidechain.GLY:
+                continue # no sidechain to model
 
+            try:
+                frame_residue = sidechain.ConstructSidechainFrameResidue(\
+                                    r.handle, rotamer_ids[i], i, rotamer_settings)
+                frame_residues.append(frame_residue)
+            except:
                 incomplete_sidechains.append(i)
+    else:
+        # no frame residues to create, just update incomplete_sidechains
+        # and cystein_indices if needed
+        for i,r in enumerate(res_list):
 
-        #let's generate the frame without any cystein sidechains
-        #this is required for the disulfid score evaluation
-        frame = sidechain.Frame(frame_residues)
-
-        #some info we have to keep track of when evaluating disulfid bonds
-        cystein_rotamers = list()
-        disulfid_indices = list()
-        cys_ca_positions = list()
-        cys_cb_positions = list()
-
-        for i in cystein_indices:
-
-            rot_grop = None
-            r = prot.residues[i]
-            ca = r.FindAtom("CA")
-            cb = r.FindAtom("CB")
-
-            if not (ca.IsValid() and cb.IsValid()):
+            if cystein_indices is not None and rotamer_ids[i] == sidechain.CYS:
+                cystein_indices.append(i)
                 continue
 
-            cys_ca_positions.append(ca.GetPos())
-            cys_cb_positions.append(cb.GetPos())
-
-            if use_frm:
-                if bbdep:
-                    rot_group = sidechain.ConstructFRMRotamerGroup(r.handle,
-                                                                   sidechain.CYD,i,
-                                                                   rotamer_library,
-                                                                   rotamer_settings,
-                                                                   phi_angles[i],
-                                                                   psi_angles[i])
-                else:
-                    rot_group = sidechain.ConstructFRMRotamerGroup(r.handle,
-                                                                   sidechain.CYD,i,
-                                                                   rotamer_library,
-                                                                   rotamer_settings)
-            else:
-                if bbdep:
-                    rot_group = sidechain.ConstructRRMRotamerGroup(r.handle,
-                                                                   sidechain.CYD,i,
-                                                                   rotamer_library,
-                                                                   rotamer_settings,
-                                                                   phi_angles[i],
-                                                                   psi_angles[i])
-                else:
-                    rot_group = sidechain.ConstructRRMRotamerGroup(r.handle,
-                                                                   sidechain.CYD,i,
-                                                                   rotamer_library,
-                                                                   rotamer_settings)
-
-            frame.AddFrameEnergy(rot_group)
-            cystein_rotamers.append(rot_group)
-
-        for i in range(len(cystein_rotamers)):
-            for j in range(i+1, len(cystein_rotamers)):
-
-                if geom.Distance(cys_ca_positions[i], cys_ca_positions[j]) > 8.0:
-                    continue #they're too far for a disulfid bond
-
-                if cystein_indices[i] in disulfid_indices or cystein_indices[j] in disulfid_indices:
-                    continue #one of them already participates in a disulfid bond!
-                             #that one might be bether but right now is first come
-                             #first served principle
-
-                min_score = float("inf")
-                min_index_k = -1
-                min_index_l = -1
-
-                for k in range(len(cystein_rotamers[i])):
-                    for l in range(len(cystein_rotamers[j])):
-                        score = sidechain.DisulfidScore(cystein_rotamers[i][k], 
-                                                        cystein_rotamers[j][l],
-                                                        cys_ca_positions[i], 
-                                                        cys_cb_positions[i],
-                                                        cys_ca_positions[j], 
-                                                        cys_cb_positions[j])
-                        if score < min_score:
-                            min_index_k = k
-                            min_index_l = l
-                            min_score = score
-
-                if min_score < 45.0:
-                    disulfid_indices.append(cystein_indices[i])
-                    disulfid_indices.append(cystein_indices[j])
-                    particle_list = list()
-                    particle_list.append(cystein_rotamers[i][min_index_k][0])
-                    new_frame_residue_i = sidechain.FrameResidue(particle_list,cystein_indices[i])
-                    particle_list = list()
-                    particle_list.append(cystein_rotamers[j][min_index_l][0])
-                    new_frame_residue_j = sidechain.FrameResidue(particle_list,cystein_indices[j])
-                    frame_residues.append(new_frame_residue_i)
-                    frame_residues.append(new_frame_residue_j)
-                    #set the position in the proteins residues
-                    cystein_rotamers[i][min_index_k].ApplyOnResidue(prot.residues[cystein_indices[i]].handle,
-                                                                    consider_hydrogens=False)
-                    cystein_rotamers[j][min_index_l].ApplyOnResidue(prot.residues[cystein_indices[j]].handle,
-                                                                    consider_hydrogens=False) 
-                    sidechain.ConnectSidechain(prot.residues[cystein_indices[i]].handle,sidechain.CYS)
-                    sidechain.ConnectSidechain(prot.residues[cystein_indices[j]].handle,sidechain.CYS)
-
-
-        #All cysteins participating in a disulfid bond have been applied to the
-        #structure and added to the frame.
-        #All remaining ones have to be handled according the given flags.
-        for idx in cystein_indices:
-            if idx in disulfid_indices:
-                continue #it's all fine
-            if keep_sidechains:
-                try:
-                    frame_residue = sidechain.ConstructSidechainFrameResidue(prot.residues[idx].handle,
-                                                                             rotamer_ids[idx],
-                                                                             idx,rotamer_settings)
-                    frame_residues.append(frame_residue)
-                except:
-                    incomplete_sidechains.append(idx) 
-            else:
-                incomplete_sidechains.append(idx)  
-
-        #we finally have to rebuild the frame if any disulfid bonds have been built
-        if len(disulfid_indices) > 0:
-            frame = sidechain.Frame(frame_residues)                 
-
-
-    else:
+            if rotamer_ids[i] == sidechain.ALA or rotamer_ids[i] == sidechain.GLY:
+                continue # no sidechain to model
+
+            incomplete_sidechains.append(i)
+
+def _AddCysteinFrameResidues(frame_residues, incomplete_sidechains,
+                             keep_sidechains, res_list, rotamer_ids,
+                             rotamer_settings, cystein_indices,
+                             disulfid_indices, disulfid_rotamers):
+    '''Update frame_residues (list) with cysteins.
+    Parameters as in _AddSidechainFrameResidues.
+    Some cysteins (in disulfid_indices) get special treatment as disulfid
+    bridges (disulfid_indices, disulfid_rotamers from _GetDisulfidBridges).
+    '''
+    # handle cysteins participating in a disulfid bond
+    for cys_idx, cys_rot in zip(disulfid_indices, disulfid_rotamers):
+        # add FrameResidue
+        frame_residue = sidechain.FrameResidue([cys_rot[0]], cys_idx)
+        frame_residues.append(frame_residue)
+        # set the position in the proteins residues
+        cys_rot.ApplyOnResidue(res_list[cys_idx].handle,
+                               consider_hydrogens=False)
+        sidechain.ConnectSidechain(res_list[cys_idx].handle, sidechain.CYS)
+
+    # add remaining ones according the given flags
+    for idx in cystein_indices:
+        if idx in disulfid_indices:
+            continue # already handled
         if keep_sidechains:
-            for i,r in enumerate(prot.residues):
-                if rotamer_ids[i] == sidechain.ALA or rotamer_ids[i] == sidechain.GLY:
-                    continue #no sidechain to model
-                try:
-                    frame_residue = sidechain.ConstructSidechainFrameResidue(r.handle,
-                                                                             rotamer_ids[i],
-                                                                             i,rotamer_settings)
-                    frame_residues.append(frame_residue)
-                except:
-                    incomplete_sidechains.append(i)
+            try:
+                frame_residue = sidechain.ConstructSidechainFrameResidue(\
+                                    res_list[idx].handle, rotamer_ids[idx],
+                                    idx, rotamer_settings)
+                frame_residues.append(frame_residue)
+            except:
+                incomplete_sidechains.append(idx) 
         else:
-            for i,r in enumerate(prot.residues):
-                if rotamer_ids[i] not in [sidechain.GLY, sidechain.ALA]:
-                    incomplete_sidechains.append(i)
-
-        frame = sidechain.Frame(frame_residues)
+            incomplete_sidechains.append(idx)
 
-
-    #build rotamers for incomplete sidechains
-    for i in incomplete_sidechains:
-
-        r = prot.residues[i]
-        rot_id = rotamer_ids[i]
+def _GetRotamerGroup(res_handle, rot_id, res_idx, rot_lib, rot_settings,
+                     phi, psi, use_frm, bbdep):
+    '''Get RotamerGroup for res_handle according to settings.'''
+    if use_frm:
+        if bbdep:
+            return sidechain.ConstructFRMRotamerGroup(res_handle, rot_id,
+                                                      res_idx, rot_lib,
+                                                      rot_settings, phi, psi)
+        else:
+            return sidechain.ConstructFRMRotamerGroup(res_handle, rot_id,
+                                                      res_idx, rot_lib,
+                                                      rot_settings)
+    else:
+        if bbdep:
+            return sidechain.ConstructRRMRotamerGroup(res_handle, rot_id,
+                                                      res_idx, rot_lib,
+                                                      rot_settings, phi, psi)
+        else:
+            return sidechain.ConstructRRMRotamerGroup(res_handle, rot_id,
+                                                      res_idx, rot_lib,
+                                                      rot_settings)
+
+def _GetRotamerGroups(res_list, rot_ids, indices, rot_lib, rot_settings,
+                      phi_angles, psi_angles, use_frm, bbdep, frame_residues):
+    '''Get list of rotamer groups from subset of res_list.
+    Residues are chosen as res_list[i] for i in indices and only if a rotamer
+    group can be created (e.g. no ALA, GLY).
+    Rotamer groups are filtered to keep only best ones (given frame).
+    Returns list of rotamer groups and list of res. indices they belong to.
+    '''
+    # res.index (res_list[i]) for each modelled sc
+    residues_with_rotamer_group = list()
+    #  linked to residue in residues_with_rotamer_group
+    rotamer_groups = list()
+    # get frame for score evaluation
+    frame = sidechain.Frame(frame_residues)
+    # build rotamers for chosen sidechains
+    for i in indices:
+        # get rotamer ID
+        r = res_list[i]
+        rot_id = rot_ids[i]
         
-        if(rot_id == sidechain.ALA or rot_id == sidechain.GLY):
+        if rot_id == sidechain.ALA or rot_id == sidechain.GLY:
             continue
 
         if rot_id == sidechain.CYS:
@@ -300,7 +215,7 @@ def Reconstruct(ent, keep_sidechains = False, build_disulfids = True,
                 omega = tor.GetAngle()
             elif i > 0:
                 # fallback computation of omega as in OST-code
-                prev = prot.residues[i-1]
+                prev = res_list[i-1]
                 if prev.IsValid() and prev.IsPeptideLinking():
                     ca_prev = prev.FindAtom("CA")
                     c_prev = prev.FindAtom("C")
@@ -319,49 +234,232 @@ def Reconstruct(ent, keep_sidechains = False, build_disulfids = True,
                 else:
                     rot_id = sidechain.TPR
 
+        # get RotamerGroup
         try:
-            if use_frm:
-                if bbdep:
-                    rot_group = sidechain.ConstructFRMRotamerGroup(r.handle,rot_id,i,
-                                                                   rotamer_library,rotamer_settings,
-                                                                   phi_angles[i],psi_angles[i])
-                else:
-                    rot_group = sidechain.ConstructFRMRotamerGroup(r.handle,rot_id,i,
-                                                                   rotamer_library,rotamer_settings)
-            else:
-                if bbdep:
-                    rot_group = sidechain.ConstructRRMRotamerGroup(r.handle,rot_id,i,
-                                                                   rotamer_library,rotamer_settings,
-                                                                   phi_angles[i],psi_angles[i])
-                else:
-                    rot_group = sidechain.ConstructRRMRotamerGroup(r.handle,rot_id,i,
-                                                                   rotamer_library,rotamer_settings)
-
+            rot_group = _GetRotamerGroup(r.handle, rot_id, i, rot_lib,
+                                         rot_settings, phi_angles[i],
+                                         psi_angles[i], use_frm, bbdep)
         except:
             continue
+        # keep best ones
         rot_group.CalculateInternalEnergies()
         frame.SetFrameEnergy(rot_group)
         rot_group.ApplySelfEnergyThresh()
         rotamer_groups.append(rot_group)
-        residues_with_rotamer_group.append(i)   
+        residues_with_rotamer_group.append(i)
+
+    return rotamer_groups, residues_with_rotamer_group
+
+def _GetDisulfidBridges(frame_residues, cystein_indices, res_list, rotamer_library,
+                        use_frm, bbdep, rotamer_settings, phi_angles, psi_angles):
+    '''Get disulfid bridges for CYS and according rotamers.
+    CYS are identified by by items in cystein_indices (into res_list).
+    Returns: disulfid_indices: list of res. index in bridge,
+             disulfid_rotamers: list of rotamers (best one for bridge).
+    '''
+    # this is required for the disulfid score evaluation
+    frame = sidechain.Frame(frame_residues)
+
+    # some info we have to keep track of when evaluating disulfid bonds
+    cystein_rotamers = list()
+    cys_ca_positions = list()
+    cys_cb_positions = list()
+
+    for i in cystein_indices:
+        # check ca, cb
+        r = res_list[i]
+        ca = r.FindAtom("CA")
+        cb = r.FindAtom("CB")
+        if not (ca.IsValid() and cb.IsValid()):
+            continue
+        cys_ca_positions.append(ca.GetPos())
+        cys_cb_positions.append(cb.GetPos())
+        # get RotamerGroup
+        rot_group = _GetRotamerGroup(r.handle, sidechain.CYD, i, rotamer_library,
+                                     rotamer_settings, phi_angles[i],
+                                     psi_angles[i], use_frm, bbdep)
+        frame.AddFrameEnergy(rot_group)
+        cystein_rotamers.append(rot_group)
+
+    # get CYS with disulfid bonds and the chosen rotamers
+    disulfid_indices = list()
+    disulfid_rotamers = list()
+    for i in range(len(cystein_rotamers)):
+        for j in range(i+1, len(cystein_rotamers)):
+
+            # too far for a disulfid bond?
+            if geom.Distance(cys_ca_positions[i], cys_ca_positions[j]) > 8.0:
+                continue 
+
+            # already done? NOTE: new one might be better, but here we do
+            #                     first come, first served
+            if cystein_indices[i] in disulfid_indices \
+               or cystein_indices[j] in disulfid_indices:
+                continue
 
+            min_score = float("inf")
+            min_index_k = -1
+            min_index_l = -1
+
+            for k in range(len(cystein_rotamers[i])):
+                for l in range(len(cystein_rotamers[j])):
+                    score = sidechain.DisulfidScore(cystein_rotamers[i][k],
+                                                    cystein_rotamers[j][l],
+                                                    cys_ca_positions[i], 
+                                                    cys_cb_positions[i],
+                                                    cys_ca_positions[j],
+                                                    cys_cb_positions[j])
+                    if score < min_score:
+                        min_index_k = k
+                        min_index_l = l
+                        min_score = score
+
+            if min_score < 45.0:
+                # update indices
+                cys_idx_i = cystein_indices[i]
+                cys_idx_j = cystein_indices[j]
+                cys_rot_i = cystein_rotamers[i][min_index_k]
+                cys_rot_j = cystein_rotamers[j][min_index_l]
+                disulfid_indices.append(cys_idx_i)
+                disulfid_indices.append(cys_idx_j)
+                disulfid_rotamers.append(cys_rot_i)
+                disulfid_rotamers.append(cys_rot_j)
+
+    return disulfid_indices, disulfid_rotamers
+
+###############################################################################
+
+def Reconstruct(ent, keep_sidechains=False, build_disulfids=True,
+                rotamer_model="frm", consider_hbonds=True,
+                consider_ligands=True, rotamer_library=None):
+    '''Reconstruct sidechains for the given structure.
+
+    :param ent:          Structure for sidechain reconstruction. Note, that the
+                         sidechain reconstruction gets directly applied on the
+                         structure itself.
+
+    :param keep_sidechains: Flag, whether complete sidechains in *ent* (i.e. 
+                            containing all required atoms) should be kept rigid
+                            and directly be added to the frame.
+
+    :param build_disulfids: Flag, whether possible disulfid bonds should be 
+                            searched. If a disulfid bond is found, the two
+                            participating cysteins are fixed and added to
+                            the frame.
+
+    :param rotamer_model: Rotamer model to be used, can either be "frm" or "rrm"
+
+    :param consider_hbonds: Flag, whether hbonds should be evaluated in the
+                            energy function. If set to False, no hydrogens will
+                            be built when building rotamers and frame and the
+                            **add_polar_hydrogens** flag won't have any
+                            consequences.
+
+    :param consider_ligands: Flag, whether to add ligands (anything in chain
+                             '_') as static objects.
+
+    :param rotamer_library: A rotamer library to extract the rotamers from. The
+                            default is the :meth:`Dunbrack <LoadDunbrackLib>`
+                            library.
+
+
+    :type ent:            :class:`ost.mol.EntityHandle`
+    :type keep_sidechains: :class:`bool`
+    :type build_disulfids: :class:`bool`
+    :type rotamer_model:   :class:`str`
+    :type consider_hbonds: :class:`bool`
+    :type consider_ligands: :class:`bool`
+    :type rotamer_library: :class:`BBDepRotamerLib` / :class:`RotamerLib`
+    '''
+
+    # setup settings
+    if rotamer_model.lower() == "frm":
+        use_frm = True
+    elif rotamer_model.lower() == "rrm":
+        use_frm = False
+    else:
+        raise RuntimeError("Only \"rrm\" and \"frm\" allowed for rotamer_model!")
 
+    rotamer_settings = sidechain.RotamerSettings()
+    rotamer_settings.consider_hbonds = consider_hbonds
+    if rotamer_library == None: 
+        rotamer_library = sidechain.LoadDunbrackLib()
+    bbdep = False
+    if type(rotamer_library) is sidechain.BBDepRotamerLib:
+        bbdep = True
+    
+    # take out ligand chain and any non-peptides
+    prot = ent.Select("peptide=true and cname!='_'")
+    
+    # parse residues (all lists of length len(prot.residues))
+    rotamer_ids = _GetRotamerIDs(prot.residues)
+    phi_angles, psi_angles = _GetDihedrals(prot.residues)
+
+    # set nter and cter (needed in _AddBackboneFrameResidues)
+    for c in prot.chains:
+        c.residues[0].SetIntProp("n_ter",1)
+        c.residues[-1].SetIntProp("c_ter",1)
+
+    # build up frame
+    frame_residues = list()         # list of frame residues connected to frame
+    incomplete_sidechains = list()  # residue indices
+    _AddBackboneFrameResidues(frame_residues, prot.residues, rotamer_ids,
+                              rotamer_settings, phi_angles)
+    
+    # add ligands?
+    if consider_ligands:
+        ligs = ent.Select("cname='_'")
+        offset = len(prot.residues)
+        _AddLigandFrameResidues(frame_residues, ligs, rotamer_settings, offset)
+
+    # check special handling of cysteins
+    if build_disulfids:
+        # residue indices of cysteins
+        cystein_indices = list()
+        # update frame_residues, incomplete_sidechains, cystein_indices
+        _AddSidechainFrameResidues(frame_residues, incomplete_sidechains,
+                                   keep_sidechains, prot.residues, rotamer_ids,
+                                   rotamer_settings, cystein_indices)
+        # update frame_residues, incomplete_sidechains with cysteins (if needed)
+        if len(cystein_indices) > 0:
+            # get disulfid bridges and according rotamers
+            disulfid_indices, disulfid_rotamers = \
+                _GetDisulfidBridges(frame_residues, cystein_indices, prot.residues,
+                                    rotamer_library, use_frm, bbdep,
+                                    rotamer_settings, phi_angles, psi_angles)
+            # update frame_residues, incomplete_sidechains
+            _AddCysteinFrameResidues(frame_residues, incomplete_sidechains,
+                                     keep_sidechains, prot.residues, rotamer_ids,
+                                     rotamer_settings, cystein_indices,
+                                     disulfid_indices, disulfid_rotamers)
+    else:
+        # update frame_residues, incomplete_sidechains
+        _AddSidechainFrameResidues(frame_residues, incomplete_sidechains,
+                                   keep_sidechains, prot.residues, rotamer_ids,
+                                   rotamer_settings)
+    
+    # get rotamer groups and residues they're linked to
+    rotamer_groups, residues_with_rotamer_group = \
+        _GetRotamerGroups(prot.residues, rotamer_ids, incomplete_sidechains,
+                          rotamer_library, rotamer_settings, phi_angles,
+                          psi_angles, use_frm, bbdep, frame_residues)
+
+    # set up graph and solve to get best rotamers
     if use_frm:
         graph = sidechain.Graph.CreateFromFRMList(rotamer_groups)
     else:
         graph = sidechain.Graph.CreateFromRRMList(rotamer_groups)
 
-
     solution = graph.Solve(100000000,0.02)
 
-
+    # update structure
     for i,rot_group,sol in zip(residues_with_rotamer_group,rotamer_groups,solution):
         try:
-            rot_group[sol].ApplyOnResidue(prot.residues[i].handle,consider_hydrogens=False)
-            sidechain.ConnectSidechain(prot.residues[i].handle,rotamer_ids[i])
+            res_handle = prot.residues[i].handle
+            rot_group[sol].ApplyOnResidue(res_handle, consider_hydrogens=False)
+            sidechain.ConnectSidechain(res_handle, rotamer_ids[i])
         except:
-            print "there is a backbone atom missing... ",prot.residues[i].GetQualifiedName()
-
+            print "there is a backbone atom missing... ", res_handle.GetQualifiedName()
 
 # these methods will be exported into module
 __all__ = ('Reconstruct',)