5.1 分子动力学1

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

从头算分子动力学（ab initio molecular dYnamics, AIMD），又称为第一性原理分子动力学，为研究电子与原子核相互耦合系统的力学进动过程的理论方法。AIMD是在早期经验力场分子动力学基础上发展起来的理论方法。经验力场分子动力学中，针对原子核之间相互作用，均采用经验性的力场，这种力场往往是对相互作用，基于这些相互作用势场好处是计算量小，能够模拟的体系较大，模拟时长也比较长，统计数据也较多。

最初的AIMD为理想的Born-Oppenheimer分子动力学（BOMD），该理论将原子核-电子耦合体系的运动问题拆分为电子结构和分子动力学两部分，用密度泛函（DFT）等方法对特定原子核构型下的电子结构进行计算，得到K-S轨道及能量，在此基础上得到理想B-O势能面上原子核感受的势能和力，进而用经典力学来考察原子核在Born-Oppenheimer势能面上的运动。

为解决理想的BOMD所遇到计算瓶颈的问题，Reberto Car和Michele Parrinello于1985年提出了一种近似的BOMD方法- Car-Parrinello MD（简称CPMD）（Car, R., & Parrinello, M，1985）。该方法首次使用了一个扩展的Lagrangian量来描述电子-原子核耦合体系的动力学问题，Lagrangian量引入了虚拟电子质量u，基于这个Lagrangian量可给出运动方程，通过原子核的经典运动来近似更新电子波函数，而不再需要每一步通过矩阵对角化对电子结构进行SCF计算，简化了电子结构计算，大幅降低了计算成本，使得AIMD在技术上首次具备了可操作性，开创了近几十年AIMD模拟的时代。

第一性原理分子动力学（AIMD）结果分析

与经典分子动力学不同，第一性原理分子动力学不需要提供力场参数，只需要提供原子初始结构，就能根据电子波函数正交化产生的虚拟力，求解牛顿运动方程。在运行优化任务时，VASP生成的XDATCAR记录的是优化步骤的离子构型；在运行AIMD任务时，记录的就是运动轨迹。而现阶段读取XDATCAR轨迹分析性质的后处理软件并不多，能读取的兼容性也并不好。__VASPKIT0.72版本之后支持了将XDATCAR转换成通用的多帧PDB文件的功能（504）以便可视化并进行后处理分析__。但是并没有提供后处理分析接口，因此我们开发了一个Python脚本XDATCAR_toolkit.py，除了实现了选择一定范围内的帧数转换成PDB文件的功能，还可以提取分子动力学模拟过程中的能量，温度并做出变化趋势图。这对判断动力学是否平衡很有帮助。另外本脚本预留了接口，可以调用读取每一帧的晶格信息和原子坐标，以便进行后续扩展编程。此脚本需要安装了numpy包的python环境，以及matplotlib包以便于画图。

在得到通用轨迹PDB文件后，就可以利用现用的分子动力学后处理软件进行处理分析，比如VMD，MDtraj，MD Analysis， Pymol等。__本教程将演示通过VMD和MD Analysis软件包分析RDF（径向分布函数）和RMSD（均方根偏差）__，前者可以用来分析结构性质，后者对判断结构是否稳定以及模拟是否平衡很有帮助。

