拼接分子的简单脚本

类别:    标签: awk 编程   阅读次数:   版权: (CC) BY-NC-SA

在分子建模, 特别是聚合物建模过程中, 有时我们需要构建具有周期性重复结构的分子, 或者是具有不同聚合度的分子. 虽然有些可视化的建模软件支持这种功能, 但使用起来不一定方便, 所以我还是写了一个简单的bash脚本来做这件事情, 放在这里备份下, 也供需要的人参考.

脚本主要使用awk语言, 利用分子单体构建聚合物. 方法很简单, 先在分子单体中定义某个原子为”头原子”, 另一原子为”尾原子”, 拼接时将第二个分子单体平移至其头原子与第一个分子单体的尾原子重合的位置, 然后将这两个原子删除. 当然, 首端的分子单体会保留其头原子, 而末端的分子单体会保留其尾原子.

使用时, 在命令行中指定分子单体的文件名称, 头原子编号, 尾原子编号, 拼接的片段数目(聚合度)即可.

为方便读取和书写, 输入文件与输出文件均采用.xyz格式: 第一行指定原子数目, 第二行是随意的标题行, 后续每行指定每个原子的元素符号与xyz坐标. 值得注意的是, 使用VMD查看分子结构时, 其中的原子编号是从0开始, 所以VMD提示编号为n的原子在从1开始的编号中为n+1号原子.

由于采用的是头原子与尾原子重合并删除的拼接方法, 所以默认情况下, 删除重合原子后, 其相邻原子间的距离可能过大, 导致可视化软件判断成键时不正常, 简单的解决方法是将头原子与尾原子所在键的长度适当缩小. 虽然这样拼接后所得聚合物首端头原子和末端尾原子的成键距离偏小, 但用于分子动力学模拟的初始构型足够了.

脚本

commol.bsh
 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
usage="\
>>>>>>>>>>>>>>>>     comMol    <<<<<<<<<<<<<<<<
>>>>>>>>>>>>>>>>    Jicun Li   <<<<<<<<<<<<<<<<
>>>>>>>>>>     2017-08-27 09:02:30     <<<<<<<<
>>Usage:  commol File.xyz #head #tail #mol"

[[ $# -lt 1 ]] && { echo "$usage"; exit; }

fxyz=$1 # 分子单体xyz文件名称
head=$2 # 头原子编号, 从1开始
tail=$3 # 尾原子编号, 从1开始
n=$4    # 拼接数目(聚合度)

fxyz=${fxyz%.xyz}.xyz    # 去除输入文件扩展名
fout=${fxyz%.xyz}~$n.xyz # 输出文件名称

awk -v head=$head -v tail=$tail -v n=$n '
NF>0 {	Natm=$1 # 原子数目
	getline Tips
	for(i=1; i<=Natm; i++) {
		getline # 单体中每个原子的原子类型, xyz坐标
		Satm[i]=$1; Xatm[i]=$2; Yatm[i]=$3; Zatm[i]=$4
	}

	for(i=1; i<=n; i++) {        # 片段数目
		for(j=1; j<=Natm; j++) { # 原子编号
			s[i,j]=Satm[j]       # i片段中j原子的原子类型
			x[i,j]=Xatm[j]       # 和初始坐标与单体相同
			y[i,j]=Yatm[j]
			z[i,j]=Zatm[j]
		}
	}
	print (Natm-2)*n+2       # 拼接后聚合物的原子数目
	print Tips, head, tail, n # 拼接信息

	for(i=1; i<=n; i++) {
		if(i>1) { # 非首端分子单体平移距离
			dx=x[i,head]-x[i-1,tail]
			dy=y[i,head]-y[i-1,tail]
			dz=z[i,head]-z[i-1,tail]
		}
		for(j=1; j<=Natm; j++) {
			x[i,j] -= dx
			y[i,j] -= dy
			z[i,j] -= dz
			# 首端分子单体不输出尾原子
			# 末端分子单体不输出头原子
			# 中间分子单体不输出头尾原子
			if( (i==1 && j!=tail) \
			  ||(i==n && j!=head) \
			  ||(i!=1 && i!=n && j!=head && j!=tail) ) \
				printf "%3s %12.6f%12.6f%12.6f\n", s[i,j], x[i,j], y[i,j], z[i,j]
		}
	}
}
' $fxyz > $fout

示例

以下面的分子为例, 头原子为3号, 尾原子为11号, 为清楚起见原子类型已改为Cl, 且它们所在键长已经适当调整.

ht.xyz
 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
24
Tips
 C   -2.25735638   -0.11936983   -0.24914958
 C   -0.73984978   -0.32271456   -0.08354569
 Cl  -2.37768477    0.47589194   -0.30646016
 H   -2.59326436   -0.63693340   -1.12333136
 H   -2.76475201   -0.50466899    0.61049816
 C   -0.43606989   -1.82550657    0.06114004
 H   -0.23245414    0.06258460   -0.94319342
 H   -0.40394179    0.19484900    0.79063609
 C    1.08143670   -2.02885130    0.22674393
 H   -0.94346553   -2.21080573    0.92078778
 Cl   1.21557328   -2.69242179    0.29063114
 H    1.58883234   -1.64355214   -0.63290380
 H    1.41734469   -1.51128773    1.10092571
 C   -0.91952625   -2.57041114   -1.19702812
 C   -1.34532335   -3.90221323   -1.10257193
 C   -0.93367443   -1.91647220   -2.43641735
 C   -1.78527002   -4.58007592   -2.24750470
 H   -1.33451760   -4.40151122   -0.15627137
 C   -1.37362358   -2.59433407   -3.58134966
 H   -0.60856921   -0.89961164   -2.50853669
 C   -1.79942144   -3.92613591   -3.48689333
 H   -2.11037260   -5.59693736   -2.17538586
 H   -1.38442919   -2.09503612   -4.52765025
 H   -2.13533239   -4.44369848   -4.36107455

视图: 投影 正交    速度:
模型: 球棍 范德华球 棍状 线框 线型   名称
左键: 转动   滚轮: 缩放   双击: 自动旋转开关   Alt+左键: 移动

Fig.1

执行下面的命令获得聚合度为10的分子

bash commol.bsh ht.xyz 3 11 10

得到聚合物结构如下


视图: 投影 正交    速度:
模型: 球棍 范德华球 棍状 线框 线型   名称
左键: 转动   滚轮: 缩放   双击: 自动旋转开关   Alt+左键: 移动

Fig.2

待完善

随意赞赏

微信

支付宝
◆本文地址: , 转载请注明◆
◆评论问题: https://jerkwin.herokuapp.com/category/3/博客, 欢迎留言◆


前一篇: 使用GROMACS计算红外光谱
后一篇: GROMACS非标准残基教程1:修改力场与增加残基

访问人次(2015年7月 9日起): | 最后更新: 2020-07-20 13:23:09 UTC | 版权所有 © 2008 - 2020 Jerkwin