2.2 AMSET

发表于 2023-05-09 更新于 2024-06-05 分类于 DFT学习

简介

https://github.com/hackingmaterials/Amset
https://hackingmaterials.lbl.gov/amset/
https://www.nature.com/articles/s41467-021-22440-5.pdf

劳伦斯·伯克利国家实验室Alex M. Ganose和Anubhav Jain开发了一种根据第一性原理输入来计算固态半导体和绝缘体载流子散射率的有效计算方法(AMSET)。本方法扩展了现有的极性和非极性电子-声子耦合、离子化杂质和基于各向同性能带结构的压电散射机制的公式以支持高度各向异性的材料。他们通过计算23种半导体的电子传输特性来测试这个公式，包括48个大原子$CH_3NH_3PbI_3$混合钙钛矿，并将结果与实验测量结果和更详细的散射模拟进行比较。相对于实验的斯皮尔曼迁移率秩系数（rs = 0.93）与使用恒定弛豫时间近似值（rs = 0.52）获得的结果相比有明显提高。他们发现，他们的方法可提供与最新技术类似的精度，而计算成本约为1/500，因此可在高通量计算工作流程中使用，以精确筛选载流子迁移率、寿命和热电功率。

Inputs obtainable from first-principles calculations. The primary input for AMSET is an uniform band structure calculation.
Scattering rates calculated in the Born approximation using common materials properties such as phonon frequencies and dielectric constants.
Transport properties calculated through the Boltzmann transport equation.
Efficient implementation that can run on a personal laptop.

安装

依赖

pymatgen==2023.3.23
scipy==1.10.1
monty==2022.9.9
matplotlib==3.7.1
BoltzTraP2==22.12.1
tqdm==4.65.0
tabulate==0.9.0
memory_profiler==0.61.0
spglib==2.0.2
click==8.1.3
sumo==2.3.6
h5py==3.8.0
pyFFTW==0.13.0
interpolation==2.2.4
numba==0.56.4

conda environment

pip install amset

Developer Installation

{.line-numbers}

1
2
3

git clone https://github.com/hackingmaterials/amset.git
cd amset
pip install -e .

准备工作

Structural relaxation:

“tight” calculation settings

{.line-numbers}

ADDGRID = True  
EDIFF = 1E-8   
EDIFFG = -5E-4  
PREC = Accurate

Dense uniform band structure and wave function coefficients:

2倍的K点密度，HSE06 is required，输出波函数，

wave ``` 得到

```wavefunction.h5 file ```
### Elastic constants：
"tight"relax   
 ```javascript  {.line-numbers}  
 NSW = 1  
 IBRION = 6

Deformation potentials：

amset deform create
single point calculation，

{.line-numbers}

1 2	NSW = 1 ICORELEVEL = 1

amset deform read undeformed def-1 def-2 def-3
产生deformation.h5 file。

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

{.line-numbers}

1
2
3

NSW = 1
IBRION = 8
LEPSILON = True

amset phonon-frequency
读取DFPT的 vasprun.xml 获得 dielectric constants 和 polar phonon frequency。

运行

有两种方式，Python脚本
或者通过读取settings.yaml和输入文件加命令行运行

脚本

cp path/example/Si/Si.py .
脚本如下：

{.line-numbers}

import warnings
from amset.core.run import Runner
from amset.plot.rates import RatesPlotter
warnings.simplefilter("ignore")
settings = {
   #general settings
   "scattering_type": ["IMP", "ADP"],
   "doping": [-1.99e14, -2.20e15, -1.72e16, -1.86e17, -1.46e18, -4.39e18],
   "temperatures": [300],
   "bandgap": 1.14,
   #electronic_structure settings
   "interpolation_factor": 50,
   #materials properties
   "deformation_potential": "deformation.h5",
   "elastic_constant": [
       [144, 53, 53, 0, 0, 0],
       [53, 144, 53, 0, 0, 0],
       [53, 53, 144, 0, 0, 0],
       [0, 0, 0, 75, 0, 0],
       [0, 0, 0, 0, 75, 0],
       [0, 0, 0, 0, 0, 75],
       [0, 0, 0, 0, 0, 75],
   ],
   "static_dielectric": [[11.7, 0, 0], [0, 11.7, 0], [0, 0, 11.7]],
   "high_frequency_dielectric": [[11.7, 0, 0], [0, 11.7, 0], [0, 0, 11.7]],
   #performance settings
   "write_mesh": True,
}

if __name__ == "__main__":
   runner = Runner.from_vasprun("vasprun.xml", settings)
   amset_data = runner.run()

   plotter = RatesPlotter(amset_data)
   plt = plotter.get_plot()
   plt.savefig("Si_rates.png", bbox_inches="tight", dpi=400)

python Si.py 即可

命令行设置

需要准备settings.yaml文件：以Si为例：

{.line-numbers}

 # general settings
 scattering_type: [IMP, ADP]
 doping: [-1.99e14, -2.20e15, -1.72e16, -1.86e17, -1.46e18, -4.39e18]
 temperatures: [300]
 bandgap: 1.14

 # electronic_structure settings
 interpolation_factor: 50

 # materials properties
 deformation_potential: deformation.h5
 elastic_constant:
  - [144,  53,  53,  0,  0,  0]
  - [ 53, 144,  53,  0,  0,  0]
  - [ 53,  53, 144,  0,  0,  0]
  - [  0,   0,   0, 75,  0,  0]
  - [  0,   0,   0,  0, 75,  0]
  - [  0,   0,   0,  0,  0, 75]
static_dielectric:
  - [11.7, 0, 0]
  - [0, 11.7, 0]
  - [0, 0, 11.7]
high_frequency_dielectric:
  - [11.7, 0, 0]
  - [0, 11.7, 0]
  - [0, 0, 11.7]
 #performance settings
 write_mesh: true

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.

设置相关参数：

General settings

scattering_type：Default: auto ；单独设置： ADP,IMP,POP 。
ADP (acoustic deformation potential scattering)
IMP (ionized impurity scattering)
PIE (piezoelectric scattering)
POP (polar optical phonon scattering)
CRT (constant relaxation time)
MFP (mean free path scattering)
doping：负值代表电子，正值代表空穴
Default: ['1.e15', '1.e16', '1.e17', '1.e18', '1.e19', '1.e20', '-1.e15', '-1.e16', '-1.e17', '-1.e18', '-1.e19', '-1.e20']
单独设置有两种： 1E13,1E14,1E15,1E16,1E17,1E18,1E19,1E20 或者 1E13:1E20:8
temperatures：Default: [300]
单独设置：300,400,500,600,700,800,900,1000 或者300:1000:8
interpolation_factor：Default: 10
控制插值密度，最终的k点等于能带计算的k点乘上插值因子。
wavefunction_coefficients：Default: wavefunction.h5
波函数系数读取的方式，产生该文件可以通过amset wave获得。
bandgap：
设置该参数会自动采用scissor修正，因此不能与scissor选项共用，且对金属不起作用。
scissor：
正值代表带隙打开，负值代表缩小，对金属无效。
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_projections：Default: False
不建议设置。该参数使用投影计算波函数交叠，往往导致差的性能，使用时设置LORBIT = 11
free_carrier_screening：Default: False
控制是否屏蔽极化光学声子和压电散射率。这可以导致在高载流子浓度下散射率的大幅度降低。

Material settings

high_frequency_dielectric：Required for: POP, PIE
static_dielectric：Required for: IMP, POP
elastic_constant：Required for: ADP, PIE
deformation_potential：Required for: ADP
piezoelectric_constant：Required for: PIE
defect_charge：Required for: IMP。Default: 1
pop_frequency：Required for: POP amset phonon-frequency获得。
compensation_factor：Required for: IMP Default: 2
mean_free_path：Required for: MFP
constant_relaxation_time：Required for: CRT 不建议和其他散射机制设置。

Performance settings：控制速度和精度，一般默认。

energy_cutoff：Default: 1.5
fd_tol：Default: 0.05. 0-1之间，越小计算K点越多。
dos_estep:Default: 0.01.越小越好越贵。
symprec：Default: 0.01。
nworkers：Default: -1。代表用所有的机器。其他设置代表几个。
cache_wavefunction：Default: True。加速计算，增加内存。如果内存不够关闭。

Output settings：输出文件设置

calculate_mobility：Default: True。是否计算n/p载流子迁移率。
separate_scattering_mobilities：Default: True。是否计算单独的散射机制的迁移率。
write_mesh：Default: False。是否将全K点结果写入磁盘。对于大的interpolation_factor可以打开。
file_format：Default: json。输出文件的格式：json, yaml, and txt.对于write_mesh=True不能用TXT。
write_input：Default: False。是否将输入设置写入amset_settings.yaml.
print_log：Default: True。是否输出日志。

Summary of scattering rates

Mechanism	Code	Requires	Type
Acoustic deformation potential scattering	ADP	n- and p-type deformation potential, elastic constant	Elastic
Piezoelectric scattering	PIE	high-frequency dielectric constant, elastic constant, piezoelectric coefficient (e)	Elastic
Polar optical phonon scattering	POP	polar optical phonon frequency, static and high-frequency dielectric constants	Inelastic
Ionized impurity scattering	IMP	static dielectric constant	Elastic

其他功能

you can use amset which relies on Fourier interpolation of the band
structure using BoltzTraP2. There are two options.

Use amset eff-mass: This calculates the conductivity effective mass as defined here: https://www.nature.com/articles/s41524-017-0013-3
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.

2.1 电子性质1

发表于 2023-05-09 更新于 2025-07-10 分类于 DFT学习

自洽计算

VASP做静态自洽的目的主要是为了得到好的电子密度和波函数文件，为了使后续的性质计算可以读取电子密度和波函数，进而增加收敛速度。
首先进行 高精度 的结构优化。然后复制文件， cp CONTCAR POSCAR 做一次自洽计算。
自洽计算要求有密的k网格点，在计算条件允许的情况下要求自洽的k网格点大致为优化时的2倍左右。
修改标签

{.line-numbers}

IBRION=-1
NSW=0
LWAVE  = .T.        (Write WAVECAR )
LCHARG = .T.        (Write CHGCAR )

这一步完成后的输出文件CHGCAR即是电荷密度，导入VESTA直接查看。

bader电荷

自洽时INCAR再加一个参数：

{.line-numbers}

1
2
3

IBRION=-1
NSW=0
LAECHG=.T.

输出文件：
AECCAR0 全电子中的芯电子部分（POTCAR）
AECCAR1 自洽计算前的交叠原子电荷（POTCAR，bader 分析不用）
AECCAR2 全电子中的价电子部分（自洽计算后）
CHGCAR 自洽计算后赝电子

输入：

{.line-numbers}

1 2	chgsum.pl AECCAR0 AECCAR2 bader CHGCAR -ref CHGCAR_sum

得到两个文件：

BCF.dat 电荷极大值序号和坐标
体积内的 Bader 电荷（积分）
距离极大值最近的原子和距离
ACF.dat 原子序号和坐标
原子电荷（极大值的电荷之后）
到零通量面最小距离
原子体积

李强博士《learn vasp the hard way》https://www.bigbrosci.com/2011/12/24/A08/ Bader电荷与原子电荷着色：
https://www.bigbrosci.com/2011/12/23/A07/ https://mp.weixin.qq.com/s/BNzjhz8SI_HXkaEnyL29Bg

成键反键

自洽时添加：

{.line-numbers}

NBANDS=96
ISYM=-1
IBRION=-1
NSW=0
LWAVE = T

vi lobsterin
原子编号改为自己的，例如

lobsterin：

COHPstartEnergy  -14
COHPendEnergy     6   
basisSet Bunge
basisfunctions Ga 4s 4p 3d
basisfunctions As 4s 4p

cohpbetween atom 1 and atom 2

前两行代表提取的能量范围；第三行代表选择的基组，如果不设置(空行)就是默认基组，如果要覆盖默认基集，可使用此关键字进行覆盖。它目前支持bunge、koga和pbeVaspFit2015；第四行代表为每个元素指定所使用的基函数；最后一行代表选择一对两个原子进行键合分析。

lobster < lobsterin

1、删除存在的WAVECAR文件，做单点计算，计算需要保存波函数文件；
2、做完单点计算在当前目录下运行lobster软件；
3、用wxDragon软件打开COHPCAR.lobster文件。
4、COHPCAR.lobster 文件包含COHP结果，第二列为pCOHP，将第二列乘-1可得-pCOHP
https://github.com/JaGeo/LobsterPy
https://zhuanlan.zhihu.com/p/470592188

注：

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

不要使用US-PP超软赝势，使用PAW赝势；请使用vasp_std版本，不要使用vasp_gam
；不支持处理对称:ISYM=0或ISYM=-1
由于需要把基于平面波的波函数投影到原子轨道的基函数上，所以计算的能带的数量要多于投影的基函数，VASP计算的时候，默认的NBANDS达不到要求，需要加大NBANDS，只要查看OUTCAR里面的NBANDS多少，grep NBANDS OUTCAR；在INCAR里面，NBANDS关键词设置成1.5倍NBANDS，足以满足计算的投影要求，在计算过程中一定要记住设置足够的NBANDS，否则会出现报错，一直不能正常算出COHP
cohpGenerator 这里通过距离或者元素类别告诉想要计算什么原子之间的COHP；cohpGenerator from 1.5 to 2.3 type Fe type N orbitalwise #输出Fe-N原子对距离在1.5到2.3的所有数据之和。
gaussianSmearingWidth 高斯展宽，如果VASP中用的ISMEAR=0，那么gaussianSmearingWidth的数值和SIGMA一样，默认值是0.2 eV，但是过大了，所以需要减小。设置格式：gaussianSmearingWidth 0.05

分析：

参考：
https://zhuanlan.zhihu.com/p/1920205994354521233

COHP（晶体轨道哈布顿居数）：量化键合强度与稳定性变化，描述的是特定能级上，两个原子轨道之间的相互作用的强度和性质

负值越大，键越稳定，可以用来分析掺杂对键能的影响（对结构稳定性的影响），负值代表成键，增强键强度
正值代表反键，削弱键强度，键易断裂
0值代表非键态

ICHOP：CHOP从-infnite到费米能级的积分，因为费米能级以下的电子是参与化学键形成的主要电子；注意进行不同体系或者不同键的ICOHP比较时，要保持积分区间一致。

ICOHP将不同能量下的原子轨道间的相互作用的贡献累加在一起，可以代表键的总强度。
ICHOP<0，表示成键，绝对值越大成键越稳定
ICHOP>0，表示反键，绝对值越大成键越不稳定

COHPCAR.lobster文件内容

第一列为能量，第2，3，4，5列为自旋向上的相互作用，其中第2，3列分别是平均相互作用的cohp和icohp，第4，5列分别是指定相互作用的chop和ichop

那么第6，7，8，9列就是对应自旋向下的。

绘制cohp-energy图像时可以选择第一列和第四列，一般绘制-cohp~energy图像。

ICOBILIST.lobster, COBICAR.lobster输出文件，包含共价键指标: 晶体轨道键指数(Crystal Orbital Bond Index, COBI)

ICOBI: 0(离子键)~1(共价键) per bond, 自旋极化体系分spin1,2，具有加和性

MadelungEnergies.lobster, SitePotentials.lobster输出文件，包含离子键指标: 马德隆能(Madelung Energy, in eV per unit cell)

ELF电子局域函数计算

电子局域函数（Electron Localization Function—ELF）是研究电子结构的手段之一。用来描述以某个位置处的电子为参考，在其附近找到与他同自旋的电子的概率，可以表征这个作为参考的电子的局域化程度，也是一种描述在多电子体系中的电子对概率的方法。
ELF的取值范围是[0，1]。
ELF=1 对应完全局域化；
ELF=1/2，对应类电子气型的成对概率；
ELF=0，对应电子完全离域化。

计算流程：做自洽的时候在INCAR中加入：

{.line-numbers}

PREC=A
LELF=T
IBRION=-1
NSW=0

运行vasp得到ELFCAR文件。

将ELFCAR拖入Vesta http://www.jp-minerals.org/vesta/en/download.html

点击 Utiltties - 2D Data Display。
点击Slice，确定切面和fractional coordinates。
导出图片：之后可以用PS工具将原子模型加入。

费米面

费米面计算，这里我们以Cu为例，见vaspkit/examples/Cu_fermi_surface。
第一步，准备好POSCAR，注意一定是原胞（primitive cell ）和用于静态计算的INCAR文件；
第二步，运行VASPKIT，输入261命令，产生用于计算计算费米面的KPOINTS和POTCAR文件，
第三步，提交VASP作业；第四步，待VASP计算完成后，再次运行VASPKIT，输入262命令，生成包含费米面数据的FERMISURFACE.bxsf文件；第五步，运行XcrySDen软件。
发现第五条能带穿过费米面，因此我们选择Band Number 5，然后点击Selected，然后得到Cu的费米面。推荐使用FermiSurfer可视化费米面。
对于半导体不存在费米面，如果想可视化导带或价带中某个能级处的等能面，新建FERMI_ENERGY.in文件然后在该文件第二行（第一行为注释行）中人为指定能级大小，注意该值一定要保证穿过价带或导带。接下来调用vaspkit-252或263命令提取可视化。

STM模拟

参考：https://zhuanlan.zhihu.com/p/538950799

https://zhuanlan.zhihu.com/p/667558468

https://www.vasp.at/wiki/index.php/STM_of_graphite

首先做一个自洽计算，输出 WAVECAR 和 CHGCAR 文件。

加入一些特殊参数计算，输出我们想要的 PARCHG 文件。

LPARD = .TRUE.
LSEPK = .FALSE.
LSEPB = .FALSE.
NBMOD = -3
EINT = -0.1 0.1

用p4vasp软件打开PARCHG文件或者vaspkit 后处理。由于 VASP 在对于远离原子表面时波函数的衰减指数描述不明确，所以模拟探针高度的时候不应该太高，否则会出现非物理的现象，比如 7 埃，就算比较远了。所以一般高度的选择要接近材料表面，这里选择了 0.5 埃。而原胞沿着 ab 平面重复的标准是 20 埃，也就是说原胞晶格 * 重复单元 = 20 埃就可以。vaspkit 里选择的是 Tersoff-Hamman 方法，现在大多数模拟 STM 的都是选择这个方法。
vaspkit处理后的输出文件包括X.grd, Y.grd, STM_*.grd三个文件，可以用以下脚本画图。
Usage: python plot.pyimport numpy as npimport matplotlib as mpl#mpl.use(‘Agg’) #silent modefrom matplotlib import pyplot as plt
x_mesh=np.loadtxt(‘X.grd’) y_mesh=np.loadtxt(‘Y.grd’) v_mesh=np.loadtxt(‘STM.grd’)
fig = plt.figure()cmap = plt.cm.get_cmap(“Greys_r”) plt.contourf(x_mesh,y_mesh,v_mesh,cmap=cmap)
plt.axis(‘equal’)plt.axis(‘off’)#plt.show()plt.savefig(“STM.jpg”,dpi=300)
或者用以下脚本处理后，得到data.txt文件，将该文件导入origin画图。Usage: sh get_data.shawk ‘{for(i=1;i<=NF;i++) print $i}’ X.grd >xawk ‘{for(i=1;i<=NF;i++) print $i}’ Y.grd >yawk ‘{for(i=1;i<=NF;i++) print $i}’ STM_0.50.grd >spaste x y s >data.txt

