diff --git a/lib/nBodySim.cpp b/lib/nBodySim.cpp
index 1701e28f15c06a3de2c810464990031203e1fb61..179a270daa45609b4e4ea0f324ff6b3edda9d530 100644
--- a/lib/nBodySim.cpp
+++ b/lib/nBodySim.cpp
@@ -81,7 +81,7 @@ void nBodySim::runSimulation(double dt, unsigned nSteps) {
 
         std::cout << std::right << "Step "    << std::setw(5)  << i 
                                 // << ": MFM = " << std::setw(12) << calculateMeanForceMagnitude() 
-                                // << ", DEV = " << std::setw(12) << verifyForces() 
+                                << ", DEV = " << std::setw(12) << verifyForces() 
                                 // << ", MCD = " << std::setw(8)  << calculateMeanCenterDistance()
                                 << ", MID = " << std::setw(8)  << calculateMeanInterparticleDistance() 
                                 << ", COM = " << "(" << COM[0] << ", " << COM[1] << ", " << COM[2] << ")"
@@ -342,4 +342,10 @@ void nBodySim::saveState2file(unsigned step, std::ofstream& file) const {
     for (unsigned i=0; i < nParticles_; ++i) {
         file << potential_[i] << std::endl;
     }
+}
+
+
+void nBodySim::saveTree2file(std::string filename) const {
+    std::ofstream file(filename);
+    tree_->save2file(file);
 }
\ No newline at end of file
diff --git a/lib/nBodySim.hpp b/lib/nBodySim.hpp
index daef9131bb22716591ff753e2aa0aab61ea9d230..4b6be7dc9ba9416045eb5966d43a5e956297a489 100644
--- a/lib/nBodySim.hpp
+++ b/lib/nBodySim.hpp
@@ -43,6 +43,8 @@ public:
     void write2file(const double* array, std::ofstream& file, unsigned dim, unsigned skip) const;
     // Write current state of simulation to file
     void saveState2file(unsigned step, std::ofstream& file) const;
+    // save location and size of all nodes in tree
+    void saveTree2file(std::string filename) const;
 
     // Functions to get private variables
     double getMass() { return mass_; }
diff --git a/lib/node.cpp b/lib/node.cpp
index 22ebc3c80131d615cf35782110419b8e80e7c30d..5f75e6ec8b51dfe4209e45635d7cb47f57fddd58 100644
--- a/lib/node.cpp
+++ b/lib/node.cpp
@@ -2,7 +2,7 @@
 
 #include <vector>
 #include <cmath>
-#include <iostream>
+#include <cassert>
 
 Node::Node(double mass, double* positions, unsigned nParticles, double softening, unsigned depth, 
         double size, Vector center, List localParticles, unsigned nLocalParticles) {
@@ -50,8 +50,8 @@ void Node::update(List localParticles, unsigned nLocalParticles) {
     for (unsigned i=0; i<8; ++i) {
         childParticles_[i] = List{};
     }
-    for (unsigned i=0; i<nLocalParticles; ++i) {
-        childParticles_[getOctant(localParticles[i])].push_back(i);
+    for (unsigned i : localParticles) {
+        childParticles_[getOctant(i)].push_back(i);
     }
 
     updateMassAndCOM();
@@ -80,13 +80,13 @@ void Node::update(List localParticles, unsigned nLocalParticles) {
             }
         } else {
             // if child wouldn't contain enough particles, don't create a new child
-            if (childParticles_[i].size() < minNParticles_) continue;
+            if (childParticles_[i].size() <= minNParticles_) continue;
             
             // if child is a nullptr and particleDivision[i] contains enough particles, create a new child
             double childSize = size_/2;
             Vector childCenter(3);
             for (unsigned j=0; j<3; ++j) {
-                childCenter[j] = center_[j] + (i & (1 << j) ? childSize : -childSize);
+                childCenter[j] = center_[j] + (i & (1 << (2-j)) ? childSize : -childSize);
             }
             children_[i] = new Node(mass_, positions_, nParticles_, softening_, depth_+1, 
                     childSize, childCenter, childParticles_[i], childParticles_[i].size());
@@ -258,4 +258,13 @@ unsigned Node::getOctant(unsigned particle) {
 
 Vector Node::getCOM() const {
     return centerOfMass_;
+}
+
+void Node::save2file(std::ofstream& file) const {
+    file << center_[0] << " " << center_[1] << " " << center_[2] << " " << size_ << " " << nLocalParticles_ << std::endl;
+    for (unsigned i=0; i<8; ++i) {
+        if (children_[i] != nullptr) {
+            children_[i]->save2file(file);
+        }
+    }
 }
\ No newline at end of file
diff --git a/lib/node.hpp b/lib/node.hpp
index b1ab8c2802bafde0a7a6df882bd673865e5fba83..1b4e732882150e55ed7c1de9662b918da5776847 100644
--- a/lib/node.hpp
+++ b/lib/node.hpp
@@ -2,6 +2,7 @@
 #define NODE_HPP
 
 #include <vector>
+#include <fstream>
 
 using Vector = std::vector<double>;
 using Matrix = std::vector<Vector>;
@@ -34,6 +35,8 @@ public:
     unsigned getOctant(unsigned particle);
     // get center of mass
     Vector getCOM() const;
+    // save node to file
+    void save2file(std::ofstream& file) const;
 
 private:
     // Node* root_;
diff --git a/lib/tree.cpp b/lib/tree.cpp
index 04bd07838772f79721b30b8677ac77abe1a3b724..c9b84e0e8a592d8c4409082e88a7ae4458da6de1 100644
--- a/lib/tree.cpp
+++ b/lib/tree.cpp
@@ -35,4 +35,8 @@ void Tree::update() {
 
 Vector Tree::getCOM() const {
     return root_->getCOM();
+}
+
+void Tree::save2file(std::ofstream& file) const {
+    root_->save2file(file);
 }
\ No newline at end of file
diff --git a/lib/tree.hpp b/lib/tree.hpp
index 90cf1a8e21f840170be4a77e80af88f32c4a35d8..dc034c029bdf1633fcd736b9f82f6fd6bd552219 100644
--- a/lib/tree.hpp
+++ b/lib/tree.hpp
@@ -4,6 +4,7 @@
 #include "node.hpp"
 
 #include <vector>
+#include <fstream>
 
 class Tree {
 public:
@@ -18,6 +19,8 @@ public:
     void update();
     // get center of mass
     Vector getCOM() const;
+    // save tree to file
+    void save2file(std::ofstream& file) const;
 
 private:
     Node* root_;
diff --git a/main.cpp b/main.cpp
index 0df677a2cb368fe1f0cf80493d3d9fbcfd000efa..0ab77d565e4b97e1512e1e9aeaa087242ffabbff 100644
--- a/main.cpp
+++ b/main.cpp
@@ -28,8 +28,11 @@ int main(int argc, char** argv) {
     // double* positions = sim.getPositions();
     unsigned N = sim.getNParticles();
     std::cout << "N = " << N << std::endl;
+    
+    sim.saveTree2file("../out/tree.txt");
 
     try {
+        std::cout << "Running simulation..." << std::endl;
         auto start = std::chrono::high_resolution_clock::now();
         sim.runSimulation(dt, nSteps);
         auto stop = std::chrono::high_resolution_clock::now();