Merge pull request #630 from electronic-structure/develop

Update master branch

CMakeLists.txt

@@ -1,7 +1,7 @@
 # We could be fine with 3.12, but 3.14 has ctest --show-only=json-v1 used in CI.
 cmake_minimum_required(VERSION 3.14)
 
-project(SIRIUS VERSION 7.2.3)
+project(SIRIUS VERSION 7.2.4)
 
 # user variables
 set(CREATE_PYTHON_MODULE    OFF CACHE BOOL "create sirius Python module")

apps/tests/test_srvo3_pwpp.f90

@@ -6,11 +6,14 @@ program test_fortran_api
 type(C_PTR) :: kset
 type(C_PTR) :: dft
 integer i, rank
-real(8) :: lat_vec(3,3), pos(3), forces(3, 5), stress(3,3), energy, scf_correction
+real(8) :: lat_vec(3,3), lat_vec1(3,3), pos(3), forces(3, 5), stress(3,3), energy, scf_correction
+real(8) :: alpha
 
 ! initialize the library
 call sirius_initialize(call_mpi_init=.true.)
 
+call MPI_COMM_RANK(MPI_COMM_WORLD, rank, i)
+
 ! create simulation context using a specified communicator
 call sirius_create_context(MPI_COMM_WORLD, handler)
 
@@ -20,7 +23,7 @@ program test_fortran_api
 
 ! atomic units are used everywhere
 ! plane-wave cutoffs are provided in a.u.^-1
-call sirius_set_parameters(handler, pw_cutoff=20.d0, gk_cutoff=7.d0)
+call sirius_set_parameters(handler, pw_cutoff=40.d0, gk_cutoff=7.d0)
 
 lat_vec = 0.d0
 do i = 1, 3
@@ -50,18 +53,60 @@ program test_fortran_api
 ! initialize the simulation handler
 call sirius_initialize_context(handler)
 
+! test different lattice updates
+do i = 1, 11, 2
+  alpha =  (1.0+(i-6)*0.3/5.0)
+  lat_vec1 = lat_vec
+  lat_vec1(1,1) = lat_vec(1,1) * alpha
+  lat_vec1(2,2) = lat_vec(2,2) / alpha
+  write(*,*)'new lattice:',lat_vec1
+  call sirius_set_lattice_vectors(handler, lat_vec1(:, 1), lat_vec1(:, 2), lat_vec1(:, 3))
+  call sirius_update_context(handler)
+enddo
+
+call sirius_set_lattice_vectors(handler, lat_vec(:, 1), lat_vec(:, 2), lat_vec(:, 3))
+call sirius_update_context(handler)
+
 ! create [2,2,2] k-grid with [0,0,0] offset and use symmetry to get irreducible set of k-points
 call sirius_create_kset_from_grid(handler, (/2, 2, 2/), (/0, 0, 0/), .true., kset)
 
 call sirius_create_ground_state(kset, dft)
 call sirius_find_ground_state(dft)
 
+if (rank.eq.0) then
+  write(*,*)'running new DFT with updated lattice vectors'
+endif
+call sirius_set_parameters(handler, iter_solver_tol=1d-2)
+lat_vec = lat_vec * 0.9d0
+call sirius_set_lattice_vectors(handler, lat_vec(:, 1), lat_vec(:, 2), lat_vec(:, 3))
+call sirius_update_context(handler)
+call sirius_find_ground_state(dft)
+
+if (rank.eq.0) then
+  write(*,*)'running new DFT with updated lattice vectors'
+endif
+call sirius_set_parameters(handler, iter_solver_tol=1d-2)
+lat_vec = lat_vec * 0.9d0
+call sirius_set_lattice_vectors(handler, lat_vec(:, 1), lat_vec(:, 2), lat_vec(:, 3))
+call sirius_update_context(handler)
+call sirius_find_ground_state(dft)
+
+if (rank.eq.0) then
+  write(*,*)'running new DFT with updated lattice vectors'
+endif
+call sirius_set_parameters(handler, iter_solver_tol=1d-2)
+lat_vec = lat_vec * 0.8d0
+call sirius_set_lattice_vectors(handler, lat_vec(:, 1), lat_vec(:, 2), lat_vec(:, 3))
+call sirius_update_context(handler)
+call sirius_find_ground_state(dft)
+
+
+
 call sirius_get_forces(dft, "total", forces)
 call sirius_get_stress_tensor(dft, "total", stress)
 call sirius_get_energy(dft, "total", energy)
 call sirius_get_energy(dft, "descf", scf_correction)
 
