分子模拟论坛 Molecular Simulation Forums

 找回密码
 立即注册

QQ登录

只需一步,快速开始

搜索
查看: 3335|回复: 10

富勒烯结构生成:CaGe程序及其改进

[复制链接]

173

主题

2919

帖子

6017

积分

论坛元老

Rank: 8Rank: 8

积分
6017
发表于 2010-9-7 10:05:00 | 显示全部楼层 |阅读模式
本文是bay 的本科论文的一部分, 两年来一直藏着
想来"保密期"已经过了, 就拿出来和大家分享分享
首发自木虫和本站共建的rcs 期刊第三期, 在这里开贴用于讨论和答疑
如果大家觉得bay 的这点工作有用, 请引用下面这篇文献
Li-Hua Gan*, Jian Liu, Qun Hui, Song-Qing Shao, Zuo-Hua Liu
General geometrical rule for stability of carbon polyhedra
Chemical Physics Letters, 2009, 472, 224–227,pdf
附件为本贴内容, 可以直接下载
也可在线查看下面诸楼层
回复

使用道具 举报

173

主题

2919

帖子

6017

积分

论坛元老

Rank: 8Rank: 8

积分
6017
 楼主| 发表于 2010-9-7 10:08:00 | 显示全部楼层
本帖最后由 bay__gulf618 于 2010-9-9 17:59 编辑
CaGe程序及其改进
Liu Jian(bay__gulf618@sina.com)
School of Chemistry and ChemicalEngineering,
Southwest
University
, Chongtsing 400715, P.R.China

摘要:本文介绍了CaGe软件包中部分程序的功能和使用方法及其输出文件的格式,在此基础上完成了任意指定的三~八元环个数的非经典富勒烯的生成程序的设计,生成格式为高斯软件可读的*.gjf格式,本文给出了程序的详尽算法和详细的使用说明,并附有C源代码。
关键词:非经典富勒烯;CaGe软件包;对偶图
1引言
碳是人类认识最早的元素,早在没有文字记录的史前时代就从烟灰、木炭等物质中了解了碳元素。它也是我们自认为研究最透彻的元素之一,石墨中的碳原子是以非常大的平面网状结构连接在一起,层与层之间相互平行地重叠在一起;金刚石的碳原子相互键合在一起形成三维网状结构。而富勒烯具有一个封闭的相对低分子量的稳定结构,并且几乎能以无穷无尽的方式进行化学反应和化学修饰。
2富勒烯的结构
第一个被发现的富勒烯是C60,它是在1985年由小罗伯特.F.柯尔(Robert F. Curl JR) ,哈罗德.W.克罗托 (Harold W. Kroto) 和理查德.E.斯莫利 (RichardE. Smalley)
等人于德克萨斯州莱斯 (Rice) 大学里检测到的[1], 并确定其结构为截角正二十面体的封闭的高对称性的球体。因为这项发现,三人
一同
分享了1996年化学诺贝尔奖金的。

在富勒
烯中,五元环相邻会造成稳定性的降低,因为沿着环面的
8 原子上的共轭π电子流的反芳香性将使结构变得不稳定[2]。最小的稳定的富勒烯是具有Ih 对称性的C60和具有D5d对称性的C70,它们也正是最小的具有分离五元环结构的富勒烯[3]。这便是大家熟知的"分离五边形(IPR)"规则。
不具备IPR规则的富勒烯是存在的,在60以下的位置上仍能找到它们的质谱信号。

图1, Kroto 获得的C原子簇质谱
在不能满足IPR规则的条件下,富勒烯可以通过两种方式闭合:
a. 有部分五元环相邻;
b. 有其他大小的环,比如四元环。
如果有两组相邻五元环,可以通过增加两个碳原子而消除五元环相邻的情形,代价是要形成一个四元环。带有两组相邻五元环的和带有一个四元环的分子在能量上比较接近。我们把只具有六元环和五元环的富勒烯称作经典富勒烯。把具有除五,六元环之外含有其他环数的称作非经典富勒烯。

