振动分析Hessian矩阵的解读问题

AutMaple

各位大佬好，

我使用cp2k/2023-gcc-9.3.0进行振动分析，目的是为了获取hessian矩阵，将其使用不同的质量矩阵转化为对应的质权笛卡尔坐标后，获取同一结构下不同同位素组成的振动频率，进而获得我们同位素地球化学所需要的RPFR值（约化配分函数比值）。

目前，我已经成功的获得了CP2K输出的hess文件。然而，与VASP和GAUSSIAN以文本形式输出的hessian矩阵不同的是，CP2K输出的hessian矩阵保存在二进制文件中，我因此不能理解该文件中的每一个数据与hessian矩阵的矩阵元之间的对应关系。


我的计算体系是含有64个原子，因而其hessian矩阵应当是192*192的矩阵。目前，我讲hess文件用od命令查看，所得到的文件是一行9个元素。

希望各位大佬能帮忙提出解决思路！

ahxb

将输入文件中&GLOBAL字段中的PRINT_LEVEL设置为MEDIUM即可在输出文件中看到hessian矩阵的矩阵元。
如果不想重新计算的话，可以做个小体系的计算，其中PRINT_LEVEL设置为MEDIUM，然后自行用二进制查看软件（如Hex Editor等）比较.hess文件和输出文件中Hessian矩阵元的对应关系，尝试类推到你需要解析的hess文件，并求其本征值与本征矢量，与输出文件对比。
另外从数学上考虑的话，可以从本征值和本征矢量计算原来的的矩阵：若A和Q分别是B的本征值和本征矢量组成的矩阵，那么A=QBQ^-1，B=QAQ^-1。应用到你的问题中，需要首先把质权坐标下的向量转化到笛卡尔坐标下，组合成正交矩阵Q；将振动频率换算为本征值（注意单位），组成对角矩阵A；然后计算QAQ^-1即可。如果你愿意，可以这么干。

sobereva

直接从输出文件里读就完了，http://sobereva.com/soft/Sobtop对周期性体系计算键/键角/二面角力常数用的Hessian矩阵也是这么读的

AutMaple

ahxb 发表于 2024-5-8 23:48
将输入文件中&GLOBAL字段中的PRINT_LEVEL设置为MEDIUM即可在输出文件中看到hessian矩阵的矩阵元。
如果不 ...

感谢，我一直以为要从hess矩阵中读取

AutMaple

sobereva 发表于 2024-5-9 01:07
直接从输出文件里读就完了，http://sobereva.com/soft/Sobtop对周期性体系计算键/键角/二面角力常数用的Hes ...

感谢社长，我找到Hessian矩阵了。

我在处理CP2K的Hessian矩阵时，由于它输出的矩阵不是实对称矩阵，所以有时候python的numpy求解出来的特征值是复数。VASP输出的Hessian矩阵是实对称矩阵，不会出现这样的问题。

在谐振近似下，能量对坐标的二阶偏导如果是连续的，那么偏导次序是不会改变偏导结果的，那么CP2K输出的矩阵不为实对称矩阵是不是来源于CP2K在处理一阶偏导和二阶偏导的方法不同。CP2K分别使用了解析方法和有限差分方法处理一阶偏导和二阶偏导。我为了使得质权坐标系下的矩阵F特征值均为实数，是否可以将矩阵F和其转置相加后平均来转变为实对称矩阵再求特征值呢？

在此我想请教社长是怎么处理非实对称矩阵求解力常数的。

sobereva

AutMaple 发表于 2024-5-14 21:54
感谢社长，我找到Hessian矩阵了。

我在处理CP2K的Hessian矩阵时，由于它输出的矩阵不是实对称矩阵，所 ...

可以按你说的处理

AutMaple

sobereva 发表于 2024-5-14 23:07
可以按你说的处理

好的，社长，感谢您的解答。

在这里我提醒后面的各位同仁，CP2K在out文件中，（在天河二号上使用CP2K时使用yhrun -N 4 -n 96 cp2k.psmp CP2K-Freq.inp | tee freq.out  可以产生out文件）
”VIB| Hessian in cartesian coordinates“输出的为非质权坐标系下的Hessian矩阵，第二列为3N个元素标识，建立质权矩阵时以此为准，顺序不要弄错。

”VIB| Frequencies after removal of the rotations and translations“下面以”VIB| Internal  Low frequencies“开头，后面跟的数据是质权坐标系下Hessian矩阵的特征值，一共有3N-3个（FULLY_PERIODIC T），或者3N-6个（FULLY_PERIODIC F），N为体系原子个数，您可以比较自己解的特征值和他们解的特征值之间的差异。