计算化学公社

标题: gaussian 输入 gjf 转成pdb文件 [打印本页]
作者
Author: meatball1982    时间: 2018-4-26 21:01
标题: gaussian 输入 gjf 转成pdb文件
本帖最后由 meatball1982 于 2018-4-26 21:03 编辑

大家好:
       我现在有许多的gaussian的输入文件，其中一个如下

wat.gjf
  1. % nproc=8
  2. % mem= 12GB

  4. # mp2/6-31+G*

  6. wat

  8. 0,1
  9. O                  -0.00000   0.00000   6.84784
  10. H                  -0.78400   0.00000   6.29384
  11. H                   0.78400   0.00000   6.29384

复制代码


我希望得到对应的pdb文件。
方法1，gaussian view,导入，导出成pdb。 因为我的输入有点多，所以手动，似乎不太现实。

方法2。我用gaussian手册中的newzmat 实现。
  1. newzmat -icom wat.gjf -opdb wat.pdb
复制代码
得到的wat.pdb没有原子的位置信息。如下。
  1. b
  2. HEADER
  3. TITLE     wat
  4. SEQRES   0      0              
  5. CONECT    0    0    0
  6. CONECT    0    0
  7. CONECT    0    0
  8. END
复制代码


方法3。组里的老师说 openbabel可以作。尝试了一下。
  1. obabel -i g03 wat.gjf -o pdb wat.pdb
复制代码


显示是
  1. 0 molecules converted
复制代码
我觉得应该还是gaussian的输入有问题。但输入是可以正常跑，并且结束的。


我的问题是：
如何才能用命令行批量实现gjf 到pdb的转换。方法2，方法3 都行，
如果不行，在只有log的情况下（没有chk),能不能通过log 转成pdb.
再不行，是不是只能自己将位置和原子信息提出来，自己再写pdb.

谢谢大家。







作者
Author: 赵云跳槽    时间: 2018-4-27 11:49
如果是我做的话，我可以先用高斯计算单点（最低计算等级，比如HF/sto-2g）得到fchk，再使用Multiwfn转换为pdb
  1. for i in *gjf
  2. do
  3. g09 $i
  4. sleep 1m   <这个按计算任务大小调整>
  5. formchk ${i%gjf}chk
  6. echo -e '100\n2\n1\n'${i%gjf}pbd'' | Multiwfn ${i%gjf}fchk > null
  7. done
  8. rm -rf null
复制代码


作者
Author: sobereva    时间: 2018-4-27 12:28
赵云跳槽 发表于 2018-4-27 11:49
如果是我做的话，我可以先用高斯计算单点（最低计算等级，比如HF/sto-2g）得到fchk，再使用Multiwfn转换为p ...


其实不如用guess(save,only)，这样只做一次初猜任务而已，不花时间
作者
Author: winterzen    时间: 2018-4-27 13:10
用newzmat转换时，Cartesian坐标的gjf的语法是：
  1. newzmat -icart xxx.gjf -opdb xxx.pdb
复制代码
作者
Author: ChemiAndy    时间: 2018-4-27 14:22
用newzmat的确很优雅。

我之前遇到了需要批量转换gjf成pdb这个问题，也是首先尝试了方法3，发现openbabel转换gjf有问题。但是用openbabel将xyz转换成pdb却非常容易，于是写了一个gjf2xyz的小script，用bash批处理很方便。后来这个gjf2xyz成了我工具箱里的常用工具。分享一下：

$ cat gjf2xyz
  1. #!/bin/bash

  3. PREFIX=`echo $1| sed 's/\..*$//'`
  4. echo "PREFIX of input file:" $PREFIX

  6. #Start line number
  7. M_START=`grep -n '^\s* $1 | head -2 | tail -1 | sed 's/://'`
  8. M_END=`grep -n '^\s*   $1 | head -3 | tail -1 | sed 's/://'`

  10. N_START=$((M_START+1))
  11. N_END=$((M_END-1))
  12. echo "N_START:" $N_START
  13. echo "N_END:" $N_END

  15. echo "$((N_END - N_START))" > ${PREFIX}.xyz
  16. cat $1 | sed -n "${N_START},${N_END}p" >> ${PREFIX}.xyz

  18. echo "Please check the xyz file:"
  19. echo "--------------------------------------------------------------"
  20. cat ${PREFIX}.xyz
  21. echo "--------------------------------------------------------------"
  22. echo "Done!"
