Lightchaser

三人行，必有我师

发表于 更新于 分类于

截断能

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

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

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

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

{.line-numbers}
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
#!/bin/bash

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

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

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

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

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

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

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

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

K点测试

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

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

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

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

cp INCAR_static INCAR

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


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

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

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

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

显示收敛情况的脚本

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

  • 简单版：

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

  • 小脚本：
    {.line-numbers}
    1
    2
    3
    4
    5
    6
    #!/usr/bin/env bash
    for i in *; do
    if [ -e $i/OUTCAR ]; then
    echo -e  $i "\t" $(grep '  without' $i/OUTCAR |tail -n 1 | awk '{print $7}')
    fi
    done
  • 复杂：
    ediff.sh如下：（安装gnuplot）
    {.line-numbers}
    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
    #!/bin/sh
    #lipai@mail.ustc.edu.cn

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

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

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

    gnuplot <<EOF 
    set term dumb
    set title 'Energy of each ionic step'
    set xlabel 'Ionic steps'
    set ylabel 'Energy(eV)'
    plot 'temp.e' u 1:5 w l  t "Energy in eV"
    set title 'Max Force of each ionic step'
    set xlabel 'Ionic steps'
    set ylabel 'Force (eV/Angstrom)'
    plot 'force.conv' w l t "Force in eV/Angstrom"
    EOF

    if [ $steps -gt 8 ] ; then
    tail -5 force.conv >temp.fff
    tail -5 temp.e >temp.ee
    gnuplot <<EOF 
    set term dumb
    set title 'Energy of each ionic step for the last few steps'
    set xlabel 'Ionic steps'
    set ylabel 'Energy(eV)'
    plot 'temp.ee' u 1:5 w l  t "Energy in eV"
    set title 'Max Force of each ionic step for the last few steps'
    set xlabel 'Ionic steps'
    set ylabel 'Force (eV/Angstrom)'
    plot 'temp.fff' w l  t "Force in eV/Angstrom"
    EOF
    rm temp.fff temp.ee
    fi

    rm temp.e temp.f temp.ff temp.fix temp.fixx

批量提交任务

适用于测试、有限位移法求力常数矩阵等情况。
简单： for i in *; do cd $i ; qsub sub4; cd $OLDPWD; done

{.line-numbers}
1
2
3
4
5
6
7
8
9
10
11
12
#!/bin/sh
root_path='pwd'

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

批量产生文件夹

