从分子动力学建模开始

发表于 2023-09-18 更新于 2024-03-25 分类于热力学与统计物理

原子建模基础

对于较大原子系统，从量子力学出发的计算是不太可能的。通过经典近似来描述原子系统的运动，能够捕捉到宏观系统的许多特征。在这样的近似中，利用牛顿运动定律和对原子间作用力的恰当描述来寻找原子的运动。这些力由给定位置的原子的势能来描述。原则上，势能函数可以从量子力学计算中确定，因为作用力和势能取决于电子的状态 --电子形成了原子间的键，是原子间作用力的来源。

然而，在经典近似中，我们将这个问题参数化：我们为相互作用构建简化模型，这些模型提供了相互作用的最重要特征。

Lennard-Jones势

此处，我们主要是用简单的模型，因为统计和热效应并没有强烈依赖系统的细节。对于两个惰性气体原子间的相互作用力主要有两个来源：

由于偶极相互作用产生的吸引力：这种相互作用的势正比于$(1/r)^6$，其中$r$为原子间的距离。即范德瓦尔斯相互作用
当两个原子被推挤在一起时，由于电子轨道重叠的可能性，存在一种排斥力，这是一种量子力学效应。用$(1/r)^n$来表示。对于氩来说，$n=12$是比较好的近似。

该组合模型称为Lennard-Jones势：

\[\begin{equation} U(r)=4 \epsilon\left(\left(\frac{\sigma}{r}\right)^{12}-\left(\frac{\sigma}{r}\right)^6\right) \end{equation}\]

当L-J势最小时：

\[\begin{equation} F(r)=-\frac{d}{d r} U(r)=0, \end{equation}\]

则有：

\[\begin{equation} F(r)=24 \epsilon_0\left((\sigma / r)^6-2(\sigma / r)^{12}\right)=0 \Rightarrow \frac{r}{\sigma}=2^{1 / 6} \end{equation}\]

L-J不仅可以模拟氩原子之间的相互作用，还可以作为需要仿真的基础模型。

初始条件

MD模拟从原子的初始构型开始，并确定所有原子的运动轨迹。初始条件由初始时刻$t_0$的所有位置$r_i(t_0)$和速度$v_i(t_0)$组成。特别地，由于大多数势随着原子间距趋于零而迅速增加，因此不要将原子放置的太近。因此，我们经常将原子有规律地放置在空间中，在一个晶格上，初始速度是随机的。

我们首先构造一个单胞，并将其复制到晶格的每一个位置，从而形成规则的单胞模式。通过上述方法可以构造出晶格。(单胞可以包含一个以上的原子)。

此处，我们使用立方晶胞。对于单胞中只有一个原子且长度为$b$的立方单胞，我们可以将原子放置在单胞中的(0,0,0)处，然后生成原子间距为$b$的晶格(具体示意图参见原书Fig. 3.1)。

然而，对于L-J体系，从理论实验等研究可知，平衡晶体结构不是简立方，而是面心立方。我们模拟的系统将由$L\times L\times L$个原胞组成，每个原胞尺寸为$b\times b\times b$，且有4个原子。

边界条件

为了利用较小的模型得到更为宏观系统的信息，我们必须采用周期性边界条件。

对运动方程进行积分

我们通过求解运动方程来确定系统的行为。对于分子动力学模拟，通常使用Celocity-Verlet算法。这与前向欧拉法近似，但非常适用于保守力。速度在$t$时刻和中间时刻$t+\Delta t/2$计算，其中$\Delta t$为时间步长。而力仅在全时间步长：$t$、$t+\Delta t$以及$t+2\Delta t$等计算。

计算中最耗时的部分是力的计算。因此，我们希望限制我们计算力的次数，并且在我们的计算中仍然具有尽可能高的精度。Velocity-Verlet算法保证了积分算法中精度和力计算次数之间的良好折衷。

在Velocity-Verlet算法中，所有原子(i = 1，..，N)的位置$r_i( t )$和速度$v_i(t)$在时间上向前传播：

\[\begin{equation} \begin{aligned} \mathbf{v}_i(t+\Delta t / 2) & =\mathbf{v}(t)+\mathbf{F}_i(t) / m_i \Delta t / 2 \\ \mathbf{r}_i(t+\Delta t) & \left.=\mathbf{r}_{(} t\right)+\mathbf{v}_i(t+\Delta t / 2) \\ \mathbf{F}_i(t+\Delta t) & =-\nabla V\left(\mathbf{r}_i(t+\Delta t)\right) \\ \mathbf{v}_i(t+\Delta t) & =\mathbf{v}(t+\Delta t / 2)+\mathbf{F}_i(t+\Delta t) / m_i \Delta t / 2 \end{aligned} \end{equation}\]

这种方法在能量守恒方面具有很好的性质，当然，它也能很好地preserve momentum。

