Fix: 1. non-collinear spin error, 2. davidson method error

Author: dyzheng

source/module_hamilt/hamilt_pw.cpp

@@ -26,7 +26,8 @@ HamiltPW::HamiltPW()
     {
         Operator* ekinetic = new Ekinetic<OperatorPW>( 
             tpiba2, 
-            gk2
+            gk2,
+            GlobalC::wfcpw->npwk_max
         );
         if(this->ops == nullptr)
         {

source/module_hamilt/ks_pw/ekinetic_pw.cpp

@@ -11,12 +11,14 @@ template class Ekinetic<OperatorPW>;
 template<>
 Ekinetic<OperatorPW>::Ekinetic(
     double tpiba2_in,
-    const double* gk2_in
+    const double* gk2_in,
+    const int gk2_dim_in
 )
 {
     this->cal_type = 10;
     this->tpiba2 = tpiba2_in;
     this->gk2 = gk2_in;
+    this->gk2_dim = gk2_dim_in;
     if( this->tpiba2 < 1e-10 || this->gk2 == nullptr
     )
     {
@@ -38,7 +40,7 @@ void Ekinetic<OperatorPW>::act
     this->max_npw = psi_in->get_nbasis() / psi_in->npol;
 
 
-    const double *gk2_ik = &(this->gk2[this->ik * this->max_npw]);
+    const double *gk2_ik = &(this->gk2[this->ik * this->gk2_dim]);
 
     for (int ib = 0; ib < n_npwx; ++ib)
     {

source/module_hamilt/ks_pw/ekinetic_pw.h

@@ -12,7 +12,8 @@ class Ekinetic : public T
     public:
     Ekinetic(
         double tpiba2_in, 
-        const double* gk2_in
+        const double* gk2_in,
+        const int gk2_dim_in
     );
 
     virtual void act
@@ -33,6 +34,8 @@ class Ekinetic : public T
 
     const double* gk2 = nullptr;
 
+    int gk2_dim = 0;
+
 };
 
 } // namespace hamilt

source/module_hamilt/ks_pw/meta_pw.cpp

@@ -49,7 +49,7 @@ void Meta<OperatorPW>::act
 
     const int npw = psi_in->get_ngk(this->ik);
     const int current_spin = this->isk[this->ik];
-    this->max_npw = psi_in->get_nbasis();
+    this->max_npw = psi_in->get_nbasis() / psi_in->npol;
     //npol == 2 case has not been considered
     this->npol = psi_in->npol;
 

source/module_hamilt/ks_pw/nonlocal_pw.cpp

@@ -75,7 +75,7 @@ void Nonlocal<OperatorPW>::act
                    &nkb,
                    &ModuleBase::ONE,
                    this->ppcell->vkb.c,
-                   &this->max_npw,
+                   &this->ppcell->vkb.nc,
                    tmpsi_in,
                    &inc,
                    &ModuleBase::ZERO,
@@ -92,7 +92,7 @@ void Nonlocal<OperatorPW>::act
                    &this->npw,
                    &ModuleBase::ONE,
                    this->ppcell->vkb.c,
-                   &this->max_npw,
+                   &this->ppcell->vkb.nc,
                    tmpsi_in,
                    &this->max_npw,
                    &ModuleBase::ZERO,
@@ -171,7 +171,7 @@ void Nonlocal<OperatorPW>::add_nonlocal_pp(std::complex<double> *hpsi_in, const
                     {
                         for (int ib = 0; ib < m; ib+=2)
                         {
-                            psind = (sum + ip2) * 2 * m + ib;
+                            psind = (sum + ip2) * m + ib;
                             becpind = ib * nkb + sum + ip;
                             becp1 = becp[becpind];
                             becp2 = becp[becpind + nkb];
@@ -201,7 +201,7 @@ void Nonlocal<OperatorPW>::add_nonlocal_pp(std::complex<double> *hpsi_in, const
                &(this->ppcell->nkb),
                &ModuleBase::ONE,
                this->ppcell->vkb.c,
-               &this->max_npw,
+               &this->ppcell->vkb.nc,
                ps,
                &inc,
                &ModuleBase::ONE,
@@ -218,7 +218,7 @@ void Nonlocal<OperatorPW>::add_nonlocal_pp(std::complex<double> *hpsi_in, const
                &(this->ppcell->nkb),
                &ModuleBase::ONE,
                this->ppcell->vkb.c,
-               &this->max_npw,
+               &this->ppcell->vkb.nc,
                ps,
                &npm,
                &ModuleBase::ONE,

source/module_hamilt/ks_pw/veff_pw.cpp

@@ -42,7 +42,7 @@ void Veff<OperatorPW>::act
     this->npol = psi_in->npol;
 
     std::complex<double> *porter = new std::complex<double>[wfcpw->nmaxgr];
-    for (int ib = 0; ib < n_npwx; ++ib)
+    for (int ib = 0; ib < n_npwx; ib += this->npol)
     {
         if (this->npol == 1)
         {

source/module_hsolver/diago_david.cpp

@@ -228,7 +228,7 @@ void DiagoDavid::cal_grad(hamilt::Hamilt* phm_in,
             }
         }
 
-        ppsi = &basis(m, 0);
+        ppsi = &basis(nbase + m, 0);
         spsi = &sp(nbase + m, 0);
         for (int ig = 0; ig < npw; ig++)
         {

source/module_psi/psi.h

@@ -53,9 +53,15 @@ class Psi
 {
 public:
     //Constructor 1: basic
-    Psi(void){};
+    Psi(void){
+        this->npol = GlobalV::NPOL;
+    };
     //Constructor 2: specify ngk only, should call resize() later
-    Psi(const int* ngk_in){this->ngk = ngk_in;}
+    Psi(const int* ngk_in)
+    {
+        this->ngk = ngk_in;
+        this->npol = GlobalV::NPOL;
+    }
     //Constructor 3: specify nk, nbands, nbasis, ngk, and do not need to call resize() later
     Psi(int nk_in, int nbd_in, int nbs_in, const int* ngk_in=nullptr)
     {
@@ -76,16 +82,20 @@ class Psi
         this->resize(nk_in, nband_in, psi_in.get_nbasis());
         this->ngk = psi_in.ngk;
         this->npol = psi_in.npol;
-        // copy from Psi from psi_in(current_k, 0, 0), 
-        // if size of k is 1, current_k in new Psi is psi_in.current_k 
-        const T* tmp = psi_in.get_pointer();
-        if(nk_in==1) for(size_t index=0; index<this->size();++index)
+
+        if(nband_in <= psi_in.get_nbands())
         {
-            psi[index] = tmp[index];
-            //current_k for this Psi only keep the spin index same as the copied Psi
-            this->current_k = psi_in.get_current_k();
-        } 
-        else for(size_t index=0; index<this->size();++index) psi[index] = tmp[index];
+            // copy from Psi from psi_in(current_k, 0, 0), 
+            // if size of k is 1, current_k in new Psi is psi_in.current_k 
+            const T* tmp = psi_in.get_pointer();
+            if(nk_in==1) for(size_t index=0; index<this->size();++index)
+            {
+                psi[index] = tmp[index];
+                //current_k for this Psi only keep the spin index same as the copied Psi
+                this->current_k = psi_in.get_current_k();
+            } 
+            else for(size_t index=0; index<this->size();++index) psi[index] = tmp[index];
+        }
     }
     // initialize the wavefunction coefficient
     // only resize and construct function now is used
@@ -235,12 +245,12 @@ class Psi
         else
         {
             const T* p = &this->psi[(range.index_1 * this->nbands + range.range_1) * this->nbasis];
-            int m = range.range_2 - range.range_1 + 1;
+            int m = (range.range_2 - range.range_1 + 1)* this->npol;
             return std::tuple<const T*, int>(p, m);
         }
     }
 
-    bool npol = 1;
+    int npol = 1;
 
  private:   
     std::vector<T> psi;