力场

短程相互作用

1-2作用

简谐键

势能形式：

$$E = \sum k (r_{ab} - r_0) ^ 2$$

$r_{ab}$：$ab$之间的距离 $k$、$r_0$：力场参数

软简谐键

势能形式：

$$E = \sum \lambda k (r_{ab} - r_0) ^ 2 / (1 + \alpha \lambda (r_{ab} - r_0) ^ 2)$$

$r_{ab}$：$ab$之间的距离 $\lambda$、$\alpha$、$k$、$r_0$：力场参数

1-3作用

简谐键角

势能形式：

$$E = \sum k (\theta_{abc} - \theta_0) ^ 2$$

$\theta_{abc}$：角$abc$ $k$、$\theta_0$：力场参数

Urey_Bradley键角

势能形式：

$$E = \sum k_1 (r_{ac} - r_0) ^ 2 + k_2 (\theta_{abc} - \theta_0) ^ 2$$

$r_{ac}$：$ac$之间的距离 $\theta_{abc}$：角abc $k_1$、$k_2$、$r_0$、$\theta_0$：力场参数

1-4作用

周期性二面角

势能形式：

$$E = \sum k (1 + cos(n \phi_{abcd} - \phi_0))$$

$\phi_{abcd}$：二面角$abcd$ $n$、$k$、$\phi_0$：力场参数

简谐二面角

势能形式：

$$E = \sum k (\phi_{abcd} - \phi_0) ^ 2$$

$\phi_{abcd}$：二面角$abcd$ $k$、$\phi_0$：力场参数

双参数1-4非键作用

势能形式：

$$E = \sum k_{LJ} * E_{LJ} + k_{\text{库伦}} * E_{\text{库伦}}$$

$E_{LJ}$、$E_{\text{库伦}}$：未修正的LJ势能和库伦势能 $k_{LJ}$、$k_{\text{库伦}}$： 力场参数

三参数1-4非键作用

势能形式：

$$E = \sum A r_{ab}^{-12} - B r_{ab}^{-6} + C q_a q_b r_{ab}^{-1}$$

$r_{ab}$：$ab$之间的距离

$q_a$、$q_b$：$a$的电荷、$b$的电荷 $A$、$B$、$C$：力场参数

1-5作用

CMAP

势能形式：

$$E = \sum \sum_{i=0}^4 \sum_{j=0}^4 k_{ij} \phi_{abcd}^i \phi_{bcde}^j$$

$\phi_{abcd}$：二面角$abcd$ $\phi_{bcde}$：二面角$bcde$ $k_{ij}$：力场参数

粒子-粒子相互作用

范德华相互作用

LJ势

势能形式：

$$E = \sum A r_{ab}^{-12} - B r_{ab}^{-6}$$

$r_{ab}$：$ab$之间的距离 $A$、$B$：力场参数

静电相互作用

库伦势

势能形式：

$$E = \sum q_a q_b r_{ab}^{-1}$$

$r_{ab}$: $ab$之间的距离 $q_a$、$q_b$: $a$的电荷、$b$的电荷

其他

虚原子

$$\begin{array}{ll} \text{种类0:} & z_0 = 2 h - z_1 \\ \text{种类1:} & \vec r_{01} = a * \vec r_{21} \\ \text{种类2:} & \vec r_{01} = a * \vec r_{21} + b \vec r_{31} \\ \text{种类3:} & \vec r_{01} = d * (\vec r_{21} + k \vec r_{31}) / |\vec r_{21} + k \vec r_{31}| \end{array}$$

其中，$0$表示虚原子，其余数字表示该虚原子依赖的其他原子的坐标 $a$、$b$、$d$、$k$、$h$为力场参数，$\vec r$表示位移，$x$、$y$、$z$为三个方向的坐标

位置限制

$$E = k \sum r_{i,ref}^2$$

$r_{i,ref}$: $i$与其参考点的距离 $k$、参考点的位置: 力场参数

控温

朗之万动力学

程序中仍保留，但尽量不要使用，因为居中朗之万效果始终比朗之万更好

实现伪代码：

    加速度 = 力 / 质量
    加速度 -= 摩擦因子 * 速度
    加速度 += 根号(2 * 摩擦因子 * 系综温度 * 玻尔兹曼因子 / dt / 质量) * 标准高斯分布随机数
    速度 += 加速度 * dt
    if (速度 > 最大速度)
    {
        速度 = 最大速度
    }
    坐标 += 速度 * dt

居中朗之万控温

Zhang, Z. et. al., Zhang Z, Liu X, Chen Z, Zheng H, Yan K, Liu J. A unified thermostat scheme for efficient configurational sampling for classical/quantum canonical ensembles via molecular dynamics. J. Chem. Phys. 2017 Jul 21;147(3):034109. doi: 10.1063/1.4991621.

