PWSCF计算电子态密度
例子，计算Cu的态密度&
一、自洽计算
calculation='scf'
restart_mode='from_scratch',
pseudo_dir = './',
outdir='./'
prefix='cu'
tstress = .true.
tprnfor = .true.
ibrav = 2, celldm(1) =6.73, nat= 1, ntyp= 1,
ecutwfc = 25.0, ecutrho = 300.0
occupations='smearing', smearing='gaussian',
degauss=0.02
&electrons diagonalization='david'&
conv_thr = 1.0e-8
mixing_beta = 0.7
ATOMIC_SPECIES
Cu 63.55 Cu.pz-d-rrkjus.UPF&
ATOMIC_POSITIONS
Cu 0.0 0.0 0.0&
K_POINTS (automatic)
8 8 8 0 0 0&
在电子自洽计算中需设置以下几个方面的参数：&
1）控制计算的部分，也就是要设置&
第一个'/'之间的关键词。
关键词calculation赋值为'scf'表示此计算是进行自洽电荷密度计算；
restart_mode表示是否是接着上一次的计算而继续的计算，赋值为'from_scratch'意味着是进行一次全新的计算开始；
pseudo_dir用来设置赝势文件所在的目录，赋值为'./'表示赝势文件放在当前计算目录；
outdir用来设置计算过程中输出文件（比如波函数、电荷密度以及势）输出到哪个目录中。赋值为'./'表示这些输出文件将放到当前计算目录中；
prefix用来定义当前计算作业的标题名，它将是一些主要输出文件的文件名。赋值为'cu'用来标记当前计算作业是对Cu进行计算；
tstress 用来设置在自洽计算过程中是否计算体系的应力，设置为 .true.表示在自洽计算过程中要计算体系的应力；
tprnfor 用来设置在自洽计算过程中是否计算体系中原子所受的力，设置为
.true.表示在自洽计算过程中要计算体系中原子所受的力;
描述所计算的体系（包括它的晶格类型、晶格常数或结构参数、原胞基矢、原胞中原子的类型数目和总的原子数目）、平面波的切断动能（也就是在展开KS轨道或
晶体波函数的平面波切断动能；另外，还包括在计算电荷密度时，展开的平面波的切断动能）、确定电子占有数的方法及相关的参数。也就是由&
..........&
第二'/'之间的关键词来设置。&
ibrav用来归属体系所属的晶格类型，赋值为2表示所计算的体系是fcc结构；
celldm(1)用来设置体系的第一个晶格常数，因为所计算的体系是fcc结构，只需设置celldm(1)，相当于指定晶格常数a的值；
nat用来指明体系的原胞中原子的总共数目，赋值为1表示所计算的原胞中只有一个原子；
ntyp用来指明体系中原子类型的数目，赋值为1表示所计算的体系只有一种类型的原子；
occupations用来设置确定电子占有数的方法，赋值为'smearing'表示采用smearing的方法来确定电子的占有数，随后须设置smearing和degauss关键词；
smearing用来指明确定电子占有数的一种具体的smearing方法，赋值为'gaussian'表示采用Gaussian函数来确定电子占有数；
degauss用来确定smearing方法中有关函数的展宽参数，赋值为0.02表示上面Gaussian函数中的展宽参数为0.02。
3）、设置电子自洽计算中本征矢量（波函数）和本征值的计算算法，自洽收敛的标准。也就是&
&electrons&
和第三个'/'之间的关键词来设置。&diagonalization用来设置在求KS方程的本征矢量和本征值时，采用具体的什么算法，赋值为'david'表示采用Davidson
iterative diagonalization with overlap
matrix方法；&
conv_thr用来设置自洽收敛标准，赋值为自洽循环过程总能的变化小于1.0e-8的化，那自洽计算就停止；
mixing_beta用来设置自洽计算过程中前后两次电荷密度混合的参数。
4）、指明体系中原子的元素名，原子量以及所采用的赝势，即ATOMIC_SPECIES
后面的设置，它们的顺序要和后面原子的坐标一一对应起来。&
Cu 63.55 Cu.pz-d-rrkjus.UPF
表示所计算的体系中原子是Cu，它的原子量为63.55，它的赝势文件为Cu.pz-d-rrkjus.UPF。
5）、给出体系原胞中原子的坐标位置，也就是ATOMIC_POSITIONS 后面的设置：
Cu 0.0 0.0 0.0
表示原胞中第一个原子是Cu，它位于原胞的原点。
6）、k点取样的设置，也就是K_POINTS 后面的设置：
K_POINTS (automatic)
表示由程序采用M-P方法自动确定k点，需给出k点取样网格的大小，以及是否在产生k点后对这些点进行平移。&
8 8 8 0 0 0
表示采用8x8x8的网格来确定k点，而且不对k点进行平移。
二、非自洽计算，增加k点，并采用四面体方法来确定电子的占有数&
calculation='nscf'
prefix='cu',
pseudo_dir
outdir='./'
celldm(1) =6.73, nat=1, ntyp=1,
= 25.0, ecutrho = 300, nbnd=8,
occupations='tetrahedra'
自洽计算中应指定nbnd，那么在自洽计算中有必要指定nbnd数量吗？如何在结果文件中查看实际使用的nband的数量。每个k点对应的能量特征值，是
不是就是band（或level, &bands consists of levels,
consists of band ；levels form
band）】&electrons
mixing_beta
ATOMIC_SPECIES
63.55 Cu.pz-d-rrkjus.UPF
ATOMIC_POSITIONS
{automatic}
三、采用dos.x计算总态密度&
outdir='./'
prefix='cu'
fildos='cu.dos',
Emin=-5.0,
Emax=25.0, DeltaE=0.1
四、采用projwfc.x来计算分波态密度 &inputpp
outdir='./'
prefix='cu'
Emin=-5.0,
Emax=25.0, DeltaE=0.1 ngauss=1,
degauss=0.02
在计算态密度的步骤就是如上面所述：a)，先进行自洽计算，保留输出的势、电荷密度和波函数；b)，然后读入上一步自洽计算得到的势或电荷密度或波函数，进行非自洽计算，其中增加k点网格，并采用四面体方法来确定电子占有数;
c)，采用dos.x计算总态密度；d)，采用projwfc.x计算分波态密度。
&【分态密度的英文原意为Atomic orbital projected density of
states，即按照不同原子往不同的轨道上投影态密度。
projwfc.x操作的英文释义：projects
wavefunctions onto orthogonalized atomic wavefunctions,
& calculates Lowdin charges, spilling parameter,
projected DOS & (separated into up and down
components for lSDA) & alternatively, computes the
local DOS(E), integrated in volumes &given in
自洽计算中occupations关键词的赋值已设置为'tetrahedra'表示采用四面体方法确定电子占有数和费米能级。另外，K_POINTS
{automatic}
下面的k点设置已增密，设置为12x12x12，为使得计算的态密度光滑，有可能需设置的更密些。
自洽计算中应指定nbnd，那么在自洽计算中有必要指定nbnd数量吗？如何在结果文件中查看实际使用的nband的数量。每个k点对应的能量特征值，是
不是就是band（或level, &bands consists of levels,
consists of band ；levels form band）】
在采用dos.x总态密度计算中，输入文件中由&inputpp
和'/'来之前的关键词来设置，它的关键词有：
outdir用来设置计算上非自洽计算输出文件的目录，设置为'./'表示是当前目录；
prefix用来标记当前所计算的体系，也确定了上一步非自洽计算输出的势或电荷密度或波函数的文件的名称，此例子中设置为'cu'，注意它们的赋值应该与上一步的非自洽计算中的一致。
fildos用来指明所计算的总态密度将写到哪个文件中，此例子中赋值为'cu.dos',表示总态密度将写到cu.dos文件；
Emin用来设置计算态密度时，能量范围的最小值，赋值为-5.0，表示将从E=-5.0
eV开始输出对应的态密度值；
Emax用来设置计算态密度时，能量范围的最大值，赋值为25.0，表示将到E=25.0eV为止输出对应的态密度；
DeltaE用来设置计算态密度时，按多大的能量间隔输出态密度，这里设为0.1eV输出态密度。
在采用projwfc.x计算态密度时，&inputpp
和'/'来之前的关键词来设置，它的关键词与dos.x的输入文件中的关键词差不多：
ngauss用来设置态密度时展宽的方法，这样是为了使得所计算的态密度看起来光滑，可以赋值：&
0，表示采用简单的高斯函数
【默认值为0】
1，表示采用一阶Methfessel-Paxton函数
-1，表示采用Marzari-Vanderbilt“冷离散“方法，其实就是一种函数形式
-99，表示采用Fermi-Dirac函数
degauss用来设置展宽函数中的展宽系数
仿照上述例子写的Sc在1atm压力下的态密度计算完整脚本为(a指压力)：
####################################################################
PW_ROOT=/home/pwscf/pwscf/espresso-4.0.5/bin/
PSEUDO_DIR=/home/workplace/pwscf-phononpy/pseudo
TMP_DIR=/home/workplace/pwscf-phononpy/tmp
#export PARA_PREFIX='mpirun -np 2'
# or & export PARA_PREFIX='mpirun' ,export
&PARA_POSTFIX= -np 3
export PW_ROOT PSEUDO_DIR TMP_DIR
export name='Sc' & &
#【定义工作环境，赝势及命令存放位置，临时文件存放位置等】
############################################################
for a in 0
cat & name.scf.in_a &&
& &calculation = 'scf'
&restart_mode='from_scratch',
& &prefix='$name',
& &pseudo_dir =
'$PSEUDO_DIR/',
& &outdir='$TMP_DIR/'
& &tstress=.t.,
& &tprnfor=.t.
&ecutwfc=90.0, ecutrho =
&occupations='smearing',
&smearing='methfessel-paxton',
&degauss=0.05
&electrons
conv_thr = 1.0d-8
mixing_beta = 0.7
#cell_dynamics= 'damp-pr' & #【自洽计算这些参数都没有必要，但为了完整性还是写出来】
# press=$a
ATOMIC_SPECIES
Sc & &44.9559 &
&Sc.pbe-nsp-van.UPF
CELL_PARAMETERS {bohr}
ATOMIC_POSITIONS {crystal}
Sc 0.4 0.0000
Sc 0.6 0.0000
K_POINTS {automatic}
8 8 4 0 0 0
PWROOT/pw.x&&/span&name.scf.in_a&name.scf.out_$a
&【#自洽计算到此结束】
##################################################################################################
cat & $name.nscf.in && EOF
& calculation='nscf',
& prefix='$name'
& pseudo_dir = '$PSEUDO_DIR/',
& outdir='$TMP_DIR/'
&system & &
& &ibrav = 0,
& nat= &2,
& ntyp= 1,
& ecutwfc = 90,
& ecutrho = 400.00,
& occupations='smearing',
& smearing='methfessel-paxton',
& degauss=0.05
&electrons
& mixing_beta = 0.7
& conv_thr = &1.0d-8
#cell_dynamics= 'damp-pr' & #【自洽计算这些参数都没有必要，但为了完整性还是写出来】
# press=$a
ATOMIC_SPECIES
Sc & &44.9559 &
&Sc.pbe-nsp-van.UPF
CELL_PARAMETERS {bohr}
ATOMIC_POSITIONS {crystal}
Sc 0.4 0.0000
Sc 0.6 0.0000
K_POINTS { automatic }
12 12 6 &0 0 0
PWROOT/pw.x&&/span&name.nscf.in
&$name.nscf.out & 【#至此非自洽部分结束。这里只是改了一个参数就是calculation='nscf',注意这里的K点可以选大点因为我们需要计算的电荷密度应该比较准。】
######################################################################################################
cat & $name.Dos.in && EOF
& prefix &= '$name'
& outdir = '$TMP_DIR/'
& fildos = '$name_dos.dat'
& DeltaE=0.01
PWROOT/dos.x&&/span&name.Dos.in
& $name.Dos.out &【#计算总的态密度】
####################################################
cat &$name.pdos.in && EOF
& outdir='$TMP_DIR/'
& prefix='$name'
& DeltaE=0.01
PWROOT/projwfc.x&&/span&name.pdos.in
&$name.pdos.out & 【#计算投影态密度】
###################################################################
已投稿到：
以上网友发言只代表其个人观点，不代表新浪网的观点或立场。小木虫 --- 500万硕博科研人员喜爱的学术科研平台
&&查看话题
【求助】quantum-espresso pwscf yambo
几者之间的关系？
最近学习quantum－espresso 和yambo 想作gw计算，看见yambo里面说用yambo作gw计算需要pwscf和abinit，请问，在yambo里面作gw原算时，需要pwscf和abinit提供什么数据？pwscf在不需要yambo的情况下能不能进行gw计算？quantum-espresso pwscf yambo 几者之间的关系？
谢谢提供信息。会不会在往后的一两个升级版本中增加计算多点GW的功能？
这是肯定的。据说他们已经搞定了。就是不知道会不会近期发布
研究生必备与500万研究生在线互动！
扫描下载送金币
浏览器进程
打开微信扫一扫
随时随地聊科研[转载]常用PWscf帮助文档及网址
已有 4393 次阅读
|个人分类:|系统分类:|关键词:网址|文章来源:转载
官方用户指南，PWscf的编译、安装、使用、常见问题等信息非常全面，通俗易懂。初学者必读。结合examples使用
PWscf安装目录下的Examples包含30多个常见实例，直接将run_example的shell脚本拖入终端即可运行。运行后在results文件夹下可以看对应输入文件和输出结果。可以参照示例写自己的输入文件。
PWscf安装目录下的Doc文件夹内有许多html文件，内容包括各种输入文件中所有变量和标签的详细定义和规范，在写输入文件的时候一定要仔细阅读。
http://www.democritos.it/mailman/listinfo/pw_users
在Google里搜索PWscf的问题排在第一的往往就是PWscf官方的邮件列表。这里汇集了全球PWscf用户的常见问题及解答。你的问题往往就在里面。对于常见的小问题这里的答案往往简单明了，事半功倍。有需要的同学可以订阅。
如何安装编译PWscf？
单机安装简易版（适用Ubuntu32位操作系统）from&
首先安装编译环境，确认联网并，在终端下输入：
sudo apt-get install gfortran
sudo apt-get install mpich-bin
然后去PWscf网站下载免费的源码包&
解压到任意目录
在终端窗口进入解压目录，方法如下：
输入：&cd&路径&（可直接把文件夹拖拽进终端自动生成路径）
./configure
等屏幕上一堆字刷完再输入
继续等刷屏 大概十分钟以后应该就好了
测试软件 进入解压目录下的examples文件夹
看到有很多example01 example02等目录
随便进一个 把run_example文件拖进终端窗口
看运行是否正常
也可以运行examples文件夹下的run_all_examples文件
运行所有的示例
计算结果在各自文件夹下的results目录里
要清除计算结果 运行examples下的make clean文件
本教程应该适用于32位Ubuntu为基础的所有衍生版
64位Ubuntu源里好像没有mpich，需要自己
一个台湾人写的安装教程（适用于Ubuntu8.10）
前幾天應數所的學弟玩 Linux 碰到問題，原來是不會裝這個...PWscf 是用來做工程及科學計算的，PWscf 的介紹如下，總之先幫他裝起來吧！上的介紹：
PWscf (Plane-Wave Self-Consistent Field) is a computer code for electronic-structure calculations within Density-Functional Theory and Density-Functional Perturbation Theory, using pseudopotentials and a plane-wave basis set. PWscf is part of the&the atomic scale. PWscf is released under the&.&distribution of codes for the quantum simulation of matter at因為學弟的電腦是 Ubuntu 8.10 ，我也不知道在其他 distribution 可不可以依照相同的步驟安裝；這裡列出一些套件，找相對應的裝應該就差不多了吧。下載PWscf
GUI 介面：
安裝 pwscf1. 預先安裝套件，請安裝下列套件，或&sudo apt-get install 套件名稱&即可
2.&，我們需要的是 ssh 自動登入，把其中設定 ssh 自動登入的步驟列舉如下
基本步驟摘要 (只用於一台）：
設置多台 ssh 可進入之免密碼機器組
3. 打開終端機，cd espresso-4.0.44.&./configure5.&make all安裝 pwgui6.&tar zxvf PWgui-4.1CVS.tgz7.&cd PWgui-4.1CVS8. 執行PWGui:&./pwgui註：後來發現 pwscf 的目錄裡本來就帶有 pwgui 了，或許可以直接執行 ./pwgui 測試看看可不可以使用。設定 pwgui設定 pwgui 可執行檔的路徑如下圖，此例中 pwscf 被解壓縮在這個目錄下：
測試 pwscfpwscf 已經提供非常多 example 在原始碼裡面，在我們的路徑下為：
以 example01 為例：
計算的結果會存在 results/ 目錄下面
by Sunforever
首先，在MS中完成建模
然后用castep进行运算，参数随便，只要run一下就可以了，开始计算后马上把任务停掉
去硬盘上你的计算任务路径下的/Documents目录，在文件夹选项里设置“显示隐藏文件”,找到.cell文件，用记事本打开
打开后就可以看到原子坐标了。这里的坐标单位是晶胞长，坐标值是晶胞尺寸的分数。在PWscf中使用时要把单位改成{crystal}
&&C & 0.9025 & 0.0661 & 0.6403
&&C & 0.9025 & 0.0661 & 0.6402
&&O & 0.6667 & 0.3333 & 0.9997
&&O & 0.3333 & 0.6667 & 0.9997
&&O & 0.6666 & 0.3333 & 0.9997
&&O & 0.3334 & 0.6669 & 0.9997
&&O & 0.6667 & 0.3333 & 0.9997
&&O & 0.3333 & 0.6667 & 0.9997
&&O & 0.6666 & 0.3334 & 0.9997
&&O & 0.3334 & 0.6667 & 0.9997
&&O & 0.6667 & 0.3333 & 0.9997
&&O & 0.3333 & 0.6667 & 0.9998
&&O & 0.6666 & 0.3333 & 0.9997
&&O & 0.3334 & 0.6669 & 0.9998
&&O & 0.6667 & 0.3333 & 0.9997
&&O & 0.3333 & 0.6667 & 0.9998
&&O & 0.6666 & 0.3334 & 0.9997
&&O & 0.3334 & 0.6667 & 0.9998
&Zn & 0.6667 & 0.3333 & 0.0001
&Zn & 0.3333 & 0.6667 & 0.0000
&Zn & 0.6666 & 0.3333 & 0.0001
&Zn & 0.3334 & 0.6669 & 0.0000
&Zn & 0.6667 & 0.3333 & 0.0001
&Zn & 0.3333 & 0.6667 & 0.0000
&Zn & 0.6666 & 0.3334 & 0.0001
&Zn & 0.3334 & 0.6667 & 0.0000
&Zn & 0.6667 & 0.3333 & 0.0001
&Zn & 0.3333 & 0.6667 & 0.0000
&Zn & 0.6666 & 0.3333 & 0.0001
&Zn & 0.3334 & 0.6669 & 0.0000
&Zn & 0.6667 & 0.3333 & 0.0001
&Zn & 0.3333 & 0.6667 & 0.0000
&Zn & 0.6666 & 0.3334 & 0.0001
&Zn & 0.3334 & 0.6667 & 0.0000
PWscf输入文件原子坐标部分示例如下：
ATOMIC_POSITIONS&{crystal}
&&C & 0.9025 & 0.0661 & 0.6403
&&C & 0.9025 & 0.0661 & 0.6402
&&O & 0.6667 & 0.3333 & 0.9997
&&O & 0.3333 & 0.6667 & 0.9997
&&O & 0.6666 & 0.3333 & 0.9997
&&O & 0.3334 & 0.6669 & 0.9997
&&O & 0.6667 & 0.3333 & 0.9997
&&O & 0.3333 & 0.6667 & 0.9997
&&O & 0.6666 & 0.3334 & 0.9997
&&O & 0.3334 & 0.6667 & 0.9997
&&O & 0.6667 & 0.3333 & 0.9997
&&O & 0.3333 & 0.6667 & 0.9998
&&O & 0.6666 & 0.3333 & 0.9997
&&O & 0.3334 & 0.6669 & 0.9998
&&O & 0.6667 & 0.3333 & 0.9997
&&O & 0.3333 & 0.6667 & 0.9998
&&O & 0.6666 & 0.3334 & 0.9997
&&O & 0.3334 & 0.6667 & 0.9998
&Zn & 0.6667 & 0.3333 & 0.0001
&Zn & 0.3333 & 0.6667 & 0.0000
&Zn & 0.6666 & 0.3333 & 0.0001
&Zn & 0.3334 & 0.6669 & 0.0000
&Zn & 0.6667 & 0.3333 & 0.0001
&Zn & 0.3333 & 0.6667 & 0.0000
&Zn & 0.6666 & 0.3334 & 0.0001
&Zn & 0.3334 & 0.6667 & 0.0000
&Zn & 0.6667 & 0.3333 & 0.0001
&Zn & 0.3333 & 0.6667 & 0.0000
&Zn & 0.6666 & 0.3333 & 0.0001
&Zn & 0.3334 & 0.6669 & 0.0000
&Zn & 0.6667 & 0.3333 & 0.0001
&Zn & 0.3333 & 0.6667 & 0.0000
&Zn & 0.6666 & 0.3334 & 0.0001
&Zn & 0.3334 & 0.6667 & 0.0000
就这么简单
运行castep那一步是为了生成cell文件，所以只要运行一下下就好
感谢Cutie同学的指导
如何进行结构优化？
在pwscf中提供了两种优化方法对原子位置进行驰豫：a) BFGS quasi-newton algorithm, b)最速下降法（quick-min Verlet）。除非初始位置很接近平衡位置，一般采用BFGS准牛顿算法比较快。在结构优化时，calculation需设置为"relax"，下面相关的参数有时也需要设置：
一、在&control .../中设置优化的收敛标准参数、步数等
nstep&：优化的最大步数；
tstress&：设置True，表示要计算体系的应力；
tprnfor：设置为True，表示要计算原子所受的力，在calculation='relax'时，默认为.True.；
etot_conv_thr：用来控制原子位置优化时，总能变化收敛的标准；默认值为1.0D-4 Ha；
forc_conv_thr：用来控制原子位置优化时，原子所受力的收敛标准，默认值为1.0D-3 Ha/a.u，只有当etot_conv_thr和forc_conv_thr的标准都满足时，优化才停止；
二、在&IONS..../中设置优化方法中的相关参数
在calculation='relax'时，ion_dynamics设置为'bfgs'表示用BFGS准牛顿算法来进行优化，设置为'damp'表示用最速下降法来进行优化。
pot_extrapolation：用来控制优化或电子迭代过程中势的混合方式，在原子位置优化时，最好设置为'second_order'，表示采用二阶方式来进行混合；
wfc_extrapolation:用来控制优化或电子迭代过程中波函数的混合方式，在原子位置优化时，最好设置为'second_order'，表示采用二阶方式来进行混合；
当设置了采用BFGS准牛顿算法来优化后，下面的参数需设置：
upscale：用来控制conv_thr在结构优化过程中最可能的缩小值，在优化快接近收敛时，conv_thr会自动减小以保证所计算的力仍然很精确，但是conv_thr并不会减小到conv_thr/upscale；
bfgs_ndim:用来控制有多少个旧的力和位移矢量会用在残差矢量的Pulay混合中，其中残差矢量是基于BFGS算法中的Hessian矩阵的逆来得到的。设置为1，就是标准的BFGS准牛顿算法；
trust_radius_max：离子优化过程中，离子每一步最大移动量； 默认值为0.8D0；
trust_radius_min：离子优化过程中，当trust_radius小于trust_radius_min&时，离子每一步最小移动量； 默认值为&1.D-3;
trust_radius_ini:默认值为&0.5D0，在原子位置优化计算中初始的离子位移量；
w_1, w_2&：用在基于Wolfe条件的线性搜索方法中的参数。
例子1：优化CO分子
calculation = "relax",
prefix = "CO",
pseudo_dir = "./",
outdir = "./",
tprnfor=.true.
forc_conv_thr=1.0D-4
ibrav = 1, celldm(1) = 12.D0, nat = 2, ntyp = 2,
ecutwfc = 24.D0, ecutrho = 144.D0,
&ELECTRONS
conv_thr = 1.D-7,
mixing_beta = 0.7D0,
ion_dynamics='bfgs',
pot_extrapolation = "second_order",
wfc_extrapolation = "second_order",
ATOMIC_SPECIES
O 1.00 O.LDA.US.RRKJ3.UPF
C 1.00 C.pz-rrkjus.UPF
ATOMIC_POSITIONS {bohr}
C 2.256 0.0 0.0
O 0.000 0.0 0.0 0 0 0
K_POINTS {Gamma}
如何进行自洽计算?
PWSCF程序包（早期的叫法），或称为ESPRESSO程序（改名后的叫法)，它包括了多几个计算模块，主要的是电子自洽计算模块pw.x，晶格动力学计算模块(ph.x, phcg.x, dynmat.x,d3.x等)，后续数据处理模块pp.x，电子输运性质计算模块pwcond.x，分子动力学模块cp.x等
一、自洽计算
例子：fcc Cu的自洽计算
calculation='scf'
restart_mode='from_scratch',
pseudo_dir = './',
outdir='./'
prefix='cu'
tstress = .true.
tprnfor = .true.
ibrav = 2, celldm(1) =6.73, nat= 1, ntyp= 1,
ecutwfc = 25.0, ecutrho = 300.0
occupations='smearing', smearing='gaussian', degauss=0.02
&electrons
diagonalization='david'
conv_thr = 1.0e-8
mixing_beta = 0.7
ATOMIC_SPECIES
Cu 63.55 Cu.pz-d-rrkjus.UPF
ATOMIC_POSITIONS
Cu 0.0 0.0 0.0
K_POINTS (automatic)
8 8 8 0 0 0
在电子自洽计算中需设置以下几个方面的参数：
1）控制计算的部分，也就是要设置
第一个'/'之间的关键词。
关键词calculation赋值为'scf'表示此计算是进行自洽电荷密度计算；
restart_mode表示是否是接着上一次的计算而继续的计算，赋值为'from_scratch'意味着是进行一次全新的计算开始；
pseudo_dir用来设置赝势文件所在的目录，赋值为'./'表示赝势文件放在当前计算目录；
outdir用来设置计算过程中输出文件（比如波函数、电荷密度以及势）输出到哪个目录中。赋值为'./'表示这些输出文件将放到当前计算目录中；
prefix用来定义当前计算作业的标题名，它将是一些主要输出文件的文件名。赋值为'cu'用来标记当前计算作业是对Cu进行计算；
tstress&用来设置在自洽计算过程中是否计算体系的应力，设置为&.true.表示在自洽计算过程中要计算体系的应力；
tprnfor&用来设置在自洽计算过程中是否计算体系中原子所受的力，设置为&.true.表示在自洽计算过程中要计算体系中原子所受的力;
2)、描述所计算的体系（包括它的晶格类型、晶格常数或结构参数、原胞基矢、原胞中原子的类型数目和总的原子数目）、平面波的切断动能（也就是在展开KS轨道或晶体波函数的平面波切断动能；另外，还包括在计算电荷密度时，展开的平面波的切断动能）、确定电子占有数的方法及相关的参数。也就是由
..........
第二'/'之间的关键词来设置。
ibrav用来归属体系所属的晶格类型，赋值为2表示所计算的体系是fcc结构；
celldm(1)用来设置体系的第一个晶格常数，因为所计算的体系是fcc结构，只需设置celldm(1)，相当于指定晶格常数a的值；
nat用来指明体系的原胞中原子的总共数目，赋值为1表示所计算的原胞中只有一个原子；
ntyp用来指明体系中原子类型的数目，赋值为1表示所计算的体系只有一种类型的原子；
occupations用来设置确定电子占有数的方法，赋值为'smearing'表示采用smearing的方法来确定电子的占有数，随后须设置smearing和degauss关键词；
smearing用来指明确定电子占有数的一种具体的smearing方法，赋值为'gaussian'表示采用Gaussian函数来确定电子占有数；
degauss用来确定smearing方法中有关函数的展宽参数，赋值为0.02表示上面Gaussian函数中的展宽参数为0.02。
3）、设置电子自洽计算中本征矢量（波函数）和本征值的计算算法，自洽收敛的标准。也就是
&electrons
和第三个'/'之间的关键词来设置。
diagonalization用来设置在求KS方程的本征矢量和本征值时，采用具体的什么算法，赋值为'david'表示采用Davidson iterative diagonalization with overlap matrix方法；
conv_thr用来设置自洽收敛标准，赋值为自洽循环过程总能的变化小于1.0e-8的化，那自洽计算就停止；
mixing_beta用来设置自洽计算过程中前后两次电荷密度混合的参数。
4）、指明体系中原子的元素名，原子量以及所采用的赝势，即ATOMIC_SPECIES&后面的设置，它们的顺序要和后面原子的坐标一一对应起来。
Cu 63.55 Cu.pz-d-rrkjus.UPF
表示所计算的体系中原子是Cu，它的原子量为63.55，它的赝势文件为Cu.pz-d-rrkjus.UPF。
5）、给出体系原胞中原子的坐标位置，也就是ATOMIC_POSITIONS&后面的设置：
Cu 0.0 0.0 0.0
表示原胞中第一个原子是Cu，它位于原胞的原点。
6）、k点取样的设置，也就是K_POINTS&后面的设置：
K_POINTS (automatic)&表示由程序采用M-P方法自动确定k点，需给出k点取样网格的大小，以及是否在产生k点后对这些点进行平移。
8 8 8 0 0 0
表示采用8x8x8的网格来确定k点，而且不对k点进行平移。&
如何计算能带结构？
计算fcc Cu的能带结构
calculation='bands'
pseudo_dir = './',
outdir='./',
prefix='cu'
ibrav = 2, celldm(1) =6.73, nat= 1, ntyp= 1,
ecutwfc = 25.0, ecutrho = 300.0, nbnd = 8
&electrons
diagonalization='david'
ATOMIC_SPECIES
Cu 63.55 Cu.pz-d-rrkjus.UPF
ATOMIC_POSITIONS
Cu 0.0 0.0 0.0
0.0 0.0 0.0 1.0
0.0 0.0 0.1 1.0
0.0 0.0 0.2 1.0
0.0 0.0 0.3 1.0
0.0 0.0 0.4 1.0
0.0 0.0 0.5 1.0
0.0 0.0 0.6 1.0
0.0 0.0 0.7 1.0
0.0 0.0 0.8 1.0
0.0 0.0 0.9 1.0
0.0 0.0 1.0 1.0
0.0 0.0 0.0 1.0
0.0 0.1 0.1 1.0
0.0 0.2 0.2 1.0
0.0 0.3 0.3 1.0
0.0 0.4 0.4 1.0
0.0 0.5 0.5 1.0
0.0 0.6 0.6 1.0
0.0 0.7 0.7 1.0
0.0 0.8 0.8 1.0
0.0 0.9 0.9 1.0
0.0 1.0 1.0 1.0
0.0 0.0 0.0 1.0
0.1 0.1 0.1 1.0
0.2 0.2 0.2 1.0
0.3 0.3 0.3 1.0
0.4 0.4 0.4 1.0
0.5 0.5 0.5 1.0
在进行能带计算时，calculation须设置为'bands'，而且在此之前须进行一次相应的自洽计算，而且要有上一步计算得到输出文件供能带计算时读入。另外最好在&system中设置nbnd，以指定计算多少条能带。在计算能带时要自己先选定一些高对称点，并产生这些高对称点之间其他点。在这个例子中，计算沿G-X-L点之间的高对称线上的能带。
在产生所要计算的特殊k点时，可以采用下面简单的f77程序来实现：
c +---------------------------------------------------------
c For generating k-points along the high-symmetry lines in
c Brillouin zone and for calculate band-structures !
c +----------------------------------------------------------
C ---------'syml'---------
c 6 : nhighk
c 20 20 20 10 20 : ndiv(i)
c X 0.5 0.0 0.5 : labhk(1),phighk(1,1),........
c G 0.0 0.0 0.0
c L 0.5 0.5 0.5
c W 0.5 0.25 0.75
c K 0.375 0.375 0.75
c G 0.0 0.0 0.0
c direct & reciprocal lattice vectors over 'emin, emax' line
C -----------------------
c max k-points = 200
program gk
implicit real*8 (a-h,o-z)
character*2 labhk
dimension tkpt(200,3),pk(200,3),phighk(10,3)
dimension disk(200),dish(10),labhk(10)
dimension ndiv(10)
open(5,file='syml',status='old')
open(7,file='inp.kpt')
read(5,*) nhighk
read(5,*) (ndiv(i),i=1,nhighk-1)
do i=1,nhighk-1
ntkp=ntkp+ndiv(i)
ntotkpt=ntkp+1
if(nhighk&10)then
write(*,*)'Number of high-symmetry k points must & 10!'
if(ntotkpt&200)then
write(*,*)'Total number of k points must &= 200!'
do i=1, nhighk
read(5,*) labhk(i),(phighk(i,j),j=1,3)
write(*,*) (labhk(i),i=1,nhighk)
c----- generating k-points along high symmetric lines --------
pk(1,1)=phighk(1,1)
pk(1,2)=phighk(1,2)
pk(1,3)=phighk(1,3)
do i = 2, nhighk
delx = (phighk(i,1) - phighk(i-1,1))/float(ndiv(i-1))
dely = (phighk(i,2) - phighk(i-1,2))/float(ndiv(i-1))
delz = (phighk(i,3) - phighk(i-1,3))/float(ndiv(i-1))
do j=1, ndiv(i-1)
ii = ii + 1
pk(ii,1) = pk(ii-1,1) + delx
pk(ii,2) = pk(ii-1,2) + dely
pk(ii,3) = pk(ii-1,3) + delz
10 format(A34)
weight=1.d0
do i=1,ntotkpt
write(7,200) pk(i,1),pk(i,2),pk(i,3),weight
200 format(3F10.6,F6.2)
c----------------------- end ---------------------------
它的输入文件为syml，输出文件为inp.kpt。其中syml输入文件的格式如下：
15 15 15 15 15 15 15
G 0.0 0.0 0.0
K -0. 0.7 0.
H -0. 0.7 0.
A 0.0 0.0 0.5
G 0.0 0.0 0.0
M 0.0 0.5 0.0
L 0.0 0.5 0.5
A 0.0 0.0 0.5
第一行用来标记有多少个特殊k点，下面是这些特殊k点之间每个要分多少个k点，接着就是这些特殊k点的坐标。
产生的inp.kpt可以之间拷贝到pw.x在计算能带时的输入文件中。
如何画能带图？
pwscf&附带了band.x和plotband.x的工具，前者是将计算出来k点坐标以及相应的本征值从out文件中收集起来（或取出）专门存储到一个文件中，以便后一个工具plotband.x进行处理。在计算能带时，先设置k点网格进行一次自洽计算，然后自己输入要计算的特殊k点并进行一次非自洽计算，得到这些特殊k点的本征值。再用band.x和ploband.x进行处理。
1). band.x的输入文件格式
prefix = 'si'
outdir = './tmp'
filband = 'sibands.dat'
spin_component=1
其中&prefix设置所计算体系的标题，以及输入文件的文件名（不包括扩展名，也就是. 'dot'后面的));
outdir用来设置上一步非自洽计算中的输出文件的目录;
filband用来设置这一步band.x处理出来的k点和本征值的输出文件(也就是将这些k点-本征值放到哪个文件中);
spin_component为1表示处理的是非自旋极化计算的本征值，如果是2表示处理的是自旋极化计算的本征值。
band.x的输出文件(fiband所设置的)的格式为：
&plot nbnd= 8, nks= 36 /
0...500000
-3.418 -0.822 5.029 5.029 7.814 9.597 9.597 13.838
0...400000
-3.891 -0.102 5.102 5.102 7.900 9.679 9.679 13.959
0...300000
-4.659 1.404 5.319 5.319 8.138 9.803 9.803 13.845
。。。。。。。。。
与在声子计算中matdyn.x得到的q点与本征值的文件的格式是一样的。
第一行中nbnd告诉了每个有多少个本征值，本征值单位是eV。nks告诉了总共有多少k点。第2行是k点的坐标，第3行是该k点对应的本征值。下面的与2、3两行类似。
2)、plotband.x的输入文件格式
sibands.dat
sibands.xmgr
sibands.ps
第一行是band.x处理得到的k点本征值文件 此例子中是sibands.
第二行是在这一步中输出的ps文件中纵坐标（本征值能量）刻度的最小值与最大值;
第三行是用来设置所出输出xmgr格式的文件的文件名;
第四行是用来设置所输出文件ps格式的文件名;
第五行是费米能级的值，这个在自洽计算的out文件中可以找到;
第六行中第一个数是用来设置输出&ps格式文件纵坐标（能量）刻度的大小，第二个数是用来设置标出费米能级的位置，它与第五行中的数相同。
同样ploband.x也可以用来处理声子本征值文件。
如何计算态密度？
例子，计算Cu的态密度
一、自洽计算（见前一个例子）
二、非自洽计算，增加k点，并采用四面体方法来确定电子的占有数
calculation='nscf'
prefix='cu',
pseudo_dir = './',
outdir='./'
ibrav=2, celldm(1) =6.73, nat=1, ntyp=1,
ecutwfc = 25.0, ecutrho = 300, nbnd=8, occupations='tetrahedra'
&electrons
conv_thr = 1.0e-8
mixing_beta = 0.7
ATOMIC_SPECIES
Cu 63.55 Cu.pz-d-rrkjus.UPF
ATOMIC_POSITIONS
Cu 0.0 0.0 0.0
K_POINTS {automatic}
12 12 12 0 0 0
三、采用dos.x计算总态密度
outdir='./'
prefix='cu'
fildos='cu.dos',
Emin=-5.0, Emax=25.0, DeltaE=0.1
四、采用projwfc.x来计算分波态密度
outdir='./'
prefix='cu'
Emin=-5.0, Emax=25.0, DeltaE=0.1 ngauss=1, degauss=0.02
在计算态密度的步骤就是如上面所述：a)，先进行自洽计算，保留输出的势、电荷密度和波函数；b)，然后读入上一步自洽计算得到的势或电荷密度或波函数，进行非自洽计算，其中增加k点网格，并采用四面体方法来确定电子占有数; c)，采用dos.x计算总态密度；d)，采用projwfc.x计算分波态密度。
在自洽计算中occupations关键词的赋值已设置为'tetrahedra'表示采用四面体方法确定电子占有数和费米能级。另外，K_POINTS {automatic}&下面的k点设置已增密，设置为12x12x12，为使得计算的态密度光滑，有可能需设置的更密些。
在采用dos.x总态密度计算中，输入文件中由&inputpp&和'/'来之前的关键词来设置，它的关键词有：
outdir用来设置计算上非自洽计算输出文件的目录，设置为'./'表示是当前目录；
prefix用来标记当前所计算的体系，也确定了上一步非自洽计算输出的势或电荷密度或波函数的文件的名称，此例子中设置为'cu'，注意它们的赋值应该与上一步的非自洽计算中的一致。
fildos用来指明所计算的总态密度将写到哪个文件中，此例子中赋值为'cu.dos',表示总态密度将写到cu.dos文件；
Emin用来设置计算态密度时，能量范围的最小值，赋值为-5.0，表示将从E=-5.0 eV开始输出对应的态密度值；
Emax用来设置计算态密度时，能量范围的最大值，赋值为25.0，表示将到E=25.0eV为止输出对应的态密度；
DeltaE用来设置计算态密度时，按多大的能量间隔输出态密度，这里设为0.1eV输出态密度。
在采用projwfc.x计算态密度时，&inputpp&和'/'来之前的关键词来设置，它的关键词与dos.x的输入文件中的关键词差不多：
ngauss用来设置态密度时展宽的方法，这样是为了使得所计算的态密度看起来光滑，可以赋值：
0，表示采用简单的高斯函数
1，表示采用一阶Methfessel-Paxton函数
-1，表示采用Marzari-Vanderbilt“冷离散“方法，其实就是一种函数形式
-99，表示采用Fermi-Dirac函数
degauss用来设置展宽函数中的展宽系数。
如何计算晶体的声子色散曲线和态密度？
pwscf是采用线性响应的方法来进行晶格动力学性质的计算。在计算晶体的声子色散和态密度时的步骤：i)用pw.x进行自洽计算; ii)用ph.x对小的q点网格进行计算，得到这些q点的动力学矩阵元;iii)用q2r.x计算出实空间中的力常数矩阵; iv)用matdyn.x计算声子色散曲线;v)用matdyn.x计算声子态密度。
下面以Sc为例子并针对pwscf的最新版本3.2.1来说明（早期版本在计算声子色散曲线较麻烦，因为它不能自动处理q点网格，然后对每个q点一次性计算，而是需要手动产生这些点，一个个计算)。
1)&用pw.x进行电子密度的自洽计算
title='Sc, hexagonal cell'
calculation = 'scf'
restart_mode='from_scratch',
prefix='sc',
pseudo_dir = './',
outdir='./tmp'
tprnfor=.true.
celldm(1)=6.05606,
celldm(3)=1.71298,
ntyp=1, nbnd= 30,
ecutwfc=30.0,
occupations ='smearing', degauss =0.01
smearing ='mp'
&electrons
diagonalization='cg'
diago_cg_maxiter= 60
mixing_mode = 'plain'
mixing_beta = 0.5
conv_thr = 1.0d-6
ATOMIC_SPECIES
Sc 44.955910 Sc.pw91-nsp-van.UPF
ATOMIC_POSITIONS (crystal)
Sc 0.4 0.0000
Sc 0.6 0.0000
K_POINTS (automatic)
8 8 6 0 0 0
注意Sc是金属，在此例子中，我们选用MP方法来确定电子的占有数(见occupations ='smearing', smearing ='mp')，这里未经测试而选用了展宽系数为0.01 Ry (见degauss=0.01)。在进行声子色散曲线的计算时，不必需对&calculation设置为'phonon'，在新版本中，直接设置为scf。这一步计算产生势以及电荷密度供一下的计算中利用到。
2)用ph.x对小的q网格点进行动力学矩阵元的计算
phonon for Sc
tr2_ph=1.0d-10,
prefix='sc',
fildvscf='scdv',
amass(1)=44.955910,
outdir='./tmp',
fildyn='sc.dyn',
elph=.false.,
trans=.true.,
ldisp=.true.
nq1=4, nq2=4, nq3=2
注意这里trans和ldisp必须设置为.true.。其中trans为.true.表示要计算声子相关的性质，ldisp设置为.true.表示要计算声子色散曲线。另外&prefix和outdir的设置尽量与上一步自洽计算中的设置一致，以能读入上一步计算得到的数据。另外nq1,nq2和nq3是用来设置q网格点的。为了得到实空间的力常数矩阵，这里采用的是先计算出q空间中小的q网格点的动力矩阵元，然后采用fft变换得到实空间的力常数矩阵。因此在这一步计算中需设置小的q网格点的网格大小。
3)用q2r.x计算实空间力常数矩阵
zasr='simple', fildyn='sc.dyn', flfrc='Sc444.fc', la2F=.false.
在q2r.x的输入文件中需指定
fildyn： 用来设置包含了q网格点的动力学矩阵元的文件，与上一步的&fildyn设置一致;
flfrc：用来设置输出力常数矩阵的文件;
la2F：用来设置是否计算出实空间中电-声耦合系数;针对计算材料的超导性质;
zasr:如何处理‘声学支求和规则“，该规则是用在处理Born有效电荷的，要求Born有效电荷的总和是零。可赋的值有:
no,表示不处理声学支求和问题
simple,&表示通过对力常数矩阵的对角元素进行修正来考虑3支声学横模的求和处理;
这里我们设置为simple。
4)用matdyn.x计算出声子色散曲线
asr='simple', amass(1)=44.955910,
flfrc='Sc444.fc', flfrq='Sc444.freq', la2F=.false., dos=.false.
这里要输入131个特殊q点的坐标。与计算能带结构时一样，需先选出要计算的高对称q点的走向以及高对称点的坐标，然后产生这些线上的q点的坐标。
计算出来的每一个q点的本征频率可按上一个blog中提到的方法处理一下后画图。
5),用matdyn.x计算声子态密度
asr='simple', amass(1)=44.955910,
flfrc='Sc444.fc', flfrq='Sc444.freq', la2F=.false., dos=.true.
fldos='phonon.dos', nk1=10, nk2=10, nk3=10, ndos=50
要计算声子态密度,dos必须设置为.true.另外fldos用来设置输出的态密度值，计算态密度时要用更密的q点网格，这需设置nk1, nk2, nk3。另外还有态密度的能量刻度上的点的数目，由ndos来设置。
如何画出声子色散曲线？
在采用pwscf中的ph.x以及matdyn.x计算得到了沿高对称q点的本征频率后，需将数据转换一下，以画出声子色散曲线。下面是计算出来的本征频率文件中的格式：
&plot nbnd= 6, nks= 91 /
0...500000
139.1 277.5 460.6
0...475000
136.6 269.2 463.2
0...450000
132.8 260.5 466.0
0...425000
128.5 252.6 468.3
0...400000
124.2 243.2 471.1
0...375000
119.9 233.3 474.8
0...350000
114.2 223.6 477.5
0...325000
.......................
第一行中的nbnds告诉了有每个q点有多少个本征振动频率，nks告诉了计算了多少个q点。下面的接着的每两行是：
每个q点对应的本征振动频率
要将q点的坐标转换为2D图中的横坐标上的长度值，下面的fortran90代码就是将这样的计算数据转换为q的长度值以及对应的本征频率值，以便画图。
! for PWSCF to plot phonon dispersion
program prog
real, allocatable :: e(:,:)
real, allocatable :: k(:,:)
real, dimension(3) ::k0,a
character(len=32):: input, output
character(len=32):: xx, yy
!read(5,*) nbands
!write(6,*) 'number of bands to be plotted'
!read(5,*) nk
write(6,*) 'name of input file, name of output file'
read(5,*) input,output
open(10,file=input, status='old')
open(11,file=output, status='new')
read(10,*) xx, yy, nbands, yy, nk
allocate(e(nk,nbands))
allocate(k(nk,3))
read(10,*)(k(i,j),j=1,3)
write(6,*)(k(i,j),j=1,3)
read(10,*) (e(i,n),n=1,nbands)
write(13,9030) (e(i,n),n=1,nbands)
9030 format (8f9.5)
do j=1,nbands
if (i.eq.1) then
a=k(i,:)-k0
dk=dk+sqrt(dot_product(a,a))
write(11,*)dk,e(i,j)
write(11,*)
end program prog
这是根据pwscf提供的一个代码所改的。
如何画电荷密度图？
在处理pwscf计算出来的电荷密度，特别是要画面电荷密度分布（等高线图）时，需借助pwscf自带的pp.x和plotrho.x工具。步骤是，先取k点网格进行自洽计算，然后采用pp.x进行数据处理得到在某一个面上的电荷密度值，最后采用plotrho.x来画图。
1)、pp.x在处理总的电荷密度要取出某个面上的电荷密度时的输入文件格式：
prefix = 'si'
outdir = './tmp'
filplot = 'sicharge'
plot_num= 0
spin_component=0
filepp(1) = 'sicharge'
weight(1) = 1.0
output_format = 2
fileout = 'si.rho.dat'
e1(1) =1.0, e1(2)=1.0, e1(3) = 0.0,
e2(1) =0.0, e2(2)=0.0, e2(3) = 1.0,
nx=56, ny=40
这里prefix设置上一步自洽计算电荷密度文件时的体系的名称;
outdir设置了上一步自洽计算中输出文件所在的目录;
filplot设置了要处理的文件的名称，这里我们要处理的是总电荷密度，它是由上一步自洽计算得到的;
plot_num设置了要处理什么类型的数据。这里因为要处理的是电荷密度，因此设置为0;
spin_component设置了是要处理总的电荷密度，还是自旋向上或向下的。设置为0表示总的; 1为自旋向上的电荷密度; 2为自旋向下的电荷密度。
nfile设置要处理的电荷密度文件有几个。这里只有一个，因此后面只设置filepp(1)和weight(1)，如果有2个以上则需指定每一个的文件名及相应的权重。
filepp(1)设置所要处理的第一个文件的文件名; weight(1)是该文件对应的权重。
iflag设置了是要画什么类型的图，因为这是是要画2D等高线图，所以设置为2，表示是要画2D图。
output_format设置处理后的数据的输出方式，这里设置为2,表示是按plotrho.x所需要的格式输出面上的电荷密度值。
fileout是用来输出面电荷密度的文件名。
e1(1), e1(2), e1(3)是用来确定面的第一个矢量，
e2(1), e2(2), e2(3)是用来确定面的第二个矢量，它们的单位是以alat为单位的。
该平面的原点由x0(1),x0(2),x0(3)来确定。
nx，ny分别用来设置平面上网格的大小，nx表示在沿第一个矢量方向的分割数，ny表示在沿第二个矢量方向上的分割数。
2)、用plotrho.x来画某个面上的电荷密度分布的输入文件格式：
si.rho.dat
第一行是用来设置要画的电荷密度的数据，它是由上一步pp.x产生的。
第二行是用来设置输出的ps文件。
第三行是用来设置在画图时，是否对数据按log或线性的来画图，设置为n表示按线性的来画，设置为y表示是按log来画图。
第四行是分别用来设置电荷密度值的最小值和最大值，以及在画等高线时的线的条数（也就是说在最小值和最大值之间要画出多少条的线)。
如何计算金属的功函数？
功函数(electronic work function)定义为使一个电子完全脱离金属表面所需要的能量，在计算中体系的真空能级减去费米能级。有关公函数计算的两篇较好详细的理论文献：
1). PhD thesis of Caspar Fall, Ab initio study of the work functions of elemental metal crystals
http://library.epfl.ch/theses/?nr=1955
2). C J Fall, N Binggeli and A Baldereschi, Deriving accurate work functions from thin-slab calculations, 1999 J. Phys.: Condens. Matter 11 .
http://www.iop.org/EJ/abstract//13/006
其中费米能级在对体系进行自洽计算可以得到。真空能级一般需对静电势（Hatree势和电子交换关联势之和）求微观平均后得到。采用pwscf计算金属的公函数的整个计算的步骤为：
i)&构造金属表面的slab模型，选取合适的原子层数以及真空层数；
ii)&采用pw.x模块对所构造的slab超原胞进行结构优化；
iii)采用pw.x模块对对优化出来的结构进行自洽计算；
iv)采用pp.x模块处理前面自洽计算出来的数据，得到静电势；
v)采用average.x模块求得静电势的微观平均值，并得到真空能级。
例子：计算理想Al(001) (1x1)表面的功函数
在pw.x进行自洽计算的输入文件al001.in:
calculation='scf'
restart_mode='from_scratch',
prefix='Al',
pseudo_dir = './',
outdir='./'
ibrav= 0, nat=11, ntyp= 1,
ecutwfc =16,
occupations='smearing', smearing='methfessel-paxton', degauss=0.01
&electrons
conv_thr = 1.0d-8
mixing_beta = 0.7
CELL_PARAMETERS cubic
0 5.41176 0
0 0 60.9909
ATOMIC_SPECIES
Al 13.867 Al.vbc.UPF
ATOMIC_POSITIONS {angstrom}
Al 1.89 2.025
Al 0 0 4.05
Al 1.89 6.075
Al 0 0 8.1
Al 1.89 10.125
Al 0 0 12.15
Al 1.89 14.175
Al 0 0 16.2
Al 1.89 18.225
Al 0 0 20.25
K_POINTS {automatic}
8 8 1 0 0 0
这里手动自己输入slab模型超原胞的晶格矢量（见CELL_PARAMETERS cubic下面的三行），共11个原子层。
在采用pw.x&计算(命令pw.x &al001.in &al001.out)完后，在该计算目录（outdir='./'所决定的)由中会产生波函数和电荷密度等文件，并产生Al.save目录（命令名由prefix='Al',所决定的）。
然后准备pp.x的输入文件,这里为pp.in，其内容如下：
prefix = 'Al'
outdir = './'
filplot = 'Al.pot'
plot_num= 11
其中prefix用于设置输出文件的前缀名，注意与pw.x的输入文件中的prefix的赋值一致。outdir用来设置输出文件的目录，注意与pw.x输入文件中的outdir的赋值一致。filpot用来设置所产生的主要输出文件的文件名。plot_num用来指定pp.x处理得到什么物理量的数据，这里设置为11，表示处理得到静电势的值。
在运行pp.x &pp.in &pp.out之后，会产生filpot所定义的Al.pot文件。下面就采用pwscf所提供的average.x计算得到静电势的微观平均值，其中average.x的输入文件，这里假设为average.in。此例子中它的输入内容如下：
第一行的值用来设置要处理几个文件，这里设置为1，表示只处理一个文件。所要处理的文件为第二行说定义的，这里为Al.pot（这是上一步由pp.x所产生得到的）。第三行用来设置每个文件所对应的权重（由于此例中，只处理一个文件，也就是一个物理量的值。如果处理多个文件时，也就是多个物理量，要把它们的值按一定的公式进行加或减的操作时，需指定每一个文件所对应的物理量在公式中的权重。比如所要处理A-B，那么A所对应的权重就是1，B所对应的权重就是-1）。
第四行的值用来设置要输出静电势微观平均值在沿某个方向上要输出多少个点的值。第五行用来设置在哪个方向上求静电势的微观平均值，这里设置为3表示在z轴方向。1和2分别对应于x和y轴。第六行用来设置计算微观平均值时的距离宽度，个人认为可以设置为原子的层间距。
通过在pw.x自洽计算得到的al001.out文件中搜'Fermi energy'字符串可以看到，本例中Al(001)(1x1)理想表面的费米能级值为3.9470 eV。通过average.x计算静电势的微观平均得到真空能级为8.376 eV，因此它的功函数为4.4287 eV。
如何用PWscf进行STM模拟？
pwscf是按Tersoff-Hamann近似来进行STM模拟的。当前版本中，不能对“晶层模型“(slab model)具有两个对称表面的情况进行处理。在进行STM模拟的步骤：1)，对所模拟的slab模型进行结构优化pw.x；2)、对优化后的结构进行自洽计算pw.x；3)、增加k点网格的大小，并进行非自洽计算pw.x；4)、将计算完后的结果采用pp.x进行处理；5)、采用plotrho.x进行画图处理。
下面介绍在4）和5）步中的输入文件及其说明：
将计算完后的结果采用pp.x进行处理以得到STM图象的数据时的输入文件：
prefix = 'AlAs110'
outdir='./tmp',
filplot = 'AlAsresm-1.0'
sample_bias=-0.0735d0,
stm_wfc_matching=.false.,
plot_num= 5
filepp(1)='AlAsresm-1.0'
weight(1)=1.0
output_format=2
e1(1)=7.0, e1(2)=0.0, e1(3)=0.0
e2(1)=0.0, e2(2)=7.07107, e2(3)=0.0
x0(1)=0.0, x0(2)=-0.18, x0(3)=3.25
nx=36 ,ny=56
fileout='AlAs110-1.0'
在用pp.x进行STM的数据处理，其输入文件中将plot_num设置5，另外需设置sample_bias，即所加的偏压大小；以及是否考虑STM模拟中波函数的匹配问题（即考虑波函数从某个位置开始衰减为零），如果把stm_wfc_matching设置为.false.，也就是不考虑匹配问题。其他的关键词的含义与对总的电荷密度进行处理（画等高线2D图）是一样的。
采用plotrho.x进行画图处理
AlAs110-1.0
AlAs110-1.0eV.ps
它的输入参数的含义与在采用plotrho.x画电荷密度等高线时的输入文件参数的含义是一样的。
第一行是用pp.x得到在某个偏压下的在某个面上的电荷密度分布；
第二行是plotrho.x画图所要得到的图的文件名；
第三行是表示所处理电荷密度是按log&（设置为y)还是线性方式（设置为n)输出画图的；
第四行是用来设置面上的电荷密度的最小值、最大值和这些极值之间的等高线的条数。
转载本文请联系原作者获取授权，同时请注明本文来自穆跃文科学网博客。链接地址：
上一篇：下一篇：
当前推荐数：0
评论 ( 个评论)
扫一扫，分享此博文
作者的其他最新博文
热门博文导读
Powered by
Copyright &