Eastsheng's Wiki

MolTwister安装与使用

2023-10-21 18:11:28

[toc]

介绍

  • 用于分子系统构建、MD模拟和MD数据的后处理软件。

环境配置

Linux: Ubuntu-22.04

1
2
3
4
5
6
7
sudo apt install git
sudo apt install gcc
sudo apt install g++
sudo apt install cmake
sudo apt install libopengl-dev
sudo apt-get install libreadline-dev
sudo apt-get install libglu1-mesa-dev freeglut3-dev mesa-common-dev

Clone、编译、安装

1
2
3
4
5
6
git clone https://github.com/MolTwister/MolTwister.git
cd MolTwister
cmake CMakeLists.txt -DINCLUDE_CUDA_COMMANDS=OFF
#cmake CMakeLists.txt -DINCLUDE_CUDA_COMMANDS=ON
make
sudo make install

安装成功输出

1
2
3
4
5
6
7
Consolidate compiler generated dependencies of target MolTwister
[100%] Built target MolTwister
Install the project...
-- Install configuration: ""
-- Installing: /usr/local/bin/moltwister
-- Installing: /usr/local/bin/moltwister-tutorials/
-- Installing: /usr/local/bin/moltwister-tutorials//water_simulation.md

教程实操

  • 终端输入moltwister打开MolTwister命令行和右侧3D View

创建一个水分子PDB

  • 在MolTwister命令行执行以下命令：

    1
    2
    3
    4
    add atom H1 at 0 0 0
    add atom O from atom 0 dist 0.96
    add atom H2 from bond 0 1 angle 104.5 0.0 dist 0.96
    genbonds

    这将在(x, y, z)=(0,0,0)处创建一个名为H1的氢原子(其中前导H产生一个氢原子)，索引为0。然后加入一个氧原子，称为O，与索引为0的原子保持0.96埃的距离。最后，在索引为0 (H1)和1 (O)的原子之间以104.5度和0.0度的二面角形成键，与O的距离为0.96埃，加入一个氢原子H2

    如下图：

  • 测量距离和角度，检查一切是否如我们所料

    1
    2
    3
    measure bondlength id 0 1
    measure bondlength id 1 
    measure angle id 0 1 2

  • measure工具还可以测量没有相连的两个原子之间距离：

    1
    measure bondlength id 0 2

  • 要创建PDB文件内容，可以这样写:

    1
    2
    print pdb # 仅打印在screen
    print pdb > my_water.pdb # 写入文件

创建liquid water

  • 打开MolTwister，确保你有my_water.pdb(在上述教程中创建的pdb文件，或类似的pdb文件)。执行如下命令:

    1
    2
    3
    4
    5
    load my_water.pdb
    sel all
    copy sel 10 10 10 3.0 3.0 3.0
    autoscale
    genbonds

    这将把PDB文件加载到MolTwister中，选择水分子的所有原子，然后在x, y和z方向上复制该选择10次，每个方向上的距离为3埃。最后，自动缩放用于放大整个系统，然后为每个水分子创建键。

创建MolTwister脚本

  • 创建一个名为gen_my_water_mol.script的文件。编写脚本，编写以下内容:

    1
    2
    3
    4
    exec("add atom H1 at 0 0 0");
    exec("add atom O from atom 0 dist 0.96");
    exec("add atom H2 from bond 0 1 angle 104.5 0.0 dist 0.96");
    exec("genbonds");

  • 在MolTwister中运行以下命令:

    1
    load gen_my_water_mol.script

    这将通过依次调用文件中给出的每个命令来创建一个水分子。

创建一个MolTwister Python脚本

  • 创建一个名为gen_cube.py的Python脚本，其中包含(required Python3):

    1
    2
    3
    4
    5
    6
    7
    8
    9
    10
    11
    import moltwister as mt

    for i in range(0, 10):
        for j in range(0, 10):
            for k in range (0, 10):
                x = i*3.0
                y = j*4.0
                z = k*1.5
                mt.mt_exec(f"add atom O at {x} {y} {z}")

    mt.mt_exec("autoscale")

  • 然后在MolTwister中执行如下命令:

    1
    load gen_cube.py

