LiPF6晶体在EC/DMC混合溶液中的溶解 | Bohrium-玻尔科研空间站

空间站广场

论文

Notebooks

比赛

课程

Apps

我的主页

我的Notebooks

我的论文库

我的足迹

我的工作空间

任务

节点

文件

数据集

镜像

项目

数据库

公开

LiPF6晶体在EC/DMC混合溶液中的溶解

分子动力学

镜影映照

Piloteye

更新于 2024-06-19

推荐镜像 :GROMACS:2022.2

推荐机型 :c3_m4_1 * NVIDIA T4

数据集

模拟文件(v2)

©️ Copyright 2023 @ Authors
作者： 镜影映照 📨
日期：2023-12-09
快速开始：点击上方的 开始连接 按钮，选择 registry.dp.tech/dptech/gromacs:2022.2镜像 和任意配置机型即可开始。

代码

文本

预置文件介绍

代码

文本

为了能顺利进行模拟，这里准备了一些文件，展示如下

代码

文本

[ ]

dataset_path = "/bohr/train-juyu/v2/"

!tree $dataset_path

代码

文本

这些文件可以分为四类，

pdb和cif文件是要模拟的体系中各分子的结构文件，定义了它们的坐标。packmol则是用于生成初始盒子的程序。
topol下的文件包含了体系的拓扑信息，定义了体系中各分子的连接方式，同时还包含了力场参数，也即原子间相互作用的信息
mdp下的文件则用于控制模拟流程，定义了模拟的参数，比如模拟的步长，温度等
对于CPU机型模拟可能较慢，因此还准备了已经模拟完成的轨迹文件md_whole.trr 和npt2.gro，用于展示溶解过程。

代码

文本

安装和导入python包，定义辅助函数

代码

文本

[ ]

# 安装必要的库

!pip install ase

!pip install mdtraj==1.9.6

!pip install nglview==3.0.3

!pip install ipywidgets==7.6.5

代码

文本

[ ]

# 导入必要的python包

from ase.io import read,write

import mdtraj as md

from mdtraj import compute_neighbors

import numpy as np

import subprocess

import nglview as nv

代码

文本

[ ]

999

100

101

102

103

104

105

106

107

108

109

110

111

112

113

114

115

116

117

118

119

120

121

122

123

124

125

126

127

128

129

130

131

132

133

134

135

136

137

138

139

140

141

142

143

144

145

146

147

148

149

150

151

152

153

154

155

156

157

158

159

160

161

162

163

164

165

166

167

168

169

170

171

172

173

174

175

176

177

178

179

180

181

›

⌄

# 定义辅助函数

def rearrange_atoms(old_gro, new_gro):

traj = md.load(old_gro)

new_index = []

Li_index = traj.topology.select("name Li")

new_index += Li_index.tolist()

P_index = traj.topology.select("name P")

for i in range(len(P_index)):

F_index = compute_neighbors(traj, 0.17, [P_index[i]])

if len(F_index[0]) != 6:

raise ValueError("P atom is not surrounded by 6 F atoms, maybe the cutoff is too small")

new_index.append(P_index[i])

new_index += F_index[0].tolist()

new_atom_names = [traj.topology.atom(i).name for i in new_index]

new_coordinates = traj.xyz[:, new_index, :]

traj.xyz = new_coordinates

for i, atom in enumerate(traj.topology.atoms):

atom.name = new_atom_names[i]

traj.save_gro(new_gro)

def gen_unwrap_top(filename,dataset_path=dataset_path):

traj = md.load(filename)

n_LIPF6 = len(traj.topology.select("element P"))

unwrap_top_str =f"""; Include forcefield parameters

#include "amber99sb.ff/forcefield.itp"

#include "{dataset_path}/topol/Atomtypes.itp"

#include "amber99sb.ff/ions.itp"

#include "{dataset_path}/topol/PF6.itp"

[ system ]

GAFF force field

[ molecules ]

; Compound nmols

LI {n_LIPF6}

PF6 {n_LIPF6}

"""

with open("unwrap.top", "w") as f:

f.write(unwrap_top_str)

def find_closest_pairs(array1, array2):

min_diffs = [np.inf, np.inf]

closest_pairs = [None, None]

for i in range(len(array2)):

