计算化学公社

标题: LDA+U体系用dynaphopy计算声子谱求助 [打印本页]

作者
Author: 高阁 时间: 2025-2-18 16:06
标题: LDA+U体系用dynaphopy计算声子谱求助
本帖最后由高阁于 2025-2-18 22:38 编辑

背景
我尝试用dynaphopy算Mg7O8Fe在2000K，50Gpa下的非简谐声子谱。目前仅在Fe上加了磁矩，在Fe的d轨道上加了U值，计算按照如下步骤进行：
1.线性响应计算U值。初始磁矩设为4，根据Calculate U for LSDA+U - VASP Wiki中方法算U。从这一步开始每一步加压。
2.结构优化；
3.phonopy用有限位移法计算0K时声子谱与FORCE_CONSTANTS。此处单点能计算时根据上一步更新磁矩，扩胞为1*2*2，共64个原子，晶格参数为7.96978、7.92590、7.92590；
4.用NVP系综训练MFLL。此处POTIM=0.7，训练了50000步。从这一步开始每一步加温；
5.用MFLL计算MD，共200000步，POTIM=0.7；
6.Dynaphopy计算非谐性声子谱。
问题
1.请问如何确定哪些物质需要设置磁矩以及初始磁矩的值？目前我只给3d以及更高轨道的物质设置了磁矩，Fe的3d轨道有6个电子，4个轨道未充满，请问设置为4是否正确？
2.请问如何确定哪些物质需要设置U值？目前我只给Fe的d轨道加了U值，请问是否正确？
3.请问训练MFLL时IVDW应该如何设置？Part 2: Machine learning force fields4.2中给出IVDW=10，但是和LDA泛函冲突，目前我使用了PBE泛函，请问是否可行？如果用LDA泛函应该如何设置？
4.请问训练MFLL时能否使用单胞，等MD时再扩胞？上述扩胞大小是否合理？
5.请问用NVP系综训练MFLL时应该如何设置LANGEVIN_GAMMA_L和LANGEVIN_GAMMA？根据线性响应和U值以及涨落、耗散的关系，我猜测LANGEVIN_GAMMA_L为涨落，应该为LANGEVIN_GAMMA（耗散）的一半左右，LANGEVIN_GAMMA具体值应该与原子质量成正比，且应该与U值大小差不多。目前我将O的LANGEVIN_GAMMA取为U值的初值，其余元素按质量取初值，LANGEVIN_GAMMA_L取为U值作为初始值再调整。请问这样做是否合理？
6.请问用MFLL算MD后应该如何将轨迹可视化？目前使用ovito和VMD都会报错。
7.请问在训练MFLL和使用MFLL进行MD计算时有少量电子步未收敛（估计<3%）是否需要重新计算？
8.请问dynaphopy算声子谱时-psm 3和 -psm 1算得结果不同，是否说明参数设置还存在问题？

感谢指导！

作者
Author: 高阁 时间: 2025-2-18 22:12
XDATCAR文件超过9.8MB，无法上传。

作者
Author: 卡开发发 时间: 2025-2-18 23:23
我似乎是从别的帖子被叫到这个帖子来，不过很多问题我不了解，只能挑了解一点的说点自己看法。
1、自由原子大概是这样，不过Fe(III)可能会是5，但一般而言按4或5猜大多数情况都能得到合理的结果。
2、Mg设置U基本没啥用，Fe(3d)考虑Hubbard U肯定是合理的，一般争议点在于是否要对O(2p)考虑Hubbard U，这点暂无定论，不同文献可能看法都有所差异。
3、我不确定你和例子分别怎样考虑，首先LDA的认可度不高，其次是例子中的IVDW=10也较为过时，（个人经验和一些文献表明）DFT-D2确实通常会比DFT-D3结果可靠性差。理论上用PBE-D3(BJ)会更好，但我不确定在MLFF过程中有什么限制。
关于MLFF的使用你可能得看看有没有更了解的人给你答复。

作者
Author: 高阁 时间: 2025-2-18 23:25

卡开发发发表于 2025-2-18 23:23
我似乎是从别的帖子被叫到这个帖子来，不过很多问题我不了解，只能挑了解一点的说点自己看法。
1、自由原 ...

是我在另一个帖子下请您过来的，非常感谢老师的指导！

作者
Author: zhy996zhy 时间: 2025-3-26 23:36

卡开发发发表于 2025-2-18 23:23
我似乎是从别的帖子被叫到这个帖子来，不过很多问题我不了解，只能挑了解一点的说点自己看法。
1、自由原 ...

