暂存

CMakeLists.txt

@@ -22,13 +22,16 @@ find_package(TBB REQUIRED)
 find_package(glad REQUIRED)
 # find_package(Matplot++ REQUIRED COMPONENTS matplot_opengl)
 find_package(Matplot++ REQUIRED)
+find_package(biu REQUIRED)
+find_package(magic_enum REQUIRED)
+find_package(Boost REQUIRED COMPONENTS headers iostreams)
+find_package(range-v3 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)
+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
-  Matplot++::matplot_opengl)
+  yaml-cpp HighFive_HighFive TBB::tbb Matplot++::matplot Matplot++::matplot_opengl biu::biu)
 # target_compile_definitions(ufo PRIVATE $<$<CONFIG:Debug>:_GLIBCXX_DEBUG>)
 
 get_property(ImportedTargets DIRECTORY "${CMAKE_SOURCE_DIR}" PROPERTY IMPORTED_TARGETS)

flake.lock

@@ -807,11 +807,11 @@
         "touchix": "touchix"
       },
       "locked": {
-        "lastModified": 1696743278,
-        "narHash": "sha256-ckBxRe2wjUiUzDQWs/WM+x0HgnRUjoZiHzO2c408dFg=",
+        "lastModified": 1697296772,
+        "narHash": "sha256-xHsRobR3rQqVrBKmLnK3QnUAdrhGQ92vbqJlUKyvl1o=",
         "owner": "CHN-beta",
         "repo": "nixos",
-        "rev": "5ffdec57c08307cc2589ec57169d72cb68cf7c8b",
+        "rev": "a56011bf6df896887256c1da7e07accab81f97d2",
         "type": "github"
       },
       "original": {
@@ -1015,11 +1015,11 @@
     },
     "nixpkgs_5": {
       "locked": {
-        "lastModified": 1696604326,
-        "narHash": "sha256-YXUNI0kLEcI5g8lqGMb0nh67fY9f2YoJsILafh6zlMo=",
+        "lastModified": 1697059129,
+        "narHash": "sha256-9NJcFF9CEYPvHJ5ckE8kvINvI84SZZ87PvqMbH6pro0=",
         "owner": "NixOS",
         "repo": "nixpkgs",
-        "rev": "87828a0e03d1418e848d3dd3f3014a632e4a4f64",
+        "rev": "5e4c2ada4fcd54b99d56d7bd62f384511a7e2593",
         "type": "github"
       },
       "original": {

flake.nix

@@ -17,7 +17,7 @@
       {
         packages = with pkgs; [ pkg-config cmake ninja ];
         buildInputs = (with pkgs; [ eigen yaml-cpp fmt highfive tbb_2021_8.dev glfw libGL range-v3 ])
-          ++ (with localPackages; [ concurrencpp matplotplusplus zpp-bits glad ]);
+          ++ (with localPackages; [ concurrencpp matplotplusplus zpp-bits glad biu ]);
         hardeningDisable = [ "all" ];
         # NIX_DEBUG = "1";
       };

include/ufo/solver.hpp

@@ -25,6 +25,7 @@
 # include <matplot/matplot.h>
 # include <matplot/backend/opengl.h>
 # include <range/v3/all.hpp>
+# include <biu.hpp>
 
 // 在相位中, 约定为使用 $\exp (2 \pi i \vec{q} \cdot \vec{r})$ 来表示原子的运动状态
 //  (而不是 $\exp (-2 \pi i \vec{q} \cdot \vec{r})$)
@@ -101,11 +102,6 @@ namespace ufo
 
   template <std::size_t N> struct ToSpanHelperType {};
   template <std::size_t N> decltype(auto) operator|(auto&& input, const ToSpanHelperType<N>&);
-
-  struct ToEigenHelper {};
-  template <typename T, std::size_t N> decltype(auto) operator|(const std::span<T, N>& input, const ToEigenHelper&);
-  template <typename T, std::size_t M, std::size_t N> decltype(auto) operator|
-    (const std::span<std::span<T, N>, M>& input, const ToEigenHelper&);
 }
 
 HIGHFIVE_REGISTER_TYPE(ufo::PhonopyComplex, ufo::create_compound_complex)

include/ufo/solver.tpp

@@ -130,26 +130,4 @@ namespace ufo
     else
       return std::span(input.data(), N);
   }