Water simulation, thermalization, RDF, MSD, VACF, and VDOS

引言

  • 在本教程中，我们将创建一个单一的水分子，然后复制它来创建一个水相。随后，我们分配了力场参数，进行了能量优化，然后进行了仿真

  • 一旦模拟完成，我们将绘制温度作为时间的函数来确定热化点。接下来是计算

    • 径向分布函数，
    • 均方位移和自扩散，
    • 速度自相关函数，
    • 以及振动态密度。
创造一个水分子
1
2
3
4
5
6
7
8
9
10
11
add atom Hw1 at 0 0 0
add atom Ow at 0.9572 0 0
add atom Hw2 from bond 0 1 angle 104.52 0 dist 0.9572
# 建立边界
mod userdefpbc -10 10 -10 10 -10 10
set userdefpbc on

autoscale
genbonds
print pdb > single_water.pdb
print xyz > single_water.xyz
分配力场参数

  • 我们现在可以为新创建的水分子分配力场参数。

  • 我们可以直接输入命令，但我们可能需要在研究过程中多次分配力场。因此，我们创建一个名为water_ff.script的文件。该文件应包含：

    1
    2
    3
    4
    5
    6
    7
    8
    9
    10
    11
    # water_ff.script-->TIP3P water

    mod mass name Hw1 to 1.00784;
    mod mass name Hw2 to 1.00784;
    mod mass name Ow to 15.999;

    add mdbond Ow Hw1 Harm all_r_less_than 1.1 1882.8 0.9572;
    add mdbond Ow Hw2 Harm all_r_less_than 1.1 1882.8 0.9572;
    add mdangle Hw1 Ow Hw2 Harm all_r_less_than 1.1 230.12 104.52;

    add mdnonbonded Ow Ow LJ 0.426768 3.188;
清除所有内容

  • 很多时候，我们希望在MolTwister中执行子任务，最终得到一个脚本、pdb文件或类似的文件，可以稍后重新加载。在创建这些之后，我们通常想要从MolTwister内存中删除我们所做的一切。这分两步完成。

  • 首先我们删除我们在3D视图中看到的:

    1
    2
    sel all
    del sel

  • 随后，我们删除所有力场参数

    1
    del ff

  • 在某些系统上，还可以打开几个MolTwister实例，这在某些情况下可能很有用。

调查一下力场

  • 注意，要使下面的Python脚本工作，必须在系统上安装MatplotLib，这可以通过调用来完成: pip install matplotlib

  • 要使MolTwister工作，必须安装Python，所以系统上应该已经安装了Python。

  • 将所需的数据再次加载到系统中:

    1
    2
    load single_water.pdb
    load water_ff.script

  • 我们可以输出数据来显示指定力场电位的图，但首先我们需要知道指定力场的指数，这可以通过调用:

    1
    2
    3
    4
    list ff
    help moldyn
    # 我们可以看看水氧之间的伦纳德-琼斯势，也就是指数为0的非成键势
    moldyn ff nonbondforceprofile 0 2.9 10.0 100 > OwOw_LJ.dat

  • 其中结果存储在OwOw_LJ.dat文件中。这可以使用以下Python脚本绘制，名称为plot_non_bonded.py:

    1
    2
    3
    4
    5
    6
    7
    8
    9
    10
    11
    12
    13
    14
    15
    16
    17
    18
    19
    20
    21
    import matplotlib.pyplot as plt

    listX = []
    listY = []
    f = open("OwOw_LJ.dat", "r")
    startOfList = False
    for line in f:
        if startOfList:
            s = line.split()
            if len(s) >= 3:
                listX.append(float(s[0])) # Extract r in Å
                listY.append(float(s[2])) # Extract U in kJ/mol

        if "Potential" in line:
            startOfList = True

    fig, ax = plt.subplots()
    ax.plot(listX, listY)

    ax.set(xlabel='r [Å]', ylabel='U [kJ/mol]', title='Ow-Ow Lennard-Jones potential')
    plt.show()

    运行脚本画出LJ势能图：

