Merge branch 'develop' into warning

Zhuxuegang2022 · Feb 29, 2024 · a3b6408 · a3b6408
2 parents b541040 + 67474ff
commit a3b6408
Show file tree

Hide file tree

Showing 76 changed files with 809 additions and 262 deletions.
diff --git a/.gitignore b/.gitignore
@@ -20,4 +20,5 @@ dist
 toolchain.tar.gz
 time.json
 *.pyc
-__pycache__
+__pycache__
+abacus.json
diff --git a/docs/advanced/input_files/input-main.md b/docs/advanced/input_files/input-main.md
@@ -17,6 +17,7 @@
     - [init\_chg](#init_chg)
     - [init\_vel](#init_vel)
     - [nelec](#nelec)
+    - [nelec_delta](#nelec_delta)
     - [nupdown](#nupdown)
     - [dft\_functional](#dft_functional)
     - [xc\_temperature](#xc_temperature)
@@ -550,6 +551,13 @@ These variables are used to control general system parameters.
   - `>0.0`: this denotes the total number of electrons in the system. Must be less than 2*nbands. 
 - **Default**: 0.0
 
+### nelec_delta
+
+- **Type**: Real
+- **Description**: 
+ the total number of electrons will be calculated by `nelec`+`nelec_delta`.
+- **Default**: 0.0
+
 ### nupdown
 
 - **Type**: Real
@@ -3524,7 +3532,8 @@ These variables are used to control the usage of QO analysis. Please note presen
   - `minimal-nodeless`: according to principle quantum number of the highest occupied state, generate only nodeless orbitals, for example Cu, only generate 1s, 2p, 3d and 4f orbitals (for Cu, 4s is occupied, thus $n_{max} = 4$)
   - `minimal-valence`: according to principle quantum number of the highest occupied state, generate only orbitals with highest principle quantum number, for example Cu, only generate 4s, 4p, 4d and 4f orbitals.
   - `full`: similarly according to the maximal principle quantum number, generate all possible orbitals, therefore for Cu, for example, will generate 1s, 2s, 2p, 3s, 3p, 3d, 4s, 4p, 4d, 4f.
-  - `energy`: will generate hydrogen-like orbitals according to Aufbau principle. For example the Cu (1s2 2s2 2p6 3s2 3p6 3d10 4s1), will generate these orbitals.
+  - `energy-full`: will generate hydrogen-like orbitals according to Aufbau principle. For example the Cu (1s2 2s2 2p6 3s2 3p6 3d10 4s1), will generate these orbitals.
+  - `energy-valence`: from the highest n (principal quantum number) layer and n-1 layer, generate all occupied and possible ls (angular momentum quantum number) for only once, for example Cu, will generate 4s, 3d and 3p orbitals.
 
   For `qo_basis pswfc`
   - `all`: use all possible pseudowavefunctions in pseudopotential file.

diff --git a/examples/spin_polarized/ATOM/run.sh b/examples/spin_polarized/ATOM/run.sh
@@ -4,7 +4,7 @@ ABACUS_PATH=$(awk -F "=" '$1=="ABACUS_PATH"{print $2}' ../../SETENV)
 ABACUS_NPROCS=$(awk -F "=" '$1=="ABACUS_NPROCS"{print $2}' ../../SETENV)
 ABACUS_THREADS=$(awk -F "=" '$1=="ABACUS_THREADS"{print $2}' ../../SETENV)
 
-OMP_NUM_THREADS=${ABACUS_THREADS} mpirun -np ${ABACUS_NPROCS} ${ABACUS_PATH} | tee scf.output
+OMP_NUM_THREADS=${ABACUS_THREADS} mpirun -np 1 ${ABACUS_PATH} | tee scf.output
 
 if [[ ! -f scf.output ]] || 
    [[ ! -f OUT.ABACUS/running_scf.log ]] ||

diff --git a/source/Makefile.Objects b/source/Makefile.Objects
@@ -452,6 +452,7 @@ OBJS_IO=input.o\
     para_json.o\
     abacusjson.o\
     general_info.o\
+    init_info.o\
 
 
 OBJS_IO_LCAO=cal_r_overlap_R.o\

diff --git a/source/driver.cpp b/source/driver.cpp
@@ -42,7 +42,8 @@ void Driver::init()
     INPUT.close_log();
 
     // (5) output the json file
-    Json::json_output();
+    //Json::create_Json(&GlobalC::ucell.symm,GlobalC::ucell.atoms,&INPUT);
+    Json::create_Json(&GlobalC::ucell,&INPUT);
     return;
 }
 

diff --git a/source/module_base/global_variable.cpp b/source/module_base/global_variable.cpp
@@ -276,6 +276,7 @@ bool psi_initializer = false;
 
 int out_chg = 0;
 double nelec = 0;
+double nelec_delta = 0;
 bool out_bandgap = false; // QO added for bandgap printing
 int out_interval = 1;    // convert from out_hsR_interval liuyu 2023-04-18
 

diff --git a/source/module_base/global_variable.h b/source/module_base/global_variable.h
@@ -308,6 +308,7 @@ extern bool psi_initializer;
 extern int out_chg;
 
 extern double nelec;
+extern double nelec_delta;
 extern bool out_bandgap;
 extern int out_interval;
 

diff --git a/source/module_basis/module_nao/hydrogen_radials.cpp b/source/module_basis/module_nao/hydrogen_radials.cpp
@@ -136,36 +136,40 @@ double HydrogenRadials::generate_hydrogen_radial_toconv(const double charge,
         if(istep == 0) printf("%10d%12.2f%14.10f%18.10e\n", istep, rmax_, norm, delta_norm);
         ++istep;
     }
+    printf("...\n");
     printf("%10d%12.2f%14.10f%18.10e\n", istep, rmax_, norm, delta_norm);
     return rmax_;
 }
 
 std::vector<std::pair<int, int>> HydrogenRadials::unzip_strategy(const int nmax,
                                                                  const std::string strategy)
 {
-    if(strategy != "energy")
+    if(strategy.substr(0, 6) != "energy")
     {
         // because for "energy", the nmax is used as the number of electrons
         assert(nmax < 8);
     }
     std::vector<std::pair<int, int>> nl_pairs;
-    if(strategy == "minimal-nodeless")
+    if(strategy.substr(0, 7) == "minimal")
     {
-        for(int n = 1; n <= nmax; n++)
+        if(strategy == "minimal-nodeless")
         {
-            std::pair<int, int> nl_pair = std::make_pair(n, n - 1);
-            nl_pairs.push_back(nl_pair);
+            for(int n = 1; n <= nmax; n++)
+            {
+                std::pair<int, int> nl_pair = std::make_pair(n, n - 1);
+                nl_pairs.push_back(nl_pair);
+            }
         }
-    }
-    else if(strategy == "minimal-valence")
-    {
-        for(int l = 0; l < nmax; l++)
+        else// if(strategy == "minimal-valence")
         {
-            std::pair<int, int> nl_pair = std::make_pair(nmax, l);
-            nl_pairs.push_back(nl_pair);
+            for(int l = 0; l < nmax; l++)
+            {
+                std::pair<int, int> nl_pair = std::make_pair(nmax, l);
+                nl_pairs.push_back(nl_pair);
+            }
         }
     }
-    else if(strategy == "energy")
+    else if(strategy.substr(0, 6) == "energy")
     {
         // 1s, -(n+1)-> 2s, -(l+1)-> 2p, 3s, -(n+1)-> 3p, 4s, -(l+1)-> 3d, 4p, 5s, -(n+1)-> 4d, 5p, 6s
         int starting_n = 1;
@@ -194,6 +198,43 @@ std::vector<std::pair<int, int>> HydrogenRadials::unzip_strategy(const int nmax,
             }
             nl_switch++;
         }
+        if(strategy == "energy-valence")
+        {
+            std::vector<int> nmax_ls;
+            std::vector<int> nmax_minus1_ls;
+            int real_nmax = 0;
+            for(auto nl_pair : nl_pairs)
+            {
+                if(nl_pair.first > real_nmax) real_nmax = nl_pair.first;
+            }
+            for(auto it = nl_pairs.begin(); it != nl_pairs.end();)
+            {
+                if(it->first == real_nmax) nmax_ls.push_back(it->second);
+                else if(it->first == real_nmax - 1) nmax_minus1_ls.push_back(it->second);
+                else
+                {
+                    it = nl_pairs.erase(it);
+                    continue;
+                }
+                ++it;
+            }
+            for(auto it = nl_pairs.begin(); it != nl_pairs.end();)
+            {
+                if(it->first == real_nmax - 1)
+                {
+                    if(std::find(nmax_ls.begin(), nmax_ls.end(), it->second) != nmax_ls.end())
+                    {
+                        it = nl_pairs.erase(it);
+                        continue;
+                    }
+                }
+                ++it;
+            }
+        }
+        for(auto nl_pair : nl_pairs)
+        {
+            std::cout << nl_pair.first << " " << nl_pair.second << std::endl;
+        }
     }
     else
     {

diff --git a/source/module_basis/module_nao/pswfc_radials.cpp b/source/module_basis/module_nao/pswfc_radials.cpp
@@ -165,6 +165,7 @@ double PswfcRadials::cut_to_convergence(const std::vector<double>& rgrid,
         if(istep == 1) printf("%10d%12.2f%14.10f%18.10e\n", istep, rgrid[ir_], norm, delta_norm);
         ++istep;
     }
+    printf("...\n");
     printf("%10d%12.2f%14.10f%18.10e\n", istep, rgrid[ir_], norm, delta_norm);
 
     rvalue = std::vector<double>(rvalue.begin() + ir_min_, rvalue.begin() + ir_ + 1);
@@ -185,7 +186,8 @@ void PswfcRadials::smooth(std::vector<double>& rgrid,
     }
 }
 
-std::vector<double> PswfcRadials::pswfc_prepossess(std::map<std::pair<int, int>, std::vector<double>>& lzeta_rvalues)
+std::vector<double> PswfcRadials::pswfc_prepossess(std::map<std::pair<int, int>, std::vector<double>>& lzeta_rvalues,
+                                                   const double conv_thr)
 {
     double nmax = 0.0;
     for(auto it = lzeta_rvalues.begin(); it != lzeta_rvalues.end(); it++)
@@ -198,7 +200,7 @@ std::vector<double> PswfcRadials::pswfc_prepossess(std::map<std::pair<int, int>,
         {
             rgrid[ir] = ir * 0.01;
         }
-        double rcut_i = cut_to_convergence(rgrid, rvalue, 1e-6);
+        double rcut_i = cut_to_convergence(rgrid, rvalue, conv_thr);
         if(rvalue.size() > nmax) nmax = rvalue.size();
         lzeta_rvalues[it->first] = rvalue; // stores in map
     }
@@ -289,7 +291,7 @@ void PswfcRadials::read_upf_pswfc(std::ifstream& ifs,
     indexing(); // build index_map_
 
     // cut rvalue to convergence and generate rgrid
-    std::vector<double> rgrid = pswfc_prepossess(result);
+    std::vector<double> rgrid = pswfc_prepossess(result, conv_thr);
     // refresh ngird value
     ngrid = rgrid.size();
 

diff --git a/source/module_basis/module_nao/pswfc_radials.h b/source/module_basis/module_nao/pswfc_radials.h
@@ -82,7 +82,8 @@ class PswfcRadials : public RadialSet {
         /// @brief call cut_to_convergence for each (l,zeta) corresponding orbital in std::map, then zero-padding to the maximal r, generate a grid
         /// @param pswfc_map a map of (l,zeta) corresponding orbital
         /// @return a vector of radial grid
-        std::vector<double> pswfc_prepossess(std::map<std::pair<int, int>, std::vector<double>>& pswfc_map);
+        std::vector<double> pswfc_prepossess(std::map<std::pair<int, int>, std::vector<double>>& pswfc_map,
+                                             const double conv_thr = 1e-6);
     private:
 
 };

diff --git a/source/module_basis/module_nao/test/hydrogen_radials_test.cpp b/source/module_basis/module_nao/test/hydrogen_radials_test.cpp
@@ -56,6 +56,58 @@ TEST_F(HydrogenRadialsTest, UnzipStrategy)
     EXPECT_EQ(nl_pairs[2].second, 2);
     EXPECT_EQ(nl_pairs[3].first, 4);
     EXPECT_EQ(nl_pairs[3].second, 3);
+    nl_pairs = hr.unzip_strategy(4, "minimal-valence");
+    EXPECT_EQ(nl_pairs.size(), 4);
+    EXPECT_EQ(nl_pairs[0].first, 4);
+    EXPECT_EQ(nl_pairs[0].second, 0);
+    EXPECT_EQ(nl_pairs[1].first, 4);
+    EXPECT_EQ(nl_pairs[1].second, 1);
+    EXPECT_EQ(nl_pairs[2].first, 4);
+    EXPECT_EQ(nl_pairs[2].second, 2);
+    EXPECT_EQ(nl_pairs[3].first, 4);
+    EXPECT_EQ(nl_pairs[3].second, 3);
+    nl_pairs = hr.unzip_strategy(1, "energy-full"); // H, 1s1 -> 1s -> [1s]
+    EXPECT_EQ(nl_pairs.size(), 1);
+    EXPECT_EQ(nl_pairs[0].first, 1);
+    EXPECT_EQ(nl_pairs[0].second, 0);
+    nl_pairs = hr.unzip_strategy(6, "energy-full"); // C, 1s1 2s2 2p2 -> 1s2s2p -> [1s][2s][2p]
+    EXPECT_EQ(nl_pairs.size(), 3);
+    EXPECT_EQ(nl_pairs[0].first, 1); // 1s
+    EXPECT_EQ(nl_pairs[0].second, 0);
+    EXPECT_EQ(nl_pairs[1].first, 2); // 2s
+    EXPECT_EQ(nl_pairs[1].second, 0);
+    EXPECT_EQ(nl_pairs[2].first, 2); // 2p
+    EXPECT_EQ(nl_pairs[2].second, 1);
+    nl_pairs = hr.unzip_strategy(29, "energy-full"); // Cu, 1s1 2s2 2p6 3s2 3p6 3d10 4s1 -> 1s2s2p3s3p4s3d -> [1s][2s][2p3s][3p4s][3d]
+    EXPECT_EQ(nl_pairs.size(), 7);
+    EXPECT_EQ(nl_pairs[0].first, 1); // 1s
+    EXPECT_EQ(nl_pairs[0].second, 0);
+    EXPECT_EQ(nl_pairs[1].first, 2); // 2s
+    EXPECT_EQ(nl_pairs[1].second, 0);
+    EXPECT_EQ(nl_pairs[2].first, 2); // 2p
+    EXPECT_EQ(nl_pairs[2].second, 1);
+    EXPECT_EQ(nl_pairs[3].first, 3); // 3s
+    EXPECT_EQ(nl_pairs[3].second, 0);
+    EXPECT_EQ(nl_pairs[4].first, 3); // 3p
+    EXPECT_EQ(nl_pairs[4].second, 1);
+    EXPECT_EQ(nl_pairs[5].first, 4); // 4s
+    EXPECT_EQ(nl_pairs[5].second, 0);
+    EXPECT_EQ(nl_pairs[6].first, 3); // 3d
+    EXPECT_EQ(nl_pairs[6].second, 2);
+    nl_pairs = hr.unzip_strategy(29, "energy-valence"); // Cu, 1s1 2s2 2p6 3s2 3p6 3d10 4s1 -> 3p4s3d
+    EXPECT_EQ(nl_pairs.size(), 3);
+    EXPECT_EQ(nl_pairs[0].first, 3); // 3p
+    EXPECT_EQ(nl_pairs[0].second, 1);
+    EXPECT_EQ(nl_pairs[1].first, 4); // 4s
+    EXPECT_EQ(nl_pairs[1].second, 0);
+    EXPECT_EQ(nl_pairs[2].first, 3); // 3d
+    EXPECT_EQ(nl_pairs[2].second, 2);
+    nl_pairs = hr.unzip_strategy(14, "energy-valence"); // Si, 1s1 2s2 2p6 3s2 3p2 -> 3s3p
+    EXPECT_EQ(nl_pairs.size(), 2);
+    EXPECT_EQ(nl_pairs[0].first, 3); // 3s
+    EXPECT_EQ(nl_pairs[0].second, 0);
+    EXPECT_EQ(nl_pairs[1].first, 3); // 3p
+    EXPECT_EQ(nl_pairs[1].second, 1);
     // full, 1s
     nl_pairs = hr.unzip_strategy(1, "full");
     EXPECT_EQ(nl_pairs.size(), 1);
@@ -84,34 +136,6 @@ TEST_F(HydrogenRadialsTest, UnzipStrategy)
     EXPECT_EQ(nl_pairs[8].second, 2);
     EXPECT_EQ(nl_pairs[9].first, 4);
     EXPECT_EQ(nl_pairs[9].second, 3);
-    nl_pairs = hr.unzip_strategy(1, "energy"); // H, 1s1 -> 1s -> [1s]
-    EXPECT_EQ(nl_pairs.size(), 1);
-    EXPECT_EQ(nl_pairs[0].first, 1);
-    EXPECT_EQ(nl_pairs[0].second, 0);
-    nl_pairs = hr.unzip_strategy(6, "energy"); // C, 1s1 2s2 2p2 -> 1s2s2p -> [1s][2s][2p]
-    EXPECT_EQ(nl_pairs.size(), 3);
-    EXPECT_EQ(nl_pairs[0].first, 1); // 1s
-    EXPECT_EQ(nl_pairs[0].second, 0);
-    EXPECT_EQ(nl_pairs[1].first, 2); // 2s
-    EXPECT_EQ(nl_pairs[1].second, 0);
-    EXPECT_EQ(nl_pairs[2].first, 2); // 2p
-    EXPECT_EQ(nl_pairs[2].second, 1);
-    nl_pairs = hr.unzip_strategy(29, "energy"); // Cu, 1s1 2s2 2p6 3s2 3p6 3d10 4s1 -> 1s2s2p3s3p4s3d -> [1s][2s][2p3s][3p4s][3d]
-    EXPECT_EQ(nl_pairs.size(), 7);
-    EXPECT_EQ(nl_pairs[0].first, 1); // 1s
-    EXPECT_EQ(nl_pairs[0].second, 0);
-    EXPECT_EQ(nl_pairs[1].first, 2); // 2s
-    EXPECT_EQ(nl_pairs[1].second, 0);
-    EXPECT_EQ(nl_pairs[2].first, 2); // 2p
-    EXPECT_EQ(nl_pairs[2].second, 1);
-    EXPECT_EQ(nl_pairs[3].first, 3); // 3s
-    EXPECT_EQ(nl_pairs[3].second, 0);
-    EXPECT_EQ(nl_pairs[4].first, 3); // 3p
-    EXPECT_EQ(nl_pairs[4].second, 1);
-    EXPECT_EQ(nl_pairs[5].first, 4); // 4s
-    EXPECT_EQ(nl_pairs[5].second, 0);
-    EXPECT_EQ(nl_pairs[6].first, 3); // 3d
-    EXPECT_EQ(nl_pairs[6].second, 2);
 }
 
 TEST_F(HydrogenRadialsTest, RadialNorm)
@@ -288,6 +312,69 @@ TEST_F(HydrogenRadialsTest, Build)
     EXPECT_EQ(hr.nzeta(0), 1);
     EXPECT_EQ(hr.nzeta_max(), 1);
     EXPECT_EQ(hr.nchi(), 3);
+    // Cu, minimal-valence, 4s 4p 4d
+    hr.build(
+        itype_,
+        charge_,
+        4,
+        rcut_,
+        dr_,
+        1e-6,
+        rank_,
+        "Cu",
+        "minimal-valence",
+        ptr_log_
+    );
+    // nmax = 4, minimal-valence, yields 4s 4p 4d 4f orbitals
+    EXPECT_EQ(hr.lmax(), 3);
+    EXPECT_EQ(hr.nzeta(0), 1);
+    EXPECT_EQ(hr.nzeta(1), 1);
+    EXPECT_EQ(hr.nzeta(2), 1);
+    EXPECT_EQ(hr.nzeta(3), 1);
+    EXPECT_EQ(hr.nzeta_max(), 1);
+    EXPECT_EQ(hr.nchi(), 4);
+    // Cu, energy-full, 1s 2s 2p 3s 3p 4s 3d
+    hr.build(
+        itype_,
+        charge_,
+        29,
+        rcut_,
+        dr_,
+        1e-6,
+        rank_,
+        "Cu",
+        "energy-full",
+        ptr_log_
+    );
+    // nmax = 29, energy-full, yields 1s 2s 2p 3s 3p 4s 3d orbitals
+    EXPECT_EQ(hr.lmax(), 2);
+    EXPECT_EQ(hr.nzeta(0), 4);
+    EXPECT_EQ(hr.nzeta(1), 2);
+    EXPECT_EQ(hr.nzeta(2), 1);
+    EXPECT_EQ(hr.nzeta(3), 0);
+    EXPECT_EQ(hr.nzeta_max(), 4);
+    EXPECT_EQ(hr.nchi(), 7);
+    // Cu, energy-valence, 3p 4s 3d
+    hr.build(
+        itype_,
+        charge_,
+        29,
+        rcut_,
+        dr_,
+        1e-6,
+        rank_,
+        "Cu",
+        "energy-valence",
+        ptr_log_
+    );
+    // nmax = 29, energy-valence, yields 3p 4s 3d orbitals
+    EXPECT_EQ(hr.lmax(), 2);
+    EXPECT_EQ(hr.nzeta(0), 1);
+    EXPECT_EQ(hr.nzeta(1), 1);
+    EXPECT_EQ(hr.nzeta(2), 1);
+    EXPECT_EQ(hr.nzeta(3), 0);
+    EXPECT_EQ(hr.nzeta_max(), 1);
+    EXPECT_EQ(hr.nchi(), 3);
     // build 1s 2s 2p 3s 3p 3d 4s 4p 4d 4f
     hr.build(
         itype_,