能带的计算以及origin绘制能带图

in 从零开始的科研之旅, VASP

终于也是把文献整理得差不多了，博客也翻新了一下，继续分享一些之前的计算。

我的计算体系是Argyrodite硫化物固态电解质体系，计算以$\rm Li_6PS_5Cl$为例吧。能带计算最好是使用primitive cell，即原胞。这是惯用晶胞或者超胞计算得到的能带是折叠的，固体物理的相关书籍中有提到，（当然计算成本也是不同的），参考公社论坛上的两篇讨论帖子：

  • http://bbs.keinsci.com/thread-11373-1-1.html
  • http://bbs.keinsci.com/thread-19361-1-1.html

原胞的变胞优化

计算使用我稍微熟一点的VASP来，首先从Materials Project网站下载晶体结构文件，下载primitive cell，如果已经有了conventional cell的cif文件，也可以使用vaspkit的602功能生成PRIMCELL.vasp原胞文件。

然后进行原胞的变胞优化，除了POSCAR的剩下三个输入文件可以使用vaspkit方便地生成，使用101功能，选择lattice relaxation进行变胞优化，102功能生成自洽计算地K点，选择Gamma撒点方案，0.03间隙，原胞计算量比较小，间隙填小一点也无妨。

笔者在Intel 5218R的集群上进行计算，分配了40个物理核心，计算时在INCAR中设置KPAR=2和NCORE=20来最大化效率，一般而言NPAR参数不需要设置，KPAR和NCORE的参数设置众说纷纭，我参考的公社某一篇帖子里的经验，前者设置为2，后者设置为物理核心数的一般，当然这是对于CPU版本的VASP而言的。

原胞的变胞优化用时很短：

1
2
3
4
5
6
7
8
9
10
11
12
 reached required accuracy - stopping structural energy minimisation
[ctan@baifq-hpc141 primitive_cell-opt]$ tail OUTCAR
                            User time (sec):       88.516
                          System time (sec):        3.010
                         Elapsed time (sec):       94.030

                   Maximum memory used (kb):      238296.
                   Average memory used (kb):          N/A

                          Minor page faults:        36126
                          Major page faults:          669
                 Voluntary context switches:         1271

在确认收敛并确认结构正常后进行能带的计算

PBE泛函计算能带

能带计算的KPOINTS与普通计算的KPOINTS不一样，通常需要第一布里渊区内的一条或几条高对称点路径来计算能带性质

VASP的计算中PBE泛函是被使用最广泛的，需要修改K点和输入文件。参考vaspkit官网的教程：http://vaspkit.cn/index.php/27.html

PBE能带的计算可以分两步进行。

首先进行自洽计算，此时也即（最常用的）静态计算，将变胞优化的CONTCAR更名为POSCAR，使用vaspkit的101功能，选择Static-Calculation，k点选择可以不变。

接着进行能带计算，需要使用vaspkit的303功能生成高对称K点，将生成的KPATH.in文件更名为KPOINTS

1
2
(base) [storm@X16 temp]$ ls
HIGH_SYMMETRY_POINTS  INCAR  KPATH.in  KPOINTS  POSCAR  POTCAR  PRIMCELL.vasp  SYMMETRY  vmd-1.9.3

可以看到，303功能也生成了原胞结构，当然前面我们已经提到过。

在进行能带计算的目录中，将上一步自洽计算的CHGCAR文件拷贝过来并在INCAR中设置电荷密度不更新，即ICHARG=11，这是必须的，为了加速计算，建议读取上一步的WAVECAR文件，虽然ISTART设置为1或者2程序都可以读取WAVECAR，但是为了保持整个计算任务中（第一步和第二步）的平面波基组一致，建议使用ISTART=2。

此种方法计算能带分两步的原因是使用高密度 k 点网格进行计算时耗时较长，因此用不耗时的普通静态计算多次迭代得到的CHGCAR和WAVECAR来加速高密度k点网格的计算。

其实对于我的小体系，直接使用高密度k点网格进行静态计算也是完全可接受的，并且SCF 和能带完全自洽，会更精确一些，用时并不长：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
[ctan@baifq-hpc141 primitive_cell-band_PBE]$ ls
4736.log  CONTCAR   HIGH_SYMMETRY_POINTS  KPOINTS  PCDAT   PRIMCELL.vasp  slurm-4736.out  vasp.out     WAVECAR
CHG       DOSCAR    INCAR                 OSZICAR  POSCAR  PROCAR         SYMMETRY        vaspout.h5   XDATCAR
CHGCAR    EIGENVAL  KPATH.in              OUTCAR   POTCAR  REPORT         vasp.err        vasprun.xml
[ctan@baifq-hpc141 primitive_cell-band_PBE]$ tail OUTCAR
                            User time (sec):      216.131
                          System time (sec):        4.733
                         Elapsed time (sec):      223.789

                   Maximum memory used (kb):      427044.
                   Average memory used (kb):          N/A

                          Minor page faults:        28812
                          Major page faults:          666
                 Voluntary context switches:         1232

HSE06杂化泛函计算能带

PBE最常被诟病的方面即能带计算，它常常低估带隙。By the way，Materials Project上的能带数据刚好能和我用PBE算出来的能带数据对上，明显小于HSE06泛函计算的数据，估计是使用PBE算的。

HSE (Heyd–Scuseria–Ernzerhof)杂化泛函在2003年首次被提出（HSE03），在2006年HSE06被提出，它能更加精确地计算能带。使用HSE06杂化泛函做计算是十分昂贵的，用时是PBE泛函的数十倍不止。

受vasp算法所限制，使用HSE06杂化泛函计算能带时只能使用自洽计算（此时为静态计算），不能进行非自洽计算。中文互联网上有一些教程是错误的，PBE可以自洽+非自洽两步，也可以自洽一步，但是HSE06杂化泛函只能一步自洽。不过可以在能带计算之前也先做一下静态计算（自洽计算），然后将产生的WAVECAR文件用于能带计算以加速。

使用vaspkit的251功能生成所需要的KPOINTS文件，这是一个“混合”的文件：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
 ------------>>
251
 -->> (01) Reading Input Parameters From INCAR File.
 ======================== K-Mesh Scheme ==========================
 1) Monkhorst-Pack Scheme
 2) Gamma Scheme

 0)   Quit
 9)   Back
 ------------->>
2
 +---------------------------- Tip ------------------------------+
 Input the K-spacing value for SCF Calculation:
 (Typical Value: 0.03-0.04 is Generally Precise Enough)
 ------------>>
0.03
 Input the K-spacing value for Band Calculation:
 (Typical Value: 0.03-0.04 for DFT and 0.04-0.06 for hybrid DFT)
 ------------>>
0.03
 +---------------------------------------------------------------+
 -->> (02) Reading K-Path From KPATH.in File.
 +-------------------------- Summary ----------------------------+
 K-Mesh for SCF Calculation:    3    3    3
 The Number of K-Points along K-Path No.1:  14
 The Number of K-Points along K-Path No.2:   6
 The Number of K-Points along K-Path No.3:  18
 The Number of K-Points along K-Path No.4:  17
 The Number of K-Points along K-Path No.5:   7
 The Number of K-Points along K-Path No.6:   7
 +---------------------------------------------------------------+
 -->> (03) Written KPOINTS File.

首先需要生成自洽（静态）计算所需要的k点网格，然后生成能带计算的高密度k点网格。INCAR文件可以沿用普通静态计算的然后改一下泛函，也可以使用101功能生成HSE06 Calculation的INCAR。

HSE06杂化泛函极为昂贵，大概在第五个电子步时，用时超过8h仍未计算完，遂放弃。之前注册某家超算平台时送了一点机时，开通了GPU版本的vasp计算，发现快了不少：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
General timing and accounting informations for this job:
========================================================
 
                 Total CPU time used (sec):    22600.775
                           User time (sec):    10925.947
                         System time (sec):    11674.828
                        Elapsed time (sec):     7781.069
 
                  Maximum memory used (kb):     6086716.
                  Average memory used (kb):          N/A
 
                         Minor page faults:      6831986
                         Major page faults:          282
                Voluntary context switches:     53436870

