首页 >> 仿真软件说明 >>ABACUS >>朱雪刚版 >> 2.2.1.1 pw中ecut测试 2024.03.12，ABACUS3.4.0测试结果（软件已更新。推荐使用最新版本）

详细内容

2.2.1.1 pw中ecut测试 2024.03.12，ABACUS3.4.0测试结果（软件已更新。推荐使用最新版本）

先附上测试脚本，后续并对脚本进行后续解释。2024.03.12

脚本执行命令是bash script.sh 即“bash”空格“脚本名称”来执行（应该是视频教学演示效果更好）

#!/bin/bash
#source activate abacus_env
OMP_NUM_THREADS=1
cpu_cores=(6 8 10 12)
upf_potential_file="Be_ONCV_PBE-1.0.upf"
path_to_potential_file="../../../../SG15_v1.0_ONCVPP_and_Orbitals"
orbital_file="../../../../SG15_v1.0_ONCVPP_and_Orbitals"

ecut_energys=(30 50 70 80 90 100 110 120 130 150 200) # the unit is Rydberg and 1 Rydberg = 13.6 eV
KPTs=(8 1 2 3 4 5 6 7 8 9 10 12 14 16 18 20)
Rcut_aus=(10 7 8 9)
num_orbs=(6s2p 5s2p 4s1p 4s1p_2016 2s)
strus=(16 2 16 54 128 250)

test_var=pw_ecut

cat > ${test_var}_energy_time << EOF
${test_var}    energy/atom    Time
EOF
echo "start to test ${test_var} with energy and time"
for ecut_energy in ${ecut_energys[@]} 
    do
        echo "strus $strus Radius cutoff ${Rcut_aus} energy cutoff ${ecut_energy} orbital ${num_orbs} KPT ${KPTs} ${KPTs} ${KPTs}"
        mkdir "${test_var}_${ecut_energy}"
        cd "${test_var}_${ecut_energy}"
            echo "now we are in the dir of $PWD "
            echo "write the INPUT"

            cat > INPUT << EOF
INPUT_PARAMETERS
#Parameters (General)

pseudo_dir      $path_to_potential_file
orbital_dir     $orbital_file
#Parameters (Accuracy)
basis_type      pw
ecutwfc         $ecut_energy
smearing_method gauss
smearing_sigma  0.01
calculation     scf
scf_nmax        100

EOF

            echo "write the KPT"
            cat > KPT << EOF
K_POINTS
0
Gamma
$KPTs $KPTs $KPTs 0 0 0
EOF
            echo "get the stru file"
            cp ../../../../STRU/Be/STRU_file/Be${strus}.STRU ./STRU
            sed -i "5s/Be_gga_[0-9]*au_100Ry_[0-9a-zA-Z]*.orb/Be_gga_${Rcut_aus}au_100Ry_${num_orbs}.orb/" STRU
            echo "run abacus in $PWD"
            OMP_NUM_THREADS=1 mpirun -n $cpu_cores abacus | tee out.log
            Time=`grep "TOTAL" out.log  | awk '{printf int($4)}'`
            energy_atom=`grep "FINAL_ETOT_IS" OUT.ABACUS/running_scf.log  | awk -v stru="${strus}" '{printf "%.8f\n",$2/stru}'`
            echo "$ecut_energy    $energy_atom    $Time" >> ../${test_var}_energy_time
        cd ..
        echo "done the $ecut_energy with energy scf calculation"
done

cat > ${test_var}_energy_time.py << EOF
import matplotlib.pyplot as plt

with open('${test_var}_energy_time', 'r') as f:
    data = f.readlines()

${test_var} = []
energy_atom = []
Time = []

for line in data[1:]:
    line = line.strip().split()
    ${test_var}.append(str(line[0]))
    energy_atom.append(float(line[1]))
    Time.append(int(line[2]))

fig, (ax1, ax2) = plt.subplots(2, 1, sharex=True)
ax1.plot(${test_var}, energy_atom, marker='o')
ax1.set_ylabel('energy_atom')
ax1.set_title('Plot of ${test_var} vs energy/atom')
ax1.set_ylim(min(energy_atom), max(energy_atom))

