计算化学公社

标题: Dynaphopy 计算高温声子谱 [打印本页]

作者
Author: Shana 时间: 2023-5-9 10:15
标题: Dynaphopy 计算高温声子谱
小弟最近在计算钙钛矿体系的声子谱，在简谐近似下，计算出的声子谱都有一定的虚频，始终无法消除，一开始花了大量的时间去提高结构优化精度，但是都以失败告终，一度都怀疑自己的结构优化靠谱性。后来从各种文献、论坛和视频中了解到钙钛矿体系有很强的非谐性，简谐近似下的声子谱往往都会带有虚频，需要考虑非简谐情形下的计算。目前有三种非简谐计算方法，可以参考这个专栏300 K下的声子谱怎么算？ - 知乎 (zhihu.com)。小弟看了教程视频之后，了解到dynaphopy计算高温声子谱的效率比较高，和vasp也十分兼容，就尝试了一下。下面是小弟的计算历程。
1. 简谐声子谱的计算。
dynaphopy的非谐修正是在phonopy的力常数基础上进行的，需要先进行0k声子谱的计算。
这是结构优化的INCAR，力的收敛标准是E-3，大部分的文献也是在这个精度范围内。结构优化的结构是单胞，10个原子。
这是phonopy计算出的声子谱，计算的结构为2*2*2的超胞，80个原子，在gmma点和X点有很强的虚频。这个虚频即使提高收敛精度和扩大超胞也是无法消除的。
2. AIMD计算
dynaphopy通过读取AIMD的计算结构进行拟合（说法不准确），需要预先准备AIMD的outcar或xdatcar。
这是小弟的AIMD的计算设置，基本参考dynaphopy的example中的设置，详细的参数涵义，都有教程介绍。
小弟这里只是采用了3000步的计算，步长是2fs，能量的收敛精度是E-6，不算太高。同样是采用80个原子的超胞。大部分做非谐计算的文献，所采用的都是在80-160原子数范围内的超胞，步长多是1-2fs，总步数大概是3000-6000，总的时长是5ps-10ps。从下面的非谐声子谱结果来看，这个设置对于非谐计算是合适的？
3. dynaphopy非谐声子谱计算。
dynaphopy的官网如下：Dynaphopy (abelcarreras.github.io)。
将结构优化的原始晶胞POSCAR，phonopy输出的力常数，和AIMD的OUTCAR放在一个文件夹内。
设置dynaphopy的输入文件input如下：分别是输入结构文件，力常数文件，原胞和超胞与POSCAR的关系设置，band是声子谱的高对称点路径，设置方式和vasp计算PBE能带一致。
然后输入dynaphopy input OUTCAR -i，进入如下的界面，选择6，进入非谐声子谱计算，再选择1，可以同时输出简谐声子谱和非简谐声子谱，此时会对peak进行拟合。拟合结束后，会输出声子谱。
这是小弟一开始使用默认设置，输出的声子谱，可以看到高温声子谱虚频加重了。从@get-it 大佬处了解到，dynaphopy默认使用最大熵方法拟合，效果可能不如fftw有效。于是小弟将拟合方式改为fftw,通过指令dynaphopy input OUTCAR -i -psm 3,采用快速傅里叶变换拟合。得到300k下的声子谱，可以看到，原来gamma点和X点的虚频，已经被成功消除了。而实验上，这个材料在室温下肯定是稳定的。

当然，如果想要得到更准确的声子谱，需要提高优化精度，采用更大的超胞计算力常数，AIMD需要更多的步数和更小的步长，比如dynaphopy的例子是采用0.7 fs,跑200000步。如果只是想得到一个比较合理的声子谱，小弟觉得适当降低精度是没问题的。比如小弟的这个计算，在24核的机器上，结构优化用时不到3小时，dfpt计算1小时，AIMD计算14小时，差不多用时一天，能够得到一个比较合适的声子谱。如果提高精度，计算用时可能会提升数倍，但是声子谱的提升可能并不会太明显。当然这是小弟的一己之见，供大家批判

