MindSponge分子动力学模拟——Constraint约束
技术背景
在前面的几篇博客中，我们已经介绍了MindSponge的基本使用方法，比如定义一个分子系统、计算分子的单点能以及迭代器的使用等。有了这些基础的教程，用户以及可以执行一些比较简单的模拟任务，比如可以跑一个能量极小化，或者是NVT过程。
当我们去执行一个模拟任务时，比较关键的一个指标是给定时间内可模拟的总时长，总时长直接决定了我们能不能看到其中的一些过渡态。当然，现在有很多改进的思路。比如可以使用硬件加速，缩减单步模拟的时长。比如可以使用增强采样，朝着给定的方向演化，更容易看到中间态。本文要介绍的方法也是其中之一，可以采用Constraint约束的方法，限制分子体系中某些不重要的自由度，这样演化步长\(dt\)就可以设置的大一些，也能够起到加速的作用。常用的Constraint约束算法有Lincs和SETTLE算法。而Shake算法由于需要迭代，因此未在MindSponge中实现。
这里还需要谈一下Constraint和Restraint的区别。一般情况下，Constraint指代一些比较强硬的约束，比如SETTLE算法，直接就固定了一个三角形的形状大小，在所有的演化过程中都不会发生变化。而Restraint更多的指代一些比较灵活的软约束，比如加一个球形谐振子势把分子系统限制在一个球体里面，就是一种Restraint约束。而Restraint也并不全是软约束，比如在MetaDynamics中会在边界处施加一个梯度非常大的Restraint约束，虽然理论上是一个软约束，但是CV一般无法越过这个边界。
不施加约束
本案例类似于上一篇博客中所介绍的NVT过程，不过上一篇博客里面用的体系是一个纯粹的多肽体系，这里我们用Molecule自带的加水的方法，模拟一个带溶剂的多肽混合体系。使用的多肽是一个简单的4个residue的体系：
REMARK Generated By Xponge (Molecule)
ATOM 1 N ALA 1 -0.095 -11.436 -0.780
ATOM 2 CA ALA 1 -0.171 -10.015 -0.507
ATOM 3 CB ALA 1 1.201 -9.359 -0.628
ATOM 4 C ALA 1 -1.107 -9.319 -1.485
ATOM 5 O ALA 1 -1.682 -9.960 -2.362
ATOM 6 N ARG 2 -1.303 -8.037 -1.397
ATOM 7 CA ARG 2 -2.194 -7.375 -2.328
ATOM 8 CB ARG 2 -3.606 -7.943 -2.235
ATOM 9 CG ARG 2 -4.510 -7.221 -3.228
ATOM 10 CD ARG 2 -5.923 -7.789 -3.136
ATOM 11 NE ARG 2 -6.831 -7.111 -4.087
ATOM 12 CZ ARG 2 -8.119 -7.421 -4.205
ATOM 13 NH1 ARG 2 -8.686 -8.371 -3.468
ATOM 14 NH2 ARG 2 -8.844 -6.747 -5.093
ATOM 15 C ARG 2 -2.273 -5.882 -2.042
ATOM 16 O ARG 2 -1.630 -5.388 -1.119
ATOM 17 N ALA 3 -3.027 -5.119 -2.777
ATOM 18 CA ALA 3 -3.103 -3.697 -2.505
ATOM 19 CB ALA 3 -1.731 -3.041 -2.625
ATOM 20 C ALA 3 -4.039 -3.001 -3.483
ATOM 21 O ALA 3 -4.614 -3.643 -4.359
ATOM 22 N ALA 4 -4.235 -1.719 -3.394
ATOM 23 CA ALA 4 -5.126 -1.057 -4.325
ATOM 24 CB ALA 4 -6.538 -1.625 -4.233
ATOM 25 C ALA 4 -5.205 0.436 -4.039
ATOM 26 O ALA 4 -4.561 0.930 -3.116
ATOM 27 OXT ALA 4 -5.915 1.166 -4.728
TER具体的使用方法这里就不多做介绍了，有相关的疑问的，可以在这个博客列表中去寻找相关问题的答案。一个不包含约束的MindSponge分子动力学模拟代码如下：
from mindspore import nn, context
from sponge import ForceField, Sponge, set_global_units, Protein, UpdaterMD
from sponge.function import VelocityGenerator
from sponge.callback import RunInfo, WriteH5MD
context.set_context(mode=context.GRAPH_MODE, device_target='GPU', device_id=1)
set_global_units('A', 'kcal/mol')
# 加载一个纯粹的多肽体系
system = Protein('../pdb/case1.pdb', template=['protein0.yaml'], rebuild_hydrogen=True)
# 在多肽体系的周边填充厚度为4A的水分子，同时系统变为周期性重复的体系
system.fill_water(edge=4, template='water.tip3p.yaml')
# 加载两个不同的力场，如果指定不使用PME算法，模拟速度会快一点，但存在能量偏移的问题，
# 而这里只做展示用，所以暂时关闭PME算法的使用
energy = ForceField(system, parameters=['AMBER.FF99SB', 'TIP3P'], use_pme=False)
# 能量极小化过程，用的是MindSpore自带的Adam优化器
min_opt = nn.Adam(system.trainable_params(), 1e-03)
md = Sponge(system, potential=energy, optimizer=min_opt)
run_info = RunInfo(200)
md.run(2000, callbacks=[run_info])
# 定义一个新的迭代器，并用change_optimizer完成切换
vgen = VelocityGenerator(300)
velocity = vgen(system.shape, system.atom_mass)
# 定义了NVT迭代器，指定300K的控温
opt = UpdaterMD(system=system,
time_step=1e-3,
velocity=velocity,
integrator='velocity_verlet',
temperature=300,
thermostat='langevin')
md.change_optimizer(opt)
nvt_h5md = WriteH5MD(system, 'test_nvt.h5md', save_freq=10, write_image=False)
md.run(2000, callbacks=[run_info, nvt_h5md])该代码的运行结果如下：
[MindSPONGE] Adding 45 hydrogen atoms for the protein molecule in 0.003 seconds.
[MindSPONGE] The box size used when filling water is: [21.184929 22.856928 15.495929]
[MindSPONGE] The edge gap along x axis is 4.0.
[MindSPONGE] The edge gap along y axis is 3.999999523162842.
[MindSPONGE] The edge gap along z axis is 3.999999761581421.
[MindSPONGE] Totally 172 waters is added to the system!
[MindSPONGE] Adding water molecule task finished!
[MindSPONGE] Started simulation at 2023-09-05 15:31:31
[MindSPONGE] Step: 0, E_pot: 6019.2305
[MindSPONGE] Step: 200, E_pot: 3667.4219
[MindSPONGE] Step: 400, E_pot: 3007.9343
[MindSPONGE] Step: 600, E_pot: 2633.8374
[MindSPONGE] Step: 800, E_pot: 2380.708
[MindSPONGE] Step: 1000, E_pot: 2233.6882
[MindSPONGE] Step: 1200, E_pot: 2129.8638
[MindSPONGE] Step: 1400, E_pot: 2050.8691
[MindSPONGE] Step: 1600, E_pot: 1993.4315
[MindSPONGE] Step: 1800, E_pot: 1949.2118
[MindSPONGE] Finished simulation at 2023-09-05 15:32:19
[MindSPONGE] Simulation time: 47.67 seconds.
--------------------------------------------------------------------------------
[MindSPONGE] Started simulation at 2023-09-05 15:32:20
[MindSPONGE] Step: 0, E_pot: 1910.4984, E_kin: 540.1826, E_tot: 2450.6812, Temperature: 316.81873, Pressure: -2254.547, Volume: 7503.4756
[MindSPONGE] Step: 200, E_pot: 2165.2405, E_kin: 421.29523, E_tot: 2586.5356, Temperature: 247.09091, Pressure: -4248.9507, Volume: 7503.4756
[MindSPONGE] Step: 400, E_pot: 2167.1772, E_kin: 479.70053, E_tot: 2646.8777, Temperature: 281.34576, Pressure: -5078.793, Volume: 7503.4756
[MindSPONGE] Step: 600, E_pot: 2196.748, E_kin: 480.39355, E_tot: 2677.1416, Temperature: 281.75226, Pressure: -2657.5989, Volume: 7503.4756
[MindSPONGE] Step: 800, E_pot: 2191.3103, E_kin: 498.1831, E_tot: 2689.4934, Temperature: 292.18588, Pressure: 506.58923, Volume: 7503.4756
[MindSPONGE] Step: 1000, E_pot: 2203.0522, E_kin: 490.20367, E_tot: 2693.2559, Temperature: 287.5059, Pressure: -3832.6213, Volume: 7503.4756
[MindSPONGE] Step: 1200, E_pot: 2205.9346, E_kin: 475.74023, E_tot: 2681.6748, Temperature: 279.02304, Pressure: -3744.449, Volume: 7503.4756
[MindSPONGE] Step: 1400, E_pot: 2191.0286, E_kin: 502.4231, E_tot: 2693.4517, Temperature: 294.67264, Pressure: -3162.2998, Volume: 7503.4756
[MindSPONGE] Step: 1600, E_pot: 2200.712, E_kin: 511.2496, E_tot: 2711.9614, Temperature: 299.84943, Pressure: -2348.196, Volume: 7503.4756
[MindSPONGE] Step: 1800, E_pot: 2209.8, E_kin: 514.364, E_tot: 2724.164, Temperature: 301.67603, Pressure: -2295.4656, Volume: 7503.4756
[MindSPONGE] Finished simulation at 2023-09-05 15:33:25
[MindSPONGE] Simulation time: 1 minutes 5.4 seconds.
--------------------------------------------------------------------------------需要注意的是，虽然我们在补水的时候加载的模板（如下所示）里面有SETTLE约束算法相关的配置信息，但是因为在UpdaterMD中未指定constraint参数，因此默认也是不添加约束条件的。
template:
base: water_3p.yaml
WAT:
atom_charge: [-0.834, 0.417, 0.417]
settle:
mandatory: false
length_unit: nm
distance:
OW-HW: 0.09572
HW-HW: 0.15139
molecule:
residue:
- WAT
length_unit: nm
coordinate:
- [0.0, 0.0, 0.0]
- [0.079079641, 0.061207927, 0.0]
- [-0.079079641, 0.061207927, 0.0]可以通过VMD插件来查看运行的轨迹，我们可以从中发现，因为未添加任何的约束条件，因此轨迹中的各个键长键角都在不断的变化。
水分子SETTLE约束
在上一个章节中，我们所运行的体系是一个多肽和水溶液的混合体系，并且没有施加任何的约束条件。这里我们先介绍第一种约束条件，对水分子施加一个SETTLE约束，简单来说，就是固定住水分子的三个原子所形成的三角形。其实这个约束更多的是在TIP4P力场中发挥作用，不过在2023年9月份的当前MindSponge版本，暂不支持TIP4P力场。
这里先简要介绍一下SETTLE约束算法的原理（如上图所示），方便大家了解SETTLE约束算法的适用范围。三角形\(A_0B_0C_0\)是原始的三角形，三角形\(A_1B_1C_1\)是在力场作用下移动了一步的三角形。我们的目标，是对新的三角形的三个顶点进行平移变换，得到一个符合约束条件（SETTLE算法的约束条件就是三个边长参数，OW-HW和HW-HW要等于给定值）的结果。而新三角形的三个点\(A_1,B_1,C_1\)，只能够在原始三角形平面的切向移动，移动后找到一个符合边长约束条件的三角形\(A_3B_3C_3\)作为最终的结果。
由此可见，SETTLE约束算法的适用范围，只能作用于三角形的体系，且必须是等腰三角形。在MindSponge中适用SETTLE约束算法的方法非常简单，就是在Updater迭代器的参数中添加一个constraint的参数即可，把构造好的SETTLE对象传进去。当然，也可以直接从模板文件里面加载SETTLE对象，这个后面章节会提到。调用SETTLE算法相关的代码实现如下：
from mindspore import nn, context
from sponge import ForceField, Sponge, set_global_units, Protein, UpdaterMD
from sponge.function import VelocityGenerator
from sponge.control import SETTLE
from sponge.callback import RunInfo, WriteH5MD
context.set_context(mode=context.GRAPH_MODE, device_target='GPU', device_id=1)
set_global_units('A', 'kcal/mol')
system = Protein('../pdb/case1.pdb', template=['protein0.yaml'], rebuild_hydrogen=True)
system.fill_water(edge=4, template='water.tip3p.yaml')
energy = ForceField(system, parameters=['AMBER.FF99SB', 'TIP3P'], use_pme=False)
min_opt = nn.Adam(system.trainable_params(), 1e-03)
md = Sponge(system, potential=energy, optimizer=min_opt)
run_info = RunInfo(200)
md.run(2000, callbacks=[run_info])
vgen = VelocityGenerator(300)
velocity = vgen(system.shape, system.atom_mass)
# 在迭代器中传入SETTLE对象
opt = UpdaterMD(system=system,
time_step=1e-3,
velocity=velocity,
integrator='velocity_verlet',
temperature=300,
thermostat='langevin',
constraint=SETTLE(system))
md.change_optimizer(opt)
nvt_h5md = WriteH5MD(system, 'test_nvt.h5md', save_freq=10, write_image=False)
md.run(2000, callbacks=[run_info, nvt_h5md])运行结果如下所示：
[MindSPONGE] Adding 45 hydrogen atoms for the protein molecule in 0.003 seconds.
[MindSPONGE] The box size used when filling water is: [21.184929 22.856928 15.495929]
[MindSPONGE] The edge gap along x axis is 4.0.
[MindSPONGE] The edge gap along y axis is 3.999999523162842.
[MindSPONGE] The edge gap along z axis is 3.999999761581421.
[MindSPONGE] Totally 172 waters is added to the system!
[MindSPONGE] Adding water molecule task finished!
[MindSPONGE] Started simulation at 2023-09-05 16:05:08
[MindSPONGE] Step: 0, E_pot: 6019.2305
[MindSPONGE] Step: 200, E_pot: 3667.421
[MindSPONGE] Step: 400, E_pot: 3007.9253
[MindSPONGE] Step: 600, E_pot: 2633.8325
[MindSPONGE] Step: 800, E_pot: 2380.7039
[MindSPONGE] Step: 1000, E_pot: 2233.6877
[MindSPONGE] Step: 1200, E_pot: 2129.8643
[MindSPONGE] Step: 1400, E_pot: 2050.8638
[MindSPONGE] Step: 1600, E_pot: 1993.42
[MindSPONGE] Step: 1800, E_pot: 1949.1956
[MindSPONGE] Finished simulation at 2023-09-05 16:05:53
[MindSPONGE] Simulation time: 45.13 seconds.
--------------------------------------------------------------------------------
[MindSPONGE] The settle constraint is used for the molecule system.
[MindSPONGE] Started simulation at 2023-09-05 16:05:54
[MindSPONGE] Step: 0, E_pot: 1910.4838, E_kin: 517.4474, E_tot: 2427.9312, Temperature: 303.48444, Pressure: -813660.7, Volume: 7503.4756
[MindSPONGE] Step: 200, E_pot: 2011.8678, E_kin: 235.43576, E_tot: 2247.3035, Temperature: 197.4598, Pressure: -13826140.0, Volume: 7503.4756
[MindSPONGE] Step: 400, E_pot: 1969.1597, E_kin: 229.6954, E_tot: 2198.855, Temperature: 192.64537, Pressure: -13705921.0, Volume: 7503.4756
[MindSPONGE] Step: 600, E_pot: 1941.4747, E_kin: 224.78912, E_tot: 2166.264, Temperature: 188.53047, Pressure: -13909719.0, Volume: 7503.4756
[MindSPONGE] Step: 800, E_pot: 1930.0374, E_kin: 207.89899, E_tot: 2137.9363, Temperature: 174.36473, Pressure: -14122749.0, Volume: 7503.4756
[MindSPONGE] Step: 1000, E_pot: 1916.6355, E_kin: 199.95618, E_tot: 2116.5918, Temperature: 167.7031, Pressure: -14373733.0, Volume: 7503.4756
[MindSPONGE] Step: 1200, E_pot: 1902.2065, E_kin: 203.48846, E_tot: 2105.695, Temperature: 170.66563, Pressure: -14358520.0, Volume: 7503.4756
[MindSPONGE] Step: 1400, E_pot: 1895.2224, E_kin: 218.30115, E_tot: 2113.5234, Temperature: 183.089, Pressure: -13955541.0, Volume: 7503.4756
[MindSPONGE] Step: 1600, E_pot: 1878.3259, E_kin: 198.46277, E_tot: 2076.7886, Temperature: 166.45056, Pressure: -14178873.0, Volume: 7503.4756
[MindSPONGE] Step: 1800, E_pot: 1880.648, E_kin: 200.18999, E_tot: 2080.838, Temperature: 167.8992, Pressure: -14218600.0, Volume: 7503.4756
[MindSPONGE] Finished simulation at 2023-09-05 16:08:00
[MindSPONGE] Simulation time: 2 minutes 6.1 seconds.
--------------------------------------------------------------------------------到这里我们还看不出什么比较特别的结论，但是我们可以从VMD可视化的结果中看看，相比于什么约束都不加，有什么样的效果。
可以看到，在上面这个动态图中，所有的水分子都在SETTLE约束算法的作用下保持键长键角不动，而中间的多肽因为不受SETTLE约束算法的影响，键长键角都在发生变化。
多肽Lincs约束
Lincs是一个只用于约束指定键的键长的约束算法，原理上相对复杂一点，可以参考下这篇博客中的介绍。
整体思路就像上面这张图一样，一根键在迭代器的作用下可能会发生偏移，此时在演化后的方向上把键长拉到指定的长度，这样就完成了约束算法。我们也可以理解为，每一个键都伸缩一下，最后再拼起来变成一个新的状态。这种约束算法，更多的适用于蛋白质的体系。一方面达不到TIP4P模型所需要的结果，另一方面Lincs算法本身是无法并行化和模块化的，性能也比不上SETTLE约束算法。因此对于水溶液体系里面一般是不用Lincs约束的。
配置Lincs约束算法的方法跟SETTLE算法是一致的，只需要在UpdaterMD中传入一个Lincs对象的constraint参数即可：
from mindspore import nn, context
from sponge import ForceField, Sponge, set_global_units, Protein, UpdaterMD
from sponge.function import VelocityGenerator
from sponge.control import Lincs
from sponge.callback import RunInfo, WriteH5MD
context.set_context(mode=context.GRAPH_MODE, device_target='GPU', device_id=1)
set_global_units('A', 'kcal/mol')
system = Protein('../pdb/case1.pdb', template=['protein0.yaml'], rebuild_hydrogen=True)
system.fill_water(edge=4, template='water.tip3p.yaml')
energy = ForceField(system, parameters=['AMBER.FF99SB', 'TIP3P'], use_pme=False)
min_opt = nn.Adam(system.trainable_params(), 1e-03)
md = Sponge(system, potential=energy, optimizer=min_opt)
run_info = RunInfo(200)
md.run(2000, callbacks=[run_info])
vgen = VelocityGenerator(300)
velocity = vgen(system.shape, system.atom_mass)
opt = UpdaterMD(system=system,
time_step=1e-3,
velocity=velocity,
integrator='velocity_verlet',
temperature=300,
thermostat='langevin',
constraint=Lincs(system, 'h-bonds'))
md.change_optimizer(opt)
nvt_h5md = WriteH5MD(system, 'test_nvt.h5md', save_freq=10, write_image=False)
md.run(2000, callbacks=[run_info, nvt_h5md])这里的Lincs约束仅仅施加在氢键上。由于这是一个多肽和水分子的混合体系，实际上这里的Lincs只施加在多肽上面的氢键（共价键）。上述代码的运行结果如下：
[MindSPONGE] Adding 45 hydrogen atoms for the protein molecule in 0.003 seconds.
[MindSPONGE] The box size used when filling water is: [21.184929 22.856928 15.495929]
[MindSPONGE] The edge gap along x axis is 4.0.
[MindSPONGE] The edge gap along y axis is 3.999999523162842.
[MindSPONGE] The edge gap along z axis is 3.999999761581421.
[MindSPONGE] Totally 172 waters is added to the system!
[MindSPONGE] Adding water molecule task finished!
[MindSPONGE] Started simulation at 2023-09-05 17:14:21
[MindSPONGE] Step: 0, E_pot: 6019.2305
[MindSPONGE] Step: 200, E_pot: 3667.4224
[MindSPONGE] Step: 400, E_pot: 3007.935
[MindSPONGE] Step: 600, E_pot: 2633.8376
[MindSPONGE] Step: 800, E_pot: 2380.7083
[MindSPONGE] Step: 1000, E_pot: 2233.6897
[MindSPONGE] Step: 1200, E_pot: 2129.8662
[MindSPONGE] Step: 1400, E_pot: 2050.8708
[MindSPONGE] Step: 1600, E_pot: 1993.4323
[MindSPONGE] Step: 1800, E_pot: 1949.2124
[MindSPONGE] Finished simulation at 2023-09-05 17:15:06
[MindSPONGE] Simulation time: 45.48 seconds.
--------------------------------------------------------------------------------
[MindSPONGE] The lincs constraint is used for the molecule system.
[MindSPONGE] Started simulation at 2023-09-05 17:15:07
[MindSPONGE] Step: 0, E_pot: 1910.498, E_kin: 515.0567, E_tot: 2425.5547, Temperature: 302.08228, Pressure: -1184091.6, Volume: 7503.4756
[MindSPONGE] Step: 200, E_pot: 2150.7756, E_kin: 402.22913, E_tot: 2553.005, Temperature: 240.10625, Pressure: -33554.688, Volume: 7503.4756
[MindSPONGE] Step: 400, E_pot: 2175.2153, E_kin: 444.4759, E_tot: 2619.6912, Temperature: 265.32498, Pressure: -80184.71, Volume: 7503.4756
[MindSPONGE] Step: 600, E_pot: 2186.1191, E_kin: 471.45215, E_tot: 2657.5713, Temperature: 281.42816, Pressure: -27345.148, Volume: 7503.4756
[MindSPONGE] Step: 800, E_pot: 2189.4526, E_kin: 504.2351, E_tot: 2693.6877, Temperature: 300.99756, Pressure: -75278.32, Volume: 7503.4756
[MindSPONGE] Step: 1000, E_pot: 2196.2683, E_kin: 522.76526, E_tot: 2719.0337, Temperature: 312.05896, Pressure: -38718.61, Volume: 7503.4756
[MindSPONGE] Step: 1200, E_pot: 2217.1887, E_kin: 489.45416, E_tot: 2706.6428, Temperature: 292.17426, Pressure: -59596.75, Volume: 7503.4756
[MindSPONGE] Step: 1400, E_pot: 2178.505, E_kin: 527.1349, E_tot: 2705.6396, Temperature: 314.66736, Pressure: -36722.918, Volume: 7503.4756
[MindSPONGE] Step: 1600, E_pot: 2161.2007, E_kin: 501.55945, E_tot: 2662.7603, Temperature: 299.4004, Pressure: -33997.46, Volume: 7503.4756
[MindSPONGE] Step: 1800, E_pot: 2185.1304, E_kin: 493.74838, E_tot: 2678.8787, Temperature: 294.73767, Pressure: -55893.01, Volume: 7503.4756
[MindSPONGE] Finished simulation at 2023-09-05 17:16:24
[MindSPONGE] Simulation time: 1 minutes 16.8 seconds.
--------------------------------------------------------------------------------从打印的输出结果中我们可以发现，这里lincs constraint被调用到了。同样的，我们可以用VMD插件对结果进行可视化：
从这个结果中我们可以看到，因为没有受到约束，水分子的键长键角都在不停的变化，多肽中的重原子共价键也在变化，只有多肽中的氢键长度保持不变。当然，如果这个案例中把代码调整为：
opt = UpdaterMD(system=system,
time_step=1e-3,
velocity=velocity,
integrator='velocity_verlet',
temperature=300,
thermostat='langevin',
constraint=Lincs(system, 'all-bonds'))那么就能够约束多肽里面所有的键，而水分子不进行约束。
全局氢键约束
上一个章节中提到了对多肽中的氢键用Lincs算法进行约束，再往上一个章节介绍了用SETTLE算法约束水分子，那么这里我们介绍如何混合使用两种不同的约束算法。还是同样的这个混合体系，我们在UpdaterMD中配置constraint的时候，不指定具体的Constraint约束算法，直接传入一个h-bonds，那么就会对全局的氢键（共价键）进行约束。这时候会出现一个问题，就是水分子和多肽分子中都存在氢键，那么哪些是使用SETTLE约束，哪些是使用Lincs约束的呢？有一个比较简单的等级区分，SETTLE约束算法的优先级要高于Lincs算法。所以，如果template配置文件中包含有SETTLE约束算法，那么在配置相应原子的约束条件的时候就会优先使用SETTLE算法。所以在这个案例中，水分子部分会被SETTLE算法约束，而多肽部分的氢键会被Lincs算法约束。
from mindspore import nn, context
from sponge import ForceField, Sponge, set_global_units, Protein, UpdaterMD
from sponge.function import VelocityGenerator
from sponge.callback import RunInfo, WriteH5MD
context.set_context(mode=context.GRAPH_MODE, device_target='GPU', device_id=1)
set_global_units('A', 'kcal/mol')
system = Protein('../pdb/case1.pdb', template=['protein0.yaml'], rebuild_hydrogen=True)
system.fill_water(edge=4, template='water.tip3p.yaml')
energy = ForceField(system, parameters=['AMBER.FF99SB', 'TIP3P'], use_pme=False)
min_opt = nn.Adam(system.trainable_params(), 1e-03)
md = Sponge(system, potential=energy, optimizer=min_opt)
run_info = RunInfo(200)
md.run(2000, callbacks=[run_info])
# 定义一个新的迭代器，并用change_optimizer完成切换
vgen = VelocityGenerator(300)
velocity = vgen(system.shape, system.atom_mass)
opt = UpdaterMD(system=system,
time_step=1e-3,
velocity=velocity,
integrator='velocity_verlet',
temperature=300,
thermostat='langevin',
constraint='h-bonds')
md.change_optimizer(opt)
nvt_h5md = WriteH5MD(system, 'test_nvt.h5md', save_freq=10, write_image=False)
md.run(2000, callbacks=[run_info, nvt_h5md])代码运行结果如下：
[MindSPONGE] Adding 45 hydrogen atoms for the protein molecule in 0.002 seconds.
[MindSPONGE] The box size used when filling water is: [21.184929 22.856928 15.495929]
[MindSPONGE] The edge gap along x axis is 4.0.
[MindSPONGE] The edge gap along y axis is 3.999999523162842.
[MindSPONGE] The edge gap along z axis is 3.999999761581421.
[MindSPONGE] Totally 172 waters is added to the system!
[MindSPONGE] Adding water molecule task finished!
[MindSPONGE] Started simulation at 2023-09-05 17:42:44
[MindSPONGE] Step: 0, E_pot: 6019.2305
[MindSPONGE] Step: 200, E_pot: 3667.4216
[MindSPONGE] Step: 400, E_pot: 3007.9355
[MindSPONGE] Step: 600, E_pot: 2633.8372
[MindSPONGE] Step: 800, E_pot: 2380.7073
[MindSPONGE] Step: 1000, E_pot: 2233.6892
[MindSPONGE] Step: 1200, E_pot: 2129.8655
[MindSPONGE] Step: 1400, E_pot: 2050.8716
[MindSPONGE] Step: 1600, E_pot: 1993.4336
[MindSPONGE] Step: 1800, E_pot: 1949.2134
[MindSPONGE] Finished simulation at 2023-09-05 17:43:30
[MindSPONGE] Simulation time: 46.60 seconds.
--------------------------------------------------------------------------------
[MindSPONGE] The settle constraint is used for the molecule system.
[MindSPONGE] The lincs constraint is used for the molecule system.
[MindSPONGE] Started simulation at 2023-09-05 17:43:31
[MindSPONGE] Step: 0, E_pot: 1910.4988, E_kin: 514.44617, E_tot: 2424.9448, Temperature: 301.7242, Pressure: -1060987.0, Volume: 7503.4756
[MindSPONGE] Step: 200, E_pot: 2022.5985, E_kin: 230.11551, E_tot: 2252.714, Temperature: 197.94637, Pressure: -13623341.0, Volume: 7503.4756
[MindSPONGE] Step: 400, E_pot: 1955.0367, E_kin: 234.73169, E_tot: 2189.7686, Temperature: 201.91722, Pressure: -14059788.0, Volume: 7503.4756
[MindSPONGE] Step: 600, E_pot: 1940.77, E_kin: 225.15894, E_tot: 2165.929, Temperature: 193.6827, Pressure: -13967141.0, Volume: 7503.4756
[MindSPONGE] Step: 800, E_pot: 1921.0624, E_kin: 215.5067, E_tot: 2136.569, Temperature: 185.3798, Pressure: -14294088.0, Volume: 7503.4756
[MindSPONGE] Step: 1000, E_pot: 1901.498, E_kin: 221.8313, E_tot: 2123.3293, Temperature: 190.82025, Pressure: -14622425.0, Volume: 7503.4756
[MindSPONGE] Step: 1200, E_pot: 1894.645, E_kin: 206.15628, E_tot: 2100.8013, Temperature: 177.33653, Pressure: -14258592.0, Volume: 7503.4756
[MindSPONGE] Step: 1400, E_pot: 1878.4142, E_kin: 205.46025, E_tot: 2083.8745, Temperature: 176.73781, Pressure: -14435189.0, Volume: 7503.4756
[MindSPONGE] Step: 1600, E_pot: 1877.7242, E_kin: 209.3739, E_tot: 2087.0981, Temperature: 180.10434, Pressure: -14059375.0, Volume: 7503.4756
[MindSPONGE] Step: 1800, E_pot: 1865.3273, E_kin: 215.8046, E_tot: 2081.1318, Temperature: 185.63606, Pressure: -14290035.0, Volume: 7503.4756
[MindSPONGE] Finished simulation at 2023-09-05 17:45:46
[MindSPONGE] Simulation time: 2 minutes 14.8 seconds.
--------------------------------------------------------------------------------可以看到，在可视化的结果中，水分子的键长键角是保持不动的，多肽中的氢键长度是保持不动的，但是多肽中的重原子共价键因为不受约束，所以也一直在动。
假如我们希望用Lincs算法来约束水分子中的氢键，那么可以在导入模板的时候，导入一个没有SETTLE参数的template，这样就可以实现用Lincs去约束水分子。但是需要注意的是，一个水分子只有两根共价键，因此Lincs算法虽然可以约束水分子的两个氢键的键长，但是无法约束其键角参数，也因此不能与TIP4P模型配合使用。
全体共价键约束
跟上面一个案例比较类似的，上一个案例是约束整个体系中所有的氢键，这里直接约束所有的共价键。只需要把h-bonds改为all-bonds即可。
from mindspore import nn, context
from sponge import ForceField, Sponge, set_global_units, Protein, UpdaterMD
from sponge.function import VelocityGenerator
from sponge.callback import RunInfo, WriteH5MD
context.set_context(mode=context.GRAPH_MODE, device_target='GPU', device_id=1)
set_global_units('A', 'kcal/mol')
system = Protein('../pdb/case1.pdb', template=['protein0.yaml'], rebuild_hydrogen=True)
system.fill_water(edge=4, template='water.tip3p.yaml')
energy = ForceField(system, parameters=['AMBER.FF99SB', 'TIP3P'], use_pme=False)
min_opt = nn.Adam(system.trainable_params(), 1e-03)
md = Sponge(system, potential=energy, optimizer=min_opt)
run_info = RunInfo(200)
md.run(2000, callbacks=[run_info])
# 定义一个新的迭代器，并用change_optimizer完成切换
vgen = VelocityGenerator(300)
velocity = vgen(system.shape, system.atom_mass)
opt = UpdaterMD(system=system,
time_step=1e-3,
velocity=velocity,
integrator='velocity_verlet',
temperature=300,
thermostat='langevin',
constraint='all-bonds')
md.change_optimizer(opt)
nvt_h5md = WriteH5MD(system, 'test_nvt.h5md', save_freq=10, write_image=False)
md.run(2000, callbacks=[run_info, nvt_h5md])运行结果如下：
[MindSPONGE] Adding 45 hydrogen atoms for the protein molecule in 0.003 seconds.
[MindSPONGE] The box size used when filling water is: [21.184929 22.856928 15.495929]
[MindSPONGE] The edge gap along x axis is 4.0.
[MindSPONGE] The edge gap along y axis is 3.999999523162842.
[MindSPONGE] The edge gap along z axis is 3.999999761581421.
[MindSPONGE] Totally 172 waters is added to the system!
[MindSPONGE] Adding water molecule task finished!
[MindSPONGE] Started simulation at 2023-09-06 09:10:50
[MindSPONGE] Step: 0, E_pot: 6019.2305
[MindSPONGE] Step: 200, E_pot: 3667.4207
[MindSPONGE] Step: 400, E_pot: 3007.9258
[MindSPONGE] Step: 600, E_pot: 2633.8313
[MindSPONGE] Step: 800, E_pot: 2380.703
[MindSPONGE] Step: 1000, E_pot: 2233.685
[MindSPONGE] Step: 1200, E_pot: 2129.8577
[MindSPONGE] Step: 1400, E_pot: 2050.8628
[MindSPONGE] Step: 1600, E_pot: 1993.4271
[MindSPONGE] Step: 1800, E_pot: 1949.2087
[MindSPONGE] Finished simulation at 2023-09-06 09:11:36
[MindSPONGE] Simulation time: 45.36 seconds.
--------------------------------------------------------------------------------
[MindSPONGE] The settle constraint is used for the molecule system.
[MindSPONGE] The lincs constraint is used for the molecule system.
[MindSPONGE] Started simulation at 2023-09-06 09:11:37
[MindSPONGE] Step: 0, E_pot: 1910.496, E_kin: 491.1751, E_tot: 2401.6711, Temperature: 288.07568, Pressure: -1245873.8, Volume: 7503.4756
[MindSPONGE] Step: 200, E_pot: 2011.3496, E_kin: 225.65253, E_tot: 2237.0022, Temperature: 198.51881, Pressure: -13732667.0, Volume: 7503.4756
[MindSPONGE] Step: 400, E_pot: 1960.8837, E_kin: 220.33536, E_tot: 2181.219, Temperature: 193.841, Pressure: -14297195.0, Volume: 7503.4756
[MindSPONGE] Step: 600, E_pot: 1923.8737, E_kin: 220.88599, E_tot: 2144.7598, Temperature: 194.32542, Pressure: -13654464.0, Volume: 7503.4756
[MindSPONGE] Step: 800, E_pot: 1913.1719, E_kin: 193.6853, E_tot: 2106.8572, Temperature: 170.39551, Pressure: -14154988.0, Volume: 7503.4756
[MindSPONGE] Step: 1000, E_pot: 1888.8213, E_kin: 197.11978, E_tot: 2085.9412, Temperature: 173.417, Pressure: -13763211.0, Volume: 7503.4756
[MindSPONGE] Step: 1200, E_pot: 1882.9513, E_kin: 208.48392, E_tot: 2091.4353, Temperature: 183.41466, Pressure: -14133313.0, Volume: 7503.4756
[MindSPONGE] Step: 1400, E_pot: 1861.6063, E_kin: 197.9073, E_tot: 2059.5137, Temperature: 174.10983, Pressure: -14121973.0, Volume: 7503.4756
[MindSPONGE] Step: 1600, E_pot: 1850.9318, E_kin: 218.16049, E_tot: 2069.0923, Temperature: 191.92766, Pressure: -14353631.0, Volume: 7503.4756
[MindSPONGE] Step: 1800, E_pot: 1856.4893, E_kin: 205.427, E_tot: 2061.9163, Temperature: 180.72531, Pressure: -14140336.0, Volume: 7503.4756
[MindSPONGE] Finished simulation at 2023-09-06 09:13:51
[MindSPONGE] Simulation time: 2 minutes 14.2 seconds.
--------------------------------------------------------------------------------在最后这个可视化结果中我们可以看到，所有的共价键的长度都被固定住了。不过这里面的水分子还是使用了SETTLE约束而不是Lincs约束。从打印的输出结果我们也可以看到，两个约束算法都有被调用到。
总结概要
本文主要介绍了在MindSponge中使用SETTLE和Lincs约束算法的方法，以及相关算法的简单原理。SETTLE约束算法主要应用于水分子体系，限制的是一个等腰三角形的拓扑结构，特点是可并行，性能较好。Lincs约束算法更多的被应用在蛋白质体系，主要限制的是每一个共价键的键长，特点是适用体系比较灵活，但总体计算量较大，且不可并行化。
版权声明
本文首发链接为：https://www.cnblogs.com/dechinphy/p/constraint.html
作者ID：DechinPhy
更多原著文章：https://www.cnblogs.com/dechinphy/
请博主喝咖啡：https://www.cnblogs.com/dechinphy/gallery/image/379634.html
手机扫一扫
移动阅读更方便
你可能感兴趣的文章