Feature: add ModuleIO::write_libxc_r()

Author: PeizeLin

source/Makefile.Objects

@@ -504,6 +504,7 @@ OBJS_IO=input_conv.o\
     write_dipole.o\
     td_current_io.o\
     write_wfc_r.o\
+    write_libxc_r.o\
     output_log.o\
     output_mat_sparse.o\
     para_json.o\

source/module_esolver/esolver_fp.cpp

@@ -8,6 +8,7 @@
 #include "module_parameter/parameter.h"
 #include "module_io/rhog_io.h"
 #include "module_io/cif_io.h"
+#include "module_io/write_libxc_r.h"
 
 namespace ModuleESolver
 {
@@ -251,7 +252,19 @@ void ESolver_FP::after_scf(const int istep)
                 &(GlobalC::ucell),
                 this->pelec->pot->get_fixed_v());
         }
-    }
+    
+		if(PARAM.inp.out_xc_r[0]>=0)
+		{
+			ModuleIO::write_libxc_r(
+				PARAM.inp.out_xc_r[0],
+				XC_Functional::get_func_id(),
+				this->pw_rhod->nrxx, // number of real-space grid
+				GlobalC::ucell.omega, // volume of cell
+				GlobalC::ucell.tpiba,
+				&this->chr,
+				this);
+		}
+	}
 }
 
 void ESolver_FP::init_after_vc(const Input_para& inp, UnitCell& cell)

source/module_io/CMakeLists.txt

