Merge pull request electronic-structure#488 from toxa81/develop

Latest updates

.github/workflows/build.yml

@@ -0,0 +1,28 @@
+name: Build
+
+on:
+  - push
+  - pull_request
+
+jobs:
+  build:
+    runs-on: ubuntu-latest
+    steps:
+      - uses: actions/checkout@v2
+        with:
+          submodules: recursive
+
+      - name: Install tools
+        run: sudo apt-get update && sudo apt-get install -y gfortran mpich libmpich-dev liblapack-dev libhdf5-dev libgsl-dev libfftw3-dev libxc-dev
+
+      - name: Build dependencies
+        run: ./ci/github_build_dependencies.sh
+
+      - name: Build SIRIUS
+        run: ./ci/github_build_sirius.sh
+
+      - name: Run unit tests
+        run: ./ci/github_run_unit_tests.sh
+
+      - name: Run verification with ReFrame
+        run: ./ci/github_run_verification.sh

.github/workflows/check-format.yml

@@ -0,0 +1,15 @@
+name: Check source code format
+
+on:
+  - push
+  - pull_request
+
+jobs:
+  check:
+    runs-on: ubuntu-latest
+    steps:
+      - uses: actions/checkout@v2
+        with:
+          submodules: recursive
+      - name: Check
+        run: ./check_format.x

CMakeLists.txt

@@ -1,6 +1,6 @@
 cmake_minimum_required(VERSION 3.12)
 
-project(SIRIUS VERSION 6.5.2)
+project(SIRIUS VERSION 6.5.3)
 
 # set language and standard
 enable_language(CXX Fortran)

apps/tests/test_davidson.cpp

