diff --git a/lib/nBodySim.cpp b/lib/nBodySim.cpp
index adb29808714655dd5a96bfaf4e461f735a4583e6..563cecd4bfe85661e0746f0b4cc174d1a7241d20 100644
--- a/lib/nBodySim.cpp
+++ b/lib/nBodySim.cpp
@@ -58,45 +58,55 @@ nBodySim::~nBodySim() {
 void nBodySim::runSimulation(double dt, unsigned nSteps) {
     for (unsigned i=0; i < nSteps; ++i) {
         calculateForces();
-        updatePositions(dt);
+        doTimeStep(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
+    #pragma omp parallel for schedule(guided)
     for (unsigned i = 0; i < nParticles_; i++) {
+        double mi = masses_[i];
+        double s = softening_[i];
+        double px = positions_[3*i];
+        double py = positions_[3*i+1];
+        double pz = positions_[3*i+2];
+        double fx = 0.0;
+        double fy = 0.0;
+        double fz = 0.0;
         for (unsigned j = i+1; j < nParticles_; j++) {
+            // if (j <= i) continue;
             // 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 dx = positions_[3*j]   - px;
+            double dy = positions_[3*j+1] - py;
+            double dz = positions_[3*j+2] - pz;
             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];
+            // calculate forces
             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;
+            double fxij = mi*mj * dx / (r2 + s*s);
+            double fyij = mi*mj * dy / (r2 + s*s);
+            double fzij = mi*mj * dz / (r2 + s*s);
+            fx += fxij;
+            fy += fyij;
+            fz += fzij;
+            forces_[3*j]   -= fxij;
+            forces_[3*j+1] -= fyij;
+            forces_[3*j+2] -= fzij;
         }
+        forces_[3*i]   += fx;
+        forces_[3*i+1] += fy;
+        forces_[3*i+2] += fz;
     }
 }
-//*/
+
 /*
 void nBodySim::calculateForces() {
     #pragma omp parallel for
@@ -129,7 +139,7 @@ void nBodySim::calculateForces() {
 }
 */
 
-void nBodySim::updatePositions(double dt) {
+void nBodySim::doTimeStep(double dt) {
     #pragma omp parallel for
     for (unsigned i=0; i < nParticles_; ++i) {
         for (unsigned j=0; j < 3; ++j) {
diff --git a/lib/nBodySim.hpp b/lib/nBodySim.hpp
index df1fa4ef1d7cc57c34fdec9774bb45bb960db45c..d296f15dd8c623d608f056afd96af31ab631c518 100644
--- a/lib/nBodySim.hpp
+++ b/lib/nBodySim.hpp
@@ -14,7 +14,7 @@ public:
     // loops over all pairs of particles and calculates forces
     void calculateForces();
     // updates positions and velocities using the calculated forces
-    void updatePositions(double dt);
+    void doTimeStep(double dt);
     // writes data to file, every row is a particle
     void write2file(const double* array, std::string filename, unsigned dim) const;