Skip to content

Commit cbaccdf

Browse files
kirk0830WHUweiqingzhou
kirk0830
and
WHUweiqingzhou
authored
Feature: Quasiatomic Orbital (QO) implement - immediate (deepmodeling#3236)
* building two center overlap * 20231120 * 20231121 * 20231124 * serial implementation done * because rcut_max is not set in previous, after addition, the corresponding unit test is forgotten to update. * full implementation of QO analysis, next is to improve accuracy of associated Laguerre Polynomial * Accurancy improved. C++-side in principle no problem for use anymore * provide python-end qo analysis into tools folder. In readme.md the usage is introduced. * recover changed INPUT files * support strategy full, use filter in postprocessing to select basis to reconstruct hamiltonian * - add new derived class pswfc_radials - extend strategy of derived class hydrogen_radials (but tested to perform ill) * implement qo_basis = pswfc * update python-end * Oops! update Makefile.Objects * add ceil() when computing supercell * update documentation * change kvec_c to kvec_d when folding * add integrated test and support output of HqoR matrix * correct CMakeLists.txt * correct CMakeLists.txt * correct CMakeLists.txt * reduce number of hydrogen-like orbitals generated with "energy" strategy * update reference value of hydrogen_radials unittest * update example and debugging the unfolding_Hk function, not completed yet * correct the method to find supercells * recover example/scf/lcao_Cu * recover examples/scf/lcao_Cu * complete necessary annotations * support the use of qo_thr for qo_basis = pswfc, flexibly controls the rcut of pswfc * add more unittest to test the symmetrical-kpoints cancellation on imaginary parts * update unittest on toqo class * complement zero_out before every so-called folding_Hk * update debugged python-end analysis package and correct CMakeLists.txt * update README.md of QO in tools/qo * Update INPUT in tools/qo/examples/abacus_input * delete redundant example files * update docs according to the newest implementation * remove commented out codes accordint to code review * provide example and update more detailed README.md * change default strategy to `minimal` for memory save concern --------- Co-authored-by: wqzhou <33364058+WHUweiqingzhou@users.noreply.github.com>
1 parent 8cf8f4b commit cbaccdf

68 files changed

Lines changed: 5478 additions & 30 deletions

Some content is hidden

Large Commits have some content hidden by default. Use the searchbox below for content that may be hidden.

‎.gitignore‎

Lines changed: 3 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -18,4 +18,6 @@ build
1818
dist
1919
.idea
2020
toolchain.tar.gz
21-
time.json
21+
time.json
22+
*.pyc
23+
__pycache__

‎docs/advanced/input_files/input-main.md‎

Lines changed: 57 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -370,6 +370,12 @@
370370
- [alpha\_trial](#alpha_trial)
371371
- [sccut](#sccut)
372372
- [sc\_file](#sc_file)
373+
- [Quasiatomic Orbital (QO) analysis](#quasiatomic-orbital-qo-analysis)
374+
- [qo\_switch](#qo_switch)
375+
- [qo\_basis](#qo_basis)
376+
- [qo\_strategy](#qo_strategy)
377+
- [qo\_screening\_coeff](#qo_screening_coeff)
378+
- [qo\_thr](#qo_thr)
373379

374380
[back to top](#full-list-of-input-keywords)
375381

@@ -3315,6 +3321,8 @@ These variables are used to control the usage of implicit solvation model. This
33153321
- **Default**: 0.00037
33163322
- **Unit**: $Bohr^{-3}$
33173323

3324+
[back to top](#full-list-of-input-keywords)
3325+
33183326
## Deltaspin
33193327

33203328
These variables are used to control the usage of deltaspin functionality.
@@ -3431,4 +3439,53 @@ for `nspin 2` case. The difference is that `lambda`, `target_mag`, and `constrai
34313439

34323440
- **Default**: none
34333441

3442+
[back to top](#full-list-of-input-keywords)
3443+
3444+
## Quasiatomic Orbital (QO) analysis
3445+
3446+
These variables are used to control the usage of QO analysis.
3447+
3448+
### qo_switch
3449+
3450+
- **Type**: Boolean
3451+
- **Description**: whether to let ABACUS output QO analysis required files
3452+
- **Default**: 0
3453+
3454+
### qo_basis
3455+
3456+
- **Type**: String
3457+
- **Description**: specify the type of atomic basis
3458+
- `pswfc`: use the pseudowavefunction in pseudopotential files as atomic basis. To use this option, please make sure in pseudopotential file there is pswfc in it.
3459+
- `hydrogen`: generate hydrogen-like atomic basis, whose charge is read from pseudopotential files presently.
3460+
3461+
*warning: to use* `pswfc` *, please use norm-conserving pseudopotentials with pseudowavefunctions, SG15 pseudopotentials cannot support this option.*
3462+
- **Default**: `hydrogen`
3463+
3464+
### qo_strategy
3465+
3466+
- **Type**: String
3467+
- **Availability**: for `qo_basis hydrogen` only.
3468+
- **Description**: specify the strategy to generate hydrogen-like orbitals
3469+
- `minimal`: according to principle quantum number of the highest occupied state, generate only nodeless orbitals, for example Cu, only generate 1s, 2p, 3d and 4f orbitals (for Cu, 4s is occupied, thus $n_{max} = 4$)
3470+
- `full`: similarly according to the maximal principle quantum number, generate all possible orbitals, therefore for Cu, for example, will generate 1s, 2s, 2p, 3s, 3p, 3d, 4s, 4p, 4d, 4f.
3471+
- `energy`: will generate hydrogen-like orbitals according to Aufbau principle. For example the Cu (1s2 2s2 2p6 3s2 3p6 3d10 4s1), will generate these orbitals.
3472+
3473+
*warning: to use* `full`, *generation strategy may cause the space spanned larger than the one spanned by numerical atomic orbitals, in this case, must filter out orbitals in some way*
3474+
- **Default**: `minimal`
3475+
3476+
### qo_screening_coeff
3477+
3478+
- **Type**: Real
3479+
- **Availability**: for `qo_basis pswfc` only.
3480+
- **Description**: a screening factor $e^{-\eta|\mathbf{r}|}$ is multiplied to the pswfc to mimic the behavior of some kind of electron. $\eta$ is the screening coefficient. Presently one scalar value can be passed to ABACUS, therefore all atom types use the same value.
3481+
- **Default**: 0.1
3482+
- **Unit**: Bohr^-1
3483+
3484+
### qo_thr
3485+
3486+
- **Type**: Real
3487+
- **Description**: the convergence threshold determining the cutoff of generated orbital. Lower threshold will yield orbital with larger cutoff radius.
3488+
- **Default**: 1.0e-6
3489+
3490+
34343491
[back to top](#full-list-of-input-keywords)

‎source/Makefile.Objects‎

Lines changed: 6 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -106,6 +106,7 @@ OBJS_MAIN=main.o\
106106
driver_run.o\
107107

108108
OBJS_BASE=abfs-vector3_order.o\
109+
assoc_laguerre.o\
109110
complexarray.o\
110111
complexmatrix.o\
111112
clebsch_gordan_coeff.o\
@@ -326,6 +327,8 @@ OBJS_ORBITAL=ORB_atomic.o\
326327
parallel_2d.o\
327328
parallel_orbitals.o\
328329
atomic_radials.o\
330+
hydrogen_radials.o\
331+
pswfc_radials.o\
329332
beta_radials.o\
330333
numerical_radial.o\
331334
radial_collection.o\
@@ -419,6 +422,8 @@ OBJS_IO=input.o\
419422
restart.o\
420423
binstream.o\
421424
to_wannier90.o\
425+
to_qo.o\
426+
to_qo_tools.o\
422427
to_wannier90_pw.o\
423428
to_wannier90_lcao_in_pw.o\
424429
to_wannier90_lcao.o\
@@ -436,6 +441,7 @@ OBJS_IO=input.o\
436441
output_rho.o\
437442
output_potential.o\
438443
output_mat_sparse.o\
444+
output_radial.o\
439445

440446
OBJS_IO_LCAO=cal_r_overlap_R.o\
441447
write_orb_info.o\

‎source/module_base/CMakeLists.txt‎

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -10,6 +10,7 @@ endif()
1010
add_library(
1111
base
1212
OBJECT
13+
assoc_laguerre.cpp
1314
clebsch_gordan_coeff.cpp
1415
complexarray.cpp
1516
complexmatrix.cpp

‎source/module_base/assoc_laguerre.cpp‎

Lines changed: 137 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,137 @@
1+
#include "module_base/assoc_laguerre.h"
2+
#include "module_base/global_function.h"
3+
//#include <tr1/cmath> // use cmath the factorial function
4+
#include <cmath>
5+
Assoc_Laguerre::Assoc_Laguerre()
6+
{
7+
}
8+
9+
Assoc_Laguerre::~Assoc_Laguerre()
10+
{
11+
}
12+
13+
void Assoc_Laguerre::generate(const int &n, const int &l, const double ns, double* const &s, double* L)
14+
{
15+
for(int i = 0; i < ns; i++)
16+
{
17+
L[i] = this->value(n, l, s[i]);
18+
}
19+
}
20+
21+
void Assoc_Laguerre::generate(const int &n, const int &l, std::vector<double> &x, std::vector<double> &y)
22+
{
23+
for(int i = 0; i < x.size(); i++)
24+
{
25+
y[i] = this->value(n, l, x[i]);
26+
}
27+
}
28+
29+
double Assoc_Laguerre::laguerre(const int &n, const double x)
30+
{
31+
if(n == 0)
32+
{
33+
return 1;
34+
}
35+
else if(n == 1)
36+
{
37+
return -x + 1;
38+
}
39+
else if(n == 2)
40+
{
41+
return 0.5 * x * x - 2 * x + 1;
42+
}
43+
else if(n == 3)
44+
{
45+
return -x * x * x / 6.0 + 3.0 * x * x / 2.0 - 3.0 * x + 1;
46+
}
47+
else if(n >= 4)
48+
{
49+
double n_ = static_cast<double>(n);
50+
double first = (2*n_ - 1 - x)/n_ * Assoc_Laguerre::laguerre(n-1, x);
51+
double second = (n_ - 1)/n_ * Assoc_Laguerre::laguerre(n-2, x);
52+
return first - second;
53+
}
54+
else
55+
{
56+
ModuleBase::WARNING_QUIT("Assoc_Laguerre::laguerre", "n is out of range");
57+
return 0;
58+
}
59+
}
60+
61+
double Assoc_Laguerre::associate_laguerre(const int &n, const double x, const int &a)
62+
{
63+
// formula from https://en.wikipedia.org/wiki/Laguerre_polynomials
64+
double n_ = static_cast<double>(n);
65+
double a_ = static_cast<double>(a);
66+
if(n == 0)
67+
{
68+
return 1;
69+
}
70+
else if(n == 1)
71+
{
72+
return -x + 1 + a_;
73+
}
74+
else if(n == 2)
75+
{
76+
return 0.5 * (x*x - 2*(a_+2)*x + (a_+1)*(a_+2));
77+
}
78+
else if(n == 3)
79+
{
80+
return -x*x*x/6.0 + (a_+3)*x*x/2.0 - (a_+2)*(a_+3)*x/2.0 + (a_+1)*(a_+2)*(a_+3)/6.0;
81+
}
82+
else if(n >= 4)
83+
{
84+
double first = (2*n_ - 1 + a_ - x)/n_ * this->associate_laguerre(n-1, x, a);
85+
double second = (n_ + a_ - 1)/n_ * this->associate_laguerre(n-2, x, a);
86+
return first - second;
87+
}
88+
else
89+
{
90+
ModuleBase::WARNING_QUIT("Assoc_Laguerre::associate_laguerre", "n is out of range");
91+
return 0;
92+
}
93+
}
94+
95+
int Assoc_Laguerre::factorial(const int &n)
96+
{
97+
if(n == 0)
98+
{
99+
return 1;
100+
}
101+
else if(n > 0)
102+
{
103+
return n * this->factorial(n-1);
104+
}
105+
else
106+
{
107+
ModuleBase::WARNING_QUIT("Assoc_Laguerre::factorial", "n is out of range");
108+
return 0;
109+
}
110+
}
111+
112+
double Assoc_Laguerre::value(const int &n, const int &l, const double &s)
113+
{
114+
int k_ = 2*l + 1;
115+
int n_ = n - l - 1;
116+
if(k_ < 0)
117+
{
118+
ModuleBase::WARNING_QUIT("Assoc_Laguerre::value", "k is out of range");
119+
return 0;
120+
}
121+
if(n_ < 0)
122+
{
123+
ModuleBase::WARNING_QUIT("Assoc_Laguerre::value", "n is out of range");
124+
return 0;
125+
}
126+
double L = 0;
127+
for(int iq = 0; iq <= n_; iq++)
128+
{
129+
L += std::pow(-s, iq) *
130+
static_cast<double>(this->factorial(n_ + k_)) /
131+
static_cast<double>(this->factorial(n_ - iq)) /
132+
static_cast<double>(this->factorial(k_ + iq)) /
133+
static_cast<double>(this->factorial(iq));
134+
}
135+
//L = std::tr1::assoc_laguerre(n_, k_, s); // use standard library
136+
return L;
137+
}

‎source/module_base/assoc_laguerre.h‎

Lines changed: 47 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,47 @@
1+
#ifndef ASSOC_LAGUEERRE_H
2+
#define ASSOC_LAGUEERRE_H
3+
4+
#include <vector>
5+
#include <string>
6+
7+
class Assoc_Laguerre
8+
{
9+
public:
10+
Assoc_Laguerre();
11+
~Assoc_Laguerre();
12+
/// @brief generate the associated Laguerre polynomial (overloaded for double*)
13+
/// @param n principal quantum number
14+
/// @param l orbital quantum number
15+
/// @param ns number of x-coordinates
16+
/// @param s x-coordinates
17+
/// @param L y-coordinates
18+
void generate(const int &n, const int &l, const double ns, double* const &s, double* L);
19+
/// @brief generate the associated Laguerre polynomial (overloaded for std::vector)
20+
/// @param n principal quantum number
21+
/// @param l orbital quantum number
22+
/// @param x x-coordinates in std::vector
23+
/// @param y y-coordinates in std::vector
24+
void generate(const int &n, const int &l, std::vector<double> &x, std::vector<double> &y);
25+
/// @brief Laguerre polynomial
26+
/// @param n degree of the polynomial
27+
/// @param x radial coordinate
28+
/// @return L_n(x)
29+
double laguerre(const int &n, const double x);
30+
/// @brief recursive relationship to find the associated Laguerre polynomial
31+
/// @param n degree of the polynomial
32+
/// @param x radial coordinate
33+
/// @param a order of the polynomial
34+
/// @return L^(a)_n(x)
35+
double associate_laguerre(const int &n, const double x, const int &a);
36+
/// @brief wrapper for associate_laguerre
37+
/// @param n principal quantum number
38+
/// @param l orbital quantum number
39+
/// @param s radial coordinate
40+
/// @return L^(2l+1)_(n-l-1)(s)
41+
double value(const int &n, const int &l, const double &s);
42+
/// @brief factorial function
43+
/// @param n
44+
/// @return n!
45+
int factorial(const int &n);
46+
};
47+
#endif // ASSOC_LAGUEERRE_H

‎source/module_base/atom_in.h‎

Lines changed: 17 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -91,7 +91,23 @@ class atom_in
9191
{"Meitnerium", 109}, {"Darmstadtium", 110}, {"Roentgenium", 111}, {"Copernicium", 112}, {"Nihonium", 113}, {"Flerovium", 114},
9292
{"Moscovium", 115}, {"Livermorium", 116}, {"Tennessine", 117}, {"Oganesson", 118}
9393
};
94-
94+
std::map<std::string, int> principle_quantum_number
95+
= {
96+
{"H", 1}, {"He", 1}, {"Li", 2}, {"Be", 2}, {"B", 2}, {"C", 2}, {"N", 2}, {"O", 2}, {"F", 2},
97+
{"Ne", 1}, {"Na", 2}, {"Mg", 3}, {"Al", 3}, {"Si", 3}, {"P", 3}, {"S", 3}, {"Cl", 3}, {"Ar", 2},
98+
{"K", 3}, {"Ca", 4}, {"Sc", 4}, {"Ti", 4}, {"V", 4}, {"Cr", 4}, {"Mn", 4}, {"Fe", 4}, {"Co", 4},
99+
{"Ni", 4}, {"Cu", 4}, {"Zn", 4}, {"Ga", 4}, {"Ge", 4}, {"As", 4}, {"Se", 4}, {"Br", 4}, {"Kr", 3},
100+
{"Rb", 4}, {"Sr", 5}, {"Y", 5}, {"Zr", 5}, {"Nb", 5}, {"Mo", 5}, {"Tc", 5}, {"Ru", 5}, {"Rh", 5},
101+
{"Pd", 5}, {"Ag", 5}, {"Cd", 5}, {"In", 5}, {"Sn", 5}, {"Sb", 5}, {"Te", 5}, {"I", 5}, {"Xe", 4},
102+
{"Cs", 5}, {"Ba", 6}, {"La", 6}, {"Ce", 6}, {"Pr", 6}, {"Nd", 6}, {"Pm", 6}, {"Sm", 6}, {"Eu", 6},
103+
{"Gd", 6}, {"Tb", 6}, {"Dy", 6}, {"Ho", 6}, {"Er", 6}, {"Tm", 6}, {"Yb", 6}, {"Lu", 6}, {"Hf", 6},
104+
{"Ta", 6}, {"W", 6}, {"Re", 6}, {"Os", 6}, {"Ir", 6}, {"Pt", 6}, {"Au", 6}, {"Hg", 6}, {"Tl", 6},
105+
{"Pb", 6}, {"Bi", 6}, {"Po", 6}, {"At", 6}, {"Rn", 6}, {"Fr", 7}, {"Ra", 7}, {"Ac", 7}, {"Th", 7},
106+
{"Pa", 7}, {"U", 7}, {"Np", 7}, {"Pu", 7}, {"Am", 7}, {"Cm", 7}, {"Bk", 7}, {"Cf", 7}, {"Es", 7},
107+
{"Fm", 7}, {"Md", 7}, {"No", 7}, {"Lr", 7}, {"Rf", 7}, {"Db", 7}, {"Sg", 7}, {"Bh", 7}, {"Hs", 7},
108+
{"Mt", 7}, {"Ds", 7}, {"Rg", 7}, {"Cn", 7}, {"Nh", 7}, {"Fl", 7}, {"Mc", 7}, {"Lv", 7}, {"Ts", 7},
109+
{"Og", 7}
110+
};
95111
};
96112

97113
#endif

‎source/module_base/formatter_contextfmt.cpp‎

Lines changed: 10 additions & 10 deletions
Original file line numberDiff line numberDiff line change
@@ -82,9 +82,9 @@ void formatter::ContextFmt::set_context(std::vector<std::string> phys_fmt) {
8282
}
8383
}
8484

85-
template<> formatter::ContextFmt& formatter::ContextFmt::operator<< <double>(double const&);
86-
template<> formatter::ContextFmt& formatter::ContextFmt::operator<< <int>(int const&);
87-
template<> formatter::ContextFmt& formatter::ContextFmt::operator<< <char>(char const&);
85+
template formatter::ContextFmt& formatter::ContextFmt::operator<< <double>(double const&);
86+
template formatter::ContextFmt& formatter::ContextFmt::operator<< <int>(int const&);
87+
template formatter::ContextFmt& formatter::ContextFmt::operator<< <char>(char const&);
8888

8989
template <>
9090
formatter::ContextFmt& formatter::ContextFmt::operator<<(const std::string& value) {
@@ -122,14 +122,14 @@ formatter::ContextFmt& formatter::ContextFmt::operator<< (char const* value)
122122
return *this << value_;
123123
}
124124

125-
template<> formatter::ContextFmt& formatter::ContextFmt::operator<< <double>(std::vector<double> const& value);
126-
template<> formatter::ContextFmt& formatter::ContextFmt::operator<< <std::string>(std::vector<std::string> const& value);
127-
template<> formatter::ContextFmt& formatter::ContextFmt::operator<< <int>(std::vector<int> const& value);
125+
template formatter::ContextFmt& formatter::ContextFmt::operator<< <double>(std::vector<double> const& value);
126+
template formatter::ContextFmt& formatter::ContextFmt::operator<< <std::string>(std::vector<std::string> const& value);
127+
template formatter::ContextFmt& formatter::ContextFmt::operator<< <int>(std::vector<int> const& value);
128128

129-
template<> formatter::ContextFmt& formatter::ContextFmt::operator<< <double>(double*& value);
130-
template<> formatter::ContextFmt& formatter::ContextFmt::operator<< <int>(int*& value);
131-
template<> formatter::ContextFmt& formatter::ContextFmt::operator<< <std::string>(std::string*& value);
132-
template<> formatter::ContextFmt& formatter::ContextFmt::operator<< <char>(char*& value);
129+
template formatter::ContextFmt& formatter::ContextFmt::operator<< <double>(double*& value);
130+
template formatter::ContextFmt& formatter::ContextFmt::operator<< <int>(int*& value);
131+
template formatter::ContextFmt& formatter::ContextFmt::operator<< <std::string>(std::string*& value);
132+
template formatter::ContextFmt& formatter::ContextFmt::operator<< <char>(char*& value);
133133

134134
void formatter::ContextFmt::reset() {
135135
Table::reset();

‎source/module_base/global_variable.cpp‎

Lines changed: 9 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -289,4 +289,13 @@ int sc_scf_nmin = 2;
289289
double alpha_trial = 0.01; // eV/uB^2
290290
double sccut = 3; // eV/uB
291291
std::string sc_file = "none";
292+
293+
//==========================================================
294+
// Quasiatomic orbital related
295+
//==========================================================
296+
bool qo_switch = false;
297+
std::string qo_basis = "hydrogen";
298+
std::string qo_strategy = "minimal";
299+
double qo_thr = 1.0e-6;
300+
std::vector<double> qo_screening_coeff = {};
292301
} // namespace GlobalV

0 commit comments

Comments
 (0)