发信人: zjliu (秋天的萝卜), 信区: Chemistry
标 题: 计算量子化学通用引擎协议(草案)
发信站: BBS 哈工大紫丁香站 (Fri Apr 1 08:19:46 2005)
(网上有彩色文档：http://www.elephantbase.net/other/ucqci.htm)
《电脑象棋和量子化学》相关资料
计算量子化学通用引擎协议(草案)
黄晨 * 2005年3月
(* 联系地址：复旦大学化学系表面化学实验室，eMail：morning_yellow@elephantbase
.net)
引言
目前计算量子化学引擎发展得不太成熟，笔者认为其主要原因是没有统一的引擎协
议。为此笔者借鉴了国际象棋通用引擎协议的一些方法(请对比《国际象棋通用引擎协议
》一文)，制定了适合计算量子化学的通用引擎协议。由于计算量子化学涉及的功能和算
法众多，笔者不可能完全了解它们的细节，所以这篇文章仅仅是协议内容的一部分，它
所能处理的问题也只是计算量子化学所涉及的最基本的问题，包括：
(1) 计算某个化学粒子(原子、分子或离子)的基本性质，例如能量、能级、光谱、
电荷分布、磁屏蔽效应等，这种计算是量子化学最基本的计算，通常称为“单点计算”
；
(2) 计算分子的振动性质，从而获得分子的红外光谱和热力学函数，通常称为“频
率计算”；
(3) 优化分子结构，即寻找势能面上的极小值点；
(4) 寻找反应过渡态，即寻找势能面上的鞍点。
一、单位制
大部分计算量子化学公式都涉及到“多电子体系的Schr?dinger方程”，为简化方程
的形式，最好采用“原子单位制”，长度、时间、质量、能量和电量的单位分别是Bohr
、Heisenberg、Thomson、Hartree和Millikan，其中Heisenberg、Thomson和Millikan是
笔者给原子单位的命名。
注：R. A. Millikan不要同R. S. Mulliken混淆，尽管他们的中文译名比较接近。
前者因油滴实验准确测定电子电量而获得1923年诺贝尔物理学奖，后者因创立分子轨道
理论而获得1966年诺贝尔化学奖。电量单位“密利根”指的是前者，原子的“密利根”
电荷指的是后者。
原子单位同国际制单位和常用单位有如下的换算关系：
1 Bohr = 5.291772108(18)×10-11m；
1 Heisenberg = 2.418884326505(16)×10-17s；
1 Thomson = 9.1093826(16)×10-31kg = 5.4857990945(24)×10-4amu；
1 Hartree = 4.35974417(75)×10-18J = 27.2113845(23)eV；
(1 Hartree = 1 Thomson·Bohr2/Heisenberg2)；
1 Millikan = 1.60217653(14)×10-19C = 4.80320441(42)×10-10esu。
常数换算参考了CODATA在2002年的推荐值，实际运算中没有必要考虑换算系数的不
确定度，请参阅笔者整理的《关于单位制、物理常数和不确定度的资料》。
二、分子构型的输入和输出格式
分子构型可以用笛卡儿坐标、自然坐标或两种坐标的混合来输入。在引擎内部，构
型由笛卡儿坐标(Cartesian)和自然坐标(Z-matrix)各保存一份，做单点计算时利用笛卡
儿坐标，做构型搜索时利用自然坐标。可以采用类似表格形式的写法，坐标必须用统一
单位Bohr(长度)和Radian(角度)。
10: add 1 atom h mass 1838
引擎自动把第一个点置于原点，所以该指令等价于：
15: add 1 atom h mass 1838 cartesian 0.000 0.000 0.000
20: add 2 atom o mass 29150 z-matrix 1 1.820
引擎自动把第二个点置于z轴正方向，所以该指令等价于：
25: add 1 atom h mass 1838 cartesian 0.000 0.000 1.820
30: add 3 atom o mass 29150 z-matrix 2 2.729 1 1.750
引擎自动把第三个点置于yOz平面的y轴正半边，所以该指令等价于：
35: add 3 atom o mass 29150 cartesian 0.000 2.685 2.036
40: add 4 atom o mass 1838 z-matrix 3 1.820 2 1.750 1 2.097
如果用自然坐标表示，键长、键角和二面角还可以使用变量，在上面H2O2的例子中
，用下面的写法可以加快分子构型搜索的速度：
110: var mh value 1838 type mass
120: var mo value 29150 type mass
130: var roh value 1.820 type bond
140: var roo value 2.729 type bond
150: var ahoh value 1.750 type angle
160: var dhooh value 2.097 type dihedral
170: add 1 atom h mass mh
180: add 2 atom o mass mo z-matrix 1 roh
190: add 3 atom o mass mo z-matrix 2 roo 1 ahoh
200: add 4 atom o mass mh z-matrix 3 roh 2 ahoh 1 dhooh
(注：这里标了行号，为后面介绍指令作例子。)
这样，分子构型就只有4个自由度了，比原先的写法减少2个自由度。帮助定义坐标
的点称为虚原子，原子名称是x，质量是0。
一个构型用笛卡儿坐标输出时，采取不同的原点和角度就可以有不同的表示，为了
让构型有可比性，输出格式中应该是唯一的。最好的做法是对四类不同的对称性构型，
给定不同的规则：
(1) 线型，原子全部在Z轴上；
(2) 非陀螺对称型，转动惯量最小的主轴在Z轴，其次是Y轴，最后是X轴；
(3) 陀螺对称型，主轴在Z轴，垂直的二次旋转轴或镜面必须和Y轴重合；
(4) 球形对称型，它包含T点群的对称性，该点群的三个相互垂直的二次旋转轴必须
和三个坐标轴重合。
一个构型也可以用若干个键长、键角和二面角参数表示，但并不是分子中任意两个
原子之间的距离都要给出的，只有成键(距离必须小于某个规定的数值)的两个原子才需
要给出键长，键角和二面角也类似。如果构型包括两个分子，那么可以假设两个分子之
间最近的两个原子是成键的。
三、输入和输出协议
输入协议就是界面向引擎发送指令的协议，协议内容通常由一系列指令集组成(用红
色表示)，输出协议则规定了引擎反馈给界面的信息，以及这些信息具体的含义(用蓝色
表示)。
1. ucqci
UCQCI是Universal Computational Quantum Chemistry Interface的缩写，这条指
令告诉引擎使用计算量子化学的通用引擎协议。
2. id {name <x> | copyright <x> | author <x>}
这是ucqci指令的反馈信息，显示引擎的版本号、版权和作者，例如：
id name Pauling 1.0 // 显示引擎的版本号是Pauling 1.0
id copyright 2005 Fudan Univ // 显示引擎的版权属于Fudan University
id author Aunty Yellow // 显示引擎的作者是Aunty Yellow
4. copyprotection {checking | ok | error}
这也是ucqci指令的反馈信息，只出现在在有版权的引擎中，显示版权检查的信息，
例如：
copyprotection checking // 表示正在做版权检查，此时引擎可能会弹出对话框，要求
用户输入版权信息
copyprotection ok // 表示通过版权检查
copyprotection error // 表示没有通过版权检查，但引擎并不中止运行，只是某些
功能受到限制
3. ucqciok
这是ucqci指令的最后一条反馈信息，表示引擎已经进入用UCQCI协议通讯的状态。
5. show {method | theory} [type <type>]
显示引擎可以提供的搜索方法、理论和基组，有时理论和基组类型很多，可以用<ty
pe>限制显示的类型。
6. method <name> [option <key> ...]]
这是show method的反馈信息，显示引擎支持的搜索方法，例如：
method freq option calcfc // 这是计算频率的方法
method berny option tight calcfc modred symm // 这是搜索稳定构型的方法
method qst2 option tight calcfc modred symm // 这是搜索过渡态的方法
...
7. theory <theory> type <type> [option <key> ...]]
这是show theory的反馈信息，显示引擎支持的理论方法，例如：
theory amber type mmforce
theory am1 type semiemp option firstderiv
theory hf type abinitio option scfbasis firstderiv secondderiv
theory pw91 type dft option dftbasis firstderiv secondderiv
theory b3lyp type mixdft option scfbasis dftbasis firstderiv secondderiv
theory g3x type accurate
...
<theory_type>通常有以下类型：
(1) mmforce：分子力场(经验方法)，有amber、dreiding、uff等；
(2) semiemp：半经验方法，有am1、mindo/3、mndo、pm3、zindo等；
(3) abinitio：从头计算方法，有hf、mp、ci、qci、cc、cas等；
(4) dft：密度泛函方法，有lyp、pw91等一系列方法；
(5) mixdft：混合密度泛函方法，有b3lyp、b3pw86等一系列方法；
(6) accurate：精确方法，有g1、g2、g3、w1等组合方法。
8. new
重新建立构型，引擎启动时并不需要，但每次重新输入构型时需要用这条指令。
9. var <var> value <value> type {mass | coor | bond | angle | dihedral}
定义变量，<var>为变量名，<value>为变量值，单位是原子单位，type后面指定为
变量类型，参见前一节H2O2例子的第110至160行。
10. add <no> atom <atom> [mass <mass>] [<coordinate>]
在构型中加入一个原子，<no>为编号，<atom>为原子名称(用原子序数或元素符号表
示，0或X表示虚原子)，<mass>为原子质量(注意单位是Thomson而不是amu)，<coordinat
e>是坐标，可以是以下形式：
(1) cartesian <x> <y> <z>：直角坐标，<x>、<y>和<z>的单位是Bohr；
(2) z-matrix <no1> <bond> [<no2> <angle> [<no3> <dihedral>]]]：自然坐标，
参见前一节H2O2例子的第170至200行。注：<no>为1时不用z-matrix关键字，<no>为2时
只用<bond>，<no>为3时只用<bond>和<angle>，<bond>的单位是Bohr，<angle>和<dihed
ral>的单位是radian。
要注意的是：编号<no>必须从1开始逐一递增，如果编号在new之后已经用过，那么
这个编号将赋予新的意义，其他和该编号有关的原子坐标可能会改变。另外，用z-matri
x关键字时，<no1>、<no2>和<no3>都必须小于<no>，而且互不相等。
11. del <no>
在构型中删除一个原子，<no>为add过的编号。原子删除后，编号仍然保留，相当于
留下一个虚原子。
12. set {var <var> | atom <no> | bond <no1> <no2> | angle <no1> <no2> <no3>
|
dihedral <no1> <no2> <no3> <no4>} <operation> [<operation> [...]]
这只在构型搜索(不包括频率计算)时起作用，用该指令来冻结或释放某个构型参数
，从而达到减少或增加自由度的目的。使用关键字atom、bond、angle和dihedral时要非
常小心，它们必须是addatom设定过的参数。<operation>可以是以下几种情况：
(1) noincrease | nodecrease：规定参数不能增加或减小，如果两个都选，那么参
数被冻结；
(2) increase | decrease：允许参数增加或减小，用来解除以上两种限定；
(3) minimun | maximum <value>：限定参数的最小值或最大值，但是最小值不能大
于当前值，但是最大值不能小于当前值。<value>可以是数值或infinite(不限定，即用
来解除限定)。
13. isready
询问引擎是否处于“就绪”状态，无论引擎是否正在运算，只要可以接收命令，就
是就绪状态。
14. readyok
这是isready的唯一反馈信息，说明引擎处于“就绪”状态。
15. pop <property> [<option> [...]]
让引擎输出计算结果或分子的其他性质。<property>可以是以下内容：
(1) atomnum，输出分子中原子的个数；
(2) var <var>，输出变量<var>的值；
(3) cartesian <atom>，输出第<atom>个原子的坐标；
(4) {bond <atom1> <atom2> | angle <atom1> <atom2> <atom3> | dihedral
<atom1> <atom2> <atom3> <atom4>}，输出由某几个原子组成的键长、键角和二面角的
数值，原子可以是任意的，只要不超过atomnum的返回值即可。
以上信息会在做构型搜索时改变，因此构型搜索的这些结果必须用pop指令来显示。
(5) energy <type>，输出构型的能量，<type>指计算能量时使用的方法，例如从头
计算的hf、mp2、mp4sdq、qcisd等或者密度泛函的becke3、lyp等。
(6) firstderiv <atom> {x | y | z}，输出能量关于某个原子坐标的一阶导数。
(7) secondderiv <atom1> {x | y | z} <atom2> {x | y | z}，输出能量关于某两
个原子坐标的二阶导数。
(8) charge <type> <atom>，输出第<atom>个原子电荷，<type>指电荷的类型，通
常是mulliken电荷。
(9) dipole | quadrupole | octapole | hexadecapole，输出分子的偶极矩、四极
矩、八极矩、十六极矩等电荷性质。
(10) base <atom> [<options>]，输出加载在第<atom>个原子上的基函数信息，例
如Gaussian函数的系数等等。在不同的基组下可以用不同的<options>选项内容，只输出
需要的信息。
(11) orbital <orbital>，输出轨道信息，轨道可以用1、2、3、4从低到高表示，
也可以用a1g1、a1g2、a1u1、a2g1等带有对称性信息的高低次序表示，还可以用homo2、
homo1、homo、lumo等前线轨道符号表示。
以上信息会在做单点计算时更新，因此单点计算的这些结果也必须用pop指令来显示
，此外对于不同的引擎，pop还可以显示其他信息，这里就不一一列举。
(12) freq <freq> [<options>]，输出频率信息，<freq> 是按照频率由低到高排列
所得的序号，注意n个原子共有3n-3个频率(包括3个转动)。<options>的内容繁多，例如
ir、raman等等。
以上信息会在做频率计算时更新，因此频率计算的这些结果也必须用pop指令来显示
。
16. info <property> <contents>
输出pop指令所要求的信息，有的信息也根据搜索过程的需要直接输出。info 信息
里的<property>对应pop指令中的<property>，这里不一一列举。
17. calc <theory> [<option> [...]]
让引擎进行单点计算(包括一阶导数和二阶导数的分析运算，如果有的话)。单点计
算结果被引擎储存起来。<theory>必须是theory输出信息所指定的理论，<option> 也必
须是，例如：
calc theory b3lyp scfbasis 6-31g* firstderiv
18. calcdone {energy <energy> | error <error_info>}
单点计算结束，error表示运算过程中出现错误(例如自洽场不收敛)，输出错误信息
。energy则说明计算成功，但只输出能量。如果需要其他信息，例如轨道、能量梯度、
电荷密度等其他信息，可以用pop来输出。
19. search <method> [<options>] <calc command>
让引擎进行构型搜索，可以是频率、能量最低构型或过渡态构型。搜索的最终构型
会被引擎存储起来。参数的含义：
<method>：搜索方法，可以是freq(频率)、berny(Berny法搜索能量最低构型)、qst
2(QST2法搜索过渡态构型)等等；
<option>：搜索过程中的选项，如followsymm(保持对称性)、nosymm(取消对称性)
、calcfc (检测频率)、loose(松散收敛)、tight(苛刻收敛)等等；
<calc command>：指定搜索过程中用到的单点计算方法，是一个完整的calc指令，
例如：
search berny tight calc theory hf scfbasis 6-31g* secondderiv
注意，calc指令的选项会影响到构型搜索的效率，例如下面指令的效率就会打很大
的折扣：
search berny tight calc theory hf scfbasis 6-31g* firstderiv
原因是berny搜索需要知道二阶导数，如果只给出一阶导数，它就会用数值方法再做
多个单点算出求得二阶导数。
20. searchdone {failed | energy <energy>}
构型搜索结束，error表示运算过程中出现错误(例如构型不收敛)，输出错误信息。
energy则说明计算成功，但只输出能量。如果需要其他信息，例如键长、键角等，可以
用pop来输出。如果是做频率分析，则必须用pop freq来输出结果。
21. stop [nextstep]
让引擎中断运算，不管是单点计算还是构型搜索。如果加上nextstep选项，则当前
正在进行的独立(不可中断)的任务结束后，再中断运算。否则在continue后，该任务就
必须重新进行这一步运算。
不同的引擎对于“独立任务”有不同的界定，但至少会把每个单点运算都当成独立
任务。
22. continue
让引擎继续运算，
23. quit
让引擎结束工作，程序终止运行。
24. bye
说明引擎已经完成了终止的准备工作(如释放内存等)，以后不会再输出任何信息，
通知界面可以为程序释放句柄了。
--
我的友情测试更新了，欢迎测试!
╔═══════════════════╗
║★★★★★友谊第一 比赛第二★★★★★║
╚═══════════════════╝
※ 来源:·哈工大紫丁香 bbs.hit.edu.cn·[FROM: 天外飞仙]
※ 来源:·哈工大紫丁香 http://bbs.hit.edu.cn·[FROM: 202.118.229.*]