Fix: force error in DFT+U

Author: dyzheng

source/module_esolver/esolver_ks_pw.cpp

@@ -393,6 +393,13 @@ void ESolver_KS_PW<T, Device>::before_scf(const int istep) {
     if (GlobalC::ucell.cell_parameter_updated) {
         this->init_after_vc(PARAM.inp, GlobalC::ucell);
     }
+    
+    // if uramping > 0, and someone want to test the pbe first in each relax step, use this line below:
+    //ModuleBase::GlobalFunc::ZEROS(GlobalC::dftu.U, GlobalC::ucell.ntype);
+
+    // if someone want to test relax from atomic charge but not extrapolated charge, use these two lines below rather than the if(){} block:
+    //this->CE.update_all_dis(GlobalC::ucell);
+    //this->sf.setup_structure_factor(&GlobalC::ucell, this->pelec->charge->rhopw);
     if (GlobalC::ucell.ionic_position_updated) {
         this->CE.update_all_dis(GlobalC::ucell);
         this->CE.extrapolate_charge(
@@ -601,10 +608,6 @@ void ESolver_KS_PW<T, Device>::iter_init(const int istep, const int iter) {
         this->p_chgmix->mixing_restart_count++;
         if (GlobalV::dft_plus_u)
         {
-            if (GlobalC::dftu.uramping > 0.01 && !GlobalC::dftu.u_converged())
-            {
-                this->p_chgmix->mixing_restart_step = GlobalV::SCF_NMAX + 1;
-            }
             if (GlobalC::dftu.uramping > 0.01)
             {
                 bool do_uramping = true;
@@ -627,6 +630,10 @@ void ESolver_KS_PW<T, Device>::iter_init(const int istep, const int iter) {
                     std::cout << " eV " << std::endl;
                 }
             }
+            if (GlobalC::dftu.uramping > 0.01 && !GlobalC::dftu.u_converged())
+            {
+                this->p_chgmix->mixing_restart_step = GlobalV::SCF_NMAX + 1;
+            }
         }
     }
     // mohan move harris functional to here, 2012-06-05

source/module_hamilt_pw/hamilt_pwdft/kernels/force_op.cpp

@@ -274,6 +274,7 @@ struct cal_force_nl_op<FPTYPE, base_device::DEVICE_CPU>
             const int nproj = atom_nh[it];
             if(orbital_l == -1)
             {
+                iat0 += atom_na[it];
                 sum0 += nproj * atom_na[it];
                 continue;
             }

source/module_hamilt_pw/hamilt_pwdft/kernels/stress_op.cpp

@@ -227,7 +227,7 @@ struct cal_stress_nl_op<FPTYPE, base_device::DEVICE_CPU>
                     FPTYPE* stress)
     {
         FPTYPE local_stress = 0;
-        int iat = 0, sum = 0;
+        int sum = 0;
         for (int it = 0; it < ntype; it++)
         {
             const int orbital_l = orbital_corr[it];
@@ -294,7 +294,6 @@ struct cal_stress_nl_op<FPTYPE, base_device::DEVICE_CPU>
                 vu += npol * npol * tlp1_2;// step for vu
             }// ia
             sum += atom_na[it] * nproj;
-            iat += atom_na[it];
         } // end it
         *stress += local_stress;
     };

source/module_hamilt_pw/hamilt_pwdft/onsite_projector.cpp

@@ -108,23 +108,25 @@ void projectors::OnsiteProjector<T, Device>::init(const std::string& orbital_dir
                                                   const ModuleBase::matrix& ekb)
 {
     this->device = base_device::get_device_type<Device>(this->ctx);
-    if(!this->initialed)
-    {
-        this->ucell = ucell_in;
-        this->ntype = ucell_in->ntype;
+    
+    this->ucell = ucell_in;
+    this->ntype = ucell_in->ntype;
 
-        this->pw_basis_ = &pw_basis;
-        this->sf_ = &sf;
+    this->pw_basis_ = &pw_basis;
+    this->sf_ = &sf;
 
-        std::vector<std::string> orb_files(ntype);
-        std::vector<int> nproj(ntype);
-        int sum_nproj = 0;
-        for(int it=0;it<ntype;++it)
-        {
-            orb_files[it] = ucell->orbital_fn[it];
-            nproj[it] = ucell->atoms[it].nwl;
-            sum_nproj += nproj[it];
-        }
+    std::vector<std::string> orb_files(ntype);
+    std::vector<int> nproj(ntype);
+    int sum_nproj = 0;
+    for(int it=0;it<ntype;++it)
+    {
+        orb_files[it] = ucell->orbital_fn[it];
+        nproj[it] = ucell->atoms[it].nwl;
+        sum_nproj += nproj[it];
+    }
+    // do orbital modulate only once
+    if(!this->initialed)
+    {
         this->lproj.resize(sum_nproj);
         int index = 0;
         for(int it=0;it<ntype;++it)
@@ -157,7 +159,7 @@ void projectors::OnsiteProjector<T, Device>::init(const std::string& orbital_dir
                         iproj, 
                         onsite_r);
 
-        ModuleBase::timer::tick("OnsiteProj", "cubspl_tabulate");
+        
         // STAGE 0 - making the interpolation table
         // CACHE 0 - if cache the irow2it, irow2iproj, irow2m, itiaiprojm2irow, <G+k|p> can be reused for 
         //           SCF, RELAX and CELL-RELAX calculation
@@ -166,24 +168,26 @@ void projectors::OnsiteProjector<T, Device>::init(const std::string& orbital_dir
         RadialProjection::RadialProjector::_build_forward_map(it2ia, it2iproj, lproj, itiaiprojm2irow_);
         //rp_._build_sbt_tab(rgrid, projs, lproj, nq, dq);
         rp_._build_sbt_tab(nproj, rgrid, projs, lproj, nq, dq, ucell_in->omega, psi.npol, tab, nhtol);
-        // For being compatible with present cal_force and cal_stress framework  
-        // uncomment the following code block if you want to use the FS_Nonlocal_tools
-        if(this->tab_atomic_ == nullptr)
-        {
-            this->tot_nproj = itiaiprojm2irow_.size();
-            this->npwx_ = this->pw_basis_->npwk_max;
-            this->size_vproj = this->tot_nproj * this->npwx_;
-            resmem_complex_op()(this->ctx, this->tab_atomic_, this->size_vproj, "OnsiteP::tab_atomic_");
-        }
-
-        delete this->fs_tools; // it is okay to delete nullptr
-        this->fs_tools = new hamilt::FS_Nonlocal_tools<T, Device>(
-            nproj, lproj, tab, nhtol, this->tab_atomic_, ucell_in, &psi, &kv, &pw_basis, &sf, wg, ekb);      
 
         ModuleBase::timer::tick("OnsiteProj", "cubspl_tabulate");
 
         this->initialed = true;
     }
+
+    // For being compatible with present cal_force and cal_stress framework  
+    // uncomment the following code block if you want to use the FS_Nonlocal_tools
+    ModuleBase::timer::tick("OnsiteProj", "cubspl_tabulate");
+    {
+        this->tot_nproj = itiaiprojm2irow_.size();
+        this->npwx_ = this->pw_basis_->npwk_max;
+        this->size_vproj = this->tot_nproj * this->npwx_;
+        resmem_complex_op()(this->ctx, this->tab_atomic_, this->size_vproj, "OnsiteP::tab_atomic_");
+    }
+
+    delete this->fs_tools; // it is okay to delete nullptr
+    this->fs_tools = new hamilt::FS_Nonlocal_tools<T, Device>(
+        nproj, lproj, tab, nhtol, this->tab_atomic_, ucell_in, &psi, &kv, &pw_basis, &sf, wg, ekb);      
+    ModuleBase::timer::tick("OnsiteProj", "cubspl_tabulate");
 }
 
 template<typename T, typename Device>