feat(input, kpoints): support shifted kspacing mesh with koffset and kmesh_type (#7247)

* feat(input, kpoints): support shifted kspacing mesh with koffset and kmesh_type

* revert input-main.md

* Fix headers

* Use sync_string instead of sync_string to be case insensitive

* add test case for kspacing koffset kmesh_type

* add test case for kspacing koffset kmesh_type

* clean test case for kspacing koffset kmesh_type

* Update CASES_CPU

---------

Co-authored-by: Mohan Chen <mohanchen@pku.edu.cn>

Author: Cstandardlib

docs/parameters.yaml

@@ -224,6 +224,22 @@ parameters:
     default_value: "0.0"
     unit: ""
     availability: ""
+  - name: koffset
+    category: System variables
+    type: Vector of Real (3 values)
+    description: |
+      Set offsets for automatic k-point mesh generated by kspacing, in each reciprocal direction. This parameter is only effective when kspacing > 0.0 and gamma_only is false.
+    default_value: 0.0 0.0 0.0
+    unit: ""
+    availability: ""
+  - name: kmesh_type
+    category: System variables
+    type: String
+    description: |
+      Set mesh type used for automatic k-point mesh generated by kspacing. Available options are gamma and mp. This parameter is only effective when kspacing > 0.0 and gamma_only is false.
+    default_value: gamma
+    unit: ""
+    availability: ""
   - name: min_dist_coef
     category: System variables
     type: Real
@@ -1690,10 +1706,14 @@ parameters:
     availability: ""
   - name: cal_syns
     category: Molecular dynamics
-    type: Boolean
+    type: "Boolean [Integer](optional)"
     description: |
       Whether to calculate and output asynchronous overlap matrix for Hefei-NAMD interface. When enabled, calculates <phi(t-1)|phi(t)> by computing overlap between basis functions at atomic positions from previous time step and current time step. The overlap is calculated by shifting atom positions backward by velocity x md_dt. Output file: OUT.*/syns_nao.csr in CSR format.
 
+      * 0 or false: disable
+      * 1 or true: enable with default precision (8 digits)
+      * 1 5: enable with custom precision (5 digits)
+
       [NOTE] Only works with LCAO basis and molecular dynamics calculations. Requires atomic velocities. Output starts from the second MD step (istep > 0).
     default_value: "False"
     unit: ""
@@ -1715,6 +1735,7 @@ parameters:
       * vw: von Weizsacker (vW) functional
       * tf+: TF + vW functional
       * wt: Wang-Teter (WT) functional
+      * ext-wt: Extended Wang-Teter functional
       * xwm: XWM functional
       * lkt: Luo-Karasiev-Trickey (LKT) functional
       * ml: Machine learning KEDF
@@ -1768,31 +1789,39 @@ parameters:
       Weight of TF KEDF (kinetic energy density functional).
     default_value: "1.0"
     unit: ""
-    availability: "OFDFT with of_kinetic=tf, tf+, wt, xwm"
+    availability: "OFDFT with of_kinetic=tf, tf+, wt, ext-wt, xwm"
   - name: of_vw_weight
     category: "OFDFT: orbital free density functional theory"
     type: Real
     description: |
       Weight of vW KEDF (kinetic energy density functional).
     default_value: "1.0"
     unit: ""
-    availability: "OFDFT with of_kinetic=vw, tf+, wt, lkt, xwm"
+    availability: "OFDFT with of_kinetic=vw, tf+, wt, ext-wt, lkt, xwm"
   - name: of_wt_alpha
     category: "OFDFT: orbital free density functional theory"
     type: Real
     description: |
       Parameter alpha of WT KEDF (kinetic energy density functional).
     default_value: ""
     unit: ""
-    availability: OFDFT with of_kinetic=wt
+    availability: "OFDFT with of_kinetic=wt, ext-wt"
   - name: of_wt_beta
     category: "OFDFT: orbital free density functional theory"
     type: Real
     description: |
       Parameter beta of WT KEDF (kinetic energy density functional).
     default_value: ""
     unit: ""
-    availability: OFDFT with of_kinetic=wt
+    availability: "OFDFT with of_kinetic=wt, ext-wt"
+  - name: of_extwt_kappa
+    category: "OFDFT: orbital free density functional theory"
+    type: Real
+    description: |
+      Parameter kappa for EXT-WT KEDF.
+    default_value: "1.0 / (2.0 * std::pow(4./3., 1./3.) - 1.0)"
+    unit: ""
+    availability: OFDFT with of_kinetic=ext-wt
   - name: of_wt_rho0
     category: "OFDFT: orbital free density functional theory"
     type: Real

source/source_cell/klist.cpp

@@ -247,8 +247,16 @@ bool K_Vectors::read_kpoints(const UnitCell& ucell,
         std::ofstream ofs(fn.c_str());
         ofs << "K_POINTS" << std::endl;
         ofs << "0" << std::endl;
-        ofs << "Gamma" << std::endl;
-        ofs << nk1 << " " << nk2 << " " << nk3 << " 0 0 0" << std::endl;
+        if (PARAM.inp.kmesh_type == "mp")
+        {
+            ofs << "Monkhorst-Pack" << std::endl;
+        }
+        else
+        {
+            ofs << "Gamma" << std::endl;
+        }
+        ofs << nk1 << " " << nk2 << " " << nk3 << " " << PARAM.inp.koffset[0] << " " << PARAM.inp.koffset[1] << " "
+            << PARAM.inp.koffset[2] << std::endl;
         ofs.close();
     }
 

source/source_cell/test/klist_test.cpp

@@ -313,6 +313,10 @@ TEST_F(KlistTest, ReadKpointsKspacing)
     PARAM.input.kspacing[0] = 0.052918; // 0.52918/Bohr = 1/A
     PARAM.input.kspacing[1] = 0.052918; // 0.52918/Bohr = 1/A
     PARAM.input.kspacing[2] = 0.052918; // 0.52918/Bohr = 1/A
+    PARAM.input.kmesh_type = "gamma";
+    PARAM.input.koffset[0] = 0.0;
+    PARAM.input.koffset[1] = 0.0;
+    PARAM.input.koffset[2] = 0.0;
     setucell();
     std::string k_file = "./support/KPT3";
     kv->read_kpoints(ucell,k_file);
@@ -328,6 +332,10 @@ TEST_F(KlistTest, ReadKpointsKspacing3values)
     PARAM.input.kspacing[0] = 0.052918; // 0.52918/Bohr = 1/A
     PARAM.input.kspacing[1] = 0.06;     // 0.52918/Bohr = 1/A
     PARAM.input.kspacing[2] = 0.07;     // 0.52918/Bohr = 1/A
+    PARAM.input.kmesh_type = "gamma";
+    PARAM.input.koffset[0] = 0.0;
+    PARAM.input.koffset[1] = 0.0;
+    PARAM.input.koffset[2] = 0.0;
     setucell();
     std::string k_file = "./support/KPT3";
     kv->read_kpoints(ucell,k_file);
@@ -343,6 +351,10 @@ TEST_F(KlistTest, ReadKpointsInvalidKspacing3values)
     PARAM.input.kspacing[0] = 0.052918; // 0.52918/Bohr = 1/A
     PARAM.input.kspacing[1] = 0;        // 0.52918/Bohr = 1/A
     PARAM.input.kspacing[2] = 0.07;     // 0.52918/Bohr = 1/A
+    PARAM.input.kmesh_type = "gamma";
+    PARAM.input.koffset[0] = 0.0;
+    PARAM.input.koffset[1] = 0.0;
+    PARAM.input.koffset[2] = 0.0;
     std::string k_file = "./support/KPT3";
     testing::internal::CaptureStdout();
     EXPECT_EXIT(kv->read_kpoints(ucell,k_file), ::testing::ExitedWithCode(1), "");
@@ -352,6 +364,72 @@ TEST_F(KlistTest, ReadKpointsInvalidKspacing3values)
     PARAM.input.kspacing[2] = 0.0;
 }
 