PROFILE, used timers:     405
=============================

数据后处理与Origin绘图

能带计算结束后，将目录拷贝到可以使用vaspkit的机器上。

对于PBE计算的能带，使用211功能处理数据，得到带隙、CBM、VBM、Fermi Energy等信息：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
 ------------>>
21
 ============================ Band Options =======================
 211) Band-Structure
 212) Projected Band-Structure of Only-One-Selected Atom
 213) Projected Band-Structure of Each Element
 214) Projected Band-Structure of Selected Atoms
 215) Projected Band-Structure by Element-Weights
 216) The Sum of Projected Band for Selected Atoms and Orbitals

 0)   Quit
 9)   Back
 ------------>>
211
 -->> (01) Reading Input Parameters From INCAR File.
 +---------------------------------------------------------------+
 |       >>> The Fermi Energy will be set to zero eV <<<         |
 +---------------------------------------------------------------+
 -->> (02) Reading Fermi-Energy from DOSCAR File.
 -->> (03) Reading Energy-Levels From EIGENVAL File.
 -->> (04) Reading K-Path From KPOINTS File.
 -->> (05) Written KLABELS File.
 +---------------------------- Tip ------------------------------+
 |If You Want to Get Fine Band Structrue by Interpolating Method.|
 | You CAN set GET_INTERPOLATED_DATA = .TRUE. in ~/.vaspkit file.|
 +---------------------------------------------------------------+
 -->> (06) Written BAND.dat File.
 -->> (07) Written REFORMATTED_BAND.dat File.
 -->> (08) Written KLINES.dat File.
 -->> (09) Written BAND_GAP File.

对于HSE06计算的能带，使用252功能处理数据，得到带隙、CBM、VBM、Fermi Energy等信息：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
 ------------>>
252

 -->> (01) Reading Input Parameters From INCAR File.
 +---------------------------------------------------------------+
 |       >>> The Fermi Energy will be set to zero eV <<<         |
 +---------------------------------------------------------------+
 -->> (02) Reading Fermi-Level From FERMI_ENERGY.in File.
 -->> (03) Reading Energy-Levels From EIGENVAL File.
 -->> (04) Reading KPT-Params in the First Line of KPOINTS File.
 -->> (05) Reading K-Path From KPATH.in File.
 -->> (06) Written KLABELS File.
 +---------------------------- Tip ------------------------------+
 |If You Want to Get Fine Band Structrue by Interpolating Method.|
 | You CAN set GET_INTERPOLATED_DATA = .TRUE. in ~/.vaspkit file.|
 +---------------------------------------------------------------+
 -->> (07) Written BAND.dat File.
 -->> (08) Written REFORMATTED_BAND.dat File.
 -->> (09) Written KLINES.dat File.
 -->> (10) Written BAND_GAP File.

绘图时我习惯把价带顶归零，即纵坐标为$\rm E-E_{VBM}$，操作如下。

把FERMI_ENERGY文件就拷贝一份，改名为FERMI_ENERGY.in。将文件FERMI_ENERGY.in中的费米能级数据改为价带顶数据，然后再重新运行252功能，此时得到的REFORMATTED_BAND.dat文件中的数据是归零的。

接下来就是origin绘图了

将REFORMATTED_BAND.dat文件拖入工作表内，确认只有一个X轴，剩下全是Y轴。新建一个sheet。打开KLABELS文本文件，手动复制数据，如下图

sheet2
sheet1

接着选中sheet1表格中的所有数据，绘制折线图：

很丑、、、

设置一下横坐标，改成高对称点的样子，看图：

sheet2第二列
sheet2第一列

这样横坐标就被设置好了：

还是很丑

接下来就是美化了，放一张最后画好的图吧：

Li6PS5Cl能带图

Comment and share

np.asmatrix函数被废弃使得CIF转POSCAR脚本报错

in 从零开始的科研之旅, VASP

使用最新版的Miniconda3，python环境的版本是14.x，顺手pip install numpy安装的numpy版本过高，使用vaspkit尝试将.cif文件转换为POSCAR时报错:

1
AttributeError: `np.mat` was removed in the NumPy 2.0 release. Use `np.asmatrix` instead

错误原因在于当前的NumPy版本太新（>=2.0），而vaspkit的cif2pos.py脚本中使用了已被移除的np.mat函数，最简单的方式是降级numpy版本：

1
pip install numpy==1.26.4

Comment and share

使用Multiwfn加载ELFCAR文件绘制等值面图

in 从零开始的科研之旅, VASP

本文仅为学习记录，并非教程

Multiwfn是一个波函数分析程序，功能非常强大，也能绘制ELF等值面图，它支持多种类型的输入文件，既有量子化学计算程序的输出文件，又有第一性原理计算程序的输出文件

