Feature: add DFT-1/2 and shell DFT-1/2, only support PW esolvers_ks_pw. (Useful information to support DFT-1/2 method with pw basis set)

source/Makefile.Objects

@@ -202,6 +202,8 @@ OBJS_CELL=atom_pseudo.o\
     bcast_cell.o\
     read_stru.o\
     read_atom_species.o\
+    sep.o\
+    sep_cell.o\
 
 OBJS_DEEPKS=LCAO_deepks.o\
         deepks_basic.o\
@@ -220,7 +222,7 @@ OBJS_DEEPKS=LCAO_deepks.o\
         deepks_phialpha.o\
         LCAO_deepks_io.o\
         LCAO_deepks_interface.o\
-        
+
 
 OBJS_ELECSTAT=elecstate.o\
     elecstate_energy_terms.o\
@@ -246,6 +248,7 @@ OBJS_ELECSTAT=elecstate.o\
     cal_nelec_nband.o\
     read_pseudo.o\
     cal_wfc.o\
+    pot_sep.o\
 
 OBJS_ELECSTAT_LCAO=elecstate_lcao.o\
       elecstate_lcao_cal_tau.o\
@@ -397,7 +400,7 @@ OBJS_HSOLVER=diago_cg.o\
     diag_const_nums.o\
     diag_hs_para.o\
     diago_pxxxgvx.o\
-    
+
 OBJS_HSOLVER_LCAO=hsolver_lcao.o\
       diago_scalapack.o\
       diago_lapack.o\
@@ -415,7 +418,7 @@ OBJS_HSOLVER_PEXSI=diago_pexsi.o\
       dist_bcd_matrix.o\
       dist_ccs_matrix.o\
       dist_matrix_transformer.o\
-      
+
 OBJS_MD=fire.o\
     langevin.o\
     md_base.o\
@@ -490,7 +493,7 @@ OBJS_RELAXATION=bfgs_basic.o\
     lbfgs.o\
     matrix_methods.o\
     line_search.o\
-    
+
 
 OBJS_SURCHEM=surchem.o\
     H_correction_pw.o\
@@ -690,7 +693,7 @@ OBJS_PARALLEL=parallel_common.o\
     parallel_kpoints.o\
     parallel_reduce.o\
     parallel_device.o
-      
+
 OBJS_SRCPW=H_Ewald_pw.o\
     dnrm2.o\
     VL_in_pw.o\
@@ -754,7 +757,8 @@ OBJS_SRCPW=H_Ewald_pw.o\
     sto_elecond.o\
     sto_dos.o\
     onsite_projector.o\
-    onsite_proj_tools.o
+    onsite_proj_tools.o\
+    VSep_in_pw.o
 
 OBJS_VDW=vdw.o\
     vdwd2_parameters.o\

source/source_cell/CMakeLists.txt

@@ -26,6 +26,8 @@ add_library(
     print_cell.cpp
     read_atom_species.cpp
     k_vector_utils.cpp
+    sep.cpp
+    sep_cell.cpp
 )
 
 if(ENABLE_COVERAGE)

source/source_cell/module_symmetry/test/symmetry_test_analysis.cpp

@@ -6,15 +6,15 @@
  * 1-1. test if the bravis lattice analysis is right;
  * 1-2. check if matrix-type and vector3-type;
  *     input and optimized lattice vectors are right;
