批量输出ESP分析信息的脚本

pika02

1. #!/bin/bash
   
2. # define a function to get value, usage: 'string' number (defaut 1)
   
3. getvalue(){
   
4.     grep -m 1 "${1}" $raw |\
   
5.     grep -Eo ' -?[0-9]+\.*[0-9]+' |\
   
6.     sed -n "${2:-1}p" |\
   
7.     sed 's/ //g'
   
8.     }
   
9. echo "Script for generating ESP relative descriptors from multiple gbw files."
   
10. tmp=$(mktemp)
    
11. list=$(mktemp)
    
12. cat << EOF > $tmp
    
13. 12
    
14. 0
    
15. -1
    
16. -1
    
17. q
    
18. EOF
    
19. # for gbws in *.gbw
    
20. #     do
    
21. #     echo "Convert from ${gbws} to wfn..."
    
22. #     $HOME/apps/orca4/orca_2aim ${gbws%.gbw}
    
23. #     rm -f *.wfx
    
24. #     echo "done"
    
25. #     done
    
26. echo "\
    
27. file,\
    
28. Volume (Angstrom^3),\
    
29. density (g/cm^3),\
    
30. Minimal value (kcal/mol),\
    
31. Maximal value (kcal/mol),\
    
32. Overall surface area (Angstrom^2),\
    
33. Positive surface area (Angstrom^2),\
    
34. Negative surface area (Angstrom^2),\
    
35. Overall average value (kcal/mol),\
    
36. Positive average value (kcal/mol),\
    
37. Negative average value (kcal/mol),\
    
38. Overall variance (sigma^2_tot) (kcal^2/mol^2),\
    
39. Positive variance (kcal^2/mol^2),\
    
40. Negative variance (kcal^2/mol^2),\
    
41. Balance of charges (nu),\
    
42. Product of sigma^2_tot and nu (kcal/mol),\
    
43. Internal charge separation (Pi) (kcal/mol),\
    
44. Molecular polarity index (MPI) (kcal/mol),\
    
45. Nonpolar surface area (|ESP| <= 10 kcal/mol) (Angstrom^2),\
    
46. Nonpolar surface area (%),\
    
47. Polar surface area (|ESP| > 10 kcal/mol) (Angstrom^2),\
    
48. Polar surface area (%)\
    
49. " > "$list"
    
50. for file in *.wfn
    
51.     do
    
52.     echo "Calculating descriptors from $file ..."
    
53.     raw=$(mktemp)
    
54.     Multiwfn $file < $tmp > $raw
    
55.     volume=$(getvalue 'Volume:' 2)
    
56.     density=$(getvalue 'Estimated density')
    
57.     minimal=$(getvalue 'Minimal' 1)
    
58.     maximal=$(getvalue 'Maximal' 2)
    
59.     osa=$(getvalue 'Overall surface area' 2)
    
60.     psa=$(getvalue 'Positive surface area' 2)
    
61.     nsa=$(getvalue 'Negative surface area' 2)
    
62.     oav=$(getvalue 'Overall average value' 2)
    
63.     pav=$(getvalue 'Positive average value' 2)
    
64.     nav=$(getvalue 'Negative average value' 2)
    
65.     ov=$(getvalue 'Overall variance' 2)
    
66.     pv=$(getvalue 'Positive variance' 2)
    
67.     nv=$(getvalue 'Negative variance' 2)
    
68.     bc=$(getvalue 'Balance of charges')
    
69.     psn=$(getvalue 'Product of sigma^2_tot and nu' 2)
    
70.     ics=$(getvalue 'Internal charge separation' 2)
    
71.     mpi=$(getvalue 'Molecular polarity index' 2)
    
72.     ns=$(getvalue 'Nonpolar surface area' 2)
    
73.     nsp=$(getvalue 'Nonpolar surface area' 3)
    
74.     ps=$(getvalue 'Polar surface area' 2)
    
75.     psp=$(getvalue 'Polar surface area' 3)
    
76.     echo "\
    
77. $file,\
    
78. $volume,\
    
79. $density,\
    
80. $minimal,\
    
81. $maximal,\
    
82. $osa,\
    
83. $psa,\
    
84. $nsa,\
    
85. $oav,\
    
86. $pav,\
    
87. $nav,\
    
88. $ov,\
    
89. $pv,\
    
90. $nv,\
    
91. $bc,\
    
92. $psn,\
    
93. $ics,\
    
94. $mpi,\
    
95. $ns,\
    
96. $nsp,\
    
97. $ps,\
    
98. $psp\
    
99. " >> $list
    
100.     done
     
101. rm $tmp
     
102. mv $list descriptors.csv
     
103. 
     
104. echo "Finished"

复制代码

改自http://bbs.keinsci.com/thread-12066-1-1.html，注释掉了批量转换gbw的功能

参见http://sobereva.com/601

pika02

20221004 修了一个bug：输出Polar surface area (|ESP| <= 10 kcal/mol)应为|ESP| > 10 kcal/mol

pika02

其实是自己用来搞化学信息学的东西当场写的

假设你稍微有点linux基础，我来废话一下这个脚本是怎么写的

3行定义了一个getvalue函数，调用grep根据输入的字符串从Multiwfn的文字输出（stdout）里找到数据所在行，并用正则表达式把里面的数字（包括前置空格，负号）取出来。因为Multiwfn的ESP输出信息里有些数字是单位的一部分（如Bohr^3），真正用得到的数字前面一定有空格（如果存在数值前面没带空格的情况，请跟我说），可能有负号和小数点。getvalue的第一个参数是要找数据用的字符串，第二个参数是想要的值在第几号位，默认是1。
不过有一点恼火的是，Nonpolar surface area (|ESP| <= 10 kcal/mol):    114.51 Angstrom^2  ( 19.19 %)这样的行，匹配正则表达式1号位的数字是10

10行，mktemp是一个很有用的生成临时文件的命令，默认会在/tmp/目录下生成一个叫“tmp.一串随机字符”的临时文件，并输出临时文件的绝对路径。一个变量名=$(mktemp)这样的操作在写shell脚本时非常好用，能瞬间生成一个不碍眼的临时文件，并在脚本以后的地方直接写 "$一个变量" 来调用。由于许多linux的/tmp都在内存里，关机就消失了，因此不需要再管它。这里的tmp是控制Multiwfn功能的神秘文本文件。11行的list也是一样的操作，先往临时文件里面写数据，最后把它移到当前目录，能有效减少硬盘读写操作。如果担心算到一半断电，那就乖乖让list指向硬盘上的文件吧。

12行生成一个神秘文本文件扔到刚刚创建的临时文件里，下面会用到。26-49行是写入表头，有两个单位的数值我喜欢拿第二个单位来用，最后nonpolar和polar surface area的数值和占比我全都要。

54行就是调用Multiwfn算东西了。轮流把当前目录下的.wfn文件喂给Multiwfn，并用刚刚的神秘tmp文件告诉Multiwfn该干什么，最后输出到一个新的临时文件raw里。然后调用一开始定义的函数从raw里面拿数据，再写进表里。当然，如果Multiwfn不在你的PATH里，请自行更改它为绝对路径。

七尺贱

pika02 发表于 2022-9-14 05:16
其实是自己用来搞化学信息学的东西当场写的

假设你稍微有点linux基础，我来废话一下这个脚本是怎么 ...

ESP分析完之后会有正的极小值点，负的极大值点，这些点一般没啥用，可以用后处理菜单的3和4，输入d来删掉这些点，而且我感觉算的东西放前面，输出的东西放后面这样可读性好一点？我看了半天才知道后面怎么算