diff --git a/PW/ewald.f90 b/PW/ewald.f90
index 2e0df75c4..7a8ff9bb4 100644
--- a/PW/ewald.f90
+++ b/PW/ewald.f90
@@ -61,7 +61,7 @@ function ewald (alat, nat, ntyp, ityp, zv, at, bg, tau, omega, g, &
   ! counter on polarization
 
   real(DP) :: charge, tpiba2, ewaldg, ewaldr, dtau (3), alpha, &
-       r (3, mxr), r2 (mxr), rmax, rr, upperbound, fact
+       r (3, mxr), r2 (mxr), rmax, rr, upperbound, fact, ds(3)
   ! total ionic charge in the cell
   ! length in reciprocal space
   ! ewald energy computed in reciprocal space
@@ -134,10 +134,10 @@ function ewald (alat, nat, ntyp, ityp, zv, at, bg, tau, omega, g, &
      !
      do na = 1, nat
         do nb = 1, nat
-           do ipol = 1, 3
-              dtau (ipol) = tau (ipol, na) - tau (ipol, nb)
-              dtau (ipol) = dtau(ipol) - anint(dtau(ipol))
-           enddo
+              dtau (:) = tau (:, na) - tau (:, nb)
+              ds(:) = MATMUL( dtau(:), bg(:,:) )
+              ds(:) = ds(:) - anint(ds(:))
+              dtau(:) = MATMUL( at(:,:), ds(:) )
            !
            ! generates nearest-neighbors shells
            !
diff --git a/PW/force_ew.f90 b/PW/force_ew.f90
index 3ee633aca..3473e0235 100644
--- a/PW/force_ew.f90
+++ b/PW/force_ew.f90
@@ -61,7 +61,7 @@ subroutine force_ew (alat, nat, ntyp, ityp, zv, at, bg, tau, &
   ! counter on polarization
 
   real(DP) :: sumnb, arg, tpiba2, alpha, dtau (3), r (3, mxr), &
-       r2 (mxr), erfc, rmax, rr, charge, upperbound, fact
+       r2 (mxr), erfc, rmax, rr, charge, upperbound, fact, ds(3)
   ! auxiliary variable for speed
   ! the argument of the exponential
   ! 2 pi /alat
@@ -141,10 +141,10 @@ subroutine force_ew (alat, nat, ntyp, ityp, zv, at, bg, tau, &
   do na = 1, nat
      do nb = 1, nat
         if (nb.eq.na) goto 50
-        do ipol = 1, 3
-           dtau (ipol) = tau (ipol, na) - tau (ipol, nb)
-           dtau (ipol) = dtau(ipol) - anint(dtau(ipol))
-        enddo
+              dtau (:) = tau (:, na) - tau (:, nb)
+              ds(:) = MATMUL( dtau(:), bg(:,:) )
+              ds(:) = ds(:) - anint(ds(:))
+              dtau(:) = MATMUL( at(:,:), ds(:) )
         !
         ! generates nearest-neighbors shells r(i)=R(i)-dtau(i)
         !
diff --git a/PW/stres_ewa.f90 b/PW/stres_ewa.f90
index 2acd83fbc..d1ae1c37a 100644
--- a/PW/stres_ewa.f90
+++ b/PW/stres_ewa.f90
@@ -55,7 +55,7 @@ subroutine stres_ewa (alat, nat, ntyp, ityp, zv, at, bg, tau, &
   ! counter on atoms
   ! number of R vectors included in r sum
 
-  real(DP) :: charge, arg, tpiba2, dtau (3), alpha, r (3, mxr), &
+  real(DP) :: charge, arg, tpiba2, dtau (3), alpha, r (3, mxr), ds(3), &
        r2 (mxr), rmax, rr, upperbound, fact, fac, g2, g2a, sdewald, sewald
   ! total ionic charge in the cell
   ! the argument of the phase
@@ -139,10 +139,10 @@ subroutine stres_ewa (alat, nat, ntyp, ityp, zv, at, bg, tau, &
      !
      do na = 1, nat
         do nb = 1, nat
-           do l = 1, 3
-              dtau (l) = tau (l, na) - tau (l, nb)
-              dtau (l) = dtau(l) - anint(dtau(l))
-           enddo
+              dtau (:) = tau (:, na) - tau (:, nb)
+              ds(:) = MATMUL( dtau(:), bg(:,:) )
+              ds(:) = ds(:) - anint(ds(:))
+              dtau(:) = MATMUL( at(:,:), ds(:) )
            !
            !     generates nearest-neighbors shells r(i)=R(i)-dtau(i)
            !