diff --git a/CMakeLists.txt b/CMakeLists.txt index 8fbf50c..2f05f91 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -14,14 +14,14 @@ find_package(Matplot++ REQUIRED) find_package(biu REQUIRED) find_package(Threads REQUIRED) -add_executable(ufo src/fold.cpp src/unfold.cpp src/plot.cpp src/raman-create-displacement.cpp src/raman-extract.cpp - src/raman-apply-contribution.cpp src/main.cpp) +add_executable(ufo src/fold.cpp src/unfold.cpp src/plot.cpp src/raman.cpp src/project-to-atom.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) if (CMAKE_BUILD_TYPE STREQUAL "Debug") target_compile_definitions(ufo PRIVATE BIU_LOGGER_DEBUG BIU_LOGGER_SOURCE_ROOT="${CMAKE_SOURCE_DIR}") + target_compile_options(ufo PRIVATE -ftemplate-backtrace-limit=0) endif() install(TARGETS ufo RUNTIME DESTINATION ${CMAKE_INSTALL_BINDIR}) diff --git a/include/ufo.hpp b/include/ufo.hpp index eaa4827..68424f6 100644 --- a/include/ufo.hpp +++ b/include/ufo.hpp @@ -11,82 +11,92 @@ namespace ufo using namespace biu::literals; using namespace biu::stream_operators; + // 许多函数通用的一个结构体 + struct CommonData + { + // 原胞相关信息，其中存储的信息是直接读入的，而不是反折叠后的结果 + struct + { + // 晶胞，单位为埃，每行代表一个格矢 + Eigen::Matrix3d Cell; + // 原子种类及其个数 + std::vector> AtomType; + // 原子位置，单位为原胞格矢 + Eigen::MatrixX3d AtomPosition; + // 单胞中每个 q 点的信息 + struct QpointType + { + // q 点的坐标，单位为倒格矢 + Eigen::Vector3d Qpoint; + struct ModeType + { + // 振动频率，单位为 THz + double Frequency; + // 振动模式，单位为埃 + Eigen::MatrixX3cd EigenVector; + // 拉曼张量 + // TODO: 单位是什么？ + std::optional WeightOnRaman; + // work around for fpr + std::optional, 3>> RamanTensor; + using serialize = zpp::bits::members<4>; + }; + // 各个模式的的信息 + std::vector Mode; + using serialize = zpp::bits::members<2>; + }; + std::vector Qpoint; + using serialize = zpp::bits::members<4>; + } Primative; + + // 超胞相关信息，包括反折叠后的信息 + struct + { + // 超胞的晶胞，单位为埃，每行代表一个格矢 + Eigen::Matrix3d Cell; + // 单胞与到超胞的变换 + Eigen::Vector3i CellMultiplier; + Eigen::Matrix3d CellDeformation; + std::vector> AtomType; + Eigen::MatrixX3d AtomPosition; + struct QpointType + { + Eigen::Vector3d Qpoint; + std::vector SubQpoint; + struct ModeType + { + double Frequency; + Eigen::MatrixX3cd EigenVector; + std::vector WeightOnUnfold; + std::optional WeightOnSelectedAtom; + std::optional WeightOnRaman; + std::optional, 3>> RamanTensor; + using serialize = zpp::bits::members<6>; + }; + std::vector Mode; + using serialize = zpp::bits::members<3>; + }; + std::vector Qpoint; + using serialize = zpp::bits::members<6>; + } Super; + + // 选择的原子 + std::optional> SelectedAtom; + + // 拉曼散射选择的偏振方向 + std::optional> RamanPolarization; + + // 各种原子的质量 + std::map AtomMass; + using serialize = zpp::bits::members<5>; + }; + void fold(std::string config_file); void unfold(std::string config_file); + void project_to_atom(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_extract(std::vector files); void raman_apply_contribution(std::string config_file); - - // 许多函数都需要用到这个，所以写到头文件中 - struct UnfoldOutput - { - Eigen::Matrix3d PrimativeCell; - Eigen::Matrix3i SuperCellTransformation; - Eigen::Vector3i SuperCellMultiplier; - Eigen::Matrix3d SuperCellDeformation; - std::optional> 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 ModeData; - }; - std::vector 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 ModeData; - }; - std::vector MetaQpointData; - - using serialize = zpp::bits::members<7>; - }; - struct DisplacementOutput - { - struct ModeData_t - { - std::size_t MetaQpointIndex; - std::size_t ModeIndex; - // 每个原子的位移，单位为埃 - Eigen::MatrixX3d AtomMovement; - // 为了使得位移最大的原子的位移恰好是 input.MaxDisplacement, 在位移上乘以了多大的系数 - double Ratio; - }; - std::vector ModeData; - Eigen::Vector3d AtomMasses; - double MaxDisplacement; - std::vector QpointIndices; - using serialize = zpp::bits::members<4>; - }; } diff --git a/src/fold.cpp b/src/fold.cpp index bc189bc..0283c89 100644 --- a/src/fold.cpp +++ b/src/fold.cpp @@ -2,17 +2,13 @@ void ufo::fold(std::string config_file) { - struct Input + struct Config { Eigen::Matrix3d SuperCellDeformation; Eigen::Vector3i SuperCellMultiplier; std::vector Qpoints; - std::optional OutputFile; - }; - struct Output - { - std::vector Qpoints; }; + auto fold = [] ( Eigen::Vector3d qpoint_in_reciprocal_primitive_cell_by_reciprocal_primitive_cell, @@ -44,15 +40,13 @@ void ufo::fold(std::string config_file) */ return (qpoint_by_reciprocal_super_cell.array() - qpoint_by_reciprocal_super_cell.array().floor()).matrix(); }; - auto input = YAML::LoadFile(config_file).as(); - Output output; - output.Qpoints = input.Qpoints - | ranges::views::transform([&](auto& qpoint) - { - return fold(qpoint, input.SuperCellDeformation * input.SuperCellMultiplier.cast().asDiagonal()); - }) - | ranges::to_vector; - - // 默认的输出太丑了，但是不想手动写了，忍一下 - std::ofstream(input.OutputFile.value_or("output.yaml")) << YAML::Node(output); + + biu::Logger::Guard log(config_file); + auto input = YAML::LoadFile(config_file).as(); + for (const auto& qpoint : input.Qpoints) log.info("{} -> {}"_f + ( + qpoint, + fold(qpoint, + input.SuperCellDeformation * input.SuperCellMultiplier.cast().asDiagonal()) + )); } diff --git a/src/main.cpp b/src/main.cpp index d082956..bb84c94 100644 --- a/src/main.cpp +++ b/src/main.cpp @@ -5,6 +5,9 @@ int main(int argc, const char** argv) using namespace biu::literals; if (argv[1] == "fold"s) ufo::fold(argv[2]); else if (argv[1] == "unfold"s) ufo::unfold(argv[2]); + else if (argv[1] == "zpp-to-yaml"s) std::cout << YAML::Node(biu::deserialize + (biu::read(std::cin))); + else if (argv[1] == "project-to-atom"s) ufo::project_to_atom(argv[2]); else if (argv[1] == "raman-create-displacement"s) ufo::raman_create_displacement(argv[2]); else if (argv[1] == "raman-extract"s) ufo::raman_extract ( diff --git a/src/plot.cpp b/src/plot.cpp index c756cdc..8f490d6 100644 --- a/src/plot.cpp +++ b/src/plot.cpp @@ -4,9 +4,9 @@ void ufo::plot_band(std::string config_file) { - struct Input + struct Config { - std::string UnfoldedDataFile; + std::string InputDataFile; // 要画图的 q 点路径列表 // 内层表示一个路径上的 q 点，外层表示不同的路径 // 单位为倒格矢 @@ -29,26 +29,26 @@ void ufo::plot_band(std::string config_file) std::optional OutputDataFile; }; - // 根据 q 点路径, 搜索要使用的 q 点，返回的是 q 点在 QpointData 中的索引以及到路径起点的距离，以及这段路径的总长度 + // 根据 q 点路径, 搜索要使用的 q 点，返回的是超胞中 Qpint 的索引和 SubQpoint 的索引，以及到路径起点的距离，以及这段路径的总长度 auto search_qpoints = [] ( const Eigen::Matrix3d& primative_cell, const std::pair& path, - const std::vector& qpoints, + const std::vector>& qpoints, double threshold, bool exclude_endpoint = false ) { // 对于 output 中的每一个点, 检查这个点是否在路径上. 如果在, 把它加入到 selected_qpoints 中 // 键为这个点到起点的距离 - boost::container::flat_map selected_qpoints; + boost::container::flat_map> 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 (std::size_t i = 0; i < qpoints.size(); i++) for (std::size_t j = 0; j < qpoints[i].size(); j++) for (auto cell_shift : biu::sequence(Eigen::Vector3i(-1, -1, -1), Eigen::Vector3i(2, 2, 2))) { - auto qpoint - = ((qpoints[i] + cell_shift.first.cast()).transpose() * primative_cell.reverse()).transpose().eval(); + Eigen::Vector3d qpoint + = primative_cell.reverse().transpose() * (qpoints[i][j] + cell_shift.first.cast()); // 计算这个点到前两个点所在直线的距离 auto distance = (end - begin).cross(qpoint - begin).norm() / (path.second - path.first).norm(); @@ -61,7 +61,8 @@ void ufo::plot_band(std::string config_file) auto distance3 = (end - begin).norm(); if (distance1 < distance3 + threshold && distance2 < distance3 + threshold) // 如果这个点不在终点处, 或者不排除终点, 则加入 - if (distance2 > threshold || !exclude_endpoint) selected_qpoints.emplace(distance1, i); + if (distance2 > threshold || !exclude_endpoint) + selected_qpoints.emplace(distance1, std::pair{i, j}); } } // 去除非常接近的点 @@ -80,11 +81,13 @@ void ufo::plot_band(std::string config_file) auto calculate_values = [] ( // search_qpoints 的第一个返回值 - const boost::container::flat_map& path, + const boost::container::flat_map>& path, // 每一条连续路径的第一个 q 点的索引 const std::set& path_begin, - // 所有 q 点的数据（需要用到它的频率和权重） - const std::vector& qpoints, + // 各个模式的频率和权重，几个维度分别为：MetaQpointIndex, ModeIndex + const std::vector>& frequency, + // 各个模式的权重，几个维度分别为：MetaQpointIndex, ModeIndex, SubQpointIndex + const std::vector>>& weight, // 用于插值的分辨率和范围 const std::array& resolution, const std::array& frequency_range, @@ -96,7 +99,7 @@ void ufo::plot_band(std::string config_file) auto blend = [&] ( // 两个点的索引 - std::size_t a, std::size_t b, + std::pair a, std::pair b, // 按照连续路径混合还是按照断开的路径混合 bool continuous, // 第一个点占的比例 @@ -105,38 +108,36 @@ void ufo::plot_band(std::string config_file) ) -> std::vector { // 混合得到的频率和权重 - std::vector frequency, weight; + std::vector frequency_result, weight_result; // 如果是连续路径，将每个模式的频率和权重按照比例混合 if (continuous) { - assert(qpoints[a].ModeData.size() == qpoints[b].ModeData.size()); - for (std::size_t i = 0; i < qpoints[a].ModeData.size(); i++) + for (std::size_t i = 0; i < frequency[a.first].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)); + frequency_result.push_back(frequency[a.first][i] * ratio + frequency[b.first][i] * (1 - ratio)); + weight_result.push_back(weight[a.first][i][a.second] * ratio + weight[b.first][i][b.second] * (1 - ratio)); } } // 如果是不连续路径，将每个模式的权重乘以比例，最后相加 else { - for (std::size_t i = 0; i < qpoints[a].ModeData.size(); i++) + for (std::size_t i = 0; i < frequency[a.first].size(); i++) { - frequency.push_back(qpoints[a].ModeData[i].Frequency); - weight.push_back(qpoints[a].ModeData[i].Weight * ratio); + frequency_result.push_back(frequency[a.first][i]); + weight_result.push_back(weight[a.first][i][a.second] * ratio); } - for (std::size_t i = 0; i < qpoints[b].ModeData.size(); i++) + for (std::size_t i = 0; i < frequency[b.first].size(); i++) { - frequency.push_back(qpoints[b].ModeData[i].Frequency); - weight.push_back(qpoints[b].ModeData[i].Weight * (1 - ratio)); + frequency_result.push_back(frequency[b.first][i]); + weight_result.push_back(weight[b.first][i][b.second] * (1 - ratio)); } } std::vector result(resolution); - for (std::size_t i = 0; i < frequency.size(); i++) + for (std::size_t i = 0; i < frequency_result.size(); i++) { - std::ptrdiff_t index = (frequency[i] - frequency_range[0]) / (frequency_range[1] - frequency_range[0]) + std::ptrdiff_t index = (frequency_result[i] - frequency_range[0]) / (frequency_range[1] - frequency_range[0]) * resolution; - if (index >= 0 && index < static_cast(resolution)) result[index] += weight[i]; + if (index >= 0 && index < static_cast(resolution)) result[index] += weight_result[i]; } return result; }; @@ -155,7 +156,8 @@ void ufo::plot_band(std::string config_file) )); else values.push_back(blend ( - std::prev(it)->second, it->second, !path_begin.contains(it->second), + std::prev(it)->second, it->second, + !path_begin.contains(std::distance(path.begin(), it)), (it->first - current_distance) / (it->first - std::prev(it)->first), resolution[1], frequency_range )); @@ -216,19 +218,19 @@ void ufo::plot_band(std::string config_file) f->save(filename, "png"); }; - auto input = YAML::LoadFile(config_file).as(); - auto unfolded_data = biu::deserialize - (biu::read(input.UnfoldedDataFile)); + auto config = YAML::LoadFile(config_file).as(); + auto input = biu::deserialize + (biu::read(config.InputDataFile)); // 搜索画图需要用到的 q 点 // key 到起点的距离，value 为 q 点在 QpointData 中的索引 - boost::container::flat_map path; + boost::container::flat_map> path; // 每一条连续路径的第一个 q 点在 path 中的索引 std::set path_begin; // x 轴的刻度，为 path 中的索引 std::set x_ticks_index; double total_distance = 0; - for (auto& line : input.Qpoints) + for (auto& line : config.Qpoints) { assert(line.size() >= 2); path_begin.insert(path.size()); @@ -237,11 +239,11 @@ void ufo::plot_band(std::string config_file) 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) + input.Primative.Cell, {line[i], line[i + 1]}, + input.Super.Qpoint + | ranges::views::transform([](auto& qpoint) { return qpoint.SubQpoint; }) | ranges::to_vector, - input.ThresholdWhenSearchingQpoints.value_or(0.001), + config.ThresholdWhenSearchingQpoints.value_or(0.001), i != line.size() - 2 ); path.merge @@ -258,39 +260,59 @@ void ufo::plot_band(std::string config_file) // 计算画图的数据 auto values = calculate_values ( - path, path_begin, unfolded_data.QpointData, input.InterpolationResolution, - input.FrequencyRange, total_distance + path, path_begin, + input.Super.Qpoint + | ranges::views::transform([](auto& qpoint) + { + return qpoint.Mode | ranges::views::transform([](auto&& mode) { return mode.Frequency; }) | ranges::to_vector; + }) + | ranges::to_vector, + input.Super.Qpoint + | ranges::views::transform([](auto& qpoint) + { + return qpoint.Mode + | ranges::views::transform([](auto&& mode) + { + return mode.WeightOnUnfold + | ranges::views::transform([&](auto&& weight) + { return weight * mode.WeightOnSelectedAtom.value_or(1) * mode.WeightOnRaman.value_or(1); }) + | ranges::to_vector; + }) + | ranges::to_vector; + }) + | ranges::to_vector, + config.InterpolationResolution, config.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; - auto y_ticks = input.YTicks.value_or(std::vector>{}) + { return path.nth(i)->first / total_distance * config.InterpolationResolution[0]; }) | ranges::to_vector; + auto y_ticks = config.YTicks.value_or(std::vector>{}) | biu::toLvalue | ranges::views::keys | ranges::views::transform([&](auto i) { - return (i - input.FrequencyRange[0]) / (input.FrequencyRange[1] - input.FrequencyRange[0]) - * input.InterpolationResolution[1]; + return (i - config.FrequencyRange[0]) / (config.FrequencyRange[1] - config.FrequencyRange[0]) + * config.InterpolationResolution[1]; }) | ranges::to_vector; - auto y_ticklabels = input.YTicks.value_or(std::vector>{}) + auto y_ticklabels = config.YTicks.value_or(std::vector>{}) | biu::toLvalue | ranges::views::values | ranges::to_vector; - if (input.OutputPictureFile) plot + if (config.OutputPictureFile) plot ( - values, input.OutputPictureFile.value(), - x_ticks, y_ticks, y_ticklabels, input.AspectRatio, input.PictureResolution + values, config.OutputPictureFile.value(), + x_ticks, y_ticks, y_ticklabels, config.AspectRatio, config.PictureResolution ); - if (input.OutputDataFile) - biu::Hdf5file(input.OutputDataFile.value(), true) + if (config.OutputDataFile) + biu::Hdf5file(config.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); + .write("InterpolationResolution", config.InterpolationResolution) + .write("FrequencyRange", config.FrequencyRange); } void ufo::plot_point(std::string config_file) { - struct Input + struct Config { std::string UnfoldedDataFile; // 要画图的 q 点 @@ -317,19 +339,20 @@ void ufo::plot_point(std::string config_file) auto search_qpoints = [] ( const Eigen::Matrix3d& primative_cell, - const Eigen::Vector3d& qpoint, const std::vector& qpoints, + const Eigen::Vector3d& qpoint, const std::vector>& qpoints, double threshold ) { biu::Logger::Guard log(qpoint); // 对于 output 中的每一个点, 检查这个点是否与所寻找的点足够近，如果足够近则返回 - for (std::size_t i = 0; i < qpoints.size(); i++) + for (std::size_t i = 0; i < qpoints.size(); i++) for (std::size_t j = 0; j < qpoints[i].size(); j++) 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())).eval(); - if ((this_qpoint - primative_cell.reverse().transpose() * qpoint).norm() < threshold) return log.rtn(i); + = (primative_cell.reverse().transpose() * (qpoints[i][j] + cell_shift.first.cast())).eval(); + if ((this_qpoint - primative_cell.reverse().transpose() * qpoint).norm() < threshold) + return log.rtn(std::pair(i, j)); } throw std::runtime_error("No q points found"); }; @@ -338,7 +361,8 @@ void ufo::plot_point(std::string config_file) auto calculate_values = [] ( // q 点的数据（需要用到它的频率和权重） - const UnfoldOutput::QpointDataType& qpoint, + const std::vector& frequency, + const std::vector& weight, // 用于插值的分辨率和范围 std::size_t resolution, const std::array& frequency_range @@ -346,12 +370,12 @@ void ufo::plot_point(std::string config_file) { biu::Logger::Guard log; std::vector result(resolution); - for (auto& mode : qpoint.ModeData) + for (std::size_t i = 0; i < frequency.size(); i++) { - double index_double = (mode.Frequency - frequency_range[0]) / (frequency_range[1] - frequency_range[0]) + double index_double = (frequency[i] - 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(resolution)) result[index] += mode.Weight; + if (index >= 0 && index < static_cast(resolution)) result[index] += weight[i]; } return log.rtn(result); }; @@ -387,45 +411,55 @@ void ufo::plot_point(std::string config_file) }; biu::Logger::Guard log; - auto input = YAML::LoadFile(config_file).as(); - auto unfolded_data = biu::deserialize - (biu::read(input.UnfoldedDataFile)); + auto config = YAML::LoadFile(config_file).as(); + auto input = biu::deserialize + (biu::read(config.UnfoldedDataFile)); auto qpoint_index = search_qpoints ( - unfolded_data.PrimativeCell, input.Qpoint, - unfolded_data.QpointData - | ranges::views::transform(&UnfoldOutput::QpointDataType::Qpoint) + input.Primative.Cell, config.Qpoint, + input.Super.Qpoint + | ranges::views::transform([](auto& qpoint) { return qpoint.SubQpoint; }) | ranges::to_vector, - input.ThresholdWhenSearchingQpoints.value_or(0.001) + config.ThresholdWhenSearchingQpoints.value_or(0.001) ); auto values = calculate_values ( - unfolded_data.QpointData[qpoint_index], - input.InterpolationResolution, input.FrequencyRange + input.Super.Qpoint[qpoint_index.first].Mode + | ranges::views::transform([](auto&& mode) { return mode.Frequency; }) + | ranges::to_vector, + input.Super.Qpoint[qpoint_index.first].Mode + | ranges::views::transform([&](auto&& mode) + { + return mode.WeightOnUnfold[qpoint_index.second] + * mode.WeightOnSelectedAtom.value_or(1) * mode.WeightOnRaman.value_or(1); + }) + | ranges::to_vector, + config.InterpolationResolution, config.FrequencyRange ); - auto x_ticks = input.XTicks.value_or(std::vector>{}) + auto x_ticks = config.XTicks.value_or(std::vector>{}) | biu::toLvalue | ranges::views::keys | ranges::views::transform([&](auto i) { - return (i - input.FrequencyRange[0]) / (input.FrequencyRange[1] - input.FrequencyRange[0]) - * input.InterpolationResolution; + return (i - config.FrequencyRange[0]) / (config.FrequencyRange[1] - config.FrequencyRange[0]) + * config.InterpolationResolution; }) | ranges::to_vector; - auto x_ticklabels = input.XTicks.value_or(std::vector>{}) + auto x_ticklabels = config.XTicks.value_or(std::vector>{}) | biu::toLvalue | ranges::views::values | ranges::to_vector; - if (input.OutputPictureFile) plot + if (config.OutputPictureFile) plot ( - values, input.OutputPictureFile.value(), - x_ticks, x_ticklabels, input.AspectRatio, input.PictureResolution + values, config.OutputPictureFile.value(), + x_ticks, x_ticklabels, config.AspectRatio, config.PictureResolution ); - if (input.OutputDataFile) - biu::Hdf5file(input.OutputDataFile.value(), true) + if (config.OutputDataFile) + biu::Hdf5file(config.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]); + .write("InterpolationResolution", config.InterpolationResolution) + .write("FrequencyRange", config.FrequencyRange); + // TODO: pfr + optional + Eigen Matrix = failed + if (config.OutputRawDataFile) + std::ofstream(*config.OutputRawDataFile) << YAML::Node(values); } diff --git a/src/project-to-atom.cpp b/src/project-to-atom.cpp new file mode 100644 index 0000000..e7b82d8 --- /dev/null +++ b/src/project-to-atom.cpp @@ -0,0 +1,37 @@ +# include +# include + +void ufo::project_to_atom(std::string config_file) +{ + struct Config + { + std::set SelectedAtom; + std::string InputFile; + std::string OutputFile; + }; + + biu::Logger::Guard log(config_file); + auto config = YAML::LoadFile(config_file).as(); + auto input = biu::deserialize(biu::read(config.InputFile)); + input.SelectedAtom = config.SelectedAtom; + for (auto& qpoint : input.Super.Qpoint) for (auto& mode : qpoint.Mode) + mode.WeightOnSelectedAtom = mode.EigenVector.rowwise().norm().cwiseProduct + ( + ranges::views::iota(0ul, input.Super.AtomPosition.rows() * 1ul) + | ranges::views::transform([&](auto i) + { return config.SelectedAtom.contains(i) ? 1. : 0.; }) + | ranges::to_vector + | biu::toEigen<> + ).sum() * (input.Super.AtomPosition.rows() * 1. / config.SelectedAtom.size()); + std::ofstream(config.OutputFile, std::ios::binary) << biu::serialize(input); + log.info("Summary:"); + for (auto i_of_qpoint : std::views::iota(0ul, input.Super.Qpoint.size())) + for (auto i_of_mode : std::views::iota(0ul, input.Super.Qpoint[i_of_qpoint].Mode.size())) + if (*input.Super.Qpoint[i_of_qpoint].Mode[i_of_mode].WeightOnSelectedAtom > 0.1) + log.info("{}:{:.2f}:{}:{:.2f}"_f + ( + i_of_qpoint, fmt::join(input.Super.Qpoint[i_of_qpoint].Qpoint, ", "), + i_of_mode, input.Super.Qpoint[i_of_qpoint].Mode[i_of_mode].Frequency, + *input.Super.Qpoint[i_of_qpoint].Mode[i_of_mode].WeightOnSelectedAtom + )); +} diff --git a/src/raman-apply-contribution.cpp b/src/raman-apply-contribution.cpp deleted file mode 100644 index 3fda163..0000000 --- a/src/raman-apply-contribution.cpp +++ /dev/null @@ -1,104 +0,0 @@ -# include - -void ufo::raman_apply_contribution(std::string config_file) -{ - struct Input - { - std::string UnfoldedDataFile; - std::string DisplacementDataFile; - Eigen::Matrix3d OriginalSusceptibility; - std::vector Susceptibilities; - std::string OutputDataFile; - std::array Polarization; - }; - - biu::Logger::Guard log(config_file); - auto input = YAML::LoadFile(config_file).as(); - auto unfolded_data = biu::deserialize - (biu::read(input.UnfoldedDataFile)); - auto displacement_data = biu::deserialize - (biu::read(input.DisplacementDataFile)); - UnfoldOutput output; - - // 整理得到的数据，作放缩 - struct mode_t - { - std::size_t MetaQpointIndex; - std::size_t ModeIndex; - Eigen::Matrix3d RamanTensor; - double Strength; - }; - auto modes = ranges::views::zip(displacement_data.ModeData, input.Susceptibilities) - | ranges::views::transform([&](const auto& data) - { - auto& [mode, susceptibility] = data; - Eigen::Matrix3d raman_tensor = susceptibility / mode.Ratio / displacement_data.MaxDisplacement; - double intensity = (input.Polarization[0].transpose() * raman_tensor * input.Polarization[1]).norm(); - log.debug("get raman tensor: {}"_f(raman_tensor)); - log.debug("get intensity: {}"_f(intensity)); - return mode_t{ mode.MetaQpointIndex, mode.ModeIndex, raman_tensor, intensity }; - }) - | ranges::to_vector; - - // 整理输出文件 - biu::for_each - ( - [&](auto&& i) { output.*i = unfolded_data.*i; }, - std::tuple(&UnfoldOutput::PrimativeCell, &UnfoldOutput::SuperCellTransformation, - &UnfoldOutput::SuperCellMultiplier, &UnfoldOutput::SuperCellDeformation, &UnfoldOutput::SelectedAtoms) - ); - // 将以前的 meta qpoint 和 mode 的索引转换到新的 - std::map meta_qpoint_map; - std::vector> mode_map; - for (auto mode : displacement_data.ModeData) - { - auto old_meta_qpoint_index = mode.MetaQpointIndex; - auto old_mode_index = mode.ModeIndex; - auto new_meta_qpoint_index = meta_qpoint_map.contains(old_meta_qpoint_index) - ? meta_qpoint_map[mode.MetaQpointIndex] - : [&] - { - log.debug("map meta qpoint {} to {}"_f(old_meta_qpoint_index, meta_qpoint_map.size())); - meta_qpoint_map[mode.MetaQpointIndex] = meta_qpoint_map.size(); - output.MetaQpointData.push_back({unfolded_data.MetaQpointData[old_mode_index].Qpoint, {}}); - return meta_qpoint_map.size() - 1; - }(); - auto new_mode_index = mode_map[new_meta_qpoint_index].size(); - log.debug("map mode {} {} to {}"_f(old_meta_qpoint_index, old_mode_index, new_mode_index)); - mode_map[new_meta_qpoint_index][old_mode_index] = new_mode_index; - output.MetaQpointData[new_meta_qpoint_index].ModeData.push_back - (unfolded_data.MetaQpointData[old_meta_qpoint_index].ModeData[old_mode_index]); - } - for (auto old_qpoint_index : displacement_data.QpointIndices) - if (meta_qpoint_map.contains(unfolded_data.QpointData[old_qpoint_index].SourceIndex)) - { - log.debug("map qpoint {} {}"_f(old_qpoint_index, unfolded_data.QpointData[old_qpoint_index].Qpoint)); - auto new_meta_qpoint_index - = meta_qpoint_map[unfolded_data.QpointData[old_qpoint_index].SourceIndex]; - output.QpointData.push_back - ({ - unfolded_data.QpointData[old_qpoint_index].Qpoint, - unfolded_data.QpointData[old_qpoint_index].Source, - new_meta_qpoint_index, - {} - }); - for (std::size_t i = 0; i < unfolded_data.QpointData[old_qpoint_index].ModeData.size(); i++) - if (mode_map[new_meta_qpoint_index].contains(i)) - { - auto new_mode_index = mode_map[new_meta_qpoint_index][i]; - log.debug("map mode {} to {}"_f(i, new_mode_index)); - if (output.QpointData.back().ModeData.size() <= new_mode_index) - { - log.debug("resize to {}"_f(new_mode_index + 1)); - output.QpointData.back().ModeData.resize(new_mode_index + 1); - } - output.QpointData.back().ModeData[new_mode_index] = - { - unfolded_data.QpointData[old_qpoint_index].ModeData[i].Frequency, - unfolded_data.QpointData[old_qpoint_index].ModeData[i].Weight * modes[new_mode_index].Strength - }; - } - } - - std::ofstream(input.OutputDataFile, std::ios::binary) << biu::serialize(output); -} diff --git a/src/raman-create-displacement.cpp b/src/raman-create-displacement.cpp deleted file mode 100644 index 7102778..0000000 --- a/src/raman-create-displacement.cpp +++ /dev/null @@ -1,128 +0,0 @@ -# include - -void ufo::raman_create_displacement(std::string config_file) -{ - struct Input - { - std::string UnfoldedDataFile; - // 搜索位于 Gamma 点的 q 点时，使用的阈值，单位为埃^-1，默认为 0.01 - std::optional ThresholdWhenSearchingQpoints; - // 搜索权重非零的模式时，使用的阈值，默认为 0.01 - std::optional ThresholdWhenSearchingModes; - // 所有原子的符号 - std::vector AtomSymbols; - // 各种原子的质量，单位为原子质量 - std::map AtomMasses; - // 原子最大位移大小，单位为埃 - double MaxDisplacement; - // 超胞，单位为埃 - Eigen::Matrix3d SuperCell; - // 超胞中各个原子的坐标，单位为超胞的格矢 - Eigen::MatrixX3d AtomPositions; - // 输出的POSCAR所在的目录 - std::string OutputPoscarDirectory; - // 输出的数据文件名 - std::string OutputDataFile; - }; - - // 假定同类型的原子一定写在一起 - auto generate_poscar = [] - ( - Eigen::Matrix3d SuperCell, Eigen::MatrixX3d AtomPositions, std::vector 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(); - }; - - biu::Logger::Guard log(config_file); - auto input = YAML::LoadFile(config_file).as(); - auto unfolded_data = biu::deserialize - (biu::read(input.UnfoldedDataFile)); - DisplacementOutput output; - - // 搜索满足条件的模式，找到满足条件的模式后，就将 MetaQpoint 的索引加入到 output 中 - // 之所以使用 MetaQpoint 的索引而不是 Qpoint 的索引，是因为 Qpoint 中可能有指向同一个 MetaQpoint 中模式的不同模式都满足要求 - // 如果写入 Qpoint 的索引，就会重复而增加之后的计算量 - std::set> 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; - log.debug("Found mode {} {} frequency {}"_f(qpoint.Qpoint, i, qpoint.ModeData[i].Frequency)); - selected_modes.insert({qpoint.SourceIndex, i}); - } - } - - // 构造输出数据 - output.AtomMasses = input.AtomSymbols - | ranges::views::transform([&](const auto& symbol) - { return input.AtomMasses.at(symbol); }) - | ranges::to_vector - | biu::toEigen<>; - output.MaxDisplacement = input.MaxDisplacement; - 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(output.AtomMasses.cwiseSqrt().cwiseInverse().rowwise().replicate(3)).eval(); - // 归一化 - mode_data.Ratio = input.MaxDisplacement / atom_movement.rowwise().norm().maxCoeff(); - mode_data.AtomMovement = atom_movement * mode_data.Ratio; - log.debug("Write mode {} {} {}"_f(i, j, mode_data.AtomMovement.rowwise().norm().transpose())); - } - - // 输出 - std::ofstream(input.OutputDataFile, std::ios::binary) << biu::serialize(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 - ); -} diff --git a/src/raman-extract.cpp b/src/raman-extract.cpp deleted file mode 100644 index 1b80857..0000000 --- a/src/raman-extract.cpp +++ /dev/null @@ -1,42 +0,0 @@ -# include - -void ufo::raman_extract(std::vector files) -{ - biu::Logger::Guard log(files); - std::vector electricity_tensors; - - for (const auto& file : files) - { - auto in_stream = std::ifstream(file); - std::string line; - // search line containing: MACROSCOPIC STATIC DIELECTRIC TENSOR - while (std::getline(in_stream, line)) - { - if (line.find("MACROSCOPIC STATIC DIELECTRIC TENSOR") != std::string::npos) - { - log.debug("find in {}"_f(file)); - // skip 1 line - std::getline(in_stream, line); - electricity_tensors.emplace_back(); - for (std::size_t i = 0; i < 3; i++) for (std::size_t j = 0; j < 3; j++) - in_stream >> electricity_tensors.back()(i, j); - log.debug("read tensor {}"_f(electricity_tensors.back())); - break; - } - } - // test if the file is read correctly - if (!in_stream) throw std::runtime_error("Error reading file: {}"_f(file)); - } - - // output the result - for (auto e: electricity_tensors) std::cout << -R"(- - [ {}, {}, {} ] - - [ {}, {}, {} ] - - [ {}, {}, {} ] -)"_f - ( - e(0, 0), e(0, 1), e(0, 2), - e(1, 0), e(1, 1), e(1, 2), - e(2, 0), e(2, 1), e(2, 2) - ); -} diff --git a/src/raman.cpp b/src/raman.cpp new file mode 100644 index 0000000..d833981 --- /dev/null +++ b/src/raman.cpp @@ -0,0 +1,199 @@ + +# include + +namespace ufo +{ + struct RamanData + { + struct ModeType + { + std::size_t QpointIndex; + std::size_t ModeIndex; + double Ratio; + using serialize = zpp::bits::members<3>; + }; + std::vector Mode; + double MaxDisplacement; + enum { Primative, Super } Cell; + using serialize = zpp::bits::members<3>; + }; + + void raman_create_displacement(std::string config_file) + { + struct Config + { + // 要计算的是原胞还是超胞 + decltype(RamanData::Cell) Cell; + // 要计算的模式，总是假定是 gamma 点的模式 + std::map> SelectedModes; + // 原子最大位移大小，单位为埃 + double MaxDisplacement; + // 输出的POSCAR所在的目录 + std::string OutputPoscarDirectory; + // 输入的数据文件名 + std::string InputDataFile; + std::string OutputDataFile; + }; + + auto generate_poscar = [] + ( + Eigen::Matrix3d SuperCell, Eigen::MatrixX3d AtomPosition, + std::vector> AtomType + ) + { + 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; + } + ss << "{}\n"_f(ranges::accumulate + ( + AtomType | ranges::views::transform([](auto&& atom) { return atom.first; }), + ""s, [](auto&& a, auto&& b) { return a + " " + b; } + )); + ss << "{}\n"_f(ranges::accumulate + ( + AtomType | ranges::views::transform([](auto&& atom) { return atom.second; }), + ""s, [](auto&& a, auto&& b) { return a + " " + std::to_string(b); } + )); + ss << "Direct\n"; + for (const auto& position : AtomPosition.rowwise()) + { + for (std::size_t i = 0; i < 3; i++) ss << position(i) << " "; + ss << std::endl; + } + return ss.str(); + }; + + biu::Logger::Guard log(config_file); + auto config = YAML::LoadFile(config_file).as(); + auto input = biu::deserialize + (biu::read(config.InputDataFile)); + RamanData output; + output.MaxDisplacement = config.MaxDisplacement; + output.Cell = config.Cell; + + auto process = [&](auto& cell) + { + std::size_t i_of_poscar = 0; + auto mass = cell.AtomType + | ranges::views::transform([&](auto&& atom) + { return ranges::views::repeat_n(input.AtomMass[atom.first], atom.second); }) + | ranges::views::join | ranges::to_vector | biu::toEigen<>; + for (auto i_of_qpoint : config.SelectedModes | ranges::views::keys) + for (auto i_of_mode : config.SelectedModes[i_of_qpoint]) + { + // 假定虚部总是为零 + auto atom_movement = cell.Qpoint[i_of_qpoint].Mode[i_of_mode].EigenVector.real() + .cwiseProduct(mass.cwiseSqrt().cwiseInverse().rowwise().replicate(3)).eval(); + // 归一化 + auto ratio = config.MaxDisplacement / atom_movement.rowwise().norm().maxCoeff(); + atom_movement *= ratio; + // 输出 + auto path = "{}/{}"_f(config.OutputPoscarDirectory, i_of_poscar); + std::filesystem::create_directories(path); + std::ofstream("{}/POSCAR"_f(path)) << generate_poscar + ( + cell.Cell, + cell.AtomPosition + atom_movement * cell.Cell.inverse(), + cell.AtomType + ); + output.Mode.push_back({i_of_qpoint, i_of_mode, ratio}); + log.debug("Write mode {} {} {}"_f(i_of_qpoint, i_of_mode, atom_movement.rowwise().norm().eval())); + i_of_poscar++; + } + }; + if (config.Cell == RamanData::Primative) process(input.Primative); + else process(input.Super); + + std::filesystem::create_directories("{}/{}"_f(config.OutputPoscarDirectory, "origin")); + std::ofstream("{}/origin/POSCAR"_f(config.OutputDataFile)) << generate_poscar + (input.Super.Cell, input.Super.AtomPosition, input.Super.AtomType); + std::ofstream(config.OutputDataFile, std::ios::binary) << biu::serialize(output); + } + + void raman_extract(std::vector files) + { + biu::Logger::Guard log(files); + std::vector electricity_tensors; + + for (const auto& file : files) + { + auto in_stream = std::ifstream(file); + std::string line; + // search line containing: MACROSCOPIC STATIC DIELECTRIC TENSOR + while (std::getline(in_stream, line)) + { + if (line.find("MACROSCOPIC STATIC DIELECTRIC TENSOR") != std::string::npos) + { + log.debug("find in {}"_f(file)); + // skip 1 line + std::getline(in_stream, line); + electricity_tensors.emplace_back(); + for (std::size_t i = 0; i < 3; i++) for (std::size_t j = 0; j < 3; j++) + in_stream >> electricity_tensors.back()(i, j); + log.debug("read tensor {}"_f(electricity_tensors.back())); + break; + } + } + // test if the file is read correctly + if (!in_stream) throw std::runtime_error("Error reading file: {}"_f(file)); + } + + // output the result + for (auto e: electricity_tensors) std::cout << +R"(- - [ {}, {}, {} ] + - [ {}, {}, {} ] + - [ {}, {}, {} ] +)"_f + ( + e(0, 0), e(0, 1), e(0, 2), + e(1, 0), e(1, 1), e(1, 2), + e(2, 0), e(2, 1), e(2, 2) + ); + } + + void raman_apply_contribution(std::string config_file) + { + struct Config + { + Eigen::Matrix3d OriginalSusceptibility; + std::vector Susceptibilities; + std::array Polarization; + std::string InputDataFile; + std::string RamanInputDataFile; + std::string OutputDataFile; + }; + + biu::Logger::Guard log(config_file); + auto config = YAML::LoadFile(config_file).as(); + auto input = biu::deserialize + (biu::read(config.InputDataFile)); + auto raman_input = biu::deserialize + (biu::read(config.RamanInputDataFile)); + + input.RamanPolarization = config.Polarization; + auto process = [&](auto& cell) + { + for (auto&& [i_of_mode, mode] : ranges::views::enumerate(raman_input.Mode)) + { + auto&& _ = cell.Qpoint[mode.QpointIndex].Mode[mode.ModeIndex]; + Eigen::Matrix3d raman_tensor = (config.Susceptibilities[i_of_mode] - config.OriginalSusceptibility) + / mode.Ratio / raman_input.MaxDisplacement; + _.RamanTensor = raman_tensor | biu::fromEigen; + _.WeightOnRaman = config.Polarization[0].transpose() * raman_tensor * config.Polarization[1]; + log.info("{}:{:.2f}:{}:{:.2f} {}"_f + ( + mode.QpointIndex, fmt::join(cell.Qpoint[mode.QpointIndex].Qpoint, ", "), + mode.ModeIndex, _.Frequency, *_.WeightOnRaman + )); + } + }; + if (raman_input.Cell == RamanData::Primative) process(input.Primative); + else process(input.Super); + + std::ofstream(config.OutputDataFile, std::ios::binary) << biu::serialize(input); + } +} diff --git a/src/unfold.cpp b/src/unfold.cpp index de20a03..bd3128c 100644 --- a/src/unfold.cpp +++ b/src/unfold.cpp @@ -2,6 +2,7 @@ # include # include # include +# include void ufo::unfold(std::string config_file) { @@ -13,11 +14,8 @@ void ufo::unfold(std::string config_file) // 只要第一部分取得足够多, 那么单胞中原子的状态就可以完全被这些平面波描述. // 将超胞中原子的运动状态投影到这些基矢上, 计算出投影的系数, 就可以将超胞的原子运动状态分解到单胞中的多个 q 点上. - struct Input + struct Config { - // 单胞的三个格矢，每行表示一个格矢的坐标，单位为埃 - Eigen::Matrix3d PrimativeCell; - // 单胞到超胞的格矢转换时用到的矩阵 // SuperCellMultiplier 是一个三维列向量且各个元素都是整数，表示单胞在各个方向扩大到多少倍之后，可以得到和超胞一样的体积 // SuperCellDeformation 是一个行列式为 1 的矩阵，它表示经过 SuperCellMultiplier 扩大后，还需要怎样的变换才能得到超胞 @@ -35,40 +33,32 @@ void ufo::unfold(std::string config_file) // 在单胞内取几个平面波的基矢 Eigen::Vector PrimativeCellBasisNumber; - // 超胞中原子的坐标，每行表示一个原子的坐标，单位为超胞的格矢 - Eigen::MatrixX3d AtomPositionBySuperCell; + // 单胞的 phonopy 输出的 phonopy.yaml，用来读入单胞的晶格、原子坐标、原子类型、原子质量 + std::string PrimativePhonopy; + // 单胞的 phonopy 输出的 band.hdf5 或 qpoint.hdf5，用来读入 q 点和振动模式 + std::string PrimativeQpoint; + // 同上，但是是超胞里的结果 + std::string SuperPhonopy; + std::string SuperQpoint; - // 从 band.hdf5 读入 QpointData - std::optional QpointDataInputFile; - - // 输出到哪些文件 - struct QpointDataOutputFileType - { - std::string Filename; - - // 如果指定，则将结果投影到那些原子上 - std::optional> SelectedAtoms; - - // 默认输出为 zpp 文件，如果指定为 true，则输出为 yaml 文件 - std::optional OutputAsYaml; - }; - std::vector QpointDataOutputFile; + std::string OutputFile; }; - // 从文件中读取 QpointData - auto read_qpoint_data = [](std::string filename) + // 从文件中读取 q 点数据 + // data 为 CommonData::Primative 或 CommonData::Super + // 返回值为原子类型和原子质量的对应关系 + auto read_qpoint = [](std::string phonopy_file, std::string qpoint_file, auto& data) { - // 读入原始数据 - // phonopy 的输出有两种可能 - // 直接指定计算的 q 点时，frequency 是 2 维，这时第一个维度是 q 点，第二个维度是不同模式 - // 计算能带时，frequency 是 3 维，相比于二维的情况多了第一个维度，表示 q 点所在路径 - // qpoint 或 path，以及 eigenvector 也有类似的变化 - // eigenvector 是三维或四维的数组，后两个维度分别表示原子运动和模式（而不是模式和原子）， - // 因为后两个维度的尺寸总是一样的（模式个数等于原子坐标个数），非常容易搞错 + // phonopy 的输出有两种可能。 + // 直接指定计算的 q 点时，frequency 是 2 维，这时第一个维度是 q 点，第二个维度是不同模式。 + // 计算能带时，frequency 是 3 维，相比于二维的情况多了第一个维度，表示 q 点所在路径。 + // qpoint 或 path，以及 eigenvector 也有类似的变化。 + // 除此以外，eigenvector 后两个维度分别指示模式的特征向量的各个维度和各个模式（而不是各个模式和特征向量的各个维度）， + // 因为后两个维度的尺寸总是一样的（模式个数等于原子坐标个数），非常容易搞错。 std::vector> qpoint; std::vector> frequency; - std::vector>> eigenvector_vector; - auto file = biu::Hdf5file(filename); + std::vector>>> eigenvector_vector; + auto file = biu::Hdf5file(qpoint_file); if (file.File.getDataSet("/frequency").getDimensions().size() == 2) file.read("/frequency", frequency) @@ -78,7 +68,7 @@ void ufo::unfold(std::string config_file) { std::vector>> temp_path; std::vector>> temp_frequency; - std::vector>>> temp_eigenvector_vector; + std::vector>>>> temp_eigenvector_vector; file.read("/frequency", temp_frequency) .read("/eigenvector", temp_eigenvector_vector) .read("/path", temp_path); @@ -87,27 +77,47 @@ void ufo::unfold(std::string config_file) 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 qpoint_data(number_of_qpoints); - for (std::size_t i = 0; i < number_of_qpoints; i++) + // 整理并写入得到结果 + auto number_of_qpoints = frequency.size(), number_of_modes = frequency[0].size(); + data.Qpoint.resize(number_of_qpoints); + for (auto i : std::views::iota(0u, number_of_qpoints)) { - 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++) + data.Qpoint[i].Qpoint = qpoint[i] | biu::toEigen<>; + data.Qpoint[i].Mode.resize(number_of_modes); + for (auto j : std::views::iota(0u, number_of_modes)) { - 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; + data.Qpoint[i].Mode[j].Frequency = frequency[i][j]; + auto number_of_atoms = number_of_modes / 3; + Eigen::MatrixX3cd eigenvector(number_of_atoms, 3); + for (auto k : std::views::iota(0u, number_of_atoms)) + for (auto l : std::views::iota(0u, 3u)) + eigenvector(k, l) = eigenvector_vector[i][k * 3 + l][j]; // 原则上讲，需要对读入的原子运动状态作相位转换, 使得它们与我们的约定一致(对超胞周期性重复)，但这个转换 phonopy 已经做了 // 这里还要需要做归一化处理 (指将数据简单地作为向量处理的归一化) - qpoint_data[i].ModeData[j].AtomMovement = eigenvectors / eigenvectors.norm(); + data.Qpoint[i].Mode[j].EigenVector = eigenvector / eigenvector.norm(); } } - return qpoint_data; + + // 读取并写入其它数据 + YAML::Node phonopy = YAML::LoadFile(phonopy_file); + data.Cell = phonopy["unit_cell"]["lattice"].as, 3>>() | biu::toEigen<>; + auto points = phonopy["points"].as>(); + data.AtomType = points + | ranges::views::transform([](const YAML::Node& point) { return point["symbol"].as(); }) + | ranges::views::chunk_by(std::ranges::equal_to{}) + | ranges::views::transform([](const auto& chunk) { return std::pair{chunk[0], chunk.size()}; }) + | ranges::to_vector; + data.AtomPosition = points + | ranges::views::transform([](const YAML::Node& point) + { return point["coordinates"].as>(); }) + | ranges::to_vector + | biu::toEigen<>; + return points + | ranges::views::transform([](const YAML::Node& point) + { return std::pair(point["symbol"].as(), point["mass"].as()); }) + | ranges::views::chunk_by(std::ranges::equal_to{}) + | ranges::views::transform([](const auto& chunk) { return chunk[0]; }) + | ranges::to>; }; // 构建基 @@ -151,31 +161,25 @@ void ufo::unfold(std::string config_file) auto construct_projection_coefficient = [] ( const std::vector>& basis, - // 实际上只需要其中的 AtomMovement - const std::vector& qpoint_data, + const std::vector>& modes, std::atomic& number_of_finished_modes ) { - // 将所有的模式取出，组成一个一维数组，稍后并行计算 - std::vector> mode_data; - for (auto& qpoint : qpoint_data) for (auto& mode : qpoint.ModeData) - mode_data.emplace_back(mode.AtomMovement); // 第一层下标对应不同模式, 第二层下标对应这个模式在反折叠后的 q 点(sub qpoint) - std::vector> projection_coefficient(mode_data.size()); + std::vector> projection_coefficient(modes.size()); // 对每个模式并行 std::transform ( - std::execution::par_unseq, mode_data.begin(), mode_data.end(), + std::execution::par, modes.begin(), modes.end(), projection_coefficient.begin(), [&](const auto& mode_data) { // 这里, mode_data 和 projection_coefficient 均指对应于一个模式的数据 std::vector projection_coefficient(basis.size()); - for (std::size_t i_of_sub_qpoint = 0; i_of_sub_qpoint < basis.size(); i_of_sub_qpoint++) + for (auto i_of_sub_qpoint : std::views::iota(0u, basis.size())) // 对于 basis 中, 对应于单胞倒格子的部分, 以及对应于不同方向的部分, 分别求内积, 然后求模方和 - for (std::size_t i_of_basis = 0; i_of_basis < basis[i_of_sub_qpoint].size(); i_of_basis++) + for (auto i_of_basis : std::views::iota(0u, basis[i_of_sub_qpoint].size())) projection_coefficient[i_of_sub_qpoint] += - (basis[i_of_sub_qpoint][i_of_basis].transpose().conjugate() * mode_data.get()) - .array().abs2().sum(); + (basis[i_of_sub_qpoint][i_of_basis].transpose().conjugate() * mode_data.get()).array().abs2().sum(); // 如果是严格地将向量分解到一组完备的基矢上, 那么不需要对计算得到的权重再做归一化处理 // 但这里并不是这样一个严格的概念. 因此对分解到各个 sub qpoint 上的权重做归一化处理 auto sum = ranges::accumulate(projection_coefficient, 0.); @@ -184,162 +188,81 @@ void ufo::unfold(std::string config_file) return projection_coefficient; } ); - // 将计算得到的投影系数重新组装成三维数组 - // 第一维是 meta qpoint，第二维是模式，第三维是 sub qpoint - std::vector>> projection_coefficient_output; - for - ( - std::size_t i_of_meta_qpoint = 0, num_of_mode_manipulated = 0; - i_of_meta_qpoint < qpoint_data.size(); - i_of_meta_qpoint++, num_of_mode_manipulated += qpoint_data[i_of_meta_qpoint].ModeData.size() - ) - projection_coefficient_output.emplace_back - ( - projection_coefficient.begin() + num_of_mode_manipulated, - projection_coefficient.begin() + num_of_mode_manipulated + qpoint_data[i_of_meta_qpoint].ModeData.size() - ); - return projection_coefficient_output; - }; - - // 组装输出，即将投影系数应用到原始数据上 - auto construct_output = [] - ( - const Input& input, - const std::vector>>& projection_coefficient, - const std::vector& qpoint_data, - const std::optional>& selected_atoms - ) - { - UnfoldOutput output; - output.PrimativeCell = input.PrimativeCell; - output.SuperCellMultiplier = input.SuperCellMultiplier; - output.SuperCellDeformation = input.SuperCellDeformation; - output.SelectedAtoms = selected_atoms; - output.MetaQpointData = qpoint_data; - for (std::size_t i_of_meta_qpoint = 0; i_of_meta_qpoint < qpoint_data.size(); i_of_meta_qpoint++) - { - // 如果需要投影到特定的原子上，需要先计算当前 meta qpoint 的不同模式的投影系数 - std::optional> 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(qpoint_data[i_of_meta_qpoint].ModeData[i_of_mode].AtomMovement.rows()) - / selected_atoms->size(); - } - } - - for - ( - auto [diff_of_sub_qpoint_by_reciprocal_modified_super_cell, i_of_sub_qpoint] - : biu::sequence(input.SuperCellMultiplier) - ) - { - auto& _ = output.QpointData.emplace_back(); - /* - SubQpointByReciprocalModifiedSuperCell = XyzOfDiffOfSubQpointByReciprocalModifiedSuperCell + - MetaQpointByReciprocalModifiedSuperCell; - SubQpoint = SubQpointByReciprocalModifiedSuperCell.transpose() * ReciprocalModifiedSuperCell; - SubQpoint = SubQpointByReciprocalPrimativeCell.transpose() * ReciprocalPrimativeCell; - ReciprocalModifiedSuperCell = ModifiedSuperCell.inverse().transpose(); - ReciprocalPrimativeCell = PrimativeCell.inverse().transpose(); - ModifiedSuperCell = SuperCellMultiplier.asDiagonal() * PrimativeCell; - MetaQpoint = MetaQpointByReciprocalModifiedSuperCell.transpose() * ReciprocalModifiedSuperCell; - MetaQpoint = MetaQpointByReciprocalSuperCell.transpose() * ReciprocalSuperCell; - ReciprocalSuperCell = SuperCell.inverse().transpose(); - ModifiedSuperCell = SuperCellDeformation * SuperCell; - SuperCell = SuperCellMultiplier.asDiagonal() * PrimativeCell; - 整理可以得到: - SubQpointByReciprocalPrimativeCell = SuperCellMultiplier.asDiagonal().inverse() * - (XyzOfDiffOfSubQpointByReciprocalModifiedSuperCell + - SuperCellDeformation.inverse() * MetaQpointByReciprocalSuperCell); - 但注意到, 这样得到的 SubQpoint 可能不在 ReciprocalPrimativeCell 中 - (当 SuperCellDeformation 不是单位矩阵时, 边界附近的一两条 SubQpoint 会出现这种情况). - 解决办法是, 在赋值时, 仅取 SubQpointByReciprocalPrimativeCell 的小数部分. - */ - auto sub_qpoint_by_reciprocal_primative_cell = - ( - input.SuperCellMultiplier.cast().cwiseInverse().asDiagonal() - * ( - diff_of_sub_qpoint_by_reciprocal_modified_super_cell.cast() - + input.SuperCellDeformation.inverse() * qpoint_data[i_of_meta_qpoint].Qpoint - ) - ).eval(); - _.Qpoint = sub_qpoint_by_reciprocal_primative_cell.array() - - sub_qpoint_by_reciprocal_primative_cell.array().floor(); - _.Source = qpoint_data[i_of_meta_qpoint].Qpoint; - _.SourceIndex = i_of_meta_qpoint; - - for (std::size_t i_of_mode = 0; i_of_mode < qpoint_data[i_of_meta_qpoint].ModeData.size(); i_of_mode++) - { - auto& __ = _.ModeData.emplace_back(); - __.Frequency = qpoint_data[i_of_meta_qpoint].ModeData[i_of_mode].Frequency; - __.Weight = projection_coefficient[i_of_meta_qpoint][i_of_mode][i_of_sub_qpoint]; - if (selected_atoms) - __.Weight *= projection_coefficient_on_atoms.value()[i_of_mode]; - } - } - } - return output; + return projection_coefficient; }; biu::Logger::Guard log; - log.info("Reading input file... "); - auto input = YAML::LoadFile(config_file).as(); - auto qpoint_data = read_qpoint_data(input.QpointDataInputFile.value_or("band.hdf5")); + + log.info("Reading input file..."); + auto config = YAML::LoadFile(config_file).as(); + CommonData output; + output.AtomMass = read_qpoint + (config.PrimativePhonopy, config.PrimativeQpoint, output.Primative); + output.AtomMass.merge(read_qpoint + (config.SuperPhonopy, config.SuperQpoint, output.Super)); + output.Super.CellDeformation = config.SuperCellDeformation; + output.Super.CellMultiplier = config.SuperCellMultiplier; log.info("Done."); - std::clog << "Constructing basis... " << std::flush; - + log.info("Constructing basis..."); auto basis = construct_basis ( - input.PrimativeCell, input.SuperCellMultiplier, - input.PrimativeCellBasisNumber, - input.AtomPositionBySuperCell - * (input.SuperCellDeformation * input.SuperCellMultiplier.cast().asDiagonal() * input.PrimativeCell) + output.Primative.Cell, output.Super.CellMultiplier, + config.PrimativeCellBasisNumber, + output.Super.AtomPosition + * (output.Super.CellDeformation * output.Super.CellMultiplier.cast().asDiagonal() * output.Primative.Cell) ); - std::clog << "Done." << std::endl; + log.info("Done."); std::clog << "Calculating projection coefficient... " << std::flush; - // 用来在屏幕上输出进度的计数器和线程 - std::atomic 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 finished; - std::thread print_thread([&] + // 将所有模式放到一维来处理 { - while (true) + auto modes = output.Super.Qpoint + | ranges::views::transform([](const auto& qpoint) + { return qpoint.Mode | ranges::views::transform([](const auto& mode) + { return std::cref(mode.EigenVector); }); }) + | ranges::views::join + | ranges::to_vector; + std::atomic number_of_finished_modes(0); + std::thread print_thread([&] { - 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(); + while (true) + { + std::osyncstream(std::clog) + << "\rCalculating projection coefficient... ({}/{})"_f(number_of_finished_modes, modes.size()) + << std::flush; + std::this_thread::sleep_for(100ms); + if (number_of_finished_modes == modes.size()) break; + } + }); + auto projection_coefficient = construct_projection_coefficient + (basis, modes, number_of_finished_modes); + std::size_t i_of_modes = 0; + for (auto& qpoint : output.Super.Qpoint) for (auto& mode : qpoint.Mode) + mode.WeightOnUnfold = projection_coefficient[i_of_modes++]; + 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(output); - } - std::clog << "Done." << std::endl; + log.info("Writing data... "); + std::ofstream(config.OutputFile, std::ios::binary) << biu::serialize(output); + log.info("Done."); + + log.info("Summary:"); + + for (auto i_of_qpoint : std::views::iota(0u, output.Super.Qpoint.size())) + for (auto i_of_mode : std::views::iota(0u, output.Super.Qpoint[i_of_qpoint].Mode.size())) + for + ( + auto i_of_sub_qpoint + : std::views::iota(0u, output.Super.Qpoint[i_of_qpoint].SubQpoint.size()) + ) + if (output.Super.Qpoint[i_of_qpoint].Mode[i_of_mode].WeightOnUnfold[i_of_sub_qpoint] > 0.01) + log.info("{}:{:.2f}:{}:{:.2f} -> {:.2f} {:.2f}"_f + ( + i_of_qpoint, fmt::join(output.Super.Qpoint[i_of_qpoint].Qpoint, ", "), + i_of_mode, output.Super.Qpoint[i_of_qpoint].Mode[i_of_mode].Frequency, + fmt::join(output.Super.Qpoint[i_of_qpoint].SubQpoint[i_of_sub_qpoint], ", "), + output.Super.Qpoint[i_of_qpoint].Mode[i_of_mode].WeightOnUnfold[i_of_sub_qpoint] + )); }