态密度

自洽计算之后，读取产生的CHGCAR、WAVECAR
上面例如bader电荷、ELF标签都关闭。
INCAR中修改：

{.line-numbers}

ISTART=1
ICHARG=11
ISMEAR=-5
NEDOS=2001
LORBIT=11

其他参数和自洽时保持一致。
运行后vaspkit的111/112/113/114命令处理即可得到。

分波态密度和局域态密度

态密度很多分析和能带的分析结果可以一一对应，很多术语也和能带分析相通。但是因为它更直观，因此在结果讨论中用得比能带分析更广泛一些。在电子能级为准连续分布的情况下，单位能量间隔内的电子态数目，即为态密度。能态密度与能带结构密切相关，是一个重要的基本函数。固体的许多特性，如电子比热、光和X射线的吸收和发射等，都与能态密度有关。在VASP中，LORBIT=10计算的就是LDOS，也就是每个原子的局域态密度 (local DOS)，是分解到每个原子上面的s,p,d；LORBIT=11，计算的就是PDOS,投影态密度(projected DOS)或分波态密度(partial DOS)，不仅分解到每个原子的s,p,d,而且还进行px,py,pz分解。

在 /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，显而易懂。

电荷密度差

Vaspkit 可以方便的处理 VASP 输出的电荷密度文件 (CHGCAR)，做电荷密度差分处理，VTST 脚本或者 VESTA 也可以做，但是不如 vaspkit方便灵活。
注：此文中讲述的功能同样适用于任何的实空间函数形式的文件处理，比如静电势差值图 (LOCPOT)，电子定域化函数差值图 (ELFCAR)，能带分解的电荷密度插值图 (PARCHG) 等等。
电荷密度差分（charge density difference）是研究电子结构的重要手段之一。可以直观的得到两个片段相互作用后的电子流向，原子形成分子过程中电子密度的变化、探究化学键的本质。电荷密度差分有以下几种形式：
（1）体系的电荷密度减去组成它的两个或几个片段的密度，也是最常见的电荷密度差分形式：

（2）自洽计算收敛以后体系的电荷密度减去该原子构型下每个原子的球对称的电荷密度(即初猜电荷密度)，也称为变形电荷密度(Deformation charge density)

（3）在某个状态的密度减去这个体系在另外一个状态的密度。比如：外加电场作用下的电荷密度减去没有外势场的电荷密度。再比如：激发态的密度减去它在基态时的密度。

以上三种电荷密度差，无论哪一种计算都需要保证原子坐标坐标必须一致！并且保证晶胞参数和 FFT-mesh 的格点密度 (NGXF, NGYF, NGZF) 一致！

(1) 的计算只需要优化 AB 的结构，A、B片段保持 AB 结构中原子坐标，不能再对片段进行单独优化。
(2) 只要优化自洽计算时候的几何结构，生成原子的球对称的电荷密度时候 (ICHARG = 12) 保持坐标不变。
(3) 也需要共用其中一个状态的几何构型（比如激发态或基态的几何构型）。

CHGCAR文件格式

CHGCAR 是包含电子密度信息的格点文件，对于自旋非极化体系（ISPIN = 1）计算只包含电荷密度，对于自旋极化体系（ISPIN = 2）计算还包含自旋电子密度。可以用 VESTA，Jmol 等程序打开。CHGCAR 的第一部分和 POSCAR, CONTCAR 的格式是完全一样的，它包含了最终结构的晶格矢量，原子核坐标等信息。紧接着是三个数字，这三个数是实空间函数的网格密度，对应于 NGXF，NGYF，NGZF 三个变量。然后是电荷密度信息 ρ(r) * Vcell。一共有 NGXF * NGYF * NGZF 个数值。比如：下面是一个 Al2O3 晶胞的 CHGCAR 文件

Al2O3 Cell opt                          
  1.00000000000000     
    6.138452   0.000551   -0.011568
    5.434311   2.854620   -0.011568
  -1.324042   -0.326661   5.501689
  O   Al
    6     4
Direct
0.160314 0.160314 0.108934
0.839686 0.839686 0.891066
0.495028 0.495028 0.257399
0.504972 0.504972 0.742601
0.826429 0.826429 0.432705
0.173571 0.173571 0.567295
0.090130 0.090130 0.795615
0.909870 0.909870 0.204385
0.342109 0.342109 0.683236
0.657891 0.657891 0.316764
100 100   96
0.43432567054E+01 0.44027037758E+01 0.45822813500E+01 0.48861003640E+01 0.53223812855E+01
0.59054783755E+01 0.66583226816E+01 0.76145085758E+01 0.88190771181E+01 0.10327304760E+02
...

Vaspkit计算几个片段的电荷密度差

以 CO 分子吸附在 Ni(100) 表面为例：
步骤一：优化 CO/Ni(100) 的结构；
步骤二：分别计算 Ni(100) 和 CO 的单点能 CO 和 Ni(100) 片段的坐标从 CO/Ni(100) 的 CONTCAR 里直接截取，不要再结构优化！！计算时也要保证三次自洽计算所采用的 FFT mesh 一致（NGXF,NGYF,NGZF）
步骤三：用 vaspkit 314 功能做电荷差分。依次输入三个片段的 CHGCAR 路径:

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 路径。

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!

步骤四：得到 CHGDIFF.vasp 文件，导入到 VESTA 里做图。

Vaspkit计算外加电场下的电荷密度差

外加电场作用下的电荷密度减去没有外势场的电荷密度。
以 InSe 二维单层材料为例：
步骤一：先优化没有外加电场的结构。
步骤二：在同样结构下计算外加电场下做单点自洽计算。添加关键词，EFIELD 控制电场力的大小（eV/Angstrom）F=qE，对应电场的单位是 E=F/q，场强单位是 V/Angstrom

1
2
3

EFIELD = 0.05
IDIPOL = 3
LDIPOL = .TRUE.

步骤三：运行 Vaspkit 314 功能做电荷差分。依次输入二个 CHGCAR 路径。

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 里做图。

电荷密度差作图

得到的电荷密度差文件必须作成图才便于被直观地考察，常用的图分为三类：
（1）3D 等值面图，isosurface.
（2）2D 切面图.
（3）1D，平面平均的数值，planar averaged charge.
通常这几种图可以配合使用讨论。等值面图和切面图用 VESTA 画最方便，平面平均图用 vaspkit 生成方便。得到的 CHGDIFF.vasp 用 VESTA 打开，isosurface 青色部分电荷密度减小，黄色密度电荷密度增加。注意 VESTA 对电荷密度的默认单位是 e/Bohr3。
vaspkit 319 功能可以把VASP的格点文件转化成 .cube的方法，.cube 在 VMD 里作图可以得到更好看的效果：

能带计算

参考链接：
https://yh-phys.github.io/2019/10/04/vasp-band-structure-dos/
自洽计算之后，读取产生的CHGCAR、WAVECAR

{.line-numbers}

1 2	ISTART=1 ICHARG=11

其他参数和自洽时保持一致。

能带计算的KPOINTS文件不同于前面，称为高对称点。
学术之友微信公众号历史记录搜索文章《五种方法生成能带结构计算的高对称点》
1、http://cryst.ehu.es/
2、https://www.materialscloud.org/work/tools/seekpath
3、http://www.aflowlib.org/aflow-online (已更新)
4、http://pymatgen.org/ (pymatgen)
5、http://www.xcrysden.org/ (XCrySDen)
Accurate DOS and Band-structure calculations：
http://cms.mpi.univie.ac.at/vasp/vasp/Accurate_DOS_Band_structure_calculations.html

可以通过vaspkit的303命令产生KPATH.in文件，cp KPATH.in KPOINTS即可。
计算完成后通过vaspkit的211命令提取能带。
注意能带路径生成操作需在自洽计算前完成。

1.准备POSCAR，调用vaspkit- 303（体相材料）或-302（二维材料）得到KPATH.in和PRIMCELL.in文件；对于二维体系，需要检查POSCAR文件的真空层是否沿z方向，如果没有，可调用vaspkit-923或vaspkit-407强制真空层沿z方向。另外，如果需要改变识别结构对称性精度（symprec默认值为1E-5），可通过vaspkit -task 302 -symprec 1E-6实现或者在~/.vaspkit中修改SYMPREC默认值。
2.执行cp PRIMCELL.vasp POSCAR命名。（如果PRIMCELL.vasp 和POSCAR不一样，并且要读取电荷密度和波函数，要从SCF计算保持结构一致。）KPATH.in推荐能带路径只针对于PRIMITIVE CELL，缺少这一步，你可能得到错误的结果。如果有必要，比较KPATH.in文件中的能带路径是否与在线能带路径产生工具SeeK-Path产生的一致，包括比较PRIMCELL.vasp和HIGH_SYMMETRY_POINTS文件。需要指出的是SeeK-Path只用于体相结构能带路径的产生。自1.3.0版本，对于二维体系不需要执行cp PRIMCELL.vasp POSCAR。

投影电荷密度

复制能带结构文件，
vi INCAR

{.line-numbers}

1
2
3

增加标签 LPARD=.T.
增加标签 IBAND=16 17
增加标签 LSEPB=.T.

LPARD：它的取值为.TRUE.或.FALSE.，它的默认值是.FALSE.，当为.TRUE.时，表示读入自
洽收敛的 CHGCAR 和 WAVECAR 并进行 partial charge density计算。
IBAND：它的取值就是你想要计算的第几条能带或哪几条能带(比如要计算第 4、 5、 6 条能带，那么就设置 IBAND = 4 5 6)，此时 NBMOD的值就是所要计算的能带的条数。
EINT：另一种方式指定所要计算的能带，它是指定计算某能量范围的partial charge density，一般是设置为两个实数，比如EINT= 4.00 5.00，此时NBOMD的值设置为-2。如果只设置了一个数，那么表示计算从 EINT 到费米能级这个范围内的，此时NBOMD 的值为-3。
KPUSE：指定所要计算的 k 点(哪个或哪几个)。
LSEPB：指定是不是要把计算的partial charge density按每个带分别写到各自对应的文件
PARCHG.nb (设置为.TRUE.)中，还是把它们合并写到一个文件中(相当于把各个带对应的partial charge density 加起来，设置为.FALSE.)，默认值为.FALSE.。
LSEPK：指定是不是把要计算的partial charge density按每个k分别写到各自对应的文件PARCHG.nk (设置为.TRUE.).

三个方法：

1、读取自恰计算的波函数，KPOINTS文件采用特定的k点

{.line-numbers}

#partial charge density
LPARD = .TRUE. 
IBAND = 17 18 # 价带顶和导带底
KPUSE = 1
LSEPB = .TRUE.
LSEPK = .TRUE.

{.line-numbers}

CsPbI3-cubic
1
r
0.500000 0.500000 0.500000 1.0

http://muchong.com/html/200906/1404880.html
2、读取自恰计算的波函数，KPOINTS文件采用自恰计算的k点

{.line-numbers}

#partial charge density
LPARD = .TRUE. 
IBAND = 17 18
KPUSE = 35
LSEPB = .TRUE.
LSEPK = .TRUE.
NELM = 200
ICHARG = 1

在scf自恰计算OUTCAR中
找到对应导带底和价带顶
对应的k-point编号

{.line-numbers}

CsPbI3-cubic
0
Gamma
8 8 8
0.0 0.0 0.0

3、读取能带计算的波函数，KPOINTS文件采用能带计算的k点

{.line-numbers}

#partial charge density
LPARD = .TRUE. 
IBAND = 17 18
KPUSE = 80 81
LSEPB = .TRUE.
LSEPK = .TRUE.
NELM = 200
ICHARG = 1

KPOINTS和能带计算一样。在band能带计算OUTCAR中找到对应导带底和价带顶对应的k-point编号。

轨道投影态密度和能带。

{.line-numbers}

1 2	LORBIT=11 NEDOS=2001

提取DOS：vaspkit -task 111/113
提取band：211/212/213/214
如果3D能带图不够平滑，可通过能带插值进行改善。
VASPKIT升级到1.3.0或更新版本，在~/.vaspkit文件中设置以下三个参数：

{.line-numbers}

1
2
3

GET_INTERPOLATED_DATA      .TRUE.     是否插值
INTERPOLATION_SPACING        0.04     插值间距，决定平滑程度
INTERPOLATION_METHOD        'cubic'   插值方法

HSE06计算电子结构：