图论中的欧拉定律表明单纯六元环的平面是不能闭合的,只有向其中加入一些较小的环才能闭合。在多面体中,设顶点(原子)数v,棱(键)数e,i元环环数为fi,根据如下几条定理:
a)
欧拉图论第一定理:3v = 2e
b)
欧拉多面体定理:∑fi+ v = e + 2
c)
对偶图内的欧拉第一定理:i fi = 2e
我们有:
(6-i)fi= 12
只有满足此关系式的多面体才能闭合。
对于只含有三,四,五,六,七,八元环的多面体,此式为
3f3 + 2f4 + f5 – f7 – 2f8 = 12

因f5 ≥ 0,此式又可表示为
3f3 + 2f4 – f7 – 2f8 ≤12

各种结构的富勒烯坐标可以使用手工构造,但对于对称性不高的分子这个构造过程十分繁琐,且当异构体数目较多时不能全部构造出来。现有一些可以生成任

C原子个数的经典富勒烯坐标的算法,如螺旋算法(spiralalgorithm)[4,5],增长变换算法[6]。但尚没有可以直接生成带有指定(三~八元)环个数的非经典富勒烯坐标的程序。本文在CaGe 包的基础上增加了此功能,且作向了*.gjf 的格式转换。下面将在介绍一些必要的数学知识和CaGe 包的部分功能后给出程序设计的算法和使用方法,并在文章后附有全部源代码。文中涉及到化学工作者不熟悉的计算机使用方法和数学理论,本文尽量给出详细的说明,使得在阅读本文时不需要参阅其他书籍。
回复

使用道具 举报

173

主题

2919

帖子

6017

积分

论坛元老

Rank: 8Rank: 8

积分
6017
 楼主| 发表于 2010-9-7 10:11:00 | 显示全部楼层
本帖最后由 bay__gulf618 于 2010-9-7 10:12 编辑
3 图论基础
由几何的方式生成的原子坐标以及共价键间的键长、键角,由于没有经过精确的量子力学计算而没有意义。完美表达一个分子结构的数学工具是图,而不是坐标。
图是这样一类由点和线组合而成的一个数学模型,它表达的仅仅是顶点集合的一种二元关系,与顶点的位置,边的长短曲直无关。它并非几何图形、工程图或美术图画,而是一种拓扑性的数学结构。
研究这类图的数学分支叫图论,由伟大的
数学家欧拉于
1736年创建。现在是离散数学的一个重要而热门的分支。
一个顶点连接的边的数目称作顶点的度数。
图中顶点的个数称作图的阶。
如果两个顶点间有边相连,称这两个顶点相邻接。
3.1 距离
图中两点间最短的路径间的边数称作两点间的距离。这条路径沿用几何学的名词,亦称作测地线。

3.2 邻接表和邻接矩阵
邻接表是图的一种表示方法,它给出与每个顶点相邻接的全部其他顶点,一般来说没有顺序。(但下面将要介绍到的embed 程序要求这个顺序)
邻接矩阵是图的另一种表示方法,当我们用一定次序把全部的n个顶点排列后,既可以构造这样一个n×n 的二进制矩阵A,满足:
aij= 1(ij邻接) / 0(ij非邻接或ij)
邻接表和邻接矩阵都可以完美表达一个图,且可以方便的相互转换而不影响点的编号。

对于邻接矩阵,有如下几个精彩有趣的性质。
定理1:An(G) 中ij号元素是图G中从vi到vj 的长n 的路径的数目。
证明:对n用数学归纳法证明。n= 1时,An= A,由邻接表定义知成立。
假设n= k 时定理成立,考虑Ak+1= AkA, Ak+1中第ij号元素为
a(k+1)ij = a(k)i1 a1j + a(k)i2 a2j + ... + a(k)in anj;
在a(k+1)ij的右端,a(k)il alj 这一项中a(k)ij
由归纳法假设,表示从vi到vl的长k 的路径的条数,而alj=1 时表示从vl到vj有一条边,这时a(k)il alj
表示从vi走到vl(走了k步)再从vl走一步到vj的路径的条数,1 当alj=0时
表示从vi走到vl,从vl“中转”,再走一步到vj的道路不存在。可见a(k+1)ij
是从vi到vj的长k+1 的路径的总条数。证毕。
推论1:经过点vi 的三元环的数目等于a(3)ii/2
推论2: 当图中没有三元环时,经过点vi 的五元环的数目等于a(5)ii/2
只给出推论2的证明:在A5中,a(5)ii表示vi经过5条边到自身的路径数,路径可以是一个五元环,亦可以是经过三元环后再来回走一条边然后到原处。当图中没有三元环时,可以排除后一种情况。在每个五元环共有两种走法,环数的2倍即为总的从vi到vi的走法,也即a(5)ii/2。这一条推论可以用在判断富勒烯是否有相邻的五元环上,且可以确定五元环的位置,并请注意它对于其他大小的环不适用。