ax2.plot(${test_var}, Time, marker='o')
ax2.set_xlabel('${test_var}')
ax2.set_ylabel('Time')
ax2.set_title('Plot of ${test_var} vs Time')
ax2.set_ylim(min(Time), max(Time))
plt.savefig('${test_var}_energy_time.png')

EOF
python ${test_var}_energy_time.py
echo "all finish"

bash脚本具体解释 2024.03.12

现在假设直接在linux下的电脑上运行ABACUS是可行的，即已经安装好ABACUS，且输入abacus会出现相应的版本信息等（注释:如果通过conda安装ABACUS则需要先用conda activate abacus_env激活对应的环境，其中abacus_env是安装时设定的环境名称）。

#!/bin/bash

是告诉机器所使用的是bash命令。（如果此处改成“#!/usr/bin/env python”则告诉脚本使用python）。

#source activate abacus_env

使用#开头表示该行是注释行即可以忽略，由于刚开始使用的是conda安装的ABACUS脚本，因此也可以在bash中激活对应安装着ABACUS的环境，在脚本中使用source。

upf_potential_file="Be_ONCV_PBE-1.0.upf"

#用变量先指定使用的势函数名称。

path_to_potential_file="../../../../SG15_v1.0_ONCVPP_and_Orbitals"

#用变量指定势函数所在的文件夹路径，这里使用的相对路劲，可以改成绝对路劲（linux下文件的位置用路劲表示，例如绝对路径/home/Be/a.txt，表示在home文件夹下有个Be文件夹，Be文件夹下有个a.txt文件。另外../表示上一层文件夹，../../表示上一层的上一层文件夹）

orbital_file="../../../../SG15_v1.0_ONCVPP_and_Orbitals"

#用变量指定原子轨道势函数所在的文件夹路径

ecut_energys=(30 50 70 80 90 100 110 120 130 150 200) # the unit is Rydberg and 1 Rydberg = 13.6 eV

此处设定ecut_energys变量，变量后紧跟等号然后紧跟小括号，括号内是要测试的截断能，同时也注释了在ABACUS中阶段能的单位是Rydberg，且1 Rydberg = 13.6 eV

KPTs=(8 1 2 3 4 5 6 7 8 9 10 12 14 16 18 20)

该变量是设置KPT文件中K点设置的变量，在这里第一个数值是8，即当调用KPTs变量时获取的是第一个数值，如果要获取后续数值需要使用bash中的变量规则，后续用到再提及。

Rcut_aus=(10 7 8 9)

该变量是对截断半径设置变量，截断半径越长计算所需要的时间也越长。

num_orbs=(6s2p 5s2p 4s1p 4s1p_2016 2s)

该变量是对原子轨道基组lcao计算相关轨道的选择做变量，数值越大，字母越多表示使用的原子轨道越详细，即算计的越准确，但时间也越长

strus=(16 2 16 54 128 250)

该变量是针对Be结构进行设置，不同体系包含不同的原子数目，这里把原子数目单独拿出来。

test_var=pw_ecut

先定义一个变量名，后续使用该变量名生成相关文件

cat > ${test_var}_energy_time << EOF
${test_var}    energy/atom    Time
EOF

这几行是先定义一个变量名，后续使用该变量名生成相关文件

然后用cat把EOF到第二个EOF之间的文字写到${test_var}_energy_time文件中。其中${test_var}，就是上边定义的名称。

即bash可通过cat > a来创建一个名称为a的新文件，<<EOF…………EOF是把之间的数据写到a文件夹内。

echo "start to test ${test_var} with energy and time"

由于脚本执行过程并不输出内容，为了确保脚本是执行的，这里在屏幕上用echo输出一行信息。

for ecut_energy in ${ecut_energys[@]} 
    do
        echo "strus $strus Radius cutoff ${Rcut_aus} energy cutoff ${ecut_energy} orbital ${num_orbs} KPT ${KPTs} ${KPTs} ${KPTs}"
        mkdir "${test_var}_${ecut_energy}"
        cd "${test_var}_${ecut_energy}"
            echo "now we are in the dir of $PWD "
            #(注意 for循环这里还没结束)

这几行开始使用变量进行循环“for ecut_energy in ${ecut_energys[@]}”是做一个for循环，即用一个新的变了ecut_energy（这个变量比前边的少个s，这里多一个s表示变量内包含多个值）循环读取ecut_energys中的值。

