cndaqiang Web Linux DFT

Learn QE1 电子结构计算

2020-03-06
cndaqiang
 
RSS

本学习教程内容主要来自互联网,个人学习记录,仅供参考。
代码仓库LearnQE@cndaqiang

参考

Summer School on Advanced Materials and Molecular Modelling
GitLab for qe2019

注意

  • 每个模块下面都有Example和Manual,如PP/Doc/user_guide.pdf
  • 默认赝势路径设置export ESPRESSO_PSEUDO=/home/cndaqiang/ONCVPSP/abinit
  • 默认outdirexport ESPRESSO_TMPDIR=
  • USPP, it is typically needed to set ecutrho > 4×ecutwfc (it should be at least 8 to 12 times larger)

PWscf

  • Molecule/Cluster 体系要设置assume_isolated参数
  • k点较多(>100)时,设置verbosity='high'输出更多信息
  • 同一元素有多种磁矩时,要使用不同的标签(Fe1,Fe2),因为初始磁矩和ntyp对应,starting_magnetization(1:ntyp)

PP模块

可执行程序

cndaqiang@girl:~/code/tdpw2.0/q-e-qe-6.4.1-devel-1.0/PP/src$ ls *.x
average.x         fs.x             plotband.x  projwfc.x       wannier_ham.x
bands.x           initial_state.x  plotproj.x  pw2bgw.x        wannier_plot.x
dos.x             molecularpdos.x  plotrho.x   pw2critic.x     wfck2r.x
epsilon.x         open_grid.x      pmw.x       pw2gw.x
fermi_proj.x      pawplot.x        ppacf.x     pw2wannier90.x
fermi_velocity.x  plan_avg.x       pp.x        sumpdos.x

pp.x,后处理pw.x和cp.x的计算结果,计算一些物理性质用于画图,图形有1D,2D,3D类型,可以指定输出格式(如xsf),输出可以被下列程序处理

  • Gnuplot (1D or 2D plots),
  • by code plotrho.x that comes with PostProc and produces PostScript 2D plots,
  • by advanced plotting software XCrySDen and gOpenMol (3D plots).

也可以输出其他运算(如pw.x -> pp.x -> average.x)需要的输入
处理的性质有

  • charge density
  • spin polarization
  • various potentials
  • local density of states at EF
  • local density of electronic entropy
  • STM images
  • selected squared wavefunction
  • ELF (electron localization function)
  • RDG (reduced density gradient)
  • integrated local density of states

pwscf.save文件夹

  • scf计算,输出文件pwscf.save:
    1. wfc波函数, dw自旋down,up自旋up,若用scf的结果计算pdos有用
    wfcdw[1-nkstot].dat wfcup[1-nkstot].dat
    2. 电荷密度,很有用,仅在scf计算,nscf不计算,nscf,pdos等后处理需要
    charge-density.dat
    3. 计算参数和输出,nscf需要从这里读入原子结构
    data-file-schema.xml
    4. 赝势
    Fe_ONCV_PBE_sr.upf

  • nscf计算时,.in文件中的原子位置会被忽略,但是得有这个选项卡

    Specified atomic positions will be IGNORED and those from the previous scf calculation will be used instead !!!

  • nscf计算,输入文件pwscf.save:
    charge-density.dat 来自scf计算
    data-file-schema.xml 来自nscf或scf计算,读入结构参数等信息,算完覆盖
    修改其中的<output>选项卡,便可调节原子结构等参数进行测试(瞎算)
  • nscf计算,输出文件pwscf.save:
    1. pdos 需要读入,有用
    wfcdw[1-nkstot].dat wfcup[1-nkstot].dat
    2. 计算参数和输出,nscf需要从这里读入原子结构
    data-file-schema.xml
    3. 赝势
    Fe_ONCV_PBE_sr.upf

  • pdos输入需要pwscf.save:
    charge-density.dat 来自scf计算
    wfcdw[1-nkstot].dat wfcup[1-nkstot].dat 来自nscf或scf计算(即scf和nscf的结果都可以用于计算pdos)
    data-file-schema.xml 来自nscf或scf计算,读入结构参数等信息

  • pdos输出pwscf.save:
    atomic_proj.xml

