Merge branch 'develop' into subspace-dug

Author: haozhihan

docs/advanced/input_files/input-main.md

@@ -11,7 +11,6 @@
     - [kpar](#kpar)
     - [bndpar](#bndpar)
     - [latname](#latname)
-    - [psi\_initializer](#psi_initializer)
     - [init\_wfc](#init_wfc)
     - [init\_chg](#init_chg)
     - [init\_vel](#init_vel)
@@ -93,6 +92,7 @@
     - [scf\_os\_stop](#scf_os_stop)
     - [scf\_os\_thr](#scf_os_thr)
     - [scf\_os\_ndim](#scf_os_ndim)
+    - [sc\_os\_ndim](#sc_os_ndim)
     - [chg\_extrap](#chg_extrap)
     - [lspinorb](#lspinorb)
     - [noncolin](#noncolin)
@@ -467,7 +467,7 @@
     - [abs\_broadening](#abs_broadening)
     - [ri\_hartree\_benchmark](#ri_hartree_benchmark)
     - [aims\_nbasis](#aims_nbasis)
-  - [Reduced Density Matrix Functional Theory](#Reduced-Density-Matrix-Functional-Theory)
+  - [Reduced Density Matrix Functional Theory](#reduced-density-matrix-functional-theory)
     - [rdmft](#rdmft)
     - [rdmft\_power\_alpha](#rdmft_power_alpha)
 
@@ -580,17 +580,6 @@ These variables are used to control general system parameters.
   - triclinic: triclinic (14)
 - **Default**: none
 
-### psi_initializer
-
-- **Type**: Integer
-- **Description**: enable the experimental feature psi_initializer, to support use numerical atomic orbitals initialize wavefunction (`basis_type pw` case).
-
-  NOTE: this feature is not well-implemented for `nspin 4` case (closed presently), and cannot use with `calculation nscf`/`esolver_type sdft` cases.
-  Available options are:
-  - 0: disable psi_initializer
-  - 1: enable psi_initializer
-- **Default**: 0
-
 ### init_wfc
 
 - **Type**: String
@@ -602,8 +591,6 @@ These variables are used to control general system parameters.
   - atomic+random: add small random numbers on atomic pseudo-wavefunctions
   - file: from binary files `WAVEFUNC*.dat`, which are output by setting [out_wfc_pw](#out_wfc_pw) to `2`.
   - random: random numbers
-
-  with `psi_initializer 1`, two more options are supported:
   - nao: from numerical atomic orbitals. If they are not enough, other wave functions are initialized with random numbers.
   - nao+random: add small random numbers on numerical atomic orbitals
 - **Default**: atomic

source/Makefile.Objects

@@ -398,6 +398,7 @@ OBJS_PSI=psi.o\
 
 OBJS_PSI_INITIALIZER=psi_initializer.o\
                      psi_initializer_random.o\
+                     psi_initializer_file.o\
                      psi_initializer_atomic.o\
                      psi_initializer_atomic_random.o\
                      psi_initializer_nao.o\
@@ -494,6 +495,7 @@ OBJS_IO=input_conv.o\
     to_wannier90_lcao.o\
     fR_overlap.o\
     unk_overlap_pw.o\
+    write_pao.o\
     write_wfc_pw.o\
     winput.o\
     write_cube.o\
@@ -669,8 +671,6 @@ OBJS_SRCPW=H_Ewald_pw.o\
     of_stress_pw.o\
     symmetry_rho.o\
     symmetry_rhog.o\
-    wavefunc.o\
-    wf_atomic.o\
     psi_init.o\
     elecond.o\
     sto_tool.o\

source/driver_run.cpp

@@ -40,10 +40,10 @@ void Driver::driver_run()
     // the life of ucell should begin here, mohan 2024-05-12
     UnitCell ucell;
     ucell.setup(PARAM.inp.latname,
-                    PARAM.inp.ntype,
-                    PARAM.inp.lmaxmax,
-                    PARAM.inp.init_vel,
-                    PARAM.inp.fixed_axes);
+                PARAM.inp.ntype,
+                PARAM.inp.lmaxmax,
+                PARAM.inp.init_vel,
+                PARAM.inp.fixed_axes);
 
     ucell.setup_cell(PARAM.globalv.global_in_stru, GlobalV::ofs_running);
     Check_Atomic_Stru::check_atomic_stru(ucell, PARAM.inp.min_dist_coef);

source/module_basis/module_ao/test/ORB_unittest.cpp

@@ -28,7 +28,6 @@ void test_orb::TearDown()
 }
     ooo.clear_after_ions(OGT, ORB, 0, nproj);
     delete[] nproj;
-    delete[] orbital_fn;
     return;
 }
 
@@ -75,7 +74,7 @@ void test_orb::set_orbs()
     ORB.init(ofs_running,
                        ntype_read,
                        "./",
-                       orbital_fn,
+                       orbital_fn.data(),
                        descriptor_file,
                        lmax,
                        lcao_ecut,
@@ -114,7 +113,7 @@ void test_orb::set_files()
 
     ModuleBase::GlobalFunc::SCAN_BEGIN(ifs, "NUMERICAL_ORBITAL");
 
-    orbital_fn = new std::string[ntype_read];
+    orbital_fn.resize(ntype_read);
 
     for (int it = 0; it < ntype_read; it++)
     {

source/module_basis/module_ao/test/ORB_unittest.h

@@ -51,7 +51,7 @@ class test_orb : public testing::Test
     double randr(double Rmax);
     void gen_table_center2();
 
-    bool force_flag = 0;
+    bool force_flag = false;
     int my_rank = 0;
     int ntype_read;
 
@@ -66,7 +66,7 @@ class test_orb : public testing::Test
     int lmax = 1;
     double lat0 = 1.0;
     std::string case_dir = "./GaAs/";
-    std::string* orbital_fn;
+    std::vector<std::string> orbital_fn;
     std::string descriptor_file;
 };
 #endif

source/module_basis/module_pw/pw_basis_k.cpp

@@ -22,7 +22,6 @@ PW_Basis_K::~PW_Basis_K()
     delete[] igl2isz_k;
     delete[] igl2ig_k;
     delete[] gk2;
-    delete[] ig2ixyz_k_;
 #if defined(__CUDA) || defined(__ROCM)
     if (this->device == "gpu") {
         if (this->precision == "single") {
@@ -169,6 +168,7 @@ void PW_Basis_K::setupIndGk()
         syncmem_int_h2d_op()(gpu_ctx, cpu_ctx, this->d_igl2isz_k, this->igl2isz_k, this->npwk_max * this->nks);
     }
 #endif
+    this->get_ig2ixyz_k();
     return;
 }
 
@@ -334,8 +334,12 @@ int& PW_Basis_K::getigl2ig(const int ik, const int igl) const
 
 void PW_Basis_K::get_ig2ixyz_k()
 {
-    delete[] this->ig2ixyz_k_;
-    this->ig2ixyz_k_ = new int [this->npwk_max * this->nks];
+    if (this->device != "gpu")
+    {
+        //only GPU need to get ig2ixyz_k
+        return;
+    }
+    int * ig2ixyz_k_cpu = new int [this->npwk_max * this->nks];
     ModuleBase::Memory::record("PW_B_K::ig2ixyz", sizeof(int) * this->npwk_max * this->nks);
     assert(gamma_only == false); //We only finish non-gamma_only fft on GPU temperarily.
     for(int ik = 0; ik < this->nks; ++ik)
@@ -348,15 +352,12 @@ void PW_Basis_K::get_ig2ixyz_k()
             int ixy = this->is2fftixy[is];
             int iy = ixy % this->ny;
             int ix = ixy / this->ny;
-            ig2ixyz_k_[igl + ik * npwk_max] = iz + iy * nz + ix * ny * nz;
+            ig2ixyz_k_cpu[igl + ik * npwk_max] = iz + iy * nz + ix * ny * nz;
         }
     }
-#if defined(__CUDA) || defined(__ROCM)
-    if (this->device == "gpu") {
-        resmem_int_op()(gpu_ctx, ig2ixyz_k, this->npwk_max * this->nks);
-        syncmem_int_h2d_op()(gpu_ctx, cpu_ctx, this->ig2ixyz_k, this->ig2ixyz_k_, this->npwk_max * this->nks);
-    }
-#endif
+    resmem_int_op()(gpu_ctx, ig2ixyz_k, this->npwk_max * this->nks);
+    syncmem_int_h2d_op()(gpu_ctx, cpu_ctx, this->ig2ixyz_k, ig2ixyz_k_cpu, this->npwk_max * this->nks);
+    delete[] ig2ixyz_k_cpu;
 }
 
 std::vector<int> PW_Basis_K::get_ig2ix(const int ik) const

source/module_basis/module_pw/pw_basis_k.h

@@ -71,8 +71,6 @@ class PW_Basis_K : public PW_Basis
         const bool xprime_in = true
     );
 
-    void get_ig2ixyz_k();
-
   public:
     int nks=0;//number of k points in this pool
     ModuleBase::Vector3<double> *kvec_d=nullptr; // Direct coordinates of k points
@@ -88,8 +86,7 @@ class PW_Basis_K : public PW_Basis
 
     int *igl2isz_k=nullptr, * d_igl2isz_k = nullptr; //[npwk_max*nks] map (igl,ik) to (is,iz)
     int *igl2ig_k=nullptr;//[npwk_max*nks] map (igl,ik) to ig
-    int *ig2ixyz_k=nullptr;
-    int *ig2ixyz_k_=nullptr;
+    int *ig2ixyz_k=nullptr; ///< [npw] map ig to ixyz
 
     double *gk2=nullptr; // modulus (G+K)^2 of G vectors [npwk_max*nks]
 
@@ -108,6 +105,8 @@ class PW_Basis_K : public PW_Basis
     double * d_gk2 = nullptr; // modulus (G+K)^2 of G vectors [npwk_max*nks]
     //create igl2isz_k map array for fft
     void setupIndGk();
+    // get ig2ixyz_k
+    void get_ig2ixyz_k();
     //calculate G+K, it is a private function
     ModuleBase::Vector3<double> cal_GplusK_cartesian(const int ik, const int ig) const;
 

source/module_basis/module_pw/test/test4-4.cpp

@@ -213,13 +213,6 @@ TEST_F(PWTEST,test4_4)
             }
         }
 
-        //check getig2ixyz_k
-        pwtest.get_ig2ixyz_k();
-        for(int igl = 0; igl < npwk ; ++igl)
-        { 
-            EXPECT_GE(pwtest.ig2ixyz_k_[igl + ik * pwtest.npwk_max], 0);
-        }
-
     }
     delete []tmp; 
     delete [] rhor;

source/module_cell/bcast_cell.cpp

@@ -1,5 +1,10 @@
 #include "unitcell.h"   
-
+#include "module_base/parallel_common.h"
+#include "module_parameter/parameter.h"
+#ifdef __EXX
+#include "module_ri/serialization_cereal.h"
+#include "module_hamilt_pw/hamilt_pwdft/global.h"
+#endif
 namespace unitcell
 {
     void bcast_atoms_tau(Atom* atoms,
@@ -12,4 +17,105 @@ namespace unitcell
         }
     #endif
     }
+    
+    void bcast_atoms_pseudo(Atom* atoms,
+                                 const int ntype)
+    {
+    #ifdef __MPI
+        MPI_Barrier(MPI_COMM_WORLD);
+        for (int i = 0; i < ntype; i++) 
+        {
+            atoms[i].bcast_atom2();
+        }
+    #endif
+    }
+
+    void bcast_Lattice(Lattice& lat)
+    {
+    #ifdef __MPI
+        MPI_Barrier(MPI_COMM_WORLD);
+        // distribute lattice parameters.
+        ModuleBase::Matrix3& latvec = lat.latvec;
+        ModuleBase::Matrix3& latvec_supercell = lat.latvec_supercell;
+        Parallel_Common::bcast_string(lat.Coordinate);
+        Parallel_Common::bcast_double(lat.lat0);
+        Parallel_Common::bcast_double(lat.lat0_angstrom);
+        Parallel_Common::bcast_double(lat.tpiba);
+        Parallel_Common::bcast_double(lat.tpiba2);
+        Parallel_Common::bcast_double(lat.omega);
+        Parallel_Common::bcast_string(lat.latName);
+
+        // distribute lattice vectors.
+        Parallel_Common::bcast_double(latvec.e11);
+        Parallel_Common::bcast_double(latvec.e12);
+        Parallel_Common::bcast_double(latvec.e13);
+        Parallel_Common::bcast_double(latvec.e21);
+        Parallel_Common::bcast_double(latvec.e22);
+        Parallel_Common::bcast_double(latvec.e23);
+        Parallel_Common::bcast_double(latvec.e31);
+        Parallel_Common::bcast_double(latvec.e32);
+        Parallel_Common::bcast_double(latvec.e33);
+
+         // distribute lattice vectors.
+        for (int i = 0; i < 3; i++)
+        {
+            Parallel_Common::bcast_double(lat.a1[i]);
+            Parallel_Common::bcast_double(lat.a2[i]);
+            Parallel_Common::bcast_double(lat.a3[i]);
+            Parallel_Common::bcast_double(lat.latcenter[i]);
+            Parallel_Common::bcast_int(lat.lc[i]);
+        }
+
+        // distribute superlattice vectors.
+        Parallel_Common::bcast_double(latvec_supercell.e11);
+        Parallel_Common::bcast_double(latvec_supercell.e12);
+        Parallel_Common::bcast_double(latvec_supercell.e13);
+        Parallel_Common::bcast_double(latvec_supercell.e21);
+        Parallel_Common::bcast_double(latvec_supercell.e22);
+        Parallel_Common::bcast_double(latvec_supercell.e23);
+        Parallel_Common::bcast_double(latvec_supercell.e31);
+        Parallel_Common::bcast_double(latvec_supercell.e32);
+        Parallel_Common::bcast_double(latvec_supercell.e33);
+
+        // distribute Change the lattice vectors or not
+    #endif
+    }
+    
+    void bcast_magnetism(Magnetism& magnet, const int ntype)
+    {
+    #ifdef __MPI
+        MPI_Barrier(MPI_COMM_WORLD);
+        Parallel_Common::bcast_double(magnet.start_magnetization, ntype);
+        if (PARAM.inp.nspin == 4) 
+        {
+            Parallel_Common::bcast_double(magnet.ux_[0]);
+            Parallel_Common::bcast_double(magnet.ux_[1]);
+            Parallel_Common::bcast_double(magnet.ux_[2]);
+        }
+    #endif
+    }
+
+    void bcast_unitcell(UnitCell& ucell)
+    {
+    #ifdef __MPI
+        const int ntype = ucell.ntype;
+        Parallel_Common::bcast_int(ucell.nat);
+
+        bcast_Lattice(ucell.lat);
+        bcast_magnetism(ucell.magnet,ntype);
+        bcast_atoms_tau(ucell.atoms,ntype);
+
+        for (int i = 0; i < ntype; i++)
+        {
+            Parallel_Common::bcast_string(ucell.orbital_fn[i]);
+        }
+
+        #ifdef __EXX
+        ModuleBase::bcast_data_cereal(GlobalC::exx_info.info_ri.files_abfs,
+                                    MPI_COMM_WORLD,
+                                    0);
+        #endif
+        return;
+    #endif
+    }
 }

source/module_cell/bcast_cell.h

@@ -1,10 +1,50 @@
 #ifndef BCAST_CELL_H
 #define BCAST_CELL_H
 
+#include "module_cell/unitcell.h"
 namespace unitcell
 {
+    /**
+     * @brief broadcast the tau array of the atoms
+     * 
+     * @param atoms: the atoms to be broadcasted [in/out]
+     * @param ntype: the number of types of the atoms [in]
+     */
     void bcast_atoms_tau(Atom* atoms,
+                     const int ntype);
+                     
+    /**
+     * @brief broadcast the pseduo of the atoms
+     * 
+     * @param atoms: the atoms to be broadcasted [in/out]
+     * @param ntype:  the number of types of the atoms [in]
+     */
+    void bcast_atoms_pseudo(Atom* atoms,
+                            const int ntype);
+    /**
+     * @brief broadcast the lattice
+     * 
+     * @param lat: the lattice to be broadcasted [in/out]
+     */
+    void bcast_Lattice(Lattice& lat);
+
+    /**
+     * @brief broadcast the magnetism
+     * 
+     * @param magnet: the magnetism to be broadcasted [in/out]
+     * @param nytpe: the number of types of the atoms [in]
+     */
+    void bcast_magnetism(Magnetism& magnet,
                          const int ntype);
+    
+    /**
+     * @brief broadcast the unitcell
+     * 
+     * @param ucell: the unitcell to be broadcasted [in/out]
+     */
+    void bcast_unitcell(UnitCell& ucell);
+
+
 }
 
 #endif // BCAST_CELL_H