老师我有一个问题，目前在计算一个本征材料，他的原胞就很大，晶格常数分别是4-22-28，所以体系原子很多，这样的话计算声子谱还需要扩胞吗，如果扩胞的话可以只在x方向扩胞使得其大于10，尽可能减少计算量吗？

作者
Author: 卡开发发 时间: 2025-3-27 12:26

zhy996zhy 发表于 2025-3-26 23:36
老师我有一个问题，目前在计算一个本征材料，他的原胞就很大，晶格常数分别是4-22-28，所以体系原子很多 ...

经验上在4那个方向必须扩展2倍以上。实际严格的情况应该同样进行收敛性测试，只是代价太大了。

作者
Author: zhy996zhy 时间: 2025-3-27 13:11

卡开发发发表于 2025-3-27 12:26
经验上在4那个方向必须扩展2倍以上。实际严格的情况应该同样进行收敛性测试，只是代价太大了。

谢谢老师，也就是说起码X方向得扩一次，Y,Z可以不扩展吗？目前组里有超算，我担心的是vasp能不能算得动，单胞就有92个原子，扩展2倍就是184，理论上2倍可行吗

作者
Author: 卡开发发 时间: 2025-3-27 19:00

zhy996zhy 发表于 2025-3-27 13:11
谢谢老师，也就是说起码X方向得扩一次，Y,Z可以不扩展吗？目前组里有超算，我担心的是vasp能不能算得动， ...

X扩一次以上，YZ不扩展，这只是经验上这样说，我不得不说这样的做法没有什么严格的理论依据。如果实在很难算得动或许可以试试别的程序，即便对于LCAO程序在一些细节上的问题可能也并非设想那样好，只能说大体系的声子计算确实过于昂贵了。

作者
Author: zhy996zhy 时间: 2025-3-27 20:13

卡开发发发表于 2025-3-27 19:00
X扩一次以上，YZ不扩展，这只是经验上这样说，我不得不说这样的做法没有什么严格的理论依据。如果实在很 ...

好的非常感谢老师，我算一下试试

作者
Author: zhy996zhy 时间: 2025-3-27 22:54

卡开发发发表于 2025-3-27 19:00
X扩一次以上，YZ不扩展，这只是经验上这样说，我不得不说这样的做法没有什么严格的理论依据。如果实在很 ...

老师我试了x方向扩胞两倍，计算结果也收敛了，但是vasprun中没有出现二阶力常数矩阵，outcar也没有，这是为什么呢，用的DFPT方法
这是我的incar设置：# Basic param
ISTART=0
ICHARG=2
ISPIN=1
PREC = Accurate !precision , can set: Low | Medium | High | Normal | Single | Accurate

#Electronic Relaxation
ENCUT = 520 !Cutoff energy for the planewave basis set in eV
NELM = 80 !Maximum number of electronic SC (selfconsistency) steps which may be performed
NELMIN = 6 !Minimum number of electronic SCF steps.
NELMDL = -3
EDIFF = 1E-7 !The global break condition for the electronic SC-loop.
LREAL = False
LALGO=38
ADDGRID=True

# Ionic relaxation
EDIFFG = -1E-5
NSW = 1 !Maximum number of ionic steps.
IBRION = 8 !How the ions are updated and moved.2-GGA optimization method; -1-scf,band

ISMEAR = 0 !0-Gaussian smearing, -5-tetrahedron smearing
SIGMA = 0.01 !Width of smearing
ALGO = All !Electronic minimization method. Fast combines speed and robustness; VeryFast and Normal can also be used.

作者
Author: 卡开发发 时间: 2025-3-27 22:59

zhy996zhy 发表于 2025-3-27 22:54
老师我试了x方向扩胞两倍，计算结果也收敛了，但是vasprun中没有出现二阶力常数矩阵，outcar也没有，这是 ...

暂时不好说，你可以把输出文件传来看看。

作者
Author: zhy996zhy 时间: 2025-3-27 23:38
老师您需要哪个输出文件，这是四个输出文件

作者
Author: zhy996zhy 时间: 2025-3-27 23:38

卡开发发发表于 2025-3-27 22:59
暂时不好说，你可以把输出文件传来看看。

麻烦老师了，非常感谢

作者
Author: 卡开发发 时间: 2025-3-28 10:47

zhy996zhy 发表于 2025-3-27 23:38
麻烦老师了，非常感谢

你的OUTCAR当中显示的INCAR部分：

