24小时热门版块排行榜

返回列表

本帖产生 1 个 1ST强帖，点击这里进行查看

当前只显示满足指定条件的回帖，点击这里查看本话题的所有回帖

zhangguangping

木虫 (著名写手)

1ST强帖: 23
应助: 71 (初中生)
贵宾: 0.031
金币: 1973.1
散金: 9059
红花: 79
帖子: 2681
在线: 2128.8小时
虫号: 529624
注册: 2008-03-20
性别: GG
专业: 原子和分子物理

[交流] 关于SIESTA中如何计算功函数的一种做法探讨已有3人参与

看到小木虫上用SIESTA的虫友多次讨论怎么用SIESTA计算金属的功函数。
http://muchong.com/bbs/viewthread.php?tid=2506450

功函数对我来说是一个比较陌生的概念。虽然知道定义，但是平时的学习中和它打交道不多。但是平时老师经常问起我计算中电极的费米能级是多少，开始我看到SIESTA计算出来的金属电极的费米能级和实验测得的功函数有比较大的差别，归结为能量零点的不同导致的。也就是说在SIESTA计算中能量是有一个参考点的。参考点的能量值就是真空中的能量值。

总之，在SIESTA中求功函数，有两个量需要知道，一个是Fermi能，另一个是真空的能量。第一个很容易知道，关键是第二个怎么获得。

在上面的帖子中，老外的回信提供了一种在2.0.1版本中求功函数的方法。

这两天在看文献的时候看到文章中提到用SIESTA计算功函数的一种做法。顿时激发起我想试试的好奇心。下面把我尝试的结果分享如下。

首先：文章是Phys. Rev. B 79, 075441,已经放到附件里面了,

文章中是这样叙述的功函数的做法：
We compute the work function by plotting the electrostatic potential profile along the direction perpendicular to the surface and taking its value in the vacuum region with respect to the Fermi level.
(PS:All the surface calculations are done using a slab geometry with 18 atomic layers and a vacuum of 14 Angstrom thickness.)
文中计算的结果如下：

