Dfpt csj

CMakeLists.txt

@@ -830,6 +830,7 @@ target_link_libraries(
   module_pwdft
   module_ofdft
   module_stodft
+  module_dfpt
   psi
   psi_initializer
   psi_overall_init

examples/18_md/02_lcao_output/INPUT

@@ -0,0 +1,38 @@
+INPUT_PARAMETERS
+#Parameters (1.General)
+suffix                 Si_nhc_nvt
+calculation            md
+nbands                 20
+symmetry               0
+pseudo_dir             ../../../tests/PP_ORB
+orbital_dir            ../../../tests/PP_ORB
+
+#Parameters (2.Iteration)
+ecutwfc                10
+scf_thr                1e-5
+scf_nmax               100
+
+#Parameters (3.Basis)
+basis_type             lcao
+ks_solver              genelpa
+gamma_only             1
+
+#Parameters (4.Smearing)
+smearing_method        gaussian
+smearing_sigma         0.001
+
+#Parameters (5.Mixing)
+mixing_type            broyden
+mixing_beta            0.3
+chg_extrap             second-order
+
+#Parameters (6.MD)
+md_type                nvt
+md_nstep               10
+md_dt                  1
+md_tfirst              300
+md_tfreq               1
+md_tchain              1
+
+dft_functional       hse
+restart_save           1

examples/18_md/02_lcao_output/KPT

@@ -0,0 +1,4 @@
+K_POINTS
+0
+Gamma
+1 1 1 0 0 0

examples/18_md/02_lcao_output/STRU

@@ -0,0 +1,22 @@
+ATOMIC_SPECIES
+Si  28.085  Si_ONCV_PBE-1.0.upf
+
+NUMERICAL_ORBITAL
+Si_gga_8au_60Ry_2s2p1d.orb
+
+LATTICE_CONSTANT
+1.8897270    # 1 Angstrom = 1.8897270 bohr
+
+LATTICE_VECTORS
+5.43090      0.00000      0.00000
+0.00000      5.43090      0.00000
+0.00000      0.00000      5.43090
+
+ATOMIC_POSITIONS
+Direct 
+
+Si
+0.0
+2
+0.000	0.000	0.000	1	1	1
+0.5 0.5 0.5 1 1 1

source/source_cell/CMakeLists.txt

@@ -30,6 +30,8 @@ add_library(
     k_vector_utils.cpp
     sep.cpp
     sep_cell.cpp
+    qlist.cpp
+    qlist.h
 )
 
 if(ENABLE_COVERAGE)

source/source_cell/k_vector_utils.cpp

@@ -1,9 +1,11 @@
 //
 // Created by rhx on 25-6-3.
+// Added by Shengjun Chen on 2026-05-26 with Q_Vectors class and related functions.
 //
 #include "k_vector_utils.h"
 
 #include "klist.h"
+#include "qlist.h"
 #include "source_base/global_variable.h"
 #include "source_base/matrix3.h"
 
@@ -59,6 +61,47 @@ void kvec_d2c(K_Vectors& kv, const ModuleBase::Matrix3& reciprocal_vec)
         }
     }
 }
+
+void kvec_d2c(Q_Vectors& qv, const ModuleBase::Matrix3& reciprocal_vec)
+{
+    if (qv.qvec_d.size() != qv.qvec_c.size())
+    {
+        qv.qvec_c.resize(qv.qvec_d.size());
+    }
+    int nqs = qv.qvec_d.size(); // always convert all q vectors
+
+    for (int i = 0; i < nqs; i++)
+    {
+        if (std::abs(qv.qvec_d[i].x) < 1.0e-10)
+        {
+            qv.qvec_d[i].x = 0.0;
+        }
+        if (std::abs(qv.qvec_d[i].y) < 1.0e-10)
+        {
+            qv.qvec_d[i].y = 0.0;
+        }
+        if (std::abs(qv.qvec_d[i].z) < 1.0e-10)
+        {
+            qv.qvec_d[i].z = 0.0;
+        }
+
+        qv.qvec_c[i] = qv.qvec_d[i] * reciprocal_vec;
+
+        if (std::abs(qv.qvec_c[i].x) < 1.0e-10)
+        {
+            qv.qvec_c[i].x = 0.0;
+        }
+        if (std::abs(qv.qvec_c[i].y) < 1.0e-10)
+        {
+            qv.qvec_c[i].y = 0.0;
+        }
+        if (std::abs(qv.qvec_c[i].z) < 1.0e-10)
+        {
+            qv.qvec_c[i].z = 0.0;
+        }
+    }
+}
+
 void kvec_c2d(K_Vectors& kv, const ModuleBase::Matrix3& latvec)
 {
     if (kv.kvec_d.size() != kv.kvec_c.size())
@@ -80,6 +123,21 @@ void kvec_c2d(K_Vectors& kv, const ModuleBase::Matrix3& latvec)
     }
 }
 
