diff --git a/.gitignore b/.gitignore index 72be43c..3e739e9 100644 --- a/.gitignore +++ b/.gitignore @@ -1,9 +1,10 @@ .vscode .direnv -test/* main plot out.yaml .ccls-cache .cache build/ +test/* +!test/14.2.2.yaml diff --git a/flake.nix b/flake.nix index 46d76b3..2ab5b6e 100644 --- a/flake.nix +++ b/flake.nix @@ -4,12 +4,12 @@ outputs = inputs: let - pkgs = inputs.nixpkgs.legacyPackages.x86_64-linux; - # pkgs = import inputs.nixpkgs - # { - # localSystem = { system = "x86_64-linux"; gcc = { arch = "alderlake"; tune = "alderlake"; }; }; - # config.allowUnfree = true; - # }; + # pkgs = inputs.nixpkgs.legacyPackages.x86_64-linux; + pkgs = import inputs.nixpkgs + { + localSystem = { system = "x86_64-linux"; gcc = { arch = "alderlake"; tune = "alderlake"; }; }; + config.allowUnfree = true; + }; localPackages = import "${inputs.nixos}/local/pkgs" { inherit pkgs; inherit (inputs.nixpkgs) lib; }; in { diff --git a/include/ufo/ufo.impl.hpp b/include/ufo/ufo.impl.hpp index 77ab49c..99cf124 100644 --- a/include/ufo/ufo.impl.hpp +++ b/include/ufo/ufo.impl.hpp @@ -7,10 +7,10 @@ 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(); + 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(); @@ -26,28 +26,61 @@ inline Input::Input(std::string filename) 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(); + 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 format: {}, should be \"yaml\".", AtomPositionInputFile.Format)); - - QPointDataInputFile.FileName = node["QPointDataInputFile"]["FileName"].as(); - QPointDataInputFile.Format = node["QPointDataInputFile"]["Format"].as(); + ("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 format: {}, should be \"yaml\" or \"hdf5\".", QPointDataInputFile.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++) { - auto& _ = QPointDataOutputFile.emplace_back(); - _.FileName = value[i]["FileName"].as(); - _.Format = value[i]["Format"].as(); - if (!std::set{"yaml", "yaml-human-readable", "zpp"}.contains(_.Format)) + 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 format: {}, should be \"yaml\", \"yaml-human-readable\" or \"zpp\".", _.Format)); + ( + "Unknown QPointDataOutputFile[{}].Format: {}, should be \"yaml\", \"yaml-human-readable\" or \"zpp\".", + i, QPointDataOutputFile[i].Format + )); } + } } if (AtomPositionInputFile.Format == "yaml") diff --git a/src/plot.cpp b/src/plot.cpp index 6c2d4af..e04e117 100644 --- a/src/plot.cpp +++ b/src/plot.cpp @@ -1,32 +1,5 @@ # include -class FloatVector -{ - protected: - std::vector Data_; - double LowerBound_, UpperBound_, Step_; - public: - FloatVector(double LowerBound, double UpperBound, double Step) - : LowerBound_(LowerBound), UpperBound_(UpperBound), Step_(Step), Data_((UpperBound - LowerBound) / Step + 1) {} - FloatVector(const FloatVector&) = default; - FloatVector(FloatVector&&) = default; - FloatVector& operator=(const FloatVector&) = default; - FloatVector& operator=(FloatVector&&) = default; - double& operator[](double i) { return Data_[static_cast((i - LowerBound_) / Step_)]; } - double operator[](double i) const { return Data_[static_cast((i - LowerBound_) / Step_)]; } - double lower_bound() const { return LowerBound_; } - double upper_bound() const { return UpperBound_; } - double step() const { return Step_; } - int size() const { return Data_.size(); } - std::map to_map() const - { - std::map result; - for (int i = 0; i < Data_.size(); i++) - result[LowerBound_ + i * Step_] = Data_[i]; - return result; - } -}; - // 要被用来画图的路径 std::vector Qpoints = { @@ -41,6 +14,13 @@ std::vector Qpoints = }; double Threshold = 0.001; +struct Point +{ + Eigen::Vector3d QPoint; + Eigen::VectorXd Frequency, Weight; + double Distance; +}; + int main(int argc, char** argv) { if (argc != 2) @@ -48,17 +28,9 @@ int main(int argc, char** argv) Output output(argv[1]); - // - - struct Point - { - Eigen::Vector3d QPoint; - FloatVector Weight; - double Distance; - }; std::vector Points; - double current_distance = 0; + double total_distance = 0; // 对于每一条路径进行搜索 for (unsigned i = 0; i < Qpoints.size() - 1; i++) { @@ -81,10 +53,16 @@ int main(int argc, char** argv) // 如果这个点在终点处, 且这条路径不是最后一条, 则不加入 if (distance2 > Threshold || i == Qpoints.size() - 2) { - auto& _ = point_of_this_path.emplace_back - (qpoint.QPoint, FloatVector{-5, 35, 0.1}, distance1); - for (auto& mode : qpoint.ModeData) - _.Weight[mode.Frequency] += mode.Weight; + auto& _ = point_of_this_path.emplace_back(); + _.QPoint = qpoint.QPoint; + _.Distance = distance1; + _.Frequency.resize(qpoint.ModeData.size()); + _.Weight.resize(qpoint.ModeData.size()); + for (unsigned j = 0; j < qpoint.ModeData.size(); j++) + { + _.Frequency(j) = qpoint.ModeData[j].Frequency; + _.Weight(j) = qpoint.ModeData[j].Weight; + } } } } @@ -101,7 +79,36 @@ int main(int argc, char** argv) point_of_this_path.erase(point_of_this_path.begin() + j); // 将结果加入 for (auto& point : point_of_this_path) - Points.emplace_back(point.QPoint, point.Weight, point.Distance + current_distance); - current_distance += (Qpoints[i + 1] - Qpoints[i]).norm(); + Points.emplace_back(point.QPoint, point.Frequency, point.Weight, point.Distance + total_distance); + total_distance += (Qpoints[i + 1] - Qpoints[i]).norm(); + } + + // 对结果插值 + std::vector interpolated_points; + for (unsigned i = 0; i < 1024; i++) + { + auto current_distance = i * total_distance / 1024; + auto& _ = interpolated_points.emplace_back(); + _.Distance = current_distance; + // 如果是开头或者结尾, 直接赋值, 否则插值 + if (current_distance < Points.front().Distance) + { + _.Frequency = Points.front().Frequency; + _.Weight = Points.front().Weight; + } + else if (current_distance > Points.back().Distance) + { + _.Frequency = Points.back().Frequency; + _.Weight = Points.back().Weight; + } + else + { + auto it = std::lower_bound(Points.begin(), Points.end(), current_distance, + [](const Point& a, double b) { return a.Distance < b; }); + _.Frequency = (it->Frequency * (it->Distance - current_distance) + + (it - 1)->Frequency * (current_distance - (it - 1)->Distance)) / (it->Distance - (it - 1)->Distance); + _.Weight = (it->Weight * (it->Distance - current_distance) + + (it - 1)->Weight * (current_distance - (it - 1)->Distance)) / (it->Distance - (it - 1)->Distance); + } } } diff --git a/src/ufo.cpp b/src/ufo.cpp index 286309b..37bde3a 100644 --- a/src/ufo.cpp +++ b/src/ufo.cpp @@ -90,7 +90,7 @@ int main(int argc, const char** argv) { // 当 SuperCellDeformation 不是单位矩阵时, input.QPointData[i_of_qpoint].QPoint 不一定在 reciprocal_primative_cell 中 // 需要首先将 q 点平移数个周期, 进入不包含 SuperCellDeformation 的超胞的倒格子中 - auto qpoint_by_reciprocal_super_cell_in_modified_reciprocal_super_cell + 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; @@ -139,16 +139,12 @@ int main(int argc, const char** argv) } current_qpoint = min_score_qpoint; } - return current_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(); - // 这一步推导过程在计算 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 diff --git a/test/14.2.2.yaml b/test/14.2.2.yaml new file mode 100644 index 0000000..95b59b6 --- /dev/null +++ b/test/14.2.2.yaml @@ -0,0 +1,12 @@ +PrimativeCell: +- [ 1.548010265167714, -2.681232429908652, 0.000000000000000 ] # a +- [ 1.548010265167714, 2.681232429908652, 0.000000000000000 ] # b +- [ 0.000000000000000, 0.000000000000000, 5.061224103556595 ] # c +SuperCellMultiplier: [ 3, 1, 1 ] +PrimativeCellBasisNumber: [ 8, 8, 8 ] +AtomPositionInputFile: { FileName: "/home/chn/Documents/lammps-SiC/14/14.2/14.2.2/14.2.2.4/band.yaml", Format: "yaml" } +QPointDataInputFile: { FileName: "/home/chn/Documents/lammps-SiC/14/14.2/14.2.2/14.2.2.4/band.yaml", Format: "yaml" } +QPointDataOutputFile: +- { FileName: "test/14.2.2.result.yaml", Format: "yaml-human-readable" } +- { FileName: "test/14.2.2.result.human-readable.yaml", Format: "yaml-human-readable" } +- { FileName: "test/14.2.2.result.zpp", Format: "zpp" }