计算化学公社
标题: 将多帧xyz文件转化成量子化学输入文件的工具：xyz2QC [打印本页]
作者Author: sobereva 时间: 2019-3-17 12:16
标题: 将多帧xyz文件转化成量子化学输入文件的工具：xyz2QC
将多帧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文件中读取初猜波函数，就可以确保没这个问题。
作者Author: happyknighthawk 时间: 2019-3-26 16:25
Sob老师，斗胆提一个请求，能否增加一个功能，就是在生成独立的gjf或者inp文件的时候，能够把这些文件统一命名成“前缀（用户指定）+编号.gjf”这种模式。当然，如果您没时间或没兴趣做这件事，就当我没说，非常感谢您
作者Author: sobereva 时间: 2019-3-27 04:33
我感觉一般用户用不着这点，而且还会额外增加一个操作步骤
如果你需要这个功能，我可以给你源码自行添加，也就几行语句的事
作者Author: happyknighthawk 时间: 2019-3-27 09:55
谢谢Sob老师回复，我的邮箱是hukun@mail.kib.ac.cn，麻烦您了
作者Author: chemqzx 时间: 2019-5-6 18:47
采用sob老师提供的GauIRC2xyz.exe将IRC中每个点的坐标转为.xyz文件后，通过本文提供的方法将xyz转为QC时，遇到所有转换的结构与template.gjf的结构相同，请问大家有遇到这样的问题吗？请问是如何解决的呢？谢谢！麻烦啦！
作者Author: sobereva 时间: 2019-5-7 08:12
template.gjf里根本就没有结构信息，结构信息部分被[GEOMETRY]所代替
作者Author: 莫落 时间: 2019-5-7 11:18
社长能否开发个小软件，将xyz轨迹文件转换成晶体文件（如cif文件），每次CP2K结果出来，手动建立晶胞可麻烦了
作者Author: sobereva 时间: 2019-5-8 09:48
暂时不在考虑范畴之内，等我什么时候研究涉及到相关问题时可能才会写
作者Author: 卡开发发 时间: 2019-5-8 12:04
我不知道你的晶格是否固定尺寸，xyz文件我不确定是否会丢失晶格信息。对于vesta的情况，cif只需要alpha、beta、gamma、a、b、c、原子的符号和坐标即可，我不知道cp2k读取的cif只要这些信息就能读取结构。如不涉及对称性问题，实现起来不是很难，我可以用python帮你试试看。
作者Author: 莫落 时间: 2019-5-8 15:25
CP2K可以直接读取cif格式的晶体文件，晶格参数在inp文件里。CP2K优化后的结构输出到.xyz文件里，而晶格参数输出到.out文件中，所以要结合两个文件提供的信息重新建立晶胞，比较麻烦。
作者Author: 卡开发发 时间: 2019-5-8 16:21
能够提取到晶格参数的话这个也不是问题，这个可以讨论下。
作者Author: 卡开发发 时间: 2019-5-13 14:06
本帖最后由 卡开发发 于 2019-5-16 11:53 编辑
这里以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的格式。
[attach]19033[/attach]
作者Author: hzamber 时间: 2019-5-14 16:59
批处理脚本，满足你的要求。
http://bbs.keinsci.com/thread-13238-1-1.html
作者Author: highlight 时间: 2019-5-14 17:39
你有没有尝试通过 MOTION / PRINT / CELL 把晶格参数单独输出出来？
输出的文件结合 xyz一起读入，应该不再需要 cif 了吧。
作者Author: sobereva 时间: 2019-5-15 20:42
刚更新了Molclus 1.8.2版，这次的xyz2QC组件产生独立的Gaussian和ORCA输入文件时已经允许用户输入前缀了。
作者Author: happyknighthawk 时间: 2019-5-15 22:39
谢谢，谢谢！
作者Author: happyknighthawk 时间: 2019-5-15 22:40
非常感谢Sob老师，这就去下载新版本！
作者Author: chemqzx 时间: 2019-6-19 15:43
感谢sob老师的回复，谢谢您！
作者Author: litangcheng 时间: 2019-10-12 16:08
想请教下楼主，cluster.xyz中有三个团簇结构，为什么用xyz2QC只生成了一个Gaussian.gif，而且Gaussian.gif中也只有一个团簇结构。
作者Author: litangcheng 时间: 2019-10-12 16:25
我知道了，看了下Viewfile，明白了
作者Author: wuzhiyi 时间: 2019-11-5 17:20
想请问一下sob老师，在molclus套件里有没有xyz2QC的反向操作QC2xyz啊？
就是用处理xyz2QC运行完的结果 然后把他变成可以用isostat处理的isomer.xyz文件。
作者Author: sobereva 时间: 2019-11-11 11:29
想法有意思，计划加入下一个版本
作者Author: wuzhiyi 时间: 2019-11-11 18:28
谢谢sob老师
作者Author: wuzhiyi 时间: 2019-11-19 17:17
发现如果第一行大于80个字符会被强行截断到80个字符
比如我有一个orca的模板
第一行是! DLPNO-CCSD(T) tightPNO aug-cc-pVTZ aug-cc-pVTZ/C tightSCF noautostart miniprint nopop
xyz2QC生产的结果就是前80个字符
! DLPNO-CCSD(T) tightPNO aug-cc-pVTZ aug-cc-pVTZ/C tightSCF noautostart miniprin
作者Author: sobereva 时间: 2019-11-20 01:26
我用的2019-Sep-23版的xyz2QC没发现这个问题
模板文件和产生的文件都附上了。当xyz2QC里选择的是ORCA而不是Gaussian的情况，读取和写入的字符上限都是200个字符
作者Author: wuzhiyi 时间: 2019-11-20 04:52
谢谢sob老师 发现是1.88版的问题 用最新版就解决了
作者Author: Evanwill 时间: 2020-6-4 14:13
提个小建议，单独QC文件的命名是按帧号来的，这样的话只有很少构象的情况，比如5个构象，文件名也可能会是00005.gjf。有没有方法，比如输出100个文件，文件名自动从001到100，输出前1000个话，就是0001到1000，等等
作者Author: sobereva 时间: 2020-6-6 10:53
我觉得统一格式更好点，之后写额外的处理的程序的时候更方便
作者Author: snljty 时间: 2020-6-6 16:11
卢老师，请问您的QC2xyz程序开发还有打算么？感觉从molclus中提取一部分代码改一下就好了，但是应该是个很实用的工具。谢谢您！
作者Author: sobereva 时间: 2020-6-6 18:43
等2.0的时候一并提供。最近太忙，有更重要的事需要做
作者Author: snljty 时间: 2020-6-6 18:52
谢谢老师！
作者Author: Evanwill 时间: 2020-6-10 13:45
再提个建议，不知道是不是只有我，我觉得，生成的gjf等文件最好能够存储在xyz文件所在的文件夹，而不是默认的molclus文件夹。或者，能够在任何文件夹，通过xyz2QC或者isostat等命令直接调用molclus的组件
作者Author: sobereva 时间: 2020-6-11 23:53
在其它文件夹调用isostat、xyz2QC，把molclus加入PATH环境变量就完了
作者Author: lanthanum 时间: 2021-8-26 01:48
本帖最后由 lanthanum 于 2021-8-26 01:59 编辑
斗胆提个请求，能否改进一下，使功能1和2识别类似如下样式的多步模板。
%chk=[FILENAME]_step01.chk
#p B3LYP/6-311G** em=GD3BJ opt freq scrf(read,SMD,solvent=n,n-DiMethylAcetamide)
Template file
0 1
[GEOMETRY]
eps=19.1
epsinf=2.090916
--link1--
%oldchk=[FILENAME]_step01.chk
%chk=[FILENAME]_step02.chk
#p M052X/6-31G* scrf=check geom=allcheck
--link1--
%oldchk=[FILENAME]_step01.chk
%chk=[FILENAME]_step03.chk
#p M052X/6-31G* geom=allcheck
--link1--
%oldchk=[FILENAME]_step01.chk
%chk=[FILENAME]_step04.chk
#p B2PLYPD3/def2TZVP geom=allcheck
其中的 [FILENAME] 将被替换成 "前缀+帧号" ，最后功能1产生一个4N步的文件，或者功能2产生N个4步的文件。
谢谢社长。
作者Author: sobereva 时间: 2021-8-26 22:38
刚才更新了官网上的molclus 1.9.9.5，并且在此文里添加了以下内容
Gaussian模板文件中可以使用[FILENAME]占位，xyz2QC产生输入文件时它将被替换为"前缀 帧号"，因此可以在模板文件里写比如%chk=/sob/[FILENAME].chk从而使得计算完毕后chk文件以恰当的名字留存在/sob目录下。这对于通过Multiwfn程序批量做波函数分析尤其有用。当有多帧的xyz文件，想对里面每个结构都用Multiwfn做波函数分析时，可以拆分完了后批量执行，得到对应的chk文件，然后再批量转换成fch，最后通过批处理脚本调用Multiwfn批量分析。怎么靠脚本实现批量转换和执行见《使用Gaussian时的几个实用脚本和命令》（http://sobereva.com/258）。
作者Author: lanthanum 时间: 2021-8-26 23:09
本帖最后由 lanthanum 于 2021-8-26 23:11 编辑
太棒了，谢谢社长。
有个小bug。
如果模板和以前一样，末尾有两个空行，那功能1生成的多步文件中，4N步与4N+1步之间是两个空行。
如果模板末尾留一个空行，那功能1和2生成的全部文件的末尾都只有一个空行。
作者Author: zjxitcc 时间: 2021-8-27 00:18
QC2xyz你想具体指啥呢？从各种量化程序的输出文件中提取出坐标？
作者Author: sobereva 时间: 2021-8-27 06:34
--link1--前头空行有多少无所谓，不是bug
作者Author: lanthanum 时间: 2021-8-27 11:27
本帖最后由 lanthanum 于 2021-8-27 11:34 编辑
谢谢社长。
试了，--link1--前头空行有多少确实无所谓，这不是bug。
另外有个小bug，功能2，在第一步会产生两个%chk=字段，一个是替换产生的%chk=C:\work\000001_step01.chk，一个是另添上的%chk=000001.chk。
作者Author: snljty 时间: 2021-8-27 11:50
其实就是molclus里早就实现了的那些功能，按照settings.ini的设置读几何结构（优化的最后一帧，频率或单点任务的机构），以及能量（电子能量（ORCA的电子能量，Gaussian的SCF能量，MP2级别能量等），频率任务的热力学能量）等写在xyz文件的注释行。只不过现在的版本是每次读计算临时文件的gau.out，ORCA.out里面的数据，如果能指定读一批别的地方计算好的ORCA0001.out,ORCA.0002.out并写到一个xyz文件里就好了。感觉这个只要做十几行代码简单的修改就可以在现有molclus的基础上实现的，所以斗胆问了一下卢老师。
作者Author: lanthanum 时间: 2021-8-27 14:56
这有个变通的办法，使用GV的分子组，把全部一批文件读到同一个分子组里（打开文件时，target选择 single new molecule group for all files），在分子组上，点鼠标右键菜单，results，molecule group table，可以看到能量，点击表头，可以排序，右下角有copy功能，左上角可以选择显示3D结构。
作者Author: snljty 时间: 2021-8-27 15:20
GaussView提取的能量经常不对，比如双杂化泛函显示的SCF部分能量。比如卢老师这里提的。http://bbs.keinsci.com/thread-13450-1-1.html
作者Author: sobereva 时间: 2021-8-28 11:22
你把你的模板文件上传我看看，并且告诉我你在xyz2QC里的完整操作。之前测试过，应该不会自动添上%chk=000001.chk
作者Author: lanthanum 时间: 2021-8-28 13:54
xyz2qc
2
cluster.xyz
回车
回车
模板和第一个结果文件，见附件。cluster.xyz是本贴1楼的文件。
作者Author: sobereva 时间: 2021-8-28 14:14
我重新更新了，现在没问题了
作者Author: lanthanum 时间: 2021-8-28 19:29
没问题了，谢谢社长。
作者Author: ananbaby2309 时间: 2022-4-25 08:32
请问脚本在哪里下载呀？（我可能瞎了，一直找不到……
作者Author: sobereva 时间: 2022-4-25 09:48
一个字一个字读本文第一段
作者Author: ReviewReview 时间: 2022-6-8 10:03
本帖最后由 ReviewReview 于 2022-6-8 10:16 编辑
"xyz2QC还允许将.xyz里的指定帧拆分成产生单独的.gjf文件，这有何意义？一般来说，像上面这样把多个结构弄到一个.gjf里运行起来最方便，但有一个缺点是如果其中一个任务失败了（如SCF不收敛、几何优化不收敛），那么整个任务就会停了，后面的步骤也不算了。如果拆成单独的.gjf文件，那就没这个问题了。"，Sober老师，这个具体怎么做呢？
作者Author: ReviewReview 时间: 2022-6-8 10:06
本帖最后由 ReviewReview 于 2022-6-8 11:17 编辑
"类似地，xyz2QC还可以产生ORCA的输入文件，也是可以作为多步任务构成单一文件，或拆分成不同的文件。对于多步任务，强烈建议在模板文件template.inp里加入noautostart关键词，要不然下一个任务的初猜会自动用前一个任务收敛的波函数，然而你选定的那些结构的几何差异可能不小，因此这样得到的初猜可能会很糟糕。noautostart代表不自动从与当前文件同名的.gbw文件中读取初猜波函数，就可以确保没这个问题。"，Sober老师，假如我选择的是Gaussian软件而非ORCA软件来使用多步任务做Opt+Freq和Energy，我也要在template.gjf文件(Gaussian软件对应的模板文件)中加入noautostart关键词吗？
作者Author: ReviewReview 时间: 2022-6-8 10:26
对不起Sober老师，是我瞎准备把不用的眼睛捐掉哈哈，刚点开xyz2QC.exe看到窗口的选项了，我把"1 //产生含多步的单一Gaussian输入文件"里的1改成2即可
。
作者Author: WOOOOWOOOO 时间: 2023-11-23 22:41
请问我用xyz2QC生成的输入文件计算后，怎么弄一个能用isostat.exe处理的xyz文件？自己用手复制来复制去搞不好
，而且我还是opt任务，能量和坐标不好找。。。。。
作者Author: sobereva 时间: 2023-11-24 02:00
目前没办法直接弄
作者Author: WOOOOWOOOO 时间: 2023-11-24 10:05
那格式的关键点是什么呀？要整理成完全一致吗，空格之类的都要一致吗，从Gaussian计算后的输出文件里摘出来的坐标和xyz文件坐标的空格这些都不一样。
作者Author: wzkchem5 时间: 2023-11-24 19:29
xyz文件的空格数，不是xyz文件格式的一部分。所以但凡是一个写得足够好的程序，读xyz文件的时候都不会对xyz文件的空格数有要求。
如果是初学编程的人写的读xyz文件的程序，那另当别论
作者Author: sobereva 时间: 2023-11-25 01:13
谈谈记录化学体系结构的xyz文件
http://sobereva.com/477
坐标部分是自由格式
注释部分尽量按照molclus产生的xyz文件的格式写，否则不保证isostat能认出里面的能量等信息
作者Author: 0325-ddd 时间: 2024-2-3 15:15
您好！斗胆请教一下您，算无机物的时候，cif文件转换为mol再转化为gjf这样计算的结果比较准确嘛？
作者Author: sobereva 时间: 2024-2-3 15:25
这种问题没意义，你只需要关注文件里记录了什么信息，如果记录的信息是等价的，并且程序在读取时都会读到，显然跟格式没关系。转化成mol完全莫名其妙，本来mol格式也不是专门记录周期性体系用的
作者Author: 0325-ddd 时间: 2024-2-3 15:30
好滴，谢谢老师（不太懂无机物，怕给老板搞错了
作者Author: 卡开发发 时间: 2024-2-3 15:44
不好说。一看记录何种数据和转换精度，例如有些按照晶胞参数（三个长度三个角度）有些则按照3x3的晶胞矢量，但两者在转换时也可能导致一些微小差异，或是分数坐标到笛卡尔坐标转换也有类似问题，二是每个格式或者写入该格式的程序的数值精度如何。
事实上相关的格式转换有现成工具，并不需要先向mol格式转换，python生态下ase就能做的很好，但转换或输出精度是否能达到你想要的样子你可以自行验证。
作者Author: 0325-ddd 时间: 2024-2-3 16:26
好的好的，非常感谢
作者Author: Whitedwarf 时间: 2024-4-18 10:17
Sob老师，请问现在还可以要到xyz2QC的源码嘛？我想在MD增强采样到不同cv路径的帧数据中按照数据点密度或者其他筛选算法，筛选出一些数据去Gaussian做能量计算，想参考xyz2QC的写法或者在源码上加一些功能，不知道能不能实现
作者Author: sobereva 时间: 2024-4-19 05:13
可以通过电子邮件向我索取，注明个人身份
作者Author: wilely 时间: 5 day ago
能不能反过来用？即把多个如100个xyz文件合并到一个xyz文件里面，好用xtb一起优化。
作者Author: Uus/pMeC6H4-/キ 时间: 5 day ago
本帖最后由 Uus/pMeC6H4-/キ 于 2024-11-18 17:59 编辑
xtb支持多帧xyz文件结构逐个读取、优化并把结果产生到同一个输出文件吗？看文档似乎提到的.xyz输入文件都是单个结构的。
既然xtb是通过命令行运行，那在shell里用个for循环遍历所有.xyz文件就能用xtb批量处理吧
作者Author: sobereva 时间: 4 day ago
cat *.xyz > all.xyz，就把当前目录下所有xyz文件合并成多帧的all.xyz了
作者Author: wilely 时间: 4 day ago
好的，谢谢
| 欢迎光临 计算化学公社 (http://bbs.keinsci.com/) |
Powered by Discuz! X3.3 |