siesta赝势生成

作者: 多米尼克2049 | 来源:发表于2020-03-10 15:00 被阅读0次

首先我们要明确为什么要生成赝势。

img122.gif
首先看这个图，回忆我们计算电子波函数的时候，哈密顿量里面的势能项分为库伦势项、Hatree项和交换关联能项。如果我们直接使用库伦势的话，也就是对应图中的Z/r势，由于这个势很深，所以波函数里会包含很多内层电子的成分，这些电子速度通常很高，对应的截断能就非常高，然后为了保证互相之间的正交性，波函数会有很多震荡，与x轴的交点叫做波节。如果我们把内层电子加原子核看成一个离子，也就是在计算的时候用一个‘赝势’

V_{pseudo}

来代替真实的势能函数，相当于不考虑内层电子对计算结果的影响，这样波函数就只包含价电子，全电子波函数就变成没有波节的‘赝波函数’。就可以大大减小计算量。这里的分界线是自己选的一个

r_c

,所以生成赝势之后需要check有没有产生非物理的态，就是所谓的Ghost State。
生成赝势的时候首先解一个薛定谔方程

全电子波函数的角标

l

就是轨道角动量，对应的全电子波函数就是球贝塞尔函数的叠加，相当于以轨道波函数作为基矢。

然后通过前面的系数来反过来拟合势函数V中的几个参数，从而得到赝势。
决定这些参数的时候往往有几个限制条件，比如电荷密度相等，这样可以保证电子的散射性质相同。
赝波函数必须连续可导
本征值必须相同。
激发态波函数也要包含进去，当使用赝势来处理激发态问题的时候。
这一部分来自The Pseudopotential Approximation

记下来再次遇到比较方便看。
新版siesta不再包含atom需要自己编译，这里记录一下。
首先在siesta根目录新建文件夹

mkdir atom

然后下载atom解压

cd atom
wget https://departments.icmab.es/leem/SIESTA_MATERIAL/Pseudos/Code/atom-4.2.7-100.tgz
tar -xvzf atom-4.2.7-100.tgz

因为atom要依赖两个库libgridxc和xmlf90，所以还需要下载这两个库编译

wget https://launchpad.net/libgridxc/trunk/0.8/+download/libgridxc-0.8.0.tgz
wget https://launchpad.net/xmlf90/trunk/1.5/+download/xmlf90-1.5.3.tar.gz
tar -xvzf xmlf90-1.5.3.tar.gz
tar -xvzf libgridxc-0.8.0.tgz

解压以后查看INSTALL文件，里面写的有编译测试方式。
先编译xmlf90，这个比较简单。
首先因为我们是在服务器上编译所以没法安装到/local/bin文件夹下面，所以我们需要提前指定安装位置，先创建一个文件夹

cd xmlf90-1.5.3/
mkdir build
./configure --prefix=/public/home/xbliu/apps/siesta/siesta-4.0.2/ATOM/xmlf90-1.5.3/build
make
make install

xmlf90就编译完了，然后是libgridxc，这个稍微麻烦一点，需要改一下fortran.mk

cd libgridxc-0.8.0/
mkdir build
cp extra/fortran.mk build/
vi fortran.mk
加一句
FC=gfortran
然后
../src/config.sh
make
make install
然后到Testers文件夹下面测试一下
cd Testers
make
会生成4个可执行文件，在当前文件夹下面执行一下看结果和Reference里面是否相同。
./test1
./test2
./test3
./test4

然后就可以编译atom了，因为atom是依赖于上面两个库的，所以需要在atom的arch.make文件里指定上面两个库的位置，这样它才能找到那两个库。

vi arch.make.sample
XMLF90_ROOT=/public/home/xbliu/apps/siesta/siesta-4.0.2/ATOM/xmlf90-1.5.3/build
GRIDXC_ROOT=/public/home/xbliu/apps/siesta/siesta-4.0.2/ATOM/libgridxc-0.8.0/build
cp arch.make.sample arch.make
make

