From e3d6225064c492e9bf929d5f22c15002b8b9c40a Mon Sep 17 00:00:00 2001 From: AsTonyshment Date: Fri, 3 Jul 2026 16:25:48 +0800 Subject: [PATCH 1/2] Test: Add snap_psibeta checks --- source/Makefile.Objects | 1 + .../source_io/module_hs/cal_r_overlap_R.cpp | 36 ++- .../module_rt/snap_projector_half_tddft.cpp | 9 +- .../module_rt/snap_projector_half_tddft.h | 41 ++- .../module_rt/snap_psibeta_half_tddft.cpp | 4 +- .../module_rt/snap_psibeta_half_tddft.h | 1 + .../source_lcao/module_rt/test/CMakeLists.txt | 17 ++ .../test/snap_psibeta_half_tddft_test.cpp | 248 ++++++++++-------- 8 files changed, 219 insertions(+), 138 deletions(-) diff --git a/source/Makefile.Objects b/source/Makefile.Objects index cd1d3b1f4cb..2240189b206 100644 --- a/source/Makefile.Objects +++ b/source/Makefile.Objects @@ -635,6 +635,7 @@ OBJS_LCAO=evolve_elec.o\ td_folding.o\ td_info.o\ velocity_op.o\ + snap_projector_half_tddft.o\ snap_psibeta_half_tddft.o\ solve_propagation.o\ boundary_fix.o\ diff --git a/source/source_io/module_hs/cal_r_overlap_R.cpp b/source/source_io/module_hs/cal_r_overlap_R.cpp index 57d6eb22971..a1e821e7268 100644 --- a/source/source_io/module_hs/cal_r_overlap_R.cpp +++ b/source/source_io/module_hs/cal_r_overlap_R.cpp @@ -187,11 +187,19 @@ void cal_r_overlap_R::construct_orbs_and_orb_r(const UnitCell& ucell, } } - iw2it.resize(PARAM.globalv.nlocal); - iw2ia.resize(PARAM.globalv.nlocal); - iw2iL.resize(PARAM.globalv.nlocal); - iw2iN.resize(PARAM.globalv.nlocal); - iw2im.resize(PARAM.globalv.nlocal); + int map_size = PARAM.globalv.nlocal; + int required_orbitals = 0; + for (int it = 0; it < ucell.ntype; ++it) + { + required_orbitals += ucell.atoms[it].nw * ucell.atoms[it].na; + } + map_size = std::max(map_size, required_orbitals); + + iw2it.resize(map_size); + iw2ia.resize(map_size); + iw2iL.resize(map_size); + iw2iN.resize(map_size); + iw2im.resize(map_size); int iw = 0; for (int it = 0; it < ucell.ntype; it++) @@ -419,11 +427,19 @@ void cal_r_overlap_R::construct_orbs_and_nonlocal_and_orb_r(const UnitCell& ucel } } - iw2it.resize(PARAM.globalv.nlocal); - iw2ia.resize(PARAM.globalv.nlocal); - iw2iL.resize(PARAM.globalv.nlocal); - iw2iN.resize(PARAM.globalv.nlocal); - iw2im.resize(PARAM.globalv.nlocal); + int map_size = PARAM.globalv.nlocal; + int required_orbitals = 0; + for (int it = 0; it < ucell.ntype; ++it) + { + required_orbitals += ucell.atoms[it].nw * ucell.atoms[it].na; + } + map_size = std::max(map_size, required_orbitals); + + iw2it.resize(map_size); + iw2ia.resize(map_size); + iw2iL.resize(map_size); + iw2iN.resize(map_size); + iw2im.resize(map_size); int iw = 0; for (int it = 0; it < ucell.ntype; it++) diff --git a/source/source_lcao/module_rt/snap_projector_half_tddft.cpp b/source/source_lcao/module_rt/snap_projector_half_tddft.cpp index a2c48584e66..6dcf4f67c83 100644 --- a/source/source_lcao/module_rt/snap_projector_half_tddft.cpp +++ b/source/source_lcao/module_rt/snap_projector_half_tddft.cpp @@ -104,8 +104,7 @@ AngularGridView angular_grid(const int ngrid) { if (!is_supported_lebedev_grid(ngrid)) { - ModuleBase::WARNING_QUIT("snap_projector_half_tddft", - "Unsupported Lebedev-Laikov grid size: " + std::to_string(ngrid)); + ModuleBase::WARNING_QUIT("snap_projector_half_tddft", "Unsupported Lebedev-Laikov grid size: " + std::to_string(ngrid)); } if (ngrid == default_lebedev_grid_points) @@ -142,9 +141,7 @@ AngularGridView angular_grid(const int ngrid) double radial_factor(const ProjectorChannel& channel, const double r, const double w_radial) { - const double projector_val - = ModuleBase::PolyInt::Polynomial_Interpolation(channel.radial_values, channel.mesh, channel.dk, r); - + const double projector_val = ModuleBase::PolyInt::Polynomial_Interpolation(channel.radial_times_r, channel.mesh, channel.dk, r); return projector_val * r * w_radial; } } // namespace @@ -264,7 +261,7 @@ void snap_projector_half_tddft(const LCAO_Orbitals& orb, } assert(channel.mesh > 0); - assert(channel.radial_values != nullptr); + assert(channel.radial_times_r != nullptr); assert(channel.radial_grid != nullptr); const double r_min = channel.radial_grid[0]; diff --git a/source/source_lcao/module_rt/snap_projector_half_tddft.h b/source/source_lcao/module_rt/snap_projector_half_tddft.h index 15baf98854e..c67b8b458fc 100644 --- a/source/source_lcao/module_rt/snap_projector_half_tddft.h +++ b/source/source_lcao/module_rt/snap_projector_half_tddft.h @@ -10,8 +10,23 @@ namespace module_rt { +/** + * @brief Numerical quadrature settings for projector snapshots. + * + * The default values reproduce the production RT-TDDFT path. + */ +struct SnapIntegrationOptions +{ + int radial_grid_num = 140; + int lebedev_grid_points = 110; +}; + /** * @brief Radial projector channel integrated against one LCAO orbital. + * + * radial_times_r stores r * p_l(r), where p_l(r) is the radial projector. + * The radial part of the volume integral is therefore evaluated as + * (r * p_l(r)) * r dr = p_l(r) * r^2 dr. */ struct ProjectorChannel { @@ -19,23 +34,20 @@ struct ProjectorChannel int mesh = 0; double dk = 0.0; double rcut = 0.0; - const double* radial_values = nullptr; + const double* radial_times_r = nullptr; const double* radial_grid = nullptr; }; /** - * @brief Numerical quadrature settings for projector snapshots. + * @brief Compute projector overlaps with default quadrature settings. * - * The default values reproduce the production RT-TDDFT path. - */ -struct SnapIntegrationOptions -{ - int radial_grid_num = 140; - int lebedev_grid_points = 110; -}; - -/** - * @brief Compute with default quadrature settings. + * The shared integral is + * I_m(A) = . + * + * The phase A is given in Cartesian coordinates. The returned nlm[0] stores + * I_m(A) for all projector magnetic components. If calc_r is true, nlm[1..3] + * store + * R_a,m(A) = . */ void snap_projector_half_tddft(const LCAO_Orbitals& orb, const std::vector& projector_channels, @@ -51,9 +63,10 @@ void snap_projector_half_tddft(const LCAO_Orbitals& orb, const char* timer_name); /** - * @brief Compute with explicit quadrature settings. + * @brief Compute projector overlaps with explicit quadrature settings. * - * If calc_r is true, nlm[1..3] also store the Cartesian position moments. + * The ProjectorChannel radial convention is always r * p_l(r). Callers that + * own different physical projectors are responsible for passing that form. */ void snap_projector_half_tddft(const LCAO_Orbitals& orb, const std::vector& projector_channels, diff --git a/source/source_lcao/module_rt/snap_psibeta_half_tddft.cpp b/source/source_lcao/module_rt/snap_psibeta_half_tddft.cpp index 2b1cf62728f..db934c81ecc 100644 --- a/source/source_lcao/module_rt/snap_psibeta_half_tddft.cpp +++ b/source/source_lcao/module_rt/snap_psibeta_half_tddft.cpp @@ -37,7 +37,7 @@ void snap_psibeta_half_tddft(const LCAO_Orbitals& orb, std::vector channels; channels.reserve(infoNL_.nproj[T0]); - // Convert nonlocal pseudopotential beta projectors to the shared grid integrator input. + // UPF nonlocal beta projectors already follow the r * beta_l(r) convention. for (int ip = 0; ip < infoNL_.nproj[T0]; ++ip) { const auto& proj = infoNL_.Beta[T0].Proj[ip]; @@ -46,7 +46,7 @@ void snap_psibeta_half_tddft(const LCAO_Orbitals& orb, channel.mesh = proj.getNr(); channel.dk = proj.getDk(); channel.rcut = proj.getRcut(); - channel.radial_values = proj.getBeta_r(); + channel.radial_times_r = proj.getBeta_r(); channel.radial_grid = proj.getRadial(); channels.push_back(channel); } diff --git a/source/source_lcao/module_rt/snap_psibeta_half_tddft.h b/source/source_lcao/module_rt/snap_psibeta_half_tddft.h index 2644fbe6baf..164c40f79d6 100644 --- a/source/source_lcao/module_rt/snap_psibeta_half_tddft.h +++ b/source/source_lcao/module_rt/snap_psibeta_half_tddft.h @@ -14,6 +14,7 @@ namespace module_rt /** * @brief Compute RT-TDDFT velocity-gauge beta-projector overlaps. * + * UPF beta projectors are stored as r * beta_l(r) and selected by atom type T0. * This overload uses the production quadrature settings. */ void snap_psibeta_half_tddft(const LCAO_Orbitals& orb, diff --git a/source/source_lcao/module_rt/test/CMakeLists.txt b/source/source_lcao/module_rt/test/CMakeLists.txt index 7a2fb16c08d..cba91874465 100644 --- a/source/source_lcao/module_rt/test/CMakeLists.txt +++ b/source/source_lcao/module_rt/test/CMakeLists.txt @@ -36,4 +36,21 @@ AddTest( TARGET MODULE_LCAO_tddft_snap_psibeta_half_test LIBS parameter ${math_libs} base device orb numerical_atomic_orbitals tddft_test_lib SOURCES snap_psibeta_half_tddft_test.cpp ../snap_projector_half_tddft.cpp ../snap_psibeta_half_tddft.cpp + ../../center2_orb.cpp + ../../center2_orb-orb11.cpp + ../../center2_orb-orb21.cpp + ../../../source_cell/setup_nonlocal.cpp + ../../../source_cell/atom_spec.cpp + ../../../source_cell/atom_pseudo.cpp + ../../../source_cell/pseudo.cpp + ../../../source_cell/read_pp.cpp + ../../../source_cell/read_pp_complete.cpp + ../../../source_cell/read_pp_upf201.cpp + ../../../source_cell/read_pp_upf100.cpp + ../../../source_cell/read_pp_vwr.cpp + ../../../source_cell/read_pp_blps.cpp + ../../../source_io/module_hs/cal_r_overlap_R.cpp + ../../../source_io/module_hs/single_R_io.cpp + ../../../source_io/module_hs/rr_sparse_writer.cpp + ../../../source_pw/module_pwdft/soc.cpp ) diff --git a/source/source_lcao/module_rt/test/snap_psibeta_half_tddft_test.cpp b/source/source_lcao/module_rt/test/snap_psibeta_half_tddft_test.cpp index 3fc7f2b5d02..f8b2f8d48f4 100644 --- a/source/source_lcao/module_rt/test/snap_psibeta_half_tddft_test.cpp +++ b/source/source_lcao/module_rt/test/snap_psibeta_half_tddft_test.cpp @@ -1,9 +1,10 @@ #include "source_lcao/module_rt/snap_psibeta_half_tddft.h" #include "source_base/ylm.h" -#include "source_basis/module_nao/radial_collection.h" -#include "source_basis/module_nao/two_center_integrator.h" +#include "source_cell/read_pp.h" #include "source_cell/setup_nonlocal.h" +#include "source_cell/unitcell.h" +#include "source_io/module_hs/cal_r_overlap_R.h" #include #include @@ -12,25 +13,50 @@ #include #include -InfoNonlocal::InfoNonlocal() +SepPot::SepPot() = default; + +SepPot::~SepPot() { - this->Beta = new Numerical_Nonlocal[1]; - this->nproj = nullptr; - this->nprojmax = 0; - this->rcutmax_Beta = 0.0; + delete[] r; + delete[] rv; } -InfoNonlocal::~InfoNonlocal() +Sep_Cell::Sep_Cell() noexcept : ntype(0), omega(0.0), tpiba2(0.0) +{ +} + +Sep_Cell::~Sep_Cell() noexcept = default; + +Magnetism::Magnetism() { - delete[] this->Beta; - delete[] this->nproj; + tot_mag = 0.0; + abs_mag = 0.0; +} + +Magnetism::~Magnetism() +{ + delete[] start_mag; +} + +UnitCell::UnitCell() +{ + itia2iat.create(1, 1); +} + +UnitCell::~UnitCell() +{ + if (set_atom_flag) + { + delete[] atoms; + } } namespace { struct ComparisonStats { - double max_real_diff = 0.0; + double max_overlap_diff = 0.0; + double max_position_diff = 0.0; double max_imag_abs = 0.0; double max_reference_abs = 0.0; }; @@ -43,66 +69,78 @@ class SnapPsibetaHalfTddftTest : public ::testing::Test ModuleBase::Ylm::set_coefficients(); const std::string root = "../../../../../"; - const std::string orb_file = "tests/PP_ORB/C_gga_8au_100Ry_2s2p1d.orb"; - const std::string full_orb_file = root + orb_file; + const std::string orb_file = "tests/PP_ORB/Ti_gga_10au_100Ry_4s2p2d1f.orb"; const std::string orbital_files[1] = {orb_file}; std::ofstream ofs("snap_psibeta_half_tddft_test.log"); - orb.init(ofs, 1, root, orbital_files, "", 2, 100.0, 0.01, 0.01, 30.0, false, 0, false, false, 0); - - build_fake_beta_projectors(); + orb.init(ofs, 1, root, orbital_files, "", 3, 100.0, 0.01, 0.01, 30.0, false, 0, false, false, 0); - orb_radials.build(1, &full_orb_file, 'o'); - beta_radials.build(1, info_nl.Beta); + ASSERT_EQ(orb.Phi[0].getLmax(), 3); + ASSERT_EQ(orb.Phi[0].getNchi(0), 4); + ASSERT_EQ(orb.Phi[0].getNchi(1), 2); + ASSERT_EQ(orb.Phi[0].getNchi(2), 2); + ASSERT_EQ(orb.Phi[0].getNchi(3), 1); - const double rmax = std::max(orb_radials.rcut_max(), beta_radials.rcut_max()); - const double cutoff = 2.0 * rmax; - const int nr = static_cast(rmax / 0.01) + 1; - - orb_radials.set_uniform_grid(true, nr, cutoff, 'i', true); - beta_radials.set_uniform_grid(true, nr, cutoff, 'i', true); - overlap_orb_beta.tabulate(orb_radials, beta_radials, 'S', nr, cutoff); + build_ti_beta_projectors(root); + initialize_r_overlap_reference(); } - void build_fake_beta_projectors() + void build_ti_beta_projectors(const std::string& root) { - const int nproj = 2; - std::vector beta_lm(nproj); - - for (int iproj = 0; iproj < nproj; ++iproj) + ucell.ntype = 1; + ucell.nat = 1; + ucell.atoms = new Atom[1]; + ucell.set_atom_flag = true; + + Atom& atom = ucell.atoms[0]; + atom.label = "Ti"; + atom.type = 0; + atom.na = 1; + atom.nwl = orb.Phi[0].getLmax(); + atom.l_nchi.resize(atom.nwl + 1); + atom.nw = 0; + for (int L = 0; L <= atom.nwl; ++L) { - const int l = iproj; - const auto& phi_ln = orb.Phi[0].PhiLN(l, 0); - beta_lm[iproj].set_NL_proj("C", - 0, - l, - phi_ln.getNr(), - phi_ln.getRab(), - phi_ln.getRadial(), - phi_ln.getPsi_r(), - orb.get_kmesh(), - orb.get_dk(), - orb.get_dr_uniform()); + atom.l_nchi[L] = orb.Phi[0].getNchi(L); + atom.nw += (2 * L + 1) * atom.l_nchi[L]; } - - info_nl.nproj = new int[1]; - info_nl.nproj[0] = nproj; - info_nl.nprojmax = nproj; - info_nl.Beta[0].set_type_info(0, "C", "NC", 1, nproj, beta_lm.data()); - info_nl.rcutmax_Beta = info_nl.Beta[0].get_rcut_max(); + atom.tau.resize(1); + atom.tau[0] = ModuleBase::Vector3(0.0, 0.0, 0.0); + + Pseudopot_upf pseudo_reader; + std::string pseudo_type = "auto"; + const int pseudo_error = pseudo_reader.init_pseudo_reader(root + "tests/PP_ORB/Ti_ONCV_PBE-1.0.upf", pseudo_type, atom.ncpp); + ASSERT_EQ(pseudo_error, 0); + ASSERT_EQ(pseudo_type, "upf201"); + ASSERT_EQ(atom.ncpp.psd, "Ti"); + ASSERT_EQ(atom.ncpp.pp_type, "NC"); + ASSERT_EQ(atom.ncpp.nbeta, 6); + ASSERT_EQ(atom.ncpp.lll, std::vector({0, 0, 1, 1, 2, 2})); + pseudo_reader.complete_default(atom.ncpp); + ASSERT_EQ(atom.ncpp.nh, 18); + ASSERT_EQ(atom.ncpp.jjj.size(), 6); + + ucell.infoNL.nproj = new int[1]; + std::ofstream log("snap_psibeta_half_tddft_nonlocal.log"); + ucell.infoNL.Set_NonLocal(0, &atom, ucell.infoNL.nproj[0], orb.get_kmesh(), orb.get_dk(), orb.get_dr_uniform(), log); + + ASSERT_EQ(ucell.infoNL.nproj[0], 6); + ucell.infoNL.nprojmax = ucell.infoNL.nproj[0]; + ucell.infoNL.rcutmax_Beta = ucell.infoNL.Beta[0].get_rcut_max(); } - static int abacus_m_to_m(const int m) + void initialize_r_overlap_reference() { - return (m % 2 == 0) ? -m / 2 : (m + 1) / 2; + r_calculator.init_nonlocal(ucell, pv, orb); } - ComparisonStats compare_zero_vector_potential(const int lebedev_grid_points) + ComparisonStats compare_zero_vector_potential(const int radial_grid_num, const int lebedev_grid_points) { const ModuleBase::Vector3 R0(0.1, -0.2, 0.3); const ModuleBase::Vector3 R1(0.4, 0.2, -0.1); const ModuleBase::Vector3 zero_A(0.0, 0.0, 0.0); module_rt::SnapIntegrationOptions options; + options.radial_grid_num = radial_grid_num; options.lebedev_grid_points = lebedev_grid_points; ComparisonStats stats; @@ -114,41 +152,45 @@ class SnapPsibetaHalfTddftTest : public ::testing::Test for (int m1 = 0; m1 < 2 * L1 + 1; ++m1) { std::vector>> grid_nlm; - module_rt::snap_psibeta_half_tddft(orb, - info_nl, - grid_nlm, - R1, - 0, - L1, - m1, - N1, - R0, - 0, - zero_A, - false, - options); - - std::vector> tci_nlm; - overlap_orb_beta.snap(0, L1, N1, abacus_m_to_m(m1), 0, R0 - R1, false, tci_nlm); - - EXPECT_FALSE(grid_nlm.empty()); - EXPECT_FALSE(tci_nlm.empty()); - if (grid_nlm.empty() || tci_nlm.empty()) + module_rt::snap_psibeta_half_tddft(orb, ucell.infoNL, grid_nlm, R1, 0, L1, m1, N1, R0, 0, zero_A, true, options); + + std::vector> reference_nlm; + r_calculator.get_psi_r_beta(ucell, reference_nlm, R1, 0, L1, m1, N1, R0, 0); + + EXPECT_EQ(grid_nlm.size(), 4); + EXPECT_EQ(reference_nlm.size(), 4); + if (grid_nlm.size() != 4 || reference_nlm.size() != 4) { continue; } - EXPECT_EQ(grid_nlm[0].size(), tci_nlm[0].size()); - if (grid_nlm[0].size() != tci_nlm[0].size()) + + bool sizes_match = true; + for (size_t dim = 0; dim < grid_nlm.size(); ++dim) + { + EXPECT_EQ(grid_nlm[dim].size(), reference_nlm[dim].size()); + sizes_match = sizes_match && (grid_nlm[dim].size() == reference_nlm[dim].size()); + } + if (!sizes_match) { continue; } - for (size_t i = 0; i < grid_nlm[0].size(); ++i) + for (size_t dim = 0; dim < grid_nlm.size(); ++dim) { - stats.max_real_diff - = std::max(stats.max_real_diff, std::abs(grid_nlm[0][i].real() - tci_nlm[0][i])); - stats.max_imag_abs = std::max(stats.max_imag_abs, std::abs(grid_nlm[0][i].imag())); - stats.max_reference_abs = std::max(stats.max_reference_abs, std::abs(tci_nlm[0][i])); + for (size_t i = 0; i < grid_nlm[dim].size(); ++i) + { + const double real_diff = std::abs(grid_nlm[dim][i].real() - reference_nlm[dim][i]); + if (dim == 0) + { + stats.max_overlap_diff = std::max(stats.max_overlap_diff, real_diff); + } + else + { + stats.max_position_diff = std::max(stats.max_position_diff, real_diff); + } + stats.max_imag_abs = std::max(stats.max_imag_abs, std::abs(grid_nlm[dim][i].imag())); + stats.max_reference_abs = std::max(stats.max_reference_abs, std::abs(reference_nlm[dim][i])); + } } } } @@ -158,50 +200,44 @@ class SnapPsibetaHalfTddftTest : public ::testing::Test } LCAO_Orbitals orb; - InfoNonlocal info_nl; - RadialCollection orb_radials; - RadialCollection beta_radials; - TwoCenterIntegrator overlap_orb_beta; + UnitCell ucell; + Parallel_Orbitals pv; + cal_r_overlap_R r_calculator; }; } // namespace TEST_F(SnapPsibetaHalfTddftTest, ZeroVectorPotentialMatchesTwoCenterIntegral) { - const double real_tolerance = 5.0e-8; + const double overlap_tolerance = 4.0e-7; + const double position_tolerance = 6.0e-7; const double imag_tolerance = 1.0e-12; - const ComparisonStats stats = compare_zero_vector_potential(110); + const ComparisonStats stats = compare_zero_vector_potential(140, 110); - EXPECT_LT(stats.max_real_diff, real_tolerance) << "max reference abs = " << stats.max_reference_abs; + EXPECT_LT(stats.max_overlap_diff, overlap_tolerance) << "max reference abs = " << stats.max_reference_abs; + EXPECT_LT(stats.max_position_diff, position_tolerance) << "max reference abs = " << stats.max_reference_abs; EXPECT_LT(stats.max_imag_abs, imag_tolerance) << "max reference abs = " << stats.max_reference_abs; } -TEST_F(SnapPsibetaHalfTddftTest, ZeroVectorPotentialHighOrderGridMatchesTwoCenterIntegral) +TEST_F(SnapPsibetaHalfTddftTest, ZeroVectorPotentialDenseRadialGridMatchesTwoCenterIntegral) { - const double real_tolerance = 5.0e-8; + const double overlap_tolerance = 3.0e-7; + const double position_tolerance = 5.0e-7; const double imag_tolerance = 1.0e-12; - const ComparisonStats stats = compare_zero_vector_potential(590); + const ComparisonStats stats = compare_zero_vector_potential(280, 110); - EXPECT_LT(stats.max_real_diff, real_tolerance) << "max reference abs = " << stats.max_reference_abs; + EXPECT_LT(stats.max_overlap_diff, overlap_tolerance) << "max reference abs = " << stats.max_reference_abs; + EXPECT_LT(stats.max_position_diff, position_tolerance) << "max reference abs = " << stats.max_reference_abs; EXPECT_LT(stats.max_imag_abs, imag_tolerance) << "max reference abs = " << stats.max_reference_abs; } -TEST_F(SnapPsibetaHalfTddftTest, ZeroVectorPotentialPositionMomentsAreReal) +TEST_F(SnapPsibetaHalfTddftTest, ZeroVectorPotentialHighOrderGridMatchesTwoCenterIntegral) { - const ModuleBase::Vector3 R0(-0.3, 0.2, 0.1); - const ModuleBase::Vector3 R1(0.2, -0.1, 0.4); - const ModuleBase::Vector3 zero_A(0.0, 0.0, 0.0); - const double tolerance = 1.0e-12; - - std::vector>> nlm; - module_rt::snap_psibeta_half_tddft(orb, info_nl, nlm, R1, 0, 1, 1, 0, R0, 0, zero_A, true); + const double overlap_tolerance = 4.0e-7; + const double position_tolerance = 6.0e-7; + const double imag_tolerance = 1.0e-12; + const ComparisonStats stats = compare_zero_vector_potential(140, 590); - ASSERT_EQ(nlm.size(), 4); - for (const auto& dim: nlm) - { - ASSERT_EQ(dim.size(), 4); - for (const std::complex& value: dim) - { - EXPECT_NEAR(value.imag(), 0.0, tolerance); - } - } + EXPECT_LT(stats.max_overlap_diff, overlap_tolerance) << "max reference abs = " << stats.max_reference_abs; + EXPECT_LT(stats.max_position_diff, position_tolerance) << "max reference abs = " << stats.max_reference_abs; + EXPECT_LT(stats.max_imag_abs, imag_tolerance) << "max reference abs = " << stats.max_reference_abs; } From 358403eada70ceff61ac579e1f5143b8b7862e5c Mon Sep 17 00:00:00 2001 From: AsTonyshment Date: Fri, 3 Jul 2026 16:46:32 +0800 Subject: [PATCH 2/2] Style: Format cal_r_overlap_R --- .../source_io/module_hs/cal_r_overlap_R.cpp | 320 ++++++++---------- source/source_io/module_hs/cal_r_overlap_R.h | 96 +++--- 2 files changed, 180 insertions(+), 236 deletions(-) diff --git a/source/source_io/module_hs/cal_r_overlap_R.cpp b/source/source_io/module_hs/cal_r_overlap_R.cpp index a1e821e7268..052a6760c0a 100644 --- a/source/source_io/module_hs/cal_r_overlap_R.cpp +++ b/source/source_io/module_hs/cal_r_overlap_R.cpp @@ -2,12 +2,12 @@ #include "rr_sparse_writer.h" #include "single_R_io.h" -#include "source_io/module_parameter/parameter.h" +#include "source_base/mathzone_add1.h" #include "source_base/parallel_reduce.h" #include "source_base/timer.h" #include "source_base/tool_quit.h" #include "source_cell/module_neighbor/sltk_grid_driver.h" -#include "source_base/mathzone_add1.h" +#include "source_io/module_parameter/parameter.h" cal_r_overlap_R::cal_r_overlap_R() { @@ -17,8 +17,7 @@ cal_r_overlap_R::~cal_r_overlap_R() { } -void cal_r_overlap_R::initialize_orb_table(const UnitCell& ucell, - const LCAO_Orbitals& orb) +void cal_r_overlap_R::initialize_orb_table(const UnitCell& ucell, const LCAO_Orbitals& orb) { const int ntype = orb.get_ntype(); int lmax_orb = -1; @@ -34,19 +33,13 @@ void cal_r_overlap_R::initialize_orb_table(const UnitCell& ucell, const int Lmax = lmax_orb + 1; const int Lmax_used = 2 * lmax_orb + 1; - Center2_Orb::init_Table_Spherical_Bessel(Lmax_used, - dr, - dk, - kmesh, - Rmesh, - psb_); + Center2_Orb::init_Table_Spherical_Bessel(Lmax_used, dr, dk, kmesh, Rmesh, psb_); ModuleBase::Ylm::set_coefficients(); MGT.init_Gaunt_CH(Lmax); MGT.init_Gaunt(Lmax); } -void cal_r_overlap_R::construct_orbs_and_orb_r(const UnitCell& ucell, - const LCAO_Orbitals& orb) +void cal_r_overlap_R::construct_orbs_and_orb_r(const UnitCell& ucell, const LCAO_Orbitals& orb) { int orb_r_ntype = 0; int mat_Nr = orb.Phi[0].PhiLN(0, 0).getNr(); @@ -137,9 +130,8 @@ void cal_r_overlap_R::construct_orbs_and_orb_r(const UnitCell& ucell, { for (int NB = 0; NB < orb.Phi[TB].getNchi(LB); ++NB) { - center2_orb21_r[TA][TB][LA][NA][LB].insert(std::make_pair( - NB, - Center2_Orb::Orb21(orbs[TA][LA][NA], orb_r, orbs[TB][LB][NB], psb_, MGT))); + center2_orb21_r[TA][TB][LA][NA][LB].insert( + std::make_pair(NB, Center2_Orb::Orb21(orbs[TA][LA][NA], orb_r, orbs[TB][LB][NB], psb_, MGT))); } } } @@ -225,7 +217,7 @@ void cal_r_overlap_R::construct_orbs_and_orb_r(const UnitCell& ucell, } } -void cal_r_overlap_R::construct_orbs_and_nonlocal_and_orb_r(const UnitCell& ucell,const LCAO_Orbitals& orb) +void cal_r_overlap_R::construct_orbs_and_nonlocal_and_orb_r(const UnitCell& ucell, const LCAO_Orbitals& orb) { const InfoNonlocal& infoNL_ = ucell.infoNL; @@ -294,37 +286,42 @@ void cal_r_overlap_R::construct_orbs_and_nonlocal_and_orb_r(const UnitCell& ucel { int nr = infoNL_.Beta[T].Proj[ip].getNr(); double dr_uniform = 0.01; - int nr_uniform = static_cast((infoNL_.Beta[T].Proj[ip].getRadial(nr-1) - infoNL_.Beta[T].Proj[ip].getRadial(0))/dr_uniform) + 1; + int nr_uniform + = static_cast((infoNL_.Beta[T].Proj[ip].getRadial(nr - 1) - infoNL_.Beta[T].Proj[ip].getRadial(0)) / dr_uniform) + 1; double* rad = new double[nr_uniform]; double* rab = new double[nr_uniform]; for (int ir = 0; ir < nr_uniform; ir++) { - rad[ir] = ir*dr_uniform; + rad[ir] = ir * dr_uniform; rab[ir] = dr_uniform; } double* y2 = new double[nr]; double* Beta_r_uniform = new double[nr_uniform]; double* dbeta_uniform = new double[nr_uniform]; - ModuleBase::Mathzone_Add1::SplineD2(infoNL_.Beta[T].Proj[ip].getRadial(), infoNL_.Beta[T].Proj[ip].getBeta_r(), nr, 0.0, 0.0, y2); - ModuleBase::Mathzone_Add1::Cubic_Spline_Interpolation( - infoNL_.Beta[T].Proj[ip].getRadial(), - infoNL_.Beta[T].Proj[ip].getBeta_r(), - y2, - nr, - rad, - nr_uniform, - Beta_r_uniform, - dbeta_uniform - ); + ModuleBase::Mathzone_Add1::SplineD2(infoNL_.Beta[T].Proj[ip].getRadial(), + infoNL_.Beta[T].Proj[ip].getBeta_r(), + nr, + 0.0, + 0.0, + y2); + ModuleBase::Mathzone_Add1::Cubic_Spline_Interpolation(infoNL_.Beta[T].Proj[ip].getRadial(), + infoNL_.Beta[T].Proj[ip].getBeta_r(), + y2, + nr, + rad, + nr_uniform, + Beta_r_uniform, + dbeta_uniform); // linear extrapolation at the zero point if (infoNL_.Beta[T].Proj[ip].getRadial(0) > 1e-10) { - double slope = (infoNL_.Beta[T].Proj[ip].getBeta_r(1) - infoNL_.Beta[T].Proj[ip].getBeta_r(0)) / (infoNL_.Beta[T].Proj[ip].getRadial(1) - infoNL_.Beta[T].Proj[ip].getRadial(0)); + double slope = (infoNL_.Beta[T].Proj[ip].getBeta_r(1) - infoNL_.Beta[T].Proj[ip].getBeta_r(0)) + / (infoNL_.Beta[T].Proj[ip].getRadial(1) - infoNL_.Beta[T].Proj[ip].getRadial(0)); Beta_r_uniform[0] = infoNL_.Beta[T].Proj[ip].getBeta_r(0) - slope * infoNL_.Beta[T].Proj[ip].getRadial(0); } - // Here, the operation beta_r / r is performed. To avoid divergence at r=0, beta_r(0) is set to beta_r(1). + // Here, the operation beta_r / r is performed. To avoid divergence at r=0, beta_r(0) is set to beta_r(1). // However, this may introduce issues, so caution is needed. for (int ir = 1; ir < nr_uniform; ir++) { @@ -348,11 +345,11 @@ void cal_r_overlap_R::construct_orbs_and_nonlocal_and_orb_r(const UnitCell& ucel true, PARAM.inp.cal_force); - delete [] rad; - delete [] rab; - delete [] y2; - delete [] Beta_r_uniform; - delete [] dbeta_uniform; + delete[] rad; + delete[] rab; + delete[] y2; + delete[] Beta_r_uniform; + delete[] dbeta_uniform; } } @@ -384,9 +381,8 @@ void cal_r_overlap_R::construct_orbs_and_nonlocal_and_orb_r(const UnitCell& ucel { for (int ip = 0; ip < infoNL_.nproj[TB]; ip++) { - center2_orb21_r_nonlocal[TA][TB][LA][NA].insert(std::make_pair( - ip, - Center2_Orb::Orb21(orbs[TA][LA][NA], orb_r, orbs_nonlocal[TB][ip], psb_, MGT))); + center2_orb21_r_nonlocal[TA][TB][LA][NA].insert( + std::make_pair(ip, Center2_Orb::Orb21(orbs[TA][LA][NA], orb_r, orbs_nonlocal[TB][ip], psb_, MGT))); } } } @@ -465,27 +461,27 @@ void cal_r_overlap_R::construct_orbs_and_nonlocal_and_orb_r(const UnitCell& ucel } } -void cal_r_overlap_R::init(const UnitCell& ucell,const Parallel_Orbitals& pv, const LCAO_Orbitals& orb) +void cal_r_overlap_R::init(const UnitCell& ucell, const Parallel_Orbitals& pv, const LCAO_Orbitals& orb) { ModuleBase::TITLE("cal_r_overlap_R", "init"); ModuleBase::timer::start("cal_r_overlap_R", "init"); this->ParaV = &pv; - initialize_orb_table(ucell,orb); - construct_orbs_and_orb_r(ucell,orb); + initialize_orb_table(ucell, orb); + construct_orbs_and_orb_r(ucell, orb); ModuleBase::timer::end("cal_r_overlap_R", "init"); return; } -void cal_r_overlap_R::init_nonlocal(const UnitCell& ucell,const Parallel_Orbitals& pv, const LCAO_Orbitals& orb) +void cal_r_overlap_R::init_nonlocal(const UnitCell& ucell, const Parallel_Orbitals& pv, const LCAO_Orbitals& orb) { ModuleBase::TITLE("cal_r_overlap_R", "init_nonlocal"); ModuleBase::timer::start("cal_r_overlap_R", "init_nonlocal"); this->ParaV = &pv; - initialize_orb_table(ucell,orb); - construct_orbs_and_nonlocal_and_orb_r(ucell,orb); + initialize_orb_table(ucell, orb); + construct_orbs_and_nonlocal_and_orb_r(ucell, orb); ModuleBase::timer::end("cal_r_overlap_R", "init_nonlocal"); return; @@ -508,44 +504,31 @@ ModuleBase::Vector3 cal_r_overlap_R::get_psi_r_psi(const ModuleBase::Vec double overlap_o = center2_orb11[T1][T2][L1][N1][L2].at(N2).cal_overlap(origin_point, distance, m1, m2); - double overlap_x = -1 * factor - * center2_orb21_r[T1][T2][L1][N1][L2].at(N2).cal_overlap(origin_point, - distance, - m1, - 1, - m2); // m = 1 - - double overlap_y = -1 * factor - * center2_orb21_r[T1][T2][L1][N1][L2].at(N2).cal_overlap(origin_point, - distance, - m1, - 2, - m2); // m = -1 - - double overlap_z = factor - * center2_orb21_r[T1][T2][L1][N1][L2].at(N2).cal_overlap(origin_point, - distance, - m1, - 0, - m2); // m = 0 - - ModuleBase::Vector3 temp_prp - = ModuleBase::Vector3(overlap_x, overlap_y, overlap_z) + R1 * overlap_o; + double overlap_x = -1 * factor * center2_orb21_r[T1][T2][L1][N1][L2].at(N2).cal_overlap(origin_point, distance, m1, 1, + m2); // m = 1 + + double overlap_y = -1 * factor * center2_orb21_r[T1][T2][L1][N1][L2].at(N2).cal_overlap(origin_point, distance, m1, 2, + m2); // m = -1 + + double overlap_z = factor * center2_orb21_r[T1][T2][L1][N1][L2].at(N2).cal_overlap(origin_point, distance, m1, 0, + m2); // m = 0 + + ModuleBase::Vector3 temp_prp = ModuleBase::Vector3(overlap_x, overlap_y, overlap_z) + R1 * overlap_o; return temp_prp; } ModuleBase::Vector3 cal_r_overlap_R::get_psi_r_gradpsi(const ModuleBase::Vector3& R1, - const int& T1, - const int& L1, - const int& m1, - const int& N1, - const ModuleBase::Vector3& R2, - const int& T2, - const int& L2, - const int& m2, - const int& N2, - const ModuleBase::Vector3& Efield, - const ModuleBase::Vector3& dR) + const int& T1, + const int& L1, + const int& m1, + const int& N1, + const ModuleBase::Vector3& R2, + const int& T2, + const int& L2, + const int& m2, + const int& N2, + const ModuleBase::Vector3& Efield, + const ModuleBase::Vector3& dR) { ModuleBase::Vector3 origin_point(0.0, 0.0, 0.0); double factor = sqrt(ModuleBase::FOUR_PI / 3.0); @@ -553,25 +536,28 @@ ModuleBase::Vector3 cal_r_overlap_R::get_psi_r_gradpsi(const ModuleBase: ModuleBase::Vector3 grad_o = center2_orb11[T1][T2][L1][N1][L2].at(N2).cal_grad_overlap(origin_point, distance, m1, m2); - ModuleBase::Vector3 grad_rx = -1 * factor * center2_orb21_r[T1][T2][L1][N1][L2].at(N2).cal_grad_overlap(origin_point, + ModuleBase::Vector3 grad_rx = -1 * factor + * center2_orb21_r[T1][T2][L1][N1][L2].at(N2).cal_grad_overlap(origin_point, distance, m1, 1, m2); // m = 1 - ModuleBase::Vector3 grad_ry = -1 * factor * center2_orb21_r[T1][T2][L1][N1][L2].at(N2).cal_grad_overlap(origin_point, + ModuleBase::Vector3 grad_ry = -1 * factor + * center2_orb21_r[T1][T2][L1][N1][L2].at(N2).cal_grad_overlap(origin_point, distance, m1, 2, m2); // m = -1 - ModuleBase::Vector3 grad_rz = factor * center2_orb21_r[T1][T2][L1][N1][L2].at(N2).cal_grad_overlap(origin_point, + ModuleBase::Vector3 grad_rz = factor + * center2_orb21_r[T1][T2][L1][N1][L2].at(N2).cal_grad_overlap(origin_point, distance, m1, 0, m2); // m = 0 - ModuleBase::Vector3 temp_prp = Efield[0] * grad_rx + Efield[1] * grad_ry + Efield[2] * grad_rz + (Efield*(R1-dR)) * grad_o; + ModuleBase::Vector3 temp_prp = Efield[0] * grad_rx + Efield[1] * grad_ry + Efield[2] * grad_rz + (Efield * (R1 - dR)) * grad_o; return temp_prp; } @@ -593,7 +579,7 @@ void cal_r_overlap_R::get_psi_r_beta(const UnitCell& ucell, nlm.resize(4); if (nproj == 0) { - for(int i = 0;i < 4;i++) + for (int i = 0; i < 4; i++) { nlm[i].resize(1); } @@ -606,7 +592,7 @@ void cal_r_overlap_R::get_psi_r_beta(const UnitCell& ucell, const int L2 = infoNL_.Beta[T2].Proj[ip].getL(); // mohan add 2021-05-07 natomwfc += 2 * L2 + 1; } - for(int i = 0;i < 4;i++) + for (int i = 0; i < 4; i++) { nlm[i].resize(natomwfc); } @@ -616,8 +602,7 @@ void cal_r_overlap_R::get_psi_r_beta(const UnitCell& ucell, const int L2 = infoNL_.Beta[T2].Proj[ip].getL(); for (int m2 = 0; m2 < 2 * L2 + 1; m2++) { - double overlap_o - = center2_orb11_nonlocal[T1][T2][L1][N1].at(ip).cal_overlap(origin_point, distance, m1, m2); + double overlap_o = center2_orb11_nonlocal[T1][T2][L1][N1].at(ip).cal_overlap(origin_point, distance, m1, m2); double overlap_x = -1 * factor * center2_orb21_r_nonlocal[T1][T2][L1][N1].at(ip).cal_overlap(origin_point, @@ -640,9 +625,9 @@ void cal_r_overlap_R::get_psi_r_beta(const UnitCell& ucell, 0, m2); // m = 0 - //nlm[index] = ModuleBase::Vector3(overlap_x, overlap_y, overlap_z) + R1 * overlap_o; + // nlm[index] = ModuleBase::Vector3(overlap_x, overlap_y, overlap_z) + R1 * overlap_o; - //nlm[index] = ModuleBase::Vector3(overlap_o, overlap_y, overlap_z);// + R1 * overlap_o; + // nlm[index] = ModuleBase::Vector3(overlap_o, overlap_y, overlap_z);// + R1 * overlap_o; nlm[0][index] = overlap_o; nlm[1][index] = overlap_x + (R1 * overlap_o).x; nlm[2][index] = overlap_y + (R1 * overlap_o).y; @@ -652,7 +637,6 @@ void cal_r_overlap_R::get_psi_r_beta(const UnitCell& ucell, } } - void cal_r_overlap_R::out_rR(const UnitCell& ucell, const Grid_Driver& gd, const int& istep, const int precision) { ModuleBase::TITLE("cal_r_overlap_R", "out_rR"); @@ -709,8 +693,7 @@ void cal_r_overlap_R::out_rR(const UnitCell& ucell, const Grid_Driver& gd, const } if (!ofs_tem1.is_open()) { - ModuleBase::WARNING_QUIT("cal_r_overlap_R::out_rR", - "Cannot open temporary sparse matrix file: " + tem1.str()); + ModuleBase::WARNING_QUIT("cal_r_overlap_R::out_rR", "Cannot open temporary sparse matrix file: " + tem1.str()); } } @@ -739,8 +722,8 @@ void cal_r_overlap_R::out_rR(const UnitCell& ucell, const Grid_Driver& gd, const int orb_index_col = iw2 / PARAM.globalv.npol; // The off-diagonal term in SOC calculaiton is zero, and the two diagonal terms are the same - int new_index = iw1 - PARAM.globalv.npol * orb_index_row - + (iw2 - PARAM.globalv.npol * orb_index_col) * PARAM.globalv.npol; + int new_index + = iw1 - PARAM.globalv.npol * orb_index_row + (iw2 - PARAM.globalv.npol * orb_index_col) * PARAM.globalv.npol; if (new_index == 0 || new_index == 3) { @@ -757,41 +740,34 @@ void cal_r_overlap_R::out_rR(const UnitCell& ucell, const Grid_Driver& gd, const int im2 = iw2im[orb_index_col]; ModuleBase::Vector3 r_distance - = (ucell.atoms[it2].tau[ia2] - ucell.atoms[it1].tau[ia1] + R_car) - * ucell.lat0; - - double overlap_o = center2_orb11[it1][it2][iL1][iN1][iL2].at(iN2).cal_overlap(origin_point, - r_distance, - im1, - im2); - - double overlap_x - = -1 * factor - * center2_orb21_r[it1][it2][iL1][iN1][iL2].at(iN2).cal_overlap(origin_point, - r_distance, - im1, - 1, - im2); // m = 1 - - double overlap_y - = -1 * factor - * center2_orb21_r[it1][it2][iL1][iN1][iL2].at(iN2).cal_overlap(origin_point, - r_distance, - im1, - 2, - im2); // m = -1 - - double overlap_z - = factor - * center2_orb21_r[it1][it2][iL1][iN1][iL2].at(iN2).cal_overlap(origin_point, - r_distance, - im1, - 0, - im2); // m = 0 - - ModuleBase::Vector3 temp_prp - = ModuleBase::Vector3(overlap_x, overlap_y, overlap_z) - + ucell.atoms[it1].tau[ia1] * ucell.lat0 * overlap_o; + = (ucell.atoms[it2].tau[ia2] - ucell.atoms[it1].tau[ia1] + R_car) * ucell.lat0; + + double overlap_o + = center2_orb11[it1][it2][iL1][iN1][iL2].at(iN2).cal_overlap(origin_point, r_distance, im1, im2); + + double overlap_x = -1 * factor + * center2_orb21_r[it1][it2][iL1][iN1][iL2].at(iN2).cal_overlap(origin_point, + r_distance, + im1, + 1, + im2); // m = 1 + + double overlap_y = -1 * factor + * center2_orb21_r[it1][it2][iL1][iN1][iL2].at(iN2).cal_overlap(origin_point, + r_distance, + im1, + 2, + im2); // m = -1 + + double overlap_z = factor + * center2_orb21_r[it1][it2][iL1][iN1][iL2].at(iN2).cal_overlap(origin_point, + r_distance, + im1, + 0, + im2); // m = 0 + + ModuleBase::Vector3 temp_prp = ModuleBase::Vector3(overlap_x, overlap_y, overlap_z) + + ucell.atoms[it1].tau[ia1] * ucell.lat0 * overlap_o; if (std::abs(temp_prp.x) > sparse_threshold) { @@ -858,10 +834,7 @@ void cal_r_overlap_R::out_rR(const UnitCell& ucell, const Grid_Driver& gd, const if (rR_nonzero_num[direction]) { - ModuleIO::output_single_R(ofs_tem1, - psi_r_psi_sparse[direction], - *(this->ParaV), - single_R_options); + ModuleIO::output_single_R(ofs_tem1, psi_r_psi_sparse[direction], *(this->ParaV), single_R_options); } else { @@ -876,8 +849,7 @@ void cal_r_overlap_R::out_rR(const UnitCell& ucell, const Grid_Driver& gd, const std::stringstream ssr; if (PARAM.inp.calculation == "md" && !PARAM.inp.out_app_flag) { - ssr << PARAM.globalv.global_matrix_dir - << "rrg" << step << ".csr"; + ssr << PARAM.globalv.global_matrix_dir << "rrg" << step << ".csr"; } else { @@ -937,16 +909,14 @@ void cal_r_overlap_R::out_rR_other(const UnitCell& ucell, } if (!ofs_tem1.is_open()) { - ModuleBase::WARNING_QUIT("cal_r_overlap_R::out_rR_other", - "Cannot open temporary sparse matrix file: " + tem1.str()); + ModuleBase::WARNING_QUIT("cal_r_overlap_R::out_rR_other", "Cannot open temporary sparse matrix file: " + tem1.str()); } } std::stringstream ssr; if (PARAM.inp.calculation == "md" && !PARAM.inp.out_app_flag) { - ssr << PARAM.globalv.global_matrix_dir - << "rrg" << step << ".csr"; + ssr << PARAM.globalv.global_matrix_dir << "rrg" << step << ".csr"; } else { @@ -979,8 +949,8 @@ void cal_r_overlap_R::out_rR_other(const UnitCell& ucell, int orb_index_col = iw2 / PARAM.globalv.npol; // The off-diagonal term in SOC calculaiton is zero, and the two diagonal terms are the same - int new_index = iw1 - PARAM.globalv.npol * orb_index_row - + (iw2 - PARAM.globalv.npol * orb_index_col) * PARAM.globalv.npol; + int new_index + = iw1 - PARAM.globalv.npol * orb_index_row + (iw2 - PARAM.globalv.npol * orb_index_col) * PARAM.globalv.npol; if (new_index == 0 || new_index == 3) { @@ -997,41 +967,34 @@ void cal_r_overlap_R::out_rR_other(const UnitCell& ucell, int im2 = iw2im[orb_index_col]; ModuleBase::Vector3 r_distance - = (ucell.atoms[it2].tau[ia2] - ucell.atoms[it1].tau[ia1] + R_car) - * ucell.lat0; - - double overlap_o = center2_orb11[it1][it2][iL1][iN1][iL2].at(iN2).cal_overlap(origin_point, - r_distance, - im1, - im2); - - double overlap_x - = -1 * factor - * center2_orb21_r[it1][it2][iL1][iN1][iL2].at(iN2).cal_overlap(origin_point, - r_distance, - im1, - 1, - im2); // m = 1 - - double overlap_y - = -1 * factor - * center2_orb21_r[it1][it2][iL1][iN1][iL2].at(iN2).cal_overlap(origin_point, - r_distance, - im1, - 2, - im2); // m = -1 - - double overlap_z - = factor - * center2_orb21_r[it1][it2][iL1][iN1][iL2].at(iN2).cal_overlap(origin_point, - r_distance, - im1, - 0, - im2); // m = 0 - - ModuleBase::Vector3 temp_prp - = ModuleBase::Vector3(overlap_x, overlap_y, overlap_z) - + ucell.atoms[it1].tau[ia1] * ucell.lat0 * overlap_o; + = (ucell.atoms[it2].tau[ia2] - ucell.atoms[it1].tau[ia1] + R_car) * ucell.lat0; + + double overlap_o + = center2_orb11[it1][it2][iL1][iN1][iL2].at(iN2).cal_overlap(origin_point, r_distance, im1, im2); + + double overlap_x = -1 * factor + * center2_orb21_r[it1][it2][iL1][iN1][iL2].at(iN2).cal_overlap(origin_point, + r_distance, + im1, + 1, + im2); // m = 1 + + double overlap_y = -1 * factor + * center2_orb21_r[it1][it2][iL1][iN1][iL2].at(iN2).cal_overlap(origin_point, + r_distance, + im1, + 2, + im2); // m = -1 + + double overlap_z = factor + * center2_orb21_r[it1][it2][iL1][iN1][iL2].at(iN2).cal_overlap(origin_point, + r_distance, + im1, + 0, + im2); // m = 0 + + ModuleBase::Vector3 temp_prp = ModuleBase::Vector3(overlap_x, overlap_y, overlap_z) + + ucell.atoms[it1].tau[ia1] * ucell.lat0 * overlap_o; if (std::abs(temp_prp.x) > sparse_threshold) { @@ -1100,10 +1063,7 @@ void cal_r_overlap_R::out_rR_other(const UnitCell& ucell, if (rR_nonzero_num[direction]) { - ModuleIO::output_single_R(ofs_tem1, - psi_r_psi_sparse[direction], - *(this->ParaV), - single_R_options); + ModuleIO::output_single_R(ofs_tem1, psi_r_psi_sparse[direction], *(this->ParaV), single_R_options); } else { diff --git a/source/source_io/module_hs/cal_r_overlap_R.h b/source/source_io/module_hs/cal_r_overlap_R.h index 543bf248fa7..1ac5c0f47f1 100644 --- a/source/source_io/module_hs/cal_r_overlap_R.h +++ b/source/source_io/module_hs/cal_r_overlap_R.h @@ -1,7 +1,6 @@ #ifndef CAL_R_OVERLAP_R_H #define CAL_R_OVERLAP_R_H -#include "source_lcao/module_ri/abfs-vector3_order.h" #include "source_base/sph_bessel_recursive.h" #include "source_base/vector3.h" #include "source_base/ylm.h" @@ -14,6 +13,7 @@ #include "source_lcao/center2_orb-orb11.h" #include "source_lcao/center2_orb-orb21.h" #include "source_lcao/center2_orb.h" +#include "source_lcao/module_ri/abfs-vector3_order.h" #include #include @@ -31,45 +31,39 @@ class cal_r_overlap_R double sparse_threshold = 1e-10; bool binary = false; - void init(const UnitCell& ucell,const Parallel_Orbitals& pv, const LCAO_Orbitals& orb); - void init_nonlocal(const UnitCell& ucell,const Parallel_Orbitals& pv, const LCAO_Orbitals& orb); - ModuleBase::Vector3 get_psi_r_psi( - const ModuleBase::Vector3& R1, - const int& T1, - const int& L1, - const int& m1, - const int& N1, - const ModuleBase::Vector3& R2, - const int& T2, - const int& L2, - const int& m2, - const int& N2 - ); - ModuleBase::Vector3 get_psi_r_gradpsi( - const ModuleBase::Vector3& R1, - const int& T1, - const int& L1, - const int& m1, - const int& N1, - const ModuleBase::Vector3& R2, - const int& T2, - const int& L2, - const int& m2, - const int& N2, - const ModuleBase::Vector3& Efield, - const ModuleBase::Vector3& dR - ); - void get_psi_r_beta( - const UnitCell& ucell, - std::vector>& nlm, - const ModuleBase::Vector3& R1, - const int& T1, - const int& L1, - const int& m1, - const int& N1, - const ModuleBase::Vector3& R2, - const int& T2 - ); + void init(const UnitCell& ucell, const Parallel_Orbitals& pv, const LCAO_Orbitals& orb); + void init_nonlocal(const UnitCell& ucell, const Parallel_Orbitals& pv, const LCAO_Orbitals& orb); + ModuleBase::Vector3 get_psi_r_psi(const ModuleBase::Vector3& R1, + const int& T1, + const int& L1, + const int& m1, + const int& N1, + const ModuleBase::Vector3& R2, + const int& T2, + const int& L2, + const int& m2, + const int& N2); + ModuleBase::Vector3 get_psi_r_gradpsi(const ModuleBase::Vector3& R1, + const int& T1, + const int& L1, + const int& m1, + const int& N1, + const ModuleBase::Vector3& R2, + const int& T2, + const int& L2, + const int& m2, + const int& N2, + const ModuleBase::Vector3& Efield, + const ModuleBase::Vector3& dR); + void get_psi_r_beta(const UnitCell& ucell, + std::vector>& nlm, + const ModuleBase::Vector3& R1, + const int& T1, + const int& L1, + const int& m1, + const int& N1, + const ModuleBase::Vector3& R2, + const int& T2); void out_rR(const UnitCell& ucell, const Grid_Driver& gd, const int& istep, const int precision = 16); void out_rR_other(const UnitCell& ucell, const int& istep, @@ -78,8 +72,8 @@ class cal_r_overlap_R private: void initialize_orb_table(const UnitCell& ucell, const LCAO_Orbitals& orb); - void construct_orbs_and_orb_r(const UnitCell& ucell,const LCAO_Orbitals& orb); - void construct_orbs_and_nonlocal_and_orb_r(const UnitCell& ucell,const LCAO_Orbitals& orb); + void construct_orbs_and_orb_r(const UnitCell& ucell, const LCAO_Orbitals& orb); + void construct_orbs_and_nonlocal_and_orb_r(const UnitCell& ucell, const LCAO_Orbitals& orb); std::vector iw2ia; std::vector iw2iL; @@ -94,25 +88,15 @@ class cal_r_overlap_R std::vector>> orbs; std::vector> orbs_nonlocal; - std::map< - size_t, - std::map>>>>> + std::map>>>>> center2_orb11; - std::map< - size_t, - std::map>>>>> + std::map>>>>> center2_orb21_r; - std::map< - size_t, - std::map>>>> - center2_orb11_nonlocal; + std::map>>>> center2_orb11_nonlocal; - std::map< - size_t, - std::map>>>> - center2_orb21_r_nonlocal; + std::map>>>> center2_orb21_r_nonlocal; const Parallel_Orbitals* ParaV = nullptr; };