@@ -171,7 +171,7 @@ void test_davidson(cmd_args const& args__)
 
     for (int r = 0; r < 2; r++) {
         double vk[] = {0.1, 0.1, 0.1};
-        K_point kp(ctx, vk, 1.0);
+        K_point kp(ctx, vk, 1.0, 0);
         kp.initialize();
         std::cout << "num_gkvec=" << kp.num_gkvec() << "\n";
         for (int i = 0; i < ctx.num_bands(); i++) {

ci/github_build_dependencies.sh

@@ -0,0 +1,4 @@
+#!/bin/bash
+export LIBSPGROOT=$HOME/local
+python3 prerequisite.py $HOME/local spfft spg
+

ci/github_build_sirius.sh

@@ -0,0 +1,7 @@
+#!/bin/bash
+
+mkdir build
+cd build
+cmake ../ -DCMAKE_BUILD_TYPE=RelWithDebInfo -DSpFFT_DIR=$HOME/local/lib/cmake/SpFFT -DCMAKE_INSTALL_PREFIX=$HOME/local -DBUILD_TESTS=1
+make -j 2 install
+

ci/github_run_unit_tests.sh

@@ -0,0 +1,3 @@
+#!/bin/bash
+cd build/apps/unit_tests
+./unit_tests.x

ci/github_run_verification.sh

@@ -0,0 +1,6 @@
+#!/bin/bash
+
+export PATH=$HOME/local/bin:$HOME/reframe/bin:$PATH
+git clone https://github.com/eth-cscs/reframe.git $HOME/reframe
+reframe -C ./reframe/config.py --system=linux:cpu -c ./reframe/checks -R -r --tag serial --exec-policy=serial
+

examples/pp-pw/QE/LiF/pw.nc.gga.in

@@ -5,8 +5,8 @@
     outdir='./',
     prefix = 'LiF_',
     verbosity='high',
-    tstress = false
-    tprnfor = false
+    tstress = true
+    tprnfor = true
     nstep=2
  /
  &system

examples/pp-pw/QE/LiF/pw.paw.lda.in

@@ -5,8 +5,8 @@
     outdir='./',
     prefix = 'LiF_',
     verbosity='high',
-    tstress = false
-    tprnfor = false
+    tstress = true
+    tprnfor = true
     nstep=2
  /
  &system
@@ -15,10 +15,10 @@
     occupations = 'smearing', smearing = 'gauss', degauss = 0.02
  /
  &electrons
-    conv_thr =  1.0d-8,
+    conv_thr =  1.0d-11,
     mixing_beta = 0.7,
-    diago_thr_init=1e-2
     electron_maxstep=100
+    diago_full_acc=.true.
  /
 &IONS
  ion_dynamics='bfgs',

src/SDDK/gvec.hpp

@@ -34,7 +34,7 @@
 #include "geometry3d.hpp"
 #include "serializer.hpp"
 #include "splindex.hpp"
-#include "../utils/profiler.hpp"
+#include "utils/profiler.hpp"
 
 using namespace geometry3d;
 

src/band/band.cpp

@@ -204,43 +204,46 @@ void Band::initialize_subspace(Hamiltonian_k& Hk__, int num_ao__) const
         phi.pw_coeffs(ispn).prime().zero();
     }
 
-    PROFILE_START("sirius::Band::initialize_subspace|kp|wf");
-
     /* generate the initial atomic wavefunctions */
     std::vector<int> atoms(ctx_.unit_cell().num_atoms());
     std::iota(atoms.begin(), atoms.end(), 0);
     Hk__.kp().generate_atomic_wave_functions(atoms, [&](int iat){return &ctx_.unit_cell().atom_type(iat).indexb_wfs();},
                                              ctx_.atomic_wf_ri(), phi);
 
+    /* generate some random noise */
+    std::vector<double> tmp(4096);
+    for (int i = 0; i < 4096; i++) {
+        tmp[i] = 1e-5 * utils::random<double>();
+    }
+    PROFILE_START("sirius::Band::initialize_subspace|kp|wf");
     /* fill remaining wave-functions with pseudo-random guess */
     assert(Hk__.kp().num_gkvec() > num_phi + 10);
-    #pragma omp parallel for schedule(static)
-    for (int i = 0; i < num_phi - num_ao__; i++) {
-        for (int igk_loc = 0; igk_loc < Hk__.kp().num_gkvec_loc(); igk_loc++) {
-            /* global index of G+k vector */
-            int igk = Hk__.kp().idxgk(igk_loc);
-            if (igk == i + 1) {
-                phi.pw_coeffs(0).prime(igk_loc, num_ao__ + i) = 1.0;
-            }
-            if (igk == i + 2) {
-                phi.pw_coeffs(0).prime(igk_loc, num_ao__ + i) = 0.5;
-            }
-            if (igk == i + 3) {
-                phi.pw_coeffs(0).prime(igk_loc, num_ao__ + i) = 0.25;
+    #pragma omp parallel
+    {
+        for (int i = 0; i < num_phi - num_ao__; i++) {
+            #pragma omp for schedule(static) nowait
+            for (int igk_loc = 0; igk_loc < Hk__.kp().num_gkvec_loc(); igk_loc++) {
+                /* global index of G+k vector */
+                int igk = Hk__.kp().idxgk(igk_loc);
+                if (igk == i + 1) {
+                    phi.pw_coeffs(0).prime(igk_loc, num_ao__ + i) = 1.0;
+                }
+                if (igk == i + 2) {
+                    phi.pw_coeffs(0).prime(igk_loc, num_ao__ + i) = 0.5;
+                }
+                if (igk == i + 3) {
+                    phi.pw_coeffs(0).prime(igk_loc, num_ao__ + i) = 0.25;
+                }
             }
         }
-    }
-
-    std::vector<double> tmp(4096);
-    for (int i = 0; i < 4096; i++) {
-        tmp[i] = utils::random<double>();
-    }
-    #pragma omp parallel for schedule(static)
-    for (int i = 0; i < num_phi; i++) {
-        for (int igk_loc = Hk__.kp().gkvec().skip_g0(); igk_loc < Hk__.kp().num_gkvec_loc(); igk_loc++) {
-            /* global index of G+k vector */
-            int igk = Hk__.kp().idxgk(igk_loc);
-            phi.pw_coeffs(0).prime(igk_loc, i) += tmp[igk & 0xFFF] * 1e-5;
+        /* add random noise */
+        for (int i = 0; i < num_phi; i++) {
+            #pragma omp for schedule(static) nowait
+            for (int igk_loc = Hk__.kp().gkvec().skip_g0(); igk_loc < Hk__.kp().num_gkvec_loc(); igk_loc++) {
+                /* global index of G+k vector */
+                int igk = Hk__.kp().idxgk(igk_loc);
+                phi.pw_coeffs(0).prime(igk_loc, i) += tmp[igk & 0xFFF];
+            }
         }
     }
 

src/band/diag_pseudo_potential.cpp

@@ -511,6 +511,13 @@ Band::diag_pseudo_potential_davidson(Hamiltonian_k& Hk__) const
                     kp.message(2, __function_name__, "%s", "wave-functions are not recomputed\n");
                 }
 
+                if (k == (itso.num_steps_ - 1) && n > itso.min_num_res_) {
+                    std::stringstream s;
+                    s << "[sirius::Band::diag_pseudo_potential_davidson] maximum number of iterations reached, but " <<
+                         n << " residual(s) did not converge for k-point " << kp.vk();
+                    WARNING(s);
+                }
+
                 /* exit the loop if the eigen-vectors are converged or this is a last iteration */
                 if (n <= itso.min_num_res_ || k == (itso.num_steps_ - 1)) {
                     kp.message(3, __function_name__, "end of iterative diagonalization; n=%i, k=%i\n", n, k);

src/k_point/k_point.cpp

@@ -977,10 +977,9 @@ K_point::generate_atomic_wave_functions(std::vector<int> atoms__,
         }
     }
 
+    PROFILE_START("sirius::K_point::generate_atomic_wave_functions|wft");
     #pragma omp parallel for schedule(static)
     for (int igk_loc = 0; igk_loc < this->num_gkvec_loc(); igk_loc++) {
-        /* global index of G+k vector */
-        //int igk = this->idxgk(igk_loc);
         /* vs = {r, theta, phi} */
         auto vs = geometry3d::spherical_coordinates(this->gkvec().gkvec_cart<index_domain_t::local>(igk_loc));
 
@@ -1014,12 +1013,14 @@ K_point::generate_atomic_wave_functions(std::vector<int> atoms__,
             }
         }
     }
+    PROFILE_STOP("sirius::K_point::generate_atomic_wave_functions|wft");
 
     for (int ia: atoms__) {
 
         double phase = twopi * dot(gkvec().vk(), unit_cell_.atom(ia).position());
         double_complex phase_k = std::exp(double_complex(0.0, phase));
 
+        PROFILE_START("sirius::K_point::generate_atomic_wave_functions|1");
         /* quickly compute phase factors without calling exp() function */
         std::vector<double_complex> phase_gk(num_gkvec_loc());
         #pragma omp parallel for schedule(static)
@@ -1030,13 +1031,18 @@ K_point::generate_atomic_wave_functions(std::vector<int> atoms__,
             /* total phase e^{-i(G+k)r_{\alpha}} */
             phase_gk[igk_loc] = std::conj(ctx_.gvec_phase_factor(G, ia) * phase_k);
         }
+        PROFILE_STOP("sirius::K_point::generate_atomic_wave_functions|1");
+
+        PROFILE_START("sirius::K_point::generate_atomic_wave_functions|2");
         int iat = unit_cell_.atom(ia).type_id();
-        #pragma omp parallel for
+        #pragma omp parallel
         for (int xi = 0; xi < indexb__(iat)->size(); xi++) {
+            #pragma omp for schedule(static) nowait
             for (int igk_loc = 0; igk_loc < num_gkvec_loc(); igk_loc++) {
                 wf__.pw_coeffs(0).prime(igk_loc, offset[ia] + xi) = wf_t[iat](igk_loc, xi) * phase_gk[igk_loc];
             }
         }
+        PROFILE_STOP("sirius::K_point::generate_atomic_wave_functions|2");
     }
 }
 

src/k_point/k_point.hpp

@@ -43,8 +43,11 @@ class K_point
     /// Unit cell object.
     Unit_cell const& unit_cell_;
 
+    /// K-point id.
+    int id_{-1};
+
     /// Weight of k-point.
-    double weight_;
+    double weight_{1.0};
 
     /// Fractional k-point coordinates.
     vector3d<double> vk_;
@@ -195,9 +198,10 @@ class K_point
 
   public:
     /// Constructor
-    K_point(Simulation_context& ctx__, double const* vk__, double weight__)
+    K_point(Simulation_context& ctx__, double const* vk__, double weight__, int id__)
         : ctx_(ctx__)
         , unit_cell_(ctx_.unit_cell())
+        , id_(id__)
         , weight_(weight__)
         , comm_(ctx_.comm_band())
         , comm_row_(ctx_.blacs_grid().comm_row())

src/k_point/k_point_set.cpp

@@ -18,8 +18,10 @@
 // OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
 
 #include <limits>
+#include "dft/smearing.hpp"
 #include "k_point/k_point.hpp"
 #include "k_point/k_point_set.hpp"
+#include "symmetry/get_irreducible_reciprocal_mesh.hpp"
 
 namespace sirius {
 
@@ -58,7 +60,7 @@ void K_point_set::create_k_mesh(vector3d<int> k_grid__, vector3d<int> k_shift__,
     mdarray<double, 2> kp;
     std::vector<double> wk;
     if (use_symmetry__) {
-        auto result = get_irreducible_reciprocal_mesh(unit_cell_.symmetry(), k_grid__, k_shift__);
+        auto result = get_irreducible_reciprocal_mesh(ctx_.unit_cell().symmetry(), k_grid__, k_shift__);
         nk          = std::get<0>(result);
         wk          = std::get<1>(result);
         auto tmp    = std::get<2>(result);
@@ -174,7 +176,7 @@ void K_point_set::find_band_occupancies()
     double ne{0};
 
     /* target number of electrons */
-    double ne_target = unit_cell_.num_valence_electrons() - ctx_.parameters_input().extra_charge_;
+    double ne_target = ctx_.unit_cell().num_valence_electrons() - ctx_.parameters_input().extra_charge_;
 
     if (std::abs(ctx_.num_fv_states() * double(ctx_.max_occupancy()) - ne_target) < 1e-10) {
         // this is an insulator, skip search for band occupancies