重构

CMakeLists.txt

@@ -14,14 +14,14 @@ 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-create-displacement.cpp src/raman-extract.cpp
+add_executable(ufo src/fold.cpp src/unfold.cpp src/plot.cpp src/raman.cpp src/project-to-atom.cpp src/main.cpp)
-  src/raman-apply-contribution.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)
 target_compile_options(ufo PRIVATE -fexperimental-library)
 if (CMAKE_BUILD_TYPE STREQUAL "Debug")
   target_compile_definitions(ufo PRIVATE BIU_LOGGER_DEBUG BIU_LOGGER_SOURCE_ROOT="${CMAKE_SOURCE_DIR}")
+  target_compile_options(ufo PRIVATE -ftemplate-backtrace-limit=0)
 endif()
 
 install(TARGETS ufo RUNTIME DESTINATION ${CMAKE_INSTALL_BINDIR})

include/ufo.hpp

@@ -11,82 +11,92 @@ namespace ufo
   using namespace biu::literals;
   using namespace biu::stream_operators;
 
+  // 许多函数通用的一个结构体
+  struct CommonData
+  {
+    // 原胞相关信息，其中存储的信息是直接读入的，而不是反折叠后的结果
+    struct
+    {
+      // 晶胞，单位为埃，每行代表一个格矢
+      Eigen::Matrix3d Cell;
+      // 原子种类及其个数
+      std::vector<std::pair<std::string, std::size_t>> AtomType;
+      // 原子位置，单位为原胞格矢
+      Eigen::MatrixX3d AtomPosition;
+      // 单胞中每个 q 点的信息
+      struct QpointType
+      {
+        // q 点的坐标，单位为倒格矢
+        Eigen::Vector3d Qpoint;
+        struct ModeType
+        {
+          // 振动频率，单位为 THz
+          double Frequency;
+          // 振动模式，单位为埃
+          Eigen::MatrixX3cd EigenVector;
+          // 拉曼张量
+          // TODO: 单位是什么？
+          std::optional<double> WeightOnRaman;
+          // work around for fpr
+          std::optional<std::array<std::array<double, 3>, 3>> RamanTensor;
+          using serialize = zpp::bits::members<4>;
+        };
+        // 各个模式的的信息
+        std::vector<ModeType> Mode;
+        using serialize = zpp::bits::members<2>;
+      };
+      std::vector<QpointType> Qpoint;
+      using serialize = zpp::bits::members<4>;
+    } Primative;
+
+    // 超胞相关信息，包括反折叠后的信息
+    struct
+    {
+      // 超胞的晶胞，单位为埃，每行代表一个格矢
+      Eigen::Matrix3d Cell;
+      // 单胞与到超胞的变换
+      Eigen::Vector3i CellMultiplier;
+      Eigen::Matrix3d CellDeformation;
+      std::vector<std::pair<std::string, std::size_t>> AtomType;
+      Eigen::MatrixX3d AtomPosition;
+      struct QpointType
+      {
+        Eigen::Vector3d Qpoint;
+        std::vector<Eigen::Vector3d> SubQpoint;
+        struct ModeType
+        {
+          double Frequency;
+          Eigen::MatrixX3cd EigenVector;
+          std::vector<double> WeightOnUnfold;
+          std::optional<double> WeightOnSelectedAtom;
+          std::optional<double> WeightOnRaman;
+          std::optional<std::array<std::array<double, 3>, 3>> RamanTensor;
+          using serialize = zpp::bits::members<6>;
+        };
+        std::vector<ModeType> Mode;
+        using serialize = zpp::bits::members<3>;
+      };
+      std::vector<QpointType> Qpoint;
+      using serialize = zpp::bits::members<6>;
+    } Super;
+
+    // 选择的原子
+    std::optional<std::set<std::size_t>> SelectedAtom;
+
+    // 拉曼散射选择的偏振方向
+    std::optional<std::array<Eigen::Vector3d, 2>> RamanPolarization;
+
+    // 各种原子的质量
+    std::map<std::string, double> AtomMass;
+    using serialize = zpp::bits::members<5>;
+  };
+
   void fold(std::string config_file);
   void unfold(std::string config_file);
+  void project_to_atom(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);
   void raman_extract(std::vector<std::string> files);
   void raman_apply_contribution(std::string config_file);
-
-  // 许多函数都需要用到这个，所以写到头文件中
-  struct UnfoldOutput
-  {
-    Eigen::Matrix3d PrimativeCell;
-    Eigen::Matrix3i SuperCellTransformation;
-    Eigen::Vector3i SuperCellMultiplier;
-    Eigen::Matrix3d SuperCellDeformation;
-    std::optional<std::vector<std::size_t>> SelectedAtoms;
-
-    // 关于各个 Q 点的数据
-    struct QpointDataType
-    {
-      // Q 点的坐标，单位为单胞的倒格矢
-      Eigen::Vector3d Qpoint;
-
-      // 来源于哪个 Q 点, 单位为超胞的倒格矢
-      Eigen::Vector3d Source;
-      std::size_t SourceIndex;
-
-      // 关于这个 Q 点上各个模式的数据
-      struct ModeDataType
-      {
-        // 模式的频率，单位为 THz
-        double Frequency;
-        // 模式的权重
-        double Weight;
-      };
-      std::vector<ModeDataType> ModeData;
-    };
-    std::vector<QpointDataType> QpointData;
-
-    struct MetaQpointDataType
-    {
-      // Q 点的坐标，单位为单胞的倒格矢
-      Eigen::Vector3d Qpoint;
-
-      // 关于这个 Q 点上各个模式的数据
-      struct ModeDataType
-      {
-        // 模式的频率，单位为 THz
-        double Frequency;
-        // 模式中各个原子的运动状态
-        // 这个数据应当是这样得到的：动态矩阵的 eigenvector 乘以 $\exp(-2 \pi i \vec q \cdot \vec r)$
-        // 这个数据可以认为是原子位移中, 关于超胞有周期性的那一部分, 再乘以原子质量的开方.
-        // 这个数据会在 unfold 时被归一化
-        Eigen::MatrixX3cd AtomMovement;
-      };
-      std::vector<ModeDataType> ModeData;
-    };
-    std::vector<MetaQpointDataType> MetaQpointData;
-
-    using serialize = zpp::bits::members<7>;
-  };
-  struct DisplacementOutput
-  {
-    struct ModeData_t
-    {
-      std::size_t MetaQpointIndex;
-      std::size_t ModeIndex;
-      // 每个原子的位移，单位为埃
-      Eigen::MatrixX3d AtomMovement;
-      // 为了使得位移最大的原子的位移恰好是 input.MaxDisplacement, 在位移上乘以了多大的系数
-      double Ratio;
-    };
-    std::vector<ModeData_t> ModeData;
-    Eigen::Vector3d AtomMasses;
-    double MaxDisplacement;
-    std::vector<std::size_t> QpointIndices;
-    using serialize = zpp::bits::members<4>;
-  };
 }

src/fold.cpp

@@ -2,17 +2,13 @@
 
 void ufo::fold(std::string config_file)
 {
-  struct Input
+  struct Config
   {
     Eigen::Matrix3d SuperCellDeformation;
     Eigen::Vector3i SuperCellMultiplier;
     std::vector<Eigen::Vector3d> Qpoints;
-    std::optional<std::string> OutputFile;
-  };
-  struct Output
-  {
-    std::vector<Eigen::Vector3d> Qpoints;
   };
+
   auto fold = []
   (
     Eigen::Vector3d qpoint_in_reciprocal_primitive_cell_by_reciprocal_primitive_cell,
@@ -44,15 +40,13 @@ void ufo::fold(std::string config_file)
     */
     return (qpoint_by_reciprocal_super_cell.array() - qpoint_by_reciprocal_super_cell.array().floor()).matrix();
   };
-  auto input = YAML::LoadFile(config_file).as<Input>();
+
-  Output output;
+  biu::Logger::Guard log(config_file);
-  output.Qpoints = input.Qpoints
+  auto input = YAML::LoadFile(config_file).as<Config>();
-    | ranges::views::transform([&](auto& qpoint)
+  for (const auto& qpoint : input.Qpoints) log.info("{} -> {}"_f
-    {
+    (
-      return fold(qpoint, input.SuperCellDeformation * input.SuperCellMultiplier.cast<double>().asDiagonal());
+      qpoint,
-    })
+      fold(qpoint,
-    | ranges::to_vector;
+        input.SuperCellDeformation * input.SuperCellMultiplier.cast<double>().asDiagonal())
-  
+    ));
-  // 默认的输出太丑了，但是不想手动写了，忍一下
-  std::ofstream(input.OutputFile.value_or("output.yaml")) << YAML::Node(output);
 }

src/main.cpp

@@ -5,6 +5,9 @@ int main(int argc, const char** argv)
   using namespace biu::literals;
   if (argv[1] == "fold"s) ufo::fold(argv[2]);
   else if (argv[1] == "unfold"s) ufo::unfold(argv[2]);
