diff --git a/CMakeLists.txt b/CMakeLists.txt index c7472f3..8fbf50c 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -14,17 +14,22 @@ find_package(Matplot++ REQUIRED) find_package(biu REQUIRED) find_package(Threads REQUIRED) -add_executable(ufo src/fold.cpp src/unfold.cpp src/plot.cpp src/raman-create-displacement.cpp src/main.cpp) +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) 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}") +endif() install(TARGETS ufo RUNTIME DESTINATION ${CMAKE_INSTALL_BINDIR}) get_property(ImportedTargets DIRECTORY "${CMAKE_SOURCE_DIR}" PROPERTY IMPORTED_TARGETS) message("Imported targets: ${ImportedTargets}") message("List of compile features: ${CMAKE_CXX_COMPILE_FEATURES}") +message("CMake build type: ${CMAKE_BUILD_TYPE}") include(CTest) add_test(NAME fold COMMAND ufo fold ${PROJECT_SOURCE_DIR}/test/fold/config.yaml) diff --git a/include/ufo.hpp b/include/ufo.hpp index 4e48475..eaa4827 100644 --- a/include/ufo.hpp +++ b/include/ufo.hpp @@ -16,6 +16,7 @@ namespace ufo 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); // 许多函数都需要用到这个，所以写到头文件中 @@ -71,4 +72,21 @@ namespace ufo 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/main.cpp b/src/main.cpp index 1291b32..d082956 100644 --- a/src/main.cpp +++ b/src/main.cpp @@ -3,11 +3,17 @@ int main(int argc, const char** argv) { using namespace biu::literals; - if (argc != 3) throw std::runtime_error("Usage: {} task config.yaml"_f(argv[0])); if (argv[1] == "fold"s) ufo::fold(argv[2]); else if (argv[1] == "unfold"s) ufo::unfold(argv[2]); else if (argv[1] == "raman-create-displacement"s) ufo::raman_create_displacement(argv[2]); - else if (argv[1] == "raman-apply-contribution"s); + else if (argv[1] == "raman-extract"s) ufo::raman_extract + ( + std::vector(argv + 2, argv + argc) + | biu::toLvalue + | ranges::views::transform([](const char* arg) { return std::string(arg); }) + | ranges::to_vector + ); + else if (argv[1] == "raman-apply-contribution"s) ufo::raman_apply_contribution(argv[2]); else if (argv[1] == "plot-band"s) ufo::plot_band(argv[2]); else if (argv[1] == "plot-point"s) ufo::plot_point(argv[2]); else throw std::runtime_error("Unknown task: {}"_f(argv[1])); diff --git a/src/raman-apply-contribution.cpp b/src/raman-apply-contribution.cpp new file mode 100644 index 0000000..3fda163 --- /dev/null +++ b/src/raman-apply-contribution.cpp @@ -0,0 +1,104 @@ +# 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 index 226bdf8..7102778 100644 --- a/src/raman-create-displacement.cpp +++ b/src/raman-create-displacement.cpp @@ -24,18 +24,6 @@ void ufo::raman_create_displacement(std::string config_file) // 输出的数据文件名 std::string OutputDataFile; }; - struct Output - { - struct ModeData_t - { - std::size_t MetaQpointIndex; - std::size_t ModeIndex; - // 每个原子的位移，单位为埃 - Eigen::MatrixX3d AtomMovement; - }; - std::vector ModeData; - using serialize = zpp::bits::members<1>; - }; // 假定同类型的原子一定写在一起 auto generate_poscar = [] @@ -70,10 +58,11 @@ void ufo::raman_create_displacement(std::string config_file) 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)); - Output output; + DisplacementOutput output; // 搜索满足条件的模式，找到满足条件的模式后，就将 MetaQpoint 的索引加入到 output 中 // 之所以使用 MetaQpoint 的索引而不是 Qpoint 的索引，是因为 Qpoint 中可能有指向同一个 MetaQpoint 中模式的不同模式都满足要求 @@ -90,11 +79,18 @@ void ufo::raman_create_displacement(std::string config_file) 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(); @@ -102,18 +98,12 @@ void ufo::raman_create_displacement(std::string config_file) mode_data.ModeIndex = j; // 未归一化的位移, 假定虚部总是为零 auto atom_movement = - unfolded_data.MetaQpointData[i].ModeData[j].AtomMovement.real().cwiseProduct - ( - ( - input.AtomSymbols - | ranges::views::transform([&](const auto& symbol) - { return input.AtomMasses.at(symbol); }) - | ranges::to_vector - | biu::toEigen<> - ).cwiseSqrt().cwiseInverse().rowwise().replicate(3) - ).eval(); + unfolded_data.MetaQpointData[i].ModeData[j].AtomMovement.real() + .cwiseProduct(output.AtomMasses.cwiseSqrt().cwiseInverse().rowwise().replicate(3)).eval(); // 归一化 - mode_data.AtomMovement = atom_movement / atom_movement.rowwise().norm().maxCoeff() * input.MaxDisplacement; + 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())); } // 输出 diff --git a/src/raman-extract.cpp b/src/raman-extract.cpp new file mode 100644 index 0000000..1b80857 --- /dev/null +++ b/src/raman-extract.cpp @@ -0,0 +1,42 @@ +# 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) + ); +}