结构优化（PBE）
静态自洽scf计算（PBE）KPOINTS和POTCAR文件不变，写波函数和电荷密度
(可选）静态自洽scf计算（HSE06）KPOINTS，POSCAR和POTCAR文件不变
DOS/band计算（HSE06）vaspkit -task 303
vaspkit→251→2→0.03→0.03生成KPOINTS，运行VASP
vaspkit-252
http://vaspkit.cn/index.php/29.html
PBE优化结构
HSE自洽计算(得到波函数)

HSE能带计算/态密度计算

{.line-numbers}

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}

 ISPIN  =  2
 NELMIN =  6            (Min electronic SCF steps) 
 LSORBIT    = .TRUE.    (Activate SOC)
 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    =  0 0 1      (Direction of the magnetic field)
 #MAGMOM   =  0 0 3    (Set this parameters manually, Local magnetic moment parallel to SAXIS, 3*NIONS*1.0 for non-collinear magnetic systems)
 NBANDS   =  32         (Set this 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

代码片段：可切换语言，无法单独设置文字格式

注意，需要使用重新编译的vasp_ncl来运行计算。另外，最好比较0d0L●S，1d0L●S和vasp_std，vasp_ncl计算结果。

如果只想给一些元素加SOC，而其他元素不加该怎么办呢？拓扑材料计算互助群的Eurus-清华给出了VASP中给特定元素加自旋轨道耦合相互作用的方法（感谢Eurus-清华）。在原来VASP，src文件夹下的relativistic.F文件中注释掉46行的

! REAL(q), PARAMETER :: INVMC2=7.45596E-6

然后53行位置加入如下6-17行的内容，然后重新编译。

  INTEGER, PARAMETER :: LMAX=3, MMAX=LMAX*2+1      COMPLEX(q) DUMMY(MMAX,MMAX,3,LMAX)      COMPLEX(q) LS(MMAX,MMAX,4,LMAX)      REAL(q) SUM, SCALE

! add parameter REAL(q), PARAMETER :: INVMC2_ORIG=7.45596E-6 REAL(q) INVMC2! add end ! add control IF(ABS(Z-83).LT.1.0E-5)THEN INVMC2=INVMC2_ORIG ELSE INVMC2=INVMC2_ORIG*0.0001 END IF! add end
LS=(0._q,0._q) CALL SETUP_LS(1,THETA,PHI,DUMMY(1:3,1:3,1:3,1),LS(1:3,1:3,1:4,1)) CALL SETUP_LS(2,THETA,PHI,DUMMY(1:5,1:5,1:3,2),LS(1:5,1:5,1:4,2)) CALL SETUP_LS(3,THETA,PHI,DUMMY(1:7,1:7,1:3,3),LS(1:7,1:7,1:4,3))
原理：这样做的原理就是将除Z-83的元素外的其他元素缩小一万倍。注意：这个脚本的第12行, IF(ABS(Z-83))，表明除了83号元素打开SOC之外，其他所有的元素SOC都是相当于关闭的。如果想打开某个特定元素的SOC，就把83改成你想要的元素序号就行。
另外，其他元素SOC关闭方式就是将其缩小一万倍。所以使用这种方法，任然是非共线计算，如果要进行wannier90计算，只要按照正常的wannier90设置就可以。

HSE06+SOC band structure calculation

Step-I: HSE06 SCF calculation with SOC (INCAR):

{.line-numbers}

#SOC-related:
LSORBIT=.TRUE.                
SAXIS= 0 0 1  
ISYM=0                   

#HSE06-related:
LHFCALC = .TRUE.
HFSCREEN = 0.2
ALGO = Damped         #Damped/Normal/All
TIME = 0.4

Step-II: HSE06 NSCF calculation with SOC:

{.line-numbers}

#SOC-related:
LSORBIT=.TRUE.                
SAXIS= 0 0 1  
ISYM=0                  

#HSE06 related:
LHFCALC = .TRUE.               
HFSCREEN = 0.2
ALGO = Normal                  !use Normal
TIME = 0.4

能带反折叠

vaspkit:

准备POSCAR
运行vaspkit -task 302/303生成KPATH.in文件
cp PRIMCELL.vasp POSCAR
运行vasp优化结构，cp CONTCAR POSCAR
运行vaspkit -task 400/400生成超胞结构POSCAR，并执行cp TRANSMAT TRANSMAT.in
运行vasp优化含有缺陷的超胞结构，cp CONTCAR POSCAR
运行vaspkit -task 281生成KPOINTS文件
准备INCAR（静态计算，设置LWAVE = T，适当增加NBANDS值）并执行VASP计算
运行vaspkit -task 282提取effective band structure(画法类似于投影能带)

如果不会计算变换矩阵（TRANSMAT.in文件不存在时），可以执行cp PRIMCELL.vasp PRIMCELL.in，281 ，这时vaspkit分别读入超胞基矢（POSCAR）和原胞基矢（PRIMCELL.in）计算变换矩阵并生成TRANSMAT.in文件。
vaspkit/examples/band_unfolding/ebs.py画图。生成的ENERGY.grd，WEIGHT.grd和MOMENTUM.grd三个文件，可使用vaspkit/examples/band_unfolding/ebs_k_resolved_dos_plot_matlalb.m脚本画图。

bandup:

第一步 得到超胞CHGCAR:输入文件参考BANDUP的example的第一步，提交任务之后得到CHGCAR。
第二步 主要是为了得到超胞非自洽计算的K点路径，分别给出单胞和超胞的晶格坐标文件prim_cell_lattice.in与supercell_lattice.in。还有单胞的高对称K点路径。
下一步命令是bandup kpts-sc-get -no_symm_avg,这一步会生成KPOINTS_supercell.out文件后面会用到。-no_symm_avg这个标签是为了减少需要计算的k点个数，减小计算量和WAVECAR大小。
第三步 利用前面第一步生成 CHGCAR; POSCAR POTCAR与第一步的一样，K点从第二步那个复制过来,然后提交任务得到波函数WAVECAR
第四步就要开始做能带反折叠了，一句命令，bandup unfold -emin -13 -emax 6 -dE 0.050 -no_symm_avg
会生成一个unfolded_EBS_not-symmetry_averaged.dat文件，画出来就行了。
https://github.com/band-unfolding/bandup

http://staff.ustc.edu.cn/~zqj/posts/Band-unfolding-tutorial/
可以处理能带反折叠的软件汇总：
https://mp.weixin.qq.com/s/-5lN2lwxsALA-rhNQxK1dw
https://vaspkit.com/tutorials.html#band-unfolding
https://github.com/QuantumMaterialsModelling/bands4vasp
https://github.com/QuantumMaterialsModelling/UnfoldingPatch4vasp

3D能带

__第一步__，准备好石墨烯的POSCAR和用于静态计算的INCAR文件；
__第二步__：运行VASPKIT并选择231生成用于3D能带计算的KPOINTS文件。为了获得平滑的3D能带面，用于产生KPOINTS文件的分辨率值应该设置在0.005左右；
__第三步__，提交VASP作业；
__第四步__，待VASP计算完成后，再次运行VASPKIT，输入232或233命令。233命令可一次性输入包含价带顶的BAND.HOMO.grd和导带底的BAND.LUMO.grd文件的能带；232可以得到其它任意能带的计算数据，这里我们选择233。
__第五步__：运行python how_to_visual.py，绘制3D能带图(python2.7环境)。
如果3D能带图不够平滑，可通过能带插值进行改善。

有效质量：

除了通过能带，也可以通过调用VASPKIT-911命令来得到带边位置，接下来我们计算K点处沿着K -> Γ和K -> M方向的电子及空穴载流子的有效质量；
第一步：准备POSCAR文件以及VPKIT.in文件，VPKIT.in文件包含以下内容；

{.line-numbers}

1                  设置为1将会生成计算有效质量的KPOINTS文件，2则计算有效质量
6                  6表示我们在K附近取6个离散点用于多项式拟合
0.015              k-cutoff, 截断半径，典型值0.015 1/Å
2                  2 表示有两个有效质量任务
0.333333 0.3333333  0.000  0.000 0.000 0.000  K->Γ     计算K点处有效质量，沿着K指向Γ方向    
0.333333 0.3333333  0.000  0.500 0.000 0.000  K->M     计算K点处有效质量，沿着K指向M方向

第二步：运行vaspkit并选择912或者913产生KPOINTS，POTCAR文件及静态计算INCAR文件；
第三步：提交VASP作业；
第四步：计算完成后把VPKIT.in文件中第一行的1修改为2，然后再次运行vaspkit并选择913，得到结果：
为了得到可靠的有效质量计算值，括号里的精度至少要在1E-7这个量级，如果误差比较大，可以适当调整k-cutoff值，一般减少该值但不能太小。
手动vaspkit处理的能带：$m^∗=3.815/B$

GW方法

参考：https://www.vasp.at/wiki/index.php/Bandgap_of_Si_in_GW

步骤：

PBE基态计算
获得虚拟轨道（50-100/atom）

INCAR：

ALGO = Exact
NBANDS  = 64
LOPTICS = .TRUE. #得到WAVEDER文件
CSHIFT = 0.1
NEDOS = 2000
# you might try
#LPEAD = .TRUE.
   
ISMEAR =  0
SIGMA  =  0.05
EDIFF  = 1E-8

完成后：

1 2	cp WAVECAR WAVECAR.DIAG cp WAVEDER WAVEDER.DIAG

GW计算
读取上一步
INCAR：

# 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

在OUTCAR中quasi-particle (QP) energies.

QP shifts <psi_nk| G(iteration)W_0 |psi_nk>: iteration 1
 for sc-GW calculations column KS-energies equals QP-energies in previous step
 and V_xc(KS)=  KS-energies - (<T + V_ion + V_H > + <T+V_H+V_ion>^1  + <V_x>^1)

 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：

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.

wannier90.win:

num_wann=8
num_bands=8

exclude_bands = 9-64

Begin Projections
Si:sp3
End Projections

dis_froz_max=9
dis_num_iter=1000

guiding_centres=true

# Bandstructure plot 
#restart         =  plot
#bands_plot      =  true
#begin kpoint_path
#L 0.50000  0.50000 0.5000 G 0.00000  0.00000 0.0000
#G 0.00000  0.00000 0.0000 X 0.50000  0.00000 0.5000
#X 0.50000  0.00000 0.5000 K 0.37500 -0.37500 0.0000
#K 0.37500 -0.37500 0.0000 G 0.00000  0.00000 0.0000
#end kpoint_path
#bands_num_points 40
#bands_plot_format gnuplot xmgrace

begin unit_cell_cart
     2.7150000     2.7150000     0.0000000
     0.0000000     2.7150000     2.7150000
     2.7150000     0.0000000     2.7150000
end unit_cell_cart

begin atoms_cart
Si       0.0000000     0.0000000     0.0000000
Si       1.3575000     1.3575000     1.3575000
end atoms_cart

mp_grid =     4     4     4

begin kpoints
     0.0000000     0.0000000     0.0000000
     0.2500000     0.0000000     0.0000000
     0.5000000     0.0000000     0.0000000
     0.2500000     0.2500000     0.0000000
     0.5000000     0.2500000     0.0000000
    -0.2500000     0.2500000     0.0000000
     0.5000000     0.5000000     0.0000000
    -0.2500000     0.5000000     0.2500000
     0.0000000     0.2500000     0.0000000
     0.0000000     0.0000000     0.2500000
    -0.2500000    -0.2500000    -0.2500000
    -0.2500000     0.0000000     0.0000000
     0.0000000    -0.2500000     0.0000000
     0.0000000     0.0000000    -0.2500000
     0.2500000     0.2500000     0.2500000
     0.0000000     0.5000000     0.0000000
     0.0000000     0.0000000     0.5000000
    -0.5000000    -0.5000000    -0.5000000
     0.0000000     0.2500000     0.2500000
     0.2500000     0.0000000     0.2500000
    -0.2500000    -0.2500000     0.0000000
    -0.2500000     0.0000000    -0.2500000
     0.0000000    -0.2500000    -0.2500000
     0.0000000     0.5000000     0.2500000
     0.2500000     0.0000000     0.5000000
    -0.2500000    -0.2500000     0.2500000
    -0.5000000    -0.2500000    -0.5000000
     0.2500000     0.5000000     0.0000000
     0.2500000    -0.2500000    -0.2500000
    -0.5000000    -0.5000000    -0.2500000
     0.0000000     0.2500000     0.5000000
    -0.2500000     0.2500000    -0.2500000
    -0.2500000    -0.5000000    -0.5000000
     0.5000000     0.0000000     0.2500000
    -0.5000000    -0.2500000     0.0000000
     0.0000000    -0.5000000    -0.2500000
    -0.2500000     0.0000000    -0.5000000
     0.2500000     0.2500000    -0.2500000
     0.5000000     0.2500000     0.5000000
    -0.2500000    -0.5000000     0.0000000
    -0.2500000     0.2500000     0.2500000
     0.5000000     0.5000000     0.2500000
     0.0000000    -0.2500000    -0.5000000
     0.2500000    -0.2500000     0.2500000
     0.2500000     0.5000000     0.5000000
    -0.5000000     0.0000000    -0.2500000
     0.0000000    -0.2500000     0.2500000
     0.2500000     0.0000000    -0.2500000
    -0.2500000    -0.2500000    -0.5000000
     0.2500000     0.5000000     0.2500000
     0.2500000    -0.2500000     0.0000000
    -0.5000000    -0.2500000    -0.2500000
     0.2500000     0.2500000     0.5000000
     0.0000000     0.2500000    -0.2500000
    -0.2500000    -0.5000000    -0.2500000
     0.5000000     0.2500000     0.2500000
    -0.2500000     0.0000000     0.2500000
     0.0000000     0.5000000     0.5000000
     0.5000000     0.0000000     0.5000000
     0.2500000    -0.2500000     0.5000000
     0.5000000     0.2500000    -0.2500000
    -0.5000000    -0.2500000    -0.7500000
     0.2500000    -0.5000000    -0.2500000
    -0.2500000     0.2500000    -0.5000000
end kpoints

Compute Wannier functions
run wannier90:

wannier90.x wannier90

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:

./gap_GW.sh OUTCAR

1.5 提高效率与准确性

发表于 2023-05-09 更新于 2025-04-03 分类于 DFT学习

截断能

截断能设置在INCAR文件的ENCUT参数。截断能（energy cutoff）指定了展开波函数用到的平面波基组（planewave basis set）所对应的截断能量。
参考： https://zhuanlan.zhihu.com/p/348826693
截断能越大，代表用来描述波函数的平面波基组越多，计算精度越高，同时计算耗时也越长。
如果没有指定ENCUT，程序将默认ENCUT=ENMAX（注意：当PREC = High时, ENCUT = 1.3*ENMAX）。扩胞并不会影响ENCUT的设置。

测试标准
ENMAX在保证定性正确的基础上大幅度降低计算耗时，如果想要定量计算某些性质就需要增大截断能了。所以ENCUT测试通常从ENMAX开始测试，一般截止到 1.5*ENMAX 即可（注意：比ENMAX小的测试可能会得到错误的能量）。当体系中每个原子的能量差收敛至0.001 eV/atom时，此ENCUT即为最终取值。

简单测试：
for i in {1..8}; do cp 400 $((400+$i*50)); sed -i "s/400/$((400+$i*50))/g" $((400+$i*50))/INCAR ; done

encut.sh脚本： 参考 https://zhuanlan.zhihu.com/p/115235135

{.line-numbers}

#!/bin/bash

#生成KPOINTS文件，KPOINTS不能与INCAR中的KSPACING参数共用，二选一
cat > KPOINTS <<!
A
0
M
9  9  9
0  0  0
!

#产生计算所需POSCAR
cat > POSCAR<<!
Si8
1.0000000000
    5.4687280000        0.0000000000        0.0000000000
    0.0000000000        5.4687280000        0.0000000000
    0.0000000000        0.0000000000        5.4687280000
Si
8
Direct
    0.0000000000        0.0000000000        0.0000000000
    0.2500000000        0.7500000000        0.7500000000
    0.5000000000        0.0000000000        0.5000000000
    0.0000000000        0.5000000000        0.5000000000
    0.5000000000        0.5000000000        0.0000000000
    0.7500000000        0.2500000000        0.7500000000
    0.7500000000        0.7500000000        0.2500000000
    0.2500000000        0.2500000000        0.2500000000
!

for i in $(seq 400 50 700) #截断能从450-800，步数为50
do

#生成vasp静态计算的INCAR
cat > INCAR_static <<!
Global Parameters
  ISTART =  0            (Read existing wavefunction; if there)
  ICHARG =  2            (Non-self-consistent: GGA/LDA band structures)
  LREAL  = .FALSE.       (Projection operators: automatic)
  ENCUT  =  $i           (Cut-off energy for plane wave basis set, in eV)
  PREC   =  Accurate     (Precision level)  
  LWAVE  = .FALSE.       (Write WAVECAR or not)
  LCHARG = .FALSE.       (Write CHGCAR or not) 
  ADDGRID= .TRUE.        (Increase grid; helps GGA convergence)

Electronic Relaxation
  ISMEAR =  -5           (Gaussian smearing; metals:1)
  #SIGMA  =  0.05        (Smearing value in eV; metals:0.2)
  NELM   =  60           (Max electronic SCF steps)  
  NELMIN =  4            (Min electronic SCF steps)
  EDIFF  =  1E-06        (SCF energy convergence; in eV)
  GGA  =  PE             (PBEsol exchange-correlation)

Ionic Relaxation
  ISIF   =  2            (Stress/relaxation: 2-Ions, 3-Shape/Ions/V，4-Shape/Ions)
  EDIFFG = -0.001        (Ionic convergence; eV/AA)
  #KSPACING = 0.10
!

cp INCAR_static INCAR
echo "ENCUT = $i eV" ; time mpirun -np 16 vasp_std #vasp并行运行命令，根据系统自行修改

rm INCAR_static
rm INCAR
#提取计算得到的能量
E=$(grep "TOTEN" OUTCAR | tail -1 | awk '{printf "%12.9f \n", $5 }')
echo $i $E >>encut_energy.out
done

K点测试

是倒空间（动量空间）的基本构成点，只取在一个倒空间晶格向量的范围内来描述 k 。总能量计算是对布里渊区内的波函数进行积分，而积分是通过对部分特殊K点的求和完成的。因此，计算一个体系时，需要确定所用KPOINTS的大小。
k.sh如下：

{.line-numbers}

#!/bin/bash
#生成vasp静态计算的INCAR
cat > INCAR_static <<!
Global Parameters
  ISTART =  0            (Read existing wavefunction; if there)
  ICHARG =  2            (Non-self-consistent: GGA/LDA band structures)
  LREAL  = .FALSE.       (Projection operators: automatic)
  ENCUT  =  500           (Cut-off energy for plane wave basis set, in eV)
  PREC   =  Accurate     (Precision level)  
  LWAVE  = .FALSE.       (Write WAVECAR or not)
  LCHARG = .FALSE.       (Write CHGCAR or not) 
  ADDGRID= .TRUE.        (Increase grid; helps GGA convergence)

Electronic Relaxation
  ISMEAR =  -5           (Gaussian smearing; metals:1)
  #SIGMA  =  0.05        (Smearing value in eV; metals:0.2)
  NELM   =  60           (Max electronic SCF steps)  
  NELMIN =  4            (Min electronic SCF steps)
  EDIFF  =  1E-06        (SCF energy convergence; in eV)
  GGA  =  PE             (PBEsol exchange-correlation)

Ionic Relaxation
  ISIF   =  2            (Stress/relaxation: 2-Ions, 3-Shape/Ions/V，4-Shape/Ions)
  EDIFFG = -0.001        (Ionic convergence; eV/AA)
  #KSPACING = 0.10
!

cp INCAR_static INCAR

#产生计算所需POSCAR
cat > POSCAR<<!
Si8
1.0000000000
    5.4687280000        0.0000000000        0.0000000000
    0.0000000000        5.4687280000        0.0000000000
    0.0000000000        0.0000000000        5.4687280000
Si
8
Direct
    0.0000000000        0.0000000000        0.0000000000
    0.2500000000        0.7500000000        0.7500000000
    0.5000000000        0.0000000000        0.5000000000
    0.0000000000        0.5000000000        0.5000000000
    0.5000000000        0.5000000000        0.0000000000
    0.7500000000        0.2500000000        0.7500000000
    0.7500000000        0.7500000000        0.2500000000
    0.2500000000        0.2500000000        0.2500000000
!


for i in $(seq 6 3 15) #K点从6-24，步数为3
do

#生成KPOINTS文件，K点类型：Monkhorst-Pack (M) 或者 Gamma (G)
cat > KPOINTS <<!
Automatic mesh
0
Monkhorst-Pack
$i $i $i
0  0  0
!

echo "KPOINTS = $i" ; time mpirun -n 16 vasp_std #vasp并行运行命令，根据系统自行修改

rm KPOINTS
#提取计算得到的能量
E=$(grep "TOTEN" OUTCAR | tail -1 | awk '{printf "%12.9f \n", $5 }')
echo $i $E >>kpoints_energy.out
done

显示收敛情况的脚本

参考：https://blog.shishiruqi.com/2019/08/19/paili-energyforce/

简单版：

grep F= OSZICAR |awk '{print $1,$5}'
grep F= OSZICAR |awk '{print $1,$5}' > conv.dat

小脚本：

{.line-numbers}

#!/usr/bin/env bash
for i in *; do
if [ -e $i/OUTCAR ]; then
echo -e  $i "\t" $(grep '  without' $i/OUTCAR |tail -n 1 | awk '{print $7}')
fi
done

复杂：
ediff.sh如下：（安装gnuplot）

{.line-numbers}

#!/bin/sh
#lipai@mail.ustc.edu.cn

begin=$1
if [ $begin =='' ]; then begin=0; fi    
awk -v begin=$begin '/E0/{if ( i<begin  ) i++;else print $0 }' OSZICAR >temp.e
    
awk '/POSITION/,/drift/{
    if(NF==6) print $4,$5,$6;
    else if($1=="total") print $1 }' OUTCAR >temp.f

awk '{if($4=="F"||$4=="T") print $4,$5,$6}' CONTCAR >temp.fix
flag=`wc temp.fix|awk '{print $1}'`
steps=`grep E0 OSZICAR |tail -1 |awk '{print $1}'`
if [ flag != '0' ] ; then
    if [ -f temp.fixx ] ; then rm temp.fixx ; fi
    for i in `seq $steps`;do
        cat temp.fix >>temp.fixx
        echo >>temp.fixx
    done
    paste temp.f temp.fixx >temp.ff
fi

awk  '{ if($1=="total") {print ++i,a;a=0}
        else {
            if($4=="F") x=0; else x=$1;
            if($5=="F") y=0; else y=$2;
            if($6=="F") z=0; else z=$3;
            force=sqrt(x^2+y^2+z^2);
            if(a<force) a=force} }' temp.ff >force.conv

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}

#!/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

```sed '3s/0.01/0.02/g' INCAR > INCAR_new```  #将INCAR中的第三行的0.01全部替换成0.02并输出位INCAR_new文件。
```sed -i '3s/0.01/0.02/g' INCAR``` # -i表示直接对源文件进行编辑，也就是说编辑之后源文件被新文件替换掉。
```for i in *; do sed -i "3s/0.05/$i/g" $i/INCAR ; done``` #将文件夹中所有的0.01替换成了文件夹序号。

我们主要对KPOINTS的文件的第四行进行批量操作，将1 1 1改成 2 2 2， 3 3 3 等。复制文件夹：
```for i in {1..6}; do cp ../ex03/0.01/ ${i}${i}${i} -r ; done```
显示行号：
```cat 333/KPOINTS  -n ```
修改文件;
```for i in {1..6}; do sed -i '4s/1 1 1/$i $i $i/g' ${i}${i}${i}/KPOINTS ; done```

