Skip to content

Commit 42ca798

Browse files
committed
1. PCA_thr element
1 parent bba6408 commit 42ca798

7 files changed

Lines changed: 30 additions & 15 deletions

File tree

source/source_hamilt/module_xc/exx_info.h

Lines changed: 2 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -56,7 +56,7 @@ struct Exx_Info
5656
bool coul_moment = false;
5757
bool rotate_abfs = false;
5858

59-
double pca_threshold = 0;
59+
std::vector<double> pca_threshold = {0};
6060
std::vector<std::string> files_abfs;
6161
std::vector<std::string> files_shrink_abfs;
6262
double C_threshold = 0;
@@ -92,7 +92,7 @@ struct Exx_Info
9292
double tolerence = 1E-12;
9393
std::vector<std::string> files_jles;
9494

95-
double pca_threshold = 0;
95+
std::vector<double> pca_threshold = {0};
9696
std::vector<std::string> files_abfs;
9797

9898
double kmesh_times = 4;

source/source_io/module_parameter/input_parameter.h

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -545,7 +545,7 @@ struct Input_para
545545
///< the evaluation of Fock exchange
546546
double exx_mixing_beta = 1.0; ///< mixing_beta for outer-loop when exx_separate_loop=1
547547
std::string exx_real_number = "default"; ///< exx calculated in real or complex
548-
double exx_pca_threshold = 0.0001; ///< threshold to screen on-site ABFs in exx
548+
std::vector<double> exx_pca_threshold = {0.0001}; ///< threshold to screen on-site ABFs in exx
549549
double exx_c_threshold = 0.0001; ///< threshold to screen C matrix in exx
550550
double exx_v_threshold = 0.1; ///< threshold to screen C matrix in exx
551551
double exx_dm_threshold = 0.0001; ///< threshold to screen density matrix in exx

source/source_io/module_parameter/read_input_item_exx_dftu.cpp

