势能形式：
：之间的距离
、：力场参数
势能形式：
：之间的距离
、、、：力场参数
势能形式：
：角
、：力场参数
势能形式：
：之间的距离
：角abc
、、、：力场参数
势能形式：
：二面角
、、：力场参数
势能形式：
：二面角
、：力场参数
势能形式：
、：未修正的LJ势能和库伦势能
、： 力场参数
势能形式：
：之间的距离
、：的电荷、的电荷
、、：力场参数
势能形式：
：二面角
：二面角
：力场参数
势能形式：
：之间的距离
、：力场参数
势能形式：
: 之间的距离
、: 的电荷、的电荷
其中，表示虚原子，其余数字表示该虚原子依赖的其他原子的坐标
、、、、为力场参数，表示位移，、、为三个方向的坐标
: 与其参考点的距离
、参考点的位置: 力场参数
程序中仍保留，但尽量不要使用，因为居中朗之万效果始终比朗之万更好
实现伪代码：
加速度 = 力 / 质量
加速度 -= 摩擦因子 * 速度
加速度 += 根号(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
实现伪代码：
因子 = 1 + dt / 时间参数 * (系综温度 / 温度 - 1)
if (随机修正)
{
因子 += 2 * 根号(系综温度 / 温度 / 自由度 / 时间参数) * 标准高斯分布随机数);
}
因子 = 根号(因子)
if (因子 > 1.25)
{
因子 = 1.25
}
else if (因子 < 0.8)
{
因子 = 0.8
}
加速度 = 力 / 质量
速度 += 加速度 * dt
if (速度 > 最大速度)
{
速度 = 最大速度
}
坐标 += 速度 * dt
速度 *= 因子
实现伪代码：
if ((步数 - 1) % 更新间隔 == 0)
{
加速度 = 力 / 质量
速度 += 加速度 * dt
坐标 += 速度 * dt / 2
速度 = 玻尔兹曼速度分布的随机速度
if (速度 > 最大速度)
{
速度 = 最大速度
}
坐标 += 速度 * dt / 2
}
else
{
加速度 = 力 / 质量
速度 += 加速度 * dt
if (速度 > 最大速度)
{
速度 = 最大速度
}
坐标 += 速度 * dt
}
实现伪代码：
拓展维度质量 = 时间常数 * 时间常数 / 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迭代
约束迭代
伪代码实现：
NVT迭代
约束迭代
if (步数 % 变化间隔 == 0)
{
变化因子 = 1 - 变化间隔 * 压缩系数 * dt / 时间常数 / 3 * (系综压强 - 压强)
if (启用随机修正)
{
变化因子 += 根号(2 * 玻尔兹曼因子 * 系综温度 * 压缩系数 / 时间常数 / 新体积) / 3 * 标准正态分布随机数
速度 /= 变化因子
}
坐标 *= 变化因子
}
伪代码实现：
NVT迭代
约束迭代
if (步数 % 变化间隔 == 0)
{
体积速度 -= 变化间隔 * 压缩系数 * dt / 时间常数 / 3 * (系综压强 - 压强) / 体积质量
体积速度根据温度约化
变化因子 += 体积速度 * dt
坐标 *= 变化因子
}