[feature] Implementation of the interface to geomtry optimization. (#778

)

Add interface to vc-sqnm code to perform lattice relaxation. Input schema is extended with `vcsqnm` section and correcponding input parameters are added.

CMakeLists.txt

@@ -26,9 +26,9 @@ set(BUILD_TESTING           OFF CACHE BOOL "build test executables") # override
 set(CUDA_ARCH               OFF CACHE STRING "the target GPU architectures; 60 (Pascal), 70 (Volta), etc. Multiple archs are supported too.")
 set(GPU_MEMORY_ALIGMENT "512" CACHE STRING "memory alignment of the GPU")
 set(USE_FP32 "AUTO" CACHE STRING "Enable single precision support.")
-set_property(CACHE USE_FP32 PROPERTY STRINGS
-  "AUTO" "ON" "OFF"
-  )
+set(USE_VCSQNM              OFF CACHE BOOL "use variable cell stabilized quasi Newton method")
+
+set_property(CACHE USE_FP32 PROPERTY STRINGS "AUTO" "ON" "OFF")
 
 # set language and standard
 

apps/dft_loop/CMakeLists.txt

@@ -1,5 +1,11 @@
+if(USE_VCSQNM)
+  find_package (Eigen3 3.3 REQUIRED NO_MODULE)
+endif()
 add_executable(sirius.scf sirius.scf.cpp)
 target_link_libraries(sirius.scf PRIVATE sirius sirius::filesystem)
+if(USE_VCSQNM)
+  target_link_libraries (sirius.scf PRIVATE Eigen3::Eigen sirius)
+endif()
 install(TARGETS sirius.scf RUNTIME DESTINATION "${CMAKE_INSTALL_PREFIX}/bin")
 set_property(TARGET sirius.scf PROPERTY POSITION_INDEPENDENT_CODE OFF)
 

apps/dft_loop/sirius.scf.cpp

@@ -1,29 +1,31 @@
-#include "utils/profiler.hpp"
 #include <sirius.hpp>
-#include "utils/filesystem.hpp"
-#include <utils/json.hpp>
 #include <cfenv>
-#include <fenv.h>
+#include "utils/profiler.hpp"
+#include "utils/filesystem.hpp"
+#include "utils/json.hpp"
+#include "dft/lattice_relaxation.hpp"
 
 using namespace sirius;
 using json = nlohmann::json;
 using namespace sddk;
 
 const std::string aiida_output_file = "output_aiida.json";
 
-enum class task_t : int
+struct task_t
 {
-    ground_state_new     = 0,
-    ground_state_restart = 1,
-    k_point_path         = 2,
-    eos                  = 3,
-    read_config          = 4
+    static const int ground_state_new         = 0;
+    static const int ground_state_restart     = 1;
+    static const int k_point_path             = 2;
+    static const int eos                      = 3;
+    static const int read_config              = 4;
+    static const int ground_state_new_relax   = 5;
+    static const int ground_state_new_vcrelax = 6;
 };
 
 void json_output_common(json& dict__)
 {
     dict__["git_hash"] = sirius::git_hash();
-    dict__["comm_world_size"] = Communicator::world().size();
+    dict__["comm_world_size"] = sddk::Communicator::world().size();
     dict__["threads_per_rank"] = omp_get_max_threads();
 }
 
@@ -84,10 +86,36 @@ create_sim_ctx(std::string fname__, cmd_args const& args__)
 }
 
 double
