空间站广场

论文

Notebooks

比赛

课程

Apps

我的主页

我的Notebooks

我的论文库

我的足迹

我的工作空间

任务

节点

文件

数据集

镜像

项目

数据库

公开

DPA-2与化学反应过渡态搜索

过渡态

DeePMD-kit

ASE

过渡态DeePMD-kitASETS

李博文

发布于 2024-02-24

推荐镜像 :DeePMD-kit:pytorch-dpa-2.0

推荐机型 :c2_m4_cpu

数据集

DPA2_model(v2)

©️ Copyright 2024 @ Authors
作者：李博文📨
日期：2024-03-01
共享协议：本作品采用知识共享署名-非商业性使用-相同方式共享 4.0 国际许可协议进行许可。
快速开始：点击上方的开始连接按钮，选择 DeePMD-kit:pytorch-dpa-2.0 镜像及 c2_m4_cpu 节点配置，选取DPA2_model(v2) 数据集，稍等片刻即可运行。

代码

文本

在完成这个教程之后，你将能够：

学习过渡态理论与过渡态搜索中的基础知识
学会通过ASE调用DPA-2模型，通过CI-NEB方法进行过渡态搜索

代码

文本

1. 过渡态与过渡态搜索

过渡态搜索这一部分的理论已经有小伙伴讲解的非常详细，这边就不重复造轮子了~感兴趣的同学可以参考：ATST-Tools | ABACUS，ASE与过渡态理论的碰撞

代码

文本

关于DPA-2模型，主要的介绍在这篇文章： Notebook for DPA-2: Towards a universal large atomic model for molecular and material simulation

代码

文本

2. 通过ASE调用DPA-2模型进行过渡态搜索

代码

文本

2.1 导入模型与数据

代码

文本

首先将模型和反应数据复制到我们所在的目录下。

代码

文本

[7]

%%bash

# 进入root目录

cd /root/

# 将模型和数据复制过来

