gaussian 输入 gjf 转成pdb文件

meatball1982

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

wat.gjf

1. % nproc=8
   
2. % mem= 12GB
   
3. 
   
4. # mp2/6-31+G*
   
5. 
   
6. wat
   
7. 
   
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
    
12. 
    

复制代码

我希望得到对应的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.

谢谢大家。

赵云跳槽

如果是我做的话，我可以先用高斯计算单点（最低计算等级，比如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

复制代码

sobereva

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

其实不如用guess(save,only)，这样只做一次初猜任务而已，不花时间

winterzen

用newzmat转换时，Cartesian坐标的gjf的语法是：

1. newzmat -icart xxx.gjf -opdb xxx.pdb

复制代码

ChemiAndy

用newzmat的确很优雅。

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

$ cat gjf2xyz

1. #!/bin/bash
   
2. 
   
3. PREFIX=`echo $1| sed 's/\..*$//'`
   
4. echo "PREFIX of input file:" $PREFIX
   
5. 
   
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/://'`
   
9. 
   
10. N_START=$((M_START+1))
    
11. N_END=$((M_END-1))
    
12. echo "N_START:" $N_START
    
13. echo "N_END:" $N_END
    
14. 
    
15. echo "$((N_END - N_START))" > ${PREFIX}.xyz
    
16. cat $1 | sed -n "${N_START},${N_END}p" >> ${PREFIX}.xyz
    
17. 
    
18. echo "Please check the xyz file:"
    
19. echo "--------------------------------------------------------------"
    
20. cat ${PREFIX}.xyz
    
21. echo "--------------------------------------------------------------"
    
22. echo "Done!"
    

复制代码

meatball1982

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

这一思路倒是想过，只计算很小的基组和方法，我的体系有点大，关键是有点多。所有当时放弃了。还是谢谢您的思路。

meatball1982

sobereva 发表于 2018-4-27 12:28
其实不如用guess(save,only)，这样只做一次初猜任务而已，不花时间

准备尝试一下。

meatball1982

winterzen 发表于 2018-4-27 13:10
用newzmat转换时，Cartesian坐标的gjf的语法是：

针对一楼给出的wat.gjf
似乎还是不行。

winterzen

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

meatball1982

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那部分有点小问题，还是说输入应该满足一定的形式？

meatball1982

winterzen 发表于 2018-4-27 15:12
我这里没有问题，转换以后的结果是

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

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

ChemiAndy

meatball1982 发表于 2018-4-27 15:16
非常感谢ChemiAndy

我将代码存成　

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

1. #!/bin/bash
   
2. 
   
3. PREFIX=`echo $1| sed 's/\..*$//'`
   
4. 
   
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/://')
   
8. 
   
9. N_START=$((M_START+1))
   
10. N_END=$((M_END-1))
    
11. 
    
12. 
    
13. echo "$((N_END - N_START))" > ${PREFIX}.xyz
    
14. sed -n '/^#/,$p' $1 | sed -n "${N_START},${N_END}p" >> ${PREFIX}.xyz
    
15. 
    
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脚本是十分强大的，值得花点时间摸索一下。

meatball1982

ChemiAndy 发表于 2018-4-28 17:10
我所有的gjf在计算设置行（#开始的行）之前都没有空行，所以我的脚本没有考虑这种意外。我针对你给出的gj ...

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


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

1. #!/bin/bash
   
2. #
   
3. 
   
4. ## outline #################################################
   
5. # ./sh_gjf2xyz.sh wat.gjf
   
6. #
   
7. 
   
8. 
   
9. ## main ###################################################
   
10. PREFIX=`echo $1| sed 's/\..*$//'`
    
11. 
    
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/://')
    
15. 
    
16. N_START=$((M_START+1))
    
17. N_END=$((M_END-1))
    
18. 
    
19. echo "$((N_END - N_START))" > ${PREFIX}.xyz
    
20. sed -n '/^#/,$p' $1 | sed -n "${N_START},${N_END}p" >> ${PREFIX}.xyz
    
21. 
    
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

复制代码

Daniel_Arndt

以前的课程设计作业中涉及类似的文件格式转换，自学了一下午awk写了个脚本完成的。真是233

ChemiAndy

meatball1982 发表于 2018-5-1 19:45
非常感谢您的帮助。
我尝试一下。
另，已经在认真学习sed和grep了。