将XDATCAR转换成PDB文件
以VASP官网中单个水分子的AIMD模拟为例。模拟的输入文件如下，模拟的步长是0.5fs，模拟步数1000步，模拟时间500fs。脚本和测试例子可以在Github仓库(https://github.com/tamaswells/VASP_script/tree/master/XDATCAR_tookit)下载。

图1. 模拟的盒子
INCAR

{.line-numbers}

PREC = Normal    ! standard precision 
ENMAX = 400      ! cutoff should be set manually
ISMEAR = 0 ; SIGMA = 0.1

ISYM = 0         ! strongly recommened for MD
IBRION = 0       ! molecular dynamics
NSW = 1000       ! 1000 steps
POTIM = 0.5      ! timestep 0.5 fs

SMASS = -3       ! Nose Hoover thermostat
TEBEG =  2000 ; TEEND = 2000 ! temperature

NBANDS = 8

POSCAR

{.line-numbers}

H2O _2
0.52918   ! scaling parameter
 12 0 0
 0 12 0
 0 0 12
1 2
select
cart
      0.00     0.00     0.00 T T F
      1.10    -1.43     0.00 T T F
      1.10     1.43     0.00 T T F

KPOINTS

{.line-numbers}

Gamma-point only
 1        ! one k-point
rec       ! in units of the reciprocal lattice vector
 0 0 0 1  ! 3 coordinates and weight

模拟完成后将XDATCAR_toolkit.py上传到文件夹中（或者置于环境变量的路径文件夹中并赋予可执行权限即可直接调用命令XDATCAR_toolkit.py运行脚本），在shell环境中运行以下命令:

python XDATCAR_toolkit.py -p -t 0.5 --pbc
即可将XDATCAR的全部帧也就是0~499.5fs的轨迹转化成PDB格式。其中-p用于开启PDB转换功能，-t 0.5用于指定时间步为0.5fs，–pbc用于获取基于第一帧演变的连续轨迹。

Now reading vasp MD energies and temperature.
Now reading vasp XDATCAR.
Total frames 1000, NpT is False
Finish reading XDATCAR.
Selected time-range:0.0~499.5fs
[debug] Now entering function plotfigure.....

运行完成后，将会在文件夹内生成Temperature.dat，Energy.dat，ENERGY.png和XDATCAR.pdb四个文件，前面两个分别为温度和能量随着模拟时间的的变化数据，第三个是使用matplotlib绘制的趋势图（如下图），最后一个是转换得到的轨迹PDB文件，可以用于可视化轨迹，亦可用于后处理分析。

-b参数用于指定转换从哪一帧开始，-e参数用于指定转换到哪一帧结束。经刘锦程博士建议，增加一个–pbc的选项，用于处理周期性获取连续的轨迹。当分子穿过盒子边界时，记录真实的位置坐标（尽管它出了边界）而不是从盒子另一边穿入的ghost原子的坐标。这对于分析与时间相关性的量（比如RMSD）很有帮助。所谓连续指的是后面的轨迹都是从第一帧演变得到的真实坐标，但是并不能保证第一帧的分子是完整的，由于周期性的缘故，第一帧内摆放的分子可能分处于盒子两侧。李继存老师有篇博文(http://jerkwin.github.io/2016/05/31/GROMACS%E8%BD%A8%E8%BF%B9%E5%91%A8%E6%9C%9F%E6%80%A7%E8%BE%B9%E7%95%8C%E6%9D%A1%E4%BB%B6%E7%9A%84%E5%A4%84%E7%90%86/)讲的很明白，可以参考。如果发现第一帧内分子不完整，可以通过添加`-i 1参数将分子向第一个原子靠近平移以获得完整的分子。如果发现不理想，可以通过调整-i`的参数获得完整的分子。

图2. 温度和系统能量的变化趋势图
RDF径向分布函数分析

得到PDB文件后，可以使用VMD，MD Analysis等分子动力学后处理软件进行分析。

使用VMD分析工具分析

打开VMD，将PDB文件拖入显示窗口，在主菜单VMD Main中选择

Extensions-Analysis-Radial Pair Distribution Function g(r)，选择分析H(type H)在O(type O)周围的概率分布。值得注意的是分析RDF时,横坐标也就是max r不能超过盒子最小边长的一半，也就是得满足最小映像约定。如图4所示，在计算RDF时，如果max r的取值大于盒子最小边长的一半，就有可能重复算到一个粒子和它的映像粒子，这使得程序的周期性判断失准。将生成的dat文件的第一列和第二列作图即可得到RDF图。

图3. VMD中计算RDF

图4. 最小映像约定示意图

使用 MD Analysis分析 RDF

MD Analysis是一个成熟的分子动力学后处理软件，使用Python编写，开源。其教程不仅步骤详细还会给出背景理论知识。可以通过conda或者pip工具在线安装。

conda config --add channels conda-forge
conda install mdanalysis
#or
pip install --upgrade MDAnalysis

RDF分析的介绍和使用方法在网页(https://www.mdanalysis.org/docs/documentation_pages/analysis/rdf.html#radial-distribution-functions-mdanalysis-analysis-rdf)上查看。使用以下的脚本得到在O原子周围找到H原子的概率，并调用`matplotlib`绘制RDF图。在1.0 Å
处出现一个尖峰，也就是对应了O-H键的平衡键长（0.96$Å$）。

{.line-numbers}

import MDAnalysis
import MDAnalysis.analysis.rdf
import matplotlib.pyplot as plt

u = MDAnalysis.Universe('XDATCAR.pdb', permissive=True)
g1= u.select_atoms('type O')
g2= u.select_atoms('type H')
rdf = MDAnalysis.analysis.rdf.InterRDF(g1,g2,nbins=75, range=(0.0, min(u.dimensions[:3])/2.0))
           
rdf.run()

fig = plt.figure(figsize=(5,4))
ax = fig.add_subplot(111)
ax.plot(rdf.bins, rdf.rdf, 'k-',  label="rdf")

ax.legend(loc="best")
ax.set_xlabel(r"Distance ($\AA$)")
ax.set_ylabel(r"RDF")
fig.savefig("RDF_all.png")
#plt.show()

图5. RDF_O_H
RMSD均方根偏差分析

VMD分析RMSD

确保使用了-pbc参数以获取连续的轨迹，将生成的XDATCAR.pdb文件拖入显示窗口。如图6右所示，第一帧内水的三个原子不在同一个镜像内，分子不完整。在进行RMSD分析时，尽管轨迹是连续的，但是在对齐分子时就会出现问题。因此在本例中需要选择第一个原子作为中心将分子平移完整，在图6左中，分子已经在同一个镜像中了。

python XDATCAR_toolkit.py -p -t 0.5 --pbc -i 1

图6. 完整和不完整的水分子
将重新生成的PDB文件拖入显示窗口，在主菜单VMD Main中选择Extensions-Analysis-Analysis-RMSD Trajectory Tool，在计算RMSD前必须先做Align（对齐），这会使得每一帧结构进行平移、旋转来与参考帧的结构尽可能贴近，从而使得RMSD最小化。刘锦程提到研究生物法分子的RMSD时需要对齐操作，而研究小分子时不需要对齐分子。

图7. VMD中计算RMSD

把左上角文本框里的默认的Protein改成all（代表所有原子都纳入考虑），然后把noh复选框的勾去掉（否则将忽略氢原子）。然后点右上角的ALIGN按钮，此时所有帧的结构就已经对齐了。本例中演示以模拟的第一帧为参考，分析氧原子位置的均方根偏差。因此在Reference mol那里选top作为参考结构，左上角文本框由all改为type O（代表计算O原子的RMSD），然后勾上Plot复选框，最后点击RMSD按钮即可得到O原子的RMSD图。在File菜单栏可以选择导出dat数据。

图8. VMD中未对齐轨迹计算的RMSD

图9. VMD中对齐了轨迹后计算的RMSD

使用 MD Analysis分析 RMSD

RMSD分析的介绍和使用方法在网页（https://www.mdanalysis.org/docs/documentation_pages/analysis/rms.html?highlight=average）上查看。

使用以下的脚本可以分别得到所有原子，氢原子，氧原子的RMSD，并调用matplotlib绘制RMSD图。网页中有一段话（Note If you use trajectory data from simulations performed under periodic boundary conditions then you must make your molecules whole before performing RMSD calculations so that the centers of mass of the selected and reference structure are properly superimposed.）也就是在计算RMSD的时候选择的分子必须是完整的，不能分处于盒子的两边。这与我们之前的描述是一致的。MD Analysis默认对齐了分子。

使用以下脚本可以绘制对齐了轨迹后所有原子，氧原子和氢原子的RMSD。

{.line-numbers}

import MDAnalysis
import MDAnalysis.analysis.rms
import matplotlib.pyplot as plt

u = MDAnalysis.Universe('XDATCAR.pdb', permissive=True)
ref = MDAnalysis.Universe('XDATCAR.pdb', permissive=True)     # reference (with the default ref_frame=0)
ref.trajectory[0] #use first frame as reference
R = MDAnalysis.analysis.rms.RMSD(u, ref,
           select="all",         # superimpose on whole backbone of all atoms # align based on all atoms
           groupselections=["type H","type O"],
           filename="rmsd_all.dat",center=True)#,   # CORE
timestep=0.0005  #0.5fs from fs to ps as Reader has no dt information, set to 1.0 ps          
R.run()
rmsd = R.rmsd.T   # transpose makes it easier for plotting
time = rmsd[1]*timestep

fig = plt.figure(figsize=(5,4))
ax = fig.add_subplot(111)
ax.plot(time, rmsd[2], 'k-',  label="all")
ax.plot(time, rmsd[3], 'r--', label="type H")
ax.plot(time, rmsd[4], 'b--', label="type O")
ax.legend(loc="best")
ax.set_xlabel("time (ps)")
ax.set_ylabel(r"RMSD ($\AA$)")
fig.savefig("rmsd_md_analysis.png")

应用python工具pymatgen，

可以直接得到MSD，diffusion，donductivity。

得到XDATCAR后，需要首先配置pymatgen。
需要用到conda命令，windows：可以安装anaconda3，通过spyder启动。
进入anaconda prompt中启动中断，配置conda environment:

1 2	conda create –name my_pymatgen python activate my_pymatgen

linux环境下一个比较好的选择是安装miniconda3
wget https://repo.anaconda.com/minico … -Windows-x86_64.exe

bash Miniconda3-latest-linux-x86_64.sh
如果原来在bashrc中配置了anaconda，则注释掉，重启terminal。
与windows下相同，需要设置一个pymatgen环境。

1 2	conda create –name my_pymatgen source activate my_pymatgen

在my_pymatgen环境下安装 pymatgen：

conda install --channel conda-forge pymatgen

安装好之后运行建立如下test.py脚本, 可以得到结果：

{.line-numbers}

import os
from pymatgen.core.trajectory import Trajectory
from pymatgen.io.vasp.outputs import Xdatcar
from pymatgen import Structure
from pymatgen.analysis.diffusion_analyzer import DiffusionAnalyzer
import numpy as np
import pickle

# 这一步是读取 XDATCAR，得到一系列结构信息
traj = Trajectory.from_file('XDATCAR')

# 这一步是实例化 DiffusionAnalyzer 的类
# 并用 from_structures 方法初始化这个类； 300 是温度，1是POTIM 的time step，1000是间隔步数
# 间隔步数（step_skip）不太容易理解，但是根据官方教程(这里具体怎么回事我不太清楚，好像potim*step_skip需要小于1000，1000为NSW值，这是我没彻底弄清楚的地方):
# dt = timesteps * self.time_step * self.step_skip

diff = DiffusionAnalyzer.from_structures(traj,'Li',300,1,1000)

# 可以用内置的 plot_msd 方法画出 MSD 图像
# 有些终端不能显示图像，这时候可以调用 export_msdt() 方法，得到数据后再自己作图
diff.plot_msd()
diff.export_msdt("write_msd")
# 接下来直接得到 离子迁移率， 单位是 mS/cm，diffusity单位是 cm^2/S

C = diff.conductivity
D = diff.diffusivity
with open('result.dat','w') as f:
    f.write('# AIMD result for Li-ion\n')
    f.write('temp\conductivity\diffusivity\n')
    f.write('%d\t%.2f  %.10f' %(300,C,D))

在1.dat中是msd，conductivity和diffusivity会直接输出在result.dat中，模拟石墨烯表面Li的MD（excessive state）结果diffusivity为2*10-7 cm^2/S，我觉得算出diffusivity后自己求conductivity比较好，请问应该怎么求？与文件对比，基本吻合（J. Phys. Chem. Lett. 2010, 1, 1176–1180；https://pubs.acs.org/doi/10.1021/jp910134u）。
其它不对的地方，欢迎批评指正。
参考：
https://pymatgen.org/installation.html
https://www.bigbrosci.com/2020/09/08/A18/

VTST的脚本里有一个xdat2xyz.pl，可以直接把XDATCAR转化成movie.xyz文件。xdatcar可以通过ase转换成ms的xtd，也可以转ext-xyz，理论上支持xyz格式的程序也能打开。ase convert xdatcar xxx.xtd（注意会生成隐藏文件xxx.arc，如果要移动xtd，应该要使两个文件在相同的目录）。

事实上，有xtd格式的话，如果你有MS的版权，有些针对xtd格式的分析其实可以直接在MS里面做。arc文件跟着xtd一块拷贝到相同目录。

离子的电导率

Pymatgen 是 python materials genomics 的缩写，它是一款基于 python 的、开源的、强大的材料分析软件（https://pymatgen.org/）。

Pymatgen 包含一系列能够表示元素（Element）、位点（Site）、分子（Molecule）、和结构（Structure）的类（Class）。它具有为很多计算软件提供前处理和后处理的能力。这些计算软件包括VASP，ABINIT，exciting，FEFF，QCHEM，LAMMPS，ADF，AIIDA，ASE，Gaussian，Lobster，Phonopy，Shengbte，Pwscf，和Zeo++等等。它能实现科研狗的众多后处理需求，包括生成相图（Phase diagram）和布拜图（Pourbaix diagrams），分析态密度和能带等等。

Pymatgen 还提供了很多数据库（Materials Project REST API，Crystallography Open Database，and other external data sources）的接口，方便大家从数据库中查询结构和其他数据。

以下是Pymatgen官网提供的后处理的例子：

Top: (left) Phase and (right) Pourbaix diagram from the Materials API. Bottom left: Calculated bandstructure plot using pymatgen’s parsing and plotting utilities. Bottom right: Arrhenius plot using pymatgen’s DiffusionAnalyzer.

本文就介绍一下如何使用 Pymatgen 的 DiffusionAnalyzer 类去计算锂离子固态电解质中锂离子电导率。

计算离子电导率的理论与公式

目前，比较准确的计算离子电导率的方法是先用NVT系综第一性原理分子动力学（AIMD，ab initio molecular dynamics）模拟材料中离子在不同温度下的运动，然后计算出离子的平均（average）均方位移（MSD，mean square displacement），再计算出自扩散系数（D
，self-diffusion coefficient），最后求得离子在某温度下的电导率（
，conductivity）。

如何进行AIMD计算

AIMD计算通常非常耗时，所以，为了减少计算成本，我们可以适当放宽计算精度。如果用 VASP 进行计算，具体的，大家可以

采用较小的截断能。氧化物用 400 eV，硫化物用 280 eV，硒化物用 270 eV
采用Gamma点作为K点设置，并使用gam版本的 VASP 进行计算
采用单胞计算，如果材料的单胞包含比较多的原子
采用合适的步长，比如2 fs，即 POTIM = 2
后处理的基本公式
一旦AIMD计算完成，大家就可以着手计算离子电导率了。本文首先先介绍以下计算过程中使用的公式，方便有兴趣的同学自己开发脚本。

平均均方位移（average MSD）可以通过以下公式计算：

是第
个离子在
时刻的位移。

自扩散系数（
）可以通过以下公式计算：

是离子在材料中的扩散维度（一般地，
），
是离子扩散的时间。

最后，离子电导率（
）可以这样计算：

是材料中的离子密度，
是元电荷，
是离子的价态，
是玻尔兹曼常数，
是温度。

电导率计算的例子
现在我们通过一个 Li_Sn_S 材料的例子来详细了解一下整个计算和处理的过程。该材料的结构显示如下：

本例中采用单胞做计算，INCAR 设置如下：

[test@ln0%tianhe2 li_sn_s]$ vi INCAR

ISTART = 0
ICHARG = 2
IBRION = 0
ISIF = 2
NPAR = 8
NSW = 30000
TEBEG = 900 #还要设置成 1500K 等等
PREC = N
POTIM = 2
SMASS = 0.0
NELMIN = 4
LWAVE = F
LCHARG = F
IALGO = 48
LREAL = A

AIMD 计算结束之后会得到 XDATCAR 文件。很多时候，由于超算的时间限制，一个完整的AIMD计算需要提交两三次，从而产生两三个 XDATCAR 文件，这时，我们只要把它们按顺序通过 cat 命令合并在一起就行。例如我们有三个 XDATCAR 文件，分别命名成 XDATCAR01，XDATCAR02，和 XDATCAR03。

[test@ln0%tianhe2 li_sn_s]$ cat XDATCAR01 XDATCAR02 XDATCAR03 > XDATCAR
新得到的XDATCAR文件，注意删掉重复的与晶格信息相关的行，一般续算的次数也不多，在使用上面命令的时候，手动把XDATCAR02, XDATCAR03 中的删除即可。

Pymatgen 大显身手
安装pymatgen
首先让我们安装 pyamtgen，推荐大家参考官网，使用 anaconda 安装，否则会出现问题。安装好了anaconda之后，不管是 linux 还是 windows, 安装 pyamtgen 的指令是一样的。下面以吕梁天河超算为例：

[test@ln0%tianhe2 li_sn_s]$ conda install –channel conda-forge pymatgen
安装完成后，我们可以试着运行 python，导入 Pyamtgen 模块，如果像下面一样没有出错，就是安装成功了。

[test@ln0%tianhe2 li_sn_s]$ python
Python 3.7.3 (default, Mar 27 2019, 22:11:17) 
[GCC 7.3.0] :: Anaconda, Inc. on linux
Type "help", "copyright", "credits" or "license" for more information.
>>> import pymatgen
>>>

查看 DiffusionAnalyzer 的类
大家可以通过官方文档（https://pymatgen.org/pymatgen.analysis.diffusion_analyzer.html）查看接下来要使用的类，熟悉一下代码的用法。

class DiffusionAnalyzer(MSONable):
    def __init__(self, structure, displacements, specie, temperature,
                 time_step, step_skip, smoothed="max", min_obs=30,
                 avg_nsteps=1000, lattices=None):

这段代码显示，运行这个类需要一系列的输入信息，包括材料结构（structure），位移（displacements），要研究的离子（specie），温度（temperature）等等。

但是这个类提供了很多方法让大家可以通过读取 XDATCAR 或者 vasprun 文件的方式来实例化，例如

@classmethod
    def from_structures(cls, structures, specie, temperature,
                        time_step, step_skip, initial_disp=None,
                        initial_structure=None, **kwargs):
        """
        Convenient constructor that takes in a list of Structure objects to
        perform diffusion analysis.
        Args:
            structures ([Structure]): list of Structure objects (must be
                ordered in sequence of run). E.g., you may have performed
                sequential VASP runs to obtain sufficient statistics.
        ... ...
        """

好了，废话不多说，直接上代码，开始进行后处理。

代码示例
新建一个文件，名字为li_conductivity.py

‘’’
分析AIMD结果，计算MSD 和 conductivity
‘’’

import os
from pymatgen.core.trajectory import Trajectory
from pymatgen.io.vasp.outputs import Xdatcar
from pymatgen import Structure
from pymatgen.analysis.diffusion_analyzer import DiffusionAnalyzer
import numpy as np
import pickle

# 这一步是读取 XDATCAR，得到一系列结构信息
traj = Trajectory.from_file('XDATCAR')

# 这一步是实例化 DiffusionAnalyzer 的类
# 并用 from_structures 方法初始化这个类； 900 是温度，2 是POTIM 的值，1是间隔步数
# 间隔步数（step_skip）不太容易理解，但是根据官方教程:
# dt = timesteps * self.time_step * self.step_skip

diff = DiffusionAnalyzer.from_structures(traj,'Li',900,2,1)

# 可以用内置的 plot_msd 方法画出 MSD 图像
# 有些终端不能显示图像，这时候可以调用 export_msdt() 方法，得到数据后再自己作图
diff.plot_msd()

# 接下来直接得到 离子迁移率， 单位是 mS/cm
C = diff.conductivity

with open('result.dat','w') as f:
    f.write('# AIMD result for Li-ion\n')
    f.write('temp\tconductivity\n')
    f.write('%d\t%.2f\n' %(900,C))

在终端运行该文件

1	[test@ln0%tianhe2 li_sn_s]$ python li_conductivity.py

一段时间后就会得到MSD图像和离子电导率

[test@ln0%tianhe2 li_sn_s]$ vi result.dat

# AIMD result for Li-ion
temp	conductivity
900	884.05

可见，该材料在 900K 时的锂离子电导率为 884.05 mS/cm。

例子下载：
链接：https://pan.baidu.com/s/1WGzOVJBoe6Ym8mvR1uWanA
提取码：jhc5

思考
简短几行代码就可以计算出离子电导率，那么如何得出材料在300K下的电导率呢？
如何计算离子在材料中的迁移势垒？
如何可视化离子在材料中的扩散路径？

如何使用 Pymatgen 可视化离子的迁移概率密度。

先举个例子，

在“Design principles for solid-state lithium superionic conductors”一文中（Wang et al., Nature Materials 2015, 14 , 1026–1031. ），作者用Ab Initio Molecular Dynamic （AIMD）计算了Li 离子在Li$\mathrm{_1}\mathrm{_2}，
\mathrm{_7}
\mathrm{_3}
\mathrm{_1}$$\mathrm{_1}，
\mathrm{_2}，和
\mathrm{_4}
\mathrm{_4}$ 四种材料中的迁移概率密度（Probability Density），结果如下图所示：

从图中可以看出Li离子在图a所示材料中主要沿c轴方向的通道迁移，而且由于这个通道连通得比较好，Li离子的迁移势垒会比较低（0.220.25 eV）。
Li离子在图b所示的材料的迁移路径形成了一个三维网格，而且由于这个概率密度比图b中的概率密度分布得更加均匀，Li离子的迁移势垒更低（0.180.19 eV）。
图b所示的材料就完全不行了，因为Li离子的概率密度仅分布在特定的位点附近，说明离子不能有效地移动。
Li离子在图d所示材料中也存在迁移局域化的行为。
作者总结说 “A general principle for the design of Li-ion conductors with low activation energy can be distilled from the above findings: all of the sites within the diffusion network should be energetically close to equivalent, with large channels connecting them.”
那么我们如何在自己的计算中画出这样的图呢？Pymatgen 举手说，它可以帮忙！

但是在开始之前，我们要安装Pymatgen的插件：Pymatgen-diffusion（https://github.com/materialsvirtuallab/pymatgen-diffusion）。

安装 Pymatgen-diffusion
推荐大家使用最新版的Anaconda安装Pymatgen及其插件。点击上面的链接，进入官网后，点击最新版本链接，

我们可以下载.zip文件，

下载完成后，大家可以解压这个文件，得到 pymatgen-diffusion-2019.8.18文件夹。

我们把其中的 pymatgen_diffusion 文件夹放到 Anaconda的site-packages文件夹下，路径是 Windows 系统：……\Anaconda\Lib\site-packages；Linux系统：……/anaconda3/lib/pythonx.x/site-packages，就算安装好了。

接下来我们可以启动python，导入这个模块，如果不报错就没有问题了。

[test@ln0%tianhe2 li_sn_s]$ python
Python 3.8.3 (default, Jul  2 2020, 16:21:59) 
[GCC 7.3.0] :: Anaconda, Inc. on linux
Type "help", "copyright", "credits" or "license" for more information.
>>> import pymatgen_diffusion
>>>

学习用法
我们可以在其github网站上通过例子学习这个模块的用法。

点击打开 probbility_analysis.ipynb 文件。

其内容如下（有所删减）：如果不想看的话可直接查看开始作图部分

from pymatgen.analysis.diffusion_analyzer import DiffusionAnalyzer
from pymatgen_diffusion.aimd.pathway import ProbabilityDensityAnalysis
import json

#ProbabilityDensityAnalysis object
filename="/Users/iekhengchu/repos/pymatgen-diffusion/pymatgen_diffusion/aimd/tests/cNa3PS4_pda.json"

data = json.load(open("../pymatgen_diffusion/aimd/tests/cNa3PS4_pda.json", "r"))
diff_analyzer = DiffusionAnalyzer.from_dict(data) # 初始化DiffusionAnalyzer类

pda = ProbabilityDensityAnalysis.from_diffusion_analyzer(diff_analyzer, interval=0.5, 
                                                         species=("Na", "Li")) #可以指定离子
#Save probability distribution to a CHGCAR-like file
pda.to_chgcar(filename="CHGCAR_new2.vasp") #保存概率密度文件

开始作图
代码（test.py）如下：

from pymatgen_diffusion.aimd.pathway import ProbabilityDensityAnalysis
from pymatgen.core.trajectory import Trajectory
from pymatgen.io.vasp.outputs import Xdatcar
from pymatgen.analysis.diffusion_analyzer import DiffusionAnalyzer

traj = Trajectory.from_file('XDATCAR')
diff = DiffusionAnalyzer.from_structures(traj,'Li',900,2,1)
pda = ProbabilityDensityAnalysis.from_diffusion_analyzer(diff,interval=0.5,species=("Li"))
pda.to_chgcar(filename="pda.vasp") #保存概率密度文件

此处理过程大概耗时8分钟，因机器而异。

在VESTA中可视化

4.3 约束DFT

发表于 2023-05-11 更新于 2025-06-07 分类于 DFT学习

约束DFT

vasp中实现

https://wiki.kfki.hu/nano/Easy_manual_occupancy_of_Kohn-Sham_levels_with_FERWE_and_FERDO

The vasp wiki tells that one can set manual uccupancy with FERWE and FERDO:

One has to converge a WAVECAR before proceeding further.

Let say your system contain 2042 electrons. Assume the ground state is a closed shell singlet. It looks enough to do the calculation without any spin polarization with ISPIN=1 since both the ground and the excited state must have the same spin singlet wavefunction. But in order to move an electron from the HOMO to the LUMO one has to preconverge the ground state with ISPIN=2.

To do an excitation, one has to promote an electron from spatial orbital #1021 to #1022. To do so, one has to

ISTART = 1

ISTART=2 also might work, but ISTART =1 is not a simple vasp continuation. ISTART=1 will transform the old WAVECAR into the new job’s bases, or truncate the bands, or create new ones if the continuation job have more NBANDS setting than the original. Keep in mind NBANDS sighly varies with the number of MPI threads on which VASP is being executed, since NBANDS/THREADS must be an integer. VASP will automatically increase this value accordingly regardless your NBANDS is fixed.

ISMEAR = -2

ISPIN = 2 # 1020 1021 1022 1023 just a commend line

FERWE = 10201.0 10.0 11.0 10000.0 #The syntax is BANDS*OCC

FERDO = 10201.0 11.0 10.0 10000.0

The 1000*0.0 means zero electron for the next bands. Note this will do not override the NBANDS = 1150 setting is this will dont calculate 2022 bands, the bands after 1150 defined with FERWE and FERDO will be simply ignored.

The following SCF cycle alterations also required for HSE06 deltaSCF calculations:

ALGO = All TIME = 0.4 LDIAG = .FALSE. LSUBROT = .FALSE.

shf3
the shellpack have the option to copy the WAVECAR along the with the inputs. OTHER=”diamSiV.WAVECAR.gz” Just uncomment the OTHER in your *.guide file!

4.2 光学性质

发表于 2023-05-11 更新于 2025-07-29 分类于 DFT学习

介电常数

首先进行结构优化,参考官网： https://www.vasp.at/wiki/index.php/Dielectric_properties_of_SiC

静态介电常数：

方法1：DFPT方法;

{.line-numbers}

IBRION = 8
NSW=1
LEPSILON = .TRUE.
LRPA=.TRUE.#RPA考虑局域场效应，默认关闭，即IP近似
LPEAD=.TRUE.#有限差分法

方法2：Response to finite electric fields

{.line-numbers}

LCALCEPS = .TRUE.
IBRION=6
NFREE=2
ISMEAR =  0
SIGMA  =  0.01

在OUTCAR中抓取：

STATIC DIELECTRIC TENSOR (including local field effects in DFT)``` 电子贡献或者

```MACROSCOPIC STATIC DIELECTRIC TENSOR (including local field effects in RPA (Hartree))```
```MACROSCOPIC STATIC DIELECTRIC TENSOR IONIC CONTRIBUTION``` 离子贡献
```BORN EFFECTIVE CHARGES (in e, cummulative output)```
```ELASTIC MODULI IONIC CONTR (kBar)```
```PIEZOELECTRIC TENSOR IONIC CONTR  for field in x, y, z        (C/m^2)```
### 频率依赖介电函数：读取WAVECAR
IPA：
 ```javascript  {.line-numbers}
 ALGO = Exact
 NBANDS  = 64
 LOPTICS = .TRUE. 
 CSHIFT = 0.100
 NEDOS = 2000  
 #and you might try with the following
 #LPEAD = .TRUE.

frequency dependent IMAGINARY DIELECTRIC FUNCTION (independent particle, no local field effects)
and
frequency dependent REAL DIELECTRIC FUNCTION (independent particle, no local field effects)

Note：对于光学性质的计算，也就是计算材料的介电函数，需要足够多的空带和致密的K网格点，使其达到非常好的收敛状态，我们才可以得到合理的光学性质；因此通常计算中，一般设置NBANDS为默认值（可在自洽的OUTCAR中以关键字NBANDS查找到）的2-3倍，K网格点为自洽值或适当增加。
使用optics.sh脚本得到介电函数的实部和虚部

{.line-numbers}

# extract image and real parts of dielectric function from vasprun.xml
awk 'BEGIN{i=1} /<imag>/,\
               /<\/imag>/ \
                {a[i]=$2 ; b[i]=$3 ; c[i]=$4; d[i]=$5 ; e[i]=$6 ; f[i]=$7; g[i]=$8; i=i+1} \
    END{for (j=12;j<i-3;j++) print a[j],b[j],c[j],d[j],e[j],f[j],g[j]}' vasprun.xml > IMAG.in
#awk 'BEGIN{i=1} /imag/,\
#                /\/imag/ \
#                 {a[i]=$2 ; b[i]=$3 ; c[i]=$4; d[i]=$5 ; e[i]=$6 ; f[i]=$7; g[i]=$8; i=i+1} \
#     END{for (j=12;j<i-3;j++) print a[j],b[j],c[j],d[j],e[j],f[j],g[j]}' vasprun.xml > IMAG.in
awk 'BEGIN{i=1} /<real>/,\
               /<\/real>/ \
                {a[i]=$2 ; b[i]=$3 ; c[i]=$4; d[i]=$5 ; e[i]=$6 ; f[i]=$7; g[i]=$8; i=i+1} \
    END{for (j=12;j<i-3;j++) print a[j],b[j],c[j],d[j],e[j],f[j],g[j]}' vasprun.xml > REAL.in

vaspkit处理：
vaspkit-task 711

压电常数

上面计算完介电常数后，一般正交晶系。
grep "PIEZOELECTRIC TENSOR IONIC CONTR for field in x, y, z (C/m^2)" OUTCAR -A10
&
grep "PIEZOELECTRIC TENSOR for field in x, y, z (C/m^2)" OUTCAR -A10

极化强度

极化强度计算
首先选取极化相FE（PhaseB, P, 铁电相）和参考相NP（PhaseA, P=0, 中心对称相），分别以下列INCAR计算，

{.line-numbers}

 SYSTEM = FE
 ENCUT = 600  #铁电极化计算建议设置高精度
 EDIFF = 1E-6
 IBRION = -1
 POTIM = 0.3
 NSW = 0
 EDIFFG = -5E-3
 ISMEAR = 0
 SIGMA = 0.05
 PREC = Accurate
 ISIF = 2
 #NPAR = 4
LWAVE = .FALSE.
LCHARG = .FALSE.
LREAL = .FALSE. 
LDIPOL = .T.
DIPOL = 0.2 0.2 0.2 #不要设置在原子及迁移路径上，设置在真空层的一边/质心
LCALCPOL = .TRUE. #计算铁电极化开关 
IPEAD=4
LPEAD=.True.

收集结果：grep “Total electronic dipole momen” OUTCAR & grep “Ionic dipole moment” OUTCAR，离子与电子相加即可，然后用NP相-铁电相即可得到极化强度P，注意单位换算。
首先要理清数据单位。VASP 计算得到的 dipole moment 的单位是 e*Å，它与库仑之间换算为:
$1e = 1.602176634 * 10^{-19} C$
$1Å = 10^{-10} m$
三维体系的极化强度：极化值除以体积。单位为 $e/Å^2$；二维体系的极化强度：极化值除以面积。单位为 e/Å。

Note：

在铁电极化计算过程中经常会出现参考相为非半导体的情况，这种情况下可以：①镜像法：以FE以NP为参照中心，建立一个-FE相（PhaseB’, 即铁电极化方向相反）。然后按照上述方法计算极化，再以FE相的极化-（-FE相）的极化值，然后除以2即得到极化强度。可以理解为1-（-1）=2，2/2=1；②线性插值法：在FE和NP相中间插入一系列中间过渡相（0%（FE）,10%,20%,30%….100%(NP)），计算它们的极化，然后可以用Origin做拟合，100%时即为NP相的极化。（100%-0%）则是极化强度P。

线性插值法计算极化强度

原文链接:

https://www.cnblogs.com/ghzhan/articles/16305679.html

BaTiO3 是钙钛矿结构，它的铁电相结构和中心对称结构如图所示 0401 ，属于四方晶系。

先构造 BaTiO3 两种极化方向的晶格结构，并用 VASP 进行结构优化得到 CONTCAR；
将上一步优化后的两个结构分别放入创建好的 ini, fin 文件夹。利用 NEB 的 nebmake.pl 命令产生这两种极化方向的中间过渡结构 (vtst下载地址: Download — Transition State Tools for VASP (utexas.edu))，具体命令为：

nebmake.pl ini/CONTCAR fin/CONTCAR 32

这里的 ‘32’ 是表示产生中间过渡的 32 种结构。执行上述命令后，当前文件夹下会产生 00, 01, 02, …, 33 个文件夹，每个文件夹下有一个 POSCAR 文件。
对每个文件夹的结构进行一次 VASP 自洽运算，INCAR 文件里面需要额外设置 DIPOL 和 LCALPOL 参数。

DIPOL 参数可以选取任一坐标，但是保证同一体系采用相同值。

LCALPOL 参数是打开极化运算。

INCAR 文件:
1
2
3
4
5
6
7
8
9
10
11
12
SYSTEM=BaTiO3
ISTART =0
ICHARG =2
PREC =Accurate
ENCUT =520
EDIFF =0.1E-07
ISMEAR = 0
LWAVE=.FALSE.
LCHARG = .FALSE.
NELM = 200
DIPOL = 0.5 0.5 0.5
LCALCPOL=.TRUE
计算完成后运行命令：
1
2
3
grep 'dipole moment' OUTCAR|rev| awk '{printf ("%s ", $4)}'|rev|awk '{printf ("%f\n", $1+$2)}' >> ../dipole_c.dat

grep 'free energy TOTEN' OUTCAR | awk '{printf ("%s\n",$5)}' >> ../energy.dat
这里通过 grep 命令产生的 dipole_c.dat 文件记录的是沿着 c 方向的极化值，这是因为 BaTiO3 是沿着 c 方向极化的。对于具体的情况需要自行修改。
处理和分析数据。

首先要理清数据单位。VASP 计算得到的 dipole moment 的单位是 e*Å，它与库仑之间换算为:

$1e = 1.602176634 * 10^{-19} C$

$1Å = 10^{-10} m$

三维体系的极化强度：极化值除以体积。单位为 e/Å2；二维体系的极化强度：极化值除以面积。单位为 e/Å。

在这个例子中，BaTiO3 的体积为 64.3586 Å3。接下来，我们来处理数据。用energy.dat 文件 (单位为 eV)，这样，我们就得到每个中间结构 (image) 的总能：
用dipole_c.dat 文件 (单位: e*Å):不同 image 的极化值如图所示
可以发现，该极化强度不是连续的，这是和选取的原胞有关，需要考虑极化量子的影响。BaTiO3 沿着 c 方向极化，所以需要对该极化值加减整数倍的 c 方向的晶格常数。这样我们得到下图：

选取最靠近中间极化强度为 0 的那条曲线，即为 BaTiO3 的极化强度曲线：

除以体积并进行简单的单位换算后为:
另外，我们也可以绘制极化值 P 与能量 E 曲线
因此， BaTiO3 的极化强度大约为 0.248 C/m2，与文献中的实验结果 0.26 C/m2 吻合。（Physical Review, 99(4), 1161–1165, 1955 ）

拟合 Landau-Ginzburg 公式

$E= \sum_{i} A/2 * P_i^2 + B/4 * P_i^4 + C/6 * P_i^6 + D/2 * \sum_{i,j}(P_i-P_j)^2$

拟合该多项式曲线，得到

$E = 6062.434883706029P^6 + 2437.756598493882P^4 + -340.5806845162749*P^2 + 10.339600058315137$
其中参数 A, B, C 分别为

$A = -0.68116 eV *(m^2/C)^2
B = 9.75103 eV *(m^2/C)^4
C = 36.37461 eV *(m^2/C)^6$
0408

绘制能量曲线的脚本

{.line-numbers}

with open('energy.dat','r') as f:
    content = f.readlines()
res = [float(i.strip('\n')) for i in content]
minres = min(res)
res = [1e3*(i-minres) for i in res]
import matplotlib.pyplot as plt
plt.figure()
plt.plot(res,'b.-')
plt.xlabel('displacement')
plt.ylabel('Energy ($meV$)')
plt.show()

绘制dipole曲线，以及考虑极化量子的脚本:

{.line-numbers}

with open('dipole_c.dat','r') as f:
    content = f.readlines()
res = [float(i.strip('\n')) for i in content]

import matplotlib.pyplot as plt
plt.figure()
plt.xlabel('displacement')
plt.ylabel('Polarization ($e*{\AA}$)')
plt.xlim([0,32])
plt.plot(res,'r.')


### 考虑极化量子
c = 4.0335000000000001
tmp = []
N = 20
for j in range(N):
    start = -N/2+j
    tmp1 = [i+start*c for i in res]
    # print(tmp1[0])
    tmp.append(tmp1)
plt.figure()
for i in tmp:
    plt.plot(i,'.')

#plt.ylim([-2,2])
plt.xlabel('displacement')
plt.ylabel('Polarization ($e*{\AA}$)')
plt.xlim([0,32])

### 单位换算
P = [-1.1965400000000006, -1.1288099999999996, -1.060039999999999, -0.990219999999999, -0.9193599999999993, -0.8474299999999992, -0.7744599999999995, -0.7004900000000003, -0.6255299999999995, -0.5496499999999997, -0.4728999999999992, -0.3953799999999994, -0.31716999999999906, -0.2384000000000004, -0.1591799999999992,-0.11827, 0.0, 0.11827, 0.1591799999999992, 0.2384000000000004, 0.31716999999999906, 0.3953799999999994, 0.4728999999999992, 0.5496499999999997, 0.6255299999999995, 0.7004900000000003, 0.7744599999999995, 0.8474299999999992, 0.9193599999999993, 0.990219999999999, 1.060039999999999, 1.1288099999999996, 1.1965400000000006]
plt.figure()
P = [i*0.248945227832799346 for i in P]
plt.xlabel('displacement')
plt.ylabel('Polarization ($C/m^2$)')
plt.xlim([0,32])
plt.plot(P,'r.')
plt.show()

绘制极化值 P 与能量 E 曲线的脚本:

{.line-numbers}

E = [0.0035557599999975764, 0.0016252999999935014, 0.0004821100000000911, 0.0, 5.983999999870093e-05, 0.0005505999999968481, 0.0013690299999993272, 0.0024182999999950994, 0.003608839999998281, 0.004859239999994713, 0.006096989999996083, 0.007258279999994954, 0.008287759999994648, 0.009138139999997463, 0.00977216999999797, 0.010162939999993625, 0.010294929999993485, 0.010162939999993625, 0.00977216999999797, 0.009138139999997463, 0.008287759999994648, 0.007258279999994954, 0.006096989999996083, 0.004859239999994713, 0.003608839999998281, 0.0024182999999950994, 0.0013690299999993272, 0.0005505999999968481, 5.983999999870093e-05, 0.0, 0.0004821100000000911, 0.0016252999999935014, 0.0035557599999975764]
P = [-1.1965400000000006, -1.1288099999999996, -1.060039999999999, -0.990219999999999, -0.9193599999999993, -0.8474299999999992, -0.7744599999999995, -0.7004900000000003, -0.6255299999999995, -0.5496499999999997, -0.4728999999999992, -0.3953799999999994, -0.31716999999999906, -0.2384000000000004, -0.1591799999999992,-0.11827, 0.0, 0.11827, 0.1591799999999992, 0.2384000000000004, 0.31716999999999906, 0.3953799999999994, 0.4728999999999992, 0.5496499999999997, 0.6255299999999995, 0.7004900000000003, 0.7744599999999995, 0.8474299999999992, 0.9193599999999993, 0.990219999999999, 1.060039999999999, 1.1288099999999996, 1.1965400000000006]
import matplotlib.pyplot as plt
import numpy as np

E = [1e3*i for i in E]
P = [i*0.248945227832799346 for i in P]

plt.figure()
plt.plot(P,E,'b.-')
plt.xlabel('Polarization ($C/m^2$)')
plt.ylabel('Energy ($meV$)')

f1 = np.polyfit(P, E, 6)
poly = ''
for i in range(len(f1)):
    poly += '{}*x^{}  +  '.format(f1[i],len(f1)-i-1)
print(poly)

x = [ -0.35+0.7*i/99 for i in range(100)]
Eval = np.polyval(f1,x)
plt.figure()
plt.plot(x,Eval,'b')
plt.plot(P,E,'r*')
plt.xlabel('Polarization ($C/m^2$)')
plt.ylabel('Energy ($meV$)')
plt.show()

拉曼计算

VASP中计算拉曼谱主要使用vasp_raman.py这个脚本，见github：
https://github.com/raman-sc/VASP

使用raman-sc有几个注意的地方：https://zhuanlan.zhihu.com/p/354651066

作者这里提供了一个算例。供大家参考链接：https://pan.baidu.com/s/1dTCP_0bRcRo0uIy0rBGttQ
提取码：4vp8

这个代码已经年代久远，使用python2写的，对大部分新学vasp的同学来说可能装的anaconda3，既使用的是python3，与python2代码是不兼容的，作者这里对vasp_raman.py进行了改写，方便python3环境运行这段代码。下载见附件。
raman.sub既提交代码的脚本，需要根据自己的软件环境进行更改，我这里给出我的提交脚本供大家参考:
export VASP_RAMAN_RUN=’mpirun -np 64 /public/home/wangxiaohui/vasp/vasp.5.4.4/bin/vasp_std &> job.out ‘意思是运行vasp的命令，找一下vasp_std所在路径

python /public/home/wangxiaohui/WORKSPACE/Ti3C2Li2/raman/zhenraman/vasp_raman.py > vasp_raman.out这里是运行raman.py代码的意思，既raman.py代码所在路径

查看raman.py代码内容得知，POSCAR.phon和OUTCAR.phon是前一步计算的结果文件，不要被作者的例子给坑了。并且POSCAR.phon，OUTCAR.phon和INCAR等文件放在同一个文件夹下。

按照脚本VASP需要计算两步：

第一步，计算Gamma点声子。用有限位移方法或者DFPT方法。

第一步结束后，执行，

cp POSCAR POSCAR.phon

cp OUTCAR OUTCAR.phon

第二步，计算宏观介电常数的导数。用vasp_raman.py 脚本。

DFPT: LEPSILON=.TRUE.

频率依赖的介电矩阵: LOPTICS=.TRUE.

{.line-numbers}

#!/bin/bash
#PBS -l select=1:ncpus=32:mpiprocs=32
#PBS -l walltime=01:00:00
#PBS -q debug
#PBS -j oe
#PBS -N Example
#PBS -V

cd $PBS_O_WORKDIR

ulimit -s unlimited  # remove limit on stack size

export VASP_RAMAN_RUN='mpirun -n 12 /u/afonari/vasp.5.3.2/vasp.5.3/vasp &> job.out'   #a按照自己路径更改
export VASP_RAMAN_PARAMS='01_10_2_0.01'  #b按照自己材料更改

python27 vasp_raman.py > vasp_raman.out

a. VASP_RAMAN_RUN设置调用VASP的命令

b. VASP_RAMAN_PARAMS=01_10_2_0.01

FIRST_MODE - 整数, 第一个计算的振动模
LAST-MODE - 整数, 最后一个计算的振动模
NDERIV - 整数, 有限位移的策略选择，目前只有2这一种选择
STEPSIZE - 浮点数, 有限位移步长，单位是Angstroms，一般选默认就可以

注意：以上这种设置格式允许我们并行计算大体系拉曼，我们可以在计算大体系的时候将Γ点的振动频率分为很多段分别计算，这样会大大提高我们计算大体系的效率。

在多段计算的时候。具体细节上要把每一段计算的OUTCAR.000n.-1.out和OUTCAR.000n.-1.out文件拷贝到同一个文件夹下，然后注释掉vasp_raman.py脚本中的运行VASP的命令。运行python2.7 vasp_raman.py > vasp_raman.out 即可得到整段的数据。

vasp_raman.py文件中执行VASP的行为331行：

1
2
3

else: # run VASP here
                  print "[__main__]: Running VASP..."
                  os.system(VASP_RAMAN_RUN) #注释掉这一行

加展宽

使用脚本broaden.py给出展宽谱，可以选择洛伦兹线型或者高斯线型。

$ python2 broaden.py result.txt
生成result.txt-broaden.dat文件，画图即可。

需要注意的是，broaden.py的第23和24行分别代表频率和强度所在列。其设置文件中对应的列。

cm1 = hw[:,1]
int1 = hw[:,2]

红外谱计算

红外谱主要是因为振动模在振动前后偶极矩发生变化，从而产生后外吸收。在振动过程中引起偶极矩改变的振动都具有红外活性。
根据上述原理，我们首先需要在计算过程中计算Gamma点的振动模，其次，我们需要计算伯恩有效电荷。

有限位移方法INCAR关键参数：

IBRION = 5
NWRITE = 3
POTIM = 0.015
LEPSILON = .TRUE
DFPT方法INCAR关键参数：

IBRION = 7
NWRITE = 3
NSW = 1
LEPSILON = .TRUE.
NELMIN=6

另外，建议进行k点测试直到伯恩有效常数收敛。

后处理

计算得到OUTCAR，然后文件夹下运行脚本。运行结果如下：

sh IR.sh

上面的振动模和活性是归一化以后的结果，也就是运行结果intensities/results/results.txt的内容。文件intensities/results/exact.res.txt是归一化之前的内容。

展宽

使用脚本broaden.py给出展宽谱，可以选择洛伦兹线型或者高斯线型。

$ python2 broaden.py result.txt
生成result.txt-broaden.dat文件，画图即可。

4.1 力学性质

发表于 2023-05-10 更新于 2023-12-17 分类于 DFT学习

弹性常数

基本概念

弹性常数描述了晶体对外加应变的响应的刚度。在材料的线性变形范围内（应变较小的情况下），体系的应力与应变满足胡克定律。也就是说，对于足够的小的变形，应力与应变成正比，即应力分量(S)是应变分量(E)的线性函数，三维材料的弹性刚度常数矩阵是6×6的。公式中Cij就是我们通常所说的弹性常数。因为刚度矩阵是对称矩阵，因此，弹性常数的独立张量元数目至多只有21个。对不同的晶系的晶体，因为对称性的关系，其独立的弹性常数是确定的。因此，晶系的对称性越高，独立的张量元数目越少。
立方晶系 ——只有3个独立矩阵元（C11，C12，C44）
六角晶系 ——有5个独立矩阵元（C11，C12，C13，C33，C44）
三角晶系
a) 32，3m，-32/m——有6个独立矩阵元（C11，C12，C13，C14，C33，C44）
b) 3，-3,——有8个独立矩阵元（C11，C12，C13，C14，C15，C33，C44，C45）
四方晶系
a) 422，4mm，-42m，4/mmm——有6个独立矩阵元（C11，C12，C13，C33，C44，C66）
b) 4，-4，4/m——有7个独立矩阵元（C11，C12，C13，C16，C33，C44，C66）
正交晶系 ——有9个独立矩阵元（C11，C12，C13，C22，C23，C33，C44，C55，C66）
单斜晶系 ——有13个独立矩阵元
三斜晶系 ——有21个独立矩阵元

