[WIP] generate wannier90 input

CMakeLists.txt

@@ -34,6 +34,7 @@ option(SIRIUS_USE_PROFILER            "measure execution of functions with timer
 option(SIRIUS_USE_MEMORY_POOL         "use memory pool" ON)
 option(SIRIUS_USE_PUGIXML             "enable reading of UPF v2 files with pugixml" OFF)
 option(SIRIUS_USE_POWER_COUNTER       "measure energy consumption with power counters" OFF)
+option(SIRIUS_USE_WANNIER90           "use Wannier90 library" OFF)
 option(BUILD_TESTING                  "build test executables" OFF) # override default setting in CTest module
 option(SIRIUS_USE_VCSQNM              "use variable cell stabilized quasi Newton method" OFF)
 
@@ -192,6 +193,10 @@ endif(SIRIUS_USE_VDWXC)
 
 find_package(costa CONFIG REQUIRED)
 
+if(SIRIUS_USE_WANNIER90)
+  find_package(Wannier90 REQUIRED)
+endif()
+
 if(SIRIUS_USE_CUDA)
   enable_language(CUDA)
   # note find cudatoolkit is called inside the include file. the
@@ -309,5 +314,6 @@ if(SIRIUS_BUILD_APPS)
   endif()
   add_subdirectory(apps/upf)
   add_subdirectory(apps/utils)
+  add_subdirectory(apps/sirius2wannier)
 endif(SIRIUS_BUILD_APPS)
 add_subdirectory(doc)

apps/mini_app/sirius.scf.cpp

@@ -82,19 +82,24 @@ preprocess_json_input(std::string fname__)
 std::unique_ptr<Simulation_context>
 create_sim_ctx(std::string fname__, cmd_args const& args__)
 {
-    auto json = preprocess_json_input(fname__);
+    std::string config_string;
+    if (isHDF5(fname__)) {
+        config_string = fname__;
+    } else {
+        auto json = preprocess_json_input(fname__);
+        config_string = json.dump();
+    }
 
-    auto ctx_ptr            = std::make_unique<Simulation_context>(json.dump(), mpi::Communicator::world());
-    Simulation_context& ctx = *ctx_ptr;
+    auto ctx = std::make_unique<Simulation_context>(config_string);
 
-    auto& inp = ctx.cfg().parameters();
+    auto& inp = ctx->cfg().parameters();
     if (inp.gamma_point() && !(inp.ngridk()[0] * inp.ngridk()[1] * inp.ngridk()[2] == 1)) {
         RTE_THROW("this is not a Gamma-point calculation")
     }
 
-    ctx.import(args__);
+    ctx->import(args__);
 
-    return ctx_ptr;
+    return ctx;
 }
 
 auto
@@ -110,6 +115,12 @@ ground_state(Simulation_context& ctx, int task_id, cmd_args const& args, int wri
                           << "+----------------------+" << std::endl;
                 break;
             }
+            case task_t::ground_state_restart: {
+                ctx.out() << "+--------------------------+" << std::endl
+                          << "| restart SCF ground state |" << std::endl
+                          << "+--------------------------+" << std::endl;
+                break;
+            }
             case task_t::ground_state_new_relax: {
                 ctx.out() << "+---------------------------------------------+" << std::endl
                           << "| new SCF ground state with atomic relaxation |" << std::endl
@@ -145,11 +156,19 @@ ground_state(Simulation_context& ctx, int task_id, cmd_args const& args, int wri
     auto& density   = dft.density();
 
     if (task_id == task_t::ground_state_restart) {
-        if (!file_exists(storage_file_name)) {
+        auto fname = args.value<fs::path>("input", storage_file_name);
+        if (!isHDF5(fname)) {
+            fname = storage_file_name;
+        }
+        if (!file_exists(fname)) {
             RTE_THROW("storage file is not found");
         }
-        density.load(storage_file_name);
-        potential.load(storage_file_name);
+        density.load(fname);
+        density.generate_paw_density();
+        //potential.load(fname);
+        potential.generate(density, ctx.use_symmetry(), true);
+        Hamiltonian0<double> H0(potential, true);
+        initialize_subspace(kset, H0);
     } else {
         dft.initial_state();
     }
@@ -222,6 +241,10 @@ ground_state(Simulation_context& ctx, int task_id, cmd_args const& args, int wri
     }
 
     if (write_state && write_output) {
+        std::cout << "Saving wf.." << std::endl;
+        kset.save("sirius.h5");
+        std::cout << " wf saved.." << std::endl;
+
         json dict;
         json_output_common(dict);
 
@@ -239,6 +262,11 @@ ground_state(Simulation_context& ctx, int task_id, cmd_args const& args, int wri
             write_json_to_file(dict, output_file);
         }
 
+///#if defined(SIRIUS_WANNIER90)
+    if (ctx.cfg().parameters().wannier()) {
+        kset.generate_w90_coeffs();
+    }
+////#endif
         // if (args.exist("aiida_output")) {
         //     json dict;
         //     json_output_common(dict);
@@ -264,6 +292,8 @@ ground_state(Simulation_context& ctx, int task_id, cmd_args const& args, int wri
         dft.check_scf_density();
     }
 
+
+
     auto repeat_update = args.value<int>("repeat_update", 0);
     if (repeat_update) {
         auto lv = ctx.unit_cell().lattice_vectors();
@@ -579,6 +609,8 @@ main(int argn, char** argv)
                    {"iterative_solver.orthogonalize=", ""},
                    {"iterative_solver.early_restart=",
                     "{double} value between 0 and 1 to control the early restart ratio in Davidson"},
+                   {"iterative_solver.energy_tolerance=", "{double} starting tolerance of iterative solver"},
+                   {"iterative_solver.num_steps=", "{int} number of steps in iterative solver"},
                    {"mixer.type=", "{string} mixer name (anderson, anderson_stable, broyden2, linear)"},
                    {"mixer.beta=", "{double} mixing parameter"},
                    {"volume_scale0=", "{double} starting volume scale for EOS calculation"},

apps/sirius2wannier/CMakeLists.txt

@@ -0,0 +1,3 @@
+add_executable(sirius2wannier sirius2wannier.cpp)
+target_link_libraries(sirius2wannier PRIVATE sirius_cxx)
+install(TARGETS sirius2wannier RUNTIME DESTINATION "${CMAKE_INSTALL_PREFIX}/bin")

apps/sirius2wannier/sirius2wannier.cpp

@@ -0,0 +1,137 @@
+#include "apps.hpp"
+#include "k_point/generate_w90_coeffs.hpp"
+#include "nlcglib/apply_hamiltonian.hpp"
+
+std::vector<std::array<double, 3>>  
+load_coordinates( const std::string& fname__ )
+{
+    std::vector<std::array<double, 3>> kp;
+
+    //read k coordinates from hdf5
+    HDF5_tree fin(fname__, hdf5_access_t::read_only);
+    std::cout << "read num_kpoints" << std::endl;
+    int num_kpoints;
+    fin["K_point_set"].read("num_kpoints", &num_kpoints, 1);
+    std::cout << "num_kpoints: " << num_kpoints << std::endl;
+    kp.resize(num_kpoints);
+    for( int ik = 0; ik < num_kpoints; ik++ ) {
+        fin["K_point_set"][ik].read("vk", &kp[ik][0], 3);
+        std::cout << "ik = " << ik << " kp = { " << kp[ik][0] << " " << kp[ik][1] << " " << kp[ik][2] << std::endl;
+    }
+    return kp;
+}
+
+
+int
+main(int argn, char** argv)
+{
+    cmd_args args(argn, argv,
+                  {{"input=", "{string} input file name"}});
+
+    sirius::initialize(1);
+
+
+    /* get the input file name */
+    auto fpath = args.value<fs::path>("input", "state.h5");
+
+    if (fs::is_directory(fpath)) {
+        fpath /= "sirius.h5";
+    }
+
+    if (!fs::exists(fpath)) {
+        if (mpi::Communicator::world().rank() == 0) {
+            std::cout << "input file does not exist" << std::endl;
+        }
+        exit(1);
+    }
+    auto fname = fpath.string();
+
+    /* create simulation context */
+    auto ctx = create_sim_ctx(fname, args);
+    ctx->initialize();
+
+    /* read the wf */
+    auto kp = load_coordinates(fname);
+    K_point_set kset(*ctx, kp);
+    std::cout << "kset initialized.\n";
+    kset.load(fname);
+
+    /* initialize the ground state */
+    DFT_ground_state dft(kset);
+    auto& potential = dft.potential();
+    auto& density   = dft.density();
+    density.load(fname);
+    density.generate_paw_density();
+    //potential.load(fname);
+    potential.generate(density, ctx->use_symmetry(), true);
+    Hamiltonian0<double> H0(potential, true);
+
+    /* checksum over wavefunctions */
+    //for (auto it : kset.spl_num_kpoints()) {
+    //    int ik = it.i;
+    //    auto Hk = H0(*kset.get<double>(ik));
+    //    for (auto is=0; is< ctx->num_spins(); is++) {
+    //      std::cout << "ik: " << ik << " ispn : "<< is << " " ; 
+    //      std::cout << kset.get<double>(ik)->spinor_wave_functions().checksum(memory_t::host, wf::spin_index(is), wf::band_range(0, ctx->num_bands())) << std::endl;
+    //    }
+    //}
+    /* check if the wfs diagonalize the hamiltonian and if the eigenvalues are correct */
+
+    la::dmatrix<std::complex<double>> psiHpsi(ctx->num_bands(), ctx->num_bands());
+
+    std::cout << "num_spins " << kset.ctx().num_spins() << std::endl;
+
+    for (auto it : kset.spl_num_kpoints()) {
+        int ik = it.i;
+        std::cout << "ik = " << ik << std::endl;
+        auto kp = kset.get<double>(ik);
+        auto Hk = H0(*kp);
+
+        bool nc_mag     = (kset.ctx().num_mag_dims() == 3);
+        int num_spinors = (kset.ctx().num_mag_dims() == 1) ? 2 : 1;
+        int num_sc      = nc_mag ? 2 : 1;
+
+        auto hpsi = std::make_unique<wf::Wave_functions<double>>(kp->gkvec_sptr(), 
+                                                                 wf::num_mag_dims(kset.ctx().num_mag_dims()),
+                                                                 wf::num_bands(ctx->num_bands()), 
+                                                                 memory_t::host);
+        auto spsi = std::make_unique<wf::Wave_functions<double>>(kp->gkvec_sptr(), 
+                                                                 wf::num_mag_dims(kset.ctx().num_mag_dims()),
+                                                                 wf::num_bands(ctx->num_bands()), 
+                                                                 memory_t::host);
+
+        for (int ispin_step = 0; ispin_step < num_spinors; ispin_step++) {
+            auto sr = nc_mag ? wf::spin_range(0, 2) : wf::spin_range(ispin_step);
+            std::cout << "ik= " << ik << " ispin_step = " << ispin_step; 
+            /* get H|psi> */
+            Hk.apply_h_s<std::complex<double>>(sr, wf::band_range(0, ctx->num_bands()), kp->spinor_wave_functions(), hpsi.get(), spsi.get());
+
+            /* get <psi|H|psi> */
+            wf::inner(kset.ctx().spla_context(), memory_t::host, sr, kp->spinor_wave_functions(), wf::band_range(0, ctx->num_bands()),
+                      *hpsi, wf::band_range(0, ctx->num_bands()), psiHpsi, 0, 0);
+
+            /* check elements that are large compared to the threshold */
+            std::vector<std::pair<int,int>> indices;
+            std::vector<double> comp;
+            for( int ibnd = 0; ibnd < ctx->num_bands(); ++ibnd ) {
+                for( int jbnd = 0; jbnd < ctx->num_bands(); ++jbnd ) {
+                    auto comparison = (ibnd == jbnd) ? kset.get<double>(ik)->band_energies(ispin_step)[ibnd]  : 0.;
+
+                    if( std::abs ( psiHpsi(jbnd,ibnd) - comparison ) > 1.e-07  ) {
+                        indices.push_back(std::pair<int,int>(jbnd, ibnd));
+                        comp.push_back(comparison);
+                    }
+                }
+            }
+            for( auto i = 0; i < indices.size(); ++i ) {
+                auto& i_ = indices[i];
+                std::cout << "      element = (";
+                std::cout << i_.first << ", " << i_.second << ") = " << psiHpsi(i_.first, i_.second) << " comparison : " << comp[i] << std::endl;
+            }
+        }//ispin_step
+    }//kpoint
+
+exit(0);
+
+    kset.generate_w90_coeffs();
+}

apps/utils/unit_cell_tools.cpp

@@ -34,7 +34,7 @@ create_supercell(cmd_args const& args__)
         std::cout << std::endl;
     }
 
-    Simulation_context ctx("sirius.json", mpi::Communicator::self());
+    Simulation_context ctx(std::string("sirius.json"), mpi::Communicator::self());
 
     auto scell_lattice_vectors = dot(ctx.unit_cell().lattice_vectors(), r3::matrix<double>(scell));
 
@@ -113,7 +113,7 @@ create_supercell(cmd_args const& args__)
 void
 find_primitive()
 {
-    Simulation_context ctx("sirius.json", mpi::Communicator::self());
+    Simulation_context ctx(std::string("sirius.json"), mpi::Communicator::self());
 
     double lattice[3][3];
     for (int i : {0, 1, 2}) {

cmake/modules/FindWannier90.cmake

@@ -0,0 +1,25 @@
+include(FindPackageHandleStandardArgs)
+find_package(PkgConfig REQUIRED)
+
+find_library(SIRIUS_WANNIER90_LIBRARIES NAMES wannier wannier90
+  PATH_SUFFIXES lib
+  HINTS
+  ENV EBROOTWANNIER90
+  ENV WANNIER90_ROOT
+)
+
+find_path(SIRIUS_WANNIER90_INCLUDE_DIR NAMES w90_wannierise.mod
+  PATH_SUFFIXES modules
+  HINTS
+  ENV EBROOTWANNIER90
+  ENV WANNIER90_ROOT
+  )
+
+find_package_handle_standard_args(Wannier90 DEFAULT_MSG SIRIUS_WANNIER90_LIBRARIES SIRIUS_WANNIER90_INCLUDE_DIR)
+
+if(Wannier90_FOUND AND NOT TARGET sirius::wannier90)
+  add_library(sirius::wannier90 INTERFACE IMPORTED)
+  set_target_properties(sirius::wannier90 PROPERTIES
+    INTERFACE_INCLUDE_DIRECTORIES "${SIRIUS_WANNIER90_INCLUDE_DIR}"
+    INTERFACE_LINK_LIBRARIES "${SIRIUS_WANNIER90_LIBRARIES}")
+endif()

examples/pp-pw/Si7Ge/sirius.json

@@ -16,7 +16,7 @@
         "gk_cutoff" : 5.0,
         "pw_cutoff" : 20.00,
 
-        "energy_tol" : 1e-8,
+        "energy_tol" : 1e-7,
         "density_tol" : 1e-7,
 
         "num_dft_iter" : 100,
@@ -33,7 +33,8 @@
     "mixer" : {
       "beta" : 0.8,
       "type" : "anderson",
-      "max_history" : 8
+      "max_history" : 8,
+      "use_hartree" : true
     },
 
     "unit_cell": {

src/CMakeLists.txt

@@ -59,6 +59,7 @@ set(_SOURCES
   "k_point/generate_spinor_wave_functions.cpp"
   "k_point/k_point.cpp"
   "k_point/k_point_set.cpp"
+  "k_point/generate_w90_coeffs.cpp"
   "radial/radial_integrals.cpp"
   "mixer/mixer_functions.cpp"
   "nlcglib/adaptor.cpp"
@@ -113,6 +114,7 @@ target_link_libraries(sirius_cxx PUBLIC ${GSL_LIBRARY}
   $<TARGET_NAME_IF_EXISTS:sirius::libvdwxc>
   $<TARGET_NAME_IF_EXISTS:nlcglib::nlcglib>
   $<TARGET_NAME_IF_EXISTS:kokkos::kokkos>
+  $<TARGET_NAME_IF_EXISTS:sirius::wannier90>
   SpFFT::spfft
   SPLA::spla
   "${SIRIUS_LINALG_LIB}"
@@ -127,7 +129,9 @@ target_link_libraries(sirius_cxx PUBLIC ${GSL_LIBRARY}
 )
 
 target_include_directories(sirius_cxx PUBLIC $<BUILD_INTERFACE:${PROJECT_SOURCE_DIR}/src>
-                                             $<BUILD_INTERFACE:${PROJECT_BINARY_DIR}/src>)
+                                         $<BUILD_INTERFACE:${PROJECT_BINARY_DIR}/src>
+                                         $<BUILD_INTERFACE:${PROJECT_BINARY_DIR}/src/mod_files>)
+
 target_compile_definitions(sirius_cxx PUBLIC
   $<$<BOOL:${SIRIUS_USE_PROFILER}>:SIRIUS_PROFILE>
   $<$<BOOL:${SIRIUS_USE_SCALAPACK}>:SIRIUS_SCALAPACK>
@@ -146,6 +150,7 @@ target_compile_definitions(sirius_cxx PUBLIC
   $<$<BOOL:${SIRIUS_USE_VCSQNM}>:SIRIUS_VCSQNM>
   $<$<BOOL:${SIRIUS_HAVE_LIBVDW_WITH_MPI}>:SIRIUS_HAVE_VDWXC_MPI>
   $<$<AND:$<BOOL:${SIRIUS_USE_MAGMA}>,$<BOOL:${SIRIUS_USE_ROCM}>>:HAVE_HIP> # Required for magma headers
+  $<$<BOOL:${SIRIUS_USE_WANNIER90}>:SIRIUS_WANNIER90>
 )
 
 if(SIRIUS_CREATE_FORTRAN_BINDINGS)