Multiwfn的输入文件(http://sobereva.com/379)

以下是Multiwfn加载ELFCAR绘制等值面图的过程

如果将Multiwfn的根目录添加到PATH环境变量以方便调用，那么画出来的图可能会遇到如下如题

在非根目录下调用Multiwfn

如果直接在根目录下双击.exe文件运行的话就能正常绘图

正常的ELF等值面图

Comment and share

掺杂质量计算的python脚本

in 从零开始的科研之旅, VASP

起因是我想合成各式各样的材料，需要用到各种药品，计算起来比较麻烦，尤其是非化学计量比的化合物比如Li5.3PS4.3Cl0.7Br，遂用AI写了一个python脚本，个人觉得还是比较实用，脚本文件和相对原子质量的文本文件需要在同一文件夹内

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
from decimal import Decimal, getcontext
import re
from sympy import symbols, Matrix, linsolve, Rational
from collections import defaultdict

# Improve Decimal precision
getcontext().prec = 100

class FormulaParser:
    @staticmethod
    def parse(formula: str) -> dict:
        elements = {}
        pattern = re.compile(r"([A-Z][a-z]?)([\d.]*)?")
        idx = 0
        while idx < len(formula):
            m = pattern.match(formula, idx)
            if not m:
                raise ValueError(f"Invalid formula at: {formula[idx:]}")
            elem, num = m.group(1), m.group(2)
            coeff = Decimal(num) if num else Decimal(1)
            elements[elem] = elements.get(elem, Decimal(0)) + coeff
            idx = m.end()
        return elements

class AtomicMassReader:
    @staticmethod
    def read(filename="atomic_masses.txt") -> dict:
        masses = {}
        try:
            with open(filename, 'r') as f:
                for line in f:
                    line = line.strip()
                    if not line or line.startswith('#'):
                        continue
                    parts = line.split()
                    if len(parts) != 2:
                        raise ValueError(f"Line format error: {line}")
                    masses[parts[0]] = Decimal(parts[1])
        except FileNotFoundError:
            raise FileNotFoundError(f"File {filename} not found")
        if not masses:
            raise ValueError("Atomic masses file empty")
        return masses

class SynthesisCalculator:
    def __init__(self, atomic_masses: dict):
        self.atomic_masses = atomic_masses
        self.target = {}
        self.target_mm = Decimal(0)
        self.reagents = []
        self.reagent_formulas = []
        self.molar_masses = []

    def input_target(self):
        f = input("Enter target formula: ").strip()
        self.target = FormulaParser.parse(f)
        for e, c in self.target.items():
            if e not in self.atomic_masses:
                raise ValueError(f"Mass of element {e} not defined")
            self.target_mm += c * self.atomic_masses[e]
        print(f"Parsed: {self.target}")
        print(f"Target molar mass = {self.target_mm.normalize()} g/mol")

    def input_reagents(self):
        n = int(input("Enter number of reagents: "))
        for i in range(n):
            f = input(f"Reagent {i+1} formula: ").strip()
            self.reagent_formulas.append(f)
            parsed = FormulaParser.parse(f)
            mm = sum(parsed[e] * self.atomic_masses[e] for e in parsed)
            self.reagents.append(parsed)
            self.molar_masses.append(mm)
            print(f"{f} molar mass = {mm.normalize()} g/mol")

    def validate(self):
        tgt = set(self.target)
        prov = set().union(*(r.keys() for r in self.reagents))
        miss = tgt - prov
        extra = prov - tgt
        if miss:
            print(f"Error: Missing elements: {miss}")
            exit()
        if extra:
            print(f"Warning: Extra elements: {extra}")

    def solve_basis(self):
        # Construct A, b corresponding to 1 mol 
        elems = sorted(set(self.target) | set().union(*(r.keys() for r in self.reagents)))
        A = [[r.get(e, Decimal(0)) for r in self.reagents] for e in elems]
        b = [self.target.get(e, Decimal(0)) for e in elems]
        # Convert to rational number to avoid floating point error
        A_rat = [[Rational(str(val)) for val in row] for row in A]
        b_rat = [Rational(str(val)) for val in b]
        vars = symbols(f'x0:{len(self.reagents)}')
        sol = linsolve((Matrix(A_rat), Matrix(b_rat)), vars)
        if not sol:
            print("Error: No exact solution for 1 mol target.")
            exit()
        tup = next(iter(sol))
        # Convert rational to decimal exact representation
        base = {}
        for i, r in enumerate(tup):
            base[vars[i]] = Decimal(r.p) / Decimal(r.q)
        return base, elems

    def calculate(self):
        t = input("Input type (0 mass g, 1 moles): ").strip()
        if t not in ['0','1']:
            print("Invalid type")
            exit()
        v = Decimal(input("Amount: ").strip())
        n = v / self.target_mm if t=='0' else v
        print(f"Target moles = {n.normalize()}")

        base, elems = self.solve_basis()
        print("\nReagent requirements:\n" + "="*40)
        total_mass = Decimal(0)
        for i, var in enumerate(sorted(base.keys(), key=lambda x: int(str(x)[1:]))):
            moles = base[var] * n
            mass = moles * self.molar_masses[i]
            total_mass += mass
            print(f"Reagent {i+1} ({self.reagent_formulas[i]}):")
            print(f"  Moles: {moles.normalize():.10f} mol")
            print(f"  Mass:  {mass.normalize():.10f} g\n")
        print(f"Total mass = {total_mass.normalize():.10f} g")

        #Verify element balance
        print("Verification:")
        for e in elems:
            actual = sum((base[symbols(f'x{i}')] * n) * r.get(e, Decimal(0))
                         for i, r in enumerate(self.reagents))
            target_amt = self.target.get(e, Decimal(0)) * n
            print(f"{e}: target {target_amt.normalize():.10f}, actual {actual.normalize():.10f}")


def main():
    print("High-Precision Compound Synthesis Calculator")
    print("="*60)
    try:
        masses = AtomicMassReader.read()
        print(f"Loaded elements: {', '.join(masses.keys())}")
    except Exception as e:
        print(e)
        exit()
    calc = SynthesisCalculator(masses)
    calc.input_target()
    calc.input_reagents()
    calc.validate()
    calc.calculate()

if __name__ == "__main__":
    main()

附一张使用截图，后面也会放到Github上面

使用截图

Comment and share

初学ELF电子定域化函数的计算与分析

in 从零开始的科研之旅, VASP

本文仅为学习记录，并非教程

这次记录文主要侧重于应用，参考了Sob老师的博文和网络上的各大教程

以FCC Cu为例进行ELF计算与分析

采用FCC面心立方最密堆积的铜单胞，使用vasp计算，先进行结构优化，确认收敛后进行ELF计算

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
Global Parameters
ISTART =  1            (Read existing wavefunction, if there)
ISPIN  =  2            (Non-Spin polarised DFT)
MAGMOM = 4*0.6
# ICHARG =  11         (Non-self-consistent: GGA/LDA band structures)
LREAL  = .FALSE.       (Projection operators: automatic)
 ENCUT  =  530        (Cut-off energy for plane wave basis set, in eV)
# PREC   =  Accurate   (Precision level: Normal or Accurate, set Accurate when perform structure lattice relaxation calculation)
LWAVE  = .TRUE.        (Write WAVECAR or not)
LCHARG = .TRUE.        (Write CHGCAR or not)
ADDGRID= .TRUE.        (Increase grid, helps GGA convergence)
LASPH  = .TRUE.        (Give more accurate total energies and band structure calculations)
PREC   = Accurate      (Accurate strictly avoids any aliasing or wrap around errors)
# LVTOT  = .TRUE.      (Write total electrostatic potential into LOCPOT or not)
# LVHAR  = .TRUE.      (Write ionic + Hartree electrostatic potential into LOCPOT or not)
# NELECT =             (No. of electrons: charged cells, be careful)
# LPLANE = .TRUE.      (Real space distribution, supercells)
# NWRITE = 2           (Medium-level output)
# KPAR   = 2           (Divides k-grid into separate groups)
# NGXF    = 300        (FFT grid mesh density for nice charge/potential plots)
# NGYF    = 300        (FFT grid mesh density for nice charge/potential plots)
# NGZF    = 300        (FFT grid mesh density for nice charge/potential plots)

Electronic Relaxation
ISMEAR =  1            (Gaussian smearing, metals:1)
SIGMA  =  0.2         (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)
# GGA  =  PS           (PBEsol exchange-correlation)

Ionic Relaxation
NSW    =  0         (Max ionic steps)
IBRION =  -1            (Algorithm: 0-MD, 1-Quasi-New, 2-CG)
ISIF   =  2            (Stress/relaxation: 2-Ions, 3-Shape/Ions/V, 4-Shape/Ions)
EDIFFG = -2E-02        (Ionic convergence, eV/AA)
# ISYM =  2            (Symmetry: 0=none, 2=GGA, 3=hybrids)
LELF = .TRUE.
NPAR = 1

在优化完的结构上进行自洽计算，添加LELF = .TRUE.这个参数，确保ELFCAR文件输出

官方wiki还提到必须显示开启NPAR = 1，然而网上的各大教程都没有提到这一点

If LELF is set, NPAR=1 has to be set explicitely in the INCAR file in addition

且有一篇讨论帖：https://wiki.vasp.at/forum/viewtopic.php?t=20020

不过我自己尝试的时候发现超算平台上添加NPAR = 1会报错，而在自己电脑上则不会报错，有点迷、、、

用Vesta载入ELFCAR文件

先设置一下等值面数值

Isosurface Level（等值面水平）是一个关键参数，用于控制可视化电子局域化程度的阈值，VESTA会绘制出ELF等于该值的等值面（isosurface），直观展示电子局域化程度高于该阈值的区

Properties-Isosurfaces-Isosurfaces level，设置完后可以按以下Tab键就会自动保存应用

等值面数值设置

这里如果设置0.75的话啥也看不到，因为铜单胞中存在着金属键，并没有很强的局域化电子，而是均匀分布的电子气

高定域性通常出现在共价作用区域、孤对电子区域、原子内核及壳层结构区域

然后就可以查看ELF图了

首先是3D的图

左上角-Edit-Lattice planes，新建一个图，根据想查看的面的来设置米勒指数hlk，以及设置距离d

3D的ELF图

然后是2D的图

左上角-Utilities-2D Data Display-Slice即可

同一个截面的2D图和3D图
001,1xd的二维截面ELF图

可以看到，顶点Cu与顶点Cu原子之间、顶点Cu和面心Cu之间的电子的局域化程度比较高并且没有明显的偏向偏离

接下来看一下一维时的电子定域化情况

左上角-Utilities-Line profile即可

  • 首先棱上的情况
边长上的ELF图

Cu和Cu两原子中间存在均匀的非零的电子局域函数，虽然不足0.5但分布平坦，这是金属键合的特征

两个铜原子的ELF(r)值相近

  • 面对角线上的情况
面对角线

三个铜原子的ELF(r)值相近

离子晶体CsF

再看一下离子键的情况，以CsF为例

CsF的ELF图

氟和铯之间出现深蓝色，这说明在这个区域电子的分布较为离域化，即氟和铯之间不存在共用电子，键的类型偏向离子键

再看看对角线上的电荷分布

对角线上的一维ELF图

两端为铯中间为氟，可以看到氟和铯ELF(r)值差别较大

共价晶体——硅

101截面的ELF图

如箭头所指，硅原子和硅原子之间有较强的电子定域化作用，键的类型为共价键

如下图，体对角线上的ELF图也能很好证明，Si-Si之间的区域电子化定域程度极高，呈现明显的尖峰

对角线上的ELF图

共价晶体——灰锡

我仍记得高中时课本讲的一个稍微有些特殊的例子——灰锡，灰锡是一个共价晶体

灰锡101晶面的ELF图

可以看到，灰锡与晶体硅结构上类似，Sn($\frac{1}{4},\frac{1}{4},\frac{1}{4}$)与Sn(0,0,0)与以及Sn(0.5,0.5,0)成键，Sn($\frac{3}{4},\frac{3}{4},\frac{1}{4}$)与Sn(0.5,0.5,0)和Sn(1,1,0)成键。成键电子局域化程度较高，为共价电子

Comment and share

2025年5月Intel平台VASP.6.4.3的编译安装

in 从零开始的科研之旅, VASP

2025年5月Intel平台VASP.6.4.3的编译安装

WSL的安装

选择Arch Linux，个人比较熟悉，仅供学习使用不作为稳定的生产力环境，配置如下，两颗E5 2676 v3

洋垃圾配置

打开powershell或者终端，输入wsl --install archlinux以自动安装，这一步需要使用代理(MAGIC)，主机上使用Clash-Verge进行代理（运气好也能连上），终端中输入：

1
2
set http_proxy=http://127.0.0.1:7897 #自己使用的端口是7897，http代理
set https_proxy=http://127.0.0.1:7897 #自己使用的端口是7897，https代理

如果没有条件使用代理，也可以手动安装，参考Arch Wiki，https://wiki.archlinuxcn.org/wiki/%E5%9C%A8_WSL_%E4%B8%8A%E5%AE%89%E8%A3%85_Arch_Linux，手动下载.wsl镜像，然后安装

笔者计算机C盘性能太差，于是想要将wsl2迁移到另一个分区

安装成功

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
[root@NAS ~]# fastfetch
                  -`                     root@NAS
                 .o+`                    --------
                `ooo/                    OS: Arch Linux x86_64
               `+oooo:                   Host: Windows Subsystem for Linux - archlinux (2.4.12)
              `+oooooo:                  Kernel: Linux 5.15.167.4-microsoft-standard-WSL2
              -+oooooo+:                 Uptime: 5 hours, 16 mins
            `/:-:++oooo+:                Packages: 155 (pacman)
           `/++++/+++++++:               Shell: bash 5.2.37
          `/++++++++++++++:              WM: WSLg 1.0.65 (Wayland)
         `/+++ooooooooooooo/`            Terminal: xterm-256color
        ooosssso++osssssso+`           CPU: Intel(R) Xeon(R) E5-2676 v3 (48) @ 2.39 GHz
       .oossssso-````/ossssss+`          Memory: 1.09 GiB / 15.57 GiB (7%)
      -osssssso.      :ssssssso.         Swap: 304.95 MiB / 4.00 GiB (7%)
     :osssssss/        osssso+++.        Disk (/): 18.67 GiB / 1006.85 GiB (2%) - ext4
    /ossssssss/        +ssssooo/-        Disk (/mnt/c): 159.05 GiB / 237.01 GiB (67%) - 9p
  `/ossssso+/:-        -:/+osssso+-      Disk (/mnt/d): 13.16 TiB / 14.55 TiB (90%) - 9p
 `+sso+:-`                 `.-/+oso:     Disk (/mnt/e): 1.70 TiB / 3.64 TiB (47%) - 9p
`++:.                           `-/+/    Disk (/mnt/f): 49.45 GiB / 476.94 GiB (10%) - 9p
.`                                 `/    Disk (/mnt/g): 25.16 TiB / 29.10 TiB (86%) - 9p
                                         Local IP (eth0): 172.31.36.99/20
                                         Locale: en_US.UTF-8

编译器的安装

1
2
3
4
5
6
7
8
9
10
11
pacman -S make
pacman -S cmake
pacman -S gcc
pacman -S gcc-fortran
####发现gcc和gfortran编译器太新，会编报错(C23)，故降级软件包
#添加中科大源来安装yay
pacman -Syy yay
yay -S downgrade #安装降级工具
downgrade -S gcc-libs #选择12.1.0版本，安装gcc的依赖
downgrade -S gcc #安装gcc，选择12.1.0版本
downgrade -S gcc-fortran #安装gforran，选择12.1.0版本

前往Intel官网下载oneAPI BaseTool kithttps://www.intel.com/content/www/us/en/developer/tools/oneapi/base-toolkit-download.html?packages=oneapi-toolkit&oneapi-toolkit-os=linux&oneapi-lin=online

1
2
chmod +x xxxxxxxxxx.sh #赋予下载下来的sh文件可执行权限
xxxx.sh #运行sh文件

oneAPI BaseTool kit不需要全部安装，只需安装MKL部分，

HPC也只需安装一部分

此处博主踩了较多的坑：

  • 编译器采用2022版，不可采用过新编译器
    • 2024年10月后，Intel删除了ifort编译器
  • gcc编译器不可太新，否则编译某些fortran模块会报错，博主从14降级到了12.1.0

VASP.6.4.3的编译安装

源码包可以去各大私密性较强的论坛或者微信公众号获取，仅供个人学习使用，不用于学术、商业用途

选择arch/目录中的模板makefile.include.intel，复制到.(上一层目录)，并修改名称为makefile.include，修改目录如下

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
# Default precompiler options
CPP_OPTIONS = -DHOST=\"LinuxIFC\" \
              -DMPI -DMPI_BLOCK=8000 -Duse_collective \
              -DscaLAPACK \
              -DCACHE_SIZE=4000 \
              -Davoidalloc \
              -Dvasp6 \
              -Duse_bse_te \
              -Dtbdyn \
              -Dfock_dblbuf

CPP         = fpp -f_com=no -free -w0  $*$(FUFFIX) $*$(SUFFIX) $(CPP_OPTIONS)

FC          = mpiifort
FCL         = mpiifort

FREE        = -free -names lowercase

FFLAGS      = -assume byterecl -w

OFLAG       = -O2 -xhost  #添加-xhost标记
OFLAG_IN    = $(OFLAG)
DEBUG       = -O0

OBJECTS     = fftmpiw.o fftmpi_map.o fftw3d.o fft3dlib.o
OBJECTS_O1 += fftw3d.o fftmpi.o fftmpiw.o
OBJECTS_O2 += fft3dlib.o

# For what used to be vasp.5.lib
CPP_LIB     = $(CPP)
FC_LIB      = $(FC)
CC_LIB      = icc
CFLAGS_LIB  = -O
FFLAGS_LIB  = -O1
FREE_LIB    = $(FREE)

OBJECTS_LIB = linpack_double.o

# For the parser library
CXX_PARS    = icpc
LLIBS       = -lstdc++

##
## Customize as of this point! Of course you may change the preceding
## part of this file as well if you like, but it should rarely be
## necessary ...
##

# When compiling on the target machine itself, change this to the
# relevant target when cross-compiling for another architecture
VASP_TARGET_CPU ?= -xHOST
FFLAGS     += $(VASP_TARGET_CPU)

# Intel MKL (FFTW, BLAS, LAPACK, and scaLAPACK)
# (Note: for Intel Parallel Studio's MKL use -mkl instead of -qmkl)
FCL        += -qmkl=sequential
MKLROOT    ?= /opt/intel/oneapi/mkl/2022.1.0 #一定要添加MKLROOT的目录，否则会报错


/home/storm/intel/oneapi/mkl/2022.1.0
LLIBS      += -L$(MKLROOT)/lib/intel64 -lmkl_scalapack_lp64 -lmkl_blacs_intelmpi_lp64
INCS        =-I$(MKLROOT)/include/fftw

# HDF5-support (optional but strongly recommended)
#CPP_OPTIONS+= -DVASP_HDF5
#HDF5_ROOT  ?= /path/to/your/hdf5/installation
#LLIBS      += -L$(HDF5_ROOT)/lib -lhdf5_fortran
#INCS       += -I$(HDF5_ROOT)/include

# For the VASP-2-Wannier90 interface (optional)
#CPP_OPTIONS    += -DVASP2WANNIER90
#WANNIER90_ROOT ?= /path/to/your/wannier90/installation
#LLIBS          += -L$(WANNIER90_ROOT)/lib -lwannier

然后开始编译

发现使用命令make all -j48会编译报错

于是改为官方wiki的命令make DEPS=1 -j48，编译成功

测试

make test进行测试，有一些报错，懒得管了、、、能跑就行

Comment and share

Vaspkit的使用

in 从零开始的科研之旅, VASP

安装、配置好vaspkit之后输入vaspkit可启动软件包

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
(base) storm@DESKTOP-HE4FQ8Q:~/my-learn/trash-vaspkit$ vaspkit
            \\\///
           / _  _ \         Hey, you must know what you are doing.
         (| (o)(o) |)       Otherwise you might get wrong results.
 o-----.OOOo--()--oOOO.------------------------------------------o
 |         VASPKIT Standard Edition 1.5.1 (27 Jan. 2024)         |
 |         Lead Developer: Vei WANG (wangvei@icloud.com)         |
 |      Main Contributors: Gang TANG, Nan XU & Jin-Cheng LIU     |
 |  Online Tutorials Available on Website: https://vaspkit.com   |
 o-----.oooO-----------------------------------------------------o
        (   )   Oooo.                          VASPKIT Made Simple
         \ (    (   )
          \_)    ) /
                (_/
 ===================== Structural Utilities ======================
 01) VASP Input-Files Generator    02) Mechanical Properties
 03) K-Path for Band-Structure     04) Structure Editor
 05) Catalysis-ElectroChem Kit     06) Symmetry Analysis
 07) Materials Databases           08) Advanced Structure Models
 ===================== Electronic Utilities ======================
 11) Density-of-States             21) Band-Structure
 23) 3D Band-Structure             25) Hybrid-DFT Band-Structure
 26) Fermi-Surface                 28) Band-Structure Unfolding
 31) Charge-Density Analysis       42) Potential Analysis
 44) Piezoelectric Properties      51) Wave-Function Analysis
 62) Magnetic Analysis             65) Spin-Texture
 68) Transport Properties
 ======================== Misc Utilities =========================
 71) Optical Properties            72) Molecular-Dynamics Kit
 74) User Interface                78) VASP2other Interface
 84) ABACUS Interface              91) Semiconductor Kit
 92) 2D-Material Kit               95) Phonon Analysis
 0)  Quit
 ------------>>