- * 1-3. double-check for  if `gtrans_convert` 
+ * 1-3. double-check for  if `gtrans_convert`
  *     gives the same results as `veccon`;
  * 1-4. check `invmap` function gives the right result;
  * 1-5 test if `gmatrix_convert` and `gmatrix_convert_int`
  *     gives the right result;
  * 2. function: `atom_ordering_new:
  * test the new atom-sort algorithm gives the right result;
  *3. function: `pricell`:
- * test if the number of primitive cells are right, 
+ * test if the number of primitive cells are right,
  * using cases whose space group
  * is different from its point group.
  ***********************************************/
@@ -34,6 +34,10 @@ UnitCell::UnitCell(){}
 UnitCell::~UnitCell(){}
 Magnetism::Magnetism(){}
 Magnetism::~Magnetism() {}
+SepPot::SepPot(){}
+SepPot::~SepPot(){}
+Sep_Cell::Sep_Cell() noexcept {}
+Sep_Cell::~Sep_Cell() noexcept {}
 
 TEST_F(SymmetryTest, AnalySys)
 {
@@ -50,8 +54,8 @@ TEST_F(SymmetryTest, AnalySys)
         int cal_ibrav = symm.real_brav;
         EXPECT_EQ(cal_ibrav, ref_ibrav);
         EXPECT_EQ(cal_point_group, ref_point_group) << "ibrav=" << stru_lib[stru].ibrav;
-        
-        //2. input and optimized lattice, gtrans_convert and veccon 
+
+        //2. input and optimized lattice, gtrans_convert and veccon
         //input lattice
         EXPECT_EQ(symm.s1, ucell.a1);
         EXPECT_EQ(symm.s2, ucell.a2);
@@ -73,7 +77,7 @@ TEST_F(SymmetryTest, AnalySys)
         symm.gtrans_convert(symm.gtrans, gtrans_optconf.data(), symm.nrotk, ucell.latvec, symm.optlat);
         symm.veccon(gtrans_veccon, gtrans_optconf_veccon, symm.nrotk, symm.s1, symm.s2, symm.s3, symm.a1, symm.a2, symm.a3);
         for(int i=0;i<symm.nrotk;++i)
-            EXPECT_EQ(gtrans_optconf[i], ModuleBase::Vector3<double>(gtrans_optconf_veccon[i*3], 
+            EXPECT_EQ(gtrans_optconf[i], ModuleBase::Vector3<double>(gtrans_optconf_veccon[i*3],
             gtrans_optconf_veccon[i*3+1], gtrans_optconf_veccon[i*3+2]));
         delete[] gtrans_veccon;
         delete[] gtrans_optconf_veccon;
@@ -108,7 +112,7 @@ TEST_F(SymmetryTest, AnalySys)
         symm.gmatrix_convert_int(symm.gmatrix, gmatrix_opt, symm.nrotk, ucell.latvec, symm.optlat); //1->2
         symm.gmatrix_convert_int(gmatrix_opt, gmatrix_input_back, symm.nrotk, symm.optlat, ucell.latvec);   //2->3
         symm.gmatrix_convert_int(gmatrix_input_back, gmatrix_opt_back, symm.nrotk, ucell.latvec, symm.optlat); //3->4
-        
+
         symm.gmatrix_convert(symm.gmatrix, kgmatrix_nonint, symm.nrotk, ucell.latvec, ucell.G);
         for (int i=0;i<symm.nrotk;++i)
         {
@@ -130,7 +134,7 @@ TEST_F(SymmetryTest, AnalySys)
             EXPECT_NEAR(gmatrix_opt[i].e31, gmatrix_opt_back[i].e31, DOUBLETHRESHOLD);
             EXPECT_NEAR(gmatrix_opt[i].e23, gmatrix_opt_back[i].e23, DOUBLETHRESHOLD);
             EXPECT_NEAR(gmatrix_opt[i].e32, gmatrix_opt_back[i].e32, DOUBLETHRESHOLD);
-            
+
             ModuleBase::Matrix3 tmpA=symm.optlat.Inverse()*gmatrix_opt[i]*symm.optlat; //A^-1*SA*A
             ModuleBase::Matrix3 tmpB=ucell.latvec.Inverse()*symm.gmatrix[i]*ucell.latvec;//B^-1*SB*B
             ModuleBase::Matrix3 tmpG_int=ucell.G.Inverse()*symm.kgmatrix[i]*ucell.G;//G^-1*SG*G
@@ -168,7 +172,7 @@ TEST_F(SymmetryTest, AnalySys)
             EXPECT_NEAR(tmpA.e23, tmpG.e23, DOUBLETHRESHOLD);
             EXPECT_NEAR(tmpA.e32, tmpG.e32, DOUBLETHRESHOLD);
         }
-            
+
         delete[] gmatrix_input_back;
         delete[] gmatrix_opt;
         delete[] gmatrix_opt_back;
@@ -180,7 +184,7 @@ TEST_F(SymmetryTest, AnalySys)
 
 TEST_F(SymmetryTest, AtomOrderingNew)
 {
-    // the old function `atom_ordering` has bugs 
+    // the old function `atom_ordering` has bugs
     // so here I do not compare with its results
     ModuleSymmetry::Symmetry symm;
     symm.epsilon=1e-5;
@@ -199,7 +203,7 @@ TEST_F(SymmetryTest, AtomOrderingNew)
     }
     //ordering
     symm.test_atom_ordering(new_pos, nat, subindex);
-    //check 
+    //check
     for (int i=0;i<nat-1;++i)
     {
         //x[i]<=x[i+1]
@@ -250,7 +254,7 @@ TEST_F(SymmetryTest, SG_Pricell)
         std::string cal_point_group = symm.pgname;
         std::string ref_space_group = supercell_lib[stru].space_group;
         std::string cal_space_group = symm.spgname;
-        
+
         int ref_ncells = supercell_lib[stru].ibrav;
         EXPECT_EQ(symm.ncell, ref_ncells);
         EXPECT_EQ(cal_point_group, ref_point_group);

source/source_cell/module_symmetry/test/symmetry_test_symtrz.cpp

@@ -19,6 +19,10 @@ UnitCell::UnitCell() {}
 UnitCell::~UnitCell() {}
 Magnetism::Magnetism() {}
 Magnetism::~Magnetism() {}
+SepPot::SepPot(){}
+SepPot::~SepPot(){}
+Sep_Cell::Sep_Cell() noexcept {}
+Sep_Cell::~Sep_Cell() noexcept {}
 
 inline std::vector<double> allocate_pos(ModuleSymmetry::Symmetry& symm, UnitCell& ucell)
 {
@@ -46,7 +50,7 @@ TEST_F(SymmetryTest, ForceSymmetry)
 {
     auto check_force = [](stru_& conf, ModuleBase::matrix& force)
     {
-        // 1. check zeros  
+        // 1. check zeros
         for (auto iat : conf.force_zero_iat)
             for (int j = 0; j < 3; ++j)
                 EXPECT_NEAR(force(iat, j), 0.0, DOUBLETHRESHOLD);

source/source_cell/sep.cpp

@@ -0,0 +1,122 @@
+#include "sep.h"
+
+#include "source_base/global_variable.h"
+#include "source_base/parallel_common.h"
+#include "source_base/tool_title.h"
+#include "source_io/output.h"
+
+#include <fstream>
+#include <sstream>
+#include <string>
+
+SepPot::SepPot()
+{
+}
+
+SepPot::~SepPot()
+{
+    delete[] r;
+    r = nullptr;
+    delete[] rv;
+    rv = nullptr;
+}
+
+int SepPot::read_sep(std::ifstream& ifs)
+{
+    std::string line;
+    while (std::getline(ifs, line))
+    {
+        std::istringstream iss(line);
+        std::string key;
+        iss >> key;
+
+        if (key == "Sep.Element")
+        {
+            iss >> label;
+        }
+        else if (key == "Sep.XcType")
+        {
+            iss >> xc_type;
+        }
+        else if (key == "Sep.Orbital")
+        {
+            iss >> orbital;
+        }
+        else if (key == "Sep.Points")
+        {
+            iss >> mesh;
+            delete[] r;
+            r = new double[mesh];
+            delete[] rv;
+            rv = new double[mesh];
+        }
+        else if (key == "Sep.StripAmount")
+        {
+            iss >> strip_elec;
+        }
+        else if (key == "<Sep.Potential")
+        {
+            double r_val, rv_val;
+            int idx = 0;
+            while (std::getline(ifs, line) && line != "Sep.Potential>")
+            {
+                std::istringstream data_line(line);
+                if (data_line >> r_val >> rv_val)
+                {
+                    r[idx] = r_val;
+                    rv[idx] = rv_val;
+                    idx++;
+                }
+            }
+            break;
+        }
+    }
+    return 0;
+}
+
+void SepPot::print_sep_info(std::ofstream& ofs)
+{
+    ofs << "\n sep_vl:";
+    ofs << "\n sep_info:";
+    ofs << "\n label         " << label;
+    ofs << "\n xc            " << xc_type;
+    ofs << "\n orbital       " << orbital;
+    ofs << "\n strip electron" << strip_elec;
+}
+
+void SepPot::print_sep_vsep(std::ofstream& ofs)
+{
+    ofs << "\n mesh  " << mesh;
+    output::printr1_d(ofs, " r : ", r, mesh);
+    output::printr1_d(ofs, " vsep : ", rv, mesh);
+    ofs << "\n -----------------------------";
+}
+
+#ifdef __MPI
+
+void SepPot::bcast_sep()
+{
+    ModuleBase::TITLE("SepPot", "bcast_sep");
+    Parallel_Common::bcast_bool(is_enable);
+    Parallel_Common::bcast_double(r_in);
+    Parallel_Common::bcast_double(r_out);
+    Parallel_Common::bcast_double(r_power);
+    Parallel_Common::bcast_double(enhence_a);
+    Parallel_Common::bcast_string(label);
+    Parallel_Common::bcast_string(xc_type);
+    Parallel_Common::bcast_string(orbital);
+    Parallel_Common::bcast_int(strip_elec);
+    Parallel_Common::bcast_int(mesh);
+
+    if (GlobalV::MY_RANK != 0 && mesh > 0)
+    {
+        r = new double[mesh];
+        rv = new double[mesh];
+    }
+
+    Parallel_Common::bcast_double(r, mesh);
+    Parallel_Common::bcast_double(rv, mesh);
+
+    return;
+}
+#endif // __MPI

source/source_cell/sep.h

@@ -0,0 +1,39 @@
+#ifndef SEP_H
+#define SEP_H
+
+#include <fstream>
+#include <string>
+
+/**
+ * Sep Potential for DFT-1/2 etc.
+ *
+ * Sep Potential
+ */
+class SepPot
+{
+  public:
+    SepPot();
+    ~SepPot();
+
+    bool is_enable = false;
+    double r_in = 0.0;      /**< cut-off radius inner */
+    double r_out = 0.0;     /**< cut-off radius outter */
+    double r_power = 20.0;  /**< shell function exp factor */
+    double enhence_a = 1.0; /**< scale sep potential */
+    std::string label;      /**< element nameof sep  */
+    std::string xc_type;    /**< Exch-Corr type */
+    std::string orbital;    /** atomic angular moment s,p,d,f */
+    int mesh = 0;           /**< number of points in radial mesh */
+    int strip_elec = 0;     /**< strip electron amount 1->0.01 50->0.5 */
+    double* r = nullptr;    /**< ridial mesh */
+    double* rv = nullptr;   /**< sep potential, but rV, unit: Ry */
+
+    int read_sep(std::ifstream& is);
+    void print_sep_info(std::ofstream& ofs);
+    void print_sep_vsep(std::ofstream& ofs);
+#ifdef __MPI
+    void bcast_sep();
+#endif /* ifdef __MPI */
+};
+
+#endif /* ifndef SEP_H */