GROMACS模拟实例--两水球碰撞

简单介绍一下如何使用GROMACS进行两水球的碰撞模拟。

模拟的世界总是很有乐趣。接下来我将介绍一下如何采用GROMACS进行两水球碰撞的模拟过程。

准备

  • 软件GROMACS
  • awk脚本

操作步骤

  • 创建水球
    这个很简单,直接使用GROMACS内置命令即可完成。首先我们来建立一个6 x 6 x 6 $nm^3$的立方体水盒子(TIP3P):

    1
    gmx solvate -cs spc216.gro -box 6 6 6 -o box.gro

    然后从中间挖取一个半径为2 nm的水球,这里直接用gmx select命令即可:

    1
    gmx select -s box.gro -ofpdb -pdbatoms selected -select 'same residue as ((x-3)^2 + (y-3)^2 + (z-3)^2 <= 2^2)'

    这里为了后面方便写awk脚本处理坐标。我使用gmx editconf命令转换上一步产生的水球坐标occupancy.pdb文件:

    1
    gmx editconf -f occupancy.pdb -o ball.gro
  • 创建碰撞初始体系
    上一步已经完成了创建水球的工作了。但是你需要两个水球来进行碰撞,下面将使用awk脚本进行水球复制和位置平移。当然如果你擅长vmd,那使用它也是极方便的。
    这一步包括复制水球坐标以及各自更改两个水球分子的原子速度,并构建一个初始形心距离为10 nm的两水球体系。下面是我针对该体系给出的一个awk脚本:

    1
    2
    3
    4
    5
    6
    7
    8
    9
    10
    11
    12
    13
    14
    15
    16
    17
    18
    19
    20
    21
    22
    23
    24
    25
    26
    27
    28
    29
    30
    31
    32
    33
    34
    35
    36
    37
    38
    39
    40
    41
    42
    43
    44
    45
    46
    47
    48
    49
    50
    51
    52
    53
    54
    55
    56
    57
    58
    #!/bin/bash
    # A awk program for generating initial system with two water ball from a water ball .gro file

    ############################################################################
    # Author : Written by Yujie Liu #
    # Email : 1051690424@qq.com #
    # Version : 1.0 #
    # Date : 2019.04.01 #
    ############################################################################

    #=========================HELP INFORMATION==============================#
    if [[ -z $1 ]] || [[ ! -e $1 ]]; then
    echo -e "\nError! No input file!!"
    echo -e "
    Usage:

    $0 <single_waterball.gro>
    "
    exit 1
    fi
    OutputFile="new_ball.gro"

    awk -v filename=$1 '
    BEGIN {
    ChangGro(filename, "left", 1)
    ChangGro(filename, "right", 0)
    }

    #++++++++++++++++++++++++++++MAIN FUNCTION++++++++++++++++++++++++++#
    # read this single .gro file, change the velocity of first water ball,
    # change the velocity and coordinate in X axis of the second water ball,
    # here the distance between two wate balls is 10 nm.
    # You can change it for your means

    function ChangGro(Filename, position, isfirst, _ARGVEND_, i) {
    getline < Filename; if(isfirst == 1) {print};
    getline < Filename; Natoms = $1+0; if(isfirst == 1) {print 2*Natoms}
    for(i = 0; i < Natoms; i++) {
    getline < Filename
    if(position == "left") {
    novelocity = substr($0, 1, 44)
    # give velocity (1, 0, 0) for left water ball
    printf("%s%8.4f%8.4f%8.4f\n", novelocity, 1, 0, 0)
    }

    else if(position == "right"){
    # give velocity (-1, 0, 0) for another water ball
    x = substr($0, 21, 8)+10 # plus 10 nm in X axis
    frontx = substr($0, 1, 20)
    afterx_and_frontvx = substr($0, 29, 44)
    printf("%s%8.3f%s%8.4f%8.4f%8.4f\n", frontx, x, afterx_and_frontvx, -1, 0, 0)
    }
    }
    getline < Filename; lx = $1+0; ly = $2+0; lz = $3+0
    if(isfirst == 0) {printf("%12.6f %12.6f %12.6f \n", lx+10, ly, lz)}
    close(Filename)
    }
    ' > ${OutputFile}

    执行上面脚本bash awk_genConfVelocity.bsh ball.gro即可得到名为new_ball.gro的初始体系,此时两水球的x方向相距10 nm,两水球各自的原子速度矢量被设置为(1, 0, 0)和(-1, 0, 0) [$nm/ps$]。表示两水球在x方向以1.0000 $nm/ps$的速度相撞。

  • 进行测试模拟
    接下来就可以进行撞击模拟了,由于这个纯水体系使用的平衡好的spc216.gro模型进行搭建的,因此不需要能量极小化这一步了。我们直接来进行NVT模拟即可,使用以下MDP文件:

    1
    2
    3
    4
    5
    6
    7
    8
    9
    10
    11
    12
    13
    14
    15
    16
    17
    18
    19
    20
    21
    22
    integrator = md
    dt = 0.001 ; ps
    nsteps = 10000
    comm-grps = system

    nstxout = 10
    nstvout = 10
    nstfout = 0
    nstlog = 1000
    nstenergy = 10
    nstxout-compressed = 0
    compressed-x-grps = system

    pbc = xyz
    rlist = 1.2
    nstlist = 1
    cutoff-scheme = group
    coulombtype = cut-off
    rcoulomb = 1.2
    vdwtype = cut-off
    rvdw = 1.2
    DispCorr = no

    ???忘记说top文件怎么产生了,直接命令:

    1
    gmx pdb2gmx -f new_ball.gro

    生成一个就行了,记得模拟初始结构还是使用new_ball.gro

    生成tpr然后跑MD即可,这里我们跑一个10 ps

    1
    2
    gmx grompp -f nvt.mdp -c new_ball.gro -o md.tpr
    gmx mdrun -deffnm md -v

    很快就完成了!

  • HAPPY PLAY
    接下来让我们看看运动轨迹情况,使用vmd md.gro md.trr打开轨迹查看一下。Interesting...

    动画视频1

    https://www.bilibili.com/video/BV1Wz4y1R72H/

拓展操作

上面的例子是在一条直线上的碰撞,那么可不可以以V字型在尖端进行碰撞呢??

答案是可以的,就直接更改上面的脚本中速度矢量分别为(1, 0, 1)和(-1, 0, 1),Z方向的盒子长度设置16 nm,表示在xz平面成V字型碰撞。同时也要关掉MDP中的消除质心移动选项,即使用:

1
comm-mode  =  None

动画视频2

https://www.bilibili.com/video/BV1DK41157i3/

总结

从动画视频中我们可以看出水球碰撞以后会发生破裂,如果是速度大小相等,方向相反的两水球碰撞,最终会形成一个有一定厚度的水层。
如果不是直线等速碰撞,沿V字型进行碰撞,结果会形成不均质水层,由于另一个方向速度叠加,该水层会继续速度叠加方向运行。

题外话

简单体系自己写写脚本进行建模还是挺方便的,不一定要依赖于现成的软件。