脚本：
 ```javascript  {.line-numbers} 
for i in 3RD.POSCAR.*
do
s=$(echo $i| cut -d "." -f 3) #代表.的3个字符
d=job-$s
mkdir $d
cp $i $d/POSCAR
cp ./INCAR ./sub.sh ./KPOINTS ./POTCAR $d
pushd $d
popd
done

节约机时

grep User OUTCAR 后得出的结果被空格分成了4部分，时间信息在第4部分里面。User 或者 Elapsed都可以。

一个合理的初始结构，可以避免很多意外的错误以及快速得到正确的结果。
如果你设置的收敛标准比较严格，精度比较高，则每一个离子步需要更多的电子步数，需要的时间也会随着离子步的增加成线性关系增长。
ENCUT和KPOINTS文件需要测试，二者越大时间越久。
这里我们就会需要一些其他的参数，NCORE和NPAR。
NCORE：控制多少个核同时计算；
NPAR：如何把计算任务分配到计算资源上面计算。
它们之间的关系是：NCORE= 计算使用的核数 / NPAR。注意：这两个参数只能选取一个来使用:
NCORE=单节点的核数，单节点的核数/2，单节点的核数/4…….个人使用经验是：NCORE = 单个节点核数 / 2 的时候，运行最省时间，设置也最方便。
当然NCORE也需要测试。

消虚频

结构优化消虚频

X, Y, Z为体系中原子的坐标，也就是POSCAR中的内容。 dx, dy,dz为虚频对应的原子振动方向上的位移。我们将（dx, dy,dz）这个位移乘以一个介于0-1之间的校正因子，然后跟POSCAR中的坐标加起来即可。以z方向为例，校正因子设为0.1， 微调后的z方向坐标为： 10.709980 + （-0.676634）0.1。
脚本vib_correct.py如下：https://www.bigbrosci.com/2022/04/29/A29/

{.line-numbers}

# -*- coding: utf-8 -*-
"""
Qiang Li
Command: python3 vib_correct.py 
This is a temporary script file.
"""
import numpy as np
from ase.io import read, write
import os

#获取虚频开始行
l_position = 0  #虚频振动方向部分在OUTCAR中的起始行数
with open('OUTCAR') as f_in:
    lines = f_in.readlines()
    wave_num = 0.0
    for num, line in enumerate(lines):
        if 'f/i' in line:
            wave_tem = float(line.rstrip().split()[6])
            if wave_tem > wave_num: #获取最大的虚频
                wave_num = wave_tem
                l_position = num+2
# ASE读POSCAR
model = read('POSCAR')
model_positions = model.get_positions()
num_atoms = len(model)
#print(model_positions)

# 获取虚频对应的OUTCAR部分 
vib_lines = lines[l_position:l_position + num_atoms] #振动部分 7222到7248行
#print(vib_lines)
vib_dis = []
for line in vib_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中。

# 微调结构
new_positions = model_positions + vib_dis * 0.4 # 0.4是微调的校正因子，即虚频对应振动位移的0.4，具体多大自己根据经验调。
model.positions = new_positions

###保存结构
write('POSCAR_new', model, vasp5=True) # POSCAR_new是微调后的结构，用于下一步的计算（别忘了把POSCAR_new改成POSCAR）。

过渡态虚频

如果过渡态计算的时候出现好几个虚频，自己直接写个python脚本也能实现。不过使用ASE可以很方便地获取POSCAR中原子数目，保存原子固定的信息。参考如下。

{.line-numbers}

#!/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

#获取虚频开始行
l_position = 0  #虚频振动方向部分在OUTCAR中的起始行数
with open('OUTCAR') as f_in:
    lines = f_in.readlines()
    wave_num = 0.0
    for num, line in enumerate(lines):
        if 'f/i' in line:
            wave_tem = float(line.rstrip().split()[6])
            if wave_tem > wave_num: #获取最大的虚频
                wave_num = wave_tem
                l_position = num+2
# 读POSCAR
with open('POSCAR') as f_pos:
    lines_pos = f_pos.readlines()

# 获取虚频对应的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 in vib_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)

# 微调结构
model_positions = np.array(model_positions)
vib_dis = np.array(vib_dis) #将振动部分写到一个array中。
new_positions = model_positions + vib_dis * 0.4 # 0.4是微调的校正因子，即虚频对应振动位移的0.4，具体多大自己根据经验调。

###保存结构
f_out = open('POSCAR_new','w')
f_out.writelines(lines_pos[:8])
f_out.write('Cartesian\n')
for i in new_positions:
    f_out.write(' '.join([str(coord) for coord in i]) + '  F  F  F\n')
f_out.close()

滑移堆叠能计算

vaspkit 926

脚本

该脚本是计算双层的堆叠能，滑移能类似。准备优化好的POSCAR和其他输入文件。
stack.sh如下：

#!/bin/bash

# 参数设置
slip_step_x=$(echo "scale=5; 1/12" | bc)  # x方向滑移步长 (1/12) scale=5 设置浮点数精度为 5 位小数
slip_step_y=$(echo "scale=5; 1/12" | bc)  # y方向滑移步长 (1/12)
slip_steps=11                             # 滑移步数（从0到11，共12步）

# 循环生成滑移结构
for i in $(seq 0 $slip_steps); do
    for j in $(seq 0 $slip_steps); do
        # 格式化文件夹名称（固定宽度，两位数字）
        folder_name=$(printf "slip_%02d_%02d" "$i" "$j")
        mkdir -p "$folder_name"
        
        # 复制输入文件（排除目录）
        find . -maxdepth 1 -type f -name "[IKPs]*" -exec cp {} "$folder_name/" \;   #这里排除复制文件夹
        cp ./sub.sh "$folder_name/"   # 复制 sub.sh提交任务的脚本
        find . -maxdepth 1 -type f -name "vdw-*" -exec cp {} "$folder_name/" \;  # 这里使用opt86b范德华修正，自行选择
        
        # 计算滑移量
        slip_x=$(echo "scale=5; $i * $slip_step_x" | bc)
        slip_y=$(echo "scale=5; $j * $slip_step_y" | bc)
        
        # 修改POSCAR文件
        awk -v slip_x="$slip_x" -v slip_y="$slip_y" '
        {
            if (NR > 9 && $3 > 0.5) {
                print $1 + slip_x, $2 + slip_y, $3, $4, $5, $6
            } else {
                print $0
            }
        }' POSCAR > "$folder_name/POSCAR"
    done
done

echo "滑移结构生成完成。"

解释：由于是六方结构，所有滑移步长1/12,可自行更改。if (NR > 9 && $3 > 0.5)里的0.5在两层之间，大于0.5的移动，小于等于不移动。要固定晶格 ISIF = 2，离子要优化，要固定一些原子的xy方向的坐标，防止平面间的滑移，优化z方向。

批量提交任务：for dir in slip_*/; do cd “$dir”；qsub sub.sh； cd .. ；done

取消所有任务(slurm)：scancel -u $USER

画图：

import numpy as np
import matplotlib.pyplot as plt

def load_data(filename):
    """加载数据文件并返回 numpy 数组"""
    try:
        data = np.loadtxt(filename)
        return data
    except Exception as e:
        print(f"Error loading file {filename}: {e}")
        return None

def reshape_data(x, y, z, shape=(11, 11)):
    """将一维数组重塑为二维数组"""
    x2 = x.reshape(shape)
    y2 = y.reshape(shape)
    z2 = z.reshape(shape)
    return x2, y2, z2

def calculate_coordinates(x, y):
    """计算新的坐标 X 和 Y"""
    X = x - 0.5 * y
    Y = 0.5 * np.sqrt(3) * y
    return X, Y

def plot_contour(X, Y, Z, title, filename, levels=100, cmap='viridis'):
    """绘制等高线图并保存"""
    fig, ax = plt.subplots(figsize=(6, 6))
    q = ax.contourf(X, Y, Z, levels=levels, cmap=cmap)
    fig.colorbar(q, label='Energy (eV)')
    ax.set_xlabel('X')
    ax.set_ylabel('Y')
    ax.set_title(title)
    ax.set_aspect(1)
    plt.savefig(filename, dpi=300, bbox_inches='tight')
    plt.close()

def main():
    # 加载数据
    fm = load_data('fm_Energy.dat')
    afm = load_data('afm_Energy.dat')

    if fm is None or afm is None:
        return

    # 提取数据
    x = fm[:, 0]
    y = fm[:, 1]
    zf = fm[:, 2]
    za = afm[:, 2]

    # 重塑数据
    x2, y2, Zf = reshape_data(x, y, zf)
    _, _, Za = reshape_data(x, y, za)

    # 计算新坐标
    X, Y = calculate_coordinates(x2, y2)

    # 计算能量差
    Zd = Zf - Za

    # 绘制图形
    plot_contour(X, Y, Zf, 'Ferromagnetic (FM) Energy', 'FM.png')
    plot_contour(X, Y, Za, 'Antiferromagnetic (AFM) Energy', 'AFM.png')
    plot_contour(X, Y, Zd, 'FM - AFM Energy Difference', 'FM-AFM.png')

if __name__ == "__main__":
    main()

常用命令行

grep "33 81" -A3 FORCE_CONSTANTS | xargs echo -n| awk '{print $1,$2,$3+$7+$11}' >> trace.dat

计算原子间二阶力常数矩阵的迹

grep 'F=' OSZICAR|awk '{print $1, "\t"$3,"\t"$9}' >AIMD.dat

提取分子动力学温度能量数值

for i in 01 02 03 04 05;do cd $i && cp CONTCAR POSCAR && cd ..;done

CINEB计算重新计算

awk '{print $1,$2,$3+$4+$5,$6+$7+$8+$9+$10,$11}' PDOS > pdos.dat

pdos轨道加和

grep LOOP+ OUTCAR |awk '{sum+=$4} END {print "Average=",sum/NR}'

grep LOOP+ OUTCAR |awk '{sum+=$4} END {print "Sum=",sum}'

grep LOOP+ OUTCAR |wc -l

提取平均离子步时间，总时间，总步数

echo -e "4\n403\n49-96 121-144\n1\n"|vaspkit

vaspkit命令行

for i in *; do mkdir $i && cp INCAR KPOINTS POSCAR POTCAR sub.sh $i && sed -i "3s/0.05/$i/g" $i/INCAR ; done

批量修改INCAR

sed -i '5s/\<0\>/0.5/' filename
请将 filename 替换为您要修改的文件名。
5s/<0>/0.5/：
5：指定要修改的行号（第 5 行）。
s：表示替换操作。
<0>：匹配完整的数字 0，确保只替换单独的 0。
0.5：要替换成的新值。

1.4 过渡态

发表于 2023-05-09 更新于 2024-06-05 分类于 DFT学习

插入的image只需要看一下有没有原子互相挨得很近，有的话调一下就行。
一开始计算先用IOPT=7，比IOPT=3快一点，然后LCLIMB=.FALSE.;算个一二十步之后再换IOPT=1或2，LCLIMB=.TRUE.
NEB计算不需要收敛，从经验来看也极少能碰到收敛的情况。主要是找到切线力尽可能等于0的点，然后用dimer算，基本都能找到过渡态。
nebresult.pl没用过，看结果“forces and energy:”下面的数据的第四列应该就是切线力，正常的计算切线力是由正值变到负值，或者由负值变到正值，然后找接近于0的点。你图里的数据只有正值，因此要么初末态构型有问题，要么你要找的点在0和1两个image之间。
不过1号image的切线力已经很接近0了，可以尝试直接做过dimer
跑个一二十步之后停止计算，‘mv CONTCAR POSCAR’之后修改INCAR再继续算。脚本都是线性插值，因此初始构型一般都很不合理，上来就用IOPT=2的话收敛慢，甚至可能出错。跑个几十步等构型相对合理一些之后，再换IOPT=2。
dimer是vasp里面自带的算法，详细用法手册里有。

计算步骤

初末态结构优化
建立两个文件夹ini(初态的意思),fin(末态)，每个文件夹放入vasp计算必备的四个文件（INCAR，POSCAR，KPOINTS，POTCAR），其中的两个POSCAR对应未优化的初末态。确保两个文件夹里面除POSCAR外，其他文件完全一样。

对于两个POSCAR，注意每一行的原子一一对应。若有固定位置的原子（比方做表面计算，要固定底部1-2层原子），最好让两个POSCAR里面的固定原子位置完全相同。
用dist.pl检查两个优化后结构(两个CONTCAR)的相似程度（每个对应原子的初末态距离的平方和，再开根号）。

具体命令为： dist.pl ini/CONTCAR fin/CONTCAR
若返回值小于5A，一般可以进行下一步。
2. 插点
插点的数目取决于前面dist.pl的返回值，一般插点数目可取（dist.pl返回值/0.8）。

具体插点的命令： nebmake.pl ini/CONTCAR fin/CONTCAR N

最后的N表示插点数目。插入点的算法为线性插值，详情请进前面给的vtst脚本链接。我这里设置N为3，执行完命令后出现文件夹00，01，02，03，04.五个文件夹里面都有nebmake.pl产生的POSCAR。

00表示初态，里面放的是ini/CONTCAR,04表示末态，放的是fin/CONTCAR, 01、02、03是插入的点。

输出中提示讲初末态对应的OUTCAR复制到对应的文件夹中，以便后续数据分析。
检查插入点的合理性：

输入命令 nebmovie.pl 0

参数0表示用POSCAR生成xyz文件；还可取1，为用CONTCAR生成。

可以看到每个文件夹内多了一个后缀名为xyz的文件。将这几个xyz文件合成一个movie.xyz文件即可在VMD中动画演示。我更倾向于使用jmol(或者ovito)来查看xyz文件，不过VMD可能大家更熟悉。

用nebmovie.pl,没有直接生成movie.xyz是因为从官方主页下载的脚本默认当使用nebmovie.pl后面不带参数或者参数为0的时候(即使用POSCAR产生xyz)，不输出movie.xyz。你想让脚本自动生成movie.xyz而不是自己去合并各个文件夹里面的xyz文件，需要自行修改nebmovie.pl文件。很简单，把倒数第二个if语句整个用#注释掉或者直接删掉即可。

也可以直接用VESTA查看01的POSCAR。

进一步可用命令nebavoid.pl 1,确保中间插入的点每一个原子间距都大于1A。该命令的参数1表示最小允许间距，可取小数。
原子间距太小说明结构有问题。要么是因为初末态POSCAR里的原子没有对应好，要么是因为nebmake.pl线性插值的方法不适用于估计反应路径路径。如果是前者，调整好对应原子，也就是CONTCAR里面相应行进行调换；如果是后者，用Materials studio之类的软件调整好结构，再放入对应文件夹。
3. 其它文件的准备
在当前目录下面放入KPOINTS，POTCAR及INCAR文件。要求KPOINTS和POTCAR与ini、fin文件夹中的一样；INCAR中的基本参数也与初末态计算的INCAR一样，另外再加入NEB计算的相关参数。
INCAR再多说两句。优化算法若想使用vtst的，需要如上设置IBRION=3,POTIM=0,IOPT设置见本文第四个图。若使用vasp自带的优化器，可以使用IBRION=1或者3，不要取2；POTIM取合适的值（0.01-0.5范围内去尝试）。

这样就可以提交NEB任务了。
4. 后处理
计算过程中，可以用nebef.pl命令（不带参数）查看计算收敛情况。
输出中，第二列即为最大受力（force of images in the neb），第三列为相应结构的能量。

前面INCAR中EDIFFG设置为-0.02（一般结构优化可取-0.02，NEB计算取-0.03），表示当所有结构最大力小于0.02eV/A时，结束计算。

另一个可以用来观察收敛情况的命令是nebbarrier.pl（同样不带参数）。该命令没有输出到屏幕的内容，而是生成neb.dat文件。
neb.dat文件第二列表示距离（即临近两结构的dist.pl的计算结果），第三列表示能量（以初态能量为参考值），第四列为力（forces along the neb）。

EDIFFG参数对应的力是nebef.pl输出中的force of images in the neb.

计算完成后，使用命令nebresult.pl.

nebresult.pl做的事情如其所输出表明的，完成了nebbarrier.pl, nebspline.pl, nebef.pl, nebmovie.pl, nebjmovie.pl, nebconverge.pl还有对各文件夹中的OUTCAR打包压缩。生成了很多文件。其中mep.eps是以dist.pl距离为横坐标，能量为纵坐标画出的能势垒图，如下：
该文件一般查看图片的软件打不开，推荐使用EPS/PS viewer. https://epsviewer.org/download.aspx

生成的spline.dat文件是对上面几个点的拟合曲线数据，已经有了上图，这些数据就意义不大了。

生成的vaspgr文件夹内是各个插点结构的收敛图（同样可用EPS/PS viewer打开），如下vaspout1.eps是01结构收敛信息：

1.3 结构优化3

发表于 2023-05-09 更新于 2025-07-08 分类于 DFT学习

加压

可通过施加静水压强、调整晶格常数的方法实现结构压缩

vi INCAR 增加标签 PSTRESS=50
（单位为 0.1 GPa）

其他不变

加电场

在vasp中根据自己所需要添加电场的场强更改INCAR中的如下参数：

1
2
3

LDIPOL=.TRUE.
IDIPOL=1/2/3/4
EFIELD= ？

对于IDIPOL参数的设置，需注意：对于IDIPOL=1-3，偶极矩将分别平行于第一、第二或第三晶格矢量的方向计算。总能量的修正计算为当前超单元中单极子/偶极子和四极子与放置在超单元中的同一偶极子之间的能量差，对应的晶格矢量接近无穷大。此标志应用于板计算。对于IDIPOL=4，将计算所有方向上的完整偶极矩，并将总能量的修正计算为当前超级单体中的单极/偶极/四极和放置在真空中的同一单极/偶极/四极之间的能量差，使用该标志计算孤立分子。

此外，EFIELD 的设置需要注意单位的换算。eV/Å（电子伏特/埃）计算完毕之后通过后处理软件分析所关注的性质即可。

状态方程拟合

我们需要对结构进行优化才可以获取稳定的晶格参数信息。有两个方法可以实现：

1 Birch-Murnaghan状态方程拟合,
2 VASP计算中通过调节ISIF参数直接优化Bulk。
直接优化前面已经介绍，下面我们先讨论一下第一个方法：BM方程拟合。

参考：大师兄科研网：https://www.bigbrosci.com/2018/02/02/ex33/

科学网耿华运的博文：https://blog.sciencenet.cn/blog-3469498-1280203.html

科学网赵志强的博文，关于拟合：https://wap.sciencenet.cn/blog-3388193-1152915.html

vaspkit：http://vaspkit.cn/index.php/48.html

大致流程（以大师兄提供脚本为例）

获得一定缩放系数的结构：for i in 0.95 0.96 0.97 0.98 0.99 1.01 1.02 1.03 1.04 1.05; do cp 1.00 $i ; sed -i "2s/1.0/$i/g" $i/POSCAR ; done
分别进行结构优化:for i in *; do cd $i; 提交任务命令; cd $OLDPWD; done
提取能量-体积以及能量-压力、压力-体积的数据:for i in *; do echo -e $i "\t" $(grep ' without' $i/OUTCAR | tail -n 1| awk '{print $7}'); done > data
拟合能量体积曲线，获得平衡晶胞参数

这个流程可通过批量计算脚本实现，参考上面链接。另外还可以获得能量随压力的变化来模拟相变。

分数占据

生成disorder或者掺杂构型的软件汇总：

https://mp.weixin.qq.com/s/vYar9UcYWAZgazBTgqUGEA
https://github.com/jichunlian/disorder
supercell： https://orex.github.io/supercell/
ATAT： https://www.brown.edu/Departments/Engineering/Labs/avdw/atat/
vaspkit：https://vaspkit.com/

disorder

安装

$ 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

扩建超胞，保证整除

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文件

&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.
/

官方解释：

nsub-tag:nsub = 2 ~ 5; Default: nsub = 2.该位置取代原子种类。本例中是2种。
subs-tag：subs = integer, integer, …; Default: None。对 应的原子数目。总和加起来应该等于该位点在SPOSCAR中的数目。
symb-tag：symb = character, character, …; Default: None。元素化学符号。空位用Kw,即汉语kong wei。本例是Sn Pb单引号。
site-tag：site = integer; Default: 1。所取代的位置在SPOSCAR中的位置。本例是第2.
prec-tag：prec = real; Default: 1D-5。当软件产生空间群与其他软件识别不一致时调整该参数，一般不用修改。
lpro-tag：lpro = logical; Default: .false.是否消除重复参数。对于大体系可以打开。
lpos、leqa、lspg、lcfg：分别控制输出文件POSCAR-X、EQAMAT、SPGMAT、CFGMAT。

运行disorder。输入disorder即可。
输出文件：

ATAT软件

简介

什么是ATAT？

其全称为 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.

如何安装ATAT：

http://brown.edu/Departments/Engineering/Labs/avdw//atat/atat3_36.tar.gz下载ATAT源码，解压后，make all开始编译；之后make install进行程序和脚本定向安装；之后在home中打开./bashrc添加ATAT安装路径于环境变量当中（添加 export PATH=$PATH:/xxxx），并source ./bashrc即可。

进入安装包的glue目录，里面有与各种dft软件的算例的接口设置。例如进入vasp可以看到它的接口设置脚本ezvasp，由于make install已经将它拷贝到了安装目录中，所以在bin中打开ezvasp，并设置相关可执行程序，包括赝势文件的名称和路径即可。

配置~/.ezvasp.rc

1
2
3

#!/bin/csh
set VASPCMD="/path/vasp_std# 设置vasp可执行程序
set POTLDA/PBE="/path/to/pot_lad" # 设置为赝势实际的路径

ATAT能做啥？

ATAT 是Walle及其合作者开发的一个利用第一性原理结合MonteCarlo模拟计算相图及热力学性能的软件包。该软件包首先采用团簇展开(cluster expansion)，用第一原理方法计算一系列构型的形成能，然后用这些能量作为Monte Carlo模拟的哈密顿量，计算材料的热力学性能和高熵合金的mcsqs建模。

MAPs

参考：
https://github.com/xu-xi/clusterexpansion/blob/master/atat.md

软件运行

一）运行MAPS：

准备lat.in文件

2）maps -d & #maps命令运行，等待后面的指令

3）touch ready

#等待十几秒钟，出现“0”文件夹，是合金中比例为0的纯相结构。

4）此时需要利用第一性原理软件计算该结构的能量

#（参见下部分和VASP软件的交互,可无需手动输入），现在手动输入能量，理解一下计算过程，假设该结构的能量已经知道，在0这个文件夹下，echo 1.1 > energy, 并删除wait。继续touch ready，就会生成1这个文件夹，里面是比例为1的另一个纯相结构，echo 2.2 > energy, 删除wait。touch ready，生成2这个文件夹，一直这样循环下去，直至达到相应的精度，最终停止输入touch stop。

二）运行MAPS与VASP交互：

1）需要lat.in，vasp.wrap文件（ #vasp的INCAR参数）

2）其余与上面运行MAPS步骤相同，直到第4步骤计算能量，进入该文件夹，运行runstruct_vasp。

大概率会报错，原因有几个：

a. 常见为~/.ezvasp.rc文件配置问题，赝势位置，vasp运行路径出错。

可以看运行命令后，是否有POTCAR文件产生，如果没有，基本就是赝势位置问题，如果POTCAR文件正常，就是vasp没有运行。

b. vasp不运行，在集群上面运行的软件，ATAT软件给了集群上面的运行办法, pollmach命令，将runstruct_vasp命令写在脚本里是可以的。

最终正确版本运行步骤：

1）准备lat.in，vasp.wrap文件，脚本文件。

lat.in文件

vasp.wrap文件

sub.sh文件

2）提交脚本文件即可

说明：

缺少.machines.rc文件的提示可忽略，该文件用于跨节点的并行
maps不断产生结构和构建团簇展开模型，-d表示使用默认参数，直接输入maps可查看可用参数
runstruct_vasp调用vasp计算所产生的结构的能量，每个目录包含一个结构及其计算结果，每个目录下的energy文件即从vasp的计算中读取的能量
maps.log中记载了现有的团簇展开模型的状态，包括对基态构型的预测和LOOCV等

什么时候停止计算：

团簇展开模型能够成功预测基态 (查看maps.log)
LOOCV足够小 (<0.025 eV suggested by ATAT manual)
有效团簇相互作用随着团簇距离的增大而衰减

以上的三个判据并不是绝对的，判断一个团簇展开模型是合理的还可以有其它方式。

停止计算：

1 2	$ touch stop $ touch stoppoll

重启计算: 与正常开始一个计算相同，如果之前的计算不是正常结束则需要删除一些文件

1
2
3

rm maps_is_running pollmach_is_running ready
maps -d &
pollmach runstruct_vasp

常用命令

corrdump -2= -3= -4= 产生团簇
getclus 查看团簇
- 第一列为团簇的阶数，第二列为团簇直径，第三列为团簇的多重度
- clusters.out中有团簇格点位置的详细信息
clusterexpand -e -cv energy 得到团簇展开系数
- 此时产生的ECI文件为energy.eci，而maps产生的ECI为eci.out，后者是形成能的团簇展开结果，前者是直接对vasp能量(即energy文件)的团簇展开
- 如果需要对其它性质做团簇展开，例如带隙，则可以在每个结构的目录下创建一个bandgap文件,其内容即为带隙的大小，然后运行clusterexpand -e -cv -pa bandgap以得到带隙团簇展开的系数bandgap.eci，其中-pa表示所展开的量已经是一个平均的量

检查计算结果

VASP是否正常结束
- ATAT会自动检查一些错误(checkerr_vasp)
- 但如果结构优化的离子步数不够大，不会有提示，因此建议在vasp.wrap中将NSW设置大一些
晶胞的应变是否很大
- checkrelax (<0.1 suggested by ATAT）

结果可视化

mapsrep 会调用gnuplot作图并产生相应的脚本mapsrep.gnu

输出文件

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
gs_str.out：基态的结构，格式和str.out一样。每一个结构后面会有end和一个空行隔开。
eci.out：eci值对应的cluster在clusters.out文件。
clusters.out：对于每一个团簇，第一行是多重度 multiplicity，第二行是团簇维度cluster diameter，第三行是团簇的格点数，下一行是格点对应的坐标。最后空行隔开每一个团簇。
str.out：在每个数字文件夹下面，和晶格输入文件相似，不同是坐标写成3*3矩阵，一个位置只写一个原子。
ref_energy.out：用于计算形成能的参考能量。

后处理

一系列用于分析ATAT计算结果的python脚本，GitHub地址为：https://github.com/xu-xi/clusterexpansion可以通过 git clone 下载到本地服务器。用“-h”命令可以查看每个脚本的帮助文档。最常用的几个脚本是：

1）cv_scan.py，测试CE模型对团簇截断半径的收敛性。

2）lc.py，测试CE模型对训练集大小的收敛性。通常~80个构型就可以得到比较精确的团簇展开模型。

3）ce_plot.py，这个脚本绘制的是拟合能量-计算能量分布图，默认-p=energy（调用的是每个数字目录下的energy文件）。energy文件是ATAT完成一个构型的DFT (relax+static) 计算后自动生成的能量文件。在实际应用中，拟合能量-计算能量分布图需要的往往不是energy而是formation energy，所以就需要编写脚本手动在每一个数字目录下生成一个形成能文件formation_energy，该值需要将目标目录下energy的值与0文件夹的energy和1文件夹的energy作线性变换得到，再写入到每一个数字目录中，读者可自行尝试。下面的脚本仅供参考，

#!/usr/bin/env python
import os,subprocess,numpy,argparse,sys,subprocess,math
energy_0 = float(subprocess.check_output('cat ./0/energy', shell = True))
energy_1 = float(subprocess.check_output('cat ./1/energy', shell = True))
for item in os.listdir(os.getcwd()):
    fullpath = os.path.join(os.getcwd(),item)
    if os.path.isdir(fullpath):
        os.chdir(fullpath)
        if not os.path.isfile('error'):
            sc_atom_number = len(file('str.out').readlines())-6
            # print sc_atom_number
            energy = float(subprocess.check_output('cat energy', shell = True))
            energy /= sc_atom_number
            #energy *= 5 normalize to unit cell,5 atoms in CsPbI3
            sc_x = float(subprocess.check_output('corrdump -pc -l=../lat.in',shell=True).strip().split()[0])
            formation_energy = energy - (1 - sc_x) * energy_0 - sc_x * energy_1
            # print formation_energy
            with open('formation energy','w') as f:
                f.write('%.7f\n' % formation_energy)
        os.chdir('../')

注：这里如果报错：NameError: name ‘file’ is not defined. Did you mean: ‘filter’?将file改成open

得到formation_energy文件后，执行命令ce_plot.py -p=formation_energy -s=”formation energy”，还可以用Origin手动绘制出形成能-掺杂浓度分布图；在得到了比较精确的CE模型后，可以进行蒙特卡罗模拟，得到团簇关联函数、系综平均能量和方差随温度变化示意图。

热力学计算

必要的输入文件：

lat.in
clusters.out
eci.out
gs_str.out

`emc2`使用示例

蒙特卡洛模拟

1	emc2 -T0=300 -T1=3000 -dT=100 -cm -x=0 -keV -gs=-1 -er=30 -aq=0 -dx=1e-5 -sigdig=12

-cm表示正则系综
-x其实是点团簇函数，若c表示格点的占据浓度，则x=2c-1
-sigdig=12是为了表示有足够的位数输出能量的方差以判断相变
-aq和-dx是自动确定MC平衡步数和取样步数的算法，其原理可参考文献Ref.1，也可以用-eq和-n来直接指定步数，需要做收敛性测试
主要输出文件为mc.out，输入emc2 -h可查看每一列所表示的内容
输出文件mcsnapshot.out为最后一个模拟温度下的MC快照。如果想得到某个温度下的snapshot，可以只设置T0而不设置T1和dT。
1
2
3
$ cp mcsnapshot.out str.out
$ runstruct_vasp -nr
$ vesta POSCAR

相图计算

基本原理可参考文献Ref.1，Martin Baker的笔记提供了更多的例子和细节

例子：

1	phb -dT=10 -ltep=1e-3 -er=30 -gs1=0 -gs2=1 -o=phb.out -keV -dx=1e-5

构建SQS/SQoS

SQS/SQoS是对无序体系的代表性结构，其中SQS假定原子完全随机的占据，而SQoS是SQS的扩展，可以考虑一定的短程有序性，其基本原理可参考相关文献Ref.2

必要文件：

rndstr.in
clusters.out

rndstr.in的一个例子：

4.1000 0.000000 0.000000
0.000000 4.1000 0.000000
0.000000 0.000000 4.1000
1.000000 0.000000 0.000000
0.000000 1.000000 0.000000
0.000000 0.000000 1.000000
0.000000 0.500000 0.500000 O=0.6667,N=0.3333
0.500000 0.000000 0.500000 O=0.6667,N=0.3333
0.500000 0.500000 0.000000 O=0.6667,N=0.3333
0.000000 0.000000 0.000000 Ba
0.500000 0.500000 0.500000 Ta

产生团簇文件clusters.out：

1	corrdump -l=rndstr.in -clus -noe -nop -2= -3= -4=

说明：

-noe和-nop表示忽略空团簇和点团簇
用getclus查看所产生的团簇。如果团簇数目过少，则后续产生SQS有可能很容易完全匹配团簇函数，导致SQS实际上不够随机；而如果团簇过多，则有可能导致近邻的团簇匹配得不够好，而通常近邻的团簇是更为重要的，可以通过-wd加上权重以避免这种情况

运行mcsqs -n= &
说明：

-n设定SQS晶胞的原子数
可通过-tcf设定目标团簇函数（可从mc.out中提取），即SQoS方法
mcsqs不会自动停下来，可用touch stopsqs停止
输出文件mcsqs.log中有对团簇函数匹配的优化记录
输出文件bestsqs.out为目前找到的最佳的SQS，可将其复制为str.out，然后运行runstruct_vasp -nr生成vasp的输入文件用于后续计算

MCSQS

一、原理

无序态的关联函数已知，因此只需要备选态的关联函数和无序态关联函数差别最小就可以，即让下方公式中的目标函数最小。

二、运行步骤

step1. 准备rndstr.in文件

1 1 1 90 90 90

1 0 0

0 1 0

0 0 1

.0 .0 .0 Cr

.0 .5 .5 Ni=0.333333,Fe=.333333, Co=.333334

.5 .0 .5 Ni=0.333333,Fe=.333333, Co=.333334

.5 .5 .0 Ni=0.333333,Fe=.333333, Co=.333334

#第一行也可以改成3*3矩阵，类似于POSCAR中晶格常数格式。这里一般不用实际的晶格常数，用1替代。

用1对后面的参数设置更容易。

step2. 运行命令
corrdump -l=rndstr.in -ro -noe -nop -clus -2=1.1
该命令通过读取rndstr.in文件中的结构，进而计算结构的对称性和团簇等信息，输出clusters.out文件。

#需要注意的是‘-2’这个参数，是用于计算关联函数的截断距离。如对于FCC构型而言，晶格常数为1，第一最近距离为2^0.5/2=0.71，第二距离为1，第三距离为1.5^0.5=1.2。所以这里设置阶段距离在1和1.2之间。一般考虑‘-2’参数即可。（注：更高的-2，-3参数会产生更无序的结构，但是需要更长的时间，一般晶格常数为1，这里选取1.1就可以）

step3. 运行命令 mcsqs -n XX
该命令主要输出sqscell.out文件。XX表示体系的原子个数，该数值一定是上述结构原子的倍数，否则会出错。

#比如上述构型，自身CrNi有4个原子，其中Ni元素与Fe,Co构成合金，因此XX至少为12个原子。该步骤运行一下，就可以进行停止。输出文件需要关注sqscell.out文件，该文件给出了XX个原子扩胞后形成的不同晶格构型。可通过挑选，删除一些不太可能的结构，比如x,y,z三个方向长度差别过大，形成的角度较小等对VASP或后续计算不友好的结构。之后用新的sqscell.out文件运行第四步骤。

step4. 运行命令mcsqs -rc
输出bestsqs.out文件和bestcorr.out文件。其中bestsqs.out就是无序结构。

#通常而言没有最完美的无序结构，这个命令会一直运行下去，可以肯定的是，再一次生成的结构比上一次生成的结构要优异，所以什么时候停止取决于自身研究的精度。

层间结合能、结合力

贴一篇文章：https://news.sciencenet.cn/htmlpaper/2022/4/202242914585844872440.shtm

Graphite interlayer distance

vasp上的教程：

准备graphite结构，

INCAR：

IVDW = 20           
LVDW_EWALD =.TRUE. 
NSW = 1 
IBRION = 2
ISIF = 4
PREC = Accurate
EDIFFG = 1e-5
LWAVE = .FALSE.
LCHARG = .FALSE.
ISMEAR = -5
SIGMA = 0.01
EDIFF = 1e-6
ALGO = Fast
NPAR = 2

提交脚本：
To run this example, execute the run.sh bash-script:

#
# 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

for d in 6.5 6.6 6.65 6.7 6.75 6.8 6.9 7.0
do
  cat>POSCAR<<!
graphite
1.0
1.22800000 -2.12695839  0.00000000
1.22800000  2.12695839  0.00000000
0.00000000  0.00000000  $d
4
direct
   0.00000000  0.00000000  0.25000000
   0.00000000  0.00000000  0.75000000
   0.33333333  0.66666667  0.25000000
   0.66666667  0.33333333  0.75000000
!
  $vasp_std
  cp OUTCAR OUTCAR.$d
  energy=$(grep "free  ene" OUTCAR.$d|awk '{print $5}')
  echo $d $energy >> results.dat 
done

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:

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 Å).

Graphite MBD binding energy

IVDW = 202           
LVDWEXPANSION =.TRUE. 
NSW = 1 
IBRION = 2
ISIF = 4
PREC = Accurate
EDIFFG = 1e-5
LWAVE = .FALSE.
LCHARG = .FALSE.
ISMEAR = -5
SIGMA = 0.01
EDIFF = 1e-6
ALGO = Fast
NPAR = 2

KPOINTS

Graphite:

Monkhorst Pack
0
gamma
16 16 8
0 0 0

Graphene:

Monkhorst Pack
0
gamma
16 16 1
0 0 0

Running this example
To run this example, execute the run.sh bash-script:

#
# 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}')

# compute interlayer binding energy (eV/atom)
deltaE=$(echo print $en2/4 - $en1/2 |python)

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].

Graphite TS binding energy

IVDW = 20            
LVDW_EWALD =.TRUE. 
NSW = 1 
IBRION = 2
ISIF = 4
PREC = Accurate
EDIFFG = 1e-5
LWAVE = .FALSE.
LCHARG = .FALSE.
ISMEAR = -5
SIGMA = 0.01
EDIFF = 1e-6
ALGO = Fast
NPAR = 2

To run this example, execute the run.sh bash-script:

#
# 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}')

# compute interlayer binding energy (eV/atom)
deltaE=$(echo print $en2/4 - $en1/2 |python)

echo "Binding energy (eV/atom): " $deltaE >results.dat

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).

1.2 结构优化2

发表于 2023-05-09 更新于 2024-06-21 分类于 DFT学习

结构优化完成后，我们已经可以从输出文件得到很多性质。这里我们先关注能量。

VASP输出的能量是什么

在优化完的文件夹下输入

{.line-numbers}

1 2	[scf1355@ln151%nc-l 2h]$ tail -1 OSZICAR 3 F= -.44829871E+02 E0= -.44829871E+02 d E =-.169526E-03 mag= -0.0000

这里的F和下面free energy含义相同，E0代表体系内能。
另外在OUTCAR中还可以找到这些

{.line-numbers}

1
2
3

free  energy   TOTEN  =       -44.82987077 eV
energy  without entropy=      -44.82987077  energy(sigma->0) =      -44.82987077
atomic energy  EATOM  =      1542.31592482

free energy指的是亥姆霍兹自由能(Helmholtz free energy)，在VASP手册中被称为自由电子能(free electronic energy)，它可以用F表示： F = E-TS+PV .而F和E对应上面的F和E0。
TOTEN和without entropy分别代表考虑和未考虑熵的贡献时自由能的值。而energy(sigma->0)为展宽趋近于零时的自由能，可以外推得到： energy(sigma->0) = E0 =1/2（F+E）

其中上面的展宽对应于INCAR中设置的ISMEAR参数，详情可以查看手册或者参考 https://zhuanlan.zhihu.com/p/384626631。
https://www.bilibili.com/video/BV19L4y1i7mC/?vd_source=17084cc867ff3bb86398ff1a79b443f2
简单理解费米面内电子全占据，费米面外没有电子。这种阶梯函数就使得，随着K点数目变化，K点处积分值变化很大，不易收敛。于是我们需要一个连续函数，代替该阶梯函数进行积分，这种替代方法就叫smear。在VASP计算中，只有在ISMEAR=-5时，才能严格满足费米能级的定义；而在其他情况下，允许电子在费米能级附近的宽度内“波动”，即将阶梯函数替代为一个平滑的函数，从而在不破坏和的精度的情况下获得更快的收敛速度。

明白非常重要的一点，从OUTCAR中直接提取的能量在不同赝势下不具备重复性，是不能够直接写到论文中去的。一般计算相对能量的时候只要保持取一样的值即可比较绝对能量没啥意义。

形成能

形成能很简单，参考高中化学形成焓的计算可以知道对于反应A+B=C+D，那么形成能=E(C)+E(D)-E(A)-E(B)
同样对于MoS2的生成焓我们可以知道： $\Delta$ H=E(MoS2)-E(Mo)-2E(S)
因此还需要计算Mo单质和S单质的单原子能量。计算过程和1.1完全一致，只需要获得结构，利用vaspkit生成KPOINTS和POTCAR文件，INCAR可以不用改。能量用F或者E0或者上面说的哪一项都可以，但要保持一致。代入计算即可。

除此之外还可以计算产生一个Mo的空位的形成能： $\Delta$ H=E(空位，MoS2)+E(Mo)-E（MoS2）
但是此时没有考虑缺陷的电荷。
可以计算分解能，剥离能，滑移能，表面能，偏聚能，吸附能等等。但是具体结构的优化会有一些处理上的不同。

单点能SP

参考： https://www.bilibili.com/read/cv10064345

单点计算是指对结果不做优化，直接输出能量。因此优化只会进行一步。

设置NSW = 0 也可以不设置，因为默认值就是0。 其他文件和参数不需要改变。提交即可。

对于一个结构在结构优化（NSW不等于0）时，第一步输出的1F和单点能的结果是一样的。

自洽计算

自洽计算其实也是单点计算，静态自洽计算前需要提供 较稳定体系 的晶格结构信息（即结构优化完的CONTCAR），从而通过电子自洽计算，通过自洽迭代求解薛定谔方程（微观中描述电子的状态，相当于宏观的运动方程））完整地计算出体系基态下费米能级（准确值，The fermi level location is accurate only in self-consistent calculation.）、电子的波函数、电荷密度等信息，可以直接分析原子间的键合作用，也可以在非自洽之后进一步分析晶体的电子结构和材料的相关性质。

因此对于电子结构计算之前必须进行一步自洽计算优化电子分布，以便加速后续计算。

{.line-numbers}

1
2
3

NSW=0  不再进行原子迟豫
LWAVE=.TRUE  （需要在WAVECAR中输出波函数）
LCHARGE=.TRUE    （需要在CHG和CHGCAR中输出电荷密度）

气体优化

在涉及单质结构优化时，可能会涉及到气体分子的优化。

一般情况下选择建一个相对大的格子，里面放入一个气体分子，格子的大小保证分子间不会影响。其他的和前面结构优化保持一致即可。
另一个办法就是找数据库， ChemScpider。网址：http://www.chemspider.com

IVDW参数

有两大类的范德华斯作用修正：

基于半经验的，包括了D2, D3, D3-BJ, TS, TS+SCS等等，这些都是在常用的交换关联泛函比如PBE的基础上，考虑色散力的作用，在PBE计算出来的总能基础上增加了额外半经验项，这一项需要设置一些参数，至于是哪种半经验公式，由IVDW的设置来决定，具体的参考VASP的手册。
范德华力计算方法，在DFT能量计算基础上增加范德华力修正。IVDW=10：DFT-D2方法；IVDW=11：DFT-D3方法。推荐首选更新的DFT-D3方法。
另一类就是vdW-DF，也就是修改交换关联项来考虑范德华修正。包括你所提到的optPBE-vdW, optB88-vdW, and optB86b-vdW等。它们的设置由修改GGA赋值，以及额外相关的一些参数来设置。这些在vasp手册上有明确的说明。
乱谈DFT-D：http://sobereva.com/83
DFT-D色散校正的使用：http://sobereva.com/210
大体系弱相互作用计算的解决之道：http://sobereva.com/214
LUSE_VDW：https://blog.csdn.net/Yaris_liu/article/details/124683272 https://zhuanlan.zhihu.com/p/47619708
包含vdW-DF，optPBE-vdW, optB88-vdW, optB88b-vdW, vdw-DF2。需要在所在的文件夹下放置一个名叫vdw_kernel.bindat的文件。

vdW-DF：采用该范德瓦尔斯修正需在INCAR中加上以下开关：
{.line-numbers}
1
2
3
GGA=RE
LUSE_VDW=.TRUE.
AGGAC=0.0000
optB88-vdW：采用该范德瓦尔斯修正需在INCAR中加上以下开关：

{.line-numbers} GGA

PARAM1=0.1833333333
PARAM2=0.2200000000
LUSE_VDW=.TRUE.
AGGAC=0.0000

optB86b-vdW：采用该范德瓦尔斯修正需在INCAR中加上以下开关：

{.line-numbers} GGA

PARAM1=0.1234
PARAM2=1.0000
LUSE_VDW=.TRUE.
AGGAC=0.0000

vdW-DF2：采用该范德瓦尔斯修正需在INCAR中加上以下开关：

{.line-numbers} GGA

1
2
3

LUSE_VDW=.TRUE.
Zab_vdW=-1.8867
AGGAC=0.0000

在考虑DFT-D2或者DFT-D3计算的过程中可同时考虑自旋轨道耦合（SOC），但是在vdW-DF、optB88-vdW 、optB86b-vdW和vdW-DF2中无法同时考虑自旋轨道耦合进行计算。

DFT+U

为什么加U?
对于含有d、f轨道电子的强关联体系(Hubbard模型)，电子间存在强烈的在位库仑相互作用，而交换相关泛函中的局域密度近似（LDA）或广义梯度近似（GGA）对电子之间的强在位库仑相互作用描述不准，此时可将研究体系的交换相关泛函分成两部分计算，也就是用DFT+U的方法进行求解，公式如下：
其中一部分电子用DFT算法（如LDA，GGA）等可以比较准确地描述；另外的d、f轨道电子通过引入Hubbard项得到正确的描述（U值就是考虑了同一个原子自旋相反的局域电子之间的库仑排斥）

如何加U?
DFT计算时输入文件的U值由以下参数确定：

{.line-numbers}

LDAU= .TRUE.|.FALSE.：开启/关闭+U功能，默认值为.FALSE.；
LDAUTYPE=1|2|4指的是+U的类型，默认值是2；其中1为Liechtenstein等提出的旋转不变LSDA+U方法；2为 Dudarev等提出的简化 LSDA+U方法；4与1类似, 但不考虑LSDA交换劈裂；
LDAUL=-1|1|2|3分别对应不加U、p、d、f轨道加U；
LDAUU、LDAUJ分别指定电子库伦相互作用项和交换相互作用项（U和J值）；
LMAXMIX =2/4/6：默认为2，加U计算时该值需大于轨道量子数，对于含有d轨道或f轨道电子的体系需对应增加至4或6。

注意：
1）必须为每一类原子指定一个U值，其顺序与结构文件POSCAR中元素排序一致；
2）同一种元素在不同的体系中U值不同；
3）U在所有电子中都存在，不管自旋是否相同，J只存在于自旋相同的电子上；当LDAUTYPE =2时，U-J的差值才有意义，即有效的U参数。
该参数需要测试或者参考文献。

磁性体系

很多体系都有磁性，磁性体系的计算需要打开自旋。以FeO为例：

{.line-numbers}

1 2	ISPIN=2 MAGMOM=原子数目初始磁矩，eg 32 45 （注：号前后不要有空格，且对应顺序）

对于一些简单的磁性体系，我们可以直接使用ISPIN=2, MAGMOM不必进行设置。对于复杂的体系则需要设置。初始值不需要与实验值完全一致，一般设置1.5倍即可。

有4个Fe原子，如果考虑铁磁态，则可以设置 MAGMON= 4*3 4*0 （Fe初设的磁矩为3，O元素为0）
有4个Fe原子，如果考虑反铁磁态，则可以设置 MAGMON=3 -3 3 -3 0 0 0 0或者MAGMON=3 3 -3 -3 0 0 0 0

表面优化

1. Bulk计算，获得稳定构型。

注意LREAL、PREC、ENCUT、EDIFF、ISMEAR保持一致。

2. 考虑不同晶面、层数。

对于表面结构，有以下几个需要注意的:

xy 方向上表面的大小（15 Å或者100原子左右）;
A 这个影响表面吸附物种的覆盖度;
B 影响体系的尺寸大小和计算时间;
C 不同的大小需要选取对应的K点；回想一下我们前面提到的经验规则。
不同的晶面，(111)， (100)，(110);
A 这取决于你研究的方向;低指数面、解离面。
B 不同晶面的表面能;
C 不同晶面的表面结构，反应活性等。
表面结构的层数（3+3以上，6Å+6Å）
A 层数多了，原子变多，体系在 z 方向的尺寸增加，也会影响计算速度；
B 计算中需要弛豫的层数；
C 不同层数对你要计算的性质会产生影响，比如表面能；
D 不同晶面需要的层数也不尽相同，一般开放的表面需要更多的层数；
E 根据自己感兴趣的性质，选择合适的层数，也就是需要你去测试一番。
Slab模型有两种，一种是上下表面对称的，一种是非对称的。对称性结构往往需要很多层，体系较大。非对称的结构体系较小，但存在偶极矩的影响，要注意加LDIPOL 和IDIPOL这两个参数来消除。1、2、3分表代表在x、y和z 方向上进行校正，4代表在所有方向上。真空层15-20埃。

3. 工具：MS、ASE、P4vasp、vaspkit。

MS操作。

一般用conventional cell。
MS操作。Windows
VESTA将CONTCAR转成cif,导入MS。
build-surface-cleave surface 切面
Surface mesh-surface vectors 可调整vector。

Surface box-top/thickness，调整终端和厚度。
Build-crystals-build vaccum slab-图2-build。

Lattice parameters-advanced-re-orient to startard把真空层沿Z轴图3。

转换成POSCAR，可优化后进一步扩胞。

ASE操作。python

块体优化后
module load python/3.7.3
python

from ase.io import read
from ase.io import write
from ase.build import surface
a=read(”CONTCAR”,format=’vasp’)
b=surface(a,(1,1,1),6,vacuum=7.5) 6层，真空层15Å
write(”POSCAR”,b,format=’vasp’)
Ctrl+D （退出 python）

大师兄ASE脚本：切金属

python3 cssm.py

vaspkit

块体优化后复制成POSCAR
vaspkit –task 803
生成SLAB***.vasp

4 表面弛豫计算

单点能

转换成POSCAR单点计算。或表面优化时第一个离子步的能量，但要保证这一步收敛。
KPOINTS设置，20-30。
INCAR设置，NSW=0或默认0.偶极矩矫正 LDIPOL = .TRUE. ; IDIPOL = 3。
POTCAR准备，提交任务脚本。完成后 grep without OUTCAR 获取单点能。

表面优化

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)

1ev/$Å^2$ =16.02J/$m^2$
$E_{surf}$是未弛豫的表面能量,A是表面积a*b。

吸附能扫描

算吸附能时一般只需要算分子在高对称点的吸附能，比如算Li在铜表面的吸附能，我们只需算Li在Cu(111）面的Top，HCP，FCC和Bridge位的吸附能即可。但是这些离散的点没法帮你构建一个“势能面”，也就是无法得到吸附能在表面究竟是怎么分布的. 下面的文献里DOI：10.1063/1.4901055有一张图非常生动，讲的分别是Li在Li(001)以及Mg在Mg(0001)表面的吸附势能分布，我们不仅可以知道哪些位置有利于吸附，还可以根据吸附能分布，描绘出分子在表面的扩散路径，这有助于我们使用NEB方法算扩散能垒。

图片

其实思路很简单，就是将表面网格化，算分子在格点的吸附能，再画出等高线图。但是大部分吸附位点其实是不稳定的，所以我们采取固定吸附分子x,y方向，在z方向弛豫以达到吸附平衡的策略。

其实手动撒点，再采集数据也是可行的，但是会比较麻烦，因此作者根据实际需要开发了一款脚本scan_adsorption_energy用于自动完成这个过程。脚本使用Python编写，需要numpy和matplotlib第三方库。我们首先算好一个吸附例子得到CONTCAR,这个可以让我们得到吸附分子的元素信息和理想的吸附高度。如下图，吸附的Li是第97号，也是最后一个原子，接着就可以变换这个结构得到不同的吸附结构。

图片

首先我们需要在目录下准备以下几个文件CONTCAR,INCAR,KPOINTS,POTCAR以及排队脚本文件(我的是vasp.pbs)，然后使用脚本产生不同的结构：

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里面包含了超胞数，选择移动的原子以及晶格常数信息。在所有的计算完成后，同样使用该脚本提取数据。

python scan_adsorption_energy.py -e
脚本会自动从每个结构文件夹内读取能量数据和吸附原子的坐标。生成两个dat文件new_direct.dat和new_dat-3D.dat。前者是原子的分数坐标和总能量，后者有4列，分别为原子的笛卡尔x分量坐标，y分量坐标，z分量坐标以及总能量。我们需要使用第1,2,4列绘制等高线图。如果需要绘制吸附能，需要自己手动将第四列的数据减去吸附原子和表面的能量。脚本会尝试使用matplotlib绘制等高线图，如下图，你也可以使用Origin Pro绘图。

图片

脚本可以在我的Github仓库下载(https://github.com/tamaswells/VASP_script/blob/master/scan_adsorption_energy.py).

固定基矢优化结构

在第一性原理计算时，二维材料的建模会采用在某一方向（一般是z轴方向）添加15-20 Å的真空层的方式来抵消层与层之间的相互作用。但是，当我们优化晶格常数的时候会遇到一个问题：那就是使用ISIF=3来优化平面内晶格常数（一般是a和b）时，c轴方向的晶格长度总是随着优化变的越来越小。

采用ISIF=4的固定晶胞体积的方法来计算。这样采用实验参数建的晶胞，由于a和b方向的变化，c轴也会变化。解决了越是优化c轴越小的问题，但是计算量还是很大。
编译vasp，就自己编译了固定z轴的脚本，专门计算二维材料。这种方法优点是直接换脚本就能解决上面优化时c轴缩小的问题，但是我想固定两个轴就需要再编译一份脚本。比较麻烦。

修改constr_cell_relax.F实现固定基矢优化结构的原理：

FCELL contains the forces on the basis vectors. These forces are used to modify the basis vectors according to the following equations:

FCELL 是个3 * 3的矩阵，包含了晶格矢量的受力，当ISIF = 3时，晶格矢量会按照如下的方法更新：

A_OLD(1:3,1:3)=A(1:3,1:3) ! F90 style
DO J=1,3
DO I=1,3
DO K=1,3
A(I,J)=A(I,J) + FCELL(I,K)*A_OLD(K,J)*STEP_SIZE
ENDDO
ENDDO
ENDDO
当通过OPTCELL文件把晶格矢量在某个方向的受力改成0以后，

IF (ICELL(I,J)==0) FCELL(I,J)=0.0

在结构更新的时候，该方向就不会改变。完成编译后，在 INCAR 文件中设置 ISIF=3，便可实现固定 Z 轴的结构优化。
对于其他方向晶格基矢的修改同理：对于a方向基矢，将 FCELL(3,I) 和 FCELL(I,3) 分别改为 FCELL(1,I) 和 FCELL(I,1) ；对于b方向基矢，则分别改为 FCELL(2,I) 和 FCELL(I,2) ；固定两个基矢则应该同时将两个方向对应的矩阵元设置为0。或参考 VASP并行可执行软件包，可对晶胞参数进行部分优化可实现添加参数文件实现不同基矢方向固定。（未做测试）

肖承诚博士的补丁，按照肖博士的方法安装以后，在INCAR中添加如下参数即可。

IOPTCELL = xx xy xz yx yy yz zx zy zz

上面的九个参数对应晶格常数矩阵。

1 xx xy xz

2 yx yy yz

3 zx zy zz

1代表驰豫，0代表固定

此时固定c轴的设置就是如下：

IOPTCELL = 1 1 0 1 1 0 0 0 0

下载地址：https://github.com/Chengcheng-Xiao/VASP_OPT_AXIS

滑移能和剥离能

材料的塑性行为主要是由于材料的结构和化学键决定的。材料发生塑性形
变在原子尺度上可以分解成两个过程。

一个是由于原子、位错或者晶界沿着某
一特定方向自由运动；第二个是在这种运动过程中存在束缚力保证在其它方向
不会产生裂纹，导致材料的断裂。因此任何一个塑性材料必须满足两个要求，
一是对某些原子层存在势垒相对低的滑动方式，以确保原子、位错或者晶界可
以沿着这些方向滑动。另一个是在这些滑动过程中，滑移平面两侧的原子依然
存在某种相对强的作用力，以确保滑移过程中材料的完整性。如金属键，其自
由电子容易发生变形和转移，在滑移过程中金属键不断的断裂和产生，从而导致了金属的滑移能量很小。而金属键的键强度很大，保证的金属滑移过程的完
整性。而离子键和共价键由于配位数和方向性等要求，在滑移时必须先发生键
的断裂，从而阻碍滑移的进行。由范德华力连接的层状材料，虽然层与层之间
容易发生滑移，但其层间作用力较弱，材料在滑移过程中容易发生破坏。

为了研究 α-Ag2S 的塑性行为，我们通过第一性原理对 α-Ag2S 的滑移过程进
行了计算。其中，我们计算的滑移面为（100）方向，滑移方向为[001]方向。其
滑移过程如图 6.10 所示，α-Ag2S 的滑移面沿着能量最低点进行滑移，滑移面不
仅发生相对移动，其面间距也会发生。因此，整个滑移路径成波浪形。我们将滑移一个周期（滑移一个晶格常数）均分 12 步。每一步我们计算系
统总能量与滑移面之间距离的关系，对于塑性材料而
言，它需要同时具有足够小的滑移能来保证材料的滑移可以进行；以及大的解离
能以至于材料不会发生断裂。

1.1 结构优化1

发表于 2023-05-09 更新于 2024-06-05 分类于 DFT学习

这里主要分享怎样快速开始第一个结构优化过程。

Linux常用命令

如果是个人服务器，需要安装Linux系统或者Windows+虚拟机+Linux；如果租用超算，安装对应超算客户端或者SSH。
这部分内容需要自己熟练，对于命令的意义在操作中就会明白。https://www.runoob.com/linux/linux-command-manual.html
常用命令如下：
原文链接：https://blog.csdn.net/m0_46422300/article/details/104645072

目录操作

1 切换目录（cd）

{.line-numbers}

cd /                 //切换到根目录
cd /bin              //切换到根目录下的bin目录
cd ../               //切换到上一级目录 或者使用命令：cd ..
cd ~                 //切换到home目录
cd -                 //切换到上次访问的目录
cd xx(文件夹名)       //切换到本目录下的名为xx的文件目录，如果目录不存在报错
cd /xxx/xx/x         //可以输入完整的路径，直接切换到目标目录，输入过程中可以使用tab键快速补全

2 查看目录（ls）

{.line-numbers}

ls                   //查看当前目录下的所有目录和文件
ls -a                //查看当前目录下的所有目录和文件（包括隐藏的文件）
ls -l                //列表查看当前目录下的所有目录和文件（列表查看，显示更多信息），与命令"ll"效果一样
ls /bin              //查看指定目录下的所有目录和文件

3 创建目录（mkdir）

{.line-numbers}

1 2	mkdir tools //在当前目录下创建一个名为tools的目录 mkdir /bin/tools //在指定目录下创建一个名为tools的目录

4 删除目录与文件（rm）

{.line-numbers}

rm 文件名              //删除当前目录下的文件
rm -f 文件名           //删除当前目录的的文件（不询问）
rm -r 文件夹名         //递归删除当前目录下此名的目录
rm -rf 文件夹名        //递归删除当前目录下此名的目录（不询问）
rm -rf *              //将当前目录下的所有目录和文件全部删除
rm -rf /*             //将根目录下的所有文件全部删除【慎用！相当于格式化系统】

5 修改目录（mv）

{.line-numbers}

1
2
3

mv 当前目录名 新目录名        //修改目录名，同样适用与文件操作
mv /usr/tmp/tool /opt       //将/usr/tmp目录下的tool目录剪切到 /opt目录下面
mv -r /usr/tmp/tool /opt    //递归剪切目录中所有文件和文件夹

6 拷贝目录（cp）

{.line-numbers}

1 2	cp /usr/tmp/tool /opt //将/usr/tmp目录下的tool目录复制到 /opt目录下面 cp -r /usr/tmp/tool /opt //递归剪复制目录中所有文件和文件夹

7 搜索目录（find）

{.line-numbers}

1	find /bin -name 'a*' //查找/bin目录下的所有以a开头的文件或者目录

8 查看当前目录（pwd）

{.line-numbers}

1	pwd //显示当前位置路径

文件操作

1 新增文件（touch）

{.line-numbers}

1	touch a.txt //在当前目录下创建名为a的txt文件（文件不存在），如果文件存在，将文件时间属性修改为当前系统时间

2 删除文件（rm）

{.line-numbers}rm 文件名link

1	rm -f 文件名 //删除当前目录的的文件（不询问）

3 编辑文件（vi、vim）

{.line-numbers}

1	vi 文件名 //打开需要编辑的文件

–进入后，操作界面有三种模式：命令模式（command mode）、插入模式（Insert mode）和底行模式（last line mode）
命令模式
-刚进入文件就是命令模式，通过方向键控制光标位置，
-使用命令”dd”删除当前整行
-使用命令”/字段”进行查找
-按”i”在光标所在字符前开始插入
-按”a”在光标所在字符后开始插入
-按”o”在光标所在行的下面另起一新行插入
-按”：”进入底行模式
插入模式
-此时可以对文件内容进行编辑，左下角会显示 “– 插入 –””
-按”ESC”进入底行模式
底行模式
-退出编辑： :q
-强制退出： :q!
-保存并退出： :wq

操作步骤示例

1.保存文件：按”ESC” -> 输入”:” -> 输入”wq”,回车 //保存并退出编辑
2.取消操作：按”ESC” -> 输入”:” -> 输入”q!”,回车 //撤销本次修改并退出编辑

补充

{.line-numbers}

1 2	vim +10 filename.txt //打开文件并跳到第10行 vim -R /etc/passwd //以只读模式打开文件

4 查看文件

{.line-numbers}

cat a.txt          //查看文件最后一屏内容
less a.txt         //PgUp向上翻页，PgDn向下翻页，"q"退出查看
more a.txt         //显示百分比，回车查看下一行，空格查看下一页，"q"退出查看
tail -100 a.txt    //查看文件的后100行，"Ctrl+C"退出查看

5 grep

{.line-numbers}

1
2
3

grep -i "the" demo_file              //在文件中查找字符串(不区分大小写)
grep -A 3 -i "example" demo_text     //输出成功匹配的行，以及该行之后的三行
grep -r "ramesh" *                   //在一个文件夹中递归查询包含指定字符串的文

文件权限

1 权限说明

文件权限简介：’r’ 代表可读（4），’w’ 代表可写（2），’x’ 代表执行权限（1），括号内代表”8421法”
##文件权限信息示例：-rwxrw-r–
-第一位：’-‘就代表是文件，’d’代表是文件夹
-第一组三位：拥有者的权限
-第二组三位：拥有者所在的组，组员的权限
-第三组三位：代表的是其他用户的权限

2 文件权限

普通授权 chmod +x a.txt
8421法 chmod 777 a.txt //1+2+4=7，"7"说明授予所有权限

打包与解压

1 说明

{.line-numbers}

.zip、.rar        //windows系统中压缩文件的扩展名
.tar              //Linux中打包文件的扩展名
.gz               //Linux中压缩文件的扩展名
.tar.gz           //Linux中打包并压缩文件的扩展名

2 打包文件

tar -zcvf 打包压缩后的文件名要打包的文件
参数说明：z：调用gzip压缩命令进行压缩; c：打包文件; v：显示运行过程; f：指定文件名;
示例：

{.line-numbers}

1	tar -zcvf a.tar file1 file2,... //多个文件压缩打包

3 解压文件

{.line-numbers}

tar -zxvf a.tar                      //解包至当前目录
tar -zxvf a.tar -C /usr------        //指定解压的位置
unzip test.zip             //解压*.zip文件 
unzip -l test.zip          //查看*.zip文件的内容

原文链接：https://blog.csdn.net/m0_46422300/article/details/104645072

其他准备

安装好VASP、vaspkit，个人电脑安装好VESTA、MS或者qvasp等建模可视化软件中的一个。如果是租用超算，还需要准备好提交任务的脚本文件，一般工程师会帮我们编好，软件安装只需要给他们提供安装包或官网即可。

vaspkit运行的五种方式

方法一 ：交换式菜单模式，直接输入vaspkit；

方法二 ：命令+参数模式，适合批处理。如vaspkit -task 102 -kpr 0.03，详细介绍见vaspki-help, 注意部分功能不支持这种模式；

方法三 ：利用bash管道传递参数,如echo -e “102\n2\n0.04\n”| vaspkit

方法四 ：同样利用bash管道传递参数, (echo 102; echo 2; echo 0.04)|vaspkit

方法五 ：vi cmd.in（文件名任意，非必须是cmd.in）包含以下三行内容

{.line-numbers}

1
2
3

102
2
0.04

然后运行vaspkit < cmd.in，产生KPOINTS文件，其他命令格式类似。

输入文件

{.line-numbers}

INCAR      in    **
STOPCAR    in
stout      out
POTCAR     in     **
KPOINTS    in     **
IBZKPT     out
POSCAR     in     **
CONTCAR    out
CHGCAR     in/out
CHG        out
WAVECAR    in/out
TMPCAR     in/out
EIGENVAL   out
DOSCAR     out
PROCAR     out
OSZICAR    out
PCDAT      out
XDATCAR    out
LOCPOT     out
ELFCAR     out
PROOUT     out

必要的输入文件有四个：INCAR，KPOINTS，POSCAR，POTCAR 。
- INCAR 告诉VASP算什么，怎么算。
- KPOINTS 包含计算的倒空间K点信息。
- POSCAR 是结构信息，包括元素、晶胞参数及各个原子的坐标信息。
- POTCAR 元素的赝势文件，描述体系中对应的原子核和电子的相关信息。
- 提交任务的脚本或者命令，需要你自己准备。

下面由简单到复杂介绍

POSCAR

结构文件信息格式如下：

{.line-numbers}

Mo S2                  #结构名称，可不写或随便写，不能没有这一行
1.0                   # universal scaling parameters，缩放系数
    3.1500000954     0.0000000000     0.0000000000 # 晶格矢量 a(1) 
    -1.5750000477    2.7279801045     0.0000000000 # 晶格矢量 a(2)
     0.0000000000    0.0000000000    12.3000001907 # 晶格矢量 a(3)
   Mo    S                                         # 元素名称，种类
    2    4                                         #对应元素的原子个数，顺序要对应
Direct                                             #
     0.333333345       0.666666689       0.250000000
     0.666666624       0.333333323       0.750000019
     0.333333345       0.666666689       0.620999986
     0.666666624       0.333333323       0.379000014
     0.666666624       0.333333323       0.121000053
     0.333333345       0.666666689       0.878999976

对于自己建模的结构，按照该格式自己写；
对于Materials projiect数据库可以直接下载POSCAR格式；
对于.cif文件，可以导入 VEASTA 等输出为POSCAR的格式，或者用 vaspkit 的411/412功能。

Mind：

POSCAR文件非常重要，大部分报错都是结构不合理导致的。一个好的初始结构不仅可以加快收敛速度，也会增加准确性。文件中元素的顺序非常重要。

KPOINTS

文件一般格式如下：

{.line-numbers}

K-POINTS      #  第一行随便写或者空着不写，但不能没有这一行
0             # 零，表示自动生成
Gamma         #  gamma点centered 
1 1 1         #  1*1*1格子，K点密度
0 0 0         #  S1 S2 S3，shift的值，一般保持 0 0 0 不变。

大家需要修改的是第三行和第四行，其他可不改。

第三行：VASP只认第一个字母，大小写均可。当然这一行也可以直接写字母G或者g。G表示的是以gamma点为中心生成网格。
另外一种是原始的Monkhorst-Pack 网格，两者的区别是 M 或者 m 在 G的基础上在三个方向上平移了1/（2N）个单位。G，也叫 gamma centered Monkhorst-Pack Grid；所以，gamma centered 只是MP网格的一种特殊情况。
因此这一行也可以一直用G不变。
第四行，在xyz三个方向上生成对应数目的K点；__这个一般需要测试__。如果开始不想测试，官网也给出了经验方法，
- d区金属，k*a ~ 30 Å
- 普通金属，k*a ~ 25 Å
- 半导体，k*a ~ 20 Å
- 绝缘体，k*a ~ 15 Å
其中a对应三个方向的晶胞参数，也就是对应方向的k点与该方向晶胞参数乘积最好在20-30之间。

如果已经安装了vaspkit程序，可以通过102功能自动生成KPOINTS文件。

POTCAR

POTCAR赝势文件相对来说最复杂也最简单。其中有很多信息，大部分信息我们不需要关注，一般只需要关注前几行。特别是赝势的种类信息。以Mo为例：

{.line-numbers}

PAW_PBE Mo 08Apr2002                   
   6.00000000000000     
 parameters from PSCTR are:
   VRHFIN =Mo: 4p5s4d
   LEXCH  = PE
   EATOM  =   217.5176 eV,   15.9871 Ry

   TITEL  = PAW_PBE Mo 08Apr2002
   LULTRA =        F    use ultrasoft PP ?
   IUNSCR =        1    unscreen: 0-lin 1-nonlin 2-no
   RPACOR =    2.200    partial core radius
   POMASS =   95.940; ZVAL   =    6.000    mass and valenz
   RCORE  =    2.750    outmost cutoff radius
   RWIGS  =    2.750; RWIGS  =    1.455    wigner-seitz radius (au A)
   ENMAX  =  224.584; ENMIN  =  168.438 eV

除此之外，一种元素也有很多赝势，例如Mo：

{.line-numbers}

[scf1355@ln151%nc-l potpaw_PBE.52]$ ls Mo*
Mo:
POTCAR  PSCTR

Mo_pv:
POTCAR  PSCTR

Mo_sv:
POTCAR  PSCTR

Mo_sv_GW:
POTCAR  PSCTR

如果没有特别的需求，直接采用VASP官网推荐的即可。参考链接：
https://cms.mpi.univie.ac.at/vasp/vasp/Recommended_PAW_potentials_DFT_calculations_using_vasp_5_2.html
不知道选哪一个就选第一个，使用方法如下：

{.line-numbers}

1	cat ~/pot/Mo/POTCAR ~/pot/S/POTCAR > POTCAR

或者用vaspkit的103命令自动生成。

Mind：

该文件看似复杂，使用时只需要查看元素种类和POSCAR是否对应即可。
可以通过grep PBE POTCAR查看顺序是否正确。

{.line-numbers}

[scf1355@ln151%nc-l 2h]$ grep PBE POTCAR 
  PAW_PBE Mo 08Apr2002                   
   TITEL  = PAW_PBE Mo 08Apr2002
  PAW_PBE S 06Sep2000                    
   TITEL  = PAW_PBE S 06Sep2000

INCAR

这个控制计算过程的文件最为复杂，手册上相关参数多达三百多个。因此前期我们记住一点：不知道含义的参数不要加。不加对于有些参数意味着用默认值，一般是合理的。INCAR参数查找参考官方：https://www.vasp.at/wiki/index.php/Category:INCAR_tag

贴一个我的结构优化的INCAR：

{.line-numbers}

ISTART =  0            (Read existing wavefunction; if there)
ISPIN  =  2           (Non-Spin polarised DFT) 
ICHARG =  2         (Non-self-consistent: GGA/LDA band structures)
ENCUT  =  500        (Cut-off energy for plane wave basis set, in eV)
PREC   =  A       (Precision level)
LWAVE  = .F.        (Write WAVECAR or not)
LCHARG = .F.        (Write CHGCAR or not)
ISMEAR =  0            (Gaussian smearing; metals:1)
SIGMA  =  0.05         (Smearing value in eV; metals:0.2
EDIFF  =  1E-08        (SCF energy convergence; in eV)
EDIFFG = -1E-02        (Ionic convergence; eV/AA) 
NSW    =  100          (Max ionic steps)
IBRION =  2            (Algorithm: 0-MD; 1-Quasi-New; 2-CG)
POTIM = 0.1
ISIF  =  3            (Stress/relaxation: 2-Ions, 3-Shape/Ions/V, 4-Shape/Ions)

上面参数详细介绍

ISTART ：初始波函数。WAVECAR存在默认1，不存在默认0。
0：随机生成波函数
1：从 WAVECAR 读取波函数，若读取失败则随机生成。
ICHARG ：初始电荷密度。当ISTART=0默认2，其他默认1。
0：通过 WAVECAR 生成电荷密度，若读取失败则通过原子电荷密度叠加生成
1：从 CHGCAR 读取电荷密度
2：通过原子电荷密度叠加生成电荷密度
11：从 CHGCAR 读取电荷密度，且在（非自洽）循环中保持不变
ISPIN ：自旋极化。默认1，关闭。不知道是否打开，可以打开，结果不会有影响。
1：非自旋极化（每个轨道上自旋向上和自旋向下的电子数量相等，适合非磁性体系）
2：自旋极化（适合铁磁、反铁磁体系）
ENCUT ：平面波截断能（单位 eV）
默认元素中最大的ENMAX。一般为 POTCAR 中最大ENMAX 参数值的 1.0 到 1.3 倍。500适用大部分体系。
PREC ：总体计算精度，默认Normal。可设为 Low、Med、High、Normal、Single、Accurate。具体参照手册，设为A或者Accurate足够精确。
LWAVE LCHARG ：控制WAVECAR和CHGCAR输出，默认输出，在结构优化可以关闭，节省储存空间。
ISMEAR ：轨道分数占据
−5：四面体方法（适合半导体和绝缘体）
0：高斯方法（适合导体、半导体和绝缘体），展宽由 SIGMA 确定
N：Methfessel-Paxton 方法（适合导体）
EDIFF ：自洽循环收敛标准（系统能量变化），单位为 eV。默认1E-4，可以设的更精确一些。一般设置为 1E-6或更高。
__EDIFFG__：离子位置优化收敛标准。默认EDIFF*10。正值为系统能量变化，负值为原子上残余力。随着体系维数下降，可适当增加原子残余力。
NSW ：离子位置优化最大步数（IBRION=1、2）。默认0。设为100足够大，也可以继续加大。
分子动力学开启时代表分子动力学模拟步数（IBRION=0）
IBRION ：离子位置优化算法。默认NSW=0/1时为-1，其他默认0。
-1：离子不移动
0：分子动力学模拟
1：准牛顿法
2：共轭梯度法
5：振动频率计算
6：弹性常数计算
7/8：DFPT方法
POTIM ：有限差分法步长。当IBRION=1/2/3默认值 0.5。当IBRION=0（MD）时必须指定。
ISIF ：需要结构优化变量。默认IBRION=0时为0，其他默认2。结构优化一般设为3，全部优化。

Mind

上面四个输入文件都是纯文本文件，且文件名不能更改，格式不能错误。
vaspkit可以同时产生INCAR、POTCAR、KPOINTS文件。

提交任务

在超算上提交或者集群上提交，只需要把四个输入文件和提交任务的脚本文件放在同一文件夹下，运行脚本即可。
在个人服务器、虚拟机提交可以通过nohup mpirun -np 核数 vasp_std > std.out &提交

Mind

提交后经常查看作业收敛情况，有问题及时暂停，避免浪费时间

查看任务完成

实时查看任务日志tail -f 日志文件
任务停止后，判断任务是否正常结束，输入grep reached 日志文件或者tail 日志文件
如果输出reached required accuracy - stopping structural energy minimisation代表正常结束，否则很可能没有收敛或者出错中断。

输出文件

任务一旦提交，会立刻生成输出文件，不过此时可能是空的。运行片刻大部分文件会产生内容。
https://www.vasp.at/wiki/index.php/Category:Output_files
其中我们先关注几个输出文件。

CONTCAR

优化后的结构文件，与POSCAR相对应。为最后一个离子步的结构。当然如果要提高精度或者续算，可以cp CONTCAR POSCAR继续计算。

OSZICAR

包含自洽计算中能量收敛等信息。

OUTCAR

这个是VASP最主要的输出文件，包含计算过程中大量信息，依次主要包括：
VASP版本；
计算开始时间和并行性CPU数；
赝势信息；
最近邻列表；
对称性信息；
晶格信息和k点坐标;
INCAR中读入的参数和其他大部分的默认参数值；
平面波个数和FFT信息；
每一步离子步数和其中每一个电子自洽的时间、内存、能量等信息；
自洽完成后的费米能和能量本征值；
应力；力；
电荷数和磁矩；
程序运行时间。
总而言之想找而不知道去哪里找，就在这里。

IBZKPT

由VASP自动生成，包含所有不等价k点的坐标和权重以及可能的四面体链接情况，它的格式与KPOINTS文件的格式是完全一样的。

WAVECAR

波函数文件，二进制文件，不能直接用文本编辑器直接打开。
波函数文件较大，可以通过(LWAVE)来控制输出。结构优化不需要输出。
它也可以作为输入文件，为后续计算提供初始波函数。

CHGCAR和CHG

电荷密度文件，包含晶格矢量，原子坐标，总电荷密度以及PAW的单中心占据情况等信息。该文件既是输出文件(由LCHARG控制)，也可以是输入文件。结构优化不需要输出。
在计算能带或者态密度时，需要读入电荷密度进行非自洽计算。

总结

输入文件KPOINTS比较简单，可以手动，注意格式；INCAR参数众多，直接参考，不知道的参数不加；POTCAR文件比较大，我们只关注赝势对应的元素顺序与POSCAR是否一致；POSCAR结构文件VESTA直接输出。提交任务脚本由超算代理提供，不同超算不同。
任务提交后关注是否正常结束，另外关注计算的时间等。

0.5 声子软件安装

发表于 2023-05-08 更新于 2023-12-19 分类于软件安装

Anaconda安装

Anaconda-Linux版本安装包【包含了phonopy安装所需的依赖包】
https://mirrors.tuna.tsinghua.edu.cn/help/anaconda/
https://www.anaconda.com/download/
将下载的准备文件上传到Linux系统
cd 到Anaconda文件的目录下，赋予权限，运行下列命令
sh [Anaconda文件名]
例如：
sh Anaconda3-2019.10-Linux-x86_64.sh

没什么特殊要求的话，之后一路回车即可安装成功。

输入conda list 如果显示已安装的包名和版本号便已安装成功。
修改配置:
输入vi ~/.bashrc （配置python环境）
在配置文件末尾中加入下列语句（引号内的路径更换为自己的相对应路径即可）
export PATH=”/路径/anaconda3/bin:$PATH”

添加后保存退出，运行 source ~/.bashrc 使配置生效

命令终端输入 python 如果出现python版本号则配置成功，输入 exit() 即可退出python命令模式。

adaconda使用介绍：

1 管理conda
(1) 检查conda
conda –version Conda会返回你安装Anaconda软件的版本。

(2) 升级conda
conda update conda Conda会检查可升级版本，并同时显示可升级的包。

2 管理运行环境
(1) 创建新运行环境
conda create –name snowflakes biopython 创建名为snowflakes的新运行环境，其中包含包biopython

(2) 激活新运行环境
Windows: activate snowflakes
Linux and macOS: source activate snowflakes
(3) 查看所有运行环境
conda info –envs 查看所有运行环境

base /home/username/Anaconda3
snowflakes * /home/username/Anaconda3/envs/snowflakes
当前激活的运行环境有且只有一个，且前面有星号(*)

(4) 查看当前运行环境
conda info -envis
当前运行环境会显示在括号中：(snowflakes)

(5) 修改当前运行环境为非激活状态
Windows: deactivate
Linux and macOS: source deactivate
(6) 复制运行环境
conda create -n flowers –clone snowflakes

(7)删除运行环境
conda remove -n flowers –all
删除运行环境之后可以通过命令查看运行环境是否删除： conda info -e

3 管理Python
(1) 创建包含Python运行的新运行环境
conda create –name snakes python=3.5 创建名为snakes并包含Python3.5的新运行环境

(2) 激活新运行环境
Windows: activate snakes
Linux and macOS: source activate snakes
(3) 查看所有运行环境
conda info –envs 查看所有运行环境

base /home/username/Anaconda3
snakes * /home/username/anaconda3/envs/snakes
snowflakes /home/username/Anaconda3/envs/snowflakes
当前激活的运行环境有且只有一个，且前面有星号(*)

(4) 查看当前运行环境
conda info -envis 或者 conda info -e
当前运行环境会显示在括号中：(snakes)

(5) 查看Python版本
pyhton –version

(6) 修改当前运行环境为非激活状态
Windows: deactivate
Linux and macOS: source deactivate
4 管理包
(1) 查看已安装包
conda list

(2) 安装指定包到当前运行环境
conda install beautifulsoup4

(3) 查找可安装包
conda search beautifulsoup4

(4) 移除包
conda remove -n beautifulsoup4 –all

作者：ZHOUZAIHUI
链接：https://www.jianshu.com/p/6e4c045a12a1

Phonopy安装

下载phonopy：https://pypi.org/project/phonopy/#files
上传解压：tar -xzvf phonopy-2.17.0.tar.gz //解压phonopy文件
cd phonopy-2.17.0
python setup.py install 安装phonopy
安装完成后在命令端输入： phonopy 如果出现如下图案则安装成功

Installed /home/hptem001/intel/oneapi/intelpython/python3.7/lib/python3.7/site-packages/phonopy-2.21.0-py3.7-linux-x86_64.egg

Processing dependencies for phonopy==2.21.0
Searching for spglib>=2.0
Reading https://pypi.org/simple/spglib/

出现报错：查看phonopy2.21说明依赖要求

python>=3.8
numpy>=1.17.0
PyYAML>=5.3
matplotlib>=2.2.2
h5py>=3.0
spglib>=2.0
scipy (optional)
seekpath (optional)

需要h5py>=3.0 、spglib>=2.0，
在这里对于类似安装都有三种：

可以联网：pip install spglib或者spglib-2.2.0
不能联网已安装python：
下载：https://pypi.org/project/spglib/#modal-close 二进制文件，由于python3.7，所以选择CP37。
pip install spglib-2.2.0-cp37-cp37m-manylinux_2_17_x86_64.whl
不能联网已安装python：
下载压缩包解压，同样安装python setup.py install。
我在安装h5py老是报错，

安装cached-property。

https://pypi.org/project/cached-property/1.5.2/#files

再安装h5py即可。

thirdorder安装

下载：https://drive.google.com/file/d/1N7bwTFwsP75Dw4FH0O-2DEodE7TVw2ub/view
解压完之后进入文件夹修改setup.py文件
主要给出spglib的INCLUDE_DIRS和LIBRARY_DIRS即include和lib64文件夹的路径.
修改完成后运行python setup.py install

ShengBTE安装

下载：https://drive.google.com/file/d/1Kur3BviQjEby34j0yxbkMF45ADJf2U0J/view
解压完进入文件夹cp arch.make.example arch.make
修改前：
修改后：
在Src文件夹内执行命令make
测试：mpirun -np 16 ShengBTE 2>BTE.err >BTE.out
添加环境变量：

vim ~/.bashrc

#加入如下命令
export PATH=/$dir/src:$PATH #$dir：ShengBTE文件夹路径

#保存退出

source ~/.bashrc

0.4 VASP6服务器安装

发表于 2023-05-08 更新于 2023-12-19 分类于软件安装

系统安装

一、下载centos镜像

centos官方：https://vault.centos.org/ 清华镜像源下载https://mirrors.tuna.tsinghua.edu.cn/
常用的镜像文件类型介绍：
DVD ISO：普通光盘完整安装版镜像，可离线安装到计算机硬盘上，包含大量的常用软件（一般选择这种jing）。
Everything ISO：包含了完整安装版的内容，并对其进行补充，集成了所有软件。
Minimal ISO：这个版本为精简版的镜像，可以安装一个基本的CentOS系统，包含了可启动系统基本所需的最小安装包。
Netinstal：在线安装版本，启动后需要联网边下载边安装。

二、制作centos安装盘-u盘

这里推荐rufus:https://rufus.ie/downloads/
无需安装，如图

“引导类型选择”选择镜像文件，选择下载好的ISO文件，开始即可。
用软碟通打开镜像，然后选择 “ 启动”再选择 “写入硬盘映像” ，再选择自己的u盘之后，再选写入。注意：Windows系统限制了LABEL的长度为11，多出的部分被截断了，所以导致U盘的LABEL只有“CentOS 7 x8”11位(实际为CentOS 7 x86_64)。
u盘格式

1、FAT32格式：

这一种格式是任何USB存储设备都会预装的文件系统，属Windows平台的传统文件格式，兼容性很好。但是呢，“上帝”还是公平的，给你开了这一扇窗，必然会关上另一扇门。即便FAT32格式兼容性好，但它不支持4GB以上的文件，因此对于很多镜像文件或大视频文件之类的也会有无可奈何的状况。
一句话概括：FAT32格式兼容性好，但不支持4GB以上的文件
2、NTFS格式：
前面说到FAT32格式是Windows平台的传统文件格式，而NTFS格式却是Windows平台应用最广泛的文件格式。它的优点在于能够支持大容量文件和超大分区，且集合了很多高级的技术，其中包括长文件名、压缩分区、数据保护和恢复等等的功能。但人无完人，任何事物也相同，那NTFS格式又差在哪了呢？差就差在它会减短闪存的寿命。为什么？大家是否知道NTFS格式是针对机械硬盘设计的，它会对硬盘的读写操作做详细的记录，而闪存储存芯片的读写次数是有限的，若使用该格式就会让闪存造成很大的负担和伤害。易受伤的人心碎了还能长长久久吗？更何况是U盘呢！
一句话概括：NTFS格式支持大容量文件和超大分区，但对闪盘芯片有伤害
3、exFAT格式
虽然说到FAT32是传统文件格式，NTFS又是最广泛的，但老毛桃可以告诉大家：exFAT格式才是最适合U盘的文件格式，它是微软为了闪存设备特地设计的文件系统，是U盘等移动设备最好的选择之一。注意：SSD和U盘同为闪存，但SSD还是用NTFS格式为好！
exFAT格式的优点在于它能够增强台式机或笔记本和移动设备之间的互操作能力、同目录下最多65536个文件，且支持访问控制，还有的就是剩余空间分配表改善空间分配行。一句话概括：exFAT格式是最适合U盘的
在linux系统下centos支持VFAT格式和FAT32格式的U盘，而且VFAT是Windows和LINUX系统都可以使用的格式。如果U盘是NTFS格式的，而又不想切换格式，那么就需装个NTFS-3G的包就可以支持NTFS。

在U盘重装系统时，我们首要的便是制作一个U盘启动盘，经常捣鼓系统的伙伴们大都知道，制作前我们需要选择U盘的模式和格式，其中模式分为2种：USB-HDD和USB-ZIP;格式分为3种，也就是上述说到的。
USB-HDD模式： USB-HDD硬盘仿真模式，此模式兼容性很高，适用于新机型，但对于一些只支持USB-ZIP模式的电脑则无法启动，一般制作U盘启动盘选择该模式。
USB-ZIP模式：大容量软盘仿真模式，此模式在一些比较老的电脑上是唯一可选的模式，但对大部分新电脑来说兼容性不好，特别是2GB以上的大容量U盘。

三、u盘安装操作系统

BIOS设置问题：

U盘启动盘有分两种的，一个是legacy启动U盘，这种启动盘不能在UEFI模式下启动。另一种是UEFI启动盘，这种就可以直接在UEFI模式下识别到。
解决方法：
1、开机不停按Del键（也可能是其他：F2、F12、esc等），进入BIOS。
2、进入BIOS后，通过左右方向键切换到boot选项，设置Secure Boot为Disabled（禁用），表示关闭安全引导。
3、将U盘启动改为第一个，方法是将其上面的选项键盘输入-号，一直到U盘启动到第一个；选中enter。
4、最后按下F10保存设置退出，开机时再按U盘启动快捷键就可以识别到U盘启动盘了。

其他原因：

1、 U盘启动盘制作失败了，需要重新制作一遍启动盘。用软碟通容易失败，推荐Rufus。

2、 USB接口接触不良，可以换几个接口尝试。

3、安装dracut报错解决方法：安装后进入了dracut:/#模式
（1）.查看你的U盘名称:
dracut:/# cd /dev
dracut:/# ls
ls然后找到sda 是硬盘对应的文件名，sdb是U盘对应的文件名，可以看到是sdb4（”4”也有可能是别的数字,视情况而定）。
（2）.修改路径:
至此重启一下，回到U盘启动第一界面安装行处，然后按下Tab键(如果不行按’e’键)，将
vmlinuz initrd=initrd.img inst.stage2=hd:LABEL=RHEL-server-7.0-x86_64-LinuxProbe.Com.iso quiet改为(注意最后面的“quiet”切勿删除):
vmlinuz initrd=initrd.img inst.stage2=hd:/dev/sdb4 quiet然后，如果是Tab按下enter键进入安装界面，就可以顺利开始安装了。如果上面按的是’e’键就需要按ctrl+x进入安装界面，就可以顺利开始安装了

安装系统

如果重启重新进入安装选项，只需进入boot将hard设为第一即可。

1、首先会进入语言地区设置，按照习惯即可。后面进入安装界面。

2、选择软件-软件选择(S)，选择如图。

安装centos不进入桌面

1
2
3

yum  groupinstall  -y  "GNOME Desktop"
reboot now
Startx

3、选择安装位置(D)进入分区。
Linux分区：顺序（主要是boot swap /）

boot 分区（因为boot是引导启动的分区，所以分区的时候必须先分boot，通常设置大小为200M 空间足够300 - 500M）
*swap（缓存分区，通常设置大小为1G 通常是物理内存大小的2倍，比如你电脑是4G的物理内存，swap分区可以是8G）
/ （根分区）
/home分区（最大）安装完成后登录，通过命令：df -h 查看分区
也可以使用lsblk查看，更直观.
下面是基本的命名方案:
第一个软盘驱动器命名为/dev/fd0第二个软驱命名为/dev/fd1.
检测到的第一个硬盘名为/dev/sda。检测到的第二块硬盘命名为/dev/sdb，以此类推。
第一个SCSI CD-ROM命名为/dev/scd0，也称为/dev/sr0。
每个SCSI磁盘上的分区通过在磁盘名后面加上一个十进制数字表示：sda1和sda2表示系统中第一个SCSI磁盘驱动器的第一个和第二个分区。
假设你有一个有2个SCSI磁盘的系统，一个在SCSI地址2，另一个在SCSI地址4。第一个磁盘(地址2)被命名为sda，第二个磁盘被命名为sdb。如果sda驱动器上有3个分区，这些分区将被命名为sda1、sda2和sda3。这同样适用于sdb磁盘及其分区。
主分区：也叫引导分区，最多可能创建4个，当创建四个主分区时候，就无法再创建扩展分区了，当然也就没有逻辑分区了。主分区是独立的，对应磁盘上的第一个分区，“一般”就是C盘。在Windows系统把所有的主分区和逻辑分区都叫做“盘”或者“驱动器”，并且把所有的可存储介质都显示为操作系统的“盘”。因此，从“盘”的概念上无法区分主分区和逻辑分区。并且盘符可以在操作系统中修改，这就是要加上“一般”二字的原因。
扩展分区：除了主分区外，剩余的磁盘空间就是扩展分区了，扩展分区可以没有，最多1个。严格地讲它不是一个实际意义的分区，它仅仅是一个指向下一个分区的指针，这种指针结构将形成一个单向链表。这样在主引导扇区中除了主分区外，仅需要存储一个被称为扩展分区的分区数据，通过这个扩展分区的数据可以找到下一个分区(实际上也就是下一个逻辑磁盘)的起始位置，以此起始位置类推可以找到所有的分区。无论系统中建立多少个逻辑磁盘，在主引导扇区中通过一个扩展分区的参数就可以逐个找到每一个逻辑磁盘。
逻辑分区：在扩展分区上面，可以创建多个逻辑分区。逻辑分区相当于一块存储截止，和操作系统还有别的逻辑分区、主分区没有什么关系，是“独立的”。
硬盘的容量＝主分区的容量＋扩展分区的容量（硬盘=C盘+X）
扩展分区的容量＝各个逻辑分区的容量之和（X=D盘+E盘+F盘）

4、开始安装。此时会进入这个界面。同时完成1/2步骤，牢记密码。等待安装完成重启。

5、登录账户，设置网络即可。

/usr/bin/xauth: file /home/user/.Xauthority does not exist
错误原因:
是因为添加用户时没有授权对应的目录，仅仅执行了useradd user而没有授权对应的家目录

chown username:username -R /home/user_dir
例如本机：

然后继续登录
发现可以登录了
但是不显示用户名和路径
查看**/etc/passwd文件后发现，新建的用户未指定shell。我们只需将其指定为/bin/bash**即可。
Solution：在root用户下
usermod -s /bin/bash name
然后继续登录，就ok了。

总结：

1.首先，备份您的重要数据。
2.下载 CentOS 7 的安装镜像。您可以从 CentOs 官方网站上下载适用于您的计算机架构的镜像文件。
3.创建一个启动盘，将下载的镜像文件写入到 USB 设备或 DVD 上，以便用于安装过程，您可以使用工具如 Rufus 或 Etcher 来创建启动盘。UltraIso软件。
“Ultraiso根据镜像名为u盘命名盘符名称，默认的镜像名太长，将会导致Centos在装机时，无法定位到镜像，从而导致安装失败，出现drauct错误”。
4.将启动盘插入计算机并重启。在计算机启动时，进入 BIOS 设置界面，并将启动优先级设置为 USB 设备或 DVD.
5.选择安装选项。在 CentoS 安装程序启动后，按照屏幕上的指示选择合适的语言、键盘布局和时区。
6.分区和磁盘设置。根据需要配置磁盘分区和文件系统。您可以选择手动分区或使用默认选项。boot分区：是引导分区；作用：系统启动，在boot分区存放着grub，内核文件等，一般200M就够了。swap交换分区：内存扩展分区；交换分区给多大？ 一般做多：8G，16G，如果系统使用到了swap分区，就直接添加物理内存或排查一下服务器有没有被黑。/ 根：所有文件的根，绝对路径的开始标志。
7.设置网络和主机名。为您的计算机设置网络连接和主机名。
8.选择软件包组件。根据您的需求选择要安装的软件包组件。可以选择最小化安装或自定义安装。安装桌面勾选GUI和GNOME两个。关闭kdump。安全策略使用默认。
9.设置管理员密码。为 root 用户设置安全容码
10.开始安装。确认所有设置后，开始安装 CentOS 7 系统
11.完成安装。安装完成后，您可以重新启动计算机，并开始使用新安装的 CentOS 7 系统请注意，重新安装系统可能会导致数据丢失，请谨慎操作并确保提前进行数据备份

安装编译器

下载官方Get the Toolkits: https://www.intel.com/content/www/us/en/developer/tools/oneapi/toolkits.html
下载base和HPC，上传到服务器上。
赋予权限：
chmod +x l_*.sh
选择IDE即可，避免出错。需要预留30g左右的磁盘空间，否则会安装失败。sh ./l_Base*.sh sh ./l_HPC*.sh
默认安装位置：
如果使用sudo安装的话，默认路径是/opt/intel/oneapi/
如果使用普通用户权限安装的话，默认路径是~/intel/oneapi/
添加环境：vim ~/.bashrc,添加
如果使用sudo安装的话，是source /opt/intel/oneapi/setvars.sh >/dev/null
如果使用普通用户权限安装的话，默认路径是source ~/intel/oneapi/setvars.sh >/dev/null

保存退出

安装vasp

tar -xzvf vasp.6.3.0.tgz
cd vasp.6.3.0
cp arch/makefile.include.intel makefile.include
vim makefile.include
注释掉 MKLROOT    ?= /path/to/your/mkl/installation
保存退出
make all或者make std gam ncl

成功后添加环境变量：vim ~/.bashrc
export PATH=/home/hptem001/bin/vasp.6.3.0/bin:$PATH
保存退出

安装vaspkit1.4

tar -xzvf vaspkit.1.4.1.linux.x64.tar.gz
cd vaspkit.1.4.1
sh setup.sh
查看cat ~/.bashrc
有没有export PATH=/home/hptem001/bin/vaspkit.1.4.1/bin:${PATH}
或者直接终端输入：vaspkit看是否有输出。
安装成功。
vim ~/.vaspkit
更改赝势路径以自动输出POTCAR
更改其他可参考：
https://zhuanlan.zhihu.com/p/649806633 自动画图
http://vaspkit.cn/index.php/17.html 能带插值

保存退出。

0.3 VASP5单机安装

发表于 2023-05-08 更新于 2023-12-19 分类于软件安装

虚拟机VMware安装（Windows用户）

VMware Workstation Pro 17.0是一款专业功能最强大的虚拟机软件，用户可以在虚拟机同时运行各种操作系统，进行开发、测试、演示和部署软件，虚拟机中复制服务器、台式机和平板环境，每个虚拟机可分配多个处理器核心、主内存和显存。VMware Workstation Pro 版延续了VMware的传统，即提供专业技术人员每天在使用虚拟机时所依赖的领先功能和性能。借助对最新版本的Windows和客户机操作系统版本、最新的处理器和硬件的支持以及连接到VMware vSphere和vCloud Air的能力，让它成为提高工作效率、节省时间和征服云计算的完美工具。
具体安装方式微信公众号都有，不再赘述。
需要提醒的一点：可以把内存和磁盘空间（80G）留大一点用于安装系统和软件，以及计算数据保存。另外注意版本和操作系统版本是否对应。

https://www.haah.net/archives/1069.html

无论虚拟机还是服务器，都需要操作系统。Linux操作系统很多，以centos为例：

https://zhuanlan.zhihu.com/p/145102034

整个过程基本不会出错，按照步骤即可。centos是开源系统，不存在正版盗版的区别。

视频操作。

https://www.bilibili.com/video/BV1tP4y1x7Pz/?spm_id_from=333.337.search-card.all.click&vd_source=17084cc867ff3bb86398ff1a79b443f2

虚拟机快速安装VMware Tool

https://zhuanlan.zhihu.com/p/279241493#:~

Intel编译器

过程很简单，一路默认。

从intel官网下载：Intel® oneAPI Base Toolkit和Intel oneAPI HPC Toolkit。
加权限:

chmod +x l_HPCKit_p_2022.3.0.8751_offline.sh

chmod +x l_BaseKit_p_2022.3.0.8767_offline.sh.

然后运行：sudo ./l_HPCKit_p_2022.3.0.8751_offline.sh 和 sudo ./l_BaseKit_p_2022.3.0.8767_offline.sh. 按照提示一步步的完成即可。

sudo vim ~/.bashrc 中加入 source /opt/intel/oneapi/setvars.sh
执行which icc; which icpc; which ifort; which mpirun 和echo $MKLROOT 检验是否安装正确。

VASP编译

解压，进入文件夹。选择arch文件下的makefile.inclue.intel 复制到vasp6.3文件夹下，修改为makefile.inclue.

修改makefile.inclue中的内容为：

检验oneMKL with VASP是否连接成功，回到vasp6.3路径下面， ldd bin/vasp_std

vi ~/.bashrc 末尾加上：
export PATH=$PATH:/opt/vasp.6.3.0/bin

source ~/.bashrc

试运行
mpirun -np 4 vasp_std INCAR