Merge https://github.com/ouqi0711/abacus-develop into develop

Author: ouqi0711

source/src_lcao/CMakeLists.txt

@@ -12,6 +12,7 @@ list(APPEND objects
     FORCE_gamma_edm.cpp
     FORCE_gamma_tvnl.cpp
     FORCE_gamma_vl.cpp
+    FORCE_gamma_vl_new.cpp
     FORCE_k.cpp
     LCAO_diago.cpp
     LCAO_evolve.cpp
@@ -33,6 +34,7 @@ list(APPEND objects
     gint_gamma_common.cpp
     gint_gamma_env.cpp
     gint_gamma_fvl.cpp
+    gint_gamma_fvl_new.cpp
     gint_gamma_mull.cpp
     gint_gamma_rho.cpp
     gint_gamma_vl.cpp

source/src_lcao/FORCE_gamma.cpp

@@ -52,6 +52,8 @@ void Force_LCAO_gamma::ftable_gamma (
     this->calFvnlDbeta(loc.dm_gamma, isforce, isstress, fvnl_dbeta, svnl_dbeta, GlobalV::vnl_method);
     this->cal_fvl_dphi(loc.dm_gamma, isforce, isstress, fvl_dphi, svl_dphi);
 
+    //this->cal_fvl_dphi_new(loc.DM, isforce, isstress, fvl_dphi, svl_dphi);
+  
     //caoyu add for DeePKS
 #ifdef __DEEPKS
     if (GlobalV::deepks_scf)

source/src_lcao/FORCE_gamma.h

@@ -128,6 +128,13 @@ class Force_LCAO_gamma
         ModuleBase::matrix& fvl_dphi,
 		ModuleBase::matrix& svl_dphi);
 
+	void cal_fvl_dphi_new(
+		double*** DM_in,
+		const bool isforce, 
+        const bool isstress,
+        ModuleBase::matrix& fvl_dphi,
+		ModuleBase::matrix& svl_dphi);
+
 	void calFvnlDbeta(
 		const std::vector<ModuleBase::matrix> &dm2d, 
 		const bool &isforce, 

source/src_lcao/FORCE_gamma_vl.cpp

@@ -40,9 +40,6 @@ void Force_LCAO_gamma::cal_fvl_dphi(
         tmpDHx[i] = this->UHM->LM->DHloc_fixed_x[i];
         tmpDHy[i] = this->UHM->LM->DHloc_fixed_y[i];
         tmpDHz[i] = this->UHM->LM->DHloc_fixed_z[i];
-        //std::cout << "  this->UHM->LM->DHloc_fixed_x=" <<  this->UHM->LM->DHloc_fixed_x[i] << std::endl;
-        //std::cout << "  this->UHM->LM->DHloc_fixed_y=" <<  this->UHM->LM->DHloc_fixed_y[i] << std::endl;
-        //std::cout << "  this->UHM->LM->DHloc_fixed_z=" <<  this->UHM->LM->DHloc_fixed_z[i] << std::endl;
     }
 
     //calculate dVL
@@ -65,13 +62,8 @@ void Force_LCAO_gamma::cal_fvl_dphi(
             GlobalC::pot.vr_eff1[ir] = GlobalC::pot.vr_eff(GlobalV::CURRENT_SPIN, ir);
         }
 
-        //  should not be set zero if VNA is used.
-        //  ModuleBase::GlobalFunc::ZEROS(this->UHM->LM->DHloc_fixed_x,this->ParaV->nloc);
-        //  ModuleBase::GlobalFunc::ZEROS(this->UHM->LM->DHloc_fixed_y,this->ParaV->nloc);
-        //  ModuleBase::GlobalFunc::ZEROS(this->UHM->LM->DHloc_fixed_z,this->ParaV->nloc);
         this->UHM->GG.cal_force(GlobalC::pot.vr_eff1);
 
-
         for(int i=0; i<GlobalV::NLOCAL; i++)
         {
             const int iat = GlobalC::ucell.iwt2iat[i];
@@ -102,26 +94,10 @@ void Force_LCAO_gamma::cal_fvl_dphi(
                         svl_dphi(1,2) += dm2d2 * this->UHM->LM->DHloc_fixed_23[index];
                         svl_dphi(2,2) += dm2d2 * this->UHM->LM->DHloc_fixed_33[index];
                     }
-                    //  std::cout << std::setw(5) << iat << std::setw(5) << iat2 
-                    //  << std::setw(5) << mu << std::setw(5) << nu
-                    //  << std::setw(15) << this->UHM->LM->DHloc_fixed_z[index] << std::endl;
                 }
             }
         }
-
-//          std::cout << "fvl_dphi:" << std::endl;
-//          for(int iat=0; iat<GlobalC::ucell.nat; ++iat)
-//          {
-//              std::cout << std::setw(5) << iat << std::setw(15) << fvl_dphi[iat][0] 
-//              << std::setw(15) << fvl_dphi[iat][1]
-//              << std::setw(15) << fvl_dphi[iat][2] << std::endl;
-//          }
-
-
     } // end spin
