diff --git a/src/examples/CMakeLists.txt b/src/examples/CMakeLists.txt
index e6cec0c77eb..8eab274c8b2 100644
--- a/src/examples/CMakeLists.txt
+++ b/src/examples/CMakeLists.txt
@@ -9,7 +9,7 @@ set(EXAMPLE_SOURCES
     dataloadbal hatom_1d binaryop dielectric hehf 3dharmonic testsolver
     testspectralprop dielectric_external_field tiny h2dynamic newsolver testcomplexfunctionsolver
     helium_exact density_smoothing siam_example ac_corr dirac-hatom
-    derivatives array_worldobject writecoeff writecoeff2)
+    derivatives array_worldobject)
  
 if(LIBXC_FOUND)
   list(APPEND EXAMPLE_SOURCES hefxc)
diff --git a/src/examples/writecoeffs/CMakeLists.txt b/src/examples/writecoeffs/CMakeLists.txt
index 6bdc4274792..c3d28d472c3 100644
--- a/src/examples/writecoeffs/CMakeLists.txt
+++ b/src/examples/writecoeffs/CMakeLists.txt
@@ -1,7 +1,7 @@
 # src/examples
 
 set(EXAMPLE_SOURCES
-    writecoeff writecoeff2 writecoeff3 h2_write h2_write_json h2_write_json_2 writecoeff_hdf5)
+    writecoeff writecoeff2 writecoeff3 writecoeff_json h2_write h2_write_json h2_write_json_2 writecoeff_hdf5)
 
 # Create executables for example applications
 foreach(example ${EXAMPLE_SOURCES})
diff --git a/src/examples/writecoeffs/FunctionIO2.h b/src/examples/writecoeffs/FunctionIO2.h
index 9bcd61e14a2..495383d53ea 100644
--- a/src/examples/writecoeffs/FunctionIO2.h
+++ b/src/examples/writecoeffs/FunctionIO2.h
@@ -188,6 +188,7 @@ template <typename T, std::size_t NDIM> struct FunctionIOData {
     f.get_impl()->world.gop.fence();
   }
 
+  // write function coordinates for dim 1, 2, 3
   void initialize_func_coeffs(const Function<T, NDIM> &f,
                               const Key<NDIM> &key) {
     const Tensor<double> &cell_width = FunctionDefaults<NDIM>::get_cell_width();
@@ -221,25 +222,29 @@ template <typename T, std::size_t NDIM> struct FunctionIOData {
         nl.push_back(key_i);
         coords.push_back(std::vector<coordT>());
 
-        for (long i = 0; i < npts_per_box; i++) {
-          for (long j = 0; j < NDIM; j++) {
-            c[j]=cell(j,0)+h*cell_width[j] *(l[j]+qx)
+        if (NDIM == 1) {
+          for (long i = 0; i < k; ++i) {
+            c[0] = cell(0, 0) + h * cell_width[0] * (l[0] + qx(i)); // x
+            coords.back().push_back(c);
           }
-        }
-
-        if (NDIM == 3) {
+        } else if (NDIM == 2) {
           for (long i = 0; i < k; ++i) {
             c[0] = cell(0, 0) + h * cell_width[0] * (l[0] + qx(i)); // x
             for (long j = 0; j < k; ++j) {
               c[1] = cell(1, 0) + h * cell_width[1] * (l[1] + qx(j)); // y
-              for (long m = 0; m < k; ++m) {
-                c[2] = cell(2, 0) + h * cell_width[2] * (l[2] + qx(m)); // z
-                coords.back().push_back(c);
-              }
+              coords.back().push_back(c);
+            }
+          }
+        } else if (NDIM == 3) {
+          for (long j = 0; j < k; ++j) {
+            c[1] = cell(1, 0) + h * cell_width[1] * (l[1] + qx(j)); // y
+            for (long m = 0; m < k; ++m) {
+              c[2] = cell(2, 0) + h * cell_width[2] * (l[2] + qx(m)); // z
+              coords.back().push_back(c);
             }
           }
         } else {
-          MADNESS_EXCEPTION("only NDIM=3 in print_grid", 0);
+          MADNESS_EXCEPTION("only NDIM <= 3 in print_grid", 0);
         }
 
         std::vector<double> values_i(npts_per_box);
@@ -307,12 +312,12 @@ template <typename T, std::size_t NDIM> struct FunctionIOData {
 template <typename T, std::size_t NDIM>
 void to_json(json &j, const FunctionIOData<T, NDIM> &p) {
   j = json{{"npts_per_box", p.npts_per_box},
+           {"num_leaf_nodes", p.num_leaf_nodes},
            {"k", p.k},
            {"cell", p.cell},
-           {"num_leaf_nodes", p.num_leaf_nodes},
-           {"coords", p.coords},
            {"nl", p.nl},
            {"ndim", p.ndim},
+           {"coords", p.coords},
            {"values", p.values}};
 }