分子动力学模拟器的简单实现

我们需要设置原子的初始构型，对给定的若干时间步的运动进行积分，然后将结果输出为标准可视化程序可以读取的文件。

无量纲运动方程

对于伦纳德-琼斯模型，我们通常使用模型的固有长度($\sigma$)和能量尺度($\epsilon_0$)作为长度和能量的基本单位。因此模拟矢量$\mathbf{r}'_i$与真实世界的$\mathbf{r}_i$有如下关系：

\[\begin{equation} \mathbf{r}_i=\sigma \mathbf{r}_i^{\prime} \Leftrightarrow \mathbf{r}_i^{\prime}=\mathbf{r}_i / \sigma \end{equation}\]

类似的，可以引入L-J时间，$\tau = \sigma\sqrt{m/\epsilon}$。其中，$m$ 是原子质量，以及L-J温度 $T_0 = \epsilon / k_B$。

使用这些符号，我们可以使用无量纲位置和时间重写L-J系统的运动方程。 $\mathbf{r}_i^{\prime}=\mathbf{r}_i/\sigma$和$t^{\prime}=t/\tau$:

\[\begin{equation} m \frac{d^2}{d t^2} \mathbf{r}_i=\sum_j 24 \epsilon_0\left(\left(\sigma / r_{i j}\right)^6-2\left(\sigma / r_{i j}\right)^{12}\right)\left(\mathbf{r}_{i j} / r_{i j}^2\right) \end{equation}\]

变成：

\[\begin{equation} \frac{d^2}{d\left(t^{\prime}\right)^2} \mathbf{r}_i^{\prime}=\sum_j 24\left(r_{i j}^{-6}-2 r_{i j}^{-12}\right)\left(\mathbf{r}_{i j}^{\prime} / r_{i j}^{\prime 2}\right) \end{equation}\]

对于氩气，$\sigma$ $ = 0.3405$ $\mu$ m，以及 $\epsilon$ $= 1.0318^{-2} $ eV。对于其它原子，你需要找到相应的参数值。

量	方程	转换系数	氩的值
长度	$x'=x/L_0$	$L_0=\sigma$	0.3405 $\mu m$
时间	$t'=t/\tau$	$\tau=\sigma \sqrt{m/\epsilon}$	$2.1569\times 10^3$ fs
力	$F'=F/F_0$	$F_0=m \sigma / \tau^2=\epsilon / \sigma$	$3.0303 \times 10^{-1}$
能量	$E'=E/E_0$	$E_0= $	$ 1.0318 ^{-2} $ eV
温度	$T'=T/T_0$	$T_0 = / k_B $	119.74 K

系统初始化

初始化系统：生成一组具有随机速度的原子晶格。该晶格由尺寸为$b$ (单位为$\sigma$)的单胞组成。若单胞的起始为$r_0$，则单胞中的四个原子的坐标为$r_0$、 $r_0 + (b/2, b/2, 0)$、$r_0 + (b/2, 0, b/2)$、以及 $r_0 + (0, b/2, b/2)$。

import numpy as np

L = 5  # 晶格尺寸
b = 2.0  # 晶胞尺寸 (units of sigma)
v0 = 1.0  # 初始化动能尺度
N = 4 * L**3  # Nr of atoms

r = np.zeros((N, 3))
v = np.zeros((N, 3))

bvec = np.array([[0, 0, 0], [b/2, b/2, 0], [b/2, 0, b/2], [0, b/2, b/2]])
ip = 0

# Generate positions
for ix in range(0, L):
    for iy in range(0, L):
        for iz in range(0, L):
            r0 = b * np.array([ix, iy, iz])  # Unit cell base position
            for k in range(0, 4):
                r[ip, :] = r0 + bvec[k, :]
                ip += 1  # Add particle

# Generate velocities
for i in range(0, ip):
    v[i, :] = v0 * np.random.randn(3)

writelammps('mymdinit.lammpstrj', L*b, L*b, L*b, r, v)

将数据写入文件的函数定义如下：

import numpy as np

def writelammps(filename, Lx, Ly, Lz, r, v):
    with open(filename, 'w') as fp:
        ip = r.shape[0]
        fp.write('ITEM: TIMESTEP\n')
        fp.write('0\n')
        fp.write('ITEM: NUMBER OF ATOMS\n')
        fp.write(f'{ip}\n')
        fp.write('ITEM: BOX BOUNDS pp pp pp\n')
        fp.write(f'0.0 {Lx}\n')
        fp.write(f'0.0 {Ly}\n')
        fp.write(f'0.0 {Lz}\n')
        fp.write('ITEM: ATOMS id type x y z vx vy vz\n')
        for i in range(ip):
            fp.write(f'{i+1} 1 {r[i, 0]} {r[i, 1]} {r[i, 2]} {v[i, 0]} {v[i, 1]} {v[i, 2]}\n')

