Merge branch 'develop' of https://github.com/deepmodeling/abacus-develop

 into stop

source/Makefile.Objects

@@ -460,9 +460,7 @@ OBJS_XC=xc_functional.o\
     xc_funct_corr_gga.o\
     xc_funct_hcth.o\
 
-OBJS_IO=input.o\
-    input_conv_tmp.o\
-    input_conv.o\
+OBJS_IO=input_conv.o\
     berryphase.o\
     bessel_basis.o\
     cal_test.o\
@@ -628,6 +626,7 @@ OBJS_SRCPW=H_Ewald_pw.o\
     forces_us.o\
     forces_nl.o\
     forces_cc.o\
+    forces_scc.o\
     fs_nonlocal_tools.o\
     force_op.o\
     stress_op.o\

source/driver.cpp

@@ -6,7 +6,6 @@
 #include "module_esolver/esolver.h"
 #include "module_hamilt_pw/hamilt_pwdft/global.h"
 #include "module_io/cal_test.h"
-#include "module_io/input.h"
 #include "module_io/input_conv.h"
 #include "module_io/para_json.h"
 #include "module_io/print_info.h"
@@ -132,7 +131,6 @@ void Driver::reading()
     read_input.write_parameters(PARAM, ss1.str());
 
     // (*temp*) copy the variables from INPUT to each class
-    Input_Conv::tmp_convert();
     Input_Conv::Convert();
 
     // (4) define the 'DIAGONALIZATION' world in MPI

source/driver_run.cpp

@@ -2,7 +2,6 @@
 #include "module_cell/check_atomic_stru.h"
 #include "module_cell/module_neighbor/sltk_atom_arrange.h"
 #include "module_hamilt_pw/hamilt_pwdft/global.h"
-#include "module_io/input.h"
 #include "module_parameter/parameter.h"
 #include "module_io/para_json.h"
 #include "module_io/print_info.h"

source/module_cell/module_symmetry/symmetry.cpp

@@ -1807,54 +1807,13 @@ void Symmetry::set_atom_map(const Atom* atoms)
                     if (equal(diff1, 0.0) && equal(diff2, 0.0) && equal(diff3, 0.0))
                     {
                         this->isym_rotiat_[k][ia] = ja;
+
                         break;
                     }
                 }
             }
         }
     }
-    this->reset_atom_map(atoms);
-}
-
-inline bool check_full_atom_map(const std::vector<std::vector<int>>& m)
-{
-    if (m.empty()) { return false;
-}
-    for (const auto& v1 : m) { for (const auto& v2 : v1) { if (v2 == -1) { return false; } } }
-    return true;
-}
-
-void Symmetry::reset_atom_map(const Atom* atoms)
-{
-    const double eps_start = this->epsilon;
-    while (!check_full_atom_map(this->isym_rotiat_) && this->epsilon < 100 * eps_start)
-    {
-        this->epsilon *= 2;
-        this->isym_rotiat_.clear();
-        this->set_atom_map(atoms);
-
-        std::cout << "try eps=" << this->epsilon << std::endl;
-        std::cout << check_full_atom_map(this->isym_rotiat_) << std::endl;
-    }
-    if (this->epsilon > 100 * eps_start)
-    {
-        ModuleBase::WARNING("Symmetry::reset_atom_map", "Cannot find all the equivalent atoms even at 100*symmetry_prec.");
-        if (ModuleSymmetry::Symmetry::symm_autoclose)
-        {
-            ModuleBase::WARNING("K_Vectors::ibz_kpoint", "Automatically set symmetry to 0 and continue ...");
-            std::cout << "Automatically set symmetry to 0 and continue ..." << std::endl;
-            ModuleSymmetry::Symmetry::symm_flag = 0;
-        }
-        else {
-            ModuleBase::WARNING_QUIT("Symmetry::reset_atom_map",
-                "Possible solutions: \n \
-1. Refine the lattice parameters in STRU;\n \
-2. Use a different`symmetry_prec`.  \n \
-3. Close symemtry: set `symmetry` to 0 in INPUT. \n \
-4. Set `symmetry_autoclose` to 1 in INPUT to automatically close symmetry when this error occurs.");
-        }
-    }
-    this->epsilon = eps_start;  //recover symmetry_prec after a successful reset
 }
 
 /// @brief return a map that is inequivalent atom index to its symmetry multiplicity

source/module_cell/module_symmetry/symmetry.h

@@ -130,8 +130,6 @@ class Symmetry : public Symmetry_Basic
 
     /// @brief  set atom map for each symmetry operation
     void set_atom_map(const Atom* atoms);
