PrimativeCellBasisNumber 使用 gamma center 的取法

AtomTranslation 改为以超胞晶格矢为单位，而不是以埃为单位
一些小修改，应该不影响结果
2026-01-11 14:53:12 +08:00 · 2026-01-08 12:36:14 +08:00 · 2026-01-08 12:36:14 +08:00 · 2026-01-05 15:21:35 +08:00 · 2025-11-05 12:05:12 +08:00 · 2025-11-05 12:02:52 +08:00
122 changed files with 1732 additions and 3180 deletions
--- a/.envrc
+++ b/.envrc
@@ -1 +1 @@
-use flake
+use flake nixos#ufo
--- a/.gitattributes
+++ b/.gitattributes
@@ -1,4 +1,2 @@
 *.hdf5 filter=lfs diff=lfs merge=lfs -text
 *.zpp filter=lfs diff=lfs merge=lfs -text
-*POSCAR* filter=lfs diff=lfs merge=lfs -text
-*.png filter=lfs diff=lfs merge=lfs -text
--- a/.gitignore
+++ b/.gitignore
@@ -1,18 +1,4 @@
-.vscode
 .direnv
-main
-plot
-out.yaml
-.ccls-cache
+build
+.vscode
 .cache
-build/
-test/*
-!test/14.2.2.4.unfold.yaml
-!test/14.2.2.5.unfold.yaml
-!test/14.2.2.5.plot.yaml
-!test/14.2.4.4.unfold.yaml
-!test/14.2.4.4.plot.yaml
-!test/14.2.5.5.unfold.yaml
-!test/14.2.5.5.plot.yaml
-!test/14.2.6.4.unfold.yaml
-!test/14.2.6.4.plot.yaml
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -1,36 +1,51 @@
 cmake_minimum_required(VERSION 3.14)
-project(ufo VERSION 0.0.0 LANGUAGES CXX)
+project(ufo VERSION 0 LANGUAGES CXX)
 enable_testing()
 include(GNUInstallDirs)

-set_property(GLOBAL PROPERTY CXX_STANDARD 23)
-set_property(GLOBAL PROPERTY CXX_STANDARD_REQUIRED ON)
-set_property(GLOBAL PROPERTY CXX_EXTENSIONS OFF)
-
 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 COMPONENTS matplot_opengl)
 find_package(Matplot++ REQUIRED)
-find_path(ZPP_BITS_INCLUDE_DIR zpp_bits.h REQUIRED)
+find_package(biu REQUIRED)
+find_package(Threads REQUIRED)
+message("biu dir: ${biu_DIR}")

-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_definitions(ufo PRIVATE $<$<CONFIG:Debug>:_GLIBCXX_DEBUG>)
+add_executable(ufo src/fold.cpp src/unfold.cpp src/project-to-mode.cpp
+  src/plot.cpp src/raman.cpp src/project-to-atom.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_definitions(ufo PRIVATE BIU_LOGGER_SOURCE_ROOT="${CMAKE_SOURCE_DIR}")
+if (CMAKE_BUILD_TYPE STREQUAL "Debug")
+  target_compile_definitions(ufo PRIVATE BIU_LOGGER_DEBUG)
+  target_compile_options(ufo PRIVATE -ftemplate-backtrace-limit=0)
+endif()
+
+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}")
+message("CMake build type: ${CMAKE_BUILD_TYPE}")
+
+include(CTest)
+
+add_test(NAME fold WORKING_DIRECTORY ${PROJECT_SOURCE_DIR}/test/fold COMMAND ufo fold config.yaml)
+add_test(NAME unfold WORKING_DIRECTORY ${PROJECT_SOURCE_DIR}/test/unfold COMMAND ufo unfold config.yaml)
+add_test(NAME zpp-to-yaml WORKING_DIRECTORY ${PROJECT_SOURCE_DIR}/test/zpp-to-yaml
+  COMMAND sh -c "${CMAKE_BINARY_DIR}/ufo zpp-to-yaml < data.zpp > data.yaml")
+add_test(NAME project-to-mode WORKING_DIRECTORY ${PROJECT_SOURCE_DIR}/test/project-to-mode
+  COMMAND ufo project-to-mode config.yaml)
+add_test(NAME project-to-atom WORKING_DIRECTORY ${PROJECT_SOURCE_DIR}/test/project-to-atom
+  COMMAND ufo project-to-atom config.yaml)
+add_test(NAME raman-create-displacement WORKING_DIRECTORY ${PROJECT_SOURCE_DIR}/test/raman-create-displacement
+  COMMAND ufo raman-create-displacement config.yaml)
+add_test(NAME raman-extract WORKING_DIRECTORY ${PROJECT_SOURCE_DIR}/test/raman-extract
+  COMMAND ufo raman-extract .)
+add_test(NAME raman-apply-contribution WORKING_DIRECTORY ${PROJECT_SOURCE_DIR}/test/raman-apply-contribution
+  COMMAND ufo raman-apply-contribution config.yaml)
--- a/README.md
+++ b/README.md
@@ -1,17 +0,0 @@
-This is a program to approximately unfold phonon spectra.
-
-Currently, there is no documentation, and some not core features are in other branches
-  and implemented in a temporary way, e.g., directly writing data in source file.
-
-I will finish the documentation and merge these features into main branch ASAP.
-
-The program is packaged in a standard way of cmake, so that you can build it just like other cmake projects:
-
-```bash
-mkdir build
-cd build
-cmake ..
-make
-```
-
-and view source files for input and output format.
--- a/default.nix
+++ b/default.nix
@@ -0,0 +1,11 @@
+{
+  stdenv, cmake, pkg-config, version ? null,
+  tbb, matplotplusplus, biu
+}: stdenv.mkDerivation
+{
+  name = "ufo";
+  src = ./.;
+  buildInputs = [ tbb matplotplusplus biu ];
+  nativeBuildInputs = [ cmake pkg-config ];
+  doCheck = true;
+}
--- a/doc/README.md
+++ b/doc/README.md
@@ -0,0 +1,45 @@
+分为几个功能：
+
+* fold：根据要计算的单胞的 q 点路径，计算超胞中对应的 q 点路径，生成的路径再交给 phonopy 计算。
+* unfold：根据 phonopy 计算的结果，将超胞的结果展开到单胞中。
+* plot：对计算结果画图。
+
+主要的输入输出格式均为 yaml。对于数据特别大的情况，也可以从 hdf5 中读取一部分数据或者将一部分数据写入到 hdf5 文件中。
+
+# fold
+
+## 输入
+
+```yaml
+# 三个整数组成的向量，表示从单胞到超胞，三个晶格矢量的倍数
+# 必写
+SuperCellMultiplier: [2, 2, 2]
+# 一个变换矩阵，表明超胞经历了怎样的扭曲。
+# 可选，默认值为单位矩阵
+SuperCellDeformation: [[1, 0, 0], [0, 1, 0], [0, 0, 1]]
+# 一个由三个浮点数组成的向量，表示考虑的 q 点
+# 必写
+Qpoints:
+  - [0, 0, 0]
+  - [0.1, 0, 0]
+  - [0.2, 0, 0]
+  - [0.3, 0, 0]
+  - [0.4, 0, 0]
+  - [0.5, 0, 0]
+# 一个 DataFile 类型的对象，表明输出结果到哪个文件
+# 必写
+OutputFile:
+```
+
+## 输出
+
+```yaml
+# 得到的 q 点坐标
+Qpoints:
+  - [0, 0, 0]
+  - [0.1, 0, 0]
+  - [0.2, 0, 0]
+  - [0.3, 0, 0]
+  - [0.4, 0, 0]
+  - [0.5, 0, 0]
+```
--- a/example/README.md
+++ b/example/README.md
@@ -1,139 +0,0 @@
-本例子演示一个将 4H-SiC 超胞中声子反折叠的例子。
-超胞中包含 95 个原子，带有一个 C 空位；我们将把它近似反折叠到仅包含 8 个原子的原胞上。
-每一步用到的文件都将给出。简单起见，我们将仅仅计算 Gamma-M 线上的声子谱。
-
-# 第一步：建立模型
-
-首先需要准备好已经已经充分弛豫的原胞模型，和在这个原胞模型基础上构造的、未弛豫的、带一个 C 空位的超胞模型。
-
-原胞和超胞的格矢矩阵分别为：
-```yaml
-PrimitiveCell:
-  - [ 3.0799, 0, 0 ]
-  - [ -1.53995, 2.66727, 0 ]
-  - [ 0, 0, 10.0822 ]
-SuperCell:
-  - [ 9.2397, 0, 0 ]
-  - [ -4.61985, 8.00181, 0 ]
-  - [ 0, 0, 10.0822 ]
-```
-
-这里的每一行代表一个格矢（与 POSCAR 相同），单位为埃。
-
-超胞的格矢矩阵可以看作是原胞的格矢矩阵左乘一个变换矩阵得到的：
-```yaml
-Transformation:
-  - [ 3, 0, 0 ]
-  - [ 2, 4, 0 ]
-  - [ 0, 0, 1 ]
-```
-
-变换矩阵的一定都是整数（否则将无法保持周期性）。
-
-得到超胞模型后，对它进行弛豫。稍后的声子反折叠步骤中，需要用到弛豫后超胞模型中各个原子的位置。
-
-注意：在稍后的声子反折叠步骤中，使用的超胞格矢矩阵是弛豫前的数据，而不是弛豫后的数据。
-弛豫过程中超胞的格矢矩阵的变化通常很小，在反折叠算法中不予考虑。
-
-注意，在从原胞模型构造超胞模型时，不能整体平移或者旋转坐标系，也不能整体平移或者旋转原子，
-  否则最终计算出来的结果将不正确。
-缺陷导致的少数原子的位移是可以接受。
-
-# 第二步：计算超胞中的 Q 点坐标
-
-思考一个问题：原胞中 Gamma-M 线上的声子模式，对应于超胞中哪些 Q 点的模式？
-这个问题经过仔细思考后并不难回答，尤其是对于那些变换矩阵为对角矩阵的超胞。
-但对于一般的超胞，手动计算时容易出错，因此我写了一个功能来处理这个问题。
-
-假定我们要求从 Gamma 到 M 点的路径上，共 11 个 Q 点的声子模式。
-这 11 个 Q 点的坐标为 $(0, 0, 0)$ 到 $(0.5, 0, 0)$。
-我们可以使用以下配置文件来计算这些 Q 点在超胞中的坐标：
-
-```yaml
-Transformation:
-  - [ 3, 0, 0 ]
-  - [ 2, 4, 0 ]
-  - [ 0, 0, 1 ]
-Qpoints:
-  - [0, 0, 0]
-  - [0.05, 0, 0]
-  - [0.1, 0, 0]
-  - [0.15, 0, 0]
-  - [0.2, 0, 0]
-  - [0.25, 0, 0]
-  - [0.3, 0, 0]
-  - [0.35, 0, 0]
-  - [0.4, 0, 0]
-  - [0.45, 0, 0]
-  - [0.5, 0, 0]
-OutputFile:
-  Filename: fold-output.yaml
-  Format: yaml
-```
-
-运行命令：
-```bash
-ufo fold fold.yaml
-```
-
-得到：
-
-```yaml
-Qpoints:
-  - [ 0.00000000, 0.00000000, 0.00000000 ]
-  - [ 0.15000000, 0.09999999, 0.00000000 ]
-  - [ 0.30000000, 0.19999998, 0.00000000 ]
-  - [ 0.45000000, 0.29999997, 0.00000000 ]
-  - [ 0.60000000, 0.39999996, 0.00000000 ]
-  - [ 0.75000000, 0.49999995, 0.00000000 ]
-  - [ 0.90000000, 0.59999994, 0.00000000 ]
-  - [ 0.05000000, 0.69999993, 0.00000000 ]
-  - [ 0.20000000, 0.79999992, 0.00000000 ]
-  - [ 0.35000000, 0.89999991, 0.00000000 ]
-  - [ 0.50000000, 0.99999990, 0.00000000 ]
-```
-
-结果表明，原胞中 Gamma-M 线经过格矢变换后分裂成了两条，
-  一条是 $(0, 0, 0)$ 点到 $(1, \frac23, 0)$ 点，另一条是 $(0, \frac23, 0)$ 点到 $(0.5, 1, 0)$ 点。
-
-大致在草稿纸上画个图，可以确认这个结果是合理的：第二段在下一个周期中与第一段相连。
-
-有时可以据此略微调整原胞中所求点的位置，使得在超胞中此处的声子可以严格求出而不需插值。
-此例子不适合演示这个方法。
-
-# 第三步：计算声子谱
-
-使用 Phonopy 计算超胞中，上一步中得到的几条路径上的声子。
-
-你可以使用任意你认为合适的参数来计算，只要能得到声子谱，之后的计算都可以进行下去。
-
-注意：这一步需要的计算量往往比较大，常常需要计算几百或几千个模型中的原子受力。
-  建议先计算其中一个并输出 WAVECAR，然后将 WAVECAR 软连接到其他计算中，以加速计算。
-
-注意：如果 phonopy 输出的数据太大，可以使用 hdf5 格式的输出以加速处理。
-
-# 第四步：声子谱反折叠
-
-反折叠需要的配置比较多且复杂。
-
-```yaml
-PrimativeCell:
-  - [ 3.0799, 0, 0 ]
-  - [ -1.53995, 2.66727, 0 ]
-  - [ 0, 0, 10.0822 ]
-Transformation:
-  - [ 3, 0, 0 ]
-  - [ 2, 4, 0 ]
-  - [ 0, 0, 1 ]
-PrimativeCellBasisNumber: [ 8, 8, 8 ]
-AtomPositionInputFile:
-  Format: yaml
-  Filename: phonopy.yaml # 或 phonopy_disp.yaml
-QpointDataInputFile:
-  Format: yaml
-  Filename: phonopy.yaml
-QpointDataOutputFile:
-  - Format: yaml
-    Filename: phonopy_folded.yaml
-```
-
--- a/example/step1/primitive_cell/POSCAR
+++ b/example/step1/primitive_cell/POSCAR
--- a/example/step1/super_cell/CONTCAR
+++ b/example/step1/super_cell/CONTCAR
@@ -1,103 +0,0 @@
-POSCAR file written by Ovito 3.7.11     
-    1.000000000000000     
-     9.2673899685017957    0.0000000000000000    0.0000000000000000
-     0.0000000000000000   10.6774221488510186   -0.0181004282034319
-     0.0000000000000000   -0.0171993991549284   10.1116441082647377
-  Si             C            
-              48              47
-Direct
- -0.0000557863121398  0.1668131534569660  0.2495178577809962
-  0.1666666570000004  0.0821651768612523  0.7494996261803192
-  0.0000785931726128 -0.0008014364573324 -0.0019171777870414
- -0.0001332313405060  0.0010323603084980  0.4990245675767412
-  0.8333163913323451  0.4170889259273483  0.2509448387496841
-  0.0010593515672643  0.3339499564673934  0.7504013581078407
-  0.8333234221246683  0.2494366873144543 -0.0001913591111189
-  0.8335033186823803  0.2508241081996696  0.5000431170293115
-  0.6666666659999976  0.6667639132151945  0.2493250038919751
-  0.8328430132025012  0.5824561773357575  0.7502934292698500
-  0.6666666670000012  0.4993068721873999  0.0001383770014351
-  0.6666666670000012  0.5005472005797680  0.5002523804663008
-  0.4999087106685396  0.9169692013808204  0.2492033506084066
-  0.6666666570000004  0.8337032368714793  0.7504363269561648
-  0.4995637733746955  0.7493807863038214 -0.0008564098823073
-  0.5003725611872856  0.7503863239197331  0.4992410940474822
-  0.3333891193121384  0.1668131534569660  0.2495178577809962
-  0.4992733589909116  0.0818752444180838  0.7497937354374519
-  0.3332547398273861 -0.0008014364573324 -0.0019171777870414
-  0.3334665643405045  0.0010323603084980  0.4990245675767412
-  0.1666666659999976  0.4169480526102884  0.2506028916346738
-  0.3322739714327337  0.3339499564673934  0.7504013581078407
-  0.1666666670000012  0.2497373319249146 -0.0000499395300293
-  0.1666666670000012  0.2505366718958373  0.4997757704091679
-  0.0007343296568028  0.6675813991651439  0.2507455593024008
-  0.1666666570000004  0.5882597899202560  0.7580329988559029
- -0.0006769397216459  0.4989432733446046  0.0018414165511462
-  0.0010383420437462  0.5016621987383314  0.5008853395289279
-  0.8334246223314590  0.9169692013808204  0.2492033506084066
-  0.0000682038488344  0.8265801375506314  0.7450704126627182
-  0.8337695596253030  0.7493807863038214 -0.0008564098823073
-  0.8329607718127133  0.7503863239197331  0.4992410940474822
-  0.6666666659999976  0.1670868510386173  0.2499476084075379
-  0.8340599550090892  0.0818752444180838  0.7497937354374519
-  0.6666666670000012 -0.0005225484838031 -0.0002518136909808
-  0.6666666670000012  0.0005098915308569  0.4992736940964735
-  0.5000169416676533  0.4170889259273483  0.2509448387496841
-  0.6666666570000004  0.3320791225644938  0.7499985648946081
-  0.5000099108753309  0.2494366873144543 -0.0001913591111189
-  0.4998300143176183  0.2508241081996696  0.5000431170293115
-  0.3325990033431961  0.6675813991651439  0.2507455593024008
-  0.5004902997974959  0.5824561773357575  0.7502934292698500
-  0.3340102727216447  0.4989432733446046  0.0018414165511462
-  0.3322949909562527  0.5016621987383314  0.5008853395289279
-  0.1666666659999976  0.9180544950396436  0.2478666411127934
-  0.3332651191511635  0.8265801375506314  0.7450704126627182
-  0.1666666670000012  0.7492229363186385  0.0001611773946995
-  0.1666666670000012  0.7586221208203656  0.5025559614474544
- -0.0003543855561368  0.1676511882651485  0.4376401280043795
-  0.1666666570000004  0.0828996116757371  0.9371848303263415
- -0.0003591309494861  0.0000023801845120  0.1856212094249273
-  0.0015663355615017 -0.0026327726480436  0.6863958225036811
-  0.8337845559136262  0.4178115825074205  0.4389674348305436
- -0.0001023038646378  0.3327868228709865  0.9381885158054116
-  0.8333373011703021  0.2502176615540846  0.1873187173259448
-  0.8342772552129253  0.2495543040285077  0.6873480333936619
-  0.6666666659999976  0.6669417550050203  0.4369673987476388
-  0.8333807305605515  0.5824145944486363  0.9381784129487265
-  0.6666666670000012  0.4997745444284661  0.1874992047550911
-  0.6666666670000012  0.4986742269313283  0.6876164226427359
-  0.4999339071178296  0.9171919757061141  0.4373995923163460
-  0.6666666570000004  0.8325907836040870  0.9377955619721441
-  0.4994945329888917  0.7500616410188495  0.1869506467429709
-  0.5026656836615087  0.7480725752499672  0.6860625197385234
-  0.3336877185561356  0.1676511882651485  0.4376401280043795
-  0.5000447115614929  0.0824669311545484  0.9373067134583003
-  0.3336924639494850  0.0000023801845120  0.1856212094249273
-  0.3317669974384970 -0.0026327726480436  0.6863958225036811
-  0.1666666659999976  0.4174122540170955  0.4385569558692831
-  0.3334356268646358  0.3327868228709865  0.9381885158054116
-  0.1666666670000012  0.2500772537096505  0.1873263641549840
-  0.1666666670000012  0.2496312870965579  0.6872470991449243
-  0.0013188200046441  0.6685626427893232  0.4398333944722485
-  0.1666666570000004  0.5809667172472394  0.9434672115062485
- -0.0002064199071260  0.5002787858188761  0.1889409420456360
-  0.0015430640609559  0.5024439923185628  0.6881853863336864
-  0.8333994258821685  0.9171919757061141  0.4373995923163460
-  0.0009625351257081  0.8317674754963278  0.9354088923528614
-  0.8338388000111076  0.7500616410188495  0.1869506467429709
-  0.8306676493384901  0.7480725752499672  0.6860625197385234
-  0.6666666659999976  0.1675291153356092  0.4379848675305987
-  0.8332886024385080  0.0824669311545484  0.9373067134583003
-  0.6666666670000012  0.0001944358287058  0.1870538575713719
-  0.6666666670000012 -0.0002360120450371  0.6867185860979695
-  0.4995487770863726  0.4178115825074205  0.4389674348305436
-  0.6666666570000004  0.3324461382685149  0.9378848515350188
-  0.4999960318296968  0.2502176615540846  0.1873187173259448
-  0.4990560777870729  0.2495543040285077  0.6873480333936619
-  0.3320145129953548  0.6685626427893232  0.4398333944722485
-  0.4999525824394457  0.5824145944486363  0.9381784129487265
-  0.3335397529071246  0.5002787858188761  0.1889409420456360
-  0.3317902689390430  0.5024439923185628  0.6881853863336864
-  0.1666666659999976  0.9216121818138033  0.4343691564472280
-  0.3323707878742898  0.8317674754963278  0.9354088923528614
-  0.1666666670000012  0.7508040359125291  0.1873373990293362
--- a/example/step1/super_cell/INCAR
+++ b/example/step1/super_cell/INCAR
@@ -1,33 +0,0 @@
-# write many output files
-LWAVE = T
-LCHARG = F
-LELF = F
-LEPSILON = F
-LCALCEPS = F
-
-# use high precision
-PREC = Accurate
-EDIFF = 1E-8
-EDIFFG = -1E-4
-
-# enough ionic steps
-NSW = 100
-
-# calculate forces and stress tensor, and allow all degrees of freedom (atom positions, cell shape, and cell volume)
-ISIF = 3
-
-# use conjugate gradient algorithm to do ionic relaxation
-IBRION = 2
-
-# use Gaussian smearing to set partial occupancies for each orbital, commonly suggested
-ISMEAR = 0
-
-# use Vosko-Wilk-Nusair interpolation to enchange the estimation of magnetic moments and energies
-VOSKOWN = 1
-
-# set symmetry
-ISYM = 2
-
-# projection done in reciprocal space for small unit cells,
-# in real space for large unit cells
-LREAL = F
--- a/example/step1/super_cell/KPOINTS
+++ b/example/step1/super_cell/KPOINTS
@@ -1,5 +0,0 @@
-3X3X3   
- 0
-Gamma    
- 3  3  3
- 0  0  0
--- a/example/step1/super_cell/POSCAR
+++ b/example/step1/super_cell/POSCAR
--- a/example/step2/fold-output.yaml
+++ b/example/step2/fold-output.yaml
@@ -1,12 +0,0 @@
-Qpoints:
-  - [ 0.00000000, 0.00000000, 0.00000000 ]
-  - [ 0.15000000, 0.09999999, 0.00000000 ]
-  - [ 0.30000000, 0.19999998, 0.00000000 ]
-  - [ 0.45000000, 0.29999997, 0.00000000 ]
-  - [ 0.60000000, 0.39999996, 0.00000000 ]
-  - [ 0.75000000, 0.49999995, 0.00000000 ]
-  - [ 0.90000000, 0.59999994, 0.00000000 ]
-  - [ 0.05000000, 0.69999993, 0.00000000 ]
-  - [ 0.20000000, 0.79999992, 0.00000000 ]
-  - [ 0.35000000, 0.89999991, 0.00000000 ]
-  - [ 0.50000000, 0.99999990, 0.00000000 ]
--- a/example/step3/1_build_models/build.zsh
+++ b/example/step3/1_build_models/build.zsh
@@ -1,4 +0,0 @@
-#!/usr/bin/env zsh
-
-cp ../../step1/super_cell/CONTCAR POSCAR
-phonopy ./phonopy.conf
--- a/example/step3/1_build_models/phonopy.conf
+++ b/example/step3/1_build_models/phonopy.conf
@@ -1,3 +0,0 @@
-CREATE_DISPLACEMENTS = .TRUE.
-DIM = 1 1 1
-SYMMETRY_TOLERANCE = 1e-4
--- a/example/step3/2_run_selfconsistence/build.zsh
+++ b/example/step3/2_run_selfconsistence/build.zsh
@@ -1,10 +0,0 @@
-#!/usr/bin/env zsh
-
-for i in {001..316}; do
-  mkdir -p job/$i
-  pushd job/$i
-    cp ../../../1_build_models/POSCAR-$i POSCAR
-    cp ../../{INCAR,POTCAR,KPOINTS} .
-    ln -s ../../WAVECAR .
-  popd
-done
--- a/example/step4/band.hdf5
+++ b/example/step4/band.hdf5
--- a/example/step4/config.yaml
+++ b/example/step4/config.yaml
@@ -1,108 +0,0 @@
-PrimativeCell:
-  - [  3.0799000263, 0.0000000000,  0.0000000000 ]
-  - [ -1.5399500132, 2.6672716639,  0.0000000000 ]
-  - [  0.0000000000, 0.0000000000, 10.0822000504 ]
-SuperCellTransformation:
-  - [ 3, 0, 0 ]
-  - [ 2, 4, 0 ]
-  - [ 0, 0, 1 ]
-PrimativeCellBasisNumber: [ 8, 8, 8 ]
-AtomPositionBySuperCell:
-  - [ 0.9999442136878602, 0.1668131534569660, 0.2495178577809962 ]
-  - [ 0.1666666570000004, 0.0821651768612523, 0.7494996261803192 ]
-  - [ 0.0000785931726128, 0.9991985635426676, 0.9980828222129586 ]
-  - [ 0.9998667686594940, 0.0010323603084980, 0.4990245675767412 ]
-  - [ 0.8333163913323451, 0.4170889259273483, 0.2509448387496841 ]
-  - [ 0.0010593515672643, 0.3339499564673934, 0.7504013581078407 ]
-  - [ 0.8333234221246683, 0.2494366873144543, 0.9998086408888811 ]
-  - [ 0.8335033186823803, 0.2508241081996696, 0.5000431170293115 ]
-  - [ 0.6666666659999976, 0.6667639132151945, 0.2493250038919751 ]
-  - [ 0.8328430132025012, 0.5824561773357575, 0.7502934292698500 ]
-  - [ 0.6666666670000012, 0.4993068721873999, 0.0001383770014351 ]
-  - [ 0.6666666670000012, 0.5005472005797680, 0.5002523804663008 ]
-  - [ 0.4999087106685396, 0.9169692013808204, 0.2492033506084066 ]
-  - [ 0.6666666570000004, 0.8337032368714793, 0.7504363269561648 ]
-  - [ 0.4995637733746955, 0.7493807863038214, 0.9991435901176927 ]
-  - [ 0.5003725611872856, 0.7503863239197331, 0.4992410940474822 ]
-  - [ 0.3333891193121384, 0.1668131534569660, 0.2495178577809962 ]
-  - [ 0.4992733589909116, 0.0818752444180838, 0.7497937354374519 ]
-  - [ 0.3332547398273861, 0.9991985635426676, 0.9980828222129586 ]
-  - [ 0.3334665643405045, 0.0010323603084980, 0.4990245675767412 ]
-  - [ 0.1666666659999976, 0.4169480526102884, 0.2506028916346738 ]
-  - [ 0.3322739714327337, 0.3339499564673934, 0.7504013581078407 ]
-  - [ 0.1666666670000012, 0.2497373319249146, 0.9999500604699707 ]
-  - [ 0.1666666670000012, 0.2505366718958373, 0.4997757704091679 ]
-  - [ 0.0007343296568028, 0.6675813991651439, 0.2507455593024008 ]
-  - [ 0.1666666570000004, 0.5882597899202560, 0.7580329988559029 ]
-  - [ 0.9993230602783540, 0.4989432733446046, 0.0018414165511462 ]
-  - [ 0.0010383420437462, 0.5016621987383314, 0.5008853395289279 ]
-  - [ 0.8334246223314590, 0.9169692013808204, 0.2492033506084066 ]
-  - [ 0.0000682038488344, 0.8265801375506314, 0.7450704126627182 ]
-  - [ 0.8337695596253030, 0.7493807863038214, 0.9991435901176927 ]
-  - [ 0.8329607718127133, 0.7503863239197331, 0.4992410940474822 ]
-  - [ 0.6666666659999976, 0.1670868510386173, 0.2499476084075379 ]
-  - [ 0.8340599550090892, 0.0818752444180838, 0.7497937354374519 ]
-  - [ 0.6666666670000012, 0.9994774515161969, 0.9997481863090192 ]
-  - [ 0.6666666670000012, 0.0005098915308569, 0.4992736940964735 ]
-  - [ 0.5000169416676533, 0.4170889259273483, 0.2509448387496841 ]
-  - [ 0.6666666570000004, 0.3320791225644938, 0.7499985648946081 ]
-  - [ 0.5000099108753309, 0.2494366873144543, 0.9998086408888811 ]
-  - [ 0.4998300143176183, 0.2508241081996696, 0.5000431170293115 ]
-  - [ 0.3325990033431961, 0.6675813991651439, 0.2507455593024008 ]
-  - [ 0.5004902997974959, 0.5824561773357575, 0.7502934292698500 ]
-  - [ 0.3340102727216447, 0.4989432733446046, 0.0018414165511462 ]
-  - [ 0.3322949909562527, 0.5016621987383314, 0.5008853395289279 ]
-  - [ 0.1666666659999976, 0.9180544950396436, 0.2478666411127934 ]
-  - [ 0.3332651191511635, 0.8265801375506314, 0.7450704126627182 ]
-  - [ 0.1666666670000012, 0.7492229363186385, 0.0001611773946995 ]
-  - [ 0.1666666670000012, 0.7586221208203656, 0.5025559614474544 ]
-  - [ 0.9996456144438632, 0.1676511882651485, 0.4376401280043795 ]
-  - [ 0.1666666570000004, 0.0828996116757371, 0.9371848303263415 ]
-  - [ 0.9996408690505139, 0.0000023801845120, 0.1856212094249273 ]
-  - [ 0.0015663355615017, 0.9973672273519564, 0.6863958225036811 ]
-  - [ 0.8337845559136262, 0.4178115825074205, 0.4389674348305436 ]
-  - [ 0.9998976961353622, 0.3327868228709865, 0.9381885158054116 ]
-  - [ 0.8333373011703021, 0.2502176615540846, 0.1873187173259448 ]
-  - [ 0.8342772552129253, 0.2495543040285077, 0.6873480333936619 ]
-  - [ 0.6666666659999976, 0.6669417550050203, 0.4369673987476388 ]
-  - [ 0.8333807305605515, 0.5824145944486363, 0.9381784129487265 ]
-  - [ 0.6666666670000012, 0.4997745444284661, 0.1874992047550911 ]
-  - [ 0.6666666670000012, 0.4986742269313283, 0.6876164226427359 ]
-  - [ 0.4999339071178296, 0.9171919757061141, 0.4373995923163460 ]
-  - [ 0.6666666570000004, 0.8325907836040870, 0.9377955619721441 ]
-  - [ 0.4994945329888917, 0.7500616410188495, 0.1869506467429709 ]
-  - [ 0.5026656836615087, 0.7480725752499672, 0.6860625197385234 ]
-  - [ 0.3336877185561356, 0.1676511882651485, 0.4376401280043795 ]
-  - [ 0.5000447115614929, 0.0824669311545484, 0.9373067134583003 ]
-  - [ 0.3336924639494850, 0.0000023801845120, 0.1856212094249273 ]
-  - [ 0.3317669974384970, 0.9973672273519564, 0.6863958225036811 ]
-  - [ 0.1666666659999976, 0.4174122540170955, 0.4385569558692831 ]
-  - [ 0.3334356268646358, 0.3327868228709865, 0.9381885158054116 ]
-  - [ 0.1666666670000012, 0.2500772537096505, 0.1873263641549840 ]
-  - [ 0.1666666670000012, 0.2496312870965579, 0.6872470991449243 ]
-  - [ 0.0013188200046441, 0.6685626427893232, 0.4398333944722485 ]
-  - [ 0.1666666570000004, 0.5809667172472394, 0.9434672115062485 ]
-  - [ 0.9997935800928740, 0.5002787858188761, 0.1889409420456360 ]
-  - [ 0.0015430640609559, 0.5024439923185628, 0.6881853863336864 ]
-  - [ 0.8333994258821685, 0.9171919757061141, 0.4373995923163460 ]
-  - [ 0.0009625351257081, 0.8317674754963278, 0.9354088923528614 ]
-  - [ 0.8338388000111076, 0.7500616410188495, 0.1869506467429709 ]
-  - [ 0.8306676493384901, 0.7480725752499672, 0.6860625197385234 ]
-  - [ 0.6666666659999976, 0.1675291153356092, 0.4379848675305987 ]
-  - [ 0.8332886024385080, 0.0824669311545484, 0.9373067134583003 ]
-  - [ 0.6666666670000012, 0.0001944358287058, 0.1870538575713719 ]
-  - [ 0.6666666670000012, 0.9997639879549629, 0.6867185860979695 ]
-  - [ 0.4995487770863726, 0.4178115825074205, 0.4389674348305436 ]
-  - [ 0.6666666570000004, 0.3324461382685149, 0.9378848515350188 ]
-  - [ 0.4999960318296968, 0.2502176615540846, 0.1873187173259448 ]
-  - [ 0.4990560777870729, 0.2495543040285077, 0.6873480333936619 ]
-  - [ 0.3320145129953548, 0.6685626427893232, 0.4398333944722485 ]
-  - [ 0.4999525824394457, 0.5824145944486363, 0.9381784129487265 ]
-  - [ 0.3335397529071246, 0.5002787858188761, 0.1889409420456360 ]
-  - [ 0.3317902689390430, 0.5024439923185628, 0.6881853863336864 ]
-  - [ 0.1666666659999976, 0.9216121818138033, 0.4343691564472280 ]
-  - [ 0.3323707878742898, 0.8317674754963278, 0.9354088923528614 ]
-  - [ 0.1666666670000012, 0.7508040359125291, 0.1873373990293362 ]
-QpointDataInputFile: band.hdf5
-QpointDataOutputFile:
-  - Filename: unfold-output.zpp
--- a/example/step4/unfold-output.zpp
+++ b/example/step4/unfold-output.zpp
--- a/example/step5/band.hdf5
+++ b/example/step5/band.hdf5
--- a/example/step5/config.yaml
+++ b/example/step5/config.yaml
@@ -1,13 +0,0 @@
-UnfoldedDataFile: unfold-output.zpp
-Qpoints:
-  -
-    - [0.0, 0.0, 0.0]
-    - [0.5, 0.0, 0.0]
-Resolution:
-  X: 64
-  Y: 256
-FrequencyRange:
-  Min: 0
-  Max: 40
-OutputPictureFile: unfold-output.png
-OutputDataFile: plot-output.hdf5
--- a/example/step5/unfold-output.hdf5
+++ b/example/step5/unfold-output.hdf5
--- a/example/step5/unfold-output.png
+++ b/example/step5/unfold-output.png
--- a/example/step5/unfold-output.zpp
+++ b/example/step5/unfold-output.zpp
--- a/flake.lock
+++ b/flake.lock
--- a/flake.nix
+++ b/flake.nix
@@ -1,25 +0,0 @@
-{
-  inputs.nixos.url = "github:CHN-beta/nixos";
-  inputs.nixpkgs.url = "github:NixOS/nixpkgs/nixos-unstable";
-
-  outputs = inputs:
-    let
-      # pkgs = inputs.nixpkgs.legacyPackages.x86_64-linux;
-      pkgs = import inputs.nixpkgs
-      {
-        localSystem = { system = "x86_64-linux"; gcc = { arch = "alderlake"; tune = "alderlake"; }; };
-        config.allowUnfree = true;
-      };
-      localPackages = import "${inputs.nixos}/local/pkgs" { inherit pkgs; inherit (inputs.nixpkgs) lib; };
-    in
-    {
-      devShell.x86_64-linux = pkgs.mkShell.override { stdenv = pkgs.gcc13Stdenv; }
-      {
-        packages = with pkgs; [ pkg-config cmake ninja ];
-        buildInputs = (with pkgs; [ eigen yaml-cpp fmt highfive tbb_2021_8.dev glfw libGL ])
-          ++ (with localPackages; [ concurrencpp matplotplusplus zpp-bits glad ]);
-        hardeningDisable = [ "all" ];
-        # NIX_DEBUG = "1";
-      };
-    };
-}
--- a/include/ufo.hpp
+++ b/include/ufo.hpp
@@ -0,0 +1,104 @@
+# pragma once
+# include <biu.hpp>
+
+namespace ufo
+{
+  // 在相位中, 约定为使用 $\exp (2 \pi i \vec{q} \cdot \vec{r})$ 来表示原子的运动状态
+  //  (而不是 $\exp (-2 \pi i \vec{q} \cdot \vec{r})$)
+  // 一些书定义的倒格矢中包含了 $2 \pi$ 的部分, 我们这里约定不包含这部分.
+  //  也就是说, 正格子与倒格子的转置相乘, 得到单位矩阵.
+
+  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;
+      std::optional<std::array<double, 3>> AtomTranslation;
+      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<7>;
+    } 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 project_to_mode(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::string path);
+  void raman_apply_contribution(std::string config_file);
+}
--- a/include/ufo/fold.hpp
+++ b/include/ufo/fold.hpp
@@ -1,37 +0,0 @@
-# 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
-      );
-  };
-}
--- a/include/ufo/plot.hpp
+++ b/include/ufo/plot.hpp
@@ -1,77 +0,0 @@
-# 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
-      );
-  };
-}
--- a/include/ufo/solver.hpp
+++ b/include/ufo/solver.hpp
@@ -1,141 +0,0 @@
-# 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)
--- a/include/ufo/unfold.hpp
+++ b/include/ufo/unfold.hpp
@@ -1,164 +0,0 @@
-# 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
-      );
-  };
-}
--- a/src/fold.cpp
+++ b/src/fold.cpp
@@ -1,86 +1,56 @@
-# include <ufo/fold.hpp>
+# include <ufo.hpp>

-namespace ufo
+void ufo::fold(std::string config_file)
 {
-  FoldSolver::InputType::InputType(std::string config_file)
+  struct Config
  {
-    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();
-    }();
-  }
+    Eigen::Matrix3d SuperCellDeformation;
+    Eigen::Vector3i SuperCellMultiplier;
+    std::vector<Eigen::Vector3d> Qpoints;
+  };

-  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
+  auto 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
-  )
+    Eigen::Matrix3d super_cell_deformation,
+    Eigen::Vector3i super_cell_multiplier,
+    Eigen::Vector3d qpoint_in_reciprocal_primitive_cell_by_reciprocal_primitive_cell
+  ) -> Eigen::Vector3d
  {
-    // 首先需要将 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
+      首先需要将 q 点坐标的单位转换为 ModifiedSuperCell 的格矢，可知：
+      QpointByReciprocalModifiedSuperCell =
+        SuperCellMultiplier.cast<double>().asDiagonal() * QpointByReciprocalPrimitiveCell;
+      接下来考虑将 q 点坐标的单位转换为 SuperCell 的格矢
+      ModifiedSuperCell = SuperCellMultiplier.transpose() * PrimativeCell;
+      SuperCell = SuperCellDeformation * ModifiedSuperCell;
+      ReciprocalModifiedSuperCell = ModifiedSuperCell.inverse().transpose();
+      ReciprocalSuperCell = SuperCell.inverse().transpose();
+      Qpoint = QpointByReciprocalModifiedSuperCell.transpose() * ReciprocalModifiedSuperCell;
+      Qpoint = QpointByReciprocalSuperCell.transpose() * ReciprocalSuperCell;
      整理可以得到:
-      QpointByReciprocalSuperCell = SuperCellDeformation * QpointByReciprocalModifiedSuperCell
+      QpointByReciprocalSuperCell = SuperCellDeformation * QpointByReciprocalModifiedSuperCell;
+      两个式子结合，可以得到：
+      QpointByReciprocalSuperCell =
+        SuperCellDeformation * SuperCellMultiplier.cast<double>().asDiagonal() * QpointByReciprocalPrimitiveCell;
    */