for j in range(i+1, len(array2)):

avg = (array2[i] + array2[j]) / 2

diff = np.linalg.norm(avg - array1)

if diff < max(min_diffs):

index = min_diffs.index(max(min_diffs))

min_diffs[index] = diff

closest_pairs[index] = [i, j]

if max(min_diffs) == 0:

break

if max(min_diffs) == 0:

break

if max(min_diffs) > 0.1:

raise ValueError("Something went wrong, the closest pairs are too far apart")

remaining_indices = list(set(range(len(array2))) - set(np.array(closest_pairs).flatten()))

closest_pairs.append(remaining_indices)

return closest_pairs

def ajust_F_order(old_gro, new_gro):

traj = md.load(old_gro)

new_index = []

Li_index = traj.topology.select("name LI")

new_index += Li_index.tolist()

P_index = traj.topology.select("name P")

for i in P_index:

new_index.append(i)

F_index = np.linspace(i+1, i+6, 6, dtype=int)

P_atom = traj.xyz[0, i, :]

F_atoms = traj.xyz[0, F_index, :]

closest_pairs = find_closest_pairs(P_atom, F_atoms)

F_array = closest_pairs[0]

for pair in closest_pairs[1:]:

mid = len(F_array) // 2

F_array = F_array[:mid] + pair + F_array[mid:]

if len(F_array) != 6:

raise ValueError("Something went wrong, the new array is not of length 6")

new_index += F_index[F_array].tolist()

traj.xyz = traj.xyz[:, new_index, :]

traj.save_pdb(new_gro)

def get_system_info(filename):

traj = md.load(filename)

com = md.compute_center_of_mass(traj)[0]

n_LIPF6 = len(traj.topology.select("element P"))

box_len = np.array([1.6*(n_LIPF6)**(1/3)]*3)

n_EC = 7*n_LIPF6

n_DMC = 7*n_LIPF6

LIPF6_pos = box_len/2 - com

return box_len, n_LIPF6, n_EC, n_DMC, LIPF6_pos

def gen_packmol_input(filename,dataset_path=dataset_path):

box_len, n_LIPF6, n_EC, n_DMC, LIPF6_pos = get_system_info(filename)

box_len_str = " ".join([str(i) for i in box_len*10])

LIPF6_pos_str = " ".join([str(i) for i in LIPF6_pos*10])

packmol_input =f"""tolerance 2.0

filetype pdb

output system.pdb

structure {filename}

number 1

fixed {LIPF6_pos_str} 0.0 0.0 0.0

end structure

structure {dataset_path}/EC.pdb

number {n_EC}

inside box 0. 0. 0. {box_len_str} # Box dimensions

end structure

structure {dataset_path}/DMC.pdb

number {n_DMC}

inside box 0. 0. 0. {box_len_str} # Box dimensions

end structure"""

with open("packmol.inp", "w") as f:

f.write(packmol_input)

return packmol_input

def gen_md_top(filename,dataset_path=dataset_path):

box_len, n_LIPF6, n_EC, n_DMC, LIPF6_pos = get_system_info(filename)

top_str =f"""; Include forcefield parameters

#include "amber99sb.ff/forcefield.itp"

#include "{dataset_path}/topol/Atomtypes.itp"

#include "amber99sb.ff/ions.itp"

#include "{dataset_path}/topol/PF6.itp"

#include "{dataset_path}/topol/EC.itp"

#include "{dataset_path}/topol/DMC.itp"

[ system ]

GAFF force field

[ molecules ]

; Compound nmols

LI {n_LIPF6}

PF6 {n_LIPF6}

EC {n_EC}

DMC {n_DMC}

"""

with open("topol.top", "w") as f:

f.write(top_str)

return top_str

代码

文本

准备输入体系的结构文件

代码

文本

结构文件的准备可分为以下几个步骤

首先需要下载LiPF6的晶体文件；
然后将晶体文件扩胞为超胞；
这之后，需要对超胞的结构进行一些调整；
最后再将调整后的LiPF6放置于EC/DMC的混合溶液中。

代码

文本

准备LiPF6晶体

代码

文本

首先从材料数据库下载LiPF6的晶体cif格式文件

代码

文本

[ ]

# 查看LiPF6的结构