-    /// @brief deal with the rarely happen incomplete atom map due to a subtle symmetry_prec
-    void reset_atom_map(const Atom* atoms);
     /// @brief check if all the atoms are movable
     ///  delta_pos symmetrization in relax is only meaningful when all the atoms are movable in all the directions.
     bool is_all_movable(const Atom* atoms, const Statistics& st)const;

source/module_elecstate/potentials/H_TDDFT_pw.cpp

@@ -4,7 +4,6 @@
 #include "module_base/timer.h"
 #include "module_hamilt_lcao/module_tddft/evolve_elec.h"
 #include "module_hamilt_pw/hamilt_pwdft/global.h"
-#include "module_io/input.h"
 #include "module_io/input_conv.h"
 
 namespace elecstate
@@ -80,7 +79,8 @@ void H_TDDFT_pw::cal_fixed_v(double *vl_pseudo)
     ModuleBase::TITLE("H_TDDFT_pw", "cal_fixed_v");
 
     //skip if velocity_gague
-    if(stype==1)return;
+    if(stype==1) {return;
+}
 
     // time evolve
     H_TDDFT_pw::istep++;
@@ -233,7 +233,7 @@ int H_TDDFT_pw::check_ncut(int t_type)
     }
     return ncut;
 }
-void H_TDDFT_pw::update_At(void)
+void H_TDDFT_pw::update_At()
 {
     std::cout << "calculate electric potential" << std::endl;
     // time evolve
@@ -345,7 +345,8 @@ double H_TDDFT_pw::cal_v_time_Gauss(const bool last)
 
     double gauss_t = (istep_int - t0 * ncut) * dt_int;
     vext_time = cos(omega * gauss_t + phase) * exp(-gauss_t * gauss_t * 0.5 / (sigma * sigma)) * amp;
-    if(last)gauss_count++;
+    if(last) {gauss_count++;
+}
 
     return vext_time;
 }
@@ -375,7 +376,8 @@ double H_TDDFT_pw::cal_v_time_trapezoid(const bool last)
     }
 
     vext_time = vext_time * amp * cos(omega * istep_int * dt_int + phase);
-    if(last)trape_count++;
+    if(last) {trape_count++;
+}
 
     return vext_time;
 }
@@ -392,7 +394,8 @@ double H_TDDFT_pw::cal_v_time_trigonometric(const bool last)
     const double timenow = istep_int * dt_int;
 
     vext_time = amp * cos(omega1 * timenow + phase1) * sin(omega2 * timenow + phase2) * sin(omega2 * timenow + phase2);
-    if(last)trigo_count++;
+    if(last) {trigo_count++;
+}
 
     return vext_time;
 }
@@ -402,10 +405,11 @@ double H_TDDFT_pw::cal_v_time_heaviside()
     double t0 = *(heavi_t0.begin() + heavi_count);
     double amp = *(heavi_amp.begin() + heavi_count);
     double vext_time = 0.0;
-    if (istep < t0)
+    if (istep < t0) {
         vext_time = amp;
-    else if (istep >= t0)
+    } else if (istep >= t0) {
         vext_time = 0.0;
+}
     heavi_count++;
 
     return vext_time;

source/module_elecstate/potentials/H_TDDFT_pw.h

@@ -1,7 +1,6 @@
 #ifndef H_TDDFT_PW_H
 #define H_TDDFT_PW_H
 
-#include "module_io/input.h"
 #include "module_io/input_conv.h"
 #include "pot_base.h"
 

source/module_esolver/esolver.h

@@ -3,7 +3,6 @@
 
 #include "module_base/matrix.h"
 #include "module_cell/unitcell.h"
-#include "module_io/input.h"
 #include "module_parameter/parameter.h"
 
 namespace ModuleESolver

source/module_esolver/esolver_fp.cpp

@@ -2,7 +2,6 @@
 
 #include "module_base/global_variable.h"
 #include "module_hamilt_pw/hamilt_pwdft/global.h"
-#include "module_io/input.h"
 #include "module_parameter/parameter.h"
 namespace ModuleESolver
 {

source/module_esolver/esolver_ks.cpp

@@ -9,7 +9,6 @@
 #include <iostream>
 
 #include "module_base/timer.h"
-#include "module_io/input.h"
 #include "module_io/json_output/init_info.h"
 #include "module_io/print_info.h"
 #include "module_parameter/parameter.h"

source/module_esolver/esolver_sdft_pw.cpp

@@ -315,15 +315,13 @@ void ESolver_SDFT_PW::after_all_runners()
                                 this->p_hamilt,
                                 (hsolver::HSolverPW_SDFT*)phsol,
                                 &stowf);
