PWSCF教程从入门到精通.docx

《PWSCF教程从入门到精通.docx》由会员分享，可在线阅读，更多相关《PWSCF教程从入门到精通.docx（41页珍藏版）》请在taowenge.com淘文阁网|工程机械CAD图纸|机械工程制图|CAD装配图下载|SolidWorks_CaTia_CAD_UG_PROE_设计图分享下载上搜索。

1、PWSCF计算实例分析（1） 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. tprnf

2、or = .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

3、 K_POINTS (automatic) 8 8 8 0 0 0 解释: 在电子自洽计算中需设置以下几个方面的参数： 1）控制计算的部分，也就是要设置 &control . 第一个/之间的关键词。 关键词calculation赋值为scf表示此计算是进行自洽电荷密度计算； restart_mode表示是否是接着上一次的计算而继续的计算，赋值为from_scratch意味着是进行一次全新的计算开始； pseudo_dir用来设置赝势文件所在的目录，赋值为./表示赝势文件放在当前计算目录； outdir用来设置计算过程中输出文件（比如波函数、电荷密度以及势）输出到哪个目录中。赋值为./表示这些输

4、出文件将放到当前计算目录中； prefix用来定义当前计算作业的标题名，它将是一些主要输出文件的文件名。赋值为cu用来标记当前计算作业是对Cu进行计算； tstress 用来设置在自洽计算过程中是否计算体系的应力，设置为 .true.表示在自洽计算过程中要计算体系的应力； tprnfor 用来设置在自洽计算过程中是否计算体系中原子所受的力，设置为 .true.表示在自洽计算过程中要计算体系中原子所受的力; 2)、描述所计算的体系（包括它的晶格类型、晶格常数或结构参数、原胞基矢、原胞中原子的类型数目和总的原子数目）、平面波的切断动能（也就是在展开KS轨道或晶体波函数的平面波切断动能；另外，还包括

5、在计算电荷密度时，展开的平面波的切断动能）、确定电子占有数的方法及相关的参数。也就是由 &system . 第二/之间的关键词来设置。 ibrav用来归属体系所属的晶格类型，赋值为2表示所计算的体系是fcc结构； celldm(1)用来设置体系的第一个晶格常数，因为所计算的体系是fcc结构，只需设置celldm(1)，相当于指定晶格常数a的值； nat用来指明体系的原胞中原子的总共数目，赋值为1表示所计算的原胞中只有一个原子； ntyp用来指明体系中原子类型的数目，赋值为1表示所计算的体系只有一种类型的原子； occupations用来设置确定电子占有数的方法，赋值为smearing表示采用s

6、mearing的方法来确定电子的占有数，随后须设置smearing和degauss关键词； smearing用来指明确定电子占有数的一种具体的smearing方法，赋值为gaussian表示采用Gaussian函数来确定电子占有数； degauss用来确定smearing方法中有关函数的展宽参数，赋值为0.02表示上面Gaussian函数中的展宽参数为0.02。 3）、设置电子自洽计算中本征矢量（波函数）和本征值的计算算法，自洽收敛的标准。也就是 &electrons . 和第三个/之间的关键词来设置。 diagonalization用来设置在求KS方程的本征矢量和本征值时，采用具体的什么算法

7、，赋值为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

8、-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点进行平移。在PWSCF计算实例分析(2)-能带结构计算 计算fcc Cu的能带结构&controlcalculation=b

9、andspseudo_dir = ./,outdir=./,prefix=cu/&systemibrav = 2, celldm(1) =6.73, nat= 1, ntyp= 1,ecutwfc = 25.0, ecutrho = 300.0, nbnd = 8/&electronsdiagonalization=david/ATOMIC_SPECIESCu 63.55 Cu.pz-d-rrkjus.UPFATOMIC_POSITIONSCu 0.0 0.0 0.0K_POINTS280.0 0.0 0.0 1.00.0 0.0 0.1 1.00.0 0.0 0.2 1.00.0 0.0 0

10、.3 1.00.0 0.0 0.4 1.00.0 0.0 0.5 1.00.0 0.0 0.6 1.00.0 0.0 0.7 1.00.0 0.0 0.8 1.00.0 0.0 0.9 1.00.0 0.0 1.0 1.00.0 0.0 0.0 1.00.0 0.1 0.1 1.00.0 0.2 0.2 1.00.0 0.3 0.3 1.00.0 0.4 0.4 1.00.0 0.5 0.5 1.00.0 0.6 0.6 1.00.0 0.7 0.7 1.00.0 0.8 0.8 1.00.0 0.9 0.9 1.00.0 1.0 1.0 1.00.0 0.0 0.0 1.00.1 0.1 0