+  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] == "raman-create-displacement"s) ufo::raman_create_displacement(argv[2]);
   else if (argv[1] == "raman-extract"s) ufo::raman_extract
   (

src/plot.cpp

@@ -4,9 +4,9 @@
 
 void ufo::plot_band(std::string config_file)
 {
-  struct Input
+  struct Config
   {
-    std::string UnfoldedDataFile;
+    std::string InputDataFile;
     // 要画图的 q 点路径列表
     // 内层表示一个路径上的 q 点，外层表示不同的路径
     // 单位为倒格矢
@@ -29,26 +29,26 @@ void ufo::plot_band(std::string config_file)
     std::optional<std::string> OutputDataFile;
   };
 
-  // 根据 q 点路径, 搜索要使用的 q 点，返回的是 q 点在 QpointData 中的索引以及到路径起点的距离，以及这段路径的总长度
+  // 根据 q 点路径, 搜索要使用的 q 点，返回的是超胞中 Qpint 的索引和 SubQpoint 的索引，以及到路径起点的距离，以及这段路径的总长度
   auto search_qpoints = []
   (
     const Eigen::Matrix3d& primative_cell,
     const std::pair<Eigen::Vector3d, Eigen::Vector3d>& path,
-    const std::vector<Eigen::Vector3d>& qpoints,
+    const std::vector<std::vector<Eigen::Vector3d>>& qpoints,
     double threshold, bool exclude_endpoint = false
   )
   {
     // 对于 output 中的每一个点, 检查这个点是否在路径上. 如果在, 把它加入到 selected_qpoints 中
     // 键为这个点到起点的距离
-    boost::container::flat_map<double, std::size_t> selected_qpoints;
+    boost::container::flat_map<double, std::pair<std::size_t, std::size_t>> selected_qpoints;
     auto begin = (path.first.transpose() * primative_cell.reverse()).transpose().eval();
     auto end = (path.second.transpose() * primative_cell.reverse()).transpose().eval();
-    for (std::size_t i = 0; i < qpoints.size(); i++)
+    for (std::size_t i = 0; i < qpoints.size(); i++) for (std::size_t j = 0; j < qpoints[i].size(); j++)
       for (auto cell_shift
         : biu::sequence(Eigen::Vector3i(-1, -1, -1), Eigen::Vector3i(2, 2, 2)))
       {
-        auto qpoint
+        Eigen::Vector3d qpoint
-          = ((qpoints[i] + cell_shift.first.cast<double>()).transpose() * primative_cell.reverse()).transpose().eval();
+          = primative_cell.reverse().transpose() * (qpoints[i][j] + cell_shift.first.cast<double>());
         // 计算这个点到前两个点所在直线的距离
         auto distance = (end - begin).cross(qpoint - begin).norm()
           / (path.second - path.first).norm();
@@ -61,7 +61,8 @@ void ufo::plot_band(std::string config_file)
           auto distance3 = (end - begin).norm();
           if (distance1 < distance3 + threshold && distance2 < distance3 + threshold)
             // 如果这个点不在终点处, 或者不排除终点, 则加入
-            if (distance2 > threshold || !exclude_endpoint) selected_qpoints.emplace(distance1, i);
+            if (distance2 > threshold || !exclude_endpoint)
+              selected_qpoints.emplace(distance1, std::pair{i, j});
         }
       }
     // 去除非常接近的点
@@ -80,11 +81,13 @@ void ufo::plot_band(std::string config_file)
   auto calculate_values = []
   (
     // search_qpoints 的第一个返回值
-    const boost::container::flat_map<double, std::size_t>& path,
+    const boost::container::flat_map<double, std::pair<std::size_t, std::size_t>>& path,
     // 每一条连续路径的第一个 q 点的索引
     const std::set<std::size_t>& path_begin,
-    // 所有 q 点的数据（需要用到它的频率和权重）
+    // 各个模式的频率和权重，几个维度分别为：MetaQpointIndex, ModeIndex
-    const std::vector<UnfoldOutput::QpointDataType>& qpoints,
+    const std::vector<std::vector<double>>& frequency,
+    // 各个模式的权重，几个维度分别为：MetaQpointIndex, ModeIndex, SubQpointIndex
+    const std::vector<std::vector<std::vector<double>>>& weight,
     // 用于插值的分辨率和范围
     const std::array<std::size_t, 2>& resolution,
     const std::array<double, 2>& frequency_range,
@@ -96,7 +99,7 @@ void ufo::plot_band(std::string config_file)
     auto blend = [&]
     (
       // 两个点的索引
-      std::size_t a, std::size_t b,
+      std::pair<std::size_t, std::size_t> a, std::pair<std::size_t, std::size_t> b,
       // 按照连续路径混合还是按照断开的路径混合
       bool continuous,
       // 第一个点占的比例
@@ -105,38 +108,36 @@ void ufo::plot_band(std::string config_file)
     ) -> std::vector<double>
     {
       // 混合得到的频率和权重
-      std::vector<double> frequency, weight;
+      std::vector<double> frequency_result, weight_result;
       // 如果是连续路径，将每个模式的频率和权重按照比例混合
       if (continuous)
       {
-        assert(qpoints[a].ModeData.size() == qpoints[b].ModeData.size());
+        for (std::size_t i = 0; i < frequency[a.first].size(); i++)
-        for (std::size_t i = 0; i < qpoints[a].ModeData.size(); i++)
         {
-          frequency.push_back
+          frequency_result.push_back(frequency[a.first][i] * ratio + frequency[b.first][i] * (1 - ratio));
-            (qpoints[a].ModeData[i].Frequency * ratio + qpoints[b].ModeData[i].Frequency * (1 - ratio));
+          weight_result.push_back(weight[a.first][i][a.second] * ratio + weight[b.first][i][b.second] * (1 - ratio));
-          weight.push_back(qpoints[a].ModeData[i].Weight * ratio + qpoints[b].ModeData[i].Weight * (1 - ratio));
         }
       }
       // 如果是不连续路径，将每个模式的权重乘以比例，最后相加
       else
       {
-        for (std::size_t i = 0; i < qpoints[a].ModeData.size(); i++)
+        for (std::size_t i = 0; i < frequency[a.first].size(); i++)
         {
-          frequency.push_back(qpoints[a].ModeData[i].Frequency);
+          frequency_result.push_back(frequency[a.first][i]);
-          weight.push_back(qpoints[a].ModeData[i].Weight * ratio);
+          weight_result.push_back(weight[a.first][i][a.second] * ratio);
         }
-        for (std::size_t i = 0; i < qpoints[b].ModeData.size(); i++)
+        for (std::size_t i = 0; i < frequency[b.first].size(); i++)
         {
-          frequency.push_back(qpoints[b].ModeData[i].Frequency);
+          frequency_result.push_back(frequency[b.first][i]);
-          weight.push_back(qpoints[b].ModeData[i].Weight * (1 - ratio));
+          weight_result.push_back(weight[b.first][i][b.second] * (1 - ratio));
         }
       }
       std::vector<double> result(resolution);
-      for (std::size_t i = 0; i < frequency.size(); i++)
+      for (std::size_t i = 0; i < frequency_result.size(); i++)
       {
-        std::ptrdiff_t index = (frequency[i] - frequency_range[0]) / (frequency_range[1] - frequency_range[0])
+        std::ptrdiff_t index = (frequency_result[i] - frequency_range[0]) / (frequency_range[1] - frequency_range[0])
           * resolution;
-        if (index >= 0 && index < static_cast<std::ptrdiff_t>(resolution)) result[index] += weight[i];
+        if (index >= 0 && index < static_cast<std::ptrdiff_t>(resolution)) result[index] += weight_result[i];
       }
       return result;
     };
@@ -155,7 +156,8 @@ void ufo::plot_band(std::string config_file)
       ));
       else values.push_back(blend
       (
-        std::prev(it)->second, it->second, !path_begin.contains(it->second),
+        std::prev(it)->second, it->second,
+        !path_begin.contains(std::distance(path.begin(), it)),
         (it->first - current_distance) / (it->first - std::prev(it)->first),
         resolution[1], frequency_range
       ));
@@ -216,19 +218,19 @@ void ufo::plot_band(std::string config_file)
     f->save(filename, "png");
   };
 
-  auto input = YAML::LoadFile(config_file).as<Input>();
+  auto config = YAML::LoadFile(config_file).as<Config>();
-  auto unfolded_data = biu::deserialize<UnfoldOutput>
+  auto input = biu::deserialize<CommonData>
-    (biu::read<std::byte>(input.UnfoldedDataFile));
+    (biu::read<std::byte>(config.InputDataFile));
   
   // 搜索画图需要用到的 q 点
   // key 到起点的距离，value 为 q 点在 QpointData 中的索引
-  boost::container::flat_map<double, std::size_t> path;
+  boost::container::flat_map<double, std::pair<std::size_t, std::size_t>> path;
   // 每一条连续路径的第一个 q 点在 path 中的索引
   std::set<std::size_t> path_begin;
   // x 轴的刻度，为 path 中的索引
   std::set<std::size_t> x_ticks_index;
   double total_distance = 0;
