2014-05-09 17:45:40
无论是做何种类型的计算研究, 首要的工作就是建模. 对分子动力学模拟MD而言, 还要加上体系的力场化, 明确指出体系中的各种相互作用. 石墨烯的MD也不例外.
建模
只要熟悉晶体方面的知识, 石墨烯的建模不算复杂, 一种简单的方法可参考 建立石墨烯（Graphene）的模型. 这种方法构造出来的是六方结构, 用作MD不是很方便. 在六方结构上进行增删原子可得到四方结构, 但手动做起来有点麻烦. 所以还是写个简单的脚本来实现吧.
脚本实现的原理如下: 将石墨烯分解为含有四个C原子的基本单元, 再将基本单元在二维平面中排布.
若以C-C键长为单位, 基本单元的长宽分别为 $\sqrt 3, 3$, 其中的C原子坐标为 $(0,1/2) (0,5/2) (\sqrt 3/2, 1) (\sqrt 3/2, 2)$
几何性质
力场化之前需要清楚石墨烯的几何性质, 主要是原子个数, 键数, 键角数, 二面角个数之间的关系.
| 化学 | 几何 | 数目 | |
|---|---|---|---|
| 碳环数 | 六边形数 | 0.5N | |
| 原子数 | 顶点数 V | N | |
| 键数 | 棱数/边数E | 1.5N | |
| 键角数 | 角数 | 3N | |
| 二面角数 | 6N | ||
| 1-3相邻数 | 3N | ||
数算方法
- 每个原子有3条键, 每条键隶属于2个原子, 总键数为3N/2=1.5N
- 每个原子周围3个键角, 总键角数3N
- 每条键对应4个二面角, 总二面角数3N*2=6N
- 1-3相邻数目与键角数相同
从拓扑角度来说, 对闭曲面, 其顶点数 $V$, 棱数 $E$, 面数 $F$ 与 欧拉示性数 $\c$ 和 亏格 $g, k$ 之间存在下面的关系
\[V-E+F=\c =\begin{cases} 2(1-g) &可定向曲面 \\ \\ 2-k &不可定向曲面 \end{cases}\]石墨烯周期性体系与闭合的碳环面拓扑等价, 为可定向曲面, 其亏格 $g=1$, 因此 $V-E+F=0$.根据前面已知的关系
$F=E-V=V/2$
因此, 也可以根据六边形个数计算键数, 键角数和二面角数
- 每个六边形有6条棱, 每条棱隶属于2个六边形, 总棱数 $E=6F/3=3F=3V/2$
- 每个六边形6个角, 总角数 $6V/2=3V$
- 每个六边形相应于12个二面角, 总二面角数 $12V/2=6V$
力场化
- 准备力场参数可以根据原子之间的连接关系推算出所有的键, 键角和二面角, 还要注意周期性边界条件PBC的使用.
- 使用Gromacs做MD时, 对这种无限体系需要使用
periodic_molecules = yes选项, 否则计算有误. - 由于此选项, 体系很小时, 并行使用核数不能过多, 否则dd出错
