空间站广场

论文

Notebooks

比赛

课程

Apps

我的主页

我的Notebooks

我的论文库

我的足迹

我的工作空间

任务

节点

文件

数据集

镜像

项目

数据库

公开

樊哲勇《分子动力学模拟》python实例 | 2.使用近邻列表的分子动力学模拟程序

樊哲勇《分子动力学模拟》程序代码python实现

mosey

yuxiangc22

更新于 2024-07-04

推荐镜像 :Basic Image:bohrium-notebook:2023-03-26

推荐机型 :c2_m4_cpu

数据集

arMD(v1)

©️ Copyright 2023 @ Authors
作者：ChenYuxiang LiDenan(AI4S1O1-材料小队学习志愿者) 📨
日期：2024-6-16
共享协议：本作品采用知识共享署名-非商业性使用-相同方式共享 4.0 国际许可协议进行许可。
快速开始：点击上方的开始连接按钮，选择 bohrium-notebook:2023-03-26 镜像及 c2_m4_cpu 节点配置即可开始运行。

代码

文本

樊哲勇《分子动力学模拟》python实例 | 2.使用近邻列表的分子动力学模拟程序

樊哲勇老师最近出版的书籍《分子动力学模拟》中，使用C++语言进行分子模拟程序的实现。未了更好地理解分子动力学模拟，AI4S1O1-材料小队学习志愿者决定针对樊老师书籍中的内容，进行python语言的程序实现以及注释说明。

本节内容对应樊老师书中的第三章

樊哲勇《分子动力学模拟》python实例 | 概述

代码

文本

本章目标

了解如何在 三斜盒子 内进行坐标系变换、施加周期性边界条件和最小镜像约定
学习 近邻列表 的构造与更新
用python实现上述内容

代码

文本

一、三斜盒子

1.1 定义

在simpleMD中，我们只考虑了正交盒子，正交盒子只有 $L_{x}$ 、 $L_{y}$ 和 $L_{z}$ 三个边长的自由度，三个边长是相互垂直正交的， $L_{x}$ 沿 $x$ 方向， $L_{y}$ 沿 $y$ 方向， $L_{z}$ 沿 $z$ 方向，因此没有角度上的自由度。本节我们将引入三斜盒子（triclinic box）。下图中左边就是正交盒子，右边是三斜盒子

alt

由上图可知，盒子的三条有向线段可以表示为矢量，记为 $a$ 、 $b$ 和 $c$ 。 $e_{x}$ 、 $e_{y}$ 和、 $e_{z}$ 是一组标准正交基，而模拟盒子的三条边 $a$ 、 $b$ 和 $c$ 可以表示为：

$⎩ ⎨ ⎧ a b c = = = a_{x} e_{x} + a_{y} e_{y} + a_{z} e_{z} b_{x} e_{x} + b_{y} e_{y} + b_{z} e_{z} c_{x} e_{x} + c_{y} e_{y} + c_{z} e_{z}$

上述的系数可以组合成一个矩阵，记为盒子矩阵 $H$ ：

$H = a_{x} a_{y} a_{z} b_{x} b_{y} b_{z} c_{x} c_{y} c_{z}$

上一章讨论的正交盒子可以看成是一种特殊的三斜盒子： $⎩ ⎨ ⎧ a b c = = = L_{x} e_{x} + 0 e_{y} + 0 e_{z} 0 e_{x} + L_{y} e_{y} + 0 e_{z} 0 e_{x} + 0 e_{y} + L_{z} e_{z}$ $H = L_{x} 00 0 L_{y} 0 00 L_{z}$