-
-  inline constexpr ToEigenHelper to_eigen;
-  template <typename T, std::size_t N> inline decltype(auto) operator|
-    (const std::span<T, N>& input, const ToEigenHelper&)
-  {
-    if constexpr (N == std::dynamic_extent)
-      return Eigen::VectorX<T>(input.data(), input.size());
-    else
-      return Eigen::Vector<T, N>(input.data());
-  }
-  template <typename T, std::size_t M, std::size_t N> inline decltype(auto) operator|
-    (const std::span<std::span<T, N>, M>& input, const ToEigenHelper&)
-  {
-    if constexpr (M == std::dynamic_extent)
-      if constexpr (N == std::dynamic_extent)
-      {
-        Eigen::Matrix
-        return Eigen::MatrixX<T>(input.data()->data(), input.size(), input.data()->size());
-      }
-      else
-        return Eigen::Matrix<T, Eigen::Dynamic, N>(input.data()->data(), input.size());
-  }
 }

include/ufo/unfold.hpp

@@ -54,7 +54,7 @@ namespace ufo
             // 这个数据是这样得到的: phonopy 输出的动态矩阵的 eigenvector 乘以 $\exp(-2 \pi i \vec q \cdot \vec r)$
             // 这个数据可以认为是原子位移中, 关于超胞有周期性的那一部分, 再乘以原子质量的开方.
             // 这个数据在读入后会被立即归一化.
-            Eigen::MatrixX3cd AtomMovement;
+            Eigen::MatrixX3cd EigenVector;
           };
           std::vector<ModeDataType> ModeData;
         };

include/ufo/unfold.impl.hpp