作者
Author: waitingseven 时间: 2023-5-9 11:21
请问这个对AIMD采用的系综有要求吗，NVT或者NPT

作者
Author: leeru 时间: 2023-5-9 14:27
采用alamode也可以达到相同的效果

作者
Author: Shana 时间: 2023-5-9 16:24

waitingseven 发表于 2023-5-9 11:21
请问这个对AIMD采用的系综有要求吗，NVT或者NPT

我用的是NVT，没有试过NPT哎

作者
Author: Shana 时间: 2023-5-9 16:24

leeru 发表于 2023-5-9 14:27
采用alamode也可以达到相同的效果

alamode有点复杂，没有dynaphopy这么简单

作者
Author: sobereva 时间: 2023-5-9 23:51
CUTOFF足够大是消除虚频的极其关键要素，绝不仅仅是看优化收敛精度

谐振近似得到的频率若有明显虚频总有办法消除（在静态计算框架下，与MD和温度无关），这是原理上注定的，否则相应物质根本无法稳定存在，这和非谐振效应强弱是两码事。PS：不过即便计算精度非常高，gamma点的acoustic模式的频率（原理为0）也可能为非常轻微的虚频，可以当做0看待。

作者
Author: Shana 时间: 2023-5-10 12:53

sobereva 发表于 2023-5-9 23:51
CUTOFF足够大是消除虚频的极其关键要素，绝不仅仅是看优化收敛精度

谐振近似得到的频率若有明显虚频总有 ...

sob 老师，您说的消除虚频方法，是不是改变原子位置，然后继续优化？但是这种改变原子位置的方法，往往伴随着对称性的变化，我认为这样即使没有虚频，也是在另外一种对称性没有虚频。sob老师，我这样理解是否合适？

作者
Author: sobereva 时间: 2023-5-10 22:20

Shana 发表于 2023-5-10 12:53
sob 老师，您说的消除虚频方法，是不是改变原子位置，然后继续优化？但是这种改变原子位置的方法，往往伴 ...

我这里是强调平面波截断能设置必须足够高，没有说改变原子位置的事。如果这个不够高，体系能量的平移不变性会较差，由此会在有限差分计算时引入虚假的能量改变并很可能造成虚频。

作者
Author: Shana 时间: 2023-5-11 10:24

sobereva 发表于 2023-5-10 22:20
我这里是强调平面波截断能设置必须足够高，没有说改变原子位置的事。如果这个不够高，体系能量的平移不变 ...

又从sob老师这学到了一点，我的截断能都是2倍推荐的截断能，这个应该是足够了吧?

作者
Author: sobereva 时间: 2023-5-11 12:36

Shana 发表于 2023-5-11 10:24
又从sob老师这学到了一点，我的截断能都是2倍推荐的截断能，这个应该是足够了吧?

遇到消灭不了的虚频，考虑尝试更高的

作者
Author: Shana 时间: 2023-5-12 09:46

sobereva 发表于 2023-5-11 12:36
遇到消灭不了的虚频，考虑尝试更高的

好的，谢谢sob老师

作者
Author: 乐平 时间: 2023-7-12 13:59
您好，您在帖子里提到

目前有三种非简谐计算方法，可以参考这个专栏300 K下的声子谱怎么算？ - 知乎 (zhihu.com)。小弟看了教程视频之后，了解到dynaphopy计算高温声子谱的效率比较高，

请问哪里有dynaphopy教程视频？谢谢！

作者
Author: Shana 时间: 2023-7-16 15:32

乐平发表于 2023-7-12 13:59
您好，您在帖子里提到

请问哪里有dynaphopy教程视频？谢谢！

官方网站学的，http://abelcarreras.github.io/DynaPhoPy/

