无机材料的基本操作

简单介绍了有关于无机材料(碳等)相关的知识。

序言

鄙人不才,在研究生阶段也研究过这方面的东西。看了许多文献,也实际重复过许多实验过程。但总感觉在这方面多为“业余”。故以此文记之,仅当乐趣罢了。

简介

碳材料现如今研究的非常多,不管在哪个领域都离不开它。我仅在分子动力学领域谈一谈。
首先,碳材料种类十分多,基本的碳材料包括石墨烯,碳纳米管和富勒烯,其他衍生出来的物质那就数不过来了。
但是在实际研究中基本上是那么几种。下面我就从最基本的建模开始谈起。

分子动力学

建模

现如今已经有许多的工具能够建立常见的碳材料模型,如下:

  • Avogadro,可以用于构建碳纳米管。
  • VMD,可以用来构建石墨烯和碳纳米管。
  • Nanotube Modeler,可以用来构建石墨烯,碳纳米管,富勒烯。
  • Materials studio,可以用来构建石墨烯,碳纳米管。
  • TubeASP,可以用于构建碳纳米管。
  • TubeGen,可以用于构建碳纳米管。
  • 石墨烯在线创建工具,可以用于构建石墨烯。
  • 纳米管在线创建工具,可以用于构建碳纳米管。
  • buildcstruct1.1,可以用于构建碳纳米管,石墨。
  • buildCstruct1.2,可以用于构建碳纳米管,石墨,以及在两端添加OH, COOH, COO-功能基团。相关文献点击
  • make_nanotube1.cc ,一段C++代码,可以产生碳纳米管(.gro)。
  • creatCNT.m,一段Matlab代码,可以产生碳纳米管(.xyz和.pdb)。
  • 合并多个pdb
    如果在vmd中构建了多个分开的分子模型,合并成一个pdb文件时可以采取以下命令:

    1
    ::TopoTools::mergemols [list 0 1]  //合并ID号为0,1的分子,返回合并后分子的ID

单独保存即可。

拓扑生成

  • gmxtop,基于opls力场的小分子拓扑生成工具,可以用于部分碳材料拓扑生成,需要检查整体电荷。
  • gmx x2top工具Gromacs模拟软件自带的拓扑生成工具,但是要注意首先需要在atomname2type.n2t文件中定义原子成键信息,ffnonbonded.itp定义LJ参数等等,可以参考链接
  • MKTOP,可以生成基于OPLS力场的小分子拓扑,不过需要自己提供原子电荷。
  • AmberTools+ACPYPE+Gaussian,用这种方式也可以创建碳材料拓扑。
  • 常用的方法举例如下:

    • gmx x2top
      建立atomname2type.n2t文件,如下:

      1
      2
      3
      4
      5
      C       CR1  0.0   12.011  3  C     0.142 C 0.142    C 0.142
      C CR1 0.0 12.011 2 C 0.142 C 0.142
      C CR1 0.0 12.011 1 C 0.142
      H HR1 0.115 1.008 1 C 0.109
      C CR1 0.0 12.011 3 C 0.142 C 0.142 H 0.108

      将n2t文件放入需要的力场文件夹下面。

    • 运行命令

      1
      gmx x2top -f CNT.pdb -o CNT.top -ff select -nopbc -name CNT -kb 400000 -kt 600 -kp 150
    • 主要问题
      gmx solvate给体系加水的时候有一部分的水分子会进入到碳管内部,这将影响整个模拟体系。因此需要通过Visual Molecular Dynamics(VMD)软件的选择语法对碳管周围以及内部水分子进行删除。

    • 语句说明
      • 选择残基名为SOL(即水分子)并且这些水分子在CNT周围3.5埃米的范围内
        set mov [atomselect top "same resid as (resname SOL and within 3.5 of resname CNT)"]
      • 移动选择的水分子向z轴移动5000埃米
        $mov moveby {0 0 5000}
      • 选择剩余的在z轴为1000埃米以下的区域部分,记为lef
        set lef [atomselect top "z<1000"]
      • 保存所需要部分的坐标为pdb格式
        $lef writepdb lef.pdb

周期性的考虑

因为有时候在模拟中要将单个石墨烯进行周期性的考虑,即考虑成无限大的石墨烯片。此时我们应该考虑的是模拟盒子的大小,类似于晶胞。可能一般人认为在GROMACS模拟软件中只需要将mdp中参数设置为periodic_molecules = yespbc=xyz不就好了。如果你这么考虑那就想的太简单了。虽然gmx x2top能够通过指定-pbc来对分子进行周期性的处理,即用于模拟无限大的碳材料(石墨烯的长宽和碳纳米管的长),但是如果你在建模的时候设置的盒子大小不能很好的匹配材料尺寸,那么你所谓的无限大或长是不正确的。比如通过vmd产生的Nanotube或者Graphene就存在这样的问题,具体参考链接1链接2链接3链接4
通过测试上面提及的建模软件,就石墨烯而言,只有Materials studio石墨烯在线创建工具这两项能够自动处理盒子大小,即它们所给出的尺寸参数正好是能够让周期性成像具有合适的距离,示意图如下:

通过Materials studio构建的石墨烯如下,分别在成像A、B、C与盒子中的石墨烯之间的原子距离恰好完整的构成正六边形,即看上去就像无限大的石墨烯片。

总而言之,如果建模软件不能给出合适的尺寸,你就应该自己通过模型去计算应该设置的盒子大小和位置,而不能仅仅通过gmx editconf -f -o -box x x x去随意设置大小,这种设置等同于上面给出的那个实例一样,没有考虑无限大或者长的体系。

其他

金纳米材料

