diff --git a/.gitignore b/.gitignore index 4ec538e..1098acc 100644 --- a/.gitignore +++ b/.gitignore @@ -5,3 +5,5 @@ main plot out.yaml .ccls-cache +.xmake/ +build/ diff --git a/flake.lock b/flake.lock index c871281..a62e708 100644 --- a/flake.lock +++ b/flake.lock @@ -807,11 +807,11 @@ "touchix": "touchix" }, "locked": { - "lastModified": 1695465859, - "narHash": "sha256-nauD552M00PmcpDkKqoYC0ql8QFaFyfKpIYm2qyRKEQ=", + "lastModified": 1695648380, + "narHash": "sha256-8zQct7OdcAtzwsvyPy07wImEtO6v0Rz//grrDtol+Jo=", "owner": "CHN-beta", "repo": "nixos", - "rev": "ad7be5bc2b22a47baf34eebb913b93c11482fb37", + "rev": "e21c7a916a66de6894c3e785a3a0d5f621458272", "type": "github" }, "original": { @@ -951,11 +951,11 @@ }, "nixpkgs-unstable": { "locked": { - "lastModified": 1695275987, - "narHash": "sha256-DjfwcdkKY7Hom5Nk8fR2fJInRf93uyeAvX+GoA1rqKo=", + "lastModified": 1695562663, + "narHash": "sha256-cnPZFTKQGldgFHWC9LHt8mDtJVTmj4WJxvKHF5uw5qA=", "owner": "CHN-beta", "repo": "nixpkgs", - "rev": "1523dd436e6f9baf1ab789131b4394b3502de789", + "rev": "82e02a1546fc59bc496f859e446ddcd412f09849", "type": "github" }, "original": { @@ -1013,6 +1013,22 @@ "type": "github" } }, + "nixpkgs_5": { + "locked": { + "lastModified": 1695360818, + "narHash": "sha256-JlkN3R/SSoMTa+CasbxS1gq+GpGxXQlNZRUh9+LIy/0=", + "owner": "NixOS", + "repo": "nixpkgs", + "rev": "e35dcc04a3853da485a396bdd332217d0ac9054f", + "type": "github" + }, + "original": { + "owner": "NixOS", + "ref": "nixos-unstable", + "repo": "nixpkgs", + "type": "github" + } + }, "nur": { "locked": { "lastModified": 1695048039, @@ -1143,7 +1159,8 @@ }, "root": { "inputs": { - "nixos": "nixos" + "nixos": "nixos", + "nixpkgs": "nixpkgs_5" } }, "sops-nix": { diff --git a/flake.nix b/flake.nix index 5b1070c..eb57543 100644 --- a/flake.nix +++ b/flake.nix @@ -1,16 +1,26 @@ { inputs.nixos.url = "github:CHN-beta/nixos"; + inputs.nixpkgs.url = "github:NixOS/nixpkgs/nixos-unstable"; - outputs = inputs: let pkgs = inputs.nixos.nixosConfigurations.pc.pkgs; in - { - devShell.x86_64-linux = pkgs.mkShell.override { stdenv = pkgs.genericPackages.gcc13Stdenv; } + outputs = inputs: + let + pkgs = inputs.nixpkgs.legacyPackages.x86_64-linux; + localPackages = inputs.nixos.nixosConfigurations.pc.pkgs.localPackages; + in { - buildInputs = with pkgs; - [ - yaml-cpp eigen fmt (localPackages.concurrencpp.override { stdenv = genericPackages.gcc13Stdenv; }) highfive - hdf5.dev tbb (localPackages.matplotplusplus.override { stdenv = genericPackages.gcc13Stdenv; }) - ]; - hardeningDisable = [ "all" ]; + devShell.x86_64-linux = pkgs.mkShell.override { stdenv = pkgs.stdenvNoCC; } + { + packages = with pkgs; [ xmake gcc13 pkg-config cmake ]; + inputsFrom = with pkgs; + [ + yaml-cpp eigen fmt localPackages.concurrencpp highfive + hdf5.dev tbb_2021_8.dev localPackages.matplotplusplus + localPackages.zpp-bits + ]; + PKG_CONFIG_PATH = "${pkgs.tbb_2021_8.dev}/lib/pkgconfig:${pkgs.yaml-cpp}/share/pkgconfig"; + yaml-cpp_DIR = "${pkgs.yaml-cpp}/share/cmake/yaml-cpp"; + hardeningDisable = [ "all" ]; + NIX_DEBUG = "1"; + }; }; - }; } diff --git a/include/ufo/ufo.hpp b/include/ufo/ufo.hpp new file mode 100644 index 0000000..21e3463 --- /dev/null +++ b/include/ufo/ufo.hpp @@ -0,0 +1,142 @@ +# 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(), 9)); } +} + +// 在相位中, 约定为使用 $\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; + // 超胞中原子的坐标，每行表示一个原子的坐标，单位为埃 + 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; + + struct InputOutputFile_ + { + std::string FileName; + std::string Format; + std::map ExtraParameters; + }; + + // 从哪个文件读入 AtomPosition, 以及这个文件的格式, 格式可选值包括 "yaml" + InputOutputFile_ AtomPositionInputFile; + // 从哪个文件读入 QPointData, 以及这个文件的格式, 格式可选值包括 "yaml" 和 "hdf5" + InputOutputFile_ QPointDataInputFile; + + // 输出到哪些文件, 以及使用怎样的格式, 格式可选值包括: + // 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); +}; + +concurrencpp::generator, unsigned>> + triplet_sequence(Eigen::Vector range); diff --git a/include/ufo/ufo.impl.hpp b/include/ufo/ufo.impl.hpp new file mode 100644 index 0000000..43bc9d4 --- /dev/null +++ b/include/ufo/ufo.impl.hpp @@ -0,0 +1,213 @@ +# 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) = node["SuperCellDeformation"][i][j].as(); + } + + for (unsigned i = 0; i < 3; i++) + PrimativeCellBasisNumber(i) = node["PrimativeCellBasisNumber"][i].as(); + + AtomPositionInputFile.FileName = node["AtomPositionInputFile"]["FileName"].as(); + AtomPositionInputFile.Format = node["AtomPositionInputFile"]["Format"].as(); + if (!std::set{"yaml"}.contains(AtomPositionInputFile.Format)) + throw std::runtime_error(fmt::format + ("Unknown format: {}, should be \"yaml\".", AtomPositionInputFile.Format)); + + QPointDataInputFile.FileName = node["QPointDataInputFile"]["FileName"].as(); + QPointDataInputFile.Format = node["QPointDataInputFile"]["Format"].as(); + if (!std::set{"yaml", "hdf5"}.contains(QPointDataInputFile.Format)) + throw std::runtime_error(fmt::format + ("Unknown format: {}, should be \"yaml\" or \"hdf5\".", QPointDataInputFile.Format)); + + if (auto value = node["QPointDataOutputFile"]) + for (unsigned i = 0; i < value.size(); i++) + { + auto& _ = QPointDataOutputFile.emplace_back(); + _.FileName = value[i]["FileName"].as(); + _.Format = value[i]["Format"].as(); + if (!std::set{"yaml", "hdf5"}.contains(_.Format)) + throw std::runtime_error(fmt::format + ("Unknown format: {}, should be \"yaml\" or \"hdf5\".", _.Format)); + } + } + + if (AtomPositionInputFile.Format == "yaml") + { + auto node = YAML::LoadFile(AtomPositionInputFile.FileName); + auto points = node["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; + 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++) + { + QPointData[i].QPoint.resize(3); + 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) +{ + 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(); + std::basic_ofstream(filename, std::ios::binary) + << std::basic_string_view{data.data(), data.size()}; + } +} + +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/main.cpp b/main.cpp deleted file mode 100644 index b659fdf..0000000 --- a/main.cpp +++ /dev/null @@ -1,440 +0,0 @@ -# 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() - { return {{"r", HighFive::AtomicType{}}, {"i", HighFive::AtomicType{}}}; } -HIGHFIVE_REGISTER_TYPE(PhonopyComplex, create_compound_complex) - -// 在相位中, 约定为使用 $\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; - Eigen::Matrix SuperCellDeformation; - // 在单胞内取几个平面波的基矢 - Eigen::Vector PrimativeCellBasisNumber; - // 超胞中原子的坐标，每行表示一个原子的坐标，单位为埃 - Eigen::MatrixX3d AtomPosition; - - // 是否调整输出结果, 使得结果中的模式适合人类阅读. 默认为 true. - // 这包括合并相近的模式, 去除权重过小的模式, 限制输出的小数位数. - // 如果想用结果来进一步画图, 则建议关闭. - std::optional Filter; - - // 关于各个 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; - - Input(std::string yaml_file, std::optional hdf5_file); -}; - -struct Output -{ - // 关于各个 Q 点的数据 - struct QPointDataType_ - { - // Q 点的坐标，单位为单胞的倒格矢 - Eigen::Vector3d QPoint; - - // 来源于哪个 Q 点, 单位为超胞的倒格矢 - Eigen::Vector3d Source; - - // 关于这个 Q 点上各个模式的数据 - struct ModeDataType_ - { - // 模式的频率，单位为 THz - double Frequency; - // 模式的权重 - double Weight; - }; - std::vector ModeData; - }; - std::vector QPointData; -}; - -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 - }; -} - -int main(int argc, const char** argv) -{ - if (argc < 3 || argc > 5) - throw std::runtime_error("Usage: " + std::string(argv[0]) + " input.yaml [input.hdf5] output.yaml"); - - std::cerr << "Reading input file..." << std::flush; - Input input(argv[1], argc > 3 ? std::make_optional(argv[2]) : std::nullopt); - 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_super_cell_in_modified_reciprocal_super_cell = [&] - { - 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()).eval(); - 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 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(); - // 这一步推导过程在计算 score 的函数中 - auto qpoint_by_reciprocal_modified_super_cell_in_modified_reciprocal_super_cell = - input.SuperCellDeformation.inverse() * qpoint_by_reciprocal_super_cell_in_modified_reciprocal_super_cell; - 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; - if (input.Filter.value_or(true)) - { - // 从小到大枚举所有的模式，并将相近的模式（相差小于 0.1 THz）合并 - std::map frequency_to_weight; - for (unsigned i_of_mode = 0; i_of_mode < input.QPointData[i_of_qpoint].ModeData.size(); i_of_mode++) - { - auto frequency = input.QPointData[i_of_qpoint].ModeData[i_of_mode].Frequency; - auto weight = projection_coefficient[i_of_qpoint][i_of_mode][i_of_sub_qpoint]; - 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; - } - } - // 仅保留权重大于 0.1 的模式 - for (auto& mode : frequency_to_weight) - if (mode.second > 0.1) - { - auto& __ = _.ModeData.emplace_back(); - __.Frequency = mode.first; - __.Weight = mode.second; - } - } - else - 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; - - // YAML 输出得太丑了，我来自己写 - std::cerr << "Writing output file..." << std::flush; - std::ofstream(argc > 3 ? argv[3] : argv[2]) << [&] - { - std::stringstream print; - auto format = input.Filter.value_or(true) ? 3 : 10; - print << "QPointData:\n"; - for (auto& qpoint: output.QPointData) - { - print << fmt::format(" - QPoint: [ {1:.{0}f}, {2:.{0}f}, {3:.{0}f} ]\n", - format, qpoint.QPoint[0], qpoint.QPoint[1], qpoint.QPoint[2]); - print << fmt::format(" Source: [ {1:.{0}f}, {2:.{0}f}, {3:.{0}f} ]\n", - format, 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", - format, mode.Frequency, mode.Weight); - } - return print.str(); - }(); - std::cerr << "Done." << std::endl; -} - -// 从文件中读取输入, 文件中应当包含: (大多数据可以直接从 phonopy 的输出中复制) -// 单胞的格矢: lattice 单位为埃 直接从 phonopy 的输出中复制 -// 超胞的倍数: SuperCellMultiplier 手动输入, 为一个包含三个整数的数组 -// 平面波的基矢个数: PrimativeCellBasisNumber 手动输入, 为一个包含三个整数的数组 -// 超胞中原子的坐标: points[*].coordinates 单位为超胞的格矢 直接从 phonopy 的输出中复制 -// 各个 Q 点的坐标: phonon[*].q-position 单位为超胞的倒格子的格矢 直接从 phonopy 的输出中复制 -// 各个模式的频率: phonon[*].band[*].frequency 单位为 THz 直接从 phonopy 的输出中复制 -// 各个模式的原子运动状态: phonon[*].band[*].eigenvector 直接从 phonopy 的输出中复制 -// 文件中可以有多余的项目, 多余的项目不管. -Input::Input(std::string yaml_file, std::optional hdf5_file) -{ - auto node = YAML::LoadFile(yaml_file); - for (unsigned i = 0; i < 3; i++) - for (unsigned j = 0; j < 3; j++) - PrimativeCell(i, j) = node["lattice"][i][j].as(); - - for (unsigned i = 0; i < 3; i++) - SuperCellMultiplier(i) = node["SuperCellMultiplier"][i].as(); - - for (unsigned i = 0; i < 3; i++) - for (unsigned j = 0; j < 3; j++) - SuperCellDeformation(i, j) = node["SuperCellDeformation"][i][j].as(); - - for (unsigned i = 0; i < 3; i++) - PrimativeCellBasisNumber(i) = node["PrimativeCellBasisNumber"][i].as(); - - if (auto value = node["Filter"]) - Filter = value.as(); - - auto points = node["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(); - AtomPosition = atom_position_to_super_cell - * (SuperCellDeformation * SuperCellMultiplier.cast().asDiagonal() * PrimativeCell); - - if (hdf5_file) - { - HighFive::File file(*hdf5_file, 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(); - } - } - } - else - { - auto phonon = node["phonon"].as>(); - QPointData.resize(phonon.size()); - for (unsigned i = 0; i < phonon.size(); i++) - { - QPointData[i].QPoint.resize(3); - 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(); - } - } - } -} diff --git a/src/ufo.cpp b/src/ufo.cpp new file mode 100644 index 0000000..286309b --- /dev/null +++ b/src/ufo.cpp @@ -0,0 +1,179 @@ +# 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_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()).eval(); + 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 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(); + // 这一步推导过程在计算 score 的函数中 + auto qpoint_by_reciprocal_modified_super_cell_in_modified_reciprocal_super_cell = + input.SuperCellDeformation.value_or(Eigen::Matrix3d::Identity()).inverse() + * qpoint_by_reciprocal_super_cell_in_modified_reciprocal_super_cell; + 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/xmake.lua b/xmake.lua new file mode 100644 index 0000000..365a08a --- /dev/null +++ b/xmake.lua @@ -0,0 +1,13 @@ +add_rules("mode.debug", "mode.release") +set_languages("cxx23") + +add_requires("pkgconfig::tbb") +add_requires("pkgconfig::yaml-cpp") + +target("ufo") + set_kind("binary") + add_includedirs("include") + add_files("src/ufo.cpp") + add_packages("tbb") + add_packages("yaml-cpp") +