mirror of
https://github.com/CHN-beta/nixos.git
synced 2026-01-12 07:29:23 +08:00
packages.ufo: 重构
This commit is contained in:
@@ -186,7 +186,7 @@
|
||||
};
|
||||
ufo = pkgs.mkShell.override { stdenv = pkgs.clang18Stdenv; }
|
||||
{
|
||||
inputsFrom = [ (pkgs.localPackages.ufo.override { version = null; }) ];
|
||||
inputsFrom = [ pkgs.localPackages.ufo ];
|
||||
packages = [ pkgs.clang-tools_18 ];
|
||||
CMAKE_EXPORT_COMPILE_COMMANDS = "1";
|
||||
};
|
||||
|
||||
@@ -26,13 +26,15 @@ find_library(LIBBACKTRACE_LIBRARY NAMES backtrace REQUIRED)
|
||||
find_package(HDF5 REQUIRED)
|
||||
find_package(concurrencpp REQUIRED)
|
||||
find_path(POCKETFFT_INCLUDE_DIR pocketfft.h REQUIRED)
|
||||
find_package(yaml-cpp REQUIRED)
|
||||
|
||||
add_library(biu src/common.cpp src/hdf5.cpp src/logger.cpp src/string.cpp)
|
||||
target_include_directories(biu PUBLIC $<BUILD_INTERFACE:${CMAKE_CURRENT_SOURCE_DIR}/include>
|
||||
$<INSTALL_INTERFACE:${CMAKE_INSTALL_INCLUDEDIR}> ${NAMEOF_INCLUDE_DIR} ${ZPP_BITS_INCLUDE_DIR}
|
||||
${LIBBACKTRACE_INCLUDE_DIR} ${POCKETFFT_INCLUDE_DIR})
|
||||
target_link_libraries(biu PUBLIC magic_enum::magic_enum fmt::fmt Boost::headers Boost::iostreams Boost::filesystem
|
||||
range-v3::range-v3 Eigen3::Eigen HighFive TgBot::TgBot ${LIBBACKTRACE_LIBRARY} hdf5::hdf5 concurrencpp::concurrencpp)
|
||||
range-v3::range-v3 Eigen3::Eigen HighFive TgBot::TgBot ${LIBBACKTRACE_LIBRARY} hdf5::hdf5 concurrencpp::concurrencpp
|
||||
yaml-cpp::yaml-cpp)
|
||||
target_compile_features(biu PUBLIC cxx_std_23)
|
||||
target_compile_options(biu PUBLIC -Wno-gnu-string-literal-operator-template)
|
||||
install(TARGETS biu EXPORT biuTargets LIBRARY DESTINATION ${CMAKE_INSTALL_LIBDIR}
|
||||
@@ -75,3 +77,7 @@ add_executable(test-hdf5 test/hdf5.cpp)
|
||||
target_link_libraries(test-hdf5 PRIVATE biu)
|
||||
set_property(TARGET test-hdf5 PROPERTY CXX_STANDARD 23 CXX_STANDARD_REQUIRED ON)
|
||||
add_test(NAME test-hdf5 COMMAND test-hdf5)
|
||||
add_executable(test-yaml test/yaml.cpp)
|
||||
target_link_libraries(test-yaml PRIVATE biu)
|
||||
set_property(TARGET test-yaml PROPERTY CXX_STANDARD 23 CXX_STANDARD_REQUIRED ON)
|
||||
add_test(NAME test-yaml COMMAND test-yaml)
|
||||
|
||||
@@ -12,3 +12,5 @@ find_path(LIBBACKTRACE_INCLUDE_DIR backtrace.h REQUIRED)
|
||||
find_library(LIBBACKTRACE_LIBRARY NAMES backtrace REQUIRED)
|
||||
find_package(HDF5 REQUIRED)
|
||||
find_package(concurrencpp REQUIRED)
|
||||
find_path(POCKETFFT_INCLUDE_DIR pocketfft.h REQUIRED)
|
||||
find_package(yaml-cpp REQUIRED)
|
||||
|
||||
@@ -1,7 +1,7 @@
|
||||
{
|
||||
stdenv, cmake, lib,
|
||||
magic-enum, fmt, boost, eigen, range-v3, nameof, zpp-bits, highfive, tgbot-cpp, libbacktrace, hdf5, concurrencpp,
|
||||
pocketfft
|
||||
pocketfft, yaml-cpp
|
||||
}: stdenv.mkDerivation rec
|
||||
{
|
||||
name = "biu";
|
||||
@@ -9,7 +9,7 @@
|
||||
buildInputs =
|
||||
[
|
||||
magic-enum fmt boost range-v3 nameof zpp-bits eigen highfive tgbot-cpp libbacktrace hdf5
|
||||
concurrencpp pocketfft
|
||||
concurrencpp pocketfft yaml-cpp
|
||||
];
|
||||
propagatedBuildInputs = buildInputs;
|
||||
nativeBuildInputs = [ cmake ];
|
||||
|
||||
@@ -10,5 +10,5 @@
|
||||
# include <biu/logger.tpp>
|
||||
# include <biu/smartref.tpp>
|
||||
# include <biu/fft.tpp>
|
||||
|
||||
# include <biu/yaml.tpp>
|
||||
# include <range/v3/all.hpp>
|
||||
|
||||
@@ -1,8 +1,6 @@
|
||||
# pragma once
|
||||
# include <mutex>
|
||||
# include <optional>
|
||||
# include <condition_variable>
|
||||
# include <cstddef>
|
||||
# include <experimental/memory>
|
||||
# include <biu/common.hpp>
|
||||
# include <biu/concepts.hpp>
|
||||
|
||||
@@ -108,17 +108,19 @@ namespace biu
|
||||
|
||||
namespace detail_
|
||||
{
|
||||
template <typename Byte> struct ReadFileReturnType;
|
||||
template <> struct ReadFileReturnType<std::byte> { using Type = std::vector<std::byte>; };
|
||||
template <> struct ReadFileReturnType<char> { using Type = std::string; };
|
||||
template <typename Byte> struct ReadReturnType;
|
||||
template <> struct ReadReturnType<std::byte> { using Type = std::vector<std::byte>; };
|
||||
template <> struct ReadReturnType<char> { using Type = std::string; };
|
||||
}
|
||||
template <typename Byte = std::byte> detail_::ReadFileReturnType<Byte>::Type
|
||||
read_file(const std::filesystem::path& path);
|
||||
template<> std::vector<std::byte> read_file<std::byte>(const std::filesystem::path& path);
|
||||
template<> std::string read_file<char>(const std::filesystem::path& path);
|
||||
template <typename Byte> detail_::ReadReturnType<Byte>::Type read(const std::filesystem::path& path);
|
||||
template <typename Byte> detail_::ReadReturnType<Byte>::Type read(std::istream& input);
|
||||
template<> std::vector<std::byte> read<std::byte>(const std::filesystem::path& path);
|
||||
template<> std::string read<char>(const std::filesystem::path& path);
|
||||
template<> std::vector<std::byte> read<std::byte>(std::istream& input);
|
||||
template<> std::string read<char>(std::istream& input);
|
||||
}
|
||||
using common::hash, common::unused, common::block_forever, common::is_interactive, common::env, common::int128_t,
|
||||
common::uint128_t, common::Empty, common::CaseInsensitiveStringLessComparator, common::RemoveMemberPointer,
|
||||
common::MoveQualifiers, common::FallbackIfNoTypeDeclared, common::exec, common::serialize, common::deserialize,
|
||||
common::sequence, common::read_file;
|
||||
common::sequence, common::read;
|
||||
}
|
||||
|
||||
@@ -1,5 +1,4 @@
|
||||
# pragma once
|
||||
# include <fstream>
|
||||
# include <boost/functional/hash.hpp>
|
||||
# include <zpp_bits.h>
|
||||
# include <biu/common.hpp>
|
||||
|
||||
@@ -38,6 +38,7 @@ namespace biu
|
||||
template <std::size_t Row, std::size_t Col> struct ToEigenHelper {};
|
||||
template <std::size_t Size> struct FromEigenVectorHelper {};
|
||||
template <std::size_t Row, std::size_t Col> struct FromEigenMatrixHelper {};
|
||||
struct FromEigenHelper {};
|
||||
|
||||
// convert 1D standard container to Eigen::Matrix, the second argument should always be unspecified
|
||||
template <typename From, std::size_t ToSize> auto operator|
|
||||
@@ -61,6 +62,10 @@ namespace biu
|
||||
// convert 2D Eigen matrix to std::vector or std::array
|
||||
template <typename Matrix, std::size_t ToRow, std::size_t ToCol> auto operator|
|
||||
(const Matrix&, const detail_::FromEigenMatrixHelper<ToRow, ToCol>&);
|
||||
|
||||
// 自动选择 fromEigenVector 还是 fromEigenMatrix
|
||||
// 如果行数或者列数 <= 1，那么选择 fromEigenVector，否则选择 fromEigenMatrix
|
||||
template <typename Matrix> auto operator|(const Matrix&, const detail_::FromEigenHelper&);
|
||||
}
|
||||
|
||||
// usage: some_value | toEigen<Row, Col>
|
||||
@@ -70,6 +75,7 @@ namespace biu
|
||||
inline constexpr detail_::FromEigenVectorHelper<Size> fromEigenVector;
|
||||
template <std::size_t Row = detail_::unspecifiedSize, std::size_t Col = detail_::unspecifiedSize>
|
||||
inline constexpr detail_::FromEigenMatrixHelper<Row, Col> fromEigenMatrix;
|
||||
inline constexpr detail_::FromEigenHelper fromEigen;
|
||||
|
||||
// test if a class is an eigen matrix
|
||||
namespace detail_
|
||||
@@ -80,7 +86,7 @@ namespace biu
|
||||
}
|
||||
template <typename Matrix> concept EigenMatrix = detail_::EigenMatrix<Matrix>::value;
|
||||
}
|
||||
using eigen::toEigen, eigen::fromEigenVector, eigen::fromEigenMatrix, eigen::EigenMatrix;
|
||||
using eigen::toEigen, eigen::fromEigenVector, eigen::fromEigenMatrix, eigen::fromEigen, eigen::EigenMatrix;
|
||||
}
|
||||
|
||||
// archive a matrix
|
||||
|
||||
@@ -254,11 +254,25 @@ namespace biu::eigen
|
||||
for (int i = 0; i < matrix.rows(); i++)
|
||||
{
|
||||
using RowVector = Eigen::RowVector
|
||||
<typename Matrix::Scalar, size.first.second == dynamicSize ? Eigen::Dynamic : size.first.second>;
|
||||
<
|
||||
typename Matrix::Scalar,
|
||||
size.first.second == dynamicSize ? Eigen::Dynamic
|
||||
: static_cast<decltype(Eigen::Dynamic)>(size.first.second)
|
||||
>;
|
||||
Eigen::Map<RowVector>(result[i].data(), 1, matrix.cols()) = matrix.row(i);
|
||||
}
|
||||
return result;
|
||||
}
|
||||
template <typename Matrix> auto detail_::operator|(const Matrix& matrix, const detail_::FromEigenHelper&)
|
||||
{
|
||||
constexpr auto
|
||||
ncols = Matrix::CompileTimeTraits::ColsAtCompileTime, nrows = Matrix::CompileTimeTraits::RowsAtCompileTime;
|
||||
if constexpr ((ncols <= 1 && ncols != Eigen::Dynamic) || (nrows <= 1 && nrows != Eigen::Dynamic))
|
||||
return matrix | fromEigenVector<>;
|
||||
else
|
||||
return matrix | fromEigenMatrix<>;
|
||||
}
|
||||
|
||||
}
|
||||
template <typename Matrix> constexpr auto Eigen::serialize(auto & archive, Matrix& matrix)
|
||||
requires biu::EigenMatrix<std::remove_cvref_t<Matrix>>
|
||||
|
||||
@@ -4,6 +4,8 @@
|
||||
# include <experimental/memory>
|
||||
# include <biu/string.hpp>
|
||||
# include <biu/concepts.hpp>
|
||||
# include <fmt/format.h>
|
||||
# include <fmt/ostream.h>
|
||||
|
||||
namespace biu
|
||||
{
|
||||
|
||||
@@ -11,6 +11,7 @@ namespace biu
|
||||
protected: T& Ref_;
|
||||
|
||||
public: SmartRef(T& val);
|
||||
public: SmartRef(T&& val);
|
||||
public: template <typename... Us> requires (std::is_constructible_v<T, Us...>) SmartRef(Us&&... val);
|
||||
public: T& operator*();
|
||||
public: T* operator->();
|
||||
|
||||
@@ -5,6 +5,8 @@
|
||||
namespace biu
|
||||
{
|
||||
template <typename T> inline SmartRef<T>::SmartRef(T& val) : Ref_(val) {}
|
||||
template <typename T> inline SmartRef<T>::SmartRef(T&& val)
|
||||
: Ptr_(std::make_unique<T>(std::move(val))), Ref_(*Ptr_) {}
|
||||
template <typename T> template <typename... Us> requires (std::is_constructible_v<T, Us...>)
|
||||
inline SmartRef<T>::SmartRef(Us&&... val)
|
||||
: Ptr_(std::make_unique<T>(std::forward<Us>(val)...)), Ref_(*Ptr_) {}
|
||||
|
||||
32
packages/biu/include/biu/yaml.hpp
Normal file
32
packages/biu/include/biu/yaml.hpp
Normal file
@@ -0,0 +1,32 @@
|
||||
# pragma once
|
||||
# include <yaml-cpp/yaml.h>
|
||||
# include <biu/eigen.hpp>
|
||||
|
||||
namespace biu
|
||||
{
|
||||
struct YamlParsable {};
|
||||
}
|
||||
|
||||
namespace YAML
|
||||
{
|
||||
template <biu::EigenMatrix Matrix> struct convert<Matrix>
|
||||
{
|
||||
static Node encode(const Matrix&);
|
||||
static bool decode(const Node& node, Matrix&);
|
||||
};
|
||||
template <biu::SpecializationOf<std::complex> Complex> struct convert<Complex>
|
||||
{
|
||||
static Node encode(const Complex&);
|
||||
static bool decode(const Node& node, Complex&);
|
||||
};
|
||||
template <biu::SpecializationOf<std::optional> Optional> struct convert<Optional>
|
||||
{
|
||||
static Node encode(const Optional&);
|
||||
static bool decode(const Node& node, Optional&);
|
||||
};
|
||||
template <typename T> struct convert
|
||||
{
|
||||
static Node encode(const T&);
|
||||
static bool decode(const Node& node, T&);
|
||||
};
|
||||
}
|
||||
75
packages/biu/include/biu/yaml.tpp
Normal file
75
packages/biu/include/biu/yaml.tpp
Normal file
@@ -0,0 +1,75 @@
|
||||
# pragma once
|
||||
# include <biu/yaml.hpp>
|
||||
# include <biu/concepts.hpp>
|
||||
# include <biu/format.hpp>
|
||||
# include <biu/eigen.hpp>
|
||||
# include <boost/pfr.hpp>
|
||||
# include <boost/pfr/core_name.hpp>
|
||||
|
||||
namespace YAML
|
||||
{
|
||||
template <biu::EigenMatrix Matrix> Node convert<Matrix>::encode(const Matrix& matrix)
|
||||
{
|
||||
auto std_matrix = matrix | biu::fromEigen;
|
||||
return convert<decltype(std_matrix)>::encode(std_matrix);
|
||||
}
|
||||
template <biu::EigenMatrix Matrix> bool convert<Matrix>::decode(const Node& node, Matrix& matrix)
|
||||
{
|
||||
using std_matrix = decltype(matrix | biu::fromEigen);
|
||||
std_matrix value;
|
||||
if (!convert<std_matrix>::decode(node, value)) return false;
|
||||
matrix = value | biu::toEigen<>;
|
||||
return true;
|
||||
}
|
||||
template <biu::SpecializationOf<std::complex> Complex> Node convert<Complex>::encode(const Complex& complex)
|
||||
{
|
||||
return convert<std::array<typename Complex::value_type, 2>>::encode({ complex.real(), complex.imag() });
|
||||
}
|
||||
template <biu::SpecializationOf<std::complex> Complex> bool convert<Complex>::decode
|
||||
(const Node& node, Complex& complex)
|
||||
{
|
||||
std::array<typename Complex::value_type, 2> arr;
|
||||
if (!convert<std::array<typename Complex::value_type, 2>>::decode(node, arr)) return false;
|
||||
complex = Complex{ arr[0], arr[1] };
|
||||
return true;
|
||||
}
|
||||
template <biu::SpecializationOf<std::optional> Optional> Node convert<Optional>::encode(const Optional& optional)
|
||||
{
|
||||
if (optional) return convert<typename Optional::value_type>::encode(*optional);
|
||||
else return YAML::Node{};
|
||||
}
|
||||
template <biu::SpecializationOf<std::optional> Optional> bool convert<Optional>::decode
|
||||
(const Node& node, Optional& optional)
|
||||
{
|
||||
if (node.IsNull() || !node.IsDefined()) optional = std::nullopt;
|
||||
else
|
||||
{
|
||||
typename Optional::value_type value;
|
||||
if (!convert<typename Optional::value_type>::decode(node, value)) return false;
|
||||
optional = value;
|
||||
}
|
||||
return true;
|
||||
}
|
||||
template <typename T> Node convert<T>::encode(const T& t)
|
||||
{
|
||||
YAML::Node node;
|
||||
boost::pfr::for_each_field(t, [&](const auto& field, auto index)
|
||||
{
|
||||
using type = std::remove_cvref_t<decltype(field)>;
|
||||
auto name = boost::pfr::get_name<decltype(index)::value, T>();
|
||||
node[name] = convert<type>::encode(field);
|
||||
});
|
||||
return node;
|
||||
}
|
||||
template <typename T> bool convert<T>::decode(const Node& node, T& t)
|
||||
{
|
||||
bool result = true;
|
||||
boost::pfr::for_each_field(t, [&](auto& field, auto index)
|
||||
{
|
||||
using type = std::remove_cvref_t<decltype(field)>;
|
||||
auto name = boost::pfr::get_name<decltype(index)::value, T>();
|
||||
result = result && convert<type>::decode(node[name], field);
|
||||
});
|
||||
return result;
|
||||
}
|
||||
}
|
||||
@@ -67,7 +67,7 @@ namespace biu
|
||||
template detail_::ExecResult<detail_::ExecMode##i> \
|
||||
exec<detail_::ExecMode##i>(detail_::ExecInput<detail_::ExecMode##i>);
|
||||
BOOST_PP_FOR(0, BIU_EXEC_PRED, BIU_EXEC_OP, BIU_EXEC_MACRO)
|
||||
template<> std::vector<std::byte> read_file<std::byte>(const std::filesystem::path& path)
|
||||
template<> std::vector<std::byte> read<std::byte>(const std::filesystem::path& path)
|
||||
{
|
||||
auto length = std::filesystem::file_size(path);
|
||||
std::vector<std::byte> buffer(length);
|
||||
@@ -75,10 +75,20 @@ namespace biu
|
||||
in.read(reinterpret_cast<char*>(buffer.data()), length);
|
||||
return buffer;
|
||||
}
|
||||
template<> std::string read_file<char>(const std::filesystem::path& path)
|
||||
template<> std::string read<char>(const std::filesystem::path& path)
|
||||
{
|
||||
auto buffer = read_file<std::byte>(path);
|
||||
auto buffer = read<std::byte>(path);
|
||||
return std::string{reinterpret_cast<char*>(buffer.data()), buffer.size()};
|
||||
}
|
||||
template<> std::vector<std::byte> read<std::byte>(std::istream& input)
|
||||
{
|
||||
static_assert(sizeof(std::byte) == sizeof(char));
|
||||
auto buffer = read<char>(input);
|
||||
return buffer
|
||||
| ranges::views::transform([](char c){ return std::byte(static_cast<unsigned char>(c)); })
|
||||
| ranges::to<std::vector<std::byte>>;
|
||||
}
|
||||
template<> std::string read<char>(std::istream& input)
|
||||
{ return std::string{std::istreambuf_iterator<char>{input}, {}}; }
|
||||
}
|
||||
}
|
||||
|
||||
@@ -14,4 +14,13 @@ int main()
|
||||
auto serialized_bob = biu::serialize(bob);
|
||||
auto bob2 = biu::deserialize<student>(serialized_bob);
|
||||
assert(bob == bob2);
|
||||
struct A
|
||||
{
|
||||
int x;
|
||||
std::string y;
|
||||
};
|
||||
A a{ 123, "abc" };
|
||||
auto b = biu::deserialize<A>(biu::serialize(a));
|
||||
assert(a.x == b.x);
|
||||
assert(a.y == b.y);
|
||||
}
|
||||
|
||||
38
packages/biu/test/yaml.cpp
Normal file
38
packages/biu/test/yaml.cpp
Normal file
@@ -0,0 +1,38 @@
|
||||
# include <biu.hpp>
|
||||
# include <yaml-cpp/yaml.h>
|
||||
|
||||
int main()
|
||||
{
|
||||
using namespace biu::literals;
|
||||
std::string data = R"(
|
||||
a: [ 1, 2, 3 ]
|
||||
b: [ [ 1, 2, 3 ], [ 4, 5, 6 ], [ 7, 8, 9 ] ]
|
||||
c: [ 1, 2 ]
|
||||
)";
|
||||
auto node = YAML::Load(data);
|
||||
auto a = node["a"].as<Eigen::Vector3i>();
|
||||
auto a2 = node["a"].as<Eigen::VectorXi>();
|
||||
auto b = node["b"].as<Eigen::Matrix3i>();
|
||||
auto b2 = node["b"].as<Eigen::MatrixXi>();
|
||||
assert(a == a2);
|
||||
assert(a(0) == 1);
|
||||
assert(a(1) == 2);
|
||||
assert(a(2) == 3);
|
||||
assert(b == b2);
|
||||
auto c = node["c"].as<std::complex<double>>();
|
||||
assert(c == 1. + 2i);
|
||||
auto d = node["d"].as<std::optional<int>>();
|
||||
auto c3 = node["c"].as<std::optional<std::complex<double>>>();
|
||||
assert(d == std::nullopt);
|
||||
assert(c3 == 1. + 2i);
|
||||
struct A
|
||||
{
|
||||
Eigen::Vector3i a;
|
||||
Eigen::Matrix3i b;
|
||||
std::complex<double> c;
|
||||
};
|
||||
auto a3 = node.as<A>();
|
||||
assert(a3.a == a);
|
||||
assert(a3.b == b);
|
||||
assert(a3.c == c);
|
||||
}
|
||||
@@ -13,6 +13,7 @@ inputs: rec
|
||||
{
|
||||
inherit glad nodesoup;
|
||||
src = inputs.topInputs.matplotplusplus;
|
||||
stdenv = inputs.pkgs.clang18Stdenv;
|
||||
};
|
||||
zpp-bits = inputs.pkgs.callPackage ./zpp-bits.nix { src = inputs.topInputs.zpp-bits; };
|
||||
eigen = inputs.pkgs.callPackage ./eigen.nix { src = inputs.topInputs.eigen; };
|
||||
@@ -66,6 +67,7 @@ inputs: rec
|
||||
{
|
||||
inherit nameof zpp-bits tgbot-cpp concurrencpp pocketfft;
|
||||
stdenv = inputs.pkgs.clang18Stdenv;
|
||||
boost = inputs.pkgs.boost186;
|
||||
fmt = inputs.pkgs.fmt_11.overrideAttrs (prev: { patches = prev.patches or [] ++ [ ./biu/fmt.patch ]; });
|
||||
};
|
||||
zxorm = inputs.pkgs.callPackage ./zxorm.nix { src = inputs.topInputs.zxorm; };
|
||||
@@ -77,7 +79,7 @@ inputs: rec
|
||||
sbatch-tui = inputs.pkgs.callPackage ./sbatch-tui { inherit biu; stdenv = inputs.pkgs.clang18Stdenv; };
|
||||
ufo = inputs.pkgs.callPackage ./ufo
|
||||
{
|
||||
inherit concurrencpp biu matplotplusplus zpp-bits;
|
||||
inherit biu matplotplusplus;
|
||||
tbb = inputs.pkgs.tbb_2021_11;
|
||||
stdenv = inputs.pkgs.clang18Stdenv;
|
||||
};
|
||||
|
||||
@@ -7,7 +7,7 @@
|
||||
inherit src;
|
||||
cmakeFlags =
|
||||
[
|
||||
"-DBUILD_SHARED_LIBS=ON" "-DMATPLOTPP_BUILD_SHARED_LIBS=ON" "-DMATPLOTPP_BUILD_EXAMPLES=OFF"
|
||||
"-DMATPLOTPP_BUILD_EXAMPLES=OFF"
|
||||
"-DMATPLOTPP_WITH_SYSTEM_NODESOUP=ON" "-DMATPLOTPP_WITH_SYSTEM_CIMG=ON"
|
||||
"-DMATPLOTPP_BUILD_EXPERIMENTAL_OPENGL_BACKEND=ON" "-DGLAD_REPRODUCIBLE=ON"
|
||||
];
|
||||
|
||||
@@ -9,16 +9,16 @@ if(NOT CMAKE_BUILD_TYPE AND NOT CMAKE_CONFIGURATION_TYPES)
|
||||
set_property(CACHE CMAKE_BUILD_TYPE PROPERTY STRINGS "Debug" "Release" "MinSizeRel" "RelWithDebInfo")
|
||||
endif()
|
||||
|
||||
find_package(yaml-cpp REQUIRED)
|
||||
find_package(TBB REQUIRED)
|
||||
find_package(Matplot++ 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)
|
||||
add_executable(ufo src/fold.cpp src/unfold.cpp src/plot.cpp src/main.cpp)
|
||||
target_include_directories(ufo PRIVATE ${PROJECT_SOURCE_DIR}/include)
|
||||
target_link_libraries(ufo PRIVATE yaml-cpp TBB::tbb Matplot++::matplot biu::biu)
|
||||
target_link_libraries(ufo PRIVATE TBB::tbb Matplot++::matplot Matplot++::matplot_opengl 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})
|
||||
|
||||
|
||||
@@ -1,11 +1,11 @@
|
||||
{
|
||||
stdenv, cmake, pkg-config, version ? null,
|
||||
yaml-cpp, eigen, fmt, concurrencpp, highfive, tbb, matplotplusplus, biu, zpp-bits
|
||||
tbb, matplotplusplus, biu
|
||||
}: stdenv.mkDerivation
|
||||
{
|
||||
name = "ufo";
|
||||
src = ./.;
|
||||
buildInputs = [ yaml-cpp eigen fmt concurrencpp highfive tbb matplotplusplus biu zpp-bits ];
|
||||
buildInputs = [ tbb matplotplusplus biu ];
|
||||
nativeBuildInputs = [ cmake pkg-config ];
|
||||
doCheck = true;
|
||||
}
|
||||
|
||||
44
packages/ufo/include/ufo.hpp
Normal file
44
packages/ufo/include/ufo.hpp
Normal file
@@ -0,0 +1,44 @@
|
||||
# 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;
|
||||
|
||||
void fold(std::string config_file);
|
||||
void unfold(std::string config_file);
|
||||
void plot(std::string config_file);
|
||||
|
||||
// unfold 和 plot 都需要用到这个，所以写出来
|
||||
struct UnfoldOutput
|
||||
{
|
||||
Eigen::Matrix3d PrimativeCell;
|
||||
|
||||
// 关于各个 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;
|
||||
};
|
||||
}
|
||||
@@ -1,36 +0,0 @@
|
||||
# pragma once
|
||||
# include <ufo/solver.hpp>
|
||||
|
||||
namespace ufo
|
||||
{
|
||||
class FoldSolver : public Solver
|
||||
{
|
||||
public:
|
||||
struct InputType
|
||||
{
|
||||
Eigen::Matrix<int, 3, 3> SuperCellTransformation;
|
||||
std::vector<Eigen::Vector3d> Qpoints;
|
||||
DataFile OutputFile;
|
||||
|
||||
InputType(std::string config_file);
|
||||
};
|
||||
struct OutputType
|
||||
{
|
||||
std::vector<Eigen::Vector3d> Qpoints;
|
||||
using serialize = zpp::bits::members<1>;
|
||||
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::Matrix<int, 3, 3> super_cell_transformation
|
||||
);
|
||||
};
|
||||
}
|
||||
@@ -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
|
||||
);
|
||||
};
|
||||
}
|
||||
@@ -1,61 +0,0 @@
|
||||
# pragma once
|
||||
# include <any>
|
||||
# include <yaml-cpp/yaml.h>
|
||||
# include <matplot/matplot.h>
|
||||
# include <biu.hpp>
|
||||
|
||||
// 在相位中, 约定为使用 $\exp (2 \pi i \vec{q} \cdot \vec{r})$ 来表示原子的运动状态
|
||||
// (而不是 $\exp (-2 \pi i \vec{q} \cdot \vec{r})$)
|
||||
// 一些书定义的倒格矢中包含了 $2 \pi$ 的部分, 我们这里约定不包含这部分.
|
||||
// 也就是说, 正格子与倒格子的转置相乘, 得到单位矩阵.
|
||||
|
||||
namespace ufo
|
||||
{
|
||||
using namespace biu::literals;
|
||||
|
||||
class Solver
|
||||
{
|
||||
public:
|
||||
virtual Solver& operator()() = 0;
|
||||
virtual ~Solver() = default;
|
||||
|
||||
inline static void zpp_write(const auto& object, std::string filename)
|
||||
{
|
||||
auto data = biu::serialize(object);
|
||||
std::ofstream file(filename, std::ios::binary | std::ios::out);
|
||||
file.exceptions(std::ios::badbit | std::ios::failbit);
|
||||
static_assert(sizeof(std::byte) == sizeof(char));
|
||||
file.write(reinterpret_cast<const char*>(data.data()), data.size());
|
||||
}
|
||||
template <typename 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())
|
||||
);
|
||||
}
|
||||
return biu::deserialize<T>(data);
|
||||
}
|
||||
|
||||
struct DataFile
|
||||
{
|
||||
std::string Filename;
|
||||
std::string Format;
|
||||
std::map<std::string, std::any> ExtraParameters;
|
||||
DataFile() = default;
|
||||
DataFile
|
||||
(
|
||||
YAML::Node node, std::set<std::string> supported_format,
|
||||
std::string config_file, bool allow_same_as_config_file = false
|
||||
);
|
||||
};
|
||||
|
||||
};
|
||||
}
|
||||
@@ -1,169 +0,0 @@
|
||||
# pragma once
|
||||
# include <ufo/solver.hpp>
|
||||
|
||||
namespace ufo
|
||||
{
|
||||
// 反折叠的原理: 将超胞中的原子运动状态, 投影到一组平面波构成的基矢中.
|
||||
// 每一个平面波的波矢由两部分相加得到: 一部分是单胞倒格子的整数倍, 所取的个数有一定任意性, 论文中建议取大约单胞中原子个数那么多个;
|
||||
// 对于没有缺陷的情况, 取一个应该就足够了.
|
||||
// 这些平面波以原胞为周期。
|
||||
// 另一部分是超胞倒格子的整数倍, 取 n 个, n 为超胞对应的单胞的倍数, 其实也就是倒空间中单胞对应倒格子中超胞的格点.
|
||||
// 只要第一部分取得足够多, 那么单胞中原子的状态就可以完全被这些平面波描述.
|
||||
// 将超胞中原子的运动状态投影到这些基矢上, 计算出投影的系数, 就可以将超胞的原子运动状态分解到单胞中的多个 q 点上.
|
||||
class UnfoldSolver : public Solver
|
||||
{
|
||||
public:
|
||||
struct InputType
|
||||
{
|
||||
// 单胞的三个格矢，每行表示一个格矢的坐标，单位为埃
|
||||
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::Matrix3i SuperCellTransformation;
|
||||
Eigen::Matrix3d SuperCellDeformation;
|
||||
Eigen::Vector3i SuperCellMultiplier;
|
||||
// 在单胞内取几个平面波的基矢
|
||||
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;
|
||||
|
||||
// 将 Transformation 矩阵分解为一个正交矩阵乘以一个整数对角矩阵
|
||||
static std::pair<Eigen::Matrix3d, Eigen::Vector3i> decompose_transformation(Eigen::Matrix3i transformation);
|
||||
|
||||
// 构建基
|
||||
// 每个 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
|
||||
);
|
||||
};
|
||||
}
|
||||
@@ -1,45 +1,22 @@
|
||||
# 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);
|
||||
SuperCellTransformation = input["SuperCellTransformation"].as<std::array<std::array<int, 3>, 3>>() | biu::toEigen<>;
|
||||
for (auto& qpoint : input["Qpoints"].as<std::vector<std::array<double, 3>>>())
|
||||
Qpoints.push_back(qpoint | biu::toEigen<>);
|
||||
OutputFile = DataFile(input["OutputFile"], {"yaml"}, config_file);
|
||||
}
|
||||
void FoldSolver::OutputType::write(std::string filename) const
|
||||
Eigen::Matrix<int, 3, 3> SuperCellTransformation;
|
||||
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 << (" - [ {:.8f}, {:.8f}, {:.8f} ]\n"_f(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_.Transformation));
|
||||
}
|
||||
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::Matrix<int, 3, 3> transformation
|
||||
)
|
||||
Eigen::Matrix<int, 3, 3> super_cell_transformation
|
||||
) -> Eigen::Vector3d
|
||||
{
|
||||
/*
|
||||
首先需要将 q 点坐标的单位转换为 ModifiedSuperCell 的格矢，可知：
|
||||
@@ -58,7 +35,7 @@ namespace ufo
|
||||
*/
|
||||
auto qpoint_by_reciprocal_super_cell =
|
||||
(
|
||||
transformation.cast<double>()
|
||||
super_cell_transformation.cast<double>()
|
||||
* qpoint_in_reciprocal_primitive_cell_by_reciprocal_primitive_cell
|
||||
).eval();
|
||||
/*
|
||||
@@ -66,5 +43,16 @@ namespace ufo
|
||||
这等价于直接取 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.SuperCellTransformation);
|
||||
})
|
||||
| ranges::to_vector;
|
||||
|
||||
// 默认的输出太丑了，但是不想手动写了，忍一下
|
||||
std::ofstream(input.OutputFile.value_or("output.yaml")) << YAML::Node(output);
|
||||
}
|
||||
|
||||
@@ -1,17 +1,15 @@
|
||||
# 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]}();
|
||||
ufo::fold(argv[2]);
|
||||
else if (argv[1] == std::string("unfold"))
|
||||
ufo::UnfoldSolver{argv[2]}();
|
||||
else if (argv[1] == std::string("plot"))
|
||||
ufo::PlotSolver{argv[2]}();
|
||||
ufo::unfold(argv[2]);
|
||||
// else if (argv[1] == std::string("plot"))
|
||||
// ufo::PlotSolver{argv[2]}();
|
||||
else
|
||||
throw std::runtime_error(fmt::format("Unknown task: {}", argv[1]));
|
||||
}
|
||||
|
||||
@@ -1,227 +1,166 @@
|
||||
# include <ufo/plot.hpp>
|
||||
# include <ufo.hpp>
|
||||
# include <matplot/matplot.h>
|
||||
# include <matplot/backend/opengl.h>
|
||||
# include <boost/container/flat_map.hpp>
|
||||
|
||||
namespace ufo
|
||||
void ufo::plot(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;
|
||||
struct { std::size_t X, Y; } Resolution;
|
||||
// 画图的频率范围
|
||||
struct { double Min, Max; } FrequencyRange;
|
||||
// 搜索 q 点时的阈值，单位为埃^-1
|
||||
std::optional<double> ThresholdWhenSearchingQpoints;
|
||||
// 是否要在 z 轴上作一些标记
|
||||
std::optional<std::vector<double>> 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);
|
||||
PrimativeCell = input["PrimativeCell"].as<std::array<std::array<double, 3>, 3>>() | biu::toEigen<>;
|
||||
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 };
|
||||
biu::Hdf5file(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::pair<unsigned, unsigned>& resolution,
|
||||
const std::pair<double, double>& 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,
|
||||
unsigned resolution, std::pair<double, double> 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 (unsigned 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 (unsigned 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 (unsigned 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++)
|
||||
{
|
||||
int index = (frequency[i] - range.first) / (range.second - range.first) * resolution;
|
||||
if (index >= 0 && index < static_cast<int>(resolution))
|
||||
result[index] += weight[i];
|
||||
int index = (frequency[i] - frequency_range.first) / (frequency_range.second - frequency_range.first)
|
||||
* resolution;
|
||||
if (index >= 0 && index < static_cast<int>(resolution)) result[index] += weight[i];
|
||||
}
|
||||
return result;
|
||||
};
|
||||
|
||||
std::vector<std::vector<double>> values;
|
||||
for (unsigned i = 0; i < resolution.first; 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 it = path.lower_bound(current_distance);
|
||||
if (it == path.begin()) values.push_back(blend
|
||||
(it->second, it->second, true, 1, resolution.second, frequency_range));
|
||||
else if (it == path.end()) values.push_back(blend
|
||||
(
|
||||
std::prev(it)->second, std::prev(it)->second, true, 1,
|
||||
resolution.second, 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.second, frequency_range
|
||||
));
|
||||
}
|
||||
return {values, x_ticks};
|
||||
}
|
||||
void PlotSolver::plot
|
||||
return values;
|
||||
};
|
||||
|
||||
// 根据数值, 画图
|
||||
auto plot = []
|
||||
(
|
||||
const std::vector<std::vector<double>>& values,
|
||||
const std::string& filename,
|
||||
@@ -237,20 +176,16 @@ namespace ufo
|
||||
for (unsigned 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<matplot::backend::opengl>(true);
|
||||
auto ax = f->current_axes();
|
||||
auto image = ax->image(std::tie(r, g, b));
|
||||
image->matrix_a(a);
|
||||
@@ -260,5 +195,67 @@ namespace ufo
|
||||
ax->y_axis().tick_values(y_ticks);
|
||||
ax->y_axis().tick_length(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.Resolution.X, input.Resolution.Y},
|
||||
{input.FrequencyRange.Min, input.FrequencyRange.Max}, total_distance
|
||||
);
|
||||
auto x_ticks = x_ticks_index | ranges::views::transform([&](auto i)
|
||||
{ return path.nth(i)->first / total_distance * input.Resolution.X; }) | ranges::to<std::vector>;
|
||||
std::vector<double> y_ticks;
|
||||
if (input.YTicks) y_ticks = input.YTicks.value()
|
||||
| ranges::views::transform([&](auto i)
|
||||
{ return (i - input.FrequencyRange.Min) / (input.FrequencyRange.Max - input.FrequencyRange.Min)
|
||||
* input.Resolution.Y; })
|
||||
| ranges::to<std::vector>;
|
||||
if (input.OutputPictureFile) plot(values, input.OutputPictureFile.value(), x_ticks, y_ticks);
|
||||
if (input.OutputDataFile)
|
||||
biu::Hdf5file(input.OutputDataFile.value(), true)
|
||||
.write("Values", values)
|
||||
.write("XTicks", x_ticks)
|
||||
.write("YTicks", y_ticks)
|
||||
.write("Resolution", std::vector{input.Resolution.X, input.Resolution.Y})
|
||||
.write("Range", std::vector{input.FrequencyRange.Min, input.FrequencyRange.Max});
|
||||
}
|
||||
|
||||
@@ -1,34 +0,0 @@
|
||||
# include <ufo/solver.hpp>
|
||||
|
||||
namespace ufo
|
||||
{
|
||||
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("Unsupported format: \"{}\""_f(Format));
|
||||
if (auto _ = node["RelativeToConfigFile"])
|
||||
{
|
||||
auto __ = _.as<bool>();
|
||||
ExtraParameters["RelativeToConfigFile"] = __;
|
||||
if (__)
|
||||
Filename = std::filesystem::path(config_file).parent_path() / Filename;
|
||||
}
|
||||
};
|
||||
|
||||
}
|
||||
@@ -1,289 +1,121 @@
|
||||
# include <ufo.hpp>
|
||||
# include <thread>
|
||||
# include <syncstream>
|
||||
# include <execution>
|
||||
# include <ufo/unfold.hpp>
|
||||
|
||||
namespace ufo
|
||||
void ufo::unfold(std::string config_file)
|
||||
{
|
||||
UnfoldSolver::InputType::InputType(std::string filename)
|
||||
{
|
||||
// read main input file
|
||||
{
|
||||
auto node = YAML::LoadFile(filename);
|
||||
PrimativeCell = node["PrimativeCell"].as<std::array<std::array<double, 3>, 3>>() | biu::toEigen<>;
|
||||
SuperCellTransformation =
|
||||
node["SuperCellTransformation"].as<std::array<std::array<int, 3>, 3>>() | biu::toEigen<>;
|
||||
std::tie(SuperCellDeformation, SuperCellMultiplier) =
|
||||
decompose_transformation(SuperCellTransformation);
|
||||
PrimativeCellBasisNumber = node["PrimativeCellBasisNumber"].as<std::array<unsigned, 3>>() | biu::toEigen<>;
|
||||
// 反折叠的原理: 将超胞中的原子运动状态, 投影到一组平面波构成的基矢中.
|
||||
// 每一个平面波的波矢由两部分相加得到: 一部分是单胞倒格子的整数倍, 所取的个数有一定任意性, 论文中建议取大约单胞中原子个数那么多个;
|
||||
// 对于没有缺陷的情况, 取一个应该就足够了.
|
||||
// 这些平面波以原胞为周期。
|
||||
// 另一部分是超胞倒格子的整数倍, 取 n 个, n 为超胞对应的单胞的倍数, 其实也就是倒空间中单胞对应倒格子中超胞的格点.
|
||||
// 只要第一部分取得足够多, 那么单胞中原子的状态就可以完全被这些平面波描述.
|
||||
// 将超胞中原子的运动状态投影到这些基矢上, 计算出投影的系数, 就可以将超胞的原子运动状态分解到单胞中的多个 q 点上.
|
||||
|
||||
AtomPositionInputFile = DataFile
|
||||
(
|
||||
node["AtomPositionInputFile"], {"yaml"},
|
||||
filename, true
|
||||
);
|
||||
QpointDataInputFile = DataFile
|
||||
(
|
||||
node["QpointDataInputFile"], {"yaml", "hdf5"},
|
||||
filename, true
|
||||
);
|
||||
if (auto value = node["QpointDataOutputFile"])
|
||||
struct Input
|
||||
{
|
||||
// 单胞的三个格矢，每行表示一个格矢的坐标，单位为埃
|
||||
Eigen::Matrix3d PrimativeCell;
|
||||
|
||||
// 单胞到超胞的格矢转换时用到的矩阵
|
||||
// SuperCellMultiplier 是一个三维列向量且各个元素都是整数，表示单胞在各个方向扩大到多少倍之后，可以得到和超胞一样的体积
|
||||
// SuperCellDeformation 是一个行列式为 1 的矩阵，它表示经过 SuperCellMultiplier 扩大后，还需要怎样的变换才能得到超胞
|
||||
// SuperCell = (SuperCellDeformation * SuperCellMultiplier.asDiagonal()) * PrimativeCell
|
||||
// ReciprocalPrimativeCell = (SuperCellDeformation * SuperCellMultiplier.asDiagonal()).transpose()
|
||||
// * ReciprocalSuperCell
|
||||
// Position = PositionToCell(line vector) * Cell
|
||||
// InversePosition = InversePositionToCell(line vector) * ReciprocalCell
|
||||
// PositionToSuperCell(line vector) * SuperCell = PositionToPrimativeCell(line vector) * PrimativeCell
|
||||
// ReciprocalPositionToSuperCell(line vector) * ReciprocalSuperCell
|
||||
// = ReciprocalPositionToPrimativeCell(line vector) * ReciprocalPrimativeCell
|
||||
Eigen::Matrix3i SuperCellTransformation;
|
||||
|
||||
// 在单胞内取几个平面波的基矢
|
||||
Eigen::Vector<unsigned, 3> PrimativeCellBasisNumber;
|
||||
|
||||
// 超胞中原子的坐标，每行表示一个原子的坐标，单位为超胞的格矢
|
||||
Eigen::MatrixX3d AtomPositionBySuperCell;
|
||||
|
||||
// 从 band.hdf5 读入 QpointData
|
||||
std::optional<std::string> QpointDataInputFile;
|
||||
|
||||
// 输出到哪些文件
|
||||
struct QpointDataOutputFileType
|
||||
{
|
||||
std::string Filename;
|
||||
|
||||
// 如果指定，则将结果投影到那些原子上
|
||||
std::optional<std::vector<std::size_t>> SelectedAtoms;
|
||||
|
||||
// 默认输出为 zpp 文件，如果指定为 true，则输出为 yaml 文件
|
||||
std::optional<bool> OutputAsYaml;
|
||||
};
|
||||
std::vector<QpointDataOutputFileType> QpointDataOutputFile;
|
||||
};
|
||||
|
||||
// 关于各个 Q 点的数据
|
||||
struct QpointData
|
||||
{
|
||||
// Q 点的坐标，单位为超胞的倒格矢
|
||||
Eigen::Vector3d Qpoint;
|
||||
|
||||
// 关于这个 Q 点上各个模式的数据
|
||||
struct ModeDataType
|
||||
{
|
||||
// 模式的频率，单位为 THz
|
||||
double Frequency;
|
||||
// 模式中各个原子的运动状态
|
||||
// 这个数据应当是这样得到的：动态矩阵的 eigenvector 乘以 $\exp(-2 \pi i \vec q \cdot \vec r)$
|
||||
// 这个数据可以认为是原子位移中, 关于超胞有周期性的那一部分, 再乘以原子质量的开方.
|
||||
// 这个数据在读入后会被立即归一化.
|
||||
Eigen::MatrixX3cd AtomMovement;
|
||||
};
|
||||
std::vector<ModeDataType> ModeData;
|
||||
};
|
||||
|
||||
// 从文件中读取 QpointData
|
||||
auto read_qpoint_data = [](std::string filename)
|
||||
{
|
||||
// 读入原始数据
|
||||
std::vector<std::vector<std::vector<double>>> frequency, path;
|
||||
std::vector<std::vector<std::vector<std::vector<biu::PhonopyComplex>>>> eigenvector_vector;
|
||||
biu::Hdf5file(filename).read("/frequency", frequency)
|
||||
.read("/eigenvector", eigenvector_vector)
|
||||
.read("/path", path);
|
||||
|
||||
// 整理得到结果
|
||||
std::vector size = { frequency.size(), frequency[0].size(), frequency[0][0].size() };
|
||||
std::vector<QpointData> qpoint(size[0] * size[1]);
|
||||
for (unsigned i = 0; i < size[0]; i++) for (unsigned j = 0; j < size[1]; j++)
|
||||
{
|
||||
qpoint[i * size[1] + j].Qpoint = Eigen::Vector3d(path[i][j].data());
|
||||
qpoint[i * size[1] + j].ModeData.resize(size[2]);
|
||||
for (unsigned k = 0; k < size[2]; k++)
|
||||
{
|
||||
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[i * size[1] + j].ModeData[k].Frequency = frequency[i][j][k];
|
||||
auto n_modes = eigenvector_vector[i][j].size() / 3;
|
||||
Eigen::MatrixX3cd eigenvectors(n_modes, 3);
|
||||
for (unsigned l = 0; l < n_modes; 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;
|
||||
// 原则上讲，需要对读入的原子运动状态作相位转换, 使得它们与我们的约定一致(对超胞周期性重复)，但这个转换 phonopy 已经做了
|
||||
// 这里还要需要做归一化处理 (指将数据简单地作为向量处理的归一化)
|
||||
qpoint[i * size[1] + j].ModeData[k].AtomMovement = eigenvectors / eigenvectors.norm();
|
||||
}
|
||||
}
|
||||
return qpoint;
|
||||
};
|
||||
|
||||
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++)
|
||||
atom_position_to_super_cell.row(i) = points[i]["coordinates"].as<std::array<double, 3>>() | biu::toEigen<>;
|
||||
auto super_cell = (SuperCellTransformation.cast<double>() * 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++)
|
||||
{
|
||||
QpointData[i].Qpoint = phonon[i]["q-position"].as<std::array<double, 3>>() | biu::toEigen<>;
|
||||
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<biu::PhonopyComplex>>>> eigenvector_vector;
|
||||
biu::Hdf5file(QpointDataInputFile.Filename, true).read("/frequency", frequency)
|
||||
.read("/eigenvector", eigenvector_vector)
|
||||
.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 << " - Qpoint: [ {1:.{0}f}, {2:.{0}f}, {3:.{0}f} ]\n"_f
|
||||
(percision, qpoint.Qpoint[0], qpoint.Qpoint[1], qpoint.Qpoint[2]);
|
||||
print << " Source: [ {1:.{0}f}, {2:.{0}f}, {3:.{0}f} ]\n"_f
|
||||
(percision, qpoint.Source[0], qpoint.Source[1], qpoint.Source[2]);
|
||||
print << " ModeData:\n";
|
||||
for (auto& mode: qpoint.ModeData)
|
||||
print << " - {{ Frequency: {1:.{0}f}, Weight: {2:.{0}f} }}\n"_f
|
||||
(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);
|
||||
}
|
||||
}
|
||||
biu::Hdf5file(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;
|
||||
}
|
||||
|
||||
std::pair<Eigen::Matrix3d, Eigen::Vector3i> UnfoldSolver::decompose_transformation(Eigen::Matrix3i transformation)
|
||||
// 将 SuperCellTransformation 矩阵分解为 SuperCellDeformation 和 SuperCellMultiplier
|
||||
auto decompose_transformation = [](Eigen::Matrix3i transformation)
|
||||
-> std::pair<Eigen::Matrix3d, Eigen::Vector3i>
|
||||
{
|
||||
// 三种整数基础矩阵
|
||||
struct Multiply { unsigned at; int multiply; };
|
||||
struct Add { unsigned from, to; int multiply; };
|
||||
struct Exchange { unsigned from, to; };
|
||||
|
||||
// 将 SuperCellTransformation 分解为一系列基础矩阵的乘积，并将这些基础矩阵从左到右输出
|
||||
auto decompose = [](Eigen::Matrix3i matrix)
|
||||
-> concurrencpp::generator<std::variant<Add, Exchange, Multiply>>
|
||||
{
|
||||
@@ -390,7 +222,7 @@ namespace ufo
|
||||
{
|
||||
auto [from, to, multiply] = std::get<Add>(i);
|
||||
auto transform = Eigen::Matrix3d::Identity().eval();
|
||||
transform(to, from) = multiply / multiplier(from) * multiplier(to);
|
||||
transform(to, from) = static_cast<double>(multiply) / multiplier(from) * multiplier(to);
|
||||
deformatin = deformatin * transform;
|
||||
}
|
||||
else if (std::holds_alternative<Exchange>(i))
|
||||
@@ -401,48 +233,65 @@ namespace ufo
|
||||
}
|
||||
}
|
||||
return {deformatin, multiplier};
|
||||
}
|
||||
};
|
||||
|
||||
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<unsigned, 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]
|
||||
std::vector<std::vector<Eigen::VectorXcd>> basis(super_cell_multiplier.prod());
|
||||
// diff_of_sub_qpoint 表示 sub qpoint 与 qpoint 的相对位置，单位为超胞的倒格矢
|
||||
for (auto [diff_of_sub_qpoint_by_reciprocal_modified_super_cell, i_of_sub_qpoint]
|
||||
: biu::sequence(super_cell_multiplier))
|
||||
{
|
||||
basis[i_of_sub_qpoint].resize(primative_cell_basis_number.prod());
|
||||
for (auto [xyz_of_basis, i_of_basis] : biu::sequence(primative_cell_basis_number))
|
||||
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>();
|
||||
// DiffOfSubQpoint
|
||||
// = (DiffOfSubQpointByReciprocalPrimativeCell.transpose() * ReciprocalPrimativeCell).transpose()
|
||||
// ReciprocalPrimativeCell = PrimativeCell.transpose().inverse()
|
||||
// SuperCell = SuperCellTransformation * PrimativeCell
|
||||
// ReciprocalSuperCell = SuperCell.transpose().inverse()
|
||||
// AtomPosition = AtomPositionBySuperCell * SuperCell.transpose()
|
||||
// 整理得到：
|
||||
//
|
||||
// 将单位转换为埃^-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,
|
||||
const std::vector<std::vector<Eigen::VectorXcd>>& basis,
|
||||
// 实际上只需要其中的 AtomMovement
|
||||
const std::vector<QpointData>& qpoint_data,
|
||||
std::atomic<unsigned>& 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());
|
||||
// 对每个模式并行
|
||||
@@ -463,36 +312,61 @@ namespace ufo
|
||||
// 但这里并不是这样一个严格的概念. 因此对分解到各个 sub qpoint 上的权重做归一化处理
|
||||
auto sum = std::accumulate
|
||||
(projection_coefficient.begin(), projection_coefficient.end(), 0.);
|
||||
for (auto& _ : projection_coefficient)
|
||||
_ /= sum;
|
||||
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++
|
||||
i_of_meta_qpoint < qpoint_data.size(); i_of_meta_qpoint++
|
||||
)
|
||||
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 Eigen::Matrix3d& primative_cell,
|
||||
const Eigen::Vector3i& super_cell_multiplier,
|
||||
const Eigen::Matrix3d& super_cell_deformation,
|
||||
const std::vector<std::vector<std::vector<double>>>& projection_coefficient,
|
||||
const std::vector<QpointData>& qpoint_data,
|
||||
const std::optional<std::vector<std::size_t>>& selected_atoms
|
||||
)
|
||||
{
|
||||
UnfoldOutput output;
|
||||
output.PrimativeCell = primative_cell;
|
||||
for (unsigned 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]
|
||||
: biu::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 = 0; i < qpoint_data[i_of_meta_qpoint].ModeData.size(); i++)
|
||||
{
|
||||
projection_coefficient_on_atoms.value().emplace_back(0);
|
||||
for (auto atom : *selected_atoms)
|
||||
projection_coefficient_on_atoms.value().back()
|
||||
+= qpoint_data[i_of_meta_qpoint].ModeData[i].AtomMovement.row(atom).array().abs2().sum();
|
||||
}
|
||||
}
|
||||
|
||||
for
|
||||
(
|
||||
auto [diff_of_sub_qpoint_by_reciprocal_modified_super_cell, i_of_sub_qpoint]
|
||||
: biu::sequence(super_cell_multiplier)
|
||||
)
|
||||
{
|
||||
auto& _ = output.QpointData.emplace_back();
|
||||
/*
|
||||
@@ -520,24 +394,84 @@ namespace ufo
|
||||
(
|
||||
super_cell_multiplier.cast<double>().cwiseInverse().asDiagonal()
|
||||
* (
|
||||
xyz_of_diff_of_sub_qpoint_by_reciprocal_modified_super_cell.cast<double>()
|
||||
+ super_cell_deformation.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>()
|
||||
+ super_cell_deformation.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 (unsigned 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;
|
||||
};
|
||||
|
||||
std::clog << "Reading input file... " << std::flush;
|
||||
auto input = YAML::LoadFile(config_file).as<Input>();
|
||||
auto qpoint_data
|
||||
= read_qpoint_data(input.QpointDataInputFile.value_or("band.hdf5"));
|
||||
std::clog << "Done." << std::endl;
|
||||
|
||||
std::clog << "Constructing basis... " << std::flush;
|
||||
auto [super_cell_deformation, super_cell_multiplier]
|
||||
= decompose_transformation(input.SuperCellTransformation);
|
||||
auto basis = construct_basis
|
||||
(
|
||||
input.PrimativeCell, super_cell_multiplier,
|
||||
input.PrimativeCellBasisNumber,
|
||||
input.AtomPositionBySuperCell
|
||||
* (input.SuperCellTransformation.cast<double>() * input.PrimativeCell).transpose()
|
||||
);
|
||||
std::clog << "Done." << std::endl;
|
||||
|
||||
std::clog << "Calculating projection coefficient... " << std::flush;
|
||||
// 用来在屏幕上输出进度的计数器和线程
|
||||
std::atomic<unsigned> 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.PrimativeCell, super_cell_multiplier, super_cell_deformation,
|
||||
projection_coefficient, qpoint_data, output_file.SelectedAtoms
|
||||
);
|
||||
if (output_file.OutputAsYaml.value_or(false))
|
||||
std::ofstream(output_file.Filename, std::ios::binary) << biu::serialize<char>(output);
|
||||
else
|
||||
std::ofstream(output_file.Filename) << YAML::Node(output);
|
||||
}
|
||||
std::clog << "Done." << std::endl;
|
||||
}
|
||||
|
||||
@@ -1,4 +1,4 @@
|
||||
Transformation:
|
||||
SuperCellTransformation:
|
||||
- [ 3, 0, 0 ]
|
||||
- [ 2, 4, 0 ]
|
||||
- [ 0, 0, 1 ]
|
||||
@@ -14,6 +14,4 @@ Qpoints:
|
||||
- [0.4, 0, 0]
|
||||
- [0.45, 0, 0]
|
||||
- [0.5, 0, 0]
|
||||
OutputFile:
|
||||
Filename: fold-output.yaml
|
||||
Format: yaml
|
||||
OutputFile: fold-output.yaml
|
||||
|
||||
@@ -1,12 +1,45 @@
|
||||
Qpoints:
|
||||
- [ 0.00000000, 0.00000000, 0.00000000 ]
|
||||
- [ 0.15000000, 0.10000000, 0.00000000 ]
|
||||
- [ 0.30000000, 0.20000000, 0.00000000 ]
|
||||
- [ 0.45000000, 0.30000000, 0.00000000 ]
|
||||
- [ 0.60000000, 0.40000000, 0.00000000 ]
|
||||
- [ 0.75000000, 0.50000000, 0.00000000 ]
|
||||
- [ 0.90000000, 0.60000000, 0.00000000 ]
|
||||
- [ 0.05000000, 0.70000000, 0.00000000 ]
|
||||
- [ 0.20000000, 0.80000000, 0.00000000 ]
|
||||
- [ 0.35000000, 0.90000000, 0.00000000 ]
|
||||
- [ 0.50000000, 0.00000000, 0.00000000 ]
|
||||
-
|
||||
- 0
|
||||
- 0
|
||||
- 0
|
||||
-
|
||||
- 0.15000000000000002
|
||||
- 0.10000000000000001
|
||||
- 0
|
||||
-
|
||||
- 0.30000000000000004
|
||||
- 0.20000000000000001
|
||||
- 0
|
||||
-
|
||||
- 0.44999999999999996
|
||||
- 0.29999999999999999
|
||||
- 0
|
||||
-
|
||||
- 0.60000000000000009
|
||||
- 0.40000000000000002
|
||||
- 0
|
||||
-
|
||||
- 0.75
|
||||
- 0.5
|
||||
- 0
|
||||
-
|
||||
- 0.89999999999999991
|
||||
- 0.59999999999999998
|
||||
- 0
|
||||
-
|
||||
- 0.049999999999999822
|
||||
- 0.69999999999999996
|
||||
- 0
|
||||
-
|
||||
- 0.20000000000000018
|
||||
- 0.80000000000000004
|
||||
- 0
|
||||
-
|
||||
- 0.35000000000000009
|
||||
- 0.90000000000000002
|
||||
- 0
|
||||
-
|
||||
- 0.5
|
||||
- 0
|
||||
- 0
|
||||
Reference in New Issue
Block a user