Merge pull request #387 from toxa81/develop

add extra charge

VERSION

@@ -1 +1 @@
-6.2.1
+6.3.0

apps/tests/test_davidson.cpp

@@ -38,7 +38,7 @@ void init_wf(K_point* kp__, Wave_functions& phi__, int num_bands__, int num_mag_
     }
 }
 
-void test_davidson()
+void test_davidson(device_t pu__, double pw_cutoff__, double gk_cutoff__, int N__)
 {
     utils::timer t1("test_davidson|setup");
 
@@ -98,31 +98,32 @@ void test_davidson()
     }
     atype.ps_total_charge_density(arho);
 
-    int N{2};
     /* lattice constant */
     double a{5};
     /* set lattice vectors */
-    ctx.unit_cell().set_lattice_vectors({{a * N, 0, 0}, 
-                                         {0, a * N, 0}, 
-                                         {0, 0, a * N}});
+    ctx.unit_cell().set_lattice_vectors({{a * N__, 0, 0}, 
+                                         {0, a * N__, 0}, 
+                                         {0, 0, a * N__}});
     /* add atoms */
-    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++) {
+    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++) {
                 ctx.unit_cell().add_atom("Cu", {i * p, j * p, k * p});
             }
         }
     }
 
     /* initialize the context */
     ctx.set_verbosity(1);
-    ctx.pw_cutoff(40);
-    ctx.gk_cutoff(10);
+    ctx.pw_cutoff(pw_cutoff__);
+    ctx.gk_cutoff(gk_cutoff__);
+    ctx.set_processing_unit(pu__);
     t1.stop();
-    ctx.set_verbosity(5);
+
+    ctx.set_verbosity(1);
     ctx.iterative_solver_tolerance(1e-12);
-    ctx.set_iterative_solver_type("exact");
+    //ctx.set_iterative_solver_type("exact");
     ctx.initialize();
 
     Density rho(ctx);
@@ -134,23 +135,34 @@ void test_davidson()
     pot.zero();
 
 
-    double vk[] = {0.1, 0.1, 0.1};
-    K_point kp(ctx, vk, 1.0);
-    kp.initialize();
-    std::cout << "num_gkvec=" << kp.num_gkvec() << "\n";
-    for (int i = 0; i < ctx.num_bands(); i++) {
-        kp.band_occupancy(i, 0, 2);
-    }
-    init_wf(&kp, kp.spinor_wave_functions(), ctx.num_bands(), 0);
+    for (int r = 0; r < 2; r++) {
+        double vk[] = {0.1, 0.1, 0.1};
+        K_point kp(ctx, vk, 1.0);
+        kp.initialize();
+        std::cout << "num_gkvec=" << kp.num_gkvec() << "\n";
+        for (int i = 0; i < ctx.num_bands(); i++) {
+            kp.band_occupancy(i, 0, 2);
+        }
+        init_wf(&kp, kp.spinor_wave_functions(), ctx.num_bands(), 0);
 
 
 
-    Hamiltonian H(ctx, pot);
-    H.prepare();
-    //Band(ctx).solve_pseudo_potential<double_complex>(kp, H);
-    //for (int i = 0; i < ctx.num_bands(); i++) {
-    //    std::cout << "energy[" << i << "]=" << kp.band_energy(i, 0) << "\n";
-    //}
+        Hamiltonian H(ctx, pot);
+        H.prepare();
+        Band(ctx).solve_pseudo_potential<double_complex>(kp, H);
+        //for (int i = 0; i < ctx.num_bands(); i++) {
+        //    std::cout << "energy[" << i << "]=" << kp.band_energy(i, 0) << "\n";
+        //}
+        std::vector<double> ekin(kp.num_gkvec());
+        for (int i = 0; i < kp.num_gkvec(); i++) {
+            ekin[i] = 0.5 * kp.gkvec().gkvec_cart<index_domain_t::global>(i).length2();
+        }
+        std::sort(ekin.begin(), ekin.end());
+
+        for (int i = 0; i < ctx.num_bands(); i++) {
+            printf("%20.16f %20.16f %20.16e\n", ekin[i], kp.band_energy(i, 0), std::abs(ekin[i] - kp.band_energy(i, 0)));
+        }
+    }
 
     //ctx.set_iterative_solver_type("davidson");
     //Band(ctx).solve_pseudo_potential<double_complex>(kp, H);