示例

苯分子的scf计算和轨道|psi(ik,iband)|^2展示

工作路径

LearnQE/DOC/material-for-ljubljana-qe-summer-school-master/Day-1/example1.benzene

SCF计算

mpirun -np 4 pw.x -i pw.benzene.scf.in | tee pw.out

PP处理

pp.x -i pp.benzene.psi2.in | tee pp.out

PP仅需要outdir/prefix.save文件夹中的内容,其他SCF计算输出不需要存在,
cp AAA.save BBB.save后,在pp.in中设置prefix = 'BBB'也可以

 &INPUTPP
    !输入数据outdir/prefix.save
    prefix    =  'test',
    !outdir    =  '/tmp',

    !处理保存到filplot_XXX,用于后续处理
    filplot   =  'pp.psi2',

    !plot_num,画图的物理量,更多选项INPUT_PP.html
    plot_num  =  7,
    !不同的plot_num对应不同参数,下面为plot_num = 7,画psi^2时的选项
    !plot_num = 7 输出 filplot_KXXX_BXXX, K k点, B band.
    kpoint    =  1,  !自旋
    kband(1)  =  1,  !画band [kband(1),kband(2)]
    kband(2)  =  16,  
    lsign     =  .true.,  !if true and k point is Gamma, plot |psi|^2 sign(psi)
 /

!可以不指定PLOT namelist,就不会输出额外的图形数据
!从INPUTPP的输出后续处理
&PLOT
   fileout       =  '.xsf', !拓展名
   iflag         =  3, !图形维度
   nfile         =  1, !输入文件数量
   output_format =  5, !输出文件格式,5 => format suitable for XCRYSDEN  (2D or user-supplied 3D region)
   weight(1)     =  1.0, !各输入文件权重
/

使用xcrysden打开xsf格式

xcrysden --xsf pp.psi2_K001_B001.xsf

选定绘图数据:
选择xcrysden的Tools/Data Grid选项,可以查看输入文件包含的数据量,可以分别指定权重,然后点OK
(此处只有一个数据)
下面开始调节画图,设置Isovalue,点击Submit 批量绘图
File/Save Current Statectrl+u,保存操作到state.xcrysden
使用已保存的操作打开新文件

xcrysden --xsf pp.psi2_K001_B002.xsf --script state.xcrysden 

直接保存成图片格式,下面脚本会自动打开xcrysden画图,可以用于批处理

xcrysden --xsf pp.psi2_K001_B002.xsf --script state.xcrysden --print 2.png

笔记本上显示的居然是这样的,应该是对xcrysden窗口的截图,分辨率等原因使得调节窗口覆盖了显示窗口,截图出错 等待后续解决方案
xsf格式也可以用VESTA打开,更美观

石墨烯DOS和能带

工作路径

LearnQE/DOC/material-for-ljubljana-qe-summer-school-master/Day-1/example2.graphene

DOS计算scf(pw.x)->nscf(pw.x)->dos.x

SCF计算

mpirun -np 4 pw.x -i pw.graphene.scf.in | tee scf.out

nscf计算/DOS计算 calculation = 'nscf'
设置occupations='tetrahedra',这同VASP等DFT程序一样,算DOS采用Blochl修正的四面体方法,但是此方法不能用来算力学性质

mpirun -np 4 pw.x -i pw.graphene.nscf.in | tee nscf.out

DOS后处理,可视化dos.x

dos.x -i dos.graphene.in | tee dos.out
./dos.py

同pp.x一样,默认会从otudir/prefix.save中读取数据进行处理pw.x
scf的计算结果也可以用来处理

 &DOS
    prefix = 'pwscf',
    ! otudir = '/tmp'
    !这几个保持默认就好了,算法上的选择
    !bz_sum
    !ngauss
    !degauss
    !能量区间和步长
    !Emin
    !Emax
    !DeltaE
    !保存到fildos
    fildos = 'graphene.dos'
 /