-        sto_elecond
-            .decide_nche(PARAM.inp.cond_dt, INPUT.cond_dtbatch, 1e-8, this->nche_sto, PARAM.inp.emin_sto, PARAM.inp.emax_sto);
+        sto_elecond.decide_nche(PARAM.inp.cond_dt, 1e-8, this->nche_sto, PARAM.inp.emin_sto, PARAM.inp.emax_sto);
         sto_elecond.sKG(PARAM.inp.cond_smear,
                         PARAM.inp.cond_fwhm,
                         PARAM.inp.cond_wcut,
                         PARAM.inp.cond_dw,
                         PARAM.inp.cond_dt,
                         PARAM.inp.cond_nonlocal,
-                        INPUT.cond_dtbatch,
                         PARAM.inp.npart_sto);
     }
 }
@@ -353,7 +351,7 @@ void ESolver_SDFT_PW::nscf()
 
     const int iter = 1;
 
-    const double diag_thr = std::max(std::min(1e-5, 0.1 *PARAM.inp.scf_thr / std::max(1.0, GlobalV::nelec)), 1e-12);
+    const double diag_thr = std::max(std::min(1e-5, 0.1 * PARAM.inp.scf_thr / std::max(1.0, GlobalV::nelec)), 1e-12);
 
     std::cout << " DIGA_THR          : " << diag_thr << std::endl;
 

source/module_esolver/print_funcs.h

@@ -1,7 +1,7 @@
 #ifndef PRINT_FUNCTIONS_H
 #define PRINT_FUNCTIONS_H
 
-#include "module_io/input.h"
+#include "module_parameter/input_parameter.h"
 #include "module_basis/module_pw/pw_basis_k.h"
 
 namespace Print_functions

source/module_hamilt_general/module_surchem/test/setcell.h

@@ -7,7 +7,6 @@
 #include "module_cell/module_neighbor/sltk_atom_arrange.h"
 #include "module_cell/module_neighbor/sltk_grid_driver.h"
 #include "module_cell/unitcell.h"
-#include "module_io/input.h"
 #include "module_hamilt_pw/hamilt_pwdft/structure_factor.h"
 
 namespace GlobalC

source/module_hamilt_general/module_vdw/test/vdw_test.cpp

@@ -1,6 +1,5 @@
 #include "module_cell/unitcell.h"
 #include "module_cell/setup_nonlocal.h"
-#include "module_io/input.h"
 #include "module_base/mathzone.h"
 #include "module_base/vector3.h"
 #include"gtest/gtest.h"

source/module_hamilt_lcao/hamilt_lcaodft/FORCE_STRESS.cpp

@@ -764,7 +764,7 @@ void Force_Stress_LCAO<T>::calForcePwPart(ModuleBase::matrix& fvl_dvl,
     //--------------------------------------------------------
     // force due to self-consistent charge.
     //--------------------------------------------------------
-    f_pw.cal_force_scc(fscc, rhopw, vnew, vnew_exist);
+    f_pw.cal_force_scc(fscc, rhopw, vnew, vnew_exist,GlobalC::ucell);
     return;
 }
 

source/module_hamilt_lcao/hamilt_lcaodft/fvnl_dbeta_gamma.cpp