@@ -26,6 +26,7 @@ list(APPEND objects
     write_dipole.cpp
     td_current_io.cpp
     write_wfc_r.cpp
+    write_libxc_r.cpp
     output_log.cpp
     para_json.cpp
     parse_args.cpp

source/module_io/cube_io.h

@@ -49,6 +49,10 @@ extern void write_cube(
     const int out_fermi = 1); // mohan add 2007-10-17
 
 
+// when MPI:
+//      read file as order (ixy,iz) to data[ixy*nz+iz]
+// when serial:
+//      read file as order (ixy,iz) to data[iz*nxy+ixy]
 extern void read_cube_core_match(
     std::ifstream &ifs,
 #ifdef __MPI
@@ -73,18 +77,30 @@ extern void read_cube_core_mismatch(
     const int ny_read,
     const int nz_read);
 
+#ifdef __MPI
+// when MPI:
+//      write data[ixy*nplane+iz*nld] to file as order (ixy,iz)
 extern void write_cube_core(
     std::ofstream &ofs_cube,
-#ifdef __MPI
     const int bz,
     const int nbz,
     const int nplane,
     const int startz_current,
-#endif
     const double*const data,
     const int nxy,
     const int nz,
+    const int nld,
     const int n_data_newline);
+#else
+// when serial:
+//      write data[iz*nxy+ixy] to file as order (ixy,iz)
+extern void write_cube_core(
+    std::ofstream &ofs_cube,
+    const double*const data,
+    const int nxy,
+    const int nz,
+    const int n_data_newline);
+#endif
 
     /**
      * @brief The trilinear interpolation method

source/module_io/read_input_item_output.cpp

@@ -76,6 +76,21 @@ void ReadInput::item_output()
         read_sync_bool(input.out_wfc_r);
         this->add_item(item);
     }
+    {
+        Input_Item item("out_xc_r");
+        item.annotation = "if >=0, output the derivatives of exchange correlation in realspace, second parameter controls the precision";
+        item.read_value = [](const Input_Item& item, Parameter& para) {
+            size_t count = item.get_size();
+            std::vector<int> out_xc_r(count); // create a placeholder vector
+            std::transform(item.str_values.begin(), item.str_values.end(), out_xc_r.begin(), [](std::string s) { return std::stoi(s); });
+            std::cout<<out_xc_r[0]<<"\t"<<out_xc_r[1]<<std::endl;
+            // assign non-negative values to para.input.out_xc_r
+            std::copy(out_xc_r.begin(), out_xc_r.end(), para.input.out_xc_r.begin());
+            std::cout<<para.input.out_xc_r[0]<<"\t"<<para.input.out_xc_r[1]<<std::endl;
+        };
+        sync_intvec(input.out_xc_r, 2, -1);
+        this->add_item(item);
+    }
     {
         Input_Item item("printe");
         item.annotation = "Print out energy for each band for every printe steps";

source/module_io/write_cube.cpp

@@ -127,7 +127,7 @@ void ModuleIO::write_cube(
     }
 
 #ifdef __MPI
-    ModuleIO::write_cube_core(ofs_cube, bz, nbz, nplane, startz_current, data, nx*ny, nz, 6);
+    ModuleIO::write_cube_core(ofs_cube, bz, nbz, nplane, startz_current, data, nx*ny, nz, 1, 6);
 #else
     ModuleIO::write_cube_core(ofs_cube, data, nx*ny, nz, 6);
 #endif
@@ -140,28 +140,25 @@ void ModuleIO::write_cube(
         /// for cube file
         ofs_cube.close();
     }
-
-    return;
 }
 
 
+#ifdef __MPI
+
 void ModuleIO::write_cube_core(
     std::ofstream &ofs_cube,
-#ifdef __MPI
     const int bz,
     const int nbz,
     const int nplane,
     const int startz_current,
-#endif
     const double*const data,
     const int nxy,
     const int nz,
+    const int nld,
     const int n_data_newline)
 {
     ModuleBase::TITLE("ModuleIO", "write_cube_core");
 
-#ifdef __MPI
-
     const int my_rank = GlobalV::MY_RANK;
     const int my_pool = GlobalV::MY_POOL;
     const int rank_in_pool = GlobalV::RANK_IN_POOL;
@@ -176,7 +173,6 @@ void ModuleIO::write_cube_core(
 
         // num_z: how many planes on processor 'ip'
         std::vector<int> num_z(nproc_in_pool, 0);
-
         for (int iz = 0; iz < nbz; iz++)
         {
             const int ip = iz % nproc_in_pool;
@@ -230,7 +226,7 @@ void ModuleIO::write_cube_core(
                     // mohan change to rho_save on 2012-02-10
                     // because this can make our next restart calculation lead
                     // to the same scf_thr as the one saved.
-                    zpiece[ixy] = data[ixy * nplane + iz - startz_current];
+                    zpiece[ixy] = data[ixy * nplane + (iz - startz_current) * nld];
                 }
             }
             // case 2: > first part rho: send the rho to
@@ -239,7 +235,7 @@ void ModuleIO::write_cube_core(
             {
                 for (int ixy = 0; ixy < nxy; ixy++)
                 {
-                    zpiece[ixy] = data[ixy * nplane + iz - startz_current];
+                    zpiece[ixy] = data[ixy * nplane + (iz - startz_current) * nld];
                 }
                 MPI_Send(zpiece.data(), nxy, MPI_DOUBLE, 0, tag, POOL_WORLD);
             }
@@ -281,7 +277,18 @@ void ModuleIO::write_cube_core(
         /// for cube file
     }
     MPI_Barrier(MPI_COMM_WORLD);
-#else
+}
+
+#else   // #ifdef __MPI
+
+void ModuleIO::write_cube_core(
+    std::ofstream &ofs_cube,
+    const double*const data,
+    const int nxy,
+    const int nz,
+    const int n_data_newline)
+{
+    ModuleBase::TITLE("ModuleIO", "write_cube_core");
     for (int ixy = 0; ixy < nxy; ixy++)
     {
         for (int iz = 0; iz < nz; iz++)
@@ -295,5 +302,6 @@ void ModuleIO::write_cube_core(
         }
         ofs_cube << "\n";
     }
-#endif
-}
+}
+
+#endif  // #ifdef __MPI