创建liquid water

  • 现在我们要使用single_water.pdb文件，创建并构建了一个包含512个水分子的水相，随机放置，没有碰撞。假设密度为$$997 kg/m^3$$，两个氢的质量为1.0078，一个氧的质量为15.999，我们先计算一下盒子的大概尺寸:

    1
    2
    3
    4
    5
    6
    7
    8
    9
    10
    11
    calculate volumefromdensity 997 1.0078,1.0078,15.999 512

    # output:
    #        Mass of atom 0 of single molecule = 1.0078 g/mol
    #        Mass of atom 1 of single molecule = 1.0078 g/mol
    #        Mass of atom 2 of single molecule = 15.9990 g/mol
    #
    #        Target density = 997.0000 kg/m^3
    #        Number of molecules = 512
    #        Destination volume = 15362.3861 Angstrom^3
    #        Destination length for cubic system = 24.8591 Angstrom

    对于一个立方模拟盒，这产生的边长约为25 Å

  • 我们现在应该确保我们从一个干净的MolTwister开始:

    1
    2
    3
    sel all
    del sel
    del ff

  • 首先我们加载水分子:

    1
    load single_water.pdb

  • 然后说del来删除之前的东西。随后，选择所有原子，将其随机放置在一个边长为25 Å，分子原子之间最小距离为1.5 Å的立方盒子中，并制作512份副本:

    1
    2
    sel all
    copy sel random 25 25 25 512 1.5

  • 检查我们是否通过选择所有的水氧(Ow)并计算它们来创建512个分子:

    1
    2
    3
    sel none
    sel atomname Ow
    measure count sel

  • MD模拟器只接受几何中心位于原点的系统。因此，我们需要通过调用将系统替换到这个位置

    1
    2
    mod atompos all geomcent 0 0 0
    autoscale

  • 现在将系统保存到PDB文件中

    1
    print pdb > water.pdb
建立周期边界条件

  • 现在我们要为模拟设置适当的周期边界条件(PBC)。

  • 从上面的构造我们知道所有的原子都在25个Å的盒子里，但是我们也可以通过:

    1
    2
    3
    4
    5
    6
    7
    8
    measure pbc

    # output: 
    #        x = [-10.0000, 10.0000]
    #        y = [-10.0000, 10.0000]
    #        z = [-10.0000, 10.0000]
    #
    #        User defined PBC!

  • 然后，我们在PBC的每一边添加足够的量，以确保分子不会靠近PBC另一边的分子。因此，我们调用

    1
    mod userdefpbc -12.5 12.5 -12.5 12.5 -12.5 12.5

  • 设置一个25 × 25 × 25 Å边界框，我们调用

    1
    set userdefpbc on

  • 为了激活它。聪明的做法是创建一个名为pbc.script的文件。并插入最后两个命令，这样我们就不必每次都手动调整PBC，而是只需编写load pbc.script