从而可以得到：

ITEM: TIMESTEP
0
ITEM: NUMBER OF ATOMS
500
ITEM: BOX BOUNDS pp pp pp
0.0 10.0
0.0 10.0
0.0 10.0
ITEM: ATOMS id type x y z vx vy vz
1 1 0.0 0.0 0.0 -0.37806747920980305 0.6198930431434015 -0.37486891843336545
2 1 1.0 1.0 0.0 -0.9735100703336169 0.11617643190081425 -1.7552597516341568
3 1 1.0 0.0 1.0 0.27186380165078294 -0.2889571390770045 0.7686946959643555
4 1 0.0 1.0 1.0 -0.8620072632510575 0.18294851178046748 -0.21705159596758267
...

运动方程的积分

为了获取初始状态，需要对运动方程进行积分。首先是读取lammps文件，

from lammpResultsRead import getSimulationSet # 自定义函数

# 获取模拟尺寸

_, bounds_list = getSimulationSet('mymdinit.lammpstrj')

[Lx,Ly,Lz] = bounds_list[:, 1] - bounds_list[:, 0]

t = 3.0; dt = 0.001; n = math.ceil(t/dt);

a = np.zeros((N, 3)) # 存储计算出的加速度

import numpy as np

@njit
def cal(a,N,r,L,v,dt):
    a = np.zeros_like(a)
    for i1 in range(N):
        for i2 in range(i1+1, N):
            dr = r[i1,:] - r[i2,:]
            for k in range(3):     # 周期性边界条件
                while dr[k] > L[k] / 2:
                    dr[k] = dr[k] - L[k]
                    
                while dr[k] < -L[k] / 2:
                    dr[k] = dr[k] + L[k]
            rr = np.dot(dr, dr)
            #print(rr)
            aa = -24 * (2 * (1/rr)**6 - (1/rr)**3) * dr / rr
            a[i1, :] = a[i1, :] + aa[0]  # from i2 on i1
            a[i2, :] = a[i2, :] - aa[1]  # from i1 on i2
    
    v = v + a * dt
    r = r + v * dt
    
    # 周期性边界条件
    
    for i1 in range(N):
        for k in range(3):
            while r[i1, k] > L[k]:
                r[i1, k] = r[i1, k] - L[k]
            while r[i1, k] < 0:
                r[i1, k] = r[i1, k] + L[k]
    return r, v



for i in range(n-1): 
    
    r,v = cal(a,N,r,L,v,dt)
    
    if i%10==0:
        print(i)
    writelammps('mymdinit.lammpstrj',i, L[0], L[1], L[2], r, v)

在上述代码中，我们重新定义了writelammps函数为：

def writelammps(filename,timestep, Lx, Ly, Lz, r, v):
    with open(filename, 'a') as fp:
        ip = r.shape[0]
        fp.write('ITEM: TIMESTEP\n')
        fp.write(f'{timestep}\n')
        fp.write('ITEM: NUMBER OF ATOMS\n')
        fp.write(f'{ip}\n')
        fp.write('ITEM: BOX BOUNDS pp pp pp\n')
        fp.write(f'0.0 {Lx}\n')
        fp.write(f'0.0 {Ly}\n')
        fp.write(f'0.0 {Lz}\n')
        fp.write('ITEM: ATOMS id type x y z vx vy vz\n')
        for i in range(ip):
            fp.write(f'{i+1} 1 {r[i, 0]} {r[i, 1]} {r[i, 2]} {v[i, 0]} {v[i, 1]} {v[i, 2]}\n')

需要注意的是上述代码速度在140步以后会发散，原作者写的matlab代码只是为了展示分子动力学程序的结构和原理，以及如何处理数据和将数据存储在文件中。如果没有加速，源代码是异常慢的。因此是没有必要完全运行该代码的。

进行了分子动力学模拟

安装Lammps

由于学习的时候，我的电脑本身就安装了Lammps，因此该部分略去。

在LAMMPS中进行仿真

# 2d L-J gas
units lj # 该命令选择L-J 单位，长度单位为\sigma，能量单位为\epsilon_0，
# 时间单位为\tao = \sigma \sqrt(m/\epsilon)，温度项T_0 = \epsilon / k_B
# 对于氩气 \sigma = 0.3405 \mu m，\epsilon = 1.0318 \cdot 10^{-2} eV
dimension 2 # 指定维度，此处是2维模拟
boundary p p p # 指定边界条件，此处为周期性边界条件
atom_style atomic # 该命令描述了每个原子/粒子的复杂性，此处我们使用了最简单的描述，atomic，用于惰性气体和粗粒度模拟。

