增加更多测试数据

统一 QPoint 为 Qpoint
fold 不再需要输入 PrimativeCell
This commit is contained in:
陈浩南 2023-10-07 15:03:56 +08:00
parent e60054579e
commit f74258e8fc
11 changed files with 166 additions and 95 deletions

4
.gitignore vendored
View File

@ -12,3 +12,7 @@ test/*
!test/14.2.2.5.plot.yaml !test/14.2.2.5.plot.yaml
!test/14.2.4.4.unfold.yaml !test/14.2.4.4.unfold.yaml
!test/14.2.4.4.plot.yaml !test/14.2.4.4.plot.yaml
!test/14.2.5.5.unfold.yaml
!test/14.2.5.5.plot.yaml
!test/14.2.6.4.unfold.yaml
!test/14.2.6.4.plot.yaml

View File

@ -8,7 +8,6 @@ namespace ufo
public: public:
struct InputType struct InputType
{ {
Eigen::Matrix3d PrimativeCell;
Eigen::Vector<unsigned, 3> SuperCellMultiplier; Eigen::Vector<unsigned, 3> SuperCellMultiplier;
std::optional<Eigen::Matrix<double, 3, 3>> SuperCellDeformation; std::optional<Eigen::Matrix<double, 3, 3>> SuperCellDeformation;
std::vector<Eigen::Vector3d> Qpoints; std::vector<Eigen::Vector3d> Qpoints;
@ -18,7 +17,7 @@ namespace ufo
}; };
struct OutputType struct OutputType
{ {
std::vector<Eigen::Vector3d> QPoints; std::vector<Eigen::Vector3d> Qpoints;
void write(std::string filename) const; void write(std::string filename) const;
}; };
protected: protected:
@ -27,7 +26,12 @@ namespace ufo
public: public:
FoldSolver(std::string config_file); FoldSolver(std::string config_file);
FoldSolver& operator()() override; FoldSolver& operator()() override;
static Eigen::Vector3d fold(Eigen::Vector3d qpoint, Eigen::Vector<unsigned, 3> super_cell_multiplier, // return value: QpointInReciprocalSuperCellByReciprocalSuperCell
std::optional<Eigen::Matrix<double, 3, 3>> super_cell_deformation); static Eigen::Vector3d fold
(
Eigen::Vector3d qpoint_in_reciprocal_primitive_cell_by_reciprocal_primitive_cell,
Eigen::Vector<unsigned, 3> super_cell_multiplier,
std::optional<Eigen::Matrix<double, 3, 3>> super_cell_deformation
);
}; };
} }

View File

@ -36,17 +36,17 @@ namespace ufo
PlotSolver& operator()() override; PlotSolver& operator()() override;
// 根据 q 点路径, 搜索要使用的 q 点 // 根据 q 点路径, 搜索要使用的 q 点
static std::vector<std::reference_wrapper<const UnfoldSolver::OutputType::QPointDataType>> search_qpoints static std::vector<std::reference_wrapper<const UnfoldSolver::OutputType::QpointDataType>> search_qpoints
( (
const std::pair<Eigen::Vector3d, Eigen::Vector3d>& path, const std::pair<Eigen::Vector3d, Eigen::Vector3d>& path,
const decltype(InputType::SourceType::QPointData)& available_qpoints, const decltype(InputType::SourceType::QpointData)& available_qpoints,
double threshold, bool exclude_endpoint = false double threshold, bool exclude_endpoint = false
); );
// 根据搜索到的 q 点, 计算每个点的数值 // 根据搜索到的 q 点, 计算每个点的数值
static std::vector<std::vector<double>> calculate_values static std::vector<std::vector<double>> calculate_values
( (
const std::vector<std::pair<Eigen::Vector3d, Eigen::Vector3d>>& path, const std::vector<std::pair<Eigen::Vector3d, Eigen::Vector3d>>& path,
const std::vector<std::vector<std::reference_wrapper<const UnfoldSolver::OutputType::QPointDataType>>>& qpoints, const std::vector<std::vector<std::reference_wrapper<const UnfoldSolver::OutputType::QpointDataType>>>& qpoints,
const decltype(InputType::FigureConfigType::Resolution)& resolution, const decltype(InputType::FigureConfigType::Resolution)& resolution,
const decltype(InputType::FigureConfigType::Range)& range const decltype(InputType::FigureConfigType::Range)& range
); );

View File

@ -41,16 +41,16 @@ namespace ufo
// 从哪个文件读入 AtomPosition, 以及这个文件的格式, 格式可选值包括 "yaml" // 从哪个文件读入 AtomPosition, 以及这个文件的格式, 格式可选值包括 "yaml"
InputOutputFile AtomPositionInputFile; InputOutputFile AtomPositionInputFile;
// 从哪个文件读入 QPointData, 以及这个文件的格式, 格式可选值包括 "yaml" 和 "hdf5" // 从哪个文件读入 QpointData, 以及这个文件的格式, 格式可选值包括 "yaml" 和 "hdf5"
InputOutputFile QPointDataInputFile; InputOutputFile QpointDataInputFile;
// 超胞中原子的坐标,每行表示一个原子的坐标,单位为埃 // 超胞中原子的坐标,每行表示一个原子的坐标,单位为埃
Eigen::MatrixX3d AtomPosition; Eigen::MatrixX3d AtomPosition;
// 关于各个 Q 点的数据 // 关于各个 Q 点的数据
struct QPointDataType struct QpointDataType
{ {
// Q 点的坐标,单位为超胞的倒格矢 // Q 点的坐标,单位为超胞的倒格矢
Eigen::Vector3d QPoint; Eigen::Vector3d Qpoint;
// 关于这个 Q 点上各个模式的数据 // 关于这个 Q 点上各个模式的数据
struct ModeDataType struct ModeDataType
@ -65,14 +65,14 @@ namespace ufo
}; };
std::vector<ModeDataType> ModeData; std::vector<ModeDataType> ModeData;
}; };
std::vector<QPointDataType> QPointData; std::vector<QpointDataType> QpointData;
// 输出到哪些文件, 以及使用怎样的格式, 格式可选值包括: // 输出到哪些文件, 以及使用怎样的格式, 格式可选值包括:
// yaml: 使用 yaml 格式输出 // yaml: 使用 yaml 格式输出
// yaml-human-readable: 使用 yaml 格式输出, 但是输出的结果更适合人类阅读, // yaml-human-readable: 使用 yaml 格式输出, 但是输出的结果更适合人类阅读,
// 包括合并相近的模式, 去除权重过小的模式, 限制输出的小数位数. // 包括合并相近的模式, 去除权重过小的模式, 限制输出的小数位数.
// zpp: 使用 zpp-bits 序列化, 可以直接被 plot.cpp 读取 // zpp: 使用 zpp-bits 序列化, 可以直接被 plot.cpp 读取
std::vector<InputOutputFile> QPointDataOutputFile; std::vector<InputOutputFile> QpointDataOutputFile;
// 从文件中读取输入 (包括一个较小的配置文件, 和一个 hdf5 或者一个 yaml 文件), 文件中应当包含: // 从文件中读取输入 (包括一个较小的配置文件, 和一个 hdf5 或者一个 yaml 文件), 文件中应当包含:
// 单胞的格矢: PrimativeCell 单位为埃 直接从 phonopy 的输出中复制 // 单胞的格矢: PrimativeCell 单位为埃 直接从 phonopy 的输出中复制
@ -90,10 +90,10 @@ namespace ufo
struct OutputType struct OutputType
{ {
// 关于各个 Q 点的数据 // 关于各个 Q 点的数据
struct QPointDataType struct QpointDataType
{ {
// Q 点的坐标,单位为单胞的倒格矢 // Q 点的坐标,单位为单胞的倒格矢
Eigen::Vector3d QPoint; Eigen::Vector3d Qpoint;
// 来源于哪个 Q 点, 单位为超胞的倒格矢 // 来源于哪个 Q 点, 单位为超胞的倒格矢
Eigen::Vector3d Source; Eigen::Vector3d Source;
@ -109,9 +109,9 @@ namespace ufo
}; };
std::vector<ModeDataType> ModeData; std::vector<ModeDataType> ModeData;
}; };
std::vector<QPointDataType> QPointData; std::vector<QpointDataType> QpointData;
void write(decltype(InputType::QPointDataOutputFile) output_files) const; void write(decltype(InputType::QpointDataOutputFile) output_files) const;
void write(std::string filename, std::string format, unsigned percision = 10) const; void write(std::string filename, std::string format, unsigned percision = 10) const;
using serialize = zpp::bits::members<1>; using serialize = zpp::bits::members<1>;
@ -153,7 +153,7 @@ namespace ufo
( (
const BasisType& basis, const BasisType& basis,
const std::vector<std::reference_wrapper<const decltype const std::vector<std::reference_wrapper<const decltype
(InputType::QPointDataType::ModeDataType::AtomMovement)>>& mode_data, (InputType::QpointDataType::ModeDataType::AtomMovement)>>& mode_data,
std::atomic<unsigned>& number_of_finished_modes std::atomic<unsigned>& number_of_finished_modes
); );
@ -162,9 +162,9 @@ namespace ufo
const decltype(InputType::SuperCellMultiplier)& super_cell_multiplier, const decltype(InputType::SuperCellMultiplier)& super_cell_multiplier,
const decltype(InputType::SuperCellDeformation)& super_cell_deformation, const decltype(InputType::SuperCellDeformation)& super_cell_deformation,
const std::vector<std::reference_wrapper<const decltype const std::vector<std::reference_wrapper<const decltype
(InputType::QPointDataType::QPoint)>>& qpoint, (InputType::QpointDataType::Qpoint)>>& qpoint,
const std::vector<std::vector<std::reference_wrapper<const decltype const std::vector<std::vector<std::reference_wrapper<const decltype
(InputType::QPointDataType::ModeDataType::Frequency)>>>& frequency, (InputType::QpointDataType::ModeDataType::Frequency)>>>& frequency,
const ProjectionCoefficientType_& projection_coefficient const ProjectionCoefficientType_& projection_coefficient
); );
}; };

View File

@ -5,9 +5,6 @@ namespace ufo
FoldSolver::InputType::InputType(std::string config_file) FoldSolver::InputType::InputType(std::string config_file)
{ {
auto input = YAML::LoadFile(config_file); auto input = YAML::LoadFile(config_file);
for (unsigned i = 0; i < 3; i++)
for (unsigned j = 0; j < 3; j++)
PrimativeCell(i, j) = input["PrimativeCell"][i][j].as<double>();
for (unsigned i = 0; i < 3; i++) for (unsigned i = 0; i < 3; i++)
SuperCellMultiplier(i) = input["SuperCellMultiplier"][i].as<unsigned>(); SuperCellMultiplier(i) = input["SuperCellMultiplier"][i].as<unsigned>();
if (input["SuperCellDeformation"]) if (input["SuperCellDeformation"])
@ -27,8 +24,8 @@ namespace ufo
std::ofstream(filename) << [&] std::ofstream(filename) << [&]
{ {
std::stringstream print; std::stringstream print;
print << "QPoints:\n"; print << "Qpoints:\n";
for (auto& qpoint : QPoints) for (auto& qpoint : Qpoints)
print << fmt::format(" - [ {:.8f}, {:.8f}, {:.8f} ]\n", qpoint(0), qpoint(1), qpoint(2)); print << fmt::format(" - [ {:.8f}, {:.8f}, {:.8f} ]\n", qpoint(0), qpoint(1), qpoint(2));
return print.str(); return print.str();
}(); }();
@ -41,7 +38,7 @@ namespace ufo
{ {
Output_.emplace(); Output_.emplace();
for (auto& qpoint : Input_.Qpoints) for (auto& qpoint : Input_.Qpoints)
Output_->QPoints.push_back(fold Output_->Qpoints.push_back(fold
( (
qpoint, Input_.SuperCellMultiplier, qpoint, Input_.SuperCellMultiplier,
Input_.SuperCellDeformation Input_.SuperCellDeformation
@ -53,13 +50,15 @@ namespace ufo
Eigen::Vector3d FoldSolver::fold Eigen::Vector3d FoldSolver::fold
( (
Eigen::Vector3d qpoint, Eigen::Vector<unsigned, 3> super_cell_multiplier, Eigen::Vector3d qpoint_in_reciprocal_primitive_cell_by_reciprocal_primitive_cell,
Eigen::Vector<unsigned, 3> super_cell_multiplier,
std::optional<Eigen::Matrix<double, 3, 3>> super_cell_deformation std::optional<Eigen::Matrix<double, 3, 3>> super_cell_deformation
) )
{ {
// 首先需要将 q 点转移到超胞的倒格子中 // 首先需要将 q 点转移到 ModifiedSuperCell 的倒格子中
// 将 q 点坐标扩大, 然后取小数部分, 就可以了 // 将 q 点坐标扩大, 然后取小数部分, 就可以了
auto qpoint_by_reciprocal_modified_super_cell = super_cell_multiplier.cast<double>().asDiagonal() * qpoint; auto qpoint_by_reciprocal_modified_super_cell = super_cell_multiplier.cast<double>().asDiagonal()
* qpoint_in_reciprocal_primitive_cell_by_reciprocal_primitive_cell;
auto qpoint_in_reciprocal_modified_super_cell_by_reciprocal_modified_super_cell = auto qpoint_in_reciprocal_modified_super_cell_by_reciprocal_modified_super_cell =
(qpoint_by_reciprocal_modified_super_cell.array() - qpoint_by_reciprocal_modified_super_cell.array().floor()) (qpoint_by_reciprocal_modified_super_cell.array() - qpoint_by_reciprocal_modified_super_cell.array().floor())
.matrix(); .matrix();

View File

@ -58,7 +58,7 @@ namespace ufo
for (auto& figure : Input_.Figures) for (auto& figure : Input_.Figures)
{ {
// 外层表示不同的线段的端点,内层表示这个线段上的 q 点 // 外层表示不同的线段的端点,内层表示这个线段上的 q 点
std::vector<std::vector<std::reference_wrapper<const UnfoldSolver::OutputType::QPointDataType>>> qpoints; std::vector<std::vector<std::reference_wrapper<const UnfoldSolver::OutputType::QpointDataType>>> qpoints;
std::vector<std::pair<Eigen::Vector3d, Eigen::Vector3d>> lines; std::vector<std::pair<Eigen::Vector3d, Eigen::Vector3d>> lines;
for (auto& path : figure.Qpoints) for (auto& path : figure.Qpoints)
for (unsigned i = 0; i < path.size() - 1; i++) for (unsigned i = 0; i < path.size() - 1; i++)
@ -66,7 +66,7 @@ namespace ufo
lines.emplace_back(path[i], path[i + 1]); lines.emplace_back(path[i], path[i + 1]);
qpoints.push_back(search_qpoints qpoints.push_back(search_qpoints
( (
lines.back(), Input_.Source.QPointData, lines.back(), Input_.Source.QpointData,
0.001, i != path.size() - 2 0.001, i != path.size() - 2
)); ));
} }
@ -76,27 +76,27 @@ namespace ufo
return *this; return *this;
} }
std::vector<std::reference_wrapper<const UnfoldSolver::OutputType::QPointDataType>> PlotSolver::search_qpoints std::vector<std::reference_wrapper<const UnfoldSolver::OutputType::QpointDataType>> PlotSolver::search_qpoints
( (
const std::pair<Eigen::Vector3d, Eigen::Vector3d>& path, const std::pair<Eigen::Vector3d, Eigen::Vector3d>& path,
const decltype(InputType::SourceType::QPointData)& available_qpoints, const decltype(InputType::SourceType::QpointData)& available_qpoints,
double threshold, bool exclude_endpoint double threshold, bool exclude_endpoint
) )
{ {
std::multimap<double, std::reference_wrapper<const UnfoldSolver::OutputType::QPointDataType>> selected_qpoints; std::multimap<double, std::reference_wrapper<const UnfoldSolver::OutputType::QpointDataType>> selected_qpoints;
// 对于 output 中的每一个点, 检查这个点是否在路径上. 如果在, 把它加入到 selected_qpoints 中 // 对于 output 中的每一个点, 检查这个点是否在路径上. 如果在, 把它加入到 selected_qpoints 中
for (auto& qpoint : available_qpoints) for (auto& qpoint : available_qpoints)
{ {
// 计算三点围成的三角形的面积的两倍 // 计算三点围成的三角形的面积的两倍
auto area = (path.second - path.first).cross(qpoint.QPoint - path.first).norm(); auto area = (path.second - path.first).cross(qpoint.Qpoint - path.first).norm();
// 计算这个点到前两个点所在直线的距离 // 计算这个点到前两个点所在直线的距离
auto distance = area / (path.second - path.first).norm(); auto distance = area / (path.second - path.first).norm();
// 如果这个点到前两个点所在直线的距离小于阈值, 则认为这个点在路径上 // 如果这个点到前两个点所在直线的距离小于阈值, 则认为这个点在路径上
if (distance < threshold) if (distance < threshold)
{ {
// 计算这个点到前两个点的距离, 两个距离都应该小于两点之间的距离 // 计算这个点到前两个点的距离, 两个距离都应该小于两点之间的距离
auto distance1 = (qpoint.QPoint - path.first).norm(); auto distance1 = (qpoint.Qpoint - path.first).norm();
auto distance2 = (qpoint.QPoint - path.second).norm(); auto distance2 = (qpoint.Qpoint - path.second).norm();
auto distance3 = (path.second - path.first).norm(); auto distance3 = (path.second - path.first).norm();
if (distance1 < distance3 + threshold && distance2 < distance3 + threshold) if (distance1 < distance3 + threshold && distance2 < distance3 + threshold)
// 如果这个点不在终点处, 或者不排除终点, 则加入 // 如果这个点不在终点处, 或者不排除终点, 则加入
@ -117,7 +117,7 @@ namespace ufo
} }
if (selected_qpoints.empty()) if (selected_qpoints.empty())
throw std::runtime_error("No q points found"); throw std::runtime_error("No q points found");
std::vector<std::reference_wrapper<const UnfoldSolver::OutputType::QPointDataType>> result; std::vector<std::reference_wrapper<const UnfoldSolver::OutputType::QpointDataType>> result;
for (auto& qpoint : selected_qpoints) for (auto& qpoint : selected_qpoints)
result.push_back(qpoint.second); result.push_back(qpoint.second);
return result; return result;
@ -126,18 +126,18 @@ namespace ufo
std::vector<std::vector<double>> PlotSolver::calculate_values std::vector<std::vector<double>> PlotSolver::calculate_values
( (
const std::vector<std::pair<Eigen::Vector3d, Eigen::Vector3d>>& path, const std::vector<std::pair<Eigen::Vector3d, Eigen::Vector3d>>& path,
const std::vector<std::vector<std::reference_wrapper<const UnfoldSolver::OutputType::QPointDataType>>>& qpoints, const std::vector<std::vector<std::reference_wrapper<const UnfoldSolver::OutputType::QpointDataType>>>& qpoints,
const decltype(InputType::FigureConfigType::Resolution)& resolution, const decltype(InputType::FigureConfigType::Resolution)& resolution,
const decltype(InputType::FigureConfigType::Range)& range const decltype(InputType::FigureConfigType::Range)& range
) )
{ {
// 整理输入 // 整理输入
std::map<double, std::reference_wrapper<const UnfoldSolver::OutputType::QPointDataType>> qpoints_with_distance; std::map<double, std::reference_wrapper<const UnfoldSolver::OutputType::QpointDataType>> qpoints_with_distance;
double total_distance = 0; double total_distance = 0;
for (unsigned i = 0; i < path.size(); i++) for (unsigned i = 0; i < path.size(); i++)
{ {
for (auto& _ : qpoints[i]) for (auto& _ : qpoints[i])
qpoints_with_distance.emplace(total_distance + (_.get().QPoint - path[i].first).norm(), _); qpoints_with_distance.emplace(total_distance + (_.get().Qpoint - path[i].first).norm(), _);
total_distance += (path[i].second - path[i].first).norm(); total_distance += (path[i].second - path[i].first).norm();
} }
@ -145,8 +145,8 @@ namespace ufo
std::vector<std::vector<double>> values; std::vector<std::vector<double>> values;
auto blend = [] auto blend = []
( (
const UnfoldSolver::OutputType::QPointDataType& a, const UnfoldSolver::OutputType::QpointDataType& a,
const UnfoldSolver::OutputType::QPointDataType& b, const UnfoldSolver::OutputType::QpointDataType& b,
double ratio, unsigned resolution, std::pair<double, double> range double ratio, unsigned resolution, std::pair<double, double> range
) -> std::vector<double> ) -> std::vector<double>
{ {

View File

@ -52,32 +52,32 @@ namespace ufo
if (!std::set<std::string>{"yaml"}.contains(AtomPositionInputFile.Format)) if (!std::set<std::string>{"yaml"}.contains(AtomPositionInputFile.Format))
throw std::runtime_error(fmt::format throw std::runtime_error(fmt::format
("Unknown AtomPositionInputFile.Format: {}, should be \"yaml\".", AtomPositionInputFile.Format)); ("Unknown AtomPositionInputFile.Format: {}, should be \"yaml\".", AtomPositionInputFile.Format));
read_file_config(node["QPointDataInputFile"], QPointDataInputFile); read_file_config(node["QpointDataInputFile"], QpointDataInputFile);
if (!std::set<std::string>{"yaml", "hdf5"}.contains(QPointDataInputFile.Format)) if (!std::set<std::string>{"yaml", "hdf5"}.contains(QpointDataInputFile.Format))
throw std::runtime_error(fmt::format throw std::runtime_error(fmt::format
("Unknown QPointDataInputFile.Format: {}, should be \"yaml\" or \"hdf5\".", QPointDataInputFile.Format)); ("Unknown QpointDataInputFile.Format: {}, should be \"yaml\" or \"hdf5\".", QpointDataInputFile.Format));
if (auto value = node["QPointDataOutputFile"]) if (auto value = node["QpointDataOutputFile"])
{ {
QPointDataOutputFile.resize(value.size()); QpointDataOutputFile.resize(value.size());
for (unsigned i = 0; i < value.size(); i++) for (unsigned i = 0; i < value.size(); i++)
{ {
read_file_config(value[i], QPointDataOutputFile[i]); read_file_config(value[i], QpointDataOutputFile[i]);
if if
( (
QPointDataOutputFile[i].ExtraParameters.contains("SameAsConfigFile") QpointDataOutputFile[i].ExtraParameters.contains("SameAsConfigFile")
&& std::any_cast<bool>(QPointDataOutputFile[i].ExtraParameters["SameAsConfigFile"]) && std::any_cast<bool>(QpointDataOutputFile[i].ExtraParameters["SameAsConfigFile"])
) )
throw std::runtime_error("QPointDataOutputFile.SameAsConfigFile should not be set."); throw std::runtime_error("QpointDataOutputFile.SameAsConfigFile should not be set.");
if if
( (
!std::set<std::string>{"yaml", "yaml-human-readable", "zpp", "hdf5"} !std::set<std::string>{"yaml", "yaml-human-readable", "zpp", "hdf5"}
.contains(QPointDataOutputFile[i].Format) .contains(QpointDataOutputFile[i].Format)
) )
throw std::runtime_error(fmt::format throw std::runtime_error(fmt::format
( (
"Unknown QPointDataOutputFile[{}].Format: {}, should be " "Unknown QpointDataOutputFile[{}].Format: {}, should be "
"\"yaml\", \"yaml-human-readable\", \"zpp\" or \"hdf5\".", "\"yaml\", \"yaml-human-readable\", \"zpp\" or \"hdf5\".",
i, QPointDataOutputFile[i].Format i, QpointDataOutputFile[i].Format
)); ));
} }
} }
@ -99,20 +99,20 @@ namespace ufo
* SuperCellMultiplier.cast<double>().asDiagonal() * PrimativeCell).eval(); * SuperCellMultiplier.cast<double>().asDiagonal() * PrimativeCell).eval();
AtomPosition = atom_position_to_super_cell * super_cell; AtomPosition = atom_position_to_super_cell * super_cell;
} }
if (QPointDataInputFile.Format == "yaml") if (QpointDataInputFile.Format == "yaml")
{ {
auto node = YAML::LoadFile(QPointDataInputFile.FileName); auto node = YAML::LoadFile(QpointDataInputFile.FileName);
auto phonon = node["phonon"].as<std::vector<YAML::Node>>(); auto phonon = node["phonon"].as<std::vector<YAML::Node>>();
QPointData.resize(phonon.size()); QpointData.resize(phonon.size());
for (unsigned i = 0; i < phonon.size(); i++) for (unsigned i = 0; i < phonon.size(); i++)
{ {
for (unsigned j = 0; j < 3; j++) for (unsigned j = 0; j < 3; j++)
QPointData[i].QPoint(j) = phonon[i]["q-position"][j].as<double>(); QpointData[i].Qpoint(j) = phonon[i]["q-position"][j].as<double>();
auto band = phonon[i]["band"].as<std::vector<YAML::Node>>(); auto band = phonon[i]["band"].as<std::vector<YAML::Node>>();
QPointData[i].ModeData.resize(band.size()); QpointData[i].ModeData.resize(band.size());
for (unsigned j = 0; j < band.size(); j++) for (unsigned j = 0; j < band.size(); j++)
{ {
QPointData[i].ModeData[j].Frequency = band[j]["frequency"].as<double>(); QpointData[i].ModeData[j].Frequency = band[j]["frequency"].as<double>();
auto eigenvector_vectors = band[j]["eigenvector"] auto eigenvector_vectors = band[j]["eigenvector"]
.as<std::vector<std::vector<std::vector<double>>>>(); .as<std::vector<std::vector<std::vector<double>>>>();
Eigen::MatrixX3cd eigenvectors(AtomPosition.rows(), 3); Eigen::MatrixX3cd eigenvectors(AtomPosition.rows(), 3);
@ -122,18 +122,18 @@ namespace ufo
= eigenvector_vectors[k][l][0] + 1i * eigenvector_vectors[k][l][1]; = eigenvector_vectors[k][l][0] + 1i * eigenvector_vectors[k][l][1];
// 需要对读入的原子运动状态作相位转换, 使得它们与我们的约定一致(对超胞周期性重复) // 需要对读入的原子运动状态作相位转换, 使得它们与我们的约定一致(对超胞周期性重复)
// 这里还要需要做归一化处理 (指将数据简单地作为向量处理的归一化) // 这里还要需要做归一化处理 (指将数据简单地作为向量处理的归一化)
auto& AtomMovement = QPointData[i].ModeData[j].AtomMovement; auto& AtomMovement = QpointData[i].ModeData[j].AtomMovement;
// AtomMovement = eigenvectors.array().colwise() * (-2 * std::numbers::pi_v<double> * 1i // AtomMovement = eigenvectors.array().colwise() * (-2 * std::numbers::pi_v<double> * 1i
// * (atom_position_to_super_cell * input.QPointData[i].QPoint)).array().exp(); // * (atom_position_to_super_cell * input.QpointData[i].Qpoint)).array().exp();
// AtomMovement /= AtomMovement.norm(); // AtomMovement /= AtomMovement.norm();
// phonopy 似乎已经进行了相位的转换!为什么? // phonopy 似乎已经进行了相位的转换!为什么?
AtomMovement = eigenvectors / eigenvectors.norm(); AtomMovement = eigenvectors / eigenvectors.norm();
} }
} }
} }
else if (QPointDataInputFile.Format == "hdf5") else if (QpointDataInputFile.Format == "hdf5")
{ {
HighFive::File file(QPointDataInputFile.FileName, HighFive::File::ReadOnly); HighFive::File file(QpointDataInputFile.FileName, HighFive::File::ReadOnly);
auto size = file.getDataSet("/frequency").getDimensions(); auto size = file.getDataSet("/frequency").getDimensions();
auto frequency = file.getDataSet("/frequency") auto frequency = file.getDataSet("/frequency")
.read<std::vector<std::vector<std::vector<double>>>>(); .read<std::vector<std::vector<std::vector<double>>>>();
@ -141,28 +141,28 @@ namespace ufo
.read<std::vector<std::vector<std::vector<std::vector<PhonopyComplex>>>>>(); .read<std::vector<std::vector<std::vector<std::vector<PhonopyComplex>>>>>();
auto path = file.getDataSet("/path") auto path = file.getDataSet("/path")
.read<std::vector<std::vector<std::vector<double>>>>(); .read<std::vector<std::vector<std::vector<double>>>>();
QPointData.resize(size[0] * size[1]); QpointData.resize(size[0] * size[1]);
for (unsigned i = 0; i < size[0]; i++) for (unsigned i = 0; i < size[0]; i++)
for (unsigned j = 0; j < size[1]; j++) 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].Qpoint = Eigen::Vector3d(path[i][j].data());
QPointData[i * size[1] + j].ModeData.resize(size[2]); QpointData[i * size[1] + j].ModeData.resize(size[2]);
for (unsigned k = 0; k < size[2]; k++) for (unsigned k = 0; k < size[2]; k++)
{ {
QPointData[i * size[1] + j].ModeData[k].Frequency = frequency[i][j][k]; QpointData[i * size[1] + j].ModeData[k].Frequency = frequency[i][j][k];
Eigen::MatrixX3cd eigenvectors(AtomPosition.rows(), 3); Eigen::MatrixX3cd eigenvectors(AtomPosition.rows(), 3);
for (unsigned l = 0; l < AtomPosition.rows(); l++) for (unsigned l = 0; l < AtomPosition.rows(); l++)
for (unsigned m = 0; m < 3; m++) for (unsigned m = 0; m < 3; m++)
eigenvectors(l, m) eigenvectors(l, m)
= eigenvector_vector[i][j][l * 3 + m][k].r + eigenvector_vector[i][j][l * 3 + m][k].i * 1i; = 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(); QpointData[i * size[1] + j].ModeData[k].AtomMovement = eigenvectors / eigenvectors.norm();
} }
} }
} }
} }
void UnfoldSolver::OutputType::write void UnfoldSolver::OutputType::write
(decltype(InputType::QPointDataOutputFile) output_files) const (decltype(InputType::QpointDataOutputFile) output_files) const
{ {
for (auto& output_file : output_files) for (auto& output_file : output_files)
write(output_file.FileName, output_file.Format); write(output_file.FileName, output_file.Format);
@ -173,11 +173,11 @@ namespace ufo
std::ofstream(filename) << [&] std::ofstream(filename) << [&]
{ {
std::stringstream print; std::stringstream print;
print << "QPointData:\n"; print << "QpointData:\n";
for (auto& qpoint: QPointData) for (auto& qpoint: QpointData)
{ {
print << fmt::format(" - QPoint: [ {1:.{0}f}, {2:.{0}f}, {3:.{0}f} ]\n", print << fmt::format(" - Qpoint: [ {1:.{0}f}, {2:.{0}f}, {3:.{0}f} ]\n",
percision, qpoint.QPoint[0], qpoint.QPoint[1], qpoint.QPoint[2]); percision, qpoint.Qpoint[0], qpoint.Qpoint[1], qpoint.Qpoint[2]);
print << fmt::format(" Source: [ {1:.{0}f}, {2:.{0}f}, {3:.{0}f} ]\n", print << fmt::format(" Source: [ {1:.{0}f}, {2:.{0}f}, {3:.{0}f} ]\n",
percision, qpoint.Source[0], qpoint.Source[1], qpoint.Source[2]); percision, qpoint.Source[0], qpoint.Source[1], qpoint.Source[2]);
print << " ModeData:\n"; print << " ModeData:\n";
@ -190,9 +190,9 @@ namespace ufo
else if (format == "yaml-human-readable") else if (format == "yaml-human-readable")
{ {
std::remove_cvref_t<decltype(*this)> output; std::remove_cvref_t<decltype(*this)> output;
std::map<unsigned, std::vector<decltype(QPointData)::const_iterator>> std::map<unsigned, std::vector<decltype(QpointData)::const_iterator>>
meta_qpoint_to_sub_qpoint_iterators; meta_qpoint_to_sub_qpoint_iterators;
for (auto it = QPointData.begin(); it != QPointData.end(); it++) for (auto it = QpointData.begin(); it != QpointData.end(); it++)
meta_qpoint_to_sub_qpoint_iterators[it->SourceIndex_].push_back(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 [meta_qpoint_index, sub_qpoint_iterators] : meta_qpoint_to_sub_qpoint_iterators)
for (auto& qpoint : sub_qpoint_iterators) for (auto& qpoint : sub_qpoint_iterators)
@ -218,8 +218,8 @@ namespace ufo
frequency_to_weight[frequency_sum / weight_sum] = weight_sum; frequency_to_weight[frequency_sum / weight_sum] = weight_sum;
} }
} }
auto& _ = output.QPointData.emplace_back(); auto& _ = output.QpointData.emplace_back();
_.QPoint = qpoint->QPoint; _.Qpoint = qpoint->Qpoint;
_.Source = qpoint->Source; _.Source = qpoint->Source;
_.SourceIndex_ = qpoint->SourceIndex_; _.SourceIndex_ = qpoint->SourceIndex_;
for (auto [frequency, weight] : frequency_to_weight) for (auto [frequency, weight] : frequency_to_weight)
@ -247,9 +247,9 @@ namespace ufo
std::vector<std::vector<double>> Source; std::vector<std::vector<double>> Source;
std::vector<std::vector<double>> Frequency; std::vector<std::vector<double>> Frequency;
std::vector<std::vector<double>> Weight; std::vector<std::vector<double>> Weight;
for (auto& qpoint : QPointData) for (auto& qpoint : QpointData)
{ {
Qpoint.emplace_back(qpoint.QPoint.data(), qpoint.QPoint.data() + 3); Qpoint.emplace_back(qpoint.Qpoint.data(), qpoint.Qpoint.data() + 3);
Source.emplace_back(qpoint.Source.data(), qpoint.Source.data() + 3); Source.emplace_back(qpoint.Source.data(), qpoint.Source.data() + 3);
Frequency.emplace_back(); Frequency.emplace_back();
Weight.emplace_back(); Weight.emplace_back();
@ -262,7 +262,7 @@ namespace ufo
HighFive::File file(filename, HighFive::File file(filename,
HighFive::File::ReadWrite | HighFive::File::Create | HighFive::File::Truncate); HighFive::File::ReadWrite | HighFive::File::Create | HighFive::File::Truncate);
file.createDataSet<double>("QPoint", file.createDataSet<double>("Qpoint",
HighFive::DataSpace::From(Qpoint)).write(Qpoint); HighFive::DataSpace::From(Qpoint)).write(Qpoint);
file.createDataSet<double>("Source", file.createDataSet<double>("Source",
HighFive::DataSpace::From(Source)).write(Source); HighFive::DataSpace::From(Source)).write(Source);
@ -298,8 +298,8 @@ namespace ufo
{ {
std::clog << "Calculating projection coefficient... " << std::flush; std::clog << "Calculating projection coefficient... " << std::flush;
std::vector<std::reference_wrapper<const decltype std::vector<std::reference_wrapper<const decltype
(InputType::QPointDataType::ModeDataType::AtomMovement)>> mode_data; (InputType::QpointDataType::ModeDataType::AtomMovement)>> mode_data;
for (auto& qpoint : Input_.QPointData) for (auto& qpoint : Input_.QpointData)
for (auto& mode : qpoint.ModeData) for (auto& mode : qpoint.ModeData)
mode_data.emplace_back(mode.AtomMovement); mode_data.emplace_back(mode.AtomMovement);
std::atomic<unsigned> number_of_finished_modes(0); std::atomic<unsigned> number_of_finished_modes(0);
@ -321,12 +321,12 @@ namespace ufo
std::clog << "\33[2K\rCalculating projection coefficient... Done." << std::endl; std::clog << "\33[2K\rCalculating projection coefficient... Done." << std::endl;
std::clog << "Constructing output... " << std::flush; std::clog << "Constructing output... " << std::flush;
std::vector<std::reference_wrapper<const decltype(InputType::QPointDataType::QPoint)>> qpoint; std::vector<std::reference_wrapper<const decltype(InputType::QpointDataType::Qpoint)>> qpoint;
std::vector<std::vector<std::reference_wrapper<const std::vector<std::vector<std::reference_wrapper<const
decltype(InputType::QPointDataType::ModeDataType::Frequency)>>> frequency; decltype(InputType::QpointDataType::ModeDataType::Frequency)>>> frequency;
for (auto& qpoint_data : Input_.QPointData) for (auto& qpoint_data : Input_.QpointData)
{ {
qpoint.emplace_back(qpoint_data.QPoint); qpoint.emplace_back(qpoint_data.Qpoint);
frequency.emplace_back(); frequency.emplace_back();
for (auto& mode_data : qpoint_data.ModeData) for (auto& mode_data : qpoint_data.ModeData)
frequency.back().emplace_back(mode_data.Frequency); frequency.back().emplace_back(mode_data.Frequency);
@ -339,7 +339,7 @@ namespace ufo
std::clog << "Done." << std::endl; std::clog << "Done." << std::endl;
} }
std::clog << "Writing output... " << std::flush; std::clog << "Writing output... " << std::flush;
Output_->write(Input_.QPointDataOutputFile); Output_->write(Input_.QpointDataOutputFile);
std::clog << "Done." << std::endl; std::clog << "Done." << std::endl;
return *this; return *this;
} }
@ -380,7 +380,7 @@ namespace ufo
( (
const BasisType& basis, const BasisType& basis,
const std::vector<std::reference_wrapper<const decltype const std::vector<std::reference_wrapper<const decltype
(InputType::QPointDataType::ModeDataType::AtomMovement)>>& mode_data, (InputType::QpointDataType::ModeDataType::AtomMovement)>>& mode_data,
std::atomic<unsigned>& number_of_finished_modes std::atomic<unsigned>& number_of_finished_modes
) )
{ {
@ -418,16 +418,16 @@ namespace ufo
const decltype(InputType::SuperCellMultiplier)& super_cell_multiplier, const decltype(InputType::SuperCellMultiplier)& super_cell_multiplier,
const decltype(InputType::SuperCellDeformation)& super_cell_deformation, const decltype(InputType::SuperCellDeformation)& super_cell_deformation,
const std::vector<std::reference_wrapper<const decltype const std::vector<std::reference_wrapper<const decltype
(InputType::QPointDataType::QPoint)>>& qpoint, (InputType::QpointDataType::Qpoint)>>& qpoint,
const std::vector<std::vector<std::reference_wrapper<const decltype const std::vector<std::vector<std::reference_wrapper<const decltype
(InputType::QPointDataType::ModeDataType::Frequency)>>>& frequency, (InputType::QpointDataType::ModeDataType::Frequency)>>>& frequency,
const ProjectionCoefficientType_& projection_coefficient const ProjectionCoefficientType_& projection_coefficient
) )
{ {
OutputType output; OutputType output;
for (unsigned i_of_qpoint = 0, num_of_mode_manipulated = 0; i_of_qpoint < qpoint.size(); i_of_qpoint++) for (unsigned i_of_qpoint = 0, num_of_mode_manipulated = 0; i_of_qpoint < qpoint.size(); i_of_qpoint++)
{ {
// 当 SuperCellDeformation 不是单位矩阵时, input.QPointData[i_of_qpoint].QPoint 不一定在 reciprocal_primative_cell 中 // 当 SuperCellDeformation 不是单位矩阵时, input.QpointData[i_of_qpoint].Qpoint 不一定在 reciprocal_primative_cell 中
// 需要首先将 q 点平移数个周期, 进入不包含 SuperCellDeformation 的超胞 (称为 ModifiedSupreCell) 的倒格子中 // 需要首先将 q 点平移数个周期, 进入不包含 SuperCellDeformation 的超胞 (称为 ModifiedSupreCell) 的倒格子中
auto qpoint_by_reciprocal_modified_super_cell_in_reciprocal_modified_super_cell auto qpoint_by_reciprocal_modified_super_cell_in_reciprocal_modified_super_cell
= !super_cell_deformation ? qpoint[i_of_qpoint].get() : [&] = !super_cell_deformation ? qpoint[i_of_qpoint].get() : [&]
@ -483,7 +483,7 @@ namespace ufo
for (auto [xyz_of_diff_of_sub_qpoint_by_reciprocal_modified_super_cell, i_of_sub_qpoint] for (auto [xyz_of_diff_of_sub_qpoint_by_reciprocal_modified_super_cell, i_of_sub_qpoint]
: triplet_sequence(super_cell_multiplier)) : triplet_sequence(super_cell_multiplier))
{ {
auto& _ = output.QPointData.emplace_back(); auto& _ = output.QpointData.emplace_back();
/* /*
SubQpointByReciprocalModifiedSuperCell = XyzOfDiffOfSubQpointByReciprocalModifiedSuperCell + SubQpointByReciprocalModifiedSuperCell = XyzOfDiffOfSubQpointByReciprocalModifiedSuperCell +
QpointInReciprocalModifiedSuperCellByReciprocalModifiedSuperCell; QpointInReciprocalModifiedSuperCellByReciprocalModifiedSuperCell;
@ -497,7 +497,7 @@ namespace ufo
(XyzOfDiffOfSubQpointByReciprocalModifiedSuperCell + (XyzOfDiffOfSubQpointByReciprocalModifiedSuperCell +
QpointInReciprocalModifiedSuperCellByReciprocalModifiedSuperCell); QpointInReciprocalModifiedSuperCellByReciprocalModifiedSuperCell);
*/ */
_.QPoint = super_cell_multiplier.cast<double>().cwiseInverse().asDiagonal() _.Qpoint = super_cell_multiplier.cast<double>().cwiseInverse().asDiagonal()
* (xyz_of_diff_of_sub_qpoint_by_reciprocal_modified_super_cell.cast<double>() * (xyz_of_diff_of_sub_qpoint_by_reciprocal_modified_super_cell.cast<double>()
+ qpoint_by_reciprocal_modified_super_cell_in_reciprocal_modified_super_cell); + qpoint_by_reciprocal_modified_super_cell_in_reciprocal_modified_super_cell);
_.Source = qpoint[i_of_qpoint]; _.Source = qpoint[i_of_qpoint];

