Sto dft

ABACUS.develop/see.txt

@@ -0,0 +1,17 @@
+diff -r source/src_pw/unitcell_pseudo.cpp /media/qianrui/work/LQR/test/abacus-develop/ABACUS.develop/source/src_pw/unitcell_pseudo.cpp
+1264,1265c1264
+< 	if ( CALCULATION!="scf-sto" && CALCULATION!="relax-sto" && CALCULATION!="md-sto" ) //qianrui 2021-2-20
+< 	{
+---
+> 	
+1291,1295c1290,1293
+< 		if(NBANDS < mag.get_neldw() )
+<         {
+<             WARNING_QUIT("unitcell","Too few spin down bands!");
+<         }
+< 	}
+---
+> 		if(NBANDS < mag.get_neldw() ) 
+> 		{
+> 			WARNING_QUIT("unitcell","Too few spin down bands!");
+> 		}

ABACUS.develop/source/Makefile.vars

@@ -28,7 +28,7 @@ ELPA_DIR = /home/mohan/1_Software/impi_elpa-16.05.005
 #ELPA_DIR = /opt/elpa/intel_2017_update4
 
 CEREAL_DIR = /home/mohan/1_Software/ABACUS_Github/cereal/cereal/
+#CEREAL_DIR = /home/qianrui/intelcompile/cereal/
 
 OBJ_DIR = obj
 NP      = 14
-

ABACUS.develop/source/input.cpp

@@ -2384,6 +2384,14 @@ void Input::Check(void)
 */
                 nstep = 1;
 
