diff --git a/src/cltools/pesmd.cpp b/src/cltools/pesmd.cpp
index 9c2bbb2056db7efc4c344908db0fbb44a60e09b9..a7ef9b91715e53a6cce311e7fc2e01ad318e949c 100644
--- a/src/cltools/pesmd.cpp
+++ b/src/cltools/pesmd.cpp
@@ -212,7 +212,7 @@ public:
     double tke=0;
     for(int i=0; i<nat; ++i) {
       for(int j=0; j<3; ++j) {
-        if( 3*i+j>dim ) break;
+        if( 3*i+j>dim-1 ) break;
         tke += 0.5*velocities[i][j]*velocities[i][j];
       }
     }
@@ -241,7 +241,7 @@ public:
       double lrand=sqrt((1.-lscale*lscale)*temp);
       for(int j=0; j<nat; ++j) {
         for(int k=0; k<3; ++k) {
-          if( 3*j+k>dim ) break;
+          if( 3*j+k>dim-1 ) break;
           therm_eng=therm_eng+0.5*velocities[j][k]*velocities[j][k];
           velocities[j][k]=lscale*velocities[j][k]+lrand*random.Gaussian();
           therm_eng=therm_eng-0.5*velocities[j][k]*velocities[0][k];
@@ -251,7 +251,7 @@ public:
       // First step of velocity verlet
       for(int j=0; j<nat; ++j) {
         for(int k=0; k<3; ++k) {
-          if( 3*j+k>dim ) break;
+          if( 3*j+k>dim-1 ) break;
           velocities[j][k] = velocities[j][k] + 0.5*tstep*forces[1+j][k];
           positions[1+j][k] = positions[1+j][k] + tstep*velocities[j][k];
         }
@@ -274,7 +274,7 @@ public:
       // Second step of velocity verlet
       for(int j=0; j<nat; ++j) {
         for(int k=0; k<3; ++k) {
-          if( 3*j+k>dim ) break;
+          if( 3*j+k>dim-1 ) break;
           velocities[j][k] = velocities[j][k] + 0.5*tstep*forces[1+j][k];
         }
       }
@@ -284,7 +284,7 @@ public:
       lrand=sqrt((1.-lscale*lscale)*temp);
       for(int j=0; j<nat; ++j) {
         for(int k=0; k<3; ++k) {
-          if( 3*j+k>dim ) break;
+          if( 3*j+k>dim-1 ) break;
           therm_eng=therm_eng+0.5*velocities[j][k]*velocities[j][k];
           velocities[j][k]=lscale*velocities[j][k]+lrand*random.Gaussian();
           therm_eng=therm_eng-0.5*velocities[j][k]*velocities[j][k];
@@ -294,7 +294,7 @@ public:
       tke=0;
       for(int i=0; i<nat; ++i) {
         for(int j=0; j<3; ++j) {
-          if( 3*i+j>dim ) break;
+          if( 3*i+j>dim-1 ) break;
           tke += 0.5*velocities[i][j]*velocities[i][j];
         }
       }