-    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;
-  }
+    auto qpoint_in_reciprocal_primitive_cell_by_reciprocal_super_cell =
+    (
+      super_cell_deformation * super_cell_multiplier.cast<double>().asDiagonal()
+        * qpoint_in_reciprocal_primitive_cell_by_reciprocal_primitive_cell
+    ).eval();
+    /*
+      到目前为止，我们还没有移动过 q 点的坐标。现在，我们将它移动整数个 ReciprocalSuperCell，直到它落在超胞的倒格子中。
+      这等价于直接取 QpointByReciprocalSuperCell - QpointByReciprocalSuperCell.floor()。
+    */
+    return (qpoint_in_reciprocal_primitive_cell_by_reciprocal_super_cell.array()
+      - qpoint_in_reciprocal_primitive_cell_by_reciprocal_super_cell.array().floor()).matrix();
+  };
+
+  biu::Logger::Guard log(config_file);
+  auto input = YAML::LoadFile(config_file).as<Config>();
+  for (const auto& qpoint : input.Qpoints) log.info("{} -> {}"_f
+    (
+      qpoint,
+      fold(input.SuperCellDeformation, input.SuperCellMultiplier, qpoint)
+    ));
 }
--- a/src/main.cpp
+++ b/src/main.cpp
@@ -1,17 +1,18 @@
-# include <ufo/fold.hpp>
-# include <ufo/unfold.hpp>
-# include <ufo/plot.hpp>
+# include <ufo.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]));
+  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] == "project-to-mode"s) ufo::project_to_mode(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(argv[2]);
+  else if (argv[1] == "raman-apply-contribution"s) ufo::raman_apply_contribution(argv[2]);
+  else if (argv[1] == "plot-band"s) ufo::plot_band(argv[2]);
+  else if (argv[1] == "plot-point"s) ufo::plot_point(argv[2]);
+  else throw std::runtime_error("Unknown task: {}"_f(argv[1]));
 }
