Fix: DFT+U nscf calculation of nspin=1 and output onsite.dm with out_chg>=0 (#6550)

* Fix: DFT+U nscf calculation of nspin=1 and output onsite.dm with out_chg>=0

* Fix: nscf force calculation with DFT+U

---------

Co-authored-by: dyzheng <[email protected]>

Author: dyzheng

source/module_esolver/esolver_ks_lcao.cpp

@@ -669,7 +669,10 @@ void ESolver_KS_LCAO<TK, TR>::iter_init(UnitCell& ucell, const int istep, const
 
     if (PARAM.inp.dft_plus_u)
     {
-        if (istep != 0 || iter != 1)
+        // set_dmr to calculate the onsite-matrix and the energy correction for DFT+U,
+        // in SCF calculation, the DMR is updated from the second iteration
+        // in NSCF calculation, the DMR is not updated, so always set_dmr here
+        if (istep != 0 || iter != 1 || PARAM.inp.calculation == "nscf")
         {
             GlobalC::dftu.set_dmr(dynamic_cast<elecstate::ElecStateLCAO<TK>*>(this->pelec)->get_DM());
         }

source/module_hamilt_lcao/module_dftu/dftu_io.cpp

@@ -55,7 +55,7 @@ void DFTU::output(const UnitCell &ucell)
 
     //Write onsite.dm
     std::ofstream ofdftu;
-    if(PARAM.inp.out_chg[0]){
+    if(PARAM.inp.out_chg[0] != -1){
       if(GlobalV::MY_RANK == 0){
         ofdftu.open(PARAM.globalv.global_out_dir + "onsite.dm");
       }
@@ -309,7 +309,7 @@ void DFTU::read_occup_m(const UnitCell& ucell,
 
                             if (PARAM.inp.nspin == 1 || PARAM.inp.nspin == 2)
                             {
-                                for (int is = 0; is < 2; is++)
+                                for (int is = 0; is < PARAM.inp.nspin; is++)
                                 {
                                     ifdftu >> word;
                                     if (strcmp("spin", word) == 0)