Restore boundary-shift propagation for LCAO RT-TDDFT MD

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

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.

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<double> frac(packed.frac[0], packed.frac[1], packed.frac[2]);
     const ModuleBase::Vector3<double> 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<double> cart = image_frac * latvec_;
     const ModuleBase::Vector3<double> vel(packed.vel[0], packed.vel[1], packed.vel[2]);
     const ModuleBase::Vector3<double> force(packed.force[0], packed.force[1], packed.force[2]);
     const ModuleBase::Vector3<int> 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<double> vel(packed.vel[0], packed.vel[1], packed.vel[2]);
     const ModuleBase::Vector3<double> force(packed.force[0], packed.force[1], packed.force[2]);
     const ModuleBase::Vector3<int> 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;
 }

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];

source/source_cell/mdcell.cpp

@@ -177,14 +177,18 @@ 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)
     {
         return;
     }
 
     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<double> cart(3 * nat_, 0.0);
-        std::vector<double> frac(3 * nat_, 0.0);
-        std::vector<double> vel(3 * nat_, 0.0);
-        std::vector<int> mbl(3 * nat_, 0);
+        const int coordinate_count = 3 * nat_;
+        std::vector<double> cart(coordinate_count, 0.0);
+        std::vector<double> frac(coordinate_count, 0.0);
+        std::vector<double> vel(coordinate_count, 0.0);
+        std::vector<double> displacement(coordinate_count, 0.0);
+        std::vector<int> boundary_shift(coordinate_count, 0);
+        std::vector<int> mbl(coordinate_count, 0);
         std::vector<int> 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,52 +245,40 @@ 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<double> 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
 
     for (std::size_t i = 0; i < owned_atoms_.size(); ++i)
     {
         const LocalAtom& atom = owned_atoms_[i];
-        ModuleBase::Vector3<double> 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

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<LocalAtom>& owned_atoms() { return owned_atoms_; }
     const std::vector<LocalAtom>& owned_atoms() const { return owned_atoms_; }

source/source_cell/module_neighlist/local_atom.h

@@ -20,6 +20,8 @@ struct LocalAtom
     ModuleBase::Vector3<double> frac;
     ModuleBase::Vector3<double> vel;
     ModuleBase::Vector3<double> force;
+    ModuleBase::Vector3<double> displacement; ///< Unwrapped fractional displacement in the current MD step.
+    ModuleBase::Vector3<int> boundary_shift; ///< Lattice translation applied in the current MD step.
     ModuleBase::Vector3<int> 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),

source/source_esolver/esolver_ks_lcao_tddft.cpp

@@ -155,6 +155,12 @@ void ESolver_KS_LCAO_TDDFT<TR, Device>::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<double> 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<TR, Device>::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++)
     {