chn/ufo
1
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

暂存

Browse Source
This commit is contained in:
陈浩南
2023-09-11 19:54:50 +08:00
parent ec16476176
commit 78cc7ed3d3
3 changed files with 1470 additions and 144 deletions
Download Patch File Download Diff File Expand all files Collapse all files

1333
flake.lock generated
View File

File diff suppressed because it is too large Load Diff

9
flake.nix
View File

@@ -1,11 +1,12 @@
{
inputs.nixpkgs.url = "github:NixOS/nixpkgs/nixos-23.05";
inputs.nixos.url = "github:CHN-beta/nixos";
outputs = inputs: let pkgs = inputs.nixpkgs.legacyPackages.x86_64-linux; in
outputs = inputs: let pkgs = inputs.nixos.nixosConfigurations.pc.pkgs; in
{
devShell.x86_64-linux = pkgs.mkShell.override { stdenv = pkgs.gcc13Stdenv; }
devShell.x86_64-linux = pkgs.mkShell.override { stdenv = pkgs.genericPackages.gcc13Stdenv; }
{
buildInputs = with pkgs; [ yaml-cpp eigen fmt ];
buildInputs = with pkgs;
[ yaml-cpp eigen fmt (localPackages.concurrencpp.override { stdenv = genericPackages.gcc13Stdenv; }) ];
nativeBuildInputs = with pkgs; [ gdb ];
};
};

272
main.cpp
View File

@@ -3,8 +3,12 @@
# include <numbers>
# include <numeric>
# include <fstream>
# include <optional>
# include <array>
# include <utility>
# include <yaml-cpp/yaml.h>
# include <eigen3/Eigen/Dense>
# include <concurrencpp/concurrencpp.h>
using namespace std::literals;
@@ -18,10 +22,10 @@ struct Input
// 单胞的三个格矢，每行表示一个格矢的坐标，单位为埃
Eigen::Matrix3d PrimativeCell;
// 超胞在各个方向上是单胞的多少倍，这是一个对角矩阵
Eigen::Matrix3i SuperCellMultiplier;
Eigen::Matrix<unsigned, 3, 3> SuperCellMultiplier;
// 在单胞内取几个平面波的基矢
// 在 debug 阶段, 先仅取一个
Eigen::Vector3i PrimativeCellBasisNumber;
Eigen::Vector<unsigned, 3> PrimativeCellBasisNumber;
// 超胞中原子的坐标，每行表示一个原子的坐标，单位为埃
Eigen::MatrixX3d AtomPosition;
@@ -39,7 +43,7 @@ struct Input
// 模式中各个原子的运动状态
// 这个数据是这样得到的: phonopy 输出的动态矩阵的 eigenvector 乘以 $\exp(-2 \pi i \vec q \cdot \vec r)$
// 这个数据可以认为是原子位移中, 关于超胞有周期性的那一部分, 再乘以原子质量的开方.
// 这个数据在读入后会立即被归一化处理.
// 这个数据在读入后不会被立即归一化.
Eigen::MatrixX3cd AtomMovement;
};
std::vector<ModeDataType_> ModeData;
@@ -55,8 +59,6 @@ struct Input
// 各个模式的频率: phonon[*].band[*].frequency 单位为 THz 直接从 phonopy 的输出中复制
// 各个模式的原子运动状态: phonon[*].band[*].eigenvector 直接从 phonopy 的输出中复制
// 文件中可以有多余的项目, 多余的项目不管.
Input(std::string filename) : Input(YAML::LoadFile(filename)) {}
Input(YAML::Node node);
};
struct Output
@@ -74,21 +76,36 @@ struct Output
double Frequency;
// 模式的权重
double Weight;
// 来源于哪个 Q 点, 单位为超胞的倒格矢
Eigen::Vector3d Source;
};
std::vector<ModeDataType_> ModeData;
};
std::vector<QPointDataType_> QPointData;
// 将数据写入文件
void write(std::string filename) const;
};
template<> struct YAML::convert<Input> { static bool decode(const Node& node, Input& input); };
template<> struct YAML::convert<Output> { static Node encode(const Output& output); };
concurrencpp::generator<std::pair<Eigen::Vector<unsigned, 3>, unsigned>>
triplet_sequence(Eigen::Vector<unsigned, 3> range)
{
for (unsigned x = 0; x < range[0]; x++)
for (unsigned y = 0; y < range[1]; y++)
for (unsigned z = 0; z < range[2]; z++)
co_yield
{
Eigen::Vector<unsigned, 3>{{x}, {y}, {z}},
x * range[1] * range[2] + y * range[2] + z
};
}
int main(int argc, const char** argv)
{
if (argc != 3)
throw std::runtime_error("Usage: " + std::string(argv[0]) + " input.yaml output.yaml");
Input input(argv[1]);
auto input = YAML::LoadFile(argv[1]).as<Input>();
// 反折叠的原理: 将超胞中的原子运动状态, 投影到一组平面波构成的基矢中.
// 每一个平面波的波矢由两部分相加得到: 一部分是单胞倒格子的整数倍, 所取的个数有一定任意性, 论文中建议取大约单胞中原子个数那么多个;
@@ -98,176 +115,159 @@ int main(int argc, const char** argv)
// 将超胞中原子的运动状态投影到这些基矢上, 计算出投影的系数, 就可以将超胞的原子运动状态分解到单胞中的多个 q 点上.
// 构建基
// 外层下标对应超胞倒格子的整数倍那部分(第二部分), 也就是对应不同反折叠后的 q 点
// 外层下标对应超胞倒格子的整数倍那部分(第二部分), 也就是对应不同反折叠后的 q 点(sub qpoint)
// 内层下标对应单胞倒格子的整数倍那部分(第一部分), 也就是对应同一个反折叠后的 q 点上的不同平面波
std::vector<std::vector<Eigen::VectorXcd>> basis;
basis.resize(input.SuperCellMultiplier.determinant());
for (int i = 0; i < input.SuperCellMultiplier(0, 0); i++)
for (int j = 0; j < input.SuperCellMultiplier(1, 1); j++)
for (int k = 0; k < input.SuperCellMultiplier(2, 2); k++)
{
// 反折叠后的某个 q 点 对应的所有的基矢
auto& basis_at_unfolded_qpoint
= basis[i * input.SuperCellMultiplier(1, 1) * input.SuperCellMultiplier(2, 2)
+ j * input.SuperCellMultiplier(2, 2) + k];
basis_at_unfolded_qpoint.resize(input.PrimativeCellBasisNumber.prod());
for (int x = 0; x < input.PrimativeCellBasisNumber(0); x++)
for (int y = 0; y < input.PrimativeCellBasisNumber(1); y++)
for (int z = 0; z < input.PrimativeCellBasisNumber(2); z++)
{
// 计算 q 点的坐标, 单位为相对于超胞的倒格矢
auto qpoint_relative_to_super_cell =
Eigen::Vector3i({{i}, {j}, {k}})
+ input.SuperCellMultiplier * Eigen::Vector3i({{x}, {y}, {z}});
// 将 q 点坐标转换为埃^-1
auto qpoint = (qpoint_relative_to_super_cell.transpose().cast<double>()
* (input.SuperCellMultiplier.cast<double>().inverse().transpose())).transpose();
// 计算基矢
auto& single_basis = basis_at_unfolded_qpoint
[x * input.PrimativeCellBasisNumber(1) * input.PrimativeCellBasisNumber(2)
+ y * input.PrimativeCellBasisNumber(2) + z];
single_basis = (-2 * std::numbers::pi_v<double> * 1i * (input.AtomPosition * qpoint)).array().exp();
}
}
basis.resize(input.SuperCellMultiplier.diagonal().prod());
for (auto [xyz_of_sub_qpoint, i_of_sub_qpoint]
: triplet_sequence(input.SuperCellMultiplier.diagonal()))
{
basis[i_of_sub_qpoint].resize(input.PrimativeCellBasisNumber.prod());
for (auto [xyz_of_basis, i_of_basis] : triplet_sequence(input.PrimativeCellBasisNumber))
{
// 计算 q 点的坐标, 单位为单胞的倒格矢
auto qpoint_relative_to_primative_cell = xyz_of_basis.cast<double>()
+ input.SuperCellMultiplier.cast<double>().inverse() * xyz_of_sub_qpoint.cast<double>();
// 将 q 点坐标转换为埃^-1
auto qpoint = (qpoint_relative_to_primative_cell.transpose() * (input.PrimativeCell.transpose().inverse()))
.transpose();
// 计算基矢
basis[i_of_sub_qpoint][i_of_basis]
= (-2 * std::numbers::pi_v<double> * 1i * (input.AtomPosition * qpoint)).array().exp();
}
}
// 计算投影的结果
// 最外层下标对应反折叠前的 q 点, 第二层下标对应不同模式, 第三层下标对应这个模式在反折叠后的 q 点
std::vector<std::vector<std::vector<double>>> projection_coefficient;
projection_coefficient.resize(input.QPointData.size());
for (int i_of_folded_qpoint = 0; i_of_folded_qpoint < input.QPointData.size(); i_of_folded_qpoint++)
for (unsigned i_of_folded_qpoint = 0; i_of_folded_qpoint < input.QPointData.size(); i_of_folded_qpoint++)
{
projection_coefficient[i_of_folded_qpoint].resize(input.QPointData[i_of_folded_qpoint].ModeData.size());
for (int i_of_mode = 0; i_of_mode < projection_coefficient[i_of_folded_qpoint].size(); i_of_mode++)
auto num_of_mode = input.QPointData[i_of_folded_qpoint].ModeData.size();
projection_coefficient[i_of_folded_qpoint].resize(num_of_mode);
for (unsigned i_of_mode = 0; i_of_mode < num_of_mode; i_of_mode++)
{
auto& coefficient_at_mode = projection_coefficient[i_of_folded_qpoint][i_of_mode];
coefficient_at_mode.resize(input.SuperCellMultiplier.determinant());
for
(int i_of_unfolded_qpoint = 0; i_of_unfolded_qpoint < coefficient_at_mode.size(); i_of_unfolded_qpoint++)
auto num_of_sub_qpoint = input.SuperCellMultiplier.diagonal().prod();
projection_coefficient[i_of_folded_qpoint][i_of_mode].resize(num_of_sub_qpoint);
for (unsigned i_of_sub_qpoint = 0; i_of_sub_qpoint < num_of_sub_qpoint; i_of_sub_qpoint++)
// 对于 basis 中, 对应于单胞倒格子的部分, 以及对应于不同方向的部分, 分别求内积, 然后求绝对值, 然后求和
for
(
int i_of_basis_in_primary_cell = 0;
i_of_basis_in_primary_cell < basis[i_of_unfolded_qpoint].size();
i_of_basis_in_primary_cell++
)
coefficient_at_mode[i_of_unfolded_qpoint] +=
// 最后, 还应该除以原子数
for (unsigned i_of_basis = 0; i_of_basis < input.PrimativeCellBasisNumber.prod(); i_of_basis++)
projection_coefficient[i_of_folded_qpoint][i_of_mode][i_of_sub_qpoint] +=
(
basis[i_of_unfolded_qpoint][i_of_basis_in_primary_cell].transpose()
* input.QPointData[i_of_folded_qpoint].ModeData[i_of_mode].AtomMovement
).array().abs2().sum();
// 归一化
auto sum = std::accumulate(coefficient_at_mode.begin(), coefficient_at_mode.end(), 0.);
for (auto& coefficient : coefficient_at_mode)
coefficient /= sum;
basis[i_of_sub_qpoint][i_of_basis].transpose()
* input.QPointData[i_of_folded_qpoint].ModeData[i_of_mode].AtomMovement
).array().abs().sum() / input.AtomPosition.rows();
}
}
// 填充输出对象
Output output;
for (int i_of_folded_qpoint = 0; i_of_folded_qpoint < input.QPointData.size(); i_of_folded_qpoint++)
// for each unfolded q point corresponding to this folded q point
for (int x = 0; x < input.SuperCellMultiplier(0, 0); x++)
for (int y = 0; y < input.SuperCellMultiplier(1, 1); y++)
for (int z = 0; z < input.SuperCellMultiplier(2, 2); z++)
{
auto& unfolded_qpoint = output.QPointData.emplace_back();
unfolded_qpoint.QPoint = input.SuperCellMultiplier.cast<double>().inverse() *
(input.QPointData[i_of_folded_qpoint].QPoint
+ Eigen::Vector3i({{x}, {y}, {z}}).cast<double>());
for (int i_of_mode = 0; i_of_mode < input.QPointData[i_of_folded_qpoint].ModeData.size(); i_of_mode++)
if (projection_coefficient[i_of_folded_qpoint][i_of_mode]
[x * input.SuperCellMultiplier(1, 1) * input.SuperCellMultiplier(2, 2)
+ y * input.SuperCellMultiplier(2, 2) + z] > 1e-3)
{
auto& mode = unfolded_qpoint.ModeData.emplace_back();
mode.Frequency = input.QPointData[i_of_folded_qpoint].ModeData[i_of_mode].Frequency;
mode.Weight = projection_coefficient[i_of_folded_qpoint][i_of_mode]
[x * input.SuperCellMultiplier(1, 1) * input.SuperCellMultiplier(2, 2)
+ y * input.SuperCellMultiplier(2, 2) + z];
}
}
for (unsigned i_of_folded_qpoint = 0; i_of_folded_qpoint < input.QPointData.size(); i_of_folded_qpoint++)
for (auto [xyz_of_sub_qpoint, i_of_sub_qpoint]
: triplet_sequence(input.SuperCellMultiplier.diagonal()))
{
auto& sub_qpoint = output.QPointData.emplace_back();
sub_qpoint.QPoint = input.SuperCellMultiplier.cast<double>().inverse() *
(input.QPointData[i_of_folded_qpoint].QPoint + xyz_of_sub_qpoint.cast<double>());
for (unsigned i_of_mode = 0; i_of_mode < input.QPointData[i_of_folded_qpoint].ModeData.size(); i_of_mode++)
{
auto& mode = sub_qpoint.ModeData.emplace_back();
mode.Frequency = input.QPointData[i_of_folded_qpoint].ModeData[i_of_mode].Frequency;
mode.Weight = projection_coefficient[i_of_folded_qpoint][i_of_mode][i_of_sub_qpoint];
mode.Source = input.QPointData[i_of_folded_qpoint].QPoint;
}
}
output.write(argv[2]);
std::ofstream(argv[2]) << YAML::Node(output);
}
Input::Input(YAML::Node root)
bool YAML::convert<Input>::decode(const Node& node, Input& input)
{
for (int i = 0; i < 3; i++)
for (int j = 0; j < 3; j++)
PrimativeCell(i, j) = root["lattice"][i][j].as<double>();
for (unsigned i = 0; i < 3; i++)
for (unsigned j = 0; j < 3; j++)
input.PrimativeCell(i, j) = node["lattice"][i][j].as<double>();
SuperCellMultiplier.setZero();
for (int i = 0; i < 3; i++)
SuperCellMultiplier(i, i) = root["SuperCellMultiplier"][i].as<int>();
input.SuperCellMultiplier.setZero();
for (unsigned i = 0; i < 3; i++)
input.SuperCellMultiplier(i, i) = node["SuperCellMultiplier"][i].as<int>();
for (int i = 0; i < 3; i++)
PrimativeCellBasisNumber(i) = root["PrimativeCellBasisNumber"][i].as<int>();
for (unsigned i = 0; i < 3; i++)
input.PrimativeCellBasisNumber(i) = node["PrimativeCellBasisNumber"][i].as<int>();
auto points = root["points"].as<std::vector<YAML::Node>>();
auto points = node["points"].as<std::vector<YAML::Node>>();
auto atom_position_to_super_cell = Eigen::MatrixX3d(points.size(), 3);
for (int i = 0; i < points.size(); i++)
for (int j = 0; j < 3; j++)
for (unsigned i = 0; i < points.size(); i++)
for (unsigned j = 0; j < 3; j++)
atom_position_to_super_cell(i, j) = points[i]["coordinates"][j].as<double>();
AtomPosition = atom_position_to_super_cell * (SuperCellMultiplier.cast<double>() * PrimativeCell);
input.AtomPosition = atom_position_to_super_cell * (input.SuperCellMultiplier.cast<double>() * input.PrimativeCell);
auto phonon = root["phonon"].as<std::vector<YAML::Node>>();
QPointData.resize(phonon.size());
for (int i = 0; i < phonon.size(); i++)
auto phonon = node["phonon"].as<std::vector<YAML::Node>>();
input.QPointData.resize(phonon.size());
for (unsigned i = 0; i < phonon.size(); i++)
{
QPointData[i].QPoint.resize(3);
for (int j = 0; j < 3; j++)
QPointData[i].QPoint(j) = phonon[i]["q-position"][j].as<double>();
input.QPointData[i].QPoint.resize(3);
for (unsigned j = 0; j < 3; j++)
input.QPointData[i].QPoint(j) = phonon[i]["q-position"][j].as<double>();
auto band = phonon[i]["band"].as<std::vector<YAML::Node>>();
QPointData[i].ModeData.resize(band.size());
for (int j = 0; j < band.size(); j++)
input.QPointData[i].ModeData.resize(band.size());
for (unsigned j = 0; j < band.size(); j++)
{
QPointData[i].ModeData[j].Frequency = band[j]["frequency"].as<double>();
auto eigenvectors = Eigen::MatrixX3cd(AtomPosition.rows(), 3);
input.QPointData[i].ModeData[j].Frequency = band[j]["frequency"].as<double>();
auto eigenvectors = Eigen::MatrixX3cd(input.AtomPosition.rows(), 3);
auto eigenvector_vectors = band[j]["eigenvector"]
.as<std::vector<std::vector<std::vector<double>>>>();
for (int k = 0; k < AtomPosition.rows(); k++)
for (int l = 0; l < 3; l++)
for (unsigned k = 0; k < input.AtomPosition.rows(); k++)
for (unsigned l = 0; l < 3; l++)
eigenvectors(k, l)
= eigenvector_vectors[k][l][0] + 1i * eigenvector_vectors[k][l][1];
// 需要对读入的原子运动状态作相位转换, 使得它们与我们的约定一致(对超胞周期性重复)
// QPointData[i].ModeData[j].AtomMovement
// = eigenvectors.array().colwise() * (-2 * std::numbers::pi_v<double> * 1i
// * (AtomPosition * QPointData[i].QPoint)).array().exp();
QPointData[i].ModeData[j].AtomMovement.resize(AtomPosition.rows(), 3);
for (int k = 0; k < AtomPosition.rows(); k++)
QPointData[i].ModeData[j].AtomMovement.row(k) = eigenvectors.row(k)
* std::exp(-2 * std::numbers::pi_v<double> * 1i
* (atom_position_to_super_cell.row(k).dot(QPointData[i].QPoint)));
// print AtomMovement
// std::cout << "AtomMovement" << std::endl;
// std::cout << QPointData[i].ModeData[j].AtomMovement << std::endl;
input.QPointData[i].ModeData[j].AtomMovement
= eigenvectors.array().colwise() * (-2 * std::numbers::pi_v<double> * 1i
* (atom_position_to_super_cell * input.QPointData[i].QPoint)).array().exp();
// input.QPointData[i].ModeData[j].AtomMovement.resize(input.AtomPosition.rows(), 3);
// for (unsigned k = 0; k < input.AtomPosition.rows(); k++)
// input.QPointData[i].ModeData[j].AtomMovement.row(k) = eigenvectors.row(k)
// * std::exp(-2 * std::numbers::pi_v<double> * 1i
// * (atom_position_to_super_cell.row(k).dot(input.QPointData[i].QPoint)));
// 这里还要需要做归一化处理
// phonopy 的文档似乎没有保证一定做了归一化处理, 再来做一遍吧.
auto sum = QPointData[i].ModeData[j].AtomMovement.cwiseAbs2().sum();
QPointData[i].ModeData[j].AtomMovement /= std::sqrt(sum);
// 这里的归一化其实并不是必须的, 因为在计算投影系数的时候, 还会进行一次归一化.
// 因为计算量并不大, 因此这里还是做一次归一化, 达到形式上的正确.
auto average =
({
auto& _ = input.QPointData[i].ModeData[j].AtomMovement;
_.cwiseAbs2().rowwise().sum().cwiseSqrt().sum() / _.rows();
});
input.QPointData[i].ModeData[j].AtomMovement /= average;
}
}
return true;
}
void Output::write(std::string filename) const
auto YAML::convert<Output>::encode(const Output& output) -> Node
{
YAML::Node root;
root["QPointData"] = YAML::Node(YAML::NodeType::Sequence);
for (int i = 0; i < QPointData.size(); i++)
Node node;
node["QPointData"] = Node(NodeType::Sequence);
for (unsigned i = 0; i < output.QPointData.size(); i++)
{
root["QPointData"][i]["QPoint"] = YAML::Node(YAML::NodeType::Sequence);
for (int j = 0; j < 3; j++)
root["QPointData"][i]["QPoint"][j] = QPointData[i].QPoint(j);
root["QPointData"][i]["ModeData"] = YAML::Node(YAML::NodeType::Sequence);
for (int j = 0; j < QPointData[i].ModeData.size(); j++)
node["QPointData"][i]["QPoint"] =
({
auto& _ = output.QPointData[i].QPoint;
std::vector<double>(_.data(), _.data() + _.size());
});
node["QPointData"][i]["ModeData"] = Node(NodeType::Sequence);
for (unsigned j = 0; j < output.QPointData[i].ModeData.size(); j++)
{
root["QPointData"][i]["ModeData"][j]["Frequency"] = QPointData[i].ModeData[j].Frequency;
root["QPointData"][i]["ModeData"][j]["Weight"] = QPointData[i].ModeData[j].Weight;
node["QPointData"][i]["ModeData"][j]["Frequency"] = output.QPointData[i].ModeData[j].Frequency;
node["QPointData"][i]["ModeData"][j]["Weight"] = output.QPointData[i].ModeData[j].Weight;
node["QPointData"][i]["ModeData"][j]["Source"] =
({
auto& _ = output.QPointData[i].ModeData[j].Source;
std::vector<double>(_.data(), _.data() + _.size());
});
}
}
std::ofstream(filename) << root;
return node;
}