进行能量优化

  • 分子被随机放置在我们创造的液相中，所以通过执行能量优化来最小化原子之间的力是一个好主意。

  • 我们将使用默认的最大优化步数和默认的学习率来实现这一点，但我们将设置当能量变化小于10 kJ/mol时算法将停止。

  • 要查看默认值，请调用:

    1
    moldyn cfg get

  • 为了设置10 kJ/mol的精度，调用

    1
    moldyn cfg set optaccuracy 10

  • 我们现在从一个新的MolTwister开始:

    1
    2
    3
    4
    5
    6
    sel all
    del sel
    del ff
    load water.pdb
    load pbc.script
    load water_ff.script

  • 然后，开始能量优化:

    1
    2
    moldyn cfg set outstride 10
    moldyn optimizeenergy

    直到它终止。这将生成一个名为traj.dcd的文件 ‘和out.txt。

  • 通过创建一个名为plot_energy.py的Python脚本，包含:

    1
    2
    3
    4
    5
    6
    7
    8
    9
    10
    11
    12
    13
    14
    15
    16
    17
    18
    19
    20
    21
    import matplotlib.pyplot as plt

    listX = []
    listY = []
    f = open("out.txt", "r")
    startOfList = False
    for line in f:
        if startOfList and "Searching" not in line:
            s = line.split()
            if len(s) >= 6:
                listX.append(float(s[0])) # Extract time step column
                listY.append(float(s[5])) # Extract energy column

        if "Step," not in line and "Step" in line:
            startOfList = True

    fig, ax = plt.subplots()
    ax.plot(listX, listY)

    ax.set(xlabel='Step', ylabel='U [kJ/mol]', title='Energy optimization')
    plt.show()

    我们可以看到优化过程的能量图:

  • 如果对结果感到满意，可以将优化过程的最后一帧导出到PDB文件。首先加载轨迹

    1
    load traj.dcd

  • 您可以通过单击3D视图并使用左右箭头键逐步查看优化过程。您还可以启用fog ，以更好的深度可视化

    1
    set fog on

  • 我们现在通过导出优化的帧(在本例中为帧499)

    1
    print pdb frame 499 > initial.pdb

    我们将使用它作为模拟的起始框架。

  • 如果您愿意，可以重命名该轨迹。DCD和out.txt文件从能源优化安全保管。它们将在后面的步骤中被删除。

执行分子动力学模拟

  • 我们重新开始清除MolTwister内存:

    1
    2
    3
    sel all
    del sel
    del ff

  • 现在删除所有剩余traj.dcd和out.txt文件

  • 我们重新加载初始系统，PBC和力场参数

    1
    2
    3
    load initial.pdb
    load pbc.script
    load water_ff.script

  • 我们配置为以0.5 fs的速度运行100000个时间步，总共产生50 ps，并且我们将输出步长设置为10帧:

    1
    2
    3
    moldyn cfg set timestep 0.5
    moldyn cfg set timesteps 100000
    moldyn cfg set outstride 10

  • 现在我们可以运行模拟

    1
    moldyn run

    在相对较新的CPU上，此模拟将花费几个小时并产生一个traj.dcd轨迹文件和一个out.txt输出文件。

    如果模拟失败，例如产生NaN (Not a Number)，那么可以尝试减少时间步长。例如，从0.5到0.2秒，以获得更好的平衡。

  • 一旦模拟完成，我们可以加载轨迹文件通过

    1
    load traj.dcd

  • 并利用键盘的左右箭头进行仿真分析。调用帮助查看如何调整每次按键的步数，以及如何显示/隐藏跨边界的键，以及在透视图和正射影视图之间切换。

绘制温度随模拟时间的函数

  • 温度作为时间步长的函数可以在out.txt中找到，我们将使用Python绘制结果。创建文件plot_temp.py，并输入以下脚本:

    1
    2
    3
    4
    5
    6
    7
    8
    9
    10
    11
    12
    13
    14
    15
    16
    17
    18
    19
    20
    21
    import matplotlib.pyplot as plt

    listX = []
    listY = []
    f = open("out.txt", "r")
    startOfList = False
    for line in f:
        if startOfList and "Searching" not in line:
            s = line.split()
            if len(s) >= 3:
                listX.append(float(s[0])) # Extract time step column
                listY.append(float(s[2])) # Extract temperature column

        if "Timestep," not in line and "Timestep" in line:
            startOfList = True

    fig, ax = plt.subplots()
    ax.plot(listX, listY)

    ax.set(xlabel='Step', ylabel='T [K]', title='Temperature')
    plt.show()

    温度vs时间：

    这应该能清楚地显示系统的热化周期。这对于确定在模拟中的哪个点找到可以执行计算的热化数据非常重要。在本教程中，我们将假设系统在8000个时间步长之后热化，这对应于DCD文件中的帧索引800，因为我们每10个时间步长输出一次。

Tags: MolTwister
← 论文笔记2023-冬 论文笔记2024-A →