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

init

Browse Source
This commit is contained in:
陈浩南
2024-11-21 16:00:37 +08:00
parent d1432c7591
commit 28e4d29f2c
29 changed files with 985 additions and 2677 deletions
Download Patch File Download Diff File Expand all files Collapse all files

2
.envrc
View File

@@ -1 +1 @@
use flake
use flake nixos#ufo

18
.gitignore vendored
View File

@@ -1,18 +1,4 @@
.vscode
.direnv
main
plot
out.yaml
.ccls-cache
build
.vscode
.cache
build/
test/*
!test/14.2.2.4.unfold.yaml
!test/14.2.2.5.unfold.yaml
!test/14.2.2.5.plot.yaml
!test/14.2.4.4.unfold.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

33
CMakeLists.txt
View File

@@ -1,36 +1,31 @@
cmake_minimum_required(VERSION 3.14)
project(ufo VERSION 0.0.0 LANGUAGES CXX)
project(ufo VERSION 0 LANGUAGES CXX)
enable_testing()
include(GNUInstallDirs)
set_property(GLOBAL PROPERTY CXX_STANDARD 23)
set_property(GLOBAL PROPERTY CXX_STANDARD_REQUIRED ON)
set_property(GLOBAL PROPERTY CXX_EXTENSIONS OFF)
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(yaml-cpp REQUIRED)
find_package(Eigen3 REQUIRED)
find_package(fmt REQUIRED)
find_package(concurrencpp REQUIRED)
find_package(HighFive REQUIRED)
find_package(TBB REQUIRED)
find_package(glad REQUIRED)
# find_package(Matplot++ REQUIRED COMPONENTS matplot_opengl)
find_package(Matplot++ REQUIRED)
find_path(ZPP_BITS_INCLUDE_DIR zpp_bits.h REQUIRED)
find_package(biu REQUIRED)
find_package(Threads REQUIRED)
add_executable(ufo src/solver.cpp src/fold.cpp src/unfold.cpp src/plot.cpp src/main.cpp)
target_include_directories(ufo PRIVATE ${PROJECT_SOURCE_DIR}/include ${ZPP_BITS_INCLUDE_DIR})
target_link_libraries(ufo PRIVATE
yaml-cpp Eigen3::Eigen fmt::fmt concurrencpp::concurrencpp HighFive_HighFive TBB::tbb Matplot++::matplot
Matplot++::matplot_opengl)
# target_compile_definitions(ufo PRIVATE $<$<CONFIG:Debug>:_GLIBCXX_DEBUG>)
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)

17
README.md
View File

@@ -1,17 +0,0 @@
This is a program to approximately unfold phonon spectra.
Currently, there is no documentation, and some not core features are in other branches
and implemented in a temporary way, e.g., directly writing data in source file.
I will finish the documentation and merge these features into main branch ASAP.
The program is packaged in a standard way of cmake, so that you can build it just like other cmake projects:
```bash
mkdir build
cd build
cmake ..
make
```
and view source files for input and output format.

11
default.nix Normal file
View File

@@ -0,0 +1,11 @@
{
stdenv, cmake, pkg-config, version ? null,
tbb, matplotplusplus, biu
}: stdenv.mkDerivation
{
name = "ufo";
src = ./.;
buildInputs = [ tbb matplotplusplus biu ];
nativeBuildInputs = [ cmake pkg-config ];
doCheck = true;
}

45
doc/README.md Normal file
View File

@@ -0,0 +1,45 @@
分为几个功能：
* 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]
```

1370
flake.lock generated
View File

File diff suppressed because it is too large Load Diff

25
flake.nix
View File

@@ -1,25 +0,0 @@
{
inputs.nixos.url = "github:CHN-beta/nixos";
inputs.nixpkgs.url = "github:NixOS/nixpkgs/nixos-unstable";
outputs = inputs:
let
# pkgs = inputs.nixpkgs.legacyPackages.x86_64-linux;
pkgs = import inputs.nixpkgs
{
localSystem = { system = "x86_64-linux"; gcc = { arch = "alderlake"; tune = "alderlake"; }; };
config.allowUnfree = true;
};
localPackages = import "${inputs.nixos}/local/pkgs" { inherit pkgs; inherit (inputs.nixpkgs) lib; };
in
{
devShell.x86_64-linux = pkgs.mkShell.override { stdenv = pkgs.gcc13Stdenv; }
{
packages = with pkgs; [ pkg-config cmake ninja ];
buildInputs = (with pkgs; [ eigen yaml-cpp fmt highfive tbb_2021_8.dev glfw libGL ])
++ (with localPackages; [ concurrencpp matplotplusplus zpp-bits glad ]);
hardeningDisable = [ "all" ];
# NIX_DEBUG = "1";
};
};
}

74
include/ufo.hpp Normal file
View File

@@ -0,0 +1,74 @@
# 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>;
};
}

37
include/ufo/fold.hpp
View File

@@ -1,37 +0,0 @@
# pragma once
# include <ufo/solver.hpp>
namespace ufo
{
class FoldSolver : public Solver
{
public:
struct InputType
{
Eigen::Vector<unsigned, 3> SuperCellMultiplier;
std::optional<Eigen::Matrix<double, 3, 3>> SuperCellDeformation;
std::vector<Eigen::Vector3d> Qpoints;
DataFile OutputFile;
InputType(std::string config_file);
};
struct OutputType
{
std::vector<Eigen::Vector3d> Qpoints;
void write(std::string filename) const;
};
protected:
InputType Input_;
std::optional<OutputType> Output_;
public:
FoldSolver(std::string config_file);
FoldSolver& operator()() override;
// return value: QpointInReciprocalSuperCellByReciprocalSuperCell
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
);
};
}

77
include/ufo/plot.hpp
View File

@@ -1,77 +0,0 @@
# pragma once
# include <ufo/unfold.hpp>
namespace ufo
{
class PlotSolver : public Solver
{
public:
struct InputType
{
Eigen::Matrix3d PrimativeCell;
struct FigureConfigType
{
std::vector<std::vector<Eigen::Vector3d>> Qpoints;
std::pair<unsigned, unsigned> Resolution;
std::pair<double, double> Range;
std::optional<std::vector<double>> YTicks;
DataFile PictureFile;
std::optional<std::vector<DataFile>> DataFiles;
};
std::vector<FigureConfigType> Figures;
struct UnfoldedDataType : public UnfoldSolver::OutputType
{
UnfoldedDataType(std::string filename);
UnfoldedDataType() = default;
};
DataFile UnfoldedDataFile;
UnfoldedDataType UnfoldedData;
InputType(std::string config_file);
};
struct OutputType
{
std::vector<std::vector<double>> Values;
std::vector<double> XTicks;
std::vector<double> YTicks;
std::pair<unsigned, unsigned> Resolution;
std::pair<double, double> Range;
OutputType() = default;
const OutputType& write(std::string filename, std::string format) const;
using serialize = zpp::bits::members<5>;
};
protected:
InputType Input_;
std::optional<std::vector<OutputType>> Output_;
public:
PlotSolver(std::string config_file);
PlotSolver& operator()() override;
// 根据 q 点路径, 搜索要使用的 q 点
static std::vector<std::reference_wrapper<const UnfoldSolver::OutputType::QpointDataType>> search_qpoints
(
const std::pair<Eigen::Vector3d, Eigen::Vector3d>& path,
const decltype(InputType::UnfoldedDataType::QpointData)& available_qpoints,
double threshold, bool exclude_endpoint = false
);
// 根据搜索到的 q 点, 计算每个点的数值
static std::tuple<std::vector<std::vector<double>>, std::vector<double>> calculate_values
(
const Eigen::Matrix3d primative_cell,
const std::vector<std::pair<Eigen::Vector3d, Eigen::Vector3d>>& path,
const std::vector<std::vector<std::reference_wrapper<const UnfoldSolver::OutputType::QpointDataType>>>& qpoints,
const decltype(InputType::FigureConfigType::Resolution)& resolution,
const decltype(InputType::FigureConfigType::Range)& range
);
// 根据数值, 画图
static void plot
(
const std::vector<std::vector<double>>& values,
const std::string& filename,
const std::vector<double>& x_ticks, const std::vector<double>& y_ticks
);
};
}

141
include/ufo/solver.hpp
View File

@@ -1,141 +0,0 @@
# pragma once
# include <iostream>
# include <array>
# include <numbers>
# include <numeric>
# include <fstream>
# include <optional>
# include <array>
# include <utility>
# include <execution>
# include <syncstream>
# include <any>
# include <map>
# include <vector>
# include <span>
# include <yaml-cpp/yaml.h>
# include <Eigen/Dense>
# include <concurrencpp/concurrencpp.h>
# include <fmt/format.h>
# include <fmt/std.h>
# include <fmt/ranges.h>
# include <highfive/H5File.hpp>
# include <zpp_bits.h>
# include <matplot/matplot.h>
# include <matplot/backend/opengl.h>
// 在相位中, 约定为使用 $\exp (2 \pi i \vec{q} \cdot \vec{r})$ 来表示原子的运动状态
// (而不是 $\exp (-2 \pi i \vec{q} \cdot \vec{r})$)
// 一些书定义的倒格矢中包含了 $2 \pi$ 的部分, 我们这里约定不包含这部分.
// 也就是说, 正格子与倒格子的转置相乘, 得到单位矩阵.
namespace Eigen
{
constexpr inline auto serialize(auto & archive, Eigen::Matrix3d& matrix)
{ return archive(std::span(matrix.data(), matrix.size())); }
constexpr inline auto serialize(auto & archive, const Eigen::Matrix3d& matrix)
{ return archive(std::span(matrix.data(), matrix.size())); }
constexpr inline auto serialize(auto & archive, Eigen::Vector3d& vector)
{ return archive(std::span(vector.data(), vector.size())); }
constexpr inline auto serialize(auto & archive, const Eigen::Vector3d& vector)
{ return archive(std::span(vector.data(), vector.size())); }
}
namespace ufo
{
using namespace std::literals;
struct PhonopyComplex { double r, i; };
inline HighFive::CompoundType create_compound_complex()
{ return {{ "r", HighFive::AtomicType<double>{}}, {"i", HighFive::AtomicType<double>{}}}; }
namespace detail_
{
template <typename T> struct SpecializationOfBitsMembersHelper : std::false_type {};
template <std::size_t N> struct SpecializationOfBitsMembersHelper<zpp::bits::members<N>> : std::true_type {};
}
template <typename T> concept ZppSerializable
= requires() { detail_::SpecializationOfBitsMembersHelper<T>::value == true; };
class Solver
{
public:
virtual Solver& operator()() = 0;
virtual ~Solver() = default;
static concurrencpp::generator<std::pair<Eigen::Vector<unsigned, 3>, unsigned>>
triplet_sequence(Eigen::Vector<unsigned, 3> range);
template <ZppSerializable T> inline static void zpp_write(const T& object, std::string filename)
{
auto [data, out] = zpp::bits::data_out();
out(object).or_throw();
static_assert(sizeof(char) == sizeof(std::byte));
std::ofstream file(filename, std::ios::binary | std::ios::out);
file.exceptions(std::ios::badbit | std::ios::failbit);
file.write(reinterpret_cast<const char*>(data.data()), data.size());
}
template <ZppSerializable T> inline static T zpp_read(std::string filename)
{
auto input = std::ifstream(filename, std::ios::binary | std::ios::in);
input.exceptions(std::ios::badbit | std::ios::failbit);
static_assert(sizeof(std::byte) == sizeof(char));
std::vector<std::byte> data;
{
std::vector<char> string(std::istreambuf_iterator<char>(input), {});
data.assign
(
reinterpret_cast<std::byte*>(string.data()),
reinterpret_cast<std::byte*>(string.data() + string.size())
);
}
auto in = zpp::bits::in(data);
T output;
in(output).or_throw();
return output;
}
class Hdf5file
{
public:
inline Hdf5file& open_for_read(std::string filename)
{
File_ = HighFive::File(filename, HighFive::File::ReadOnly);
return *this;
}
inline Hdf5file& open_for_write(std::string filename)
{
File_ = HighFive::File(filename, HighFive::File::ReadWrite | HighFive::File::Create
| HighFive::File::Truncate);
return *this;
}
template <typename T> inline Hdf5file& read(T& object, std::string name)
{
object = File_->getDataSet(name).read<std::remove_cvref_t<decltype(object)>>();
return *this;
}
template <typename T> inline Hdf5file& write(const T& object, std::string name)
{
File_->createDataSet(name, object);
return *this;
}
protected:
std::optional<HighFive::File> File_;
};
struct DataFile
{
std::string Filename;
std::string Format;
std::map<std::string, std::any> ExtraParameters;
inline DataFile() = default;
DataFile
(
YAML::Node node, std::set<std::string> supported_format,
std::string config_file, bool allow_same_as_config_file = false
);
};
};
}
HIGHFIVE_REGISTER_TYPE(ufo::PhonopyComplex, ufo::create_compound_complex)

164
include/ufo/unfold.hpp
View File

@@ -1,164 +0,0 @@
# pragma once
# include <ufo/solver.hpp>
namespace ufo
{
// 反折叠的原理: 将超胞中的原子运动状态, 投影到一组平面波构成的基矢中.
// 每一个平面波的波矢由两部分相加得到: 一部分是单胞倒格子的整数倍, 所取的个数有一定任意性, 论文中建议取大约单胞中原子个数那么多个;
// 对于没有缺陷的情况, 取一个应该就足够了.
// 另一部分是超胞倒格子的整数倍, 取 n 个, n 为超胞对应的单胞的倍数, 其实也就是倒空间中单胞对应倒格子中超胞的格点.
// 只要第一部分取得足够多, 那么单胞中原子的状态就可以完全被这些平面波描述.
// 将超胞中原子的运动状态投影到这些基矢上, 计算出投影的系数, 就可以将超胞的原子运动状态分解到单胞中的多个 q 点上.
class UnfoldSolver : public Solver
{
public:
struct InputType
{
// 单胞的三个格矢，每行表示一个格矢的坐标，单位为埃
Eigen::Matrix3d PrimativeCell;
// 单胞到超胞的格矢转换时用到的矩阵
// SuperCellMultiplier 是一个三维列向量且各个元素都是整数，表示单胞在各个方向扩大到多少倍之后，可以得到和超胞一样的体积
// SuperCsolver.hpp>mation 是一个行列式为 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::Vector<unsigned, 3> SuperCellMultiplier;
std::optional<Eigen::Matrix<double, 3, 3>> SuperCellDeformation;
// 在单胞内取几个平面波的基矢
Eigen::Vector<unsigned, 3> PrimativeCellBasisNumber;
// 从哪个文件读入 AtomPosition, 以及这个文件的格式, 格式可选值包括 "yaml"
DataFile AtomPositionInputFile;
// 从哪个文件读入 QpointData, 以及这个文件的格式, 格式可选值包括 "yaml" 和 "hdf5"
DataFile QpointDataInputFile;
// 超胞中原子的坐标，每行表示一个原子的坐标，单位为埃
Eigen::MatrixX3d AtomPosition;
// 关于各个 Q 点的数据
struct QpointDataType
{
// Q 点的坐标，单位为超胞的倒格矢
Eigen::Vector3d Qpoint;
// 关于这个 Q 点上各个模式的数据
struct ModeDataType
{
// 模式的频率，单位为 THz
double Frequency;
// 模式中各个原子的运动状态
// 这个数据是这样得到的: phonopy 输出的动态矩阵的 eigenvector 乘以 $\exp(-2 \pi i \vec q \cdot \vec r)$
// 这个数据可以认为是原子位移中, 关于超胞有周期性的那一部分, 再乘以原子质量的开方.
// 这个数据在读入后会被立即归一化.
Eigen::MatrixX3cd AtomMovement;
};
std::vector<ModeDataType> ModeData;
};
std::vector<QpointDataType> QpointData;
// 输出到哪些文件, 以及使用怎样的格式, 格式可选值包括:
// yaml: 使用 yaml 格式输出
// yaml-human-readable: 使用 yaml 格式输出, 但是输出的结果更适合人类阅读,
// 包括合并相近的模式, 去除权重过小的模式, 限制输出的小数位数.
// zpp: 使用 zpp-bits 序列化, 可以直接被 plot.cpp 读取
std::vector<DataFile> QpointDataOutputFile;
// 从文件中读取输入 (包括一个较小的配置文件, 和一个 hdf5 或者一个 yaml 文件), 文件中应当包含:
// 单胞的格矢: PrimativeCell 单位为埃 直接从 phonopy 的输出中复制
// 超胞的倍数: SuperCellMultiplier 手动输入, 为一个包含三个整数的数组
// 超胞的变形: SuperCellDeformation 手动输入, 为一个三阶方阵
// 平面波的基矢个数: PrimativeCellBasisNumber 手动输入, 为一个包含三个整数的数组
// 另外还有一个文件, 直接将 phonopy 的输出复制过来即可, 如果是 yaml, 应该包含下面的内容:
// 超胞中原子的坐标: points[*].coordinates 单位为超胞的格矢 直接从 phonopy 的输出中复制
// 各个 Q 点的坐标: phonon[*].q-position 单位为超胞的倒格子的格矢 直接从 phonopy 的输出中复制
// 各个模式的频率: phonon[*].band[*].frequency 单位为 THz 直接从 phonopy 的输出中复制
// 各个模式的原子运动状态: phonon[*].band[*].eigenvector 直接从 phonopy 的输出中复制
// 文件中可以有多余的项目, 多余的项目不管.
InputType(std::string filename);
};
struct OutputType
{
// 关于各个 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;
void write(decltype(InputType::QpointDataOutputFile) output_files) const;
void write(std::string filename, std::string format, unsigned percision = 10) const;
using serialize = zpp::bits::members<1>;
virtual ~OutputType() = default;
};
// 第一层是不同的 sub qpoint, 第二层是单胞内不同的平面波
using BasisType = std::vector<std::vector<Eigen::VectorXcd>>;
protected:
InputType Input_;
std::optional<OutputType> Output_;
std::optional<BasisType> Basis_;
// 第一层是不同的模式, 第二层是不同的 sub qpoint
using ProjectionCoefficientType_ = std::vector<std::vector<double>>;
public:
UnfoldSolver(std::string config_file);
UnfoldSolver& operator()() override;
// 构建基
// 每个 q 点对应的一组 sub qpoint。不同的 q 点所对应的 sub qpoint 是不一样的，但 sub qpoint 与 q 点的相对位置一致。
// 这里 xyz_of_diff_of_sub_qpoint 即表示这个相对位置。
// 由于基只与这个相对位置有关（也就是说，不同 q 点的基是一样的），因此可以先计算出所有的基，这样降低计算量。
// 外层下标对应超胞倒格子的整数倍那部分(第二部分), 也就是不同的 sub qpoint
// 内层下标对应单胞倒格子的整数倍那部分(第一部分), 也就是 sub qpoint 上的不同平面波（取的数量越多，结果越精确）
static BasisType construct_basis
(
const decltype(InputType::PrimativeCell)& primative_cell,
const decltype(InputType::SuperCellMultiplier)& super_cell_multiplier,
const decltype(InputType::PrimativeCellBasisNumber)&
primative_cell_basis_number,
const decltype(InputType::AtomPosition)& atom_position
);
// 计算投影系数, 是反折叠的核心步骤
ProjectionCoefficientType_ construct_projection_coefficient
(
const BasisType& basis,
const std::vector<std::reference_wrapper<const decltype
(InputType::QpointDataType::ModeDataType::AtomMovement)>>& mode_data,
std::atomic<unsigned>& number_of_finished_modes
);
OutputType construct_output
(
const decltype(InputType::SuperCellMultiplier)& super_cell_multiplier,
const decltype(InputType::SuperCellDeformation)& super_cell_deformation,
const std::vector<std::reference_wrapper<const decltype
(InputType::QpointDataType::Qpoint)>>& meta_qpoint_by_reciprocal_super_cell,
const std::vector<std::vector<std::reference_wrapper<const decltype
(InputType::QpointDataType::ModeDataType::Frequency)>>>& frequency,
const ProjectionCoefficientType_& projection_coefficient
);
};
}

122
src/fold.cpp
View File

@@ -1,86 +1,58 @@
# include <ufo/fold.hpp>
# include <ufo.hpp>
namespace ufo
void ufo::fold(std::string config_file)
{
FoldSolver::InputType::InputType(std::string config_file)
struct Input
{
auto input = YAML::LoadFile(config_file);
for (unsigned i = 0; i < 3; i++)
SuperCellMultiplier(i) = input["SuperCellMultiplier"][i].as<unsigned>();
if (input["SuperCellDeformation"])
{
SuperCellDeformation.emplace();
for (unsigned i = 0; i < 3; i++)
for (unsigned j = 0; j < 3; j++)
(*SuperCellDeformation)(i, j) = input["SuperCellDeformation"][i][j].as<double>();
}
for (auto& qpoint : input["Qpoints"].as<std::vector<std::vector<double>>>())
Qpoints.push_back(Eigen::Vector3d
{{qpoint.at(0)}, {qpoint.at(1)}, {qpoint.at(2)}});
OutputFile = DataFile(input["OutputFile"], {"yaml"}, config_file);
}
void FoldSolver::OutputType::write(std::string filename) const
Eigen::Matrix3d SuperCellDeformation;
Eigen::Vector3i SuperCellMultiplier;
std::vector<Eigen::Vector3d> Qpoints;
std::optional<std::string> OutputFile;
};
struct Output
{
std::ofstream(filename) << [&]
{
std::stringstream print;
print << "Qpoints:\n";
for (auto& qpoint : Qpoints)
print << fmt::format(" - [ {:.8f}, {:.8f}, {:.8f} ]\n", qpoint(0), qpoint(1), qpoint(2));
return print.str();
}();
}
FoldSolver::FoldSolver(std::string config_file) : Input_(config_file) {}
FoldSolver& FoldSolver::operator()()
{
if (!Output_)
{
Output_.emplace();
for (auto& qpoint : Input_.Qpoints)
Output_->Qpoints.push_back(fold
(
qpoint, Input_.SuperCellMultiplier,
Input_.SuperCellDeformation
));
}
Output_->write(Input_.OutputFile.Filename);
return *this;
}
Eigen::Vector3d FoldSolver::fold
std::vector<Eigen::Vector3d> Qpoints;
};
auto 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
)
Eigen::Matrix3d super_cell_transformation
) -> Eigen::Vector3d
{
// 首先需要将 q 点转移到 ModifiedSuperCell 的倒格子中
// 将 q 点坐标扩大, 然后取小数部分, 就可以了
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 =
(qpoint_by_reciprocal_modified_super_cell.array() - qpoint_by_reciprocal_modified_super_cell.array().floor())
.matrix();
if (!super_cell_deformation)
return qpoint_in_reciprocal_modified_super_cell_by_reciprocal_modified_super_cell;
/*
对 q 点平移数个 SupreCell, 直到它落在超胞的倒格子中
这等价于直接将 q 点坐标用 SuperCell 的倒格子表示, 然后取小数部分.
ModifiedSuperCell = SuperCellMultiplier * PrimativeCell
SuperCell = SuperCellDeformation * ModifiedSuperCell
ReciprocalModifiedSuperCell = ModifiedSuperCell.inverse().transpose()
ReciprocalSuperCell = SuperCell.inverse().transpose()
Qpoint = QpointByReciprocalModifiedSuperCell.transpose() * ReciprocalModifiedSuperCell
Qpoint = QpointByReciprocalSuperCell.transpose() * ReciprocalSuperCell
首先需要将 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 * QpointByReciprocalModifiedSuperCell;
两个式子结合，可以得到：
QpointByReciprocalSuperCell = SuperCellDeformation * SuperCellMultiplier * QpointByReciprocalPrimitiveCell;
*/
auto qpoint_in_reciprocal_modified_super_cell_by_reciprocal_super_cell =
(*super_cell_deformation * qpoint_in_reciprocal_modified_super_cell_by_reciprocal_modified_super_cell).eval();
auto qpoint_in_reciprocal_super_cell_by_reciprocal_super_cell =
qpoint_in_reciprocal_modified_super_cell_by_reciprocal_super_cell.array()
- qpoint_in_reciprocal_modified_super_cell_by_reciprocal_super_cell.array().floor();
return qpoint_in_reciprocal_super_cell_by_reciprocal_super_cell;
}
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);
}

