localPackages.ufo: init

.gitignore

@@ -1,5 +1,5 @@
 result
-result-man
+result-*
 outputs
 .direnv
 build

flake.nix

@@ -100,6 +100,7 @@
               duc = pkgs.pkgsStatic.duc.override { enableCairo = false; cairo = null; pango = null; };
             in pkgs.pkgsStatic.localPackages.hpcstat.override
               { inherit openssh duc; standalone = true; version = inputs.self.rev or "dirty"; };
+          ufo = pkgs.pkgsStatic.localPackages.ufo.override { version = inputs.self.rev or "dirty"; };
           nixpkgs = pkgs;
         }
         // (
@@ -200,6 +201,12 @@
             buildInputs = [ pkgs.clang-tools_18 ];
             CMAKE_EXPORT_COMPILE_COMMANDS = "1";
           };
+          ufo = pkgs.mkShell
+          {
+            inputsFrom = [ (inputs.self.packages.x86_64-linux.ufo.override { version = null; }) ];
+            packages = [ pkgs.clang-tools_18 ];
+            CMAKE_EXPORT_COMPILE_COMMANDS = "1";
+          };
         };
     };
 }

local/pkgs/default.nix

@@ -79,6 +79,8 @@ inputs: rec
   mkPnpmPackage = inputs.pkgs.callPackage ./mkPnpmPackage.nix {};
   nodejs-with-pnpm9 = inputs.pkgs.callPackage ./nodejs-with-pnpm9.nix {};
   sbatch-tui = inputs.pkgs.callPackage ./sbatch-tui { inherit biu; };
