diff --git a/lib/nBodySim.cpp b/lib/nBodySim.cpp
index b2eca9d3f01ce0f5294ee822529e0344429bc97f..adb29808714655dd5a96bfaf4e461f735a4583e6 100644
--- a/lib/nBodySim.cpp
+++ b/lib/nBodySim.cpp
@@ -61,7 +61,43 @@ void nBodySim::runSimulation(double dt, unsigned nSteps) {
         updatePositions(dt);
     }
 }
-
+///*
+void nBodySim::calculateForces() {
+    #pragma omp parallel for
+    for (unsigned i=0; i < 3*nParticles_; ++i) {
+        forces_[i] = 0.0;
+    }
+    // loop over each pair of particles and calculate the force between them
+    #pragma omp parallel for
+    for (unsigned i = 0; i < nParticles_; i++) {
+        for (unsigned j = i+1; j < nParticles_; j++) {
+            // calculate distance between particles
+            double dx = positions_[3*j]   - positions_[3*i];
+            double dy = positions_[3*j+1] - positions_[3*i+1];
+            double dz = positions_[3*j+2] - positions_[3*i+2];
+            double r2 = dx*dx + dy*dy + dz*dz;
+            double r = std::sqrt(r2);
+            dx /= r;
+            dy /= r;
+            dz /= r;
+            // calculate force
+            double mi = masses_[i];
+            double mj = masses_[j];
+            double s = softening_[i];
+            double fx = mi*mj * dx / (r2 + s*s);
+            double fy = mi*mj * dy / (r2 + s*s);
+            double fz = mi*mj * dz / (r2 + s*s);
+            forces_[3*i]   += fx;
+            forces_[3*i+1] += fy;
+            forces_[3*i+2] += fz;
+            forces_[3*j]   -= fx;
+            forces_[3*j+1] -= fy;
+            forces_[3*j+2] -= fz;
+        }
+    }
+}
+//*/
+/*
 void nBodySim::calculateForces() {
     #pragma omp parallel for
     for (unsigned i=0; i < 3*nParticles_; ++i) {
@@ -91,6 +127,7 @@ void nBodySim::calculateForces() {
         }
     }
 }
+*/
 
 void nBodySim::updatePositions(double dt) {
     #pragma omp parallel for