复制代码
作者
Author: meatball1982    时间: 2018-4-27 14:56
赵云跳槽 发表于 2018-4-27 11:49
如果是我做的话，我可以先用高斯计算单点（最低计算等级，比如HF/sto-2g）得到fchk，再使用Multiwfn转换为p ...

这一思路倒是想过，只计算很小的基组和方法，我的体系有点大，关键是有点多。所有当时放弃了。还是谢谢您的思路。
作者
Author: meatball1982    时间: 2018-4-27 14:57
sobereva 发表于 2018-4-27 12:28
其实不如用guess(save,only)，这样只做一次初猜任务而已，不花时间

准备尝试一下。
作者
Author: meatball1982    时间: 2018-4-27 14:58
winterzen 发表于 2018-4-27 13:10
用newzmat转换时，Cartesian坐标的gjf的语法是：

针对一楼给出的wat.gjf
似乎还是不行。
作者
Author: winterzen    时间: 2018-4-27 15:12
meatball1982 发表于 2018-4-27 14:58
针对一楼给出的wat.gjf
似乎还是不行。

我这里没有问题，转换以后的结果是
  1. HEADER
  2. REMARK wat
  3. HETATM    1  O           1       0.000   0.000   6.848                       O
  4. HETATM    2  H           1      -0.784   0.000   6.294                       H
  5. HETATM    3  H           1       0.784   0.000   6.294                       H
  6. CONECT    1    2    3
  7. CONECT    2    1
  8. CONECT    3    1
  9. END
复制代码


用的Gaussian 09 D.01版中的newzmat
作者
Author: meatball1982    时间: 2018-4-27 15:16
ChemiAndy 发表于 2018-4-27 14:22
用newzmat的确很优雅。

我之前遇到了需要批量转换gjf成pdb这个问题，也是首先尝试了方法3，发现openbabe ...

非常感谢ChemiAndy

我将代码存成　
gjf2xyz
./gjf2xyz wat.gjf

wat.gjf 为１Ｌ的输入。

得到
  1. PREFIX of input file: wat
  2. ./gjf2xyz: command substitution: line 7: unexpected EOF while looking for matching `''
  3. ./gjf2xyz: command substitution: line 8: syntax error: unexpected end of file
  4. ./gjf2xyz: command substitution: line 8: unexpected EOF while looking for matching `''
  5. ./gjf2xyz: command substitution: line 9: syntax error: unexpected end of file
  6. N_START: 1
  7. N_END: -1
  8. sed: -e expression #1, char 3: unexpected `,'
  9. Please check the xyz file:
  10. --------------------------------------------------------------
  11. -2
  12. --------------------------------------------------------------
  13. Done
复制代码


是不是在找M_START和M_END那部分有点小问题，还是说输入应该满足一定的形式？
作者
Author: meatball1982    时间: 2018-4-27 16:23
本帖最后由 meatball1982 于 2018-5-2 09:57 编辑
winterzen 发表于 2018-4-27 15:12
我这里没有问题，转换以后的结果是

好的。我再找一个新一点Gaussian试试。非常感谢。

更新： 2018年 05月 02日 星期三 09:56:52 CST
换了个新版本的gaussian, g16.
可行。

作者
Author: ChemiAndy    时间: 2018-4-28 17:10
本帖最后由 ChemiAndy 于 2018-4-28 17:25 编辑
meatball1982 发表于 2018-4-27 15:16
非常感谢ChemiAndy

我将代码存成　

