首页/文章/ 详情

Martini 3.0教程:参数化一个新的小分子

4月前浏览5084

0. Introduction

本教程将讨论如何为新的小分子构建Martini 3.0拓扑。该教程基于Martini 2教程,但加入了Martini 3.0中的一些更改。建议阅读文献[1]

原文用的例子是1-乙基萘,并使用Gromacs 2019.x或更高版本。这里我们以2-萘酰胺为例,因为我要做的体系是2-萘酰胺(叉腰)

注意:本教程基于Martini 3的公开测试版(v.3.0.b.3.2),将在Martini 3的最终参数可用时立即进行更新。

1)生成原子参考数据

我们将需要原子参考数据来提取CG模型的成键参数。请注意,在本教程中,我们将需要所有氢原子来提取键长,因此请确保您的原子结构包含所有氢。

原文使用了LigParGen服务器[2]来获得原子级(全原子,AA)结构和力场拓扑的方法,当然也可以随意使用自己喜欢的原子力场。这里我改用我习惯的acpype程序(对应Amber力场)或ATB(对应GROMOS力场)来生成拓扑。

通过acpype程序,我们已经获得了2-萘酰胺的结构文件(gro)和Amber拓扑,分别名为NAP_GMX.gro,NAP_GMX.top和NAP_GMX.itp。习惯上将NAP_GMX.top重命名为topol.top。

进入该文件夹,将刚才生成的itp、top、gro文件复 制过来,并执行:

