diff --git a/cpp/field.cpp b/cpp/field.cpp
index c4cb192a8a543733d015d43c793104c7b5969868..b92fce63ecc9ddedd7a093a3c68c410604d2159c 100644
--- a/cpp/field.cpp
+++ b/cpp/field.cpp
@@ -1455,9 +1455,9 @@ void field<rnumber, be, fc>::Hermitian_reflect()
 template <typename rnumber,
           field_backend be,
           field_components fc>
-void field<rnumber, be, fc>::symmetrize()
+void field<rnumber, be, fc>::symmetrize_FFTW()
 {
-    TIMEZONE("field::symmetrize");
+    TIMEZONE("field::symmetrize_FFTW");
     assert(!this->real_space_representation);
     typename fftw_interface<rnumber>::complex *cdata = this->get_cdata();
 
@@ -1506,6 +1506,57 @@ void field<rnumber, be, fc>::symmetrize()
     this->dft();
     this->normalize();
     return;
+}
+
+template <typename rnumber,
+          field_backend be,
+          field_components fc>
+void field<rnumber, be, fc>::symmetrize()
+{
+    TIMEZONE("field::symmetrize");
+    assert(!this->real_space_representation);
+    typename fftw_interface<rnumber>::complex *cdata = this->get_cdata();
+
+    // make 0 mode real
+    if (this->myrank == this->clayout->rank[0][0])
+    {
+        for (ptrdiff_t cc = 0; cc < ncomp(fc); cc++)
+            cdata[cc][1] = 0.0;
+    }
+    // put kx = nx/2 modes to 0
+    for (ptrdiff_t iy = 0; iy < ptrdiff_t(this->clayout->subsizes[0]); iy++)
+    for (ptrdiff_t iz = 0; iz < ptrdiff_t(this->clayout->subsizes[1]); iz++)
+    {
+        ptrdiff_t cindex = this->get_cindex(this->clayout->sizes[2]-1, iy, iz);
+        for (int cc = 0; cc < int(ncomp(fc)); cc++) {
+            (*(cdata + cc + ncomp(fc)*cindex))[0] = 0.0;
+            (*(cdata + cc + ncomp(fc)*cindex))[1] = 0.0;
+        }
+    }
+    // put ky = ny/2 modes to 0
+    if (this->clayout->myrank == this->clayout->rank[0][this->clayout->sizes[0]/2])
+    {
+        for (ptrdiff_t iz = 0; iz < ptrdiff_t(this->clayout->subsizes[1]); iz++)
+        for (ptrdiff_t ix = 0; ix < ptrdiff_t(this->clayout->subsizes[2]); ix++)
+        {
+            ptrdiff_t cindex = this->get_cindex(ix, this->clayout->sizes[0]/2-this->clayout->starts[0], iz);
+            for (int cc = 0; cc < int(ncomp(fc)); cc++) {
+                (*(cdata + cc + ncomp(fc)*cindex))[0] = 0.0;
+                (*(cdata + cc + ncomp(fc)*cindex))[1] = 0.0;
+            }
+        }
+    }
+    // put kz = nz/2 modes to 0
+    for (ptrdiff_t iy = 0; iy < ptrdiff_t(this->clayout->subsizes[0]); iy++)
+    for (ptrdiff_t ix = 0; ix < ptrdiff_t(this->clayout->subsizes[2]); ix++)
+    {
+        ptrdiff_t cindex = this->get_cindex(ix, iy, this->clayout->sizes[1]/2);
+        for (int cc = 0; cc < int(ncomp(fc)); cc++) {
+            (*(cdata + cc + ncomp(fc)*cindex))[0] = 0.0;
+            (*(cdata + cc + ncomp(fc)*cindex))[1] = 0.0;
+        }
+    }
+
     // symmetrize kx = 0 plane, line by line, for ky != 0
     MPI_Status *mpistatus = new MPI_Status;
     // bufferp will hold data from "plus", i.e. iy
diff --git a/cpp/field.hpp b/cpp/field.hpp
index cd7a06d74fd4bc670920121976bd28008b40003c..f9c88bcdafd5b0bfc7f7b954b5e8980086a7ac6d 100644
--- a/cpp/field.hpp
+++ b/cpp/field.hpp
@@ -122,6 +122,7 @@ class field
         void ift();
         void normalize();
         void symmetrize();
+        void symmetrize_FFTW();
         void Hermitian_reflect();
 
         /* stats */