INCAR:
ISTART = 1
ICHARG = 2
ISPIN = 1
PREC = Accurate
ENCUT = 520
NELM = 80
NELMIN = 6
NELMDL = -3
EDIFF = 1E-7
LREAL = False
LALGO = 38
ADDGRID = True
EDIFFG = -1E-5
NSW = 0
IBRION = 8
ISMEAR = 0
SIGMA = 0.01
ALGO = All
注意，这里的NSW是0，你可能需要找找原因。

作者
Author: zhy996zhy 时间: 2025-3-28 11:12
本帖最后由 zhy996zhy 于 2025-3-28 11:14 编辑

卡开发发发表于 2025-3-28 10:47
你的OUTCAR当中显示的INCAR部分：

INCAR:

谢谢卡老师，我ALGO用的是all，我觉得可能和这个有关，也就是说他没有启动微扰吗，昨晚我重新跑了一个，用ALGO等于Normal，在log中多输出了以下内容：Linear response reoptimize wavefunctions to high precision
DAV: 1 -0.731358150770E+03 -0.22605E-06 -0.44333E-09  4320 0.844E-05
DAV: 2 -0.731358150771E+03 -0.29831E-09 -0.30888E-09  4320 0.502E-05
DAV: 3 -0.731358150771E+03 -0.50204E-09 -0.16288E-09  4320 0.268E-05
Linear response DOF=       72
Linear response progress:
  Degree of freedom: 1/ 72
   N    E                   dE          d eps    ncg    rms       rms(c)
RMM: 1 -0.248471924828E-01 -0.24847E-01 0.19780E-02 50396 0.157E-01
RMM: 2 -0.230871519360E-01 0.17600E-02 -0.33716E-03 36088 0.451E-02 0.706E-02
RMM: 3 -0.233513375550E-01 -0.26419E-03 -0.79871E-04 41575 0.155E-02 0.385E-02
RMM: 4 -0.236379932416E-01 -0.28666E-03 -0.15534E-04 42126 0.717E-03 0.200E-02
RMM: 5 -0.238251768152E-01 -0.18718E-03 -0.64102E-05 40426 0.493E-03 0.501E-03
RMM: 6 -0.238101256490E-01 0.15051E-04 -0.35625E-05 53008 0.147E-03 0.306E-03
RMM: 7 -0.237959064699E-01 0.14219E-04 -0.21891E-05 52194 0.939E-04 0.621E-04
RMM: 8 -0.237979485303E-01 -0.20421E-05 -0.21052E-05 59351 0.562E-04 0.262E-04
RMM: 9 -0.237976896918E-01 0.25884E-06 -0.20747E-05 63408 0.507E-04 0.865E-05
RMM:  10 -0.237974339776E-01 0.25571E-06 -0.18351E-05 66205 0.492E-04 0.358E-05
RMM:  11 -0.237974000625E-01 0.33915E-07 -0.18584E-05 67839 0.485E-04 0.146E-05
RMM:  12 -0.237973528784E-01 0.47184E-07 -0.16917E-05 68205 0.486E-04 0.125E-05
first order change of Fermi-energy -6.075729177101863E-003
这是不意味着有了微扰
另外这是还没跑完的outcar文件

作者
Author: 卡开发发 时间: 2025-3-28 11:54

zhy996zhy 发表于 2025-3-28 11:12
谢谢卡老师，我ALGO用的是all，我觉得可能和这个有关，也就是说他没有启动微扰吗，昨晚我重新跑了一个， ...

我不排除vasp内部可能ALGO等于特定值会强制NSW=0，我有怀疑过这件事，但搜了下源码中NSW和IBRION等被赋值0的情况没有找到。

当然按照道理说OUTCAR开头的INCAR信息应该和你输入是一致的，否则OUTCAT开头INCAR信息的IBRION=8按照正常逻辑也会被打印为IBRION=-1，但实际上这里并没有，而是在后续参数check的部分显示IBRION=-1，这里应该是在检查过NSW=0后被改为了IBRION=-1。

所以或许你应该再检查一下那个任务的INCAR，或者重新提交试试，具体机制我不明确。

作者
Author: zhy996zhy 时间: 2025-3-28 12:16

卡开发发发表于 2025-3-28 11:54
我不排除vasp内部可能ALGO等于特定值会强制NSW=0，我有怀疑过这件事，但搜了下源码中NSW和IBRION等被赋值 ...