实现伪代码：

    加速度 = 力 / 质量
    速度 += 加速度 * dt
    if (速度 > 最大速度)
    {
        速度 = 最大速度
    }
    坐标 += 速度 * dt / 2
    衰减系数 = exp(-摩擦因子 * dt)
    速度 *= 衰减系数
    速度 += 根号((1 - 衰减系数 * 衰减系数) * 系综温度 * 玻尔兹曼因子 / 质量) * 标准高斯分布随机数
    坐标 += 速度 * dt / 2

Berendsen控温

实现伪代码：

 因子 = 1 + dt / 时间参数 * (系综温度 / 温度 - 1)
    if (随机修正)
    {
        因子 += 2 * 根号(系综温度 / 温度 / 自由度 / 时间参数) * 标准高斯分布随机数);
    }
    因子 = 根号(因子)
    if (因子 > 1.25)
    {
        因子 = 1.25
    }
    else if (因子 < 0.8)
    {
        因子 = 0.8
    }
    加速度 = 力 / 质量
    速度 += 加速度 * dt
    if (速度 > 最大速度)
    {
        速度 = 最大速度
    }
    坐标 += 速度 * dt
    速度 *= 因子

Andersen控温

实现伪代码：

if ((步数 - 1) % 更新间隔 == 0)
    {
        加速度 = 力 / 质量
        速度 += 加速度 * dt
        坐标 += 速度 * dt / 2
        速度 = 玻尔兹曼速度分布的随机速度
        if (速度 > 最大速度)
        {
            速度 = 最大速度
        }
        坐标 += 速度 * dt / 2
    }
    else
    {
        加速度 = 力 / 质量
        速度 += 加速度 * dt
        if (速度 > 最大速度)
        {
            速度 = 最大速度
        }
        坐标 += 速度 * dt
    }

Nose-Hoover控温

实现伪代码：

    拓展维度质量 = 时间常数 * 时间常数 / 4 / π / π
    拓展维度速度[0] = (2 * 总动能 - 自由度 * 玻尔兹曼因子 * 系综温度) / 拓展维度质量
    if (链长度 > 1)
    {
        拓展维度速度[0] -= 拓展维度速度[0] * 拓展维度速度[1] * dt
    }
    for (int i = 1; i < 链长度 - 1; i += 1)
    {
        拓展维度速度[i] = 拓展维度速度[i - 1] * 拓展维度速度[i - 1] - 玻尔兹曼因子 * 系综温度) / 拓展维度质量
        拓展维度速度[i] *= dt
        拓展维度速度[i] -= 拓展维度速度[0] * 拓展维度速度[1] * dt
    }
    拓展维度速度[链长度 - 1] = 拓展维度速度[链长度 - 2] * 拓展维度速度[链长度 - 2] - 玻尔兹曼因子 * 系综温度 / 拓展维度质量
    拓展维度速度[链长度 - 1] *= dt
    加速度 = 力 / 质量 - 速度 * 拓展维度速度[0]
    速度 += 加速度 * dt
    if (速度 > 最大速度)
    {
        速度 = 最大速度
    }
    坐标 += 速度 * dt

控压

蒙特卡洛控压

伪代码实现：

    if (步数 % 变化间隔 == 0)
    {
        记录旧的能量和力，随机改变体积，计算新的能量和力
        额外项 = 系综压强* (新体积 - 旧体积)
        额外项 -= 粒子数 * 玻尔兹曼因子 * 系综温度 * ln(新体积 / 旧体积)
        额外项 -= 表面数量 * 表面张力 * (新表面积 - 旧表面积)
        接受概率 = exp(- (新能量 - 旧能量 + 额外项) / 玻尔兹曼因子 / 系综温度)
        if (不接受)
        {
            恢复旧体积、能量和力
        }
        if (尝试次数 % 参数检查间隔 == 0)
        {
            if (接受率 < 最小接受率)
            {
                体积变化最大值 *= 0.9
            }
            if (接受率 > 最大接受率)
            {
                体积变化最大值 *= 1.1
            }
        }
    }
    NVT迭代
    约束迭代

Berendsen控压

伪代码实现：

    NVT迭代
    约束迭代
    if (步数 % 变化间隔 == 0)
    {
        变化因子 = 1 - 变化间隔 * 压缩系数 * dt / 时间常数 / 3 * (系综压强 - 压强)
        if (启用随机修正)
        {
            变化因子 += 根号(2 * 玻尔兹曼因子 * 系综温度 * 压缩系数 / 时间常数 / 新体积) / 3 * 标准正态分布随机数
            速度 /= 变化因子
        }
        坐标 *= 变化因子
    }

Andersen控压

伪代码实现：

    NVT迭代
    约束迭代
    if (步数 % 变化间隔 == 0)
    {
        体积速度 -= 变化间隔 * 压缩系数 * dt / 时间常数 / 3 * (系综压强 - 压强) / 体积质量
        体积速度根据温度约化
        变化因子 += 体积速度 * dt
        坐标 *= 变化因子
    }