Merge branch 'develop' into fft5

Author: A-006

docs/advanced/input_files/input-main.md

@@ -989,7 +989,7 @@ calculations.
 
 - **Type**: String
 - **Description**: In our package, the XC functional can either be set explicitly using the `dft_functional` keyword in `INPUT` file. If `dft_functional` is not specified, ABACUS will use the xc functional indicated in the pseudopotential file.
-  On the other hand, if dft_functional is specified, it will overwrite the functional from pseudopotentials and performs calculation with whichever functional the user prefers. We further offer two ways of supplying exchange-correlation functional. The first is using 'short-hand' names such as 'LDA', 'PBE', 'SCAN'. A complete list of 'short-hand' expressions can be found in [the source code](../../../source/module_hamilt_general/module_xc/xc_functional.cpp). The other way is only available when ***compiling with LIBXC***, and it allows for supplying exchange-correlation functionals as combinations of LIBXC keywords for functional components, joined by a plus sign, for example, 'dft_functional='LDA_X_1D_EXPONENTIAL+LDA_C_1D_CSC'. The list of LIBXC keywords can be found on its [website](https://www.tddft.org/programs/libxc/functionals/). In this way, **we support all the LDA,GGA and mGGA functionals provided by LIBXC**.
+  On the other hand, if dft_functional is specified, it will overwrite the functional from pseudopotentials and performs calculation with whichever functional the user prefers. We further offer two ways of supplying exchange-correlation functional. The first is using 'short-hand' names such as 'LDA', 'PBE', 'SCAN'. A complete list of 'short-hand' expressions can be found in [the source code](../../../source/module_hamilt_general/module_xc/xc_functional.cpp). The other way is only available when ***compiling with LIBXC***, and it allows for supplying exchange-correlation functionals as combinations of LIBXC keywords for functional components, joined by a plus sign, for example, dft_functional='LDA_X_1D_EXPONENTIAL+LDA_C_1D_CSC'. The list of LIBXC keywords can be found on its [website](https://libxc.gitlab.io/functionals/). In this way, **we support all the LDA,GGA and mGGA functionals provided by LIBXC**.
 
   Furthermore, the old INPUT parameter exx_hybrid_type for hybrid functionals has been absorbed into dft_functional. Options are `hf` (pure Hartree-Fock), `pbe0`(PBE0), `hse` (Note: in order to use HSE functional, LIBXC is required). Note also that HSE has been tested while PBE0 has NOT been fully tested yet, and the maximum CPU cores for running exx in parallel is $N(N+1)/2$, with N being the number of atoms. And forces for hybrid functionals are not supported yet.
 
@@ -1760,14 +1760,14 @@ The band (KS orbital) energy for each (k-point, spin, band) will be printed in t
 
 - **Type**: Boolean
 - **Availability**: Numerical atomic orbital basis
-- **Description**: Whether to print Hamiltonian matrices H(R)/density matrics DM(R) in npz format. This feature does not work for gamma-only calculations. Currently only intended for internal usage.
+- **Description**: Whether to print Hamiltonian matrices $H(R)$/density matrics $DM(R)$ in npz format. This feature does not work for gamma-only calculations. Currently only intended for internal usage.
 - **Default**: False
 
 ### dm_to_rho
 
 - **Type**: Boolean
 - **Availability**: Numerical atomic orbital basis
-- **Description**: Reads density matrix DM(R) in npz format and creates electron density on grids. This feature does not work for gamma-only calculations. Only supports serial calculations. Currently only intended for internal usage.
+- **Description**: Reads density matrix $DM(R)$ in npz format and creates electron density on grids. This feature does not work for gamma-only calculations. Only supports serial calculations. Currently only intended for internal usage.
 - **Default**: False
 
 ### out_app_flag

source/driver.cpp

@@ -183,7 +183,7 @@ void Driver::atomic_world()
     //--------------------------------------------------
 
     // where the actual stuff is done
-    this->driver_run();
+    this->driver_run(GlobalC::ucell);
 
     ModuleBase::timer::finish(GlobalV::ofs_running);
     ModuleBase::Memory::print_all(GlobalV::ofs_running);

source/driver.h

@@ -1,6 +1,8 @@
 #ifndef DRIVER_H
 #define DRIVER_H
 
+#include "module_cell/unitcell.h"
+
 class Driver
 {
   public:
@@ -34,7 +36,7 @@ class Driver
     void atomic_world();
 
     // the actual calculations
-    void driver_run();
+    void driver_run(UnitCell& ucell);
 };
 
 #endif

source/driver_run.cpp

@@ -24,7 +24,8 @@
  * the configuration-changing subroutine takes force and stress and updates the
  * configuration
  */
-void Driver::driver_run() {
+void Driver::driver_run(UnitCell& ucell)
+{
     ModuleBase::TITLE("Driver", "driver_line");
     ModuleBase::timer::tick("Driver", "driver_line");
 
@@ -39,37 +40,35 @@ void Driver::driver_run() {
 #endif
 
     // the life of ucell should begin here, mohan 2024-05-12
-    // delete ucell as a GlobalC in near future
-    GlobalC::ucell.setup_cell(PARAM.globalv.global_in_stru, GlobalV::ofs_running);
-    Check_Atomic_Stru::check_atomic_stru(GlobalC::ucell,
-                                         PARAM.inp.min_dist_coef);
+    ucell.setup_cell(PARAM.globalv.global_in_stru, GlobalV::ofs_running);
+    Check_Atomic_Stru::check_atomic_stru(ucell, PARAM.inp.min_dist_coef);
 
     //! 2: initialize the ESolver (depends on a set-up ucell after `setup_cell`)
-    ModuleESolver::ESolver* p_esolver = ModuleESolver::init_esolver(PARAM.inp, GlobalC::ucell);
+    ModuleESolver::ESolver* p_esolver = ModuleESolver::init_esolver(PARAM.inp, ucell);
 
     //! 3: initialize Esolver and fill json-structure
-    p_esolver->before_all_runners(PARAM.inp, GlobalC::ucell);
+    p_esolver->before_all_runners(ucell, PARAM.inp);
 
     // this Json part should be moved to before_all_runners, mohan 2024-05-12
 #ifdef __RAPIDJSON
-    Json::gen_stru_wrapper(&GlobalC::ucell);
+    Json::gen_stru_wrapper(&ucell);
 #endif
 
     const std::string cal_type = PARAM.inp.calculation;
 
     //! 4: different types of calculations
     if (cal_type == "md")
     {
-        Run_MD::md_line(GlobalC::ucell, p_esolver, PARAM);
+        Run_MD::md_line(ucell, p_esolver, PARAM);
     }
     else if (cal_type == "scf" || cal_type == "relax" || cal_type == "cell-relax" || cal_type == "nscf")
     {
         Relax_Driver rl_driver;
-        rl_driver.relax_driver(p_esolver);
+        rl_driver.relax_driver(p_esolver, ucell);
     }
     else if (cal_type == "get_S")
     {
-        p_esolver->runner(0, GlobalC::ucell);
+        p_esolver->runner(ucell, 0);
     }
     else
     {
@@ -79,11 +78,11 @@ void Driver::driver_run() {
         //! test_neighbour(LCAO),
         //! gen_bessel(PW), et al.
         const int istep = 0;
-        p_esolver->others(istep);
+        p_esolver->others(ucell, istep);
     }
 
     //! 5: clean up esolver
-    p_esolver->after_all_runners();
+    p_esolver->after_all_runners(ucell);
 
     ModuleESolver::clean_esolver(p_esolver);
 

source/module_base/opt_TN.hpp

@@ -1,9 +1,9 @@
 #ifndef OPT_TN_H
 #define OPT_TN_H
 
-#include <limits>
+#include "opt_CG.h"
 
-#include "./opt_CG.h"
+#include <limits>
 
 namespace ModuleBase
 {
@@ -25,7 +25,7 @@ class Opt_TN
     {
         this->mach_prec_ = std::numeric_limits<double>::epsilon(); // get machine precise
     }
-    ~Opt_TN(){};
+    ~Opt_TN() {};
 
     /**
      * @brief Allocate the space for the arrays in cg_.
@@ -54,7 +54,9 @@ class Opt_TN
     {
         this->iter_ = 0;
         if (nx_new != 0)
+        {
             this->nx_ = nx_new;
+        }
         this->cg_.refresh(nx_new);
     }
 
@@ -167,17 +169,23 @@ void Opt_TN::next_direct(double* px,
         epsilon = this->get_epsilon(px, cg_direct);
         // epsilon = 1e-9;
         for (int i = 0; i < this->nx_; ++i)
+        {
             temp_x[i] = px[i] + epsilon * cg_direct[i];
+        }
         (t->*p_calGradient)(temp_x, temp_gradient);
         for (int i = 0; i < this->nx_; ++i)
+        {
             temp_Hcgd[i] = (temp_gradient[i] - pgradient[i]) / epsilon;
+        }
 
         // get CG step length and update rdirect
         cg_alpha = cg_.step_length(temp_Hcgd, cg_direct, cg_ifPD);
         if (cg_ifPD == -1) // Hessian is not positive definite, and cgiter = 1.
         {
             for (int i = 0; i < this->nx_; ++i)
+            {
                 rdirect[i] += cg_alpha * cg_direct[i];
+            }
             flag = -1;
             break;
         }
@@ -188,14 +196,18 @@ void Opt_TN::next_direct(double* px,
         }
 
         for (int i = 0; i < this->nx_; ++i)
+        {
             rdirect[i] += cg_alpha * cg_direct[i];
+        }
 
         // store residuals used in truncated conditions
         last_residual = curr_residual;
         curr_residual = cg_.get_residual();
         cg_iter = cg_.get_iter();
         if (cg_iter == 1)
+        {
             init_residual = curr_residual;
+        }
 
         // check truncated conditions
         // if (curr_residual < 1e-12)

source/module_cell/unitcell.cpp

@@ -551,6 +551,9 @@ void UnitCell::setup_cell(const std::string& fn, std::ofstream& log) {
     this->atoms = new Atom[this->ntype]; // atom species.
     this->set_atom_flag = true;
 
+    this->symm.epsilon = PARAM.inp.symmetry_prec;
+    this->symm.epsilon_input = PARAM.inp.symmetry_prec;
+
     bool ok = true;
     bool ok2 = true;
 

source/module_elecstate/module_charge/charge_init.cpp

@@ -157,7 +157,25 @@ void Charge::init_rho(elecstate::efermi& eferm_iout,
     {
         for (int is = 0; is < PARAM.inp.nspin; ++is)
         {
-            GlobalC::restart.load_disk("charge", is, this->nrxx, rho[is]);
+            try
+            {
+                GlobalC::restart.load_disk("charge", is, this->nrxx, rho[is]);
+            }
+            catch (const std::exception& e)
+            {
+                // try to load from the output of `out_chg` 
+                std::stringstream ssc;
+                ssc << PARAM.globalv.global_readin_dir << "SPIN" << is + 1 << "_CHG.cube";
+                if (ModuleIO::read_vdata_palgrid(GlobalC::Pgrid,
+                    (PARAM.inp.esolver_type == "sdft" ? GlobalV::RANK_IN_STOGROUP : GlobalV::MY_RANK),
+                    GlobalV::ofs_running,
+                    ssc.str(),
+                    this->rho[is],
+                    GlobalC::ucell.nat))
+                {
+                    GlobalV::ofs_running << " Read in the charge density: " << ssc.str() << std::endl;
+                }
+            }
         }
         GlobalC::restart.info_load.load_charge_finish = true;
     }

source/module_esolver/esolver.cpp

@@ -247,8 +247,8 @@ ESolver* init_esolver(const Input_para& inp, UnitCell& ucell)
         {
             p_esolver = new ESolver_KS_LCAO<std::complex<double>, std::complex<double>>();
         }
-        p_esolver->before_all_runners(inp, ucell);
-        p_esolver->runner(0, ucell); // scf-only
+        p_esolver->before_all_runners(ucell, inp);
+        p_esolver->runner(ucell, 0); // scf-only
         // force and stress is not needed currently,
         // they will be supported after the analytical gradient
         // of LR-TDDFT is implemented.

source/module_esolver/esolver.h

@@ -20,26 +20,26 @@ class ESolver
     }
 
     //! initialize the energy solver by using input parameters and cell modules
-    virtual void before_all_runners(const Input_para& inp, UnitCell& cell) = 0;
+    virtual void before_all_runners(UnitCell& ucell, const Input_para& inp) = 0;
 
     //! run energy solver
-    virtual void runner(const int istep, UnitCell& cell) = 0;
+    virtual void runner(UnitCell& cell, const int istep) = 0;
 
     //! perform post processing calculations
-    virtual void after_all_runners(){};
+    virtual void after_all_runners(UnitCell& ucell){};
 
     //! deal with exx and other calculation than scf/md/relax/cell-relax:
     //! such as nscf, get_wf and get_pchg
-    virtual void others(const int istep){};
+    virtual void others(UnitCell& ucell, const int istep) {};
 
     //! calculate total energy of a given system
     virtual double cal_energy() = 0;
 
     //! calcualte forces for the atoms in the given cell
-    virtual void cal_force(ModuleBase::matrix& force) = 0;
+    virtual void cal_force(UnitCell& ucell, ModuleBase::matrix& force) = 0;
 
     //! calcualte stress of given cell
-    virtual void cal_stress(ModuleBase::matrix& stress) = 0;
+    virtual void cal_stress(UnitCell& ucell, ModuleBase::matrix& stress) = 0;
 
     bool conv_esolver = true; // whether esolver is converged
 

source/module_esolver/esolver_dp.cpp

@@ -31,7 +31,7 @@
 namespace ModuleESolver
 {
 
-void ESolver_DP::before_all_runners(const Input_para& inp, UnitCell& ucell)
+void ESolver_DP::before_all_runners(UnitCell& ucell, const Input_para& inp)
 {
     ucell_ = &ucell;
     dp_potential = 0;
@@ -57,7 +57,7 @@ void ESolver_DP::before_all_runners(const Input_para& inp, UnitCell& ucell)
 #endif
 }
 
-void ESolver_DP::runner(const int istep, UnitCell& ucell)
+void ESolver_DP::runner(UnitCell& ucell, const int istep)
 {
     ModuleBase::TITLE("ESolver_DP", "runner");
     ModuleBase::timer::tick("ESolver_DP", "runner");
@@ -127,13 +127,13 @@ double ESolver_DP::cal_energy()
     return dp_potential;
 }
 
-void ESolver_DP::cal_force(ModuleBase::matrix& force)
+void ESolver_DP::cal_force(UnitCell& ucell, ModuleBase::matrix& force)
 {
     force = dp_force;
     ModuleIO::print_force(GlobalV::ofs_running, *ucell_, "TOTAL-FORCE (eV/Angstrom)", force, false);
 }
 
-void ESolver_DP::cal_stress(ModuleBase::matrix& stress)
+void ESolver_DP::cal_stress(UnitCell& ucell, ModuleBase::matrix& stress)
 {
     stress = dp_virial;
 
@@ -148,7 +148,7 @@ void ESolver_DP::cal_stress(ModuleBase::matrix& stress)
     ModuleIO::print_stress("TOTAL-STRESS", stress, true, false);
 }
 
-void ESolver_DP::after_all_runners()
+void ESolver_DP::after_all_runners(UnitCell& ucell)
 {
     GlobalV::ofs_running << "\n\n --------------------------------------------" << std::endl;
     GlobalV::ofs_running << std::setprecision(16);