+  ufo = inputs.pkgs.callPackage ./ufo
+    { inherit concurrencpp biu glad matplotplusplus zpp-bits; tbb = inputs.pkgs.tbb_2021_11; };
 
   fromYaml = content: builtins.fromJSON (builtins.readFile
     (inputs.pkgs.runCommand "toJSON" {}

local/pkgs/ufo/.gitignore

@@ -0,0 +1 @@
+test

local/pkgs/ufo/CMakeLists.txt

@@ -0,0 +1,33 @@
+cmake_minimum_required(VERSION 3.14)
+project(ufo VERSION 0 LANGUAGES CXX)
+enable_testing()
+include(GNUInstallDirs)
+
+if(NOT CMAKE_BUILD_TYPE AND NOT CMAKE_CONFIGURATION_TYPES)
+  message("Setting build type to 'Release' as none was specified.")
+  set(CMAKE_BUILD_TYPE Release CACHE STRING "Choose the type of build." FORCE)
+  set_property(CACHE CMAKE_BUILD_TYPE PROPERTY STRINGS "Debug" "Release" "MinSizeRel" "RelWithDebInfo")
+endif()
+
+find_package(yaml-cpp REQUIRED)
+find_package(Eigen3 REQUIRED)
+find_package(fmt REQUIRED)
+find_package(concurrencpp REQUIRED)
+find_package(HighFive REQUIRED)
+find_package(TBB REQUIRED)
+find_package(glad REQUIRED)
+find_package(Matplot++ REQUIRED)
+find_path(ZPP_BITS_INCLUDE_DIR zpp_bits.h REQUIRED)
+
+add_executable(ufo src/solver.cpp src/fold.cpp src/unfold.cpp src/plot.cpp src/main.cpp)
+target_include_directories(ufo PRIVATE ${PROJECT_SOURCE_DIR}/include ${ZPP_BITS_INCLUDE_DIR})
+target_link_libraries(ufo PRIVATE
+  yaml-cpp Eigen3::Eigen fmt::fmt concurrencpp::concurrencpp HighFive_HighFive TBB::tbb Matplot++::matplot
+  Matplot++::matplot_opengl)
+target_compile_features(ufo PRIVATE cxx_std_23)
+
+install(TARGETS ufo RUNTIME DESTINATION ${CMAKE_INSTALL_BINDIR})
+
+get_property(ImportedTargets DIRECTORY "${CMAKE_SOURCE_DIR}" PROPERTY IMPORTED_TARGETS)
+message("Imported targets: ${ImportedTargets}")
+message("List of compile features: ${CMAKE_CXX_COMPILE_FEATURES}")

local/pkgs/ufo/default.nix

@@ -0,0 +1,10 @@
+{
+  stdenv, cmake, pkg-config,
+  yaml-cpp, eigen, fmt, concurrencpp, highfive, tbb, glad, matplotplusplus, biu, zpp-bits
+}: stdenv.mkDerivation
+{
+  name = "ufo";
+  src = ./.;
+  buildInputs = [ yaml-cpp eigen fmt concurrencpp highfive tbb glad matplotplusplus biu zpp-bits ];
+  nativeBuildInputs = [ cmake pkg-config ];
+}

local/pkgs/ufo/include/ufo/fold.hpp

@@ -0,0 +1,37 @@
+# pragma once
+# include <ufo/solver.hpp>
+
+namespace ufo
+{
+  class FoldSolver : public Solver
+  {
+    public:
+      struct InputType
+      {
+        Eigen::Vector<unsigned, 3> SuperCellMultiplier;
+        std::optional<Eigen::Matrix<double, 3, 3>> SuperCellDeformation;
+        std::vector<Eigen::Vector3d> Qpoints;
+        DataFile OutputFile;
+
+        InputType(std::string config_file);
+      };
+      struct OutputType
+      {
+        std::vector<Eigen::Vector3d> Qpoints;
+        void write(std::string filename) const;
+      };
+    protected:
+      InputType Input_;
+      std::optional<OutputType> Output_;
+    public:
+      FoldSolver(std::string config_file);
+      FoldSolver& operator()() override;
+      // return value: QpointInReciprocalSuperCellByReciprocalSuperCell
+      static Eigen::Vector3d fold
+      (
+        Eigen::Vector3d qpoint_in_reciprocal_primitive_cell_by_reciprocal_primitive_cell,
+        Eigen::Vector<unsigned, 3> super_cell_multiplier,
+        std::optional<Eigen::Matrix<double, 3, 3>> super_cell_deformation
+      );
+  };
+}

local/pkgs/ufo/include/ufo/plot.hpp

@@ -0,0 +1,77 @@
+# pragma once
+# include <ufo/unfold.hpp>
+
+namespace ufo
+{
+  class PlotSolver : public Solver
+  {
+    public:
+      struct InputType
+      {
+        Eigen::Matrix3d PrimativeCell;
+
+        struct FigureConfigType
+        {
+          std::vector<std::vector<Eigen::Vector3d>> Qpoints;
+          std::pair<unsigned, unsigned> Resolution;
+          std::pair<double, double> Range;
+          std::optional<std::vector<double>> YTicks;
+          DataFile PictureFile;
+          std::optional<std::vector<DataFile>> DataFiles;
+        };
+        std::vector<FigureConfigType> Figures;
+
+        struct UnfoldedDataType : public UnfoldSolver::OutputType
+        {
+          UnfoldedDataType(std::string filename);
+          UnfoldedDataType() = default;
+        };
+        DataFile UnfoldedDataFile;
+        UnfoldedDataType UnfoldedData;
+
+        InputType(std::string config_file);
+      };
+      struct OutputType
+      {
+        std::vector<std::vector<double>> Values;
+        std::vector<double> XTicks;
+        std::vector<double> YTicks;
+        std::pair<unsigned, unsigned> Resolution;
+        std::pair<double, double> Range;
+
+        OutputType() = default;
+        const OutputType& write(std::string filename, std::string format) const;
+        using serialize = zpp::bits::members<5>;
+      };
+    protected:
+      InputType Input_;
+      std::optional<std::vector<OutputType>> Output_;
+    public:
+      PlotSolver(std::string config_file);
+      PlotSolver& operator()() override;
+
+      // 根据 q 点路径, 搜索要使用的 q 点
+      static std::vector<std::reference_wrapper<const UnfoldSolver::OutputType::QpointDataType>> search_qpoints
+      (
+        const std::pair<Eigen::Vector3d, Eigen::Vector3d>& path,
+        const decltype(InputType::UnfoldedDataType::QpointData)& available_qpoints,
+        double threshold, bool exclude_endpoint = false
+      );
+      // 根据搜索到的 q 点, 计算每个点的数值
+      static std::tuple<std::vector<std::vector<double>>, std::vector<double>> calculate_values
+      (
+        const Eigen::Matrix3d primative_cell,
+        const std::vector<std::pair<Eigen::Vector3d, Eigen::Vector3d>>& path,
+        const std::vector<std::vector<std::reference_wrapper<const UnfoldSolver::OutputType::QpointDataType>>>& qpoints,
+        const decltype(InputType::FigureConfigType::Resolution)& resolution,
+        const decltype(InputType::FigureConfigType::Range)& range
+      );
+      // 根据数值, 画图
+      static void plot
+      (
+        const std::vector<std::vector<double>>& values,
+        const std::string& filename,
+        const std::vector<double>& x_ticks, const std::vector<double>& y_ticks
+      );
+  };
+}

local/pkgs/ufo/include/ufo/solver.hpp

@@ -0,0 +1,141 @@
+# pragma once
+# include <iostream>
+# include <array>
+# include <numbers>
+# include <numeric>
+# include <fstream>
+# include <optional>
+# include <array>
+# include <utility>
+# include <execution>
+# include <syncstream>
+# include <any>
+# include <map>
+# include <vector>
+# include <span>
+# include <yaml-cpp/yaml.h>
+# include <Eigen/Dense>
+# include <concurrencpp/concurrencpp.h>
+# include <fmt/format.h>
+# include <fmt/std.h>
+# include <fmt/ranges.h>
+# include <highfive/H5File.hpp>
+# include <zpp_bits.h>
+# include <matplot/matplot.h>
+# include <matplot/backend/opengl.h>
+
+// 在相位中, 约定为使用 $\exp (2 \pi i \vec{q} \cdot \vec{r})$ 来表示原子的运动状态
+//  (而不是 $\exp (-2 \pi i \vec{q} \cdot \vec{r})$)
+// 一些书定义的倒格矢中包含了 $2 \pi$ 的部分, 我们这里约定不包含这部分.
+//  也就是说, 正格子与倒格子的转置相乘, 得到单位矩阵.
+
+namespace Eigen
+{
+  constexpr inline auto serialize(auto & archive, Eigen::Matrix3d& matrix)
+    { return archive(std::span(matrix.data(), matrix.size())); }
+  constexpr inline auto serialize(auto & archive, const Eigen::Matrix3d& matrix)
+    { return archive(std::span(matrix.data(), matrix.size())); }
+  constexpr inline auto serialize(auto & archive, Eigen::Vector3d& vector)
+    { return archive(std::span(vector.data(), vector.size())); }
+  constexpr inline auto serialize(auto & archive, const Eigen::Vector3d& vector)
+    { return archive(std::span(vector.data(), vector.size())); }
+}
+
+namespace ufo
+{
+  using namespace std::literals;
+  struct PhonopyComplex { double r, i; };
+  inline HighFive::CompoundType create_compound_complex()
+    { return {{ "r", HighFive::AtomicType<double>{}}, {"i", HighFive::AtomicType<double>{}}}; }
+
+	namespace detail_
+	{
+		template <typename T> struct SpecializationOfBitsMembersHelper : std::false_type {};
+		template <std::size_t N> struct SpecializationOfBitsMembersHelper<zpp::bits::members<N>> : std::true_type {};
+	}
+	template <typename T> concept ZppSerializable
+    = requires() { detail_::SpecializationOfBitsMembersHelper<T>::value == true; };
+
+  class Solver
+  {
+    public:
+      virtual Solver& operator()() = 0;
+      virtual ~Solver() = default;
+
+      static concurrencpp::generator<std::pair<Eigen::Vector<unsigned, 3>, unsigned>>
+        triplet_sequence(Eigen::Vector<unsigned, 3> range);
+
+      template <ZppSerializable T> inline static void zpp_write(const T& object, std::string filename)
+      {
+        auto [data, out] = zpp::bits::data_out();
+        out(object).or_throw();
+        static_assert(sizeof(char) == sizeof(std::byte));
+        std::ofstream file(filename, std::ios::binary | std::ios::out);
+        file.exceptions(std::ios::badbit | std::ios::failbit);
+        file.write(reinterpret_cast<const char*>(data.data()), data.size());
+      }
+      template <ZppSerializable T> inline static T zpp_read(std::string filename)
+      {
+        auto input = std::ifstream(filename, std::ios::binary | std::ios::in);
+        input.exceptions(std::ios::badbit | std::ios::failbit);
+        static_assert(sizeof(std::byte) == sizeof(char));
+        std::vector<std::byte> data;
+        {
+          std::vector<char> string(std::istreambuf_iterator<char>(input), {});
+          data.assign
+          (
+            reinterpret_cast<std::byte*>(string.data()),
+            reinterpret_cast<std::byte*>(string.data() + string.size())
+          );
+        }
+        auto in = zpp::bits::in(data);
+        T output;
+        in(output).or_throw();
+        return output;
+      }
+
+      class Hdf5file
+      {
+        public:
+          inline Hdf5file& open_for_read(std::string filename)
+          {
+            File_ = HighFive::File(filename, HighFive::File::ReadOnly);
+            return *this;
+          }
+          inline Hdf5file& open_for_write(std::string filename)
+          {
+            File_ = HighFive::File(filename, HighFive::File::ReadWrite | HighFive::File::Create
+              | HighFive::File::Truncate);
+            return *this;
+          }
+          template <typename T> inline Hdf5file& read(T& object, std::string name)
+          {
+            object = File_->getDataSet(name).read<std::remove_cvref_t<decltype(object)>>();
+            return *this;
+          }
+          template <typename T> inline Hdf5file& write(const T& object, std::string name)
+          {
+            File_->createDataSet(name, object);
+            return *this;
+          }
+        protected:
+          std::optional<HighFive::File> File_;
+      };
+
+      struct DataFile
+      {
+        std::string Filename;
+        std::string Format;
+        std::map<std::string, std::any> ExtraParameters;
+        inline DataFile() = default;
+        DataFile
+        (
+          YAML::Node node, std::set<std::string> supported_format,
+          std::string config_file, bool allow_same_as_config_file = false
+        );
+      };
+
+  };
+}
+
+HIGHFIVE_REGISTER_TYPE(ufo::PhonopyComplex, ufo::create_compound_complex)

local/pkgs/ufo/include/ufo/unfold.hpp

@@ -0,0 +1,164 @@
+# pragma once
+# include <ufo/solver.hpp>
+
+namespace ufo
+{
+  // 反折叠的原理: 将超胞中的原子运动状态, 投影到一组平面波构成的基矢中.
+  // 每一个平面波的波矢由两部分相加得到: 一部分是单胞倒格子的整数倍, 所取的个数有一定任意性, 论文中建议取大约单胞中原子个数那么多个;
+  //  对于没有缺陷的情况, 取一个应该就足够了.
+  // 另一部分是超胞倒格子的整数倍, 取 n 个, n 为超胞对应的单胞的倍数, 其实也就是倒空间中单胞对应倒格子中超胞的格点.
+  // 只要第一部分取得足够多, 那么单胞中原子的状态就可以完全被这些平面波描述.
+  // 将超胞中原子的运动状态投影到这些基矢上, 计算出投影的系数, 就可以将超胞的原子运动状态分解到单胞中的多个 q 点上.
+  class UnfoldSolver : public Solver
+  {
+    public:
+      struct InputType
+      {
+        // 单胞的三个格矢，每行表示一个格矢的坐标，单位为埃
+        Eigen::Matrix3d PrimativeCell;
+        // 单胞到超胞的格矢转换时用到的矩阵
+        // SuperCellMultiplier 是一个三维列向量且各个元素都是整数，表示单胞在各个方向扩大到多少倍之后，可以得到和超胞一样的体积
+        // SuperCsolver.hpp>mation 是一个行列式为 1 的矩阵，它表示经过 SuperCellMultiplier 扩大后，还需要怎样的变换才能得到超胞
+        // SuperCell = (SuperCellDeformation * SuperCellMultiplier.asDiagonal()) * PrimativeCell
+        // ReciprocalPrimativeCell = (SuperCellDeformation * SuperCellMultiplier.asDiagonal()).transpose()
+        //  * ReciprocalSuperCell
+        // Position = PositionToCell(line vector) * Cell
+        // InversePosition = InversePositionToCell(line vector) * ReciprocalCell
+        // PositionToSuperCell(line vector) * SuperCell = PositionToPrimativeCell(line vector) * PrimativeCell
+        // ReciprocalPositionToSuperCell(line vector) * ReciprocalSuperCell
+        //  = ReciprocalPositionToPrimativeCell(line vector) * ReciprocalPrimativeCell
+        Eigen::Vector<unsigned, 3> SuperCellMultiplier;
+        std::optional<Eigen::Matrix<double, 3, 3>> SuperCellDeformation;
+        // 在单胞内取几个平面波的基矢
+        Eigen::Vector<unsigned, 3> PrimativeCellBasisNumber;
+
+        // 从哪个文件读入 AtomPosition, 以及这个文件的格式, 格式可选值包括 "yaml"
+        DataFile AtomPositionInputFile;
+        // 从哪个文件读入 QpointData, 以及这个文件的格式, 格式可选值包括 "yaml" 和 "hdf5"
+        DataFile QpointDataInputFile;
+
+        // 超胞中原子的坐标，每行表示一个原子的坐标，单位为埃
+        Eigen::MatrixX3d AtomPosition;
+        // 关于各个 Q 点的数据
+        struct QpointDataType
+        {
+          // Q 点的坐标，单位为超胞的倒格矢
+          Eigen::Vector3d Qpoint;
+
+          // 关于这个 Q 点上各个模式的数据
+          struct ModeDataType
+          {
+            // 模式的频率，单位为 THz
+            double Frequency;
+            // 模式中各个原子的运动状态
+            // 这个数据是这样得到的: phonopy 输出的动态矩阵的 eigenvector 乘以 $\exp(-2 \pi i \vec q \cdot \vec r)$
+            // 这个数据可以认为是原子位移中, 关于超胞有周期性的那一部分, 再乘以原子质量的开方.
+            // 这个数据在读入后会被立即归一化.
+            Eigen::MatrixX3cd AtomMovement;
+          };
+          std::vector<ModeDataType> ModeData;
+        };
+        std::vector<QpointDataType> QpointData;
+
+        // 输出到哪些文件, 以及使用怎样的格式, 格式可选值包括:
+        // yaml: 使用 yaml 格式输出
+        // yaml-human-readable: 使用 yaml 格式输出, 但是输出的结果更适合人类阅读,
+        //  包括合并相近的模式, 去除权重过小的模式, 限制输出的小数位数.
+        // zpp: 使用 zpp-bits 序列化, 可以直接被 plot.cpp 读取
+        std::vector<DataFile> QpointDataOutputFile;
+
+        // 从文件中读取输入 (包括一个较小的配置文件, 和一个 hdf5 或者一个 yaml 文件), 文件中应当包含:
+        // 单胞的格矢: PrimativeCell 单位为埃 直接从 phonopy 的输出中复制
+        // 超胞的倍数: SuperCellMultiplier 手动输入, 为一个包含三个整数的数组
+        // 超胞的变形: SuperCellDeformation 手动输入, 为一个三阶方阵
+        // 平面波的基矢个数: PrimativeCellBasisNumber 手动输入, 为一个包含三个整数的数组
+        // 另外还有一个文件, 直接将 phonopy 的输出复制过来即可, 如果是 yaml, 应该包含下面的内容:
+        // 超胞中原子的坐标: points[*].coordinates 单位为超胞的格矢 直接从 phonopy 的输出中复制
+        // 各个 Q 点的坐标: phonon[*].q-position 单位为超胞的倒格子的格矢 直接从 phonopy 的输出中复制
+        // 各个模式的频率: phonon[*].band[*].frequency 单位为 THz 直接从 phonopy 的输出中复制
+        // 各个模式的原子运动状态: phonon[*].band[*].eigenvector 直接从 phonopy 的输出中复制
+        // 文件中可以有多余的项目, 多余的项目不管.
+        InputType(std::string filename);
+      };
+      struct OutputType
+      {
+        // 关于各个 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;
+
+        void write(decltype(InputType::QpointDataOutputFile) output_files) const;
+        void write(std::string filename, std::string format, unsigned percision = 10) const;
+
+        using serialize = zpp::bits::members<1>;
+
+        virtual ~OutputType() = default;
+      };
+
+      // 第一层是不同的 sub qpoint, 第二层是单胞内不同的平面波
+      using BasisType = std::vector<std::vector<Eigen::VectorXcd>>;
+    protected:
+      InputType Input_;
+      std::optional<OutputType> Output_;
+      std::optional<BasisType> Basis_;
+
+      // 第一层是不同的模式, 第二层是不同的 sub qpoint
+      using ProjectionCoefficientType_ = std::vector<std::vector<double>>;
+
+    public:
+      UnfoldSolver(std::string config_file);
+      UnfoldSolver& operator()() override;
+
+      // 构建基
+      // 每个 q 点对应的一组 sub qpoint。不同的 q 点所对应的 sub qpoint 是不一样的，但 sub qpoint 与 q 点的相对位置一致。
+      // 这里 xyz_of_diff_of_sub_qpoint 即表示这个相对位置。
+      // 由于基只与这个相对位置有关（也就是说，不同 q 点的基是一样的），因此可以先计算出所有的基，这样降低计算量。
+      // 外层下标对应超胞倒格子的整数倍那部分(第二部分), 也就是不同的 sub qpoint
+      // 内层下标对应单胞倒格子的整数倍那部分(第一部分), 也就是 sub qpoint 上的不同平面波（取的数量越多，结果越精确）
+      static BasisType construct_basis
+      (
+        const decltype(InputType::PrimativeCell)& primative_cell,
+        const decltype(InputType::SuperCellMultiplier)& super_cell_multiplier,
+        const decltype(InputType::PrimativeCellBasisNumber)&
+          primative_cell_basis_number,
+        const decltype(InputType::AtomPosition)& atom_position
+      );
+
+      // 计算投影系数, 是反折叠的核心步骤
+      ProjectionCoefficientType_ construct_projection_coefficient
+      (
+        const BasisType& basis,
+        const std::vector<std::reference_wrapper<const decltype
+          (InputType::QpointDataType::ModeDataType::AtomMovement)>>& mode_data,
+        std::atomic<unsigned>& number_of_finished_modes
+      );
+
+      OutputType construct_output
+      (
+        const decltype(InputType::SuperCellMultiplier)& super_cell_multiplier,
+        const decltype(InputType::SuperCellDeformation)& super_cell_deformation,
+        const std::vector<std::reference_wrapper<const decltype
+          (InputType::QpointDataType::Qpoint)>>& meta_qpoint_by_reciprocal_super_cell,
+        const std::vector<std::vector<std::reference_wrapper<const decltype
+          (InputType::QpointDataType::ModeDataType::Frequency)>>>& frequency,
+        const ProjectionCoefficientType_& projection_coefficient
+      );
+  };
+}

local/pkgs/ufo/src/fold.cpp

@@ -0,0 +1,86 @@
+# include <ufo/fold.hpp>
+
+namespace ufo
+{
+  FoldSolver::InputType::InputType(std::string config_file)
+  {
+    auto input = YAML::LoadFile(config_file);
+    for (unsigned i = 0; i < 3; i++)
+      SuperCellMultiplier(i) = input["SuperCellMultiplier"][i].as<unsigned>();
+    if (input["SuperCellDeformation"])
+    {
+      SuperCellDeformation.emplace();
+      for (unsigned i = 0; i < 3; i++)
+        for (unsigned j = 0; j < 3; j++)
+          (*SuperCellDeformation)(i, j) = input["SuperCellDeformation"][i][j].as<double>();
+    }
+    for (auto& qpoint : input["Qpoints"].as<std::vector<std::vector<double>>>())
+      Qpoints.push_back(Eigen::Vector3d
+        {{qpoint.at(0)}, {qpoint.at(1)}, {qpoint.at(2)}});
+    OutputFile = DataFile(input["OutputFile"], {"yaml"}, config_file);
+  }
+  void FoldSolver::OutputType::write(std::string filename) const
+  {
+    std::ofstream(filename) << [&]
+    {
+      std::stringstream print;
+      print << "Qpoints:\n";
+      for (auto& qpoint : Qpoints)
+        print << fmt::format("  - [ {:.8f}, {:.8f}, {:.8f} ]\n", qpoint(0), qpoint(1), qpoint(2));
+      return print.str();
+    }();
+  }
+
+  FoldSolver::FoldSolver(std::string config_file) : Input_(config_file) {}
+  FoldSolver& FoldSolver::operator()()
+  {
+    if (!Output_)
+    {
+      Output_.emplace();
+      for (auto& qpoint : Input_.Qpoints)
+        Output_->Qpoints.push_back(fold
+        (
+          qpoint, Input_.SuperCellMultiplier,
+          Input_.SuperCellDeformation
+        ));
+    }
+    Output_->write(Input_.OutputFile.Filename);
+    return *this;
+  }
+
+  Eigen::Vector3d FoldSolver::fold
+  (
+    Eigen::Vector3d qpoint_in_reciprocal_primitive_cell_by_reciprocal_primitive_cell,
+    Eigen::Vector<unsigned, 3> super_cell_multiplier,
+    std::optional<Eigen::Matrix<double, 3, 3>> super_cell_deformation
+  )
+  {
+    // 首先需要将 q 点转移到 ModifiedSuperCell 的倒格子中
+    // 将 q 点坐标扩大, 然后取小数部分, 就可以了
+    auto qpoint_by_reciprocal_modified_super_cell = super_cell_multiplier.cast<double>().asDiagonal()
+      * qpoint_in_reciprocal_primitive_cell_by_reciprocal_primitive_cell;
+    auto qpoint_in_reciprocal_modified_super_cell_by_reciprocal_modified_super_cell =
+      (qpoint_by_reciprocal_modified_super_cell.array() - qpoint_by_reciprocal_modified_super_cell.array().floor())
+        .matrix();
+    if (!super_cell_deformation)
+      return qpoint_in_reciprocal_modified_super_cell_by_reciprocal_modified_super_cell;
+    /*
+      对 q 点平移数个 SupreCell, 直到它落在超胞的倒格子中
+      这等价于直接将 q 点坐标用 SuperCell 的倒格子表示, 然后取小数部分.
+      ModifiedSuperCell = SuperCellMultiplier * PrimativeCell
+      SuperCell = SuperCellDeformation * ModifiedSuperCell
+      ReciprocalModifiedSuperCell = ModifiedSuperCell.inverse().transpose()
+      ReciprocalSuperCell = SuperCell.inverse().transpose()
+      Qpoint = QpointByReciprocalModifiedSuperCell.transpose() * ReciprocalModifiedSuperCell
+      Qpoint = QpointByReciprocalSuperCell.transpose() * ReciprocalSuperCell
+      整理可以得到:
+      QpointByReciprocalSuperCell = SuperCellDeformation * QpointByReciprocalModifiedSuperCell
+    */
+    auto qpoint_in_reciprocal_modified_super_cell_by_reciprocal_super_cell =
+      (*super_cell_deformation * qpoint_in_reciprocal_modified_super_cell_by_reciprocal_modified_super_cell).eval();
+    auto qpoint_in_reciprocal_super_cell_by_reciprocal_super_cell =
+      qpoint_in_reciprocal_modified_super_cell_by_reciprocal_super_cell.array()
+      - qpoint_in_reciprocal_modified_super_cell_by_reciprocal_super_cell.array().floor();
+    return qpoint_in_reciprocal_super_cell_by_reciprocal_super_cell;
+  }
+}

local/pkgs/ufo/src/main.cpp

@@ -0,0 +1,17 @@
+# include <ufo/fold.hpp>
+# include <ufo/unfold.hpp>
+# include <ufo/plot.hpp>
+
+int main(int argc, const char** argv)
+{
+  if (argc != 3)
+    throw std::runtime_error(fmt::format("Usage: {} task config.yaml", argv[0]));
+  if (argv[1] == std::string("fold"))
+    ufo::FoldSolver{argv[2]}();
+  else if (argv[1] == std::string("unfold"))
+    ufo::UnfoldSolver{argv[2]}();
+  else if (argv[1] == std::string("plot"))
+    ufo::PlotSolver{argv[2]}();
+  else
+    throw std::runtime_error(fmt::format("Unknown task: {}", argv[1]));
+}

local/pkgs/ufo/src/plot.cpp

@@ -0,0 +1,266 @@
+# include <ufo/plot.hpp>
+
+namespace ufo
+{
+  PlotSolver::InputType::UnfoldedDataType::UnfoldedDataType(std::string filename)
+  {
+    static_cast<UnfoldSolver::OutputType&>(*this) = zpp_read<UnfoldSolver::OutputType>(filename);
+  }
+
+  PlotSolver::InputType::InputType(std::string config_file)
+  {
+    auto input = YAML::LoadFile(config_file);
+    for (unsigned i = 0; i < 3; i++)
+      for (unsigned j = 0; j < 3; j++)
+        PrimativeCell(i, j) = input["PrimativeCell"][i][j].as<double>();
+    for (auto& figure : input["Figures"].as<std::vector<YAML::Node>>())
+    {
+      Figures.emplace_back();
+      auto qpoints = figure["Qpoints"]
+        .as<std::vector<std::vector<std::vector<double>>>>();
+      for (auto& line : qpoints)
+      {
+        Figures.back().Qpoints.emplace_back();
+        for (auto& point : line)
+          Figures.back().Qpoints.back().emplace_back(point.at(0), point.at(1), point.at(2));
+        if (Figures.back().Qpoints.back().size() < 2)
+          throw std::runtime_error("Not enough points in a line");
+      }
+      if (Figures.back().Qpoints.size() < 1)
+        throw std::runtime_error("Not enough lines in a figure");
+      Figures.back().Resolution = figure["Resolution"].as<std::pair<unsigned, unsigned>>();
+      Figures.back().Range = figure["Range"].as<std::pair<double, double>>();
+      Figures.back().PictureFile
+        = DataFile(figure["PictureFile"], {"png"}, config_file);
+      if (figure["YTicks"])
+        Figures.back().YTicks = figure["YTicks"].as<std::vector<double>>();
+      if (figure["DataFiles"])
+      {
+        Figures.back().DataFiles.emplace();
+        for (auto& data_file : figure["DataFiles"].as<std::vector<YAML::Node>>())
+          Figures.back().DataFiles->emplace_back()
+            = DataFile(data_file, {"hdf5", "zpp"}, config_file);
+      }
+    }
+    UnfoldedDataFile = DataFile(input["UnfoldedDataFile"], {"zpp"}, config_file);
+    UnfoldedData = UnfoldedDataType(UnfoldedDataFile.Filename);
+  }
+  const PlotSolver::OutputType& PlotSolver::OutputType::write(std::string filename, std::string format) const
+  {
+    if (format == "zpp")
+      zpp_write(*this, filename);
+    else if (format == "hdf5")
+    {
+      std::vector resolution{ Resolution.first, Resolution.second };
+      std::vector range{ Range.first, Range.second };
+      Hdf5file{}.open_for_write(filename).write(Values, "Values")
+        .write(XTicks, "XTicks")
+        .write(YTicks, "YTicks")
+        .write(resolution, "Resolution")
+        .write(range, "Range");
+    }
+    return *this;
+  }
+
+  PlotSolver::PlotSolver(std::string config_file) : Input_(config_file) {}
+
+  PlotSolver& PlotSolver::operator()()
+  {
+    Output_.emplace();
+    for (auto& figure : Input_.Figures)
+    {
+      // 外层表示不同的线段的端点，内层表示这个线段上的 q 点
+      std::vector<std::vector<std::reference_wrapper<const UnfoldSolver::OutputType::QpointDataType>>> qpoints;
+      std::vector<std::pair<Eigen::Vector3d, Eigen::Vector3d>> lines;
+      for (auto& path : figure.Qpoints)
+        for (unsigned i = 0; i < path.size() - 1; i++)
+        {
+          lines.emplace_back(path[i], path[i + 1]);
+          qpoints.push_back(search_qpoints
+          (
+            lines.back(), Input_.UnfoldedData.QpointData,
+            0.001, i != path.size() - 2
+          ));
+        }
+      auto [values, x_ticks] = calculate_values
+      (
+        Input_.PrimativeCell, lines, qpoints, figure.Resolution, figure.Range
+      );
+      auto y_ticks = figure.YTicks.value_or(std::vector<double>{});
+      for (auto& _ : y_ticks)
+        _ = (_ - figure.Range.first) / (figure.Range.second - figure.Range.first) * figure.Resolution.second;
+      plot(values, figure.PictureFile.Filename, x_ticks, y_ticks);
+      Output_->emplace_back();
+      Output_->back().Values = std::move(values);
+      Output_->back().XTicks = std::move(x_ticks);
+      Output_->back().YTicks = std::move(y_ticks);
+      Output_->back().Resolution = figure.Resolution;
+      Output_->back().Range = figure.Range;
+      if (figure.DataFiles)
+        for (auto& data_file : *figure.DataFiles)
+          Output_->back().write(data_file.Filename, data_file.Format);
+    }
+    return *this;
+  }
+
+  std::vector<std::reference_wrapper<const UnfoldSolver::OutputType::QpointDataType>> PlotSolver::search_qpoints
+  (
+    const std::pair<Eigen::Vector3d, Eigen::Vector3d>& path,
+    const decltype(InputType::UnfoldedDataType::QpointData)& available_qpoints,
+    double threshold, bool exclude_endpoint
+  )
+  {
+    std::multimap<double, std::reference_wrapper<const UnfoldSolver::OutputType::QpointDataType>> selected_qpoints;
+    // 对于 output 中的每一个点, 检查这个点是否在路径上. 如果在, 把它加入到 selected_qpoints 中
+    for (auto& qpoint : available_qpoints)
+    {
+      // 计算三点围成的三角形的面积的两倍
+      auto area = (path.second - path.first).cross(qpoint.Qpoint - path.first).norm();
+      // 计算这个点到前两个点所在直线的距离
+      auto distance = area / (path.second - path.first).norm();
+      // 如果这个点到前两个点所在直线的距离小于阈值, 则认为这个点在路径上
+      if (distance < threshold)
+      {
+        // 计算这个点到前两个点的距离, 两个距离都应该小于两点之间的距离
+        auto distance1 = (qpoint.Qpoint - path.first).norm();
+        auto distance2 = (qpoint.Qpoint - path.second).norm();
+        auto distance3 = (path.second - path.first).norm();
+        if (distance1 < distance3 + threshold && distance2 < distance3 + threshold)
+          // 如果这个点不在终点处, 或者不排除终点, 则加入
+          if (distance2 > threshold || !exclude_endpoint)
+            selected_qpoints.emplace(distance1, std::ref(qpoint));
+      }
+    }
+    // 去除非常接近的点
+    for (auto it = selected_qpoints.begin(); it != selected_qpoints.end();)
+    {
+      auto next = std::next(it);
+      if (next == selected_qpoints.end())
+        break;
+      else if (next->first - it->first < threshold)
+        selected_qpoints.erase(next);
+      else
+        it = next;
+    }
+    if (selected_qpoints.empty())
+      throw std::runtime_error("No q points found");
+    std::vector<std::reference_wrapper<const UnfoldSolver::OutputType::QpointDataType>> result;
+    for (auto& qpoint : selected_qpoints)
+      result.push_back(qpoint.second);
+    return result;
+  }
+
+  std::tuple<std::vector<std::vector<double>>, std::vector<double>> PlotSolver::calculate_values
+  (
+    const Eigen::Matrix3d primative_cell,
+    const std::vector<std::pair<Eigen::Vector3d, Eigen::Vector3d>>& path,
+    const std::vector<std::vector<std::reference_wrapper<const UnfoldSolver::OutputType::QpointDataType>>>& qpoints,
+    const decltype(InputType::FigureConfigType::Resolution)& resolution,
+    const decltype(InputType::FigureConfigType::Range)& range
+  )
+  {
+    // 整理输入
+    std::map<double, std::reference_wrapper<const UnfoldSolver::OutputType::QpointDataType>> qpoints_with_distance;
+    double total_distance = 0;
+    std::vector<double> x_ticks;
+    for (unsigned i = 0; i < path.size(); i++)
+    {
+      for (auto& _ : qpoints[i])
+        qpoints_with_distance.emplace
+        (
+          total_distance
+            + ((_.get().Qpoint - path[i].first).transpose() * primative_cell.inverse().transpose()).norm(),
+          _
+        );
+      total_distance += ((path[i].second - path[i].first).transpose() * primative_cell.inverse().transpose()).norm();
+      if (i != path.size() - 1)
+        x_ticks.push_back(total_distance);
+    }
+    for (auto& _ : x_ticks)
+      _ = _ / total_distance * resolution.first;
+
+    // 插值
+    std::vector<std::vector<double>> values;
+    auto blend = []
+    (
+      const UnfoldSolver::OutputType::QpointDataType& a,
+      const UnfoldSolver::OutputType::QpointDataType& b,
+      double ratio, unsigned resolution, std::pair<double, double> range
+    ) -> std::vector<double>
+    {
+      // 计算插值结果
+      std::vector<double> frequency, weight;
+      for (unsigned i = 0; i < a.ModeData.size(); i++)
+      {
+        frequency.push_back(a.ModeData[i].Frequency * ratio + b.ModeData[i].Frequency * (1 - ratio));
+        weight.push_back(a.ModeData[i].Weight * ratio + b.ModeData[i].Weight * (1 - ratio));
+      }
+      std::vector<double> result(resolution);
+      for (unsigned i = 0; i < frequency.size(); i++)
+      {
+        int index = (frequency[i] - range.first) / (range.second - range.first) * resolution;
+        if (index >= 0 && index < static_cast<int>(resolution))
+          result[index] += weight[i];
+      }
+      return result;
+    };
+    for (unsigned i = 0; i < resolution.first; i++)
+    {
+      auto current_distance = total_distance * i / resolution.first;
+      auto it = qpoints_with_distance.lower_bound(current_distance);
+      if (it == qpoints_with_distance.begin())
+        values.push_back(blend(it->second.get(), it->second.get(), 1, resolution.second, range));
+      else if (it == qpoints_with_distance.end())
+        values.push_back(blend(std::prev(it)->second.get(), std::prev(it)->second.get(), 1, resolution.second,
+          range));
+      else
+        values.push_back(blend
+        (
+          std::prev(it)->second.get(), it->second.get(),
+          (it->first - current_distance) / (it->first - std::prev(it)->first),
+          resolution.second, range)
+        );
+    }
+    return {values, x_ticks};
+  }
+  void PlotSolver::plot
+  (
+    const std::vector<std::vector<double>>& values,
+    const std::string& filename,
+    const std::vector<double>& x_ticks, const std::vector<double>& y_ticks
+  )
+  {
+    std::vector<std::vector<double>>
+      r(values[0].size(), std::vector<double>(values.size(), 0)),
+      g(values[0].size(), std::vector<double>(values.size(), 0)),
+      b(values[0].size(), std::vector<double>(values.size(), 0)),
+      a(values[0].size(), std::vector<double>(values.size(), 0));
+    for (unsigned i = 0; i < values[0].size(); i++)
+      for (unsigned j = 0; j < values.size(); j++)
+      {
+        auto v = values[j][i];
+        if (v < 0.05)
+          v = 0;
+        a[i][j] = v * 100 * 255;
+        if (a[i][j] > 255)
+          a[i][j] = 255;
+        r[i][j] = 255 - v * 2 * 255;
+        if (r[i][j] < 0)
+          r[i][j] = 0;
+        g[i][j] = 255 - v * 2 * 255;
+        if (g[i][j] < 0)
+          g[i][j] = 0;
+        b[i][j] = 255;
+      }
+    auto f = matplot::figure<matplot::backend::gnuplot>(true);
+    auto ax = f->current_axes();
+    auto image = ax->image(std::tie(r, g, b));
+    image->matrix_a(a);
+    ax->y_axis().reverse(false);
+    ax->x_axis().tick_values(x_ticks);
+    ax->x_axis().tick_length(1);
+    ax->y_axis().tick_values(y_ticks);
+    ax->y_axis().tick_length(1);
+    f->save(filename, "png");
+  }
+}

local/pkgs/ufo/src/solver.cpp

@@ -0,0 +1,47 @@
+# include <ufo/solver.hpp>
+
+namespace ufo
+{
+  concurrencpp::generator<std::pair<Eigen::Vector<unsigned, 3>, unsigned>> Solver::triplet_sequence
+    (Eigen::Vector<unsigned, 3> range)
+  {
+    for (unsigned x = 0; x < range[0]; x++)
+      for (unsigned y = 0; y < range[1]; y++)
+        for (unsigned z = 0; z < range[2]; z++)
+          co_yield
+          {
+            Eigen::Vector<unsigned, 3>{{x}, {y}, {z}},
+            x * range[1] * range[2] + y * range[2] + z
+          };
+  }
+
+  Solver::DataFile::DataFile
+    (YAML::Node node, std::set<std::string> supported_format, std::string config_file, bool allow_same_as_config_file)
+  {
+    if (auto _ = node["SameAsConfigFile"])
+    {
+      auto __ = _.as<bool>();
+      if (__ && !allow_same_as_config_file)
+        throw std::runtime_error("\"SameAsConfigFile: true\" is not allowed here.");
+      ExtraParameters["SameAsConfigFile"] = __;
+      if (__)
+      {
+        Filename = config_file;
+        Format = "yaml";
+        return;
+      }
+    }
+    Filename = node["Filename"].as<std::string>();
+    Format = node["Format"].as<std::string>();
+    if (!supported_format.contains(Format))
+      throw std::runtime_error(fmt::format("Unsupported format: \"{}\"", Format));
+    if (auto _ = node["RelativeToConfigFile"])
+    {
+      auto __ = _.as<bool>();
+      ExtraParameters["RelativeToConfigFile"] = __;
+      if (__)
+        Filename = std::filesystem::path(config_file).parent_path() / Filename;
+    }
+  };
+
+}

local/pkgs/ufo/src/unfold.cpp

@@ -0,0 +1,429 @@
+# include <ufo/unfold.hpp>
+
+namespace ufo
+{
+  UnfoldSolver::InputType::InputType(std::string filename)
+  {
+    // read main input file
+    {
+      auto node = YAML::LoadFile(filename);
+        for (unsigned i = 0; i < 3; i++)
+          for (unsigned j = 0; j < 3; j++)
+            PrimativeCell(i, j) = node["PrimativeCell"][i][j].as<double>();
+
+      for (unsigned i = 0; i < 3; i++)
+        SuperCellMultiplier(i) = node["SuperCellMultiplier"][i].as<int>();
+
+      if (auto value = node["SuperCellDeformation"])
+      {
+        SuperCellDeformation.emplace();
+        for (unsigned i = 0; i < 3; i++)
+          for (unsigned j = 0; j < 3; j++)
+            (*SuperCellDeformation)(i, j) = value[i][j].as<double>();
+      }
+
+      for (unsigned i = 0; i < 3; i++)
+        PrimativeCellBasisNumber(i) = node["PrimativeCellBasisNumber"][i].as<int>();
+
+      AtomPositionInputFile = DataFile
+      (
+        node["AtomPositionInputFile"], {"yaml"},
+        filename, true
+      );
+      QpointDataInputFile = DataFile
+      (
+        node["QpointDataInputFile"], {"yaml", "hdf5"},
+        filename, true
+      );
+      if (auto value = node["QpointDataOutputFile"])
+      {
+        QpointDataOutputFile.resize(value.size());
+        for (unsigned i = 0; i < value.size(); i++)
+          QpointDataOutputFile[i] = DataFile
+          (
+            value[i], {"yaml", "yaml-human-readable", "zpp", "hdf5"},
+            filename, false
+          );
+      }
+    }
+
+    if (AtomPositionInputFile.Format == "yaml")
+    {
+      auto node = YAML::LoadFile(AtomPositionInputFile.Filename);
+      std::vector<YAML::Node> points;
+      if (auto _ = node["points"])
+        points = _.as<std::vector<YAML::Node>>();
+      else
+        points = node["unit_cell"]["points"].as<std::vector<YAML::Node>>();
+      auto atom_position_to_super_cell = Eigen::MatrixX3d(points.size(), 3);
+      for (unsigned i = 0; i < points.size(); i++)
+        for (unsigned j = 0; j < 3; j++)
+          atom_position_to_super_cell(i, j) = points[i]["coordinates"][j].as<double>();
+      auto super_cell = (SuperCellDeformation.value_or(Eigen::Matrix3d::Identity())
+        * SuperCellMultiplier.cast<double>().asDiagonal() * PrimativeCell).eval();
+      AtomPosition = atom_position_to_super_cell * super_cell;
+    }
+    if (QpointDataInputFile.Format == "yaml")
+    {
+      auto node = YAML::LoadFile(QpointDataInputFile.Filename);
+      auto phonon = node["phonon"].as<std::vector<YAML::Node>>();
+      QpointData.resize(phonon.size());
+      for (unsigned i = 0; i < phonon.size(); i++)
+      {
+        for (unsigned j = 0; j < 3; j++)
+          QpointData[i].Qpoint(j) = phonon[i]["q-position"][j].as<double>();
+        auto band = phonon[i]["band"].as<std::vector<YAML::Node>>();
+        QpointData[i].ModeData.resize(band.size());
+        for (unsigned j = 0; j < band.size(); j++)
+        {
+          QpointData[i].ModeData[j].Frequency = band[j]["frequency"].as<double>();
+          auto eigenvector_vectors = band[j]["eigenvector"]
+            .as<std::vector<std::vector<std::vector<double>>>>();
+          Eigen::MatrixX3cd eigenvectors(AtomPosition.rows(), 3);
+          for (unsigned k = 0; k < AtomPosition.rows(); k++)
+            for (unsigned l = 0; l < 3; l++)
+              eigenvectors(k, l)
+                = eigenvector_vectors[k][l][0] + 1i * eigenvector_vectors[k][l][1];
+          // 需要对读入的原子运动状态作相位转换, 使得它们与我们的约定一致(对超胞周期性重复)
+          // 这里还要需要做归一化处理 (指将数据简单地作为向量处理的归一化)
+          auto& AtomMovement = QpointData[i].ModeData[j].AtomMovement;
+          // AtomMovement = eigenvectors.array().colwise() * (-2 * std::numbers::pi_v<double> * 1i
+          //   * (atom_position_to_super_cell * input.QpointData[i].Qpoint)).array().exp();
+          // AtomMovement /= AtomMovement.norm();
+          // phonopy 似乎已经进行了相位的转换！为什么？
+          AtomMovement = eigenvectors / eigenvectors.norm();
+        }
+      }
+    }
+    else if (QpointDataInputFile.Format == "hdf5")
+    {
+      std::vector<std::vector<std::vector<double>>> frequency, path;
+      std::vector<std::vector<std::vector<std::vector<PhonopyComplex>>>> eigenvector_vector;
+      Hdf5file{}.open_for_read(QpointDataInputFile.Filename).read(frequency, "/frequency")
+        .read(eigenvector_vector, "/eigenvector")
+        .read(path, "/path");
+      std::vector size = { frequency.size(), frequency[0].size(), frequency[0][0].size() };
+      QpointData.resize(size[0] * size[1]);
+      for (unsigned i = 0; i < size[0]; i++)
+        for (unsigned j = 0; j < size[1]; j++)
+        {
+          QpointData[i * size[1] + j].Qpoint = Eigen::Vector3d(path[i][j].data());
+          QpointData[i * size[1] + j].ModeData.resize(size[2]);
+          for (unsigned k = 0; k < size[2]; k++)
+          {
+            QpointData[i * size[1] + j].ModeData[k].Frequency = frequency[i][j][k];
+            Eigen::MatrixX3cd eigenvectors(AtomPosition.rows(), 3);
+            for (unsigned l = 0; l < AtomPosition.rows(); l++)
+              for (unsigned m = 0; m < 3; m++)
+                eigenvectors(l, m)
+                  = eigenvector_vector[i][j][l * 3 + m][k].r + eigenvector_vector[i][j][l * 3 + m][k].i * 1i;
+            QpointData[i * size[1] + j].ModeData[k].AtomMovement = eigenvectors / eigenvectors.norm();
+          }
+        }
+    }
+  }
+
+  void UnfoldSolver::OutputType::write
+    (decltype(InputType::QpointDataOutputFile) output_files) const
+  {
+    for (auto& output_file : output_files)
+      write(output_file.Filename, output_file.Format);
+  }
+  void UnfoldSolver::OutputType::write(std::string filename, std::string format, unsigned percision) const
+  {
+    if (format == "yaml")
+      std::ofstream(filename) << [&]
+      {
+        std::stringstream print;
+        print << "QpointData:\n";
+        for (auto& qpoint: QpointData)
+        {
+          print << fmt::format("  - Qpoint: [ {1:.{0}f}, {2:.{0}f}, {3:.{0}f} ]\n",
+            percision, qpoint.Qpoint[0], qpoint.Qpoint[1], qpoint.Qpoint[2]);
+          print << fmt::format("    Source: [ {1:.{0}f}, {2:.{0}f}, {3:.{0}f} ]\n",
+            percision, qpoint.Source[0], qpoint.Source[1], qpoint.Source[2]);
+          print << "    ModeData:\n";
+          for (auto& mode: qpoint.ModeData)
+            print << fmt::format("      - {{ Frequency: {1:.{0}f}, Weight: {2:.{0}f} }}\n",
+              percision, mode.Frequency, mode.Weight);
+        }
+        return print.str();
+      }();
+    else if (format == "yaml-human-readable")
+    {
+      std::remove_cvref_t<decltype(*this)> output;
+      std::map<unsigned, std::vector<decltype(QpointData)::const_iterator>>
+        meta_qpoint_to_sub_qpoint_iterators;
+      for (auto it = QpointData.begin(); it != QpointData.end(); it++)
+        meta_qpoint_to_sub_qpoint_iterators[it->SourceIndex_].push_back(it);
+      for (auto [meta_qpoint_index, sub_qpoint_iterators] : meta_qpoint_to_sub_qpoint_iterators)
+        for (auto& qpoint : sub_qpoint_iterators)
+        {
+          std::map<double, double> frequency_to_weight;
+          for (unsigned i_of_mode = 0; i_of_mode < qpoint->ModeData.size(); i_of_mode++)
+          {
+            auto frequency = qpoint->ModeData[i_of_mode].Frequency;
+            auto weight = qpoint->ModeData[i_of_mode].Weight;
+            auto it_lower = frequency_to_weight.lower_bound(frequency - 0.1);
+            auto it_upper = frequency_to_weight.upper_bound(frequency + 0.1);
+            if (it_lower == it_upper)
+              frequency_to_weight[frequency] = weight;
+            else
+            {
+              auto frequency_sum = std::accumulate(it_lower, it_upper, 0.,
+                [](const auto& a, const auto& b) { return a + b.first * b.second; });
+              auto weight_sum = std::accumulate(it_lower, it_upper, 0.,
+                [](const auto& a, const auto& b) { return a + b.second; });
+              frequency_sum += frequency * weight;
+              weight_sum += weight;
+              frequency_to_weight.erase(it_lower, it_upper);
+              frequency_to_weight[frequency_sum / weight_sum] = weight_sum;
+            }
+          }
+          auto& _ = output.QpointData.emplace_back();
+          _.Qpoint = qpoint->Qpoint;
+          _.Source = qpoint->Source;
+          _.SourceIndex_ = qpoint->SourceIndex_;
+          for (auto [frequency, weight] : frequency_to_weight)
+            if (weight > 0.1)
+            {
+              auto& __ = _.ModeData.emplace_back();
+              __.Frequency = frequency;
+              __.Weight = weight;
+            }
+        }
+      output.write(filename, "yaml", 3);
+    }
+    else if (format == "zpp")
+      zpp_write(*this, filename);
+    else if (format == "hdf5")
+    {
+      std::vector<std::vector<double>> Qpoint, Source, Frequency, Weight;
+      for (auto& qpoint : QpointData)
+      {
+        Qpoint.emplace_back(qpoint.Qpoint.data(), qpoint.Qpoint.data() + 3);
+        Source.emplace_back(qpoint.Source.data(), qpoint.Source.data() + 3);
+        Frequency.emplace_back();
+        Weight.emplace_back();
+        for (auto& mode : qpoint.ModeData)
+        {
+          Frequency.back().push_back(mode.Frequency);
+          Weight.back().push_back(mode.Weight);
+        }
+      }
+      Hdf5file{}.open_for_write(filename).write(Qpoint, "/Qpoint")
+        .write(Source, "/Source")
+        .write(Frequency, "/Frequency")
+        .write(Weight, "/Weight");
+    }
+  }
+
+  UnfoldSolver::UnfoldSolver(std::string config_file) : Input_([&]
+  {
+    std::clog << "Reading input file... " << std::flush;
+    return config_file;
+  }())
+  {
+    std::clog << "Done." << std::endl;
+  }
+
+  UnfoldSolver& UnfoldSolver::operator()()
+  {
+    if (!Basis_)
+    {
+      std::clog << "Constructing basis... " << std::flush;
+      Basis_ = construct_basis
+      (
+        Input_.PrimativeCell, Input_.SuperCellMultiplier,
+        Input_.PrimativeCellBasisNumber, Input_.AtomPosition
+      );
+      std::clog << "Done." << std::endl;
+    }
+    if (!Output_)
+    {
+      std::clog << "Calculating projection coefficient... " << std::flush;
+      std::vector<std::reference_wrapper<const decltype
+        (InputType::QpointDataType::ModeDataType::AtomMovement)>> mode_data;
+      for (auto& qpoint : Input_.QpointData)
+        for (auto& mode : qpoint.ModeData)
+          mode_data.emplace_back(mode.AtomMovement);
+      std::atomic<unsigned> number_of_finished_modes(0);
+      std::thread print_thread([&]
+      {
+        unsigned n;
+        while ((n = number_of_finished_modes) < mode_data.size())
+        {
+          std::osyncstream(std::cerr) << fmt::format("\rCalculating projection coefficient... ({}/{})",
+            number_of_finished_modes, mode_data.size()) << std::flush;
+          std::this_thread::sleep_for(100ms);
+          number_of_finished_modes.wait(n);
+        }
+      });
+      auto projection_coefficient = construct_projection_coefficient
+        (*Basis_, mode_data, number_of_finished_modes);
+      number_of_finished_modes = mode_data.size();
+      print_thread.join();
+      std::clog << "\33[2K\rCalculating projection coefficient... Done." << std::endl;
+
+      std::clog << "Constructing output... " << std::flush;
+      std::vector<std::reference_wrapper<const decltype(InputType::QpointDataType::Qpoint)>> qpoint;
+      std::vector<std::vector<std::reference_wrapper<const
+        decltype(InputType::QpointDataType::ModeDataType::Frequency)>>> frequency;
+      for (auto& qpoint_data : Input_.QpointData)
+      {
+        qpoint.emplace_back(qpoint_data.Qpoint);
+        frequency.emplace_back();
+        for (auto& mode_data : qpoint_data.ModeData)
+          frequency.back().emplace_back(mode_data.Frequency);
+      }
+      Output_ = construct_output
+      (
+        Input_.SuperCellMultiplier,
+        Input_.SuperCellDeformation, qpoint, frequency, projection_coefficient
+      );
+      std::clog << "Done." << std::endl;
+    }
+    std::clog << "Writing output... " << std::flush;
+    Output_->write(Input_.QpointDataOutputFile);
+    std::clog << "Done." << std::endl;
+    return *this;
+  }
+
+  UnfoldSolver::BasisType UnfoldSolver::construct_basis
+  (
+    const decltype(InputType::PrimativeCell)& primative_cell,
+    const decltype(InputType::SuperCellMultiplier)& super_cell_multiplier,
+    const decltype(InputType::PrimativeCellBasisNumber)& primative_cell_basis_number,
+    const decltype(InputType::AtomPosition)& atom_position
+  )
+  {
+    BasisType basis(super_cell_multiplier.prod());
+    // 每个 q 点对应的一组 sub qpoint。不同的 q 点所对应的 sub qpoint 是不一样的，但 sub qpoint 与 q 点的相对位置一致。
+    // 这里 xyz_of_diff_of_sub_qpoint 即表示这个相对位置，单位为超胞的倒格矢
+    for (auto [xyz_of_diff_of_sub_qpoint_by_reciprocal_modified_super_cell, i_of_sub_qpoint]
+      : triplet_sequence(super_cell_multiplier))
+    {
+      basis[i_of_sub_qpoint].resize(primative_cell_basis_number.prod());
+      for (auto [xyz_of_basis, i_of_basis] : triplet_sequence(primative_cell_basis_number))
+      {
+        // 计算 q 点的坐标, 单位为单胞的倒格矢
+        auto diff_of_sub_qpoint_by_reciprocal_primative_cell = xyz_of_basis.cast<double>()
+          + super_cell_multiplier.cast<double>().cwiseInverse().asDiagonal()
+          * xyz_of_diff_of_sub_qpoint_by_reciprocal_modified_super_cell.cast<double>();
+        // 将 q 点坐标转换为埃^-1
+        auto qpoint = (diff_of_sub_qpoint_by_reciprocal_primative_cell.transpose()
+          * (primative_cell.transpose().inverse())).transpose();
+        // 计算基矢
+        basis[i_of_sub_qpoint][i_of_basis]
+          = (2i * std::numbers::pi_v<double> * (atom_position * qpoint)).array().exp();
+      }
+    }
+    return basis;
+  }
+
+  std::vector<std::vector<double>> UnfoldSolver::construct_projection_coefficient
+  (
+    const BasisType& basis,
+    const std::vector<std::reference_wrapper<const decltype
+      (InputType::QpointDataType::ModeDataType::AtomMovement)>>& mode_data,
+    std::atomic<unsigned>& number_of_finished_modes
+  )
+  {
+    // 第一层下标对应不同模式, 第二层下标对应这个模式在反折叠后的 q 点(sub qpoint)
+    std::vector<std::vector<double>> projection_coefficient(mode_data.size());
+    // 对每个模式并行
+    std::transform
+    (
+      std::execution::par, mode_data.begin(), mode_data.end(),
+      projection_coefficient.begin(), [&](const auto& mode_data)
+      {
+        // 这里, mode_data 和 projection_coefficient 均指对应于一个模式的数据
+        std::vector<double> projection_coefficient(basis.size());
+        for (unsigned i_of_sub_qpoint = 0; i_of_sub_qpoint < basis.size(); i_of_sub_qpoint++)
+          // 对于 basis 中, 对应于单胞倒格子的部分, 以及对应于不同方向的部分, 分别求内积, 然后求模方和
+          for (unsigned i_of_basis = 0; i_of_basis < basis[i_of_sub_qpoint].size(); i_of_basis++)
+            projection_coefficient[i_of_sub_qpoint] +=
+              (basis[i_of_sub_qpoint][i_of_basis].transpose().conjugate() * mode_data.get())
+                .array().abs2().sum();
+        // 如果是严格地将向量分解到一组完备的基矢上, 那么不需要对计算得到的权重再做归一化处理
+        // 但这里并不是这样一个严格的概念. 因此对分解到各个 sub qpoint 上的权重做归一化处理
+        auto sum = std::accumulate
+          (projection_coefficient.begin(), projection_coefficient.end(), 0.);
+        for (auto& _ : projection_coefficient)
+          _ /= sum;
+        number_of_finished_modes++;
+        return projection_coefficient;
+      }
+    );
+    return projection_coefficient;
+  }
+
+  UnfoldSolver::OutputType UnfoldSolver::construct_output
+  (
+    const decltype(InputType::SuperCellMultiplier)& super_cell_multiplier,
+    const decltype(InputType::SuperCellDeformation)& super_cell_deformation,
+    const std::vector<std::reference_wrapper<const decltype
+      (InputType::QpointDataType::Qpoint)>>& meta_qpoint_by_reciprocal_super_cell,
+    const std::vector<std::vector<std::reference_wrapper<const decltype
+      (InputType::QpointDataType::ModeDataType::Frequency)>>>& frequency,
+    const ProjectionCoefficientType_& projection_coefficient
+  )
+  {
+    OutputType output;
+    for
+    (
+      unsigned i_of_meta_qpoint = 0, num_of_mode_manipulated = 0;
+      i_of_meta_qpoint < meta_qpoint_by_reciprocal_super_cell.size();
+      i_of_meta_qpoint++
+    )
+    {
+      for (auto [xyz_of_diff_of_sub_qpoint_by_reciprocal_modified_super_cell, i_of_sub_qpoint]
+        : triplet_sequence(super_cell_multiplier))
+      {
+        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 =
+        (
+          super_cell_multiplier.cast<double>().cwiseInverse().asDiagonal()
+          * (
+            xyz_of_diff_of_sub_qpoint_by_reciprocal_modified_super_cell.cast<double>()
+            + super_cell_deformation.value_or(Eigen::Matrix3d::Identity()).inverse()
+              * meta_qpoint_by_reciprocal_super_cell[i_of_meta_qpoint].get().cast<double>()
+          )
+        ).eval();
+        _.Qpoint = sub_qpoint_by_reciprocal_primative_cell.array()
+          - sub_qpoint_by_reciprocal_primative_cell.array().floor();
+        _.Source = meta_qpoint_by_reciprocal_super_cell[i_of_meta_qpoint];
+        _.SourceIndex_ = i_of_meta_qpoint;
+        for (unsigned i_of_mode = 0; i_of_mode < frequency[i_of_meta_qpoint].size(); i_of_mode++)
+        {
+          auto& __ = _.ModeData.emplace_back();
+          __.Frequency = frequency[i_of_meta_qpoint][i_of_mode];
+          __.Weight = projection_coefficient[num_of_mode_manipulated + i_of_mode][i_of_sub_qpoint];
+        }
+      }
+      num_of_mode_manipulated += frequency[i_of_meta_qpoint].size();
+    }
+    return output;
+  }
+}