-  for (auto& line : input.Qpoints)
+  for (auto& line : config.Qpoints)
   {
     assert(line.size() >= 2);
     path_begin.insert(path.size());
@@ -237,11 +239,11 @@ void ufo::plot_band(std::string config_file)
       x_ticks_index.insert(path.size());
       auto [this_path, this_distance] = search_qpoints
       (
-        unfolded_data.PrimativeCell, {line[i], line[i + 1]},
+        input.Primative.Cell, {line[i], line[i + 1]},
-        unfolded_data.QpointData
+        input.Super.Qpoint
-          | ranges::views::transform(&UnfoldOutput::QpointDataType::Qpoint)
+          | ranges::views::transform([](auto& qpoint) { return qpoint.SubQpoint; })
           | ranges::to_vector,
-        input.ThresholdWhenSearchingQpoints.value_or(0.001),
+        config.ThresholdWhenSearchingQpoints.value_or(0.001),
         i != line.size() - 2
       );
       path.merge
@@ -258,39 +260,59 @@ void ufo::plot_band(std::string config_file)
   // 计算画图的数据
   auto values = calculate_values
   (
-    path, path_begin, unfolded_data.QpointData, input.InterpolationResolution,
+    path, path_begin,
-    input.FrequencyRange, total_distance
+    input.Super.Qpoint
+      | ranges::views::transform([](auto& qpoint)
+      {
+        return qpoint.Mode | ranges::views::transform([](auto&& mode) { return mode.Frequency; }) | ranges::to_vector;
+      })
+      | ranges::to_vector,
+    input.Super.Qpoint
+      | ranges::views::transform([](auto& qpoint)
+      {
+        return qpoint.Mode
+          | ranges::views::transform([](auto&& mode)
+          {
+            return mode.WeightOnUnfold
+              | ranges::views::transform([&](auto&& weight)
+                { return weight * mode.WeightOnSelectedAtom.value_or(1) * mode.WeightOnRaman.value_or(1); })
+              | ranges::to_vector;
+          })
+          | ranges::to_vector;
+      })
+      | ranges::to_vector,
+    config.InterpolationResolution, config.FrequencyRange, total_distance
   );
   auto x_ticks = x_ticks_index | ranges::views::transform([&](auto i)
-    { return path.nth(i)->first / total_distance * input.InterpolationResolution[0]; }) | ranges::to<std::vector>;
+    { return path.nth(i)->first / total_distance * config.InterpolationResolution[0]; }) | ranges::to_vector;
-  auto y_ticks = input.YTicks.value_or(std::vector<std::pair<double, std::string>>{})
+  auto y_ticks = config.YTicks.value_or(std::vector<std::pair<double, std::string>>{})
     | biu::toLvalue | ranges::views::keys
     | ranges::views::transform([&](auto i)
     {
-      return (i - input.FrequencyRange[0]) / (input.FrequencyRange[1] - input.FrequencyRange[0])
+      return (i - config.FrequencyRange[0]) / (config.FrequencyRange[1] - config.FrequencyRange[0])
-        * input.InterpolationResolution[1];
+        * config.InterpolationResolution[1];
     })
     | ranges::to_vector;
-  auto y_ticklabels = input.YTicks.value_or(std::vector<std::pair<double, std::string>>{})
+  auto y_ticklabels = config.YTicks.value_or(std::vector<std::pair<double, std::string>>{})
     | biu::toLvalue | ranges::views::values | ranges::to_vector;
-  if (input.OutputPictureFile) plot
+  if (config.OutputPictureFile) plot
   (
-    values, input.OutputPictureFile.value(),
+    values, config.OutputPictureFile.value(),
-    x_ticks, y_ticks, y_ticklabels, input.AspectRatio, input.PictureResolution
+    x_ticks, y_ticks, y_ticklabels, config.AspectRatio, config.PictureResolution
   );
-  if (input.OutputDataFile)
+  if (config.OutputDataFile)
-    biu::Hdf5file(input.OutputDataFile.value(), true)
+    biu::Hdf5file(config.OutputDataFile.value(), true)
       .write("Values", values)
       .write("XTicks", x_ticks)
       .write("YTicks", y_ticks)
       .write("YTickLabels", y_ticklabels)