Lines changed: 10 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -245,12 +245,20 @@ void ReadInput::item_exx()
245245
Input_Item item("exx_pca_threshold");
246246
item.annotation = "threshold to screen on-site ABFs in exx";
247247
item.category = "Exact Exchange (LCAO)";
248-
item.type = "Real";
248+
item.type = "Vector of Real (n values where n is the number of atomic types)";
249249
item.description = "To accelerate the evaluation of four-center integrals (), the product of atomic orbitals are expanded in the basis of auxiliary basis functions (ABF): . The size of the ABF (i.e. number of ) is reduced using principal component analysis. When a large PCA threshold is used, the number of ABF will be reduced, hence the calculation becomes faster. However, this comes at the cost of computational accuracy. A relatively safe choice of the value is 1e-4.";
250250
item.default_value = "1E-4";
251251
item.unit = "";
252252
item.availability = "";
253-
read_sync_double(input.exx_pca_threshold);
253+
item.read_value = [](const Input_Item& item, Parameter& para) {
254+
para.input.exx_pca_threshold.clear();
255+
size_t count = item.get_size();
256+
for (int i = 0; i < count; i++)
257+
{
258+
para.input.exx_pca_threshold.push_back(std::stod(item.str_values[i]));
259+
}
260+
};
261+
sync_doublevec(input.exx_pca_threshold, para.input.ntype, 0.0001);
254262
this->add_item(item);
255263
}
256264
{

source/source_io/test/read_input_ptest.cpp

Lines changed: 2 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -288,7 +288,8 @@ TEST_F(InputParaTest, ParaRead)
288288
EXPECT_EQ(param.inp.exx_hybrid_step, 100);
289289
EXPECT_DOUBLE_EQ(std::stod(param.inp.exx_fock_lambda[0]), 0.3);
290290
EXPECT_DOUBLE_EQ(param.inp.exx_mixing_beta, 1.0);
291-
EXPECT_DOUBLE_EQ(param.inp.exx_pca_threshold, 0);
291+
ASSERT_EQ(param.inp.exx_pca_threshold.size(), 1);
292+
EXPECT_DOUBLE_EQ(param.inp.exx_pca_threshold[0], 0);
292293
EXPECT_DOUBLE_EQ(param.inp.exx_c_threshold, 0);
293294
EXPECT_DOUBLE_EQ(param.inp.exx_v_threshold, 0);
294295
EXPECT_DOUBLE_EQ(param.inp.exx_dm_threshold, 0);

source/source_lcao/module_ri/RPA_LRI.hpp

Lines changed: 7 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -52,7 +52,12 @@ void RPA_LRI<T, Tdata>::postSCF(const UnitCell& ucell,
5252
this->out_eigen_vector(parav, psi);
5353
this->out_struc(ucell);
5454

55-
std::cout << "rpa_pca_threshold: " << this->info.pca_threshold << std::endl;
55+
std::cout << "rpa_pca_threshold:";
56+
for (const double value: this->info.pca_threshold)
57+
{
58+
std::cout << " " << value;
59+
}
60+
std::cout << std::endl;
5661
std::cout << "rpa_ccp_rmesh_times: " << this->info.ccp_rmesh_times << std::endl;
5762
std::cout << "rpa_lcao_exx(Ha): " << std::fixed << std::setprecision(15) << exx_cut_coulomb->Eexx / 2.0 << std::endl;
5863

@@ -162,7 +167,7 @@ void RPA_LRI<T, Tdata>::cal_postSCF_exx(const elecstate::DensityMatrix<T, Tdata>
162167
orb,
163168
this->lcaos,
164169
this->info.kmesh_times,
165-
this->info.shrink_abfs_pca_thr);
170+
std::vector<double>(ucell.ntype, this->info.shrink_abfs_pca_thr));
166171
if (this->info.files_shrink_abfs.empty())
167172
{
168173
this->abfs_shrink = abfs_same_atom;

source/source_lcao/module_ri/exx_abfs-construct_orbs.cpp

Lines changed: 6 additions & 5 deletions
Original file line numberDiff line numberDiff line change
@@ -82,10 +82,10 @@ std::vector<std::vector<std::vector<Numerical_Orbital_Lm>>> Exx_Abfs::Construct_
8282
const LCAO_Orbitals& orb,
8383
const std::vector<std::vector<std::vector<Numerical_Orbital_Lm>>> &orbs,
8484
const double kmesh_times_mot,
85-
const double times_threshold )
85+
const std::vector<double>& times_threshold )
8686
{
8787
ModuleBase::TITLE("Exx_Abfs::Construct_Orbs::abfs_same_atom");
88-
if(times_threshold>1)
88+
if(std::all_of(times_threshold.begin(), times_threshold.end(), [](const double value) { return value > 1; }))
8989
{ return std::vector<std::vector<std::vector<Numerical_Orbital_Lm>>>(orb.get_ntype()); }
9090

9191
const std::vector<std::vector<std::vector<std::vector<double>>>>
@@ -265,9 +265,9 @@ std::vector<std::vector<std::vector<std::vector<double>>>> Exx_Abfs::Construct_O
265265
const std::vector<std::vector<std::vector<Numerical_Orbital_Lm>>> &abfs,
266266
const std::vector<std::vector<std::vector<Numerical_Orbital_Lm>>> &orbs,
267267
const double kmesh_times_mot,
268-
const double times_threshold )
268+
const std::vector<double>& times_threshold )
269269
{
270-
if(times_threshold>1)
270+
if(std::all_of(times_threshold.begin(), times_threshold.end(), [](const double value) { return value > 1; }))
271271
return std::vector<std::vector<std::vector<std::vector<double>>>>(abfs.size());
272272

273273
const std::vector<std::vector<std::pair<std::vector<double>,RI::Tensor<double>>>>
@@ -276,6 +276,7 @@ std::vector<std::vector<std::vector<std::vector<double>>>> Exx_Abfs::Construct_O
276276
const std::vector<std::vector<std::vector<std::vector<double>>>> psis = get_psi( abfs );
277277
std::vector<std::vector<std::vector<std::vector<double>>>> psis_new( psis.size() );
278278

279+
assert(times_threshold.size() == eig.size());
279280
for( size_t T=0; T!=eig.size(); ++T )
280281
{
281282
double eig_value_max = 0;
@@ -285,7 +286,7 @@ std::vector<std::vector<std::vector<std::vector<double>>>> Exx_Abfs::Construct_O
285286
//ofs<<T<<"\t"<<L<<"\t"<<M<<"\t"<<eig[T][L].first[M]<<std::endl;
286287
eig_value_max = std::max( eig_value_max, eig[T][L].first[M] );
287288
}
288-
const double eig_value_threshold = eig_value_max * times_threshold;
289+
const double eig_value_threshold = eig_value_max * times_threshold[T];
289290

290291
//ofs<<"eig_value_max:\t"<<eig_value_max<<std::endl;
291292
//ofs<<"eig_value_threshold:\t"<<eig_value_threshold<<std::endl;

source/source_lcao/module_ri/exx_abfs-construct_orbs.h

Lines changed: 2 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -25,7 +25,7 @@ class Exx_Abfs::Construct_Orbs
2525
const LCAO_Orbitals& orb,
2626
const std::vector<std::vector<std::vector<Numerical_Orbital_Lm>>> &lcaos,
2727
const double kmesh_times_mot,
28-
const double times_threshold=0);
28+
const std::vector<double>& times_threshold);
2929

3030
static void print_orbs_size(
3131
const UnitCell& ucell,
@@ -84,7 +84,7 @@ class Exx_Abfs::Construct_Orbs
8484
const std::vector<std::vector<std::vector<Numerical_Orbital_Lm>>> &abfs,
8585
const std::vector<std::vector<std::vector<Numerical_Orbital_Lm>>> &orbs,
8686
const double kmesh_times_mot,
87-
const double times_threshold );
87+
const std::vector<double>& times_threshold );
8888

8989
static std::vector<std::vector<std::vector<std::vector<double>>>> div_r(
9090
const std::vector<std::vector<std::vector<std::vector<double>>>> &psirs,

0 commit comments

Comments
 (0)