任何命令的参数解释可使用command --help或man command）
==cd /share/apps/== 下面有超算可以使用的软件
==cd==到家目录，==ls -al== 可看到有shell配置文件.bashrc，==vim .bashrc== 进入编辑
按i或o进入编辑模式，输入以下内容（以使用vasp软件为例）
#!/bin/bash
export PATH=/share/apps/vasp/vasp.5.4.4/bin/:$PATH #bin下为可执行文件
export LD_LIBRARY=/share/apps/vasp/vasp.5.4.4/vasp.5.lib/:$LD_LIBRARY #lib下为库文件Esc返回命令模式，输入：wq保存退出
重新登陆服务器或者==source ~/.bashrc== 环境变量生效
（比使用ftp方便）
rz：把window文件传送到linux服务器
rz 不会覆盖同名文件 rz -bye 覆盖同名文件
sz：把linux上的文件传送到window
==sz file1 file2 dir/*==
（节点为并行节点）
新建一个sh脚本，==vim sub_task.sh==输入以下内容：
#!/bin/bash
#BSUB -q mpi #使用mpi队列
#BSUB -n 24 #指定进程数
#BSUB -o %J.out #指定标准输出文件，%J指用作业的ID作为文件名
#BSUB -e %J.err #指定标准错误文件
export OMP_NUM_THREADS 2
export MP_TASK_AFFINITY=core:$OMP_NUM_THREADS
mpirun /share/apps/vasp/vasp.5.4.4/bin/vasp_std #运行可执行文件
退出后添加可执行权限
==chmod +x sub_track.sh==
提交任务：==bsub < sub_track.sh==
检查任务状态：==bjobs==
结束任务进程：==bkill jobid== ==bkill -r jobid== （-r为强制结束，谨慎使用）
查看队列：==bqueues==
详细命令介绍看手册
==rm -rf !(file1 | file2)==
在家目录：==source .bashrc== 或 ==. .bashrc==
不在家目录：==source /.bashrc== 或 ==/.bashrc==
用于设置指令的别名（在.bashrc中设置）
利用脚本处理文本文件（很有用）
sed -e xxxxx或者 sed -f xxxxx 只能将内容输出到屏幕或文件，不能对文件进行操作
sed -i ‘4d’ ./INCAR 删除文件第四行
sed -i ‘4cISPIN = 2’ ./INCAR 将文件第四行替换为ISPIN = 2
sed -i ‘4iISPIN = 2’ ./INCAR 在文件第四行后面插入ISPIN = 2
sed -i ‘4aISPIN = 2’ ./INCAR 在文件第四行前面插入ISPIN =2
sed -i ‘$aENCUT = 400’ ./INCAR 在文件当前输入位置$符后插入ENCUT = 400
sed -i ‘3,4s/that/this/g’ ./INCAR 在文件第三、四行寻找that替换为this（不要引号也可）
不使用单（双）引号则用\反斜杠分割，用反斜杠需要保证插入的语句是连续的不能有空格，如：sed -i 4i\ISPIN=2 ./INCAR。
或者不加反斜杠也行，如：sed -i 4iISPIN=2 ./INCAR若有空格则需要引号
${}：对变量的替换，同$var
$()：对命令的替换，同``（反引号，波浪号下面）
eg：
==for i in {1..8}; do cp 400 $((400+50*****$i)); sed -i “s/400/$((400+50*****$i))/g” $((400+50*****$i))/INCAR; done==
==for i in *; do cd $i; sub < sub_task.sh; cd $OLDPW; done==
（这里需进入vasp输入文件的目录，这样执行vasp才知道将输出文件储存在当前目录，否则会报错）
（处理文本文件）
awk ‘{print $1,$4}’ log.txt 默认每行按空格或TAB分割，出现空格或TAB就当作一项，输出每一行的第1、4
awk ‘{printf “%-8s %-10s\n”,$1,$4}’ log.txt 格式化输出
awk -F, ‘{print $1 $2}’ log.txt -F指定分割字符为,逗号，输入每一行的1、2项
awk ‘BEGIN{FS=”,”} {print $1,$2}’ log.txt 用内建变量表示分割
awk -F ‘[ ,]’ ‘{print $1,$2,$5}’ log.txt 使用多个分割符，先使用空格分割，然后对分割结果再用“,”分割
awk -va=1 ‘{print $1,$1+a}’ log.txt -v设置变量，此处设置变量a等于1，输出每一行第一行第1项和第1项+1（若第1项为数字则进行运算，若为字符则在其后加上数字1）
awk -va=1 -vb=s ‘{print $1,$1+a,$1+b}’ log.txt
eg:
for i in *; do echo -e $i “\t” $(grep User $i/OUTCAR | awk ‘{print $4}’); done
考虑ENCUT值的选择对计算时间的影响：遍历每个文件夹，输出文件夹名、制表符以及OUTCAR中任务运行时间
例1：
通过awk命令得到计算时间将数据保存为data文件
==for i in *; do echo -e
使用bash脚本: ==vim plt.py==
#!/usr/bin/env python #指定脚本解释程序
import matplotlib.pyplot as plt
import numpy as np
x,y = np.load_text(‘data.txt’,delimeter = ‘,’,usecols = (0,1),unpack = True) #分隔符为,逗号
plt.xlabel(‘ENCUT/eV’)
plt.ylabel(‘Time/s’)
plt.plot(x,y, linewidth = 2)
plt.show()==chmod +x plt.py== 添加执行权限或使用python解释器打开 ==python plt.py==
例2：
通过awk命令得到电子步优化最终能量与ENCUT关系
==for i in *; do echo -e
（1）os模块的system方法
注意python的变量需要经过转化才能作为Linux的变量 os.environ['S']=str(S)
其它参考下面的链接或本人编写的批量初始化代码即可
（2）os模块popen方法
（3）subprocess模块
引用：https://blog.csdn.net/ronon77/article/details/84774575
引用：https://wenwen.sogou.com/z/q778923295.htm
引用：http://www.zzvips.com/article/67036.html
引用：https://www.jb51.net/article/63897.htm
引用：https://www.cnblogs.com/momoyan/p/9145992.html
处理很多数据的时候，感觉没有直接用python，然后利用Shell与Python交互方便
引用：https://blog.51cto.com/nanwangting/932095
（1）使用$0,$1,$2...
$0 —当前脚本文件名
$n —传递给脚本的参数，$1表示第一个参数，$2表示第二个参数 ... 超过10个需要使用${10},${11}...
#!/usr/bin/bash
echo "脚本$0"
echo "脚本第一个参数$1"
echo "脚本第二个参数$2"执行：==./test.sh 1 2==
输出：脚本./test.sh
脚本第一个参数1
脚本第二个参数2
（2）getops
一般不会出现，本服务器查看硬盘限额命令为 ==mmlsquota==
引用：https://blog.51cto.com/13438667/2082594
==cp ./{INCAR,POSCAR,KPOINTS,POTCAR} ./data==
==du -h==
（1）方法一：
==echo -n "不换行输出"==
（2）方法二：
==echo -e "字符串\c"==
==echo -e 处理特殊字符;==
可接的特殊字符有：
\c 最后不加上换行符号；
\f 换行但光标仍旧停留在原来的位置；
\n 换行且光标移至行首；
\r 光标移至行首，但不换行；
\t 插入tab；
\v 与\f相同；
\ 插入\字符；
引用：https://blog.csdn.net/wteruiycbqqvwt/article/details/98988462
电子存在四个量子数，n、l、m、n，其中m为自旋量子数，取值为+1/2和-1/2。量子化学中用自旋多重度来区分一组简并的波函数，这些波函数之间只存在自旋角动量的不同。
自旋多重度定义为2S+1，其中S是自旋角动量，S与体系内单电子数N有关，S=N/2，即自旋多重度=N+1。
1）如果体系所有的占据轨道的两个自旋轨道都是充满的，根据泡利不相容原理，两电子自旋相反，体系没有自旋角动量。故只存在一种自旋状态（无自旋），自旋多重度为1。
2）如果体系有一个占据轨道只占据了一个自旋轨道，另一个自旋轨道为空，体系分为两部分，一部分是由没有自旋的全部充满的轨道构成，另一部分就是这个单电子产生的+1/2或-1/2的自旋，自旋多重度为2。
3）如果体系有两个以上的单电子，根据电子排布规则，成单电子优先占据与其他单电子自旋轨道方向一致的自旋轨道，自旋不能抵消，故体系一部分为充满的轨道，一部分为单电子占据的轨道，自旋多重度为N+1。
4）每个电子产生1/2的自旋，由于自旋多重度=N+1=2*****[N*****1/2]+1=2S+1。
5）引用：http://blog.sciencenet.cn/blog-671981-639926.html
1）在VASP中，控制自旋的参数是INCAR中的ISPIN，ISPIN=1进行非自旋极化计算，ISPIN=2进行自旋极化计算，为了进一步得到我们想要的组态，可通过NUPDOWN指定自旋向上和自旋向下的电子数差。
（在进行固定自旋多重态的计算中，如果初始的电荷密度的自旋极矩和INCAR中设置的NUPDOWN不同的话，收敛会很慢。如果是从读入的初始波函数开始计算，则没有类似问题，且结果一直准确。ICHARGE=2表示从电荷密度开始计算，VASP 会自动改设定从MAGMOMT to NUPDOWN/NIONS；ICHARGE=1表示从CHGCAR文件中的电荷密度开始计算，但初始自旋极矩经常不准确。）
2）通过ISMEAR=-2控制，设置FERWE、FERDO实现。引用：
https://wiki.kfki.hu/nano/Easy_manual_occupancy_of_Kohn-Sham_levels_with_FERWE_and_FERDO
（1）raw文件
需要原子种类、模拟的盒子、原子坐标和原子受力、系统能量和维里应力的信息
单位：（不同软件可能采用的单位体系不一样，需要转换，这与Lammps中设置不同的units需转换一个道理）
| Property | Time | Length | Energy | Force | Pressure |
|---|---|---|---|---|---|
| Unit | 皮秒 | 埃 | eV | eV/埃 | Bar |
对应需要准备五个文件：
box.raw, coord.raw, force.raw, energy.raw and virial.raw
以力文件为例：
$ cat force.raw
-0.724 2.039 -0.951 0.841 -0.464 0.363
6.737 1.554 -5.587 -2.803 0.062 2.222
-1.968 -0.163 1.020 -0.225 -0.789 0.343
#一行代表一个frame(即一个构型)，例子中有2个原子，每个原子3个方向受力，共6列
#其它raw软件同理，coord.raw(有frame行，atom*3列) box.raw(frame行，盒子三个方向9列) virial.raw(frame行，维里应力三个方向9列)还需要原子种类文件：type.raw
$ cat type.raw
0 1
#一行内原子种类用整数依次表示（保证与原子坐标准备文件中的原子顺序一致）如Si02即为：0 1 1 （2）set文件
把raw文件变成deepmd训练需要的set文件
使用命令：$deepmd_source_dir/data/raw/raw_to_set.sh
$ ls
box.raw coord.raw energy.raw force.raw type.raw virial.raw
$ $deepmd_source_dir/data/raw/raw_to_set.sh 2000
nframe is 6000 #共6000个frame，分成3个集
nline per set is 2000
will make 3 sets
making set 0 ...
making set 1 ...
making set 2 ...
$ ls
box.raw coord.raw energy.raw force.raw set.000 set.001 set.002 type.raw virial.raw #set.002是测试集参数文件（json格式）：
四个部分（model、learning_rate、loss、training）
以水分子为例：
"model": {
"type_map": ["O", "H"],
"descriptor" :{ #设置描述符D
"type": "se_a", #smooth-edition angular infomation
"rcut_smth": 5.80, #周围原子对中心原子的贡献从5.80开始降为0
"rcut": 6.00, #周围原子对中心原子的贡献到6.00降为0 (截断半径)
"sel": [46, 92], #截断半径内所能容纳的最大的原子数（46个氧原子，92个氢原子）
"neuron": [25, 50, 100], #神经网络大小
"axis_neuron": 16, #embedding matrix的宽度
"resnet_dt": false, #描述符的残参网络设置 根据经验设置
"seed": 1, #随机数种子
"_comment": " that's all"
},
"fitting_net" : {
"neuron": [240, 240, 240],
"resnet_dt": true,
"seed": 1,
"_comment": " that's all"
},
"_comment": " that's all"
}
"loss" : {
"start_pref_e": 0.02, #开始能量误差的权重系数
"limit_pref_e": 1, #最终能量误差的权重系数
"start_pref_f": 1000, #开始力的权重系数 一开始几乎全部以力的误差作为损失函数降低的贡献
"limit_pref_f": 1,
"start_pref_v": 0, #维里应力 这里简单不考虑维里应力，计算固体考虑弹性常数要打开
"limit_pref_v": 0,
"_comment": " that's all"
}
"learning_rate" :{
"type": "exp", #指数方式衰减
"start_lr": 0.005, #一开始的lr
"decay_steps": 5000, #每5000步衰减一次
"decay_rate": 0.95, #衰减到原来的95% 新的版本里面有stop_lr,可自己计算衰减率
"_comment": "that's all"
}
"training" : {
"systems": ["../data/"], #训练文件夹的位置
"set_prefix": "set", #
"stop_batch": 1000000, #训练步数
"batch_size": 1, #同时训练的数量
"seed": 1,
"_comment": " display and restart",
"_comment": " frequencies counted in batch",
"disp_file": "lcurve.out", #指定文件记录损失函数变化
"disp_freq": 100, #隔多少步输出损失函数变化
"numb_test": 10, #测试数量
"save_freq": 1000, #每1000步输出临时文件，以防止意外结束从头计算
"save_ckpt": "model.ckpt", #
"load_ckpt": "model.ckpt",
"disp_training":true,
"time_training":true,
"profiling": false,
"profiling_file":"timeline.json",
"_comment": "that's all"
}Window键+x 菜单选择Windows Powershell（管理员）输入notepad会跳出记事本
启动不自动进入conda环境：
==conda config --set auto_activate_base false==
（deepmd需要在conda环境下，而vasp不能在conda环境下，否则不能正常计算）
激活环境：==conda activate base== 取消激活：==conda deactivate==
安装软件AutoHotKey（较简单）
AutoHotKey是一款著名的windows系统快捷键设置的软件，轻便小巧。
官方下载: https://autohotkey.com/download/ahk-install.exe
（1）先安装AutoHotKey
（2）打开记事本，把如下内容复制粘贴进去：
; Typora
; 快捷增加字体颜色
; SendInput {Text} 解决中文输入法问题
#IfWinActive ahk_exe Typora.exe
{
; Ctrl+Alt+O 橙色
^!o::addFontColor("orange")
; Ctrl+Alt+R 红色
^!r::addFontColor("red")
; Ctrl+Alt+B 浅蓝色
^!b::addFontColor("cornflowerblue")
}
; 快捷增加字体颜色
addFontColor(color){
clipboard := "" ; 清空剪切板
Send {ctrl down}c{ctrl up} ; 复制
SendInput {TEXT}<font color='%color%'>
SendInput {ctrl down}v{ctrl up} ; 粘贴
If(clipboard = ""){
SendInput {TEXT}</font> ; Typora 在这不会自动补充
}else{
SendInput {TEXT}</ ; Typora中自动补全标签
}
}（3）将文件保存为ahk后缀的文件，如TyporaHotKey.ahk
（4）双击运行
（5）在Typora软件里就可以使用快捷键：
如按Ctrl+Alt+O添加橙色，Ctrl+Alt+R 红色，按Ctrl+\取消样式！
也可以右键 MyHotkeyScript.ahk 脚本文件，点击Compile Script编译脚本成exe程序，就可以不用下载Autohotkey在其他电脑上运行了。
上面脚本只写了橙色、红色、浅蓝三种颜色，你可以按需照例增加其他颜色或快捷方式！
引用：https://blog.csdn.net/superit401/article/details/106344453/
Typora 有一个「高亮」的格式，类似于荧光笔，但是感觉默认的颜色偏亮，看久了不舒服，所以利用修改主题文件的方式来自定义颜色。
操作很简单，先找到主题文件：「文件」 ==> 「偏好设置」（或者直接 Ctrl + 逗号），在右边「外观」栏中找到「打开主题文件」打开：
mark {
background: #a9d18e;
border-bottom: 0px solid #ffffff;
padding: 0.0px;
margin: 0 0px;
}如果只是想要单纯改变颜色，也可以只写 background 的这一行：
mark {
background: #a9d18e;
}像这样：
其中 backgraoud 后面的十六进制数为所需要的颜色。border-bottom 是下划线的大小和颜色。padding 就是上下左右的边框大小。margin 就是所标记文字离左右文字的距离。
这里做一个各个参数的对比：
Ctrl+逗号 打开「偏好设置」，在「通用」里最下面打开高级设置，找到下图位置，添加自己需要的快捷键：
引用：https://blog.csdn.net/weixin_39679367/article/details/104391070
from scipy.optimize import minimize
import numpy as np
from matplotlib import pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
from matplotlib import cm as cm
import pandas as pd
from scipy import signal
class GPR:
def __init__(self, optimize=True):
self.is_fit = False
self.train_X, self.train_y = None, None
self.params = {"l": 0.4, "sigma_f": 0.1}
self.optimize = optimize
def fit(self, X, y):
# store train data
self.train_X = np.asarray(X)
self.train_y = np.asarray(y)
# hyper parameters optimization
def negative_log_likelihood_loss(params):
self.params["l"], self.params["sigma_f"] = params[0], params[1]
Kyy = self.kernel(self.train_X, self.train_X) + 1e-8 * np.eye(len(self.train_X))
return 0.5 * self.train_y.T.dot(np.linalg.inv(Kyy)).dot(self.train_y) + 0.5 * np.linalg.slogdet(Kyy)[
1] + 0.5 * len(self.train_X) * np.log(2 * np.pi)
if self.optimize:
res = minimize(negative_log_likelihood_loss, [self.params["l"], self.params["sigma_f"]],
bounds=((1e-4, 1e4), (1e-4, 1e4)),
method='L-BFGS-B')
self.params["l"], self.params["sigma_f"] = res.x[0], res.x[1]
self.is_fit = True
def predict(self, X):
if not self.is_fit:
print("GPR Model not fit yet.")
return
X = np.asarray(X)
Kff = self.kernel(self.train_X, self.train_X) # (N, N)
Kyy = self.kernel(X, X) # (k, k)
Kfy = self.kernel(self.train_X, X) # (N, k)
Kff_inv = np.linalg.inv(Kff + 1e-8 * np.eye(len(self.train_X))) # (N, N)
mu = Kfy.T.dot(Kff_inv).dot(self.train_y)
cov = Kyy - Kfy.T.dot(Kff_inv).dot(Kfy)
return mu, cov
def kernel(self, x1, x2):
dist_matrix = np.sum(x1 ** 2, 1).reshape(-1, 1) + np.sum(x2 ** 2, 1) - 2 * np.dot(x1, x2.T)
return self.params["sigma_f"] ** 2 * np.exp(-0.5 / self.params["l"] ** 2 * dist_matrix)
def rotate(angle):
ax.view_init(azim=angle)
IO = "C:\\Users\\haoji\\Desktop\\config_data.xlsx"
data = pd.read_excel(IO,header=None)
train_X = data.values[:,0:2].tolist()
train_y = data.values[:,2].tolist()
gpr = GPR(optimize=True)
gpr.fit(train_X,train_y)
test_d1 = np.arange(1,5, 0.1)
test_d2 = np.arange(1,5, 0.1)
test_d1, test_d2 = np.meshgrid(test_d1, test_d2)
test_X = [[d1, d2] for d1, d2 in zip(test_d1.ravel(), test_d2.ravel())]
mu, cov = gpr.predict(test_X)
z = mu.reshape(test_d1.shape)
fig = plt.figure(figsize=(14, 10))
ax = Axes3D(fig)
_ = ax.plot_surface(test_d1, test_d2, z, cmap=cm.coolwarm, linewidth=0, alpha=0.4, antialiased=False)
ax.scatter(np.asarray(train_X)[:,0], np.asarray(train_X)[:,1], train_y, c=train_y, cmap=cm.coolwarm)
#ax.contourf(test_d1, test_d2, z, zdir='z', offset=0, cmap=cm.coolwarm, alpha=0.7)
ax.set_title("l=%.2f sigma_f=%.2f" % (gpr.params["l"], gpr.params["sigma_f"]))
ax.set_xlabel("R1/Angstrom")
ax.set_ylabel("R2/Angstrom")
ax.set_zlabel("Energy/eV")
#elev,azim = 90,0
#ax.view_init(elev,azim)（该条目中相关的bash脚本仅适用于本人实际操作时的情况，应根据自己的需要和实际文件格式等进行调整）
（1）初次实现批量构建准备文件（产生的Si-O1-O2三个原子在同一直线的构型）
（2）批量提交任务
批量提交的命令：
提交任务的脚本:
（3）提取数据
得到的result文件用第3小节的代码处理
（4）批量构建准备文件（利用与python交互实现R1_R2_theta三个变量的准备文件构建）
（5）利用脚本参数传递分批提交任务
（1）此代码为处理服务器批量导出的数据result
result_list = list()
count = 0
with open('C:\\Users\\haoji\\Desktop\\result') as result:
for content in result.readlines():
content_array = content.split(" ")
for unit in content_array:
if unit.endswith('\n'):
print(unit[0:-1],end=' ')
else:
print(unit,end=' ')
count += 1
if count == 2:
print()
count = 0 利用bash脚本初步实现批量构建准备文件及批量提交任务、批量提取数据
（1）批量构建准备文件
bash脚本如下：
得到文件：
输入文件有：
INCAR (指示计算内容和参数)
INCAR 中没有考虑自旋 ISPIN参数没打开，而我们选的体系SiO2不是根据晶格生成的 (Si原子配位为4)，而是直接计算盒子中三个孤立原子的相互作用，按道理这里的单电子数不为零，需要考虑自旋。我不是很清楚这里究竟需不需要考虑自旋，我开了ISPIN后计算的能量不一样，具体还需要老师咨询下DFT专家，同时输出文件OUTCAR中关于每个自旋轨道电子的排布究竟怎么样才是合理的我也不是很清楚。
为了只得到三个原子的相互作用，盒子要足够大，尽量减少周期性的影响。我在POSCAR中测试了合适的盒子尺寸参数，且由于盒子中原子距离变化后存在不收敛的问题，导致部分任务失败，我测试了电子迭代步数NELM保证所有任务收敛。
POSCAR (坐标信息)
这是计算的另一个关键，需要调整坐标的位置得到不同的构型。目前我只是通过bash脚本实现了直线型的，bash脚本有很多限制 (如之前遇到命令不支持浮点数运算等)，所以要再引入$\theta$参数 (三个原子相对坐标信息由二维->三维) 用bash脚本可能不太行。后面可能要尝试用Python来写构造数据，看看能不能实现Python和shell的交互。
POTCAR (赝势文件)
KPOINTS (k文件)
（2）批量任务提交和数据处理
采用以下命令进行批量任务提交
任务提交脚本如下:
计算任务过程中出现了并行效率低的提示，暂时没找到处理的方法，不影响计算。
输出文件中重点是OUTCAR文件，我通过探索实现用以下命令批量获取能量数据。过程中两次批量任务的经命令得到的数据格式不知为何不统一，通过相关命令实现了统一。
最后得到数据通过相关python代码进行再处理导入excel文件中：
利用之前编写的高斯回归代码尝试对数据进行拟合，绘制包含两个自变量的R1、R2的势能面，只是进行训练，没有进行测试。
代码如下：
from scipy.optimize import minimize
import numpy as np
from matplotlib import pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
from matplotlib import cm as cm
import pandas as pd
from scipy import signal
class GPR:
def __init__(self, optimize=True):
self.is_fit = False
self.train_X, self.train_y = None, None
self.params = {"l": 0.4, "sigma_f": 0.1}
self.optimize = optimize
def fit(self, X, y):
# store train data
self.train_X = np.asarray(X)
self.train_y = np.asarray(y)
# hyper parameters optimization
def negative_log_likelihood_loss(params):
self.params["l"], self.params["sigma_f"] = params[0], params[1]
Kyy = self.kernel(self.train_X, self.train_X) + 1e-8 * np.eye(len(self.train_X))
return 0.5 * self.train_y.T.dot(np.linalg.inv(Kyy)).dot(self.train_y) + 0.5 * np.linalg.slogdet(Kyy)[
1] + 0.5 * len(self.train_X) * np.log(2 * np.pi)
if self.optimize:
res = minimize(negative_log_likelihood_loss, [self.params["l"], self.params["sigma_f"]],
bounds=((1e-4, 1e4), (1e-4, 1e4)),
method='L-BFGS-B')
self.params["l"], self.params["sigma_f"] = res.x[0], res.x[1]
self.is_fit = True
def predict(self, X):
if not self.is_fit:
print("GPR Model not fit yet.")
return
X = np.asarray(X)
Kff = self.kernel(self.train_X, self.train_X) # (N, N)
Kyy = self.kernel(X, X) # (k, k)
Kfy = self.kernel(self.train_X, X) # (N, k)
Kff_inv = np.linalg.inv(Kff + 1e-8 * np.eye(len(self.train_X))) # (N, N)
mu = Kfy.T.dot(Kff_inv).dot(self.train_y)
cov = Kyy - Kfy.T.dot(Kff_inv).dot(Kfy)
return mu, cov
def kernel(self, x1, x2):
dist_matrix = np.sum(x1 ** 2, 1).reshape(-1, 1) + np.sum(x2 ** 2, 1) - 2 * np.dot(x1, x2.T)
return self.params["sigma_f"] ** 2 * np.exp(-0.5 / self.params["l"] ** 2 * dist_matrix)
def rotate(angle):
ax.view_init(azim=angle)
IO = "C:\\Users\\haoji\\Desktop\\config_data.xlsx"
data = pd.read_excel(IO,header=None)
train_X = data.values[:,0:2].tolist()
train_y = data.values[:,2].tolist()
gpr = GPR(optimize=True)
gpr.fit(train_X,train_y)
test_d1 = np.arange(1,5, 0.1)
test_d2 = np.arange(1,5, 0.1)
test_d1, test_d2 = np.meshgrid(test_d1, test_d2)
test_X = [[d1, d2] for d1, d2 in zip(test_d1.ravel(), test_d2.ravel())]
mu, cov = gpr.predict(test_X)
z = mu.reshape(test_d1.shape)
fig = plt.figure(figsize=(14, 10))
ax = Axes3D(fig)
_ = ax.plot_surface(test_d1, test_d2, z, cmap=cm.coolwarm, linewidth=0, alpha=0.4, antialiased=False)
ax.scatter(np.asarray(train_X)[:,0], np.asarray(train_X)[:,1], train_y, c=train_y, cmap=cm.coolwarm)
#ax.contourf(test_d1, test_d2, z, zdir='z', offset=0, cmap=cm.coolwarm, alpha=0.7)
ax.set_title("l=%.2f sigma_f=%.2f" % (gpr.params["l"], gpr.params["sigma_f"]))
ax.set_xlabel("R1/Angstrom")
ax.set_ylabel("R2/Angstrom")
ax.set_zlabel("Energy/eV")
#elev,azim = 90,0
#ax.view_init(elev,azim) 得到如下图像：
目前的问题是不知道如何调用函数去获取曲面的切面，得到二维的图像便于分析。其次目前只是采用高斯回归进行简单的测试，后面还需学习神经网络算法。
目前在学习李航的《统计学习方法》，看网上的相关机器学习视频加强自己的理论基础，并尝试编写相关代码巩固。我想对理论有个较好的理解后才能更好的利用deepmd的神经网络。
- 利用python实现R_1、R_2、\theta三个参数的POSCAR文件构建，由直线型拓展到平面并进行计算。
- 开始尝试DeepMD平台的操作流程
- 继续机器学习相关理论的学习
- 下周有许多课程汇报任务，可能时间会比较紧张，请老师谅解
利用python与shell脚本的交互，实现了R1、R2、theta三个参数的POSCAR文件等的构建。其中R1、R2选取的范围为**[0.6,2.8] Angstrom**，间隔为0.2Angstrom；theta选取的范围为**[30°,180°]，间隔为5°**。同时==考虑到交换对称性，保证R2>=R1，避免重复计算==
三个原子确定一个平面（代码中取的是yz平面），因此令原子位于：
Si(10,10,10)，O1(10,10-R1,10)，O2(10,10+R1cos(180-θ),10-R1sin(180-θ))
示意图如下：
代码如下：
得到文件：
由于服务器一个节点只能==一次提交150~200个任务==，故要分批提交任务。利用下述脚本实现参数传递，分批提交任务。
下述命令提交R1=0.6，R2={0.6,0.8,1.0,1.2,1.4}的所有构型计算任务 5*[(180-30)/5]=150个任务
同时计算过程中由于文件过大，经询问同学后知可==设置L.WAVE=.FALSE，LCHARGE=.FALSE==，因为只是单点计算，不需要保留波函数和电荷密度信息。
raw文件需要==原子种类、模拟的盒子、原子坐标和原子受力、系统能量和维里应力==的信息（维里应力针对固体体系，本体系不需要）
对应需要准备五个文件：
box.raw, coord.raw, force.raw, energy.raw and virial.raw
以力文件为例：（==上述文件的行数与frame数一致==）
$ cat force.raw
-0.724 2.039 -0.951 0.841 -0.464 0.363
6.737 1.554 -5.587 -2.803 0.062 2.222
-1.968 -0.163 1.020 -0.225 -0.789 0.343
#一行代表一个frame(即一个构型)，例子中有2个原子，每个原子3个方向受力，共6列
#其它raw软件同理，coord.raw(有frame行，atom*3列) box.raw(frame行，盒子三个方向的张量共9列) virial.raw(frame行，维里应力三个方向的张量9列)，energy.raw(frame行，1列)还需要原子种类文件：type.raw（假设每个frame原子种类一样，故==仅需要1行==）
$ cat type.raw
0 1
#一行内原子种类用整数依次表示（保证与原子坐标准备文件中的原子顺序一致）如Si02即为：0 1 1 由于deepmd提供的用于生成raw文件的dpdata程序我看着有点乱，我就直接照它的源码按照文件需要的格式修改了代码，用python提取每个构型的所有的相关数据。
代码如下：
# -*- coding: utf-8 -*-
"""
@author: haojie
"""
import numpy as np
def EXTRACT_POSCAR(file_path,cartesian = True):
'''
提取POSCAR里面的盒子坐标、原子坐标、原子种类信息，返回system用于main调用得到raw文件
'''
f = open(file_path,mode="r")
content = f.readlines()
system = {}
system['atom_names'] = [str(x) for x in content[5].split()]
system['atom_numbs'] = [int(x) for x in content[6].split()]
scale = float(content[1])
box = []
for i in range(2,5):
boxv = [float(ii) for ii in content[i].split()]
boxv = np.array(boxv)*scale
box.append(boxv)
system['box'] = [np.array(box)]
natoms = sum(system['atom_numbs'])
system['natoms'] = natoms
coord = []
for i in range(8,8+natoms):
coordv = [float(ii) for ii in content[i].split()]
coordv = np.array(coordv)*scale
coord.append(coordv)
system['coord'] = [np.array(coord)]
return system
def EXTRACT_OUTCAR(file_path):
'''
提取OUTCAR里面的能量和受力信息，返回Property用于main调用得到raw文件
'''
f = open(file_path,mode="r")
count = 0
for line in f.readlines():
count +=1
f.close()
Property = {}
search_string1 = "TOTAL-FORCE"
f = open(file_path,mode="r")
num = 0
flag = False
force_xyz = []
for line in f.readlines():
if search_string1 in line:
flag = True
if flag and num < 6:
num +=1
if num >2:
force = [str(x) for x in line.split()][3:]
force = [float(x) for x in force]
force = np.array(force)
force_xyz.append(force)
f.close()
Property['force'] = [force_xyz]
search_string2 = "entropy="
f = open(file_path,mode="r")
energy_whole = []
for line in f.readlines():
if search_string2 in line:
energy = line.split()[-1:]
energy_whole.append(energy)
f.close()
Property['energy']=[energy_whole]
return Property
def main():
################################################################################
#POSCAR文件得到box.raw coord.raw type.raw
file_path1 = "C://Users//haoji//Desktop//POSCAR"
System = EXTRACT_POSCAR(file_path1)
print(System)
#box.raw
box_raw = open("box.raw",mode="a")
for i in range(0,3):
axis = [str(x) for x in System['box'][0][i]]
axis = axis[0] + " " + axis[1] + " "+ axis[2] + " "
box_raw.write(axis)
box_raw.write("\n")
box_raw.close()
#coord.raw
coord_raw = open("coord.raw",mode="a")
for i in range(0,System['natoms']):
coord_per = [str(x) for x in System['coord'][0][i]]
coord_per = coord_per[0] + " " + coord_per[1] + " " + coord_per[2] + " "
coord_raw.write(coord_per)
coord_raw.write('\n')
coord_raw.close()
#type_raw
type_raw = open("type.raw",mode="a")
for i in range(0,len(System['atom_numbs'])):
type_per = str(i)
type_per = (type_per + " ")*System['atom_numbs'][i]
type_raw.write(type_per)
type_raw.close()
################################################################################
#OUTCAR文件得到force.raw energy.raw
file_path2 = "C://Users//haoji//Desktop//OUTCAR"
Property = EXTRACT_OUTCAR(file_path2)
print(Property)
#force.raw
force_raw = open("force.raw",mode="a")
for i in range(0,3):
force = [str(x) for x in Property['force'][0][i]]
force = force[0] + " " + force[1]+ " " + force[2] + " "
force_raw.write(force)
force_raw.write("\n")
force_raw.close()
#energy.raw
energy_raw = open("energy.raw",mode="a")
energy = [str(x) for x in Property['energy'][0][0]]
energy = energy[0]
energy_raw.write(energy)
energy_raw.write("\n")
energy_raw.close()
if __name__ == "__main__":
main()用以下命令将每个构型的数据整合到一起：（以box.raw为例）
得到文件：
（1）box.raw
（2）coord.raw
（3）force.raw
（4）energy.raw
（5）type.raw
训练需要参数文件指定训练信息（json格式）：（相当于lammps的in文件）
四个部分（==model、learning_rate、loss、training==）
以水分子为例：
"model": {
"type_map": ["O", "H"],
"descriptor" :{ #设置描述符D 感觉描述符相当于对数据的前处理？使得数据处理后满足各种对称性
"type": "se_a", #smooth-edition angular infomation
"rcut_smth": 5.80, #周围原子对中心原子的贡献从5.80开始降为0
"rcut": 6.00, #周围原子对中心原子的贡献到6.00降为0 (截断半径)
"sel": [46, 92], #截断半径内所能容纳的最大的原子数（46个氧原子，92个氢原子）
"neuron": [25, 50, 100], #神经网络大小
"axis_neuron": 16, #embedding matrix的宽度
"resnet_dt": false, #描述符的残参网络设置 根据经验设置
"seed": 1, #随机数种子
"_comment": " that's all"
},
"fitting_net" : {
"neuron": [240, 240, 240], #拟合的神经网络
"resnet_dt": true,
"seed": 1,
"_comment": " that's all"
},
"_comment": " that's all"
}
"loss" : {
"start_pref_e": 0.02, #开始能量误差的权重系数
"limit_pref_e": 1, #最终能量误差的权重系数
"start_pref_f": 1000, #开始力的权重系数 一开始几乎全部以力的误差作为损失函数降低的贡献
"limit_pref_f": 1,
"start_pref_v": 0, #维里应力 这里简单不考虑维里应力，计算固体考虑弹性常数要打开
"limit_pref_v": 0,
"_comment": " that's all"
}
"learning_rate" :{
"type": "exp", #指数方式衰减
"start_lr": 0.005, #一开始的lr
"decay_steps": 5000, #每5000步衰减一次
"decay_rate": 0.95, #衰减到原来的95% 新的版本用的是stop_lr,系统会自己计算衰减率
"_comment": "that's all"
}
"training" : {
"systems": ["../data/"], #训练文件夹的位置
"set_prefix": "set", #
"stop_batch": 1000000, #训练步数
"batch_size": 1, #同时训练的数量
"seed": 1,
"_comment": " display and restart",
"_comment": " frequencies counted in batch",
"disp_file": "lcurve.out", #指定文件记录损失函数变化
"disp_freq": 100, #隔多少步输出损失函数变化
"numb_test": 10, #测试数量
"save_freq": 1000, #每1000步输出临时文件，以防止意外结束从头计算
"save_ckpt": "model.ckpt", #
"load_ckpt": "model.ckpt",
"disp_training":true,
"time_training":true,
"profiling": false,
"profiling_file":"timeline.json",
"_comment": "that's all"
}下面为官网的训练水分子的流程：（==SiO2对照着就行==）
（1）分割训练集和测试集（下面的为例子的数据6000个frame分割成3组，最后一组作测试集）
利用raw_to_set.sh分割
$ ls
box.raw coord.raw energy.raw force.raw type.raw virial.raw
$ $deepmd_source_dir/data/raw/raw_to_set.sh 2000
nframe is 6000
nline per set is 2000
will make 3 sets
making set 0 ...
making set 1 ...
making set 2 ...
$ ls
box.raw coord.raw energy.raw force.raw set.000 set.001 set.002 type.raw virial.raw
本体系共2418个构型，相应的==每个set 806个frame==
（2）模型训练
$ cd $deepmd_source_dir/examples/water/train/
$ dp train water_se_a.json #执行参数文件
SiO2体系的json文件如下：（基本没改，照水分子那个改的）
{
"_comment": "SiO2",
"model": {
"type_map": ["Si", "O"],
"descriptor" :{
"type": "se_a",
"sel": [46, 92],
"rcut_smth": 5.80,
"rcut": 6.00,
"neuron": [25, 50, 100],
"resnet_dt": false,
"axis_neuron": 16,
"seed": 1,
"_comment": " that's all"
},
"fitting_net" : {
"neuron": [240, 240, 240],
"resnet_dt": true,
"seed": 1,
"_comment": " that's all"
},
"_comment": " that's all"
},
"learning_rate" :{
"type": "exp",
"decay_steps": 5000,
"start_lr": 0.001,
"stop_lr": 3.51e-8,
"_comment": "that's all"
},
"loss" :{
"start_pref_e": 0.02,
"limit_pref_e": 1,
"start_pref_f": 1000,
"limit_pref_f": 1,
"start_pref_v": 0,
"limit_pref_v": 0,
"_comment": " that's all"
},
"_comment": " traing controls",
"training" : {
"systems": ["./"],
"set_prefix": "set",
"stop_batch": 1000000,
"batch_size": 1,
"seed": 1,
"_comment": " display and restart",
"_comment": " frequencies counted in batch",
"disp_file": "lcurve.out",
"disp_freq": 100,
"numb_test": 10,
"save_freq": 1000,
"save_ckpt": "model.ckpt",
"load_ckpt": "model.ckpt",
"disp_training":true,
"time_training":true,
"profiling": false,
"profiling_file":"timeline.json",
"_comment": "that's all"
},
"_comment": "that's all"
}（3）训练结果文件
其中有个lcurve.out文件记录损失函数变化（2~7列分别为总测试误差、总训练误差、能量测试误差、能量训练误差、力测试误差、力训练误差，8列为学习率变化）
# batch l2_tst l2_trn l2_e_tst l2_e_trn l2_f_tst l2_f_trn lr
0 2.67e+01 2.57e+01 2.21e-01 2.22e-01 8.44e-01 8.12e-01 1.0e-03
100 6.14e+00 5.40e+00 3.01e-01 2.99e-01 1.93e-01 1.70e-01 1.0e-03
200 5.02e+00 4.49e+00 1.53e-01 1.53e-01 1.58e-01 1.42e-01 1.0e-03
300 4.36e+00 3.71e+00 7.32e-02 7.27e-02 1.38e-01 1.17e-01 1.0e-03
400 4.04e+00 3.29e+00 3.16e-02 3.22e-02 1.28e-01 1.04e-01 1.0e-03
SiO2体系==训练初期的==lcurve.out：
SiO2体系==训练中期的==lcurve.out：（每20W步列了一个）
SiO2体系==训练结束的==lcurve.out：
精度明显不够，视频中水分子的训练结果精度达到了10-4~10-5，且训练集和测试集的误差一般要在一个量级。
整个训练过程的精度只在10-1~10-2附近波动
- 进一步学习DeepMD软件相关参数的设置，测试并分析计算失败的原因、
- 用当前拟合出来的模型，利用提供的python接口预测一些新构型的能量值，然后再用VASP跑对应的构型，对比一下测试的精度偏差了多少（相当于手动算误差可视化一下）
- SRTP校级答辩需要准备一下答辩材料
- 我对神经网络还不是很懂，可能需要点时间学习，但是没有机器学习的知识不知道能不能上手深度学习。