作者
Author: 乐平 时间: 2023-7-17 10:02

Shana 发表于 2023-7-16 15:32
官方网站学的，http://abelcarreras.github.io/DynaPhoPy/

谢谢，但是这不是您说的“视频”……

另外，我理解是是，还是需要先跑 phonopy 结合 VASP 计算谐性声子谱，然后再跑 AIMD，用 DynaPhoPy从轨迹中提取信息得到非谐性声子谱。是这样的流程吧？

作者
Author: Shana 时间: 2023-7-17 20:13

乐平发表于 2023-7-17 10:02
谢谢，但是这不是您说的“视频”……

另外，我理解是是，还是需要先跑 phonopy 结合 VASP 计算谐性声 ...

是的，这个软件很小众，基本没有视频。

作者
Author: 李江山 时间: 2023-11-1 10:09
你好，请问我得到了quasiparticles_data.yaml文件如何画出声子频率，因为我到了画图那一步报错了，只生成了这个文件，可以手动作图吗

作者
Author: Iridium-191 时间: 2024-1-19 15:10

sobereva 发表于 2023-5-9 23:51
CUTOFF足够大是消除虚频的极其关键要素，绝不仅仅是看优化收敛精度

谐振近似得到的频率若有明显虚频总有 ...

老师好，那对于比如铁电体或者STO这种量子顺电体来说，或者其他高温相，室温相在0K下都是非稳态的相，这种是不是只能通过有限温度下的计算来解决了？谢谢老师

作者
Author: 小付小氟 时间: 2024-4-8 16:08
将拟合方式改为fftw,通过指令dynaphopy input OUTCAR -i -psm 3,采用快速傅里叶变换拟合。我想问用端口之间输出了图片。输入什么样的指令导出图片对应的数值呢？使用这种方法拟合，是对应的力常数文件发生改变吗

作者
Author: Shana 时间: 2024-4-17 09:42

小付小氟发表于 2024-4-8 16:08
将拟合方式改为fftw,通过指令dynaphopy input OUTCAR -i -psm 3,采用快速傅里叶变换拟合。我想问用端口之间 ...

能够直接输出力常数的

作者
Author: 九点一定起 时间: 2024-6-12 11:03
您好，我想问下那个input文件中的晶胞与原胞的转换矩阵是怎么得到的呀，如何得到这个矩阵呢

作者
Author: Shana 时间: 2024-6-12 15:14

九点一定起发表于 2024-6-12 11:03
您好，我想问下那个input文件中的晶胞与原胞的转换矩阵是怎么得到的呀，如何得到这个矩阵呢

你可以将晶胞转换成原胞，然后用原胞扩胞算声子和AIMD。这样转换矩阵就是对角矩阵

作者
Author: 九点一定起 时间: 2024-6-18 11:11

Shana 发表于 2024-6-12 15:14
你可以将晶胞转换成原胞，然后用原胞扩胞算声子和AIMD。这样转换矩阵就是对角矩阵

您好我这个声子和AIMD都算完了但是这个矩阵具体是怎么来的呢在哪里看呢

作者
Author: 九点一定起 时间: 2024-6-18 12:37

九点一定起发表于 2024-6-18 11:11
您好我这个声子和AIMD都算完了但是这个矩阵具体是怎么来的呢在哪里看呢

就是dynaphopy 中input中需要输入的转换矩阵这些我都算完了不知道这个怎么设置是怎么算出来的还是怎么得到的呢麻烦您了

作者
Author: Shana 时间: 2024-6-28 10:42

九点一定起发表于 2024-6-18 12:37
就是dynaphopy 中input中需要输入的转换矩阵这些我都算完了不知道这个怎么设置是怎么算出来的还是怎么 ...

