Merge pull request #913 from wenfei-li/develop

gint : removed more unused codes; unify the interfaces a little bit

Author: mohanchen

source/module_elecstate/elecstate_lcao.cpp

@@ -94,7 +94,7 @@ void ElecStateLCAO::psiToRho(const psi::Psi<std::complex<double>>& psi)
     //------------------------------------------------------------
 
     ModuleBase::GlobalFunc::NOTE("Calculate the charge on real space grid!");
-    this->uhm->GK.cal_rho_k(this->loc->DM_R);
+    this->uhm->GK.cal_rho_k(this->loc->DM_R, this->charge);
 
     this->charge->renormalize_rho();
 
@@ -142,7 +142,7 @@ void ElecStateLCAO::psiToRho(const psi::Psi<double>& psi)
     // calculate the charge density on real space grid.
     //------------------------------------------------------------
     ModuleBase::GlobalFunc::NOTE("Calculate the charge on real space grid!");
-    this->uhm->GG.cal_rho(this->loc->DM);
+    this->uhm->GG.cal_rho(this->loc->DM, this->charge);
 
     this->charge->renormalize_rho();
 

source/module_esolver/esolver_ks_lcao_elec.cpp

@@ -124,11 +124,11 @@ namespace ModuleESolver
             // calculate the charge density
             if (GlobalV::GAMMA_ONLY_LOCAL)
             {
-                this->UHM.GG.cal_rho(this->LOC.DM);
+                this->UHM.GG.cal_rho(this->LOC.DM, (Charge*)(&GlobalC::CHR));
             }
             else
             {
-                this->UHM.GK.cal_rho_k(this->LOC.DM_R);
+                this->UHM.GK.cal_rho_k(this->LOC.DM_R, (Charge*)(&GlobalC::CHR));
             }
 
             // renormalize the charge density

source/src_io/istate_charge.cpp

@@ -136,7 +136,7 @@ void IState_Charge::begin(Gint_Gamma &gg)
 
 			// (4) calculate charge density for a particular 
 			// band.
-   			gg.cal_rho(this->loc->DM);
+   			gg.cal_rho(this->loc->DM, (Charge*)(&GlobalC::CHR));
 			GlobalC::CHR.save_rho_before_sum_band(); //xiaohui add 2014-12-09
 			std::stringstream ss;
 			std::stringstream ss1;

source/src_lcao/CMakeLists.txt

@@ -36,6 +36,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/LCAO_hamilt.cpp

@@ -34,11 +34,7 @@ void LCAO_Hamilt::set_lcao_matrices(void)
     ModuleBase::timer::tick("LCAO_Hamilt","set_lcao_matrices");
 
     if(GlobalV::GAMMA_ONLY_LOCAL)
