From 5c0f0d5f691d1234f49718b99a6370d1835d8d34 Mon Sep 17 00:00:00 2001 From: Yu Liu <77716030+YuLiu98@users.noreply.github.com> Date: Sat, 27 Jul 2024 20:54:56 +0800 Subject: [PATCH] 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 | 153 +++++++++++++++++++++++ source/module_io/output_log.h | 20 +++ source/module_io/test/CMakeLists.txt | 2 +- source/module_io/test/outputlog_test.cpp | 73 +++++++++++ source/module_io/write_pot.cpp | 15 +++ 5 files changed, 262 insertions(+), 1 deletion(-) diff --git a/source/module_io/output_log.cpp b/source/module_io/output_log.cpp index b0402abba7..bbdcf88999 100644 --- a/source/module_io/output_log.cpp +++ b/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 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 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 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, diff --git a/source/module_io/output_log.h b/source/module_io/output_log.h index afea63dc3a..76bb6d293d 100644 --- a/source/module_io/output_log.h +++ b/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 diff --git a/source/module_io/test/CMakeLists.txt b/source/module_io/test/CMakeLists.txt index ec877ad419..c065e06f14 100644 --- a/source/module_io/test/CMakeLists.txt +++ b/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( diff --git a/source/module_io/test/outputlog_test.cpp b/source/module_io/test/outputlog_test.cpp index 903e6be9b7..adda80ee67 100644 --- a/source/module_io/test/outputlog_test.cpp +++ b/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[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; } \ No newline at end of file diff --git a/source/module_io/write_pot.cpp b/source/module_io/write_pot.cpp index 133680ccbc..d309117365 100644 --- a/source/module_io/write_pot.cpp +++ b/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 //-------------------------------------------