-      .write("InterpolationResolution", input.InterpolationResolution)
+      .write("InterpolationResolution", config.InterpolationResolution)
-      .write("FrequencyRange", input.FrequencyRange);
+      .write("FrequencyRange", config.FrequencyRange);
 }
 
 void ufo::plot_point(std::string config_file)
 {
-  struct Input
+  struct Config
   {
     std::string UnfoldedDataFile;
     // 要画图的 q 点
@@ -317,19 +339,20 @@ void ufo::plot_point(std::string config_file)
   auto search_qpoints = []
   (
     const Eigen::Matrix3d& primative_cell,
-    const Eigen::Vector3d& qpoint, const std::vector<Eigen::Vector3d>& qpoints,
+    const Eigen::Vector3d& qpoint, const std::vector<std::vector<Eigen::Vector3d>>& qpoints,
     double threshold
   )
   {
     biu::Logger::Guard log(qpoint);
     // 对于 output 中的每一个点, 检查这个点是否与所寻找的点足够近，如果足够近则返回
-    for (std::size_t i = 0; i < qpoints.size(); i++)
+    for (std::size_t i = 0; i < qpoints.size(); i++) for (std::size_t j = 0; j < qpoints[i].size(); j++)
       for (auto cell_shift
         : biu::sequence(Eigen::Vector3i(-1, -1, -1), Eigen::Vector3i(2, 2, 2)))
       {
         auto this_qpoint
-          = (primative_cell.reverse().transpose() * (qpoints[i] + cell_shift.first.cast<double>())).eval();
+          = (primative_cell.reverse().transpose() * (qpoints[i][j] + cell_shift.first.cast<double>())).eval();
-        if ((this_qpoint - primative_cell.reverse().transpose() * qpoint).norm() < threshold) return log.rtn(i);
+        if ((this_qpoint - primative_cell.reverse().transpose() * qpoint).norm() < threshold)
+          return log.rtn(std::pair(i, j));
       }
     throw std::runtime_error("No q points found");
   };
@@ -338,7 +361,8 @@ void ufo::plot_point(std::string config_file)
   auto calculate_values = []
   (
     // q 点的数据（需要用到它的频率和权重）
-    const UnfoldOutput::QpointDataType& qpoint,
+    const std::vector<double>& frequency,
+    const std::vector<double>& weight,
     // 用于插值的分辨率和范围
     std::size_t resolution,
     const std::array<double, 2>& frequency_range
@@ -346,12 +370,12 @@ void ufo::plot_point(std::string config_file)
   {
     biu::Logger::Guard log;
     std::vector<double> result(resolution);
-    for (auto& mode : qpoint.ModeData)
+    for (std::size_t i = 0; i < frequency.size(); i++)
     {
-      double index_double = (mode.Frequency - frequency_range[0]) / (frequency_range[1] - frequency_range[0])
+      double index_double = (frequency[i] - frequency_range[0]) / (frequency_range[1] - frequency_range[0])
         * (resolution - 1);
       std::ptrdiff_t index = std::round(index_double);
-      if (index >= 0 && index < static_cast<std::ptrdiff_t>(resolution)) result[index] += mode.Weight;
+      if (index >= 0 && index < static_cast<std::ptrdiff_t>(resolution)) result[index] += weight[i];
     }
     return log.rtn(result);
   };
@@ -387,45 +411,55 @@ void ufo::plot_point(std::string config_file)
   };
 
   biu::Logger::Guard log;
-  auto input = YAML::LoadFile(config_file).as<Input>();
+  auto config = YAML::LoadFile(config_file).as<Config>();
-  auto unfolded_data = biu::deserialize<UnfoldOutput>
+  auto input = biu::deserialize<CommonData>
-    (biu::read<std::byte>(input.UnfoldedDataFile));
+    (biu::read<std::byte>(config.UnfoldedDataFile));
   
   auto qpoint_index = search_qpoints
   (
-    unfolded_data.PrimativeCell, input.Qpoint,
+    input.Primative.Cell, config.Qpoint,
-    unfolded_data.QpointData
+    input.Super.Qpoint
-      | ranges::views::transform(&UnfoldOutput::QpointDataType::Qpoint)
+      | ranges::views::transform([](auto& qpoint) { return qpoint.SubQpoint; })
       | ranges::to_vector,
-    input.ThresholdWhenSearchingQpoints.value_or(0.001)
+    config.ThresholdWhenSearchingQpoints.value_or(0.001)
   );
   auto values = calculate_values
   (
-    unfolded_data.QpointData[qpoint_index],
+    input.Super.Qpoint[qpoint_index.first].Mode
-    input.InterpolationResolution, input.FrequencyRange
+      | ranges::views::transform([](auto&& mode) { return mode.Frequency; })
+      | ranges::to_vector,
+    input.Super.Qpoint[qpoint_index.first].Mode
+      | ranges::views::transform([&](auto&& mode)
+      {
+        return mode.WeightOnUnfold[qpoint_index.second]
+          * mode.WeightOnSelectedAtom.value_or(1) * mode.WeightOnRaman.value_or(1);
+      })
+      | ranges::to_vector,
+    config.InterpolationResolution, config.FrequencyRange
   );
-  auto x_ticks = input.XTicks.value_or(std::vector<std::pair<double, std::string>>{})
+  auto x_ticks = config.XTicks.value_or(std::vector<std::pair<double, std::string>>{})
     | biu::toLvalue | ranges::views::keys
     | ranges::views::transform([&](auto i)
     {
-      return (i - input.FrequencyRange[0]) / (input.FrequencyRange[1] - input.FrequencyRange[0])
+      return (i - config.FrequencyRange[0]) / (config.FrequencyRange[1] - config.FrequencyRange[0])
-        * input.InterpolationResolution;
+        * config.InterpolationResolution;
     })
     | ranges::to_vector;
-  auto x_ticklabels = input.XTicks.value_or(std::vector<std::pair<double, std::string>>{})
+  auto x_ticklabels = config.XTicks.value_or(std::vector<std::pair<double, std::string>>{})
     | biu::toLvalue | ranges::views::values | ranges::to_vector;
-  if (input.OutputPictureFile) plot
+  if (config.OutputPictureFile) plot
   (
-    values, input.OutputPictureFile.value(),
+    values, config.OutputPictureFile.value(),
-    x_ticks, x_ticklabels, input.AspectRatio, input.PictureResolution
+    x_ticks, x_ticklabels, config.AspectRatio, config.PictureResolution
   );
-  if (input.OutputDataFile)
+  if (config.OutputDataFile)
-    biu::Hdf5file(input.OutputDataFile.value(), true)
+    biu::Hdf5file(config.OutputDataFile.value(), true)
       .write("Values", values)
       .write("XTicks", x_ticks)
       .write("XTickLabels", x_ticklabels)
-      .write("InterpolationResolution", input.InterpolationResolution)
+      .write("InterpolationResolution", config.InterpolationResolution)
-      .write("FrequencyRange", input.FrequencyRange);
+      .write("FrequencyRange", config.FrequencyRange);
-  if (input.OutputRawDataFile)
+  // TODO: pfr + optional + Eigen Matrix = failed
-    std::ofstream(*input.OutputRawDataFile) << YAML::Node(unfolded_data.QpointData[qpoint_index]);
+  if (config.OutputRawDataFile)
+    std::ofstream(*config.OutputRawDataFile) << YAML::Node(values);
 }

src/project-to-atom.cpp

@@ -0,0 +1,37 @@
+# include <ufo.hpp>
+# include <ranges>
+
+void ufo::project_to_atom(std::string config_file)
+{
+  struct Config
+  {
+    std::set<std::size_t> SelectedAtom;
+    std::string InputFile;
+    std::string OutputFile;
+  };
+
+  biu::Logger::Guard log(config_file);
+  auto config = YAML::LoadFile(config_file).as<Config>();
+  auto input = biu::deserialize<CommonData>(biu::read<std::byte>(config.InputFile));
+  input.SelectedAtom = config.SelectedAtom;
+  for (auto& qpoint : input.Super.Qpoint) for (auto& mode : qpoint.Mode)
+    mode.WeightOnSelectedAtom = mode.EigenVector.rowwise().norm().cwiseProduct
+    (
+      ranges::views::iota(0ul, input.Super.AtomPosition.rows() * 1ul)
+        | ranges::views::transform([&](auto i)
+          { return config.SelectedAtom.contains(i) ? 1. : 0.; })
+        | ranges::to_vector
+        | biu::toEigen<>
+    ).sum() * (input.Super.AtomPosition.rows() * 1. / config.SelectedAtom.size());
+  std::ofstream(config.OutputFile, std::ios::binary) << biu::serialize<char>(input);
+  log.info("Summary:");
+  for (auto i_of_qpoint : std::views::iota(0ul, input.Super.Qpoint.size()))
+    for (auto i_of_mode : std::views::iota(0ul, input.Super.Qpoint[i_of_qpoint].Mode.size()))
+      if (*input.Super.Qpoint[i_of_qpoint].Mode[i_of_mode].WeightOnSelectedAtom > 0.1)
+        log.info("{}:{:.2f}:{}:{:.2f}"_f
+        (
+          i_of_qpoint, fmt::join(input.Super.Qpoint[i_of_qpoint].Qpoint, ", "),
+          i_of_mode, input.Super.Qpoint[i_of_qpoint].Mode[i_of_mode].Frequency,
+          *input.Super.Qpoint[i_of_qpoint].Mode[i_of_mode].WeightOnSelectedAtom
+        ));
+}

src/raman-apply-contribution.cpp

@@ -1,104 +0,0 @@
-# include <ufo.hpp>
-
-void ufo::raman_apply_contribution(std::string config_file)
-{
-  struct Input
-  {
-    std::string UnfoldedDataFile;
-    std::string DisplacementDataFile;
-    Eigen::Matrix3d OriginalSusceptibility;
-    std::vector<Eigen::Matrix3d> Susceptibilities;
-    std::string OutputDataFile;
-    std::array<Eigen::Vector3d, 2> Polarization;
-  };
-
-  biu::Logger::Guard log(config_file);
-  auto input = YAML::LoadFile(config_file).as<Input>();
-  auto unfolded_data = biu::deserialize<UnfoldOutput>
-    (biu::read<std::byte>(input.UnfoldedDataFile));
-  auto displacement_data = biu::deserialize<DisplacementOutput>
-    (biu::read<std::byte>(input.DisplacementDataFile));
-  UnfoldOutput output;
-
-  // 整理得到的数据，作放缩
-  struct mode_t
-  {
-    std::size_t MetaQpointIndex;
-    std::size_t ModeIndex;
-    Eigen::Matrix3d RamanTensor;
-    double Strength;
-  };
-  auto modes = ranges::views::zip(displacement_data.ModeData, input.Susceptibilities)
-    | ranges::views::transform([&](const auto& data)
-      {
-        auto& [mode, susceptibility] = data;
-        Eigen::Matrix3d raman_tensor = susceptibility / mode.Ratio / displacement_data.MaxDisplacement;
-        double intensity = (input.Polarization[0].transpose() * raman_tensor * input.Polarization[1]).norm();
-        log.debug("get raman tensor: {}"_f(raman_tensor));
-        log.debug("get intensity: {}"_f(intensity));
-        return mode_t{ mode.MetaQpointIndex, mode.ModeIndex, raman_tensor, intensity };
-      })
-    | ranges::to_vector;
-
-  // 整理输出文件
-  biu::for_each
-  (
-    [&](auto&& i) { output.*i = unfolded_data.*i; },
-    std::tuple(&UnfoldOutput::PrimativeCell, &UnfoldOutput::SuperCellTransformation,
-      &UnfoldOutput::SuperCellMultiplier, &UnfoldOutput::SuperCellDeformation, &UnfoldOutput::SelectedAtoms)
-  );
-  // 将以前的 meta qpoint 和 mode 的索引转换到新的
-  std::map<std::size_t, std::size_t> meta_qpoint_map;
-  std::vector<std::map<std::size_t, std::size_t>> mode_map;
-  for (auto mode : displacement_data.ModeData)
-  {
-    auto old_meta_qpoint_index = mode.MetaQpointIndex;
-    auto old_mode_index = mode.ModeIndex;
-    auto new_meta_qpoint_index = meta_qpoint_map.contains(old_meta_qpoint_index)
-      ? meta_qpoint_map[mode.MetaQpointIndex]
-      : [&]
-      {
-        log.debug("map meta qpoint {} to {}"_f(old_meta_qpoint_index, meta_qpoint_map.size()));
-        meta_qpoint_map[mode.MetaQpointIndex] = meta_qpoint_map.size();
-        output.MetaQpointData.push_back({unfolded_data.MetaQpointData[old_mode_index].Qpoint, {}});
-        return meta_qpoint_map.size() - 1;
-      }();
-    auto new_mode_index = mode_map[new_meta_qpoint_index].size();
-    log.debug("map mode {} {} to {}"_f(old_meta_qpoint_index, old_mode_index, new_mode_index));
-    mode_map[new_meta_qpoint_index][old_mode_index] = new_mode_index;
-    output.MetaQpointData[new_meta_qpoint_index].ModeData.push_back
-      (unfolded_data.MetaQpointData[old_meta_qpoint_index].ModeData[old_mode_index]);
-  }
-  for (auto old_qpoint_index : displacement_data.QpointIndices)
-    if (meta_qpoint_map.contains(unfolded_data.QpointData[old_qpoint_index].SourceIndex))
-    {
-      log.debug("map qpoint {} {}"_f(old_qpoint_index, unfolded_data.QpointData[old_qpoint_index].Qpoint));
-      auto new_meta_qpoint_index
-        = meta_qpoint_map[unfolded_data.QpointData[old_qpoint_index].SourceIndex];
-      output.QpointData.push_back
-      ({
-        unfolded_data.QpointData[old_qpoint_index].Qpoint,
-        unfolded_data.QpointData[old_qpoint_index].Source,
-        new_meta_qpoint_index,
-        {}
-      });
-      for (std::size_t i = 0; i < unfolded_data.QpointData[old_qpoint_index].ModeData.size(); i++)
-        if (mode_map[new_meta_qpoint_index].contains(i))
-        {
-          auto new_mode_index = mode_map[new_meta_qpoint_index][i];
-          log.debug("map mode {} to {}"_f(i, new_mode_index));
-          if (output.QpointData.back().ModeData.size() <= new_mode_index)
-          {
-            log.debug("resize to {}"_f(new_mode_index + 1));
-            output.QpointData.back().ModeData.resize(new_mode_index + 1);
-          }
-          output.QpointData.back().ModeData[new_mode_index] =
-          {
-            unfolded_data.QpointData[old_qpoint_index].ModeData[i].Frequency,
-            unfolded_data.QpointData[old_qpoint_index].ModeData[i].Weight * modes[new_mode_index].Strength
-          };
-        }
-    }
-
-  std::ofstream(input.OutputDataFile, std::ios::binary) << biu::serialize<char>(output);
-}

src/raman-create-displacement.cpp

@@ -1,128 +0,0 @@
-# include <ufo.hpp>
-
-void ufo::raman_create_displacement(std::string config_file)
-{
-  struct Input
-  {
-    std::string UnfoldedDataFile;
-    // 搜索位于 Gamma 点的 q 点时，使用的阈值，单位为埃^-1，默认为 0.01
-    std::optional<double> ThresholdWhenSearchingQpoints;
-    // 搜索权重非零的模式时，使用的阈值，默认为 0.01
-    std::optional<double> ThresholdWhenSearchingModes;
-    // 所有原子的符号
-    std::vector<std::string> AtomSymbols;
-    // 各种原子的质量，单位为原子质量
-    std::map<std::string, double> AtomMasses;
-    // 原子最大位移大小，单位为埃
-    double MaxDisplacement;
-    // 超胞，单位为埃
-    Eigen::Matrix3d SuperCell;
-    // 超胞中各个原子的坐标，单位为超胞的格矢
-    Eigen::MatrixX3d AtomPositions;
-    // 输出的POSCAR所在的目录
-    std::string OutputPoscarDirectory;
-    // 输出的数据文件名
-    std::string OutputDataFile;
-  };
-
-  // 假定同类型的原子一定写在一起
-  auto generate_poscar = []
-  (
-    Eigen::Matrix3d SuperCell, Eigen::MatrixX3d AtomPositions, std::vector<std::string> AtomSymbols
-  )
-  {
-    std::stringstream ss;
-    ss << "some random comment to make VASP happy\n1.0\n";
-    for (std::size_t i = 0; i < 3; i++)
-    {
-      for (std::size_t j = 0; j < 3; j++) ss << SuperCell(i, j) << " ";
-      ss << std::endl;
-    }
-    auto atom_symbols = AtomSymbols | ranges::views::chunk_by(std::ranges::equal_to{});
-    ss << "{}\n"_f(ranges::accumulate
-    (
-      atom_symbols | ranges::views::transform([](auto&& chunk) { return chunk[0]; }),
-      ""s, [](auto&& a, auto&& b) { return a + " " + b; }
-    ));
-    ss << "{}\n"_f(ranges::accumulate
-    (
-      atom_symbols | ranges::views::transform([](auto&& chunk) { return chunk.size(); }),
-      ""s, [](auto&& a, auto&& b) { return a + " " + std::to_string(b); }
-    ));
-    ss << "Direct\n";
-    for (const auto& position : AtomPositions.rowwise())
-    {
-      for (std::size_t i = 0; i < 3; i++) ss << position(i) << " ";
-      ss << std::endl;
-    }
-    return ss.str();
-  };
-
-  biu::Logger::Guard log(config_file);
-  auto input = YAML::LoadFile(config_file).as<Input>();
-  auto unfolded_data = biu::deserialize<UnfoldOutput>
-    (biu::read<std::byte>(input.UnfoldedDataFile));
-  DisplacementOutput output;
-
-  // 搜索满足条件的模式，找到满足条件的模式后，就将 MetaQpoint 的索引加入到 output 中
-  // 之所以使用 MetaQpoint 的索引而不是 Qpoint 的索引，是因为 Qpoint 中可能有指向同一个 MetaQpoint 中模式的不同模式都满足要求
-  // 如果写入 Qpoint 的索引，就会重复而增加之后的计算量
-  std::set<std::pair<std::size_t, std::size_t>> selected_modes;
-  for (const auto& qpoint : unfolded_data.QpointData)
-  {
-    if
-    (
-      (unfolded_data.PrimativeCell.reverse().transpose() * qpoint.Qpoint).norm()
-        > input.ThresholdWhenSearchingQpoints.value_or(0.01)
-    )
-      continue;
-    for (std::size_t i = 0; i < qpoint.ModeData.size(); i++)
-    {
-      if (qpoint.ModeData[i].Weight < input.ThresholdWhenSearchingModes.value_or(0.01)) continue;
-      log.debug("Found mode {} {} frequency {}"_f(qpoint.Qpoint, i, qpoint.ModeData[i].Frequency));
-      selected_modes.insert({qpoint.SourceIndex, i});
-    }
-  }
-
-  // 构造输出数据
-  output.AtomMasses = input.AtomSymbols
-    | ranges::views::transform([&](const auto& symbol)
-      { return input.AtomMasses.at(symbol); })
-    | ranges::to_vector
-    | biu::toEigen<>;
-  output.MaxDisplacement = input.MaxDisplacement;
-  for (auto [i, j] : selected_modes)
-  {
-    auto& mode_data = output.ModeData.emplace_back();
-    mode_data.MetaQpointIndex = i;
-    mode_data.ModeIndex = j;
-    // 未归一化的位移, 假定虚部总是为零
-    auto atom_movement =
-      unfolded_data.MetaQpointData[i].ModeData[j].AtomMovement.real()
-        .cwiseProduct(output.AtomMasses.cwiseSqrt().cwiseInverse().rowwise().replicate(3)).eval();
-    // 归一化
-    mode_data.Ratio = input.MaxDisplacement / atom_movement.rowwise().norm().maxCoeff();
-    mode_data.AtomMovement = atom_movement * mode_data.Ratio;
-    log.debug("Write mode {} {} {}"_f(i, j, mode_data.AtomMovement.rowwise().norm().transpose()));
-  }
-
-  // 输出
-  std::ofstream(input.OutputDataFile, std::ios::binary) << biu::serialize<char>(output);
-  for (std::size_t i = 0; i < output.ModeData.size(); i++)
-  {
-    std::filesystem::create_directories(input.OutputPoscarDirectory + "/" + std::to_string(i));
-    std::ofstream(input.OutputPoscarDirectory + "/" + std::to_string(i) + "/POSCAR") << generate_poscar
-    (
-      input.SuperCell,
-      input.AtomPositions + output.ModeData[i].AtomMovement * input.SuperCell.inverse(),
-      input.AtomSymbols
-    );
-  }
-  std::filesystem::create_directories(input.OutputPoscarDirectory + "/original");
-  std::ofstream(input.OutputPoscarDirectory + "/original/POSCAR") << generate_poscar
-  (
-    input.SuperCell,
-    input.AtomPositions,
-    input.AtomSymbols
-  );
-}

src/raman-extract.cpp

@@ -1,42 +0,0 @@
-# include <ufo.hpp>
-
-void ufo::raman_extract(std::vector<std::string> files)
-{
-  biu::Logger::Guard log(files);
-  std::vector<Eigen::Matrix3d> electricity_tensors;
-
-  for (const auto& file : files)
-  {
-    auto in_stream = std::ifstream(file);
-    std::string line;
-    // search line containing: MACROSCOPIC STATIC DIELECTRIC TENSOR
-    while (std::getline(in_stream, line))
-    {
-      if (line.find("MACROSCOPIC STATIC DIELECTRIC TENSOR") != std::string::npos)
-      {
-        log.debug("find in {}"_f(file));
-        // skip 1 line
-        std::getline(in_stream, line);
-        electricity_tensors.emplace_back();
-        for (std::size_t i = 0; i < 3; i++) for (std::size_t j = 0; j < 3; j++)
-          in_stream >> electricity_tensors.back()(i, j);
-        log.debug("read tensor {}"_f(electricity_tensors.back()));
-        break;
-      }
-    }
-    // test if the file is read correctly
-    if (!in_stream) throw std::runtime_error("Error reading file: {}"_f(file));
-  }
-
-  // output the result
-  for (auto e: electricity_tensors) std::cout <<
-R"(- - [ {}, {}, {} ]
-  - [ {}, {}, {} ]
-  - [ {}, {}, {} ]
-)"_f
-    (
-      e(0, 0), e(0, 1), e(0, 2),
-      e(1, 0), e(1, 1), e(1, 2),
-      e(2, 0), e(2, 1), e(2, 2)
-    );
-}