lattice hex 0.75 # 该命令生成晶格原子的坐标，并不生成原子，仅当我们生成原子的时候才有原子。
# hex 二维六方晶格，因此每个原子有六个最近邻原子。
# 当我们选择LJ单元进行模拟时，0.75被称为尺度，是约化密度\rho'。
# (注意到如果我们不使用LJ单位,量表的解释是不同的,参见LAMMPS文档以获得更多的信息)。
region simbox block 0 20 0 10 -0.1 0.1 # 该命令定义了范围，0 < x < 20, 0 < y < 10, -0.1 < z < 0.1； 并将该区域命名为simbox。
create_box 1 simbox # 创建区域，模拟盒子只包含1 ( 1 )种原子，因此数量为1。
create_atoms 1 box  # 向模拟box中填充我们定义的类型为1的原子。

mass 1 1.0 # 该命令定义类型1的原子相对于L-J模型的质量为1。这意味着所有原子的质量都与L-J模型无量纲化时使用的质量m相同。
velocity all create 2.5 87287 # 该命令产生随机速度(高斯分布)，使系统中所有原子类型的初始温度为2.5 (L-J无量纲单位)。
# 87287 是随机数生成器的seed，只要不改变seed值，就会产生与本次模拟相同的初始速度分布。
pair_style lj/cut 2.5 # 该命令制定了我们想要使用的长度为2.5的L-J势，也就是说当距离大于2.5时，相互作用近似为0。
pair_coeff 1 1 1.0 1.0 2.5 # 指定L-J参数，前两个数字表示我们使用的是第1类原子和第1类原子之间的相互作用。
# L-J参数为1.0、1.0以及2.5。这意味着两原子之间的相互作用具有的\sigma值对应于1.0乘以一般的\sigma值；
# \epsilon_0的值对应于1.0乘以总体的\epsilon的值；最后的2.5是cut-off，与上面的2.5一样。

neighbor 0.3 bin  
neigh_modify every 20 delay 0 check no
# 这个模块包含了仿真器是如何使用定期更新的邻居列表来计算原子之间的相互作用的信息。
# 不需要改变这些参数，但改变它们通常不会对模拟结果产生任何影响，只会对模拟的效率产生影响。

fix 1 all nve # 指定在原子上进行何种模拟，这可以通过一个或多个可以应用于原子区域的fix命令来完成。
# 此处，我们称fix为1，被应用于所有粒子(all)。指定模拟在常数nve下运行，固定数量的粒子(n)，定容体积(v，表示模拟box不在过程中发生变化)，以及固定的能量(e)。
# 你可能会惊讶于能量不变的部分。积分算法是否保证能量恒定？是的。
# 然而，也可能存在我们想要给系统的某一部分添加能量的情况，在这种情况下，基本的迭代算法仍然保持能量，但是我们添加了可能改变系统总能量的额外项。

dump 1 all custom 10 dump.lammpstrj id type x y z vx vy vz
# 该命令告诉模拟器输出状态。1是我们给该转存(dump)的名字；
# 输出all原子使用custom输出形式，且每10个时间步输出到文件dump.lammpstrj文件中，后面规定了每个原子输出的是什么。
thermo 100 # 该命令规定了每100个时间步向log.lammps日志文件和终端输出日志。
run 5000 # 该命令开始模拟，并指定它将运行5000个时间步。

Advanced models

Coarsening

对于上述的模型，如果仍以原子均匀分布的状态启动系统个，但初始能量较低，会发生什么呢？此外，我们保持系统的平均动能大致恒定(这意味着我们大约保持系统内温度恒定)。

units lj
dimension 2
boundary p p p
atom_style atomic

lattice hex 0.5
region simbox block 0 80 0 40 -0.1 0.1
create_box 1 simbox
create_atoms 1 box
 
mass 1 1.0
veloctity all create 0.05 87287  # 此处

pair_style lj/cut 2.5
pair_coeff 1 1 1.0 1.0 2.5

neighbor 0.3 bin
neigh_modifty every 20 delay 0 check no

fix 1 all nvt temp 0.25 0.25 1.0 # 此处
dump 1 all atom 1000 dump.lammpstrj
thermo 100
run 50000

这会导致相分离：体系分离成了液体和气体相，展示一种名为“准晶分解”的现象。该现象后续会讨论。

硅晶体的断裂

原书待续

练习

暂时不看

量	方程	转换系数	氩的值
长度	\(x'=x/L_0\)	\(L_0=\sigma\)	0.3405 \(\mu m\)
时间	\(t'=t/\tau\)	\(\tau=\sigma \sqrt{m/\epsilon}\)	\(2.1569\times 10^3\) fs
力	\(F'=F/F_0\)	\(F_0=m \sigma / \tau^2=\epsilon / \sigma\)	\(3.0303 \times 10^{-1}\)
能量	\(E'=E/E_0\)	$E_0= $	$ 1.0318 ^{-2} $ eV
温度	\(T'=T/T_0\)	$T_0 = / k_B $	119.74 K