Fix: Support pseudopotential with dense grid (deepmodeling#3678)

* <Fix> Fix the "Error: xhegvd failed!" when running KSDFT with al_lps.cpi * <TEST> Add two integrate tests for LPS 101_PW_lps_pseudopots and 907_OF_LPS * <Fix> Correct the atom label read from BLPS * <Test> Update BLPS part of read_pp_test.cpp
Zhuxuegang2022 · Mar 8, 2024 · 544dc05 · 544dc05
1 parent 9c01ed4
commit 544dc05
Show file tree

Hide file tree

Showing 17 changed files with 10,190 additions and 5 deletions.
diff --git a/source/module_cell/pseudo.cpp b/source/module_cell/pseudo.cpp
@@ -100,7 +100,7 @@ void pseudo::set_pseudo_h(const Pseudopot_upf &upf)
 
 	// mohan update 2021-02-22
 	//  max number of points in the atomic radial mesh
-	int ndmx = 2000; 
+	int ndmx = 200000; 
 	if (this->mesh > ndmx)
 	{
 		std::cout << "\n set_pseudo_h, too many grid points,";

diff --git a/source/module_cell/read_pp_blps.cpp b/source/module_cell/read_pp_blps.cpp
@@ -1,4 +1,6 @@
 #include "read_pp.h"
+#include "module_base/atom_in.h"
+#include "module_base/element_name.h"
 
 int Pseudopot_upf::read_pseudo_blps(std::ifstream &ifs)
 {
@@ -37,6 +39,16 @@ int Pseudopot_upf::read_pseudo_blps(std::ifstream &ifs)
     this->zp = static_cast<int>(zion);
     ifs.ignore(300, '\n');
 
+    atom_in ai;
+    for (auto each_type:  ModuleBase::element_name)
+    {
+        if (zatom == ai.atom_Z[each_type])
+        {
+            this->psd = each_type;
+            break;
+        }
+    }
+
     int pspcod, pspxc, lloc, r2well;
     ifs >> pspcod >> pspxc >> this->lmax >> lloc >> this->mesh >> r2well;
 

diff --git a/source/module_cell/test/read_pp_test.cpp b/source/module_cell/test/read_pp_test.cpp
@@ -494,7 +494,7 @@ TEST_F(ReadPPTest, BLPS)
 	EXPECT_FALSE(upf->tvanp);
 	EXPECT_FALSE(upf->has_so);
 	EXPECT_EQ(upf->nbeta,0);
-	EXPECT_EQ(upf->psd,"Silicon");
+	EXPECT_EQ(upf->psd,"Si");
 	EXPECT_EQ(upf->zp,4);
 	EXPECT_EQ(upf->lmax,0);
 	EXPECT_EQ(upf->mesh,1601);

diff --git a/source/module_elecstate/module_charge/charge.cpp b/source/module_elecstate/module_charge/charge.cpp
@@ -360,9 +360,15 @@ void Charge::atomic_rho(const int spin_number_need,
                                     rhoatm[ir] = atom->ncpp.rho_at[ir] / ModuleBase::FOUR_PI / r2;
                                 }
                                 rhoatm[0]
-                                    = pow((rhoatm[2] / rhoatm[1]), 1. / (atom->ncpp.r[2] - atom->ncpp.r[1])); // zws add
-                                rhoatm[0] = pow(rhoatm[0], atom->ncpp.r[1]);
-                                rhoatm[0] = rhoatm[1] / rhoatm[0];
+                                    = pow((rhoatm[2] / rhoatm[1]), atom->ncpp.r[1] / (atom->ncpp.r[2] - atom->ncpp.r[1])); // zws add, sunliang updated 2024-03-04
+                                if (rhoatm[0] < 1e-12)
+                                {
+                                    rhoatm[0] = rhoatm[1];
+                                }
+                                else
+                                {
+                                    rhoatm[0] = rhoatm[1] / rhoatm[0];
+                                }
 
                                 double charge = 0.0;
                                 ModuleBase::Integral::Simpson_Integral(atom->ncpp.msh,