Update mater branch (electronic-structure#799)

* [ELPA] Fit latest elpa into the code (electronic-structure#798)

Switch to latest elpa-2022.11.001.rc2

* [refactor] Change geometry3d -> r3 namespace (electronic-structure#793)

No change in the code logic, just renaming namespace and classes.

CMakeLists.txt

@@ -1,6 +1,6 @@
 cmake_minimum_required(VERSION 3.18)
 
-project(SIRIUS VERSION 7.4.0)
+project(SIRIUS VERSION 7.4.1)
 
 # user variables
 set(CREATE_PYTHON_MODULE    OFF CACHE BOOL "create sirius Python module")

apps/dft_loop/sirius.scf.cpp

@@ -476,13 +476,13 @@ void run_tasks(cmd_args const& args)
 
         double t{0};
         for (size_t i = 0; i < vertex.size() - 1; i++) {
-            vector3d<double> v0 = vector3d<double>(vertex[i].second);
-            vector3d<double> v1 = vector3d<double>(vertex[i + 1].second);
-            vector3d<double> dv = v1 - v0;
-            vector3d<double> dv_cart = dot(ctx->unit_cell().reciprocal_lattice_vectors(), dv);
+            r3::vector<double> v0 = r3::vector<double>(vertex[i].second);
+            r3::vector<double> v1 = r3::vector<double>(vertex[i + 1].second);
+            r3::vector<double> dv = v1 - v0;
+            r3::vector<double> dv_cart = dot(ctx->unit_cell().reciprocal_lattice_vectors(), dv);
             int np = std::max(10, static_cast<int>(30 * dv_cart.length()));
             for (int j = 1; j <= np; j++) {
-                vector3d<double> v = v0 + dv * static_cast<double>(j) / np;
+                r3::vector<double> v = v0 + dv * static_cast<double>(j) / np;
                 ks.add_kpoint(&v[0], 1.0);
                 t += dv_cart.length() / np;
                 x_axis.push_back(t);

apps/nlcg/sirius.nlcg.cpp

@@ -1,8 +1,8 @@
-#include "utils/profiler.hpp"
 #include <sirius.hpp>
 #include <utils/json.hpp>
-#include "nlcglib/adaptor.hpp"
+#include <nlcglib/adaptor.hpp>
 #include <nlcglib/nlcglib.hpp>
+#include "utils/profiler.hpp"
 
 using namespace sirius;
 using json = nlohmann::json;
@@ -30,15 +30,14 @@ std::unique_ptr<Simulation_context> create_sim_ctx(std::string     fname__,
 
     auto& inp = ctx.cfg().parameters();
     if (inp.gamma_point() && !(inp.ngridk()[0] * inp.ngridk()[1] * inp.ngridk()[2] == 1)) {
-        TERMINATE("this is not a Gamma-point calculation")
+        RTE_THROW("this is not a Gamma-point calculation")
     }
 
     ctx.import(args__);
 
     return ctx_ptr;
 }
 
-
 double ground_state(Simulation_context& ctx,
                     task_t              task,
                     cmd_args const&     args,

apps/tests/test_davidson.cpp

@@ -134,12 +134,12 @@ void test_davidson(cmd_args const& args__)
         json_conf["parameters"]["num_bands"] = num_bands;
     }
 
-    std::vector<geometry3d::vector3d<double>> coord;
+    std::vector<r3::vector<double>> coord;
     double p = 1.0 / N;
     for (int i = 0; i < N; i++) {
         for (int j = 0; j < N; j++) {
             for (int k = 0; k < N; k++) {
-                coord.push_back(geometry3d::vector3d<double>(i * p, j * p, k * p));
+                coord.push_back(r3::vector<double>(i * p, j * p, k * p));
             }
         }
     }

apps/tests/test_eigen.cpp

@@ -48,7 +48,7 @@ test_diag(sddk::BLACS_grid const& blacs_grid__, int N__, int n__, int nev__, int
         if (std::is_same<T, double>::value) {
             printf("real data type\n");
         }
-        if (std::is_same<T, double_complex>::value) {
+        if (std::is_same<T, std::complex<double>>::value) {
             printf("complex data type\n");
         }
     }
@@ -138,7 +138,7 @@ void test_diag2(sddk::BLACS_grid const& blacs_grid__,
 {
     auto solver = Eigensolver_factory(name__);
 
-    sddk::matrix<double_complex> full_mtrx;
+    sddk::matrix<std::complex<double>> full_mtrx;
     int n;
     if (blacs_grid__.comm().rank() == 0) {
         sddk::HDF5_tree h5(fname__, sddk::hdf5_access_t::read_only);
@@ -148,22 +148,22 @@ void test_diag2(sddk::BLACS_grid const& blacs_grid__,
         if (n != m) {
             TERMINATE("not a square matrix");
         }
-        full_mtrx = sddk::matrix<double_complex>(n, n);
+        full_mtrx = sddk::matrix<std::complex<double>>(n, n);
         h5.read("/mtrx", full_mtrx);
         blacs_grid__.comm().bcast(&n, 1, 0);
         blacs_grid__.comm().bcast(full_mtrx.at(sddk::memory_t::host), static_cast<int>(full_mtrx.size()), 0);
     } else {
         blacs_grid__.comm().bcast(&n, 1, 0);
-        full_mtrx = sddk::matrix<double_complex>(n, n);
+        full_mtrx = sddk::matrix<std::complex<double>>(n, n);
         blacs_grid__.comm().bcast(full_mtrx.at(sddk::memory_t::host), static_cast<int>(full_mtrx.size()), 0);
     }
     if (blacs_grid__.comm().rank() == 0) {
         printf("matrix size: %i\n", n);
     }
 
     std::vector<double> eval(n);
-    sddk::dmatrix<double_complex> A(n, n, blacs_grid__, bs__, bs__);
-    sddk::dmatrix<double_complex> Z(n, n, blacs_grid__, bs__, bs__);
+    sddk::dmatrix<std::complex<double>> A(n, n, blacs_grid__, bs__, bs__);
+    sddk::dmatrix<std::complex<double>> Z(n, n, blacs_grid__, bs__, bs__);
 
     for (int j = 0; j < A.num_cols_local(); j++) {
         for (int i = 0; i < A.num_rows_local(); i++) {
@@ -199,7 +199,7 @@ void call_test(std::vector<int> mpi_grid__,
             if (type__ == 0) {
                 t = test_diag<double>(blacs_grid, N__, n__, nev__, bs__, test_gen__, name__, *solver);
             } else {
-                t = test_diag<double_complex>(blacs_grid, N__, n__, nev__, bs__, test_gen__, name__, *solver);
+                t = test_diag<std::complex<double>>(blacs_grid, N__, n__, nev__, bs__, test_gen__, name__, *solver);
             }
             /* skip first "warmup" measurment */
             if (i) {

apps/tests/test_examples.cpp

@@ -12,7 +12,7 @@ if (Communicator::world().rank() == 0) {
 }
 /* reciprocal lattice vectors in
    inverse atomic units */
-matrix3d<double> M = {{1, 0, 0},
+r3::matrix<double> M = {{1, 0, 0},
                       {0, 1, 0},
                       {0, 0, 1}};
 /* G-vector cutoff radius in
@@ -47,7 +47,7 @@ if (Communicator::world().rank() == 0) {
 }
 /* reciprocal lattice vectors in 
     inverse atomic units */
-matrix3d<double> M = {{1, 0, 0},
+r3::matrix<double> M = {{1, 0, 0},
                       {0, 1, 0},
                       {0, 0, 1}};
 /* G-vector cutoff radius in
@@ -86,9 +86,9 @@ spfft_transform_type<double> spfft(
 
 /* create data buffer with local number of G-vectors
    and fill with random numbers */
-mdarray<double_complex, 1> f(gvp.gvec_count_fft());
+mdarray<std::complex<double>, 1> f(gvp.gvec_count_fft());
 f = [](int64_t){
-  return utils::random<double_complex>();
+  return utils::random<std::complex<double>>();
 };
 /* transform to real space */
 spfft.backward(reinterpret_cast<double const*>(
@@ -119,7 +119,7 @@ spla::Context spla_ctx(SPLA_PU_HOST);
 
 /* reciprocal lattice vectors in 
    inverse atomic units */
-matrix3d<double> M = {{1, 0, 0},
+r3::matrix<double> M = {{1, 0, 0},
                       {0, 1, 0},
                       {0, 0, 1}};
 /* G-vector cutoff radius in
@@ -138,15 +138,15 @@ wf::Wave_functions<double> wf(gvec, wf::num_mag_dims(0), wf::num_bands(N), memor
 for (int i = 0; i < N; i++) {
     for (int j = 0; j < gvec->count(); j++) {
         wf.pw_coeffs(j, wf::spin_index(0), wf::band_index(i)) =
-            utils::random<double_complex>();
+            utils::random<std::complex<double>>();
     }
 }
 /* create a 2x2 BLACS grid */
 BLACS_grid grid(Communicator::world(), 2, 2);
 /* cyclic block size */
 int bs = 16;
 /* create a distributed overlap matrix */
-dmatrix<double_complex> o(N, N, grid, bs, bs);
+dmatrix<std::complex<double>> o(N, N, grid, bs, bs);
 /* create temporary wave-functions */
 wf::Wave_functions<double> tmp(gvec, wf::num_mag_dims(0), wf::num_bands(N), memory_t::host);
 /* orthogonalize wave-functions */

apps/tests/test_gemm.cpp

@@ -5,7 +5,7 @@
 typedef double gemm_type;
 int const nop_gemm = 2;
 #else
-typedef double_complex gemm_type;
+typedef std::complex<double> gemm_type;
 int const nop_gemm = 8;
 #endif
 

apps/tests/test_hdf5.cpp

@@ -36,7 +36,7 @@ void test2()
     f["node1"].write("md1", md1);
     f["node1"].write(0, md1);
 
-    sddk::mdarray<double_complex, 2> md2(2, 4);
+    sddk::mdarray<std::complex<double>, 2> md2(2, 4);
     md2.zero();
     f["node1"].write("md2", md2);
     f["node1"].write(1, md2);
@@ -55,7 +55,7 @@ void test3()
     f["node1"].read("md1", md1);
     f["node1"].read(0, md1);
 
-    sddk::mdarray<double_complex, 2> md2(2, 4);
+    sddk::mdarray<std::complex<double>, 2> md2(2, 4);
     md2.zero();
     f["node1"].read("md2", md2);
     f["node1"].read(1, md2);

apps/tests/test_hloc.cpp

@@ -117,7 +117,7 @@ int main(int argn, char** argv)
         json_conf["parameters"]["gamma_point"] = reduce_gvec;
 
         auto ctx = sirius::create_simulation_context(json_conf, {{10, 0, 0}, {0, 10, 0}, {0, 0, 10}}, 0,
-            std::vector<geometry3d::vector3d<double>>(), false, false);
+            std::vector<r3::vector<double>>(), false, false);
         for (int i = 0; i < repeat; i++) {
             if (fp32) {
 #if defined(USE_FP32)

apps/tests/test_mdarray.cpp

@@ -52,8 +52,8 @@ void f3()
         #pragma omp parallel
         {
             int tid = omp_get_thread_num();
-            sddk::mdarray<double_complex, 2> a(100, 100);
-            a(0, 0) = double_complex(tid, tid);
+            sddk::mdarray<std::complex<double>, 2> a(100, 100);
+            a(0, 0) = std::complex<double>(tid, tid);
         }
     }
 }

apps/tests/test_phase.cpp

@@ -5,28 +5,28 @@ using namespace sirius;
 void test1()
 {
     int N = 30000000;
-    std::vector<vector3d<double>> a(N);
+    std::vector<r3::vector<double>> a(N);
 
     for (int i = 0; i < N; i++) {
         double r = utils::random<double>();
         a[i] = {r, r, r};
     }
-    std::vector<double_complex> phase(N, 0);
+    std::vector<std::complex<double>> phase(N, 0);
     double t1{0};
     double t2{0};
 
     for (int i = 0; i < 10; i++) {
         double t = -omp_get_wtime();
         #pragma omp parallel for
         for (int i = 0; i < N; i++) {
-            phase[i] = std::exp(double_complex(0, dot(a[i], a[i])));
+            phase[i] = std::exp(std::complex<double>(0, dot(a[i], a[i])));
         }
         t1 += (t + omp_get_wtime());
 
         t = -omp_get_wtime();
         #pragma omp parallel for schedule(static)
         for (int i = 0; i < N; i++) {
-            phase[i] = std::exp(double_complex(0, dot(a[i], a[i])));
+            phase[i] = std::exp(std::complex<double>(0, dot(a[i], a[i])));
         }
         t2 += (t + omp_get_wtime());
     }

apps/tests/test_pppw_xc.cpp

@@ -32,11 +32,11 @@ void test_davidson(cmd_args const& args__)
     json_conf["control"]["mpi_grid_dims"] = mpi_grid;
 
     double p = 1.0 / N;
-    std::vector<geometry3d::vector3d<double>> coord;
+    std::vector<r3::vector<double>> coord;
     for (int i = 0; i < N; i++) {
         for (int j = 0; j < N; j++) {
             for (int k = 0; k < N; k++) {
-                coord.push_back(geometry3d::vector3d<double>(i * p, j * p, k * p));
+                coord.push_back(r3::vector<double>(i * p, j * p, k * p));
             }
         }
     }

apps/tests/test_reduce.cpp

@@ -9,8 +9,8 @@ int main(int argn, char** argv)
 
     bool pr = (Communicator::world().rank() == 0);
 
-    sddk::mdarray<double_complex, 2> tmp_d(100, 200);
-    sddk::mdarray<double_complex, 2> tmp(100, 200);
+    sddk::mdarray<std::complex<double>, 2> tmp_d(100, 200);
+    sddk::mdarray<std::complex<double>, 2> tmp(100, 200);
 
     std::cout << Communicator::world().rank() << " tmp_d " << tmp_d.at(memory_t::host) << std::endl;
     std::cout << Communicator::world().rank() << " tmp " << tmp.at(memory_t::host) << std::endl;
@@ -34,22 +34,22 @@ int main(int argn, char** argv)
     }
 
     if (pr) {
-        std::cout << "allreduce<reinterpret_cast<double_complex>> " << std::endl;
+        std::cout << "allreduce<reinterpret_cast<std::complex<double>>> " << std::endl;
     }
     Communicator::world().allreduce(reinterpret_cast<double*>(tmp.at(memory_t::host)), 2 * static_cast<int>(tmp.size()));
 
     if (pr) {
-        std::cout << "reduce<reinterpret_cast<double_complex>> " << std::endl;
+        std::cout << "reduce<reinterpret_cast<std::complex<double>>> " << std::endl;
     }
     Communicator::world().reduce(reinterpret_cast<double*>(tmp.at(memory_t::host)), 2 * static_cast<int>(tmp.size()), 1);
 
     if (pr) {
-        std::cout << "allreduce<double_complex> " << std::endl;
+        std::cout << "allreduce<std::complex<double>> " << std::endl;
     }
     Communicator::world().allreduce(tmp.at(memory_t::host), static_cast<int>(tmp.size()));
 
     if (pr) {
-        std::cout << "reduce<double_complex> " << std::endl;
+        std::cout << "reduce<std::complex<double>> " << std::endl;
     }
     Communicator::world().reduce(tmp.at(memory_t::host), static_cast<int>(tmp.size()), 1);
 

apps/tests/test_sym.cpp

@@ -20,9 +20,9 @@ void test_sym(cmd_args const& args__)
     json_conf["parameters"]["num_bands"] = 1;
     //ctx.so_correction(true);
 
-    std::vector<geometry3d::vector3d<double>> coord;
-    coord.push_back(geometry3d::vector3d<double>({1.0/2, 1.0/2, 0}));
-    coord.push_back(geometry3d::vector3d<double>({1.0/3, 1.0/3, 0}));
+    std::vector<r3::vector<double>> coord;
+    coord.push_back(r3::vector<double>({1.0/2, 1.0/2, 0}));
+    coord.push_back(r3::vector<double>({1.0/3, 1.0/3, 0}));
     //ctx.unit_cell().add_atom("Cu", {0.113500, 0.613500, 0.886500});
     //ctx.unit_cell().add_atom("Cu", {0.613500, 0.886500, 0.113500});
     //ctx.unit_cell().add_atom("Cu", {0.886500, 0.113500, 0.613500});
@@ -33,8 +33,8 @@ void test_sym(cmd_args const& args__)
             coord, false, false);
     Simulation_context& ctx = *ctx_ptr;
 
-    vector3d<int> k_grid(4, 4, 1);
-    vector3d<int> k_shift(0, 0, 0);
+    r3::vector<int> k_grid(4, 4, 1);
+    r3::vector<int> k_shift(0, 0, 0);
 
     K_point_set kset_sym(ctx, k_grid, k_shift, true);
     K_point_set kset_nosym(ctx, k_grid, k_shift, false);
@@ -79,12 +79,12 @@ void test_sym(cmd_args const& args__)
 
             auto rotm = sht::rotation_matrix<double>(2, eang, pr);
 
-            auto vk1 = geometry3d::reduce_coordinates(dot(R, kset_sym.get<double>(ik)->vk())).first;
+            auto vk1 = r3::reduce_coordinates(dot(R, kset_sym.get<double>(ik)->vk())).first;
 
             std::cout << "isym: " << isym << " k: " << kset_sym.get<double>(ik)->vk() << " k1: " << vk1 << std::endl;
 
             /* compute <phi|G+k>w<G+k|phi> using k1 from the irreducible set */
-            sddk::mdarray<double_complex, 3> dm(5, 5, na);
+            sddk::mdarray<std::complex<double>, 3> dm(5, 5, na);
             dm.zero();
 
             int ik1 = kset_nosym.find_kpoint(vk1);
@@ -112,10 +112,10 @@ void test_sym(cmd_args const& args__)
                 int ja = sym[isym].spg_op.sym_atom[ia];
 
                 double phase = twopi * dot(kset_sym.get<double>(ik)->vk(), ctx.unit_cell().atom(ia).position());
-                auto dephase_k = std::exp(double_complex(0.0, phase));
+                auto dephase_k = std::exp(std::complex<double>(0.0, phase));
 
                 phase = twopi * dot(kset_sym.get<double>(ik)->vk(), ctx.unit_cell().atom(ja).position());
-                auto phase_k = std::exp(double_complex(0.0, phase));
+                auto phase_k = std::exp(std::complex<double>(0.0, phase));
 
                 std::cout << "ia : " << ia << " -> " << ja << std::endl;
 
@@ -126,7 +126,7 @@ void test_sym(cmd_args const& args__)
                 int jb = type_j.indexb_wfs().offset(2) + nawf.second[ja];
 
                 for (int ig = 0; ig < kset_sym.get<double>(ik)->num_gkvec(); ig++) {
-                    sddk::mdarray<double_complex, 1> v1(5);
+                    sddk::mdarray<std::complex<double>, 1> v1(5);
                     v1.zero();
                     for (int m = 0; m < 5; m++) {
                         for (int mp = 0; mp < 5; mp++) {
@@ -139,7 +139,7 @@ void test_sym(cmd_args const& args__)
                 }
             }
 
-            sddk::mdarray<double_complex, 3> dm1(5, 5, na);
+            sddk::mdarray<std::complex<double>, 3> dm1(5, 5, na);
             dm1.zero();
 
             for (int ia = 0; ia < na; ia++) {

apps/tests/test_wf_inner.cpp

@@ -66,7 +66,7 @@ void test_wf_inner(std::vector<int> mpi_grid_dims__, double cutoff__, int num_ba
         for (int i = 0; i < ovlp.num_rows_local(); i++) {
             auto irow = ovlp.irow(i);
             /* 2 is accumulated from two spins */
-            double_complex z = ovlp(i, j) - 2 * static_cast<double>(irow + 1) / (jcol + 1);
+            std::complex<double> z = ovlp(i, j) - 2 * static_cast<double>(irow + 1) / (jcol + 1);
             max_diff = std::max(max_diff, std::abs(z));
         }
     }