From d8f2fdcf7533376d09fc273f4b4c9460f25f654f Mon Sep 17 00:00:00 2001 From: chn Date: Wed, 20 Sep 2023 16:37:42 +0800 Subject: [PATCH] =?UTF-8?q?add=20fold=20man=20=E5=A2=9E=E5=8A=A0=E8=BE=93?= =?UTF-8?q?=E5=87=BA?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit --- fold.cpp | 98 ++++++++++++++++++++++++++++++++++++++++++++++++++++++++ main.cpp | 14 +++++++- 2 files changed, 111 insertions(+), 1 deletion(-) create mode 100644 fold.cpp diff --git a/fold.cpp b/fold.cpp new file mode 100644 index 0000000..c39db2a --- /dev/null +++ b/fold.cpp @@ -0,0 +1,98 @@ +# include +# include +# include +# include +# include +# include +# include +# include +# include +# include +# include +# include + +using namespace std::literals; + +// 计算单胞中的 q 点在超胞中的对应 + +struct Input +{ + // 单胞的三个格矢，每行表示一个格矢的坐标，单位为埃 + Eigen::Matrix3d PrimativeCell; + // 单胞到超胞的格矢转换时用到的矩阵 + Eigen::Vector SuperCellMultiplier; + Eigen::Matrix SuperCellDeformation; + + // Q 点的坐标，单位为单胞的倒格矢 + std::vector QPointData; +}; +template<> struct YAML::convert { static bool decode(const Node& node, Input& input); }; + +struct Output +{ + // Q 点的坐标，单位为超胞的倒格矢 + std::vector QPointData; +}; + +int main(int argc, char** argv) +{ + if (argc != 3) + throw std::runtime_error("Usage: " + std::string(argv[0]) + " input.yaml output.yaml"); + + auto input = YAML::LoadFile(argv[1]).as(); + Output output; + + for (auto qpoint_by_reciprocal_primative_cell : input.QPointData) + { + // 计算出这个 q 点的绝对坐标, 再计算出它相对于超胞倒格子的相对坐标. 将这个结果取小数部分, 就得到了 meta qpoint 的坐标(相对于超胞倒格子) + std::cout << "PrimativeCell:\n" << input.PrimativeCell << "\n"; + auto reciprocal_primative_cell = input.PrimativeCell.inverse().transpose(); + std::cout << "reciprocal_primative_cell:\n" << reciprocal_primative_cell << "\n"; + auto qpoint = (qpoint_by_reciprocal_primative_cell.transpose() * reciprocal_primative_cell).transpose(); + std::cout << "qpoint:\n" << qpoint << "\n"; + auto reciprocal_super_cell = (input.SuperCellDeformation * input.SuperCellMultiplier.cast().asDiagonal() + * input.PrimativeCell).inverse().transpose(); + std::cout << "reciprocal_super_cell:\n" << reciprocal_super_cell << "\n"; + auto qpoint_by_reciprocal_super_cell = (qpoint.transpose() * reciprocal_super_cell.inverse()).transpose(); + std::cout << "qpoint_by_reciprocal_super_cell:\n" << qpoint_by_reciprocal_super_cell << "\n"; + auto meta_qpoint_by_reciprocal_super_cell = [&] + { auto _ = qpoint_by_reciprocal_super_cell.array(); return _ - _.floor(); }(); + std::cout << "meta_qpoint_by_reciprocal_super_cell:\n" << meta_qpoint_by_reciprocal_super_cell << "\n"; + output.QPointData.push_back(meta_qpoint_by_reciprocal_super_cell); + } + + std::ofstream(argv[2]) << [&] + { + std::stringstream print; + print << "QPointData:\n"; + for (auto& qpoint: output.QPointData) + print << fmt::format(" - [ {:.3f}, {:.3f}, {:.3f} ]\n", qpoint(0), qpoint(1), qpoint(2)); + return print.str(); + }(); +} + +bool YAML::convert::decode(const Node& node, Input& input) +{ + for (unsigned i = 0; i < 3; i++) + for (unsigned j = 0; j < 3; j++) + input.PrimativeCell(i, j) = node["lattice"][i][j].as(); + + input.SuperCellMultiplier.setZero(); + for (unsigned i = 0; i < 3; i++) + input.SuperCellMultiplier(i) = node["SuperCellMultiplier"][i].as(); + + for (unsigned i = 0; i < 3; i++) + for (unsigned j = 0; j < 3; j++) + input.SuperCellDeformation(i, j) = node["SuperCellDeformation"][i][j].as(); + + auto points = node["points"].as>>(); + + input.QPointData.resize(points.size()); + for (unsigned i = 0; i < points.size(); i++) + { + for (unsigned j = 0; j < 3; j++) + input.QPointData[i](j) = points[i][j]; + } + + return true; +} diff --git a/main.cpp b/main.cpp index 642922d..89eabe8 100644 --- a/main.cpp +++ b/main.cpp @@ -10,6 +10,7 @@ # include # include # include +# include using namespace std::literals; @@ -70,7 +71,7 @@ struct Output // 关于各个 Q 点的数据 struct QPointDataType_ { - // Q 点的坐标，单位为超胞的倒格矢 + // Q 点的坐标，单位为单胞的倒格矢 Eigen::Vector3d QPoint; // 来源于哪个 Q 点, 单位为超胞的倒格矢 @@ -110,7 +111,9 @@ int main(int argc, const char** argv) if (argc != 3) throw std::runtime_error("Usage: " + std::string(argv[0]) + " input.yaml output.yaml"); + std::cerr << "Reading input file..." << std::flush; auto input = YAML::LoadFile(argv[1]).as(); + std::cerr << "Done." << std::endl; // 反折叠的原理: 将超胞中的原子运动状态, 投影到一组平面波构成的基矢中. // 每一个平面波的波矢由两部分相加得到: 一部分是单胞倒格子的整数倍, 所取的个数有一定任意性, 论文中建议取大约单胞中原子个数那么多个; @@ -125,6 +128,7 @@ int main(int argc, const char** argv) // 由于基只与这个相对位置有关（也就是说，不同 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 即表示这个相对位置，单位为超胞的倒格矢 @@ -146,12 +150,15 @@ int main(int argc, const char** argv) = (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()); for (unsigned i_of_qpoint = 0; i_of_qpoint < input.QPointData.size(); i_of_qpoint++) { + std::cerr << fmt::format("\rCalculating projection coefficient for qpoint {}/{}...", + i_of_qpoint, input.QPointData.size()) << std::flush; projection_coefficient[i_of_qpoint].resize(input.QPointData[i_of_qpoint].ModeData.size()); for (unsigned i_of_mode = 0; i_of_mode < input.QPointData[i_of_qpoint].ModeData.size(); i_of_mode++) { @@ -173,8 +180,10 @@ int main(int argc, const char** argv) __ /= sum; } } + 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++) for (auto [xyz_of_diff_of_sub_qpoint_by_reciprocal_super_cell, i_of_sub_qpoint] @@ -235,9 +244,11 @@ int main(int argc, const char** argv) __.Weight = projection_coefficient[i_of_qpoint][i_of_mode][i_of_sub_qpoint]; } } + std::cerr << "Done." << std::endl; // std::ofstream(argv[2]) << YAML::Node(output); // YAML 输出得太丑了，我来自己写 + std::cerr << "Writing output file..." << std::flush; std::ofstream(argv[2]) << [&] { std::stringstream print; @@ -256,6 +267,7 @@ int main(int argc, const char** argv) } return print.str(); }(); + std::cerr << "Done." << std::endl; } // 从文件中读取输入, 文件中应当包含: (大多数据可以直接从 phonopy 的输出中复制)