diff --git a/CMakeLists.txt b/CMakeLists.txt index 4fa2902..e026da1 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -22,7 +22,7 @@ find_package(TBB REQUIRED) find_package(Matplot++ REQUIRED) find_path(ZPP_BITS_INCLUDE_DIR zpp_bits.h REQUIRED) -add_executable(ufo src/ufo.cpp) +add_executable(ufo src/unfold.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) diff --git a/include/ufo/solver.hpp b/include/ufo/solver.hpp new file mode 100644 index 0000000..2f72654 --- /dev/null +++ b/include/ufo/solver.hpp @@ -0,0 +1,70 @@ +# pragma once +# include +# include +# include +# include +# include +# include +# include +# include +# include +# include +# include +# include +# include +# include +# include +# include +# include +# include +# include +# include +# include +# include +# include + +// 在相位中, 约定为使用 $\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{}}, {"i", HighFive::AtomicType{}}}; } + class Solver + { + public: + virtual Solver& operator()() = 0; + ~Solver() = default; + + inline static concurrencpp::generator, unsigned>> + triplet_sequence(Eigen::Vector 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{{x}, {y}, {z}}, + x * range[1] * range[2] + y * range[2] + z + }; + } + }; +} + +HIGHFIVE_REGISTER_TYPE(ufo::PhonopyComplex, ufo::create_compound_complex) diff --git a/include/ufo/ufo.hpp b/include/ufo/ufo.hpp deleted file mode 100644 index 0fe548e..0000000 --- a/include/ufo/ufo.hpp +++ /dev/null @@ -1,155 +0,0 @@ -# include -# include -# include -# include -# include -# include -# include -# include -# include -# include -# include -# include -# include -# include -# include -# include -# include -# include -# include -# include -# include -# include -# include - -using namespace std::literals; - -struct PhonopyComplex { double r, i; }; -HighFive::CompoundType create_compound_complex(); -HIGHFIVE_REGISTER_TYPE(PhonopyComplex, create_compound_complex) - -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())); } -} - -// 在相位中, 约定为使用 $\exp (2 \pi i \vec{q} \cdot \vec{r})$ 来表示原子的运动状态 -// (而不是 $\exp (-2 \pi i \vec{q} \cdot \vec{r})$) -// 一些书定义的倒格矢中包含了 $2 \pi$ 的部分, 我们这里约定不包含这部分. -// 也就是说, 正格子与倒格子的转置相乘, 得到单位矩阵. - -struct Input -{ - // 单胞的三个格矢，每行表示一个格矢的坐标，单位为埃 - Eigen::Matrix3d PrimativeCell; - // 单胞到超胞的格矢转换时用到的矩阵 - // 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::Vector SuperCellMultiplier; - std::optional> SuperCellDeformation; - // 在单胞内取几个平面波的基矢 - Eigen::Vector PrimativeCellBasisNumber; - - struct InputOutputFile_ - { - std::string FileName; - std::string Format; - std::map ExtraParameters; - }; - - // 从哪个文件读入 AtomPosition, 以及这个文件的格式, 格式可选值包括 "yaml" - InputOutputFile_ AtomPositionInputFile; - // 从哪个文件读入 QPointData, 以及这个文件的格式, 格式可选值包括 "yaml" 和 "hdf5" - InputOutputFile_ 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 ModeData; - }; - std::vector QPointData; - - // 输出到哪些文件, 以及使用怎样的格式, 格式可选值包括: - // yaml: 使用 yaml 格式输出 - // yaml-human-readable: 使用 yaml 格式输出, 但是输出的结果更适合人类阅读, 包括合并相近的模式, 去除权重过小的模式, 限制输出的小数位数. - // zpp: 使用 zpp-bits 序列化, 可以直接被 plot.cpp 读取 - std::vector 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 的输出中复制 - // 文件中可以有多余的项目, 多余的项目不管. - Input(std::string filename); -}; - -struct Output -{ - // 关于各个 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 ModeData; - }; - std::vector QPointData; - - void write(std::string filename, std::string format, unsigned percision = 10) const; - Output() = default; - Output(std::string filename); - - using serialize = zpp::bits::members<1>; -}; - -concurrencpp::generator, unsigned>> - triplet_sequence(Eigen::Vector range); diff --git a/include/ufo/ufo.impl.hpp b/include/ufo/ufo.impl.hpp deleted file mode 100644 index 79f93e8..0000000 --- a/include/ufo/ufo.impl.hpp +++ /dev/null @@ -1,268 +0,0 @@ -# include - -inline HighFive::CompoundType create_compound_complex() - { return {{"r", HighFive::AtomicType{}}, {"i", HighFive::AtomicType{}}}; } - -inline Input::Input(std::string filename) -{ - // read main input file - { - auto node = YAML::LoadFile(filename); - for (unsigned i = 0; i < 3; i++) - for (unsigned j = 0; j < 3; j++) - PrimativeCell(i, j) = node["PrimativeCell"][i][j].as(); - - for (unsigned i = 0; i < 3; i++) - SuperCellMultiplier(i) = node["SuperCellMultiplier"][i].as(); - - if (auto value = node["SuperCellDeformation"]) - { - SuperCellDeformation.emplace(); - for (unsigned i = 0; i < 3; i++) - for (unsigned j = 0; j < 3; j++) - (*SuperCellDeformation)(i, j) = value[i][j].as(); - } - - for (unsigned i = 0; i < 3; i++) - PrimativeCellBasisNumber(i) = node["PrimativeCellBasisNumber"][i].as(); - - auto read_file_config = [filename](YAML::Node source, InputOutputFile_& config) - { - if (auto _ = source["SameAsConfigFile"]) - { - auto __ = _.as(); - config.ExtraParameters["SameAsConfigFile"] = __; - if (__) - { - config.FileName = filename; - config.Format = "yaml"; - return; - } - } - config.FileName = source["FileName"].as(); - config.Format = source["Format"].as(); - if (auto _ = source["RelativeToConfigFile"]) - { - auto __ = _.as(); - config.ExtraParameters["RelativeToConfigFile"] = __; - if (__) - config.FileName = std::filesystem::path(filename).parent_path() / config.FileName; - } - }; - read_file_config(node["AtomPositionInputFile"], AtomPositionInputFile); - if (!std::set{"yaml"}.contains(AtomPositionInputFile.Format)) - throw std::runtime_error(fmt::format - ("Unknown AtomPositionInputFile.Format: {}, should be \"yaml\".", AtomPositionInputFile.Format)); - read_file_config(node["QPointDataInputFile"], QPointDataInputFile); - if (!std::set{"yaml", "hdf5"}.contains(QPointDataInputFile.Format)) - throw std::runtime_error(fmt::format - ("Unknown QPointDataInputFile.Format: {}, should be \"yaml\" or \"hdf5\".", QPointDataInputFile.Format)); - if (auto value = node["QPointDataOutputFile"]) - { - QPointDataOutputFile.resize(value.size()); - for (unsigned i = 0; i < value.size(); i++) - { - read_file_config(value[i], QPointDataOutputFile[i]); - if - ( - QPointDataOutputFile[i].ExtraParameters.contains("SameAsConfigFile") - && std::any_cast(QPointDataOutputFile[i].ExtraParameters["SameAsConfigFile"]) - ) - throw std::runtime_error("QPointDataOutputFile.SameAsConfigFile should not be set."); - if - ( - !std::set{"yaml", "yaml-human-readable", "zpp"} - .contains(QPointDataOutputFile[i].Format) - ) - throw std::runtime_error(fmt::format - ( - "Unknown QPointDataOutputFile[{}].Format: {}, should be \"yaml\", \"yaml-human-readable\" or \"zpp\".", - i, QPointDataOutputFile[i].Format - )); - } - } - } - - if (AtomPositionInputFile.Format == "yaml") - { - auto node = YAML::LoadFile(AtomPositionInputFile.FileName); - std::vector points; - if (auto _ = node["points"]) - points = _.as>(); - else - points = node["unit_cell"]["points"].as>(); - 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(); - auto super_cell = (SuperCellDeformation.value_or(Eigen::Matrix3d::Identity()) - * SuperCellMultiplier.cast().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>(); - 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(); - auto band = phonon[i]["band"].as>(); - QPointData[i].ModeData.resize(band.size()); - for (unsigned j = 0; j < band.size(); j++) - { - QPointData[i].ModeData[j].Frequency = band[j]["frequency"].as(); - auto eigenvector_vectors = band[j]["eigenvector"] - .as>>>(); - 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 * 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") - { - HighFive::File file(QPointDataInputFile.FileName, HighFive::File::ReadOnly); - auto size = file.getDataSet("/frequency").getDimensions(); - auto frequency = file.getDataSet("/frequency") - .read>>>(); - auto eigenvector_vector = file.getDataSet("/eigenvector") - .read>>>>(); - auto path = file.getDataSet("/path") - .read>>>(); - 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(); - } - } - } -} - -inline void Output::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") - { - Output output; - std::map> - 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 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") - { - auto [data, out] = zpp::bits::data_out(); - out(*this).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(data.data()), data.size()); - } -} - -inline Output::Output(std::string filename) -{ - auto input = std::ifstream(filename, std::ios::binary | std::ios::in); - input.exceptions(std::ios::badbit | std::ios::failbit); - std::vector data; - { - std::vector string(std::istreambuf_iterator(input), {}); - data.assign - ( - reinterpret_cast(string.data()), - reinterpret_cast(string.data() + string.size()) - ); - } - auto in = zpp::bits::in(data); - in(*this).or_throw(); -} - -inline concurrencpp::generator, unsigned>> - triplet_sequence(Eigen::Vector 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{{x}, {y}, {z}}, - x * range[1] * range[2] + y * range[2] + z - }; -} diff --git a/include/ufo/unfold.hpp b/include/ufo/unfold.hpp new file mode 100644 index 0000000..6c21285 --- /dev/null +++ b/include/ufo/unfold.hpp @@ -0,0 +1,172 @@ +# pragma once +# include + +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 SuperCellMultiplier; + std::optional> SuperCellDeformation; + // 在单胞内取几个平面波的基矢 + Eigen::Vector PrimativeCellBasisNumber; + + struct InputOutputFile + { + std::string FileName; + std::string Format; + std::map ExtraParameters; + }; + + // 从哪个文件读入 AtomPosition, 以及这个文件的格式, 格式可选值包括 "yaml" + InputOutputFile AtomPositionInputFile; + // 从哪个文件读入 QPointData, 以及这个文件的格式, 格式可选值包括 "yaml" 和 "hdf5" + InputOutputFile 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 ModeData; + }; + std::vector QPointData; + + // 输出到哪些文件, 以及使用怎样的格式, 格式可选值包括: + // yaml: 使用 yaml 格式输出 + // yaml-human-readable: 使用 yaml 格式输出, 但是输出的结果更适合人类阅读, + // 包括合并相近的模式, 去除权重过小的模式, 限制输出的小数位数. + // zpp: 使用 zpp-bits 序列化, 可以直接被 plot.cpp 读取 + std::vector 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 ModeData; + }; + std::vector 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>; + protected: + InputType Input_; + std::optional Output_; + std::optional Basis_; + + // 第一层是不同的模式, 第二层是不同的 sub qpoint + using ProjectionCoefficientType_ = std::vector>; + + 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>& mode_data, + std::atomic& number_of_finished_modes + ); + + OutputType construct_output + ( + const decltype(InputType::PrimativeCell)& primative_cell, + const decltype(InputType::SuperCellMultiplier)& super_cell_multiplier, + const decltype(InputType::SuperCellDeformation)& super_cell_deformation, + const std::vector>& qpoint, + const std::vector>>& frequency, + const ProjectionCoefficientType_& projection_coefficient + ); + }; +} diff --git a/src/main.cpp b/src/main.cpp new file mode 100644 index 0000000..638944f --- /dev/null +++ b/src/main.cpp @@ -0,0 +1,9 @@ +# include + +int main(int argc, const char** argv) +{ + if (argc != 2) + throw std::runtime_error(fmt::format("Usage: {} config.yaml", argv[0])); + + ufo::UnfoldSolver{argv[1]}(); +} diff --git a/src/ufo.cpp b/src/ufo.cpp deleted file mode 100644 index aec0a48..0000000 --- a/src/ufo.cpp +++ /dev/null @@ -1,175 +0,0 @@ -# include - -int main(int argc, const char** argv) -{ - if (argc != 2) - throw std::runtime_error(fmt::format("Usage: {} config.yaml", argv[0])); - - std::cerr << "Reading input file..." << std::flush; - Input input(argv[1]); - std::cerr << " Done." << std::endl; - - // 反折叠的原理: 将超胞中的原子运动状态, 投影到一组平面波构成的基矢中. - // 每一个平面波的波矢由两部分相加得到: 一部分是单胞倒格子的整数倍, 所取的个数有一定任意性, 论文中建议取大约单胞中原子个数那么多个; - // 对于没有缺陷的情况, 取一个应该就足够了. - // 另一部分是超胞倒格子的整数倍, 取 n 个, n 为超胞对应的单胞的倍数, 其实也就是倒空间中单胞对应倒格子中超胞的格点. - // 只要第一部分取得足够多, 那么单胞中原子的状态就可以完全被这些平面波描述. - // 将超胞中原子的运动状态投影到这些基矢上, 计算出投影的系数, 就可以将超胞的原子运动状态分解到单胞中的多个 q 点上. - - // 构建基 - // 每个 q 点对应的一组 sub qpoint。不同的 q 点所对应的 sub qpoint 是不一样的，但 sub qpoint 与 q 点的相对位置一致。 - // 这里 xyz_of_diff_of_sub_qpoint 即表示这个相对位置。 - // 由于基只与这个相对位置有关（也就是说，不同 q 点的基是一样的），因此可以先计算出所有的基，这样降低计算量。 - // 外层下标对应超胞倒格子的整数倍那部分(第二部分), 也就是不同的 sub qpoint - // 内层下标对应单胞倒格子的整数倍那部分(第一部分), 也就是 sub qpoint 上的不同平面波（取的数量越多，结果越精确） - std::cerr << "Calculating basis..." << std::flush; - std::vector> basis(input.SuperCellMultiplier.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(input.SuperCellMultiplier)) - { - basis[i_of_sub_qpoint].resize(input.PrimativeCellBasisNumber.prod()); - for (auto [xyz_of_basis, i_of_basis] : triplet_sequence(input.PrimativeCellBasisNumber)) - { - // 计算 q 点的坐标, 单位为单胞的倒格矢 - auto diff_of_sub_qpoint_by_reciprocal_primative_cell = xyz_of_basis.cast() - + input.SuperCellMultiplier.cast().cwiseInverse().asDiagonal() - * xyz_of_diff_of_sub_qpoint_by_reciprocal_modified_super_cell.cast(); - // 将 q 点坐标转换为埃^-1 - auto qpoint = (diff_of_sub_qpoint_by_reciprocal_primative_cell.transpose() - * (input.PrimativeCell.transpose().inverse())).transpose(); - // 计算基矢 - basis[i_of_sub_qpoint][i_of_basis] - = (2i * std::numbers::pi_v * (input.AtomPosition * qpoint)).array().exp(); - } - } - std::cerr << "Done." << std::endl; - - // 计算投影的结果 - // 最外层下标对应反折叠前的 q 点, 第二层下标对应不同模式, 第三层下标对应这个模式在反折叠后的 q 点(sub qpoint) - std::vector>> projection_coefficient(input.QPointData.size()); - std::atomic finished_qpoint(0); - // 对每个 q 点并行 - std::transform - ( - std::execution::par, input.QPointData.begin(), input.QPointData.end(), - projection_coefficient.begin(), [&](const auto& qpoint_data) - { - std::osyncstream(std::cerr) << fmt::format("\rCalculating projection coefficient... ({}/{})", - finished_qpoint, input.QPointData.size()) << std::flush; - std::vector> projection_coefficient(qpoint_data.ModeData.size()); - // 这里, qpoint_data 和 projection_coefficient 均指对应于一个 q 点的数据 - for (unsigned i_of_mode = 0; i_of_mode < qpoint_data.ModeData.size(); i_of_mode++) - { - auto& _ = projection_coefficient[i_of_mode]; - _.resize(input.SuperCellMultiplier.prod()); - for (unsigned i_of_sub_qpoint = 0; i_of_sub_qpoint < input.SuperCellMultiplier.prod(); i_of_sub_qpoint++) - // 对于 basis 中, 对应于单胞倒格子的部分, 以及对应于不同方向的部分, 分别求内积, 然后求模方和 - for (unsigned i_of_basis = 0; i_of_basis < input.PrimativeCellBasisNumber.prod(); i_of_basis++) - _[i_of_sub_qpoint] += - (basis[i_of_sub_qpoint][i_of_basis].transpose().conjugate() * qpoint_data.ModeData[i_of_mode].AtomMovement) - .array().abs2().sum(); - // 如果是严格地将向量分解到一组完备的基矢上, 那么不需要对计算得到的权重再做归一化处理 - // 但这里并不是这样一个严格的概念. 因此对分解到各个 sub qpoint 上的权重做归一化处理 - auto sum = std::accumulate(_.begin(), _.end(), 0.); - for (auto& __ : _) - __ /= sum; - } - finished_qpoint++; - std::osyncstream(std::cerr) << fmt::format("\rCalculating projection coefficient...({}/{})", - finished_qpoint, input.QPointData.size()) << std::flush; - return projection_coefficient; - }); - std::cerr << " Done." << std::endl; - - // 填充输出对象 - std::cerr << "Filling output object..." << std::flush; - Output output; - for (unsigned i_of_qpoint = 0; i_of_qpoint < input.QPointData.size(); i_of_qpoint++) - { - // 当 SuperCellDeformation 不是单位矩阵时, input.QPointData[i_of_qpoint].QPoint 不一定在 reciprocal_primative_cell 中 - // 需要首先将 q 点平移数个周期, 进入不包含 SuperCellDeformation 的超胞的倒格子中 - auto qpoint_by_reciprocal_modified_super_cell_in_modified_reciprocal_super_cell - = !input.SuperCellDeformation ? input.QPointData[i_of_qpoint].QPoint : [&] - { - auto current_qpoint = input.QPointData[i_of_qpoint].QPoint; - // 给一个 q 点打分 - // 计算这个 q 点以 modified_reciprocal_supre_cell 为单位的坐标, 依次考虑每个维度, 总分为每个维度之和. - // 如果这个坐标大于 0 小于 1, 则打 0 分. - // 如果这个坐标小于 0, 则打这个坐标的相反数分. - // 如果这个坐标大于 1, 则打这个坐标减去 1 的分. - auto score = [&](Eigen::Vector3d qpoint_by_reciprocal_super_cell) - { - // SuperCell = SuperCellDeformation * SuperCellMultiplier.asDiagonal() * PrimativeCell - // ModifiedSuperCell = SuperCellMultiplier.asDiagonal() * PrimativeCell - // ReciprocalSuperCell = SuperCell.inverse().transpose() - // ReciprocalModifiedSuperCell = ModifiedSuperCell.inverse().transpose() - // qpoint.transpose() = qpoint_by_reciprocal_super_cell.transpose() * ReciprocalSuperCell - // qpoint.transpose() = qpoint_by_reciprocal_modified_super_cell.transpose() * ReciprocalModifiedSuperCell - auto qpoint_by_reciprocal_modified_super_cell = - (input.SuperCellDeformation->inverse() * qpoint_by_reciprocal_super_cell).eval(); - double score = 0; - for (unsigned i = 0; i < 3; i++) - { - auto coordinate = qpoint_by_reciprocal_modified_super_cell[i]; - if (coordinate < 0) - score -= coordinate; - else if (coordinate > 1) - score += coordinate - 1; - } - return score; - }; - while (score(current_qpoint) > 0) - { - double min_score = std::numeric_limits::max(); - Eigen::Vector3d min_score_qpoint; - for (int x = -1; x <= 1; x++) - for (int y = -1; y <= 1; y++) - for (int z = -1; z <= 1; z++) - { - auto this_qpoint = current_qpoint - + Eigen::Matrix{{x}, {y}, {z}}.cast(); - auto this_score = score(this_qpoint); - if (this_score < min_score) - { - min_score = this_score; - min_score_qpoint = this_qpoint; - } - } - current_qpoint = min_score_qpoint; - } - return input.SuperCellDeformation->inverse() * current_qpoint; - }(); - for (auto [xyz_of_diff_of_sub_qpoint_by_reciprocal_modified_super_cell, i_of_sub_qpoint] - : triplet_sequence(input.SuperCellMultiplier)) - { - auto& _ = output.QPointData.emplace_back(); - auto reciprocal_modified_super_cell = - (input.SuperCellMultiplier.cast().asDiagonal() * input.PrimativeCell).inverse().transpose(); - // sub qpoint 的坐标，单位为埃^-1 - auto sub_qpoint = ((xyz_of_diff_of_sub_qpoint_by_reciprocal_modified_super_cell.cast() - + qpoint_by_reciprocal_modified_super_cell_in_modified_reciprocal_super_cell) - .transpose() * reciprocal_modified_super_cell).transpose(); - // 将坐标转换为相对于单胞的倒格矢的坐标并写入 - // 由 sub_qpoint.transpose() = sub_qpoint_by_reciprocal_primative_cell.transpose() - // * PrimativeCell.transpose().inverse() - // 得到 sub_qpoint_by_reciprocal_primative_cell = PrimativeCell * sub_qpoint - _.QPoint = input.PrimativeCell * sub_qpoint; - _.Source = input.QPointData[i_of_qpoint].QPoint; - _.SourceIndex_ = i_of_qpoint; - for (unsigned i_of_mode = 0; i_of_mode < input.QPointData[i_of_qpoint].ModeData.size(); i_of_mode++) - { - auto& __ = _.ModeData.emplace_back(); - __.Frequency = input.QPointData[i_of_qpoint].ModeData[i_of_mode].Frequency; - __.Weight = projection_coefficient[i_of_qpoint][i_of_mode][i_of_sub_qpoint]; - } - } - } - std::cerr << " Done." << std::endl; - - std::cerr << "Writing output file..." << std::flush; - for (auto& output_file : input.QPointDataOutputFile) - output.write(output_file.FileName, output_file.Format); - std::cerr << " Done." << std::endl; -} diff --git a/src/unfold.cpp b/src/unfold.cpp new file mode 100644 index 0000000..af60b50 --- /dev/null +++ b/src/unfold.cpp @@ -0,0 +1,498 @@ +# include + +namespace ufo +{ + UnfoldSolver::InputType::InputType(std::string filename) + { + // read main input file + { + auto node = YAML::LoadFile(filename); + for (unsigned i = 0; i < 3; i++) + for (unsigned j = 0; j < 3; j++) + PrimativeCell(i, j) = node["PrimativeCell"][i][j].as(); + + for (unsigned i = 0; i < 3; i++) + SuperCellMultiplier(i) = node["SuperCellMultiplier"][i].as(); + + if (auto value = node["SuperCellDeformation"]) + { + SuperCellDeformation.emplace(); + for (unsigned i = 0; i < 3; i++) + for (unsigned j = 0; j < 3; j++) + (*SuperCellDeformation)(i, j) = value[i][j].as(); + } + + for (unsigned i = 0; i < 3; i++) + PrimativeCellBasisNumber(i) = node["PrimativeCellBasisNumber"][i].as(); + + auto read_file_config = [filename](YAML::Node source, InputOutputFile& config) + { + if (auto _ = source["SameAsConfigFile"]) + { + auto __ = _.as(); + config.ExtraParameters["SameAsConfigFile"] = __; + if (__) + { + config.FileName = filename; + config.Format = "yaml"; + return; + } + } + config.FileName = source["FileName"].as(); + config.Format = source["Format"].as(); + if (auto _ = source["RelativeToConfigFile"]) + { + auto __ = _.as(); + config.ExtraParameters["RelativeToConfigFile"] = __; + if (__) + config.FileName = std::filesystem::path(filename).parent_path() / config.FileName; + } + }; + read_file_config(node["AtomPositionInputFile"], AtomPositionInputFile); + if (!std::set{"yaml"}.contains(AtomPositionInputFile.Format)) + throw std::runtime_error(fmt::format + ("Unknown AtomPositionInputFile.Format: {}, should be \"yaml\".", AtomPositionInputFile.Format)); + read_file_config(node["QPointDataInputFile"], QPointDataInputFile); + if (!std::set{"yaml", "hdf5"}.contains(QPointDataInputFile.Format)) + throw std::runtime_error(fmt::format + ("Unknown QPointDataInputFile.Format: {}, should be \"yaml\" or \"hdf5\".", QPointDataInputFile.Format)); + if (auto value = node["QPointDataOutputFile"]) + { + QPointDataOutputFile.resize(value.size()); + for (unsigned i = 0; i < value.size(); i++) + { + read_file_config(value[i], QPointDataOutputFile[i]); + if + ( + QPointDataOutputFile[i].ExtraParameters.contains("SameAsConfigFile") + && std::any_cast(QPointDataOutputFile[i].ExtraParameters["SameAsConfigFile"]) + ) + throw std::runtime_error("QPointDataOutputFile.SameAsConfigFile should not be set."); + if + ( + !std::set{"yaml", "yaml-human-readable", "zpp"} + .contains(QPointDataOutputFile[i].Format) + ) + throw std::runtime_error(fmt::format + ( + "Unknown QPointDataOutputFile[{}].Format: {}, should be \"yaml\", \"yaml-human-readable\" or \"zpp\".", + i, QPointDataOutputFile[i].Format + )); + } + } + } + + if (AtomPositionInputFile.Format == "yaml") + { + auto node = YAML::LoadFile(AtomPositionInputFile.FileName); + std::vector points; + if (auto _ = node["points"]) + points = _.as>(); + else + points = node["unit_cell"]["points"].as>(); + 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(); + auto super_cell = (SuperCellDeformation.value_or(Eigen::Matrix3d::Identity()) + * SuperCellMultiplier.cast().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>(); + 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(); + auto band = phonon[i]["band"].as>(); + QPointData[i].ModeData.resize(band.size()); + for (unsigned j = 0; j < band.size(); j++) + { + QPointData[i].ModeData[j].Frequency = band[j]["frequency"].as(); + auto eigenvector_vectors = band[j]["eigenvector"] + .as>>>(); + 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 * 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") + { + HighFive::File file(QPointDataInputFile.FileName, HighFive::File::ReadOnly); + auto size = file.getDataSet("/frequency").getDimensions(); + auto frequency = file.getDataSet("/frequency") + .read>>>(); + auto eigenvector_vector = file.getDataSet("/eigenvector") + .read>>>>(); + auto path = file.getDataSet("/path") + .read>>>(); + QPointData.resize(size[0] * size[1]); + for (unsigned i = 0; i < size[0]; i++) + for (unsigned j = 0; j < size[1]; j++) + { + QPointData[i * size[1] + j].QPoint = Eigen::Vector3d(path[i][j].data()); + QPointData[i * size[1] + j].ModeData.resize(size[2]); + for (unsigned k = 0; k < size[2]; k++) + { + QPointData[i * size[1] + j].ModeData[k].Frequency = frequency[i][j][k]; + Eigen::MatrixX3cd eigenvectors(AtomPosition.rows(), 3); + for (unsigned l = 0; l < AtomPosition.rows(); l++) + for (unsigned m = 0; m < 3; m++) + eigenvectors(l, m) + = eigenvector_vector[i][j][l * 3 + m][k].r + eigenvector_vector[i][j][l * 3 + m][k].i * 1i; + QPointData[i * size[1] + j].ModeData[k].AtomMovement = eigenvectors / eigenvectors.norm(); + } + } + } + } + + void UnfoldSolver::OutputType::write + (decltype(InputType::QPointDataOutputFile) output_files) const + { + for (auto& output_file : output_files) + write(output_file.FileName, output_file.Format); + } + void UnfoldSolver::OutputType::write(std::string filename, std::string format, unsigned percision) const + { + if (format == "yaml") + std::ofstream(filename) << [&] + { + std::stringstream print; + print << "QPointData:\n"; + for (auto& qpoint: QPointData) + { + print << fmt::format(" - QPoint: [ {1:.{0}f}, {2:.{0}f}, {3:.{0}f} ]\n", + percision, qpoint.QPoint[0], qpoint.QPoint[1], qpoint.QPoint[2]); + print << fmt::format(" Source: [ {1:.{0}f}, {2:.{0}f}, {3:.{0}f} ]\n", + percision, qpoint.Source[0], qpoint.Source[1], qpoint.Source[2]); + print << " ModeData:\n"; + for (auto& mode: qpoint.ModeData) + print << fmt::format(" - {{ Frequency: {1:.{0}f}, Weight: {2:.{0}f} }}\n", + percision, mode.Frequency, mode.Weight); + } + return print.str(); + }(); + else if (format == "yaml-human-readable") + { + std::remove_cvref_t output; + std::map> + 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 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") + { + auto [data, out] = zpp::bits::data_out(); + out(*this).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(data.data()), data.size()); + } + } + + 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> mode_data; + for (auto& qpoint : Input_.QPointData) + for (auto& mode : qpoint.ModeData) + mode_data.emplace_back(mode.AtomMovement); + std::atomic 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 << "\rCalculating projection coefficient... Done." << std::endl; + + std::clog << "Constructing output... " << std::flush; + std::vector> qpoint; + std::vector>> 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_.PrimativeCell, Input_.SuperCellMultiplier, + Input_.SuperCellDeformation, qpoint, frequency, projection_coefficient + ); + std::clog << "Done." << std::endl; + } + std::clog << "Writing output... " << std::flush; + Output_->write(Input_.QPointDataOutputFile); + std::clog << "Done." << std::endl; + return *this; + } + + UnfoldSolver::BasisType UnfoldSolver::construct_basis + ( + const decltype(InputType::PrimativeCell)& primative_cell, + const decltype(InputType::SuperCellMultiplier)& super_cell_multiplier, + const decltype(InputType::PrimativeCellBasisNumber)& primative_cell_basis_number, + const decltype(InputType::AtomPosition)& atom_position + ) + { + BasisType basis(super_cell_multiplier.prod()); + // 每个 q 点对应的一组 sub qpoint。不同的 q 点所对应的 sub qpoint 是不一样的，但 sub qpoint 与 q 点的相对位置一致。 + // 这里 xyz_of_diff_of_sub_qpoint 即表示这个相对位置，单位为超胞的倒格矢 + for (auto [xyz_of_diff_of_sub_qpoint_by_reciprocal_modified_super_cell, i_of_sub_qpoint] + : triplet_sequence(super_cell_multiplier)) + { + basis[i_of_sub_qpoint].resize(primative_cell_basis_number.prod()); + for (auto [xyz_of_basis, i_of_basis] : triplet_sequence(primative_cell_basis_number)) + { + // 计算 q 点的坐标, 单位为单胞的倒格矢 + auto diff_of_sub_qpoint_by_reciprocal_primative_cell = xyz_of_basis.cast() + + super_cell_multiplier.cast().cwiseInverse().asDiagonal() + * xyz_of_diff_of_sub_qpoint_by_reciprocal_modified_super_cell.cast(); + // 将 q 点坐标转换为埃^-1 + auto qpoint = (diff_of_sub_qpoint_by_reciprocal_primative_cell.transpose() + * (primative_cell.transpose().inverse())).transpose(); + // 计算基矢 + basis[i_of_sub_qpoint][i_of_basis] + = (2i * std::numbers::pi_v * (atom_position * qpoint)).array().exp(); + } + } + return basis; + } + + std::vector> UnfoldSolver::construct_projection_coefficient + ( + const BasisType& basis, + const std::vector>& mode_data, + std::atomic& number_of_finished_modes + ) + { + // 第一层下标对应不同模式, 第二层下标对应这个模式在反折叠后的 q 点(sub qpoint) + std::vector> projection_coefficient(mode_data.size()); + // 对每个模式并行 + std::transform + ( + std::execution::par, mode_data.begin(), mode_data.end(), projection_coefficient.begin(), + [&](const auto& mode_data) + { + // 这里, mode_data 和 projection_coefficient 均指对应于一个模式的数据 + std::vector projection_coefficient(basis.size()); + for (unsigned i_of_sub_qpoint = 0; i_of_sub_qpoint < basis.size(); i_of_sub_qpoint++) + // 对于 basis 中, 对应于单胞倒格子的部分, 以及对应于不同方向的部分, 分别求内积, 然后求模方和 + for (unsigned i_of_basis = 0; i_of_basis < basis[i_of_sub_qpoint].size(); i_of_basis++) + projection_coefficient[i_of_sub_qpoint] += + (basis[i_of_sub_qpoint][i_of_basis].transpose().conjugate() * mode_data.get()) + .array().abs2().sum(); + // 如果是严格地将向量分解到一组完备的基矢上, 那么不需要对计算得到的权重再做归一化处理 + // 但这里并不是这样一个严格的概念. 因此对分解到各个 sub qpoint 上的权重做归一化处理 + auto sum = std::accumulate + (projection_coefficient.begin(), projection_coefficient.end(), 0.); + for (auto& _ : projection_coefficient) + _ /= sum; + number_of_finished_modes++; + return projection_coefficient; + } + ); + return projection_coefficient; + } + + UnfoldSolver::OutputType UnfoldSolver::construct_output + ( + const decltype(InputType::PrimativeCell)& primative_cell, + const decltype(InputType::SuperCellMultiplier)& super_cell_multiplier, + const decltype(InputType::SuperCellDeformation)& super_cell_deformation, + const std::vector>& qpoint, + const std::vector>>& frequency, + const ProjectionCoefficientType_& projection_coefficient + ) + { + OutputType output; + for (unsigned i_of_qpoint = 0, num_of_mode_in_all_qpoint = 0; i_of_qpoint < qpoint.size(); i_of_qpoint++) + { + // 当 SuperCellDeformation 不是单位矩阵时, input.QPointData[i_of_qpoint].QPoint 不一定在 reciprocal_primative_cell 中 + // 需要首先将 q 点平移数个周期, 进入不包含 SuperCellDeformation 的超胞的倒格子中 + auto qpoint_by_reciprocal_modified_super_cell_in_modified_reciprocal_super_cell + = !super_cell_deformation ? qpoint[i_of_qpoint].get() : [&] + { + auto current_qpoint = qpoint[i_of_qpoint].get(); + // 给一个 q 点打分 + // 计算这个 q 点以 modified_reciprocal_supre_cell 为单位的坐标, 依次考虑每个维度, 总分为每个维度之和. + // 如果这个坐标大于 0 小于 1, 则打 0 分. + // 如果这个坐标小于 0, 则打这个坐标的相反数分. + // 如果这个坐标大于 1, 则打这个坐标减去 1 的分. + auto score = [&](Eigen::Vector3d qpoint_by_reciprocal_super_cell) + { + // SuperCell = SuperCellDeformation * SuperCellMultiplier.asDiagonal() * PrimativeCell + // ModifiedSuperCell = SuperCellMultiplier.asDiagonal() * PrimativeCell + // ReciprocalSuperCell = SuperCell.inverse().transpose() + // ReciprocalModifiedSuperCell = ModifiedSuperCell.inverse().transpose() + // qpoint.transpose() = qpoint_by_reciprocal_super_cell.transpose() * ReciprocalSuperCell + // qpoint.transpose() = qpoint_by_reciprocal_modified_super_cell.transpose() * ReciprocalModifiedSuperCell + auto qpoint_by_reciprocal_modified_super_cell = + (super_cell_deformation->inverse() * qpoint_by_reciprocal_super_cell).eval(); + double score = 0; + for (unsigned i = 0; i < 3; i++) + { + auto coordinate = qpoint_by_reciprocal_modified_super_cell[i]; + if (coordinate < 0) + score -= coordinate; + else if (coordinate > 1) + score += coordinate - 1; + } + return score; + }; + while (score(current_qpoint) > 0) + { + double min_score = std::numeric_limits::max(); + Eigen::Vector3d min_score_qpoint; + for (int x = -1; x <= 1; x++) + for (int y = -1; y <= 1; y++) + for (int z = -1; z <= 1; z++) + { + auto this_qpoint = current_qpoint + + Eigen::Matrix{{x}, {y}, {z}}.cast(); + auto this_score = score(this_qpoint); + if (this_score < min_score) + { + min_score = this_score; + min_score_qpoint = this_qpoint; + } + } + current_qpoint = min_score_qpoint; + } + return super_cell_deformation->inverse() * current_qpoint; + }(); + for (auto [xyz_of_diff_of_sub_qpoint_by_reciprocal_modified_super_cell, i_of_sub_qpoint] + : triplet_sequence(super_cell_multiplier)) + { + auto& _ = output.QPointData.emplace_back(); + auto reciprocal_modified_super_cell = + (super_cell_multiplier.cast().asDiagonal() * primative_cell).inverse().transpose(); + // sub qpoint 的坐标，单位为埃^-1 + auto sub_qpoint = ((xyz_of_diff_of_sub_qpoint_by_reciprocal_modified_super_cell.cast() + + qpoint_by_reciprocal_modified_super_cell_in_modified_reciprocal_super_cell) + .transpose() * reciprocal_modified_super_cell).transpose(); + // 将坐标转换为相对于单胞的倒格矢的坐标并写入 + // 由 sub_qpoint.transpose() = sub_qpoint_by_reciprocal_primative_cell.transpose() + // * PrimativeCell.transpose().inverse() + // 得到 sub_qpoint_by_reciprocal_primative_cell = PrimativeCell * sub_qpoint + _.QPoint = primative_cell * sub_qpoint; + _.Source = qpoint[i_of_qpoint]; + _.SourceIndex_ = i_of_qpoint; + for (unsigned i_of_mode = 0; i_of_mode < frequency[i_of_qpoint].size(); i_of_mode++, num_of_mode_in_all_qpoint++) + { + auto& __ = _.ModeData.emplace_back(); + __.Frequency = frequency[i_of_qpoint][i_of_mode]; + __.Weight = projection_coefficient[num_of_mode_in_all_qpoint][i_of_sub_qpoint]; + } + } + } + return output; + } +} + +// inline Output::Output(std::string filename) +// { +// auto input = std::ifstream(filename, std::ios::binary | std::ios::in); +// input.exceptions(std::ios::badbit | std::ios::failbit); +// std::vector data; +// { +// std::vector string(std::istreambuf_iterator(input), {}); +// data.assign +// ( +// reinterpret_cast(string.data()), +// reinterpret_cast(string.data() + string.size()) +// ); +// } +// auto in = zpp::bits::in(data); +// in(*this).or_throw(); +// } +