merge fold into main

CMakeLists.txt

@@ -22,7 +22,7 @@ find_package(TBB REQUIRED)
 find_package(Matplot++ REQUIRED)
 find_path(ZPP_BITS_INCLUDE_DIR zpp_bits.h REQUIRED)
 
-add_executable(ufo src/unfold.cpp src/plot.cpp src/main.cpp)
+add_executable(ufo 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)

fold.cpp

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

include/ufo/fold.hpp

@@ -0,0 +1,33 @@
+# pragma once
+# include <ufo/solver.hpp>
+
+namespace ufo
+{
+  class FoldSolver : public Solver
+  {
+    public:
+      struct InputType
+      {
+        Eigen::Matrix3d PrimativeCell;
+        Eigen::Vector<unsigned, 3> SuperCellMultiplier;
+        std::optional<Eigen::Matrix<double, 3, 3>> SuperCellDeformation;
+        std::vector<Eigen::Vector3d> Qpoints;
+        std::string OutputFilename;
+
+        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;
+      static Eigen::Vector3d fold(Eigen::Vector3d qpoint, Eigen::Vector<unsigned, 3> super_cell_multiplier,
+        std::optional<Eigen::Matrix<double, 3, 3>> super_cell_deformation);
+  };
+}

src/fold.cpp

@@ -0,0 +1,87 @@
+# 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++)
+      for (unsigned j = 0; j < 3; j++)
+        PrimativeCell(i, j) = input["PrimativeCell"][i][j].as<double>();
+    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)}});
+    OutputFilename = input["OutputFilename"].as<std::string>();
+  }
+  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_.OutputFilename);
+    return *this;
+  }
+
+  Eigen::Vector3d FoldSolver::fold
+  (
+    Eigen::Vector3d qpoint, Eigen::Vector<unsigned, 3> super_cell_multiplier,
+    std::optional<Eigen::Matrix<double, 3, 3>> super_cell_deformation
+  )
+  {
+    // 首先需要将 q 点转移到超胞的倒格子中
+    // 将 q 点坐标扩大, 然后取小数部分, 就可以了
+    auto qpoint_by_reciprocal_modified_super_cell = super_cell_multiplier.cast<double>().asDiagonal() * qpoint;
+    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;
+  }
+}

src/main.cpp

@@ -1,10 +1,14 @@
+# 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("unfold"))
+  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]}();

test.cpp

@@ -1,27 +0,0 @@
-# include <eigen3/Eigen/Dense>
-# include <zpp_bits.h>
-# include <span>
-# include <iostream>
-# include <vector>
-
-namespace Eigen
-{
-  constexpr auto serialize(auto & archive, Eigen::Matrix3d& matrix)
-  {
-    return archive(std::span<double, 9>(matrix.data(), 9));
-  }
-  constexpr auto serialize(auto & archive, const Eigen::Matrix3d& matrix)
-  {
-    return archive(std::span<const double, 9>(matrix.data(), 9));
-  }
-}
-
-int main()
-{
-  auto [data, in, out] = zpp::bits::data_in_out();
-  std::vector<Eigen::Matrix3d> PrimativeCell = { Eigen::Matrix3d::Identity() };
-  out(PrimativeCell).or_throw();
-  std::vector<Eigen::Matrix3d> PrimativeCell2;
-  in(PrimativeCell2).or_throw();
-  std::cout << PrimativeCell2[0] << std::endl;
-}