本帖最后由 丁越 于 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，第一列是构象编号，第二列是相对应能量。 - 47 -2257.94540220
- 58 -2257.94637200
- 52 -2257.94742030
- 48 -2257.94822210
- 55 -2257.94907700
- 8 -2257.95001480
- 4 -2257.95100720
- 25 -2257.95219110
- 21 -2257.95303130
- 38 -2257.95437110
- 56 -2257.95793080
输入./Boltzmann.sh test.txt，然后提示输入温度，这里回车使用默认298.15K。Boltzmann分布的结果便显示在下面了。显示的数据是按照能量由低到高排列的，第二列顺便计算出了每个构象相对于最低能量的构象的差值，单位转化为了Kcal/mol。第三列是由Boltzmann公式Ni/N_ref=exp{(E_ref - E_i)/k_b*T}得到的值；不同构象的分布比例列在第四列了，默认保留了两位小数。 - The Temperature has been set as 298.15 K.
- ################################### Boltzmann Distribution ###################################
- Index deltaE(Kcal/mol) Q_i(relat) Percent(%)
- 56 0.000000 1.000000 96.91
- 38 2.233712 0.023004 2.23
- 21 3.074436 0.005562 0.54
- 25 3.601662 0.002283 0.22
- 4 4.344559 0.000651 0.06
- 8 4.967290 0.000228 0.02
- 55 5.555759 0.000084 0.01
- 48 6.092209 0.000034 0.00
- 52 6.595339 0.000015 0.00
- 58 7.253147 0.000005 0.00
- 47 7.861697 0.000002 0.00
- ##############################################################################################
- **The calcualted Boltzmann distribution has been exported to result.txt in current folder** ^o^
接下来说说脚本代码的构建： - #!/bin/bash
- #Author: Yue Ding
- #This script is used to calculate the Boltzmann distribution of different conformations.The energy unit is defaulted as Hartree (a.u.).
- #Note:Check the input .txt file to make sure it dosen't contain any blank lines!!!!!
- #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
- read -p 'Input the Temperature (K), default 298.15 K by press Enter: ' T
- if [[ -z ${T} ]]; then
- echo "The Temperature has been set as 298.15 K."
- T=298.15
- sleep 0.5
- else
- sleep 0.5
- fi
- echo '################################### Boltzmann Distribution ###################################'
- cat $1 |sort -nk 2 > ./temp1.txt
- E_ref=`cat temp1.txt |awk 'NR==1{print $2}'`
- 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
- cat temp2.txt |awk -v var1="${T}" '{printf "%d\t %f\t %f\n", $1, $2, exp($3/var1)}' > ./temp3.txt
- sum=`cut -f 3 'temp3.txt'|awk '{all += $1}END {print all}'`
- 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
- echo '##############################################################################################'
- echo '**The calcualted Boltzmann distribution has been exported to result.txt in current folder** ^o^ '
- 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文件。
| 
发表于 Post on 2021-12-4 22:52:59
楼主 Author
