Chapt 6: Molecule Dynamics

Author
Affiliation

Feng Wan

B832@仲英楼 or 3133@泓理楼 // School of Physcis @ Xian Jiaoton University

Show the code 
import numpy as np
import matplotlib.pyplot as plt
from sympy import *
from IPython.display import Markdown
from sympy.solvers.ode.systems import dsolve_system
from sympy import *

1 Final project

  • 自选大作业格式要求：
  • 选题背景介绍
  • 物理模型与理论
  • 程序求解
  • 结果分析 +（发现抄袭，一律挂科）

2 Molecule Dynamics - 引言

  • Monte Carlo 方法 是一种随机模拟方法，或统计试验方法。通过不断产生随机数序列来模拟过程。
  • 分子动力学（Molecular Dynamics，简称MD） 是一种确定性模拟方法。通过数值求解一个个的粒子运动方程来模拟整个系统的行为。
    • 按系统内部的内禀动力学规律来计算并确定位形的转变。
    • 分子动力学方法中不存在任何随机因素。可以看作是体系在一段时间内的发展过程的模拟。
    • 计算量大，占内存多

适用的体系： 少 体 —— 多体

              点粒子  ——  具有内部结构的体系

处理的微观客体：分 子 —— 其它微观粒子

2.1 分子动力学方法发展历史

  • Alder和Wainght在1957年至1959年间应用于理想“硬球”液体模型，很多简单液体中分子之间的相互作用的重要性质在两人的研究中被发现 ；
  • Rahman于1964年应用一种更接近的液体模型模拟了液氦；
  • Verlet于1967~1971年，Levesque于1983年做了进一步的研究，模拟了更复杂的液体——水、溶盐、聚合物和蛋白质等；
  • 使用分子动力学模拟研究真实系统在1974年由Rahman和Stilliger完成，研究对象为液体水；
  • 第一个蛋白质的分子动力学研究报道是McCammon在1977年，研究对象是牛胰岛素抑制剂。

2.2 分子动力学实例

蛋白质的折叠

通过分子动力学模拟 Barnase （芽孢杆菌RNA酶 ）的解折叠过程，反过来研究其可能的折叠路径。

  • 水化（溶剂化）

  • 描述生物大分子周围的溶剂（水）分子分布。

  • 抹香鲸肌红蛋白的水化

  • 水化位置（蓝，实验）

  • 水化数密度最大值（黄，MD 模拟）

材料的脆性断裂

  • 分子动力学是在原子、分子水平上求解多体问题的重要的计算机模拟方法。

  • 通过求解所有粒子的运动方程，分子动力学方法可以用于模拟与粒子运动路径相关的基本过程。

  • 在经典分子动力学中，粒子的运动行为是通过经典的Newton运动方程所描述。系统的所有粒子服从经典力学的运动规律，它的动力学方程就是从经典力学的运动方程——拉格朗日(lagrange)方程和哈密顿(Hamilton)方程导出。

  • 分子动力学得天独厚的优势是能够模拟分子的运动轨道，即能够提供分子运动和变化的最为微观的、细致的信息。这对传统的力学、热力学、光谱学等实验方法都是难以办到的。
  • 即使是现代分析理论和实验测量手段，要想准确的描述物质的非平衡态性质都是非常困难的。

  • 如果能够获得分子动力学计算得到的运行轨道，就能通过非常直接的计算公式给出感兴趣的非平衡态物理量。

  • 这种特点对现今所有的理论、实验、计算手段来讲都是独一无二的。

  • 正是这种不可替代的作用，使得分子动力学在当今自然科学研究中获得了如此广泛和深入的应用。

  • 实际上，分子动力学模拟方法和随机模拟方法一样都面临着两个基本限制：
  • 一个是有限观测时间的限制；
    • 另一个是有限系统大小的限制。
    • 通常人们感兴趣的是体系在热力学极限下（即粒子数目趋于无穷时）的性质。
    • 但是计算机模拟允许的体系大小要比热力学极限小得多，因此可能会出现有限尺寸效应。
    • 为了减小有限尺寸效应，人们往往引入周期性、全反射、漫反射等边界条件。
    • 当然边界条件的引入显然会影响体系的某些性质。
  • 对体系的分子运动方程组采用计算机进行数值求解时，需要将运动方程离散化为有限差分方程。
    • 常用的求解方法有欧拉法、龙格-库塔法、Verlet算法等。
    • 数值计算的误差阶数显然取决于所采用的数值求解方法的近似阶数。
    • 原则上，只要计算机计算速度足够大，内存足够多，我们可以使计算误差足够小。

2.3 分子动力学运动方程数值求解

  1. 系统的动力学机制决定运动方程的形式。
  2. 在分子动力学方法处理过程中，方程组的建立是通过对物理体系的微观数学描述给出的。
  3. 在这个微观的物理体系中，每个分子都各自服从经典的牛顿力学。
  4. 每个分子运动的内禀动力学是用理论力学上的哈密顿量或者拉格朗日量来描述，也可以直接用牛顿运动方程来描述。
  5. 这种方法可以处理与时间有关的过程，因而可以处理非平衡态问题。

2.3.1 运动方程

  1. 采用分子动力学方法时，必须对一组分子运动微分方程做数值求解。
  2. 从计算数学的角度来看，这是个求一个初值问题的微分方程的解。
  3. 实际上计算数学为了求解这种问题已经发展了许多的算法。
  4. 但是并不是所有的这些算法都可以用来解决物理问题。

2.3.2 空间描述

  • 在空间描述物体的运动，如果其本身的大小可以忽略时，就可以将其看作是粒子（或质点）。
  • 粒子描述：
    • 空间位置：\(r\)
    • 速度：\(v = dr/dt\)
    • 加速度：\(a = dv/dt = d^2r/dt^2\)
  • 若一个系统由 \(N\) 个粒子组成，则粒子描述：
    • 空间位置：\(r^1, r^2, ..., r^N\)
    • 笛卡尔坐标系，粒子有 \(3N\) 个自由度
  • 设系统有 \(s\) 个自由度
    • 广义坐标： \(q_1, q_2, ..., q_N\)
    • 广义速度： \(\dot{q}_1, \dot{q}_2, ..., \dot{q}_N\)

2.4 最小作用量原理

  • 莫培督1744年提出最小作用量原理：
    • 保守的、完整的力学系统，由某一初位形转变到另一位形的一切具有相同能量的可能运动中，真实的运动是其作用量具有极小值的那种运动。
  • 力学系统中，构造能量函数 \(L\) 及其作用量 \(S\) \[ \begin{aligned}[l] L & = & L\left[ q(t), \dot{q}(t), t\right] \\ S\left[q(t)\right] &=& \int_{t_1}^{t_2} L\left[q(t), \dot{q}(t), t\right]dt \end{aligned} \]
  • 作用量的积分式叫做泛函(functional)，作用量取极值的方法就是求其变分 \(\delta S = 0\) \[ \delta S = \delta \int_{t_1}^{t_2} L\left[q(t), \dot{q}(t), t\right]dt = 0 \]