3.3
可平面图
如果能够将图在平面(或球面)上画出,使它的全部边互不相交,那么这个图就是可平面图。由于富勒烯都是球面的,下面的讨论如果没有特别指明,皆指可平面图。

3.4 对偶图 (dual graph)
给定平面图G,如图1所示,其面分别为f1,f2,...,fk, G 的对偶图G*
的构造如下:图G*
的结点集为 {f1,f2,fk}, 对于G 中每一条边,与其关联的面为fi和fj,我们在G*
中画一条边fifj。

通俗地讲,就是把G中的面对应成G*
中的顶点,把G 中的顶点对应成G*
中的面。G*
称作G 的对偶图,G 称作原始图 (primalgraph) 。
定理2:如果G 中面,边,顶点数分别为f,e,v,则在G*
中分别为v,e,f。
定理3:如果一个图中没有度数为一的顶点(称作叶子)或在同一个顶点连接两次的边(称作伪环),则G与G*
一一对应,且有(G*)*≌ G。
如果图中所有顶点度数相同,那么这个图称作正则图。R-正则图指图中所有顶点度数皆为r。 在我们要研究的富勒烯分子中,每个顶点的度数皆为三,是一个3-正则图,它的对偶图恰是一个全部由三角形构成的网。
定理4:在3-正则图中,f= v/2 + 2,f为面数,v为顶点数。
证明: 欧拉图论第一定理:3v = 2e
欧拉多面体定理:∑fi + v = e +2
推得:∑fi= v/2 + 2
回复

使用道具 举报

173

主题

2919

帖子

6017

积分

论坛元老

Rank: 8Rank: 8

积分
6017
 楼主| 发表于 2010-9-7 10:13:00 | 显示全部楼层
本帖最后由 bay__gulf618 于 2010-9-7 10:20 编辑
4 CaGe 软件包
CaGe 是一个专业的用于生成经典富勒烯和碳纳米管坐标的开源软件包 (OpenSource software package)[7]
,主体部分用c 和java 语言写成。目前尚只有Unix/Linux 和mac 版本。下载地址为http://www.mathematik.uni-bielefeld.de/~CaGe/.
包中的软件主要由用以生成邻接表的生成器(generator)和把邻接表转换为笛卡儿坐标的转换器组成。前者有fullgen,
tubetype, plantri等,后者是embed。除此之外还有一些简单的文本处理器如cat,tee等,那是Unix系统的标准配置,在每台机器上都有它们的完整版本,在这里用处不大。
下面介绍一种在无java 且无gui(图形界面) 的Unix/Linux 上安装使用CaGe 软件包的方法。

4.1
安装CaGe
从CaGe的网站上下载这个软件包CaGe-sources.tar.gz。
新建一个文件夹如cage,
解压CaGe-sources.tar.gz:
mkdir cage && cd cage && tar zxvf
../CaGe-sources.tar.gz

进入Generators 文件夹:cd Generators
将它们编译成可执行文件: make
把Generators目录放入PATH变量中:输入PATH=${PATH}:${PWD}
或将此语句加入~/.bashrc文件中使之可以在任意目录下使用。
如果没有报错,说明安装成功。安装时不需要使用root用户权限。

4.2 fullgen

