diff --git a/src/eri_mme/eri_mme_lattice_summation.F b/src/eri_mme/eri_mme_lattice_summation.F
index 6d3f26fed5..8eeb360500 100644
--- a/src/eri_mme/eri_mme_lattice_summation.F
+++ b/src/eri_mme/eri_mme_lattice_summation.F
@@ -759,23 +759,23 @@ SUBROUTINE pgf_sum_3c_1d(S_G, RA, RB, RC, zeta, zetb, zetc, a_mm, lgth, R_bounds
          CALL pgf_sum_3c_rspace_1d_generic(S_G, RA, RB, RC, zeta, zetb, zetc, a_mm, lgth, R_bounds_3)
       ELSE
          #:for l_max in range(0,lmax_unroll+1)
-         IF (l_max == ${l_max}$) THEN
-            #:for m_max in range(0,lmax_unroll+1)
-               IF (m_max == ${m_max}$) THEN
-                  #:for n_max in range(0, lmax_unroll+1)
-                     IF (n_max == ${n_max}$) THEN
-                        #:for prop_exp in range(0,2)
-                           IF (prop_exp == ${prop_exp}$) THEN
-                              CALL pgf_sum_3c_rspace_1d_${l_max}$_${m_max}$_${n_max}$_exp_${prop_exp}$ (S_G, RA, RB, RC, &
+            IF (l_max == ${l_max}$) THEN
+               #:for m_max in range(0,lmax_unroll+1)
+                  IF (m_max == ${m_max}$) THEN
+                     #:for n_max in range(0, lmax_unroll+1)
+                        IF (n_max == ${n_max}$) THEN
+                           #:for prop_exp in range(0,2)
+                              IF (prop_exp == ${prop_exp}$) THEN
+                                 CALL pgf_sum_3c_rspace_1d_${l_max}$_${m_max}$_${n_max}$_exp_${prop_exp}$ (S_G, RA, RB, RC, &
                                                                                            zeta, zetb, zetc, a_mm, lgth, R_bounds_3)
-                              RETURN
-                           END IF
-                        #:endfor
-                     END IF
-                  #:endfor
-               END IF
-            #:endfor
-         END IF
+                                 RETURN
+                              END IF
+                           #:endfor
+                        END IF
+                     #:endfor
+                  END IF
+               #:endfor
+            END IF
          #:endfor
       END IF
 
diff --git a/src/ewalds_multipole_sr.fypp b/src/ewalds_multipole_sr.fypp
index b37aae3542..a544ea1c28 100644
--- a/src/ewalds_multipole_sr.fypp
+++ b/src/ewalds_multipole_sr.fypp
@@ -1108,14 +1108,14 @@
          energy = energy + eloc
       #:endif
       #:if store_forces
-      IF (do_forces) THEN
-         forces(1, atom_a) = forces(1, atom_a) - fr(1)
-         forces(2, atom_a) = forces(2, atom_a) - fr(2)
-         forces(3, atom_a) = forces(3, atom_a) - fr(3)
-         forces(1, atom_b) = forces(1, atom_b) + fr(1)
-         forces(2, atom_b) = forces(2, atom_b) + fr(2)
-         forces(3, atom_b) = forces(3, atom_b) + fr(3)
-      END IF
+         IF (do_forces) THEN
+            forces(1, atom_a) = forces(1, atom_a) - fr(1)
+            forces(2, atom_a) = forces(2, atom_a) - fr(2)
+            forces(3, atom_a) = forces(3, atom_a) - fr(3)
+            forces(1, atom_b) = forces(1, atom_b) + fr(1)
+            forces(2, atom_b) = forces(2, atom_b) + fr(2)
+            forces(3, atom_b) = forces(3, atom_b) + fr(3)
+         END IF
       #:endif
       ! Electric fields
       IF (do_efield) THEN
diff --git a/src/grid/dgemm/grid_dgemm_tensor_local.h b/src/grid/dgemm/grid_dgemm_tensor_local.h
index 905e41cf93..03b08836ed 100644
--- a/src/grid/dgemm/grid_dgemm_tensor_local.h
+++ b/src/grid/dgemm/grid_dgemm_tensor_local.h
@@ -46,9 +46,7 @@ static inline void initialize_tensor(struct tensor_ *a, const int dim,
                                      const int *const sizes) {
   if (a == NULL)
     return;
-#if !defined(NDEBUG)
-  memset(a, 0, sizeof(*a));
-#endif
+
   a->dim_ = dim;
   for (int d = 0; d < dim; d++)
     a->size[d] = sizes[d];
diff --git a/src/grid/ref/grid_ref_collint.h b/src/grid/ref/grid_ref_collint.h
index 88aec98b8c..8ae0e8245e 100644
--- a/src/grid/ref/grid_ref_collint.h
+++ b/src/grid/ref/grid_ref_collint.h
@@ -66,7 +66,6 @@ ortho_cx_to_grid(const int lp, const int k1, const int k2, const int j1,
 #if (GRID_DO_COLLOCATE)
     // collocate
     double reg[4] = {0.0, 0.0, 0.0, 0.0};
-#pragma omp simd reduction(+ : reg)
     for (int lxp = 0; lxp <= lp; lxp++) {
       const double p = pol[0][lxp][i + cmax];
       reg[0] += cx[lxp * 4 + 0] * p;
@@ -82,7 +81,6 @@ ortho_cx_to_grid(const int lp, const int k1, const int k2, const int j1,
 #else
     // integrate
     const double reg[4] = {*grid_0, *grid_1, *grid_2, *grid_3};
-#pragma omp simd
     for (int lxp = 0; lxp <= lp; lxp++) {
       const double p = pol[0][lxp][i + cmax];
       cx[lxp * 4 + 0] += reg[0] * p;
diff --git a/src/pw/realspace_grid_types.F b/src/pw/realspace_grid_types.F
index 2e43d56d2f..fc21f14b46 100644
--- a/src/pw/realspace_grid_types.F
+++ b/src/pw/realspace_grid_types.F
@@ -1916,6 +1916,9 @@ SUBROUTINE rs_grid_zero(rs)
       CALL timeset(routineN, handle)
       l(1) = LBOUND(rs%r, 1); l(2) = LBOUND(rs%r, 2); l(3) = LBOUND(rs%r, 3)
       u(1) = UBOUND(rs%r, 1); u(2) = UBOUND(rs%r, 2); u(3) = UBOUND(rs%r, 3)
+!$OMP PARALLEL DO DEFAULT(NONE) COLLAPSE(3) &
+!$OMP             PRIVATE(i,j,k) &
+!$OMP             SHARED(rs,l,u)
       DO k = l(3), u(3)
       DO j = l(2), u(2)
       DO i = l(1), u(1)
@@ -1923,6 +1926,7 @@ SUBROUTINE rs_grid_zero(rs)
       END DO
       END DO
       END DO
+!$OMP END PARALLEL DO
       CALL timestop(handle)
 
    END SUBROUTINE rs_grid_zero