2.5 拉格朗日(Lagrange)方程

  • 由最小作用量原理可导出拉格朗日方程 \[ \frac{d}{dt} \frac{\partial L}{\partial \dot{q}_i} - \frac{\partial L}{\partial q_i} = 0 \]
  • 对于孤立的保守系统，每个粒子在势场U中运动，则 \[ L(q_i, \dot{q}_i, t) = \frac{1}{2}m_i \dot{q}_i^2 - U_i \]
  • 系统整体的Lagrange函数是 \[ L(q_i, \dot{q}_i, t) = \sum_i L(q_i, \dot{q}_i, t) = \frac{1}{2}\sum_i m_i \dot{q}_i^2 - \sum_i U_i \]
  • 带入拉格朗日方程可以得到第i个粒子的牛顿运动方程（ \(\alpha\) 指每个粒子的自由度） \[ m_i \ddot{q}_{i\alpha} = - \frac{\partial U_i}{\partial q_{i\alpha}} = -\Delta U_i \]

2.6 哈密顿(Hamilton)方程

哈密顿函数 \(H\) 是动量和坐标的函数，是动能和势能之和 \[ H = \sum_i \frac{p_i^2}{2m} + \sum_{i < j}u(r_j) \] 变量为动量 \(p\) 和坐标 \(r\) 的Hamilton方程 \[ \dot{p}_i = - \frac{\partial H}{\partial q_i} ~~~~ \dot{q}_i = \frac{\partial H}{\partial p_i} \] 如果系统的哈密顿函数不显含时间，就有 \(dH/dt=0\)，即得到能量守恒定律。

2.7 粒子运动方程的数值解法

设粒子的坐标、动量及其作用力分别用 \(x(t)\)，\(v(t)\)，\(p(t)\)，\(f(x,t)\)表示，其速度初始值为 \(x(0)\)，\(v(0)\)，\(p(0)\)，\(f(0)\)。则决定粒子运动的牛顿方程是 \[ m\ddot{r} = - \frac{\partial U }{\partial r} \] 同时 \[ \left\lbrace \begin{aligned}[l] v(t) &=& \dot{r} = \frac{\partial H}{\partial p} = \frac{p}{m} \\ m\ddot{r} &=& \frac{dp}{dt} = -\frac{\partial H}{\partial r} = - \frac{\partial U}{\partial r} \end{aligned} \right. ~~~~~ \text{IVP of ODE} \]

2.8 差分计算的二步法与Verlet算法

用更高阶的泰勒展开 \[ \left\lbrace \begin{aligned}[l] f(t+h) &=&f(t) + hf_t + \frac{h^2}{2}f_{tt} + O(h^3) \\ f(t-h) &=&f(t) - hf_t + \frac{h^2}{2}f_{tt} + O(h^3) \end{aligned} \right. \]

\[ \Rightarrow \left\lbrace \begin{aligned}[l] f_{tt} &=& \frac{1}{h^2}[f(t+h) - 2f(t) + f(t-h)] + O(h^2) \\ f_{t} &=& \frac{f(t+h) - f(t-h)}{2h} + O(h^2) \end{aligned} \right. \]

Verlet method

\[ \xrightarrow[\ddot{r}(t) = \frac{F(t)}{m}]{r(t)\equiv f(t), ~~ p(t) = m\dot{r}(t)} \left\lbrace \begin{aligned}[l] r(t+h) &=& 2r(t) - r(t - h) + h^2 \frac{F(t)}{m} + O(h^4) \\ p(t) &=& \frac{m}{2h}[r(t+h) -r(t - h)] + O(h^2) \end{aligned} \right. \] \(p\) 的求解要比 \(x\) 落后一步。

2.8.1 Verlet算法的速度形式

\[ \left\lbrace \begin{aligned}[l] f(t+h) &=&f(t) + hf_t + \frac{h^2}{2}f_{tt} + O(h^3) \\ f(t-h) &=&f(t) - hf_t + \frac{h^2}{2}f_{tt} + O(h^3) \end{aligned} \right. \]