@@ -161,22 +173,22 @@ void test_davidson()
     //    K_point kp(ctx, vk, 1.0);
     //    kp.initialize();
 
-    Hamiltonian0 h0(ctx, pot);
-    auto hk = h0(kp);
+    //Hamiltonian0 h0(ctx, pot);
+    //auto hk = h0(kp);
 
-    //    init_wf(&kp, kp.spinor_wave_functions(), ctx.num_bands(), 0);
+    ////    init_wf(&kp, kp.spinor_wave_functions(), ctx.num_bands(), 0);
 
-    auto eval = davidson<double_complex>(hk, kp.spinor_wave_functions(), 0, 4, 40, 1e-12, 0, 1e-7, [](int i, int ispn){return 1.0;}, false);
+    //auto eval = davidson<double_complex>(hk, kp.spinor_wave_functions(), 0, 4, 40, 1e-12, 0, 1e-7, [](int i, int ispn){return 1.0;}, false);
 
-    std::vector<double> ekin(kp.num_gkvec());
-    for (int i = 0; i < kp.num_gkvec(); i++) {
-        ekin[i] = 0.5 * kp.gkvec().gkvec_cart<index_domain_t::global>(i).length2();
-    }
-    std::sort(ekin.begin(), ekin.end());
+    //std::vector<double> ekin(kp.num_gkvec());
+    //for (int i = 0; i < kp.num_gkvec(); i++) {
+    //    ekin[i] = 0.5 * kp.gkvec().gkvec_cart<index_domain_t::global>(i).length2();
+    //}
+    //std::sort(ekin.begin(), ekin.end());
 
-    for (int i = 0; i < ctx.num_bands(); i++) {
-        printf("%20.16f %20.16f %20.16e\n", ekin[i], eval[i], std::abs(ekin[i] - eval[i]));
-    }
+    //for (int i = 0; i < ctx.num_bands(); i++) {
+    //    printf("%20.16f %20.16f %20.16e\n", ekin[i], eval[i], std::abs(ekin[i] - eval[i]));
+    //}
 
     //}
 