--- a/src/plot.cpp
+++ b/src/plot.cpp
@@ -1,233 +1,179 @@
-# include <ufo/plot.hpp>
+# include <ufo.hpp>
+# include <matplot/matplot.h>
+# include <boost/container/flat_map.hpp>

-namespace ufo
+void ufo::plot_band(std::string config_file)
 {
-  PlotSolver::InputType::UnfoldedDataType::UnfoldedDataType(std::string filename)
+  struct Config
  {
-    static_cast<UnfoldSolver::OutputType&>(*this) = zpp_read<UnfoldSolver::OutputType>(filename);
-  }
+    std::string InputDataFile;
+    // 要画图的 q 点路径列表
+    // 内层表示一个路径上的 q 点，外层表示不同的路径
+    // 单位为倒格矢
+    std::vector<std::vector<Eigen::Vector3d>> Qpoints;
+    // 插值时使用的分辨率（不影响画出来图片的分辨率和横纵比）
+    std::array<std::size_t, 2> InterpolationResolution;
+    // 画图区域的y轴和x轴的比例。如果不指定，则由matplot++自动调整（通常调整为正方形，即 1）
+    std::optional<double> AspectRatio;
+    // 整张图片的分辨率
+    std::optional<std::array<std::size_t, 2>> PictureResolution;
+    // 画图的频率范围
+    std::array<double, 2> FrequencyRange;
+    // 搜索 q 点时的阈值，单位为埃^-1
+    std::optional<double> ThresholdWhenSearchingQpoints;
+    // 是否要在 y 轴上作一些标记
+    std::optional<std::vector<std::pair<double, std::string>>> YTicks;
+    // 是否输出图片
+    std::optional<std::string> OutputPictureFile;
+    // 是否输出数据，可以进一步使用 matplotlib 画图
+    std::optional<std::string> OutputDataFile;
+  };