到这里编译的工作就完成了，下面是赝势的生成。
直接进入Tutorial文件夹，里面有三个文件夹All_electron,PS_generation和Utils
All_electron里面是全电子的计算，PS_generation是一些示例输入文件，Utils是生成、测试赝势以及做全电子计算的脚本。在根目录下还有一个文件夹Contrib里面有一个输入参数的建议值在atom_table.txt文件里面，我们可以从建议值开始。
我们以La元素为例来演示一下如何生成赝势。

cd PS_generation
mkdir La
vi La.pbe.inp

然后从atom_table.txt里面得到La的建议值，注意这个值需要测试

   pe      Lanthanum         Guess    Ecut ~ 40Ry   l=0 as local
        tm2
 n=La c=pbr
       0.0       0.0       0.0       0.0       0.0       0.0
   11    4
    6    0     2.00      0.00
    6    1     0.00      0.00
    5    2     1.00      0.00
    4    3     0.00      0.00
   3.50     4.10     3.50     3.50

这里面最重要的是4个轨道的电子占据和截断半径，第一行前面的那个pe是用来指定这是一个什么类型的计算。pg就是生成赝势，pe就是生成赝势的时候做原子核修正（core correction），pt就是做赝势生成测试，而ae是做全电子计算。第二行是指定用TM方法生成赝势，这个相关文章里面都有，然后n=La是元素名，c=pbr是指定pbe的交换关联势，然后下面指定作为离子核的轨道的数目，相当于atomic里面的config='[Kr] 5s2 5p1 4d10'的[Kr]的轨道个数，后面一个数是需要计算的价带的数目，再下面4行就是轨道的主量子数、角量子数和占据数，最后一行是核半径（core radii $r_c$ ），单位是bohr，这个参数主要是用来做非线性交换关联修正的，详情见S. G. Louie, S. Froyen, and M. L. Cohen, Phys. Rev. B 26, 1738 (1982)..简单介绍一下子就是在传统方法中，赝核的电荷密度等于给定半径 $r_{pc}$ 的电荷密度，并且半径内有如下形式
$\rho_{pc}(r) = Ar \sin(br)$
然后程序会拟合出合适的A和b，新的用来生成GGA赝势的方法里这个函数换成了
$\rho_{pc}(r) = r^2 e^{a+br^2+cr^4}$
做核修正的时候注意第一行应该是pe而不是pg。
这个修正主要是为了加强赝势的可移植性。
设置好输入文件，然后只需要执行一下脚本就可以了

../../Utils/pg.sh La.pbe.inp

赝势就生成好了。

error in gamfind (aka v0pp) - delta not found
stop parameter =861

这个我搜到说是由于格式标准没有对齐，由于ATOM对输入格式要求十分严格，所以在写输入文件的时候最好是在给的例子上改，但是改来改去还是不行，最后还是改的半径 $r_c$ 解决了这个问题，生成赝势的时候 $r_c$ 的选取非常具有技巧性，最好是照着前面说的table.txt来选

KBgen: WARNING: Ghost states have been detected
KBgen: WARNING: Some parameter should be changed in the
KBgen: WARNING: pseudopotential generation procedure.
Stopping Program from Node:    0

这个也是由于 $r_c$ 选取的问题，改了 $r_c$ 之后，就没有这个问题了。

这个上面也有介绍ATOM的使用
Siesta 赝势atom程序学习笔记

在你要生成一个原子的赝势的时候，要首先做全电子计算，看每个轨道的本征能量。
这个时候输入文件就是这样

   ae Si Ground state all-electron
   Si   ca
       0.0
    3    2
    3    0      2.00
    3    1      2.00

12345678901234567890123456789012345678901234567890      Ruler

首先指定计算类型是ae，全电子计算。然后只需要指定交换关联势的类型，然后指定价电子轨道的主量子数和角量子数还有电子占据就可以了，最后一行是一个标尺，用来对齐上面的输入格式的，格式要求是非常严格的。
然后charge.gplot和charge.gps是gnuplot脚本，用来画电荷密度的，pseudo.gplot,pseudo.gps是画波函数，pots是画势能在实空间分布的，scpots是对比屏蔽势函数与非屏蔽势函数的对比。gplot结尾的脚本需要Xmanager，但是我装了Xmanager之后它就是一闪而过，不知道为什么。gps脚本可以通过添加gnuplot的输出格式命令让它输出为pdf或者png格式的图。