有个PRIMITIVE MATRIX和SUPERCELL MATRIX，比如你将晶胞先约化成原胞，那么PRIMITIVE MATRIX的对角元就是1，1，1.然后2*2*2扩胞，SUPERCELL MATRIX的对角元就是2，2，2
PRIMITIVE MATRIX
1 0 0
0 1 0
0 0 1
SUPERCELL MATRIX
2 0 0
0 2 0
0 0 2

作者
Author: Shana 时间: 2024-6-28 10:43

九点一定起发表于 2024-6-18 12:37
就是dynaphopy 中input中需要输入的转换矩阵这些我都算完了不知道这个怎么设置是怎么算出来的还是怎么 ...

phonopy是可以输出PRIMITIVE MATRIX的，可以参考输出界面

作者
Author: 17873772769 时间: 2024-10-21 17:07
Shana老师，您好，在dynaphopy的参数设置中写了使用fftw方法时，psm可设置为2、3、4，我在使用dynaphopy input OUTCAR -i -psm 3时，会存在报错“ModuleNotFoundError: No module named 'pyfftw'”，但我使用2时可以正常运行，这结果应该是一样的么（使用后确实无虚频了，且看起来结果合理）

作者
Author: Shana 时间: 2024-10-28 16:27

17873772769 发表于 2024-10-21 17:07
Shana老师，您好，在dynaphopy的参数设置中写了使用fftw方法时，psm可设置为2、3、4，我在使用dynaphopy in ...

从我的使用感受上来看，dynaphopy的几个算法处理的结果是不大一样的，无虚频的话，消除的很成功

作者
Author: etaipxe 时间: 2024-12-26 19:19
您好，请问0K下声子谱没有虚频，但是300K反而有虚频了，这样正常吗？有什么解决方法呢

作者
Author: Shana 时间: 2024-12-27 16:16

etaipxe 发表于 2024-12-26 19:19
您好，请问0K下声子谱没有虚频，但是300K反而有虚频了，这样正常吗？有什么解决方法呢

这个是有可能的，能找到很多参考文献的，请搜索anharmonicity/非谐性。这个可能是材料本身的特性，不一定是错误的结果。

作者
Author: mly 时间: 2025-1-24 23:14
想请教您一下，我想计算300K下的吉布斯自由能，我需要在AIMD的INCAR中设置300K，在使用phonopy提取热力学性质时，也输入300K就可以吗？我看phonopy也可以计算300K-1000K（或者是其他温度）这俩有什么区别吗？感谢

作者
Author: Shana 时间: 2025-2-3 21:12

mly 发表于 2025-1-24 23:14
想请教您一下，我想计算300K下的吉布斯自由能，我需要在AIMD的INCAR中设置300K，在使用phonopy提取热力学性 ...

phonopy的300k自由能可以通过自由能的定义去算，即考虑T*S部分，S 是熵，和声子有关，phonopy可以根据声子谱计算出熵。当然，考虑到温度对声子的影响，phonopy支持QHA方法来计算有限温的自由能。建议看看phonopy的例子部分https://phonopy.github.io/phonopy/examples.html，以及理论部分https://phonopy.github.io/phonop ... perties-expressions。

作者
Author: mly 时间: 2025-2-14 09:46

Shana 发表于 2025-2-3 21:12
phonopy的300k自由能可以通过自由能的定义去算，即考虑T*S部分，S 是熵，和声子有关，phonopy可以根据声 ...

感谢

作者
Author: 高阁 时间: 2025-2-17 17:48
请问非谐性声子谱的数据怎么像phonopy一样导出呢？我需要用它们计算后续内容，感谢！还有我觉得这个声子谱似乎有些奇怪，请问它正常吗？这是2000k下的

作者
Author: Shana 时间: 2025-2-22 12:02

高阁发表于 2025-2-17 17:48
请问非谐性声子谱的数据怎么像phonopy一样导出呢？我需要用它们计算后续内容，感谢！还有我觉得这个声子谱 ...

