mirror of
https://github.com/CHN-beta/nixos.git
synced 2026-01-12 04:39:23 +08:00
packages.ufo: split repository
This commit is contained in:
17
flake.lock
generated
17
flake.lock
generated
@@ -1309,6 +1309,7 @@
|
||||
"spectroscopy": "spectroscopy",
|
||||
"sqlite-orm": "sqlite-orm",
|
||||
"tgbot-cpp": "tgbot-cpp",
|
||||
"ufo": "ufo",
|
||||
"v-sim": "v-sim",
|
||||
"vaspberry": "vaspberry",
|
||||
"winapps": "winapps",
|
||||
@@ -1623,6 +1624,22 @@
|
||||
"type": "github"
|
||||
}
|
||||
},
|
||||
"ufo": {
|
||||
"flake": false,
|
||||
"locked": {
|
||||
"lastModified": 1732177086,
|
||||
"narHash": "sha256-zmrzTQGXkR54igJUhYp0pFqS2RdV69Wi/wgyFME/K+E=",
|
||||
"ref": "refs/heads/main",
|
||||
"rev": "28e4d29f2c70d1f3b80a092b75b81a4793455980",
|
||||
"revCount": 65,
|
||||
"type": "git",
|
||||
"url": "https://git.chn.moe/chn/ufo.git"
|
||||
},
|
||||
"original": {
|
||||
"type": "git",
|
||||
"url": "https://git.chn.moe/chn/ufo.git"
|
||||
}
|
||||
},
|
||||
"v-sim": {
|
||||
"flake": false,
|
||||
"locked": {
|
||||
|
||||
@@ -69,6 +69,7 @@
|
||||
nixos-wallpaper = { url = "git+https://git.chn.moe/chn/nixos-wallpaper.git"; flake = false; };
|
||||
spectroscopy = { url = "github:skelton-group/Phonopy-Spectroscopy"; flake = false; };
|
||||
vaspberry = { url = "github:Infant83/VASPBERRY"; flake = false; };
|
||||
ufo = { url = "git+https://git.chn.moe/chn/ufo.git"; flake = false; };
|
||||
};
|
||||
|
||||
outputs = inputs: let localLib = import ./flake/lib.nix inputs.nixpkgs.lib; in
|
||||
|
||||
@@ -76,7 +76,7 @@ inputs: rec
|
||||
sqlite-orm = inputs.pkgs.callPackage ./sqlite-orm.nix { src = inputs.topInputs.sqlite-orm; };
|
||||
mkPnpmPackage = inputs.pkgs.callPackage ./mkPnpmPackage.nix {};
|
||||
sbatch-tui = inputs.pkgs.callPackage ./sbatch-tui { inherit biu; stdenv = inputs.pkgs.clang18Stdenv; };
|
||||
ufo = inputs.pkgs.callPackage ./ufo
|
||||
ufo = inputs.pkgs.callPackage inputs.topInputs.ufo
|
||||
{
|
||||
inherit biu matplotplusplus;
|
||||
tbb = inputs.pkgs.tbb_2021_11;
|
||||
|
||||
@@ -1 +0,0 @@
|
||||
use flake .#ufo
|
||||
@@ -1,31 +0,0 @@
|
||||
cmake_minimum_required(VERSION 3.14)
|
||||
project(ufo VERSION 0 LANGUAGES CXX)
|
||||
enable_testing()
|
||||
include(GNUInstallDirs)
|
||||
|
||||
if(NOT CMAKE_BUILD_TYPE AND NOT CMAKE_CONFIGURATION_TYPES)
|
||||
message("Setting build type to 'Release' as none was specified.")
|
||||
set(CMAKE_BUILD_TYPE Release CACHE STRING "Choose the type of build." FORCE)
|
||||
set_property(CACHE CMAKE_BUILD_TYPE PROPERTY STRINGS "Debug" "Release" "MinSizeRel" "RelWithDebInfo")
|
||||
endif()
|
||||
|
||||
find_package(TBB REQUIRED)
|
||||
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-create-displacement.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)
|
||||
target_compile_options(ufo PRIVATE -fexperimental-library)
|
||||
|
||||
install(TARGETS ufo RUNTIME DESTINATION ${CMAKE_INSTALL_BINDIR})
|
||||
|
||||
get_property(ImportedTargets DIRECTORY "${CMAKE_SOURCE_DIR}" PROPERTY IMPORTED_TARGETS)
|
||||
message("Imported targets: ${ImportedTargets}")
|
||||
message("List of compile features: ${CMAKE_CXX_COMPILE_FEATURES}")
|
||||
|
||||
include(CTest)
|
||||
add_test(NAME fold COMMAND ufo fold ${PROJECT_SOURCE_DIR}/test/fold/config.yaml)
|
||||
|
||||
@@ -1,11 +0,0 @@
|
||||
{
|
||||
stdenv, cmake, pkg-config, version ? null,
|
||||
tbb, matplotplusplus, biu
|
||||
}: stdenv.mkDerivation
|
||||
{
|
||||
name = "ufo";
|
||||
src = ./.;
|
||||
buildInputs = [ tbb matplotplusplus biu ];
|
||||
nativeBuildInputs = [ cmake pkg-config ];
|
||||
doCheck = true;
|
||||
}
|
||||
@@ -1,45 +0,0 @@
|
||||
分为几个功能：
|
||||
|
||||
* fold：根据要计算的单胞的 q 点路径，计算超胞中对应的 q 点路径，生成的路径再交给 phonopy 计算。
|
||||
* unfold：根据 phonopy 计算的结果，将超胞的结果展开到单胞中。
|
||||
* plot：对计算结果画图。
|
||||
|
||||
主要的输入输出格式均为 yaml。对于数据特别大的情况，也可以从 hdf5 中读取一部分数据或者将一部分数据写入到 hdf5 文件中。
|
||||
|
||||
# fold
|
||||
|
||||
## 输入
|
||||
|
||||
```yaml
|
||||
# 三个整数组成的向量，表示从单胞到超胞，三个晶格矢量的倍数
|
||||
# 必写
|
||||
SuperCellMultiplier: [2, 2, 2]
|
||||
# 一个变换矩阵，表明超胞经历了怎样的扭曲。
|
||||
# 可选，默认值为单位矩阵
|
||||
SuperCellDeformation: [[1, 0, 0], [0, 1, 0], [0, 0, 1]]
|
||||
# 一个由三个浮点数组成的向量，表示考虑的 q 点
|
||||
# 必写
|
||||
Qpoints:
|
||||
- [0, 0, 0]
|
||||
- [0.1, 0, 0]
|
||||
- [0.2, 0, 0]
|
||||
- [0.3, 0, 0]
|
||||
- [0.4, 0, 0]
|
||||
- [0.5, 0, 0]
|
||||
# 一个 DataFile 类型的对象，表明输出结果到哪个文件
|
||||
# 必写
|
||||
OutputFile:
|
||||
```
|
||||
|
||||
## 输出
|
||||
|
||||
```yaml
|
||||
# 得到的 q 点坐标
|
||||
Qpoints:
|
||||
- [0, 0, 0]
|
||||
- [0.1, 0, 0]
|
||||
- [0.2, 0, 0]
|
||||
- [0.3, 0, 0]
|
||||
- [0.4, 0, 0]
|
||||
- [0.5, 0, 0]
|
||||
```
|
||||
@@ -1,74 +0,0 @@
|
||||
# pragma once
|
||||
# include <biu.hpp>
|
||||
|
||||
namespace ufo
|
||||
{
|
||||
// 在相位中, 约定为使用 $\exp (2 \pi i \vec{q} \cdot \vec{r})$ 来表示原子的运动状态
|
||||
// (而不是 $\exp (-2 \pi i \vec{q} \cdot \vec{r})$)
|
||||
// 一些书定义的倒格矢中包含了 $2 \pi$ 的部分, 我们这里约定不包含这部分.
|
||||
// 也就是说, 正格子与倒格子的转置相乘, 得到单位矩阵.
|
||||
|
||||
using namespace biu::literals;
|
||||
using namespace biu::stream_operators;
|
||||
|
||||
void fold(std::string config_file);
|
||||
void unfold(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);
|
||||
void raman_apply_contribution(std::string config_file);
|
||||
|
||||
// 许多函数都需要用到这个，所以写到头文件中
|
||||
struct UnfoldOutput
|
||||
{
|
||||
Eigen::Matrix3d PrimativeCell;
|
||||
Eigen::Matrix3i SuperCellTransformation;
|
||||
Eigen::Vector3i SuperCellMultiplier;
|
||||
Eigen::Matrix3d SuperCellDeformation;
|
||||
std::optional<std::vector<std::size_t>> SelectedAtoms;
|
||||
|
||||
// 关于各个 Q 点的数据
|
||||
struct QpointDataType
|
||||
{
|
||||
// Q 点的坐标，单位为单胞的倒格矢
|
||||
Eigen::Vector3d Qpoint;
|
||||
|
||||
// 来源于哪个 Q 点, 单位为超胞的倒格矢
|
||||
Eigen::Vector3d Source;
|
||||
std::size_t SourceIndex;
|
||||
|
||||
// 关于这个 Q 点上各个模式的数据
|
||||
struct ModeDataType
|
||||
{
|
||||
// 模式的频率，单位为 THz
|
||||
double Frequency;
|
||||
// 模式的权重
|
||||
double Weight;
|
||||
};
|
||||
std::vector<ModeDataType> ModeData;
|
||||
};
|
||||
std::vector<QpointDataType> QpointData;
|
||||
|
||||
struct MetaQpointDataType
|
||||
{
|
||||
// Q 点的坐标，单位为单胞的倒格矢
|
||||
Eigen::Vector3d Qpoint;
|
||||
|
||||
// 关于这个 Q 点上各个模式的数据
|
||||
struct ModeDataType
|
||||
{
|
||||
// 模式的频率，单位为 THz
|
||||
double Frequency;
|
||||
// 模式中各个原子的运动状态
|
||||
// 这个数据应当是这样得到的：动态矩阵的 eigenvector 乘以 $\exp(-2 \pi i \vec q \cdot \vec r)$
|
||||
// 这个数据可以认为是原子位移中, 关于超胞有周期性的那一部分, 再乘以原子质量的开方.
|
||||
// 这个数据会在 unfold 时被归一化
|
||||
Eigen::MatrixX3cd AtomMovement;
|
||||
};
|
||||
std::vector<ModeDataType> ModeData;
|
||||
};
|
||||
std::vector<MetaQpointDataType> MetaQpointData;
|
||||
|
||||
using serialize = zpp::bits::members<7>;
|
||||
};
|
||||
}
|
||||
@@ -1,58 +0,0 @@
|
||||
# include <ufo.hpp>
|
||||
|
||||
void ufo::fold(std::string config_file)
|
||||
{
|
||||
struct Input
|
||||
{
|
||||
Eigen::Matrix3d SuperCellDeformation;
|
||||
Eigen::Vector3i SuperCellMultiplier;
|
||||
std::vector<Eigen::Vector3d> Qpoints;
|
||||
std::optional<std::string> OutputFile;
|
||||
};
|
||||
struct Output
|
||||
{
|
||||
std::vector<Eigen::Vector3d> Qpoints;
|
||||
};
|
||||
auto fold = []
|
||||
(
|
||||
Eigen::Vector3d qpoint_in_reciprocal_primitive_cell_by_reciprocal_primitive_cell,
|
||||
Eigen::Matrix3d super_cell_transformation
|
||||
) -> Eigen::Vector3d
|
||||
{
|
||||
/*
|
||||
首先需要将 q 点坐标的单位转换为 ModifiedSuperCell 的格矢，可知：
|
||||
QpointByReciprocalModifiedSuperCell = SuperCellMultiplier * QpointByReciprocalPrimitiveCell;
|
||||
接下来考虑将 q 点坐标的单位转换为 SuperCell 的格矢
|
||||
ModifiedSuperCell = SuperCellMultiplier * PrimativeCell;
|
||||
SuperCell = SuperCellDeformation * ModifiedSuperCell;
|
||||
ReciprocalModifiedSuperCell = ModifiedSuperCell.inverse().transpose();
|
||||
ReciprocalSuperCell = SuperCell.inverse().transpose();
|
||||
Qpoint = QpointByReciprocalModifiedSuperCell.transpose() * ReciprocalModifiedSuperCell;
|
||||
Qpoint = QpointByReciprocalSuperCell.transpose() * ReciprocalSuperCell;
|
||||
整理可以得到:
|
||||
QpointByReciprocalSuperCell = SuperCellDeformation * QpointByReciprocalModifiedSuperCell;
|
||||
两个式子结合，可以得到：
|
||||
QpointByReciprocalSuperCell = SuperCellDeformation * SuperCellMultiplier * QpointByReciprocalPrimitiveCell;
|
||||
*/
|
||||
auto qpoint_by_reciprocal_super_cell =
|
||||
(
|
||||
super_cell_transformation * qpoint_in_reciprocal_primitive_cell_by_reciprocal_primitive_cell
|
||||
).eval();
|
||||
/*
|
||||
到目前为止，我们还没有移动过 q 点的坐标。现在，我们将它移动整数个 ReciprocalSuperCell，直到它落在超胞的倒格子中。
|
||||
这等价于直接取 QpointByReciprocalSuperCell - QpointByReciprocalSuperCell.floor()。
|
||||
*/
|
||||
return (qpoint_by_reciprocal_super_cell.array() - qpoint_by_reciprocal_super_cell.array().floor()).matrix();
|
||||
};
|
||||
auto input = YAML::LoadFile(config_file).as<Input>();
|
||||
Output output;
|
||||
output.Qpoints = input.Qpoints
|
||||
| ranges::views::transform([&](auto& qpoint)
|
||||
{
|
||||
return fold(qpoint, input.SuperCellDeformation * input.SuperCellMultiplier.cast<double>().asDiagonal());
|
||||
})
|
||||
| ranges::to_vector;
|
||||
|
||||
// 默认的输出太丑了，但是不想手动写了，忍一下
|
||||
std::ofstream(input.OutputFile.value_or("output.yaml")) << YAML::Node(output);
|
||||
}
|
||||
@@ -1,14 +0,0 @@
|
||||
# include <ufo.hpp>
|
||||
|
||||
int main(int argc, const char** argv)
|
||||
{
|
||||
using namespace biu::literals;
|
||||
if (argc != 3) throw std::runtime_error("Usage: {} task config.yaml"_f(argv[0]));
|
||||
if (argv[1] == "fold"s) ufo::fold(argv[2]);
|
||||
else if (argv[1] == "unfold"s) ufo::unfold(argv[2]);
|
||||
else if (argv[1] == "raman-create-displacement"s) ufo::raman_create_displacement(argv[2]);
|
||||
else if (argv[1] == "raman-apply-contribution"s);
|
||||
else if (argv[1] == "plot-band"s) ufo::plot_band(argv[2]);
|
||||
else if (argv[1] == "plot-point"s) ufo::plot_point(argv[2]);
|
||||
else throw std::runtime_error("Unknown task: {}"_f(argv[1]));
|
||||
}
|
||||
@@ -1,431 +0,0 @@
|
||||
# include <ufo.hpp>
|
||||
# include <matplot/matplot.h>
|
||||
# include <boost/container/flat_map.hpp>
|
||||
|
||||
void ufo::plot_band(std::string config_file)
|
||||
{
|
||||
struct Input
|
||||
{
|
||||
std::string UnfoldedDataFile;
|
||||
// 要画图的 q 点路径列表
|
||||
// 内层表示一个路径上的 q 点，外层表示不同的路径
|
||||
// 单位为倒格矢
|
||||
std::vector<std::vector<Eigen::Vector3d>> Qpoints;
|
||||
// 插值时使用的分辨率（不影响画出来图片的分辨率和横纵比）
|
||||
std::array<std::size_t, 2> InterpolationResolution;
|
||||
// 画图区域的y轴和x轴的比例。如果不指定，则由matplot++自动调整（通常调整为正方形，即 1）
|
||||
std::optional<double> AspectRatio;
|
||||
// 整张图片的分辨率
|
||||
std::optional<std::array<std::size_t, 2>> PictureResolution;
|
||||
// 画图的频率范围
|
||||
std::array<double, 2> FrequencyRange;
|
||||
// 搜索 q 点时的阈值，单位为埃^-1
|
||||
std::optional<double> ThresholdWhenSearchingQpoints;
|
||||
// 是否要在 y 轴上作一些标记
|
||||
std::optional<std::vector<std::pair<double, std::string>>> YTicks;
|
||||
// 是否输出图片
|
||||
std::optional<std::string> OutputPictureFile;
|
||||
// 是否输出数据，可以进一步使用 matplotlib 画图
|
||||
std::optional<std::string> OutputDataFile;
|
||||
};
|
||||
|
||||
// 根据 q 点路径, 搜索要使用的 q 点，返回的是 q 点在 QpointData 中的索引以及到路径起点的距离，以及这段路径的总长度
|
||||
auto search_qpoints = []
|
||||
(
|
||||
const Eigen::Matrix3d& primative_cell,
|
||||
const std::pair<Eigen::Vector3d, Eigen::Vector3d>& path,
|
||||
const std::vector<Eigen::Vector3d>& qpoints,
|
||||
double threshold, bool exclude_endpoint = false
|
||||
)
|
||||
{
|
||||
// 对于 output 中的每一个点, 检查这个点是否在路径上. 如果在, 把它加入到 selected_qpoints 中
|
||||
// 键为这个点到起点的距离
|
||||
boost::container::flat_map<double, std::size_t> selected_qpoints;
|
||||
auto begin = (path.first.transpose() * primative_cell.reverse()).transpose().eval();
|
||||
auto end = (path.second.transpose() * primative_cell.reverse()).transpose().eval();
|
||||
for (std::size_t i = 0; i < qpoints.size(); i++)
|
||||
for (auto cell_shift
|
||||
: biu::sequence(Eigen::Vector3i(-1, -1, -1), Eigen::Vector3i(2, 2, 2)))
|
||||
{
|
||||
auto qpoint
|
||||
= ((qpoints[i] + cell_shift.first.cast<double>()).transpose() * primative_cell.reverse()).transpose().eval();
|
||||
// 计算这个点到前两个点所在直线的距离
|
||||
auto distance = (end - begin).cross(qpoint - begin).norm()
|
||||
/ (path.second - path.first).norm();
|
||||
// 如果这个点到前两个点所在直线的距离小于阈值, 则认为这个点在这条直线上，但不一定在这两个点之间
|
||||
if (distance < threshold)
|
||||
{
|
||||
// 计算这个点到前两个点的距离, 两个距离都应该小于两点之间的距离
|
||||
auto distance1 = (qpoint - begin).norm();
|
||||
auto distance2 = (qpoint - end).norm();
|
||||
auto distance3 = (end - begin).norm();
|
||||
if (distance1 < distance3 + threshold && distance2 < distance3 + threshold)
|
||||
// 如果这个点不在终点处, 或者不排除终点, 则加入
|
||||
if (distance2 > threshold || !exclude_endpoint) selected_qpoints.emplace(distance1, i);
|
||||
}
|
||||
}
|
||||
// 去除非常接近的点
|
||||
for (auto it = selected_qpoints.begin(); it != selected_qpoints.end();)
|
||||
{
|
||||
auto next = std::next(it);
|
||||
if (next == selected_qpoints.end()) break;
|
||||
else if (next->first - it->first < threshold) selected_qpoints.erase(next);
|
||||
else it = next;
|
||||
}
|
||||
if (selected_qpoints.empty()) throw std::runtime_error("No q points found");
|
||||
return std::make_pair(selected_qpoints, (end - begin).norm());
|
||||
};
|
||||
|
||||
// 根据搜索到的 q 点, 计算图中每个点的值
|
||||
auto calculate_values = []
|
||||
(
|
||||
// search_qpoints 的第一个返回值
|
||||
const boost::container::flat_map<double, std::size_t>& path,
|
||||
// 每一条连续路径的第一个 q 点的索引
|
||||
const std::set<std::size_t>& path_begin,
|
||||
// 所有 q 点的数据（需要用到它的频率和权重）
|
||||
const std::vector<UnfoldOutput::QpointDataType>& qpoints,
|
||||
// 用于插值的分辨率和范围
|
||||
const std::array<std::size_t, 2>& resolution,
|
||||
const std::array<double, 2>& frequency_range,
|
||||
// 路径的总长度
|
||||
double total_distance
|
||||
)
|
||||
{
|
||||
// 按比例混合两个 q 点的结果，得到可以用于画图的那一列数据
|
||||
auto blend = [&]
|
||||
(
|
||||
// 两个点的索引
|
||||
std::size_t a, std::size_t b,
|
||||
// 按照连续路径混合还是按照断开的路径混合
|
||||
bool continuous,
|
||||
// 第一个点占的比例
|
||||
double ratio,
|
||||
std::size_t resolution, std::array<double, 2> frequency_range
|
||||
) -> std::vector<double>
|
||||
{
|
||||
// 混合得到的频率和权重
|
||||
std::vector<double> frequency, weight;
|
||||
// 如果是连续路径，将每个模式的频率和权重按照比例混合
|
||||
if (continuous)
|
||||
{
|
||||
assert(qpoints[a].ModeData.size() == qpoints[b].ModeData.size());
|
||||
for (std::size_t i = 0; i < qpoints[a].ModeData.size(); i++)
|
||||
{
|
||||
frequency.push_back
|
||||
(qpoints[a].ModeData[i].Frequency * ratio + qpoints[b].ModeData[i].Frequency * (1 - ratio));
|
||||
weight.push_back(qpoints[a].ModeData[i].Weight * ratio + qpoints[b].ModeData[i].Weight * (1 - ratio));
|
||||
}
|
||||
}
|
||||
// 如果是不连续路径，将每个模式的权重乘以比例，最后相加
|
||||
else
|
||||
{
|
||||
for (std::size_t i = 0; i < qpoints[a].ModeData.size(); i++)
|
||||
{
|
||||
frequency.push_back(qpoints[a].ModeData[i].Frequency);
|
||||
weight.push_back(qpoints[a].ModeData[i].Weight * ratio);
|
||||
}
|
||||
for (std::size_t i = 0; i < qpoints[b].ModeData.size(); i++)
|
||||
{
|
||||
frequency.push_back(qpoints[b].ModeData[i].Frequency);
|
||||
weight.push_back(qpoints[b].ModeData[i].Weight * (1 - ratio));
|
||||
}
|
||||
}
|
||||
std::vector<double> result(resolution);
|
||||
for (std::size_t i = 0; i < frequency.size(); i++)
|
||||
{
|
||||
std::ptrdiff_t index = (frequency[i] - frequency_range[0]) / (frequency_range[1] - frequency_range[0])
|
||||
* resolution;
|
||||
if (index >= 0 && index < static_cast<std::ptrdiff_t>(resolution)) result[index] += weight[i];
|
||||
}
|
||||
return result;
|
||||
};
|
||||
|
||||
std::vector<std::vector<double>> values;
|
||||
for (std::size_t i = 0; i < resolution[0]; i++)
|
||||
{
|
||||
auto current_distance = total_distance * i / resolution[0];
|
||||
auto it = path.lower_bound(current_distance);
|
||||
if (it == path.begin()) values.push_back(blend
|
||||
(it->second, it->second, true, 1, resolution[1], frequency_range));
|
||||
else if (it == path.end()) values.push_back(blend
|
||||
(
|
||||
std::prev(it)->second, std::prev(it)->second, true, 1,
|
||||
resolution[1], frequency_range
|
||||
));
|
||||
else values.push_back(blend
|
||||
(
|
||||
std::prev(it)->second, it->second, !path_begin.contains(it->second),
|
||||
(it->first - current_distance) / (it->first - std::prev(it)->first),
|
||||
resolution[1], frequency_range
|
||||
));
|
||||
}
|
||||
return values;
|
||||
};
|
||||
|
||||
// 根据数值, 画图
|
||||
auto plot = []
|
||||
(
|
||||
const std::vector<std::vector<double>>& values,
|
||||
const std::string& filename,
|
||||
const std::vector<double>& x_ticks, const std::vector<double>& y_ticks,
|
||||
const std::vector<std::string>& y_ticklabels,
|
||||
const std::optional<double>& aspect_ratio,
|
||||
const std::optional<std::array<std::size_t, 2>>& resolution
|
||||
)
|
||||
{
|
||||
std::vector<std::vector<double>>
|
||||
r(values[0].size(), std::vector<double>(values.size(), 0)),
|
||||
g(values[0].size(), std::vector<double>(values.size(), 0)),
|
||||
b(values[0].size(), std::vector<double>(values.size(), 0)),
|
||||
a(values[0].size(), std::vector<double>(values.size(), 0));
|
||||
for (std::size_t i = 0; i < values[0].size(); i++)
|
||||
for (std::size_t j = 0; j < values.size(); j++)
|
||||
{
|
||||
auto v = values[j][i];
|
||||
if (v < 0.05) v = 0;
|
||||
a[i][j] = v * 100 * 255;
|
||||
if (a[i][j] > 255) a[i][j] = 255;
|
||||
r[i][j] = 255 - v * 2 * 255;
|
||||
if (r[i][j] < 0) r[i][j] = 0;
|
||||
g[i][j] = 255 - v * 2 * 255;
|
||||
if (g[i][j] < 0) g[i][j] = 0;
|
||||
b[i][j] = 255;
|
||||
}
|
||||
auto f = matplot::figure(true);
|
||||
auto ax = f->current_axes();
|
||||
auto image = ax->image(std::tie(r, g, b));
|
||||
image->matrix_a(a);
|
||||
ax->y_axis().reverse(false);
|
||||
ax->x_axis().tick_values(x_ticks);
|
||||
ax->x_axis().tick_length(1);
|
||||
ax->x_axis().ticklabels(std::vector<std::string>(x_ticks.size()));
|
||||
ax->y_axis().tick_values(y_ticks);
|
||||
ax->y_axis().tick_length(1);
|
||||
ax->y_axis().ticklabels(y_ticklabels);
|
||||
if (aspect_ratio)
|
||||
{
|
||||
ax->axes_aspect_ratio_auto(false);
|
||||
ax->axes_aspect_ratio(*aspect_ratio);
|
||||
}
|
||||
if (resolution)
|
||||
{
|
||||
f->width((*resolution)[0]);
|
||||
f->height((*resolution)[1]);
|
||||
}
|
||||
f->save(filename, "png");
|
||||
};
|
||||
|
||||
auto input = YAML::LoadFile(config_file).as<Input>();
|
||||
auto unfolded_data = biu::deserialize<UnfoldOutput>
|
||||
(biu::read<std::byte>(input.UnfoldedDataFile));
|
||||
|
||||
// 搜索画图需要用到的 q 点
|
||||
// key 到起点的距离，value 为 q 点在 QpointData 中的索引
|
||||
boost::container::flat_map<double, std::size_t> path;
|
||||
// 每一条连续路径的第一个 q 点在 path 中的索引
|
||||
std::set<std::size_t> path_begin;
|
||||
// x 轴的刻度，为 path 中的索引
|
||||
std::set<std::size_t> x_ticks_index;
|
||||
double total_distance = 0;
|
||||
for (auto& line : input.Qpoints)
|
||||
{
|
||||
assert(line.size() >= 2);
|
||||
path_begin.insert(path.size());
|
||||
for (std::size_t i = 0; i < line.size() - 1; i++)
|
||||
{
|
||||
x_ticks_index.insert(path.size());
|
||||
auto [this_path, this_distance] = search_qpoints
|
||||
(
|
||||
unfolded_data.PrimativeCell, {line[i], line[i + 1]},
|
||||
unfolded_data.QpointData
|
||||
| ranges::views::transform(&UnfoldOutput::QpointDataType::Qpoint)
|
||||
| ranges::to_vector,
|
||||
input.ThresholdWhenSearchingQpoints.value_or(0.001),
|
||||
i != line.size() - 2
|
||||
);
|
||||
path.merge
|
||||
(
|
||||
this_path
|
||||
| ranges::views::transform([&](auto& p)
|
||||
{ return std::make_pair(p.first + total_distance, p.second); })
|
||||
| ranges::to<boost::container::flat_map>
|
||||
);
|
||||
total_distance += this_distance;
|
||||
}
|
||||
}
|
||||
|
||||
// 计算画图的数据
|
||||
auto values = calculate_values
|
||||
(
|
||||
path, path_begin, unfolded_data.QpointData, input.InterpolationResolution,
|
||||
input.FrequencyRange, total_distance
|
||||
);
|
||||
auto x_ticks = x_ticks_index | ranges::views::transform([&](auto i)
|
||||
{ return path.nth(i)->first / total_distance * input.InterpolationResolution[0]; }) | ranges::to<std::vector>;
|
||||
auto y_ticks = input.YTicks.value_or(std::vector<std::pair<double, std::string>>{})
|
||||
| biu::toLvalue | ranges::views::keys
|
||||
| ranges::views::transform([&](auto i)
|
||||
{
|
||||
return (i - input.FrequencyRange[0]) / (input.FrequencyRange[1] - input.FrequencyRange[0])
|
||||
* input.InterpolationResolution[1];
|
||||
})
|
||||
| ranges::to_vector;
|
||||
auto y_ticklabels = input.YTicks.value_or(std::vector<std::pair<double, std::string>>{})
|
||||
| biu::toLvalue | ranges::views::values | ranges::to_vector;
|
||||
if (input.OutputPictureFile) plot
|
||||
(
|
||||
values, input.OutputPictureFile.value(),
|
||||
x_ticks, y_ticks, y_ticklabels, input.AspectRatio, input.PictureResolution
|
||||
);
|
||||
if (input.OutputDataFile)
|
||||
biu::Hdf5file(input.OutputDataFile.value(), true)
|
||||
.write("Values", values)
|
||||
.write("XTicks", x_ticks)
|
||||
.write("YTicks", y_ticks)
|
||||
.write("YTickLabels", y_ticklabels)
|
||||
.write("InterpolationResolution", input.InterpolationResolution)
|
||||
.write("FrequencyRange", input.FrequencyRange);
|
||||
}
|
||||
|
||||
void ufo::plot_point(std::string config_file)
|
||||
{
|
||||
struct Input
|
||||
{
|
||||
std::string UnfoldedDataFile;
|
||||
// 要画图的 q 点
|
||||
Eigen::Vector3d Qpoint;
|
||||
// 插值的分辨率
|
||||
std::size_t InterpolationResolution;
|
||||
std::optional<double> AspectRatio;
|
||||
std::optional<std::array<std::size_t, 2>> PictureResolution;
|
||||
// 画图的频率范围
|
||||
std::array<double, 2> FrequencyRange;
|
||||
// 搜索 q 点时的阈值，单位为埃^-1
|
||||
std::optional<double> ThresholdWhenSearchingQpoints;
|
||||
// 是否要在 z 轴上作一些标记
|
||||
std::optional<std::vector<std::pair<double, std::string>>> XTicks;
|
||||
// 是否输出图片
|
||||
std::optional<std::string> OutputPictureFile;
|
||||
// 是否输出插值后数据，可以进一步使用 matplotlib 画图
|
||||
std::optional<std::string> OutputDataFile;
|
||||
// 是否输出插值前数据，可以配合 phonopy 结果深入研究
|
||||
std::optional<std::string> OutputRawDataFile;
|
||||
};
|
||||
|
||||
// 根据 q 点路径, 搜索要使用的 q 点，返回的是 q 点在 QpointData 中的索引
|
||||
auto search_qpoints = []
|
||||
(
|
||||
const Eigen::Matrix3d& primative_cell,
|
||||
const Eigen::Vector3d& qpoint, const std::vector<Eigen::Vector3d>& qpoints,
|
||||
double threshold
|
||||
)
|
||||
{
|
||||
biu::Logger::Guard log(qpoint);
|
||||
// 对于 output 中的每一个点, 检查这个点是否与所寻找的点足够近，如果足够近则返回
|
||||
for (std::size_t i = 0; i < qpoints.size(); i++)
|
||||
for (auto cell_shift
|
||||
: biu::sequence(Eigen::Vector3i(-1, -1, -1), Eigen::Vector3i(2, 2, 2)))
|
||||
{
|
||||
auto this_qpoint
|
||||
= (primative_cell.reverse().transpose() * (qpoints[i] + cell_shift.first.cast<double>())).eval();
|
||||
if ((this_qpoint - primative_cell.reverse().transpose() * qpoint).norm() < threshold) return log.rtn(i);
|
||||
}
|
||||
throw std::runtime_error("No q points found");
|
||||
};
|
||||
|
||||
// 根据搜索到的 q 点, 计算图中每个点的值
|
||||
auto calculate_values = []
|
||||
(
|
||||
// q 点的数据（需要用到它的频率和权重）
|
||||
const UnfoldOutput::QpointDataType& qpoint,
|
||||
// 用于插值的分辨率和范围
|
||||
std::size_t resolution,
|
||||
const std::array<double, 2>& frequency_range
|
||||
)
|
||||
{
|
||||
biu::Logger::Guard log;
|
||||
std::vector<double> result(resolution);
|
||||
for (auto& mode : qpoint.ModeData)
|
||||
{
|
||||
double index_double = (mode.Frequency - frequency_range[0]) / (frequency_range[1] - frequency_range[0])
|
||||
* (resolution - 1);
|
||||
std::ptrdiff_t index = std::round(index_double);
|
||||
if (index >= 0 && index < static_cast<std::ptrdiff_t>(resolution)) result[index] += mode.Weight;
|
||||
}
|
||||
return log.rtn(result);
|
||||
};
|
||||
|
||||
// 根据数值, 画图
|
||||
auto plot = []
|
||||
(
|
||||
const std::vector<double>& values, const std::string& filename,
|
||||
const std::vector<double>& x_ticks, const std::vector<std::string>& x_ticklabels,
|
||||
const std::optional<double>& aspect_ratio, const std::optional<std::array<std::size_t, 2>>& resolution
|
||||
)
|
||||
{
|
||||
biu::Logger::Guard log;
|
||||
auto f = matplot::figure(true);
|
||||
auto ax = f->current_axes();
|
||||
auto image = ax->area(values, 0, false, "");
|
||||
ax->y_axis().reverse(false);
|
||||
ax->x_axis().tick_values(x_ticks);
|
||||
ax->x_axis().tick_length(1);
|
||||
ax->x_axis().ticklabels(x_ticklabels);
|
||||
ax->y_axis().tick_values({});
|
||||
if (aspect_ratio)
|
||||
{
|
||||
ax->axes_aspect_ratio_auto(false);
|
||||
ax->axes_aspect_ratio(*aspect_ratio);
|
||||
}
|
||||
if (resolution)
|
||||
{
|
||||
f->width((*resolution)[0]);
|
||||
f->height((*resolution)[1]);
|
||||
}
|
||||
f->save(filename, "png");
|
||||
};
|
||||
|
||||
biu::Logger::Guard log;
|
||||
auto input = YAML::LoadFile(config_file).as<Input>();
|
||||
auto unfolded_data = biu::deserialize<UnfoldOutput>
|
||||
(biu::read<std::byte>(input.UnfoldedDataFile));
|
||||
|
||||
auto qpoint_index = search_qpoints
|
||||
(
|
||||
unfolded_data.PrimativeCell, input.Qpoint,
|
||||
unfolded_data.QpointData
|
||||
| ranges::views::transform(&UnfoldOutput::QpointDataType::Qpoint)
|
||||
| ranges::to_vector,
|
||||
input.ThresholdWhenSearchingQpoints.value_or(0.001)
|
||||
);
|
||||
auto values = calculate_values
|
||||
(
|
||||
unfolded_data.QpointData[qpoint_index],
|
||||
input.InterpolationResolution, input.FrequencyRange
|
||||
);
|
||||
auto x_ticks = input.XTicks.value_or(std::vector<std::pair<double, std::string>>{})
|
||||
| biu::toLvalue | ranges::views::keys
|
||||
| ranges::views::transform([&](auto i)
|
||||
{
|
||||
return (i - input.FrequencyRange[0]) / (input.FrequencyRange[1] - input.FrequencyRange[0])
|
||||
* input.InterpolationResolution;
|
||||
})
|
||||
| ranges::to_vector;
|
||||
auto x_ticklabels = input.XTicks.value_or(std::vector<std::pair<double, std::string>>{})
|
||||
| biu::toLvalue | ranges::views::values | ranges::to_vector;
|
||||
if (input.OutputPictureFile) plot
|
||||
(
|
||||
values, input.OutputPictureFile.value(),
|
||||
x_ticks, x_ticklabels, input.AspectRatio, input.PictureResolution
|
||||
);
|
||||
if (input.OutputDataFile)
|
||||
biu::Hdf5file(input.OutputDataFile.value(), true)
|
||||
.write("Values", values)
|
||||
.write("XTicks", x_ticks)
|
||||
.write("XTickLabels", x_ticklabels)
|
||||
.write("InterpolationResolution", input.InterpolationResolution)
|
||||
.write("FrequencyRange", input.FrequencyRange);
|
||||
if (input.OutputRawDataFile)
|
||||
std::ofstream(*input.OutputRawDataFile) << YAML::Node(unfolded_data.QpointData[qpoint_index]);
|
||||
}
|
||||
@@ -1,138 +0,0 @@
|
||||
# include <ufo.hpp>
|
||||
|
||||
void ufo::raman_create_displacement(std::string config_file)
|
||||
{
|
||||
struct Input
|
||||
{
|
||||
std::string UnfoldedDataFile;
|
||||
// 搜索位于 Gamma 点的 q 点时，使用的阈值，单位为埃^-1，默认为 0.01
|
||||
std::optional<double> ThresholdWhenSearchingQpoints;
|
||||
// 搜索权重非零的模式时，使用的阈值，默认为 0.01
|
||||
std::optional<double> ThresholdWhenSearchingModes;
|
||||
// 所有原子的符号
|
||||
std::vector<std::string> AtomSymbols;
|
||||
// 各种原子的质量，单位为原子质量
|
||||
std::map<std::string, double> AtomMasses;
|
||||
// 原子最大位移大小，单位为埃
|
||||
double MaxDisplacement;
|
||||
// 超胞，单位为埃
|
||||
Eigen::Matrix3d SuperCell;
|
||||
// 超胞中各个原子的坐标，单位为超胞的格矢
|
||||
Eigen::MatrixX3d AtomPositions;
|
||||
// 输出的POSCAR所在的目录
|
||||
std::string OutputPoscarDirectory;
|
||||
// 输出的数据文件名
|
||||
std::string OutputDataFile;
|
||||
};
|
||||
struct Output
|
||||
{
|
||||
struct ModeData_t
|
||||
{
|
||||
std::size_t MetaQpointIndex;
|
||||
std::size_t ModeIndex;
|
||||
// 每个原子的位移，单位为埃
|
||||
Eigen::MatrixX3d AtomMovement;
|
||||
};
|
||||
std::vector<ModeData_t> ModeData;
|
||||
using serialize = zpp::bits::members<1>;
|
||||
};
|
||||
|
||||
// 假定同类型的原子一定写在一起
|
||||
auto generate_poscar = []
|
||||
(
|
||||
Eigen::Matrix3d SuperCell, Eigen::MatrixX3d AtomPositions, std::vector<std::string> AtomSymbols
|
||||
)
|
||||
{
|
||||
std::stringstream ss;
|
||||
ss << "some random comment to make VASP happy\n1.0\n";
|
||||
for (std::size_t i = 0; i < 3; i++)
|
||||
{
|
||||
for (std::size_t j = 0; j < 3; j++) ss << SuperCell(i, j) << " ";
|
||||
ss << std::endl;
|
||||
}
|
||||
auto atom_symbols = AtomSymbols | ranges::views::chunk_by(std::ranges::equal_to{});
|
||||
ss << "{}\n"_f(ranges::accumulate
|
||||
(
|
||||
atom_symbols | ranges::views::transform([](auto&& chunk) { return chunk[0]; }),
|
||||
""s, [](auto&& a, auto&& b) { return a + " " + b; }
|
||||
));
|
||||
ss << "{}\n"_f(ranges::accumulate
|
||||
(
|
||||
atom_symbols | ranges::views::transform([](auto&& chunk) { return chunk.size(); }),
|
||||
""s, [](auto&& a, auto&& b) { return a + " " + std::to_string(b); }
|
||||
));
|
||||
ss << "Direct\n";
|
||||
for (const auto& position : AtomPositions.rowwise())
|
||||
{
|
||||
for (std::size_t i = 0; i < 3; i++) ss << position(i) << " ";
|
||||
ss << std::endl;
|
||||
}
|
||||
return ss.str();
|
||||
};
|
||||
|
||||
auto input = YAML::LoadFile(config_file).as<Input>();
|
||||
auto unfolded_data = biu::deserialize<UnfoldOutput>
|
||||
(biu::read<std::byte>(input.UnfoldedDataFile));
|
||||
Output output;
|
||||
|
||||
// 搜索满足条件的模式，找到满足条件的模式后，就将 MetaQpoint 的索引加入到 output 中
|
||||
// 之所以使用 MetaQpoint 的索引而不是 Qpoint 的索引，是因为 Qpoint 中可能有指向同一个 MetaQpoint 中模式的不同模式都满足要求
|
||||
// 如果写入 Qpoint 的索引，就会重复而增加之后的计算量
|
||||
std::set<std::pair<std::size_t, std::size_t>> selected_modes;
|
||||
for (const auto& qpoint : unfolded_data.QpointData)
|
||||
{
|
||||
if
|
||||
(
|
||||
(unfolded_data.PrimativeCell.reverse().transpose() * qpoint.Qpoint).norm()
|
||||
> input.ThresholdWhenSearchingQpoints.value_or(0.01)
|
||||
)
|
||||
continue;
|
||||
for (std::size_t i = 0; i < qpoint.ModeData.size(); i++)
|
||||
{
|
||||
if (qpoint.ModeData[i].Weight < input.ThresholdWhenSearchingModes.value_or(0.01)) continue;
|
||||
selected_modes.insert({qpoint.SourceIndex, i});
|
||||
}
|
||||
}
|
||||
|
||||
// 构造输出数据
|
||||
for (auto [i, j] : selected_modes)
|
||||
{
|
||||
auto& mode_data = output.ModeData.emplace_back();
|
||||
mode_data.MetaQpointIndex = i;
|
||||
mode_data.ModeIndex = j;
|
||||
// 未归一化的位移, 假定虚部总是为零
|
||||
auto atom_movement =
|
||||
unfolded_data.MetaQpointData[i].ModeData[j].AtomMovement.real().cwiseProduct
|
||||
(
|
||||
(
|
||||
input.AtomSymbols
|
||||
| ranges::views::transform([&](const auto& symbol)
|
||||
{ return input.AtomMasses.at(symbol); })
|
||||
| ranges::to_vector
|
||||
| biu::toEigen<>
|
||||
).cwiseSqrt().cwiseInverse().rowwise().replicate(3)
|
||||
).eval();
|
||||
// 归一化
|
||||
mode_data.AtomMovement = atom_movement / atom_movement.rowwise().norm().maxCoeff() * input.MaxDisplacement;
|
||||
}
|
||||
|
||||
// 输出
|
||||
std::ofstream(input.OutputDataFile, std::ios::binary) << biu::serialize<char>(output);
|
||||
for (std::size_t i = 0; i < output.ModeData.size(); i++)
|
||||
{
|
||||
std::filesystem::create_directories(input.OutputPoscarDirectory + "/" + std::to_string(i));
|
||||
std::ofstream(input.OutputPoscarDirectory + "/" + std::to_string(i) + "/POSCAR") << generate_poscar
|
||||
(
|
||||
input.SuperCell,
|
||||
input.AtomPositions + output.ModeData[i].AtomMovement * input.SuperCell.inverse(),
|
||||
input.AtomSymbols
|
||||
);
|
||||
}
|
||||
std::filesystem::create_directories(input.OutputPoscarDirectory + "/original");
|
||||
std::ofstream(input.OutputPoscarDirectory + "/original/POSCAR") << generate_poscar
|
||||
(
|
||||
input.SuperCell,
|
||||
input.AtomPositions,
|
||||
input.AtomSymbols
|
||||
);
|
||||
}
|
||||
@@ -1,345 +0,0 @@
|
||||
# include <ufo.hpp>
|
||||
# include <thread>
|
||||
# include <syncstream>
|
||||
# include <execution>
|
||||
|
||||
void ufo::unfold(std::string config_file)
|
||||
{
|
||||
// 反折叠的原理: 将超胞中的原子运动状态, 投影到一组平面波构成的基矢中.
|
||||
// 每一个平面波的波矢由两部分相加得到: 一部分是单胞倒格子的整数倍, 所取的个数有一定任意性, 论文中建议取大约单胞中原子个数那么多个;
|
||||
// 对于没有缺陷的情况, 取一个应该就足够了.
|
||||
// 这些平面波以原胞为周期。
|
||||
// 另一部分是超胞倒格子的整数倍, 取 n 个, n 为超胞对应的单胞的倍数, 其实也就是倒空间中单胞对应倒格子中超胞的格点.
|
||||
// 只要第一部分取得足够多, 那么单胞中原子的状态就可以完全被这些平面波描述.
|
||||
// 将超胞中原子的运动状态投影到这些基矢上, 计算出投影的系数, 就可以将超胞的原子运动状态分解到单胞中的多个 q 点上.
|
||||
|
||||
struct Input
|
||||
{
|
||||
// 单胞的三个格矢，每行表示一个格矢的坐标，单位为埃
|
||||
Eigen::Matrix3d PrimativeCell;
|
||||
|
||||
// 单胞到超胞的格矢转换时用到的矩阵
|
||||
// SuperCellMultiplier 是一个三维列向量且各个元素都是整数，表示单胞在各个方向扩大到多少倍之后，可以得到和超胞一样的体积
|
||||
// SuperCellDeformation 是一个行列式为 1 的矩阵，它表示经过 SuperCellMultiplier 扩大后，还需要怎样的变换才能得到超胞
|
||||
// SuperCell = (SuperCellDeformation * SuperCellMultiplier.asDiagonal()) * PrimativeCell
|
||||
// ReciprocalPrimativeCell = (SuperCellDeformation * SuperCellMultiplier.asDiagonal()).transpose()
|
||||
// * ReciprocalSuperCell
|
||||
// Position = PositionToCell(line vector) * Cell
|
||||
// InversePosition = InversePositionToCell(line vector) * ReciprocalCell
|
||||
// PositionToSuperCell(line vector) * SuperCell = PositionToPrimativeCell(line vector) * PrimativeCell
|
||||
// ReciprocalPositionToSuperCell(line vector) * ReciprocalSuperCell
|
||||
// = ReciprocalPositionToPrimativeCell(line vector) * ReciprocalPrimativeCell
|
||||
Eigen::Matrix3d SuperCellDeformation;
|
||||
Eigen::Vector3i SuperCellMultiplier;
|
||||
|
||||
// 在单胞内取几个平面波的基矢
|
||||
Eigen::Vector<std::size_t, 3> PrimativeCellBasisNumber;
|
||||
|
||||
// 超胞中原子的坐标，每行表示一个原子的坐标，单位为超胞的格矢
|
||||
Eigen::MatrixX3d AtomPositionBySuperCell;
|
||||
|
||||
// 从 band.hdf5 读入 QpointData
|
||||
std::optional<std::string> QpointDataInputFile;
|
||||
|
||||
// 输出到哪些文件
|
||||
struct QpointDataOutputFileType
|
||||
{
|
||||
std::string Filename;
|
||||
|
||||
// 如果指定，则将结果投影到那些原子上
|
||||
std::optional<std::vector<std::size_t>> SelectedAtoms;
|
||||
|
||||
// 默认输出为 zpp 文件，如果指定为 true，则输出为 yaml 文件
|
||||
std::optional<bool> OutputAsYaml;
|
||||
};
|
||||
std::vector<QpointDataOutputFileType> QpointDataOutputFile;
|
||||
};
|
||||
|
||||
// 从文件中读取 QpointData
|
||||
auto read_qpoint_data = [](std::string filename)
|
||||
{
|
||||
// 读入原始数据
|
||||
// phonopy 的输出有两种可能
|
||||
// 直接指定计算的 q 点时，frequency 是 2 维，这时第一个维度是 q 点，第二个维度是不同模式
|
||||
// 计算能带时，frequency 是 3 维，相比于二维的情况多了第一个维度，表示 q 点所在路径
|
||||
// qpoint 或 path，以及 eigenvector 也有类似的变化
|
||||
// eigenvector 是三维或四维的数组，后两个维度分别表示原子运动和模式（而不是模式和原子），
|
||||
// 因为后两个维度的尺寸总是一样的（模式个数等于原子坐标个数），非常容易搞错
|
||||
std::vector<std::array<double, 3>> qpoint;
|
||||
std::vector<std::vector<double>> frequency;
|
||||
std::vector<std::vector<std::vector<biu::PhonopyComplex>>> eigenvector_vector;
|
||||
auto file = biu::Hdf5file(filename);
|
||||
|
||||
if (file.File.getDataSet("/frequency").getDimensions().size() == 2)
|
||||
file.read("/frequency", frequency)
|
||||
.read("/eigenvector", eigenvector_vector)
|
||||
.read("/qpoint", qpoint);
|
||||
else
|
||||
{
|
||||
std::vector<std::vector<std::array<double, 3>>> temp_path;
|
||||
std::vector<std::vector<std::vector<double>>> temp_frequency;
|
||||
std::vector<std::vector<std::vector<std::vector<biu::PhonopyComplex>>>> temp_eigenvector_vector;
|
||||
file.read("/frequency", temp_frequency)
|
||||
.read("/eigenvector", temp_eigenvector_vector)
|
||||
.read("/path", temp_path);
|
||||
frequency = temp_frequency | ranges::views::join | ranges::to_vector;
|
||||
qpoint = temp_path | ranges::views::join | ranges::to_vector;
|
||||
eigenvector_vector = temp_eigenvector_vector | ranges::views::join | ranges::to_vector;
|
||||
}
|
||||
|
||||
// 整理得到结果
|
||||
auto number_of_qpoints = frequency.size(), num_of_modes = frequency[0].size();
|
||||
std::vector<UnfoldOutput::MetaQpointDataType> qpoint_data(number_of_qpoints);
|
||||
for (std::size_t i = 0; i < number_of_qpoints; i++)
|
||||
{
|
||||
qpoint_data[i].Qpoint = qpoint[i] | biu::toEigen<>;
|
||||
qpoint_data[i].ModeData.resize(num_of_modes);
|
||||
for (std::size_t j = 0; j < num_of_modes; j++)
|
||||
{
|
||||
qpoint_data[i].ModeData[j].Frequency = frequency[i][j];
|
||||
auto number_of_atoms = eigenvector_vector[i].size() / 3;
|
||||
Eigen::MatrixX3cd eigenvectors(number_of_atoms, 3);
|
||||
for (std::size_t k = 0; k < number_of_atoms; k++) for (std::size_t l = 0; l < 3; l++)
|
||||
eigenvectors(k, l)
|
||||
= eigenvector_vector[i][k * 3 + l][j].r + eigenvector_vector[i][k * 3 + l][j].i * 1i;
|
||||
// 原则上讲，需要对读入的原子运动状态作相位转换, 使得它们与我们的约定一致(对超胞周期性重复)，但这个转换 phonopy 已经做了
|
||||
// 这里还要需要做归一化处理 (指将数据简单地作为向量处理的归一化)
|
||||
qpoint_data[i].ModeData[j].AtomMovement = eigenvectors / eigenvectors.norm();
|
||||
}
|
||||
}
|
||||
return qpoint_data;
|
||||
};
|
||||
|
||||
// 构建基
|
||||
// 每个 q 点对应一组 sub qpoint。不同的 q 点所对应的 sub qpoint 是不一样的，但 sub qpoint 与 q 点的相对位移在不同 q 点之间是相同的。
|
||||
// 由于基只与这个相对位置有关（也就是说，不同 q 点的基是一样的），因此可以先计算出所有的基，这样降低计算量。
|
||||
// 外层下标对应超胞倒格子的整数倍那部分(第二部分), 也就是不同的 sub qpoint
|
||||
// 内层下标对应单胞倒格子的整数倍那部分(第一部分), 也就是 sub qpoint 上的不同平面波（取的数量越多，结果越精确）
|
||||
auto construct_basis = []
|
||||
(
|
||||
Eigen::Matrix3d primative_cell, Eigen::Vector3i super_cell_multiplier,
|
||||
Eigen::Vector<std::size_t, 3> primative_cell_basis_number, Eigen::MatrixX3d atom_position
|
||||
)
|
||||
{
|
||||
biu::Logger::Guard log;
|
||||
std::vector<std::vector<Eigen::VectorXcd>> basis(super_cell_multiplier.prod());
|
||||
// diff_of_sub_qpoint 表示 sub qpoint 与 qpoint 的相对位置，单位为超胞的倒格矢
|
||||
for (auto [diff_of_sub_qpoint_by_reciprocal_modified_super_cell, i_of_sub_qpoint]
|
||||
: biu::sequence(super_cell_multiplier))
|
||||
{
|
||||
basis[i_of_sub_qpoint].resize(primative_cell_basis_number.prod());
|
||||
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<double>()
|
||||
+ super_cell_multiplier.cast<double>().cwiseInverse().asDiagonal()
|
||||
* diff_of_sub_qpoint_by_reciprocal_modified_super_cell.cast<double>();
|
||||
// 将单位转换为埃^-1
|
||||
auto diff_of_sub_qpoint = (diff_of_sub_qpoint_by_reciprocal_primative_cell.transpose()
|
||||
* (primative_cell.transpose().inverse())).transpose();
|
||||
// 计算基矢
|
||||
basis[i_of_sub_qpoint][i_of_basis]
|
||||
= (2i * std::numbers::pi_v<double> * (atom_position * diff_of_sub_qpoint)).array().exp();
|
||||
}
|
||||
}
|
||||
return basis;
|
||||
};
|
||||
|
||||
// 计算从超胞到原胞的投影系数（不是分原子的投影系数），是反折叠的核心步骤
|
||||
// 返回的投影系数是一个三维数组，第一维对应不同的 q 点，第二维对应不同的模式，第三维对应不同的 sub qpoint
|
||||
auto construct_projection_coefficient = []
|
||||
(
|
||||
const std::vector<std::vector<Eigen::VectorXcd>>& basis,
|
||||
// 实际上只需要其中的 AtomMovement
|
||||
const std::vector<UnfoldOutput::MetaQpointDataType>& qpoint_data,
|
||||
std::atomic<std::size_t>& number_of_finished_modes
|
||||
)
|
||||
{
|
||||
// 将所有的模式取出，组成一个一维数组，稍后并行计算
|
||||
std::vector<std::reference_wrapper<const Eigen::MatrixX3cd>> mode_data;
|
||||
for (auto& qpoint : qpoint_data) for (auto& mode : qpoint.ModeData)
|
||||
mode_data.emplace_back(mode.AtomMovement);
|
||||
// 第一层下标对应不同模式, 第二层下标对应这个模式在反折叠后的 q 点(sub qpoint)
|
||||
std::vector<std::vector<double>> projection_coefficient(mode_data.size());
|
||||
// 对每个模式并行
|
||||
std::transform
|
||||
(
|
||||
std::execution::par_unseq, mode_data.begin(), mode_data.end(),
|
||||
projection_coefficient.begin(), [&](const auto& mode_data)
|
||||
{
|
||||
// 这里, mode_data 和 projection_coefficient 均指对应于一个模式的数据
|
||||
std::vector<double> projection_coefficient(basis.size());
|
||||
for (std::size_t i_of_sub_qpoint = 0; i_of_sub_qpoint < basis.size(); i_of_sub_qpoint++)
|
||||
// 对于 basis 中, 对应于单胞倒格子的部分, 以及对应于不同方向的部分, 分别求内积, 然后求模方和
|
||||
for (std::size_t i_of_basis = 0; i_of_basis < basis[i_of_sub_qpoint].size(); i_of_basis++)
|
||||
projection_coefficient[i_of_sub_qpoint] +=
|
||||
(basis[i_of_sub_qpoint][i_of_basis].transpose().conjugate() * mode_data.get())
|
||||
.array().abs2().sum();
|
||||
// 如果是严格地将向量分解到一组完备的基矢上, 那么不需要对计算得到的权重再做归一化处理
|
||||
// 但这里并不是这样一个严格的概念. 因此对分解到各个 sub qpoint 上的权重做归一化处理
|
||||
auto sum = ranges::accumulate(projection_coefficient, 0.);
|
||||
for (auto& _ : projection_coefficient) _ /= sum;
|
||||
number_of_finished_modes++;
|
||||
return projection_coefficient;
|
||||
}
|
||||
);
|
||||
// 将计算得到的投影系数重新组装成三维数组
|
||||
// 第一维是 meta qpoint，第二维是模式，第三维是 sub qpoint
|
||||
std::vector<std::vector<std::vector<double>>> projection_coefficient_output;
|
||||
for
|
||||
(
|
||||
std::size_t i_of_meta_qpoint = 0, num_of_mode_manipulated = 0;
|
||||
i_of_meta_qpoint < qpoint_data.size();
|
||||
i_of_meta_qpoint++, num_of_mode_manipulated += qpoint_data[i_of_meta_qpoint].ModeData.size()
|
||||
)
|
||||
projection_coefficient_output.emplace_back
|
||||
(
|
||||
projection_coefficient.begin() + num_of_mode_manipulated,
|
||||
projection_coefficient.begin() + num_of_mode_manipulated + qpoint_data[i_of_meta_qpoint].ModeData.size()
|
||||
);
|
||||
return projection_coefficient_output;
|
||||
};
|
||||
|
||||
// 组装输出，即将投影系数应用到原始数据上
|
||||
auto construct_output = []
|
||||
(
|
||||
const Input& input,
|
||||
const std::vector<std::vector<std::vector<double>>>& projection_coefficient,
|
||||
const std::vector<UnfoldOutput::MetaQpointDataType>& qpoint_data,
|
||||
const std::optional<std::vector<std::size_t>>& selected_atoms
|
||||
)
|
||||
{
|
||||
UnfoldOutput output;
|
||||
output.PrimativeCell = input.PrimativeCell;
|
||||
output.SuperCellMultiplier = input.SuperCellMultiplier;
|
||||
output.SuperCellDeformation = input.SuperCellDeformation;
|
||||
output.SelectedAtoms = selected_atoms;
|
||||
output.MetaQpointData = qpoint_data;
|
||||
for (std::size_t i_of_meta_qpoint = 0; i_of_meta_qpoint < qpoint_data.size(); i_of_meta_qpoint++)
|
||||
{
|
||||
// 如果需要投影到特定的原子上，需要先计算当前 meta qpoint 的不同模式的投影系数
|
||||
std::optional<std::vector<double>> projection_coefficient_on_atoms;
|
||||
if (selected_atoms)
|
||||
{
|
||||
projection_coefficient_on_atoms.emplace();
|
||||
for (std::size_t i_of_mode = 0; i_of_mode < qpoint_data[i_of_meta_qpoint].ModeData.size(); i_of_mode++)
|
||||
{
|
||||
projection_coefficient_on_atoms->emplace_back(0);
|
||||
for (auto atom : *selected_atoms)
|
||||
projection_coefficient_on_atoms->back()
|
||||
+= qpoint_data[i_of_meta_qpoint].ModeData[i_of_mode].AtomMovement.row(atom).array().abs2().sum();
|
||||
projection_coefficient_on_atoms->back() *=
|
||||
static_cast<double>(qpoint_data[i_of_meta_qpoint].ModeData[i_of_mode].AtomMovement.rows())
|
||||
/ selected_atoms->size();
|
||||
}
|
||||
}
|
||||
|
||||
for
|
||||
(
|
||||
auto [diff_of_sub_qpoint_by_reciprocal_modified_super_cell, i_of_sub_qpoint]
|
||||
: biu::sequence(input.SuperCellMultiplier)
|
||||
)
|
||||
{
|
||||
auto& _ = output.QpointData.emplace_back();
|
||||
/*
|
||||
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 =
|
||||
(
|
||||
input.SuperCellMultiplier.cast<double>().cwiseInverse().asDiagonal()
|
||||
* (
|
||||
diff_of_sub_qpoint_by_reciprocal_modified_super_cell.cast<double>()
|
||||
+ input.SuperCellDeformation.inverse() * qpoint_data[i_of_meta_qpoint].Qpoint
|
||||
)
|
||||
).eval();
|
||||
_.Qpoint = sub_qpoint_by_reciprocal_primative_cell.array()
|
||||
- sub_qpoint_by_reciprocal_primative_cell.array().floor();
|
||||
_.Source = qpoint_data[i_of_meta_qpoint].Qpoint;
|
||||
_.SourceIndex = i_of_meta_qpoint;
|
||||
|
||||
for (std::size_t i_of_mode = 0; i_of_mode < qpoint_data[i_of_meta_qpoint].ModeData.size(); i_of_mode++)
|
||||
{
|
||||
auto& __ = _.ModeData.emplace_back();
|
||||
__.Frequency = qpoint_data[i_of_meta_qpoint].ModeData[i_of_mode].Frequency;
|
||||
__.Weight = projection_coefficient[i_of_meta_qpoint][i_of_mode][i_of_sub_qpoint];
|
||||
if (selected_atoms)
|
||||
__.Weight *= projection_coefficient_on_atoms.value()[i_of_mode];
|
||||
}
|
||||
}
|
||||
}
|
||||
return output;
|
||||
};
|
||||
|
||||
biu::Logger::Guard log;
|
||||
log.info("Reading input file... ");
|
||||
auto input = YAML::LoadFile(config_file).as<Input>();
|
||||
auto qpoint_data = read_qpoint_data(input.QpointDataInputFile.value_or("band.hdf5"));
|
||||
log.info("Done.");
|
||||
|
||||
std::clog << "Constructing basis... " << std::flush;
|
||||
|
||||
auto basis = construct_basis
|
||||
(
|
||||
input.PrimativeCell, input.SuperCellMultiplier,
|
||||
input.PrimativeCellBasisNumber,
|
||||
input.AtomPositionBySuperCell
|
||||
* (input.SuperCellDeformation * input.SuperCellMultiplier.cast<double>().asDiagonal() * input.PrimativeCell)
|
||||
);
|
||||
std::clog << "Done." << std::endl;
|
||||
|
||||
std::clog << "Calculating projection coefficient... " << std::flush;
|
||||
// 用来在屏幕上输出进度的计数器和线程
|
||||
std::atomic<std::size_t> number_of_finished_modes(0);
|
||||
auto number_of_modes = ranges::accumulate
|
||||
(
|
||||
qpoint_data
|
||||
| ranges::views::transform([](const auto& qpoint)
|
||||
{ return qpoint.ModeData.size(); }),
|
||||
0ul
|
||||
);
|
||||
std::atomic<bool> finished;
|
||||
std::thread print_thread([&]
|
||||
{
|
||||
while (true)
|
||||
{
|
||||
std::osyncstream(std::clog)
|
||||
<< "\rCalculating projection coefficient... ({}/{})"_f(number_of_finished_modes, number_of_modes)
|
||||
<< std::flush;
|
||||
std::this_thread::sleep_for(100ms);
|
||||
if (finished) break;
|
||||
}
|
||||
});
|
||||
auto projection_coefficient = construct_projection_coefficient(basis, qpoint_data, number_of_finished_modes);
|
||||
finished = true;
|
||||
print_thread.join();
|
||||
std::clog << "\33[2K\rCalculating projection coefficient... Done." << std::endl;
|
||||
|
||||
std::clog << "Writing data... " << std::flush;
|
||||
for (auto& output_file : input.QpointDataOutputFile)
|
||||
{
|
||||
auto output = construct_output
|
||||
(input, projection_coefficient, qpoint_data, output_file.SelectedAtoms);
|
||||
if (output_file.OutputAsYaml.value_or(false)) std::ofstream(output_file.Filename) << YAML::Node(output);
|
||||
else std::ofstream(output_file.Filename, std::ios::binary) << biu::serialize<char>(output);
|
||||
}
|
||||
std::clog << "Done." << std::endl;
|
||||
}
|
||||
@@ -1,18 +0,0 @@
|
||||
SuperCellMultiplier: [3, 4, 1]
|
||||
SuperCellDeformation:
|
||||
- [ 1, 0, 0 ]
|
||||
- [ 0.6666, 1, 0 ]
|
||||
- [ 0, 0, 1 ]
|
||||
Qpoints:
|
||||
- [0, 0, 0]
|
||||
- [0.05, 0, 0]
|
||||
- [0.1, 0, 0]
|
||||
- [0.15, 0, 0]
|
||||
- [0.2, 0, 0]
|
||||
- [0.25, 0, 0]
|
||||
- [0.3, 0, 0]
|
||||
- [0.35, 0, 0]
|
||||
- [0.4, 0, 0]
|
||||
- [0.45, 0, 0]
|
||||
- [0.5, 0, 0]
|
||||
OutputFile: fold-output.yaml
|
||||
Reference in New Issue
Block a user