将NAP_GMX.itp的[ atomtypes ]字段剪切到amber99s b-ildn.ff/ffnonbond.itp的原子类型中,并删除topol.top中的[ default ]字段,添加两行:

    #include "amber99s b-ildn.ff/forcefield.itp"#include "amber99s b-ildn.ff/spc.itp"

    运行如下命令:

    gmx insert-molecules -ci NAP_GMX.gro -nmol 1 -box 3.8 3.8 3.8 -rot xyz -seed 0 -o NAP_box.gro
    gmx solvate -cp NAP_box.gro -p topol.top -o NAP_sol.gro
    gmx grompp -f em.mdp -c NAP_sol.gro -p topol.top -o em
    gmx mdrun -v -deffnm em
    gmx grompp -f eq.mdp -c em.gro -p topol.top -o eq
    gmx mdrun -v -deffnm eq
    gmx grompp -f run.mdp -c eq.gro -p topol.top -o md
    gmx mdrun -v -deffnm md

    上述命令运行能量最小化的命令,然后执行250 ps的NPT平衡,并执行50 ns的MD运行(检查脚本和各mdp文件以了解更多信息)。在这种情况下,使用的溶剂是水;但是,也可输入其它平衡好的溶剂盒子。应该选择一种代表分子将花费大部分时间的环境的溶剂。

    2)原子到CG粒子的映射

    映射(即,将分子拆分为由CG粒子描述的构建单元)是粗粒化的核心,并且依赖于经验、化学知识和反复的试验。将分子映射到Martini 3.0模型时,应遵循以下准则:

    • 仅考虑非氢原子来定义映射;

    • 避免在两个小粒子之间切断特定的化学基团(例如酰胺或羧酸根);

    • 保证分子的对称性;此外,希望尽可能保留AA结构的体积和形状;

    • 4:1、3:1和2:1映射的默认选项是Regular(R),Small(S)和Tiny(T)粒子;它们是线性片段的默认选项,例如,辛烷包含的两个4:1片段;

    • 就计算性能而言,R型粒子是最佳选择,其粒子大小可以合理地代表4:1线性分子。

    • T型粒子特别适合代表芳环的平面度;

    • S型粒子通常更好地模拟脂族环的“大块”形状。

    • 应优化粒子的数量,以使映射中的最大错配为原子结构中每10个非氢原子±1个非氢原子;

    • 完全分支的碎片通常应使用较小尺寸的珠子(合理的理由是分支基团的中心原子被掩埋,也就是说,它不暴露于环境中,从而减少了其对相互作用的影响);例如,一个新戊烷基团包含5个非氢原子,但是由于它是完全支链的,因此可以安全地将其建模为规则的珠子。

    注意,在Martini 3.0中,共轭的、厚原子的结构最好由Tiny(T)粒子描述。这样可以确保与堆积相关的特性与全原子数据一致[1]。因此,在这种情况下,萘部分的10个碳原子映射到5个T珠,如下图所示:

    剩下的酰胺基团,映射为S粒子,因为S粒子的大小适合描述3个非氢原子,且酰胺不宜被切断。粒子在图中也已编号,以供进一步参考。

    3)从原子模拟生成CG映射轨迹

    使用刚才创建的映射,现在将在1)中进行的模拟转换为CG分辨率。一种方法是创建一个Gromacs(AA到CG)索引文件,其中每个索引组代表一个珠子,并包含映射的原子数。

    不用从头开始手动创建索引文件,而是可以使用CGbuilder工具[3]获得AA到CG的索引文件。只需加载分子的原子pdb / gro文件,单击要第一个粒子包含的原子。如果改变主意,则再次单击以将其删除,单击“新珠子”按钮创建下一个珠子;完成后下载文件。实际上,该工具还允许根据所选的映射获得珠子的初始CG配置(gro文件)和CG-AA映射文件(map文件)。

    注意,CGbuilder要求原子不能以不同于1的权重贡献给特定的粒子,而在将原子结构映射到Martini时有时会需要这样做。在这种情况下,索引和/或映射文件应随后手动完善。

    在开始之前:关于Martini 2.x的一个重要变化是,在将原子结构映射到CG分辨率时,现在要考虑氢原子以确定珠子的相对位置。这应该反映在AA到CG索引文件中,也就是说,索引还应该包含氢(对CGbuilder而言,就是也要单击氢!)。一般规则是将某个氢原子映射到含有与其相连的非氢原子的珠子上。

    现在,可以尝试通过CGbuilder映射NAP.gro。

    映射后的结构文件如下,重命名为NAP_CG.gro。

    Generated with cgbuilder
    6
       1NAP     BB    1  -0.358  -0.002  -0.046
       1NAP    SC1    2  -0.088  -0.095  -0.003
       1NAP    SC2    3  -0.052   0.206   0.027
       1NAP    SC3    4   0.111   0.012   0.002
       1NAP    SC4    5   0.273  -0.183  -0.021
       1NAP    SC5    6   0.349   0.099   0.004
    10 10 10

    索引文件如下,重命名为mapping.ndx。

    [ BB ]
    21 22 20 17 18
    [ SC1 ]
    10 5 19
    [ SC2 ]
    2 1 11 12
    [ SC3 ]
    4 3
    [ SC4 ]
    14 6 7 15
    [ SC5 ]
    16 8 9 13

    将mapping.ndx文件复 制到当前目录。

    现在,我们考虑了氢原子,因为对AA结构的基于几何中心(COG)的映射,是在Martini 3.0中获得成键参数的标准程序[1]。因此,在将AA结构映射到CG分辨率时,我们需要考虑氢原子。由于gmx traj的异常行为(潜在的错误,请参见注释[4]),如果我们要坚持使用gmx traj(替代方法包括使用MDAnalysis Python库),我们需要在运行gmx traj之前进行一些修改,即我们首先需要创建一个原子结构相同,但原子质量全部相等(e.g.,全为1)的全原子的tpr文件。

    复 制NAP_GMX.itp为NAP_COG.itp,并将原子质量全部改为1。

    NAP_COG.itp文件内容(部分):

    [ atoms ]
    ;   nr  type  resi  res  atom  cgnr     charge      mass       ; qtot   bond_type
        1   ca     1   UNK    C4    1    -0.110000     1.00000 ; qtot -0.110
        2   ca     1   UNK    C5    2    -0.122000     1.00000 ; qtot -0.232
        3   ca     1   UNK    C6    3    -0.023000     1.00000 ; qtot -0.255
        4   ca     1   UNK   C11    4    -0.050000     1.00000 ; qtot -0.305
        5   ca     1   UNK   C12    5    -0.054000     1.00000 ; qtot -0.359
        6   ca     1   UNK   C10    6    -0.108000     1.00000 ; qtot -0.467
        7   ca     1   UNK    C9    7    -0.129000     1.00000 ; qtot -0.596
        8   ca     1   UNK    C8    8    -0.117000     1.00000 ; qtot -0.713
        9   ca     1   UNK    C7    9    -0.124000     1.00000 ; qtot -0.837
       10   ha     1   UNK   H12   10     0.157000     1.00000 ; qtot -0.680
       11   ha     1   UNK    H4   11     0.134000     1.00000 ; qtot -0.546
       12   ha     1   UNK    H5   12     0.136000     1.00000 ; qtot -0.410
       13   ha     1   UNK    H7   13     0.135000     1.00000 ; qtot -0.275
       14   ha     1   UNK   H10   14     0.138000     1.00000 ; qtot -0.137
       15   ha     1   UNK    H9   15     0.136000     1.00000 ; qtot -0.001
       16   ha     1   UNK    H8   16     0.134000     1.00000 ; qtot 0.133
       17    c     1   UNK     C   17     0.670701     1.00000 ; qtot 0.804
       18    o     1   UNK     O   18    -0.611101     1.00000 ; qtot 0.193
       19   ca     1   UNK    C3   19    -0.144600     1.00000 ; qtot 0.048
       20    n     1   UNK     N   20    -0.673001     1.00000 ; qtot -0.625
       21   hn     1   UNK    H1   21     0.312500     1.00000 ; qtot -0.313
       22   hn     1   UNK    H2   22     0.312500     1.00000 ; qtot 0.000

    复 制topol.itp为topol_COG.itp,并改为引入NAP_COG.itp。

    topol_COG.itp文件内容(部分):

    ; NAP_GMX.top created by acpype (Rev: 0) on Fri Jan  8 22:17:10 2021
    #include "amber99s b-ildn.ff/forcefield.itp"#include "amber99s b-ildn.ff/spc.itp"
    ; Include NAP_GMX.itp topology
    #include "NAP_COG.itp"
    [ system ]
    NAP in water
    [ molecules ]
    ; Compound        nmols
    NAP              1    
    SOL              1692

    准备好上述文件后,执行命令:

    echo 0 | gmx trjconv -f md.xtc -o AA-traj.whole.xtc -s md.tpr -pbc whole
    gmx grompp -f run.mdp -c eq.gro -p topol_COG.top -o AA-COG.tpr
    rm mdout.mdp # clean-up
    # 因为要划分为6个粒子, 所以用seq 0 5生成0到5的序列, 且指定-ng 6
    seq 0 5 | gmx traj -f AA-traj.whole.xtc -s AA-COG.tpr -oxt mapped.xtc -n mapping.ndx -ng 6 -com

    这将保证:

    1. 首先,确保AA轨迹是完整的,即目标分子在轨迹文件的一个或多个帧中未被周期性边界条件分割(gmx trjconv -pbc whole ...命令);

    2. 随后创建一个AA-COG.tpr,将在接下来的步骤(gmx grompp -p ...命令)中用于COG映射;

    3. 最后,将AA轨迹映射到CG分辨率:gmx traj -f ...命令将进行COG映射,因为它使用的是AA-COG.tpr

    4)创建初始CGitptpr文件

    CGitp文件的创建必须手动完成,尽管从现有itp文件中进行一些复 制粘贴可能有助于正确设置格式。您(可能)需要的不同指令是:

    [ moleculetype ]:一行包含分子名称和许多排除项。对于Martini,排除的标准数量为1。

    [ atoms ]每个原子一行,含有原子编号,粒子类型,残基编号,残基名,原子名,电荷组,电荷,质量。小分子的残基编号和残基名称将相同。原子名称可以自由选择。在Martini,每个粒子都有自己的电荷组。通常不指定质量,在这种情况下,质量取自粒子定义。

    [ bonds ]每个键一行,包含原子1,原子2,函数类型,平均长度,作用力常数。Martini的函数形式为1。可以将键长设置为在步骤5中获得的平均长度。应将力常数根据所获得的直方图的宽度来调整:较小的宽度表示较高的力常数,反之亦然(另请参见以下各节)。

    [ constraints ][ angles ][ dihedrals ][ exclusions ] 的正确格式见Gromacs手册。

    构建初步的CG itp文件,名为NAP_initial.itp。实际上,您甚至不需要真正关心成键项。查看第5步的开头,以了解需要哪些成键项。

    NAP-initial.itp的文件内容如下:

    ;;;;;; 2-NAPHTHAMIDE
    [ moleculetype ]
    ; molname    nrexcl
     NAP          1
    [ atoms ]
    ; nr type resnr residue atom cgnr charge mass
      1  SC2   0    NAP    BB    1    0
      2  TC4   0    NAP   SC1    2    0     45
      3  TC4   0    NAP   SC2    3    0     45
      4  TC4   0    NAP   SC3    4    0      0
      5  TC4   0    NAP   SC4    5    0     45
      6  TC4   0    NAP   SC5    6    0     45
    [ virtual_sitesn ]
    ; site funct  constructing atom indices
      4     1     2 3 5 6

    在继续下一步之前,我们需要一个CG tpr文件来从映射的轨迹生成键,角度和二面角的分布。为此,应再准备一个topol_CG.top,文件内容如下:

    #include "martini_v3.0.b.3.2.itp"
    #include "NAP_initial.itp"
    [ system ]
    One molecule
    [ molecules ]
    NAP                    1

    编写用于CG的mdp文件,名为martini.mdp:

    ;
    ; STANDARD MD INPUT OPTIONS FOR MARTINI 2.x
    ; Updated 15 Jul 2015 by DdJ
    ; Adaptation to em/NVT/etc. by RA
    ;
    ; for use with GROMACS 5
    ; For a thorough comparison of different mdp options in combination with the Martini force field, see:
    ; D.H. de Jong et al., Martini straight: boosting performance using a shorter cutoff and GPUs, submitted.
    title                    = Martini ('new' parameters for) NPT run
    ; TIMESTEP IN MARTINI
    ; Most simulations are numerically stable with dt=40 fs,
    ; however better energy conservation is achieved using a
    ; 20-30 fs timestep.
    ; Time steps smaller than 20 fs are not required unless specifically stated in the itp file.
    integrator               = md
    dt                       = 0.02  
    nsteps                   = 1000000 ; 30 ns
    nstcomm                  = 100
    comm-grps       =
    nstxout                  = 0
    nstvout                  = 0
    nstfout                  = 0
    nstlog                   = 1000
    nstenergy                = 100
    nstxout-compressed       = 1000
    compressed-x-precision   = 100
    compressed-x-grps        = System
    energygrps               = System
    ; NEIGHBOURLIST and MARTINI
    ; To achieve faster simulations in combination with the Verlet-neighborlist
    ; scheme, Martini can be simulated with a straight cutoff. In order to
    ; do so, the cutoff distance is reduced 1.1 nm.
    ; Neighborlist length should be optimized depending on your hardware setup:
    ; updating ever 20 steps should be fine for classic systems, while updating
    ; every 30-40 steps might be better for GPU based systems.
    ; The Verlet neighborlist scheme will automatically choose a proper neighborlist
    ; length, based on a energy drift tolerance.
    ;
    ; Coulomb interactions can alternatively be treated using a reaction-field,
    ; giving slightly better properties.
    ; Please realize that electrostVatic interactions in the Martini model are
    ; not considered to be very accurate to begin with, especially as the
    ; screening in the system is set to be uniform across the system with
    ; a screening constant of 15. When using PME, please make sure your
    ; system properties are still reasonable.
    ;
    ; With the polarizable water model, the relative electrostatic screening
    ; (epsilon_r) should have a value of 2.5, representative of a low-dielectric
    ; apolar solvent. The polarizable water itself will perform the explicit screening
    ; in aqueous environment. In this case, the use of PME is more realistic.
    cutoff-scheme            = Verlet
    nstlist                  = 20
    ns_type                  = grid
    pbc                      = xyz
    verlet-buffer-tolerance  = 0.005
    coulombtype              = cutoff
    coulomb-modifier         = Potential-shift-verlet
    rcoulomb                 = 1.1
    epsilon_r                = 15; 2.5 (with polarizable water)
    vdw_type                 = cutoff  
    vdw-modifier             = Potential-shift-verlet
    rvdw                     = 1.1
    ; MARTINI and TEMPERATURE/PRESSURE
    ; normal temperature and pressure coupling schemes can be used.
    ; It is recommended to couple individual groups in your system separately.
    ; Good temperature control can be achieved with the velocity rescale (V-rescale)
    ; thermostat using a coupling constant of the order of 1 ps. Even better
    ; temperature control can be achieved by reducing the temperature coupling
    ; constant to 0.1 ps, although with such tight coupling (approaching
    ; the time step) one can no longer speak of a weak-coupling scheme.
    ; We therefore recommend a coupling time constant of at least 0.5 ps.
    ; The Berendsen thermostat is less suited since it does not give
    ; a well described thermodynamic ensemble.
    ;
    ; Pressure can be controlled with the Parrinello-Rahman barostat,
    ; with a coupling constant in the range 4-8 ps and typical compressibility
    ; in the order of 10e-4 - 10e-5 bar-1. Note that, for equilibration purposes,
    ; the Berendsen barostat probably gives better results, as the Parrinello-
    ; Rahman is prone to oscillating behaviour. For bilayer systems the pressure
    ; coupling should be done semiisotropic.
    tcoupl                   = v-rescale
    tc-grps                  = System
    tau_t                    = 1.0  
    ref_t                    = 298  
    Pcoupl                   = parrinello-rahman
    Pcoupltype               = isotropic
    tau_p                    = 12.0  ;parrinello-rahman is more stable with larger tau-p, DdJ, 20130422
    compressibility          = 3e-4
    ref_p                    =    1
    gen_vel                  = no
    gen_temp                 = 298
    gen_seed                 = 473529
    ; MARTINI and CONSTRAINTS
    ; for ring systems and stiff bonds constraints are defined
    ; which are best handled using Lincs.
    constraints              = none
    constraint_algorithm     = Lincs

    最后在官网下载martini.itp(当前名为martini_v3.0.b.3.2.itp),然后运行脚本:

    seq 0 5 | gmx traj -f AA-traj.whole.xtc -s AA-COG.tpr -n mapping.ndx -oxt molecule.gro -ng 6 -com -e 0
    gmx grompp -f martini_v2.x_new_run.mdp -c molecule.gro -p topol_CG.top -o CG.tpr -maxwarn 1

    这将:

    1. 从轨迹中提取一帧(当然将其映射);

    2. 使用该帧,以及top和mdp(请参阅网站上的示例)文件,为分子创建CG.tpr文件。

    遵循《Martini 3.0 圣经》[1],您还可以根据化学构造单元指定粒子类型。有关粒子选择的详细讨论,请参阅本教程的最后一节。

    5)从CG映射轨迹生成目标CG分布

    我们需要从步骤3映射到CG的原子模拟中,测量我们想要在CG模型中使用的成键相互作用的长度/角度(bonds, constraints, angles, proper and improper dihedrals)。但是,我们需要定义哪些成键项?让我们回到绘图表中,确定应该在哪些粒子之间定义成键相互作用。

    5.1)关于CG模型的成键项的选择

    连接2-萘酰胺内的T粒子的键很可能非常刚性,也就是说,它们的分布将非常集中。这要求使用constraints[1]。将粒子模型保持在一起的“幼稚”方法是把所有粒子都约束起来(下图A):

    但是,由于约束算法要满足越来越多的连接约束,这种模型容易出现数值不稳定性。另一种选择是构建“铰链”模型,其中通过5个约束将4个外部粒子(图中的粒子2、3、5和6)连接起来,以形成“铰链”构造,而中心粒子(编号 4)被描述为虚拟位点(Virtual Site, VS),即位置完全由其构造粒子定义的粒子(上图B)。使用虚拟粒子不仅可以提高模型的数值稳定性,还可以提高性能。由于VS没有质量,4号粒子的质量应分摊给四个构造粒子,这样,粒子2、3、5和6的质量应分别为45(= 36 + 36/4)。

    然后可以通过两个键即1-2和1-3连接粒子1。需要两个improper项(1-3-6-2和2-1-3-4),以将粒子1保持在萘环上的正确位置。最后还需要使用improper项 2-3-6-5来保持萘环平坦。在这种情况下,所有T粒子之间的排斥也应适用,因为分子非常硬并且在这种情况下不需要分子内相互作用。

    5.2)索引文件和目标分布的生成

    一旦确定了成键项,就为每个键创建一个带有索引[bondX]的键索引文件,其中包含成对的CG粒子,例如:

    [bond1]
     1  2
    [bond2]
     1  4
    ...

    角度(带有三个CG粒子)和二面角(带有四个)的情况类似。编写脚本,为感兴趣的所有键,角度和二面体生成分布。对于2-萘酰胺,有七个键(5个约束和2个键)和三个二面角,如前所述。

    教程提供了5_generate_target_distr.sh脚本,可通过运行

    bash 5_generate_target_distr.sh

    自动调用轨迹,生成分布。如果要用于新的体系,记得修改其中的gmx命令。

    检查文件夹bond_mapped和dihedrals_mapped的结果,会发现每个成键分布为bonds_mapped / distr_bond_X.xvg,而bonds_mapped / data_bonds.txt 给出了每个键的均值和标准偏差。

    对于每个键,脚本使用以下命令(在此示例中,该命令应用于索引为0的第一个键):

    echo 0 | gmx distance -f mapped.xtc -n bonds.ndx -s CG.tpr -oall bonds_mapped/bond_0.xvg -xvg none
    gmx analyze -f bonds_mapped/bond_0.xvg -dist bonds_mapped/distr_bond_0.xvg -xvg none -bw 0.001

    对二面角也类似:

    echo 0 | gmx angle -type dihedral -f mapped.xtc -n dihedrals.ndx -ov dihedrals_mapped/dih_0.xvg
    gmx analyze -f dihedrals_mapped/dih_0.xvg -dist dihedrals_mapped/distr_dih_0.xvg -xvg none -bw 1.0

    6)创建CG模拟

    现在,我们可以完成CG模型的第一阶段NAP_take1.itp的确定,在该模型中,我们可以使用data_bonds.txt和data_dihedrals.txt文件中包含的信息来更好地猜测结合参数。

    复 制NAP_initial.itp为NAP_take1.itp。使用上一步的输入来调整ENAP_take1.itp。

    ;;;;;; 2-NAPHTHAMIDE
    [ moleculetype ]
    ; molname    nrexcl
     NAP          1
    [ atoms ]
    ; nr type resnr residue atom cgnr charge mass
      1  SC2   0    NAP    BB    1    0    
      2  TC4   0    NAP   SC1    2    0     45  
      3  TC4   0    NAP   SC2    3    0     45  
      4  TC4   0    NAP   SC3    4    0      0
      5  TC4   0    NAP   SC4    5    0     45  
      6  TC4   0    NAP   SC5    6    0     45  
    [bonds]
    ; i  j  funct length
     1  2   1     0.28079   25000 ; cog
     1  3   1     0.38768   25000 ; cog
    #ifndef FLEXIBLE
    [constraints] ; For minimizations, you may use define=-DFLEXIBLE to use a stiff-bond version of the topology that is more amenable to minimization.
    #endif
    ; i  j  funct length
     2  3   1     0.30060 1000000 ; cog
     2  5   1     0.36661 1000000 ; cog
     2  6   1     0.47223 1000000 ; cog
     3  6   1     0.41003 1000000 ; cog
     5  6   1     0.29058 1000000 ; cog
    [ dihedrals ]
    ; improper
    ; i j k l  funct  ref.angle force_k
     1 3 6 2    2        0        50  ; controls planarity of amide group
     2 1 3 4    2        0        50  ; keeps amide on position 1 of naphthalene
     2 3 6 5    2        0       200  ; keeps naphthalene planar
    [ virtual_sitesn ]
    ; site funct  constructing atom indices
      4     1     2 3 5 6
    [ exclusions ]
    1 2 3 4 5 6
    2 3 4 5 6
    3 4 5 6
    4 5 6
    5 6

    同时修改topol_CG.top,如下:

    #include "martini_v3.0.b.3.2.itp"
    #include "martini_v3.0_solvents.itp"
    #include "NAP_takel.itp"
    [ system ]
    One molecule
    [ molecules ]
    NAP                    1

    执行:

    gmx solvate -cp NAP_CG.gro -cs box_CG_WN_eq.gro -p topol_CG.top -o NAP_sol.gro -box 4.5 4.5 4.5
    gmx grompp -f martini_em.mdp -c NAP_sol.gro -p topol_CG.top -o 1-min.tpr -po 1-min.mdp  -maxwarn 1
    gmx mdrun -v -deffnm 1-min
    gmx grompp -p topol_CG.top -c 1-min.gro -f martini_eq.mdp -o 2-eq.tpr -po 2-eq.mdp
    gmx mdrun -v -deffnm 2-eq
    gmx grompp -p topol_CG.top -c 2-eq.gro -f martini_run.mdp -o 3-run.tpr -po 3-run.mdp
    gmx mdrun -v -deffnm 3-run

    将对水中的Martini系统进行能量最小化,然后进行NPT平衡,并执行50 ns的MD产生相模拟(检查脚本和各种mdp文件以了解更多信息)。

    检查所生成的结构文件和轨迹文件。

    一旦运行了MD,就可以使用为映射轨迹生成的索引文件来生成CG轨迹的分布:

    cp ../5_target-distr/bonds.ndx .
    cp ../5_target-distr/dihedrals.ndx .
    bash 6_generate_CG_distr.sh

    它将生成一系列文件,与上一步中的5_generate_target_distr.sh所做的相似,但现在针对CG轨迹。

    现在,可以绘制使用脚本draw_bond_distribution.py和draw_dihedral_distribution.py绘制二者的分布并进行比较。

    发现萘环上的键长分布与全原子吻合较好,而1-3和1-4(前两张子图)分布不太合理。

    观察模拟结构,发现酰胺基团和3号原子(SC2)距离较近,不合理。因此在NAP_take1.itp中添加键角约束。

    添加键角约束后重新跑一次能量最小化、平衡相和产生相模拟,再计算键长分布,如下图。

    可以看到两者分布比上一次模拟合理一些。

    同理计算二面角分布,如下图。

    注意,有时分布会出现两个峰,而CG模型无法考虑双峰性。如果在第一次迭代中不令人满意(很可能会发生),则应使用CG itp中的平衡值和力常数进行迭代,直到达成令人满意的协议。

    8)与实验结果进行比较,进一步完善和最终考虑

    粗粒化后的模拟结果应从划分自由能、纯液体或有机晶体的密度等进行合理性验证,以及其它更具体的参考材料。

    (注:本文章来源于网络,如有侵权,请联系我们删除。)

    来源:模拟之家
    ACTACPSystemSTEPS化学UG材料试验
    著作权归作者所有,欢迎分享,未经许可,不得转载
    首次发布时间:2024-01-13
    最近编辑:4月前
    刘十三613
    博士 分子动力学、GROMACS
    获赞 95粉丝 39文章 3课程 27
    点赞
    收藏
    未登录
    还没有评论

    课程
    培训
    服务
    行家

    VIP会员 学习 福利任务 兑换礼品
    下载APP
    联系我们
    帮助与反馈