23
src/main.cpp
View File

@@ -1,17 +1,14 @@
# include <ufo/fold.hpp>
# include <ufo/unfold.hpp>
# include <ufo/plot.hpp>
# include <ufo.hpp>
int main(int argc, const char** argv)
{
if (argc != 3)
throw std::runtime_error(fmt::format("Usage: {} task config.yaml", argv[0]));
if (argv[1] == std::string("fold"))
ufo::FoldSolver{argv[2]}();
else if (argv[1] == std::string("unfold"))
ufo::UnfoldSolver{argv[2]}();
else if (argv[1] == std::string("plot"))
ufo::PlotSolver{argv[2]}();
else
throw std::runtime_error(fmt::format("Unknown task: {}", argv[1]));
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]));
}

575
src/plot.cpp
View File

@@ -1,233 +1,177 @@
# include <ufo/plot.hpp>
# include <ufo.hpp>
# include <matplot/matplot.h>
# include <boost/container/flat_map.hpp>
namespace ufo
void ufo::plot_band(std::string config_file)
{
PlotSolver::InputType::UnfoldedDataType::UnfoldedDataType(std::string filename)
struct Input
{
static_cast<UnfoldSolver::OutputType&>(*this) = zpp_read<UnfoldSolver::OutputType>(filename);
}
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;
};
PlotSolver::InputType::InputType(std::string 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 (auto& figure : input["Figures"].as<std::vector<YAML::Node>>())
{
Figures.emplace_back();
auto qpoints = figure["Qpoints"]
.as<std::vector<std::vector<std::vector<double>>>>();
for (auto& line : qpoints)
{
Figures.back().Qpoints.emplace_back();
for (auto& point : line)
Figures.back().Qpoints.back().emplace_back(point.at(0), point.at(1), point.at(2));
if (Figures.back().Qpoints.back().size() < 2)
throw std::runtime_error("Not enough points in a line");
}
if (Figures.back().Qpoints.size() < 1)
throw std::runtime_error("Not enough lines in a figure");
Figures.back().Resolution = figure["Resolution"].as<std::pair<unsigned, unsigned>>();
Figures.back().Range = figure["Range"].as<std::pair<double, double>>();
Figures.back().PictureFile
= DataFile(figure["PictureFile"], {"png"}, config_file);
if (figure["YTicks"])
Figures.back().YTicks = figure["YTicks"].as<std::vector<double>>();
if (figure["DataFiles"])
{
Figures.back().DataFiles.emplace();
for (auto& data_file : figure["DataFiles"].as<std::vector<YAML::Node>>())
Figures.back().DataFiles->emplace_back()
= DataFile(data_file, {"hdf5", "zpp"}, config_file);
}
}
UnfoldedDataFile = DataFile(input["UnfoldedDataFile"], {"zpp"}, config_file);
UnfoldedData = UnfoldedDataType(UnfoldedDataFile.Filename);
}
const PlotSolver::OutputType& PlotSolver::OutputType::write(std::string filename, std::string format) const
{
if (format == "zpp")
zpp_write(*this, filename);
else if (format == "hdf5")
{
std::vector resolution{ Resolution.first, Resolution.second };
std::vector range{ Range.first, Range.second };
Hdf5file{}.open_for_write(filename).write(Values, "Values")
.write(XTicks, "XTicks")
.write(YTicks, "YTicks")
.write(resolution, "Resolution")
.write(range, "Range");
}
return *this;
}
PlotSolver::PlotSolver(std::string config_file) : Input_(config_file) {}
PlotSolver& PlotSolver::operator()()
{
Output_.emplace();
for (auto& figure : Input_.Figures)
{
// 外层表示不同的线段的端点，内层表示这个线段上的 q 点
std::vector<std::vector<std::reference_wrapper<const UnfoldSolver::OutputType::QpointDataType>>> qpoints;
std::vector<std::pair<Eigen::Vector3d, Eigen::Vector3d>> lines;
for (auto& path : figure.Qpoints)
for (unsigned i = 0; i < path.size() - 1; i++)
{
lines.emplace_back(path[i], path[i + 1]);
qpoints.push_back(search_qpoints
(
lines.back(), Input_.UnfoldedData.QpointData,
0.001, i != path.size() - 2
));
}
auto [values, x_ticks] = calculate_values
(
Input_.PrimativeCell, lines, qpoints, figure.Resolution, figure.Range
);
auto y_ticks = figure.YTicks.value_or(std::vector<double>{});
for (auto& _ : y_ticks)
_ = (_ - figure.Range.first) / (figure.Range.second - figure.Range.first) * figure.Resolution.second;
plot(values, figure.PictureFile.Filename, x_ticks, y_ticks);
Output_->emplace_back();
Output_->back().Values = std::move(values);
Output_->back().XTicks = std::move(x_ticks);
Output_->back().YTicks = std::move(y_ticks);
Output_->back().Resolution = figure.Resolution;
Output_->back().Range = figure.Range;
if (figure.DataFiles)
for (auto& data_file : *figure.DataFiles)
Output_->back().write(data_file.Filename, data_file.Format);
}
return *this;
}
std::vector<std::reference_wrapper<const UnfoldSolver::OutputType::QpointDataType>> PlotSolver::search_qpoints
// 根据 q 点路径, 搜索要使用的 q 点，返回的是 q 点在 QpointData 中的索引以及到路径起点的距离，以及这段路径的总长度
auto search_qpoints = []
(
const Eigen::Matrix3d& primative_cell,
const std::pair<Eigen::Vector3d, Eigen::Vector3d>& path,
const decltype(InputType::UnfoldedDataType::QpointData)& available_qpoints,
double threshold, bool exclude_endpoint
const std::vector<Eigen::Vector3d>& qpoints,
double threshold, bool exclude_endpoint = false
)
{
std::multimap<double, std::reference_wrapper<const UnfoldSolver::OutputType::QpointDataType>> selected_qpoints;
// 对于 output 中的每一个点, 检查这个点是否在路径上. 如果在, 把它加入到 selected_qpoints 中
for (auto& qpoint : available_qpoints)
{
// 计算三点围成的三角形的面积的两倍
auto area = (path.second - path.first).cross(qpoint.Qpoint - path.first).norm();
// 计算这个点到前两个点所在直线的距离
auto distance = area / (path.second - path.first).norm();
// 如果这个点到前两个点所在直线的距离小于阈值, 则认为这个点在路径上
if (distance < threshold)
// 键为这个点到起点的距离
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 distance1 = (qpoint.Qpoint - path.first).norm();
auto distance2 = (qpoint.Qpoint - path.second).norm();
auto distance3 = (path.second - path.first).norm();
if (distance1 < distance3 + threshold && distance2 < distance3 + threshold)
// 如果这个点不在终点处, 或者不排除终点, 则加入
if (distance2 > threshold || !exclude_endpoint)
selected_qpoints.emplace(distance1, std::ref(qpoint));
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 (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");
std::vector<std::reference_wrapper<const UnfoldSolver::OutputType::QpointDataType>> result;
for (auto& qpoint : selected_qpoints)
result.push_back(qpoint.second);
return result;
}
if (selected_qpoints.empty()) throw std::runtime_error("No q points found");
return std::make_pair(selected_qpoints, (end - begin).norm());
};
std::tuple<std::vector<std::vector<double>>, std::vector<double>> PlotSolver::calculate_values
// 根据搜索到的 q 点, 计算图中每个点的值
auto calculate_values = []
(
const Eigen::Matrix3d primative_cell,
const std::vector<std::pair<Eigen::Vector3d, Eigen::Vector3d>>& path,
const std::vector<std::vector<std::reference_wrapper<const UnfoldSolver::OutputType::QpointDataType>>>& qpoints,
const decltype(InputType::FigureConfigType::Resolution)& resolution,
const decltype(InputType::FigureConfigType::Range)& range
// 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
)
{
// 整理输入
std::map<double, std::reference_wrapper<const UnfoldSolver::OutputType::QpointDataType>> qpoints_with_distance;
double total_distance = 0;
std::vector<double> x_ticks;
for (unsigned i = 0; i < path.size(); i++)
{
for (auto& _ : qpoints[i])
qpoints_with_distance.emplace
(
total_distance
+ ((_.get().Qpoint - path[i].first).transpose() * primative_cell.inverse().transpose()).norm(),
_
);
total_distance += ((path[i].second - path[i].first).transpose() * primative_cell.inverse().transpose()).norm();
if (i != path.size() - 1)
x_ticks.push_back(total_distance);
}
for (auto& _ : x_ticks)
_ = _ / total_distance * resolution.first;
// 插值
std::vector<std::vector<double>> values;
auto blend = []
// 按比例混合两个 q 点的结果，得到可以用于画图的那一列数据
auto blend = [&]
(
const UnfoldSolver::OutputType::QpointDataType& a,
const UnfoldSolver::OutputType::QpointDataType& b,
double ratio, unsigned resolution, std::pair<double, double> range
// 两个点的索引
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;
for (unsigned i = 0; i < a.ModeData.size(); i++)
// 如果是连续路径，将每个模式的频率和权重按照比例混合
if (continuous)
{
frequency.push_back(a.ModeData[i].Frequency * ratio + b.ModeData[i].Frequency * (1 - ratio));
weight.push_back(a.ModeData[i].Weight * ratio + b.ModeData[i].Weight * (1 - ratio));
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 (unsigned i = 0; i < frequency.size(); i++)
for (std::size_t i = 0; i < frequency.size(); i++)
{
int index = (frequency[i] - range.first) / (range.second - range.first) * resolution;
if (index >= 0 && index < static_cast<int>(resolution))
result[index] += weight[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;
};
for (unsigned i = 0; i < resolution.first; i++)
std::vector<std::vector<double>> values;
for (std::size_t i = 0; i < resolution[0]; i++)
{
auto current_distance = total_distance * i / resolution.first;
auto it = qpoints_with_distance.lower_bound(current_distance);
if (it == qpoints_with_distance.begin())
values.push_back(blend(it->second.get(), it->second.get(), 1, resolution.second, range));
else if (it == qpoints_with_distance.end())
values.push_back(blend(std::prev(it)->second.get(), std::prev(it)->second.get(), 1, resolution.second,
range));
else
values.push_back(blend
(
std::prev(it)->second.get(), it->second.get(),
(it->first - current_distance) / (it->first - std::prev(it)->first),
resolution.second, range)
);
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, x_ticks};
}
void PlotSolver::plot
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<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>>
@@ -235,32 +179,253 @@ namespace ufo
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 (unsigned i = 0; i < values[0].size(); i++)
for (unsigned j = 0; j < values.size(); j++)
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;
if (v < 0.05) v = 0;
a[i][j] = v * 100 * 255;
if (a[i][j] > 255)
a[i][j] = 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;
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;
if (g[i][j] < 0) g[i][j] = 0;
b[i][j] = 255;
}
auto f = matplot::figure<matplot::backend::gnuplot>(true);
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]);
}