-    // test mohan tmp
-//  test_gamma(this->UHM->LM->DHloc_fixed_x,"this->UHM->LM->DHloc_fixed_x");
-
     if(isstress)
     {
         for(int i=0;i<3;i++)

source/src_lcao/FORCE_gamma_vl_new.cpp

@@ -0,0 +1,28 @@
+#include "FORCE_gamma.h"
+#include "../src_pw/global.h"
+#include "../module_base/timer.h"
+
+void Force_LCAO_gamma::cal_fvl_dphi_new(
+	double*** DM_in,
+	const bool isforce, 
+    const bool isstress,
+    ModuleBase::matrix& fvl_dphi,
+	ModuleBase::matrix& svl_dphi)
+{   
+    ModuleBase::TITLE("Force_LCAO_gamma","cal_fvl_dphi_new");
+    ModuleBase::timer::tick("Force_LCAO_gamma","cal_fvl_dphi_new");
+    int istep = 1;
+    GlobalC::pot.init_pot(istep, GlobalC::pw.strucFac);
+    fvl_dphi.zero_out();
+    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);
+        }
+
+        this->UHM->GG.cal_force_new(DM_in[is], GlobalC::pot.vr_eff1, fvl_dphi);
+    }
+    ModuleBase::timer::tick("Force_LCAO_gamma","cal_fvl_dphi_new");
+}

source/src_lcao/LCAO_nnr.cpp

@@ -86,52 +86,7 @@ void Grid_Technique::cal_nnrg()
 							this->nlocdimg[iat] += nelement; 
 							this->nad[iat]++;
 							++count;
-
-							/*
-							   ofs << std::setw(10) << iat << std::setw(10) << iat2
-							   << std::setw(10) << GlobalC::GridD.getBox(ad).x 
-							   << std::setw(10) << GlobalC::GridD.getBox(ad).y 
-							   << std::setw(10) << GlobalC::GridD.getBox(ad).z 
-							   << std::setw(20) << distance << std::endl;
-							 */
 						}
-						// there is another possibility that i and j are adjacent atoms.
-						// which is that <i|beta> are adjacents while <beta|j> are also
-						// adjacents, these considerations are only considered in k-point
-						// algorithm, 
-						// mohan fix bug 2012-07-03
-						/*
-						else if(distance >= rcut)
-						{
-							for (int ad0 = 0; ad0 < GlobalC::GridD.getAdjacentNum()+1; ++ad0)
-							{
-								const int T0 = GlobalC::GridD.getType(ad0);
-								const int I0 = GlobalC::GridD.getNatom(ad0);
-								const int iat0 = GlobalC::ucell.itia2iat(T0, I0);
-								const int start0 = GlobalC::ucell.itiaiw2iwt(T0, I0, 0);
-
-								tau0 = GlobalC::GridD.getAdjacentTau(ad0);	
-								dtau1 = tau0 - tau1; 
-								dtau2 = tau0 - tau2;
-
-								double distance1 = dtau1.norm() * GlobalC::ucell.lat0;
-								double distance2 = dtau2.norm() * GlobalC::ucell.lat0;
-
-								double rcut1 = GlobalC::ORB.Phi[T1].getRcut() + GlobalC::ucell.infoNL.Beta[T0].get_rcut_max();
-								double rcut2 = GlobalC::ORB.Phi[T2].getRcut() + GlobalC::ucell.infoNL.Beta[T0].get_rcut_max();
-
-								if( distance1 < rcut1 && distance2 < rcut2 )
-								{
-									const int nelement = atom1->nw * atom2->nw;
-									this->nnrg += nelement;
-									this->nlocdimg[iat] += nelement; 
-									this->nad[iat]++;
-									++count;
-									break;	
-								} // dis1, dis2
-							}//ad0
-						}//distance
-						*/
 					}// end iat2
 				}// end ad
 //				GlobalV::ofs_running << " iat=" << iat << " nlocstartg=" << nlocstartg[iat] << " nad=" << nad[iat] << std::endl;