cp -r /bohr/dpa2model-msob/v2/* ./

# 查看文件夹内容

tree -L 2

.
├── model
│   └── model.pt
├── rxn0
│   ├── neb.traj
│   ├── p.xyz
│   ├── r.xyz
│   └── ts.xyz
├── rxn1
│   ├── p.xyz
│   ├── r.xyz

│   └── ts.xyz
└── rxn2
    ├── p.xyz
    ├── r.xyz
    └── ts.xyz

4 directories, 11 files

代码

文本

可以看到有一个训练好的模型和若干实例反应，接下来先安装一些后面数据处理所需的库。

代码

文本

[8]

# 安装数据分析所需的库

! pip install py3Dmol dpdata ase rmsd

Requirement already satisfied: py3Dmol in /opt/mamba/lib/python3.10/site-packages (2.0.4)
Requirement already satisfied: dpdata in /opt/mamba/lib/python3.10/site-packages (0.2.17)
Requirement already satisfied: ase in /opt/mamba/lib/python3.10/site-packages (3.22.1)
Requirement already satisfied: rmsd in /opt/mamba/lib/python3.10/site-packages (1.5.1)
Requirement already satisfied: h5py in /opt/mamba/lib/python3.10/site-packages (from dpdata) (3.10.0)
Requirement already satisfied: monty in /opt/mamba/lib/python3.10/site-packages (from dpdata) (2024.2.26)
Requirement already satisfied: wcmatch in /opt/mamba/lib/python3.10/site-packages (from dpdata) (8.5)
Requirement already satisfied: numpy>=1.14.3 in /opt/mamba/lib/python3.10/site-packages (from dpdata) (1.26.2)
Requirement already satisfied: scipy in /opt/mamba/lib/python3.10/site-packages (from dpdata) (1.11.4)
Requirement already satisfied: matplotlib>=3.1.0 in /opt/mamba/lib/python3.10/site-packages (from ase) (3.4.3)
Requirement already satisfied: pyparsing>=2.2.1 in /opt/mamba/lib/python3.10/site-packages (from matplotlib>=3.1.0->ase) (3.1.1)
Requirement already satisfied: cycler>=0.10 in /opt/mamba/lib/python3.10/site-packages (from matplotlib>=3.1.0->ase) (0.12.1)
Requirement already satisfied: kiwisolver>=1.0.1 in /opt/mamba/lib/python3.10/site-packages (from matplotlib>=3.1.0->ase) (1.4.5)
Requirement already satisfied: python-dateutil>=2.7 in /opt/mamba/lib/python3.10/site-packages (from matplotlib>=3.1.0->ase) (2.8.2)
Requirement already satisfied: pillow>=6.2.0 in /opt/mamba/lib/python3.10/site-packages (from matplotlib>=3.1.0->ase) (10.1.0)
Requirement already satisfied: bracex>=2.1.1 in /opt/mamba/lib/python3.10/site-packages (from wcmatch->dpdata) (2.4)
Requirement already satisfied: six>=1.5 in /opt/mamba/lib/python3.10/site-packages (from python-dateutil>=2.7->matplotlib>=3.1.0->ase) (1.16.0)
WARNING: Running pip as the 'root' user can result in broken permissions and conflicting behaviour with the system package manager. It is recommended to use a virtual environment instead: https://pip.pypa.io/warnings/venv

代码

文本

导入训练好的DPA-2模型，模型路径为：/root/model/model.pt

代码

文本

[9]

# 导入DPA-2的ASE接口

from deepmd_pt.utils.ase_calc import DPCalculator

# 读取模型，创建DPA-2计算器实例

calc = DPCalculator(model = '/root/model/model.pt')

代码

文本

2.2 计算分子的能量与力

代码

文本

用ase读取一个xyz文件，并且尝试计算它的能量和力，注意ase默认输出的单位是电子伏特（ev）和埃（Å）。

代码

文本

[10]

# 导入ase的读取和写入模块

from ase.io import read,write

# 随便读取一个xyz文件，这里可以修改为任何已有的xyz文件，如果你的个人文件夹（/personal）中有结构也可以尝试读取

data = read('root/rxn0/r.xyz')

# DPA-2模型计算能量和力需要cell信息，这里设置的值对能量计算无影响

data.set_cell([10,10,10])

# 为读取的xyz结构设置一个计算器

data.calc = calc

# 分别打印该结构的能量和力信息

print(f"energy：{data.get_potential_energy()} (ev)")

print(f"forces：{data.get_forces()} (ev/Å)")

energy：-8636.717467628849 (ev)
forces：[[-1.16979954e-02 -8.11188200e-03 -9.46256525e-05]
 [ 3.64224220e-02  4.81465264e-02  8.15033011e-03]
 [ 6.83969372e-03  4.54776803e-02 -3.78356003e-04]
 [-7.19191241e-02 -1.76068434e-02  1.63492994e-03]
 [-1.95175025e-02 -3.02620462e-02 -8.43866504e-03]
 [-3.32095354e-02 -2.73219947e-02  1.32785867e-03]
 [ 5.96679969e-02 -1.22404825e-02 -2.29667202e-03]
 [ 3.29042145e-02 -3.23853076e-02 -8.74523257e-03]
 [ 5.09830203e-04  3.43043498e-02  8.84043258e-03]] (ev/Å)

代码

文本

2.3 进行过渡态的搜索

代码

文本

读取一个反应的反应物和产物，通过CI-NEB方法计算路径上的鞍点，这个点就是我们想要的过渡态。（这一步可能需要1-5分钟，请耐心等待）

代码

文本

[11]

#导入ase中的NEB模块，LBFGS优化器和FIRE优化器

from ase.neb import NEB

from ase.optimize import LBFGS,FIRE

# 选择你想要计算的反应（rxn0/rxn1/rxn2均可）

rxn_name='/root/rxn0'

initial_structure = read(f'{rxn_name}/r.xyz') # 初始结构

final_structure = read(f'{rxn_name}/p.xyz') # 最终结构

images = [initial_structure]

images += [initial_structure.copy() for _ in range(10)] # 在中间插入10个图像

images += [final_structure]

# 设置NEB计算所需参数，这里采用攀爬图像方法，用于寻找反应路径上的最高点，

neb = NEB(images, climb=True, method='improvedtangent', allow_shared_calculator=True)

# 采用IDPP插值对中间的图像进行计算，生成初始路径

neb.interpolate()

# 对每个图像设置DPA-2计算器

for image in images:

image.set_cell([10,10,10])

image.set_calculator(calc)

# 使用LBFGS优化器进行优化，轨迹保存为neb.traj

optimizer = LBFGS(neb, trajectory=f'{rxn_name}/neb.traj')

#开始运行

optimizer.run(fmax=0.05)

       Step     Time          Energy         fmax
*Force-consistent energies used in optimization.
LBFGS:    0 13:19:57    -8633.801599*       5.6959
LBFGS:    1 13:20:08    -8634.239104*       4.0642
LBFGS:    2 13:20:20    -8634.516158*       3.4608
LBFGS:    3 13:20:31    -8634.592346*       1.5398
LBFGS:    4 13:20:42    -8634.639983*       1.2405
LBFGS:    5 13:20:53    -8634.662912*       0.7734
LBFGS:    6 13:21:04    -8634.678395*       0.6354
LBFGS:    7 13:21:15    -8634.692497*       0.9626
LBFGS:    8 13:21:26    -8634.694729*       0.8859
LBFGS:    9 13:21:38    -8634.693292*       0.5709
LBFGS:   10 13:21:49    -8634.695464*       0.5142
LBFGS:   11 13:22:00    -8634.696080*       0.3178
LBFGS:   12 13:22:11    -8634.695762*       0.3383
LBFGS:   13 13:22:22    -8634.696360*       0.2134
LBFGS:   14 13:22:33    -8634.696628*       0.1762
LBFGS:   15 13:22:45    -8634.696702*       0.1438
LBFGS:   16 13:22:56    -8634.696782*       0.1188
LBFGS:   17 13:23:07    -8634.696765*       0.0643
LBFGS:   18 13:23:18    -8634.696781*       0.0520
LBFGS:   19 13:23:29    -8634.696781*       0.0505
LBFGS:   20 13:23:40    -8634.696771*       0.0363

True

代码

文本

绘制能量变化的折线图，将路径上能量最高点的结构保存为xyz文件。

代码

文本

[12]

# 导入绘图的第三方库

import matplotlib.pyplot as plt

# 导入kcal和mol单位，用于单位转换

from ase.units import kcal,mol

# 计算最终路径上每个点的能量

energies = [image.get_potential_energy() for image in images]

# 绘制能量变化折线图

energies_kcal = [(i-energies[0])*mol/kcal for i in energies] #将能量单位从ew转化为kcal/mol，并且以第一个结构为0点

plt.plot(energies_kcal)

plt.xlabel('Image Number') # 设置x轴标题

plt.ylabel('Energy (kcal/mol)') # 设置y轴标题

plt.show() # 展示图像

# 将能量最高点以及整条路径保存为xyz文件

write(f'{rxn_name}/neb_ci.xyz', images[energies.index(max(energies))])

write(f'{rxn_name}/neb_path.xyz', images)

代码

文本

2.4 分析NEB路径最高点与真实过渡态的差异

代码

文本

看下标准过渡态结构和NEB路径上的能量最高点结构是否一致 (按住鼠标左键拖动，可以让分子在空间中旋转)

代码

文本

[13]

import py3Dmol # 导入绘图所需模块

# 分别画出标准过渡态结构以及路径上能量最高点结构

for xyz in ['ts.xyz','neb_ci.xyz']:

view = py3Dmol.view(width=300, height=300)

# 向视图中添加XYZ格式的分子，并设置显示样式为球棍模型

print(f'{rxn_name}/{xyz}')

view.addModel(open(f'{rxn_name}/{xyz}').read(), "xyz")

view.setStyle({'stick':{'radius': 0.1}, 'sphere':{'radius': 0.4}}) # 设置球棍模型参数

view.zoomTo() # 调整视图以适应分子模型

# 显示分子模型

view.show()

/root/rxn0/ts.xyz

/root/rxn0/neb_ci.xyz

代码

文本

计算NEB最高点的结构与标准过渡态的几何误差（单位：埃），对于小分子而言RMSD小于0.1就已经非常相似了。

代码

文本

[14]

import os

# 计算两个结构之间的RMSD

rmsd = os.system(f'calculate_rmsd {rxn_name}/ts.xyz {rxn_name}/neb_ci.xyz')

0.00417101261925496

代码

文本

查看NEB计算得到的反应路径，感受化学反应的变化过程。

代码

文本

[15]

def read_xyz_file(filepath):

"""读取包含多个XYZ结构的文件，并返回一个包含所有结构的列表"""

structures = []

with open(filepath, 'r') as file:

content = file.read()

# 以空行分割不同的XYZ结构

for block in content.strip().split('\n\n'):

structures.append(block)

return structures

# 假设已经通过之前的read_xyz_file函数读取了XYZ结构到xyz_structures列表中

xyz_structures = read_xyz_file(f'{rxn_name}/neb_path.xyz')

# 创建一个py3Dmol视图

view = py3Dmol.view(width=400, height=400)

# 将所有结构作为帧添加到视图中

for xyz in xyz_structures:

view.addModelsAsFrames(xyz, "xyz")

# 设置球棍模型样式

view.setStyle({'model': -1}, {'stick':{'radius': 0.1}, 'sphere':{'radius': 0.4}}) # '-1'表示对所有模型应用这个样式

# 设置动画参数

view.animate({'loop': "forward", 'reps': 0,'interval': 500}) # 'reps': 0 表示无限循环

# 调整视图以适应所有分子模型

view.zoomTo()

# 显示分子模型动画

view.show()

代码

文本

计算得到的路径与标准路径的能垒误差

rxn0/1/2的实际能垒：47.35/46.78/45.68kcal/mol

代码

文本

[16]

# 输出能量最高点与反应物的能量差

print(f"{rxn_name}的预测能垒：{round((max(energies_kcal)-energies_kcal[0]),3)} (kcal/mol)")

/root/rxn0的预测能垒：46.598 (kcal/mol)

代码

文本

过渡态

DeePMD-kit

ASE

过渡态DeePMD-kitASETS

已赞12

本文被以下合集收录

DPA-2

yuxiangc22

更新于 2024-09-03

5 篇5 人关注

测试合集文章列表100篇

xingyanshi@dp.tech很长的名字xingyanshi@dp.tech很长的名字很长的

更新于 2024-08-04

104 篇2 人关注