首先我们要明确为什么要生成赝势。
img122.gif
首先看这个图，回忆我们计算电子波函数的时候，哈密顿量里面的势能项分为库伦势项、Hatree项和交换关联能项。如果我们直接使用库伦势的话，也就是对应图中的Z/r势，由于这个势很深，所以波函数里会包含很多内层电子的成分，这些电子速度通常很高，对应的截断能就非常高，然后为了保证互相之间的正交性，波函数会有很多震荡，与x轴的交点叫做波节。如果我们把内层电子加原子核看成一个离子，也就是在计算的时候用一个‘赝势’
生成赝势的时候首先解一个薛定谔方程
全电子波函数的角标
然后通过前面的系数来反过来拟合势函数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-charges.png
In-pots.png
这部分主要是参考ATOM的官方手册
ATOM User Manual
网友评论