计算过程

采用第一性原理计算弹性常数有两种方法，第一种方法是应力-应变方法， 即通过给结构施加不同应变，分别计算出所对应的应力大小，然后利用公式（1）拟合得到一次项系数，从而得到。第二种方法是能量-应变方法， 即通过给结构施加不同应变后计算出体系总能相对于基态能量变化（应变能），再利用公式（2）进行二次多项式拟合，其中二次多项式系数是晶体的某个弹性常数或者弹性常数组合。第一种方法优点是每进行一次计算可一次得到六个独立分量，缺点是为了得到准确的应力大小，必须选择更高的截断能和更密集的K格点。在同样的精度下，能量-应变法相比应力-应变方法要求的截断能和K点数目相对较少，缺点是计算应变数目有所增加。VASPKIT中弹性常数计算基于能量-应变法，目前不支持三斜晶系。
一般计算弹性常数我们可以采用应变应力或者应变能量关系，进行拟合。或者可以直接通过VASP直接计算。

方法一：

首先介绍直接用VASP在输入文件中添加参数进行计算的方法。

当然首先要进行结构优化，

在INCAR中添加关键参数，进行弹性常数的计算

{.line-numbers}

NSW=1
EDIFFG=-1e-3
ISIF=3
IBRION=6 #计算弹性常数
POTIM=0.015
NFREE = 4

计算完成后可以查看OUTCAR得到弹性刚度矩阵。或者可以利用VASPKIT直接得到弹性常数。
请注意OUTCAR中的刚度矩阵是未经过处理的刚度矩阵，二维材料需要乘以C轴长度的0.01，还需要了解Voigt简标，11>1 22>2 33>3 23>4 13>5 12>6。然后我们通过Origin等工具绘图即可得到材料的杨氏模量，泊松比等力学性质。

方法二：

VASPKIT和VASP计算晶体的弹性常数具体计算步骤分为：

准备优化彻底的POSCAR文件，注意通常采用标准的惯用原胞计算弹性常数，如果不确信POSCAR文件中是否是标准的惯用原胞，可以用vaspkit-603/604生成标准结构；
运行vaspkit 选择102生成KPOINTS，由于计算弹性常数对K-mesh要求很高，因此对于半导体（金属体）体系，生成K点的精度应不小于 0.03(0.02)×2π {\AA}^{-1}

INCAR参考设置如下；

{.line-numbers}

Global Parameters
 ISTART =  0        
 LREAL  =  F      
 PREC   =  High  （截断能设置默认值1.5-2倍）
 LWAVE  = F       
 LCHARG = F     
 ADDGRID= .TRUE.    
Electronic Relaxation
 ISMEAR =  0          
 SIGMA  =  0.05       
 NELM   =  40           
 NELMIN =  4           
 EDIFF  =  1E-08       
Ionic Relaxation
 NELMIN =  6          
 NSW    =  100         
 IBRION =  2           
 ISIF   =  2    （切记选择2,如果选择3会把施加应变后原胞重新优化成平衡原胞）
 EDIFFG = -1E-02

准备VPKIT.in文件并设置第一行为1（预处理）；运行vaspkit并选择201, 将生成用于计算弹性常数文件；

{.line-numbers}

1                    ! 设置1将产生计算弹性常数的输入文件，2则计算弹性常数
3D                   ! 2D为二维体系，3D为三维体系
7                    ! 7个应变
-0.015 -0.010 -0.005 0.000 0.005 0.010 0.015  ! 应变变化范围

批量提交vasp作业；
总计有9个大文件夹，分别以该空间群所需要的独立弹性常数来命名；在每个大文件夹内则是以应变强度命名的实际计算执行的文件夹，文件夹的数量可根据预处理的VPKIT.in文件的最后一行的应变设置来更改。在每一个应变文件夹内则是施加了应变的POSCAR和事先准备好的INCAR、KPOINTS和POTCAR文件的链接，可直接进行vasp计算。
这里通过一个简单的两重for循环命令依次序完成所需要的计算任务。
脚本 :

{.line-numbers}

#!/bin/bash
path=`pwd`
#记录当前文件夹位置并定义到变量 path中，称为初始目录
pwd
for cij in `ls -F | grep /$`
#ls -F| grep /$ 命令为返回当前文件夹中所有子目录
#将所有子目录添加到变量cij中，并准备进行for循环
do
cd ${path}/$cij
#进入cij变量所代表的子目录中
pwd
for s in strain_*
#将cij代表的子目录中所有strain文件夹添加到变量s中
#此部添加变量s的方式可更改为添加变量cij的方式，但不建议这么做
do
cd ${path}/$cij/$s
#进入初始目录下变量cij和变量s组成的目录，由于for循环的设置，将遍历所有可计算的strain文件夹
pwd
 mpirun -np 16 vasp_std
 #执行vasp计算
cd ${path}
#返回初始目录
pwd
done
done

再次修改VPKIT.in文件中第一行为2（后处理），然后再次运行vaspkit并选择201，得到以下结果；
如果在计算弹性常数时希望强制固定某个体系的空间群，可在SYMMETRY.in文件第二行设置空间群，这样晶系将通过指定空间群判断并计算相应的独立弹性常数，这对于掺杂或者合金体系比较有用。

根据弹性常数计算和三维可视化材料力学量

首先利用VASPKIT-204命令得到体弹性模量、杨氏模量、剪切模量及泊松比随角度的依赖关系，具体算法可参考文献【J. Phys. Condens. Matter, 28, 275201 (2016)和Comput. Phys. Commun. 181, 2102–2115 (2010)】。接下来我们以某单斜体系的弹性常数为例来演示如何进行材料力学量三维可视化。

第一步新建ELASTIC_TENSOR.in并输入6x6弹性常数矩阵，注意第一行是注释行不可省略。如果是二维体系，则输入文件为ELASTIC_TENSOR_2D.in，注意二维体系弹性常数矩阵大小为3x3。

{.line-numbers}

# comment line (in GPa)
   228.38   85.741   81.503        0   -0.737        0
   85.741   217.47   94.201        0  -20.213        0
   81.503   94.201   178.81        0   -9.472        0
        0        0        0   35.094        0  -17.851
   -0.737  -20.213   -9.472        0   37.778        0
        0        0        0  -17.851        0   42.708

第二步运行VASPKIT 204命令得到MECHANICS_3D.dat
打开MECHANICS_3D.dat可看到第一行给出了每一列数据所代表的物理量，第二行给出了球坐标系下分别划分仰角θ∈[0,180]和方位角ϕ∈[0,360]的格点大小。
第三步从vaspkit/examples/angular_dependent_mechanics (VASPKIT ver. >= 1.3.2)文件夹中拷贝mechanics_3d_plot_matlab.m到当前目录，调用Matlab软件运行该脚本，得到以下信息.
如果我们想可视化杨氏模量，则输入2回车即可得到。
另外，以下两个软件也可以实现材料力学量三维可视化。
Elate
ElasticPOST

3.4 非谐声子

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

介绍

ALAMODE中非谐声子计算部分采用了自洽声子计算方法。

具体教程：https://alamode.readthedocs.io/en/latest/anphondir/formalism_anphon.html#self-consistent-phonon-scph-calculation

https://www.bilibili.com/video/BV1mv411M7W7/?spm_id_from=333.999.0.0&vd_source=17084cc867ff3bb86398ff1a79b443f2

TDEP 是其采用技术Temperature dependent effective potential的缩写。
下面是开发者的网站，但是并不能下载到源代码。http://ollehellman.github.io/page/index.html好消息是第一性原理计算软件ABINIT中集成了TDEP方法，官网user guide>a-TDEP。需要研究的可以自行研究。ABINIT官网：https://www.abinit.org/

Dynaphopy 是采用了一种基于AIMD（第一性原理分子动力学）的杂化方法的python程序包：http://abelcarreras.github.io/DynaPhoPy/index.html
文中给出了参考文献，需要的自行下载。

ALAMODE

安装

https://alamode.readthedocs.io/en/latest/index.html#

ALAMODE计算Si简谐声子谱

1、利用ALM计算displacement patterns

优化后的超胞结构重命名POSCAR.orig，准备输入文件si_alm.in：

&general
  PREFIX = si222
  MODE = suggest
  NAT = 64; NKD = 1
  KD = Si
/

&interaction
  NORDER = 1  # 1: harmonic, 2: cubic, ..
/

&cell
  20.406 # factor in Bohr unit
  1.0 0.0 0.0 # a1
  0.0 1.0 0.0 # a2
  0.0 0.0 1.0 # a3
/

&cutoff 
  Si-Si None
/

&position
  1 0.0000000000000000 0.0000000000000000 0.0000000000000000   
  1 0.0000000000000000 0.0000000000000000 0.5000000000000000
  1 0.0000000000000000 0.2500000000000000 0.2500000000000000
  1 0.0000000000000000 0.2500000000000000 0.7500000000000000
  1 0.0000000000000000 0.5000000000000000 0.0000000000000000
  1 0.0000000000000000 0.5000000000000000 0.5000000000000000
  1 0.0000000000000000 0.7500000000000000 0.2500000000000000

晶格常数替代为自己的超胞结构。

运行ALM: alm si_alm.in > si_alm.log1，将会生成si222.pattern_HARMONIC文件，包含坐标信息。对于Si只有一个displacement patterns。

2、计算原子间力

在此之前要确定位移的大小，一般简谐0.01-0.04，0.01大多可以。运行python displace.py --VASP=POSCAR.orig --mag=0.01 -pf si222.pattern_HARMONIC

之后对每个结构提交计算，脚本：

#!/bin/bash

# Assume we have 20 displaced configurations for QE [disp01.pw.in,..., disp20.pw.in].

for ((i=1;i<=20;i++))
do
    num=`echo $i | awk '{printf("%02d",$1)}'`
    mkdir ${num}
    cd ${num}
    cp ../disp${num}.pw.in .
    pw.x < disp${num}.pw.in > disp${num}.pw.out
    cd ../
done

计算完成后，提取信息：python extract.py --VASP=POSCAR.orig vasprun*.xml > DFSET_harmonic，生成DFSET_harmonic文件。

3、拟合力常数

修改1中的 si_alm.in ：

&general
  PREFIX = si222
  MODE = optimize   # <-- here
  NAT = 64; NKD = 1
  KD = Si
/

&optimize
  DFSET = DFSET_harmonic
/

运行alm si_alm.in > si_alm.log2，生成si222.fcs and si222.xml

4、计算声子谱和DOS

准备si_phband.in文件，结构中用晶胞结构，注意&cell中单位是Bohr，具体可设置如下

&general
  PREFIX = si
  MODE   = phonons
  FCSXML = si.xml

  NKD = 1; KD = Si
  MASS = 28.085
/

&cell
  10.184247326051628 
  0.0000000000000000    0.5071343999939496    0.5071343999939496
  0.5071343999939496    0.0000000000000000    0.5071343999939496
  0.5071343999939496    0.5071343999939496    0.0000000000000000
/

&kpoint
  1   # KPMODE = 1: line mode
  G 0.0 0.0 0.0 X 0.5 0.5 0.0 51
  X 0.5 0.5 1.0 G 0.0 0.0 0.0 51
  G 0.0 0.0 0.0 L 0.5 0.5 0.5 51
/

执行anphon si_phband.in > si_phband.log会产生si.bands文件，可以用plotband.py脚本画图。执行python plotband.py si.bands得到

准备si_phdos.in文件，结构中用晶胞结构，注意&cell中单位是Bohr，具体可设置如下

&general
  PREFIX = si
  MODE   = phonons
  FCSXML = si.xml

  NKD = 1; KD = Si
  MASS = 28.085
/

&cell
  10.184247326051628 
  0.0000000000000000    0.5071343999939496    0.5071343999939496
  0.5071343999939496    0.0000000000000000    0.5071343999939496
  0.5071343999939496    0.5071343999939496    0.0000000000000000
/

&kpoint
  2   # KPMODE = 2: uniform mesh mode
  20 20 20
/

执行anphon si_phdos.in > si_phdos.log会产生si.dos文件，可以用plotdos.py脚本画图。执行python plotdos.py –emax 550 –nokey si.dos得到

5. 估算对热导率的立方IFCs(非简谐IFCs)

这里计算非简谐IFCs，我们把si_alm1.in复制为si_alm3.in，将NORDER = 1改为NORDER = 2.我们还要注意这里考虑计算三阶力常数，所以&cutoff下的参数我们要再指定个截止半径，一般要求略大于第二近邻的距离（可以从si_alm1.log查看，如下图）。当然你接着用None也行，这样可能大大增加参数的数量，从而增加计算成本。

运行alm si_alm3.in > si_alm3.log，可以看到多了si_cubic.相关文件。运行python displace.py –OpenMX=Si.dat –mag=0.04 -pf si_cubic.pattern_ANHARM3，可见产生了disp01.dat到disp11.dat11个OpenMX输入文件，我们准备shell脚本批量提交，shell脚本命名为pbs.sh，具体为

#!/bin/bash

for ((i=1;i<=11;i++))
do
    num=`echo $i | awk '{printf("%02d",$1)}'`
    mkdir ${num}
    cd ${num}
    cp ../disp${num}.dat .
    cp disp${num}.dat Si.dat
    cp ../openmx.pbs .
    qsub openmx.pbs
    cd ../
done

总结

简谐计算步骤

输入文件：优化后的超胞结构重命名POSCAR.orig，准备输入文件si_alm.in（其中mode：suggest）；输出：si222.pattern_HARMONIC文件，包含坐标信息。
displace.py对上面的.pattern文件处理，得到随机位移结构；对随机位移机构进行力的计算，生成不同vasprun.xml；extract.py处理上面vasprun.xml，生成DFSET_harmonic文件。
力常数拟合输入文件：DFSET_harmonic文件，修改1中的 si_alm.in，其中mode改成opt, 添加下面；生成si222.fcs and si222.xml。
1
2
3
&optimize
DFSET = DFSET_harmonic
/
后处理输入文件：准备phband.in文件，修改MODE = phonons
FCSXML = si.xml；修改kmesh格式可得到声子谱、态密度、热性质等。

非谐计算

准备简谐计算得到的si222.fcs and si222.xml，AIMD得到的vasprun.xml，超胞结构。
displace.py 对AIMD的vasprun.xml处理，得到随机结构；对随机位移机构进行力的计算，生成不同vasprun.xml；extract.py处理上面vasprun.xml，生成DFSET_AIMD_random文件。
Cross validation (CV)：alm_cv.in文件中添加简谐力常数 FC2XML = ../harmonic.xml，添加高阶力常数部分标签和cv相关标签，mode = cv；得到*cvset文件；python3 cvscore.py *cvset > si222.cvscore得到Minimum cvscore
拟合力常数输入文件：alm_opt.in中mode=opt；修改cv相关标签；生成高阶力常数cv.fcs and cv.xml.在这一步alm 的输入文件&general 部分加入 FC3_SHENGBTE=1、 FC4_SHENGBTE=1可输出ShengBTE格式的输入文件。

SCPH计算：scph.in文件中修改MODE=SCPH和温度标签；添加&scph模块，见上面例子。anphon scph.in > scph.log提交，查看收敛：grep “conv” scph.log。生成输出文件：.scph_dymat、.scph_dfc2、.scph_bands。对.scph_dfc2处理得到有效二阶力常数文件，用于后续计算。

>dfc2
DFC2--a generator of renormalized harmonic FCs from SCPH outputs.
XML file containing original Fc2 : STo222.xml 
Output xml filename with anharmonic correction : ST0222_SCPH2-2_300K.xml
FC2 correction file fromSCPH calculation :STo_scph2-2.scph_dfc2
Target temperature : 300


New XML file STo222_SCPH2-2_300K.xml was created successfully.

Alamode计算高温声子色散

Anharmonic Phonons (self-consistent phonon theory)
主要使用alamode软件。版本：alamode-1.1.0

输入文件

vasprun.xml #分子动力学输出文件
Dfile  #包含原子的力和位移
extract_disp.py #选取结构，并计算位移
extract_force.py #选取结构，并获得受力
PPOSCAR #初始超胞
alm.in #提取力常数的输入文件
anphono.in #计算声子谱的输入文件

从分子动力学数据中选取一定数量的结构，并提取它们的力和计算与初始状态相比的位移，输出到Dfile文件当中

python extract_disp.py
#运行后根据提示输入开始的步数和间隔步数
#例如总共1000步，从第一步开始提取，每间隔五步提取一次
#得到disp.dat

python extract_force.py

方法同上得到force.dat
将disp.dat和force.dat文件按如下形式合并到一起

disp1x  disp1y  disp1z   force1x  force1y  force1z
       ......
       ......
dispNx  dispNy  dispNz   forceNx  forceNy  forceNz

获取力常数

提取四阶力常数的alm.in文件如下（使用LASSO方法）：

二阶力常数可以通过有限位移法提取(详见官网手册)

&general
  PREFIX = ThO2_ENET
  MODE = optimize
  NAT = 81; NKD = 2
  KD = Th O
#  PRINTSYM=1
/

&interaction
  NORDER = 5  # 1: harmonic, 2: cubic, ..
  NBODY = 2 3 3 2 2
/

&cell
  1 # factor in Bohr unit
  0.0000000000000   15.919498944104   15.919498944104
  15.919498944104   0.0000000000000   15.919498944104
  15.919498944104   15.919498944104   0.0000000000000
/

&optimize
  LMODEL= enet
  FC2XML = ThO2_2nd.xml
  DFSET=ThO2_dffile

  CV = 0
  L1_RATIO = 1.0
  L1_ALPHA = 1.0e-06
  CV_MINALPHA = 1.0e-6
  CV_MAXALPHA = 0.02
  CV_NALPHA = 100

  STANDARDIZE = 1
  CONV_TOL = 1.0e-9
/

&cutoff 
  *-* None None 12.0 12.0 12.0
/

&position
1  0.0000000000000000  0.0000000000000000  0.0000000000000000
1  0.3333333333333333  0.0000000000000000  0.0000000000000000
..................

注意晶格常数单位的转换 $1 Bohr= 0.529177208 Angstrom$

在终端运行：

bash
alm alm.in > alm.log
得到Nb.xml

Anharmonic Phonon

anphono.in输入文件如下：

&general
 PREFIX =scph
 MODE = SCPH
 FCSXML =Nb.xml #forceconstant file

 NKD = 1; KD = Nb
 MASS = 92.90638

 #NONANALYTIC = 0; BORNINFO = PbTe.born

 TMIN = 0; TMAX = 300; DT = 100   #Temperature and step
/

&scph
 KMESH_SCPH        = 6 6 6
 KMESH_INTERPOLATE = 6 6 6

 SELF_OFFDIAG = 0
 RESTART_SCPH = 0

 MIXALPHA = 0.1
 MAXITER = 500
 TOL_SCPH = 1.0e-10
/

&cell   ## primitive cell is best
  6.2507550151332
  -0.5  0.5  0.5
   0.5 -0.5  0.5
   0.5  0.5 -0.5
/

&kpoint
  1  # KPMODE = 1: get phonon band  2: get others properties
  G 0.0 0.0 0.0 H 0.5 -0.5 0.5 201
  H 0.5 -0.5 0.5 P 0.25 0.25 0.25 201
  P 0.25 0.25 0.25 G 0.0 0.0 0.0 201
  G 0.0 0.0 0.0  N 0.0 0.0 0.5 201
/

在终端运行：

anphon anphono.in > anphono.log
得到scph.scph_bands，里面有不同温度下的声子谱

进入安装目录下的tools文件夹：make Makefile，编译df2c。

[scf1355@ln151%nc-l ah]$ ~/soft/alamode/tools/dfc2 
 DFC2 -- a generator of renormalized harmonic FCs from SCPH outputs.
 XML file containing original FC2 : crte2-442.xml
 Output xml filename with anharmonic correction : crte2-new.xml
 FC2 correction file from SCPH calculation : crte2-442.scph_dfc2
 Target temperature : 300

输出crte2-new.xml.利用https://github.com/ttadano/alamode/discussions/192 中提到的脚本转换成phonopy格式。此时的力常数文件由于对称性数量不匹配，利用phonopy读取并写出FULL_FORCE_CONSTANTS。

Alamode_to_ShengBTE

Alamode计算热导率步骤

