From 65d533bd80bf327f6cca6aad8964255f27c94d76 Mon Sep 17 00:00:00 2001
From: Chen Nuo <49788094+Cstandardlib@users.noreply.github.com>
Date: Sun, 2 Aug 2026 01:28:08 +0800
Subject: [PATCH] Fix redundant density symmetrization in OFDFT

---
 source/source_esolver/esolver_of.cpp | 24 ++++++++++++++----------
 1 file changed, 14 insertions(+), 10 deletions(-)

diff --git a/source/source_esolver/esolver_of.cpp b/source/source_esolver/esolver_of.cpp
index 6437ace50c..034bb16503 100644
--- a/source/source_esolver/esolver_of.cpp
+++ b/source/source_esolver/esolver_of.cpp
@@ -236,30 +236,34 @@ void ESolver_OF::before_opt(const int istep, UnitCell& ucell)
 
     elecstate::init_scf(ucell, Pgrid, sf.strucFac, locpp.numeric, istep, PARAM.globalv.global_out_dir, PARAM.inp, this->pelec);
 
-    Symmetry_rho::symmetrize_rho(PARAM.inp.nspin, this->chr, this->pw_rho, ucell.symm);
-
-    for (int is = 0; is < PARAM.inp.nspin; ++is)
+    const int nspin = PARAM.inp.nspin;
+    if (PARAM.inp.init_chg == "file")
     {
-        if (PARAM.inp.init_chg != "file")
+        Symmetry_rho::symmetrize_rho(nspin, this->chr, this->pw_rho, ucell.symm);
+        for (int is = 0; is < nspin; ++is)
         {
             for (int ibs = 0; ibs < this->pw_rho->nrxx; ++ibs)
             {
-                // Here we initialize rho to be uniform,
-                // because the rho got by pot.init_pot -> Charge::atomic_rho may contain minus elements.
-                this->chr.rho[is][ibs] = this->nelec_[is] / ucell.omega;
                 this->pphi_[is][ibs] = sqrt(this->chr.rho[is][ibs]);
             }
         }
-        else
+    }
+    else
+    {
+        // Non-file densities are replaced with a uniform density, so
+        // symmetrizing them would only add an unnecessary FFT round trip.
+        for (int is = 0; is < nspin; ++is)
         {
             for (int ibs = 0; ibs < this->pw_rho->nrxx; ++ibs)
             {
+                // The density from pot.init_pot -> Charge::atomic_rho may contain negative elements.
+                this->chr.rho[is][ibs] = this->nelec_[is] / ucell.omega;
                 this->pphi_[is][ibs] = sqrt(this->chr.rho[is][ibs]);
             }
         }
     }
 
-    for (int is = 0; is < PARAM.inp.nspin; ++is)
+    for (int is = 0; is < nspin; ++is)
     {
         this->pelec->eferm.set_efval(is, 0);
         this->theta_[is] = 0.;
@@ -267,7 +271,7 @@ void ESolver_OF::before_opt(const int istep, UnitCell& ucell)
         ModuleBase::GlobalFunc::ZEROS(this->pdEdphi_[is], this->pw_rho->nrxx);
         ModuleBase::GlobalFunc::ZEROS(this->pdirect_[is], this->pw_rho->nrxx);
     }
-    if (PARAM.inp.nspin == 1)
+    if (nspin == 1)
     {
         this->theta_[0] = 0.2;
     }