-  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
+  // 根据 q 点路径, 搜索要使用的 q 点，返回的是超胞中 Qpint 的索引和 SubQpoint 的索引，以及到路径起点的距离，以及这段路径的总长度
+  auto search_qpoints = []
  (
+    const Eigen::Matrix3d& primative_cell,
    const std::pair<Eigen::Vector3d, Eigen::Vector3d>& path,
-    const decltype(InputType::UnfoldedDataType::QpointData)& available_qpoints,
-    double threshold, bool exclude_endpoint
+    const std::vector<std::vector<Eigen::Vector3d>>& qpoints,
+    double threshold, bool exclude_endpoint = false
  )
  {
-    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)
+    // 键为这个点到起点的距离
+    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 j = 0; j < qpoints[i].size(); j++)
+      for (auto cell_shift
+        : biu::sequence(Eigen::Vector3i(-1, -1, -1), Eigen::Vector3i(2, 2, 2)))
      {
-        // 计算这个点到前两个点的距离, 两个距离都应该小于两点之间的距离
-        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));
+        Eigen::Vector3d qpoint
+          = 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();
+        // 如果这个点到前两个点所在直线的距离小于阈值, 则认为这个点在这条直线上，但不一定在这两个点之间
+        if (distance < threshold)
+        {
+          // 计算这个点到前两个点的距离, 两个距离都应该小于两点之间的距离
+          auto distance1 = (qpoint - begin).norm();
+          auto distance2 = (qpoint - end).norm();
+          auto distance3 = (end - begin).norm();
+          if (distance1 < distance3 + threshold && distance2 < distance3 + threshold)
+            // 如果这个点不在终点处, 或者不排除终点, 则加入
+            if (distance2 > threshold || !exclude_endpoint)
+              selected_qpoints.emplace(distance1, std::pair{i, j});
+        }
      }
-    }
    // 去除非常接近的点
    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 (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;
-  }
+    if (selected_qpoints.empty()) throw std::runtime_error("No q points found");
+    return std::make_pair(selected_qpoints, (end - begin).norm());
+  };

-  std::tuple<std::vector<std::vector<double>>, std::vector<double>> PlotSolver::calculate_values
+  // 根据搜索到的 q 点, 计算图中每个点的值
+  auto 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
+    // search_qpoints 的第一个返回值
+    const boost::container::flat_map<double, std::pair<std::size_t, std::size_t>>& path,
+    // 每一条连续路径的第一个 q 点的索引
+    const std::set<std::size_t>& path_begin,
+    // 各个模式的频率和权重，几个维度分别为：MetaQpointIndex, ModeIndex
+    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,
+    // 路径的总长度
+    double total_distance
  )
  {
-    // 整理输入
-    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 = []
+    // 按比例混合两个 q 点的结果，得到可以用于画图的那一列数据
+    auto blend = [&]
    (
-      const UnfoldSolver::OutputType::QpointDataType& a,
-      const UnfoldSolver::OutputType::QpointDataType& b,
-      double ratio, unsigned resolution, std::pair<double, double> range
+      // 两个点的索引
+      std::pair<std::size_t, std::size_t> a, std::pair<std::size_t, std::size_t> b,
+      // 按照连续路径混合还是按照断开的路径混合
+      bool continuous,
+      // 第一个点占的比例
+      double ratio,
+      std::size_t resolution, std::array<double, 2> frequency_range
    ) -> std::vector<double>
    {
-      // 计算插值结果
-      std::vector<double> frequency, weight;
-      for (unsigned i = 0; i < a.ModeData.size(); i++)
+      // 混合得到的频率和权重
+      std::vector<double> frequency_result, weight_result;
+      // 如果是连续路径，将每个模式的频率和权重按照比例混合
+      if (continuous)
      {
-        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));
+        for (std::size_t i = 0; i < frequency[a.first].size(); i++)
+        {
+          frequency_result.push_back(frequency[a.first][i] * ratio + frequency[b.first][i] * (1 - ratio));
+          weight_result.push_back(weight[a.first][i][a.second] * ratio + weight[b.first][i][b.second] * (1 - ratio));
+        }
+      }
+      // 如果是不连续路径，将每个模式的权重乘以比例，最后相加
+      else
+      {
+        for (std::size_t i = 0; i < frequency[a.first].size(); i++)
+        {
+          frequency_result.push_back(frequency[a.first][i]);
+          weight_result.push_back(weight[a.first][i][a.second] * ratio);
+        }
+        for (std::size_t i = 0; i < frequency[b.first].size(); i++)
+        {
+          frequency_result.push_back(frequency[b.first][i]);
+          weight_result.push_back(weight[b.first][i][b.second] * (1 - ratio));
+        }
      }
      std::vector<double> result(resolution);
-      for (unsigned i = 0; i < frequency.size(); i++)
+      for (std::size_t i = 0; i < frequency_result.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];
+        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_result[i];
      }
      return result;
    };
-    for (unsigned i = 0; i < resolution.first; i++)
+
+    std::vector<std::vector<double>> values;
+    for (std::size_t i = 0; i < resolution[0]; 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)
-        );
+      auto current_distance = total_distance * i / resolution[0];
+      auto it = path.lower_bound(current_distance);
+      if (it == path.begin()) values.push_back(blend
+          (it->second, it->second, true, 1, resolution[1], frequency_range));
+      else if (it == path.end()) values.push_back(blend
+      (
+        std::prev(it)->second, std::prev(it)->second, true, 1,
+        resolution[1], frequency_range
+      ));
+      else values.push_back(blend
+      (
+        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
+      ));
    }
-    return {values, x_ticks};
-  }
-  void PlotSolver::plot
+    return values;
+  };
+
+  // 根据数值, 画图
+  auto plot = []
  (
    const std::vector<std::vector<double>>& values,
    const std::string& filename,
-    const std::vector<double>& x_ticks, const std::vector<double>& y_ticks
+    const std::vector<double>& x_ticks, const std::vector<double>& y_ticks,
+    const std::vector<std::string>& y_ticklabels,
+    const std::optional<double>& aspect_ratio,
+    const std::optional<std::array<std::size_t, 2>>& resolution
  )
  {
    std::vector<std::vector<double>>
@@ -235,32 +181,285 @@ namespace ufo
      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++)
+    for (std::size_t i = 0; i < values[0].size(); i++)
+      for (std::size_t j = 0; j < values.size(); j++)
      {
        auto v = values[j][i];
-        if (v < 0.05)
-          v = 0;
+        if (v < 0.05) v = 0;
        a[i][j] = v * 100 * 255;
-        if (a[i][j] > 255)
-          a[i][j] = 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;
+        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;
+        if (g[i][j] < 0) g[i][j] = 0;
        b[i][j] = 255;
      }
-    auto f = matplot::figure<matplot::backend::gnuplot>(true);
+    auto f = matplot::figure(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->x_axis().ticklabels(std::vector<std::string>(x_ticks.size()));
    ax->y_axis().tick_values(y_ticks);
    ax->y_axis().tick_length(1);
+    ax->y_axis().ticklabels(y_ticklabels);
+    if (aspect_ratio)
+    {
+      ax->axes_aspect_ratio_auto(false);
+      ax->axes_aspect_ratio(*aspect_ratio);
+    }
+    if (resolution)
+    {
+      f->width((*resolution)[0]);
+      f->height((*resolution)[1]);
+    }
    f->save(filename, "png");
+  };
+
+  auto config = YAML::LoadFile(config_file).as<Config>();
+  auto input = biu::deserialize<CommonData>
+    (biu::read<std::byte>(config.InputDataFile));
+  
+  // 搜索画图需要用到的 q 点
+  // key 到起点的距离，value 为 q 点在 QpointData 中的索引
+  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 : config.Qpoints)
+  {
+    assert(line.size() >= 2);
+    path_begin.insert(path.size());
+    for (std::size_t i = 0; i < line.size() - 1; i++)
+    {
+      x_ticks_index.insert(path.size());
+      auto [this_path, this_distance] = search_qpoints
+      (
+        input.Primative.Cell, {line[i], line[i + 1]},
+        input.Super.Qpoint
+          | ranges::views::transform([](auto& qpoint) { return qpoint.SubQpoint; })
+          | ranges::to_vector,
+        config.ThresholdWhenSearchingQpoints.value_or(0.001),
+        i != line.size() - 2
+      );
+      path.merge
+      (
+        this_path
+        | ranges::views::transform([&](auto& p)
+          { return std::make_pair(p.first + total_distance, p.second); })
+        | ranges::to<boost::container::flat_map>
+      );
+      total_distance += this_distance;
+    }
  }
+
+  // 计算画图的数据
+  auto values = calculate_values
+  (
+    path, path_begin,
+    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 * config.InterpolationResolution[0]; }) | ranges::to_vector;
+  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 - config.FrequencyRange[0]) / (config.FrequencyRange[1] - config.FrequencyRange[0])
+        * config.InterpolationResolution[1];
+    })
+    | ranges::to_vector;
+  auto y_ticklabels = config.YTicks.value_or(std::vector<std::pair<double, std::string>>{})
+    | biu::toLvalue | ranges::views::values | ranges::to_vector;
+  if (config.OutputPictureFile) plot
+  (
+    values, config.OutputPictureFile.value(),
+    x_ticks, y_ticks, y_ticklabels, config.AspectRatio, config.PictureResolution
+  );
+  if (config.OutputDataFile)
+    biu::Hdf5file(config.OutputDataFile.value(), true)
+      .write("Values", values)
+      .write("XTicks", x_ticks)
+      .write("YTicks", y_ticks)
+      .write("YTickLabels", y_ticklabels)
+      .write("InterpolationResolution", config.InterpolationResolution)
+      .write("FrequencyRange", config.FrequencyRange);
+}
+
+void ufo::plot_point(std::string config_file)
+{
+  struct Config
+  {
+    std::string UnfoldedDataFile;
+    // 要画图的 q 点
+    Eigen::Vector3d Qpoint;
+    // 插值的分辨率
+    std::size_t InterpolationResolution;
+    std::optional<double> AspectRatio;
+    std::optional<std::array<std::size_t, 2>> PictureResolution;
+    // 画图的频率范围
+    std::array<double, 2> FrequencyRange;
+    // 搜索 q 点时的阈值，单位为埃^-1
+    std::optional<double> ThresholdWhenSearchingQpoints;
+    // 是否要在 z 轴上作一些标记
+    std::optional<std::vector<std::pair<double, std::string>>> XTicks;
+    // 是否输出图片
+    std::optional<std::string> OutputPictureFile;
+    // 是否输出插值后数据，可以进一步使用 matplotlib 画图
+    std::optional<std::string> OutputDataFile;
+    // 是否输出插值前数据，可以配合 phonopy 结果深入研究
+    std::optional<std::string> OutputRawDataFile;
+  };
+
+  // 根据 q 点路径, 搜索要使用的 q 点，返回的是 q 点在 QpointData 中的索引
+  auto search_qpoints = []
+  (
+    const Eigen::Matrix3d& primative_cell,
+    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 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][j] + cell_shift.first.cast<double>())).eval();
+        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");
+  };
+
+  // 根据搜索到的 q 点, 计算图中每个点的值
+  auto calculate_values = []
+  (
+    // q 点的数据（需要用到它的频率和权重）
+    const std::vector<double>& frequency,
+    const std::vector<double>& weight,
+    // 用于插值的分辨率和范围
+    std::size_t resolution,
+    const std::array<double, 2>& frequency_range
+  )
+  {
+    biu::Logger::Guard log;
+    std::vector<double> result(resolution);
+    for (std::size_t i = 0; i < frequency.size(); i++)
+    {
+      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] += weight[i];
+    }
+    return log.rtn(result);
+  };
+
+  // 根据数值, 画图
+  auto plot = []
+  (
+    const std::vector<double>& values, const std::string& filename,
+    const std::vector<double>& x_ticks, const std::vector<std::string>& x_ticklabels,
+    const std::optional<double>& aspect_ratio, const std::optional<std::array<std::size_t, 2>>& resolution
+  )
+  {
+    biu::Logger::Guard log;
+    auto f = matplot::figure(true);
+    auto ax = f->current_axes();
+    auto image = ax->area(values, 0, false, "");
+    ax->y_axis().reverse(false);
+    ax->x_axis().tick_values(x_ticks);
+    ax->x_axis().tick_length(1);
+    ax->x_axis().ticklabels(x_ticklabels);
+    ax->y_axis().tick_values({});
+    if (aspect_ratio)
+    {
+      ax->axes_aspect_ratio_auto(false);
+      ax->axes_aspect_ratio(*aspect_ratio);
+    }
+    if (resolution)
+    {
+      f->width((*resolution)[0]);
+      f->height((*resolution)[1]);
+    }
+    f->save(filename, "png");
+  };
+
+  biu::Logger::Guard log;
+  auto config = YAML::LoadFile(config_file).as<Config>();
+  auto input = biu::deserialize<CommonData>
+    (biu::read<std::byte>(config.UnfoldedDataFile));
+  
+  auto qpoint_index = search_qpoints
+  (
+    input.Primative.Cell, config.Qpoint,
+    input.Super.Qpoint
+      | ranges::views::transform([](auto& qpoint) { return qpoint.SubQpoint; })
+      | ranges::to_vector,
+    config.ThresholdWhenSearchingQpoints.value_or(0.001)
+  );
+  auto values = calculate_values
+  (
+    input.Super.Qpoint[qpoint_index.first].Mode
+      | 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 = config.XTicks.value_or(std::vector<std::pair<double, std::string>>{})
+    | biu::toLvalue | ranges::views::keys
+    | ranges::views::transform([&](auto i)
+    {
+      return (i - config.FrequencyRange[0]) / (config.FrequencyRange[1] - config.FrequencyRange[0])
+        * config.InterpolationResolution;
+    })
+    | ranges::to_vector;
+  auto x_ticklabels = config.XTicks.value_or(std::vector<std::pair<double, std::string>>{})
+    | biu::toLvalue | ranges::views::values | ranges::to_vector;
+  if (config.OutputPictureFile) plot
+  (
+    values, config.OutputPictureFile.value(),
+    x_ticks, x_ticklabels, config.AspectRatio, config.PictureResolution
+  );
+  if (config.OutputDataFile)
+    biu::Hdf5file(config.OutputDataFile.value(), true)
+      .write("Values", values)
+      .write("XTicks", x_ticks)
+      .write("XTickLabels", x_ticklabels)
+      .write("InterpolationResolution", config.InterpolationResolution)
+      .write("FrequencyRange", config.FrequencyRange);
+  // TODO: pfr + optional + Eigen Matrix = failed
+  if (config.OutputRawDataFile)
+    std::ofstream(*config.OutputRawDataFile) << YAML::Node(values);
 }