以下是中文翻译对照：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
===================== 结构工具 ======================
01) VASP 输入文件生成器
02) 力学性质
03) 能带结构的 K 路径
04) 结构编辑器
05) 催化 - 电化学套件
06) 对称性分析
07) 材料数据库
08) 高级结构模型
==================== 电子工具 ======================
11) 密度态
21) 能带结构
23) 三维能带结构
25) 混合 DFT 能带结构
26) 费米面
28) 能带结构展开
31) 电荷密度分析
42) 势分析
44) 压电性质
51) 波函数分析
62) 磁性分析
65) 自旋纹理
68) 输运性质
====================== 杂项工具 ========================
71) 光学性质
72) 分子动力学套件
74) 用户界面
78) VASP 转其他格式接口
84) ABACUS 接口
91) 半导体套件
92) 二维材料套件
95) 声子分析
0) 退出

我目前还是初学者小白，使用得最多的功能还是生成输入文件功能尤其是生成INCAR文件，

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
 ------------>>
01
 ==================== VASP Input Files Options ===================
 101) Customize INCAR File
 102) Generate KPOINTS File for SCF Calculation
 103) Generate POTCAR File with Default Setting
 104) Generate POTCAR File with User Specified Potential
 105) Generate POSCAR File from cif (no fractional occupations)
 106) Generate POSCAR File from Material Studio xsd (retain fixes)
 107) Reformat POSCAR File in Specified Order of Elements
 108) Successive Procedure to Generate VASP Files and Check
 109) Submit Job Queue

 0)   Quit
 9)   Back
 ------------>>
  ------------>>