下一行do是执行的意思，do的结束是done（在后边的脚本中）。

再下一行echo是输出包含变量的信息，即结构、截断、能量、原子轨道、KPT设置等。

mkdir是在当前目录下创建文件，此处使用变量"${test_var}_${ecut_energy}"来命名，即首先创建pw_ecut_30文件夹，后续循环过程分别创建pw_ecut_50、pw_ecut_70、pw_ecut_80等文件夹。

cd是打开文件夹，即打开刚创建的pw_ecut_30文件夹

echo输出现在所做文件夹的路径，其中 "now we are in the dir of $PWD "中的“$PWD”是linux下默认的变量，会给出当前文件夹的路径。

下边这段bash命令还在for循环里边主要是创建INPUT文件

#(注，接上边的for循环内容)
            echo "write the INPUT"
            cat > INPUT << EOF
INPUT_PARAMETERS
#Parameters (General)
pseudo_dir      $path_to_potential_file
orbital_dir     $orbital_file
#Parameters (Accuracy)
basis_type      pw
ecutwfc         $ecut_energy
smearing_method gauss
smearing_sigma  0.01
calculation     scf
scf_nmax        100
EOF
#（注，for循环还没结束，需要继续往下看）

echo 在屏幕输出，要写INPUT文件了。

cat > INPUT << EOF…………EOF，跟上边解释相同，把EOF和EOF之间的内容写到INPUT文件夹下

INPUT_PARAMETERS是ABACUS中INPUT文件的提示符（该项我不是很确定）

#Parameters (General) #注释行，表示开始写一些基本的命令

pseudo_dir $path_to_potential_file #使用变量指定upf文件夹所在的文件夹

orbital_dir $orbital_file #使用变量指定lcao文件夹所在的文件夹

#Parameters (Accuracy) # 注释行，表示下边对准确性进行设置，主要是阶段能设置

basis_type pw #使用的势函数形式，pw表示平面波基组，lcao表示原子轨道基组，这里使用pw。

ecutwfc $ecut_energy #截断能设置，这里使用变量$ecut_energy即第一次循环取30，然后再取50、70等

smearing_method gauss #参数我不是很多先不解释，对于Be晶体来说，或者金属来说加上之后更容易收敛

smearing_sigma 0.01 #该处的0.01是参考的ABACUS文件夹下exzample中的设置，具体怎么设置准确我不清楚。

calculation scf #设定计算的方法，scf表示自洽性计算，即结构文件中原子位置不动，对波函数进行初猜，同时优化波函数至当前构型的能量最低点。除了scf还有md、relax、cell-relax等，md是进行第一性原理动力学模拟、relax是进行原子位置移动并寻找能量最低点、cell-relax是对原子坐标和晶格进行移动和优化

scf_nmax 100 #在自洽迭代时，迭代的步数，即优化波函数的最大次数，如果十几步就优化好了则只优化十几步，但如果一致优化不到程序默认的或我们设定的标准（此处用程序默认值1e-7）则一直优化，直到达到100步时退出，并在输出文件夹给出不收敛的提示。

写完INPUT之后开始准备KPT文件

#(注继续上边的for循环)
            echo "write the KPT"
            cat > KPT << EOF
K_POINTS
0
Gamma
$KPTs $KPTs $KPTs 0 0 0
EOF
#(for循环还未结束)

这几行时准备KPT文件，如果INPUT里边设置了kspacing则不需要准备KPT文件

echo输出信息，开始写KPT

cat > KPT << EOF…………EOF把相关内容写到KPT文件中

K_POINTS #KPT内的具体格式，不熟悉先不解释

0 #KPT内的具体格式，不熟悉先不解释

Gamma #KPT内的具体格式，不熟悉先不解释

$KPTs $KPTs $KPTs 0 0 0 #使用变量KPTs中的第一个值。

下边是获取结构文件，同时对结构文件进行相应的修改

#(注继续上边的for循环)
            echo "get the stru file"
            cp ../../../../STRU/Be/STRU_file/Be${strus}.STRU ./STRU
            sed -i "5s/Be_gga_[0-9]*au_100Ry_[0-9a-zA-Z]*.orb/Be_gga_${Rcut_aus}au_100Ry_${num_orbs}.orb/" STRU
            #(for循环还未结束)

