diff --git a/modelling/pymod/_closegaps.py b/modelling/pymod/_closegaps.py
index 1fc457c4addf565df206c4ee4b9d74fc2ca024fa..fd6a10da1765e180150ebb7d62b056c640df115e 100644
--- a/modelling/pymod/_closegaps.py
+++ b/modelling/pymod/_closegaps.py
@@ -53,18 +53,11 @@ def _GetBestLC(mhandle, loop_candidates, start_resnum, chain_idx,
     all_scores = db_scores.Copy()  # copy to avoid changing db_scores
 
     # add backbone scores
-    for key in ScoringWeights.GetBackboneScoringKeys():
-        scores = loop_candidates.CalculateScores(mhandle.backbone_scorer,
-                                                 key, start_resnum, chain_idx)
-        all_scores.Set(key, scores)
+    loop_candidates.CalculateBackboneScores(all_scores, mhandle.backbone_scorer,
+                                            start_resnum, chain_idx)
 
     # add all atom scores?
     if with_aa:
-        # get stuff from mhandle
-        aa_scorer = mhandle.all_atom_scorer
-        aa_score_env = mhandle.all_atom_scorer_env
-        aa_sc_env = mhandle.all_atom_sidechain_env
-        sc_rec = mhandle.sidechain_reconstructor
         # get best ones based on previous scores
         non_aa_weights = ScoringWeights.GetWeights(with_db, False)
         non_aa_scores = all_scores.LinearCombine(non_aa_weights)
@@ -73,34 +66,9 @@ def _GetBestLC(mhandle, loop_candidates, start_resnum, chain_idx,
         # extract relevant loop candidates and scores
         aa_loop_candidates = loop_candidates.Extract(min_indices)
         all_scores = all_scores.Extract(min_indices)
-        
-        # TODO: use LoopCandidates feature to do all of the below...
-
-        # keep backup of loop of interest
-        loop_length = len(loop_candidates[0])
-        loop_backup = aa_sc_env.GetEnvironment(start_resnum, loop_length,
-                                               chain_idx)
-        # setup empty lists for AA results
-        aa_scoring_keys = ScoringWeights.GetAllAtomScoringKeys()
-        aa_scores = {key: list() for key in aa_scoring_keys}
-        # score them w/o relaxation
-        for lc_bb_list in aa_loop_candidates:
-            # setup sc env. and reconstruct
-            aa_sc_env.SetEnvironment(lc_bb_list, start_resnum, chain_idx)
-            sc_result = sc_rec.Reconstruct(start_resnum, loop_length, chain_idx)
-            # fix env. for scoring (always contains full loop too!)
-            aa_score_env.SetEnvironment(sc_result.env_pos)
-            for key in aa_scoring_keys:
-                score = aa_scorer[key].CalculateScore(start_resnum,
-                                                      loop_length, chain_idx)
-                aa_scores[key].append(score)
-        # reset aa_sc_env
-        aa_sc_env.SetEnvironment(loop_backup)
-        # NOTE: we leave the aa_score_env in a modified state
-        # -> this is ok for aa_score_env as we always set proper surroundings
-        # -> this is not ok for SC (due to frame approach), so we back it up
-        for key in aa_scoring_keys:
-            all_scores.Set(key, aa_scores[key])
+        # add all atom scores
+        aa_loop_candidates.CalculateAllAtomScores(all_scores, mhandle,
+                                                  start_resnum, chain_idx)
 
     # get / return best
     scores = all_scores.LinearCombine(weights)
@@ -768,8 +736,7 @@ def FillLoopsByDatabase(mhandle, fragment_db, structure_db,
             # check for stem rmsd before ApplyCCD
             if add_db_features:
                 db_scores = ScoreContainer()
-                db_scores.Set(ScoringWeights.GetStemRMSDsKey(),
-                              candidates.CalculateStemRMSDs(n_stem, c_stem))
+                candidates.CalculateStemRMSDs(db_scores, n_stem, c_stem)
             # try to close loops
             try:
                 #pylint: disable=broad-except
@@ -796,21 +763,18 @@ def FillLoopsByDatabase(mhandle, fragment_db, structure_db,
 
             # deal with DB features
             if add_db_features:
-                # setup
-                assert(len(orig_indices) == len(candidates))
-                seq_prof_key = ScoringWeights.GetSequenceProfileScoresKey()
-                str_prof_key = ScoringWeights.GetStructureProfileScoresKey()
                 # get subset of stem rmsd scores
+                assert(len(orig_indices) == len(candidates))
                 db_scores = db_scores.Extract(orig_indices)
                 # add profile scores
                 prof = mhandle.profiles[actual_gap.GetChainIndex()]
                 start_pos = n_stem.GetNumber().GetNum() - 1
-                scores = candidates.CalculateSequenceProfileScores(\
-                               structure_db, prof, start_pos)
-                db_scores.Set(seq_prof_key, scores)
-                scores = candidates.CalculateStructureProfileScores(\
-                               structure_db, prof, start_pos)
-                db_scores.Set(str_prof_key, scores)
+                candidates.CalculateSequenceProfileScores(db_scores,
+                                                          structure_db,
+                                                          prof, start_pos)
+                candidates.CalculateStructureProfileScores(db_scores,
+                                                           structure_db,
+                                                           prof, start_pos)
                 # update list
                 actual_db_scores.append(db_scores)
 
diff --git a/modelling/pymod/export_loop_candidate.cc b/modelling/pymod/export_loop_candidate.cc
index 15c7e81ba0620b8d76ef351c47dce2cd3c6bb40a..16292170844dd9a9fd68113bb50263bba2e53727 100644
--- a/modelling/pymod/export_loop_candidate.cc
+++ b/modelling/pymod/export_loop_candidate.cc
@@ -161,7 +161,40 @@ boost::python::list WrapCalculateLinearScores(LoopCandidatesPtr p,
   return return_list;
 }
 
-boost::python::list WrapCalculateSequenceProfileScores(LoopCandidatesPtr p,
+void WrapCalculateBackboneScores(LoopCandidatesPtr p,
+                                 ScoreContainer& score_container,
+                                 scoring::BackboneOverallScorerPtr scorer,
+                                 uint start_resnum, uint chain_idx) {
+  p->CalculateBackboneScores(score_container, scorer, start_resnum, chain_idx);
+}
+void WrapCalculateBackboneScoresK(LoopCandidatesPtr p,
+                                  ScoreContainer& score_container,
+                                  scoring::BackboneOverallScorerPtr scorer,
+                                  const boost::python::list& keys,
+                                  uint start_resnum, uint chain_idx) {
+  std::vector<String> v_keys;
+  core::ConvertListToVector(keys, v_keys);
+  p->CalculateBackboneScores(score_container, scorer, v_keys, start_resnum,
+                             chain_idx);
+}
+void WrapCalculateAllAtomScores(LoopCandidatesPtr p,
+                                ScoreContainer& score_container,
+                                const ModellingHandle& mhandle,
+                                uint start_resnum, uint chain_idx) {
+  p->CalculateAllAtomScores(score_container, mhandle, start_resnum, chain_idx);
+}
+void WrapCalculateAllAtomScoresK(LoopCandidatesPtr p,
+                                 ScoreContainer& score_container,
+                                 const ModellingHandle& mhandle,
+                                 const boost::python::list& keys,
+                                 uint start_resnum, uint chain_idx) {
+  std::vector<String> v_keys;
+  core::ConvertListToVector(keys, v_keys);
+  p->CalculateAllAtomScores(score_container, mhandle, v_keys, start_resnum,
+                            chain_idx);
+}
+
+boost::python::list WrapCalcSequenceProfScores(LoopCandidatesPtr p,
                                loop::StructureDBPtr structure_db,
                                const ost::seq::ProfileHandle& prof,
                                uint offset) {
@@ -173,8 +206,16 @@ boost::python::list WrapCalculateSequenceProfileScores(LoopCandidatesPtr p,
   core::AppendVectorToList(scores, return_list);
   return return_list;
 }
+void WrapCalcSequenceProfScoresSC(LoopCandidatesPtr p,
+                                  ScoreContainer& score_container,
+                                  loop::StructureDBPtr structure_db,
+                                  const ost::seq::ProfileHandle& prof,
+                                  uint offset) {
+  p->CalculateSequenceProfileScores(score_container, structure_db, prof,
+                                    offset);
+}
 
-boost::python::list WrapCalculateStructureProfileScores(LoopCandidatesPtr p,
+boost::python::list WrapCalcStructureProfScores(LoopCandidatesPtr p,
                                loop::StructureDBPtr structure_db,
                                const ost::seq::ProfileHandle& prof,
                                uint offset) {
@@ -186,6 +227,14 @@ boost::python::list WrapCalculateStructureProfileScores(LoopCandidatesPtr p,
   core::AppendVectorToList(scores, return_list);
   return return_list;
 }
+void WrapCalcStructureProfScoresSC(LoopCandidatesPtr p,
+                                   ScoreContainer& score_container,
+                                   loop::StructureDBPtr structure_db,
+                                   const ost::seq::ProfileHandle& prof,
+                                   uint offset) {
+  p->CalculateStructureProfileScores(score_container, structure_db, prof,
+                                     offset);
+}
 
 boost::python::list WrapCalculateStemRMSDs(LoopCandidatesPtr p,
                                const ost::mol::ResidueHandle& n_stem,
@@ -197,6 +246,12 @@ boost::python::list WrapCalculateStemRMSDs(LoopCandidatesPtr p,
   core::AppendVectorToList(scores, return_list);
   return return_list;
 }
+void WrapCalculateStemRMSDsSC(LoopCandidatesPtr p,
+                              ScoreContainer& score_container,
+                              const ost::mol::ResidueHandle& n_stem,
+                              const ost::mol::ResidueHandle& c_stem) {
+  p->CalculateStemRMSDs(score_container, n_stem, c_stem);
+}
 
 LoopCandidatesPtr WrapExtractList(LoopCandidatesPtr p,
                                   const boost::python::list& idx_list) {
@@ -243,12 +298,32 @@ void export_loop_candidate() {
     .def("CalculateLinearScores", &WrapCalculateLinearScores,
          (arg("scorer"), arg("linear_weights"), arg("start_resnum"), 
           arg("chain_index")=0))
-    .def("CalculateSequenceProfileScores", &WrapCalculateSequenceProfileScores,
+    .def("CalculateBackboneScores", &WrapCalculateBackboneScores,
+         (arg("score_container"), arg("scorer"),
+          arg("start_resnum"), arg("chain_index")=0))
+    .def("CalculateBackboneScores", &WrapCalculateBackboneScoresK,
+         (arg("score_container"), arg("scorer"), arg("keys"),
+          arg("start_resnum"), arg("chain_index")=0))
+    .def("CalculateAllAtomScores", &WrapCalculateAllAtomScores,
+         (arg("score_container"), arg("mhandle"),
+          arg("start_resnum"), arg("chain_index")=0))
+    .def("CalculateAllAtomScores", &WrapCalculateAllAtomScoresK,
+         (arg("score_container"), arg("mhandle"), arg("keys"),
+          arg("start_resnum"), arg("chain_index")=0))
+    .def("CalculateSequenceProfileScores", &WrapCalcSequenceProfScores,
          (arg("structure_db"), arg("prof"), arg("offset")=0))
-    .def("CalculateStructureProfileScores", &WrapCalculateStructureProfileScores,
+    .def("CalculateSequenceProfileScores", &WrapCalcSequenceProfScoresSC,
+         (arg("score_container"), arg("structure_db"), arg("prof"),
+          arg("offset")=0))
+    .def("CalculateStructureProfileScores", &WrapCalcStructureProfScores,
          (arg("structure_db"), arg("prof"), arg("offset")=0))
+    .def("CalculateStructureProfileScores", &WrapCalcStructureProfScoresSC,
+         (arg("score_container"), arg("structure_db"), arg("prof"),
+          arg("offset")=0))
     .def("CalculateStemRMSDs", &WrapCalculateStemRMSDs,
          (arg("n_stem"), arg("c_stem")))
+    .def("CalculateStemRMSDs", &WrapCalculateStemRMSDsSC,
+         (arg("score_container"), arg("n_stem"), arg("c_stem")))
     .def("__iter__", iterator<LoopCandidates>())
     .def("__len__", &LoopCandidates::size)
     .def("__getitem__", lc_getitem, 
diff --git a/modelling/src/loop_candidate.cc b/modelling/src/loop_candidate.cc
index ccbe9628e4648def3bc4fab7c5ffe280c83ad17a..1be73cad4c97461437b86c3060ac37f48533e877 100644
--- a/modelling/src/loop_candidate.cc
+++ b/modelling/src/loop_candidate.cc
@@ -1,5 +1,6 @@
 #include <promod3/modelling/loop_candidate.hh>
 #include <promod3/modelling/kic.hh>
+#include <promod3/modelling/scoring_weights.hh>
 #include <promod3/core/message.hh>
 #include <promod3/core/runtime_profiling.hh>
 #include <promod3/core/geom_stems.hh>
@@ -438,12 +439,95 @@ std::vector<Real> LoopCandidates::CalculateLinearScores(
                                             start_resnum, chain_idx);
 }
 
+void LoopCandidates::CalculateBackboneScores(
+                    ScoreContainer& score_container,
+                    scoring::BackboneOverallScorerPtr scorer,
+                    uint start_resnum, uint chain_idx) const {
+  const std::vector<String>&
+  keys = ScoringWeights::GetInstance().GetBackboneScoringKeys();
+  CalculateBackboneScores(score_container, scorer, keys, start_resnum,
+                          chain_idx);
+}
+void LoopCandidates::CalculateBackboneScores(
+                    ScoreContainer& score_container,
+                    scoring::BackboneOverallScorerPtr scorer,
+                    const std::vector<String>& keys,
+                    uint start_resnum, uint chain_idx) const {
+  for (uint i = 0; i < keys.size(); ++i) {
+    std::vector<Real> scores = CalculateScores(scorer, keys[i], start_resnum,
+                                               chain_idx);
+    score_container.Set(keys[i], scores);
+  }
+}
+
+void LoopCandidates::CalculateAllAtomScores(
+                    ScoreContainer& score_container,
+                    const ModellingHandle& mhandle,
+                    uint start_resnum, uint chain_idx) const {
+  const std::vector<String>&
+  keys = ScoringWeights::GetInstance().GetAllAtomScoringKeys();
+  CalculateAllAtomScores(score_container, mhandle, keys, start_resnum,
+                         chain_idx);
+}
+void LoopCandidates::CalculateAllAtomScores(
+                    ScoreContainer& score_container,
+                    const ModellingHandle& mhandle,
+                    const std::vector<String>& keys,
+                    uint start_resnum, uint chain_idx) const {
+
+  core::ScopedTimerPtr timer = core::StaticRuntimeProfiler::StartScoped(
+                          "LoopCandidates::CalculateAllAtomScores", 2);
+
+  // get stuff from mhandle (short names for convenience)
+  loop::AllAtomEnvPtr aa_score_env = mhandle.all_atom_scorer_env;
+  scoring::AllAtomOverallScorerPtr aa_scorer = mhandle.all_atom_scorer;
+  loop::AllAtomEnvPtr aa_sc_env = mhandle.all_atom_sidechain_env;
+  sidechain::SidechainReconstructorPtr sc_rec = mhandle.sidechain_reconstructor;
+  sidechain::SidechainReconstructionDataPtr sc_result;
+
+  // prepare storage for all scores
+  std::vector< std::vector<Real> > aa_scores(keys.size());
+  for (uint i_k = 0; i_k < keys.size(); ++i_k) {
+    aa_scores[i_k].resize(size());
+  }
+
+  // keep backup of loop of interest
+  const uint loop_length = sequence_.size();
+  loop::AllAtomEnvPositionsPtr loop_backup;
+  loop_backup = aa_sc_env->GetEnvironment(start_resnum, loop_length, chain_idx);
+
+  // for each loop: reconstruct sidechains and get all scores
+  for (uint i_lc = 0; i_lc < size(); ++i_lc) {
+    // setup sidechain env. and reconstruct
+    aa_sc_env->SetEnvironment(candidates_[i_lc], start_resnum, chain_idx);
+    sc_result = sc_rec->Reconstruct(start_resnum, loop_length, chain_idx);
+    // fix env. for scoring (always contains loop and relevant surrounding!)
+    aa_score_env->SetEnvironment(*(sc_result->env_pos));
+    for (uint i_k = 0; i_k < keys.size(); ++i_k) {
+      scoring::AllAtomScorerPtr scorer = aa_scorer->Get(keys[i_k]);
+      aa_scores[i_k][i_lc] = scorer->CalculateScore(start_resnum, loop_length,
+                                                    chain_idx);
+    }
+  }
+
+  // update score container
+  for (uint i_k = 0; i_k < keys.size(); ++i_k) {
+    score_container.Set(keys[i_k], aa_scores[i_k]);
+  }
+
+  // reset aa_sc_env
+  aa_sc_env->SetEnvironment(*loop_backup);
+  // NOTE: we leave the aa_score_env in a modified state
+  // -> this is ok for aa_score_env as we always set proper surroundings
+  // -> this is not ok for SC (due to frame approach), so we back it up
+}
+
 std::vector<Real> LoopCandidates::CalculateSequenceProfileScores(
                     loop::StructureDBPtr structure_db,
                     const ost::seq::ProfileHandle& prof, uint offset) const {
 
   core::ScopedTimerPtr timer = core::StaticRuntimeProfiler::StartScoped(
-                                "LoopCandidates::CalculateSequenceProfileScores", 2);
+                          "LoopCandidates::CalculateSequenceProfileScores", 2);
   // check
   if (!HasFragmentInfos()) {
     throw promod3::Error("FragmentInfos are invalid in LoopCandidates!");
@@ -459,13 +543,22 @@ std::vector<Real> LoopCandidates::CalculateSequenceProfileScores(
   return scores;
 
 }
+void LoopCandidates::CalculateSequenceProfileScores(
+                    ScoreContainer& score_container,
+                    loop::StructureDBPtr structure_db,
+                    const ost::seq::ProfileHandle& prof, uint offset) const {
+  // use default key
+  score_container.Set(
+            ScoringWeights::GetInstance().GetSequenceProfileScoresKey(),
+            CalculateSequenceProfileScores(structure_db, prof, offset));
+}
 
 std::vector<Real> LoopCandidates::CalculateStructureProfileScores(
                     loop::StructureDBPtr structure_db,
                     const ost::seq::ProfileHandle& prof, uint offset) const {
 
   core::ScopedTimerPtr timer = core::StaticRuntimeProfiler::StartScoped(
-                                "LoopCandidates::CalculateStructureProfileScores", 2);
+                          "LoopCandidates::CalculateStructureProfileScores", 2);
   // check
   if (!HasFragmentInfos()) {
     throw promod3::Error("FragmentInfos are invalid in LoopCandidates!");
@@ -480,6 +573,15 @@ std::vector<Real> LoopCandidates::CalculateStructureProfileScores(
   }
   return scores;
 }
+void LoopCandidates::CalculateStructureProfileScores(
+                    ScoreContainer& score_container,
+                    loop::StructureDBPtr structure_db,
+                    const ost::seq::ProfileHandle& prof, uint offset) const {
+  // use default key
+  score_container.Set(
+            ScoringWeights::GetInstance().GetStructureProfileScoresKey(),
+            CalculateStructureProfileScores(structure_db, prof, offset));
+}
 
 std::vector<Real> LoopCandidates::CalculateStemRMSDs(
                     const ost::mol::ResidueHandle& n_stem, 
@@ -514,5 +616,13 @@ std::vector<Real> LoopCandidates::CalculateStemRMSDs(
   }
   return scores;
 }
+void LoopCandidates::CalculateStemRMSDs(
+                    ScoreContainer& score_container,
+                    const ost::mol::ResidueHandle& n_stem, 
+                    const ost::mol::ResidueHandle& c_stem) const {
+  // use default key
+  score_container.Set(ScoringWeights::GetInstance().GetStemRMSDsKey(),
+                      CalculateStemRMSDs(n_stem, c_stem));
+}
 
 }}
diff --git a/modelling/src/loop_candidate.hh b/modelling/src/loop_candidate.hh
index f10819e1a214486e03b402383fd26c4ce4674a02..1504b0b4c1a003ce2b069262c84f46b662aba9af 100644
--- a/modelling/src/loop_candidate.hh
+++ b/modelling/src/loop_candidate.hh
@@ -6,7 +6,10 @@
 #include <promod3/loop/torsion_sampler.hh>
 #include <promod3/scoring/backbone_overall_scorer.hh>
 #include <promod3/modelling/ccd.hh>
+#include <promod3/modelling/model.hh>
 #include <promod3/modelling/monte_carlo.hh>
+#include <promod3/modelling/score_container.hh>
+
 #include <ost/seq/profile_handle.hh>
 #include <vector>
 
@@ -84,24 +87,52 @@ public:
 
   std::vector<Real>
   CalculateScores(scoring::BackboneOverallScorerPtr scorer,
-                  const String& key, uint start_resnum, uint chain_idx) const;
+                  const String& key, uint start_resnum,
+                  uint chain_idx = 0) const;
   std::vector<Real>
   CalculateLinearScores(scoring::BackboneOverallScorerPtr scorer,
                         const std::map<String, Real>& linear_weights,
-                        uint start_resnum, uint chain_idx) const;
+                        uint start_resnum, uint chain_idx = 0) const;
+  
+  void CalculateBackboneScores(ScoreContainer& score_container,
+                               scoring::BackboneOverallScorerPtr scorer,
+                               uint start_resnum, uint chain_idx = 0) const;
+  void CalculateBackboneScores(ScoreContainer& score_container,
+                               scoring::BackboneOverallScorerPtr scorer,
+                               const std::vector<String>& keys,
+                               uint start_resnum, uint chain_idx = 0) const;
+  void CalculateAllAtomScores(ScoreContainer& score_container,
+                              const ModellingHandle& mhandle,
+                              uint start_resnum, uint chain_idx = 0) const;
+  void CalculateAllAtomScores(ScoreContainer& score_container,
+                              const ModellingHandle& mhandle,
+                              const std::vector<String>& keys,
+                              uint start_resnum, uint chain_idx = 0) const;
 
   std::vector<Real>
   CalculateSequenceProfileScores(loop::StructureDBPtr structure_db,
                                  const ost::seq::ProfileHandle& prof,
                                  uint offset = 0) const;
+  void CalculateSequenceProfileScores(ScoreContainer& score_container,
+                                      loop::StructureDBPtr structure_db,
+                                      const ost::seq::ProfileHandle& prof,
+                                      uint offset = 0) const;
+
   std::vector<Real>
   CalculateStructureProfileScores(loop::StructureDBPtr structure_db,
                                   const ost::seq::ProfileHandle& prof,
                                   uint offset = 0) const;
+  void CalculateStructureProfileScores(ScoreContainer& score_container,
+                                       loop::StructureDBPtr structure_db,
+                                       const ost::seq::ProfileHandle& prof,
+                                       uint offset = 0) const;
 
   std::vector<Real>
   CalculateStemRMSDs(const ost::mol::ResidueHandle& n_stem, 
                      const ost::mol::ResidueHandle& c_stem) const;
+  void CalculateStemRMSDs(ScoreContainer& score_container,
+                          const ost::mol::ResidueHandle& n_stem, 
+                          const ost::mol::ResidueHandle& c_stem) const;
 
   // Add/Remove candidates
   void Add(const loop::BackboneList& candidate);
diff --git a/modelling/src/model.cc b/modelling/src/model.cc
index 9e2b8b26abab6ad078186f4b3929efdcc47d5f08..57d1601730ad8f51c09795df5e2320a339161a5a 100644
--- a/modelling/src/model.cc
+++ b/modelling/src/model.cc
@@ -95,7 +95,7 @@ ModellingHandle ModellingHandle::Copy() const{
 }
 
 
-int CountEnclosedGaps(ModellingHandle& mhandle, const StructuralGap& gap,
+int CountEnclosedGaps(const ModellingHandle& mhandle, const StructuralGap& gap,
                       bool insertions_only)
 {
   core::ScopedTimerPtr prof = core::StaticRuntimeProfiler::StartScoped(
@@ -109,7 +109,7 @@ int CountEnclosedGaps(ModellingHandle& mhandle, const StructuralGap& gap,
   if (!gap.IsNTerminal()) my_b_res = gap.before.GetNumber();
   ResNum my_a_res = my_b_res + gap.GetLength() + 1;
 
-  for (StructuralGapList::iterator i = mhandle.gaps.begin();
+  for (StructuralGapList::const_iterator i = mhandle.gaps.begin();
        i != mhandle.gaps.end(); ++i) {
     
     // skip different chains
diff --git a/modelling/src/model.hh b/modelling/src/model.hh
index 9637b1763415802139f6c215290145c7bb786cff..b1bd18e45d22bd142f2292e5954e4bcdd2e8e86e 100644
--- a/modelling/src/model.hh
+++ b/modelling/src/model.hh
@@ -49,7 +49,7 @@ struct ModellingHandle {
 };
 
 // see Python doc
-int CountEnclosedGaps(ModellingHandle& mhandle, const StructuralGap& gap,
+int CountEnclosedGaps(const ModellingHandle& mhandle, const StructuralGap& gap,
                       bool insertions_only=false);
 
 // see Python doc