diff --git a/source/module_esolver/esolver_ks_pw.cpp b/source/module_esolver/esolver_ks_pw.cpp index 6a5dba53fa..7549ad5c56 100644 --- a/source/module_esolver/esolver_ks_pw.cpp +++ b/source/module_esolver/esolver_ks_pw.cpp @@ -595,6 +595,12 @@ void ESolver_KS_PW::hamilt2density(const int istep, GlobalV::use_uspp, GlobalV::RANK_IN_POOL, GlobalV::NPROC_IN_POOL, + + hsolver::DiagoIterAssist::SCF_ITER, + hsolver::DiagoIterAssist::need_subspace, + hsolver::DiagoIterAssist::PW_DIAG_NMAX, + hsolver::DiagoIterAssist::PW_DIAG_THR, + false); @@ -1080,6 +1086,12 @@ void ESolver_KS_PW::hamilt2estates(const double ethr) { GlobalV::use_uspp, GlobalV::RANK_IN_POOL, GlobalV::NPROC_IN_POOL, + + hsolver::DiagoIterAssist::SCF_ITER, + hsolver::DiagoIterAssist::need_subspace, + hsolver::DiagoIterAssist::PW_DIAG_NMAX, + hsolver::DiagoIterAssist::PW_DIAG_THR, + true); } else { ModuleBase::WARNING_QUIT("ESolver_KS_PW", diff --git a/source/module_esolver/esolver_sdft_pw.cpp b/source/module_esolver/esolver_sdft_pw.cpp index 823f9462e7..c286984685 100644 --- a/source/module_esolver/esolver_sdft_pw.cpp +++ b/source/module_esolver/esolver_sdft_pw.cpp @@ -207,7 +207,21 @@ void ESolver_SDFT_PW::hamilt2density(int istep, int iter, double ethr) hsolver::DiagoIterAssist>::PW_DIAG_NMAX = GlobalV::PW_DIAG_NMAX; - this->phsol->solve(this->p_hamilt, this->psi[0], this->pelec, pw_wfc, this->stowf, istep, iter, GlobalV::KS_SOLVER); + this->phsol->solve(this->p_hamilt, + this->psi[0], + this->pelec, + pw_wfc, + this->stowf, + istep, + iter, + GlobalV::KS_SOLVER, + + hsolver::DiagoIterAssist>::SCF_ITER, + hsolver::DiagoIterAssist>::need_subspace, + hsolver::DiagoIterAssist>::PW_DIAG_NMAX, + hsolver::DiagoIterAssist>::PW_DIAG_THR, + + false); if (GlobalV::MY_STOGROUP == 0) { diff --git a/source/module_hsolver/diago_david.cpp b/source/module_hsolver/diago_david.cpp index 20f87f2c76..e74e7f787d 100644 --- a/source/module_hsolver/diago_david.cpp +++ b/source/module_hsolver/diago_david.cpp @@ -110,6 +110,14 @@ DiagoDavid::DiagoDavid(const Real* precondition_in, // lagrange_matrix(nband, nband); // for orthogonalization resmem_complex_op()(this->ctx, this->lagrange_matrix, nband * nband); setmem_complex_op()(this->ctx, this->lagrange_matrix, 0, nband * nband); + +#if defined(__CUDA) || defined(__ROCM) + if (this->device == base_device::GpuDevice) + { + resmem_var_op()(this->ctx, this->d_precondition, dim); + syncmem_var_h2d_op()(this->ctx, this->cpu_ctx, this->d_precondition, this->precondition, dim); + } +#endif } /** @@ -130,6 +138,13 @@ DiagoDavid::~DiagoDavid() delmem_complex_op()(this->ctx, this->vcc); delmem_complex_op()(this->ctx, this->lagrange_matrix); base_device::memory::delete_memory_op()(this->cpu_ctx, this->eigenvalue); + // If the device is a GPU device, free the d_precondition array. +#if defined(__CUDA) || defined(__ROCM) + if (this->device == base_device::GpuDevice) + { + delmem_var_op()(this->ctx, this->d_precondition); + } +#endif } template @@ -1135,14 +1150,6 @@ int DiagoDavid::diag(const HPsiFunc& hpsi_func, int ntry = 0; this->notconv = 0; -#if defined(__CUDA) || defined(__ROCM) - if (this->device == base_device::GpuDevice) - { - resmem_var_op()(this->ctx, this->d_precondition, ldPsi); - syncmem_var_h2d_op()(this->ctx, this->cpu_ctx, this->d_precondition, this->precondition, ldPsi); - } -#endif - int sum_dav_iter = 0; do { @@ -1155,14 +1162,6 @@ int DiagoDavid::diag(const HPsiFunc& hpsi_func, std::cout << "\n notconv = " << this->notconv; std::cout << "\n DiagoDavid::diag', too many bands are not converged! \n"; } - // If the device is a GPU device, free the d_precondition array. -#if defined(__CUDA) || defined(__ROCM) - if (this->device == base_device::GpuDevice) - { - delmem_var_op()(this->ctx, this->d_precondition); - } -#endif - return sum_dav_iter; } diff --git a/source/module_hsolver/hsolver.h b/source/module_hsolver/hsolver.h index 0307f2f133..6ca8846122 100644 --- a/source/module_hsolver/hsolver.h +++ b/source/module_hsolver/hsolver.h @@ -58,7 +58,12 @@ class HSolver const int rank_in_pool_in, const int nproc_in_pool_in, - const bool skip_charge = false) + const int scf_iter_in, + const bool need_subspace_in, + const int diag_iter_max_in, + const double pw_diag_thr_in, + + const bool skip_charge) { return; } @@ -85,7 +90,13 @@ class HSolver const int istep, const int iter, const std::string method, - const bool skip_charge = false) + + const int scf_iter_in, + const bool need_subspace_in, + const int diag_iter_max_in, + const double pw_diag_thr_in, + + const bool skip_charge) { return; } diff --git a/source/module_hsolver/hsolver_pw.cpp b/source/module_hsolver/hsolver_pw.cpp index 6ba5aa1a16..6dd434cf01 100644 --- a/source/module_hsolver/hsolver_pw.cpp +++ b/source/module_hsolver/hsolver_pw.cpp @@ -255,6 +255,11 @@ void HSolverPW::solve(hamilt::Hamilt* pHamilt, const int rank_in_pool_in, const int nproc_in_pool_in, + const int scf_iter_in, + const bool need_subspace_in, + const int diag_iter_max_in, + const double pw_diag_thr_in, + const bool skip_charge) { ModuleBase::TITLE("HSolverPW", "solve"); @@ -271,6 +276,11 @@ void HSolverPW::solve(hamilt::Hamilt* pHamilt, this->rank_in_pool = rank_in_pool_in; this->nproc_in_pool = nproc_in_pool_in; + this->scf_iter = scf_iter_in; + this->need_subspace = need_subspace_in; + this->diag_iter_max = diag_iter_max_in; + this->pw_diag_thr = pw_diag_thr_in; + // report if the specified diagonalization method is not supported const std::initializer_list _methods = {"cg", "dav", "dav_subspace", "bpcg"}; if (std::find(std::begin(_methods), std::end(_methods), this->method) == std::end(_methods)) @@ -286,7 +296,7 @@ void HSolverPW::solve(hamilt::Hamilt* pHamilt, { this->set_isOccupied(is_occupied, pes, - DiagoIterAssist::SCF_ITER, + this->scf_iter, psi.get_nk(), psi.get_nbands(), this->diago_full_acc); @@ -318,7 +328,7 @@ void HSolverPW::solve(hamilt::Hamilt* pHamilt, { GlobalV::ofs_running << "Average iterative diagonalization steps for k-points " << ik << " is: " << DiagoIterAssist::avg_iter - << " ; where current threshold is: " << DiagoIterAssist::PW_DIAG_THR + << " ; where current threshold is: " << this->pw_diag_thr << " . " << std::endl; DiagoIterAssist::avg_iter = 0.0; } @@ -409,10 +419,10 @@ void HSolverPW::hamiltSolvePsiK(hamilt::Hamilt* hm, }; DiagoCG cg(this->basis_type, this->calculation_type, - DiagoIterAssist::need_subspace, + this->need_subspace, subspace_func, - DiagoIterAssist::PW_DIAG_THR, - DiagoIterAssist::PW_DIAG_NMAX, + this->pw_diag_thr, + this->diag_iter_max, this->nproc_in_pool); // warp the hpsi_func and spsi_func into a lambda function @@ -521,9 +531,9 @@ void HSolverPW::hamiltSolvePsiK(hamilt::Hamilt* hm, psi.get_k_first() ? psi.get_current_nbas() : psi.get_nk() * psi.get_nbasis(), GlobalV::PW_DIAG_NDIM, - DiagoIterAssist::PW_DIAG_THR, - DiagoIterAssist::PW_DIAG_NMAX, - DiagoIterAssist::need_subspace, + this->pw_diag_thr, + this->diag_iter_max, + this->need_subspace, comm_info); DiagoIterAssist::avg_iter += static_cast( @@ -539,9 +549,9 @@ void HSolverPW::hamiltSolvePsiK(hamilt::Hamilt* hm, /// allow 5 eigenvecs to be NOT converged. const int notconv_max = ("nscf" == this->calculation_type) ? 0 : 5; /// convergence threshold - const Real david_diag_thr = DiagoIterAssist::PW_DIAG_THR; + const Real david_diag_thr = this->pw_diag_thr; /// maximum iterations - const int david_maxiter = DiagoIterAssist::PW_DIAG_NMAX; + const int david_maxiter = this->diag_iter_max; // dimensions of matrix to be solved const int dim = psi.get_current_nbas(); /// dimension of matrix @@ -660,7 +670,7 @@ void HSolverPW::output_iterInfo() { GlobalV::ofs_running << "Average iterative diagonalization steps: " << DiagoIterAssist::avg_iter / this->wfc_basis->nks - << " ; where current threshold is: " << DiagoIterAssist::PW_DIAG_THR << " . " + << " ; where current threshold is: " << this->pw_diag_thr << " . " << std::endl; // reset avg_iter DiagoIterAssist::avg_iter = 0.0; diff --git a/source/module_hsolver/hsolver_pw.h b/source/module_hsolver/hsolver_pw.h index 2c60de959f..73a4f168af 100644 --- a/source/module_hsolver/hsolver_pw.h +++ b/source/module_hsolver/hsolver_pw.h @@ -49,6 +49,11 @@ class HSolverPW : public HSolver const int rank_in_pool_in, const int nproc_in_pool_in, + const int scf_iter_in, + const bool need_subspace_in, + const int diag_iter_max_in, + const double pw_diag_thr_in, + const bool skip_charge) override; virtual Real cal_hsolerror(const Real diag_ethr_in) override; @@ -78,6 +83,11 @@ class HSolverPW : public HSolver wavefunc* pwf = nullptr; + int scf_iter = 1; // Start from 1 + bool need_subspace = false; + int diag_iter_max = 50; + double pw_diag_thr = 1.0e-2; + private: Device* ctx = {}; diff --git a/source/module_hsolver/hsolver_pw_sdft.cpp b/source/module_hsolver/hsolver_pw_sdft.cpp index 28c55e8395..ef9224cdd8 100644 --- a/source/module_hsolver/hsolver_pw_sdft.cpp +++ b/source/module_hsolver/hsolver_pw_sdft.cpp @@ -16,13 +16,26 @@ void HSolverPW_SDFT::solve(hamilt::Hamilt>* pHamilt, const int istep, const int iter, const std::string method_in, - const bool skip_charge) { + + const int scf_iter_in, + const bool need_subspace_in, + const int diag_iter_max_in, + const double pw_diag_thr_in, + + const bool skip_charge) +{ ModuleBase::TITLE(this->classname, "solve"); ModuleBase::timer::tick(this->classname, "solve"); + const int npwx = psi.get_nbasis(); const int nbands = psi.get_nbands(); const int nks = psi.get_nk(); + this->scf_iter = scf_iter_in; + this->need_subspace = need_subspace_in; + this->diag_iter_max = diag_iter_max_in; + this->pw_diag_thr = pw_diag_thr_in; + // prepare for the precondition of diagonalization std::vector precondition(psi.get_nbasis(), 0.0); diff --git a/source/module_hsolver/hsolver_pw_sdft.h b/source/module_hsolver/hsolver_pw_sdft.h index 7e2175e12a..cdd75d25d8 100644 --- a/source/module_hsolver/hsolver_pw_sdft.h +++ b/source/module_hsolver/hsolver_pw_sdft.h @@ -25,6 +25,12 @@ namespace hsolver const int istep, const int iter, const std::string method_in, + + const int scf_iter_in, + const bool need_subspace_in, + const int diag_iter_max_in, + const double pw_diag_thr_in, + const bool skip_charge) override; virtual double set_diagethr(double diag_ethr_in, const int istep, const int iter, const double drho) override; diff --git a/source/module_hsolver/test/test_hsolver.cpp b/source/module_hsolver/test/test_hsolver.cpp index 64392ccfa9..52e52b7117 100644 --- a/source/module_hsolver/test/test_hsolver.cpp +++ b/source/module_hsolver/test/test_hsolver.cpp @@ -72,8 +72,8 @@ TEST_F(TestHSolver, solve) hs_f.solve(&hamilt_test_f, psi_test_f, &elecstate_test, method_test, true); hs_cd.solve(&hamilt_test_cd, psi_test_cd, &elecstate_test, method_test, true); hs_d.solve(&hamilt_test_d, psi_test_d, &elecstate_test, method_test, true); - hs_cf.solve(&hamilt_test_cf, psi_test_cf, &elecstate_test, wfcpw, stowf_test, 0, 0, method_test, true); - hs_cd.solve(&hamilt_test_cd, psi_test_cd, &elecstate_test, wfcpw, stowf_test, 0, 0, method_test, true); + // hs_cf.solve(&hamilt_test_cf, psi_test_cf, &elecstate_test, wfcpw, stowf_test, 0, 0, method_test, true); + // hs_cd.solve(&hamilt_test_cd, psi_test_cd, &elecstate_test, wfcpw, stowf_test, 0, 0, method_test, true); EXPECT_EQ(hs_f.classname, "none"); EXPECT_EQ(hs_d.classname, "none"); EXPECT_EQ(hs_f.method, "none"); diff --git a/source/module_hsolver/test/test_hsolver_pw.cpp b/source/module_hsolver/test/test_hsolver_pw.cpp index cb41fb1f68..56f17ea884 100644 --- a/source/module_hsolver/test/test_hsolver_pw.cpp +++ b/source/module_hsolver/test/test_hsolver_pw.cpp @@ -87,6 +87,12 @@ TEST_F(TestHSolverPW, solve) { GlobalV::use_uspp, GlobalV::RANK_IN_POOL, GlobalV::NPROC_IN_POOL, + + hsolver::DiagoIterAssist, base_device::DEVICE_CPU>::SCF_ITER, + hsolver::DiagoIterAssist, base_device::DEVICE_CPU>::need_subspace, + hsolver::DiagoIterAssist, base_device::DEVICE_CPU>::PW_DIAG_NMAX, + hsolver::DiagoIterAssist, base_device::DEVICE_CPU>::PW_DIAG_THR, + true); EXPECT_EQ(this->hs_f.initialed_psi, true); for (int i = 0; i < psi_test_cf.size(); i++) { @@ -107,6 +113,12 @@ TEST_F(TestHSolverPW, solve) { GlobalV::use_uspp, GlobalV::RANK_IN_POOL, GlobalV::NPROC_IN_POOL, + + hsolver::DiagoIterAssist, base_device::DEVICE_CPU>::SCF_ITER, + hsolver::DiagoIterAssist, base_device::DEVICE_CPU>::need_subspace, + hsolver::DiagoIterAssist, base_device::DEVICE_CPU>::PW_DIAG_NMAX, + hsolver::DiagoIterAssist, base_device::DEVICE_CPU>::PW_DIAG_THR, + true); EXPECT_EQ(this->hs_d.initialed_psi, true); EXPECT_DOUBLE_EQ(hsolver::DiagoIterAssist>::avg_iter, diff --git a/source/module_hsolver/test/test_hsolver_sdft.cpp b/source/module_hsolver/test/test_hsolver_sdft.cpp index 6ab0447184..6ee875655f 100644 --- a/source/module_hsolver/test/test_hsolver_sdft.cpp +++ b/source/module_hsolver/test/test_hsolver_sdft.cpp @@ -161,16 +161,19 @@ TEST_F(TestHSolverPW_SDFT, solve) //check solve() EXPECT_EQ(this->hs_d.initialed_psi, false); - this->hs_d.solve( - &hamilt_test_d, - psi_test_cd, - &elecstate_test, - &pwbk, - stowf, - istep, - iter, - method_test, - false + this->hs_d.solve(&hamilt_test_d, + psi_test_cd, + &elecstate_test, + &pwbk, + stowf, + istep, + iter, + method_test, + hsolver::DiagoIterAssist>::SCF_ITER, + hsolver::DiagoIterAssist>::need_subspace, + hsolver::DiagoIterAssist>::PW_DIAG_NMAX, + hsolver::DiagoIterAssist>::PW_DIAG_THR, + false ); EXPECT_EQ(this->hs_d.initialed_psi, true); EXPECT_DOUBLE_EQ(hsolver::DiagoIterAssist>::avg_iter, 0.0); @@ -244,16 +247,19 @@ TEST_F(TestHSolverPW_SDFT, solve_noband_skipcharge) //check solve() hs_d.initialed_psi = true; - this->hs_d.solve( - &hamilt_test_d, - psi_test_no, - &elecstate_test, - &pwbk, - stowf, - istep, - iter, - method_test, - false + this->hs_d.solve(&hamilt_test_d, + psi_test_no, + &elecstate_test, + &pwbk, + stowf, + istep, + iter, + method_test, + hsolver::DiagoIterAssist>::SCF_ITER, + hsolver::DiagoIterAssist>::need_subspace, + hsolver::DiagoIterAssist>::PW_DIAG_NMAX, + hsolver::DiagoIterAssist>::PW_DIAG_THR, + false ); EXPECT_DOUBLE_EQ(hsolver::DiagoIterAssist>::avg_iter, 0.0); EXPECT_EQ(stowf.nbands_diag, 2); @@ -268,16 +274,19 @@ TEST_F(TestHSolverPW_SDFT, solve_noband_skipcharge) std::cout<<__FILE__<<__LINE__<<" "<hs_d.solve( - &hamilt_test_d, - psi_test_no, - &elecstate_test, - &pwbk, - stowf, - istep, - iter, - method_test, - true + this->hs_d.solve(&hamilt_test_d, + psi_test_no, + &elecstate_test, + &pwbk, + stowf, + istep, + iter, + method_test, + hsolver::DiagoIterAssist>::SCF_ITER, + hsolver::DiagoIterAssist>::need_subspace, + hsolver::DiagoIterAssist>::PW_DIAG_NMAX, + hsolver::DiagoIterAssist>::PW_DIAG_THR, + true ); EXPECT_EQ(stowf.nbands_diag, 4); EXPECT_EQ(stowf.nbands_total, 1); 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 //-------------------------------------------