Fix 1. get_pchg with nspin=4 and 2. out_mul with DFT+U and 3. fix updates in develop branch

source/module_cell/module_neighbor/sltk_grid.cpp

@@ -152,12 +152,19 @@ void Grid::setMemberVariables(std::ofstream& ofs_in, //  output data to ofs
     this->clear_atoms();
 
     // random selection, in order to estimate again.
-    this->x_min = ucell.atoms[0].tau[0].x;
-    this->y_min = ucell.atoms[0].tau[0].y;
-    this->z_min = ucell.atoms[0].tau[0].z;
-    this->x_max = ucell.atoms[0].tau[0].x;
-    this->y_max = ucell.atoms[0].tau[0].y;
-    this->z_max = ucell.atoms[0].tau[0].z;
+    for (int it = 0; it < ucell.ntype; it++)
+    {
+        if (ucell.atoms[it].na > 0)
+        {
+            this->x_min = ucell.atoms[it].tau[0].x;
+            this->y_min = ucell.atoms[it].tau[0].y;
+            this->z_min = ucell.atoms[it].tau[0].z;
+            this->x_max = ucell.atoms[it].tau[0].x;
+            this->y_max = ucell.atoms[it].tau[0].y;
+            this->z_max = ucell.atoms[it].tau[0].z;
+            break;
+        }
+    }
 
     ModuleBase::Vector3<double> vec1(ucell.latvec.e11, ucell.latvec.e12, ucell.latvec.e13);
     ModuleBase::Vector3<double> vec2(ucell.latvec.e21, ucell.latvec.e22, ucell.latvec.e23);

source/module_cell/update_cell.cpp

@@ -384,41 +384,21 @@ void periodic_boundary_adjustment(Atom* atoms,
     // first adjust direct coordinates,
     // then update them into cartesian coordinates,
     //----------------------------------------------
-    for (int i=0;i<ntype;i++) {
-        Atom* atom = &atoms[i];
-        for (int j=0;j<atom->na;j++) {
-            printf("the taud is %f %f %f\n",atom->taud[j].x,atom->taud[j].y,atom->taud[j].z);
-        }
-    }
     for (int it = 0; it < ntype; it++) {
         Atom* atom = &atoms[it];
         for (int ia = 0; ia < atom->na; ia++) {
             // mohan update 2011-03-21
-            if (atom->taud[ia].x < 0)
+            for (int ik = 0; ik < 3; ik++) 
             {
-                atom->taud[ia].x += 1.0;
-            }
-            if (atom->taud[ia].y < 0)
-            {
-                atom->taud[ia].y += 1.0;
-            }
-            if (atom->taud[ia].z < 0)
-            {
-                atom->taud[ia].z += 1.0;
-            }
-            if (atom->taud[ia].x >= 1.0)
-            {
-                atom->taud[ia].x -= 1.0;
-            }
-            if (atom->taud[ia].y >= 1.0)
-            {
-                atom->taud[ia].y -= 1.0;
-            }
-            if (atom->taud[ia].z >= 1.0)
-            {
-                atom->taud[ia].z -= 1.0;
+                if (atom->taud[ia][ik] < 0) 
+                {
+                    atom->taud[ia][ik] += 1.0;
+                }
+                if (atom->taud[ia][ik] >= 1.0) 
+                {
+                    atom->taud[ia][ik] -= 1.0;
+                }
             }
-
             if (atom->taud[ia].x < 0 
                 || atom->taud[ia].y < 0
                 || atom->taud[ia].z < 0 

source/module_elecstate/module_charge/charge_init.cpp

@@ -37,6 +37,7 @@ void Charge::init_rho(elecstate::efermi& eferm_iout,
     this->pgrid = &pgrid;
 
     bool read_error = false;
+    bool read_kin_error = false;
     if (PARAM.inp.init_chg == "file" || PARAM.inp.init_chg == "auto")
     {
         GlobalV::ofs_running << " try to read charge from file" << std::endl;
@@ -101,72 +102,100 @@ void Charge::init_rho(elecstate::efermi& eferm_iout,
             }
         }
 
-        if (XC_Functional::get_func_type() == 3 || XC_Functional::get_func_type() == 5)
+        if (read_error)
         {
-            GlobalV::ofs_running << " try to read kinetic energy density from file" << std::endl;
-            // try to read charge from binary file first, which is the same as QE
-            std::vector<std::complex<double>> kin_g_space(PARAM.inp.nspin * this->ngmc, {0.0, 0.0});
-            std::vector<std::complex<double>*> kin_g;
-            for (int is = 0; is < PARAM.inp.nspin; is++)
+            const std::string warn_msg
+                = " WARNING: \"init_chg\" is enabled but ABACUS failed to read charge density from file.\n"
+                  " Please check if there is SPINX_CHG.cube (X=1,...) or {suffix}-CHARGE-DENSITY.restart in the "
+                  "directory.\n";
+            std::cout << std::endl << warn_msg;
+            if (PARAM.inp.init_chg == "file")
             {
-                kin_g.push_back(kin_g_space.data() + is * this->ngmc);
+                ModuleBase::WARNING_QUIT("Charge::init_rho",
+                                         "Failed to read in charge density from file.\nIf you want to use atomic "
+                                         "charge initialization, \nplease set init_chg to atomic in INPUT.");
             }
+        }
 
-            std::stringstream binary;
-            binary << PARAM.globalv.global_readin_dir << PARAM.inp.suffix + "-TAU-DENSITY.restart";
-            if (ModuleIO::read_rhog(binary.str(), rhopw, kin_g.data()))
+        if (XC_Functional::get_func_type() == 3 || XC_Functional::get_func_type() == 5)
+        {
+            // If the charge density is not read in, then the kinetic energy density is not read in either
+            if (!read_error)
             {
-                GlobalV::ofs_running << " Read in the kinetic energy density: " << binary.str() << std::endl;
-                for (int is = 0; is < PARAM.inp.nspin; ++is)
+                GlobalV::ofs_running << " try to read kinetic energy density from file" << std::endl;
+                // try to read charge from binary file first, which is the same as QE
+                std::vector<std::complex<double>> kin_g_space(PARAM.inp.nspin * this->ngmc, {0.0, 0.0});
+                std::vector<std::complex<double>*> kin_g;
+                for (int is = 0; is < PARAM.inp.nspin; is++)
                 {
-                    rhopw->recip2real(kin_g[is], this->kin_r[is]);
+                    kin_g.push_back(kin_g_space.data() + is * this->ngmc);
                 }
-            }
-            else
-            {
-                for (int is = 0; is < PARAM.inp.nspin; is++)
+
+                std::stringstream binary;
+                binary << PARAM.globalv.global_readin_dir << PARAM.inp.suffix + "-TAU-DENSITY.restart";
+                if (ModuleIO::read_rhog(binary.str(), rhopw, kin_g.data()))
                 {
-                    std::stringstream ssc;
-                    ssc << PARAM.globalv.global_readin_dir << "SPIN" << is + 1 << "_TAU.cube";
-                    // mohan update 2012-02-10, sunliang update 2023-03-09
-                    if (ModuleIO::read_vdata_palgrid(
-                            pgrid,
-                            (PARAM.inp.esolver_type == "sdft" ? GlobalV::RANK_IN_STOGROUP : GlobalV::MY_RANK),
-                            GlobalV::ofs_running,
-                            ssc.str(),
-                            this->kin_r[is],
-                            ucell.nat))
+                    GlobalV::ofs_running << " Read in the kinetic energy density: " << binary.str() << std::endl;
+                    for (int is = 0; is < PARAM.inp.nspin; ++is)
                     {
-                        GlobalV::ofs_running << " Read in the kinetic energy density: " << ssc.str() << std::endl;
+                        rhopw->recip2real(kin_g[is], this->kin_r[is]);
                     }
                 }
+                else
+                {
+                    for (int is = 0; is < PARAM.inp.nspin; is++)
+                    {
+                        std::stringstream ssc;
+                        ssc << PARAM.globalv.global_readin_dir << "SPIN" << is + 1 << "_TAU.cube";
+                        // mohan update 2012-02-10, sunliang update 2023-03-09
+                        if (ModuleIO::read_vdata_palgrid(
+                                pgrid,
+                                (PARAM.inp.esolver_type == "sdft" ? GlobalV::RANK_IN_STOGROUP : GlobalV::MY_RANK),
+                                GlobalV::ofs_running,
+                                ssc.str(),
+                                this->kin_r[is],
+                                ucell.nat))
+                        {
+                            GlobalV::ofs_running << " Read in the kinetic energy density: " << ssc.str() << std::endl;
+                        }
+                        else
+                        {
+                            read_kin_error = true;
+                            std::cout << " WARNING: \"init_chg\" is enabled but ABACUS failed to read kinetic energy "
+                                         "density from file.\n"
+                                         " Please check if there is SPINX_TAU.cube (X=1,...) or "
+                                         "{suffix}-TAU-DENSITY.restart in the directory.\n"
+                                      << std::endl;
+                            break;
+                        }
+                    }
+                }
+            }
+            else
+            {
+                read_kin_error = true;
             }
         }
     }
-    if (read_error)
+
+    if (PARAM.inp.init_chg == "atomic" || read_error) // mohan add 2007-10-17
     {
-        const std::string warn_msg = " WARNING: \"init_chg\" is enabled but ABACUS failed to read charge density from file.\n"
-                                     " Please check if there is SPINX_CHG.cube (X=1,...) or {suffix}-CHARGE-DENSITY.restart in the directory.\n";
-        std::cout << std::endl << warn_msg;
-        if (PARAM.inp.init_chg == "auto")
+        if (read_error)
         {
-            std::cout << " Charge::init_rho: use atomic initialization instead." << std::endl << std::endl;
-        }
-        else if (PARAM.inp.init_chg == "file")
-        {
-            ModuleBase::WARNING_QUIT("Charge::init_rho", "Failed to read in charge density from file.\nIf you want to use atomic charge initialization, \nplease set init_chg to atomic in INPUT.");
+            std::cout << " Charge::init_rho: use atomic initialization instead." << std::endl;
         }
+        this->atomic_rho(PARAM.inp.nspin, ucell.omega, rho, strucFac, ucell);
     }
 
-    if (PARAM.inp.init_chg == "atomic" || 
-        (PARAM.inp.init_chg == "auto" && read_error)) // mohan add 2007-10-17
+    // wenfei 2021-7-29 : initial tau = 3/5 rho^2/3, Thomas-Fermi
+    if (XC_Functional::get_func_type() == 3 || XC_Functional::get_func_type() == 5)
     {
-        this->atomic_rho(PARAM.inp.nspin, ucell.omega, rho, strucFac, ucell);
-
-        // liuyu 2023-06-29 : move here from atomic_rho(), which will be called several times in charge extrapolation
-        // wenfei 2021-7-29 : initial tau = 3/5 rho^2/3, Thomas-Fermi
-        if (XC_Functional::get_func_type() == 3 || XC_Functional::get_func_type() == 5)
+        if (PARAM.inp.init_chg == "atomic" || read_kin_error)
         {
+            if (read_kin_error)
+            {
+                std::cout << " Charge::init_rho: init kinetic energy density from rho." << std::endl;
+            }
             const double pi = 3.141592653589790;
             const double fact = (3.0 / 5.0) * pow(3.0 * pi * pi, 2.0 / 3.0);
             for (int is = 0; is < PARAM.inp.nspin; ++is)

source/module_esolver/esolver_gets.cpp

@@ -147,7 +147,7 @@ void ESolver_GetS::runner(UnitCell& ucell, const int istep)
 }
 
 void ESolver_GetS::after_all_runners(UnitCell& ucell) {};
-double ESolver_GetS::cal_energy() {};
+double ESolver_GetS::cal_energy() { return 0.0; };
 void ESolver_GetS::cal_force(UnitCell& ucell, ModuleBase::matrix& force) {};
 void ESolver_GetS::cal_stress(UnitCell& ucell, ModuleBase::matrix& stress) {};
 

source/module_esolver/esolver_ks.h

@@ -43,7 +43,7 @@ class ESolver_KS : public ESolver_FP
     virtual void iter_init(UnitCell& ucell, const int istep, const int iter);
 
     //! Something to do after hamilt2density function in each iter loop.
-    virtual void iter_finish(UnitCell& ucell, const int istep, int& iter);
+    virtual void iter_finish(UnitCell& ucell, const int istep, int& iter) override;
 
     // calculate electron density from a specific Hamiltonian with ethr
     virtual void hamilt2density_single(UnitCell& ucell, const int istep, const int iter, const double ethr);

source/module_hamilt_lcao/module_deepks/LCAO_deepks_interface.cpp

@@ -69,7 +69,7 @@ void LCAO_Deepks_Interface<TK, TR>::out_deepks_labels(const double& etot,
 
         if (PARAM.inp.deepks_bandgap)
         {
-            const int nocc = PARAM.inp.nelec / 2;
+            const int nocc = (PARAM.inp.nelec+1) / 2;
             std::vector<double> o_tot(nks);
             for (int iks = 0; iks < nks; ++iks)
             {

source/module_hsolver/diago_iter_assist.cpp

@@ -118,6 +118,9 @@ void DiagoIterAssist<T, Device>::diagH_subspace(const hamilt::Hamilt<T, Device>*
     // after generation of H and S matrix, diag them
     DiagoIterAssist::diagH_LAPACK(nstart, n_band, hcc, scc, nstart, en, vcc);
 
+
+    const int ld_temp = in_place ? dmax : dmin;
+
     { // code block to calculate evc
         gemm_op<T, Device>()(ctx,
                              'N',
@@ -132,12 +135,12 @@ void DiagoIterAssist<T, Device>::diagH_subspace(const hamilt::Hamilt<T, Device>*
                              nstart,
                              &zero,
                              temp,
-                             dmin);
+                             ld_temp);
     }
 
     if (!in_place)
     {
-        matrixSetToAnother<T, Device>()(ctx, n_band, temp, dmin, evc.get_pointer(), dmax);
+        matrixSetToAnother<T, Device>()(ctx, n_band, temp, ld_temp, evc.get_pointer(), dmax);
         delmem_complex_op()(ctx, temp);
     }
     delmem_complex_op()(ctx, hcc);

source/module_io/get_pchg_lcao.cpp

@@ -200,7 +200,9 @@ void IState_Charge::begin(Gint_k& gk,
         mode = 3;
     }
 
-    const int fermi_band = static_cast<int>((nelec + 1) / 2 + 1.0e-8);
+    int fermi_band = static_cast<int>((nelec + 1) / 2 + 1.0e-8);
+    // for nspin = 4 case, fermi_band is the number of electrons because of each band has 1 electron
+    if(PARAM.inp.nspin == 4) fermi_band = nelec;
     std::cout << " number of electrons = " << nelec << std::endl;
     std::cout << " number of occupied bands = " << fermi_band << std::endl;
 
@@ -216,15 +218,15 @@ void IState_Charge::begin(Gint_k& gk,
             elecstate::DensityMatrix<std::complex<double>, double> DM(this->ParaV, nspin_dm, kv.kvec_d, kv.get_nks() / nspin_dm);
 
 #ifdef __MPI
-            this->idmatrix(ib, nspin, nelec, nlocal, wg, DM, kv, if_separate_k);
+            this->idmatrix(ib, nspin_dm, nelec, nlocal, wg, DM, kv, if_separate_k);
 #else
             ModuleBase::WARNING_QUIT("IState_Charge::begin", "The `pchg` calculation is only available for MPI now!");
 #endif
             // If contribution from different k-points need to be output separately
             if (if_separate_k)
             {
                 // For multi-k, loop over all real k-points
-                for (int ik = 0; ik < kv.get_nks() / nspin; ++ik)
+                for (int ik = 0; ik < kv.get_nks() / nspin_dm; ++ik)
                 {
                     for (int is = 0; is < nspin; ++is)
                     {
@@ -506,16 +508,8 @@ void IState_Charge::idmatrix(const int& ib,
     ModuleBase::TITLE("IState_Charge", "idmatrix");
     assert(wg.nr == kv.get_nks());
 
-    const int fermi_band = static_cast<int>((nelec + 1) / 2 + 1.0e-8);
-
     // To ensure the normalization of charge density in multi-k calculation (if if_separate_k is true)
-    double wg_sum_k = 0;
-    double wg_sum_k_homo = 0;
-    for (int ik = 0; ik < kv.get_nks() / nspin; ++ik)
-    {
-        wg_sum_k += wg(ik, ib);
-        wg_sum_k_homo += wg(ik, fermi_band - 1);
-    }
+    double wg_sum_k = 1.0;
 
     for (int ik = 0; ik < kv.get_nks(); ++ik)
     {
@@ -530,11 +524,11 @@ void IState_Charge::idmatrix(const int& ib,
             double wg_value;
             if (if_separate_k)
             {
-                wg_value = (ib < fermi_band) ? wg_sum_k : wg_sum_k_homo;
+                wg_value = wg_sum_k;
             }
             else
             {
-                wg_value = (ib < fermi_band) ? wg(ik, ib) : wg(ik, fermi_band - 1);
+                wg_value = wg(ik, 0);
             }
             wg_local[ib_local] = wg_value;
         }

source/module_io/input_conv.h

@@ -12,8 +12,8 @@
 #include <iomanip>
 #include <iostream>
 #include <regex.h>
-#include <stdio.h>
-#include <string.h>
+#include <cstdio>
+#include <cstring>
 #include <string>
 #include <vector>
 #include <algorithm>
@@ -31,7 +31,7 @@ void tmp_convert();
  * @brief Pass the data members from the INPUT instance(defined in
  * module_io/input.cpp) to GlobalV and GlobalC.
  */
-void Convert(void);
+void Convert();
 
 /**
  * @brief To parse input parameters as expressions into vectors
@@ -47,7 +47,7 @@ void parse_expression(const std::string& fn, std::vector<T>& vec)
 {
     ModuleBase::TITLE("Input_Conv", "parse_expression");
     int count = 0;
-    std::string pattern("([0-9]+\\*[0-9.]+|[0-9,.]+)");
+    std::string pattern("([-+]?[0-9]+\\*[-+]?[0-9.]+|[-+]?[0-9,.]+)");
     std::vector<std::string> str;
     std::stringstream ss(fn);
     std::string section;
@@ -103,6 +103,7 @@ void parse_expression(const std::string& fn, std::vector<T>& vec)
         {
             size_t pos = sub_str.find("*");
             int num = stoi(sub_str.substr(0, pos));
+            assert(num>=0);
             T occ = stof(sub_str.substr(pos + 1, sub_str.size()));
             // std::vector<double> ocp_temp(num, occ);
             // const std::vector<double>::iterator dest = vec.begin() + count;

source/module_io/output_mulliken.h

@@ -138,7 +138,7 @@ void cal_mag(Parallel_Orbitals* pv,
             auto sc_lambda
                 = new hamilt::DeltaSpin<hamilt::OperatorLCAO<TK, double>>(nullptr,
                                                                           kv.kvec_d,
-                                                                          nullptr,
+                                                                          dynamic_cast<hamilt::HamiltLCAO<TK, double>*>(p_ham)->getHR(),
                                                                           ucell,
                                                                           &gd,
                                                                           two_center_bundle.overlap_orb_onsite.get(),
@@ -164,7 +164,7 @@ void cal_mag(Parallel_Orbitals* pv,
             auto sc_lambda = new hamilt::DeltaSpin<hamilt::OperatorLCAO<std::complex<double>, std::complex<double>>>(
                 nullptr,
                 kv.kvec_d,
-                nullptr,
+                dynamic_cast<hamilt::HamiltLCAO<std::complex<double>, std::complex<double>>*>(p_ham)->getHR(),
                 ucell,
                 &gd,
                 two_center_bundle.overlap_orb_onsite.get(),