Feature : CG diagonalization & mixing of dij in PAW

source/module_cell/module_paw/paw_cell.cpp

@@ -76,6 +76,29 @@ void Paw_Cell::init_paw_cell(
     }
 
     this -> init_rhoij();
+    //this -> init_mix_dij();
+}
+
+void Paw_Cell::init_mix_dij()
+{
+    first_iter = true;
+    count = 0;
+
+    if(GlobalV::RANK_IN_POOL == 0)
+    {
+        dij_save.resize(natom);
+        for(int iat = 0; iat < natom; iat ++)
+        {
+            const int it = atom_type[iat];
+            const int nproj = paw_element_list[it].get_mstates();
+            const int size_dij = nproj * (nproj+1) / 2;
+            dij_save[iat].resize(size_dij * nspden);
+            for(int i = 0; i < size_dij * nspden; i ++)
+            {
+                dij_save[iat][i] = 0.0;
+            }
+        }
+    }
 }
 
 void Paw_Cell::init_rhoij()

source/module_cell/module_paw/paw_cell.h

@@ -236,6 +236,15 @@ class Paw_Cell
     double calculate_ecore();
 
     private:
+
+    bool first_iter;
+    int count;
+
+    void init_mix_dij();
+    void mix_dij(const int iat, double*dij_libpaw);
+
+    std::vector<std::vector<double>> dij_save;
+
 // Info to be passed to libpaw_interface:
 // 1. ecut, ecutpaw : kinetic energy cutoff of the planewave basis set
 // there will be one coarse grid for density/potential, and a fine grid for PAW

source/module_cell/module_paw/paw_cell_libpaw.cpp

@@ -192,6 +192,28 @@ void Paw_Cell::set_libpaw_atom(const int natom_in, const int ntypat_in, const in
     }
 }
 
+void Paw_Cell::mix_dij(const int iat, double*dij_paw)
+{
+    double mixing_beta = 0.1;
+
+    std::cout << "mixing_beta : " << mixing_beta << std::endl;
+
+    const int it = atom_type[iat];
+    const int nproj = paw_element_list[it].get_mstates();
+    const int size_dij = nproj * (nproj+1) / 2;
+    for(int i = 0; i < size_dij * nspden; i ++)
+    {  
+        if(!first_iter) dij_paw[i] = dij_save[iat][i] * (1.0 - mixing_beta) + dij_paw[i] * mixing_beta;
+
+        if(count > 30) dij_paw[i] = dij_save[iat][i];
+
+        dij_save[iat][i] = dij_paw[i];
+    }
+
+    first_iter = false;
+    count ++;
+}
+
 // Sets filename_list
 // I'm going to read directly from STRU file
 void Paw_Cell::set_libpaw_files()
@@ -609,13 +631,16 @@ void Paw_Cell::set_dij()
            dij[is] = new double[nproj * nproj];
         }
 
-
-
 #ifdef __MPI
-        if(GlobalV::RANK_IN_POOL == 0) extract_dij(iat,size_dij,dij_libpaw);
+        if(GlobalV::RANK_IN_POOL == 0)
+        {
+            extract_dij(iat,size_dij,dij_libpaw);
+            //mix_dij(iat,dij_libpaw);
+        }
         Parallel_Common::bcast_double(dij_libpaw,size_dij*nspden);
 #else
         extract_dij(iat,size_dij,dij_libpaw);
+        //mix_dij(iat,dij_libpaw);
 #endif
 
         for(int is = 0; is < nspden; is ++)

source/module_elecstate/module_charge/charge_mixing.cpp

@@ -99,6 +99,7 @@ void Charge_Mixing::auto_set(const double& bandgap_in, const UnitCell& ucell_)
     {
         this->mixing->mixing_beta = this->mixing_beta = 0.7;
     }
+    GlobalV::MIXING_BETA = mixing_beta;
     GlobalV::ofs_running << "      Autoset mixing_beta to " << this->mixing_beta << std::endl;
 
     bool has_trans_metal = false;

source/module_hamilt_pw/hamilt_pwdft/hamilt_pw.cpp

@@ -10,6 +10,10 @@
 #include "operator_pw/meta_pw.h"
 #include "operator_pw/nonlocal_pw.h"
 
+#ifdef USE_PAW
+#include "module_cell/module_paw/paw_cell.h"
+#endif
+
 namespace hamilt
 {
 
@@ -207,6 +211,21 @@ void HamiltPW<T, Device>::sPsi(const T* psi_in, // psi
     const T one{1, 0};
     const T zero{0, 0};
 
+    if(GlobalV::use_paw)
+    {
+#ifdef USE_PAW
+#ifdef __DEBUG
+        assert(psi.get_k_first());
+#endif
+        for(int m = 0; m < nbands; m ++)
+        {
+            GlobalC::paw_cell.paw_nl_psi(1, reinterpret_cast<const std::complex<double>*> (&psi_in[m * npw]),
+                reinterpret_cast<std::complex<double>*>(&spsi[m * nrow]));
+        }
+#endif
+        return;
+    }
+
     syncmem_op()(this->ctx, this->ctx, spsi, psi_in, static_cast<size_t>(nbands * nrow));
     if (GlobalV::use_uspp)
     {

source/module_io/input.cpp

@@ -4026,7 +4026,7 @@ int Input::count_ntype(const std::string &fn)
         while (true)
         {
             ModuleBase::GlobalFunc::READ_VALUE(ifa, temp);
-            if (temp == "LATTICE_CONSTANT" || temp == "NUMERICAL_ORBITAL" || temp == "NUMERICAL_DESCRIPTOR"
+            if (temp == "LATTICE_CONSTANT" || temp == "NUMERICAL_ORBITAL" || temp == "NUMERICAL_DESCRIPTOR" || temp == "PAW_FILES"
                 || ifa.eof())
             {
                 break;

source/module_io/parameter_pool.cpp

@@ -46,7 +46,7 @@ int count_ntype(const std::string& fn)
         while (true)
         {
             ModuleBase::GlobalFunc::READ_VALUE(ifa, temp);
-            if (temp == "LATTICE_CONSTANT" || temp == "NUMERICAL_ORBITAL" || temp == "NUMERICAL_DESCRIPTOR"
+            if (temp == "LATTICE_CONSTANT" || temp == "NUMERICAL_ORBITAL" || temp == "NUMERICAL_DESCRIPTOR" || temp == "PAW_FILES"
                 || ifa.eof())
             {
                 break;