gint_k : add interface for new force implementation (not done yet)

Author: wenfei-li

source/src_lcao/CMakeLists.txt

@@ -37,6 +37,7 @@ list(APPEND objects
     gint_gamma_vl.cpp
     gint_k.cpp
     gint_k_fvl.cpp
+    gint_k_fvl_new.cpp
     gint_k_init.cpp
     gint_k_rho.cpp
     gint_k_vl.cpp

source/src_lcao/FORCE_k.cpp

@@ -70,6 +70,7 @@ void Force_LCAO_k::ftable_k (
     // doing on the real space grid.
     // ---------------------------------------
     this->cal_fvl_dphi_k(dm2d, isforce, isstress, fvl_dphi, svl_dphi);
+	//this->cal_fvl_dphi_k_new(isforce, isstress, fvl_dphi, svl_dphi);
 
     this->calFvnlDbeta(dm2d, isforce, isstress, fvnl_dbeta, svnl_dbeta, GlobalV::vnl_method);
 
@@ -1152,6 +1153,50 @@ void Force_LCAO_k::cal_fvl_dphi_k(
 	return;
 }
 
+// calculate the force due to < phi | Vlocal | dphi >
+void Force_LCAO_k::cal_fvl_dphi_k_new(
+	const bool isforce, 
+    const bool isstress,
+    ModuleBase::matrix& fvl_dphi,
+	ModuleBase::matrix& svl_dphi)
+{
+	ModuleBase::TITLE("Force_LCAO_k","cal_fvl_dphi_k_new");
+	ModuleBase::timer::tick("Force_LCAO_k","cal_fvl_dphi_k_new");
+
+	if(!isforce&&!isstress) 
+	{
+		ModuleBase::timer::tick("Force_LCAO_k","cal_fvl_dphi_k_new");
+		return;
+	}
+
+	int istep = 1;
+
+	// if Vna potential is not used.
+	GlobalC::pot.init_pot(istep, GlobalC::pw.strucFac);
+
+
+	for(int is=0; is<GlobalV::NSPIN; ++is)
+	{
+		GlobalV::CURRENT_SPIN = is;
+		for(int ir=0; ir<GlobalC::pw.nrxx; ir++)
+		{
+			GlobalC::pot.vr_eff1[ir] = GlobalC::pot.vr_eff(GlobalV::CURRENT_SPIN, ir);
+		}
+
+		//--------------------------------
+		// Grid integration here.
+		//--------------------------------
+		// fvl_dphi can not be set to zero here if Vna is used
+		if(isstress||isforce) 
+		{
+			this->UHM->GK.cal_force_k(isforce, isstress, fvl_dphi,svl_dphi,GlobalC::pot.vr_eff1);
+		}
+	}
+
+	ModuleBase::timer::tick("Force_LCAO_k","cal_fvl_dphi_k_new");
+	return;
+}
+
 void Force_LCAO_k::calFvnlDbeta(
 	double** dm2d, 
 	const bool &isforce, 

source/src_lcao/FORCE_k.h

@@ -60,6 +60,7 @@ class Force_LCAO_k : public Force_LCAO_gamma
 
 	// calculate the force due to < phi | Vlocal | dphi >
 	void cal_fvl_dphi_k(double** dm2d, const bool isforce, const bool isstress, ModuleBase::matrix& fvl_dphi, ModuleBase::matrix& svl_dphi);
+	void cal_fvl_dphi_k_new(const bool isforce, const bool isstress, ModuleBase::matrix& fvl_dphi, ModuleBase::matrix& svl_dphi);
 
 	// old method to calculate the force due to < phi | dbeta > < beta | phi >
 	void cal_fvnl_dbeta_k(double** dm2d, const bool isforce, const bool isstress, ModuleBase::matrix& fvnl_dbeta, ModuleBase::matrix& svnl_dbeta);

source/src_lcao/gint_k.h

@@ -109,6 +109,13 @@ class Gint_k : public Gint_k_init
         //mohan add 2011-06-19
         //zhengdy add 2016-10-18
 
+    void cal_force_k(
+        const bool isforce,
+        const bool isstress,
+        ModuleBase::matrix& fvl_dphi, 
+        ModuleBase::matrix& svl_dphi, 
+        const double* vl);
+
     private:
 
     //------------------------------------------------------

source/src_lcao/gint_k_fvl_new.cpp

@@ -0,0 +1,110 @@
+#include "gint_k.h"
+#include "../src_pw/global.h"
+#include "global_fp.h" // mohan add 2021-01-30
+
+#include "../module_base/ylm.h"
+#include "../module_base/timer.h"
+
+void Gint_k::cal_force_k(
+	const bool isforce,
+	const bool isstress,
+	ModuleBase::matrix& fvl_dphi, 
+	ModuleBase::matrix& svl_dphi, 
+	const double *vl)
+{
+	ModuleBase::TITLE("Gint_k","cal_force_k");
+	ModuleBase::timer::tick("Gint_k","cal_force_k");
+
+    const int max_size = GlobalC::GridT.max_atom;
+
+	if(max_size)
+    {
+        const int nbx = GlobalC::GridT.nbx;
+        const int nby = GlobalC::GridT.nby;
+        const int nbz_start = GlobalC::GridT.nbzp_start;
+        const int nbz = GlobalC::GridT.nbzp;
+        const double dv = GlobalC::ucell.omega/GlobalC::pw.ncxyz;
+    
+        const int ncyz = GlobalC::pw.ncy*GlobalC::pw.nczp; // mohan add 2012-03-25
+
+        for (int i=0; i<nbx; i++)
+        {
+            const int ibx = i*GlobalC::pw.bx;
+            for (int j=0; j<nby; j++)
+            {
+                const int jby = j*GlobalC::pw.by;
+                for (int k=nbz_start; k<nbz_start+nbz; k++)
+                {
+                    const int kbz = k*GlobalC::pw.bz-GlobalC::pw.nczp_start;
+    
+                    const int grid_index = (k-nbz_start) + j * nbz + i * nby * nbz;
+    
+                    // get the value: how many atoms has orbital value on this grid.
+                    const int na_grid = GlobalC::GridT.how_many_atoms[ grid_index ];
+                    if(na_grid==0) continue;
+
+                    // it's a uniform grid to save orbital values, so the delta_r is a constant.
+                    const double delta_r = GlobalC::ORB.dr_uniform;						
+                    
+                    // here vindex refers to local potentials
+                    int* vindex = Gint_Tools::get_vindex(ncyz, ibx, jby, kbz);	
+                    
+                    int * block_iw, * block_index, * block_size;
+                    Gint_Tools::get_block_info(na_grid, grid_index, block_iw, block_index, block_size);
+
+                    //------------------------------------------------------
+                    // whether the atom-grid distance is larger than cutoff
+                    //------------------------------------------------------
+                    bool **cal_flag = Gint_Tools::get_cal_flag(na_grid, grid_index);
+
+                    // set up band matrix psir_ylm and psir_DM
+                    const int LD_pool = max_size*GlobalC::ucell.nwmax;
+
+                    Gint_Tools::Array_Pool<double> psir_ylm(GlobalC::pw.bxyz, LD_pool);
+                    Gint_Tools::Array_Pool<double> dpsir_ylm_x(GlobalC::pw.bxyz, LD_pool);
+                    Gint_Tools::Array_Pool<double> dpsir_ylm_y(GlobalC::pw.bxyz, LD_pool);
+                    Gint_Tools::Array_Pool<double> dpsir_ylm_z(GlobalC::pw.bxyz, LD_pool);
+
+                    Gint_Tools::cal_dpsir_ylm(
+                        na_grid, grid_index, delta_r,
+                        block_index, block_size, 
+                        cal_flag,
+                        psir_ylm.ptr_2D,
+                        dpsir_ylm_x.ptr_2D,
+                        dpsir_ylm_y.ptr_2D,
+                        dpsir_ylm_z.ptr_2D
+                    );
+
+                    double *vldr3 = Gint_Tools::get_vldr3(vl, ncyz, ibx, jby, kbz, dv);
+                    const Gint_Tools::Array_Pool<double> psir_vlbr3    = Gint_Tools::get_psir_vlbr3(
+						na_grid, LD_pool,
+						block_index, cal_flag,
+						vldr3, psir_ylm.ptr_2D);
+                    
+					Gint_Tools::Array_Pool<double> psir_vlbr3_DM(GlobalC::pw.bxyz, LD_pool);
+                    if(isforce)
+                    {
+                    }
+                    if(isstress)
+                    {
+                    }
+
+                    free(vldr3);		vldr3=nullptr;
+                    delete[] block_iw;
+                    delete[] block_index;
+                    delete[] block_size;
+
+                    for(int ib=0; ib<GlobalC::pw.bxyz; ++ib)
+                        free(cal_flag[ib]);
+                    free(cal_flag);			cal_flag=nullptr;
+                }//k
+            }//j
+        }//i
+    }//max_size
+
+
+	ModuleBase::timer::tick("Gint_k","cal_force_k");
+	return;
+}
+
+