diff --git a/src/raman.cpp b/src/raman.cpp index 7ea2789..15ab759 100644 --- a/src/raman.cpp +++ b/src/raman.cpp @@ -69,6 +69,10 @@ namespace ufo auto process = [&](auto& cell) { + // 我们需要做两件事。 + // 第一件事,是将原子振动幅度进行缩放,使得最大位移为指定值。这是为了保证计算时用到的位移幅度不是太大也不是太小,有个令人满意的精确度。 + // 第二件事是,计算这个振幅下对应的粒子数(声子数),并将它的倒数存储为 ratio, + // 这样在后续计算 Raman 张量时就可以直接用这个 ratio 来除以位移幅度,从而得到正确的 Raman 张量。 auto mass = cell.AtomType | ranges::views::transform([&](auto&& atom) { return ranges::views::repeat_n(input.AtomMass[atom.first], atom.second); }) @@ -76,12 +80,15 @@ namespace ufo for (auto i_of_qpoint : config.SelectedModes | ranges::views::keys) for (auto i_of_mode : config.SelectedModes[i_of_qpoint]) { - // 假定虚部总是为零 - auto move_eigenvector = cell.Qpoint[i_of_qpoint].Mode[i_of_mode].EigenVector.real() + // 假定虚部总是为零,将特征向量的实部取出,除以 sqrt(mass) 得到位移,这时位移还没有经过放缩 + auto movement_before_scale = cell.Qpoint[i_of_qpoint].Mode[i_of_mode].EigenVector.real() .cwiseProduct(mass.cwiseSqrt().cwiseInverse().rowwise().replicate(3)).eval(); - // 归一化 - auto ratio = config.MaxDisplacement / move_eigenvector.rowwise().norm().maxCoeff(); - auto atom_movement = (move_eigenvector * ratio).eval(); + // 对位移进行放缩 + auto movement_scale_ratio = config.MaxDisplacement / movement_before_scale.rowwise().norm().maxCoeff(); + auto movement_after_scale = (movement_before_scale * movement_scale_ratio).eval(); + // 计算这个位移下的粒子数,只需要是一个比例关系即可 + auto phonon_number = movement_after_scale.rowwise().squaredNorm().cwiseProduct(mass).sum() + * cell.Qpoint[i_of_qpoint].Mode[i_of_mode].Frequency; // 输出 for (auto direction : {1, -1}) { @@ -89,10 +96,11 @@ namespace ufo (config.OutputPoscarDirectory, i_of_qpoint, i_of_mode, direction > 0 ? "+" : "-"); std::filesystem::create_directories(path); std::ofstream("{}/POSCAR"_f(path)) << generate_poscar - (cell.Cell, cell.AtomPosition + atom_movement * cell.Cell.inverse() * direction, cell.AtomType); + (cell.Cell, cell.AtomPosition + movement_after_scale * cell.Cell.inverse() * direction, cell.AtomType); } - output.ModeRatio[{i_of_qpoint, i_of_mode}] = ratio; - log.debug("Write mode {} {} {}"_f(i_of_qpoint, i_of_mode, atom_movement.rowwise().norm().eval())); + output.ModeRatio[{i_of_qpoint, i_of_mode}] = 1 / phonon_number; + log.debug("Write mode {} {} {}"_f + (i_of_qpoint, i_of_mode, movement_after_scale.rowwise().norm().eval(), 1 / phonon_number)); } }; if (config.Cell == RamanData::Primative) process(input.Primative); @@ -179,8 +187,7 @@ namespace ufo auto&& [i_of_qpoint, i_of_mode] = i; auto&& _ = cell.Qpoint[i_of_qpoint].Mode[i_of_mode]; auto&& susceptibility = config.Susceptibility[i_of_qpoint][i_of_mode]; - Eigen::Matrix3d raman_tensor = (susceptibility["+"] - susceptibility["-"]) / 2 - / ratio / raman_input.MaxDisplacement; + Eigen::Matrix3d raman_tensor = (susceptibility["+"] - susceptibility["-"]) / 2 * ratio; _.RamanTensor = raman_tensor | biu::fromEigen; _.WeightOnRaman = config.Polarization[0].transpose() * raman_tensor * config.Polarization[1]; log.info("{}:{:.2f}:{}:{:.2f} {}"_f