--- a/src/project-to-atom.cpp
+++ b/src/project-to-atom.cpp
@@ -0,0 +1,36 @@
+# 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().squaredNorm().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
+        ));
+}
--- a/src/project-to-mode.cpp
+++ b/src/project-to-mode.cpp
@@ -0,0 +1,130 @@
+# include <ufo.hpp>
+
+void ufo::project_to_mode(std::string config_file)
+{
+  // 将选定的超胞中的某一个 q 点上的每一个模式到反折叠后的某个 q 点的投影，投射到单胞中的对应 q 点的模式上
+  // 与 unfold 的方法类似，只有一个重要的区别：
+  // 组成基矢的平面波的波矢的两个部分中，第二个部分不再需要。因为这个部分的结果已经算出来了，只需要将结果乘上 WeightOnUnfold 即可。
+  struct Config
+  {
+    // 在单胞内取几个平面波的基矢
+    Eigen::Vector3i PrimativeCellBasisNumber;
+    // 要选定的超胞中的 q 点，sub qpoint，以及单胞中的 q 点
+    std::size_t SuperQpointIndex, SubQpointIndex, PrimativeQpointIndex;
+    // 输入输出文件的路径，直接输出 yaml，因为这些数据不再需要后处理
+    std::string InputDataFile, OutputFile;
+  };
+  struct Output
+  {
+    Eigen::Vector3d SuperQpoint, SubQpoint, PrimativeQpoint;
+    std::vector<std::vector<double>> Coefficient;
+  };
+
+  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));
+  Output output
+  {
+    .SuperQpoint = input.Super.Qpoint[config.SuperQpointIndex].Qpoint,
+    .SubQpoint = input.Super.Qpoint[config.SuperQpointIndex].SubQpoint[config.SubQpointIndex],
+    .PrimativeQpoint = input.Primative.Qpoint[config.PrimativeQpointIndex].Qpoint,
+    .Coefficient = std::vector<std::vector<double>>(input.Super.Qpoint[config.SuperQpointIndex].Mode.size(),
+      std::vector<double>(input.Primative.Qpoint[config.PrimativeQpointIndex].Mode.size()))
+  };
+  if ((output.SubQpoint - output.PrimativeQpoint).norm() > 0.01)
+    log.error("sub qpoint {} != primative qpoint {}"_f(output.SubQpoint, output.PrimativeQpoint));
+
+  // 基在原胞中的表示
+  auto primative_basis = [&]
+  {
+    biu::Logger::Guard log;
+    std::vector<Eigen::VectorXcd> basis((config.PrimativeCellBasisNumber.array() * 2 - 1).prod());
+    for (auto [xyz_of_basis, i_of_basis] : biu::sequence
+      ((-config.PrimativeCellBasisNumber + Eigen::Vector3i::Ones()).eval(), config.PrimativeCellBasisNumber))
+    {
+      // 波矢就是 xyz_of_basis, 单位为单胞的倒格矢
+      // 将它的单位转换成埃^-1
+      auto wavevector =
+        (xyz_of_basis.transpose().cast<double>() * input.Primative.Cell.inverse().transpose()).transpose();
+      // 将原子坐标单位也转换为埃
+      auto atom_position = input.Primative.AtomPosition * input.Primative.Cell;
+      // 计算基矢
+      basis[i_of_basis] = (2i * std::numbers::pi_v<double> * (atom_position * wavevector)).array().exp();
+    }
+    return basis;
+  }();
+  // 基在超胞中的表示
+  auto super_basis = [&]
+  {
+    biu::Logger::Guard log;
+    std::vector<Eigen::VectorXcd> basis((config.PrimativeCellBasisNumber.array() * 2 - 1).prod());
+    auto diff_of_sub_qpoint_by_reciprocal_modified_super_cell = [&]
+    {
+      for
+      (
+        auto [diff_of_sub_qpoint_by_reciprocal_modified_super_cell, i_of_sub_qpoint]
+          : biu::sequence(input.Super.CellMultiplier)
+      )
+        if (i_of_sub_qpoint == config.SubQpointIndex) return diff_of_sub_qpoint_by_reciprocal_modified_super_cell;
+      std::unreachable();
+    }();
+    for (auto [xyz_of_basis, i_of_basis] : biu::sequence
+      ((-config.PrimativeCellBasisNumber + Eigen::Vector3i::Ones()).eval(), config.PrimativeCellBasisNumber))
+    {
+      // 计算波矢, 单位为单胞的倒格矢
+      auto wavevector_by_reciprocal_primative_cell = xyz_of_basis.cast<double>()
+        + input.Super.CellMultiplier.cast<double>().cwiseInverse().asDiagonal()
+        * diff_of_sub_qpoint_by_reciprocal_modified_super_cell.cast<double>();
+      // 将单位转换为埃^-1
+      auto wavevector =
+        (wavevector_by_reciprocal_primative_cell.transpose() * input.Primative.Cell.inverse().transpose()).transpose();
+      // 将原子坐标单位也转换为埃
+      auto atom_translation =
+        (input.Super.AtomTranslation.value_or(std::array{0., 0., 0.}) | biu::toEigen<>)
+        .transpose().replicate(input.Super.AtomPosition.rows(), 1).eval();
+      auto atom_position = (input.Super.AtomPosition - atom_translation)
+        * (input.Super.CellDeformation * input.Super.CellMultiplier.cast<double>().asDiagonal() * input.Primative.Cell);
+      // 计算基矢
+      basis[i_of_basis] = (2i * std::numbers::pi_v<double> * (atom_position * wavevector)).array().exp();
+    }
+    return basis;
+  }();
+  // 将所有模式投影到基矢上
+  auto primative_projection_coefficient =
+    input.Primative.Qpoint[config.PrimativeQpointIndex].Mode
+      | ranges::views::transform([&](const auto& mode)
+        {
+          return primative_basis
+            | ranges::views::transform([&](const auto& basis)
+              { return (basis.transpose().conjugate() * mode.EigenVector).eval() | biu::fromEigen; })
+            | ranges::to_vector
+            | biu::toEigen<>;
+        })
+      | ranges::to_vector;
+  auto super_projection_coefficient =
+    input.Super.Qpoint[config.SuperQpointIndex].Mode
+      | ranges::views::transform([&](const auto& mode)
+        {
+          return super_basis
+            | ranges::views::transform([&](const auto& basis)
+              { return (basis.transpose().conjugate() * mode.EigenVector).eval() | biu::fromEigen; })
+            | ranges::to_vector
+            | biu::toEigen<>;
+        })
+      | ranges::to_vector;
+  // 将它们互相投影，得到结果
+  for (auto [i_of_super_mode, super_mode] : ranges::views::enumerate(super_projection_coefficient))
+  {
+    for (auto [i_of_primative_mode, primative_mode]
+      : ranges::views::enumerate(primative_projection_coefficient))
+      output.Coefficient[i_of_super_mode][i_of_primative_mode] =
+        std::norm((super_mode.transpose().conjugate() * primative_mode).trace());
+    auto sum = ranges::accumulate(output.Coefficient[i_of_super_mode], 0.);
+    for (auto& coefficient : output.Coefficient[i_of_super_mode])
+      coefficient *=
+        input.Super.Qpoint[config.SuperQpointIndex].Mode[i_of_super_mode].WeightOnUnfold[config.SubQpointIndex] / sum;
+  }
+  // 输出
+  std::ofstream(config.OutputFile) << YAML::Node(output);
+}
--- a/src/raman.cpp
+++ b/src/raman.cpp
@@ -0,0 +1,197 @@
+# include <ufo.hpp>
+
+namespace ufo
+{
+  struct RamanData
+  {
+    std::map<std::pair<std::size_t, std::size_t>, double> ModeRatio;
+    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;
+      // 要计算哪个 q 点的哪些模式，总是假定是 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)
+    {
+      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 move_eigenvector = cell.Qpoint[i_of_qpoint].Mode[i_of_mode].EigenVector.real()
+            .cwiseProduct(mass.cwiseSqrt().cwiseInverse().rowwise().replicate(3)).eval();
+          // 归一化
+          auto ratio = config.MaxDisplacement / move_eigenvector.rowwise().norm().maxCoeff();
+          auto atom_movement = (move_eigenvector * ratio).eval();
+          // 输出
+          for (auto direction : {1, -1})
+          {
+            auto path = "{}/{}/{}/{}"_f
+              (config.OutputPoscarDirectory, i_of_qpoint, i_of_mode, direction > 0 ? "+" : "-");
+            std::filesystem::create_directories(path);
+            std::ofstream("{}/POSCAR"_f(path)) << generate_poscar
+              (cell.Cell, cell.AtomPosition + atom_movement * cell.Cell.inverse() * direction, cell.AtomType);
+          }
+          output.ModeRatio[{i_of_qpoint, i_of_mode}] = ratio;
+          log.debug("Write mode {} {} {}"_f(i_of_qpoint, i_of_mode, atom_movement.rowwise().norm().eval()));
+        }
+    };
+    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.OutputPoscarDirectory)) << 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::string path)
+  {
+    biu::Logger::Guard log(path);
+    auto get_dielectric_tensor = [](std::filesystem::path file)
+    {
+      biu::Logger::Guard log(file);
+      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);
+          Eigen::Matrix3d dielectric_tensor;
+          for (std::size_t i = 0; i < 3; i++) for (std::size_t j = 0; j < 3; j++)
+            in_stream >> dielectric_tensor(i, j);
+          log.debug("read tensor {}"_f(dielectric_tensor));
+          return dielectric_tensor;
+        }
+      }
+      // test if the file is read correctly
+      throw std::runtime_error("Error reading file: {}"_f(file));
+    };
+    auto search_path = [&](this auto self, std::size_t indent, std::string path) -> void
+    {
+      if (std::filesystem::exists(path + "/OUTCAR"))
+      {
+        auto e = get_dielectric_tensor(path + "/OUTCAR");
+        for (std::size_t i = 0; i < 3; i++) std::cout << "{}- [ {}, {}, {} ]\n"_f
+          (std::string(indent * 2, ' '), e(i, 0), e(i, 1), e(i, 2));
+      }
+      else for (auto entry : std::filesystem::directory_iterator(path)) if (entry.is_directory())
+      {
+        std::cout << "{}{}:\n"_f
+        (
+          std::string(indent * 2, ' '),
+          (std::set{"+"s, "-"s}.contains(entry.path().filename()))
+            ? "\"{}\""_f(entry.path().filename()) : entry.path().filename().string()
+        );
+        self(indent + 1, entry.path());
+      }
+    };
+    search_path(0, path);
+  }
+
+  void raman_apply_contribution(std::string config_file)
+  {
+    struct Config
+    {
+      Eigen::Matrix3d OriginalSusceptibility;
+      std::map<std::size_t, std::map<std::size_t, std::map<std::string, Eigen::Matrix3d>>> Susceptibility;
+      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, ratio] : raman_input.ModeRatio)
+      {
+        auto&& [i_of_qpoint, i_of_mode] = i;
+        auto&& _ = cell.Qpoint[i_of_qpoint].Mode[i_of_mode];
+        auto&& susceptibility = config.Susceptibility[i_of_qpoint][i_of_mode];
+        Eigen::Matrix3d raman_tensor = (susceptibility["+"] - susceptibility["-"]) / 2 / ratio;
+        _.RamanTensor = raman_tensor | biu::fromEigen;
+        _.WeightOnRaman = config.Polarization[0].transpose() * raman_tensor * config.Polarization[1];
+        log.info("{}:{:.2f}:{}:{:.2f} {}"_f
+        (
+          i_of_qpoint, fmt::join(cell.Qpoint[i_of_qpoint].Qpoint, ", "),
+          i_of_mode, _.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);
+  }
+}
--- a/src/solver.cpp
+++ b/src/solver.cpp
@@ -1,47 +0,0 @@
-# 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;
-    }
-  };
-
-}
--- a/src/unfold.cpp
+++ b/src/unfold.cpp
@@ -1,429 +1,311 @@
-# include <ufo/unfold.hpp>
+# include <ufo.hpp>
+# include <thread>
+# include <syncstream>
+# include <execution>
+# include <ranges>

