计算化学公社

 找回密码 Forget password
 注册 Register
Views: 450|回复 Reply: 3

[Multiwfn资源与经验] 批量输出ESP分析信息的脚本

[复制链接 Copy URL]

80

帖子

2

威望

556

eV
积分
676

Level 4 (黑子)

发表于 Post on 2022-9-14 04:13:27 | 显示全部楼层 Show all |阅读模式 Reading model
本帖最后由 pika02 于 2022-10-4 00:12 编辑
  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. echo "Finished"
复制代码

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

参见http://sobereva.com/601

评分 Rate

参与人数
Participants 3
eV +20 收起 理由
Reason
sobereva + 10
七尺贱 + 5 好物!
丁越 + 5 赞!

查看全部评分 View all ratings

80

帖子

2

威望

556

eV
积分
676

Level 4 (黑子)

 楼主 Author| 发表于 Post on 2022-10-4 00:14:29 | 显示全部楼层 Show all
20221004 修了一个bug:输出Polar surface area (|ESP| <= 10 kcal/mol)应为|ESP| > 10 kcal/mol

80

帖子

2

威望

556

eV
积分
676

Level 4 (黑子)

 楼主 Author| 发表于 Post on 2022-9-14 05:16:15 | 显示全部楼层 Show all
其实是自己用来搞化学信息学的东西当场写的

假设你稍微有点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里,请自行更改它为绝对路径。

174

帖子

0

威望

836

eV
积分
1010

Level 4 (黑子)

发表于 Post on 2022-9-14 14:29:57 | 显示全部楼层 Show all
本帖最后由 七尺贱 于 2022-9-14 14:37 编辑
pika02 发表于 2022-9-14 05:16
其实是自己用来搞化学信息学的东西当场写的

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

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

本版积分规则 Credits rule

手机版 Mobile version|北京科音自然科学研究中心 Beijing Kein Research Center for Natural Sciences|京公网安备 11010502035419号|计算化学公社 — 北京科音旗下高水平计算化学交流论坛 ( 京ICP备14038949号-1 )|网站地图

GMT+8, 2023-2-2 23:44 , Processed in 0.905245 second(s), 24 queries .

快速回复 返回顶部 返回列表 Return to list