\[ \rightarrow \left\lbrace \begin{aligned}[l] r(t+h) &=& r(t) + hv(t) + h^2 \frac{F(t + h)}{2m} ~~~~~~ (3) \\ r(t-h) &=& r(t) - hv(t) + h^2 \frac{F(t + h)}{2m} ~~~~~~ (3') \end{aligned} \right. \] 匀加速直线运动的位移公式

(3’) push forward with step \(h\): \[ r(t) = r(t + h) - hv(t + h) + h^2 \frac{F(t + h)}{2m} ~~~~~~(4) \] plug (3) into (4) \[ v(t + h) = v(t) + h \frac{F(t+h)+F(t)}{2m} \]

Note 1: Verlet算法的速度形式

\[ \left\lbrace \begin{aligned}[l] r_i(t + h) &=& r_i(t) + hv_i(t) + h^2\frac{F_i(t)}{2m_i} \\ v_i(t + h) &=& v_i(t) + h\frac{F_i(t+h) + F_i(t)}{2m_i} \\ F_i(t) &=& -\sum \nabla_j V(r_{ij}) \end{aligned} \right. \] \(x\) 和 \(v\) （其实是动量 \(p\) ）的求解同步，且算法的数值计算的稳定性也得到了很大的加强。

Example 1 (对于一维谐振子) \[ H = \frac{p^2}{2m} + \frac{1}{2}kx^2; ~~~~~~ F(t) = -kx \] \[ \rightarrow \left\lbrace \begin{aligned}[l] x(t + h) &=& x(t) + hv(t) + h^2\frac{F(t)}{2m} \\ v(t + h) &=& v(t) + h\frac{F(t+h) + F(t)}{2m} \\ \end{aligned} \right. \]

Code:

%filename: chapter6_example_1_Euler.m
%一维谐振子运动方程的数值求解；
clear all
close all
clc
tic
m=1; %谐振子质量
k=1;%弹簧劲度系数


h=0.001; %时间步长,0.5,0.2,0.05
x0=0;   %初始位置,可以任意设置
p0=100;   %初始动量,可以任意设置
a=0;    %时间起点
b=100;   %时间终点
t=a:h:b;
N=(b-a)/h;
x=zeros(1,N+1);
p=zeros(1,N+1);
x(1)=x0;
p(1)=p0;

for n=1:N %差分算法的单步法，Euler法迭代公式
    x(n+1)=x(n)+h*p(n)/m;
    p(n+1)=p(n)-h*k*x(n);
end

A=sqrt(x0^2+p0^2/(m*k)); %谐振子的振幅
pha=-acos(x0/A);          %谐振子的初始相位，若x0与p0同号，pha为负，否则pha为正
xx=A*cos(sqrt(k/m)*t+pha); %谐振子的精确解
pp=-m*A*sqrt(k/m)*sin(sqrt(k/m)*t+pha); %动量的精确解
EE=pp.^2/(2*m)+0.5*k*xx.^2; %能量的精确值
E=p.^2/(2*m)+0.5*k*x.^2; %能量的数值解
subplot(3,1,1),plot(t,x,'r',t,xx,'b')
title('谐振子的x(t)曲线')
xlabel('t');ylabel('x');
subplot(3,1,2),plot(x,p,'r',xx,pp,'b')
title('谐振子的x-p相图')
xlabel('x');ylabel('P');
subplot(3,1,3),plot(t,E,'r',t,EE,'b')
title('谐振子的能量随时间变化关系')
xlabel('t');ylabel('E');
max(x-xx)
toc
%filename: chapter6_example_2_RK.m
%利用四阶RK算法求解一维谐振子的运动方程;
function two_variables_RK
clear all
close all
clc

tic 
m=1; %谐振子质量
k=1;%弹簧劲度系数
h=0.01; %时间步长,0.5,0.2,0.05
x0=0;   %初始位置,可以任意设置
p0=100;   %初始动量,可以任意设置
a=0;    %时间起点
b=100;   %时间终点
t=a:h:b;
N=(b-a)/h;
xx=zeros(1,N+1);
pp=zeros(1,N+1);
%--------------RK算法------------------------
f=@(t,x)[x(2)/m, -k*x(1)]; % 定义 f，x(2)=pp,x(1)=xx
z0 = [x0, p0]; 
z = myrungekutta(f, t, z0);
xx = z(:,1); 
pp = z(:,2); 
A=sqrt(x0^2+p0^2/(m*k));  %谐振子的振幅
pha=-acos(x0/A);          %谐振子的初始相位，若x0与p0同号，pha为负，否则pha为正
xxx=A*cos(sqrt(k/m)*t+pha); %谐振子的精确解
ppp=-m*A*sqrt(k/m)*sin(sqrt(k/m)*t+pha); %动量的精确解
EE=ppp.^2/(2*m)+0.5*k*xxx.^2; %能量的精确值
E=pp.^2/(2*m)+0.5*k*xx.^2; %能量的数值解
toc
subplot(2,1,1),plot(t,xx,'r',t,xxx,'b')
title('位置随时间的演化关系')
legend('数值解','精确解')
subplot(2,1,2),plot(xx,pp,'r')
title('位置-动量相图')
figure;subplot(2,1,1),plot(t,E,'r',t,EE,'b')
xlabel('时间')
ylabel('能量')
title('谐振子的能量随时间变化关系')
subplot(2,1,2);plot(t,E-EE')
xlabel('时间')
ylabel('能量误差')
title('谐振子能量的数值解与解析解的误差随时间变化关系')


function z = myrungekutta(f, t, z0)
n = length(t);
h = t(2)-t(1);
z=zeros(n,2);
z(1,:)=z0;
for i = 1:n-1
    k1 = h*f(t(i), z(i,:));
    k2=h*f(t(i)+0.5*h, z(i,:) + k1/2);
    k3=h*f(t(i)+0.5*h, z(i,:) + k2/2);
    k4=h*f(t(i)+h, z(i,:) + k3);
    z(i+1,:)=z(i,:) + (k1 + 2*k2 + 2*k3 + k4)/6;
end


%filename: chapter6_example_2_VVerlet.m
%一维谐振子运动方程的数值求解-velocity-verlet算法;
clear all
close all
clc
tic  %获得CPU时间
m=1; %谐振子质量
k=1; %弹簧劲度系数
h=0.01; %时间步长,0.5,0.2,0.05
x0=0;   %初始位置,可以任意设置
p0=100;   %初始动量,可以任意设置
a=0;    %时间起点
b=100;   %时间终点
t=a:h:b;
N=(b-a)/h;
x=zeros(1,N+1);
p=zeros(1,N+1);
x(1)=x0;
p(1)=p0;

for n=1:N %速度Verlet算法
    x(n+1)=x(n)+h*p(n)/m-h^2*k*x(n)/(2*m);
    p(n+1)=p(n)+h*(-k*x(n+1)-k*x(n))/2;
end

A=sqrt(x0^2+p0^2/(m*k)); %谐振子的振幅
pha=-acos(x0/A);          %谐振子的初始相位，若x0与p0同号，pha为负，否则pha为正
xx=A*cos(sqrt(k/m)*t+pha); %谐振子的精确解
pp=-m*A*sqrt(k/m)*sin(sqrt(k/m)*t+pha); %动量的精确解
EE=pp.^2/(2*m)+0.5*k*xx.^2; %能量的精确值
E=p.^2/(2*m)+0.5*k*x.^2; %能量的数值解
toc

subplot(2,1,1),plot(t,x,'r',t,xx,'b')
title('谐振子的x(t)曲线')
subplot(2,1,2),plot(x,p,'r',xx,pp,'b')
title('谐振子的x-p相图')
figure;subplot(2,1,1),plot(t,E,'r',t,EE,'b')
xlabel('时间')
ylabel('能量')
title('谐振子的能量随时间变化关系')
subplot(2,1,2);plot(t,E-EE)
xlabel('时间')
ylabel('能量误差')
title('谐振子能量的数值解与解析解的误差随时间变化关系')

2.9 Python version

import numpy as np
def euler(ode, x0, y0, h, xp):
    steps = int((xp - x0) / h)
    xi = np.zeros(steps + 1)
    yi = np.zeros((steps + 1, len(y0)))
    xi[0] = x0
    yi[0] = y0
    for i in np.arange(steps):
        xi[i + 1] = xi[i] + h
        yi[i + 1] = yi[i] + h * ode(xi[i], yi[i])
    return xi, yi


def adams2step(ode, x0, y0, h, xp):
    steps = int((xp - x0) / h)
    xi = np.zeros(steps + 1)
    yi = np.zeros((steps + 1, len(y0)))
    xi[0] = x0
    yi[0] = y0
    
    ym = y0
    xm = x0
    
    xi[1] = x0 + h
    yi[1] = y0 + h * ode(x0, y0)
    for i in np.arange(1, steps):
        xi[i + 1] = xi[i] + h
        yp = yi[i] + h * (1.5 * ode(xi[i], yi[i]) - 0.5 * ode(xm, ym))
        xm = xi[i]
        ym = yi[i]
        yi[i + 1] = yp
    return xi, yi 


def rk4(func, x0, y0, h, xp):
    steps = int((xp - x0) / h)
    xi = np.zeros(steps + 1)
    yi = np.zeros((steps + 1, len(y0)))
    xi[0] = x0
    yi[0] = y0
    for i in range(steps):
        k1 = func(xi[i], yi[i])
        k2 = func(xi[i] + 0.5 * h, yi[i] + 0.5 * h * k1)
        k3 = func(xi[i] + 0.5 * h, yi[i] + 0.5 * h * k2)
        k4 = func(xi[i] + h, yi[i] + h * k3)
        xi[i + 1] = xi[i] + h
        yi[i + 1] = yi[i] + h / 6.0 * (k1 + 2.0 * k2 + 2.0 * k3 + k4)
    return xi, yi
Show the code 
from python.odefunctions import *

def ode(x, y):
    m = 1
    k = 1
    res = y.copy()
    res[0] = y[1] / m
    res[1] = -k * y[0]
    return res

def verlet(x0, y0, h, xp):
    m = 1
    k = 1
    steps = int((xp - x0) / h)
    xi = np.zeros(steps + 1)
    yi = np.zeros((steps + 1, len(y0)))
    xi[0] = x0
    yi[0, :] = y0
    for i in np.arange(steps):
        xi[i + 1] = xi[i] + h
        yi[i + 1, 0] = yi[i, 0] + h * yi[i, 1] + 0.5 * h * h / m * (-k * yi[i, 0])
        yi[i + 1, 1] = yi[i, 1] + h * 0.5 / m * (-k * (yi[i, 0] + yi[i + 1, 0]))
    return xi, yi

x0 = 0
xp = 100.0
h = 0.1
y0 = np.array([0, 100.0])

x_e, y_e = euler(ode, x0, y0, h, xp)
x_rk, y_rk = rk4(ode, x0, y0, h, xp)
x_v, y_v = verlet(x0, y0, h, xp)

plt.figure(figsize = (12, 4))
plt.subplot(121)
##plt.plot(x_e, y_e[:, 0], 'r-', label = 'euler')
plt.plot(x_rk, y_rk[:, 0], 'b-.', label = 'rk4')
plt.plot(x_v, y_v[:, 0], 'c-', label = 'verlet')
plt.legend()
plt.xlabel(r'$t$')
plt.ylabel(r'$x$')
plt.subplot(122)
#plt.plot(y_e[:, 0], y_e[:, 1], 'r-', label = 'euler')
plt.plot(y_rk[:, 0], y_rk[:, 1], 'b-.', label = 'rk4')
plt.plot(y_v[:, 0], y_v[:, 1], 'c-', label = 'verlet')
plt.legend()
plt.xlabel(r'$x$')
plt.ylabel(r'$p$')
plt.show()

plt.figure(figsize = (12, 4))
plt.subplot(121)
plt.plot(x_rk, y_rk[:, 0] - np.sin(x_rk) * 100.0, 'b-.', label = 'rk4')
plt.plot(x_v, y_v[:, 0] - np.sin(x_v) * 100.0, 'c-', label = 'verlet')
plt.legend()
plt.xlabel(r'$t$')
plt.ylabel(r'$x$')
plt.subplot(122)
plt.plot(x_rk, y_rk[:, 1] - np.cos(x_rk) * 100.0, 'b-.', label = 'rk4')
plt.plot(x_v, y_v[:, 1] - np.cos(x_v) * 100.0, 'c-', label = 'verlet')
plt.legend()
plt.xlabel(r'$x$')
plt.ylabel(r'$p$')
plt.show()

ek_rk = 0.5 * (y_rk[:, 0]**2.0 + y_rk[:, 1]**2.0)
ek_v = 0.5 * (y_v[:, 0]**2.0 + y_v[:, 1]**2.0)


plt.plot(x_rk, ek_rk, 'r-', label = 'rk4')
#plt.plot(x_v, ek_v, 'c-.', label = 'verlet')
plt.plot(x_rk, np.ones(x_rk.shape) * 5000.0, 'k-')
plt.show()

与单步法相比，显然在相同的步长下，Velocity Verlet算法具有更高的精度。

比较Velocity Verlet算法与四阶RK算法的耗时

Example 2 (利用分子动力学方法研究哈雷彗星的运动)  

哈雷彗星与太阳间的势函数为 \[V(r) = - G\frac{Mm}{r}\]

在质心坐标系下 \[\mu \frac{d^2 \vec{r}}{dt^2} = \vec{f} = -GMm \frac{\vec{r}}{r^3}\]

其中 \(\mu = \frac{Mm}{M + m} \approx m\)

以太阳位置为原点，远日点为初始位置，则 \(x_0 = r_\max, v_{x0} = 0, y_0 = 0\), and \(v_{y0} = v_\min.\) \(r_\max = 5.28 \times 10^{12}~\text{m}\), \(v_\min = 9.13 \times 10^2~\text{m/s}\).

分别利用Velocity Verlet算法和四阶RK算法研究哈雷彗星的运动

Show the code 
def plotcircle(ax, r, mystr, pos):
    theta = np.linspace(0, 2.0 * np.pi, 200)
    xx = np.sin(theta) * r
    yy = np.cos(theta) * r
    ax.plot(xx, yy, 'k-')
    ax.text(pos, yy[0] + 0.02, mystr, fontsize = 14)

def plotelipse(ax, a, b):
    theta = np.linspace(0, 2.0 * np.pi, 200)
    xx = np.sin(theta) * a
    yy = np.cos(theta) * b + b * 0.9
    ax.plot(xx, yy, 'k-.')
    ax.plot(xx[0], yy[0], 'bo', markersize = 15)
    ax.text(-0.12, yy[0] + 0.05, "Halley", fontsize = 14)

fig = plt.figure(figsize = (3, 9))
ax = plt.axes()
ax.plot(0, 0, 'ro', markersize = 15)
plotcircle(ax, 0.2, 'Earth', -0.1)
plotcircle(ax, 0.3, 'Mars', -0.1)
plotcircle(ax, 0.5, 'Mercury', -0.15)
plotcircle(ax, 0.8, 'Saturn', -0.13)
plotcircle(ax, 1.1, 'Uranus', -0.13)
plotcircle(ax, 1.6, 'Neptunus', -0.16)
plotelipse(ax, 0.2, 1.0);
ax.set_xlim([-0.4, 0.4]);
ax.set_ylim([-0.2, 2.2]);
plt.axis('off');

\[ \begin{eqnarray} (x, y): ~~x^{k+1} &=& x^k + v_x^k \tau + \frac{\tau^2}{2} g_x^k, \\ (v_x, v_y):~~~v_x^{k+1} &=& v_x^k + \frac{\tau}{2} (g_x^{k+1} + g_x^k), \\ (g_x, g_y):~~~g_x &=& -\kappa \frac{x}{r^3}, \\ r = \sqrt{x^2 + y^2},& & \kappa = GM, \end{eqnarray} \]

Code:

%filename: chapter6_example_3_Halley_VVerlet.m
%哈雷彗星的运动轨道求解-velocity-verlet算法
clear all
close all
clc
%长度单位我们用2.68*10^12为标度，则远日电rmax=1.97
%时间单位我们用哈雷彗星的周期76年为标度
tic
k=39.478428; %GM乘积

h=0.00001; %时间步长,0.0001时误差还较大,注意：时间单位为76年,步长为6.66小时
a=0;    %时间起点
b=3;   %时间终点
t=a:h:b;
N=(b-a)/h;
x0=1.966843;   %初始位置,远日点
y0=0;
vx0=0;     %初始速度
vy0=0.815795; 

x=zeros(1,N+1);
y=zeros(1,N+1);
vx=zeros(1,N+1);
vy=zeros(1,N+1);
x(1,1)=x0;
y(1,1)=y0;
vx(1,1)=vx0;
vy(1,1)=vy0;
%axis([-0.5 2 -0.3 0.3])
for n=1:N %速度Verlet算法
     %hold on
     %plot(x(n),y(n),x(1),y(1),'r*','markersize',8);
     %pause(0.001);
    x(1,n+1)=x(1,n)+h*vx(1,n)-h^2/2*k*x(1,n)/sqrt(x(1,n)^2+y(1,n)^2)^3;
    y(1,n+1)=y(1,n)+h*vy(1,n)-h^2/2*k*y(1,n)/sqrt(x(1,n)^2+y(1,n)^2)^3;
    vx(1,n+1)=vx(1,n)+h/2*(-k*x(1,n+1)/(x(1,n+1)^2+y(1,n+1)^2)^1.5...
    -k*x(1,n)/(x(1,n)^2+y(1,n)^2)^1.5);
    vy(1,n+1)=vy(1,n)+h/2*(-k*y(1,n+1)/(x(1,n+1)^2+y(1,n+1)^2)^1.5...
    -k*y(1,n)/(x(1,n)^2+y(1,n)^2)^1.5);
end
  %hold off
toc
subplot(2,1,1),plot(t,sqrt(x.^2+y.^2),'r')
xlabel('t')
ylabel('r')
title('哈雷彗星的r(t)曲线')
subplot(2,1,2),plot(x,y,'b',x(1),y(1),'r*','markersize',8)
hold on
plot(0,0,'ro','markersize',8)
hold off
title('哈雷彗星的轨道曲线')

RK4:

%filename: chapter6_example_3_Halley_RK.m
%利用四阶RK算法求解哈雷彗星的运动轨道求解
function four_variables_RK
clear all
close all
clc

tic 
k=39.478428; %GM乘积

h=0.00001; %时间步长,0.0001时误差还较大,注意：时间单位为76年,步长为6.66小时
a=0;    %时间起点
b=3;   %时间终点
t=a:h:b;
N=(b-a)/h;
x0=1.966843;   %初始位置,远日点
y0=0;
vx0=0;     %初始速度
vy0=0.815795; 

x=zeros(1,N+1);
y=zeros(1,N+1);
vx=zeros(1,N+1);
vy=zeros(1,N+1);
%--------------RK算法------------------------
f=@(t,z)[z(3),z(4),-k*z(1)/sqrt(z(1)^2+z(2)^2)^3,-k*z(2)/sqrt(z(1)^2+z(2)^2)^3]; % 定义 f，y(2)=pp,y(1)=xx
z0 = [x0, y0, vx0, vy0]; 

z = myrungekutta(f, t, z0);
x  = z(:,1); 
y  = z(:,2); 
vx = z(:,3); 
vy = z(:,4); 
toc
subplot(2,1,1),plot(t,sqrt(x.^2+y.^2),'r')
title('哈雷彗星的r(t)曲线')
subplot(2,1,2),plot(x,y,'b',x(1),y(1),'r*','markersize',8)
title('哈雷彗星的轨道曲线')
hold on
plot(0,0,'ro','markersize',8)
hold off



function z = myrungekutta(f, t, z0)
n = length(t);
h = t(2)-t(1);
z=zeros(n,4);
z(1,:)=z0;
for i = 1:n-1
    k1 = h*f(t(i), z(i,:));
    k2 = h*f(t(i)+0.5*h, z(i,:) + k1/2);
    k3 = h*f(t(i)+0.5*h, z(i,:) + k2/2);
    k4 = h*f(t(i)+h, z(i,:) + k3);
    z(i+1,:)=z(i,:) + (k1 + 2*k2 + 2*k3 + k4)/6;
end

3 原胞与边界条件

  • 当系统粒子数目很少时，可以对所有粒子进行MD模拟；
  • 通常来说，对于一小块样品，其内部粒子数也是很多的。
  • 要想精确求出系综中粒子间相互作用和运动状态是不可能的，只有通过少量粒子的运动性质计算来求出宏观性质。
  • 由于粒子数过少，在构建的计算原胞边缘的粒子较多，其受力状态与内部粒子不同。
  • 为了解决这一问题，分子动力学方法需要对计算原胞设置适当的边界条件。

3.1 原胞

  • 分子动力学模拟方法往往用于研究大块物质在给定密度下的性质，而实际计算模拟不可能在几乎是无穷大的系统中进行。
  • 所以必须引进一个叫做分子动力学原胞的体积元，以维持一个恒定的密度。
  • 对气体和液体，如果所占体积足够大，并且系统处于热平衡状态的情况下，那么这个体积的形状是无关紧要的。
  • 对于晶态的系统，原胞的形状是有影响的。
  • 为了计算简便，取一个立方形的体积为分子动力学原胞。设原胞的线度大小为 \(L\)，则体积为 \(L^3\)。
  • 由于引进这样的立方体箱子，将产生六个我们不希望出现的表面。
  • 模拟中碰撞这些箱子的表面的粒子应当被反射回到原胞内部，特别是对粒子数目很少的系统。
  • 这些表面的存在对系统的任何一种性质都会有重大的影响。

3.2 边界条件

  • 采用分子动力学方法，必需对被计算的粒子系统给定适当的边界条件。这些边界条件大致可分成四种。

    • 自由表面边界(free-surface boundary)这种边界条件常用于大型的自由分子的模拟。
    • 固定边界(rigid boundary)在所有要计算到粒子晶胞之外还要包上几层结构相同的位置不变的粒子，包层的厚度必须大于粒子间相互作用的力程范围。包层部分代表了与运动粒子起作用的宏观晶体的那一部分。
    • 柔性边界 (flexible boundary)这种边界比固定边界更接近实际。它允许边界上的粒子有微小的移动以反映内层粒子的作用力施加到它们身上时的情况。
    • 周期性边界(periodic boundary)在模拟较大的系统时，为了消除表面效应或边界效应，常采用周期性边界条件。就是让原胞上下、左右、前后对边上的粒子间有相互作用。

3.3 周期性边界条件

  • 为了将分子动力学原胞有限立方体内的模拟，扩展到真实大系统的模拟，通常采用周期性边界条件。
  • 采用周期性边界条件，就可以消除引入原胞后的表面效应，构造出一个准无穷大的体积来更精确地代表宏观系统。
  • 实际上，这里我们做了一个假定，即让这个小体积原胞镶嵌在一个无穷大的大块物质之中。
  • 做一个假设：让这个小体积原胞镶嵌在一个无穷大的大块物质中，周期性边界条件的数学表示形式为： \[ A(\vec{x}) = A(\vec{x} + \vec{n}L), ~~~ \vec{n} = (n_1, n_2, n_3) \]
  • 其中 \(A\) 为任意可观测量，\(n_1\), \(n_2\), \(n_3\) 为任意整数。
  • 周期性边界条件就是命令基本原胞完全等同地重复无穷多次。
  • 实现周期性边界条件的具体操作：当有一个粒子穿过基本原胞的六方体表面时；就让这个粒子也以相同的速度穿过此表面对面的表面，重新进入该原胞内。

3.4 截断距离

  • 出于计算上的考虑，力场必须截断，即力场在某一范围内有效。
  • 通常选取一个适当的力程截断距离(cut-off distance)将位势截断，以减少计算势能所耗用的时间。
  • 如果预先不采取特殊措施，一个分子动力学步，求和运算遍及所有相邻粒子，若只进行一次求和运算，则不存在任何困难。
  • 问题是这种求和运算将反复不断地进行多次，工作量就大得难以忍受，所需的总运行时间有可能 \(99％\) 是用来寻找相邻粒子、计算位势及作用在粒子上的力，时间浪费极大。

  • 势函数直接截断：

  • 一般分子势能函数在大于某个距离时，其值便迅速趋近于零。

  • 因此可以定义势函数 \[ V^s(r_{ij}) = \left\lbrace \begin{aligned}[l] V(r_{ij}) - V_c, ~~~ r_{ij} \leq r_c \\ 0, ~~~ r_{ij} > r_c \end{aligned} \right. \] 其中 \(V(r)\) 为一般势能，\(V_c\) 为截断势能，\(r_c\) 为截断半径。这样一来，便可减少许多不必要的计算量。

  • 对于不同分子动力学原胞内粒子间的相互作用，如果相互作用是短程力，可以在长度 \(r_c\) 处截断，显然 \(r_c\) 越大，所用时间就越多。

  • \(V_c\) 必须要足够小，以使截断不会显著地影响模拟结果。

  • 原胞尺度 \(L\) 的数值应选得很大，\(L\) 通常选得比 \(r_c\) 大很多。

  • 一般选择原胞尺度满足不等式条件 \(L/2>r_c\)，使得距离大于 \(L/2\) 的粒子的相互作用可以忽略，以避免有限尺寸效应。

3.5 最小像力约定

  • 在考虑粒子间的相互作用时，通常采用最小像力约定。

  • 最小像力约定是：设原胞内有 \(N\) 个粒子，在重复的分子动力学原胞中，每一个粒子只同它所在的原胞内的另外 \(N-1\) 个中的每个粒子或其最邻近的影像粒子发生相互作用。

  • 如 \(r_i\) 处的粒子 \(i\) 同 \(r_j\) 处的粒子 \(j\) 之间的距离为 \[ r_{ij} = \min\left(\left|\vec{r}_i - \vec{r}_j + \vec{n}L\right|\right) \]

  • 实际上这个约定就是通过满足不等式条件 \(r_c < L/2\) 来截断位势。

  • 采用最小像力约定后，原胞内第 \(i\) 个粒子与周围粒子的相互作用势和相互作用力为： \[ U_i(\vec{R}) = \sum_{j = 1, N}^{i \neq j} u(r_{ij}); ~~~ \vec{F}_i(\vec{R}) = \sum_{j = 1, N}^{i \neq j} F(r_{ij})\hat{r}_{ij} \] \(\vec{R} = \{\vec{r}_1, \vec{r}_2, ... \vec{r}_N,\}\)表示元胞内所有粒子的坐标。 \(\hat{r}_{ij}\) 是沿 \(r_j-r_i\) 方向的单位矢量。

  • 采用最小像力约定会使得在截断处粒子的受力有一个 \(\delta\) -函数的奇异性，这会给模拟计算带来误差。

  • 为减小这种误差，总可以将截断处粒子的相互作用势能换算成 \(V(r)- V_c\) ，以保证在截断处相互作用为零。

3.6 势函数

  • 分子动力学模拟的核心问题是求解牛顿运动方程组，而关键问题是原子间作用势的确定。
  • 通常是通过实验拟合或半经验解法得到原子间作用势，然后求得系统能量。
  • 固体间的性质，依赖于分子的本性与分子间的相互作用力，由统计力学的系统理论可知，只要知道决定分子间作用力的势能函数形态，就可以透过计算求得各种力学性质。

3.7 分子间的相互作用力

  • 以力的本性来区分有： 万有引力、强核力、弱核力、电磁力。
    • 其中万有引力太小，例如在典型的分子间几乎接触的距离约 \(0.2－0.3\) nm，对于能量的贡献约为 \(－10^{-30}\) J左右；
    • 强核力和弱核力的作用距离都在 \(10^{-4}\) nm以内，分子间的距离通常比此值还要大得多；
    • 只有电磁力能够正确解释分子间的相互作用。
  • 分子间作用力主要有：
    • 长程的引力，包括静电作用、诱导作用和色散作用；
    • 短程的斥力，它是由于电子云重叠所产生的。
  • 将分子间作用力的理论透过量子力学的计算与光谱等实验的方式，可求得分子间势能函数。
    • 两原子核的运动受一个只依赖于两原子核之间距离 \(r\) 的位势 \(V(r)\) 支配。
    • 在远距离处是吸引势（Van der Waals相互作用），而在近距离处是排斥势（原子核的Coulomb相互作用和电子的Pauli斥力）。
  • 势函数 \(U\) 看作是原子间相互作用势 \(V\) 和外势 \(U_e\) 之和 \[ U(r_i) = V(r_{ij}) + U_e(r_i) \]

3.8 势函数发展简介

  • 1924年，Lennard-Jones发表了著名的负幂函数式的Lennard-Jones势函数的解析式。
  • 1929年，Morse发表了指数式的Morse势。1931年，Born和Mayer发表了描述离子晶体中两离子间相互作用的Born-Mayer势函数。
  • 1983年至1984年间，Daw和Baskes提出了金属材料中的嵌入原子势(EAM)的概念和算法。
  • 几乎同时，Finnis和Sinclair根据密度函数二次矩理论提出形式上与EAM基本相同的经验F-S模型。

3.8.1 Lennard-Jones势

  • 早期的分子动力学主要用于模拟惰性气体和液体，发展了一些对势模型。
  • 其中应用非常普遍的是Lennard-Jones势。Lennard-Jones势由John Edward Lennard-Jones在1924年提出。
  • Lennard-Jones势函数表示为： \[ V(r_{ij}) = 4\varepsilon_{\alpha\beta} \left[\left(\frac{\sigma_{\alpha\beta}}{r_{ij}}\right)^{12} - \left(\frac{\sigma_{\alpha\beta}}{r_{ij}}\right)^{6}\right] \]
  • 式中，第一项代表短程排斥力项，这一项对 \(r_{ij}\) 的导数的负值为正，即为排斥，幂次高，在 \(r_{ij}\) 小的时候起主导作用；
  • 第二项代表远程吸引力项，这一项对 \(r_{ij}\) 的导数的负值为负，即为吸引，幕次低，在 \(r_{ij}\) 大的时候起主导作用。
    • 其中：\(-\varepsilon\) 为位势的最小值，这个最小值出现在距离 \(r\) 等于 \(2^{1/6}\sigma\) 的地方。\(\alpha\) 和 \(\beta\) 分别代表粒子 \(A\) 和 \(B\)。 \(r=\sigma\) 时位势为零。
    • 在分子动力学中，一般使用约化单位。当选取原子的质量 \(=1\) 时，原子的动量和速度以及力和加速度在数值上相等。
  • 在Lennard-Jones势作用下，第 \(i\) 个粒子与元胞内其它 \(N-1\) 个粒子或其最邻近的影像粒子的相互作用力在 \(x\) 方向的分量为： \[ F_{i,x} = 48\left(\frac{\varepsilon_{\alpha\beta}}{\sigma_{\alpha\beta}^2}\right)\sum_{j \neq i}^N(x_i - x_j)\left[\left(\frac{\sigma_{\alpha\beta}}{r_{ij}}\right)^{14}-\frac{1}{2}\left(\frac{\sigma_{\alpha\beta}}{r_{ij}}\right)^{8}\right] \]
  • Lennard-Jones势函数与力场函数曲线
Show the code 
x = np.linspace(0.6, 2.6, 100)
var('r')
var('V_r F_r', cls = Function)
vr = 4 * (r**(-12) - r**(-6))
fr = -vr.diff(r)
vrr = lambdify((r), vr)
frr = lambdify((r), fr)
Markdown(f'$${latex(Eq(V_r(r), vr))}$$')

\[V_{r}{\left(r \right)} = - \frac{4}{r^{6}} + \frac{4}{r^{12}}\]

Show the code 
Markdown(f'$${latex(Eq(F_r(r), fr))}$$')

\[F_{r}{\left(r \right)} = - \frac{24}{r^{7}} + \frac{48}{r^{13}}\]

Show the code 
plt.plot(x, vrr(x), lw = 4, label = r'$V(r)$')
plt.plot(x, frr(x), lw = 4, label = r'$F(r)$')
result = np.where(frr(x) >= 0)
xx = x[result[0][-1]]
yy = np.array([-4, 5])
plt.fill_betweenx(yy, 0, xx * np.ones(yy.shape), color = 'red', alpha = 0.3)
plt.fill_betweenx(yy, xx * np.ones(yy.shape), 2.6 * np.ones(yy.shape), color = 'blue', alpha = 0.3)
plt.plot(x, np.zeros(x.shape), 'k-.')
plt.ylim([-4, 5])
plt.xlim([0.6, 2.6])
plt.legend()
plt.text(0.62, -2, 'Repulsive force', fontsize = 14)
plt.text(1.6, 2, 'Attractive force', fontsize = 14)
plt.xlabel(r'$r_{ij}$', fontsize = 14)
plt.show()

3.9 分子动力学模拟的基本步骤

  • 平衡态分子动力学模拟的实际步骤可以划分为四步：
    1. 设定模拟所采用的模型
    • 模型的设定，就是选取原胞与势函数。势函数的研究和物理系统上对物质的描述研究息息相关。
    1. 给定初始条件
    • 运动方程的求解需要知道粒子的初始位置和速度，不同的算法要求不同的初始条件。
    • 如：速度verlet算法需要零时刻的坐标以及零时刻的速度值。
    1. 趋于平衡的计算过程
    • 为使得系统平衡，模拟中设计一个趋衡过程，即在这个过程中，增加或者从系统中移出能量，直到持续给出确定的能量值。
    • 称这时的系统已经达到平衡。
    1. 宏观物理量的计算
    • 实际计算宏观物理量往往是在分子动力学模拟的最后阶段进行的。
    • 它是沿着相空间轨迹求平均来计算得到的。

3.9.1 设定模拟模型

  • 设定模型是分子动力学模拟的第一步工作。
  • 例如在一个分子系统中，假定两个分子间的相互作用势为硬球势，其势函数表示为 \[ U(r) = \left\lbrace \begin{aligned}[l] +\infty, ~~~ (r < \sigma) \\ 0, ~~~~ (r \geq \sigma) \end{aligned} \right. \]
  • 常用的是Lennard-Jones 势函数、EAM势函数等，不同的物质状态描述用不同的势函数。
  • 为了计算方便，取分子动力学元胞为立方形，元胞的线度大小为 \(L\)，则其体积为 \(L^3\)。
  • 根据给定密度 \(\rho\) 和指定的单个元胞中的粒子数 \(N\)，确定出元胞的 \(L\) 值。采用周期性边界条件和最小像力约定。
  • 根据实际需要解决的问题，构造初始原胞中粒子物理结构模型。
  • 模型势函数一旦确定，就可以根据经典物理学规律求得模拟中的守恒量。例如对在微正则系综的模拟中能量、动量和角动量均为守恒量。

3.9.2 给定初始条件

  • 分子动力学模拟的过程进入对系统微分方程组做数值求解时，需要知道粒子的初始位置和速度的数值。

  • 不同的算法要求不同的初始条件。

  • 一般讲，系统的初始条件不可能知道，实际上也不需要精确选择代求系统的初始条件，因为模拟时间足够长时，系统就会忘掉初始条件。

  • 当然，合理的初始条件可以加快系统趋于平衡的时间和步伐，获得好的精度。

  • 常用的初始条件可以选择为：

    • 令初始位置在差分划分网格的格子上，初始速度则从玻尔兹曼分布随机抽样得到；
    • 令初始位置随机的偏离差分划分网格的格子上，初始速度为零；
    • 令初始位置随机的偏离差分划分网格的格子上，初始速度也是从玻尔兹曼分布随机抽样得到。 在我们后面的程序中，选择了第一种初始条件。

3.9.3 趋于平衡计算

  • 按照给出的运动方程、边界条件和初始条件，就可以进行分子动力学模拟计算。

  • 但是，这样计算出的系统不会具有所要求的系统能量，并且这个状态本身也还不是一个平衡态。

  • 为了使系统达到平衡，模拟中需要一个趋衡过程。在这个过程中，增加或从系统中移出能量，直到系统具有所要求的能量。

  • 然后，再对运动方程中的时间向前积分若干步，使系统持续给出确定能量值。 称这时系统已经达到平衡态。这段达到平衡所需的时间称为弛豫时间。

  • 称这时系统已经达到平衡态。这段达到平衡所需的时间称为弛豫时间。

  • 分子动力学计算的基本思想是：赋予分子体系初始运动状态之后，利用分子的自然运动在相空间中抽取样本进行统计分子动力学计算，时间步长就是抽样的间隔，因而时间步长的选取对动力学模拟非常重要。它决定了模拟所需要的时间。

  • 时间步长 \(h\) 太长会造成分子间激烈碰撞，体系数据溢出；

  • 为了减小误差，步长 \(h\) 必须取得小一些；

  • 但是取得太小，系统模拟的弛豫时间就很长，太短的时间步长还会降低模拟过程搜索相空间的能力。

  • 因此一般选取的时间步长为体系各个自由度中最短运动周期的十分之一。

  • 这里需要积累一定的模拟经验，选择适当的时间步长 \(h\) 。

3.9.4 宏观物理量的计算

  • 在模拟过程中计算出的动能值是在不连续的路径上的值，在时间的各个间断点 \(\mu\) 上计算动能的平均值 \[ \bar{E}_k = \frac{1}{n-n_0} \sigma_{\mu > n_0}^n \sigma_{i = 1}^N \frac{(p_i^2)^\mu}{2m} \] 在分子动力学模拟过程中，温度是需要加以监测的物理量，特别是在模拟的起始阶段。根据能量均分定理，我们可以从平均动能值计算得到温度值， \[ T = \frac{\bar{E}_k}{\frac{d}{2}Nk_B} \] 其中 \(d\) 为每个粒子的自由度，如果不考虑系统所受的约束，则 \(d＝3\)。

3.9.5 平衡态分子动力学模拟

  • 在经典分子动力学模拟方法的应用当中，存在着对两种系统状态的分子动力学模拟。
    • 一种是对平衡态的MD模拟；
    • 另一类是对非平衡态的分子动力学模拟。
  • 对平衡态系综分子动力学模拟又可以分为如下类型：
    • 微正则系综的分子动力学 \((N,V,E)\) 模拟
    • 正则系综的分子动力学 \((N,V,T)\) 模拟
    • 等温等压系综分子动力学 \((N,P,T)\) 模拟
    • 等焓等压系综分子动力学 \((N,P,H)\) 模拟
    • 巨正则系综分子动力学 \((V,T, \mu)\) 模拟
  • 这里仅对平衡态的分子动力学方法中几类模拟做介绍
3.9.5.1 系综
  • 系综(ensemble) 是指在一定的宏观条件下(约束条件)，大量性质和结构完全相同的、处于各种运动状态的、各自独立的系统的集合。全称为统计系综。
  • 系综是用统计方法描述热力学系统的统计规律性时引入的一个基本概念；
  • 系综是统计理论的一种表述方式，系综理论使统计物理成为普遍的微观统计理论 ；
  • 系综并不是实际的物体，构成系综的系统才是实际物体。
  • 约束条件是由一组外加宏观参量来表示。在平衡统计力学范畴下，可以用来处理稳定系综。
3.9.5.2 常用系综分类
  • 根据宏观约束条件，系综被分为以下几种：
    1. 正则系综(canonical ensemble)，全称应为“宏观正则系综”，简写为\(NVT\)，即表示具有确定的粒子数 \((N)\)、体积 \((V)\)、温度 \((T)\)。 正则系综是蒙特卡罗方法模拟处理的典型代表。假定 \(N\) 个粒子处在体积为 \(V\) 的盒子内，将其埋入温度恒为 \(T\) 的热浴中。此时，总能量 \((E)\) 和系统压强 \((P)\) 可能在某一平均值附近起伏变化。 平衡体系代表封闭系统，与大热源热接触平衡的恒温系统。正则系综的特征函数是亥姆霍兹自由能 \(F(N,V,T)\)。
    2. 微正则系综(micro-canonical ensemble)，简写为NVE，即表示具有确定的粒子数 \((N)\)、体积 \((V)\)、总能量 \((E)\)。微正则系综广泛被应用在分子动力学模拟中。假定 \(N\) 个粒子处在体积为 \(V\) 的盒子内，并固定总能量 \((E)\)。此时，系综的温度 \((T)\) 和系统压强 \((P)\) 可能在某一平均值附近起伏变化。平衡体系为孤立系统，与外界即无能量交换，也无粒子交换。微正则系综的特征函数是熵 \(S(N,V,E)\)。
    3. 等温等压(constant-pressure,constant-temperature)，简写为NPT，即表示具有确定的粒子数 \((N)\)、压强 \((P)\)、温度 \((T)\)。一般是在蒙特卡罗模拟中实现。其总能量 \((E)\) 和系统体积 \((V)\) 可能存在起伏。体系是可移动系统壁情况下的恒温热浴。特征函数是吉布斯自由能 \(G(N,P,T)\)。
    4. 等压等焓(contant-pressure,constant- enthalpy)，简写为\(NPH\)，即表示具有确定的粒子数 \((N)\)、压强 \((P)\)、焓 \((H)\)。由于由于 \(H =E+PV\)，故在该系综下进行模拟时要保持压力与焓值为固定，其调节技术的实现也有一定的难度，这种系综在实际的分子动力学模拟中已经很少遇到了。
    5. 巨正则系综(grand canonical ensemble)，简写为 \(VT\mu\)，即表示具有确定的粒体积 \((V)\)、温度 \((T)\) 和化学势 \((\mu)\)。巨正则系综通常是蒙特卡罗模拟的对象和手段。此时、系统能量 \((E)\)、压强 \((P)\) 和粒子数 \((N)\) 会在某一平均值附近有一个起伏。体系是一个开放系统，与大热源大粒子源热接触平衡而具有恒定的 \(T\)。特征函数是马休(Massieu)函数 \(J(\mu，V，T)\)。

3.9.6 微正则系综的分子动力学模拟

  • 粒子数恒定、体积恒定、能量恒定、整个系统的总动量恒等于零。 分子动力学模拟步骤如下（Verlet算法）：

    1. 给定初始空间位置 \(\{r_i(0)\}\) 。
    2. 计算在第n步时粒子所受的力\(\{F_i(n)\}, F_i(n)=F_i(t_n)\)。
    3. 利用如下公式，计算在第 \(n+1\) 步时所有粒子所处的空间位置。 \[ \vec{r}_i^{n+1} = 2\vec{r}_i^n - \vec{r}_i^{n-1} + \vec{F}_i^n h^2/m \]
    4. 计算第 \(n\) 步的速度 \[ \vec{v}_i^n = (\vec{r}_i^{n+1} - \vec{r}_i^{n-1})/2h \]
    5. 返回到步骤 ii)，开始下一步的模拟计算。

  • Verlet算法的速度形式：

    1. 给定初始空间位置 \(\{\vec{r}_i^1\}\) 。
    2. 给定初始速度 \(\{\vec{v}_i^1\}\) 。
    3. 利用公式： \[ \vec{v}_i^{n+1} = \vec{r}_i^n + h\vec{v}_i^n + \vec{F}_i^n h^2 / 2m \] 计算在第 \(n+1\) 步时所有粒子所处的空间位置 \(\{ \vec{r}_i^{n+1} \}\) 。
    4. 计算在第 \(n+1\) 步时所有粒子的速度： \[ \vec{v}_i^{n+1} = \vec{v}_i^n + h(\vec{F}_i^{n+1} + \vec{F}_i^n) / 2m \]
    5. 返回到步骤 iii) ，开始第 \(n+2\) 步的模拟计算。

  • Verlet速度形式的算法比前一种算法好些。

  • 它不仅可以在计算中得到同一时间步上的空间位置和速度，而且数值计算的稳定性也提高了。

  • 一般情况下，对于能量的确定的系统不可能给出精确的初始条件。

  • 这时需要先给出一个合理的初始条件，然后在模拟过程中逐渐调节系统能量达到给定值。

  • 其步骤为：

    • 首先将运动方程组解出若干步的结果；
    • 然后计算出动能和位能；
    • 假如总能量不等于给定恒定值，则通过调整速度来实现总能量不变的目的，具体做法就是将速度乘以一个标度(scaling)因子 \(\beta\)，即： \[ \{ \vec{v}_i^{n+1} \beta \} \rightarrow \{ \vec{v}_i^{n+1} \} \]

  • 速度标度因子 \(\beta\) : \[ \beta = \left[ \frac{T^*(N - 1)}{16\sum_i v_i^2} \right]^{1/2} \]

Back to top