diff --git a/cpp/field.cpp b/cpp/field.cpp
index d22aa180eacd8ddd15e16c54a24af27c77b091eb..2cbada82083a709798d49503b0394983ba0d639c 100644
--- a/cpp/field.cpp
+++ b/cpp/field.cpp
@@ -1581,6 +1581,7 @@ void field<rnumber, be, fc>::symmetrize()
             if (this->clayout->myrank == rankp)
             {
                 ptrdiff_t iyy = iy - this->clayout->starts[0];
+                #pragma omp simd
                 for (ptrdiff_t iz = zstart; iz < zend; iz++)
                 {
                     ptrdiff_t cindex = this->get_cindex(0, iyy, iz);
@@ -1594,6 +1595,7 @@ void field<rnumber, be, fc>::symmetrize()
             if (this->clayout->myrank == rankm)
             {
                 ptrdiff_t iyy = (this->clayout->sizes[0] - iy) - this->clayout->starts[0];
+                #pragma omp simd
                 for (ptrdiff_t iz = zstart; iz < zend; iz++)
                 {
                     ptrdiff_t cindex = this->get_cindex(0, iyy, iz);
@@ -1648,11 +1650,11 @@ void field<rnumber, be, fc>::symmetrize()
             if (this->clayout->myrank == rankp)
             {
                 ptrdiff_t iyy = iy - this->clayout->starts[0];
+                #pragma omp simd
                 for (ptrdiff_t iz = zstart; iz < zend; iz++)
                 {
-                    if (iz == 0)
-                        continue;
-                    ptrdiff_t izz = (this->clayout->sizes[1] - iz);
+                    // modulo operation makes sense only for slab decomposition
+                    ptrdiff_t izz = (this->clayout->sizes[1] - iz) % this->clayout->sizes[1];
                     ptrdiff_t cindex = this->get_cindex(0, iyy, iz);
                     for (int cc = 0; cc < int(ncomp(fc)); cc++)
                     {
@@ -1660,26 +1662,16 @@ void field<rnumber, be, fc>::symmetrize()
                         (*(cdata + ncomp(fc)*cindex + cc))[1] =  ((*(bufferp + ncomp(fc)*iz+cc))[1] - (*(bufferm + ncomp(fc)*izz+cc))[1])/2;
                     }
                 }
-                // zero value treated separately because of izz
-                if (omp_get_thread_num() == 0)
-                {
-                    ptrdiff_t cindex = this->get_cindex(0, iyy, 0);
-                    for (int cc = 0; cc < int(ncomp(fc)); cc++)
-                    {
-                        (*(cdata + cc + ncomp(fc)*cindex))[0] =  ((*(bufferp + cc))[0] + (*(bufferm + cc))[0])/2;
-                        (*(cdata + cc + ncomp(fc)*cindex))[1] =  ((*(bufferp + cc))[1] - (*(bufferm + cc))[1])/2;
-                    }
-                }
             }
             // if I my rank is either plus or minus, I should update my slice
             if (this->clayout->myrank == rankm)
             {
                 ptrdiff_t iyy = (this->clayout->sizes[0] - iy) - this->clayout->starts[0];
+                #pragma omp simd
                 for (ptrdiff_t iz = zstart; iz < zend; iz++)
                 {
-                    if (iz == 0)
-                        continue;
-                    ptrdiff_t izz = (this->clayout->sizes[1] - iz);
+                    // modulo operation makes sense only for slab decomposition
+                    ptrdiff_t izz = (this->clayout->sizes[1] - iz) % this->clayout->sizes[1];
                     ptrdiff_t cindex = this->get_cindex(0, iyy, izz);
                     for (int cc = 0; cc < int(ncomp(fc)); cc++)
                     {
@@ -1687,16 +1679,6 @@ void field<rnumber, be, fc>::symmetrize()
                         (*(cdata + ncomp(fc)*cindex + cc))[1] = -((*(bufferp + ncomp(fc)*iz+cc))[1] - (*(bufferm + ncomp(fc)*izz+cc))[1])/2;
                     }
                 }
-                // zero value treated separately because of izz
-                if (omp_get_thread_num() == 0)
-                {
-                    ptrdiff_t cindex = this->get_cindex(0, iyy, 0);
-                    for (int cc = 0; cc < int(ncomp(fc)); cc++)
-                    {
-                        (*(cdata + cc + ncomp(fc)*cindex))[0] =  ((*(bufferp + cc))[0] + (*(bufferm + cc))[0])/2;
-                        (*(cdata + cc + ncomp(fc)*cindex))[1] = -((*(bufferp + cc))[1] - (*(bufferm + cc))[1])/2;
-                    }
-                }
             }
         }// end omp parallel region
         }// end if