@@ -1612,6 +1612,7 @@ class colvar_grid_gradient : public colvar_grid<cvm::real>
16121612 std::vector<int > bin (nd, 0 );
16131613 // We process points where exp(-d^2 / 2sigma^2) > 10^-3
16141614 // This implies distance < 3.72 * sigma
1615+ cvm::real inv_squared_smooth = 1 / (smoothing*smoothing);
16151616 cvm::real const cutoff_factor = 3.72 ;
16161617 // TODO: make sure that this is not > min nx /2
16171618 cvm::real const cutoff = cutoff_factor * smoothing; // take like floor()
@@ -1621,8 +1622,8 @@ class colvar_grid_gradient : public colvar_grid<cvm::real>
16211622 std::vector<int >ix_max (nd, 0 );
16221623 std::vector<int >periodic_offset (nd,0 );
16231624 for (size_t i =0 ; i < nd; i++) {
1624- ixmin = cv_value[0 ] - cutoff;
1625- ixmax = cv_value[0 ] + cutoff;
1625+ ixmin = cv_value[i ] - cutoff;
1626+ ixmax = cv_value[i ] + cutoff;
16261627 int cut_ixmin = static_cast <int >(cvm::floor (std::max (ixmin,static_cast <cvm::real>(0 ))));
16271628 int cut_ixmax = static_cast <int >(cvm::floor (std::min (ixmax, static_cast <cvm::real>(nx[0 ]) - 1 )));
16281629 if (periodic[i]) {
@@ -1645,7 +1646,7 @@ class colvar_grid_gradient : public colvar_grid<cvm::real>
16451646 }
16461647 dist += std::pow (ix_copy[dim] +0.5 - cv_value[dim], 2 );
16471648 }
1648- acc_force (ix_copy, force, true , cvm::exp (-dist/smoothing ));
1649+ acc_force (ix_copy, force, true , cvm::exp (-dist * inv_squared_smooth ));
16491650 }
16501651 } else {
16511652 // cvm::log("we update the force with acc_force, force = " + cvm::to_str(force[0]) + " " + cvm::to_str(force[1]));
0 commit comments