-ground_state(Simulation_context& ctx, task_t task, cmd_args const& args, int write_output)
+ground_state(Simulation_context& ctx, int task_id, cmd_args const& args, int write_output)
 {
     print_memory_usage(ctx.out(), FILE_LINE);
 
+    if (ctx.comm().rank() == 0) {
+        switch (task_id) {
+            case task_t::ground_state_new: {
+                ctx.out() << "+----------------------+" << std::endl
+                          << "| new 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
+                          << "+---------------------------------------------+" << std::endl;
+                break;
+            }
+            case task_t::ground_state_new_vcrelax: {
+                ctx.out() << "+---------------------------------------------------------+" << std::endl
+                          << "| new SCF ground state with atomic and lattice relaxation |" << std::endl
+                          << "+---------------------------------------------------------+" << std::endl;
+                break;
+            }
+            default: {
+                break;
+            }
+        }
+    }
+
     auto& inp = ctx.cfg().parameters();
 
     std::string ref_file = args.value<std::string>("test_against", "");
@@ -96,55 +124,95 @@ ground_state(Simulation_context& ctx, task_t task, cmd_args const& args, int wri
 
     bool const reduce_kp = ctx.use_symmetry() && ctx.cfg().parameters().use_ibz();
     K_point_set kset(ctx, ctx.cfg().parameters().ngridk(), ctx.cfg().parameters().shiftk(), reduce_kp);
+
     DFT_ground_state dft(kset);
 
     print_memory_usage(ctx.out(), FILE_LINE);
 
     auto& potential = dft.potential();
     auto& density = dft.density();
 
-    if (task == task_t::ground_state_restart) {
+    if (task_id == task_t::ground_state_restart) {
         if (!utils::file_exists(storage_file_name)) {
-            TERMINATE("storage file is not found");
+            RTE_THROW("storage file is not found");
         }
         density.load();
         potential.load();
     } else {
         dft.initial_state();
     }
 
-    /* launch the calculation */
-    auto result = dft.find(inp.density_tol(), inp.energy_tol(), ctx.cfg().iterative_solver().energy_tolerance(),
-            inp.num_dft_iter(), write_state);
-    /* compute forces and stress */
+    bool compute_stress{false};
+    bool compute_forces{false};
     if (ctx.cfg().control().print_stress() && !ctx.full_potential()) {
-        Stress& s       = dft.stress();
-        auto stress_tot = s.calc_stress_total();
-        s.print_info();
-        result["stress"] = std::vector<std::vector<double>>(3, std::vector<double>(3));
-        for (int i = 0; i < 3; i++) {
-            for (int j = 0; j < 3; j++) {
-                result["stress"][i][j] = stress_tot(j, i);
-            }
-        }
+        compute_stress = true;
     }
     if (ctx.cfg().control().print_forces()) {
-        Force& f         = dft.forces();
-        auto& forces_tot = f.calc_forces_total();
-        f.print_info();
-        result["forces"] = std::vector<std::vector<double>>(ctx.unit_cell().num_atoms(), std::vector<double>(3));
-        for (int i = 0; i < ctx.unit_cell().num_atoms(); i++) {
-            for (int j = 0; j < 3; j++) {
-                result["forces"][i][j] = forces_tot(j, i);
+        compute_forces = true;
+    }
+
+    nlohmann::json result;
+
+    Lattice_relaxation lr(dft);
+
+    switch (task_id) {
+        case task_t::ground_state_new:
+        case task_t::ground_state_restart: {
+            /* launch the calculation */
+            result = dft.find(inp.density_tol(), inp.energy_tol(), ctx.cfg().iterative_solver().energy_tolerance(),
+                    inp.num_dft_iter(), write_state);
+
+            if (compute_stress) {
+                dft.stress().calc_stress_total();
             }
+            if (compute_forces) {
+                dft.forces().calc_forces_total();
+            }
+            /* compute forces and stress */
+            if (ctx.cfg().control().print_stress() && !ctx.full_potential()) {
+                rte::ostream out(dft.ctx().out(), __func__);
+                dft.stress().print_info(out, dft.ctx().verbosity());
+                result["stress"] = std::vector<std::vector<double>>(3, std::vector<double>(3));
+                auto st = dft.stress().stress_total();
+                for (int i = 0; i < 3; i++) {
+                    for (int j = 0; j < 3; j++) {
+                        result["stress"][i][j] = st(j, i);
+                    }
+                }
+            }
+            if (ctx.cfg().control().print_forces()) {
+                rte::ostream out(dft.ctx().out(), __func__);
+                dft.forces().print_info(out, dft.ctx().verbosity());
+                result["forces"] = std::vector<std::vector<double>>(ctx.unit_cell().num_atoms(), std::vector<double>(3));
+                auto& ft = dft.forces().forces_total();
+                for (int i = 0; i < ctx.unit_cell().num_atoms(); i++) {
+                    for (int j = 0; j < 3; j++) {
+                        result["forces"][i][j] = ft(j, i);
+                    }
+                }
+            }
+            break;
+        }
+        case task_t::ground_state_new_relax: {
+            result = lr.find(ctx.cfg().vcsqnm().num_steps(), ctx.cfg().vcsqnm().forces_tol());
+            break;
+        }
+        case task_t::ground_state_new_vcrelax: {
+            result = lr.find(ctx.cfg().vcsqnm().num_steps(), ctx.cfg().vcsqnm().forces_tol(),
+                    ctx.cfg().vcsqnm().stress_tol());
+            break;
+        }
+        default: {
+            RTE_OUT(ctx.out()) << "task " << task_id << " is not handeled" << std::endl;
+            break;
         }
     }
 
     if (write_state && write_output) {
         json dict;
         json_output_common(dict);
 
-        dict["task"] = static_cast<int>(task);
+        dict["task"] = task_id;
         dict["context"] = ctx.serialize();
         dict["ground_state"] = result;
         //dict["timers"] = utils::timer::serialize();
@@ -282,7 +350,7 @@ ground_state(Simulation_context& ctx, task_t task, cmd_args const& args, int wri
 void run_tasks(cmd_args const& args)
 {
     /* get the task id */
-    task_t task = static_cast<task_t>(args.value<int>("task", 0));
+    int task_id = args.value<int>("task", 0);
 
     /* get the input file name */
     auto fpath = args.value<fs::path>("input", "sirius.json");
@@ -299,24 +367,26 @@ void run_tasks(cmd_args const& args)
     }
 
     auto fname = fpath.string();
-
-    if (task == task_t::ground_state_new || task == task_t::ground_state_restart) {
+    if (task_id == task_t::ground_state_new ||
+        task_id == task_t::ground_state_restart ||
+        task_id == task_t::ground_state_new_relax ||
+        task_id == task_t::ground_state_new_vcrelax) {
         auto ctx = create_sim_ctx(fname, args);
         ctx->initialize();
         //if (ctx->comm().rank() == 0) {
         //    auto dict = ctx->serialize();
         //    std::ofstream ofs("setup.json", std::ofstream::out | std::ofstream::trunc);
         //    ofs << dict.dump(4);
         //}
-        //if (ctx->full_potential()) {
-        //    ctx->gk_cutoff(ctx->aw_cutoff() / ctx->unit_cell().min_mt_radius());
-        //}
-        ground_state(*ctx, task, args, 1);
+        int write_output{1};
+        ground_state(*ctx, task_id, args, write_output);
     }
-    if (task == task_t::eos) {
+    if (task_id == task_t::eos) {
         auto s0 = std::pow(args.value<double>("volume_scale0"), 1.0 / 3);
         auto s1 = std::pow(args.value<double>("volume_scale1"), 1.0 / 3);
 
+        int write_output{0};
+
         int rank{0};
         int num_steps{10};
         std::vector<double> volume;
@@ -329,7 +399,7 @@ void run_tasks(cmd_args const& args)
             auto lv = ctx->unit_cell().lattice_vectors() * s;
             ctx->unit_cell().set_lattice_vectors(lv);
             ctx->initialize();
-            auto e = ground_state(*ctx, task_t::ground_state_new, args, 0);
+            auto e = ground_state(*ctx, task_t::ground_state_new, args, write_output);
             volume.push_back(ctx->unit_cell().omega());
             energy.push_back(e);
         }
@@ -340,7 +410,7 @@ void run_tasks(cmd_args const& args)
             }
         }
     }
-    if (task == task_t::read_config) {
+    if (task_id == task_t::read_config) {
         //int count{0};
         //while (true) {
         //    std::stringstream s;
@@ -367,10 +437,9 @@ void run_tasks(cmd_args const& args)
             ctx1.unit_cell().atom(ia).set_position(ctx2.unit_cell().atom(ia).position());
         }
         ctx1.update();
-
     }
 
-    if (task == task_t::k_point_path) {
+    if (task_id == task_t::k_point_path) {
         auto ctx = create_sim_ctx(fname, args);
         ctx->cfg().iterative_solver().energy_tolerance(1e-12);
         ctx->gamma_point(false);
@@ -431,7 +500,7 @@ void run_tasks(cmd_args const& args)
         if (!ctx->full_potential()) {
             band.initialize_subspace(ks, H0);
             if (ctx->hubbard_correction()) {
-                TERMINATE("fix me");
+                RTE_THROW("fix me");
                 //potential.U().compute_occupation_matrix(ks); // TODO: this is wrong; U matrix should come form the saved file
                 //potential.U().calculate_hubbard_potential_and_energy(potential.U().occupation_matrix());
             }
@@ -507,11 +576,6 @@ int main(int argn, char** argv)
     args.register_key("--volume_scale1=", "{double} final volume scale for EOS calculation");
 
     args.parse_args(argn, argv);
-    if (args.exist("help")) {
-        std::printf("Usage: %s [options]\n", argv[0]);
-        args.print_help();
-        return 0;
-    }
 
 #if defined(_GNU_SOURCE)
     if (args.exist("fpe")) {

apps/nlcg/sirius.nlcg.cpp

@@ -130,7 +130,7 @@ double ground_state(Simulation_context& ctx,
     if (ctx.cfg().control().print_stress() && !ctx.full_potential()) {
         Stress& s       = dft.stress();
         auto stress_tot = s.calc_stress_total();
-        s.print_info();
+        s.print_info(dft.ctx().out(), dft.ctx().verbosity());
         result["stress"] = std::vector<std::vector<double>>(3, std::vector<double>(3));
         for (int i = 0; i < 3; i++) {
             for (int j = 0; j < 3; j++) {
@@ -141,7 +141,7 @@ double ground_state(Simulation_context& ctx,
     if (ctx.cfg().control().print_forces()) {
         Force& f         = dft.forces();
         auto& forces_tot = f.calc_forces_total();
-        f.print_info();
+        f.print_info(dft.ctx().out(), dft.ctx().verbosity());
         result["forces"] = std::vector<std::vector<double>>(ctx.unit_cell().num_atoms(), std::vector<double>(3));
         for (int i = 0; i < ctx.unit_cell().num_atoms(); i++) {
             for (int j = 0; j < 3; j++) {

apps/nlcg/sirius.test.nlcg.cpp

@@ -101,7 +101,7 @@ double ground_state(Simulation_context& ctx,
     if (ctx.control().print_stress_ && !ctx.full_potential()) {
         Stress& s       = dft.stress();
         auto stress_tot = s.calc_stress_total();
-        s.print_info();
+        s.print_info(dft.ctx().out(), dft.ctx().verbosity());
         result["stress"] = std::vector<std::vector<double>>(3, std::vector<double>(3));
         for (int i = 0; i < 3; i++) {
             for (int j = 0; j < 3; j++) {
@@ -112,7 +112,7 @@ double ground_state(Simulation_context& ctx,
     if (ctx.control().print_forces_) {
         Force& f         = dft.forces();
         auto& forces_tot = f.calc_forces_total();
-        f.print_info();
+        f.print_info(dft.ctx().out(), dft.ctx().verbosity());
         result["forces"] = std::vector<std::vector<double>>(ctx.unit_cell().num_atoms(), std::vector<double>(3));
         for (int i = 0; i < ctx.unit_cell().num_atoms(); i++) {
             for (int j = 0; j < 3; j++) {

apps/tests/test_fortran_api.f90

@@ -122,6 +122,12 @@ program test_fortran_api
 enddo
 
 
+call sirius_option_get("control", "output", SIRIUS_STRING_TYPE, C_LOC(str_val), len(str_val))
+write(*,*)trim(adjustl(str_val))
+str_val = "file:output.txt"
+call sirius_option_set(handler, "control", "output", SIRIUS_STRING_TYPE, C_LOC(str_val), len(trim(adjustl(str_val))))
+
+
 call sirius_initialize_context(handler)
 !call sirius_print_info(handler)
 !

src/CMakeLists.txt

@@ -134,6 +134,7 @@ target_compile_definitions(sirius PUBLIC
   $<$<BOOL:${USE_ROCM}>:SIRIUS_GPU SIRIUS_ROCM>
   $<$<BOOL:${USE_VDWXC}>:SIRIUS_USE_VDWXC>
   $<$<BOOL:${USE_FP32_BOOL}>:USE_FP32>
+  $<$<BOOL:${USE_VCSQNM}>:SIRIUS_VCSQNM>
   $<$<BOOL:${HAVE_LIBVDW_WITH_MPI}>:SIRIUS_HAVE_VDWXC_MPI>
   $<$<AND:$<BOOL:${USE_MAGMA}>,$<BOOL:${USE_ROCM}>>:HAVE_HIP> # Required for magma headers
 )

src/README.md

@@ -5,6 +5,7 @@ SIRIUS_PRINT_MEMORY_USAGE
 SIRIUS_PRINT_PERFORMANCE
 SIRIUS_PRINT_CHECKSUM
 SIRIUS_PRINT_MPI_LAYOUT
+SIRIUS_PRINT_TIMING
 SIRIUS_EV_SOLVER
 SIRIUS_VERBOSITY
 SIRIUS_SAVE_CONFIG