view = nv.show_ase(read(f"{dataset_path}/LiPF6.cif"))

view.clear_representations()

view.add_representation("ball+stick", selection=range(0,3),color = "white")

view.add_representation("ball+stick", selection=range(3,24))

view

代码

文本

cif文件只包含了LiPF6的一个单胞，有3个LiPF6。将鼠标移动到图中原子上，可以看见白色的是Li，黄色的为P，绿色的为F。由于晶体的周期性，边缘的PF6是不完整的。

为了得到较大的晶体，首先需要对其进行扩胞。这里在三个方向上，分别让晶体扩展2、4、1倍。

代码

文本

[ ]

# 指定x,y,z方向上的扩胞倍数

cell_repaet = [2,4,1]

# 进行扩胞，并将扩胞后的结构写入gro文件

bulk_rutile_LIPF6 = read(f"{dataset_path}/LiPF6.cif")

lattice_LIPF6 = bulk_rutile_LIPF6.repeat(cell_repaet)

write("LiPF6_supercell.gro",lattice_LIPF6)

# 查看超胞的结构

traj = md.load("LiPF6_supercell.gro")

view = nv.show_mdtraj(traj)

view.clear_representations()

view.add_representation("ball+stick", selection=traj.topology.select("name Li"),color = "white")

view.add_representation("ball+stick", selection=traj.topology.select("name P F"))

view

代码

文本

如此一来，超胞的结构便生成好了。但是超胞边缘PF6的结构由于周期性并不完整，因此还需要对结构文件做一些调整。这里我们使用gromacs将PF6的结构变得完整。

代码

文本

[ ]

# 重新排列原子顺序，使得LiPF6的结构符合GROMACS的要求

rearrange_atoms("LiPF6_supercell.gro", "LiPF6_rearranged.gro")

# 生成unwrap.top文件

gen_unwrap_top("LiPF6_rearranged.gro")

# 生成gromacs输入文件并运行

gen_unwrap_tpr = f"gmx grompp -f {dataset_path}/mdp/unwrap.mdp -p unwrap.top -o unwrap.tpr -c LiPF6_rearranged.gro -maxwarn 1"

gen_unwrap_structure = "echo 0 | gmx trjconv -f LiPF6_rearranged.gro -s unwrap.tpr -pbc whole -o LiPF6_unwrap.gro"

subprocess.run(gen_unwrap_tpr, shell=True)

subprocess.run(gen_unwrap_structure, shell=True)

代码

文本

[ ]

# 查看新的结构

traj = md.load("LiPF6_unwrap.gro")

view = nv.show_mdtraj(traj)

view.clear_representations()

view.add_representation("ball+stick", selection="Li")

view.add_representation("ball+stick", selection="PF6")

view

代码

文本

这样PF6的结构就完整了，最后还需要简单调整每个PF6中F原子的顺序，使得其与力场文件中一致。

代码

文本

[ ]

# 调整F原子的顺序以与top文件中的顺序一致

ajust_F_order("LiPF6_unwrap.gro", "LiPF6_ajusted_F.pdb")

代码

文本

组装各组分生成模拟盒子

代码

文本

调整了LiPF6的结构之后，就可以通过packmol将其放置于EC/DMC的混合溶液中了。

packmol生成的结构没有盒子大小信息，通过gromacs简单添加一下即可。

代码

文本

[ ]

filename = "LiPF6_ajusted_F.pdb"

box_len, n_LIPF6, n_EC, n_DMC, LIPF6_pos = get_system_info(filename)

# 生成模拟盒子

gen_packmol_input(filename)

packmol = f"chmod +x {dataset_path}/packmol &{dataset_path}/packmol < packmol.inp"

subprocess.run(packmol, shell=True)

# 添加模拟盒子大小信息

editconf = f"gmx editconf -f system.pdb -o system.gro -box {box_len[0]} {box_len[1]} {box_len[2]}"

subprocess.run(editconf, shell=True)

代码

文本

[ ]

# 查看体系结构

view = nv.show_mdtraj(md.load("system.gro"))

view.clear_representations()

view.add_representation("ball+stick", selection="Li")

view.add_representation("ball+stick", selection="PF6")

view.add_representation("line", selection="EC DMC",linewidth=0.5)

view