src/raman.cpp

@@ -0,0 +1,199 @@
+
+# include <ufo.hpp>
+
+namespace ufo
+{
+  struct RamanData
+  {
+    struct ModeType
+    {
+      std::size_t QpointIndex;
+      std::size_t ModeIndex;
+      double Ratio;
+      using serialize = zpp::bits::members<3>;
+    };
+    std::vector<ModeType> Mode;
+    double MaxDisplacement;
+    enum { Primative, Super } Cell;
+    using serialize = zpp::bits::members<3>;
+  };
+
+  void raman_create_displacement(std::string config_file)
+  {
+    struct Config
+    {
+      // 要计算的是原胞还是超胞
+      decltype(RamanData::Cell) Cell;
+      // 要计算的模式，总是假定是 gamma 点的模式
+      std::map<std::size_t, std::set<std::size_t>> SelectedModes;
+      // 原子最大位移大小，单位为埃
+      double MaxDisplacement;
+      // 输出的POSCAR所在的目录
+      std::string OutputPoscarDirectory;
+      // 输入的数据文件名
+      std::string InputDataFile;
+      std::string OutputDataFile;
+    };
+
+    auto generate_poscar = []
+    (
+      Eigen::Matrix3d SuperCell, Eigen::MatrixX3d AtomPosition,
+      std::vector<std::pair<std::string, std::size_t>> AtomType
+    )
+    {
+      std::stringstream ss;
+      ss << "some random comment to make VASP happy\n1.0\n";
+      for (std::size_t i = 0; i < 3; i++)
+      {
+        for (std::size_t j = 0; j < 3; j++) ss << SuperCell(i, j) << " ";
+        ss << std::endl;
+      }
+      ss << "{}\n"_f(ranges::accumulate
+      (
+        AtomType | ranges::views::transform([](auto&& atom) { return atom.first; }),
+        ""s, [](auto&& a, auto&& b) { return a + " " + b; }
+      ));
+      ss << "{}\n"_f(ranges::accumulate
+      (
+        AtomType | ranges::views::transform([](auto&& atom) { return atom.second; }),
+        ""s, [](auto&& a, auto&& b) { return a + " " + std::to_string(b); }
+      ));
+      ss << "Direct\n";
+      for (const auto& position : AtomPosition.rowwise())
+      {
+        for (std::size_t i = 0; i < 3; i++) ss << position(i) << " ";
+        ss << std::endl;
+      }
+      return ss.str();
+    };
+
+    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));
+    RamanData output;
+    output.MaxDisplacement = config.MaxDisplacement;
+    output.Cell = config.Cell;
+
+    auto process = [&](auto& cell)
+    {
+      std::size_t i_of_poscar = 0;
+      auto mass = cell.AtomType
+        | ranges::views::transform([&](auto&& atom)
+          { return ranges::views::repeat_n(input.AtomMass[atom.first], atom.second); })
+        | ranges::views::join | ranges::to_vector | biu::toEigen<>;
+      for (auto i_of_qpoint : config.SelectedModes | ranges::views::keys)
+        for (auto i_of_mode : config.SelectedModes[i_of_qpoint])
+        {
+          // 假定虚部总是为零
+          auto atom_movement = cell.Qpoint[i_of_qpoint].Mode[i_of_mode].EigenVector.real()
+            .cwiseProduct(mass.cwiseSqrt().cwiseInverse().rowwise().replicate(3)).eval();
+          // 归一化
+          auto ratio = config.MaxDisplacement / atom_movement.rowwise().norm().maxCoeff();
+          atom_movement *= ratio;
+          // 输出
+          auto path = "{}/{}"_f(config.OutputPoscarDirectory, i_of_poscar);
+          std::filesystem::create_directories(path);
+          std::ofstream("{}/POSCAR"_f(path)) << generate_poscar
+          (
+            cell.Cell,
+            cell.AtomPosition + atom_movement * cell.Cell.inverse(),
+            cell.AtomType
+          );
+          output.Mode.push_back({i_of_qpoint, i_of_mode, ratio});
+          log.debug("Write mode {} {} {}"_f(i_of_qpoint, i_of_mode, atom_movement.rowwise().norm().eval()));
+          i_of_poscar++;
+        }
+    };
+    if (config.Cell == RamanData::Primative) process(input.Primative);
+    else process(input.Super);
+
+    std::filesystem::create_directories("{}/{}"_f(config.OutputPoscarDirectory, "origin"));
+    std::ofstream("{}/origin/POSCAR"_f(config.OutputDataFile)) << generate_poscar
+      (input.Super.Cell, input.Super.AtomPosition, input.Super.AtomType);
+    std::ofstream(config.OutputDataFile, std::ios::binary) << biu::serialize<char>(output);
+  }
+
+  void raman_extract(std::vector<std::string> files)
+  {
+    biu::Logger::Guard log(files);
+    std::vector<Eigen::Matrix3d> electricity_tensors;
+
+    for (const auto& file : files)
+    {
+      auto in_stream = std::ifstream(file);
+      std::string line;
+      // search line containing: MACROSCOPIC STATIC DIELECTRIC TENSOR
+      while (std::getline(in_stream, line))
+      {
+        if (line.find("MACROSCOPIC STATIC DIELECTRIC TENSOR") != std::string::npos)
+        {
+          log.debug("find in {}"_f(file));
+          // skip 1 line
+          std::getline(in_stream, line);
+          electricity_tensors.emplace_back();
+          for (std::size_t i = 0; i < 3; i++) for (std::size_t j = 0; j < 3; j++)
+            in_stream >> electricity_tensors.back()(i, j);
+          log.debug("read tensor {}"_f(electricity_tensors.back()));
+          break;
+        }
+      }
+      // test if the file is read correctly
+      if (!in_stream) throw std::runtime_error("Error reading file: {}"_f(file));
+    }
+
+    // output the result
+    for (auto e: electricity_tensors) std::cout <<
+R"(- - [ {}, {}, {} ]
+  - [ {}, {}, {} ]
+  - [ {}, {}, {} ]
+)"_f
+      (
+        e(0, 0), e(0, 1), e(0, 2),
+        e(1, 0), e(1, 1), e(1, 2),
+        e(2, 0), e(2, 1), e(2, 2)
+      );
+  }
+
+  void raman_apply_contribution(std::string config_file)
+  {
+    struct Config
+    {
+      Eigen::Matrix3d OriginalSusceptibility;
+      std::vector<Eigen::Matrix3d> Susceptibilities;
+      std::array<Eigen::Vector3d, 2> Polarization;
+      std::string InputDataFile;
+      std::string RamanInputDataFile;
+      std::string OutputDataFile;
+    };
+
+    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));
+    auto raman_input = biu::deserialize<RamanData>
+      (biu::read<std::byte>(config.RamanInputDataFile));
+
+    input.RamanPolarization = config.Polarization;
+    auto process = [&](auto& cell)
+    {
+      for (auto&& [i_of_mode, mode] : ranges::views::enumerate(raman_input.Mode))
+      {
+        auto&& _ = cell.Qpoint[mode.QpointIndex].Mode[mode.ModeIndex];
+        Eigen::Matrix3d raman_tensor = (config.Susceptibilities[i_of_mode] - config.OriginalSusceptibility)
+          / mode.Ratio / raman_input.MaxDisplacement;
+        _.RamanTensor = raman_tensor | biu::fromEigen;
+        _.WeightOnRaman = config.Polarization[0].transpose() * raman_tensor * config.Polarization[1];
+        log.info("{}:{:.2f}:{}:{:.2f} {}"_f
+        (
+          mode.QpointIndex, fmt::join(cell.Qpoint[mode.QpointIndex].Qpoint, ", "),
+          mode.ModeIndex, _.Frequency, *_.WeightOnRaman
+        ));
+      }
+    };
+    if (raman_input.Cell == RamanData::Primative) process(input.Primative);
+    else process(input.Super);
+
+    std::ofstream(config.OutputDataFile, std::ios::binary) << biu::serialize<char>(input);
+  }
+}

