Gaussian的Link、IOp与非标准计算路径
Gaussian的Link、IOp与非标准计算路径
Link, IOp and non-standard calculation route of Gaussian
文/Sobereva @北京科音 2010-Mar-9
Link, IOp and non-standard calculation route of Gaussian
文/Sobereva @北京科音 2010-Mar-9
1. Gaussian程序的基本结构
高斯是由很多功能不同的独立可执行程序组成的，它们被称为Link，在它们之间通过临时文件来传递数据，最主要的是rwf文件，每个Link的功能在高斯手册里都有简单介绍。Link间的调用顺序实际上是随意的，但是只有某些固定顺序的调用（还包括循环、分支、跳转方式的调用）才有意义，才能完成特定的任务。每个Link对应的可执行文件名的通式为lyyxx.exe，例如L314对应的可执行程序l314.exe就是yy=3、xx=14。yy代表的是这个Link所属的overlay，也称为层，功能相近的子程序被归在一层，如overlay 6主要是分析波函数，所以包括了L601（布居数和相关的分析）、L607（NBO分析）、L609（AIM分析）等子程序。同一层的Link也经常连在一起被调用，如L301、L302、L303、L311、L314它们总是连在一起调用来完成基组的设定和电子积分。L0、L1、L9999相对其它Link较为特殊，它们并不用来执行实际的计算，而且不论何种任务，它们都一定会执行。
L0：在程序最开始执行，用来初始化运行环境。与L0对应的程序并不叫L0.exe，若在高斯03中，dos版本它就对应g03.exe，windows GUI版本就是g03W.exe，unix版本就对应g03。L0会创建L1进程，然后在全部Link调用结束前会一直保持睡眠状态。L0还会把输入文件复制成一份“内部输入文件”到用户设的scratch目录中，其名字在windows版本下就是gxx.inp，在unix版本下为“进程id号.inp”，任务正常结束后会被自动删掉，其内容相对于输入文件没有了注释行（开头为!的行），并且用@来引用的外部文件会在其中被展开。
L1：被Link 0所调用。用来初始化临时文件、解析输入文件中route section段落的关键词、生成调用以后Link所需的命令行指令。L1和L0共同组成overlay 0。
L9999：最后一个被调用的Link，进行任务的扫尾工作。它确保最重要的信息已写入chk，输出第三方程序所需的文件（如.wfn波函数文件），生成档案条目（是指在格言之前输出的那段计算数据的紧凑描述），输出古怪的格言，最后结束整个任务。
2. 内部选项(IOp,Internal Option)
每个Link都有其IOp，各个选项的选项值决定了这个Link如何运行，每个IOp的解释可以在高斯官方网站查到，遇见个别查不到的，需用《察看Gaussian全部IOP的方法》（http://sobereva.com/30）文中的方法。在route section段落所指定的计算任务和方法被L1所解析后，哪些Link会被调用，并且以什么参数来调用就确定了。例如# HF/STO-3G会被解析成下面这样，称为非标准路径（不同高斯版本解析出的结果可能有些不同）：
1/38=1/1;
2/17=6,18=5,40=1/2;
3/6=3,11=9,16=1,25=1,30=1/1,2,3;
4//1;
5/5=2,32=1,38=5/2;
6/7=2,8=2,9=2,10=2,28=1/1;
99/5=1,9=1/99;
上面每行内容格式是：yy/选项=数值,选项=数值.../xx1,xx2,xx3...;。例如第5行，yy=5，xx=2说明要调用L502，并且令其选项5的值为2、选项32的值为1、选项38的值为5。再例如第3行，说明要调用L301、L302、L303，并且将选项6的值为3、选项11的值为9...选项30的值为1这样的设定传递给它们，属于同一层且在同一行的Link会共享IOp设定。
解析出的内容并不是把所有选项和赋值都列出来，选项的值如果是默认的就不会被列出，一般各个IOp默认的值都是0。例如第四行，说明以各个选项都用默认值的方式来调用L401。
可以用IOp关键字修改各层选项的值，格式为IOp(层/选项=值,层/选项=值...)，比如route section多填上IOp(3/40=2,3/11=1,4/9=2)，则第三行成了3/6=3,11=1,16=1,25=1,30=1,40=2/1,2,3;，第四行成了4/9=2/1;。
给某个Link传递它本身并没有的选项设定等于没有传递。例如freq任务被解析出的内容会有7/8=1,10=1,25=1/1,2,3,16;这样一行，而L701/702/703并没有选项8，所以虽然写了，但实际上等价于只给L716传递了选项8=1这个设定。同一层的不同Link往往同时具有相同的选项，选项的设定传递给它们时所代表的含义有时相同（必定这几个Link功能有相关性），但也有时不同（这种情况下这些Link不会有机会在一起使用）。例如IOp(1/13)控制L103/113/114（用于几何优化）及L115（计算IRC）用什么方法更新Hessian矩阵，也控制L121（ADMP计算）的多时间步长参数。对这两种情况1/13的功能是不相关的，却只能给这个选项设一个值，看似会有矛盾，但由于ADMP不可能与优化或IRC在一起执行，所以不必担心这个问题。
想查看route section被解析的结果，可以进入高斯自带的testrt.exe程序，然后输入比如HF/STO-3G，按两次回车就可得到上述结果。或者在windows版高斯的Existing File Job Edit窗口中点对勾按钮。也可以在任务的输出文件的开始部分看到（使用#T精简输出时除外）。
3. Link的调用顺序
来看一个复杂点的例子，# HF/STO-3G opt的非标准路径如下：1/18=20,38=1/1,3;
2/9=110,17=6,18=5,40=1/2;
3/6=3,11=1,16=1,25=1,30=1/1,2,3;
4/7=1/1;
5/5=2,38=5/2;
6/7=2,8=2,9=2,10=2,28=1/1;
7//1,2,3,16;
1/18=20/3(3);
2/9=110/2;
6/7=2,8=2,9=2,10=2,19=2,28=1/1;
99//99;
2/9=110/2; 重新定位坐标，计算对称性，检查变量
3/6=3,11=1,16=1,25=1,30=1/1,2,3; 产生基组信息，计算积分
4/5=5,7=1,16=3/1; MO初猜
5/5=2,38=5/2; SCF迭代
7//1,2,3,16; 计算导数，处理优化和频率的信息
1/18=20/3(-5); Berny优化（计算新坐标，判断是否位移和力已收敛）
2/9=110/2;
6/7=2,8=2,9=2,10=2,19=2,28=1/1;
99/9=1/99;
Link会按照非标准路径的书写顺序从上往下执行，即L101->L103->L202->L301->L302...，利用#P输出的文件内的Enter...和Leave Link...提示就能验证。在某些行末（但在分号之前）会看到括号标识，它表明完成当前层后如何跳转，若值为n，当n>0表明要往后跳转n+1次，当n<0表明要往前跳转n次，默认是0，即按顺序执行下一行而不跳转。比如第8行1/18=20/3(3);说明这行执行结束后（当前位置仍在这行）向后跳4次，执行2/9=110/2;。倒数第四行1/18=20/3(-5);说明运行完L103之后往前跳5次，即接下来执行2/9=110/2;，这样就形成了一个循环，反复计算新的电子积分、计算导数并移动坐标，构成了几何优化的最主要过程。
跳转指令一般都会被执行，但对于某些Link是否跳转也取决于其执行结果，否则上述循环将无限进行下去。第一次调用L103后遇到的跳转命令(3)只要当L103正常结束就会执行，所以一般都会跳过L9999，而不会因此直接结束任务。进入循环后，当循环内的L103发现位移和受力都已经收敛，就不再执行(-5)所示的跳转命令，而是接着向下执行L202输出当前坐标变量，然后运行L601输出布居、多极矩等信息，最终到L9999结束任务。
4. 自定义计算路径
上面的两个非标准路径的例子虽然名字含有“非标准”，但实际上Link的调用顺序是标准的，是完成程序内设功能的默认流程。用户也可以自己定制计算路径，成为真正意义的“非标准”路径。方法是使用nonstd关键字，然后在后面写上非标准路径。也可以结合#T、#P、oldconstant等设定，因为它们并不属于非标准路径涉及的设定，与nonstd没有冲突。例如这是一个以HF/STO-3计算水分子的任务所对应的自定义计算路径输入文件：%chk=C:\ltwd\h2o.chk
#P nonstd
1/38=1/1;
2/17=6,18=5,40=1/2;
3/11=9,16=1,25=1,30=1/1,2,3;
4//1;
5/5=2,32=1,38=5/2;
6/7=2,8=2,9=2,10=2,28=1/1;
99/5=1,9=1/99;
Title Card Required
0 1
O 0.00000000 0.00000000 -0.11094313
H -0.00000000 -0.78383672 0.44331313
H -0.00000000 0.78383672 0.44331313
我们可以以此为基础修改。例如只想获得距离矩阵及分子点群信息，不想干别的，只需要执行L202就可以了，于是非标准路径只保留下面这三行：
1/38=1/1; 读取标题和分子说明部分
2/17=6,18=5,40=1/2; 重新定位坐标，计算对称性，检查变量
99/5=1,9=1/99;
5. 自定义计算路径应用实例
自定义计算路径有一定实用价值。假设我们想获取几何优化中每步的Mulliken电荷，然而高斯默认的是只在内循环前后各调用一次L601，而不会每步都输出一次电荷。想达成这个目的可以写脚本提取优化任务的输出文件中的结构，然后生成一堆单点输入文件再批量执行，显然这比较麻烦。利用自定义非标准路径则容易得多，只需要把6/7=2,8=2,9=2,10=2,28=1/1;这行挪到内循环里面即可，上面的例子即改为：前略...
7//1,2,3,16;
6/7=2,8=2,9=2,10=2,28=1/1; <---添加此行
1/18=20/3(-6); <---由于前面多加了一行，为了跳回原处，所以(-5)也改成了(-6)。
2/9=110/2;
后略...
当然也可以把此行写在7//1,2,3,16;的前面，因为它们并不影响密度矩阵，所以结果一样。但不应写在5/5=2,38=5/2;之前，因为当前结构收敛的波函数此时尚未生成。也不能写成：
...
1/18=20/3;
6/7=2,8=2,9=2,10=2,28=1/1(-6);
...
这样的话会陷入无限循环，因为L601不能判断什么时候不做跳转。
再举一个例子，在scan任务中获得每一步的AIM键临界点。先用testrt获知scan关键字对应的是：
...前略
2/29=3/2;
3/11=9,16=1,25=1,30=1/1,2,3;
4/5=5,16=3/1;
5/5=2,38=5/2;
1/60=1/8(-4);
99/9=1/99;
再将aim关键字输入testrt，可获知在非标准路径中与aim关键字有关的是6/7=2,8=2,9=2,10=2,28=1,35=4/1,9;，由于不打算运行L601，而且只有选项35对L609有意义，所以可以只写6/35=4/9;。故改成下面这样，并写在# nonstd之后即可。
...前略
2/29=3/2;
3/11=9,16=1,25=1,30=1/1,2,3;
4/5=5,16=3/1;
5/5=2,38=5/2;
6/35=4/9;
1/60=1/8(-5);
99/9=1/99;
其中L108用来做势能面扫描，当步数达到所设值时会忽略跳转命令往下执行。
同理，利用上述方法结合IRC计算，还可以一次得到IRC路径各点的Mayer键级(6/80=1)，用脚本提取出来后就能绘制反应路径上的键级变化曲线。
6. 与计算路径相关的几个选项
有两个Link 0命令（写在%行的设定）值得一提，可以结合非标准路径使用。%KJob Lxxx M 说明执行Lxxx M次之后停止任务，此选项方便直接地人为控制Link的调用。比如%kjob L502 3说明执行L502三次后立刻中断任务，最后不会经过L9999。
%Subst Lxxx dir 从dir目录下取得Lxxx.exe的执行文件，代替默认路径下的同名Link执行文件。如果某个子程序自行修改过，有多个版本放在不同文件夹，可以用这个功能来指定调用哪个。如果需要同时指定多个Link的路径，就写多次这个条目依次指定。
还有几个route section关键字也与执行路径相关，虽然也可以直接用nonstd的方法实现相同功能，但直接用这些关键字有时更方便，免得把非标准路径全写出来。
Skip=M 跳过前M层。比如第一个例子若写成# HF/STO-3G skip=5，则非标准路径就只剩下了6/7=2,8=2,9=2,10=2,28=1/1; 和 99/5=1,9=1/99;。
Skip=Ovn 跳过overlay=n前面的层。如使用Skip=Ov4则第一个例子从4//1开始执行。
ExtraLinks=yyxx 每次yy层执行完毕后再执行一次Lyyxx。例如第三节的例子，如果写成# HF/STO-3G opt ExtraLinks=607，则执行完L602就会接着执行L607(NBO分析)，在非标准路径上也能看到发生了改变，三处6/.../1;都变成了6/.../1,7。
ExtraOverlays 用了这个关键词后，空一行后写上执行内容（并且与标题行之间留一空行），则这行会被插入到overlay 9之前。例如第一个例子的输入文件写成：
# SP ExtraOverlays
2/9=1/2;
Title Card Required
...后略
从非标准路径上就会看到，2/9=1/2;被插入到了99/5=1,9=1/99;的前面，所以任务临结束前会调用一次L202。