CASTEP计算bader电荷的结果有点异常

minagami

从http://emuch.net/html/201411/8144608.html 《Bader 电荷分析，结果如何解读》一帖中我了解到原来能用CASTEP的输出结果计算bader电荷。计算的方法和代码主要参照了这个页面http://theory.cm.utexas.edu/henkelman/code/bader/ 《Code: Bader Charge Analysis》
计算步骤如下：
1.在.param文件添加WRITE_FORMATTED_DENSITY : TRUE，执行单点能计算，得到.den_fmt文件；

2.先执行gfortran den2cube.f90 -o den2cube，生成二进制执行文件den2cube。再把.den_fmt和.castep放在同一目录，执行den2cube seedname，把.den_fmt文件转变成.cube文件，这一步开始好像已经出问题了。.cube文件的前几行如下：
PdO
Extracted from a CASTEP density formatted file using den2cube
    4    0.000000    0.000000    0.000000
   31    0.190925    0.000000    0.000000
   31    0.000000    0.190925    0.000000
   55    0.000000    0.000000    0.186593
    8    6.000000    0.000000    0.000000    2.519005
    8    6.000000    0.000000    0.000000    7.557014
   46   10.000000    2.863880    0.000000    0.000000
   46   10.000000    0.000000    2.863880    5.038009
  0.20870E-03  0.22442E-03  0.25416E-03  0.32745E-03  0.44478E-03  0.66455E-03
  0.10105E-02  0.15474E-02  0.24120E-02  0.37542E-02  0.53537E-02  0.65546E-02
  0.68665E-02  0.67029E-02  0.67029E-02  0.68665E-02  0.65546E-02  0.53537E-02
  0.37542E-02  0.24120E-02  0.15474E-02  0.10105E-02  0.66455E-03  0.44478E-03
  0.32745E-03  0.25416E-03  0.22442E-03  0.20870E-03  0.22442E-03  0.25416E-03
和用另外一个工具check2xsf从.check文件提取出的电荷密度相比好像小了三个数量级。
在check2xsf里我使用的命令是check2xsf -c --cube PdO.check PdO.cube，得到的.cude文件的前几行如下:
DENSITY:

4 0.0 0.0 0.0
30 0.190925 0.000000 0.000000
30 0.000000 0.190925 0.000000
54 0.000000 0.000000 0.186593
8 0.0 0.000000 0.000000 2.519005
8 0.0 0.000000 0.000000 7.557015
46 0.0 2.863880 0.000000 0.000000
46 0.0 0.000000 2.863880 5.038010
0.207061
0.222656
0.252166
0.324876
下半部分的一列或几列数字对应的应该都是电荷密度呀，不知为何差异会这么大？

3. 执行bader PdO.cube，用cube文件得到bader电荷，在输出文件中ACF.dat的结果如下：
    #         X           Y           Z        CHARGE     MIN DIST    ATOMIC VOL
--------------------------------------------------------------------------------
    1   0.0000000   0.0000000   2.5190050   0.0394056   0.0000000  14.6373755
    2   0.0000000   0.0000000   7.5570140   0.0394056   0.0000000  14.6373755
    3   2.8638800   0.0000000   0.0000000   0.0477797   0.0000000  20.1604000
    4   0.0000000   2.8638800   5.0380090   0.0430494   0.0000000  18.4599614
--------------------------------------------------------------------------------
    VACUUM CHARGE:               0.0832
    VACUUM VOLUME:             291.6116
    NUMBER OF ELECTRONS:         0.2528
NUMBER OF ELECTRONS好像应该是整个体系的电子数，我算的体系是PdO晶胞，包含两个Pd原子和两个O原子，算出这么小的数值明显不合理呀？
而用另外一种方法得到的.cube文件算出的电荷数也有异常：
    #         X           Y           Z        CHARGE     MIN DIST    ATOMIC VOL
--------------------------------------------------------------------------------
    1   0.0000000   0.0000000   2.5190050  47.0103477   1.6792117  66.3783211
    2   0.0000000   0.0000000   7.5570150  47.0103477   1.6792110  66.3783211
    3   2.8638800   0.0000000   0.0000000  60.9626573   1.6017751  98.9043106
    4   0.0000000   2.8638800   5.0380100  60.9626573   1.6017744  98.9043106