@@ -181,9 +136,6 @@ void Grid_Technique::cal_nnrg()
 	}
 	allocate_find_R2 = true;
 
-//	GlobalV::ofs_running << std::setw(5) << "b1" << std::setw(5) << "b2" << std::setw(5) << "b3"
-//	<< std::setw(8) << "iat" << std::setw(8) << "ad" << std::setw(8) << "iat2"
-//	<< std::setw(8) << "find_R2" << std::setw(8) << "find_R2st" << std::setw(8) << "dis" << std::endl;
 	for (int T1 = 0; T1 < GlobalC::ucell.ntype; T1++)
 	{
 		for (int I1 = 0; I1 < GlobalC::ucell.atoms[T1].na; I1++)
@@ -203,7 +155,6 @@ void Grid_Technique::cal_nnrg()
 				const int I2 = GlobalC::GridD.getNatom(ad);
 				const int iat2 = GlobalC::ucell.itia2iat(T2,I2);
 
-				
 				// if this atom is in this processor.
 				if(this->in_this_processor[iat])
 				{
@@ -216,28 +167,6 @@ void Grid_Technique::cal_nnrg()
 						const int b1 = GlobalC::GridD.getBox(ad).x;
 						const int b2 = GlobalC::GridD.getBox(ad).y;
 						const int b3 = GlobalC::GridD.getBox(ad).z;
-                        
-						// for test
-						/*
-						if( this->cal_RindexAtom(b1, b2, b3, iat2) == 232 )
-						{
-							std::cout << " ====== nnrg =========" << std::endl;
-							std::cout << " index=" << cal_RindexAtom(b1, b2, b3, iat2) << std::endl;
-							std::cout << " iat=" << iat << " iat2=" << iat2 << std::endl;
-							std::cout << " R1 = " << tau1.x << " " << tau1.y << " " << tau1.z << std::endl;
-							std::cout << " R2 = " << GlobalC::GridD.getAdjacentTau(ad).x 
-							<< " " << GlobalC::GridD.getAdjacentTau(ad).y 
-							<< " " << GlobalC::GridD.getAdjacentTau(ad).z << std::endl;
-							std::cout << std::setprecision(25);
-							std::cout << " distance = " << distance << std::endl;
-							std::cout << " box = " << b1 << " " << b2 << " " << b3 << std::endl;
-							std::cout << " rcut = " << rcut << std::endl;
-						}
-						*/
-						
-
-						//std::cout << " iat=" << iat << " find_R2=" << this->cal_RindexAtom(b1, b2, b3, iat2) <<
-						// " b1=" << b1 << " b2=" << b2 << " b3=" << b3 << " iat2=" << iat2 << " distance=" << distance << std::endl;
 
 						// mohan fix bug 2011-06-26, should be '<', not '<='	
 						//			if(distance < rcut)
@@ -258,19 +187,6 @@ void Grid_Technique::cal_nnrg()
 							find_R2[iat][count] = this->cal_RindexAtom(b1, b2, b3, iat2);
 
 
-							/*if(iat==50 && iat2==96)
-							{
-								GlobalV::ofs_running << " ************** iat=" << iat << " count=" << count << " find_R2=" << find_R2[iat][count] << 
-								" b1=" << b1 << " b2=" << b2 << " b3=" << b3 << " iat2=" << iat2 << " distance=" << distance 
-								<< " rcut=" << rcut <<std::endl;
-							}
-							else if(find_R2[iat][count]==10536)
-							{
-								GlobalV::ofs_running << " ************** iat=" << iat << " count=" << count << " find_R2=" << find_R2[iat][count] << 
-								" b1=" << b1 << " b2=" << b2 << " b3=" << b3 << " iat2=" << iat2 << " distance=" << distance 
-								<< " rcut=" << rcut <<std::endl;
-							}*/
-
 							// find_R2st
 							// note: the first must be zero.
 							// find_R2st: start position of each adjacen atom.
@@ -287,21 +203,6 @@ void Grid_Technique::cal_nnrg()
 		}
 	}
 
