Bash

Bash—Bourne shell

一种很棒的脚本语言,处理文本文件十分方便。故记之,以备。
注意点awk<强大的脚本编程工具,Linux三剑客中的王者>:

  • 数组默认从1开始,0开始需手动指定
  • 自定义函数时候,函数内部参数尽量与外部参数设置不同
  • 函数内部出现的参数,且未在函数参数列表中的均为全局变量。这点尤其要注意!
  • 需要通过函数返回多个值,请用数组。
  • awk中引入数值记得+0,比如substr()+0"'$a'"+0
  • 采用输出重定向>>写入文件以后,如果紧接着后面还需要用getline打开读取,一定要记得先close前面的文件。
  • 写入文件内容较多时,可以用fflush(filename)来刷新缓存,便于及时查看输出文件内容是否正确。
  • 自定义函数中尽量减少中间变量的使用。特别是未初始化的变量,将其加入到函数参数列表中以防变量污染。

简单碳纳米管构建脚本

一个简单的bash脚本用于构建CNT,支持输出文件格式xyzgro。如果使用gro,则盒子Z方向为周期性。可以用gmx editconf -c进行居中处理即可。下载链接

叠合并计算两结构之间的RMSD值

直接读取含有多帧结构构型的xyz或者gro文件,然后进行叠合并以第一帧结构为参考结构,计算出其他帧相对于参考结构的RMSD值,单位为或者纳米,取决于使用的输入文件类型。
需要注意几点:

  • 目前只能够计算相同结构的RMSD,原子数和元素必须对应
  • 对所有原子进行拟合,包括H,主流原子CHONS
  • 支持原子质量权重

使用实例;
对结构进行质量权重

1
bash awk_fitmol.bsh -m y -i fit.xyz


对结构不用质量权重

1
bash awk_fitmol.bsh -m n -i fit.xyz

该结果与VMD一致。脚本和实例下载链接

自相关函数计算脚本