echo先输出该部分的操作

cp复制“前、前、前、前一个文件夹下，STRU文件夹下，Be文件夹下，STRU_file文件夹下，的Be16.STRU文件到当前文件夹下，并命名成STRU

sed是对文件夹下的内容进行修改，即对第五行的内容进行对应原子轨道的修改，改行命令解释可咨询chatgpt。

然后运行软件ABACUS，并提取关心的物理量，由于也是第一次写脚本，因此脚本只提取了能量做对比。

#(注继续上边的for循环)
            echo "run abacus in $PWD"
            OMP_NUM_THREADS=1 mpirun -n $cpu_cores abacus | tee out.log
            Time=`grep "TOTAL" out.log  | awk '{printf int($4)}'`
            energy_atom=`grep "FINAL_ETOT_IS" OUT.ABACUS/running_scf.log  | awk -v stru="${strus}" '{printf "%.8f\n",$2/stru}'`
            echo "$ecut_energy    $energy_atom    $Time" >> ../${test_var}_energy_time
        cd ..
        echo "done the $ecut_energy with energy scf calculation"
done
#(for循环结束)

echo输出当前的文件夹路径，以便我们确保所在位置是准确的。

OMP_NUM_THREADS=1 mpirun -n $cpu_cores abacus | tee out.log#运行ABACUS的命令，一般是“OMP_NUM_THREADS=1 mpirun -n 16 abacus”这里使用变量是为了后续测试cpu核数对计算的影响。从我后续测试来看，cpu核数为物理核数时计算最快，例如E52680V2有两颗cpu，单个cpu是10物理核心2个线程，即总共20核心，但此时不是我把所有核心用完最快，反而我使用“OMP_NUM_THREADS=1 mpirun -n 10 abacus”是最快的。

|tee out.log是把bash脚本中运行的命令从屏幕输出出来，该命令是我咨询chatgpt获得的。

Time=‘……’用grep命令从out.log文件夹下获取包含TOTAL行的内容，此行第五个是程序有运行时间，用awk获取对应时间

energy_atom=‘……’用grep获取OUT.ABACUS文件夹下running_scf.log文件中包含“FINAL_ETOT_IS”的行，改行有结构的总能量，然后用awk获取第二个字符数值，同时把总能量除以原子个数换算成单个原子的能量。

echo "$ecut_energy $energy_atom $Time" >> ../${test_var}_energy_time#这里有个>>重定向，即把上边获取的原子能量，原子数目，和模拟时间等数值输出到前边创建的${test_var}_energy_time文件下。

cd ..进入上一层目录，因为循环时我们进入了一个创建的目录，这里返回去，再做其他阶段能的文件的计算。

echo " 输出完成了30能量截断的测试。

done 表示上边使用的for循环结束。

脚本运行至此，测试所需的数据已存储在文件pw_ecut_energy_time中，下边采用python进行作图。python脚本也可以通过bash进行添加：

cat > ${test_var}_energy_time.py << EOF
import matplotlib.pyplot as plt
with open('${test_var}_energy_time', 'r') as f:
    data = f.readlines()
${test_var} = []
energy_atom = []
Time = []

for line in data[1:]:
    line = line.strip().split()
    ${test_var}.append(str(line[0]))
    energy_atom.append(float(line[1]))
    Time.append(int(line[2]))

fig, (ax1, ax2) = plt.subplots(2, 1, sharex=True)
ax1.plot(${test_var}, energy_atom, marker='o')
ax1.set_ylabel('energy_atom')
ax1.set_title('Plot of ${test_var} vs energy/atom')
ax1.set_ylim(min(energy_atom), max(energy_atom))

ax2.plot(${test_var}, Time, marker='o')
ax2.set_xlabel('${test_var}')
ax2.set_ylabel('Time')
ax2.set_title('Plot of ${test_var} vs Time')
ax2.set_ylim(min(Time), max(Time))
plt.savefig('${test_var}_energy_time.png')

EOF

同样使用cat> filename << EOF …………EOF进行文件的创建和内容的添加

import matplotlib.pyplot as plt#画图使用的python库

with open('${test_var}_energy_time', 'r') as f:

