AIMD输出电荷信息

zhl@sy

大家好！ 我现在用AIMD计算Al/Fe2O3反应过程，想输出不同时间下体系原子的电荷变化，都有哪些方法可以实现呢？

卡开发发

若是Bader电荷我暂时还没啥办法，需要了解Bader程序本身以及电荷网格存储的子程序才行。但基于Wigner-Seiz半径积分的原子电荷还是可以简单实现的，可以修改主程序main.F(位置在7260 Format上面一行)，我手头的版本是5.3.2的：

1.         io_end                           ！ vasp原始程序
   
2. ！======插入如下代码
   
3.      IF (JOBPAR/=0 .OR. IO%LORBIT>=10  ) THEN
   
4.          
   
5.          IF (IO%LORBIT==11 .OR. IO%LORBIT==1 .OR. IO%LORBIT==12 .OR. IO%LORBIT==2) THEN
   
6.             LPAR=LMDIMP
   
7.          ELSE
   
8.             LPAR=LDIMP
   
9.          ENDIF
   
10. 
    
11.          DEALLOCATE(PAR)
    
12.          ALLOCATE(PAR(WDES%NB_TOT,WDES%NKPTS,LPAR,T_INFO%NIONP,WDES%NCDIJ))
    
13.          
    
14.          IF (IO%LORBIT>=10) THEN
    
15.             CALL SPHPRO_FAST( &
    
16.                  GRID,LATT_CUR, P,T_INFO,W, WDES, 71,IO%IU0,&
    
17.                  INFO%LOVERL,LMDIM,CQIJ, LPAR, LDIMP, LMDIMP, .TRUE., IO%LORBIT,PAR)
    
18.          ELSE
    
19.             CALL SPHPRO( &
    
20.                  GRID,LATT_CUR, P,T_INFO,W, WDES, 71,IO%IU0,&
    
21.                  INFO%LOVERL,LMDIM,CQIJ, LPAR, LDIMP, LMDIMP, LTRUNC, IO%LORBIT,PAR)
    
22.          ENDIF      
    
23. 
    
24.      ENDIF
    
25. ！======结束插入
    