-namespace ufo
+void ufo::unfold(std::string config_file)
 {
-  UnfoldSolver::InputType::InputType(std::string filename)
+  // 反折叠的原理: 将超胞中的原子运动状态, 投影到一组平面波构成的基矢中.
+  // 每一个平面波的波矢由两部分相加得到: 一部分是单胞倒格子的整数倍, 所取的个数有一定任意性, 论文中建议取大约单胞中原子个数那么多个，
+  // 但我觉得可以多取一些；这些平面波以原胞为周期，它用于尽可能描述
+  // 另一部分是超胞倒格子的整数倍, 取 n 个, n 为超胞对应的单胞的倍数, 其实也就是倒空间中单胞对应倒格子中超胞的格点.
+  // 只要第一部分取得足够多, 那么单胞中原子的状态就可以完全被这些平面波描述.
+  // 将超胞中原子的运动状态投影到这些基矢上, 计算出投影的系数, 就可以将超胞的原子运动状态分解到单胞中的多个 q 点上.
+
+  struct Config
  {
-    // read main input file
+    // 单胞到超胞的格矢转换时用到的矩阵
+    // SuperCellMultiplier 是一个三维列向量且各个元素都是整数，表示单胞在各个方向扩大到多少倍之后，可以得到和超胞一样的体积
+    // SuperCellDeformation 是一个行列式为 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::Vector3i SuperCellMultiplier;
+    Eigen::Matrix3d SuperCellDeformation;
+    std::optional<std::array<double, 3>> AtomTranslation;
+
+    // 在单胞内取几个平面波的基矢
+    Eigen::Vector3i PrimativeCellBasisNumber;
+
+    // 单胞的 phonopy 输出的 phonopy.yaml，用来读入单胞的晶格、原子坐标、原子类型、原子质量
+    std::string PrimativePhonopy;
+    // 单胞的 phonopy 输出的 band.hdf5 或 qpoint.hdf5，用来读入 q 点和振动模式
+    std::string PrimativeQpoint;
+    // 同上，但是是超胞里的结果
+    std::string SuperPhonopy;
+    std::string SuperQpoint;
+
+    std::string OutputFile;
+  };
+
+  // 从文件中读取 q 点数据
+  // data 为 CommonData::Primative 或 CommonData::Super
+  // 返回值为原子类型和原子质量的对应关系
+  auto read_qpoint = [](std::string phonopy_file, std::string qpoint_file, auto& data)
+  {
+    biu::Logger::Guard log(phonopy_file, qpoint_file, data);
+
+    // phonopy 的输出有两种可能。
+    // 直接指定计算的 q 点时，frequency 是 2 维，这时第一个维度是 q 点，第二个维度是不同模式。
+    // 计算能带时，frequency 是 3 维，相比于二维的情况多了第一个维度，表示 q 点所在路径。
+    // qpoint 或 path，以及 eigenvector 也有类似的变化。
+    // 除此以外，eigenvector 后两个维度分别指示模式的特征向量的各个维度和各个模式（而不是各个模式和特征向量的各个维度），
+    // 因为后两个维度的尺寸总是一样的（模式个数等于原子坐标个数），非常容易搞错。
+    std::vector<std::array<double, 3>> qpoint;
+    std::vector<std::vector<double>> frequency;
+    std::vector<std::vector<std::vector<std::complex<double>>>> eigenvector_vector;
+    auto file = biu::Hdf5file(qpoint_file);
+
+    if (file.File.getDataSet("/frequency").getDimensions().size() == 2)
+      file.read("/frequency", frequency)
+        .read("/eigenvector", eigenvector_vector)
+        .read("/qpoint", qpoint);
+    else
    {
-      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>();
+      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<std::complex<double>>>>> temp_eigenvector_vector;
+      file.read("/frequency", temp_frequency)
+        .read("/eigenvector", temp_eigenvector_vector)
+        .read("/path", temp_path);
+      frequency = temp_frequency | ranges::views::join | ranges::to_vector;
+      qpoint = temp_path | ranges::views::join | ranges::to_vector;
+      eigenvector_vector = temp_eigenvector_vector | ranges::views::join | ranges::to_vector;
+    }

-      for (unsigned i = 0; i < 3; i++)
-        SuperCellMultiplier(i) = node["SuperCellMultiplier"][i].as<int>();
-
-      if (auto value = node["SuperCellDeformation"])
+    // 整理并写入得到结果
+    auto number_of_qpoints = frequency.size(), number_of_modes = frequency[0].size();
+    data.Qpoint.resize(number_of_qpoints);
+    for (auto i : std::views::iota(0u, number_of_qpoints)) 
+    {
+      data.Qpoint[i].Qpoint = qpoint[i] | biu::toEigen<>;
+      data.Qpoint[i].Mode.resize(number_of_modes);
+      for (auto j : std::views::iota(0u, number_of_modes)) 
      {
-        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
-          );
+        data.Qpoint[i].Mode[j].Frequency = frequency[i][j];
+        auto number_of_atoms = number_of_modes / 3;
+        Eigen::MatrixX3cd eigenvector(number_of_atoms, 3);
+        for (auto k : std::views::iota(0u, number_of_atoms))
+          for (auto l : std::views::iota(0u, 3u))
+            eigenvector(k, l) = eigenvector_vector[i][k * 3 + l][j];
+        // 原则上讲，需要对读入的原子运动状态作相位转换, 使得它们与我们的约定一致(对超胞周期性重复)，但这个转换 phonopy 已经做了
+        // 这里还要需要做归一化处理 (指将数据简单地作为向量处理的归一化)，但似乎 phonopy 也已经做了，但我们也再做一次，也没什么坏处
+        data.Qpoint[i].Mode[j].EigenVector = eigenvector / eigenvector.norm();
      }
    }

-    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();
-          }
-        }
-    }
-  }
+    // 读取并写入其它数据
+    YAML::Node phonopy = YAML::LoadFile(phonopy_file);
+    data.Cell = phonopy["unit_cell"]["lattice"].as<Eigen::Matrix3d>();
+    auto points = phonopy["unit_cell"]["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>>;
+  };

-  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
+  // 构建基
+  // 每个 q 点对应一组 sub qpoint。不同的 q 点所对应的 sub qpoint 是不一样的，但 sub qpoint 与 q 点的相对位移在不同 q 点之间是相同的。
+  // 由于基只与这个相对位置有关（也就是说，不同 q 点的基是一样的），因此可以先计算出所有的基，这样降低计算量。
+  // 外层下标对应超胞倒格子的整数倍那部分(第二部分), 也就是不同的 sub qpoint
+  // 内层下标对应单胞倒格子的整数倍那部分(第一部分), 也就是 sub qpoint 上的不同平面波（取的数量越多，结果越精确）
+  auto 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
+    Eigen::Matrix3d primative_cell, Eigen::Vector3i super_cell_multiplier,
+    Eigen::Vector3i primative_cell_basis_number, Eigen::MatrixX3d 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))
+    biu::Logger::Guard log;
+    std::vector<std::vector<Eigen::VectorXcd>> basis(super_cell_multiplier.prod());
+    // diff_of_sub_qpoint 表示 sub qpoint 与 qpoint 的相对位置，单位为超胞的倒格矢
+    for (auto [diff_of_sub_qpoint_by_reciprocal_modified_super_cell, i_of_sub_qpoint]
+      : biu::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))
+      basis[i_of_sub_qpoint].resize((primative_cell_basis_number.array() * 2 - 1).prod());
+      for (auto [xyz_of_basis, i_of_basis] : biu::sequence
+        ((-primative_cell_basis_number + Eigen::Vector3i::Ones()).eval(), primative_cell_basis_number))
      {
-        // 计算 q 点的坐标, 单位为单胞的倒格矢
-        auto diff_of_sub_qpoint_by_reciprocal_primative_cell = xyz_of_basis.cast<double>()
+        // 计算波矢, 单位为单胞的倒格矢
+        auto wavevector_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();
+          * diff_of_sub_qpoint_by_reciprocal_modified_super_cell.cast<double>();
+        // 将单位转换为埃^-1
+        auto wavevector =
+          (wavevector_by_reciprocal_primative_cell.transpose() * primative_cell.inverse().transpose()).transpose();
        // 计算基矢
        basis[i_of_sub_qpoint][i_of_basis]
-          = (2i * std::numbers::pi_v<double> * (atom_position * qpoint)).array().exp();
+          = (2i * std::numbers::pi_v<double> * (atom_position * wavevector)).array().exp();
      }
    }
    return basis;
-  }
+  };

