计算化学公社

 找回密码 Forget password
 注册 Register
Views: 1856|回复 Reply: 5

[算法与编程] 一键计算Boltzmann分布的脚本Boltzmann.sh

[复制链接 Copy URL]

327

帖子

9

威望

1992

eV
积分
2499

Level 5 (御坂)

发表于 Post on 2021-12-4 22:52:59 | 显示全部楼层 Show all |阅读模式 Reading model
本帖最后由 丁越 于 2021-12-5 14:01 编辑

一键计算Boltzmann分布的脚本Boltzmann.sh

    今天在计算一个有多重构象小分子的Raman光谱时,感觉用Excel一个个的求不同构象的Boltzmann分布比例太费劲了,遂花了点时间写了小脚本,使之可以一键求得某温度下的分布比例。由于目前常用的Gaussian、CP2K的能量单位默认为Hartree,所以使用脚本时切记不要再去转化为Kcal/mol之类的单位。下面就来说说脚本的使用:
    首先建立一个.txt后缀的文本文件,文件第一列写入构象的Index,这个Index数值是任意的,目的只是为了标识不同的构象; 然后tab空一列,第二列写入不同构象的能量,粗糙计算的话可以直接用电子能量,讲究的话最好使用自由能。非常重要的一点:这个文本文件不要留有任何空行!!!!
    使用:./Boltzmann xxx.txt

举个例子:下面是我计算的那个分子能量最低的11个构象,输入文件命名为test.txt,第一列是构象编号,第二列是相对应能量。
  1. 47       -2257.94540220
  2. 58       -2257.94637200
  3. 52       -2257.94742030
  4. 48       -2257.94822210
  5. 55       -2257.94907700
  6. 8        -2257.95001480
  7. 4        -2257.95100720
  8. 25       -2257.95219110
  9. 21       -2257.95303130
  10. 38       -2257.95437110
  11. 56       -2257.95793080
复制代码
       输入./Boltzmann.sh test.txt,然后提示输入温度,这里回车使用默认298.15K。Boltzmann分布的结果便显示在下面了。显示的数据是按照能量由低到高排列的,第二列顺便计算出了每个构象相对于最低能量的构象的差值,单位转化为了Kcal/mol。第三列是由Boltzmann公式Ni/N_ref=exp{(E_ref - E_i)/k_b*T}得到的值;不同构象的分布比例列在第四列了,默认保留了两位小数。
  1. The Temperature has been set as 298.15 K.
  2. ################################### Boltzmann Distribution ###################################
  3.           Index    deltaE(Kcal/mol)            Q_i(relat)            Percent(%)
  4.              56         0.000000                1.000000              96.91
  5.              38         2.233712                0.023004               2.23
  6.              21         3.074436                0.005562               0.54
  7.              25         3.601662                0.002283               0.22
  8.               4         4.344559                0.000651               0.06
  9.               8         4.967290                0.000228               0.02
  10.              55         5.555759                0.000084               0.01
  11.              48         6.092209                0.000034               0.00
  12.              52         6.595339                0.000015               0.00
  13.              58         7.253147                0.000005               0.00
  14.              47         7.861697                0.000002               0.00
  15. ##############################################################################################
  16. **The calcualted Boltzmann distribution has been exported to result.txt in current folder** ^o^
复制代码

接下来说说脚本代码的构建:
  1. #!/bin/bash
  2. #Author: Yue Ding
  3. #This script is used to calculate the Boltzmann distribution of different conformations.The energy unit is defaulted as Hartree (a.u.).
  4. #Note:Check the input .txt file to make sure it dosen't contain any blank lines!!!!!
  5. #Usage: Creating a .txt file that contains the index and hartree energy of conformations, and then you can calculate it just by inputing ./Boltzmann xxx.txt

  6. read -p 'Input the Temperature (K), default 298.15 K by press Enter: ' T
  7. if [[ -z ${T} ]]; then
  8.   echo "The Temperature has been set as 298.15 K."
  9.   T=298.15
  10.   sleep 0.5
  11. else
  12.         sleep 0.5
  13. fi

  14. echo '################################### Boltzmann Distribution ###################################'
  15. cat $1 |sort -nk 2 > ./temp1.txt
  16. E_ref=`cat temp1.txt |awk 'NR==1{print $2}'`
  17. cat temp1.txt |awk -v var="${E_ref}" '{printf "%d\t %f\t %f\n", $1, ($2-var)*627.5, (var-$2)*315940}' > ./temp2.txt
  18. cat temp2.txt |awk -v var1="${T}" '{printf "%d\t %f\t %f\n", $1, $2, exp($3/var1)}' > ./temp3.txt
  19. sum=`cut -f 3 'temp3.txt'|awk '{all += $1}END {print all}'`
  20. cat temp3.txt |awk -v var2="${sum}" 'NR==1{printf "%15s\t %18s\t %16s\t %14s\n", "Index","deltaE(Kcal/mol)","Q_i(relat)","Percent(%)"} NR>=1{printf "%15d\t %15.6f\t %15.6f\t %10.2f\n", $1, $2, $3, ($3/var2)*100}'|tee reslut.txt
  21. echo '##############################################################################################'
  22. echo '**The calcualted Boltzmann distribution has been exported to result.txt in current folder** ^o^ '
  23. rm -f temp*.txt