三斜盒子的体积 $V$ 等于上述矩阵行列式 $det (H)$ 的绝对值： ${V det (H) = ∣ det (H) ∣ = a_{x} (b_{y} c_{z} - c_{y} b_{z}) + b_{x} (c_{y} a_{z} - a_{y} c_{z}) + c_{x} (a_{y} b_{z} - b_{y} a_{z})$

为方便讨论，我们记盒子的矩阵 $H$ 的逆为 $G = H^{- 1}$ 。

对于三斜盒子，三个面的面积 $A_{b c}$ , $A_{a c}$ 和 $A_{ab}$ 为: $⎩ ⎨ ⎧ A_{b c} A_{a c} A_{ab} = = = ∣∣ b \times c ∣∣ ∣∣ a \times c ∣∣ ∣∣ a \times b ∣∣$

则垂直于平面的厚度 $L_{a}$ 、 $L_{b}$ 和 $L_{c}$ 为： $⎩ ⎨ ⎧ L_{a} L_{b} L_{c} = = = V / A_{b c} V / A_{a c} V / A_{ab}$

1.2 盒子里的原子坐标

对于一个盒子内的原子，其坐标 $r = (x, y, z)$ 可以表示为:

$r = s_{a} a + s_{b} b + s_{c} c$

$(s_{a}, s_{b}, s_{c})$ 表示原子在盒子内的相对位置，因为原子只能在盒子内部或其表面，故定义域为 $[0, 1]$ ，称为分数坐标。上式可以表示为如下矩阵形式：

$x y z = a_{x} a_{y} a_{z} b_{x} b_{y} b_{z} c_{x} c_{y} c_{z} s_{a} s_{b} s_{c}$

上式可用 $H$ 简写为:

$x y z = H s_{a} s_{b} s_{c}$

将盒子矩阵的逆 $G$ 作用在上述等式两边，可得：

$s_{a} s_{b} s_{c} = G x y z$

通过上述2个等式，我们便可以很方便的将盒子内的原子坐标在分数坐标和笛卡尔坐标之间相互转化。

1.3 周期性边界条件

在得到原子的分数坐标 $(s_{a}, s_{b}, s_{c})$ 后，可很方便的施加周期性边界条件：

将原坐标 $r$ 转换为分数坐标 $s$ 。
对于分数坐标 $s$ 中任意一个分量 $s_{i}$ ：若 $s_{i} > 1$ ，则 $s_{i} = s_{i} - 1$ ；若 $s_{i} < 0$ ，则 $s_{i} = s_{i} + 1$ 。
将分数坐标 $s$ 转换为笛卡尔坐标 $r$

1.4 最小镜像约定

如果有粒子 $i$ （坐标 $r_{i} = (x_{i}, y_{i}, z_{i})$ ）和粒子 $j$ （坐标 $r_{j} = (x_{j}, y_{j}, z_{j})$ ）,在三斜盒子内计算两个粒子的相对坐标时，也可以很方便的施加最小镜像约定：

计算粒子 $i$ 和粒子 $j$ 的相对坐标： $r_{ij} = (x_{ij}, y_{ij}, z_{ij})$ 。
用盒子逆矩阵 $G$ 可以将 $r_{ij}$ 变为分数坐标 $s_{ij} = (ξ_{ij}, η_{ij}, ζ_{ij})$ ，变换关系为：

$ξ_{ij} η_{ij} ζ_{ij} = G x_{ij} y_{ij} z_{ij}$

实施最小镜像约定：对分数坐标 $s_{ij}$ 中任意一个坐标分量 $s$ ，若 $s < - \frac{1}{2}$ ，则 $s = s + 1$ ；若 $s > \frac{1}{2}$ ，则 $s = s - 1$ 。
将变换后的分数坐标转换为普通的笛卡尔坐标：

$x_{ij} y_{ij} z_{ij} = H ξ_{ij} η_{ij} ζ_{ij}$

1.5 三斜盒子的代码实现

大家可以思考如何用python3实现以下功能：

计算三斜晶胞的体积，边长
分数坐标和笛卡尔坐标相互转化
施加最小镜像约定和周期性边界条件

下面给出了一个参考代码，需要注意的是：代码中的盒子矩阵是上述 $H$ 的转置

代码

文本

[1]

import numpy as np

from ase.io import read

def calculate_volume_length(box: np.ndarray) -> None:

'''

计算晶胞的体积和晶胞的三个晶向的长度

Args:

box: 盒子矩阵

'''

volume = np.abs(np.linalg.det(box))

area_ab = np.linalg.norm(np.cross(box[0], box[1]))

area_bc = np.linalg.norm(np.cross(box[1], box[2]))

area_ac = np.linalg.norm(np.cross(box[2], box[0]))

l_c = volume / area_ab

l_a = volume / area_bc

l_b = volume / area_ac

print(f'volume: {volume}')

print(f'l_a: {l_a}, l_b: {l_b}, l_c: {l_c}')

def coord_conversion(coords: np.ndarray, box: np.ndarray) -> None:

'''

坐标系转换

Args:

coords: 笛卡尔坐标

box: 盒子矩阵

'''

box_inv = np.linalg.inv(box)

# 笛卡尔坐标转为分数坐标

frac_coords = np.dot(coords, box_inv)

# 分数坐标转为笛卡尔坐标

cart_coords = np.dot(frac_coords, box)

print(f'frac_coords: \n{frac_coords}')

print(f'cart_coords: \n{cart_coords}')

def apply_pbc(coords: np.ndarray, box: np.ndarray) -> None:

'''

应用周期性边界条件

Args:

coords: 笛卡尔坐标

box: 盒子矩阵

'''

box_inv = np.linalg.inv(box)

frac_coords = np.dot(coords, box_inv)

frac_coords[frac_coords < 0] += 1

frac_coords[frac_coords >= 1] -= 1

new_coords = np.dot(frac_coords, box)

print(f'origin coords: \n{coords}')

print(f'new coords: \n{new_coords}')

def apply_mic(rij: np.ndarray, box: np.ndarray) -> None:

'''

应用最小镜像约定

Args:

rij: 原子间距

box: 盒子矩阵

'''

box_inv = np.linalg.inv(box)

frac_coords = np.dot(rij, box_inv)

frac_coords[frac_coords < -0.5] += 1

frac_coords[frac_coords >= 0.5] -= 1

new_rij = np.dot(frac_coords, box)

print(f'origin rij: \n{rij}')

print(f'new rij: \n{new_rij}')

def main():

atom = read("/bohr/md-simulation-linearMD-gk7w/v1/dataset/Si.vasp")

box = atom.get_cell()

coords = atom.get_positions()

calculate_volume_length(box)

coord_conversion(coords, box)

apply_pbc(coords, box)

apply_mic(coords[0] - coords[1], box)

main()

volume: 161.31810738968053
l_a: 5.443702372939453, l_b: 5.443702372939453, l_c: 5.443702372939453
frac_coords: 
[[ 7.50000000e-01  7.50000000e-01  2.50000000e-01]
 [-4.45929298e-33  5.00000000e-01  5.00000000e-01]
 [ 7.50000000e-01  2.50000000e-01  7.50000000e-01]
 [ 0.00000000e+00  0.00000000e+00  0.00000000e+00]
 [ 2.50000000e-01  7.50000000e-01  7.50000000e-01]
 [ 5.00000000e-01  5.00000000e-01  2.45444455e-33]
 [ 2.50000000e-01  2.50000000e-01  2.50000000e-01]
 [ 5.00000000e-01  0.00000000e+00  5.00000000e-01]]
cart_coords: 
[[4.08277678e+00 4.08277678e+00 1.36092559e+00]
 [4.50000000e-16 2.72185119e+00 2.72185119e+00]
 [4.08277678e+00 1.36092559e+00 4.08277678e+00]
 [0.00000000e+00 0.00000000e+00 0.00000000e+00]
 [1.36092559e+00 4.08277678e+00 4.08277678e+00]
 [2.72185119e+00 2.72185119e+00 3.00000000e-16]
 [1.36092559e+00 1.36092559e+00 1.36092559e+00]
 [2.72185119e+00 0.00000000e+00 2.72185119e+00]]
origin coords: 
[[4.08277678e+00 4.08277678e+00 1.36092559e+00]
 [4.50000000e-16 2.72185119e+00 2.72185119e+00]
 [4.08277678e+00 1.36092559e+00 4.08277678e+00]
 [0.00000000e+00 0.00000000e+00 0.00000000e+00]
 [1.36092559e+00 4.08277678e+00 4.08277678e+00]
 [2.72185119e+00 2.72185119e+00 3.00000000e-16]
 [1.36092559e+00 1.36092559e+00 1.36092559e+00]
 [2.72185119e+00 0.00000000e+00 2.72185119e+00]]
new coords: 
[[4.08277678e+00 4.08277678e+00 1.36092559e+00]
 [4.50000000e-16 2.72185119e+00 2.72185119e+00]
 [4.08277678e+00 1.36092559e+00 4.08277678e+00]
 [0.00000000e+00 0.00000000e+00 0.00000000e+00]
 [1.36092559e+00 4.08277678e+00 4.08277678e+00]
 [2.72185119e+00 2.72185119e+00 3.00000000e-16]
 [1.36092559e+00 1.36092559e+00 1.36092559e+00]
 [2.72185119e+00 0.00000000e+00 2.72185119e+00]]
origin rij: 
[ 4.08277678  1.36092559 -1.36092559]
new rij: 
[-1.36092559  1.36092559 -1.36092559]

代码

文本

二、近邻列表的构建和维护

2.1 为什么要使用近邻列表

在上一章中，我们使用LJ势计算了体系的能量和每个原子的受力，并根据原子受力更新了原子的位置和速度。对于一个有 $N$ 个原子的体系，求能量和力时需要遍历原子对，其计算量与 $N^{2}$ 成正比。当更新原子的位置和速度时，需要遍历每个原子，其计算量与 $N$ 成正比。因此可以写出总计算量 $y = c_{0} + c_{1} N + c_{2} N^{2}$ ,其中 $c_{0}$ 、 $c_{1}$ 和 $c_{2}$ 为常数。当 $N$ 很大时， $N^{2}$ 占主导地位，该程序的计算复杂度为 $O (N^{2})$ ，具有平方标度。当粒子数 $N$ 很大时，计算量也会变得很大。

解决上述问题的一个方式是构建近邻列表（neighbor list），在其中提前储存好每一个原子周围的近邻原子，在需要计算时直接使用。近邻列表中有几个关键参数：

近邻列表的截断半径 $R_{c}$ ：在构建近邻列表时，我们只会考虑距离在 $R_{c}$ 以内的原子为近邻原子， $R_{c}$ 应该大于势函数的截断半径 $r_{c}$ 。 $R_{c}$ 越大，近邻原子数量越多，但内存消耗会更大。
缓冲半径 $δ = R_{c} - r_{c}$ ：一般而言，刚构建好的近邻列表在若干步内都是安全的，不需要每一步都更新。若 $δ$ 较小，则每次构建近邻列表所需时间较小，但更新近邻列表的频率较高。一般而言 $δ = 1 \overset{˚}{A}$ 是不错的选择。

alt

2.2 何时更新近邻列表

原子的近邻不是一成不变的，因此构建好的近邻列表需要每隔一段时间更新一下。更新的时机，既可以按照某种固定频率，也可以根据原子坐标变化来判断。下面介绍一种决定何时更新近邻列表的算法：

定义一套额外的参考坐标 $r_{re f}$ ，每个原子的坐标都初始化为0.
在更新原子位置时，对每个原子 $i$ 计算其位移矢量： $d_{i} = ∣ r_{i} - r_{re f} ∣$
若 $2 d_{max} > δ$ ，则更新近邻列表，同时更新 $r_{re f} = r_{i}$

这样，当原子移动了一个较大距离后，程序就会自动更新近邻列表。下面是该算法的一个python实现：

代码

文本

[2]

import numpy as np

def checkIfNeedUpdate(coord_now: np.ndarray, coord_ref: np.ndarray, delta: float=1) -> bool:

'''

检查是否需要更新近邻列表

Parameters

----------

coord_now: np.ndarray, 当前的原子坐标

coord_ref: np.ndarray, 参考坐标

delta: float, 阈值(单位: Angstrom)

Returns

-------

bool, 是否需要更新近邻列表

'''

diff = coord_now - coord_ref

max_displacement = np.max(np.linalg.norm(diff, axis=1))

update = 2 * max_displacement > delta

return update

def main():

# 随机生成原子坐标，用于测试

natom = 100

coord_now = np.random.rand(natom, 3)

coord_ref = np.zeros((natom, 3))

if checkIfNeedUpdate(coord_now, coord_ref):

print("Need update!")

# 更新参考坐标

coord_ref = np.copy(coord_now)

main()

Need update!

代码

文本

2.3 平方标度算法构建近邻列表

近邻列表的平方标度构建方法较为简单：

对于一个原子 $i$ ，遍历该原子和所有原子间的距离，若其间距小于 $R_{c}$ ，则将对应原子记为 $i$ 的近邻。

该算法仍具有 $O (N^{2})$ 的复杂度，但和原来相比，该算法不用在每一步都遍历所有原子对，只需要间隔一段时间更新近邻列表即可，从而节省时间。

下面给出了该算法的一个python实现：

代码

文本

[3]

import numpy as np

from ase.io import read

def find_neighbor_on2(box: np.ndarray, coords: np.ndarray, cutoff: float=3.0):

'''

构造近邻列表，复杂度O(n^2)

Args:

box: 晶胞矩阵

coords: 原子坐标

cutoff: 近邻截断半径

Returns:

neighbor_list: 近邻列表，每一行储存了对应原子的近邻原子的索引

neighbor_number: 每个原子的近邻数目

'''

# 初始化

natoms = len(coords)

max_neighbors = 10

neighbor_list = np.zeros((natoms,max_neighbors), dtype=int)

neighbor_number = np.zeros(natoms, dtype=int)

# 遍历上三角矩阵

for i in range(natoms-1):

for j in range(i+1, natoms):

rij = coords[j] - coords[i]

rij_frac = np.dot(rij, np.linalg.inv(box))

rij_frac[rij_frac > 0.5] -= 1

rij_frac[rij_frac < -0.5] += 1

rij = np.dot(rij_frac, box)

dij = np.linalg.norm(rij)

if dij < cutoff:

neighbor_list[i, neighbor_number[i]] = j

neighbor_list[j, neighbor_number[j]] = i

neighbor_number[i] += 1

neighbor_number[j] += 1

return neighbor_list, neighbor_number

atoms = read("/bohr/md-simulation-linearMD-gk7w/v1/dataset/Si.vasp")

box = atoms.get_cell()

coords = atoms.get_positions()

neighbor_list, neighbor_number = find_neighbor_on2(box, coords)

代码

文本

2.4 线性标度算法构建近邻列表

在MD模拟中最早提出线性标度算法的可能是Quentrec和Brot。该算法的主要思想有以下两点：

整个模拟盒子被划分为一系列小胞，每个小胞的任何厚度均不小于近邻列表的截断半径 $R_{c}$ 。
对于任何一个原子，只需要在27个小胞（一个是该原子所在的小胞，另外26个是与该胞挨着的）中寻找近邻原子。(下图展示了二维空间中的情形，对于任意原子，只需要在该原子本身所在的1个小胞+8个近邻小胞，共9个小胞内寻找近邻原子即可)

这样，我们可以先计算原子在哪个小胞内。然后遍历每个原子，并只在对应的小胞内寻找近邻，避免了遍历所有的原子对，因此提高了效率。

alt

在实现该算法时，首先我们需要确定将整个模拟盒子划分成多少个小胞，可以根据1.1节中的公式来计算盒子在每个方向上的厚度 $L_{a}$ 、 $L_{b}$ 和 $L_{c}$ ，将厚度除以近邻列表的截断半径，并向下取整，即可得到每个方向上的小胞个数： $⎩ ⎨ ⎧ N_{a} N_{b} N_{c} = = = [L_{a} / R_{c}] [L_{b} / R_{c}] [L_{c} / R_{c}]$

下面给出了该算法的一个python实现：

代码

文本

[4]

import numpy as np

from ase.io import read

def getThickness(box: np.ndarray) -> np.ndarray:

"""

Calculate the thickness of the box along the three axes.

The thickness is defined as the volume of the box divided by the area of the face perpendicular to the axis.

Parameters:

box (np.ndarray): 3x3 matrix representing the box vectors.

Returns:

np.ndarray: The thickness along the three axes.

"""

volume = np.abs(np.linalg.det(box))

a = box[0]

b = box[1]

c = box[2]

thickness_a = volume / np.linalg.norm(np.cross(b, c))

thickness_b = volume / np.linalg.norm(np.cross(a, c))

thickness_c = volume / np.linalg.norm(np.cross(a, b))

return np.array([thickness_a, thickness_b, thickness_c])

def getCell(thickness: np.ndarray, box: np.ndarray, coords: np.ndarray, numCells: np.ndarray, cutoff: float) -> np.ndarray:

'''

Calculate the atom in which cell.

Parameters:

thickness (np.ndarray): The thickness of the box along the three axes.

box (np.ndarray): 3x3 matrix representing the box vectors.

coords (np.ndarray): The coordinates of the atoms.

numCells (np.ndarray): The number of cells along the three axes.

cutoff (float): The cutoff radius.

Returns:

cell_index (np.ndarray): The index of the cell in which the atom is located.

'''

coord_frac = np.dot(coords, np.linalg.inv(box))

# 计算晶胞索引，并保证晶胞索引在晶胞范围内

cell_index = np.floor(coord_frac * thickness / cutoff).astype(int)

cell_index = np.mod(cell_index, numCells)

return cell_index

def find_neighbor_on1(coords: np.ndarray, box: np.ndarray, cutoff: float) -> np.ndarray:

'''

Find the neighbors of each atom in the system.

Parameters:

coords (np.ndarray): The coordinates of the atoms.

box (np.ndarray): 3x3 matrix representing the box vectors.

Returns:

neighbor_list (np.ndarray): The list of neighbors for each atom.

neighbor_number (np.ndarray): The number of neighbors for each atom.

'''

# 初始化

natoms = len(coords)

max_neighbor = 1000

neighbor_list = np.full((natoms, max_neighbor), -1, dtype=int)

neighbor_number = np.zeros(natoms, dtype=int)

# 计算晶胞厚度

thickness = getThickness(box)

# 计算整个模拟盒子可以划分的小胞

numCells = np.floor(thickness / cutoff).astype(int)

# 计算原子所在的小胞

cell_index = getCell(thickness, box, coords, numCells, cutoff)

# 建立字典，键为小胞索引，值为该小胞中的原子索引

cell_atoms = {tuple(idx): [] for idx in np.ndindex(*numCells)}

for atom_idx, cell_idx in enumerate(cell_index):

cell_atoms[tuple(cell_idx)].append(atom_idx)

# 计算邻居小胞的偏移

offsets = np.array([[-1,0,1]]*3)

neighbor_offsets = np.array(np.meshgrid(*offsets)).T.reshape(-1, 3)

# 遍历每一个原子

for n in range(natoms):

# 当前原子所在的小胞

current_cell = tuple(cell_index[n])

# 邻居小胞 = 当前小胞 + 偏移

neighbor_cells = (current_cell + neighbor_offsets) % numCells

neighbor_cells = set([tuple(cell) for cell in neighbor_cells])

# 遍历所有邻居小胞中的原子

all_neighbors = np.array([m for cell in neighbor_cells if cell in cell_atoms for m in cell_atoms[cell] if n != m])

if len(all_neighbors) == 0:

continue

# 计算原子之间的距离

rij = coords[all_neighbors] - coords[n]

rij_frac = np.dot(rij, np.linalg.inv(box))

rij_frac[rij_frac > 0.5] -= 1

rij_frac[rij_frac < -0.5] += 1

rij = np.dot(rij_frac, box)

rij_length = np.linalg.norm(rij, axis=1)

# 选取距离小于cutoff的原子

valid_neighbors = all_neighbors[rij_length < cutoff]

if len(valid_neighbors) > max_neighbor:

raise ValueError(f"Atom {n} has {len(valid_neighbors)} neighbors: {valid_neighbors}")

# 保存结果

neighbor_number[n] = len(valid_neighbors)

neighbor_list[n, :len(valid_neighbors)] = valid_neighbors

return neighbor_list, neighbor_number

atoms = read("/bohr/md-simulation-linearMD-gk7w/v1/dataset/Si.vasp")

box = atoms.get_cell()

coords = atoms.get_positions()

neighbor_list, neighbor_number = find_neighbor_on1(coords, box, 3.0)

print(neighbor_number)

[4 4 4 4 4 4 4 4]

代码

文本

三、用ASE或pymatgen快速构建近邻列表(拓展阅读)

Python语言的一大优势是具有丰富的第三方库，可以很方便的帮助我们快速实现某些功能。在材料模拟领域，较常用的库有ASE(全称：Atomic Simulation Environment, 官网：https://wiki.fysik.dtu.dk/ase/) 和Pymatgen（全称：Python Materials Genomics，官网：https://pymatgen.org/）。这里展示如何基于ASE或Pymatgen快速构建近邻列表。这里只展示基本用法，高阶用法可以前往官网阅读手册。

代码

文本

3.1 用ASE构建近邻列表

代码

文本

[7]

# ASE官方手册：https://wiki.fysik.dtu.dk/ase/ase/neighborlist.html

from ase.io import read

from ase.neighborlist import build_neighbor_list

# Read from file

atoms = read('/bohr/md-simulation-linearMD-gk7w/v1/dataset/Si.vasp')

# Build neighbor list

neighbor_list = build_neighbor_list(atoms, self_interaction=False)

print(neighbor_list.get_neighbors(0)[0])

[5 7 1 3]

代码

文本

3.2 用Pymatgen构建近邻列表

代码

文本

[6]

from pymatgen.core import Structure

import numpy as np

# 读取结构文件

structure = Structure.from_file("/bohr/md-simulation-linearMD-gk7w/v1/dataset/Si.vasp")

# 构造近邻列表

center_atom_index, neighbor_atom_index, offset_vectors, distances = structure.get_neighbor_list(r = 3.0)

print(neighbor_atom_index[center_atom_index == 0])

[5 7 1 3]

代码

文本

四、课后作业

通过本节的学习，你应该已经掌握了三斜盒子、近邻列表的相关基础知识，并能通过Python实现相关功能。课后希望各位同学能运用本节所学知识，使用python实现linearMD程序，并完成相应测试。C++版的代码可以参考樊老师的github仓库：https://github.com/brucefan1983/Molecular-Dynamics-Simulation ，Python版的参考代码详见本notebook的数据集。

代码

文本

樊哲勇《分子动力学模拟》程序代码python实现

点个赞吧