+    }
+	else if (calculation == "scf-sto")  // qianrui 2021-2-20
+    {
+                if(mem_saver == 1)
+                {
+                        mem_saver = 0;
+                        AUTO_SET("mem_savre","0");
+                }
     }
     else if (calculation == "relax")  // pengfei 2014-10-13
     {

ABACUS.develop/source/run_pw.cpp

@@ -75,7 +75,14 @@ void Run_pw::plane_wave_line(void)
     //=====================================
     CHR.allocate(NSPIN, pw.nrxx, pw.ngmc);
     pot.allocate(pw.nrxx);
-    wf.allocate(kv.nks);
+    if ( NBANDS != 0 || (CALCULATION!="scf-sto" && CALCULATION!="relax-sto" && CALCULATION!="md-sto") )//qianrui add 
+	{
+        wf.allocate(kv.nks);
+    }
+    else
+    {
+        wf.npwx = wf.setupIndGk(pw, kv.nks);
+    }
 	UFFT.allocate();
 
     //=======================
@@ -116,7 +123,12 @@ void Run_pw::plane_wave_line(void)
     //================================
     // Initial start wave functions
     //================================
-   	wf.wfcinit();
+    if ( NBANDS != 0 || (CALCULATION!="scf-sto" && CALCULATION!="relax-sto" && CALCULATION!="md-sto") )//qianrui add 
+	{
+    	wf.wfcinit();
+    }
+
+
 
 	switch(exx_global.info.hybrid_type)				// Peize Lin add 2019-03-09
 	{

ABACUS.develop/source/src_global/global_function.h

@@ -184,7 +184,7 @@ void SCAN_END(ifstream &ifs, const string &TargetName);
 template<class T>
 static inline void DCOPY( const T &a, T &b, const int &dim)
 {
-    for (int i=0; i<dim; i++) b[i] = a[i];
+    for (int i=0; i<dim; ++i) b[i] = a[i];
 }
 
 void BLOCK_HERE( const string &description );

ABACUS.develop/source/src_pw/energy.cpp

@@ -87,6 +87,18 @@ void energy::calculate_etot(void)
 	}			
 
     //Quxin adds for DFT+U energy correction on 20201029
+
+	/*cout << resetiosflags(ios::scientific) << endl;
+	cout << setprecision(6) << endl;
+	cout << " eband=" << eband << endl;
+	cout << " deband=" << deband << endl;
+	cout << " etxc-etxcc=" << etxc-etxcc << endl;
+	cout << " ewld=" << ewld << endl;
+	cout << " ehart=" << ehart << endl;
+	cout << " demet=" << demet << endl;
+	cout << " descf=" << descf << endl;
+	cout << " efiled=" << Efield::etotefield << endl;
+	cout << " fermienergy= "<<ef<<endl;*/
     if(INPUT.dft_plus_u) 
 	{
 		this->etot += dftu.EU;																	  

ABACUS.develop/source/src_pw/pseudopot_cell_vnl.h

@@ -29,6 +29,7 @@ class pseudopot_cell_vnl: public pseudopot_cell_vl
 	friend class Hamilt_PW;
 	friend class WF_atomic;
 	friend class wavefunc;
+	friend class Stochastic_hchi;
 
 	void init(const int ntype, const bool allocate_vkb=1);
 

ABACUS.develop/source/src_pw/sto_che.cpp

@@ -1,12 +1,155 @@
-#include "tools.h"
-#include "global.h"
 #include "sto_che.h" 
+#include "global.h"
 
 
 Stochastic_Chebychev::Stochastic_Chebychev()
 {
+    initplan = false;
+    initcoef = false;
+    getcoef = false;
+    getpolyval = false;
+    extend = 10;
+    norder = 5;
+    ccoef = new complex<double> [1];
+    coef = new double [1];
+    polyvalue = new complex<double> [1];
 }
 
 Stochastic_Chebychev::~Stochastic_Chebychev()
 {
+    if(initplan)
+    {
+        fftw_destroy_plan(plancoef);
+    }
+    delete [] coef;
+    delete [] ccoef;
+    delete [] polyvalue;
+}
+void Stochastic_Chebychev:: init()
+{
+    norder = STO_WF.nche_sto;
+    assert(norder > 5);
+    assert(extend >= 1);
+    if(norder != 0)
+    {
+        norder2 = 2 * norder * extend;
+        delete[] coef;
+        delete[] ccoef;
+        delete[] polyvalue;
+        ccoef = new complex<double> [norder2];
+        coef = new double [norder2];
+        polyvalue = new complex<double> [norder];
+        initcoef = true;
+    }
+    else
+    {
+        WARNING_QUIT("Stochastic_Chebychev", "The Chebychev expansion order should be at least one!");
+    }
+
+}
+
+void Stochastic_Chebychev:: calcoef(double fun(double))
+{
+    if(!initcoef) WARNING_QUIT("Stochastic_Chebychev", "Please init coef first!");
+    for(int i = 0; i < norder2; ++i)
+    {
+        coef[i]=fun(cos((i+0.5)*TWO_PI/norder2));
+    }
+     if(!initplan)
+    {
+        initplan = true;
+        plancoef = fftw_plan_dft_r2c_1d(norder2, coef, (fftw_complex *) ccoef, FFTW_MEASURE);
+    }
+    fftw_execute(plancoef);
+    complex<double> ui(0,1);
+    for(int i = 0; i<norder; ++i)
+    {
+        if(i == 0)
+        {
+            coef[i] = real(exp(ui*i*PI/norder2) * ccoef[i]) / norder2;
+        }
+        else
+        {
+            coef[i] = real(exp(ui*i*PI/norder2) * ccoef[i]) * 2 / norder2;
+        }
+    }
+    getcoef = true;
+}
+
+void Stochastic_Chebychev:: recurs(double&tnp1, double& tn, double& tn_1, double &t)
+{   
+    tnp1 = 2*t*tn-tn_1;
+}
+
+complex<double> Stochastic_Chebychev:: calresult()
+{
+    if(!getcoef||!getpolyval) WARNING_QUIT("Stochastic_Chebychev", "Please calculate coef or polyval first!");
+    complex<double> result = 0;
+    for(int ior = 0; ior < norder; ++ior)
+    {
+        result += coef[ior] * polyvalue[ior];
+    }
+    return result;   
+}
+
+void Stochastic_Chebychev:: calresult(double &t, double& result)
+{
+    if(!getcoef) WARNING_QUIT("Stochastic_Chebychev", "Please calculate coef first!");
+    double tnp1, tn, tn_1;
+    tn_1 = 1;
+    tn = tn_1 * t;
+    //0- & 1-st order
+    result = coef[0] * tn_1 + coef[1] * tn;
+
+    //more than 1-st orders
+    for(int ior = 2; ior < norder; ++ior)
+    {
+        recurs(tnp1, tn, tn_1, t);
+        result += coef[ior] * tnp1;
+        tn_1 = tn;
+        tn = tnp1; 
+    }
+    return;
+}
+
+
+void Stochastic_Chebychev:: calpolyval(void tfun(complex<double> *in, complex<double> *out), int& ndim, complex<double> *wavein)
+{
+
+    complex<double> *arraynp1, *arrayn, *arrayn_1;
+    arraynp1 = new complex<double> [ndim];
+    arrayn = new complex<double> [ndim];
+    arrayn_1 = new complex<double> [ndim];
+
+    DCOPY(wavein, arrayn_1, ndim);
+    tfun(arrayn_1, arrayn);
+
+    ZEROS(polyvalue,norder);
+    //0- & 1-st order
+    for(int i = 0; i < ndim; ++i) 
+    {
+        polyvalue[0] += conj(wavein[i]) * wavein[i];// 0-th order : <wavein | wavein> = ndim
+        polyvalue[1] += conj(wavein[i]) * arrayn[i];// 1-st order : <wavein | tfun | wavein>
+    }
+
+    //more than 1-st orders
+    for(int ior = 2; ior < norder; ++ior)
+    {
+        recurs(arraynp1, arrayn, arrayn_1, tfun, ndim);
+        for(int i = 0; i < ndim; ++i) // n-th order : <wavein | T_n(tfun) | wavein>
+        {
+            polyvalue[ior] += conj(wavein[i]) * arraynp1[i];
+        }
+
+        complex<double>* tem = arrayn_1;
+        arrayn_1 = arrayn;
+        arrayn = arraynp1;
+        arraynp1 = tem; 
+    }
+
+    delete [] arraynp1;
+    delete [] arrayn;
+    delete [] arrayn_1;
+    getpolyval = true;
+    return;
 }

ABACUS.develop/source/src_pw/sto_che.h

@@ -2,7 +2,7 @@
 #define INCLUDE_STO_CHEBYCHEV_H
 
 #include "tools.h"
-
+#include "../src_global/mymath.h"
 //----------------------------------------------
 // Chebychev Filtering
 //----------------------------------------------
@@ -15,12 +15,87 @@ class Stochastic_Chebychev
     // constructor and deconstructor
     Stochastic_Chebychev();
     ~Stochastic_Chebychev();
+    void init();
+
+    void calcoef(double fun(double));
+    complex<double> calresult();
+    void calresult(double &t, double &result);
+
+    template<class T>
+    void calresult(void fun(T *in, T *out), int& ndim, T *wavein, T *waveout);
+
+
+
+    void calpolyval(void fun(complex<double> *in, complex<double> *out), int& ndim, complex<double> *wavein);
+
     int norder;
+    int extend;
+    int norder2;  // 2 * norder
+    double* coef;  //[norder2] expansion coefficient of each order, only first norder coefficients are usefull
+    complex<double> *ccoef;  //[norder2] temporary complex expansion coefficient of each order, only first norder coefficients are usefull.
+    complex<double> *polyvalue; //
+    fftw_plan plancoef;
+    bool initplan, initcoef, getcoef, getpolyval;
 
 	private:
-
+    void recurs(double&tnp1, double &tn, double &tn_1, double& t); //tnp1: T_(n+1), tn: T_n, tn_1: T_(n-1)
+    template<class T>
+    void recurs(T *arraynp1, T* arrayn, T *arrayn_1, void fun(T *in,T *out), int& ndim);
 
 
 };
 
+template<class T>
+void Stochastic_Chebychev:: recurs(T *arraynp1, T* arrayn, T *arrayn_1, void fun(T *in,T *out), int& ndim)
+{
+    fun(arrayn,arraynp1);
+    for(int i = 0; i < ndim; ++i)
+    {
+        arraynp1[i]=2*arraynp1[i]-arrayn_1[i];
+    }
+}
+
+template<class T>
+void Stochastic_Chebychev:: calresult(void tfun(T *in, T *out), int &ndim, T *wavein, T *waveout)
+{
+    if(!getcoef) WARNING_QUIT("Stochastic_Chebychev", "Please calculate coef first!");
+
+    T *arraynp1, *arrayn, *arrayn_1;
+    arraynp1 = new T [ndim];
+    arrayn = new T [ndim];
+    arrayn_1 = new T [ndim];
+    DCOPY(wavein, arrayn_1, ndim);
+    tfun(arrayn_1, arrayn);
+
+    //0- & 1-st order
+    for(int i = 0; i < ndim; ++i)
+    {
+        waveout[i] = coef[0] * arrayn_1[i] + coef[1] * arrayn[i];
+    }
+
+
+
+    //more than 1-st orders
+    for(int ior = 2; ior < norder; ++ior)
+    {
+        recurs(arraynp1, arrayn, arrayn_1, tfun, ndim);
+        for(int i = 0; i < ndim; ++i)
+        {
+            waveout[i] += coef[ior] * arraynp1[i];
+        }
+        T * tem = arrayn_1;
+        arrayn_1 = arrayn;
+        arrayn = arraynp1;
+        arraynp1 = tem; 
+    }
+    delete [] arraynp1;
+    delete [] arrayn;
+    delete [] arrayn_1;
+    return;
+}
+
+
+
+
+
 #endif// Eelectrons_Chebychev