Merge branch 'develop' of github.com:deepmodeling/abacus-develop into develop

Author: pxlxingliang

source/module_elecstate/test/updaterhok_pw_test.cpp

@@ -299,15 +299,15 @@ TEST_F(EState,RhoPW)
     GlobalC::CHR.allocate(GlobalV::NSPIN, GlobalC::pw.nrxx, GlobalC::pw.ngmc);
     // GlobalC::pot.allocate(GlobalC::pw.nrxx);
     // we need to supply NBANDS here
-    GlobalC::wf.allocate(GlobalC::kv.nks);
+    psi::Psi<std::complex<double>>* psi = GlobalC::wf.allocate(GlobalC::kv.nks);
     // std::cout<<"npwx "<<GlobalC::wf.npwx<<std::endl;
     GlobalC::UFFT.allocate();
 
     //====== read wavefunction ==========================================
     std::stringstream ssw;
     ssw <<GlobalV::global_out_dir<< "WAVEFUNC";
     // we need to supply out_wfc_pw here
-    read_wfc2(ssw.str(), GlobalC::wf.psi[0], GlobalC::pw.gcar);
+    read_wfc2(ssw.str(), psi[0], GlobalC::pw.gcar);
 
     // copy data from old wf.evc to new evc(an object of Psi)
     evc.resize(GlobalC::kv.nks,GlobalV::NBANDS,GlobalC::wf.npwx);
@@ -316,9 +316,10 @@ TEST_F(EState,RhoPW)
     {
 	    for(int j=0;j<GlobalC::wf.npwx;j++)
             {
-		    evc(i,j)=GlobalC::wf.psi[0](i,j);
+		    evc(i,j) = psi[0](i,j);
 	    }
     }
+    delete psi;
     // using class ElecStatePW to calculate rho
     elecstate::MockElecStatePW* kk;
     kk = new elecstate::MockElecStatePW(&GlobalC::pw,&GlobalC::CHR,GlobalV::NBANDS);

source/module_elecstate/test/updaterhok_pw_test.h

@@ -132,17 +132,17 @@ void Use_FFT::allocate()
     return;
 }
 
-void wavefunc::allocate(const int nks)
+psi::Psi<complex<double>>* wavefunc::allocate(const int nks)
 {
 	this->npwx = GlobalC::pw.setupIndGk(this->igk, GlobalC::kv.ngk);
 	this->wg.create(nks,GlobalV::NBANDS);
 	this->ekb = new double*[nks];
-	this->psi = new psi::Psi<std::complex<double>>(nks, GlobalV::NBANDS,npwx, nullptr);
+	psi::Psi<std::complex<double>>* psi = new psi::Psi<std::complex<double>>(nks, GlobalV::NBANDS,npwx, nullptr);
 	for (int ik=0;ik<nks;ik++)
 	{
 		this->ekb[ik] = new double[GlobalV::NBANDS];
 	}
-	return;
+	return psi;
 }
 
 bool Charge::read_rho(const int &is, const std::string &fn, double* rho) //add by dwan

source/module_esolver/esolver.h

@@ -5,6 +5,8 @@
 #include "../module_cell/unitcell_pseudo.h"
 #include "../src_pw/energy.h"
 #include "../module_base/matrix.h"
+//--------------temporary----------------------------
+#include "module_psi/psi.h"
 
 namespace ModuleESolver
 {
@@ -16,7 +18,15 @@ namespace ModuleESolver
         ESolver() {
             classname = "ESolver";
         }
-        virtual ~ESolver() {};
+        
+        virtual ~ESolver() 
+        {
+            //--------------temporary----------------------------
+            if(this->psi != nullptr)
+            {
+                delete psi;
+            }
+        }
 
         //virtual void Init(Input_EnSolver &inp, matrix &lattice_v)=0
         virtual void Init(Input& inp, UnitCell_pseudo& cell) = 0;
@@ -44,6 +54,13 @@ namespace ModuleESolver
         //get iterstep used in current scf
         virtual int getniter() { return 0; }
         string classname;
+
+        //--------------temporary----------------------------
+        // this is the interface of non-self-consistant calculation
+        virtual void nscf() {};
+        
+        //wavefunction coefficients
+        psi::Psi<std::complex<double>>* psi=nullptr;
     };
 
     void init_esolver(ESolver*& p_esolver, const string use_esol);

source/module_esolver/esolver_ks.h

@@ -28,7 +28,11 @@ namespace ModuleESolver
 
         // calculate electron density from a specific Hamiltonian
         virtual void hamilt2density(const int istep, const int iter, const double ethr);
-        // get
+
+        // calculate electron states from a specific Hamiltonian
+        virtual void hamilt2estates(const double ethr){};
+
+        // get current step of Ionic simulation
         virtual int getniter() override;
 
     protected:

source/module_esolver/esolver_ks_lcao_elec.cpp

@@ -26,6 +26,10 @@
 #endif
 #include "../src_pw/H_Ewald_pw.h"
 
+#include "module_hamilt/hamilt_lcao.h"
+#include "module_elecstate/elecstate_lcao.h"
+#include "module_hsolver/hsolver_lcao.h"
+
 namespace ModuleESolver
 {
 
@@ -188,6 +192,75 @@ namespace ModuleESolver
             srho.begin(is, GlobalC::CHR, GlobalC::pw, GlobalC::Pgrid, GlobalC::symm);
         }
 
+        //init Psi, HSolver, ElecState, Hamilt
+        if(this->phsol != nullptr)
+        {
+            if(this->phsol->classname != "HSolverLCAO")
+            {
+                delete this->phsol;
+                this->phsol = nullptr;
+            }
+        }
+        else
+        {
+            this->phsol = new hsolver::HSolverLCAO();
+            this->phsol->method = GlobalV::KS_SOLVER;
+        }
+        if(this->pelec != nullptr)
+        {
+            if(this->pelec->classname != "ElecStateLCAO")
+            {
+                delete this->pelec;
+                this->pelec = nullptr;
+            }
+        }
+        else
+        {
+            this->pelec = new elecstate::ElecStateLCAO(
+                    (Charge*)(&(GlobalC::CHR)),
+                    &(GlobalC::kv), 
+                    GlobalC::kv.nks,
+                    GlobalV::NBANDS,
+                    &(this->LOC),
+                    &(this->UHM),
+                    &(this->LOWF));
+        }
+        if(this->phami != nullptr)
+        {
+            if(this->phami->classname != "HamiltLCAO")
+            {
+                delete this->phami;
+                this->phami = nullptr;
+            }
+        }
+        else
+        {
+            // three cases for hamilt class
+            if(GlobalV::GAMMA_ONLY_LOCAL)
+            {
+                this->phami = new hamilt::HamiltLCAO<double, double>(
+                    &(this->UHM.GG),
+                    &(this->UHM.genH),
+                    &(this->LM) );
+            }
+            // non-collinear spin case would not use the second template now, 
+            // would add this feature in the future
+            /*else if(GlobalV::NSPIN==4)
+            {
+                this->phami = new hamilt::HamiltLCAO<std::complex<double>, std::complex<double>>(
+                    &(this->UHM.GK),
+                    &(this->UHM.genH),
+                    &(this->LM) );
+            }*/ 
+            else
+            {
+                this->phami = new hamilt::HamiltLCAO<std::complex<double>, double>(
+                    &(this->UHM.GK),
+                    &(this->UHM.genH),
+                    &(this->LM) );
+            }
+        }
+
         ModuleBase::timer::tick("ESolver_KS_LCAO", "beforescf");
         return;
     }