if [ ! -d "$DIR" ];then echo Please specify the source directory. exit fi
forNAMEin 3RD.POSCAR.* do NUM="${NAME#3RD.POSCAR.}" # the number of the 3nd POSCAR file CODE=$(head -1 $NAME) # the code of the POSCAR file # find the source POSCAR file in the source directory INFO=$(grep $CODE $DIR/3RD.POSCAR.* | head -1) [ x"$INFO" == x ] && echo "NOT FOUND for job-$NUM" && continue SOURCE_NUM="${INFO%:*}" # the number of the source POSCAR file SOURCE_NUM="${SOURCE_NUM##*.}"
#module load spglib/master #thirdorder_vasp.py sow 331 -10 #生成2*2*2的超胞 0.6表示截断半径为6埃。 #sow a b c m abc表示生成的超胞大小，m为正数的时候表示截断半径（单位：nm）（例如0.6，表示截断半径为6.0埃）；为负数时表示第几近邻截断（此时应取负整数，例如-3表示第三近邻）
for i in 3RD.POSCAR.* do s=$(echo $i| cut -d "." -f 3) d=job-$s mkdir $d cp $i $d/POSCAR cp ./INCAR ./sub.sh ./KPOINTS ./POTCAR $d pushd $d popd done
批量提交任务
{.line-numbers}
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17
#!/bin/sh root_path='pwd'
#在每个job文件夹里进行静态计算。 for i in {001..1050} #按顺序从job-001到1050做静态vasp计算,根据生成的文件夹设定 do cd job-$i echo "job-$i", time sbatch sub.sh wait cd .. done
Besides the routine output files from previous ShengBTE program, FourPhonon generates these output files:
BTE.Numprocess_4ph: number of allowed four-phonon scattering processes, for each irreducible q point and phonon band
BTE.P4: phase space available for four-phonon processes, for each irreducible q point and phonon band
BTE.P4_total: total volume in phase space available for four-phonon processes
BTE.P4_plusplus*, BTE.P4_plusminus*, BTE.P4_minusminus*: similar to BTE.P4 but only includes contributions from ++/+-/- - processes
*Note for four-phonon scatterings, there are three different channels: recombination (++), redistribution (+-) and splitting (- -) processes.
Under temperature-dependent directories: (the unit of scattering rates is ps−1 if not specified)
BTE.WP4: weighted phase space available for four-phonon processes ([rad/ps]$^{-5}$, 2nd column) vs angular frequency (rad/ps, 1st column) for each irreducible q point and phonon band
BTE.WP4_plusplus*, BTE.WP4_plusminus*, BTE.WP4_minusminus*: similar to BTE.WP4 but only includes contributions from ++/+-/- - processes
BTE.w_3ph: three-phonon scattering rates for each irreducible q point and phonon band, this file replaces the original output BTE.w_anharmonic. Absorption and emission processes are written out into BTE.w_3ph_plus and BTE.w_3ph_minus
BTE.w_4ph: four-phonon scattering rates for each irreducible q point and phonon band. Similarly, we provide the contributions from different channels: BTE.w_4ph_plusplus, BTE.w_4ph_plusminus and BTE.w_4ph_minusminus
BTE.w_4ph_normal: four-phonon scattering rates of normal processes, for each irreducible q point and phonon band. This file has 5 columns and each one represents: angular frequency in rad/ps, recombination channel, redistribution channel, splitting channel, overall scattering rates from normal processes
BTE.w_4ph_Umklapp: four-phonon scattering rates of Umklapp processes, for each irreducible q point and phonon band. The format is the same as BTE.w_4ph_normal
BTE.kappa*: thermal conductivity results are written out as usual
1987年，Baroni、Giannozzi和Testa提出了DFPT方法。DFPT通过计算系统能量对外场微扰的响应来求出晶格动力学性质。该方法最大的优势在于它不限定微扰的波矢与原胞边界正交，不需要超原胞也可以对任意波矢求解。Castep、Vasp等采用的是一种linear response theory 的方法（或者称为 density perturbation functional theory，DFPT），直接计算出原子的移动而导致的势场变化，再进一步构造出动力学矩阵。进而计算出声子谱。
# -*- coding: utf-8 -*- """ Qiang Li Command: python3 vib_correct.py This is a temporary script file. """ import numpy as np from ase.ioimport read, write import os
# 获取虚频对应的OUTCAR部分 vib_lines = lines[l_position:l_position + num_atoms] #振动部分 7222到7248行 #print(vib_lines) vib_dis = [] for line invib_lines: #model_positions = [float(i) for i in line.rstrip().split()[:3]] vib_infor = [float(i) for i in line.rstrip().split()[3:]] # dx， dy， dz对应的那三列 vib_dis.append(vib_infor) vib_dis = np.array(vib_dis) #将振动部分写到一个array中。
#!/usr/bin/env python # -*- coding: utf-8 -*- """ Qiang Li Command: python3 vib_correct.py This is a temporary script file. """ import numpy as np import os
# 获取虚频对应的OUTCAR部分 num_atoms = sum([int(i) for i in lines_pos[6].rstrip().split()]) vib_lines = lines[l_position:l_position + num_atoms] #振动部分 7222到7248行
model_positions = [] vib_dis = [] for line invib_lines: position = [float(i) for i in line.rstrip().split()[:3]] vib_infor = [float(i) for i in line.rstrip().split()[3:]] # dx， dy， dz对应的那三列 model_positions.append(position) vib_dis.append(vib_infor)
###保存结构 f_out = open('POSCAR_new','w') f_out.writelines(lines_pos[:8]) f_out.write('Cartesian\n') for i innew_positions: f_out.write(' '.join([str(coord) for coord in i]) + ' F F F\n') f_out.close()
随后进行数据提取，以vaspkit1.12版本为例，打开vaspkit，输入73，打开VASP2other Interface 子功能，选择739 Sort Phonon Band Structure for Phononpy，运行完成即可获得group velocity相关数据以及各原子的声子谱数据；在运行过程中我们可以看到软件说明find group velocity并写入输出文件中。 1.5版本为789功能。
generate 1000 supercells with displacements at 300K.
New in v2.25 Symfc (symfc/symfc) is based on fitting approach and any displacements set of atoms in supercell can be handled. For example, random displacements generated by RANDOM_DISPLACEMENTS can be used to compute force constants. To use symfc, its python module has to be installed via conda-forge or pip.
New in v2.3 ALM (ttadano/ALM) is based on fitting approach and any displacements set of atoms in supercell can be handled. For example, random displacements generated by RANDOM_DISPLACEMENTS can be used to compute force constants. To use ALM, its python module has to be installed via conda-forge or building it. The installation instruction is found here.
When ALM is used, please cite the paper: T. Tadano and S. Tsuneyuki, J. Phys. Soc. Jpn. 87, 041015 (2018).
Dielectric constants, piezoelectric constants and polar-phonon frequency：
DFPT方法，Note, DFPT cannot be used with hybrid exchange-correlation functionals. In these cases the LCALCEPS flag should be used in combination with IBRION = 6
amset run –directory path/to/files 运行即可 其他设置可以通过命令行加入： amset run –interpolation-factor 20 画图： amset plot rates mesh_99x99x99.h5
Output filesThe transport results will be written to transport_99x99x99.json. The 99x99x99 is thefinal interpolated mesh upon which scattering rates and the transport properties arecomputed.The write_mesh option is set to True, meaning scattering rates for all doping levelsand temperatures will be written to the mesh_99x99x99.h5 file.
zero_weighted_kpoints： 设置方式： keep: Keep zero-weighted k-points in the band structure. drop: Drop zero-weighted k-points, keeping only the weighted k-points. prefer: Drop weighted-kpoints if zero-weighted k-points are present in the calculation (useful for cheap hybrid calculations).
Use amset plot band --stats vasprun.xml: This calculates the effective mass in the same way as sumo but using a band structure interpolated from a uniform calculation.
orbitalWise 这个参数一般加在cohpBetween或者cohpGenerator之后，用以把原子之间的相互作用具体分到s p d 这些轨道上，比如计算s-s, p-p, d-d之间的pCOOP, pCOHP； 比如：cohpBetween atom 1 atom 2 orbitalWise
cohpGenerator from 1.4 to 1.5 type Ga type Sr orbitalWise
在 /vaspkit.0.72/examples/LDOS_PDOS/Partial_DOS_of_CO_on_Ni_111_surface 目录下启动 vaspkit ，输入命令 115 选择子功能 The Sum of Projected DOS for Selected Atoms and orbitals 。我们需要提取的是O的s，p轨道，C的s，p轨道以及表面的dxy和d3z2-r2(dz2)轨道。首先提示你选择元素（累加），第一次输入元素O，回车后提示你输入提取的轨道名（累加），第一次输入轨道s。回车后重复以上两个操作。如果想结束输入，在元素选择行直接按回车键结束输入。
元素行接受自由格式输入，你可以输入以下内容元素行接受自由格式输入，你可以输入以下内容1-3 4 Ni表示选择元素1,2,3,4和Ni元素的PDOS进行累加。轨道行只支持标准输入，如果使用了LORBIT=10不投影轨道，你只能选择 s p d f中的一个或多个，如果使用了LORBIT=11投影了轨道，你可以从s p d f和s py pz px dxy dyz dz2 dxz dx2 f-3 f-2 f-1 f0 f1 f2 f3中选择一个或多个输入。轨道行还支持输入all计算所有轨道的DOS之和。
比如元素行输入C O ,对应的轨道选择s px，那么提取的就是所有C和O元素的s轨道和px轨道的PDOS之和。 PDOS_USER.dat中的轨道名 为#Energy C&O_s&px，显而易懂。
以 CO 分子吸附在 Ni(100) 表面为例： 步骤一：优化 CO/Ni(100) 的结构； 步骤二：分别计算 Ni(100) 和 CO 的单点能 CO 和 Ni(100) 片段的坐标从 CO/Ni(100) 的 CONTCAR 里直接截取，不要再结构优化！！计算时也要保证三次自洽计算所采用的 FFT mesh 一致（NGXF,NGYF,NGZF） 步骤三：用 vaspkit 314 功能做电荷差分。依次输入三个片段的 CHGCAR 路径:
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16
314 ======== File Options ======== Input the Names of Charge/Potential Files with Space: (e.g., to get AB-A-B, type: ~/AB/CHGCAR ./A/CHGCAR ../B/CHGCAR) (e.g., to get A-B, type: ~/A/CHGCAR ./B/CHGCAR)
------------>> ./CHGCAR ./co/CHGCAR ./slab/CHGCAR
-->> (01) Reading Structural Parameters from ./CHGCAR File... -->> (02) Reading Charge Density From ./CHGCAR File... -->> (03) Reading Structural Parameters from ./co/CHGCAR File... -->> (04) Reading Charge Density From ./co/CHGCAR File... -->> (05) Reading Structural Parameters from ./slab/CHGCAR File... -->> (06) Reading Charge Density From ./slab/CHGCAR File... -->> (07) Written CHGDIFF.vasp File!
生成 CHGDIFF.vasp 包含电荷密度差的信息，可以直接导入到 VESTA 里作图。 图片
Vaspkit计算变形电荷密度差
变形电荷密度差：自洽计算收敛以后体系的电荷密度减去该原子构型下每个原子的球对称的电荷密度(即初猜电荷密度)。 以 CO 分子为例： 步骤一：先自洽计算优化 CO 分子。 步骤二：新键文件夹，在优化好的结构基础上用 ICHARG = 12 做非自洽计算。 步骤三：用 vaspkit 314 功能做电荷差分。依次输入两个 CHGCAR 路径。
1 2 3 4 5 6 7 8 9 10 11 12 13 14
314 ======================= File Options ============================ Input the Names of Charge/Potential Files with Space: (e.g., to get AB-A-B, type: ~/AB/CHGCAR ./A/CHGCAR ../B/CHGCAR) (e.g., to get A-B, type: ~/A/CHGCAR ./B/CHGCAR)
------------>> ./CHGCAR ./deform/CHGCAR
-->> (01) Reading Structural Parameters from ./CHGCAR File... -->> (02) Reading Charge Density From ./CHGCAR File... -->> (03) Reading Structural Parameters from ./deform/CHGCAR File... -->> (04) Reading Charge Density From ./deform/CHGCAR File... -->> (05) Written CHGDIFF.vasp File!
314 ======================= File Options ============================ Input the Names of Charge/Potential Files with Space: (e.g., to get AB-A-B, type: ~/AB/CHGCAR ./A/CHGCAR ../B/CHGCAR) (e.g., to get A-B, type: ~/A/CHGCAR ./B/CHGCAR)
------------>> ./CHGCAR ./noele/CHGCAR
-->> (01) Reading Structural Parameters from ./CHGCAR File... -->> (02) Reading Charge Density From ./CHGCAR File... -->> (03) Reading Structural Parameters from ./noele/CHGCAR File... -->> (04) Reading Charge Density From ./noele/CHGCAR File... -->> (05) Written CHGDIFF.vasp File! 步骤四：得到 CHGDIFF.vasp 文件，导入到 VESTA 里做图。
ALGO = Damped HFSCREEN = 0.2 LHFCALC = .True. LMAXFOCK = 4 PRECFOCK = Fast TIME = 0.4
+SOC计算电子结构：
结构优化（PBE）
静态自洽scf计算（PBE+SOC）KPOINTS和POTCAR文件不变，写波函数和电荷密度
非自洽计算（PBE+SOC） Spin-Orbit Coupling Calculation
{.line-numbers}
1 2 3 4 5 6 7 8 9
ISPIN = 2 NELMIN = 6 (Min electronic SCF steps) LSORBIT = .TRUE. (ActivateSOC) GGA_COMPAT = .FALSE. (Apply spherical cutoff on gradient field) #VOSKOWN = 1 (Enhances the magnetic moments and the magnetic energies) LMAXMIX = 2 (For d elements increase LMAXMIX to 4, f: LMAXMIX = 6) SAXIS = 001 (Directionof the magnetic field) #MAGMOM = 003 (Setthis parameters manually, Local magnetic moment parallel to SAXIS, 3*NIONS*1.0for non-collinear magnetic systems) NBANDS = 32 (Setthis parameters manually, 2 * number of bands of collinear-run)
调节SOC强度
调节自旋轨道耦合作用的大小
要调优VASP中的自旋轨道耦合强度，需要修改源文件vasp_source_code_path/src/relativistic.F 并重新编译。如果要将强度减少一半，可以在relativistic.F的129行L●S项乘以0.5d0，如下： DO I=0,1DO J=0,1DO M =1,2LL+1DO MP=1,2LL+1 DLLMM(LMP+MP-1,LM+M-1,J+2I+1)=DLLMM(LMP+MP-1,LM+M-1,J+2I+1)+ & 0.5d0SUMLS(M,MP,I+2*J+1,LL) !!!line 129 relativistic.F fileEND DOEND DOEND DOEND DO
# Frequency dependent dielectric tensor including # local field effects within the RPA (default) or # including changes in the DFT xc-potential (LRPA=.FALSE.). # N.B.: beware one first has to have done a # calculation with ALGO=Exact, LOPTICS=.TRUE. # and a reasonable number of virtual states (see above) ALGO = GW0 ; LSPECTRAL = .TRUE. ; NOMEGA = 50 # be sure to take the same number of bands as for # the LOPTICS=.TRUE. calculation, otherwise the # WAVEDER file is not read correctly NBANDS = 64 # Add this to update the quasiparticle energies # in the Green's function (GW0) #NELM = 4 ISMEAR = 0 SIGMA = 0.05 EDIFF = 1E-8
k-point 1 : 0.0000 0.0000 0.0000 band No. KS-energies QP-energies sigma(KS) V_xc(KS) V^pw_x(r,r') Z occupation .........
LWANNIER90计算能带结构
重新开始第三步
INCAR：
1 2 3 4 5 6 7 8
ALGO = GW0 ; LSPECTRAL = .TRUE. ; NOMEGA = 50 #LRPA = .FALSE. ## be sure to take the same number of bands as for ## the LOPTICS=.TRUE. calculation, otherwise the ## WAVEDER file is not read correctly NBANDS = 64 ##VASP2WANNIER90 LWANNIER90=.TRUE.
This run generates the wannier90 standard output (wannier90.wout) and the file wannier90.chk needed for the wannier interpolation (next step)
Obtain bandstructure (Wannier interpolation) Uncomment the bandstructure plot flags in wannier90.win and rerun (restart) wannier90:
wannier90.x wannier90
This run generates the following bandstructure files which can be visualized using xmgrace or gnuplot:
wannier90_band.agr
wannier90_band.dat
wannier90_band.gnu
to plot the band structure using gnuplot: gnuplot -persist wannier90_band.gnu
Beyond the random-phase-approximation
To include local field effects beyond the random-phase-approximation in the description of the frequency dependent dielectric response function (local field effects in DFT) add the following line to your INCAR file:
LRPA = .FALSE. and again restart from the WAVECAR and WAVEDER files from step 2.
Beyond G0W0: GW0
The most usual step beyond single-shot GW (G0W0) is to iterate the quasi-particle energies in the Greens functions. This is the socalled GW0 approximation. To have VASP do, for instance, 4 iterations of the QP-energies in G, add the following line to the INCAR file:
NELM = 4 and again restart from the WAVECAR and WAVEDER files from step 2.
To quickly find the QP-energy of the highest lying occupied state after 4 iterations of the QP energies in G, type:
gnuplot <<EOF set term dumb set title 'Energy of each ionic step' set xlabel 'Ionic steps' set ylabel 'Energy(eV)' plot 'temp.e' u 1:5 w l t "Energy in eV" set title 'Max Force of each ionic step' set xlabel 'Ionic steps' set ylabel 'Force (eV/Angstrom)' plot 'force.conv' w l t "Force in eV/Angstrom" EOF
if [ $steps -gt 8 ] ; then tail -5 force.conv >temp.fff tail -5 temp.e >temp.ee gnuplot <<EOF set term dumb set title 'Energy of each ionic step for the last few steps' set xlabel 'Ionic steps' set ylabel 'Energy(eV)' plot 'temp.ee' u 1:5 w l t "Energy in eV" set title 'Max Force of each ionic step for the last few steps' set xlabel 'Ionic steps' set ylabel 'Force (eV/Angstrom)' plot 'temp.fff' w l t "Force in eV/Angstrom" EOF rm temp.fff temp.ee fi
rm temp.e temp.f temp.ff temp.fix temp.fixx
批量提交任务
适用于测试、有限位移法求力常数矩阵等情况。 简单： for i in *; do cd $i ; qsub sub4; cd $OLDPWD; done
{.line-numbers}
1 2 3 4 5 6 7 8 9 10 11 12
#!/bin/sh root_path='pwd'
#在每个job文件夹里进行静态计算。 for i in {001..1316} #按顺序从job-001到116做静态vasp计算,根据生成的文件夹设定 do cd job-$i echo "job-$i", time sbatch sub.sh wait cd .. done
批量产生文件夹
i in {2..9}; do cp 0.01 0.0$i -r ; done ``` #批量复制0.01文件夹复制成0.02-0.09
# -*- coding: utf-8 -*- """ Qiang Li Command: python3 vib_correct.py This is a temporary script file. """ import numpy as np from ase.ioimport read, write import os
# 获取虚频对应的OUTCAR部分 vib_lines = lines[l_position:l_position + num_atoms] #振动部分 7222到7248行 #print(vib_lines) vib_dis = [] for line invib_lines: #model_positions = [float(i) for i in line.rstrip().split()[:3]] vib_infor = [float(i) for i in line.rstrip().split()[3:]] # dx， dy， dz对应的那三列 vib_dis.append(vib_infor) vib_dis = np.array(vib_dis) #将振动部分写到一个array中。
#!/usr/bin/env python # -*- coding: utf-8 -*- """ Qiang Li Command: python3 vib_correct.py This is a temporary script file. """ import numpy as np import os
# 获取虚频对应的OUTCAR部分 num_atoms = sum([int(i) for i in lines_pos[6].rstrip().split()]) vib_lines = lines[l_position:l_position + num_atoms] #振动部分 7222到7248行
model_positions = [] vib_dis = [] for line invib_lines: position = [float(i) for i in line.rstrip().split()[:3]] vib_infor = [float(i) for i in line.rstrip().split()[3:]] # dx， dy， dz对应的那三列 model_positions.append(position) vib_dis.append(vib_infor)
###保存结构 f_out = open('POSCAR_new','w') f_out.writelines(lines_pos[:8]) f_out.write('Cartesian\n') for i innew_positions: f_out.write(' '.join([str(coord) for coord in i]) + ' F F F\n') f_out.close()
nebresult.pl没用过，看结果“forces and energy:”下面的数据的第四列应该就是切线力，正常的计算切线力是由正值变到负值，或者由负值变到正值，然后找接近于0的点。你图里的数据只有正值，因此要么初末态构型有问题，要么你要找的点在0和1两个image之间。 不过1号image的切线力已经很接近0了，可以尝试直接做过dimer
另一个可以用来观察收敛情况的命令是nebbarrier.pl（同样不带参数）。该命令没有输出到屏幕的内容，而是生成neb.dat文件。 neb.dat文件第二列表示距离（即临近两结构的dist.pl的计算结果），第三列表示能量（以初态能量为参考值），第四列为力（forces along the neb）。
EDIFFG参数对应的力是nebef.pl输出中的force of images in the neb.
$ unzip disorder-main $ cd disorder-main $ make disorder (Only disorder is compilied) $ make supercell (Only supercell is compilied) $ make or make all (Both disorder and supercell are compilied) $ make clean (rm *.o *.mod disorder supercell) export PATH=$PATH:disorder_installation_path/disorder-main/bin (Add this to the ~/.bashrc file as you wish)
Note: The default compiler is ifort, and you can modify the Makefile in the src directory to use another compiler (e.g. gfortran).
需要准备两个文件SPOSCAR和INDSOD,以examples的1和2为例：
例子1：SnxPb1-xTe
扩建超胞，保证整除
1 2 3 4 5 6 7 8 9
Te Pb 1.0 13.1323995589999996 0.0000000000000000 0.0000000000000000 0.0000000000000000 13.1323995589999996 0.0000000000000000 0.0000000000000000 0.0000000000000000 13.1323995589999996 Te Pb 32 32 Direct * * *
准备INDSOD文件
1 2 3 4 5 6 7 8 9 10 11 12
&input nsub = 2 subs = 16,16 symb = 'Sn','Pb' ! The quotes is unnecessary for the ifort compiler site = 2 prec = 1D-5 ! lpro = .true. ! lpos = .true. ! leqa = .true. ! lspg = .true. lcfg = .false. /
其全称为 Alloy Theoretic Automated Toolkit (ATAT) 引用官网的一句话，ATAT is a generic name that refers to a collection of alloy theory tools developped by Axel van de Walle, in collaboration with various research groups and with various sources of financial support.
maps.log：开始可能的warnings: ’Not enough known energies to fit CE’ ’True ground states not = fitted ground states’ ’New ground states predicted, see predstr.out’ 随着结构增多: ’Among structures of known energy, true and predicted ground states agree.’ ’No other ground states of xx atoms/unit cell or less exist.’ 最后一行有crossvalidation score (in eV/atom)，该值要小于0.025。精确小于0.001.
fit.out:拟合的结构，一个结构一行，共6列，分别代表concentration、 energy 、fitted_energy (energy-fitted_energy)、 weight 、index
predstr.out：预测的能量，一个结构一行，共5列，分别代表concentration、 energy 、predicted_energy、index、status.其中status包括status is either b for busy (being calculated), e for error, u for unknown (not yet calculated) or g if that structure is predicted to be a ground state (g can be combined with the above).
gs.out：基态能量，一个结构一行，四列，分别是concentration、 energy、 fitted_energy、 index
# # To run VASP this script calls $vasp_std # (or posibly $vasp_gam and/or $vasp_ncl). # These variables can be defined by sourcing vaspcmd . vaspcmd 2> /dev/null
# # When vaspcmd is not available and $vasp_std, # $vasp_gam, and/or $vasp_ncl are not set as environment # variables, you can specify them here [ -z "`echo $vasp_std`" ] && vasp_std="mpirun -np 8 /path-to-your-vasp/vasp_std" [ -z "`echo $vasp_gam`" ] && vasp_gam="mpirun -np 8 /path-to-your-vasp/vasp_gam" [ -z "`echo $vasp_ncl`" ] && vasp_ncl="mpirun -np 8 /path-to-your-vasp/vasp_ncl"
The optimal length of the lattice vector c normal to the stacking direction is determined in a series of single point calculations with varied value of c (all other degrees of freedom are fixed at their experimental values).
The computed c vs. energy dependence is written in the file results.dat and can be visualized e.g. using xmgrace. The optimal value can be obtained using the attached utility (python with numpy or Numeric is needed):
./utilities/fit.py results.dat This should yield:
1 2 3 4 5 6 7
200 iterations performed Ch-square: 4.30305519481e-09 ---------
E0(eV): -37.433456779 d0(A): 6.65603352689 The computed value of 6.66 Å agrees well with experiment (6.71 Å).
# # To run VASP this script calls $vasp_std # (or posibly $vasp_gam and/or $vasp_ncl). # These variables can be defined by sourcing vaspcmd . vaspcmd 2> /dev/null
# # When vaspcmd is not available and $vasp_std, # $vasp_gam, and/or $vasp_ncl are not set as environment # variables, you can specify them here [ -z "`echo $vasp_std`" ] && vasp_std="mpirun -np 8 /path-to-your-vasp/vasp_std" [ -z "`echo $vasp_gam`" ] && vasp_gam="mpirun -np 8 /path-to-your-vasp/vasp_gam" [ -z "`echo $vasp_ncl`" ] && vasp_ncl="mpirun -np 8 /path-to-your-vasp/vasp_ncl"
# # The real work starts here #
# Here the work starts rm results.dat
drct=$(pwd)
for i in graphene graphite do cd $drct/$i $vasp_std done
cd $drct
# obtain total energy for graphite en2=$(grep "free ene" graphite/OUTCAR |tail -1|awk '{print $5}')
# obtain total energy for graphene en1=$(grep "free ene" graphene/OUTCAR |tail -1|awk '{print $5}')
echo "Binding energy (eV/atom): " $deltaE >results.dat
Note that the calculation is performed in two steps (two separate single-point calculations) in which the energy for bulk graphite and for graphene are obtained. The binding energy is computed automatically and it is written in the file results.dat. (N.B.: for the latter python needs to be available.)
The computed value of 0.050 eV/A is now fairly close to the RPA reference of 0.048 eV/atom [1].
# # To run VASP this script calls $vasp_std # (or posibly $vasp_gam and/or $vasp_ncl). # These variables can be defined by sourcing vaspcmd . vaspcmd 2> /dev/null
# # When vaspcmd is not available and $vasp_std, # $vasp_gam, and/or $vasp_ncl are not set as environment # variables, you can specify them here [ -z "`echo $vasp_std`" ] && vasp_std="mpirun -np 8 /path-to-your-vasp/vasp_std" [ -z "`echo $vasp_gam`" ] && vasp_gam="mpirun -np 8 /path-to-your-vasp/vasp_gam" [ -z "`echo $vasp_ncl`" ] && vasp_ncl="mpirun -np 8 /path-to-your-vasp/vasp_ncl"
# # The real work starts here #
rm results.dat
drct=$(pwd)
for i in graphene graphite do cd $drct/$i $vasp_std done
cd $drct
# obtain total energy for graphite en2=$(grep "free ene" graphite/OUTCAR |tail -1|awk '{print $5}')
# obtain total energy for graphene en1=$(grep "free ene" graphene/OUTCAR |tail -1|awk '{print $5}')
echo "Binding energy (eV/atom): " $deltaE >results.dat
Note that the calculation is performed in two steps (two separate single-point calculations) in which the energy for bulk graphite and for graphene are obtained. The binding energy is computed automatically and it is written in the file results.dat. (N.B.: for the latter python needs to be available.)
Even though the TS method predicts a reasonable geometry (see the Graphite interlayer distance example) it overestimates the energetics strongly: the computed binding energy of -0.083 eV/atom is too large compared to the RPA reference of 0.048 eV/atom. This overestimation is - at least in part - due to neglecting the many-body interactions (see example Graphite MBD binding energy).
自洽计算其实也是单点计算，静态自洽计算前需要提供 较稳定体系 的晶格结构信息（即结构优化完的CONTCAR），从而通过电子自洽计算，通过自洽迭代求解薛定谔方程（微观中描述电子的状态，相当于宏观的运动方程）） 完整地计算出体系基态下费米能级（准确值，The fermi level location is accurate only in self-consistent calculation.）、 电子的波函数、电荷密度等信息，可以直接分析原子间的键合作用，也可以在非自洽之后进一步分析晶体的电子结构和材料的相关性质。
xy 方向上表面的大小（15 Å或者100原子左右）; A 这个影响表面吸附物种的覆盖度; B 影响体系的尺寸大小和计算时间; C 不同的大小需要选取对应的K点；回想一下我们前面提到的经验规则。
不同的晶面，(111)， (100)，(110); A 这取决于你研究的方向;低指数面、解离面。 B 不同晶面的表面能; C 不同晶面的表面结构，反应活性等。
表面结构的层数（3+3以上，6Å+6Å） A 层数多了，原子变多，体系在 z 方向的尺寸增加，也会影响计算速度； B 计算中需要弛豫的层数； C 不同层数对你要计算的性质会产生影响，比如表面能； D 不同晶面需要的层数也不尽相同，一般开放的表面需要更多的层数； E 根据自己感兴趣的性质，选择合适的层数，也就是需要你去测试一番。
1）POSCAR第7行后面加一行加S， Sed –i ‘7a S’ POSCAR 原子坐标加F/T。 Sed –i ’10,11s/$/ F F F/g’ POSCAR Sed –i ’12,13,14s/$/ T T T/g’ POSCAR vaspkit –task 402/403 或者ctrl+v,选中列-大写A-输入 T T T - ESC. 2）POTCAR和KPOINTS不变。 INCAR设置：IBRION=1、2，EDIFF=1E-6，EDIFFG=-0.03（比块体小一点），ISIF=2， 提交任务。
5 表面能计算
表面能是创造物质表面时，破坏分子间化学键所需消耗的能量。在固体物理理论中，表面原子比物质内部的原子具有更多的能量，因此，根据能量最低原理，原子会自发的趋于物质内部而不是表面。表面能的另一种定义是，材料表面相对于材料内部所多出的能量。把一个固体材料分解成小块需要破坏它内部的化学键，所以需要消耗能量。如果这个分解的过程是可逆的，那么把材料分解成小块所需要的能量就和小块材料表面所增加的能量相等。但事实上，只有在真空中刚刚形成的表面才符合上述能量守恒。因为新形成的表面是非常不稳定的，它们通过表面原子重组和相互间的反应，或者对周围其他分子或原子的吸附，从而使表面能量降低。 γ=($E_{surf}$ - N * $E_{bulkatom}$)/(2A)
tamas@tamas-PC:~/Desktop/scan$ python scan_adsorption_energy.py -g Input which atom you want to move to scan adsorption energy–>-1 脚本提示我们移动第几个原子产生不同的吸附结构，我们可以输入97或者-1，移动最后一个Li。这个Cu(111)面是(331)的超胞，因此我们最需要取左下角的1/9的小区域撒点即可，然后程序会将小区域内的数据点进行平移得到整个表面的数据。如何根据不同的表面调整这个参数以及撒点的密度呢？可以使用-u参数进行自定义设置。
Input how many supercells in x axis–>3 Input how many supercells in y axis–>3 Input interporate how many points for each unit cell in x axis,even number are suggested–>2 Input interporate how many points for each unit cell in y axis,even number are suggested–>2 command to submit jobs,e.g. qsub vasp.pbs–>qsub vasp.pbs Input which atom you want to move to scan adsorption energy–>-1 手动选择x方向的超胞数3，y方向的超胞数3,在小区域内的x方向的撒点数，y方向的撒点数（为了包含中心点，建议选择偶数）和提交排队脚本的命令（我的是qsub vasp.pbs）。脚本会自动产生4个结构（2*2），并会将每个计算所需的文件都复制到相应的文件夹内，并提交排队脚本进行计算。同时脚本会写出一个grid.info里面包含了超胞数，选择移动的原子以及晶格常数信息。 在所有的计算完成后，同样使用该脚本提取数据。
cd / //切换到根目录 cd /bin //切换到根目录下的bin目录 cd ../ //切换到上一级目录 或者使用命令：cd .. cd ~ //切换到home目录 cd - //切换到上次访问的目录 cd xx(文件夹名) //切换到本目录下的名为xx的文件目录，如果目录不存在报错 cd /xxx/xx/x //可以输入完整的路径，直接切换到目标目录，输入过程中可以使用tab键快速补全
2 查看目录（ls）
{.line-numbers}
1 2 3 4
ls //查看当前目录下的所有目录和文件 ls -a //查看当前目录下的所有目录和文件（包括隐藏的文件） ls -l //列表查看当前目录下的所有目录和文件（列表查看，显示更多信息），与命令"ll"效果一样 ls /bin //查看指定目录下的所有目录和文件
INCARin ** STOPCARin stout out POTCARin ** KPOINTSin ** IBZKPT out POSCARin ** CONTCAR out CHGCARin/out CHG out WAVECARin/out TMPCARin/out EIGENVAL out DOSCAR out PROCAR out OSZICAR out PCDAT out XDATCAR out LOCPOT out ELFCAR out PROOUT out