虽然ncsf和scf的结果都可以用来画DOS,但区别还是有的,如图
scf(kpoints 9x9x1 smearing)
nscf(kpoints 12x12x1 tetrahedra)

能带计算scf(pw.x)->bands(pw.x)->bands.x

SCF计算

mpirun -np 4 pw.x -i pw.graphene.scf.in | tee scf.out

能带计算calculation = 'nscf'
k点路径的选择
xcrystal打开pwscf结构,Tools/k-path Selection
左键旋转,右键选中高对称点,OK,保存格式选择pwscf
pw.x计算bands时的k点设置示例

K_POINTS {crystal_b}
3!高对称点数
   0.0000000000     0.0000000000     0.0000000000     40 !此k点到下一k点的点数
   0.6666666667    -0.3333333333     0.0000000000     20
   0.5000000000    -0.5000000000     0.0000000000     0

计算

mpirun -np 4 pw.x -i pw.graphene.bands.in | tee pw.bands.out

能带后处理,可视化bands.x

bands.x -i bands.graphene.in | tee bands.out

输出文件

  • graphene.bands.dat.gnu # 距离 EIG \n …
  • graphene.bands.dat # kx,ky,kz \n EIG[1:nband]

使用py对graphene.bands.dat.gnu画图

./bands.py

bands.graphene.in

 &BANDS
    !读取outdir/pwscf
    prefix  = 'pwscf'
    ! otudir  = '/tmp/'
    !输出文件filband
    filband = 'graphene.bands.dat'
    !spin_component = 1/2 !1 up, 2 down
    !lsigma #noncollinear case
    !
    !If lp = .true. matrix elements of the momentum operator p between conduction and valence bands
    !              are computed and written to file specified in filp
    !lp
    !filp
    
    lsym = .true.,
    !no_overlap 
    ! plot_2d  2D格式 kx,ky,Energy
    !firstk
    !lastk
 /

另外也可以使用plotband.x处理graphene.bands.dat得到的结果和graphene.bands.dat.gnu类似

cndaqiang@girl:example2.graphene$ plotband.x 
     Input file > graphene.bands.dat
Reading    8 bands at     48 k-points
Range:  -18.9560   14.6190eV  Emin, Emax > -18.9560   14.6190
high-symmetry point:  0.0000 0.0000 0.0000   x coordinate   0.0000
high-symmetry point:  0.3333 0.5774 0.0000   x coordinate   0.6667
high-symmetry point:  0.0000 0.5774 0.0000   x coordinate   1.0000
high-symmetry point:  0.0000 0.0000 0.0000   x coordinate   1.5773
output file (gnuplot/xmgr) > plotband.dat.gnu        
bands in gnuplot/xmgr format written to file plotband.dat.gnu   
output file (ps) > 回车
stopping ...

Si状态方程

工作路径

LearnQE/DOC/material-for-ljubljana-qe-summer-school-master/Day-1/example3.Si

状态方程ev.x

使用PW/tools/ev.x可以计算状态方程,直接运行,输入相应参数,数据文件两列[a,E]

cndaqiang@girl:example3.Si/ex3.alat/reference$ ev.x
     Lattice parameter or Volume are in (au, Ang) > au
     Enter type of bravais lattice (fcc, bcc, sc, noncubic) > fcc
     Enter type of equation of state :
     1=birch1, 2=birch2, 3=keane, 4=murnaghan > 1
     Input file > Etot-vs-alat.dat
     Output file > out.dat
     Beware: file out.dat will be overwritten
cndaqiang@girl:example3.Si/ex3.alat/reference$ head out.dat
# equation of state: birch 1st order.  chisq =   0.5913D-07
# a0 = 10.2786 a.u., k0 =  915 kbar, dk0 =  3.96 d2k0 =  0.000 emin =  -15.80790
# a0 =  5.43921 Ang, k0 =  91.5 GPa,  V0 =   271.48 (a.u.)^3,  V0 =   40.23 A^3 

#########################################################################
# Lat.Par       E_calc        E_fit       E_diff    Pressure      Enthalpy
# a.u.            Ry           Ry            Ry        GPa           Ry
#########################################################################
  9.70000     -15.77943     -15.77927    -0.00016      22.46      -15.43108
  9.80000     -15.78874     -15.78890     0.00016      17.39      -15.51061

