From eb6b9034175ce07a7f2398a3998b08964428cfa1 Mon Sep 17 00:00:00 2001 From: chn Date: Tue, 3 Dec 2024 16:25:17 +0800 Subject: [PATCH] =?UTF-8?q?=E5=86=99=E4=BA=86=20project=5Fto=5Fmode=20?= =?UTF-8?q?=E4=BD=86=E8=BF=98=E6=B2=A1=E8=B0=83=E8=AF=95=E5=A5=BD?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit --- CMakeLists.txt | 5 +- include/ufo.hpp | 1 + src/main.cpp | 1 + src/project-to-mode.cpp | 151 +++++++++++++++++++++++++++++++ src/unfold.cpp | 54 +++++++++-- test/project-to-mode/.gitignore | 1 + test/project-to-mode/config.yaml | 6 ++ test/project-to-mode/data.zpp | 3 + test/zpp-to-yaml/data.zpp | 4 +- 9 files changed, 215 insertions(+), 11 deletions(-) create mode 100644 src/project-to-mode.cpp create mode 100644 test/project-to-mode/.gitignore create mode 100644 test/project-to-mode/config.yaml create mode 100644 test/project-to-mode/data.zpp diff --git a/CMakeLists.txt b/CMakeLists.txt index 9ea57ee..3008a16 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -14,7 +14,8 @@ find_package(Matplot++ REQUIRED) find_package(biu REQUIRED) find_package(Threads REQUIRED) -add_executable(ufo src/fold.cpp src/unfold.cpp src/plot.cpp src/raman.cpp src/project-to-atom.cpp src/main.cpp) +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) @@ -37,3 +38,5 @@ add_test(NAME fold WORKING_DIRECTORY ${PROJECT_SOURCE_DIR}/test/fold COMMAND ufo 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) diff --git a/include/ufo.hpp b/include/ufo.hpp index 68424f6..7bde046 100644 --- a/include/ufo.hpp +++ b/include/ufo.hpp @@ -94,6 +94,7 @@ namespace ufo 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); diff --git a/src/main.cpp b/src/main.cpp index bb84c94..4afe529 100644 --- a/src/main.cpp +++ b/src/main.cpp @@ -8,6 +8,7 @@ int main(int argc, const char** argv) else if (argv[1] == "zpp-to-yaml"s) std::cout << YAML::Node(biu::deserialize (biu::read(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 ( diff --git a/src/project-to-mode.cpp b/src/project-to-mode.cpp new file mode 100644 index 0000000..0737b71 --- /dev/null +++ b/src/project-to-mode.cpp @@ -0,0 +1,151 @@ +# include + +void ufo::project_to_mode(std::string config_file) +{ + // 将选定的超胞中的某一个 q 点上的每一个模式到反折叠后的某个 q 点的投影，投射到单胞中的对应 q 点的模式上 + + struct Config + { + // 在单胞内取几个平面波的基矢 + Eigen::Vector 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> Coefficient; + }; + + biu::Logger::Guard log(config_file); + auto config = YAML::LoadFile(config_file).as(); + auto input = biu::deserialize + (biu::read(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>(input.Super.Qpoint[config.SuperQpointIndex].Mode.size(), + std::vector(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 basis(config.PrimativeCellBasisNumber.prod()); + for (auto [xyz_of_basis, i_of_basis] + : biu::sequence(config.PrimativeCellBasisNumber)) + { + // 波矢就是 xyz_of_basis, 单位为单胞的倒格矢 + // 将它的单位转换成埃^-1 + auto wavevector = input.Primative.Cell.inverse() * xyz_of_basis.cast(); + log.debug("xyz_of_basis: {}"_f(xyz_of_basis)); + log.debug("wavevector: {:.3f}"_f(fmt::join(wavevector, ", "))); + // 将原子坐标单位也转换为埃 + auto atom_position = input.Primative.AtomPosition * input.Primative.Cell; + // 计算基矢 + basis[i_of_basis] = (2i * std::numbers::pi_v * (atom_position * wavevector)).array().exp(); + } + return basis; + }(); + // 基在超胞中的表示 + auto super_basis = [&] + { + biu::Logger::Guard log; + std::vector basis(config.PrimativeCellBasisNumber.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(); + }(); + log.debug("diff_of_sub_qpoint_by_reciprocal_modified_super_cell: {}"_f + (diff_of_sub_qpoint_by_reciprocal_modified_super_cell)); + for (auto [xyz_of_basis, i_of_basis] + : biu::sequence(config.PrimativeCellBasisNumber)) + { + // 计算波矢, 单位为单胞的倒格矢 + auto wavevector_by_reciprocal_primative_cell = xyz_of_basis.cast() + + input.Super.CellMultiplier.cast().cwiseInverse().asDiagonal() + * diff_of_sub_qpoint_by_reciprocal_modified_super_cell.cast(); + // 将单位转换为埃^-1 + auto wavevector = input.Primative.Cell.inverse() * wavevector_by_reciprocal_primative_cell; + // 将原子坐标单位也转换为埃 + auto atom_position = input.Super.AtomPosition + * (input.Super.CellDeformation * input.Super.CellMultiplier.cast().asDiagonal() * input.Primative.Cell); + // 计算基矢 + basis[i_of_basis] = (2i * std::numbers::pi_v * (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; + for (auto i_of_primative_mode : std::views::iota(0u, primative_projection_coefficient.size())) + for (auto i_of_direction : std::array{0, 1, 2}) + { + if (i_of_primative_mode != 0 && i_of_primative_mode != 11) continue; + log.debug("primitive mode: {} {} {:.3f}"_f + (i_of_primative_mode, i_of_direction, fmt::join(primative_projection_coefficient[i_of_primative_mode].col(i_of_direction).array().abs(), ", "))); + log.debug("primitive mode: {} {} {:.3f}"_f + (i_of_primative_mode, i_of_direction, fmt::join(primative_projection_coefficient[i_of_primative_mode].col(i_of_direction).array().arg(), ", "))); + } + log.debug("primitive 0 {}"_f(primative_projection_coefficient[0].col(2).transpose().conjugate() * primative_projection_coefficient[0].col(2))); + log.debug("primitive 11 {}"_f(primative_projection_coefficient[11].col(2).transpose().conjugate() * primative_projection_coefficient[11].col(2))); + log.debug("0 x 11 {}"_f(primative_projection_coefficient[11].col(2).transpose().conjugate() * primative_projection_coefficient[0].col(2))); + 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; + log.debug("super 2 x primitive 0 {}"_f(super_projection_coefficient[2].col(2).transpose().conjugate() * primative_projection_coefficient[0].col(2))); + log.debug("super 2 x primitive 11 {}"_f(super_projection_coefficient[2].col(2).transpose().conjugate() * primative_projection_coefficient[11].col(2))); + log.debug("super 2 x super 2 {}"_f(super_projection_coefficient[2].col(2).transpose().conjugate() * super_projection_coefficient[2].col(2))); + for (auto i_of_direction : std::array{0, 1, 2}) + { + log.debug("super mode: 2 {} {:.3f}"_f + (i_of_direction, fmt::join(super_projection_coefficient[2].col(i_of_direction).array().abs(), ", "))); + log.debug("super mode: 2 {} {:.3f}"_f + (i_of_direction, fmt::join(super_projection_coefficient[2].col(i_of_direction).array().arg(), ", "))); + } + // 将它们互相投影，得到结果 + 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)) + for (auto i_of_direction : std::array{0, 1, 2}) + output.Coefficient[i_of_super_mode][i_of_primative_mode] += + (super_mode.col(i_of_direction).transpose().conjugate() * primative_mode.col(i_of_direction)).norm(); + 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); +} diff --git a/src/unfold.cpp b/src/unfold.cpp index f0767df..d8ffcbd 100644 --- a/src/unfold.cpp +++ b/src/unfold.cpp @@ -7,9 +7,8 @@ void ufo::unfold(std::string config_file) { // 反折叠的原理: 将超胞中的原子运动状态, 投影到一组平面波构成的基矢中. - // 每一个平面波的波矢由两部分相加得到: 一部分是单胞倒格子的整数倍, 所取的个数有一定任意性, 论文中建议取大约单胞中原子个数那么多个; - // 对于没有缺陷的情况, 取一个应该就足够了. - // 这些平面波以原胞为周期。 + // 每一个平面波的波矢由两部分相加得到: 一部分是单胞倒格子的整数倍, 所取的个数有一定任意性, 论文中建议取大约单胞中原子个数那么多个， + // 但我觉得可以多取一些；这些平面波以原胞为周期，它用于尽可能描述 // 另一部分是超胞倒格子的整数倍, 取 n 个, n 为超胞对应的单胞的倍数, 其实也就是倒空间中单胞对应倒格子中超胞的格点. // 只要第一部分取得足够多, 那么单胞中原子的状态就可以完全被这些平面波描述. // 将超胞中原子的运动状态投影到这些基矢上, 计算出投影的系数, 就可以将超胞的原子运动状态分解到单胞中的多个 q 点上. @@ -143,16 +142,15 @@ void ufo::unfold(std::string config_file) for (auto [xyz_of_basis, i_of_basis] : biu::sequence(primative_cell_basis_number)) { - // 计算 q 点的坐标, 单位为单胞的倒格矢 - auto diff_of_sub_qpoint_by_reciprocal_primative_cell = xyz_of_basis.cast() + // 计算波矢, 单位为单胞的倒格矢 + auto wavevector_by_reciprocal_primative_cell = xyz_of_basis.cast() + super_cell_multiplier.cast().cwiseInverse().asDiagonal() * diff_of_sub_qpoint_by_reciprocal_modified_super_cell.cast(); // 将单位转换为埃^-1 - auto diff_of_sub_qpoint = (diff_of_sub_qpoint_by_reciprocal_primative_cell.transpose() - * (primative_cell.transpose().inverse())).transpose(); + auto wavevector = primative_cell.inverse() * wavevector_by_reciprocal_primative_cell; // 计算基矢 basis[i_of_sub_qpoint][i_of_basis] - = (2i * std::numbers::pi_v * (atom_position * diff_of_sub_qpoint)).array().exp(); + = (2i * std::numbers::pi_v * (atom_position * wavevector)).array().exp(); } } return basis; @@ -204,6 +202,46 @@ void ufo::unfold(std::string config_file) (config.SuperPhonopy, config.SuperQpoint, output.Super)); output.Super.CellDeformation = config.SuperCellDeformation; output.Super.CellMultiplier = config.SuperCellMultiplier; + // 填充 SubQpoint + for (auto i_of_super_qpoint : std::views::iota(0u, output.Super.Qpoint.size())) + for + ( + auto [diff_of_sub_qpoint_by_reciprocal_modified_super_cell, i_of_sub_qpoint] + : biu::sequence(config.SuperCellMultiplier) + ) + { + /* + 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().cwiseInverse().asDiagonal() + * ( + diff_of_sub_qpoint_by_reciprocal_modified_super_cell.cast() + + 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()); + } log.info("Done."); log.info("Constructing basis..."); diff --git a/test/project-to-mode/.gitignore b/test/project-to-mode/.gitignore new file mode 100644 index 0000000..fce0015 --- /dev/null +++ b/test/project-to-mode/.gitignore @@ -0,0 +1 @@ +output.yaml \ No newline at end of file diff --git a/test/project-to-mode/config.yaml b/test/project-to-mode/config.yaml new file mode 100644 index 0000000..f7842e1 --- /dev/null +++ b/test/project-to-mode/config.yaml @@ -0,0 +1,6 @@ +PrimativeCellBasisNumber: [ 16, 16, 16 ] +SuperQpointIndex: 0 +SubQpointIndex: 0 +PrimativeQpointIndex: 0 +InputDataFile: data.zpp +OutputFile: output.yaml diff --git a/test/project-to-mode/data.zpp b/test/project-to-mode/data.zpp new file mode 100644 index 0000000..65ce6cd --- /dev/null +++ b/test/project-to-mode/data.zpp @@ -0,0 +1,3 @@ +version https://git-lfs.github.com/spec/v1 +oid sha256:b1d66dad782f3070af08e8eb1357aa81b4d37404a3ac189df8d92cf6fba8587d +size 2921511 diff --git a/test/zpp-to-yaml/data.zpp b/test/zpp-to-yaml/data.zpp index bf725e6..65ce6cd 100644 --- a/test/zpp-to-yaml/data.zpp +++ b/test/zpp-to-yaml/data.zpp @@ -1,3 +1,3 @@ version https://git-lfs.github.com/spec/v1 -oid sha256:136f41ca83e6829a9bb3c58025bb80cde7f7fdb0cef34f88f352cf6f36f972de -size 2921063 +oid sha256:b1d66dad782f3070af08e8eb1357aa81b4d37404a3ac189df8d92cf6fba8587d +size 2921511