26. 7260  FORMAT(/ &                  ！ vasp原始程序

复制代码

然后重新编译。设置完LORBIT=11之后每做完一个SCF loop都会把电荷打在屏幕上。有问题可以及时联系。

zhl@sy

卡开发发 发表于 2017-5-20 21:33
若是Bader电荷我暂时还没啥办法，需要了解Bader程序本身以及电荷网格存储的子程序才行。但基于Wigner-Seiz ...

我用的vasp版本是5.2.2，不知是否可以进行编译，另外想问一下这种电荷分析的方法与bader电荷分析有什么差异吗？该方法给出的电荷信息可信度怎样？谢谢！

卡开发发

zhl@sy 发表于 2017-5-21 13:03
我用的vasp版本是5.2.2，不知是否可以进行编译，另外想问一下这种电荷分析的方法与bader电荷分析有什么差 ...

我用的vasp版本是5.2.2，不知是否可以进行编译

可以，我实测的时候采用的就是5.2.x的。

这种电荷分析的方法与bader电荷分析有什么差异吗？该方法给出的电荷信息可信度怎样？

这些方法都有一定的任意性，比如WS半径选取的方式和Bader电荷中zero-flux这些对于某些特定体系都有可能违背化学常识而失效。
WS球的方法实际上用的并不算多，Bader的认可度高一些。但技术上有不少难点：
1、Bader电荷的精度依赖计算中电荷格点的密度，计算和调试的复杂度会大不少，因为实际处理总是得检查一下积分得到的总电子数目是否合理；
2、得专门独立写一个subroutine去代替原来vtst中提供的checksum.pl来合并价电子和芯电子电荷的格点，并且还得将Bader那个程序的代码都嵌入vasp内部，如果不是vasp内部人士，这样的改动愿意做的人应该不多；
3、如果不进行上述内嵌处理，就得储存每一步的价电子和芯电子的密度网格，之后也得写专门的脚本或程序来处理，如果要跑1000+帧的话，储存的代价实在太大了。

zhl@sy

卡开发发 发表于 2017-5-21 14:19
可以，我实测的时候采用的就是5.2.x的。

谢谢您的详细解答。我之前查的文献有给出bader电荷信息的，但vasp在给出反应过程中bader电荷似乎很困难，之前我看有用MS计算Milliken电荷的，难道在这方面vasp无能为力了？

卡开发发

zhl@sy 发表于 2017-5-21 16:30
谢谢您的详细解答。我之前查的文献有给出bader电荷信息的，但vasp在给出反应过程中bader电荷似乎很困难， ...

难道在这方面vasp无能为力了？

也不见得完全无能为力，比如你可以考虑用ASE做分子动力学的引擎调用VASP做SCF，也许对python代码做适当修改能够实现。

因为ASE本身只是调用VASP做单点能和力，这样可以在VASP的输入文件指定单点能后输出电荷密度，然后在ASE的分子动力学部分进行修改，要求每一步分子动力学在VASP做SCF后都对电荷密度做一次Bader分析，并储存到指定文件。

这样的修改工作量要比直接修改VASP本身小一些，当然运行的效率自然也会差一些。不过我对ASE和python不熟悉，你可能需要自己摸索一下。无论采用哪一种方法都对原理和编程的了解有一定的要求。

之前我看有用MS计算Milliken电荷的

没记错的话CASTEP或者DMol3应当也没有办法实现每一步离子步之后都做布居分析，加上源码封闭，连修改的机会都没有了。

zhl@sy

卡开发发 发表于 2017-5-21 19:06
也不见得完全无能为力，比如你可以考虑用ASE做分子动力学的引擎调用VASP做SCF，也许对python代码做适当 ...

谢谢您，学习了。ASE我是第一次听到，我再摸索一下。另外，我看到一篇文献中讲述有Milliken电荷计算的（见附件），但不知是否用CASTEP计算的。

卡开发发

zhl@sy 发表于 2017-5-21 20:47
谢谢您，学习了。ASE我是第一次听到，我再摸索一下。另外，我看到一篇文献中讲述有Milliken电 ...

ASE≡Atomic Simulation Environment

应该不是CASTEP，CASTEP并不支持共轭梯度极小化能量的方式做SCF。我觉得CMPD的可能性会比较大。但是否需要修改程序来完成每一步都做布居分析，这点我不清楚。

zhl@sy

卡开发发 发表于 2017-5-21 22:22
ASE≡Atomic Simulation Environment

应该不是CASTEP，CASTEP并不支持共轭梯度极小化能量的方式做SCF ...

是CPMD吧，Car-Parinello Molecular Dynamics，它相比BOMD谁在电荷及布局分析上更加方便有效一些呢？

卡开发发

zhl@sy 发表于 2017-5-22 08:54
是CPMD吧，Car-Parinello Molecular Dynamics，它相比BOMD谁在电荷及布局分析上更加方便有效一些呢？

不好意思，笔误，是CPMD，我指的也许是采用是叫CPMD的那个程序，方法还是BOMD。

zhl@sy

卡开发发 发表于 2017-5-22 10:57
不好意思，笔误，是CPMD，我指的也许是采用是叫CPMD的那个程序，方法还是BOMD。

老师您好！
之前问过您关于Bader电荷的问题。现在有个新问题向您请教：
我用AIMD计算了一个反应体系，设置的步数为1E4，由于文件较大所以没有设置电荷输出，但要求每200steps做一次SCF输出原子构型。由于我想分析该反应过程中原子之间电荷的转移规律，因此考虑计算bader电荷。采用的方法是从轨迹文件中导出特定步数下的原子轨迹文件，随后再去做SCF并计算bader电荷，并没有从AIMD的CHGCAR中提取分析。想问一下这种处理方式是否合适？是否能达到我想要的目的？
谢谢您！

卡开发发

zhl@sy 发表于 2017-9-10 23:20
老师您好！
之前问过您关于Bader电荷的问题。现在有个新问题向您请教：
我用AIMD计算了一个反应体系， ...

合适，而且控制起来更方便。动力学或者构型优化当中从第二步开始波函数和电荷密度就会采用上一步的外推，FFT网格也并不像CASTEP每一步都做修正，如果SCF收敛标准稍微不严格，以及晶格发生变化等，这些都会来带误差，如果拿这些结构出来做单点自洽反而结果可靠些。

建议最好写个脚本自动去提取结构生成输入文件，每个Bader电荷做完就把对应CHGCAR和CHG这些都删了，否则对硬盘是负担。

zhl@sy

卡开发发 发表于 2017-9-11 03:14
合适，而且控制起来更方便。动力学或者构型优化当中从第二步开始波函数和电荷密度就会采用上一步的外推， ...

好的，我利用这种方法计算了某一step下体系的Bader电荷，如附件图片所示。这个体系中间包含一个界面，其中红色的为Fe原子，我在赝势中采用的是Fed7s1，由于上下基体部分Fe原子没有反应，而只有中间部分有反应发生，因此其Bader电荷应在8.0附近，但从结果上来看Bader电荷比较分散，应该如何看待这种数据的分散性呢？谢谢老师！

卡开发发

zhl@sy 发表于 2017-9-11 09:07
好的，我利用这种方法计算了某一step下体系的Bader电荷，如附件图片所示。这个体系中间包含一个界面，其 ...

金属体系的电子本来就不是很局域，所以有一些扰动就会导致电荷看起来比较分散。但到底是不是Bader电荷反应出来的结果一定正确，这个我不好说。

zhl@sy

卡开发发 发表于 2017-9-11 13:31
金属体系的电子本来就不是很局域，所以有一些扰动就会导致电荷看起来比较分散。但到底是不是Bader电荷反 ...

谢谢老师！可能我没有说清楚，图中计算的结果是基于AIMD过程中某一时刻原子构型做电子的DFT得到的。此时Fe原子的分布比较混乱，没有了规则的原子结构排列，导致电荷的重新排布，使得Bader电荷显得较为分散，但整体上仍围绕8变化，是否可这样理解？另外，我想通过反应区域附近Fe原子Bader电荷的总体变化体现电荷的转移，利用bader是否能够实现呢？