138
src/raman-create-displacement.cpp Normal file
View File

@@ -0,0 +1,138 @@
# 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
);
}

47
src/solver.cpp
View File

@@ -1,47 +0,0 @@
# include <ufo/solver.hpp>
namespace ufo
{
concurrencpp::generator<std::pair<Eigen::Vector<unsigned, 3>, unsigned>> Solver::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
};
}
Solver::DataFile::DataFile
(YAML::Node node, std::set<std::string> supported_format, std::string config_file, bool allow_same_as_config_file)
{
if (auto _ = node["SameAsConfigFile"])
{
auto __ = _.as<bool>();
if (__ && !allow_same_as_config_file)
throw std::runtime_error("\"SameAsConfigFile: true\" is not allowed here.");
ExtraParameters["SameAsConfigFile"] = __;
if (__)
{
Filename = config_file;
Format = "yaml";
return;
}
}
Filename = node["Filename"].as<std::string>();
Format = node["Format"].as<std::string>();
if (!supported_format.contains(Format))
throw std::runtime_error(fmt::format("Unsupported format: \"{}\"", Format));
if (auto _ = node["RelativeToConfigFile"])
{
auto __ = _.as<bool>();
ExtraParameters["RelativeToConfigFile"] = __;
if (__)
Filename = std::filesystem::path(config_file).parent_path() / Filename;
}
};
}