set term pngcairo lw 2 font "AR PL UKai CN,14"
set output "precipitation.png"
replot
set output

gnuplot的命令我还不太懂，用起来就感觉比较混乱。
然后做完全电子计算之后就可以生成赝势，这个时候输入文件最好换一个名字，因为它输出文件都是在一个文件夹里，这个文件夹的名字就是输入文件的名字，这里我们以Si.pbe.inp为例

#
#  Pseudopotential generation for Silicon
#  pg: simple generation
#
   pg      Silicon
        tm2      3.0             # PS flavor, logder R
 n=Si c=car                      # Symbol, XC flavor,{ |r|s}
       0.0       0.0       0.0       0.0       0.0       0.0
    3    4                       # norbs_core, norbs_valence
    3    0      2.00      0.00   # 3s2
    3    1      2.00      0.00   # 3p2
    3    2      0.00      0.00   # 3d0
    4    3      0.00      0.00   # 4f0
      1.90      1.90      1.90      1.90      0.00      0.00
#
# Last line (above): 
#    rc(s)     rc(p)     rc(d)     rc(f)   rcore_flag  rcore
#
#23456789012345678901234567890123456789012345678901234567890

这个输入文件多了很多注释，与ae的计算相比，计算的标签换成了pg,然后就是多了一个每个轨道的截止半径，半径以内，电荷密度通过一个解析函数来拟合，半径以外，与全电子计算结果相同。这个半径的选择其实并不唯一，所以需要测试，包括看输出文件的本征能量和计算电子结构以及晶格参数进行测试。通常来说半径越小，赝势越‘硬’，需要的截断能就越高，但是同时可移植性就越好，半径越大则是反过来。
做完生成赝势的计算，下面就进行赝势的测试，也就是pt的计算
，首先pt的输入格式需要ae和pt两个部分。
所以输入文件是这样的。

#
# All-electron calculations for a series of Si configurations
#
   ae Si Test -- GS 3s2 3p2
   Si   ca
       0.0
    3    2
    3    0      2.00
    3    1      2.00
   ae Si Test -- 3s2 3p1 3d1
   Si   ca
       0.0
    3    3
    3    0      2.00
    3    1      1.00
    3    2      1.00
   ae Si Test -- 3s1 3p3
   Si   ca
       0.0
    3    2
    3    0      1.00
    3    1      3.00
   ae Si Test -- 3s1 3p2 3d1
   Si   ca
       0.0
    3    3
    3    0      1.00
    3    1      2.00
    3    2      1.00
   ae Si Test -- 3s0 3p3 3d1
   Si   ca
       0.0
    3    3
    3    0      0.00
    3    1      3.00
    3    2      1.00

#
# Pseudopotential test calculations
#
   pt Si Test -- GS 3s2 3p2
   Si   ca
       0.0
    3    2
    3    0      2.00
    3    1      2.00
   pt Si Test -- 3s2 3p1 3d1
   Si   ca
       0.0
    3    3
    3    0      2.00
    3    1      1.00
    3    2      1.00
   pt Si Test -- 3s1 3p3
   Si   ca
       0.0
    3    2
    3    0      1.00
    3    1      3.00
   pt Si Test -- 3s1 3p2 3d1
   Si   ca
       0.0
    3    3
    3    0      1.00
    3    1      2.00
    3    2      1.00
   pt Si Test -- 3s0 3p3 3d1
   Si   ca
       0.0
    3    3
    3    0      0.00
    3    1      3.00
    3    2      1.00

12345678901234567890123456789012345678901234567890      Ruler

官方给的例子中是给了不同电子占据下的赝势的对比。
然后通过脚本画出对应势函数和电荷密度的对比，通过对比输出文件的本征值来看赝势的生成质量。

In-charges.png

In-pots.png
这部分主要是参考ATOM的官方手册
ATOM User Manual

siesta赝势生成

相关文章

网友评论

延伸阅读

深度阅读

栏目导航

热点阅读