diff --git a/lib/nBodySim.cpp b/lib/nBodySim.cpp
index 9088f7967c32aab3ab4adf040c7a1401eb5a9b10..03243f6e3600cb3d81a1f2d071c15ff981cbef40 100644
--- a/lib/nBodySim.cpp
+++ b/lib/nBodySim.cpp
@@ -68,25 +68,27 @@ void nBodySim::runSimulation(double dt, unsigned nSteps) {
     for (unsigned i=0; i < nSteps; ++i) {
 
         write2file(positions_, file, 3, skip);
+        Vector COM = getCenterOfMass();
 
         // Integrate with leapfrog
         auto start = std::chrono::high_resolution_clock::now();
         kick(dt/2.0);
         drift(dt);
-        updateForces();
+        treeUpdateForces();
         kick(dt/2.0);
         auto stop = std::chrono::high_resolution_clock::now();
         auto duration = std::chrono::duration_cast<std::chrono::milliseconds>(stop - start);
 
         std::cout << std::right << "Step "    << std::setw(5)  << i 
-                                << ": MFM = " << std::setw(12) << calculateMeanForceMagnitude() 
+                                // << ": MFM = " << std::setw(12) << calculateMeanForceMagnitude() 
                                 // << ", DEV = " << std::setw(12) << verifyForces() 
-                                << ", MCD = " << std::setw(8)  << calculateMeanCenterDistance()
-                                // << ", MID = " << std::setw(8)  << calculateMeanCenterDistance() 
+                                // << ", MCD = " << std::setw(8)  << calculateMeanCenterDistance()
+                                << ", MID = " << std::setw(8)  << calculateMeanInterparticleDistance() 
+                                << ", COM = " << std::setw(4)  << "(" << COM[0] << ", " << COM[1] << ", " << COM[2] << ") "
                                 << ", E = "   << std::setw(12) << calculateTotalEnergy()
                                 << ", T = "   << std::setw(8)  << duration.count()/1000.0 << " s"
                                 << std::endl;
-        // updateTree();
+        updateTree();
     }
     file.close();
 }
@@ -289,6 +291,22 @@ double nBodySim::calculateTotalEnergy() const {
     return totalEnergy;
 }
 
+Vector nBodySim::getCenterOfMass() const {
+    /*
+    Vector com = {0.0, 0.0, 0.0};
+    for (unsigned i = 0; i < nParticles_; i++) {
+        com[0] += positions_[3*i];
+        com[1] += positions_[3*i+1];
+        com[2] += positions_[3*i+2];
+    }
+    com[0] /= nParticles_;
+    com[1] /= nParticles_;
+    com[2] /= nParticles_;
+    return com;
+    */
+    return tree_->getCOM(); 
+}
+
 void nBodySim::write2file(const double* array, std::ofstream& file, unsigned dim, unsigned skip) const {
     // writes array to file, with N rows and dim columns
     for (unsigned i = 0; i < nParticles_; i+=skip) {
diff --git a/lib/nBodySim.hpp b/lib/nBodySim.hpp
index fdf41202018ec4cc0055e6c17bcab96fa095cc4b..daef9131bb22716591ff753e2aa0aab61ea9d230 100644
--- a/lib/nBodySim.hpp
+++ b/lib/nBodySim.hpp
@@ -36,6 +36,8 @@ public:
     double calculateMeanCenterDistance() const;
     // Loop over all particles and calculate the total energy
     double calculateTotalEnergy() const;
+    // get center of mass
+    Vector getCenterOfMass() const;
 
     // Writes data to file, every row is a particle
     void write2file(const double* array, std::ofstream& file, unsigned dim, unsigned skip) const;
diff --git a/lib/node.cpp b/lib/node.cpp
index b1bf6c67716814e540b4fd4da6e464e754f52649..22ebc3c80131d615cf35782110419b8e80e7c30d 100644
--- a/lib/node.cpp
+++ b/lib/node.cpp
@@ -254,4 +254,8 @@ unsigned Node::getOctant(unsigned particle) {
             }
         }
     }
+}
+
+Vector Node::getCOM() const {
+    return centerOfMass_;
 }
\ No newline at end of file
diff --git a/lib/node.hpp b/lib/node.hpp
index 244e530be40cf68a61f1c8f91b35c6204c109b51..b1ab8c2802bafde0a7a6df882bd673865e5fba83 100644
--- a/lib/node.hpp
+++ b/lib/node.hpp
@@ -32,6 +32,8 @@ public:
     // Function that takes a particle index in the global particles_ array and returns the octant it is in relative to the center fo the current node 
     // this function doesn't check if the particle is actually in the node.
     unsigned getOctant(unsigned particle);
+    // get center of mass
+    Vector getCOM() const;
 
 private:
     // Node* root_;
diff --git a/lib/tree.cpp b/lib/tree.cpp
index 3d20168e7beb376b507a338a8bab7a29f65fe29c..04bd07838772f79721b30b8677ac77abe1a3b724 100644
--- a/lib/tree.cpp
+++ b/lib/tree.cpp
@@ -31,4 +31,8 @@ void Tree::updateForce(unsigned particle) {
 
 void Tree::update() {
     root_->update(allParticles_, nParticles_);
+}
+
+Vector Tree::getCOM() const {
+    return root_->getCOM();
 }
\ No newline at end of file
diff --git a/lib/tree.hpp b/lib/tree.hpp
index 86cac25b74508a1cb47211296c739546c38c772f..90cf1a8e21f840170be4a77e80af88f32c4a35d8 100644
--- a/lib/tree.hpp
+++ b/lib/tree.hpp
@@ -16,6 +16,8 @@ public:
     void updateForce(unsigned particle);
     // Update tree
     void update();
+    // get center of mass
+    Vector getCOM() const;
 
 private:
     Node* root_;