-    {
-        // mohan add 2012-03-29
-        // calculate the grid integration of 'Vl' matrix for gamma algorithms.
-        this->GG.prepare(GlobalC::ucell.latvec, GlobalC::ucell.lat0);
-    
+    {   
         // calulate the 'S', 'T' and 'Vnl' matrix for gamma algorithms.
         this->calculate_STNR_gamma();	
 
@@ -50,7 +46,6 @@ void LCAO_Hamilt::set_lcao_matrices(void)
 
         // calculate the grid integration of 'Vl' matrix for l-points algorithms.
         this->GK.init(GlobalC::pw.nbx, GlobalC::pw.nby, GlobalC::pw.nbzp, GlobalC::pw.nbzp_start, GlobalC::pw.ncxyz, this->LM);
-
     }
 
     // initial the overlap matrix is done.	

source/src_lcao/gint_gamma.cpp

@@ -5,16 +5,7 @@
 
 Gint_Gamma::Gint_Gamma()
 {
-    ylm1 = new double[100];
-    ylm2 = new double[100]; // can used for L=9
-    iq = new int[1];
-    x0 = new double[1];
-    x1 = new double[1];
-    x2 = new double[1];
-    x3 = new double[1];
-    x12 = new double[1];
-    x03 = new double[1];
-    
+   
     sender_index_size = 1;
 	sender_local_index = new int[1];
     sender_size_process = new int[1];
@@ -32,15 +23,6 @@ Gint_Gamma::Gint_Gamma()
 
 Gint_Gamma::~Gint_Gamma()
 {
-    delete[] ylm1;
-    delete[] ylm2;
-    delete[] iq;
-    delete[] x0;
-    delete[] x1;
-    delete[] x2;
-    delete[] x3;
-    delete[] x12;
-    delete[] x03;
 
     delete[] sender_local_index;
     delete[] sender_size_process;
@@ -51,61 +33,4 @@ Gint_Gamma::~Gint_Gamma()
     delete[] receiver_size_process;
     delete[] receiver_displacement_process;
     delete[] receiver_buffer;
-}
-
-
-void Gint_Gamma::save_atoms_on_grid(const Grid_Technique &gt)
-{
-    ModuleBase::TITLE("Grid_Integral","save_atoms_on_grid");
-
-    // mohan change.
-    max_size = gt.max_atom;
-
-	if(max_size == 0)
-	{
-		// mohan add return 2011-03-15
-		GlobalV::ofs_warning << " processor " << GlobalV::MY_RANK << ": no atom on sub-fft-grid." << std::endl;
-		return;
-	}
-
-    delete[] iq;
-    delete[] x0;
-    delete[] x1;
-    delete[] x2;
-    delete[] x3;
-    delete[] x12;
-    delete[] x03;
-    this->iq = new int[max_size];
-    this->x0 = new double[max_size];
-    this->x1 = new double[max_size];
-    this->x2 = new double[max_size];
-    this->x3 = new double[max_size];
-    this->x12 = new double[max_size];
-    this->x03 = new double[max_size];
-
-	ModuleBase::GlobalFunc::ZEROS(iq, max_size);
-	ModuleBase::GlobalFunc::ZEROS(x0, max_size);
-	ModuleBase::GlobalFunc::ZEROS(x1, max_size);
-	ModuleBase::GlobalFunc::ZEROS(x2, max_size);
-	ModuleBase::GlobalFunc::ZEROS(x3, max_size);
-	ModuleBase::GlobalFunc::ZEROS(x12, max_size);
-	ModuleBase::GlobalFunc::ZEROS(x03, max_size);
-
-	this->vfactor = std::abs(this->latvec0.Det())/gt.ncxyz;
-
-    return;
-}
-
-void Gint_Gamma::prepare(
-    const ModuleBase::Matrix3 &latvec_in,
-    const double& lat0_in)
-{
-	ModuleBase::TITLE("Grid_Base_Beta","prepare");
-
-	this->lat0 = lat0_in;
-
-	this->latvec0 = latvec_in;
-	this->latvec0 *= this->lat0;
-	
-	return;
 }

source/src_lcao/gint_gamma.h

@@ -11,6 +11,7 @@
 #include "../module_base/global_variable.h"
 #include "grid_technique.h"
 #include "LCAO_matrix.h"
+#include "../src_pw/charge.h"
 #include <omp.h>
 
 //=========================================================
@@ -30,7 +31,7 @@ class Gint_Gamma
 	void cal_vlocal( const double*const vlocal, LCAO_Matrix &lm);
 
 	// (2) calculate charge density
-	double cal_rho(double*** DM_in);
+	void cal_rho(double*** DM_in, Charge* chr);
 
 	// (3) calcualte the forces related to grid
 	void cal_force(double*** DM_in, const double*const vlocal,
@@ -42,36 +43,13 @@ class Gint_Gamma
 
 	// (5) calculate the Mulliken charge
 	void cal_mulliken(double** mulliken);
-	
-	void prepare( 
-		const ModuleBase::Matrix3 &latvec_in, 
-        const double& lat0_in);
 
 private:
 
-	double vfactor;
-	ModuleBase::Matrix3 latvec0;
-	double lat0;
     double***  DM;   //pointer to LOC.DM
-    double* transformer;
-	double psiv1;
-	double psiv2;
-	double* ylm1;
-	double* ylm2;
 
 	int grid_index;			// may delete?
 	int max_size;
-	
-	// these parameters are for interpolation.
-	// we store these parameters at first to speed
-	// up the calculation.
-	double *x0;
-	double *x1;
-	double *x2;
-	double *x3;
-	double* x12;
-	double* x03;
-	int *iq;
 
     ///===============================
     /// Use MPI_Alltoallv to convert a grid distributed matrix
@@ -90,32 +68,17 @@ class Gint_Gamma
     int *receiver_size_process;
     int *receiver_displacement_process;
     double* receiver_buffer;
-    
-    void save_atoms_on_grid(const Grid_Technique& gt);
 
 	// for calculation of < phi_i | Vlocal | phi_j >
 	// Input:	vlocal[ir]
 	// Output:	GridVlocal.ptr_2D[iw1_lo][iw2_lo]
 	Gint_Tools::Array_Pool<double> gamma_vlocal(const double*const vlocal) const;  
-
-	// for calculation of charege 
-	// Input:	DM[is][iw1_lo][iw2_lo]
-	// Output:	rho.ptr_2D[is][ir]
-	Gint_Tools::Array_Pool<double> gamma_charge(const double*const*const*const DM) const;
-
 	// for calculation of Mulliken charge.
 	void gamma_mulliken(double** mulliken);
 
 	// for calculation of envelope functions.
 	void gamma_envelope(const double* wfc, double* rho);// mohan add 2011-07-01
 
-
-	// for calculatin of < dphi_i | Vlocal | phi_j > for foce calculation
-	// on regular FFT real space grid.
-	void gamma_force(const double*const*const*const DM, const double*const vlocal,
-			ModuleBase::matrix& force, ModuleBase::matrix& stress,
-			const bool is_force, const bool is_stress);
-
 	void cal_meshball_vlocal(
 		const int na_grid,  						// how many atoms on this (i,j,k) grid
 		const int LD_pool,
@@ -139,7 +102,7 @@ class Gint_Gamma
 		const double*const*const psir_ylm,				// psir_ylm[GlobalC::pw.bxyz][LD_pool]
 		const int*const vindex,							// vindex[GlobalC::pw.bxyz]
 		const double*const*const*const DM,				// DM[GlobalV::NSPIN][lgd_now][lgd_now]
-		Gint_Tools::Array_Pool<double> &rho) const;		// rho[GlobalV::NSPIN][GlobalC::pw.nrxx]
+		Charge* chr) const;		// rho[GlobalV::NSPIN][GlobalC::pw.nrxx]
 
 	void cal_meshball_force(
 		const int grid_index,

source/src_lcao/gint_gamma_env.cpp

@@ -10,7 +10,7 @@ void Gint_Gamma::cal_env(const double* wfc, double* rho)
     ModuleBase::TITLE("Gint_Gamma","cal_env");
     ModuleBase::timer::tick("Gint_Gamma","cal_env");
 
-    this->save_atoms_on_grid(GlobalC::GridT);
+	this->max_size = GlobalC::GridT.max_atom;
     this->gamma_envelope(wfc, rho);
 
     ModuleBase::timer::tick("Gint_Gamma","cal_env");
@@ -141,15 +141,7 @@ void Gint_Gamma::gamma_envelope(const double* wfc, double* rho)
 						// we can get the parameters we need to do interpolation and
 						// store them first!! these can save a lot of effort.
 						const double position = distance[ib][id] / delta_r;
-						/*
-						   this->iq[id] = static_cast<int>(position);
-						   this->x0[id] = position - static_cast<double>(iq[id]);
-						   this->x1[id] = 1.0 - x0[id];
-						   this->x2[id] = 2.0 - x0[id];
-						   this->x3[id] = 3.0 - x0[id];
-						   this->x12[id] = x1[id]*x2[id] / 6.0;
-						   this->x03[id] = x0[id]*x3[id] / 2.0;
-						 */
+
 						int ip;
 						double dx, dx2, dx3;
 						double c1, c2, c3, c4;