594
src/unfold.cpp
View File

@@ -1,384 +1,245 @@
# include <ufo/unfold.hpp>
# include <ufo.hpp>
# include <thread>
# include <syncstream>
# include <execution>
namespace ufo
void ufo::unfold(std::string config_file)
{
UnfoldSolver::InputType::InputType(std::string filename)
// 反折叠的原理: 将超胞中的原子运动状态, 投影到一组平面波构成的基矢中.
// 每一个平面波的波矢由两部分相加得到: 一部分是单胞倒格子的整数倍, 所取的个数有一定任意性, 论文中建议取大约单胞中原子个数那么多个;
// 对于没有缺陷的情况, 取一个应该就足够了.
// 这些平面波以原胞为周期。
// 另一部分是超胞倒格子的整数倍, 取 n 个, n 为超胞对应的单胞的倍数, 其实也就是倒空间中单胞对应倒格子中超胞的格点.
// 只要第一部分取得足够多, 那么单胞中原子的状态就可以完全被这些平面波描述.
// 将超胞中原子的运动状态投影到这些基矢上, 计算出投影的系数, 就可以将超胞的原子运动状态分解到单胞中的多个 q 点上.
struct Input
{
// read main input file
// 单胞的三个格矢，每行表示一个格矢的坐标，单位为埃
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
{
auto node = YAML::LoadFile(filename);
for (unsigned i = 0; i < 3; i++)
for (unsigned j = 0; j < 3; j++)
PrimativeCell(i, j) = node["PrimativeCell"][i][j].as<double>();
std::string Filename;
for (unsigned i = 0; i < 3; i++)
SuperCellMultiplier(i) = node["SuperCellMultiplier"][i].as<int>();
// 如果指定，则将结果投影到那些原子上
std::optional<std::vector<std::size_t>> SelectedAtoms;
if (auto value = node["SuperCellDeformation"])
// 默认输出为 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++)
{
SuperCellDeformation.emplace();
for (unsigned i = 0; i < 3; i++)
for (unsigned j = 0; j < 3; j++)
(*SuperCellDeformation)(i, j) = value[i][j].as<double>();
}
for (unsigned i = 0; i < 3; i++)
PrimativeCellBasisNumber(i) = node["PrimativeCellBasisNumber"][i].as<int>();
AtomPositionInputFile = DataFile
(
node["AtomPositionInputFile"], {"yaml"},
filename, true
);
QpointDataInputFile = DataFile
(
node["QpointDataInputFile"], {"yaml", "hdf5"},
filename, true
);
if (auto value = node["QpointDataOutputFile"])
{
QpointDataOutputFile.resize(value.size());
for (unsigned i = 0; i < value.size(); i++)
QpointDataOutputFile[i] = DataFile
(
value[i], {"yaml", "yaml-human-readable", "zpp", "hdf5"},
filename, false
);
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;
};
if (AtomPositionInputFile.Format == "yaml")
{
auto node = YAML::LoadFile(AtomPositionInputFile.Filename);
std::vector<YAML::Node> points;
if (auto _ = node["points"])
points = _.as<std::vector<YAML::Node>>();
else
points = node["unit_cell"]["points"].as<std::vector<YAML::Node>>();
auto atom_position_to_super_cell = Eigen::MatrixX3d(points.size(), 3);
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>();
auto super_cell = (SuperCellDeformation.value_or(Eigen::Matrix3d::Identity())
* SuperCellMultiplier.cast<double>().asDiagonal() * PrimativeCell).eval();
AtomPosition = atom_position_to_super_cell * super_cell;
}
if (QpointDataInputFile.Format == "yaml")
{
auto node = YAML::LoadFile(QpointDataInputFile.Filename);
auto phonon = node["phonon"].as<std::vector<YAML::Node>>();
QpointData.resize(phonon.size());
for (unsigned i = 0; i < phonon.size(); i++)
{
for (unsigned j = 0; j < 3; j++)
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 (unsigned j = 0; j < band.size(); j++)
{
QpointData[i].ModeData[j].Frequency = band[j]["frequency"].as<double>();
auto eigenvector_vectors = band[j]["eigenvector"]
.as<std::vector<std::vector<std::vector<double>>>>();
Eigen::MatrixX3cd eigenvectors(AtomPosition.rows(), 3);
for (unsigned k = 0; k < AtomPosition.rows(); k++)
for (unsigned l = 0; l < 3; l++)
eigenvectors(k, l)
= eigenvector_vectors[k][l][0] + 1i * eigenvector_vectors[k][l][1];
// 需要对读入的原子运动状态作相位转换, 使得它们与我们的约定一致(对超胞周期性重复)
// 这里还要需要做归一化处理 (指将数据简单地作为向量处理的归一化)
auto& AtomMovement = QpointData[i].ModeData[j].AtomMovement;
// AtomMovement = eigenvectors.array().colwise() * (-2 * std::numbers::pi_v<double> * 1i
// * (atom_position_to_super_cell * input.QpointData[i].Qpoint)).array().exp();
// AtomMovement /= AtomMovement.norm();
// phonopy 似乎已经进行了相位的转换！为什么？
AtomMovement = eigenvectors / eigenvectors.norm();
}
}
}
else if (QpointDataInputFile.Format == "hdf5")
{
std::vector<std::vector<std::vector<double>>> frequency, path;
std::vector<std::vector<std::vector<std::vector<PhonopyComplex>>>> eigenvector_vector;
Hdf5file{}.open_for_read(QpointDataInputFile.Filename).read(frequency, "/frequency")
.read(eigenvector_vector, "/eigenvector")
.read(path, "/path");
std::vector size = { frequency.size(), frequency[0].size(), frequency[0][0].size() };
QpointData.resize(size[0] * size[1]);
for (unsigned i = 0; i < size[0]; i++)
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].ModeData.resize(size[2]);
for (unsigned k = 0; k < size[2]; k++)
{
QpointData[i * size[1] + j].ModeData[k].Frequency = frequency[i][j][k];
Eigen::MatrixX3cd eigenvectors(AtomPosition.rows(), 3);
for (unsigned l = 0; l < AtomPosition.rows(); l++)
for (unsigned m = 0; m < 3; m++)
eigenvectors(l, m)
= 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();
}
}
}
}
void UnfoldSolver::OutputType::write
(decltype(InputType::QpointDataOutputFile) output_files) const
{
for (auto& output_file : output_files)
write(output_file.Filename, output_file.Format);
}
void UnfoldSolver::OutputType::write(std::string filename, std::string format, unsigned percision) const
{
if (format == "yaml")
std::ofstream(filename) << [&]
{
std::stringstream print;
print << "QpointData:\n";
for (auto& qpoint: QpointData)
{
print << fmt::format(" - Qpoint: [ {1:.{0}f}, {2:.{0}f}, {3:.{0}f} ]\n",
percision, qpoint.Qpoint[0], qpoint.Qpoint[1], qpoint.Qpoint[2]);
print << fmt::format(" Source: [ {1:.{0}f}, {2:.{0}f}, {3:.{0}f} ]\n",
percision, qpoint.Source[0], qpoint.Source[1], qpoint.Source[2]);
print << " ModeData:\n";
for (auto& mode: qpoint.ModeData)
print << fmt::format(" - {{ Frequency: {1:.{0}f}, Weight: {2:.{0}f} }}\n",
percision, mode.Frequency, mode.Weight);
}
return print.str();
}();
else if (format == "yaml-human-readable")
{
std::remove_cvref_t<decltype(*this)> output;
std::map<unsigned, std::vector<decltype(QpointData)::const_iterator>>
meta_qpoint_to_sub_qpoint_iterators;
for (auto it = QpointData.begin(); it != QpointData.end(); 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& qpoint : sub_qpoint_iterators)
{
std::map<double, double> frequency_to_weight;
for (unsigned i_of_mode = 0; i_of_mode < qpoint->ModeData.size(); i_of_mode++)
{
auto frequency = qpoint->ModeData[i_of_mode].Frequency;
auto weight = qpoint->ModeData[i_of_mode].Weight;
auto it_lower = frequency_to_weight.lower_bound(frequency - 0.1);
auto it_upper = frequency_to_weight.upper_bound(frequency + 0.1);
if (it_lower == it_upper)
frequency_to_weight[frequency] = weight;
else
{
auto frequency_sum = std::accumulate(it_lower, it_upper, 0.,
[](const auto& a, const auto& b) { return a + b.first * b.second; });
auto weight_sum = std::accumulate(it_lower, it_upper, 0.,
[](const auto& a, const auto& b) { return a + b.second; });
frequency_sum += frequency * weight;
weight_sum += weight;
frequency_to_weight.erase(it_lower, it_upper);
frequency_to_weight[frequency_sum / weight_sum] = weight_sum;
}
}
auto& _ = output.QpointData.emplace_back();
_.Qpoint = qpoint->Qpoint;
_.Source = qpoint->Source;
_.SourceIndex_ = qpoint->SourceIndex_;
for (auto [frequency, weight] : frequency_to_weight)
if (weight > 0.1)
{
auto& __ = _.ModeData.emplace_back();
__.Frequency = frequency;
__.Weight = weight;
}
}
output.write(filename, "yaml", 3);
}
else if (format == "zpp")
zpp_write(*this, filename);
else if (format == "hdf5")
{
std::vector<std::vector<double>> Qpoint, Source, Frequency, Weight;
for (auto& qpoint : QpointData)
{
Qpoint.emplace_back(qpoint.Qpoint.data(), qpoint.Qpoint.data() + 3);
Source.emplace_back(qpoint.Source.data(), qpoint.Source.data() + 3);
Frequency.emplace_back();
Weight.emplace_back();
for (auto& mode : qpoint.ModeData)
{
Frequency.back().push_back(mode.Frequency);
Weight.back().push_back(mode.Weight);
}
}
Hdf5file{}.open_for_write(filename).write(Qpoint, "/Qpoint")
.write(Source, "/Source")
.write(Frequency, "/Frequency")
.write(Weight, "/Weight");
}
}
UnfoldSolver::UnfoldSolver(std::string config_file) : Input_([&]
{
std::clog << "Reading input file... " << std::flush;
return config_file;
}())
{
std::clog << "Done." << std::endl;
}
UnfoldSolver& UnfoldSolver::operator()()
{
if (!Basis_)
{
std::clog << "Constructing basis... " << std::flush;
Basis_ = construct_basis
(
Input_.PrimativeCell, Input_.SuperCellMultiplier,
Input_.PrimativeCellBasisNumber, Input_.AtomPosition
);
std::clog << "Done." << std::endl;
}
if (!Output_)
{
std::clog << "Calculating projection coefficient... " << std::flush;
std::vector<std::reference_wrapper<const decltype
(InputType::QpointDataType::ModeDataType::AtomMovement)>> mode_data;
for (auto& qpoint : Input_.QpointData)
for (auto& mode : qpoint.ModeData)
mode_data.emplace_back(mode.AtomMovement);
std::atomic<unsigned> number_of_finished_modes(0);
std::thread print_thread([&]
{
unsigned n;
while ((n = number_of_finished_modes) < mode_data.size())
{
std::osyncstream(std::cerr) << fmt::format("\rCalculating projection coefficient... ({}/{})",
number_of_finished_modes, mode_data.size()) << std::flush;
std::this_thread::sleep_for(100ms);
number_of_finished_modes.wait(n);
}
});
auto projection_coefficient = construct_projection_coefficient
(*Basis_, mode_data, number_of_finished_modes);
number_of_finished_modes = mode_data.size();
print_thread.join();
std::clog << "\33[2K\rCalculating projection coefficient... Done." << std::endl;
std::clog << "Constructing output... " << std::flush;
std::vector<std::reference_wrapper<const decltype(InputType::QpointDataType::Qpoint)>> qpoint;
std::vector<std::vector<std::reference_wrapper<const
decltype(InputType::QpointDataType::ModeDataType::Frequency)>>> frequency;
for (auto& qpoint_data : Input_.QpointData)
{
qpoint.emplace_back(qpoint_data.Qpoint);
frequency.emplace_back();
for (auto& mode_data : qpoint_data.ModeData)
frequency.back().emplace_back(mode_data.Frequency);
}
Output_ = construct_output
(
Input_.SuperCellMultiplier,
Input_.SuperCellDeformation, qpoint, frequency, projection_coefficient
);
std::clog << "Done." << std::endl;
}
std::clog << "Writing output... " << std::flush;
Output_->write(Input_.QpointDataOutputFile);
std::clog << "Done." << std::endl;
return *this;
}
UnfoldSolver::BasisType UnfoldSolver::construct_basis
// 构建基
// 每个 q 点对应一组 sub qpoint。不同的 q 点所对应的 sub qpoint 是不一样的，但 sub qpoint 与 q 点的相对位移在不同 q 点之间是相同的。
// 由于基只与这个相对位置有关（也就是说，不同 q 点的基是一样的），因此可以先计算出所有的基，这样降低计算量。
// 外层下标对应超胞倒格子的整数倍那部分(第二部分), 也就是不同的 sub qpoint
// 内层下标对应单胞倒格子的整数倍那部分(第一部分), 也就是 sub qpoint 上的不同平面波（取的数量越多，结果越精确）
auto construct_basis = []
(
const decltype(InputType::PrimativeCell)& primative_cell,
const decltype(InputType::SuperCellMultiplier)& super_cell_multiplier,
const decltype(InputType::PrimativeCellBasisNumber)& primative_cell_basis_number,
const decltype(InputType::AtomPosition)& atom_position
Eigen::Matrix3d primative_cell, Eigen::Vector3i super_cell_multiplier,
Eigen::Vector<std::size_t, 3> primative_cell_basis_number, Eigen::MatrixX3d atom_position
)
{
BasisType basis(super_cell_multiplier.prod());
// 每个 q 点对应的一组 sub qpoint。不同的 q 点所对应的 sub qpoint 是不一样的，但 sub qpoint 与 q 点的相对位置一致。
// 这里 xyz_of_diff_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))
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] : triplet_sequence(primative_cell_basis_number))
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()
* xyz_of_diff_of_sub_qpoint_by_reciprocal_modified_super_cell.cast<double>();
// 将 q 点坐标转换为埃^-1
auto qpoint = (diff_of_sub_qpoint_by_reciprocal_primative_cell.transpose()
* 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 * qpoint)).array().exp();
= (2i * std::numbers::pi_v<double> * (atom_position * diff_of_sub_qpoint)).array().exp();
}
}
return basis;
}
};
std::vector<std::vector<double>> UnfoldSolver::construct_projection_coefficient
// 计算从超胞到原胞的投影系数（不是分原子的投影系数），是反折叠的核心步骤
// 返回的投影系数是一个三维数组，第一维对应不同的 q 点，第二维对应不同的模式，第三维对应不同的 sub qpoint
auto construct_projection_coefficient = []
(
const BasisType& basis,
const std::vector<std::reference_wrapper<const decltype
(InputType::QpointDataType::ModeDataType::AtomMovement)>>& mode_data,
std::atomic<unsigned>& number_of_finished_modes
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, mode_data.begin(), mode_data.end(),
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 (unsigned i_of_sub_qpoint = 0; i_of_sub_qpoint < basis.size(); i_of_sub_qpoint++)
for (std::size_t i_of_sub_qpoint = 0; i_of_sub_qpoint < basis.size(); i_of_sub_qpoint++)
// 对于 basis 中, 对应于单胞倒格子的部分, 以及对应于不同方向的部分, 分别求内积, 然后求模方和
for (unsigned i_of_basis = 0; i_of_basis < basis[i_of_sub_qpoint].size(); i_of_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 = std::accumulate
(projection_coefficient.begin(), projection_coefficient.end(), 0.);
for (auto& _ : projection_coefficient)
_ /= sum;
auto sum = ranges::accumulate(projection_coefficient, 0.);
for (auto& _ : projection_coefficient) _ /= sum;
number_of_finished_modes++;
return projection_coefficient;
}
);
return projection_coefficient;
}
UnfoldSolver::OutputType UnfoldSolver::construct_output
(
const decltype(InputType::SuperCellMultiplier)& super_cell_multiplier,
const decltype(InputType::SuperCellDeformation)& super_cell_deformation,
const std::vector<std::reference_wrapper<const decltype
(InputType::QpointDataType::Qpoint)>>& meta_qpoint_by_reciprocal_super_cell,
const std::vector<std::vector<std::reference_wrapper<const decltype
(InputType::QpointDataType::ModeDataType::Frequency)>>>& frequency,
const ProjectionCoefficientType_& projection_coefficient
)
{
OutputType output;
// 将计算得到的投影系数重新组装成三维数组
// 第一维是 meta qpoint，第二维是模式，第三维是 sub qpoint
std::vector<std::vector<std::vector<double>>> projection_coefficient_output;
for
(
unsigned i_of_meta_qpoint = 0, num_of_mode_manipulated = 0;
i_of_meta_qpoint < meta_qpoint_by_reciprocal_super_cell.size();
i_of_meta_qpoint++
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++)
{
for (auto [xyz_of_diff_of_sub_qpoint_by_reciprocal_modified_super_cell, i_of_sub_qpoint]
: triplet_sequence(super_cell_multiplier))
// 如果需要投影到特定的原子上，需要先计算当前 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();
/*
@@ -404,26 +265,81 @@ namespace ufo
*/
auto sub_qpoint_by_reciprocal_primative_cell =
(
super_cell_multiplier.cast<double>().cwiseInverse().asDiagonal()
input.SuperCellMultiplier.cast<double>().cwiseInverse().asDiagonal()
* (
xyz_of_diff_of_sub_qpoint_by_reciprocal_modified_super_cell.cast<double>()
+ super_cell_deformation.value_or(Eigen::Matrix3d::Identity()).inverse()
* meta_qpoint_by_reciprocal_super_cell[i_of_meta_qpoint].get().cast<double>()
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 = meta_qpoint_by_reciprocal_super_cell[i_of_meta_qpoint];
_.SourceIndex_ = i_of_meta_qpoint;
for (unsigned i_of_mode = 0; i_of_mode < frequency[i_of_meta_qpoint].size(); i_of_mode++)
_.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 = frequency[i_of_meta_qpoint][i_of_mode];
__.Weight = projection_coefficient[num_of_mode_manipulated + i_of_mode][i_of_sub_qpoint];
__.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];
}
}
num_of_mode_manipulated += frequency[i_of_meta_qpoint].size();
}
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;
}