Al金属体系

工作路径

LearnQE/DOC/material-for-ljubljana-qe-summer-school-master/Day-1/example4.Al

scf计算

金属体系因为费米面的原因,要计算得准确,相比半导体需要增加K点密度,并选择合适的degauss,smearing

    occupations = 'smearing',
    smearing = 'm-v',
    degauss = 0.01

Summer School on Advanced Materials and Molecular Modelling的结果

电荷密度

scf计算

cndaqiang@girl:example4.Al/ex2.chdens$ mpirun -np 4 pw.x -i pw.al.scf.in | tee scf.out

PP计算,xcrysden查看,查看方法与前面的查看xsf格式相同

pp.x -i pp.charge.in 
xcrysden --xsf chargedensity.xsf 

PP输出文件

&INPUTPP
    prefix="pwscf"
    !filplot = "chargedensity"  !输出用于后续处理,不必须输出
    plot_num = 17
    !plot_num = 17 价电子 all-electron valence charge density
    !can be performed for PAW calculations only
    !requires a very dense real-space grid!
    !plot_num = 21 全电子 all-electron charge density (valence+core).
    !spin_component = 0  !0 All, 1 up, 2 down
/
&PLOT
    nfile = 1 
    weight(1)     =  1.0, !各输入文件权重
    iflag  = 3 
    weight(1) = 1.0
    output_format = 5 !6 cube, 5xsf
    !fileout = 'chargedensity.cube' 
    fileout = 'chargedensity.xsf' 
/

铝原子较大可能会覆盖电荷的显示,把原子半径调小就看见了

这次vesta打开chargedensity.xsf没有自动显示电荷密度,不知道是不是数值太小,可以Utilities/2D Data Display查看 把PP的&PLOT调整可绘制cube格式,vesta都支持xsf和cube格式

    output_format = 6
    fileout = 'chargedensity.cube' 

Fe磁性材料

Iron has two remarkable features: it is magnetic and it requires an Ultrasoft PP (USPP) since its localized 3d atomic states are very hard. 若不用USPP,ectoff要取很高
Fe是金属材料,也要设置smearing
工作路径

LearnQE/DOC/material-for-ljubljana-qe-summer-school-master/Day-1/example5.Fe

铁磁FM,反铁磁AFM计算

  • 计算反铁磁,要构建超胞,并指定不同的初始磁矩,还要使用不同的标签(Fe1,Fe2),
    因为初始磁矩和ntyp对应,starting_magnetization(1:ntyp)
    &SYSTEM
    
      ntyp = 2, !设置两种类型
      starting_magnetization(1) =  0.6
      starting_magnetization(2) = -0.6
     !...
    /
    ATOMIC_SPECIES
     Fe1  1.  Fe_ONCV_PBE_sr.upf
     Fe2  1.  Fe_ONCV_PBE_sr.upf
    !...
    

    分别运行

    mpirun -np 4 pw.x -i pw.fe_fm.scf.in | tee fm.out
    mpirun -np 4 pw.x -i pw.fe_afm.scf.in | tee afm.out
    

    检查磁性和能量,可以看到fm总能量/原子数更低

    cndaqiang@girl:example5.Fe$ grep magnetization fm.out | tail -2
       total magnetization       =      1.853765 Bohr mag/cell
       absolute magnetization    =      1.863940 Bohr mag/cell
    cndaqiang@girl:example5.Fe$ grep ! fm.out 
    !    total energy              =    -240.27909971 Ry
    cndaqiang@girl:example5.Fe$ grep magnetization afm.out | tail -2
       total magnetization       =     -0.000295 Bohr mag/cell
       absolute magnetization    =      0.008622 Bohr mag/cell
    cndaqiang@girl:example5.Fe$ grep ! afm.out 
    !    total energy              =    -480.47700464 Ry
    

DOS

nscf算DOS时,占据方式为occupations = 'tetrahedra',不用设置smearing等参数

