博文
拼接分子的简单脚本
|
2017-08-27 05:51:29 感谢 李传玺 整理说明
在分子建模, 特别是聚合物建模过程中, 有时我们需要构建具有周期性重复结构的分子, 或者是具有不同聚合度的分子. 虽然有些可视化的建模软件支持这种功能, 但使用起来不一定方便, 所以我还是写了一个简单的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.6fn", 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.377684770.47589194-0.30646016
H -2.59326436-0.63693340-1.12333136
H -2.76475201-0.504668990.61049816
C -0.43606989-1.825506570.06114004
H -0.232454140.06258460-0.94319342
H -0.403941790.194849000.79063609
C 1.08143670-2.028851300.22674393
H -0.94346553-2.210805730.92078778
Cl 1.21557328-2.692421790.29063114
H 1.58883234-1.64355214-0.63290380
H 1.41734469-1.511287731.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 |
Fig.1
执行下面的命令获得聚合度为10的分子
得到聚合物结构如下
Fig.2待完善
输入与输出最好使用pdb格式, 但由于pdb文件格式严格, 且版本众多, 所以暂时使用xyz格式
自动根据成键信息缩短头尾原子距离, 但代码会变得过于复杂, 暂时没有必要
做成在线小工具, 类似 构建表面修饰长链分子的纳米颗粒模型
拼接时可随机旋转分子, 避免线性结构
◆本文地址: https://jerkwin.github.io/2017/08/27/拼接分子的简单脚本/, 转载请注明◆
◆评论问题: https://jerkwin.herokuapp.com/category/3/博客, 欢迎留言◆
https://blog.sciencenet.cn/blog-548663-1072914.html
上一篇：使用GROMACS计算红外光谱
下一篇：GROMACS非标准残基教程1：修改力场与增加残基