101
 +---------------------------- Tip ------------------------------+
 |          WARNNING: You MUST know what wou are doing!          |
 |Some Parameters in INCAR file need to be set/adjusted manually.|
 +---------------------------------------------------------------+
 ======================== INCAR Options ==========================
 ST) Static-Calculation            SR) Standard Relaxation
 MG) Magnetic Properties           SO) Spin-Orbit Coupling
 D3) DFT-D3 no-damping Correction  H6) HSE06 Calculation
 PU) DFT+U Calculation             MD) Molecular Dynamics
 GW) GW0 Calculation               BS) BSE Calculation
 DC) Elastic Constant              EL) ELF Calculation
 BD) Bader Charge Analysis         OP) Optical Properties
 EC) Static Dielectric Constant    PC) Decomposed Charge Density
 PH) Phonon-Calculation            PY) Phonon with Phononpy
 NE) Nudged Elastic Band (NEB)     DM) The Dimer Method
 FQ) Frequence Calculation         LR) Lattice Relaxation
 MT) Meta-GGA Calculation          PZ) Piezoelectric Calculation

 0)   Quit
 9)   Back
 ------------>>
 Input Key-Parameters (STH6D3 means HSE06-D3 Static-Calcualtion)

以下是中文对照：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
+---------------------------- 提示 ------------------------------+
|          警告：你必须清楚自己在做什么！          |
|INCAR 文件中的一些参数需要手动设置/调整。|
+---------------------------------------------------------------+
======================= INCAR 选项 ==========================
ST) 静态计算            SR) 标准松弛
MG) 磁性性质            SO) 自旋轨道耦合
D3) DFT-D3 无阻尼校正  H6) HSE06 计算
PU) DFT+U 计算          MD) 分子动力学
GW) GW0 计算            BS) BSE 计算
DC) 弹性常数            EL) ELF 计算
BD) Bader 电荷分析      OP) 光学性质
EC) 静态介电常数        PC) 分解电荷密度
PH) 声子计算            PY) 使用 Phonopy 的声子计算
NE) 强制弹性带（NEB）   DM) 二聚体方法
FQ) 频率计算            LR) 晶格松弛
MT) Meta-GGA 计算       PZ) 压电计算

LR晶格弛豫就对应着结构优化，当然生成的INCAR文件可能仍需要进行一些修改，比如ENCUT可能需要调整为赝势中ENMAX的1.3倍

生成KPOINTS文件

对于非能带的计算，选择102) Generate KPOINTS File for SCF Calculation使用程序自动撒点即可，但是需要用户选择撒点方式和撒点密度，选择Gamma Scheme，使用0.04精度

关于生成 KPOINTS 文件

  • 在 VASP 计算中，KPOINTS 文件用于指定 k 点网格的设置。对于非能带计算，比如进行能量计算、优化结构等，通常不需要像计算能带结构那样沿特定的高对称路径采样 k 点，而是可以在整个布里渊区均匀地分布 k 点，这就是所说的 “自动撒点”。

撒点方式

  • Monkhorst-Pack 方式 ：这是最常用的自动撒点方式。它通过指定一组三个整数（nx, ny, nz），在三个倒易晶格方向上均匀地分布 k 点，形成的网格可以系统地逼近布里渊区的积分。例如，对于简单晶格，若指定 k 点网格为 3×3×3，程序就会按照 Monkhorst-Pack 方式在这三个方向各取 3 个点，生成一个均匀分布的 k 点网格。
  • Gamma 方式 ：这种方式会在 k 点网格中包含 Γ 点（布里渊区的中心点）。对于金属等具有高对称性的材料，在 Γ 点附近采样很重要，这时候 Gamma 方式撒点更合适。比如指定 k 点网格为 4×4×4 并采用 Gamma 方式撒点，生成的网格会以 Γ 点为中心，向周围均匀分布其他 k 点。

K 点密度

  • K 点密度指的是在布里渊区单位体积内分布的 k 点数量，它决定了 k 点网格的疏密程度。较高的 K 点密度意味着更精细地采样布里渊区，能够更精确地计算系统的物理性质，但也会增加计算量和资源消耗。
  • 例如，对于一个计算资源充足的大型超算任务，为了获得高精度结果，可以选择较大的 K 点密度，如 8×8×8 甚至更高；而对于初步探索或计算资源有限的情况，可以适当降低 K 点密度，如 4×4×4。