所有的计算都是基于密度泛函软件 VASP 来进行的。
首先利用第一性原理的密度泛函对晶体结构进行优化，弛豫出稳定的晶格结构，然后利用密度泛函微扰理论计算出波恩有效电荷和静态介电张量，以便用于后续求解动力学矩阵的非解析部分。
随后为了获得更加精准的二阶力常数，通过差分位移的方法生成微扰结构，后利用静态密度泛函理论计算出力和位移信息，再利用最小二乘法得到所需要的二阶力常数。
对于高阶力常数，本文利用分子动力学模拟得到 80 个准随机结构，对每个准随机结构的原子施加随机方向的微小位移，以便得到更加随机的结构，最后利用最小绝对收缩和选择算子的方法得到高阶力常数。
结合上述已经得到的介电张量、波恩有效电荷、二阶力常数和三阶力常数作为输入，最后通过ShengBTE声子玻尔兹曼输运方程求得简谐近似下晶格热导率、声子群速度、声子散射率等热输运系数。
结合介电张量、波恩有效电荷、二阶力常数和四阶力常数作为输入进行自洽声子理论的计算，这将得到与温度相关的非谐动力学矩阵，对其进行傅里叶变换得到与温度相关的有效二阶力常数，用有效二阶力常数代替上述的二阶力常数，在此基础上求解声子玻尔兹曼输运方程求得自洽声子理论下晶格热导率、声子群速度、声子散射率等热输运系数。
在简谐近似下存在虚频，与实验观察到的稳定性不一致。因此，传统的基于简谐原子间力常数的玻尔兹曼输运方程解的结果不再有效，并且非谐效应在足够大的位移振幅内变得显著。在 ALAMODE 软件包中，通过有限位移法在 2×2×2 超胞中获得简谐原子间力常数 (IFCs) ，位移大小设置为 0.01Å。
随后，通过压缩感知晶格动力学方法训练三阶和四阶力常数。具体来说，本文首先模拟了在 300 K、时间步长为 2fs 的 4000 步的从头算分子动力学，得到了 80 个快照。为了得到准随机结构，本文将到得到快照中的原子向随机方向移动 0.1 Å。随后，本文将获得的准随机结构在 4×4×4 的 k 点网格中进行了静态 DFT 计算，以计算 Hellmann–Feynman 力。最终，利用得到的 80 个准随机结构的位移和力组成的数据集通过最小绝对收缩和选择算子 (least absolute shrinkage and selection operator，简称 LASSO)方法训练出非谐力常数。
在计算中本文考虑了超胞内的所有原子间最近邻相互作用来计算简谐和三阶力常数，四阶力常数只考虑了第五最近邻内的相互作用，由于五阶和六阶力常数对结果的影响很小，因此只考虑了次近邻相互作用。
本文根据 ALAMODE 软件包中采用的 SCP 理论，利用获得的简谐和非谐力常数来计算非谐声子频率。
声子玻尔兹曼输运方程的解是通过 ShengBTE 的修订版 FourPhonon 软件包来实现的，并建立了具有高斯涂抹 (scale broad 参数为 0.05) 的 18×18×18 的 q 点网格来模拟相关积分和声子波矢量。此外，在计算中得到了大约 2.3×105 个三声子过程和 1.2×109 个四声子过程。因此，对于三声子散射，本文采用迭代方案来求解玻尔兹曼输运方程，而通过单模弛豫时间近似(SMRTA) 来处理四声子散射过程，这是因为计算成本巨大。

步骤分析

上面第3步，利用alm来算大概分为三步，mode设为suggest产生随机位移，vasp求力和能量，mode设为opt拟合二阶力常数即可。alm模块。
上面第4步，利用AIMD产生随机位移，vasp求力和能量，压缩感知晶格动力学方法训练三阶和四阶力常数mode设为cv，随后mode设为opt拟合高阶力常数，alm模块。最后mode设为SCPH获得变温声子谱和有效二阶力常数,anphon模块。后续运行dfc2模块生成指定温度的有效二阶力常数，用https://github.com/ttadano/alamode/discussions/192 中提到的脚本转换成phonopy格式。
上面第5步，为了产生ShengBTE的力常数输入文件，在拟合高阶力常数时，alm 的输入文件&general 部分加入 FC3_SHENGBTE=1、 FC4_SHENGBTE=1即可。如果简谐声子谱有虚频，上面的有效二阶力常数也需要，即第7步所述。
上面第6步见ShengBTE教程3.2节。
Alamode可以得到的结果包括格林艾森常数、声子群速度、三声子相空间、热导率原子贡献率和元素、变温声子谱和态密度、MSD、可视化振动模式、QHA等,anphon模块。

Alamode 连接 shengBTE

首先，github上安装最新版本的develop branches 或 feature/shengbte_4ph branches

[3rd 4th IFCs]

拟合力常数时，alm 的输入文件&general 部分加入 FC3_SHENGBTE=1、 FC4_SHENGBTE=1即可。

&general

PREFIX = BN-1000

MODE = opt

NAT = 128

NKD = 2; KD = B N

TOLERANCE = 1.0e-3

FC3_SHENGBTE=1

FC4_SHENGBTE=1

Note：输出的文件名为：FORCE_CONSTANT_3RD、FORCE_CONSTANT_4TH，和ShengBTE的文件名差一个S，FORCE_CONSTANTS_3RD、FORCE_CONSTANTS_4TH

[2nd IFCs]

首先通过QE的DFPT计算一个不含SCP的二阶力常数，BN.ifc2。

做SCPH计算，会生成计算温度点的二阶力常数SCP等效修正部分，例如BN.scph_dfc2（仅考虑四声子散射对重整化的影响），若进一步考虑三声子散射的重整化作用，文件名为BN.scph+bubble(0)_dfc2。

alamode/tools/scph_to_qefc.py 可以将SCP的修正部分与QE计算的二阶力常数相加，得到指定温度重整化后的二阶力常数。

Usage:

python scph_to_qefc.py original_QEfc2 scph_fc2_correction temperature > scp_fc2

python scph_to_qefc.py BN.ifc2 BN.scph_dfc2 1000 > espresso.ifc2

Authors and references

Weizhe Yuan yuanwz@stu.hit.edn.cn

W. Yuan et al., Physical Review B 110, 165412 (2024).

注意的问题

The accuracy of the anharmonic force constants can be checked by increasing the size of the training dataset. If the physical properties, such as thermal conductivity and free energy, change significantly with the size of the training dataset, it indicates that the training is insufficient. If the change is minor, your force constants may be OK.
报错：Warning: The following force constant doesn’t exist in the original file:
You can neglect this warning as long as the correction term (the last column) is small.
This warning is raised when the absolute value of the correction is larger than 1.0e-10 and the corresponding force constant does not exist in the original FCSXML file. In your case, the correction term is only slightly above the tolerance. So, it should not cause any problems hopefully.
When the correction term is much larger than 1.0e-10, this warning indicates a possible error in the inputs. Likely errors are the input structure is not fully symmetrized, or the numerical digits of lattice vectors and fractional coordinates are insufficient.
the given KMESH_INTERPOLATE is not commensurate with the supercell size of the original FCSXML.
In the current implementation of alamode, only the former choice,
renormalized 2FCs + conventional (old) 3FCs,
is supported. Another possible choice may be
renormalized 2FCs + renormalized 3FCs,
but it is not supported within alamode yet.
When MODE = RTA, the anphon code only reads the 2FCs and 3FCs from given XML files. By contrast, when MODE = SCPH, the code additionally reads the 4FCs because they are necessary to compute the renormalized 2FCs based on the self-consistent phonon theory.
Therefore, the second combination of
renormalized 2FCs + high order FCs using SCPH
is not possible.
I performed the SCPH calculation. It is not converged at low temperatures (0K, 100K, 200K), while it achieved convergence at high temperatures (300~1200K).A potential solution is to use a smaller MIXALPHA and a smaller DT.
I would like to use second-order, third-order, and fourth-order force constants to obtain anharmonic potentials and calculate the potential energy as a function of atomic displacement. My goal is to generate a plot similar to the one in the paper “Orbitally driven giant phonon anharmonicity in SnSe” (Nature Phys. 2015).Please use https://github.com/ttadano/alamode/blob/develop/tools/calcpes.py.
now that I have obtained the renormalized force constant file xxx.xml at a finite temperature, I now want to convert the renormalized force constant file directly to FORCE_CONSTANS
How to convert a second-order force constant matrix obtained by almode to FORCE_CONSTANTS in phonopy format？

python convert.py alamode.xml > FORCE_CONSTANTS

Could you give me some suggestion? How shoud I covert this with phonopy?

FULL_FORCE_CONSTANTS = .TRUE.
WRITE_FORCE_CONSTANTS = .TRUE.
Note: the old FORCE_CONSTANTS will be overwrited

How small the fitting error should be?

It depends on the Taylor expansion potential and displacement magnitude you choose.
In the standard harmonic calculation where –mag=0.01 is used in displace.py and the all harmonic interactions are considered (cutoff = None), the fitting error is usually less than 5% (~1–2% in most cases).
In the calculation of third-order force constants with –mag=0.04, the fitting error should be small as well. Indeed, in many cases, we obtain much smaller fitting errors (< 1%) than the harmonic case.
In the temperature-dependent effective potential method, we try to fit the harmonic potential to the displacement-force datasets sampled by ab initio molecular dynamics at finite temperature. Therefore, the fitting error tends to be much larger (> 10%).

dynaphopy

简谐声子谱的计算。
dynaphopy的非谐修正是在phonopy的力常数基础上进行的，需要先进行0k声子谱的计算。
结构优化的INCAR，力的收敛标准是E-3，大部分的文献也是在这个精度范围内。结构优化的结构是单胞，10个原子。
这是phonopy计算的声子谱，计算的结构为222的超胞，80个原子，在gmma点和X点有很强的虚频。这个虚频即使提高收敛精度和扩大超胞也是无法消除的。
AIMD计算
dynaphopy通过读取AIMD的计算结构进行拟合（说法不准确），需要预先准备AIMD的outcar或xdatcar。
AIMD的计算设置，基本参考dynaphopy的example中的设置，详细的参数涵义，都有教程介绍。
这里只是采用了3000步的计算，步长是2fs，能量的收敛精度是E-6，不算太高。同样是采用80个原子的超胞。大部分做非谐计算的文献，所采用的都是在80-160原子数范围内的超胞，步长多是1-2fs，总步数大概是3000-6000，总的时长是5ps-10ps。从下面的非谐声子谱结果来看，这个设置对于非谐计算是合适的？
dynaphopy非谐声子谱计算。
dynaphopy的官网如下：Dynaphopy (abelcarreras.github.io)。
将结构优化的原始晶胞POSCAR，phonopy输出的力常数，和AIMD的OUTCAR放在一个文件夹内。
设置dynaphopy的输入文件input如下：分别是输入结构文件，力常数文件，原胞和超胞与POSCAR的关系设置，band是声子谱的高对称点路径，设置方式和vasp计算PBE能带一致。
然后输入dynaphopy input OUTCAR -i，进入下的界面，选择6，进入非谐声子谱计算，再选择1，可以同时输出简谐声子谱和非简谐声子谱，此时会对peak进行拟合。拟合结束后，会输出声子谱。
一开始使用默认设置，输出的声子谱，可以看到高温声子谱虚频加重了。从@get-it 大佬处了解到，dynaphopy默认使用最大熵方法拟合，效果可能不如fftw有效。于是小弟将拟合方式改为fftw,通过指令dynaphopy input OUTCAR -i -psm 3,采用快速傅里叶变换拟合。得到300k下的声子谱，可以看到，原来gamma点和X点的虚频，已经被成功消除了。而实验上，这个材料在室温下肯定是稳定的。

当然，如果想要得到更准确的声子谱，需要提高优化精度，采用更大的超胞计算力常数，AIMD需要更多的步数和更小的步长，比如dynaphopy的例子是采用0.7 fs,跑200000步。如果只是想得到一个比较合理的声子谱，适当降低精度是没问题的。比如这个计算，在24核的机器上，结构优化用时不到3小时，dfpt计算1小时，AIMD计算14小时，差不多用时一天，能够得到一个比较合适的声子谱。如果提高精度，计算用时可能会提升数倍，但是声子谱的提升可能并不会太明显。

3.3 电声耦合

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

MC-电声耦合

需要vasp6.0之后的版本,INCAR开启

PHON_LMC = .TRUE.
IBRION = 6

全MC采样

标签PHON_NSTRUCT设置由于MC抽样而生成的结构的数量。应该监视观测值相对于这个数量的收敛。
标签TEBEG=0也是需要的，用于选择运行抽样的温度。
另外，PHON_LBOSE可以设置为.TRUE.或.FALSE. (默认PHON_LBOSE=.TRUE.)，分别选择玻色-爱因斯坦统计或麦克斯韦-玻尔兹曼统计。

0K的一个INCAR文件如下:

System = DEFAULT
PREC = Accurate
ISMEAR = 0; SIGMA = 0.1;
IBRION = 6

PHON_LMC = .TRUE.
PHON_NSTRUCT = 100
TEBEG = 0.0

随后生成Wycoff positions畸变的许多POSCAR，输出为POSCAR.TEBEG.NUMBER

ZG configuration即时采样（one-shot）

这种方法仅使用单一的扭曲结构，因此比完全的MC抽样快上几个数量级，同时保持了与收敛的超晶胞尺寸的完全MC抽样非常接近的准确性。例如，我们展示了对于带隙的零点重整化，准确度在ZG配置和完全MC抽样之间在5毫电子伏特之内。因此，我们建议优先使用这种方法，当很难实现超晶胞尺寸的收敛或者5毫电子伏特的精度已经足够时。

输入
要选择ZG配置，必须在INCAR文件中设置PHON_NSTRUCT=0。不同温度的数量和温度列表（以K为单位）必须使用标签PHON_NTLIST和PHON_TLIST，在INCAR文件中分别提供。一个示例如下：
1
2
PHON_NTLIST = 4
PHON_TLIST = 0.0 100.0 200.0 350.0

给出了一个温度范围为0-700K（步长为100K）的示例INCAR文件：

System = DEFAULT
PREC = Accurate
ISMEAR = 0; SIGMA = 0.1;
IBRION = 6
PHON_NTLIST = 8
PHON_TLIST = 0.0 100.0 200.0 300.0 400.0 500.0 600.0 700.0
PHON_NSTRUCT = 0
PHON_LMC = .TRUE.

输出

类似于MC抽样，ZG配置方法会产生多个具有不同扭曲Wycoff位置但不变的布里渊矩阵的POSCAR文件。这些文件被标记为

POSCAR.TEMP
其中TEMP遍历由PHON_TLIST定义的所有温度。

one-shot例子（金刚石带隙重整化）

输入文件

INCAR：

general:
 System = cd-C
 PREC = Accurate
 ALGO = FAST
 ISMEAR = 0 
 SIGMA = 0.1;
 IBRION = 6         #获得动力学矩阵声子计算
 PHON_LMC = .TRUE.  #开启电声耦合计算
 PHON_NSTRUCT = 0   #选择one-shot方法
 PHON_NTLIST = 1    #温度点个数
 PHON_TLIST = 0.0   #对应温度

准备优化后的结构，需要扩胞，以保证声子准确。KPOINTS、POTCAR根据POSCAR调整。

计算

一、获得结构

按照超胞结构和上面的INCAR提交任务，结束后获得OUTCAR。
执行： cp OUTCAR OUTCAR.init
新的畸变后的结构输出为POSCAR.T=0.

二、计算特定位移畸变结构的电子能级

cp POSCAR.T=0. POSCAR

INCAR 修改为

 System = cd-C
PREC = Accurate
ALGO = FAST
ISMEAR = 0 
SIGMA = 0.1  #去掉PHON_标签

执行计算，完成后cp OUTCAR OUTCAR.T=0.

三、提取ZPR

具有特殊位移的带隙减去没有特殊位移的带隙，就可以得到零点带隙重整的值。带隙保存在OUTCAR文件里。

执行例子中的脚本：bash extract_zpr.sh,
脚本内容：

#!/bin/bash 

i="OUTCAR.T=0."
j="OUTCAR.init"

homo1=`awk '/NELECT/ {print $3/2}' $i`
homo2=`awk '/NELECT/ {print $3/2-1}' $i`
homo3=`awk '/NELECT/ {print $3/2-2}' $i`
lumo1=`awk '/NELECT/ {print $3/2+var+1}' $i`
lumo2=`awk '/NELECT/ {print $3/2+var+2}' $i`
lumo3=`awk '/NELECT/ {print $3/2+var+3}' $i`
lumo4=`awk '/NELECT/ {print $3/2+var+4}' $i`
lumo5=`awk '/NELECT/ {print $3/2+var+5}' $i`
lumo6=`awk '/NELECT/ {print $3/2+var+6}' $i`
e1a=`grep "   $homo1  " $i | head -$nkpt | sort -n -k 2 | tail -1 | awk '{print $2}'`
e1b=`grep "   $homo2  " $i | head -$nkpt | sort -n -k 2 | tail -1 | awk '{print $2}'`
e1c=`grep "   $homo3  " $i | head -$nkpt | sort -n -k 2 | tail -1 | awk '{print $2}'`
e2a=`grep "   $lumo1  " $i | head -$nkpt | sort -n -k 2 | head -1 | awk '{print $2}'`
e2b=`grep "   $lumo2  " $i | head -$nkpt | sort -n -k 2 | head -1 | awk '{print $2}'`
e2c=`grep "   $lumo3  " $i | head -$nkpt | sort -n -k 2 | head -1 | awk '{print $2}'`
e2d=`grep "   $lumo4  " $i | head -$nkpt | sort -n -k 2 | head -1 | awk '{print $2}'`
e2e=`grep "   $lumo5  " $i | head -$nkpt | sort -n -k 2 | head -1 | awk '{print $2}'`
e2f=`grep "   $lumo6  " $i | head -$nkpt | sort -n -k 2 | head -1 | awk '{print $2}'`   
 
homo_ref=`awk '/NELECT/ {print $3/2}' $j`
lumo_ref=`awk '/NELECT/ {print $3/2+var+1}' $j`

h_ref=`grep "   $homo_ref  " $j | head -$nkpt | sort -n -k 2 | tail -1 | awk '{print $2}'`
l_ref=`grep "   $lumo_ref  " $j | head -$nkpt | sort -n -k 2 | head -1 | awk '{print $2}'`  
 
echo "The band gap (in eV) without zero-point vibrations is:"
echo "$h_ref $l_ref" |awk '{print ($2-$1)}'
echo "The band gap (in eV) including zero-point vibrations is:"
echo "$e1a $e1b $e1c $e2a $e2b $e2c $e2d $e2e $e2f" |awk '{print (($4+$5+$6+$7+$8+$9)/6.0-($1+$2+$3)/3.0)}'
echo "The zero-point renormalization of the band gap (in eV) is:" 
echo "$e1a $e1b $e1c $e2a $e2b $e2c $e2d $e2e $e2f $h_ref $l_ref" |awk '{print (($4+$5+$6+$7+$8+$9)/6.0-($1+$2+$3)/3.0)-($11-$10)}'

输出如下：

The band gap (in eV) without zero-point vibrations is:
4.4049
The band gap (in eV) including zero-point vibrations is:
4.05102
The zero-point renormalization of the band gap (in eV) is:
-0.353883

精度

精度取决于超胞大小。超胞越大越准确。

温度依赖的带隙计算

输入文件

INCAR：

 System = cd-C
 PREC = Accurate
 ALGO = FAST
 ISMEAR = 0 
 SIGMA = 0.1;
 IBRION = 6         #获得动力学矩阵声子计算
 PHON_LMC = .TRUE.  #开启电声耦合计算
 PHON_NSTRUCT = 0   #选择one-shot方法
PHON_NTLIST = 8     #温度个数
PHON_TLIST = 0.0 100.0 200.0 300.0 400.0 500.0 600.0 700.0

准备优化后的结构，需要扩胞，以保证声子准确。KPOINTS、POTCAR根据POSCAR调整。

计算

一、获取特定位移的结构

按照超胞结构和上面的INCAR提交任务，结束后获得OUTCAR。
执行： cp OUTCAR OUTCAR.init
新的畸变后的结构输出为POSCAR.T=0. ~ POSCAR.T=700.

二、计算特定位移畸变结构的电子能级
cp POSCAR.T=0. POSCAR

INCAR 修改为

 System = cd-C
PREC = Accurate
ALGO = FAST
ISMEAR = 0 
SIGMA = 0.1  #去掉PHON_标签

对每个温度下的结构执行计算。可通过脚本run_temperature.sh批量提交。
脚本内容：

#!/bin/bash

# please enter your executable path here
vasp_exec=./vasp_gam

#please enter the number of processors used for VASP here
np=8


for i in 0 100 200 300 400 500 600 700
do
   cp POSCAR.T\=$i. POSCAR
   mpirun -np $np $vasp_exec
   mv OUTCAR OUTCAR.T\=$i
done

三、分析OUTCAR

利用脚本extract_temp.sh

!/bin/bash 

if [ -f gap_vs_temp.dat ]
then
   rm gap_vs_temp.dat
fi

touch gap_vs_temp.dat
counter=0

for temp in 0 100 200 300 400 500 600 700
do
   i="OUTCAR.T=$temp"

   homo1=`awk '/NELECT/ {print $3/2}' $i`
   homo2=`awk '/NELECT/ {print $3/2-1}' $i`
   homo3=`awk '/NELECT/ {print $3/2-2}' $i`
   lumo1=`awk '/NELECT/ {print $3/2+var+1}' $i`
   lumo2=`awk '/NELECT/ {print $3/2+var+2}' $i`
   lumo3=`awk '/NELECT/ {print $3/2+var+3}' $i`
   lumo4=`awk '/NELECT/ {print $3/2+var+4}' $i`
   lumo5=`awk '/NELECT/ {print $3/2+var+5}' $i`
   lumo6=`awk '/NELECT/ {print $3/2+var+6}' $i`
   e1a=`grep "^    $homo1   " $i | head -$nkpt | sort -n -k 2 | tail -1 | awk '{print $2}'`
   e1b=`grep "^    $homo2   " $i | head -$nkpt | sort -n -k 2 | tail -1 | awk '{print $2}'`
   e1c=`grep "^    $homo3   " $i | head -$nkpt | sort -n -k 2 | tail -1 | awk '{print $2}'`
   e2a=`grep "^    $lumo1   " $i | head -$nkpt | sort -n -k 2 | head -1 | awk '{print $2}'`
   e2b=`grep "^    $lumo2   " $i | head -$nkpt | sort -n -k 2 | head -1 | awk '{print $2}'`
   e2c=`grep "^    $lumo3   " $i | head -$nkpt | sort -n -k 2 | head -1 | awk '{print $2}'`
   e2d=`grep "^    $lumo4   " $i | head -$nkpt | sort -n -k 2 | head -1 | awk '{print $2}'`
   e2e=`grep "^    $lumo5   " $i | head -$nkpt | sort -n -k 2 | head -1 | awk '{print $2}'`
   e2f=`grep "^    $lumo6   " $i | head -$nkpt | sort -n -k 2 | head -1 | awk '{print $2}'`

   if [ $temp -eq "0" ]
   then
      ref=`echo "$e1a $e1b $e1c $e2a $e2b $e2c $e2d $e2e $e2f" |awk '{print (($4+$5+$6+$7+$8+$9)/6.0-($1+$2+$3)/3.0)}'`
   fi

   echo "$e1a $e1b $e1c $e2a $e2b $e2c $e2d $e2e $e2f $temp $ref" |awk '{print $10,(($4+$5+$6+$7+$8+$9)/6.0-($1+$2+$3)/3.0)-$11}' >> gap_vs_temp.dat
done

运行脚本后生成gap_vs_temp.dat文件，画图即可。

gnuplot画图：
gnuplot -e "set terminal jpeg;set xlabel 'T (K)'; set ylabel 'band gap'; set style data lines; plot 'gap_vs_temp.dat', 'C_exp_points_offset0.dat' w circles, 'C_exp_fit_offset0.dat'" > gap_vs_temp.jpg

包含体积效应的温度依赖带隙

仅考虑温度效应会低估带隙，要考虑体积效应。

一、准简谐近似方法获得体积变化

三个脚本：

1
2
3

quasi_harm_4x4x4_diamond_create_pos_and_run_vasp.sh
quasi_harm_4x4x4_diamond_make_energy_vs_volume_plots.sh
quasi_harm_4x4x4_diamond_obtain_fitting.sh

对超胞结构创建不同体积的结构进行体积-能量计算：

运行脚本：quasi_harm_4x4x4_diamond_create_pos_and_run_vasp.sh

#! /bin/bash

# please enter your executable path here
vasp_exec=./vasp_gam

#please enter the number of processors used for VASP here
np=8

cp INCAR.qh INCAR

for i in 6.13521592 6.27789536 6.4205748 6.56325424 6.70593368 6.84861312 6.99129256 7.133972 7.27665144 7.41933088 7.56201032 7.70468976 7.8473692 7.99004864 8.13272808
do
sed "s/7.13397200/${i}/g" POSCAR.4x4x4 > POSCAR_$i
cp POSCAR_$i POSCAR
mpirun -np 8 $vasp_exec

mv OUTCAR OUTCAR_$i
done

该脚本创建15个POSCAR，体积每步变化2%,运行脚本quasi_harm_4x4x4_diamond_create_pos_and_run_vasp.sh提交计算。

二、计算完成后运行quasi_harm_4x4x4_diamond_make_energy_vs_volume_plots.sh数据处理。
脚本内容：

#!/bin/bash