src/unfold.cpp

@@ -2,6 +2,7 @@
 # include <thread>
 # include <syncstream>
 # include <execution>
+# include <ranges>
 
 void ufo::unfold(std::string config_file)
 {
@@ -13,11 +14,8 @@ void ufo::unfold(std::string config_file)
   // 只要第一部分取得足够多, 那么单胞中原子的状态就可以完全被这些平面波描述.
   // 将超胞中原子的运动状态投影到这些基矢上, 计算出投影的系数, 就可以将超胞的原子运动状态分解到单胞中的多个 q 点上.
 
-  struct Input
+  struct Config
   {
-    // 单胞的三个格矢，每行表示一个格矢的坐标，单位为埃
-    Eigen::Matrix3d PrimativeCell;
-
     // 单胞到超胞的格矢转换时用到的矩阵
     // SuperCellMultiplier 是一个三维列向量且各个元素都是整数，表示单胞在各个方向扩大到多少倍之后，可以得到和超胞一样的体积
     // SuperCellDeformation 是一个行列式为 1 的矩阵，它表示经过 SuperCellMultiplier 扩大后，还需要怎样的变换才能得到超胞
@@ -35,40 +33,32 @@ void ufo::unfold(std::string config_file)
     // 在单胞内取几个平面波的基矢
     Eigen::Vector<std::size_t, 3> PrimativeCellBasisNumber;
 
-    // 超胞中原子的坐标，每行表示一个原子的坐标，单位为超胞的格矢
+    // 单胞的 phonopy 输出的 phonopy.yaml，用来读入单胞的晶格、原子坐标、原子类型、原子质量
-    Eigen::MatrixX3d AtomPositionBySuperCell;
+    std::string PrimativePhonopy;
+    // 单胞的 phonopy 输出的 band.hdf5 或 qpoint.hdf5，用来读入 q 点和振动模式
+    std::string PrimativeQpoint;
+    // 同上，但是是超胞里的结果
+    std::string SuperPhonopy;
+    std::string SuperQpoint;
 
-    // 从 band.hdf5 读入 QpointData
+    std::string OutputFile;
-    std::optional<std::string> QpointDataInputFile;
-
-    // 输出到哪些文件
-    struct QpointDataOutputFileType
-    {
-      std::string Filename;
-
-      // 如果指定，则将结果投影到那些原子上
-      std::optional<std::vector<std::size_t>> SelectedAtoms;
-
-      // 默认输出为 zpp 文件，如果指定为 true，则输出为 yaml 文件
-      std::optional<bool> OutputAsYaml;
-    };
-    std::vector<QpointDataOutputFileType> QpointDataOutputFile;
   };
 
-  // 从文件中读取 QpointData
+  // 从文件中读取 q 点数据
-  auto read_qpoint_data = [](std::string filename)
+  // data 为 CommonData::Primative 或 CommonData::Super
+  // 返回值为原子类型和原子质量的对应关系
+  auto read_qpoint = [](std::string phonopy_file, std::string qpoint_file, auto& data)
   {
-    // 读入原始数据
+    // phonopy 的输出有两种可能。
-    // phonopy 的输出有两种可能
+    // 直接指定计算的 q 点时，frequency 是 2 维，这时第一个维度是 q 点，第二个维度是不同模式。
-    // 直接指定计算的 q 点时，frequency 是 2 维，这时第一个维度是 q 点，第二个维度是不同模式
+    // 计算能带时，frequency 是 3 维，相比于二维的情况多了第一个维度，表示 q 点所在路径。
-    // 计算能带时，frequency 是 3 维，相比于二维的情况多了第一个维度，表示 q 点所在路径
+    // qpoint 或 path，以及 eigenvector 也有类似的变化。
-    // qpoint 或 path，以及 eigenvector 也有类似的变化
+    // 除此以外，eigenvector 后两个维度分别指示模式的特征向量的各个维度和各个模式（而不是各个模式和特征向量的各个维度），
-    // eigenvector 是三维或四维的数组，后两个维度分别表示原子运动和模式（而不是模式和原子），
+    // 因为后两个维度的尺寸总是一样的（模式个数等于原子坐标个数），非常容易搞错。
-    //  因为后两个维度的尺寸总是一样的（模式个数等于原子坐标个数），非常容易搞错
     std::vector<std::array<double, 3>> qpoint;
     std::vector<std::vector<double>> frequency;
-    std::vector<std::vector<std::vector<biu::PhonopyComplex>>> eigenvector_vector;
+    std::vector<std::vector<std::vector<std::complex<double>>>> eigenvector_vector;
-    auto file = biu::Hdf5file(filename);
+    auto file = biu::Hdf5file(qpoint_file);
 
     if (file.File.getDataSet("/frequency").getDimensions().size() == 2)
       file.read("/frequency", frequency)
@@ -78,7 +68,7 @@ void ufo::unfold(std::string config_file)
     {
       std::vector<std::vector<std::array<double, 3>>> temp_path;
       std::vector<std::vector<std::vector<double>>> temp_frequency;
-      std::vector<std::vector<std::vector<std::vector<biu::PhonopyComplex>>>> temp_eigenvector_vector;
+      std::vector<std::vector<std::vector<std::vector<std::complex<double>>>>> temp_eigenvector_vector;
       file.read("/frequency", temp_frequency)
         .read("/eigenvector", temp_eigenvector_vector)
         .read("/path", temp_path);
@@ -87,27 +77,47 @@ void ufo::unfold(std::string config_file)
       eigenvector_vector = temp_eigenvector_vector | ranges::views::join | ranges::to_vector;
     }
 
-    // 整理得到结果
+    // 整理并写入得到结果
-    auto number_of_qpoints = frequency.size(), num_of_modes = frequency[0].size();
+    auto number_of_qpoints = frequency.size(), number_of_modes = frequency[0].size();
-    std::vector<UnfoldOutput::MetaQpointDataType> qpoint_data(number_of_qpoints);
+    data.Qpoint.resize(number_of_qpoints);
-    for (std::size_t i = 0; i < number_of_qpoints; i++)
+    for (auto i : std::views::iota(0u, number_of_qpoints)) 
     {
-      qpoint_data[i].Qpoint = qpoint[i] | biu::toEigen<>;
+      data.Qpoint[i].Qpoint = qpoint[i] | biu::toEigen<>;
-      qpoint_data[i].ModeData.resize(num_of_modes);
+      data.Qpoint[i].Mode.resize(number_of_modes);
-      for (std::size_t j = 0; j < num_of_modes; j++)
+      for (auto j : std::views::iota(0u, number_of_modes)) 
       {
-        qpoint_data[i].ModeData[j].Frequency = frequency[i][j];
+        data.Qpoint[i].Mode[j].Frequency = frequency[i][j];
-        auto number_of_atoms = eigenvector_vector[i].size() / 3;
+        auto number_of_atoms = number_of_modes / 3;
-        Eigen::MatrixX3cd eigenvectors(number_of_atoms, 3);
+        Eigen::MatrixX3cd eigenvector(number_of_atoms, 3);
-        for (std::size_t k = 0; k < number_of_atoms; k++) for (std::size_t l = 0; l < 3; l++)
+        for (auto k : std::views::iota(0u, number_of_atoms))
-          eigenvectors(k, l)
+          for (auto l : std::views::iota(0u, 3u))
-            = eigenvector_vector[i][k * 3 + l][j].r + eigenvector_vector[i][k * 3 + l][j].i * 1i;
+            eigenvector(k, l) = eigenvector_vector[i][k * 3 + l][j];
         // 原则上讲，需要对读入的原子运动状态作相位转换, 使得它们与我们的约定一致(对超胞周期性重复)，但这个转换 phonopy 已经做了
         // 这里还要需要做归一化处理 (指将数据简单地作为向量处理的归一化)
-        qpoint_data[i].ModeData[j].AtomMovement = eigenvectors / eigenvectors.norm();
+        data.Qpoint[i].Mode[j].EigenVector = eigenvector / eigenvector.norm();
       }
     }
-    return qpoint_data;
+
+    // 读取并写入其它数据
+    YAML::Node phonopy = YAML::LoadFile(phonopy_file);
+    data.Cell = phonopy["unit_cell"]["lattice"].as<std::array<std::array<double, 3>, 3>>() | biu::toEigen<>;
+    auto points = phonopy["points"].as<std::vector<YAML::Node>>();
+    data.AtomType = points
+      | ranges::views::transform([](const YAML::Node& point) { return point["symbol"].as<std::string>(); })
+      | ranges::views::chunk_by(std::ranges::equal_to{})
+      | ranges::views::transform([](const auto& chunk) { return std::pair{chunk[0], chunk.size()}; })
+      | ranges::to_vector;
+    data.AtomPosition = points
+      | ranges::views::transform([](const YAML::Node& point)
+        { return point["coordinates"].as<std::array<double, 3>>(); })
+      | ranges::to_vector
+      | biu::toEigen<>;
+    return points
+      | ranges::views::transform([](const YAML::Node& point)
+        { return std::pair(point["symbol"].as<std::string>(), point["mass"].as<double>()); })
+      | ranges::views::chunk_by(std::ranges::equal_to{})
+      | ranges::views::transform([](const auto& chunk) { return chunk[0]; })
+      | ranges::to<std::map<std::string, double>>;
   };
 
   // 构建基
@@ -151,31 +161,25 @@ void ufo::unfold(std::string config_file)
   auto construct_projection_coefficient = []
   (
     const std::vector<std::vector<Eigen::VectorXcd>>& basis,
-    // 实际上只需要其中的 AtomMovement
+    const std::vector<std::reference_wrapper<const Eigen::MatrixX3cd>>& modes,
-    const std::vector<UnfoldOutput::MetaQpointDataType>& qpoint_data,
     std::atomic<std::size_t>& number_of_finished_modes
   )
   {
-    // 将所有的模式取出，组成一个一维数组，稍后并行计算
-    std::vector<std::reference_wrapper<const Eigen::MatrixX3cd>> mode_data;
-    for (auto& qpoint : qpoint_data) for (auto& mode : qpoint.ModeData)
-      mode_data.emplace_back(mode.AtomMovement);
     // 第一层下标对应不同模式, 第二层下标对应这个模式在反折叠后的 q 点(sub qpoint)
-    std::vector<std::vector<double>> projection_coefficient(mode_data.size());
+    std::vector<std::vector<double>> projection_coefficient(modes.size());
     // 对每个模式并行
     std::transform
     (
-      std::execution::par_unseq, mode_data.begin(), mode_data.end(),
+      std::execution::par, modes.begin(), modes.end(),
       projection_coefficient.begin(), [&](const auto& mode_data)
       {
         // 这里, mode_data 和 projection_coefficient 均指对应于一个模式的数据
         std::vector<double> projection_coefficient(basis.size());
-        for (std::size_t i_of_sub_qpoint = 0; i_of_sub_qpoint < basis.size(); i_of_sub_qpoint++)
+        for (auto i_of_sub_qpoint : std::views::iota(0u, basis.size()))
           // 对于 basis 中, 对应于单胞倒格子的部分, 以及对应于不同方向的部分, 分别求内积, 然后求模方和
-          for (std::size_t i_of_basis = 0; i_of_basis < basis[i_of_sub_qpoint].size(); i_of_basis++)
+          for (auto i_of_basis : std::views::iota(0u, basis[i_of_sub_qpoint].size()))
             projection_coefficient[i_of_sub_qpoint] +=
-              (basis[i_of_sub_qpoint][i_of_basis].transpose().conjugate() * mode_data.get())
+              (basis[i_of_sub_qpoint][i_of_basis].transpose().conjugate() * mode_data.get()).array().abs2().sum();
-                .array().abs2().sum();
         // 如果是严格地将向量分解到一组完备的基矢上, 那么不需要对计算得到的权重再做归一化处理
         // 但这里并不是这样一个严格的概念. 因此对分解到各个 sub qpoint 上的权重做归一化处理
         auto sum = ranges::accumulate(projection_coefficient, 0.);
@@ -184,162 +188,81 @@ void ufo::unfold(std::string config_file)
         return projection_coefficient;
       }
     );
-    // 将计算得到的投影系数重新组装成三维数组
+    return projection_coefficient;
-    // 第一维是 meta qpoint，第二维是模式，第三维是 sub qpoint
-    std::vector<std::vector<std::vector<double>>> projection_coefficient_output;
-    for
-    (
-      std::size_t i_of_meta_qpoint = 0, num_of_mode_manipulated = 0;
-      i_of_meta_qpoint < qpoint_data.size();
-      i_of_meta_qpoint++, num_of_mode_manipulated += qpoint_data[i_of_meta_qpoint].ModeData.size()
-    )
-      projection_coefficient_output.emplace_back
-      (
-        projection_coefficient.begin() + num_of_mode_manipulated,
-        projection_coefficient.begin() + num_of_mode_manipulated + qpoint_data[i_of_meta_qpoint].ModeData.size()
-      );
-    return projection_coefficient_output;
-  };
-
-  // 组装输出，即将投影系数应用到原始数据上
-  auto construct_output = []
-  (
-    const Input& input,
-    const std::vector<std::vector<std::vector<double>>>& projection_coefficient,
-    const std::vector<UnfoldOutput::MetaQpointDataType>& qpoint_data,
-    const std::optional<std::vector<std::size_t>>& selected_atoms
-  )
-  {
-    UnfoldOutput output;
-    output.PrimativeCell = input.PrimativeCell;
-    output.SuperCellMultiplier = input.SuperCellMultiplier;
-    output.SuperCellDeformation = input.SuperCellDeformation;
-    output.SelectedAtoms = selected_atoms;
-    output.MetaQpointData = qpoint_data;
-    for (std::size_t i_of_meta_qpoint = 0; i_of_meta_qpoint < qpoint_data.size(); i_of_meta_qpoint++)
-    {
-      // 如果需要投影到特定的原子上，需要先计算当前 meta qpoint 的不同模式的投影系数
-      std::optional<std::vector<double>> projection_coefficient_on_atoms;
-      if (selected_atoms)
-      {
-        projection_coefficient_on_atoms.emplace();
-        for (std::size_t i_of_mode = 0; i_of_mode < qpoint_data[i_of_meta_qpoint].ModeData.size(); i_of_mode++)
-        {
-          projection_coefficient_on_atoms->emplace_back(0);
-          for (auto atom : *selected_atoms)
-            projection_coefficient_on_atoms->back()
-              += qpoint_data[i_of_meta_qpoint].ModeData[i_of_mode].AtomMovement.row(atom).array().abs2().sum();
-          projection_coefficient_on_atoms->back() *=
-            static_cast<double>(qpoint_data[i_of_meta_qpoint].ModeData[i_of_mode].AtomMovement.rows())
-              / selected_atoms->size();
-        }
-      }
-
-      for
-      (
-        auto [diff_of_sub_qpoint_by_reciprocal_modified_super_cell, i_of_sub_qpoint]
-          : biu::sequence(input.SuperCellMultiplier)
-      )
-      {
-        auto& _ = output.QpointData.emplace_back();
-        /*
-          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 =
-        (
-          input.SuperCellMultiplier.cast<double>().cwiseInverse().asDiagonal()
-          * (
-            diff_of_sub_qpoint_by_reciprocal_modified_super_cell.cast<double>()
-              + input.SuperCellDeformation.inverse() * qpoint_data[i_of_meta_qpoint].Qpoint
-          )
-        ).eval();
-        _.Qpoint = sub_qpoint_by_reciprocal_primative_cell.array()
-          - sub_qpoint_by_reciprocal_primative_cell.array().floor();
-        _.Source = qpoint_data[i_of_meta_qpoint].Qpoint;
-        _.SourceIndex = i_of_meta_qpoint;
-
-        for (std::size_t i_of_mode = 0; i_of_mode < qpoint_data[i_of_meta_qpoint].ModeData.size(); i_of_mode++)
-        {
-          auto& __ = _.ModeData.emplace_back();
-          __.Frequency = qpoint_data[i_of_meta_qpoint].ModeData[i_of_mode].Frequency;
-          __.Weight = projection_coefficient[i_of_meta_qpoint][i_of_mode][i_of_sub_qpoint];
-          if (selected_atoms)
-            __.Weight *= projection_coefficient_on_atoms.value()[i_of_mode];
-        }
-      }
-    }
-    return output;
   };
 
   biu::Logger::Guard log;
-  log.info("Reading input file... ");
+
-  auto input = YAML::LoadFile(config_file).as<Input>();
+  log.info("Reading input file...");
-  auto qpoint_data = read_qpoint_data(input.QpointDataInputFile.value_or("band.hdf5"));
+  auto config = YAML::LoadFile(config_file).as<Config>();
+  CommonData output;
+  output.AtomMass = read_qpoint
+    (config.PrimativePhonopy, config.PrimativeQpoint, output.Primative);
+  output.AtomMass.merge(read_qpoint
+    (config.SuperPhonopy, config.SuperQpoint, output.Super));
+  output.Super.CellDeformation = config.SuperCellDeformation;
+  output.Super.CellMultiplier = config.SuperCellMultiplier;
   log.info("Done.");
 
-  std::clog << "Constructing basis... " << std::flush;
+  log.info("Constructing basis...");
-
   auto basis = construct_basis
   (
-    input.PrimativeCell, input.SuperCellMultiplier,
+    output.Primative.Cell, output.Super.CellMultiplier,
-    input.PrimativeCellBasisNumber,
+    config.PrimativeCellBasisNumber,
-    input.AtomPositionBySuperCell
+    output.Super.AtomPosition
-      * (input.SuperCellDeformation * input.SuperCellMultiplier.cast<double>().asDiagonal() * input.PrimativeCell)
+      * (output.Super.CellDeformation * output.Super.CellMultiplier.cast<double>().asDiagonal() * output.Primative.Cell)
   );
-  std::clog << "Done." << std::endl;
+  log.info("Done.");
 
   std::clog << "Calculating projection coefficient... " << std::flush;
-  // 用来在屏幕上输出进度的计数器和线程
+  // 将所有模式放到一维来处理
-  std::atomic<std::size_t> number_of_finished_modes(0);
-  auto number_of_modes = ranges::accumulate
-  (
-    qpoint_data
-      | ranges::views::transform([](const auto& qpoint)
-        { return qpoint.ModeData.size(); }),
-    0ul
-  );
-  std::atomic<bool> finished;
-  std::thread print_thread([&]
   {
-    while (true)
+    auto modes = output.Super.Qpoint
+      | ranges::views::transform([](const auto& qpoint)
+        { return qpoint.Mode | ranges::views::transform([](const auto& mode)
+          { return std::cref(mode.EigenVector); }); })
+      | ranges::views::join
+      | ranges::to_vector;
+    std::atomic<std::size_t> number_of_finished_modes(0);
+    std::thread print_thread([&]
     {
-      std::osyncstream(std::clog)
+      while (true)
-        << "\rCalculating projection coefficient... ({}/{})"_f(number_of_finished_modes, number_of_modes)
+      {
-        << std::flush;
+        std::osyncstream(std::clog)
-      std::this_thread::sleep_for(100ms);
+          << "\rCalculating projection coefficient... ({}/{})"_f(number_of_finished_modes, modes.size())
-      if (finished) break;
+          << std::flush;
-    }
+        std::this_thread::sleep_for(100ms);
-  });
+        if (number_of_finished_modes == modes.size()) break;
-  auto projection_coefficient = construct_projection_coefficient(basis, qpoint_data, number_of_finished_modes);
+      }
-  finished = true;
+    });
-  print_thread.join();
+    auto projection_coefficient = construct_projection_coefficient
+      (basis, modes, number_of_finished_modes);
+    std::size_t i_of_modes = 0;
+    for (auto& qpoint : output.Super.Qpoint) for (auto& mode : qpoint.Mode)
+      mode.WeightOnUnfold = projection_coefficient[i_of_modes++];
+    print_thread.join();
+  }
   std::clog << "\33[2K\rCalculating projection coefficient... Done." << std::endl;
 
-  std::clog << "Writing data... " << std::flush;
+  log.info("Writing data... ");
-  for (auto& output_file : input.QpointDataOutputFile)
+  std::ofstream(config.OutputFile, std::ios::binary) << biu::serialize<char>(output);
-  {
+  log.info("Done.");
-    auto output = construct_output
+
-      (input, projection_coefficient, qpoint_data, output_file.SelectedAtoms);
+  log.info("Summary:");
-    if (output_file.OutputAsYaml.value_or(false)) std::ofstream(output_file.Filename) << YAML::Node(output);
+
-    else std::ofstream(output_file.Filename, std::ios::binary) << biu::serialize<char>(output);
+  for (auto i_of_qpoint : std::views::iota(0u, output.Super.Qpoint.size()))
-  }
+    for (auto i_of_mode : std::views::iota(0u, output.Super.Qpoint[i_of_qpoint].Mode.size()))
-  std::clog << "Done." << std::endl;
+      for
+      (
+        auto i_of_sub_qpoint
+          : std::views::iota(0u, output.Super.Qpoint[i_of_qpoint].SubQpoint.size())
+      )
+      if (output.Super.Qpoint[i_of_qpoint].Mode[i_of_mode].WeightOnUnfold[i_of_sub_qpoint] > 0.01)
+        log.info("{}:{:.2f}:{}:{:.2f} -> {:.2f} {:.2f}"_f
+        (
+          i_of_qpoint, fmt::join(output.Super.Qpoint[i_of_qpoint].Qpoint, ", "),
+          i_of_mode, output.Super.Qpoint[i_of_qpoint].Mode[i_of_mode].Frequency,
+          fmt::join(output.Super.Qpoint[i_of_qpoint].SubQpoint[i_of_sub_qpoint], ", "),
+          output.Super.Qpoint[i_of_qpoint].Mode[i_of_mode].WeightOnUnfold[i_of_sub_qpoint]
+        ));
 }