写了 project_to_mode 但还没调试好

CMakeLists.txt

@@ -14,7 +14,8 @@ find_package(Matplot++ REQUIRED)
 find_package(biu REQUIRED)
 find_package(Threads REQUIRED)
 
-add_executable(ufo src/fold.cpp src/unfold.cpp src/plot.cpp src/raman.cpp src/project-to-atom.cpp src/main.cpp)
+add_executable(ufo src/fold.cpp src/unfold.cpp src/project-to-mode.cpp
+  src/plot.cpp src/raman.cpp src/project-to-atom.cpp src/main.cpp)
 target_include_directories(ufo PRIVATE ${PROJECT_SOURCE_DIR}/include)
 target_link_libraries(ufo PRIVATE TBB::tbb Matplot++::matplot biu::biu)
 target_compile_features(ufo PRIVATE cxx_std_23)
@@ -37,3 +38,5 @@ add_test(NAME fold WORKING_DIRECTORY ${PROJECT_SOURCE_DIR}/test/fold COMMAND ufo
 add_test(NAME unfold WORKING_DIRECTORY ${PROJECT_SOURCE_DIR}/test/unfold COMMAND ufo unfold config.yaml)
 add_test(NAME zpp-to-yaml WORKING_DIRECTORY ${PROJECT_SOURCE_DIR}/test/zpp-to-yaml
   COMMAND sh -c "${CMAKE_BINARY_DIR}/ufo zpp-to-yaml < data.zpp > data.yaml")
+add_test(NAME project-to-mode WORKING_DIRECTORY ${PROJECT_SOURCE_DIR}/test/project-to-mode
+  COMMAND ufo project-to-mode config.yaml)

include/ufo.hpp

@@ -94,6 +94,7 @@ namespace ufo
   void fold(std::string config_file);
   void unfold(std::string config_file);
   void project_to_atom(std::string config_file);
+  void project_to_mode(std::string config_file);
   void plot_band(std::string config_file);
   void plot_point(std::string config_file);
   void raman_create_displacement(std::string config_file);

src/main.cpp

@@ -8,6 +8,7 @@ int main(int argc, const char** argv)
   else if (argv[1] == "zpp-to-yaml"s) std::cout << YAML::Node(biu::deserialize<ufo::CommonData>
     (biu::read<std::byte>(std::cin)));
   else if (argv[1] == "project-to-atom"s) ufo::project_to_atom(argv[2]);
+  else if (argv[1] == "project-to-mode"s) ufo::project_to_mode(argv[2]);
   else if (argv[1] == "raman-create-displacement"s) ufo::raman_create_displacement(argv[2]);
   else if (argv[1] == "raman-extract"s) ufo::raman_extract
   (

src/project-to-mode.cpp

@@ -0,0 +1,151 @@
+# include <ufo.hpp>
+
+void ufo::project_to_mode(std::string config_file)
+{
+  // 将选定的超胞中的某一个 q 点上的每一个模式到反折叠后的某个 q 点的投影，投射到单胞中的对应 q 点的模式上
+
+  struct Config
+  {
+    // 在单胞内取几个平面波的基矢
+    Eigen::Vector<std::size_t, 3> PrimativeCellBasisNumber;
+    // 要选定的超胞中的 q 点，sub qpoint，以及单胞中的 q 点
+    std::size_t SuperQpointIndex, SubQpointIndex, PrimativeQpointIndex;
+    // 输入输出文件的路径，直接输出 yaml，因为这些数据不再需要后处理
+    std::string InputDataFile, OutputFile;
+  };
+  struct Output
+  {
+    Eigen::Vector3d SuperQpoint, SubQpoint, PrimativeQpoint;
+    std::vector<std::vector<double>> Coefficient;
+  };
+
+  biu::Logger::Guard log(config_file);
+  auto config = YAML::LoadFile(config_file).as<Config>();
+  auto input = biu::deserialize<CommonData>
+    (biu::read<std::byte>(config.InputDataFile));
+  Output output
+  {
+    .SuperQpoint = input.Super.Qpoint[config.SuperQpointIndex].Qpoint,
+    .SubQpoint = input.Super.Qpoint[config.SuperQpointIndex].SubQpoint[config.SubQpointIndex],
+    .PrimativeQpoint = input.Primative.Qpoint[config.PrimativeQpointIndex].Qpoint,
+    .Coefficient = std::vector<std::vector<double>>(input.Super.Qpoint[config.SuperQpointIndex].Mode.size(),
+      std::vector<double>(input.Primative.Qpoint[config.PrimativeQpointIndex].Mode.size()))
+  };
+  if ((output.SubQpoint - output.PrimativeQpoint).norm() > 0.01)
+    log.error("sub qpoint {} != primative qpoint {}"_f(output.SubQpoint, output.PrimativeQpoint));
+
+  // 基在原胞中的表示
+  auto primative_basis = [&]
+  {
+    biu::Logger::Guard log;
+    std::vector<Eigen::VectorXcd> basis(config.PrimativeCellBasisNumber.prod());
+    for (auto [xyz_of_basis, i_of_basis]
+      : biu::sequence(config.PrimativeCellBasisNumber))
+    {
+      // 波矢就是 xyz_of_basis, 单位为单胞的倒格矢
+      // 将它的单位转换成埃^-1
+      auto wavevector = input.Primative.Cell.inverse() * xyz_of_basis.cast<double>();
+      log.debug("xyz_of_basis: {}"_f(xyz_of_basis));
+      log.debug("wavevector: {:.3f}"_f(fmt::join(wavevector, ", ")));
+      // 将原子坐标单位也转换为埃
+      auto atom_position = input.Primative.AtomPosition * input.Primative.Cell;
+      // 计算基矢
+      basis[i_of_basis] = (2i * std::numbers::pi_v<double> * (atom_position * wavevector)).array().exp();
+    }
+    return basis;
+  }();
+  // 基在超胞中的表示
+  auto super_basis = [&]
+  {
+    biu::Logger::Guard log;
+    std::vector<Eigen::VectorXcd> basis(config.PrimativeCellBasisNumber.prod());
+    auto diff_of_sub_qpoint_by_reciprocal_modified_super_cell = [&]
+    {
+      for
+      (
+        auto [diff_of_sub_qpoint_by_reciprocal_modified_super_cell, i_of_sub_qpoint]
+          : biu::sequence(input.Super.CellMultiplier)
+      )
+        if (i_of_sub_qpoint == config.SubQpointIndex) return diff_of_sub_qpoint_by_reciprocal_modified_super_cell;
+      std::unreachable();
+    }();
+    log.debug("diff_of_sub_qpoint_by_reciprocal_modified_super_cell: {}"_f
+      (diff_of_sub_qpoint_by_reciprocal_modified_super_cell));
+    for (auto [xyz_of_basis, i_of_basis]
+      : biu::sequence(config.PrimativeCellBasisNumber))
+    {
+      // 计算波矢, 单位为单胞的倒格矢
+      auto wavevector_by_reciprocal_primative_cell = xyz_of_basis.cast<double>()
+        + input.Super.CellMultiplier.cast<double>().cwiseInverse().asDiagonal()
+        * diff_of_sub_qpoint_by_reciprocal_modified_super_cell.cast<double>();
+      // 将单位转换为埃^-1
+      auto wavevector = input.Primative.Cell.inverse() * wavevector_by_reciprocal_primative_cell;
+      // 将原子坐标单位也转换为埃
+      auto atom_position = input.Super.AtomPosition
+        * (input.Super.CellDeformation * input.Super.CellMultiplier.cast<double>().asDiagonal() * input.Primative.Cell);
+      // 计算基矢
+      basis[i_of_basis] = (2i * std::numbers::pi_v<double> * (atom_position * wavevector)).array().exp();
+    }
+    return basis;
+  }();
+  // 将所有模式投影到基矢上
+  auto primative_projection_coefficient =
+    input.Primative.Qpoint[config.PrimativeQpointIndex].Mode
+      | ranges::views::transform([&](const auto& mode)
+        {
+          return primative_basis
+            | ranges::views::transform([&](const auto& basis)
+              { return (basis.transpose().conjugate() * mode.EigenVector).eval() | biu::fromEigen; })
+            | ranges::to_vector
+            | biu::toEigen<>;
+        })
+      | ranges::to_vector;
+  for (auto i_of_primative_mode : std::views::iota(0u, primative_projection_coefficient.size()))
+    for (auto i_of_direction : std::array{0, 1, 2})
+    {
+      if (i_of_primative_mode != 0 && i_of_primative_mode != 11) continue;
+      log.debug("primitive mode: {} {} {:.3f}"_f
+        (i_of_primative_mode, i_of_direction, fmt::join(primative_projection_coefficient[i_of_primative_mode].col(i_of_direction).array().abs(), ", ")));
+      log.debug("primitive mode: {} {} {:.3f}"_f
+        (i_of_primative_mode, i_of_direction, fmt::join(primative_projection_coefficient[i_of_primative_mode].col(i_of_direction).array().arg(), ", ")));
+    }
+  log.debug("primitive 0 {}"_f(primative_projection_coefficient[0].col(2).transpose().conjugate() * primative_projection_coefficient[0].col(2)));
+  log.debug("primitive 11 {}"_f(primative_projection_coefficient[11].col(2).transpose().conjugate() * primative_projection_coefficient[11].col(2)));
+  log.debug("0 x 11 {}"_f(primative_projection_coefficient[11].col(2).transpose().conjugate() * primative_projection_coefficient[0].col(2)));
+  auto super_projection_coefficient =
+    input.Super.Qpoint[config.SuperQpointIndex].Mode
+      | ranges::views::transform([&](const auto& mode)
+        {
+          return super_basis
+            | ranges::views::transform([&](const auto& basis)
+              { return (basis.transpose().conjugate() * mode.EigenVector).eval() | biu::fromEigen; })
+            | ranges::to_vector
+            | biu::toEigen<>;
+        })
+      | ranges::to_vector;
+  log.debug("super 2 x primitive 0 {}"_f(super_projection_coefficient[2].col(2).transpose().conjugate() * primative_projection_coefficient[0].col(2)));
+  log.debug("super 2 x primitive 11 {}"_f(super_projection_coefficient[2].col(2).transpose().conjugate() * primative_projection_coefficient[11].col(2)));
+  log.debug("super 2 x super 2 {}"_f(super_projection_coefficient[2].col(2).transpose().conjugate() * super_projection_coefficient[2].col(2)));
+  for (auto i_of_direction : std::array{0, 1, 2})
+  {
+    log.debug("super mode: 2 {} {:.3f}"_f
+      (i_of_direction, fmt::join(super_projection_coefficient[2].col(i_of_direction).array().abs(), ", ")));
+    log.debug("super mode: 2 {} {:.3f}"_f
+      (i_of_direction, fmt::join(super_projection_coefficient[2].col(i_of_direction).array().arg(), ", ")));
+  }
+  // 将它们互相投影，得到结果
+  for (auto [i_of_super_mode, super_mode] : ranges::views::enumerate(super_projection_coefficient))
+  {
+    for (auto [i_of_primative_mode, primative_mode]
+      : ranges::views::enumerate(primative_projection_coefficient))
+      for (auto i_of_direction : std::array{0, 1, 2})
+        output.Coefficient[i_of_super_mode][i_of_primative_mode] +=
+          (super_mode.col(i_of_direction).transpose().conjugate() * primative_mode.col(i_of_direction)).norm();
+    auto sum = ranges::accumulate(output.Coefficient[i_of_super_mode], 0.);
+    for (auto& coefficient : output.Coefficient[i_of_super_mode])
+      coefficient *=
+        input.Super.Qpoint[config.SuperQpointIndex].Mode[i_of_super_mode].WeightOnUnfold[config.SubQpointIndex] / sum;
+  }
+  // 输出
+  std::ofstream(config.OutputFile) << YAML::Node(output);
+}

src/unfold.cpp

@@ -7,9 +7,8 @@
 void ufo::unfold(std::string config_file)
 {
   // 反折叠的原理: 将超胞中的原子运动状态, 投影到一组平面波构成的基矢中.
-  // 每一个平面波的波矢由两部分相加得到: 一部分是单胞倒格子的整数倍, 所取的个数有一定任意性, 论文中建议取大约单胞中原子个数那么多个;
-  //  对于没有缺陷的情况, 取一个应该就足够了.
-  // 这些平面波以原胞为周期。
+  // 每一个平面波的波矢由两部分相加得到: 一部分是单胞倒格子的整数倍, 所取的个数有一定任意性, 论文中建议取大约单胞中原子个数那么多个，
+  // 但我觉得可以多取一些；这些平面波以原胞为周期，它用于尽可能描述
   // 另一部分是超胞倒格子的整数倍, 取 n 个, n 为超胞对应的单胞的倍数, 其实也就是倒空间中单胞对应倒格子中超胞的格点.
   // 只要第一部分取得足够多, 那么单胞中原子的状态就可以完全被这些平面波描述.
   // 将超胞中原子的运动状态投影到这些基矢上, 计算出投影的系数, 就可以将超胞的原子运动状态分解到单胞中的多个 q 点上.
@@ -143,16 +142,15 @@ void ufo::unfold(std::string config_file)
       for (auto [xyz_of_basis, i_of_basis]
         : biu::sequence(primative_cell_basis_number))
       {
-        // 计算 q 点的坐标, 单位为单胞的倒格矢
-        auto diff_of_sub_qpoint_by_reciprocal_primative_cell = xyz_of_basis.cast<double>()
+        // 计算波矢, 单位为单胞的倒格矢
+        auto wavevector_by_reciprocal_primative_cell = xyz_of_basis.cast<double>()
           + super_cell_multiplier.cast<double>().cwiseInverse().asDiagonal()
           * diff_of_sub_qpoint_by_reciprocal_modified_super_cell.cast<double>();
         // 将单位转换为埃^-1
-        auto diff_of_sub_qpoint = (diff_of_sub_qpoint_by_reciprocal_primative_cell.transpose()
-          * (primative_cell.transpose().inverse())).transpose();
+        auto wavevector = primative_cell.inverse() * wavevector_by_reciprocal_primative_cell;
         // 计算基矢
         basis[i_of_sub_qpoint][i_of_basis]
-          = (2i * std::numbers::pi_v<double> * (atom_position * diff_of_sub_qpoint)).array().exp();
+          = (2i * std::numbers::pi_v<double> * (atom_position * wavevector)).array().exp();
       }
     }
     return basis;
@@ -204,6 +202,46 @@ void ufo::unfold(std::string config_file)
     (config.SuperPhonopy, config.SuperQpoint, output.Super));
   output.Super.CellDeformation = config.SuperCellDeformation;
   output.Super.CellMultiplier = config.SuperCellMultiplier;
+  // 填充 SubQpoint
+  for (auto i_of_super_qpoint : std::views::iota(0u, output.Super.Qpoint.size()))
+    for
+    (
+      auto [diff_of_sub_qpoint_by_reciprocal_modified_super_cell, i_of_sub_qpoint]
+        : biu::sequence(config.SuperCellMultiplier)
+    )
+    {
+      /*
+        SubQpointByReciprocalModifiedSuperCell = XyzOfDiffOfSubQpointByReciprocalModifiedSuperCell +
+          MetaQpointByReciprocalModifiedSuperCell;
+        SubQpoint = SubQpointByReciprocalModifiedSuperCell.transpose() * ReciprocalModifiedSuperCell;
+        SubQpoint = SubQpointByReciprocalPrimativeCell.transpose() * ReciprocalPrimativeCell;
+        ReciprocalModifiedSuperCell = ModifiedSuperCell.inverse().transpose();
+        ReciprocalPrimativeCell = PrimativeCell.inverse().transpose();
+        ModifiedSuperCell = SuperCellMultiplier.asDiagonal() * PrimativeCell;
+        MetaQpoint = MetaQpointByReciprocalModifiedSuperCell.transpose() * ReciprocalModifiedSuperCell;
+        MetaQpoint = MetaQpointByReciprocalSuperCell.transpose() * ReciprocalSuperCell;
+        ReciprocalSuperCell = SuperCell.inverse().transpose();
+        ModifiedSuperCell = SuperCellDeformation * SuperCell;
+        SuperCell = SuperCellMultiplier.asDiagonal() * PrimativeCell;
+        整理可以得到:
+        SubQpointByReciprocalPrimativeCell = SuperCellMultiplier.asDiagonal().inverse() *
+          (XyzOfDiffOfSubQpointByReciprocalModifiedSuperCell +
+            SuperCellDeformation.inverse() * MetaQpointByReciprocalSuperCell);
+        但注意到, 这样得到的 SubQpoint 可能不在 ReciprocalPrimativeCell 中
+          (当 SuperCellDeformation 不是单位矩阵时, 边界附近的一两条 SubQpoint 会出现这种情况).
+        解决办法是, 在赋值时, 仅取 SubQpointByReciprocalPrimativeCell 的小数部分.
+      */
+      auto sub_qpoint_by_reciprocal_primative_cell =
+      (
+        config.SuperCellMultiplier.cast<double>().cwiseInverse().asDiagonal()
+        * (
+          diff_of_sub_qpoint_by_reciprocal_modified_super_cell.cast<double>()
+            + config.SuperCellDeformation.inverse() * output.Super.Qpoint[i_of_super_qpoint].Qpoint
+        )
+      ).eval();
+      output.Super.Qpoint[i_of_super_qpoint].SubQpoint.push_back(sub_qpoint_by_reciprocal_primative_cell.array()
+        - sub_qpoint_by_reciprocal_primative_cell.array().floor());
+    }
   log.info("Done.");
 
   log.info("Constructing basis...");

test/project-to-mode/.gitignore

@@ -0,0 +1 @@
+output.yaml

test/project-to-mode/config.yaml

@@ -0,0 +1,6 @@
+PrimativeCellBasisNumber: [ 16, 16, 16 ]
+SuperQpointIndex: 0
+SubQpointIndex: 0
+PrimativeQpointIndex: 0
+InputDataFile: data.zpp
+OutputFile: output.yaml