diff --git a/include/ufo/unfold.hpp b/include/ufo/unfold.hpp index a652e29..cef8c10 100644 --- a/include/ufo/unfold.hpp +++ b/include/ufo/unfold.hpp @@ -57,11 +57,13 @@ namespace ufo { // 模式的频率，单位为 THz double Frequency; - // 模式中各个原子的运动状态 - // 这个数据是这样得到的: phonopy 输出的动态矩阵的 eigenvector 乘以 $\exp(-2 \pi i \vec q \cdot \vec r)$ - // 这个数据可以认为是原子位移中, 关于超胞有周期性的那一部分, 再乘以原子质量的开方. - // 这个数据在读入后会被立即归一化. - Eigen::MatrixX3cd AtomMovement; + /* + phonopy 中定义的动态矩阵的特征向量, 代表模式中各个原子的运动状态 + (也即原子位移中, 关于超胞有周期性的那一部分, 再乘以原子质量的开方). + phonopy 的动态矩阵的定义与一般的定义不同, 因此这里的这个特征向量不是一般意义上的动态矩阵的特征向量. + 尽管在计算时, 确实需要将这个特征向量转化为标准的动态矩阵的特征向量, 但为了节省计算量, 在反折叠时一起计算, 而不是在输入时计算. + */ + Eigen::MatrixX3cd PhonopyEigenvector; }; std::vector ModeData; }; @@ -118,54 +120,41 @@ namespace ufo 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 - ); + // 计算单个 MetaQpoint 到所有 SubQpoint 的投影系数, 是反折叠的核心步骤 + struct ConstructProjectionCoefficientInputType + { + struct + { + const Eigen::Matrix3d& PrimativeCell; + const Eigen::Vector& SuperCellMultiplier; + const std::optional>& SuperCellDeformation; + const Eigen::Vector& PrimativeCellBasisNumber; + const Eigen::MatrixX3d& AtomPosition; + } GlobalData; + struct + { + const Eigen::Vector3d& MetaQpointByReciprocalSuperCell; + const std::vector& ModeData; + } MetaQpointData; + }; + struct ConstructProjectionCoefficientOutputType + { + struct SubQpointDataType + { + Eigen::Vector3d SubQpointByReciprocalPrimativeCell; + std::vector ModeData; + }; + std::vector SubQpointData; + }; + static ConstructProjectionCoefficientOutputType construct_projection_coefficient + (ConstructProjectionCoefficientInputType Input); }; } diff --git a/src/unfold.cpp b/src/unfold.cpp index 1f4add4..b0e6d9c 100644 --- a/src/unfold.cpp +++ b/src/unfold.cpp @@ -284,16 +284,6 @@ namespace ufo 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; @@ -344,6 +334,25 @@ namespace ufo return *this; } + auto UnfoldSolver::construct_projection_coefficient(ConstructProjectionCoefficientInputType Input) + -> ConstructProjectionCoefficientOutputType + { + /* + Eigenvector = PhonopyEigenvector * exp(2i * pi * MetaQpoint * AtomPosition); + Basis = exp(2i * pi * (SubQpoint + ABunchOfPlanWaveWhichIsTimesOfReciprocalPrimativeCell) * AtomPosition); + ProjectionCoefficient = sum_over_all_basis_about_single_sub_qpoint + ((Basis.conjugate().transpose() * Eigenvector).abs2().sum()); + SubQpoint = XyzOfDiffOfSubQpoint + MetaQpoint; + 整理可以得到: + ProjectionCoefficient = sum_over_all_basis_about_single_sub_qpoint + ((exp(2i * pi * (XyzOfDiffOfSubQpoint + ABunchOfPlanWaveWhichIsTimesOfReciprocalPrimativeCell) * AtomPosition) + .conjugate().transpose() * PhonopyEigenvector).abs2().sum()); + */ + + + } + + UnfoldSolver::BasisType UnfoldSolver::construct_basis ( const decltype(InputType::PrimativeCell)& primative_cell, @@ -446,14 +455,13 @@ namespace ufo ModifiedSuperCell = SuperCellMultiplier.asDiagonal() * PrimativeCell; MetaQpoint = MetaQpointByReciprocalModifiedSuperCell.transpose() * ReciprocalModifiedSuperCell; MetaQpoint = MetaQpointByReciprocalSuperCell.transpose() * ReciprocalSuperCell; - ReciprocalModifiedSuperCell = ModifiedSuperCell.inverse().transpose(); ReciprocalSuperCell = SuperCell.inverse().transpose(); ModifiedSuperCell = SuperCellDeformation * SuperCell; SuperCell = SuperCellMultiplier.asDiagonal() * PrimativeCell; 整理可以得到: SubQpointByReciprocalPrimativeCell = SuperCellMultiplier.asDiagonal().inverse() * (XyzOfDiffOfSubQpointByReciprocalModifiedSuperCell + - SuperCellDeformation * MetaQpointByReciprocalSuperCell); + SuperCellDeformation.inverse() * MetaQpointByReciprocalSuperCell); 但注意到, 这样得到的 SubQpoint 可能不在 ReciprocalPrimativeCell 中 (当 SuperCellDeformation 不是单位矩阵时, 边界附近的一两条 SubQpoint 会出现这种情况). 解决办法是, 在赋值时, 仅取 SubQpointByReciprocalPrimativeCell 的小数部分. @@ -463,7 +471,7 @@ namespace ufo super_cell_multiplier.cast().cwiseInverse().asDiagonal() * ( xyz_of_diff_of_sub_qpoint_by_reciprocal_modified_super_cell.cast() - + super_cell_deformation.value_or(Eigen::Matrix3d::Identity()) + + super_cell_deformation.value_or(Eigen::Matrix3d::Identity()).inverse() * meta_qpoint_by_reciprocal_super_cell[i_of_meta_qpoint].get().cast() ) ).eval();