diff --git a/core/src/cluster.cc b/core/src/cluster.cc
index 4ccbee05011a57720d975d38a52439ef0188064b..3c579730429acd631911538e2df574101ff955cf 100644
--- a/core/src/cluster.cc
+++ b/core/src/cluster.cc
@@ -19,12 +19,12 @@ struct Cluster {
 // gives avg distance between two clusters
 Real ClusterDistanceAVG(const Cluster& one, 
                         const Cluster& two, 
-                        Real** dist_matrix) {
+                        const ost::TriMatrix<Real>& dist_matrix) {
   typedef std::vector<uint>::const_iterator MyIter;
   Real sum  = 0.0;
   for (MyIter i = one.members.begin(); i != one.members.end(); ++i) {
     for (MyIter j = two.members.begin(); j != two.members.end(); ++j) {
-      sum += dist_matrix[*i][*j];
+      sum += dist_matrix.Get(*i, *j);
     }
   }
   return sum / (one.members.size() * two.members.size());
@@ -33,13 +33,13 @@ Real ClusterDistanceAVG(const Cluster& one,
 // gives max distance between two clusters
 Real ClusterDistanceMAX(const Cluster& one, 
                         const Cluster& two, 
-                        Real** dist_matrix) {
+                        const ost::TriMatrix<Real>& dist_matrix) {
   typedef std::vector<uint>::const_iterator MyIter;
   Real max_dist = -std::numeric_limits<Real>::max();
   Real current_dist;
   for (MyIter i = one.members.begin(); i != one.members.end(); ++i) {
     for (MyIter j = two.members.begin(); j != two.members.end(); ++j) {
-      current_dist = dist_matrix[*i][*j];
+      current_dist = dist_matrix.Get(*i, *j);
       if (current_dist > max_dist) max_dist = current_dist;
     }
   }
@@ -49,13 +49,13 @@ Real ClusterDistanceMAX(const Cluster& one,
 // gives min distance between two clusters
 Real ClusterDistanceMIN(const Cluster& one, 
                         const Cluster& two, 
-                        Real** dist_matrix) {
+                        const ost::TriMatrix<Real>& dist_matrix) {
   typedef std::vector<uint>::const_iterator MyIter;
   Real min_dist = std::numeric_limits<Real>::max();
   Real current_dist;
   for (MyIter i = one.members.begin(); i != one.members.end(); ++i) {
     for (MyIter j = two.members.begin(); j != two.members.end(); ++j) {
-      current_dist = dist_matrix[*i][*j];
+      current_dist = dist_matrix.Get(*i, *j);
       if (current_dist < min_dist) min_dist = current_dist;
     }
   }
@@ -98,12 +98,15 @@ bool PositiveCompareQueueEntries(const QueueEntry& a, const QueueEntry& b){
 
 namespace promod3 { namespace core {
 
-void DoClustering(Real** dist_matrix, uint num, DistanceMetric metric,
+void DoClustering(const ost::TriMatrix<Real>& dist_matrix, 
+                  DistanceMetric metric,
                   Real thresh, char direction, 
                   std::vector<std::vector<uint> >& output){
 
   output.clear();
 
+  uint num = dist_matrix.GetSize();
+
   // do the stupid cases
   if(num == 0){
     return;
@@ -129,7 +132,7 @@ void DoClustering(Real** dist_matrix, uint num, DistanceMetric metric,
   //let's get a function ptr to the right metric function
   typedef Real (*metric_f_ptr)(const Cluster& one, 
                                const Cluster& two, 
-                               Real** dist_matrix);
+                               const ost::TriMatrix<Real>& dist_matrix);
   metric_f_ptr metric_ptr;
 
   switch (metric) {
@@ -176,7 +179,7 @@ void DoClustering(Real** dist_matrix, uint num, DistanceMetric metric,
 
   for(uint i = 0; i < num; ++i){
     for(uint j = i + 1; j < num; ++j){
-      distance_queue.push(QueueEntry(i, j, dist_matrix[i][j]));
+      distance_queue.push(QueueEntry(i, j, dist_matrix.Get(i, j)));
     }
   }
 
diff --git a/core/src/cluster.hh b/core/src/cluster.hh
index 032ed7886163d5d539dd424354963e30ae94c973..5c86bf02796e6566ba7fbcdbc3fef71ec655a46d 100644
--- a/core/src/cluster.hh
+++ b/core/src/cluster.hh
@@ -4,6 +4,7 @@
 #include <vector>
 #include <limits>
 #include <ost/base.hh>
+#include <ost/tri_matrix.hh>
 
 namespace promod3 { namespace core {
 
@@ -14,9 +15,9 @@ enum DistanceMetric {
 };
   
 
-void DoClustering(Real** dist_matrix, uint num, DistanceMetric metric,
-	                Real thresh, char direction, 
-	                std::vector<std::vector<uint> >& output);
+void DoClustering(const ost::TriMatrix<Real>& dist_matrix, 
+                  DistanceMetric metric, Real thresh, char direction, 
+                  std::vector<std::vector<uint> >& output);
 
 }} //ns
 
diff --git a/core/src/superpose.cc b/core/src/superpose.cc
index 7eb504029aace3c2305ae4a438afccbcf0e1e55a..597fa69c5b814409e923ce3a60b7e1f8bfd0cd2c 100644
--- a/core/src/superpose.cc
+++ b/core/src/superpose.cc
@@ -412,7 +412,7 @@ std::pair<geom::Mat4,Real> Superposition(EMatX3& pos_one, EMatX3& pos_two,
   return std::make_pair(transform, rmsd); 
 }
 
-void FillRMSDMatrix(EMatX3List& position_list, Real** data){
+void FillRMSDMatrix(EMatX3List& position_list, ost::TriMatrix<Real>& data){
 
   if(position_list.empty()){
     throw promod3::Error("Empty position list in FillRMSDMatrix!");
@@ -422,6 +422,11 @@ void FillRMSDMatrix(EMatX3List& position_list, Real** data){
   int num_pos = position_list.size();
   int pos_size = position_list[0].rows(); 
 
+  if(num_pos != data.GetSize()) {
+    throw promod3::Error("Number of positions and size of rmsd matrix is "
+                         "inconsistent in FillRMSDMatrix function!");
+  }
+
   if(pos_size < 3) {
     throw promod3::Error("Observed superposition with < 3 positions to "
                          "superpose!");
@@ -442,11 +447,10 @@ void FillRMSDMatrix(EMatX3List& position_list, Real** data){
   Real rmsd;
   EMat3 rot;
   for(int i = 0; i < num_pos; ++i){
-    data[i][i] = 0.0;
+    data.Set(i,i,Real(0.0));
     for(int j = i + 1; j < num_pos; ++j){
       TheobaldRMSD(position_list[i], position_list[j], rmsd, rot);
-      data[i][j] = rmsd;
-      data[j][i] = rmsd;
+      data.Set(i, j, rmsd);
     }
   }
 }
diff --git a/core/src/superpose.hh b/core/src/superpose.hh
index 8a169a784d23389fe8244108e22e905f12ecd367..8e9839c9575695d1fdc8520250cba9f9b02223e1 100644
--- a/core/src/superpose.hh
+++ b/core/src/superpose.hh
@@ -3,6 +3,7 @@
 
 #include <ost/geom/vecmat3_op.hh>
 #include <ost/geom/mat4.hh>
+#include <ost/tri_matrix.hh>
 #include <promod3/core/message.hh>
 
 #include <Eigen/SVD>
@@ -43,7 +44,7 @@ std::pair<geom::Mat4,Real> Superposition(EMatX3& pos_one, EMatX3& pos_two,
 
 // Fill matrix with pairwise RMSD values
 // The positions you pass will be centered, their mean will be the origin...
-void FillRMSDMatrix(EMatX3List& position_list, Real** data);
+void FillRMSDMatrix(EMatX3List& position_list, ost::TriMatrix<Real>& data);
 
 // ITERATIVE MIN RMSD SUPERPOSITION ALGORITHMS
 // You'll always know what indices are actually used for the final
diff --git a/modelling/src/loop_candidate.cc b/modelling/src/loop_candidate.cc
index 49885a31c854a2fa4edec7558196c8ac21f07ac9..82fda8923c1cab9c909f45903756889f6de2f577 100644
--- a/modelling/src/loop_candidate.cc
+++ b/modelling/src/loop_candidate.cc
@@ -288,10 +288,7 @@ LoopCandidates::GetClusters(Real max_dist,
     clusters[0].push_back(0);
   } else if (num_loops > 0) {
 
-    Real** rmsd_matrix = new Real*[num_loops];
-    for (uint i= 0; i < num_loops; ++i) {
-      rmsd_matrix[i] = new Real[num_loops];
-    }
+    ost::TriMatrix<Real> rmsd_matrix(num_loops, Real(0.0));
     
     if(superposed_rmsd){
       // prepare the Eigen positions
@@ -311,24 +308,16 @@ LoopCandidates::GetClusters(Real max_dist,
     }
     else{
       for(uint i = 0; i < this->size(); ++i){
-        rmsd_matrix[i][i] = 0.0;
         for(uint j = i + 1; j < this->size(); ++j){
           Real ca_rmsd = candidates_[i].CARMSD(candidates_[j], false);
-          rmsd_matrix[i][j] = ca_rmsd;
-          rmsd_matrix[j][i] = ca_rmsd;
+          rmsd_matrix.Set(i, j, ca_rmsd);
         }
       }
     }
 
     // get clusters
-    promod3::core::DoClustering(rmsd_matrix, num_loops,
-                                promod3::core::AVG_DIST_METRIC,
+    promod3::core::DoClustering(rmsd_matrix, promod3::core::AVG_DIST_METRIC,
                                 max_dist, '-', clusters);
-    // clean up
-    for (uint i = 0; i < num_loops; ++i) {
-      delete [] rmsd_matrix[i];
-    }
-    delete [] rmsd_matrix;
   }
 }
 
@@ -368,10 +357,7 @@ LoopCandidatesPtr LoopCandidates::GetLargestCluster(Real max_dist,
     return Extract(indices);
   } else {
 
-    Real** rmsd_matrix = new Real*[num_loops];
-    for (uint i= 0; i < num_loops; ++i) {
-      rmsd_matrix[i] = new Real[num_loops];
-    }
+    ost::TriMatrix<Real> rmsd_matrix(num_loops, Real(0.0));
 
     if(superposed_rmsd){
       // prepare the Eigen positions
@@ -391,11 +377,9 @@ LoopCandidatesPtr LoopCandidates::GetLargestCluster(Real max_dist,
     }
     else{
       for(uint i = 0; i < this->size(); ++i){
-        rmsd_matrix[i][i] = 0.0;
         for(uint j = i + 1; j < this->size(); ++j){
           Real ca_rmsd = candidates_[i].CARMSD(candidates_[j], false);
-          rmsd_matrix[i][j] = ca_rmsd;
-          rmsd_matrix[j][i] = ca_rmsd;
+          rmsd_matrix.Set(i, j, ca_rmsd);
         }
       }
     }
@@ -404,11 +388,10 @@ LoopCandidatesPtr LoopCandidates::GetLargestCluster(Real max_dist,
     int max_num = 0;
     int max_idx = 0;
 
-    // going triangular might be more efficient
     for(uint i = 0; i < num_loops; ++i){
       int num = 0;
       for(uint j = 0; j < num_loops; ++j){
-        if(rmsd_matrix[i][j] < max_dist){
+        if(rmsd_matrix.Get(i, j) < max_dist) {
           ++num;
         }
       }
@@ -420,17 +403,11 @@ LoopCandidatesPtr LoopCandidates::GetLargestCluster(Real max_dist,
 
     std::vector<uint> indices;
     for(uint i = 0; i < num_loops; ++i){
-      if(rmsd_matrix[max_idx][i] < max_dist){
+      if(rmsd_matrix.Get(max_idx, i) < max_dist) {
         indices.push_back(i);
       }
     }
 
-    // clean up
-    for (uint i = 0; i < num_loops; ++i) {
-      delete [] rmsd_matrix[i];
-    }
-    delete [] rmsd_matrix;
-
     return Extract(indices);
   }
 }
diff --git a/modelling/src/rigid_blocks.cc b/modelling/src/rigid_blocks.cc
index 4a4fb3eccbec6809d74db5385d4337017d2a6a6d..ef24abeb7a846f4039e1fd5331a5fc13b954691b 100644
--- a/modelling/src/rigid_blocks.cc
+++ b/modelling/src/rigid_blocks.cc
@@ -51,10 +51,7 @@ void DoIt(promod3::core::EMatX3& pos_mat_one,
                              initial_transformations);
 
   uint num_blocks = initial_indices.size();
-  Real** dist_mat = new Real*[num_blocks];
-  for(uint i= 0; i < num_blocks; ++i){
-    dist_mat[i] = new Real[num_blocks];
-  }
+  ost::TriMatrix<Real> dist_mat(num_blocks, Real(0.0));
 
   //calculate the pairwise distances
   for(uint i = 0; i < num_blocks; ++i){
@@ -62,22 +59,14 @@ void DoIt(promod3::core::EMatX3& pos_mat_one,
       Real dist = IndexListSimilarity(initial_indices[i],
                                       initial_indices[j],
                                       num_rows);
-      dist_mat[i][j] = dist;
-      dist_mat[j][i] = dist;
+      dist_mat.Set(i, j, dist);
     }
   }  
 
   std::vector<std::vector<uint> > clusters;
-  promod3::core::DoClustering(dist_mat, num_blocks,
-                              promod3::core::MIN_DIST_METRIC,
+  promod3::core::DoClustering(dist_mat, promod3::core::MIN_DIST_METRIC,
                               cluster_thresh, '+', clusters);
 
-  for(uint i = 0; i < num_blocks; ++i){
-    delete [] dist_mat[i];
-  }
-  delete [] dist_mat;
-
-
   //map over the final indices
   indices.clear();
   std::vector<uint>::iterator it;
diff --git a/scoring/src/constraint_constructor.cc b/scoring/src/constraint_constructor.cc
index dec98385908c87acd1ceed41d2c4174db6629582..15e200177219085de069322553b78a6615bd9ba3 100644
--- a/scoring/src/constraint_constructor.cc
+++ b/scoring/src/constraint_constructor.cc
@@ -24,32 +24,22 @@ void ClusterSequences(const ost::seq::AlignmentHandle& aln,
   //clusters all sequences, neglecting the first
   int num_sequences = aln.GetCount() - 1;
 
-  //pairwise distances come in here
-  Real** dist_matrix = new Real*[num_sequences];
-  for(int i = 0; i < num_sequences; ++i){
-    dist_matrix[i] = new Real[num_sequences];
-    memset(dist_matrix[i], 0, num_sequences*sizeof(Real));
-  }
+  ost::TriMatrix<Real> dist_matrix(num_sequences, Real(0.0));
 
   //sequence similarity serves as distance for clustering
   for(int i = 0; i < num_sequences; ++i){
     for(int j = i+1; j < num_sequences; ++j){
       Real seq_sim = ost::seq::alg::SequenceSimilarity(aln, subst, true, 
                                                        i+1, j+1);
-      dist_matrix[i][j] = seq_sim;
-      dist_matrix[j][i] = seq_sim;
+      dist_matrix.Set(i,j, seq_sim);
     }
   }
 
   //lets cluster the stuff
-  promod3::core::DoClustering(dist_matrix, num_sequences, metric,
+  promod3::core::DoClustering(dist_matrix, metric,
                               cluster_thresh, '+', clusters); 
 
-  //clean up
-  for(int i = 0; i < num_sequences; ++i){
-    delete [] dist_matrix[i];
-  }
-  delete [] dist_matrix;
+
 }
 
 void GenerateClusterWeights(const ost::seq::AlignmentHandle& aln,