这个可以见文后参考资料3实例。模拟这种元素体系要看你的研究内容了,如果你仅仅只是想看看这种材料与其他分子之间的物理相互作用,比如吸附等相关的性质,那么倒是可以通过GROMACS去做。但前提是有文献依据,即有相关的模拟参数。
如果你将金看做单个原子的组成,想要研究一下金纳米球或者其他各种形态与其他分子相互作用的关系,你可以这样做:

  • 构建金纳米球。首先通过MS构建一个足够大的超晶胞,导出pdb结构文件,然后通过VMD的选择语法截取一个球体。比如我想从超晶胞中截取一个球心为(3, 3, 3),半径为3 nm的金球,并保存为pdb文件用于GROMACS,可以通过如下命令:
    1
    2
    3
    4
    set mov [atomselect top "(x-30)^2+(y-30)^2+(z-30)^2 > 30^2"]
    $mov moveby {0 0 5000}
    set lef [atomselect top "z<1000"]
    $lef writepdb lef.pdb

得到金球lef.pdb文件。当然也可以直接选择小于半径的区域保存即可。

  • 添加溶剂。gmx solvate即可,如果你想加入其它分子啥的也是可以的。后面只需要改改拓扑就好了。
  • 然后需要构建拓扑文件,比如下面的:

    1
    2
    3
    4
    5
    6
    7
    8
    9
    10
    11
    12
    13
    14
    15
    16
    17
    18
    19
    20
    21
    22
    23
    #include "gromos53a6.ff/forcefield.itp"

    [ atomtypes ]
    AU 79 0.000 0.000 A 0.029247582 9.657102e-06

    [ moleculetype ]
    ; molname nrexcl
    AU 1

    [ atoms ]
    ; id at type res nr residu name at name cg nr charge mass
    1 AU 1 AU AU 1 0.0 196.9665

    #include "gromos53a6.ff/spce.itp"

    [ system ]
    ; Name
    Au sphere in water

    [ molecules ]
    ; Compound #mols
    AU 6642
    SOL 12956
  • 最后模拟即可。EMNVTNPT等等。。

    可以看到即使你没有对金添加位置限制或者固定,模拟过后整个金球也不会垮掉。这可能是因为金原子的相对分子质量大,半径大,使得金原子之间的范德华力十分强,因此可以维持原有形态。

二氧化硅界面吸附

这方面已经有非常多的文献研究了,不外乎先建立模型,然后进行拓扑文件的生成。此做法与石墨烯体系相似。比如我想研究一下水溶液中单层SiO2(0 1 0 晶面)对乙醇分子的吸附情况:

  • 建模型。最方便的就是用MS构建单层二氧化硅周期性模型。
  • 系统组合。可以用gmx内置工具,结合VMD,构建复合体系,比如像下面这样的模型:
  • 生成拓扑文件。有机小分子的拓扑就不讲了。如果我想研究周期性无限大的二氧化硅,那么首先还是像石墨烯那样建模,然后弄一个atomname2tpye.n2t文件,比如下面这样的(不同模型不一样):
    1
    2
    3
    4
    5
    6
    SI   SIR         0.31   28.0855 4    O  0.160   O 0.160   O 0.160   O 0.160  
    SI SIR 0 28.0855 4 O 0.160 O 0.160 H 0.147 H 0.147
    O OR 0 15.9994 2 SI 0.160 SI 0.160
    O OR -0.71 15.9994 2 SI 0.172 H 0.098
    H HR 0 1.008 1 SI 0.147
    H HR 0.4 1.008 1 O 0.098

然后通过gmx x2top生成拓扑文件即可,注意检查正确性。同时也需要将一些非键参数添加到ffnonbonded.itp中。这些参数一般可以参考文献来源。注意二氧化硅采用位置限制或者冻结(mdp参数),和无限大石墨烯相似。

  • 最后就是模拟了,然后分析结论。

实用命令

  • 如果你需要计算刚性石墨烯或者苯酚子平面法线与z坐标轴(或者分子平面与XY轴构成平面之间)的夹角,研究该分子在模拟过程中夹角的变化,可以通过gromacs自带的分析工具:gmx gangle
    程序通过分子中的三个原子坐标确定所在平面,首先应该将该三原子对添加到index.ndx文件中,这里记作[plane],然后继续下面操作:
    计算夹角随模拟时间的变化采用一下命令:

    1
    gmx gangle -f md_nopbc.xtc -s md.tpr -n index.ndx -g1 plane -group1 'group plane' -g2 z -oall angle.xvg -binw 1.0

    计算夹角在模拟过程中的占有率可采用一下命令:

    1
    gmx gangle -f md_nopbc.xtc -s md.tpr -n index.ndx -g1 plane -group1 'group plane' -g2 z -oh angle.xvg -binw 1.0

-g1后面选择平面,并通过三原子对确定分子所在平面的法线。
-g2后面选择的是z轴。
-binw指定角度输出间隔。

  • 选择碳管内部水分子命令(比如研究碳管中的受限水分子的行为)
    1
    2
    3
    gmx select -f -s -on <<EOD
    "Inside" same residue as ((x-1.5)^2+(y-1.5)^2 < 0.678^2 and z>1 and z<4 and resname SOL)
    EOD

表示碳管半径0.678 nm,圆柱体在xy平面圆心坐标(1.5,1.5);z方向选择1<z<2范围的分子;残基名称选择SOL;最终取交集,并通过same residue as保持分子完整,最终命名为组Inside。注意在gmx中长度单位均为nm,与vmd的埃米注意区分。

参考资料

[1] 使用Gromacs模拟碳纳米管的一个简单例子
[2] 基本VMD操作及碳纳米管和水体系示例
[3] SiO2界面吸附的分子动力学模拟_熊勇
[4] 在Gromacs中模拟金纳米线拉伸过程