mpirun -np 4 pw.x -i pw.fe_fm.scf.in | tee fm.out
mpirun -np 4 pw.x -i pw.fe_fm.nscf.in | tee fmnscf.out
dos.x -i dos.in 
./dos_nspin.py #画图

可以对比同样的k点(12x12x12)使用scf和nscf的结果分别画的dos,nscf的结果更平滑

PDOS projwfc.x

mpirun -np 4 pw.x -i pw.fe_fm.scf.in | tee fm.out
mpirun -np 4 pw.x -i pw.fe_fm.nscf.in | tee fmnscf.out
projwfc.x -i pdos.in 
./pdos_nspin.py  #画图

生成类似fe.pdos.pdos_atm#1(Fe)_wfc#1(s)文件名的数据文件
文件结构(nspin=2时加倍):

  • s
    ldos(S) pdos(S)
  • p
    ldos(p) pdos: px py pz
  • d
    ldos(d) pdos:dz2 dxz dyz dxy dx2-y2

LDOS: s, p, d, f轨道的态密度
PDOS: s, p, d, f轨道的投影态密度
projwfc.x的输入文件和dos.x输入类似

&PROJWFC
    filpdos = 'fe.pdos' !pdos是filpdos,dos是fildos
    Emin = 8.0
    Emax = 25.0,
    DeltaE = 0.1
/

FeO DFT+U

QE仅支持特定元素+U
/Modules/set_hubbard_n.f90我们可以看到支持的元素和加U轨道的主量子数n
/Modules/set_hubbard_l.f90我们可以看到支持的元素和加U轨道的轨道量子数l
对于Fe,n=3,l=2,即3d轨道加U
如果所要加U的元素不再QE支持范围内,可以改代码支持

DFT+U (formerly known as LDA+U) currently works only for a few selected elements. Modify Modules/set_hubbard_l.f90 and PW/src/tabd.f90 if you plan to use DFT+U with an element that is not configured there.

赝势的选取很影响结果,如左边Fe_ONCV_PBE_fr.upf算得FeO为导体,右边Fe.pbesol-spn-kjpaw_psl.0.2.1.UPF算的FeO为半导体
而使用Fe.pbesol-spn-kjpaw_psl.1.0.0.UPF,直接报错

     axis vectors are left-handed
               file Fe.pbesol-spn-kjpaw_psl.1.0.0.UPF: wavefunction(s)  3P 3D renormalized
               file Fe.pbesol-spn-kjpaw_psl.1.0.0.UPF: wavefunction(s)  3P 3D renormalized
           3          10

 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
     Error in routine set_dft_from_name (1):
      conflicting values for igcx
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

     stopping ...

因此之后的计算采用推荐的赝势

Fe1  55.845  Fe.pbesol-spn-kjpaw_psl.0.2.1.UPF
Fe2  55.845  Fe.pbesol-spn-kjpaw_psl.0.2.1.UPF
O     16.0   O.pbesol-n-kjpaw_psl.0.1.UPF

工作路径

LearnQE/DOC/material-for-ljubljana-qe-summer-school-master/Day-2/example4.functionals/ex1.DFT+U

不加U

虽然FeO实验是半导体,但是计算是导体,所以SCF计算要设置smearing
反铁磁材料,所以建超胞,Fe标签不同,用于指定不同的初始磁矩

mpirun -np 4 pw.x -i pw.FeO.scf.in | tee scf.out
mpirun -np 4 pw.x -i pw.FeO.nscf.in | tee nscf.out
projwfc.x -i projwfc.FeO.in
./pdos_nspin.py

加U

添加U值在scf和nscf的计算中,参数相同

    lda_plus_u = .true.,
    !+U的类型
    lda_plus_u_kind = 0, 
    !   0   simplified version of Cococcioni and de Gironcoli,
    !   PRB 71, 035105 (2005), using Hubbard_U
    !   1   rotationally invariant scheme of Liechtenstein et al.,
    !   using Hubbard_U and Hubbard_J

    !U_projection_type轨道的选择方式
    U_projection_type = 'atomic',
    !Hubbard_U(i), i=1,ntyp 仅对需要的元素+U值, 
    !U值可根据Timrov, N. Marzari, M. Cococcioni, PRB 98, 085127 (2018)计算
    Hubbard_U(1) = 5.2
    Hubbard_U(2) = 5.2