11、.1 1.00.2 0.2 0.2 1.00.3 0.3 0.3 1.00.4 0.4 0.4 1.00.5 0.5 0.5 1.0解释：在进行能带计算时，calculation须设置为bands，而且在此之前须进行一次相应的自洽计算，而且要有上一步计算得到输出文件供能带计算时读入。另外最好在&system中设置nbnd，以指定计算多少条能带。在计算能带时要自己先选定一些高对称点，并产生这些高对称点之间其他点。在这个例子中，计算沿G-X-L点之间的高对称线上的能带。在产生所要计算的特殊k点时，可以采用下面简单的f77程序来实现：c +-c For generating k-points alo

12、ng the high-symmetry lines inc Brillouin zone and for calculate band-structures !c +-C -syml-c 6 : nhighkc 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.0c L 0.5 0.5 0.5c W 0.5 0.25 0.75c K 0.375 0.375 0.75c G 0.0 0.0 0.0c direct & reciprocal lattice vectors over emin

13、, emax lineC -c max k-points = 200program gkimplicit real*8 (a-h,o-z)character*2 labhkdimension tkpt(200,3),pk(200,3),phighk(10,3)dimension disk(200),dish(10),labhk(10)dimension ndiv(10)copen(5,file=syml,status=old)open(7,file=inp.kpt)cread(5,*) nhighkread(5,*) (ndiv(i),i=1,nhighk-1)do i=1,nhighk-1n

14、tkp=ntkp+ndiv(i)enddontotkpt=ntkp+1 if(nhighk10)thenwrite(*,*)Number of high-symmetry k points must 200)thenwrite(*,*)Total number of k points must kvecs_FS.in kvecs_FS.in kvecs_FS.out kvecs_FS.x kvecs_FS.out注释：这里$Sysname用来标记所计算的体系，这个例子中给它赋值为cu，$Calc_Type用来标记所进行的是Fermi surface的计算，可以赋值为FS，$nabc用来设置在计

15、算费米面时所要计算的k点网格，这里设置为12x12x12的k点网格。$n_start用来标记从哪个能带开始计算起，这里设置为3，$n_last用来设置要计算到哪根能带为止，这里设置为6，一般它们两个要设置在费米能级所在的带的上下，也就是它们的范围之内要包括了费米能级。$a1, $a2,$a3用来标记所得到k点网格产生的k点坐标。$E_Fermi用来标记在自洽计算中得到的费米能级值。 二、进行非自洽的计算，得到这些k点的本征值：cat cu.fs.in cu.fs.inpw.x cu.fs.out 三、准备将计算的本征值和对应的k点写成xcrysden的bxsf格式mv cu.fs.out Ba

16、nds.outcat input_FS bands_fs.outmv Bands.bxsf cu.fs.bxsfPWSCF计算实例分析计算Gamma点声子频率 PWSCF的ph.x模块是基于linear response理论或密度泛函微扰理论来计算声子频率。计算步骤：1、先用pw.x进行自洽计算2、采用ph.x计算Gamma点声子频率一、自洽计算&control calculation=scf, restart_mode=from_scratch, prefix=si pseudo_dir = ./, outdir=./&system ibrav = 2, celldm(1) =10.20,

17、nat= 2, ntyp= 1, ecutwfc = 18.0/&electrons mixing_beta = 0.7 conv_thr = 1.0d-8/ATOMIC_SPECIESSi 28.086 Si.vbc.UPFATOMIC_POSITIONSSi 0.00 0.00 0.00 Si 0.25 0.25 0.25K_POINTS (automatic)8 8 8 0 0 0 pw.x si.scf.out上面是对fcc Si进行自洽的电子结构计算，输入文件与前面对Cu的计算类似。二、用ph.x进行Gamma点声子频率的计算输入文件Phonon of Si at Gamma&inp

18、utph tr2_ph=1.0d-14, prefix=si, epsil=.true., amass(1)=28.08, outdir=./,fildyn=si.dynG,/0.0 0.0 0.0ph.x si.phG.out上面第一行是注释行，关键词由&inputph ./给出，在/之后给出的q点Gamma点的坐标。tr2_ph用来设置声子计算时自洽收敛的标准；prefix与前面自洽计算中的关键词prefix一致，用来标记所计算的体系，注意它的值要与自洽计算中的一致；epsil在q=0（也就是Gamma点时），如果体系半导体，而且epsil设置为.true.，则表示计算半导体的宏观介电常数