bandshift=0
gwrun=-1
dgbd=-1
val=-1
con=-1
gwldadiff=-1
test=-1
while [[ $# -gt 0 ]]
do
   key="$1"
   case $key in
   esac
   shift
done

if [ -f "helpscript.perl" ]; then
   rm helpscript.perl
fi


cat > helpscript.perl  <<EOF
#!/bin/perl

use strict;
use warnings;

my \$zahler=0;
my @entropy=0;
my \$ezp=0;
my \$fhelmholtz=0;
my \$uenergy;
my \$kboltzmann=8.6173303*10**(-5.0);
my \$ntemp=8;
my \$tmax=700;
my \$tmin=0;
my \$tstep=100;
for (my \$itemp=1;\$itemp<=\$ntemp;\$itemp++)
{
   \$entropy[\$itemp]=0;
}
while (<>)
{
   chomp;
   \$_=~s/^/ /;
   my @help=split(/[\t,\s]+/);
   \$zahler=\$zahler+1;
   if (\$zahler == 1) {\$uenergy=\$help[1];}
   else 
   {
      my \$homega=\$help[2]/1000;
      \$ezp=\$ezp+\$homega*0.5;
      for (my \$itemp=1;\$itemp<=\$ntemp;\$itemp++)
      {
         my \$temp=(\$itemp-1)*\$tstep+\$tmin;
         my \$kbt=\$kboltzmann*\$temp;
         if (\$temp < 0.0000001) 
         {
         }
         else
         {
              \$entropy[\$itemp]=\$entropy[\$itemp]-\$kboltzmann*log(1-exp(-\$homega/\$kbt));
              \$entropy[\$itemp]=\$entropy[\$itemp]+\$kboltzmann*\$homega/\$kbt*(1/(exp(\$homega/\$kbt)-1));
         }
      }
   }
   last if eof;
}

\$ezp=\$ezp; #/ \$zahler;

for (my \$itemp=1;\$itemp<=\$ntemp;\$itemp++)
{
    my \$temp=(\$itemp-1)*\$tstep+\$tmin;
        
    \$entropy[\$itemp]=\$entropy[\$itemp]; # / \$zahler; 
    \$fhelmholtz=\$uenergy+\$ezp-\$temp*\$entropy[\$itemp];
    printf ("%15.8e %15.8e\n",\$temp,\$fhelmholtz);
}
printf ("#NOTEMP %15.8e\n",\$uenergy);
EOF

rm OUTTEMP*
first=0
for i in OUTCAR_*; do

   echo "Starting $i"
   v2=${i##OUTCAR_}
   if [ -f "helpfile.help" ]; then
      rm helpfile.help
   fi
   touch helpfile.help
   cp OUTCAR_$v2 OUTCAR
   awk '/free energy/' OUTCAR | tail -n 1| awk '{print $5}' >> helpfile.help
   awk '/[0-9]* f .* THz/ {print $1,$10}' OUTCAR >> helpfile.help
   awk '/[0-9]* f.i.*THz/ {print $1,$9}' OUTCAR >> helpfile.help
   volume=`awk '/volume of cell/ {print $5}' OUTCAR | tail -n 1`

   perl helpscript.perl helpfile.help > hhhhelp.txt

   runcount=0
   while read line; do
      runcount=$((runcount+1))
      if [[ $first -eq 0 ]]; then
         echo $line | awk -v var="$volume" '{print var,$2,$1}' > OUTTEMP_$runcount
      else
         echo $line | awk -v var="$volume" '{print var,$2}' >> OUTTEMP_$runcount
      fi
   done < ./hhhhelp.txt

   first=1
done

for i in OUTTEMP_*; do
   v2=${i##OUTTEMP_}
   mv OUTTEMP_$v2 OUTHELP
   sort -n -k 1 OUTHELP > OUTTEMP_$v2
done

rm helpfile.help
rm helpscript.perl
rm hhhhelp.txt
rm OUTHELP

自由能曲线保存在OUTTEMP_*.

运行脚本quasi_harm_4x4x4_diamond_obtain_fitting.sh获得温度-体积关系。
脚本内容：

#!/bin/bash

for i in OUTTEMP_*
do
   cp $i OUTTEMP.current
   #extract temperature
   temp=`head -n 1 OUTTEMP.current|awk '{print $3}'`
   #do fitting
   gnuplot -e "E(V)=E0+9.0/8.0*B0*V0*((V0/V)**(2.0/3.0)-1)**2 + 9.0/16.0*B0*\
      (B0P-4)*V0*((V0/V)**(2.0/3.0)-1.0)**3.0 + R*((V0/V)**(2.0/3.0)-1.0)**4.0;\
      B0P = 1;B0 = 1;V0 = 720;E0 = -1150;R  = -1.0;fit E(x) 'OUTTEMP.current'  u 1:2 via B0P,B0,V0,E0,R" &> suppress_output
   #extract volume from fit
   a=`grep "V0" fit.log|grep "=" |tail -n 1|awk '{print ($3/2.0)**(1.0/3.0)}'`
   #print temperature and volume to 
   echo "temperature: $temp, a_latt: $a" 
done

rm suppress_output
rm OUTTEMP.current

输出如下：

temperature: 0.00000000e+00, a_latt: 7.18012
temperature: 1.00000000e+02, a_latt: 7.18014
temperature: 2.00000000e+02, a_latt: 7.18064
temperature: 3.00000000e+02, a_latt: 7.18267
temperature: 4.00000000e+02, a_latt: 7.18613
temperature: 5.00000000e+02, a_latt: 7.19037
temperature: 6.00000000e+02, a_latt: 7.19493
temperature: 7.00000000e+02, a_latt: 7.19959
temperature: #NOTEMP, a_latt: 7.15218

三、最后将QHA获得的晶胞参数，替换掉POSCAR.T=* files中的晶胞参数。再次运行run_temperature.sh 和extract_temp.sh.获得数据画图。

通过添加体积效应，与实验获得了更好的一致性。在这个教程中，我们只使用了一个4x4x4的单元格，因为更大的单元格已经相当耗时，但对于收敛的5x5x5单元格，两条曲线都应该与实验相比略有变差。由于在本例中使用的PBE无法充分描述电子交换和相关性，因此预计实验与理论之间会有差异。为了获得真正的优异一致性，需要使用GW近似方法。
严格来说，将体积效应添加到电子-声子相互作用的正确方法是首先为每个温度更改体积，然后为该温度计算电子-声子相互作用。在这个教程和参考文献[5]中，我们是反过来做的。因此，电子-声子相互作用只需要计算一次。在参考文献[5]中，我们观察到这两种方法给出了非常相似的结果。

3.2 热输运2

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

三阶力常数矩阵

需要用到thirdorder程序。

使用方法：

准备好POSCAR,INCAR,KPOINTS,POTCAR和提交任务脚本（以sub.sh为例）；
其中POSCAR与二阶力常数矩阵求解时的POSCAR-unit一致，KPOINTS也与二阶一致。INCAR文件与二阶力常数求解时的有限位移法一致，其他不变。完成后运行thirdorder实现扩胞。

~/安装目录/thirdorder_vasp.py sow 2 2 2 -12
2 2 2是扩胞倍数，要和二阶力常数矩阵保持一致或者略小。-12指的是考虑多少个近邻原子，需要测试，但是如果计算资源紧缺的时候取个-10就行了，但有风险；正值是指截断半径。之后生成成功以后会出现一堆3RD开头的文件，那些都是POSCAR，一般有几百上千个。

参考phonopy计算声子谱的有限位移法，需要把这些POSCAR单独计算。也就是有多少POSCAR，建立多少文件夹，分别进行多少次计算。
find job* -name vasprun.xml|sort -n|/ thirdorder_vasp.py reap a b c -d a b c和-d与前面一致。该命令可获得FORCE_CONSTANT_3RD也即是三阶力常数矩阵。

脚本方法

测试-d参数（秦光照老师）

{.line-numbers}

#!/bin/bash
#
# Construct the 3nd files based on the finished calculations
#
# Written BY QIN, GuangZhao <qin.phys@gmail.com>
# 2016-11-22

DIR="/vol-th/home/xhaiyan/WangN/Au6S2/third241-14"

if [ ! -d "$DIR" ];then
    echo Please specify the source directory.
    exit
fi

for NAME in 3RD.POSCAR.*
do
    NUM="${NAME#3RD.POSCAR.}"   # the number of the 3nd POSCAR file
    CODE=$(head -1 $NAME)       # the code of the POSCAR file
# find the source POSCAR file in the source directory
    INFO=$(grep $CODE $DIR/3RD.POSCAR.* | head -1)
    [ x"$INFO" == x ] && echo "NOT FOUND for job-$NUM" && continue
    SOURCE_NUM="${INFO%:*}"     # the number of the source POSCAR file
    SOURCE_NUM="${SOURCE_NUM##*.}"

    echo "$DIR/job-$SOURCE_NUM/ --> job-$NUM/"
#    diff 3RD.POSCAR.$NUM $DIR/3RD.POSCAR.$SOURCE_NUM
    rsync -a $DIR/job-$SOURCE_NUM/vasprun.xml ./job-$NUM/
done

针对大量POSCAR批量建立文件夹

{.line-numbers}

#module load spglib/master
#thirdorder_vasp.py sow 3 3 1 -10
#生成2*2*2的超胞 0.6表示截断半径为6埃。
#sow a b c m  abc表示生成的超胞大小，m为正数的时候表示截断半径（单位：nm）（例如0.6，表示截断半径为6.0埃）；为负数时表示第几近邻截断（此时应取负整数，例如-3表示第三近邻）

for i in 3RD.POSCAR.*
do
s=$(echo $i| cut -d "." -f 3)
d=job-$s
mkdir $d
cp $i $d/POSCAR
cp ./INCAR ./sub.sh ./KPOINTS ./POTCAR $d
pushd $d
popd
done

批量提交任务

{.line-numbers}

#!/bin/sh
root_path='pwd'

#在每个job文件夹里进行静态计算。
for i in {001..1050} #按顺序从job-001到1050做静态vasp计算,根据生成的文件夹设定
do
cd job-$i
echo "job-$i", time 
sbatch sub.sh
wait
cd ..
done

#上述静态计算完成后提取三阶力常数文件
#cd ${root_path}

#find job* -name vasprun.xml| sort -n| python thirdorder_vasp.py reap 2  2  1 -12 #最后四个数字跟你扩胞时一致

脚本赋予权限：chmod +x .sh;运行脚本：./sh

ShengBTE获得热输运性质

运行ShengBTE计算热输运性质

获得二阶力常数计算：
计算流程与声子谱计算流程相同，用DFPT和有限位移法都可以。将得到的FORCE_CONSTANTS重命名为FORCE_CONSTANTS_2ND。
获得三阶力常数计算：见上，输出为FORCE_CONSTANTS_3RD文件。

准备CONTROL文件：格式如下:

{.line-numbers}

#CONTROL for ShengBTE
&allocations
        nelements=2              #元素种类
        natoms=6                 #原子个数
        ngrid(:)=20 20 5         #计算网格密度
&end
&crystal
        lfactor=0.1              #软件长度单位换算
        lattvec(:,1)=3.1613337552503711    0.0000000000000025    0.0000000000000000 #原胞晶胞参数
        lattvec(:,2)=-1.5806668776251651    2.7377953418477841    0.0000000000000000
        lattvec(:,3)=-0.0000000000000000    0.0000000000000000   12.3318234981499426
        elements= "Mo" "S"  #原子顺序
        types= 1 1 2 2 2 2  #原子个数，2代表第二种元素，4个2代表有4个第二类原子
       positions(:,1)= 0.3333333449999998  0.6666666890000030  0.2500000000000000 #原子坐标，按顺序
       positions(:,2)=0.6666666240000012  0.3333333229999980  0.7500000189999980
       positions(:,3)= 0.3333333449999998  0.6666666890000030  0.6230719397502427
       positions(:,4)=0.6666666240000012  0.3333333229999980  0.3769280602499874
       positions(:,5)=0.6666666240000012  0.3333333229999980  0.1230720067501656
       positions(:,6)=0.3333333449999998  0.6666666890000030  0.8769280222497613
       scell(:)=4 4 1 #扩胞倍数，和前面一致
&end
&parameters
       T=300  #温度
&end
&flags
      nonanalytic=.TRUE.
      nanowires=.FALSE.
      convergence=.TRUE.  #迭代法，默认最大1000步
&end

详细参数可参考手册。

提交任务
获得不同温度的数据，只需要依次修改CONTROL文件的T tag。或者通过脚本：

{.line-numbers}

for i in {350..1350..100} # 对i赋值，区间为300到900，每搁10取一个数，即（300 310 320 330 ....890 900），可以根据自身情况自行设置
do  #开始循环
mkdir $i #建立单一任务文件夹
cp sub.sh CONTROL FORCE_CONSTANTS_2ND FORCE_CONSTANTS_3RD $i  
#将二阶三阶力学常数文件和CONTROL文件复制到单一任务文件夹
#注意，这里的CONTROL文件夹是删减过的
cd $i #进入单一文件夹内
echo \&parameters >> CONTROL      #补充CONTROL文件
  echo  T= $i  >> CONTROL    #补充CONTROL文件 
  #温度设置 如果同时设置多个温度，同一行只会算第一个并报错。如果多次设定一个温度到同一个CONTROL文件，只会运算最后一个。
#echo   scalebroad=0.1    >> CONTROL    #补充CONTROL文件
echo \&end    >> CONTROL    #补充CONTROL文件
echo \&flags  >> CONTROL    #补充CONTROL文件
echo    nonanalytic=.TRUE.   >> CONTROL    #补充CONTROL文件
echo    nanowires=.FALSE.   >> CONTROL    #补充CONTROL文件
echo    convergence=.True. >> CONTROL
echo \&end  >> CONTROL    #补充CONTROL文件
 #这里补充CONTROL文件用的是笨方法，如果有相应改进方法可以自行更改，简洁代码
sbatch sub.sh  #执行ShengBTE计算，这里的执行命令应视情况而定。
cd .. #返回初始目录
done  #结束任务

数据处理

计算完成后，先检查收敛情况，找到T300K文件夹下的BTE.kappa,第一列代表迭代步数，一般未达到默认1000步可认为收敛正常。
计算完成后在当前目录下会生成T300K文件夹和一些输出文件。文件处理参见天玑算b站视频：https://www.bilibili.com/video/BV1Yh411n7tx/?spm_id_from=333.788&vd_source=17084cc867ff3bb86398ff1a79b443f2

其他参考：https://www.bilibili.com/read/cv22090733/

注意事项

计算得到的热导率和文献报道有较大偏差，可能偏大或者偏小，可能的原因有：

1.采用的2阶力常数FORCE_CONSTANT中超胞内原子排序不是按照ucatomsc_zsc_y*sc_x排列；力常数文件与control文件中原胞内原子次序是否一致；

2.shengBTE的control文件中晶格参数单位默认nm，有转换系数的参数；

3.使用QE利用phonopy计算得到的FORCE_CONSTANT文件中，单位需要转换，不能直接利用shengBTE计算，vasp利用phonopy则单位正确；
1. 3阶力常数的计算量较大，但精度要求很高，如0.0003Rh/bohr^2量级的力的误差，也会导致最终热导率的数值出现大的偏离；如使用cp2k计算力常数，通常使用QE等平面波软件，对于超胞体系，随着扩胞，在较大的截断半径下，所需计算量非常大，使用CP2K中的QS模块也能进行第一性原理计算，具有相对较高的计算效率，但是CP2K不同于平面波的QE等，计算中仅有gamma点计算，需要采用实空间较大超胞，如果使用和QE等相同的超胞，则会导致原子受力精度不高，一般在~0.0003Rh/bohr^2量级的误差，也会对某些体系，特别是二维体系的热导率有很大的影响。在使用cp2k进行3阶力常数大量计算之前，先进行超胞收敛测试，可以在相同的位移构型下，不断扩胞，判断力计算的变化情况和误差。用QE等计算也建议先进行k点等计算参数的收敛判断，否则会浪费大量时间，热导率计算结果却有很大偏离。
2. 使用经验势能模型，在通过lammps进行热输运计算之前，可以先利用分子静力学软件GULP结合shengBTE模块，预判势能应用于热输运计算的准确性。但经验势能用于热输运精度和前景一般，机器势能表现更好。
原文链接：https://blog.csdn.net/odin_linux/article/details/130084321
静态计算可以各个job文件夹同时分别计算，不用按顺序来。等所有job文件夹计算完后，自行运行“find job…” 提取三阶力常数的命令。严格来说，胞的大小和三阶力常数计算的截断半径都应该进行测试，测试以所计算得到的热导率的收敛进行判断。计算资源有限无法测试的话，尽可能取大一点。Born有效电荷和介电常数的计算应该不是必要的，如果你要考虑在Gamma点附近LO-TO劈裂的话(主要针对离子晶体)，则需要进行计算。可以加上，感觉要好点。
二阶力常数和三阶力常数所用的超胞大小可以不统一，我看文献中有这种情形，不过一般作者采用的二阶力常数的超胞都是大于等于三阶力常数超胞的大小。
原文链接：https://zhuanlan.zhihu.com/p/147810223
官方网站：https://www.shengbte.org/

四声子过程

除了ShengBTE的CONTROL文件中的常规输入，FourPhonon还需要四阶力常数和一些新的标签在CONTROL文件中。请检查ShengBTE网站以获取其他标签的定义。

并行环境

1.1和1.0版本是为MPI并行性编写的。从1.2版本开始，支持四声子散射的迭代求解器，我们已经迁移到OpenMP，以处理该迭代求解器所需的大内存。未来我们可能会支持MPI+OpenMP混合并行性，以允许更多灵活性。确保在编译时添加-qopenmp，或者：

export FFLAGS=-qopenmp -traceback -debug -O2 -static_intel

四阶IFCs文件：FORCE_CONSTANTS_4TH

该文件包含四阶原子间力常数矩阵，并使用稀疏描述来节省空间。要构建四阶力常数，可以参考Fourthorder的python脚本。该力常数的格式是对三阶力常数的直接扩展，包含nb个这样的矩阵块：

一行空白
一个从1开始的顺序索引（从1到nb）
一行表示第二个晶胞的笛卡尔坐标（单位：Å）
一行表示第三个晶胞的笛卡尔坐标（单位：Å）
一行表示第四个晶胞的笛卡尔坐标（单位：Å）
一行表示四个相关原子的从1开始的索引，每个原子索引从1到natoms
81行的力常数矩阵（单位：eV Å⁴）。开头的索引用于标记笛卡尔坐标轴。
该文件的一个示例：


1
0.0000000000e+00 0.0000000000e+00 0.0000000000e+00
0.0000000000e+00 0.0000000000e+00 0.0000000000e+00
0.0000000000e+00 0.0000000000e+00 0.0000000000e+00
     1      1      1      1
 1  1  1  1       -42.0002584218
 1  1  1  2         0.0000000000
 1  1  1  3         0.0000000000
 1  1  2  1         0.0000000000
 1  1  2  2        40.3875689003
 1  1  2  3         0.0000000000
 1  1  3  1         0.0000000000
 1  1  3  2         0.0000000000
 1  1  3  3        40.3875689003
 1  2  1  1         0.0000000000
 1  2  1  2        40.3875689003
 1  2  1  3         0.0000000000
 1  2  2  1        40.3875689003
 1  2  2  2         0.0000000000
 1  2  2  3         0.0000000000
 1  2  3  1         0.0000000000
 1  2  3  2         0.0000000000
 1  2  3  3         0.0000000000
 1  3  1  1         0.0000000000
 1  3  1  2         0.0000000000
 1  3  1  3        40.3875689003
 1  3  2  1         0.0000000000
 1  3  2  2         0.0000000000
 1  3  2  3         0.0000000000
 1  3  3  1        40.3875689003
 1  3  3  2         0.0000000000
 1  3  3  3         0.0000000000
 2  1  1  1         0.0000000000
 2  1  1  2        40.3875689003
 2  1  1  3         0.0000000000
 2  1  2  1        40.3875689003
 2  1  2  2         0.0000000000
 2  1  2  3         0.0000000000
 2  1  3  1         0.0000000000
 2  1  3  2         0.0000000000
 2  1  3  3         0.0000000000
 2  2  1  1        40.3875689003
 2  2  1  2         0.0000000000
 2  2  1  3         0.0000000000
 2  2  2  1         0.0000000000
 2  2  2  2       -42.0002584218
 2  2  2  3         0.0000000000
 2  2  3  1         0.0000000000
 2  2  3  2         0.0000000000
 2  2  3  3        40.3875689003
 2  3  1  1         0.0000000000
 2  3  1  2         0.0000000000
 2  3  1  3         0.0000000000
 2  3  2  1         0.0000000000
 2  3  2  2         0.0000000000
 2  3  2  3        40.3875689003
 2  3  3  1         0.0000000000
 2  3  3  2        40.3875689003
 2  3  3  3         0.0000000000
 3  1  1  1         0.0000000000
 3  1  1  2         0.0000000000
 3  1  1  3        40.3875689003
 3  1  2  1         0.0000000000
 3  1  2  2         0.0000000000
 3  1  2  3         0.0000000000
 3  1  3  1        40.3875689003
 3  1  3  2         0.0000000000
 3  1  3  3         0.0000000000
 3  2  1  1         0.0000000000
 3  2  1  2         0.0000000000
 3  2  1  3         0.0000000000
 3  2  2  1         0.0000000000
 3  2  2  2         0.0000000000
 3  2  2  3        40.3875689003
 3  2  3  1         0.0000000000
 3  2  3  2        40.3875689003
 3  2  3  3         0.0000000000
 3  3  1  1        40.3875689003
 3  3  1  2         0.0000000000
 3  3  1  3         0.0000000000
 3  3  2  1         0.0000000000
 3  3  2  2        40.3875689003
 3  3  2  3         0.0000000000
 3  3  3  1         0.0000000000
 3  3  3  2         0.0000000000
 3  3  3  3       -42.0002584218

CONTROL 文件设置

此 CONTROL 文件包含所有用户指定的设置和参数，包括晶体结构信息、展宽因子、q 网格、温度、函数等。

一般用法

要启用 FourPhonon 功能，您需要添加一个新的 &flags 名称列表：

1
2

four_phonon （逻辑，默认为 .false.）：计算四声子相空间和四声子散射率。
four_phonon_iteration （逻辑，默认为 .false.）：精确计算四声子散射并求解玻尔兹曼输运方程（BTE）。这仅在 four_phonon=.true. 时有效。

标志的实际用法
以下是这些标志的一些实用组合：

onlyharmonic=.true. 和 four_phonon=.true.：仅计算四声子相空间。
convergence=.false. 和 four_phonon=.true.：在 RTA 层面计算三声子和四声子散射的热导率；此设置将输出不同四声子通道的贡献。
four_phonon=.true. （收敛默认值为 .true.）：使用三声子迭代方案计算热导率，但将四声子散射处理为 RTA 层面。
four_phonon=.true. 和 four_phonon_iteration=.true.：利用 BTE 的精确解计算三声子和四声子散射的热导率。

注意： CONTROL 文件中的所有其他参数，例如温度或 q 网格，适用于三声子和四声子过程。FourPhonon 包不支持纳米线功能。

加速 RTA 解的采样方法

为了从散射过程样本中估计散射率，请引用 [Z. Guo 等人，npj Comput. Mater. 10, 31 (2024).] 并使用相关标签。

参数设置

num_sample_process_3ph 和 num_sample_process_3ph_phase_space

num_sample_process_3ph 和 num_sample_process_3ph_phase_space （整数，默认值为 -1）：表示从每个模式中提取的样本数，以用于估计三声子相空间和散射率。 -1 表示不进行取样，执行严格计算。注意，当 convergence=.true.，即使用迭代方案计算三声子散射时，num_sample_process_3ph 必须为 -1，因为取样方法基于弛豫时间近似。

num_sample_process_4ph 和 num_sample_process_4ph_phase_space

num_sample_process_4ph 和 num_sample_process_4ph_phase_space （整数，默认值为 -1）：表示从每个模式中提取的样本数，以用于估计四声子相空间和散射率。 -1 表示不进行取样，执行严格计算。注意，当 four_phonon=.false. 时，num_sample_process_4ph 和 num_sample_process_4ph_phase_space 必须均为 -1。

示例

计算三声子散射率（精确计算，迭代），并在 RTA 层面采用 100,000 个样本过程计算四声子散射：

&parameters
    T=300
    scalebroad=1
    num_sample_process_3ph_phase_space = -1
    num_sample_process_3ph = -1
    num_sample_process_4ph_phase_space = 100000
    num_sample_process_4ph = 100000         
&end

TDEP IFCs 的输入

对于 TDEP 格式的谐波力常数，请使用以下描述的标签/脚本：

tdep=.true. 启用 TDEP 格式谐波 IFCs 的处理，该输入文件名应始终为 infile.forceconstant。如果系统是极性材料，除了常规参数（如 Born 有效电荷和介电常数）外，我们还需要用户手动输入耦合参数 (\lambda_{Ewald})，该参数在 TDEP IFCs 的末尾近似写出，如下所示：

1	0.675614929199219 # Coupling parameter in Ewald summation

注意：输入的谐波 IFCs 是未带极性修正的裸 IFCs，以便 FourPhonon 能够处理并在内部添加修正。该模块在很大程度上依赖 TDEP 源代码，我们无法保证成功。建议用户检查处理后的频率和三声子散射率。

对于非极性材料，我们还开发了一种脚本，将 TDEP 格式转换为 Phonopy 格式，名为 tdep2phonopy.f90。这样您无需激活 tdep 标志。编译 Fortran 后，运行可执行文件（以下是示例）：

1	./tdep2phonopy outfile.forceconstant 10 10 1 2

这将 TDEP IFCs （outfile.forceconstant）转换为 Phonopy 格式，构建一个来自双原子单位胞的 (10 \times 10 \times 1) 超胞。请检查转换后的声子色散。我们建议使用足够大的超胞来存储所有的第二阶 IFC 中的成对相互作用。

示例

计算三声子散射率（精确计算，迭代），并在 RTA 层面计算四声子散射。使用的 TDEP 格式谐波 IFCs 并应用非解析修正，Ewald 参数为 (\lambda =0.675614929199219)。同样需要 Born 有效电荷（此处未列出），与 ShengBTE 对极性系统的处理相似。

&parameters
    T=300
    scalebroad=0.1
    Ewald=0.675614929199219
&end
&flags
    tdep=.TRUE.
    nonanalytic=.TRUE. ! 此标签现在需要明确指定
    convergence=.TRUE.
    four_phonon=.TRUE.
&end

TDEP 格式的 3阶/4阶 IFCs

对于 TDEP 格式的 3阶/4阶 IFCs，我们开发了一个独立的脚本，将其转换为 ShengBTE 格式，名为 tdep2ShengBTE.f90。我们需要在同一目录下有 infile.ucposcar。编译 Fortran 后，运行可执行文件：

示例

1	./tdep2ShengBTE 3 2 outfile.forceconstant_thirdorder

这将 TDEP 的 3 阶 IFCs（outfile.forceconstant_thirdorder）转换为 ShengBTE 格式。

Output files

Besides the routine output files from previous ShengBTE program, FourPhonon generates these output files:

BTE.Numprocess_4ph: number of allowed four-phonon scattering processes, for each irreducible q point and phonon band
BTE.P4: phase space available for four-phonon processes, for each irreducible q point and phonon band
BTE.P4_total: total volume in phase space available for four-phonon processes
BTE.P4_plusplus*, BTE.P4_plusminus*, BTE.P4_minusminus*: similar to BTE.P4 but only includes contributions from ++/+-/- - processes

*Note for four-phonon scatterings, there are three different channels: recombination (++), redistribution (+-) and splitting (- -) processes.

Under temperature-dependent directories: (the unit of scattering rates is ps−1 if not specified)

BTE.WP4: weighted phase space available for four-phonon processes ([rad/ps]$^{-5}$, 2nd column) vs angular frequency (rad/ps, 1st column) for each irreducible q point and phonon band
BTE.WP4_plusplus*, BTE.WP4_plusminus*, BTE.WP4_minusminus*: similar to BTE.WP4 but only includes contributions from ++/+-/- - processes
BTE.w_3ph: three-phonon scattering rates for each irreducible q point and phonon band, this file replaces the original output BTE.w_anharmonic. Absorption and emission processes are written out into BTE.w_3ph_plus and BTE.w_3ph_minus
BTE.w_4ph: four-phonon scattering rates for each irreducible q point and phonon band. Similarly, we provide the contributions from different channels: BTE.w_4ph_plusplus, BTE.w_4ph_plusminus and BTE.w_4ph_minusminus
BTE.w_4ph_normal: four-phonon scattering rates of normal processes, for each irreducible q point and phonon band. This file has 5 columns and each one represents: angular frequency in rad/ps, recombination channel, redistribution channel, splitting channel, overall scattering rates from normal processes
BTE.w_4ph_Umklapp: four-phonon scattering rates of Umklapp processes, for each irreducible q point and phonon band. The format is the same as BTE.w_4ph_normal
BTE.kappa*: thermal conductivity results are written out as usual

3.1 热输运1 声子谱

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

声子谱

声子谱的计算主要有两种方法，一种是直接法，另一种是微扰密度泛函方法（DFPT）。

直接法，或称frozen-phonon方法，是通过在优化后的平衡结构中引入原子位移，计算作用在原子上的Hellmann-Feynman力，进而由动力学矩阵算出声子色散曲线。直接法的缺陷在于它要求声子波矢与原胞边界supersize
正交，或者原胞足够大使得Hellmann-Feynman力在原胞外可以忽略不计。这使得对于复杂系统，如对称性高的晶体、合金、超晶格等材料需要采用超原胞。超原胞的采用使计算量急剧增加，极大的限制了该方法的使用。
1987年，Baroni、Giannozzi和Testa提出了DFPT方法。DFPT通过计算系统能量对外场微扰的响应来求出晶格动力学性质。该方法最大的优势在于它不限定微扰的波矢与原胞边界正交，不需要超原胞也可以对任意波矢求解。Castep、Vasp等采用的是一种linear response theory 的方法（或者称为 density perturbation functional theory，DFPT），直接计算出原子的移动而导致的势场变化，再进一步构造出动力学矩阵。进而计算出声子谱。

DFPT计算步骤

结构优化
声子谱的计算需要对原胞结构做高精度的充分优化，否则很容易出现虚频；
{.line-numbers}
1
2
EDIFF = 1E-08
EDIFFG = -1E-06
可以采用分步优化，直到达到精度。然后优化后的结果，cp CONTCAR POSCAR：

$ phonopy -d --dim="4 4 1"

之后会生成一堆文件，其中有 SPOSCAR 和POSCAR-00* ，分别对应不同的方法密度泛函微扰法 (DFPT) 以及有限位移法。

{.line-numbers}

$ ls
phonopy_disp.yaml  POSCAR-002  POSCAR-005  POSCAR-008  POSCAR-011  POSCAR-014  SPOSCAR
POSCAR             POSCAR-003  POSCAR-006  POSCAR-009  POSCAR-012  POSCAR-015
POSCAR-001         POSCAR-004  POSCAR-007  POSCAR-010  POSCAR-013  POSCAR-016
$ cp POSCAR POSCAR-unitcell
$ cp SPOSCAR POSCAR

KPOINTS需适当减小，可以的话最好再进行一次收敛测试
INCAR

{.line-numbers}

IBRION =  8            (determines the Hessian matrix using DFPT)
EDIFF  =  1E-08        (SCF energy convergence; in eV) 
PREC   =  Accurate     (Precision level)  
ENCUT  =  500          (Cut-off energy for plane wave basis set, in eV) 
IALGO  =  38           (Davidson block iteration scheme)
LREAL  = .FALSE.       (Projection operators: false)
LWAVE  = .FLASE.       (Write WAVECAR or not)
LCHARG = .FLASE.       (Write CHGCAR or not) 
ADDGRID= .TRUE.        (Increase grid; helps GGA convergence) 
NSW    = 1

上面参数 IBRION = 8 NSW = 1是注意修改的。
提交VASP计算

后处理
phonopy --fc vasprun.xml

读取vasprun.xml，生成FORCE_CONSTANTS文件。

编辑band.conf文件，该文件给出了高对称点路径的信息
Note:

{.line-numbers}

ATOM_NAME = * *  #元素名称，顺序与POSCAR保持一致
DIM = 4 4 1      #建立超胞的大小，--dim后面的
BAND = 0.5 0.0 0.0  0.0 0.0 0.0  0.333333 0.333333 0.0  0.5 0.0 0.0 #高对称点路径，每组高对称点之间用两个空格隔开
FORCE_CONSTANTS = READ #表示读取力常数文件FORCE_CONSTANTS
FC_SYMMETRY = .TRUE.

准备mesh.conf文件，如下所示：

{.line-numbers}

1
2
3

ATOM_NAME = Si
DIM = 2 2 2
MP = 24 24 24

{.line-numbers}

phonopy --dim="4 4 1" -c POSCAR-unitcell -p -s band.conf
(单位是THz)或者
phonopy --dim="4 4 1" --factor=524.147 -c POSCAR-unitcell -p -s band.conf
(单位是cm-1)

（读取band.conf、FORCE_CONSTANTS、POSCAR-unitcell，输出band.yaml、phonopy.yaml以及pdf图）

1
2
3

phonopy-bandplot  --gnuplot> phono.dat
如果上述命令不行，可以用旧版本phonopy 
bandplot --gnuplot> phonon.dat

读取band.yaml，就可以获得二阶力常数矩阵声子谱的数据，其中第一行为高对称点的坐标。

准备mesh.conf以计算态密度和PDOS：

{.line-numbers}

ATOM_NAME = 
DIM = 4 4 1
MP=12 12 1
PDOS=1 2, 3 4, 5 6
FORCE_CONSTANTS = READ
BAND_POINTS = 71

运行命令：phonopy –dim='4 4 1' -c POSCAR-unitcell mesh.conf或者乘上因子修改单位。phonopy –dim='4 4 1' --factor=524.147 -c POSCAR-unitcell mesh.conf
得到projected_dos.dat文件，拖到Origin中画图，得到声子态密度图像（从第二列数据开始，每一列对应每个原子的态密度），对比着之前计算得到的声子谱可以分析声学支与光学支之间的相互作用关系

有限位移方法

这里需要准备和 POSCAR-00* 数目一样的文件夹，并把 POSCAR-00* 复制到每个文件夹中并命名为 POSCAR 。每个文件夹内都需要准备 POTCAR 、 POSCAR 、 INCAR 、 KPOINTS 。并且分别计算。
可用以下sh脚本快速进行：手动修改第一行文件夹数目和vasp运行文件名称。

{.line-numbers}

for i in {001..006}
do 
mkdir $i
cd $i
cp ../POSCAR-$i POSCAR
cp ../{INCAR,KPOINTS,POTCAR,sub.sh} .
运行提交任务命令
cd ../
done

INCAR 可以设置为如下内容：
INCAR设置如下(静态计算)：

{.line-numbers}

PREC = Accurate
IBRION = -1
NSW = 1
ENCUT = 500
EDIFF = 1.0e-08
EDIFFG = -0.001
ISMEAR = 0
SIGMA = 0.05
IALGO = 38
LREAL = .FALSE.
LWAVE = .FALSE.
LCHARG = .FALSE.

里面IBRION = -1 NSW = 1注意修改

KPOINTS需适当减小，可以的话最好再进行一次收敛测试
提交VASP计算

后处理

{.line-numbers}

1	phonopy -f ./disp-*/vasprun.xml

准备band.conf

{.line-numbers}

ATOM_NAME = In Se
DIM = 4 4 1
BAND = 0.5 0.0 0.0  0.0 0.0 0.0  0.333333 0.333333 0.0  0.5 0.0 0.0
FORCE_CONSTANTS = WRITE  #注意这里是 WRITE
FC_SYMMETRY = .TRUE.

{.line-numbers}

phonopy -c POSCAR-unitcell band.conf -p -s
phonopy --factor=524.147 -c POSCAR-unitcell -p -s band.conf
#旧版本phonopy
bandplot  --gnuplot> phonon.dat
就可以获得二阶力常数矩阵 FORCE_CONSTANTS 和声子谱的数据
#新版本phonopy
phonopy-bandplot --gnuplot > phonon.dat

Mind

如果计算完的声子谱其他地方都好，就是在Γ点有一点点虚频，那么这个材料很可能是稳定的，只是你优化做得不够好，进一步提高优化的精度可消除这一点点虚频。
对于二维材料，如果在Γ点出现很小的虚频，基本可以认为这个材料是稳定的，大部分二维材料都会有此现象；尤其是VASP结合phonopy code计算二维材料的声子谱在Γ点更是容易出现虚频；使用Quantum-ESPRESSO的PWSCF和PH模块计算声子谱对于内存的需求较小一些，且对于二维材料的声子谱计算更友好一些，后期会详细介绍Quantum-ESPRESSO计算声子谱的具体步骤。
声子谱是表示组成材料原子的集体振动模式。如果材料的原胞包含n个原子，那么声子谱总共有3n支，其中有3条声学支，3n-3条光学支。声学支表示原胞的整体振动，光学支表示原胞内原子间的相对振动。

计算出的声子谱有虚频，往往表示该材料不稳定。 有些情况下，我们可以利用虚频信息使不稳定的材料变得稳定。如果声子谱的一条声学支存在虚频，主要位于Γ点和M点1/2处（对应倒格矢的1/4位置）。倒格矢的1/4，对应晶格长度的4倍。我们可能需要将原胞沿上述倒格矢方向扩大四倍，进一步优化原子位置，才可能得到比较稳定的晶胞。

消虚频

DFPT方法与有限位移法互换尝试。或者更换标准胞扩胞计算。
扩胞，在虚频很小时可以尝试扩大扩胞倍数。
提高精度优化后计算声子谱。声子谱对力的精度要求特别高，所以提高精度多次优化后再计算声子谱也是可行的。
修改POTCAR，尝试不同的赝势方法。（不太建议，因为需改之后所有计算都需要统一成稳定构型的赝势。不过有效）
修改POSCAR原子坐标，再优化，再计算声子谱。
修改方法：对于DFPT计算的声子谱，使用grep cm-1 OUTCAR查看频率，寻找最大的虚频频率。然后在OUTCAR中从最后往前找到最大的虚频频率对应的坐标。
dx dy dz对应虚频的位移，把这些坐标与POSCAR中的坐标相加或者0.1-0.9之间的倍数，形成新的POSCAR。进行优化，再计算声子谱。
不排除真的算不出来没虚频的，据我所知，SnTe就没见到谁大大方方发文章谈声子谱的（算出来虚频消除不了不好意思讲）。
INCAR参数对虚频的影响：ADDGRID影响很小，LREAL有影响。截断能和K点密度在收敛情况下影响小。

结构优化消虚频

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

phonopy声子群速度

操作步骤

首先对结构高精度优化后扩胞，计算得到没有虚频的声子谱；
对band.conf文件进行修改，并添加群速度计算相关参数GROUP_VELOCITY = .TRUE.和GV_DELTA_Q = 0.01，程序需要读取本征值和本征矢，所以必须设置 EIGENVECTORS =.TRUE.，此外为了防止vaspkit在导出数据时报错，需要删除BAND_LABELS高对称路径名称，修改后band.conf文件内容如下；

而且具体使用过程中，关于执行命令还有执行文件的tag内容还有一些可替换的关系,–gv, –group_velocity 可在执行计算命令时替代执行文件中的GROUP_VELOCITY = .TRUE.内容， –gv_delta_q 可在执行计算命令时替代执行文件中的(GV_DELTA_Q)内容

对声子谱重新进行计算，打开phonopy.yaml文件可以看到，phonopy在计算声子谱的过程中考虑到了群速度；打开band.yaml文件中也保存了有关群速度的相关信息；同时在计算声子谱中phonopy所储存的数据文件mesh.yaml和band.yaml文件中也保存了有关群速度的相关信息。
随后进行数据提取，以vaspkit1.12版本为例，打开vaspkit，输入73，打开VASP2other Interface 子功能，选择739 Sort Phonon Band Structure for Phononpy，运行完成即可获得group velocity相关数据以及各原子的声子谱数据；在运行过程中我们可以看到软件说明find group velocity并写入输出文件中。
1.5版本为789功能。

格林艾森常数（Grüneisen parameters）

为表示声子非简谐效应的参数。phonopy计算模式格林艾森参数方式见官网。
以NaCl为例，首先需要进行三个声子谱计算，分别是平衡体积（orig），稍小体积（minus）和稍大体积（plus）的结构，计算时改变晶格参数使得体积获得改变。

计算完成后获得不同体积下的原子间力常数文件FORCE_SETS 或 FORCE_CONSTANTS （后续命令中需要添加–readfc参数）。

随后通过phonopy-gruneisen（旧版gruneisen）命令获得对应于声子带结构中倒易空间中的路径的模式格林艾森常数（具体设定可根据自身计算情况更改）
phonopy-gruneisen orig plus minus --dim="2 2 2" --pa="0 1/2 1/2 1/2 0 1/2 1/2 1/2 0" --band="1/2 1/4 3/4 0 0 0 1/2 1/2 1/2 1/2 0.0 1/2" -p -c POSCAR-unitcell
在倒易网格中的模式格林艾森常数（具体设定可根据自身计算情况更改）：
phonopy-gruneisen orig plus minus --dim="2 2 2" --pa="0 1/2 1/2 1/2 0 1/2 1/2 1/2 0" --mesh="20 20 20" -p -c POSCAR-unitcell --color="RB"
phonopy-gruneisen 命令选项包括–dim，–mp, –mesh，–band，–pa, –primitive_axis，–readfc，–band_points，–nac，–factor，–nomeshsym，-p，-c，-s, –save，-o等，功能上均与phonopy命令相同，但除上述列出的选项外，可以在phonopy命令上使用的其它选项是不能在phonopy-gruneisen 上使用的。如有必要，可以自行修改phonopy-gruneisen代码。

如查看帮助，可输入
phonopy-gruneisen -h
通过 phonopy-gruneisenplot 命令导出计算得到的数据，使用方式类似于phonopy-bandplot。

原文链接：https://blog.csdn.net/icehoqion/article/details/130550807

热学性质

参考：
http://supersunsir.com/2020/08/01/vasp/vasp_phonopy/

热膨胀系数

我们通过准简谐近似（QHA）的方法计算热膨胀系数和等压热容等性质，

首先，我们要把优化好的POSCAR加应变（比如0.96 0.97 0.98 0.99 1 1.01 1.02 1.03 1.04的缩放系数，这样就有9个结构文件），
然后把每一个POSCAR都去重复步骤（2）也就是去计算每一个缩放系数声子谱。大概就是这个效果
最后我们通过phonopy里面的QHA模块去得到此材料的非偕性质phonopy-qha -p -s e-v.dat *00/thermal_properties.yaml就可以得到很多热力学性质比如热膨胀系数，格林爱森参数，自由能，等压热容等等。

先进行高精度的晶格优化。
修改缩放系数为0.95-1.05新建10个POSCAR.
对11个POSCAR晶体进行声子谱计算
thermo.conf文件后处理生成11个对应的thermal_properties.yaml
QHA计算热膨胀系数

另一个写的步骤：

高精度的原胞优化
在此基础上施加形变（梯度可以设置为千分之一）
优化多个结构
使用DFPT方法计算每个结构的声子谱（如果结构有问题，使用phonopy –symmetry）
在每个文件夹中得到thermal_properties.yaml
运行脚本得到e-v.dat
运行phonopy-qha的相关命令得到各种dat文件
Over

可以参考作者的thermo.conf文件：

{.line-numbers}

DIM = 2 2 2
ATOM_NAME = sio
MP = 31 31 31
TPROP = .TRUE.
TMAX = 2000
TSTEP = 2
FORCE_CONSTANTS = READ

然后是如何写e-v.dat文件

大概意思就是找到之前计算的那11个POSCAR的vasprun.xml结果文件，然后记事本打开搜索最后一个e_wo_entrp。
下载phonopy官方文件中Si-QHA的example，
解压vasprun_xmls.tar.lzma我们可以得到算例中用到的vasprun.xml文件。

不同温度下热力学性质计算（自由能、热容、振动熵）

在DFPT或者Finite displacement方法计算声子谱的基础上，得到mesh.yaml文件
新建一个文件夹，将FORCE_CONSTANTS、POSCAR、mesh.conf放入该文件夹中（POSCAR为原胞）

得到mesh.yaml文件

{.line-numbers}

1	phonopy mesh.conf -c POSCAR

mesh.conf的内容

{.line-numbers}

ATOM_NAME = Hf Si O
DIM = 2 2 2
MP = 13 13 13
FORCE_CONSTANTS = READ

得到thermal_properties.yaml文件

phonopy -c POSCAR mesh.conf -t
得到thermal_properties.dat文件
phonopy-propplot --gnuplot thermal_properties.yaml >thermal_properties.dat
将dat文件导入到Origin中画图

dat文件的5列数据分别为：

{.line-numbers}

temperature[K]
Free energy [kJ/mol]
Entropy [J/K/mol]
Heat capacity [J/K/mol]
Energy [kJ/mol]

下面是能用到的各种命令

{.line-numbers}

扩胞用
phonopy -d --dim="2 2 2" -c POSCAR_UNIT
查找晶格对称性
phonopy --symmetry
收集力常数
phonopy --fc vasprun.xml
得到band.yalm文件（好像不需要这一步）
phonopy --dim="2 2 2" -c POSCAR-UNIT -p band.conf
得到thermal_properties.yaml文件
phonopy -c POSCAR mesh.conf  -t

下面是产生e-v.dat文件的脚本，需要根据自己的目录名进行修改：

{.line-numbers}

for i in -0.005 -0.004 -0.003 -0.002 -0.001 0.001 0.002 0.003 0.004 0.005
do
 cd AC_Uni$i;
 pwd;
 volume=`grep volume vasprun.xml | tail -1| awk '{printf("%20.13f", $3)}'`;
 energy=`grep 'E0' OSZICAR | tail -1 | awk '{printf "%12.4f \t", $5}'`;
 echo ${volume} ${energy} >> ../e-v.dat;
 cd ../; #Exit the Strain folder
done

声子谱振动可视化

上面生成band.yaml文件，利用JMOL打开。https://jmol.sourceforge.net/
或者网站：https://henriquemiranda.github.io/phononwebsite/phonon.html
band.conf添加标签ANIME=0 0 0将000这个Q点的振动模式输出到anime.ascii,用v_sim软件打开即可。https://l_sim.gitlab.io/v_sim/index.en.html
运行命令 phonopy -f vasprun.xml-{001,002}得到FORSETS
运行phonopy –dim=”3 3 3” –ct=”0 0 0”得到gamma点的振动模式密度
DFPT方法得到的OUTCAR可以利用PyVibMS、JMOL
从band.yaml生成vesta读取：https://github.com/AdityaRoy-1996/Phonopy_VESTA

注意不要输出声子群速度，写出声子本征矢量。需要POSCAR.vesta格式。

从OUTCAR生成vesta读取：https://github.com/Stanford-MCTG/VASP-plot-modes

其他参考：https://www.bilibili.com/read/cv11848434/

声子谱反折叠

https://github.com/yuzie007/upho

热位移

热学性质是根据倒易空间采样网格上的声子频率计算得出的，必须设置Mesh tag进行使用，并且必须检查它们相对于取样网格的收敛性。通常这种计算对计算量要求不高，因此随着采样网格密度的增加，收敛性很容易达到。此外，计算时Phonopy对频率的默认单位为Thz。
在DFPT和冷冻声子的计算基础上，在thermal_properties.conf 中设置：

TPROP = .TRUE. #进行热学性质计算
DIM = 1 1 1 #扩胞情况
MP = 16 16 8 #Mesh网格密度
TMAX = 2000
TMIN = 0
TSTEP = 10 #温度上下限默认为100,0，step默认值为10

计算得到thermal_properties.yaml文件：

phonopy -c POSCAR thermal_properties.conf -t

得到thermal_properties.dat文件：

phonopy-propplot --gnuplot thermal_properties.yaml >thermal_properties.dat

五列数据分别为：temperature[K]、Free energy [kJ/mol]、Entropy [J/K/mol]、Heat capacity [J/K/mol]、Energy [kJ/mol]。

Phonopy官网还给出了两个参数：
PRETEND_REAL = .TRUE.
这里将虚频声子当成实声子进入计算，仅在有小虚频实作测试用；
CUTOFF_FREQUENCY = 0.1
默认值为0，将虚频声子频率视为负值，计算高于截止频率的声子。

热位移计算

计算投影到笛卡尔轴上的均方位移与温度的函数关系。声子频率以 THz 为单位（phononopy 的默认设置），用于获得均方位移。同样需要与Mesh tag共同使用，温度设置与前述相同。新建thermal_displacements.conf用于输入tag。

TDISP = .TRUE. #进行均方位移计算
DIM = 1 1 1 #扩胞情况
MP = 16 16 8 #Mesh网格密度
TMAX = 2000
TMIN = 0
TSTEP = 10 #温度上下限默认为100,0，step默认值为10

XYZ_PROJECTION = .TRUE.
打开该tag，将对特征向量按实空间进行投影，得到不同原子在a,b,c三个轴向上的位移投影，数据存储为frequency atom1_x atom1_y atom1_z atom2_x atom2_y atom2_z…形式。

PROJECTION_DIRECTION = 1 1 0
打开该tag，将原子位移在指定实空间轴向上进行投影。
phonopy -c POSCAR thermal_displacements.conf -p -s
执行后，输出结果将保存在thermal_displacements.yaml中。注意，计算的原子位移为均方位移。

此外的两个tag：

TDISPMAT = .TRUE.
打开该tag时，将以3*3矩阵形式储存原子在六个方向上的均方位移xx, yy, zz, yz, xz, xy，结果写入thermal_displacement_matrices.yaml中。

TDISPMAT_CIF = 1273.0
该tag用于指定计算在某温度（1273K）下，原子的均方位移矩阵，结果输出在tdispmat.cif中。

Phonopy其他功能

phonopy --symmetry -c POSCAR-unit

BPOSCAR和PPOSCAR生成，分别是标准胞和原胞，vasp格式。

phonopy --rd 10 --dim 2 2 2 --pa auto --amplitude 0.03 -c POSCAR-unitcell

–rd 10 代表生成10个随机位移结构，–dim 2 2 2扩胞倍数，–pa auto沿三个坐标轴投影方向，–amplitude 0.03位移幅度，-c指定对象。

输出”phonopy_disp.yaml” and supercells

phonopy-load phonopy_disp_orig.yaml --rd 1000 --temperature 300

generate 1000 supercells with displacements at 300K.

New in v2.25 Symfc (symfc/symfc) is based on fitting approach and any displacements set of atoms in supercell can be handled. For example, random displacements generated by RANDOM_DISPLACEMENTS can be used to compute force constants. To use symfc, its python module has to be installed via conda-forge or pip.

New in v2.3 ALM (ttadano/ALM) is based on fitting approach and any displacements set of atoms in supercell can be handled. For example, random displacements generated by RANDOM_DISPLACEMENTS can be used to compute force constants. To use ALM, its python module has to be installed via conda-forge or building it. The installation instruction is found here.

When ALM is used, please cite the paper: T. Tadano and S. Tsuneyuki, J. Phys. Soc. Jpn. 87, 041015 (2018).

FC_CALCULATOR = ALM

2.4 3D能带

发表于 2023-05-09 更新于 2025-12-18 分类于 DFT学习

vaspkit处理

开始之前可以做好正常的scf和band计算，确定需要第几条带。在二维能带里，我们只看沿着高对称路径的点上的能带。在三维能带图中，我们需要计算特定带的布里渊区的所有点。

新建文件夹用于3D计算，可以复制scf。vaspkit的231功能用于产生布里渊区撒点的KPOINTS。所生成的KP点包含1加权KP点和0加权KP点(注意，为了获得平滑的3D能带，用于生成KP点文件的分辨率值应该设置<0.01)。
提交计算，INCAR同scf。

后处理

vaspkit

vaspkit 232/233都可用于后处理，232可选择特定的带，233仅适用于非磁性半导体体系。完成之后会出现包含价带顶的BAND.HOMO.grd和导带底的BAND.LUMO.grd文件的能带。grep ^band PROCAR | head -n30
执行python how_to_visual.py，绘制3D能带图(python2.7环境)；这个文件在VASPKIT的安装包中，解压之后路径为：vaspkit/vaspkit.1.3.5/examples/3D_band；
特别的，如果能带不够美观，可通过能带插值进行改善。VASPKIT升级到1.3.0或更新版本，在~/.vaspkit文件中设置以下三个参数
how_to_visual.m对应于MATLAB代码， how_to_visual.py 对应于python代码，均可进行对BAND_HOMO.grd BAND_LUMO.grd KX.grd KY.grd 数据文件的读取和绘图。使用how_to_visual.m代码在MATLAB中进行绘制可实现角度拖动和放大缩放和自定义颜色等操作。

详细参考：https://mp.weixin.qq.com/s/CHnphDzUgeZA3hD0cENx5Q

https://mp.weixin.qq.com/s/4KWdZh6H3utVkF77DhuHFg

https://mp.weixin.qq.com/s/DxLl-mEZsdxq5gbft2cpzg

qvasp后处理

参考：https://mp.weixin.qq.com/s/A_Kb5IaSwY2U0UgpcsghiQ

qvasp -3dband处理提取VBM和CBM能带，
下载数据，放于origin中画图。步骤是全选数据，然后Plot->3D Surface-> Color Map Surface,点击ok
origin中精修（例如去除网格线之类的），并把价带顶和导带底merge到一起。并调整色彩美感。

脚本处理

绘图contour：就是类似等高线，如想看两个能带是否有nodal-line，可以将两个能带做差投影到一个平面。

from mpl_toolkits.mplot3d.axes3d import Axes3D
from matplotlib import cm
import matplotlib.pyplot as plt
import numpy as np
import sys
num=201
X = np.arange(0, num, 1)
Y = np.arange(0, num, 1)
X, Y = np.meshgrid(X, Y)

file=open("14u.dat",'r+')
Z=file.readlines()
Z=[Z[i].strip().split() for i in range(len(Z))]
Z=[float(Z[i][0]) for i in range(len(Z))]
Z=np.array(Z)
Z.resize(num,num)

file2=open("15u.dat",'r+')
Z2=file2.readlines()
Z2=[Z2[i].strip().split() for i in range(len(Z2))]
Z2=[float(Z2[i][0]) for i in range(len(Z2))]
Z2=np.array(Z2)
Z2.resize(num,num)
Z3 = Z- Z2print(np.min(np.abs(Z3)))

fig = plt.figure(figsize=(10,10))
ax = fig.add_subplot(1, 1, 1)
c=plt.contourf(Y,X,Z3,150,alpha=1,cmap=plt.cm.coolwarm)
plt.contour(Y,X,Z3,3,colors='black',linewidths=0.3)
plt.xlim(0,num)
plt.ylim(0,num)
plt.colorbar(c)
#plt.clabel(C, inline=True,fontsize=10)
plt.savefig("contour-3d",dpi=600)

脚本2

利用vasplit-232处理得到的数据绘制，输入对应vaspkit输出的能带index，生成2D和3D的指定能带，能带差值以及差值的绝对值，能带劈裂能的统计以及txt格式的数据格式。

#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
绘制VASPKIT 232功能输出的自旋劈裂图
修改版：调整3D图z轴标签位置，避免与标题重叠
"""

import numpy as np
import matplotlib.pyplot as plt
from matplotlib import cm
from mpl_toolkits.mplot3d import Axes3D
import os
import sys
import glob

# =============================================
# 1. 读取grd文件 - 修改版
# =============================================

def read_grd_file_simple(filename):
    """
    简单的.grd文件读取函数
    假设文件包含多个数值，每行可能有多个值
    """
    print(f"正在读取文件: {filename}")
    
    with open(filename, 'r') as f:
        lines = f.readlines()
    
    # 收集所有数值
    all_values = []
    for line in lines:
        line = line.strip()
        if line:
            # 分割数值，支持空格和制表符分隔
            parts = line.split()
            for part in parts:
                try:
                    all_values.append(float(part))
                except ValueError:
                    print(f"警告: 跳过非数值: {part}")
    
    # 转换为numpy数组
    data = np.array(all_values)
    
    # 尝试确定网格尺寸
    n_total = len(data)
    
    # 如果是KX或KY文件，我们可能需要从文件名推断
    if 'kx' in filename.lower() or 'ky' in filename.lower():
        # 对于KX/KY文件，尝试找到匹配的能带文件大小
        # 这里我们假设能带文件是正确格式的
        return data, 0, 0  # 暂时返回0尺寸，后面再处理
    else:
        # 对于能带文件，尝试推断尺寸
        nx = int(np.sqrt(n_total))
        ny = nx
        
        while nx * ny < n_total and nx <= ny:
            nx += 1
        
        while nx * ny > n_total and nx > 1:
            nx -= 1
        
        if nx * ny != n_total:
            ny = int(np.ceil(n_total / nx))
        
        print(f"简单读取: 数据点总数 = {n_total}, 推断网格 = {ny}×{nx}")
        
        # 重塑数据
        data_2d = data.reshape((ny, nx))
        
        return data_2d, nx, ny

# =============================================
# 2. 获取用户输入的能带编号
# =============================================

def get_band_index():
    """
    获取用户输入的能带编号
    """
    print("=" * 60)
    print("VASPKIT 232功能数据可视化")
    print("=" * 60)
    
    # 获取用户输入
    while True:
        try:
            band_input = input("请输入要分析的能带编号 (例如: 23): ").strip()
            band_index = int(band_input)
            if band_index > 0:
                return band_index
            else:
                print("错误: 能带编号必须是正整数!")
        except ValueError:
            print("错误: 请输入有效的数字!")
        except KeyboardInterrupt:
            print("\n程序已终止")
            sys.exit(0)

# =============================================
# 3. 查找文件
# =============================================

def find_files(band_index):
    """
    根据能带编号查找对应的文件
    """
    # 定义可能的文件名模式
    band_patterns = [
        f"BAND_B{band_index}_UP.grd",  # 默认格式
        f"BAND_B{band_index}_UP.GRD",  # 大写扩展名
        f"band_b{band_index}_up.grd",  # 全小写
        f"B{band_index}_UP.grd",       # 可能没有BAND前缀
        f"BAND_{band_index}_UP.grd",   # 可能没有B前缀
    ]
    
    dw_patterns = [
        f"BAND_B{band_index}_DW.grd",
        f"BAND_B{band_index}_DW.GRD",
        f"band_b{band_index}_dw.grd",
        f"B{band_index}_DW.grd",
        f"BAND_{band_index}_DW.grd",
    ]
    
    # 查找UP文件
    up_file = None
    for pattern in band_patterns:
        matches = glob.glob(pattern)
        if matches:
            up_file = matches[0]
            break
    
    # 查找DW文件
    dw_file = None
    for pattern in dw_patterns:
        matches = glob.glob(pattern)
        if matches:
            dw_file = matches[0]
            break
    
    # 查找KX和KY文件
    kx_patterns = ["KX.grd", "kx.grd", "KX.GRD", "KX.dat", "kx.dat"]
    ky_patterns = ["KY.grd", "ky.grd", "KY.GRD", "KY.dat", "ky.dat"]
    
    kx_file = None
    for pattern in kx_patterns:
        matches = glob.glob(pattern)
        if matches:
            kx_file = matches[0]
            break
    
    ky_file = None
    for pattern in ky_patterns:
        matches = glob.glob(pattern)
        if matches:
            ky_file = matches[0]
            break
    
    return kx_file, ky_file, up_file, dw_file

# =============================================
# 4. 主程序
# =============================================

def main():
    # 获取能带编号
    band_index = get_band_index()
    print(f"分析能带编号: B{band_index}")
    
    # 查找文件
    kx_file, ky_file, up_file, dw_file = find_files(band_index)
    
    # 检查文件是否存在
    files_found = []
    missing_files = []
    
    for fname, ftype in [(kx_file, "KX"), (ky_file, "KY"), (up_file, "UP"), (dw_file, "DW")]:
        if fname:
            files_found.append((fname, ftype))
            print(f"找到{ftype}文件: {fname}")
        else:
            missing_files.append(ftype)
    
    if missing_files:
        print(f"\n错误: 以下文件未找到: {', '.join(missing_files)}")
        print("请检查:")
        print(f"1. 当前目录下是否有KX.grd/KY.grd文件")
        print(f"2. 是否有BAND_B{band_index}_UP.grd和BAND_B{band_index}_DW.grd文件")
        print(f"3. 尝试不同的文件名格式 (如B{band_index}_UP.grd)")
        sys.exit(1)
    
    print(f"\n成功找到所有文件!")
    
    # 读取数据
    print("\n读取k点网格数据...")
    kx_data = read_grd_file_simple(kx_file)
    ky_data = read_grd_file_simple(ky_file)
    
    print("\n读取能带数据...")
    band_up_data = read_grd_file_simple(up_file)
    band_dw_data = read_grd_file_simple(dw_file)
    
    # 提取数据
    kx_2d, nx_kx, ny_kx = kx_data
    ky_2d, nx_ky, ny_ky = ky_data
    band_up_2d, nx_up, ny_up = band_up_data
    band_dw_2d, nx_dw, ny_dw = band_dw_data
    
    # 如果KX/KY是1D数组，尝试重塑
    if len(kx_2d.shape) == 1:
        # 推断网格尺寸
        n_total = len(kx_2d)
        nx = int(np.sqrt(n_total))
        ny = nx
        
        while nx * ny < n_total and nx <= ny:
            nx += 1
        
        while nx * ny > n_total and nx > 1:
            nx -= 1
        
        if nx * ny != n_total:
            ny = int(np.ceil(n_total / nx))
        
        print(f"重塑KX: {n_total}点 -> {ny}×{nx}网格")
        kx_2d = kx_2d.reshape((ny, nx))
        ky_2d = ky_2d.reshape((ny, nx))
    
    # 确保所有数组形状一致
    target_shape = kx_2d.shape
    print(f"\n目标网格形状: {target_shape}")
    
    # 调整能带数据形状
    if band_up_2d.shape != target_shape:
        print(f"调整UP能带形状: {band_up_2d.shape} -> {target_shape}")
        if band_up_2d.size == target_shape[0] * target_shape[1]:
            band_up_2d = band_up_2d.reshape(target_shape)
        else:
            print("警告: UP能带数据点数量不匹配!")
    
    if band_dw_2d.shape != target_shape:
        print(f"调整DW能带形状: {band_dw_2d.shape} -> {target_shape}")
        if band_dw_2d.size == target_shape[0] * target_shape[1]:
            band_dw_2d = band_dw_2d.reshape(target_shape)
        else:
            print("警告: DW能带数据点数量不匹配!")
    
    # 计算自旋劈裂
    spin_splitting = band_up_2d - band_dw_2d
    
    # 统计信息
    print("\n" + "=" * 60)
    print(f"能带 B{band_index} 自旋劈裂统计信息:")
    print(f"  最小值: {np.nanmin(spin_splitting):.6f} eV")
    print(f"  最大值: {np.nanmax(spin_splitting):.6f} eV")
    print(f"  平均值: {np.nanmean(spin_splitting):.6f} eV")
    print(f"  标准差: {np.nanstd(spin_splitting):.6f} eV")
    print(f"  绝对平均值: {np.nanmean(np.abs(spin_splitting)):.6f} eV")
    
    # =============================================
    # 5. 绘制2D映射图
    # =============================================
    
    print("\n" + "=" * 60)
    print("绘制2D映射图...")
    
    # 创建子图
    fig_2d, axes_2d = plt.subplots(2, 2, figsize=(14, 12))
    fig_2d.suptitle(f'2D Spin Splitting Maps (UP - DW) - Band B{band_index}', fontsize=16, fontweight='bold')
    
    # 定义颜色范围
    vmax = max(abs(np.nanmin(spin_splitting)), abs(np.nanmax(spin_splitting)))
    if vmax == 0:
        vmax = 0.01
    
    # 1. Spin-up 能带
    ax1 = axes_2d[0, 0]
    im1 = ax1.pcolormesh(kx_2d, ky_2d, band_up_2d, 
                        cmap='viridis', shading='auto')
    ax1.set_title(f'Spin-up Band Energy (B{band_index})', fontsize=12)
    ax1.set_xlabel('$k_x$ (1/Å)', fontsize=10)
    ax1.set_ylabel('$k_y$ (1/Å)', fontsize=10)
    ax1.set_aspect('equal')
    plt.colorbar(im1, ax=ax1, label='Energy (eV)')
    
    # 2. Spin-down 能带
    ax2 = axes_2d[0, 1]
    im2 = ax2.pcolormesh(kx_2d, ky_2d, band_dw_2d, 
                        cmap='viridis', shading='auto')
    ax2.set_title(f'Spin-down Band Energy (B{band_index})', fontsize=12)
    ax2.set_xlabel('$k_x$ (1/Å)', fontsize=10)
    ax2.set_ylabel('$k_y$ (1/Å)', fontsize=10)
    ax2.set_aspect('equal')
    plt.colorbar(im2, ax=ax2, label='Energy (eV)')
    
    # 3. 自旋劈裂 (UP - DW)
    ax3 = axes_2d[1, 0]
    im3 = ax3.pcolormesh(kx_2d, ky_2d, spin_splitting, 
                        cmap='RdBu_r', shading='auto',
                        vmin=-vmax, vmax=vmax)
    ax3.set_title(f'Spin Splitting (UP - DW) - B{band_index}', fontsize=12)
    ax3.set_xlabel('$k_x$ (1/Å)', fontsize=10)
    ax3.set_ylabel('$k_y$ (1/Å)', fontsize=10)
    ax3.set_aspect('equal')
    cbar3 = plt.colorbar(im3, ax=ax3, label='ΔE (eV)')
    
    # 4. 自旋劈裂绝对值
    ax4 = axes_2d[1, 1]
    im4 = ax4.pcolormesh(kx_2d, ky_2d, np.abs(spin_splitting), 
                        cmap='hot', shading='auto')
    ax4.set_title(f'Absolute Spin Splitting - B{band_index}', fontsize=12)
    ax4.set_xlabel('$k_x$ (1/Å)', fontsize=10)
    ax4.set_ylabel('$k_y$ (1/Å)', fontsize=10)
    ax4.set_aspect('equal')
    plt.colorbar(im4, ax=ax4, label='|ΔE| (eV)')
    
    plt.tight_layout()
    # 生成带有能带编号的输出文件名
    output_2d = f'spin_splitting_2d_maps_B{band_index}.png'
    plt.savefig(output_2d, dpi=600, bbox_inches='tight')
    print(f"2D映射图已保存为: {output_2d}")
    
    # =============================================
    # 6. 绘制3D曲面图 - 修改z轴标签位置
    # =============================================
    
    print("\n" + "=" * 60)
    print("绘制3D曲面图...")
    
    fig_3d = plt.figure(figsize=(14, 10))
    fig_3d.suptitle(f'3D Spin Splitting Surface (UP - DW) - Band B{band_index}', fontsize=16, fontweight='bold')
    
    # 创建3D子图
    # 1. 自旋劈裂3D曲面
    ax3d_1 = fig_3d.add_subplot(221, projection='3d')
    surf1 = ax3d_1.plot_surface(kx_2d, ky_2d, spin_splitting, 
                               cmap='RdBu_r', linewidth=0, 
                               antialiased=True, alpha=0.9, rstride=5, cstride=5)
    ax3d_1.set_xlabel('$k_x$ (1/Å)', fontsize=10)
    ax3d_1.set_ylabel('$k_y$ (1/Å)', fontsize=10)
    # 修改：增加z轴标签与坐标轴的距离，避免与数字标签重叠
    ax3d_1.set_zlabel('ΔE (eV)', fontsize=10, labelpad=15)
    ax3d_1.set_title(f'Spin Splitting Surface - B{band_index}', fontsize=12)
    fig_3d.colorbar(surf1, ax=ax3d_1, shrink=0.5, aspect=10, label='ΔE (eV)', pad=0.12)
    
    # 2. 自旋劈裂3D线框图
    ax3d_2 = fig_3d.add_subplot(222, projection='3d')
    wire1 = ax3d_2.plot_wireframe(kx_2d, ky_2d, spin_splitting, 
                                 rstride=10, cstride=10, 
                                 color='blue', linewidth=0.5, alpha=0.7)
    ax3d_2.set_xlabel('$k_x$ (1/Å)', fontsize=10)
    ax3d_2.set_ylabel('$k_y$ (1/Å)', fontsize=10)
    # 修改：增加z轴标签与坐标轴的距离
    ax3d_2.set_zlabel('ΔE (eV)', fontsize=10, labelpad=15)
    ax3d_2.set_title(f'Spin Splitting Wireframe - B{band_index}', fontsize=12)
    
    # 3. UP能带3D曲面
    ax3d_3 = fig_3d.add_subplot(223, projection='3d')
    surf3 = ax3d_3.plot_surface(kx_2d, ky_2d, band_up_2d, 
                               cmap='viridis', linewidth=0, 
                               antialiased=True, alpha=0.9, rstride=5, cstride=5)
    ax3d_3.set_xlabel('$k_x$ (1/Å)', fontsize=10)
    ax3d_3.set_ylabel('$k_y$ (1/Å)', fontsize=10)
    # 修改：增加z轴标签与坐标轴的距离
    ax3d_3.set_zlabel('E (eV)', fontsize=10, labelpad=15)
    ax3d_3.set_title(f'Spin-up Band Surface - B{band_index}', fontsize=12)
    fig_3d.colorbar(surf3, ax=ax3d_3, shrink=0.5, aspect=10, label='Energy (eV)', pad=0.12)
    
    # 4. DW能带3D曲面
    ax3d_4 = fig_3d.add_subplot(224, projection='3d')
    surf4 = ax3d_4.plot_surface(kx_2d, ky_2d, band_dw_2d, 
                               cmap='viridis', linewidth=0, 
                               antialiased=True, alpha=0.9, rstride=5, cstride=5)
    ax3d_4.set_xlabel('$k_x$ (1/Å)', fontsize=10)
    ax3d_4.set_ylabel('$k_y$ (1/Å)', fontsize=10)
    # 修改：增加z轴标签与坐标轴的距离
    ax3d_4.set_zlabel('E (eV)', fontsize=10, labelpad=15)
    ax3d_4.set_title(f'Spin-down Band Surface - B{band_index}', fontsize=12)
    fig_3d.colorbar(surf4, ax=ax3d_4, shrink=0.5, aspect=10, label='Energy (eV)', pad=0.12)
    
    plt.tight_layout()
    # 生成带有能带编号的输出文件名
    output_3d = f'spin_splitting_3d_surfaces_B{band_index}.png'
    plt.savefig(output_3d, dpi=600, bbox_inches='tight')
    print(f"3D曲面图已保存为: {output_3d}")
    
    # =============================================
    # 7. 绘制统计分布图
    # =============================================
    
    print("\n" + "=" * 60)
    print("绘制统计分布图...")
    
    fig_stats, (ax_hist, ax_scatter) = plt.subplots(1, 2, figsize=(14, 6))
    
    # 直方图
    spin_splitting_flat = spin_splitting.flatten()
    spin_splitting_flat = spin_splitting_flat[~np.isnan(spin_splitting_flat)]
    
    ax_hist.hist(spin_splitting_flat, bins=100, density=True, 
                alpha=0.7, color='steelblue', edgecolor='black')
    ax_hist.axvline(x=0, color='red', linestyle='--', linewidth=1.5, label='ΔE=0')
    ax_hist.set_xlabel('Spin Splitting ΔE (eV)', fontsize=12)
    ax_hist.set_ylabel('Probability Density', fontsize=12)
    ax_hist.set_title(f'Distribution of Spin Splitting Values - Band B{band_index}', fontsize=14)
    ax_hist.legend()
    ax_hist.grid(True, alpha=0.3)
    
    # 散点密度图
    scatter = ax_scatter.scatter(kx_2d.flatten(), ky_2d.flatten(), 
                               c=spin_splitting.flatten(), 
                               cmap='RdBu_r', s=1, alpha=0.7,
                               vmin=-vmax, vmax=vmax)
    ax_scatter.set_xlabel('$k_x$ (1/Å)', fontsize=12)
    ax_scatter.set_ylabel('$k_y$ (1/Å)', fontsize=12)
    ax_scatter.set_title(f'Scatter Density of Spin Splitting - Band B{band_index}', fontsize=14)
    ax_scatter.set_aspect('equal')
    plt.colorbar(scatter, ax=ax_scatter, label='ΔE (eV)')
    
    plt.tight_layout()
    # 生成带有能带编号的输出文件名
    output_stats = f'spin_splitting_statistics_B{band_index}.png'
    plt.savefig(output_stats, dpi=600, bbox_inches='tight')
    print(f"统计分布图已保存为: {output_stats}")
    
    # =============================================
    # 8. 保存数据
    # =============================================
    
    print("\n" + "=" * 60)
    print("保存数据文件...")
    
    # 保存自旋劈裂数据
    output_data = np.column_stack((
        kx_2d.flatten(),
        ky_2d.flatten(),
        band_up_2d.flatten(),
        band_dw_2d.flatten(),
        spin_splitting.flatten()
    ))
    
    # 生成带有能带编号的输出文件名
    data_filename = f'spin_splitting_data_B{band_index}.txt'
    np.savetxt(data_filename, output_data,
               header='kx(1/A) ky(1/A) E_up(eV) E_dw(eV) Delta_E(eV)',
               fmt='%.8f', comments='')
    
    print(f"数据已保存为: {data_filename}")
    
    # 保存统计摘要
    summary_filename = f'spin_splitting_summary_B{band_index}.txt'
    with open(summary_filename, 'w') as f:
        f.write("=" * 60 + "\n")
        f.write(f"Spin Splitting Analysis Summary - Band B{band_index}\n")
        f.write("=" * 60 + "\n\n")
        f.write(f"Grid dimensions: {kx_2d.shape}\n")
        f.write(f"Total k-points: {kx_2d.size}\n\n")
        f.write("Statistics:\n")
        f.write(f"  Minimum ΔE: {np.nanmin(spin_splitting):.6f} eV\n")
        f.write(f"  Maximum ΔE: {np.nanmax(spin_splitting):.6f} eV\n")
        f.write(f"  Mean ΔE: {np.nanmean(spin_splitting):.6f} eV\n")
        f.write(f"  Std Dev ΔE: {np.nanstd(spin_splitting):.6f} eV\n")
        f.write(f"  Mean |ΔE|: {np.nanmean(np.abs(spin_splitting)):.6f} eV\n")
        f.write(f"\nFiles generated:\n")
        f.write(f"  1. {output_2d}\n")
        f.write(f"  2. {output_3d}\n")
        f.write(f"  3. {output_stats}\n")
        f.write(f"  4. {data_filename}\n")
        f.write(f"  5. {summary_filename}\n")
    
    print(f"统计摘要已保存为: {summary_filename}")
    
    # =============================================
    # 9. 显示图形
    # =============================================
    
    print("\n" + "=" * 60)
    print("绘制完成! 显示所有图形...")
    print("=" * 60)
    
    plt.show()

if __name__ == "__main__":
    main()

2.3 wannier90+wanniertools

发表于 2023-05-09 更新于 2026-01-09 分类于 DFT学习

参考：

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

https://mp.weixin.qq.com/s?__biz=MzI3MzIyMjI2Nw==&mid=2247488523&idx=1&sn=90a544712f166e248993e1198eb40981&chksm=eb27cb1fdc50420921bb11d3bc8f4833c14418b0d0c0d755af6687fca024fc9d15d999ae7813&scene=21#wechat_redirect

https://blog.sciencenet.cn/blog-2909108-1154273.html

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

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

一、安装编译

wannier90安装

下载: https://wannier.org/download/

2.1版本是可以与VASP结合的最新版本。但是2.1版本默认安装与VASP接口并不好，主要是借助肖承诚博士写了一个Fortran的接口（https://github.com/Chengcheng-Xiao/VASP2WAN90_v2_fix），需要注意的是这个接口是针对VASP 5.4.4版本的。1.2版本与VASP的接口是好的，所以1.2版本默认安装就好。对比两个版本的不同，1.2版本的主要缺点就在于不能构建向上自旋和向下自旋的能带，也就是没有自旋轨道耦合作用的铁磁和反铁磁体系。 0201

1.1 wannier90 2.1安装

首先，解压安装包

tar -zxvf wannier90-2.1.0.tar.gz

其次，进入文件夹

cd wannier90-2.1.0/

然后，准备编译文件（这里老王用的是ifort，注意要检查ifort和mpiifort执行命令）

cp config/make.inc.ifort make.inc

接着，编译

make

完成后，编译库，得到libwannier.a文件

make lib

1.2 VASP编译（注意以上接口VASP2WAN90_v2_fix只针对VASP 5.4.4）

首先拷贝VASP2WAN90_v2_fix接口文件中的mlwf.patch 到VASP代码编译目录下。然后执行如下命令。

patch -p0 < mlwf.patch

接着在VASP makefile.include 文件中加入下面两行。注意路径。

1
2
3

CPP_OPTIONS+=-DVASP2WANNIER90v2

LLIBS+=/path/to/your/wannier90_distro/libwannier.a

最后编译即可

make all

wannier-tools安装

开源软件包WannierTools，是一个用于研究新型拓扑材料的软件。此代码在紧束缚模型中工作，可以由另一个软件包Wannier90生成。

它可以通过计算Wilson loop来帮助对给定材料的拓扑相进行分类，通过角分辨光电发射(ARPES)和红外光谱技术，可以得到表面态光谱扫描隧道显微镜(STM)实验。它还可以识别Weyl/Dirac点和节点线结构的位置，计算闭合动量环周围的贝里相位和部分布里渊区(BZ)的贝里曲率。
此外，WannierTools还可以计算非磁性金属和半金属的普通磁电阻利用玻尔兹曼输运理论，计算给定磁场方向和强度下的朗道能级谱;并从超胞计算中得到展开的能谱。

wanniertools官方安装教程

http://www.wanniertools.com/installation.html

安装依赖环境

1、Intel编译器即oneapi（Fortran compiler (Gfortran or ifort)，MPICH version higher than 2.1.5，Lapack and Blas library, (Intel MKL recommended)）

2、Arpack-ng

安装步骤

下载wanniertools和arpack安装包并分别解压

1 2	unzip wannier_tools-master.zip tar vzxf arpack96.tar.gz

先进入arpack，修改编译文件并编译
在ARmake.inc中修改编译路径，这里以实际安装路径为准

1	home = $(HOME)/ARPACK

修改编译器，参考VASP安装时的makefile.include

1 2	FC=f77 FFLAGS= -O -cg89

改为

1
2
3

FC         = mpiifort
FCL        = mpiifort -mkl=sequential
FFLAGS  = -O -cg89

然后执行make编译

make lib

直至产生libarpack_SUN4.a文件，该文件会因为不同机器不同环境而不同名称，只要确定好，然后在编译wanniertools的时候指定好就行了

在wanniertools的src文件夹的makefile文件中该默认值为

1
2
3

#ARPACK LIBRARY

ARPACK=/Users/quan/work/workplace/ARPACK/libarpack MAC.a

修改为刚刚已经生成的文件路径，这里使用的是Makefile.intel-mpi，重命名为makefile

1	cp Makefile.intel-mpi Makefile

然后编译即可

make

最后生成wt.x文件即为安装成功

1	cp -f wt.x../bin

二、计算过程(以石墨烯为例)

1、结构优化加自洽计算

POSCAR为

{.line-numbers}

Surf-mp-48                             
   1.00000000000000     
     1.2344083195740001   -2.1380573715510001    0.0000000000000000
     1.2344083195740001    2.1380573715510001    0.0000000000000000
     0.0000000000000000    0.0000000000000000   12.0000000000000000
   C 
     2
Direct
  0.0000000000000000  0.0000000000000000  0.5000000000000000
  0.3333329999917112  0.6666670000099160  0.5000000000000000

INCAR为

{.line-numbers}

SYSTEM=Graphene
Global Parameters
ISTART =  0            (Read existing wavefunction; if there)
ISPIN  =  1            (Non-Spin polarised DFT)
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   =  A            (Precision level)
ADDGRID= .TRUE.        (Increase grid; helps GGA convergence)
 
Electronic Relaxation
ALGO = N
ISMEAR =  0            (Gaussian smearing; metals:1)
SIGMA  =  0.05         (Smearing value in eV; metals:0.2)
NELM   =  90           (Max electronic SCF steps)
NELMIN =  6            (Min electronic SCF steps)
EDIFF  =  1E-08        (SCF energy convergence; in eV)
 
Ionic Relaxation
NSW    =  0            (Max ionic steps)
IBRION = -1            (Algorithm: 0-MD; 1-Quasi-New; 2-CG)
ISIF   =  3            (Stress/relaxation: 2-Ions, 3-Shape/Ions/V, 4-Shape/Ions)
EDIFFG = -2E-02        (Ionic convergence; eV/AA)

加SOC在自洽计算加上
LSORBIT=T

2、非自洽计算投影能带

在自洽计算基础上，cp -r scf band
修改INCAR

1
2
3

ISTART =  1         
ICHARG =  11 
LORBIT =  10/11

KPOINTS 改为高对称点路径。

画一下能带结构，然后再做一下Fatband分析，了解费米面附近能带是有哪些成分构成的。以C的pz轨道。

3、非自洽vasp产生.win文件

在自洽计算基础上 cp -r scf wan1
修改INCAR

ISTART =  1         
ICHARG =  11 
NBANDS = 40
LWANNIER90 = .TRUE.
NUM_WANN = 2

注意，INCAR中使用NPAR参数会出错

KPOINTS、POSCAR、POTCAR不改。

准备好wannier90.win文件，同时在INCAR里头添加一行LWANNIER90=.TRUE.和一行ICHARG=11，然后再运行一遍VASP。顺利的话，会产生 wannier90.amn, wannier90.mmn, wannier90.eig, wannier90.chk 等4个文件。wannier90.win里头的内容也有可能会被改变。包括结构信息，k点信息。有了这三个文件就可以进一步调节参数，构造wannier函数了。请注意wannier90.win里头的num_bands要和OUTCAR里头的NBANDS一致。.mmn积分矩阵和.amn布洛赫态在轨道上的尝试投影和UNK0000x.1存储实空间单胞的布洛赫态相关信息，这些文件是运行链接wannier90的DFT计算输出的结果，而gaas.win是计算参数设置文件

这一步的wannier90.win文件只需要加入几个关键参数，其他会自动生成。

# 解纠缠窗口
dis_win_max = 
dis_win_min=
# 冻结窗口
dis_froz_max = 
dis_forz_min = #(default = dis_win_min)

num_bands =    #（default = nbnd - exclude_bands)
num_wann =      #(num_wann <= num_bands, 等于投影总轨道数，此处为32）
exclude_bands =    

begin projections   ## block
...
end projections

mp_grid = 8 8 8  ## 与自洽的KPOINTS一致

4、wannier90.x计算

运行wannier90主程序。mpirun -np 40 wannier90.x wannier90 &

注意在第3步打开hr_plot = T，关闭write_hr = true；在第4步打开write_hr = true，关闭hr_plot = T。

这里后面一个winner90代表wannier90.win的文件.win前面的部分。如果你的win文件是name.win运行的地方就可以是

wannier90.x name &

运行过程中当num_bands > num_wann 有个迭代过程，控制迭代步数可以通过

1 2	dis_num_iter = 200 # 迭代1，默认200步 num_iter = 200 # 迭代2，默认200步

最后会生成wannier90.chk文件，restart计算会读取这个文件。

加入能带部分的话第一次运行就会出来，wannier90_band.dat文件，这个可以用来和DFT结果对比。
构造好Wannier函数以后，仔细检查一下wannier90.wout，
搜索”Final State”，在这里可以看到每条Wannier轨道的展宽，由此可以检查你构造的Wannier函数是否足够局域。如果有某个Wannier函数的展宽
比晶格常数大许多，就表明你的Wannier函数构造不是很理想，需要不断调节投影轨道，解纠缠窗口以及Frozen窗口等。从这些轨道的中心和展宽也能
看出轨道之间的简并性如何，由此也可以判断所构造的Wannier函数对称性如何，这个也会影响后续的分析，特别是表面态的计算。查看wannier90.wout末尾的Final State的最后一列小于3，其结果是可以的。

5、

紧接着，你就可以开始书写WannierTools的主输入文件了，
目前叫input.dat，以后会改成http://wt.in (WannierTools v2.2之后)。这里不具体介绍每一个参数，在说明文档中有详细的介绍，中文版的说明文档正在翻译中。
input.dat文件有模板，使用的时候直接拷贝然后加以修改。很多内容都可以直接从wannier90.win和wannier90.wout中拷
贝。 Bi2Se3的主输入文件input.dat可以在附件中下载得到。我们可以通过控制输入文件中的一些参数来实现我们所需要的功能。

常见报错

首先跑接口会生成3个文件.amm .mmm .eig通过这三个文件构造wannier函数。也可以提前把.win文件准备好放进去。以上三个文件生成后运行wannier90.x会成.chk文件
Restar=T是需要读取.chk文件的

Error: restart requested but wannier90.chkfile not found

原因：加了restart = plot 这个命令，但却没有wannier90.chk这个文件。Restart是需要读取wannier90.chk的。

解决办法：注释掉restart = plot，然后重新运行wannier90.x wannier90 &得到wannier90.chk.

问题：Wannier90 is running in LIBRARY MODE

Setting up k-point neighbours…

Ignoring in input file

解决办法：尽量串行，不要并行。请再次运行 $postw90 wannier90的为输出的数据

Spinors与文件INCAR中的ISPIN的设置有关，同时存在或者同时去掉，否则影响.amn文件的形成。（自旋轨道耦合）
谨慎用use_bloch_phases=T一般是F
问题：param_read: mismatch in wannier90.eig

解决办法：每一步计算的核数要相等

问题：dis_windows: More states in the forze window than target WFs

解决办法：把dis_windows改大点

Fewer projections defined than the number of wannier functions requested.

投影的轨道太少

解决办法：增加投影的轨道

Kmesh_get_bvector: Not enongh bevectors found

解决办法：添加kemsh_tol

Wanted band:1 found band:17

Wanted kpoint:2 found kpoint:1

解决办法：常见原因是num_bands设置错误，返回band计算文件夹查看NBANDS

Dis_windows: Energy window contains fewer states than number sf target WFs

解决办法：把dis_windows改大点

Dis_windows: more states in the frozen window than WFs

解决办法：调节frozen window, 投影太多或者这个窗口太大涉及到别的轨道

Wannier 运行有两种模式一是library mode，一种是单机模式stand-alone mode

解决办法：library mode是需要进行串行。

param_get_projection: too many projections defined

解决办法：begin projections random

原因：num_wann 与设置投影数不匹配。

解决办法：不加soc：num_wann 等于每种原子投影轨道数之和

soc：num_wann 等于两倍的每种原子投影轨道数之和

Gnuplot的报错信息”arc_bulk.gnu”, line 26: warning: No unsable data in this plot to auto-scale axis range
param_get_projection: Problem reading m state into string

解决办法：已解决Fe: s;px;py;pz;dz2;dxz;dyz;dx2-y2;dxy

不是Fe:s;p;d

能带拟合的不好可能是因为，投影的不对，还有核数和能带数的关系

解决办法：太马虎了，这样的问题肯定是K点写错了

能带布里渊区的宽度不一样有可能是wannier里规定的K点数并没有完全计算。

num_iter = 2000 #wannier center 迭代次数

num_print_cycles = 40 #wannier center 迭代每40步输出一次

dis_num_iter = 1200 #解纠缠的迭代步数#
num_iter最小化过程中的总迭代次数。如果您希望生成投影wf而不是最大本地化的wf，则设置num_iter=0(请参阅教程中的示例8)？默认值为100
dis_mix_ratio = 1.d0 # 0.5<dis_mix_ratio<1, 默认值0.5

#在解纠缠过程中，要利用混合参数进行收敛，值为0.5是“安全”的选择。使用1.0(即不混合)通常会得到更快的收敛，但可能会导致最小化，在某些情况下是不稳定的。

dis_conv_tol = 1.0E-10 # Delta < 1.000E-10 over 3 iterations

dis_win_max = 15

dis_froz_max = 14

num_wan = 12

exclude_bands : 1-11,26 #在使用这个参数时，num_bands = 总能带数减去排除在外的能带数

num_bands = 14

#如果使用能量窗口规定想要拟合的能带区间，需要同时使用这个几个参数

wannier_plot_supercell =2 #在一个对应于“超胞”的网格上生成WFs
wannier_plot_list =2 #列举要画出WFs的序号，与wannier center 的序号是一致的
运行接口的情况有三种，一是选择投影能带数发生了改变，二是，需要写入wannier90.win的输入文件；三是需要在INCAR添加关键参数，输出wannier.x需要的文件。（仅仅是调节能量窗口是不需要重新运行接口的）
无没有projection是肯定没有amn的

wannier90.win参数介绍

系统参数

num_wann： I， Number of WF 
num_bands： I， Number of bands passed to the code
unit_cell_cart： P， Unit cell vectors in Cartesian coordinates，不用设默认即可
atoms_cart*： P， Positions of atoms in Cartesian coordinates，会自动补全 
atoms_frac*： R， Positions of atoms in fractional coordinates with respect to the lattice vectors会自动补全
mp_grid： I，Dimensions of the Monkhorst-Pack grid of k-points与自洽KPOINTS一致
spinors： L， WF are spinors考虑SOC或者不考虑，默认是关闭不考虑

任务控制

exclude_bands： I， List of bands to exclude from the calculation不考虑那几条能带
restart： S， Restart from checkpoint file 
iprint： I， Output verbosity level
spin： S， Which spin channel to read
translate_home_cell： L， To translate final Wannier centres to home unit cell when writing xyz file 
write_xyz： L， To write atomic positions and final centres in xyz file format

画图参数

wannier_plot： L， Plot the WF
bands_plot： L， Plot interpolated band structure 
kpoint_path： P， K-point path for the interpolated band structure 
bands_num_points： I， Number of points along the first section of the k-point path 
bands_plot_format： S， File format in which to plot the interpolated 
bands bands_plot_project： I， WF to project the band structure onto 
bands_plot_mode： S， Slater-Koster type interpolation or Hamiltonian cut-off
fermi_surface_plot： L， Plot the Fermi surface
write_hr： L， Write the Hamiltonian in the WF basis
dist_cutoff： P， Cut-off for the distance between WF

解纠缠参数

dis_win_min： P， Bottom of the outer energy window 
dis_win_max： P， Top of the outer energy window 
dis_froz_min： P， Bottom of the inner (frozen) energy window 
dis_froz_max： P， Top of the inner (frozen) energy window 
dis_num_iter： I， Number of iterations for the minimisation of ΩI

wannier化参数

1
2
3

num_iter： I， Number of iterations for the minimisation of Ω
num_print_cycles： I， Control frequency of printing
guiding_centres： L， Use guiding centres

wt.in参数介绍

可参考：https://mp.weixin.qq.com/s/HUDobUio4BdQfVLImZfYmg

先给一个模版：

&TB_FILE
Hrfile = 'wannier90_hr.dat'  !default
Package = 'VASP' !default
/


&CONTROL
AHC_calc = T
ANE_calc = F
BulkBand_calc         = T
BulkFS_calc           = F
BulkGap_cube_calc     = F
BulkGap_plane_calc    = F
SlabBand_calc         = F
WireBand_calc         = F
SlabSS_calc           = F
SlabArc_calc          = F
SlabQPI_calc          = F
Z2_3D_calc            = F
Chern_3D_calc         = F
SlabSpintexture_calc  = F
Wanniercenter_calc    = T
BerryPhase_calc       = T
BerryCurvature_calc   = F
EffectiveMass_calc    = F
/

&SYSTEM
#NSLAB = 10              ! for thin film system default=10
#NSLAB1= 4               ! nanowire system
#NSLAB2= 4               ! nanowire system
NumOccupied = 18        ! NumOccupied=num_wan
SOC = 1                 ! soc on
E_FERMI = 4.4195        ! e-fermi
#Bx= 0, By= 0, Bz= 0     
#surf_onsite= 0.0        
/

&PARAMETERS
Eta_Arc = 0.0001     ! infinite small value, like brodening
#E_arc = -0.12         ! energy for calculate Fermi Arc
OmegaNum = 1001      ! omega number step=1meV
OmegaMin = -1     ! energy interval
OmegaMax =  1     ! energy interval
Nk1 = 151           ! number k points  odd number would be better
Nk2 = 151           ! number k points  odd number would be better
Nk3 = 1           ! number k points  odd number would be better test
#NP = 1              ! number of principle layers default=2
#Gap_threshold = 1.0 ! threshold for GapCube output
TMin = 100     ! energy interval 注意 Tmin ≠0
TMax = 300    ! energy interval
NumT = 100  ! omega number  （300-100）/100 从100开始间隔100K输出一次
/

LATTICE
Angstrom
-2.069  -3.583614  0.000000     ! crystal lattice information
 2.069  -3.583614  0.000000
 0.000   2.389075  9.546667

ATOM_POSITIONS
5                               ! number of atoms for projectors
Direct                          ! Direct or Cartisen coordinate
 Bi 0.3990    0.3990    0.6970
 Bi 0.6010    0.6010    0.3030
 Se 0     0     0.5
 Se 0.2060    0.2060    0.1180
 Se 0.7940    0.7940    0.8820

PROJECTORS
 3 3 3 3 3          ! number of projectors
Bi px py pz         ! projectors
Bi px py pz
Se px py pz
Se px py pz
Se px py pz

SURFACE            ! from version v2.2.6 on, you can specify a surface with only two lattice vectors
 1  0  0
 0  1  0

KPATH_BULK            ! BulkBand_calc=T
4              ! number of k line only for bulk band
G 0.00000 0.00000 0.0000 Z 0.00000 0.00000 0.5000
Z 0.00000 0.00000 0.5000 F 0.50000 0.50000 0.0000
F 0.50000 0.50000 0.0000 G 0.00000 0.00000 0.0000
G 0.00000 0.00000 0.0000 L 0.50000 0.00000 0.0000

KPATH_SLAB       !SlabBand_calc=T or SlabSS_calc=T
2        ! numker of k line 
K 0.33 0.67 G 0.0 0.0  
G 0.0 0.0 M 0.5 0.5

KPLANE_SLAB      !SlabArc_calc=T, SlabSpintexture_calc= T
-0.1 -0.1      ! Original point for 2D k plane
 0.2  0.0      ! The first vector to define 2D k plane
 0.0  0.2      ! The second vector to define 2D k plane  for arc plots

KPLANE_BULK     !BulkGap_plane_calc=T, BerryCurvature_calc=T, wanniercenter_calc=T
 0.00  0.00  0.50   ! Original point for 3D k plane  k3=0.5, bar{a}, along k1
 1.00  0.00  0.00   ! The first vector to define 3d k space plane
 0.00  0.50  0.00   ! The second vector to define 3d k space plane

!> The following 5 matrices are for backup using, will not affect the main input for WannierTools
 0.00  0.00  0.00   ! Original point for 3D k plane  k3=0.0, bar{a}, along k1
 1.00  0.00  0.00   ! The first vector to define 3d k space plane
 0.00  0.50  0.00   ! The second vector to define 3d k space plane

 0.00  0.50  0.00   ! Original point for 3D k plane  k2=0.5, bar{c}, along ka
 0.00  0.00  1.00   ! The first vector to define 3d k space plane
 0.50  0.00  0.00   ! The second vector to define 3d k space plane

 0.00  0.00  0.00   ! Original point for 3D k plane  k2=0, bar{c}, along ka
 0.00  0.00  1.00   ! The first vector to define 3d k space plane
 0.50  0.00  0.00   ! The second vector to define 3d k space plane

 0.50  0.00  0.00   ! Original point for 3D k plane  k1=0.5, bar{c}, along kb
 0.00  0.00  1.00   ! The first vector to define 3d k space plane
 0.00  0.50  0.00   ! The second vector to define 3d k space plane

 0.00  0.00  0.00   ! Original point for 3D k plane  k1=0, bar{c}, along kb
 0.00  0.00  1.00   ! The first vector to define 3d k space plane
 0.00  0.50  0.00   ! The second vector to define 3d k space plane



KCUBE_BULK      !BulkGap_cube_calc=T
-0.50 -0.50 -0.50   ! Original point for 3D k plane
 1.00  0.00  0.00   ! The first vector to define 3d k space plane
 0.00  1.00  0.00   ! The second vector to define 3d k space plane
 0.00  0.00  1.00   ! The third vector to define 3d k cube


EFFECTIVE_MASS      ! optional
2                   ! The i'th band to be calculated
0.01                ! k step in unit of (1/Angstrom)
0.0 0.0 0.0         ! k point where the effective mass calculated.


WANNIER_CENTRES     ! copy from wannier90.wout
Cartesian
 -0.000040  -1.194745   6.638646
  0.000038  -1.196699   6.640059
 -0.000032  -1.192363   6.640243
 -0.000086  -3.583414   2.908040
  0.000047  -3.581457   2.906587
 -0.000033  -3.585864   2.906443
 -0.000001   1.194527   4.773338
  0.000003   1.194538   4.773336
 -0.000037   1.194536   4.773327
  0.000006  -1.194384   1.130261
 -0.000018  -1.216986   1.140267
  0.000007  -1.172216   1.140684
  0.000011  -3.583770   8.416406
 -0.000002  -3.561169   8.406398
 -0.000007  -3.605960   8.405979
  0.000086  -1.194737   6.638626
 -0.000047  -1.196693   6.640080
  0.000033  -1.192286   6.640223
  0.000040  -3.583406   2.908021
 -0.000038  -3.581452   2.906608
  0.000032  -3.585788   2.906424
  0.000001   1.194548   4.773330
 -0.000003   1.194537   4.773332
  0.000037   1.194539   4.773340
 -0.000011  -1.194381   1.130260
  0.000002  -1.216981   1.140268
  0.000007  -1.172191   1.140687
 -0.000006  -3.583766   8.416405
  0.000018  -3.561165   8.406400
 -0.000007  -3.605935   8.405982

Tips:

2D材料就按照bulk去算
需要测试Nk1 、NK2、NK3对结果的影响
必须设置的模块：&TB_FILE、&CONTROL、&SYSTEM、&PARAMETERS、LATTICE、ATOM_POSITIONS、PROJECTORS、WANNIER_CENTRES ；选择设置SURFACE、KPATH_BULK、 KPATH_SLAB、 KPLANE_SLAB、KPLANE_BULK 、KCUBE_BULK、EFFECTIVE_MASS 等。