--------------------------------------------------------------------------------
    VACUUM CHARGE:               0.0000
    VACUUM VOLUME:               0.0000
    NUMBER OF ELECTRONS:       215.9460
即使CHARGE一栏表示的是总电子数，O原子对应的总电子数也太多了吧。

不知道这里有没有人用过CASTEP来进行Bader电荷分析？欢迎各位一起来讨论~

zyj19831206

兄弟，你也是用CASTEP的，能加扣扣吗？

minagami

zyj19831206 发表于 2015-12-23 20:26
兄弟，你也是用CASTEP的，能加扣扣吗？

qq号已pm

zyj19831206

minagami 发表于 2015-12-23 20:28
qq号已pm

用castep搜索过渡态还是我的弱项，做了好几个反应，对于虚频的处理还是糊里糊涂。

minagami

zyj19831206 发表于 2015-12-23 20:30
用castep搜索过渡态还是我的弱项，做了好几个反应，对于虚频的处理还是糊里糊涂。

对于虚频我也没什么好办法，听说能够用Dmol3的EF法来消除，不过我实际上我也没有消除成功过

zyj19831206

minagami 发表于 2015-12-23 20:34
对于虚频我也没什么好办法，听说能够用Dmol3的EF法来消除，不过我实际上我也没有消除成功过

dmol3的EF方法不就是我发的那个帖子？就是寻找到合适的极小点，继续做ts search而已吧。

minagami

zyj19831206 发表于 2015-12-23 20:46
dmol3的EF方法不就是我发的那个帖子？就是寻找到合适的极小点，继续做ts search而已吧。

你的帖子《请问在dmol3中，这种消除虚频的办法是如何实施的？ 》里说的只是反应物和产物有虚频的情况吧。EF法即本征值跟踪法，对应的操作的Dmol3里的TS Optimization。

zyj19831206

minagami 发表于 2015-12-23 20:51
你的帖子《请问在dmol3中，这种消除虚频的办法是如何实施的？ 》里说的只是反应物和产物有虚频的情况吧。 ...

那个不怎么靠谱，有时候貌似把所有的虚频全消灭光了，我觉得还是要先ts conformation,然后寻找极小点，但是我不知道特定配位和xtd轨迹中结构之间的关系。

minagami

zyj19831206 发表于 2015-12-23 20:52
那个不怎么靠谱，有时候貌似把所有的虚频全消灭光了，我觉得还是要先ts conformation,然后寻找极小点，但 ...

你在能量走势图里选择能量最低点，然后xtd文件里就会出现对应的结构。接着复制粘贴到另一个空白文件里就行了。

zyj19831206

minagami 发表于 2015-12-23 20:55
你在能量走势图里选择能量最低点，然后xtd文件里就会出现对应的结构。接着复制粘贴到另一个空白文件里就 ...

当选择的那个点变成黄色就可以了？还是要show crosshair?

minagami

选择的那个点变成黄色显示选中就可以了，不用show crosshair

zyj19831206

minagami 发表于 2015-12-23 20:59
选择的那个点变成黄色显示选中就可以了，不用show crosshair

点变黄，xtd文件就自动变换了，这个还真没想过，有时候就卡在这种简单的近乎弱智的问题上，哎

minagami

zyj19831206 发表于 2015-12-23 21:03
点变黄，xtd文件就自动变换了，这个还真没想过，有时候就卡在这种简单的近乎弱智的问题上，哎

我也感觉经常会遇到类似的情况。唉，现实中没人指导，只能靠自己找资料自学和在网上多交流了。

卡开发发

检查一下，这两个程序可能其中有不能正确识别晶格的，另外这些程序可能没法正确对应CASTEP中超软赝势的芯电子态，需要对源码进行修改（当时只是无意间看了一下程序怎么count电子数的，发觉根本无法对应CASTEP中的价电子），这样再使用Bader程序分析肯定有问题。如果仅仅只是总电子数count错误导致，按说按照比例换算一下应该也成。

minagami

多谢你的回复！我还以为是我操作或者是设置上出了问题，没想到竟然是源代码本身有问题，如果是这样的话那就没办法了。我还不会读fortran的代码，如果不用超软而是用模守恒赝势的话不知道能不能行？