12
test/14.2.2.4.unfold.yaml
View File

@@ -1,12 +0,0 @@
PrimativeCell:
- [ 1.548010265167714, -2.681232429908652, 0.000000000000000 ] # a
- [ 1.548010265167714, 2.681232429908652, 0.000000000000000 ] # b
- [ 0.000000000000000, 0.000000000000000, 5.061224103556595 ] # c
SuperCellMultiplier: [ 3, 1, 1 ]
PrimativeCellBasisNumber: [ 8, 8, 8 ]
AtomPositionInputFile: { FileName: "/home/chn/Documents/lammps-SiC/14/14.2/14.2.2/14.2.2.4/band.yaml", Format: "yaml" }
QpointDataInputFile: { FileName: "/home/chn/Documents/lammps-SiC/14/14.2/14.2.2/14.2.2.4/band.yaml", Format: "yaml" }
QpointDataOutputFile:
- { FileName: "test/14.2.2.4.result.yaml", Format: "yaml" }
- { FileName: "test/14.2.2.4.result.human-readable.yaml", Format: "yaml-human-readable" }
- { FileName: "test/14.2.2.4.result.zpp", Format: "zpp" }

14
test/14.2.2.5.plot.yaml
View File

@@ -1,14 +0,0 @@
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.2.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.2.5.png