我所有的gjf在计算设置行（#开始的行）之前都没有空行，所以我的脚本没有考虑这种意外。我针对你给出的gjf文件修改了脚本，考虑了存在空行的可能性。新脚本
  1. #!/bin/bash

  3. PREFIX=`echo $1| sed 's/\..*$//'`

  5. #Start line number
  6. M_START=$(sed -n '/^#/,$p' $1 | grep -n '^\s* | head -2 | tail -1 | sed 's/://')
  7. M_END=$(sed -n '/^#/,$p' $1   | grep -n '^\s*  | head -3 | tail -1 | sed 's/://')

  9. N_START=$((M_START+1))
  10. N_END=$((M_END-1))


  13. echo "$((N_END - N_START))" > ${PREFIX}.xyz
  14. sed -n '/^#/,$p' $1 | sed -n "${N_START},${N_END}p" >> ${PREFIX}.xyz

  16. echo "Please check the generated xyz file:"
  17. echo "--------------------------------------------------------------"
  18. cat ${PREFIX}.xyz
  19. echo "--------------------------------------------------------------"
  20. echo "Done!"
复制代码


运行结果：
  1. ./gjf2xyz wat.gjf
  2. Please check the generated xyz file:
  3. --------------------------------------------------------------
  4. 3
  5. 0,1
  6. O                  -0.00000   0.00000   6.84784
  7. H                  -0.78400   0.00000   6.29384
  8. H                   0.78400   0.00000   6.29384
  9. --------------------------------------------------------------
复制代码


解释：
1. 这个脚本就是提取#行之后的所有行，然后根据gaussian输入文件的空行规则，从第二个空行的第二行开始，到第三个空行之前一行结束，为gjf的xyz坐标行（对Z坐标无效），直接输出。
2. xyz文件的前两行：一般第一行第一个数字为总原子数，第二行无规则，我这里直接把总电荷和自旋数写上，以便输出成其他软件的输入文件时使用；

你试下，有错误告诉我。
正则表达式和sed/grep结合使用编写shell脚本是十分强大的，值得花点时间摸索一下。
作者
Author: meatball1982    时间: 2018-5-1 19:45
本帖最后由 meatball1982 于 2018-5-2 14:51 编辑
ChemiAndy 发表于 2018-4-28 17:10
我所有的gjf在计算设置行（#开始的行）之前都没有空行，所以我的脚本没有考虑这种意外。我针对你给出的gj ...

非常感谢您的帮助。
我尝试一下。
另，已经在认真学习sed和grep了。


更新：
根据您提出的思路，我修改了一个简单版本的。
通过脚本 gjf -> xyz, xyz -> pdb

  1. #!/bin/bash
  2. #

  4. ## outline #################################################
  5. # ./sh_gjf2xyz.sh wat.gjf
  6. #


  9. ## main ###################################################
  10. PREFIX=`echo $1| sed 's/\..*$//'`

  12. #Start line number
  13. M_START=$(sed -n '/^#/,$p' $1 | grep -n '^ | head -2 | tail -1 | sed 's/://')
  14. M_END=$(sed   -n '/^#/,$p' $1 | grep -n '^   | head -3 | tail -1 | sed 's/://')

  16. N_START=$((M_START+1))
  17. N_END=$((M_END-1))

  19. echo "$((N_END - N_START))" > ${PREFIX}.xyz
  20. sed -n '/^#/,$p' $1 | sed -n "${N_START},${N_END}p" >> ${PREFIX}.xyz

  22. ## logs ###################################################
  23. #
  24. # http://bbs.keinsci.com/thread-9826-1-1.html#pid67064
  25. # typed by : ChemiAndy
  26. # mod   by : meatball1982
  27. # mod   on : 2018年 05月 02日 星期三 11:48:02 CST
复制代码

  1. babel -ixyz wat.xyz  -opdb wat.pdb
复制代码




作者
Author: Daniel_Arndt    时间: 2018-5-2 13:07
以前的课程设计作业中涉及类似的文件格式转换，自学了一下午awk写了个脚本完成的。真是233
作者
Author: ChemiAndy    时间: 2018-5-4 14:21
meatball1982 发表于 2018-5-1 19:45
非常感谢您的帮助。
我尝试一下。
另，已经在认真学习sed和grep了。

作者
Author: cokoy    时间: 2019-3-7 11:00
obabel -i g03 wat.gjf -o pdb wat.pdb
命令错了，应该是
obabel -i g03 wat.gjf -o pdb -O wat.pdb
作者
Author: meatball1982    时间: 2019-3-7 21:58
cokoy 发表于 2019-3-7 11:00
obabel -i g03 wat.gjf -o pdb wat.pdb
命令错了，应该是
obabel -i g03 wat.gjf -o pdb -O wat.pdb

好的。谢谢，我研究一下。
好长时间的帖子了。得回忆一下。



欢迎光临 计算化学公社 (http://bbs.keinsci.com/) Powered by Discuz! X3.3