生成KPOINTS的同时，POTCAR也会被自动生成，前提是在vaspkit的配置文件中正确设置了赝势的路径和赝势稳健的种类（默认type为PBE）

Comment and share

离子电导率的计算

in 从零开始的科研之旅, VASP

使用VASP计算离子电导率，似乎比较麻烦

一个能量较低、结构合理的起始三维构型是进行MD模拟的第一步

分子动力学模拟的基本步骤

  • 确定起始构型
  • 选用适当的力场和模拟软件
  • 构建体系和能量最小化
  • 平衡过程
  • 数据采集过程
  • 结果分析

根据系统的传播方式，分子模拟方法可以分为分子动力学MD和蒙特卡洛MC两类方法

  • 分子动力学通过数值积分求解牛顿运动方程，生成系统轨迹，可研究结构、动力学和热力学性质。
  • 蒙特卡洛方法基于概率规则生成新构型，侧重采样平衡态结构和热力学性质，不涉及时间动力学。

力场是一组描述分子中原子间相互作用的势能函数和参数的集合，用于计算系统的总能量和原子受力，从而驱动分子动力学（MD）模拟或蒙特卡罗（MC）采样：

  • 将量子力学描述的分子相互作用简化为经典力学模型，降低计算成本。
  • 通过参数化的势能函数，预测分子的结构、动力学和热力学性质（如键长、自由能、扩散系数等）。

固定电荷力场：固定原子电荷，不考虑环境引起的电子云变形（极化），计算效率高

  • 生物分子力场
    • AMBER
      • 广泛用于蛋白质、核酸模拟，参数优化针对溶液环境，如 ff14SB、ff19SB
    • CHARMM
      • 广泛用于蛋白质、核酸模拟，参数优化针对溶液环境，如 ff14SB、ff19SB
    • GROMOS
      • 用于膜蛋白和聚合物，参数简洁，计算速度快
  • 小分子力场
    • GAFF
      • 可参数化多种有机小分子，支持自定义分子
    • UFF
      • 覆盖周期表多数元素，适合无机材料和有机分子

极化力场：显示考虑原子极化率（电子云随电场变化），通过附加自由度（如 Drude 振子、诱导偶极）模拟极化效应，精度更高但计算成本显著增加

  • AMOEBA：基于原子极化率和多极矩，适用于极性体系（如水、蛋白质 - 配体相互作用）。
  • Drude 模型：每个原子附加虚拟振子模拟极化，如 CHARMM-Drude 力场。
  • Polarizable AMBER (PAMBS)：在 AMBER 框架中引入极化参数。

粗粒化力场（Coarse-Grained Force Fields）：将多个原子简化为单个 “珠子”（如蛋白质的每个残基、脂质的疏水尾），降低计算复杂度，适用于大尺度系统（如细胞膜、聚合物）

  • Martini：广泛用于膜蛋白、病毒颗粒和软物质，4:1 粗粒化（4 个原子→1 个珠子）。
  • MARTINI 3.0：支持多尺度模拟，兼容全原子模型。

使用pymatgen提取数据

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
'''
分析AIMD结果，计算MSD 和 conductivity
'''
import os
from pymatgen.core.trajectory import Trajectory
from pymatgen.io.vasp.outputs import Xdatcar
from pymatgen.core import Structure
from pymatgen.analysis.diffusion.analyzer import DiffusionAnalyzer
import numpy as np
import pickle
import matplotlib.pyplot as plt

# 这一步是读取 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() 方法，得到数据后再自己作图
ax = diff.export_msdt("msd.dat")
#plt.show()

# 接下来直接得到 离子迁移率， 单位是 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
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
import numpy as np
import matplotlib.pyplot as plt

# 读取数据
data = np.loadtxt('msd.txt', comments='#')  # 假设你的文件名为 msd.txt

# 提取时间步长和 MSD 数据
time = data[:, 0]  # 第一列是时间步长
msd = data[:, 1]   # 第二列是总的 MSD

# 绘制 MSD 曲线
plt.figure(figsize=(10, 6))
plt.plot(time, msd, label='MSD', color='blue', linewidth=2)
plt.xlabel('Time (fs)', fontsize=14)
plt.ylabel('Mean Square Displacement (Ų)', fontsize=14)
plt.title('MSD vs Time', fontsize=16)
plt.legend(fontsize=12)
plt.grid(True)
plt.show()
  1. System preparation
  2. Minimization/Relaxation
  3. Equilibration
  4. Production

系统准备可以说是模拟中最关键的阶段，但在许多情况下却得到了最少的关注；如果你的系统准备有缺陷，这些缺陷可能会证明是致命的。最坏的可能结果是，准备好的系统并非你想要的（例如，它包含不正确的分子或质子化状态），但在化学上却是有效的，并且可以被你的力场很好地描述，因此没有错误地通过了剩余的步骤——实际上，这是系统准备问题中常见的后果。如果系统在随后的平衡步骤中行为良好，不应假设系统已正确准备；这里应该格外小心。

生产流程如下：

生产流程图

为什么需要有恒温器？

分子动力学模拟用于观察和获取某个研究系统的有趣性质。在很多情况下，为了模拟实验室条件下进行的实验，希望从正则（恒温）配分函数中进行抽样。一般来说，如果系统在模拟过程中必须保持温度不变，将采用某种温控算法。

分子动力学模拟的温度通常使用根据等分配定理定义的动能来测量： $$ \frac{3}{2}Nk_BT = <\sum_{i=1}^{N}\frac{1}{2}m_iv_i^2> $$ 带角标的括号表示温度是一个时间平均值。如果我们使用等分配定理来计算分子动力学模拟单时间点的温度快照，而不是时间平均，这个量称为瞬时温度。瞬时温度不一定等于目标温度；实际上，在正则配分函数中，瞬时温度应该围绕目标温度波动。

温控算法通过改变本质上微观可逆的牛顿运动方程来工作。因此，如果希望计算扩散系数等动力学性质，则最好不使用温控器；相反，在系统达到所需温度后，应关闭温控器。然而，虽然所有温控器都会给出非物理动力学，但有些研究发现它们对特定动力学性质的计算影响很小，并且它们在生产模拟中也常被使用。

一些恒温器

  • Gaussian
    • 使用了高斯最小约束原理来确定维持瞬时温度所需的最小扰动力。这种温控器不采样正则分布；相反，它采样的是等速（恒动能）配分函数。但是，等速配分函数采样与正则配分函数相同的配置相空间，因此，使用任一配分函数都可以等效地获得位置依赖（结构）平衡性质。然而，速度依赖（动力学）性质在两个配分函数之间不会等效。这种温控器仅适用于某些高级应用。
  • Simple Velocity Rescaling
    • 实施起来最容易的温控器之一,也是最不物理的温控器之一。这种温控器依赖于重新标粒子的动量，使模拟的瞬时温度恰好与目标温度相匹配，不建议使用
  • Berendsen
    • 既不采样正则分布也不采样等速分布，它的使用可能导致模拟伪象，不建议使用
  • Bussi-Donadio-Parrinello (Canonical Sampling Velocity Rescaling)
    • 正确采样了正则配分函数，也被称为“V-rescale”热浴方法，非常常用的热浴方法
  • Andersen
    • 只能用于采样结构性质，因为动力学性质可能会受到突然碰撞的很大影响，不建议使用
  • Langevin
    • 不了解
  • Noés-Hoover
    • 不了解