14
test/14.2.5.5.plot.yaml Normal file
View File

@ -0,0 +1,14 @@
PrimativeCell:
- [ 1.548010265167714, -2.681232429908652, 0.000000000000000 ] # a
- [ 1.548010265167714, 2.681232429908652, 0.000000000000000 ] # b
- [ 0.000000000000000, 0.000000000000000, 5.061224103556595 ] # c
SourceFilename: ./test/14.2.5.5.result.zpp
Figures:
- Qpoints:
- - [ 0, 0, 0 ]
- [ 0.5, 0, 0 ]
- [ 0.3333333, 0.3333333, 0 ]
- [ 0, 0, 0 ]
Resolution: [ 1024, 400 ]
Range: [ -5, 35 ]
Filename: ./test/14.2.5.5.png

18
test/14.2.5.5.unfold.yaml Normal file
View File

@ -0,0 +1,18 @@
PrimativeCell:
- [ 1.548010265167714, -2.681232429908652, 0.000000000000000 ] # a
- [ 1.548010265167714, 2.681232429908652, 0.000000000000000 ] # b
- [ 0.000000000000000, 0.000000000000000, 5.061224103556595 ] # c
SuperCellMultiplier: [ 1, 2, 1 ]
SuperCellDeformation:
- [ 1, 0, 0]
- [ 1, 1, 0]
- [ 0, 0, 1]
PrimativeCellBasisNumber: [ 8, 8, 8 ]
AtomPositionInputFile:
FileName: "/home/chn/Documents/lammps-SiC/14/14.2/14.2.5/14.2.5.5/phonopy.yaml"
Format: "yaml"
QPointDataInputFile: { FileName: "/home/chn/Documents/lammps-SiC/14/14.2/14.2.5/14.2.5.5/band.hdf5", Format: "hdf5" }
QPointDataOutputFile:
- { FileName: "test/14.2.5.5.result.yaml", Format: "yaml" }
- { FileName: "test/14.2.5.5.result.human-readable.yaml", Format: "yaml-human-readable" }
- { FileName: "test/14.2.5.5.result.zpp", Format: "zpp" }

