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..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) { @@ -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,22 @@ 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]); } } + if (positions_updated) + { + backing_unitcell_->ionic_position_updated = true; + } return; } #endif @@ -263,24 +268,17 @@ 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; } + 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_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 bc4c5bbc448..00e7e9bf7d0 100644 --- a/source/source_esolver/esolver_ks_lcao_tddft.cpp +++ b/source/source_esolver/esolver_ks_lcao_tddft.cpp @@ -155,6 +155,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) @@ -210,7 +216,16 @@ 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, + boundary_A, + &(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..641c3ff9964 100644 --- a/source/source_lcao/module_rt/boundary_fix.cpp +++ b/source/source_lcao/module_rt/boundary_fix.cpp @@ -1,12 +1,15 @@ #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 ModuleBase::Vector3& cart_At_last, const Parallel_Orbitals* pv, ct::Tensor& hk_last, ct::Tensor& sk_last, @@ -15,6 +18,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 +43,13 @@ 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; + // 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 @@ -37,10 +57,20 @@ 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); - psi_last->fix_k(ik); - module_rt::boundary_shift_c(phase, psi_last[0].get_pointer(), 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); + } + 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); + } } } } @@ -96,4 +126,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..0b61b5d0cfb 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,9 @@ 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] 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 @@ -27,7 +29,9 @@ 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 ModuleBase::Vector3& cart_At_last, const Parallel_Orbitals* pv, ct::Tensor& hk_last, ct::Tensor& sk_last, @@ -61,4 +65,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/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/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..4b84963538f 100644 --- a/source/source_md/unittests/CMakeLists.txt +++ b/source/source_md/unittests/CMakeLists.txt @@ -13,7 +13,40 @@ 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_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 + ../../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) + 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 new file mode 100644 index 00000000000..c0fa9644d74 --- /dev/null +++ b/source/source_md/unittests/test_md_base.cpp @@ -0,0 +1,570 @@ +#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" +#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 "source_lcao/module_rt/td_folding.h" + +#include + +#include +#include +#include +#include +#include + +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: + 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 bool positions_updated) + { + const Atom& actual = backing_types[0]; + const Atom& reference = reference_types[0]; + EXPECT_EQ(backing.ionic_position_updated, positions_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 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); + } + + 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); + } + + 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; + } + 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]; + } + 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); + 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] = initial; + } + } + 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); + 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 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 = 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); + } + } + } + } + + 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(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({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); + for (int kpoint = 0; kpoint < 4; ++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); + 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] = initial; + } + } + 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) + { + 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); + 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 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); + } + } + } + } +}; + +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(true); + expect_backing_state(true); + 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) + { + 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); + 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, 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(); + 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); + } + 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 + +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/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/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