新建
ABACUS对比CP2K精度和效率测试 | Si的状态方程（EOS）
John Yan
kirk0830
weiqingzhou@whu.edu.cn
推荐镜像 :abacus-cp2k-elk:v1.1.1
推荐机型 :c16_m32_cpu
赞 2
12
目录
©️ Copyright 2024 @ Authors
作者：John Yan (AISI电子结构团队实习生) 📨
日期：2024-7-12
共享协议：本作品采用知识共享署名-非商业性使用-相同方式共享 4.0 国际许可协议进行许可。
快速开始：点击上方的 开始连接 按钮，选择 abacus-cp2k-elk:v1.1.1 镜像及 c16_m32_cpu 节点配置，使用其他节点配置可能会产生计算结果和时间无法复现，并在连接后点击右上方的切换内核，选择Python (ipykernel)内核稍等片刻即可运行。
代码
文本
目录：
1. 对比工作流程介绍
2. 准备文件介绍
3. 使用ELK进行Si的EOS的全电子计算
4. 使用CP2K基于GTH-PBE赝势和DZVP-GTH-q4基组进行Si的EOS计算
5. 使用ABACUS基于由CP2K的GTH-PBE赝势生成的赝势进行Si的EOS计算
6. ABACUS和CP2K在Si体系EOS上的精度和效率对比总结
代码
文本
0. 写在前面
Kohn-Sham密度泛函理论计算的关键一步就是选择合适的基组对KS方程进行参数化，将其转换为本征值求解问题。不同的基组选择，将极大影响计算的精度和效率。最为常见的选择是平面波基组，优点是体系精度可以随着平面波数量系统性增加，但是巨大的基组数量也限制了可计算体系的尺寸大小，一般限制在1000个原子以内。选择局域基组可以显著的增加可计算体系的大小，在一定精度的保证下，大幅提升计算的效率。计算软件 CP2K支持局域的高斯基组GTOs来进行计算，广泛地应用于计算化学领域。而国产开源密度泛函理论计算软件ABACUS除了平面波基组之中，还良好支持数值原子轨道基组。这一系列的教程通过使用ABACUS的LCAO基组和CP2K的GTOs基组的结果进行系统的测试，我们相信这可以帮助熟悉CP2K的用户快速建立起对ABACUS计算结果的认知和信心。
代码
文本
1. 对比工作流程介绍
代码
文本
对比的精度指标：以全电子计算的Si体系的EOS作为精度对标参考，将cp2k和abacus计算的eos与其进行比较，计算delta数值，数值越大表明精度越低。
对比的效率指标：完成计算的时间，cp2k用时读取自output.log中Timing-Total time，abacus用时读取自log中的TIME STATISTICS-TIME-total
Step1:构建精度对标参考：使用ELK进行Si的EOS计算（包括k点收敛性测试）
Step2:使用cp2k常用的赝势和基组进行Si的EOS计算（包括k点、cutoff、rel_cutoff收敛性测试）
Step3:基于cp2k使用的赝势生成abacus可以读取的upf赝势、并基于此赝势生成轨道文件（包括ecutwfc收敛性测试）
Step4:基于以上赝势和轨道文件分别进行ABACUS的Si的eos的pw和lcao计算（包括k点收敛性测试）
Step5:总结精度和效率指标进行比较
代码
文本
2. 准备文件介绍
代码
文本
为了确保ABACUS、CP2K和ELK计算的结果具有可比较性，需要对参数进行对齐，目前总结了以下关键参数需要对齐：
任务类型 scf计算对应的参数
elk：tasks 0
cp2k：RUNTYPE ENERGY
abacus：calculation scf
晶胞信息： 统一通过ATOMKIT工具由相应的cif文件生成
泛函类型： pbe泛函对应的参数
elk：xctype 20
cp2k：XC_FUNCTIONAL PBE
abacus：dft_functional PBE
k点： 根据k点收敛性测试结果进行对齐
elk：ngirdk 7, 7, 7
cp2k：KPOINT 7, 7, 7
abacus：kspacing 0.10 (对应的k点为7, 7, 7)
新旧电荷密度混合方法： broydpm方法对应的参数
elk：mixtype 3
cp2k：&MIXING METHOD BROYDEN_MIXING
abacus：mixing_type broyden
此外还有自旋极化、自旋轨道耦合等参数，可以参考官网的关键词手册，这里不进行详细介绍。
代码
文本
进入项目目录，使用tree命令查看文件结构，当前目录包括四个子目录，分别包含abacus部分的工作、cp2k部分的工作、elk部分的工作以及summary总结工作。其中需要准备的输入文件是elk的elk.in(类似INPUT文件)、Si.in(描述元素性质的文件)，cp2k的input.inp文件，abacus的INPUT文件、STRU文件。
代码
文本
[1]
# 进入工作目录，使用tree命令查看文件结构
!cd /root/abacus_benchmark_demo && tree -L 2
. ├── abacus │ ├── abacus_eos │ ├── abacus_kconv │ ├── abacus_orb_gen │ ├── abacus_pw_ecutconv │ ├── abacus_pw_eos │ └── abacus_upf_gen ├── cp2k │ ├── cp2k_cutoffconv │ ├── cp2k_eos │ ├── cp2k_kconv │ └── cp2k_relcutoffconv ├── elk │ ├── elk_eos │ ├── elk_kconv │ └── elk_rgkmax_conv └── summary ├── summary.csv ├── summary.py ├── summary_eos.png └── summary_table.png 17 directories, 4 files
代码
文本
3. 使用ELK进行Si的EOS的全电子计算
代码
文本
我们需要确定精度的标定标准。对于使用赝势的DFT软件，最好的选择是使用全电子计算软件的计算结果作为参考。这里我们使用ELK进行Si的EOS计算（另外的例子也可见ABINIT：https://docs.abinit.org/tutorial/paw3/ ），首先对ELK进行k点收敛性测试。
代码
文本
3.1 ELK的K点收敛性测试
代码
文本
进入工作目录，查看文件
代码
文本
[2]
!cd /root/abacus_benchmark_demo/elk/elk_kconv && ls
Si.in elk_kconv.png elk_si_k10 elk_si_k4 elk_si_k7 elk.in elk_kconv.py elk_si_k2 elk_si_k5 elk_si_k8 elk_kconv.csv elk_si_k1 elk_si_k3 elk_si_k6 elk_si_k9
代码
文本
查看输入文件Si.in模版，以下任务需要修改其中的ngridk字段
代码
文本
[3]
!cat /root/abacus_benchmark_demo/elk/elk_kconv/elk.in
tasks
0
xctype
20
sppath
'./'
avec
1.8897260000 0.0000000000 0.0000000000
0.0000000000 1.8897260000 0.0000000000
0.0000000000 0.0000000000 1.8897260000
scale
5.43070000
atoms
1 : nspecies
'Si.in' : spfname
8
0.0000000000 0.0000000000 0.0000000000 0.000 0.000 0.000
0.0000000000 0.5000000000 0.5000000000 0.000 0.000 0.000
0.5000000000 0.0000000000 0.5000000000 0.000 0.000 0.000
0.5000000000 0.5000000000 0.0000000000 0.000 0.000 0.000
0.7500000000 0.7500000000 0.2500000000 0.000 0.000 0.000
0.7500000000 0.2500000000 0.7500000000 0.000 0.000 0.000
0.2500000000 0.7500000000 0.7500000000 0.000 0.000 0.000
0.2500000000 0.2500000000 0.2500000000 0.000 0.000 0.000
! Tolerance on conv. of total energy (absolute)
epsengy
1.e-6
! Use adaptive linear mixing of densities
! 1 - Adaptive linear
! 2 - Pulay mixing
mixtype
3
stype
3
tempk
300
ngridk
1 1 1代码
文本
查看自动修改k点创建子目录、运行elk计算、回收elk结果并可视化的python脚本
代码
文本
[4]
!cat /root/abacus_benchmark_demo/elk/elk_kconv/elk_kconv.py
import os
import shutil
import pandas as pd
import matplotlib.pyplot as plt
# 获取当前目录
current_dir = os.getcwd()
# 要复制的文件列表
files_to_copy = ['elk.in', 'Si.in']
# 初始化列表来存储数据
data = []
# 创建子目录并复制文件，修改elk.in文件中的ngridk值
for i in range(10):
kpoints_value = 1 + i # 步长调整为1
sub_dir_name = f'elk_si_k{kpoints_value}'
sub_dir_path = os.path.join(current_dir, sub_dir_name)
# 创建子目录
os.makedirs(sub_dir_path, exist_ok=True)
# 复制文件到子目录
for file_name in files_to_copy:
shutil.copy(os.path.join(current_dir, file_name), sub_dir_path)
# 修改elk.in文件中的ngridk下一行的值
elk_in_path = os.path.join(sub_dir_path, 'elk.in')
with open(elk_in_path, 'r+') as file:
lines = file.readlines()
# 寻找ngridk关键字并修改下一行的值
for j in range(len(lines)):
if 'ngridk' in lines[j].lower():
lines[j + 1] = f'{kpoints_value} {kpoints_value} {kpoints_value}\n'
break
file.seek(0) # 回到文件开头
file.writelines(lines) # 写入修改后的内容
file.truncate() # 删除未写入的部分
# 可取消注释运行，大约需要几分钟时间
# # 切换到子目录并运行计算任务
os.chdir(sub_dir_path)
os.system('mpirun -np 16 /root/elk-9.6.8/src/elk | tee out.log')
os.chdir(current_dir) # 返回到上一级目录
# 读取计算结果
totenergy_file = os.path.join(sub_dir_path, 'TOTENERGY.OUT')
# 检查TOTENERGY.OUT文件是否存在
if os.path.isfile(totenergy_file):
# 读取TOTENERGY.OUT文件的最后一行
with open(totenergy_file, 'r') as f:
lines = f.readlines()
if lines:
try:
total_energy = float(lines[-1].strip())
# 将k点和总能存储到列表中
data.append([kpoints_value, total_energy])
except ValueError:
print(f"Error reading total energy in {totenergy_file}")
# 将数据转换为DataFrame并保存为CSV文件
df = pd.DataFrame(data, columns=['K Points', 'Total Energy'])
df.to_csv('elk_kconv.csv', index=False)
# 绘制图表
plt.figure(figsize=(10, 6)) # 设置图形大小
plt.plot(df['K Points'], df['Total Energy'], marker='o') # 绘制折线图
plt.title('Kpoints vs Total Energy') # 设置图表标题
plt.xlabel('Kpoints') # 设置x轴标签
plt.ylabel('Total Energy (Hartree)') # 设置y轴标签
plt.grid(True) # 显示网格
plt.savefig('elk_kconv.png') # 保存图形
plt.show() # 显示图形
print("所有子目录创建、文件修改、计算任务和数据处理完成。")
代码
文本
运行这行代码会自动进行具有不同k点的计算
代码
文本
[5]
# 取消注释运行
# !cd /root/abacus_benchmark_demo/elk/elk_kconv && python /root/abacus_benchmark_demo/elk/elk_kconv/elk_kconv.py
代码
文本
显示k点收敛曲线
代码
文本
[6]
from PIL import Image
img=Image.open('/root/abacus_benchmark_demo/elk/elk_kconv/elk_kconv.png')
img
代码
文本
根据收敛曲线以及下文中abacus和cp2k的k点收敛结果，选定k点为（7，7，7）
代码
文本
3.2 ELK的rgkmax收敛性测试
代码
文本
rgkmax是elk所采用的FP-LAPW方法中muffin-tin potential实空间的半径（记作rmt）与平面波G+k（记作gk）的乘积。因此测试rgkmax的收敛等同于测试平面波基组的数量。进入工作目录，查看文件
代码
文本
[7]
!cd /root/abacus_benchmark_demo/elk/elk_rgkmax_conv && ls
Si.in elk_rgkmax_conv.png elk_si_rgk3 elk_si_rgk6 elk_si_rgk9 elk.in elk_rgkmax_conv.py elk_si_rgk4 elk_si_rgk7 elk_rgkmax_conv.csv elk_si_rgk10 elk_si_rgk5 elk_si_rgk8
代码
文本
取消注释可运行elk的rgkmax收敛性测试计算和结果回收脚本
代码
文本
[8]
# !cd /root/abacus_benchmark_demo/elk/elk_rgkmax_conv && python /root/abacus_benchmark_demo/elk/elk_rgkmax_conv/elk_rgkmax_conv.py
代码
文本
[9]
!cat /root/abacus_benchmark_demo/elk/elk_rgkmax_conv/elk_rgkmax_conv.py
import os
import shutil
import pandas as pd
import matplotlib.pyplot as plt
# 获取当前目录
current_dir = os.getcwd()
# 要复制的文件列表
files_to_copy = ['elk.in', 'Si.in']
# 初始化列表来存储数据
data = []
# 创建子目录并复制文件，修改elk.in文件中的rgkmax值
for i in range(8):
rgkmax_value = 3 + i # 步长调整为1
sub_dir_name = f'elk_si_rgk{rgkmax_value}'
sub_dir_path = os.path.join(current_dir, sub_dir_name)
# 创建子目录
os.makedirs(sub_dir_path, exist_ok=True)
# 复制文件到子目录
for file_name in files_to_copy:
shutil.copy(os.path.join(current_dir, file_name), sub_dir_path)
# 修改elk.in文件中的rgkmax值
elk_in_path = os.path.join(sub_dir_path, 'elk.in')
with open(elk_in_path, 'r+') as file:
lines = file.readlines()
# 寻找rgkmax关键字并修改下一行的值
for j in range(len(lines)):
if 'rgkmax' in lines[j].lower():
lines[j + 1] = f'{rgkmax_value}\n'
break
file.seek(0) # 回到文件开头
file.writelines(lines) # 写入修改后的内容
file.truncate() # 删除未写入的部分
# 可取消注释运行，大约需要几分钟时间
# 切换到子目录并运行计算任务
os.chdir(sub_dir_path)
os.system('mpirun -np 16 /root/elk-9.6.8/src/elk | tee out.log')
os.chdir(current_dir) # 返回到上一级目录
# 读取计算结果
totenergy_file = os.path.join(sub_dir_path, 'TOTENERGY.OUT')
# 检查TOTENERGY.OUT文件是否存在
if os.path.isfile(totenergy_file):
# 读取TOTENERGY.OUT文件的最后一行
with open(totenergy_file, 'r') as f:
lines = f.readlines()
if lines:
try:
total_energy = float(lines[-1].strip())
# 将rgkmax和总能存储到列表中
data.append([rgkmax_value, total_energy])
except ValueError:
print(f"Error reading total energy in {totenergy_file}")
# 将数据转换为DataFrame并保存为CSV文件
df = pd.DataFrame(data, columns=['RGKMAX', 'Total Energy'])
df.to_csv('elk_rgkmax_conv.csv', index=False)
# 绘制图表
plt.figure(figsize=(10, 6)) # 设置图形大小
plt.plot(df['RGKMAX'], df['Total Energy'], marker='o') # 绘制折线图
plt.title('RGKMAX vs Total Energy') # 设置图表标题
plt.xlabel('RGKMAX') # 设置x轴标签
plt.ylabel('Total Energy (Hartree)') # 设置y轴标签
plt.grid(True) # 显示网格
plt.savefig('elk_rgkmax_conv.png') # 保存图形
plt.show() # 显示图形
print("所有子目录创建、文件修改、计算任务和数据处理完成。")
代码
文本
查看收敛性结果
代码
文本
[10]
from PIL import Image
img=Image.open('/root/abacus_benchmark_demo/elk/elk_rgkmax_conv/elk_rgkmax_conv.png')
img
代码
文本
3.3 ELK的Si的EOS计算
代码
文本
进入工作目录，查看文件
代码
文本
[11]
!cd /root/abacus_benchmark_demo/elk/elk_eos && ls
Si.in Si_eos_5.40 Si_eos_5.70 elk.in elk_eos.py Si_eos_5.20 Si_eos_5.50 Si_eos_5.80 elk_eos.csv elk_time.csv Si_eos_5.30 Si_eos_5.60 Si_eos_5.90 elk_eos.png
代码
文本
查看自动化eos计算脚本
代码
文本
[12]
!cat /root/abacus_benchmark_demo/elk/elk_eos/elk_eos.py
import os
import shutil
import pandas as pd
import matplotlib.pyplot as plt
import re
import csv
# 定义提取总时间的函数
def extract_total_time(file_path):
with open(file_path, 'r') as file:
content = file.read()
# 定义正则表达式来匹配 Timings (CPU seconds) 部分的 total 值
match = re.search(r'total\s*:\s*([\d\.]+)', content)
if match:
return float(match.group(1)) # 返回 total 对应的数值
else:
return None
# 保存INFO.OUT提取的数据到CSV文件
def save_to_csv(data, output_file):
with open(output_file, 'w', newline='') as csvfile:
writer = csv.writer(csvfile)
writer.writerow(['Directory', 'Total Time (CPU seconds)'])
writer.writerows(data)
# 计算平均值并添加到最后一行
if data:
avg_time = sum(row[1] for row in data) / len(data)
writer.writerow(['Average', avg_time])
# 获取当前目录
current_dir = os.getcwd()
# 要复制的文件列表
files_to_copy = ['elk.in', 'Si.in']
# 初始化列表来存储数据
data = []
timing_data = []
# 创建子目录并复制文件，修改elk.in文件中的scale值
for i in range(8):
scale_value = 5.2 + i * 0.1 # 步长为0.1
sub_dir_name = f'Si_eos_{scale_value:.2f}'
sub_dir_path = os.path.join(current_dir, sub_dir_name)
# 创建子目录
os.makedirs(sub_dir_path, exist_ok=True)
# 复制文件到子目录
for file_name in files_to_copy:
shutil.copy(os.path.join(current_dir, file_name), sub_dir_path)
# 修改elk.in文件中的scale下一行的值
elk_in_path = os.path.join(sub_dir_path, 'elk.in')
with open(elk_in_path, 'r+') as file:
lines = file.readlines()
# 寻找scale关键字并修改下一行的值
for j in range(len(lines)):
if 'scale' in lines[j].lower():
lines[j + 1] = f'{scale_value:.2f}\n' # 保留两位小数
break
file.seek(0) # 回到文件开头
file.writelines(lines) # 写入修改后的内容
file.truncate() # 删除未写入的部分
# 可取消注释运行，大约需要几分钟时间
# 切换到子目录并运行计算任务
os.chdir(sub_dir_path)
os.system('mpirun -np 16 /root/elk-9.6.8/src/elk | tee out.log')
os.chdir(current_dir) # 返回到上一级目录
# 读取计算结果
totenergy_file = os.path.join(sub_dir_path, 'TOTENERGY.OUT')
info_out_file = os.path.join(sub_dir_path, 'INFO.OUT')
# 检查TOTENERGY.OUT文件是否存在
if os.path.isfile(totenergy_file):
# 读取TOTENERGY.OUT文件的最后一行
with open(totenergy_file, 'r') as f:
lines = f.readlines()
if lines:
try:
total_energy = float(lines[-1].strip())
# 将scale值和总能存储到列表中
data.append([scale_value, total_energy])
except ValueError:
print(f"Error reading total energy in {totenergy_file}")
# 检查INFO.OUT文件是否存在并提取总时间
if os.path.isfile(info_out_file):
total_time = extract_total_time(info_out_file)
if total_time is not None:
timing_data.append([sub_dir_name, total_time])
# 处理数据，将第一列转化为体积，第二列转化为eV单位
processed_data = [[scale**3, energy*27.2107219] for scale, energy in data]
# 将数据转换为DataFrame并保存为CSV文件
df = pd.DataFrame(processed_data, columns=['Volume (Å^3)', 'Total Energy (eV)'])
df.to_csv('elk_eos.csv', index=False)
# 绘制图表
plt.figure(figsize=(10, 6)) # 设置图形大小
plt.plot(df['Volume (Å^3)'], df['Total Energy (eV)'], marker='o') # 绘制折线图
plt.title('Volume vs Total Energy') # 设置图表标题
plt.xlabel('Volume (Å^3)') # 设置x轴标签
plt.ylabel('Total Energy (eV)') # 设置y轴标签
plt.grid(True) # 显示网格
plt.savefig('elk_eos.png') # 保存图形
plt.show() # 显示图形
# 保存提取的总时间数据到CSV
save_to_csv(timing_data, 'elk_time.csv')
print("所有子目录创建、文件修改、计算任务和数据处理完成。")
print(f'计算时间数据已保存到 elk_time.csv')
代码
文本
取消注释可以运行自动化eos计算脚本，进行eos计算和结果回收（需要十分钟左右）
代码
文本
[56]
# !cd /root/abacus_benchmark_demo/elk/elk_eos && python /root/abacus_benchmark_demo/elk/elk_eos/elk_eos.py
代码
文本
可视化eos计算结果
代码
文本
[14]
from PIL import Image
img=Image.open('/root/abacus_benchmark_demo/elk/elk_eos/elk_eos.png')
img
代码
文本
4. 使用CP2K基于GTH-PBE赝势和DZVP-GTH-q4基组进行Si的EOS计算
代码
文本
选定cp2k常用的GTH-PBE赝势和DZVP-GTH-q4基组进行Si的EOS计算，在EOS计算前需要进行k点、CUTOFF、REL_CUTOFF收敛性测试。
代码
文本
4.1 CP2K的K点收敛性测试
代码
文本
进入工作目录，查看文件
代码
文本
[15]
! cd /root/abacus_benchmark_demo/cp2k/cp2k_kconv && ls
cp2k_Si_kpoints_1 cp2k_Si_kpoints_4 cp2k_Si_kpoints_8 cp2k_kconv.py cp2k_Si_kpoints_10 cp2k_Si_kpoints_5 cp2k_Si_kpoints_9 input.inp cp2k_Si_kpoints_2 cp2k_Si_kpoints_6 cp2k_kconv.csv cp2k_Si_kpoints_3 cp2k_Si_kpoints_7 cp2k_kconv.png
代码
文本
查看输入文件模版input.inp
代码
文本
[16]
!cat /root/abacus_benchmark_demo/cp2k/cp2k_kconv/input.inp
&GLOBAL
PROJECT Si
PRINT_LEVEL MEDIUM
RUN_TYPE ENERGY
&END GLOBAL
&FORCE_EVAL
METHOD Quickstep
&SUBSYS
&CELL
A [angstrom] 5.43070000 0.00000000 0.00000000
B [angstrom] 0.00000000 5.43070000 0.00000000
C [angstrom] 0.00000000 0.00000000 5.43070000
PERIODIC XYZ
&END CELL
&COORD
SCALED
Si 0.00000000 0.00000000 0.00000000
Si 0.00000000 0.50000000 0.50000000
Si 0.50000000 0.00000000 0.50000000
Si 0.50000000 0.50000000 0.00000000
Si 0.75000000 0.75000000 0.25000000
Si 0.75000000 0.25000000 0.75000000
Si 0.25000000 0.75000000 0.75000000
Si 0.25000000 0.25000000 0.25000000
&END COORD
&KIND Si
ELEMENT Si
BASIS_SET DZVP-GTH-q4
POTENTIAL GTH-PBE
&END KIND
&END SUBSYS
&DFT
BASIS_SET_FILE_NAME GTH_BASIS_SETS
POTENTIAL_FILE_NAME POTENTIAL
# WFN_RESTART_FILE_NAME Si2-RESTART.wfn
CHARGE 0 #Net charge
MULTIPLICITY 1 #Spin multiplicity
&QS
EPS_DEFAULT 1.0E-11 #Set all EPS_xxx to values such that the energy will be correct up to this value
&END QS
&POISSON
PERIODIC XYZ #Direction(s) of PBC for calculating electrostatics
PSOLVER PERIODIC #The way to solve Poisson equation
&END POISSON
&XC
&XC_FUNCTIONAL PBE
&END XC_FUNCTIONAL
&END XC
&MGRID
CUTOFF 500
REL_CUTOFF 60
&END MGRID
&KPOINTS
SCHEME MONKHORST-PACK 1 1 1
&END KPOINTS
&SCF
MAX_SCF 128
EPS_SCF 1E-06 #Convergence threshold of density matrix of inner SCF
# SCF_GUESS RESTART #Use wavefunction from WFN_RESTART_FILE_NAME file as initial guess
# IGNORE_CONVERGENCE_FAILURE #Continue calculation even if SCF not converged, works for version >= 2024.1
&DIAGONALIZATION
ALGORITHM STANDARD #Algorithm for diagonalization
&END DIAGONALIZATION
&MIXING #How to mix old and new density matrices
METHOD BROYDEN_MIXING #PULAY_MIXING is also a good alternative
ALPHA 0.4 #Default. Mixing 40% of new density matrix with the old one
NBROYDEN 8 #Default is 4. Number of previous steps stored for the actual mixing scheme
&END MIXING
&PRINT
&RESTART #Note: Use "&RESTART OFF" can prevent generating .wfn file
BACKUP_COPIES 0 #Maximum number of backup copies of wfn file. 0 means never
&END RESTART
&END PRINT
&END SCF
&END DFT
&END FORCE_EVAL
代码
文本
查看cp2k的k点收敛性测试脚本
代码
文本
[17]
!cat /root/abacus_benchmark_demo/cp2k/cp2k_kconv/cp2k_kconv.py
import os
import shutil
import re
import csv
import pandas as pd
import matplotlib.pyplot as plt
# 获取当前目录
current_dir = os.getcwd()
# 要复制的文件列表
files_to_copy = ['input.inp']
# 初始化列表来存储数据
data = []
# 创建子目录并复制文件，修改input.inp文件中的KPOINTS字段的值
for i in range(10):
kpoints_value = i + 1 # KPOINTS的递增值
sub_dir_name = f'cp2k_Si_kpoints_{kpoints_value}' # 根据kpoints_value命名子目录
sub_dir_path = os.path.join(current_dir, sub_dir_name)
# 创建子目录
os.makedirs(sub_dir_path, exist_ok=True)
# 复制文件到子目录
for file_name in files_to_copy:
shutil.copy(os.path.join(current_dir, file_name), sub_dir_path)
# 修改input.inp文件中的KPOINTS字段的值
input_inp_path = os.path.join(sub_dir_path, 'input.inp')
with open(input_inp_path, 'r+') as file:
content = file.read()
# 修改KPOINTS字段的值
lines = content.splitlines()
for j, line in enumerate(lines):
if line.strip() == '&KPOINTS':
# 修改下一行的SCHEME MONKHORST-PACK的值
lines[j + 1] = f' SCHEME MONKHORST-PACK {kpoints_value} {kpoints_value} {kpoints_value}'
break
# 将修改后的内容写回文件
content = '\n'.join(lines)
file.seek(0) # 回到文件开头
file.write(content) # 写入修改后的内容
file.truncate() # 删除未写入的部分
# 去掉以下三行代码的注释可以进行计算，所需时长约几分钟
# 切换到子目录并运行计算任务
os.chdir(sub_dir_path)
os.system('mpirun -n 16 cp2k.popt -i input.inp -o output.log')
os.chdir(current_dir) # 返回到上一级目录
# 遍历所有子目录，提取input.inp和output.log中的数据
for root, dirs, files in os.walk(current_dir):
for dir_name in dirs:
dir_path = os.path.join(root, dir_name)
# 构建 input.inp 和 output.log 文件路径
input_file_path = os.path.join(dir_path, 'input.inp')
log_file_path = os.path.join(dir_path, 'output.log')
if os.path.isfile(input_file_path) and os.path.isfile(log_file_path):
try:
# 从 input.inp 文件中读取 SCHEME MONKHORST-PACK 后的第一个数字
with open(input_file_path, 'r') as f:
for line in f:
if 'SCHEME MONKHORST-PACK' in line:
match = re.search(r'SCHEME MONKHORST-PACK\s+(\d+)', line)
if match:
scheme_number = int(match.group(1))
break
# 从 output.log 文件中读取 Total energy: 后的数字
with open(log_file_path, 'r') as f:
for line in f:
if 'Total energy:' in line:
total_energy = float(line.split()[2])
break
# 将读取的数据添加到列表中
data.append([scheme_number, total_energy])
except ValueError:
# 跳过无法转换为浮点数的目录名
print(f"无法从目录 {dir_name} 中提取数据")
continue
# 将数据按照Scheme Number排序
data.sort(key=lambda x: x[0])
# 将数据写入 CSV 文件
with open('cp2k_kconv.csv', 'w', newline='') as csvfile:
writer = csv.writer(csvfile)
writer.writerow(['Scheme Number', 'Total Energy'])
writer.writerows(data)
print("数据已保存到 cp2k_kconv.csv 文件中。")
# 读取CSV文件并绘制图表
df = pd.read_csv('cp2k_kconv.csv')
plt.figure(figsize=(10, 6)) # 设置图形大小
plt.plot(df['Scheme Number'], df['Total Energy'], marker='o') # 绘制折线图
plt.title('Kpoints vs Total Energy') # 设置图表标题
plt.xlabel('Kpoints') # 设置x轴标签
plt.ylabel('Total Energy (Hartree)') # 设置y轴标签
plt.grid(True) # 显示网格
plt.savefig('cp2k_kconv.png') # 保存图形
plt.show() # 显示图形
print("所有子目录创建、文件修改、计算任务和数据处理完成。")
代码
文本
取消注释即可运行cp2k的计算并回收结果
代码
文本
[18]
# !cd /root/abacus_benchmark_demo/cp2k/cp2k_kconv && python /root/abacus_benchmark_demo/cp2k/cp2k_kconv/cp2k_kconv.py
代码
文本
可视化计算结果
代码
文本
[19]
from PIL import Image
img=Image.open('/root/abacus_benchmark_demo/cp2k/cp2k_kconv/cp2k_kconv.png')
img
代码
文本
根据k点收敛曲线，选定k点为（7，7，7）
代码
文本
4.2 CP2K的CUTOFF收敛性测试
代码
文本
进入工作目录，查看文件
代码
文本
[20]
! cd /root/abacus_benchmark_demo/cp2k/cp2k_cutoffconv && ls
cp2k_cutoff_100 cp2k_cutoff_350 cp2k_cutoff_600 cp2k_cutoffconv.png cp2k_cutoff_150 cp2k_cutoff_400 cp2k_cutoff_650 cp2k_cutoffconv.py cp2k_cutoff_200 cp2k_cutoff_450 cp2k_cutoff_700 input.inp cp2k_cutoff_250 cp2k_cutoff_500 cp2k_cutoff_750 cp2k_cutoff_300 cp2k_cutoff_550 cp2k_cutoffconv.csv
代码
文本
取消注释可运行计算并回收结果
代码
文本
[21]
# !cd /root/abacus_benchmark_demo/cp2k/cp2k_cutoffconv && python /root/abacus_benchmark_demo/cp2k/cp2k_cutoffconv/cp2k_cutoffconv.py
代码
文本
可视化计算结果
代码
文本
[22]
from PIL import Image
img=Image.open('/root/abacus_benchmark_demo/cp2k/cp2k_cutoffconv/cp2k_cutoffconv.png')
img
代码
文本
根据收敛性曲线，选定CUTOFF为300 Hartree。
代码
文本
4.3 CP2K的REL_CUTOFF收敛性测试
代码
文本
进入工作目录，查看文件
代码
文本
[23]
! cd /root/abacus_benchmark_demo/cp2k/cp2k_relcutoffconv && ls
cp2k_relcutconv.png cp2k_relcutoff_40 cp2k_relcutoff_70 input.inp cp2k_relcutoff_20 cp2k_relcutoff_50 cp2k_relcutoffconv.csv cp2k_relcutoff_30 cp2k_relcutoff_60 cp2k_relcutoffconv.py
代码
文本
取消注释可运行计算并回收结果
代码
文本
[24]
# !cd /root/abacus_benchmark_demo/cp2k/cp2k_relcutoffconv && python /root/abacus_benchmark_demo/cp2k/cp2k_relcutoffconv/cp2k_relcutoffconv.py
代码
文本
可视化结果
代码
文本
[25]
from PIL import Image
img=Image.open('/root/abacus_benchmark_demo/cp2k/cp2k_relcutoffconv/cp2k_relcutconv.png')
img
代码
文本
根据收敛曲线，选定收敛REL_CUTOFF为50 Hartree。
代码
文本
4.4 CP2K的Si的EOS计算
代码
文本
进入工作目录，查看文件
代码
文本
[26]
! cd /root/abacus_benchmark_demo/cp2k/cp2k_eos && ls
Si_eos_5.20 Si_eos_5.50 Si_eos_5.80 cp2k_eos.png input.inp Si_eos_5.30 Si_eos_5.60 Si_eos_5.90 cp2k_eos.py Si_eos_5.40 Si_eos_5.70 cp2k_eos.csv cp2k_time.csv
代码
文本
取消注释可运行计算并回收结果
代码
文本
[27]
# !python /root/abacus_benchmark_demo/cp2k/cp2k_eos/cp2k_eos.py
代码
文本
可视化EOS计算结果
代码
文本
[28]
from PIL import Image
img=Image.open('/root/abacus_benchmark_demo/cp2k/cp2k_eos/cp2k_eos.png')
img
代码
文本
5. 使用ABACUS基于由CP2K的GTH-PBE赝势生成的赝势进行Si的EOS计算
代码
文本
为了对齐赝势，可以将cp2k使用的GTH-PBE赝势通过ATOM模块转换成abacus可读取的upf赝势文件进行直接比较，并且可以基于此生成的赝势进一步生成轨道文件，进行lcao计算。
代码
文本
5.1 基于CP2K的GTH-PBE赝势生成ABACUS可读取的upf赝势文件
代码
文本
进入工作目录，查看文件
代码
文本
[29]
!cd /root/abacus_benchmark_demo/abacus/abacus_upf_gen && ls
PROJECT-Si.UPF-1.upf Si.inp Si.out
代码
文本
查看生成upf赝势任务的输入文件
代码
文本
[30]
!cat /root/abacus_benchmark_demo/abacus/abacus_upf_gen/Si.inp
&GLOBAL
PROGRAM_NAME ATOM
&END GLOBAL
&ATOM
ELEMENT Si
ELECTRON_CONFIGURATION CORE 3s2 3p2
CORE [Ne]
&METHOD
METHOD_TYPE KOHN-SHAM
&XC
&XC_FUNCTIONAL PBE
&END XC_FUNCTIONAL
&END XC
&END METHOD
&PP_BASIS
BASIS_TYPE GEOMETRICAL_GTO
&END PP_BASIS
&POTENTIAL
PSEUDO_TYPE GTH
>H_POTENTIAL
2 2
0.44000000 1 -6.26928833
2
0.43563383 2 8.95174150 -2.70627082
3.49378060
0.49794218 1 2.43127673
&END
&END POTENTIAL
&PRINT
&ANALYZE_BASIS
OVERLAP_CONDITION_NUMBER T
COMPLETENESS T
&END ANALYZE_BASIS
&UPF_FILE
FILENAME Si.UPF
&END
&END
&END ATOM代码
文本
取消注释可运行代码并得到对应的赝势upf文件
代码
文本
[31]
#!cd /root/abacus_benchmark_demo/abacus/abacus_upf_gen && mpirun -np 8 /root/cp2k/exe/local/cp2k.popt -i Si.inp -o Si.out
代码
文本
查看生成的赝势的前15行
代码
文本
[32]
!head -15 /root/abacus_benchmark_demo/abacus/abacus_upf_gen/PROJECT-Si.UPF-1.upf
<UPF version="2.0.1">
<PP_INFO>
Converted from CP2K GTH format
<PP_INPUTFILE>
Si
2 2 0 0
0.44000000000000 1 -6.26928833000000
2
0.43563383000000 2 8.95174150000000 -2.70627082000000
3.49378060000000
0.49794218000000 1 2.43127673000000
</PP_INPUTFILE>
</PP_INFO>
<PP_HEADER
generated="Generated in analytical, separable form"
代码
文本
5.2 基于生成的赝势文件进行平面波计算确定生成轨道的截断能
代码
文本
进入工作目录，查看文件
代码
文本
[33]
!cd /root/abacus_benchmark_demo/abacus/abacus_pw_ecutconv && ls
INPUT abacus_si_ecutwfc130.00 abacus_si_ecutwfc210.00 PROJECT-Si.UPF-1.upf abacus_si_ecutwfc140.00 abacus_si_ecutwfc30.00 STRU abacus_si_ecutwfc150.00 abacus_si_ecutwfc40.00 abacus_pw_ecutconv.csv abacus_si_ecutwfc160.00 abacus_si_ecutwfc50.00 abacus_pw_ecutconv.png abacus_si_ecutwfc170.00 abacus_si_ecutwfc60.00 abacus_pw_ecutconv.py abacus_si_ecutwfc180.00 abacus_si_ecutwfc70.00 abacus_si_ecutwfc100.00 abacus_si_ecutwfc190.00 abacus_si_ecutwfc80.00 abacus_si_ecutwfc110.00 abacus_si_ecutwfc20.00 abacus_si_ecutwfc90.00 abacus_si_ecutwfc120.00 abacus_si_ecutwfc200.00 cp2k_gen_gth_pbe_si.upf
代码
文本
取消注释运行abacus的pw计算并回收计算结果
代码
文本
[3]
!cd /root/abacus_benchmark_demo/abacus/abacus_pw_ecutconv && python /root/abacus_benchmark_demo/abacus/abacus_pw_ecutconv/abacus_pw_ecutconv.py
^C
代码
文本
[35]
from PIL import Image
img=Image.open('/root/abacus_benchmark_demo/abacus/abacus_pw_ecutconv/abacus_pw_ecutconv.png')
img
代码
文本
根据能量收敛曲线选定Ecutwfc为100 Ry。
代码
文本
5.3 基于生成的upf赝势文件生成相应的orb轨道文件
代码
文本
由于cp2k生成的赝势文件的
代码
文本
[36]
!cd /root/abacus_benchmark_demo/abacus/abacus_orb_gen && tree -L 1
. ├── PROJECT-Si.UPF-1.upf ├── SIAB ├── SIAB_INPUT.json ├── Si-dimer-1.75 ├── Si-dimer-2.00 ├── Si-dimer-2.25 ├── Si-dimer-2.75 ├── Si-dimer-3.75 ├── Si-monomer ├── Si-trimer-1.90 ├── Si-trimer-2.10 ├── Si-trimer-2.60 ├── Si_1s1p ├── Si_2s2p1d ├── Si_3s3p2d ├── Si_ONCV_PBE-1.0.upf ├── cp2k_gen_gth_pbe_si.upf └── pp_info_modify.py 13 directories, 5 files
代码
文本
[37]
!head -15 /root/abacus_benchmark_demo/abacus/abacus_orb_gen/PROJECT-Si.UPF-1.upf
<UPF version="2.0.1">
<PP_INFO>
Converted from CP2K GTH format
<PP_INPUTFILE>
Si
2 2 0 0
0.44000000000000 1 -6.26928833000000
2
0.43563383000000 2 8.95174150000000 -2.70627082000000
3.49378060000000
0.49794218000000 1 2.43127673000000
</PP_INPUTFILE>
</PP_INFO>
<PP_HEADER
generated="Generated in analytical, separable form"
代码
文本
运行以下脚本可以自动将原upf文件的更新，并保存为cp2k_gen_gth_pbe_si.upf
代码
文本
[38]
!cd /root/abacus_benchmark_demo/abacus/abacus_orb_gen && python /root/abacus_benchmark_demo/abacus/abacus_orb_gen/pp_info_modify.py
已将Si_ONCV_PBE-1.0.upf中的<PP_INFO>部分替换到PROJECT-Si.UPF-1.upf中，并保存到cp2k_gen_gth_pbe_si.upf。
代码
文本
查看修改后的upf文件的前15行
代码
文本
[39]
!head -15 /root/abacus_benchmark_demo/abacus/abacus_orb_gen/cp2k_gen_gth_pbe_si.upf
<UPF version="2.0.1"> <PP_INFO> This pseudopotential file has been produced using the code ONCVPSP (Optimized Norm-Conservinng Vanderbilt PSeudopotential) scalar-relativistic version 2.1.1, 03/26/2014 by D. R. Hamann The code is available through a link at URL www.mat-simresearch.com. Documentation with the package provides a full discription of the input data below. While it is not required under the terms of the GNU GPL, it is suggested that you cite D. R. Hamann, Phys. Rev. B 88, 085117 (2013) in any publication using these pseudopotentials.
代码
文本
查看使用upf文件生成lcao轨道文件的配置文件
代码
文本
[40]
!cat /root/abacus_benchmark_demo/abacus/abacus_orb_gen/SIAB_INPUT.json
{
"pseudo_dir": "./",
"ecutwfc": 100.0,
"pseudo_name": "cp2k_gen_gth_pbe_si.upf",
"bessel_nao_rcut": [
6,
7,
8,
9,
10
],
"smearing_sigma": 0.01,
"optimizer": "bfgs",
"max_steps": 5000,
"spill_coefs": [
0.0,
1.0
],
"spill_guess": "atomic",
"nthreads_rcut": 4,
"jY_type": "reduced",
"reference_systems": [
{
"shape": "dimer",
"nbands": "auto",
"bond_lengths": "auto",
"nspin": 1
},
{
"shape": "trimer",
"nbands": "auto",
"bond_lengths": "auto",
"nspin": 1
}
],
"orbitals": [
{
"zeta_notation": "Z",
"shape": "dimer",
"orb_ref": "none",
"nbands_ref": "occ"
},
{
"zeta_notation": "DZP",
"shape": "dimer",
"orb_ref": "Z",
"nbands_ref": "occ"
},
{
"zeta_notation": "TZDP",
"shape": "trimer",
"orb_ref": "DZP",
"nbands_ref": "all"
}
],
"environment": "",
"mpi_command": "mpirun -np 16",
"abacus_command": "abacus"
}代码
文本
取消注释，运行以下脚本可实现轨道的生成，c16m32 cpu节点运行大约需要1小时时间，不建议取消注释运行。 注：如遇到无法执行情况，可通过下方SSH连接方式直接运行轨道生成文件。
代码
文本
[1]
# !cd /root/abacus_benchmark_demo/abacus/abacus_orb_gen && /root/miniconda3/envs/orbgen/bin/python SIAB/SIAB_nouvelle.py -i SIAB_INPUT.json
===================================================================================================
.oooooo..o ooooo .o. oooooooooo.
d8P' `Y8 `888' .888. `888' `Y8b
Y88bo. 888 .8"888. 888 888
`"Y8888o. 888 .8' `888. 888oooo888'
`"Y88b 888 .88ooo8888. 888 `88b
oo .d8P 888 .8' `888. 888 .88P
8""88888P' o888o o88o o8888o o888bood8P'
New version of Systematically Improvable Atomic-orbital Basis (SIAB) method for generating
numerical atomic orbitals (NAOs) for Linar Combinations of Atomic Orbitals (LCAO) based electronic
structure calculations.
This version is refactored from PTG_DPSI, by ABACUS-AISI developers.
Github repo: https://github.com/abacusmodeling/ABACUS-orbitals
Tutorials: https://mcresearch.github.io/abacus-user-guide/abacus-nac1.html
https://mcresearch.github.io/abacus-user-guide/abacus-nac2.html
https://mcresearch.github.io/abacus-user-guide/abacus-nac3.html
See reference for more information.
===================================================================================================
Parsing SIAB input file SIAB_INPUT.json with version 0.1.0
Autoset: adding orbitals... will add electrons to l = 2 orbitals according to Hund's rule
Autoset: strategy 'fill up to' lmax = 2 subshell
Autoset: minimal number of electrons to fill up to l = 2 is 20.0
WARNING: since SIAB version 2.1(2024.6.3), the original functionality invoked by value "auto" is replaced by
"scan", and for dimer the "auto" now will directly use in-built dimer database if available, otherwise will
fall back to "scan". This warning will be print everytime if "auto" is used. To disable this warning, specify
directly the "bond_lengths" in any one of following ways:
1. a list of floats, e.g. [2.0, 2.5, 3.0]
2. a string "default", which will use default bond length for dimer, and scan for other shapes, for other shapes, will
fall back to "scan".
3. a string "scan", which will scan bond lengths for present shape.
DUPLICATE CHECK-1 pass: folder Si-dimer-1.75 exists
DUPLICATE CHECK-2 pass: INPUT and INPUTw exist
KEYWORD "pseudo_dir" has different values. Original: /personal/upf_prep/Si_orb_gen3, new: /root/abacus_benchmark_demo/abacus/abacus_orb_gen
Difference detected, start a new job.
Run ABACUS calculation on reference structure.
Reference structure: dimer
Bond length: 1.75
present directory: /root/abacus_benchmark_demo/abacus/abacus_orb_gen/Si-dimer-1.75
OMP_NUM_THREADS: 1
run with command: mpirun -np 16 abacus
ABACUS v3.7.1
Atomic-orbital Based Ab-initio Computation at UStc
Website: http://abacus.ustc.edu.cn/
Documentation: https://abacus.deepmodeling.com/
Repository: https://github.com/abacusmodeling/abacus-develop
https://github.com/deepmodeling/abacus-develop
Commit: f2d31f56c (Wed Jul 17 13:06:46 2024 +0800)
Fri Jul 19 16:39:52 2024
MAKE THE DIR : OUT.Si-dimer-1.75/
RUNNING WITH DEVICE : CPU / Intel(R)%Xeon(R)%Platinum
UNIFORM GRID DIM : 128 * 128 * 128
UNIFORM GRID DIM(BIG) : 128 * 128 * 128
DONE(1.02586 SEC) : SETUP UNITCELL
DONE(1.027 SEC) : INIT K-POINTS
---------------------------------------------------------
Self-consistent calculations for electrons
---------------------------------------------------------
SPIN KPOINTS PROCESSORS
1 1 16
---------------------------------------------------------
Use plane wave basis
---------------------------------------------------------
ELEMENT NATOM XC
Si 2
---------------------------------------------------------
Initial plane wave basis and FFT box
---------------------------------------------------------
DONE(1.05301 SEC) : INIT PLANEWAVE
DONE(1.0981 SEC) : LOCAL POTENTIAL
DONE(1.13074 SEC) : NON-LOCAL POTENTIAL
MEMORY FOR PSI (MB) : 1.03027
DONE(1.15437 SEC) : INIT BASIS
-------------------------------------------
SELF-CONSISTENT :
-------------------------------------------
START CHARGE : atomic
DONE(1.35787 SEC) : INIT SCF
* * * * * *
<< Start SCF iteration.
ITER ETOT/eV EDIFF/eV DRHO TIME/s
DA1 -2.03925398e+02 0.00000000e+00 4.4363e-01 0.59
DA2 -2.04672968e+02 -7.47569318e-01 6.0855e-01 0.44
DA3 -2.06380661e+02 -1.70769373e+00 2.9608e-02 0.38
DA4 -2.06342101e+02 3.85600124e-02 1.2667e-02 0.40
DA5 -2.06349589e+02 -7.48799588e-03 1.2559e-03 0.35
DA6 -2.06349462e+02 1.27467329e-04 2.7170e-04 0.39
DA7 -2.06350322e+02 -8.59784738e-04 8.0993e-05 0.41
DA8 -2.06350455e+02 -1.33038586e-04 1.5287e-06 0.37
DA9 -2.06350470e+02 -1.51174506e-05 1.0499e-06 0.69
^C
[mpiexec@bohrium-36094-1168648] Sending Ctrl-C to processes as requested
[mpiexec@bohrium-36094-1168648] Press Ctrl-C again to force abort
===================================================================================
= BAD TERMINATION OF ONE OF YOUR APPLICATION PROCESSES
= RANK 0 PID 939 RUNNING AT bohrium-36094-1168648
= KILLED BY SIGNAL: 2 (Interrupt)
===================================================================================
===================================================================================
= BAD TERMINATION OF ONE OF YOUR APPLICATION PROCESSES
= RANK 1 PID 940 RUNNING AT bohrium-36094-1168648
= KILLED BY SIGNAL: 2 (Interrupt)
===================================================================================
===================================================================================
= BAD TERMINATION OF ONE OF YOUR APPLICATION PROCESSES
= RANK 2 PID 941 RUNNING AT bohrium-36094-1168648
= KILLED BY SIGNAL: 2 (Interrupt)
===================================================================================
===================================================================================
= BAD TERMINATION OF ONE OF YOUR APPLICATION PROCESSES
= RANK 3 PID 942 RUNNING AT bohrium-36094-1168648
= KILLED BY SIGNAL: 2 (Interrupt)
===================================================================================
===================================================================================
= BAD TERMINATION OF ONE OF YOUR APPLICATION PROCESSES
= RANK 4 PID 943 RUNNING AT bohrium-36094-1168648
= KILLED BY SIGNAL: 2 (Interrupt)
===================================================================================
===================================================================================
= BAD TERMINATION OF ONE OF YOUR APPLICATION PROCESSES
= RANK 5 PID 944 RUNNING AT bohrium-36094-1168648
= KILLED BY SIGNAL: 2 (Interrupt)
===================================================================================
代码
文本
[ ]
# !python /root/abacus_benchmark_demo/abacus/abacus_orb_gen/SIAB/SIAB_nouvelle.py -i /root/abacus_benchmark_demo/abacus/abacus_orb_gen/SIAB_INPUT.json
代码
文本
SSH连接方式运行轨道生成代码：点击上方绿色的已连接按钮，点击通过SSH连接到节点，在底部出现的SSH命令行界面依次输入以下命令：
# 激活轨道生成环境
conda activate orbgen
#进入工作目录
cd /root/abacus_benchmark_demo/abacus/abacus_orb_gen
#运行生成轨道脚本
python SIAB/SIAB_nouvelle.py -i SIAB_INPUT.json
即可进行轨道的生成，使用C16M32的CPU节点约运行1h左右
代码
文本
查看输出的轨道文件的前20行，以Si_gga_7au_100.0Ry_2s2p1d.orb为例
代码
文本
[42]
! head -20 /root/abacus_benchmark_demo/abacus/abacus_orb_gen/Si_2s2p1d/7au_100.0Ry/Si_gga_7au_100.0Ry_2s2p1d.orb
---------------------------------------------------------------------------
Element Si
Energy Cutoff(Ry) 100.0
Radius Cutoff(a.u.) 7
Lmax 2
Number of Sorbital--> 2
Number of Porbital--> 2
Number of Dorbital--> 1
---------------------------------------------------------------------------
SUMMARY END
Mesh 701
dr 0.01
Type L N
0 0 0
-9.94327772020867e-03 -9.85918654232421e-03 -9.60697963972621e-03 -9.18685682932890e-03
-8.59915088530649e-03 -7.84432715863854e-03 -6.92298304561070e-03 -5.83584730639540e-03
-4.58377923515861e-03 -3.16776768343987e-03 -1.58892993885413e-03 1.51489538554255e-04
2.05212051968615e-03 4.11146854661615e-03 6.32791660654468e-03 8.69972693610516e-03
1.12250429524361e-02 1.39018913072218e-02 1.67281840597182e-02 1.97017209646062e-02
代码
文本
5.3 ABACUS的K点收敛性测试
代码
文本
进入工作目录，查看文件
代码
文本
[43]
!cd /root/abacus_benchmark_demo/abacus/abacus_kconv && ls
INPUT abacus_si_kspacing0.02 abacus_si_kspacing0.14 STRU abacus_si_kspacing0.04 abacus_si_kspacing0.16 Si_gga_7au_100.0Ry_2s2p1d.orb abacus_si_kspacing0.06 abacus_si_kspacing0.18 abacus_kconv.csv abacus_si_kspacing0.08 abacus_si_kspacing0.20 abacus_kconv.png abacus_si_kspacing0.10 abacus_si_kspacing0.22 abacus_kconv.py abacus_si_kspacing0.12 cp2k_gen_gth_pbe_si.upf
代码
文本
取消注释运行代码可进行abacus计算和结果回收
代码
文本
[44]
#!cd /root/abacus_benchmark_demo/abacus/abacus_kconv && python /root/abacus_benchmark_demo/abacus/abacus_kconv/abacus_kconv.py
代码
文本
可视化计算结果
代码
文本
[45]
from PIL import Image
img=Image.open('/root/abacus_benchmark_demo/abacus/abacus_kconv/abacus_kconv.png')
img
代码
文本
根据收敛曲线，参考elk和cp2k的设置，以下计算中kspacing选取0.10，对应k点为（7，7，7）
代码
文本
5.4 ABACUS的Si的EOS计算（lcao）
代码
文本
进入工作目录，查看文件
代码
文本
[46]
!cd /root/abacus_benchmark_demo/abacus/abacus_eos && ls
INPUT abacus_Si_eos_5.50 abacus_eos.png STRU abacus_Si_eos_5.60 abacus_eos.py Si_gga_7au_100.0Ry_2s2p1d.orb abacus_Si_eos_5.70 abacus_time.csv abacus_Si_eos_5.20 abacus_Si_eos_5.80 cp2k_gen_gth_pbe_si.upf abacus_Si_eos_5.30 abacus_Si_eos_5.90 abacus_Si_eos_5.40 abacus_eos.csv
代码
文本
取消注释可以运行abacus的eos计算并回收结果
代码
文本
[2]
!cd /root/abacus_benchmark_demo/abacus/abacus_eos && python /root/abacus_benchmark_demo/abacus/abacus_eos/abacus_eos.py
^C
代码
文本
可视化计算结果
代码
文本
[48]
from PIL import Image
img=Image.open('/root/abacus_benchmark_demo/abacus/abacus_eos/abacus_eos.png')
img
代码
文本
5.4 ABACUS的Si的EOS计算（pw）
代码
文本
进入工作目录，查看文件
代码
文本
[49]
!cd /root/abacus_benchmark_demo/abacus/abacus_pw_eos && ls
INPUT abacus_Si_eos_5.60 abacus_pw_eos.png STRU abacus_Si_eos_5.70 abacus_pw_eos.py abacus_Si_eos_5.20 abacus_Si_eos_5.80 abacus_pw_time.csv abacus_Si_eos_5.30 abacus_Si_eos_5.90 cp2k_gen_gth_pbe_si.upf abacus_Si_eos_5.40 abacus_eos.csv abacus_Si_eos_5.50 abacus_pw_eos.csv
代码
文本
取消注释可以运行abacus的pw的eos计算并回收结果
代码
文本
[50]
# !cd /root/abacus_benchmark_demo/abacus/abacus_pw_eos && python /root/abacus_benchmark_demo/abacus/abacus_pw_eos/abacus_pw_eos.py
代码
文本
可视化计算结果
代码
文本
[51]
from PIL import Image
img=Image.open('/root/abacus_benchmark_demo/abacus/abacus_pw_eos/abacus_pw_eos.png')
img
代码
文本
6. ABACUS和CP2K在Si体系EOS上的精度和效率对比总结
代码
文本
6.1 EOS的Birch-Murnaghan方程拟合
代码
文本
很多实际的应用，如测定相变、膨胀系数、弹性、热电效应、组分变化、应力等，都要求精确地测定晶格常数。理论计算平衡晶格常数多采用固体的状态方程（equation of state, EOS），该方程表示晶体体系的总能量E随体系体积V的变化，即具有E(V)的形式。晶体的状态方程对于基础科学和应用科学具有重要的意义，比如平衡体积（V0）、体弹性模量（B0）以及它的一阶导数（B0'），这些可测的物理量与晶体的状态方程直接相关。在高压下状态方程可用好几种不同的函数形式来描述，比如Murnaghan方程、Birch-Murnaghan(BM)方程和普适方程等。
三阶Birch-Murnaghan方程如下：
代码
文本
代码
文本
由该方程可以看出，当V=V0时，E=E0，这就分别是平衡晶格常数及其对应的能量。在实际操作中，我们就可以在平衡晶格常数a0附近计算得到若干E-V数据点，经曲线拟合得到体系最低能量对应的晶格常数。采用最小二乘法拟合可以得到V0, E0, B0, B0'。
代码
文本
6.2 EOS的比较指标
代码
文本
Delta指标用于比较通过两种不同的DFT计算方法(a)和(b)计算的EOS。这里，Delta = Delta(a, b)定义为：
其中，Ea(V)和Eb(V))分别是通过方法(a)和(b)得到的数据点的Birch-Murnaghan拟合，两条EOS曲线已经根据它们的最低能量进行了对齐，如前所述，积分范围覆盖了以中心体积V0为中心的6%的体积范围,使得最小体积为0.94 V0，最大体积为1.06 V0
代码
文本
6.3 精度和效率指标对比总结
代码
文本
进入工作目录,查看文件
代码
文本
[52]
!cd /root/abacus_benchmark_demo/summary && ls
summary.csv summary.py summary_eos.png summary_table.png
代码
文本
取消注释运行代码可实现此前计算EOS的拟合和对标elk的delta数值
代码
文本
[53]
#!python /root/abacus_benchmark_demo/summary/summary.py
代码
文本
查看EOS拟合结果
代码
文本
[54]
from PIL import Image
img=Image.open('/root/abacus_benchmark_demo/summary/summary_eos.png')
img
代码
文本
查看结果对比表格
代码
文本
[55]
from PIL import Image
img=Image.open('/root/abacus_benchmark_demo/summary/summary_table.png')
img
代码
文本
总结：在Si体系下，ABACUS通过使用和cp2k相同的赝势，以及由此赝势生成的轨道文件进行EOS计算。在其他参数尽可能对齐的情况下，ABACUS的平面波方法（pw）计算精度最高，但耗时最长；ABACUS的轨道方法（lcao）计算精度略低于平面波（pw）方法但优于CP2K，同时耗时约为cp2k的60%，在精度和效率上均优于cp2k。
代码
文本
已赞2
本文被以下合集收录
ML
微信用户Sj6o
更新于 2024-04-29
10 篇0 人关注
推荐阅读
公开
ABACUS对比CP2K精度和效率测试 | 基于phono3py计算Si的晶格热导率
John Yan
更新于 2024-08-07
公开
DPA-2与化学反应过渡态搜索
李博文
发布于 2024-02-24
12 赞27 转存文件4 评论