博文
构建表面修饰长链分子的纳米颗粒模型
|
2016-02-21 00:01:09 初稿
2016-08-16 22:58:21 修订, 增加随机选择, 指定法向功能
在进行与纳米颗粒有关的分子动力学模拟时, 有时候需要对纳米颗粒的表面进行修饰, 例如, 将一些长链的表面活性剂分子连接到颗粒的表面, 以改变纳米颗粒的亲水性. 这种模型的构建牵涉到两个问题, 一个是如何得到纳米颗粒的结构, 另一个是有了纳米颗粒的结构后, 如何将一些长链分子连接到颗粒的表面.
构建球形纳米颗粒一般而言, 我们可以根据原子的晶胞结构来构建相应的晶体, 然后根据一定的条件去除不需要的原子, 就得到了一定形状的纳米颗粒. 对于球形的纳米颗粒, 这可以利用我编写的一个小工具来完成. 只要选择好晶胞类型(一般使用FCC或HCP, 因为这两种是最密堆积), 团簇类型球体, 设定好晶胞参数和球体半径, 点击创建就可以得到纳米颗粒的坐标文件了. 在坐标文件中, 正常位点以C表示, 面心位点以F表示, 六方位点以H表示, 表面原子则以S表示. 程序中表面原子的确定方法比较粗略, 你可以根据自己的需要自行决定将哪些原子指定为表面原子, 只要将相应的元素符号改为S即可.
作为示例, 我使用HCP晶胞类型, 其他参数默认, 得到的结构文件如下(中间部分省略):
| HCP.xyz | |
|---|---|
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17 | 751
HCP Sphere Radius=5
H -7.500000-11.258330-2.449490
H -7.500000-11.2583302.449490
S -9.000000-8.660254-7.348469
S -9.000000-10.392305-4.898979
H -9.000000-8.660254-2.449490
C -9.000000-10.3923050.000000
H -9.000000-8.6602542.449490
S -9.000000-10.3923054.898979
S -9.000000-8.6602547.348469
H -10.500000-6.062178-7.348469
C -10.500000-7.794229-4.898979
.....(省略)
S 9.00000010.392305-4.898979
C 9.00000010.3923050.000000
S 9.00000010.3923054.898979 |
为避免后面添加长链分子后, 分子其中的氢原子与颗粒中的H位点相混淆, 我们先将H原子类型替换为C. 如果颗粒的原子类型不是C的话, 你可能也需要将C改为你需要的原子类型, 以免与长链分子中的碳原子混淆.
将这个文件保存为HCP.xyz备用.
构建长链分子长链分子的构建可采用你熟悉的分子编辑软件, 只要能输出xyz格式的结构文件即可.
我们以下面的分子作为示例:
| C9.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
27
28 | 26
C9COOH
C 8.790538241.01170907-0.03365831
C 9.280804121.708268141.24932447
C 8.787212313.166967101.26214971
H 10.350347471.693390921.27681426
H 8.894850191.193606622.10434675
C 9.277478183.863526172.54513250
H 9.173166233.681628620.40712744
H 7.717668963.181844321.23465992
C 8.783886375.322225132.55795773
H 8.891524263.348864653.40015477
H 10.347021533.848648952.57262229
H 9.169840295.836886651.70293546
H 7.714343025.337102352.53046794
C 9.274152256.018784203.84094052
C 8.780560447.477483163.85376576
H 10.343695606.003906983.86843031
H 8.888198335.504122684.69596279
C 9.270826318.174042235.13674854
H 9.166514367.992144682.99874348
H 7.711017097.492360383.82627597
H 8.884872397.659380715.99177081
H 10.340369668.159165015.16423833
C 8.777234509.632741195.14957378
O 8.6038300910.247864664.06554214
O 8.5159497210.288243416.39333706
H 7.8215593210.939918936.27192456 |
Fig.1
我们可以利用C1和C23来指定这个分子的长轴方向. 由于C1会连接到纳米颗粒表面, 所以需要将本来与它相连的氢原子删除.
将结构保存为C9.xyz备用.
将长链分子连接到纳米颗粒表面在进行表面修饰时, 我们要在指定为S的表面原子上连接长链分子. 连接时, 要使长链分子的长轴方向沿着表面的法向(对球体其法向即径向), 这样各个分子之间的重叠和交叉才最小. 手动进行这些操作也可以, 但当分子数较多时就变得很麻烦了. 因此, 我写了一个简单的bash脚本apdMol来完成这件事. 这个脚本在运行时需要两个结构文件, 一个是纳米颗粒的结构文件, 一个是长链分子的结构文件, 此外, 还需要给出长链分子中的两个原子编号来指定长链分子的长轴方向. 有了这些信息后, 脚本会自动将纳米颗粒中标为S的表面原子替换为长链分子, 并保证长链分子的长轴方向沿颗粒的法向.
原始脚本发布后, 有人希望能增加些功能
随机选取部分表面原子进行修饰
能够使用其他法向, 这样就能用于修饰表面
第二个功能实现起来并不复杂, 但对于第一个功能, 必须有些考虑, 想明白这里所说的”随机”是什么意思. 对于球形颗粒来说, 其表面随机分布的点看起来可能并不是均匀的, 有些地方疏, 有些地方密, 这可能并不是你需要的结果. 此外, 有些人所说的随机很可能指的是随机均匀分布, 需要把一定数量的点均匀分布在表面上, 使点与点之间的间距尽量相等且最大. 这两种随机方法实现起来都不简单, 就不详细说明了. 目前代码中实现的随机是最简单随机方法, 就是随机从表面点中选取一定数目的点. 严格说来这是一种很差的随机方法(甚至可说是错误的), 但容易实现, 对一般应用来说也不会引起多大问题. 在使用代码时, 请不要忘记这一点.
| apdMol.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
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144 | usage=">>>>>>>>>>>>>>>> apdMol <<<<<<<<<<<<<<<<>>>>>>>>>>>>>>>> Jerkwin@gmail.com <<<<<<<<<<<<<<<<>>>>>>>>>> 2016-02-20 23:20:36 初版 <<<<<<<<<>>>>>>>>>> 2016-08-16 22:58:21 修订 <<<<<<<<<>> 用法: bash $0 <MolA.xyz> <MolB.xyz> {Iatm} {Jatm} {Lsrf} {Apd#%} {Norm}>> 默认: bash $0 <MolA.xyz> <MolB.xyz> 1 2 S 50% COG>> MolA: 待修饰的团簇或表面>> MolB: 修饰到MolA的分子>> Iatm: MolB的#Iatm原子将重合到MolA中的某个表面原子>> Jatm: #Iatm->#Jatm原子定义MolB的主轴方向>> Lsrf: MolA中表面原子的标识符号>> Apd#%: MolA中要修饰的表面原子的个数或百分比>> Norm: MolA表面原子的法向 对球形团簇可使用COG(到几何中心法向) 对表面可使用X|Y|Z 也可直接指定{x,y,z}(尚未实现)"
[[ $# -lt 2 ]] && { echo "$usage"; exit; }
MolA=$1; MolB=$2
Iatm=1; [[ x$3!= x ]] && { Iatm=$3; }
Jatm=2; [[ x$4!= x ]] && { Jatm=$4; }
Lsrf="S"; [[ x$5!= x ]] && { Lsrf=$5; }
Napd="50%"; [[ x$6!= x ]] && { Napd=$6; }
Norm="COG"; [[ x$7!= x ]] && { Norm=$7; }
awk -v file=$MolB -v Iatm=$Iatm -v Jatm=$Jatm \-v Lsrf=$Lsrf -v Napd=$Napd -v Norm=$Norm 'BEGIN{ Pi=4*atan2(1,1)
getline NatmB < file
getline< file
for(i=1; i<=NatmB; i++) {
getline< file
SatmB[i]=$1; XatmB[i]=$2; YatmB[i]=$3; ZatmB[i]=$4
}
# 使Iatm置于原点
x=XatmB[Iatm]; y=YatmB[Iatm]; z=ZatmB[Iatm]
for(i=1; i<=NatmB; i++) {
XatmB[i] -= x; YatmB[i] -= y; ZatmB[i] -= z
}
# 获取分子主轴, 归一化
Vij[1]=XatmB[Jatm]
Vij[2]=YatmB[Jatm]
Vij[3]=ZatmB[Jatm]
UniVect(Vij)
# 表面原子法向
Vatm[1]=0; Vatm[2]=0; Vatm[3]=0if(Norm=="X") Vatm[1]=1elseif(Norm=="Y") Vatm[2]=1elseif(Norm=="Z") Vatm[3]=1
}
{ NatmA=$1getline
Nsrf=0; Xcom=0; Ycom=0; Zcom=0for(i=1; i<=NatmA; i++) {
getline
SatmA[i]=$1; XatmA[i]=$2; YatmA[i]=$3; ZatmA[i]=$4
Xcom +=$2; Ycom +=$3; Zcom +=$4if(SatmA[i]==Lsrf) { Nsrf++; Isrf[Nsrf]=i }
}
Xcom /= NatmA; Ycom /= NatmA; Zcom /= NatmA
# 随机混洗
FisherYates(Isrf, Nsrf)
if(index(Napd, "%")) { sub("%","",Napd); Napd=int(Nsrf*Napd/100) }
if(Napd>Nsrf) Napd=Nsrf
print NatmA+(NatmB-1)*Napd
print Napd"/"Nsrf " Surface Atoms Appended with Another Molecule"for(i=1; i<=NatmA; i++)
printf"%4s %15.9f %15.9f %15.9fn", SatmA[i], XatmA[i], YatmA[i], ZatmA[i]
for(i=1; i<=Napd; i++) {
Itmp=Isrf[i]
if(Norm=="COG") { # 待转至的轴, 相对于分子几何中心的径向
Vatm[1]=XatmA[Itmp]-Xcom
Vatm[2]=YatmA[Itmp]-Ycom
Vatm[3]=ZatmA[Itmp]-Zcom
}
UniVect(Vatm)
# 绕两轴中线旋转180度
x=(Vatm[1]+Vij[1])/2
y=(Vatm[2]+Vij[2])/2
z=(Vatm[3]+Vij[3])/2
AxiRotMat(Pi, x, y, z)
# 追加分子for(j=2; j<=NatmB; j++) {
V[1] = XatmB[j]
V[2] = YatmB[j]
V[3] = ZatmB[j]
RotVect(RotMat, V)
printf"%4s %15.9f %15.9f %15.9f #%iMolBn", SatmB[j],
V[1]+XatmA[Itmp], V[2]+YatmA[Itmp], V[3]+ZatmA[Itmp], i
}
}
}
function Acos(x) { returnatan2(sqrt(1-x*x),x) }
function Dot(A, B) { return A[1]*B[1]+A[2]*B[2]+A[3]*B[3] }
function VectAng(A, B) { return Acos(Dot(A, B)/sqrt(Dot(A, A)*Dot(B, B))) }
function UniVect(V, Rtmp) {
Rtmp=Dot(V, V)
if (Rtmp!=0) {
Rtmp=1/sqrt(Rtmp)
V[1] *= Rtmp; V[2] *= Rtmp; V[3] *= Rtmp
}
}
function RotVect(RotMat, V, Vx,Vy,Vz) { # 按照矩阵RotMat旋转矢量V
Vx = V[1]; Vy = V[2]; Vz = V[3]
V[1]= RotMat[1,1]*Vx + RotMat[1,2]*Vy + RotMat[1,3]*Vz
V[2]= RotMat[2,1]*Vx + RotMat[2,2]*Vy + RotMat[2,3]*Vz
V[3]= RotMat[3,1]*Vx + RotMat[3,2]*Vy + RotMat[3,3]*Vz
}
function AxiRotMat(tht, x, y, z, Rtmp, w) { # 绕矢量(x,y,z)旋转tht弧度, 四元数方法
Rtmp=x*x+y*y+z*z
if(Rtmp>1E-6) { Rtmp=1/sqrt(Rtmp); x *= Rtmp; y *= Rtmp; z *= Rtmp }
Rtmp=sin(tht/2); w=cos(tht/2); x *= Rtmp; y *= Rtmp; z *= Rtmp
RotMat[1,1]=1-2*(y*y+z*z); RotMat[1,2]=2*(x*y-w*z); RotMat[1,3]=2*(x*z+w*y)
RotMat[2,1]=2*(x*y+w*z); RotMat[2,2]=1-2*(x*x+z*z); RotMat[2,3]=2*(y*z-w*x)
RotMat[3,1]=2*(x*z-w*y); RotMat[3,2]=2*(y*z+w*x); RotMat[3,3]=1-2*(x*x+y*y)
}
function FisherYates(arr, len, i,j,tmp) {
srand()
for(i=len; i>1; i--) {
j =int(rand()*i)+1
tmp = arr[i]
arr[i] = arr[j]
arr[j] = tmp
}
}
'$MolA >${MolA%.*}+${MolB%.*}.xyz
|
使用bash apdMol.bsh HCP.xyz C9.xyz 1 23运行脚本, 会得到HCP+C9.xyz文件. 默认修饰了50%的表面原子.
Fig.2
将得到的构型稍加修改, 就可以用于模拟了.
代码下载点击这里下载脚本及示例文件及代码.
◆本文地址: http://jerkwin.github.io/2016/02/21/构建表面修饰长链分子的纳米颗粒模型/, 转载请注明◆
https://blog.sciencenet.cn/blog-548663-957611.html
上一篇：GaussView使用cube文件作图时的成键问题
下一篇：X射线安检仪图像的颜色映射
