Feature: output vacuum level when output electrostatic potential (#4799)

* Feature: output vacuum level

* Tests: add unit tests for vacuum level

* Test: update CMakeLists.txt for unittests

* Refactor: remove rho_basis dependency in vacuum level

source/module_io/output_log.cpp

@@ -29,6 +29,159 @@ void output_efermi(bool& convergence, double& efermi, std::ofstream& ofs_running
     }
 }
 
+void output_vacuum_level(const UnitCell* ucell,
+                         const double* const* rho,
+                         const double* v_elecstat,
+                         const int& nx,
+                         const int& ny,
+                         const int& nz,
+                         const int& nxyz,
+                         const int& nrxx,
+                         const int& nplane,
+                         const int& startz_current,
+                         std::ofstream& ofs_running)
+{
+    // determine the vacuum direction
+    double vacuum[3] = {0.0};
+    for (int dir = 0; dir < 3; dir++)
+    {
+        std::vector<double> pos;
+        for (int it = 0; it < ucell->ntype; ++it)
+        {
+            for (int ia = 0; ia < ucell->atoms[it].na; ++ia)
+            {
+                pos.push_back(ucell->atoms[it].taud[ia][dir]);
+            }
+        }
+
+        std::sort(pos.begin(), pos.end());
+        for (int i = 1; i < pos.size(); i++)
+        {
+            vacuum[dir] = std::max(vacuum[dir], pos[i] - pos[i - 1]);
+        }
+
+        // consider the periodic boundary condition
+        vacuum[dir] = std::max(vacuum[dir], pos[0] + 1 - pos[pos.size() - 1]);
+    }
+
+    // we assume that the cell is a cuboid
+    // get the direction with the largest vacuum
+    int direction = 2;
+    vacuum[0] *= ucell->latvec.e11;
+    vacuum[1] *= ucell->latvec.e22;
+    vacuum[2] *= ucell->latvec.e33;
+    if (vacuum[0] > vacuum[2])
+    {
+        direction = 0;
+    }
+    if (vacuum[1] > vacuum[direction])
+    {
+        direction = 1;
+    }
+
+    int length = nz;
+    if (direction == 0)
+    {
+        length = nx;
+    }
+    else if (direction == 1)
+    {
+        length = ny;
+    }
+
+    // get the average along the direction in real space
+    auto average = [](const int& ny,
+                      const int& nxyz,
+                      const int& nrxx,
+                      const int& nplane,
+                      const int& startz_current,
+                      const int& direction,
+                      const int& length,
+                      const double* v,
+                      double* ave,
+                      bool abs) {
+        for (int ir = 0; ir < nrxx; ++ir)
+        {
+            int index = 0;
+            if (direction == 0)
+            {
+                index = ir / (ny * nplane);
+            }
+            else if (direction == 1)
+            {
+                int i = ir / (ny * nplane);
+                index = ir / nplane - i * ny;
+            }
+            else if (direction == 2)
+            {
+                index = ir % nplane + startz_current;
+            }
+
+            double value = abs ? std::fabs(v[ir]) : v[ir];
+
+            ave[index] += value;
+        }
+
+#ifdef __MPI
+        MPI_Allreduce(MPI_IN_PLACE, ave, length, MPI_DOUBLE, MPI_SUM, POOL_WORLD);
+#endif
+
+        int surface = nxyz / length;
+        for (int i = 0; i < length; ++i)
+        {
+            ave[i] /= surface;
+        }
+    };
+
+    // average charge density along direction
+    std::vector<double> totchg(nrxx, 0.0);
+    for (int ir = 0; ir < nrxx; ++ir)
+    {
+        totchg[ir] = rho[0][ir];
+    }
+    if (GlobalV::NSPIN == 2)
+    {
+        for (int ir = 0; ir < nrxx; ++ir)
+        {
+            totchg[ir] += rho[1][ir];
+        }
+    }
+
+    std::vector<double> ave(length, 0.0);
+    average(ny, nxyz, nrxx, nplane, startz_current, direction, length, totchg.data(), ave.data(), true);
+
+    // set vacuum to be the point in space where the electronic charge density is the minimum
+    // get the index corresponding to the minimum charge density
+    int min_index = 0;
+    double min_value = 1e9;
+    double windows[7] = {0.1, 0.2, 0.3, 0.4, 0.3, 0.2, 0.1};
+    for (int i = 0; i < length; i++)
+    {
+        double sum = 0;
+        int temp = i - 3 + length;
+        // use a sliding average to smoothen in charge density
+        for (int win = 0; win < 7; win++)
+        {
+            int index = (temp + win) % length;
+            sum += ave[index] * windows[win];
+        }
+
+        if (sum < min_value)
+        {
+            min_value = sum;
+            min_index = i;
+        }
+    }
+
+    // average electrostatic potential along direction
+    ave.assign(ave.size(), 0.0);
+    average(ny, nxyz, nrxx, nplane, startz_current, direction, length, v_elecstat, ave.data(), false);
+
+    // get the vacuum level
+    double vacuum_level = ave[min_index] * ModuleBase::Ry_to_eV;
+    ofs_running << "The vacuum level is " << vacuum_level << " eV" << std::endl;
+}
+
 void print_force(std::ofstream& ofs_running,
                  const UnitCell& cell,
                  const std::string& name,

source/module_io/output_log.h

@@ -22,6 +22,26 @@ void output_convergence_after_scf(bool& convergence, double& energy, std::ofstre
 /// @param ofs_running the output stream
 void output_efermi(bool& convergence, double& efermi, std::ofstream& ofs_running = GlobalV::ofs_running);
 
+/// @brief calculate and output the vacuum level
+/// We first determine the vacuum direction, then get the vacuum position based on the minimum of charge density,
+/// finally output the vacuum level, i.e., the electrostatic potential at the vacuum position.
+///
+/// @param ucell the unitcell
+/// @param rho charge density
+/// @param v_elecstat electrostatic potential
+/// @param ofs_running the output stream
+void output_vacuum_level(const UnitCell* ucell,
+                         const double* const* rho,
+                         const double* v_elecstat,
+                         const int& nx,
+                         const int& ny,
+                         const int& nz,
+                         const int& nxyz,
+                         const int& nrxx,
+                         const int& nplane,
+                         const int& startz_current,
+                         std::ofstream& ofs_running = GlobalV::ofs_running);
+
 /// @brief output atomic forces
 /// @param ofs the output stream
 /// @param cell the unitcell

source/module_io/test/CMakeLists.txt

@@ -103,7 +103,7 @@ AddTest(
 AddTest(
   TARGET io_output_log_test
   LIBS base ${math_libs} device
-  SOURCES ../output_log.cpp outputlog_test.cpp
+  SOURCES ../output_log.cpp outputlog_test.cpp ../../module_basis/module_pw/test/test_tool.cpp
 )
 
 AddTest(

source/module_io/test/outputlog_test.cpp

@@ -9,6 +9,11 @@
 #include "module_base/constants.h"
 #include "module_base/global_variable.h"
 #include "module_io/output_log.h"
+
+#ifdef __MPI
+#include "module_basis/module_pw/test/test_tool.h"
+#include "mpi.h"
+#endif
 /**
  * - Tested Functions:
  *  - output_convergence_after_scf()
@@ -123,9 +128,19 @@ UnitCell::UnitCell()
     ntype = 1;
     nat = 2;
     atoms = new Atom[ntype];
+
+    latvec.e11 = latvec.e22 = latvec.e33 = 10;
+    latvec.e12 = latvec.e13 = latvec.e23 = 0;
+    latvec.e21 = latvec.e31 = latvec.e32 = 0;
+
+    atoms[0].taud = new ModuleBase::Vector3<double>[2];
+    atoms[0].taud[0].set(0.5456, 0, 0.54275);
+    atoms[0].taud[1].set(0.54, 0.8495, 0.34175);
 }
 UnitCell::~UnitCell()
 {
+    if (atoms != nullptr)
+        delete[] atoms;
 }
 InfoNonlocal::InfoNonlocal()
 {
@@ -146,6 +161,8 @@ Atom::Atom()
 }
 Atom::~Atom()
 {
+    if (taud != nullptr)
+        delete[] taud;
 }
 Atom_pseudo::Atom_pseudo()
 {
@@ -160,6 +177,41 @@ pseudo::~pseudo()
 {
 }
 
+TEST(OutputVacuumLevelTest, OutputVacuumLevel)
+{
+    GlobalV::NSPIN = 1;
+    UnitCell ucell;
+    const int nx = 50, ny = 50, nz = 50, nxyz = 125000, nrxx = 125000, nplane = 50, startz_current = 0;
+
+    double** rho = new double*[1];
+    rho[0] = new double[nrxx];
+    double* v_elecstat = new double[nrxx];
+    for (int ir = 0; ir < nrxx; ++ir)
+    {
+        rho[0][ir] = 0.01 * ir;
+        v_elecstat[ir] = 0.02 * ir;
+    }
+
+    std::ofstream ofs_running("test_output_vacuum_level.txt");
+    ModuleIO::output_vacuum_level(&ucell, rho, v_elecstat, nx, ny, nz, nxyz, nrxx, nplane, startz_current, ofs_running);
+    ofs_running.close();
+
+    std::ifstream ifs_running("test_output_vacuum_level.txt");
+    std::stringstream ss;
+    ss << ifs_running.rdbuf();
+    std::string file_content = ss.str();
+    ifs_running.close();
+
+    std::string expected_content = "The vacuum level is 2380.86 eV\n";
+
+    EXPECT_EQ(file_content, expected_content);
+    std::remove("test_output_vacuum_level.txt");
+
+    delete[] rho[0];
+    delete[] rho;
+    delete[] v_elecstat;
+}
+
 TEST(PrintForce, PrintForce)
 {
     UnitCell ucell;
@@ -237,4 +289,25 @@ TEST(PrintStress, PrintStress)
     EXPECT_THAT(output_str, testing::HasSubstr(" TOTAL-PRESSURE: 49035.075992 KBAR"));
     ifs.close();
     std::remove("test.txt");
+}
+
+int main(int argc, char** argv)
+{
+#ifdef __MPI
+    setupmpi(argc, argv, GlobalV::NPROC, GlobalV::MY_RANK);
+    divide_pools(GlobalV::NPROC,
+                 GlobalV::MY_RANK,
+                 GlobalV::NPROC_IN_POOL,
+                 GlobalV::KPAR,
+                 GlobalV::MY_POOL,
+                 GlobalV::RANK_IN_POOL);
+#endif
+
+    testing::InitGoogleTest(&argc, argv);
+    int result = RUN_ALL_TESTS();
+
+#ifdef __MPI
+    finishmpi();
+#endif
+    return result;
 }

source/module_io/write_pot.cpp

@@ -6,6 +6,7 @@
 #include "module_elecstate/potentials/efield.h"
 #include "module_hamilt_pw/hamilt_pwdft/global.h"
 #include "module_io/cube_io.h"
+#include "module_io/output_log.h"
 
 namespace ModuleIO
 {
@@ -178,6 +179,20 @@ void write_elecstat_pot(
         }
     }
 
+    //-------------------------------------------
+    //! Get the vacuum level of the system
+    //-------------------------------------------
+    ModuleIO::output_vacuum_level(ucell,
+                                  chr->rho,
+                                  v_elecstat.data(),
+                                  rho_basis->nx,
+                                  rho_basis->ny,
+                                  rho_basis->nz,
+                                  rho_basis->nxyz,
+                                  rho_basis->nrxx,
+                                  rho_basis->nplane,
+                                  rho_basis->startz_current);
+
     //-------------------------------------------
     //! Write down the electrostatic potential
     //-------------------------------------------