@@ -9,29 +9,12 @@ namespace ufo
     // read main input file
     {
       auto node = YAML::LoadFile(filename);
-      PrimativeCell = [](YAML::Node node)
-      {
-        Eigen::Matrix3d matrix;
-        for (unsigned i = 0; i < 3; i++)
-          for (unsigned j = 0; j < 3; j++)
-            matrix(i, j) = node[i][j].as<double>();
-        return matrix;
-      }(node["PrimativeCell"]);
-
-      for (unsigned i = 0; i < 3; i++)
-        SuperCellMultiplier(i) = node["SuperCellMultiplier"][i].as<int>();
-
+      PrimativeCell = node["PrimativeCell"].as<std::vector<std::vector<double>>>() | biu::toEigen<3, 3>;
+      SuperCellMultiplier = node["SuperCellMultiplier"].as<std::vector<unsigned>>() | toEigenVector<3>;
       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>();
-
+        SuperCellDeformation = value.as<std::vector<std::vector<double>>>() | toEigenMatrix<3, 3>;
+      PrimativeCellBasisNumber = node["PrimativeCellBasisNumber"].as<std::vector<unsigned>>()
+        | toEigenVector<3>;
       AtomPositionInputFile = DataFile
       (
         node["AtomPositionInputFile"], {"yaml"},
@@ -43,29 +26,30 @@ namespace ufo
         filename, true
       );
       if (auto value = node["QpointDataOutputFile"])
-      {
-        QpointDataOutputFile.resize(value.size());
-        for (unsigned i = 0; i < value.size(); i++)
-          QpointDataOutputFile[i] = DataFile
+        for (auto&& v : value)
+          QpointDataOutputFile.push_back(DataFile
           (
-            value[i], {"yaml", "yaml-human-readable", "zpp", "hdf5"},
+            v, {"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>();
+      YAML::Node source = node["points"];
+      if (!source)
+        source = node["unit_cell"]["points"];
+      auto atom_position_to_super_cell = source
+        | std::views::transform([](YAML::Node&& v)
+          { return v["coordinates"].as<std::vector<double>>(); })
+        | std::
+        | toEigenMatrix<std::dynamic_extent, 3>
+
+      auto atom_position_to_super_cell = Eigen::MatrixX3d(source, 3);
+      for (unsigned i = 0; i < source.size(); i++)
+        atom_position_to_super_cell.row(i) = source[i]["coordinates"].as<std::vector<double>>()
+          | toEigenVector<3>;
       auto super_cell = (SuperCellDeformation.value_or(Eigen::Matrix3d::Identity())
         * SuperCellMultiplier.cast<double>().asDiagonal() * PrimativeCell).eval();
       AtomPosition = atom_position_to_super_cell * super_cell;
@@ -73,7 +57,83 @@ namespace ufo
     if (QpointDataInputFile.Format == "yaml")
     {
       auto node = YAML::LoadFile(QpointDataInputFile.Filename);
+      QpointData = node["phonon"].as<std::vector<YAML::Node>>()
+        | std::views::transform([](YAML::Node&& v)
+          {
+            return QpointDataType
+            {
+              .Qpoint = v["q-position"].as<std::vector<double>>() | toEigenVector<3>,
+              .ModeData = std::vector(
+                v["band"].as<std::vector<YAML::Node>>()
+                | std::views::transform([](YAML::Node&& v)
+                  {
+                    return QpointDataType::ModeDataType
+                    {
+                      .Frequency = v["frequency"].as<double>(),
+                      .EigenVector = 
+                      v["frequency"].as<double>(),
+                      Eigen::MatrixX3cd(v["eigenvector"]
+                        .as<std::vector<std::vector<std::vector<double>>>>()
+                        | std::views::transform([](auto&& v)
+                          {
+                            return Eigen::Vector3cd(v | toEigenVector<2>);
+                          })
+                        | toEigenMatrix)
+                    };
+                  })
+                | toStdVector
+              )
+
+              )
+
+              v["q-position"].as<std::vector<double>>() | toEigenVector<3>,
+              v["band"].as<std::vector<YAML::Node>>()
+              | std::views::transform([](auto&& v)
+                {
+                  return ModeDataType
+                  {
+                    v["frequency"].as<double>(),
+                    Eigen::MatrixX3cd(v["eigenvector"]
+                      .as<std::vector<std::vector<std::vector<double>>>>()
+                      | std::views::transform([](auto&& v)
+                        {
+                          return Eigen::Vector3cd(v | toEigenVector<2>);
+                        })
+                      | toEigenMatrix)
+                  };
+                })
+              | toStdVector
+            };
+
+            QpointDataType output;
+            for (unsigned i = 0; i < 3; i++)
+              output.Qpoint(i) = v["q-position"][i].as<double>();
+            output.ModeData = v["band"].as<std::vector<YAML::Node>>()
+              | std::views::transform([](auto&& v)
+                {
+                  ModeDataType output;
+                  output.Frequency = v["frequency"].as<double>();
+                  auto eigenvector_vectors = v["eigenvector"]
+                    .as<std::vector<std::vector<std::vector<double>>>>();
+                  Eigen::MatrixX3cd eigenvectors(AtomPosition.rows(), 3);
+                  for (unsigned i = 0; i < AtomPosition.rows(); i++)
+                    for (unsigned j = 0; j < 3; j++)
+                      eigenvectors(i, j)
+                        = eigenvector_vectors[i][j][0] + 1i * eigenvector_vectors[i][j][1];
+                  output.EigenVector = eigenvectors / eigenvectors.norm();
+                  return output;
+                })
+              | toStdVector;
+            return output;
+          })
+
+
       auto phonon = node["phonon"].as<std::vector<YAML::Node>>();
+
+
+
+
+      QpointData = std::vector<QpointDataType>(phonon.size());
       QpointData.resize(phonon.size());
       for (unsigned i = 0; i < phonon.size(); i++)
       {