压浴类型 原理 特点 适用性 局限性
Simple volume rescaling 每次修改系统体积使瞬时压强等于目标压强 操作简单直接 仅可用于模拟初始阶段，帮助系统初步接近目标压强 不能采样正确配分函数，不能用于生产采样；不平滑地接近目标压强，可能导致系统出现不真实的物理问题
Berendsen 将系统耦合到弱相互作用的压力浴，通过缩放因子定期缩放体积 比简单体积缩放更真实地接近目标压强，产生更真实的压力波动 对于平衡开始阶段可能有用 采样配分函数定义不佳，不能保证是 NPT 或 NPH，不应用于生产采样
Andersen 在运动方程中增加额外自由度，将系统耦合到假想压力浴（类似各向同性活塞作用） 采样正确的配分函数 适用于对配分函数采样要求高的场景 是各向同性的，不能对系统的一部分施加各向异性压力
Parrinello - Rahman 对 Andersen 压缩计的扩展，支持模拟盒大小和形状的各向异性缩放 具有 Andersen 压缩计的性质，增加了各向异性支持，适用于固体模拟（可模拟晶体格子中相变形状变化）等 在固体模拟等领域应用优势明显 并非所有模拟器都已实现该压缩计，限制了用户选择
Martyna - Tuckerman - Tobias - Klein (MTTK) 引入替代运动方程，正确采样较小系统配分函数 通常被视为对 Parrinello - Rahman 的改进，适用于小系统 为小系统的模拟提供了更准确的配分函数采样方法 未提及明显独特的局限性（与其他扩展自由度压缩计类似，受模拟器实现情况等限制）

Comment and share

截断能ENCUT的选取

in 从零开始的科研之旅, VASP

ENCUT是用来指定平面波基组能量截止值的参数

以上是wiki对ENCUT的描述，POTCAR文件中已经包含了元素的ENMAX和ENMIN，这也是互联网上对ENCUT的选择方案对通常是1.3倍的ENMAX

侯柱峰老师提供了另外一种方法并指出通过该方法选择的ENCUT通常能满足1.3倍ENMAX，供参考使用

截断能指定了用于波函数展开的平面波基组的截断能量，此能量越大则用来描述波函数的平面波基组越多，精度越高，但计算也越耗时。

可以通过计算测试来选择合适的ENCUT。

以下是用于测试的脚本：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
#!/bin/sh
rm WAVECAR
for i in 150 200 250 300 350 400
do
cat > INCAR <<!
SYSTEM = Si-Diamond
ENCUT = $i
ISTART = 0 ; ICHARG = 2
ISMEAR = -5
PREC = Accurate
!
echo "ENCUT = $i eV" ; time vasp
E=`grep "TOTEN" OUTCAR | tail -1 | awk '{printf "%12.6f \n", $5 }'`
echo $i $E >>comment
done

计算后会得到comment文件，它列出了在不同ENCUT下计算得到的总能量

1
2
3
4
5
6
150 -11.900655
200 -11.938864
250 -11.944599
300 -11.945248
350 -11.945503
400 -11.945622

==总能量变化在0.001 eV即可==，因此在这个例子中ENCUT可以选择250

Comment and share

表面能的计算

in 从零开始的科研之旅, VASP

关于表面能

根据VASP官方手册Hands on Session III，表面能的定义如下，大师兄的文章也提到了这个定义( Nanoscale, 2017, 9, 13089–13094) $$ \sigma=\sigma^{unrelax}+E^{relax}=\dfrac{1}{2}(E_{surf}-N_{atoms}\cdot E_{bulk})+E^{relax} $$

$$ \gamma_s=\dfrac{1}{2A}(E_{surf}^{unrelax}-N\cdot E_{bulk})+\dfrac{1}{A}(E_{surf}^{relax}-E_{surf}^{unrelax}) $$

其中，弛豫能量被定义为： Erelax = Esurfrelax − Esurfunrelax

  • 手册的公式和大师兄文章中公式的区别在于表面积A，表面能应当除以表面积，所以后者更合理一些
  • Ebulk为体相中单个原子的能量，所以N ⋅ Ebulk为这个体相的能量，即为OUTCAR直接读出的能量

关于表面积A的计算

这部分有些难理解，涉及到一点高等数学和线性代数🤕

表面积可通过晶格基矢量的叉积确定，在POSCAR或CONTCAR文件中，晶格基矢量通常以三行的形式定义

对于面积的计算，可设一组不共线的基矢量为a1 = (a11, a12, 0), a2 = (a21, a22, 0)

下面这张图应该挺形象

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

a⃗ × b⃗ = (0, 0, a11 × a22 − a12 × a21)

面积为外积模乘以缩放银子 实际面积 = s2 ⋅ |a11a22 − a12a21| PS：对于三维的体积 体积 = s3 ⋅ |a1 ⋅ (a2 × a3)|

p(1x1)Cu(111)的表面能

前面已经算得差不多了

slab模型的能量信息

1
2
3
4
5
6
[ctan@baifq-hpc141 primitive-slab-opt-1]$ grep '  without' OUTCAR
  energy  without entropy=      -13.96941183  energy(sigma->0) =      -13.97673168
  energy  without entropy=      -13.96945528  energy(sigma->0) =      -13.97676958
  energy  without entropy=      -13.97006413  energy(sigma->0) =      -13.97726241
  energy  without entropy=      -13.97020643  energy(sigma->0) =      -13.97735554
  energy  without entropy=      -13.97017969  energy(sigma->0) =      -13.97735013

slab模型优化后的结构信息

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
CONTCAR\(1\1\1)
   1.00000000000000
     2.5701000690000000    0.0000000000000000    0.0000000000000000
    -1.2850500345000000    2.2257719501000000    0.0000000000000000
     0.0000000000000000    0.0000000000000000   21.2954006195000005
  Cu/0e71e558f37
               4
Selective dynamics
Direct
  0.0000000000000000  0.0000000000000000  0.0000000000000000   F   F   F
  0.6666700000000034  0.3333299999999966  0.0985399999999998   F   F   F
  0.3333311322281505  0.6666688677718494  0.1972481653092927   T   T   T
 -0.0000004145140272  0.0000004145140272  0.2951613337135311   T   T   T

  0.00000000E+00  0.00000000E+00  0.00000000E+00
  0.00000000E+00  0.00000000E+00  0.00000000E+00
  0.00000000E+00  0.00000000E+00  0.00000000E+00
  0.00000000E+00  0.00000000E+00  0.00000000E+00

bulk模型的能量信息

1
2
3
4
[ctan@baifq-hpc141 primitive-slab-opt-1]$ grep '  without' .bulk-calc/OUTCAR
  energy  without entropy=      -14.90258332  energy(sigma->0) =      -14.90537746
  energy  without entropy=      -14.90568039  energy(sigma->0) =      -14.90844732
  energy  without entropy=      -14.90859467  energy(sigma->0) =      -14.91129762
  • A = 2.570100069 * 2.2257719501 − 0.0000000000000000 * (−1.2850500345) = 5.720456642530275 Å2
  • Esurfunrelax = −13.9767 eV
  • Ebulk ⋅ N = −14.9112 eV
  • Esurfacerelax = −13.9773 eV

所以表面能γs = 0.08158 eV/Å2 ,又1 eV/Å2=16.0218 J/m2，所以表面能γs= 1.3070 J/m2

c(1x1)Cu(111)的表面能

之前大师兄提到过，对于面心立方最密堆积的形式，slab模型的建立可以使用原胞也可以使用单胞，那么接下来看看用单胞计算表面能吧

c(1x1)Cu(111)优化计算的能量信息：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
[ctan@baifq-hpc141 conventional-slab-opt]$ grep '  without' OUTCAR
  energy  without entropy=      -55.94972046  energy(sigma->0) =      -55.94434499
  energy  without entropy=      -55.94997517  energy(sigma->0) =      -55.94458201
  energy  without entropy=      -55.95396000  energy(sigma->0) =      -55.94820485
  energy  without entropy=      -55.95519341  energy(sigma->0) =      -55.94928561
  energy  without entropy=      -55.95554732  energy(sigma->0) =      -55.94959634
  energy  without entropy=      -55.95713778  energy(sigma->0) =      -55.95095507
  energy  without entropy=      -55.95788514  energy(sigma->0) =      -55.95158140
  energy  without entropy=      -55.95893404  energy(sigma->0) =      -55.95242747
  energy  without entropy=      -55.95963672  energy(sigma->0) =      -55.95296199
  energy  without entropy=      -55.96031694  energy(sigma->0) =      -55.95340564
  energy  without entropy=      -55.96050719  energy(sigma->0) =      -55.95349958
  energy  without entropy=      -55.96053048  energy(sigma->0) =      -55.95316488
  energy  without entropy=      -55.96067995  energy(sigma->0) =      -55.95346716
  energy  without entropy=      -55.96069855  energy(sigma->0) =      -55.95356576
  energy  without entropy=      -55.96051601  energy(sigma->0) =      -55.95360396
  energy  without entropy=      -55.96050802  energy(sigma->0) =      -55.95361124
  energy  without entropy=      -55.96041586  energy(sigma->0) =      -55.95360050
  energy  without entropy=      -55.96018279  energy(sigma->0) =      -55.95347398