dynaphopy input OUTCAR -sfc FORCE_CONSTANTS 这个指令会给出修正后的力常数，然后用phonopy画出修正后的声子谱

作者
Author: 高阁 时间: 2025-2-23 11:09

Shana 发表于 2025-2-22 12:02
dynaphopy input OUTCAR -sfc FORCE_CONSTANTS 这个指令会给出修正后的力常数，然后用phonopy画出修正 ...

谢谢

作者
Author: yuanzhe 时间: 2025-2-24 21:11
本帖最后由 yuanzhe 于 2025-2-25 14:56 编辑

shana老师，用vasp跑aimd需要用phonopy扩胞的POSCAR进行计算吗？kpoints是直接vaspkit还是手动设置吗？

作者
Author: Shana 时间: 2025-2-25 17:28

yuanzhe 发表于 2025-2-24 21:11
shana老师，用vasp跑aimd需要用phonopy扩胞的POSCAR进行计算吗？kpoints是直接vaspkit还是手动设置吗？

扩胞系数要和phonopy的系数一样，不然结果会有问题。晶胞比较大的话（扩胞系数比较大），用单k点就很好了。

作者
Author: yuanzhe 时间: 2025-3-3 15:47

Shana 发表于 2025-2-25 17:28
扩胞系数要和phonopy的系数一样，不然结果会有问题。晶胞比较大的话（扩胞系数比较大），用单k点就很好了 ...

谢谢shana老师解惑，我一般扩2*2*2，用vaspkit生成的就可以了嘛

作者
Author: Shana 时间: 2025-3-5 15:22

yuanzhe 发表于 2025-3-3 15:47
谢谢shana老师解惑，我一般扩2*2*2，用vaspkit生成的就可以了嘛

可以，phonopy也可以扩胞。phonopy扩完胞后，会有一个SPOSCAR，这个是未做位移的supercell。

作者
Author: yuanzhe 时间: 2025-3-11 15:19

Shana 发表于 2025-3-5 15:22
可以，phonopy也可以扩胞。phonopy扩完胞后，会有一个SPOSCAR，这个是未做位移的supercell。

谢谢shana老师！！！

作者
Author: Autumn899 时间: 2025-4-8 22:56
请问什么命令可以导出能在origin里画图的声子数据文件呀

作者
Author: Shana 时间: 2025-4-11 09:32

Autumn899 发表于 2025-4-8 22:56
请问什么命令可以导出能在origin里画图的声子数据文件呀

先保存重整后的力常数，然后用phonopy计算声子谱，导出dat数据

作者
Author: Autumn899 时间: 2025-4-11 17:01

Shana 发表于 2025-4-11 09:32
先保存重整后的力常数，然后用phonopy计算声子谱，导出dat数据

好的！

作者
Author: cdx 时间: 2025-6-16 08:59
本帖最后由 cdx 于 2025-6-16 09:00 编辑

老师请问利用Dynaphopy可以得到声子寿命吗？已完成非谐计算，需要如何导出声子寿命数据，求指导

作者
Author: cdx 时间: 2025-6-16 15:03

Shana 发表于 2025-2-22 12:02
dynaphopy input OUTCAR -sfc FORCE_CONSTANTS 这个指令会给出修正后的力常数，然后用phonopy画出修正 ...

老师你好，请问可以正常获得非谐修正后的力常数，可以利用phonopy导出声子寿命数据吗

作者
Author: Shana 时间: 2025-6-17 17:36

cdx 发表于 2025-6-16 15:03
老师你好，请问可以正常获得非谐修正后的力常数，可以利用phonopy导出声子寿命数据吗

这个我不知道哎。没用过

作者
Author: cdx 时间: 2025-6-27 22:12

小付小氟发表于 2024-4-8 16:08
**** 作者被禁止或删除内容自动屏蔽 ****

你好，请问您有找到导出图片中数据的方法吗

欢迎光临计算化学公社 (http://ccc.keinsci.com/)