本脚本能够计算速度自相关函数等,需要提供一个输入xvg文件,该文件能够通过gmx traj -ov生成。自动忽略含有[@#&]字符的行,第一列数据必须是时间,后面每三列为同类粒子的速度分量。脚本未采用FFT变换,速度欠佳。处理小几万帧还可以。下载链接。结果展示,如下图分别用该脚本与gmx velacc工具比较:

特殊构型创建脚本

主要是结合gmx select命令来构建一些特殊构型的溶剂体系,比如球体,椭球体,圆柱体,板状等。可以参考博文构建特殊构型的GROMACS脚本。脚本和实例下载链接

氧化石墨烯构建脚本

目前有很多的工具可以用来创建石墨烯和碳纳米管模型,参看博文无机材料的基本操作,个别还可以创建氧化石墨烯【包括-OH,-COOH和-O-基团】。但是毕竟是网页版的,可能也不太十分的方便,因此这里我也提供一个创建氧化石墨烯的awk脚本,可以添加-OH-COOH-COO--O-和边界-H基团,输出结构文件格式为.gro,便于GROMACS使用。功能基团默认的键长和键角来自使用xtb优化的氧化石墨烯模型。下载链接

转换gmx anaeig -2d输出xvg为概率分布或者自由能

这个源程序来自于博文浅谈PCA与g_covar+g_anaeig+ddtdp+sigmaplot做自由能面图的方法,原来是个Fortran程序,但是我看了一下,觉得处理这样的数据用awk完成也是十分方便的,因此使用awk重写了一下这个程序。下载链接。但是我觉得这个工具使用的思路不仅限于这个用途,你自己也可以改一改用于其他用途,比如有人问到的如何将拉氏图的两个二面角($\phi$, $\psi$)与自由能(G)联系起来(文献中有这样的图)。

拉氏图绘图脚本

这本是个python脚本,暂且放到这里吧。用于gmx rama输出作图,需要先将输出的rama.xvg使用egrep命令提取出单个氨基酸所有帧的$\phi$和$\psi$,然后进行MD热图的绘制。下载链接

STRIDE后处理统计工具

STRIDEDSSP相似,都是用来计算蛋白二级结构的,这里也给出一个处理STRIDE输出文件的脚本,用于统计各种二级结构在总氨基酸中所占比例。下载链接链接

DSSP后处理统计工具

一个简单的脚本,读取DSSP输出,统计各种二级结构在总氨基酸中所占比例,下载链接。适用于一帧pdb结构文件。完整流程:

1
2
dssp -i file.pdb -o ss.dssp
bash awk_dssp_ss_count.bsh ss.dssp

输出:

1
2
3
4
5
6
7
8
9
10
11
This Protein has  2 Chains,
Total Number of Residue is 196,
Each Structure percentage as follows:

B-Sheet is: 48.469%
B-Bridge is: 2.041%
Bend is: 13.776%
Turn is: 7.143%
A-Helix is: 5.612%
3-Helix is: 1.531%
Coil is: 21.429%

平均滞留时间脚本

这个awk脚本用来处理gmx select -om命令产生的mask.xvg,原始博文参考GROMACS分析教程:使用g_select计算平均滞留时间 ,此方法是基于该文的最后部分Matlab处理方式更改的,测试速度要比博文中提及的awk代码要快,但是两者的结果有一点点差距,当然我也写过C版本的,那个最快。链接

简单的石墨烯结构产生脚本

一个简单的建模石墨烯脚本。下载链接

单个线性分子组成圆柱体体系

一个简单的建模脚本,能够将线性分子堆积成圆柱体体系,用于模拟。可调参数输入-h查看,参考了构建表面修饰长链分子的纳米颗粒模型一文中的函数代码。代码见链接,例子下载例子

简单突变脚本

一个简单的突变蛋白单个氨基酸为丙氨酸的脚本,参考GMXPBSA 2.1: A GROMACS tool to perform MM/PBSA and computational alanine scanning一文源码。代码见链接

分子线性拼接脚本

参看原始博文拼接分子的简单脚本,修改成了读取gro结构文件(gromacs常用)。代码见链接

拉氏图转彩色脚本

gmx rama -o生成的rama.xvg文件转为彩色rama.xpm图片。代码见链接.

熵计算脚本

简单的脚本,通过协方差矩阵的特征值(gmx covar -o)产生的eigenval.xvg来计算熵值。采用quasi-Harmonic approximation,原始文献参看J Chem Phys 115:6289等式13。代码见链接.
相互作用熵脚本,原始文献参看Interaction Entropy: A New Paradigm for Highly Efficient and Reliable Computation of Protein–Ligand Binding Free Energy

能量分解脚本

结合g_mmpbsa使用的脚本,思路源于g_mmpbsa脚本,此处是Python脚本,可能用起来不太方便,我将它改写成了awk脚本。代码见链接1,链接2

氢键分析

思路来源于readHBmap.tar.gz,此处是个Python脚本,我将它改写成了bash脚本。代码见链接

构建两水球碰撞体系

构建两个水球的碰撞体系,并更改两水球原子沿X轴速度大小相等方向相反。原始博文参看GROMACS模拟实例–两水球碰撞。下载链接

轨迹pdb或者gro文件拆分为单个pdb或者gro结构

代码见链接

二分法求解单调函数零点

估计区域[a, b],然后利用二分求解方程在该区域的零点值,代码见链接

Levenberg-Marquardt算法求解拟合参数

参照别人的Python代码,硬着头皮写完了,awk自定义矩阵函数属实恶心,看来要转战Py了^-^。有问题,但是大体是这样,代码看链接

在球面上生成均匀分布的数

参考文章见Generating uniformly distributed numbers on a sphereawk代码见链接

温度对二元铜锌合金系统属性影响

采用二分法,有限差分方式确定温度对合金体系不同属性的影响。awk代码看链接

使用速度Verlet算法模拟铝原子体系。

改自C-Example2:使用速度Verlet算法计算铝原子体系能量一文。采用awk编程进行小规模模拟,速度的确和C或者Fortran等专业的编程差距太大。仅供学习之用。代码较长,就不贴出来了,看链接

高斯-约旦(Gauss-Jordan)消元法求N阶逆矩阵

最近研究了一下Levenberg-Marquard算法来求解非线性方程参数,但是基本上能找到的要么是C或者Fortran程序,要么是PythonMatlab。本来想参考这些现成的来写个awk的脚本,无奈awk支持的内置函数实在是太少了。特别是矩阵运算这一块完全都需要自己编写函数才行,太麻烦了。其中有个求逆矩阵的函数,挺有趣,记之。看代码

线性、二次方程及Freundlich等温吸附拟合

代码

Vasp轨迹坐标(XDATCAR)文件转换为.gro格式轨迹文件

原文参看博文新版VASP轨迹文件转换为GROMACS格式。但感觉有些毛病,稍微修改了一下以供使用。看代码

F3/AOP序参数计算实例

看过水分子序参数的计算一博文的都知道通过调用GROMACS自带分析工具计算虽然方便,但是效率并不高。前面的相关博文中我也给出了更具效率的C程序。折中一下,你可以直接通过bash独自去处理gro文件,不调用GROMACS的内置工具,这样效率也会大幅度提高。注意:提取gro中的数据不要用域来判断,用substr函数来得到固定位置的字符或者值。比如代码

bash下的MD程序?

尝试用bash下的awk编程跑MD程序?代码改自C-Example1:铝原子体系势能与体积的变化关系一文。看代码

MS轨迹坐标提取

代码.

格栅参数计算

简单的格栅参数计算脚本,只需更改部分初始参数,即可一次性得到所求参数。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
#!/bin/bash
# A simple bash script for grating design
# Author: liuyujie

Q=10 # m^3/s 流量
B1=5 # 栅前渠宽,一般由经验公式:B1=sqrt(2*Q/v)
S=0.01 # m 栅条宽度
g=9.81 # m/(s^2) 重力加速度
Ktotal=1.5 # 总变化系数
h=0.4 # m 栅前水深
h1=0.3 # m 栅前渠道超高,一般取0.3m
alpha1=20 # degree 进水渠道展开角度,一般取20度
v=0.9 # m/s 过栅流速,一般取0.8~1.0 (m^2/s)
b=0.02 # m 粗格栅50-100mm,中格栅10-40mm,细格栅3-10mm
alpha=60 # degree 格栅倾角,机械清除60~90度,人工清除为30~60度
k=3 # 格栅受污染堵塞后,水头损失增大数倍,一般取k=3
beta=2.42 # 栅条形状系数,用于求阻力系数,锐边矩形=2.42,迎水面为半圆形的
# 矩形=1.83,圆形=1.79,迎水面和背水面均为半圆形的矩形=1.67,
# 梯形=2.00,正方形=0.64(收缩系数)
W1=0.07 # m^3/(1000m^3) 单位体积污水栅渣量,一般取0.1~0.01,粗格栅取小值,细格栅取大值

awk -v Q=$Q -v S=$S -v g=$g -v Ktotal=$Ktotal -v h=$h -v h1=$h1 -v alpha1=$alpha1 \
-v v=$v -v b=$b -v alpha=$alpha -v k=$k -v beta=$beta -v W1=$W1 -v B1=$B1 \
-v PI=3.141592653 '
BEGIN {
n = Q * sqrt(sin(alpha*PI/180.0))/(b * h * v); n = round(n); B = S*(n-1)+b*n
printf("格栅间隙数n=%d(个) 格栅槽宽度B=%.1f(m)\n", n, B)
if(beta!=0.64) {
Rc = beta * (S/b)^(4/3)
}
else {
Rc = ((b+S)/(beta*b)-1)^2
}
h0 = Rc * v^2/(2*g) * sin(alpha*PI/180.0)
h2 = k * h0
printf("计算水头损失h0=%.6f(m) 过栅水头损失h2=%.6f(m)\n", h0, h2)
printf("栅后槽总高度H=%.6f(m)\n", h+h1+h2)
L1 = (B-B1)/(2*sin(alpha1*PI/180.0)/cos(alpha1*PI/180.0))
L2 = 0.5*L1
H1 = h + h1
L = L1 + L2 + 0.5 + 1.0 + H1/(sin(alpha*PI/180.0)/cos(alpha*PI/180.0))
printf("进水渠道肩宽L1=%.6f(m) 渐宽部分长度L2=%.6f(m)\n栅前槽高H1=%.1f(m) " \
"栅槽总长度L=%.6f(m)\n", L1, L2, H1, L)
W = Q*W1*86400/(Ktotal*1000)
printf("每日栅渣量W=%.2f(m^3/d)\n", W)
}

function round(kl) {
b5 = int(kl)
cc = kl * 10
df = b5 * 10 + 5
if(cc >= df) {
e4 = b5 + 1
}
else {
e4 = b5
}
return e4
}
'