+TEST_F(KlistTest, ReadKpointsKspacingShiftedGamma)
+{
+    kv->nspin = 1;
+    PARAM.input.kspacing[0] = 0.052918; // 0.52918/Bohr = 1/A
+    PARAM.input.kspacing[1] = 0.052918;
+    PARAM.input.kspacing[2] = 0.052918;
+    PARAM.input.kmesh_type = "gamma";
+    PARAM.input.koffset[0] = 0.5;
+    PARAM.input.koffset[1] = 0.5;
+    PARAM.input.koffset[2] = 0.5;
+    setucell();
+
+    std::string k_file = "./support/KPT3";
+    kv->read_kpoints(ucell, k_file);
+
+    EXPECT_EQ(kv->get_nkstot(), 343);
+    EXPECT_EQ(kv->get_k_kword(), "Gamma");
+    EXPECT_DOUBLE_EQ(kv->get_koffset(0), 0.5);
+    EXPECT_DOUBLE_EQ(kv->get_koffset(1), 0.5);
+    EXPECT_DOUBLE_EQ(kv->get_koffset(2), 0.5);
+    EXPECT_NEAR(kv->kvec_d[0].x, 1.0 / 14.0, 1e-12);
+    EXPECT_NEAR(kv->kvec_d[0].y, 1.0 / 14.0, 1e-12);
+    EXPECT_NEAR(kv->kvec_d[0].z, 1.0 / 14.0, 1e-12);
+
+    PARAM.input.kspacing[0] = 0.0;
+    PARAM.input.kspacing[1] = 0.0;
+    PARAM.input.kspacing[2] = 0.0;
+    PARAM.input.koffset[0] = 0.0;
+    PARAM.input.koffset[1] = 0.0;
+    PARAM.input.koffset[2] = 0.0;
+    PARAM.input.kmesh_type = "gamma";
+}
+
+TEST_F(KlistTest, ReadKpointsKspacingShiftedMP)
+{
+    kv->nspin = 1;
+    PARAM.input.kspacing[0] = 0.052918; // 0.52918/Bohr = 1/A
+    PARAM.input.kspacing[1] = 0.052918;
+    PARAM.input.kspacing[2] = 0.052918;
+    PARAM.input.kmesh_type = "mp";
+    PARAM.input.koffset[0] = 0.5;
+    PARAM.input.koffset[1] = 0.5;
+    PARAM.input.koffset[2] = 0.5;
+    setucell();
+
+    std::string k_file = "./support/KPT3";
+    kv->read_kpoints(ucell, k_file);
+
+    EXPECT_EQ(kv->get_nkstot(), 343);
+    EXPECT_EQ(kv->get_k_kword(), "Monkhorst-Pack");
+    EXPECT_DOUBLE_EQ(kv->get_koffset(0), 0.5);
+    EXPECT_DOUBLE_EQ(kv->get_koffset(1), 0.5);
+    EXPECT_DOUBLE_EQ(kv->get_koffset(2), 0.5);
+    EXPECT_NEAR(kv->kvec_d[0].x, -5.5 / 14.0, 1e-12);
+    EXPECT_NEAR(kv->kvec_d[0].y, -5.5 / 14.0, 1e-12);
+    EXPECT_NEAR(kv->kvec_d[0].z, -5.5 / 14.0, 1e-12);
+
+    PARAM.input.kspacing[0] = 0.0;
+    PARAM.input.kspacing[1] = 0.0;
+    PARAM.input.kspacing[2] = 0.0;
+    PARAM.input.koffset[0] = 0.0;
+    PARAM.input.koffset[1] = 0.0;
+    PARAM.input.koffset[2] = 0.0;
+    PARAM.input.kmesh_type = "gamma";
+}
+
 TEST_F(KlistTest, ReadKpointsGamma)
 {
     std::string k_file = "./support/KPT";

source/source_io/module_parameter/input_parameter.h

@@ -64,6 +64,8 @@ struct Input_para
     int diago_proc = 0;                             ///< the number of procs used to diag. mohan add 2012-01-13
     int nbspline = -1;                              ///< the order of B-spline basis(>=0) if it is -1 (default)
     std::vector<double> kspacing = {0.0, 0.0, 0.0}; ///< kspacing for k-point generation
+    std::vector<double> koffset = {0.0, 0.0, 0.0};  ///< koffset for kspacing-generated k-point mesh
+    std::string kmesh_type = "gamma";               ///< k-point mesh type for kspacing-generated k-point mesh: gamma or mp
     double min_dist_coef = 0.2;                     ///< allowed minimum distance between two atoms
 
     std::string device = "auto";

source/source_io/module_parameter/read_input_item_system.cpp

@@ -4,8 +4,8 @@
 #include "read_input_tool.h"
 #include "source_base/module_device/device.h"
 
-#include <fstream>
 #include <unistd.h>
+#include <algorithm>
 
 namespace ModuleIO
 {
@@ -613,6 +613,53 @@ Available options are:
         };
         this->add_item(item);
     }