-	//---------
-	// for test
-	//---------
-	/*
-	GlobalV::ofs_running << " print find_R2 " << std::endl;
-	for(int i=0; i<GlobalC::ucell.nat; i++)
-	{
-		for(int j=0; j<nad[i]; j++)
-		{
-			GlobalV::ofs_running << " i=" << i << " j=" << j << " find_R2=" << find_R2[i][j] << std::endl;
-		}
-	}
-	GlobalV::ofs_running << std::endl;
-	*/
-
 	return;
 }
 
@@ -330,25 +231,10 @@ void Grid_Technique::cal_max_box_index(void)
 		}
 	}
 
-	/*
-	ModuleBase::GlobalFunc::OUT(GlobalV::ofs_running,"maxB1",maxB1);
-	ModuleBase::GlobalFunc::OUT(GlobalV::ofs_running,"maxB2",maxB2);
-	ModuleBase::GlobalFunc::OUT(GlobalV::ofs_running,"maxB3",maxB3);
-	ModuleBase::GlobalFunc::OUT(GlobalV::ofs_running,"minB1",minB1);
-	ModuleBase::GlobalFunc::OUT(GlobalV::ofs_running,"minB2",minB2);
-	ModuleBase::GlobalFunc::OUT(GlobalV::ofs_running,"minB3",minB3);
-	*/
-
 	nB1 = maxB1-minB1+1;
 	nB2 = maxB2-minB2+1;
 	nB3 = maxB3-minB3+1;
 
-	/*
-	ModuleBase::GlobalFunc::OUT(GlobalV::ofs_running,"nB1",nB1);
-	ModuleBase::GlobalFunc::OUT(GlobalV::ofs_running,"nB2",nB2);
-	ModuleBase::GlobalFunc::OUT(GlobalV::ofs_running,"nB3",nB3);
-	*/
-
 	nbox = nB1 * nB2 * nB3;
 
 	//ModuleBase::GlobalFunc::OUT(GlobalV::ofs_running,"nbox",nbox);

source/src_lcao/gint_gamma.h

@@ -35,6 +35,7 @@ class Gint_Gamma : public Grid_Base_Beta
 
 	// (3) calcualte the forces related to grid
 	void cal_force( const double*const vlocal);
+	void cal_force_new(double** DM_in, const double*const vlocal, ModuleBase::matrix& force);
 
 	// (4) calcualte the envelope function
 	void cal_env(const double* wfc, double* rho);
@@ -106,6 +107,7 @@ class Gint_Gamma : public Grid_Base_Beta
 	// for calculatin of < dphi_i | Vlocal | phi_j > for foce calculation
 	// on regular FFT real space grid.
 	void gamma_force(const double*const vlocal) const;
+	void gamma_force_new(const double*const*const DM, const double*const vlocal, ModuleBase::matrix& force);
 
 	void cal_meshball_vlocal(
 		const int na_grid,  						// how many atoms on this (i,j,k) grid
@@ -131,7 +133,22 @@ class Gint_Gamma : public Grid_Base_Beta
 		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]
-	
+
+	void cal_meshball_force(
+		const int grid_index,
+		const int na_grid,  					    // how many atoms on this (i,j,k) grid
+		const int LD_pool,
+		const int*const block_iw,				    // block_iw[na_grid],	index of wave functions for each block
+		const int*const block_size, 			    // block_size[na_grid],	number of columns of a band
+		const int*const block_index,		    	// block_index[na_grid+1], count total number of atomis orbitals
+		const bool*const*const cal_flag,	    	// cal_flag[GlobalC::pw.bxyz][na_grid],	whether the atom-grid distance is larger than cutoff
+		const double*const*const psir_vlbr3,	    // psir_vlbr3[GlobalC::pw.bxyz][LD_pool]
+		const double*const*const dpsir_x,	    // psir_vlbr3[GlobalC::pw.bxyz][LD_pool]
+		const double*const*const dpsir_y,	    // psir_vlbr3[GlobalC::pw.bxyz][LD_pool]
+		const double*const*const dpsir_z,	    // psir_vlbr3[GlobalC::pw.bxyz][LD_pool]
+		const double*const*const DM,
+		ModuleBase::matrix &force);
+
 	// extract the local potentials.
 	// vldr3[GlobalC::pw.bxyz]
     double* get_vldr3(const double* const vlocal, const int ncyz, const int ibx, const int jby, const int kbz) const;