-  std::vector<std::vector<double>> UnfoldSolver::construct_projection_coefficient
+  // 计算从超胞到原胞的投影系数（不是分原子的投影系数），是反折叠的核心步骤
+  // 返回的投影系数是一个三维数组，第一维对应不同的 q 点，第二维对应不同的模式，第三维对应不同的 sub qpoint
+  auto 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
+    const std::vector<std::vector<Eigen::VectorXcd>>& basis,
+    const std::vector<std::reference_wrapper<const Eigen::MatrixX3cd>>& modes,
+    std::atomic<std::size_t>& number_of_finished_modes
  )
  {
    // 第一层下标对应不同模式, 第二层下标对应这个模式在反折叠后的 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, 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 (unsigned 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 (unsigned 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())
-                .array().abs2().sum();
+              (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;
+        auto sum = ranges::accumulate(projection_coefficient, 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;
+  biu::Logger::Guard log;
+
+  log.info("Reading input file...");
+  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;
+  output.Super.AtomTranslation = config.AtomTranslation;
+  // 填充 SubQpoint
+  for (auto i_of_super_qpoint : std::views::iota(0u, output.Super.Qpoint.size()))
    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++
+      auto [diff_of_sub_qpoint_by_reciprocal_modified_super_cell, i_of_sub_qpoint]
+        : biu::sequence(config.SuperCellMultiplier)
    )
    {
-      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();
+      /*
+        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 =
+      (
+        config.SuperCellMultiplier.cast<double>().cwiseInverse().asDiagonal()
+        * (
+          diff_of_sub_qpoint_by_reciprocal_modified_super_cell.cast<double>()
+            + config.SuperCellDeformation.inverse() * output.Super.Qpoint[i_of_super_qpoint].Qpoint
+        )
+      ).eval();
+      output.Super.Qpoint[i_of_super_qpoint].SubQpoint.push_back(sub_qpoint_by_reciprocal_primative_cell.array()
+        - sub_qpoint_by_reciprocal_primative_cell.array().floor());
    }
-    return output;
+  log.info("Done.");
+
+  log.info("Constructing basis...");
+  auto basis = construct_basis
+  (
+    output.Primative.Cell, output.Super.CellMultiplier,
+    config.PrimativeCellBasisNumber,
+    output.Super.AtomPosition
+      * (output.Super.CellDeformation * output.Super.CellMultiplier.cast<double>().asDiagonal() * output.Primative.Cell)
+  );
+  log.info("Done.");
+
+  std::clog << "Calculating projection coefficient... " << std::flush;
+  // 将所有q点的模式flatten一下，放到一个一维数组中，方便并行计算
+  {
+    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([&]
+    {
+      while (true)
+      {
+        std::osyncstream(std::clog)
+          << "\rCalculating projection coefficient... ({}/{})"_f(number_of_finished_modes, modes.size())
+          << std::flush;
+        std::this_thread::sleep_for(100ms);
+        if (number_of_finished_modes == modes.size()) break;
+      }
+    });
+    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;
+
+  log.info("Writing data... ");
+  std::ofstream(config.OutputFile, std::ios::binary) << biu::serialize<char>(output);
+  log.info("Done.");
+
+  log.info("Summary:");
+
+  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()))
+      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]
+        ));
 }
--- a/test/14.2.2.4.unfold.yaml
+++ b/test/14.2.2.4.unfold.yaml
@@ -1,12 +0,0 @@
-PrimativeCell:
- [     1.548010265167714,    -2.681232429908652,     0.000000000000000 ] # a
- [     1.548010265167714,     2.681232429908652,     0.000000000000000 ] # b
- [     0.000000000000000,     0.000000000000000,     5.061224103556595 ] # c
-SuperCellMultiplier: [ 3, 1, 1 ]
-PrimativeCellBasisNumber: [ 8, 8, 8 ]
-AtomPositionInputFile: { FileName: "/home/chn/Documents/lammps-SiC/14/14.2/14.2.2/14.2.2.4/band.yaml", Format: "yaml" }
-QpointDataInputFile: { FileName: "/home/chn/Documents/lammps-SiC/14/14.2/14.2.2/14.2.2.4/band.yaml", Format: "yaml" }
-QpointDataOutputFile:
- { FileName: "test/14.2.2.4.result.yaml", Format: "yaml" }
- { FileName: "test/14.2.2.4.result.human-readable.yaml", Format: "yaml-human-readable" }
- { FileName: "test/14.2.2.4.result.zpp", Format: "zpp" }
--- a/test/14.2.2.5.plot.yaml
+++ b/test/14.2.2.5.plot.yaml
@@ -1,14 +0,0 @@
-PrimativeCell:
- [     1.548010265167714,    -2.681232429908652,     0.000000000000000 ] # a
- [     1.548010265167714,     2.681232429908652,     0.000000000000000 ] # b
- [     0.000000000000000,     0.000000000000000,     5.061224103556595 ] # c
-SourceFilename: ./test/14.2.2.5.result.zpp
-Figures:
- Qpoints:
-  - - [ 0, 0, 0 ]
-    - [ 0.5, 0, 0 ]
-    - [ 0.3333333, 0.3333333, 0 ]
-    - [ 0, 0, 0 ]
-  Resolution: [ 1024, 400 ]
-  Range: [ -5, 35 ]
-  Filename: ./test/14.2.2.5.png
--- a/test/14.2.2.5.unfold.yaml
+++ b/test/14.2.2.5.unfold.yaml
@@ -1,12 +0,0 @@
-PrimativeCell:
- [     1.548010265167714,    -2.681232429908652,     0.000000000000000 ] # a
- [     1.548010265167714,     2.681232429908652,     0.000000000000000 ] # b
- [     0.000000000000000,     0.000000000000000,     5.061224103556595 ] # c
-SuperCellMultiplier: [ 3, 1, 1 ]
-PrimativeCellBasisNumber: [ 8, 8, 8 ]
-AtomPositionInputFile: { FileName: "/home/chn/Documents/lammps-SiC/14/14.2/14.2.2/14.2.2.5/band.yaml", Format: "yaml" }
-QpointDataInputFile: { FileName: "/home/chn/Documents/lammps-SiC/14/14.2/14.2.2/14.2.2.5/band.yaml", Format: "yaml" }
-QpointDataOutputFile:
- { FileName: "test/14.2.2.5.result.yaml", Format: "yaml" }
- { FileName: "test/14.2.2.5.result.human-readable.yaml", Format: "yaml-human-readable" }
- { FileName: "test/14.2.2.5.result.zpp", Format: "zpp" }
--- a/test/14.2.4.4.plot.yaml
+++ b/test/14.2.4.4.plot.yaml
@@ -1,14 +0,0 @@
-PrimativeCell:
- [     1.548010265167714,    -2.681232429908652,     0.000000000000000 ] # a
- [     1.548010265167714,     2.681232429908652,     0.000000000000000 ] # b
- [     0.000000000000000,     0.000000000000000,     5.061224103556595 ] # c
-SourceFilename: ./test/14.2.4.4.result.zpp
-Figures:
- Qpoints:
-  - - [ 0, 0, 0 ]
-    - [ 0.5, 0, 0 ]
-    - [ 0.3333333, 0.3333333, 0 ]
-    - [ 0, 0, 0 ]
-  Resolution: [ 1024, 400 ]
-  Range: [ -5, 35 ]
-  Filename: ./test/14.2.4.4.png
--- a/test/14.2.4.4.unfold.yaml
+++ b/test/14.2.4.4.unfold.yaml
@@ -1,14 +0,0 @@
-PrimativeCell:
- [     1.548010265167714,    -2.681232429908652,     0.000000000000000 ] # a
- [     1.548010265167714,     2.681232429908652,     0.000000000000000 ] # b
- [     0.000000000000000,     0.000000000000000,     5.061224103556595 ] # c
-SuperCellMultiplier: [ 2, 2, 1 ]
-PrimativeCellBasisNumber: [ 8, 8, 8 ]
-AtomPositionInputFile: 
-  FileName: "/home/chn/Documents/lammps-SiC/14/14.2/14.2.4/14.2.4.4/phonopy.yaml"
-  Format: "yaml"
-QpointDataInputFile: { FileName: "/home/chn/Documents/lammps-SiC/14/14.2/14.2.4/14.2.4.4/band.hdf5", Format: "hdf5" }
-QpointDataOutputFile:
- { FileName: "test/14.2.4.4.result.yaml", Format: "yaml" }
- { FileName: "test/14.2.4.4.result.human-readable.yaml", Format: "yaml-human-readable" }
- { FileName: "test/14.2.4.4.result.zpp", Format: "zpp" }
--- a/test/14.2.5.5.plot.yaml
+++ b/test/14.2.5.5.plot.yaml
@@ -1,14 +0,0 @@
-PrimativeCell:
- [     1.548010265167714,    -2.681232429908652,     0.000000000000000 ] # a
- [     1.548010265167714,     2.681232429908652,     0.000000000000000 ] # b
- [     0.000000000000000,     0.000000000000000,     5.061224103556595 ] # c
-SourceFilename: ./test/14.2.5.5.result.zpp
-Figures:
- Qpoints:
-  - - [ 0, 0, 0 ]
-    - [ 0.5, 0, 0 ]
-    - [ 0.3333333, 0.3333333, 0 ]
-    - [ 0, 0, 0 ]
-  Resolution: [ 1024, 400 ]
-  Range: [ -5, 35 ]
-  Filename: ./test/14.2.5.5.png
--- a/test/14.2.5.5.unfold.yaml
+++ b/test/14.2.5.5.unfold.yaml
@@ -1,18 +0,0 @@
-PrimativeCell:
- [     1.548010265167714,    -2.681232429908652,     0.000000000000000 ] # a
- [     1.548010265167714,     2.681232429908652,     0.000000000000000 ] # b
- [     0.000000000000000,     0.000000000000000,     5.061224103556595 ] # c
-SuperCellMultiplier: [ 1, 2, 1 ]
-SuperCellDeformation:
- [ 1, 0, 0]
- [ 1, 1, 0]
- [ 0, 0, 1]
-PrimativeCellBasisNumber: [ 8, 8, 8 ]
-AtomPositionInputFile:
-  FileName: "/home/chn/Documents/lammps-SiC/14/14.2/14.2.5/14.2.5.5/phonopy.yaml"
-  Format: "yaml"
-QpointDataInputFile: { FileName: "/home/chn/Documents/lammps-SiC/14/14.2/14.2.5/14.2.5.5/band.hdf5", Format: "hdf5" }
-QpointDataOutputFile:
- { FileName: "test/14.2.5.5.result.yaml", Format: "yaml" }
- { FileName: "test/14.2.5.5.result.human-readable.yaml", Format: "yaml-human-readable" }
- { FileName: "test/14.2.5.5.result.zpp", Format: "zpp" }
--- a/test/14.2.6.4.plot.yaml
+++ b/test/14.2.6.4.plot.yaml
@@ -1,15 +0,0 @@
-PrimativeCell:
- [     1.548010265167714,    -2.681232429908652,     0.000000000000000 ] # a
- [     1.548010265167714,     2.681232429908652,     0.000000000000000 ] # b
- [     0.000000000000000,     0.000000000000000,     5.061224103556595 ] # c
-SourceFilename: ./test/14.2.6.4.result.zpp
-Figures:
- Qpoints:
-  - - [ 0, 0, 0 ]
-    - [ 0.5, 0, 0 ]
-    - [ 0.3333333, 0.3333333, 0 ]
-    - [ 0, 0, 0 ]
-  Resolution: [ 1024, 1024 ]
-  Range: [ -5, 35 ]
-  Filename: ./test/14.2.6.4.png
-  YTicks: [ 0, 5, 10, 15, 20, 25, 30 ]
--- a/test/14.2.6.4.unfold.yaml
+++ b/test/14.2.6.4.unfold.yaml
@@ -1,18 +0,0 @@
-PrimativeCell:
- [     1.548010265167714,    -2.681232429908652,     0.000000000000000 ] # a
- [     1.548010265167714,     2.681232429908652,     0.000000000000000 ] # b
- [     0.000000000000000,     0.000000000000000,     5.061224103556595 ] # c
-SuperCellMultiplier: [ 4, 8, 1 ]
-SuperCellDeformation:
- [ 1, 0, 0]
- [ 1, 1, 0]
- [ 0, 0, 1]
-PrimativeCellBasisNumber: [ 8, 8, 8 ]
-AtomPositionInputFile:
-  FileName: "/home/chn/Documents/lammps-SiC/14/14.2/14.2.6/14.2.6.4/phonopy.yaml"
-  Format: "yaml"
-QpointDataInputFile: { FileName: "/home/chn/Documents/lammps-SiC/14/14.2/14.2.6/14.2.6.4/band.hdf5", Format: "hdf5" }
-QpointDataOutputFile:
- { FileName: "test/14.2.6.4.result.yaml", Format: "yaml" }
- { FileName: "test/14.2.6.4.result.human-readable.yaml", Format: "yaml-human-readable" }
- { FileName: "test/14.2.6.4.result.zpp", Format: "zpp" }
--- a/example/step2/fold.yaml
+++ b/example/step2/fold.yaml
@@ -1,7 +1,7 @@
-SuperCellMultiplier: [ 3, 4, 1 ]
+SuperCellMultiplier: [3, 4, 1]
 SuperCellDeformation:
  - [ 1, 0, 0 ]
-  - [ 0.6666666, 1, 0 ]
+  - [ 0.6666, 1, 0 ]
  - [ 0, 0, 1 ]
 Qpoints:
  - [0, 0, 0]
@@ -15,6 +15,4 @@ Qpoints:
  - [0.4, 0, 0]
  - [0.45, 0, 0]
  - [0.5, 0, 0]
-OutputFile:
-  Filename: fold-output.yaml
-  Format: yaml
+OutputFile: fold-output.yaml
--- a/test/project-to-atom/.gitignore
+++ b/test/project-to-atom/.gitignore
@@ -0,0 +1 @@
+output.zpp
--- a/test/project-to-atom/config.yaml
+++ b/test/project-to-atom/config.yaml
@@ -0,0 +1,3 @@
+SelectedAtom: [0, 2, 4]
+InputFile: input.zpp
+OutputFile: output.zpp
--- a/test/project-to-atom/input.zpp
+++ b/test/project-to-atom/input.zpp
--- a/test/project-to-mode/.gitignore
+++ b/test/project-to-mode/.gitignore
@@ -0,0 +1 @@
+output.yaml
--- a/test/project-to-mode/config.yaml
+++ b/test/project-to-mode/config.yaml
@@ -0,0 +1,6 @@
+PrimativeCellBasisNumber: [ 8, 8, 8 ]
+SuperQpointIndex: 0
+SubQpointIndex: 0
+PrimativeQpointIndex: 0
+InputDataFile: data.zpp
+OutputFile: output.yaml
--- a/test/project-to-mode/data.zpp
+++ b/test/project-to-mode/data.zpp
--- a/test/raman-apply-contribution/.gitignore
+++ b/test/raman-apply-contribution/.gitignore
@@ -0,0 +1 @@
+output.zpp
--- a/test/raman-apply-contribution/config.yaml
+++ b/test/raman-apply-contribution/config.yaml
@@ -0,0 +1,228 @@
+OriginalSusceptibility:
+  - [ 7.006328, -3.3e-05, -2e-06 ]
+  - [ -3.4e-05, 7.002539, -2.4e-05 ]
+  - [ -3e-06, -1.2e-05, 7.303458 ]
+Susceptibility:
+  0:
+    0:
+      "+":
+        - [ 7.006379, -5e-06, -2e-06 ]
+        - [ -1.6e-05, 7.002547, -3e-05 ]
+        - [ -1.1e-05, -6e-06, 7.303627 ]
+      "-":
+        - [ 7.006463, 4e-06, -4.5e-05 ]
+        - [ 1e-06, 7.002729, -1.7e-05 ]
+        - [ -4.9e-05, 6e-06, 7.303617 ]
+    1:
+      "+":
+        - [ 7.006331, 7e-06, 1.6e-05 ]
+        - [ 2e-06, 7.002539, -1.9e-05 ]
+        - [ 1.7e-05, -3e-06, 7.303612 ]
+      "-":
+        - [ 7.006366, -8e-06, 1.5e-05 ]
+        - [ -1.1e-05, 7.002548, -8e-06 ]
+        - [ 1.3e-05, 1.1e-05, 7.303606 ]
+    10:
+      "+":
+        - [ 7.00768, 0.004054, 1e-06 ]
+        - [ 0.004046, 7.003291, -4e-05 ]
+        - [ -1e-06, -1.7e-05, 7.305554 ]
+      "-":
+        - [ 7.007955, -0.004016, -1e-06 ]
+        - [ -0.004015, 7.003451, -9e-06 ]
+        - [ -4e-06, 7e-06, 7.305848 ]
+    107:
+      "+":
+        - [ 7.007532, -4.7e-05, -6e-06 ]
+        - [ -4.8e-05, 7.003828, -9e-06 ]
+        - [ -4e-06, 4e-06, 7.307759 ]
+      "-":
+        - [ 7.0076, -5e-06, -1.2e-05 ]
+        - [ -5e-06, 7.003875, -8e-06 ]
+        - [ -1e-05, 2e-06, 7.307945 ]
+    115:
+      "+":
+        - [ 7.007683, 5e-06, -5.3e-05 ]
+        - [ -0, 7.003913, 5e-06 ]
+        - [ -4e-05, 2.8e-05, 7.307388 ]
+      "-":
+        - [ 7.007917, 5e-06, -3e-06 ]
+        - [ 4e-06, 7.004128, -3e-06 ]
+        - [ -2e-06, 1.1e-05, 7.307671 ]
+    15:
+      "+":
+        - [ 7.010948, -8.3e-05, -0.007718 ]
+        - [ -8e-05, 7.004814, -3.4e-05 ]
+        - [ -0.00772, -1.5e-05, 7.3101 ]
+      "-":
+        - [ 7.011174, 3e-06, 0.00773 ]
+        - [ -5e-06, 7.005176, -7e-06 ]
+        - [ 0.007728, 1e-05, 7.3104 ]
+    16:
+      "+":
+        - [ 7.008787, 7.2e-05, -1.8e-05 ]
+        - [ 7.4e-05, 7.007313, 0.007694 ]
+        - [ -2.4e-05, 0.007696, 7.310191 ]
+      "-":
+        - [ 7.008947, -1.3e-05, -7e-06 ]
+        - [ -2.1e-05, 7.007438, -0.007742 ]
+        - [ -7e-06, -0.007727, 7.310434 ]
+    185:
+      "+":
+        - [ 7.018886, 1.8e-05, -5e-05 ]
+        - [ 2.7e-05, 7.015141, -3e-05 ]
+        - [ -5e-05, -1.7e-05, 7.307471 ]
+      "-":
+        - [ 7.003657, 1e-06, 9e-06 ]
+        - [ -1e-06, 6.999828, -1.4e-05 ]
+        - [ 1e-05, 5e-06, 7.326057 ]
+    2:
+      "+":
+        - [ 7.006421, -8e-06, -4e-05 ]
+        - [ -1.5e-05, 7.002682, -6.7e-05 ]
+        - [ -3.9e-05, -4.5e-05, 7.303477 ]
+      "-":
+        - [ 7.00633, 2.4e-05, -7.4e-05 ]
+        - [ 2.2e-05, 7.00253, -9.1e-05 ]
+        - [ -7.7e-05, -7.2e-05, 7.303417 ]
+    251:
+      "+":
+        - [ 7.008862, 4.9e-05, -8.8e-05 ]
+        - [ 4.8e-05, 7.011172, -0.00086 ]
+        - [ -7.8e-05, -0.000845, 7.305162 ]
+      "-":
+        - [ 7.009107, -5e-06, -2e-06 ]
+        - [ -8e-06, 7.011587, 0.000791 ]
+        - [ 3e-06, 0.000809, 7.305493 ]
+    252:
+      "+":
+        - [ 7.015066, -4.3e-05, -0.000838 ]
+        - [ -4.5e-05, 7.005147, -2.9e-05 ]
+        - [ -0.000841, -9e-06, 7.305197 ]
+      "-":
+        - [ 7.015285, -4e-06, 0.000833 ]
+        - [ -2e-06, 7.005368, -4e-06 ]
+        - [ 0.000824, 1.2e-05, 7.305474 ]
+    281:
+      "+":
+        - [ 6.981143, 1.4e-05, -5.8e-05 ]
+        - [ 1.1e-05, 7.038724, 4.3e-05 ]
+        - [ -5.7e-05, 5.9e-05, 7.306158 ]
+      "-":
+        - [ 7.037206, 2e-06, 0 ]
+        - [ 3e-06, 6.983191, -1e-05 ]
+        - [ 2e-06, 0, 7.306576 ]
+    282:
+      "+":
+        - [ 7.014582, -0.027613, 1.6e-05 ]
+        - [ -0.027632, 7.005279, -7.1e-05 ]
+        - [ 4e-06, -6.7e-05, 7.306158 ]
+      "-":
+        - [ 7.014711, 0.027632, -2e-06 ]
+        - [ 0.027633, 7.005478, -1.1e-05 ]
+        - [ -3e-06, 4e-06, 7.306578 ]
+    287:
+      "+":
+        - [ 7.032791, 2.6e-05, -4.1e-05 ]
+        - [ 3.2e-05, 7.028882, -1.7e-05 ]
+        - [ -3.6e-05, -2e-06, 7.250185 ]
+      "-":
+        - [ 6.985897, -3e-06, 4e-06 ]
+        - [ 5e-06, 6.982421, -1.4e-05 ]
+        - [ 2e-06, -3e-06, 7.382705 ]
+    292:
+      "+":
+        - [ 7.007262, 1.3e-05, 1.3e-05 ]
+        - [ 9e-06, 7.013971, -5.9e-05 ]
+        - [ 7e-06, -4.6e-05, 7.305084 ]
+      "-":
+        - [ 7.011688, 0, 2e-06 ]
+        - [ -2e-06, 7.009318, -7e-06 ]
+        - [ 4e-06, 9e-06, 7.306455 ]
+    293:
+      "+":
+        - [ 7.015281, 0.002384, -0 ]
+        - [ 0.002373, 7.005685, -4.4e-05 ]
+        - [ -4e-06, -2.8e-05, 7.3058 ]
+      "-":
+        - [ 7.015423, -0.002375, -9e-06 ]
+        - [ -0.00237, 7.005799, -7e-06 ]
+        - [ -8e-06, 1.1e-05, 7.305927 ]
+    311:
+      "+":
+        - [ 7.017108, 0.000122, -0.023832 ]
+        - [ 0.00013, 7.006382, -5.6e-05 ]
+        - [ -0.02384, -3.4e-05, 7.306937 ]
+      "-":
+        - [ 7.017506, 5e-06, 0.023741 ]
+        - [ -6e-06, 7.006811, -2e-06 ]
+        - [ 0.023733, 1.6e-05, 7.307537 ]
+    312:
+      "+":
+        - [ 7.010335, 7.8e-05, -4.6e-05 ]
+        - [ 8.4e-05, 7.013598, 0.023664 ]
+        - [ -4.1e-05, 0.023683, 7.307236 ]
+      "-":
+        - [ 7.010567, -1e-06, -3e-06 ]
+        - [ -9e-06, 7.013715, -0.023726 ]
+        - [ -4e-06, -0.023709, 7.307591 ]
+    327:
+      "+":
+        - [ 7.008517, 2.9e-05, -1.6e-05 ]
+        - [ 2.5e-05, 7.004613, -9.3e-05 ]
+        - [ -1.8e-05, -7.1e-05, 7.324773 ]
+      "-":
+        - [ 7.012664, -1e-05, 1e-06 ]
+        - [ -8e-06, 7.008951, -1.2e-05 ]
+        - [ 5e-06, 7e-06, 7.306262 ]
+    378:
+      "+":
+        - [ 7.012128, 1e-05, 1.8e-05 ]
+        - [ 1e-06, 7.008401, -1e-05 ]
+        - [ 1.3e-05, -1e-06, 7.322152 ]
+      "-":
+        - [ 7.012153, 1e-06, -1.3e-05 ]
+        - [ 9e-06, 7.008396, -1.4e-05 ]
+        - [ -1.9e-05, 3e-06, 7.322317 ]
+    379:
+      "+":
+        - [ 7.010873, -5.9e-05, -3.8e-05 ]
+        - [ -7.8e-05, 7.007166, -3.6e-05 ]
+        - [ -4.3e-05, -2.5e-05, 7.319279 ]
+      "-":
+        - [ 7.011186, -5e-06, -1.4e-05 ]
+        - [ 6e-06, 7.007443, -1e-06 ]
+        - [ -1.9e-05, 1.4e-05, 7.319705 ]
+    5:
+      "+":
+        - [ 7.007563, -0.001516, 5e-06 ]
+        - [ -0.001515, 7.00316, -5.5e-05 ]
+        - [ 0, -2.9e-05, 7.305665 ]
+      "-":
+        - [ 7.007833, 0.001555, -9e-06 ]
+        - [ 0.001554, 7.003468, -1.1e-05 ]
+        - [ -1.2e-05, 1e-05, 7.306111 ]
+    6:
+      "+":
+        - [ 7.005528, 1.8e-05, 4.1e-05 ]
+        - [ 1.2e-05, 7.005511, -3.7e-05 ]
+        - [ 3.7e-05, -1.5e-05, 7.30606 ]
+      "-":
+        - [ 7.008811, -2e-06, 4e-06 ]
+        - [ -5e-06, 7.002567, -3e-06 ]
+        - [ 7e-06, 1.7e-05, 7.306146 ]
+    9:
+      "+":
+        - [ 7.011048, 8e-05, -2.7e-05 ]
+        - [ 8.6e-05, 7.000001, 6e-06 ]
+        - [ -2.6e-05, 4.1e-05, 7.30555 ]
+      "-":
+        - [ 7.003087, -6e-06, 3e-06 ]
+        - [ -3e-06, 7.008199, -1e-06 ]
+        - [ 5e-06, 1.8e-05, 7.305867 ]
+Polarization:
+  - [ 0, 1, 0 ]
+  - [ 0, 1, 0 ]
+InputDataFile: data.zpp
+RamanInputDataFile: raman.zpp
+OutputDataFile: output.zpp
--- a/test/raman-apply-contribution/data.zpp
+++ b/test/raman-apply-contribution/data.zpp
--- a/test/raman-apply-contribution/raman.zpp
+++ b/test/raman-apply-contribution/raman.zpp
--- a/test/raman-create-displacement/.gitignore
+++ b/test/raman-create-displacement/.gitignore
@@ -0,0 +1,2 @@
+POSCAR
+output.zpp
--- a/test/raman-create-displacement/config.yaml
+++ b/test/raman-create-displacement/config.yaml
@@ -0,0 +1,7 @@
+Cell: Super
+SelectedModes: 
+  0: [ 0, 1, 2, 5, 6, 9, 10, 15, 16, 107, 115, 185, 251, 252, 281, 282, 287, 292, 293, 311, 312, 327, 378, 379 ]
+MaxDisplacement: 0.02
+OutputPoscarDirectory: POSCAR
+InputDataFile: data.zpp
+OutputDataFile: output.zpp
--- a/test/raman-create-displacement/data.zpp
+++ b/test/raman-create-displacement/data.zpp
--- a/test/raman-extract/.gitattributes
+++ b/test/raman-extract/.gitattributes
@@ -0,0 +1 @@
+OUTCAR filter=lfs diff=lfs merge=lfs -text
--- a/test/raman-extract/0/0/+/OUTCAR
+++ b/test/raman-extract/0/0/+/OUTCAR
--- a/test/raman-extract/0/0/-/OUTCAR
+++ b/test/raman-extract/0/0/-/OUTCAR
--- a/test/raman-extract/0/1/+/OUTCAR
+++ b/test/raman-extract/0/1/+/OUTCAR
--- a/test/raman-extract/0/1/-/OUTCAR
+++ b/test/raman-extract/0/1/-/OUTCAR
--- a/test/raman-extract/0/10/+/OUTCAR
+++ b/test/raman-extract/0/10/+/OUTCAR
--- a/test/raman-extract/0/10/-/OUTCAR
+++ b/test/raman-extract/0/10/-/OUTCAR
--- a/test/raman-extract/0/107/+/OUTCAR
+++ b/test/raman-extract/0/107/+/OUTCAR
--- a/test/raman-extract/0/107/-/OUTCAR
+++ b/test/raman-extract/0/107/-/OUTCAR
--- a/test/raman-extract/0/115/+/OUTCAR
+++ b/test/raman-extract/0/115/+/OUTCAR
--- a/test/raman-extract/0/115/-/OUTCAR
+++ b/test/raman-extract/0/115/-/OUTCAR
--- a/test/raman-extract/0/15/+/OUTCAR
+++ b/test/raman-extract/0/15/+/OUTCAR
--- a/test/raman-extract/0/15/-/OUTCAR
+++ b/test/raman-extract/0/15/-/OUTCAR
--- a/test/raman-extract/0/16/+/OUTCAR
+++ b/test/raman-extract/0/16/+/OUTCAR
--- a/test/raman-extract/0/16/-/OUTCAR
+++ b/test/raman-extract/0/16/-/OUTCAR
--- a/test/raman-extract/0/185/+/OUTCAR
+++ b/test/raman-extract/0/185/+/OUTCAR
--- a/test/raman-extract/0/185/-/OUTCAR
+++ b/test/raman-extract/0/185/-/OUTCAR
--- a/test/raman-extract/0/2/+/OUTCAR
+++ b/test/raman-extract/0/2/+/OUTCAR
--- a/test/raman-extract/0/2/-/OUTCAR
+++ b/test/raman-extract/0/2/-/OUTCAR
--- a/test/raman-extract/0/251/+/OUTCAR
+++ b/test/raman-extract/0/251/+/OUTCAR
--- a/test/raman-extract/0/251/-/OUTCAR
+++ b/test/raman-extract/0/251/-/OUTCAR
--- a/test/raman-extract/0/252/+/OUTCAR
+++ b/test/raman-extract/0/252/+/OUTCAR
--- a/test/raman-extract/0/252/-/OUTCAR
+++ b/test/raman-extract/0/252/-/OUTCAR
--- a/test/raman-extract/0/281/+/OUTCAR
+++ b/test/raman-extract/0/281/+/OUTCAR
--- a/test/raman-extract/0/281/-/OUTCAR
+++ b/test/raman-extract/0/281/-/OUTCAR
--- a/test/raman-extract/0/282/+/OUTCAR
+++ b/test/raman-extract/0/282/+/OUTCAR
--- a/test/raman-extract/0/282/-/OUTCAR
+++ b/test/raman-extract/0/282/-/OUTCAR
--- a/test/raman-extract/0/287/+/OUTCAR
+++ b/test/raman-extract/0/287/+/OUTCAR
--- a/test/raman-extract/0/287/-/OUTCAR
+++ b/test/raman-extract/0/287/-/OUTCAR
--- a/test/raman-extract/0/292/+/OUTCAR
+++ b/test/raman-extract/0/292/+/OUTCAR
--- a/test/raman-extract/0/292/-/OUTCAR
+++ b/test/raman-extract/0/292/-/OUTCAR
--- a/test/raman-extract/0/293/+/OUTCAR
+++ b/test/raman-extract/0/293/+/OUTCAR
--- a/test/raman-extract/0/293/-/OUTCAR
+++ b/test/raman-extract/0/293/-/OUTCAR
--- a/test/raman-extract/0/311/+/OUTCAR
+++ b/test/raman-extract/0/311/+/OUTCAR
--- a/test/raman-extract/0/311/-/OUTCAR
+++ b/test/raman-extract/0/311/-/OUTCAR
--- a/test/raman-extract/0/312/+/OUTCAR
+++ b/test/raman-extract/0/312/+/OUTCAR
--- a/test/raman-extract/0/312/-/OUTCAR
+++ b/test/raman-extract/0/312/-/OUTCAR
--- a/Show More
+++ b/Show More
Author	SHA1	Message	Date
Haonan Chen	b130375ff8	PrimativeCellBasisNumber 使用 gamma center 的取法	2026-01-11 14:53:12 +08:00
Haonan Chen	82ad565448	AtomTranslation 改为以超胞晶格矢为单位，而不是以埃为单位	2026-01-08 12:36:14 +08:00
Haonan Chen	cb7dd54e7f	一些小修改，应该不影响结果	2026-01-08 12:36:14 +08:00
Haonan Chen	19736049c9	fix build	2026-01-05 15:21:35 +08:00
chn	dc6b431bf9	fix raman tensor scaling	2025-11-05 12:05:12 +08:00
chn	32febb2843	Revert "fix raman tensor calculation scaling" This reverts commit `a54e769b24`.	2025-11-05 12:02:52 +08:00
chn	a54e769b24	fix raman tensor calculation scaling	2025-11-05 11:37:44 +08:00
chn	8a137d5342	a little improvement in code (should not change result)	2025-10-11 14:45:59 +08:00
chn	b21a8fb198	Revert "fix raman" This reverts commit `2ddfca800a`.	2025-10-11 10:18:55 +08:00
chn	2ddfca800a	fix raman	2025-10-10 16:29:16 +08:00
chn	745353d089	project-to-atom: fix	2025-04-01 18:15:43 +08:00
chn	c4c549c669	project-to-atom: add test	2025-04-01 17:30:01 +08:00
chn	388c77db5f	remove unused config	2025-04-01 17:22:26 +08:00
chn	d4578c7991	fix	2024-12-24 17:47:02 +08:00
chn	f871c28768	修改拉曼输出文件布局	2024-12-22 15:53:11 +08:00
chn	d584951fd7	fix raman, add test	2024-12-04 15:21:09 +08:00
chn	6e64054ffc	fix	2024-12-03 19:11:48 +08:00
chn	0c073b73b5	clean up	2024-12-03 18:15:28 +08:00
chn	bf3c8430ed	解决bug	2024-12-03 17:56:22 +08:00
chn	eb6b903417	写了 project_to_mode 但还没调试好	2024-12-03 16:25:17 +08:00
chn	3f98ce8378	add zpp-to-yaml test	2024-12-02 18:51:17 +08:00
chn	7b4ea3bbf3	add unfold test	2024-12-02 18:27:40 +08:00
chn	db1d3befd6	重构	2024-12-02 15:47:33 +08:00
chn	b94823f8aa	add raman support	2024-11-24 22:13:29 +08:00
chn	28e4d29f2c	init	2024-11-21 16:18:06 +08:00