执行

mpirun -np 4 pw.x -i pw.FeO.scfU.in | tee scfU.out
mpirun -np 4 pw.x -i pw.FeO.nscfU.in | tee nscfU.out
projwfc.x -i projwfc.FeO.in
./pdos_nspin.py

加U后算出来的仍是导体,根据Summer School on Advanced Materials and Molecular Modelling,优化到了局域最小值
通过设置starting_ns_eigenvalue(m,ispin,ityp),计算出全局最小值

加U & starting_ns_eigenvalue(m,ispin,ityp)

在scf时添加,尽在第一个scf初始迭代时生效

    starting_ns_eigenvalue(5,2,1) = 1.0
    starting_ns_eigenvalue(5,1,2) = 1.0

选择5,2,1与5,1,2的原因是,根据scf+U的屏幕输出,在第一个scf内,atom1的spin2的第5个本征值(0.296)是非简并的(前面四个都是2重简并)
atom2同理

Why 5-th state? Because it is the one which is non-degenerate and if occupied fully could lead to an insulating result.

     iteration #  1     ecut=    30.00 Ry     beta= 0.30
......
atom    1   Tr[ns(na)] (up, down, total) =   5.00636  1.09500  6.10136
   spin  1
    eigenvalues: 
  1.000  1.000  1.002  1.002  1.002
......
   spin  2
    eigenvalues: 
  0.129  0.129  0.270  0.270  0.296
......

在输出会有这样的语句提示 Modify starting ns matrices according to input values ,然后修改values 其余不变

mpirun -np 4 pw.x -i pw.FeO.scfUns.in | tee scfUns.out
mpirun -np 4 pw.x -i pw.FeO.nscfU.in | tee nscfU.out
projwfc.x -i projwfc.FeO.in
./pdos_nspin.py

三种情况的DOS,+U与starting_ns_eigenvalue可以算得半导体的情况

通过比较能量,调整之后能量是更低了一些 -735.172(+U & starting_ns_eigenvalue ) < -734.99424361(+U)

光学吸收谱计算

参考QE实践详解

手册/PP/Doc/eps_man.pdf

输入参数示例,下面的内容是epsilon.x的输入

&inputpp
  !outdir='./tmp'
  calculation='eps'
/
&energy_grid
  smeartype='gauss'
  intersmear=0.50 !这里0.5已经很大了,太大会出现峰特别宽,就像fft的dumpling取大了一样,用0.136好一些
  intrasmear=0.0
  wmin=0.0  !能量区间,单位eV
  wmax=60.0 !能量区间,单位eV
  nbndmin=1
  nbndmax=0
  nw=2000
  shift=0.0
/

calculation='eps'时输出epsr.dat, epsi.dat, eels.dat and ieps.dat

  • epsr.dat, epsi.dat介电函数张量的实部虚部$\epsilon_{1,\alpha,\alpha}(\omega),\epsilon_{2,\alpha,\alpha}(\omega)$
    虚部对应吸收谱
  • eels.dat电子能量损失谱
  • ieps.dat the diagonal components of dielectric tensor calculated on the imaginary axis of frequency (via London transformation) $\epsilon_{\alpha,\alpha}(i\omega)$

运行

epsilon.x -i epsilon.in | tee epsilon.out
epsi_pwscf2gpaw.py epsi_pwscf.dat
spectrum.gpaw.py ./gpaw.epsi_pwscf.dat
vs ./gpaw.epsi_pwscf.dat.png

所以nbnd要取多一些,不然能量的区间太小了

报错

     non uniform kpt grid

则要关闭对称性,使所有k点权重相同

noinv=T
nosym=T

本文首发于我的博客@cndaqiang.
本博客所有文章除特别声明外,均采用 CC BY-SA 4.0 协议 ,转载请注明出处!



评论


广告

目录

广告
访客数据