data = f.readlines() #这两行主要是读取文件数据，并按行读取

${test_var} = [] #定义保存测试量的数组（变量存储相关数据）

energy_atom = [] #定义单个原子能量的变量

Time = [] #定义程序执行时间的变量

for line in data[1:]:

line = line.strip().split() #使用for 循环利用line变量读取data中的数据，[1:]是python中的切片，具体理解可以问chatgpt。

${test_var}.append(str(line[0])) #把截断半径的数据加入到相关变量中。

energy_atom.append(float(line[1])) #把原子能量加入相关变量中

Time.append(int(line[2])) #把使用时间存储在相关变量中

fig, (ax1, ax2) = plt.subplots(2, 1, sharex=True) #python中画图，沿着x轴包含两个子图

ax1.plot(${test_var}, energy_atom, marker='o') #第一个图是测试截断和能量的关系

ax1.set_ylabel('energy_atom') #y轴标签是energy_atom 注：这里少了能量单位，应该是 energy （ev/atom），刚开始写的脚本出点儿错误是正常的，望理解。

ax1.set_title('Plot of ${test_var} vs energy/atom') #图的标题，即x轴标签

ax1.set_ylim(min(energy_atom), max(energy_atom)) #y轴的范围，上下分别取能量的最大和最小值。

ax2.plot(${test_var}, Time, marker='o') #第二个图画测试量和时间关系

ax2.set_xlabel('${test_var}') #x轴的标签

ax2.set_ylabel('Time')# y轴的标签

ax2.set_title('Plot of ${test_var} vs Time') #图片标题

ax2.set_ylim(min(Time), max(Time)) #y轴上下限

plt.savefig('${test_var}_energy_time.png') #保存图片

写完python脚本之后，利用在bash脚本中利用python *.py可执行python脚本

python ${test_var}_energy_time.py #执行python脚本
echo "all finish"

并最终输出完成测试“all finish”。

预期的效果

运行脚本用

bash ecut_convergence.sh

如果一切顺利，即可在服务器上相应文件夹下看到用python做的图，以及对应的数据信息。如下：

测试结果

结构文件Be16 atoms

pw_ecut energy/atom Time

30 -361.60480453 96

50 -362.05282097 183

70 -362.05631545 316

80 -362.05648695 426

90 -362.05654260 580

100 -362.05664713 734

110 -362.05673714 894

120 -362.05677717 931

130 -362.05679495 1080

150 -362.05681318 1483

200 -362.05684192 2396

图中由于30截断较低是能量跟50之后的数据偏差太大，因此作图时最好把30的数据去掉，然后再作图对比，注可问chatgpt修改python脚本（for line in data[1:]:）修改为（for line in data[2:]: ）。

从测试结果可以看出当能量截断ecut设定在100以上时能量的偏差就降低到了0.1meV以上，从我听说的标准上来看是可行的，因此后续我的计算能量截断选择为“ecutwfc 100” 。

存在问题

由于本人也是初次使用bash脚本，因此在做测试时只测试了能量。其实还需要测试力，EOS等（目前EOS可以用DPGEN来完成），当然可能还需要测试其他的量。如果知道测试那些量，就可以尝试写相关脚本，目前我的问题是不知道应该测试那些量，另外编写脚本比较费时间，如果大家能把使用的脚本分享出来，可将降低其他人的入门时间，也必定会从他人的分享中获益。

截止2024年6月29日，ABACUS开发者已经开始出收敛性测试的处理python脚本《ABACUS中的ecutwfc和kspacing的收敛性测试》，不过该脚本目前也只有能量测试，需要使用的研究人员可以参考修改。

简化版本的bash脚本

由于以上脚本是经过多次测试之后写的脚本，变量较多，不利于新手理解，因此提供一个简单只测ecut的脚本，但该脚本不再解释。由于脚本中没有复制STRU和*.upf文件,因此需要先把STRU文件和势函数文件放到bash脚本的相同文件夹下。然后运行bash pw_ecut_test_siplife.sh

pw_ecut_test_siplife.sh如下

#!/bin/bash
ecut_energys=(30 50 70 80 90 100 110 120 130 150 200) # the unit is Rydberg and 1 Rydberg = 13.6 eV