14
test/14.2.6.4.plot.yaml Normal file
View File

@ -0,0 +1,14 @@
PrimativeCell:
- [ 1.548010265167714, -2.681232429908652, 0.000000000000000 ] # a
- [ 1.548010265167714, 2.681232429908652, 0.000000000000000 ] # b
- [ 0.000000000000000, 0.000000000000000, 5.061224103556595 ] # c
SourceFilename: ./test/14.2.6.4.result.zpp
Figures:
- Qpoints:
- - [ 0, 0, 0 ]
- [ 0.5, 0, 0 ]
- [ 0.3333333, 0.3333333, 0 ]
- [ 0, 0, 0 ]
Resolution: [ 1024, 400 ]
Range: [ -5, 35 ]
Filename: ./test/14.2.6.4.png

18
test/14.2.6.4.unfold.yaml Normal file
View File

@ -0,0 +1,18 @@
PrimativeCell:
- [ 1.548010265167714, -2.681232429908652, 0.000000000000000 ] # a
- [ 1.548010265167714, 2.681232429908652, 0.000000000000000 ] # b
- [ 0.000000000000000, 0.000000000000000, 5.061224103556595 ] # c
SuperCellMultiplier: [ 4, 8, 1 ]
SuperCellDeformation:
- [ 1, 0, 0]
- [ 1, 1, 0]
- [ 0, 0, 1]
PrimativeCellBasisNumber: [ 8, 8, 8 ]
AtomPositionInputFile:
FileName: "/home/chn/Documents/lammps-SiC/14/14.2/14.2.6/14.2.6.4/phonopy.yaml"
Format: "yaml"
QPointDataInputFile: { FileName: "/home/chn/Documents/lammps-SiC/14/14.2/14.2.6/14.2.6.4/band.hdf5", Format: "hdf5" }
QPointDataOutputFile:
- { FileName: "test/14.2.6.4.result.yaml", Format: "yaml" }
- { FileName: "test/14.2.6.4.result.human-readable.yaml", Format: "yaml-human-readable" }
- { FileName: "test/14.2.6.4.result.zpp", Format: "zpp" }