复制代码

首先需要把test.txt文件交给Boltzmann脚本读取,那么$1对应着脚本接收到的第一个参数。紧接着需要将第中第二列数据由小到大排列,所以使用了cat $1 |sort -nk 2,注意对数字排列大小时一定要带-n(我当时没带-n,快要被离谱的排序搞到吐血。。);-k指对第几列数据进行排序。
接下来就是定义能量最低值作为参考值,即E_ref=`cat temp1.txt |awk 'NR==1{print $2}'`
这个temp2.txt文件的作用就是记录三列数据,分别是构象Index、转化单位为Kcal/mol后得到的能量差值deltaE(Kcal/mol)、以及Boltzmann公式Ni/N_ref=exp{(E_ref - E_i)/k_b*T}中的(E_ref - E_i)/k_b的值(都是单位换算问题)。awk -v var="${E_ref}" 中-v(variable的简写) 的作用是可以在awk中读取外部定义的变量,我觉得非常有用,这里我们将外部定义的E_ref读取进来然后赋值给变量var。
第三个temp3.txt的作用就是记录构象Index、转化单位得到的能量差值deltaE(Kcal/mol)、以及运算exp{(E_ref - E_i)/k_b*T}。
sum=`cut -f 3 'temp3.txt'|awk '{all += $1}END {print all}'` 这一列就是求exp{(E_ref - E_i)/k_b*T}的加和,也即 Q_i(relat)的加和,并且由此可以计算不同构象的比例了。
接下来的一大长串就是调整每列数据的间隔并且给他加个标题。
最后删除无用的临时temp文件。

Boltzmann.sh

1.47 KB, 下载次数 Times of downloads: 38

评分 Rate

参与人数
Participants 5
eV +27 收起 理由
Reason
Jiongci-21 + 2 谢谢
Para1lel + 5 谢谢分享
ggdh + 5 GJ!
zsu007 + 5
sobereva + 10

查看全部评分 View all ratings

自由发挥,野蛮生长

877

帖子

36

威望

4803

eV
积分
6400

Level 6 (一方通行)

发表于 Post on 2021-12-5 12:14:04 | 显示全部楼层 Show all
我又来提意见了。。。
因为默认温度就是298.15,而且一般算Boltzmann都用这个温度,所以建议设一个默认值,就是如果用户直接回车,用默认温度,而不是强迫用户输入298.15。

327

帖子

9

威望

1992

eV
积分
2499

Level 5 (御坂)

 楼主 Author| 发表于 Post on 2021-12-5 13:35:57 | 显示全部楼层 Show all
ggdh 发表于 2021-12-5 12:14
我又来提意见了。。。
因为默认温度就是298.15,而且一般算Boltzmann都用这个温度,所以建议设一个默认值 ...

欧了钟叔,我改下脚本
自由发挥,野蛮生长

1

帖子

0

威望

7

eV
积分
8

Level 1 能力者

发表于 Post on 2022-11-17 10:38:23 | 显示全部楼层 Show all
大神好,怎么运行您这个Boltzmann.sh文件呢?怎么调用?

6691

帖子

0

威望

4012

eV
积分
10703

Level 6 (一方通行)

发表于 Post on 2022-11-17 15:21:35 | 显示全部楼层 Show all
gca1983 发表于 2022-11-17 03:38
大神好,怎么运行您这个Boltzmann.sh文件呢?怎么调用?

自己用搜索引擎搜shell脚本的调用方法,一步一步照着做
BDF(https://bdf-manual.readthedocs.io/zh_CN/latest/Introduction.html)、ORCA(https://orcaforum.kofo.mpg.de/index.php)开发团队成员

124

帖子

0

威望

950

eV
积分
1074

Level 4 (黑子)

发表于 Post on 2023-1-7 22:06:48 | 显示全部楼层 Show all
老师,您好,请教您,我运行了脚本,可是没有任何反应是什么原因呢?如图,麻烦您啦
202301072205404624..png

本版积分规则 Credits rule

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

GMT+8, 2023-2-6 04:26 , Processed in 0.358494 second(s), 25 queries .

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