@@ -216,17 +228,27 @@ void test_davidson()
 
 int main(int argn, char** argv)
 {
-    cmd_args args;
+    cmd_args args(argn, argv, {{"device=", "(string) CPU or GPU"},
+                               {"pw_cutoff=", "(double) plane-wave cutoff for density and potential"},
+                               {"gk_cutoff=", "(double) plane-wave cutoff for wave-functions"},
+                               {"N=", "(int) cell multiplicity"},
+                               {"mpi_grid=", "(int[2]) dimensions of the MPI grid for band diagonalization"}
+                              });
 
-    args.parse_args(argn, argv);
     if (args.exist("help")) {
         printf("Usage: %s [options]\n", argv[0]);
         args.print_help();
         return 0;
     }
 
+    auto pu        = get_device_t(args.value<std::string>("device", "CPU"));
+    auto pw_cutoff = args.value<double>("pw_cutoff", 30);
+    auto gk_cutoff = args.value<double>("gk_cutoff", 10);
+    auto N         = args.value<int>("N", 1);
+    auto mpi_grid  = args.value<std::vector<int>>("mpi_grid", {1, 1});
+
     sirius::initialize(1);
-    test_davidson();
+    test_davidson(pu, pw_cutoff, gk_cutoff, N);
     int rank = Communicator::world().rank();
     sirius::finalize();
     if (!rank) {

doc/doxygen.cfg

@@ -38,7 +38,7 @@ PROJECT_NAME           = "SIRIUS"
 # could be handy for archiving the generated documentation or if some version
 # control system is used.
 
-PROJECT_NUMBER         = "6.2.1"
+PROJECT_NUMBER         = "6.3.0"
 
 # Using the PROJECT_BRIEF tag one can provide an optional one line description
 # for a project that appears at the top of each page and should give viewer a

python_module/CMakeLists.txt

@@ -2,6 +2,8 @@ if(CREATE_PYTHON_MODULE)
   set(CMAKE_CXX_STANDARD 14)
   find_package(mpi4py REQUIRED)
 
+  add_definitions("-DPYBIND11_CPP14")
+
   set(pb11_src_dir "${PROJECT_SOURCE_DIR}/python_module/pybind11")
   check_git_submodule(pybind11 "${pb11_src_dir}")
   if(NOT pybind11_avail)

python_module/magnetization.hpp

@@ -12,11 +12,9 @@ std::tuple<std::list<std::vector<double>>, std::list<double>> magnetization(Dens
     std::list<double> lm;
 
     for (int j=0; j < ctx.num_mag_dims(); ++j) {
-        std::vector<double> vec(3);
-        double m;
-        density.magnetization(j).integrate(vec, m);
-        lvec.push_back(vec);
-        lm.push_back(m);
+        auto result = density.magnetization(j).integrate();
+        lvec.push_back(std::get<2>(result));
+        lm.push_back(std::get<1>(result));
     }
 
     return std::make_tuple(lvec, lm);

src/Band/diag_pseudo_potential.hpp

@@ -456,14 +456,15 @@ inline int Band::diag_pseudo_potential_davidson(K_point*       kp__,
         /* current subspace size */
         int N = num_bands;
 
-        utils::timer t1("sirius::Band::diag_pseudo_potential_davidson|evp");
+        auto p1 = std::unique_ptr<utils::profiler>(new utils::profiler(__function_name__, __FILE__, __LINE__,
+            "sirius::Band::diag_pseudo_potential_davidson|evp"));
         /* solve generalized eigen-value problem with the size N and get lowest num_bands eigen-vectors */
         if (std_solver.solve(N, num_bands, hmlt, eval.data(), evec)) {
             std::stringstream s;
             s << "error in diagonalziation";
             TERMINATE(s);
         }
-        t1.stop();
+        p1 = nullptr;
 
         evp_work_count() += 1;
 
@@ -585,7 +586,8 @@ inline int Band::diag_pseudo_potential_davidson(K_point*       kp__,
 
             eval_old = eval;
 
-            utils::timer t1("sirius::Band::diag_pseudo_potential_davidson|evp");
+            auto p1 = std::unique_ptr<utils::profiler>(new utils::profiler(__function_name__, __FILE__, __LINE__,
+                "sirius::Band::diag_pseudo_potential_davidson|evp"));
             if (itso.orthogonalize_) {
                 /* solve standard eigen-value problem with the size N */
                 if (std_solver.solve(N, num_bands, hmlt, eval.data(), evec)) {
@@ -601,7 +603,7 @@ inline int Band::diag_pseudo_potential_davidson(K_point*       kp__,
                     TERMINATE(s);
                 }
             }
-            t1.stop();
+            p1 = nullptr;
 
             evp_work_count() += std::pow(static_cast<double>(N) / num_bands, 3);
 

src/CMakeLists.txt

@@ -17,10 +17,10 @@ endif()
 set(_SOURSES "dft_ground_state.cpp"
   "Hubbard/apply_hubbard_potential.cpp"
   "Hubbard/Hubbard.cpp"
-"Hubbard/hubbard_generate_atomic_orbitals.cpp"
-"Hubbard/hubbard_occupancies_derivatives.cpp"
-"Hubbard/hubbard_occupancy.cpp"
-"Hubbard/hubbard_potential_energy.cpp")
+  "Hubbard/hubbard_generate_atomic_orbitals.cpp"
+  "Hubbard/hubbard_occupancies_derivatives.cpp"
+  "Hubbard/hubbard_occupancy.cpp"
+  "Hubbard/hubbard_potential_energy.cpp")
 
 # create library with .cpp, .cu and .f90 sources
 add_library(sirius STATIC "${_SOURSES};${_CUSOURCES};${_FSOURCES}")
@@ -41,13 +41,6 @@ target_link_libraries(sirius PUBLIC OpenMP::OpenMP_CXX)
 
 if(CREATE_FORTRAN_BINDINGS)
   set_target_properties(sirius PROPERTIES Fortran_MODULE_DIRECTORY mod_files)
-endif()
-
-#if(USE_CUDA)
-#  add_library(sirius STATIC "${_CUSOURCES}")
-#endif()
-
-if(CREATE_FORTRAN_BINDINGS)
   INSTALL ( CODE
     "EXECUTE_PROCESS (COMMAND \"${CMAKE_COMMAND}\" -E copy_directory \"${PROJECT_BINARY_DIR}/src/mod_files\" \"${CMAKE_INSTALL_PREFIX}/include/sirius/\")"
     )

src/Density/density.hpp

@@ -139,7 +139,7 @@ class Density : public Field4D
 
         int ia{-1};
 
-        /// ae and ps local unified densities+magnetization
+        /* AE and PS local unified densities + magnetization */
         std::vector<Spheric_function<function_domain_t::spectral, double>> ae_density_;
         std::vector<Spheric_function<function_domain_t::spectral, double>> ps_density_;
     };
@@ -385,9 +385,7 @@ class Density : public Field4D
 
     inline void normalize()
     {
-        std::vector<double> nel_mt;
-        double              nel_it;
-        double nel   = rho().integrate(nel_mt, nel_it);
+        double nel = std::get<0>(rho().integrate());
         double scale = unit_cell_.num_electrons() / nel;
 
         /* renormalize interstitial part */
@@ -410,9 +408,7 @@ class Density : public Field4D
     {
         double nel{0};
         if (ctx_.full_potential()) {
-            std::vector<double> nel_mt;
-            double              nel_it;
-            nel = rho().integrate(nel_mt, nel_it);
+            nel = std::get<0>(rho().integrate());
         } else {
             nel = rho().f_0().real() * unit_cell_.omega();
         }

src/Density/generate_valence.hpp

@@ -45,12 +45,13 @@ inline void Density::generate_valence(K_point_set const& ks__)
         TERMINATE(s);
     }
 
-    if (std::abs(occ_val - unit_cell_.num_valence_electrons()) > 1e-8 && ctx_.comm().rank() == 0) {
+    if (std::abs(occ_val - unit_cell_.num_valence_electrons() + ctx_.parameters_input().extra_charge_) > 1e-8 &&
+        ctx_.comm().rank() == 0) {
         std::stringstream s;
         s << "wrong band occupancies" << std::endl
           << "  computed : " << occ_val << std::endl
-          << "  required : " << unit_cell_.num_valence_electrons() << std::endl
-          << "  difference : " << std::abs(occ_val - unit_cell_.num_valence_electrons());
+          << "  required : " << unit_cell_.num_valence_electrons() - ctx_.parameters_input().extra_charge_ << std::endl
+          << "  difference : " << std::abs(occ_val - unit_cell_.num_valence_electrons() + ctx_.parameters_input().extra_charge_);
         WARNING(s);
     }
 
@@ -131,6 +132,12 @@ inline void Density::generate_valence(K_point_set const& ks__)
     if (!ctx_.full_potential()) {
         augment();
 
+        /* remove extra chanrge */
+        if (ctx_.gvec().comm().rank() == 0) {
+            rho().f_pw_local(0) += ctx_.parameters_input().extra_charge_ / ctx_.unit_cell().omega();
+
+        }
+
         if (ctx_.control().print_hash_ && ctx_.comm().rank() == 0) {
             auto h = mdarray<double_complex, 1>(&rho().f_pw_local(0), ctx_.gvec().count()).hash();
             utils::print_hash("rho", h);

src/Hamiltonian/non_local_operator.hpp

@@ -703,9 +703,10 @@ class Q_operator : public Non_local_operator
 //};
 
 /// Apply non-local part of the Hamiltonian and S operators.
-/** This operations must be combined becuase of the expensive inner product between wave-functions and beta
- *  projectors, which is computed only once. 
- *  \param [in]  spins   Range of the spin index 
+/** This operations must be combined because of the expensive inner product between wave-functions and beta
+ *  projectors, which is computed only once.
+ *
+ *  \param [in]  spins   Range of the spin index.
  *  \param [in]  N       Starting index of the wave-functions.
  *  \param [in]  n       Number of wave-functions to which D and Q are applied.
  *  \param [in]  beta    Beta-projectors.

src/K_point/generate_atomic_wave_functions.hpp

@@ -40,7 +40,7 @@ inline void K_point::generate_atomic_wave_functions(const basis_functions_index
     lmax = std::max(lmax, unit_cell_.lmax());
 
     auto& atom_type   = unit_cell_.atom(atom).type();
-#pragma omp parallel for schedule(static)
+    #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);
@@ -131,7 +131,7 @@ inline void K_point::compute_gradient_wave_functions(Wave_functions& phi,
         qalpha[igk_loc] = double_complex(0.0, -G[direction]);
     }
 
-#pragma omp parallel for schedule(static)
+    #pragma omp parallel for schedule(static)
     for (int nphi = 0; nphi < num_wf; nphi++) {
         for (int ispn = 0; ispn < phi.num_sc(); ispn++) {
             for (int igk_loc = 0; igk_loc < this->num_gkvec_loc(); igk_loc++) {