+void kvec_c2d(Q_Vectors& qv, const ModuleBase::Matrix3& latvec)
+{
+    if (qv.qvec_d.size() != qv.qvec_c.size())
+    {
+        qv.qvec_d.resize(qv.qvec_c.size());
+    }
+    int nqs = qv.qvec_d.size(); // always convert all q vectors
+
+    ModuleBase::Matrix3 RT = latvec.Transpose();
+    for (int i = 0; i < nqs; i++)
+    {
+        qv.qvec_d[i] = qv.qvec_c[i] * RT;
+    }
+}
+
 void set_both_kvec(K_Vectors& kv, const ModuleBase::Matrix3& G, const ModuleBase::Matrix3& R, std::string& skpt)
 {
     if (true) // Originally GlobalV::FINAL_SCF, but we don't have this variable in the new code.
@@ -145,6 +203,71 @@ void set_both_kvec(K_Vectors& kv, const ModuleBase::Matrix3& G, const ModuleBase
     return;
 }
 
+void set_both_kvec(Q_Vectors& qv, const ModuleBase::Matrix3& G, const ModuleBase::Matrix3& R, std::string& sqpt)
+{
+    if (true) // Originally GlobalV::FINAL_SCF, but we don't have this variable in the new code.
+    {
+        if (qv.get_q_nqstot() == 0)
+        {
+            qv.qd_done = true;
+            qv.qc_done = false;
+        }
+        else
+        {
+            if (qv.get_q_qword() == "Cartesian" || qv.get_q_qword() == "C")
+            {
+                qv.qc_done = true;
+                qv.qd_done = false;
+            }
+            else if (qv.get_q_qword() == "Direct" || qv.get_q_qword() == "D")
+            {
+                qv.qd_done = true;
+                qv.qc_done = false;
+            }
+            else
+            {
+                GlobalV::ofs_warning << " Error : neither Cartesian nor Direct qpoint." << std::endl;
+            }
+        }
+    }
+
+    // set cartesian q vectors.
+    if (!qv.qc_done && qv.qd_done)
+    {
+        KVectorUtils::kvec_d2c(qv, G);
+        qv.qc_done = true;
+    }
+
+    // set direct q vectors
+    else if (qv.qc_done && !qv.qd_done)
+    {
+        KVectorUtils::kvec_c2d(qv, R);
+        qv.qd_done = true;
+    }
+    std::string table;
+    table += " Q-POINTS DIRECT COORDINATES\n";
+    table += FmtCore::format("%8s%12s%12s%12s%8s\n", "QPOINTS", "DIRECT_X", "DIRECT_Y", "DIRECT_Z", "WEIGHT");
+    for (int i = 0; i < qv.get_q_nqstot(); i++)
+    {
+        table += FmtCore::format("%8d%12.8f%12.8f%12.8f%8.4f\n",
+                                 i + 1,
+                                 qv.qvec_d[i].x,
+                                 qv.qvec_d[i].y,
+                                 qv.qvec_d[i].z,
+                                 qv.wq[i]);
+    }
+    GlobalV::ofs_running << table << std::endl;
+    if (GlobalV::MY_RANK == 0)
+    {
+        std::stringstream ss;
+        ss << " " << std::setw(40) << "nqstot now"
+           << " = " << qv.get_q_nqstot() << std::endl;
+        ss << table << std::endl;
+        sqpt = ss.str();
+    }
+    return;
+}
+
 void set_after_vc(K_Vectors& kv, const int& nspin_in, const ModuleBase::Matrix3& reciprocal_vec)
 {
     GlobalV::ofs_running << "\n SETUP K-POINTS" << std::endl;
@@ -175,6 +298,9 @@ void set_after_vc(K_Vectors& kv, const int& nspin_in, const ModuleBase::Matrix3&
     print_klists(kv, GlobalV::ofs_running);
 }
 
+// Don't need to set_after_vc for Q_Vectors, since the q-points are only
+// used in the calculation of the wavefunctions.
+
 void print_klists(const K_Vectors& kv, std::ofstream& ofs)
 {
     ModuleBase::TITLE("KVectorUtils", "print_klists");
@@ -217,6 +343,47 @@ void print_klists(const K_Vectors& kv, std::ofstream& ofs)
     return;
 }
 
+void print_klists(const Q_Vectors& qv, std::ofstream& ofs)
+{
+    ModuleBase::TITLE("KVectorUtils", "print_qlists");
+    int nqs = qv.get_nqs();
+    int nqstot = qv.get_q_nqstot();
+
+    if (nqstot < nqs)
+    {
+        std::cout << "\n nqstot=" << nqstot;
+        std::cout << "\n nqs=" << nqs;
+        ModuleBase::WARNING_QUIT("print_qlists", "nqstot < nqs");
+    }
+    std::string table;
+    table += " Q-POINTS CARTESIAN COORDINATES\n";
+    table += FmtCore::format("%8s%12s%12s%12s%8s\n", "QPOINTS", "CARTESIAN_X", "CARTESIAN_Y", "CARTESIAN_Z", "WEIGHT");
+    for (int i = 0; i < nqs; i++)
+    {
+        table += FmtCore::format("%8d%12.8f%12.8f%12.8f%8.4f\n",
+                                 i + 1,
+                                 qv.qvec_c[i].x,
+                                 qv.qvec_c[i].y,
+                                 qv.qvec_c[i].z,
+                                 qv.wq[i]);
+    }
+    GlobalV::ofs_running << "\n" << table << std::endl;
+
+    table.clear();
+    table += " Q-POINTS DIRECT COORDINATES\n";
+    table += FmtCore::format("%8s%12s%12s%12s%8s\n", "QPOINTS", "DIRECT_X", "DIRECT_Y", "DIRECT_Z", "WEIGHT");
+    for (int i = 0; i < nqs; i++)
+    {
+        table += FmtCore::format("%8d%12.8f%12.8f%12.8f%8.4f\n",
+                                 i + 1,
+                                 qv.qvec_d[i].x,
+                                 qv.qvec_d[i].y,
+                                 qv.qvec_d[i].z,
+                                 qv.wq[i]);
+    }
+    GlobalV::ofs_running << "\n" << table << std::endl;
+    return
+
 #ifdef __MPI
 void kvec_mpi_k(K_Vectors& kv)
 {
@@ -349,6 +516,9 @@ void kvec_mpi_k(K_Vectors& kv)
 } // END SUBROUTINE
 #endif
 
+#ifdef __MPI
+// THIS Should be completed
+#endif
 
 void kvec_ibz_kpoint(K_Vectors& kv,
                      const ModuleSymmetry::Symmetry& symm,
@@ -807,4 +977,15 @@ void kvec_ibz_kpoint(K_Vectors& kv,
 
     return;
 }
+
+void kvec_ibz_kpoint(Q_Vectors& qv,
+                     const ModuleSymmetry::Symmetry& symm,
+                     bool use_symm,
+                     std::string& skpt,
+                     const UnitCell& ucell,
+                     bool& match)
+{
+// 大修特修
+}
+
 } // namespace KVectorUtils

source/source_cell/k_vector_utils.h

@@ -1,5 +1,6 @@
 //
 // Created by rhx on 25-6-3.
+// Added by Shengjun Chen on 2026-05-26 with Q_Vectors class and related functions.
 //
 
 #ifndef K_VECTOR_UTILS_H
@@ -9,6 +10,7 @@
 #include "source_cell/unitcell.h"
 
 class K_Vectors;
+class Q_Vectors;
 
 namespace KVectorUtils
 {
@@ -122,6 +124,26 @@ void kvec_ibz_kpoint(K_Vectors& kv,
                      std::string& skpt,
                      const UnitCell& ucell,
                      bool& match);
+
+
+// Function Overloading for Q_Vectors
+// Added for Q_Vectors, by Shengjun Chen on 2026-05-26.
+void kvec_d2c(Q_Vectors& qv, const ModuleBase::Matrix3& reciprocal_vec);
+void kvec_c2d(Q_Vectors& qv, const ModuleBase::Matrix3& latvec);
+void set_both_kvec(Q_Vectors& qv, const ModuleBase::Matrix3& G, const ModuleBase::Matrix3& R, std::string& sqpt);
+// Don't need to set_after_vc for Q_Vectors, since the q-points are only
+// used in phonon calculation and won't change after volume change.
+void print_klists(const Q_Vectors& qv, std::ofstream& ofs);
+#ifdef __MPI
+void kvec_mpi_k(Q_Vectors& qv);
+#endif // __MPI
+void kvec_ibz_kpoint(Q_Vectors& qv,
+                     const ModuleSymmetry::Symmetry& symm,
+                     bool use_symm,
+                     std::string& sqpt,
+                     const UnitCell& ucell,
+                     bool& match);
+
 } // namespace KVectorUtils
 
 #endif // K_VECTOR_UTILS_H