diff --git a/.envrc b/.envrc index 8392d15..2a81941 100644 --- a/.envrc +++ b/.envrc @@ -1 +1 @@ -use flake \ No newline at end of file +use flake nixos#ufo diff --git a/.gitignore b/.gitignore index 9ae3fb0..a9037bf 100644 --- a/.gitignore +++ b/.gitignore @@ -1,18 +1,4 @@ -.vscode .direnv -main -plot -out.yaml -.ccls-cache +build +.vscode .cache -build/ -test/* -!test/14.2.2.4.unfold.yaml -!test/14.2.2.5.unfold.yaml -!test/14.2.2.5.plot.yaml -!test/14.2.4.4.unfold.yaml -!test/14.2.4.4.plot.yaml -!test/14.2.5.5.unfold.yaml -!test/14.2.5.5.plot.yaml -!test/14.2.6.4.unfold.yaml -!test/14.2.6.4.plot.yaml diff --git a/CMakeLists.txt b/CMakeLists.txt index 0d84b3f..c7472f3 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -1,36 +1,31 @@ cmake_minimum_required(VERSION 3.14) -project(ufo VERSION 0.0.0 LANGUAGES CXX) +project(ufo VERSION 0 LANGUAGES CXX) enable_testing() include(GNUInstallDirs) -set_property(GLOBAL PROPERTY CXX_STANDARD 23) -set_property(GLOBAL PROPERTY CXX_STANDARD_REQUIRED ON) -set_property(GLOBAL PROPERTY CXX_EXTENSIONS OFF) - if(NOT CMAKE_BUILD_TYPE AND NOT CMAKE_CONFIGURATION_TYPES) message("Setting build type to 'Release' as none was specified.") set(CMAKE_BUILD_TYPE Release CACHE STRING "Choose the type of build." FORCE) set_property(CACHE CMAKE_BUILD_TYPE PROPERTY STRINGS "Debug" "Release" "MinSizeRel" "RelWithDebInfo") endif() -find_package(yaml-cpp REQUIRED) -find_package(Eigen3 REQUIRED) -find_package(fmt REQUIRED) -find_package(concurrencpp REQUIRED) -find_package(HighFive REQUIRED) find_package(TBB REQUIRED) -find_package(glad REQUIRED) -# find_package(Matplot++ REQUIRED COMPONENTS matplot_opengl) find_package(Matplot++ REQUIRED) -find_path(ZPP_BITS_INCLUDE_DIR zpp_bits.h REQUIRED) +find_package(biu REQUIRED) +find_package(Threads REQUIRED) -add_executable(ufo src/solver.cpp src/fold.cpp src/unfold.cpp src/plot.cpp src/main.cpp) -target_include_directories(ufo PRIVATE ${PROJECT_SOURCE_DIR}/include ${ZPP_BITS_INCLUDE_DIR}) -target_link_libraries(ufo PRIVATE - yaml-cpp Eigen3::Eigen fmt::fmt concurrencpp::concurrencpp HighFive_HighFive TBB::tbb Matplot++::matplot - Matplot++::matplot_opengl) -# target_compile_definitions(ufo PRIVATE $<$:_GLIBCXX_DEBUG>) +add_executable(ufo src/fold.cpp src/unfold.cpp src/plot.cpp src/raman-create-displacement.cpp src/main.cpp) +target_include_directories(ufo PRIVATE ${PROJECT_SOURCE_DIR}/include) +target_link_libraries(ufo PRIVATE TBB::tbb Matplot++::matplot biu::biu) +target_compile_features(ufo PRIVATE cxx_std_23) +target_compile_options(ufo PRIVATE -fexperimental-library) + +install(TARGETS ufo RUNTIME DESTINATION ${CMAKE_INSTALL_BINDIR}) get_property(ImportedTargets DIRECTORY "${CMAKE_SOURCE_DIR}" PROPERTY IMPORTED_TARGETS) message("Imported targets: ${ImportedTargets}") message("List of compile features: ${CMAKE_CXX_COMPILE_FEATURES}") + +include(CTest) +add_test(NAME fold COMMAND ufo fold ${PROJECT_SOURCE_DIR}/test/fold/config.yaml) + diff --git a/README.md b/README.md deleted file mode 100644 index 6b2a8c7..0000000 --- a/README.md +++ /dev/null @@ -1,17 +0,0 @@ -This is a program to approximately unfold phonon spectra. - -Currently, there is no documentation, and some not core features are in other branches - and implemented in a temporary way, e.g., directly writing data in source file. - -I will finish the documentation and merge these features into main branch ASAP. - -The program is packaged in a standard way of cmake, so that you can build it just like other cmake projects: - -```bash -mkdir build -cd build -cmake .. -make -``` - -and view source files for input and output format. diff --git a/default.nix b/default.nix new file mode 100644 index 0000000..3d2364d --- /dev/null +++ b/default.nix @@ -0,0 +1,11 @@ +{ + stdenv, cmake, pkg-config, version ? null, + tbb, matplotplusplus, biu +}: stdenv.mkDerivation +{ + name = "ufo"; + src = ./.; + buildInputs = [ tbb matplotplusplus biu ]; + nativeBuildInputs = [ cmake pkg-config ]; + doCheck = true; +} diff --git a/doc/README.md b/doc/README.md new file mode 100644 index 0000000..aa5d6bd --- /dev/null +++ b/doc/README.md @@ -0,0 +1,45 @@ +分为几个功能： + +* fold：根据要计算的单胞的 q 点路径，计算超胞中对应的 q 点路径，生成的路径再交给 phonopy 计算。 +* unfold：根据 phonopy 计算的结果，将超胞的结果展开到单胞中。 +* plot：对计算结果画图。 + +主要的输入输出格式均为 yaml。对于数据特别大的情况，也可以从 hdf5 中读取一部分数据或者将一部分数据写入到 hdf5 文件中。 + +# fold + +## 输入 + +```yaml +# 三个整数组成的向量，表示从单胞到超胞，三个晶格矢量的倍数 +# 必写 +SuperCellMultiplier: [2, 2, 2] +# 一个变换矩阵，表明超胞经历了怎样的扭曲。 +# 可选，默认值为单位矩阵 +SuperCellDeformation: [[1, 0, 0], [0, 1, 0], [0, 0, 1]] +# 一个由三个浮点数组成的向量，表示考虑的 q 点 +# 必写 +Qpoints: + - [0, 0, 0] + - [0.1, 0, 0] + - [0.2, 0, 0] + - [0.3, 0, 0] + - [0.4, 0, 0] + - [0.5, 0, 0] +# 一个 DataFile 类型的对象，表明输出结果到哪个文件 +# 必写 +OutputFile: +``` + +## 输出 + +```yaml +# 得到的 q 点坐标 +Qpoints: + - [0, 0, 0] + - [0.1, 0, 0] + - [0.2, 0, 0] + - [0.3, 0, 0] + - [0.4, 0, 0] + - [0.5, 0, 0] +``` diff --git a/flake.lock b/flake.lock deleted file mode 100644 index 93c15c3..0000000 --- a/flake.lock +++ /dev/null @@ -1,1370 +0,0 @@ -{ - "nodes": { - "aagl": { - "inputs": { - "flake-compat": "flake-compat", - "nixpkgs": [ - "nixos", - "nixpkgs" - ] - }, - "locked": { - "lastModified": 1696252780, - "narHash": "sha256-sQEjVzzstiaNLyiFJ19EMwwbDSSNDyQZIbPiLonlDCQ=", - "owner": "ezKEa", - "repo": "aagl-gtk-on-nix", - "rev": "0c9d93bdb311f7948f9fb0e98d869316d78eec12", - "type": "github" - }, - "original": { - "owner": "ezKEa", - "repo": "aagl-gtk-on-nix", - "type": "github" - } - }, - "deploy-rs": { - "inputs": { - "flake-compat": "flake-compat_2", - "nixpkgs": [ - "nixos", - "nixpkgs" - ], - "utils": "utils" - }, - "locked": { - "lastModified": 1695052866, - "narHash": "sha256-agn7F9Oww4oU6nPiw+YiYI9Xb4vOOE73w8PAoBRP4AA=", - "owner": "serokell", - "repo": "deploy-rs", - "rev": "e3f41832680801d0ee9e2ed33eb63af398b090e9", - "type": "github" - }, - "original": { - "owner": "serokell", - "repo": "deploy-rs", - "type": "github" - } - }, - "dguibert-nur-packages": { - "inputs": { - "dwl-src": "dwl-src", - "dwm-src": "dwm-src", - "emacs-overlay": "emacs-overlay", - "emacs-src": "emacs-src", - "flake-parts": "flake-parts", - "flake-utils": "flake-utils_2", - "mako-src": "mako-src", - "nix": "nix", - "nixpkgs": [ - "nixos", - "nixpkgs" - ], - "st-src": "st-src" - }, - "locked": { - "lastModified": 1693549169, - "narHash": "sha256-1mgWe0BH63pEOI2MaclIxGyQbyh0e/k+PbjCasuNfS4=", - "owner": "CHN-beta", - "repo": "dguibert-nur-packages", - "rev": "3d4419864595ee05e82a391206d9f8499e321265", - "type": "github" - }, - "original": { - "owner": "CHN-beta", - "repo": "dguibert-nur-packages", - "type": "github" - } - }, - "dwl-src": { - "flake": false, - "locked": { - "lastModified": 1680697480, - "narHash": "sha256-vaxCJB4ykNVwZNKau42Xu+D39A5mXgjfyKgOE2Frsa4=", - "owner": "dguibert", - "repo": "dwl", - "rev": "844b04d23db2de161a058bec477ef7ff772ce298", - "type": "github" - }, - "original": { - "owner": "dguibert", - "ref": "pu", - "repo": "dwl", - "type": "github" - } - }, - "dwm-src": { - "flake": false, - "locked": { - "lastModified": 1649608053, - "narHash": "sha256-+v2wof46XG/n85zVY82GIuKgzr0K9s1apPfaDSyuqos=", - "owner": "dguibert", - "repo": "dwm", - "rev": "306e1790145849abd30ef6c81051755a6813ff30", - "type": "github" - }, - "original": { - "owner": "dguibert", - "ref": "pu", - "repo": "dwm", - "type": "github" - } - }, - "emacs-overlay": { - "inputs": { - "flake-utils": "flake-utils", - "nixpkgs": [ - "nixos", - "dguibert-nur-packages", - "nixpkgs" - ], - "nixpkgs-stable": "nixpkgs-stable" - }, - "locked": { - "lastModified": 1691144360, - "narHash": "sha256-DPaalQfXWXbyiuSqbk0dWKusniMWQSzAu0e36xU6rbA=", - "owner": "nix-community", - "repo": "emacs-overlay", - "rev": "4ba15d6f4310459e6da08dcd4d3df7f4d102bdf0", - "type": "github" - }, - "original": { - "owner": "nix-community", - "repo": "emacs-overlay", - "type": "github" - } - }, - "emacs-src": { - "flake": false, - "locked": { - "lastModified": 1685570069, - "narHash": "sha256-lp6HgksMTfUboBm/BBLHFwiZ1aQnzHEAOhRL8I9RHOg=", - "owner": "emacs-mirror", - "repo": "emacs", - "rev": "17c7915ab947ebeec6ea5ad3eb4cad1f24d5d4fc", - "type": "github" - }, - "original": { - "owner": "emacs-mirror", - "ref": "emacs-29", - "repo": "emacs", - "type": "github" - } - }, - "flake-compat": { - "flake": false, - "locked": { - "lastModified": 1673956053, - "narHash": "sha256-4gtG9iQuiKITOjNQQeQIpoIB6b16fm+504Ch3sNKLd8=", - "owner": "edolstra", - "repo": "flake-compat", - "rev": "35bb57c0c8d8b62bbfd284272c928ceb64ddbde9", - "type": "github" - }, - "original": { - "owner": "edolstra", - "repo": "flake-compat", - "type": "github" - } - }, - "flake-compat_2": { - "flake": false, - "locked": { - "lastModified": 1668681692, - "narHash": "sha256-Ht91NGdewz8IQLtWZ9LCeNXMSXHUss+9COoqu6JLmXU=", - "owner": "edolstra", - "repo": "flake-compat", - "rev": "009399224d5e398d03b22badca40a37ac85412a1", - "type": "github" - }, - "original": { - "owner": "edolstra", - "repo": "flake-compat", - "type": "github" - } - }, - "flake-compat_3": { - "flake": false, - "locked": { - "lastModified": 1673956053, - "narHash": "sha256-4gtG9iQuiKITOjNQQeQIpoIB6b16fm+504Ch3sNKLd8=", - "owner": "edolstra", - "repo": "flake-compat", - "rev": "35bb57c0c8d8b62bbfd284272c928ceb64ddbde9", - "type": "github" - }, - "original": { - "owner": "edolstra", - "repo": "flake-compat", - "type": "github" - } - }, - "flake-compat_4": { - "flake": false, - "locked": { - "lastModified": 1673956053, - "narHash": "sha256-4gtG9iQuiKITOjNQQeQIpoIB6b16fm+504Ch3sNKLd8=", - "owner": "edolstra", - "repo": "flake-compat", - "rev": "35bb57c0c8d8b62bbfd284272c928ceb64ddbde9", - "type": "github" - }, - "original": { - "owner": "edolstra", - "repo": "flake-compat", - "type": "github" - } - }, - "flake-compat_5": { - "flake": false, - "locked": { - "lastModified": 1673956053, - "narHash": "sha256-4gtG9iQuiKITOjNQQeQIpoIB6b16fm+504Ch3sNKLd8=", - "owner": "edolstra", - "repo": "flake-compat", - "rev": "35bb57c0c8d8b62bbfd284272c928ceb64ddbde9", - "type": "github" - }, - "original": { - "owner": "edolstra", - "repo": "flake-compat", - "type": "github" - } - }, - "flake-compat_6": { - "flake": false, - "locked": { - "lastModified": 1673956053, - "narHash": "sha256-4gtG9iQuiKITOjNQQeQIpoIB6b16fm+504Ch3sNKLd8=", - "owner": "edolstra", - "repo": "flake-compat", - "rev": "35bb57c0c8d8b62bbfd284272c928ceb64ddbde9", - "type": "github" - }, - "original": { - "owner": "edolstra", - "repo": "flake-compat", - "type": "github" - } - }, - "flake-parts": { - "inputs": { - "nixpkgs-lib": "nixpkgs-lib" - }, - "locked": { - "lastModified": 1682984683, - "narHash": "sha256-fSMthG+tp60AHhNmaHc4StT3ltfHkQsJtN8GhfLWmtI=", - "owner": "hercules-ci", - "repo": "flake-parts", - "rev": "86684881e184f41aa322e653880e497b66429f3e", - "type": "github" - }, - "original": { - "id": "flake-parts", - "type": "indirect" - } - }, - "flake-parts_2": { - "inputs": { - "nixpkgs-lib": "nixpkgs-lib_2" - }, - "locked": { - "lastModified": 1685662779, - "narHash": "sha256-cKDDciXGpMEjP1n6HlzKinN0H+oLmNpgeCTzYnsA2po=", - "owner": "hercules-ci", - "repo": "flake-parts", - "rev": "71fb97f0d875fd4de4994dfb849f2c75e17eb6c3", - "type": "github" - }, - "original": { - "owner": "hercules-ci", - "repo": "flake-parts", - "type": "github" - } - }, - "flake-parts_3": { - "inputs": { - "nixpkgs-lib": [ - "nixos", - "nixpak", - "nixpkgs" - ] - }, - "locked": { - "lastModified": 1696343447, - "narHash": "sha256-B2xAZKLkkeRFG5XcHHSXXcP7To9Xzr59KXeZiRf4vdQ=", - "owner": "hercules-ci", - "repo": "flake-parts", - "rev": "c9afaba3dfa4085dbd2ccb38dfade5141e33d9d4", - "type": "github" - }, - "original": { - "owner": "hercules-ci", - "repo": "flake-parts", - "type": "github" - } - }, - "flake-parts_4": { - "inputs": { - "nixpkgs-lib": [ - "nixos", - "nixpak", - "hercules-ci-effects", - "hercules-ci-agent", - "nixpkgs" - ] - }, - "locked": { - "lastModified": 1688466019, - "narHash": "sha256-VeM2akYrBYMsb4W/MmBo1zmaMfgbL4cH3Pu8PGyIwJ0=", - "owner": "hercules-ci", - "repo": "flake-parts", - "rev": "8e8d955c22df93dbe24f19ea04f47a74adbdc5ec", - "type": "github" - }, - "original": { - "owner": "hercules-ci", - "repo": "flake-parts", - "type": "github" - } - }, - "flake-utils": { - "inputs": { - "systems": "systems" - }, - "locked": { - "lastModified": 1689068808, - "narHash": "sha256-6ixXo3wt24N/melDWjq70UuHQLxGV8jZvooRanIHXw0=", - "owner": "numtide", - "repo": "flake-utils", - "rev": "919d646de7be200f3bf08cb76ae1f09402b6f9b4", - "type": "github" - }, - "original": { - "owner": "numtide", - "repo": "flake-utils", - "type": "github" - } - }, - "flake-utils-plus": { - "inputs": { - "flake-utils": [ - "nixos", - "nur-xddxdd", - "flake-utils" - ] - }, - "locked": { - "lastModified": 1696331477, - "narHash": "sha256-YkbRa/1wQWdWkVJ01JvV+75KIdM37UErqKgTf0L54Fk=", - "owner": "gytis-ivaskevicius", - "repo": "flake-utils-plus", - "rev": "bfc53579db89de750b25b0c5e7af299e0c06d7d3", - "type": "github" - }, - "original": { - "owner": "gytis-ivaskevicius", - "repo": "flake-utils-plus", - "type": "github" - } - }, - "flake-utils_2": { - "inputs": { - "systems": "systems_2" - }, - "locked": { - "lastModified": 1685518550, - "narHash": "sha256-o2d0KcvaXzTrPRIo0kOLV0/QXHhDQ5DTi+OxcjO8xqY=", - "owner": "numtide", - "repo": "flake-utils", - "rev": "a1720a10a6cfe8234c0e93907ffe81be440f4cef", - "type": "github" - }, - "original": { - "owner": "numtide", - "repo": "flake-utils", - "type": "github" - } - }, - "flake-utils_3": { - "locked": { - "lastModified": 1659877975, - "narHash": "sha256-zllb8aq3YO3h8B/U0/J1WBgAL8EX5yWf5pMj3G0NAmc=", - "owner": "numtide", - "repo": "flake-utils", - "rev": "c0e246b9b83f637f4681389ecabcb2681b4f3af0", - "type": "github" - }, - "original": { - "owner": "numtide", - "repo": "flake-utils", - "type": "github" - } - }, - "flake-utils_4": { - "inputs": { - "systems": "systems_4" - }, - "locked": { - "lastModified": 1689068808, - "narHash": "sha256-6ixXo3wt24N/melDWjq70UuHQLxGV8jZvooRanIHXw0=", - "owner": "numtide", - "repo": "flake-utils", - "rev": "919d646de7be200f3bf08cb76ae1f09402b6f9b4", - "type": "github" - }, - "original": { - "owner": "numtide", - "repo": "flake-utils", - "type": "github" - } - }, - "flake-utils_5": { - "inputs": { - "systems": "systems_5" - }, - "locked": { - "lastModified": 1681202837, - "narHash": "sha256-H+Rh19JDwRtpVPAWp64F+rlEtxUWBAQW28eAi3SRSzg=", - "owner": "numtide", - "repo": "flake-utils", - "rev": "cfacdce06f30d2b68473a46042957675eebb3401", - "type": "github" - }, - "original": { - "owner": "numtide", - "repo": "flake-utils", - "type": "github" - } - }, - "flake-utils_6": { - "locked": { - "lastModified": 1638122382, - "narHash": "sha256-sQzZzAbvKEqN9s0bzWuYmRaA03v40gaJ4+iL1LXjaeI=", - "owner": "numtide", - "repo": "flake-utils", - "rev": "74f7e4319258e287b0f9cb95426c9853b282730b", - "type": "github" - }, - "original": { - "owner": "numtide", - "repo": "flake-utils", - "type": "github" - } - }, - "flake-utils_7": { - "inputs": { - "systems": "systems_6" - }, - "locked": { - "lastModified": 1694529238, - "narHash": "sha256-zsNZZGTGnMOf9YpHKJqMSsa0dXbfmxeoJ7xHlrt+xmY=", - "owner": "numtide", - "repo": "flake-utils", - "rev": "ff7b65b44d01cf9ba6a71320833626af21126384", - "type": "github" - }, - "original": { - "owner": "numtide", - "repo": "flake-utils", - "type": "github" - } - }, - "flake-utils_8": { - "inputs": { - "systems": "systems_7" - }, - "locked": { - "lastModified": 1685518550, - "narHash": "sha256-o2d0KcvaXzTrPRIo0kOLV0/QXHhDQ5DTi+OxcjO8xqY=", - "owner": "numtide", - "repo": "flake-utils", - "rev": "a1720a10a6cfe8234c0e93907ffe81be440f4cef", - "type": "github" - }, - "original": { - "owner": "numtide", - "repo": "flake-utils", - "type": "github" - } - }, - "haskell-flake": { - "locked": { - "lastModified": 1684780604, - "narHash": "sha256-2uMZsewmRn7rRtAnnQNw1lj0uZBMh4m6Cs/7dV5YF08=", - "owner": "srid", - "repo": "haskell-flake", - "rev": "74210fa80a49f1b6f67223debdbf1494596ff9f2", - "type": "github" - }, - "original": { - "owner": "srid", - "ref": "0.3.0", - "repo": "haskell-flake", - "type": "github" - } - }, - "hercules-ci-agent": { - "inputs": { - "flake-parts": "flake-parts_4", - "haskell-flake": "haskell-flake", - "nixpkgs": "nixpkgs_3" - }, - "locked": { - "lastModified": 1688568579, - "narHash": "sha256-ON0M56wtY/TIIGPkXDlJboAmuYwc73Hi8X9iJGtxOhM=", - "owner": "hercules-ci", - "repo": "hercules-ci-agent", - "rev": "367dd8cd649b57009a6502e878005a1e54ad78c5", - "type": "github" - }, - "original": { - "id": "hercules-ci-agent", - "type": "indirect" - } - }, - "hercules-ci-effects": { - "inputs": { - "flake-parts": [ - "nixos", - "nixpak", - "flake-parts" - ], - "hercules-ci-agent": "hercules-ci-agent", - "nixpkgs": [ - "nixos", - "nixpak", - "nixpkgs" - ] - }, - "locked": { - "lastModified": 1695684520, - "narHash": "sha256-yORqGB0i1OtEf9MOCCT2BIbOd8txPZn216CM+ylMmhY=", - "owner": "hercules-ci", - "repo": "hercules-ci-effects", - "rev": "91fae5824f5f1199f61693c6590b4a89abaed9d7", - "type": "github" - }, - "original": { - "owner": "hercules-ci", - "repo": "hercules-ci-effects", - "type": "github" - } - }, - "home-manager": { - "inputs": { - "nixpkgs": [ - "nixos", - "nixpkgs" - ] - }, - "locked": { - "lastModified": 1695108154, - "narHash": "sha256-gSg7UTVtls2yO9lKtP0yb66XBHT1Fx5qZSZbGMpSn2c=", - "owner": "nix-community", - "repo": "home-manager", - "rev": "07682fff75d41f18327a871088d20af2710d4744", - "type": "github" - }, - "original": { - "owner": "nix-community", - "ref": "release-23.05", - "repo": "home-manager", - "type": "github" - } - }, - "impermanence": { - "locked": { - "lastModified": 1694622745, - "narHash": "sha256-z397+eDhKx9c2qNafL1xv75lC0Q4nOaFlhaU1TINqb8=", - "owner": "nix-community", - "repo": "impermanence", - "rev": "e9643d08d0d193a2e074a19d4d90c67a874d932e", - "type": "github" - }, - "original": { - "owner": "nix-community", - "repo": "impermanence", - "type": "github" - } - }, - "lmix": { - "inputs": { - "nixpkgs": [ - "nixos", - "nixpkgs" - ], - "nurl": "nurl", - "utils": "utils_2" - }, - "locked": { - "lastModified": 1693749173, - "narHash": "sha256-+ACUbkspARXLzgxxgiAbINOPYa2+oiR4JhNimpVP9wU=", - "owner": "CHN-beta", - "repo": "lmix", - "rev": "73b61e4318432e50c4b742edbeae82f6ad2fe69a", - "type": "github" - }, - "original": { - "owner": "CHN-beta", - "repo": "lmix", - "type": "github" - } - }, - "lowdown-src": { - "flake": false, - "locked": { - "lastModified": 1653759997, - "narHash": "sha256-mvWCdd79cNh8VK9wfHsMMF2XFYaHd9DhdlwSa33uW+M=", - "owner": "kristapsdz", - "repo": "lowdown", - "rev": "edca6ce6d5336efb147321a43c47a698de41bb7c", - "type": "github" - }, - "original": { - "owner": "kristapsdz", - "repo": "lowdown", - "type": "github" - } - }, - "mako-src": { - "flake": false, - "locked": { - "lastModified": 1684537028, - "narHash": "sha256-y0gp33fCX/qYdexGceDu3ZEUM/hCU2if9BFCSHqZ11g=", - "owner": "emersion", - "repo": "mako", - "rev": "0e111547b17dd8d2733c48fc4082cdb419d7688d", - "type": "github" - }, - "original": { - "owner": "emersion", - "ref": "master", - "repo": "mako", - "type": "github" - } - }, - "napalm": { - "inputs": { - "flake-utils": "flake-utils_3", - "nixpkgs": [ - "nixos", - "nixpkgs" - ] - }, - "locked": { - "lastModified": 1693989153, - "narHash": "sha256-gx39Y3opGB25+44OjM+h1bdJyzgLD963va8ULGYlbhM=", - "owner": "nix-community", - "repo": "napalm", - "rev": "a8215ccf1c80070f51a92771f3bc637dd9b9f7ee", - "type": "github" - }, - "original": { - "owner": "nix-community", - "repo": "napalm", - "type": "github" - } - }, - "nix": { - "inputs": { - "flake-compat": "flake-compat_3", - "lowdown-src": "lowdown-src", - "nixpkgs": [ - "nixos", - "dguibert-nur-packages", - "nixpkgs" - ], - "nixpkgs-regression": "nixpkgs-regression" - }, - "locked": { - "lastModified": 1692264617, - "narHash": "sha256-ZOf28XhU8Fvir3biY58zmIaShiUP/s4GwB4ruuQkMlc=", - "owner": "dguibert", - "repo": "nix", - "rev": "fa563827bdc7f1c939173b1036fb121e43f947c7", - "type": "github" - }, - "original": { - "owner": "dguibert", - "ref": "pu", - "repo": "nix", - "type": "github" - } - }, - "nix-alien": { - "inputs": { - "flake-compat": "flake-compat_4", - "flake-utils": "flake-utils_4", - "nix-index-database": [ - "nixos", - "nix-index-database" - ], - "nixpkgs": "nixpkgs_2" - }, - "locked": { - "lastModified": 1695714965, - "narHash": "sha256-uukcDCyFOIMo5vJWJbLJk2phHZtJ1DE7YrypSV48gII=", - "owner": "thiagokokada", - "repo": "nix-alien", - "rev": "a948cf76e084f4ac770793c6ff9c57ad8b8c099f", - "type": "github" - }, - "original": { - "owner": "thiagokokada", - "repo": "nix-alien", - "type": "github" - } - }, - "nix-index-database": { - "inputs": { - "nixpkgs": [ - "nixos", - "nixpkgs" - ] - }, - "locked": { - "lastModified": 1696131323, - "narHash": "sha256-Y47r8Jo+9rs+XUWHcDPZtkQs6wFeZ24L4CQTfVwE+vY=", - "owner": "Mic92", - "repo": "nix-index-database", - "rev": "031d4b22505fdea47bd53bfafad517cd03c26a4f", - "type": "github" - }, - "original": { - "owner": "Mic92", - "repo": "nix-index-database", - "type": "github" - } - }, - "nix-vscode-extensions": { - "inputs": { - "flake-compat": "flake-compat_5", - "flake-utils": "flake-utils_5", - "nixpkgs": [ - "nixos", - "nixpkgs" - ] - }, - "locked": { - "lastModified": 1693358717, - "narHash": "sha256-OYGe2Yay1QoodZZmvPYBFGAoTrRfyKLzFs2vON4gRek=", - "owner": "nix-community", - "repo": "nix-vscode-extensions", - "rev": "50c4bce16b93e7ca8565d51fafabc05e9f0515da", - "type": "github" - }, - "original": { - "owner": "nix-community", - "repo": "nix-vscode-extensions", - "rev": "50c4bce16b93e7ca8565d51fafabc05e9f0515da", - "type": "github" - } - }, - "nixd": { - "inputs": { - "flake-parts": "flake-parts_2", - "nixpkgs": [ - "nixos", - "nixpkgs" - ] - }, - "locked": { - "lastModified": 1695137077, - "narHash": "sha256-wJ8EpYjsqrR4GFAF67wJKmZd4q86KuODWAag4acQL5Q=", - "owner": "nix-community", - "repo": "nixd", - "rev": "e8f144ca50fe71e74d247e5308ae7ce122f0a0e6", - "type": "github" - }, - "original": { - "owner": "nix-community", - "repo": "nixd", - "type": "github" - } - }, - "nixos": { - "inputs": { - "aagl": "aagl", - "deploy-rs": "deploy-rs", - "dguibert-nur-packages": "dguibert-nur-packages", - "home-manager": "home-manager", - "impermanence": "impermanence", - "lmix": "lmix", - "napalm": "napalm", - "nix-alien": "nix-alien", - "nix-index-database": "nix-index-database", - "nix-vscode-extensions": "nix-vscode-extensions", - "nixd": "nixd", - "nixos-cn": "nixos-cn", - "nixpak": "nixpak", - "nixpkgs": "nixpkgs_4", - "nixpkgs-unstable": "nixpkgs-unstable", - "nur": "nur", - "nur-xddxdd": "nur-xddxdd", - "pnpm2nix-nzbr": "pnpm2nix-nzbr", - "qchem": "qchem", - "sops-nix": "sops-nix", - "touchix": "touchix" - }, - "locked": { - "lastModified": 1696743278, - "narHash": "sha256-ckBxRe2wjUiUzDQWs/WM+x0HgnRUjoZiHzO2c408dFg=", - "owner": "CHN-beta", - "repo": "nixos", - "rev": "5ffdec57c08307cc2589ec57169d72cb68cf7c8b", - "type": "github" - }, - "original": { - "owner": "CHN-beta", - "repo": "nixos", - "type": "github" - } - }, - "nixos-cn": { - "inputs": { - "flake-utils": "flake-utils_6", - "nixpkgs": [ - "nixos", - "nixpkgs" - ] - }, - "locked": { - "lastModified": 1682818384, - "narHash": "sha256-l8jh9BQj6nfjPDYGyrZkZwX1GaOqBX+pBHU+7fFZU3w=", - "owner": "nixos-cn", - "repo": "flakes", - "rev": "2d475ec68cca251ef6c6c69a9224db5c264c5e5b", - "type": "github" - }, - "original": { - "owner": "nixos-cn", - "repo": "flakes", - "type": "github" - } - }, - "nixpak": { - "inputs": { - "flake-parts": "flake-parts_3", - "hercules-ci-effects": "hercules-ci-effects", - "nixpkgs": [ - "nixos", - "nixpkgs" - ] - }, - "locked": { - "lastModified": 1696478570, - "narHash": "sha256-Zqktub0f4M8K0jDHFYaTwsGUddkH3UqHU0NNfGJmIKY=", - "owner": "nixpak", - "repo": "nixpak", - "rev": "271e01d3912c5c622ca7fa99d63d790bea980de0", - "type": "github" - }, - "original": { - "owner": "nixpak", - "repo": "nixpak", - "type": "github" - } - }, - "nixpkgs": { - "locked": { - "lastModified": 1689850295, - "narHash": "sha256-fUYf6WdQlhd2H+3aR8jST5dhFH1d0eE22aes8fNIfyk=", - "owner": "nixos", - "repo": "nixpkgs", - "rev": "5df4d78d54f7a34e9ea1f84a22b4fd9baebc68d0", - "type": "github" - }, - "original": { - "owner": "nixos", - "ref": "nixos-unstable", - "repo": "nixpkgs", - "type": "github" - } - }, - "nixpkgs-lib": { - "locked": { - "dir": "lib", - "lastModified": 1682879489, - "narHash": "sha256-sASwo8gBt7JDnOOstnps90K1wxmVfyhsTPPNTGBPjjg=", - "owner": "NixOS", - "repo": "nixpkgs", - "rev": "da45bf6ec7bbcc5d1e14d3795c025199f28e0de0", - "type": "github" - }, - "original": { - "dir": "lib", - "owner": "NixOS", - "ref": "nixos-unstable", - "repo": "nixpkgs", - "type": "github" - } - }, - "nixpkgs-lib_2": { - "locked": { - "dir": "lib", - "lastModified": 1685564631, - "narHash": "sha256-8ywr3AkblY4++3lIVxmrWZFzac7+f32ZEhH/A8pNscI=", - "owner": "NixOS", - "repo": "nixpkgs", - "rev": "4f53efe34b3a8877ac923b9350c874e3dcd5dc0a", - "type": "github" - }, - "original": { - "dir": "lib", - "owner": "NixOS", - "ref": "nixos-unstable", - "repo": "nixpkgs", - "type": "github" - } - }, - "nixpkgs-regression": { - "locked": { - "lastModified": 1643052045, - "narHash": "sha256-uGJ0VXIhWKGXxkeNnq4TvV3CIOkUJ3PAoLZ3HMzNVMw=", - "owner": "NixOS", - "repo": "nixpkgs", - "rev": "215d4d0fd80ca5163643b03a33fde804a29cc1e2", - "type": "github" - }, - "original": { - "owner": "NixOS", - "repo": "nixpkgs", - "rev": "215d4d0fd80ca5163643b03a33fde804a29cc1e2", - "type": "github" - } - }, - "nixpkgs-stable": { - "locked": { - "lastModified": 1690927903, - "narHash": "sha256-D5gCaCROnjEKDOel//8TO/pOP87pAEtT0uT8X+0Bj/U=", - "owner": "NixOS", - "repo": "nixpkgs", - "rev": "bd836ac5e5a7358dea73cb74a013ca32864ccb86", - "type": "github" - }, - "original": { - "owner": "NixOS", - "ref": "nixos-23.05", - "repo": "nixpkgs", - "type": "github" - } - }, - "nixpkgs-unstable": { - "locked": { - "lastModified": 1696511470, - "narHash": "sha256-eSIw2JSXKHjZwkZGM2d6Cj2/ega6+cqoWZ+EM6lWY04=", - "owner": "CHN-beta", - "repo": "nixpkgs", - "rev": "c3cf1cffb3e119bd166c51ca86605c1cf12988d9", - "type": "github" - }, - "original": { - "owner": "CHN-beta", - "ref": "nixos-unstable", - "repo": "nixpkgs", - "type": "github" - } - }, - "nixpkgs_2": { - "locked": { - "lastModified": 1692007866, - "narHash": "sha256-X8w0vPZjZxMm68VCwh/BHDoKRGp+BgzQ6w7Nkif6IVM=", - "owner": "NixOS", - "repo": "nixpkgs", - "rev": "de2b8ddf94d6cc6161b7659649594c79bd66c13b", - "type": "github" - }, - "original": { - "owner": "NixOS", - "ref": "nixpkgs-unstable", - "repo": "nixpkgs", - "type": "github" - } - }, - "nixpkgs_3": { - "locked": { - "lastModified": 1688322751, - "narHash": "sha256-eW62dC5f33oKZL7VWlomttbUnOTHrAbte9yNUNW8rbk=", - "owner": "NixOS", - "repo": "nixpkgs", - "rev": "0fbe93c5a7cac99f90b60bdf5f149383daaa615f", - "type": "github" - }, - "original": { - "owner": "NixOS", - "ref": "nixos-unstable", - "repo": "nixpkgs", - "type": "github" - } - }, - "nixpkgs_4": { - "locked": { - "lastModified": 1696511131, - "narHash": "sha256-IIhn6F8D26Kix77guTW/4KdpwBzpSHJ3mjG1C8FAwHc=", - "owner": "CHN-beta", - "repo": "nixpkgs", - "rev": "1bac8e4beb5b30458994710236b9db265829327b", - "type": "github" - }, - "original": { - "owner": "CHN-beta", - "ref": "nixos-23.05", - "repo": "nixpkgs", - "type": "github" - } - }, - "nixpkgs_5": { - "locked": { - "lastModified": 1696604326, - "narHash": "sha256-YXUNI0kLEcI5g8lqGMb0nh67fY9f2YoJsILafh6zlMo=", - "owner": "NixOS", - "repo": "nixpkgs", - "rev": "87828a0e03d1418e848d3dd3f3014a632e4a4f64", - "type": "github" - }, - "original": { - "owner": "NixOS", - "ref": "nixos-unstable", - "repo": "nixpkgs", - "type": "github" - } - }, - "nur": { - "locked": { - "lastModified": 1696506445, - "narHash": "sha256-ozu7YxmHsvxSyQazVlkajF8A8U7TaXz3asCL5hFxgNk=", - "owner": "nix-community", - "repo": "NUR", - "rev": "0178289e0bd913fe9847605b01d6e15b7d076f6e", - "type": "github" - }, - "original": { - "owner": "nix-community", - "repo": "NUR", - "type": "github" - } - }, - "nur-xddxdd": { - "inputs": { - "flake-utils": "flake-utils_7", - "flake-utils-plus": "flake-utils-plus", - "nixpkgs": [ - "nixos", - "nixpkgs" - ], - "nvfetcher": "nvfetcher" - }, - "locked": { - "lastModified": 1696487499, - "narHash": "sha256-wvrBwhLpdF+oK5v3Lzgb1Yhz3vT1DHzIL3HKST/tCwU=", - "owner": "xddxdd", - "repo": "nur-packages", - "rev": "9e53a952689cacfd88987c55466450e3076ced05", - "type": "github" - }, - "original": { - "owner": "xddxdd", - "repo": "nur-packages", - "type": "github" - } - }, - "nurl": { - "inputs": { - "nixpkgs": "nixpkgs" - }, - "locked": { - "lastModified": 1692200072, - "narHash": "sha256-s1GdyzB9ACSqQGxJ8Syc0hFn5w8YXhOzDVdNJxQHNpg=", - "owner": "nix-community", - "repo": "nurl", - "rev": "647bc872cbb7927af0ce0a18735749ed1399f85d", - "type": "github" - }, - "original": { - "owner": "nix-community", - "repo": "nurl", - "type": "github" - } - }, - "nvfetcher": { - "inputs": { - "flake-compat": "flake-compat_6", - "flake-utils": [ - "nixos", - "nur-xddxdd", - "flake-utils" - ], - "nixpkgs": [ - "nixos", - "nur-xddxdd", - "nixpkgs" - ] - }, - "locked": { - "lastModified": 1693539235, - "narHash": "sha256-ACmCq1+RnVq+EB7yeN6fThUR3cCJZb6lKEfv937WG84=", - "owner": "berberman", - "repo": "nvfetcher", - "rev": "2bcf73dea96497ac9c36ed320b457caa705f9485", - "type": "github" - }, - "original": { - "owner": "berberman", - "repo": "nvfetcher", - "type": "github" - } - }, - "pnpm2nix-nzbr": { - "inputs": { - "flake-utils": "flake-utils_8", - "nixpkgs": [ - "nixos", - "nixpkgs" - ] - }, - "locked": { - "lastModified": 1692918416, - "narHash": "sha256-SPbFwuIJ5HIhoZxI6snQ3QGyrusqB8eqX/3aOtx+afA=", - "owner": "CHN-beta", - "repo": "pnpm2nix-nzbr", - "rev": "7472a7749fc9f1ef7b75b618d6645f84ead2764d", - "type": "github" - }, - "original": { - "owner": "CHN-beta", - "repo": "pnpm2nix-nzbr", - "type": "github" - } - }, - "qchem": { - "inputs": { - "nixpkgs": [ - "nixos", - "nixpkgs" - ] - }, - "locked": { - "lastModified": 1696260682, - "narHash": "sha256-iccjl57qw6aEe9nsCYFbF2bl7NEI/3Y4cn1U+QYvrFk=", - "owner": "Nix-QChem", - "repo": "NixOS-QChem", - "rev": "7324cb54b7687718ed7b05581998f105fe2fd3e3", - "type": "github" - }, - "original": { - "owner": "Nix-QChem", - "repo": "NixOS-QChem", - "type": "github" - } - }, - "root": { - "inputs": { - "nixos": "nixos", - "nixpkgs": "nixpkgs_5" - } - }, - "sops-nix": { - "inputs": { - "nixpkgs": [ - "nixos", - "nixpkgs" - ], - "nixpkgs-stable": [ - "nixos", - "nixpkgs" - ] - }, - "locked": { - "lastModified": 1696320910, - "narHash": "sha256-fbuEc6wylH+0VxG48lhPBK+SQJHfo2lusUwWHZNipIM=", - "owner": "Mic92", - "repo": "sops-nix", - "rev": "746c7fa1a64c1671a4bf287737c27fdc7101c4c2", - "type": "github" - }, - "original": { - "owner": "Mic92", - "repo": "sops-nix", - "type": "github" - } - }, - "st-src": { - "flake": false, - "locked": { - "lastModified": 1619005886, - "narHash": "sha256-PGD1oVKzaFARUJnn+IhgcDwqx8xuvnouLhUjaGwsMto=", - "owner": "dguibert", - "repo": "st", - "rev": "d7faa1e059e1cdb7ba0d2383c1f0fc81d6d00e15", - "type": "github" - }, - "original": { - "owner": "dguibert", - "ref": "pu", - "repo": "st", - "type": "github" - } - }, - "systems": { - "locked": { - "lastModified": 1681028828, - "narHash": "sha256-Vy1rq5AaRuLzOxct8nz4T6wlgyUR7zLU309k9mBC768=", - "owner": "nix-systems", - "repo": "default", - "rev": "da67096a3b9bf56a91d16901293e51ba5b49a27e", - "type": "github" - }, - "original": { - "owner": "nix-systems", - "repo": "default", - "type": "github" - } - }, - "systems_2": { - "locked": { - "lastModified": 1681028828, - "narHash": "sha256-Vy1rq5AaRuLzOxct8nz4T6wlgyUR7zLU309k9mBC768=", - "owner": "nix-systems", - "repo": "default", - "rev": "da67096a3b9bf56a91d16901293e51ba5b49a27e", - "type": "github" - }, - "original": { - "owner": "nix-systems", - "repo": "default", - "type": "github" - } - }, - "systems_3": { - "locked": { - "lastModified": 1681028828, - "narHash": "sha256-Vy1rq5AaRuLzOxct8nz4T6wlgyUR7zLU309k9mBC768=", - "owner": "nix-systems", - "repo": "default", - "rev": "da67096a3b9bf56a91d16901293e51ba5b49a27e", - "type": "github" - }, - "original": { - "owner": "nix-systems", - "repo": "default", - "type": "github" - } - }, - "systems_4": { - "locked": { - "lastModified": 1681028828, - "narHash": "sha256-Vy1rq5AaRuLzOxct8nz4T6wlgyUR7zLU309k9mBC768=", - "owner": "nix-systems", - "repo": "default", - "rev": "da67096a3b9bf56a91d16901293e51ba5b49a27e", - "type": "github" - }, - "original": { - "owner": "nix-systems", - "repo": "default", - "type": "github" - } - }, - "systems_5": { - "locked": { - "lastModified": 1681028828, - "narHash": "sha256-Vy1rq5AaRuLzOxct8nz4T6wlgyUR7zLU309k9mBC768=", - "owner": "nix-systems", - "repo": "default", - "rev": "da67096a3b9bf56a91d16901293e51ba5b49a27e", - "type": "github" - }, - "original": { - "owner": "nix-systems", - "repo": "default", - "type": "github" - } - }, - "systems_6": { - "locked": { - "lastModified": 1681028828, - "narHash": "sha256-Vy1rq5AaRuLzOxct8nz4T6wlgyUR7zLU309k9mBC768=", - "owner": "nix-systems", - "repo": "default", - "rev": "da67096a3b9bf56a91d16901293e51ba5b49a27e", - "type": "github" - }, - "original": { - "owner": "nix-systems", - "repo": "default", - "type": "github" - } - }, - "systems_7": { - "locked": { - "lastModified": 1681028828, - "narHash": "sha256-Vy1rq5AaRuLzOxct8nz4T6wlgyUR7zLU309k9mBC768=", - "owner": "nix-systems", - "repo": "default", - "rev": "da67096a3b9bf56a91d16901293e51ba5b49a27e", - "type": "github" - }, - "original": { - "owner": "nix-systems", - "repo": "default", - "type": "github" - } - }, - "touchix": { - "inputs": { - "nixpkgs": [ - "nixos", - "nixpkgs" - ] - }, - "locked": { - "lastModified": 1685202181, - "narHash": "sha256-yRj9Vh3T1pL+cNxhiQCFyH67Ys1h7pcrOfSZM10Xb+g=", - "owner": "CHN-beta", - "repo": "touchix", - "rev": "9cedc2f3dc007875525af480976d91b2990847de", - "type": "github" - }, - "original": { - "owner": "CHN-beta", - "repo": "touchix", - "type": "github" - } - }, - "utils": { - "locked": { - "lastModified": 1667395993, - "narHash": "sha256-nuEHfE/LcWyuSWnS8t12N1wc105Qtau+/OdUAjtQ0rA=", - "owner": "numtide", - "repo": "flake-utils", - "rev": "5aed5285a952e0b949eb3ba02c12fa4fcfef535f", - "type": "github" - }, - "original": { - "owner": "numtide", - "repo": "flake-utils", - "type": "github" - } - }, - "utils_2": { - "inputs": { - "systems": "systems_3" - }, - "locked": { - "lastModified": 1689068808, - "narHash": "sha256-6ixXo3wt24N/melDWjq70UuHQLxGV8jZvooRanIHXw0=", - "owner": "numtide", - "repo": "flake-utils", - "rev": "919d646de7be200f3bf08cb76ae1f09402b6f9b4", - "type": "github" - }, - "original": { - "owner": "numtide", - "repo": "flake-utils", - "type": "github" - } - } - }, - "root": "root", - "version": 7 -} diff --git a/flake.nix b/flake.nix deleted file mode 100644 index eef8c04..0000000 --- a/flake.nix +++ /dev/null @@ -1,25 +0,0 @@ -{ - inputs.nixos.url = "github:CHN-beta/nixos"; - inputs.nixpkgs.url = "github:NixOS/nixpkgs/nixos-unstable"; - - outputs = inputs: - let - # pkgs = inputs.nixpkgs.legacyPackages.x86_64-linux; - pkgs = import inputs.nixpkgs - { - localSystem = { system = "x86_64-linux"; gcc = { arch = "alderlake"; tune = "alderlake"; }; }; - config.allowUnfree = true; - }; - localPackages = import "${inputs.nixos}/local/pkgs" { inherit pkgs; inherit (inputs.nixpkgs) lib; }; - in - { - devShell.x86_64-linux = pkgs.mkShell.override { stdenv = pkgs.gcc13Stdenv; } - { - packages = with pkgs; [ pkg-config cmake ninja ]; - buildInputs = (with pkgs; [ eigen yaml-cpp fmt highfive tbb_2021_8.dev glfw libGL ]) - ++ (with localPackages; [ concurrencpp matplotplusplus zpp-bits glad ]); - hardeningDisable = [ "all" ]; - # NIX_DEBUG = "1"; - }; - }; -} diff --git a/include/ufo.hpp b/include/ufo.hpp new file mode 100644 index 0000000..4e48475 --- /dev/null +++ b/include/ufo.hpp @@ -0,0 +1,74 @@ +# pragma once +# include + +namespace ufo +{ + // 在相位中, 约定为使用 $\exp (2 \pi i \vec{q} \cdot \vec{r})$ 来表示原子的运动状态 + // (而不是 $\exp (-2 \pi i \vec{q} \cdot \vec{r})$) + // 一些书定义的倒格矢中包含了 $2 \pi$ 的部分, 我们这里约定不包含这部分. + // 也就是说, 正格子与倒格子的转置相乘, 得到单位矩阵. + + using namespace biu::literals; + using namespace biu::stream_operators; + + void fold(std::string config_file); + void unfold(std::string config_file); + void plot_band(std::string config_file); + void plot_point(std::string config_file); + void raman_create_displacement(std::string config_file); + void raman_apply_contribution(std::string config_file); + + // 许多函数都需要用到这个，所以写到头文件中 + struct UnfoldOutput + { + Eigen::Matrix3d PrimativeCell; + Eigen::Matrix3i SuperCellTransformation; + Eigen::Vector3i SuperCellMultiplier; + Eigen::Matrix3d SuperCellDeformation; + std::optional> 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>; + }; +} diff --git a/include/ufo/fold.hpp b/include/ufo/fold.hpp deleted file mode 100644 index 6ce66dc..0000000 --- a/include/ufo/fold.hpp +++ /dev/null @@ -1,37 +0,0 @@ -# pragma once -# include - -namespace ufo -{ - class FoldSolver : public Solver - { - public: - struct InputType - { - Eigen::Vector SuperCellMultiplier; - std::optional> SuperCellDeformation; - std::vector Qpoints; - DataFile OutputFile; - - InputType(std::string config_file); - }; - struct OutputType - { - std::vector Qpoints; - void write(std::string filename) const; - }; - protected: - InputType Input_; - std::optional Output_; - public: - FoldSolver(std::string config_file); - FoldSolver& operator()() override; - // return value: QpointInReciprocalSuperCellByReciprocalSuperCell - static Eigen::Vector3d fold - ( - Eigen::Vector3d qpoint_in_reciprocal_primitive_cell_by_reciprocal_primitive_cell, - Eigen::Vector super_cell_multiplier, - std::optional> super_cell_deformation - ); - }; -} \ No newline at end of file diff --git a/include/ufo/plot.hpp b/include/ufo/plot.hpp deleted file mode 100644 index 8c203ee..0000000 --- a/include/ufo/plot.hpp +++ /dev/null @@ -1,77 +0,0 @@ -# pragma once -# include - -namespace ufo -{ - class PlotSolver : public Solver - { - public: - struct InputType - { - Eigen::Matrix3d PrimativeCell; - - struct FigureConfigType - { - std::vector> Qpoints; - std::pair Resolution; - std::pair Range; - std::optional> YTicks; - DataFile PictureFile; - std::optional> DataFiles; - }; - std::vector 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> Values; - std::vector XTicks; - std::vector YTicks; - std::pair Resolution; - std::pair Range; - - OutputType() = default; - const OutputType& write(std::string filename, std::string format) const; - using serialize = zpp::bits::members<5>; - }; - protected: - InputType Input_; - std::optional> Output_; - public: - PlotSolver(std::string config_file); - PlotSolver& operator()() override; - - // 根据 q 点路径, 搜索要使用的 q 点 - static std::vector> search_qpoints - ( - const std::pair& path, - const decltype(InputType::UnfoldedDataType::QpointData)& available_qpoints, - double threshold, bool exclude_endpoint = false - ); - // 根据搜索到的 q 点, 计算每个点的数值 - static std::tuple>, std::vector> calculate_values - ( - const Eigen::Matrix3d primative_cell, - const std::vector>& path, - const std::vector>>& qpoints, - const decltype(InputType::FigureConfigType::Resolution)& resolution, - const decltype(InputType::FigureConfigType::Range)& range - ); - // 根据数值, 画图 - static void plot - ( - const std::vector>& values, - const std::string& filename, - const std::vector& x_ticks, const std::vector& y_ticks - ); - }; -} diff --git a/include/ufo/solver.hpp b/include/ufo/solver.hpp deleted file mode 100644 index e1e2c36..0000000 --- a/include/ufo/solver.hpp +++ /dev/null @@ -1,141 +0,0 @@ -# pragma once -# include -# include -# include -# include -# include -# include -# include -# include -# include -# include -# include -# include -# include -# include -# include -# include -# include -# include -# include -# include -# include -# include -# include -# include - -// 在相位中, 约定为使用 $\exp (2 \pi i \vec{q} \cdot \vec{r})$ 来表示原子的运动状态 -// (而不是 $\exp (-2 \pi i \vec{q} \cdot \vec{r})$) -// 一些书定义的倒格矢中包含了 $2 \pi$ 的部分, 我们这里约定不包含这部分. -// 也就是说, 正格子与倒格子的转置相乘, 得到单位矩阵. - -namespace Eigen -{ - constexpr inline auto serialize(auto & archive, Eigen::Matrix3d& matrix) - { return archive(std::span(matrix.data(), matrix.size())); } - constexpr inline auto serialize(auto & archive, const Eigen::Matrix3d& matrix) - { return archive(std::span(matrix.data(), matrix.size())); } - constexpr inline auto serialize(auto & archive, Eigen::Vector3d& vector) - { return archive(std::span(vector.data(), vector.size())); } - constexpr inline auto serialize(auto & archive, const Eigen::Vector3d& vector) - { return archive(std::span(vector.data(), vector.size())); } -} - -namespace ufo -{ - using namespace std::literals; - struct PhonopyComplex { double r, i; }; - inline HighFive::CompoundType create_compound_complex() - { return {{ "r", HighFive::AtomicType{}}, {"i", HighFive::AtomicType{}}}; } - - namespace detail_ - { - template struct SpecializationOfBitsMembersHelper : std::false_type {}; - template struct SpecializationOfBitsMembersHelper> : std::true_type {}; - } - template concept ZppSerializable - = requires() { detail_::SpecializationOfBitsMembersHelper::value == true; }; - - class Solver - { - public: - virtual Solver& operator()() = 0; - virtual ~Solver() = default; - - static concurrencpp::generator, unsigned>> - triplet_sequence(Eigen::Vector range); - - template inline static void zpp_write(const T& object, std::string filename) - { - auto [data, out] = zpp::bits::data_out(); - out(object).or_throw(); - static_assert(sizeof(char) == sizeof(std::byte)); - std::ofstream file(filename, std::ios::binary | std::ios::out); - file.exceptions(std::ios::badbit | std::ios::failbit); - file.write(reinterpret_cast(data.data()), data.size()); - } - template 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 data; - { - std::vector string(std::istreambuf_iterator(input), {}); - data.assign - ( - reinterpret_cast(string.data()), - reinterpret_cast(string.data() + string.size()) - ); - } - auto in = zpp::bits::in(data); - T output; - in(output).or_throw(); - return output; - } - - class Hdf5file - { - public: - inline Hdf5file& open_for_read(std::string filename) - { - File_ = HighFive::File(filename, HighFive::File::ReadOnly); - return *this; - } - inline Hdf5file& open_for_write(std::string filename) - { - File_ = HighFive::File(filename, HighFive::File::ReadWrite | HighFive::File::Create - | HighFive::File::Truncate); - return *this; - } - template inline Hdf5file& read(T& object, std::string name) - { - object = File_->getDataSet(name).read>(); - return *this; - } - template inline Hdf5file& write(const T& object, std::string name) - { - File_->createDataSet(name, object); - return *this; - } - protected: - std::optional File_; - }; - - struct DataFile - { - std::string Filename; - std::string Format; - std::map ExtraParameters; - inline DataFile() = default; - DataFile - ( - YAML::Node node, std::set supported_format, - std::string config_file, bool allow_same_as_config_file = false - ); - }; - - }; -} - -HIGHFIVE_REGISTER_TYPE(ufo::PhonopyComplex, ufo::create_compound_complex) diff --git a/include/ufo/unfold.hpp b/include/ufo/unfold.hpp deleted file mode 100644 index c53dc00..0000000 --- a/include/ufo/unfold.hpp +++ /dev/null @@ -1,164 +0,0 @@ -# pragma once -# include - -namespace ufo -{ - // 反折叠的原理: 将超胞中的原子运动状态, 投影到一组平面波构成的基矢中. - // 每一个平面波的波矢由两部分相加得到: 一部分是单胞倒格子的整数倍, 所取的个数有一定任意性, 论文中建议取大约单胞中原子个数那么多个; - // 对于没有缺陷的情况, 取一个应该就足够了. - // 另一部分是超胞倒格子的整数倍, 取 n 个, n 为超胞对应的单胞的倍数, 其实也就是倒空间中单胞对应倒格子中超胞的格点. - // 只要第一部分取得足够多, 那么单胞中原子的状态就可以完全被这些平面波描述. - // 将超胞中原子的运动状态投影到这些基矢上, 计算出投影的系数, 就可以将超胞的原子运动状态分解到单胞中的多个 q 点上. - class UnfoldSolver : public Solver - { - public: - struct InputType - { - // 单胞的三个格矢，每行表示一个格矢的坐标，单位为埃 - Eigen::Matrix3d PrimativeCell; - // 单胞到超胞的格矢转换时用到的矩阵 - // SuperCellMultiplier 是一个三维列向量且各个元素都是整数，表示单胞在各个方向扩大到多少倍之后，可以得到和超胞一样的体积 - // SuperCsolver.hpp>mation 是一个行列式为 1 的矩阵，它表示经过 SuperCellMultiplier 扩大后，还需要怎样的变换才能得到超胞 - // SuperCell = (SuperCellDeformation * SuperCellMultiplier.asDiagonal()) * PrimativeCell - // ReciprocalPrimativeCell = (SuperCellDeformation * SuperCellMultiplier.asDiagonal()).transpose() - // * ReciprocalSuperCell - // Position = PositionToCell(line vector) * Cell - // InversePosition = InversePositionToCell(line vector) * ReciprocalCell - // PositionToSuperCell(line vector) * SuperCell = PositionToPrimativeCell(line vector) * PrimativeCell - // ReciprocalPositionToSuperCell(line vector) * ReciprocalSuperCell - // = ReciprocalPositionToPrimativeCell(line vector) * ReciprocalPrimativeCell - Eigen::Vector SuperCellMultiplier; - std::optional> SuperCellDeformation; - // 在单胞内取几个平面波的基矢 - Eigen::Vector 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 ModeData; - }; - std::vector QpointData; - - // 输出到哪些文件, 以及使用怎样的格式, 格式可选值包括: - // yaml: 使用 yaml 格式输出 - // yaml-human-readable: 使用 yaml 格式输出, 但是输出的结果更适合人类阅读, - // 包括合并相近的模式, 去除权重过小的模式, 限制输出的小数位数. - // zpp: 使用 zpp-bits 序列化, 可以直接被 plot.cpp 读取 - std::vector 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 ModeData; - }; - std::vector 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>; - protected: - InputType Input_; - std::optional Output_; - std::optional Basis_; - - // 第一层是不同的模式, 第二层是不同的 sub qpoint - using ProjectionCoefficientType_ = std::vector>; - - public: - UnfoldSolver(std::string config_file); - UnfoldSolver& operator()() override; - - // 构建基 - // 每个 q 点对应的一组 sub qpoint。不同的 q 点所对应的 sub qpoint 是不一样的，但 sub qpoint 与 q 点的相对位置一致。 - // 这里 xyz_of_diff_of_sub_qpoint 即表示这个相对位置。 - // 由于基只与这个相对位置有关（也就是说，不同 q 点的基是一样的），因此可以先计算出所有的基，这样降低计算量。 - // 外层下标对应超胞倒格子的整数倍那部分(第二部分), 也就是不同的 sub qpoint - // 内层下标对应单胞倒格子的整数倍那部分(第一部分), 也就是 sub qpoint 上的不同平面波（取的数量越多，结果越精确） - static BasisType construct_basis - ( - const decltype(InputType::PrimativeCell)& primative_cell, - const decltype(InputType::SuperCellMultiplier)& super_cell_multiplier, - const decltype(InputType::PrimativeCellBasisNumber)& - primative_cell_basis_number, - const decltype(InputType::AtomPosition)& atom_position - ); - - // 计算投影系数, 是反折叠的核心步骤 - ProjectionCoefficientType_ construct_projection_coefficient - ( - const BasisType& basis, - const std::vector>& mode_data, - std::atomic& number_of_finished_modes - ); - - OutputType construct_output - ( - const decltype(InputType::SuperCellMultiplier)& super_cell_multiplier, - const decltype(InputType::SuperCellDeformation)& super_cell_deformation, - const std::vector>& meta_qpoint_by_reciprocal_super_cell, - const std::vector>>& frequency, - const ProjectionCoefficientType_& projection_coefficient - ); - }; -} diff --git a/src/fold.cpp b/src/fold.cpp index 9732be9..bc189bc 100644 --- a/src/fold.cpp +++ b/src/fold.cpp @@ -1,86 +1,58 @@ -# include +# include -namespace ufo +void ufo::fold(std::string config_file) { - FoldSolver::InputType::InputType(std::string config_file) + struct Input { - auto input = YAML::LoadFile(config_file); - for (unsigned i = 0; i < 3; i++) - SuperCellMultiplier(i) = input["SuperCellMultiplier"][i].as(); - if (input["SuperCellDeformation"]) - { - SuperCellDeformation.emplace(); - for (unsigned i = 0; i < 3; i++) - for (unsigned j = 0; j < 3; j++) - (*SuperCellDeformation)(i, j) = input["SuperCellDeformation"][i][j].as(); - } - for (auto& qpoint : input["Qpoints"].as>>()) - Qpoints.push_back(Eigen::Vector3d - {{qpoint.at(0)}, {qpoint.at(1)}, {qpoint.at(2)}}); - OutputFile = DataFile(input["OutputFile"], {"yaml"}, config_file); - } - void FoldSolver::OutputType::write(std::string filename) const + Eigen::Matrix3d SuperCellDeformation; + Eigen::Vector3i SuperCellMultiplier; + std::vector Qpoints; + std::optional OutputFile; + }; + struct Output { - std::ofstream(filename) << [&] - { - std::stringstream print; - print << "Qpoints:\n"; - for (auto& qpoint : Qpoints) - print << fmt::format(" - [ {:.8f}, {:.8f}, {:.8f} ]\n", qpoint(0), qpoint(1), qpoint(2)); - return print.str(); - }(); - } - - FoldSolver::FoldSolver(std::string config_file) : Input_(config_file) {} - FoldSolver& FoldSolver::operator()() - { - if (!Output_) - { - Output_.emplace(); - for (auto& qpoint : Input_.Qpoints) - Output_->Qpoints.push_back(fold - ( - qpoint, Input_.SuperCellMultiplier, - Input_.SuperCellDeformation - )); - } - Output_->write(Input_.OutputFile.Filename); - return *this; - } - - Eigen::Vector3d FoldSolver::fold + std::vector Qpoints; + }; + auto fold = [] ( Eigen::Vector3d qpoint_in_reciprocal_primitive_cell_by_reciprocal_primitive_cell, - Eigen::Vector super_cell_multiplier, - std::optional> super_cell_deformation - ) + Eigen::Matrix3d super_cell_transformation + ) -> Eigen::Vector3d { - // 首先需要将 q 点转移到 ModifiedSuperCell 的倒格子中 - // 将 q 点坐标扩大, 然后取小数部分, 就可以了 - auto qpoint_by_reciprocal_modified_super_cell = super_cell_multiplier.cast().asDiagonal() - * qpoint_in_reciprocal_primitive_cell_by_reciprocal_primitive_cell; - auto qpoint_in_reciprocal_modified_super_cell_by_reciprocal_modified_super_cell = - (qpoint_by_reciprocal_modified_super_cell.array() - qpoint_by_reciprocal_modified_super_cell.array().floor()) - .matrix(); - if (!super_cell_deformation) - return qpoint_in_reciprocal_modified_super_cell_by_reciprocal_modified_super_cell; /* - 对 q 点平移数个 SupreCell, 直到它落在超胞的倒格子中 - 这等价于直接将 q 点坐标用 SuperCell 的倒格子表示, 然后取小数部分. - ModifiedSuperCell = SuperCellMultiplier * PrimativeCell - SuperCell = SuperCellDeformation * ModifiedSuperCell - ReciprocalModifiedSuperCell = ModifiedSuperCell.inverse().transpose() - ReciprocalSuperCell = SuperCell.inverse().transpose() - Qpoint = QpointByReciprocalModifiedSuperCell.transpose() * ReciprocalModifiedSuperCell - Qpoint = QpointByReciprocalSuperCell.transpose() * ReciprocalSuperCell + 首先需要将 q 点坐标的单位转换为 ModifiedSuperCell 的格矢，可知： + QpointByReciprocalModifiedSuperCell = SuperCellMultiplier * QpointByReciprocalPrimitiveCell; + 接下来考虑将 q 点坐标的单位转换为 SuperCell 的格矢 + ModifiedSuperCell = SuperCellMultiplier * PrimativeCell; + SuperCell = SuperCellDeformation * ModifiedSuperCell; + ReciprocalModifiedSuperCell = ModifiedSuperCell.inverse().transpose(); + ReciprocalSuperCell = SuperCell.inverse().transpose(); + Qpoint = QpointByReciprocalModifiedSuperCell.transpose() * ReciprocalModifiedSuperCell; + Qpoint = QpointByReciprocalSuperCell.transpose() * ReciprocalSuperCell; 整理可以得到: - QpointByReciprocalSuperCell = SuperCellDeformation * QpointByReciprocalModifiedSuperCell + QpointByReciprocalSuperCell = SuperCellDeformation * QpointByReciprocalModifiedSuperCell; + 两个式子结合，可以得到： + QpointByReciprocalSuperCell = SuperCellDeformation * SuperCellMultiplier * QpointByReciprocalPrimitiveCell; */ - auto qpoint_in_reciprocal_modified_super_cell_by_reciprocal_super_cell = - (*super_cell_deformation * qpoint_in_reciprocal_modified_super_cell_by_reciprocal_modified_super_cell).eval(); - auto qpoint_in_reciprocal_super_cell_by_reciprocal_super_cell = - qpoint_in_reciprocal_modified_super_cell_by_reciprocal_super_cell.array() - - qpoint_in_reciprocal_modified_super_cell_by_reciprocal_super_cell.array().floor(); - return qpoint_in_reciprocal_super_cell_by_reciprocal_super_cell; - } + auto qpoint_by_reciprocal_super_cell = + ( + super_cell_transformation * qpoint_in_reciprocal_primitive_cell_by_reciprocal_primitive_cell + ).eval(); + /* + 到目前为止，我们还没有移动过 q 点的坐标。现在，我们将它移动整数个 ReciprocalSuperCell，直到它落在超胞的倒格子中。 + 这等价于直接取 QpointByReciprocalSuperCell - QpointByReciprocalSuperCell.floor()。 + */ + return (qpoint_by_reciprocal_super_cell.array() - qpoint_by_reciprocal_super_cell.array().floor()).matrix(); + }; + auto input = YAML::LoadFile(config_file).as(); + 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); } diff --git a/src/main.cpp b/src/main.cpp index 8bed4aa..1291b32 100644 --- a/src/main.cpp +++ b/src/main.cpp @@ -1,17 +1,14 @@ -# include -# include -# include +# include int main(int argc, const char** argv) { - if (argc != 3) - throw std::runtime_error(fmt::format("Usage: {} task config.yaml", argv[0])); - if (argv[1] == std::string("fold")) - ufo::FoldSolver{argv[2]}(); - else if (argv[1] == std::string("unfold")) - ufo::UnfoldSolver{argv[2]}(); - else if (argv[1] == std::string("plot")) - ufo::PlotSolver{argv[2]}(); - else - throw std::runtime_error(fmt::format("Unknown task: {}", argv[1])); + using namespace biu::literals; + if (argc != 3) throw std::runtime_error("Usage: {} task config.yaml"_f(argv[0])); + if (argv[1] == "fold"s) ufo::fold(argv[2]); + else if (argv[1] == "unfold"s) ufo::unfold(argv[2]); + else if (argv[1] == "raman-create-displacement"s) ufo::raman_create_displacement(argv[2]); + else if (argv[1] == "raman-apply-contribution"s); + else if (argv[1] == "plot-band"s) ufo::plot_band(argv[2]); + else if (argv[1] == "plot-point"s) ufo::plot_point(argv[2]); + else throw std::runtime_error("Unknown task: {}"_f(argv[1])); } diff --git a/src/plot.cpp b/src/plot.cpp index f607807..c756cdc 100644 --- a/src/plot.cpp +++ b/src/plot.cpp @@ -1,233 +1,177 @@ -# include +# include +# include +# include -namespace ufo +void ufo::plot_band(std::string config_file) { - PlotSolver::InputType::UnfoldedDataType::UnfoldedDataType(std::string filename) + struct Input { - static_cast(*this) = zpp_read(filename); - } + std::string UnfoldedDataFile; + // 要画图的 q 点路径列表 + // 内层表示一个路径上的 q 点，外层表示不同的路径 + // 单位为倒格矢 + std::vector> Qpoints; + // 插值时使用的分辨率（不影响画出来图片的分辨率和横纵比） + std::array InterpolationResolution; + // 画图区域的y轴和x轴的比例。如果不指定，则由matplot++自动调整（通常调整为正方形，即 1） + std::optional AspectRatio; + // 整张图片的分辨率 + std::optional> PictureResolution; + // 画图的频率范围 + std::array FrequencyRange; + // 搜索 q 点时的阈值，单位为埃^-1 + std::optional ThresholdWhenSearchingQpoints; + // 是否要在 y 轴上作一些标记 + std::optional>> YTicks; + // 是否输出图片 + std::optional OutputPictureFile; + // 是否输出数据，可以进一步使用 matplotlib 画图 + std::optional OutputDataFile; + }; - PlotSolver::InputType::InputType(std::string config_file) - { - auto input = YAML::LoadFile(config_file); - for (unsigned i = 0; i < 3; i++) - for (unsigned j = 0; j < 3; j++) - PrimativeCell(i, j) = input["PrimativeCell"][i][j].as(); - for (auto& figure : input["Figures"].as>()) - { - Figures.emplace_back(); - auto qpoints = figure["Qpoints"] - .as>>>(); - 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>(); - Figures.back().Range = figure["Range"].as>(); - Figures.back().PictureFile - = DataFile(figure["PictureFile"], {"png"}, config_file); - if (figure["YTicks"]) - Figures.back().YTicks = figure["YTicks"].as>(); - if (figure["DataFiles"]) - { - Figures.back().DataFiles.emplace(); - for (auto& data_file : figure["DataFiles"].as>()) - Figures.back().DataFiles->emplace_back() - = DataFile(data_file, {"hdf5", "zpp"}, config_file); - } - } - UnfoldedDataFile = DataFile(input["UnfoldedDataFile"], {"zpp"}, config_file); - UnfoldedData = UnfoldedDataType(UnfoldedDataFile.Filename); - } - const PlotSolver::OutputType& PlotSolver::OutputType::write(std::string filename, std::string format) const - { - if (format == "zpp") - zpp_write(*this, filename); - else if (format == "hdf5") - { - std::vector resolution{ Resolution.first, Resolution.second }; - std::vector range{ Range.first, Range.second }; - Hdf5file{}.open_for_write(filename).write(Values, "Values") - .write(XTicks, "XTicks") - .write(YTicks, "YTicks") - .write(resolution, "Resolution") - .write(range, "Range"); - } - return *this; - } - - PlotSolver::PlotSolver(std::string config_file) : Input_(config_file) {} - - PlotSolver& PlotSolver::operator()() - { - Output_.emplace(); - for (auto& figure : Input_.Figures) - { - // 外层表示不同的线段的端点，内层表示这个线段上的 q 点 - std::vector>> qpoints; - std::vector> 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{}); - 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> PlotSolver::search_qpoints + // 根据 q 点路径, 搜索要使用的 q 点，返回的是 q 点在 QpointData 中的索引以及到路径起点的距离，以及这段路径的总长度 + auto search_qpoints = [] ( + const Eigen::Matrix3d& primative_cell, const std::pair& path, - const decltype(InputType::UnfoldedDataType::QpointData)& available_qpoints, - double threshold, bool exclude_endpoint + const std::vector& qpoints, + double threshold, bool exclude_endpoint = false ) { - std::multimap> 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 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()).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> 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> PlotSolver::calculate_values + // 根据搜索到的 q 点, 计算图中每个点的值 + auto calculate_values = [] ( - const Eigen::Matrix3d primative_cell, - const std::vector>& path, - const std::vector>>& qpoints, - const decltype(InputType::FigureConfigType::Resolution)& resolution, - const decltype(InputType::FigureConfigType::Range)& range + // search_qpoints 的第一个返回值 + const boost::container::flat_map& path, + // 每一条连续路径的第一个 q 点的索引 + const std::set& path_begin, + // 所有 q 点的数据（需要用到它的频率和权重） + const std::vector& qpoints, + // 用于插值的分辨率和范围 + const std::array& resolution, + const std::array& frequency_range, + // 路径的总长度 + double total_distance ) { - // 整理输入 - std::map> qpoints_with_distance; - double total_distance = 0; - std::vector 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> values; - auto blend = [] + // 按比例混合两个 q 点的结果，得到可以用于画图的那一列数据 + auto blend = [&] ( - const UnfoldSolver::OutputType::QpointDataType& a, - const UnfoldSolver::OutputType::QpointDataType& b, - double ratio, unsigned resolution, std::pair range + // 两个点的索引 + std::size_t a, std::size_t b, + // 按照连续路径混合还是按照断开的路径混合 + bool continuous, + // 第一个点占的比例 + double ratio, + std::size_t resolution, std::array frequency_range ) -> std::vector { - // 计算插值结果 + // 混合得到的频率和权重 std::vector frequency, weight; - for (unsigned i = 0; i < a.ModeData.size(); i++) + // 如果是连续路径，将每个模式的频率和权重按照比例混合 + if (continuous) { - frequency.push_back(a.ModeData[i].Frequency * ratio + b.ModeData[i].Frequency * (1 - ratio)); - weight.push_back(a.ModeData[i].Weight * ratio + b.ModeData[i].Weight * (1 - ratio)); + assert(qpoints[a].ModeData.size() == qpoints[b].ModeData.size()); + for (std::size_t i = 0; i < qpoints[a].ModeData.size(); i++) + { + frequency.push_back + (qpoints[a].ModeData[i].Frequency * ratio + qpoints[b].ModeData[i].Frequency * (1 - ratio)); + weight.push_back(qpoints[a].ModeData[i].Weight * ratio + qpoints[b].ModeData[i].Weight * (1 - ratio)); + } + } + // 如果是不连续路径，将每个模式的权重乘以比例，最后相加 + else + { + for (std::size_t i = 0; i < qpoints[a].ModeData.size(); i++) + { + frequency.push_back(qpoints[a].ModeData[i].Frequency); + weight.push_back(qpoints[a].ModeData[i].Weight * ratio); + } + for (std::size_t i = 0; i < qpoints[b].ModeData.size(); i++) + { + frequency.push_back(qpoints[b].ModeData[i].Frequency); + weight.push_back(qpoints[b].ModeData[i].Weight * (1 - ratio)); + } } std::vector result(resolution); - for (unsigned i = 0; i < frequency.size(); i++) + for (std::size_t i = 0; i < frequency.size(); i++) { - int index = (frequency[i] - range.first) / (range.second - range.first) * resolution; - if (index >= 0 && index < static_cast(resolution)) - result[index] += weight[i]; + std::ptrdiff_t index = (frequency[i] - frequency_range[0]) / (frequency_range[1] - frequency_range[0]) + * resolution; + if (index >= 0 && index < static_cast(resolution)) result[index] += weight[i]; } return result; }; - for (unsigned i = 0; i < resolution.first; i++) + + std::vector> values; + for (std::size_t i = 0; i < resolution[0]; i++) { - auto current_distance = total_distance * i / resolution.first; - auto it = qpoints_with_distance.lower_bound(current_distance); - if (it == qpoints_with_distance.begin()) - values.push_back(blend(it->second.get(), it->second.get(), 1, resolution.second, range)); - else if (it == qpoints_with_distance.end()) - values.push_back(blend(std::prev(it)->second.get(), std::prev(it)->second.get(), 1, resolution.second, - range)); - else - values.push_back(blend - ( - std::prev(it)->second.get(), it->second.get(), - (it->first - current_distance) / (it->first - std::prev(it)->first), - resolution.second, range) - ); + auto current_distance = total_distance * i / resolution[0]; + auto it = path.lower_bound(current_distance); + if (it == path.begin()) values.push_back(blend + (it->second, it->second, true, 1, resolution[1], frequency_range)); + else if (it == path.end()) values.push_back(blend + ( + std::prev(it)->second, std::prev(it)->second, true, 1, + resolution[1], frequency_range + )); + else values.push_back(blend + ( + std::prev(it)->second, it->second, !path_begin.contains(it->second), + (it->first - current_distance) / (it->first - std::prev(it)->first), + resolution[1], frequency_range + )); } - return {values, x_ticks}; - } - void PlotSolver::plot + return values; + }; + + // 根据数值, 画图 + auto plot = [] ( const std::vector>& values, const std::string& filename, - const std::vector& x_ticks, const std::vector& y_ticks + const std::vector& x_ticks, const std::vector& y_ticks, + const std::vector& y_ticklabels, + const std::optional& aspect_ratio, + const std::optional>& resolution ) { std::vector> @@ -235,32 +179,253 @@ namespace ufo g(values[0].size(), std::vector(values.size(), 0)), b(values[0].size(), std::vector(values.size(), 0)), a(values[0].size(), std::vector(values.size(), 0)); - for (unsigned i = 0; i < values[0].size(); i++) - for (unsigned j = 0; j < values.size(); j++) + for (std::size_t i = 0; i < values[0].size(); i++) + for (std::size_t j = 0; j < values.size(); j++) { auto v = values[j][i]; - if (v < 0.05) - v = 0; + if (v < 0.05) v = 0; a[i][j] = v * 100 * 255; - if (a[i][j] > 255) - a[i][j] = 255; + if (a[i][j] > 255) a[i][j] = 255; r[i][j] = 255 - v * 2 * 255; - if (r[i][j] < 0) - r[i][j] = 0; + if (r[i][j] < 0) r[i][j] = 0; g[i][j] = 255 - v * 2 * 255; - if (g[i][j] < 0) - g[i][j] = 0; + if (g[i][j] < 0) g[i][j] = 0; b[i][j] = 255; } - auto f = matplot::figure(true); + auto f = matplot::figure(true); auto ax = f->current_axes(); auto image = ax->image(std::tie(r, g, b)); image->matrix_a(a); ax->y_axis().reverse(false); ax->x_axis().tick_values(x_ticks); ax->x_axis().tick_length(1); + ax->x_axis().ticklabels(std::vector(x_ticks.size())); ax->y_axis().tick_values(y_ticks); ax->y_axis().tick_length(1); + ax->y_axis().ticklabels(y_ticklabels); + if (aspect_ratio) + { + ax->axes_aspect_ratio_auto(false); + ax->axes_aspect_ratio(*aspect_ratio); + } + if (resolution) + { + f->width((*resolution)[0]); + f->height((*resolution)[1]); + } f->save(filename, "png"); + }; + + auto input = YAML::LoadFile(config_file).as(); + auto unfolded_data = biu::deserialize + (biu::read(input.UnfoldedDataFile)); + + // 搜索画图需要用到的 q 点 + // key 到起点的距离，value 为 q 点在 QpointData 中的索引 + 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) + { + 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 + ); + total_distance += this_distance; + } } + + // 计算画图的数据 + auto values = calculate_values + ( + path, path_begin, unfolded_data.QpointData, input.InterpolationResolution, + input.FrequencyRange, total_distance + ); + auto x_ticks = x_ticks_index | ranges::views::transform([&](auto i) + { return path.nth(i)->first / total_distance * input.InterpolationResolution[0]; }) | ranges::to; + auto y_ticks = input.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]; + }) + | ranges::to_vector; + auto y_ticklabels = input.YTicks.value_or(std::vector>{}) + | biu::toLvalue | ranges::views::values | ranges::to_vector; + if (input.OutputPictureFile) plot + ( + values, input.OutputPictureFile.value(), + x_ticks, y_ticks, y_ticklabels, input.AspectRatio, input.PictureResolution + ); + if (input.OutputDataFile) + biu::Hdf5file(input.OutputDataFile.value(), true) + .write("Values", values) + .write("XTicks", x_ticks) + .write("YTicks", y_ticks) + .write("YTickLabels", y_ticklabels) + .write("InterpolationResolution", input.InterpolationResolution) + .write("FrequencyRange", input.FrequencyRange); +} + +void ufo::plot_point(std::string config_file) +{ + struct Input + { + std::string UnfoldedDataFile; + // 要画图的 q 点 + Eigen::Vector3d Qpoint; + // 插值的分辨率 + std::size_t InterpolationResolution; + std::optional AspectRatio; + std::optional> PictureResolution; + // 画图的频率范围 + std::array FrequencyRange; + // 搜索 q 点时的阈值，单位为埃^-1 + std::optional ThresholdWhenSearchingQpoints; + // 是否要在 z 轴上作一些标记 + std::optional>> XTicks; + // 是否输出图片 + std::optional OutputPictureFile; + // 是否输出插值后数据，可以进一步使用 matplotlib 画图 + std::optional OutputDataFile; + // 是否输出插值前数据，可以配合 phonopy 结果深入研究 + std::optional OutputRawDataFile; + }; + + // 根据 q 点路径, 搜索要使用的 q 点，返回的是 q 点在 QpointData 中的索引 + auto search_qpoints = [] + ( + const Eigen::Matrix3d& primative_cell, + 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 (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); + } + throw std::runtime_error("No q points found"); + }; + + // 根据搜索到的 q 点, 计算图中每个点的值 + auto calculate_values = [] + ( + // q 点的数据（需要用到它的频率和权重） + const UnfoldOutput::QpointDataType& qpoint, + // 用于插值的分辨率和范围 + std::size_t resolution, + const std::array& frequency_range + ) + { + biu::Logger::Guard log; + std::vector result(resolution); + for (auto& mode : qpoint.ModeData) + { + double index_double = (mode.Frequency - frequency_range[0]) / (frequency_range[1] - frequency_range[0]) + * (resolution - 1); + std::ptrdiff_t index = std::round(index_double); + if (index >= 0 && index < static_cast(resolution)) result[index] += mode.Weight; + } + return log.rtn(result); + }; + + // 根据数值, 画图 + auto plot = [] + ( + const std::vector& values, const std::string& filename, + const std::vector& x_ticks, const std::vector& x_ticklabels, + const std::optional& aspect_ratio, const std::optional>& resolution + ) + { + biu::Logger::Guard log; + auto f = matplot::figure(true); + auto ax = f->current_axes(); + auto image = ax->area(values, 0, false, ""); + ax->y_axis().reverse(false); + ax->x_axis().tick_values(x_ticks); + ax->x_axis().tick_length(1); + ax->x_axis().ticklabels(x_ticklabels); + ax->y_axis().tick_values({}); + if (aspect_ratio) + { + ax->axes_aspect_ratio_auto(false); + ax->axes_aspect_ratio(*aspect_ratio); + } + if (resolution) + { + f->width((*resolution)[0]); + f->height((*resolution)[1]); + } + f->save(filename, "png"); + }; + + biu::Logger::Guard log; + auto input = YAML::LoadFile(config_file).as(); + auto unfolded_data = biu::deserialize + (biu::read(input.UnfoldedDataFile)); + + auto qpoint_index = search_qpoints + ( + unfolded_data.PrimativeCell, input.Qpoint, + unfolded_data.QpointData + | ranges::views::transform(&UnfoldOutput::QpointDataType::Qpoint) + | ranges::to_vector, + input.ThresholdWhenSearchingQpoints.value_or(0.001) + ); + auto values = calculate_values + ( + unfolded_data.QpointData[qpoint_index], + input.InterpolationResolution, input.FrequencyRange + ); + auto x_ticks = input.XTicks.value_or(std::vector>{}) + | biu::toLvalue | ranges::views::keys + | ranges::views::transform([&](auto i) + { + return (i - input.FrequencyRange[0]) / (input.FrequencyRange[1] - input.FrequencyRange[0]) + * input.InterpolationResolution; + }) + | ranges::to_vector; + auto x_ticklabels = input.XTicks.value_or(std::vector>{}) + | biu::toLvalue | ranges::views::values | ranges::to_vector; + if (input.OutputPictureFile) plot + ( + values, input.OutputPictureFile.value(), + x_ticks, x_ticklabels, input.AspectRatio, input.PictureResolution + ); + if (input.OutputDataFile) + biu::Hdf5file(input.OutputDataFile.value(), true) + .write("Values", values) + .write("XTicks", x_ticks) + .write("XTickLabels", x_ticklabels) + .write("InterpolationResolution", input.InterpolationResolution) + .write("FrequencyRange", input.FrequencyRange); + if (input.OutputRawDataFile) + std::ofstream(*input.OutputRawDataFile) << YAML::Node(unfolded_data.QpointData[qpoint_index]); } diff --git a/src/raman-create-displacement.cpp b/src/raman-create-displacement.cpp new file mode 100644 index 0000000..226bdf8 --- /dev/null +++ b/src/raman-create-displacement.cpp @@ -0,0 +1,138 @@ +# 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; + }; + 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 = [] + ( + 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(); + }; + + auto input = YAML::LoadFile(config_file).as(); + auto unfolded_data = biu::deserialize + (biu::read(input.UnfoldedDataFile)); + Output 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; + selected_modes.insert({qpoint.SourceIndex, i}); + } + } + + // 构造输出数据 + for (auto [i, j] : selected_modes) + { + auto& mode_data = output.ModeData.emplace_back(); + mode_data.MetaQpointIndex = i; + mode_data.ModeIndex = j; + // 未归一化的位移, 假定虚部总是为零 + auto atom_movement = + unfolded_data.MetaQpointData[i].ModeData[j].AtomMovement.real().cwiseProduct + ( + ( + input.AtomSymbols + | ranges::views::transform([&](const auto& symbol) + { return input.AtomMasses.at(symbol); }) + | ranges::to_vector + | biu::toEigen<> + ).cwiseSqrt().cwiseInverse().rowwise().replicate(3) + ).eval(); + // 归一化 + mode_data.AtomMovement = atom_movement / atom_movement.rowwise().norm().maxCoeff() * input.MaxDisplacement; + } + + // 输出 + std::ofstream(input.OutputDataFile, std::ios::binary) << biu::serialize(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/solver.cpp b/src/solver.cpp deleted file mode 100644 index 4d6ce5e..0000000 --- a/src/solver.cpp +++ /dev/null @@ -1,47 +0,0 @@ -# include - -namespace ufo -{ - concurrencpp::generator, unsigned>> Solver::triplet_sequence - (Eigen::Vector range) - { - for (unsigned x = 0; x < range[0]; x++) - for (unsigned y = 0; y < range[1]; y++) - for (unsigned z = 0; z < range[2]; z++) - co_yield - { - Eigen::Vector{{x}, {y}, {z}}, - x * range[1] * range[2] + y * range[2] + z - }; - } - - Solver::DataFile::DataFile - (YAML::Node node, std::set supported_format, std::string config_file, bool allow_same_as_config_file) - { - if (auto _ = node["SameAsConfigFile"]) - { - auto __ = _.as(); - 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(); - Format = node["Format"].as(); - if (!supported_format.contains(Format)) - throw std::runtime_error(fmt::format("Unsupported format: \"{}\"", Format)); - if (auto _ = node["RelativeToConfigFile"]) - { - auto __ = _.as(); - ExtraParameters["RelativeToConfigFile"] = __; - if (__) - Filename = std::filesystem::path(config_file).parent_path() / Filename; - } - }; - -} diff --git a/src/unfold.cpp b/src/unfold.cpp index e3eea4d..de20a03 100644 --- a/src/unfold.cpp +++ b/src/unfold.cpp @@ -1,384 +1,245 @@ -# include +# include +# include +# include +# include -namespace ufo +void ufo::unfold(std::string config_file) { - UnfoldSolver::InputType::InputType(std::string filename) + // 反折叠的原理: 将超胞中的原子运动状态, 投影到一组平面波构成的基矢中. + // 每一个平面波的波矢由两部分相加得到: 一部分是单胞倒格子的整数倍, 所取的个数有一定任意性, 论文中建议取大约单胞中原子个数那么多个; + // 对于没有缺陷的情况, 取一个应该就足够了. + // 这些平面波以原胞为周期。 + // 另一部分是超胞倒格子的整数倍, 取 n 个, n 为超胞对应的单胞的倍数, 其实也就是倒空间中单胞对应倒格子中超胞的格点. + // 只要第一部分取得足够多, 那么单胞中原子的状态就可以完全被这些平面波描述. + // 将超胞中原子的运动状态投影到这些基矢上, 计算出投影的系数, 就可以将超胞的原子运动状态分解到单胞中的多个 q 点上. + + struct Input { - // read main input file + // 单胞的三个格矢，每行表示一个格矢的坐标，单位为埃 + Eigen::Matrix3d PrimativeCell; + + // 单胞到超胞的格矢转换时用到的矩阵 + // SuperCellMultiplier 是一个三维列向量且各个元素都是整数，表示单胞在各个方向扩大到多少倍之后，可以得到和超胞一样的体积 + // SuperCellDeformation 是一个行列式为 1 的矩阵，它表示经过 SuperCellMultiplier 扩大后，还需要怎样的变换才能得到超胞 + // SuperCell = (SuperCellDeformation * SuperCellMultiplier.asDiagonal()) * PrimativeCell + // ReciprocalPrimativeCell = (SuperCellDeformation * SuperCellMultiplier.asDiagonal()).transpose() + // * ReciprocalSuperCell + // Position = PositionToCell(line vector) * Cell + // InversePosition = InversePositionToCell(line vector) * ReciprocalCell + // PositionToSuperCell(line vector) * SuperCell = PositionToPrimativeCell(line vector) * PrimativeCell + // ReciprocalPositionToSuperCell(line vector) * ReciprocalSuperCell + // = ReciprocalPositionToPrimativeCell(line vector) * ReciprocalPrimativeCell + Eigen::Matrix3d SuperCellDeformation; + Eigen::Vector3i SuperCellMultiplier; + + // 在单胞内取几个平面波的基矢 + Eigen::Vector PrimativeCellBasisNumber; + + // 超胞中原子的坐标，每行表示一个原子的坐标，单位为超胞的格矢 + Eigen::MatrixX3d AtomPositionBySuperCell; + + // 从 band.hdf5 读入 QpointData + std::optional QpointDataInputFile; + + // 输出到哪些文件 + struct QpointDataOutputFileType { - auto node = YAML::LoadFile(filename); - for (unsigned i = 0; i < 3; i++) - for (unsigned j = 0; j < 3; j++) - PrimativeCell(i, j) = node["PrimativeCell"][i][j].as(); + std::string Filename; - for (unsigned i = 0; i < 3; i++) - SuperCellMultiplier(i) = node["SuperCellMultiplier"][i].as(); + // 如果指定，则将结果投影到那些原子上 + std::optional> SelectedAtoms; - if (auto value = node["SuperCellDeformation"]) + // 默认输出为 zpp 文件，如果指定为 true，则输出为 yaml 文件 + std::optional OutputAsYaml; + }; + std::vector QpointDataOutputFile; + }; + + // 从文件中读取 QpointData + auto read_qpoint_data = [](std::string filename) + { + // 读入原始数据 + // 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); + + if (file.File.getDataSet("/frequency").getDimensions().size() == 2) + file.read("/frequency", frequency) + .read("/eigenvector", eigenvector_vector) + .read("/qpoint", qpoint); + else + { + std::vector>> temp_path; + std::vector>> temp_frequency; + std::vector>>> temp_eigenvector_vector; + file.read("/frequency", temp_frequency) + .read("/eigenvector", temp_eigenvector_vector) + .read("/path", temp_path); + frequency = temp_frequency | ranges::views::join | ranges::to_vector; + qpoint = temp_path | ranges::views::join | ranges::to_vector; + eigenvector_vector = temp_eigenvector_vector | ranges::views::join | ranges::to_vector; + } + + // 整理得到结果 + auto number_of_qpoints = frequency.size(), num_of_modes = frequency[0].size(); + std::vector qpoint_data(number_of_qpoints); + for (std::size_t i = 0; i < number_of_qpoints; i++) + { + qpoint_data[i].Qpoint = qpoint[i] | biu::toEigen<>; + qpoint_data[i].ModeData.resize(num_of_modes); + for (std::size_t j = 0; j < num_of_modes; j++) { - SuperCellDeformation.emplace(); - for (unsigned i = 0; i < 3; i++) - for (unsigned j = 0; j < 3; j++) - (*SuperCellDeformation)(i, j) = value[i][j].as(); - } - - for (unsigned i = 0; i < 3; i++) - PrimativeCellBasisNumber(i) = node["PrimativeCellBasisNumber"][i].as(); - - AtomPositionInputFile = DataFile - ( - node["AtomPositionInputFile"], {"yaml"}, - filename, true - ); - QpointDataInputFile = DataFile - ( - node["QpointDataInputFile"], {"yaml", "hdf5"}, - filename, true - ); - if (auto value = node["QpointDataOutputFile"]) - { - QpointDataOutputFile.resize(value.size()); - for (unsigned i = 0; i < value.size(); i++) - QpointDataOutputFile[i] = DataFile - ( - value[i], {"yaml", "yaml-human-readable", "zpp", "hdf5"}, - filename, false - ); + qpoint_data[i].ModeData[j].Frequency = frequency[i][j]; + auto number_of_atoms = eigenvector_vector[i].size() / 3; + Eigen::MatrixX3cd eigenvectors(number_of_atoms, 3); + for (std::size_t k = 0; k < number_of_atoms; k++) for (std::size_t l = 0; l < 3; l++) + eigenvectors(k, l) + = eigenvector_vector[i][k * 3 + l][j].r + eigenvector_vector[i][k * 3 + l][j].i * 1i; + // 原则上讲，需要对读入的原子运动状态作相位转换, 使得它们与我们的约定一致(对超胞周期性重复)，但这个转换 phonopy 已经做了 + // 这里还要需要做归一化处理 (指将数据简单地作为向量处理的归一化) + qpoint_data[i].ModeData[j].AtomMovement = eigenvectors / eigenvectors.norm(); } } + return qpoint_data; + }; - if (AtomPositionInputFile.Format == "yaml") - { - auto node = YAML::LoadFile(AtomPositionInputFile.Filename); - std::vector points; - if (auto _ = node["points"]) - points = _.as>(); - else - points = node["unit_cell"]["points"].as>(); - auto atom_position_to_super_cell = Eigen::MatrixX3d(points.size(), 3); - for (unsigned i = 0; i < points.size(); i++) - for (unsigned j = 0; j < 3; j++) - atom_position_to_super_cell(i, j) = points[i]["coordinates"][j].as(); - auto super_cell = (SuperCellDeformation.value_or(Eigen::Matrix3d::Identity()) - * SuperCellMultiplier.cast().asDiagonal() * PrimativeCell).eval(); - AtomPosition = atom_position_to_super_cell * super_cell; - } - if (QpointDataInputFile.Format == "yaml") - { - auto node = YAML::LoadFile(QpointDataInputFile.Filename); - auto phonon = node["phonon"].as>(); - QpointData.resize(phonon.size()); - for (unsigned i = 0; i < phonon.size(); i++) - { - for (unsigned j = 0; j < 3; j++) - QpointData[i].Qpoint(j) = phonon[i]["q-position"][j].as(); - auto band = phonon[i]["band"].as>(); - QpointData[i].ModeData.resize(band.size()); - for (unsigned j = 0; j < band.size(); j++) - { - QpointData[i].ModeData[j].Frequency = band[j]["frequency"].as(); - auto eigenvector_vectors = band[j]["eigenvector"] - .as>>>(); - 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 * 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>> frequency, path; - std::vector>>> eigenvector_vector; - Hdf5file{}.open_for_read(QpointDataInputFile.Filename).read(frequency, "/frequency") - .read(eigenvector_vector, "/eigenvector") - .read(path, "/path"); - std::vector size = { frequency.size(), frequency[0].size(), frequency[0][0].size() }; - QpointData.resize(size[0] * size[1]); - for (unsigned i = 0; i < size[0]; i++) - for (unsigned j = 0; j < size[1]; j++) - { - QpointData[i * size[1] + j].Qpoint = Eigen::Vector3d(path[i][j].data()); - QpointData[i * size[1] + j].ModeData.resize(size[2]); - for (unsigned k = 0; k < size[2]; k++) - { - QpointData[i * size[1] + j].ModeData[k].Frequency = frequency[i][j][k]; - Eigen::MatrixX3cd eigenvectors(AtomPosition.rows(), 3); - for (unsigned l = 0; l < AtomPosition.rows(); l++) - for (unsigned m = 0; m < 3; m++) - eigenvectors(l, m) - = eigenvector_vector[i][j][l * 3 + m][k].r + eigenvector_vector[i][j][l * 3 + m][k].i * 1i; - QpointData[i * size[1] + j].ModeData[k].AtomMovement = eigenvectors / eigenvectors.norm(); - } - } - } - } - - void UnfoldSolver::OutputType::write - (decltype(InputType::QpointDataOutputFile) output_files) const - { - for (auto& output_file : output_files) - write(output_file.Filename, output_file.Format); - } - void UnfoldSolver::OutputType::write(std::string filename, std::string format, unsigned percision) const - { - if (format == "yaml") - std::ofstream(filename) << [&] - { - std::stringstream print; - print << "QpointData:\n"; - for (auto& qpoint: QpointData) - { - print << fmt::format(" - Qpoint: [ {1:.{0}f}, {2:.{0}f}, {3:.{0}f} ]\n", - percision, qpoint.Qpoint[0], qpoint.Qpoint[1], qpoint.Qpoint[2]); - print << fmt::format(" Source: [ {1:.{0}f}, {2:.{0}f}, {3:.{0}f} ]\n", - percision, qpoint.Source[0], qpoint.Source[1], qpoint.Source[2]); - print << " ModeData:\n"; - for (auto& mode: qpoint.ModeData) - print << fmt::format(" - {{ Frequency: {1:.{0}f}, Weight: {2:.{0}f} }}\n", - percision, mode.Frequency, mode.Weight); - } - return print.str(); - }(); - else if (format == "yaml-human-readable") - { - std::remove_cvref_t output; - std::map> - 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 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> Qpoint, Source, Frequency, Weight; - for (auto& qpoint : QpointData) - { - Qpoint.emplace_back(qpoint.Qpoint.data(), qpoint.Qpoint.data() + 3); - Source.emplace_back(qpoint.Source.data(), qpoint.Source.data() + 3); - Frequency.emplace_back(); - Weight.emplace_back(); - for (auto& mode : qpoint.ModeData) - { - Frequency.back().push_back(mode.Frequency); - Weight.back().push_back(mode.Weight); - } - } - Hdf5file{}.open_for_write(filename).write(Qpoint, "/Qpoint") - .write(Source, "/Source") - .write(Frequency, "/Frequency") - .write(Weight, "/Weight"); - } - } - - UnfoldSolver::UnfoldSolver(std::string config_file) : Input_([&] - { - std::clog << "Reading input file... " << std::flush; - return config_file; - }()) - { - std::clog << "Done." << std::endl; - } - - UnfoldSolver& UnfoldSolver::operator()() - { - if (!Basis_) - { - std::clog << "Constructing basis... " << std::flush; - Basis_ = construct_basis - ( - Input_.PrimativeCell, Input_.SuperCellMultiplier, - Input_.PrimativeCellBasisNumber, Input_.AtomPosition - ); - std::clog << "Done." << std::endl; - } - if (!Output_) - { - std::clog << "Calculating projection coefficient... " << std::flush; - std::vector> mode_data; - for (auto& qpoint : Input_.QpointData) - for (auto& mode : qpoint.ModeData) - mode_data.emplace_back(mode.AtomMovement); - std::atomic 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> qpoint; - std::vector>> frequency; - for (auto& qpoint_data : Input_.QpointData) - { - qpoint.emplace_back(qpoint_data.Qpoint); - frequency.emplace_back(); - for (auto& mode_data : qpoint_data.ModeData) - frequency.back().emplace_back(mode_data.Frequency); - } - Output_ = construct_output - ( - Input_.SuperCellMultiplier, - Input_.SuperCellDeformation, qpoint, frequency, projection_coefficient - ); - std::clog << "Done." << std::endl; - } - std::clog << "Writing output... " << std::flush; - Output_->write(Input_.QpointDataOutputFile); - std::clog << "Done." << std::endl; - return *this; - } - - UnfoldSolver::BasisType UnfoldSolver::construct_basis + // 构建基 + // 每个 q 点对应一组 sub qpoint。不同的 q 点所对应的 sub qpoint 是不一样的，但 sub qpoint 与 q 点的相对位移在不同 q 点之间是相同的。 + // 由于基只与这个相对位置有关（也就是说，不同 q 点的基是一样的），因此可以先计算出所有的基，这样降低计算量。 + // 外层下标对应超胞倒格子的整数倍那部分(第二部分), 也就是不同的 sub qpoint + // 内层下标对应单胞倒格子的整数倍那部分(第一部分), 也就是 sub qpoint 上的不同平面波（取的数量越多，结果越精确） + auto construct_basis = [] ( - const decltype(InputType::PrimativeCell)& primative_cell, - const decltype(InputType::SuperCellMultiplier)& super_cell_multiplier, - const decltype(InputType::PrimativeCellBasisNumber)& primative_cell_basis_number, - const decltype(InputType::AtomPosition)& atom_position + Eigen::Matrix3d primative_cell, Eigen::Vector3i super_cell_multiplier, + Eigen::Vector primative_cell_basis_number, Eigen::MatrixX3d atom_position ) { - BasisType basis(super_cell_multiplier.prod()); - // 每个 q 点对应的一组 sub qpoint。不同的 q 点所对应的 sub qpoint 是不一样的，但 sub qpoint 与 q 点的相对位置一致。 - // 这里 xyz_of_diff_of_sub_qpoint 即表示这个相对位置，单位为超胞的倒格矢 - for (auto [xyz_of_diff_of_sub_qpoint_by_reciprocal_modified_super_cell, i_of_sub_qpoint] - : triplet_sequence(super_cell_multiplier)) + biu::Logger::Guard log; + std::vector> basis(super_cell_multiplier.prod()); + // diff_of_sub_qpoint 表示 sub qpoint 与 qpoint 的相对位置，单位为超胞的倒格矢 + for (auto [diff_of_sub_qpoint_by_reciprocal_modified_super_cell, i_of_sub_qpoint] + : biu::sequence(super_cell_multiplier)) { basis[i_of_sub_qpoint].resize(primative_cell_basis_number.prod()); - for (auto [xyz_of_basis, i_of_basis] : triplet_sequence(primative_cell_basis_number)) + for (auto [xyz_of_basis, i_of_basis] + : biu::sequence(primative_cell_basis_number)) { // 计算 q 点的坐标, 单位为单胞的倒格矢 auto diff_of_sub_qpoint_by_reciprocal_primative_cell = xyz_of_basis.cast() + super_cell_multiplier.cast().cwiseInverse().asDiagonal() - * xyz_of_diff_of_sub_qpoint_by_reciprocal_modified_super_cell.cast(); - // 将 q 点坐标转换为埃^-1 - auto qpoint = (diff_of_sub_qpoint_by_reciprocal_primative_cell.transpose() + * diff_of_sub_qpoint_by_reciprocal_modified_super_cell.cast(); + // 将单位转换为埃^-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 * (atom_position * qpoint)).array().exp(); + = (2i * std::numbers::pi_v * (atom_position * diff_of_sub_qpoint)).array().exp(); } } return basis; - } + }; - std::vector> UnfoldSolver::construct_projection_coefficient + // 计算从超胞到原胞的投影系数（不是分原子的投影系数），是反折叠的核心步骤 + // 返回的投影系数是一个三维数组，第一维对应不同的 q 点，第二维对应不同的模式，第三维对应不同的 sub qpoint + auto construct_projection_coefficient = [] ( - const BasisType& basis, - const std::vector>& mode_data, - std::atomic& number_of_finished_modes + const std::vector>& basis, + // 实际上只需要其中的 AtomMovement + const std::vector& qpoint_data, + 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::transform ( - std::execution::par, mode_data.begin(), mode_data.end(), + std::execution::par_unseq, mode_data.begin(), mode_data.end(), projection_coefficient.begin(), [&](const auto& mode_data) { // 这里, mode_data 和 projection_coefficient 均指对应于一个模式的数据 std::vector projection_coefficient(basis.size()); - for (unsigned i_of_sub_qpoint = 0; i_of_sub_qpoint < basis.size(); i_of_sub_qpoint++) + for (std::size_t i_of_sub_qpoint = 0; i_of_sub_qpoint < basis.size(); i_of_sub_qpoint++) // 对于 basis 中, 对应于单胞倒格子的部分, 以及对应于不同方向的部分, 分别求内积, 然后求模方和 - for (unsigned i_of_basis = 0; i_of_basis < basis[i_of_sub_qpoint].size(); i_of_basis++) + for (std::size_t i_of_basis = 0; i_of_basis < basis[i_of_sub_qpoint].size(); i_of_basis++) projection_coefficient[i_of_sub_qpoint] += (basis[i_of_sub_qpoint][i_of_basis].transpose().conjugate() * mode_data.get()) .array().abs2().sum(); // 如果是严格地将向量分解到一组完备的基矢上, 那么不需要对计算得到的权重再做归一化处理 // 但这里并不是这样一个严格的概念. 因此对分解到各个 sub qpoint 上的权重做归一化处理 - auto sum = std::accumulate - (projection_coefficient.begin(), projection_coefficient.end(), 0.); - for (auto& _ : projection_coefficient) - _ /= sum; + auto sum = ranges::accumulate(projection_coefficient, 0.); + for (auto& _ : projection_coefficient) _ /= sum; number_of_finished_modes++; return projection_coefficient; } ); - return projection_coefficient; - } - - UnfoldSolver::OutputType UnfoldSolver::construct_output - ( - const decltype(InputType::SuperCellMultiplier)& super_cell_multiplier, - const decltype(InputType::SuperCellDeformation)& super_cell_deformation, - const std::vector>& meta_qpoint_by_reciprocal_super_cell, - const std::vector>>& frequency, - const ProjectionCoefficientType_& projection_coefficient - ) - { - OutputType output; + // 将计算得到的投影系数重新组装成三维数组 + // 第一维是 meta qpoint，第二维是模式，第三维是 sub qpoint + std::vector>> projection_coefficient_output; for ( - unsigned i_of_meta_qpoint = 0, num_of_mode_manipulated = 0; - i_of_meta_qpoint < meta_qpoint_by_reciprocal_super_cell.size(); - i_of_meta_qpoint++ + std::size_t i_of_meta_qpoint = 0, num_of_mode_manipulated = 0; + i_of_meta_qpoint < qpoint_data.size(); + i_of_meta_qpoint++, num_of_mode_manipulated += qpoint_data[i_of_meta_qpoint].ModeData.size() ) + projection_coefficient_output.emplace_back + ( + projection_coefficient.begin() + num_of_mode_manipulated, + projection_coefficient.begin() + num_of_mode_manipulated + qpoint_data[i_of_meta_qpoint].ModeData.size() + ); + return projection_coefficient_output; + }; + + // 组装输出，即将投影系数应用到原始数据上 + auto construct_output = [] + ( + const Input& input, + const std::vector>>& 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++) { - for (auto [xyz_of_diff_of_sub_qpoint_by_reciprocal_modified_super_cell, i_of_sub_qpoint] - : triplet_sequence(super_cell_multiplier)) + // 如果需要投影到特定的原子上，需要先计算当前 meta qpoint 的不同模式的投影系数 + std::optional> 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(); /* @@ -404,26 +265,81 @@ namespace ufo */ auto sub_qpoint_by_reciprocal_primative_cell = ( - super_cell_multiplier.cast().cwiseInverse().asDiagonal() + input.SuperCellMultiplier.cast().cwiseInverse().asDiagonal() * ( - xyz_of_diff_of_sub_qpoint_by_reciprocal_modified_super_cell.cast() - + super_cell_deformation.value_or(Eigen::Matrix3d::Identity()).inverse() - * meta_qpoint_by_reciprocal_super_cell[i_of_meta_qpoint].get().cast() + 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 = meta_qpoint_by_reciprocal_super_cell[i_of_meta_qpoint]; - _.SourceIndex_ = i_of_meta_qpoint; - for (unsigned i_of_mode = 0; i_of_mode < frequency[i_of_meta_qpoint].size(); i_of_mode++) + _.Source = qpoint_data[i_of_meta_qpoint].Qpoint; + _.SourceIndex = i_of_meta_qpoint; + + for (std::size_t i_of_mode = 0; i_of_mode < qpoint_data[i_of_meta_qpoint].ModeData.size(); i_of_mode++) { auto& __ = _.ModeData.emplace_back(); - __.Frequency = frequency[i_of_meta_qpoint][i_of_mode]; - __.Weight = projection_coefficient[num_of_mode_manipulated + i_of_mode][i_of_sub_qpoint]; + __.Frequency = qpoint_data[i_of_meta_qpoint].ModeData[i_of_mode].Frequency; + __.Weight = projection_coefficient[i_of_meta_qpoint][i_of_mode][i_of_sub_qpoint]; + if (selected_atoms) + __.Weight *= projection_coefficient_on_atoms.value()[i_of_mode]; } } - num_of_mode_manipulated += frequency[i_of_meta_qpoint].size(); } return output; + }; + + biu::Logger::Guard log; + log.info("Reading input file... "); + auto input = YAML::LoadFile(config_file).as(); + auto qpoint_data = read_qpoint_data(input.QpointDataInputFile.value_or("band.hdf5")); + log.info("Done."); + + std::clog << "Constructing basis... " << std::flush; + + auto basis = construct_basis + ( + input.PrimativeCell, input.SuperCellMultiplier, + input.PrimativeCellBasisNumber, + input.AtomPositionBySuperCell + * (input.SuperCellDeformation * input.SuperCellMultiplier.cast().asDiagonal() * input.PrimativeCell) + ); + std::clog << "Done." << std::endl; + + 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) + { + std::osyncstream(std::clog) + << "\rCalculating projection coefficient... ({}/{})"_f(number_of_finished_modes, number_of_modes) + << std::flush; + std::this_thread::sleep_for(100ms); + if (finished) break; + } + }); + auto projection_coefficient = construct_projection_coefficient(basis, qpoint_data, number_of_finished_modes); + finished = true; + print_thread.join(); + std::clog << "\33[2K\rCalculating projection coefficient... Done." << std::endl; + + std::clog << "Writing data... " << std::flush; + for (auto& output_file : input.QpointDataOutputFile) + { + auto output = construct_output + (input, projection_coefficient, qpoint_data, output_file.SelectedAtoms); + if (output_file.OutputAsYaml.value_or(false)) std::ofstream(output_file.Filename) << YAML::Node(output); + else std::ofstream(output_file.Filename, std::ios::binary) << biu::serialize(output); } + std::clog << "Done." << std::endl; } diff --git a/test/14.2.2.4.unfold.yaml b/test/14.2.2.4.unfold.yaml deleted file mode 100644 index 672aa0e..0000000 --- a/test/14.2.2.4.unfold.yaml +++ /dev/null @@ -1,12 +0,0 @@ -PrimativeCell: -- [ 1.548010265167714, -2.681232429908652, 0.000000000000000 ] # a -- [ 1.548010265167714, 2.681232429908652, 0.000000000000000 ] # b -- [ 0.000000000000000, 0.000000000000000, 5.061224103556595 ] # c -SuperCellMultiplier: [ 3, 1, 1 ] -PrimativeCellBasisNumber: [ 8, 8, 8 ] -AtomPositionInputFile: { FileName: "/home/chn/Documents/lammps-SiC/14/14.2/14.2.2/14.2.2.4/band.yaml", Format: "yaml" } -QpointDataInputFile: { FileName: "/home/chn/Documents/lammps-SiC/14/14.2/14.2.2/14.2.2.4/band.yaml", Format: "yaml" } -QpointDataOutputFile: -- { FileName: "test/14.2.2.4.result.yaml", Format: "yaml" } -- { FileName: "test/14.2.2.4.result.human-readable.yaml", Format: "yaml-human-readable" } -- { FileName: "test/14.2.2.4.result.zpp", Format: "zpp" } diff --git a/test/14.2.2.5.plot.yaml b/test/14.2.2.5.plot.yaml deleted file mode 100644 index 6b445a2..0000000 --- a/test/14.2.2.5.plot.yaml +++ /dev/null @@ -1,14 +0,0 @@ -PrimativeCell: -- [ 1.548010265167714, -2.681232429908652, 0.000000000000000 ] # a -- [ 1.548010265167714, 2.681232429908652, 0.000000000000000 ] # b -- [ 0.000000000000000, 0.000000000000000, 5.061224103556595 ] # c -SourceFilename: ./test/14.2.2.5.result.zpp -Figures: -- Qpoints: - - - [ 0, 0, 0 ] - - [ 0.5, 0, 0 ] - - [ 0.3333333, 0.3333333, 0 ] - - [ 0, 0, 0 ] - Resolution: [ 1024, 400 ] - Range: [ -5, 35 ] - Filename: ./test/14.2.2.5.png diff --git a/test/14.2.2.5.unfold.yaml b/test/14.2.2.5.unfold.yaml deleted file mode 100644 index e90577c..0000000 --- a/test/14.2.2.5.unfold.yaml +++ /dev/null @@ -1,12 +0,0 @@ -PrimativeCell: -- [ 1.548010265167714, -2.681232429908652, 0.000000000000000 ] # a -- [ 1.548010265167714, 2.681232429908652, 0.000000000000000 ] # b -- [ 0.000000000000000, 0.000000000000000, 5.061224103556595 ] # c -SuperCellMultiplier: [ 3, 1, 1 ] -PrimativeCellBasisNumber: [ 8, 8, 8 ] -AtomPositionInputFile: { FileName: "/home/chn/Documents/lammps-SiC/14/14.2/14.2.2/14.2.2.5/band.yaml", Format: "yaml" } -QpointDataInputFile: { FileName: "/home/chn/Documents/lammps-SiC/14/14.2/14.2.2/14.2.2.5/band.yaml", Format: "yaml" } -QpointDataOutputFile: -- { FileName: "test/14.2.2.5.result.yaml", Format: "yaml" } -- { FileName: "test/14.2.2.5.result.human-readable.yaml", Format: "yaml-human-readable" } -- { FileName: "test/14.2.2.5.result.zpp", Format: "zpp" } diff --git a/test/14.2.4.4.plot.yaml b/test/14.2.4.4.plot.yaml deleted file mode 100644 index 6c812ae..0000000 --- a/test/14.2.4.4.plot.yaml +++ /dev/null @@ -1,14 +0,0 @@ -PrimativeCell: -- [ 1.548010265167714, -2.681232429908652, 0.000000000000000 ] # a -- [ 1.548010265167714, 2.681232429908652, 0.000000000000000 ] # b -- [ 0.000000000000000, 0.000000000000000, 5.061224103556595 ] # c -SourceFilename: ./test/14.2.4.4.result.zpp -Figures: -- Qpoints: - - - [ 0, 0, 0 ] - - [ 0.5, 0, 0 ] - - [ 0.3333333, 0.3333333, 0 ] - - [ 0, 0, 0 ] - Resolution: [ 1024, 400 ] - Range: [ -5, 35 ] - Filename: ./test/14.2.4.4.png diff --git a/test/14.2.4.4.unfold.yaml b/test/14.2.4.4.unfold.yaml deleted file mode 100644 index 4c21443..0000000 --- a/test/14.2.4.4.unfold.yaml +++ /dev/null @@ -1,14 +0,0 @@ -PrimativeCell: -- [ 1.548010265167714, -2.681232429908652, 0.000000000000000 ] # a -- [ 1.548010265167714, 2.681232429908652, 0.000000000000000 ] # b -- [ 0.000000000000000, 0.000000000000000, 5.061224103556595 ] # c -SuperCellMultiplier: [ 2, 2, 1 ] -PrimativeCellBasisNumber: [ 8, 8, 8 ] -AtomPositionInputFile: - FileName: "/home/chn/Documents/lammps-SiC/14/14.2/14.2.4/14.2.4.4/phonopy.yaml" - Format: "yaml" -QpointDataInputFile: { FileName: "/home/chn/Documents/lammps-SiC/14/14.2/14.2.4/14.2.4.4/band.hdf5", Format: "hdf5" } -QpointDataOutputFile: -- { FileName: "test/14.2.4.4.result.yaml", Format: "yaml" } -- { FileName: "test/14.2.4.4.result.human-readable.yaml", Format: "yaml-human-readable" } -- { FileName: "test/14.2.4.4.result.zpp", Format: "zpp" } diff --git a/test/14.2.5.5.plot.yaml b/test/14.2.5.5.plot.yaml deleted file mode 100644 index ad219e2..0000000 --- a/test/14.2.5.5.plot.yaml +++ /dev/null @@ -1,14 +0,0 @@ -PrimativeCell: -- [ 1.548010265167714, -2.681232429908652, 0.000000000000000 ] # a -- [ 1.548010265167714, 2.681232429908652, 0.000000000000000 ] # b -- [ 0.000000000000000, 0.000000000000000, 5.061224103556595 ] # c -SourceFilename: ./test/14.2.5.5.result.zpp -Figures: -- Qpoints: - - - [ 0, 0, 0 ] - - [ 0.5, 0, 0 ] - - [ 0.3333333, 0.3333333, 0 ] - - [ 0, 0, 0 ] - Resolution: [ 1024, 400 ] - Range: [ -5, 35 ] - Filename: ./test/14.2.5.5.png diff --git a/test/14.2.5.5.unfold.yaml b/test/14.2.5.5.unfold.yaml deleted file mode 100644 index c19c01b..0000000 --- a/test/14.2.5.5.unfold.yaml +++ /dev/null @@ -1,18 +0,0 @@ -PrimativeCell: -- [ 1.548010265167714, -2.681232429908652, 0.000000000000000 ] # a -- [ 1.548010265167714, 2.681232429908652, 0.000000000000000 ] # b -- [ 0.000000000000000, 0.000000000000000, 5.061224103556595 ] # c -SuperCellMultiplier: [ 1, 2, 1 ] -SuperCellDeformation: -- [ 1, 0, 0] -- [ 1, 1, 0] -- [ 0, 0, 1] -PrimativeCellBasisNumber: [ 8, 8, 8 ] -AtomPositionInputFile: - FileName: "/home/chn/Documents/lammps-SiC/14/14.2/14.2.5/14.2.5.5/phonopy.yaml" - Format: "yaml" -QpointDataInputFile: { FileName: "/home/chn/Documents/lammps-SiC/14/14.2/14.2.5/14.2.5.5/band.hdf5", Format: "hdf5" } -QpointDataOutputFile: -- { FileName: "test/14.2.5.5.result.yaml", Format: "yaml" } -- { FileName: "test/14.2.5.5.result.human-readable.yaml", Format: "yaml-human-readable" } -- { FileName: "test/14.2.5.5.result.zpp", Format: "zpp" } diff --git a/test/14.2.6.4.plot.yaml b/test/14.2.6.4.plot.yaml deleted file mode 100644 index 810c926..0000000 --- a/test/14.2.6.4.plot.yaml +++ /dev/null @@ -1,15 +0,0 @@ -PrimativeCell: -- [ 1.548010265167714, -2.681232429908652, 0.000000000000000 ] # a -- [ 1.548010265167714, 2.681232429908652, 0.000000000000000 ] # b -- [ 0.000000000000000, 0.000000000000000, 5.061224103556595 ] # c -SourceFilename: ./test/14.2.6.4.result.zpp -Figures: -- Qpoints: - - - [ 0, 0, 0 ] - - [ 0.5, 0, 0 ] - - [ 0.3333333, 0.3333333, 0 ] - - [ 0, 0, 0 ] - Resolution: [ 1024, 1024 ] - Range: [ -5, 35 ] - Filename: ./test/14.2.6.4.png - YTicks: [ 0, 5, 10, 15, 20, 25, 30 ] diff --git a/test/14.2.6.4.unfold.yaml b/test/14.2.6.4.unfold.yaml deleted file mode 100644 index 18b8695..0000000 --- a/test/14.2.6.4.unfold.yaml +++ /dev/null @@ -1,18 +0,0 @@ -PrimativeCell: -- [ 1.548010265167714, -2.681232429908652, 0.000000000000000 ] # a -- [ 1.548010265167714, 2.681232429908652, 0.000000000000000 ] # b -- [ 0.000000000000000, 0.000000000000000, 5.061224103556595 ] # c -SuperCellMultiplier: [ 4, 8, 1 ] -SuperCellDeformation: -- [ 1, 0, 0] -- [ 1, 1, 0] -- [ 0, 0, 1] -PrimativeCellBasisNumber: [ 8, 8, 8 ] -AtomPositionInputFile: - FileName: "/home/chn/Documents/lammps-SiC/14/14.2/14.2.6/14.2.6.4/phonopy.yaml" - Format: "yaml" -QpointDataInputFile: { FileName: "/home/chn/Documents/lammps-SiC/14/14.2/14.2.6/14.2.6.4/band.hdf5", Format: "hdf5" } -QpointDataOutputFile: -- { FileName: "test/14.2.6.4.result.yaml", Format: "yaml" } -- { FileName: "test/14.2.6.4.result.human-readable.yaml", Format: "yaml-human-readable" } -- { FileName: "test/14.2.6.4.result.zpp", Format: "zpp" } diff --git a/test/fold/config.yaml b/test/fold/config.yaml new file mode 100644 index 0000000..3bc8ff4 --- /dev/null +++ b/test/fold/config.yaml @@ -0,0 +1,18 @@ +SuperCellMultiplier: [3, 4, 1] +SuperCellDeformation: + - [ 1, 0, 0 ] + - [ 0.6666, 1, 0 ] + - [ 0, 0, 1 ] +Qpoints: + - [0, 0, 0] + - [0.05, 0, 0] + - [0.1, 0, 0] + - [0.15, 0, 0] + - [0.2, 0, 0] + - [0.25, 0, 0] + - [0.3, 0, 0] + - [0.35, 0, 0] + - [0.4, 0, 0] + - [0.45, 0, 0] + - [0.5, 0, 0] +OutputFile: fold-output.yaml