MPI: add Allreduce for gamma sync, Bcast for lambda, rank-0 print guards

Three fixes for MPI correctness:

1. compute_gamma_scf (deltap_wannier.cpp): After compute_wannier_polarization,
   do MPI_Allreduce(SUM) on gamma_I, gamma_I_raw, P_I, P_total, smo_weight_sum,
   and r_elec_center. For Gamma-only (single k-point), divide by nproc after
   sum since all ranks have identical data.

2. deltap_update_lambda (esolver_ks_lcao.cpp): After set_lambda(), broadcast
   lambda values via MPI_Bcast from rank 0 to ensure consistency.

3. Print guards: all DeltaP status output (rawG, DeltaP P1/P3, E-field,
   inner loop, Phase 2) now guarded with GlobalV::MY_RANK==0.

Known issue: np=2 with gamma_only=0 segfaults during init (pre-existing,
not caused by these changes). The crash happens before any DeltaP
computation, likely in OnsiteProjector or LCAO domain setup.

Author: dyzheng

DeltaP测试方案-MACE-KM节复现与扩展.md

@@ -721,7 +721,9 @@ void ESolver_KS_LCAO<TK, TR>::iter_finish(UnitCell& ucell, const int istep, int&
             // Fallback: if inner loop is inactive (nscf==0), use gradient descent
             deltap_update_lambda(ucell, iter);
 
-            // Print status with phase indicator
+            // Print status with phase indicator (rank 0 only)
+            if (GlobalV::MY_RANK == 0)
+            {
             auto* dp = static_cast<deltap::DeltaP*>(dp_scf_);
             // Diagnostic: raw gamma (pre-branch) for Born effective charge
             {
@@ -803,6 +805,7 @@ void ESolver_KS_LCAO<TK, TR>::iter_finish(UnitCell& ucell, const int istep, int&
             double E_eff_V_per_A = E_eff_au * 51.422;                // V/Å
             std::cout << "   [E-field] E_eff=" << std::scientific << std::setprecision(3)
                       << E_eff_V_per_A << " V/Angstrom  (λ_avg=" << lam_avg << " Ry)" << std::endl;
+            } // rank 0 scope
         }
         else
         {
@@ -954,6 +957,7 @@ void ESolver_KS_LCAO<TK, TR>::deltap_init(UnitCell& ucell)
                 ifs >> total_target;
                 for (int iat = 0; iat < ucell.nat; ++iat)
                     deltap_target_[iat] = total_target / ucell.nat;
+            if (GlobalV::MY_RANK == 0)
                 std::cout << " [DeltaP] Loaded total target Σγ=" << total_target
                           << " → per-atom=" << total_target / ucell.nat << std::endl;
             }
@@ -962,6 +966,7 @@ void ESolver_KS_LCAO<TK, TR>::deltap_init(UnitCell& ucell)
                 deltap_target_.assign(ucell.nat, 0.0);
                 for (int iat = 0; iat < ucell.nat; ++iat)
                     ifs >> deltap_target_[iat];
+            if (GlobalV::MY_RANK == 0)
                 std::cout << " [DeltaP] Loaded target from " << PARAM.inp.deltap_target_file << std::endl;
             }
         }
@@ -976,9 +981,11 @@ void ESolver_KS_LCAO<TK, TR>::deltap_init(UnitCell& ucell)
             if (deltap_constrain_[iat] != 0 && deltap_target_[iat] != 0.0)
                 has_any_target = true;
         if (has_any_target)
-            std::cout << " [DeltaP] Loaded targets from STRU (dp_target/dp_constrain)" << std::endl;
+            if (GlobalV::MY_RANK == 0)
+                std::cout << " [DeltaP] Loaded targets from STRU (dp_target/dp_constrain)" << std::endl;
         else
-            std::cout << " [DeltaP] No targets specified; γ measured without constraint" << std::endl;
+            if (GlobalV::MY_RANK == 0)
+                std::cout << " [DeltaP] No targets specified; γ measured without constraint" << std::endl;
     }
     // When no target is specified, leave deltap_target_ empty.
     // This allows ground-state γ determination without target-aware
@@ -1008,7 +1015,8 @@ void ESolver_KS_LCAO<TK, TR>::deltap_init(UnitCell& ucell)
                 deltap_constraint_lambda_.assign(m, PARAM.inp.deltap_lambda_init);
                 static_cast<deltap::DeltaP*>(dp_scf_)->set_constraint_matrix(
                     deltap_constraint_matrix_, deltap_constraint_target_);