cat > ecut_energy_time << EOF
ecut    energy/atom    Time
EOF
echo "start to test ecut with energy and time"
for ecut_energy in ${ecut_energys[@]} 
    do
        echo "do the energy cutoff ${ecut_energy} test"
        mkdir "ecut_${ecut_energy}"
        cd "ecut_${ecut_energy}"
            echo "now we are in the dir of $PWD "
            echo "write the INPUT"

            cat > INPUT << EOF
INPUT_PARAMETERS
#Parameters (General)

pseudo_dir      ../
orbital_dir     ../
#Parameters (Accuracy)
basis_type      pw
ecutwfc         $ecut_energy
smearing_method gauss
smearing_sigma  0.01
calculation     scf
scf_nmax        100

EOF

            echo "write the KPT"
            cat > KPT << EOF
K_POINTS
0
Gamma
6 6 6 0 0 0
EOF

            echo "run abacus in $PWD"
            cp ../STRU ./STRU
            OMP_NUM_THREADS=1 mpirun -n 6 abacus | tee out.log
            Time=`grep "TOTAL" out.log  | awk '{printf int($4)}'`
            energy_atom=`grep "FINAL_ETOT_IS" OUT.ABACUS/running_scf.log  | awk -v stru="16" '{printf "%.8f\n",$2/stru}'`
            echo "$ecut_energy    $energy_atom    $Time" >> ../ecut_energy_time
        cd ..
        echo "done the $ecut_energy with energy scf calculation"
done

cat > ecut_energy_time.py << EOF
import matplotlib.pyplot as plt

with open('ecut_energy_time', 'r') as f:
    data = f.readlines()

ecut = []
energy_atom = []
Time = []

for line in data[2:]:
    line = line.strip().split()
    ecut.append(str(line[0]))
    energy_atom.append(float(line[1]))
    Time.append(int(line[2]))

fig, (ax1, ax2) = plt.subplots(2, 1, sharex=True)
ax1.plot(ecut, energy_atom, marker='o')
ax1.set_ylabel('energy_atom')
ax1.set_title('Plot of ecut vs energy/atom')
ax1.set_ylim(min(energy_atom), max(energy_atom))

ax2.plot(ecut, Time, marker='o')
ax2.set_xlabel('ecut')
ax2.set_ylabel('Time')
ax2.set_title('Plot of ecut vs Time')
ax2.set_ylim(min(Time), max(Time))
plt.savefig('ecut_energy_time.png')

EOF
python ecut_energy_time.py
echo "all finish"

作者：朱雪刚邮箱:xuegangzhu@qq.com；工作单位：石家庄学院理学院/北京科学智能研究院（AISI）访问学者2023.07-2024.09，访问导师北京大学陈默涵; 徐张满仓邮箱: xuzhangmancang@dp.tech

截止2024.07.21录制视频教程已上传至Bohrium的课程《DeePMD应用案例讲解：铜原子掺杂铍晶体的机器学习势函数拟合过程演示》网址： https://bohrium.dp.tech/courses/1075495070?tab=courses 后续会在Bohrium平台更新

注意：后续的更新，大部分会在Bohrium课程平台进行，请看教程入门的同学加入课程进行学习；且录制的视频课程会把一些个人观点给加入，从个人观点来看比文字教程的内容更多

写教程内容讨论QQ群:143276924 DPGEN+ABACUS教程准备；

ABACUS软件的QQ群：759914681，群内有专职开发人员，目前ABACUS提问问题以github的issue为主，群内可作为辅助提问。

欢迎大家推广本教程，让更多的dp入门学习者有个参考，目前2024.09.16我联系的微信公众号推广是 lammps加油站的小马老师（我也报名了小马老师的一对一辅导）。当然也看到lammps爱好者在转发推广。感谢大家公众号的推广。

上一篇2.2.1.2 KPT测试（做DPGEN要用上kspaci下一篇2.2.1 ABACUS对pw 基组进行测试。2024.0

2.2.1.1 pw中ecut测试 2024.03.12，ABACUS3.4.0测试结果（软件已更新。推荐使用最新版本）

先附上测试脚本，后续并对脚本进行后续解释。2024.03.12

bash脚本具体解释 2024.03.12

预期的效果

存在问题

简化版本的bash脚本