首先我们要明确为什么要生成赝势。
首先看这个图,回忆我们计算电子波函数的时候,哈密顿量里面的势能项分为库伦势项、Hatree项和交换关联能项。如果我们直接使用库伦势的话,也就是对应图中的Z/r势,由于这个势很深,所以波函数里会包含很多内层电子的成分,这些电子速度通常很高,对应的截断能就非常高,然后为了保证互相之间的正交性,波函数会有很多震荡,与x轴的交点叫做波节。如果我们把内层电子加原子核看成一个离子,也就是在计算的时候用一个‘赝势’来代替真实的势能函数,相当于不考虑内层电子对计算结果的影响,这样波函数就只包含价电子,全电子波函数就变成没有波节的‘赝波函数’。就可以大大减小计算量。这里的分界线是自己选的一个,所以生成赝势之后需要check有没有产生非物理的态,就是所谓的Ghost State。
生成赝势的时候首先解一个薛定谔方程
全电子波函数的角标就是轨道角动量,对应的全电子波函数就是球贝塞尔函数的叠加,相当于以轨道波函数作为基矢。
然后通过前面的系数来反过来拟合势函数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 ),单位是bohr,这个参数主要是用来做非线性交换关联修正的,详情见S. G. Louie, S. Froyen, and M. L. Cohen, Phys. Rev. B 26, 1738 (1982)..简单介绍一下子就是在传统方法中,赝核的电荷密度等于给定半径的电荷密度,并且半径内有如下形式
然后程序会拟合出合适的A和b,新的用来生成GGA赝势的方法里这个函数换成了
做核修正的时候注意第一行应该是pe
而不是pg
。
这个修正主要是为了加强赝势的可移植性。
设置好输入文件,然后只需要执行一下脚本就可以了
../../Utils/pg.sh La.pbe.inp
赝势就生成好了。
error in gamfind (aka v0pp) - delta not found
stop parameter =861
这个我搜到说是由于格式标准没有对齐,由于ATOM对输入格式要求十分严格,所以在写输入文件的时候最好是在给的例子上改,但是改来改去还是不行,最后还是改的半径解决了这个问题,生成赝势的时候的选取非常具有技巧性,最好是照着前面说的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
这个也是由于选取的问题,改了之后,就没有这个问题了。
这个上面也有介绍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-pots.png
这部分主要是参考ATOM的官方手册
ATOM User Manual
网友评论