将多帧xyz文件转化成量子化学输入文件的工具：xyz2QC

sobereva

将多帧xyz文件转化成量子化学输入文件的工具：xyz2QC

xyz2QC: A tool to convert multi-frame xyz files to input files of quantum chemistry program

文/Sobereva@北京科音

First release: 2019-Mar-17  Last update: 2021-Aug-26

1 前言

笔者开发的免费的做构象搜索和团簇构型搜索的Molclus程序是如今最为流行的构型和构象搜索程序之一，介绍和下载见官网http://www.keinsci.com/research/molclus.html。Molclus做搜索过程一般流程是先用十分廉价但粗糙的级别（如调用Openbabel跑MMFF94、调用xtb跑GFN-xTB）做结构优化和初步筛选，将其中能量最低的、为数不多的一批体系保留，最后再调用量子化学程序用更准确的方法算能量或做进一步优化。在整个过程中，初筛的耗时很低，在个人计算机上跑也没有任何压力，而最后对为数不多的筛出来的体系进一步做DFT/后HF量子化学计算才是占整个搜索过程的大头，才是真正有必要弄到超算上跑的（对于那些没有自己的像样的服务器的人而言）。然而，在超算上，计算任务一般是以提交方式进行的，而molclus这样自动调用其它程序去运行的方式在超算上不方便使用。为了解决这个矛盾，本文介绍笔者开发的xyz2QC程序。xyz2QC作为Molclus程序（1.7版及之后）的一个子程序发布，在Molclus压缩包里就可以找到。


2 xyz2QC介绍

xyz2QC可以基于Gaussian和ORCA的模板文件，将含有多帧的.xyz文件中指定的一批帧转化成Gaussian和ORCA输入文件。

模板文件的格式要求和molclus的完全一样，Gaussian的必须叫template.gjf，ORCA的必须叫template.inp，运行xyz2QC前必须将之放在当前目录下。

读进xyz2QC的多帧的.xyz可以是genmer、gentor直接产生的，也可以是molclus跑完后产生的isomers.xyz，也可以是isostat筛选、排序、归簇后产生的cluster.xyz。

转化出的Gaussian/ORCA的输入文件，既可以是所有被选择的结构作为多步任务出现在同一个输入文件里的，也可以要求拆分成一个个独立的文件，比如000001.gjf、000002.gjf...文件名是.xyz里的帧号，文件名数字前头可以自行定义前缀。


3 例子

xyz2QC使用极其简单，下面是一个简单例子。

假设我们之前用molclus对某个有机分子在粗糙的半经验方法PM7下做了批量优化，并用isostat做了处理，得到了按照能量排序后的cluster.xyz文件（可以这里下载：http://sobereva.com/attach/472/cluster.xyz）。这里我们想把其中能量最低的三个转化成.gjf文件，之后弄到超算/服务器上算更准确的能量，于是我们准备一个Gaussian模板文件template.gjf，内容为比如

%nproc=4
# M062X/TZVP

Template file

0 1
[GEOMETRY]
<---此处有空行
<---此处有空行

我们将cluster.xyz、template.gjf都放到xyz2QC目录下，然后进入此目录，启动xyz2QC，按照屏幕操作，依次输入：
1  //产生含多步的单一Gaussian输入文件
cluster.xyz  //输入.xyz的实际路径，因为此文件就在当前目录，因此这里不用输绝对路径
1-3  //我们只需要算能量最低的三个结构，由于cluster.xyz里已经按照能量排序，所以这里我们选择xyz文件里的第1~3帧

此时程序提示Gaussian.gjf已经产生在了当前目录，按回车退出即可。这个文件的内容是这样的：
%nproc=4
# M062X/TZVP

Template file

0 1
O          -2.30907426      1.90383755      0.11387296
...略
H          -3.80099226     -1.36391720      0.45491869

--link1--
%nproc=4
# M062X/TZVP

Template file

0 1
O          -2.47455224      1.86708010     -0.16509991
...略
H          -3.70656780     -1.42641602      0.64335011

