常用PWscf帮助文档及网址
User Guide
官方用户指南,PWscf的编译、安装、使用、常见问题等信息非常全面,通俗易懂。初学者必读。结合examples使用
Examples
PWscf安装目录下的Examples包含30多个常见实例,直接将run_example的shell脚本拖入终端即可运行。运行后在results文件夹下可以看对应输入文件和输出结果。可以参照示例写自己的输入文件。
Doc
PWscf安装目录下的Doc文件夹内有许多html文件,内容包括各种输入文件中所有变量和标签的详细定义和规范,在写输入文件的时候一定要仔细阅读。
PWscf邮件列表
在Google里搜索PWscf的问题排在第一的往往就是PWscf官方的邮件列表。这里汇集了全球PWscf用户的常见问题及解答。你的问题往往就在里面。对于常见的小问题这里的答案往往简单明了,事半功倍。有需要的同学可以订阅。
大学及研究所的课件、讲义、ppt等 (以下译自PWscf的wiki )
关于理论背景和计算实例的一些ppt
伊朗某大学关于计算材料学的课程的一些lectures
,
MIT的课件
,
北卡州立大学的课件
UIUC计算材料中心2006年夏令营的课件和视频
康奈尔大学纳米科技中心2006年秋季研讨会的讲座和指南,包括平面波赝势法、密度泛函理论、分子动力学、纳米光学等方面的资料
密度泛函理论和实例
如何安装编译PWscf?
单机安装简易版(适用Ubuntu32位操作系统)from
首先安装编译环境,确认联网并,在终端下输入:
sudo apt-get install gfortran
sudo apt-get install mpich-bin
然后去PWscf网站下载免费的源码包
解压到任意目录
在终端窗口进入解压目录,方法如下:
输入: cd 路径 (可直接把文件夹拖拽进终端自动生成路径)
然后输入
./configure
等屏幕上一堆字刷完再输入
make all
继续等刷屏 大概十分钟以后应该就好了
测试软件 进入解压目录下的examples文件夹
看到有很多example01 example02等目录
随便进一个 把run_example文件拖进终端窗口
看运行是否正常
也可以运行examples文件夹下的run_all_examples文件
运行所有的示例
计算结果在各自文件夹下的results目录里
要清除计算结果 运行examples下的make clean文件
本教程应该适用于32位Ubuntu为基础的所有衍生版
64位Ubuntu源里好像没有mpich,需要自己
一个台湾人写的安装教程(适用于Ubuntu8.10)
pwscf - ubuntu 8.10 安裝筆記
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因為學弟的電腦是 Ubuntu 8.10 ,我也不知道在其他 distribution 可不可以依照相同的步驟安裝;這裡列出一些套件,找相對應的裝應該就差不多了吧。
the atomic scale. PWscf is released under the . distribution of codes for the quantum simulation of matter at
下載
PWscf
- PWscf:
- PWGui:
安裝 pwscf
1. 預先安裝套件,請安裝下列套件,或 sudo apt-get install 套件名稱 即可
- 編譯環境:build-essential
- FFT: fftw-double-dev fftw-single-dev
- MPI: libmpich1.0-dev
- SSH: ssh
- Fortran: gfortran
基本步驟摘要 (只用於一台):設置多台 ssh 可進入之免密碼機器組
- 打 cd ~/.ssh(注意 ssh 前面有一個點)進入使用者個人隱藏的 ssh 設定目錄(在 UNIX 下凡是名稱以一點作為開始的檔案及目錄,都會隱藏,只用 ls 看不見,要用 ls -a 才看得見)。
- 打 ssh-keygen 產生 ssh 的鑰匙,它會問要輸入 passphrase,給不給都可以,之後你在目錄下看會到 id_rsa 及 id_rsa.pub 兩個檔案。
- 把 id_rsa.pub 拷貝為檔名叫作 authorized_keys 的檔案。
- ~/.ssh 底下的 id_rsa、id_rsa.pub、authorized_keys 等檔案,都要 chmod 成為群組 (g) 與其他 (o) 皆不可讀與寫。chmod g-w [filename] 、 chmod g-r [filename] ,比較快的方法是 chmod 500 [filename]。
- 現在可以試著以 ssh username@hostname (遠端)或 ssh localhost (本機)試著用 ssh 以不需要密碼的方式登入遠端或本機系統,如果不需要密碼可以登入,則表示設定成功。
觀念是,遠端的那台機器上之使用者的 .ssh 子目錄下之 authorized_keys 檔案之中要有記載想免密碼近端機器之公鑰匙。
按上面的方法,先在近端機器的使用者目錄下產生 .ssh 子目錄以及其下相關檔案。
把其中 id_rsa.pub 傳送到遠端,並利用 cat id_rsa.pub >> authorized_keys(注意是兩個大於號,代表追加)來把其內容加入到遠端機器的 authorized_keys 檔案中。
測試從近端到遠端以 ssh 登入,看是否需要密碼。
3. 打開終端機,cd espresso-4.0.4
4. ./configure
5. make all
安裝 pwgui
6. tar zxvf PWgui-4.1CVS.tgz
7. cd PWgui-4.1CVS
8. 執行PWGui: ./pwgui
註:後來發現 pwscf 的目錄裡本來就帶有 pwgui 了,或許可以直接執行 ./pwgui 測試看看可不可以使用。
設定 pwgui
設定 pwgui 可執行檔的路徑如下圖,此例中 pwscf 被解壓縮在這個目錄下:
- /home/labpc/PWgui/espreso-4.0.4/
測試 pwscf
pwscf 已經提供非常多 example 在原始碼裡面,在我們的路徑下為:
- /home/labpc/PWgui/espreso-4.0.4/examples/
- cd /home/labpc/PWgui/espreso-4.0.4/ (此為解開之後的路徑)
- cd example01
- ./run_example
by Sunforever如何从MS中导出原子坐标?
如何进行结构优化?
转自:http://old.blog.edu.cn/user1/11542/archives/2006/1578700.shtml
在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分子
&CONTROL
calculation = "relax",
prefix = "CO",
pseudo_dir = "./",
outdir = "./",
tprnfor=.true.
forc_conv_thr=1.0D-4
/
&SYSTEM
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,
/
&IONS
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的自洽计算
&control
calculation='scf'
restart_mode='from_scratch',
pseudo_dir = './',
outdir='./'
prefix='cu'
tstress = .true.
tprnfor = .true.
/
&system
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)控制计算的部分,也就是要设置
&control
.....
第一个'/'之间的关键词。
关键词calculation赋值为'scf'表示此计算是进行自洽电荷密度计算;
restart_mode表示是否是接着上一次的计算而继续的计算,赋值为'from_scratch'意味着是进行一次全新的计算开始;
pseudo_dir用来设置赝势文件所在的目录,赋值为'./'表示赝势文件放在当前计算目录;
outdir用来设置计算过程中输出文件(比如波函数、电荷密度以及势)输出到哪个目录中。赋值为'./'表示这些输出文件将放到当前计算目录中;
prefix用来定义当前计算作业的标题名,它将是一些主要输出文件的文件名。赋值为'cu'用来标记当前计算作业是对Cu进行计算;
tstress 用来设置在自洽计算过程中是否计算体系的应力,设置为 .true.表示在自洽计算过程中要计算体系的应力;
tprnfor 用来设置在自洽计算过程中是否计算体系中原子所受的力,设置为 .true.表示在自洽计算过程中要计算体系中原子所受的力;
2)、描述所计算的体系(包括它的晶格类型、晶格常数或结构参数、原胞基矢、原胞中原子的类型数目和总的原子数目)、平面波的切断动能(也就是在展开KS轨道或晶体波函数的平面波切断动能;另外,还包括在计算电荷密度时,展开的平面波的切断动能)、确定电子占有数的方法及相关的参数。也就是由
&system
..........
第二'/'之间的关键词来设置。
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点进行平移。
如何计算能带结构?
转自:http://old.blog.edu.cn/user1/11542/archives/2006/1573345.shtml
计算fcc Cu的能带结构
&control
calculation='bands'
pseudo_dir = './',
outdir='./',
prefix='cu'
/
&system
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
K_POINTS
28
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)
c
open(5,file='syml',status='old')
open(7,file='inp.kpt')
c
read(5,*) nhighk
read(5,*) (ndiv(i),i=1,nhighk-1)
do i=1,nhighk-1
ntkp=ntkp+ndiv(i)
enddo
ntotkpt=ntkp+1
if(nhighk>10)then
write(*,*)'Number of high-symmetry k points must < 10!'
STOP
endif
if(ntotkpt>200)then
write(*,*)'Total number of k points must <= 200!'
STOP
endif
do i=1, nhighk
read(5,*) labhk(i),(phighk(i,j),j=1,3)
enddo
write(*,*) (labhk(i),i=1,nhighk)
c
c----- generating k-points along high symmetric lines --------
c
c
pk(1,1)=phighk(1,1)
pk(1,2)=phighk(1,2)
pk(1,3)=phighk(1,3)
ii=1
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
enddo
enddo
c
10 format(A34)
weight=1.d0
do i=1,ntotkpt
write(7,200) pk(i,1),pk(i,2),pk(i,3),weight
enddo
200 format(3F10.6,F6.2)
stop
end
c----------------------- end ---------------------------
它的输入文件为syml,输出文件为inp.kpt。其中syml输入文件的格式如下:
8
15 15 15 15 15 15 15
G 0.0 0.0 0.0
K -0.33333333333 0.6666666666667 0.000000000
H -0.33333333333 0.6666666666667 0.500000000
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在计算能带时的输入文件中。
如何画能带图?
转自:http://old.blog.edu.cn/user1/11542/archives/2007/1651378.shtml
pwscf 附带了band.x和plotband.x的工具,前者是将计算出来k点坐标以及相应的本征值从out文件中收集起来(或取出)专门存储到一个文件中,以便后一个工具plotband.x进行处理。在计算能带时,先设置k点网格进行一次自洽计算,然后自己输入要计算的特殊k点并进行一次非自洽计算,得到这些特殊k点的本征值。再用band.x和ploband.x进行处理。
1). band.x的输入文件格式
&inputpp
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 0.500000 0.500000
-3.418 -0.822 5.029 5.029 7.814 9.597 9.597 13.838
0.400000 0.400000 0.400000
-3.891 -0.102 5.102 5.102 7.900 9.679 9.679 13.959
0.300000 0.300000 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
-6.0 10
sibands.xmgr
sibands.ps
6.255
1.0 6.255
第一行是band.x处理得到的k点本征值文件 此例子中是sibands.dat;
第二行是在这一步中输出的ps文件中纵坐标(本征值能量)刻度的最小值与最大值;
第三行是用来设置所出输出xmgr格式的文件的文件名;
第四行是用来设置所输出文件ps格式的文件名;
第五行是费米能级的值,这个在自洽计算的out文件中可以找到;
第六行中第一个数是用来设置输出 ps格式文件纵坐标(能量)刻度的大小,第二个数是用来设置标出费米能级的位置,它与第五行中的数相同。
同样ploband.x也可以用来处理声子本征值文件。
如何计算态密度?
转自:http://old.blog.edu.cn/user1/11542/archives/2006/1575092.shtml
例子,计算Cu的态密度
一、自洽计算(见前一个例子)
二、非自洽计算,增加k点,并采用四面体方法来确定电子的占有数
&control
calculation='nscf'
prefix='cu',
pseudo_dir = './',
outdir='./'
/
&system
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计算总态密度
&inputpp
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计算分波态密度。
在自洽计算中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用来设置展宽函数中的展宽系数。
如何计算晶体的声子色散曲线和态密度?
转自:http://old.blog.edu.cn/user1/11542/archives/2007/1647197.shtml
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进行电子密度的自洽计算
&control
title='Sc, hexagonal cell'
calculation = 'scf'
restart_mode='from_scratch',
prefix='sc',
pseudo_dir = './',
outdir='./tmp'
tprnfor=.true.
/
&system
ibrav=4,
celldm(1)=6.05606,
celldm(3)=1.71298,
nat=2,
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.3333333333333286 0.6666666666666714 0.2500000000000000
Sc 0.6666666666666714 0.3333333333333286 0.7500000000000000
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
&inputph
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计算实空间力常数矩阵
&input
zasr='simple', fildyn='sc.dyn', flfrc='Sc444.fc', la2F=.false.
/
在q2r.x的输入文件中需指定
fildyn: 用来设置包含了q网格点的动力学矩阵元的文件,与上一步的 fildyn设置一致;
flfrc:用来设置输出力常数矩阵的文件;
la2F:用来设置是否计算出实空间中电-声耦合系数;针对计算材料的超导性质;
zasr:如何处理‘声学支求和规则“,该规则是用在处理Born有效电荷的,要求Born有效电荷的总和是零。可赋的值有:
no,表示不处理声学支求和问题
simple, 表示通过对力常数矩阵的对角元素进行修正来考虑3支声学横模的求和处理;
crystal,
one-dim,
zero-dim
这里我们设置为simple。
4)用matdyn.x计算出声子色散曲线
&input
asr='simple', amass(1)=44.955910,
flfrc='Sc444.fc', flfrq='Sc444.freq', la2F=.false., dos=.false.
/
131
0.000000 0.000000 0.000000 0.00
0.000000 0.016667 0.000000 0.00
0.000000 0.033333 0.000000 0.00
0.000000 0.050000 0.000000 0.00
0.000000 0.066667 0.000000 0.00
0.000000 0.083333 0.000000 0.00
0.000000 0.100000 0.000000 0.00
0.000000 0.116667 0.000000 0.00
0.000000 0.133333 0.000000 0.00
0.000000 0.150000 0.000000 0.00
0.000000 0.166667 0.000000 0.00
0.000000 0.183333 0.000000 0.00
0.000000 0.200000 0.000000 0.00
0.000000 0.216667 0.000000 0.00
0.000000 0.233333 0.000000 0.00
0.000000 0.250000 0.000000 0.00
0.000000 0.266667 0.000000 0.00
....
这里要输入131个特殊q点的坐标。与计算能带结构时一样,需先选出要计算的高对称q点的走向以及高对称点的坐标,然后产生这些线上的q点的坐标。
计算出来的每一个q点的本征频率可按上一个blog中提到的方法处理一下后画图。
5),用matdyn.x计算声子态密度
&input
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来设置。
如何画出声子色散曲线?
转自:http://old.blog.edu.cn/user1/11542/archives/2007/1645444.shtml
在采用pwscf中的ph.x以及matdyn.x计算得到了沿高对称q点的本征频率后,需将数据转换一下,以画出声子色散曲线。下面是计算出来的本征频率文件中的格式:
&plot nbnd= 6, nks= 91 /
0.500000 0.000000 0.500000
139.6336 194.6691 277.0816 409.1945 460.9638 475.5246
0.475000 0.000000 0.475000
136.1839 189.6126 269.1646 418.0272 463.4473 475.7502
0.450000 0.000000 0.450000
132.4453 184.0388 260.9538 426.7985 466.0779 475.9910
0.425000 0.000000 0.425000
128.4070 177.8615 252.3483 435.3956 468.8332 476.3233
0.400000 0.000000 0.400000
124.0518 171.0092 243.2433 443.7172 471.6875 476.8231
0.375000 0.000000 0.375000
119.3569 163.4369 233.5371 451.6723 474.6123 477.5588
0.350000 0.000000 0.350000
114.2958 155.1312 223.1375 459.1796 477.5767 478.5845
0.325000 0.000000 0.325000
.......................
第一行中的nbnds告诉了有每个q点有多少个本征振动频率,nks告诉了计算了多少个q点。下面的接着的每两行是:
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))
do i=1,nk
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)
enddo
do j=1,nbands
dk=0
do i=1,nk
if (i.eq.1) then
k0=k(i,:)
endif
a=k(i,:)-k0
dk=dk+sqrt(dot_product(a,a))
write(11,*)dk,e(i,j)
k0=k(i,:)
enddo
write(11,*)
enddo
stop
end program prog
--------
这是根据pwscf提供的一个代码所改的。
如何画电荷密度图?
转自:http://old.blog.edu.cn/user1/11542/archives/2007/1649978.shtml
在处理pwscf计算出来的电荷密度,特别是要画面电荷密度分布(等高线图)时,需借助pwscf自带的pp.x和plotrho.x工具。步骤是,先取k点网格进行自洽计算,然后采用pp.x进行数据处理得到在某一个面上的电荷密度值,最后采用plotrho.x来画图。
1)、pp.x在处理总的电荷密度要取出某个面上的电荷密度时的输入文件格式:
&inputpp
prefix = 'si'
outdir = './tmp'
filplot = 'sicharge'
plot_num= 0
spin_component=0
/
&plot
nfile = 1
filepp(1) = 'sicharge'
weight(1) = 1.0
iflag = 2
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
si.rho.ps
n
0 0.09 6
第一行是用来设置要画的电荷密度的数据,它是由上一步pp.x产生的。
第二行是用来设置输出的ps文件。
第三行是用来设置在画图时,是否对数据按log或线性的来画图,设置为n表示按线性的来画,设置为y表示是按log来画图。
第四行是分别用来设置电荷密度值的最小值和最大值,以及在画等高线时的线的条数(也就是说在最小值和最大值之间要画出多少条的线)。
如何计算金属的功函数?
转自:http://old.blog.edu.cn/user1/11542/archives/2007/1637193.shtml
功函数(electronic work function)定义为使一个电子完全脱离金属表面所需要的能量,在计算中体系的真空能级减去费米能级。有关公函数计算的两篇较好详细的理论文献:
1). PhD thesis of Caspar Fall, Ab initio study of the work functions of elemental metal crystals
2). C J Fall, N Binggeli and A Baldereschi, Deriving accurate work functions from thin-slab calculations, 1999 J. Phys.: Condens. Matter 11 2689-2696.
其中费米能级在对体系进行自洽计算可以得到。真空能级一般需对静电势(Hatree势和电子交换关联势之和)求微观平均后得到。采用pwscf计算金属的公函数的整个计算的步骤为:
i) 构造金属表面的slab模型,选取合适的原子层数以及真空层数;
ii) 采用pw.x模块对所构造的slab超原胞进行结构优化;
iii)采用pw.x模块对对优化出来的结构进行自洽计算;
iv)采用pp.x模块处理前面自洽计算出来的数据,得到静电势;
v)采用average.x模块求得静电势的微观平均值,并得到真空能级。
例子:计算理想Al(001) (1x1)表面的功函数
在pw.x进行自洽计算的输入文件al001.in:
&control
calculation='scf'
restart_mode='from_scratch',
prefix='Al',
pseudo_dir = './',
outdir='./'
/
&system
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
5.41176 0 0
0 5.41176 0
0 0 60.9909
ATOMIC_SPECIES
Al 13.867 Al.vbc.UPF
ATOMIC_POSITIONS {angstrom}
Al 0 0 0
Al 1.43189 1.43189 2.025
Al 0 0 4.05
Al 1.43189 1.43189 6.075
Al 0 0 8.1
Al 1.43189 1.43189 10.125
Al 0 0 12.15
Al 1.43189 1.43189 14.175
Al 0 0 16.2
Al 1.43189 1.43189 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
然后准备pp.x的输入文件,这里为pp.in,其内容如下:
&inputpp
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
1
Al.pot
1.D0
5000
3
2.95
第一行的值用来设置要处理几个文件,这里设置为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模拟?
转自:http://old.blog.edu.cn/user1/11542/archives/2007/1650357.shtml
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图象的数据时的输入文件:
&inputpp
prefix = 'AlAs110'
outdir='./tmp',
filplot = 'AlAsresm-1.0'
sample_bias=-0.0735d0,
stm_wfc_matching=.false.,
plot_num= 5
/
&plot
nfile=1
filepp(1)='AlAsresm-1.0'
weight(1)=1.0
iflag=2
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
n
0.00005 0.0078 8
它的输入参数的含义与在采用plotrho.x画电荷密度等高线时的输入文件参数的含义是一样的。
第一行是用pp.x得到在某个偏压下的在某个面上的电荷密度分布;
第二行是plotrho.x画图所要得到的图的文件名;
第三行是表示所处理电荷密度是按log (设置为y)还是线性方式(设置为n)输出画图的;
第四行是用来设置面上的电荷密度的最小值、最大值和这些极值之间的等高线的条数。