@@ -20,33 +20,38 @@ void Force_LCAO<double>::cal_fvnl_dbeta(const elecstate::DensityMatrix<double, d
     ModuleBase::TITLE("Force_LCAO", "cal_fvnl_dbeta");
     ModuleBase::timer::tick("Force_LCAO", "cal_fvnl_dbeta");
 
-    double r0[3];
-    double r1[3];
-
+// use schedule(dynamic) for load balancing because adj_num is various
+#pragma omp parallel
+{ 
+    ModuleBase::matrix local_svnl_dbeta(3, 3);
+    #pragma omp for schedule(dynamic)
     for (int iat = 0; iat < ucell.nat; iat++)
     {
         const int it = ucell.iat2it[iat];
         const int ia = ucell.iat2ia[iat];
         const int T0 = it;
         const int I0 = ia;
+        double r0[3]{0};
+        double r1[3]{0};
 
         const ModuleBase::Vector3<double> tau0 = ucell.atoms[it].tau[ia];
         // find ajacent atom of atom ia
         // gd.Find_atom( ucell.atoms[it].tau[ia] );
-        gd.Find_atom(ucell, tau0, it, ia);
+        AdjacentAtomInfo adjs;
+        gd.Find_atom(ucell, tau0, it, ia, &adjs);
         const double Rcut_Beta = ucell.infoNL.Beta[it].get_rcut_max();
 
         std::vector<std::unordered_map<int, std::vector<std::vector<double>>>> nlm_tot;
-        nlm_tot.resize(gd.getAdjacentNum() + 1); // this saves <psi_i|beta> and <psi_i|\nabla|beta>
+        nlm_tot.resize(adjs.adj_num + 1); // this saves <psi_i|beta> and <psi_i|\nabla|beta>
 
         // Step 1 : Calculate and save <psi_i|beta> and <psi_i|\nabla|beta>
-        for (int ad1 = 0; ad1 < gd.getAdjacentNum() + 1; ad1++)
+        for (int ad1 = 0; ad1 < adjs.adj_num + 1; ad1++)
         {
-            const int T1 = gd.getType(ad1);
+            const int T1 = adjs.ntype[ad1];
             const Atom* atom1 = &ucell.atoms[T1];
-            const int I1 = gd.getNatom(ad1);
+            const int I1 = adjs.natom[ad1];
             const int start1 = ucell.itiaiw2iwt(T1, I1, 0);
-            const ModuleBase::Vector3<double> tau1 = gd.getAdjacentTau(ad1);
+            const ModuleBase::Vector3<double> tau1 = adjs.adjacent_tau[ad1];
             const double Rcut_AO1 = orb.Phi[T1].getRcut();
 
             nlm_tot[ad1].clear();
@@ -87,22 +92,22 @@ void Force_LCAO<double>::cal_fvnl_dbeta(const elecstate::DensityMatrix<double, d
         } // ad
 
         // Step 2 : sum to get beta<psi_i|beta><beta|\nabla|psi_j>
-        for (int ad1 = 0; ad1 < gd.getAdjacentNum() + 1; ++ad1)
+        for (int ad1 = 0; ad1 < adjs.adj_num + 1; ++ad1)
         {
-            const int T1 = gd.getType(ad1);
+            const int T1 = adjs.ntype[ad1];
             const Atom* atom1 = &ucell.atoms[T1];
-            const int I1 = gd.getNatom(ad1);
+            const int I1 = adjs.natom[ad1];
             const int start1 = ucell.itiaiw2iwt(T1, I1, 0);
-            const ModuleBase::Vector3<double> tau1 = gd.getAdjacentTau(ad1);
+            const ModuleBase::Vector3<double> tau1 = adjs.adjacent_tau[ad1];
             const double Rcut_AO1 = orb.Phi[T1].getRcut();
 
-            for (int ad2 = 0; ad2 < gd.getAdjacentNum() + 1; ad2++)
+            for (int ad2 = 0; ad2 < adjs.adj_num + 1; ad2++)
             {
-                const int T2 = gd.getType(ad2);
+                const int T2 = adjs.ntype[ad2];
                 const Atom* atom2 = &ucell.atoms[T2];
-                const int I2 = gd.getNatom(ad2);
+                const int I2 = adjs.natom[ad2];
                 const int start2 = ucell.itiaiw2iwt(T2, I2, 0);
-                const ModuleBase::Vector3<double> tau2 = gd.getAdjacentTau(ad2);
+                const ModuleBase::Vector3<double> tau2 = adjs.adjacent_tau[ad2];
                 const double Rcut_AO2 = orb.Phi[T2].getRcut();
 
                 const double dist1 = (tau1 - tau0).norm() * ucell.lat0;
@@ -207,8 +212,8 @@ void Force_LCAO<double>::cal_fvnl_dbeta(const elecstate::DensityMatrix<double, d
                             {
                                 for (int jpol = ipol; jpol < 3; jpol++)
                                 {
-                                    svnl_dbeta(ipol, jpol)
-                                        += sum / 2.0 * (nlm[ipol] * r0[jpol] + nlm_t[ipol] * r1[jpol]);
+                                    local_svnl_dbeta(ipol, jpol)
+                                    += sum / 2.0 * (nlm[ipol] * r0[jpol] + nlm_t[ipol] * r1[jpol]);
                                 }
                             }
                         }
@@ -218,6 +223,21 @@ void Force_LCAO<double>::cal_fvnl_dbeta(const elecstate::DensityMatrix<double, d
         }             // ad1
     }                 // iat
 
+    // sum up local_svnl_dbeta to svnl_dbeta
+    if (isstress)
+    {
+        #pragma omp critical
+        {
+            for (int ipol = 0; ipol < 3; ipol++)
+            {
+                for (int jpol = ipol; jpol < 3; jpol++)
+                {
+                    svnl_dbeta(ipol, jpol) += local_svnl_dbeta(ipol, jpol);
+                }
+            }
+        }
+    }
+}
     if (isstress)
     {
         StressTools::stress_fill(ucell.lat0, ucell.omega, svnl_dbeta);