其中PW是利用Abinit用平面波计算的收敛结果，Exp是试验结果。其他的都是在SIESTA中利用不同大小的基组计算的。(关于基组的大小，参看文献中详细叙述）

其次：使用SIESTA-3.0-rc2版本，做尝试的时候，计算的Au，并且我选择了12layers,中间的真空间隔是14 Angstrom. 计算的时候为了减小计算量采用的是SZP,(毕竟是尝试嘛)
下面是我的输入文件：
----------------------------------------------
SystemName  bulk_Au
SystemLabel bulk_Au

NumberOfSpecies 1
NumberOfAtoms 144

%block ChemicalSpeciesLabel
1 79  Au
%endblock ChemicalSpeciesLabel

PAO.EnergyShift    100 meV

PAO.BasisType split

PAO.BasisSize SZP

LatticeConstant 1.0 Ang

%block LatticeVectors
8.65128000 0.00000000 0.00000000
0.00000000 9.98964000 0.00000000
0.00000000 0.00000000 39.9003800
%endblock LatticeVectors

#kgrid_cutoff 20. Ang

%block kgrid_Monkhorst_Pack
  6 0 0 0.5
  0 6 0 0.5
  0 0 6 0.5
%endblock kgrid_Monkhorst_Pack

# SPIN options

xc.functional       GGA          # Exchange-correlation functional
xc.authors          PBE          # Exchange-correlation version
SpinPolarized       F             # Logical parameters are: yes or no
FixSpin             F
TotalSpin          0.0
NonCollinearSpin    F                # 'T', 'F'
MeshCutoff            300. Ry             # Equivalent plane wave cutoff for the grid

# SCF options

DM.Tolerance       1.d-4       # Tolerance in maximum difference between input and output DM
MaxSCFIterations    300       # Maximum number of SCF iter
DM.UseSaveDM          T          # to use continuation files
DM.MixingWeight    0.05       # New DM amount for next SCF cycle
DM.NumberPulay       5
DM.MixSCF1          F
DM.PulayOnFile    F          # Store in memory ('F') or in files ('T')
# NeglNonOverlapInt    T          # 'F'=do not neglect
SolutionMethod       Diagon       # OrderN or Diagon or transiesta if you do a transport calculation
ElectronicTemperature  300 K       # Temp. for Fermi smearing

# MD options

MD.TypeOfRun          CG          # Type of dynamics:
MD.VariableCell       F
MD.NumCGsteps       000          # Number of CG steps for coordinate optimization
ZM.UnitsLength       Ang          #the units of length used during Z-matrix input
ZM.UnitsAngle       deg          #the units of angles used during Z-matrix input
ZM.ForceTolLength    0.02 eV/Ang #controls the convergence with respect to forces on Z-matrix lengths
ZM.ForceTolAngle    0.00356549 Ry/rad    #controls the convergence with respect to forces on Z-matrix angles
ZM.MaxDisplLength    0.1  Ang    # controls the maximum change in a Zmatrix length during an optimisation step
ZM.MaxDisplAngle    0.003 rad    # controls the maximum change in a Z-matrix angle during an optimisation step
ZM.CalcAllForces    T          # Default value

AtomicCoordinatesFormat    Ang
%block AtomicCoordinatesAndAtomicSpecies
5.767520 0.000000 0.000000    1
0.000000 0.000000 0.000000    1
2.883760 0.000000 0.000000    1
7.209400 2.497409 0.000000    1
1.441880 2.497409 0.000000    1
4.325640 2.497409 0.000000    1
5.767520 4.994818 0.000000    1
0.000000 4.994818 0.000000    1
2.883760 4.994818 0.000000    1
7.209400 7.492227 0.000000    1
1.441880 7.492227 0.000000    1
4.325640 7.492227 0.000000    1
2.883760 8.324698 2.354580    1
5.767520 8.324698 2.354580    1
0.000000 8.324698 2.354580    1
4.325640 0.832467 2.354580    1
7.209400 0.832467 2.354580    1
1.441880 0.832467 2.354580    1
2.883760 3.329876 2.354580    1
5.767520 3.329876 2.354580    1
0.000000 3.329876 2.354580    1
4.325640 5.827285 2.354580    1
7.209400 5.827285 2.354580    1
1.441880 5.827285 2.354580    1
4.325640 9.157169 4.709160    1
7.209400 9.157169 4.709160    1
1.441880 9.157169 4.709160    1
5.767520 1.664938 4.709160    1
0.000000 1.664938 4.709160    1
2.883760 1.664938 4.709160    1
4.325640 4.162347 4.709160    1
7.209400 4.162347 4.709160    1
1.441880 4.162347 4.709160    1
5.767520 6.659756 4.709160    1
0.000000 6.659756 4.709160    1
2.883760 6.659756 4.709160    1
5.767520 0.000000 7.063740    1
0.000000 0.000000 7.063740    1
2.883760 0.000000 7.063740    1
7.209400 2.497409 7.063740    1
1.441880 2.497409 7.063740    1
4.325640 2.497409 7.063740    1
5.767520 4.994818 7.063740    1
0.000000 4.994818 7.063740    1
2.883760 4.994818 7.063740    1
7.209400 7.492227 7.063740    1
1.441880 7.492227 7.063740    1
4.325640 7.492227 7.063740    1
2.883760 8.324698 9.418320    1
5.767520 8.324698 9.418320    1
0.000000 8.324698 9.418320    1
4.325640 0.832467 9.418320    1
7.209400 0.832467 9.418320    1
1.441880 0.832467 9.418320    1
2.883760 3.329876 9.418320    1
5.767520 3.329876 9.418320    1
0.000000 3.329876 9.418320    1
4.325640 5.827285 9.418320    1
7.209400 5.827285 9.418320    1
1.441880 5.827285 9.418320    1
4.325640 9.157169 11.772900    1
7.209400 9.157169 11.772900    1
1.441880 9.157169 11.772900    1
5.767520 1.664938 11.772900    1
0.000000 1.664938 11.772900    1
2.883760 1.664938 11.772900    1
4.325640 4.162347 11.772900    1
7.209400 4.162347 11.772900    1
1.441880 4.162347 11.772900    1
5.767520 6.659756 11.772900    1
0.000000 6.659756 11.772900    1
2.883760 6.659756 11.772900    1
5.767520 0.000000 25.772900    1
0.000000 0.000000 25.772900    1
2.883760 0.000000 25.772900    1
7.209400 2.497409 25.772900    1
1.441880 2.497409 25.772900    1
4.325640 2.497409 25.772900    1
5.767520 4.994818 25.772900    1
0.000000 4.994818 25.772900    1
2.883760 4.994818 25.772900    1
7.209400 7.492227 25.772900    1
1.441880 7.492227 25.772900    1
4.325640 7.492227 25.772900    1
2.883760 8.324698 28.127480    1
5.767520 8.324698 28.127480    1
0.000000 8.324698 28.127480    1
4.325640 0.832467 28.127480    1
7.209400 0.832467 28.127480    1
1.441880 0.832467 28.127480    1
2.883760 3.329876 28.127480    1
5.767520 3.329876 28.127480    1
0.000000 3.329876 28.127480    1
4.325640 5.827285 28.127480    1
7.209400 5.827285 28.127480    1
1.441880 5.827285 28.127480    1
4.325640 9.157169 30.482060    1
7.209400 9.157169 30.482060    1
1.441880 9.157169 30.482060    1
5.767520 1.664938 30.482060    1
0.000000 1.664938 30.482060    1
2.883760 1.664938 30.482060    1
4.325640 4.162347 30.482060    1
7.209400 4.162347 30.482060    1
1.441880 4.162347 30.482060    1
5.767520 6.659756 30.482060    1
0.000000 6.659756 30.482060    1
2.883760 6.659756 30.482060    1
5.767520 0.000000 32.836640    1
0.000000 0.000000 32.836640    1
2.883760 0.000000 32.836640    1
7.209400 2.497409 32.836640    1
1.441880 2.497409 32.836640    1
4.325640 2.497409 32.836640    1
5.767520 4.994818 32.836640    1
0.000000 4.994818 32.836640    1
2.883760 4.994818 32.836640    1
7.209400 7.492227 32.836640    1
1.441880 7.492227 32.836640    1
4.325640 7.492227 32.836640    1
2.883760 8.324698 35.191220    1
5.767520 8.324698 35.191220    1
0.000000 8.324698 35.191220    1
4.325640 0.832467 35.191220    1
7.209400 0.832467 35.191220    1
1.441880 0.832467 35.191220    1
2.883760 3.329876 35.191220    1
5.767520 3.329876 35.191220    1
0.000000 3.329876 35.191220    1
4.325640 5.827285 35.191220    1
7.209400 5.827285 35.191220    1
1.441880 5.827285 35.191220    1
4.325640 9.157169 37.545800    1
7.209400 9.157169 37.545800    1
1.441880 9.157169 37.545800    1
5.767520 1.664938 37.545800    1
0.000000 1.664938 37.545800    1
2.883760 1.664938 37.545800    1
4.325640 4.162347 37.545800    1
7.209400 4.162347 37.545800    1
1.441880 4.162347 37.545800    1
5.767520 6.659756 37.545800    1
0.000000 6.659756 37.545800    1
2.883760 6.659756 37.545800    1
%endblock AtomicCoordinatesAndAtomicSpecies

# Output options

WriteCoorInitial       T
WriteCoorStep          T
WriteCoorXmol          T
WriteForces             T
WriteEigenvalues       F          # If .false., it writes them in the file Systemlabel.EIG
WriteMullikenPop       1          # Write Mulliken Population Analysis
#WriteKpoints          T
#WriteKbands             T
#WriteBands             T
WriteMDCoorXmol          T
WriteMDhistory          T

# Options for saving or reading information

MD.UseSaveZM             T    # Use stored positions and velocities
MD.UseSaveCG             T    # Use stored positions and velocities
SaveRho                   T    # Write valence pseudocharge at the mesh
SaveDeltaRho             T    # Write RHOscf-RHOatm at the mesh
SaveElectrostaticPotential  T    # Write the total elect. pot. at the mesh
                              # (local pseudopotential + Hartree)
SaveTotalPotential       T    # write the valence total effective local potential
                              # (local pseudopotential + Hartree + Vxc)
WriteSiestaDim             T    # Write minimum dim to siesta.h and stop
WriteDenchar             T    # Write information for DENCHAR
----------------------------------------------------------------------------------------

最后：计算完毕，在bulk_Au.EIG文件中得到Fermi能级是：-3.1042 eV
然后我对VH,或者VT文件，做处理，找到真空的能量，数值是：1.4083 eV
最后我得到的Au的功函数是4.5125 eV.和文献中的计算结果差别不是很大。(因为我用的是SZP),如果做更精确的计算或者得到收敛值，需要增大基组。

PS:其中获得真空能量的时候利用VH和VT文件是一样的。对VH处理的结果：

(计算中我把真空放到中间了，构建了一个类似平行板电容器的结构，其实做一个slab是等价的)
其中对VH的处理有很多方法，我举一例：可以利用Util下面的Macroave。上面的图就是利用Macroave做的。至于在xy面上取不取平均都无所谓，毕竟要的数据是真空处的，这个值不依赖于取不取平均而变化,因为真空处处相等。
或者用任何能从VH或者VT得到真空中能量值的方法都行。也可以利用我在http://muchong.com/bbs/viewthread.php?tid=2386188&fpage=4分享的工具，沿着z方向的某条线提取数据，作图，得到真空的能量值。下面是我利用这种方法得到的图：

得到的真空能量为：1.4082 eV
和上面方法得到的稍微有点差别。(理论上应该没有任何区别，可能是因为编写程序的时候，用的数据精度的问题，另外VH，VT文件中的数据单位是Ry)

[ Last edited by zhangguangping on 2011-5-7 at 09:10 ]

回复此楼

» 收录本帖的淘帖专辑推荐

VASP	第一性原理相关文档	第一性原理计算经验	第一性原理和电化学
@个人收集@	计算-siesta	Siesta（Transiesta、ATK）	量子化学计算与使用经验
我的小木虫	第一性原理	学习备用	Siesta