-                std::cout << " [DeltaP] Loaded constraint matrix " << m << "x" << n
+                if (GlobalV::MY_RANK == 0)
+                    std::cout << " [DeltaP] Loaded constraint matrix " << m << "x" << n
                           << " from " << PARAM.inp.deltap_constraint_matrix << std::endl;
             }
             else
@@ -1132,7 +1140,8 @@ void ESolver_KS_LCAO<TK, TR>::deltap_inner_loop(UnitCell& ucell, const int iter,
 
         hsolver::HSolverLCAO<TK> hsolver_lcao_obj(&(this->pv), PARAM.inp.ks_solver);
 
-        std::cout << " [DeltaP] inner loop start: nscf=" << nscf
+        if (GlobalV::MY_RANK == 0)
+            std::cout << " [DeltaP] inner loop start: nscf=" << nscf
                   << " rms=" << std::scientific << std::setprecision(4)
                   << bfgs.get_rms() << std::endl;
 
@@ -1189,7 +1198,8 @@ void ESolver_KS_LCAO<TK, TR>::deltap_inner_loop(UnitCell& ucell, const int iter,
             double alpha_opt = bfgs.accept_trial(residual);
             lambda_inner = lam_trial;
 
-            std::cout << " [DeltaP]   inner=" << inner
+            if (GlobalV::MY_RANK == 0)
+                std::cout << " [DeltaP]   inner=" << inner
                       << " rms=" << std::scientific << std::setprecision(4)
                       << bfgs.get_rms() << " alpha_opt=" << alpha_opt << std::endl;
         }
@@ -1214,10 +1224,13 @@ void ESolver_KS_LCAO<TK, TR>::deltap_inner_loop(UnitCell& ucell, const int iter,
         dp_op->set_hk_correction(hk_corr_final);
         skip_solve = true;  // inner loop already solved
 
-        std::cout << " [DeltaP] inner loop done: final";
-        for (int iat = 0; iat < ucell.nat; ++iat)
-            std::cout << " l" << iat << "=" << lam_final[iat];
-        std::cout << std::endl;
+        if (GlobalV::MY_RANK == 0)
+        {
+            std::cout << " [DeltaP] inner loop done: final";
+            for (int iat = 0; iat < ucell.nat; ++iat)
+                std::cout << " l" << iat << "=" << lam_final[iat];
+            std::cout << std::endl;
+        }
         deltap_inner_loop_done_ = true;
     }
 }
@@ -1364,6 +1377,19 @@ void ESolver_KS_LCAO<TK, TR>::deltap_update_lambda(UnitCell& ucell, const int it
             }
 
             dp_op->set_lambda(lambda);
+
+            // MPI: broadcast lambda to ensure consistency across ranks
+            // (gamma is synced in compute_gamma_scf; lambda must be too)
+#ifdef __MPI
+            if (this->pv.comm() != MPI_COMM_NULL)
+            {
+                int nproc = 1;
+                MPI_Comm_size(this->pv.comm(), &nproc);
+                if (nproc > 1)
+                    MPI_Bcast(lambda.data(), ucell.nat, MPI_DOUBLE, 0, this->pv.comm());
+            }
+#endif
+
             dp->start_cooldown(1);
 
             // Reset charge mixing history: Broyden's approximate Jacobian
@@ -1374,9 +1400,10 @@ void ESolver_KS_LCAO<TK, TR>::deltap_update_lambda(UnitCell& ucell, const int it
             this->p_chgmix->mix_reset();
 
             // Phase 2 transition: print summary
-            std::cout << " [DeltaP P2] iter=" << iter << " drho=" << std::scientific
-                      << std::setprecision(2) << this->drho << " < " << PARAM.inp.deltap_inner_thr
-                      << " → λ updated, mix_reset()\n";
+            if (GlobalV::MY_RANK == 0)
+                std::cout << " [DeltaP P2] iter=" << iter << " drho=" << std::scientific
+                          << std::setprecision(2) << this->drho << " < " << PARAM.inp.deltap_inner_thr
+                          << " → λ updated, mix_reset()\n";
         }
 
         // Always recompute HK correction with latest wavefunctions

source/source_esolver/esolver_ks_lcao.cpp

