Add deepks_out_freq_elec for outputting deepks labels during electronic steps

docs/advanced/input_files/input-main.md

@@ -188,6 +188,7 @@
     - [bessel\_nao\_sigma](#bessel_nao_sigma)
   - [DeePKS](#deepks)
     - [deepks\_out\_labels](#deepks_out_labels)
+    - [deepks\_out\_freq\_elec](#deepks_out_freq_elec)
     - [deepks\_scf](#deepks_scf)
     - [deepks\_equiv](#deepks_equiv)
     - [deepks\_model](#deepks_model)
@@ -2159,6 +2160,13 @@ Warning: this function is not robust enough for the current version. Please try
   This is not needed when `deepks_out_labels` equals 2. 
 - **Default**: 0
 
+### deepks_out_freq_elec
+
+- **Type**: Integer
+- **Availability**: numerical atomic orbital basis
+- **Description**: When `deepks_out_freq_elec` is greater than 0, print labels and descriptors for DeePKS in OUT.${suffix}/DeePKS_Labels_Elec per `deepks_out_freq_elec` electronic iterations, with suffix `_e*` to distinguish different steps. Often used with `deepks_out_labels` equals 1.
+- **Default**: 0
+
 ### deepks_scf
 
 - **Type**: Boolean

source/module_hamilt_lcao/module_deepks/LCAO_deepks_interface.h

@@ -30,6 +30,8 @@ class LCAO_Deepks_Interface
     /// @param[in] psid
     /// @param[in] dm
     /// @param[in] p_ham
+    /// @param[in] iter
+    /// @param[in] conv_esolver
     /// @param[in] rank
     void out_deepks_labels(const double& etot,
                            const int& nks,
@@ -44,6 +46,8 @@ class LCAO_Deepks_Interface
                            const psi::Psi<TK>& psid,
                            const elecstate::DensityMatrix<TK, double>* dm,
                            hamilt::HamiltLCAO<TK, TR>* p_ham,
+                           const int& iter,
+                           const bool& conv_esolver, 
                            const int rank,
                            std::ostream& ofs_running);
 

source/module_hamilt_lcao/module_deepks/LCAO_deepks_io.cpp

@@ -80,7 +80,7 @@ void LCAO_deepks_io::save_npy_d(const int nat,
                                 const std::vector<int>& inl2l,
                                 const bool deepks_equiv,
                                 const std::vector<torch::Tensor>& descriptor,
-                                const std::string& out_dir,
+                                const std::string& dm_eig_file,
                                 const int rank)
 {
     ModuleBase::TITLE("LCAO_deepks_io", "save_npy_d");
@@ -107,9 +107,7 @@ void LCAO_deepks_io::save_npy_d(const int nat,
         const long unsigned dshape[] = {static_cast<unsigned long>(nat), static_cast<unsigned long>(des_per_atom)};
         if (rank == 0)
         {
-            std::string file_dm_eig = out_dir + "deepks_dm_eig.npy";
-            // std::string file_dm_eig = "dm_eig.npy";
-            npy::SaveArrayAsNumpy(file_dm_eig, false, 2, dshape, npy_des);
+            npy::SaveArrayAsNumpy(dm_eig_file, false, 2, dshape, npy_des);
         }
     }
     else
@@ -127,9 +125,7 @@ void LCAO_deepks_io::save_npy_d(const int nat,
         const long unsigned dshape[] = {static_cast<unsigned long>(nat), static_cast<unsigned long>(des_per_atom)};
         if (rank == 0)
         {
-            std::string file_dm_eig = out_dir + "deepks_dm_eig.npy";
-            // std::string file_dm_eig = "dm_eig.npy";
-            npy::SaveArrayAsNumpy(file_dm_eig, false, 2, dshape, npy_des);
+            npy::SaveArrayAsNumpy(dm_eig_file, false, 2, dshape, npy_des);
         }
     }
     return;

source/module_hamilt_lcao/module_deepks/LCAO_deepks_io.h

@@ -49,7 +49,7 @@ void save_npy_d(const int nat,
                 const std::vector<int>& inl2l,
                 const bool deepks_equiv,
                 const std::vector<torch::Tensor>& descriptor,
-                const std::string& out_dir,
+                const std::string& dm_eig_file,
                 const int rank);
 
 // save energy

source/module_hamilt_lcao/module_deepks/deepks_basic.cpp

@@ -137,13 +137,14 @@ void DeePKS_domain::cal_edelta_gedm_equiv(const int nat,
     ModuleBase::TITLE("DeePKS_domain", "cal_edelta_gedm_equiv");
     ModuleBase::timer::tick("DeePKS_domain", "cal_edelta_gedm_equiv");
 
+    const std::string file_d = PARAM.globalv.global_out_dir + "deepks_dm_eig.npy";;
     LCAO_deepks_io::save_npy_d(nat,
                                des_per_atom,
                                inlmax,
                                inl2l,
                                PARAM.inp.deepks_equiv,
                                descriptor,
-                               PARAM.globalv.global_out_dir,
+                               file_d,
                                rank); // libnpy needed
 
     if (rank == 0)

source/module_io/ctrl_output_lcao.cpp

@@ -173,6 +173,8 @@ void ctrl_output_lcao(UnitCell& ucell,
 			*psi,
 			pelec->get_DM(),
 			p_ham_deepks,
+            -1, // -1 when called in after scf
+            true, // no used when after scf
 			GlobalV::MY_RANK,
             GlobalV::ofs_running);
 #endif

source/module_io/read_input_item_deepks.cpp

@@ -13,6 +13,18 @@ void ReadInput::item_deepks()
         read_sync_int(input.deepks_out_labels);
         this->add_item(item);
     }
+    {
+        Input_Item item("deepks_out_freq_elec");
+        item.annotation = ">0 frequency of electronic iteration to output descriptors and labels for deepks.";
+        read_sync_int(input.deepks_out_freq_elec);
+        item.reset_value = [](const Input_Item& item, Parameter& para) {
+            if (para.input.deepks_out_freq_elec < 0)
+            {
+                para.input.deepks_out_freq_elec = 0;
+            }
+        };
+        this->add_item(item);
+    }
     {
         Input_Item item("deepks_scf");
         item.annotation = ">0 add V_delta to Hamiltonian";

source/module_io/read_set_globalv.cpp

@@ -97,6 +97,10 @@ void ReadInput::set_global_dir(const Input_para& inp, System_para& sys)
     sys.global_mlkedf_descriptor_dir = sys.global_out_dir + "MLKEDF_Descriptors/";
     sys.global_mlkedf_descriptor_dir = to_dir(sys.global_mlkedf_descriptor_dir);
 
+    /// get the global directory for DeePKS labels during electronic steps 
+    sys.global_deepks_label_elec_dir = sys.global_out_dir + "DeePKS_Labels_Elec/";
+    sys.global_deepks_label_elec_dir = to_dir(sys.global_deepks_label_elec_dir);
+
     /// get the global readin directory
     sys.global_readin_dir = inp.read_file_dir;
     sys.global_readin_dir = to_dir(sys.global_readin_dir);

source/module_parameter/input_parameter.h

@@ -263,6 +263,8 @@ struct Input_para
     //==========================================================
     int deepks_out_labels = 0;         ///< (need libnpy) prints energy and force labels and
                                        ///< descriptors for training, wenfei 2022-1-12
+    int deepks_out_freq_elec = 0;      ///< (need libnpy) frequency of electronic iteration to output
+                                       ///< descriptors and labels, default is 0, which means no output until convergence
     bool deepks_scf = false;           ///< (need libnpy and libtorch) if set to true, a trained model
                                        ///< would be needed to calculate V_delta and F_delta
     int deepks_bandgap = 0;       ///< for bandgap label. QO added 2021-12-15

source/module_parameter/system_parameter.h

@@ -44,6 +44,7 @@ struct System_para
     std::string global_stru_dir = "";       ///< global structure directory
     std::string global_matrix_dir = "";     ///< global matrix directory
     std::string global_mlkedf_descriptor_dir = ""; ///< global ML KEDF descriptor directory
+    std::string global_deepks_label_elec_dir = ""; ///< global directory for DeePKS labels during electronic steps 
     std::string log_file = "log";           ///< log file name
 
     bool deepks_setorb = false;             ///< true if "deepks" is set

source/source_base/global_file.cpp

@@ -194,6 +194,47 @@ void ModuleBase::Global_File::make_dir_out(
 #endif
     }
 
+    if(PARAM.inp.deepks_out_freq_elec > 0)
+    {
+        int make_dir_deepks_elec = 0;
+        std::string command1 =  "test -d " + PARAM.globalv.global_deepks_label_elec_dir + " || mkdir " + PARAM.globalv.global_deepks_label_elec_dir;
+
+        times = 0;
+        while(times<GlobalV::NPROC)
+        {
+            if(rank==times)
+            {
+                if ( system( command1.c_str() ) == 0 )
+                {
+                    std::cout << " MAKE THE DEEPKS LABELS (ELEC) DIR    : " << PARAM.globalv.global_deepks_label_elec_dir << std::endl;
+                    make_dir_deepks_elec = 1;
+                }
+                else
+                {
+                    std::cout << " PROC " << rank << " CAN NOT MAKE THE DEEPKS LABELS (ELEC) DIR !!! " << std::endl;
+                    make_dir_deepks_elec = 0;
+                }
+            }
+#ifdef __MPI
+            Parallel_Reduce::reduce_all(make_dir_deepks_elec);
+#endif
+            if(make_dir_deepks_elec > 0)
+            { 
+                break;
+            }
+            ++times;
+        }
+
+#ifdef __MPI
+        if(make_dir_deepks_elec == 0)
+        {
+            std::cout << " CAN NOT MAKE THE DEEPKS LABELS (ELEC) DIR......." << std::endl;
+            ModuleBase::QUIT();
+        }
+        MPI_Barrier(MPI_COMM_WORLD);
+#endif
+    }
+
     // mohan add 2010-09-12
     if(out_alllog)
     {

source/source_esolver/esolver_ks_lcao.cpp

@@ -795,7 +795,7 @@ void ESolver_KS_LCAO<TK, TR>::iter_finish(UnitCell& ucell, const int istep, int&
         GlobalC::dftu.output(ucell);
     }
 
-    // 2) for deepks, calculate delta_e
+    // 2) for deepks, calculate delta_e, output labels during electronic steps
 #ifdef __MLALGO
     if (PARAM.inp.deepks_scf)
     {
@@ -870,6 +870,37 @@ void ESolver_KS_LCAO<TK, TR>::iter_finish(UnitCell& ucell, const int istep, int&
     {
         GlobalC::dftu.initialed_locale = true;
     }
+
+    // 9) for deepks, output labels during electronic steps (after conv_esolver is renewed)
+#ifdef __MLALGO
+    if (PARAM.inp.deepks_out_labels >0 && PARAM.inp.deepks_out_freq_elec)
+    {
+        if (iter % PARAM.inp.deepks_out_freq_elec == 0 )
+        {
+            hamilt::HamiltLCAO<TK, TR>* p_ham_deepks = dynamic_cast<hamilt::HamiltLCAO<TK, TR>*>(this->p_hamilt);
+            std::shared_ptr<LCAO_Deepks<TK>> ld_shared_ptr(&ld, [](LCAO_Deepks<TK>*) {});
+            LCAO_Deepks_Interface<TK, TR> deepks_interface(ld_shared_ptr);
+
+            deepks_interface.out_deepks_labels(this->pelec->f_en.etot,
+                                               this->kv.get_nks(),
+                                               ucell.nat,
+                                               PARAM.globalv.nlocal,
+                                               this->pelec->ekb,
+                                               this->kv.kvec_d,
+                                               ucell,
+                                               orb_,
+                                               this->gd,
+                                               &(this->pv),
+                                               *(this->psi),
+                                               dynamic_cast<const elecstate::ElecStateLCAO<TK>*>(this->pelec)->get_DM(),
+                                               p_ham_deepks,
+                                               iter,
+                                               conv_esolver,
+                                               GlobalV::MY_RANK,
+                                               GlobalV::ofs_running);
+        }
+    }
+#endif
 }
 
 template class ESolver_KS_LCAO<double, double>;

tests/09_DeePKS/103_NO_GO_deepks_out_freq_elec/INPUT

@@ -0,0 +1,39 @@
+INPUT_PARAMETERS
+#Parameters (1.General)
+suffix			    autotest
+calculation     	scf
+
+nbands			    6
+symmetry		    0
+pseudo_dir		    ../../PP_ORB
+orbital_dir		    ../../PP_ORB
+
+#Parameters (2.Iteration)
+ecutwfc			    50
+scf_thr			    5e-3 #relatively large to avoid many steps
+scf_nmax            50
+
+#Parameters (3.Basis)
+basis_type		    lcao
+gamma_only		    1
+
+#Parameters (4.Smearing)
+smearing_method		gaussian
+smearing_sigma		0.02
+
+#Parameters (5.Mixing)
+mixing_type		    broyden
+mixing_beta		    0.8
+mixing_gg0          0
+
+#Parameters (6.File)
+deepks_out_labels   1 
+
+deepks_bandgap      1
+deepks_v_delta      1
+cal_force           1
+cal_stress          1
+
+deepks_out_freq_elec 1
+
+dft_functional      lda

tests/09_DeePKS/103_NO_GO_deepks_out_freq_elec/KPT

@@ -0,0 +1,4 @@
+K_POINTS
+0
+Gamma
+1 1 1 0 0 0

tests/09_DeePKS/103_NO_GO_deepks_out_freq_elec/README

@@ -0,0 +1 @@
+test output for deepks_out_freq_elec > 0, for gamma_only H2O molecule

tests/09_DeePKS/103_NO_GO_deepks_out_freq_elec/STRU

@@ -0,0 +1,32 @@
+ATOMIC_SPECIES
+H 1.008  H_ONCV_PBE2PZ-1.0.upf
+O 15.9994 O_ONCV_PBE2PZ-1.0.upf
+
+NUMERICAL_ORBITAL
+H_gga_6au_60Ry_2s1p.orb
+O_gga_6au_60Ry_2s2p1d.orb
+
+NUMERICAL_DESCRIPTOR
+../Model_ProjOrb/2au_20Ry_jle.orb
+
+LATTICE_CONSTANT
+1
+
+LATTICE_VECTORS
+15 0 0
+0 15 0
+0 0 15
+
+ATOMIC_POSITIONS
+Cartesian
+
+H
+0
+2
+-12.081531451316582 1.463368531712373 10.304287878967891 1 1 1
+-12.056180479123837 4.25408045699522 10.010554611214044 1 1 1
+O
+0
+1
+-13.1930046246741 2.91132430713516 10.440289103003526 1 1 1
+

tests/09_DeePKS/103_NO_GO_deepks_out_freq_elec/result.ref

@@ -0,0 +1,27 @@
+etotref -463.3717558791495890
+etotperatomref -154.4572519597
+totalforceref 3.661556
+totalstressref 9.690490
+deepks_e_label 17.02859183994638
+deepks_edelta 0.0
+deepks_o_label 0.3103148500892656
+deepks_odelta 0.0
+deepks_oprec -0.42236872134755643
+deepks_h_label 49.03922791708764
+deepks_vdelta 0.0
+deepks_vdp 176.26236130240005
+deepks_f_label 0.07120588316979731
+deepks_fdelta 0.0
+deepks_fpre 19.53717375096851
+deepks_s_label 0.08886477565384877
+deepks_sdelta 0.0
+deepks_spre 19.22629462223564
+deepks_e_label_elec 51.047596216236645
+deepks_edelta_elec 0
+deepks_o_label_elec .86018933752816427
+deepks_odelta_elec 0
+deepks_oprec_elec -1.23599151361440239
+deepks_h_label_elec 147.028760179109205
+deepks_vdelta_elec 0
+deepks_vdp_elec 529.07335498441735
+totaltimeref 1.03

tests/09_DeePKS/103_NO_KP_deepks_out_freq_elec/INPUT

@@ -0,0 +1,39 @@
+INPUT_PARAMETERS
+#Parameters (1.General)
+suffix			    autotest
+calculation     	scf
+
+nbands			    6
+symmetry		    0
+pseudo_dir		    ../../PP_ORB
+orbital_dir		    ../../PP_ORB
+
+#Parameters (2.Iteration)
+ecutwfc			    50
+scf_thr			    5e-3 #relatively large to avoid many steps
+scf_nmax            50
+
+#Parameters (3.Basis)
+basis_type		    lcao
+gamma_only		    0
+
+#Parameters (4.Smearing)
+smearing_method		gaussian
+smearing_sigma		0.02
+
+#Parameters (5.Mixing)
+mixing_type		    broyden
+mixing_beta		    0.8
+mixing_gg0          0
+
+#Parameters (6.File)
+deepks_out_labels   1 
+
+deepks_bandgap      1
+deepks_v_delta      1
+cal_force           1
+cal_stress          1
+
+deepks_out_freq_elec 1
+
+dft_functional      lda