add fold

man 增加输出

fold.cpp

@@ -0,0 +1,98 @@
+# include <iostream>
+# include <array>
+# include <numbers>
+# include <numeric>
+# include <fstream>
+# include <optional>
+# include <array>
+# include <utility>
+# include <yaml-cpp/yaml.h>
+# include <eigen3/Eigen/Dense>
+# include <concurrencpp/concurrencpp.h>
+# include <fmt/format.h>
+
+using namespace std::literals;
+
+// 计算单胞中的 q 点在超胞中的对应
+
+struct Input
+{
+  // 单胞的三个格矢，每行表示一个格矢的坐标，单位为埃
+  Eigen::Matrix3d PrimativeCell;
+  // 单胞到超胞的格矢转换时用到的矩阵
+  Eigen::Vector<unsigned, 3> SuperCellMultiplier;
+  Eigen::Matrix<double, 3, 3> SuperCellDeformation;
+
+  // Q 点的坐标，单位为单胞的倒格矢
+  std::vector<Eigen::Vector3d> QPointData;
+};
+template<> struct YAML::convert<Input> { static bool decode(const Node& node, Input& input); };
+
+struct Output
+{
+  // Q 点的坐标，单位为超胞的倒格矢
+  std::vector<Eigen::Vector3d> QPointData;
+};
+
+int main(int argc, char** argv)
+{
+  if (argc != 3)
+    throw std::runtime_error("Usage: " + std::string(argv[0]) + " input.yaml output.yaml");
+
+  auto input = YAML::LoadFile(argv[1]).as<Input>();
+  Output output;
+
+  for (auto qpoint_by_reciprocal_primative_cell : input.QPointData)
+  {
+    // 计算出这个 q 点的绝对坐标, 再计算出它相对于超胞倒格子的相对坐标. 将这个结果取小数部分, 就得到了 meta qpoint 的坐标(相对于超胞倒格子)
+    std::cout << "PrimativeCell:\n" << input.PrimativeCell << "\n";
+    auto reciprocal_primative_cell = input.PrimativeCell.inverse().transpose();
+    std::cout << "reciprocal_primative_cell:\n" << reciprocal_primative_cell << "\n";
+    auto qpoint = (qpoint_by_reciprocal_primative_cell.transpose() * reciprocal_primative_cell).transpose();
+    std::cout << "qpoint:\n" << qpoint << "\n";
+    auto reciprocal_super_cell = (input.SuperCellDeformation * input.SuperCellMultiplier.cast<double>().asDiagonal()
+      * input.PrimativeCell).inverse().transpose();
+    std::cout << "reciprocal_super_cell:\n" << reciprocal_super_cell << "\n";
+    auto qpoint_by_reciprocal_super_cell = (qpoint.transpose() * reciprocal_super_cell.inverse()).transpose();
+    std::cout << "qpoint_by_reciprocal_super_cell:\n" << qpoint_by_reciprocal_super_cell << "\n";
+    auto meta_qpoint_by_reciprocal_super_cell = [&]
+      { auto _ = qpoint_by_reciprocal_super_cell.array(); return _ - _.floor(); }();
+    std::cout << "meta_qpoint_by_reciprocal_super_cell:\n" << meta_qpoint_by_reciprocal_super_cell << "\n";
+    output.QPointData.push_back(meta_qpoint_by_reciprocal_super_cell);
+  }
+
+  std::ofstream(argv[2]) << [&]
+  {
+    std::stringstream print;
+    print << "QPointData:\n";
+    for (auto& qpoint: output.QPointData)
+      print << fmt::format("  - [ {:.3f}, {:.3f}, {:.3f} ]\n", qpoint(0), qpoint(1), qpoint(2));
+    return print.str();
+  }();
+}
+
+bool YAML::convert<Input>::decode(const Node& node, Input& input)
+{
+  for (unsigned i = 0; i < 3; i++)
+    for (unsigned j = 0; j < 3; j++)
+      input.PrimativeCell(i, j) = node["lattice"][i][j].as<double>();
+
+  input.SuperCellMultiplier.setZero();
+  for (unsigned i = 0; i < 3; i++)
+    input.SuperCellMultiplier(i) = node["SuperCellMultiplier"][i].as<int>();
+
+  for (unsigned i = 0; i < 3; i++)
+    for (unsigned j = 0; j < 3; j++)
+      input.SuperCellDeformation(i, j) = node["SuperCellDeformation"][i][j].as<double>();
+
+  auto points = node["points"].as<std::vector<std::vector<double>>>();
+
+  input.QPointData.resize(points.size());
+  for (unsigned i = 0; i < points.size(); i++)
+  {
+    for (unsigned j = 0; j < 3; j++)
+      input.QPointData[i](j) = points[i][j];
+  }
+
+  return true;
+}