+    {
+        Input_Item item("koffset");
+        item.annotation = "offset for kspacing-generated automatic k-point mesh";
+        item.category = "System variables";
+        item.type = "Vector of Real (3 values)";
+        item.description = "Set offsets for automatic k-point mesh generated by kspacing, in each reciprocal direction. "
+                          "This parameter is only effective when kspacing > 0.0 and gamma_only is false.";
+        item.default_value = "0.0 0.0 0.0";
+        item.read_value = [](const Input_Item& item, Parameter& para) {
+            if (item.get_size() != 3)
+            {
+                ModuleBase::WARNING_QUIT("ReadInput", "koffset must provide three values.");
+            }
+            para.input.koffset[0] = std::stod(item.str_values[0]);
+            para.input.koffset[1] = std::stod(item.str_values[1]);
+            para.input.koffset[2] = std::stod(item.str_values[2]);
+        };
+        sync_doublevec(input.koffset, 3, 0.0);
+        this->add_item(item);
+    }
+    {
+        Input_Item item("kmesh_type");
+        item.annotation = "mesh type for kspacing-generated automatic k-point mesh";
+        item.category = "System variables";
+        item.type = "String";
+        item.description = "Set mesh type used for automatic k-point mesh generated by kspacing. "
+                          "Available options are gamma and mp. This parameter is only effective when kspacing > 0.0 "
+                          "and gamma_only is false.";
+        item.default_value = "gamma";
+        item.read_value = [](const Input_Item& item, Parameter& para) {
+            para.input.kmesh_type = strvalue;
+            std::transform(para.input.kmesh_type.begin(),
+                           para.input.kmesh_type.end(),
+                           para.input.kmesh_type.begin(),
+                           ::tolower);
+        };
+        sync_string(input.kmesh_type);
+        item.check_value = [](const Input_Item& item, const Parameter& para) {
+            std::vector<std::string> avail_list = {"gamma", "mp"};
+            if (std::find(avail_list.begin(), avail_list.end(), para.input.kmesh_type) == avail_list.end())
+            {
+                const std::string warningstr = nofound_str(avail_list, "kmesh_type");
+                ModuleBase::WARNING_QUIT("ReadInput", warningstr);
+            }
+        };
+        this->add_item(item);
+    }
     {
         Input_Item item("min_dist_coef");
         item.annotation = "factor related to the allowed minimum distance "

tests/01_PW/210_PW_kspace_shift/INPUT

@@ -0,0 +1,23 @@
+INPUT_PARAMETERS 
+suffix            autotest
+calculation       scf
+init_wfc          random
+
+symmetry          1
+ecutwfc           5
+scf_thr           1e-06
+scf_nmax          30
+basis_type        pw
+latname           bco
+
+mixing_type       broyden
+mixing_beta       0.4
+
+smearing_method   gaussian
+smearing_sigma    0.02
+
+kspacing          0.6
+koffset           0.5 0.5 0.5
+kmesh_type        gamma
+pseudo_dir        ../../PP_ORB
+pw_seed           1

tests/01_PW/210_PW_kspace_shift/README

@@ -0,0 +1,2 @@
+fast integration case for kspacing-generated shifted Gamma mesh
+covers new INPUT parameters: koffset and kmesh_type

tests/01_PW/210_PW_kspace_shift/STRU

@@ -0,0 +1,22 @@
+#This is the atom file containing all the information
+#about the lattice structure.
+
+ATOMIC_SPECIES
+H 1.0008   H_ONCV_PBE-1.0.upf
+ 
+NUMERICAL_ORBITAL
+H_gga_6au_60Ry_2s1p.orb
+
+LATTICE_CONSTANT
+10.0  			#Lattice constant
+
+LATTICE_PARAMETERS
+1.5 2.0
+
+ATOMIC_POSITIONS
+Cartesian   		#Cartesian(Unit is LATTICE_CONSTANT)
+H 			#Name of element	
+0.0			#Magnetic for this element.
+2			#Number of atoms
+0.00 0.00 -0.0661400 0 0 0	#x,y,z, move_x, move_y, move_z
+0.00 0.00  0.0661400 0 0 0	#x,y,z, move_x, move_y, move_z

tests/01_PW/210_PW_kspace_shift/result.ref

@@ -0,0 +1,6 @@
+etotref -27.08734962492505
+etotperatomref -13.5436748125
+pointgroupref D_2h
+spacegroupref D_2h
+nksibzref 2
+totaltimeref 0.28

tests/01_PW/CASES_CPU.txt

@@ -111,6 +111,7 @@ scf_out_chg_tau
 207_PW_skip
 208_PW_CG_float
 209_PW_DFTHALF
+210_PW_kspace_shift
 801_PW_LT_sc
 802_PW_LT_fcc
 803_PW_LT_bcc