结构信息：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
[ctan@baifq-hpc141 conventional-slab-opt]$ cat CONTCAR
CONTCAR\(1\1\1)
   1.00000000000000
     5.1402001381000000    0.0000000000000000    0.0000000000000000
    -2.5701000690000000    4.4515439000999999    0.0000000000000000
     0.0000000000000000    0.0000000000000000   21.2954006195000005
  Cu/0e71e558f37
              16
Direct
  0.0000006270807907 -0.0000006270807907  0.0011656342599333
  0.6666672048402060  0.3333327951597896  0.0988195745854847
  0.3333327951597896  0.6666672048402060  0.1968004254145135
 -0.0000006270807907  0.0000006270807907  0.2944543657400676
  0.0000006270807907  0.4999993729192098  0.0011656342599333
  0.6666672048402060  0.8333327951597940  0.0988195745854847
  0.3333327951597896  0.1666672048402087  0.1968004254145135
 -0.0000006270807907  0.5000006270807924  0.2944543657400676
  0.5000006270807924  0.4999993729192098  0.0011656342599333
  0.1666672048402087  0.8333327951597940  0.0988195745854847
  0.8333327951597940  0.1666672048402087  0.1968004254145135
  0.4999993729192098  0.5000006270807924  0.2944543657400676
  0.5000006270807924 -0.0000006270807907  0.0011656342599333
  0.1666672048402087  0.3333327951597896  0.0988195745854847
  0.8333327951597940  0.6666672048402060  0.1968004254145135
  0.4999993729192098  0.0000006270807907  0.2944543657400676

  0.00000000E+00  0.00000000E+00  0.00000000E+00
  0.00000000E+00  0.00000000E+00  0.00000000E+00
  0.00000000E+00  0.00000000E+00  0.00000000E+00
  0.00000000E+00  0.00000000E+00  0.00000000E+00
  0.00000000E+00  0.00000000E+00  0.00000000E+00
  0.00000000E+00  0.00000000E+00  0.00000000E+00
  0.00000000E+00  0.00000000E+00  0.00000000E+00
  0.00000000E+00  0.00000000E+00  0.00000000E+00
  0.00000000E+00  0.00000000E+00  0.00000000E+00
  0.00000000E+00  0.00000000E+00  0.00000000E+00
  0.00000000E+00  0.00000000E+00  0.00000000E+00
  0.00000000E+00  0.00000000E+00  0.00000000E+00
  0.00000000E+00  0.00000000E+00  0.00000000E+00
  0.00000000E+00  0.00000000E+00  0.00000000E+00
  0.00000000E+00  0.00000000E+00  0.00000000E+00
  0.00000000E+00  0.00000000E+00  0.00000000E+00

Cu(111)的体相能量：

1
2
3
4
[ctan@baifq-hpc141 conventional-slab-opt]$ grep '  without' .bulk-calc/OUTCAR
  energy  without entropy=      -14.90258332  energy(sigma->0) =      -14.90537746
  energy  without entropy=      -14.90568039  energy(sigma->0) =      -14.90844732
  energy  without entropy=      -14.90859467  energy(sigma->0) =      -14.91129762
  • 表面积A=5.1402001381000000 * 4.4515439000999999 = 22.8818 Å2
  • Esurfunrelax = -55.9443 eV
  • N ⋅ Ebulk = -59.6452 eV
    • N指的是slab模型中的原子个数而非bulk中的原子个数
  • Esurfrelax = -55.9535 eV

计算得表面能γs=0.08047 eV/Å2 = 1.2892 J/m2

这与p(1x1)Cu(111)计算所得的表面能有一定的差距，目前现阶段我还无法解决这个问题🤔😵‍💫

p(1x1)Ni(100)的表面能

这个例子来自于vasp官方的手册Hands on Session III

1
2
3
4
5
6
7
8
9
10
11
12
general:
SYSTEM = clean Ni(100) surface
ISTART = 0 ; ICHARG=2
ENCUT = 270
ISMEAR = 2 ; SIGMA = 0.2
spin:
ISPIN=2
MAGMOM = 5*1
dynamic:
IBRION = 1
NSW = 100
POTIM = 0.2

有个疑惑的点是手册中的INCAR没有加上偶极矫正，我在自己算的时候加上了偶极矫正

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
System = Nickel-slab-opt
ISMEAR = 1
ISTART = 0
ICHARG = 2
SIGMA = 0.05
ENCUT = 360
ALGO = Normal
ISPIN = 2
MAGMOM = 5*1
NSW = 200
NELM = 200
EDIFF = 1E-6
EDIFFG = -0.01
IDIPOL = 3
LDIPOL = .TRUE.
LWAVE = .FALSE.
LCHARG = .FALSE.
IBRION = 1
POTIM = 0.05

体相能量信息：

1
2
3
4
[ctan@baifq-hpc141 primitive-opt]$ grep '  without' .bulk/OUTCAR
  energy  without entropy=      -21.88865316  energy(sigma->0) =      -21.88759144
  energy  without entropy=      -21.88930785  energy(sigma->0) =      -21.88825056
  energy  without entropy=      -21.89056178  energy(sigma->0) =      -21.88952707

弛豫与未弛豫的slab模型的能量信息：

1
2
3
4
5
6
7
8
[ctan@baifq-hpc141 primitive-opt]$ grep '  without' OUTCAR
  energy  without entropy=      -25.97155386  energy(sigma->0) =      -25.97291371
  energy  without entropy=      -25.97233802  energy(sigma->0) =      -25.97366730
  energy  without entropy=      -25.97412653  energy(sigma->0) =      -25.97534340
  energy  without entropy=      -25.97511519  energy(sigma->0) =      -25.97602399
  energy  without entropy=      -25.97516432  energy(sigma->0) =      -25.97619350
  energy  without entropy=      -25.97516747  energy(sigma->0) =      -25.97619463
  energy  without entropy=      -25.97517140  energy(sigma->0) =      -25.97620156

优化后的结构：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
CONTCAR\(1\0\0)
   1.00000000000000
     2.4846000671000001    0.0000000000000000    0.0000000000000000
    -1.2423000336000001    2.1517267763999999    0.0000000000000000
     0.0000000000000000    0.0000000000000000   18.7047996521000002
  Ni_pv/368cd815
               5
Direct
  0.6666685456044492  0.3333314543955508  0.0005454007373249
  0.0000003831810720 -0.0000003831810720  0.1078816295897802
  0.3333305745814666  0.6666694254185335  0.2169200598414882
  0.6666681017676376  0.3333318982323624  0.3259932620996842
  0.0000023948653782 -0.0000023948653782  0.4332496477317282

  0.00000000E+00  0.00000000E+00  0.00000000E+00
  0.00000000E+00  0.00000000E+00  0.00000000E+00
  0.00000000E+00  0.00000000E+00  0.00000000E+00
  0.00000000E+00  0.00000000E+00  0.00000000E+00
  0.00000000E+00  0.00000000E+00  0.00000000E+00
  • 表面积A=2.4846000671000001*2.1517267763999999=5.346180493024307 eV
  • Esurfunrelax = -25.97291371 eV
  • Esurfrelax = -25.97620156 eV
  • N ⋅ Ebulk = -21.88952707 eV

表面能 γs = 0.1293 eV/Å2 = 2.0715 J/m2

小结

  • 本节主要是联系表面能的计算
  • Ni(100)表面能的计算选用的赝势为Ni_sv赝势

Comment and share

John Doe

author.bio


author.job


Changchun, China