经典分子动力学的准确程度依赖于其所采用的势能模型的准确程度，虽然第一性分子动力学的结果较为准确，但目前还是较难应用于大体系的计算。因此，从经典分子动力学出发，通过机器学习建立更加准确的势能模型，不仅提高了势能模型的准确程度，而且能够发挥经典分动的优势而应用于大体系的计算，是从计算准确度和计算效率出发看目前较为可行的一种方案。
一、平台选择
机器学习势能最简单的实现方式是通过人工神经网络，基于第一性原理程序制作的数据集，实现模型参数的优化，目前已有的平台已有很多：
1：通过tensorflow自己创建神经网络，但是困难之处在于创建的模型与目前通用的分动模拟程序如lammps等的接口问题，需要有一定的编程基础，时间成本较高，不推荐新手从这里入门；
2：deepMD，已经集成相应的lammps环境，随着版本的发展，结合conda，程序的安装越来越便捷，其对原子近邻环境的描述和其余平台有所区别，通过一个矩阵运算考虑了所有对称性，然后用一个网络来拟合，最后再用一个网络学习势能模型。由于机器学习势能开展分子动力学计算，近邻配位环境的转换是最耗时的一部分，所以不要期望deepmd势能用于分动计算拥有和经典势能一样的速度，但是要快于第一性原理分动计算，不过cp2k也挺快哦。deepmd可以学习能量、力以及virial，力的学习对于开展分子动力学是非常必要的，所以最早出现的aenet用于拟合势能面还可以，用于分动就比较奔溃。最后，最好的就是deepmd集成了lammps，lammps分动计算的功能十分强大，有很多好用的命令和工具，所以省掉开展分动计算没有接口和势能验证方面的问题，适合新手入手。
3：aenet The Atomic Energy Network 不推荐
一个简单的机器学习势神经网络实现，分动结构可以用基于python平台的ASE(Atomic Simulation Environment)，只对能量进行拟合，只是一个全连接层网络，原子近邻环境可以有两种不同的函数组合。没有力的拟合，则基本上分动计算没法实现，只能拟合势能面。若要实现分动，也只能通过ASE后者python的接口，不过不推荐。
4：GPUMD的NEP势能：适合热输运性质的分动计算
如果是热输运方面的分动应用，GPUMD非常推荐，只要有一高端显卡，同样的数据集，其势能分动的可靠性在经验上要优于DeepMD，并且计算速度非常快！特别是针对于热输运性质的研究和计算，GPUMD的表现非常棒。
二、数据集制作
1. 如果想要快速制作数据，推荐使用CP2K开展第一性分动计算，不过CP2K使用类似高斯函数组的基组，可能是基组选择和赝势的设置，对某些二维材料，计算结果和平面波软件的结果有较大差异，比如晶格常数、声子色散不正确等，但是对三维非奇葩材料，算起来准确度没问题，并且优点是很快，缺点是没有平面波软件的倒空间性质，如电子结构的色散等，只有能量和力等计算，作为学习数据集还是可以。
2.如果精度要求较高，可以制作随机位移(通过经典势能的分动或者第一性分动提供），然后用平面波基组的第一性原理软件如vasp,qe等逐帧计算，时间相对较慢，但是准确性可以很高。
3.用VASP等第一性软件的分动提供能量和力的数据集。
数据集的数目参考：一般8000组（超胞结构）以上，能量精度能到1e-2以下，力的精度在5e-2以下，在此精度，通过ALAMODE的3阶力常数拟合能够得到2阶较为准确的声子色散结构，但是3阶力常数的准确性可能在5%~10%，分动计算与网络中参数数目有关，基本正常。基于shengBTE计算热输运性质基本和第一性原理结果接近，但是准确性要提高，可能还需要继续扩充数据集。
三、机器学习势能训练
1、使用deepmd，可以参考新版deepmd的manual，训练起来有两个网络，一个是近邻环境表征用的，一个是全连接层网络学习势能；
应用示例:
1.数据集制作：
采用cp2k第一性分动制作结构-能量数据集(为了提高精度，也可进一步通过QE等第一性计算程序计算相应结构下的能量和力)，将数据集处理为aenet所需的xsf格式文件，下面为python处理脚本，实现cp2k计算完成后轨迹内xyz坐标的读取以及力的读取，并放入独立的文件内用于aenet的学习：
#输入1：cp2k计算结果文件夹：cp2k-pos-1.xyz文件所在路径,可以为列表，1重循环针对此，转为filepath
filepath0=[R"D:\latio3pbnm_cp2k_md_odin2_210328_pbnmXyz1uc_mdNVT300"]
filepath0.append(R"D:\latio3pbnm_cp2k_md_odin2_210328_pbnmXyz1uc_mdNVE300")
#输入2：xyz文件名前缀：latio3pbnm_cp2k_md_cutoff350_relcut60 -pos-1.xyz
file_pre0=[R'\latio3pbnm_cp2k_md_cutoff350_relcut60']
file_pre0.append(R'\latio3pbnm_cp2k_md_cutoff350_relcut60')
#输入3：计算文件2重标识，最终全部按照序号制作成独立文件夹，序号的起点后可控制不为0
file_ident1=['-pos-1']
file_ident2=['.xyz']
#自动计算：
file_num=len(filepath0)
file_ident1_num=len(file_ident1)
file_ident2_num=len(file_ident2)
#输入5：输出文件夹父目录：不存在则会创建
out_path_pre=R"D:/cp2k_latio3_md/aenet_datafiles"
#输入6；创建文件夹序号已有的总数目，即文件夹序号的起点
dir_num0=0
#输入7：aenet用总的generate.in文件:包含"aenet_gener_traindata.in" "aenet_gener_testdata.in"
train_filename=out_path_pre+"\\"+"aenet_gener_traindata.in"
test_filename=out_path_pre+"\\"+"aenet_gener_testdata.in"
#输入8：以上汇总xsf文件中总路径前缀
ryzen_filepath='..' #默认计算时在当前1级目录下创建并进入了一个新的计算目录
#输入6.2: 训练和测试的比例，即训练的总文件数目占总帧数目的比例,将训练和测试分别输出到相应汇总文件中去
#采用10个顺序分配法
train_ind_ratio=0.85 #即0-8 分隔0.85 9
#输入9: 晶格结构，同QE的cell_para格式,为字符串空格分隔，直接字符串形式输出
#5.6336 5.6156 7.9145
vec1='5.6336 0 0'
vec2='0 5.6156 0'
vec3='0 0 7.9145'
#
#输出总文件目录是否存在，不存在则创建
if not os.path.exists(out_path_pre):
os.makedirs(out_path_pre)
#aenet的训练集和测试集
train_file=open(train_filename,'w')
train_line=[]
test_file=open(test_filename,'w')
test_line=[]
#
train_data_num=0
test_data_num=0
dir_ind=0
for k0 in range(file_num):
filepath=filepath0[k0]
file_pre=file_pre0[k0]
for k1 in range(file_ident1_num):
for p1 in range(file_ident2_num):
#自动计算：cp2k的力计算输出文件名:latio3pbnm_cp2k_md_cutoff350_relcut60 -frc-1.xyz
force_filename=filepath+file_pre+'-frc-1.xyz'
#在输出总文件夹中下级目录，将所有文件夹循环和内部2重文件标识循环，组织为起点dir_num0的序列列表
out_path=out_path_pre+'\\'+repr(dir_ind+dir_num0)
file_lable=dir_ind+dir_num0
dir_ind=dir_ind+1
#将2级文件标识转为总输出文件下起始标号0的独立文件序列
if not os.path.exists(out_path):
os.makedirs(out_path)
#1打开xyz文件
filename=filepath+file_pre+file_ident1[k1]+file_ident2[p1]
xyz_file=open(filename,'r')
force_file=open(force_filename,'r')
#1.1读入原子个数
xyz_line=xyz_file.readlines()
xyz_file.close()
force_line=force_file.readlines()
force_file.close()
#总行数
line_num=len(xyz_line)
temp_list=xyz_line[0].split()
atom_num=round(float(temp_list[0]))
#计算得到每一帧占用的行数: 1atomnum+1(i,time,E)+20(La 0.0000000000 0.0000000000 0.0000000000)
#每帧第二行：i为步数，time为时间单位fs，E为pot能量，单位Ha
#原子坐标单位为angstrom
pe_line=atom_num+2
#帧数
pic_num=round(line_num/pe_line)
for k2 in range(pic_num):
#创建写文件,每帧一个原子坐标.txt文件，不新建独立文件夹用于后期QE计算；
#aenet用，每帧有独立文件夹
d_now2=out_path+'\\'+repr(k2)+'\\'+'xsf'
#linux计算汇总内路径标识
d_now3=ryzen_filepath+'/'+repr(file_lable)+'/'+repr(k2)+'/'+'xsf'
if not os.path.exists(d_now2):
os.makedirs(d_now2)
#1qe计算用坐标文本文件编号命名并全部放置在总目录下
out_file_name=out_path+'\\'+repr(k2)+'.txt'
qe_file=open(out_file_name,'w')
#坐标信息存储
atom_xyz=[] #字符串
force_xyz=[] #字符串拆分成list,转为float后单位转换
for k3 in range(atom_num):
res_line=xyz_line[k2*pe_line+2+k3].split()
force_res_line=force_line[k2*pe_line+2+k3].split()
#转为qe用 type x y z
qe_line=res_line
#字符串形式用于aenet
atom_xyz.append(xyz_line[k2*pe_line+2+k3].strip())
#力的单位转换从Ha/bohr----eV/angstrom
force_xyz.append([float(force_res_line[1])*2*13.60569/0.52918,float(force_res_line[2])*2*13.60569/0.52918,float(force_res_line[3])*2*13.60569/0.52918])
print("%s %s %s %s"%(qe_line[0],qe_line[1],qe_line[2],qe_line[3]),file=qe_file)
qe_file.close()
#能量
#2制作aenet数据集文件
#在每帧文件夹d_now2下创建帧编号命名.xsf文件
xsf_filename=d_now2+'\\'+repr(k2)+".xsf"
xsf2_filename=d_now3+'/'+repr(k2)+".xsf"
xsf_file=open(xsf_filename,'w')
#将能量Ha转为eV 13.60569
# i = 6, time = 6.000, E = -553.2237210041
ene_temp_line=xyz_line[k2*pe_line+1].split()
ene_res=float(ene_temp_line[8])*13.60569*2
#make xsf file
print('# total energy = %.12f eV'%(ene_res),file=xsf_file)
print("\nCRYSTAL\nPRIMVEC",file=xsf_file)
print(" %s"%(vec1),file=xsf_file)
print(" %s"%(vec2),file=xsf_file)
print(" %s"%(vec3),file=xsf_file)
print("PRIMCOORD\n%.0f 1"%(atom_num),file=xsf_file)
for k13 in range(atom_num):
print("%s %.14f %.14f %.14f"%(atom_xyz[k13],force_xyz[k13][0],force_xyz[k13][1],force_xyz[k13][2]),file=xsf_file)
xsf_file.close()
#归类到总的训练集或测试集中,每个总文件夹独立分配
if (k2%10)/10<=train_ind_ratio:
train_data_num=train_data_num+1
train_line.append(xsf2_filename)
else:
test_data_num=test_data_num+1
test_line.append(xsf2_filename)
#则总输出文件下/0起始(2个文件标识命名的计算次)/qe用坐标.txt+帧序号独立文件夹/xsf/帧序号.xsf
#输出generate用文件
print("FILES\n%d"%(train_data_num),file=train_file)
print("FILES\n%d"%(test_data_num),file=test_file)
for k14 in range(train_data_num):
print("%s"%(train_line[k14]),file=train_file)
for k15 in range(test_data_num):
print("%s"%(test_line[k15]),file=test_file)
train_file.close()
test_file.close()
2.linux shell计算脚本
#!/bin/bash
file_pre="aenet_latio3"
file_comp="redlinux"
file_date="210402"
######################AENET###############################
############redlinux-gnu_openmpi
export PATH="openmpi的路径:$PATH"
export LD_LIBRARY_PATH="openmpi的lib:$LD_LIBRARY_PATH"
#aenet:path
export PATH="aenet的路径aenet-2.0.4/bin:$PATH"
#aenet安装后相应的执行文件名
red_gener="generate.x-2.0.4-gfortran_openblas_mpi"
red_train="train.x-2.0.4-gfortran_openblas_mpi"
red_pred="predict.x-2.0.4-gfortran_openblas_mpi"
################################
gener=${red_gener}
train=${red_train}
pred=${red_pred}
mpi_pre="mpirun "
#stp_datafiles="aenet_stp_datafiles"
train_files="../aenet_gener_traindata.in"
test_files="../aenet_gener_testdata.in"
aenet_files=${train_files}
##########################################################
hidden1_sh=20
hidden2_sh=20
hidden_sh="2 ${hidden1_sh}:tanh ${hidden2_sh}:tanh"
###run_j
run_gener="F"
run_train="F"
run_pred="T"
###
workdir="${file_pre}_${file_comp}_${file_date}"
if test ! -d ${workdir}
then
echo "${workdir} not exist and make!!!"
mkdir ${workdir}
else
echo "${workdir} already exist and not make!!!"
fi
cd ${workdir}
#cp ../${stp_datafiles}/*.stp ./
for f_ident1 in "cp2k_nvt300_${hidden1_sh}t${hidden2_sh}t"
#for f_ident1 in 1
do
file_count="${f_ident1}"
gener_filename_in="./generate_${file_pre}_${file_comp}_${file_date}_${file_count}.in"
gener_filename_out="./generate_${file_pre}_${file_comp}_${file_date}_${file_count}.out"
train_filename_in="./train_${file_pre}_${file_comp}_${file_date}_${file_count}.in"
train_filename_out="./train_${file_pre}_${file_comp}_${file_date}_${file_count}.out"
pred_filename_in="./pred_${file_pre}_${file_comp}_${file_date}_${file_count}.in"
pred_filename_out="./pred_${file_pre}_${file_comp}_${file_date}_${file_count}.out"
#
#1:atomtype.stp
setup_sh="SETUPS"
for atom_type in "Ti" "O" "La"
do
if test ! -r ${atom_type}.fingerprint.stp
then
cat > ${atom_type}.fingerprint.stp << EOF1
DESCR
N. Artrith and A. Urban, Comput. Mater. Sci. 114 (2016) 135-150.
N. Artrith, A. Urban, and G. Ceder, Phys. Rev. B 96 (2017) 014112.
END DESCR
ATOM ${atom_type}
ENV 3
La
Ti
O
RMIN 0.55d0
BASIS type=Chebyshev
radial_Rc = 8.0 radial_N = 16 angular_Rc = 8.0 angular_N = 8
EOF1
fi
setup_sh="${setup_sh}
${atom_type} ${atom_type}.fingerprint.stp"
done
#2:
cat > ${gener_filename_in} << EOF0
OUTPUT ${file_count}_generate_res
TYPES
3
Ti -1577.203929041 | eV
La -859.121667221 | eV
O -429.356978284 | eV
${setup_sh}
EOF0
#put dataset path into .in file
cat ${aenet_files} >> ${gener_filename_in}
if [ ${run_gener} = "T" ]
then
echo "running generate.x as ${mpi_pre} ${gener} ${gener_filename_in} > ${gener_filename_out}>>>>>begin>>>>>"
${gener} ${gener_filename_in} > ${gener_filename_out}
echo "running generate.x done!!!>>>>>>>end<<<<<<<"
else
echo "Not running generate.x!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!"
fi
cat > ${train_filename_in} << EOF4
TRAININGSET ${file_count}_generate_res
TESTPERCENT 10
ITERATIONS 20
#MAXENERGY 1.0
TIMING
SAVE_ENERGIES
METHOD
bfgs
!METHOD 2 Steepest Descent
!online_sd gamma=5.0d-7 alpha=0.25d0
!online_sd gamma=1.0d-8 alpha=0.25d0
!
!METHOD 3 Extended Kalman Filter
!ekf
!
!METHOD 4 Levenberg-Marquardt
!lm batchsize=8000 learnrate=0.1d0 iter=3 conv=0.001 adjust=5.0
NETWORKS
! atom network hidden
! types file-name layers nodes:activation
Ti Ti.${file_count}.nn ${hidden_sh}
O O.${file_count}.nn ${hidden_sh}
La La.${file_count}.nn ${hidden_sh}
EOF4
if [ ${run_train} = "T" ]
then
echo "running generate.x as ${mpi_pre} ${train} ${train_filename_in} > ${train_filename_out}>>>>>begin>>>>>"
${mpi_pre} ${train} ${train_filename_in} > ${train_filename_out}
echo "running generate.x done!!!>>>>>>>end<<<<<<<"
else
echo "Not running train.x !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!"
fi
#
cat > ${pred_filename_in} << EOF5
TYPES
3
Ti
La
O
NETWORKS
Ti Ti.${file_count}.nn
O O.${file_count}.nn
La La.${file_count}.nn
FORCES
# or optimize coordinates:
#
# RELAX
# method=bfgs F_conv=1.0d-2 E_conv=1.0d-6 steps=99
#
# method: optimization method (currently only BFGS)
# F_conv: convergence thershold for the forces
# E_conv: convergence threshold for the energy
# steps: max. number of iterations
EOF5
#put dataset path into predict.in file
cat ${test_files} >> ${pred_filename_in}
if [ ${run_pred} = "T" ]
then
echo "running predict.x as ${mpi_pre} ${pred} ${pred_filename_in} > ${pred_filename_out}>>>>>begin>>>>>"
${pred} ${pred_filename_in} > ${pred_filename_out}
echo "running predict.x done!!!>>>>>>>end<<<<<<<"
else
echo "Not running predict.x !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!"
fi
cd ../
done
扫一扫
2123
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
