CudaSPONGE参考手册
1.前言及免责声明
本手册参考了GROMACS的参考手册[1.1] 进行组织，但并非对是其的中文翻译，而是根据其结构讲述CudsSPONGE中的各种设置——实际上对各种MD模拟软件的代码进行阅读以后就能够发现，几乎每一个MD软件都会对文献中的算法进行适当调整，而CudaSPONGE中也包含了这样的操作，也是本参考手册的意义所在。
本手册并不完整，并且由于撰稿人时间不足——我们的首要任务是改进SPONGE——本手册并不奢求完整。本手册还在不断完善中，信息可能也会不完全正确。
欢迎对形式和内容提出意见，可以直接在本页面的讨论中对相关信息进行补充。
2.SPONGE
分子模拟是一个非常有意思的词语，这个汉语词汇可以是英语中两个词组的翻译，分别是molecular simulation和molecular modeling。虽然这两个词组也分别有可以互相区分的翻译，例如分别翻译为分子仿真和分子建模，但在汉语语境下很多场合仍然是用分子模拟进行翻译。Molecular modeling泛指通过各种分子模型来描述复杂化学体系的一般过程，目的是根据原子尺度上的详细知识来理解和预测宏观特性，而molecular simulation则通常指代的是以分子动力学模拟和蒙特卡洛模拟为代表的具体的利用原子进行模拟的方法。
SPONGE全名为“Simulation Package tOwards Next GEneration molecular modeling”，也即面向下一代分子模拟的模拟软件，其中包含两个模拟，但是它们的意义就如上所述存在区别。从名字中也可以看出，SPONGE的中心语是Simulation Package，也即模拟软件，它的核心主要还是集中在molecular simulation，但是目的则是称为下一代分子模拟工具，也即努力方向是更全面的、更先进的molecular modeling。就目前而言，SPONGE主要还是作为一个分子动力学模拟工具，并包括一些附带的分子模拟功能。
3.数学基础定义
3.1 矩阵和矢量定义
在本文中，涉及到诸多物理量，其中矩阵和矢量均使用黑体书写，如力矢量F \mathbf{F} F Π \mathbf{\Pi} Π F F F
在未加声明的情况下，矢量均为行矢量，例如某个粒子的力矢量形状为( 1 , 3 ) (1,3) ( 1 , 3 )
F = ( F x F y F z ) (3-1) \mathbf{F} =
\left(
\begin{array}{ccc}
F_x & F_y & F_z
\end{array}
\right)
\tag{3-1} F = ( F x F y F z ) ( 3-1 ) 在CudaSPONGE的程序中，使用数据结构体VECTOR表示对应的矢量。
一个形状为( m , n ) (m,n) ( m , n ) Y \mathbf{Y} Y x x x ∂ Y ∂ x \frac{\partial\mathbf{Y}}{\partial x} ∂ x ∂ Y ( m , n ) (m,n) ( m , n )
∂ Y ∂ x = ( ∂ y 11 ∂ x ∂ y 12 ∂ x ⋯ ∂ y 1 n ∂ x ∂ y 21 ∂ x ∂ y 22 ∂ x ⋯ ∂ y 2 n ∂ x ⋯ ⋯ ⋱ ⋯ ∂ y m 1 ∂ x ∂ y m 1 ∂ x ⋯ ∂ y m n ∂ x ) (3-2) \frac{\partial\mathbf{Y}}{\partial x} =
\left(
\begin{array}{cccc}
\frac{\partial y_{11}}{\partial x} & \frac{\partial y_{12}}{\partial x} & \cdots & \frac{\partial y_{1n}}{\partial x} \\
\frac{\partial y_{21}}{\partial x} & \frac{\partial y_{22}}{\partial x} & \cdots & \frac{\partial y_{2n}}{\partial x} \\
\cdots & \cdots & \ddots & \cdots \\
\frac{\partial y_{m1}}{\partial x} & \frac{\partial y_{m1}}{\partial x} & \cdots & \frac{\partial y_{mn}}{\partial x} \\
\end{array}
\right)
\tag{3-2} ∂ x ∂ Y = ∂ x ∂ y 11 ∂ x ∂ y 21 ⋯ ∂ x ∂ y m 1 ∂ x ∂ y 12 ∂ x ∂ y 22 ⋯ ∂ x ∂ y m 1 ⋯ ⋯ ⋱ ⋯ ∂ x ∂ y 1 n ∂ x ∂ y 2 n ⋯ ∂ x ∂ y mn ( 3-2 ) 一个形状为( p , q ) (p,q) ( p , q ) Y \mathbf{Y} Y ( m , n ) (m,n) ( m , n ) X \mathbf{X} X ∂ Y ∂ X \frac{\partial\mathbf{Y}}{\partial\mathbf{X}} ∂ X ∂ Y ( m n , p q ) (mn,pq) ( mn , pq )
∂ Y ∂ X = ( ∂ y 11 ∂ x 11 ∂ y 21 ∂ x 11 ⋯ ∂ y p 1 ∂ x 11 ⋯ ∂ y p q ∂ x 11 ∂ y 11 ∂ x 21 ∂ y 21 ∂ x 21 ⋯ ∂ y p 1 ∂ x 21 ⋯ ∂ y p q ∂ x 21 ⋯ ⋯ ⋱ ⋯ ⋱ ⋯ ∂ y 11 ∂ x m n ∂ y 21 ∂ x m n ⋯ ∂ y p 1 ∂ x m n ⋯ ∂ y p q ∂ x m n ) (3-3) \frac{\partial \mathbf{Y}}{\partial \mathbf{X}} =
\left(
\begin{array}{cccccc}
\frac{\partial y_{11}}{\partial x_{11}} & \frac{\partial y_{21}}{\partial x_{11}} & \cdots & \frac{\partial y_{p1}}{\partial x_{11}} & \cdots & \frac{\partial y_{pq}}{\partial x_{11}} \\
\frac{\partial y_{11}}{\partial x_{21}} & \frac{\partial y_{21}}{\partial x_{21}} & \cdots & \frac{\partial y_{p1}}{\partial x_{21}} & \cdots & \frac{\partial y_{pq}}{\partial x_{21}} \\
\cdots & \cdots & \ddots & \cdots &\ddots & \cdots \\
\frac{\partial y_{11}}{\partial x_{mn}} & \frac{\partial y_{21}}{\partial x_{mn}} & \cdots & \frac{\partial y_{p1}}{\partial x_{mn}} & \cdots & \frac{\partial y_{pq}}{\partial x_{mn}} \\
\end{array}
\right)
\tag{3-3} ∂ X ∂ Y = ∂ x 11 ∂ y 11 ∂ x 21 ∂ y 11 ⋯ ∂ x mn ∂ y 11 ∂ x 11 ∂ y 21 ∂ x 21 ∂ y 21 ⋯ ∂ x mn ∂ y 21 ⋯ ⋯ ⋱ ⋯ ∂ x 11 ∂ y p 1 ∂ x 21 ∂ y p 1 ⋯ ∂ x mn ∂ y p 1 ⋯ ⋯ ⋱ ⋯ ∂ x 11 ∂ y pq ∂ x 21 ∂ y pq ⋯ ∂ x mn ∂ y pq ( 3-3 ) 这种先按分母展开的导数定义被称为分母布局。按分子展开的导数定义为分子布局，两者互为倒置。
一个形状为( m , n ) (m,n) ( m , n ) X \mathbf{X} X ( n , l ) (n,l) ( n , l ) Y \mathbf{Y} Y X Y \mathbf{X}\mathbf{Y} XY ( m , l ) (m,l) ( m , l )
X Y = ( ∑ i n x 1 i y i 1 ∑ i n x 1 i y i 2 ⋯ ∑ i n x 1 i y i l ∑ i n x 2 i y i 1 ∑ i n x 2 i y i 2 ⋯ ∑ i n x 2 i y i l ⋯ ⋯ ⋱ ⋯ ∑ i n x m i y m 1 ∑ i n x m i y i 2 ⋯ ∑ i n x m i y i l ) (3-4) \mathbf{X}\mathbf{Y} =
\left(
\begin{array}{cccc}
\sum_i^n x_{1i}y_{i1} & \sum_i^n x_{1i}y_{i2} & \cdots & \sum_i^n x_{1i}y_{il} \\
\sum_i^n x_{2i}y_{i1} & \sum_i^n x_{2i}y_{i2} & \cdots & \sum_i^n x_{2i}y_{il} \\
\cdots & \cdots & \ddots & \cdots \\
\sum_i^n x_{mi}y_{m1} & \sum_i^n x_{mi}y_{i2} & \cdots & \sum_i^n x_{mi}y_{il} \\
\end{array}
\right)
\tag{3-4} XY = ∑ i n x 1 i y i 1 ∑ i n x 2 i y i 1 ⋯ ∑ i n x mi y m 1 ∑ i n x 1 i y i 2 ∑ i n x 2 i y i 2 ⋯ ∑ i n x mi y i 2 ⋯ ⋯ ⋱ ⋯ ∑ i n x 1 i y i l ∑ i n x 2 i y i l ⋯ ∑ i n x mi y i l ( 3-4 ) 一个形状为( m , n ) (m,n) ( m , n ) X \mathbf{X} X ( m , n ) (m,n) ( m , n ) Y \mathbf{Y} Y t r ( X T Y ) {\rm tr}(\mathbf{X}^{\rm T}\mathbf{Y}) tr ( X T Y )
t r ( X T Y ) = ∑ i m ∑ j n x i j y i j (3-5) {\rm tr}(\mathbf{X}^{\rm T}\mathbf{Y}) = \sum_i^m\sum_j^n x_{ij}y_{ij}
\tag{3-5} tr ( X T Y ) = i ∑ m j ∑ n x ij y ij ( 3-5 ) 一个形状为( m , n ) (m,n) ( m , n ) X \mathbf{X} X ( p , q ) (p,q) ( p , q ) Y \mathbf{Y} Y X ⊗ Y \mathbf{X} \otimes \mathbf{Y} X ⊗ Y ( m p , q n ) (mp,qn) ( m p , q n )
X ⊗ Y = ( x 11 Y x 12 Y ⋯ x 1 n Y x 21 Y x 22 Y ⋯ x 2 n Y ⋯ ⋯ ⋱ ⋯ x m 1 Y x m 2 Y ⋯ x m n Y ) (3-6) \mathbf{X} \otimes \mathbf{Y} =
\left(
\begin{array}{cccc}
x_{11}\mathbf{Y} & x_{12}\mathbf{Y} & \cdots & x_{1n}\mathbf{Y} \\
x_{21}\mathbf{Y} & x_{22}\mathbf{Y} & \cdots & x_{2n}\mathbf{Y} \\
\cdots & \cdots & \ddots & \cdots \\
x_{m1}\mathbf{Y} & x_{m2}\mathbf{Y} & \cdots & x_{mn}\mathbf{Y} \\
\end{array}
\right)
\tag{3-6} X ⊗ Y = x 11 Y x 21 Y ⋯ x m 1 Y x 12 Y x 22 Y ⋯ x m 2 Y ⋯ ⋯ ⋱ ⋯ x 1 n Y x 2 n Y ⋯ x mn Y ( 3-6 ) 由上，有一些重要推论，包括求导的链式法则：
∂ Z ∂ X = ∂ Y ∂ X ∂ Z ∂ Y (3-7) \frac{\partial \mathbf{Z}}{\partial \mathbf{X}} =
\frac{\partial \mathbf{Y}}{\partial \mathbf{X}}
\frac{\partial \mathbf{Z}}{\partial \mathbf{Y}}
\tag{3-7} ∂ X ∂ Z = ∂ X ∂ Y ∂ Y ∂ Z ( 3-7 ) 对于任意Y = A X B \mathbf{Y}=\mathbf{A}\mathbf{X}\mathbf{B} Y = AXB
∂ Y ∂ X = B ⊗ A T (3-8) \frac{\partial \mathbf{Y}}{\partial \mathbf{X}} =
\mathbf{B} \otimes \mathbf{A}^{\rm T}
\tag{3-8} ∂ X ∂ Y = B ⊗ A T ( 3-8 ) 3.2 单位制
CudaSPONGE在程序内部使用一套自身的单位制，与程序外部的数据交换的单位根据情况进行设定。
内部的单位使用为：
物理量 单位 解释 长度 Å 10 − 10 m \rm 10^{-10}\ m 1 0 − 10 m 时间 20.455 fs 20.455 × 10 − 12 s \rm 20.455 \times 10^{-12}\ s 20.455 × 1 0 − 12 s 速度 Å/(20.455 fs) 长度单位除以时间单位 电荷 18.2223 e e为元电荷 质量 Da 碳12原子质量的1/12 温度 K 热力学温标 能量 kcal/mol 1 k c a l / m o l = 4.184 k J / m o l \rm 1\ kcal/mol = 4.184\ kJ/mol 1 kcal/mol = 4.184 kJ/mol 力 kcal/mol/Å 压强 1.4395 × 10 − 5 b a r \rm 1.4395\times10^{-5}\ bar 1.4395 × 1 0 − 5 bar N A / ( 4.184 × 10 28 ) N_{\rm A}/( 4.184\times10^{28}) N A / ( 4.184 × 1 0 28 )
其中，
电荷选择该单位的原因是希望在该单位制下，库伦定律中的系数为1，也即
E C o u l o m b = q 1 q 2 r 12 (3-9) E_{\rm Coulomb} = \frac{q_1q_2}{r_{12}} \tag{3-9} E Coulomb = r 12 q 1 q 2 ( 3-9 ) 时间选择该单位的原因是希望在该单位制下，动能计算中的系数为1 2 \frac{1}{2} 2 1
E k i n e t i c = 1 2 m v 2 (3-10) E_{\rm kinetic} = \frac{1}{2}mv^2 \tag{3-10} E kinetic = 2 1 m v 2 ( 3-10 ) 压强选择该单位的原因是希望在该单位制下，体积功的计算系数为1，也即
W = P V (3-11) W = PV \tag{3-11} W = P V ( 3-11 ) 4.分子动力学
分子动力学的基础是经典力学，在这里通过牛顿力学和哈密顿力学进行简单表述。
4.1 牛顿力学
对于任意体系，我们可以计算任意时刻的力F \mathbf{F} F v \mathbf{v} v r \mathbf{r} r
{ d v d t = F m d r d t = v (4-1) \left\{
\begin{array}{l}
\frac{d\mathbf{v}}{dt}=\frac{\mathbf{F}}{m} \\
\frac{d\mathbf{r}}{dt}=\mathbf{v}
\end{array}
\right.
\tag{4-1} { d t d v = m F d t d r = v ( 4-1 ) 将该微分方程组对时间t t t
{ v ( t + 0.5 Δ t ) = v ( t − 0.5 Δ t ) + F m Δ t r ( t + Δ t ) = r ( t ) + v ( t + 0.5 Δ t ) Δ t (4-2) \left\{
\begin{array}{l}
\mathbf{v}(t+0.5\Delta t)=\mathbf{v}(t-0.5\Delta t) + \frac{\mathbf{F}}{m} \Delta t\\
\mathbf{r}(t+\Delta t)=\mathbf{r}(t) + \mathbf{v}(t+0.5\Delta t)\Delta t
\end{array}
\right.
\tag{4-2} { v ( t + 0.5Δ t ) = v ( t − 0.5Δ t ) + m F Δ t r ( t + Δ t ) = r ( t ) + v ( t + 0.5Δ t ) Δ t ( 4-2 ) 4.2 哈密顿力学
在哈密顿力学框架下，使用正则坐标{ q } \{q\} { q } { p } \{p\} { p } H \mathscr{H} H T ( { p } ) T(\{p\}) T ({ p }) U ( { q } ) U(\{q\}) U ({ q })
H = T + U (4-3) \mathscr{H}=T + U \tag{4-3} H = T + U ( 4-3 ) 有哈密顿正则方程：
{ d p i d t = − ∂ H ∂ q i d q i d t = ∂ H ∂ p i (4-4) \left\{
\begin{array}{l}
\frac{dp_i}{dt} = -\frac {\partial \mathscr{H}}{\partial q_i} \\
\frac{dq_i}{dt} = \frac {\partial \mathscr{H}}{\partial p_i} \\
\end{array}
\right.
\tag{4-4} { d t d p i = − ∂ q i ∂ H d t d q i = ∂ p i ∂ H ( 4-4 ) 不难看出，实际上式(4-1)和式(4-4)是完全等价的，只是写法不同。
相空间是包含了体系所有正则坐标和正则动量的空间。在某时刻t t t d q d p d\mathbf{q}d\mathbf{p} d q d p ρ ( q , p , t ) \rho(\mathbf{q},\mathbf{p},t) ρ ( q , p , t )
在分子动力学模拟中，一个关键的方程是刘维尔方程，它认为分布函数沿着相空间的任何轨迹是常数：
d ρ ( q , p , t ) d t = 0 (4-5) \frac{d\rho(\mathbf{q},\mathbf{p},t)}{dt} = 0 \tag{4-5} d t d ρ ( q , p , t ) = 0 ( 4-5 ) 对式(4-5)展开，可得：
∂ ρ ∂ t + ∑ i ( ∂ ρ ∂ q i d q i d t + ∂ ρ ∂ p i d p i d t ) = 0 (4-6) \frac {\partial \rho}{\partial t} + \sum_i(\frac {\partial \rho}{\partial q_i} \frac {dq_i}{dt}+ \frac {\partial \rho}{\partial p_i}\frac {dp_i}{dt}) = 0 \tag{4-6} ∂ t ∂ ρ + i ∑ ( ∂ q i ∂ ρ d t d q i + ∂ p i ∂ ρ d t d p i ) = 0 ( 4-6 ) 利用哈密顿正则方程式(4-4)替换式(4-6)中的时间全微分项，并定义刘维尔算符
i L ^ = ∑ i ( ∂ H ∂ p i ∂ ∂ q i − ∂ H ∂ q i ∂ ∂ p i ) (4-7) i\hat L= \sum_i\Big(\frac{\partial \mathscr H}{\partial p_i} \frac{\partial}{\partial q_i} - \frac{\partial \mathscr H}{\partial q_i} \frac{\partial}{\partial p_i}\Big)
\tag{4-7} i L ^ = i ∑ ( ∂ p i ∂ H ∂ q i ∂ − ∂ q i ∂ H ∂ p i ∂ ) ( 4-7 ) 可得刘维尔方程：
∂ ρ ∂ t = − i L ^ ρ (4-8) \frac{\partial \rho} {\partial t} = -i\hat L \rho
\tag{4-8} ∂ t ∂ ρ = − i L ^ ρ ( 4-8 ) 该方程非常类似含时薛定谔方程，只相差一个系数ℏ \hbar ℏ
i ℏ ∂ φ ∂ t = H ^ φ (4-9) i\hbar \frac{\partial \varphi} {\partial t} = \hat H \varphi
\tag{4-9} i ℏ ∂ t ∂ φ = H ^ φ ( 4-9 ) 求解式(4-8)所示的微分方程，有
ρ ( q , p , t ) = e x p ( − i L ^ t ) ρ ( q , p , 0 ) (4-10) \rho(\mathbf{q},\mathbf{p},t) = {\rm exp}(-i\hat Lt)\rho(\mathbf{q},\mathbf{p},0) \tag{4-10} ρ ( q , p , t ) = exp ( − i L ^ t ) ρ ( q , p , 0 ) ( 4-10 ) 对于该指数形式，有Trotters-Suzuki分解：
e x p ( − i ( A ^ + B ^ ) t ) ≈ e x p ( − i A ^ t ) e x p ( − i B ^ t ) (4-11) {\rm exp}(-i(\hat A + \hat B)t) \approx {\rm exp}(-i\hat At) {\rm exp}(-i\hat Bt) \tag{4-11} exp ( − i ( A ^ + B ^ ) t ) ≈ exp ( − i A ^ t ) exp ( − i B ^ t ) ( 4-11 ) Trotters-Suzuki分解非常重要，首先是可以利用t = t Δ t Δ t t =\frac{t}{\Delta t}{\Delta t} t = Δ t t Δ t ( e x p ( − i L ^ Δ t ) ) t Δ t \big({\rm exp}(-i\hat L\Delta t)\big)^{\frac{t}{\Delta t}} ( exp ( − i L ^ Δ t ) ) Δ t t H \mathscr H H
5.周期性边界条件（PBC）
5.1 三维周期性边界条件
大部分模拟中使用三维周期性边界条件，也即认为模拟只包括一个晶胞，而该体系是该晶胞在空间中无限延展的。下图是二维PBC的示例，该图片源自GROMACS手册。
l x l_x l x l y l_y l y l z l_z l z α \alpha α β \beta β γ \gamma γ b \mathbf{b} b l x \mathbf{l}_x l x x x x l y \mathbf{l}_y l y x y xy x y l z \mathbf{l}_z l z b \mathbf{b} b
b = ( l x l y l z ) = ( b x x 0 0 b y x b y y 0 b z x b z y b z z ) (5-1) \mathbf{b} =
\left(
\begin{array}{l}
\mathbf{l}_x \\
\mathbf{l}_y \\
\mathbf{l}_z
\end{array}
\right)
=
\left(
\begin{array}{lll}
b_{xx} & 0 & 0 \\
b_{yx} & b_{yy} & 0 \\
b_{zx} & b_{zy} & b_{zz}
\end{array}
\right)
\tag{5-1} b = l x l y l z = b xx b y x b z x 0 b y y b z y 0 0 b z z ( 5-1 ) 体积V = d e t ( b ) = b x x b y y b z z V = {\rm det}(\mathbf{b}) = b_{xx}b_{yy}b_{zz} V = det ( b ) = b xx b y y b z z b − 1 \mathbf{b}^{-1} b − 1
此时，任意位置和位移的分数坐标s \mathbf{s} s r \mathbf{r} r
{ s = r b − 1 r = s b (5-2) \left\{
\begin{array}{l}
\mathbf{s} = \mathbf{r} \mathbf{b}^{-1}\\
\mathbf{r} = \mathbf{s} \mathbf{b}
\end{array}
\right.
\tag{5-2} { s = r b − 1 r = sb ( 5-2 ) 在周期性边界下i i i j j j r i j \mathbf{r}_{ij} r ij
r i j = r i − r j − f l o o r ( ( r i − r j ) b − 1 + C ) b (5-3) \mathbf{r}_{ij} = \mathbf{r}_i - \mathbf{r}_j - {\rm floor}((\mathbf{r}_i - \mathbf{r}_j)\mathbf{b}^{-1}+\mathbf{C})\mathbf{b}
\tag{5-3} r ij = r i − r j − floor (( r i − r j ) b − 1 + C ) b ( 5-3 ) 其中，C = ( 0.5 0.5 0.5 ) \mathbf{C} = (0.5\ 0.5\ 0.5) C = ( 0.5 0.5 0.5 ) f l o o r {\rm floor} floor
模拟中，通常会要求对粒子求和时只计算截断距离r c r_c r c r c r_c r c
{ l x > 2 ∗ r c l y > 2 ∗ r c l z > 2 ∗ r c (5-4) \left\{
\begin{array}{l}
l_x > 2 * r_c \\
l_y > 2 * r_c \\
l_z > 2 * r_c
\end{array}
\right.
\tag{5-4} ⎩ ⎨ ⎧ l x > 2 ∗ r c l y > 2 ∗ r c l z > 2 ∗ r c ( 5-4 ) 粒子的笛卡尔坐标平移任意个周期对于体系是等效的，本身不会影响模拟，只会影响输出轨迹中的可视化。CudaSPONGE在储存粒子坐标时，总会将其储存在距离原点最近的正数值，同时保证分子的完整性。
5.2 其他边界条件
CudaSPONGE的程序中的实现中本质只包含了三维周期性边界条件，但是可以将其改变视为其他周期性边界条件。
对于非周期性边界条件下的模拟，CudaSPONGE中会额外保存一套不进行周期性映射的坐标，用于计算非键相互作用。而对于成键相互作用仍然使用周期性边界条件下的条件，因此在输入时需要将盒子边长设大，防止出错。
对于伪2维周期性边界条件，非周期性维度只能是z轴。对于z轴方向，非键近程作用和成键作用处理仍然与周期性边界条件时相同，因此需要谨慎地设置z轴的边长。对于非键远程作用，提供额外的PMC-IZ算法[5.1] 。
6.约束
在这里，我们将constrain、constraint等翻译为约束，而restrain、restraint等翻译为限制。约束和限制的区别在这里定义为，约束是通过算法强制满足，而限制是添加能量惩罚项。
6.1 SETTLE
SETTLE是保持水形状的迭代解析解。默认情况下，如果使用约束算法，都将会对水使用SETTLE算法。SETTLE的实现在文献[6.1] 的基础上，并没有要求等腰三角形的情况，同时也添加了对速度和压强的修正。
6.2 SHAKE
SHAKE是通过迭代来计算约束。SHAKE的实现主要参考了文献[6.2] 和文献[6.3] 。
其主要思路为，对于更新前位移为r 0 \mathbf{r_0} r 0 d d d i i i r i \mathbf{r}_i r i F \mathbf{F} F
F = 1 2 r i 2 − d 2 r i r 0 r 0 (6-1) \mathbf{F}=\frac{1}{2}\frac{\mathbf{r}_i^2-d^2}{\mathbf{r}_i \mathbf{r}_0}\mathbf{r_0}
\tag{6-1} F = 2 1 r i r 0 r i 2 − d 2 r 0 ( 6-1 ) 如果需要对角度进行约束，CudaSPONGE能够将通过余弦定理转化为对距离的约束，但通常这种约束并不物理，难以收敛，且对步长的帮助很少，因此不推荐对角度进行约束。
7.控温
控温主要是对粒子的速度进行更改。上文中提到了Trotter-Suzuki分解，它指导我们能够将速度和位置先后迭代。同样地，对于控温，我们也可以将对应的粒子速度的修正迭代单独进行，但这种计算是有多种不同的方式。我们使用T ^ ( Δ t ) = e x p ( − i L ^ T Δ t ) \hat T(\Delta t) = {\rm exp}(-i\hat L_T\Delta t) T ^ ( Δ t ) = exp ( − i L ^ T Δ t ) r ^ ( Δ t ) = e x p ( − i L ^ r Δ t ) \hat r(\Delta t) = {\rm exp}(-i\hat L_r\Delta t) r ^ ( Δ t ) = exp ( − i L ^ r Δ t ) v ^ ( Δ t ) = e x p ( − i L ^ v Δ t ) \hat v(\Delta t) = {\rm exp}(-i\hat L_v\Delta t) v ^ ( Δ t ) = exp ( − i L ^ v Δ t )
e x p ( − i L ^ Δ t ) ≈ v ^ ( Δ t 2 ) T ^ ( Δ t ) r ^ ( Δ t ) v ^ ( Δ t 2 ) ≈ v ^ ( Δ t 2 ) r ^ ( Δ t 2 ) T ^ ( Δ t ) r ^ ( Δ t 2 ) v ^ ( Δ t 2 ) (7-1) \begin{array}{lll}
{\rm exp}(-i\hat L \Delta t) & \approx & \hat v(\frac {\Delta t}{2}) \hat T(\Delta t) \hat r(\Delta t)\hat v(\frac {\Delta t}{2}) \\
& \approx &\hat v(\frac {\Delta t}{2})\hat r(\frac {\Delta t}{2}) \hat T(\Delta t)\hat r(\frac {\Delta t}{2})\hat v(\frac {\Delta t}{2}) \\
\end{array}
\tag{7-1} exp ( − i L ^ Δ t ) ≈ ≈ v ^ ( 2 Δ t ) T ^ ( Δ t ) r ^ ( Δ t ) v ^ ( 2 Δ t ) v ^ ( 2 Δ t ) r ^ ( 2 Δ t ) T ^ ( Δ t ) r ^ ( 2 Δ t ) v ^ ( 2 Δ t ) ( 7-1 ) 其中第一种迭代方案被称为普通迭代方案，第二种方式被称为居中迭代方案，该方案的设计主要参考了文献[7.1] 。
7.1 Andersen控温
Andersen热浴通过周期性地令NVE 系综体系中的粒子按照Maxwell-Boltzmann分布重新分配速度，从而实现控温。在SPONGE的实现中，每经过设定的步长进行全部原子的速度重分布。速度重分布步长时不宜设置过小，否则会丧失动力学信息，速度时间自相关函数减弱。
7.2 Berendsen控温
Berendsen控温是一种变标度控温法。假设体系被放置在温度为T 0 T_0 T 0 T ( t ) T(t) T ( t )
d T d t = 1 τ [ T 0 − T ( t ) ] (7-2) \begin{array}{l}
\frac{d T}{d t} = \frac{1}{\tau} \left[T_0 - T(t)\right]
\end{array}
\tag{7-2} d t d T = τ 1 [ T 0 − T ( t ) ] ( 7-2 ) 其中τ \tau τ T ( t + Δ 1 2 ) T(t+\Delta \frac{1}{2}) T ( t + Δ 2 1 ) t − Δ t 2 t-\Delta \frac{t}{2} t − Δ 2 t
T ( t + Δ t 2 ) = T ( t − Δ t 2 ) + [ d T d t ] t − Δ t / 2 ⋅ Δ t + O [ ( Δ t ) 2 ] (7-3) \begin{array}{l}
T\left(t + \frac{\Delta t}{2}\right) = T\left(t - \frac{\Delta t}{2}\right) + \left[ \frac{dT}{dt} \right]_{t - \Delta t / 2} \cdot \Delta t + \mathcal{O}\left[(\Delta t)^2\right]
\end{array}
\tag{7-3} T ( t + 2 Δ t ) = T ( t − 2 Δ t ) + [ d t d T ] t − Δ t /2 ⋅ Δ t + O [ ( Δ t ) 2 ] ( 7-3 ) 再运用Fourier定律得到
T ( t + Δ t / 2 ) T ( t − Δ t / 2 ) ≈ 1 + 1 T ( t − Δ t / 2 ) ⋅ 1 τ [ T 0 − T ( t − Δ t / 2 ) ] ⋅ Δ t = 1 + [ T 0 T ( t − Δ t / 2 ) − 1 ] ⋅ Δ t τ ≡ λ 2 (7-4) \begin{array}{ll}
\frac{T(t + \Delta t/2)}{T(t - \Delta t/2)} &\approx 1 + \frac{1}{T(t - \Delta t/2)} \cdot \frac{1}{\tau} [T_0 - T(t - \Delta t/2)] \cdot \Delta t \\
&= 1 + \left[ \frac{T_0}{T(t - \Delta t/2)} - 1 \right] \cdot \frac{\Delta t}{\tau} \equiv \lambda^2 \\
\end{array}
\tag{7-4} T ( t − Δ t /2 ) T ( t + Δ t /2 ) ≈ 1 + T ( t − Δ t /2 ) 1 ⋅ τ 1 [ T 0 − T ( t − Δ t /2 )] ⋅ Δ t = 1 + [ T ( t − Δ t /2 ) T 0 − 1 ] ⋅ τ Δ t ≡ λ 2 ( 7-4 ) 标度因子λ \lambda λ
λ = 1 + [ T 0 T ( t − Δ t / 2 ) − 1 ] ⋅ Δ t τ = 1 + ( T 0 − T T ) ⋅ Δ t τ (7-5) \begin{array}{l}
\lambda = \sqrt{1 + \left[ \frac{T_0}{T(t - \Delta t/2)} - 1 \right] \cdot \frac{\Delta t}{\tau}} = \sqrt{1 + \left( \frac{T_0 - T}{T} \right) \cdot \frac{\Delta t}{\tau}}
\end{array}
\tag{7-5} λ = 1 + [ T ( t − Δ t /2 ) T 0 − 1 ] ⋅ τ Δ t = 1 + ( T T 0 − T ) ⋅ τ Δ t ( 7-5 ) 通过将标度因子代入粒子速度从而控制体系温度。
7.3 Bussi控温
Bussi控温在Berdensen的基础上引入了Weiner随机项，以保证能够给出正确的动能分布:
λ = 1 + ( T 0 − T T ) ⋅ Δ t τ + 2 T 0 T f τ d W (7-6) \begin{array}{l}
\lambda = \sqrt{1 + \left( \frac{T_0 - T}{T} \right) \cdot \frac{\Delta t}{\tau}} + 2 \sqrt{\frac{T_0}{T f \tau}} \mathrm{d} W
\end{array}
\tag{7-6} λ = 1 + ( T T 0 − T ) ⋅ τ Δ t + 2 T f τ T 0 d W ( 7-6 ) 其中f f f d W \mathrm{d}W d W λ \lambda λ 0.8 0.8 0.8 1.2 1.2 1.2
7.4 Langevin控温
(居中)郎之万热浴利用Langevin方程引入涨落与阻尼项，从而对系统进行控温。对于粒子i i i
m i d 2 r i d t 2 = F i − γ m i d r i d t + 2 γ m i k B T R i ( t ) (7-7) \begin{array}{l}
m_i \frac{d^2 \mathbf{r}_i}{dt^2} = \mathbf{F}_i - \gamma m_i \frac{d \mathbf{r}_i}{dt} + \sqrt{2 \gamma m_i k_B T} \, \mathbf{R}_i(t)
\end{array}
\tag{7-7} m i d t 2 d 2 r i = F i − γ m i d t d r i + 2 γ m i k B T R i ( t ) ( 7-7 ) 其中γ \gamma γ R i ( t ) \mathbf{R}_i(t) R i ( t )
7.5 Nose-Hoover Chain控温
8.控压
控压主要是对周期性边界条件进行更改。类似于控温部分介绍的式(7.1)的分解，CudaSPONGE中对于控压主要采用如下的分解方式（P ^ ( Δ t ) = e x p ( − i L ^ P Δ t ) \hat P(\Delta t) = {\rm exp}(-i\hat L_P\Delta t) P ^ ( Δ t ) = exp ( − i L ^ P Δ t )
e x p ( − i L ^ Δ t ) ≈ v ^ ( Δ t 2 ) P ^ ( Δ t ) T ^ ( Δ t ) r ^ ( Δ t ) v ^ ( Δ t 2 ) ≈ v ^ ( Δ t 2 ) P ^ ( Δ t ) r ^ ( Δ t 2 ) T ^ ( Δ t ) r ^ ( Δ t 2 ) v ^ ( Δ t 2 ) (8-1) \begin{array}{lll}
{\rm exp}(-i\hat L \Delta t) & \approx & \hat v(\frac {\Delta t}{2}) \hat P(\Delta t) \hat T(\Delta t) \hat r(\Delta t)\hat v(\frac {\Delta t}{2}) \\
& \approx & \hat v(\frac {\Delta t}{2}) \hat P(\Delta t) \hat r(\frac {\Delta t}{2}) \hat T(\Delta t)\hat r(\frac {\Delta t}{2})\hat v(\frac {\Delta t}{2}) \\
\end{array}
\tag{8-1} exp ( − i L ^ Δ t ) ≈ ≈ v ^ ( 2 Δ t ) P ^ ( Δ t ) T ^ ( Δ t ) r ^ ( Δ t ) v ^ ( 2 Δ t ) v ^ ( 2 Δ t ) P ^ ( Δ t ) r ^ ( 2 Δ t ) T ^ ( Δ t ) r ^ ( 2 Δ t ) v ^ ( 2 Δ t ) ( 8-1 ) 对于NPT系综，我们此时的正则坐标是分数坐标s i \mathbf{s}_i s i b \mathbf{b} b s i \mathbf{s}_i s i b \mathbf{b} b
L = ( T − U ) s y s t e m + ( T − U ) b o x = 1 2 ∑ i = 1 N m i s i ˙ T b T b s i ˙ − ∑ i < j N Φ ( r i j ) + 1 2 m g t r ( b ˙ T b ˙ ) − p e x t V = 1 2 ∑ i = 1 N m i s i ˙ T G s i ˙ − ∑ i < j N Φ ( r i j ) + 1 2 m g t r ( b ˙ T b ˙ ) − p e x t V (8-2) \begin{array}{l}
L = (T - U)_{system} + (T-U)_{box} \\
= \frac{1}{2} \sum_{i=1}^N m_i \dot{\bold{s_i}}^T \mathbf{b}^T \mathbf{b} \dot{\bold{s_i}} - \sum_{i<j}^N \Phi({r_{ij}} )
+ \frac{1}{2} m_g \mathrm{tr}\left( \dot{\mathbf{b}}^T \dot{\mathbf{b}} \right) - p_{ext} V \\
= \frac{1}{2} \sum_{i=1}^N m_i \dot{\bold{s_i}}^T \bold{G} \dot{\bold{s_i}} - \sum_{i<j}^N \Phi({r_{ij}})
+ \frac{1}{2} m_g \mathrm{tr}\left( \dot{\mathbf{b}}^T \dot{\mathbf{b}} \right) - p_{ext} V
\end{array}
\tag{8-2} L = ( T − U ) sy s t e m + ( T − U ) b o x = 2 1 ∑ i = 1 N m i s i ˙ T b T b s i ˙ − ∑ i < j N Φ ( r ij ) + 2 1 m g tr ( b ˙ T b ˙ ) − p e x t V = 2 1 ∑ i = 1 N m i s i ˙ T G s i ˙ − ∑ i < j N Φ ( r ij ) + 2 1 m g tr ( b ˙ T b ˙ ) − p e x t V ( 8-2 ) 其中度规矩阵G = b T b \bold{G} = \mathbf{b}^T \mathbf{b} G = b T b m g m_g m g p e x t p_{ext} p e x t
d d t ( ∂ L ∂ s ˙ i ) − ∂ L ∂ s i = 0 d d t ( ∂ L ∂ b ˙ ) − ∂ L ∂ b = 0 (8-3) \begin{array}{l}
\displaystyle \frac{d}{dt}\left(\frac{\partial L}{\partial \dot{\mathbf{s}}_i}\right)
- \frac{\partial L}{\partial \mathbf{s}_i}=0 \\
\displaystyle \frac{d}{dt}\left(\frac{\partial L}{\partial \dot{\mathbf{b}}}\right)
- \frac{\partial L}{\partial \mathbf{b}}=0
\end{array}
\tag{8-3} d t d ( ∂ s ˙ i ∂ L ) − ∂ s i ∂ L = 0 d t d ( ∂ b ˙ ∂ L ) − ∂ b ∂ L = 0 ( 8-3 ) 压强张量（又被称为应力张量、压力张量等）被定义为
Π = − ( P e x t 0 0 0 P e x t 0 0 0 P e x t + γ s / b z z ) + 1 V [ ∑ i m i ( v i x v i x v i x v i y v i x v i z v i y v i x v i y v i y v i y v i z v i z v i x v i z v i y v i z v i z ) − ( ∂ U ∂ b x x ∂ U ∂ b x y ∂ U ∂ b x z ∂ U ∂ b y x ∂ U ∂ b y y ∂ U ∂ b y z ∂ U ∂ b z x ∂ U ∂ b z y ∂ U ∂ b z z ) b T ] (8-4) \begin{aligned}
\mathbf{\Pi}=&-\left(
\begin{array}{lll}
P_{\rm ext} & 0 & 0 \\
0 & P_{\rm ext} & 0 \\
0 & 0 & P_{\rm ext} + \gamma_s/b_{zz} \\
\end{array}
\right) \\
&+ \frac{1}{V}
\left[
\sum_i m_i
\left(
\begin{array}{lll}
v_{ix} v_{ix} & v_{ix} v_{iy} & v_{ix} v_{iz} \\
v_{iy} v_{ix} & v_{iy} v_{iy} & v_{iy} v_{iz} \\
v_{iz} v_{ix} & v_{iz} v_{iy} & v_{iz} v_{iz} \\
\end{array}
\right) \right. \\
&\left.\qquad -
\left(
\begin{array}{lll}
\frac{\partial U}{\partial b_{xx}} & \frac{\partial U}{\partial b_{xy}} & \frac{\partial U}{\partial b_{xz}} \\
\frac{\partial U}{\partial b_{yx}} & \frac{\partial U}{\partial b_{yy}} & \frac{\partial U}{\partial b_{yz}} \\
\frac{\partial U}{\partial b_{zx}} & \frac{\partial U}{\partial b_{zy}} & \frac{\partial U}{\partial b_{zz}} \\
\end{array}
\right)
\mathbf{b}^{\rm T}
\right]
\end{aligned}
\tag{8-4} Π = − P ext 0 0 0 P ext 0 0 0 P ext + γ s / b z z + V 1 i ∑ m i v i x v i x v i y v i x v i z v i x v i x v i y v i y v i y v i z v i y v i x v i z v i y v i z v i z v i z − ∂ b xx ∂ U ∂ b y x ∂ U ∂ b z x ∂ U ∂ b x y ∂ U ∂ b y y ∂ U ∂ b z y ∂ U ∂ b x z ∂ U ∂ b y z ∂ U ∂ b z z ∂ U b T ( 8-4 ) 上式中，势能对盒子的偏导数中，势能的自变量是盒子b \mathbf{b} b s i \mathbf{s}_i s i P e x t P_{\rm ext} P ext γ s \gamma_s γ s
如果我们需要各向同性（Isotropic）或半各向同性（Semiisotropic）或半各向异性（Semianisotropic）变化，可以在各向异性（Anisotropic）的基础上强行施加约束，三种约束后的Π ′ \mathbf{\Pi}' Π ′ Π \mathbf{\Pi} Π
Π ′ = { 1 3 ( T r ( Π ) 0 0 0 T r ( Π ) 0 0 0 T r ( Π ) ) i s o t r o p i c ( 1 2 ( Π x x + Π y y ) 0 0 0 1 2 ( Π x x + Π y y ) 0 0 0 Π z z ) s e m i i s o t r o p i c ( Π x x 0 0 0 Π y y 0 0 0 Π z z ) s e m i a n i s o t r o p i c (8-5) \mathbf{\Pi}'=
\left\{
\begin{array}{lll}
\frac{1}{3}\left(
\begin{array}{lll}
{\rm Tr}(\mathbf{\Pi}) & 0 & 0 \\
0 & {\rm Tr}(\mathbf{\Pi}) & 0 \\
0 & 0 & {\rm Tr}(\mathbf{\Pi}) \\
\end{array}
\right) & {\rm isotropic} \\
\left(
\begin{array}{lll}
\frac{1}{2}(\mathbf{\Pi}_{xx} + \mathbf{\Pi}_{yy}) & 0 & 0 \\
0 & \frac{1}{2}(\mathbf{\Pi}_{xx} + \mathbf{\Pi}_{yy}) & 0 \\
0 & 0 & \mathbf{\Pi}_{zz} \\
\end{array} \right) & {\rm semiisotropic} \\
\left(
\begin{array}{lll}
\mathbf{\Pi}_{xx} & 0 & 0 \\
0 & \mathbf{\Pi}_{yy} & 0 \\
0 & 0 & \mathbf{\Pi}_{zz} \\
\end{array} \right) & {\rm semianisotropic} \\
\end{array}
\right.
\tag{8-5} Π ′ = ⎩ ⎨ ⎧ 3 1 Tr ( Π ) 0 0 0 Tr ( Π ) 0 0 0 Tr ( Π ) 2 1 ( Π xx + Π y y ) 0 0 0 2 1 ( Π xx + Π y y ) 0 0 0 Π z z Π xx 0 0 0 Π y y 0 0 0 Π z z isotropic semiisotropic semianisotropic ( 8-5 ) 在此基础之上，SPONGE还提供了强行约束Π \mathbf{\Pi} Π
若势能U U U d U = ∑ i j d r i j ∂ U ∂ r i j = − ∑ i j d r i j F i j T dU=\sum_{ij} d\mathbf{r}_{ij} \frac{\partial U}{\partial \mathbf{r}_{ij}}=-\sum_{ij} d\mathbf{r}_{ij} \mathbf{F}_{ij}^{\rm T} d U = ∑ ij d r ij ∂ r ij ∂ U = − ∑ ij d r ij F ij T
− ( ∂ U ∂ b x x ∂ U ∂ b x y ∂ U ∂ b x z ∂ U ∂ b y x ∂ U ∂ b y y ∂ U ∂ b y z ∂ U ∂ b z x ∂ U ∂ b z y ∂ U ∂ b z z ) b T = ∑ i j r i j T F i j (8-6) -
\left(
\begin{array}{lll}
\frac{\partial U}{\partial b_{xx}} & \frac{\partial U}{\partial b_{xy}} & \frac{\partial U}{\partial b_{xz}} \\
\frac{\partial U}{\partial b_{yx}} & \frac{\partial U}{\partial b_{yy}} & \frac{\partial U}{\partial b_{yz}} \\
\frac{\partial U}{\partial b_{zx}} & \frac{\partial U}{\partial b_{zy}} & \frac{\partial U}{\partial b_{zz}} \\
\end{array}
\right)
\mathbf{b}^{\rm T} = \sum_{ij}\mathbf{r}^{\rm T}_{ij}\mathbf{F}_{ij}
\tag{8-6} − ∂ b xx ∂ U ∂ b y x ∂ U ∂ b z x ∂ U ∂ b x y ∂ U ∂ b y y ∂ U ∂ b z y ∂ U ∂ b x z ∂ U ∂ b y z ∂ U ∂ b z z ∂ U b T = ij ∑ r ij T F ij ( 8-6 ) 同样地，如果在计算势能的某部分U 0 U_0 U 0 r i \mathbf{r}_i r i F i \mathbf{F}_i F i
− ( ∂ U ∂ b x x ∂ U ∂ b x y ∂ U ∂ b x z ∂ U ∂ b y x ∂ U ∂ b y y ∂ U ∂ b y z ∂ U ∂ b z x ∂ U ∂ b z y ∂ U ∂ b z z ) b T = ∑ i r i T F i (8-7) -
\left(
\begin{array}{lll}
\frac{\partial U}{\partial b_{xx}} & \frac{\partial U}{\partial b_{xy}} & \frac{\partial U}{\partial b_{xz}} \\
\frac{\partial U}{\partial b_{yx}} & \frac{\partial U}{\partial b_{yy}} & \frac{\partial U}{\partial b_{yz}} \\
\frac{\partial U}{\partial b_{zx}} & \frac{\partial U}{\partial b_{zy}} & \frac{\partial U}{\partial b_{zz}} \\
\end{array}
\right)
\mathbf{b}^{\rm T} = \sum_i \mathbf{r}^{\rm T}_i \mathbf{F}_i
\tag{8-7} − ∂ b xx ∂ U ∂ b y x ∂ U ∂ b z x ∂ U ∂ b x y ∂ U ∂ b y y ∂ U ∂ b z y ∂ U ∂ b x z ∂ U ∂ b y z ∂ U ∂ b z z ∂ U b T = i ∑ r i T F i ( 8-7 ) 8.1 Andersen控压
Andersen控压和Andersen控温同时在文献[8.1] 中提出，其核心思想是将体积作为一个拓展自由度进行微正则系综和/或正则系综模拟，值得注意的是这个体系中，正则坐标是原子的分数坐标f \mathbf{f} f r \mathbf{r} r [8.2] 。在SPONGE中，这类拓展系综的算法被统一视为Andersen控压的延伸，使用的积分方法与目前存在的各方法均不相同，但总体而言与Martyna、Tuckerman、Tobias和Klein等人（MTTK控压，文献[8.3] ）较为相近。MTTK提出的控压算法中，使用了Nose-Hoover算法控制盒子的速度，而在SPONGE的实现中则使用Langevin算法控制盒子的速度（文献[8.4] ），同时我们认为自由度的倒数约等于0，也即1 / N f ≈ 0 1/N_f \approx 0 1/ N f ≈ 0 i i i r i \mathbf{r}_i r i v i \mathbf{v}_i v i b \mathbf{b} b g \mathbf{g} g
{ d r i d t = r i g d v i d t = − v i g d b d t = b g d g d t = L a n g e v i n ( g ) (8-8) \left\{
\begin{array}{lll}
\frac{d\mathbf{r}_i}{dt} & = & \mathbf{r}_i \mathbf{g}\\
\frac{d\mathbf{v}_i}{dt} & = & -\mathbf{v}_i\mathbf{g} \\
\frac{d\mathbf{b}}{dt} & = & \mathbf{b}\mathbf{g} \\
\frac{d\mathbf{g}}{dt} & = & {\rm Langevin}(\mathbf{g})
\end{array}
\right.
\tag{8-8} ⎩ ⎨ ⎧ d t d r i d t d v i d t d b d t d g = = = = r i g − v i g bg Langevin ( g ) ( 8-8 ) 各向同性时，上式的Π \mathbf{\Pi} Π R \mathbf{R} R g \mathbf{g} g b \mathbf{b} b P P P R R R g g g l l l
{ r i ′ = e x p ( g Δ t ) r i v i ′ = ( 1 − g Δ t ) v i l ′ = e x p ( g Δ t ) l g ′ = e x p ( − γ l n Δ t ) ( g + P V Δ t / m g ) + ( 1 − e x p ( − 2 γ l n Δ t ) ) k B T / m g R (8-9) \left\{
\begin{array}{lll}
\mathbf{r}'_i & = & {\rm exp}(g\Delta t)\mathbf{r}_i \\
\mathbf{v}'_i & = & (1-g\Delta t)\mathbf{v}_i \\
l' & = & {\rm exp}(g\Delta t)l \\
g' & = & {\rm exp}( -\gamma_{\rm ln} \Delta t) (g + PV \Delta t/ m_g) + \sqrt{(1 - {\rm exp}(-2\gamma_{\rm ln} \Delta t)) k_{\rm B} T / m_g} R
\end{array}
\right.
\tag{8-9} ⎩ ⎨ ⎧ r i ′ v i ′ l ′ g ′ = = = = exp ( g Δ t ) r i ( 1 − g Δ t ) v i exp ( g Δ t ) l exp ( − γ ln Δ t ) ( g + P V Δ t / m g ) + ( 1 − exp ( − 2 γ ln Δ t )) k B T / m g R ( 8-9 ) 上式中，m g m_g m g γ l n \gamma_{\rm ln} γ ln R \mathbf{R} R
半各向同性时，上式的Π \mathbf{\Pi} Π R \mathbf{R} R g \mathbf{g} g b \mathbf{b} b ∀ α ∈ { x , y , z } \forall \alpha \in \{x, y, z\} ∀ α ∈ { x , y , z }
{ r i α ′ = e x p ( g α α Δ t ) r i α v i α ′ = ( 1 − g α α Δ t ) v i α b α α ′ = e x p ( g α α Δ t ) b α α g α α ′ = e x p ( − γ l n Δ t ) ( g α α + Π α α V Δ t / m g ) + ( 1 − e x p ( − 2 γ l n Δ t ) ) k B T / m g R (8-10) \left\{
\begin{array}{lll}
\mathbf{r}'_{i\alpha} & = & {\rm exp}(\mathbf{g}_{\alpha\alpha}\Delta t)\mathbf{r}_{i\alpha} \\
\mathbf{v}'_{i\alpha} & = & (1-\mathbf{g}_{\alpha\alpha}\Delta t)\mathbf{v}_{i\alpha} \\
\mathbf{b}_{\alpha\alpha}' & = & {\rm exp}(\mathbf{g}_{\alpha\alpha}\Delta t)\mathbf{b}_{\alpha\alpha} \\
\mathbf{g}_{\alpha\alpha}' & = & {\rm exp}( -\gamma_{\rm ln} \Delta t)(\mathbf{g}_{\alpha\alpha} + \Pi_{\alpha\alpha}V \Delta t/ m_g) + \sqrt{(1 - {\rm exp}(-2\gamma_{\rm ln} \Delta t)) k_{\rm B} T / m_g} R
\end{array}
\right.
\tag{8-10} ⎩ ⎨ ⎧ r i α ′ v i α ′ b α α ′ g α α ′ = = = = exp ( g α α Δ t ) r i α ( 1 − g α α Δ t ) v i α exp ( g α α Δ t ) b α α exp ( − γ ln Δ t ) ( g α α + Π α α V Δ t / m g ) + ( 1 − exp ( − 2 γ ln Δ t )) k B T / m g R ( 8-10 ) 各向异性时，将式(8.6)展开进行计算，此时对∀ α , β , γ ∈ { x , y , z } \forall \alpha,\beta,\gamma \in \{x, y, z\} ∀ α , β , γ ∈ { x , y , z } α ≠ β \alpha \neq \beta α = β
{ r i α ′ = e x p ( ∑ γ r i γ r i α g γ α Δ t ) r i α v i α ′ = v i α − ∑ γ v i γ g γ α Δ t b α α ′ = e x p ( ∑ γ b α γ b α α g γ α Δ t ) b α α b α β ′ = b α β + ∑ γ b α γ g γ β Δ t g α γ ′ = e x p ( − γ l n Δ t ) ( g α γ + Π α γ V Δ t / m g ) + ( 1 − e x p ( − 2 γ l n Δ t ) ) k B T / m g R (8-11) \left\{
\begin{array}{lll}
\mathbf{r}'_{i\alpha} & = & {\rm exp}(\sum_\gamma \frac{\mathbf{r}_{i\gamma}}{\mathbf{r}_{i\alpha}}\mathbf{g}_{\gamma\alpha}\Delta t)\mathbf{r}_{i\alpha} \\
\mathbf{v}'_{i\alpha} & = & \mathbf{v}_{i\alpha} - \sum_\gamma \mathbf{v}_{i\gamma}\mathbf{g}_{\gamma\alpha}\Delta t \\
\mathbf{b}_{\alpha\alpha}' & = & {\rm exp}(\sum_\gamma \frac{\mathbf{b}_{\alpha\gamma }}{\mathbf{b}_{\alpha\alpha}}\mathbf{g}_{\gamma\alpha}\Delta t) \mathbf{b}_{\alpha\alpha} \\
\mathbf{b}_{\alpha\beta}' & = & \mathbf{b}_{\alpha\beta} + \sum_\gamma \mathbf{b}_{\alpha\gamma}\mathbf{g}_{\gamma\beta} \Delta t \\
\mathbf{g}_{\alpha\gamma}' & = & {{\rm exp}( -\gamma_{\rm ln}\Delta t)} (\mathbf{g}_{\alpha\gamma} + \Pi_{\alpha\gamma}V \Delta t/ m_g) + \sqrt{(1 - {\rm exp}(-2\gamma_{\rm ln} \Delta t)) k_{\rm B} T / m_g} R
\end{array}
\right.
\tag{8-11} ⎩ ⎨ ⎧ r i α ′ v i α ′ b α α ′ b α β ′ g α γ ′ = = = = = exp ( ∑ γ r i α r iγ g γ α Δ t ) r i α v i α − ∑ γ v iγ g γ α Δ t exp ( ∑ γ b α α b α γ g γ α Δ t ) b α α b α β + ∑ γ b α γ g γ β Δ t exp ( − γ ln Δ t ) ( g α γ + Π α γ V Δ t / m g ) + ( 1 − exp ( − 2 γ ln Δ t )) k B T / m g R ( 8-11 ) 8.2 Berendsen控压
从Andersen控压出发，l i m γ l n → ∞ , T → 0 {\rm lim} \gamma_{\rm ln} \to \infty, T\to 0 lim γ ln → ∞ , T → 0
g α β = Π α β V Δ t / m g (8-12) \mathbf{g}_{\alpha\beta} = \Pi_{\alpha\beta}V \Delta t/ m_g
\tag{8-12} g α β = Π α β V Δ t / m g ( 8-12 ) 8.3 Bussi控压
从Andersen控压出发，l i m γ l n → ∞ {\rm lim} \gamma_{\rm ln} \to \infty lim γ ln → ∞
g α β = Π α β V Δ t / m g + k B T / m g R (8-13) \mathbf{g}_{\alpha\beta} = \Pi_{\alpha\beta}V \Delta t/ m_g + \sqrt{k_{\rm B} T / m_g} R
\tag{8-13} g α β = Π α β V Δ t / m g + k B T / m g R ( 8-13 ) 8.4 MC控压
9.虚原子
虚原子（virtual atom）也被称虚拟原子、为虚拟相互作用位点（virtual interaction site）、虚拟位点、哑原子（dump atom）等，将部分复杂的势能简化为另一个虚拟的原子的方法。例如对于势能U = r 0 − 6 + r 0 − 12 U=r_0^{-6} + r_0^{-12} U = r 0 − 6 + r 0 − 12 r v = r 0 2 r_v=r_0 ^ 2 r v = r 0 2 U = r 0 − 6 + r v − 6 U=r_0^{-6} + r_v^{-6} U = r 0 − 6 + r v − 6 U = r 0 − 6 + r 0 − 12 + r 0 − 24 U=r_0^{-6} + r_0^{-12} + r_0^{-24} U = r 0 − 6 + r 0 − 12 + r 0 − 24 r 1 = r 0 2 r_1=r_0 ^ 2 r 1 = r 0 2 r 2 = r 1 2 r_2=r_1 ^ 2 r 2 = r 1 2 U = r 0 − 6 + r 1 − 6 + r 2 − 6 U = r_0^{-6} + r_1^{-6} + r_2^{-6} U = r 0 − 6 + r 1 − 6 + r 2 − 6 n + 1 n+1 n + 1 r n \mathbf{r}_n r n n n n
对于包含虚原子的模拟，我们先将虚原子当作真实存在的原子进行普通计算，计算完以后需要重新分配力。对势能的原始形式和包含虚原子的形式分别进行微分，有：
d U = d b ( ∂ U ∂ b ) r 0 + d r 0 ( ∂ U ∂ r 0 ) b = d b ( ∂ U ∂ b ) r i + ∑ n d r n ( ∂ U ∂ r n ) r i , i ≠ n , b (9-1) \begin{array}{lllll}
dU & = & d{\mathbf{b}} \left(\frac{\partial U}{\partial \mathbf{b}}\right)_{\mathbf{r}_{0}} & + & d{\mathbf{r_0}} \left(\frac{\partial U}{\partial \mathbf{r_0}}\right)_{\mathbf{b}}
\\
& = & d{\mathbf{b}} \left(\frac{\partial U}{\partial \mathbf{b}}\right)_{\mathbf{r}_{i}} & + & \sum_n d{\mathbf{r}_n} \left(\frac{\partial U}{\partial \mathbf{r}_n}\right)_{\mathbf{r}_{i,i \neq n},\mathbf{b}}
\end{array}
\tag{9-1} d U = = d b ( ∂ b ∂ U ) r 0 d b ( ∂ b ∂ U ) r i + + d r 0 ( ∂ r 0 ∂ U ) b ∑ n d r n ( ∂ r n ∂ U ) r i , i = n , b ( 9-1 ) 上式中，偏导括号外的部分表示在求偏导时不变的量。
根据上式，对于力的计算，有
( ∂ U ∂ r 0 ) b = ∑ n d r n d r 0 ( ∂ U ∂ r n ) r i , i ≠ n , b (9-2) \left(\frac{\partial U}{\partial \mathbf{r_0}}\right)_{\mathbf{b}} =
\sum_n \frac {d{\mathbf{r}_n}}{d\mathbf{r}_0} \left(\frac{\partial U}{\partial \mathbf{r_n}}\right)_{\mathbf{r}_{i,i \neq n},\mathbf{b}}
\tag{9-2} ( ∂ r 0 ∂ U ) b = n ∑ d r 0 d r n ( ∂ r n ∂ U ) r i , i = n , b ( 9-2 ) d r n / d r 0 {d{\mathbf{r}_n}} / {d\mathbf{r}_0} d r n / d r 0
而对于维里计算，则还需要再将势能按两种情况的分数坐标f n \mathbf{f}_n f n
d U = d b ( ∂ U ∂ b ) f 0 + d f 0 ( ∂ U ∂ f 0 ) b = d b ( ∂ U ∂ b ) f i + ∑ n d f n ( ∂ U ∂ f n ) f i , i ≠ n , b (9-3) \begin{array}{lllll}
dU & = & d{\mathbf{b}} \left(\frac{\partial U}{\partial \mathbf{b}}\right)_{\mathbf{f}_{0}} & + & d{\mathbf{f_0}} \left(\frac{\partial U}{\partial \mathbf{f_0}}\right)_{\mathbf{b}}
\\
& = & d{\mathbf{b}} \left(\frac{\partial U}{\partial \mathbf{b}}\right)_{\mathbf{f}_{i}} & + & \sum_n d{\mathbf{f}_n} \left(\frac{\partial U}{\partial \mathbf{f}_n}\right)_{\mathbf{f}_{i,i \neq n},\mathbf{b}}
\end{array}
\tag{9-3} d U = = d b ( ∂ b ∂ U ) f 0 d b ( ∂ b ∂ U ) f i + + d f 0 ( ∂ f 0 ∂ U ) b ∑ n d f n ( ∂ f n ∂ U ) f i , i = n , b ( 9-3 ) 结合(9-1)和(9-3)
( ∂ U ∂ b ) f 0 = ( ∂ U ∂ b ) f i + ∑ n ( ∂ f n ∂ b ) f 0 ( ∂ U ∂ f n ) f i , i ≠ n , b = ( ∂ U ∂ b ) f i − ∑ n ( ∂ f n ∂ b ) f 0 b F n T (9-4) \begin{array}{llllll}
\left(\frac{\partial U}{\partial \mathbf{b}}\right)_{\mathbf{f}_{0}} & = & \left(\frac{\partial U}{\partial \mathbf{b}}\right)_{\mathbf{f}_{i}} & + & \sum_n & \left(\frac{\partial \mathbf{f}_n}{\partial \mathbf{b}}\right)_{\mathbf{f}_{0}} \left(\frac{\partial U}{\partial \mathbf{f}_n}\right)_{\mathbf{f}_{i,i \neq n},\mathbf{b}} \\
& = & \left(\frac{\partial U}{\partial \mathbf{b}}\right)_{\mathbf{f}_{i}} & - & \sum_n & \left(\frac{\partial \mathbf{f}_n}{\partial \mathbf{b}}\right)_{\mathbf{f}_{0}} \mathbf{b} \mathbf{F}_n^{\rm T}
\end{array}
\tag{9-4} ( ∂ b ∂ U ) f 0 = = ( ∂ b ∂ U ) f i ( ∂ b ∂ U ) f i + − ∑ n ∑ n ( ∂ b ∂ f n ) f 0 ( ∂ f n ∂ U ) f i , i = n , b ( ∂ b ∂ f n ) f 0 b F n T ( 9-4 ) 从上面的讨论中可以看到，计算最关键的就在于获得r n \mathbf{r}_n r n
SPONGE中，目前支持5种形式的虚原子，分别定义如下：（以r v \mathbf{r}_v r v
种类0：依赖一个原子r 0 \mathbf{r}_0 r 0
r v = r 0 ( 1 0 0 0 1 0 0 0 − 1 ) f v = f 0 ( 1 0 0 0 1 0 0 0 − 1 ) \begin{array}{lll}
\mathbf{r}_v = \mathbf{r}_0 \left(\begin{array}{lll} 1 & 0 & 0 \\ 0 & 1 & 0 \\ 0 & 0 & -1\end{array}\right) \\
\mathbf{f}_v = \mathbf{f}_0 \left(\begin{array}{lll} 1 & 0 & 0 \\ 0 & 1 & 0 \\ 0 & 0 & -1\end{array}\right)
\end{array} r v = r 0 1 0 0 0 1 0 0 0 − 1 f v = f 0 1 0 0 0 1 0 0 0 − 1
种类1：依赖2个原子r 1 \mathbf{r}_1 r 1 r 2 \mathbf{r}_2 r 2 a a a
r v = a r 21 f v = a f 21 \begin{array}{l}
\mathbf{r}_v = a \mathbf{r}_{21} \\
\mathbf{f}_v = a \mathbf{f}_{21}
\end{array} r v = a r 21 f v = a f 21
种类2：依赖3个原子r 1 \mathbf{r}_1 r 1 r 2 \mathbf{r}_2 r 2 r 3 \mathbf{r}_3 r 3 a a a b b b
r v = a r 21 + b r 31 f v = a f 21 + b f 31 \begin{array}{l}
\mathbf{r}_v = a \mathbf{r}_{21} + b \mathbf{r}_{31} \\
\mathbf{f}_v = a \mathbf{f}_{21} + b \mathbf{f}_{31}
\end{array} r v = a r 21 + b r 31 f v = a f 21 + b f 31
种类3：依赖3个原子r 1 \mathbf{r}_1 r 1 r 2 \mathbf{r}_2 r 2 r 3 \mathbf{r}_3 r 3 d d d k k k
r v = d ∗ r 12 + k r 23 ∣ r 12 + k r 23 ∣ \begin{array}{l}
\mathbf{r}_v = d * \frac{ \mathbf{r}_{12} + k \mathbf{r}_{23}} {|\mathbf{r}_{12} + k \mathbf{r}_{23}|}
\end{array} r v = d ∗ ∣ r 12 + k r 23 ∣ r 12 + k r 23
种类4：依赖于n n n { r i } \{\mathbf{r}_i\} { r i } n n n
r v = ∑ i w i r i f v = ∑ i w i f i \begin{array}{l}
\mathbf{r}_v = \sum_i w_i \mathbf{r}_i \\
\mathbf{f}_v = \sum_i w_i \mathbf{f}_i
\end{array} r v = ∑ i w i r i f v = ∑ i w i f i
参考文献
对于晶胞，可以使用三个晶胞的边的长度、、和夹角、、进行描述，也可以将三个边的矢量构造成一个盒矩阵进行描述。通过合理建系，我们可以让第一条边在轴上，第二条边在平面上，第三条边则为任意位置。此时盒矩阵是一个三角矩阵。这里我们以矩阵的每一行作为边的行矢量，可以表示为一个下三角矩阵：