From 18a495ac5b1a652ff36ae9b51959c99d56676927 Mon Sep 17 00:00:00 2001 From: jinzx10 Date: Tue, 6 Jun 2023 17:31:54 +0800 Subject: [PATCH 01/13] remove unnecessary "this" --- .../module_nao/numerical_radial.cpp | 78 ++++++++++--------- 1 file changed, 42 insertions(+), 36 deletions(-) diff --git a/source/module_basis/module_nao/numerical_radial.cpp b/source/module_basis/module_nao/numerical_radial.cpp index 9e0b708ef0..87b42549f0 100644 --- a/source/module_basis/module_nao/numerical_radial.cpp +++ b/source/module_basis/module_nao/numerical_radial.cpp @@ -11,8 +11,8 @@ #include "module_base/constants.h" #include "module_base/cubic_spline.h" -#include "module_base/spherical_bessel_transformer.h" #include "module_base/math_integral.h" +#include "module_base/spherical_bessel_transformer.h" using ModuleBase::PI; @@ -26,45 +26,45 @@ NumericalRadial::NumericalRadial() NumericalRadial::NumericalRadial(const NumericalRadial& other) { - this->symbol_ = other.symbol_; - this->itype_ = other.itype_; - this->izeta_ = other.izeta_; - this->l_ = other.l_; + symbol_ = other.symbol_; + itype_ = other.itype_; + izeta_ = other.izeta_; + l_ = other.l_; - this->nr_ = other.nr_; - this->nk_ = other.nk_; + nr_ = other.nr_; + nk_ = other.nk_; - this->is_fft_compliant_ = other.is_fft_compliant_; + is_fft_compliant_ = other.is_fft_compliant_; - this->pr_ = other.pr_; - this->pk_ = other.pk_; + pr_ = other.pr_; + pk_ = other.pk_; - this->use_internal_transformer_ = other.use_internal_transformer_; + use_internal_transformer_ = other.use_internal_transformer_; // deep copy if (other.ptr_rgrid()) { - this->rgrid_ = new double[nr_]; - this->rvalue_ = new double[nr_]; - std::memcpy(this->rgrid_, other.rgrid_, nr_ * sizeof(double)); - std::memcpy(this->rvalue_, other.rvalue_, nr_ * sizeof(double)); + rgrid_ = new double[nr_]; + rvalue_ = new double[nr_]; + std::memcpy(rgrid_, other.rgrid_, nr_ * sizeof(double)); + std::memcpy(rvalue_, other.rvalue_, nr_ * sizeof(double)); } if (other.ptr_kgrid()) { - this->kgrid_ = new double[nk_]; - this->kvalue_ = new double[nk_]; - std::memcpy(this->kgrid_, other.kgrid_, nk_ * sizeof(double)); - std::memcpy(this->kvalue_, other.kvalue_, nk_ * sizeof(double)); + kgrid_ = new double[nk_]; + kvalue_ = new double[nk_]; + std::memcpy(kgrid_, other.kgrid_, nk_ * sizeof(double)); + std::memcpy(kvalue_, other.kvalue_, nk_ * sizeof(double)); } if (use_internal_transformer_) { - this->sbt_ = new ModuleBase::SphericalBesselTransformer; + sbt_ = new ModuleBase::SphericalBesselTransformer; } else { - this->sbt_ = other.sbt_; + sbt_ = other.sbt_; } } @@ -116,7 +116,8 @@ NumericalRadial& NumericalRadial::operator=(const NumericalRadial& rhs) if (rhs.use_internal_transformer_) { - if (!use_internal_transformer_) { + if (!use_internal_transformer_) + { sbt_ = new ModuleBase::SphericalBesselTransformer; } } @@ -422,12 +423,12 @@ void NumericalRadial::radtab(const char op, // currently only FFT-compliant grids are supported! // FFT-based transform requires that two NumericalRadial objects have exactly the same grid - assert(this->is_fft_compliant_ && ket.is_fft_compliant_); - assert(this->nr_ == ket.nr_); - assert(this->rcut() == ket.rcut()); + assert(is_fft_compliant_ && ket.is_fft_compliant_); + assert(nr_ == ket.nr_); + assert(rcut() == ket.rcut()); double* ktmp = new double[nk_]; - std::transform(this->kvalue_, this->kvalue_ + nk_, ket.kvalue_, ktmp, std::multiplies()); + std::transform(kvalue_, kvalue_ + nk_, ket.kvalue_, ktmp, std::multiplies()); int op_pk = 0; switch (op) @@ -445,14 +446,14 @@ void NumericalRadial::radtab(const char op, { // derivative of the radial table if (l == 0) { // j'_0(x) = -j_1(x) - sbt_->radrfft(1, nk_, this->kcut(), ktmp, table, this->pk_ + ket.pk_ + op_pk - 1); + sbt_->radrfft(1, nk_, kcut(), ktmp, table, pk_ + ket.pk_ + op_pk - 1); std::for_each(table, table + nr_, [](double& x) { x *= -1; }); } else { // (2*l+1) * j'_l(x) = l * j_{l-1}(x) - (l+1) * j_{l+1}(x) double* rtmp = new double[nr_]; - sbt_->radrfft(l + 1, nk_, this->kcut(), ktmp, table, this->pk_ + ket.pk_ + op_pk - 1); - sbt_->radrfft(l - 1, nk_, this->kcut(), ktmp, rtmp, this->pk_ + ket.pk_ + op_pk - 1); + sbt_->radrfft(l + 1, nk_, kcut(), ktmp, table, pk_ + ket.pk_ + op_pk - 1); + sbt_->radrfft(l - 1, nk_, kcut(), ktmp, rtmp, pk_ + ket.pk_ + op_pk - 1); std::transform(table, table + nr_, rtmp, table, [l](double x, double y) { return (l * y - (l + 1) * x) / (2 * l + 1); }); @@ -461,33 +462,38 @@ void NumericalRadial::radtab(const char op, } else { // radial table - sbt_->radrfft(l, nk_, this->kcut(), ktmp, table, this->pk_ + ket.pk_ + op_pk); + sbt_->radrfft(l, nk_, kcut(), ktmp, table, pk_ + ket.pk_ + op_pk); } delete[] ktmp; } -void NumericalRadial::normalize(bool for_r_space) { +void NumericalRadial::normalize(bool for_r_space) +{ int& ngrid = for_r_space ? nr_ : nk_; + + // tbu stands for "to be updated" double*& grid_tbu = for_r_space ? rgrid_ : kgrid_; double*& value_tbu = for_r_space ? rvalue_ : kvalue_; - double factor = 0.0; + double factor = 0.0; double* integrand = new double[ngrid]; - double* rab = new double[ngrid-1]; + double* rab = new double[ngrid]; std::adjacent_difference(grid_tbu, grid_tbu + ngrid, rab); - std::transform(value_tbu, value_tbu+ngrid, grid_tbu, integrand, std::multiplies()); + std::transform(value_tbu, value_tbu + ngrid, grid_tbu, integrand, std::multiplies()); std::for_each(integrand, integrand + ngrid, [](double& x) { x *= x; }); + // FIXME Simpson_Integral should use only ngrid-1 rab points! + rab[ngrid - 1] = rab[ngrid - 2]; ModuleBase::Integral::Simpson_Integral(ngrid, integrand, rab, factor); - factor = 1./std::sqrt(factor); + factor = 1. / std::sqrt(factor); std::for_each(value_tbu, value_tbu + ngrid, [factor](double& x) { x *= factor; }); transform(for_r_space); delete[] rab; delete[] integrand; - //unit test TBD! + // unit test TBD! } void NumericalRadial::transform(const bool forward) From 37ad9d9eb74eb0defe10b5db45d2e1adc3f92df9 Mon Sep 17 00:00:00 2001 From: jinzx10 Date: Tue, 6 Jun 2023 17:42:36 +0800 Subject: [PATCH 02/13] allow copy of RadialSet & calculate rcut_max on demand --- source/module_basis/module_nao/radial_set.cpp | 85 ++++++++++++++++++- source/module_basis/module_nao/radial_set.h | 45 +++------- 2 files changed, 96 insertions(+), 34 deletions(-) diff --git a/source/module_basis/module_nao/radial_set.cpp b/source/module_basis/module_nao/radial_set.cpp index 4dab92210a..7be0f1885a 100644 --- a/source/module_basis/module_nao/radial_set.cpp +++ b/source/module_basis/module_nao/radial_set.cpp @@ -7,6 +7,88 @@ RadialSet::~RadialSet() delete[] index_map_; } +double RadialSet::rcut_max() const { + double rmax = 0.0; + for (int i = 0; i != nchi_; ++i) { + if (chi_[i].rcut() > rmax) { + rmax = chi_[i].rcut(); + } + } + return rmax; +} + +RadialSet::RadialSet(const RadialSet& other) { + symbol_ = other.symbol_; + itype_ = other.itype_; + lmax_ = other.lmax_; + + nzeta_ = nullptr; + if (lmax_ >= 0) { + nzeta_ = new int[lmax_ + 1]; + for (int l = 0; l <= lmax_; l++) + nzeta_[l] = other.nzeta_[l]; + } + nzeta_max_ = other.nzeta_max_; + nchi_ = other.nchi_; + + chi_ = nullptr; + if (nchi_ > 0) { + chi_ = new NumericalRadial[nchi_]; + for (int i = 0; i < nchi_; i++) { + chi_[i] = other.chi_[i]; // deep copy + } + } + + index_map_ = nullptr; + if (lmax_ >= 0 && nzeta_max_ > 0) { + index_map_ = new int[(lmax_ + 1) * nzeta_max_]; + for (int i = 0; i < (lmax_ + 1) * nzeta_max_; i++) { + index_map_[i] = other.index_map_[i]; + } + } +} + +RadialSet& RadialSet::operator=(const RadialSet& rhs) { + if (&rhs == this) { + return *this; + } + + symbol_ = rhs.symbol_; + itype_ = rhs.itype_; + lmax_ = rhs.lmax_; + + delete[] nzeta_; + nzeta_ = nullptr; + if (lmax_ >= 0) { + nzeta_ = new int[lmax_ + 1]; + for (int l = 0; l <= lmax_; l++) + nzeta_[l] = rhs.nzeta_[l]; + } + nzeta_max_ = rhs.nzeta_max_; + nchi_ = rhs.nchi_; + + //rcut_max_ = rhs.rcut_max_; + + delete[] chi_; + chi_ = nullptr; + if (nchi_ > 0) { + chi_ = new NumericalRadial[nchi_]; + for (int i = 0; i < nchi_; i++) { + chi_[i] = rhs.chi_[i]; // deep copy + } + } + + delete[] index_map_; + index_map_ = nullptr; + if (lmax_ >= 0 && nzeta_max_ > 0) { + index_map_ = new int[(lmax_ + 1) * nzeta_max_]; + for (int i = 0; i < (lmax_ + 1) * nzeta_max_; i++) { + index_map_[i] = rhs.index_map_[i]; + } + } + return *this; +} + int RadialSet::index(int l, int izeta) const { assert(l >= 0 && l <= lmax_); @@ -31,6 +113,7 @@ void RadialSet::set_grid(const bool for_r_space, const int ngrid, const double* { for (int i = 0; i < nchi_; i++) chi_[i].set_grid(for_r_space, ngrid, grid, mode); + //rcut_max_ = grid[ngrid-1]; } void RadialSet::set_uniform_grid(const bool for_r_space, const int ngrid, const double cutoff, const char mode) @@ -50,7 +133,7 @@ void RadialSet::cleanup() nzeta_max_ = 0; nchi_ = 0; - rcut_max_ = 0; + //rcut_max_ = 0; delete[] chi_; chi_ = nullptr; diff --git a/source/module_basis/module_nao/radial_set.h b/source/module_basis/module_nao/radial_set.h index acc28e9c3e..ca4b8f2669 100644 --- a/source/module_basis/module_nao/radial_set.h +++ b/source/module_basis/module_nao/radial_set.h @@ -22,7 +22,10 @@ class RadialSet { public: RadialSet(){}; - ~RadialSet(); + RadialSet(const RadialSet&); //!< deep copy + RadialSet& operator=(const RadialSet&); //!< deep copy + + virtual ~RadialSet(); //! build the set of numerical radial functions from a file virtual void build(const std::string& file, //!< orbital or pseudopotential file @@ -37,36 +40,14 @@ class RadialSet * Get access to private members. * */ //!@{ - const std::string& symbol() const - { - return symbol_; - } - int itype() const - { - return itype_; - } - int lmax() const - { - return lmax_; - } - - int nzeta(int l) const - { - return nzeta_[l]; - } - int nzeta_max() const - { - return nzeta_max_; - } - int nchi() const - { - return nchi_; - } - - double rcut_max() const - { - return rcut_max_; - } + const std::string& symbol() const { return symbol_; } + int itype() const { return itype_; } + int lmax() const { return lmax_; } + double rcut_max() const; + + int nzeta(int l) const { return nzeta_[l]; } + int nzeta_max() const { return nzeta_max_; } + int nchi() const { return nchi_; } const NumericalRadial& chi(int l, int izeta); //!@} @@ -98,8 +79,6 @@ class RadialSet int nzeta_max_ = 0; //!< maximum number of NumericalRadial objects among each angular momentum int nchi_ = 0; //!< total number of NumericalRadial objects - double rcut_max_ = 0; //!< maximum r-space cutoff radius among all NumericalRadial objects - NumericalRadial* chi_ = nullptr; //!< array of NumericalRadial objects int* index_map_ = nullptr; //!< map (l,izeta) to an index in chi_ array From 4926b19554f47f903d7811aa70fa7869d2434a1c Mon Sep 17 00:00:00 2001 From: jinzx10 Date: Tue, 6 Jun 2023 17:46:01 +0800 Subject: [PATCH 03/13] clean up --- source/module_basis/module_nao/atomic_radials.cpp | 14 +++++--------- source/module_basis/module_nao/atomic_radials.h | 11 +++-------- source/module_basis/module_nao/numerical_radial.h | 9 ++++----- .../module_nao/test/atomic_radials_test.cpp | 2 +- .../module_nao/test/numerical_radial_test.cpp | 1 - 5 files changed, 13 insertions(+), 24 deletions(-) diff --git a/source/module_basis/module_nao/atomic_radials.cpp b/source/module_basis/module_nao/atomic_radials.cpp index 73e0bb9945..e2f5f019d0 100644 --- a/source/module_basis/module_nao/atomic_radials.cpp +++ b/source/module_basis/module_nao/atomic_radials.cpp @@ -45,6 +45,11 @@ void AtomicRadials::build(const std::string& file, const int itype, std::ofstrea itype_ = itype; read_abacus_orb(ifs, ptr_log, rank); + + if (rank == 0) + { + ifs.close(); + } } void AtomicRadials::read_abacus_orb(std::ifstream& ifs, std::ofstream* ptr_log, const int rank) @@ -82,10 +87,6 @@ void AtomicRadials::read_abacus_orb(std::ifstream& ifs, std::ofstream* ptr_log, { ifs >> orb_ecut_; } - else if (tmp == "Cutoff(a.u.)") - { - ifs >> rcut_max_; // orbitals in an ABACUS orbital file have the same cutoff radius - } else if (tmp == "Lmax") { ifs >> lmax_; @@ -212,11 +213,6 @@ void AtomicRadials::read_abacus_orb(std::ifstream& ifs, std::ofstream* ptr_log, chi_[index(l, izeta)].normalize(); } - if (rank == 0) - { - ifs.close(); - } - delete[] is_read; delete[] rvalue; delete[] rgrid; diff --git a/source/module_basis/module_nao/atomic_radials.h b/source/module_basis/module_nao/atomic_radials.h index 84a4e87213..62afc3faaa 100644 --- a/source/module_basis/module_nao/atomic_radials.h +++ b/source/module_basis/module_nao/atomic_radials.h @@ -25,10 +25,8 @@ class AtomicRadials : public RadialSet { public: - AtomicRadials(){}; - ~AtomicRadials() - { - } // ~RadialSet() is called automatically + AtomicRadials() {} + ~AtomicRadials() {} // ~RadialSet() is called automatically //! Build the class from an orbital file void build(const std::string& file, //!< orbital file name @@ -38,10 +36,7 @@ class AtomicRadials : public RadialSet ); //! Get the energy cutoff as given by the orbital file - double orb_ecut() const - { - return orb_ecut_; - } + double orb_ecut() const { return orb_ecut_; } private: double orb_ecut_; //!< energy cutoff as given by the orbital file diff --git a/source/module_basis/module_nao/numerical_radial.h b/source/module_basis/module_nao/numerical_radial.h index 88dcad87e0..f539584f89 100644 --- a/source/module_basis/module_nao/numerical_radial.h +++ b/source/module_basis/module_nao/numerical_radial.h @@ -153,9 +153,6 @@ class NumericalRadial //! Saves the data to file (what data, in what format?) void save(std::string file = "" /*! file name */) const; - //! Set itype for the object. - void set_itype(const int itype) { itype_ = itype; } - //! Computes the radial table for two-center integrals. /*! * TODO add support for non-FFT-compliant grid @@ -249,13 +246,15 @@ class NumericalRadial //! gets r-space grid cutoff distance double rcut() const { - return rgrid_ ? rgrid_[nr_ - 1] : 0.0; + assert(rgrid_); + return rgrid_[nr_ - 1]; } //! gets k-space grid cutoff distance double kcut() const { - return kgrid_ ? kgrid_[nk_ - 1] : 0.0; + assert(kgrid_); + return kgrid_[nk_ - 1]; } //! gets the pointer to r-space grid points diff --git a/source/module_basis/module_nao/test/atomic_radials_test.cpp b/source/module_basis/module_nao/test/atomic_radials_test.cpp index f215b2c2a2..0fdfc048cb 100644 --- a/source/module_basis/module_nao/test/atomic_radials_test.cpp +++ b/source/module_basis/module_nao/test/atomic_radials_test.cpp @@ -14,7 +14,7 @@ /*! * Tested functions: * - * - read + * - build * - parse an orbital file and initialize the NumericalRadial objects * * - all "getters" diff --git a/source/module_basis/module_nao/test/numerical_radial_test.cpp b/source/module_basis/module_nao/test/numerical_radial_test.cpp index 2cbac16380..429b2bc38a 100644 --- a/source/module_basis/module_nao/test/numerical_radial_test.cpp +++ b/source/module_basis/module_nao/test/numerical_radial_test.cpp @@ -201,7 +201,6 @@ TEST_F(NumericalRadialTest, BuildAndGet) EXPECT_EQ(chi.nr(), sz); EXPECT_EQ(chi.nk(), 0); EXPECT_EQ(chi.rcut(), grid[sz - 1]); - EXPECT_EQ(chi.kcut(), 0); ASSERT_NE(chi.ptr_rgrid(), nullptr); ASSERT_NE(chi.ptr_rvalue(), nullptr); From 4f4453975975041552128889d24abb03bab78e87 Mon Sep 17 00:00:00 2001 From: jinzx10 Date: Tue, 6 Jun 2023 18:02:55 +0800 Subject: [PATCH 04/13] add "RadialCollection", the top-level container for radial functions --- source/module_basis/module_nao/CMakeLists.txt | 1 + .../module_nao/radial_collection.cpp | 95 ++++++++++++ .../module_nao/radial_collection.h | 95 ++++++++++++ .../module_nao/test/CMakeLists.txt | 13 +- .../test/radial_collection_test.cpp | 142 ++++++++++++++++++ 5 files changed, 345 insertions(+), 1 deletion(-) create mode 100644 source/module_basis/module_nao/radial_collection.cpp create mode 100644 source/module_basis/module_nao/radial_collection.h create mode 100644 source/module_basis/module_nao/test/radial_collection_test.cpp diff --git a/source/module_basis/module_nao/CMakeLists.txt b/source/module_basis/module_nao/CMakeLists.txt index 1268400eac..914adfeb69 100644 --- a/source/module_basis/module_nao/CMakeLists.txt +++ b/source/module_basis/module_nao/CMakeLists.txt @@ -5,6 +5,7 @@ if(ENABLE_LCAO) numerical_radial.cpp radial_set.cpp atomic_radials.cpp + radial_collection.cpp ) if(ENABLE_COVERAGE) diff --git a/source/module_basis/module_nao/radial_collection.cpp b/source/module_basis/module_nao/radial_collection.cpp new file mode 100644 index 0000000000..d237ce8bf6 --- /dev/null +++ b/source/module_basis/module_nao/radial_collection.cpp @@ -0,0 +1,95 @@ +#include "module_basis/module_nao/radial_collection.h" + +#include "module_base/spherical_bessel_transformer.h" +#include "module_basis/module_nao/atomic_radials.h" +#include "module_basis/module_nao/beta_radials.h" + +double RadialCollection::rcut_max() const +{ + double rmax = 0.0; + for (int itype = 0; itype < ntype_; ++itype) + { + rmax = std::max(rmax, radset_[itype]->rcut_max()); + } + return rmax; +} + +RadialCollection::~RadialCollection() +{ + for (int itype = 0; itype < ntype_; ++itype) + { + delete radset_[itype]; + } + delete[] radset_; +} + +void RadialCollection::cleanup() +{ + for (int itype = 0; itype < ntype_; ++itype) + { + delete radset_[itype]; + } + delete[] radset_; + radset_ = nullptr; + ntype_ = 0; + lmax_ = -1; + nchi_ = 0; +} + +void RadialCollection::build(const int nfile, const std::string* const file, const char type) +{ + assert(type == 'o' || type == 'p'); + + cleanup(); + + ntype_ = nfile; + radset_ = new RadialSet*[ntype_]; + switch (type) + { + case 'o': + for (int itype = 0; itype < ntype_; ++itype) + { + radset_[itype] = new AtomicRadials; + radset_[itype]->build(file[itype], itype); + } + break; + case 'p': + for (int itype = 0; itype < ntype_; ++itype) + { + radset_[itype] = new BetaRadials; + radset_[itype]->build(file[itype], itype); + } + break; + default:; /* not supposed to happen */ + } + + for (int itype = 0; itype < ntype_; ++itype) + { + lmax_ = std::max(lmax_, radset_[itype]->lmax()); + nchi_ += radset_[itype]->nchi(); + } +} + +void RadialCollection::set_transformer(ModuleBase::SphericalBesselTransformer* sbt, int update) +{ + for (int itype = 0; itype < ntype_; ++itype) + { + radset_[itype]->set_transformer(sbt, update); + } +} + +void RadialCollection::set_grid(const bool for_r_space, const int ngrid, const double* grid, const char mode) +{ + for (int itype = 0; itype < ntype_; ++itype) + { + radset_[itype]->set_grid(for_r_space, ngrid, grid, mode); + } +} + +void RadialCollection::set_uniform_grid(const bool for_r_space, const int ngrid, const double cutoff, const char mode) +{ + for (int itype = 0; itype < ntype_; ++itype) + { + radset_[itype]->set_uniform_grid(for_r_space, ngrid, cutoff, mode); + } +} diff --git a/source/module_basis/module_nao/radial_collection.h b/source/module_basis/module_nao/radial_collection.h new file mode 100644 index 0000000000..308a7e6277 --- /dev/null +++ b/source/module_basis/module_nao/radial_collection.h @@ -0,0 +1,95 @@ +#ifndef RADIAL_COLLECTION_H_ +#define RADIAL_COLLECTION_H_ + +#include + +#include "module_basis/module_nao/radial_set.h" + +//! A class that holds the collection of all numerical radial functions of a kind. +/*! + * This class is supposed to hold the radial functions of all numerical atomic + * orbitals or Kleinman-Bylander beta functions from all elements involved in + * calculation. The way it stores the radial functions is type(element)-wise, + * i.e. all radial functions of the same itype (element) are stored in a + * RadialSet object, and the collection holds an array of such objects. + * */ +class RadialCollection +{ + public: + RadialCollection(){}; + ~RadialCollection(); + + void build(const int nfile, const std::string* const file, const char type = 'o'); + + /*! @name Getters + * + * Get access to private members (and their properties). + * */ + //!@{ + //! element symbol of a given type + const std::string& symbol(int itype) const { return radset_[itype]->symbol(); } + + //! number of RadialSet objects in the collection + int ntype() const { return ntype_; } + + //! maximum angular momentum of all NumericalRadial objects in the collection + int lmax() const { return lmax_; } + + //! maximum cutoff radius of all NumericalRadial objects in the collection + double rcut_max() const; + + //! number of distinct radial functions of a given type and angular momentum + int nzeta(int itype, int l) const { return radset_[itype]->nzeta(l); } + + //! maximum number of distinct radial functions of a given type among all angular momentum + int nzeta_max(int itype) const { return radset_[itype]->nzeta_max(); } + + //! total number of NumericalRadial objects in the collection + int nchi() const { return nchi_; } + + //! number of NumericalRadial objects of a given type + int nchi(int itype) const { return radset_[itype]->nchi(); } + + //! get access to the NumericalRadial object with given type, angular momentum and zeta number + const NumericalRadial& operator()(const int itype, const int l, const int izeta) const + { + assert(itype >= 0 && itype < ntype_); + return radset_[itype]->chi(l, izeta); + } + + //! get access to the RadialSet object with given type + const RadialSet& operator()(const int itype) const + { + assert(itype >= 0 && itype < ntype_); + return *radset_[itype]; + } + //!@} + + /*! @name property setters for all RadialSet objects + * + * @see RadialSet + * */ + //!@{ + //! Set a spherical Bessel transformers for all RadialSet objects + //! @see RadialSet::set_transformer + void set_transformer(ModuleBase::SphericalBesselTransformer* sbt = nullptr, int update = 0); + + //! Set a common grid for all RadialSet objects + //! @see RadialSet::set_grid + void set_grid(const bool for_r_space, const int ngrid, const double* grid, const char mode = 'i'); + + //! Set a common uniform grid for all RadialSet objects + //! @see RadialSet::set_uniform_grid + void set_uniform_grid(const bool for_r_space, const int ngrid, const double cutoff, const char mode = 'i'); + //!@} + + private: + int ntype_ = 0; //!< number of RadialSet in the collection + int lmax_ = -1; //!< maximum angular momentum of all NumericalRadial objects in the collection + int nchi_ = 0; //!< total number of NumericalRadial objects in the collection + RadialSet** radset_ = nullptr; + + void cleanup(); +}; + +#endif diff --git a/source/module_basis/module_nao/test/CMakeLists.txt b/source/module_basis/module_nao/test/CMakeLists.txt index 108abc0c8b..c3effe58ba 100644 --- a/source/module_basis/module_nao/test/CMakeLists.txt +++ b/source/module_basis/module_nao/test/CMakeLists.txt @@ -10,7 +10,7 @@ AddTest( AddTest( TARGET atomic_radials - SOURCES + SOURCES atomic_radials_test.cpp ../atomic_radials.cpp ../radial_set.cpp @@ -18,3 +18,14 @@ AddTest( LIBS ${math_libs} device base ) +AddTest( + TARGET radial_collection + SOURCES + radial_collection_test.cpp + ../radial_collection.cpp + ../atomic_radials.cpp + ../beta_radials.cpp + ../radial_set.cpp + ../numerical_radial.cpp + LIBS ${math_libs} device base +) diff --git a/source/module_basis/module_nao/test/radial_collection_test.cpp b/source/module_basis/module_nao/test/radial_collection_test.cpp new file mode 100644 index 0000000000..f97fbf987e --- /dev/null +++ b/source/module_basis/module_nao/test/radial_collection_test.cpp @@ -0,0 +1,142 @@ +#include "module_basis/module_nao/radial_collection.h" +#include "gtest/gtest.h" + +#ifdef __MPI +#include +#endif + +/*********************************************************** + * Unit test of class "RadialCollection" + ***********************************************************/ +/*! + * Tested functions: + * + * - build + * - parse an orbital file and initialize the NumericalRadial objects + * + * - all "getters" + * - Get access to private members. + * + * - all "batch setters" + * - Set a property for all RadialSet objects at once + * */ +class RadialCollectionTest : public ::testing::Test +{ + protected: + void SetUp(); + void TearDown(); + + RadialCollection orb; //!< object under test + int nfile = 0; // number of orbital files + std::string* file = nullptr; //!< orbitals file to read from + std::string log_file = "./test_files/atomic_orbital.log"; //!< file for logging +}; + +void RadialCollectionTest::SetUp() { + nfile = 4; + file = new std::string[nfile]; + file[0] = "../../../../../tests/PP_ORB/C_gga_8au_100Ry_2s2p1d.orb"; + file[1] = "../../../../../tests/PP_ORB/H_gga_8au_60Ry_2s1p.orb"; + file[2] = "../../../../../tests/PP_ORB/O_gga_10au_100Ry_2s2p1d.orb"; + file[3] = "../../../../../tests/PP_ORB/Fe_gga_9au_100Ry_4s2p2d1f.orb"; +} + +void RadialCollectionTest::TearDown() { + delete[] file; +} + +TEST_F(RadialCollectionTest, BuildAndGet) { + orb.build(nfile, file, 'o'); + + EXPECT_EQ(orb.symbol(0), "C"); + EXPECT_EQ(orb.symbol(1), "H"); + EXPECT_EQ(orb.symbol(2), "O"); + EXPECT_EQ(orb.symbol(3), "Fe"); + + EXPECT_EQ(orb.ntype(), 4); + EXPECT_EQ(orb.lmax(), 3); + EXPECT_DOUBLE_EQ(orb.rcut_max(), 10.0); + + EXPECT_EQ(orb.nzeta(0,0), 2); + EXPECT_EQ(orb.nzeta(0,1), 2); + EXPECT_EQ(orb.nzeta(0,2), 1); + + EXPECT_EQ(orb.nzeta(1,0), 2); + EXPECT_EQ(orb.nzeta(1,1), 1); + + EXPECT_EQ(orb.nzeta(2,0), 2); + EXPECT_EQ(orb.nzeta(2,1), 2); + EXPECT_EQ(orb.nzeta(2,2), 1); + + EXPECT_EQ(orb.nzeta(3,0), 4); + EXPECT_EQ(orb.nzeta(3,1), 2); + EXPECT_EQ(orb.nzeta(3,2), 2); + EXPECT_EQ(orb.nzeta(3,3), 1); + + EXPECT_EQ(orb.nzeta_max(0), 2); + EXPECT_EQ(orb.nzeta_max(1), 2); + EXPECT_EQ(orb.nzeta_max(2), 2); + EXPECT_EQ(orb.nzeta_max(3), 4); + + EXPECT_EQ(orb.nchi(0), 5); + EXPECT_EQ(orb.nchi(1), 3); + EXPECT_EQ(orb.nchi(2), 5); + EXPECT_EQ(orb.nchi(3), 9); + EXPECT_EQ(orb.nchi(), 22); + + EXPECT_EQ(orb(0).itype(), 0); + EXPECT_EQ(orb(3).itype(), 3); + + EXPECT_EQ(orb(0, 1, 0).l(), 1); + EXPECT_EQ(orb(3, 3, 0).l(), 3); +} + +TEST_F(RadialCollectionTest, BatchSet) { + orb.build(nfile, file, 'o'); + + ModuleBase::SphericalBesselTransformer sbt; + orb.set_transformer(&sbt); + orb.set_uniform_grid(true, 2001, 20.0); + + EXPECT_EQ(orb.rcut_max(), 20.0); + for (int itype = 0; itype != orb.ntype(); ++itype) { + for (int l = 0; l <= orb.lmax(); ++l) { + for (int izeta = 0; izeta != orb.nzeta(itype, l); ++izeta) { + EXPECT_EQ(&sbt, orb(itype, l, izeta).ptr_sbt()); + EXPECT_DOUBLE_EQ(orb(itype, l, izeta).rcut(), 20.0); + } + } + } + + double* grid = new double[3]; + grid[0] = 0.0; + grid[1] = 1.0; + grid[2] = 3.14; + + orb.set_grid(true, 3, grid, 'i'); + for (int itype = 0; itype != orb.ntype(); ++itype) { + for (int l = 0; l <= orb.lmax(); ++l) { + for (int izeta = 0; izeta != orb.nzeta(itype, l); ++izeta) { + EXPECT_DOUBLE_EQ(orb(itype, l, izeta).rcut(), 3.14); + } + } + } + delete[] grid; +} + +int main(int argc, char** argv) +{ + +#ifdef __MPI + MPI_Init(&argc, &argv); +#endif + + testing::InitGoogleTest(&argc, argv); + int result = RUN_ALL_TESTS(); + +#ifdef __MPI + MPI_Finalize(); +#endif + + return result; +} From 36d595f92c78df45ab1f9b6753ed6601814c4b7a Mon Sep 17 00:00:00 2001 From: jinzx10 Date: Tue, 6 Jun 2023 21:04:33 +0800 Subject: [PATCH 05/13] enable MPI test for module_nao & bug fix for reading orbital file --- .../module_basis/module_nao/atomic_radials.cpp | 17 +++++++++++++++-- .../module_basis/module_nao/test/CMakeLists.txt | 2 -- .../module_nao/test/atomic_radials_test.cpp | 11 +++++++++-- 3 files changed, 24 insertions(+), 6 deletions(-) diff --git a/source/module_basis/module_nao/atomic_radials.cpp b/source/module_basis/module_nao/atomic_radials.cpp index e2f5f019d0..faa560c333 100644 --- a/source/module_basis/module_nao/atomic_radials.cpp +++ b/source/module_basis/module_nao/atomic_radials.cpp @@ -9,6 +9,7 @@ void AtomicRadials::build(const std::string& file, const int itype, std::ofstream* ptr_log, const int rank) { + // deallocates all arrays and reset variables cleanup(); std::ifstream ifs; @@ -21,6 +22,7 @@ void AtomicRadials::build(const std::string& file, const int itype, std::ofstrea } #ifdef __MPI + // NOTE Paralle_Common::bcast_xxx use GlobalV::MY_RANK. Use with care! Parallel_Common::bcast_bool(is_open); #endif @@ -144,19 +146,29 @@ void AtomicRadials::read_abacus_orb(std::ifstream& ifs, std::ofstream* ptr_log, } #ifdef __MPI + // NOTE Paralle_Common::bcast_xxx use GlobalV::MY_RANK. Use with care! Parallel_Common::bcast_string(symbol_); Parallel_Common::bcast_double(orb_ecut_); Parallel_Common::bcast_int(lmax_); - Parallel_Common::bcast_int(nzeta_, lmax_ + 1); Parallel_Common::bcast_int(nchi_); Parallel_Common::bcast_int(nzeta_max_); - Parallel_Common::bcast_int(index_map_, (lmax_ + 1) * nzeta_max_); Parallel_Common::bcast_int(ngrid); Parallel_Common::bcast_double(dr); #endif + if (rank != 0) { + nzeta_ = new int[lmax_ + 1]; + index_map_ = new int[(lmax_ + 1) * nzeta_max_]; + } + +#ifdef __MPI + // NOTE Paralle_Common::bcast_xxx use GlobalV::MY_RANK. Use with care! + Parallel_Common::bcast_int(nzeta_, lmax_ + 1); + Parallel_Common::bcast_int(index_map_, (lmax_ + 1) * nzeta_max_); +#endif + double* rvalue = new double[ngrid]; double* rgrid = new double[ngrid]; for (int ir = 0; ir != ngrid; ++ir) @@ -201,6 +213,7 @@ void AtomicRadials::read_abacus_orb(std::ifstream& ifs, std::ofstream* ptr_log, } #ifdef __MPI + // NOTE Paralle_Common::bcast_xxx use GlobalV::MY_RANK. Use with care! Parallel_Common::bcast_int(l); Parallel_Common::bcast_int(izeta); Parallel_Common::bcast_double(rvalue, ngrid); diff --git a/source/module_basis/module_nao/test/CMakeLists.txt b/source/module_basis/module_nao/test/CMakeLists.txt index c3effe58ba..56d16e6d3d 100644 --- a/source/module_basis/module_nao/test/CMakeLists.txt +++ b/source/module_basis/module_nao/test/CMakeLists.txt @@ -1,5 +1,3 @@ -remove_definitions(-D__MPI) - AddTest( TARGET numerical_radial SOURCES diff --git a/source/module_basis/module_nao/test/atomic_radials_test.cpp b/source/module_basis/module_nao/test/atomic_radials_test.cpp index 0fdfc048cb..d493dc7ed2 100644 --- a/source/module_basis/module_nao/test/atomic_radials_test.cpp +++ b/source/module_basis/module_nao/test/atomic_radials_test.cpp @@ -7,6 +7,7 @@ #endif #include "module_base/constants.h" +#include "module_base/global_variable.h" /*********************************************************** * Unit test of class "AtomicRadials" @@ -38,8 +39,11 @@ class AtomicRadialsTest : public ::testing::Test TEST_F(AtomicRadialsTest, ReadAndGet) { +#ifdef __MPI + MPI_Comm_rank(MPI_COMM_WORLD, &GlobalV::MY_RANK); +#endif - Ti_radials.build(file); + Ti_radials.build(file, 0, nullptr, GlobalV::MY_RANK); EXPECT_EQ(Ti_radials.lmax(), 3); EXPECT_EQ(Ti_radials.nzeta(0), 4); @@ -69,9 +73,12 @@ TEST_F(AtomicRadialsTest, ReadAndGet) TEST_F(AtomicRadialsTest, BatchSet) { +#ifdef __MPI + MPI_Comm_rank(MPI_COMM_WORLD, &GlobalV::MY_RANK); +#endif int itype = 5; - Ti_radials.build(file, itype); + Ti_radials.build(file, itype, nullptr, GlobalV::MY_RANK); EXPECT_EQ(Ti_radials.itype(), 5); EXPECT_EQ(Ti_radials.chi(0, 0).itype(), 5); From 825157afe15d13fb35dd46a659aa8f2c826cbaa7 Mon Sep 17 00:00:00 2001 From: jinzx10 Date: Wed, 7 Jun 2023 21:41:21 +0800 Subject: [PATCH 06/13] clean up --- .../module_nao/atomic_radials.cpp | 28 +----- .../module_basis/module_nao/atomic_radials.h | 4 +- .../module_nao/numerical_radial.h | 72 ++++------------ .../module_nao/radial_collection.cpp | 6 +- .../module_nao/radial_collection.h | 10 +++ source/module_basis/module_nao/radial_set.cpp | 85 ++++++++++++++----- source/module_basis/module_nao/radial_set.h | 22 +++-- .../test/radial_collection_test.cpp | 32 ++++--- 8 files changed, 130 insertions(+), 129 deletions(-) diff --git a/source/module_basis/module_nao/atomic_radials.cpp b/source/module_basis/module_nao/atomic_radials.cpp index faa560c333..5210770bb0 100644 --- a/source/module_basis/module_nao/atomic_radials.cpp +++ b/source/module_basis/module_nao/atomic_radials.cpp @@ -22,7 +22,6 @@ void AtomicRadials::build(const std::string& file, const int itype, std::ofstrea } #ifdef __MPI - // NOTE Paralle_Common::bcast_xxx use GlobalV::MY_RANK. Use with care! Parallel_Common::bcast_bool(is_open); #endif @@ -119,34 +118,14 @@ void AtomicRadials::read_abacus_orb(std::ifstream& ifs, std::ofstream* ptr_log, * 3. a map from (l, izeta) to 1-d array index in chi_ * */ nchi_ = 0; - nzeta_max_ = 0; for (int l = 0; l <= lmax_; ++l) { nchi_ += nzeta_[l]; - nzeta_max_ = std::max(nzeta_[l], nzeta_max_); - } - - index_map_ = new int[(lmax_ + 1) * nzeta_max_]; - int index_chi = 0; - for (int l = 0; l <= lmax_; ++l) - { - for (int izeta = 0; izeta != nzeta_max_; ++izeta) - { - if (izeta >= nzeta_[l]) - { - index_map_[l * nzeta_max_ + izeta] = -1; // -1 means no such orbital - } - else - { - index_map_[l * nzeta_max_ + izeta] = index_chi; - ++index_chi; - } - } } + indexing(); // calculate nzeta_max_ and build index_map_ } #ifdef __MPI - // NOTE Paralle_Common::bcast_xxx use GlobalV::MY_RANK. Use with care! Parallel_Common::bcast_string(symbol_); Parallel_Common::bcast_double(orb_ecut_); Parallel_Common::bcast_int(lmax_); @@ -158,13 +137,13 @@ void AtomicRadials::read_abacus_orb(std::ifstream& ifs, std::ofstream* ptr_log, Parallel_Common::bcast_double(dr); #endif - if (rank != 0) { + if (rank != 0) + { nzeta_ = new int[lmax_ + 1]; index_map_ = new int[(lmax_ + 1) * nzeta_max_]; } #ifdef __MPI - // NOTE Paralle_Common::bcast_xxx use GlobalV::MY_RANK. Use with care! Parallel_Common::bcast_int(nzeta_, lmax_ + 1); Parallel_Common::bcast_int(index_map_, (lmax_ + 1) * nzeta_max_); #endif @@ -213,7 +192,6 @@ void AtomicRadials::read_abacus_orb(std::ifstream& ifs, std::ofstream* ptr_log, } #ifdef __MPI - // NOTE Paralle_Common::bcast_xxx use GlobalV::MY_RANK. Use with care! Parallel_Common::bcast_int(l); Parallel_Common::bcast_int(izeta); Parallel_Common::bcast_double(rvalue, ngrid); diff --git a/source/module_basis/module_nao/atomic_radials.h b/source/module_basis/module_nao/atomic_radials.h index 62afc3faaa..b4aac608f5 100644 --- a/source/module_basis/module_nao/atomic_radials.h +++ b/source/module_basis/module_nao/atomic_radials.h @@ -15,11 +15,9 @@ * int element_index = 1; * std::ofstream ofs_log("/path/to/log/file"); * std::string orb_file = "/path/to/orbital/file"; - * int rank = 0; - * MPI_Comm_rank(MPI_COMM_WORLD, &rank); * * AtomicRadials O_radials; - * O_radials.build(orb_file, element_index, ofs_log, rank); + * O_radials.build(orb_file, element_index, ofs_log, GlobalV::MY_RANK); * * */ class AtomicRadials : public RadialSet diff --git a/source/module_basis/module_nao/numerical_radial.h b/source/module_basis/module_nao/numerical_radial.h index f539584f89..9beb32d86f 100644 --- a/source/module_basis/module_nao/numerical_radial.h +++ b/source/module_basis/module_nao/numerical_radial.h @@ -85,7 +85,7 @@ class NumericalRadial const int p = 0, //!< exponent of the implicit power term in input values, @see @ref group1 const int izeta = 0, //!< index for the multiplicity of radial functions of the same itype and l const std::string symbol = "", //!< chemical symbol - const int itype = 0 //!< index for the element in calculation + const int itype = 0 //!< index for the element in calculation ); //! Sets a SphericalBesselTransformer. @@ -208,40 +208,22 @@ class NumericalRadial * */ ///@{ //! gets symbol_ - std::string const& symbol() const - { - return symbol_; - } + std::string const& symbol() const { return symbol_; } //! gets itype_ - int itype() const - { - return itype_; - } + int itype() const { return itype_; } //! gets izeta_ - int izeta() const - { - return izeta_; - } + int izeta() const { return izeta_; } //! gets the angular momentum - int l() const - { - return l_; - } + int l() const { return l_; } //! gets the number of r-space grid points - int nr() const - { - return nr_; - } + int nr() const { return nr_; } //! gets the number of k-space grid points - int nk() const - { - return nk_; - } + int nk() const { return nk_; } //! gets r-space grid cutoff distance double rcut() const @@ -258,52 +240,28 @@ class NumericalRadial } //! gets the pointer to r-space grid points - const double* ptr_rgrid() const - { - return rgrid_; - } + const double* ptr_rgrid() const { return rgrid_; } //! gets the pointer to k-space grid points - const double* ptr_kgrid() const - { - return kgrid_; - } + const double* ptr_kgrid() const { return kgrid_; } //! gets the pointer to r-space values - const double* ptr_rvalue() const - { - return rvalue_; - } + const double* ptr_rvalue() const { return rvalue_; } //! gets the pointer to k-space values - const double* ptr_kvalue() const - { - return kvalue_; - } + const double* ptr_kvalue() const { return kvalue_; } //! gets the exponent of the pre-multiplied power term in rvalues_. @see pr_ - double pr() const - { - return pr_; - } + double pr() const { return pr_; } //! gets the exponent of the pre-multiplied power term in kvalues_. @see pk_ - double pk() const - { - return pk_; - } + double pk() const { return pk_; } //! gets the flag for FFT-compliancy. @see is_fft_compliant_ - bool is_fft_compliant() const - { - return is_fft_compliant_; - } + bool is_fft_compliant() const { return is_fft_compliant_; } //! gets the pointer to the SphericalBesselTransformer - const ModuleBase::SphericalBesselTransformer* ptr_sbt() const - { - return sbt_; - } + const ModuleBase::SphericalBesselTransformer* ptr_sbt() const { return sbt_; } ///@} diff --git a/source/module_basis/module_nao/radial_collection.cpp b/source/module_basis/module_nao/radial_collection.cpp index d237ce8bf6..15ed392123 100644 --- a/source/module_basis/module_nao/radial_collection.cpp +++ b/source/module_basis/module_nao/radial_collection.cpp @@ -36,15 +36,15 @@ void RadialCollection::cleanup() nchi_ = 0; } -void RadialCollection::build(const int nfile, const std::string* const file, const char type) +void RadialCollection::build(const int nfile, const std::string* const file, const char file_type) { - assert(type == 'o' || type == 'p'); + assert(file_type == 'o' || file_type == 'p'); cleanup(); ntype_ = nfile; radset_ = new RadialSet*[ntype_]; - switch (type) + switch (file_type) { case 'o': for (int itype = 0; itype < ntype_; ++itype) diff --git a/source/module_basis/module_nao/radial_collection.h b/source/module_basis/module_nao/radial_collection.h index 308a7e6277..08f59e11f8 100644 --- a/source/module_basis/module_nao/radial_collection.h +++ b/source/module_basis/module_nao/radial_collection.h @@ -87,6 +87,16 @@ class RadialCollection int ntype_ = 0; //!< number of RadialSet in the collection int lmax_ = -1; //!< maximum angular momentum of all NumericalRadial objects in the collection int nchi_ = 0; //!< total number of NumericalRadial objects in the collection + + //! array of RadialSet objects + /*! + * Each object hold a set of NumericalRadial objects that belong to a single element. + * The set could be either atomic orbitals or beta functions, in which case the + * underlying objects are AtomicRadials or BetaRadials, respectively. + * + * NOTE: AtomicRadials and BetaRadials do not necessarily have the same size as + * RadialSet. Therefore, a multilevel pointer is necessary for polymorphism. + * */ RadialSet** radset_ = nullptr; void cleanup(); diff --git a/source/module_basis/module_nao/radial_set.cpp b/source/module_basis/module_nao/radial_set.cpp index 7be0f1885a..113f8a3607 100644 --- a/source/module_basis/module_nao/radial_set.cpp +++ b/source/module_basis/module_nao/radial_set.cpp @@ -1,5 +1,7 @@ #include "module_basis/module_nao/radial_set.h" +#include + RadialSet::~RadialSet() { delete[] nzeta_; @@ -7,23 +9,28 @@ RadialSet::~RadialSet() delete[] index_map_; } -double RadialSet::rcut_max() const { +double RadialSet::rcut_max() const +{ double rmax = 0.0; - for (int i = 0; i != nchi_; ++i) { - if (chi_[i].rcut() > rmax) { + for (int i = 0; i != nchi_; ++i) + { + if (chi_[i].rcut() > rmax) + { rmax = chi_[i].rcut(); } } return rmax; } -RadialSet::RadialSet(const RadialSet& other) { +RadialSet::RadialSet(const RadialSet& other) +{ symbol_ = other.symbol_; itype_ = other.itype_; lmax_ = other.lmax_; nzeta_ = nullptr; - if (lmax_ >= 0) { + if (lmax_ >= 0) + { nzeta_ = new int[lmax_ + 1]; for (int l = 0; l <= lmax_; l++) nzeta_[l] = other.nzeta_[l]; @@ -32,24 +39,30 @@ RadialSet::RadialSet(const RadialSet& other) { nchi_ = other.nchi_; chi_ = nullptr; - if (nchi_ > 0) { + if (nchi_ > 0) + { chi_ = new NumericalRadial[nchi_]; - for (int i = 0; i < nchi_; i++) { + for (int i = 0; i < nchi_; i++) + { chi_[i] = other.chi_[i]; // deep copy } } index_map_ = nullptr; - if (lmax_ >= 0 && nzeta_max_ > 0) { + if (lmax_ >= 0 && nzeta_max_ > 0) + { index_map_ = new int[(lmax_ + 1) * nzeta_max_]; - for (int i = 0; i < (lmax_ + 1) * nzeta_max_; i++) { + for (int i = 0; i < (lmax_ + 1) * nzeta_max_; i++) + { index_map_[i] = other.index_map_[i]; } } } -RadialSet& RadialSet::operator=(const RadialSet& rhs) { - if (&rhs == this) { +RadialSet& RadialSet::operator=(const RadialSet& rhs) +{ + if (&rhs == this) + { return *this; } @@ -59,7 +72,8 @@ RadialSet& RadialSet::operator=(const RadialSet& rhs) { delete[] nzeta_; nzeta_ = nullptr; - if (lmax_ >= 0) { + if (lmax_ >= 0) + { nzeta_ = new int[lmax_ + 1]; for (int l = 0; l <= lmax_; l++) nzeta_[l] = rhs.nzeta_[l]; @@ -67,22 +81,24 @@ RadialSet& RadialSet::operator=(const RadialSet& rhs) { nzeta_max_ = rhs.nzeta_max_; nchi_ = rhs.nchi_; - //rcut_max_ = rhs.rcut_max_; - delete[] chi_; chi_ = nullptr; - if (nchi_ > 0) { + if (nchi_ > 0) + { chi_ = new NumericalRadial[nchi_]; - for (int i = 0; i < nchi_; i++) { + for (int i = 0; i < nchi_; i++) + { chi_[i] = rhs.chi_[i]; // deep copy } } delete[] index_map_; index_map_ = nullptr; - if (lmax_ >= 0 && nzeta_max_ > 0) { + if (lmax_ >= 0 && nzeta_max_ > 0) + { index_map_ = new int[(lmax_ + 1) * nzeta_max_]; - for (int i = 0; i < (lmax_ + 1) * nzeta_max_; i++) { + for (int i = 0; i < (lmax_ + 1) * nzeta_max_; i++) + { index_map_[i] = rhs.index_map_[i]; } } @@ -96,6 +112,36 @@ int RadialSet::index(int l, int izeta) const return index_map_[l * nzeta_max_ + izeta]; } +void RadialSet::indexing() +{ + if (!nzeta_) + { + return; + } + + assert(lmax_ >= 0); + nzeta_max_ = *std::max_element(nzeta_, nzeta_ + lmax_ + 1); + + delete[] index_map_; + index_map_ = new int[(lmax_ + 1) * nzeta_max_]; + int index_chi = 0; + for (int l = 0; l <= lmax_; ++l) + { + for (int izeta = 0; izeta != nzeta_max_; ++izeta) + { + if (izeta >= nzeta_[l]) + { + index_map_[l * nzeta_max_ + izeta] = -1; // -1 means no such orbital + } + else + { + index_map_[l * nzeta_max_ + izeta] = index_chi; + ++index_chi; + } + } + } +} + const NumericalRadial& RadialSet::chi(int l, int izeta) { int i = index_map_[l * nzeta_max_ + izeta]; @@ -113,7 +159,6 @@ void RadialSet::set_grid(const bool for_r_space, const int ngrid, const double* { for (int i = 0; i < nchi_; i++) chi_[i].set_grid(for_r_space, ngrid, grid, mode); - //rcut_max_ = grid[ngrid-1]; } void RadialSet::set_uniform_grid(const bool for_r_space, const int ngrid, const double cutoff, const char mode) @@ -133,8 +178,6 @@ void RadialSet::cleanup() nzeta_max_ = 0; nchi_ = 0; - //rcut_max_ = 0; - delete[] chi_; chi_ = nullptr; diff --git a/source/module_basis/module_nao/radial_set.h b/source/module_basis/module_nao/radial_set.h index ca4b8f2669..ae9aae8ec6 100644 --- a/source/module_basis/module_nao/radial_set.h +++ b/source/module_basis/module_nao/radial_set.h @@ -7,21 +7,20 @@ #include "module_basis/module_nao/numerical_radial.h" -//! An abstract class representing a related set of numerical radial functions +//! An abstract class representing a related set of numerical radial functions. /*! * This abstract class represents a set of numerical radial functions from - * a single file and of the same "itype" number. This is supposed to be the - * base class for concrete classes like AtomicRadials and BetaRadials, in - * which case "itype" is the element index in calculation, and they represent - * the set of all radial functions of the numerical atomic orbitals and - * Kleinman-Bylander beta functions of a single element respectively. + * a single file. This is supposed to be the base class for concrete classes + * like AtomicRadials and BetaRadials, in which case they represent the set + * of all radial functions of the numerical atomic orbitals and Kleinman- + * Bylander beta functions of a single element respectively. * * @see AtomicRadials BetaRadials * */ class RadialSet { public: - RadialSet(){}; + RadialSet() {} RadialSet(const RadialSet&); //!< deep copy RadialSet& operator=(const RadialSet&); //!< deep copy @@ -45,7 +44,11 @@ class RadialSet int lmax() const { return lmax_; } double rcut_max() const; - int nzeta(int l) const { return nzeta_[l]; } + int nzeta(int l) const + { + assert(l >= 0 && l <= lmax_); + return nzeta_[l]; + } int nzeta_max() const { return nzeta_max_; } int nchi() const { return nchi_; } @@ -88,6 +91,9 @@ class RadialSet //! get the index in chi_ array from (l,izeta) int index(const int l, const int izeta) const; + + //! calculate nzeta_max_ and build index_map_ from nzeta_ and lmax_ + void indexing(); }; #endif diff --git a/source/module_basis/module_nao/test/radial_collection_test.cpp b/source/module_basis/module_nao/test/radial_collection_test.cpp index f97fbf987e..c3cda85cee 100644 --- a/source/module_basis/module_nao/test/radial_collection_test.cpp +++ b/source/module_basis/module_nao/test/radial_collection_test.cpp @@ -12,7 +12,8 @@ * Tested functions: * * - build - * - parse an orbital file and initialize the NumericalRadial objects + * - parse a series of orbital or pseudopotential files and + * initialize AtomicRadials/BetaRadials objects accordingly. * * - all "getters" * - Get access to private members. @@ -27,18 +28,19 @@ class RadialCollectionTest : public ::testing::Test void TearDown(); RadialCollection orb; //!< object under test - int nfile = 0; // number of orbital files + int nfile = 0; // number of orbital/pseudopotential files std::string* file = nullptr; //!< orbitals file to read from std::string log_file = "./test_files/atomic_orbital.log"; //!< file for logging }; void RadialCollectionTest::SetUp() { + std::string dir = "../../../../../tests/PP_ORB/"; nfile = 4; file = new std::string[nfile]; - file[0] = "../../../../../tests/PP_ORB/C_gga_8au_100Ry_2s2p1d.orb"; - file[1] = "../../../../../tests/PP_ORB/H_gga_8au_60Ry_2s1p.orb"; - file[2] = "../../../../../tests/PP_ORB/O_gga_10au_100Ry_2s2p1d.orb"; - file[3] = "../../../../../tests/PP_ORB/Fe_gga_9au_100Ry_4s2p2d1f.orb"; + file[0] = dir + "C_gga_8au_100Ry_2s2p1d.orb"; + file[1] = dir + "H_gga_8au_60Ry_2s1p.orb"; + file[2] = dir + "O_gga_10au_100Ry_2s2p1d.orb"; + file[3] = dir + "Fe_gga_9au_100Ry_4s2p2d1f.orb"; } void RadialCollectionTest::TearDown() { @@ -84,11 +86,17 @@ TEST_F(RadialCollectionTest, BuildAndGet) { EXPECT_EQ(orb.nchi(3), 9); EXPECT_EQ(orb.nchi(), 22); - EXPECT_EQ(orb(0).itype(), 0); - EXPECT_EQ(orb(3).itype(), 3); + for (int itype = 0; itype <= 3; ++itype) { + EXPECT_EQ(orb(itype).itype(), itype); + } - EXPECT_EQ(orb(0, 1, 0).l(), 1); - EXPECT_EQ(orb(3, 3, 0).l(), 3); + for (int itype = 0; itype <= 3; ++itype) { + for (int l = 0; l <= orb(itype).lmax(); ++l) { + for (int izeta = 0; izeta != orb(itype).nzeta(l); ++izeta) { + EXPECT_EQ(orb(itype, l, izeta).l(), l); + } + } + } } TEST_F(RadialCollectionTest, BatchSet) { @@ -100,7 +108,7 @@ TEST_F(RadialCollectionTest, BatchSet) { EXPECT_EQ(orb.rcut_max(), 20.0); for (int itype = 0; itype != orb.ntype(); ++itype) { - for (int l = 0; l <= orb.lmax(); ++l) { + for (int l = 0; l <= orb(itype).lmax(); ++l) { for (int izeta = 0; izeta != orb.nzeta(itype, l); ++izeta) { EXPECT_EQ(&sbt, orb(itype, l, izeta).ptr_sbt()); EXPECT_DOUBLE_EQ(orb(itype, l, izeta).rcut(), 20.0); @@ -115,7 +123,7 @@ TEST_F(RadialCollectionTest, BatchSet) { orb.set_grid(true, 3, grid, 'i'); for (int itype = 0; itype != orb.ntype(); ++itype) { - for (int l = 0; l <= orb.lmax(); ++l) { + for (int l = 0; l <= orb(itype).lmax(); ++l) { for (int izeta = 0; izeta != orb.nzeta(itype, l); ++izeta) { EXPECT_DOUBLE_EQ(orb(itype, l, izeta).rcut(), 3.14); } From 02a74e72d893d24500c94766438cff3cc042875e Mon Sep 17 00:00:00 2001 From: jinzx10 Date: Wed, 7 Jun 2023 21:42:08 +0800 Subject: [PATCH 07/13] BetaRadials implemented --- source/module_basis/module_nao/CMakeLists.txt | 1 + .../module_basis/module_nao/beta_radials.cpp | 433 ++++++++++++++++++ source/module_basis/module_nao/beta_radials.h | 54 +++ .../module_nao/test/CMakeLists.txt | 10 + .../module_nao/test/beta_radials_test.cpp | 195 ++++++++ 5 files changed, 693 insertions(+) create mode 100644 source/module_basis/module_nao/beta_radials.cpp create mode 100644 source/module_basis/module_nao/beta_radials.h create mode 100644 source/module_basis/module_nao/test/beta_radials_test.cpp diff --git a/source/module_basis/module_nao/CMakeLists.txt b/source/module_basis/module_nao/CMakeLists.txt index 914adfeb69..2b848af517 100644 --- a/source/module_basis/module_nao/CMakeLists.txt +++ b/source/module_basis/module_nao/CMakeLists.txt @@ -5,6 +5,7 @@ if(ENABLE_LCAO) numerical_radial.cpp radial_set.cpp atomic_radials.cpp + beta_radials.cpp radial_collection.cpp ) diff --git a/source/module_basis/module_nao/beta_radials.cpp b/source/module_basis/module_nao/beta_radials.cpp new file mode 100644 index 0000000000..782d4e0e2f --- /dev/null +++ b/source/module_basis/module_nao/beta_radials.cpp @@ -0,0 +1,433 @@ +#include "module_basis/module_nao/beta_radials.h" + +#include "module_base/parallel_common.h" +#include "module_base/tool_quit.h" + +void BetaRadials::build(const std::string& file, const int itype, std::ofstream* ptr_log, const int rank) { + /* + * Build a BetaRadials object from beta functions read from a pseudopotential file + * + * NOTE: only the rank == 0 process reads the file + * */ + cleanup(); + + std::ifstream ifs; + bool is_open = false; + + if (rank == 0) + { + ifs.open(file); + is_open = ifs.is_open(); + } + +#ifdef __MPI + Parallel_Common::bcast_bool(is_open); +#endif + + if (!is_open) + { + ModuleBase::WARNING_QUIT("BetaRadials::read", "Couldn't open pseudopotential file: " + file); + } + + if (ptr_log) + { + (*ptr_log) << "\n\n\n\n"; + (*ptr_log) << " >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>" << std::endl; + (*ptr_log) << " | |" << std::endl; + (*ptr_log) << " | SETUP BETA PROJECTORS |" << std::endl; + (*ptr_log) << " | |" << std::endl; + (*ptr_log) << " | Projector information includes the cutoff radius, angular |" << std::endl; + (*ptr_log) << " | momentum, zeta number and numerical values on a radial grid. |" << std::endl; + (*ptr_log) << " | |" << std::endl; + (*ptr_log) << " <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<" << std::endl; + (*ptr_log) << "\n\n\n\n"; + } + + itype_ = itype; + + // identify the version of UPF file + // UPF 2.0.1 format starts with + int upf_version = 100; + if (rank == 0) { + std::string first_line; + std::getline(ifs, first_line); + if (first_line.find("2.0.1") != std::string::npos) { + upf_version = 201; + } + } + +#ifdef __MPI + Parallel_Common::bcast_int(upf_version); +#endif + + switch (upf_version) { + case 100: + read_beta_upf100(ifs, ptr_log, rank); + break; + case 201: + read_beta_upf201(ifs, ptr_log, rank); + break; + default: /* not supposed to happen */; + } + + if (rank == 0) + { + ifs.close(); + } +} + +void BetaRadials::read_beta_upf100(std::ifstream& ifs, std::ofstream* ptr_log, const int rank) { + /* + * Read the nonlocal beta functions from the pseudopotential file (in old UPF format) + * */ + std::stringstream ss; + std::string tmp; + std::string line; + int ngrid_max = 0; + bool is_good = false; + + if (rank == 0) { + /* + * Read the header, including + * + * 1. Element symbol + * 2. Maximum angular momentum + * 3. Number of radial grid points + * 4. Number of beta functions + * */ + while (std::getline(ifs, line)) { + if (line.find("Element") != std::string::npos) { + ss.str(""); + ss << line; + ss >> symbol_; + } else if (line.find("Max angular momentum component") != std::string::npos) { + ss.str(""); + ss << line; + ss >> lmax_; + } else if (line.find("Number of points in mesh") != std::string::npos) { + ss.str(""); + ss << line; + ss >> ngrid_max; + } else if (line.find("Number of Projectors") != std::string::npos) { + ss.str(""); + ss << line; + ss >> tmp >> nchi_; // nchi_ is the number of beta functions + } + + // should obtain valid nchi_, symbol_ & ngrid_max upon reaching the end of header + if (line.find("") != std::string::npos) { + // lmax could be -1 when there is no beta projectors, see, e.g., H.pz-vbc.UPF + is_good = (nchi_ >= 0) && symbol_.size() && (ngrid_max > 0); + break; + } + } + } + +#ifdef __MPI + Parallel_Common::bcast_bool(is_good); + Parallel_Common::bcast_string(symbol_); + Parallel_Common::bcast_int(lmax_); + Parallel_Common::bcast_int(nchi_); + Parallel_Common::bcast_int(ngrid_max); +#endif + + if (!is_good) + { + ModuleBase::WARNING_QUIT("BetaRadials::read_beta_upf100", "PP_HEADER error"); + } + + // In case some pseudopotential file does not have any beta function + if (nchi_ == 0) { + return; + } + + double* rgrid = new double[ngrid_max]; + if (rank == 0) { + /* + * Read the radial grid + * */ + while (ifs >> tmp) { + if (tmp == "") { + break; + } + } + assert(!ifs.eof()); + + for (int ir = 0; ir != ngrid_max; ++ir) { + ifs >> rgrid[ir]; + } + + ifs >> tmp; + assert(tmp == ""); + } + +#ifdef __MPI + Parallel_Common::bcast_double(rgrid, ngrid_max); +#endif + + assert(lmax_ >= 0); + nzeta_ = new int[lmax_ + 1]; + for (int l = 0; l <= lmax_; ++l) { + nzeta_[l] = 0; + } + + // rbeta, l & ngrid are directly read from file + double* rbeta = new double[ngrid_max]; // beta function is given as r*beta(r) + int l = 0; + int ngrid = 0; + + int l_last = -1; + int izeta = 0; + + chi_ = new NumericalRadial[nchi_]; + for (int i = 0; i != nchi_; ++i) { + if (rank == 0) { + /* + * Read the beta functions, including + * + * 1. Angular momentum + * 2. Number of actual radial grid points (should be smaller than ngrid_max) + * 3. Numerical values (r*beta(r)) on the radial grid + * + * and record the zeta number (izeta). + * */ + while (std::getline(ifs, line)) { + if (line.find("") != std::string::npos) { + break; + } + } + + ifs >> tmp >> l; + ifs >> tmp >> tmp; // skip "Beta" "L" + ifs >> ngrid; + + if (l == l_last) { + izeta += 1; + } else { + izeta = 0; + l_last = l; + } + + for (int ir = 0; ir != ngrid; ++ir) { + ifs >> rbeta[ir]; + } + + nzeta_[l] += 1; + + ifs >> tmp; + assert(tmp == ""); + } // rank == 0 + +#ifdef __MPI + Parallel_Common::bcast_int(l); + Parallel_Common::bcast_int(ngrid); + Parallel_Common::bcast_int(izeta); + Parallel_Common::bcast_double(rbeta, ngrid); +#endif + chi_[i].build(l, true, ngrid, rgrid, rbeta, 1, izeta, symbol_, itype_); + } + +#ifdef __MPI + Parallel_Common::bcast_int(nzeta_, lmax_ + 1); +#endif + + indexing(); + + delete[] rgrid; + delete[] rbeta; +} + +void BetaRadials::read_beta_upf201(std::ifstream& ifs, std::ofstream* ptr_log, const int rank) { + /* + * Read the nonlocal beta functions from the pseudopotential file (in UPF 2.0.1 format) + * */ + std::string tmp; + std::string line; + int ngrid_max = 0; + bool is_good = false; + + if (rank == 0) { + /* + * Read the header, including + * + * 1. Element symbol + * 2. Maximum angular momentum + * 3. Number of radial grid points + * 4. Number of beta functions + * */ + while (std::getline(ifs, line)) { + if (line.find("element=") != std::string::npos) { + symbol_ = trim201(line); + } else if (line.find("l_max=") != std::string::npos) { + lmax_ = std::stoi(trim201(line)); + } else if (line.find("mesh_size=") != std::string::npos) { + ngrid_max = std::stoi(trim201(line)); + } else if (line.find("number_of_proj=") != std::string::npos) { + nchi_ = std::stoi(trim201(line)); + } + + if (line.find("/>") != std::string::npos) { + is_good = (nchi_ >= 0) && (ngrid_max >= 0) && symbol_.size(); + break; + } + } + } + +#ifdef __MPI + Parallel_Common::bcast_bool(is_good); + Parallel_Common::bcast_string(symbol_); + Parallel_Common::bcast_int(lmax_); + Parallel_Common::bcast_int(nchi_); + Parallel_Common::bcast_int(ngrid_max); +#endif + + if (!is_good) + { + ModuleBase::WARNING_QUIT("BetaRadials::read_beta_upf201", "PP_HEADER error"); + } + + // In case some pseudopotential file does not have any beta function + if (nchi_ == 0) { + return; + } + + double* rgrid = new double[ngrid_max]; + if (rank == 0) { + /* + * Read the radial grid + * */ + while (std::getline(ifs, line)) { + if (line.find("> rgrid[ir]; + } + + ifs >> tmp; + assert(tmp == ""); + } + +#ifdef __MPI + Parallel_Common::bcast_double(rgrid, ngrid_max); +#endif + + nzeta_ = new int[lmax_ + 1]; + for (int l = 0; l <= lmax_; ++l) { + nzeta_[l] = 0; + } + + // rbeta, l & ngrid are directly read from file + double* rbeta = new double[ngrid_max]; + int l = 0; + int ngrid = 0; + + int l_last = -1; + int izeta = 0; + + chi_ = new NumericalRadial[nchi_]; + for (int i = 0; i != nchi_; ++i) { + if (rank == 0) { + /* + * Read the beta functions, including + * + * 1. Angular momentum + * 2. Number of actual radial grid points (should be smaller than ngrid_max) + * 3. Numerical values on the radial grid + * + * and record the zeta number (izeta). + * */ + while (std::getline(ifs, line)) { + if (line.find("") != std::string::npos) { + is_good &= (l >= 0) && (l <= lmax_) && (ngrid > 0) && (ngrid <= ngrid_max); + break; + } + } + + if (l == l_last) { + izeta += 1; + } else { + izeta = 0; + l_last = l; + } + + for (int ir = 0; ir != ngrid; ++ir) { + ifs >> rbeta[ir]; + } + + // reverse scan to determine the grid size + for (; ngrid > 0; --ngrid) { + if (std::abs(rbeta[ngrid-1]) > 1e-12) { + break; + } + } + + nzeta_[l] += 1; + + ifs >> tmp; + assert(tmp.find(" 0); + } // rank == 0 + +#ifdef __MPI + Parallel_Common::bcast_bool(is_good); + Parallel_Common::bcast_int(l); + Parallel_Common::bcast_int(ngrid); + Parallel_Common::bcast_int(izeta); + Parallel_Common::bcast_double(rbeta, ngrid); +#endif + + if (!is_good) + { + ModuleBase::WARNING_QUIT("BetaRadials::read_beta_upf201", "PP_BETA error"); + } + + chi_[i].build(l, true, ngrid, rgrid, rbeta, 1, izeta, symbol_, itype_); + } + +#ifdef __MPI + Parallel_Common::bcast_int(nzeta_, lmax_ + 1); +#endif + + indexing(); + + delete[] rgrid; + delete[] rbeta; +} + +std::string BetaRadials::trim201(std::string const& str) { + + // extract the substring between quotation marks (with whitespace trimmed) + // str MUST contain exactly a pair of quotation marks + + std::string::size_type start = str.find('"'); + std::string::size_type end = str.find_last_of('"'); + std::string tmp = str.substr(start+1, end-start-1); + + if (tmp.length() == 0) { + return tmp; + } + + start = tmp.find_first_not_of(" \t"); + end = tmp.find_last_not_of(" \t"); + return tmp.substr(start, end+1-start); +} + diff --git a/source/module_basis/module_nao/beta_radials.h b/source/module_basis/module_nao/beta_radials.h new file mode 100644 index 0000000000..4f11dd0bf9 --- /dev/null +++ b/source/module_basis/module_nao/beta_radials.h @@ -0,0 +1,54 @@ +#ifndef BETA_RADIALS_H_ +#define BETA_RADIALS_H_ + +#include "module_basis/module_nao/radial_set.h" + +//! The radial part of beta functions of a single element +/*! + * This class represents the radial part of all Kleinman-Bylander beta + * functions of a single element as read from a pseudopotential file. + * + * @see RadialSet + * + * Usage: + * + * int element_index = 1; + * std::ofstream ofs_log("/path/to/log/file"); + * std::string upf_file = "/path/to/pseudopotential/file"; + * + * BetaRadials O_beta; + * O_beta.build(orb_file, element_index, ofs_log, GlobalV::MY_RANK); + * + * */ +class BetaRadials: public RadialSet +{ + public: + BetaRadials() {} + ~BetaRadials() {} + + //! Build the class from a pseudopotential file + void build(const std::string& file, //!< pseudopotential file name + const int itype = 0, //!< element index in calculation + std::ofstream* ptr_log = nullptr, //!< output file stream for logging + const int rank = 0 //!< MPI rank + ); + + private: + + //! Read beta projectors from a pseudopotential file of UPF 1.0.0 format + void read_beta_upf100(std::ifstream& ifs, //!< input file stream from orbital file + std::ofstream* ptr_log = nullptr, //!< output file stream for logging + const int rank = 0 //!< MPI rank + ); + + //! Read beta projectors from a pseudopotential file of UPF 2.0.1 format + void read_beta_upf201(std::ifstream& ifs, //!< input file stream from orbital file + std::ofstream* ptr_log = nullptr, //!< output file stream for logging + const int rank = 0 //!< MPI rank + ); + + //! extract the substring between a pair of quotation marks (for UPF v2.0.1) + std::string trim201(std::string const& str); +}; + +#endif diff --git a/source/module_basis/module_nao/test/CMakeLists.txt b/source/module_basis/module_nao/test/CMakeLists.txt index 56d16e6d3d..a530d7fe69 100644 --- a/source/module_basis/module_nao/test/CMakeLists.txt +++ b/source/module_basis/module_nao/test/CMakeLists.txt @@ -16,6 +16,16 @@ AddTest( LIBS ${math_libs} device base ) +AddTest( + TARGET beta_radials + SOURCES + beta_radials_test.cpp + ../beta_radials.cpp + ../radial_set.cpp + ../numerical_radial.cpp + LIBS ${math_libs} device base +) + AddTest( TARGET radial_collection SOURCES diff --git a/source/module_basis/module_nao/test/beta_radials_test.cpp b/source/module_basis/module_nao/test/beta_radials_test.cpp new file mode 100644 index 0000000000..056dacf2c5 --- /dev/null +++ b/source/module_basis/module_nao/test/beta_radials_test.cpp @@ -0,0 +1,195 @@ +#include "module_basis/module_nao/beta_radials.h" + +#include "gtest/gtest.h" + +#ifdef __MPI +#include +#endif + +#include "module_base/constants.h" +#include "module_base/global_variable.h" + +/*********************************************************** + * Unit test of class "BetaRadials" + ***********************************************************/ +/*! + * Tested functions: + * + * - build + * - parse a pseudopotential file and initialize individual NumericalRadial objects + * + * - all "getters" + * - Get access to private members. + * + * - all "batch setters" + * - Set a property for all NumericalRadial objects at once + * */ +class BetaRadialsTest : public ::testing::Test +{ + protected: + void SetUp(); + void TearDown() {} + + BetaRadials beta; //!< object under test + + std::string dir = "../../../../../tests/PP_ORB/"; //!< directory with test files + std::string log_file = "./test_files/atomic_orbital.log"; //!< file for logging +}; + +void BetaRadialsTest::SetUp() { +#ifdef __MPI + MPI_Comm_rank(MPI_COMM_WORLD, &GlobalV::MY_RANK); +#endif +} + +TEST_F(BetaRadialsTest, ReadAndGet100) +{ + /* + * Read beta projectors from a UPF file of old format + * */ + std::string file100 = "Zn.LDA.UPF"; //!< a UPF 1.0.0 file to read from + beta.build(dir+file100, 314, nullptr, GlobalV::MY_RANK); + + EXPECT_EQ(beta.itype(), 314); + EXPECT_EQ(beta.symbol(), "Zn"); + EXPECT_EQ(beta.lmax(), 2); + EXPECT_EQ(beta.nzeta(0), 0); + EXPECT_EQ(beta.nzeta(1), 1); + EXPECT_EQ(beta.nzeta(2), 1); + EXPECT_EQ(beta.nzeta_max(), 1); + EXPECT_EQ(beta.nchi(), 2); + + EXPECT_DOUBLE_EQ(beta.rcut_max(), 3.36331676102); + + EXPECT_DOUBLE_EQ(beta.chi(1, 0).ptr_rgrid()[0] , 3.21829939971e-05); + EXPECT_DOUBLE_EQ(beta.chi(1, 0).ptr_rgrid()[4] , 3.39007851980e-05); + EXPECT_DOUBLE_EQ(beta.chi(1, 0).ptr_rgrid()[888], 3.31987661585e+00); + EXPECT_DOUBLE_EQ(beta.chi(1, 0).ptr_rgrid()[889], 3.36331676102e+00); + + EXPECT_DOUBLE_EQ(beta.chi(2, 0).ptr_rgrid()[0] , 3.21829939971e-05); + EXPECT_DOUBLE_EQ(beta.chi(2, 0).ptr_rgrid()[4] , 3.39007851980e-05); + EXPECT_DOUBLE_EQ(beta.chi(2, 0).ptr_rgrid()[887], 3.27699753773e+00); + EXPECT_DOUBLE_EQ(beta.chi(2, 0).ptr_rgrid()[888], 3.31987661585e+00); + + EXPECT_DOUBLE_EQ(beta.chi(1, 0).ptr_rvalue()[0] , -9.73759791529e-09); + EXPECT_DOUBLE_EQ(beta.chi(1, 0).ptr_rvalue()[4] , -1.08048430719e-08); + EXPECT_DOUBLE_EQ(beta.chi(1, 0).ptr_rvalue()[600], -5.66090228695e-02); + EXPECT_DOUBLE_EQ(beta.chi(1, 0).ptr_rvalue()[888], -1.89079314582e-11); + EXPECT_DOUBLE_EQ(beta.chi(1, 0).ptr_rvalue()[889], 1.99483148995e-12); + + EXPECT_DOUBLE_EQ(beta.chi(2, 0).ptr_rvalue()[0] , -2.84316719619e-11); + EXPECT_DOUBLE_EQ(beta.chi(2, 0).ptr_rvalue()[4] , -3.32316831459e-11); + EXPECT_DOUBLE_EQ(beta.chi(2, 0).ptr_rvalue()[600], -3.87224488590e-01); + EXPECT_DOUBLE_EQ(beta.chi(2, 0).ptr_rvalue()[887], 3.31894955664e-11); + EXPECT_DOUBLE_EQ(beta.chi(2, 0).ptr_rvalue()[888], -3.22751710952e-12); +} + +TEST_F(BetaRadialsTest, ReadAndGetVoid) { + /* + * This test reads a UPF file with no beta projectors. + * The test is to check that the code does not crash. + * */ + std::string file0 = "H.pz-vbc.UPF"; //!< a UPF file with no beta projectors + beta.build(dir+file0, 5, nullptr, GlobalV::MY_RANK); + + EXPECT_EQ(beta.nchi(), 0); + EXPECT_EQ(beta.itype(), 5); + EXPECT_EQ(beta.lmax(), -1); + EXPECT_EQ(beta.symbol(), "H"); +} + +TEST_F(BetaRadialsTest, ReadAndGet201) +{ + /* + * This test read beta projectors from a UPF file of 2.0.1 format + * */ + std::string file201 = "Pb_ONCV_PBE-1.0.upf"; //!< a UPF 2.0.1 file to read from + beta.build(dir+file201, 999, nullptr, GlobalV::MY_RANK); + + EXPECT_EQ(beta.itype(), 999); + EXPECT_EQ(beta.symbol(), "Pb"); + EXPECT_EQ(beta.lmax(), 3); + EXPECT_EQ(beta.nzeta(0), 2); + EXPECT_EQ(beta.nzeta(1), 2); + EXPECT_EQ(beta.nzeta(2), 2); + EXPECT_EQ(beta.nzeta(3), 2); + EXPECT_EQ(beta.nzeta_max(), 2); + EXPECT_EQ(beta.nchi(), 8); + + // NOTE: neither "cutoff_radius_index" nor "cutoff_radius" is reliable! + // the code reads all the values first and then reverse scan to determine the grid size + EXPECT_DOUBLE_EQ(beta.rcut_max(), 3.68); + + EXPECT_EQ(beta.chi(0,0).rcut(), 3.64); + EXPECT_EQ(beta.chi(0,0).nr(), 365); + EXPECT_EQ(beta.chi(0,0).izeta(), 0); + EXPECT_DOUBLE_EQ(beta.chi(0, 0).ptr_rgrid()[0] , 0.0000); + EXPECT_DOUBLE_EQ(beta.chi(0, 0).ptr_rgrid()[8] , 0.0800); + EXPECT_DOUBLE_EQ(beta.chi(0, 0).ptr_rgrid()[364], 3.6400); + EXPECT_DOUBLE_EQ(beta.chi(0, 0).ptr_rvalue()[0] , 0.0000000000e+00); + EXPECT_DOUBLE_EQ(beta.chi(0, 0).ptr_rvalue()[4] , 5.9689893417e-02); + EXPECT_DOUBLE_EQ(beta.chi(0, 0).ptr_rvalue()[364], -4.5888625103e-07); + + EXPECT_EQ(beta.chi(3,1).rcut(), 3.68); + EXPECT_EQ(beta.chi(3,1).nr(), 369); + EXPECT_EQ(beta.chi(3,1).izeta(), 1); + EXPECT_DOUBLE_EQ(beta.chi(3, 1).ptr_rgrid()[0] , 0.0000); + EXPECT_DOUBLE_EQ(beta.chi(3, 1).ptr_rgrid()[8] , 0.0800); + EXPECT_DOUBLE_EQ(beta.chi(3, 1).ptr_rgrid()[368], 3.6800); + EXPECT_DOUBLE_EQ(beta.chi(3, 1).ptr_rvalue()[0] , 0.0000000000e+00); + EXPECT_DOUBLE_EQ(beta.chi(3, 1).ptr_rvalue()[4] , 1.7908487484e-06); + EXPECT_DOUBLE_EQ(beta.chi(3, 1).ptr_rvalue()[368], -7.0309158570e-06); +} + +TEST_F(BetaRadialsTest, BatchSet) +{ + std::string file201 = "Pb_ONCV_PBE-1.0.upf"; //!< a UPF 2.0.1 file to read from + beta.build(dir+file201, 999, nullptr, GlobalV::MY_RANK); + + ModuleBase::SphericalBesselTransformer sbt; + beta.set_transformer(&sbt); + for (int l = 0; l != beta.lmax(); ++l) { + for (int izeta = 0; izeta != beta.nzeta(l); ++izeta) { + EXPECT_EQ(beta.chi(l, izeta).ptr_sbt(), &sbt); + } + } + + beta.set_uniform_grid(true, 2001, 20.0); + for (int l = 0; l != beta.lmax(); ++l) { + for (int izeta = 0; izeta != beta.nzeta(l); ++izeta) { + EXPECT_EQ(beta.chi(l, izeta).nr(), 2001); + EXPECT_EQ(beta.chi(l, izeta).rcut(), 20.0); + EXPECT_EQ(beta.chi(l, izeta).ptr_rgrid()[1500], 15.0); + EXPECT_EQ(beta.chi(l, izeta).ptr_rvalue()[1500], 0.0); + } + } + + double grid[5] = {0.0, 1.1, 2.2, 3.3, 4.4}; + beta.set_grid(true, 5, grid); + for (int l = 0; l != beta.lmax(); ++l) { + for (int izeta = 0; izeta != beta.nzeta(l); ++izeta) { + EXPECT_EQ(beta.chi(l, izeta).ptr_rgrid()[0], 0.0); + EXPECT_EQ(beta.chi(l, izeta).ptr_rgrid()[1], 1.1); + EXPECT_EQ(beta.chi(l, izeta).ptr_rgrid()[2], 2.2); + EXPECT_EQ(beta.chi(l, izeta).ptr_rgrid()[3], 3.3); + EXPECT_EQ(beta.chi(l, izeta).ptr_rgrid()[4], 4.4); + } + } +} + +int main(int argc, char** argv) +{ + +#ifdef __MPI + MPI_Init(&argc, &argv); +#endif + + testing::InitGoogleTest(&argc, argv); + int result = RUN_ALL_TESTS(); + +#ifdef __MPI + MPI_Finalize(); +#endif + + return result; +} From bdabee380a729600433bda358eb21349f26224e4 Mon Sep 17 00:00:00 2001 From: dyzheng Date: Thu, 8 Jun 2023 15:13:51 +0800 Subject: [PATCH 08/13] Fix: compiling warning by both Intel and GNU compilers (#2609) * Fix: warning detected by intel compiler * Fix: compiling warning detected by GNU compiler * Fix: error --------- Co-authored-by: dyzheng --- source/module_base/test/global_function_test.cpp | 6 ++++-- source/module_base/test/tool_threading_test.cpp | 2 ++ source/module_cell/test/unitcell_test_para.cpp | 6 ++++-- .../module_hamilt_general/module_xc/xc_functional_vxc.cpp | 5 +++++ .../hamilt_lcaodft/operator_lcao/operator_lcao.cpp | 5 +++++ source/module_io/restart.cpp | 6 ++++-- source/module_io/test/input_test.cpp | 2 +- 7 files changed, 25 insertions(+), 7 deletions(-) diff --git a/source/module_base/test/global_function_test.cpp b/source/module_base/test/global_function_test.cpp index f79bad0ce6..63fef746c2 100644 --- a/source/module_base/test/global_function_test.cpp +++ b/source/module_base/test/global_function_test.cpp @@ -417,8 +417,10 @@ TEST_F(GlobalFunctionTest, MakeDir) { GlobalV::MY_RANK = 0; ModuleBase::GlobalFunc::MAKE_DIR("scf"); - std::system("test -d "); - std::system("rm -r scf "); + auto error1 = std::system("test -d "); + EXPECT_EQ(error1, 0); + auto error2 = std::system("rm -r scf "); + EXPECT_EQ(error2, 0); SUCCEED(); } diff --git a/source/module_base/test/tool_threading_test.cpp b/source/module_base/test/tool_threading_test.cpp index ac2fc552c5..1624235b99 100644 --- a/source/module_base/test/tool_threading_test.cpp +++ b/source/module_base/test/tool_threading_test.cpp @@ -129,7 +129,9 @@ TEST(ToolThreadingTEST, OmpParalle) EXPECT_THAT(output1,testing::HasSubstr("1")); EXPECT_THAT(output1,testing::HasSubstr("0")); } +#ifndef _OPENMP #define _OPENMP +#endif TEST(ToolThreadingTEST, OmpParalleOpenmp) { std::string output2; diff --git a/source/module_cell/test/unitcell_test_para.cpp b/source/module_cell/test/unitcell_test_para.cpp index 20c90d26ac..11e1949a52 100644 --- a/source/module_cell/test/unitcell_test_para.cpp +++ b/source/module_cell/test/unitcell_test_para.cpp @@ -155,8 +155,10 @@ TEST_F(UcellTest,ReadPseudo) ifs.close(); std::string command1 = "test -d C && rm -rf C"; std::string command2 = "test -d H && rm -rf H"; - std::system( command1.c_str() ); - std::system( command2.c_str() ); + auto error1 = std::system( command1.c_str() ); + EXPECT_EQ(error1,0); + auto error2 = std::system( command2.c_str() ); + EXPECT_EQ(error2,0); } //read_cell_pseudopots //bcast_unitcell2 diff --git a/source/module_hamilt_general/module_xc/xc_functional_vxc.cpp b/source/module_hamilt_general/module_xc/xc_functional_vxc.cpp index 62096bfc03..44cb1751a7 100644 --- a/source/module_hamilt_general/module_xc/xc_functional_vxc.cpp +++ b/source/module_hamilt_general/module_xc/xc_functional_vxc.cpp @@ -454,6 +454,11 @@ std::tuple XC_Functional::v_xc_libxc( // Peiz ModuleBase::timer::tick("XC_Functional","v_xc_libxc"); return std::make_tuple( etxc, vtxc, std::move(v_nspin4) ); } // end if(4==GlobalV::NSPIN) + else//NSPIN != 1,2,4 is not supported + { + throw std::domain_error("GlobalV::NSPIN ="+std::to_string(GlobalV::NSPIN) + +" unfinished in "+std::string(__FILE__)+" line "+std::to_string(__LINE__)); + } } //the interface to libxc xc_mgga_exc_vxc(xc_func,n,rho,grho,laplrho,tau,e,v1,v2,v3,v4) diff --git a/source/module_hamilt_lcao/hamilt_lcaodft/operator_lcao/operator_lcao.cpp b/source/module_hamilt_lcao/hamilt_lcaodft/operator_lcao/operator_lcao.cpp index 63fccc47eb..3fc5d3a082 100644 --- a/source/module_hamilt_lcao/hamilt_lcaodft/operator_lcao/operator_lcao.cpp +++ b/source/module_hamilt_lcao/hamilt_lcaodft/operator_lcao/operator_lcao.cpp @@ -173,6 +173,11 @@ void OperatorLCAO::init(const int ik_in) break; } + default: + { + ModuleBase::WARNING_QUIT("OperatorLCAO::init", "unknown cal_type"); + break; + } } if(this->next_op != nullptr) {//it is not the last node, loop next init() function diff --git a/source/module_io/restart.cpp b/source/module_io/restart.cpp index 0b60f4d8a0..cd54e11ed2 100644 --- a/source/module_io/restart.cpp +++ b/source/module_io/restart.cpp @@ -24,7 +24,8 @@ void Restart::write_file2(const std::string &file_name, const void*const ptr, co { const int file = open(file_name.c_str(), O_WRONLY|O_CREAT|O_TRUNC, S_IRUSR|S_IWUSR); if(-1==file) throw std::runtime_error("can't open restart save file. \nerrno="+ModuleBase::GlobalFunc::TO_STRING(errno)+".\n"+ModuleBase::GlobalFunc::TO_STRING(__FILE__)+" line "+ModuleBase::GlobalFunc::TO_STRING(__LINE__)); - write(file, ptr, size); + auto error = write(file, ptr, size); + if(-1==error) throw std::runtime_error("can't write restart save file. \nerrno="+ModuleBase::GlobalFunc::TO_STRING(errno)+".\n"+ModuleBase::GlobalFunc::TO_STRING(__FILE__)+" line "+ModuleBase::GlobalFunc::TO_STRING(__LINE__) ); close(file); } @@ -32,7 +33,8 @@ void Restart::read_file2(const std::string &file_name, void*const ptr, const siz { const int file = open(file_name.c_str(), O_RDONLY); if(-1==file) throw std::runtime_error("can't open restart load file. \nerrno="+ModuleBase::GlobalFunc::TO_STRING(errno)+".\n"+ModuleBase::GlobalFunc::TO_STRING(__FILE__)+" line "+ModuleBase::GlobalFunc::TO_STRING(__LINE__)); - read(file, ptr, size); + auto error = read(file, ptr, size); + if(-1==error) throw std::runtime_error("can't read restart load file. \nerrno="+ModuleBase::GlobalFunc::TO_STRING(errno)+".\n"+ModuleBase::GlobalFunc::TO_STRING(__FILE__)+" line "+ModuleBase::GlobalFunc::TO_STRING(__LINE__)); close(file); } diff --git a/source/module_io/test/input_test.cpp b/source/module_io/test/input_test.cpp index 2bb481257c..24361d6f66 100644 --- a/source/module_io/test/input_test.cpp +++ b/source/module_io/test/input_test.cpp @@ -1456,7 +1456,7 @@ TEST_F(InputTest, Check) class ReadKSpacingTest : public ::testing::Test { protected: - void SetUp() + void SetUp() override { // create a temporary file for testing char tmpname[] = "tmpfile.tmp"; From c0cc66cbd12d62f63fc5c003c036c32fb0bbd75f Mon Sep 17 00:00:00 2001 From: dyzheng Date: Thu, 8 Jun 2023 15:14:31 +0800 Subject: [PATCH 09/13] Fix: optimized printing of start magnetization (#2605) * Fix: optimized printing of start magnetization * Fix: update refs in examples --------- Co-authored-by: dyzheng --- .../Pt-slab/running_scf.ref | 27 -------------- .../Pt-slab/running_scf.ref1 | 27 -------------- .../Pt-slab/running_scf.ref2 | 27 -------------- .../electric_field/Pt-slab/running_scf.ref | 27 -------------- examples/hse/lcao_Si2/running_scf.log_ref | 2 - source/module_base/global_function.h | 8 +--- .../module_base/test/global_function_test.cpp | 4 +- source/module_cell/read_atoms.cpp | 37 +++++++++++++++---- 8 files changed, 33 insertions(+), 126 deletions(-) diff --git a/examples/compensating_charge/Pt-slab/running_scf.ref b/examples/compensating_charge/Pt-slab/running_scf.ref index c3a70e0684..7534bfc886 100644 --- a/examples/compensating_charge/Pt-slab/running_scf.ref +++ b/examples/compensating_charge/Pt-slab/running_scf.ref @@ -55,33 +55,6 @@ L=1, number of zeta = 1 L=2, number of zeta = 1 number of atom for this type = 27 - start magnetization = FALSE - start magnetization = FALSE - start magnetization = FALSE - start magnetization = FALSE - start magnetization = FALSE - start magnetization = FALSE - start magnetization = FALSE - start magnetization = FALSE - start magnetization = FALSE - start magnetization = FALSE - start magnetization = FALSE - start magnetization = FALSE - start magnetization = FALSE - start magnetization = FALSE - start magnetization = FALSE - start magnetization = FALSE - start magnetization = FALSE - start magnetization = FALSE - start magnetization = FALSE - start magnetization = FALSE - start magnetization = FALSE - start magnetization = FALSE - start magnetization = FALSE - start magnetization = FALSE - start magnetization = FALSE - start magnetization = FALSE - start magnetization = FALSE TOTAL ATOM NUMBER = 27 diff --git a/examples/dipole_correction/Pt-slab/running_scf.ref1 b/examples/dipole_correction/Pt-slab/running_scf.ref1 index 6e739912b5..f53ed349bb 100644 --- a/examples/dipole_correction/Pt-slab/running_scf.ref1 +++ b/examples/dipole_correction/Pt-slab/running_scf.ref1 @@ -56,33 +56,6 @@ L=1, number of zeta = 1 L=2, number of zeta = 1 number of atom for this type = 27 - start magnetization = FALSE - start magnetization = FALSE - start magnetization = FALSE - start magnetization = FALSE - start magnetization = FALSE - start magnetization = FALSE - start magnetization = FALSE - start magnetization = FALSE - start magnetization = FALSE - start magnetization = FALSE - start magnetization = FALSE - start magnetization = FALSE - start magnetization = FALSE - start magnetization = FALSE - start magnetization = FALSE - start magnetization = FALSE - start magnetization = FALSE - start magnetization = FALSE - start magnetization = FALSE - start magnetization = FALSE - start magnetization = FALSE - start magnetization = FALSE - start magnetization = FALSE - start magnetization = FALSE - start magnetization = FALSE - start magnetization = FALSE - start magnetization = FALSE TOTAL ATOM NUMBER = 27 diff --git a/examples/dipole_correction/Pt-slab/running_scf.ref2 b/examples/dipole_correction/Pt-slab/running_scf.ref2 index f06d639ec1..63692a7f6b 100644 --- a/examples/dipole_correction/Pt-slab/running_scf.ref2 +++ b/examples/dipole_correction/Pt-slab/running_scf.ref2 @@ -56,33 +56,6 @@ L=1, number of zeta = 1 L=2, number of zeta = 1 number of atom for this type = 27 - start magnetization = FALSE - start magnetization = FALSE - start magnetization = FALSE - start magnetization = FALSE - start magnetization = FALSE - start magnetization = FALSE - start magnetization = FALSE - start magnetization = FALSE - start magnetization = FALSE - start magnetization = FALSE - start magnetization = FALSE - start magnetization = FALSE - start magnetization = FALSE - start magnetization = FALSE - start magnetization = FALSE - start magnetization = FALSE - start magnetization = FALSE - start magnetization = FALSE - start magnetization = FALSE - start magnetization = FALSE - start magnetization = FALSE - start magnetization = FALSE - start magnetization = FALSE - start magnetization = FALSE - start magnetization = FALSE - start magnetization = FALSE - start magnetization = FALSE TOTAL ATOM NUMBER = 27 diff --git a/examples/electric_field/Pt-slab/running_scf.ref b/examples/electric_field/Pt-slab/running_scf.ref index 9106c25e76..a06f820c69 100644 --- a/examples/electric_field/Pt-slab/running_scf.ref +++ b/examples/electric_field/Pt-slab/running_scf.ref @@ -56,33 +56,6 @@ L=1, number of zeta = 1 L=2, number of zeta = 1 number of atom for this type = 27 - start magnetization = FALSE - start magnetization = FALSE - start magnetization = FALSE - start magnetization = FALSE - start magnetization = FALSE - start magnetization = FALSE - start magnetization = FALSE - start magnetization = FALSE - start magnetization = FALSE - start magnetization = FALSE - start magnetization = FALSE - start magnetization = FALSE - start magnetization = FALSE - start magnetization = FALSE - start magnetization = FALSE - start magnetization = FALSE - start magnetization = FALSE - start magnetization = FALSE - start magnetization = FALSE - start magnetization = FALSE - start magnetization = FALSE - start magnetization = FALSE - start magnetization = FALSE - start magnetization = FALSE - start magnetization = FALSE - start magnetization = FALSE - start magnetization = FALSE TOTAL ATOM NUMBER = 27 diff --git a/examples/hse/lcao_Si2/running_scf.log_ref b/examples/hse/lcao_Si2/running_scf.log_ref index 6fe86a4053..507b7bb9a7 100644 --- a/examples/hse/lcao_Si2/running_scf.log_ref +++ b/examples/hse/lcao_Si2/running_scf.log_ref @@ -56,8 +56,6 @@ L=1, number of zeta = 2 L=2, number of zeta = 1 number of atom for this type = 2 - start magnetization = FALSE - start magnetization = FALSE TOTAL ATOM NUMBER = 2 diff --git a/source/module_base/global_function.h b/source/module_base/global_function.h index fa77c76b26..358454cc1d 100644 --- a/source/module_base/global_function.h +++ b/source/module_base/global_function.h @@ -49,9 +49,7 @@ void OUT(std::ofstream &ofs,const std::string &name,const T &a) template void OUT(std::ofstream &ofs,const std::string &name,const T &x, const T&y) { - std::stringstream name2; - name2 << "[" << name << "]"; - ofs<< " " << std::setw(40) < void OUT(std::ofstream &ofs,const std::string &name,const T &x, const T &y, const T &z) { - std::stringstream name2; - name2 << "[" << name << "]"; - ofs<< " " << std::setw(40) <0 to avoid too much output + //print when ia!=0 && mag[ia] != mag[0] to avoid too much output + if(ia==0 || (ia!=0 + && (atoms[it].m_loc_[ia].x != atoms[it].m_loc_[0].x + || atoms[it].m_loc_[ia].y != atoms[it].m_loc_[0].y + || atoms[it].m_loc_[ia].z != atoms[it].m_loc_[0].z))) + { + //use a stringstream to generate string: "concollinear magnetization of element it is:" + std::stringstream ss; + ss << "magnetization of element " << it+1; + if(ia!=0) + { + ss<<" (atom"<0 to avoid too much output + //print when ia!=0 && mag[ia] != mag[0] to avoid too much output + if(ia==0 || (ia!=0 && atoms[it].mag[ia] != atoms[it].mag[0])) + { + //use a stringstream to generate string: "cocollinear magnetization of element it is:" + std::stringstream ss; + ss << "magnetization of element " << it+1; + if(ia!=0) + { + ss<<" (atom"< Date: Thu, 8 Jun 2023 15:19:04 +0800 Subject: [PATCH 10/13] Test: UnitTest for module_tddft (#2580) * Merge develop branch to TDDFT branch (#2249) * Add the dipole and spectrum post-processing script. (#1928) * Add the Dipole and post-processing script necessary to calculate the spectrum using TDDFT. * fix td_val test bug. * Fix: bug in parsing strings mixed with " ", "\t" (#1951) * Use template to reconstruct parse_expression * Feature: output R matrix at each MD step * Modify'matrix_HS' to 'matrix' for R matrix output * Merge branches 'develop' and 'develop' of https://github.com/1041176461/abacus-develop into develop * Fix: modify index in parse_expression * Fix: bug in parsing strings mixed with " ", "\t" --------- Co-authored-by: jiyuang * solve conflicts * solve conflicts * solve conflicts * solve conflicts * solve conflicts * Fix : bug of multi-kpoint for tddft (#2023) * fix bug of electric field in tddft * fix bug of dipole and electric field * fix bug of multi-kpoint for tddft * unit conversion of out_efield * delete useless paramters for tddft (td_val_elec*) * delete useless parameters for tddft * pack all changes in tddft * solve conflicts * delete td_val --------- Co-authored-by: HeFuxiang94 <78629482+HeFuxiang94@users.noreply.github.com> Co-authored-by: jiyuyang <1041176461@qq.com> Co-authored-by: jiyuang * New propagator for better Convergence (#2253) * New propagator for better Convergence * Add files via upload * Add files via upload * Modified the code's implementation approach. I have modified the implementation approach of the code, and this one involves less changes to the original code. It also includes handling for multiple k-points. * Delete module_esolver directory * Delete module_hamilt_general directory * Delete module_hamilt_lcao directory * Add files via upload * Update dipole.py (#2270) * Simplify matrix multiplication and add new propagator method (#2272) * Add the dipole and spectrum post-processing script. (#1928) * Add the Dipole and post-processing script necessary to calculate the spectrum using TDDFT. * fix td_val test bug. * Fix: bug in parsing strings mixed with " ", "\t" (#1951) * Use template to reconstruct parse_expression * Feature: output R matrix at each MD step * Modify'matrix_HS' to 'matrix' for R matrix output * Merge branches 'develop' and 'develop' of https://github.com/1041176461/abacus-develop into develop * Fix: modify index in parse_expression * Fix: bug in parsing strings mixed with " ", "\t" --------- Co-authored-by: jiyuang * solve conflicts * solve conflicts * solve conflicts * solve conflicts * solve conflicts * Fix : bug of multi-kpoint for tddft (#2023) * fix bug of electric field in tddft * fix bug of dipole and electric field * fix bug of multi-kpoint for tddft * unit conversion of out_efield * delete useless paramters for tddft (td_val_elec*) * delete useless parameters for tddft * pack all changes in tddft * solve conflicts * delete td_val * simplify matrix multiplication and update new propagator for tddft * add autotests for new propagator * fix bug of ekb in tddft * delete temp output * update ref of new propagator * update CASES --------- Co-authored-by: HeFuxiang94 <78629482+HeFuxiang94@users.noreply.github.com> Co-authored-by: jiyuyang <1041176461@qq.com> Co-authored-by: jiyuang Co-authored-by: Qianrui <76200646+Qianruipku@users.noreply.github.com> * delete td_htype and update ref of tddft autotests * delete td_htype * refactor tddft code * move read_paramter of electric field to input_conv.cpp for tddft * sve conflicts * move definition of tag * convert read_parameter * add namespace module_tddft * remove tmp file * move read_paramters of electric field in tddft to input_conv.cpp * add force of tddft efield * add annotations * move read_parameters from H_TDDFT_pw.cpp to input_conv.cpp * add ifdef __LCAO * remove input parameter of read_td_efield * remove input parameter of read_td_efield * fix bug for input UTs * change name of parameters of compute_force in H_TDDFT_pw.h * improve annotation * solve conflicts * solve conflicts * add UT for read_td_efield * add ifdef __LCAO * add ifdef __LCAO * mistake for TEST and TEST_F * fix ReadTdEfieldTest * fix ReadTdEfieldTest * add doublenear() for EXPECT_EQ in ReadEfieldTest * add nampspace for doublenear * add nampspace for doublenear * mistake * replace doublenear with EXPECT_NEAR * TEST_F * delete nband in propagator * add UT for module tddft * comment add_test * fix bug of prepare() in H_TDDFT_pw.cpp * split tddft_test * fix bug * fix bug of UT for tddft * add tddft prefix * add UT for cblacs_gridinit * fix bug of UT for tddft * fix bug of UT for tddft --------- Co-authored-by: HeFuxiang94 <78629482+HeFuxiang94@users.noreply.github.com> Co-authored-by: jiyuyang <1041176461@qq.com> Co-authored-by: jiyuang Co-authored-by: HTZhao <104255052+ESROAMER@users.noreply.github.com> Co-authored-by: Omega <62006565+Satinelamp@users.noreply.github.com> Co-authored-by: Qianrui <76200646+Qianruipku@users.noreply.github.com> --- source/module_base/blacs_connector.h | 5 +- source/module_base/test/CMakeLists.txt | 5 + .../module_base/test/blacs_connector_test.cpp | 53 ++++++++ .../potentials/H_TDDFT_pw.cpp | 6 +- .../hamilt_lcaodft/FORCE_STRESS.cpp | 122 ++++++++--------- .../module_tddft/CMakeLists.txt | 4 + .../module_tddft/evolve_psi.cpp | 2 +- .../module_tddft/propagator.cpp | 22 ++- .../module_tddft/propagator.h | 16 +-- .../module_tddft/test/CMakeLists.txt | 34 +++++ .../module_tddft/test/bandenergy_test.cpp | 109 +++++++++++++++ .../module_tddft/test/middle_hamilt_test.cpp | 107 +++++++++++++++ .../module_tddft/test/norm_psi_test.cpp | 125 ++++++++++++++++++ .../module_tddft/test/propagator_test1.cpp | 114 ++++++++++++++++ .../module_tddft/test/propagator_test2.cpp | 119 +++++++++++++++++ .../module_tddft/test/propagator_test3.cpp | 124 +++++++++++++++++ .../module_tddft/test/tddft_test.cpp | 50 +++++++ .../module_tddft/test/tddft_test.h | 50 +++++++ .../module_tddft/test/upsi_test1.cpp | 76 +++++++++++ .../module_tddft/test/upsi_test2.cpp | 69 ++++++++++ .../module_tddft/test/upsi_test3.cpp | 86 ++++++++++++ 21 files changed, 1208 insertions(+), 90 deletions(-) create mode 100644 source/module_base/test/blacs_connector_test.cpp create mode 100644 source/module_hamilt_lcao/module_tddft/test/CMakeLists.txt create mode 100644 source/module_hamilt_lcao/module_tddft/test/bandenergy_test.cpp create mode 100644 source/module_hamilt_lcao/module_tddft/test/middle_hamilt_test.cpp create mode 100644 source/module_hamilt_lcao/module_tddft/test/norm_psi_test.cpp create mode 100644 source/module_hamilt_lcao/module_tddft/test/propagator_test1.cpp create mode 100644 source/module_hamilt_lcao/module_tddft/test/propagator_test2.cpp create mode 100644 source/module_hamilt_lcao/module_tddft/test/propagator_test3.cpp create mode 100644 source/module_hamilt_lcao/module_tddft/test/tddft_test.cpp create mode 100644 source/module_hamilt_lcao/module_tddft/test/tddft_test.h create mode 100644 source/module_hamilt_lcao/module_tddft/test/upsi_test1.cpp create mode 100644 source/module_hamilt_lcao/module_tddft/test/upsi_test2.cpp create mode 100644 source/module_hamilt_lcao/module_tddft/test/upsi_test3.cpp diff --git a/source/module_base/blacs_connector.h b/source/module_base/blacs_connector.h index d29ac47fae..98fd06fcde 100644 --- a/source/module_base/blacs_connector.h +++ b/source/module_base/blacs_connector.h @@ -33,7 +33,8 @@ extern "C" void Cblacs_gridmap(int* icontxt, int *usermap, int ldumap, int nprow, int npcol); // Informational and Miscellaneous void Cblacs_gridinfo(int icontxt, int* nprow, int *npcol, int *myprow, int *mypcol); - int Cblacs_pnum(int icontxt, int prow, int pcol); - void Cblacs_pcoord(int icontxt, int pnum, int *prow, int *pcol); + void Cblacs_gridinit(int* icontxt, char* layout, int nprow, int npcol); + int Cblacs_pnum(int icontxt, int prow, int pcol); + void Cblacs_pcoord(int icontxt, int pnum, int *prow, int *pcol); void Cblacs_exit(int icontxt); } diff --git a/source/module_base/test/CMakeLists.txt b/source/module_base/test/CMakeLists.txt index dc55ca3aed..c281e8beab 100644 --- a/source/module_base/test/CMakeLists.txt +++ b/source/module_base/test/CMakeLists.txt @@ -4,6 +4,11 @@ AddTest( LIBS ${math_libs} SOURCES blas_connector_test.cpp ) +AddTest( + TARGET base_blacs_connector + LIBS ${math_libs} + SOURCES blacs_connector_test.cpp +) AddTest( TARGET base_timer SOURCES timer_test.cpp ../timer.cpp diff --git a/source/module_base/test/blacs_connector_test.cpp b/source/module_base/test/blacs_connector_test.cpp new file mode 100644 index 0000000000..17cd92d268 --- /dev/null +++ b/source/module_base/test/blacs_connector_test.cpp @@ -0,0 +1,53 @@ +#include "../blacs_connector.h" + +#include + +#include "gtest/gtest.h" + +/************************************************ + * unit test of functions in blacs_connector.h + ***********************************************/ + +/** + * - Tested Function + * - Cblacs_gridinit + * - Initializes a grid of processors with a given number of rows and columns. + * The function creates a cartesian topology of all the processors initialized + * by the BLS library. In this topology, each processor is identified by its + * coordinates (row, col) in the grid. + */ + +TEST(blacs_connector, Cblacs_gridinit) +{ + int icontxt; + char layout[] = "ROW"; + int nprow = 1; + int npcol = 1; + + int myid, nprocs; + Cblacs_pinfo(&myid, &nprocs); + Cblacs_get(-1, 0, &icontxt); + + // Call the Cblacs_gridinit() function + Cblacs_gridinit(&icontxt, layout, nprow, npcol); + + // Check if the grid context handle is created successfully + EXPECT_EQ(icontxt, 0); +} + +int main(int argc, char** argv) +{ + int myrank; + int mysize; + + MPI_Init(&argc, &argv); + MPI_Comm_size(MPI_COMM_WORLD, &mysize); + MPI_Comm_rank(MPI_COMM_WORLD, &myrank); + + testing::InitGoogleTest(&argc, argv); + + int result = 0; + result = RUN_ALL_TESTS(); + MPI_Finalize(); + return 0; +} diff --git a/source/module_elecstate/potentials/H_TDDFT_pw.cpp b/source/module_elecstate/potentials/H_TDDFT_pw.cpp index c1b67cc64d..ef47e62b90 100644 --- a/source/module_elecstate/potentials/H_TDDFT_pw.cpp +++ b/source/module_elecstate/potentials/H_TDDFT_pw.cpp @@ -312,19 +312,19 @@ double H_TDDFT_pw::cal_v_time_heaviside() void H_TDDFT_pw::prepare(const UnitCell& cell, int& dir) { - if (dir == 0) + if (dir == 1) { bvec[0] = cell.G.e11; bvec[1] = cell.G.e12; bvec[2] = cell.G.e13; } - else if (dir == 1) + else if (dir == 2) { bvec[0] = cell.G.e21; bvec[1] = cell.G.e22; bvec[2] = cell.G.e23; } - else if (dir == 2) + else if (dir == 3) { bvec[0] = cell.G.e31; bvec[1] = cell.G.e32; diff --git a/source/module_hamilt_lcao/hamilt_lcaodft/FORCE_STRESS.cpp b/source/module_hamilt_lcao/hamilt_lcaodft/FORCE_STRESS.cpp index b72c9e26c1..29251e31da 100644 --- a/source/module_hamilt_lcao/hamilt_lcaodft/FORCE_STRESS.cpp +++ b/source/module_hamilt_lcao/hamilt_lcaodft/FORCE_STRESS.cpp @@ -16,30 +16,32 @@ double Force_Stress_LCAO::force_invalid_threshold_ev = 0.00; double Force_Stress_LCAO::output_acc = 1.0e-8; -Force_Stress_LCAO::Force_Stress_LCAO(Record_adj& ra, const int nat_in) : - RA(&ra), f_pw(nat_in), nat(nat_in){} -Force_Stress_LCAO::~Force_Stress_LCAO() {} - -void Force_Stress_LCAO::getForceStress( - const bool isforce, - const bool isstress, - const bool istestf, - const bool istests, - Local_Orbital_Charge& loc, - const elecstate::ElecState* pelec, - const psi::Psi* psid, - const psi::Psi>* psi, - LCAO_Hamilt &uhm, - ModuleBase::matrix& fcs, - ModuleBase::matrix & scs, - const Structure_Factor& sf, - const K_Vectors& kv, - ModulePW::PW_Basis* rhopw, +Force_Stress_LCAO::Force_Stress_LCAO(Record_adj& ra, const int nat_in) : RA(&ra), f_pw(nat_in), nat(nat_in) +{ +} +Force_Stress_LCAO::~Force_Stress_LCAO() +{ +} + +void Force_Stress_LCAO::getForceStress(const bool isforce, + const bool isstress, + const bool istestf, + const bool istests, + Local_Orbital_Charge& loc, + const elecstate::ElecState* pelec, + const psi::Psi* psid, + const psi::Psi>* psi, + LCAO_Hamilt& uhm, + ModuleBase::matrix& fcs, + ModuleBase::matrix& scs, + const Structure_Factor& sf, + const K_Vectors& kv, + ModulePW::PW_Basis* rhopw, #ifdef __EXX - Exx_LRI& exx_lri_double, - Exx_LRI>& exx_lri_complex, + Exx_LRI& exx_lri_double, + Exx_LRI>& exx_lri_complex, #endif - ModuleSymmetry::Symmetry* symm) + ModuleSymmetry::Symmetry* symm) { ModuleBase::TITLE("Force_Stress_LCAO", "getForceStress"); ModuleBase::timer::tick("Force_Stress_LCAO", "getForceStress"); @@ -154,13 +156,13 @@ void Force_Stress_LCAO::getForceStress( svl_dphi, #endif uhm, - kv); - //implement vdw force or stress here - // Peize Lin add 2014-04-04, update 2021-03-09 - // jiyy add 2019-05-18, update 2021-05-02 - ModuleBase::matrix force_vdw; - ModuleBase::matrix stress_vdw; - + kv); + // implement vdw force or stress here + // Peize Lin add 2014-04-04, update 2021-03-09 + // jiyy add 2019-05-18, update 2021-05-02 + ModuleBase::matrix force_vdw; + ModuleBase::matrix stress_vdw; + auto vdw_solver = vdw::make_vdw(GlobalC::ucell, INPUT); if (vdw_solver != nullptr) { @@ -219,13 +221,13 @@ void Force_Stress_LCAO::getForceStress( if (isforce) { force_dftu.create(nat, 3); - } - if(isstress) - { - stress_dftu.create(3, 3); - } - GlobalC::dftu.force_stress(loc.dm_gamma, loc.dm_k, *uhm.LM, force_dftu, stress_dftu, kv); - } + } + if (isstress) + { + stress_dftu.create(3, 3); + } + GlobalC::dftu.force_stress(loc.dm_gamma, loc.dm_k, *uhm.LM, force_dftu, stress_dftu, kv); + } if (!GlobalV::GAMMA_ONLY_LOCAL) this->flk.finish_k(); @@ -366,35 +368,37 @@ void Force_Stress_LCAO::getForceStress( GlobalC::ld.save_npy_f(fcs, "f_tot.npy", GlobalC::ucell.nat); // Ty/Bohr, F_tot if (GlobalV::deepks_scf) { - GlobalC::ld.save_npy_f(fcs - GlobalC::ld.F_delta, "f_base.npy", GlobalC::ucell.nat); //Ry/Bohr, F_base + GlobalC::ld.save_npy_f(fcs - GlobalC::ld.F_delta, "f_base.npy", GlobalC::ucell.nat); // Ry/Bohr, F_base - if(GlobalV::GAMMA_ONLY_LOCAL) + if (GlobalV::GAMMA_ONLY_LOCAL) { GlobalC::ld.cal_gdmx(loc.dm_gamma[0], - GlobalC::ucell, - GlobalC::ORB, - GlobalC::GridD, - pv->trace_loc_row, - pv->trace_loc_col, - isstress); + GlobalC::ucell, + GlobalC::ORB, + GlobalC::GridD, + pv->trace_loc_row, + pv->trace_loc_col, + isstress); } else - { - GlobalC::ld.cal_gdmx_k(loc.dm_k, - GlobalC::ucell, - GlobalC::ORB, - GlobalC::GridD, - pv->trace_loc_row, - pv->trace_loc_col, - kv.nks, - kv.kvec_d, - isstress); + { + GlobalC::ld.cal_gdmx_k(loc.dm_k, + GlobalC::ucell, + GlobalC::ORB, + GlobalC::GridD, + pv->trace_loc_row, + pv->trace_loc_col, + kv.nks, + kv.kvec_d, + isstress); } - if(GlobalV::deepks_out_unittest) GlobalC::ld.check_gdmx(GlobalC::ucell.nat); + if (GlobalV::deepks_out_unittest) + GlobalC::ld.check_gdmx(GlobalC::ucell.nat); GlobalC::ld.cal_gvx(GlobalC::ucell.nat); - if(GlobalV::deepks_out_unittest) GlobalC::ld.check_gvx(GlobalC::ucell.nat); - GlobalC::ld.save_npy_gvx(GlobalC::ucell.nat);// /Bohr, grad_vx + if (GlobalV::deepks_out_unittest) + GlobalC::ld.check_gvx(GlobalC::ucell.nat); + GlobalC::ld.save_npy_gvx(GlobalC::ucell.nat); // /Bohr, grad_vx } else { @@ -712,7 +716,7 @@ void Force_Stress_LCAO::calForceStressIntegralPart(const bool isGammaOnly, #else ModuleBase::matrix& svl_dphi, #endif - LCAO_Hamilt &uhm, + LCAO_Hamilt& uhm, const K_Vectors& kv) { if (isGammaOnly) @@ -759,7 +763,7 @@ void Force_Stress_LCAO::calForceStressIntegralPart(const bool isGammaOnly, svl_dphi, #endif uhm, - kv); + kv); } return; } diff --git a/source/module_hamilt_lcao/module_tddft/CMakeLists.txt b/source/module_hamilt_lcao/module_tddft/CMakeLists.txt index 8e0adf2b7d..f41dad565d 100644 --- a/source/module_hamilt_lcao/module_tddft/CMakeLists.txt +++ b/source/module_hamilt_lcao/module_tddft/CMakeLists.txt @@ -19,4 +19,8 @@ if(ENABLE_LCAO) add_coverage(tddft) endif() + IF (BUILD_TESTING) + add_subdirectory(test) + endif() + endif() diff --git a/source/module_hamilt_lcao/module_tddft/evolve_psi.cpp b/source/module_hamilt_lcao/module_tddft/evolve_psi.cpp index ebd578efdb..c30f9891c7 100644 --- a/source/module_hamilt_lcao/module_tddft/evolve_psi.cpp +++ b/source/module_hamilt_lcao/module_tddft/evolve_psi.cpp @@ -67,7 +67,7 @@ void evolve_psi(const int nband, /// @input Stmp, Htmp, print_matrix /// @output U_operator Propagator prop(propagator, pv); - prop.compute_propagator(nband, nlocal, Stmp, Htmp, H_laststep, U_operator, print_matrix); + prop.compute_propagator(nlocal, Stmp, Htmp, H_laststep, U_operator, print_matrix); // (3)->>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> diff --git a/source/module_hamilt_lcao/module_tddft/propagator.cpp b/source/module_hamilt_lcao/module_tddft/propagator.cpp index f2711ed2f8..c48e4142ff 100644 --- a/source/module_hamilt_lcao/module_tddft/propagator.cpp +++ b/source/module_hamilt_lcao/module_tddft/propagator.cpp @@ -23,8 +23,7 @@ inline int globalIndex(int localindex, int nblk, int nprocs, int myproc) return gIndex; } -void Propagator::compute_propagator(const int nband, - const int nlocal, +void Propagator::compute_propagator(const int nlocal, const std::complex* Stmp, const std::complex* Htmp, const std::complex* H_laststep, @@ -35,16 +34,16 @@ void Propagator::compute_propagator(const int nband, switch (ptype) { case 0: - compute_propagator_cn2(nband, nlocal, Stmp, Htmp, U_operator, print_matrix); + compute_propagator_cn2(nlocal, Stmp, Htmp, U_operator, print_matrix); break; case 1: tag = 1; - compute_propagator_taylor(nband, nlocal, Stmp, Htmp, U_operator, print_matrix, tag); + compute_propagator_taylor(nlocal, Stmp, Htmp, U_operator, print_matrix, tag); break; case 2: - compute_propagator_etrs(nband, nlocal, Stmp, Htmp, H_laststep, U_operator, print_matrix); + compute_propagator_etrs(nlocal, Stmp, Htmp, H_laststep, U_operator, print_matrix); break; @@ -54,8 +53,7 @@ void Propagator::compute_propagator(const int nband, } } -void Propagator::compute_propagator_cn2(const int nband, - const int nlocal, +void Propagator::compute_propagator_cn2(const int nlocal, const std::complex* Stmp, const std::complex* Htmp, std::complex* U_operator, @@ -259,8 +257,7 @@ void Propagator::compute_propagator_cn2(const int nband, delete[] ipiv; } -void Propagator::compute_propagator_taylor(const int nband, - const int nlocal, +void Propagator::compute_propagator_taylor(const int nlocal, const std::complex* Stmp, const std::complex* Htmp, std::complex* U_operator, @@ -580,8 +577,7 @@ void Propagator::compute_propagator_taylor(const int nband, delete[] ipiv; } -void Propagator::compute_propagator_etrs(const int nband, - const int nlocal, +void Propagator::compute_propagator_etrs(const int nlocal, const std::complex* Stmp, const std::complex* Htmp, const std::complex* H_laststep, @@ -593,8 +589,8 @@ void Propagator::compute_propagator_etrs(const int nband, ModuleBase::GlobalFunc::ZEROS(U1, this->ParaV->nloc); ModuleBase::GlobalFunc::ZEROS(U2, this->ParaV->nloc); int tag = 2; - compute_propagator_taylor(nband, nlocal, Stmp, Htmp, U1, print_matrix, tag); - compute_propagator_taylor(nband, nlocal, Stmp, H_laststep, U2, print_matrix, tag); + compute_propagator_taylor(nlocal, Stmp, Htmp, U1, print_matrix, tag); + compute_propagator_taylor(nlocal, Stmp, H_laststep, U2, print_matrix, tag); ScalapackConnector::gemm('N', 'N', nlocal, diff --git a/source/module_hamilt_lcao/module_tddft/propagator.h b/source/module_hamilt_lcao/module_tddft/propagator.h index 7d9e65b2c4..a156ca1355 100644 --- a/source/module_hamilt_lcao/module_tddft/propagator.h +++ b/source/module_hamilt_lcao/module_tddft/propagator.h @@ -24,7 +24,6 @@ class Propagator /** * @brief compute propagator * - * @param[in] nband number of bands * @param[in] nlocal number of orbitals * @param[in] Stmp overlap matrix * @param[in] Htmp H(t+dt/2) or H(t+dt) @@ -32,8 +31,7 @@ class Propagator * @param[in] print_matirx print internal matrix or not * @param[out] U_operator operator of propagator */ - void compute_propagator(const int nband, - const int nlocal, + void compute_propagator(const int nlocal, const std::complex* Stmp, const std::complex* Htmp, const std::complex* H_laststep, @@ -50,15 +48,13 @@ class Propagator /** * @brief compute propagator of method Crank-Nicolson * - * @param[in] nband number of bands * @param[in] nlocal number of orbitals * @param[in] Stmp overlap matrix * @param[in] Htmp H(t+dt/2) or H(t+dt) * @param[in] print_matirx print internal matrix or not * @param[out] U_operator operator of propagator */ - void compute_propagator_cn2(const int nband, - const int nlocal, + void compute_propagator_cn2(const int nlocal, const std::complex* Stmp, const std::complex* Htmp, std::complex* U_operator, @@ -67,7 +63,6 @@ class Propagator /** * @brief compute propagator of method 4th Taylor * - * @param[in] nband number of bands * @param[in] nlocal number of orbitals * @param[in] Stmp overlap matrix * @param[in] Htmp H(t+dt/2) or H(t+dt) @@ -75,8 +70,7 @@ class Propagator * @param[in] tag a parametre different for 4th Taylor and ETRS * @param[out] U_operator operator of propagator */ - void compute_propagator_taylor(const int nband, - const int nlocal, + void compute_propagator_taylor(const int nlocal, const std::complex* Stmp, const std::complex* Htmp, std::complex* U_operator, @@ -86,7 +80,6 @@ class Propagator /** * @brief compute propagator of method ETRS * - * @param[in] nband number of bands * @param[in] nlocal number of orbitals * @param[in] Stmp overlap matrix * @param[in] Htmp H(t+dt/2) or H(t+dt) @@ -94,8 +87,7 @@ class Propagator * @param[in] print_matirx print internal matrix or not * @param[out] U_operator operator of propagator */ - void compute_propagator_etrs(const int nband, - const int nlocal, + void compute_propagator_etrs(const int nlocal, const std::complex* Stmp, const std::complex* Htmp, const std::complex* H_laststep, diff --git a/source/module_hamilt_lcao/module_tddft/test/CMakeLists.txt b/source/module_hamilt_lcao/module_tddft/test/CMakeLists.txt new file mode 100644 index 0000000000..675806ca8d --- /dev/null +++ b/source/module_hamilt_lcao/module_tddft/test/CMakeLists.txt @@ -0,0 +1,34 @@ +remove_definitions(-D __MPI) + +add_library(tddft_test_lib tddft_test.cpp) + +AddTest( + TARGET tddft_middle_hamilt_test + LIBS ${math_libs} base device tddft_test_lib + SOURCES middle_hamilt_test.cpp ../middle_hamilt.cpp +) + +AddTest( + TARGET tddft_bandenergy_test + LIBS ${math_libs} base device tddft_test_lib + SOURCES bandenergy_test.cpp ../bandenergy.cpp +) + +AddTest( + TARGET tddft_norm_psi_test + LIBS ${math_libs} base device tddft_test_lib + SOURCES norm_psi_test.cpp ../norm_psi.cpp +) + +AddTest( + TARGET tddft_upsi_test + LIBS ${math_libs} base device tddft_test_lib + SOURCES upsi_test1.cpp upsi_test2.cpp upsi_test3.cpp ../upsi.cpp +) + +AddTest( + TARGET tddft_propagator_test + LIBS ${math_libs} base device tddft_test_lib + SOURCES propagator_test1.cpp propagator_test2.cpp propagator_test3.cpp ../propagator.cpp +) + diff --git a/source/module_hamilt_lcao/module_tddft/test/bandenergy_test.cpp b/source/module_hamilt_lcao/module_tddft/test/bandenergy_test.cpp new file mode 100644 index 0000000000..45ef2d6359 --- /dev/null +++ b/source/module_hamilt_lcao/module_tddft/test/bandenergy_test.cpp @@ -0,0 +1,109 @@ +#include "module_hamilt_lcao/module_tddft/bandenergy.h" + +#include +#include +#include + +#include "module_basis/module_ao/parallel_orbitals.h" +#include "module_hamilt_lcao/module_tddft/evolve_elec.h" +#include "tddft_test.h" + +/************************************************ + * unit test of functions in bandenergy.h + ***********************************************/ + +/** + * - Tested Function + * - compute_ekb + * - compute band energy ekb . + */ + +#define doublethreshold 1e-8 +double module_tddft::Evolve_elec::td_print_eij = -1; + +Parallel_Orbitals::Parallel_Orbitals() +{ +} +Parallel_Orbitals::~Parallel_Orbitals() +{ +} + +TEST(BandEnergyTest, testBandEnergy) +{ + std::complex* psi_k; + std::complex* Htmp; + double* ekb; + int nband = 3; + int nlocal = 4; + bool print_matrix = false; + Parallel_Orbitals* pv; + pv = new Parallel_Orbitals(); + pv->nloc = nlocal * nlocal; + pv->nloc_wfc = nlocal * nband; + pv->ncol = nlocal; + pv->nrow = nlocal; + pv->ncol_bands = nband; + pv->dim0 = 1; + pv->dim1 = 1; + pv->nb = 1; + + int dim[2]; + int period[2] = {1, 1}; + int reorder = 0; + dim[0] = nprow; + dim[1] = npcol; + + MPI_Cart_create(MPI_COMM_WORLD, 2, dim, period, reorder, &pv->comm_2D); + + // Initialize input matrices + int info; + int mb = 1, nb = 1, lda = nband, ldc = nlocal; + int irsrc = 0, icsrc = 0, lld = numroc_(&nlocal, &mb, &myprow, &irsrc, &nprow), + lld1 = numroc_(&nband, &mb, &myprow, &irsrc, &nprow); + descinit_(pv->desc, &nlocal, &nlocal, &mb, &nb, &irsrc, &icsrc, &ictxt, &lld, &info); + descinit_(pv->desc_wfc, &nlocal, &nband, &mb, &nb, &irsrc, &icsrc, &ictxt, &lld, &info); + descinit_(pv->desc_Eij, &nband, &nband, &mb, &nb, &irsrc, &icsrc, &ictxt, &lld, &info); + + // Initialize data + Htmp = new std::complex[nlocal * nlocal]; + psi_k = new std::complex[nlocal * nband]; + ekb = new double[nband]; + + for (int i = 0; i < nlocal; ++i) + { + for (int j = 0; j < nlocal; ++j) + { + if (i == j) + { + Htmp[i * nlocal + j] = std::complex(1.0, 0.0); + } + } + } + Htmp[1] = 0.5; + Htmp[4] = 0.5; + + psi_k[0] = 1.0; + psi_k[1] = 1.0; + psi_k[2] = 0.0; + psi_k[3] = 0.0; + psi_k[4] = 2.0; + psi_k[5] = 1.0; + psi_k[6] = 1.0; + psi_k[7] = 0.0; + psi_k[8] = 3.0; + psi_k[9] = 0.0; + psi_k[10] = 0.0; + psi_k[11] = 1.0; + + // Call the function + module_tddft::compute_ekb(pv, nband, nlocal, Htmp, psi_k, ekb); + + // Check the results + EXPECT_NEAR(ekb[0], 3.0, doublethreshold); + EXPECT_NEAR(ekb[1], 8.0, doublethreshold); + EXPECT_NEAR(ekb[2], 10.0, doublethreshold); + + delete[] psi_k; + delete[] Htmp; + delete[] ekb; +} diff --git a/source/module_hamilt_lcao/module_tddft/test/middle_hamilt_test.cpp b/source/module_hamilt_lcao/module_tddft/test/middle_hamilt_test.cpp new file mode 100644 index 0000000000..8387a82de7 --- /dev/null +++ b/source/module_hamilt_lcao/module_tddft/test/middle_hamilt_test.cpp @@ -0,0 +1,107 @@ +#include "module_hamilt_lcao/module_tddft/middle_hamilt.h" + +#include +#include +#include + +#include "module_basis/module_ao/parallel_orbitals.h" +#include "tddft_test.h" + +/************************************************ + * unit test of functions in middle_hamilt.h + ***********************************************/ + +/** + * - Tested Function + * - half_Hmatrix + * - compute H(t+dt/2). + */ + +#define doublethreshold 1e-8 + +Parallel_Orbitals::Parallel_Orbitals() +{ +} +Parallel_Orbitals::~Parallel_Orbitals() +{ +} + +TEST(MiddleHamiltTest, testMiddleHamilt) +{ + std::complex* Htmp; + std::complex* Hlaststep; + int nband = 3; + int nlocal = 4; + bool print_matrix = false; + Parallel_Orbitals* pv; + pv = new Parallel_Orbitals(); + pv->nloc = nlocal * nlocal; + pv->ncol = nlocal; + pv->nrow = nlocal; + + // Initialize input matrices + int info; + int mb = 1, nb = 1, lda = nlocal, ldc = nlocal; + int irsrc = 0, icsrc = 0, lld = numroc_(&nlocal, &mb, &myprow, &irsrc, &nprow); + descinit_(pv->desc, &nlocal, &nlocal, &mb, &nb, &irsrc, &icsrc, &ictxt, &lld, &info); + + // Initialize data + Htmp = new std::complex[nlocal * nlocal]; + Hlaststep = new std::complex[nlocal * nlocal]; + + for (int i = 0; i < nlocal; ++i) + { + for (int j = 0; j < nlocal; ++j) + { + if (i == j) + { + Htmp[i * nlocal + j] = std::complex(1.0, 0.0); + Hlaststep[i * nlocal + j] = std::complex(1.0 + 0.2 * (i * nlocal + j), 0.0); + } + else + { + Hlaststep[i * nlocal + j] = std::complex(0.2 * (i * nlocal + j), 0.0); + } + } + } + + // Call the function + module_tddft::half_Hmatrix(pv, nband, nlocal, Htmp, Hlaststep, print_matrix); + + // Check the results + EXPECT_NEAR(Htmp[0].real(), 1.0, doublethreshold); + EXPECT_NEAR(Htmp[0].imag(), 0.0, doublethreshold); + EXPECT_NEAR(Htmp[1].real(), 0.1, doublethreshold); + EXPECT_NEAR(Htmp[1].imag(), 0.0, doublethreshold); + EXPECT_NEAR(Htmp[2].real(), 0.2, doublethreshold); + EXPECT_NEAR(Htmp[2].imag(), 0.0, doublethreshold); + EXPECT_NEAR(Htmp[3].real(), 0.3, doublethreshold); + EXPECT_NEAR(Htmp[3].imag(), 0.0, doublethreshold); + EXPECT_NEAR(Htmp[4].real(), 0.4, doublethreshold); + EXPECT_NEAR(Htmp[4].imag(), 0.0, doublethreshold); + EXPECT_NEAR(Htmp[5].real(), 1.5, doublethreshold); + EXPECT_NEAR(Htmp[5].imag(), 0.0, doublethreshold); + EXPECT_NEAR(Htmp[6].real(), 0.6, doublethreshold); + EXPECT_NEAR(Htmp[6].imag(), 0.0, doublethreshold); + EXPECT_NEAR(Htmp[7].real(), 0.7, doublethreshold); + EXPECT_NEAR(Htmp[7].imag(), 0.0, doublethreshold); + EXPECT_NEAR(Htmp[8].real(), 0.8, doublethreshold); + EXPECT_NEAR(Htmp[8].imag(), 0.0, doublethreshold); + EXPECT_NEAR(Htmp[9].real(), 0.9, doublethreshold); + EXPECT_NEAR(Htmp[9].imag(), 0.0, doublethreshold); + EXPECT_NEAR(Htmp[10].real(), 2.0, doublethreshold); + EXPECT_NEAR(Htmp[10].imag(), 0.0, doublethreshold); + EXPECT_NEAR(Htmp[11].real(), 1.1, doublethreshold); + EXPECT_NEAR(Htmp[11].imag(), 0.0, doublethreshold); + EXPECT_NEAR(Htmp[12].real(), 1.2, doublethreshold); + EXPECT_NEAR(Htmp[12].imag(), 0.0, doublethreshold); + EXPECT_NEAR(Htmp[13].real(), 1.3, doublethreshold); + EXPECT_NEAR(Htmp[13].imag(), 0.0, doublethreshold); + EXPECT_NEAR(Htmp[14].real(), 1.4, doublethreshold); + EXPECT_NEAR(Htmp[14].imag(), 0.0, doublethreshold); + EXPECT_NEAR(Htmp[15].real(), 2.5, doublethreshold); + EXPECT_NEAR(Htmp[15].imag(), 0.0, doublethreshold); + + delete[] Htmp; + delete[] Hlaststep; +} diff --git a/source/module_hamilt_lcao/module_tddft/test/norm_psi_test.cpp b/source/module_hamilt_lcao/module_tddft/test/norm_psi_test.cpp new file mode 100644 index 0000000000..5a1bdc0361 --- /dev/null +++ b/source/module_hamilt_lcao/module_tddft/test/norm_psi_test.cpp @@ -0,0 +1,125 @@ +#include "module_hamilt_lcao/module_tddft/norm_psi.h" + +#include +#include +#include + +#include "module_basis/module_ao/parallel_orbitals.h" +#include "tddft_test.h" + +/************************************************ + * unit test of functions in norm_psi.h + ***********************************************/ + +/** + * - Tested Function + * - norm_psi + * - normalize the wave function. + */ + +#define doublethreshold 1e-8 + +Parallel_Orbitals::Parallel_Orbitals() +{ +} +Parallel_Orbitals::~Parallel_Orbitals() +{ +} + +TEST(NormPsiTest, testNormPsi) +{ + std::complex* psi_k; + std::complex* Stmp; + int nband = 3; + int nlocal = 4; + bool print_matrix = false; + Parallel_Orbitals* pv; + pv = new Parallel_Orbitals(); + pv->nloc = nlocal * nlocal; + pv->nloc_wfc = nlocal * nband; + pv->ncol = nlocal; + pv->nrow = nlocal; + pv->ncol_bands = nband; + pv->dim0 = 1; + pv->dim1 = 1; + pv->nb = 1; + + int dim[2]; + int period[2] = {1, 1}; + int reorder = 0; + dim[0] = nprow; + dim[1] = npcol; + + MPI_Cart_create(MPI_COMM_WORLD, 2, dim, period, reorder, &pv->comm_2D); + + // Initialize input matrices + int info; + int mb = 1, nb = 1, lda = nband, ldc = nlocal; + int irsrc = 0, icsrc = 0, lld = numroc_(&nlocal, &mb, &myprow, &irsrc, &nprow), + lld1 = numroc_(&nband, &mb, &myprow, &irsrc, &nprow); + descinit_(pv->desc, &nlocal, &nlocal, &mb, &nb, &irsrc, &icsrc, &ictxt, &lld, &info); + descinit_(pv->desc_wfc, &nlocal, &nband, &mb, &nb, &irsrc, &icsrc, &ictxt, &lld, &info); + descinit_(pv->desc_Eij, &nband, &nband, &mb, &nb, &irsrc, &icsrc, &ictxt, &lld, &info); + + // Initialize data + Stmp = new std::complex[nlocal * nlocal]; + psi_k = new std::complex[nlocal * nband]; + + for (int i = 0; i < nlocal; ++i) + { + for (int j = 0; j < nlocal; ++j) + { + if (i == j) + { + Stmp[i * nlocal + j] = std::complex(1.0, 0.0); + } + } + } + Stmp[1] = 0.5; + Stmp[4] = 0.5; + + psi_k[0] = 1.0; + psi_k[1] = 1.0; + psi_k[2] = 0.0; + psi_k[3] = 0.0; + psi_k[4] = 2.0; + psi_k[5] = 1.0; + psi_k[6] = 1.0; + psi_k[7] = 0.0; + psi_k[8] = 3.0; + psi_k[9] = 0.0; + psi_k[10] = 0.0; + psi_k[11] = 1.0; + + // Call the function + module_tddft::norm_psi(pv, nband, nlocal, Stmp, psi_k, print_matrix); + + // Check the results + EXPECT_NEAR(psi_k[0].real(), 0.577350269189626, doublethreshold); + EXPECT_NEAR(psi_k[0].imag(), 0.0, doublethreshold); + EXPECT_NEAR(psi_k[1].real(), 0.577350269189626, doublethreshold); + EXPECT_NEAR(psi_k[1].imag(), 0.0, doublethreshold); + EXPECT_NEAR(psi_k[2].real(), 0.0, doublethreshold); + EXPECT_NEAR(psi_k[2].imag(), 0.0, doublethreshold); + EXPECT_NEAR(psi_k[3].real(), 0.0, doublethreshold); + EXPECT_NEAR(psi_k[3].imag(), 0.0, doublethreshold); + EXPECT_NEAR(psi_k[4].real(), 0.707106781186547, doublethreshold); + EXPECT_NEAR(psi_k[4].imag(), 0.0, doublethreshold); + EXPECT_NEAR(psi_k[5].real(), 0.353553390593274, doublethreshold); + EXPECT_NEAR(psi_k[5].imag(), 0.0, doublethreshold); + EXPECT_NEAR(psi_k[6].real(), 0.353553390593274, doublethreshold); + EXPECT_NEAR(psi_k[6].imag(), 0.0, doublethreshold); + EXPECT_NEAR(psi_k[7].real(), 0.0, doublethreshold); + EXPECT_NEAR(psi_k[7].imag(), 0.0, doublethreshold); + EXPECT_NEAR(psi_k[8].real(), 0.948683298050514, doublethreshold); + EXPECT_NEAR(psi_k[8].imag(), 0.0, doublethreshold); + EXPECT_NEAR(psi_k[9].real(), 0.0, doublethreshold); + EXPECT_NEAR(psi_k[9].imag(), 0.0, doublethreshold); + EXPECT_NEAR(psi_k[10].real(), 0.0, doublethreshold); + EXPECT_NEAR(psi_k[10].imag(), 0.0, doublethreshold); + EXPECT_NEAR(psi_k[11].real(), 0.316227766016838, doublethreshold); + EXPECT_NEAR(psi_k[11].imag(), 0.0, doublethreshold); + + delete[] psi_k; + delete[] Stmp; +} diff --git a/source/module_hamilt_lcao/module_tddft/test/propagator_test1.cpp b/source/module_hamilt_lcao/module_tddft/test/propagator_test1.cpp new file mode 100644 index 0000000000..ce17a13185 --- /dev/null +++ b/source/module_hamilt_lcao/module_tddft/test/propagator_test1.cpp @@ -0,0 +1,114 @@ +#include +#include +#include + +#include "module_basis/module_ao/parallel_orbitals.h" +#include "module_hamilt_lcao/module_tddft/propagator.h" +#include "module_io/input.h" +#include "tddft_test.h" + +/************************************************ + * unit test of functions in propagator.h + ***********************************************/ + +/** + * - Tested Function + * - Propagator::compute_propagator_cn2 + * - compute propagator of method Crank-Nicolson. + */ + +Input INPUT; +#define doublethreshold 1e-8 + +Parallel_Orbitals::Parallel_Orbitals() +{ +} +Parallel_Orbitals::~Parallel_Orbitals() +{ +} + +TEST(PropagatorTest, testPropagatorCN) +{ + std::complex* U_operator; + std::complex* Stmp; + std::complex* Htmp; + int nlocal = 4; + bool print_matrix = false; + Parallel_Orbitals* pv; + pv = new Parallel_Orbitals(); + pv->nloc = nlocal * nlocal; + pv->ncol = nlocal; + INPUT.mdp.md_dt = 4; + + // Initialize input matrices + int info; + int mb = 1, nb = 1, lda = nlocal, ldc = nlocal; + int irsrc = 0, icsrc = 0, lld = numroc_(&nlocal, &mb, &myprow, &irsrc, &nprow); + descinit_(pv->desc, &nlocal, &nlocal, &mb, &nb, &irsrc, &icsrc, &ictxt, &lld, &info); + + // Initialize data + U_operator = new std::complex[nlocal * nlocal]; + Stmp = new std::complex[nlocal * nlocal]; + Htmp = new std::complex[nlocal * nlocal]; + + for (int i = 0; i < nlocal * nlocal; ++i) + { + U_operator[i] = std::complex(0.0, 0.0); + } + for (int i = 0; i < nlocal; ++i) + { + for (int j = 0; j < nlocal; ++j) + { + if (i == j) + { + Htmp[i * nlocal + j] = std::complex(1.0, 0.0); + Stmp[i * nlocal + j] = std::complex(1.0, 0.0); + } + } + } + Stmp[1] = 0.5; + Stmp[4] = 0.5; + + // Call the function + int propagator = 0; + module_tddft::Propagator prop(propagator, pv); + prop.compute_propagator(nlocal, Stmp, Htmp, nullptr, U_operator, print_matrix); + + // Check the results + EXPECT_NEAR(U_operator[0].real(), -0.107692307692308, doublethreshold); + EXPECT_NEAR(U_operator[0].imag(), -0.861538461538462, doublethreshold); + EXPECT_NEAR(U_operator[1].real(), 0.492307692307692, doublethreshold); + EXPECT_NEAR(U_operator[1].imag(), -0.0615384615384615, doublethreshold); + EXPECT_NEAR(U_operator[2].real(), 0.0, doublethreshold); + EXPECT_NEAR(U_operator[2].imag(), 0.0, doublethreshold); + EXPECT_NEAR(U_operator[3].real(), 0.0, doublethreshold); + EXPECT_NEAR(U_operator[3].imag(), 0.0, doublethreshold); + EXPECT_NEAR(U_operator[4].real(), 0.492307692307692, doublethreshold); + EXPECT_NEAR(U_operator[4].imag(), -0.0615384615384615, doublethreshold); + EXPECT_NEAR(U_operator[5].real(), -0.107692307692308, doublethreshold); + EXPECT_NEAR(U_operator[5].imag(), -0.861538461538462, doublethreshold); + EXPECT_NEAR(U_operator[6].real(), 0.0, doublethreshold); + EXPECT_NEAR(U_operator[6].imag(), 0.0, doublethreshold); + EXPECT_NEAR(U_operator[7].real(), 0.0, doublethreshold); + EXPECT_NEAR(U_operator[7].imag(), 0.0, doublethreshold); + EXPECT_NEAR(U_operator[8].real(), 0.0, doublethreshold); + EXPECT_NEAR(U_operator[8].imag(), 0.0, doublethreshold); + EXPECT_NEAR(U_operator[9].real(), 0.0, doublethreshold); + EXPECT_NEAR(U_operator[9].imag(), 0.0, doublethreshold); + EXPECT_NEAR(U_operator[10].real(), 0.0, doublethreshold); + EXPECT_NEAR(U_operator[10].imag(), -1.0, doublethreshold); + EXPECT_NEAR(U_operator[11].real(), 0.0, doublethreshold); + EXPECT_NEAR(U_operator[11].imag(), 0.0, doublethreshold); + EXPECT_NEAR(U_operator[12].real(), 0.0, doublethreshold); + EXPECT_NEAR(U_operator[12].imag(), 0.0, doublethreshold); + EXPECT_NEAR(U_operator[13].real(), 0.0, doublethreshold); + EXPECT_NEAR(U_operator[13].imag(), 0.0, doublethreshold); + EXPECT_NEAR(U_operator[14].real(), 0.0, doublethreshold); + EXPECT_NEAR(U_operator[14].imag(), 0.0, doublethreshold); + EXPECT_NEAR(U_operator[15].real(), 0.0, doublethreshold); + EXPECT_NEAR(U_operator[15].imag(), -1.0, doublethreshold); + + delete[] U_operator; + delete[] Htmp; + delete[] Stmp; +} \ No newline at end of file diff --git a/source/module_hamilt_lcao/module_tddft/test/propagator_test2.cpp b/source/module_hamilt_lcao/module_tddft/test/propagator_test2.cpp new file mode 100644 index 0000000000..e015e17f45 --- /dev/null +++ b/source/module_hamilt_lcao/module_tddft/test/propagator_test2.cpp @@ -0,0 +1,119 @@ +#include +#include +#include + +#include "module_basis/module_ao/parallel_orbitals.h" +#include "module_hamilt_lcao/module_tddft/propagator.h" +#include "module_io/input.h" +#include "tddft_test.h" + +/************************************************ + * unit test of functions in propagator.h + ***********************************************/ + +/** + * - Tested Function + * - Propagator::compute_propagator_taylor + * - compute propagator of method 4th Taylor expansion. + */ + +#define doublethreshold 1e-8 + +TEST(PropagatorTest, testPropagatorTaylor) +{ + std::complex* U_operator; + std::complex* Stmp; + std::complex* Htmp; + int nlocal = 4; + bool print_matrix = false; + Parallel_Orbitals* pv; + pv = new Parallel_Orbitals(); + pv->nloc = nlocal * nlocal; + pv->ncol = nlocal; + pv->nrow = nlocal; + pv->dim0 = 1; + pv->dim1 = 1; + pv->nb = 1; + + int dim[2]; + int period[2] = {1, 1}; + int reorder = 0; + dim[0] = nprow; + dim[1] = npcol; + + MPI_Cart_create(MPI_COMM_WORLD, 2, dim, period, reorder, &pv->comm_2D); + + INPUT.mdp.md_dt = 4; + + // Initialize input matrices + int info; + int mb = 1, nb = 1, lda = nlocal, ldc = nlocal; + int irsrc = 0, icsrc = 0, lld = numroc_(&nlocal, &mb, &myprow, &irsrc, &nprow); + descinit_(pv->desc, &nlocal, &nlocal, &mb, &nb, &irsrc, &icsrc, &ictxt, &lld, &info); + + // Initialize data + U_operator = new std::complex[nlocal * nlocal]; + Stmp = new std::complex[nlocal * nlocal]; + Htmp = new std::complex[nlocal * nlocal]; + + for (int i = 0; i < nlocal * nlocal; ++i) + { + U_operator[i] = std::complex(0.0, 0.0); + } + for (int i = 0; i < nlocal; ++i) + { + for (int j = 0; j < nlocal; ++j) + { + if (i == j) + { + Htmp[i * nlocal + j] = std::complex(1.0, 0.0); + Stmp[i * nlocal + j] = std::complex(1.0, 0.0); + } + } + } + Stmp[1] = 0.5; + Stmp[4] = 0.5; + + // Call the function + int propagator = 1; + module_tddft::Propagator prop(propagator, pv); + prop.compute_propagator(nlocal, Stmp, Htmp, nullptr, U_operator, print_matrix); + + // Check the results + EXPECT_NEAR(U_operator[0].real(), 1.95473251028807, doublethreshold); + EXPECT_NEAR(U_operator[0].imag(), 2.8641975308642, doublethreshold); + EXPECT_NEAR(U_operator[1].real(), -1.7119341563786, doublethreshold); + EXPECT_NEAR(U_operator[1].imag(), -3.80246913580247, doublethreshold); + EXPECT_NEAR(U_operator[2].real(), 0.0, doublethreshold); + EXPECT_NEAR(U_operator[2].imag(), 0.0, doublethreshold); + EXPECT_NEAR(U_operator[3].real(), 0.0, doublethreshold); + EXPECT_NEAR(U_operator[3].imag(), 0.0, doublethreshold); + EXPECT_NEAR(U_operator[4].real(), -1.7119341563786, doublethreshold); + EXPECT_NEAR(U_operator[4].imag(), -3.80246913580247, doublethreshold); + EXPECT_NEAR(U_operator[5].real(), 1.95473251028807, doublethreshold); + EXPECT_NEAR(U_operator[5].imag(), 2.8641975308642, doublethreshold); + EXPECT_NEAR(U_operator[6].real(), 0.0, doublethreshold); + EXPECT_NEAR(U_operator[6].imag(), 0.0, doublethreshold); + EXPECT_NEAR(U_operator[7].real(), 0.0, doublethreshold); + EXPECT_NEAR(U_operator[7].imag(), 0.0, doublethreshold); + EXPECT_NEAR(U_operator[8].real(), 0.0, doublethreshold); + EXPECT_NEAR(U_operator[8].imag(), 0.0, doublethreshold); + EXPECT_NEAR(U_operator[9].real(), 0.0, doublethreshold); + EXPECT_NEAR(U_operator[9].imag(), 0.0, doublethreshold); + EXPECT_NEAR(U_operator[10].real(), -0.333333333333333, doublethreshold); + EXPECT_NEAR(U_operator[10].imag(), -0.666666666666667, doublethreshold); + EXPECT_NEAR(U_operator[11].real(), 0.0, doublethreshold); + EXPECT_NEAR(U_operator[11].imag(), 0.0, doublethreshold); + EXPECT_NEAR(U_operator[12].real(), 0.0, doublethreshold); + EXPECT_NEAR(U_operator[12].imag(), 0.0, doublethreshold); + EXPECT_NEAR(U_operator[13].real(), 0.0, doublethreshold); + EXPECT_NEAR(U_operator[13].imag(), 0.0, doublethreshold); + EXPECT_NEAR(U_operator[14].real(), 0.0, doublethreshold); + EXPECT_NEAR(U_operator[14].imag(), 0.0, doublethreshold); + EXPECT_NEAR(U_operator[15].real(), -0.333333333333333, doublethreshold); + EXPECT_NEAR(U_operator[15].imag(), -0.666666666666667, doublethreshold); + + delete[] U_operator; + delete[] Htmp; + delete[] Stmp; +} \ No newline at end of file diff --git a/source/module_hamilt_lcao/module_tddft/test/propagator_test3.cpp b/source/module_hamilt_lcao/module_tddft/test/propagator_test3.cpp new file mode 100644 index 0000000000..4023b518fe --- /dev/null +++ b/source/module_hamilt_lcao/module_tddft/test/propagator_test3.cpp @@ -0,0 +1,124 @@ +#include +#include +#include + +#include "module_basis/module_ao/parallel_orbitals.h" +#include "module_hamilt_lcao/module_tddft/propagator.h" +#include "module_io/input.h" +#include "tddft_test.h" + +/************************************************ + * unit test of functions in propagator.h + ***********************************************/ + +/** + * - Tested Function + * - Propagator::compute_propagator_etrs + * - compute propagator of method ETRS. + */ + +#define doublethreshold 1e-8 + +TEST(PropagatorTest, testPropagatorETRS) +{ + std::complex* U_operator; + std::complex* Stmp; + std::complex* Htmp; + std::complex* Hlaststep; + int nlocal = 4; + bool print_matrix = false; + Parallel_Orbitals* pv; + pv = new Parallel_Orbitals(); + pv->nloc = nlocal * nlocal; + pv->ncol = nlocal; + pv->nrow = nlocal; + pv->dim0 = 1; + pv->dim1 = 1; + pv->nb = 1; + + int dim[2]; + int period[2] = {1, 1}; + int reorder = 0; + dim[0] = nprow; + dim[1] = npcol; + + MPI_Cart_create(MPI_COMM_WORLD, 2, dim, period, reorder, &pv->comm_2D); + + INPUT.mdp.md_dt = 4; + + // Initialize input matrices + int info; + int mb = 1, nb = 1, lda = nlocal, ldc = nlocal; + int irsrc = 0, icsrc = 0, lld = numroc_(&nlocal, &mb, &myprow, &irsrc, &nprow); + descinit_(pv->desc, &nlocal, &nlocal, &mb, &nb, &irsrc, &icsrc, &ictxt, &lld, &info); + + // Initialize data + U_operator = new std::complex[nlocal * nlocal]; + Stmp = new std::complex[nlocal * nlocal]; + Htmp = new std::complex[nlocal * nlocal]; + Hlaststep = new std::complex[nlocal * nlocal]; + + for (int i = 0; i < nlocal * nlocal; ++i) + { + U_operator[i] = std::complex(0.0, 0.0); + } + for (int i = 0; i < nlocal; ++i) + { + for (int j = 0; j < nlocal; ++j) + { + if (i == j) + { + Htmp[i * nlocal + j] = std::complex(1.0, 0.0); + Hlaststep[i * nlocal + j] = std::complex(1.0, 0.0); + Stmp[i * nlocal + j] = std::complex(1.0, 0.0); + } + } + } + Stmp[1] = 0.5; + Stmp[4] = 0.5; + Hlaststep[0] = 1.1; + + // Call the function + int propagator = 2; + module_tddft::Propagator prop(propagator, pv); + prop.compute_propagator(nlocal, Stmp, Htmp, Hlaststep, U_operator, print_matrix); + + // Check the results + EXPECT_NEAR(U_operator[0].real(), -0.0105865569272976, doublethreshold); + EXPECT_NEAR(U_operator[0].imag(), -0.228336412132297, doublethreshold); + EXPECT_NEAR(U_operator[1].real(), 0.195527023319616, doublethreshold); + EXPECT_NEAR(U_operator[1].imag(), -0.728701844231062, doublethreshold); + EXPECT_NEAR(U_operator[2].real(), 0.0, doublethreshold); + EXPECT_NEAR(U_operator[2].imag(), 0.0, doublethreshold); + EXPECT_NEAR(U_operator[3].real(), 0.0, doublethreshold); + EXPECT_NEAR(U_operator[3].imag(), 0.0, doublethreshold); + EXPECT_NEAR(U_operator[4].real(), 0.247662246608749, doublethreshold); + EXPECT_NEAR(U_operator[4].imag(), -0.694263679317177, doublethreshold); + EXPECT_NEAR(U_operator[5].real(), -0.0219180003048316, doublethreshold); + EXPECT_NEAR(U_operator[5].imag(), -0.300090839811004, doublethreshold); + EXPECT_NEAR(U_operator[6].real(), 0.0, doublethreshold); + EXPECT_NEAR(U_operator[6].imag(), 0.0, doublethreshold); + EXPECT_NEAR(U_operator[7].real(), 0.0, doublethreshold); + EXPECT_NEAR(U_operator[7].imag(), 0.0, doublethreshold); + EXPECT_NEAR(U_operator[8].real(), 0.0, doublethreshold); + EXPECT_NEAR(U_operator[8].imag(), 0.0, doublethreshold); + EXPECT_NEAR(U_operator[9].real(), 0.0, doublethreshold); + EXPECT_NEAR(U_operator[9].imag(), 0.0, doublethreshold); + EXPECT_NEAR(U_operator[10].real(), -0.401041666666667, doublethreshold); + EXPECT_NEAR(U_operator[10].imag(), -0.902777777777778, doublethreshold); + EXPECT_NEAR(U_operator[11].real(), 0.0, doublethreshold); + EXPECT_NEAR(U_operator[11].imag(), 0.0, doublethreshold); + EXPECT_NEAR(U_operator[12].real(), 0.0, doublethreshold); + EXPECT_NEAR(U_operator[12].imag(), 0.0, doublethreshold); + EXPECT_NEAR(U_operator[13].real(), 0.0, doublethreshold); + EXPECT_NEAR(U_operator[13].imag(), 0.0, doublethreshold); + EXPECT_NEAR(U_operator[14].real(), 0.0, doublethreshold); + EXPECT_NEAR(U_operator[14].imag(), 0.0, doublethreshold); + EXPECT_NEAR(U_operator[15].real(), -0.401041666666667, doublethreshold); + EXPECT_NEAR(U_operator[15].imag(), -0.902777777777778, doublethreshold); + + delete[] U_operator; + delete[] Htmp; + delete[] Stmp; + delete[] Hlaststep; +} \ No newline at end of file diff --git a/source/module_hamilt_lcao/module_tddft/test/tddft_test.cpp b/source/module_hamilt_lcao/module_tddft/test/tddft_test.cpp new file mode 100644 index 0000000000..fedb46a976 --- /dev/null +++ b/source/module_hamilt_lcao/module_tddft/test/tddft_test.cpp @@ -0,0 +1,50 @@ +#include "tddft_test.h" + +#include +#include + +#include "module_base/blacs_connector.h" +#include "module_base/scalapack_connector.h" +#include "module_basis/module_ao/parallel_orbitals.h" + +/************************************************ + * unit test of module_tddft + ***********************************************/ + +int myprow, nprow, ictxt, mypcol, npcol; + +void MPIInit() +{ + int myrank; + int mysize; + + MPI_Init(NULL, NULL); + MPI_Comm_size(MPI_COMM_WORLD, &mysize); + MPI_Comm_rank(MPI_COMM_WORLD, &myrank); + + // Set up BLACS context and grid + int nprocs; + nprow = 1; + npcol = 1; + Cblacs_pinfo(&myrank, &mysize); + Cblacs_get(-1, 0, &ictxt); + Cblacs_gridinit(&ictxt, "Row", nprow, npcol); + Cblacs_gridinfo(ictxt, &nprow, &npcol, &myprow, &mypcol); +} + +/************************************************ + * unit test of module_tddft + ***********************************************/ +int main(int argc, char** argv) +{ + MPIInit(); + + int result = 0; + testing::InitGoogleTest(&argc, argv); + result = RUN_ALL_TESTS(); + + Cblacs_exit(ictxt); + + // MPI_Finalize(); + return 0; +} diff --git a/source/module_hamilt_lcao/module_tddft/test/tddft_test.h b/source/module_hamilt_lcao/module_tddft/test/tddft_test.h new file mode 100644 index 0000000000..c13f0b0e4e --- /dev/null +++ b/source/module_hamilt_lcao/module_tddft/test/tddft_test.h @@ -0,0 +1,50 @@ +#ifndef __TDDFTTEST +#define __TDDFTTEST +#include + +#include "gtest/gtest.h" +#include "module_basis/module_ao/parallel_orbitals.h" + +using namespace std; +extern int myprow, nprow, ictxt, mypcol, npcol; + +class TDDFTTEST : public testing::Test +{ + public: + static void SetUpTestCase() + { + cout << "\033[32m" + << "============================" + << "\033[0m" << endl; + cout << "\033[32m" + << "= TDDFT MODULE TEST START =" + << "\033[0m" << endl; + cout << "\033[32m" + << "============================" + << "\033[0m" << endl; + } + static void TearDownTestCase() + { + cout << "\033[32m" + << "============================" + << "\033[0m" << endl; + cout << "\033[32m" + << "= TDDFT MODULE TEST END =" + << "\033[0m" << endl; + cout << "\033[32m" + << "============================" + << "\033[0m" << endl; + } + void SetUp() + { + cout << "\033[32m" + << "[ CASE ]" + << "\033[0m" + << " "; + } + void TearDown() + { + } +}; + +#endif diff --git a/source/module_hamilt_lcao/module_tddft/test/upsi_test1.cpp b/source/module_hamilt_lcao/module_tddft/test/upsi_test1.cpp new file mode 100644 index 0000000000..e73221f80b --- /dev/null +++ b/source/module_hamilt_lcao/module_tddft/test/upsi_test1.cpp @@ -0,0 +1,76 @@ +#include +#include +#include + +#include "module_basis/module_ao/parallel_orbitals.h" +#include "module_hamilt_lcao/module_tddft/upsi.h" +#include "tddft_test.h" + +#define doublethreshold 1e-8 + +/************************************************ + * unit test of functions in upsi.h + ***********************************************/ + +/** + * - Tested Function + * - upsi + * - apply U_operator to the wave function of the previous step for new wave function. + */ + +Parallel_Orbitals::Parallel_Orbitals() +{ +} +Parallel_Orbitals::~Parallel_Orbitals() +{ +} + +TEST(UpsiTest, testUpsi1) +{ + std::complex* U_operator; + std::complex* psi_k_laststep; + std::complex* psi_k; + int nband = 2; + int nlocal = 2; + bool print_matrix = false; + Parallel_Orbitals* pv; + pv = new Parallel_Orbitals(); + pv->nloc = nlocal * nband; + pv->ncol = nlocal; + pv->ncol_bands = nband; + + // Initialize input matrices + int info; + int mb = 1, nb = 1, lda = nband, ldc = nlocal; + int irsrc = 0, icsrc = 0, lld = numroc_(&nlocal, &mb, &myprow, &irsrc, &nprow); + descinit_(pv->desc, &nlocal, &nlocal, &mb, &nb, &irsrc, &icsrc, &ictxt, &lld, &info); + descinit_(pv->desc_wfc, &nlocal, &nband, &mb, &nb, &irsrc, &icsrc, &ictxt, &lld, &info); + + // Initialize data + U_operator = new std::complex[nlocal * nlocal]; + psi_k_laststep = new std::complex[nlocal * nband]; + psi_k = new std::complex[nlocal * nband]; + + for (int i = 0; i < nlocal * nlocal; ++i) + { + U_operator[i] = std::complex(1.0, 0.0); + } + for (int i = 0; i < nlocal * nband; ++i) + { + psi_k_laststep[i] = std::complex(1.0, 0.0); + psi_k[i] = std::complex(0.0, 0.0); + } + + // Call the function + module_tddft::upsi(pv, nband, nlocal, U_operator, psi_k_laststep, psi_k, print_matrix); + + // Check the results + for (int i = 0; i < nlocal * nband; ++i) + { + EXPECT_NEAR(psi_k[i].real(), nlocal, doublethreshold); + EXPECT_NEAR(psi_k[i].imag(), 0.0, doublethreshold); + } + delete[] U_operator; + delete[] psi_k; + delete[] psi_k_laststep; +} diff --git a/source/module_hamilt_lcao/module_tddft/test/upsi_test2.cpp b/source/module_hamilt_lcao/module_tddft/test/upsi_test2.cpp new file mode 100644 index 0000000000..33ca60f7a6 --- /dev/null +++ b/source/module_hamilt_lcao/module_tddft/test/upsi_test2.cpp @@ -0,0 +1,69 @@ +#include +#include +#include + +#include "module_basis/module_ao/parallel_orbitals.h" +#include "module_hamilt_lcao/module_tddft/upsi.h" +#include "tddft_test.h" + +#define doublethreshold 1e-8 + +/************************************************ + * unit test of functions in upsi.h + ***********************************************/ + +/** + * - Tested Function + * - upsi + * - apply U_operator to the wave function of the previous step for new wave function. + */ + +TEST(UpsiTest, testUpsi2) +{ + std::complex* U_operator; + std::complex* psi_k_laststep; + std::complex* psi_k; + int nband = 3; + int nlocal = 4; + bool print_matrix = false; + Parallel_Orbitals* pv; + pv = new Parallel_Orbitals(); + pv->nloc = nlocal * nband; + pv->ncol = nlocal; + pv->ncol_bands = nband; + + // Initialize input matrices + int info; + int mb = 1, nb = 1, lda = nband, ldc = nlocal; + int irsrc = 0, icsrc = 0, lld = numroc_(&nlocal, &mb, &myprow, &irsrc, &nprow); + descinit_(pv->desc, &nlocal, &nlocal, &mb, &nb, &irsrc, &icsrc, &ictxt, &lld, &info); + descinit_(pv->desc_wfc, &nlocal, &nband, &mb, &nb, &irsrc, &icsrc, &ictxt, &lld, &info); + + // Initialize data + U_operator = new std::complex[nlocal * nlocal]; + psi_k_laststep = new std::complex[nlocal * nband]; + psi_k = new std::complex[nlocal * nband]; + + for (int i = 0; i < nlocal * nlocal; ++i) + { + U_operator[i] = std::complex(1.0, 0.0); + } + for (int i = 0; i < nlocal * nband; ++i) + { + psi_k_laststep[i] = std::complex(2.0, 0.0); + psi_k[i] = std::complex(0.0, 0.0); + } + + // Call the function + module_tddft::upsi(pv, nband, nlocal, U_operator, psi_k_laststep, psi_k, print_matrix); + + // Check the results + for (int i = 0; i < nlocal * nband; ++i) + { + EXPECT_NEAR(psi_k[i].real(), 2.0 * nlocal, doublethreshold); + EXPECT_NEAR(psi_k[i].imag(), 0.0, doublethreshold); + } + delete[] U_operator; + delete[] psi_k; + delete[] psi_k_laststep; +} diff --git a/source/module_hamilt_lcao/module_tddft/test/upsi_test3.cpp b/source/module_hamilt_lcao/module_tddft/test/upsi_test3.cpp new file mode 100644 index 0000000000..576adf60ee --- /dev/null +++ b/source/module_hamilt_lcao/module_tddft/test/upsi_test3.cpp @@ -0,0 +1,86 @@ +#include +#include +#include + +#include "module_basis/module_ao/parallel_orbitals.h" +#include "module_hamilt_lcao/module_tddft/upsi.h" +#include "tddft_test.h" + +#define doublethreshold 1e-8 + +/************************************************ + * unit test of functions in upsi.h + ***********************************************/ + +/** + * - Tested Function + * - upsi + * - apply U_operator to the wave function of the previous step for new wave function. + */ + +TEST(UpsiTest, testUpsi3) +{ + std::complex* U_operator; + std::complex* psi_k_laststep; + std::complex* psi_k; + int nband = 3; + int nlocal = 4; + bool print_matrix = false; + Parallel_Orbitals* pv; + pv = new Parallel_Orbitals(); + pv->nloc = nlocal * nband; + pv->ncol = nlocal; + pv->ncol_bands = nband; + + // Initialize input matrices + int info; + int mb = 1, nb = 1, lda = nband, ldc = nlocal; + int irsrc = 0, icsrc = 0, lld = numroc_(&nlocal, &mb, &myprow, &irsrc, &nprow); + descinit_(pv->desc, &nlocal, &nlocal, &mb, &nb, &irsrc, &icsrc, &ictxt, &lld, &info); + descinit_(pv->desc_wfc, &nlocal, &nband, &mb, &nb, &irsrc, &icsrc, &ictxt, &lld, &info); + + // Initialize data + U_operator = new std::complex[nlocal * nlocal]; + psi_k_laststep = new std::complex[nlocal * nband]; + psi_k = new std::complex[nlocal * nband]; + + for (int i = 0; i < nlocal; ++i) + { + for (int j = 0; j < nlocal; ++j) + { + if (i == j) + { + U_operator[i * nlocal + j] = std::complex(1.0, 0.0); + } + else + { + U_operator[i * nlocal + j] = std::complex(0.0, 0.0); + } + } + } + + for (int i = 0; i < nlocal; ++i) + { + for (int j = 0; j < nband; ++j) + { + psi_k_laststep[i * nband + j] = std::complex(i * nband + j, 0.0); + psi_k[i * nband + j] = std::complex(0.0, 0.0); + } + } + + // Call the function + module_tddft::upsi(pv, nband, nlocal, U_operator, psi_k_laststep, psi_k, print_matrix); + + // Check the results + for (int i = 0; i < nlocal; ++i) + { + for (int j = 0; j < nband; ++j) + { + EXPECT_NEAR(psi_k[i * nband + j].real(), i * nband + j, doublethreshold); + EXPECT_NEAR(psi_k[i].imag(), 0.0, doublethreshold); + } + } + delete[] U_operator; + delete[] psi_k; + delete[] psi_k_laststep; +} From 65707408108eb4ea7c71355f111c734f0e345009 Mon Sep 17 00:00:00 2001 From: Yike Huang <67682086+kirk0830@users.noreply.github.com> Date: Fri, 9 Jun 2023 10:36:46 +0800 Subject: [PATCH 11/13] Test: add module_io/bessel_basis unit test (#2614) * Test: module_io/bessel_basis unit test still leave readin_C4 and allocate_C4 two functions untested, will complete soon. * Update bessel_basis_test.cpp add annotations to describe tested functions --- source/module_io/bessel_basis.cpp | 36 +- source/module_io/bessel_basis.h | 140 ++++- source/module_io/numerical_basis.cpp | 8 +- source/module_io/numerical_descriptor.cpp | 4 +- source/module_io/test/CMakeLists.txt | 6 + source/module_io/test/INPUTs | 3 + source/module_io/test/bessel_basis_test.cpp | 557 ++++++++++++++++++ .../BesselBasis_UnitTest_C4_AtomType0.html | 25 + 8 files changed, 728 insertions(+), 51 deletions(-) create mode 100644 source/module_io/test/INPUTs create mode 100644 source/module_io/test/bessel_basis_test.cpp create mode 100644 source/module_io/test/support/BesselBasis_UnitTest_C4_AtomType0.html diff --git a/source/module_io/bessel_basis.cpp b/source/module_io/bessel_basis.cpp index 421f8e5318..e82d39baec 100644 --- a/source/module_io/bessel_basis.cpp +++ b/source/module_io/bessel_basis.cpp @@ -1,5 +1,5 @@ -#include "bessel_basis.h" #include "module_hamilt_pw/hamilt_pwdft/global.h" +#include "bessel_basis.h" #include "module_base/parallel_common.h" #include "module_base/math_integral.h" #include "module_base/math_sphbes.h" @@ -25,9 +25,11 @@ void Bessel_Basis::init( const bool &smooth, const double &sigma, const double &rcut_in, - const double &tol_in, + const double &tol_in, + const UnitCell& ucell, const double &dk, - const double &dr) + const double &dr + ) { ModuleBase::TITLE("Bessel_Basis", "init"); this->Dk = dk; @@ -63,14 +65,13 @@ void Bessel_Basis::init( if( start_from_file ) { // setup C4 - this->allocate_C4(ntype, lmax_in, GlobalC::ucell.nmax, Ecut_number); - + this->allocate_C4(ntype, lmax_in, ucell.nmax, Ecut_number, ucell); // check tolerence - this->readin_C4("INPUTs", ntype, ecut, rcut, Ecut_number, tolerence); + this->readin_C4("INPUTs", ntype, ecut, rcut, Ecut_number, tolerence, ucell); #ifdef __MPI Parallel_Common::bcast_double( C4.ptr, C4.getSize() ); #endif - this->init_Faln(ntype, lmax_in, GlobalC::ucell.nmax, Ecut_number); + this->init_Faln(ntype, lmax_in, ucell.nmax, Ecut_number, ucell); } return; @@ -124,7 +125,8 @@ void Bessel_Basis::init_Faln( const int &ntype, const int &lmax, const int &nmax, - const int &ecut_number) + const int &ecut_number, + const UnitCell& ucell) { ModuleBase::TITLE("Bessel_Basis","init_Faln"); ModuleBase::timer::tick("Spillage","init_Faln"); @@ -135,9 +137,9 @@ void Bessel_Basis::init_Faln( this->nwfc = 0; for(int it=0; it> filec4; - for(int il=0; il< GlobalC::ucell.atoms[it].nwl+1; il++) + for(int il=0; il< ucell.atoms[it].nwl+1; il++) { - for(int in=0; in< GlobalC::ucell.atoms[it].l_nchi[il]; in++) + for(int in=0; in< ucell.atoms[it].l_nchi[il]; in++) { //for tests //std::cout << "\n" << std::setw(5) << it << std::setw(5) << il << std::setw(5) << in; @@ -448,7 +451,8 @@ void Bessel_Basis::allocate_C4( const int &ntype, const int &lmax, const int &nmax, - const int &ecut_number) + const int &ecut_number, + const UnitCell& ucell) { ModuleBase::TITLE("Bessel_Basis","allocate_C4"); @@ -456,9 +460,9 @@ void Bessel_Basis::allocate_C4( for(int it=0; it=rcut, j_l(q*r) = 0 (if not smoothed) double rcut; + /// @brief energy cutoff for determining kmesh and number of SBFs double ecut; double tolerence; - - bool smooth; // mohan add 2009-01-18 + /// @brief whether smooth SBFs around cutoff radius, resulting in non-zero values. For importance of smooth of SBFs, see J. Phys.: Condens. Matter 22 (2010) 445501, eqn 6. (mohan add 2009-01-18) + bool smooth; + /// @brief stddev of smooth function (Gaussian function, centered at rcut) double sigma; - //======================================================== - // MEMBER FUNCTIONS : - // NAME : readin ( read in parameters ) - //======================================================== + /// @brief Allocate memory for C4 matrix and initialize all elements to one. + /// @param ntype number of atom types + /// @param lmax maximal angular momentum of localized orbitals + /// @param nmax maximal principal quantum number of localized orbitals + /// @param ecut_number number of SBFs void allocate_C4( const int &ntype, const int &lmax, const int &nmax, - const int &ecut_number); - + const int &ecut_number, + const UnitCell& ucell + ); + + /// @brief Read C4 from external file. Presently an O(N^2) search algorithm is used. A HTML parser is needed in the future to improve performance. + /// @param name name of external file where C4-stored file information is contained + /// @param ntype number of atom types + /// @param ecut energy cutoff + /// @param rcut cutoff radius + /// @param ecut_number number of SBFs + /// @param tolerence accurancy of SBFs, here only used for consistency check void readin_C4( const std::string &name, const int &ntype, const int &ecut, const int &rcut, const int &ecut_number, - const double &tolerence); + const double &tolerence, + const UnitCell& ucell + ); void init_TableOne(void); - void init_Faln( const int &ntype,const int &lmax,const int &nmax,const int &ecut_number); - int nwfc; // number of localized wave functions - - // init table, used for output overlap Q. + /// @brief calculate F_{aln}(it, il, in, ik) = sum_{ie}{C4(it, il, in, ie)*TableOne(il, ie, ik)}, where TableOne is overlap integral between two spherical bessel functions (jle(r) and jlk(r)) + /// @param ntype number of atomtype + /// @param lmax maximal angular momentum + /// @param nmax maximal chi + /// @param ecut_number number of SBFs + void init_Faln( + const int &ntype, + const int &lmax, + const int &nmax, + const int &ecut_number, + const UnitCell& ucell + ); + + /// @brief number of localized wave functions + int nwfc; + + /// @brief calculate element value of TableOne matrix + /// @details (be called in Bessel_Basis::init(), used for outputing overlap Q matrix) initialize the table whose matrix element is the result of integral int{dr r^2 jle(r)*jlk(r)}, TableOne has three subscript (l, ie, ik), the first runs over orbitals' angular momentum and ie, ik run over ecut_number and kmesh SBFs + /// @param smooth_in whether jle(r) SBF is smoothed by a Gaussian function + /// @param sigma_in stddev for controlling smearing of Gaussian function for smoothing jle(r) + /// @param ecutwfc planewave kinetic energy cutoff for controlling kspace sampling + /// @param rcut cutoff radius of SBFs + /// @param dr realspace grid + /// @param dk kspace grid + /// @param lmax maximal angular momentum for SBFs + /// @param ecut_number number of SBFs + /// @param tolerence accurancy of SBFs void init_TableOne( const bool smooth_in, const double &sigma_in, diff --git a/source/module_io/numerical_basis.cpp b/source/module_io/numerical_basis.cpp index ad4bc0eaf4..36cde008e6 100644 --- a/source/module_io/numerical_basis.cpp +++ b/source/module_io/numerical_basis.cpp @@ -43,7 +43,9 @@ void Numerical_Basis::start_from_file_k(const int& ik, ModuleBase::ComplexMatrix INPUT.bessel_nao_smooth, INPUT.bessel_nao_sigma, INPUT.bessel_nao_rcut, - INPUT.bessel_nao_tolerence ); + INPUT.bessel_nao_tolerence, + GlobalC::ucell + ); this->mu_index = this->init_mu_index(); this->init_label = true; } @@ -71,7 +73,9 @@ void Numerical_Basis::output_overlap(const psi::Psi>& psi, INPUT.bessel_nao_smooth, INPUT.bessel_nao_sigma, INPUT.bessel_nao_rcut, - INPUT.bessel_nao_tolerence ); + INPUT.bessel_nao_tolerence, + GlobalC::ucell + ); this->mu_index = this->init_mu_index(); this->init_label = true; } diff --git a/source/module_io/numerical_descriptor.cpp b/source/module_io/numerical_descriptor.cpp index fe91ecd690..9e269745ac 100644 --- a/source/module_io/numerical_descriptor.cpp +++ b/source/module_io/numerical_descriptor.cpp @@ -53,7 +53,9 @@ void Numerical_Descriptor::output_descriptor(const psi::Psi INPUT.bessel_descriptor_smooth, INPUT.bessel_descriptor_sigma, rcut_in, - tol_in ); + tol_in, + GlobalC::ucell + ); this->nmax = Numerical_Descriptor::bessel_basis.get_ecut_number(); this->init_mu_index(); this->init_label = true; diff --git a/source/module_io/test/CMakeLists.txt b/source/module_io/test/CMakeLists.txt index 358dbe1c38..de2c90c48a 100644 --- a/source/module_io/test/CMakeLists.txt +++ b/source/module_io/test/CMakeLists.txt @@ -111,3 +111,9 @@ AddTest( TARGET io_parse_args SOURCES parse_args_test.cpp ) + +AddTest( + TARGET io_bessel_basis_test + LIBS ${math_libs} base device + SOURCES bessel_basis_test.cpp ../bessel_basis.cpp +) \ No newline at end of file diff --git a/source/module_io/test/INPUTs b/source/module_io/test/INPUTs new file mode 100644 index 0000000000..618420b3b5 --- /dev/null +++ b/source/module_io/test/INPUTs @@ -0,0 +1,3 @@ + +./support/BesselBasis_UnitTest_C4_AtomType0.html + diff --git a/source/module_io/test/bessel_basis_test.cpp b/source/module_io/test/bessel_basis_test.cpp new file mode 100644 index 0000000000..1472cfd630 --- /dev/null +++ b/source/module_io/test/bessel_basis_test.cpp @@ -0,0 +1,557 @@ +/// @details unit test for module_io/bessel_basis, not tested functions: readin_C4, allocate_C4 + +#include +#include + +#include + +#include +#include +#include +#include + +#include + +#include "../bessel_basis.h" +#include "../../module_cell/unitcell.h" +#include "../../module_elecstate/magnetism.h" + +#ifdef __LCAO +#include "../../module_cell/setup_nonlocal.h" +#endif +#include "gtest/gtest.h" + +/// @brief Simpson integral +/// @attention this function is a COMPLETE version, but there is no improvement in performance. +/// @param x variable stored in a vector +/// @param y function value stored in a vector +/// @return integral value +double SimpsonIntegral(const std::vector &x, const std::vector &y) +{ + double result = 0.0; + result += (y[0] + y[y.size() - 1]); + /* x and y must have the same size and their length must be a same odd number */ + assert(x.size() == y.size()); + assert(x.size() % 2 == 1); + double h = x[1] - x[0]; + for (int i = 1; i < x.size() - 1; i+=2) + { + result += (2*y[i] + 4*y[i+1]); + } + result *= h/3; + return result; +} +/// @brief recursive definition of 1st-kind Spherical Bessel function of the first kind +/// @attention this function is a COMPLETE version, can replace the one in basic mathematical library +/// @param l order of 1st-kind spherical Bessel function +/// @param x variable +/// @return value of 1st-kind spherical Bessel function +double SphericalBessel(int l, double x) { + /* c++ cannot handle limits well */ + if (x == 0) { + if (l == 0) { + return 1; + } + else { + return 0; + } + } + else { + if (l == 0) { + return std::sin(x) / x; + } + else if (l == 1) { + return (std::sin(x) / (x * x)) - std::cos(x) / x; + } + else { + return ((2 * l - 1) / x) * SphericalBessel(l - 1, x) - SphericalBessel(l - 2, x); + } + } +} +/// @brief 1st-kind Spherical Bessel function of the first kind mapped on realspace grid +/// @attention this function is a COMPLETE version, can replace the one in basic mathematical library +/// @param l order of 1st-kind spherical Bessel function +/// @param q wave vector +/// @param r realspace grid +/// @return function value on realspace grid +std::vector CalculateSphericalBessel(int l, double q, const std::vector& r) { + std::vector bessel; + for (const auto& radius : r) { + double x = q * radius; + double besselValue = SphericalBessel(l, x); + bessel.push_back(besselValue); + } + return bessel; +} +/// @brief Get the first p zeros of 1st-kind, l-ordered Spherical Bessel function from a constant table +/// @attention this function is a INCOMPLETE version, due to limited support of numerical table. +/// @param order l, angular momentum +/// @param number p, number of zeros needed +/// @return a vector of q, the q is from j_l(qr) +std::vector GetSphericalBesselZeros(int order, int number) { + std::map, double> zeros; + + zeros[{0, 1}] = 3.14159; zeros[{0, 2}] = 6.28318; zeros[{0, 3}] = 9.42477; + zeros[{0, 4}] = 12.5664; zeros[{0, 5}] = 15.708; zeros[{0, 6}] = 18.8495; + zeros[{0, 7}] = 21.9911; zeros[{0, 8}] = 25.1327; zeros[{0, 9}] = 28.2743; + zeros[{0, 10}] = 31.4159; + + zeros[{1, 1}] = 4.49341; zeros[{1, 2}] = 7.72525; zeros[{1, 3}] = 10.9041; + zeros[{1, 4}] = 14.0662; zeros[{1, 5}] = 17.2208; zeros[{1, 6}] = 20.3713; + zeros[{1, 7}] = 23.5181; zeros[{1, 8}] = 26.6617; zeros[{1, 9}] = 29.8029; + zeros[{1, 10}] = 32.9425; + + zeros[{2, 1}] = 5.76346; zeros[{2, 2}] = 9.09501; zeros[{2, 3}] = 12.3229; + zeros[{2, 4}] = 15.5146; zeros[{2, 5}] = 18.6861; zeros[{2, 6}] = 21.8457; + zeros[{2, 7}] = 24.9989; zeros[{2, 8}] = 28.1498; zeros[{2, 9}] = 31.2997; + zeros[{2, 10}] = 34.4491; + + zeros[{3, 1}] = 7.01559; zeros[{3, 2}] = 10.4013; zeros[{3, 3}] = 13.5821; + zeros[{3, 4}] = 16.7496; zeros[{3, 5}] = 19.9023; zeros[{3, 6}] = 23.0446; + zeros[{3, 7}] = 26.1799; zeros[{3, 8}] = 29.3105; zeros[{3, 9}] = 32.4377; + zeros[{3, 10}] = 35.5629; + + zeros[{4, 1}] = 8.26356; zeros[{4, 2}] = 11.6209; zeros[{4, 3}] = 14.7965; + zeros[{4, 4}] = 17.9598; zeros[{4, 5}] = 21.113; zeros[{4, 6}] = 24.2583; + zeros[{4, 7}] = 27.3979; zeros[{4, 8}] = 30.5325; zeros[{4, 9}] = 33.6635; + zeros[{4, 10}] = 36.7914; + + zeros[{5, 1}] = 9.51045; zeros[{5, 2}] = 12.8377; zeros[{5, 3}] = 16.0106; + zeros[{5, 4}] = 19.1714; zeros[{5, 5}] = 22.3224; zeros[{5, 6}] = 25.4666; + zeros[{5, 7}] = 28.6055; zeros[{5, 8}] = 31.7408; zeros[{5, 9}] = 34.873; + zeros[{5, 10}] = 38.0025; + + std::vector result; + for (int i = 1; i <= number; ++i) { + result.push_back(zeros[{order, i}]); + } + return result; +} +/// @brief Get mod of q vector of Spherical Bessel functions, all q satisfy when r=`rcut`, j_l(qr)=0. +/// @details first solve the equation j_l(x) = 0, therefore get the table (l, k) -> x, where l is the order of SBF and k is the k-th zero of j_l(x). Then let x = q*rcut, therefore q = x/rcut, return it. +/// @attention this function itself is a COMPLETE version, while the function it called GetSphericalBesselZeros may be INCOMPLETE, due to limited support of numerical table. +/// @param order the angular momentum of Spherical Bessel functions +/// @param number number of q to be returned +/// @param rcut 'cutoff radius' of Spherical Bessel functions. When r=rcut, Spherical Bessel functions are zero, and for r>rcut, they are zero as required by concept of constructing localized atomic orbital. +/// @return a vector of q, the q is from j_l(qr) +std::vector GetqList(int order, int number, double rcut) { + std::vector qList; + std::vector zerosList = GetSphericalBesselZeros(order, number); + for (const auto& zero : zerosList) { + qList.push_back(zero / rcut); + } + return qList; +} +/// @brief overload operator * for vectors, elements will be multiplied one by one +/// @param x one vector +/// @param y another vector +/// @return vector after multiplication +std::vector operator*(const std::vector& x, const std::vector& y) { + std::vector result; + for (int i = 0; i < x.size(); ++i) { + result.push_back(x[i] * y[i]); + } + return result; +} +/// @brief init_TableOne for unit test +/// @attention this is a COMPLETE version of init_TableOne, can replace the one in module_io/bessel_basis.cpp +/// @param smooth whether to smooth the function (gaussian function) +/// @param sigma stddev of gaussian function for smoothing +/// @param ecutwfc control the number of Spheical Bessel functions +/// @param rcut cutoff radius of Spherical Bessel functions, r>rcut, Spherical Bessel functions are zero +/// @param lmax maximal angular momentum of Spherical Bessel functions +/// @param dr grid spacing of r +/// @param dk grid spacing of k +/// @return `std::vector>>` TableOne[angular momentum][index of Spherical Bessel function][Arbitrary wave vector k] +std::vector>> GenerateTableOne(const bool smooth, const double sigma, const double ecutwfc, const double rcut, const int lmax, const double dr, const double dk){ + std::vector rGrid; + std::vector SmoothFactor_rGrid; + int rGridNum = static_cast(rcut/dr) + 4; + if (rGridNum % 2 == 0) + { + rGridNum += 1; + } + for (int indexr=0; indexr kGrid; + int kGridNum = static_cast(sqrt(ecutwfc)/dk) + 4 + 1; + if (kGridNum % 2 == 0) + { + kGridNum += 1; + } + for (double indexk=0; indexk(sqrt(ecutwfc)*rcut/M_PI); + + std::vector>> TableOne; + + for (int l=0; l qList = GetqList(l, NumBesselFunction, rcut); + /* should initialize TableOne first */ + if (l == 0) + { + for (int l=0; l> Zeros2d; + for (int indexq=0; indexq Zeros1d; + for (int indexk=0; indexk jle_rGrid = CalculateSphericalBessel(l, qList[indexq], rGrid); + if (smooth) + { + jle_rGrid = jle_rGrid*SmoothFactor_rGrid; + } + for (int indexk=0; indexk jlk_rGrid = CalculateSphericalBessel(l, kGrid[indexk], rGrid); + std::vector function_rGrid = jlk_rGrid*jle_rGrid*rGrid*rGrid; + TableOne[l][indexq][indexk] = SimpsonIntegral(rGrid, function_rGrid); + } + } + } + return TableOne; +} +/// @brief Improved version of module_io/bessel_basis::readin_C4 and allocate_C4 functions, for generating C4 matrix but now with a higher speed on accessing elements +/// @attention function will read in total number of chi-s both from file and from input, and assert that they are the same. It is also just a INCOMPLETE version, for a complete version, HTML parser library will be included and the other parameter, NumAtomType, will also be used for calibrating data. +/// @param FileName name of external file where C4-stored file information is contained +/// @param NumAtomType number of atom types +/// @param l maximal angular momentum of localized orbitals +/// @param NumChi number of contracted spherical bessel functions to fit one atomic orbital +/// @param NumBesselFunction number of spherical bessel functions in one contracted spherical bessel function +/// @return a map whose key is a string of "atomtype l chi q", and value is the corresponding C4 value +std::unordered_map ReadinC4(const std::string &FileName, const int &NumAtomType, const int &l, const int &NumChi, const int &NumBesselFunction){ + std::unordered_map C4Map; + std::ifstream C4File(FileName); + /* plan to add an HTML parser library in the future... */ + std::string word; + bool b_ReadC4Line = false; + int TotalNumChi = 0; + + assert(!C4File.fail()); + while (C4File.good()) + { + if (!b_ReadC4Line){ + C4File >> word; + if (word == "") + { + b_ReadC4Line = true; + C4File >> word; /* n_chi value read */ + TotalNumChi = std::stoi(word); + assert(TotalNumChi == NumChi); + continue; + } + } + else{ + for (int indexchi = 0; indexchi < TotalNumChi; indexchi++){ + + C4File >> word; + C4File >> word; + C4File >> word; /* skip title1, 2 and 3 */ + + C4File >> word; std::string key = word; key += " "; + C4File >> word; key += word; key += " "; + C4File >> word; key += word; key += " "; + + for (int indexNumBesselFunction = 0; indexNumBesselFunction < NumBesselFunction; indexNumBesselFunction++) + { + std::string keyForMap = key + std::to_string(indexNumBesselFunction); + C4File >> word; + C4Map[keyForMap] = std::stod(word); + } + } + break; + } + } + C4File.close(); + return C4Map; +} +/// @brief Generate F_{alpha, l, chi, k} matrix +/// @param FileName name of external file where C4-stored file information is contained +/// @param NumAtomType number of atom types +/// @param l maximal angular momentum of localized orbitals +/// @param NumChi number of contracted spherical bessel functions to fit one atomic orbital +/// @param NumBesselFunction number of spherical bessel functions in one contracted spherical bessel function +/// @param vvv_d_TableOne Integral table, has subscripts (l, q, k), whose element is the result of integral int{dr r^2 jle(r)*jlk(r)}. l runs over angular momentum, q runs over all spherical bessel functions, k runs over k points sampled controlled by ecutwfc and dk by sqrt(ecutwfc)/dk + 1 + 4. +/// @return F_{alpha, l, chi, k} matrix, whose element is the result of sum_{q}{C4(alpha, l, chi, q)*TableOne(l, q, k)} +std::vector>>> GenerateFaln(const std::string &FileName, const int &NumAtomType, const int &lmax, const int &NumChi, const int &NumBesselFunction, std::vector>> vvv_d_TableOne){ + std::unordered_map umap_str_d_C4Map = ReadinC4(FileName, NumAtomType, lmax, NumChi, NumBesselFunction); + std::vector>>> vvvv_d_Faln; + for (int indexAtomType = 0; indexAtomType < NumAtomType; indexAtomType++) + { + std::vector>> vvv_d_Faln; + for (int indexl = 0; indexl < lmax+1; indexl++) + { + std::vector> vv_d_Faln; + for (int indexChi = 0; indexChi < NumChi; indexChi++) + { + std::vector v_d_Faln; + for (int indexk = 0; indexk < vvv_d_TableOne[0][0].size(); indexk++) + { + double Faln = 0.0; + for (int indexBesselFunction = 0; indexBesselFunction < NumBesselFunction; indexBesselFunction++) + { + std::string key = std::to_string(indexAtomType) + " " + std::to_string(indexl) + " " + std::to_string(indexChi) + " " + std::to_string(indexBesselFunction); + Faln += umap_str_d_C4Map[key]*vvv_d_TableOne[indexl][indexBesselFunction][indexk]; + } + v_d_Faln.push_back(Faln); + } + vv_d_Faln.push_back(v_d_Faln); + } + vvv_d_Faln.push_back(vv_d_Faln); + } + vvvv_d_Faln.push_back(vvv_d_Faln); + } + return vvvv_d_Faln; +} + +/* OVERLOAD (de)constructors... */ +pseudo_nc::pseudo_nc() +{ +} +pseudo_nc::~pseudo_nc() +{ +} +Atom::Atom() +{ +} +Atom::~Atom() +{ +} +Atom_pseudo::Atom_pseudo() +{ +} +Atom_pseudo::~Atom_pseudo() +{ +} +UnitCell::UnitCell() +{ +} +UnitCell::~UnitCell() +{ +} +Magnetism::Magnetism() +{ +} +Magnetism::~Magnetism() +{ +} + +#ifdef __LCAO +InfoNonlocal::InfoNonlocal() +{ +} +InfoNonlocal::~InfoNonlocal() +{ +} +#endif +/* OVERLOAD printM3 function? */ +/* +void output::printM3(std::ofstream &ofs, const std::string &description, const ModuleBase::Matrix3 &m) +{ +} +*/ + +/************************************************ + * unit test of bessel_basis.cpp + ***********************************************/ + +/* + * - Tested Functions: + * - InitTest: Bessel_Basis::init(start_from_file, ecutwfc, ntype, lmax_in, smooth, sigma, rcut_in, tol_in, ucell, dk, dr) + * - initialize Bessel_Basis class object with following steps: + * - call init_TableOne(...): calculate TableOne (l, ie, ik) = int{dr r^2 jle(r)*jlk(r)} (jle(r) and jlk(r) are Spherical Bessel functions) + * - return, unless set start_from_file = true, will also do the following: + * - call allocate_C4(...): allocate memory for C4 4d-matrix (it, il, in, ie) and set all elements to 1.0 + * - call readin_C4(...): read C4 value from external file and store in C4 4d-matrix + * - call init_Faln(...): calculate F_{aln}(it, il, in, ik) = sum_{ie}{C4(it, il, in, ie)*TableOne(il, ie, ik)} + * - PolynomialInterpolation2Test: Bessel_Basis::Polynomial_Interpolation2(const int &l, const int &ie, const double &gnorm) + * - return (cubic spline) interpolated element value of TableOne 3d-matrix + * - PolynomialInterpolationTest: Bessel_Basis::Polynomial_Interpolation(const int &it, const int &l, const int &ic, const double &gnorm) + * - return (cubic spline) interpolated element value of Faln 4d-matrix + */ + +#define private public +class TestBesselBasis : public ::testing::Test { +protected: + UnitCell ucell; + Bessel_Basis besselBasis; + + void SetUp() override { + ucell.ntype = 1; + ucell.lmax = 0; + ucell.nmax = 1; + + ucell.atoms = new Atom[1]; + ucell.atoms[0].l_nchi = new int[1]; + + ucell.atoms[0].nwl = 0; + ucell.atoms[0].l_nchi[0] = 1; + /* setup_cell manually */ + } + void TearDown() override { + delete[] ucell.atoms->l_nchi; + delete[] ucell.atoms; + } +}; + + TEST_F(TestBesselBasis, InitTest) { + /* parameters required by Bessel_Basis::init() */ + bool b_TestFaln = false; + double d_EnergyCutoff = 4; + int i_Ntype = 1; + int i_Lmax = 0; + bool b_Smooth = false; + double d_SmoothSigma = 0.1; + double d_CutoffRadius = 2; + double d_Tolerance = 0.01; + double d_dk = 0.01; + double d_dr = 2; + /* number of Spherical Bessel functions with given order l, + used to fit one chi */ + int i_Nq = static_cast(sqrt(d_EnergyCutoff)*d_CutoffRadius/M_PI); + /* number of SBF is expected to be 1 */ + besselBasis.init( + b_TestFaln, d_EnergyCutoff, i_Ntype, i_Lmax, b_Smooth, + d_SmoothSigma, d_CutoffRadius, d_Tolerance, + ucell, d_dk, d_dr + ); + EXPECT_EQ(besselBasis.get_ecut_number(), i_Nq); + EXPECT_EQ(besselBasis.get_ecut(), d_EnergyCutoff); + EXPECT_EQ(besselBasis.get_rcut(), d_CutoffRadius); + EXPECT_EQ(besselBasis.get_tolerence(), d_Tolerance); + EXPECT_EQ(besselBasis.get_smooth(), b_Smooth); + EXPECT_EQ(besselBasis.get_sigma(), d_SmoothSigma); +} + +TEST_F(TestBesselBasis, PolynomialInterpolation2Test) { + /* function Bessel_Basis::Polynomial_Interpolation2 is to do + cubic interpolation on TableOne, this matrix has, dimension l*nq*nk */ + /* parameters required by Bessel_Basis::init() */ + bool b_TestFaln = false; + double d_EnergyCutoff = 4; + int i_Ntype = 1; + int i_Lmax = 0; + bool b_Smooth = false; + double d_SmoothSigma = 0.1; + double d_CutoffRadius = 2.0; + double d_Tolerance = 0.01; + double d_dk = 2.0; + double d_dr = 0.01; /* for d_dr will largely impact accurancy of numerical integration, will not give a toy-value */ + + /* number of angular momentum considered should be 1 (lmax = 0)->The first dimension */ + int i_Nq = static_cast(sqrt(d_EnergyCutoff)*d_CutoffRadius/M_PI); + /* number of SBF is expected to be 1->The second dimension */ + int i_kMesh = static_cast(sqrt(d_EnergyCutoff)/d_dk) + 4 + 1; + /* the third dimension is at least to be 6 */ + /* therefore the expected dimension of TableOne is 1*1*6 */ + + besselBasis.init( + b_TestFaln, d_EnergyCutoff, i_Ntype, i_Lmax, b_Smooth, + d_SmoothSigma, d_CutoffRadius, d_Tolerance, + ucell, d_dk, d_dr + ); + /* gnorm for interpolation */ + double d_Gnorm = 1.0; + double d_position = d_Gnorm/d_dk; + int i_position = static_cast(d_position); + assert(i_position < i_kMesh-4); + double d_x0 = d_position - static_cast(i_position); + double d_x1 = 1.0 - d_x0; + double d_x2 = 2.0 - d_x0; + double d_x3 = 3.0 - d_x0; + + std::vector>> vvv_d_TableOne = GenerateTableOne( + b_Smooth, d_SmoothSigma, d_EnergyCutoff, d_CutoffRadius, + i_Lmax, d_dr, d_dk + ); + double d_yExpected = vvv_d_TableOne[0][0][i_position]*d_x1*d_x2*d_x3/6.0+ + vvv_d_TableOne[0][0][i_position]*d_x0*d_x2*d_x3/2.0- + vvv_d_TableOne[0][0][i_position]*d_x1*d_x0*d_x3/2.0+ + vvv_d_TableOne[0][0][i_position]*d_x1*d_x2*d_x0/6.0; + double d_yTested = besselBasis.Polynomial_Interpolation2(0, 0, d_Gnorm); + EXPECT_NEAR(d_yExpected, d_yTested, 0.01); +} +TEST_F(TestBesselBasis, PolynomialInterpolationTest) { + /* function Bessel_Basis::Polynomial_Interpolation is to do + cubic interpolation on Faln, this matrix has, dimension atomtype*l*nchi*nk. + To obtain Faln matrix, it is needed to first get C4 from external file. + The C4 is coefficent of SBF, has dimension atomtype*l*nchi*nq. + Therefore the Faln, equals the contraction of nq between C4 and TableOne. + C4: atomtype*l*nchi*nq + TableOne: l*nq*nk -> Faln: atomtype*l*nchi*nk */ + + /* parameters required by Bessel_Basis::init() */ + bool b_TestFaln = true; + double d_EnergyCutoff = 4; + int i_Ntype = 1; + int i_Lmax = 0; + bool b_Smooth = false; + double d_SmoothSigma = 0.1; + double d_CutoffRadius = 2.0; + double d_Tolerance = 0.01; + double d_dk = 2.0; + double d_dr = 0.01; /* for d_dr will largely impact accurancy of numerical integration, will not give a toy-value */ + + int i_Nq = static_cast(sqrt(d_EnergyCutoff)*d_CutoffRadius/M_PI); + int i_kMesh = static_cast(sqrt(d_EnergyCutoff)/d_dk) + 4 + 1; + /* + manipulate Bessel_Basis::init_Faln function + because for(int it=0; it(d_position); + assert(i_position < i_kMesh-4); + double d_x0 = d_position - static_cast(i_position); + double d_x1 = 1.0 - d_x0; + double d_x2 = 2.0 - d_x0; + double d_x3 = 3.0 - d_x0; + + std::vector>> vvv_d_TableOne = GenerateTableOne( + b_Smooth, d_SmoothSigma, d_EnergyCutoff, d_CutoffRadius, + i_Lmax, d_dr, d_dk + ); + std::vector>>> vvvv_d_Faln = GenerateFaln( + "./support/BesselBasis_UnitTest_C4_AtomType0.html", i_Ntype, i_Lmax, 1, i_Nq, vvv_d_TableOne + ); + double d_yExpected = vvvv_d_Faln[0][0][0][i_position]*d_x1*d_x2*d_x3/6.0+ + vvvv_d_Faln[0][0][0][i_position]*d_x0*d_x2*d_x3/2.0- + vvvv_d_Faln[0][0][0][i_position]*d_x1*d_x0*d_x3/2.0+ + vvvv_d_Faln[0][0][0][i_position]*d_x1*d_x2*d_x0/6.0; + double d_yTested = besselBasis.Polynomial_Interpolation(0, 0, 0, d_Gnorm); + EXPECT_NEAR(d_yExpected, d_yTested, 0.01); +} +#undef private diff --git a/source/module_io/test/support/BesselBasis_UnitTest_C4_AtomType0.html b/source/module_io/test/support/BesselBasis_UnitTest_C4_AtomType0.html new file mode 100644 index 0000000000..05122a725e --- /dev/null +++ b/source/module_io/test/support/BesselBasis_UnitTest_C4_AtomType0.html @@ -0,0 +1,25 @@ + +Bessel Basis Unit Test. Present file will be read as C4 coefficients stored file of atom type 0. +INPUTS tag: +In INPUTS tag what enclosed are EnergyCutoff, CutoffRadius, Nq, Tolerance, respectively. +They should be totally identical with function Bessel_Basis::init_Faln(), or assert() will fail. +C4 tag: +The first line of this tag will be the total chi. The chi acts similarly with contracted Gaussian type +orbitals, which means for one chi, there will be Nq spherical Bessel functions, and for each angular +momentum l, there will be numbers of chi-s to form it. C4 file contains all information of chi-s, +i.e., all over l(s). For example the H, whose numerical orbitals file is named as H_gga_6au_100Ry_2s1p.orb, +which means there are 3 kinds of orbitals, 1s, 2s and 2p. The cutoff radius is 6 au, and the energy cutoff +is 100 Ry. For orbital "1s", there will be a chi_1s, then for "2s" the chi_2s, for "2p" the chi_2p, the +first line is thus the number chi_1s + chi_2s + chi_2p. +The following line has contents in sequence of: +title1 title2 title3 AtomType AngularMomentum ChiIndex C4Vaule1 C4Value2 ... C4ValueNq +Final description: +In unit test case, only one chi, one l is used, and EnergyCutoff, CutoffRadius are selected to have Nq=1. + + +4 2 1 0.01 + + +1 +TEST TEST TEST 0 0 0 1.0 + From 856320dea048c21c9a676a503592dc2a34c39c34 Mon Sep 17 00:00:00 2001 From: dyzheng Date: Fri, 9 Jun 2023 18:06:18 +0800 Subject: [PATCH 12/13] Revert "Refactor: output efermi in pw basis (#2604)" (#2616) This reverts commit a23df8019d1417d0464c3b82d4bfdfc822ccb815. --- source/module_esolver/esolver_ks_lcao.cpp | 7 ++++--- source/module_esolver/esolver_ks_pw.cpp | 11 ++++------- 2 files changed, 8 insertions(+), 10 deletions(-) diff --git a/source/module_esolver/esolver_ks_lcao.cpp b/source/module_esolver/esolver_ks_lcao.cpp index 67ea65062b..02892a3840 100644 --- a/source/module_esolver/esolver_ks_lcao.cpp +++ b/source/module_esolver/esolver_ks_lcao.cpp @@ -787,16 +787,17 @@ void ESolver_KS_LCAO::afterscf(const int istep) if (this->conv_elec) { + // xiaohui add "OUT_LEVEL", 2015-09-16 if (GlobalV::OUT_LEVEL != "m") - { GlobalV::ofs_running << std::setprecision(16); + if (GlobalV::OUT_LEVEL != "m") GlobalV::ofs_running << " EFERMI = " << this->pelec->eferm.ef * ModuleBase::Ry_to_eV << " eV" << std::endl; - } } else { GlobalV::ofs_running << " !! convergence has not been achieved @_@" << std::endl; - std::cout << " !! CONVERGENCE HAS NOT BEEN ACHIEVED !!" << std::endl; + if (GlobalV::OUT_LEVEL == "ie" || GlobalV::OUT_LEVEL == "m") // xiaohui add "m" option, 2015-09-16 + std::cout << " !! CONVERGENCE HAS NOT BEEN ACHIEVED !!" << std::endl; } #ifdef __DEEPKS diff --git a/source/module_esolver/esolver_ks_pw.cpp b/source/module_esolver/esolver_ks_pw.cpp index 03ce818b68..6f46808fcc 100644 --- a/source/module_esolver/esolver_ks_pw.cpp +++ b/source/module_esolver/esolver_ks_pw.cpp @@ -533,16 +533,13 @@ void ESolver_KS_PW::afterscf(const int istep) } if (this->conv_elec) { - if (GlobalV::OUT_LEVEL != "m") - { - GlobalV::ofs_running << std::setprecision(16); - GlobalV::ofs_running << " EFERMI = " << this->pelec->eferm.ef * ModuleBase::Ry_to_eV << " eV" << std::endl; - } + GlobalV::ofs_running << "\n charge density convergence is achieved" << std::endl; + GlobalV::ofs_running << " final etot is " << this->pelec->f_en.etot * ModuleBase::Ry_to_eV << " eV" + << std::endl; } else { - GlobalV::ofs_running << " !! convergence has not been achieved @_@" << std::endl; - std::cout << " !! CONVERGENCE HAS NOT BEEN ACHIEVED !!" << std::endl; + GlobalV::ofs_running << " convergence has NOT been achieved!" << std::endl; } if (GlobalV::OUT_LEVEL != "m") From 8822af8dd8b770030c1dca966bf4fae905a2fdb7 Mon Sep 17 00:00:00 2001 From: jinzx10 Date: Sat, 10 Jun 2023 13:10:07 +0800 Subject: [PATCH 13/13] const qualifiers appended --- source/module_basis/module_nao/radial_collection.cpp | 2 +- source/module_basis/module_nao/radial_collection.h | 10 +++++----- source/module_basis/module_nao/radial_set.cpp | 6 +++--- source/module_basis/module_nao/radial_set.h | 6 +++--- 4 files changed, 12 insertions(+), 12 deletions(-) diff --git a/source/module_basis/module_nao/radial_collection.cpp b/source/module_basis/module_nao/radial_collection.cpp index 15ed392123..7131a90a34 100644 --- a/source/module_basis/module_nao/radial_collection.cpp +++ b/source/module_basis/module_nao/radial_collection.cpp @@ -70,7 +70,7 @@ void RadialCollection::build(const int nfile, const std::string* const file, con } } -void RadialCollection::set_transformer(ModuleBase::SphericalBesselTransformer* sbt, int update) +void RadialCollection::set_transformer(ModuleBase::SphericalBesselTransformer* const sbt, const int update) { for (int itype = 0; itype < ntype_; ++itype) { diff --git a/source/module_basis/module_nao/radial_collection.h b/source/module_basis/module_nao/radial_collection.h index 08f59e11f8..db9e3be062 100644 --- a/source/module_basis/module_nao/radial_collection.h +++ b/source/module_basis/module_nao/radial_collection.h @@ -27,7 +27,7 @@ class RadialCollection * */ //!@{ //! element symbol of a given type - const std::string& symbol(int itype) const { return radset_[itype]->symbol(); } + const std::string& symbol(const int itype) const { return radset_[itype]->symbol(); } //! number of RadialSet objects in the collection int ntype() const { return ntype_; } @@ -39,16 +39,16 @@ class RadialCollection double rcut_max() const; //! number of distinct radial functions of a given type and angular momentum - int nzeta(int itype, int l) const { return radset_[itype]->nzeta(l); } + int nzeta(const int itype, const int l) const { return radset_[itype]->nzeta(l); } //! maximum number of distinct radial functions of a given type among all angular momentum - int nzeta_max(int itype) const { return radset_[itype]->nzeta_max(); } + int nzeta_max(const int itype) const { return radset_[itype]->nzeta_max(); } //! total number of NumericalRadial objects in the collection int nchi() const { return nchi_; } //! number of NumericalRadial objects of a given type - int nchi(int itype) const { return radset_[itype]->nchi(); } + int nchi(const int itype) const { return radset_[itype]->nchi(); } //! get access to the NumericalRadial object with given type, angular momentum and zeta number const NumericalRadial& operator()(const int itype, const int l, const int izeta) const @@ -72,7 +72,7 @@ class RadialCollection //!@{ //! Set a spherical Bessel transformers for all RadialSet objects //! @see RadialSet::set_transformer - void set_transformer(ModuleBase::SphericalBesselTransformer* sbt = nullptr, int update = 0); + void set_transformer(ModuleBase::SphericalBesselTransformer* const sbt = nullptr, const int update = 0); //! Set a common grid for all RadialSet objects //! @see RadialSet::set_grid diff --git a/source/module_basis/module_nao/radial_set.cpp b/source/module_basis/module_nao/radial_set.cpp index 113f8a3607..f194370175 100644 --- a/source/module_basis/module_nao/radial_set.cpp +++ b/source/module_basis/module_nao/radial_set.cpp @@ -105,7 +105,7 @@ RadialSet& RadialSet::operator=(const RadialSet& rhs) return *this; } -int RadialSet::index(int l, int izeta) const +int RadialSet::index(const int l, const int izeta) const { assert(l >= 0 && l <= lmax_); assert(izeta >= 0 && izeta < nzeta_[l]); @@ -142,14 +142,14 @@ void RadialSet::indexing() } } -const NumericalRadial& RadialSet::chi(int l, int izeta) +const NumericalRadial& RadialSet::chi(const int l, const int izeta) { int i = index_map_[l * nzeta_max_ + izeta]; assert(i >= 0 && i < nchi_); return chi_[i]; } -void RadialSet::set_transformer(ModuleBase::SphericalBesselTransformer* sbt, int update) +void RadialSet::set_transformer(ModuleBase::SphericalBesselTransformer* const sbt, const int update) { for (int i = 0; i < nchi_; i++) chi_[i].set_transformer(sbt, update); diff --git a/source/module_basis/module_nao/radial_set.h b/source/module_basis/module_nao/radial_set.h index ae9aae8ec6..7b1b1b1417 100644 --- a/source/module_basis/module_nao/radial_set.h +++ b/source/module_basis/module_nao/radial_set.h @@ -44,7 +44,7 @@ class RadialSet int lmax() const { return lmax_; } double rcut_max() const; - int nzeta(int l) const + int nzeta(const int l) const { assert(l >= 0 && l <= lmax_); return nzeta_[l]; @@ -52,7 +52,7 @@ class RadialSet int nzeta_max() const { return nzeta_max_; } int nchi() const { return nchi_; } - const NumericalRadial& chi(int l, int izeta); + const NumericalRadial& chi(const int l, const int izeta); //!@} /*! @name property setters for all NumericalRadial objects @@ -62,7 +62,7 @@ class RadialSet //!@{ //! Set a spherical Bessel transformers for all NumericalRadial objects //! @see NumericalRadial::set_transformer - void set_transformer(ModuleBase::SphericalBesselTransformer* sbt = nullptr, int update = 0); + void set_transformer(ModuleBase::SphericalBesselTransformer* const sbt = nullptr, const int update = 0); //! Set a common grid for all NumericalRadial objects //! @see NumericalRadial::set_grid