From 24fad25ce41abbd165c8ee4aa54d952e710bd83c Mon Sep 17 00:00:00 2001 From: zaotian Date: Wed, 28 Oct 2026 13:05:47 +0800 Subject: [PATCH 1/7] fix(rt): restore MD boundary shifts and phase corrections --- .ci/slurm/config.ini | 6 + source/source_base/parallel_cell.cpp | 20 ++ source/source_base/parallel_cell.h | 2 + source/source_cell/domain_decomposition.cpp | 57 +-- source/source_cell/domain_decomposition.h | 2 + source/source_cell/mdcell.cpp | 94 +++-- .../source_cell/module_neighlist/local_atom.h | 6 + .../source_esolver/esolver_ks_lcao_tddft.cpp | 10 +- source/source_lcao/module_rt/boundary_fix.cpp | 33 +- source/source_lcao/module_rt/boundary_fix.h | 10 +- source/source_md/md_base.cpp | 28 +- source/source_md/test/test_md_base.cpp | 46 ++- source/source_md/unittests/CMakeLists.txt | 16 + source/source_md/unittests/test_md_base.cpp | 340 ++++++++++++++++++ .../05_rtTDDFT/29_NO_MD_cross_boundary/INPUT | 31 ++ tests/05_rtTDDFT/29_NO_MD_cross_boundary/KPT | 4 + .../05_rtTDDFT/29_NO_MD_cross_boundary/README | 1 + tests/05_rtTDDFT/29_NO_MD_cross_boundary/STRU | 22 ++ .../29_NO_MD_cross_boundary/result.ref | 3 + tests/05_rtTDDFT/CASES_CPU.txt | 1 + tests/05_rtTDDFT/CASES_GPU.txt | 1 + 21 files changed, 647 insertions(+), 86 deletions(-) create mode 100644 source/source_md/unittests/test_md_base.cpp create mode 100644 tests/05_rtTDDFT/29_NO_MD_cross_boundary/INPUT create mode 100644 tests/05_rtTDDFT/29_NO_MD_cross_boundary/KPT create mode 100644 tests/05_rtTDDFT/29_NO_MD_cross_boundary/README create mode 100644 tests/05_rtTDDFT/29_NO_MD_cross_boundary/STRU create mode 100644 tests/05_rtTDDFT/29_NO_MD_cross_boundary/result.ref diff --git a/.ci/slurm/config.ini b/.ci/slurm/config.ini index 5ccd00630ac..78a86416f89 100644 --- a/.ci/slurm/config.ini +++ b/.ci/slurm/config.ini @@ -631,3 +631,9 @@ suite = 01_PW name = 090_PW_out_pchg_wfc_spinor resource = pw_gpu1 runner = autotest_gpu + +[case.101] +suite = 05_rtTDDFT +name = 29_NO_MD_cross_boundary +resource = gpu2 +runner = autotest_gpu diff --git a/source/source_base/parallel_cell.cpp b/source/source_base/parallel_cell.cpp index 0531e288ce6..2564cbbf1b7 100644 --- a/source/source_base/parallel_cell.cpp +++ b/source/source_base/parallel_cell.cpp @@ -61,6 +61,26 @@ int CommunicationDomain::sum(int value) const return value; } +void CommunicationDomain::sum(int* values, int count) const +{ +#ifdef __MPI + if (communicator_ != MPI_COMM_NULL) + { + MPI_Allreduce(MPI_IN_PLACE, values, count, MPI_INT, MPI_SUM, communicator_); + } +#endif +} + +void CommunicationDomain::sum(double* values, int count) const +{ +#ifdef __MPI + if (communicator_ != MPI_COMM_NULL) + { + MPI_Allreduce(MPI_IN_PLACE, values, count, MPI_DOUBLE, MPI_SUM, communicator_); + } +#endif +} + double CommunicationDomain::max(double value) const { #ifdef __MPI diff --git a/source/source_base/parallel_cell.h b/source/source_base/parallel_cell.h index bd02ccaa035..e705a3d70b8 100644 --- a/source/source_base/parallel_cell.h +++ b/source/source_base/parallel_cell.h @@ -21,6 +21,8 @@ class CommunicationDomain int size() const; int max(int value) const; int sum(int value) const; + void sum(int* values, int count) const; + void sum(double* values, int count) const; double max(double value) const; /// Gather one bounded text batch on rank zero; other ranks receive an empty string. /// maximum_count bounds both local payloads and their total in MPI count units. diff --git a/source/source_cell/domain_decomposition.cpp b/source/source_cell/domain_decomposition.cpp index 8565e8c6a10..b0bca2908fb 100644 --- a/source/source_cell/domain_decomposition.cpp +++ b/source/source_cell/domain_decomposition.cpp @@ -992,6 +992,11 @@ DomainDecomposition::PackedAtom DomainDecomposition::pack_atom( packed.mbl[0] = atom.mbl.x; packed.mbl[1] = atom.mbl.y; packed.mbl[2] = atom.mbl.z; + for (int axis = 0; axis < 3; ++axis) + { + packed.displacement[axis] = atom.displacement[axis]; + packed.boundary_shift[axis] = atom.boundary_shift[axis]; + } packed.mass = atom.mass; packed.image_shift[0] = image_shift[0]; packed.image_shift[1] = image_shift[1]; @@ -1006,21 +1011,27 @@ LocalAtom DomainDecomposition::unpack_ghost_atom(const PackedAtom& packed) const { const ModuleBase::Vector3 frac(packed.frac[0], packed.frac[1], packed.frac[2]); const ModuleBase::Vector3 image_frac(packed.frac[0] + packed.image_shift[0], - packed.frac[1] + packed.image_shift[1], - packed.frac[2] + packed.image_shift[2]); + packed.frac[1] + packed.image_shift[1], + packed.frac[2] + packed.image_shift[2]); const ModuleBase::Vector3 cart = image_frac * latvec_; const ModuleBase::Vector3 vel(packed.vel[0], packed.vel[1], packed.vel[2]); const ModuleBase::Vector3 force(packed.force[0], packed.force[1], packed.force[2]); const ModuleBase::Vector3 mbl(packed.mbl[0], packed.mbl[1], packed.mbl[2]); - return LocalAtom(cart, - frac, - vel, - force, - mbl, - packed.mass, - packed.type, - packed.type_index, - packed.owner_rank); + LocalAtom atom(cart, + frac, + vel, + force, + mbl, + packed.mass, + packed.type, + packed.type_index, + packed.owner_rank); + for (int axis = 0; axis < 3; ++axis) + { + atom.displacement[axis] = packed.displacement[axis]; + atom.boundary_shift[axis] = packed.boundary_shift[axis]; + } + return atom; } LocalAtom DomainDecomposition::unpack_owned_atom(const PackedAtom& packed) const @@ -1030,13 +1041,19 @@ LocalAtom DomainDecomposition::unpack_owned_atom(const PackedAtom& packed) const const ModuleBase::Vector3 vel(packed.vel[0], packed.vel[1], packed.vel[2]); const ModuleBase::Vector3 force(packed.force[0], packed.force[1], packed.force[2]); const ModuleBase::Vector3 mbl(packed.mbl[0], packed.mbl[1], packed.mbl[2]); - return LocalAtom(cart, - frac, - vel, - force, - mbl, - packed.mass, - packed.type, - packed.type_index, - packed.owner_rank); + LocalAtom atom(cart, + frac, + vel, + force, + mbl, + packed.mass, + packed.type, + packed.type_index, + packed.owner_rank); + for (int axis = 0; axis < 3; ++axis) + { + atom.displacement[axis] = packed.displacement[axis]; + atom.boundary_shift[axis] = packed.boundary_shift[axis]; + } + return atom; } diff --git a/source/source_cell/domain_decomposition.h b/source/source_cell/domain_decomposition.h index c6a57829f23..053ee512c5a 100644 --- a/source/source_cell/domain_decomposition.h +++ b/source/source_cell/domain_decomposition.h @@ -60,6 +60,8 @@ class DomainDecomposition double frac[3]; double vel[3]; double force[3]; + double displacement[3]; + int boundary_shift[3]; int mbl[3]; double mass; int image_shift[3]; diff --git a/source/source_cell/mdcell.cpp b/source/source_cell/mdcell.cpp index 486f539f9c7..e1abea1e9a4 100644 --- a/source/source_cell/mdcell.cpp +++ b/source/source_cell/mdcell.cpp @@ -185,6 +185,10 @@ void MDCell::sync_backing_unitcell() } sync_backing_unitcell_geometry_(); + for (int it = 0; it < backing_unitcell_->ntype; ++it) + { + backing_unitcell_->atoms[it].boundary_shift.resize(backing_unitcell_->atoms[it].na); + } #ifdef __MPI if (size_ > 1) @@ -195,36 +199,42 @@ void MDCell::sync_backing_unitcell() type_offset[it + 1] = type_offset[it] + backing_unitcell_->atoms[it].na; } - std::vector cart(3 * nat_, 0.0); - std::vector frac(3 * nat_, 0.0); - std::vector vel(3 * nat_, 0.0); - std::vector mbl(3 * nat_, 0); + const int coordinate_count = 3 * nat_; + std::vector cart(coordinate_count, 0.0); + std::vector frac(coordinate_count, 0.0); + std::vector vel(coordinate_count, 0.0); + std::vector displacement(coordinate_count, 0.0); + std::vector boundary_shift(coordinate_count, 0); + std::vector mbl(coordinate_count, 0); std::vector owner(nat_, 0); for (std::size_t i = 0; i < owned_atoms_.size(); ++i) { const LocalAtom& atom = owned_atoms_[i]; const int iat = type_offset[atom.type] + atom.type_index; - cart[3 * iat] = atom.cart.x; - cart[3 * iat + 1] = atom.cart.y; - cart[3 * iat + 2] = atom.cart.z; - frac[3 * iat] = atom.frac.x; - frac[3 * iat + 1] = atom.frac.y; - frac[3 * iat + 2] = atom.frac.z; - vel[3 * iat] = atom.vel.x; - vel[3 * iat + 1] = atom.vel.y; - vel[3 * iat + 2] = atom.vel.z; - mbl[3 * iat] = atom.mbl.x; - mbl[3 * iat + 1] = atom.mbl.y; - mbl[3 * iat + 2] = atom.mbl.z; + for (int axis = 0; axis < 3; ++axis) + { + const int component = 3 * iat + axis; + cart[component] = atom.cart[axis]; + frac[component] = atom.frac[axis]; + vel[component] = atom.vel[axis]; + displacement[component] = atom.displacement[axis]; + boundary_shift[component] = atom.boundary_shift[axis]; + mbl[component] = atom.mbl[axis]; + } owner[iat] = 1; } - MPI_Allreduce(MPI_IN_PLACE, cart.data(), 3 * nat_, MPI_DOUBLE, MPI_SUM, comm_); - MPI_Allreduce(MPI_IN_PLACE, frac.data(), 3 * nat_, MPI_DOUBLE, MPI_SUM, comm_); - MPI_Allreduce(MPI_IN_PLACE, vel.data(), 3 * nat_, MPI_DOUBLE, MPI_SUM, comm_); - MPI_Allreduce(MPI_IN_PLACE, mbl.data(), 3 * nat_, MPI_INT, MPI_SUM, comm_); - MPI_Allreduce(MPI_IN_PLACE, owner.data(), nat_, MPI_INT, MPI_SUM, comm_); + ModuleBase::CommunicationDomain domain; + domain.initialize(comm_); + domain.sum(cart.data(), coordinate_count); + domain.sum(frac.data(), coordinate_count); + domain.sum(vel.data(), coordinate_count); + domain.sum(displacement.data(), coordinate_count); + domain.sum(boundary_shift.data(), coordinate_count); + domain.sum(mbl.data(), coordinate_count); + const int atom_count = nat_; + domain.sum(owner.data(), atom_count); for (int it = 0; it < backing_unitcell_->ntype; ++it) { @@ -235,27 +245,19 @@ void MDCell::sync_backing_unitcell() { throw std::runtime_error("MDCell backing UnitCell atom ownership is invalid."); } - backing_unitcell_->atoms[it].tau[ia].set(cart[3 * iat], cart[3 * iat + 1], cart[3 * iat + 2]); - ModuleBase::Vector3 displacement(frac[3 * iat] - backing_unitcell_->atoms[it].taud[ia].x, - frac[3 * iat + 1] - backing_unitcell_->atoms[it].taud[ia].y, - frac[3 * iat + 2] - backing_unitcell_->atoms[it].taud[ia].z); - for (int k = 0; k < 3; ++k) + for (int axis = 0; axis < 3; ++axis) { - if (displacement[k] > 0.5) - { - displacement[k] -= 1.0; - } - else if (displacement[k] < -0.5) - { - displacement[k] += 1.0; - } + const int component = 3 * iat + axis; + backing_unitcell_->atoms[it].tau[ia][axis] = cart[component]; + backing_unitcell_->atoms[it].taud[ia][axis] = frac[component]; + backing_unitcell_->atoms[it].vel[ia][axis] = vel[component]; + backing_unitcell_->atoms[it].dis[ia][axis] = displacement[component]; + backing_unitcell_->atoms[it].boundary_shift[ia][axis] = boundary_shift[component]; + backing_unitcell_->atoms[it].mbl[ia][axis] = mbl[component]; } - backing_unitcell_->atoms[it].taud[ia].set(frac[3 * iat], frac[3 * iat + 1], frac[3 * iat + 2]); - backing_unitcell_->atoms[it].dis[ia] = displacement; - backing_unitcell_->atoms[it].vel[ia].set(vel[3 * iat], vel[3 * iat + 1], vel[3 * iat + 2]); - backing_unitcell_->atoms[it].mbl[ia].set(mbl[3 * iat], mbl[3 * iat + 1], mbl[3 * iat + 2]); } } + backing_unitcell_->ionic_position_updated = true; return; } #endif @@ -263,24 +265,14 @@ void MDCell::sync_backing_unitcell() for (std::size_t i = 0; i < owned_atoms_.size(); ++i) { const LocalAtom& atom = owned_atoms_[i]; - ModuleBase::Vector3 displacement = atom.frac - backing_unitcell_->atoms[atom.type].taud[atom.type_index]; - for (int k = 0; k < 3; ++k) - { - if (displacement[k] > 0.5) - { - displacement[k] -= 1.0; - } - else if (displacement[k] < -0.5) - { - displacement[k] += 1.0; - } - } backing_unitcell_->atoms[atom.type].tau[atom.type_index] = atom.cart; backing_unitcell_->atoms[atom.type].taud[atom.type_index] = atom.frac; - backing_unitcell_->atoms[atom.type].dis[atom.type_index] = displacement; + backing_unitcell_->atoms[atom.type].dis[atom.type_index] = atom.displacement; + backing_unitcell_->atoms[atom.type].boundary_shift[atom.type_index] = atom.boundary_shift; backing_unitcell_->atoms[atom.type].vel[atom.type_index] = atom.vel; backing_unitcell_->atoms[atom.type].mbl[atom.type_index] = atom.mbl; } + backing_unitcell_->ionic_position_updated = true; } #ifdef __MPI diff --git a/source/source_cell/module_neighlist/local_atom.h b/source/source_cell/module_neighlist/local_atom.h index a311081c2b1..d74a8d8c465 100644 --- a/source/source_cell/module_neighlist/local_atom.h +++ b/source/source_cell/module_neighlist/local_atom.h @@ -20,6 +20,8 @@ struct LocalAtom ModuleBase::Vector3 frac; ModuleBase::Vector3 vel; ModuleBase::Vector3 force; + ModuleBase::Vector3 displacement; ///< Unwrapped fractional displacement in the current MD step. + ModuleBase::Vector3 boundary_shift; ///< Lattice translation applied in the current MD step. ModuleBase::Vector3 mbl; double mass; int type; @@ -31,6 +33,8 @@ struct LocalAtom frac(0.0, 0.0, 0.0), vel(0.0, 0.0, 0.0), force(0.0, 0.0, 0.0), + displacement(0.0, 0.0, 0.0), + boundary_shift(0, 0, 0), mbl(1, 1, 1), mass(1.0), type(0), @@ -52,6 +56,8 @@ struct LocalAtom frac(frac_in), vel(vel_in), force(force_in), + displacement(0.0, 0.0, 0.0), + boundary_shift(0, 0, 0), mbl(mbl_in), mass(mass_in), type(type_in), diff --git a/source/source_esolver/esolver_ks_lcao_tddft.cpp b/source/source_esolver/esolver_ks_lcao_tddft.cpp index ff2b1e06e40..8ac4f81f462 100644 --- a/source/source_esolver/esolver_ks_lcao_tddft.cpp +++ b/source/source_esolver/esolver_ks_lcao_tddft.cpp @@ -202,7 +202,15 @@ void ESolver_KS_LCAO_TDDFT::runner(BaseCell& basecell, const int ist // Update the last-step matrices and wave functions if any atom crosses the boundary // Apply a phase correction to H, S, and psi to keep consistency when atoms cross periodic boundaries const size_t len_hs_ik = use_tensor && use_lapack ? PARAM.globalv.nlocal * PARAM.globalv.nlocal : this->pv.nloc; - module_rt::reset_matrix_boundary(ucell, this->kv, &(this->pv), this->Hk_laststep, this->Sk_laststep, this->psi_laststep, len_hs_ik); + const int kpoint_count = this->kv.get_nks(); + module_rt::reset_matrix_boundary(ucell, + this->kv.kvec_d, + kpoint_count, + &(this->pv), + this->Hk_laststep, + this->Sk_laststep, + this->psi_laststep, + len_hs_ik); for (int estep = 0; estep < estep_max; estep++) { diff --git a/source/source_lcao/module_rt/boundary_fix.cpp b/source/source_lcao/module_rt/boundary_fix.cpp index 7629a7177f4..570017c0387 100644 --- a/source/source_lcao/module_rt/boundary_fix.cpp +++ b/source/source_lcao/module_rt/boundary_fix.cpp @@ -1,12 +1,14 @@ #include "boundary_fix.h" #include "source_base/libm/libm.h" #include "source_base/constants.h" +#include "source_base/parallel_cell.h" #include "source_base/vector3.h" namespace module_rt{ void reset_matrix_boundary(const UnitCell& ucell, - const K_Vectors& kv, + const std::vector>& kvec_d, + const int nks, const Parallel_Orbitals* pv, ct::Tensor& hk_last, ct::Tensor& sk_last, @@ -15,6 +17,20 @@ void reset_matrix_boundary(const UnitCell& ucell, { ModuleBase::TITLE("module_rt", "reset_matrix_boundary"); ModuleBase::timer::start("module_rt", "reset_matrix_boundary"); + // Tensor-LAPACK keeps full H/S matrices on rank zero, while psi stays distributed. + const bool gathered_matrices = len_hs != static_cast(pv->nloc); + Parallel_Orbitals full_distribution; + const Parallel_Orbitals* matrix_distribution = pv; + bool update_matrices = true; + if (gathered_matrices) + { + const int basis_size = pv->get_global_row_size(); + full_distribution.set_serial(basis_size, basis_size); + full_distribution.set_atomic_trace(ucell.get_iat2iwt(), ucell.nat, basis_size); + matrix_distribution = &full_distribution; + const ModuleBase::CommunicationDomain domain = ModuleBase::world_comm_domain(); + update_matrices = domain.rank() == 0; + } const ModuleBase::Vector3 zero = {0, 0, 0}; for(size_t iat = 0; iat < ucell.nat; iat++) { @@ -26,10 +42,10 @@ void reset_matrix_boundary(const UnitCell& ucell, #ifdef _OPENMP #pragma omp parallel for schedule(dynamic) #endif - for(int ik = 0; ik < kv.get_nks(); ik++) + for(int ik = 0; ik < nks; ik++) { const ModuleBase::Vector3 tmp_rshift(rshift.x, rshift.y, rshift.z); - const double arg = -kv.kvec_d[ik] * tmp_rshift * ModuleBase::TWO_PI; + const double arg = -kvec_d[ik] * tmp_rshift * ModuleBase::TWO_PI; //skip unrelevent ik if(arg==0)continue; //calculate correction phase @@ -37,8 +53,13 @@ void reset_matrix_boundary(const UnitCell& ucell, ModuleBase::libm::sincos(arg, &sinp, &cosp); const std::complex phase = std::complex(cosp, sinp); //phase correction for Hamiltionian, overlap matrix and c vec. - module_rt::boundary_shift_mat(phase, hk_last.template data>() + ik * len_hs, pv, iat); - module_rt::boundary_shift_mat(phase, sk_last.template data>() + ik * len_hs, pv, iat); + if (update_matrices) + { + std::complex* h = hk_last.template data>() + ik * len_hs; + std::complex* s = sk_last.template data>() + ik * len_hs; + module_rt::boundary_shift_mat(phase, h, matrix_distribution, iat); + module_rt::boundary_shift_mat(phase, s, matrix_distribution, iat); + } psi_last->fix_k(ik); module_rt::boundary_shift_c(phase, psi_last[0].get_pointer(), pv, iat); } @@ -96,4 +117,4 @@ void boundary_shift_c(const std::complex& phase, } return; } -} //namespace module_rt \ No newline at end of file +} //namespace module_rt diff --git a/source/source_lcao/module_rt/boundary_fix.h b/source/source_lcao/module_rt/boundary_fix.h index 14127cf6a1f..96dbe301cef 100644 --- a/source/source_lcao/module_rt/boundary_fix.h +++ b/source/source_lcao/module_rt/boundary_fix.h @@ -6,7 +6,7 @@ #define BOUNDARY_FIX_H #include "source_cell/unitcell.h" -#include "source_cell/klist.h" +#include #include "source_basis/module_ao/parallel_orbitals.h" #include "source_psi/psi.h" #include "source_base/module_container/ATen/core/tensor.h" @@ -16,7 +16,8 @@ namespace module_rt{ * @brief Add phases to the matrix and coefficient from the previous step to correct the boundary discontinuity. * * @param[in] ucell Unitcell information -* @param[in] kv K-point vectors +* @param[in] kvec_d Direct k-point vectors +* @param[in] nks Number of local k points * @param[in] pv information of parallel * @param[in] hk_last Hamiltonian matrix from last step * @param[in] sk_last Overlap matrix from last step @@ -27,7 +28,8 @@ namespace module_rt{ * @param[out] sk_last the fixed wavefunctions */ void reset_matrix_boundary(const UnitCell& ucell, - const K_Vectors& kv, + const std::vector>& kvec_d, + const int nks, const Parallel_Orbitals* pv, ct::Tensor& hk_last, ct::Tensor& sk_last, @@ -61,4 +63,4 @@ void boundary_shift_c(const std::complex& phase, const Parallel_Orbitals* pv, const size_t iat); }// namespace module_rt -#endif // BOUNDARY_FIX_H \ No newline at end of file +#endif // BOUNDARY_FIX_H diff --git a/source/source_md/md_base.cpp b/source/source_md/md_base.cpp index 6146e25117e..dc583afb136 100644 --- a/source/source_md/md_base.cpp +++ b/source/source_md/md_base.cpp @@ -1,5 +1,6 @@ #include "md_base.h" #include "md_func.h" +#include "source_base/parallel_cell.h" #include "source_cell/unitcell.h" #ifdef __MPI #include "mpi.h" @@ -92,6 +93,8 @@ void MD_base::second_half() void MD_base::update_pos() { std::vector& atoms = mdcell.owned_atoms(); + int invalid_position = 0; + const double epsilon = 1.0e-12; for (std::size_t i = 0; i < atoms.size(); ++i) { LocalAtom& atom = atoms[i]; @@ -108,18 +111,37 @@ void MD_base::update_pos() } } pos = pos * mdcell.GT(); + atom.displacement = pos; atom.frac += pos; + atom.boundary_shift.set(0, 0, 0); for (int k = 0; k < 3; ++k) { - atom.frac[k] -= std::floor(atom.frac[k]); - // A tiny negative value can round to exactly 1.0 after adding one. + if (atom.frac[k] < 0.0) + { + ++atom.boundary_shift[k]; + atom.frac[k] += 1.0; + } + // Keep the legacy second check: a tiny negative can round to 1.0. if (atom.frac[k] >= 1.0) { - atom.frac[k] = 0.0; + --atom.boundary_shift[k]; + atom.frac[k] -= 1.0; + } + if (atom.frac[k] < -epsilon || atom.frac[k] >= 1.0 + epsilon) + { + invalid_position = 1; } } atom.cart = atom.frac * mdcell.latvec(); } + ModuleBase::CommunicationDomain domain = ModuleBase::world_comm_domain(); +#ifdef __MPI + domain.initialize(mdcell.communicator()); +#endif + if (domain.max(invalid_position) != 0) + { + ModuleBase::WARNING_QUIT("MD_base::update_pos", "Movement of atom is larger than the cell length"); + } return; } diff --git a/source/source_md/test/test_md_base.cpp b/source/source_md/test/test_md_base.cpp index 2d0f69cfdb7..35546f4a9c9 100644 --- a/source/source_md/test/test_md_base.cpp +++ b/source/source_md/test/test_md_base.cpp @@ -54,6 +54,8 @@ TEST(MdPositionTest, TinyNegativeDisplacementRemainsInsideHalfOpenCell) EXPECT_LT(moving.frac[axis], 1.0); EXPECT_DOUBLE_EQ(moving.frac[axis], 0.0); EXPECT_DOUBLE_EQ(moving.cart[axis], 0.0); + EXPECT_EQ(moving.boundary_shift[axis], 0); + EXPECT_DOUBLE_EQ(moving.displacement[axis], moving.vel[axis]); } } @@ -65,7 +67,7 @@ TEST(MdPositionTest, FiniteBoundaryCrossingsAndFrozenComponentsKeepTheirPosition parameters.md_tfirst = 0.0; parameters.md_tlast = 0.0; MD_base integrator(parameters, false, true, 0, cell); - cell.owned_atoms()[0].vel.set(-0.5, 1.75, -1.25); + cell.owned_atoms()[0].vel.set(-0.5, 0.75, -1.25); cell.owned_atoms()[1].mbl.set(0, 0, 0); cell.owned_atoms()[1].vel.set(-0.5, 1.75, -1.25); std::ofstream log; @@ -74,9 +76,51 @@ TEST(MdPositionTest, FiniteBoundaryCrossingsAndFrozenComponentsKeepTheirPosition EXPECT_DOUBLE_EQ(moving.frac.x, 0.75); EXPECT_DOUBLE_EQ(moving.frac.y, 0.25); EXPECT_DOUBLE_EQ(moving.frac.z, 0.5); + EXPECT_EQ(moving.boundary_shift.x, 1); + EXPECT_EQ(moving.boundary_shift.y, -1); + EXPECT_EQ(moving.boundary_shift.z, 1); const LocalAtom& frozen = cell.owned_atoms()[1]; EXPECT_DOUBLE_EQ(frozen.frac.x, 0.25); EXPECT_DOUBLE_EQ(frozen.frac.y, 0.5); EXPECT_DOUBLE_EQ(frozen.frac.z, 0.75); } +TEST(MdPositionTest, ExactUpperBoundaryIsWrappedAndTheNextStepClearsItsShift) +{ + MDCell cell = make_position_test_cell(); + MD_para parameters; + parameters.md_dt = ModuleBase::AU_to_FS; + parameters.md_tfirst = 0.0; + parameters.md_tlast = 0.0; + MD_base integrator(parameters, false, true, 0, cell); + LocalAtom& moving = cell.owned_atoms()[0]; + moving.vel.set(0.75, 0.0, 0.0); + std::ofstream log; + integrator.first_half(log); + EXPECT_DOUBLE_EQ(moving.frac.x, 0.0); + EXPECT_DOUBLE_EQ(moving.displacement.x, 0.75); + EXPECT_EQ(moving.boundary_shift.x, -1); + moving.vel.set(0.0, 0.0, 0.0); + integrator.first_half(log); + EXPECT_DOUBLE_EQ(moving.frac.x, 0.0); + EXPECT_DOUBLE_EQ(moving.displacement.x, 0.0); + EXPECT_EQ(moving.boundary_shift.x, 0); +} + +TEST(MdPositionTest, RejectsMovementBeyondTheLegacySingleCellRange) +{ + MDCell cell = make_position_test_cell(); + MD_para parameters; + parameters.md_dt = ModuleBase::AU_to_FS; + parameters.md_tfirst = 0.0; + parameters.md_tlast = 0.0; + MD_base integrator(parameters, false, true, 0, cell); + LocalAtom& moving = cell.owned_atoms()[0]; + std::ofstream log; + for (const double velocity : {-1.5, 2.0}) + { + moving.frac.x = 0.25; + moving.vel.x = velocity; + EXPECT_EXIT(integrator.first_half(log), testing::ExitedWithCode(1), ""); + } +} } // namespace diff --git a/source/source_md/unittests/CMakeLists.txt b/source/source_md/unittests/CMakeLists.txt index 14e6a4fe706..46970aff310 100644 --- a/source/source_md/unittests/CMakeLists.txt +++ b/source/source_md/unittests/CMakeLists.txt @@ -13,7 +13,23 @@ if(ENABLE_MPI) COMMAND ${MPIEXEC_EXECUTABLE} ${MPIEXEC_NUMPROC_FLAG} 4 ${MPIEXEC_PREFLAGS} $ ${MPIEXEC_POSTFLAGS} WORKING_DIRECTORY ${CMAKE_CURRENT_BINARY_DIR}) + addtest( + TARGET MODULE_MD_base_mpi + KEEP_FEATURE_DEFINITIONS __MPI + LIBS parameter psi device container + SOURCES test_md_base.cpp ../md_base.cpp + ../../source_lcao/module_rt/boundary_fix.cpp + ../../source_basis/module_ao/parallel_orbitals.cpp + ../../source_base/parallel_2d.cpp + ../../source_base/parallel_common.cpp + ../../source_cell/magnetism.cpp ${md_test_depend_files}) + add_test( + NAME MODULE_MD_base_mpi_np4 + COMMAND ${MPIEXEC_EXECUTABLE} ${MPIEXEC_NUMPROC_FLAG} 4 ${MPIEXEC_PREFLAGS} + $ ${MPIEXEC_POSTFLAGS} + WORKING_DIRECTORY ${CMAKE_CURRENT_BINARY_DIR}) set_tests_properties( MODULE_MD_dump_mpi MODULE_MD_dump_mpi_np4 + MODULE_MD_base_mpi MODULE_MD_base_mpi_np4 PROPERTIES ENVIRONMENT "OMP_NUM_THREADS=1;MKL_NUM_THREADS=1;I_MPI_PIN=0" TIMEOUT 60) endif() diff --git a/source/source_md/unittests/test_md_base.cpp b/source/source_md/unittests/test_md_base.cpp new file mode 100644 index 00000000000..673b4366174 --- /dev/null +++ b/source/source_md/unittests/test_md_base.cpp @@ -0,0 +1,340 @@ +#include "source_md/md_base.h" +#include "source_base/constants.h" +#include "source_base/parallel_cell.h" +#include "source_base/parallel_global.h" +#include "source_cell/domain_decomposition.h" +#include "source_cell/unitcell.h" +#include "source_cell/update_cell.h" +#include "source_lcao/module_rt/boundary_fix.h" + +#include + +#include +#include +#include +#include + +namespace +{ +class MdBoundaryTest : public testing::Test +{ +protected: + ModuleBase::CommunicationDomain domain; + ModuleBase::Matrix3 lattice; + std::array backing_types; + std::array reference_types; + UnitCell backing; + MDCell cell; + DomainDecomposition decomposition; + MD_para parameters; + + void SetUp() override + { + domain = ModuleBase::world_comm_domain(); + lattice.e11 = 1.0; + lattice.e22 = 1.0; + lattice.e33 = 1.0; + for (Atom* type : {backing_types.data(), reference_types.data()}) + { + type->na = 4; + type->nw = 1; + type->tau.resize(4); + type->taud.resize(4); + type->dis.resize(4); + type->vel.resize(4); + type->mbl.resize(4, ModuleBase::Vector3(1, 1, 1)); + type->boundary_shift.resize(4); + type->taud[0].set(0.99, 0.25, 0.75); + type->taud[1].set(0.01, 0.75, 0.25); + type->taud[2].set(0.0, 0.0, 0.0); + type->taud[3].set(0.5, 0.5, 0.5); + type->mbl[3].set(0, 0, 0); + for (int index = 0; index < 4; ++index) + { + type->tau[index] = type->taud[index] * lattice; + } + } + backing.ntype = 1; + backing.nat = 4; + backing.atoms = backing_types.data(); + backing.set_iat2itia(); + backing.set_iat2iwt(1); + std::vector owned; + const int source_rank = domain.size() - 1; + if (domain.rank() == source_rank) + { + for (int index = 0; index < 4; ++index) + { + LocalAtom atom; + atom.frac = backing_types[0].taud[index]; + atom.cart = backing_types[0].tau[index]; + atom.mbl = backing_types[0].mbl[index]; + atom.type_index = index; + atom.owner_rank = source_rank; + owned.push_back(atom); + } + } + const ModuleBase::Matrix3 inverse = lattice.Inverse(); + const std::vector labels(1, "X"); + const std::vector masses(1, 1.0); + const std::vector counts(1, 4); + cell.initialize_from_owned_atoms(lattice, inverse, 1.0, 1.0, 4, owned, + labels, masses, counts, 0.0, domain); + cell.set_backing_unitcell(backing); + decomposition.init(domain, lattice, 1.0, 0.0, 0.0); + cell.set_neighbor_cutoff(0.05); + parameters.md_dt = ModuleBase::AU_to_FS; + parameters.md_tfirst = -1.0; + parameters.md_tlast = -1.0; + } + + ModuleBase::Vector3 velocity_for(int index, int step) const + { + if (step != 0) + { + return ModuleBase::Vector3(0.0, 0.01, 0.0); + } + if (index == 0) + { + return ModuleBase::Vector3(0.02, 0.6, 0.0); + } + if (index == 1) + { + return ModuleBase::Vector3(-0.02, -0.6, 0.0); + } + return ModuleBase::Vector3(-1.0e-21, -2.0e-21, -3.0e-21); + } + + void expect_local_state(const LocalAtom& atom) + { + const int index = atom.type_index; + const Atom& reference = reference_types[0]; + SCOPED_TRACE(index); + for (int axis = 0; axis < 3; ++axis) + { + EXPECT_DOUBLE_EQ(atom.frac[axis], reference.taud[index][axis]); + EXPECT_DOUBLE_EQ(atom.displacement[axis], reference.dis[index][axis]); + EXPECT_EQ(atom.boundary_shift[axis], reference.boundary_shift[index][axis]); + } + } + + void expect_backing_state() + { + const Atom& actual = backing_types[0]; + const Atom& reference = reference_types[0]; + EXPECT_TRUE(backing.ionic_position_updated); + for (int index = 0; index < 4; ++index) + { + SCOPED_TRACE(index); + for (int axis = 0; axis < 3; ++axis) + { + EXPECT_DOUBLE_EQ(actual.taud[index][axis], reference.taud[index][axis]); + EXPECT_DOUBLE_EQ(actual.tau[index][axis], reference.tau[index][axis]); + EXPECT_DOUBLE_EQ(actual.dis[index][axis], reference.dis[index][axis]); + EXPECT_EQ(actual.boundary_shift[index][axis], reference.boundary_shift[index][axis]); + } + } + } + using Complex = std::complex; + + Complex phase_for(int index, int kpoint) const + { + const double kx = kpoint == 0 ? 0.0 : 0.25; + const double argument = -ModuleBase::TWO_PI * kx * reference_types[0].boundary_shift[index].x; + const Complex exponent(0.0, argument); + return std::exp(exponent); + } + + Complex matrix_value(int row, int column) const + { + const double real = 1.0 + row + column; + const double imaginary = 0.1 * (row - column); + return Complex(real, imaginary); + } + + void expect_phase_correction(bool gathered) + { + const int basis_size = 4; + const int band_count = 2; + Parallel_Orbitals orbitals; + const int block_size = 1; + const int init_result = orbitals.init(basis_size, basis_size, block_size, cell.communicator()); + ASSERT_EQ(init_result, 0); + std::ofstream log; + const int band_result = orbitals.set_nloc_wfc_Eij(band_count, log, log); + ASSERT_EQ(band_result, 0); + orbitals.set_atomic_trace(backing.get_iat2iwt(), backing.nat, basis_size); + std::vector> points(2); + points[1].set(0.25, 0.0, 0.0); + const int rows = gathered ? basis_size : orbitals.get_row_size(); + const int columns = gathered ? basis_size : orbitals.get_col_size(); + const int matrix_size = rows * columns; + const ct::TensorShape shape({2, matrix_size}); + ct::Tensor h(ct::DataType::DT_COMPLEX_DOUBLE, shape); + ct::Tensor s(ct::DataType::DT_COMPLEX_DOUBLE, shape); + const std::vector plane_wave_counts; + psi::Psi coefficients(2, orbitals.ncol_bands, orbitals.nrow, plane_wave_counts, true); + for (int kpoint = 0; kpoint < 2; ++kpoint) + { + for (int column = 0; column < columns; ++column) + { + const int global_column = gathered ? column : orbitals.local2global_col(column); + for (int row = 0; row < rows; ++row) + { + const int global_row = gathered ? row : orbitals.local2global_row(row); + const int element = kpoint * matrix_size + column * rows + row; + h.data()[element] = matrix_value(global_row, global_column); + s.data()[element] = 0.1 * h.data()[element]; + } + } + coefficients.fix_k(kpoint); + for (int element = 0; element < orbitals.nrow * orbitals.ncol_bands; ++element) + { + coefficients.get_pointer()[element] = Complex(1.0, 0.5); + } + } + module_rt::reset_matrix_boundary(backing, points, 2, &orbitals, h, s, &coefficients, matrix_size); + for (int kpoint = 0; kpoint < 2; ++kpoint) + { + for (int column = 0; column < columns; ++column) + { + const int global_column = gathered ? column : orbitals.local2global_col(column); + for (int row = 0; row < rows; ++row) + { + const int global_row = gathered ? row : orbitals.local2global_row(row); + const int element = kpoint * matrix_size + column * rows + row; + Complex phase(1.0, 0.0); + if (!gathered || domain.rank() == 0) + { + phase = std::conj(phase_for(global_row, kpoint)) * phase_for(global_column, kpoint); + } + const Complex expected_h = matrix_value(global_row, global_column) * phase; + const Complex expected_s = 0.1 * expected_h; + const double h_error = std::abs(h.data()[element] - expected_h); + const double s_error = std::abs(s.data()[element] - expected_s); + EXPECT_LT(h_error, 1.0e-12); + EXPECT_LT(s_error, 1.0e-12); + } + } + coefficients.fix_k(kpoint); + for (int column = 0; column < orbitals.ncol_bands; ++column) + { + for (int row = 0; row < orbitals.nrow; ++row) + { + const int global_row = orbitals.local2global_row(row); + const Complex expected = Complex(1.0, 0.5) * std::conj(phase_for(global_row, kpoint)); + const int element = column * orbitals.nrow + row; + const double error = std::abs(coefficients.get_pointer()[element] - expected); + EXPECT_LT(error, 1.0e-12); + } + } + } + } +}; + +TEST_F(MdBoundaryTest, MatchesLegacyWrappingThroughMigrationAndRepeatedSync) +{ + const int rank = domain.rank(); + MD_base integrator(parameters, false, true, rank, cell); + std::ofstream log; + for (int step = 0; step < 2; ++step) + { + SCOPED_TRACE(step); + for (int index = 0; index < 4; ++index) + { + ModuleBase::Vector3 displacement = velocity_for(index, step); + for (int axis = 0; axis < 3; ++axis) + { + if (!reference_types[0].mbl[index][axis]) + { + displacement[axis] = 0.0; + } + } + reference_types[0].taud[index] += displacement; + reference_types[0].dis[index] = displacement; + } + unitcell::periodic_boundary_adjustment(reference_types.data(), lattice, 1); + for (LocalAtom& atom : cell.owned_atoms()) + { + const int index = atom.type_index; + atom.vel = velocity_for(index, step); + } + integrator.first_half(log); + for (const LocalAtom& atom : cell.owned_atoms()) + { + expect_local_state(atom); + } + // Before migration, rank zero and the other nonowners can be empty. + backing.ionic_position_updated = false; + cell.sync_backing_unitcell(); + expect_backing_state(); + decomposition.migrate_owned_atoms(cell); + int migrated = 0; + for (const LocalAtom& atom : cell.owned_atoms()) + { + expect_local_state(atom); + const int owner = decomposition.owner_rank_from_frac(atom.frac); + EXPECT_EQ(atom.owner_rank, owner); + if (owner != domain.size() - 1) + { + ++migrated; + } + } + const int global_migrated = domain.sum(migrated); + if (domain.size() > 1) + { + EXPECT_GT(global_migrated, 0); + } + for (const LocalAtom& ghost : cell.ghost_atoms()) + { + expect_local_state(ghost); + } + // Restart output may synchronize the same MD step a second time. + for (int repeat = 0; repeat < 2; ++repeat) + { + backing.ionic_position_updated = false; + cell.sync_backing_unitcell(); + expect_backing_state(); + } + expect_phase_correction(false); + expect_phase_correction(true); + const int crossing = step == 0 ? 1 : 0; + EXPECT_EQ(backing_types[0].boundary_shift[0].x, -crossing); + EXPECT_EQ(backing_types[0].boundary_shift[1].x, crossing); + EXPECT_EQ(backing_types[0].boundary_shift[2].norm2(), 0); + EXPECT_EQ(backing_types[0].boundary_shift[3].norm2(), 0); + } +} + +TEST_F(MdBoundaryTest, InitialSyncClearsStaleBackingState) +{ + for (int index = 0; index < 4; ++index) + { + backing_types[0].boundary_shift[index].set(7, 8, 9); + backing_types[0].dis[index].set(7.0, 8.0, 9.0); + } + backing.ionic_position_updated = false; + cell.sync_backing_unitcell(); + expect_backing_state(); +} +} // namespace + +int main(int argc, char** argv) +{ + int process_count = 1; + int thread_count = 1; + int rank = 0; + Parallel_Global::read_pal_param(argc, argv, process_count, thread_count, rank); + POOL_WORLD = MPI_COMM_NULL; + KP_WORLD = MPI_COMM_NULL; + INT_BGROUP = MPI_COMM_NULL; + BP_WORLD = MPI_COMM_NULL; + GRID_WORLD = MPI_COMM_NULL; + DIAG_WORLD = MPI_COMM_NULL; + ::testing::InitGoogleTest(&argc, argv); + const int local_result = RUN_ALL_TESTS(); + const ModuleBase::CommunicationDomain domain = ModuleBase::world_comm_domain(); + const int result = domain.max(local_result); + Parallel_Global::finalize_mpi(); + return result; +} diff --git a/tests/05_rtTDDFT/29_NO_MD_cross_boundary/INPUT b/tests/05_rtTDDFT/29_NO_MD_cross_boundary/INPUT new file mode 100644 index 00000000000..9e28ec729ef --- /dev/null +++ b/tests/05_rtTDDFT/29_NO_MD_cross_boundary/INPUT @@ -0,0 +1,31 @@ +INPUT_PARAMETERS +suffix autotest +pseudo_dir ../../PP_ORB +orbital_dir ../../PP_ORB +calculation md +esolver_type tddft +basis_type lcao +gamma_only 0 +symmetry -1 +nbands 2 +ecutwfc 20 + +scf_thr 1e-8 +scf_nmax 100 +mixing_type broyden +mixing_beta 0.7 +mixing_gg0 0.0 + +md_type nve +md_nstep 3 +estep_per_md 2 +md_dt 0.05 +td_dt 0.025 +init_vel 1 +# Keep the supplied finite crossing velocities. +md_tfirst -1 +md_dumpfreq 1 +md_restartfreq 1 + +td_vext 0 +td_stype 0 diff --git a/tests/05_rtTDDFT/29_NO_MD_cross_boundary/KPT b/tests/05_rtTDDFT/29_NO_MD_cross_boundary/KPT new file mode 100644 index 00000000000..5184d98cf62 --- /dev/null +++ b/tests/05_rtTDDFT/29_NO_MD_cross_boundary/KPT @@ -0,0 +1,4 @@ +K_POINTS +1 +Direct +0.25 0.0 0.0 1.0 diff --git a/tests/05_rtTDDFT/29_NO_MD_cross_boundary/README b/tests/05_rtTDDFT/29_NO_MD_cross_boundary/README new file mode 100644 index 00000000000..1cc5fb6e741 --- /dev/null +++ b/tests/05_rtTDDFT/29_NO_MD_cross_boundary/README @@ -0,0 +1 @@ +Minimal LCAO rt-TDDFT regression for finite periodic crossings in both directions at a non-Gamma k point. diff --git a/tests/05_rtTDDFT/29_NO_MD_cross_boundary/STRU b/tests/05_rtTDDFT/29_NO_MD_cross_boundary/STRU new file mode 100644 index 00000000000..bc863b1562e --- /dev/null +++ b/tests/05_rtTDDFT/29_NO_MD_cross_boundary/STRU @@ -0,0 +1,22 @@ +ATOMIC_SPECIES +H 1.00794 H.LDA.UPF + +NUMERICAL_ORBITAL +H_lda_8.0au_100Ry_2s1p.orb + +LATTICE_CONSTANT +15 + +LATTICE_VECTORS +1 0 0 +0 1 0 +0 0 1 + +ATOMIC_POSITIONS +Direct + +H +0 +2 +0.9999 0.5 0.45 m 1 0 0 v 0.001 0 0 +0.0001 0.5 0.55 m 1 0 0 v -0.001 0 0 diff --git a/tests/05_rtTDDFT/29_NO_MD_cross_boundary/result.ref b/tests/05_rtTDDFT/29_NO_MD_cross_boundary/result.ref new file mode 100644 index 00000000000..5b36622d81b --- /dev/null +++ b/tests/05_rtTDDFT/29_NO_MD_cross_boundary/result.ref @@ -0,0 +1,3 @@ +etotref -30.89775473222166 +etotperatomref -15.44887736611083 +totaltimeref 2.5 diff --git a/tests/05_rtTDDFT/CASES_CPU.txt b/tests/05_rtTDDFT/CASES_CPU.txt index 8cbca6fadab..0cb759facdb 100644 --- a/tests/05_rtTDDFT/CASES_CPU.txt +++ b/tests/05_rtTDDFT/CASES_CPU.txt @@ -26,3 +26,4 @@ 26_PW_single_TDDFT 27_PW_gmres_recycle_TDDFT 28_PW_gmres_subspace_TDDFT +29_NO_MD_cross_boundary diff --git a/tests/05_rtTDDFT/CASES_GPU.txt b/tests/05_rtTDDFT/CASES_GPU.txt index c2cf48ba997..82bd2ae4719 100644 --- a/tests/05_rtTDDFT/CASES_GPU.txt +++ b/tests/05_rtTDDFT/CASES_GPU.txt @@ -26,3 +26,4 @@ 26_PW_single_TDDFT 27_PW_gmres_recycle_TDDFT 28_PW_gmres_subspace_TDDFT +29_NO_MD_cross_boundary From d3731877b07955c6eca7b7b06e928c488cd439f9 Mon Sep 17 00:00:00 2001 From: zaotian Date: Wed, 28 Oct 2026 13:19:35 +0800 Subject: [PATCH 2/7] test(rt): merge finite boundary regression into case 06 --- .ci/slurm/config.ini | 6 ---- tests/05_rtTDDFT/06_NO_dir_TDDFT/INPUT | 11 +++++-- tests/05_rtTDDFT/06_NO_dir_TDDFT/KPT | 2 +- tests/05_rtTDDFT/06_NO_dir_TDDFT/README | 1 + tests/05_rtTDDFT/06_NO_dir_TDDFT/STRU | 18 +++++------ tests/05_rtTDDFT/06_NO_dir_TDDFT/result.ref | 10 +++--- .../05_rtTDDFT/29_NO_MD_cross_boundary/INPUT | 31 ------------------- tests/05_rtTDDFT/29_NO_MD_cross_boundary/KPT | 4 --- .../05_rtTDDFT/29_NO_MD_cross_boundary/README | 1 - tests/05_rtTDDFT/29_NO_MD_cross_boundary/STRU | 22 ------------- .../29_NO_MD_cross_boundary/result.ref | 3 -- tests/05_rtTDDFT/CASES_CPU.txt | 1 - tests/05_rtTDDFT/CASES_GPU.txt | 1 - 13 files changed, 24 insertions(+), 87 deletions(-) create mode 100644 tests/05_rtTDDFT/06_NO_dir_TDDFT/README delete mode 100644 tests/05_rtTDDFT/29_NO_MD_cross_boundary/INPUT delete mode 100644 tests/05_rtTDDFT/29_NO_MD_cross_boundary/KPT delete mode 100644 tests/05_rtTDDFT/29_NO_MD_cross_boundary/README delete mode 100644 tests/05_rtTDDFT/29_NO_MD_cross_boundary/STRU delete mode 100644 tests/05_rtTDDFT/29_NO_MD_cross_boundary/result.ref diff --git a/.ci/slurm/config.ini b/.ci/slurm/config.ini index 78a86416f89..5ccd00630ac 100644 --- a/.ci/slurm/config.ini +++ b/.ci/slurm/config.ini @@ -631,9 +631,3 @@ suite = 01_PW name = 090_PW_out_pchg_wfc_spinor resource = pw_gpu1 runner = autotest_gpu - -[case.101] -suite = 05_rtTDDFT -name = 29_NO_MD_cross_boundary -resource = gpu2 -runner = autotest_gpu diff --git a/tests/05_rtTDDFT/06_NO_dir_TDDFT/INPUT b/tests/05_rtTDDFT/06_NO_dir_TDDFT/INPUT index be85d721db0..01e6af7a49a 100644 --- a/tests/05_rtTDDFT/06_NO_dir_TDDFT/INPUT +++ b/tests/05_rtTDDFT/06_NO_dir_TDDFT/INPUT @@ -8,16 +8,17 @@ calculation md esolver_type tddft basis_type lcao gamma_only 0 +symmetry -1 nbands 5 #Parameter (Electronic Structure) ecutwfc 20 -scf_nmax 50 +scf_nmax 100 mixing_type broyden mixing_beta 0.7 mixing_gg0 0.0 -scf_thr 1.0e-6 +scf_thr 1.0e-8 cal_stress 1 stress_thr 1e-6 @@ -27,8 +28,12 @@ force_thr_ev 1.0e-3 #Parameter (MD) md_type nve md_nstep 3 -estep_per_md 1 +estep_per_md 2 md_dt 0.05 +td_dt 0.025 +md_tfirst -1 +md_dumpfreq 1 +md_restartfreq 1 init_vel 1 #Parameter (External Field) diff --git a/tests/05_rtTDDFT/06_NO_dir_TDDFT/KPT b/tests/05_rtTDDFT/06_NO_dir_TDDFT/KPT index e769af76382..03803d6f1dd 100644 --- a/tests/05_rtTDDFT/06_NO_dir_TDDFT/KPT +++ b/tests/05_rtTDDFT/06_NO_dir_TDDFT/KPT @@ -1,4 +1,4 @@ K_POINTS 0 Gamma -2 1 1 0 0 0 +4 1 1 0 0 0 diff --git a/tests/05_rtTDDFT/06_NO_dir_TDDFT/README b/tests/05_rtTDDFT/06_NO_dir_TDDFT/README new file mode 100644 index 00000000000..0f22c0395dd --- /dev/null +++ b/tests/05_rtTDDFT/06_NO_dir_TDDFT/README @@ -0,0 +1 @@ +LCAO rt-TDDFT with a y-directed Gaussian field and finite periodic crossings in both directions. diff --git a/tests/05_rtTDDFT/06_NO_dir_TDDFT/STRU b/tests/05_rtTDDFT/06_NO_dir_TDDFT/STRU index af322391215..eb843c14106 100644 --- a/tests/05_rtTDDFT/06_NO_dir_TDDFT/STRU +++ b/tests/05_rtTDDFT/06_NO_dir_TDDFT/STRU @@ -8,15 +8,15 @@ LATTICE_CONSTANT 15 LATTICE_VECTORS -1 0 0 #latvec1 -0 1 0 #latvec2 -0 0 1 #latvec3 +1 0 0 +0 1 0 +0 0 1 ATOMIC_POSITIONS -Cartesian +Direct -H #label -0 #magnetism -2 #number of atoms -0.502095122795 0.474071446432 0.509122850526 m 1 1 1 v -1.58317459417e-05 0.000196314331857 0.000681467807833 -0.497904877205 0.525928553568 0.592277149473 m 1 1 1 v 1.58317459417e-05 -0.000196314331857 -0.000681467807833 +H +0 +2 +0.9999 0.5 0.45 m 1 1 1 v 0.001 0 0 +0.0001 0.5 0.55 m 1 1 1 v -0.001 0 0 diff --git a/tests/05_rtTDDFT/06_NO_dir_TDDFT/result.ref b/tests/05_rtTDDFT/06_NO_dir_TDDFT/result.ref index 7070f16d32b..4301ad54668 100644 --- a/tests/05_rtTDDFT/06_NO_dir_TDDFT/result.ref +++ b/tests/05_rtTDDFT/06_NO_dir_TDDFT/result.ref @@ -1,5 +1,5 @@ -etotref -30.44916398293971 -etotperatomref -15.2245819915 -totalforceref 0.593208 -totalstressref 4.354288 -totaltimeref 1.55 +etotref -30.89742338434547 +etotperatomref -15.4487116922 +totalforceref 1.490274 +totalstressref 1.942380 +totaltimeref 0.84 diff --git a/tests/05_rtTDDFT/29_NO_MD_cross_boundary/INPUT b/tests/05_rtTDDFT/29_NO_MD_cross_boundary/INPUT deleted file mode 100644 index 9e28ec729ef..00000000000 --- a/tests/05_rtTDDFT/29_NO_MD_cross_boundary/INPUT +++ /dev/null @@ -1,31 +0,0 @@ -INPUT_PARAMETERS -suffix autotest -pseudo_dir ../../PP_ORB -orbital_dir ../../PP_ORB -calculation md -esolver_type tddft -basis_type lcao -gamma_only 0 -symmetry -1 -nbands 2 -ecutwfc 20 - -scf_thr 1e-8 -scf_nmax 100 -mixing_type broyden -mixing_beta 0.7 -mixing_gg0 0.0 - -md_type nve -md_nstep 3 -estep_per_md 2 -md_dt 0.05 -td_dt 0.025 -init_vel 1 -# Keep the supplied finite crossing velocities. -md_tfirst -1 -md_dumpfreq 1 -md_restartfreq 1 - -td_vext 0 -td_stype 0 diff --git a/tests/05_rtTDDFT/29_NO_MD_cross_boundary/KPT b/tests/05_rtTDDFT/29_NO_MD_cross_boundary/KPT deleted file mode 100644 index 5184d98cf62..00000000000 --- a/tests/05_rtTDDFT/29_NO_MD_cross_boundary/KPT +++ /dev/null @@ -1,4 +0,0 @@ -K_POINTS -1 -Direct -0.25 0.0 0.0 1.0 diff --git a/tests/05_rtTDDFT/29_NO_MD_cross_boundary/README b/tests/05_rtTDDFT/29_NO_MD_cross_boundary/README deleted file mode 100644 index 1cc5fb6e741..00000000000 --- a/tests/05_rtTDDFT/29_NO_MD_cross_boundary/README +++ /dev/null @@ -1 +0,0 @@ -Minimal LCAO rt-TDDFT regression for finite periodic crossings in both directions at a non-Gamma k point. diff --git a/tests/05_rtTDDFT/29_NO_MD_cross_boundary/STRU b/tests/05_rtTDDFT/29_NO_MD_cross_boundary/STRU deleted file mode 100644 index bc863b1562e..00000000000 --- a/tests/05_rtTDDFT/29_NO_MD_cross_boundary/STRU +++ /dev/null @@ -1,22 +0,0 @@ -ATOMIC_SPECIES -H 1.00794 H.LDA.UPF - -NUMERICAL_ORBITAL -H_lda_8.0au_100Ry_2s1p.orb - -LATTICE_CONSTANT -15 - -LATTICE_VECTORS -1 0 0 -0 1 0 -0 0 1 - -ATOMIC_POSITIONS -Direct - -H -0 -2 -0.9999 0.5 0.45 m 1 0 0 v 0.001 0 0 -0.0001 0.5 0.55 m 1 0 0 v -0.001 0 0 diff --git a/tests/05_rtTDDFT/29_NO_MD_cross_boundary/result.ref b/tests/05_rtTDDFT/29_NO_MD_cross_boundary/result.ref deleted file mode 100644 index 5b36622d81b..00000000000 --- a/tests/05_rtTDDFT/29_NO_MD_cross_boundary/result.ref +++ /dev/null @@ -1,3 +0,0 @@ -etotref -30.89775473222166 -etotperatomref -15.44887736611083 -totaltimeref 2.5 diff --git a/tests/05_rtTDDFT/CASES_CPU.txt b/tests/05_rtTDDFT/CASES_CPU.txt index 0cb759facdb..8cbca6fadab 100644 --- a/tests/05_rtTDDFT/CASES_CPU.txt +++ b/tests/05_rtTDDFT/CASES_CPU.txt @@ -26,4 +26,3 @@ 26_PW_single_TDDFT 27_PW_gmres_recycle_TDDFT 28_PW_gmres_subspace_TDDFT -29_NO_MD_cross_boundary diff --git a/tests/05_rtTDDFT/CASES_GPU.txt b/tests/05_rtTDDFT/CASES_GPU.txt index 82bd2ae4719..c2cf48ba997 100644 --- a/tests/05_rtTDDFT/CASES_GPU.txt +++ b/tests/05_rtTDDFT/CASES_GPU.txt @@ -26,4 +26,3 @@ 26_PW_single_TDDFT 27_PW_gmres_recycle_TDDFT 28_PW_gmres_subspace_TDDFT -29_NO_MD_cross_boundary From 7b3a835bc70a5f1c7db9ac1fea0abf149cda6f18 Mon Sep 17 00:00:00 2001 From: zaotian Date: Wed, 28 Oct 2026 15:34:41 +0800 Subject: [PATCH 3/7] fix(md): preserve initial charge-extrapolation state Mark synchronized positions updated only after an MD drift. Initial and restart-output synchronization preserve the existing flag, so the initial force evaluation does not advance charge-extrapolation history. Add a regression for the initial and subsequent force evaluations. --- source/source_cell/mdcell.cpp | 12 +++- source/source_cell/mdcell.h | 3 +- source/source_md/md_func.cpp | 4 +- source/source_md/run_md.cpp | 3 +- source/source_md/unittests/test_md_base.cpp | 75 ++++++++++++++++++--- 5 files changed, 80 insertions(+), 17 deletions(-) diff --git a/source/source_cell/mdcell.cpp b/source/source_cell/mdcell.cpp index e1abea1e9a4..4e1870f9b77 100644 --- a/source/source_cell/mdcell.cpp +++ b/source/source_cell/mdcell.cpp @@ -177,7 +177,7 @@ void MDCell::refresh_cart_from_frac() neighbor_layout_valid_ = false; } -void MDCell::sync_backing_unitcell() +void MDCell::sync_backing_unitcell(const bool positions_updated) { if (backing_unitcell_ == nullptr) { @@ -257,7 +257,10 @@ void MDCell::sync_backing_unitcell() } } } - backing_unitcell_->ionic_position_updated = true; + if (positions_updated) + { + backing_unitcell_->ionic_position_updated = true; + } return; } #endif @@ -272,7 +275,10 @@ void MDCell::sync_backing_unitcell() backing_unitcell_->atoms[atom.type].vel[atom.type_index] = atom.vel; backing_unitcell_->atoms[atom.type].mbl[atom.type_index] = atom.mbl; } - backing_unitcell_->ionic_position_updated = true; + if (positions_updated) + { + backing_unitcell_->ionic_position_updated = true; + } } #ifdef __MPI diff --git a/source/source_cell/mdcell.h b/source/source_cell/mdcell.h index 04fa383e26d..4c3edcf7883 100644 --- a/source/source_cell/mdcell.h +++ b/source/source_cell/mdcell.h @@ -57,7 +57,8 @@ class MDCell : public BaseCell const NeighborSearch& neighbor_search() const; void set_lattice_vectors(const ModuleBase::Matrix3& latvec); void refresh_cart_from_frac(); - void sync_backing_unitcell(); + /// Copy state; mark a position update only when requested, otherwise preserve the flag. + void sync_backing_unitcell(const bool positions_updated); std::vector& owned_atoms() { return owned_atoms_; } const std::vector& owned_atoms() const { return owned_atoms_; } diff --git a/source/source_md/md_func.cpp b/source/source_md/md_func.cpp index 3536dd4a0b7..55927c345bf 100644 --- a/source/source_md/md_func.cpp +++ b/source/source_md/md_func.cpp @@ -343,7 +343,9 @@ void force_virial(ModuleESolver::ESolver* p_esolver, { backing_velocities[static_cast(it)] = ucell.atoms[it].vel; } - mdcell.sync_backing_unitcell(); + // Step zero evaluates the initial coordinates before the first MD drift. + const bool positions_updated = istep > 0; + mdcell.sync_backing_unitcell(positions_updated); for (int it = 0; it < ucell.ntype; ++it) { std::copy(backing_velocities[static_cast(it)].begin(), diff --git a/source/source_md/run_md.cpp b/source/source_md/run_md.cpp index 114f9cbb3a5..d71ac5679a1 100644 --- a/source/source_md/run_md.cpp +++ b/source/source_md/run_md.cpp @@ -285,7 +285,8 @@ void md_line(MDCell& mdcell, { if (mdcell.has_backing_unitcell()) { - mdcell.sync_backing_unitcell(); + const bool positions_updated = false; + mdcell.sync_backing_unitcell(positions_updated); } std::stringstream file; file << PARAM.globalv.global_stru_dir << "STRU_MD_" << mdrun->step_ + mdrun->step_rst_; diff --git a/source/source_md/unittests/test_md_base.cpp b/source/source_md/unittests/test_md_base.cpp index 673b4366174..2c57df6f11a 100644 --- a/source/source_md/unittests/test_md_base.cpp +++ b/source/source_md/unittests/test_md_base.cpp @@ -1,4 +1,5 @@ #include "source_md/md_base.h" +#include "source_md/md_func.h" #include "source_base/constants.h" #include "source_base/parallel_cell.h" #include "source_base/parallel_global.h" @@ -16,6 +17,23 @@ namespace { +class PositionStateSolver : public ModuleESolver::ESolver +{ +public: + std::vector observed_position_updates; + + void before_all_runners(BaseCell&, const Input_para&) override {} + void after_all_runners(BaseCell&) override {} + void runner(BaseCell& base, const int) override + { + const UnitCell& cell = static_cast(base); + observed_position_updates.push_back(cell.ionic_position_updated); + } + double cal_energy() override { return 0.0; } + void cal_force(BaseCell&, ModuleBase::matrix& force) override { force.zero_out(); } + void cal_stress(BaseCell&, ModuleBase::matrix& stress) override { stress.zero_out(); } +}; + class MdBoundaryTest : public testing::Test { protected: @@ -118,11 +136,11 @@ class MdBoundaryTest : public testing::Test } } - void expect_backing_state() + void expect_backing_state(const bool positions_updated) { const Atom& actual = backing_types[0]; const Atom& reference = reference_types[0]; - EXPECT_TRUE(backing.ionic_position_updated); + EXPECT_EQ(backing.ionic_position_updated, positions_updated); for (int index = 0; index < 4; ++index) { SCOPED_TRACE(index); @@ -266,8 +284,8 @@ TEST_F(MdBoundaryTest, MatchesLegacyWrappingThroughMigrationAndRepeatedSync) } // Before migration, rank zero and the other nonowners can be empty. backing.ionic_position_updated = false; - cell.sync_backing_unitcell(); - expect_backing_state(); + cell.sync_backing_unitcell(true); + expect_backing_state(true); decomposition.migrate_owned_atoms(cell); int migrated = 0; for (const LocalAtom& atom : cell.owned_atoms()) @@ -292,9 +310,13 @@ TEST_F(MdBoundaryTest, MatchesLegacyWrappingThroughMigrationAndRepeatedSync) // Restart output may synchronize the same MD step a second time. for (int repeat = 0; repeat < 2; ++repeat) { - backing.ionic_position_updated = false; - cell.sync_backing_unitcell(); - expect_backing_state(); + const bool positions_updated = repeat == 0; + if (positions_updated) + { + backing.ionic_position_updated = false; + } + cell.sync_backing_unitcell(positions_updated); + expect_backing_state(true); } expect_phase_correction(false); expect_phase_correction(true); @@ -306,16 +328,47 @@ TEST_F(MdBoundaryTest, MatchesLegacyWrappingThroughMigrationAndRepeatedSync) } } -TEST_F(MdBoundaryTest, InitialSyncClearsStaleBackingState) +TEST_F(MdBoundaryTest, InitialForceEvaluationDoesNotAdvancePositionUpdateState) +{ + const int rank = domain.rank(); + MD_base integrator(parameters, false, true, rank, cell); + PositionStateSolver solver; + const std::string readin_directory; + backing.ionic_position_updated = false; + integrator.setup(&solver, readin_directory, decomposition); + ASSERT_EQ(solver.observed_position_updates.size(), 1u); + EXPECT_FALSE(solver.observed_position_updates[0]); + // The legacy setup marks the state only after the first force evaluation. + EXPECT_TRUE(backing.ionic_position_updated); + + backing.ionic_position_updated = false; + for (LocalAtom& atom : cell.owned_atoms()) + { + const int index = atom.type_index; + atom.vel = velocity_for(index, 0); + } + std::ofstream log; + integrator.first_half(log); + MD_func::force_virial(&solver, 1, cell, decomposition, + integrator.potential, false, integrator.virial, false); + ASSERT_EQ(solver.observed_position_updates.size(), 2u); + EXPECT_TRUE(solver.observed_position_updates[1]); +} + +TEST_F(MdBoundaryTest, InitialSyncPreservesPositionUpdateState) { for (int index = 0; index < 4; ++index) { backing_types[0].boundary_shift[index].set(7, 8, 9); backing_types[0].dis[index].set(7.0, 8.0, 9.0); } - backing.ionic_position_updated = false; - cell.sync_backing_unitcell(); - expect_backing_state(); + const bool positions_updated = false; + for (const bool initial_state : {false, true}) + { + backing.ionic_position_updated = initial_state; + cell.sync_backing_unitcell(positions_updated); + expect_backing_state(initial_state); + } } } // namespace From 37dea1ae4ba682d476e15e67d2df54715770687c Mon Sep 17 00:00:00 2001 From: zaotian Date: Thu, 29 Oct 2026 12:59:47 +0800 Subject: [PATCH 4/7] fix(rt): handle threaded and hybrid-gauge boundary corrections --- .../source_esolver/esolver_ks_lcao_tddft.cpp | 7 + source/source_lcao/module_rt/boundary_fix.cpp | 15 +- source/source_lcao/module_rt/boundary_fix.h | 2 + source/source_md/unittests/CMakeLists.txt | 17 ++ source/source_md/unittests/test_md_base.cpp | 182 +++++++++++++++++- 5 files changed, 211 insertions(+), 12 deletions(-) diff --git a/source/source_esolver/esolver_ks_lcao_tddft.cpp b/source/source_esolver/esolver_ks_lcao_tddft.cpp index 8ac4f81f462..ef4254f3df1 100644 --- a/source/source_esolver/esolver_ks_lcao_tddft.cpp +++ b/source/source_esolver/esolver_ks_lcao_tddft.cpp @@ -147,6 +147,12 @@ void ESolver_KS_LCAO_TDDFT::runner(BaseCell& basecell, const int ist //---------------------------------------------------------------- // 1) before_scf (preparation before the electronic-step loop) //---------------------------------------------------------------- + // Preserve the field endpoint used by the saved H/S/psi before advancing the field. + ModuleBase::Vector3 boundary_A; + if (this->inp_->td_stype == 2) + { + boundary_A = td_p->endpoint_A(); + } // Prepare the endpoint field before before_scf initializes the potential and constructs TD operators. const int first_estep = this->totstep + 1; if (this->inp_->td_stype == 0 || first_estep == 0) @@ -206,6 +212,7 @@ void ESolver_KS_LCAO_TDDFT::runner(BaseCell& basecell, const int ist module_rt::reset_matrix_boundary(ucell, this->kv.kvec_d, kpoint_count, + boundary_A, &(this->pv), this->Hk_laststep, this->Sk_laststep, diff --git a/source/source_lcao/module_rt/boundary_fix.cpp b/source/source_lcao/module_rt/boundary_fix.cpp index 570017c0387..641c3ff9964 100644 --- a/source/source_lcao/module_rt/boundary_fix.cpp +++ b/source/source_lcao/module_rt/boundary_fix.cpp @@ -9,6 +9,7 @@ namespace module_rt{ void reset_matrix_boundary(const UnitCell& ucell, const std::vector>& kvec_d, const int nks, + const ModuleBase::Vector3& cart_At_last, const Parallel_Orbitals* pv, ct::Tensor& hk_last, ct::Tensor& sk_last, @@ -45,7 +46,10 @@ void reset_matrix_boundary(const UnitCell& ucell, for(int ik = 0; ik < nks; ik++) { const ModuleBase::Vector3 tmp_rshift(rshift.x, rshift.y, rshift.z); - const double arg = -kvec_d[ik] * tmp_rshift * ModuleBase::TWO_PI; + // Hybrid folding uses exp(i * (2*pi*k.R + A.R_cart)). + // Transform the saved state with its own field endpoint. + const ModuleBase::Vector3 cart_shift = tmp_rshift * ucell.latvec * ucell.lat0; + const double arg = -kvec_d[ik] * tmp_rshift * ModuleBase::TWO_PI - cart_At_last * cart_shift; //skip unrelevent ik if(arg==0)continue; //calculate correction phase @@ -60,8 +64,13 @@ void reset_matrix_boundary(const UnitCell& ucell, module_rt::boundary_shift_mat(phase, h, matrix_distribution, iat); module_rt::boundary_shift_mat(phase, s, matrix_distribution, iat); } - psi_last->fix_k(ik); - module_rt::boundary_shift_c(phase, psi_last[0].get_pointer(), pv, iat); + if (pv->nrow > 0 && pv->ncol_bands > 0) + { + // Saved LCAO coefficients use k-first storage. Index it without + // changing the shared Psi cursor in this OpenMP loop. + std::complex* psi_k = &(*psi_last)(ik, 0, 0); + module_rt::boundary_shift_c(phase, psi_k, pv, iat); + } } } } diff --git a/source/source_lcao/module_rt/boundary_fix.h b/source/source_lcao/module_rt/boundary_fix.h index 96dbe301cef..0b61b5d0cfb 100644 --- a/source/source_lcao/module_rt/boundary_fix.h +++ b/source/source_lcao/module_rt/boundary_fix.h @@ -18,6 +18,7 @@ namespace module_rt{ * @param[in] ucell Unitcell information * @param[in] kvec_d Direct k-point vectors * @param[in] nks Number of local k points +* @param[in] cart_At_last Saved hybrid-gauge Cartesian A in Hartree atomic units; zero in other gauges * @param[in] pv information of parallel * @param[in] hk_last Hamiltonian matrix from last step * @param[in] sk_last Overlap matrix from last step @@ -30,6 +31,7 @@ namespace module_rt{ void reset_matrix_boundary(const UnitCell& ucell, const std::vector>& kvec_d, const int nks, + const ModuleBase::Vector3& cart_At_last, const Parallel_Orbitals* pv, ct::Tensor& hk_last, ct::Tensor& sk_last, diff --git a/source/source_md/unittests/CMakeLists.txt b/source/source_md/unittests/CMakeLists.txt index 46970aff310..4b84963538f 100644 --- a/source/source_md/unittests/CMakeLists.txt +++ b/source/source_md/unittests/CMakeLists.txt @@ -19,6 +19,10 @@ if(ENABLE_MPI) LIBS parameter psi device container SOURCES test_md_base.cpp ../md_base.cpp ../../source_lcao/module_rt/boundary_fix.cpp + ../../source_lcao/module_rt/td_folding.cpp + ../../source_hamilt/module_hcontainer/base_matrix.cpp + ../../source_hamilt/module_hcontainer/atom_pair.cpp + ../../source_hamilt/module_hcontainer/hcontainer.cpp ../../source_basis/module_ao/parallel_orbitals.cpp ../../source_base/parallel_2d.cpp ../../source_base/parallel_common.cpp @@ -32,4 +36,17 @@ if(ENABLE_MPI) MODULE_MD_dump_mpi MODULE_MD_dump_mpi_np4 MODULE_MD_base_mpi MODULE_MD_base_mpi_np4 PROPERTIES ENVIRONMENT "OMP_NUM_THREADS=1;MKL_NUM_THREADS=1;I_MPI_PIN=0" TIMEOUT 60) + foreach(rank_count IN ITEMS 1 4) + math(EXPR processor_count "${rank_count} * 4") + add_test( + NAME MODULE_MD_base_mpi_np${rank_count}_omp4 + COMMAND ${MPIEXEC_EXECUTABLE} ${MPIEXEC_NUMPROC_FLAG} ${rank_count} + ${MPIEXEC_PREFLAGS} $ + ${MPIEXEC_POSTFLAGS} + WORKING_DIRECTORY ${CMAKE_CURRENT_BINARY_DIR}) + set_tests_properties( + MODULE_MD_base_mpi_np${rank_count}_omp4 + PROPERTIES ENVIRONMENT "OMP_NUM_THREADS=4;MKL_NUM_THREADS=1;I_MPI_PIN=0" + PROCESSORS "${processor_count}" TIMEOUT 60) + endforeach() endif() diff --git a/source/source_md/unittests/test_md_base.cpp b/source/source_md/unittests/test_md_base.cpp index 2c57df6f11a..e8c0c2d015b 100644 --- a/source/source_md/unittests/test_md_base.cpp +++ b/source/source_md/unittests/test_md_base.cpp @@ -7,12 +7,14 @@ #include "source_cell/unitcell.h" #include "source_cell/update_cell.h" #include "source_lcao/module_rt/boundary_fix.h" +#include "source_lcao/module_rt/td_folding.h" #include #include #include #include +#include #include namespace @@ -157,7 +159,8 @@ class MdBoundaryTest : public testing::Test Complex phase_for(int index, int kpoint) const { - const double kx = kpoint == 0 ? 0.0 : 0.25; + const std::array kx_values = {{0.0, 0.25, -0.25, 0.375}}; + const double kx = kx_values[kpoint]; const double argument = -ModuleBase::TWO_PI * kx * reference_types[0].boundary_shift[index].x; const Complex exponent(0.0, argument); return std::exp(exponent); @@ -170,6 +173,141 @@ class MdBoundaryTest : public testing::Test return Complex(real, imaginary); } + std::vector fold_hybrid(const Parallel_Orbitals& orbitals, + const ModuleBase::Vector3& kpoint, + const ModuleBase::Vector3& vector_potential, + const bool wrapped) const + { + hamilt::HContainer real_space(&orbitals); + for (int row = 0; row < 4; ++row) + { + if (orbitals.get_nrow_atom(row) == 0) + { + continue; + } + for (int column = 0; column < 4; ++column) + { + if (orbitals.get_ncol_atom(column) == 0) + { + continue; + } + ModuleBase::Vector3 translation(column - row, row - column, 0); + if (wrapped) + { + translation += backing_types[0].boundary_shift[row] - backing_types[0].boundary_shift[column]; + } + hamilt::AtomPair pair(row, column, translation, &orbitals); + real_space.insert_pair(pair); + } + } + real_space.allocate(nullptr, true); + std::map, Complex> phases; + for (int index = 0; index < real_space.size_atom_pairs(); ++index) + { + hamilt::AtomPair& pair = real_space.get_atom_pair(index); + const ModuleBase::Vector3 translation = pair.get_R_index(0); + const int row = pair.get_atom_i(); + const int column = pair.get_atom_j(); + pair.get_pointer(0)[0] = matrix_value(row, column); + const ModuleBase::Vector3 direct(translation.x, translation.y, translation.z); + const ModuleBase::Vector3 cartesian = direct * backing.latvec * backing.lat0; + const double angle = vector_potential * cartesian; + const Complex exponent(0.0, angle); + phases[translation] = std::exp(exponent); + } + std::vector folded(orbitals.nloc, Complex(0.0, 0.0)); + const int rows = orbitals.get_row_size(); + module_rt::folding_HR_td(real_space, folded.data(), kpoint, vector_potential, phases, rows, 1); + return folded; + } + + void expect_hybrid_fold_correction(const bool gathered) + { + const int basis_size = 4; + const int band_count = 2; + Parallel_Orbitals orbitals; + const int init_result = orbitals.init(basis_size, basis_size, 1, cell.communicator()); + ASSERT_EQ(init_result, 0); + std::ofstream log; + const int band_result = orbitals.set_nloc_wfc_Eij(band_count, log, log); + ASSERT_EQ(band_result, 0); + orbitals.set_atomic_trace(backing.get_iat2iwt(), backing.nat, basis_size); + Parallel_Orbitals full; + full.set_serial(basis_size, basis_size); + full.set_atomic_trace(backing.get_iat2iwt(), backing.nat, basis_size); + const Parallel_Orbitals& matrix_layout = gathered ? full : orbitals; + const int matrix_size = matrix_layout.nloc; + const ct::TensorShape shape({4, matrix_size}); + ct::Tensor h(ct::DataType::DT_COMPLEX_DOUBLE, shape); + ct::Tensor s(ct::DataType::DT_COMPLEX_DOUBLE, shape); + const std::vector plane_wave_counts; + psi::Psi coefficients(4, orbitals.ncol_bands, orbitals.nrow, plane_wave_counts, true); + const std::array kx_values = {{0.0, 0.25, -0.25, 0.375}}; + std::vector> points(4); + // A nonorthogonal cell and lat0 != 1 catch missing Cartesian/unit conversions. + backing.latvec = ModuleBase::Matrix3(1.0, 0.2, 0.1, 0.0, 1.0, 0.3, 0.0, 0.0, 1.0); + backing.lat0 = 3.0; + const ModuleBase::Vector3 saved_A(0.07, -0.11, 0.05); + const ModuleBase::Vector3 next_A(-0.13, 0.17, 0.09); + std::vector> before(4); + std::vector> expected(4); + for (int kpoint = 0; kpoint < 4; ++kpoint) + { + points[kpoint].set(kx_values[kpoint], 0.0, 0.0); + before[kpoint] = fold_hybrid(matrix_layout, points[kpoint], saved_A, false); + expected[kpoint] = fold_hybrid(matrix_layout, points[kpoint], saved_A, true); + const std::vector wrong_time = fold_hybrid(matrix_layout, points[kpoint], next_A, true); + double field_time_difference = 0.0; + for (int element = 0; element < matrix_size; ++element) + { + const int offset = kpoint * matrix_size + element; + h.data()[offset] = before[kpoint][element]; + s.data()[offset] = 0.1 * before[kpoint][element]; + field_time_difference += std::abs(expected[kpoint][element] - wrong_time[element]); + } + EXPECT_GT(field_time_difference, 0.01); + coefficients.fix_k(kpoint); + for (int element = 0; element < orbitals.nrow * orbitals.ncol_bands; ++element) + { + coefficients.get_pointer()[element] = Complex(1.0 + kpoint, 0.5 - kpoint); + } + } + coefficients.fix_k(0); + module_rt::reset_matrix_boundary(backing, points, 4, saved_A, + &orbitals, h, s, &coefficients, matrix_size); + EXPECT_EQ(coefficients.get_current_k(), 0); + for (int kpoint = 0; kpoint < 4; ++kpoint) + { + for (int element = 0; element < matrix_size; ++element) + { + const int offset = kpoint * matrix_size + element; + const bool corrected = !gathered || domain.rank() == 0; + const Complex target = corrected ? expected[kpoint][element] : before[kpoint][element]; + const double h_error = std::abs(h.data()[offset] - target); + const double s_error = std::abs(s.data()[offset] - 0.1 * target); + EXPECT_LT(h_error, 1.0e-12); + EXPECT_LT(s_error, 1.0e-12); + } + coefficients.fix_k(kpoint); + for (int column = 0; column < orbitals.ncol_bands; ++column) + { + for (int row = 0; row < orbitals.nrow; ++row) + { + const int global_row = orbitals.local2global_row(row); + const ModuleBase::Vector3& shift = backing_types[0].boundary_shift[global_row]; + const ModuleBase::Vector3 direct(shift.x, shift.y, shift.z); + const ModuleBase::Vector3 cartesian = direct * backing.latvec * backing.lat0; + const double angle = ModuleBase::TWO_PI * (points[kpoint] * direct) + saved_A * cartesian; + const Complex exponent(0.0, angle); + const Complex target = Complex(1.0 + kpoint, 0.5 - kpoint) * std::exp(exponent); + const int element = column * orbitals.nrow + row; + const double error = std::abs(coefficients.get_pointer()[element] - target); + EXPECT_LT(error, 1.0e-12); + } + } + } + } + void expect_phase_correction(bool gathered) { const int basis_size = 4; @@ -182,17 +320,19 @@ class MdBoundaryTest : public testing::Test const int band_result = orbitals.set_nloc_wfc_Eij(band_count, log, log); ASSERT_EQ(band_result, 0); orbitals.set_atomic_trace(backing.get_iat2iwt(), backing.nat, basis_size); - std::vector> points(2); + std::vector> points(4); points[1].set(0.25, 0.0, 0.0); + points[2].set(-0.25, 0.0, 0.0); + points[3].set(0.375, 0.0, 0.0); const int rows = gathered ? basis_size : orbitals.get_row_size(); const int columns = gathered ? basis_size : orbitals.get_col_size(); const int matrix_size = rows * columns; - const ct::TensorShape shape({2, matrix_size}); + const ct::TensorShape shape({4, matrix_size}); ct::Tensor h(ct::DataType::DT_COMPLEX_DOUBLE, shape); ct::Tensor s(ct::DataType::DT_COMPLEX_DOUBLE, shape); const std::vector plane_wave_counts; - psi::Psi coefficients(2, orbitals.ncol_bands, orbitals.nrow, plane_wave_counts, true); - for (int kpoint = 0; kpoint < 2; ++kpoint) + psi::Psi coefficients(4, orbitals.ncol_bands, orbitals.nrow, plane_wave_counts, true); + for (int kpoint = 0; kpoint < 4; ++kpoint) { for (int column = 0; column < columns; ++column) { @@ -208,11 +348,16 @@ class MdBoundaryTest : public testing::Test coefficients.fix_k(kpoint); for (int element = 0; element < orbitals.nrow * orbitals.ncol_bands; ++element) { - coefficients.get_pointer()[element] = Complex(1.0, 0.5); + coefficients.get_pointer()[element] = Complex(1.0 + kpoint, 0.5 - kpoint); } } - module_rt::reset_matrix_boundary(backing, points, 2, &orbitals, h, s, &coefficients, matrix_size); - for (int kpoint = 0; kpoint < 2; ++kpoint) + coefficients.fix_k(0); + const ModuleBase::Vector3 no_vector_potential; + module_rt::reset_matrix_boundary(backing, points, 4, no_vector_potential, + &orbitals, h, s, &coefficients, matrix_size); + // Per-k coefficient access must not mutate a shared Psi cursor. + EXPECT_EQ(coefficients.get_current_k(), 0); + for (int kpoint = 0; kpoint < 4; ++kpoint) { for (int column = 0; column < columns; ++column) { @@ -240,7 +385,8 @@ class MdBoundaryTest : public testing::Test for (int row = 0; row < orbitals.nrow; ++row) { const int global_row = orbitals.local2global_row(row); - const Complex expected = Complex(1.0, 0.5) * std::conj(phase_for(global_row, kpoint)); + const Complex initial(1.0 + kpoint, 0.5 - kpoint); + const Complex expected = initial * std::conj(phase_for(global_row, kpoint)); const int element = column * orbitals.nrow + row; const double error = std::abs(coefficients.get_pointer()[element] - expected); EXPECT_LT(error, 1.0e-12); @@ -328,6 +474,24 @@ TEST_F(MdBoundaryTest, MatchesLegacyWrappingThroughMigrationAndRepeatedSync) } } +TEST_F(MdBoundaryTest, HybridBoundaryMatchesFreshFoldingAtSavedFieldTime) +{ + const int rank = domain.rank(); + MD_base integrator(parameters, false, true, rank, cell); + for (LocalAtom& atom : cell.owned_atoms()) + { + const int index = atom.type_index; + atom.vel = velocity_for(index, 0); + } + std::ofstream log; + integrator.first_half(log); + cell.sync_backing_unitcell(true); + EXPECT_EQ(backing_types[0].boundary_shift[0].x, -1); + EXPECT_EQ(backing_types[0].boundary_shift[1].x, 1); + expect_hybrid_fold_correction(false); + expect_hybrid_fold_correction(true); +} + TEST_F(MdBoundaryTest, InitialForceEvaluationDoesNotAdvancePositionUpdateState) { const int rank = domain.rank(); From a659109b4791acfa9b85bcdde87aae56b0eb3b0a Mon Sep 17 00:00:00 2001 From: zaotian Date: Thu, 29 Oct 2026 13:02:49 +0800 Subject: [PATCH 5/7] docs(test): record boundary-regression validation and provenance --- .../05_rtTDDFT/06_NO_dir_TDDFT/VALIDATION.md | 106 ++++++++++++++++++ 1 file changed, 106 insertions(+) create mode 100644 tests/05_rtTDDFT/06_NO_dir_TDDFT/VALIDATION.md diff --git a/tests/05_rtTDDFT/06_NO_dir_TDDFT/VALIDATION.md b/tests/05_rtTDDFT/06_NO_dir_TDDFT/VALIDATION.md new file mode 100644 index 00000000000..c5bc76070d0 --- /dev/null +++ b/tests/05_rtTDDFT/06_NO_dir_TDDFT/VALIDATION.md @@ -0,0 +1,106 @@ +# Boundary regression verification + +Recorded on 2026-10-29 for PR #8229. + +## Executable and configuration + +- Fixed source patch: `37dea1ae4ba682d476e15e67d2df54715770687c`; compiled before committing + that patch on top of `7b3a835bc70a5f1c7db9ac1fea0abf149cda6f18`. + `abacus_std_para --version` reports `v3.11.0-beta11`. +- Fixed executable SHA-256: `2dde0853955017f62dca9ddb39611367f894028b701678dba0ad4f02a9680bca`. +- Negative control: unmodified develop + `e5ad53046d4f4c4502bd04657b7fbeba60d97ad4`, + executable SHA-256 `aaff3e368eb0d6b5eb8630af53c97f82598978e3b26c77a5bbb985cf256dc1f9`. +- CPU double-precision LCAO; Intel icpx 2022.1, Intel MPI 2021.6, + MKL 2022.1, ELPA 2021.11.002, Libxc 6.2.2 and LibRI enabled. + MPI, OpenMP, testing and float FFTW were enabled; CUDA was disabled. +- Compilation used at most 20 jobs. Runtime defaults were + `OMP_NUM_THREADS=1 MKL_NUM_THREADS=1 I_MPI_PIN=0`. + Explicit OpenMP tests used four threads, at most four MPI ranks, and + `MKL_NUM_THREADS=1`; builds and runtime tests did not overlap. + +The existing configured build was updated with: + +```bash +module load abacus +export OMP_NUM_THREADS=1 MKL_NUM_THREADS=1 I_MPI_PIN=0 +cmake --build build-fix --target MODULE_MD_base_mpi MODULE_MD_func \ + MODULE_MD_verlet MODULE_MD_dump_mpi MODULE_MD_run \ + MODULE_CELL_NEIGHBOR_mdcell_migrate_mpi abacus_std_para --parallel 20 +ctest --test-dir build-fix -V \ + -R '^MODULE_MD_(base_mpi(_np4|_np[14]_omp4)?|func|verlet|dump_mpi(_np4)?|run)$|^MODULE_CELL_NEIGHBOR_(mdcell_migrate_mpi|domain_decomposition_np[14])$' +``` + +All 12 CTests passed. The new tests failed with the old boundary implementation: +the shared Psi cursor changed, and hybrid H/S disagreed with fresh folding. +The fixed tests cover Gamma and three active non-Gamma k points, both crossing +directions, nonorthogonal lattice vectors, distributed and gathered H/S, and +1/4 MPI ranks with 1/4 OpenMP threads. The hybrid oracle folds real-space +atom pairs again after shifting their lattice indices, using the saved field +endpoint; a distinct next endpoint ensures this is a time-sensitive oracle. + +## Integration commands and results + +Commands were run on fresh copies of tracked case files, with shared PP_ORB +and integrate resources. Here `fixed` is the absolute executable path and +`repo` is the checkout root. + +```bash +cd "$repo/tests/05_rtTDDFT" +bash ../integrate/Autotest.sh -a "$fixed" -n 4 -o 1 -j 1 -f CASES_CPU.txt +# Also run cases 06 and 10 for each (ranks, threads) in (1,1), (1,4), (4,1), (4,4): +bash ../integrate/Autotest.sh -a "$fixed" -n "$ranks" -o "$threads" -j 1 \ + -f CASES_CPU.txt -r '^(06_NO_dir_TDDFT|10_NO_MD_boundary_shift)$' +``` + +All 28 RT cases passed (113 property checks), as did the additional 06/10 +rank/thread combinations. Existing references were unchanged. + +The previous CI failures were rerun with the same harness using these lists: + +```bash +cd "$repo/tests/01_PW" +printf '%s\n' 093_PW_MSST 094_PW_MSST2 md_nvt 102_PW_MD_1O 103_PW_MD_2O > /tmp/pr8229-pw.list +bash ../integrate/Autotest.sh -a "$fixed" -n 4 -o 1 -j 1 -f /tmp/pr8229-pw.list +cd "$repo/tests/03_NAO_multik" +printf '%s\n' md_out_wf md_chg_extra > /tmp/pr8229-nao.list +bash ../integrate/Autotest.sh -a "$fixed" -n 2 -o 1 -j 1 -f /tmp/pr8229-nao.list +``` + +All seven cases passed. CUDA runtime was not tested locally. + +## Case 06 reference and negative control + +The replacement reference was introduced in `d373187`, using the boundary +implementation in `24fad25`. It was reproduced from the fixed source above, +with the committed INPUT/STRU/KPT, by running this in a scratch copy: + +```bash +bash ../integrate/Autotest.sh -a "$fixed" -n 4 -o 1 -j 1 \ + -f CASES_CPU.txt -r '^06_NO_dir_TDDFT$' -g +# In a separate scratch copy, with the committed reference and unmodified develop: +bash ../integrate/Autotest.sh -a "$develop" -n 4 -o 1 -j 1 \ + -f CASES_CPU.txt -r '^06_NO_dir_TDDFT$' +``` + +The standard property collector generates energy/per-atom energy, force sum, +stress sum and timing. No custom validator or tolerance change is used. +Timing is hardware dependent and is not a numerical pass criterion. + +| Run | Final energy (eV) | Force sum | Stress sum | +| --- | ---: | ---: | ---: | +| Committed reference | -30.89742338434547 | 1.490274 | 1.942380 | +| Fixed regeneration, 4 ranks / 1 thread | -30.89742338437829 | 1.490274 | 1.942380 | +| Unmodified develop, same crossing input | -13.12053220665643 | 31.992600 | 50.138461 | +| Fixed noncrossing translation control | -30.89742338437838 | 1.490274 | 1.942380 | + +Develop exited normally but the harness failed, with a 17.7768911777 eV +energy error. In the independent noncrossing control, both fractional x +coordinates were translated by `20/45` modulo one (an integer shift on +the 45-point FFT grid). Its energy agrees with the fixed crossing run to +less than 1e-9 eV; force and stress sums agree at the reported precision. + +The fresh-fold hybrid unit test verifies the boundary transformation only. +A separate Gamma-point, nonzero-file-A diagnostic still showed an +origin-translation discrepancy with fully frozen atoms, equally before and +after this patch. It is not used as a passing reference for this fix. From 0277d7ca3d99b91231b76363f8a6f7d7526b5c07 Mon Sep 17 00:00:00 2001 From: zaotian Date: Thu, 29 Oct 2026 13:10:49 +0800 Subject: [PATCH 6/7] docs(test): keep PR-specific validation outside the source tree --- .../05_rtTDDFT/06_NO_dir_TDDFT/VALIDATION.md | 106 ------------------ 1 file changed, 106 deletions(-) delete mode 100644 tests/05_rtTDDFT/06_NO_dir_TDDFT/VALIDATION.md diff --git a/tests/05_rtTDDFT/06_NO_dir_TDDFT/VALIDATION.md b/tests/05_rtTDDFT/06_NO_dir_TDDFT/VALIDATION.md deleted file mode 100644 index c5bc76070d0..00000000000 --- a/tests/05_rtTDDFT/06_NO_dir_TDDFT/VALIDATION.md +++ /dev/null @@ -1,106 +0,0 @@ -# Boundary regression verification - -Recorded on 2026-10-29 for PR #8229. - -## Executable and configuration - -- Fixed source patch: `37dea1ae4ba682d476e15e67d2df54715770687c`; compiled before committing - that patch on top of `7b3a835bc70a5f1c7db9ac1fea0abf149cda6f18`. - `abacus_std_para --version` reports `v3.11.0-beta11`. -- Fixed executable SHA-256: `2dde0853955017f62dca9ddb39611367f894028b701678dba0ad4f02a9680bca`. -- Negative control: unmodified develop - `e5ad53046d4f4c4502bd04657b7fbeba60d97ad4`, - executable SHA-256 `aaff3e368eb0d6b5eb8630af53c97f82598978e3b26c77a5bbb985cf256dc1f9`. -- CPU double-precision LCAO; Intel icpx 2022.1, Intel MPI 2021.6, - MKL 2022.1, ELPA 2021.11.002, Libxc 6.2.2 and LibRI enabled. - MPI, OpenMP, testing and float FFTW were enabled; CUDA was disabled. -- Compilation used at most 20 jobs. Runtime defaults were - `OMP_NUM_THREADS=1 MKL_NUM_THREADS=1 I_MPI_PIN=0`. - Explicit OpenMP tests used four threads, at most four MPI ranks, and - `MKL_NUM_THREADS=1`; builds and runtime tests did not overlap. - -The existing configured build was updated with: - -```bash -module load abacus -export OMP_NUM_THREADS=1 MKL_NUM_THREADS=1 I_MPI_PIN=0 -cmake --build build-fix --target MODULE_MD_base_mpi MODULE_MD_func \ - MODULE_MD_verlet MODULE_MD_dump_mpi MODULE_MD_run \ - MODULE_CELL_NEIGHBOR_mdcell_migrate_mpi abacus_std_para --parallel 20 -ctest --test-dir build-fix -V \ - -R '^MODULE_MD_(base_mpi(_np4|_np[14]_omp4)?|func|verlet|dump_mpi(_np4)?|run)$|^MODULE_CELL_NEIGHBOR_(mdcell_migrate_mpi|domain_decomposition_np[14])$' -``` - -All 12 CTests passed. The new tests failed with the old boundary implementation: -the shared Psi cursor changed, and hybrid H/S disagreed with fresh folding. -The fixed tests cover Gamma and three active non-Gamma k points, both crossing -directions, nonorthogonal lattice vectors, distributed and gathered H/S, and -1/4 MPI ranks with 1/4 OpenMP threads. The hybrid oracle folds real-space -atom pairs again after shifting their lattice indices, using the saved field -endpoint; a distinct next endpoint ensures this is a time-sensitive oracle. - -## Integration commands and results - -Commands were run on fresh copies of tracked case files, with shared PP_ORB -and integrate resources. Here `fixed` is the absolute executable path and -`repo` is the checkout root. - -```bash -cd "$repo/tests/05_rtTDDFT" -bash ../integrate/Autotest.sh -a "$fixed" -n 4 -o 1 -j 1 -f CASES_CPU.txt -# Also run cases 06 and 10 for each (ranks, threads) in (1,1), (1,4), (4,1), (4,4): -bash ../integrate/Autotest.sh -a "$fixed" -n "$ranks" -o "$threads" -j 1 \ - -f CASES_CPU.txt -r '^(06_NO_dir_TDDFT|10_NO_MD_boundary_shift)$' -``` - -All 28 RT cases passed (113 property checks), as did the additional 06/10 -rank/thread combinations. Existing references were unchanged. - -The previous CI failures were rerun with the same harness using these lists: - -```bash -cd "$repo/tests/01_PW" -printf '%s\n' 093_PW_MSST 094_PW_MSST2 md_nvt 102_PW_MD_1O 103_PW_MD_2O > /tmp/pr8229-pw.list -bash ../integrate/Autotest.sh -a "$fixed" -n 4 -o 1 -j 1 -f /tmp/pr8229-pw.list -cd "$repo/tests/03_NAO_multik" -printf '%s\n' md_out_wf md_chg_extra > /tmp/pr8229-nao.list -bash ../integrate/Autotest.sh -a "$fixed" -n 2 -o 1 -j 1 -f /tmp/pr8229-nao.list -``` - -All seven cases passed. CUDA runtime was not tested locally. - -## Case 06 reference and negative control - -The replacement reference was introduced in `d373187`, using the boundary -implementation in `24fad25`. It was reproduced from the fixed source above, -with the committed INPUT/STRU/KPT, by running this in a scratch copy: - -```bash -bash ../integrate/Autotest.sh -a "$fixed" -n 4 -o 1 -j 1 \ - -f CASES_CPU.txt -r '^06_NO_dir_TDDFT$' -g -# In a separate scratch copy, with the committed reference and unmodified develop: -bash ../integrate/Autotest.sh -a "$develop" -n 4 -o 1 -j 1 \ - -f CASES_CPU.txt -r '^06_NO_dir_TDDFT$' -``` - -The standard property collector generates energy/per-atom energy, force sum, -stress sum and timing. No custom validator or tolerance change is used. -Timing is hardware dependent and is not a numerical pass criterion. - -| Run | Final energy (eV) | Force sum | Stress sum | -| --- | ---: | ---: | ---: | -| Committed reference | -30.89742338434547 | 1.490274 | 1.942380 | -| Fixed regeneration, 4 ranks / 1 thread | -30.89742338437829 | 1.490274 | 1.942380 | -| Unmodified develop, same crossing input | -13.12053220665643 | 31.992600 | 50.138461 | -| Fixed noncrossing translation control | -30.89742338437838 | 1.490274 | 1.942380 | - -Develop exited normally but the harness failed, with a 17.7768911777 eV -energy error. In the independent noncrossing control, both fractional x -coordinates were translated by `20/45` modulo one (an integer shift on -the 45-point FFT grid). Its energy agrees with the fixed crossing run to -less than 1e-9 eV; force and stress sums agree at the reported precision. - -The fresh-fold hybrid unit test verifies the boundary transformation only. -A separate Gamma-point, nonzero-file-A diagnostic still showed an -origin-translation discrepancy with fully frozen atoms, equally before and -after this patch. It is not used as a passing reference for this fix. From aeba8de56e2bf54c068a77c58792529e27d8f204 Mon Sep 17 00:00:00 2001 From: zaotian Date: Thu, 29 Oct 2026 15:49:52 +0800 Subject: [PATCH 7/7] test(rt): cover saved-field timing through the MD runner --- source/source_md/unittests/test_md_base.cpp | 23 +++++++++++++++---- tests/05_rtTDDFT/18_NO_hyb_TDDFT/INPUT | 18 ++++++++++----- tests/05_rtTDDFT/18_NO_hyb_TDDFT/README | 1 + tests/05_rtTDDFT/18_NO_hyb_TDDFT/STRU | 12 +++++----- .../18_NO_hyb_TDDFT/current_tot.txt.ref | 18 ++++++--------- tests/05_rtTDDFT/18_NO_hyb_TDDFT/result.ref | 7 +++--- .../18_NO_hyb_TDDFT/vector_pot.txt.ref | 8 +++++++ 7 files changed, 56 insertions(+), 31 deletions(-) create mode 100644 tests/05_rtTDDFT/18_NO_hyb_TDDFT/vector_pot.txt.ref diff --git a/source/source_md/unittests/test_md_base.cpp b/source/source_md/unittests/test_md_base.cpp index e8c0c2d015b..c0fa9644d74 100644 --- a/source/source_md/unittests/test_md_base.cpp +++ b/source/source_md/unittests/test_md_base.cpp @@ -191,7 +191,9 @@ class MdBoundaryTest : public testing::Test { continue; } - ModuleBase::Vector3 translation(column - row, row - column, 0); + const int translation_x = column - row; + const int translation_y = row - column; + ModuleBase::Vector3 translation(translation_x, translation_y, 0); if (wrapped) { translation += backing_types[0].boundary_shift[row] - backing_types[0].boundary_shift[column]; @@ -267,9 +269,12 @@ class MdBoundaryTest : public testing::Test } EXPECT_GT(field_time_difference, 0.01); coefficients.fix_k(kpoint); + const double coefficient_real = 1.0 + kpoint; + const double coefficient_imaginary = 0.5 - kpoint; + const Complex initial(coefficient_real, coefficient_imaginary); for (int element = 0; element < orbitals.nrow * orbitals.ncol_bands; ++element) { - coefficients.get_pointer()[element] = Complex(1.0 + kpoint, 0.5 - kpoint); + coefficients.get_pointer()[element] = initial; } } coefficients.fix_k(0); @@ -289,6 +294,9 @@ class MdBoundaryTest : public testing::Test EXPECT_LT(s_error, 1.0e-12); } coefficients.fix_k(kpoint); + const double coefficient_real = 1.0 + kpoint; + const double coefficient_imaginary = 0.5 - kpoint; + const Complex initial(coefficient_real, coefficient_imaginary); for (int column = 0; column < orbitals.ncol_bands; ++column) { for (int row = 0; row < orbitals.nrow; ++row) @@ -299,7 +307,7 @@ class MdBoundaryTest : public testing::Test const ModuleBase::Vector3 cartesian = direct * backing.latvec * backing.lat0; const double angle = ModuleBase::TWO_PI * (points[kpoint] * direct) + saved_A * cartesian; const Complex exponent(0.0, angle); - const Complex target = Complex(1.0 + kpoint, 0.5 - kpoint) * std::exp(exponent); + const Complex target = initial * std::exp(exponent); const int element = column * orbitals.nrow + row; const double error = std::abs(coefficients.get_pointer()[element] - target); EXPECT_LT(error, 1.0e-12); @@ -346,9 +354,12 @@ class MdBoundaryTest : public testing::Test } } coefficients.fix_k(kpoint); + const double coefficient_real = 1.0 + kpoint; + const double coefficient_imaginary = 0.5 - kpoint; + const Complex initial(coefficient_real, coefficient_imaginary); for (int element = 0; element < orbitals.nrow * orbitals.ncol_bands; ++element) { - coefficients.get_pointer()[element] = Complex(1.0 + kpoint, 0.5 - kpoint); + coefficients.get_pointer()[element] = initial; } } coefficients.fix_k(0); @@ -380,12 +391,14 @@ class MdBoundaryTest : public testing::Test } } coefficients.fix_k(kpoint); + const double coefficient_real = 1.0 + kpoint; + const double coefficient_imaginary = 0.5 - kpoint; + const Complex initial(coefficient_real, coefficient_imaginary); for (int column = 0; column < orbitals.ncol_bands; ++column) { for (int row = 0; row < orbitals.nrow; ++row) { const int global_row = orbitals.local2global_row(row); - const Complex initial(1.0 + kpoint, 0.5 - kpoint); const Complex expected = initial * std::conj(phase_for(global_row, kpoint)); const int element = column * orbitals.nrow + row; const double error = std::abs(coefficients.get_pointer()[element] - expected); diff --git a/tests/05_rtTDDFT/18_NO_hyb_TDDFT/INPUT b/tests/05_rtTDDFT/18_NO_hyb_TDDFT/INPUT index 492e21dc163..eef00dd0998 100644 --- a/tests/05_rtTDDFT/18_NO_hyb_TDDFT/INPUT +++ b/tests/05_rtTDDFT/18_NO_hyb_TDDFT/INPUT @@ -6,10 +6,10 @@ INPUT_PARAMETERS calculation md esolver_type tddft md_type nve -md_nstep 1 -estep_per_md 10 -td_dt 0.005 -md_tfirst 0 +md_nstep 4 +estep_per_md 2 +td_dt 0.025 +md_tfirst -1 # rt-TDDFT parameters @@ -36,11 +36,17 @@ gamma_only 0 # electronic structure calculations ecutwfc 20 -scf_nmax 50 -scf_thr 1e-6 +scf_nmax 100 +scf_thr 1e-8 # charge mixing mixing_type broyden mixing_beta 0.7 mixing_gg0 0.0 + +md_dt 0.05 +init_vel 1 +symmetry -1 +md_dumpfreq 1 +out_vecpot 1 diff --git a/tests/05_rtTDDFT/18_NO_hyb_TDDFT/README b/tests/05_rtTDDFT/18_NO_hyb_TDDFT/README index a47d3411f75..e891923a377 100644 --- a/tests/05_rtTDDFT/18_NO_hyb_TDDFT/README +++ b/tests/05_rtTDDFT/18_NO_hyb_TDDFT/README @@ -1 +1,2 @@ Test RT-TDDFT hybrid gauge output and parameter td_dt for periodic system +Also checks saved-field timing across a periodic boundary with two electronic steps per MD step. diff --git a/tests/05_rtTDDFT/18_NO_hyb_TDDFT/STRU b/tests/05_rtTDDFT/18_NO_hyb_TDDFT/STRU index 822d2c7cc68..e9eebc358ca 100644 --- a/tests/05_rtTDDFT/18_NO_hyb_TDDFT/STRU +++ b/tests/05_rtTDDFT/18_NO_hyb_TDDFT/STRU @@ -13,10 +13,10 @@ LATTICE_VECTORS 0.5 0.5 0.0 ATOMIC_POSITIONS -Cartesian //Cartesian or Direct coordinate. +Direct -Si // Element type -0.0 // magnetism -2 // number of atoms -0.00 0.00 0.00 0 0 0 -0.25 0.25 0.25 1 1 1 +Si +0.0 +2 +0.9995 0.0 0.0 m 1 1 1 v 0.0 0.001 0.001 +0.2495 0.25 0.25 m 1 1 1 v 0.0 -0.001 -0.001 diff --git a/tests/05_rtTDDFT/18_NO_hyb_TDDFT/current_tot.txt.ref b/tests/05_rtTDDFT/18_NO_hyb_TDDFT/current_tot.txt.ref index 7f49c2cedde..0ecf32cb2cc 100644 --- a/tests/05_rtTDDFT/18_NO_hyb_TDDFT/current_tot.txt.ref +++ b/tests/05_rtTDDFT/18_NO_hyb_TDDFT/current_tot.txt.ref @@ -1,11 +1,7 @@ -1 -1.2080528880846813e-19 1.2080528880847670e-19 1.2080528880844887e-19 -2 2.1944657158836502e-07 -2.1944657277845298e-07 3.1856225789490216e-06 -3 2.3687980954535415e-07 -2.3687978488285760e-07 8.0274188859243376e-06 -4 2.8641942462454050e-08 -2.8641948921304256e-08 9.6346681720893599e-06 -5 1.2615316193284446e-08 -1.2615330098784816e-08 9.5047117782733680e-06 -6 -9.9299854444834876e-09 9.9299763872357555e-09 9.3134367812586245e-06 -7 -3.7912652295188386e-08 3.7912648596369048e-08 9.0666648293070098e-06 -8 -7.0591239361918577e-08 7.0591241073221420e-08 8.7705306596311231e-06 -9 -1.0760302590107121e-07 1.0760303271470401e-07 8.4311195836446455e-06 -10 -1.4885404710402708e-07 1.4885405845888459e-07 8.0543604498738173e-06 -11 -1.9431312768314243e-07 1.9431314289794967e-07 7.6460963157211179e-06 +1 -1.2080517284530851e-19 1.2080531731065871e-19 1.2080531731065948e-19 +2 1.2451440495663431e-06 -1.2104466814246842e-06 1.5049160765224529e-05 +3 9.8538281614216140e-07 -9.4914223443828920e-07 3.6655451677813790e-05 +4 -5.2766514366548575e-07 5.9153652074363848e-07 3.9125916000482234e-05 +5 -1.6767814605672313e-06 1.7344158416499257e-06 3.0008656483392178e-05 +6 -3.2561979854203605e-06 3.2774917151898570e-06 1.8545534269939497e-05 +7 -5.0827836916354501e-06 5.0738257431279963e-06 5.8331186497094963e-06 diff --git a/tests/05_rtTDDFT/18_NO_hyb_TDDFT/result.ref b/tests/05_rtTDDFT/18_NO_hyb_TDDFT/result.ref index a31830fc496..6eeefc4b913 100644 --- a/tests/05_rtTDDFT/18_NO_hyb_TDDFT/result.ref +++ b/tests/05_rtTDDFT/18_NO_hyb_TDDFT/result.ref @@ -1,4 +1,5 @@ -etotref -202.3030374605464 -etotperatomref -101.1515187303 +etotref -202.2977241326445 +etotperatomref -101.1488620663 CompareCurrent_pass 0 -totaltimeref 5.26 +CompareVectorPot_pass 0 +totaltimeref 5.12 diff --git a/tests/05_rtTDDFT/18_NO_hyb_TDDFT/vector_pot.txt.ref b/tests/05_rtTDDFT/18_NO_hyb_TDDFT/vector_pot.txt.ref new file mode 100644 index 00000000000..7af02933c63 --- /dev/null +++ b/tests/05_rtTDDFT/18_NO_hyb_TDDFT/vector_pot.txt.ref @@ -0,0 +1,8 @@ +#istep A_x A_y A_z # endpoint values, Hartree atomic units +1 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 +2 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00 +3 0.0000000000000000e+00 0.0000000000000000e+00 -1.0049513776379543e-03 +4 0.0000000000000000e+00 0.0000000000000000e+00 -2.0099027552759087e-03 +5 0.0000000000000000e+00 0.0000000000000000e+00 -2.0099027552759087e-03 +6 0.0000000000000000e+00 0.0000000000000000e+00 -2.0099027552759087e-03 +7 0.0000000000000000e+00 0.0000000000000000e+00 -2.0099027552759087e-03