diff --git a/src/MD.cpp b/src/MD.cpp
index b060820..0d615cc 100644
--- a/src/MD.cpp
+++ b/src/MD.cpp
@@ -592,9 +592,10 @@ double PotentialAndAcceleration(double dt) {
 	double Pot = 0.;
 
 	for (int i = 0; i < N - 1; i++) {
+		double ai0 = 0., ai1 = 0., ai2 = 0.;
 		double rI[3] = {r[i][0], r[i][1], r[i][2]};
 		for (int j = i + 1; j < N; j++) {
-			double quot, rnorm, term1, term2;
+			double quot, term2;
 			//  component-by-componenent position of i relative to j
 			double rij[3];
 			//  sum of squares of the components
@@ -613,18 +614,25 @@ double PotentialAndAcceleration(double dt) {
 			double rSqd_7 = rSqd_3 * rSqd_3 * rSqd;
 			double f = (48. - (24. * rSqd_3)) / rSqd_7;
 			// from F = ma, where m = 1 in natural units!
-			for (int k = 0; k < 3; k++) {
-				double tmp = rij[k] * f;
-				a[i][k] += tmp;
-				a[j][k] -= tmp;
-			}
+			double tmp = rij[0] * f;
+			double tmp1 = rij[1] * f;
+			double tmp2 = rij[2] * f;
+			ai0 += tmp;
+			ai1 += tmp1;
+			ai2 += tmp2;
+			a[j][0] -= tmp;
+			a[j][1] -= tmp1;
+			a[j][2] -= tmp2;
 
 			quot = sigma * sigma / rSqd;
 			term2 = quot * quot * quot;
-			Pot += epsilon_8 * term2 * (term2 - 1.);
+			Pot += term2 * (term2 - 1.);
 		}
+		a[i][0] += ai0;
+		a[i][1] += ai1;
+		a[i][2] += ai2;
 	}
-	return Pot;
+	return Pot*epsilon_8;
 }
 
 // Function to calculate the potential energy of the system