我明白了老师，我再检查一下，重新跑那个OUTCAR 中IBRION是8，NSW是1，也有以下输出：Linear response reoptimize wavefunctions to high precision
DAV: 1 -0.731358150771E+03 -0.22685E-06 -0.44903E-09  4320 0.844E-05
DAV: 2 -0.731358150771E+03 0.36380E-09 -0.32204E-09  4320 0.502E-05
DAV: 3 -0.731358150770E+03 0.21100E-09 -0.18410E-09  4320 0.268E-05
Linear response DOF=       72
Linear response progress:
  Degree of freedom: 1/ 72
   N    E                   dE          d eps    ncg    rms       rms(c)
RMM: 1 -0.248471924827E-01 -0.24847E-01 0.19780E-02 50396 0.157E-01
RMM: 2 -0.230871519418E-01 0.17600E-02 -0.33716E-03 36088 0.451E-02 0.706E-02
RMM: 3 -0.233513401164E-01 -0.26419E-03 -0.79871E-04 41575 0.155E-02 0.385E-02
RMM: 4 -0.236379966169E-01 -0.28666E-03 -0.15534E-04 42126 0.717E-03 0.200E-02
RMM: 5 -0.238251825231E-01 -0.18719E-03 -0.64102E-05 40427 0.493E-03 0.501E-03
RMM: 6 -0.238101303555E-01 0.15052E-04 -0.35627E-05 53008 0.147E-03 0.306E-03
RMM: 7 -0.237959157941E-01 0.14215E-04 -0.21905E-05 52190 0.939E-04 0.621E-04
RMM: 8 -0.237979507496E-01 -0.20350E-05 -0.21008E-05 59352 0.562E-04 0.262E-04
RMM: 9 -0.237976832849E-01 0.26746E-06 -0.20606E-05 63412 0.507E-04 0.865E-05
RMM:  10 -0.237974415404E-01 0.24174E-06 -0.18336E-05 66204 0.492E-04 0.358E-05
RMM:  11 -0.237974015809E-01 0.39959E-07 -0.18550E-05 67839 0.485E-04 0.146E-05
RMM:  12 -0.237973589879E-01 0.42593E-07 -0.16976E-05 68205 0.486E-04 0.125E-05
first order change of Fermi-energy -6.076101350771523E-003

作者
Author: 卡开发发 时间: 2025-3-28 12:59

zhy996zhy 发表于 2025-3-28 12:16
我明白了老师，我再检查一下，重新跑那个OUTCAR 中IBRION是8，NSW是1，也有以下输出：Linear response re ...

是啊，要开始声子计算NSW至少是大于0的值，否则IBRION会被自动改成-1。

作者
Author: zhy996zhy 时间: 2025-3-28 14:41

卡开发发发表于 2025-3-28 12:59
是啊，要开始声子计算NSW至少是大于0的值，否则IBRION会被自动改成-1。

非常感谢卡开老师，希望这次计算顺利吧

作者
Author: zhy996zhy 时间: 2025-4-3 19:04

卡开发发发表于 2025-3-28 12:59
是啊，要开始声子计算NSW至少是大于0的值，否则IBRION会被自动改成-1。

卡开老师，抱歉想在问个问题，这是我的outcar文件声子谱计算结果有三个虚频，以及绘制的声子谱，请问这种大小的虚频有可能消除吗？如果要消除，可以往那个方向努力呢？

作者
Author: zhy996zhy 时间: 2025-4-3 19:44
卡开老师，这是高对称点：High-symmetry points (in fractional coordinates).
0.0000000000 0.0000000000 0.0000000000    GAMMA
0.0000000000 0.5000000000 0.0000000000    Z
0.0000000000 0.0000000000 0.5000000000    B
0.0000000000 0.0000000000  -0.5000000000    B_2
0.5000000000 0.0000000000 0.0000000000    Y
  -0.5000000000 0.0000000000 0.0000000000    Y_2
0.5000000000 0.5000000000 0.0000000000    C
  -0.5000000000 0.5000000000 0.0000000000    C_2
0.0000000000 0.5000000000 0.5000000000    D
0.0000000000 0.5000000000  -0.5000000000    D_2
  -0.5000000000 0.0000000000 0.5000000000    A
  -0.5000000000 0.5000000000 0.5000000000    E
  -0.5000000000 0.0000000000 0.5000000000    H
  -0.5000000000 0.0000000000 0.5000000000    H_2
  -0.5000000000 0.0000000000  -0.5000000000    H_4
  -0.5000000000 0.5000000000 0.5000000000    M
  -0.5000000000 0.5000000000 0.5000000000    M_2
  -0.5000000000 0.5000000000  -0.5000000000    M_4

作者
Author: zhy996zhy 时间: 2025-4-3 21:35

zhy996zhy 发表于 2025-3-28 14:41
非常感谢卡开老师，希望这次计算顺利吧