12
test/14.2.2.5.unfold.yaml
View File

@@ -1,12 +0,0 @@
PrimativeCell:
- [ 1.548010265167714, -2.681232429908652, 0.000000000000000 ] # a
- [ 1.548010265167714, 2.681232429908652, 0.000000000000000 ] # b
- [ 0.000000000000000, 0.000000000000000, 5.061224103556595 ] # c
SuperCellMultiplier: [ 3, 1, 1 ]
PrimativeCellBasisNumber: [ 8, 8, 8 ]
AtomPositionInputFile: { FileName: "/home/chn/Documents/lammps-SiC/14/14.2/14.2.2/14.2.2.5/band.yaml", Format: "yaml" }
QpointDataInputFile: { FileName: "/home/chn/Documents/lammps-SiC/14/14.2/14.2.2/14.2.2.5/band.yaml", Format: "yaml" }
QpointDataOutputFile:
- { FileName: "test/14.2.2.5.result.yaml", Format: "yaml" }
- { FileName: "test/14.2.2.5.result.human-readable.yaml", Format: "yaml-human-readable" }
- { FileName: "test/14.2.2.5.result.zpp", Format: "zpp" }

14
test/14.2.4.4.plot.yaml
View File

@@ -1,14 +0,0 @@
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.4.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.4.4.png

14
test/14.2.4.4.unfold.yaml
View File

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

14
test/14.2.5.5.plot.yaml
View File

@@ -1,14 +0,0 @@
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
View File

@@ -1,18 +0,0 @@
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" }

15
test/14.2.6.4.plot.yaml
View File

@@ -1,15 +0,0 @@
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, 1024 ]
Range: [ -5, 35 ]
Filename: ./test/14.2.6.4.png
YTicks: [ 0, 5, 10, 15, 20, 25, 30 ]

18
test/14.2.6.4.unfold.yaml
View File

@@ -1,18 +0,0 @@
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" }

18
test/fold/config.yaml Normal file
View File

@@ -0,0 +1,18 @@
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