--link1--
%nproc=4
# M062X/TZVP

Template file

0 1
O          -2.34444893      1.89718026      0.06136974
...略
H          -3.75360804     -1.39079507      0.53898348

有一定Gaussian基础的人都知道，在Gaussian的输入文件里可以用--link1--来将多步任务进行分隔，因此直接执行Gaussian.gjf，这三个结构的能量就都会挨个算出来了。之后我们可以用比如Ultraedit或者Linux下的grep命令直接把带有SCF Done字样的行都提取出来，就直接得到这三个结构的能量了，非常方便！

注意--link1--前头的空行数目必须对，否则多步任务将运行不成功，上面产生的Gaussian.gjf这样是没问题的。--link1--前面空行数目取决于template.gjf末尾的空行数。

xyz2QC还允许将.xyz里的指定帧拆分成产生单独的.gjf文件，这有何意义？一般来说，像上面这样把多个结构弄到一个.gjf里运行起来最方便，但有一个缺点是如果其中一个任务失败了（如SCF不收敛、几何优化不收敛），那么整个任务就会停了，后面的步骤也不算了。如果拆成单独的.gjf文件，那就没这个问题了。

Gaussian模板文件中可以使用[FILENAME]占位，xyz2QC产生输入文件时它将被替换为"前缀 帧号"，因此可以在模板文件里写比如%chk=/sob/[FILENAME].chk从而使得计算完毕后chk文件以恰当的名字留存在/sob目录下。这对于通过Multiwfn程序批量做波函数分析尤其有用。当有多帧的xyz文件，想对里面每个结构都用Multiwfn做波函数分析时，可以拆分完了后批量执行，得到对应的chk文件，然后再批量转换成fch，最后通过批处理脚本调用Multiwfn批量分析。怎么靠脚本实现批量转换和执行见《使用Gaussian时的几个实用脚本和命令》（http://sobereva.com/258）。

类似地，xyz2QC还可以产生ORCA的输入文件，也是可以作为多步任务构成单一文件，或拆分成不同的文件。对于多步任务，强烈建议在模板文件template.inp里加入noautostart关键词，要不然下一个任务的初猜会自动用前一个任务收敛的波函数，然而你选定的那些结构的几何差异可能不小，因此这样得到的初猜可能会很糟糕。noautostart代表不自动从与当前文件同名的.gbw文件中读取初猜波函数，就可以确保没这个问题。

happyknighthawk

Sob老师，斗胆提一个请求，能否增加一个功能，就是在生成独立的gjf或者inp文件的时候，能够把这些文件统一命名成“前缀（用户指定）+编号.gjf”这种模式。当然，如果您没时间或没兴趣做这件事，就当我没说，非常感谢您

sobereva

happyknighthawk 发表于 2019-3-26 16:25
Sob老师，斗胆提一个请求，能否增加一个功能，就是在生成独立的gjf或者inp文件的时候，能够把这些文件统一 ...

我感觉一般用户用不着这点，而且还会额外增加一个操作步骤
如果你需要这个功能，我可以给你源码自行添加，也就几行语句的事

happyknighthawk

sobereva 发表于 2019-3-27 04:33
我感觉一般用户用不着这点，而且还会额外增加一个操作步骤
如果你需要这个功能，我可以给你源码自行添加 ...

谢谢Sob老师回复，我的邮箱是hukun@mail.kib.ac.cn，麻烦您了

chemqzx

采用sob老师提供的GauIRC2xyz.exe将IRC中每个点的坐标转为.xyz文件后，通过本文提供的方法将xyz转为QC时，遇到所有转换的结构与template.gjf的结构相同，请问大家有遇到这样的问题吗？请问是如何解决的呢？谢谢！麻烦啦！

sobereva

chemqzx 发表于 2019-5-6 18:47
采用sob老师提供的GauIRC2xyz.exe将IRC中每个点的坐标转为.xyz文件后，通过本文提供的方法将xyz转为QC时， ...

template.gjf里根本就没有结构信息，结构信息部分被[GEOMETRY]所代替