卡开老师，我重新确定了两个虚频最低点对应坐标，一个在G点和Z点之间，另一个在高对称点D上

作者
Author: 卡开发发 时间: 2025-4-3 22:44

zhy996zhy 发表于 2025-4-3 21:35
卡开老师，我重新确定了两个虚频最低点对应坐标，一个在G点和Z点之间，另一个在高对称点D上

那你可能得考虑对精度和超胞尺寸测试一下了。

作者
Author: zhy996zhy 时间: 2025-4-4 08:58

卡开发发发表于 2025-4-3 22:44
那你可能得考虑对精度和超胞尺寸测试一下了。

卡开老师，您觉得这两个虚频问题大吗，是否可以忽略，审稿人会不会一直追问

作者
Author: 卡开发发 时间: 2025-4-4 09:21

zhy996zhy 发表于 2025-4-4 08:58
卡开老师，您觉得这两个虚频问题大吗，是否可以忽略，审稿人会不会一直追问

按道理不太会，其实按道理这个虚频足够小了。

作者
Author: zhy996zhy 时间: 2025-4-4 11:11

卡开发发发表于 2025-4-4 09:21
按道理不太会，其实按道理这个虚频足够小了。

谢谢卡开老师，我重新看了下，outcar 没给全，最大虚频在-6cm-1，如果是这个量级呢

作者
Author: 卡开发发 时间: 2025-4-5 04:25

zhy996zhy 发表于 2025-4-4 11:11
谢谢卡开老师，我重新看了下，outcar 没给全，最大虚频在-6cm-1，如果是这个量级呢

有多少虚频呢？

作者
Author: zhy996zhy 时间: 2025-4-5 09:58
数量吗，还是挺多的大概有 50 多条，就在图片中那两个下沉区域

作者
Author: zhy996zhy 时间: 2025-4-6 15:13

卡开发发发表于 2025-4-5 04:25
有多少虚频呢？

就是利用phonopy得出的band.dat绘图后有最大-0.21THZ的虚频，outcar中只有那三个，我不太清楚是否有可能是phonopy本身插值的问题

作者
Author: 卡开发发 时间: 2025-4-6 18:47

zhy996zhy 发表于 2025-4-6 15:13
就是利用phonopy得出的band.dat绘图后有最大-0.21THZ的虚频，outcar中只有那三个，我不太清楚是否有可能 ...

确实有可能，但也不排除扩更大胞能够缓解。

作者
Author: zhy996zhy 时间: 2025-4-6 20:30

卡开发发发表于 2025-4-6 18:47
确实有可能，但也不排除扩更大胞能够缓解。

好的谢谢老师，唉主要放图上不好看，只看outcar理论上够小了

作者
Author: zhy996zhy 时间: 2025-4-13 13:12

卡开发发发表于 2025-4-6 18:47
确实有可能，但也不排除扩更大胞能够缓解。

卡开老师，我重新优化结构，力收敛到-7，计算声子谱putcar中结果如下，
550 f/i= 0.000188 THz    0.001183 2PiTHz 0.006280 cm-1    0.000779 meV
551 f/i= 0.000510 THz    0.003202 2PiTHz 0.016998 cm-1    0.002107 meV
552 f/i= 0.000652 THz    0.004095 2PiTHz 0.021741 cm-1    0.002696 meV，
想询问下您觉得这个虚频够小可以忽略吗？

作者
Author: 卡开发发 时间: 2025-4-14 08:10
本帖最后由卡开发发于 2025-4-14 08:14 编辑

zhy996zhy 发表于 2025-4-13 13:12
卡开老师，我重新优化结构，力收敛到-7，计算声子谱putcar中结果如下，
550 f/i= 0.000188 THz 0 ...

原则上说，这么大的虚频肯定是可忽略的，剩下的看文章怎么写和图怎么出。但先不要急，如果这么看，一个是之前的收敛精度可能并不高，而且你之前可能没有重新优化过结构，建议多优化几次来消除一下Pulay Stress，再去算声子比较好，我没办法了解到你全部的计算流程。

作者
Author: zhy996zhy 时间: 2025-4-14 21:16

卡开发发发表于 2025-4-14 08:10
原则上说，这么大的虚频肯定是可忽略的，剩下的看文章怎么写和图怎么出。但先不要急，如果这么看，一个是 ...

好的谢谢卡开老师，我计划先再尝试下波恩修正，这个是重新优化过结构算的声子谱，力收敛到10-7之后一直波动，感觉很难下降了，您说的方式我也尝试一下

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