19、。如果是金属性的体系，或者非Gamma的计算，则不能设置为.true.；amass(1)用来设置体系中第一类原子的原子量；outdir用来设置计算输出的目录；fildyn用来设置指定要将动力学矩阵元输出到什么文件中。PWSCF计算实例分析原子位置优化在pwscf中提供了两种优化方法对原子位置进行驰豫：a) BFGS quasi-newton algorithm, b)最速下降法（quick-min Verlet）。除非初始位置很接近平衡位置，一般采用BFGS准牛顿算法比较快。在结构优化时，calculation需设置为relax，下面相关的参数有时也需要设置： 一、在&control ./中设

20、置优化的收敛标准参数、步数等 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./中设置优化方法中的相关参数 在calcul

21、ation=relax时，ion_dynamics设置为bfgs表示用BFGS准牛顿算法来进行优化，设置为damp表示用最速下降法来进行优化。 pot_extrapolation：用来控制优化或电子迭代过程中势的混合方式，在原子位置优化时，最好设置为second_order，表示采用二阶方式来进行混合； wfc_extrapolation:用来控制优化或电子迭代过程中波函数的混合方式，在原子位置优化时，最好设置为second_order，表示采用二阶方式来进行混合； 当设置了采用BFGS准牛顿算法来优化后，下面的参数需设置：upscale：用来控制conv_thr在结构优化过程中最可能的缩小值

22、，在优化快接近收敛时，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;trus

23、t_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

24、_thr = 1.D-7, mixing_beta = 0.7D0,/&IONS ion_dynamics=bfgs,pot_extrapolation = second_order, wfc_extrapolation = second_order,/ATOMIC_SPECIESO 1.00 O.LDA.US.RRKJ3.UPFC 1.00 C.pz-rrkjus.UPFATOMIC_POSITIONS bohrC 2.256 0.0 0.0O 0.000 0.0 0.0 0 0 0 K_POINTS Gamma 如何采用PWSCF计算金属的功函数 http:/library.epfl.ch

25、/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 2689-2696.http:/www.iop.org/EJ/abstract/0953-8984/11/13/006其中费米能级在对体系进行自洽计算可以得到。真空能级一般需对静电势（Hatree势和电子交换关联势之和）求微观平均后得到。采用pwscf计算金属的公函数的整个计算的步骤为

26、：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=./&

27、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 cubic5.41176 0 00 5.41176 00 0 60.9909 ATOMIC_SPECIESAl 13.867 Al.vbc.UPF ATOMIC_POSITIONS angstromAl 0 0 0Al 1.43189 1.43189 2

28、.025Al 0 0 4.05Al 1.43189 1.43189 6.075Al 0 0 8.1Al 1.43189 1.43189 10.125Al 0 0 12.15Al 1.43189 1.43189 14.175Al 0 0 16.2Al 1.43189 1.43189 18.225Al 0 0 20.25K_POINTS automatic 8 8 1 0 0 0这里手动自己输入slab模型超原胞的晶格矢量（见CELL_PARAMETERS cubic下面的三行），共11个原子层。在采用pw.x 计算(命令pw.x al001.out)完后，在该计算目录（outdir=./所决定的

29、)由中会产生波函数和电荷密度等文件，并产生Al.save目录（命令名由prefix=Al,所决定的）。然后准备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，表示处理

30、得到静电势的值。在运行pp.x pp.out之后，会产生filpot所定义的Al.pot文件。下面就采用pwscf所提供的average.x计算得到静电势的微观平均值，其中average.x的输入文件，这里假设为average.in。此例子中它的输入内容如下：1Al.pot1.D0500032.95第一行的值用来设置要处理几个文件，这里设置为1，表示只处理一个文件。所要处理的文件为第二行说定义的，这里为Al.pot（这是上一步由pp.x所产生得到的）。第三行用来设置每个文件所对应的权重（由于此例中，只处理一个文件，也就是一个物理量的值。如果处理多个文件时，也就是多个物理量，要把它们的值按一定的

31、公式进行加或减的操作时，需指定每一个文件所对应的物理量在公式中的权重。比如所要处理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，因此它的功函数为处理PWSCF中计算得到本征频率 处理PWSCF中计算得到本征频率 在采用pwscf中的ph.x