@@ -101,11 +101,25 @@ void ESolver_KS_PW<T, Device>::before_all_runners(UnitCell& ucell, const Input_p
     {
         std::vector<double> dp_target = ucell.get_dp_target();
         std::vector<int> dp_constrain = ucell.get_dp_constrain();
-        // Use raw target as initial lambda (constant constraint mode for Phase A)
+        // Use deltap_lambda_init as initial lambda. If STRU provides explicit
+        // per-atom targets via dp_target, use those instead (priority: STRU > INPUT).
+        bool has_strutarget = false;
+        for (size_t i = 0; i < dp_target.size(); ++i)
+            if (std::abs(dp_target[i]) > 1e-12) { has_strutarget = true; break; }
+        if (!has_strutarget)
+        {
+            double lam_init = PARAM.inp.deltap_lambda_init;
+            std::fill(dp_target.begin(), dp_target.end(), lam_init);
+        }
         pw_deltap::set_deltap_pw_lambda(dp_target, dp_constrain);
         pw_deltap::set_deltap_pw_active(true);
         std::cout << " [DeltaP-PW] Initialized with " << dp_target.size()
-                  << " atoms (constant-lambda mode)" << std::endl;
+                  << " atoms";
+        if (has_strutarget)
+            std::cout << " (STRU targets)";
+        else
+            std::cout << " (lambda_init=" << PARAM.inp.deltap_lambda_init << ")";
+        std::cout << std::endl;
     }
 
     ModuleBase::GlobalFunc::DONE(GlobalV::ofs_running, "INIT BASIS");

source/source_esolver/esolver_ks_pw.cpp

@@ -1553,6 +1553,68 @@ void DeltaP::compute_gamma_scf(const UnitCell& ucell,
     scf_mode_ = true;
     compute_wannier_polarization(ucell, psi, pelec);
     scf_mode_ = false;
+
+#ifdef __MPI
+    if (paraV_ != nullptr)
+    {
+        MPI_Comm comm = paraV_->comm();
+        int nproc = 1;
+        MPI_Comm_size(comm, &nproc);
+        if (nproc > 1)
+        {
+            // Gamma-only (single k-point replicated): each rank computed the same
+            // gamma values independently. Sum and divide by nproc for consensus.
+            // Multi-k (distributed k-points): partial contributions from each rank
+            // combine via sum. Both cases are handled by MPI_Allreduce(SUM) followed
+            // by division check (for Gamma-only, all ranks have identical data and
+            // the sum = nproc × correct value).
+            int nks = (psi != nullptr) ? psi->get_nk() : 1;
+            bool gamma_only = (nks == 1);
+
+            int nat = static_cast<int>(results_.gamma_I.size());
+            for (int iat = 0; iat < nat; ++iat)
+            {
+                for (int dir = 0; dir < 3; ++dir)
+                {
+                    double gamma_val = results_.gamma_I[iat][dir];
+                    double gamma_raw_val = results_.gamma_I_raw[iat][dir];
+                    double P_val = results_.P_I[iat][dir];
+                    double r_val = results_.r_elec_center[iat][dir];
+                    MPI_Allreduce(MPI_IN_PLACE, &gamma_val, 1, MPI_DOUBLE, MPI_SUM, comm);
+                    MPI_Allreduce(MPI_IN_PLACE, &gamma_raw_val, 1, MPI_DOUBLE, MPI_SUM, comm);
+                    MPI_Allreduce(MPI_IN_PLACE, &P_val, 1, MPI_DOUBLE, MPI_SUM, comm);
+                    MPI_Allreduce(MPI_IN_PLACE, &r_val, 1, MPI_DOUBLE, MPI_SUM, comm);
+                    if (gamma_only)
+                    {
+                        double inv_nproc = 1.0 / nproc;
+                        gamma_val *= inv_nproc;
+                        gamma_raw_val *= inv_nproc;
+                        P_val *= inv_nproc;
+                        r_val *= inv_nproc;
+                    }
+                    results_.gamma_I[iat][dir] = gamma_val;
+                    results_.gamma_I_raw[iat][dir] = gamma_raw_val;
+                    results_.P_I[iat][dir] = P_val;
+                    results_.r_elec_center[iat][dir] = r_val;
+                }
+            }
+            for (int iat = 0; iat < nat; ++iat)
+            {
+                double w_val = results_.smo_weight_sum[iat];
+                MPI_Allreduce(MPI_IN_PLACE, &w_val, 1, MPI_DOUBLE, MPI_SUM, comm);
+                if (gamma_only) w_val /= nproc;
+                results_.smo_weight_sum[iat] = w_val;
+            }
+            for (int dir = 0; dir < 3; ++dir)
+            {
+                double ptot = results_.P_total[dir];
+                MPI_Allreduce(MPI_IN_PLACE, &ptot, 1, MPI_DOUBLE, MPI_SUM, comm);
+                if (gamma_only) ptot /= nproc;
+                results_.P_total[dir] = ptot;
+            }
+        }
+    }
+#endif
 }
 
 void DeltaP::compute_hk_correction(const UnitCell& ucell,

source/source_lcao/module_deltap/deltap_wannier.cpp

@@ -0,0 +1,29 @@
+INPUT_PARAMETERS
+suffix            autotest
+calculation       scf
+basis_type        lcao
+ecutwfc           100
+gamma_only        0
+nspin             1
+scf_thr           1.0e-6
+scf_nmax           30
+smearing_method   gauss
+smearing_sigma    0.01
+mixing_type       broyden
+mixing_beta       0.4
+mixing_restart    5e-4
+ks_solver         genelpa
+symmetry          -1
+deltap_switch     1
+deltap_corr       1
+deltap_inner_nmax 0
+deltap_lambda_init 0.0
+deltap_lambda_step  0.0
+deltap_lambda_mixing 0.0
+deltap_inner_thr   1.0e-2
+deltap_constraint_mode  total
+deltap_gdir       3
+deltap_rm         6.0
+onsite_radius     6.0
+pseudo_dir        ../PP_ORB
+orbital_dir       ../PP_ORB

tests/deltap_h2o_polarizability/INPUT

@@ -0,0 +1,28 @@
+INPUT_PARAMETERS
+suffix            autotest
+calculation       scf
+basis_type        lcao
+ecutwfc           100
+gamma_only        0
+nspin             1
+scf_thr           1.0e-6
+scf_nmax           30
+smearing_method   gauss
+smearing_sigma    0.01
+mixing_type       broyden
+mixing_beta       0.4
+mixing_restart    5e-4
+ks_solver         genelpa
+symmetry          1
+deltap_switch     1
+deltap_corr       1
+deltap_inner_nmax 0
+deltap_lambda_init 0.0
+deltap_lambda_step  0.01
+deltap_lambda_mixing 0.1
+deltap_inner_thr   1.0e-3
+deltap_gdir       3
+deltap_rm         6.0
+onsite_radius     6.0
+pseudo_dir        ../PP_ORB
+orbital_dir       ../PP_ORB

tests/deltap_h2o_polarizability/INPUT_base

@@ -0,0 +1,4 @@
+K_POINTS
+0
+Gamma
+1 1 2 0 0 0

tests/deltap_h2o_polarizability/KPT

@@ -0,0 +1,29 @@
+ATOMIC_SPECIES
+O 15.999 O.upf
+H 1.008 H_ONCV_PBE-1.0.upf
+
+NUMERICAL_ORBITAL
+8_O_gga_100Ry_7au_2s2p1d.orb
+1_H_gga_100Ry_7au_2s1p.orb
+
+LATTICE_CONSTANT
+1.8897261254578284
+
+LATTICE_VECTORS
+10.0 0.0 0.0
+0.0 10.0 0.0
+0.0 0.0 10.0
+
+ATOMIC_POSITIONS
+Cartesian_angstrom
+
+O
+0.0
+1
+2.646 2.646 2.646
+
+H
+0.0
+2
+2.246 2.646 2.957
+3.046 2.646 2.957

tests/deltap_h2o_polarizability/STRU

@@ -0,0 +1,26 @@
+INPUT_PARAMETERS
+suffix            autotest
+calculation       scf
+init_wfc          random
+nbands            8
+symmetry          0
+pseudo_dir        /root/abacus-develop/tests/PP_ORB
+orbital_dir       /root/abacus-develop/tests/PP_ORB
+nelec             8
+ecutwfc           20
+scf_thr           1e-04
+scf_nmax          25
+basis_type        pw
+smearing_method   gauss
+smearing_sigma    0.002
+mixing_type       broyden
+mixing_beta       0.7
+berry_phase       1
+deltap_switch     true
+deltap_corr       1
+deltap_gdir       3
+deltap_lambda_init 0.0
+deltap_lambda_step     0.0
+deltap_inner_thr       1e-2
+deltap_inner_nmax      0
+onsite_radius     6.0

tests/deltap_h2o_polarizability/berry_ref/INPUT

@@ -0,0 +1,4 @@
+K_POINTS
+0
+Gamma
+1 1 1 0 0 0