莫落

社长能否开发个小软件，将xyz轨迹文件转换成晶体文件（如cif文件），每次CP2K结果出来，手动建立晶胞可麻烦了

sobereva

莫落 发表于 2019-5-7 11:18
社长能否开发个小软件，将xyz轨迹文件转换成晶体文件（如cif文件），每次CP2K结果出来，手动建立晶胞可麻烦 ...

暂时不在考虑范畴之内，等我什么时候研究涉及到相关问题时可能才会写

卡开发发

莫落 发表于 2019-5-7 11:18
社长能否开发个小软件，将xyz轨迹文件转换成晶体文件（如cif文件），每次CP2K结果出来，手动建立晶胞可麻烦 ...

我不知道你的晶格是否固定尺寸，xyz文件我不确定是否会丢失晶格信息。对于vesta的情况，cif只需要alpha、beta、gamma、a、b、c、原子的符号和坐标即可，我不知道cp2k读取的cif只要这些信息就能读取结构。如不涉及对称性问题，实现起来不是很难，我可以用python帮你试试看。

莫落

卡开发发 发表于 2019-5-8 12:04
我不知道你的晶格是否固定尺寸，xyz文件我不确定是否会丢失晶格信息。对于vesta的情况，cif只需要alpha、 ...

CP2K可以直接读取cif格式的晶体文件，晶格参数在inp文件里。CP2K优化后的结构输出到.xyz文件里，而晶格参数输出到.out文件中，所以要结合两个文件提供的信息重新建立晶胞，比较麻烦。

卡开发发

莫落 发表于 2019-5-8 15:25
CP2K可以直接读取cif格式的晶体文件，晶格参数在inp文件里。CP2K优化后的结构输出到.xyz文件里，而晶格参 ...

能够提取到晶格参数的话这个也不是问题，这个可以讨论下。

卡开发发

莫落 发表于 2019-5-8 15:25
CP2K可以直接读取cif格式的晶体文件，晶格参数在inp文件里。CP2K优化后的结构输出到.xyz文件里，而晶格参 ...

这里以xtl为例子，因为同样能够被M$识别。python实现CP2K做cif或xtl格式转换的要点在于：
1、原子坐标由直角坐标转换为分数坐标，转换过程需要倒格矢量，参考Callaway的《Quantum Theory of the Solid State》的$1.2

2、通过晶格矢量计算计算晶格参数，涉及到简单的矢量模长和夹角计算。
附件中的CrysFormat.py能够实现上述过程。

读取out、xyz和写入xtl比较简单，核心技巧在于python正则表达式的使用，附件中的CP2K_io.py能够实现。CP2K的问题在于晶格矢量保留的小数太少了，因此转换过程会有数值误差。

最后，主程序读取out和xyz的文件名和指定xtl的文件名，调用前面两个py程序即可，附件中的CP2K_xtl是一个例子，根据需求原则上也能写出转换成cif的格式。

hzamber

happyknighthawk 发表于 2019-3-26 16:25
Sob老师，斗胆提一个请求，能否增加一个功能，就是在生成独立的gjf或者inp文件的时候，能够把这些文件统一 ...

批处理脚本，满足你的要求。
http://bbs.keinsci.com/thread-13238-1-1.html

highlight

莫落 发表于 2019-5-8 15:25
CP2K可以直接读取cif格式的晶体文件，晶格参数在inp文件里。CP2K优化后的结构输出到.xyz文件里，而晶格参 ...

你有没有尝试通过 MOTION / PRINT / CELL 把晶格参数单独输出出来？
输出的文件结合 xyz一起读入，应该不再需要 cif 了吧。

sobereva

happyknighthawk 发表于 2019-3-26 16:25
Sob老师，斗胆提一个请求，能否增加一个功能，就是在生成独立的gjf或者inp文件的时候，能够把这些文件统一 ...

刚更新了Molclus 1.8.2版，这次的xyz2QC组件产生独立的Gaussian和ORCA输入文件时已经允许用户输入前缀了。