代码

文本

至此，我们已经准备好了结构文件，下面可以进行模拟了。

代码

文本

进行分子动力学模拟

代码

文本

分子动力学模拟的流程可以分为以下几步：

进行能量最小化
进行较短的nvt模拟，初步平衡体系
进行压力为1000 bar的npt模拟
进行压力为1 bar的npt模拟
进行较长的压力为1 bar的npt模拟，模拟LiPF6的溶解

需要注意的是，在2~5中，使用了较高的模拟温度（390K），以在短时间内实现LiPF6的溶解；同时在2~4中，我们对LiPF6的位置进行了约束，确保其结构的稳定。

代码

文本

能量最小化

代码

文本

[ ]

# 生成模拟所需的top文件

gen_md_top("LiPF6_ajusted_F.pdb")

# 生成能量最小化所需的输入文件并运行

gen_em_tpr = f"gmx grompp -f {dataset_path}/mdp/em.mdp -p topol.top -o em.tpr -c system.gro -maxwarn 1"

em_run = "gmx mdrun -v -deffnm em"

subprocess.run(gen_em_tpr, shell=True)

subprocess.run(em_run, shell=True)

代码

文本

nvt模拟

代码

文本

[ ]

# 进行一个较短的nvt模拟，以初步平衡体系

gen_nvt_tpr = f"gmx grompp -f {dataset_path}/mdp/nvt.mdp -p topol.top -o nvt.tpr -c em.gro -r system.gro -maxwarn 1"

nvt_run = "gmx mdrun -v -deffnm nvt"

subprocess.run(gen_nvt_tpr, shell=True)

subprocess.run(nvt_run, shell=True)

代码

文本

npt1 模拟

代码

文本

[ ]

# 在使用packmol创建盒子时，体系较大，因此进行一个压强较大的模拟快速压缩盒子

gen_npt1_tpr = f"gmx grompp -f {dataset_path}/mdp/npt1.mdp -p topol.top -o npt1.tpr -c nvt.gro -r nvt.gro -maxwarn 1"

npt1_run = "gmx mdrun -v -deffnm npt1"

subprocess.run(gen_npt1_tpr, shell=True)

subprocess.run(npt1_run, shell=True)

代码

文本

npt2 模拟

代码

文本

[ ]

# 进行常压npt模拟，以平衡体系

gen_npt2_tpr = f"gmx grompp -f {dataset_path}/mdp/npt2.mdp -p topol.top -o npt2.tpr -c npt1.gro -r npt1.gro -maxwarn 1"

npt2_run = "gmx mdrun -v -deffnm npt2"

subprocess.run(gen_npt2_tpr, shell=True)

subprocess.run(npt2_run, shell=True)

代码

文本

md 模拟

代码

文本

[ ]

# 进行较长的md模拟，观察LiPF6的溶解

gen_md_tpr = f"gmx grompp -f {dataset_path}/mdp/md.mdp -p topol.top -o md.tpr -c npt2.gro -r npt2.gro -maxwarn 1"

md_run = "gmx mdrun -v -deffnm md -rdd 0.5"

subprocess.run(gen_md_tpr, shell=True)

subprocess.run(md_run, shell=True)

# 转换轨迹使分子在跨越周期性边界时保持完整

traj_whole = "echo 0 | gmx trjconv -f md.trr -s md.tpr -pbc whole -o md_whole.trr"

subprocess.run(traj_whole, shell=True)

代码

文本

查看LiPF6的溶解过程

代码

文本

[ ]

import mdtraj as md

import nglview as nv

# 若未进行模拟，则注释第一行，并把第二行的注释去掉

traj = md.load("md_whole.trr", top="npt2.gro")

# traj = md.load(f"{dataset_path}/md_whole.trr", top=f"{dataset_path}/npt2.gro")

view = nv.show_mdtraj(traj)

view.clear_representations()

view.add_representation("ball+stick", selection="Li")

view.add_representation("ball+stick", selection="PF6")

view.add_representation("line", selection="EC DMC",linewidth= 0.5)

view.center()

view.add_unitcell()

view.camera = 'orthographic'

view

代码

文本

分子动力学

已赞3

本文被以下合集收录

bohr61096f

更新于 2024-08-27

47 篇0 人关注