-call MPI_COMM_RANK(MPI_COMM_WORLD, rank, i)
 
 if (rank.eq.0) then
   write(*,*)'Forces:'

src/context/input.cpp

@@ -27,22 +27,19 @@
 
 namespace sirius {
 
-void Hubbard_input::read(json const& parser)
+void
+Hubbard_input::read(json const& parser)
 {
     if (!parser.count("hubbard")) {
         return;
     }
-
     auto section = parser["hubbard"];
 
-    orthogonalize_hubbard_orbitals_ =
-        section.value("orthogonalize_hubbard_wave_functions", orthogonalize_hubbard_orbitals_);
+    orthogonalize_hubbard_orbitals_ = section.value("orthogonalize", orthogonalize_hubbard_orbitals_);
 
-    normalize_hubbard_orbitals_ =
-        section.value("normalize_hubbard_wave_functions", normalize_hubbard_orbitals_);
+    normalize_hubbard_orbitals_ = section.value("normalize", normalize_hubbard_orbitals_);
 
-    simplified_hubbard_correction_ =
-        section.value("simplified_hubbard_correction", simplified_hubbard_correction_);
+    simplified_hubbard_correction_ = section.value("simplified", simplified_hubbard_correction_);
 
     if (section.count("projection_method")) {
         std::string projection_method = parser["hubbard"]["projection_method"].get<std::string>();
@@ -66,63 +63,77 @@ void Hubbard_input::read(json const& parser)
 
     auto v = parser["unit_cell"]["atom_types"].get<std::vector<std::string>>();
 
-    for (auto& label : v) {
-        if (section.count(label)) {
-            if (species_with_U.count(label)) {
-                throw std::runtime_error("U-correction for atom " + label + " has already been defined");
-            }
+    auto sec = parser["hubbard"]["local"];
 
-            hubbard_correction_ = true;
+    if (sec.size() == 0) {
+        throw std::runtime_error(
+            "The Hubbard correction section is defined but contain no information about atoms with U-correction");
+    }
 
-            hubbard_orbital_t ho;
+    for (int elem = 0; elem < sec.size(); elem++) {
+        std::string label;
+        label = sec[elem].value("atom_type", label);
+        if (species_with_U.count(label)) {
+            throw std::runtime_error("U-correction for atom " + label + " has already been defined");
+        }
 
-            ho.coeff[0] = section[label].value("U", ho.coeff[0]);
-            ho.coeff[1] = section[label].value("J", ho.coeff[1]);
-            ho.coeff[2] = section[label].value("B", ho.coeff[2]);
-            ho.coeff[2] = section[label].value("E2", ho.coeff[2]);
-            ho.coeff[3] = section[label].value("E3", ho.coeff[3]);
-            ho.coeff[4] = section[label].value("alpha", ho.coeff[4]);
-            ho.coeff[5] = section[label].value("beta", ho.coeff[5]);
+        // check that the atom type is actually defined
 
-            /* now convert eV in Ha */
-            for (int s = 0; s < 6; s++) {
-                ho.coeff[s] /= ha2ev;
-            }
-            ho.l = section[label].value("l", ho.l);
-            ho.n = section[label].value("n", ho.n);
-
-            if (ho.l == -1 || ho.n == -1) {
-                std::string level;
-                section[label].value("hubbard_orbital", level);
-                std::map<char, int> const map_l = {{'s', 0},{'p', 1}, {'d', 2}, {'f', 3}};
-
-                std::istringstream iss(std::string(1, level[0]));
-                iss >> ho.n;
-                if (ho.n <= 0 || iss.fail()) {
-                    std::stringstream s;
-                    s << "wrong principal quantum number : " << std::string(1, level[0]);
-                    throw std::runtime_error(s.str());
-                }
-                ho.l = map_l.at(level[1]);
-            }
+        auto found = std::find(v.begin(), v.end(), label);
 
-            if (!section[label].count("occupancy")) {
-                throw std::runtime_error("initial occupancy of the Hubbard orbital is not set");
-            }
-            ho.occupancy = section[label].value("occupancy", ho.occupancy);
-            ho.initial_occupancy = section[label].value("initial_occupancy", ho.initial_occupancy);
+        if (found == v.end()) {
+            throw std::runtime_error("The atom type " + label + " can not be found in the unit cell declaration");
+        }
 
-            int sz = static_cast<int>(ho.initial_occupancy.size());
-            int lmmax = 2 * ho.l + 1;
+        hubbard_correction_ = true;
 
-            if (!(sz == 0 || sz == lmmax || sz == 2 * lmmax)) {
+        hubbard_orbital_t ho;
+
+        ho.coeff[0] = sec[elem].value("U", ho.coeff[0]);
+        ho.coeff[1] = sec[elem].value("J", ho.coeff[1]);
+        ho.coeff[2] = sec[elem].value("B", ho.coeff[2]);
+        ho.coeff[2] = sec[elem].value("E2", ho.coeff[2]);
+        ho.coeff[3] = sec[elem].value("E3", ho.coeff[3]);
+        ho.coeff[4] = sec[elem].value("alpha", ho.coeff[4]);
+        ho.coeff[5] = sec[elem].value("beta", ho.coeff[5]);
+
+        /* now convert eV in Ha */
+        for (int s = 0; s < 6; s++) {
+            ho.coeff[s] /= ha2ev;
+        }
+        ho.l = sec[elem].value("l", ho.l);
+        ho.n = sec[elem].value("n", ho.n);
+        if (ho.l == -1 || ho.n == -1) {
+            std::string level;
+            sec[elem].value("hubbard_orbital", level);
+            std::map<char, int> const map_l = {{'s', 0}, {'p', 1}, {'d', 2}, {'f', 3}};
+
+            std::istringstream iss(std::string(1, level[0]));
+            iss >> ho.n;
+            if (ho.n <= 0 || iss.fail()) {
                 std::stringstream s;
-                s << "wrong size of initial occupacies vector (" << sz << ") for l = " << ho.l;
+                s << "wrong principal quantum number : " << std::string(1, level[0]);
                 throw std::runtime_error(s.str());
             }
+            ho.l = map_l.at(level[1]);
+        }
 
-            species_with_U[label] = ho;
+        if (!sec[elem].count("occupancy")) {
+            throw std::runtime_error("initial occupancy of the Hubbard orbital is not set");
         }
+        ho.occupancy         = sec[elem].value("occupancy", ho.occupancy);
+        ho.initial_occupancy = sec[elem].value("initial_occupancy", ho.initial_occupancy);
+
+        int sz    = static_cast<int>(ho.initial_occupancy.size());
+        int lmmax = 2 * ho.l + 1;
+
+        if (!(sz == 0 || sz == lmmax || sz == 2 * lmmax)) {
+            std::stringstream s;
+            s << "wrong size of initial occupacies vector (" << sz << ") for l = " << ho.l;
+            throw std::runtime_error(s.str());
+        }
+
+        species_with_U[label] = ho;
     }
 }
 

src/context/input_schema.json

@@ -451,7 +451,7 @@
                     "default" : [],
                     "title" : "List of XC functionals (typically contains exchange term and correlation term).",
                     "description" : "Naming convention of LibXC is used, names should be provided in capital letters.",
-                    "example" : ["XC_LDA_X", "XC_LDA_C_PZ"] 
+                    "example" : ["XC_LDA_X", "XC_LDA_C_PZ"]
                 },
                 "core_relativity" : {
                     "type" : "string",
@@ -748,13 +748,21 @@
                                 "minItems" : 2,
                                 "maxItems" : 2
                             },
+                            "T" : {
+                                "type" : "array",
+                                "items" : {
+                                    "type" : "integer"
+                                },
+                                "minItems" : 3,
+                                "maxItems" : 3
+                            },
                             "l" : {
                                 "type" : "array",
                                 "items" : {
                                     "type" : "integer"
                                 },
                                 "minItems" : 2,
-                                "maxItems" : 2 
+                                "maxItems" : 2
                             },
                             "V" : {
                                 "type" : "number",
@@ -767,4 +775,3 @@
         }
     }
 }
-