Skip to content

Commit

Permalink
Feature: Quasiatomic Orbital (QO) implement - immediate (deepmodeling…
Browse files Browse the repository at this point in the history
…#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 <[email protected]>
  • Loading branch information
kirk0830 and WHUweiqingzhou authored Jan 2, 2024
1 parent 8cf8f4b commit cbaccdf
Show file tree
Hide file tree
Showing 68 changed files with 5,478 additions and 30 deletions.
4 changes: 3 additions & 1 deletion .gitignore
Original file line number Diff line number Diff line change
Expand Up @@ -18,4 +18,6 @@ build
dist
.idea
toolchain.tar.gz
time.json
time.json
*.pyc
__pycache__
57 changes: 57 additions & 0 deletions docs/advanced/input_files/input-main.md
Original file line number Diff line number Diff line change
Expand Up @@ -370,6 +370,12 @@
- [alpha\_trial](#alpha_trial)
- [sccut](#sccut)
- [sc\_file](#sc_file)
- [Quasiatomic Orbital (QO) analysis](#quasiatomic-orbital-qo-analysis)
- [qo\_switch](#qo_switch)
- [qo\_basis](#qo_basis)
- [qo\_strategy](#qo_strategy)
- [qo\_screening\_coeff](#qo_screening_coeff)
- [qo\_thr](#qo_thr)

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

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

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

## Deltaspin

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

- **Default**: none

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

## Quasiatomic Orbital (QO) analysis

These variables are used to control the usage of QO analysis.

### qo_switch

- **Type**: Boolean
- **Description**: whether to let ABACUS output QO analysis required files
- **Default**: 0

### qo_basis

- **Type**: String
- **Description**: specify the type of atomic basis
- `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.
- `hydrogen`: generate hydrogen-like atomic basis, whose charge is read from pseudopotential files presently.

*warning: to use* `pswfc` *, please use norm-conserving pseudopotentials with pseudowavefunctions, SG15 pseudopotentials cannot support this option.*
- **Default**: `hydrogen`

### qo_strategy

- **Type**: String
- **Availability**: for `qo_basis hydrogen` only.
- **Description**: specify the strategy to generate hydrogen-like orbitals
- `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$)
- `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.
- `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.

*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*
- **Default**: `minimal`

### qo_screening_coeff

- **Type**: Real
- **Availability**: for `qo_basis pswfc` only.
- **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.
- **Default**: 0.1
- **Unit**: Bohr^-1

### qo_thr

- **Type**: Real
- **Description**: the convergence threshold determining the cutoff of generated orbital. Lower threshold will yield orbital with larger cutoff radius.
- **Default**: 1.0e-6


[back to top](#full-list-of-input-keywords)
6 changes: 6 additions & 0 deletions source/Makefile.Objects
Original file line number Diff line number Diff line change
Expand Up @@ -106,6 +106,7 @@ OBJS_MAIN=main.o\
driver_run.o\

OBJS_BASE=abfs-vector3_order.o\
assoc_laguerre.o\
complexarray.o\
complexmatrix.o\
clebsch_gordan_coeff.o\
Expand Down Expand Up @@ -326,6 +327,8 @@ OBJS_ORBITAL=ORB_atomic.o\
parallel_2d.o\
parallel_orbitals.o\
atomic_radials.o\
hydrogen_radials.o\
pswfc_radials.o\
beta_radials.o\
numerical_radial.o\
radial_collection.o\
Expand Down Expand Up @@ -419,6 +422,8 @@ OBJS_IO=input.o\
restart.o\
binstream.o\
to_wannier90.o\
to_qo.o\
to_qo_tools.o\
to_wannier90_pw.o\
to_wannier90_lcao_in_pw.o\
to_wannier90_lcao.o\
Expand All @@ -436,6 +441,7 @@ OBJS_IO=input.o\
output_rho.o\
output_potential.o\
output_mat_sparse.o\
output_radial.o\

OBJS_IO_LCAO=cal_r_overlap_R.o\
write_orb_info.o\
Expand Down
1 change: 1 addition & 0 deletions source/module_base/CMakeLists.txt
Original file line number Diff line number Diff line change
Expand Up @@ -10,6 +10,7 @@ endif()
add_library(
base
OBJECT
assoc_laguerre.cpp
clebsch_gordan_coeff.cpp
complexarray.cpp
complexmatrix.cpp
Expand Down
137 changes: 137 additions & 0 deletions source/module_base/assoc_laguerre.cpp
Original file line number Diff line number Diff line change
@@ -0,0 +1,137 @@
#include "module_base/assoc_laguerre.h"
#include "module_base/global_function.h"
//#include <tr1/cmath> // use cmath the factorial function
#include <cmath>
Assoc_Laguerre::Assoc_Laguerre()
{
}

Assoc_Laguerre::~Assoc_Laguerre()
{
}

void Assoc_Laguerre::generate(const int &n, const int &l, const double ns, double* const &s, double* L)
{
for(int i = 0; i < ns; i++)
{
L[i] = this->value(n, l, s[i]);
}
}

void Assoc_Laguerre::generate(const int &n, const int &l, std::vector<double> &x, std::vector<double> &y)
{
for(int i = 0; i < x.size(); i++)
{
y[i] = this->value(n, l, x[i]);
}
}

double Assoc_Laguerre::laguerre(const int &n, const double x)
{
if(n == 0)
{
return 1;
}
else if(n == 1)
{
return -x + 1;
}
else if(n == 2)
{
return 0.5 * x * x - 2 * x + 1;
}
else if(n == 3)
{
return -x * x * x / 6.0 + 3.0 * x * x / 2.0 - 3.0 * x + 1;
}
else if(n >= 4)
{
double n_ = static_cast<double>(n);
double first = (2*n_ - 1 - x)/n_ * Assoc_Laguerre::laguerre(n-1, x);
double second = (n_ - 1)/n_ * Assoc_Laguerre::laguerre(n-2, x);
return first - second;
}
else
{
ModuleBase::WARNING_QUIT("Assoc_Laguerre::laguerre", "n is out of range");
return 0;
}
}

double Assoc_Laguerre::associate_laguerre(const int &n, const double x, const int &a)
{
// formula from https://en.wikipedia.org/wiki/Laguerre_polynomials
double n_ = static_cast<double>(n);
double a_ = static_cast<double>(a);
if(n == 0)
{
return 1;
}
else if(n == 1)
{
return -x + 1 + a_;
}
else if(n == 2)
{
return 0.5 * (x*x - 2*(a_+2)*x + (a_+1)*(a_+2));
}
else if(n == 3)
{
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;
}
else if(n >= 4)
{
double first = (2*n_ - 1 + a_ - x)/n_ * this->associate_laguerre(n-1, x, a);
double second = (n_ + a_ - 1)/n_ * this->associate_laguerre(n-2, x, a);
return first - second;
}
else
{
ModuleBase::WARNING_QUIT("Assoc_Laguerre::associate_laguerre", "n is out of range");
return 0;
}
}

int Assoc_Laguerre::factorial(const int &n)
{
if(n == 0)
{
return 1;
}
else if(n > 0)
{
return n * this->factorial(n-1);
}
else
{
ModuleBase::WARNING_QUIT("Assoc_Laguerre::factorial", "n is out of range");
return 0;
}
}

double Assoc_Laguerre::value(const int &n, const int &l, const double &s)
{
int k_ = 2*l + 1;
int n_ = n - l - 1;
if(k_ < 0)
{
ModuleBase::WARNING_QUIT("Assoc_Laguerre::value", "k is out of range");
return 0;
}
if(n_ < 0)
{
ModuleBase::WARNING_QUIT("Assoc_Laguerre::value", "n is out of range");
return 0;
}
double L = 0;
for(int iq = 0; iq <= n_; iq++)
{
L += std::pow(-s, iq) *
static_cast<double>(this->factorial(n_ + k_)) /
static_cast<double>(this->factorial(n_ - iq)) /
static_cast<double>(this->factorial(k_ + iq)) /
static_cast<double>(this->factorial(iq));
}
//L = std::tr1::assoc_laguerre(n_, k_, s); // use standard library
return L;
}
47 changes: 47 additions & 0 deletions source/module_base/assoc_laguerre.h
Original file line number Diff line number Diff line change
@@ -0,0 +1,47 @@
#ifndef ASSOC_LAGUEERRE_H
#define ASSOC_LAGUEERRE_H

#include <vector>
#include <string>

class Assoc_Laguerre
{
public:
Assoc_Laguerre();
~Assoc_Laguerre();
/// @brief generate the associated Laguerre polynomial (overloaded for double*)
/// @param n principal quantum number
/// @param l orbital quantum number
/// @param ns number of x-coordinates
/// @param s x-coordinates
/// @param L y-coordinates
void generate(const int &n, const int &l, const double ns, double* const &s, double* L);
/// @brief generate the associated Laguerre polynomial (overloaded for std::vector)
/// @param n principal quantum number
/// @param l orbital quantum number
/// @param x x-coordinates in std::vector
/// @param y y-coordinates in std::vector
void generate(const int &n, const int &l, std::vector<double> &x, std::vector<double> &y);
/// @brief Laguerre polynomial
/// @param n degree of the polynomial
/// @param x radial coordinate
/// @return L_n(x)
double laguerre(const int &n, const double x);
/// @brief recursive relationship to find the associated Laguerre polynomial
/// @param n degree of the polynomial
/// @param x radial coordinate
/// @param a order of the polynomial
/// @return L^(a)_n(x)
double associate_laguerre(const int &n, const double x, const int &a);
/// @brief wrapper for associate_laguerre
/// @param n principal quantum number
/// @param l orbital quantum number
/// @param s radial coordinate
/// @return L^(2l+1)_(n-l-1)(s)
double value(const int &n, const int &l, const double &s);
/// @brief factorial function
/// @param n
/// @return n!
int factorial(const int &n);
};
#endif // ASSOC_LAGUEERRE_H
18 changes: 17 additions & 1 deletion source/module_base/atom_in.h
Original file line number Diff line number Diff line change
Expand Up @@ -91,7 +91,23 @@ class atom_in
{"Meitnerium", 109}, {"Darmstadtium", 110}, {"Roentgenium", 111}, {"Copernicium", 112}, {"Nihonium", 113}, {"Flerovium", 114},
{"Moscovium", 115}, {"Livermorium", 116}, {"Tennessine", 117}, {"Oganesson", 118}
};

std::map<std::string, int> principle_quantum_number
= {
{"H", 1}, {"He", 1}, {"Li", 2}, {"Be", 2}, {"B", 2}, {"C", 2}, {"N", 2}, {"O", 2}, {"F", 2},
{"Ne", 1}, {"Na", 2}, {"Mg", 3}, {"Al", 3}, {"Si", 3}, {"P", 3}, {"S", 3}, {"Cl", 3}, {"Ar", 2},
{"K", 3}, {"Ca", 4}, {"Sc", 4}, {"Ti", 4}, {"V", 4}, {"Cr", 4}, {"Mn", 4}, {"Fe", 4}, {"Co", 4},
{"Ni", 4}, {"Cu", 4}, {"Zn", 4}, {"Ga", 4}, {"Ge", 4}, {"As", 4}, {"Se", 4}, {"Br", 4}, {"Kr", 3},
{"Rb", 4}, {"Sr", 5}, {"Y", 5}, {"Zr", 5}, {"Nb", 5}, {"Mo", 5}, {"Tc", 5}, {"Ru", 5}, {"Rh", 5},
{"Pd", 5}, {"Ag", 5}, {"Cd", 5}, {"In", 5}, {"Sn", 5}, {"Sb", 5}, {"Te", 5}, {"I", 5}, {"Xe", 4},
{"Cs", 5}, {"Ba", 6}, {"La", 6}, {"Ce", 6}, {"Pr", 6}, {"Nd", 6}, {"Pm", 6}, {"Sm", 6}, {"Eu", 6},
{"Gd", 6}, {"Tb", 6}, {"Dy", 6}, {"Ho", 6}, {"Er", 6}, {"Tm", 6}, {"Yb", 6}, {"Lu", 6}, {"Hf", 6},
{"Ta", 6}, {"W", 6}, {"Re", 6}, {"Os", 6}, {"Ir", 6}, {"Pt", 6}, {"Au", 6}, {"Hg", 6}, {"Tl", 6},
{"Pb", 6}, {"Bi", 6}, {"Po", 6}, {"At", 6}, {"Rn", 6}, {"Fr", 7}, {"Ra", 7}, {"Ac", 7}, {"Th", 7},
{"Pa", 7}, {"U", 7}, {"Np", 7}, {"Pu", 7}, {"Am", 7}, {"Cm", 7}, {"Bk", 7}, {"Cf", 7}, {"Es", 7},
{"Fm", 7}, {"Md", 7}, {"No", 7}, {"Lr", 7}, {"Rf", 7}, {"Db", 7}, {"Sg", 7}, {"Bh", 7}, {"Hs", 7},
{"Mt", 7}, {"Ds", 7}, {"Rg", 7}, {"Cn", 7}, {"Nh", 7}, {"Fl", 7}, {"Mc", 7}, {"Lv", 7}, {"Ts", 7},
{"Og", 7}
};
};

#endif
20 changes: 10 additions & 10 deletions source/module_base/formatter_contextfmt.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -82,9 +82,9 @@ void formatter::ContextFmt::set_context(std::vector<std::string> phys_fmt) {
}
}

template<> formatter::ContextFmt& formatter::ContextFmt::operator<< <double>(double const&);
template<> formatter::ContextFmt& formatter::ContextFmt::operator<< <int>(int const&);
template<> formatter::ContextFmt& formatter::ContextFmt::operator<< <char>(char const&);
template formatter::ContextFmt& formatter::ContextFmt::operator<< <double>(double const&);
template formatter::ContextFmt& formatter::ContextFmt::operator<< <int>(int const&);
template formatter::ContextFmt& formatter::ContextFmt::operator<< <char>(char const&);

template <>
formatter::ContextFmt& formatter::ContextFmt::operator<<(const std::string& value) {
Expand Down Expand Up @@ -122,14 +122,14 @@ formatter::ContextFmt& formatter::ContextFmt::operator<< (char const* value)
return *this << value_;
}

template<> formatter::ContextFmt& formatter::ContextFmt::operator<< <double>(std::vector<double> const& value);
template<> formatter::ContextFmt& formatter::ContextFmt::operator<< <std::string>(std::vector<std::string> const& value);
template<> formatter::ContextFmt& formatter::ContextFmt::operator<< <int>(std::vector<int> const& value);
template formatter::ContextFmt& formatter::ContextFmt::operator<< <double>(std::vector<double> const& value);
template formatter::ContextFmt& formatter::ContextFmt::operator<< <std::string>(std::vector<std::string> const& value);
template formatter::ContextFmt& formatter::ContextFmt::operator<< <int>(std::vector<int> const& value);

template<> formatter::ContextFmt& formatter::ContextFmt::operator<< <double>(double*& value);
template<> formatter::ContextFmt& formatter::ContextFmt::operator<< <int>(int*& value);
template<> formatter::ContextFmt& formatter::ContextFmt::operator<< <std::string>(std::string*& value);
template<> formatter::ContextFmt& formatter::ContextFmt::operator<< <char>(char*& value);
template formatter::ContextFmt& formatter::ContextFmt::operator<< <double>(double*& value);
template formatter::ContextFmt& formatter::ContextFmt::operator<< <int>(int*& value);
template formatter::ContextFmt& formatter::ContextFmt::operator<< <std::string>(std::string*& value);
template formatter::ContextFmt& formatter::ContextFmt::operator<< <char>(char*& value);

void formatter::ContextFmt::reset() {
Table::reset();
Expand Down
9 changes: 9 additions & 0 deletions source/module_base/global_variable.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -289,4 +289,13 @@ int sc_scf_nmin = 2;
double alpha_trial = 0.01; // eV/uB^2
double sccut = 3; // eV/uB
std::string sc_file = "none";

//==========================================================
// Quasiatomic orbital related
//==========================================================
bool qo_switch = false;
std::string qo_basis = "hydrogen";
std::string qo_strategy = "minimal";
double qo_thr = 1.0e-6;
std::vector<double> qo_screening_coeff = {};
} // namespace GlobalV
Loading

0 comments on commit cbaccdf

Please sign in to comment.