i in {2..9}; do cp 0.01 0.0$i -r ; done ``` #批量复制0.01文件夹复制成0.02-0.09
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
```sed '3s/0.01/0.02/g' INCAR > INCAR_new```  #将INCAR中的第三行的0.01全部替换成0.02并输出位INCAR_new文件。
```sed -i '3s/0.01/0.02/g' INCAR``` # -i表示直接对源文件进行编辑，也就是说编辑之后源文件被新文件替换掉。
```for i in *; do sed -i "3s/0.05/$i/g" $i/INCAR ; done``` #将文件夹中所有的0.01替换成了文件夹序号。

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

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

节约机时

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

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

消虚频

结构优化消虚频

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

{.line-numbers}
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
# -*- 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}
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
#!/usr/bin/env python
# -*- coding: utf-8 -*-
"""
Qiang Li
Command: python3 vib_correct.py
This is a temporary script file.
"""
import numpy as np
import os

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

# 获取虚频对应的OUTCAR部分    
num_atoms = sum([int(i) for i in lines_pos[6].rstrip().split()])
vib_lines = lines[l_position:l_position + num_atoms] #振动部分 7222到7248行

model_positions = []
vib_dis = []
for line in vib_lines:
    position = [float(i) for i in line.rstrip().split()[:3]]
    vib_infor = [float(i) for i in line.rstrip().split()[3:]] # dx， dy， dz对应的那三列
    model_positions.append(position)
    vib_dis.append(vib_infor)

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

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

滑移堆叠能计算

vaspkit 926

脚本

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

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

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

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

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

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

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

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

画图：

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
import numpy as np
import matplotlib.pyplot as plt

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

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

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

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

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

    if fm is None or afm is None:
        return

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

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

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

    # 计算能量差
    Zd = Zf - Za

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

if __name__ == "__main__":
    main()

常用命令行

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

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

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

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

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

CINEB计算重新计算

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

pdos轨道加和

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

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

grep LOOP+ OUTCAR |wc -l

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

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

vaspkit命令行

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

批量修改INCAR

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

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

计算步骤

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

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

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

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

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

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

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

输入命令 nebmovie.pl 0

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

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

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

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

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

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

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

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

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

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

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

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

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

发表于 更新于 分类于

加压

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

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

其他不变

加电场

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

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

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

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

状态方程拟合

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

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

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

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

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

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

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

  1. 获得一定缩放系数的结构：for i in 0.95 0.96 0.97 0.98 0.99 1.01 1.02 1.03 1.04 1.05; do cp 1.00 $i ; sed -i "2s/1.0/$i/g" $i/POSCAR ; done

  2. 分别进行结构优化:for i in *; do cd $i; 提交任务命令; cd $OLDPWD; done

  3. 提取能量-体积以及能量-压力、压力-体积的数据:for i in *; do echo -e $i "\t" $(grep ' without' $i/OUTCAR | tail -n 1| awk '{print $7}'); done > data

  4. 拟合能量体积曲线，获得平衡晶胞参数

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

分数占据

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

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

disorder

安装

1
2
3
4
5
6
7
8
$ unzip disorder-main
$ cd disorder-main
$ make disorder (Only disorder is compilied)
$ make supercell (Only supercell is compilied)
$ make or make all (Both disorder and supercell are compilied)
$ make clean (rm *.o *.mod disorder supercell)
export PATH=$PATH:disorder_installation_path/disorder-main/bin
(Add this to the ~/.bashrc file as you wish)

Note: The default compiler is ifort, and you can modify the Makefile
in the src directory to use another compiler (e.g. gfortran).

需要准备两个文件SPOSCAR和INDSOD,以examples的1和2为例：

例子1：SnxPb1-xTe

  1. 扩建超胞，保证整除
    1
    2
    3
    4
    5
    6
    7
    8
    9
    Te Pb
       1.0
        13.1323995589999996    0.0000000000000000    0.0000000000000000
         0.0000000000000000   13.1323995589999996    0.0000000000000000
         0.0000000000000000    0.0000000000000000   13.1323995589999996
      Te   Pb
      32   32
    Direct
    * * *
  2. 准备INDSOD文件
    1
    2
    3
    4
    5
    6
    7
    8
    9
    10
    11
    12
    &input
      nsub = 2
      subs = 16,16
      symb = 'Sn','Pb' ! The quotes is unnecessary for the ifort compiler
      site = 2
      prec = 1D-5
    ! lpro = .true.
    ! lpos = .true.
    ! leqa = .true.
    ! lspg = .true.
      lcfg = .false.
    /
    官方解释：
    1
    2
    3
    4
    5
    6
    7
    nsub-tag:nsub = 2 ~ 5; Default: nsub = 2.该位置取代原子种类。本例中是2种。
    subs-tag：subs = integer, integer, …; Default: None。对 应的原子数目。总和加起来应该等于该位点在SPOSCAR中的数目。
    symb-tag：symb = character, character, …; Default: None。元素化学符号。空位用Kw,即汉语kong wei。本例是Sn Pb单引号。
    site-tag：site = integer; Default: 1。所取代的位置在SPOSCAR中的位置。本例是第2.
    prec-tag：prec = real; Default: 1D-5。当软件产生空间群与其他软件识别不一致时调整该参数，一般不用修改。
    lpro-tag：lpro = logical; Default: .false.是否消除重复参数。对于大体系可以打开。
    lpos、leqa、lspg、lcfg：分别控制输出文件POSCAR-X、EQAMAT、SPGMAT、CFGMAT。
  3. 运行disorder。输入disorder即可。
  4. 输出文件：

ATAT软件

简介

  • 什么是ATAT？

其全称为 Alloy Theoretic Automated Toolkit (ATAT)
引用官网的一句话，ATAT is a generic name that refers to a collection of alloy theory tools developped by Axel van de Walle, in collaboration with various research groups and with various sources of financial support.

  • 如何安装ATAT：

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

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

配置~/.ezvasp.rc

1
2
3
#!/bin/csh
set VASPCMD="/path/vasp_std# 设置vasp可执行程序
set POTLDA/PBE="/path/to/pot_lad" # 设置为赝势实际的路径
  • ATAT能做啥？

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

MAPs

软件运行

一）运行MAPS：

  1. 准备lat.in文件

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

3）touch ready

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

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

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

二）运行MAPS与VASP交互：

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

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

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

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

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

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

最终正确版本运行步骤：

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

lat.in文件

vasp.wrap文件

sub.sh文件

2）提交脚本文件即可

说明：

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

什么时候停止计算：

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

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

停止计算：

1
2
$ touch stop
$ touch stoppoll

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

1
2
3
rm maps_is_running pollmach_is_running ready
maps -d &
pollmach runstruct_vasp

常用命令

  • corrdump -2= -3= -4= 产生团簇

  • getclus 查看团簇

    • 第一列为团簇的阶数，第二列为团簇直径，第三列为团簇的多重度
    • clusters.out中有团簇格点位置的详细信息

  • clusterexpand -e -cv energy 得到团簇展开系数

    • 此时产生的ECI文件为energy.eci，而maps产生的ECI为eci.out，后者是形成能的团簇展开结果，前者是直接对vasp能量(即energy文件)的团簇展开
    • 如果需要对其它性质做团簇展开，例如带隙，则可以在每个结构的目录下创建一个bandgap文件,其内容即为带隙的大小，然后运行clusterexpand -e -cv -pa bandgap以得到带隙团簇展开的系数bandgap.eci，其中-pa表示所展开的量已经是一个平均的量

检查计算结果

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

结果可视化

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

输出文件

  1. maps.log：开始可能的warnings:
    ’Not enough known energies to fit CE’
    ’True ground states not = fitted ground states’
    ’New ground states predicted, see predstr.out’
    随着结构增多:
    ’Among structures of known energy, true and predicted ground states agree.’
    ’No other ground states of xx atoms/unit cell or less exist.’
    最后一行有crossvalidation score (in eV/atom)，该值要小于0.025。精确小于0.001.

  2. fit.out:拟合的结构，一个结构一行，共6列，分别代表concentration、 energy 、fitted_energy (energy-fitted_energy)、 weight 、index

  3. predstr.out：预测的能量，一个结构一行，共5列，分别代表concentration、 energy 、predicted_energy、index、status.其中status包括status is either
    b for busy (being calculated),
    e for error,
    u for unknown (not yet calculated) or
    g if that structure is predicted to be a ground state (g can be combined with the above).

  4. gs.out：基态能量，一个结构一行，四列，分别是concentration、 energy、 fitted_energy、 index

  5. gs_str.out：基态的结构，格式和str.out一样。每一个结构后面会有end和一个空行隔开。

  6. eci.out：eci值对应的cluster在clusters.out文件。

  7. clusters.out：对于每一个团簇，第一行是多重度 multiplicity，第二行是团簇维度cluster diameter，第三行是团簇的格点数，下一行是格点对应的坐标。最后空行隔开每一个团簇。

  8. str.out：在每个数字文件夹下面，和晶格输入文件相似，不同是坐标写成3*3矩阵，一个位置只写一个原子。

  9. ref_energy.out：用于计算形成能的参考能量。

后处理

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

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

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

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

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

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

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

热力学计算

必要的输入文件：

  • lat.in
  • clusters.out
  • eci.out
  • gs_str.out
emc2使用示例

蒙特卡洛模拟

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

相图计算

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

例子：

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

构建SQS/SQoS

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

必要文件：

  • rndstr.in
  • clusters.out

rndstr.in的一个例子：

1
2
3
4
5
6
7
8
9
10
11
4.1000 0.000000 0.000000
0.000000 4.1000 0.000000
0.000000 0.000000 4.1000
1.000000 0.000000 0.000000
0.000000 1.000000 0.000000
0.000000 0.000000 1.000000
0.000000 0.500000 0.500000 O=0.6667,N=0.3333
0.500000 0.000000 0.500000 O=0.6667,N=0.3333
0.500000 0.500000 0.000000 O=0.6667,N=0.3333
0.000000 0.000000 0.000000 Ba
0.500000 0.500000 0.500000 Ta

产生团簇文件clusters.out：

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

说明：

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

运行mcsqs -n= &
说明：

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

MCSQS

一、原理

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

二、运行步骤

step1. 准备rndstr.in文件

1 1 1 90 90 90

1 0 0

0 1 0

0 0 1

.0 .0 .0 Cr

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

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

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

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

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

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

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

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

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

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

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

层间结合能、结合力

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

Graphite interlayer distance

vasp上的教程：

准备graphite结构，

INCAR：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
IVDW = 20           
LVDW_EWALD =.TRUE. 
NSW = 1 
IBRION = 2
ISIF = 4
PREC = Accurate
EDIFFG = 1e-5
LWAVE = .FALSE.
LCHARG = .FALSE.
ISMEAR = -5
SIGMA = 0.01
EDIFF = 1e-6
ALGO = Fast
NPAR = 2

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

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
#
# To run VASP this script calls $vasp_std
# (or posibly $vasp_gam and/or $vasp_ncl).
# These variables can be defined by sourcing vaspcmd
. vaspcmd 2> /dev/null

#
# When vaspcmd is not available and $vasp_std,
# $vasp_gam, and/or $vasp_ncl are not set as environment
# variables, you can specify them here
[ -z "`echo $vasp_std`" ] && vasp_std="mpirun -np 8 /path-to-your-vasp/vasp_std"
[ -z "`echo $vasp_gam`" ] && vasp_gam="mpirun -np 8 /path-to-your-vasp/vasp_gam"
[ -z "`echo $vasp_ncl`" ] && vasp_ncl="mpirun -np 8 /path-to-your-vasp/vasp_ncl"

#
# The real work starts here
#

rm results.dat

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

The optimal length of the lattice vector c normal to the stacking direction is determined in a series of single point calculations with varied value of c (all other degrees of freedom are fixed at their experimental values).

The computed c vs. energy dependence is written in the file results.dat and can be visualized e.g. using xmgrace. The optimal value can be obtained using the attached utility (python with numpy or Numeric is needed):

./utilities/fit.py results.dat
This should yield:

1
2
3
4
5
6
7
200 iterations performed
Ch-square: 4.30305519481e-09
---------

       E0(eV):         -37.433456779
       d0(A):  6.65603352689
The computed value of 6.66 Å agrees well with experiment (6.71 Å).

Graphite MBD binding energy

1
2
3
4
5
6
7
8
9
10
11
12
13
14
IVDW = 202           
LVDWEXPANSION =.TRUE. 
NSW = 1 
IBRION = 2
ISIF = 4
PREC = Accurate
EDIFFG = 1e-5
LWAVE = .FALSE.
LCHARG = .FALSE.
ISMEAR = -5
SIGMA = 0.01
EDIFF = 1e-6
ALGO = Fast
NPAR = 2

KPOINTS

Graphite:

1
2
3
4
5
Monkhorst Pack
0
gamma
16 16 8
0 0 0

Graphene:

1
2
3
4
5
Monkhorst Pack
0
gamma
16 16 1
0 0 0

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

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
#
# To run VASP this script calls $vasp_std
# (or posibly $vasp_gam and/or $vasp_ncl).
# These variables can be defined by sourcing vaspcmd
. vaspcmd 2> /dev/null

#
# When vaspcmd is not available and $vasp_std,
# $vasp_gam, and/or $vasp_ncl are not set as environment
# variables, you can specify them here
[ -z "`echo $vasp_std`" ] && vasp_std="mpirun -np 8 /path-to-your-vasp/vasp_std"
[ -z "`echo $vasp_gam`" ] && vasp_gam="mpirun -np 8 /path-to-your-vasp/vasp_gam"
[ -z "`echo $vasp_ncl`" ] && vasp_ncl="mpirun -np 8 /path-to-your-vasp/vasp_ncl"

#
# The real work starts here
#

# Here the work starts
rm results.dat

drct=$(pwd)

for i in graphene graphite
do
  cd $drct/$i
  $vasp_std
done

cd $drct

# obtain total energy for graphite 
en2=$(grep "free  ene" graphite/OUTCAR |tail -1|awk '{print $5}') 

# obtain total energy for graphene
en1=$(grep "free  ene" graphene/OUTCAR |tail -1|awk '{print $5}')

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

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

Note that the calculation is performed in two steps (two separate single-point calculations) in which the energy for bulk graphite and for graphene are obtained. The binding energy is computed automatically and it is written in the file results.dat. (N.B.: for the latter python needs to be available.)

The computed value of 0.050 eV/A is now fairly close to the RPA reference of 0.048 eV/atom [1].

Graphite TS binding energy

1
2
3
4
5
6
7
8
9
10
11
12
13
14
IVDW = 20            
LVDW_EWALD =.TRUE. 
NSW = 1 
IBRION = 2
ISIF = 4
PREC = Accurate
EDIFFG = 1e-5
LWAVE = .FALSE.
LCHARG = .FALSE.
ISMEAR = -5
SIGMA = 0.01
EDIFF = 1e-6
ALGO = Fast
NPAR = 2

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

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
#
# To run VASP this script calls $vasp_std
# (or posibly $vasp_gam and/or $vasp_ncl).
# These variables can be defined by sourcing vaspcmd
. vaspcmd 2> /dev/null

#
# When vaspcmd is not available and $vasp_std,
# $vasp_gam, and/or $vasp_ncl are not set as environment
# variables, you can specify them here
[ -z "`echo $vasp_std`" ] && vasp_std="mpirun -np 8 /path-to-your-vasp/vasp_std"
[ -z "`echo $vasp_gam`" ] && vasp_gam="mpirun -np 8 /path-to-your-vasp/vasp_gam"
[ -z "`echo $vasp_ncl`" ] && vasp_ncl="mpirun -np 8 /path-to-your-vasp/vasp_ncl"

#
# The real work starts here
#

rm results.dat

drct=$(pwd)

for i in graphene graphite
do
  cd $drct/$i
  $vasp_std
done

cd $drct

# obtain total energy for graphite 
en2=$(grep "free  ene" graphite/OUTCAR |tail -1|awk '{print $5}') 

# obtain total energy for graphene
en1=$(grep "free  ene" graphene/OUTCAR |tail -1|awk '{print $5}')

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

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

Note that the calculation is performed in two steps (two separate single-point calculations) in which the energy for bulk graphite and for graphene are obtained. The binding energy is computed automatically and it is written in the file results.dat. (N.B.: for the latter python needs to be available.)

Even though the TS method predicts a reasonable geometry (see the Graphite interlayer distance example) it overestimates the energetics strongly: the computed binding energy of -0.083 eV/atom is too large compared to the RPA reference of 0.048 eV/atom. This overestimation is - at least in part - due to neglecting the many-body interactions (see example Graphite MBD binding energy).

发表于 更新于 分类于

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

VASP输出的能量是什么

在优化完的文件夹下输入

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

这里的F和下面free energy含义相同，E0代表体系内能。
另外在OUTCAR中还可以找到这些
{.line-numbers}
1
2
3
free  energy   TOTEN  =       -44.82987077 eV
energy  without entropy=      -44.82987077  energy(sigma->0) =      -44.82987077
atomic energy  EATOM  =      1542.31592482

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

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

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

形成能

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

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

单点能SP

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

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

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

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

自洽计算

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

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

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

气体优化

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

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

IVDW参数

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

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

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

DFT+U

  1. 为什么加U?
    对于含有d、f轨道电子的强关联体系(Hubbard模型)，电子间存在强烈的在位库仑相互作用，而交换相关泛函中的局域密度近似（LDA）或广义梯度近似（GGA）对电子之间的强在位库仑相互作用描述不准， 此时可将研究体系的交换相关泛函分成两部分计算，也就是用DFT+U的方法进行求解，公式如下：
    其中一部分电子用DFT算法（如LDA，GGA）等可以比较准确地描述；另外的d、f轨道电子通过引入Hubbard项得到正确的描述（U值就是考虑了同一个原子自旋相反的局域电子之间的库仑排斥）
  2. 如何加U?
    DFT计算时输入文件的U值由以下参数确定：
    {.line-numbers}
    1
    2
    3
    4
    5
    LDAU= .TRUE.|.FALSE.：开启/关闭+U功能，默认值为.FALSE.；
    LDAUTYPE=1|2|4指的是+U的类型，默认值是2；其中1为Liechtenstein等提出的旋转不变LSDA+U方法；2为 Dudarev等提出的简化 LSDA+U方法；4与1类似, 但不考虑LSDA交换劈裂；
    LDAUL=-1|1|2|3分别对应不加U、p、d、f轨道加U；
    LDAUU、LDAUJ分别指定电子库伦相互作用项和交换相互作用项（U和J值）；
    LMAXMIX =2/4/6：默认为2，加U计算时该值需大于轨道量子数，对于含有d轨道或f轨道电子的体系需对应增加至4或6。
    注意：
    1）必须为每一类原子指定一个U值，其顺序与结构文件POSCAR中元素排序一致；
    2）同一种元素在不同的体系中U值不同；
    3）U在所有电子中都存在，不管自旋是否相同，J只存在于自旋相同的电子上；当LDAUTYPE =2时，U-J的差值才有意义，即有效的U参数。
    该参数需要测试或者参考文献。

磁性体系

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

{.line-numbers}
1
2
ISPIN=2
MAGMOM=原子数目*初始磁矩，eg 3*2 4*5 （注：*号前后不要有空格，且对应顺序）

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

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

  2. 有4个Fe原子，如果考虑反铁磁态，则可以设置 MAGMON=3 -3 3 -3 0 0 0 0或者MAGMON=3 3 -3 -3 0 0 0 0

表面优化

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

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

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

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

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

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

MS操作。

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

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

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

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

ASE操作。python

块体优化后
module load python/3.7.3
python

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

大师兄ASE脚本：切金属

python3 cssm.py

vaspkit

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

4 表面弛豫计算

单点能

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

表面优化

1）POSCAR第7行后面加一行加S， Sed –i ‘7a S’ POSCAR
原子坐标加F/T。 Sed –i ’10,11s/$/ F F F/g’ POSCAR
Sed –i ’12,13,14s/$/ T T T/g’ POSCAR
vaspkit –task 402/403 或者ctrl+v,选中列-大写A-输入 T T T - ESC.
2）POTCAR和KPOINTS不变。
INCAR设置：IBRION=1、2，EDIFF=1E-6，EDIFFG=-0.03（比块体小一点），ISIF=2，
提交任务。

5 表面能计算

表面能是创造物质表面时，破坏分子间化学键所需消耗的能量。在固体物理理论中，表面原子比物质内部的原子具有更多的能量，因此，根据能量最低原理，原子会自发的趋于物质内部而不是表面。表面能的另一种定义是，材料表面相对于材料内部所多出的能量。把一个固体材料分解成小块需要破坏它内部的化学键，所以需要消耗能量。如果这个分解的过程是可逆的，那么把材料分解成小块所需要的能量就和小块材料表面所增加的能量相等。但事实上，只有在真空中刚刚形成的表面才符合上述能量守恒。因为新形成的表面是非常不稳定的，它们通过表面原子重组和相互间的反应，或者对周围其他分子或原子的吸附，从而使表面能量降低。
γ=($E_{surf}$ - N * $E_{bulkatom}$)/(2A)

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

吸附能扫描

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

图片

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

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

图片

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

tamas@tamas-PC:~/Desktop/scan$ python scan_adsorption_energy.py -g
Input which atom you want to move to scan adsorption energy–>-1
脚本提示我们移动第几个原子产生不同的吸附结构，我们可以输入97或者-1，移动最后一个Li。这个Cu(111)面是(331)的超胞，因此我们最需要取左下角的1/9的小区域撒点即可，然后程序会将小区域内的数据点进行平移得到整个表面的数据。如何根据不同的表面调整这个参数以及撒点的密度呢？可以使用-u参数进行自定义设置。

Input how many supercells in x axis–>3
Input how many supercells in y axis–>3
Input interporate how many points for each unit cell in x axis,even number are suggested–>2
Input interporate how many points for each unit cell in y axis,even number are suggested–>2
command to submit jobs,e.g. qsub vasp.pbs–>qsub vasp.pbs
Input which atom you want to move to scan adsorption energy–>-1
手动选择x方向的超胞数3，y方向的超胞数3,在小区域内的x方向的撒点数，y方向的撒点数（为了包含中心点，建议选择偶数）和提交排队脚本的命令（我的是qsub vasp.pbs）。脚本会自动产生4个结构（2*2），并会将每个计算所需的文件都复制到相应的文件夹内，并提交排队脚本进行计算。同时脚本会写出一个grid.info里面包含了超胞数，选择移动的原子以及晶格常数信息。 在所有的计算完成后，同样使用该脚本提取数据。

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

图片

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

固定基矢优化结构

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

  1. 采用ISIF=4的固定晶胞体积的方法来计算。这样采用实验参数建的晶胞，由于a和b方向的变化，c轴也会变化。解决了越是优化c轴越小的问题，但是计算量还是很大。

  2. 编译vasp，就自己编译了固定z轴的脚本，专门计算二维材料。这种方法优点是直接换脚本就能解决上面优化时c轴缩小的问题，但是我想固定两个轴就需要再编译一份脚本。比较麻烦。

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

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

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

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

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

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

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

IOPTCELL = xx xy xz yx yy yz zx zy zz

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

1 xx xy xz

2 yx yy yz

3 zx zy zz

1代表驰豫，0代表固定

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

IOPTCELL = 1 1 0 1 1 0 0 0 0

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

滑移能和剥离能

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

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

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

发表于 更新于 分类于

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

Linux常用命令

目录操作

1 切换目录（cd）

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

2 查看目录（ls）

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

3 创建目录（mkdir）

{.line-numbers}
1
2
mkdir tools          //在当前目录下创建一个名为tools的目录
mkdir /bin/tools     //在指定目录下创建一个名为tools的目录

4 删除目录与文件（rm）

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

5 修改目录（mv）

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

6 拷贝目录（cp）

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

7 搜索目录（find）

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

8 查看当前目录（pwd）

{.line-numbers}
1
pwd                         //显示当前位置路径

文件操作

1 新增文件（touch）

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

2 删除文件（rm）

{.line-numbers}rm 文件名link
1
rm -f 文件名           //删除当前目录的的文件（不询问）

3 编辑文件（vi、vim）

{.line-numbers}
1
vi 文件名              //打开需要编辑的文件

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

操作步骤示例

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

补充

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

4 查看文件

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

5 grep

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

文件权限

1 权限说明

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

2 文件权限

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

打包与解压

1 说明

{.line-numbers}
1
2
3
4
.zip、.rar        //windows系统中压缩文件的扩展名
.tar              //Linux中打包文件的扩展名
.gz               //Linux中压缩文件的扩展名
.tar.gz           //Linux中打包并压缩文件的扩展名

2 打包文件

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

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

3 解压文件

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

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

其他准备

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

vaspkit运行的五种方式

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

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

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

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

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

{.line-numbers}
1
2
3
102
2
0.04

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

输入文件

{.line-numbers}
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
INCAR      in    **
STOPCAR    in
stout      out
POTCAR     in     **
KPOINTS    in     **
IBZKPT     out
POSCAR     in     **
CONTCAR    out
CHGCAR     in/out
CHG        out
WAVECAR    in/out
TMPCAR     in/out
EIGENVAL   out
DOSCAR     out
PROCAR     out
OSZICAR    out
PCDAT      out
XDATCAR    out
LOCPOT     out
ELFCAR     out
PROOUT     out

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

下面由简单到复杂介绍

POSCAR

结构文件信息格式如下：

{.line-numbers}
1
2
3
4
5
6
7
8
9
10
11
12
13
14
Mo S2                  #结构名称，可不写或随便写，不能没有这一行
1.0                   # universal scaling parameters，缩放系数
    3.1500000954     0.0000000000     0.0000000000 # 晶格矢量 a(1) 
    -1.5750000477    2.7279801045     0.0000000000 # 晶格矢量 a(2)
     0.0000000000    0.0000000000    12.3000001907 # 晶格矢量 a(3)
   Mo    S                                         # 元素名称，种类
    2    4                                         #对应元素的原子个数，顺序要对应
Direct                                             #
     0.333333345       0.666666689       0.250000000
     0.666666624       0.333333323       0.750000019
     0.333333345       0.666666689       0.620999986
     0.666666624       0.333333323       0.379000014
     0.666666624       0.333333323       0.121000053
     0.333333345       0.666666689       0.878999976
  • 对于自己建模的结构，按照该格式自己写；
  • 对于Materials projiect数据库可以直接下载POSCAR格式；
  • 对于.cif文件，可以导入 VEASTA 等输出为POSCAR的格式，或者用 vaspkit 的411/412功能。

Mind：

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

KPOINTS

文件一般格式如下：

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

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

  • 第三行：VASP只认第一个字母，大小写均可。当然这一行也可以直接写字母G或者g。G表示的是以gamma点为中心生成网格。
    另外一种是原始的Monkhorst-Pack 网格，两者的区别是 M 或者 m 在 G的基础上在三个方向上平移了1/（2N）个单位。G，也叫 gamma centered Monkhorst-Pack Grid；所以，gamma centered 只是MP网格的一种特殊情况。
    因此这一行也可以一直用G不变。

  • 第四行，在xyz三个方向上生成对应数目的K点；__这个一般需要测试__。如果开始不想测试，官网也给出了经验方法，

    • d区金属，k*a ~ 30 Å
    • 普通金属，k*a ~ 25 Å
    • 半导体，k*a ~ 20 Å
    • 绝缘体，k*a ~ 15 Å

    其中a对应三个方向的晶胞参数，也就是对应方向的k点与该方向晶胞参数乘积最好在20-30之间。

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

POTCAR

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

{.line-numbers}
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
PAW_PBE Mo 08Apr2002                   
   6.00000000000000     
 parameters from PSCTR are:
   VRHFIN =Mo: 4p5s4d
   LEXCH  = PE
   EATOM  =   217.5176 eV,   15.9871 Ry

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

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

{.line-numbers}
1
2
3
4
5
6
7
8
9
10
11
12
[scf1355@ln151%nc-l potpaw_PBE.52]$ ls Mo*
Mo:
POTCAR  PSCTR

Mo_pv:
POTCAR  PSCTR

Mo_sv:
POTCAR  PSCTR

Mo_sv_GW:
POTCAR  PSCTR

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

{.line-numbers}
1
cat ~/pot/Mo/POTCAR ~/pot/S/POTCAR  > POTCAR

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

Mind：

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

{.line-numbers}
1
2
3
4
5
[scf1355@ln151%nc-l 2h]$ grep PBE POTCAR 
  PAW_PBE Mo 08Apr2002                   
   TITEL  = PAW_PBE Mo 08Apr2002
  PAW_PBE S 06Sep2000                    
   TITEL  = PAW_PBE S 06Sep2000

INCAR

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

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

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

上面参数详细介绍

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

Mind

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

提交任务

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

Mind

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

查看任务完成

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

输出文件

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

CONTCAR

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

OSZICAR

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

OUTCAR

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

IBZKPT

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

WAVECAR

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

CHGCAR和CHG

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

总结

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

发表于 更新于 分类于

Anaconda安装

  1. Anaconda-Linux版本安装包 【包含了phonopy安装所需的依赖包】
    https://mirrors.tuna.tsinghua.edu.cn/help/anaconda/
    https://www.anaconda.com/download/

  2. 将下载的准备文件上传到Linux系统
    cd 到Anaconda文件的目录下，赋予权限，运行下列命令
    sh [Anaconda文件名]
    例如：
    sh Anaconda3-2019.10-Linux-x86_64.sh

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

    输入conda list 如果显示已安装的包名和版本号便已安装成功。

  3. 修改配置:
    输入vi ~/.bashrc （配置python环境）
    在配置文件末尾中加入下列语句（引号内的路径更换为自己的相对应路径即可）
    export PATH=”/路径/anaconda3/bin:$PATH”

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

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

adaconda使用介绍：

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

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

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

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

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

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

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

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

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

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

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

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

(5) 查看Python版本
pyhton –version

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

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

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

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

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

Phonopy安装

  1. 下载phonopy：https://pypi.org/project/phonopy/#files

  2. 上传解压：tar -xzvf phonopy-2.17.0.tar.gz //解压phonopy文件

  3. cd phonopy-2.17.0
    python setup.py install 安装phonopy
    安装完成后在命令端输入： phonopy 如果出现如下图案则安装成功
    0201

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

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

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

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

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

  • 可以联网：pip install spglib或者spglib-2.2.0

  • 不能联网已安装python：
    下载：https://pypi.org/project/spglib/#modal-close 二进制文件，由于python3.7，所以选择CP37。
    pip install spglib-2.2.0-cp37-cp37m-manylinux_2_17_x86_64.whl

  • 不能联网已安装python：
    下载压缩包解压，同样安装python setup.py install。

  • 我在安装h5py老是报错，
    0202

    安装cached-property。

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

再安装h5py即可。

thirdorder安装

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

ShengBTE安装

  1. 下载：https://drive.google.com/file/d/1Kur3BviQjEby34j0yxbkMF45ADJf2U0J/view

  2. 解压完进入文件夹cp arch.make.example arch.make
    修改前：0203
    修改后：0204

  3. 在Src文件夹内执行命令make

  4. 测试：mpirun -np 16 ShengBTE 2>BTE.err >BTE.out

  5. 添加环境变量：

    vim ~/.bashrc

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

#保存退出

source ~/.bashrc

发表于 更新于 分类于

系统安装

一、下载centos镜像

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

二、制作centos安装盘-u盘

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

    “引导类型选择”选择镜像文件，选择下载好的ISO文件，开始即可。

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

  • u盘格式

    1、FAT32格式：

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

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

三、u盘安装操作系统

BIOS设置问题：

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

其他原因：

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

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

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

安装系统

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

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

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

安装centos不进入桌面

1
2
3
yum  groupinstall  -y  "GNOME Desktop"
reboot now
Startx

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

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

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

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

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

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

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

总结：

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

安装编译器

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

保存退出

安装vasp

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

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

安装vaspkit1.4

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

保存退出。

发表于 更新于 分类于

虚拟机VMware安装（Windows用户）

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

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

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

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

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

视频操作。

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

虚拟机快速安装VMware Tool

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

Intel编译器

过程很简单，一路默认。

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

chmod +x l_HPCKit_p_2022.3.0.8751_offline.sh

chmod +x l_BaseKit_p_2022.3.0.8767_offline.sh.

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

  1. sudo vim ~/.bashrc 中加入 source /opt/intel/oneapi/setvars.sh

  2. 执行which icc; which icpc; which ifort; which mpirun 和echo $MKLROOT 检验是否安装正确。

VASP编译

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

修改makefile.inclue中的内容为：

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

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

source ~/.bashrc

试运行
mpirun -np 4 vasp_std INCAR

发表于 更新于 分类于

查看cpu核数

cat /proc/cpuinfo| grep “physical id”| sort| uniq| wc -l

cat /proc/cpuinfo| grep “cpu cores”| uniq

cat /proc/cpuinfo | grep ‘process’ | sort | uniq | wc -l

系统处理器基本信息记录在/proc/cpuinfo文件中;sort 可针对文本文件的内容，以行为单位来排序;uniq 命令用于检查及删除文本文件中重复出现的行列，一般与 sort 命令结合使用。wc （word count）命令常用于计算文件的行数、字数和字节数，-l , –lines : 显示行数。

查看内存总大小

cat /proc/meminfo | grep “MemTotal” | awk ‘{print $2/1024/1024}’

awk ‘{print $2/1024/1024}’ 以空格分隔，并将分割后的变量打印出来，因为文件里面是KB，转成GB，除两次1024，最后就是内存的大小。

查看硬盘大小

两种方法

第一种（安全，建议使用），所有的相加。

df -h
第二种，不安全,查看磁盘分区，相加。

fdisk -l

fdisk是磁盘分区命令，而且普通用户没法执行这个命令，如果使用出错，很容易完整删除整个系统，虚拟机玩玩就好了。

求第四列数据的平均值和总和：
grep LOOP+ OUTCAR |awk ‘{sum+=$4} END {print “Average=”,sum/NR}’
grep LOOP+ OUTCAR |awk ‘{sum+=$4} END {print “Sum=”,sum}’

发表于 更新于 分类于

常用计算相关数据库

结构查找

  1. The Materials Project https://www.materialsproject.org/
  2. ICSD the Inorganic Crystal Structure Database 无机晶体数据库 http://www2.fiz-karlsruhe.de/icsd_home.html , 需要账号密码
  3. CCDC – The Cambridge Crystallographic Data Centre 剑桥晶体数据库。https://www.ccdc.cam.ac.uk/
  4. COD – 开放晶体结构数据库（Crystallography Open Database）。http://www.crystallography.net/cod/
  5. AMCSD – American Mineralogist Crystal Structure Database。http://rruff.geo.arizona.edu/AMS/amcsd.php
  6. AFLOW：http://www.aflowlib.org/
  7. springer material : https://materials.springer.com/periodictable#
  8. 拓扑材料：（该部分内容参考：http://blog.sciencenet.cn/blog-1502061-1130705.html ）
    Materiae ：http://materiae.iphy.ac.cn/#/
    Topological Materials Arsenal：
    https://ccmp.nju.edu.cn/
    Topological Materials Database:
    http://topologicalquantumchemistry.org/#/
  9. 二维材料：
    Materials Cloud：
    https://www.materialscloud.org/discover/2dstructures/dashboard/ptable
    2D Materials：https://materialsweb.org/twodmaterials
    C2DB: https://cmr.fysik.dtu.dk/c2db/c2db.html
  10. 超导材料：
    第I类超导体：http://www.superconductors.org/Type1.htm
    第II类超导体：http://www.superconductors.org/Type2.htm
  11. 磁结构数据库：MAGNDATA:
    http://webbdcrista1.ehu.es/magndata/index.php?show_db=1
  12. ChemSpider :http://www.chemspider.com/ RSC旗下的一个数据库

物理性质

  1. NIST 数据库 ： https://physics.nist.gov/cuu/Constants/index.html

  2. CRC Handbook of physical chemistry查找热力学参数，晶格常数，entropy等等。
    http://hbcponline.com/faces/contents/ContentsSearch.xhtml

  3. http://www2.ucdsb.on.ca/tiss/stretton/database/inorganic_thermo.htm 无机化合物的物理和热力学相关数据

  4. http://www.knowledgedoor.com/2/elements_handbook/cohesive_energy.html 结合能

  5. http://www.surface-tension.de/ 表面张力

  6. http://www.wiredchemist.com/chemistry/data/thermodynamic-data

  7. https://atct.anl.gov/

  8. 物理常数
    https://physics.nist.gov/cuu/Constants/index.html
    http://wild.life.nctu.edu.tw/class/common/energy-unit-conv-table.html
    http://halas.rice.edu/conversions

  9. 空间群对称性和其布里渊区的信息：
    http://cryst.ehu.es/
    http://www.cryst.ehu.es/cryst/get_kvec.html

软件官网

  1. VASP官网 http://www.vasp.at/

  2. VASP 教程视频 http://www.nersc.gov/users/training/events/3-day-vasp-workshop/

  3. VASP官方论坛 https://cms.mpi.univie.ac.at/vasp-forum/ 需要注册才可以提问

  4. Henkelman 课题组 http://theory.cm.utexas.edu/henkelman/

  5. ASE： Atomic Simulation Environment： https://wiki.fysik.dtu.dk/ase/

  6. Pymatgen： http://pymatgen.org/, 适合材料计算的

  7. K-path ： https://www.materialscloud.org/work/tools/seekpath

  8. Chemml: https://hachmannlab.github.io/chemml/

  9. p4vasp：http://www.p4vasp.at/

  10. RDkit：https://www.rdkit.org/

  11. vaspkit ： https://sourceforge.net/projects/vaspkit/files/ VASP计算前后处理。

  12. VESTA ： http://jp-minerals.org/vesta/en/ 日本开发的可视化建模软件

  13. Xcrysden ：http://www.xcrysden.org/， 用于批量做图，选择K-path做能带图，等。

  14. Jmol: (sourceforge.net)

  15. phonopy：https://atztogo.github.io/phonopy/

  16. Lobster：http://schmeling.ac.rwth-aachen.de/cohp/index.php?menuID=6

  17. Wannier90：http://www.wannier.org/

免费的学习网站和论坛

  1. 思想家公社的门口：量子化学·分子模拟·二次元 (sobereva.com)

  2. 计算化学公社 - 高水平计算化学、理论化学交流论坛 (keinsci.com)

  3. Learn VASP The Hard Way (bigbrosci.com)大师兄科研网

  4. 世事如棋 (https://blog.shishiruqi.com/)

  5. Physics (yh-phys.github.io)

  6. https://www.guanjihuan.com/

  7. 学术之友知乎 https://zhuanlan.zhihu.com/gangtangdft

  8. 一个人就是一个叠加态（Li pai）：http://blog.sina.com.cn/lipai91

  9. http://blog.sciencenet.cn/home.php?mod=space&uid=1502061

  10. http://blog.sciencenet.cn/home.php?mod=space&uid=2909108

  11. http://blog.sciencenet.cn/home.php?mod=space&uid=548663

  12. http://blog.sciencenet.cn/home.php?mod=space&uid=3352196

  13. http://blog.sciencenet.cn/home.php?mod=space&uid=3102863

  14. http://blog.sciencenet.cn/u/zhangfrank

  15. http://blog.sciencenet.cn/home.php?mod=space&uid=1524047

  16. http://blog.sciencenet.cn/home.php?mod=space&uid=3222255

  17. http://blog.sciencenet.cn/home.php?mod=space&uid=3388193

  18. http://blog.sciencenet.cn/home.php?mod=space&uid=478347

  19. 叶小球：http://blog.sciencenet.cn/home.php?mod=space&uid=567091

  20. Nan Xu github：https://github.com/tamaswells

  21. QijingZheng github：https://github.com/QijingZheng

  22. Chenists github：https://github.com/Chenists/VASP

  23. ChengCheng Xiao：https://github.com/Chengcheng-Xiao?tab=repositories

  24. ShaoZhengjiang github：https://github.com/PytLab

公众号

B站up

程序检索（更新至2024/07/21）

1、Sym4state.jl 程序：磁性材料的高效计算包
文章题目：

MagneticTB: A package for tight-binding model of magnetic and non-magnetic materials

文章链接：

https://doi.org/10.1016/j.cpc.2024.109283

附件下载：

https://github .com /A -LOST -WAPITI /Sym4state .jl

2、Crystal Toolkit

该项目由Materials Project (https://materialsproject.org/) 员工Matthew Horton主导开发，具有很多很方便的功能，比如可视化晶体结构、原胞和单胞之间的转换、切一个slab、构建晶界、猜测元素氧化态、替换和取代、分析局域配位环境、生成XRD图谱、计算的能带结构和态密度、相图等等。这些功能中有些甚至可以替代Materials Studio对应的功能(可视化晶体结构、单胞和原胞转换、构建slab、计算XRD图谱等)，也有一些特色的功能(分析局域配位环境等)。

Crystal Toolkit主页

https://crystaltoolkit.org/

Crystal Toolkit GitHub主页

https://github.com/materialsproject/crystaltoolkit

3、吸附位点扫描

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

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

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

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

4、nMoldyn-3

We present a new implementation of the program nMoldyn, which has been developed for the computation and decomposition of neutron scattering intensities from Molecular Dynamics trajectories
./nMOLDYNStart.py

https://github.com/khinsen/nMOLDYN3/

  1. Modemap

https://github.com/JMSkelton/ModeMap/tree/master

典型的用法包括三到四个步骤，实施一系列短程序：

  • 准备一系列结构，这些结构沿着一个（1D映射）或两个（2D映射）声子模式，在一系列的幅度上“冻结”（ModeMap.py）。

  • 在这些结构上运行单点能量计算，并提取总能量。提供了一个适用于维也纳第一性原理模拟软件包（VASP）代码的基本脚本（ExtractTotalEnergies.py）。

  • 对计算结果进行后处理，以生成势能表面图（ModeMap_PostProcess.py）。

  • 如果需要，可以将计算出的势能表面拟合成多项式函数，以便进一步分析（ModeMap_PolynomialFit.py）。

  • 再次，如果需要，通过使用拟合的势能面作为输入，进行进一步的后处理，使用1D薛定谔求解器。由J. Buckeridge（UCL）编写的Fortran代码在1DSchrodingerSolver中提供。该代码使用傅里叶方法获得非简谐势中的特征值和特征向量，并确定一个有效的重正化（简谐）频率，以再现其在给定温度下对热力学分区函数的贡献。由J. M. Frost编写的TISH代码可以与多项式拟合一起工作，用于研究带隙变形势能。

限制条件

  1. Twister

https://github.com/qtm-iisc/Twister/blob/main/README.md

TWISTER 是一个 Python 包，它帮助你在引入扭曲的二维材料之间找到和构造相称的莫尔超晶格。

TWISTER 还可以通过 LAMMPS 接口研究使用经典力场的莫尔超晶格的结构重构。请查看 FFRelaxation/ 以获取示例。

  • 转角双层磷烯。

从包含输入文件 get_ang.inp 的目录中，在终端运行以下命令：
python path_to_Get_Ang/SRC/get_ang.py

输入文件的文档可以在 Documentation/ 文件夹中找到。

为了找到适当的扭转角，运行 get_ang.inp 以获取一系列扭转角。

  1. 如果生成了输出文件，它们将具有没有不匹配的适当扭转角。

  2. 会生成一个 twist.inp 文件，包含新的晶格向量和最小变形及最小超晶格面积的扭转角。运行 twister.py 获取超晶格中的原子位置。

每个角度对应的输出文件将以 “Angle_%angle%” 的形式写入 solutions/ 文件夹中。
输出的文档可以在 Documentation/ 文件夹中找到。

为了可视化生成的结构，我们提供一个工具 tovasp.py，它读取 twist.inp 和 superlattice.dat 以生成 POSCAR.vasp。可以直接在 VESTA （http://jp-minerals.org/vesta）中进行可视化。
python path_to_Twister/SRC/tovasp.py

  1. Phasego

附件下载：http://dx.doi.org/10.17632/xp7gvmrxrg.1http://cpc.cs.qub.ac.uk/summaries/AEVQ_v2_0.html

“Phasego 包从第一性原理计算获得的声子态密度中提取亥姆霍兹自由能。在状态方程拟合的帮助下，它在固定温度/压力下降低了作为压力/温度函数的吉布斯自由能。基于准谐波近似（QHA），它计算所有感兴趣结构之间可能的相边界，最后自动绘制相图。对于单相分析，Phasego 可以从数值上推导出许多属性，例如热膨胀系数、体积模量、热容、热压、Hugoniot 压力-体积-温度关系、Gruneisen 参数和 Debye 温度。为了检验它的相变分析能力，我在这里举两个例子：半导体GaN和金属Fe。对于 GaN，Phasego 自动确定并绘制了所提供的闪锌矿 (ZB)、纤锌矿 (WZ) 和岩盐 (RS) 结构之间的相边界。在 Fe 的情况下，结果表明，在高温下，电子热激发自由能校正显着改变了体心立方 (bcc)、面心立方 (fcc) 和六方密堆积 (hcp) 之间的相界结构

0%