main.cpp

@@ -10,6 +10,7 @@
 # include <eigen3/Eigen/Dense>
 # include <concurrencpp/concurrencpp.h>
 # include <fmt/format.h>
+# include <highfive/H5Easy.hpp>
 
 using namespace std::literals;
 
@@ -70,7 +71,7 @@ struct Output
   // 关于各个 Q 点的数据
   struct QPointDataType_
   {
-    // Q 点的坐标，单位为超胞的倒格矢
+    // Q 点的坐标，单位为单胞的倒格矢
     Eigen::Vector3d QPoint;
 
     // 来源于哪个 Q 点, 单位为超胞的倒格矢
@@ -110,7 +111,9 @@ int main(int argc, const char** argv)
   if (argc != 3)
     throw std::runtime_error("Usage: " + std::string(argv[0]) + " input.yaml output.yaml");
 
+  std::cerr << "Reading input file..." << std::flush;
   auto input = YAML::LoadFile(argv[1]).as<Input>();
+  std::cerr << "Done." << std::endl;
 
   // 反折叠的原理: 将超胞中的原子运动状态, 投影到一组平面波构成的基矢中.
   // 每一个平面波的波矢由两部分相加得到: 一部分是单胞倒格子的整数倍, 所取的个数有一定任意性, 论文中建议取大约单胞中原子个数那么多个;
@@ -125,6 +128,7 @@ int main(int argc, const char** argv)
   // 由于基只与这个相对位置有关（也就是说，不同 q 点的基是一样的），因此可以先计算出所有的基，这样降低计算量。
   // 外层下标对应超胞倒格子的整数倍那部分(第二部分), 也就是不同的 sub qpoint
   // 内层下标对应单胞倒格子的整数倍那部分(第一部分), 也就是 sub qpoint 上的不同平面波（取的数量越多，结果越精确）
+  std::cerr << "Calculating basis..." << std::flush;
   std::vector<std::vector<Eigen::VectorXcd>> basis(input.SuperCellMultiplier.prod());
   // 每个 q 点对应的一组 sub qpoint。不同的 q 点所对应的 sub qpoint 是不一样的，但 sub qpoint 与 q 点的相对位置一致。
   // 这里 xyz_of_diff_of_sub_qpoint 即表示这个相对位置，单位为超胞的倒格矢
@@ -146,12 +150,15 @@ int main(int argc, const char** argv)
         = (2i * std::numbers::pi_v<double> * (input.AtomPosition * qpoint)).array().exp();
     }
   }
+  std::cerr << "Done." << std::endl;
 
   // 计算投影的结果
   // 最外层下标对应反折叠前的 q 点, 第二层下标对应不同模式, 第三层下标对应这个模式在反折叠后的 q 点(sub qpoint)
   std::vector<std::vector<std::vector<double>>> projection_coefficient(input.QPointData.size());
   for (unsigned i_of_qpoint = 0; i_of_qpoint < input.QPointData.size(); i_of_qpoint++)
   {
+    std::cerr << fmt::format("\rCalculating projection coefficient for qpoint {}/{}...",
+      i_of_qpoint, input.QPointData.size()) << std::flush;
     projection_coefficient[i_of_qpoint].resize(input.QPointData[i_of_qpoint].ModeData.size());
     for (unsigned i_of_mode = 0; i_of_mode < input.QPointData[i_of_qpoint].ModeData.size(); i_of_mode++)
     {
@@ -173,8 +180,10 @@ int main(int argc, const char** argv)
         __ /= sum;
     }
   }
+  std::cerr << "Done." << std::endl;
 
   // 填充输出对象
+  std::cerr << "Filling output object..." << std::flush;
   Output output;
   for (unsigned i_of_qpoint = 0; i_of_qpoint < input.QPointData.size(); i_of_qpoint++)
     for (auto [xyz_of_diff_of_sub_qpoint_by_reciprocal_super_cell, i_of_sub_qpoint]
@@ -235,9 +244,11 @@ int main(int argc, const char** argv)
           __.Weight = projection_coefficient[i_of_qpoint][i_of_mode][i_of_sub_qpoint];
         }
     }
+  std::cerr << "Done." << std::endl;
 
   // std::ofstream(argv[2]) << YAML::Node(output);
   // YAML 输出得太丑了，我来自己写
+  std::cerr << "Writing output file..." << std::flush;
   std::ofstream(argv[2]) << [&]
   {
     std::stringstream print;
@@ -256,6 +267,7 @@ int main(int argc, const char** argv)
     }
     return print.str();
   }();
+  std::cerr << "Done." << std::endl;
 }
 
 // 从文件中读取输入, 文件中应当包含: (大多数据可以直接从 phonopy 的输出中复制)