虽然实验上Cu的不均匀扩散主要出现在晶界和空洞处，但由于拟合势函数时并未拟合表面构型和晶界构型，因此计算性质时也不能研究未拟合材料的性质。如果后续想扩展研究则需要验证势函数准确性，并添加合适的训练集进行验证和后续的计算。
由于拟合DP势函数时考虑了单原子和双原子掺杂，因此研究单原子铜在铍晶体中的扩散可以进行，而晶界和表面上的扩散则需要验证，如果验证出现较大偏差则需要补充数据集再训练势函数。
准备文件
计划参考单斌老师的书籍和Bohrium上的notebook进行准备input文件，利用atomsk搭建相关的模型文件
LAMMPS做扩散的输入文件input样例--从《DFT 到 MD｜超详细「深度势能」材料计算上手指南》中复制并修改。
label loop
variable temperature index 50 250 450 650
variable data_file index Be127_Cu1_hcp.lmp
variable path_to_dp string /home/zxg/BeCu/dpgennew/dp_potential/inter09
variable path_to_conf string /home/zxg/BeCu/STRU/Be/lmp_file
shell mkdir msd_${data_file}_temp_${temperature}
shell cd msd_${data_file}_temp_${temperature}
log log_${data_file}_temp_${temperature}.lammps
units metal
boundary p p p
atom_style atomic
read_data ${path_to_conf}/${data_file}
#replicate 2 2 2
mass 1 9.012182
mass 2 63.546
group Be type 1
group Cu type 2
pair_style deepmd ${path_to_dp}/000/frozen_model.pb ${path_to_dp}/001/frozen_model.pb ${path_to_dp}/002/frozen_model.pb ${path_to_dp}/003/frozen_model.pb out_freq 100 out_file model_devi${data_file}_temp${temperature}.out
pair_coeff * *
velocity all create ${temperature} 23456789
fix 1 all nvt temp ${temperature} ${temperature} 0.5
timestep 0.001
#rdf calculation
compute rdf all rdf 100 1 1 1 2 2 2
fix 2 all ave/time 100 1 100 c_rdf[*] file ${data_file}_temp_${temperature}.rdf mode vector
#msd calculation
compute msd1 Be msd
compute msd2 Cu msd
fix 3 all ave/time 100 1 100 c_msd1[4] c_msd2[4] file ${data_file}_temp_${temperature}.msd
thermo_style custom step temp pe ke etotal press density lx ly lz vol
thermo 100
dump 1 all custom 100 ${data_file}_temp_${temperature}.dump id type x y z
run 1000000
write_data ${data_file}_temp_${temperature}.restart_data
shell cd ..
clear
next temperature
jump SELF loop
#clear
#next data_file
#jump SELF loop
DP势函数，使用DPGEN迭代过程中的势函数，在INPUT文件夹加入相对准确的路径即可。
结构文件使用DPGEN迭代过程中的单晶Be128对应的lammpsconfs，并将最后一个原子类型由代表Be的1修改成代表Cu铜的2
运行lmp程序
运行命令执行LAMMPS进行计算
首先需要激活能运行LAMMPS的环境
conda activate deepmd
然后运行LAMMPS的可执行程序进行计算：
lmp -in in.diffusion
计算结果可视化 2024.06.21
运行完成之后会形成几个文件夹，如下
(deepmd) [root@87_gpu diffusion]# tree
.
├── ave_rdf.txt
├── in.density
├── in.diffusion_msd
├── lmp_msd.py
├── lmp_rdf.py
├── log.lammps
├── msd_Be127_Cu1_hcp.lmp_temp_250
│ ├── Be127_Cu1_hcp.lmp_temp_250.dump
│ ├── Be127_Cu1_hcp.lmp_temp_250.msd
│ ├── Be127_Cu1_hcp.lmp_temp_250.rdf
│ ├── Be127_Cu1_hcp.lmp_temp_250.restart_data
│ ├── log_Be127_Cu1_hcp.lmp_temp_250.lammps
│ └── model_deviBe127_Cu1_hcp.lmp_temp250.out
├── msd_Be127_Cu1_hcp.lmp_temp_450
│ ├── Be127_Cu1_hcp.lmp_temp_450.dump
│ ├── Be127_Cu1_hcp.lmp_temp_450.msd
│ ├── Be127_Cu1_hcp.lmp_temp_450.rdf
│ ├── Be127_Cu1_hcp.lmp_temp_450.restart_data
│ ├── log_Be127_Cu1_hcp.lmp_temp_450.lammps
│ └── model_deviBe127_Cu1_hcp.lmp_temp450.out
├── msd_Be127_Cu1_hcp.lmp_temp_50
│ ├── Be127_Cu1_hcp.lmp_temp_50.dump
│ ├── Be127_Cu1_hcp.lmp_temp_50.msd
│ ├── Be127_Cu1_hcp.lmp_temp_50.rdf
│ ├── Be127_Cu1_hcp.lmp_temp_50.restart_data
│ ├── log_Be127_Cu1_hcp.lmp_temp_50.lammps
│ └── model_deviBe127_Cu1_hcp.lmp_temp50.out
├── msd_Be127_Cu1_hcp.lmp_temp_650
│ ├── Be127_Cu1_hcp.lmp_temp_650.dump
│ ├── Be127_Cu1_hcp.lmp_temp_650.msd
│ ├── Be127_Cu1_hcp.lmp_temp_650.rdf
│ ├── Be127_Cu1_hcp.lmp_temp_650.restart_data
│ ├── log_Be127_Cu1_hcp.lmp_temp_650.lammps
│ └── model_deviBe127_Cu1_hcp.lmp_temp650.out
同样用python脚本对rdf画图
import numpy as np
import matplotlib.pyplot as plt
nbins = 100 # define the number of bins in the RDF
with open("/home/zxg/BeCu/lammps/diffusion/msd_Be127_Cu1_hcp.lmp_temp_250/Be127_Cu1_hcp.lmp_temp_250.rdf", "r") as f: # read the licl.rdf file
lines = f.readlines()
lines = lines[3:]
data = np.zeros((nbins, 7))
count = 0
for line in lines:
nums = line.split()
if len(nums) == 8:
for i in range(1, 8):
data[int(nums[0])-1, i-1] += float(nums[i]) # accumulatie data for each bin
if len(nums) == 2:
count += 1 # count the number of accumulations for each bin
ave_rdf = data / count # calculate the averaged RDF data
np.savetxt('ave_rdf.txt', ave_rdf)
labels = ['Be-Be', 'Be-Cu', 'Cu-Cu']
colors = ['orange', 'purple', 'green']
for i, label, color in zip(range(1, 7, 2), labels, colors):
plt.plot(ave_rdf[:, 0], ave_rdf[:, i], color=color, label=label)
plt.title('RDF')
plt.xlabel('r/Å')
plt.ylabel('g(r)')
plt.legend()
plt.savefig('rdf_250.png', dpi=300)
画图结果展示：
用python做msd的图形脚本为：
import numpy as np
import matplotlib.pyplot as plt
data = np.loadtxt('/home/zxg/BeCu/lammps/diffusion/msd_Be127_Cu1_hcp.lmp_temp_250/Be127_Cu1_hcp.lmp_temp_250.msd', skiprows=2)
time = data[:, 0]
msd1 = data[:, 1]
msd2 = data[:, 2]
plt.plot(time/1000, msd1, 'b-', label='$Be$') # 1fs= 1/1000ps
plt.plot(time/1000, msd2, 'r-', label='$Cu$')
plt.xlabel('time(ps)')
plt.ylabel('$MSD(A^2)$')
plt.title('msd')
plt.legend()
plt.savefig('msd_temp_250.png', dpi=300)
# 计算扩散系数
slope1, residuals = np.polyfit(time, msd1, 1) # 使用 numpy 的 polyfit 函数对 Li+ 和 Cl- 离子的均方位移数据进行一阶多项式（线性）拟合得到斜率
slope2, residuals = np.polyfit(time, msd2, 1) # 这里，我们假设 MSD 与时间成线性关系，即 MSD(t) = 6*D*t，其中 D 是扩散系数。
Diff1 = slope1/6 * 1e-5 # D=1/6*slope; 1 A^2/fs = 1e-5 m^2/s
Diff2 = slope2/6 * 1e-5
print(f"Diffusion Coefficients of Be: {Diff1} m^2/s")
print(f"Diffusion Coefficients of Cu: {Diff2} m^2/s")
结果展示：由于线条较粗，导致Cu的msd覆盖了Be，可通过平均方法降低图片的可观测度。
计算结果显示铜掺杂铍晶体几乎不扩散，即msd也在某个常数值左右来回震荡。这时可通过计算不同温度下的msd并结合相关的理论公式拟合来预测固体材料中的原子扩散。
另外，可以将计算进行一下扩展，即在铜原子附近设置一个缺陷Be空位在研究msd的变化等。
作者：朱雪刚 邮箱: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爱好者在转发推广。感谢大家公众号的推广。
