diff --git a/.github/workflows/coverage.yml b/.github/workflows/coverage.yml index 549b5d4757..c552a0b9e4 100644 --- a/.github/workflows/coverage.yml +++ b/.github/workflows/coverage.yml @@ -31,3 +31,4 @@ jobs: if: ${{ ! cancelled() }} with: gcov: true + token: ${{ secrets.CODECOV_TOKEN }} diff --git a/CMakeLists.txt b/CMakeLists.txt index dd237d5e2b..b03faf8fad 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -91,8 +91,9 @@ You can install Git first and reinstall abacus.") RESULT_VARIABLE GIT_COMMIT_DATE_RESULT ) if(GIT_COMMIT_HASH_RESULT EQUAL 0 AND GIT_COMMIT_DATE_RESULT EQUAL 0) - set(COMMIT "${GIT_COMMIT_HASH} (${GIT_COMMIT_DATE})") - add_definitions("-DCOMMIT=\"${COMMIT}\"") + add_definitions(-DCOMMIT_INFO) + file(WRITE "${CMAKE_CURRENT_BINARY_DIR}/commit.h" "#define COMMIT \"${GIT_COMMIT_HASH} (${GIT_COMMIT_DATE})\"\n") + include_directories(${CMAKE_CURRENT_BINARY_DIR}) message(STATUS "Current commit hash: ${GIT_COMMIT_HASH}") message(STATUS "Last commit date: ${GIT_COMMIT_DATE}") else() @@ -320,9 +321,6 @@ endif() # Warning: CMake add support to HIP in version 3.21. This is rather a new version. # Use cmake with AMD-ROCm: https://rocmdocs.amd.com/en/latest/Installation_Guide/Using-CMake-with-AMD-ROCm.html if(USE_ROCM) - if(COMMIT_INFO) - message(FATAL_ERROR "Commit info is not supported on ROCm. Try to recompile with cmake command -DCOMMIT_INFO=OFF.") - endif() if(NOT DEFINED ROCM_PATH ) set (ROCM_PATH "/opt/rocm" CACHE STRING "Default ROCM installation directory." ) endif () @@ -396,10 +394,21 @@ if(MKLROOT) list(APPEND math_libs -lifcore) endif() else() + # In compatibility to builtin FindLAPACK.cmake before v3.5.4 + if(DEFINED LAPACK_DIR) + string(APPEND CMAKE_PREFIX_PATH ";${LAPACK_DIR}") + endif() + if(DEFINED LAPACK_LIBRARY) + set(LAPACK_LIBRARIES ${LAPACK_LIBRARY}) + endif() + if(DEFINED BLAS_DIR) + string(APPEND CMAKE_PREFIX_PATH ";${BLAS_DIR}") + endif() + find_package(FFTW3 REQUIRED) find_package(LAPACK REQUIRED) include_directories(${FFTW3_INCLUDE_DIRS}) - list(APPEND math_libs FFTW3::FFTW3 LAPACK::LAPACK) + list(APPEND math_libs FFTW3::FFTW3 LAPACK::LAPACK BLAS::BLAS) if(ENABLE_LCAO) find_package(ScaLAPACK REQUIRED) @@ -733,4 +742,4 @@ endif() if(ENABLE_RAPIDJSON) target_link_libraries(${ABACUS_BIN_NAME} json_output) -endif() \ No newline at end of file +endif() diff --git a/cmake/FindELPA.cmake b/cmake/FindELPA.cmake index 4105e47592..904be208a5 100644 --- a/cmake/FindELPA.cmake +++ b/cmake/FindELPA.cmake @@ -9,11 +9,27 @@ find_package(PkgConfig) +# Compatible layer towards old manual routines +if(DEFINED ELPA_INCLUDE_DIR) + set(ELPA_INCLUDE_DIRS ${ELPA_INCLUDE_DIR}) +endif() +if(DEFINED ELPA_LIBRARIES) + set(ELPA_LINK_LIBRARIES ${ELPA_LIBRARIES}) +endif() + find_path(ELPA_INCLUDE_DIRS elpa/elpa.h HINTS ${ELPA_DIR} PATH_SUFFIXES "include" "include/elpa" ) +# Fix #3589 +# First if judges if ELPA dir specified +if(ELPA_INCLUDE_DIRS MATCHES "^/usr/include/elpa/.*") + # Second if judges if global visible ELPA header found + if(DEFINED ELPA_DIR OR CMAKE_PREFIX_PATH MATCHES ".*elpa.*") + unset(ELPA_INCLUDE_DIRS) + endif() +endif() if(USE_OPENMP) find_library(ELPA_LINK_LIBRARIES NAMES elpa_openmp elpa @@ -28,6 +44,9 @@ else() ) endif() +# Incompatible with ELPA earlier than 2021.11.001 +# Before ELPA 2021.11.001, its pkg-config file +# is named like "elpa-2021.05.002.pc". if(NOT ELPA_INCLUDE_DIRS AND PKG_CONFIG_FOUND) if(DEFINED ELPA_DIR) string(APPEND CMAKE_PREFIX_PATH ";${ELPA_DIR}") diff --git a/cmake/FindLAPACK.cmake b/cmake/FindLAPACK.cmake deleted file mode 100644 index c240d5facf..0000000000 --- a/cmake/FindLAPACK.cmake +++ /dev/null @@ -1,31 +0,0 @@ -# - Find LAPACK -# Find the native double precision ScaLAPACK headers and libraries. -# -# LAPACK_LIBRARIES - List of libraries when using ScaLAPACK. -# LAPACK_FOUND - True if ScaLAPACK is found. -# - -find_library(LAPACK_LIBRARY - NAMES openblas blas - HINTS ${LAPACK_DIR} - PATH_SUFFIXES "lib" -) - -# Handle the QUIET and REQUIRED arguments and -# set ScaLAPACK_FOUND to TRUE if all variables are non-zero. -include(FindPackageHandleStandardArgs) -find_package_handle_standard_args(LAPACK DEFAULT_MSG LAPACK_LIBRARY) - -# Copy the results to the output variables and target. -if(LAPACK_FOUND) - set(LAPACK_LIBRARIES ${LAPACK_LIBRARY}) - - if(NOT TARGET LAPACK::LAPACK) - add_library(LAPACK::LAPACK UNKNOWN IMPORTED) - set_target_properties(LAPACK::LAPACK PROPERTIES - IMPORTED_LINK_INTERFACE_LANGUAGES "C" - IMPORTED_LOCATION "${LAPACK_LIBRARY}") - endif() -endif() - -mark_as_advanced(LAPACK_LIBRARY) diff --git a/cmake/FindLibxc.cmake b/cmake/FindLibxc.cmake index 4a3c04cba7..2800942ca9 100644 --- a/cmake/FindLibxc.cmake +++ b/cmake/FindLibxc.cmake @@ -3,21 +3,24 @@ include(FindPackageHandleStandardArgs) if(DEFINED Libxc_DIR) string(APPEND CMAKE_PREFIX_PATH ";${Libxc_DIR}") endif() -# Using CMake interface as default. +# Using pkg-config interface as default, to +# avoid linking to wrong global visible Libxc instead of +# specified one. # NO REQUIRED here, otherwhile it would throw error # with no LibXC found. -find_package(Libxc HINTS +find_package(PkgConfig) +if(PKG_CONFIG_FOUND) + pkg_search_module(Libxc IMPORTED_TARGET GLOBAL libxc) + find_package_handle_standard_args(Libxc DEFAULT_MSG Libxc_LINK_LIBRARIES Libxc_INCLUDE_DIRS) +endif() +if(NOT Libxc_FOUND) + find_package(Libxc REQUIRED HINTS ${Libxc_DIR}/share/cmake/Libxc ${Libxc_DIR}/lib/cmake/Libxc ${Libxc_DIR}/lib64/cmake/Libxc ) -if(NOT TARGET Libxc::xc) - find_package(PkgConfig REQUIRED) - pkg_search_module(Libxc REQUIRED IMPORTED_TARGET GLOBAL libxc) - find_package_handle_standard_args(Libxc DEFAULT_MSG Libxc_LINK_LIBRARIES Libxc_INCLUDE_DIRS) endif() - # Copy the results to the output variables and target. # if find_package() above works, Libxc::xc would be present and # below would be skipped. diff --git a/docs/advanced/input_files/input-main.md b/docs/advanced/input_files/input-main.md index 97a90b08dc..517cdefd09 100644 --- a/docs/advanced/input_files/input-main.md +++ b/docs/advanced/input_files/input-main.md @@ -1019,16 +1019,16 @@ Note that `mixing_beta_mag` is not euqal to `mixing_beta` means that $\rho_{up}$ ### mixing_restart -- **Type**: Integer -- **Description**: At `mixing_restart`-th iteration, SCF will restart by using output charge density from perivos iteration as input charge density directly, and start a new mixing. `mixing_restart=0|1` means SCF starts from scratch. +- **Type**: double +- **Description**: If the density difference between input and output `drho` is smaller than `mixing_restart`, SCF will restart at next step which means SCF will restart by using output charge density from perivos iteration as input charge density directly, and start a new mixing. Notice that `mixing_restart` will only take effect once in one SCF. - **Default**: 0 ### mixing_dmr - **Type**: bool -- **Availability**: Only for `mixing_restart>=2` -- **Description**: At `mixing_restart`-th iteration, SCF will start a mixing for real-space density matrix by using the same coefficiences as the mixing of charge density. +- **Availability**: Only for `mixing_restart>=0.0` +- **Description**: At n-th iteration which is calculated by `drho using overload_cast_ = pybind11::detail::overload_cast_impl; -void bind_math_base(py::module& m) +void bind_base_math(py::module& m) { py::module module_base = m.def_submodule("ModuleBase"); @@ -219,4 +220,6 @@ void bind_math_base(py::module& m) static_cast(x_info.ptr), static_cast(w_info.ptr)); }); + py::class_(module_base, "SphericalBesselTransformer") + .def(py::init<>()); } \ No newline at end of file diff --git a/python/pyabacus/src/py_m_nao.cpp b/python/pyabacus/src/py_m_nao.cpp new file mode 100644 index 0000000000..5fd53e5e72 --- /dev/null +++ b/python/pyabacus/src/py_m_nao.cpp @@ -0,0 +1,99 @@ +#include +#include +#include +#include "module_basis/module_nao/radial_collection.h" +#include "module_basis/module_nao/two_center_integrator.h" +#include "module_base/vector3.h" + +namespace py = pybind11; +using namespace pybind11::literals; +template +using overload_cast_ = pybind11::detail::overload_cast_impl; + +void bind_m_nao(py::module& m) +{ + // Create the submodule for Module NAO + py::module m_nao = m.def_submodule("ModuleNAO"); + + // Bind the RadialCollection class + py::class_(m_nao, "RadialCollection") + .def(py::init<>()) + .def("build", [](RadialCollection& self, int nfile, const py::list &file_list, char ftype){ + std::vector files; + files.reserve(nfile); + for (auto file : file_list) + { + files.push_back(file.cast()); + } + self.build(nfile, files.data(), ftype); + }, "nfile"_a, "file_list"_a, "ftype"_a = '\0') + .def("set_transformer", &RadialCollection::set_transformer, "sbt"_a, "update"_a = 0) + .def("set_uniform_grid", &RadialCollection::set_uniform_grid, "for_r_space"_a, "ngrid"_a, "cutoff"_a, "mode"_a = 'i', "enable_fft"_a = false) + .def("set_grid", [](RadialCollection& self, const bool for_r_space, const int ngrid, py::array_t grid, const char mode = 'i'){ + py::buffer_info grid_info = grid.request(); + if (grid_info.ndim != 1) + { + throw std::runtime_error("Input array must be 1-dimensional"); + } + self.set_grid(for_r_space, ngrid, static_cast(grid_info.ptr), mode); + }, "for_r_space"_a, "ngrid"_a, "grid"_a, "mode"_a = 'i') + // Getters + .def("symbol", &RadialCollection::symbol, "itype"_a) + .def_property_readonly("ntype", &RadialCollection::ntype) + .def("lmax", overload_cast_()(&RadialCollection::lmax, py::const_), "itype"_a) + .def_property_readonly("lmax", overload_cast_<>()(&RadialCollection::lmax, py::const_)) + .def("rcut_max", overload_cast_()(&RadialCollection::rcut_max, py::const_), "itype"_a) + .def_property_readonly("rcut_max", overload_cast_<>()(&RadialCollection::rcut_max, py::const_)) + .def("nzeta", &RadialCollection::nzeta, "itype"_a, "l"_a) + .def("nzeta_max", overload_cast_()(&RadialCollection::nzeta_max, py::const_), "itype"_a) + .def_property_readonly("nzeta_max", overload_cast_<>()(&RadialCollection::nzeta_max, py::const_)) + .def("nchi", overload_cast_()(&RadialCollection::nchi, py::const_), "itype"_a) + .def_property_readonly("nchi", overload_cast_<>()(&RadialCollection::nchi, py::const_)); + //Bind the TwoCenterIntegrator class + py::class_(m_nao, "TwoCenterIntegrator") + .def(py::init<>()) + .def("tabulate", &TwoCenterIntegrator::tabulate, "bra"_a, "ket"_a, "op"_a, "nr"_a, "cutoff"_a) + .def("calculate", [](TwoCenterIntegrator& self, const int itype1, const int l1, const int izeta1, const int m1, const int itype2, const int l2, const int izeta2, const int m2, py::array_t pvR, bool cal_grad = false) + { + py::buffer_info pvR_info = pvR.request(); + if (pvR_info.size != 3) + { + throw std::runtime_error("Radial part must have 3 elements"); + } + double* cvR = static_cast(pvR_info.ptr); + ModuleBase::Vector3 vR(cvR[0], cvR[1], cvR[2]); + double out[1] = {0.0}; + double* grad_out = nullptr; + if (cal_grad) + { + grad_out = new double[3]; + } + self.calculate(itype1, l1, izeta1, m1, itype2, l2, izeta2, m2, vR, out, grad_out); + py::array_t out_array({1}, out); + if (cal_grad) + { + py::array_t grad_out_array({3}, grad_out); + return py::make_tuple(out_array, grad_out_array); + } + else + { + py::array_t grad_out_array({0}); + return py::make_tuple(out_array, grad_out_array); + } + + }, "itype1"_a, "l1"_a, "izeta1"_a, "m1"_a, "itype2"_a, "l2"_a, "izeta2"_a, "m2"_a, "pvR"_a, "cal_grad"_a = false) + .def("snap", [](TwoCenterIntegrator& self, const int itype1, const int l1, const int izeta1, const int m1, const int itype2, py::array_t pvR, const bool deriv){ + py::buffer_info pvR_info = pvR.request(); + if (pvR_info.size != 3) + { + throw std::runtime_error("Radial part must have 3 elements"); + } + double* cvR = static_cast(pvR_info.ptr); + ModuleBase::Vector3 vR(cvR[0], cvR[1], cvR[2]); + // TODO: check deriv & out memory allocation + std::vector> out; + self.snap(itype1, l1, izeta1, m1, itype2, vR, deriv, out); + return out; + }, "itype1"_a, "l1"_a, "izeta1"_a, "m1"_a, "itype2"_a, "pvR"_a, "deriv"_a = false); + +} \ No newline at end of file diff --git a/python/pyabacus/src/pyabacus/__init__.py b/python/pyabacus/src/pyabacus/__init__.py index 94d8c0d5b8..f20598fbea 100644 --- a/python/pyabacus/src/pyabacus/__init__.py +++ b/python/pyabacus/src/pyabacus/__init__.py @@ -1,4 +1,4 @@ from __future__ import annotations # from ._core import __doc__, __version__, NumericalRadial, ModuleBase -from ._core import ModuleBase -__all__ = ["ModuleBase"] \ No newline at end of file +from ._core import ModuleBase, ModuleNAO +__all__ = ["ModuleBase", "ModuleNAO"] \ No newline at end of file diff --git a/python/pyabacus/tests/indexmap.py b/python/pyabacus/tests/indexmap.py new file mode 100644 index 0000000000..f9917919d8 --- /dev/null +++ b/python/pyabacus/tests/indexmap.py @@ -0,0 +1,107 @@ +def _index_map(ntype, natom, lmax, nzeta=None): + ''' + Bijective map between the composite index (itype, iatom, l, zeta, m) + and linearized orbital index mu. + + An atomic orbital is labeled by its type & atomic index, its angular + momentum quantum numbers l & m, and possibly a zeta number. Suppose + there's a total of N orbitals, each orbital can also be assigned a + unique index mu \in [0, N-1]. + + This function returns a pair of bijective maps between the composite + index (itype, iatom, l, zeta, m) and the linearized orbital index mu. + The composite index is linearized in C-style (lexicographic order). + + Parameters + ---------- + ntype : int + Number of atom types. + natom : list of int + Number of atoms for each type. + lmax : list of int + lmax[i] specifies the maximum angular momentum of type i. + nzeta : list of list of int + nzeta[i][l] specifies the number of zeta orbitals of the + angular momentum l of type i. + If None, nzeta is assumed to be 1 for all. + + Returns + ------- + comp2mu : dict + A dict that maps a composite index (itype, iatom, l, zeta, m) + to its linearized index. + mu2comp : dict + A dict that maps a linearized index to its composite index. + + Notes + ----- + The linearized index arranges m in accordance with ABACUS: + + 0, 1, -1, 2, -2, 3, -3, ..., l, -l + + ''' + if nzeta is None: + nzeta = [[1]*(lmax[itype]+1) for itype in range(ntype)] + + assert len(natom) == ntype + assert len(lmax) == ntype + assert lmax == [len(nzeta[itype])-1 for itype in range(ntype)] + + comp2mu = {} + mu2comp = {} + mu = 0 + for itype in range(ntype): + for iatom in range(natom[itype]): + for l in range(lmax[itype]+1): + for zeta in range(nzeta[itype][l]): + ''' + In ABACUS, real spherical harmonics Ylm arranges its m + in the following order: + + 0, 1, -1, 2, -2, 3, -3, ..., l, -l + + (see module_base/ylm.cpp and module_base/math_ylmreal.cpp + for details) + + ''' + for mm in range(0, 2*l+1): + m = -mm // 2 if mm % 2 == 0 else (mm + 1) // 2 + comp2mu[(itype, iatom, l, zeta, m)] = mu + mu2comp[mu] = (itype, iatom, l, zeta, m) + mu += 1 + + return comp2mu, mu2comp + + +############################################################ +# Test +############################################################ +import unittest + +class _TestIndexMap(unittest.TestCase): + + def test_index_map(self): + ntype = 3 + natom = [2, 1, 3] + lmax = [1, 2, 4] + nzeta = [[2,3], [1,0,1], [1, 2, 2, 1, 3]] + comp2mu, mu2comp = _index_map(ntype, natom, lmax, nzeta) + + # check the total number of orbitals + nao = sum(sum( (2*l+1) * nzeta[itype][l] for l in range(lmax[itype]+1) ) \ + * natom[itype] for itype in range(ntype)) + self.assertEqual( len(mu2comp.items()), nao ) + + # check bijectivity + for mu in range(nao): + self.assertEqual( comp2mu[mu2comp[mu]], mu ) + + # check the first and the last + self.assertEqual( mu2comp[0], (0, 0, 0, 0, 0) ) + self.assertEqual( mu2comp[nao-1], \ + (ntype-1, natom[-1]-1, lmax[-1], nzeta[-1][-1]-1, -lmax[-1]) ) + + +if __name__ == '__main__': + unittest.main() + diff --git a/python/pyabacus/tests/test_base_math.py b/python/pyabacus/tests/test_base_math.py index 2dcb55d8c3..671d10edc8 100644 --- a/python/pyabacus/tests/test_base_math.py +++ b/python/pyabacus/tests/test_base_math.py @@ -1,17 +1,20 @@ from __future__ import annotations import pytest -import pyabacus as m +from pyabacus import ModuleBase as m import numpy as np def test_sphbes(): - s = m.ModuleBase.Sphbes() + s = m.Sphbes() # test for sphbesj assert s.sphbesj(1, 0.0) == 0.0 assert s.sphbesj(0, 0.0) == 1.0 +def test_sbt(): + sbt = m.SphericalBesselTransformer() + @pytest.fixture def simpson_setup(): n = 1000 @@ -22,7 +25,7 @@ def simpson_setup(): def test_simpson(simpson_setup): n, func, dx= simpson_setup - s = m.ModuleBase.Integral() + s = m.Integral() assert s.simpson(n, func, dx) == pytest.approx(0, abs=1e-10) diff --git a/python/pyabacus/tests/test_m_nao.py b/python/pyabacus/tests/test_m_nao.py new file mode 100644 index 0000000000..0890a06635 --- /dev/null +++ b/python/pyabacus/tests/test_m_nao.py @@ -0,0 +1,99 @@ +from __future__ import annotations + +# import pytest +from pyabacus import ModuleNAO as nao +from pyabacus import ModuleBase as base +import numpy as np + +import copy +from indexmap import _index_map +from wigner_real import R, wigner_d_real + +import matplotlib.pyplot as plt + +orb_dir = '../../../tests/PP_ORB/' +file_list = ['H_gga_8au_100Ry_2s1p.orb', 'O_gga_10au_100Ry_2s2p1d.orb'] + +file_list = [orb_dir + orbfile for orbfile in file_list ] +nfile = len(file_list) + +orb = nao.RadialCollection() +orb.build(nfile, file_list, 'o') +sbt = base.SphericalBesselTransformer() +orb.set_transformer(sbt) + +rmax = orb.rcut_max * 2.0 +dr = 0.01 +nr = int(rmax / dr) + 1 + +orb.set_uniform_grid(True, nr, rmax, 'i', True) +S_intor = nao.TwoCenterIntegrator() +S_intor.tabulate(orb, orb, 'S', nr, rmax) +vR = np.array([1.0,0,0]).astype(np.float64) +#out = S_intor.calculate(1, 1, 0, 1, 0, 1, 0, 1, vR)[0] +#print(out) + +comp2mu, mu2comp = _index_map(2, [2, 1], [1, 2], [[2,1], [2,2,1]]) +#print(comp2mu, mu2comp) + +nao = len(comp2mu.items()) + +coord = [ + [ + np.array([-12.081531451316582,16.463368531712373,10.304287878967891]), + np.array([ -12.056180479123837,19.25408045699522,10.010554611214044]) + ], + [ + np.array([-13.1930046246741,17.91132430713516,10.440289103003526]) + ] + ] + + +S = np.zeros((nao, nao)) +for mu in range(nao): + itype1, iatom1, l1, zeta1, m1 = mu2comp[mu] + for nu in range(nao): + itype2, iatom2, l2, zeta2, m2 = mu2comp[nu] + vR = coord[itype2][iatom2] - coord[itype1][iatom1] + S[mu, nu] = S_intor.calculate(itype1, l1, zeta1, m1, itype2, l2, zeta2, m2, vR)[0] + + +S_prerot = np.copy(S) +coord_pre = copy.deepcopy(coord) + +val_pre, vec_pre = np.linalg.eigh(S_prerot) + +# random rotation +alpha = np.random.rand() * 2 * np.pi +beta = np.random.rand() * np.pi +gamma = np.random.rand() * 2 * np.pi +Rmat = R(alpha, beta, gamma) + +for coord_type in coord: + for i in range(len(coord_type)): + coord_type[i] = Rmat @ coord_type[i] + +D = np.zeros((nao, nao)) +mu = 0 +while mu < nao: + itype, iatom, l, zeta, m = mu2comp[mu] + D[mu:mu+2*l+1, mu:mu+2*l+1] = wigner_d_real(l, alpha, beta, gamma) + mu += 2*l+1 + +S_e3 = D @ S_prerot @ D.T # overlap matrix after rotation (by angular momentum theory) + +# numerical evaluation of S_e3 by two-center integrals +S = np.zeros((nao, nao)) +for mu in range(nao): + itype1, iatom1, l1, zeta1, m1 = mu2comp[mu] + for nu in range(nao): + itype2, iatom2, l2, zeta2, m2 = mu2comp[nu] + vR = coord[itype2][iatom2] - coord[itype1][iatom1] + S[mu, nu] = S_intor.calculate(itype1, l1, zeta1, m1, itype2, l2, zeta2, m2, vR)[0] + +val, vec = np.linalg.eigh(S) + +#print('S eigval diff: ', np.linalg.norm(val-val_pre, np.inf)) +print('norm(S_e3 - S_numer) = ', np.linalg.norm(S_e3 - S, np.inf)) + + diff --git a/python/pyabacus/tests/wigner_real.py b/python/pyabacus/tests/wigner_real.py new file mode 100644 index 0000000000..656a29fba0 --- /dev/null +++ b/python/pyabacus/tests/wigner_real.py @@ -0,0 +1,200 @@ +import numpy as np +from sympy.physics.wigner import wigner_d +from sympy import matrix2numpy, Ynm, Matrix + +# rotation matrices (counter-clockwise) +def Rz(angle): + return np.array([[np.cos(angle), -np.sin(angle), 0], + [np.sin(angle), np.cos(angle), 0], + [0, 0, 1]]) + + +def Ry(angle): + return np.array([[np.cos(angle), 0, np.sin(angle)], + [0, 1, 0], + [-np.sin(angle), 0, np.cos(angle)]]) + + +def R(alpha, beta, gamma): + return Rz(alpha) @ Ry(beta) @ Rz(gamma) + + +def cart2sph(x, y, z): + ''' + Cartesian to spherical coordinates. + + ''' + r = np.sqrt(x**2 + y**2 + z**2) + polar = np.arccos(z/r) + azimuth = np.arctan2(y, x) + return r, polar, azimuth + + +def sph2cart(r, polar, azimuth): + ''' + Spherical to Cartesian coordinates. + + ''' + x = r * np.sin(polar) * np.cos(azimuth) + y = r * np.sin(polar) * np.sin(azimuth) + z = r * np.cos(polar) + return x, y, z + + +def Ylm_real(l, m, polar, azimuth): + ''' + Real spherical harmonics (with sign convention in accordance with ABACUS) + + ''' + if m > 0: + return np.sqrt(2) * complex(Ynm(l, m, polar, azimuth)).real + elif m < 0: + return np.sqrt(2) * complex(Ynm(l, -m, polar, azimuth)).imag + else: # m == 0 + return float(Ynm(l, 0, polar, azimuth)) + + +def Yl_real(l, polar, azimuth, m_order='abacus'): + ''' + An array of real spherical harmonics for a given l. + + Note + ---- + ABACUS arranges m in the following order: + + 0, 1, -1, 2, -2, ..., l, -l + + + ''' + if m_order == 'abacus': + return np.array([Ylm_real(l, (m+1)//2 * (-1)**(m+1), polar, azimuth) \ + for m in range(0, 2*l+1)]) + elif m_order == 'descend': + return np.array([Ylm_real(l, m, polar, azimuth) for m in range(l,-l-1,-1)]) + elif m_order == 'ascend': + return np.array([Ylm_real(l, m, polar, azimuth) for m in range(-l, l+1)]) + else: + raise ValueError('m_order must be "abacus", "descend" or "ascend"') + + +def ToReal(l): + ''' + Standard-to-real transform for spherical harmonics: + + |lm_real> = \sum_{m'} |lm'> U[m',m] + + where m is arranged in descending order (in order to use with sympy's wigner_d) + + ''' + U = np.zeros((2*l+1, 2*l+1), dtype=complex) + U[l, l] = 1.0 + for m in range(1, l+1): + U[l-m, l-m] = 1.0 / np.sqrt(2) + U[l-m, l+m] = -1j / np.sqrt(2) + U[l+m, l-m] = (-1)**m / np.sqrt(2) + U[l+m, l+m] = 1j * (-1)**m / np.sqrt(2) + + return U + + +def wigner_d_real(l, alpha, beta, gamma): + ''' + Wigner D-matrix for real spherical harmonics. + + This function returns the Wigner D-matrix for real spherical harmonics where + m is arranged in accordance with ABACUS, and the rotations are counter-clockwise: + + D_real[m',m] = + + Notes + ----- + sympy's wigner_d generates the Wigner D-matrix for standard spherical + harmonics where m is arranged in descending order and the rotations + are clockwise. See sympy's source code and documentation for details. + + ''' + D = matrix2numpy(wigner_d(l, -alpha, -beta, -gamma), dtype=complex) + + # transform to real spherical harmonics + U = ToReal(l) + D_real = np.real(U.T.conj() @ D @ U) + + # rearrange m to the ABACUS order from the descending order + idx = np.zeros(2*l+1, dtype=int) + idx[::2] = np.arange(l, 2*l+1) + idx[1::2] = np.arange(l-1, -1, -1) + # idx = [l, l-1, l+1, l-2, l+2, ..., 0, 2*l] + + D_real = D_real[:, idx] + D_real = D_real[idx, :] + + return D_real + + +######################################################################## +# Tests +######################################################################## + +import unittest + +class _TestWignerReal(unittest.TestCase): + + def test_yreal(self): + # random vector + v = np.random.randn(3) + r, polar, azimuth = cart2sph(*v) + + C = np.sqrt(0.75/np.pi) + YReal_ref = [C*v[2]/r, -C*v[0]/r, -C*v[1]/r] # z/r, -x/r, -y/r + YReal_1 = [Ylm_real(1, m, polar, azimuth) for m in [0, 1, -1]] + YReal_2 = Yl_real(1, polar, azimuth) + + self.assertTrue(np.allclose(YReal_ref, YReal_1, rtol=1e-12, atol=1e-12)) + self.assertTrue(np.allclose(YReal_ref, YReal_2, rtol=1e-12, atol=1e-12)) + + + def test_toreal(self): + l = 3 + + # random vector + v = np.random.randn(3) + r, polar, azimuth = cart2sph(*v) + + Yl = np.array([complex(Ynm(l, m, polar, azimuth)) for m in range(l, -l-1, -1)]) + Yl_real_ref = Yl_real(l, polar, azimuth, m_order='descend') + U = ToReal(l) + + # unitarity + self.assertTrue(np.allclose(U.T.conj() @ U, np.eye(2*l+1), rtol=1e-12, atol=1e-12)) + + # Yl_real_ref == Yl @ U + self.assertTrue(np.allclose(Yl @ U, Yl_real_ref, rtol=1e-12, atol=1e-12)) + + + def test_wigner_d_real(self): + l = 3 + + # Yl_real along some random vector v + v = np.random.randn(3) + _, polar, azimuth = cart2sph(*v) + Yl_real_v = Yl_real(l, polar, azimuth) # ABACUS order + + # Euler angles + alpha = np.random.rand() * 2 * np.pi + beta = np.random.rand() * np.pi + gamma = np.random.rand() * 2 * np.pi + + # Yl_real along the rotated vector u = R @ v + u = R(alpha, beta, gamma) @ v + _, polar, azimuth = cart2sph(*u) + Yl_real_u = Yl_real(l, polar, azimuth) # ABACUS order + + # Wigner D-matrix for real spherical harmonics with m in ABACUS order + D_real = wigner_d_real(l, alpha, beta, gamma) + + # Yl_real_v == Yl_real_u @ D_real + self.assertTrue( np.allclose(Yl_real_v, Yl_real_u @ D_real) ) + + +if __name__ == '__main__': + unittest.main() diff --git a/source/Makefile.Objects b/source/Makefile.Objects index 47172623ee..a71b0e8de7 100644 --- a/source/Makefile.Objects +++ b/source/Makefile.Objects @@ -448,12 +448,11 @@ OBJS_IO=input.o\ output_rho.o\ output_potential.o\ output_mat_sparse.o\ - output_radial.o\ para_json.o\ abacusjson.o\ general_info.o\ init_info.o\ - + readin_info.o\ OBJS_IO_LCAO=cal_r_overlap_R.o\ write_orb_info.o\ diff --git a/source/driver_run.cpp b/source/driver_run.cpp index ffa5bc7488..5222be6aa5 100644 --- a/source/driver_run.cpp +++ b/source/driver_run.cpp @@ -5,6 +5,7 @@ #include "module_io/print_info.h" #include "module_io/winput.h" #include "module_md/run_md.h" +#include "module_io/para_json.h" /** * @brief This is the driver function which defines the workflow of ABACUS calculations. @@ -44,9 +45,14 @@ void Driver::driver_run() ModuleBase::QUIT(); } - // 4. Initialize Esolver + // 4. Initialize Esolver,and fill json-structure p_esolver->Init(INPUT, GlobalC::ucell); + +#ifdef __RAPIDJSON + Json::gen_stru_wrapper(&GlobalC::ucell); +#endif + //------------------------------------------------------------ // This part onward needs to be refactored. //---------------------------MD/Relax------------------------- diff --git a/source/module_base/atom_in.h b/source/module_base/atom_in.h index 4ee66cad28..5458a28fd7 100644 --- a/source/module_base/atom_in.h +++ b/source/module_base/atom_in.h @@ -108,6 +108,57 @@ class atom_in {"Mt", 7}, {"Ds", 7}, {"Rg", 7}, {"Cn", 7}, {"Nh", 7}, {"Fl", 7}, {"Mc", 7}, {"Lv", 7}, {"Ts", 7}, {"Og", 7} }; + /// @brief ground state electron configuration, sequence of orbitals in key is in accord with the sequence of n then l + /// @note 1s2s2p3s3p 4s 3d4p5s 4d5p6s 4f5d6p7s 5f6d7p, from NIST periodic table, + /// @details see: https://www.nist.gov/system/files/documents/2019/12/10/nist_periodictable_july2019_crop.pdf + std::map> groundstate_electronconfiguration + = { + // 1s + {"H", {1}}, {"He", {2}}, // 1st period + // 1s 2s + {"Li", {2, 1}}, {"Be", {2, 2}}, // 2nd period + // 1s 2s 2p + {"B", {2, 2, 1}}, {"C", {2, 2, 2}}, {"N", {2, 2, 3}}, {"O", {2, 2, 4}}, {"F", {2, 2, 5}}, {"Ne", {2, 2, 6}}, // 2nd period + // 1s 2s 2p 3s 1s 2s 2p 3s + {"Na", {2, 2, 6, 1}}, {"Mg", {2, 2, 6, 2}}, // 3rd period + // 1s 2s 2p 3s 3p 1s 2s 2p 3s 3p 1s 2s 2p 3s 3p + {"Al", {2, 2, 6, 2, 1}}, {"Si", {2, 2, 6, 2, 2}}, {"P", {2, 2, 6, 2, 3}}, // 3rd period + {"S", {2, 2, 6, 2, 4}}, {"Cl", {2, 2, 6, 2, 5}}, {"Ar", {2, 2, 6, 2, 6}}, // 3rd period + // 1s 2s 2p 3s 3p 3d 4s 1s 2s 2p 3s 3p 3d 4s 1s 2s 2p 3s 3p 3d 4s + {"K", {2, 2, 6, 2, 6, 0, 1}}, {"Ca", {2, 2, 6, 2, 6, 0, 2}}, {"Sc", {2, 2, 6, 2, 6, 1, 2}}, // 4th period + {"Ti", {2, 2, 6, 2, 6, 2, 2}}, {"V", {2, 2, 6, 2, 6, 3, 2}}, {"Cr", {2, 2, 6, 2, 6, 4, 2}}, // 4th period + {"Mn", {2, 2, 6, 2, 6, 5, 2}}, {"Fe", {2, 2, 6, 2, 6, 6, 2}}, {"Co", {2, 2, 6, 2, 6, 7, 2}}, // 4th period + {"Ni", {2, 2, 6, 2, 6, 8, 2}}, {"Cu", {2, 2, 6, 2, 6, 10, 1}}, {"Zn", {2, 2, 6, 2, 6, 10, 2}}, // 4th period + // 1s 2s 2p 3s 3p 3d 4s 4p 1s 2s 2p 3s 3p 3d 4s 4p 1s 2s 2p 3s 3p 3d 4s 4p + {"Ga", {2, 2, 6, 2, 6, 10, 2, 1}}, {"Ge", {2, 2, 6, 2, 6, 10, 2, 2}}, {"As", {2, 2, 6, 2, 6, 10, 2, 3}}, // 4th period + {"Se", {2, 2, 6, 2, 6, 10, 2, 4}}, {"Br", {2, 2, 6, 2, 6, 10, 2, 5}}, {"Kr", {2, 2, 6, 2, 6, 10, 2, 6}}, // 4th period + // 1s 2s 2p 3s 3p 3d 4s 4p 4d 4f 5s 1s 2s 2p 3s 3p 3d 4s 4p 4d 4f 5s 1s 2s 2p 3s 3p 3d 4s 4p 4d 4f 5s + {"Rb", {2, 2, 6, 2, 6, 10, 2, 6, 0, 0, 1}}, {"Sr", {2, 2, 6, 2, 6, 10, 2, 6, 0, 0, 2}}, {"Y", {2, 2, 6, 2, 6, 10, 2, 6, 1, 0, 2}}, // 5th period + {"Zr", {2, 2, 6, 2, 6, 10, 2, 6, 2, 0, 2}}, {"Nb", {2, 2, 6, 2, 6, 10, 2, 6, 4, 0, 1}}, {"Mo", {2, 2, 6, 2, 6, 10, 2, 6, 5, 0, 1}}, // 5th period + {"Tc", {2, 2, 6, 2, 6, 10, 2, 6, 5, 0, 2}}, {"Ru", {2, 2, 6, 2, 6, 10, 2, 6, 7, 0, 1}}, {"Rh", {2, 2, 6, 2, 6, 10, 2, 6, 8, 0, 1}}, // 5th period + {"Pd", {2, 2, 6, 2, 6, 10, 2, 6, 10, 0, 0}}, {"Ag", {2, 2, 6, 2, 6, 10, 2, 6, 10, 0, 1}}, {"Cd", {2, 2, 6, 2, 6, 10, 2, 6, 10, 0, 2}}, // 5th period + // 1s 2s 2p 3s 3p 3d 4s 4p 4d 4f 5s 5p 1s 2s 2p 3s 3p 3d 4s 4p 4d 4f 5s 5p 1s 2s 2p 3s 3p 3d 4s 4p 4d 4f 5s 5p + {"In", {2, 2, 6, 2, 6, 10, 2, 6, 10, 0, 2, 1}}, {"Sn", {2, 2, 6, 2, 6, 10, 2, 6, 10, 0, 2, 2}}, {"Sb", {2, 2, 6, 2, 6, 10, 2, 6, 10, 0, 2, 3}}, // 6th period + {"Te", {2, 2, 6, 2, 6, 10, 2, 6, 10, 0, 2, 4}}, {"I", {2, 2, 6, 2, 6, 10, 2, 6, 10, 0, 2, 5}}, {"Xe", {2, 2, 6, 2, 6, 10, 2, 6, 10, 0, 2, 6}}, // 6th period + // 1s 2s 2p 3s 3p 3d 4s 4p 4d 4f 5s 5p 5d 5f 5g 6s 1s 2s 2p 3s 3p 3d 4s 4p 4d 4f 5s 5p 5d 5f 5g 6s + {"Cs", {2, 2, 6, 2, 6, 10, 2, 6, 10, 0, 2, 6, 0, 0, 0, 1}}, {"Ba", {2, 2, 6, 2, 6, 10, 2, 6, 10, 0, 2, 6, 0, 0, 0, 2}}, // 6th period + {"La", {2, 2, 6, 2, 6, 10, 2, 6, 10, 0, 2, 6, 1, 0, 0, 2}}, {"Ce", {2, 2, 6, 2, 6, 10, 2, 6, 10, 1, 2, 6, 1, 0, 0, 2}}, // 6th period + {"Pr", {2, 2, 6, 2, 6, 10, 2, 6, 10, 3, 2, 6, 0, 0, 0, 2}}, {"Nd", {2, 2, 6, 2, 6, 10, 2, 6, 10, 4, 2, 6, 0, 0, 0, 2}}, // 6th period + {"Pm", {2, 2, 6, 2, 6, 10, 2, 6, 10, 5, 2, 6, 0, 0, 0, 2}}, {"Sm", {2, 2, 6, 2, 6, 10, 2, 6, 10, 6, 2, 6, 0, 0, 0, 2}}, // 6th period + {"Eu", {2, 2, 6, 2, 6, 10, 2, 6, 10, 7, 2, 6, 0, 0, 0, 2}}, {"Gd", {2, 2, 6, 2, 6, 10, 2, 6, 10, 7, 2, 6, 1, 0, 0, 2}}, // 6th period + {"Tb", {2, 2, 6, 2, 6, 10, 2, 6, 10, 9, 2, 6, 0, 0, 0, 2}}, {"Dy", {2, 2, 6, 2, 6, 10, 2, 6, 10, 10, 2, 6, 0, 0, 0, 2}}, // 6th period + {"Ho", {2, 2, 6, 2, 6, 10, 2, 6, 10, 11, 2, 6, 0, 0, 0, 2}}, {"Er", {2, 2, 6, 2, 6, 10, 2, 6, 10, 12, 2, 6, 0, 0, 0, 2}}, // 6th period + {"Tm", {2, 2, 6, 2, 6, 10, 2, 6, 10, 13, 2, 6, 0, 0, 0, 2}}, {"Yb", {2, 2, 6, 2, 6, 10, 2, 6, 10, 14, 2, 6, 0, 0, 0, 2}}, // 6th period + {"Lu", {2, 2, 6, 2, 6, 10, 2, 6, 10, 14, 2, 6, 1, 0, 0, 2}}, {"Hf", {2, 2, 6, 2, 6, 10, 2, 6, 10, 14, 2, 6, 2, 0, 0, 2}}, // 6th period + {"Ta", {2, 2, 6, 2, 6, 10, 2, 6, 10, 14, 2, 6, 3, 0, 0, 2}}, {"W", {2, 2, 6, 2, 6, 10, 2, 6, 10, 14, 2, 6, 4, 0, 0, 2}}, // 6th period + {"Re", {2, 2, 6, 2, 6, 10, 2, 6, 10, 14, 2, 6, 5, 0, 0, 2}}, {"Os", {2, 2, 6, 2, 6, 10, 2, 6, 10, 14, 2, 6, 6, 0, 0, 2}}, // 6th period + {"Ir", {2, 2, 6, 2, 6, 10, 2, 6, 10, 14, 2, 6, 7, 0, 0, 2}}, {"Pt", {2, 2, 6, 2, 6, 10, 2, 6, 10, 14, 2, 6, 9, 0, 0, 1}}, // 6th period + {"Au", {2, 2, 6, 2, 6, 10, 2, 6, 10, 14, 2, 6, 10, 0, 0, 1}}, {"Hg", {2, 2, 6, 2, 6, 10, 2, 6, 10, 14, 2, 6, 10, 0, 0, 2}}, // 6th period + // 1s 2s 2p 3s 3p 3d 4s 4p 4d 4f 5s 5p 5d 5f 5g 6s 6p 1s 2s 2p 3s 3p 3d 4s 4p 4d 4f 5s 5p 5d 5f 5g 6s 6p + {"Tl", {2, 2, 6, 2, 6, 10, 2, 6, 10, 14, 2, 6, 10, 0, 0, 2, 1}}, {"Pb", {2, 2, 6, 2, 6, 10, 2, 6, 10, 14, 2, 6, 10, 0, 0, 2, 2}}, // 6th period + {"Bi", {2, 2, 6, 2, 6, 10, 2, 6, 10, 14, 2, 6, 10, 0, 0, 2, 3}}, {"Po", {2, 2, 6, 2, 6, 10, 2, 6, 10, 14, 2, 6, 10, 0, 0, 2, 4}}, // 6th period + {"At", {2, 2, 6, 2, 6, 10, 2, 6, 10, 14, 2, 6, 10, 0, 0, 2, 5}}, {"Rn", {2, 2, 6, 2, 6, 10, 2, 6, 10, 14, 2, 6, 10, 0, 0, 2, 6}} // 6th period + }; }; #endif \ No newline at end of file diff --git a/source/module_base/global_variable.cpp b/source/module_base/global_variable.cpp index 9a8fe56748..c1f769e93a 100644 --- a/source/module_base/global_variable.cpp +++ b/source/module_base/global_variable.cpp @@ -248,7 +248,7 @@ std::string of_kernel_file = "WTkernel.txt"; std::string MIXING_MODE = "broyden"; double MIXING_BETA = 0.7; int MIXING_NDIM = 8; -int MIXING_RESTART = 0; +double MIXING_RESTART = 0.0; double MIXING_GG0 = 1.00; double MIXING_BETA_MAG = 1.6; double MIXING_GG0_MAG = 1.00; diff --git a/source/module_base/global_variable.h b/source/module_base/global_variable.h index fddcb4f734..2df0c6f57a 100644 --- a/source/module_base/global_variable.h +++ b/source/module_base/global_variable.h @@ -277,7 +277,7 @@ extern std::string of_kernel_file; // The name of WT kernel file. extern std::string MIXING_MODE; extern double MIXING_BETA; extern int MIXING_NDIM; -extern int MIXING_RESTART; +extern double MIXING_RESTART; extern double MIXING_GG0; extern bool MIXING_TAU; extern double MIXING_BETA_MAG; diff --git a/source/module_base/parallel_global.cpp b/source/module_base/parallel_global.cpp index e67e563455..8b8e6a7bab 100644 --- a/source/module_base/parallel_global.cpp +++ b/source/module_base/parallel_global.cpp @@ -223,7 +223,8 @@ void Parallel_Global::read_mpi_parameters(int argc,char **argv) #else const char* version = "unknown"; #endif -#ifdef COMMIT +#ifdef COMMIT_INFO +#include "commit.h" const char* commit = COMMIT; #else const char* commit = "unknown"; diff --git a/source/module_base/ylm.cpp b/source/module_base/ylm.cpp index b837a7bc52..5e86b996b4 100644 --- a/source/module_base/ylm.cpp +++ b/source/module_base/ylm.cpp @@ -11,7 +11,44 @@ namespace ModuleBase { int Ylm::nlm = 0; -std::vector Ylm::ylmcoef(100); +std::vector Ylm::ylmcoef = { + 1.0 / sqrt(ModuleBase::FOUR_PI), + sqrt (3.0 / ModuleBase::FOUR_PI), + sqrt (15.0) / 2.0, + sqrt (5.0) / 2.0, + sqrt (5.0), + 1.0 / sqrt(3.0), + sqrt (5.0 / 3.0), + sqrt (35.0 / 9.0), + sqrt (7.0/3.0)/1.5, + sqrt (35.0 / 8.0), + sqrt (7.0 / 8.0), + sqrt (7.0), + 1.0 / sqrt (15.0), + sqrt (14.0 / 15.0), + sqrt (14.0 / 3.0), + sqrt(7.0)*3.0/4.0, + 9.0/4.0/sqrt(5.0), + sqrt(21.0/5.0), + sqrt(24.0/25.0), + sqrt(21.0)/2.0, + sqrt(3.0)/2.0, + 0.5/sqrt(7.0), + 1.5*sqrt(3.0/7.0), + 3.0/sqrt(2.0), + 0.6*sqrt(11.0), + 0.8*sqrt(11.0/7.0), + sqrt (33.0/8.0), + sqrt (55.0/56.0), + sqrt (33.0/7.0), + sqrt (11.0)*2.0/7.0, + sqrt (11.0)*0.75, + sqrt (11.0)*0.25, + sqrt (11.0), + 1.0/3.0/sqrt(5.0), + 2.0/3.0*sqrt(11.0/5.0), + sqrt(22.0/5.0) +}; // here Lmax == max angular momentum + 1 void Ylm::get_ylm_real( const int &Lmax, const ModuleBase::Vector3 &vec, double ylmr[] ) @@ -1287,48 +1324,7 @@ void Ylm::hes_rl_sph_harm return; } - -void Ylm::set_coefficients(void) -{ - Ylm::ylmcoef[0] = 1.0 / sqrt(ModuleBase::FOUR_PI); - Ylm::ylmcoef[1] = sqrt (3.0 / ModuleBase::FOUR_PI); - Ylm::ylmcoef[2] = sqrt (15.0) / 2.0; - Ylm::ylmcoef[3] = sqrt (5.0) / 2.0; - Ylm::ylmcoef[4] = sqrt (5.0); - Ylm::ylmcoef[5] = 1.0 / sqrt(3.0); - Ylm::ylmcoef[6] = sqrt (5.0 / 3.0); - Ylm::ylmcoef[7] = sqrt (35.0 / 9.0); - Ylm::ylmcoef[8] = sqrt (7.0/3.0)/1.5; - Ylm::ylmcoef[9] = sqrt (35.0 / 8.0); - Ylm::ylmcoef[10] = sqrt (7.0 / 8.0); - Ylm::ylmcoef[11] = sqrt (7.0); - Ylm::ylmcoef[12] = 1.0 / sqrt (15.0); - Ylm::ylmcoef[13] = sqrt (14.0 / 15.0); - Ylm::ylmcoef[14] = sqrt (14.0 / 3.0); - Ylm::ylmcoef[15] = sqrt(7.0)*3.0/4.0; - Ylm::ylmcoef[16] = 9.0/4.0/sqrt(5.0); - Ylm::ylmcoef[17] = sqrt(21.0/5.0); - Ylm::ylmcoef[18] = sqrt(24.0/25.0); - Ylm::ylmcoef[19] = sqrt(21.0)/2.0; - Ylm::ylmcoef[20] = sqrt(3.0)/2.0; - Ylm::ylmcoef[21] = 0.5/sqrt(7.0); - Ylm::ylmcoef[22] = 1.5*sqrt(3.0/7.0); - Ylm::ylmcoef[23] = 3.0/sqrt(2.0); - Ylm::ylmcoef[24] = 0.6*sqrt(11.0); - Ylm::ylmcoef[25] = 0.8*sqrt(11.0/7.0); - Ylm::ylmcoef[26] = sqrt (33.0/8.0); - Ylm::ylmcoef[27] = sqrt (55.0/56.0); - Ylm::ylmcoef[28] = sqrt (33.0/7.0); - Ylm::ylmcoef[29] = sqrt (11.0)*2.0/7.0; - Ylm::ylmcoef[30] = sqrt (11.0)*0.75; - Ylm::ylmcoef[31] = sqrt (11.0)*0.25; - Ylm::ylmcoef[32] = sqrt (11.0); - Ylm::ylmcoef[33] = 1.0/3.0/sqrt(5.0); - Ylm::ylmcoef[34] = 2.0/3.0*sqrt(11.0/5.0); - Ylm::ylmcoef[35] = sqrt(22.0/5.0); - return; -} - +void Ylm::set_coefficients (){}; /* void Ylm::test1 (void) { diff --git a/source/module_basis/module_ao/test/CMakeLists.txt b/source/module_basis/module_ao/test/CMakeLists.txt index 61a7784f76..dd3fa0f25b 100644 --- a/source/module_basis/module_ao/test/CMakeLists.txt +++ b/source/module_basis/module_ao/test/CMakeLists.txt @@ -46,7 +46,7 @@ AddTest( LIBS ${math_libs} device numerical_atomic_orbitals container formatter SOURCES 1_snap_equal_test.cpp ORB_unittest.cpp ../../../module_base/parallel_reduce.cpp ../../../module_base/parallel_global.cpp ../../../module_base/parallel_common.cpp - ../../../module_base/assoc_laguerre.cpp ../../../module_io/output_radial.cpp + ../../../module_base/assoc_laguerre.cpp ${depend_files} ) install(DIRECTORY GaAs DESTINATION ${CMAKE_CURRENT_BINARY_DIR}/../../../tests) diff --git a/source/module_basis/module_nao/hydrogen_radials.cpp b/source/module_basis/module_nao/hydrogen_radials.cpp index 2823450ccb..d21c23ccb0 100644 --- a/source/module_basis/module_nao/hydrogen_radials.cpp +++ b/source/module_basis/module_nao/hydrogen_radials.cpp @@ -1,6 +1,7 @@ #include "module_basis/module_nao/hydrogen_radials.h" #include "module_base/global_variable.h" #include "module_base/math_integral.h" +#include "module_base/atom_in.h" // for calculating slater screening constant #include #include #include @@ -13,6 +14,7 @@ HydrogenRadials& HydrogenRadials::operator=(const HydrogenRadials& rhs) void HydrogenRadials::build(const int itype, const double charge, + const bool with_slater_screening, const int nmax, const double rcut, const double dr, @@ -22,16 +24,18 @@ void HydrogenRadials::build(const int itype, const std::string strategy, std::ofstream* ptr_log) { + if(with_slater_screening) {printf("Build hydrogen_radials with Slater screening coefficients.\n");} cleanup(); itype_ = itype; symbol_ = symbol; // rcut should be determined as soon as possible... //generate_hydrogen_radials(charge, nmax, 10.0, dr, rank, ptr_log); - hydrogen(charge, nmax, dr, conv_thr, rank, strategy, ptr_log); + hydrogen(charge, with_slater_screening, nmax, dr, conv_thr, rank, strategy, ptr_log); set_rcut_max(); } std::vector HydrogenRadials::generate_hydrogen_radial_segment(const double charge, + const bool with_slater_screening, const int n, const int l, const double rmin, @@ -51,8 +55,15 @@ std::vector HydrogenRadials::generate_hydrogen_radial_segment(const doub rgrid[ir] = rmin + ir * dr; } + double screened_charge = charge; + if(with_slater_screening) + { + double sigma = slater_screening(symbol_, n, l); + screened_charge = charge - sigma; + } + double norm_factor = sqrt( - 4.0*std::pow(charge, 3)* + 4.0*std::pow(screened_charge, 3)* static_cast(this->assoc_laguerre_.factorial(n - l - 1)) / std::pow(double(n), 4) / static_cast(this->assoc_laguerre_.factorial(n + l)) / @@ -62,7 +73,7 @@ std::vector HydrogenRadials::generate_hydrogen_radial_segment(const doub for(int ir = 0; ir != ngrid; ++ir) { // Bohr radius is 1.0 - double rho = 2.0 * rgrid[ir] * charge / n / a0; + double rho = 2.0 * rgrid[ir] * screened_charge / n / a0; rvalue[ir] = norm_factor * std::pow(rho, l) * exp(-rho/2.0) * this->assoc_laguerre_.value( n, l, rho ); @@ -85,13 +96,14 @@ double HydrogenRadials::radial_norm(const std::vector rgrid, } double HydrogenRadials::generate_hydrogen_radial_toconv(const double charge, - const int n, - const int l, - const double conv_thr, - const int rank, - std::vector& rgrid, - std::vector& rvalue, - std::ofstream* ptr_log) + const bool with_slater_screening, + const int n, + const int l, + const double conv_thr, + const int rank, + std::vector& rgrid, + std::vector& rvalue, + std::ofstream* ptr_log) { double norm = 0.0; double rmax_ = 0.0; // in Bohr @@ -121,7 +133,7 @@ double HydrogenRadials::generate_hydrogen_radial_toconv(const double charge, rgrid_segment[ir] = rmin_ + ir * dr; } std::vector rvalue_segment = generate_hydrogen_radial_segment( - charge, n, l, rmin_, rmax_, dr, rank, ptr_log); + charge, with_slater_screening, n, l, rmin_, rmax_, dr, rank, ptr_log); // before push back, pop back the last element if(rgrid.size() != 0) { @@ -231,10 +243,6 @@ std::vector> HydrogenRadials::unzip_strategy(const int nmax, ++it; } } - for(auto nl_pair : nl_pairs) - { - std::cout << nl_pair.first << " " << nl_pair.second << std::endl; - } } else { @@ -266,6 +274,7 @@ void HydrogenRadials::smooth(std::vector& rgrid, std::map, std::pair, std::vector>> HydrogenRadials::generate_orb(const double charge, + const bool with_slater_screening, const int nmax, const double dr, const double conv_thr, @@ -287,6 +296,7 @@ HydrogenRadials::generate_orb(const double charge, std::vector rgrid; std::vector rvalue; double rmax_nl = generate_hydrogen_radial_toconv(charge, + with_slater_screening, n, l, conv_thr, @@ -361,6 +371,7 @@ HydrogenRadials::mapping_nl_lzeta(const int nmax, } void HydrogenRadials::hydrogen(const double charge, + const bool with_slater_screening, const int nmax, const double dr, const double conv_thr, @@ -369,7 +380,7 @@ void HydrogenRadials::hydrogen(const double charge, std::ofstream* ptr_log) { std::map, std::pair, std::vector>> orbitals = - generate_orb(charge, nmax, dr, conv_thr, rank, strategy, ptr_log); + generate_orb(charge, with_slater_screening, nmax, dr, conv_thr, rank, strategy, ptr_log); std::map, std::pair> nl_lzeta = mapping_nl_lzeta(nmax, strategy); nchi_ = orbitals.size(); @@ -389,3 +400,54 @@ void HydrogenRadials::hydrogen(const double charge, //++ichi; } } + +double HydrogenRadials::slater_screening(const std::string symbol, + const int n, + const int l) +{ + atom_in atom_db; + std::vector electron_config = atom_db.groundstate_electronconfiguration[symbol]; + int isubshell = 0; + double sigma = 0.0; + int _len = 0; + for(int n_ = 1; n_ <= n; ++n_) + { + if(n_ == n) _len += l + 1; + else _len += n_; + } + if(_len > electron_config.size()) + { + printf("Error: electron configuration is not enough for %s\n", symbol.c_str()); + printf("n = %d, l = %d\n", n, l); + exit(1); + } + // special case for 1s: for H and He, use 0.30 for 1s screening constant + if(symbol == "H") return 0.0; // only one 1s electron, no screening by "other electrons" + else if(symbol == "He") return 0.30; // only two 1s electrons, one screening the other + else if(n == 1) return 0.30; + for(int n_ = 1; n_ <= n; ++n_) + { + int lmax = (n_ == n) ? l : n_ - 1; + for(int l_ = 0; l_ <= lmax; ++l_) + { + int nelec = electron_config[isubshell]; + if(n - n_ >= 2) sigma += nelec * 1.0; + else if(n - n_ == 1) + { + double screening = (l > 1)? 1.00 : 0.85; + sigma += nelec * screening; + } + else + { + if(l_ == l) sigma += (nelec - 1) * 0.35; + else + { + double screening = (l > 1)? 1.00 : 0.35; + sigma += nelec * screening; + } + } + ++isubshell; + } + } + return sigma; +} \ No newline at end of file diff --git a/source/module_basis/module_nao/hydrogen_radials.h b/source/module_basis/module_nao/hydrogen_radials.h index 5098bb7b03..a1242fb709 100644 --- a/source/module_basis/module_nao/hydrogen_radials.h +++ b/source/module_basis/module_nao/hydrogen_radials.h @@ -73,6 +73,7 @@ class HydrogenRadials : public RadialSet /// @param ptr_log pointer to the log ofstream void build(const int itype = 0, const double charge = 1.0, + const bool with_slater_screening = false, const int nmax = 0, const double rcut = 10.0, const double dr = 0.01, @@ -107,6 +108,7 @@ class HydrogenRadials : public RadialSet /// @param ptr_log pointer to the log ofstream /// @return the rmax of present radial function double generate_hydrogen_radial_toconv(const double charge, + const bool with_slater_screening, const int n, const int l, const double conv_thr, @@ -130,6 +132,7 @@ class HydrogenRadials : public RadialSet /// @param ptr_log pointer to the log ofstream std::map, std::pair, std::vector>> generate_orb(const double charge = 1.0, + const bool with_slater_screening = false, const int nmax = 0, const double dr = 0.01, const double conv_thr = 1e-6, @@ -152,13 +155,21 @@ class HydrogenRadials : public RadialSet /// @param strategy strategy string /// @param ptr_log pointer to the log ofstream void hydrogen(const double charge = 1.0, + const bool with_slater_screening = false, const int nmax = 0, const double dr = 0.01, const double conv_thr = 1e-6, const int rank = 0, const std::string strategy = "minimal-valence", std::ofstream* ptr_log = nullptr); - + /// @brief return the Slater screening constant for calculating effective nuclear charge + /// @note hoping to get a more accurate estimation of hydrogen-like atom radial function, by including many-electron effect in this way + /// @details algorithm: https://laney.edu/pinar-alscher/wp-content/uploads/sites/219/2016/04/Slater-rules-revised.pdf + /// @param n principal quantum number + /// @param l angular momentum quantum number + double slater_screening(const std::string symbol, + const int n, + const int l); private: /// @brief generate hydrogen-like radial functions for a given n, l, in a given range [rmin, rmax] /// @param charge charge of the nucleus @@ -171,6 +182,7 @@ class HydrogenRadials : public RadialSet /// @param ptr_log pointer to the log ofstream /// @return the radial function stored in std::vector std::vector generate_hydrogen_radial_segment(const double charge = 1.0, + const bool with_slater_screening = false, const int n = 0, const int l = 0, const double rmin = 0.0, diff --git a/source/module_basis/module_nao/radial_collection.cpp b/source/module_basis/module_nao/radial_collection.cpp index cb5a40ba25..19e09f49fb 100644 --- a/source/module_basis/module_nao/radial_collection.cpp +++ b/source/module_basis/module_nao/radial_collection.cpp @@ -209,6 +209,7 @@ void RadialCollection::build(const int nfile, const std::string* const file, con void RadialCollection::build(const int ntype, const double* const charges, + const bool with_slater_screening, const int* const nmax, const std::string* symbols, const double conv_thr, @@ -223,6 +224,7 @@ void RadialCollection::build(const int ntype, radset_[itype] = new HydrogenRadials; radset_[itype]->build(itype, charges[itype], + with_slater_screening, nmax[itype], 10.0, // rcut should be determined automatically, in principle... 0.01, @@ -327,11 +329,11 @@ char RadialCollection::check_file_type(const std::string& file) const return file_type; } -void RadialCollection::to_file(const std::string& appendix) +void RadialCollection::to_file(const std::string& appendix, const std::string& format) const { for (int itype = 0; itype < ntype_; ++itype) { std::string fname = radset_[itype]->symbol() + "_" + appendix + ".orb"; - radset_[itype]->to_file(fname); + radset_[itype]->write_abacus_orb(fname); } } \ No newline at end of file diff --git a/source/module_basis/module_nao/radial_collection.h b/source/module_basis/module_nao/radial_collection.h index 8002136dec..3e1865e0f5 100644 --- a/source/module_basis/module_nao/radial_collection.h +++ b/source/module_basis/module_nao/radial_collection.h @@ -32,6 +32,7 @@ class RadialCollection /// builds the collection from quasi hydrogen radial functions void build(const int ntype, const double* const charges, + const bool with_slater_screening, const int* const nmax, const std::string* symbols = nullptr, const double conv_thr = 1e-10, @@ -123,7 +124,15 @@ class RadialCollection const bool enable_fft = false); ///@} - void to_file(const std::string&); + /** + * @brief export all RadialSet objects to a file in a given format. + * + * Supported formats: + * - "abacus_orb" (default): ABACUS Numerical atomic orbital format + */ + void to_file(const std::string& appendix, ///< file name + const std::string& format = "abacus_orb" ///< file format + ) const; private: int ntype_ = 0; ///< number of RadialSet in the collection diff --git a/source/module_basis/module_nao/radial_set.cpp b/source/module_basis/module_nao/radial_set.cpp index 1ef5174a12..259f792891 100644 --- a/source/module_basis/module_nao/radial_set.cpp +++ b/source/module_basis/module_nao/radial_set.cpp @@ -3,9 +3,9 @@ #include #include #include +#include #include "module_base/spherical_bessel_transformer.h" -#include "module_io/output_radial.h" RadialSet::~RadialSet() { @@ -193,29 +193,62 @@ void RadialSet::cleanup() index_map_ = nullptr; } -void RadialSet::to_file(const std::string& file_name, const int rank) const +void RadialSet::write_abacus_orb(const std::string& file_name, const int rank) const { - ModuleIO::OutputRadial out_radial; - out_radial.initialize(file_name); - out_radial.configure( - symbol_, - 100.0, // fake value - rcut_max_, - lmax_, - nzeta_, - int(rcut_max_/0.01) + 1, - 0.01 - ); - for(int l = 0; l <= lmax_; l++) + std::ofstream file_to; + file_to.open(file_name, std::ios::out); + + std::vector sublayers = {"S", "P", "D", "F", "G", "H", "I", "J", "K"}; + if(file_to.good()) { - for(int izeta = 0; izeta < nzeta_[l]; izeta++) + for(int i = 0; i < 75; ++i) + { + file_to << "-"; + } + file_to << std::endl; + // left aligned + file_to << std::left << std::setw(28) << "Element" << symbol_ << std::endl; + file_to << std::left << std::setw(28) << "Energy Cutoff(Ry)" << std::to_string(int(100.0)) << std::endl; + // rcut .1f, not scientific + file_to << std::left << std::setw(28) << "Radius Cutoff(a.u.)" + << std::fixed << std::setprecision(1) << rcut_max_ << std::endl; + file_to << std::left << std::setw(28) << "Lmax" << lmax_ << std::endl; + for(int l = 0; l <= lmax_; ++l) + { + std::string title = "Number of " + sublayers[l] + "orbital-->"; + file_to << std::left << std::setw(28) << title << nzeta_[l] << std::endl; + } + for(int i = 0; i < 75; ++i) + { + file_to << "-"; + } + file_to << std::endl; + file_to << "SUMMARY END\n\n"; + file_to << std::left << std::setw(28) << "Mesh" << std::setprecision(0) << int(rcut_max_/0.01) + 1 << std::endl; + file_to << std::left << std::setw(28) << "dr" << std::setprecision(2) << 0.01 << std::endl; + + for(int l = 0; l <= lmax_; l++) { - out_radial.push( - chi_[index(l, izeta)].nr(), - chi_[index(l, izeta)].rgrid(), - chi_[index(l, izeta)].rvalue() - ); + for(int izeta = 0; izeta < nzeta_[l]; izeta++) + { + file_to << std::right << std::setw(20) << "Type" + << std::right << std::setw(20) << "L" + << std::right << std::setw(20) << "N" << std::endl; + file_to << std::right << std::setw(20) << std::to_string(0) + << std::right << std::setw(20) << std::to_string(l) + << std::right << std::setw(20) << std::to_string(izeta); + for(int i = 0; i < int(rcut_max_/0.01) + 1; i++) + { + if(i % 4 == 0) + { + file_to << std::endl; + } + file_to << std::left << std::setw(22) << std::setprecision(14) << std::scientific + << chi_[index(l, izeta)].rvalue()[i]; + } + file_to << std::endl; + } } } - out_radial.finalize(); + file_to.close(); } diff --git a/source/module_basis/module_nao/radial_set.h b/source/module_basis/module_nao/radial_set.h index 9a03e58e15..b875aadc0e 100644 --- a/source/module_basis/module_nao/radial_set.h +++ b/source/module_basis/module_nao/radial_set.h @@ -64,26 +64,28 @@ class RadialSet ) {} /** - * @brief Builds from quasi hydrogen radial functions. + * @brief Builds from hydrogen-like radial functions. * - * + * Currently only HydrogenRadials objects are supposed to used this + * interface. */ - virtual void build(const int itype = 0, - const double charge = 1.0, - const int nmax = 0, - const double rcut = 10.0, - const double dr = 0.01, - const double conv_thr = 1e-6, - const int rank = 0, - const std::string symbol = "", - const std::string strategy = "minimal-valence", - std::ofstream* const ptr_log = nullptr - ) {} + virtual void build(const int = 0, ///< the element index in calculation + const double = 1.0, ///< nuclear charge + const bool = false, ///< whether to include Slater screening + const int = 0, ///< maximal principal quantum number or electrons + const double = 10.0, ///< maximal radius + const double = 0.01, ///< radial grid step + const double = 1e-6, ///< convergence threshold for norm of pseudowavefunction + const int = 0, ///< MPI rank + const std::string = "", ///< the name of the element + const std::string = "minimal-valence", ///< the strategy to generate whole set of radial functions + std::ofstream* const = nullptr ///< output file stream for logging + ) {} /** - * @brief Builds the object from pseudopotential file + * @brief Builds from pseudopotential file * - * Currently only AtomicRadials objects are supposed to used this + * Currently only PswfcRadials objects are supposed to used this * interface. */ virtual void build(const std::string&, ///< file name @@ -93,8 +95,19 @@ class RadialSet std::ofstream* const = nullptr, ///< output file stream for logging const int = 0 ///< MPI rank ) {} - - void to_file(const std::string& file_name, const int rank = 0) const; + /** + * @brief write any RadialSet object to a file in ABACUS numerical atomic orbital format. + * + * This function will write any RadialSet object to file in ABACUS numerical atomic orbital format. + * However its counterparts, + * `read_abacus_orb()` is in AtomicRadials, + * `read_beta_upf100()` and `read_beta_upf201()` are in BetaRadials, + * `read_upf_pswfc()` is in PswfcRadials, + * `read_coeff()` is in SphbesRadials, + * due to read-in procedures always vary from one to another. + */ + void write_abacus_orb(const std::string&, ///< file name + const int = 0) const; ///< MPI rank ///@} diff --git a/source/module_basis/module_nao/test/CMakeLists.txt b/source/module_basis/module_nao/test/CMakeLists.txt index 9242ea44e8..0a90a43a21 100644 --- a/source/module_basis/module_nao/test/CMakeLists.txt +++ b/source/module_basis/module_nao/test/CMakeLists.txt @@ -14,7 +14,6 @@ AddTest( ../atomic_radials.cpp ../radial_set.cpp ../numerical_radial.cpp - ../../../module_io/output_radial.cpp ../../module_ao/ORB_atomic_lm.cpp ../../module_ao/ORB_atomic.cpp LIBS ${math_libs} device base @@ -29,7 +28,6 @@ AddTest( ../numerical_radial.cpp ../../module_ao/ORB_atomic_lm.cpp ../../module_ao/ORB_atomic.cpp - ../../../module_io/output_radial.cpp LIBS ${math_libs} device base ) @@ -42,7 +40,6 @@ AddTest( ../numerical_radial.cpp ../../module_ao/ORB_atomic_lm.cpp ../../module_ao/ORB_atomic.cpp - ../../../module_io/output_radial.cpp LIBS ${math_libs} device base ) @@ -55,7 +52,6 @@ AddTest( ../numerical_radial.cpp ../../module_ao/ORB_atomic_lm.cpp ../../module_ao/ORB_atomic.cpp - ../../../module_io/output_radial.cpp LIBS ${math_libs} device base ) @@ -66,7 +62,6 @@ AddTest( ../sphbes_radials.cpp ../radial_set.cpp ../numerical_radial.cpp - ../../../module_io/output_radial.cpp ../../module_ao/ORB_atomic_lm.cpp ../../module_ao/ORB_atomic.cpp LIBS ${math_libs} device base @@ -83,7 +78,6 @@ AddTest( ../pswfc_radials.cpp ../radial_set.cpp ../numerical_radial.cpp - ../../../module_io/output_radial.cpp ../sphbes_radials.cpp ../../module_ao/ORB_atomic_lm.cpp ../../module_ao/ORB_atomic.cpp @@ -106,7 +100,6 @@ AddTest( ../two_center_bundle.cpp ../two_center_integrator.cpp ../real_gaunt_table.cpp - ../../../module_io/output_radial.cpp LIBS ${math_libs} device base container orb ) @@ -135,7 +128,6 @@ AddTest( ../radial_set.cpp ../numerical_radial.cpp ../two_center_bundle.cpp - ../../../module_io/output_radial.cpp LIBS ${math_libs} device base container orb ) @@ -155,7 +147,6 @@ AddTest( ../sphbes_radials.cpp ../radial_set.cpp ../numerical_radial.cpp - ../../../module_io/output_radial.cpp LIBS ${math_libs} device base container orb ) diff --git a/source/module_basis/module_nao/test/hydrogen_radials_test.cpp b/source/module_basis/module_nao/test/hydrogen_radials_test.cpp index d4c7b333b6..db169801d1 100644 --- a/source/module_basis/module_nao/test/hydrogen_radials_test.cpp +++ b/source/module_basis/module_nao/test/hydrogen_radials_test.cpp @@ -256,6 +256,7 @@ TEST_F(HydrogenRadialsTest, GenerateHydrogenRadialToconv) double rmax_chg1_n1l0 = hr.generate_hydrogen_radial_toconv( 1.0, + false, 1, 0, 1e-7, @@ -273,6 +274,7 @@ TEST_F(HydrogenRadialsTest, GenerateHydrogenRadialToconv) double rmax_chg4_n2l1 = hr.generate_hydrogen_radial_toconv( 4.0, + false, 2, 1, 1e-7, @@ -298,6 +300,7 @@ TEST_F(HydrogenRadialsTest, Build) hr.build( itype_, charge_, + false, nmax_, rcut_, dr_, @@ -316,6 +319,7 @@ TEST_F(HydrogenRadialsTest, Build) hr.build( itype_, charge_, + false, 4, rcut_, dr_, @@ -337,6 +341,7 @@ TEST_F(HydrogenRadialsTest, Build) hr.build( itype_, charge_, + false, 29, rcut_, dr_, @@ -355,9 +360,35 @@ TEST_F(HydrogenRadialsTest, Build) EXPECT_EQ(hr.nzeta_max(), 4); EXPECT_EQ(hr.nchi(), 7); // Cu, energy-valence, 3p 4s 3d + printf("Unittest for generating Cu energy-valence orbitals without Slater screening:\n"); hr.build( itype_, - charge_, + 29, // use the real nuclear charge for Cu + false, + 29, + rcut_, + dr_, + 1e-6, + rank_, + "Cu", + "energy-valence", + ptr_log_ + ); + // nmax = 29, energy-valence, yields 3p 4s 3d orbitals + EXPECT_EQ(hr.lmax(), 2); + EXPECT_EQ(hr.nzeta(0), 1); + EXPECT_EQ(hr.nzeta(1), 1); + EXPECT_EQ(hr.nzeta(2), 1); + EXPECT_EQ(hr.nzeta(3), 0); + EXPECT_EQ(hr.nzeta_max(), 1); + EXPECT_EQ(hr.nchi(), 3); + // test with Slater screening on Cu + // Cu, energy-valence, 3p 4s 3d + printf("Unittest for generating Cu energy-valence orbitals with Slater screening:\n"); + hr.build( + itype_, + 29, // use the real nuclear charge for Cu + true, 29, rcut_, dr_, @@ -379,6 +410,7 @@ TEST_F(HydrogenRadialsTest, Build) hr.build( itype_, charge_, + false, 4, rcut_, dr_, @@ -398,6 +430,31 @@ TEST_F(HydrogenRadialsTest, Build) EXPECT_EQ(hr.nchi(), 10); } +TEST_F(HydrogenRadialsTest, SlaterScreeningTest) +{ + HydrogenRadials hr; + double sigma = hr.slater_screening("H", 1, 0); + EXPECT_NEAR(sigma, 0.0, 1e-6); + sigma = hr.slater_screening("He", 1, 0); + EXPECT_NEAR(sigma, 0.30, 1e-6); + sigma = hr.slater_screening("F", 2, 1); + EXPECT_NEAR(sigma, 3.8, 1e-6); + sigma = hr.slater_screening("Ca", 4, 0); + EXPECT_NEAR(sigma, 17.15, 1e-6); + sigma = hr.slater_screening("Sc", 4, 0); + EXPECT_NEAR(sigma, 18, 1e-6); + sigma = hr.slater_screening("Cu", 4, 0); + EXPECT_NEAR(sigma, 25.3, 1e-6); + sigma = hr.slater_screening("Cu", 3, 2); + EXPECT_NEAR(sigma, 21.15, 1e-6); + sigma = hr.slater_screening("Pt", 6, 0); + EXPECT_NEAR(sigma, 74.45, 1e-6); + sigma = hr.slater_screening("Pt", 5, 1); + EXPECT_NEAR(sigma, 57.65, 1e-6); + sigma = hr.slater_screening("Os", 1, 0); + EXPECT_NEAR(sigma, 0.3, 1e-6); +} + int main(int argc, char** argv) { testing::InitGoogleTest(&argc, argv); diff --git a/source/module_cell/module_symmetry/symmetry.cpp b/source/module_cell/module_symmetry/symmetry.cpp index fa6c9f1eef..8b0ca16dbb 100644 --- a/source/module_cell/module_symmetry/symmetry.cpp +++ b/source/module_cell/module_symmetry/symmetry.cpp @@ -110,7 +110,45 @@ void Symmetry::analy_sys(const Lattice& lat, const Statistics& st, Atom* atoms, test_brav = true; // output the real ibrav and point group this->setgroup(this->symop, this->nop, this->real_brav); - this->getgroup(nrot_out, nrotk_out, ofs_running); + + if (GlobalV::NSPIN > 1) pricell_loop = this->magmom_same_check(atoms); + + if (!pricell_loop && GlobalV::NSPIN == 2) + {//analyze symmetry for spin-up atoms only + std::vector pos_spinup; + for (int it = 0;it < ntype;++it) + { + int na_spinup = 0; + for (int ia = 0;ia < atoms[it].na;++ia) if (atoms[it].mag[ia] > -this->epsilon) ++na_spinup; + this->na[it] = na_spinup; + //update newpos + for (int ia = 0;ia < atoms[it].na;++ia) + { + if (atoms[it].mag[ia] > -this->epsilon) + { + pos_spinup.push_back(this->newpos[3 * (istart[it] + ia)]); + pos_spinup.push_back(this->newpos[3 * (istart[it] + ia) + 1]); + pos_spinup.push_back(this->newpos[3 * (istart[it] + ia) + 2]); + } + } + // update start to spin-up configuration + if (it > 0) istart[it] = istart[it - 1] + na[it - 1]; + if (na[it] < na[itmin_type]) + { + this->itmin_type = it; + this->itmin_start = istart[it]; + } + } + this->getgroup(nrot_out, nrotk_out, ofs_running, pos_spinup.data()); + // recover na and istart + for (int it = 0;it < ntype;++it) + { + this->na[it] = atoms[it].na; + if (it > 0) istart[it] = istart[it - 1] + na[it - 1]; + } + } + else + this->getgroup(nrot_out, nrotk_out, ofs_running, this->newpos); }; if (GlobalV::CALCULATION == "cell-relax" && nrotk > 0) @@ -231,8 +269,6 @@ void Symmetry::analy_sys(const Lattice& lat, const Statistics& st, Atom* atoms, this->set_atom_map(atoms); - if (GlobalV::NSPIN > 1) pricell_loop = this->magmom_same_check(atoms); - if (GlobalV::CALCULATION == "relax") this->all_mbl = this->is_all_movable(atoms, st); delete[] newpos; @@ -765,9 +801,9 @@ void Symmetry::lattice_type( } -void Symmetry::getgroup(int &nrot, int &nrotk, std::ofstream &ofs_running) +void Symmetry::getgroup(int& nrot, int& nrotk, std::ofstream& ofs_running, double* pos) { - ModuleBase::TITLE("Symmetry","getgroup"); + ModuleBase::TITLE("Symmetry", "getgroup"); //-------------------------------------------------------------------------------- //return all possible space group operators that reproduce a lattice with basis @@ -790,7 +826,7 @@ void Symmetry::getgroup(int &nrot, int &nrotk, std::ofstream &ofs_running) for (int i = 0; i < nop; ++i) { // std::cout << "symop = " << symop[i].e11 <<" "<< symop[i].e12 <<" "<< symop[i].e13 <<" "<< symop[i].e21 <<" "<< symop[i].e22 <<" "<< symop[i].e23 <<" "<< symop[i].e31 <<" "<< symop[i].e32 <<" "<< symop[i].e33 << std::endl; - this->checksym(this->symop[i], this->gtrans[i], this->newpos); + this->checksym(this->symop[i], this->gtrans[i], pos); // std::cout << "s_flag =" <plat.Det()<optlat.Det()/this->plat.Det()); this->ncell=floor(ncell_double+0.5); - if(this->ncell != ntrans) - { - std::cout << " WARNING: PRICELL: NCELL != NTRANS !" << std::endl; - std::cout << " NCELL=" << ncell << ", NTRANS=" << ntrans << std::endl; - std::cout << " Suggest solution: Use a larger `symmetry_prec`. " << std::endl; + + auto reset_pcell = [this]() -> void { std::cout << " Now regard the structure as a primitive cell." << std::endl; this->ncell = 1; this->ptrans = std::vector >(1, ModuleBase::Vector3(0, 0, 0)); GlobalV::ofs_running << "WARNING: Original cell may have more than one primitive cells, \ but we have to treat it as a primitive cell. Use a larger `symmetry_prec`to avoid this warning." << std::endl; + }; + if (this->ncell != ntrans) + { + std::cout << " WARNING: PRICELL: NCELL != NTRANS !" << std::endl; + std::cout << " NCELL=" << ncell << ", NTRANS=" << ntrans << std::endl; + std::cout << " Suggest solution: Use a larger `symmetry_prec`. " << std::endl; + reset_pcell(); return; } if(std::abs(ncell_double-double(this->ncell)) > this->epsilon*100) { - std::cout << " ERROR: THE NUMBER OF PRIMITIVE CELL IS NOT AN INTEGER !" << std::endl; - ModuleBase::QUIT(); + std::cout << " WARNING: THE NUMBER OF PRIMITIVE CELL IS NOT AN INTEGER !" << std::endl; + std::cout << " NCELL(double)=" << ncell_double << ", NTRANS=" << ncell << std::endl; + std::cout << " Suggest solution: Use a larger `symmetry_prec`. " << std::endl; + reset_pcell(); + return; } GlobalV::ofs_running<<"Original cell was built up by "<ncell<<" primitive cells."<ncell) { - std::cout << " ERROR: Number of cells and number of vectors did not agree."; + std::cout << " WARNING: Number of cells and number of vectors did not agree."; std::cout<<"Try to change symmetry_prec in INPUT." << std::endl; - ModuleBase::QUIT(); + reset_pcell(); } return; } diff --git a/source/module_cell/module_symmetry/symmetry.h b/source/module_cell/module_symmetry/symmetry.h index bf00c758af..7dfee17663 100644 --- a/source/module_cell/module_symmetry/symmetry.h +++ b/source/module_cell/module_symmetry/symmetry.h @@ -101,7 +101,7 @@ class Symmetry : public Symmetry_Basic void change_lattice(void); - void getgroup(int &nrot, int &nrotk, std::ofstream &ofs_running); + void getgroup(int& nrot, int& nrotk, std::ofstream& ofs_running, double* pos); void checksym(ModuleBase::Matrix3& s, ModuleBase::Vector3& gtrans, double* pos); /// @brief primitive cell analysis void pricell(double* pos, const Atom* atoms); diff --git a/source/module_elecstate/module_charge/charge_mixing.h b/source/module_elecstate/module_charge/charge_mixing.h index bc6f74e4d1..98bdc277f2 100644 --- a/source/module_elecstate/module_charge/charge_mixing.h +++ b/source/module_elecstate/module_charge/charge_mixing.h @@ -99,6 +99,9 @@ class Charge_Mixing double get_mixing_gg0() const {return mixing_gg0;} Base_Mixing::Mixing* get_mixing() const {return mixing;} + // for mixing restart + int mixing_restart = 0; //which step to restart mixing during SCF + private: // mixing_data diff --git a/source/module_elecstate/module_dm/test/test_dm_R_init.cpp b/source/module_elecstate/module_dm/test/test_dm_R_init.cpp index fee15a4bd7..e3f1441898 100644 --- a/source/module_elecstate/module_dm/test/test_dm_R_init.cpp +++ b/source/module_elecstate/module_dm/test/test_dm_R_init.cpp @@ -2,6 +2,7 @@ #include "gmock/gmock.h" #include "gtest/gtest.h" +#define private public #include "module_elecstate/module_dm/density_matrix.h" #include "module_hamilt_lcao/module_hcontainer/hcontainer.h" @@ -276,6 +277,39 @@ TEST_F(DMTest, DMInit4) delete kv; } +// test for save_DMR +TEST_F(DMTest, saveDMR) +{ + // initalize a kvectors + K_Vectors* kv = nullptr; + int nspin = 2; + int nks = 4; // since nspin = 2 + kv = new K_Vectors; + kv->nks = nks; + kv->kvec_d.resize(nks); + kv->kvec_d[1].x = 0.5; + kv->kvec_d[3].x = 0.5; + // construct a DM + elecstate::DensityMatrix, double> DM(kv, paraV, nspin); + Grid_Driver gd(0, 0, 0); + DM.init_DMR(&gd, &ucell); + // construct another DM + elecstate::DensityMatrix, double> DM_test(kv, paraV, nspin); + DM_test.init_DMR(*DM.get_DMR_pointer(1)); + DM_test.save_DMR(); + EXPECT_EQ(DM_test.get_DMR_pointer(1)->get_nnr(), DM.get_DMR_pointer(1)->get_nnr()); + EXPECT_EQ(DM_test.get_DMR_pointer(1)->get_nnr(), DM_test._DMR_save[0].size()); + // add a new AtomPair, act as a relaxation + hamilt::AtomPair tmp_ap(9, 9, 1, 0, 0, paraV); + DM_test.get_DMR_pointer(1)->insert_pair(tmp_ap); + DM_test.get_DMR_pointer(1)->allocate(); + // update DMR_save + DM_test.save_DMR(); + EXPECT_EQ(DM_test.get_DMR_pointer(1)->get_nnr(), DM_test._DMR_save[0].size()); + // delete + delete kv; +} + int main(int argc, char** argv) { #ifdef __MPI diff --git a/source/module_esolver/esolver_ks.cpp b/source/module_esolver/esolver_ks.cpp index 54f430c5b0..5a978d5abc 100644 --- a/source/module_esolver/esolver_ks.cpp +++ b/source/module_esolver/esolver_ks.cpp @@ -401,8 +401,13 @@ namespace ModuleESolver hsolver_error = this->phsol->cal_hsolerror(); } } - - this->conv_elec = (drho < this->scf_thr && iter!=GlobalV::MIXING_RESTART); + // mixing will restart at this->p_chgmix->mixing_restart steps + if (drho <= GlobalV::MIXING_RESTART && GlobalV::MIXING_RESTART > 0.0 && this->p_chgmix->mixing_restart > iter) + { + this->p_chgmix->mixing_restart = iter + 1; + } + // drho will be 0 at this->p_chgmix->mixing_restart step, which is not ground state + this->conv_elec = (drho < this->scf_thr && !(iter==this->p_chgmix->mixing_restart && GlobalV::MIXING_RESTART > 0.0)); // If drho < hsolver_error in the first iter or drho < scf_thr, we do not change rho. if (drho < hsolver_error || this->conv_elec) @@ -412,25 +417,8 @@ namespace ModuleESolver else { //----------charge mixing--------------- - //before first calling mix_rho(), bandgap and cell structure can be analyzed to get better default parameters - // if(iter == 1) - // { - // double bandgap_for_autoset = 0.0; - // if (!GlobalV::TWO_EFERMI) - // { - // this->pelec->cal_bandgap(); - // bandgap_for_autoset = this->pelec->bandgap; - // } - // else - // { - // this->pelec->cal_bandgap_updw(); - // bandgap_for_autoset = std::min(this->pelec->bandgap_up, this->pelec->bandgap_dw); - // } - // p_chgmix->auto_set(bandgap_for_autoset, GlobalC::ucell); - // } - // mixing will restart after GlobalV::MIXING_RESTART steps - // So, GlobalV::MIXING_RESTART=1 means mix from scratch - if (GlobalV::MIXING_RESTART > 0 && iter == GlobalV::MIXING_RESTART - 1) + // mixing will restart after this->p_chgmix->mixing_restart steps + if (GlobalV::MIXING_RESTART > 0 && iter == this->p_chgmix->mixing_restart - 1) { // do not mix charge density } @@ -468,7 +456,7 @@ namespace ModuleESolver if(stop) break; } // notice for restart - if (GlobalV::MIXING_RESTART > 0 && iter == GlobalV::MIXING_RESTART - 1) + if (GlobalV::MIXING_RESTART > 0 && iter == this->p_chgmix->mixing_restart - 1) { std::cout<<"SCF restart after this step!"< void ESolver_KS_LCAO::eachiterinit(const int istep, const int iter) { - if (iter == 1 || iter == GlobalV::MIXING_RESTART) + if (iter == 1) + { + this->p_chgmix->init_mixing(); // init mixing + this->p_chgmix->mixing_restart = GlobalV::SCF_NMAX; + } + // for mixing restart + if (iter == this->p_chgmix->mixing_restart && GlobalV::MIXING_RESTART > 0.0) { this->p_chgmix->init_mixing(); - if (iter == GlobalV::MIXING_RESTART && GlobalV::MIXING_DMR) // for mixing_dmr + if (GlobalV::MIXING_DMR) // for mixing_dmr { // allocate memory for dmr_mdata const elecstate::DensityMatrix* dm @@ -601,7 +607,7 @@ namespace ModuleESolver // save input rho this->pelec->charge->save_rho_before_sum_band(); // save density matrix for mixing - if (GlobalV::MIXING_RESTART > 0 && GlobalV::MIXING_DMR && iter >= GlobalV::MIXING_RESTART) + if (GlobalV::MIXING_RESTART > 0 && GlobalV::MIXING_DMR && iter >= this->p_chgmix->mixing_restart) { elecstate::DensityMatrix* dm = dynamic_cast*>(this->pelec)->get_DM(); @@ -788,7 +794,7 @@ namespace ModuleESolver void ESolver_KS_LCAO::eachiterfinish(int iter) { // mix density matrix - if (GlobalV::MIXING_RESTART > 0 && iter >= GlobalV::MIXING_RESTART && GlobalV::MIXING_DMR ) + if (GlobalV::MIXING_RESTART > 0 && iter >= this->p_chgmix->mixing_restart && GlobalV::MIXING_DMR ) { elecstate::DensityMatrix* dm = dynamic_cast*>(this->pelec)->get_DM(); diff --git a/source/module_esolver/esolver_ks_pw.cpp b/source/module_esolver/esolver_ks_pw.cpp index c5c4bc684b..c19e64b08e 100644 --- a/source/module_esolver/esolver_ks_pw.cpp +++ b/source/module_esolver/esolver_ks_pw.cpp @@ -492,7 +492,13 @@ void ESolver_KS_PW::othercalculation(const int istep) template void ESolver_KS_PW::eachiterinit(const int istep, const int iter) { - if (iter == 1 || iter == GlobalV::MIXING_RESTART) + if (iter == 1) + { + this->p_chgmix->init_mixing(); + this->p_chgmix->mixing_restart = GlobalV::SCF_NMAX; + } + // for mixing restart + if (iter == this->p_chgmix->mixing_restart && GlobalV::MIXING_RESTART > 0.0) { this->p_chgmix->init_mixing(); } diff --git a/source/module_hamilt_lcao/hamilt_lcaodft/operator_lcao/test/CMakeLists.txt b/source/module_hamilt_lcao/hamilt_lcaodft/operator_lcao/test/CMakeLists.txt index 28b02effa0..5aeae5314c 100644 --- a/source/module_hamilt_lcao/hamilt_lcaodft/operator_lcao/test/CMakeLists.txt +++ b/source/module_hamilt_lcao/hamilt_lcaodft/operator_lcao/test/CMakeLists.txt @@ -9,7 +9,6 @@ AddTest( ../../../../module_basis/module_ao/parallel_2d.cpp ../../../../module_basis/module_ao/parallel_orbitals.cpp ../../../../module_basis/module_ao/ORB_atomic_lm.cpp tmp_mocks.cpp ../../../../module_hamilt_general/operator.cpp - ../../../../module_io/output_radial.cpp ) AddTest( @@ -20,7 +19,6 @@ AddTest( ../../../../module_basis/module_ao/parallel_2d.cpp ../../../../module_basis/module_ao/parallel_orbitals.cpp ../../../../module_basis/module_ao/ORB_atomic_lm.cpp tmp_mocks.cpp ../../../../module_hamilt_general/operator.cpp - ../../../../module_io/output_radial.cpp ) AddTest( @@ -31,7 +29,6 @@ AddTest( ../../../../module_basis/module_ao/parallel_2d.cpp ../../../../module_basis/module_ao/parallel_orbitals.cpp ../../../../module_basis/module_ao/ORB_atomic_lm.cpp tmp_mocks.cpp ../../../../module_hamilt_general/operator.cpp - ../../../../module_io/output_radial.cpp ) AddTest( @@ -42,7 +39,6 @@ AddTest( ../../../../module_basis/module_ao/parallel_2d.cpp ../../../../module_basis/module_ao/parallel_orbitals.cpp ../../../../module_basis/module_ao/ORB_atomic_lm.cpp tmp_mocks.cpp ../../../../module_hamilt_general/operator.cpp - ../../../../module_io/output_radial.cpp ) AddTest( @@ -53,7 +49,6 @@ AddTest( ../../../../module_basis/module_ao/parallel_2d.cpp ../../../../module_basis/module_ao/parallel_orbitals.cpp ../../../../module_basis/module_ao/ORB_atomic_lm.cpp tmp_mocks.cpp ../../../../module_hamilt_general/operator.cpp - ../../../../module_io/output_radial.cpp ) install(FILES parallel_operator_tests.sh DESTINATION ${CMAKE_CURRENT_BINARY_DIR}) diff --git a/source/module_io/CMakeLists.txt b/source/module_io/CMakeLists.txt index 618f300170..bcc98f797f 100644 --- a/source/module_io/CMakeLists.txt +++ b/source/module_io/CMakeLists.txt @@ -28,7 +28,6 @@ list(APPEND objects output_log.cpp output_rho.cpp output_potential.cpp - output_radial.cpp parameter_pool.cpp para_json.cpp ) diff --git a/source/module_io/input.cpp b/source/module_io/input.cpp index 037e4f3db9..c68311c5ae 100644 --- a/source/module_io/input.cpp +++ b/source/module_io/input.cpp @@ -75,7 +75,8 @@ void Input::Init(const std::string& fn) #else const char* version = "unknown"; #endif -#ifdef COMMIT +#ifdef COMMIT_INFO +#include "commit.h" const char* commit = COMMIT; #else const char* commit = "unknown"; @@ -305,7 +306,7 @@ void Input::Default(void) mixing_mode = "broyden"; mixing_beta = -10; mixing_ndim = 8; - mixing_restart = 0; + mixing_restart = 0.0; mixing_gg0 = 1.00; // use Kerker defaultly mixing_beta_mag = -10.0; // only set when nspin == 2 || nspin == 4 mixing_gg0_mag = 0.0; // defaultly exclude Kerker from mixing magnetic density @@ -3136,8 +3137,15 @@ void Input::Default_2(void) // jiyy add 2019-08-04 } if(qo_screening_coeff.size() != ntype) { - double default_screening_coeff = (qo_screening_coeff.size() == 1)? qo_screening_coeff[0]: 0.1; - qo_screening_coeff.resize(ntype, default_screening_coeff); + if(qo_basis == "pswfc") + { + double default_screening_coeff = (qo_screening_coeff.size() == 1)? qo_screening_coeff[0]: 0.1; + qo_screening_coeff.resize(ntype, default_screening_coeff); + } + else + { + // if length of qo_screening_coeff is not 0, turn on Slater screening + } } if(qo_strategy.size() != ntype) { @@ -3352,7 +3360,7 @@ void Input::Bcast() Parallel_Common::bcast_string(mixing_mode); Parallel_Common::bcast_double(mixing_beta); Parallel_Common::bcast_int(mixing_ndim); - Parallel_Common::bcast_int(mixing_restart); + Parallel_Common::bcast_double(mixing_restart); Parallel_Common::bcast_double(mixing_gg0); // mohan add 2014-09-27 Parallel_Common::bcast_double(mixing_beta_mag); Parallel_Common::bcast_double(mixing_gg0_mag); diff --git a/source/module_io/input.h b/source/module_io/input.h index 8629e2746b..ececbaeeaa 100644 --- a/source/module_io/input.h +++ b/source/module_io/input.h @@ -233,8 +233,8 @@ class Input std::string mixing_mode; // "plain","broyden",... double mixing_beta; // 0 : no_mixing int mixing_ndim; // used in Broyden method - int mixing_restart; - double mixing_gg0; // used in kerker method. mohan add 2014-09-27 + double mixing_restart; // mixing will restart once if drho is smaller than mixing_restart + double mixing_gg0; // used in kerker method double mixing_beta_mag; double mixing_gg0_mag; double mixing_gg0_min; diff --git a/source/module_io/json_output/CMakeLists.txt b/source/module_io/json_output/CMakeLists.txt index 7dea34b51b..b54336fc4e 100644 --- a/source/module_io/json_output/CMakeLists.txt +++ b/source/module_io/json_output/CMakeLists.txt @@ -2,6 +2,7 @@ list(APPEND objects_json abacusjson.cpp general_info.cpp init_info.cpp + readin_info.cpp ) add_library( diff --git a/source/module_io/json_output/general_info.cpp b/source/module_io/json_output/general_info.cpp index 1ac054f34b..479ff5fd16 100644 --- a/source/module_io/json_output/general_info.cpp +++ b/source/module_io/json_output/general_info.cpp @@ -19,7 +19,8 @@ void gen_general_info(Input *input) #else const std::string version = "unknown"; #endif -#ifdef COMMIT +#ifdef COMMIT_INFO +#include "commit.h" const std::string commit = COMMIT; #else const std::string commit = "unknown"; diff --git a/source/module_io/json_output/general_info.h b/source/module_io/json_output/general_info.h index afed718a06..7bee072098 100644 --- a/source/module_io/json_output/general_info.h +++ b/source/module_io/json_output/general_info.h @@ -1,6 +1,8 @@ #include "module_io/input.h" -//Add json objects to gener_info +/** +* @brief In this part of the code to complete the general_info part of the json tree. +*/ namespace Json { #ifdef __RAPIDJSON diff --git a/source/module_io/json_output/init_info.cpp b/source/module_io/json_output/init_info.cpp index ff2e7107f5..3a8e469cbe 100644 --- a/source/module_io/json_output/init_info.cpp +++ b/source/module_io/json_output/init_info.cpp @@ -11,19 +11,6 @@ namespace Json #ifdef __RAPIDJSON -// void gen_init(ModuleSymmetry::Symmetry *symm,Atom *atoms){ -// std::string pgname = symm->pgname; -// std::string spgname = symm->spgname; -// AbacusJson::add_json({"init", "pgname"}, pgname,false); -// AbacusJson::add_json({"init", "spgname"}, spgname,false); - - -// int numAtoms = atoms->na; -// AbacusJson::add_json({"init", "natom"}, numAtoms,false); -// AbacusJson::add_json({"init", "nband"}, INPUT.nbands,false); - - -// } void gen_init(UnitCell *ucell){ std::string pgname = ucell->symm.pgname; @@ -36,15 +23,19 @@ void gen_init(UnitCell *ucell){ AbacusJson::add_json({"init", "natom"}, numAtoms,false); AbacusJson::add_json({"init", "nband"}, GlobalV::NBANDS,false); - int ntype = ucell->ntype; - AbacusJson::add_json({"init", "nelectron"}, ntype,false); - + int ntype = ucell->ntype,nelec_total=0; + AbacusJson::add_json({"init", "nelectron"}, 0,false); for (int it = 0; it < ntype; it++) { std::string label = ucell->atoms[it].label; int number = ucell->atoms[it].ncpp.zv; + + nelec_total+=ucell->atoms[it].ncpp.zv * ucell->atoms[it].na; AbacusJson::add_json({"init", "nelectron_each_type",label}, number,false); } + + AbacusJson::add_json({"init", "nelectron"}, nelec_total,false); + } @@ -54,5 +45,89 @@ void add_nkstot(int nkstot,int nkstot_ibz){ Json::AbacusJson::add_json({"init", "nkstot_ibz"}, nkstot_ibz,false); } + +void gen_stru(UnitCell *ucell){ + + int ntype = ucell->ntype; + + + + //array of pseudopotential file + std::string* pseudo_fn = ucell->pseudo_fn; + + //array of orbital file + std::string* orbital_fn = ucell->orbital_fn; + + + //add atom element,orbital file and pseudopotential file + for(int i=0;iatoms[i].label; + + std::string atom_element = ucell->atoms[i].ncpp.psd; + + Json::jsonValue element_obj(JobjectType); + element_obj.JaddStringKV(atom_label,atom_element); + Json::AbacusJson::add_json({"init","element"}, element_obj,false); + + + std::string orbital_str = GlobalV::global_orbital_dir + orbital_fn[i]; + if(!orbital_str.compare("")){ + Json::jsonValue nullValue; + nullValue.SetNull(); + Json::AbacusJson::add_json({"init","orb",atom_label}, nullValue,false); + }else { + Json::AbacusJson::add_json({"init","orb",atom_label}, orbital_str,false); + } + std::string pseudo_str = pseudo_fn[i]; + Json::AbacusJson::add_json({"init","pp",atom_label}, pseudo_str,false); + + + } + + + //atom coordinate, mag and label + double lat0 = ucell->lat0; + std::string* label = ucell->atom_label; + for(int i=0;i* tau = ucell->atoms[i].tau; + int na = ucell->atoms[i].na; + for(int j=0;jatoms[i].mag[j],true); + + std::string str = label[i]; + Json::AbacusJson::add_json({"init","label"}, str,true); + } + + } + + + //cell + { + Json::jsonValue cellArray1(JarrayType); + Json::jsonValue cellArray2(JarrayType); + Json::jsonValue cellArray3(JarrayType); + cellArray1.JPushBack(ucell->latvec.e11); + cellArray1.JPushBack(ucell->latvec.e12); + cellArray1.JPushBack(ucell->latvec.e13); + cellArray2.JPushBack(ucell->latvec.e21); + cellArray2.JPushBack(ucell->latvec.e22); + cellArray2.JPushBack(ucell->latvec.e23); + cellArray3.JPushBack(ucell->latvec.e31); + cellArray3.JPushBack(ucell->latvec.e32); + cellArray3.JPushBack(ucell->latvec.e33); + Json::AbacusJson::add_json({"init","cell"}, cellArray1,true); + Json::AbacusJson::add_json({"init","cell"}, cellArray2,true); + Json::AbacusJson::add_json({"init","cell"}, cellArray3,true); + } + return; +} + + #endif } // namespace Json \ No newline at end of file diff --git a/source/module_io/json_output/init_info.h b/source/module_io/json_output/init_info.h index 90e6b62e59..324d8e498c 100644 --- a/source/module_io/json_output/init_info.h +++ b/source/module_io/json_output/init_info.h @@ -3,15 +3,27 @@ #include "module_cell/unitcell.h" -//Add json objects to init +/** +* @brief In this part of the code to complete the init part of the json tree. +*/ namespace Json { #ifdef __RAPIDJSON // void gen_init(ModuleSymmetry::Symmetry *symm,Atom *atoms); + +/** +* @param ucell: ucell for reading json parameters. +*/ void gen_init(UnitCell *ucell); +/** +* @param nkstot,nkstot_ibz: two param in json tree +*/ void add_nkstot(int nkstot,int nkstot_ibz); - +/** +* @param ucell: ucell for reading structure init in abacus. +*/ +void gen_stru(UnitCell *ucell); #endif } \ No newline at end of file diff --git a/source/module_io/json_output/readin_info.cpp b/source/module_io/json_output/readin_info.cpp new file mode 100644 index 0000000000..d43a23a1c0 --- /dev/null +++ b/source/module_io/json_output/readin_info.cpp @@ -0,0 +1,17 @@ +#include "readin_info.h" +#include "module_io/input.h" +#include "../para_json.h" +#include "abacusjson.h" + + +//Add json objects to init +namespace Json +{ + +#ifdef __RAPIDJSON + + + + +#endif +} // namespace Json \ No newline at end of file diff --git a/source/module_io/json_output/readin_info.h b/source/module_io/json_output/readin_info.h new file mode 100644 index 0000000000..8c8dfb3321 --- /dev/null +++ b/source/module_io/json_output/readin_info.h @@ -0,0 +1,18 @@ +#include "module_cell/module_symmetry/symmetry.h" +#include "module_cell/atom_spec.h" +#include "module_cell/unitcell.h" + + + +/** +* @brief In this part of the code to complete the readin part of the json tree. +* @param ucell: ucell for reading json parameters +*/ +namespace Json +{ +#ifdef __RAPIDJSON + + + +#endif +} \ No newline at end of file diff --git a/source/module_io/json_output/test/CMakeLists.txt b/source/module_io/json_output/test/CMakeLists.txt index ec09db3965..5628afca48 100644 --- a/source/module_io/json_output/test/CMakeLists.txt +++ b/source/module_io/json_output/test/CMakeLists.txt @@ -7,6 +7,6 @@ remove_definitions(-DUSE_PAW) AddTest( TARGET io_json_output_json LIBS ${math_libs} base device cell_info - SOURCES para_json_test.cpp ../general_info.cpp ../init_info.cpp + SOURCES para_json_test.cpp ../general_info.cpp ../init_info.cpp ../readin_info.cpp ../../para_json.cpp ../abacusjson.cpp ../../input.cpp ../../output.cpp ) \ No newline at end of file diff --git a/source/module_io/json_output/test/para_json_test.cpp b/source/module_io/json_output/test/para_json_test.cpp index c3baefd4cd..2e40d14654 100644 --- a/source/module_io/json_output/test/para_json_test.cpp +++ b/source/module_io/json_output/test/para_json_test.cpp @@ -6,6 +6,7 @@ #include "module_io/para_json.h" #include "../init_info.h" #include "../general_info.h" +#include "../readin_info.h" #include "version.h" /************************************************ @@ -89,6 +90,7 @@ TEST(AbacusJsonTest, AddJson) } // Validate json parameters in doc objects + ASSERT_EQ(Json::AbacusJson::doc["array"][0]["int"].GetInt(), 0); ASSERT_STREQ(Json::AbacusJson::doc["array"][0]["string"].GetString(), "0"); ASSERT_STREQ(Json::AbacusJson::doc["array"][0]["0"].GetString(), "string"); @@ -257,7 +259,6 @@ Magnetism::~Magnetism() { delete[] this->start_magnetization; } - TEST(AbacusJsonTest, InitInfo) { UnitCell ucell; @@ -271,9 +272,9 @@ TEST(AbacusJsonTest, InitInfo) ucell.ntype = 3; GlobalV::NBANDS = 10; - for(int i=0;i<3;i++){ - ucell.atoms[i].label = std::to_string(i); - ucell.atoms[i].ncpp.zv = i*3; + for(int i=0;i<1;i++){ + ucell.atoms[i].label = "Si"; + ucell.atoms[i].ncpp.zv = 3; } @@ -296,8 +297,92 @@ TEST(AbacusJsonTest, InitInfo) ASSERT_STREQ(Json::AbacusJson::doc["init"]["point_group"].GetString(), "T_d"); ASSERT_STREQ(Json::AbacusJson::doc["init"]["point_group_in_space"].GetString(), "O_h"); - ASSERT_EQ(Json::AbacusJson::doc["init"]["nelectron_each_type"]["0"].GetInt(), 0); - ASSERT_EQ(Json::AbacusJson::doc["init"]["nelectron_each_type"]["1"].GetInt(), 3); - ASSERT_EQ(Json::AbacusJson::doc["init"]["nelectron_each_type"]["2"].GetInt(), 6); - + ASSERT_EQ(Json::AbacusJson::doc["init"]["nelectron_each_type"]["Si"].GetInt(), 3); + } + + + +TEST(AbacusJsonTest, Init_stru_test){ + //init ucell + UnitCell ucell; + + Atom atomlist[1]; + std::string label[1]; + + ModuleBase::Matrix3 latvec; + latvec.e11=0.1; + latvec.e12=0.1; + latvec.e13=0.1; + + latvec.e21=0.2; + latvec.e22=0.2; + latvec.e23=0.2; + + latvec.e31=0.3; + latvec.e32=0.3; + latvec.e33=0.3; + ucell.latvec = latvec; + + + double lat0 = 10.0; + ucell.ntype = 1; + ucell.pseudo_fn = new std::string[1]; + ucell.orbital_fn = new std::string[1]; + ucell.atoms = atomlist; + ucell.atom_label = new std::string[1]; + ucell.lat0 = lat0; + + ModuleBase::Vector3 tau[2]; + + Json::AbacusJson::doc.Parse("{}"); + + double mag[2]; + //fill ucell + for(int i=0;i<1;i++){ + ucell.atom_label[i]="Si"; + atomlist[i].na = 2; + atomlist[i].label = "Fe"; + ucell.pseudo_fn[i] = "si.ufp"; + ucell.atoms[i].tau = new ModuleBase::Vector3[2]; + atomlist[i].mag = new double[2]; + for(int j=0;j sublayers = {"S", "P", "D", "F", "G", "H", "I", "J", "K"}; - for(int i = 0; i < 75; ++i) - { - file_to_ << "-"; - } - file_to_ << "\n"; - // left aligned - file_to_ << std::left << std::setw(28) << "Element" << symbol_ << "\n"; - file_to_ << std::left << std::setw(28) << "Energy Cutoff(Ry)" << std::to_string(int(ecut_)) << "\n"; - // rcut .1f, not scientific - file_to_ << std::left << std::setw(28) << "Radius Cutoff(a.u.)" << std::fixed << std::setprecision(1) << rcut_ << "\n"; - file_to_ << std::left << std::setw(28) << "Lmax" << lmax_ << "\n"; - for(int i = 0; i < lmax_ + 1; ++i) - { - std::string title = "Number of " + sublayers[i] + "orbital-->"; - file_to_ << std::left << std::setw(28) << title << l_nchi_[i] << "\n"; - } - for(int i = 0; i < 75; ++i) - { - file_to_ << "-"; - } - file_to_ << "\n"; - file_to_ << "SUMMARY END\n\n"; - file_to_ << std::left << std::setw(28) << "Mesh" << std::setprecision(0) << ngrid_ << "\n"; - file_to_ << std::left << std::setw(28) << "dr" << std::setprecision(2) << dr_ << "\n"; - } - else - { - ModuleBase::WARNING_QUIT("ModuleIO::OutputRadial::write_header", "file is not good"); - } -} - -void ModuleIO::OutputRadial::configure(const std::string& symbol, - const double& ecut, - const double& rcut, - const int& lmax, - const int* l_nchi, - const int& ngrid, - const double& dr - ) -{ - symbol_ = symbol; - ecut_ = ecut; - rcut_ = rcut; - lmax_ = lmax; - l_nchi_.resize(lmax_ + 1); - for(int i = 0; i < lmax_ + 1; ++i) - { - l_nchi_[i] = l_nchi[i]; - } - ngrid_ = ngrid; - dr_ = dr; - write_header(); -} - -void ModuleIO::OutputRadial::push(int ngrid, - const double* rgrid, - const double* chi) -{ - if(ngrid == 0) - { - ModuleBase::WARNING("ModuleIO::OutputRadial::push", "ngrid is 0, use the configured one"); - ngrid = ngrid_; - } - if(rgrid == nullptr) - { - ModuleBase::WARNING("ModuleIO::OutputRadial::push", "rgrid is nullptr, use the configured one"); - } - - file_to_ << std::right << std::setw(20) << "Type" - << std::right << std::setw(20) << "L" - << std::right << std::setw(20) << "N" << "\n"; - file_to_ << std::right << std::setw(20) - << std::to_string(0) - << std::right < lmax_ + 1) - { - ModuleBase::WARNING_QUIT("ModuleIO::OutputRadial::push", "current_l_ is out of range"); - } -} - -void ModuleIO::OutputRadial::finalize() -{ - file_to_.close(); -} \ No newline at end of file diff --git a/source/module_io/output_radial.h b/source/module_io/output_radial.h deleted file mode 100644 index a488c2700d..0000000000 --- a/source/module_io/output_radial.h +++ /dev/null @@ -1,70 +0,0 @@ -#ifndef RADIAL_IO_H_ -#define RADIAL_IO_H_ - -#include -#include - -namespace ModuleIO -{ - class OutputRadial - { - public: - OutputRadial() {}; - ~OutputRadial() {}; - /// @brief will create a ofstream object and open the file - /// @param filename name of the file to be opened - void initialize(const std::string& filename); - /// @brief configure basic information of the file - /// @param symbol element symbol - /// @param ecut energy cutoff of numerical atomic orbital - /// @param rcut radius cutoff of numerical atomic orbital - /// @param lmax maximum angular momentum - /// @param l_nchi number of radial functions for each angular momentum - /// @param ngrid number of grid points - /// @param dr grid spacing - void configure(const std::string& symbol = "H", - const double& ecut = 100.0, - const double& rcut = 10.0, - const int& lmax = 0, - const int* l_nchi = nullptr, - const int& ngrid = 0, - const double& dr = 0.01 - ); - /// @brief import radialset data in - /// @param ngrid number of grid points - /// @param rgrid radial grid - /// @param chi zeta function - void push(int ngrid = 0, - const double* rgrid = nullptr, - const double* chi = nullptr); - /// @brief will close the file and delete the ofstream object - void finalize(); - - // getters - std::string symbol() const {return symbol_;}; - double ecut() const {return ecut_;}; - double rcut() const {return rcut_;}; - int lmax() const {return lmax_;}; - std::vector l_nchi() const {return l_nchi_;}; - int ngrid() const {return ngrid_;}; - double dr() const {return dr_;}; - int current_l() const {return current_l_;}; - int current_chi() const {return current_chi_;}; - - private: - std::ofstream file_to_; - std::string symbol_ = "H"; - double ecut_ = 100.0; - double rcut_ = 10.0; - int lmax_ = 0; - std::vector l_nchi_; - int ngrid_ = 0; - double dr_ = 0.01; - - int current_l_ = 0; - int current_chi_ = 0; - - void write_header(); - }; -} // namespace ModuleIO -#endif // RADIAL_IO_H_ \ No newline at end of file diff --git a/source/module_io/para_json.cpp b/source/module_io/para_json.cpp index 5e79431069..d0af231b2b 100644 --- a/source/module_io/para_json.cpp +++ b/source/module_io/para_json.cpp @@ -11,6 +11,7 @@ #include "json_output/abacusjson.h" #include "json_output/general_info.h" #include "json_output/init_info.h" +#include "json_output/readin_info.h" #endif // __RAPIDJSON @@ -40,10 +41,22 @@ void create_Json(UnitCell *ucell,Input *input){ #ifdef __RAPIDJSON gen_general_info(input); gen_init(ucell); + // gen_stru(ucell); #endif json_output(); } +void gen_stru_wrapper(UnitCell *ucell){ +#ifdef __RAPIDJSON +#ifdef __MPI + if (GlobalV::MY_RANK == 0) + gen_stru(ucell); +#elif + gen_stru(ucell); +#endif +#endif +} + void convert_time(std::time_t time_now, std::string& time_str) { std::tm* tm = std::localtime(&time_now); diff --git a/source/module_io/para_json.h b/source/module_io/para_json.h index bba5d075c3..9fb665e587 100644 --- a/source/module_io/para_json.h +++ b/source/module_io/para_json.h @@ -18,4 +18,6 @@ void json_output(); // Convert time_t to string void convert_time(std::time_t time_now, std::string& time_str); +// generate struture wrapper function +void gen_stru_wrapper(UnitCell *ucell); } // namespace Json diff --git a/source/module_io/parameter_pool.cpp b/source/module_io/parameter_pool.cpp index 0a45e2b39d..3590c38eee 100644 --- a/source/module_io/parameter_pool.cpp +++ b/source/module_io/parameter_pool.cpp @@ -820,7 +820,7 @@ void input_parameters_set(std::map input_parameters } else if (input_parameters.count("mixing_restart") != 0) { - INPUT.mixing_restart = *static_cast(input_parameters["mixing_restart"].get()); + INPUT.mixing_restart = *static_cast(input_parameters["mixing_restart"].get()); } else if (input_parameters.count("mixing_gg0") != 0) { diff --git a/source/module_io/test/CMakeLists.txt b/source/module_io/test/CMakeLists.txt index aee5e4b47d..35cac76c0b 100644 --- a/source/module_io/test/CMakeLists.txt +++ b/source/module_io/test/CMakeLists.txt @@ -159,17 +159,9 @@ AddTest( ../../module_cell/atom_spec.cpp ../../module_cell/parallel_kpoints.cpp ../../module_cell/test/support/mock_unitcell.cpp - ../output_radial.cpp ) endif() -AddTest( - TARGET output_radial_test - LIBS base ${math_libs} device - SOURCES output_radial_test.cpp ../output_radial.cpp -) - - AddTest( TARGET read_wfc_pw_test LIBS base ${math_libs} device planewave diff --git a/source/module_io/test/input_conv_test.cpp b/source/module_io/test/input_conv_test.cpp index fcdc4b6b69..b9742c6819 100644 --- a/source/module_io/test/input_conv_test.cpp +++ b/source/module_io/test/input_conv_test.cpp @@ -184,7 +184,7 @@ TEST_F(InputConvTest, Conv) EXPECT_TRUE(GlobalV::decay_grad_switch); EXPECT_EQ(GlobalV::sc_file, "sc.json"); - EXPECT_EQ(GlobalV::MIXING_RESTART,0); + EXPECT_EQ(GlobalV::MIXING_RESTART,0.0); EXPECT_EQ(GlobalV::MIXING_DMR,false); } diff --git a/source/module_io/test/input_test_para.cpp b/source/module_io/test/input_test_para.cpp index 7d53b57311..46e9c5d729 100644 --- a/source/module_io/test/input_test_para.cpp +++ b/source/module_io/test/input_test_para.cpp @@ -382,7 +382,7 @@ TEST_F(InputParaTest, Bcast) EXPECT_TRUE(INPUT.mdp.dump_virial); EXPECT_FALSE(INPUT.mixing_tau); EXPECT_FALSE(INPUT.mixing_dftu); - EXPECT_EQ(INPUT.mixing_restart,0); + EXPECT_EQ(INPUT.mixing_restart,0.0); EXPECT_EQ(INPUT.mixing_dmr,false); EXPECT_EQ(INPUT.out_bandgap, 0); EXPECT_EQ(INPUT.out_mat_t, 0); diff --git a/source/module_io/test/output_radial_test.cpp b/source/module_io/test/output_radial_test.cpp deleted file mode 100644 index 19fabf0963..0000000000 --- a/source/module_io/test/output_radial_test.cpp +++ /dev/null @@ -1,87 +0,0 @@ -#include -#include "../output_radial.h" -#include - -class OutputRadialTest : public ::testing::Test -{ - protected: - virtual void SetUp() - { - symbol_ = "H"; - ecut_ = 100.0; - rcut_ = 10.0; - lmax_ = 0; - l_nchi_ = {1}; - ngrid_ = 101; - dr_ = 0.01; - rgrid_ = new double[ngrid_]; - chi_ = new double[ngrid_]; - for(int i = 0; i < ngrid_; ++i) - { - rgrid_[i] = i * dr_; - chi_[i] = std::exp(-rgrid_[i]); - } - } - - virtual void TearDown() - { - delete[] rgrid_; - delete[] chi_; - } - - std::string symbol_; - double ecut_; - double rcut_; - int lmax_; - std::vector l_nchi_; - int ngrid_; - double dr_; - double* rgrid_; - double* chi_; -}; - -TEST_F(OutputRadialTest, initialize) -{ - ModuleIO::OutputRadial output_radial; - output_radial.initialize("test_output_radial.dat"); - EXPECT_EQ(output_radial.symbol(), "H"); - EXPECT_EQ(output_radial.ecut(), 100.0); - EXPECT_EQ(output_radial.rcut(), 10.0); - EXPECT_EQ(output_radial.lmax(), 0); - EXPECT_EQ(output_radial.ngrid(), 0); - EXPECT_EQ(output_radial.dr(), 0.01); - EXPECT_EQ(output_radial.current_l(), 0); - EXPECT_EQ(output_radial.current_chi(), 0); - output_radial.finalize(); -} - -TEST_F(OutputRadialTest, configure) -{ - ModuleIO::OutputRadial output_radial; - output_radial.initialize("test_output_radial.dat"); - output_radial.configure(symbol_, ecut_, rcut_, lmax_, l_nchi_.data(), ngrid_, dr_); - EXPECT_EQ(output_radial.symbol(), symbol_); - EXPECT_EQ(output_radial.ecut(), ecut_); - EXPECT_EQ(output_radial.rcut(), rcut_); - EXPECT_EQ(output_radial.lmax(), lmax_); - for(int i = 0; i < lmax_ + 1; ++i) - { - EXPECT_EQ(output_radial.l_nchi()[i], l_nchi_[i]); - } - EXPECT_EQ(output_radial.ngrid(), ngrid_); - EXPECT_EQ(output_radial.dr(), dr_); - EXPECT_EQ(output_radial.current_l(), 0); - EXPECT_EQ(output_radial.current_chi(), 0); - output_radial.finalize(); -} - -TEST_F(OutputRadialTest, push) -{ - ModuleIO::OutputRadial output_radial; - output_radial.initialize("test_output_radial.dat"); - output_radial.configure(symbol_, ecut_, rcut_, lmax_, l_nchi_.data(), ngrid_, dr_); - output_radial.push(ngrid_, rgrid_, chi_); - EXPECT_EQ(output_radial.current_l(), 1); - EXPECT_EQ(output_radial.current_chi(), 0); - output_radial.finalize(); -} diff --git a/source/module_io/test/to_qo_test.cpp b/source/module_io/test/to_qo_test.cpp index 25cf5bae3c..2d792913e1 100644 --- a/source/module_io/test/to_qo_test.cpp +++ b/source/module_io/test/to_qo_test.cpp @@ -59,6 +59,8 @@ void define_fcc_cell(UnitCell& ucell) GlobalV::global_out_dir = "./"; GlobalV::qo_screening_coeff = {0.1, 0.1}; + GlobalV::qo_thr = 1e-6; + GlobalV::ofs_running = std::ofstream("unittest.log"); } void define_sc_cell(UnitCell& ucell) @@ -92,6 +94,8 @@ void define_sc_cell(UnitCell& ucell) GlobalV::global_out_dir = "./"; GlobalV::qo_screening_coeff = {0.1}; + GlobalV::qo_thr = 1e-6; + GlobalV::ofs_running = std::ofstream("unittest.log"); } class toQOTest : public testing::Test @@ -162,7 +166,7 @@ TEST_F(toQOTest, BuildNao) define_fcc_cell(ucell); toQO tqo("hydrogen", {"minimal-nodeless", "minimal-nodeless"}); tqo.unwrap_unitcell(&ucell); - tqo.build_nao(ucell.ntype, ucell.orbital_fn); + tqo.build_nao(ucell.ntype, "./", ucell.orbital_fn, 0); EXPECT_EQ(tqo.p_nao()->nchi(), 10); // not (l, m)-resoluted EXPECT_EQ(tqo.nphi(), 26); // (l, m)-resoluted } @@ -172,7 +176,13 @@ TEST_F(toQOTest, BuildHydrogenMinimal) define_fcc_cell(ucell); toQO tqo("hydrogen", {"minimal-nodeless", "minimal-nodeless"}); tqo.unwrap_unitcell(&ucell); - tqo.build_ao(ucell.ntype); + tqo.build_ao(ucell.ntype, + "./", + ucell.pseudo_fn, + {}, + GlobalV::qo_thr, + GlobalV::ofs_running, + 0); EXPECT_EQ(tqo.p_ao()->nchi(), 5); // Si: 1s, 2p, 3d, C: 1s, 2p EXPECT_EQ(tqo.nchi(), 13); // Si: 1s, 2px, 2py, 2pz, 3dz2, 3dxz, 3dyz, 3dx2-y2, 3dxy, C: 1s, 2px, 2py, 2pz tqo.p_ao()->to_file("special_use_unittest"); @@ -199,13 +209,22 @@ TEST_F(toQOTest, ScanSupercellForAtom) define_fcc_cell(ucell); toQO tqo("hydrogen", {"minimal-nodeless", "minimal-nodeless"}); tqo.unwrap_unitcell(&ucell); - tqo.build_nao(ucell.ntype, ucell.orbital_fn); + tqo.build_nao(ucell.ntype, // ntype + "./", // orbital_dir + ucell.orbital_fn,// orbital_fn + 0); // rank std::vector nmax = std::vector(ucell.ntype); for(int itype = 0; itype < ucell.ntype; itype++) { nmax[itype] = tqo.atom_database().principle_quantum_number[tqo.symbols()[itype]]; } - tqo.build_ao(ucell.ntype); + tqo.build_ao(ucell.ntype, + "./", + ucell.pseudo_fn, + {}, + GlobalV::qo_thr, + GlobalV::ofs_running, + 0); std::vector> n1n2n3 = tqo.scan_supercell_for_atom(0, 0); EXPECT_EQ(n1n2n3.size(), 13); // 13 = 3*3*3 - 2 - 3*4 } @@ -236,8 +255,17 @@ TEST_F(toQOTest, ScanSupercellFCC) define_fcc_cell(ucell); toQO tqo("hydrogen", {"minimal-nodeless", "minimal-nodeless"}); tqo.unwrap_unitcell(&ucell); - tqo.build_nao(ucell.ntype, ucell.orbital_fn); - tqo.build_ao(ucell.ntype); + tqo.build_nao(ucell.ntype, + "./", + ucell.orbital_fn, + 0); + tqo.build_ao(ucell.ntype, + "./", + ucell.pseudo_fn, + {}, + GlobalV::qo_thr, + GlobalV::ofs_running, + 0); tqo.scan_supercell(); EXPECT_EQ(tqo.nR(), 13); } @@ -247,9 +275,18 @@ TEST_F(toQOTest, ScanSupercellSC1) define_sc_cell(ucell); toQO tqo("hydrogen", {"minimal-nodeless"}); tqo.unwrap_unitcell(&ucell); - tqo.build_nao(ucell.ntype, ucell.orbital_fn); + tqo.build_nao(ucell.ntype, + "./", + ucell.orbital_fn, + 0); GlobalV::qo_thr = 1e-6; - tqo.build_ao(ucell.ntype); + tqo.build_ao(ucell.ntype, + "./", + ucell.pseudo_fn, + {}, + GlobalV::qo_thr, + GlobalV::ofs_running, + 0); tqo.scan_supercell(); EXPECT_EQ(tqo.nR(), 19); // 3*3*3 - 8 (corner 111, -1-1-1, etc) } @@ -259,13 +296,22 @@ TEST_F(toQOTest, AllocateOvlpMinimal) define_fcc_cell(ucell); toQO tqo("hydrogen", {"minimal-nodeless", "minimal-nodeless"}); tqo.unwrap_unitcell(&ucell); - tqo.build_nao(ucell.ntype, ucell.orbital_fn); + tqo.build_nao(ucell.ntype, + "./", + ucell.orbital_fn, + 0); std::vector nmax = std::vector(ucell.ntype); for(int itype = 0; itype < ucell.ntype; itype++) { nmax[itype] = tqo.atom_database().principle_quantum_number[tqo.symbols()[itype]]; } - tqo.build_ao(ucell.ntype); + tqo.build_ao(ucell.ntype, + "./", + ucell.pseudo_fn, + {}, + GlobalV::qo_thr, + GlobalV::ofs_running, + 0); tqo.scan_supercell(); tqo.allocate_ovlp(true); tqo.allocate_ovlp(false); @@ -297,6 +343,7 @@ TEST_F(toQOTest, AllocateOvlpMinimal) TEST_F(toQOTest, Initialize) { define_fcc_cell(ucell); + GlobalV::qo_screening_coeff = {}; toQO tqo("hydrogen", {"minimal-nodeless", "minimal-nodeless"}); std::vector> kvecs_c; kvecs_c.push_back(ModuleBase::Vector3(0.0, 0.0, 0.0)); // Gamma point @@ -306,6 +353,7 @@ TEST_F(toQOTest, Initialize) TEST_F(toQOTest, CalculateOvlpR) { define_fcc_cell(ucell); + GlobalV::qo_screening_coeff = {}; GlobalV::qo_thr = 1e-10; toQO tqo("hydrogen", {"minimal-nodeless", "minimal-nodeless"}); std::vector> kvecs_c; @@ -338,6 +386,7 @@ TEST_F(toQOTest, CalculateOvlpR) TEST_F(toQOTest, CalculateSelfOvlpRMinimal) { define_fcc_cell(ucell); + GlobalV::qo_screening_coeff = {}; GlobalV::qo_thr = 1e-10; toQO tqo("hydrogen", {"minimal-nodeless", "minimal-nodeless"}); std::vector> kvecs_c; @@ -369,6 +418,7 @@ TEST_F(toQOTest, CalculateSelfOvlpKSymmetrical) { define_fcc_cell(ucell); GlobalV::qo_thr = 1e-10; + GlobalV::qo_screening_coeff = {}; toQO tqo("hydrogen", {"minimal-nodeless", "minimal-nodeless"}); std::vector> kvecs_c; kvecs_c.push_back(ModuleBase::Vector3(-0.25, -0.25, -0.25)); // pair 1 @@ -503,7 +553,13 @@ TEST_F(toQOTest, BuildHydrogenFull) toQO tqo("hydrogen", {"full", "full"}); tqo.unwrap_unitcell(&ucell); GlobalV::qo_thr = 1e-10; - tqo.build_ao(ucell.ntype); + tqo.build_ao(ucell.ntype, + "./", + ucell.pseudo_fn, + {}, + GlobalV::qo_thr, + GlobalV::ofs_running, + 0); EXPECT_EQ(tqo.p_ao()->nchi(), 9); // Si: 1s, 2s, 2p, 3s, 3p, 3d, C: 1s, 2s, 2p EXPECT_EQ(tqo.nchi(), 19); tqo.p_ao()->to_file("special_use_unittest"); @@ -513,6 +569,7 @@ TEST_F(toQOTest, CalculateSelfOvlpRFull) { define_fcc_cell(ucell); GlobalV::qo_thr = 1e-10; + GlobalV::qo_screening_coeff = {}; toQO tqo("hydrogen", {"full", "full"}); std::vector> kvecs_c; kvecs_c.push_back(ModuleBase::Vector3(0.0, 0.0, 0.0)); // Gamma point @@ -555,7 +612,13 @@ TEST_F(toQOTest, BuildPswfcPartial1) define_fcc_cell(ucell); toQO tqo("pswfc", {"s", "s"}); tqo.unwrap_unitcell(&ucell); - tqo.build_ao(ucell.ntype, ucell.pseudo_fn); + tqo.build_ao(ucell.ntype, + "./", + ucell.pseudo_fn, + GlobalV::qo_screening_coeff, + GlobalV::qo_thr, + GlobalV::ofs_running, + 0); EXPECT_EQ(tqo.p_ao()->nchi(), 5); // AO will always read and import all orbitals EXPECT_EQ(tqo.nchi(), 2); } @@ -565,7 +628,13 @@ TEST_F(toQOTest, BuildPswfcPartial2) define_fcc_cell(ucell); toQO tqo("pswfc", {"ps", "s"}); tqo.unwrap_unitcell(&ucell); - tqo.build_ao(ucell.ntype, ucell.pseudo_fn); + tqo.build_ao(ucell.ntype, + "./", + ucell.pseudo_fn, + GlobalV::qo_screening_coeff, + GlobalV::qo_thr, + GlobalV::ofs_running, + 0); EXPECT_EQ(tqo.p_ao()->nchi(), 5); // AO will always read and import all orbitals EXPECT_EQ(tqo.nchi(), 8); // the first element is Si, it has two p orbitals, so 3+3+1+1 } @@ -575,7 +644,13 @@ TEST_F(toQOTest, BuildPswfcPartial3) define_fcc_cell(ucell); toQO tqo("pswfc", {"all", "p"}); tqo.unwrap_unitcell(&ucell); - tqo.build_ao(ucell.ntype, ucell.pseudo_fn); + tqo.build_ao(ucell.ntype, + "./", + ucell.pseudo_fn, + GlobalV::qo_screening_coeff, + GlobalV::qo_thr, + GlobalV::ofs_running, + 0); EXPECT_EQ(tqo.p_ao()->nchi(), 5); // AO will always read and import all orbitals EXPECT_EQ(tqo.nchi(), 10); } @@ -583,9 +658,16 @@ TEST_F(toQOTest, BuildPswfcPartial3) TEST_F(toQOTest, BuildPswfcAll) { define_fcc_cell(ucell); + GlobalV::qo_thr = 1e-10; toQO tqo("pswfc", {"all", "all"}); tqo.unwrap_unitcell(&ucell); - tqo.build_ao(ucell.ntype, ucell.pseudo_fn); + tqo.build_ao(ucell.ntype, + "./", + ucell.pseudo_fn, + GlobalV::qo_screening_coeff, + GlobalV::qo_thr, + GlobalV::ofs_running, + 0); EXPECT_EQ(tqo.p_ao()->nchi(), 5); EXPECT_EQ(tqo.nchi(), 11); tqo.p_ao()->to_file("special_use_unittest"); @@ -596,10 +678,19 @@ TEST_F(toQOTest, ScanSupercellSC2) define_sc_cell(ucell); toQO tqo("pswfc", {"all"}); tqo.unwrap_unitcell(&ucell); - tqo.build_nao(ucell.ntype, ucell.orbital_fn); + tqo.build_nao(ucell.ntype, + "./", + ucell.orbital_fn, + 0); GlobalV::qo_screening_coeff[0] = 0.1; // use this to control the tailing of radial function GlobalV::qo_thr = 1e-6; - tqo.build_ao(ucell.ntype, ucell.pseudo_fn); // radius = 13.6 Bohr + tqo.build_ao(ucell.ntype, + "./", + ucell.pseudo_fn, + GlobalV::qo_screening_coeff, + GlobalV::qo_thr, + GlobalV::ofs_running, + 0); // radius = 13.6 Bohr tqo.scan_supercell(); EXPECT_EQ(tqo.nR(), 81); // 5*5*5 - 12(edge center) - 8*4(corner) } @@ -609,10 +700,19 @@ TEST_F(toQOTest, ScanSupercellSC3) define_sc_cell(ucell); toQO tqo("pswfc", {"all"}); tqo.unwrap_unitcell(&ucell); - tqo.build_nao(ucell.ntype, ucell.orbital_fn); + tqo.build_nao(ucell.ntype, + "./", + ucell.orbital_fn, + 0); GlobalV::qo_screening_coeff[0] = 0.25; // use this to control the tailing of radial function GlobalV::qo_thr = 1e-6; - tqo.build_ao(ucell.ntype, ucell.pseudo_fn); // radius = 13.6 Bohr + tqo.build_ao(ucell.ntype, + "./", + ucell.pseudo_fn, + GlobalV::qo_screening_coeff, + GlobalV::qo_thr, + GlobalV::ofs_running, + 0); // radius = 13.6 Bohr tqo.scan_supercell(); EXPECT_EQ(tqo.nR(), 57); // 5*5*5 - 12(edge center) - 8*(8-1)(corner) = 5*5*5 - 12(edge center) - 8*(2*2*2-1)(corner) GlobalV::qo_screening_coeff[0] = 0.1; @@ -623,10 +723,19 @@ TEST_F(toQOTest, ScanSupercellSC4) define_sc_cell(ucell); toQO tqo("pswfc", {"all"}); tqo.unwrap_unitcell(&ucell); - tqo.build_nao(ucell.ntype, ucell.orbital_fn); + tqo.build_nao(ucell.ntype, + "./", + ucell.orbital_fn, + 0); GlobalV::qo_screening_coeff[0] = 0.5; // use this to control the tailing of radial function GlobalV::qo_thr = 1e-6; - tqo.build_ao(ucell.ntype, ucell.pseudo_fn); // radius = 13.6 Bohr + tqo.build_ao(ucell.ntype, + "./", + ucell.pseudo_fn, + GlobalV::qo_screening_coeff, + GlobalV::qo_thr, + GlobalV::ofs_running, + 0); // radius = 13.6 Bohr tqo.scan_supercell(); EXPECT_EQ(tqo.nR(), 33); // 3*3*3 + 6(face) GlobalV::qo_screening_coeff[0] = 0.1; @@ -668,13 +777,13 @@ TEST_F(toQOTest, CalculateSelfOvlpRPswfc) } std::remove("Si_special_use_unittest.orb"); std::remove("C_special_use_unittest.orb"); - //tqo.write_ovlp(tqo.ovlp_R()[0], "QO_self_ovlp.dat"); } TEST_F(toQOTest, CalculateOvlpKGamma) { define_fcc_cell(ucell); GlobalV::qo_thr = 1e-10; + GlobalV::qo_screening_coeff = {}; toQO tqo("hydrogen", {"minimal-nodeless", "minimal-nodeless"}); std::vector> kvecs_c; kvecs_c.push_back(ModuleBase::Vector3(0.0, 0.0, 0.0)); // Gamma point @@ -692,6 +801,32 @@ TEST_F(toQOTest, CalculateOvlpKGamma) } } } + EXPECT_TRUE(all_real); +} + +TEST_F(toQOTest, CalculateOvlpKSlaterGamma) +{ + define_fcc_cell(ucell); + GlobalV::qo_thr = 1e-10; + GlobalV::qo_screening_coeff = {0.1}; + toQO tqo("hydrogen", {"energy-full", "energy-full"}); + std::vector> kvecs_c; + kvecs_c.push_back(ModuleBase::Vector3(0.0, 0.0, 0.0)); // Gamma point + tqo.initialize(&ucell, kvecs_c); + tqo.calculate_ovlp_k(0); + // all should be real numbers at Gamma point + bool all_real = true; + for(int i = 0; i < tqo.ovlp_k().size(); i++) + { + for(int j = 0; j < tqo.ovlp_k()[i].size(); j++) + { + if(tqo.ovlp_k()[i][j].imag() != 0.0) + { + all_real = false; + } + } + } + EXPECT_TRUE(all_real); } TEST_F(toQOTest, CalculateSelfOvlpKPswfcSymmetrical) @@ -792,6 +927,7 @@ TEST_F(toQOTest, CalculateHydrogenlike) { define_fcc_cell(ucell); GlobalV::qo_thr = 1e-10; + GlobalV::qo_screening_coeff = {}; toQO tqo("hydrogen", {"minimal-nodeless", "minimal-nodeless"}); std::vector> kvecs_c; kvecs_c.push_back(ModuleBase::Vector3(0.0, 0.0, 0.0)); // Gamma point diff --git a/source/module_io/test/write_input_test.cpp b/source/module_io/test/write_input_test.cpp index 0103c0a14b..72393774ee 100644 --- a/source/module_io/test/write_input_test.cpp +++ b/source/module_io/test/write_input_test.cpp @@ -393,7 +393,7 @@ TEST_F(write_input, Mixing7) EXPECT_THAT(output, testing::HasSubstr("mixing_angle -10 #angle mixing parameter for non-colinear calculations")); EXPECT_THAT(output, testing::HasSubstr("mixing_tau 0 #whether to mix tau in mGGA calculation")); EXPECT_THAT(output, testing::HasSubstr("mixing_dftu 0 #whether to mix locale in DFT+U calculation")); - EXPECT_THAT(output, testing::HasSubstr("mixing_restart 0 #which step to restart mixing during SCF")); + EXPECT_THAT(output, testing::HasSubstr("mixing_restart 0 #threshold to restart mixing during SCF")); EXPECT_THAT(output, testing::HasSubstr("mixing_dmr 0 #whether to mix real-space density matrix")); EXPECT_THAT(output, testing::HasSubstr("")); ifs.close(); diff --git a/source/module_io/to_qo.cpp b/source/module_io/to_qo.cpp index 2a17014ff3..5e448ae169 100644 --- a/source/module_io/to_qo.cpp +++ b/source/module_io/to_qo.cpp @@ -42,14 +42,23 @@ void toQO::initialize(UnitCell* p_ucell, */ // build the numerical atomic orbital basis // PARALLELIZATION STRATEGY: use RANK-0 to read in the files, then broadcast - build_nao(p_ucell_->ntype, p_ucell_->orbital_fn); + build_nao(p_ucell_->ntype, + GlobalV::global_orbital_dir, + p_ucell_->orbital_fn, + GlobalV::MY_RANK); // build another atomic orbital // PARALLELIZATION STRATEGY: only RANK-0 works #ifdef __MPI if(GlobalV::MY_RANK == 0) { #endif - build_ao(ntype_, p_ucell_->pseudo_fn); + build_ao(ntype_, + GlobalV::global_pseudo_dir, + p_ucell_->pseudo_fn, + GlobalV::qo_screening_coeff, + GlobalV::qo_thr, + GlobalV::ofs_running, + GlobalV::MY_RANK); // neighbor list search scan_supercell(); // build grids @@ -71,7 +80,10 @@ void toQO::initialize(UnitCell* p_ucell, #endif } -void toQO::build_nao(const int ntype, const std::string* const orbital_fn) +void toQO::build_nao(const int ntype, + const std::string orbital_dir, + const std::string* const orbital_fn, + const int rank) { // build the numerical atomic orbital basis ModuleBase::SphericalBesselTransformer sbt; @@ -82,11 +94,11 @@ void toQO::build_nao(const int ntype, const std::string* const orbital_fn) Parallel_Common::bcast_int(ntype_); #endif std::string* orbital_fn_ = new std::string[ntype_]; - if(GlobalV::MY_RANK == 0) + if(rank == 0) { for(int it = 0; it < ntype_; it++) { - orbital_fn_[it] = GlobalV::global_orbital_dir + orbital_fn[it]; + orbital_fn_[it] = orbital_dir + orbital_fn[it]; } } #ifdef __MPI @@ -108,7 +120,7 @@ void toQO::build_nao(const int ntype, const std::string* const orbital_fn) nphi_ += _nphi_it*na_[it]; } #ifdef __MPI - if(GlobalV::MY_RANK == 0) + if(rank == 0) { #endif printf("Build numerical atomic orbital basis done.\n"); @@ -126,10 +138,21 @@ bool toQO::orbital_filter(const int l, const std::string spec) else return false; } -void toQO::build_hydrogen(const int ntype, const double* const charges, const int* const nmax) +void toQO::build_hydrogen(const int ntype, + const double* const charges, + const bool slater_screening, + const int* const nmax, + const double qo_thr, + const int rank) { ao_ = std::unique_ptr(new RadialCollection); - ao_->build(ntype, charges, nmax, symbols_.data(), GlobalV::qo_thr, strategies_.data()); + ao_->build(ntype, + charges, + slater_screening, + nmax, + symbols_.data(), + qo_thr, + strategies_.data()); ModuleBase::SphericalBesselTransformer sbt; ao_->set_transformer(sbt); @@ -144,24 +167,30 @@ void toQO::build_hydrogen(const int ntype, const double* const charges, const in } #ifdef __MPI - if(GlobalV::MY_RANK == 0) + if(rank == 0) { #endif - printf("Build arbitrary atomic orbital basis done.\n"); + if(nchi_ > 0) printf("Build arbitrary atomic orbital basis done.\n"); + else ModuleBase::WARNING_QUIT("toQO::initialize", "Error: no atomic orbital is built."); #ifdef __MPI } #endif } -void toQO::build_pswfc(const int ntype, const std::string* const pspot_fn, const double* const screening_coeffs) +void toQO::build_pswfc(const int ntype, + const std::string pseudo_dir, + const std::string* const pspot_fn, + const double* const screening_coeffs, + const double qo_thr, + const int rank) { ao_ = std::unique_ptr(new RadialCollection); std::string* pspot_fn_ = new std::string[ntype_]; for(int it = 0; it < ntype; it++) { - pspot_fn_[it] = GlobalV::global_pseudo_dir + pspot_fn[it]; + pspot_fn_[it] = pseudo_dir + pspot_fn[it]; } - ao_->build(ntype, pspot_fn_, screening_coeffs, GlobalV::qo_thr); + ao_->build(ntype, pspot_fn_, screening_coeffs, qo_thr); ModuleBase::SphericalBesselTransformer sbt; ao_->set_transformer(sbt); @@ -176,7 +205,7 @@ void toQO::build_pswfc(const int ntype, const std::string* const pspot_fn, const } #ifdef __MPI - if(GlobalV::MY_RANK == 0) + if(rank == 0) { #endif printf("Build arbitrary atomic orbital basis done.\n"); @@ -185,21 +214,66 @@ void toQO::build_pswfc(const int ntype, const std::string* const pspot_fn, const #endif delete[] pspot_fn_; } - -void toQO::build_ao(const int ntype, const std::string* const pspot_fn) +/* +void toQO::build_szv(const int ntype) +{ + // build the numerical atomic orbital basis + ModuleBase::SphericalBesselTransformer sbt; + ao_ = std::unique_ptr(new RadialCollection); + // add GlobalV::global_orbital_dir ahead of orbital_fn + int ntype_ = ntype; + std::vector orbital_fn_(ntype_, "szv.orb"); + ao_->build(ntype_, orbital_fn_.data(), 'o'); + ao_->set_transformer(sbt); + for(int itype = 0; itype < ntype; itype++) + { + int _nchi_it = 0; + for(int l = 0; l <= ao_->lmax(itype); l++) + { + _nchi_it += (2*l+1)*ao_->nzeta(itype, l); + } + nchi_ += _nchi_it * na_[itype]; + } +} +*/ +void toQO::build_ao(const int ntype, + const std::string pseudo_dir, + const std::string* const pspot_fn, + const std::vector screening_coeffs, + const double qo_thr, + const std::ofstream& ofs_running, + const int rank) { if(qo_basis_ == "hydrogen") { - build_hydrogen(ntype_, charges_.data(), nmax_.data()); + bool with_slater_screening = std::find_if(screening_coeffs.begin(), screening_coeffs.end(), + [](double sc) { return sc > 1e-10; }) != screening_coeffs.end(); + build_hydrogen(ntype_, + charges_.data(), + with_slater_screening, + nmax_.data(), + qo_thr, + rank); } else if(qo_basis_ == "pswfc") { - build_pswfc(ntype_, pspot_fn, GlobalV::qo_screening_coeff.data()); + build_pswfc(ntype_, + pseudo_dir, + pspot_fn, + screening_coeffs.data(), + qo_thr, + rank); + } + /* + else if(qo_basis_ == "szv") + { + build_szv(ntype_); } + */ else { #ifdef __MPI - if(GlobalV::MY_RANK == 0) + if(rank == 0) { #endif // Not implemented error @@ -270,7 +344,7 @@ void toQO::calculate_ovlp_R(const int iR) for(int mj : mjs) { // TWO ATOMIC ORBITALS ARE SPECIFIED, THEN WE NEED TO CALCULATE THE OVERLAP IN SUPERCELL - ModuleBase::Vector3 rij = p_ucell_->atoms[it].tau[ia] - p_ucell_->atoms[jt].tau[ja]; + ModuleBase::Vector3 rij = p_ucell_->atoms[jt].tau[ja] - p_ucell_->atoms[it].tau[ia]; // there is waste here, but for easy to understand, I don't optimize it. ModuleBase::Vector3 R = supercells_[iR]; ModuleBase::Vector3 Rij; diff --git a/source/module_io/to_qo.h b/source/module_io/to_qo.h index c278877093..2e98dc6274 100644 --- a/source/module_io/to_qo.h +++ b/source/module_io/to_qo.h @@ -65,21 +65,42 @@ class toQO /// @brief build RadialCollection for numerical atomic orbitals /// @param ntype number of atom types /// @param orbital_fn filenames of numerical atomic orbitals - void build_nao(const int ntype, const std::string* const orbital_fn); + void build_nao(const int ntype, + const std::string orbital_dir, + const std::string* const orbital_fn, + const int rank); /// @brief build RadialCollection for atomic orbitals /// @param ntype number of atom types /// @param charges charges of atoms /// @param nmax maximum principle quantum number of atoms - void build_hydrogen(const int ntype, const double* const charges, const int* const nmax); + void build_hydrogen(const int ntype, + const double* const charges, + const bool slater_screening, + const int* const nmax, + const double qo_thr, + const int rank); + //function might be implemented in future + //void build_szv(const int ntype); /// @brief build RadialCollection for atomic orbitals /// @param ntype number of atom types /// @param pspot_fn filenames of pseudopotentials /// @param screening_coeffs screening coefficients of pseudopotentials, appears like a factor (exp[-s*r]) scaling the pswfc - void build_pswfc(const int ntype, const std::string* const pspot_fn, const double* const screening_coeffs); + void build_pswfc(const int ntype, + const std::string pseudo_dir, + const std::string* const pspot_fn, + const double* const screening_coeffs, + const double qo_thr, + const int rank); /// @brief build RadialCollection for atomic orbitals /// @param ntype number of atom types /// @param pspot_fn filenames of pseudopotentials, if use qo_basis = hydrogen, omit this parameter - void build_ao(const int ntype, const std::string* const pspot_fn = nullptr); + void build_ao(const int ntype, + const std::string pseudo_dir, + const std::string* const pspot_fn = nullptr, + const std::vector screening_coeffs = std::vector(), + const double qo_thr = 1e-10, + const std::ofstream& ofs = std::ofstream(), + const int rank = 0); /* * Main functions */ diff --git a/source/module_io/write_input.cpp b/source/module_io/write_input.cpp index 3970efb9e5..f33bab90c3 100644 --- a/source/module_io/write_input.cpp +++ b/source/module_io/write_input.cpp @@ -249,7 +249,7 @@ ModuleBase::GlobalFunc::OUTP(ofs, "out_bandgap", out_bandgap, "if true, print ou ModuleBase::GlobalFunc::OUTP(ofs, "mixing_type", mixing_mode, "plain; pulay; broyden"); ModuleBase::GlobalFunc::OUTP(ofs, "mixing_beta", mixing_beta, "mixing parameter: 0 means no new charge"); ModuleBase::GlobalFunc::OUTP(ofs, "mixing_ndim", mixing_ndim, "mixing dimension in pulay or broyden"); - ModuleBase::GlobalFunc::OUTP(ofs, "mixing_restart", mixing_restart, "which step to restart mixing during SCF"); + ModuleBase::GlobalFunc::OUTP(ofs, "mixing_restart", mixing_restart, "threshold to restart mixing during SCF"); ModuleBase::GlobalFunc::OUTP(ofs, "mixing_gg0", mixing_gg0, "mixing parameter in kerker"); ModuleBase::GlobalFunc::OUTP(ofs, "mixing_beta_mag", mixing_beta_mag, "mixing parameter for magnetic density"); ModuleBase::GlobalFunc::OUTP(ofs, "mixing_gg0_mag", mixing_gg0_mag, "mixing parameter in kerker"); diff --git a/source/module_relax/relax_old/ions_move_bfgs.cpp b/source/module_relax/relax_old/ions_move_bfgs.cpp index e44e4f9dc5..5d3e65b4b6 100644 --- a/source/module_relax/relax_old/ions_move_bfgs.cpp +++ b/source/module_relax/relax_old/ions_move_bfgs.cpp @@ -40,7 +40,17 @@ void Ions_Move_BFGS::start(UnitCell& ucell, const ModuleBase::matrix& force, con // istep must be set eariler. // use force to setup gradient. - Ions_Move_Basic::setup_gradient(ucell, force, this->pos, this->grad); + // Only the first step needs to generate the pos from ucell. + // In the following steps, the pos is updated by BFGS methods. + if (first_step) + { + Ions_Move_Basic::setup_gradient(ucell, force, this->pos, this->grad); + first_step = false; + } + else{ + std::vector pos_tmp(3 * ucell.nat); + Ions_Move_Basic::setup_gradient(ucell, force, pos_tmp.data(), this->grad); + } // use energy_in and istep to setup etot and etot_old. Ions_Move_Basic::setup_etot(energy_in, 0); // use gradient and etot and etot_old to check diff --git a/source/module_relax/relax_old/ions_move_bfgs.h b/source/module_relax/relax_old/ions_move_bfgs.h index 5338c43c3e..42a08a2593 100644 --- a/source/module_relax/relax_old/ions_move_bfgs.h +++ b/source/module_relax/relax_old/ions_move_bfgs.h @@ -17,6 +17,7 @@ class Ions_Move_BFGS : public BFGS_Basic bool init_done; void bfgs_routine(const double& lat0); void restart_bfgs(const double& lat0); + bool first_step=true; // If it is the first step of the relaxation. The pos is only generated from ucell in the first step, and in the following steps, the pos is generated from the previous step. }; #endif diff --git a/source/version.h b/source/version.h index f7de8c919e..9465a46987 100644 --- a/source/version.h +++ b/source/version.h @@ -1,3 +1,3 @@ #ifndef VERSION -#define VERSION "v3.5.3" +#define VERSION "v3.5.4" #endif diff --git a/tools/qo/README.md b/tools/qo/README.md index f3962e1d8c..d04b7052b1 100644 --- a/tools/qo/README.md +++ b/tools/qo/README.md @@ -1,101 +1,97 @@ -# Quasiatomic Orbital Analysis for ABACUS - -- Hamiltonian matrix transformation - -## Requirements - -- [NumPy](https://numpy.org/) -- [SciPy](https://www.scipy.org/) -- [Matplotlib](https://matplotlib.org/) (optional, for plotting) - -## Installation +*Use Markdown render to convert this document into directly human-readable version* +# ABACUS LCAO2QO Hands-on +## Introduction +This document is a hands-on tutorial for the ABACUS LCAO2QO tool. For basic information about QO (Quasiatomic Orbital) method, please refer to *Quasiatomic orbitals for ab initio tight-binding analysis* by *Qian et al.*: https://journals.aps.org/prb/abstract/10.1103/PhysRevB.78.245112. +### Basic ideas +QO is for projecting as many as possible eigenstates of a system onto a atom-centered basis set, unlike MLWFs (Maximally Localized Wannier Functions) that would produce bond-like Wannier functions for some systems. The QO is defined from an optimization problem that maximal the norm of the projection of atom-centered orbitals onto the super-plane spanned by eigenstates and constructed unoccupied states: +$\max\mathcal{L}:=\max \sum_{\mu}{\left\|\left(\sum_{n\mathbf{k}}{\hat{P}_{\psi_{n\mathbf{k}}}}+\sum_{m\mathbf{k}}{\hat{P}_{c_{m\mathbf{k}}}}\right)|A_{\mu}\rangle\right\|^2}$ +where $\hat{P}_{\psi_{n\mathbf{k}}}$ and $\hat{P}_{c_{m\mathbf{k}}}$ are the projection operators for occupied and unoccupied states, respectively. $A_{\mu}$ is the atom-centered orbital. +The former part of projection is fixed once the occupied states are determined, and the latter part is determined by selection of $Nq-Nb$ from $\sum_{m\mathbf{k}}{\hat{P}_{c_{m\mathbf{k}}}}|A_{\mu}\rangle$ states, with total number of $Nq$ and number of occupied states $Nb$. +### Formulation +See Github issue #3640 for details: https://github.com/deepmodeling/abacus-develop/issues/3640 +## Prerequisites +- ABACUS +- Python3 +- Numpy +- Scipy +- Matplotlib (optional) +## Usage +### Basis set selection +ABACUS provides various kinds of atom-centered orbitals construction, such as hydrogen-like orbitals, exponentially decay scaled pseudowavefunctions, single zeta (not implemented, develop up to user actual needs). Select orbitals that fit most to the system you are interested in. +### Input file of ABACUS +To run QO analysis, in ABACUS INPUT file, following keywords are needed to be set: +```bash +qo_switch 1 +qo_basis hydrogen +qo_strategy energy-valence +qo_thr 1.0e-10 +#qo_screening_coeff 0.1 +``` +- `qo_switch` is a switch to turn on QO analysis. +- `qo_basis` is the basis set selection, parameters above define the hydrogen-like orbital construction +- `qo_strategy` is the strategy to set what orbitals needed to construct. For hydrogen-like orbitals, for example if element is Cu (1s2 2s2 2p6 3s2 3p6 3d10 4s1), `energy-valence` corresponds to construct only 3p 3d and 4s orbitals, `energy-full` corresponds to construct all orbitals, `minimal-nodeless` corresponds to construct 1s, 2p, 3d and 4f orbitals, `minimal-valence` corresponds to construct 4s 4p 4d 4f orbitals. For more information about this keyword, please refer to ABACUS manual: https://abacus.deepmodeling.com/en/latest/advanced/input_files/input-main.html#qo-strategy +- `qo_thr` is the threshold to determine the realspace tailing of orbitals, which is used to determine the cutoff radius of orbitals. Default value is 1.0e-10. +- `qo_screening_coeff`, for `qo_basis hydrogen` it controls whether use Slater screening coefficients to include many-electron effect into the orbitals. Any float number as value will turn on this feature. Default value is off. -This code runs like scripts, you can also install this by: +For some systems, pseudowavefunction works better than hydrogen-like orbitals (but not all pseudopotential has pseudowavefunction in file, check your pseudopotential file before run). To use pseudowavefunction as `qo_basis`, set as following: ```bash -python3 setup.py install +qo_switch 1 +qo_basis pswfc +qo_strategy all +qo_thr 1.0e-10 +qo_screening_coeff 0.5 ``` -, then you can use it like a package: -```python -import abacus2qo +For `qo_basis pswfc`, `qo_screening_coeff` is necessary and for Si single crystal the optimal value is tested to be `0.5`. Test on this parameter should be done for different systems. A larger value will result in more localized orbitals. +### Output file of ABACUS +Once the calculation is done, the output file will contain the QO analysis results. The QO analysis results are stored in `OUT.${suffix}` directory, where `${suffix}` is set in `INPUT` and default is `ABACUS`. You will find following files: +```bash +data-0-H ... +data-0-S ... ``` -If you would like to change code and run it, you can also install it by: +These files are Hamiltonian and overlap matrix of numerical atomic orbitals basis sets. For each kpoint a pair of these files will be generated, and different spin channels are defined as different kpoints. ```bash -python3 setup.py develop +LOWF_K_1.txt ... ``` +These files are eigenstates represented by numerical atomic orbital basis sets. For each kpoint a file will be generated. +```bash +QO_ovlp_0.dat ... +``` +This file contains information of $\langle A_\mu(\mathbf{k})|\phi_\nu(\mathbf{k})\rangle$, file name indiced by kpoint. +```bash +QO_supercells.dat +``` +This file stores information of supercell considered when calculating two-center-integrals, will be usefull in kpoint extrapolation. +```bash +kpoints +istate.info +running_scf.log +... +``` +Basic output files of ABACUS. +### Run `postprocess.py` +After the calculation is done, run `postprocess.py` to generate QO analysis results. The script will read the output files of ABACUS and generate the QO analysis results. Following defines an example of `postprocess.py`: +```python +if __name__ == "__main__": -## Usage - -1. Perform ABACUS `basis_type lcao` calculation with additional keywords: - ```Plain text - qo_switch 1 - qo_basis pswfc - ``` - , this will tell ABACUS to extract pseudowavefunction from pseudopotential file. To use this without any errors, please make sure there IS pseudowavefunction in pseudopotential file. - One can also add `qo_screening_coeff` keyword in `INPUT` to tune the behavior of pseudowavefunction (default is set to 0.1): - ```Plain text - qo_screening_coeff 0.1 - ``` - `qo_screening_coeff` always to be `0.1`, but whatever value you like, it must be a positive number to ensure a proper behavior at least at infinity, i.e. $\lim_{|\mathbf{r}|\rightarrow\infty}\tilde{\psi}(|\mathbf{r}|)\rightarrow0$. - There is also a parameter `qo_thr` controlling the tailing of pseudowavefunction, which defines as the convergence threshold of norm of pseudowavefunction when realspace cutoff increasing. To change this parameter, just add it in `INPUT`: - ```Plain text - qo_thr 1e-6 - ``` - There is also one *very experimental, very numerically instable* feature: - ```Plain text - qo_switch 1 - qo_basis hydrogen - ``` - To use hydrogen-like radial function as projector. However, because the charges cannot be set very physically presently (unless let user set them), ABACUS will read atomic charge from pseudopotential file, which may cause unphysical shrink or expansion. - For example see `examples`, find `INPUT`, `STRU` and `KPT`. -2. Copy output Hamiltonian and overlap matrices $H(\mathbf{k})$ and $S(\mathbf{k})$ files (`data-[i]-H` and `data-[i]-S`), the overlap matrices between AO and numerical atomic orbitals $S^{\chi\phi}(\mathbf{k})$ `QO_ovlp_[i].dat`, converged wavefunction `LOWF_K_[i].txt`, kpoints information summary `kpoints` from ABACUS output `OUT.*` directory to the path you like. Or in another way, just assign the `path` in `source/main.py` to the path of `OUT.*` directory. -3. Then you specify the path, number of kpoints, number of states want to introduce and the lattice vector $\mathbf{R}$ (presently it is the one set to $(0, 0, 0)$) in `source/main.py`. - ```python - import source.components.driver as driver - import numpy as np - - if __name__ == "__main__": - # example of QO - path = "./tests/integrate/220_NO_KP_QO/OUT.ABACUS" - # fully non-reduced kpoints - nkpts = 125 - # band range to reproduce - band_range = (0, 13) - - # create driver - d_ = driver.toQO_Driver() - # initialize driver - d_.initialize(path, nkpts, "scf", band_range) - # expand space spanned by selected bands - d_.space_expansion() - # calculate H(R) and S(R) in R-space in QO representation - HRs, SRs = d_.reproduce_hamiltonian(Rs=[ - np.array([0, 0, 0]) - ]) - ``` - `HRs` and `SRs` are list of `numpy.ndarray`, so it is convinient to use `numpy.save` to save them. - You can also find a file `QO_supercells.dat` in `OUT.${suffix}/`, which contains the information of supercell. You can use it as a reference to construct `Rs`. -4. run it! - -## Formulation info - -QO is simple in formulation, its central idea is to maximize the norm of QO, which is the projection of AO in extended space. The space is spanned by the selected bands concatenated with the pseudowavefunction. - -First project pseudowavefunction with NAO, which means add projection operator of NAO at left side of AO: -$\hat{P}|A_\mu(\mathbf{k})\rangle$ -, where $\mathbf{k}$ specifies the kpoint, $\mu$ distinguishes the AO, it is actually a index describing which type, which atom, which $(n, l, m)$. Projection operator $\hat{P}$ can be represented as: -$\hat{P} = \sum_{\alpha\beta}{|\phi_{\alpha}(\mathbf{k})\rangle \mathbf{S}^{-1}_{\alpha\beta}\langle\phi_{\beta}(\mathbf{k})|A_{\mu}(\mathbf{k})\rangle}$ -The matrix element $\mathbf{S}^{-1}_{\alpha\beta}$, hopefully denoted as it is here without any misleading, is the matrix element of inverse overlap matrix between NAOs. + # False to run unittest, True to run production + run_type = "production" # can be "unittest" or "production" + ndim = 7 # dimension of Monkhorst-Pack grid + qo_basis = "hydrogen" + qo_strategy = "minimal-valence" + path = f"./scf/{qo_basis}/{qo_strategy}/OUT.ABACUS-{ndim}{ndim}{ndim}/" + #fpic = f"{qo_basis}-{qo_strategy}-{ndim}{ndim}{ndim}.png" + # only for production + eig_thr = 1e-10 + ib_min, ib_max = 0, 4 -Similarly, we project AO into selected eigenstates spanned space: -$\hat{P}_{\psi}|A_\mu(\mathbf{k})\rangle$ = $\sum_{n}{|\psi_{\alpha}(\mathbf{k})\rangle\langle\psi_{\beta}(\mathbf{k})|}A_{\mu}(\mathbf{k})\rangle$ -, one can quickly rewrite it as: -$\begin{cases} - |A_{\mu}(\mathbf{k})\rangle \rightarrow \sum_{\alpha \beta}{\left( \mathbf{S}_{\alpha \beta}^{-1}\langle \phi _{\beta}(\mathbf{k})|A_{\mu}(\mathbf{k})\rangle \right) |\phi _{\alpha}(\mathbf{k})\rangle}\\ - |A_{\mu}^{\parallel}(\mathbf{k})\rangle \rightarrow \sum_{\alpha}{\left( \sum_{\beta}{\langle \phi _{\beta}(\mathbf{k})|A_{\mu}(\mathbf{k})\rangle \sum_n{c_{n\alpha}\left( \mathbf{k} \right) c_{n\beta}^{\dagger}\left( \mathbf{k} \right)}} \right) |\phi _{\alpha}(\mathbf{k})\rangle}\\ -\end{cases}$ + if run_type == "production": + nkpts = ndim*ndim*ndim + hqos_k, sqos_k, kvecs_d = production(path, nkpts, ib_min, ib_max, eig_thr, error_estimation_with_lcao=True) -Therefore the virtual states introduced by AO, which definitely is NAOs linearly combined in some manner, is: -$|A_{\mu}^{\bot}(\mathbf{k})\rangle =|A_{\mu}(\mathbf{k})\rangle -|A_{\mu}^{\parallel}(\mathbf{k})\rangle$ -Then canonical orhogonalization is performed, get only parts of eigenstates, denoted as $\{c_m(\mathbf{k})\}$. Then merge spaces spanned by $\{c_m(\mathbf{k})\}$ and $\{|\psi_{n}(\mathbf{k})\rangle\}$, denote as $\{\varphi_i(\mathbf{k})\}$. -Then the QO is: -$|Q_{\mu}(\mathbf{k})\rangle =\sum_i{|\varphi _i(\mathbf{k})\rangle \langle \varphi _i(\mathbf{k})|A_{\mu}(\mathbf{k})\rangle}$ \ No newline at end of file + elif run_type == "unittest": + unittest.main() +``` +The `production` function will read the output files of ABACUS and generate the QO analysis results. The `path` is the path to the output files of ABACUS. `nkpts` is the number of kpoints. `ib_min` and `ib_max` are the range of orbitals to be considered. `eig_thr` is the threshold to determine the eigenstates to be considered. `error_estimation_with_lcao` is a switch to turn on error estimation with LCAO. If set to `True`, the script will calculate the error estimation with LCAO. Collect information output by `production(...)` function. The `hqos_k` and `sqos_k` are the Hamiltonian and overlap matrix of QO basis sets. The `kvecs_d` is the direct coordinates of kpoints. +QO postprocess is programmed with unittest, therefore to run the script, set `run_type` to `production` and run the script (but mostly you don't need this). +## Acknowledgement +ABACUS LCAO2QO module uses two-center-integrator module refactored by @jinzx10, who also helps in code debugging. @mohanchen, @dyzheng, @WHUweiqingzhou and @QG-phys provide suggestions on code maintainance and participate discussion in technical aspects. \ No newline at end of file diff --git a/tools/qo/abacus2qo/__init__.py b/tools/qo/abacus2qo/__init__.py deleted file mode 100644 index e69de29bb2..0000000000 diff --git a/tools/qo/abacus2qo/components/__init__.py b/tools/qo/abacus2qo/components/__init__.py deleted file mode 100644 index e69de29bb2..0000000000 diff --git a/tools/qo/abacus2qo/components/basis_filter.py b/tools/qo/abacus2qo/components/basis_filter.py deleted file mode 100644 index 50e23b1840..0000000000 --- a/tools/qo/abacus2qo/components/basis_filter.py +++ /dev/null @@ -1,48 +0,0 @@ -import numpy as np -import scipy.linalg as la -class toQO_BasisFilter: - """Balance between flexible and accurate basis selection, filter basis according to practical band structure - matrix should be k-point resoluted, i.e., matrix[ik] is the matrix for kpoint ik - Raises: - ValueError: if band range is not properly configured - """ - band_range: tuple # min and max band index - overlap_matrix_filter_threshold: float # filter thresold for overlap matrix, only larger than this threshold will basis retain - denstiy_matrix_filter_threshold: float # filter thresold for density matrix, only larger than this threshold will basis retain - - def __init__(self) -> None: - pass - - def set_overlap_filter(self, overlap_matrix_filter_threshold: float) -> None: - """set overlap matrix filter threshold - - Args: - overlap_matrix_filter_threshold (float): overlap matrix filter threshold - """ - self.overlap_matrix_filter_threshold = overlap_matrix_filter_threshold - - def filter_via_overlap_matrix(self, sk: np.ndarray, saok: np.ndarray) -> np.ndarray: - """calculate AO overlap matrix in NAO representation. Those diagonal elements that are failed to be larger than threshold will be filtered out - because they cannot be represented by NAOs - - Args: - nkpts (int): number of kpoints - sk (np.ndarray): overlap matrix in k-space - saok (np.ndarray): AO-NAO overlap matrix - - Returns: - deleted_basis_indices (list): list of list, each list contains indices of selected basis for each kpoint - """ - selected_basis_indices = [] - _saok = np.zeros_like(saok) - - _x = la.solve(sk, saok.conj().T) - _ovlp = saok @ _x - - for i in range(_ovlp.shape[0]): - if _ovlp[i, i] > self.overlap_matrix_filter_threshold: - selected_basis_indices.append(i) - # filter out basis - _saok = saok[selected_basis_indices, :] - - return selected_basis_indices, _saok diff --git a/tools/qo/abacus2qo/components/calculator.py b/tools/qo/abacus2qo/components/calculator.py deleted file mode 100644 index be9ffb0bfd..0000000000 --- a/tools/qo/abacus2qo/components/calculator.py +++ /dev/null @@ -1,190 +0,0 @@ -import numpy as np -import scipy.linalg as la -import abacus2qo.components.safe_guard as sg - -class toQO_Calculator: - """python-end the Quasiatomic orbital (QO) analysis - """ - sg_: sg.toQO_SafeGuard - - def __init__(self) -> None: - self.sg_ = sg.toQO_SafeGuard() - - def folding_HR(self, matrices_R: list, supercells: list, kpoint: np.ndarray) -> np.ndarray: - """calculate fold matrix from R to k-space - - Args: - matrices_R (list): matrices in R-space - supercells (list): corresponding supercells direct vectors n1n2n3 - kpoint (np.ndarray): one specific kpoint direct coordinates - - Returns: - np.ndarray: matrix in k-space - """ - matrix_k = np.zeros_like(matrices_R[0]) - for iR in range(len(supercells)): - arg = np.exp(1j * supercells[iR] @ kpoint * 2 * np.pi) - matrix_k += arg * matrices_R[iR] - return matrix_k - - def unfolding_Hk(self, matrices_k: list, kpoints: list, supercell: np.ndarray) -> np.ndarray: - """calculate unfold matrix from k to R-space - - Args: - matrices_k (list): matrices in k-space - kpoints (list): corresponding equivalent kpoints, arranged in the same order as matrices_k - supercell (np.ndarray): one specific supercell direct coordinates - - Returns: - np.ndarray: matrix in R-space - """ - - matrix_R = np.zeros_like(matrices_k[0]) - for ik in range(len(kpoints)): - kpoint = kpoints[ik] - arg = np.exp(-1j * kpoint @ supercell * 2 * np.pi) - matrix_R += arg * matrices_k[ik]/np.sqrt(len(kpoints)) - return matrix_R - - def projto_nao(self, sk: np.ndarray, saok: np.ndarray) -> np.ndarray: - """calculate anything in nao representation with the nao projector - - Args: - sk (np.ndarray): nao overlap matrix in k-space - saok (np.ndarray): ao-nao overlap matrix in k-space - - Returns: - np.ndarray: anything in nao representation - """ - print("Calculate Atomic Orbital (AO) in NAO representation.") - psi_chi = np.linalg.solve(sk, saok.conj().T) - return psi_chi - - def projto_eigstate(self, psi_lcao: np.ndarray, saok: np.ndarray) -> np.ndarray: - """calculate any states projected onto eigenstates of Hamiltonian in k-space represented by NAO - - Args: - psi_lcao (np.ndarray): states in LCAO representation - saok (np.ndarray): ao-nao overlap matrix in k-space - - Returns: - np.ndarray: states projected onto eigenstates of Hamiltonian in k-space represented by NAO - """ - psi_chi_para = psi_lcao @ psi_lcao.conj().T @ saok.conj().T - return psi_chi_para - - def canonical_orthogonalization(self, psi_chi_orth: np.ndarray, sk: np.ndarray, m: int) -> np.ndarray: - """this is a general canonical orthogonalization procedure for any subspace - - Args: - psi_chi_orth (np.ndarray): states in one certain representation, here is nao - sk (np.ndarray): basis function overlap matrix - m (int): number of states retained - - Returns: - np.ndarray: states orthogonalized - """ - print("Perform canonical orthogonalization for newly appended subspace from AO introduction.") - wk = psi_chi_orth.conj().T @ sk @ psi_chi_orth - yk_k, vk_k = la.eigh(wk) - print("Eigenvalues of W(k) are: \n", yk_k) - # truncate m largest eigenvalues and eigenvectors - yk_k = yk_k[-m:] - vk_k = vk_k[:, -m:] - print("Get ", m, " largest eigenvalues and eigenvectors. Selected eigenvalues are: \n", yk_k) - # calculate psi_complem - yk_k = np.diag([np.sqrt(1/yk_k[i]) for i in range(m)]) - psi_complem = psi_chi_orth @ vk_k @ yk_k - print("Check orthogonalities between states in extended space.") - self.sg_.check_eigenstates(psi_complem, sk) - return psi_complem - - def merge_space(self, psi_lcao: np.ndarray, psi_complem: np.ndarray, hk: np.ndarray, sk: np.ndarray) -> np.ndarray: - """calculate psi_exten - - psi_exten = [psi_lcao, psi_complem] - - Args: - psi_lcao (np.ndarray): states in LCAO representation - psi_complem (np.ndarray): states in complementary representation - hk (np.ndarray): Hamiltonian in k-space - sk (np.ndarray): overlap in k-space - - Returns: - np.ndarray: extended states in k-space - """ - print("Combine occupied states and constructed virtual states to get extended wavefunction" - +"\n(empty states are appended)." - +"\nShape of occupied states is: ", psi_lcao.shape, "\nShape of empty states is: ", psi_complem.shape) - psi_exten = np.concatenate((psi_lcao, psi_complem), axis=1) - print("Shape of extended states is: ", psi_exten.shape) - # check orthogonality - print("Check orthogonality of psi_exten.") - sk_exten = psi_exten.conj().T @ sk @ psi_exten - # if sk_exten is not identity matrix? - error = self.sg_.check_identity(sk_exten) - while error > 1e-6: - print("Error of sk_exten not being identity is: ", error, ", unacceptable.") - exit() - - # if hk_exten is not diagonal in supspace psi_lcao? - hk_exten = psi_exten.conj().T @ hk @ psi_exten - self.sg_.check_diagonal(hk_exten[:psi_lcao.shape[1], :psi_lcao.shape[1]]) - # get extended energy spectrum - print("Get extended energy spectrum.") - eigvals_exten, eigvecs_exten = la.eigh(hk_exten, sk_exten) - print("Eigenvalues of H_exten(k) are: \n", eigvals_exten) - eigvals, eigvecs = la.eigh(hk, sk) - print("Comparatively eigenvalues of H(k) are: \n", eigvals) - return psi_exten - - def calculate_qo(self, saok: np.ndarray, psi_exten: np.ndarray, sk: np.ndarray) -> np.ndarray: - """calculate qo represented by NAO in kspace, return with nrows = nchi and ncols = nphi - - Args: - saok (np.ndarray): ao-nao overlap matrix in k-space - psi_exten (np.ndarray): extended states in k-space - sk (np.ndarray): overlap in k-space - - Returns: - np.ndarray: qo represented by NAO in kspace - """ - print("Calculate QO.") - qo = psi_exten @ psi_exten.conj().T @ saok.conj().T - # this saok is overlap between AO and NAO, line is AO, column is NAO - # then normalize qo - for i in range(qo.shape[1]): - qo[:, i] = qo[:, i] / np.sqrt(qo[:, i].conj().T @ sk @ qo[:, i]) - print("QO Normalization: after, norm of the ", i, "-th QO is: %10.8f"%((qo[:, i].conj().T @ sk @ qo[:, i]).real)) - return qo - - def calculate_hqok(self, qo: np.ndarray, hk: np.ndarray, sk: np.ndarray) -> np.ndarray: - """calculate hqok - - Args: - qo (np.ndarray): QO in k-space - hk (np.ndarray): Hamiltonian represented by QO in k-space - sk (np.ndarray): QO overlap in k-space - - Returns: - np.ndarray: hqok - """ - print("Calculate hamiltonian matrix in QO basis in k-space.") - hqok = qo.conj().T @ hk @ qo - sqok = qo.conj().T @ sk @ qo # this is overlap matrix in QO basis - - skip_diagonalization = False - eigval_s, eigvec_s = la.eigh(sqok) - print("Eigenvalues of overlap of in basis Sqo(k) are: \n", eigval_s) - if abs(eigval_s[0]) < 1e-6: - print("WARNING: BAD CONDITION NUMBER OF OVERLAP MATRIX IN QO BASIS, THIS WILL CAUSE PROBLEMS IN DIAGONALIZATION.") - return np.zeros_like(hqok), np.zeros_like(sqok) - - else: - eigvals_qo, eigvecs_qo = la.eigh(hqok, sqok) - print("Eigenvalues of Hamiltonian in QO basis Hqo(k) are: \n", eigvals_qo) - - eigvals_nao, eigvecs_nao = la.eigh(hk, sk) - print("Eigenvalues of Hamiltonian in NAO basis H(k) are: \n", eigvals_nao) - - return hqok, sqok diff --git a/tools/qo/abacus2qo/components/data_container.py b/tools/qo/abacus2qo/components/data_container.py deleted file mode 100644 index 2cf1dcde1c..0000000000 --- a/tools/qo/abacus2qo/components/data_container.py +++ /dev/null @@ -1,34 +0,0 @@ -import numpy as np - -class toQO_DataContainer: - """the one who really stores data, works as a repo - """ - # hamiltonians - hk: list # hamiltonian matrix in k-space - hqok: np.ndarray # hamiltonian matrix in QO basis in k-space - sqok: np.ndarray # overlap matrix in QO basis in k-space - # homogeneous overlaps - sk: list # overlap matrix in k-space - wk: list # overlap matrix in QO basis in k-space, list of np.ndarray - # heterogeneous overlaps/transformations - saok: list # AO-NAO overlap matrix in k-space, list of np.ndarray - # data - kpoints: list # kpoints direct coordinates - equivalent_kpoints: list # weights of kpoints - energies: list # energies of wavefunctions - # dimensions - nkpts: int # number of kpoints - nbands: int # number of bands - nchi: int # number of AOs - nphi: int # number of NAOs - # wavefunction in nao representations - psi_lcao: list # wavefunction represented by numerical atomic orbitals - psi_chi: list # ao represented by numerical atomic orbitals - psi_chi_para: list # parallel component of QO wavefunction represented by NAOs - psi_chi_orth: list # orthogonal component of QO wavefunction represented by NAOs - psi_complem: list # complementary component of QO wavefunction represented by NAOs - psi_exten: list # extended component of QO wavefunction represented by NAOs - psi_qo: list # QO represented by NAOs - - def __init__(self) -> None: - pass diff --git a/tools/qo/abacus2qo/components/data_manager.py b/tools/qo/abacus2qo/components/data_manager.py deleted file mode 100644 index f5e99553d1..0000000000 --- a/tools/qo/abacus2qo/components/data_manager.py +++ /dev/null @@ -1,163 +0,0 @@ -import numpy as np -import abacus2qo.components.data_container as dc -import abacus2qo.tools.hamiltonian as ham -import abacus2qo.tools.wavefunction as wf -import abacus2qo.tools.qo_ovlp as qov -import abacus2qo.tools.kpoints as kpt -import abacus2qo.components.safe_guard as sg -from scipy.linalg import eigh - -class toQO_DataManager: - """Manage memory storing the data - """ - data: dc.toQO_DataContainer - sg_: sg.toQO_SafeGuard - - def __init__(self) -> None: - self.data = dc.toQO_DataContainer() # Data manager borns with its data - self.sg_ = sg.toQO_SafeGuard() - - def kpoints_eq(self, kpt1: np.ndarray, kpt2: np.ndarray) -> bool: - """check if two kpoints are equal""" - return np.linalg.norm(kpt1 - kpt2) < 1e-6 - - def align_according_to_kpoints_ref(self, to_align: list, kpts: list, kpts_ref: list) -> list: - """align kpoints according to kpoints file""" - indexing = [] - # first check if number of kpoints are the same - if len(kpts) != len(kpts_ref): - raise ValueError("kpoints read from files are not the same") - nkpts = len(kpts_ref) - for ik in range(nkpts): - _kpt = kpts_ref[ik] - for jk in range(nkpts): - if self.kpoints_eq(_kpt, kpts_ref[jk]): - indexing.append(jk) - break - # indexing has the meaning: for the first kpoint in kpts_ref, the data would be found in the indexing[0]th position of to_align - # align - aligned = [] - for ik in range(nkpts): - aligned.append(to_align[indexing[ik]]) - - return aligned - - def check_kpoints_alignment(self, kpts1: list, kpts2: list) -> bool: - """check if two kpoints lists are aligned""" - nkpts = len(kpts1) - for ik in range(nkpts): - if self.kpoints_eq(kpts1[ik], kpts2[ik]): - return False - return True - - def read(self, nkpts: int, calculation: str, path: str, band_range: tuple) -> None: - """read data from files without any changes - - Args: - - nkpts (int): number of kpoints - path (str): path to the data - band_range (tuple): min and max band index - - Raises: - - ValueError: None in this level of function - """ - print(""" ---------------------------------------------------------------------------------------------------------- -$$$$$$$\ $$$$$$\ $$$$$$$$\ $$$$$$\ $$$$$$\ $$\ $$\ $$$$$$$\ $$$$$$\ $$$$$$$\ $$$$$$$$\ -$$ __$$\ $$ __$$\\__$$ __|$$ __$$\ \_$$ _|$$$\ $$$ |$$ __$$\ $$ __$$\ $$ __$$\\__$$ __| -$$ | $$ |$$ / $$ | $$ | $$ / $$ | $$ | $$$$\ $$$$ |$$ | $$ |$$ / $$ |$$ | $$ | $$ | -$$ | $$ |$$$$$$$$ | $$ | $$$$$$$$ | $$ | $$\$$\$$ $$ |$$$$$$$ |$$ | $$ |$$$$$$$ | $$ | -$$ | $$ |$$ __$$ | $$ | $$ __$$ | $$ | $$ \$$$ $$ |$$ ____/ $$ | $$ |$$ __$$< $$ | -$$ | $$ |$$ | $$ | $$ | $$ | $$ | $$ | $$ |\$ /$$ |$$ | $$ | $$ |$$ | $$ | $$ | -$$$$$$$ |$$ | $$ | $$ | $$ | $$ | $$$$$$\ $$ | \_/ $$ |$$ | $$$$$$ |$$ | $$ | $$ | -\_______/ \__| \__| \__| \__| \__| \______|\__| \__|\__| \______/ \__| \__| \__| ---------------------------------------------------------------------------------------------------------- -In this step, read data dumped from ABACUS. Please ensure you have enough precision on -printing several critical quantities, such as H(k), S(k), wavefunction and Sao(k).\n""") - # int - self.data.nkpts = nkpts - # list of np.ndarray(2d), list of np.ndarray(2d), int - self.data.hk, self.data.sk, self.data.nphi = ham.parse(self.data.nkpts, path) - # list of np.ndarray(2d), list of np.ndarray(nkpts, 3) - self.data.saok, saok_kpts = qov.parse(self.data.nkpts, path) - # list of np.ndarray(2d), list of np.ndarray(nkpts, 3), list of np.ndarray(1d) - self.data.psi_lcao, psi_kpts, self.data.energies = wf.parse(self.data.nkpts, path) - # list of np.ndarray(nkpts, 3), list of list of np.ndarray(nkpts, 3) - kpts, self.data.equivalent_kpoints = kpt.parse(path) - # presently there may be at least two sets of kpoints in ABACUS - # Hamiltonian and Overlap matrix -> kpoints file, kvec_d - # Overlap AO with NAO -> kpoints file, kvec_d - # wavefunction -> itself, kvec_c - # align all according to kpoints file - # if not self.check_kpoints_alignment(kpts, psi_kpts): - # self.data.psi_lcao = self.align_according_to_kpoints_ref(self.data.psi_lcao, psi_kpts, kpts) - # self.data.energies = self.align_according_to_kpoints_ref(self.data.energies, psi_kpts, kpts) - # if not self.check_kpoints_alignment(kpts, saok_kpts): - # self.data.saok = self.align_according_to_kpoints_ref(self.data.saok, saok_kpts, kpts) - self.data.kpoints = kpts - print(""" ------------------------------------------------------------------------------------------------ -$$$$$$$\ $$$$$$\ $$$$$$$$\ $$$$$$\ $$$$$$\ $$\ $$\ $$$$$$$$\ $$$$$$\ $$\ $$\ -$$ __$$\ $$ __$$\\__$$ __|$$ __$$\ $$ __$$\ $$ | $$ |$$ _____|$$ __$$\ $$ | $$ | -$$ | $$ |$$ / $$ | $$ | $$ / $$ | $$ / \__|$$ | $$ |$$ | $$ / \__|$$ |$$ / -$$ | $$ |$$$$$$$$ | $$ | $$$$$$$$ | $$ | $$$$$$$$ |$$$$$\ $$ | $$$$$ / -$$ | $$ |$$ __$$ | $$ | $$ __$$ | $$ | $$ __$$ |$$ __| $$ | $$ $$< -$$ | $$ |$$ | $$ | $$ | $$ | $$ | $$ | $$\ $$ | $$ |$$ | $$ | $$\ $$ |\$$\ -$$$$$$$ |$$ | $$ | $$ | $$ | $$ | \$$$$$$ |$$ | $$ |$$$$$$$$\ \$$$$$$ |$$ | \$$\ -\_______/ \__| \__| \__| \__| \__| \______/ \__| \__|\________| \______/ \__| \__| ------------------------------------------------------------------------------------------------ -Due to interface between C++ and Python will lose precision, in this step, -orthonormalities between eigenstates are needed to check.""") - for ik in range(self.data.nkpts): - print("-"*50, "\nFor k-point No.", ik) - error = self.sg_.check_eigenstates(self.data.psi_lcao[ik].T, self.data.sk[ik]) - if error > 1e-6: - print("Warning: eigenstates are not orthonormal, re-diagonalize Hamiltonian and overlap matrix") - eigvals, eigvecs = eigh(self.data.hk[ik], self.data.sk[ik]) - self.data.psi_lcao[ik] = eigvecs.T - self.data.energies[ik] = eigvals - print("-"*50, "\nOrthonormalities check finished. \nThen band selection according to user settings.") - # kpoints bands orbitals - self.data.psi_lcao = [np.array(psi_k)[band_range[0]:band_range[1], :] for psi_k in self.data.psi_lcao] - self.data.energies = [np.array(energies_k)[band_range[0]:band_range[1]] for energies_k in self.data.energies] - - # transpose psi_lcao to make it consistent with psi_chi - self.data.psi_lcao = [psi_k.T for psi_k in self.data.psi_lcao] - - def resize(self) -> None: - """resize the data container according to the data read from files or after basis filtering - - Raises: - ValueError: if AO filtered set cannot span larger space than the eigenvectors of the Hamiltonian - """ - self.data.nbands = self.data.psi_lcao[0].shape[1] - self.data.nchi = self.data.saok[0].shape[0] - self.data.nphi = self.data.hk[0].shape[0] - - _m = self.data.nchi - self.data.psi_lcao[0].shape[1] - if _m < 0: - print("number of AOs:", self.data.nchi) - print("number of bands selected in solved occupied eigenstates:", self.data.psi_lcao[0].shape[1]) - raise ValueError("current filtered AO set cannot span larger space than the eigenvectors of the Hamiltonian") - - self.data.psi_chi = [np.zeros( - (self.data.nphi, self.data.nchi)) for ik in range(self.data.nkpts)] - self.data.psi_qo = [np.zeros( - (self.data.nphi, self.data.nchi)) for ik in range(self.data.nkpts)] # how many AOs are there, how many QOs there are. - self.data.psi_chi_para = [np.zeros( - (self.data.nphi, self.data.nchi)) for ik in range(self.data.nkpts)] - self.data.psi_chi_orth = [np.zeros( - (self.data.nphi, self.data.nchi)) for ik in range(self.data.nkpts)] - self.data.wk = [np.zeros( - (self.data.nchi, self.data.nchi)) for ik in range(self.data.nkpts)] - self.data.psi_complem = [np.zeros( - (self.data.nphi, _m)) for ik in range(self.data.nkpts)] - self.data.psi_exten = [np.zeros( - (self.data.nphi, self.data.nchi)) for ik in range(self.data.nkpts)] - self.data.hqok = [np.zeros( - (self.data.nchi, self.data.nchi)) for ik in range(self.data.nkpts)] - self.data.sqok = [np.zeros( - (self.data.nchi, self.data.nchi)) for ik in range(self.data.nkpts)] - \ No newline at end of file diff --git a/tools/qo/abacus2qo/components/driver.py b/tools/qo/abacus2qo/components/driver.py deleted file mode 100644 index 7321115df7..0000000000 --- a/tools/qo/abacus2qo/components/driver.py +++ /dev/null @@ -1,195 +0,0 @@ -import abacus2qo.components.data_manager as dm -import abacus2qo.components.basis_filter as bf -import abacus2qo.components.calculator as cal -""" -1. filter out all irrelevant AOs from overlap matrix of AO in NAO representation -? maybe AO can also be normalized here? -2. should make sure still have number of AO larger than number of bands -3. calculate wk, diagonalize wk, get eigenvalues and eigenvectors, ... -4. get QO, along with the full Hamiltonian in its representation -5. multiply QO with psi_lcao, get psi (selected band range) in QO representation, calculate density matrix, get minimal set of QO -6. reproduce the band structure to check the accuracy -""" - -class toQO_Driver: - """we have too many toQO_* classes! Let it drive them - """ - dm_: dm.toQO_DataManager - bf_: bf.toQO_BasisFilter - cal_: cal.toQO_Calculator - - ao_basis_indices_: list - - def __init__(self) -> None: - pass - - def initialize(self, - path, - nkpts, - calculation, - band_range, - overlap_matrix_filter_threshold = 0.0): - - self.dm_ = dm.toQO_DataManager() - self.dm_.read(nkpts, calculation, path, band_range) - - - self.bf_ = bf.toQO_BasisFilter() - self.bf_.set_overlap_filter(overlap_matrix_filter_threshold) - - self.dm_.resize() - - self.cal_ = cal.toQO_Calculator() - - def filter_redundant_ao(self): - """filter out all irrelevant AOs from overlap matrix of AO in NAO representation - - Returns: - selected_basis_indices (list): list of list, each list contains indices of selected basis for each kpoint - - Notes: - This function is not practically used presently, because the charge for hydrogen-like orbital cannot be set very physically. - """ - print(""" ----------------------------------------------------------------------------------- -$$$$$$$$\ $$$$$$\ $$\ $$$$$$$$\ $$$$$$$$\ $$$$$$$\ $$$$$$\ $$$$$$\ -$$ _____|\_$$ _|$$ | \__$$ __|$$ _____|$$ __$$\ $$ __$$\ $$ __$$\ -$$ | $$ | $$ | $$ | $$ | $$ | $$ | $$ / $$ |$$ / $$ | -$$$$$\ $$ | $$ | $$ | $$$$$\ $$$$$$$ | $$$$$$$$ |$$ | $$ | -$$ __| $$ | $$ | $$ | $$ __| $$ __$$< $$ __$$ |$$ | $$ | -$$ | $$ | $$ | $$ | $$ | $$ | $$ | $$ | $$ |$$ | $$ | -$$ | $$$$$$\ $$$$$$$$\ $$ | $$$$$$$$\ $$ | $$ | $$ | $$ | $$$$$$ | -\__| \______|\________|\__| \________|\__| \__| \__| \__| \______/ ----------------------------------------------------------------------------------- -In this step, irrelevant AOs will be filtered out from overlap matrix of AO in NAO -representation. - """) - filtered_saok = [] - selected_basis_indices = [] - for ik in range(self.dm_.data.nkpts): - sk = self.dm_.data.sk[ik] - saok = self.dm_.data.saok[ik] - _selected_basis_indices, _saok = self.bf_.filter_via_overlap_matrix(sk, saok) - selected_basis_indices.append(_selected_basis_indices) - filtered_saok.append(_saok) - # update data - - self.dm_.data.saok = filtered_saok - self.dm_.resize() - return selected_basis_indices - - def space_expansion(self): - """Expand space to include orthogonal components to occupied states constructed from the arbitrary AO set. - """ - print(""" ----------------------------------------------------------------------------------------- -$$$$$$$$\ $$\ $$\ $$$$$$$\ $$$$$$\ $$\ $$\ $$$$$$\ $$$$$$\ $$$$$$\ $$\ $$\ -$$ _____|$$ | $$ |$$ __$$\ $$ __$$\ $$$\ $$ |$$ __$$\ \_$$ _|$$ __$$\ $$$\ $$ | -$$ | \$$\ $$ |$$ | $$ |$$ / $$ |$$$$\ $$ |$$ / \__| $$ | $$ / $$ |$$$$\ $$ | -$$$$$\ \$$$$ / $$$$$$$ |$$$$$$$$ |$$ $$\$$ |\$$$$$$\ $$ | $$ | $$ |$$ $$\$$ | -$$ __| $$ $$< $$ ____/ $$ __$$ |$$ \$$$$ | \____$$\ $$ | $$ | $$ |$$ \$$$$ | -$$ | $$ /\$$\ $$ | $$ | $$ |$$ |\$$$ |$$\ $$ | $$ | $$ | $$ |$$ |\$$$ | -$$$$$$$$\ $$ / $$ |$$ | $$ | $$ |$$ | \$$ |\$$$$$$ |$$$$$$\ $$$$$$ |$$ | \$$ | -\________|\__| \__|\__| \__| \__|\__| \__| \______/ \______| \______/ \__| \__| ----------------------------------------------------------------------------------------- -In this step, space spanned by selected bands-corresponding states will merge with -space spanned by orthogonal components of the arbitrary AO set.""") - for ik in range(self.dm_.data.nkpts): - print("-"*50, "\nFor k-point No.", ik) - self.dm_.data.psi_chi[ik] = self.cal_.projto_nao(self.dm_.data.sk[ik], self.dm_.data.saok[ik]) - self.dm_.data.psi_chi_para[ik] = self.cal_.projto_eigstate(self.dm_.data.psi_lcao[ik], self.dm_.data.saok[ik]) - self.dm_.data.psi_chi_orth[ik] = self.dm_.data.psi_chi[ik] - self.dm_.data.psi_chi_para[ik] - print("Back check the orthogonality between occupied states and constructed virtual states.") - _zero = self.dm_.data.psi_chi_orth[ik].conj().T @ self.dm_.data.sk[ik] @ self.dm_.data.psi_lcao[ik] - print("The result should be zero, and the result is: \n", _zero) - m = self.dm_.data.nchi - self.dm_.data.nbands - print("Number of empty states is: ", m, "\nNumber of bands to reproduce is: ", self.dm_.data.nbands, "\nNumber of basis functions is: ", self.dm_.data.nchi) - if m < 0: - raise ValueError("Number of AOs is smaller than number of bands selected.") - self.dm_.data.psi_complem[ik] = self.cal_.canonical_orthogonalization(self.dm_.data.psi_chi_orth[ik], self.dm_.data.sk[ik], m) - self.dm_.data.psi_exten[ik] = self.cal_.merge_space(self.dm_.data.psi_lcao[ik], - self.dm_.data.psi_complem[ik], - self.dm_.data.hk[ik], - self.dm_.data.sk[ik]) - print("-"*50) - - def reproduce_hamiltonian(self, Rs: list, test: str = "no"): - """get QO, reproduce selected pieces of energy spectrum - - Args: - Rs (list): list of R vectors - test (str, optional): test mode. Defaults to "no". - - supported test modes: - "no": no test - "lcao": return H(R) and S(R) in LCAO rep. - "w/o diag": return H(R) and S(R) in QO rep. without diagonalization to check the accuracy of the algorithm - - """ - print(""" ---------------------------------------------------------------------------------------------- -$$$$$$$$\ $$$$$$$\ $$$$$$\ $$\ $$\ $$$$$$\ $$$$$$$$\ $$$$$$\ $$$$$$$\ $$\ $$\ -\__$$ __|$$ __$$\ $$ __$$\ $$$\ $$ |$$ __$$\ $$ _____|$$ __$$\ $$ __$$\ $$$\ $$$ | - $$ | $$ | $$ |$$ / $$ |$$$$\ $$ |$$ / \__|$$ | $$ / $$ |$$ | $$ |$$$$\ $$$$ | - $$ | $$$$$$$ |$$$$$$$$ |$$ $$\$$ |\$$$$$$\ $$$$$\ $$ | $$ |$$$$$$$ |$$\$$\$$ $$ | - $$ | $$ __$$< $$ __$$ |$$ \$$$$ | \____$$\ $$ __| $$ | $$ |$$ __$$< $$ \$$$ $$ | - $$ | $$ | $$ |$$ | $$ |$$ |\$$$ |$$\ $$ |$$ | $$ | $$ |$$ | $$ |$$ |\$ /$$ | - $$ | $$ | $$ |$$ | $$ |$$ | \$$ |\$$$$$$ |$$ | $$$$$$ |$$ | $$ |$$ | \_/ $$ | - \__| \__| \__|\__| \__|\__| \__| \______/ \__| \______/ \__| \__|\__| \__| ---------------------------------------------------------------------------------------------- -In this step, H(k) (in NAO rep.) will be transformed into the one in Quasi-atomic Orbital (QO) -rep. - """) - if test == "no": - for ik in range(self.dm_.data.nkpts): - print("-"*50, "\nFor k-point No.", ik) - self.dm_.data.psi_qo[ik] = self.cal_.calculate_qo( - self.dm_.data.saok[ik], - self.dm_.data.psi_exten[ik], - self.dm_.data.sk[ik]) - self.dm_.data.hqok[ik], self.dm_.data.sqok[ik] = self.cal_.calculate_hqok( - self.dm_.data.psi_qo[ik], - self.dm_.data.hk[ik], - self.dm_.data.sk[ik]) - print("-"*50) - print(""" ----------------------------------------------------------------------------------------- -$$\ $$\ $$\ $$\ $$$$$$$$\ $$$$$$\ $$\ $$$$$$$\ $$$$$$\ $$\ $$\ $$$$$$\ -$$ | $$ |$$$\ $$ |$$ _____|$$ __$$\ $$ | $$ __$$\ \_$$ _|$$$\ $$ |$$ __$$\ -$$ | $$ |$$$$\ $$ |$$ | $$ / $$ |$$ | $$ | $$ | $$ | $$$$\ $$ |$$ / \__| -$$ | $$ |$$ $$\$$ |$$$$$\ $$ | $$ |$$ | $$ | $$ | $$ | $$ $$\$$ |$$ |$$$$\ -$$ | $$ |$$ \$$$$ |$$ __| $$ | $$ |$$ | $$ | $$ | $$ | $$ \$$$$ |$$ |\_$$ | -$$ | $$ |$$ |\$$$ |$$ | $$ | $$ |$$ | $$ | $$ | $$ | $$ |\$$$ |$$ | $$ | -\$$$$$$ |$$ | \$$ |$$ | $$$$$$ |$$$$$$$$\ $$$$$$$ |$$$$$$\ $$ | \$$ |\$$$$$$ | - \______/ \__| \__|\__| \______/ \________|\_______/ \______|\__| \__| \______/ ----------------------------------------------------------------------------------------- -In this step, H(k) will be converted into H(R). However, one should be sure the number of -k points defined in ABACUS KPT is larger than number of supercells searched in R space to -avoid information loss. ----------------------------------------------------------------------------------------- -Note1: present QO algorithm is numerically unstable, reproduce of band structure is not - guaranteed. It is because the linear dependence of QO, yields eigenvalue of their - overlap matrix to be zero, which brings numerical instability. ----------------------------------------------------------------------------------------- -Note2: the sequence of output matrix in QO representation is, (it, ia, l, zeta, m), "it" - is the index of atom type, "ia" is the index of atom of present type, "l" is the - angular momentum, "zeta" is the multiplicities of present angular momentum, "m" is - magnetic quantum number, arranged in the order in accordance with ABACUS: - Y00, Y10, Y11, Y1-1, Y20, Y21, Y2-1, Y22, Y2-2, ... ----------------------------------------------------------------------------------------- - """) - result_HR = [] - result_SR = [] - for R in Rs: - - if test == "no": - hqoR = self.cal_.unfolding_Hk(self.dm_.data.hqok, self.dm_.data.kpoints, R) - sqoR = self.cal_.unfolding_Hk(self.dm_.data.sqok, self.dm_.data.kpoints, R) - result_HR.append(hqoR) - result_SR.append(sqoR) - if test == "lcao": - hR = self.cal_.unfolding_Hk(self.dm_.data.hk, self.dm_.data.kpoints, R) - sR = self.cal_.unfolding_Hk(self.dm_.data.sk, self.dm_.data.kpoints, R) - result_HR.append(hR) - result_SR.append(sR) - return result_HR, result_SR \ No newline at end of file diff --git a/tools/qo/abacus2qo/components/safe_guard.py b/tools/qo/abacus2qo/components/safe_guard.py deleted file mode 100644 index 07973afc10..0000000000 --- a/tools/qo/abacus2qo/components/safe_guard.py +++ /dev/null @@ -1,92 +0,0 @@ -import numpy as np - -class toQO_SafeGuard: - """I am not that soap. Use to check hermiticity, orthonormality, rank, linear dependence of matrices - """ - def __init__(self) -> None: - pass - - def initialize(self) -> None: - """initialize the class - """ - pass - - def check_hermiticity(self, mat: np.ndarray) -> float: - """check if a matrix is Hermite - - Args: - mat (np.ndarray): matrix to be checked - """ - error = np.linalg.norm(mat - mat.conj().T) - print("Error of matrix being Hermite is: %14.10e"%error) - return error - - def check_orthonormality(self, mat: np.ndarray) -> float: - """check if a matrix is orthonormal - - Args: - mat (np.ndarray): matrix to be checked - """ - error = np.linalg.norm(mat - np.eye(mat.shape[-1])) - print("Error of matrix being orthonormal is: %14.10e"%error) - return error - - def check_identity(self, mat: np.ndarray) -> float: - """check if a matrix is identity - - Args: - mat (np.ndarray): matrix to be checked - """ - error = np.linalg.norm(mat - np.eye(mat.shape[-1])) - print("Error of matrix being identity is: %14.10e"%error) - return error - - def check_diagonal(self, mat: np.ndarray) -> float: - """check if a matrix is diagonal - - Args: - mat (np.ndarray): matrix to be checked - """ - error = np.linalg.norm(mat - np.diag(np.diag(mat))) - print("Error of matrix being diagonal is: %14.10e"%error) - return error - - def check_rank(self, mat: np.ndarray, tol = -1) -> int: - """check rank of a matrix - - Args: - mat (np.ndarray): matrix to be checked - """ - if tol < 0: - rank = np.linalg.matrix_rank(mat) - else: - rank = np.linalg.matrix_rank(mat, tol = tol) - print("Rank of matrix is: ", rank) - return rank - - def check_linear_dependence(self, mat: np.ndarray) -> int: - """check if a matrix is linearly dependent - - Args: - mat (np.ndarray): matrix to be checked - """ - rank = np.linalg.matrix_rank(mat) - if rank < mat.shape[-1]: - print("Matrix is linearly dependent!") - exit(1) - else: - return rank - - def check_eigenstates(self, eigenvec: np.ndarray, ovlp: np.ndarray) -> float: - """check if eigenstates are orthonormal - - Args: - eigenvec (np.ndarray): eigenstates - ovlp (np.ndarray): overlap matrix - """ - mat = eigenvec.conj().T @ ovlp @ eigenvec - print("Analyzing matrix of shape: ", mat.shape) - error = np.linalg.norm(mat - np.eye(eigenvec.shape[-1])) - print("Error of eigenstates being orthonormal is: %14.10e"%error) - - return error \ No newline at end of file diff --git a/tools/qo/abacus2qo/tools/__init__.py b/tools/qo/abacus2qo/tools/__init__.py deleted file mode 100644 index e69de29bb2..0000000000 diff --git a/tools/qo/abacus2qo/tools/basic_functions.py b/tools/qo/abacus2qo/tools/basic_functions.py deleted file mode 100644 index db0f94b0b9..0000000000 --- a/tools/qo/abacus2qo/tools/basic_functions.py +++ /dev/null @@ -1,43 +0,0 @@ -import numpy as np -def make_complex(number: str) -> complex: - """make complex number from string - - Args: - number (str): string of number - - Returns: - complex: complex number - """ - - read_real = False - read_imag = False - real = "" - imag = "" - for letter in number: - if letter == "(": - read_real = True - continue - elif letter == ",": - read_real = False - read_imag = True - continue - elif letter == ")": - read_imag = False - continue - else: - if read_real: - real += letter - elif read_imag: - imag += letter - try: - real = float(real) - imag = float(imag) - result = real+imag*1j - return result - except ValueError: - print(f"Error: {number}") - raise ValueError - -if __name__ == "__main__": - - print(make_complex("(4.32385486762811e-02,-1.59403849290982e-02)")) \ No newline at end of file diff --git a/tools/qo/abacus2qo/tools/hamiltonian.py b/tools/qo/abacus2qo/tools/hamiltonian.py deleted file mode 100644 index acc33a8366..0000000000 --- a/tools/qo/abacus2qo/tools/hamiltonian.py +++ /dev/null @@ -1,84 +0,0 @@ -import numpy as np -from abacus2qo.tools.basic_functions import make_complex - -def recover_from_upper_triangle(matrix: list) -> np.ndarray: - """recover full matrix from upper-right triangle - - Args: - matrix (list): matrix of upper-right triangle - - Returns: - np.ndarray: full matrix - """ - matrix = np.array(matrix) - matrix = matrix + matrix.conj().T - for i in range(len(matrix)): - matrix[i][i] = matrix[i][i] / 2 - return matrix - -def parse(nkpts: int, path = "./") -> tuple[list, list, int]: - """parse hamiltonian and overlap matrix H(k) and S(k) - - Args: - nkpts (int): number of kpoints - """ - hamiltonian = [] - overlap = [] - num_bands = 0 - if path[-1] != "/": - path += "/" - for ik in range(nkpts): - h_fname = path + f"data-{ik}-H" - s_fname = path + f"data-{ik}-S" - - hamiltonian_k = [] - overlap_k = [] - # first read Hamiltonian H(k) - with open(h_fname, "r") as f: - lines = f.readlines() - for il, line in enumerate(lines): - if len(line.strip()) == 0: - continue - hamiltonian_k_b = [] - words = line.strip().split() - if il == 0: - num_bands = int(words[0]) - # delete the first element - words.pop(0) - for icol in range(il): - hamiltonian_k_b.append(complex(0, 0)) - for word in words: - hamiltonian_k_b.append(make_complex(word)) - - hamiltonian_k.append(hamiltonian_k_b) - - hamiltonian_k = recover_from_upper_triangle(hamiltonian_k) - hamiltonian.append(hamiltonian_k) - - # then read overlap matrix S(k) - with open(s_fname, "r") as f: - lines = f.readlines() - for il, line in enumerate(lines): - if len(line.strip()) == 0: - continue - overlap_k_b = [] - words = line.strip().split() - if il == 0: - num_bands = int(words[0]) - # delete the first element - words.pop(0) - for icol in range(il): - overlap_k_b.append(complex(0, 0)) - for word in words: - overlap_k_b.append(make_complex(word)) - - overlap_k.append(overlap_k_b) - - overlap_k = recover_from_upper_triangle(overlap_k) - overlap.append(overlap_k) - - #return recover_from_upper_triangle(np.array(hamiltonian)), recover_from_upper_triangle(np.array(overlap)), num_bands - return hamiltonian, overlap, num_bands - -if __name__ == "__main__": - pass \ No newline at end of file diff --git a/tools/qo/abacus2qo/tools/kpoints.py b/tools/qo/abacus2qo/tools/kpoints.py deleted file mode 100644 index 5824998f98..0000000000 --- a/tools/qo/abacus2qo/tools/kpoints.py +++ /dev/null @@ -1,89 +0,0 @@ -import numpy as np - -def check_if_reduce(path = "./"): - - if path[-1] != "/": - path += "/" - log_fname = path + f"kpoints" - line = "start" - - with open(log_fname, "r") as f: - contents = f.readlines() - for line in contents: - line = line.strip() - if line.startswith("KPT DIRECT_X"): - return True - return False - -def parse(path = "./"): - if check_if_reduce(path): - return parse_symmetry_true(path) - else: - return parse_symmetry_false(path) - -def parse_symmetry_true(path = "./"): - - if path[-1] != "/": - path += "/" - log_fname = path + f"kpoints" - line = "start" - - kpoints = [] - equiv_kpoints = [] - - read_kpoint_reduction_information = False - with open(log_fname, "r") as f: - while len(line) > 0: - line = f.readline().strip() - - if line.startswith("KPT DIRECT_X"): - read_kpoint_reduction_information = True - continue - if line.startswith("nkstot = "): - nkpts = line.split()[2] - - if read_kpoint_reduction_information: - - kpt_ind, x, y, z, irr_kpt_ind, _x, _y, _z = line.split() - if int(irr_kpt_ind) > len(kpoints): - kpoint = np.array([float(_x), float(_y), float(_z)]) - kpoints.append(kpoint) - equiv_kpoints.append([kpoint]) - else: - kpoint = np.array([float(x), float(y), float(z)]) - equiv_kpoints[int(irr_kpt_ind)-1].append(kpoint) - - if int(kpt_ind) == int(nkpts): - break - continue - - - return kpoints, equiv_kpoints - -def parse_symmetry_false(path = "./"): - - if path[-1] != "/": - path += "/" - log_fname = path + f"kpoints" - line = "start" - - kpoints = [] - nkpts = 0 - with open(log_fname, "r") as f: - contents = f.readlines() - for line in contents: - line = line.strip() - if line.startswith("nkstot now = "): - nkpts = int(line.split()[-1]) - if line[0].isdigit(): - kpoint = np.array([float(x) for x in line.split()[1:4]]) - kpoints.append(kpoint) - nkpts -= 1 - if nkpts == 0: - break - return kpoints, kpoints - -if __name__ == "__main__": - kpoints, equiv_kpoints = parse("./work") - print(kpoints) - print(equiv_kpoints) \ No newline at end of file diff --git a/tools/qo/abacus2qo/tools/qo_ovlp.py b/tools/qo/abacus2qo/tools/qo_ovlp.py deleted file mode 100644 index 547897c864..0000000000 --- a/tools/qo/abacus2qo/tools/qo_ovlp.py +++ /dev/null @@ -1,40 +0,0 @@ -import numpy as np -from abacus2qo.tools.basic_functions import make_complex - -def parse(nkpts: int, path = "./") -> tuple[list, list]: - """read QO overlap matrix S(k) from file - - Args: - nkpts (int): number of kpoints - """ - qo_ovlp = [] - kpoints = [] - - if path[-1] != "/": - path += "/" - for ik in range(nkpts): - qo_fname = path + f"QO_ovlp_{ik}.dat" - qo_ovlp_k = [] - with open(qo_fname, "r") as f: - lines = f.readlines() - for line in lines: - if line.startswith("KPOINT_COORDINATE:"): - _words = line.split(":") - _words = _words[-1].split() - kpoints.append(np.array([float(number) for number in _words])) - else: - qo_ovlp_k.append([make_complex(number) for number in line.split()]) - qo_ovlp.append(np.array(qo_ovlp_k)) - - _qo_ovlp_R0 = np.zeros_like(qo_ovlp[0]) - for ik in range(nkpts): - _qo_ovlp_R0 += qo_ovlp[ik] - # sum over all element imaginary parts and check if near zero - if np.sum(np.abs(np.imag(_qo_ovlp_R0))) > 1e-6: - raise ValueError("kpoint symmetry error, fail to cancel imaginary part at R = 0.") - return qo_ovlp, kpoints - -if __name__ == "__main__": - qo_ovlp = parse(2) - print(qo_ovlp) - print(qo_ovlp.shape) \ No newline at end of file diff --git a/tools/qo/abacus2qo/tools/wavefunction.py b/tools/qo/abacus2qo/tools/wavefunction.py deleted file mode 100644 index 2adfa86277..0000000000 --- a/tools/qo/abacus2qo/tools/wavefunction.py +++ /dev/null @@ -1,65 +0,0 @@ -import numpy as np - -def parse(nkpts: int, path: str = "./") -> tuple[list, list, list]: - """parse LOWF_K_*.txt file generated by ABACUS with input parameter out_wfc_lcao 1. - - Args: - nkpts (int): number of k-points - """ - wfc = [] - kpoints = [] - energies = [] - if path[-1] != "/": - path += "/" - for ik in range(1, nkpts+1): - fname = path + f"LOWF_K_{ik}.txt" - wfc_k = [] - energies_k = [] - wfc_k_b = [] - with open(fname, "r") as f: - lines = f.readlines() - for line in lines: - line = line.strip() - if line.endswith("(index of k points)"): - #print("kpoint index:", line) - pass - elif line.endswith("(number of bands)"): - #print("number of bands:", line) - pass - elif line.endswith("(number of orbitals)"): - #print("number of orbitals:", line) - pass - elif line.endswith("(band)"): - #print("band index:", line) - pass - elif line.endswith("(Ry)"): - energies_k.append(float(line.split()[0])) - pass - elif line.endswith("(Occupations)"): - occ = float(line.split()[0]) - if len(wfc_k_b) > 0: - # group every two elements into a complex number - wfc_k_b = np.array(wfc_k_b) - wfc_k_b = wfc_k_b[::2] + 1j * wfc_k_b[1::2] - wfc_k.append(wfc_k_b) - wfc_k_b = [] - else: - num_valid_words = len(line.split()) - if num_valid_words == 3: - #print("kpoint coordinates:", line) - kpoints.append(np.array([float(word) for word in line.split()])) - else: - for coeff in line.strip().split(): # remove leading and trailing whitespaces - wfc_k_b.append(float(coeff)) - # last band - wfc_k_b = np.array(wfc_k_b) - wfc_k_b = wfc_k_b[::2] + 1j * wfc_k_b[1::2] - wfc_k.append(wfc_k_b) - - wfc.append(np.array(wfc_k)) - energies.append(np.array(energies_k)) - - return wfc, kpoints, energies - -if __name__ == "__main__": - print(parse(8, path = "./examples/input/")) \ No newline at end of file diff --git a/tools/qo/examples/INPUT b/tools/qo/examples/INPUT index f86fca2cb0..3a93f615d6 100644 --- a/tools/qo/examples/INPUT +++ b/tools/qo/examples/INPUT @@ -1,13 +1,15 @@ INPUT_PARAMETERS -pseudo_dir ../../../tests/PP_ORB -orbital_dir ../../../tests/PP_ORB +pseudo_dir /abacus-develop/tests/PP_ORB/ +orbital_dir /abacus-develop/tests/PP_ORB/ -ecutwfc 50 +ecutwfc 100 scf_nmax 100 -scf_thr 1e-6 +scf_thr 1e-7 basis_type lcao qo_switch 1 # turn on QO analysis output qo_basis pswfc # use pseudowavefunction read from pseudopotentials qo_screening_coeff 0.5 # controls the shrink of AO basis, larger value means tighter basis -qo_thr 1e-6 # controls the accuracy of QO basis +qo_thr 1e-10 # controls the accuracy of QO basis +qo_strategy all # use all pswfc to form QO basis +symmetry -1 # to ensure matrix elements of symmetry kpoints can be cancelled diff --git a/tools/qo/examples/KPT b/tools/qo/examples/KPT index de4006c594..7c59c733af 100644 --- a/tools/qo/examples/KPT +++ b/tools/qo/examples/KPT @@ -1,4 +1,4 @@ K_POINTS 0 MP -5 5 5 0 0 0 +9 9 9 0 0 0 diff --git a/tools/qo/examples/STRU b/tools/qo/examples/STRU index 37805fca79..069387c455 100644 --- a/tools/qo/examples/STRU +++ b/tools/qo/examples/STRU @@ -1,23 +1,23 @@ #This test is from example/scf/lcao_Si but change to Pseudopotential with pswfc ATOMIC_SPECIES -Sn 1.000 Sn_dojo_nsoc.upf #MUST USE PSEUDOPOTENTIAL WITH PSWFC! +Si 1.000 Si_dojo_nsoc.upf #MUST USE PSEUDOPOTENTIAL WITH PSWFC! NUMERICAL_ORBITAL -Sn_dojo_6au.orb +./Si_gga_8au_100Ry_2s2p1d.orb LATTICE_CONSTANT -1.889726877 +1.8897261258369282 #Lattice constant LATTICE_VECTORS - 4.64567349 0.00000000 0.00000000 - 2.32283620 4.02327092 0.00000000 - 2.32283597 1.34109040 3.79317605 +0.0000000000 2.7218510000 2.7218510000 +2.7218510000 0.0000000000 2.7218510000 +2.7218510000 2.7218510000 0.0000000000 ATOMIC_POSITIONS -Direct -Sn -0.00 -2 - 0.87500000 0.87500000 0.87500000 m 1 1 1 - 0.12500000 0.12500000 0.12500000 m 1 1 1 +Direct #Cartesian(Unit is LATTICE_CONSTANT) +Si #Name of element +0.0 #Magnetic for this element. +2 #Number of atoms +0.00 0.00 0.00 1 1 1 #x,y,z, move_x, move_y, move_z +0.25 0.25 0.25 1 1 1 diff --git a/tools/qo/examples/reference/band/INPUT b/tools/qo/examples/reference/band/INPUT new file mode 100644 index 0000000000..8471a978c6 --- /dev/null +++ b/tools/qo/examples/reference/band/INPUT @@ -0,0 +1,11 @@ +INPUT_PARAMETERS +# Created by Atomic Simulation Enviroment +pseudo_dir /abacus-develop/tests/PP_ORB/ +orbital_dir /abacus-develop/tests/PP_ORB/ +ntype 1 +ecutwfc 100 +basis_type lcao +symmetry 0 +calculation nscf +out_band 1 +init_chg file diff --git a/tools/qo/examples/reference/band/KPT b/tools/qo/examples/reference/band/KPT new file mode 100644 index 0000000000..f763f1cde8 --- /dev/null +++ b/tools/qo/examples/reference/band/KPT @@ -0,0 +1,11 @@ +K_POINTS +8 +Line +0.0000000000 0.0000000000 0.0000000000 30 +0.5000000000 0.0000000000 0.5000000000 30 +0.6250000000 0.2500000000 0.6250000000 1 +0.3750000000 0.3750000000 0.7500000000 30 +0.0000000000 0.0000000000 0.0000000000 30 +0.5000000000 0.5000000000 0.5000000000 30 +0.5000000000 0.2500000000 0.7500000000 30 +0.5000000000 0.0000000000 0.5000000000 1 diff --git a/tools/qo/examples/reference/band/STRU b/tools/qo/examples/reference/band/STRU new file mode 100644 index 0000000000..069387c455 --- /dev/null +++ b/tools/qo/examples/reference/band/STRU @@ -0,0 +1,23 @@ +#This test is from example/scf/lcao_Si but change to Pseudopotential with pswfc + +ATOMIC_SPECIES +Si 1.000 Si_dojo_nsoc.upf #MUST USE PSEUDOPOTENTIAL WITH PSWFC! + +NUMERICAL_ORBITAL +./Si_gga_8au_100Ry_2s2p1d.orb + +LATTICE_CONSTANT +1.8897261258369282 #Lattice constant + +LATTICE_VECTORS +0.0000000000 2.7218510000 2.7218510000 +2.7218510000 0.0000000000 2.7218510000 +2.7218510000 2.7218510000 0.0000000000 + +ATOMIC_POSITIONS +Direct #Cartesian(Unit is LATTICE_CONSTANT) +Si #Name of element +0.0 #Magnetic for this element. +2 #Number of atoms +0.00 0.00 0.00 1 1 1 #x,y,z, move_x, move_y, move_z +0.25 0.25 0.25 1 1 1 diff --git a/tools/qo/examples/reference/scf/INPUT b/tools/qo/examples/reference/scf/INPUT new file mode 100644 index 0000000000..c2badae774 --- /dev/null +++ b/tools/qo/examples/reference/scf/INPUT @@ -0,0 +1,9 @@ +INPUT_PARAMETERS +pseudo_dir /abacus-develop/tests/PP_ORB/ +orbital_dir /abacus-develop/tests/PP_ORB/ + +ecutwfc 100 +scf_nmax 100 +scf_thr 1e-7 +basis_type lcao +out_chg 1 diff --git a/tools/qo/examples/reference/scf/KPT b/tools/qo/examples/reference/scf/KPT new file mode 100644 index 0000000000..7c59c733af --- /dev/null +++ b/tools/qo/examples/reference/scf/KPT @@ -0,0 +1,4 @@ +K_POINTS +0 +MP +9 9 9 0 0 0 diff --git a/tools/qo/examples/reference/scf/STRU b/tools/qo/examples/reference/scf/STRU new file mode 100644 index 0000000000..069387c455 --- /dev/null +++ b/tools/qo/examples/reference/scf/STRU @@ -0,0 +1,23 @@ +#This test is from example/scf/lcao_Si but change to Pseudopotential with pswfc + +ATOMIC_SPECIES +Si 1.000 Si_dojo_nsoc.upf #MUST USE PSEUDOPOTENTIAL WITH PSWFC! + +NUMERICAL_ORBITAL +./Si_gga_8au_100Ry_2s2p1d.orb + +LATTICE_CONSTANT +1.8897261258369282 #Lattice constant + +LATTICE_VECTORS +0.0000000000 2.7218510000 2.7218510000 +2.7218510000 0.0000000000 2.7218510000 +2.7218510000 2.7218510000 0.0000000000 + +ATOMIC_POSITIONS +Direct #Cartesian(Unit is LATTICE_CONSTANT) +Si #Name of element +0.0 #Magnetic for this element. +2 #Number of atoms +0.00 0.00 0.00 1 1 1 #x,y,z, move_x, move_y, move_z +0.25 0.25 0.25 1 1 1 diff --git a/tools/qo/main.py b/tools/qo/main.py deleted file mode 100644 index 8162daa55d..0000000000 --- a/tools/qo/main.py +++ /dev/null @@ -1,38 +0,0 @@ -"""Main program of Quasiatomic Orbital (QO) analysis postprocessing package - -Author: ykhuang -Date: 2023/12/25, Merry Christmas - -Instruction ------------ -This program is used to postprocess the output of ABACUS calculation. To -obtain files needed, please set the following parameters in ABACUS input file: -``` -# Quasiatomic Orbital Analysis -qo_switch 1 # turn on quasiatomic orbital analysis -qo_basis pswfc # use pseudowavefunction as QO corresponding AO -qo_thr 1e-6 # controls the realspace spreading of AO -qo_screening_coeff 0.5 # controls the exponential decay-manner of AO -``` -""" -import abacus2qo.components.driver as driver -import numpy as np - -if __name__ == "__main__": - # example of QO - path = "./examples/OUT.ABACUS/" - # fully non-reduced kpoints - nkpts = 125 - # band range to reproduce - band_range = (0, 13) - - # create driver - d_ = driver.toQO_Driver() - # initialize driver - d_.initialize(path, nkpts, "scf", band_range) - # expand space spanned by selected bands - d_.space_expansion() - # calculate H(R) and S(R) in R-space in QO representation - HRs, SRs = d_.reproduce_hamiltonian(Rs=[ - np.array([0, 0, 0]) - ]) \ No newline at end of file diff --git a/tools/qo/postprocess.py b/tools/qo/postprocess.py new file mode 100644 index 0000000000..4db5b48476 --- /dev/null +++ b/tools/qo/postprocess.py @@ -0,0 +1,690 @@ +""" +---------------------------------------------------------------------------------- +| ABACUS LCAO2QO(Quasiatomic Orbital) Instruction | +---------------------------------------------------------------------------------- +| Author: Yike HUANG +| Email: huangyk@aisi.ac.cn +| Date: 2024-03-02 +| +| TO USE: +| Can either use this script in-place or by importing it as a module. +| 0. Run ABACUS: +| - Run ABACUS with the following settings: +| - `calculation scf` +| - `qo_switch 1` +| - `qo_basis hydrogen` # can be hydrogen, pswfc and nao +| - `qo_strategy energy-valence` # for more strategies, see ABACUS manual +| - `qo_thr 1e-10` # threshold for determing range of integral in realspace +| If want to turn on Slater Screening to mimic many-electron behavior (sort +| of), add the following settings: +| - `qo_screening_coeff 0.1` +| More details please refer to ABACUS manual. +| 1. In-place: +| - Put this script in the same directory as the output files of ABACUS. +| - Run the script with the command `python all_in_one.py` after modifying +| the parameters in the `if __name__ == "__main__":` block. +| 2. As a module: +| - Import this script as a module in another script. +| - Call the function `production` with the appropriate parameters, the +| function `production` returns the Hamiltonian and overlap matrices in +| QO representation along with valid kpoints direct coordinates. +| 9999. Unittest: +| - Run the script with the command `python all_in_one.py` to perform +| unittests after modifying the parameters in the `if __name__ == "__main__":` +| block. Unittest will read data stored in the same folder when performing +| production runs, so it is necessary targeting the correct path first. +---------------------------------------------------------------------------------- +""" + +import numpy as np +# TESTED +def cxx_topycomplex(num: str): + """cxx prints complex numbers in the form (a,b)""" + num = num.replace('(', '').replace(')', '') + num = num.split(',') + return complex(float(num[0]), float(num[1])) +# TESTED +def read_mat_hs(fhs): + """""" + with open(fhs, 'r') as f: + data = "".join(f.readlines()).replace('\n', '').split() + size = int(data[0]) + indexing = [] + for i in range(size): + indexing += [(i, j) for j in range(i, size)] + data = data[1:] + mat = np.zeros((size, size), dtype=np.complex128) + for inum, number in enumerate(data): + mat[indexing[inum]] = cxx_topycomplex(number) + mat = mat + mat.conj().T + for i in range(size): + mat[i, i] = mat[i, i] / 2 + return mat +# TESTED +def read_lowf(flowf): + """read the lowf matrix in k-space, and the k-vector. + + Args: + flowf (str): path to the file containing the lowf matrix + + Returns: + np.ndarray: lowf matrix in k-space + np.ndarray: k-vector + """ + with open(flowf, 'r') as f: + lines = [line.strip() for line in f.readlines()] + for line in lines: + if line.endswith("(number of bands)"): + nband = int(line.split()[0]) + elif line.endswith("(number of orbitals)"): + nlocal = int(line.split()[0]) + else: + continue + lowf_k = np.zeros((nlocal, nband), dtype=np.complex128) + kvec_c = None + ib, ilocal = 0, 0 + for line in lines: + if line.endswith("(band)"): + ib = int(line.split()[0]) - 1 + ilocal = 0 + continue + if not line.endswith(")"): + if line.count(" ") != 2: + # it means it is a band rather than coordinates of kpoint + nums = line.split() + c = [complex(float(nums[i]), float(nums[i+1])) for i in range(0, len(nums), 2)] + for i in range(len(c)): + lowf_k[ilocal, ib] = c[i] + ilocal += 1 + else: + kvec_c = np.array([float(x) for x in line.split()]) + return lowf_k, kvec_c +# TESTED +def read_ao_proj(fao_proj): + """Read the atomic orbital projection matrix in k-space, and the k-vector. + + Args: + fao_proj (str): path to the file containing the atomic orbital projection matrix + + Returns: + np.ndarray: atomic orbital projection matrix in k-space + np.ndarray: k-vector + """ + with open(fao_proj, 'r') as f: + lines = [line.strip() for line in f.readlines()] + kvec_d = None + if lines[0].startswith("KPOINT_COORDINATE"): + kvec_d = np.array([float(x) for x in lines[0].split()[-3:]]) + lines = lines[1:] + nlocal = len(lines[0].split()) + nao = len(lines) + ao_proj = np.zeros((nao, nlocal), dtype=np.complex128) + for i, line in enumerate(lines): + nums = [cxx_topycomplex(num) for num in line.split()] + for j in range(nlocal): + ao_proj[i, j] = nums[j] + # because in ao_proj, the row index is AO, we need to dagger it + return ao_proj.conj().T, kvec_d +# TESTED +import os +def read(path: str, + nk: int, + read_H: bool = True, + read_S: bool = True, + read_LOWF: bool = True, + read_QO: bool = True) -> tuple[list, list, list, list, list]: + """Read H(k), S(k), LOWF(k), and <\phi(k)|A(k)> from ABACUS output files. + + Args: + path (str): path to the directory containing the output files + nk (int): number of kpoints + read_H (bool, optional): whether to read H(k). Defaults to True. + read_S (bool, optional): whether to read S(k). Defaults to True. + read_LOWF (bool, optional): whether to read LOWF(k). Defaults to True. + read_QO (bool, optional): whether to read <\phi(k)|A(k)>. Defaults to True. + + Returns: + tuple[list, list, list, list, list]: H(k), S(k), LOWF(k), <\phi(k)|A(k)>, k-vectors + """ + hs_k, ss_k, lowfs_k, aos_proj_k, kvecs_c, kvecs_d = [], [], [], [], [], [] + if read_H: + for ik in range(nk): + fh = os.path.join(path, f"data-{ik}-H") + hs_k.append(read_mat_hs(fh)) + if read_S: + for ik in range(nk): + fs = os.path.join(path, f"data-{ik}-S") + ss_k.append(read_mat_hs(fs)) + if read_LOWF: + for ik in range(nk): + flowf = os.path.join(path, f"LOWF_K_{ik+1}.txt") + lowf_k, kvec_c = read_lowf(flowf) + lowfs_k.append(lowf_k) + kvecs_c.append(kvec_c) + if read_QO: + for ik in range(nk): + fao_proj = os.path.join(path, f"QO_ovlp_{ik}.dat") + ao_proj_k, kvec_d = read_ao_proj(fao_proj) + aos_proj_k.append(ao_proj_k) + kvecs_d.append(kvec_d) + + return hs_k, ss_k, lowfs_k, aos_proj_k, kvecs_d +# TESTED +def cal_denmat(lowf_k: np.ndarray) -> np.ndarray: + """calculate the density matrix from the lowf matrix. + + Args: + lowf_k (np.ndarray): lowf matrix in k-space + + Returns: + np.ndarray: density matrix in k-space + """ + return lowf_k@lowf_k.conj().T +# TESTED +import scipy as sp +def lowdin_onW(s_k: np.ndarray, denmat_k: np.ndarray, ao_proj_k: np.ndarray, m: int): + """Perform Lowdin orthogonalization on the W matrix, which is defined as W = A^+ * (S^-1 - D) * A. + + Args: + s_k (np.ndarray): overlap matrix in k-space + denmat_k (np.ndarray): density matrix in k-space + ao_proj_k (np.ndarray): atomic orbital projection matrix in k-space + m (int): number of bands wanted + + Returns: + np.ndarray: the lowdin orthogonalized matrix in the quasi-atomic orbital basis, can check + orthonormality with overlap matrix + """ + nrow = s_k.shape[0] + ncol = s_k.shape[1] + assert nrow == ncol + sinv_k = np.linalg.solve(s_k, np.eye(N=nrow, M=ncol, dtype=np.complex128)) + W_mat = ao_proj_k.conj().T@(sinv_k - denmat_k)@ao_proj_k + # forcely let diagonal elements to be real + W_mat = (W_mat + W_mat.conj().T) / 2 + # diagonalize W and get only m largest eigenvalues and corresponding eigenvectors + eigvals, eigvecs = sp.linalg.eigh(W_mat) + eigvals = eigvals[-m:] + eigvecs = eigvecs[:, -m:] + lambdas = np.diag(1.0/np.sqrt(eigvals)) + + return (sinv_k - denmat_k)@ao_proj_k@eigvecs@lambdas + # call return as lowf_k_bar +# TESTED +def cal_tb_hs(h_k: np.ndarray, s_k: np.ndarray, denmat_k_tilde: np.ndarray, ao_proj_k: np.ndarray, norm_thr: float = 1e-10): + """calculate the tight-binding Hamiltonian and overlap matrix in the quasi-atomic orbital basis. + + Args: + h_k (np.ndarray): Hamiltonian matrix in k-space + s_k (np.ndarray): overlap matrix in k-space + denmat_k_tilde (np.ndarray): density matrix in k-space + ao_proj_k (np.ndarray): atomic orbital projection matrix in k-space + norm_thr (float, optional): threshold for the norm of quasi-atomic orbitals. Defaults to 1e-10. + + Returns: + tuple: tight-binding Hamiltonian and overlap matrix in the quasi-atomic orbital basis. + """ + qo_k = denmat_k_tilde@ao_proj_k + # normalize qo_k + nqo = qo_k.shape[1] + for iqo in range(nqo): + norm_qo = np.sqrt(qo_k[:, iqo].conj().T@s_k@qo_k[:, iqo]) + if norm_qo < norm_thr: + raise ValueError(f"Norm of qo_k is too small: norm < {norm_thr}.") + qo_k[:, iqo] = qo_k[:, iqo] / norm_qo + + hqo_k = qo_k.conj().T@h_k@qo_k + sqo_k = qo_k.conj().T@s_k@qo_k + + return hqo_k, sqo_k +# TESTED +def kR_transform(dst, mats_in: list, srcs: list, direction: str = "R->k"): + """single-side transformation from k-space to R-space, or from R-space to k-space. + + Args: + dst (np.ndarray): destination vector + mats_in (list): input matrices + srcs (list): source vectors + direction (str, optional): "R->k" or "k->R". Defaults to "R->k". + + Returns: + np.ndarray: transformed matrix + """ + mat_out = np.zeros_like(mats_in[0], dtype=np.complex128) + if len(mats_in) != len(srcs): + print(f"Number of input matrices and source vectors are not equal: {len(mats_in)} (num of mats) != {len(srcs)} (num of srcs).") + print(f"Shape of first matrix: {mats_in[0].shape}.") + assert False + phase = 0.0+1.0j if direction == "R->k" else 0.0-1.0j + for mat_in, src in zip(mats_in, srcs): + arg = np.exp(phase * 2.0 * np.pi * np.dot(dst, src)) + mat_out += arg * mat_in + return mat_out / np.sqrt(len(mats_in)) +# TESTED +def k_extrapolation(mats_k: list, kvecs_in: list, Rs: list, kvecs_out: list): + """for kpoint extrapolation, from k-space to R-space, and then back to k-space. + + Args: + mats_k (list): input matrices in k-space + kvecs_in (list): input k-vectors + Rs (list): supercell vectors + kvecs_out (list): output k-vectors + + Returns: + list: output matrices in k-space + """ + mats_R = [kR_transform(Rs[i], mats_k, kvecs_in, direction="k->R") for i in range(len(Rs))] + mats_kpath = [kR_transform(kvecs_out[i], mats_R, Rs, direction="R->k") for i in range(len(kvecs_out))] + return mats_kpath + +def norm_filter_k(denmat_k: np.ndarray, s_k: np.ndarray, ao_proj_k: np.ndarray, norm_thr: float = 1e-10): + """filter out unbounded states represented by norms of atomic orbitals. + + Args: + denmat_k (np.ndarray): density matrix in k-space + s_k (np.ndarray): overlap matrix in k-space + ao_proj_k (np.ndarray): atomic orbital projection matrix in k-space + norm_thr (float, optional): threshold for the norm of atomic orbitals. Defaults to 1e-10. + + Returns: + list: indices of unbounded states + """ + ao_k = denmat_k@ao_proj_k + norms = [np.sqrt(ao_k[:, iao].conj().T@s_k@ao_k[:, iao]) for iao in range(ao_k.shape[1])] + unbounded_states = [iqo for iqo, norm in enumerate(norms) if norm < norm_thr] + return unbounded_states + +def nao_space_reduce(selected_indices: list, s_k: np.ndarray): + return s_k[:, selected_indices] + +"""Expect a file contain information describing row and column indices, like + + +0,0,0,0,0 0,0,0,1,0 0,0,0,2,0 0,0,1,0,0 + + +[it],[ia],[l],[izeta],[m] [it],[ia],[l],[izeta],[m] [it],[ia],[l],[izeta],[m] [it],[ia],[l],[izeta],[m] ... + + +""" +def read_matrix_index(findex: str): + with open(findex, 'r') as f: + lines = [line.strip() for line in f.readlines()] + # return bidirectional maps, tuple to index and index to tuple + # row index + row_index = lines[2].split() + row_index = [tuple([int(x) for x in index.split(",")]) for index in row_index] + row_index = {index: i for i, index in enumerate(row_index)} + row_index_inv = {i: index for i, index in enumerate(row_index)} + # col index + col_index = lines[4].split() + col_index = [tuple([int(x) for x in index.split(",")]) for index in col_index] + col_index = {index: i for i, index in enumerate(col_index)} + col_index_inv = {i: index for i, index in enumerate(col_index)} + return row_index, row_index_inv, col_index, col_index_inv + +import itertools as it +def mats_hs_k(path: str, nk: int, ib_min: int, ib_max: int = -1, qo: bool = True): + print("Number of kpoints: ", nk) + hs_k, ss_k, lowfs_k, aos_proj_k, kvecs_d = read(path, nk) + lowfs_k = lowfs_k if ib_max == -1 else [lowfs_k[ik][:, ib_min:ib_max] for ik in range(nk)] + + if not qo: + return hs_k, ss_k, kvecs_d + else: + # 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 + # present is TZDP, 3s3p2d, arrange like [1s, 2s, 3s, 1p-1, 1p0, 1p1, 2p-1, 2p0, 2p1, 3p-1, 3p0, 3p1, 1d-2, 1d-1, 1d0, 1d1, 1d2, 2d-2, 2d-1, 2d0, 2d1, 2d2] + # 0 1 2 3 4 5 6 7 8 9 10 11 12 + # for DZP, 2s2p1d, arrange like [1s, 2s, 1p-1, 1p0, 1p1, 2p-1, 2p0, 2p1, 1d-2, 1d-1, 1d0, 1d1, 1d2] + # only want 1s1p + # therefore the szv_indices from TZDP: + #szv_indices = [0, 3, 4, 5, 22, 25, 26, 27] + # therefore the szv_indices from DZP: + #szv_indices = [0, 2, 3, 4, 13, 15, 16, 17] + #dzp_indices = [0, 1, 3, 4, 5, 6, 7, 8, 12, 13, 14, 15, 16, 22, 23, 25, 26, 27, 28, 29, 30, 34, 35, 36, 37, 38] + #aos_proj_k = [ss_k[ik][:, szv_indices] for ik in range(nk)] + #print("Shape of ao_proj_k: ", aos_proj_k[0].shape) # (8, 44), 8 is the number of selected atomic orbitals, 44 is the number of numerical atomic orbitals + #aos_proj_k = ss_k + + nband = lowfs_k[0].shape[1] + print("Number of bands selected: ", nband) + m = aos_proj_k[0].shape[1] - nband + print("Calculating quasi-atomic orbitals, number of bands unoccupied: ", m) + denmats_k = [cal_denmat(lowf_k) for lowf_k in lowfs_k] + ubsts_k = [norm_filter_k(denmat_k=denmats_k[ik], s_k=ss_k[ik], ao_proj_k=aos_proj_k[ik], norm_thr=1e-8) for ik in range(nk)] + iks_ubsts = [ik for ik, ubsts in enumerate(ubsts_k) if len(ubsts) > 0] + print("Kpoints with unbounded states: ", iks_ubsts) + lowfs_bar_k = [lowdin_onW(s_k=ss_k[ik], + denmat_k=denmats_k[ik], + ao_proj_k=aos_proj_k[ik], + m=m) for ik in range(nk)] + lowfs_tilde_k = [np.concatenate((lowfs_k[ik], lowfs_bar_k[ik]), axis=1) for ik in range(nk)] + _one = [lowf_tilde_k.conj().T@ss_k[ik]@lowf_tilde_k for ik, lowf_tilde_k in enumerate(lowfs_tilde_k)] + # immediate check of orthonormality + for ik in range(nk): + for i, j in it.product(range(nband+m), repeat=2): + if i != j: + #print(_one[ik][i, j]) + assert abs(_one[ik][i, j].real) < 1e-6 and abs(_one[ik][i, j].imag) < 1e-6 + else: + #print(_one[ik][i, j]) + assert abs(_one[ik][i, j].real - 1) < 1e-5 and abs(_one[ik][i, j].imag) < 1e-6 + # go on + denmats_tilde_k = [cal_denmat(lowf_tilde_k) for lowf_tilde_k in lowfs_tilde_k] + hqos_k, sqos_k = [], [] + for ik in range(nk): + hqo_k, sqo_k = cal_tb_hs(hs_k[ik], ss_k[ik], denmats_tilde_k[ik], aos_proj_k[ik]) + hqos_k.append(hqo_k) + sqos_k.append(sqo_k) + return hqos_k, sqos_k, kvecs_d + +def monkhorst_pack(size): + """Construct a uniform sampling of k-space of given size.""" + if np.less_equal(size, 0).any(): + raise ValueError(f'Illegal size: {list(size)}') + kpts = np.indices(size).transpose((1, 2, 3, 0)).reshape((-1, 3)) + return (kpts - np.array(size)//2) + +def bad_condition_diag(Hk, Sk, eig_thr = 1e-10): + """diagonalize the generalized eigenvalue problem Hk * x = e * Sk * x, and return the eigenvalues. + + Args: + Hk (np.ndarray): Hamiltonian matrix + Sk (np.ndarray): overlap matrix + eig_thr (float, optional): for filtering out linear dependent components. Defaults to 1e-10. + + Returns: + eigvals: eigenvalues of the generalized eigenvalue problem, sorted in ascending order, padded with a large number if necessary. + """ + ndim = Hk.shape[0] + s_eigvals, s_eigvecs = np.linalg.eigh(Sk) + #print("eigenvalues of sk:\n", s_eigvals) + s_eigvecs = s_eigvecs[:, s_eigvals > eig_thr] + s_eigvals = s_eigvals[s_eigvals > eig_thr] + s_eigvals = np.diag(1.0/np.sqrt(s_eigvals)) + s_eigvecs = np.dot(s_eigvecs, s_eigvals) + Hk = np.dot(np.conj(s_eigvecs.T), np.dot(Hk, s_eigvecs)) + eigs = np.linalg.eigvalsh(Hk) + # fill in the missing eigenvalues with a large number + eigs = np.concatenate((eigs, np.ones(ndim-len(eigs))*1e6)) + return eigs + +def production(path: str, # path of OUT.${suffix} folder + nkpts: int, # number of kpoints in KPT file + ib_min: int = 0, # minimum index of bands selected to reproduce + ib_max: int = -1, # maximum index of bands selected to reproduce + diag_thr: float = 1e-10, # threshold for filtering out linear dependent components + error_estimation_with_lcao: bool = False): # whether to estimate errors with LCAO + """Production code for the whole process of calculating Hamiltonian and overlap matrices in QO basis. + + Args: + path (str): path to the directory containing the output files + nkpts (int): number of kpoints + ib_min (int, optional): minimum index of bands selected to reproduce. Defaults to 0. + ib_max (int, optional): maximum index of bands selected to reproduce. Defaults to -1. + diag_thr (float, optional): threshold for filtering out linear dependent components. Defaults to 1e-10. + error_estimation_with_lcao (bool, optional): whether to estimate errors with LCAO. Defaults to False. + + Returns: + list: Hamiltonian matrices in QO basis + list: overlap matrices in QO basis + list: k-vectors + """ + # first calculate things repsented by QO basis + hqos_k, sqos_k, kvecs_d = mats_hs_k(path, nkpts, ib_min, ib_max, qo=True) + eigs_qo = [bad_condition_diag(hqos_k[ik], sqos_k[ik], diag_thr) for ik in range(nkpts)] + eigs_qo = np.asarray(eigs_qo) + # then calculate things represented by AO basis, if necessary + if error_estimation_with_lcao: + haos_k, saos_k, kvecs_d = mats_hs_k(path, nkpts, ib_min, ib_max, qo=False) + eigs_ao = [bad_condition_diag(haos_k[ik], saos_k[ik], diag_thr) for ik in range(nkpts)] + eigs_ao = np.asarray(eigs_ao) + # calculate errors + e_errors_k = eigs_qo[:, ib_min:ib_max] - eigs_ao[:, ib_min:ib_max] + e_errors_k = [np.linalg.norm(e_errors_k[ik]) for ik in range(nkpts)] + exclude_kvecs = [ik for ik, e_error_k in enumerate(e_errors_k) if e_error_k > 1e-8] + print("Number of kpoints excluded due to large error:", len(exclude_kvecs), "\nTheir indices:", exclude_kvecs) + for ik in exclude_kvecs: + print("kpoint: ", kvecs_d[ik], "has an error of: ", e_errors_k[ik]) + print("eigvals of sqos_k: ", sp.linalg.eigvalsh(sqos_k[ik])) + else: + exclude_kvecs = [] + # exclude kpoints with large errors + hqos_k = [hqos_k[ik] for ik in range(nkpts) if ik not in exclude_kvecs] + sqos_k = [sqos_k[ik] for ik in range(nkpts) if ik not in exclude_kvecs] + kvecs_d = [kvecs_d[ik] for ik in range(nkpts) if ik not in exclude_kvecs] + + return hqos_k, sqos_k, kvecs_d + # kpoints extrapolation: QO validation by reproducing band structure + # Rlist = np.loadtxt(path+'/QO_supercells.dat', dtype=np.int32) + # kpath = np.loadtxt('./band/kpoints')[:,1:4] + # print("There are", len(Rlist), "supercells") + + # hk, sk = k_extrapolation(hqos_k, kvecs_d, Rlist, kpath), k_extrapolation(sqos_k, kvecs_d, Rlist, kpath) + + # eigs = [] + # for ik in range(len(kpath)): + # eigs.append(bad_condition_diag(hk[ik], sk[ik], diag_thr)) + # eigs = np.asarray(eigs) + + # band = np.loadtxt('./band/OUT.ABACUS/BANDS_1.dat') + # figure_title = 'Band Structure of QO Hamiltonian: KPT' + str(ndim) + 'x' + str(ndim) + 'x' + str(ndim) + # import matplotlib.pyplot as plt + # plt.plot(band[:, 1], eigs[:,:]*13.6, 'o', markerfacecolor='none', markeredgecolor='b') + # plt.plot(band[:, 1], band[:, 2:]) + # plt.ylim([-7.5,25]) + # plt.title(figure_title) + # plt.savefig(fpic) + +import unittest +class QOUnittest(unittest.TestCase): + def test_read_mat_hs(self): + for ik in range(ndim**3): + print("Unittest test_read_mat_hs on kpoint", ik) + fhs = path + f'/data-{ik}-S' + mat = read_mat_hs(fhs) + nrows, ncols = mat.shape + # check square matrix + self.assertEqual(nrows, ncols) + # check hermitian + self.assertListEqual(mat.tolist(), mat.conj().T.tolist()) + # check not-all-zero + self.assertNotEqual(np.sum(np.abs(mat)), 0) + + def test_read_lowf_k(self): + for ik in range(ndim**3): + print("Unittest test_read_lowf_k on kpoint", ik) + flowf = path + f"/LOWF_K_{ik+1}.txt" + lowf_k, kvec_d = read_lowf(flowf) + nlocal, nband = lowf_k.shape + # check number of bands + self.assertGreater(nband, 0) + # check number of local orbitals + self.assertGreater(nlocal, 0) + # check dimension of space + self.assertGreater(nlocal, nband) + + fs = path + f'/data-{ik}-S' + s_k = read_mat_hs(fs) + nrows, ncols = s_k.shape + self.assertEqual(nrows, ncols) # should be a square matrix + self.assertEqual(nrows, nlocal) # should be overlap matrix of basis functions + _one = lowf_k.conj().T@s_k@lowf_k + # orthonormality + for i in range(nband): + for j in range(nband): + if i != j: + self.assertAlmostEqual(abs(_one[i, j].real), 0, delta=1e-5) + self.assertAlmostEqual(abs(_one[i, j].imag), 0, delta=1e-5) + else: + self.assertAlmostEqual(abs(_one[i, j].real), 1, delta=1e-5) + self.assertAlmostEqual(abs(_one[i, j].imag), 0, delta=1e-5) + + def test_read_ao_proj(self): + for ik in range(ndim**3): + print("Unittest test_read_ao_proj on kpoint", ik) + fao_proj = path + f'/QO_ovlp_{ik}.dat' + ao_proj, kvec_d = read_ao_proj(fao_proj) + nrows, ncols = ao_proj.shape + self.assertGreater(nrows, 0) + self.assertGreater(ncols, 0) + self.assertGreater(nrows, ncols) # space dimension comparision + + def test_cal_denmat(self): + for ik in range(ndim**3): + print("Unittest test_cal_denmat on kpoint", ik) + flowf = path + f"/LOWF_K_{ik+1}.txt" + lowf_k, kvec_d = read_lowf(flowf) + denmat_k = cal_denmat(lowf_k) + + nrows, ncols = denmat_k.shape + self.assertEqual(nrows, ncols) + nlocal, nband = lowf_k.shape + self.assertEqual(nrows, nlocal) + self.assertEqual(ncols, nlocal) + + denmat_dagger_k = denmat_k.conj().T + # D = D^+ + for i in range(nrows): + for j in range(ncols): + self.assertAlmostEqual(denmat_k[i, j], denmat_dagger_k[i, j], delta=1e-5) + + fs = path + f'/data-{ik}-S' + s_k = read_mat_hs(fs) + D_k = denmat_k.tolist() + DSD_k = (denmat_k@s_k@denmat_k).tolist() + self.assertEqual(len(D_k), len(DSD_k)) + self.assertEqual(len(D_k[0]), len(DSD_k[0])) + # D = DSD + for i in range(len(D_k)): + for j in range(len(D_k[0])): + self.assertAlmostEqual(D_k[i][j], DSD_k[i][j], delta=1e-3) + + def test_lowdin_onW(self): + for ik in range(ndim**3): + print("Unittest test_lowdin_onW on kpoint", ik) + ndim_test = 2 + flowf = path + f"/LOWF_K_{ik+1}.txt" + lowf_k, kvec_d = read_lowf(flowf) + denmat_k = cal_denmat(lowf_k) + fs = path + f'/data-{ik}-S' + s_k = read_mat_hs(fs) + fao_proj = path + f'/QO_ovlp_{ik}.dat' + ao_proj_k, kvec_d = read_ao_proj(fao_proj) + + lowf_bar = lowdin_onW(s_k=s_k, + denmat_k=denmat_k, + ao_proj_k=ao_proj_k, + m = ndim_test) + _one = lowf_bar.conj().T@s_k@lowf_bar + self.assertEqual(_one.shape, (ndim_test, ndim_test)) + # orthonormality + for i in range(ndim_test): + for j in range(ndim_test): + if i != j: + self.assertAlmostEqual(_one[i, j], 0, delta=1e-6) + else: + self.assertAlmostEqual(_one[i, j], 1, delta=1e-6) + # D = DSD + lowf_tilde_k = np.concatenate((lowf_k, lowf_bar), axis=1) + denmat_tilde_k = cal_denmat(lowf_tilde_k) + D_tilde_k = denmat_tilde_k.tolist() + DSD_tilde_k = (denmat_tilde_k@s_k@denmat_tilde_k).tolist() + for i in range(len(D_tilde_k)): + for j in range(len(D_tilde_k[0])): + self.assertAlmostEqual(D_tilde_k[i][j], DSD_tilde_k[i][j], delta=1e-5) + + def test_cal_tb_hs(self): + print("Unittest test_cal_tb_hs") + # eye + nbasis = 16 + I = np.eye(nbasis, dtype=np.complex128) + # CASE 1: the most ideal case + # create an Hermite matrix as model Hamiltonian + H = np.random.rand(nbasis, nbasis) + 1.0j*np.random.rand(nbasis, nbasis) + H = H + H.conj().T + # eye as model overlap matrix + S = I + # eye as model denmat + D = I + # eye as model ao_proj + A = I + # test + Htb, Stb = cal_tb_hs(H, S, D, A) + # check shape + self.assertEqual(Htb.shape, (nbasis, nbasis)) + self.assertEqual(Stb.shape, (nbasis, nbasis)) + # check hermitian + self.assertListEqual(Htb.tolist(), Htb.conj().T.tolist()) + self.assertListEqual(Stb.tolist(), Stb.conj().T.tolist()) + # check not-all-zero + self.assertNotEqual(np.sum(np.abs(Htb)), 0) + self.assertNotEqual(np.sum(np.abs(Stb)), 0) + # CASE 2: less ideal case + # directly truncate A, for the meaning the AO basis is not as many as original basis + nAO = 3 + A = A[:, :nAO] + D = I[:, :nAO]@I[:nAO, :] # usually D is not a full-rank matrix + Htb, Stb = cal_tb_hs(H, S, D, A) + # check shape + self.assertEqual(Htb.shape, (nAO, nAO)) + self.assertEqual(Stb.shape, (nAO, nAO)) + # check hermitian + self.assertListEqual(Htb.tolist(), Htb.conj().T.tolist()) + self.assertListEqual(Stb.tolist(), Stb.conj().T.tolist()) + # check not-all-zero + self.assertNotEqual(np.sum(np.abs(Htb)), 0) + self.assertNotEqual(np.sum(np.abs(Stb)), 0) + + def test_kR_transform(self): + for nk in [1, 3, 5]: + kpoints = it.product(np.arange(-0.5+0.5/nk, 0.5, 1.0/nk), repeat=3) + kpoints = np.array(list(kpoints)) + Rs = it.product(range(-int((nk-1)/2), int((nk+1)/2)), repeat=3) + Rs = np.array(list(Rs)) + assert len(kpoints) == len(Rs) + ndim = 20 # square matrix dimension + mats_0 = [np.random.rand(ndim, ndim) for i in range(len(kpoints))] + mats_k = [kR_transform(kpoints[i], mats_0, Rs, "R->k") for i in range(len(kpoints))] + mats_kR = [kR_transform(Rs[i], mats_k, kpoints, "k->R") for i in range(len(kpoints))] + for ik in range(len(kpoints)): + for i in range(ndim): + for j in range(ndim): + self.assertAlmostEqual(mats_kR[ik][i, j].real, mats_0[ik][i, j].real, delta=1e-12) + self.assertAlmostEqual(mats_kR[ik][i, j].imag, mats_0[ik][i, j].imag, delta=1e-12) + + def test_k_extrapolation(self): + for nk in [1, 3, 5]: + kpoints = it.product(np.arange(-0.5+0.5/nk, 0.5, 1.0/nk), repeat=3) + kpoints = np.array(list(kpoints)) + Rs = it.product(range(-int((nk-1)/2), int((nk+1)/2)), repeat=3) + Rs = np.array(list(Rs)) + assert len(kpoints) == len(Rs) + ndim = 20 # square matrix dimension + mats_0 = [np.random.rand(ndim, ndim) for i in range(len(kpoints))] + mats_k2R2k = k_extrapolation(mats_0, kpoints, Rs, kpoints) + for ik in range(len(kpoints)): + for i in range(ndim): + for j in range(ndim): + self.assertAlmostEqual(mats_k2R2k[ik][i, j].real, mats_0[ik][i, j].real, delta=1e-12) + self.assertAlmostEqual(mats_k2R2k[ik][i, j].imag, mats_0[ik][i, j].imag, delta=1e-12) + +if __name__ == "__main__": + + # False to run unittest, True to run production + run_type = "production" # can be "unittest" or "production" + ndim = 7 # dimension of Monkhorst-Pack grid + qo_basis = "hydrogen" + qo_strategy = "minimal-valence" + path = f"./scf/{qo_basis}/{qo_strategy}/OUT.ABACUS-{ndim}{ndim}{ndim}/" + #fpic = f"{qo_basis}-{qo_strategy}-{ndim}{ndim}{ndim}.png" + # only for production + eig_thr = 1e-10 + ib_min, ib_max = 0, 4 + + if run_type == "production": + nkpts = ndim*ndim*ndim + hqos_k, sqos_k, kvecs_d = production(path, nkpts, ib_min, ib_max, eig_thr, error_estimation_with_lcao=True) + + elif run_type == "unittest": + unittest.main() diff --git a/tools/qo/setup.py b/tools/qo/setup.py deleted file mode 100644 index 53dd047d5b..0000000000 --- a/tools/qo/setup.py +++ /dev/null @@ -1,26 +0,0 @@ -import pathlib -import setuptools - - -here = pathlib.Path(__file__).parent.resolve() -readme = (here / 'README.md').read_text(encoding='utf-8') - -# did not include torch and pyscf here -install_requires=["scipy", "numpy"] - -setuptools.setup( - name="abacus2qo", - author="Yike HUANG", - author_email="huangyk@aici.ac.cn", - description="Quasiatomic Orbital (QO) analysis", - long_description=readme, - long_description_content_type="text/markdown", - packages=setuptools.find_packages(include=['abacus2qo', 'abacus2qo.*']), - classifiers=[ - 'Development Status :: 3 - Alpha', - "Programming Language :: Python :: 3.7", - ], - install_requires=install_requires, - python_requires=">=3.7", - version="0.0.2" -) \ No newline at end of file