add raman support

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)

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<std::string> 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_t> ModeData;
+    Eigen::Vector3d AtomMasses;
+    double MaxDisplacement;
+    std::vector<std::size_t> QpointIndices;
+    using serialize = zpp::bits::members<4>;
+  };
 }

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<const char*>(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]));

src/raman-apply-contribution.cpp

@@ -0,0 +1,104 @@
+# include <ufo.hpp>
+
+void ufo::raman_apply_contribution(std::string config_file)
+{
+  struct Input
+  {
+    std::string UnfoldedDataFile;
+    std::string DisplacementDataFile;
+    Eigen::Matrix3d OriginalSusceptibility;
+    std::vector<Eigen::Matrix3d> Susceptibilities;
+    std::string OutputDataFile;
+    std::array<Eigen::Vector3d, 2> Polarization;
+  };
+
+  biu::Logger::Guard log(config_file);
+  auto input = YAML::LoadFile(config_file).as<Input>();
+  auto unfolded_data = biu::deserialize<UnfoldOutput>
+    (biu::read<std::byte>(input.UnfoldedDataFile));
+  auto displacement_data = biu::deserialize<DisplacementOutput>
+    (biu::read<std::byte>(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<std::size_t, std::size_t> meta_qpoint_map;
+  std::vector<std::map<std::size_t, std::size_t>> 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<char>(output);
+}

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_t> 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<Input>();
   auto unfolded_data = biu::deserialize<UnfoldOutput>
     (biu::read<std::byte>(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
+      unfolded_data.MetaQpointData[i].ModeData[j].AtomMovement.real()
-      (
+        .cwiseProduct(output.AtomMasses.cwiseSqrt().cwiseInverse().rowwise().replicate(3)).eval();
-        (
-          input.AtomSymbols
-            | ranges::views::transform([&](const auto& symbol)
-              { return input.AtomMasses.at(symbol); })
-            | ranges::to_vector
-            | biu::toEigen<>
-        ).cwiseSqrt().cwiseInverse().rowwise().replicate(3)
-      ).eval();
     // 归一化
-    mode_data.AtomMovement = atom_movement / atom_movement.rowwise().norm().maxCoeff() * input.MaxDisplacement;
+    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()));
   }
 
   // 输出

src/raman-extract.cpp

@@ -0,0 +1,42 @@
+# include <ufo.hpp>
+
+void ufo::raman_extract(std::vector<std::string> files)
+{
+  biu::Logger::Guard log(files);
+  std::vector<Eigen::Matrix3d> 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)
+    );
+}