mirror of
https://github.com/CHN-beta/ufo.git
synced 2024-10-22 19:58:44 +08:00
plot 完成模式的查找
This commit is contained in:
parent
312dbd0529
commit
e042ae4c15
@ -27,6 +27,11 @@ target_include_directories(ufo PRIVATE ${PROJECT_SOURCE_DIR}/include ${ZPP_BITS_
|
||||
target_link_libraries(ufo PRIVATE
|
||||
yaml-cpp Eigen3::Eigen fmt::fmt concurrencpp::concurrencpp HighFive_HighFive TBB::tbb Matplot++::matplot)
|
||||
|
||||
add_executable(plot src/plot.cpp)
|
||||
target_include_directories(plot PRIVATE ${PROJECT_SOURCE_DIR}/include ${ZPP_BITS_INCLUDE_DIR})
|
||||
target_link_libraries(plot PRIVATE
|
||||
yaml-cpp Eigen3::Eigen fmt::fmt concurrencpp::concurrencpp HighFive_HighFive TBB::tbb Matplot++::matplot)
|
||||
|
||||
get_property(ImportedTargets DIRECTORY "${CMAKE_SOURCE_DIR}" PROPERTY IMPORTED_TARGETS)
|
||||
message("Imported targets: ${ImportedTargets}")
|
||||
message("List of compile features: ${CMAKE_CXX_COMPILE_FEATURES}")
|
||||
|
||||
@ -29,7 +29,9 @@ HIGHFIVE_REGISTER_TYPE(PhonopyComplex, create_compound_complex)
|
||||
namespace Eigen
|
||||
{
|
||||
constexpr inline auto serialize(auto & archive, Eigen::Matrix3d& matrix)
|
||||
{ return archive(std::span<double, 9>(matrix.data(), 9)); }
|
||||
{ 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())); }
|
||||
}
|
||||
|
||||
// 在相位中, 约定为使用 $\exp (2 \pi i \vec{q} \cdot \vec{r})$ 来表示原子的运动状态
|
||||
@ -56,6 +58,19 @@ struct Input
|
||||
std::optional<Eigen::Matrix<double, 3, 3>> SuperCellDeformation;
|
||||
// 在单胞内取几个平面波的基矢
|
||||
Eigen::Vector<unsigned, 3> PrimativeCellBasisNumber;
|
||||
|
||||
struct InputOutputFile_
|
||||
{
|
||||
std::string FileName;
|
||||
std::string Format;
|
||||
std::map<std::string, std::any> ExtraParameters;
|
||||
};
|
||||
|
||||
// 从哪个文件读入 AtomPosition, 以及这个文件的格式, 格式可选值包括 "yaml"
|
||||
InputOutputFile_ AtomPositionInputFile;
|
||||
// 从哪个文件读入 QPointData, 以及这个文件的格式, 格式可选值包括 "yaml" 和 "hdf5"
|
||||
InputOutputFile_ QPointDataInputFile;
|
||||
|
||||
// 超胞中原子的坐标，每行表示一个原子的坐标，单位为埃
|
||||
Eigen::MatrixX3d AtomPosition;
|
||||
// 关于各个 Q 点的数据
|
||||
@ -79,18 +94,6 @@ struct Input
|
||||
};
|
||||
std::vector<QPointDataType_> QPointData;
|
||||
|
||||
struct InputOutputFile_
|
||||
{
|
||||
std::string FileName;
|
||||
std::string Format;
|
||||
std::map<std::string, std::any> ExtraParameters;
|
||||
};
|
||||
|
||||
// 从哪个文件读入 AtomPosition, 以及这个文件的格式, 格式可选值包括 "yaml"
|
||||
InputOutputFile_ AtomPositionInputFile;
|
||||
// 从哪个文件读入 QPointData, 以及这个文件的格式, 格式可选值包括 "yaml" 和 "hdf5"
|
||||
InputOutputFile_ QPointDataInputFile;
|
||||
|
||||
// 输出到哪些文件, 以及使用怎样的格式, 格式可选值包括:
|
||||
// yaml: 使用 yaml 格式输出
|
||||
// yaml-human-readable: 使用 yaml 格式输出, 但是输出的结果更适合人类阅读, 包括合并相近的模式, 去除权重过小的模式, 限制输出的小数位数.
|
||||
@ -138,6 +141,8 @@ struct Output
|
||||
void write(std::string filename, std::string format, unsigned percision = 10);
|
||||
Output() = default;
|
||||
Output(std::string filename);
|
||||
|
||||
using serialize = zpp::bits::members<1>;
|
||||
};
|
||||
|
||||
concurrencpp::generator<std::pair<Eigen::Vector<unsigned, 3>, unsigned>>
|
||||
|
||||
@ -193,11 +193,20 @@ inline void Output::write(std::string filename, std::string format, unsigned per
|
||||
else if (format == "zpp")
|
||||
{
|
||||
auto [data, out] = zpp::bits::data_out();
|
||||
out(*this).or_throw();
|
||||
std::basic_ofstream<std::byte>(filename, std::ios::binary)
|
||||
<< std::basic_string_view{data.data(), data.size()};
|
||||
}
|
||||
}
|
||||
|
||||
inline Output::Output(std::string filename)
|
||||
{
|
||||
auto [data, in] = zpp::bits::data_in();
|
||||
auto input = std::basic_ifstream<std::byte>(filename, std::ios::binary);
|
||||
data.assign(std::istreambuf_iterator<std::byte>(input), {});
|
||||
in(*this).or_throw();
|
||||
}
|
||||
|
||||
inline concurrencpp::generator<std::pair<Eigen::Vector<unsigned, 3>, unsigned>>
|
||||
triplet_sequence(Eigen::Vector<unsigned, 3> range)
|
||||
{
|
||||
|
||||
164
src/plot.cpp
164
src/plot.cpp
@ -1,93 +1,107 @@
|
||||
# include <ufo/ufo.impl.hpp>
|
||||
|
||||
// 一个临时的程序, 用于将数据导出画图
|
||||
|
||||
struct Output
|
||||
class FloatVector
|
||||
{
|
||||
struct QPointDataType_
|
||||
{
|
||||
std::vector<double> QPoint;
|
||||
struct ModeDataType_
|
||||
protected:
|
||||
std::vector<double> 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<int>((i - LowerBound_) / Step_)]; }
|
||||
double operator[](double i) const { return Data_[static_cast<int>((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<double, double> to_map() const
|
||||
{
|
||||
double Frequency;
|
||||
double Weight;
|
||||
};
|
||||
std::vector<ModeDataType_> ModeData;
|
||||
};
|
||||
std::vector<QPointDataType_> QPointData;
|
||||
std::map<double, double> result;
|
||||
for (int i = 0; i < Data_.size(); i++)
|
||||
result[LowerBound_ + i * Step_] = Data_[i];
|
||||
return result;
|
||||
}
|
||||
};
|
||||
|
||||
// BAND = 0.000000 0.000000 0.000000 0.500000 0.000000 0.000000 0.333333 0.333333 0.000000 0.000000 0.000000 0.000000 0 0 0.5 0.5 0 0.5 0.333333 0.333333 0.5 0 0 0.5
|
||||
// BAND_LABELS = $\Gamma$ M K $\Gamma$ A L H A
|
||||
|
||||
// Gamma 0 0 0
|
||||
// M 1/2 0 0
|
||||
// K 1/3 1/3 0
|
||||
// A 0 0 1/2
|
||||
// L 1/2 0 1/2
|
||||
// H 1/3 1/3 1/2
|
||||
|
||||
std::vector<std::vector<double>> Qpoints =
|
||||
// 要被用来画图的路径
|
||||
std::vector<Eigen::Vector3d> Qpoints =
|
||||
{
|
||||
{0, 0, 0},
|
||||
{0.025, 0, 0},
|
||||
{0.05, 0, 0},
|
||||
{0.075, 0, 0},
|
||||
{0.1, 0, 0},
|
||||
{0.125, 0, 0},
|
||||
{0.15, 0, 0},
|
||||
{0.175, 0, 0},
|
||||
{0.2, 0, 0},
|
||||
{0.225, 0, 0},
|
||||
{0.25, 0, 0},
|
||||
{0.275, 0, 0},
|
||||
{0.3, 0, 0},
|
||||
{0.325, 0, 0},
|
||||
{0.35, 0, 0},
|
||||
{0.375, 0, 0},
|
||||
{0.4, 0, 0},
|
||||
{0.425, 0, 0},
|
||||
{0.45, 0, 0},
|
||||
{0.475, 0, 0},
|
||||
{0.5, 0, 0}
|
||||
{ 0, 0, 0 },
|
||||
{ 0.5, 0, 0 },
|
||||
{ 1. / 3, 1. / 3, 0 },
|
||||
{ 0, 0, 0 },
|
||||
{ 0, 0, 0.5 },
|
||||
{ 0.5, 0, 0.5 },
|
||||
{ 1. / 3, 1. / 3, 0.5 },
|
||||
{ 0, 0, 0.5 }
|
||||
};
|
||||
double Threshold = 0.001;
|
||||
|
||||
int main(int argc, char** argv)
|
||||
{
|
||||
YAML::Node root = YAML::LoadFile(argv[1]);
|
||||
Output output;
|
||||
output.QPointData.resize(root["QPointData"].size());
|
||||
for (int i = 0; i < root["QPointData"].size(); i++)
|
||||
{
|
||||
output.QPointData[i].QPoint = root["QPointData"][i]["QPoint"].as<std::vector<double>>();
|
||||
output.QPointData[i].ModeData.resize(root["QPointData"][i]["ModeData"].size());
|
||||
for (int j = 0; j < root["QPointData"][i]["ModeData"].size(); j++)
|
||||
{
|
||||
output.QPointData[i].ModeData[j].Frequency
|
||||
= root["QPointData"][i]["ModeData"][j]["Frequency"].as<double>();
|
||||
output.QPointData[i].ModeData[j].Weight
|
||||
= root["QPointData"][i]["ModeData"][j]["Weight"].as<double>();
|
||||
}
|
||||
}
|
||||
if (argc != 2)
|
||||
throw std::runtime_error("Usage: plot output.zpp");
|
||||
|
||||
// 外层表示 q 点坐标, 内层表示频率
|
||||
// 频率取 -5 到 30 THz, 每 0.1 THz 一个点
|
||||
std::vector<std::vector<double>> data(21);
|
||||
for (int i = 0; i < 21; i++)
|
||||
Output output(argv[1]);
|
||||
|
||||
//
|
||||
|
||||
struct Point
|
||||
{
|
||||
data[i].resize(351);
|
||||
Eigen::Vector3d QPoint;
|
||||
FloatVector Weight;
|
||||
double Distance;
|
||||
};
|
||||
std::vector<Point> Points;
|
||||
|
||||
double current_distance = 0;
|
||||
// 对于每一条路径进行搜索
|
||||
for (unsigned i = 0; i < Qpoints.size() - 1; i++)
|
||||
{
|
||||
std::vector<Point> point_of_this_path;
|
||||
// 对于 output 中的每一个点, 检查这个点是否在路径上. 如果在, 把它加入到 point_of_this_path 中
|
||||
for (auto& qpoint : output.QPointData)
|
||||
if (std::abs(qpoint.QPoint[0] - Qpoints[i][0]) < 1e-3 &&
|
||||
std::abs(qpoint.QPoint[1] - Qpoints[i][1]) < 1e-3 &&
|
||||
std::abs(qpoint.QPoint[2] - Qpoints[i][2]) < 1e-3)
|
||||
{
|
||||
// 计算三点围成的三角形的面积的两倍
|
||||
auto area = (Qpoints[i + 1] - Qpoints[i]).cross(qpoint.QPoint - Qpoints[i]).norm();
|
||||
// 计算这个点到前两个点所在直线的距离
|
||||
auto distance = area / (Qpoints[i + 1] - Qpoints[i]).norm();
|
||||
// 如果这个点到前两个点所在直线的距离小于阈值, 则认为这个点在路径上
|
||||
if (distance < Threshold)
|
||||
{
|
||||
for (auto& mode : qpoint.ModeData)
|
||||
data[i][static_cast<int>((mode.Frequency + 5) * 10)] += mode.Weight;
|
||||
break;
|
||||
// 计算这个点到前两个点的距离, 两个距离都应该小于两点之间的距离
|
||||
auto distance1 = (qpoint.QPoint - Qpoints[i]).norm();
|
||||
auto distance2 = (qpoint.QPoint - Qpoints[i + 1]).norm();
|
||||
auto distance3 = (Qpoints[i + 1] - Qpoints[i]).norm();
|
||||
if (distance1 < distance3 + Threshold && distance2 < distance3 + Threshold)
|
||||
// 如果这个点在终点处, 且这条路径不是最后一条, 则不加入
|
||||
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;
|
||||
}
|
||||
}
|
||||
}
|
||||
// 对筛选结果排序
|
||||
std::sort(point_of_this_path.begin(), point_of_this_path.end(),
|
||||
[](const Point& a, const Point& b) { return a.Distance < b.Distance; });
|
||||
// 去除非常接近的点
|
||||
for (unsigned j = 1; j < point_of_this_path.size(); j++)
|
||||
while
|
||||
(
|
||||
j < point_of_this_path.size()
|
||||
&& point_of_this_path[j].Distance - point_of_this_path[j - 1].Distance < Threshold
|
||||
)
|
||||
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();
|
||||
}
|
||||
|
||||
matplot::image(data, true);
|
||||
matplot::colorbar().limits({0, 0.5});
|
||||
matplot::show();
|
||||
}
|
||||
|
||||
11
test.cpp
11
test.cpp
@ -2,6 +2,7 @@
|
||||
# include <zpp_bits.h>
|
||||
# include <span>
|
||||
# include <iostream>
|
||||
# include <vector>
|
||||
|
||||
namespace Eigen
|
||||
{
|
||||
@ -18,9 +19,9 @@ namespace Eigen
|
||||
int main()
|
||||
{
|
||||
auto [data, in, out] = zpp::bits::data_in_out();
|
||||
Eigen::Matrix3d PrimativeCell = Eigen::Matrix3d::Identity();
|
||||
out(PrimativeCell);
|
||||
Eigen::Matrix3d PrimativeCell2;
|
||||
in(PrimativeCell2);
|
||||
std::cout << PrimativeCell2 << std::endl;
|
||||
std::vector<Eigen::Matrix3d> PrimativeCell = { Eigen::Matrix3d::Identity() };
|
||||
out(PrimativeCell).or_throw();
|
||||
std::vector<Eigen::Matrix3d> PrimativeCell2;
|
||||
in(PrimativeCell2).or_throw();
|
||||
std::cout << PrimativeCell2[0] << std::endl;
|
||||
}
|
||||
Loading…
Reference in New Issue
Block a user