fullgen 是fullerenegenerator (富勒烯生成器) 的简写,名如其意。
它的基本的语法是:
fullgen nipr code # > outfile
n为将要生成的富勒烯的碳原子数。IPR选项表示只生成满足“分立五边形”规则(IPR)的结构。">"符号为输出重定向。code 后的#表示接1-8数字,表示生成文件的类型。其中6 为生成人类可读的格式(也即ASCII格式)。其他类型说明不在CaGe 包中,而在另一个名为plantri的软件包中,见plantri 包(下面介绍)中fullgen-guide.txt文件,本文在此不作说明。
例子,C60富勒烯的结构。
首先,我们用fullgen生成邻接文件: fullgen 60ipr code 6 > codes
其中,生成文件为codes。文件格式为:
>>writegraph3dplanar
100018192
20001340
30002324
40003541
50004306
.....
除首行外,第1列为C原子编号,2-4列均为0,5-7列为与第1列相邻的原子编号。当分子图以平面图方式画出时,这三个数字代表的C原子按照相同的顺逆时针方向旋转。这个文件提供了富勒烯分子的邻接信息。除此之外,fullgen 还可以处理对称性要求,设定顶点数的起止范围等。
4.3 embed
embed是由邻接表生成笛卡尔坐标的小程序。embed 的本意是嵌入,这里指把邻接表表示的图形嵌入到一个平面或球面上。
它的基本语法为
embed -d # outfile
#为2或3, 其中为2 生成二维坐标(嵌入平面),为3 生成三维坐标(嵌入球面)。
接上例,我们用embed程序将图代码文件转成笛卡尔坐标:
embed –d 3 coord
生成的坐标输出在名叫“coord”的文件中。 其格式如下:
>>writegraph3d
11.01272542-3.31656263-0.0993744118192
20.59064728-3.11044795-1.418255921340
31.32384299-2.27018974-2.264608202324
40.41757863-1.56673538-3.067007933541
50.66659436-0.22913073-3.396866844306
......
除首行外,2-4列为原子笛卡儿坐标,其他各行/列的意义同上。
因fullgen 生成的所有结构都存储在同一个文件中,而embed 只能读取其中的第一个。当需要embed计算全部的结构时,需要对输入文件手工编辑或编写一个脚本或程序使这一过程自动化,在下一节中介绍的a2gjf 文件就实现了了这一功能。
程序生成的坐标仅有3位小数,不能满足某些软件的精度要求。本文经过修改使之小数位数增至8位。
4.4 tubetype
tubetype是一个生成碳纳米管的小程序。
使用语法为
tubetypeoffset1 offset2 [tube m] > outfile
offset 为表征碳纳米管结构的位移量,m为管径方向上上六元环的个数,默认值为4
例子:输出长度为10的SWNT(5,5)临接表
tubetype 55 tube 10 > codes
遗憾的是本程序的输出格式为二进制,不能输出ascii格式。
格式如下(十六进制):
3E 3E 70 6C 61 6E61 72 5F 63 6F 64 65 20 6C 65 3C3C 00
2E 00 05 00 02 00 06 00
00 00 09 00 01 00 03 00
00 00 24 00 02 00 04 00
......
除文件开头的说明(前18个字节:">>planar_codele外,其他使用两个字节存放一个整数。存储格式为planar_code le,见本节对plantri的介绍。在x86框架下,每一个整数的存储方式为littleendian(是指低地址存放最低有效字节)。
此程序的生成的结构中,某些C 原子仅与两个C 原子相邻接。其输出文件经处理后可以被embed 读取,但图形显示有错误。
4.5 plantri
plantri 本身是一个没有化学背景的纯数学软件,它的开发和维护者也不是cage 的作者。但由于这个程序的作用同cage有相近之处,被放入cage包中。它的最新版本可以在这里下载。
http://cs.anu.edu.au/~bdm/plantri/
有趣的是,在plantri包中,也包含有fullgen 文件,甚至关于fullgen
的详细的说明在这里而不是在它自己的包中。这里使用的plantri程序来自于cage 包。
plantri是 planetriangulations的缩写,它的最基本的作用是生成一个全部由三角形构成的平面网状图(triangulations),这种图形被在本文中称作三角形网。在上节关于对偶图的讨论可知,这个图形的对偶图就是一个3-正则的平面网,嵌入到球面上即是富勒烯的数学模型。
它的常用语法为
plantri n-m#c# [-a][-d] > outfile
n为原始图的顶点数,#为正整数,-m#为图中最小的顶点度数δ,-c#为图的连通度κ,Whitney证明出[8]:κ≤δ。
当处理带四元环的富勒烯时,使用-m4即可。-a表示输出文件为ascii格式,否则输出planar code格式。-d表示输出图形的对偶图,值得注意的是当指定这个选项时,所有其他选项针对的是它的原始图而不是它本身。
例子:比如输出不带三度结点的阶数为10 的三角形网
plantri 10-m4 –a > outfile
结果为
10 bcdef,afghc,abhid,acie,adijf,aejgb,bfjh,bgjic,chjed,eihgf
10 bcdef,afghc,abhd,achije,adjf,aejgb,bfjih,bgidc,dhgj,digfe
10 bcde,aefghic,abid,acije,adjfb,bejg,bfjh,bgji,bhjdc,dihgfe
......
a-j表示十个顶点,第一列为阶数,下面每个逗号间依次表示与a-j
相邻接的顶点编号。每一行表示一个图形。
另一个例子:输出不带三元环的阶数为20 的3-正则图,它的面数=20/2+
2,它的对偶图是没有三度结点的阶数为12 的三角形网
plantri 12 -m4 –d > outfile
格式如下(十六进制)
3E 3E 70 6C 61 6E61 72 5F 63 6F 64 65 3C 3C
14 02 03 04
00 01 05 06
00 01 07 08
......
这是一种名为planar_code的存储格式,是plantri默认的,也是cage包中最常用的格式。在此笔者详细介绍一下这种格式,因为在本文提供的程序将直接使用到它。
除文件开头的说明(前15个字节:">>planar_code外,其他每一个字节存放一个整数。顶点数的编号从1开始最大到255。文件中有一些字节的数值为0。在它之间的三个数字表示与给定编号的顶点相邻接的其他顶点的编号,值为0的字节表示间隔。间隔三个数字在图展成平面时按照同一个顺逆时针方向旋转。每个图形间亦使用值为0的字节作间隔。
表一:planar_code和 ascii 两种格式比较
Table1: A compare between the format “planar_code” and “ascii”



planar_code(默认)
ascii(带选项 -a)
文件格式
二进制
文本
每个顶点的表示方式
1个字节
一个字符+空格
顶点名称
1,2,3...
字母 'a','b',...
顶点间间隔
0值字节
逗号 ','
图形间间隔
0值字节
新行 '\n'
回复

使用道具 举报

173

主题

2919

帖子

6017

积分

论坛元老

Rank: 8Rank: 8

积分
6017
 楼主| 发表于 2010-9-7 10:16:00 | 显示全部楼层
5 程序改造
新的软件包是在CaGe 包中的allowed_deg.c
embed.c
plantri.c 三个文件的基础上,增加了四个文件:class.cpp, e2gjf.cpp, run.sh, Makefile。使它可以按照linux 下的标准方式发行
,安
装和使用。新的软件包是一个完整的,独立的软件包。
##############################################################
此部分暂不提供,如果需要请向作者索取
(bay__gulf618@sina.com)
##############################################################
回复

使用道具 举报

173

主题

2919

帖子

6017

积分

论坛元老

Rank: 8Rank: 8

积分
6017
 楼主| 发表于 2010-9-7 10:21:00 | 显示全部楼层
本帖最后由 bay__gulf618 于 2010-9-7 10:22 编辑
6 测试
可以通过得出的异构体数和查看每个结构的空间图形两个方面验证本程序的正确性。
(a)使用time run.shn 0 -1 0 0 > /dev/null命令,生成最小含四元环(或不含)的C原子数为n 的结构数目,并可显示运行时间。在n=30时在屏幕上的输出如下
plantri_ad 17 -F4F5F6
230 triangulationswritten to stdout; cpu=0.11 sec
plantri_ad 17 -F4F5F6-d
230 cubic graphswritten to stdout; cpu=0.10 sec
there are 230 in totel
real
0m0.299s

user
0m0.234s

sys
0m0.009s

由此方法得C30
含有四元环(或不含)的所有异构体数目为230,运行时间0.299s 。
表二: 不同大小的四元环富勒烯的异构体数目
Table 2: the number of isomers of various size of fullerene
C原子数
异构体数
运行时间/s
C原子数
异构体数
运行时间/s
20
23
0.023
34
536
1.636
22
32
0.026
36
820
4.613
24
59
0.035
38
1175
13.02
26
93
0.055
40
1735
37.15
28
152
0.105
42
2418
105.53
30
230
0.237
44
3443
299.35
32
374
0.607
46
4711
686.24
回复

使用道具 举报

173

主题

2919

帖子

6017

积分

论坛元老

Rank: 8Rank: 8

积分
6017
 楼主| 发表于 2010-9-7 10:23:00 | 显示全部楼层
7程序的不足
程序尚有几点不足:
(1) 不能处理对称性问题;
(2) 当原子数较大时,需要大量的机时;
(3)
不能处理环的邻接问题。但在2 的限制条件下,小分子的富勒烯只有极少除六元环外,各个大小的环间各不相邻的要求。


8版权说明
在新的软件包中:
embed.c 的维护者是lisken@Mathematik.Uni-Bielefeld.de;
allowed_deg.c的维护者是 GunnarBrinkmann(University of Bielefeld, gunnar@mathematik.uni-bielefeld.de)和 Ulrike vonNathusius;
plantri.c的维护者是GunnarBrinkmann 和Brendan McKay( Australian National University
bdm@cs.anu.edu.au
);
软件包其他部分由笔者(bay__gulf618@sina.com.cn)维护.
该软件包为自由软件(free software),并遵照GPL v2 (GNUGeneral Public License
v2)协议发放。在满足下列条件下可以自由阅读,修改和发布源代码[10]:
(1)保留对原作作者说明;
(2)不得侵犯他人自由阅读,修改和发布源代码的权利;
(3)代码本身不得用于商业用途。
其他请参见GPL v2 的详则。


参考文献

[1]H.W.Kroto,J.R.Heath, S.C.O'Brien, R.F.Curl, R.E.Smalley. Nature, 1985, 18,162-163.
[2]A.D.J.Havmet, J.Am.Chem.Soc.
1986, 108:319-321.
[3]H.W.Kroto, Nature1987, 329,529-31.
[4]D.E.Manolopoulos,J.C.May, S.E.Down, Chem.Phys.Lett.1991,181,105.
[5]G.Brinkmann, D. Franceus, P. W. Fowler , Jack E. Graver, Chem.Phys.Lett. 2006, 428, 386.
[6]G.Brinkmann, A. W. M. Dress, J. Algori.1997, 23, 345.
[7]G. Brinkmann, O. Delgado Friedrichs, S.Lisken, A. Peeters, N. Van Cleemput, MATCH Commun. Math. Comput. Chem.,63(3), pp. 533-552, 2010,
[8]Whitney.H.,Congruent graph and theconnectivity of graph. Am. J. Math.1932 54 150-168.
[9]P.W.Fowler, T.heine,D.E,Manolopous,D.Mitchell,G.Orlandi,R.Schmit,G.Seifertand F.Zerbetto. J
.Phys.Chem.

1996,100,6984-6991
.
[1
0
]
R. M. Stallman, http://www.gnu.org/licenses/gpl.html,2007
回复

使用道具 举报

39

主题

363

帖子

771

积分

高级会员

Rank: 4

积分
771
发表于 2010-9-7 10:54:00 | 显示全部楼层
不错,怪不得bay这么牛,本科时基础就这么好了!
回复

使用道具 举报

31

主题

152

帖子

341

积分

中级会员

Rank: 3Rank: 3

积分
341
发表于 2010-11-3 21:14:00 | 显示全部楼层
好啊,谢谢分享
回复

使用道具 举报

108

主题

385

帖子

884

积分

高级会员

Rank: 4

积分
884
发表于 2011-3-4 12:57:00 | 显示全部楼层
你太牛了,我看都看不懂
回复

使用道具 举报

您需要登录后才可以回帖 登录 | 立即注册

本版积分规则

Archiver|手机版|小黑屋|分子模拟论坛

GMT+8, 2021-11-28 04:41 , Processed in 0.077469 second(s), 26 queries .

Powered by Discuz! X3.4

© 2001-2017 Comsenz Inc.

快速回复 返回顶部 返回列表