Fix : dftu force and stress (#1868)

* Fix : force matrix not allocated but used

* fix : calling finish_k after dftu force & stress

---------

Co-authored-by: wenfei-li <[email protected]>

Author: wenfei-li

source/module_hamilt_lcao/hamilt_lcaodft/FORCE_STRESS.cpp

@@ -196,6 +196,8 @@ void Force_Stress_LCAO::getForceStress(
 		}
 		GlobalC::dftu.force_stress(loc.dm_gamma, loc.dm_k, *uhm.LM, force_dftu, stress_dftu);
 	}
+
+	if(!GlobalV::GAMMA_ONLY_LOCAL) this->flk.finish_k();
 #ifdef __EXX
 	//Force and Stress contribution from exx
 	ModuleBase::matrix force_exx;
@@ -583,7 +585,7 @@ void Force_Stress_LCAO::getForceStress(
 			}
 			if(GlobalV::dft_plus_u)
 			{
-				sc_pw.print_stress("DFTU     STRESS",sigmaxc,GlobalV::TEST_STRESS,ry);
+				sc_pw.print_stress("DFTU     STRESS",stress_dftu,GlobalV::TEST_STRESS,ry);
 			}
 			sc_pw.print_stress("TOTAL    STRESS",scs,GlobalV::TEST_STRESS,ry);
 

source/module_hamilt_lcao/hamilt_lcaodft/FORCE_k.cpp

@@ -177,18 +177,6 @@ void Force_LCAO_k::ftable_k(const bool isforce,
         Parallel_Reduce::reduce_double_pool(svl_dphi.c, svl_dphi.nr * svl_dphi.nc);
     }
 
-    // test the force.
-    /*
-    std::cout << " overlap force" << std::endl;
-    for(int iat=0; iat<GlobalC::ucell.nat; ++iat)
-    {
-        const double fac = ModuleBase::Ry_to_eV / 0.529177;
-        std::cout << std::setw(5) << iat+1 << std::setw(15) << foverlap[iat][0] *fac<< std::setw(15) <<
-    foverlap[iat][1]*fac << std::setw(15) << foverlap[iat][2]*fac << std::endl;
-    }
-    */
-
-    this->finish_k();
 
     ModuleBase::timer::tick("Force_LCAO_k", "ftable_k");
     return;

source/module_hamilt_lcao/module_dftu/dftu_folding.cpp

@@ -189,15 +189,13 @@ void DFTU::folding_matrix_k(const int ik, const int dim1, const int dim2, std::c
                         // the index of orbitals in this processor
                         const int iw1_all = start1 + ii;
                         const int mu = this->LM->ParaV->trace_loc_row[iw1_all];
-                        if (mu < 0)
-                            continue;
+                        if (mu < 0) continue;
 
                         for (int jj = 0; jj < atom2->nw * GlobalV::NPOL; jj++)
                         {
                             int iw2_all = start2 + jj;
                             const int nu = this->LM->ParaV->trace_loc_col[iw2_all];
-                            if (nu < 0)
-                                continue;
+                            if (nu < 0) continue;
 
                             int iic;
                             if (ModuleBase::GlobalFunc::IS_COLUMN_MAJOR_KS_SOLVER())