- 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.258330 2.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.392305 0.000000
H -9.000000 -8.660254 2.449490
S -9.000000 -10.392305 4.898979
S -9.000000 -8.660254 7.348469
H -10.500000 -6.062178 -7.348469
C -10.500000 -7.794229 -4.898979
.....(省略)
S 9.000000 10.392305 -4.898979
C 9.000000 10.392305 0.000000
S 9.000000 10.392305 4.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.79053824 1.01170907 -0.03365831
C 9.28080412 1.70826814 1.24932447
C 8.78721231 3.16696710 1.26214971
H 10.35034747 1.69339092 1.27681426
H 8.89485019 1.19360662 2.10434675
C 9.27747818 3.86352617 2.54513250
H 9.17316623 3.68162862 0.40712744
H 7.71766896 3.18184432 1.23465992
C 8.78388637 5.32222513 2.55795773
H 8.89152426 3.34886465 3.40015477
H 10.34702153 3.84864895 2.57262229
H 9.16984029 5.83688665 1.70293546
H 7.71434302 5.33710235 2.53046794
C 9.27415225 6.01878420 3.84094052
C 8.78056044 7.47748316 3.85376576
H 10.34369560 6.00390698 3.86843031
H 8.88819833 5.50412268 4.69596279
C 9.27082631 8.17404223 5.13674854
H 9.16651436 7.99214468 2.99874348
H 7.71101709 7.49236038 3.82627597
H 8.88487239 7.65938071 5.99177081
H 10.34036966 8.15916501 5.16423833
C 8.77723450 9.63274119 5.14957378
O 8.60383009 10.24786466 4.06554214
O 8.51594972 10.28824341 6.39333706
H 7.82155932 10.93991893 6.27192456
|
我们可以利用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]=0
if(Norm=="X") Vatm[1]=1
else if(Norm=="Y") Vatm[2]=1
else if(Norm=="Z") Vatm[3]=1
}
{ NatmA=$1
getline
Nsrf=0; Xcom=0; Ycom=0; Zcom=0
for(i=1; i<=NatmA; i++) {
getline
SatmA[i]=$1; XatmA[i]=$2; YatmA[i]=$3; ZatmA[i]=$4
Xcom += $2; Ycom += $3; Zcom += $4
if(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.9f\n", 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 #%iMolB\n", SatmB[j], \
V[1]+XatmA[Itmp], V[2]+YatmA[Itmp], V[3]+ZatmA[Itmp], i
}
}
}
function Acos(x) { return atan2(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%的表面原子.
将得到的构型稍加修改, 就可以用于模拟了.
代码下载
点击这里下载脚本及示例文件及代码.
网络资料
评论
- 2016-03-10 08:19:06
刘小康好厉害！！！ -
2016-03-11 15:21:43
Jerkwin谢谢鼓励. - 2016-03-11 21:27:51
赞赞 -
2016-03-14 23:30:35
Jerkwin谢. - 2016-07-26 09:48:24
小和尚老师您好，我需要对表面进行部分接枝，能否实现随机选取表面原子进行部分接枝？ -
2016-07-26 10:35:39
Jerkwin可以, 只要指定相应的原子即可. 至于如何随机, 那要看你的定义. - 2016-08-03 13:13:45
张蕴瀚老师您好，我需要对一个金表面进行修饰巯基分子，但只需要一层金原子，这样是不是没了法向？这个脚本应该如何修改？ -
2016-08-04 10:35:09
Jerkwin表面的话更简单, 手动做也不麻烦. 使用这个脚本的话, 你要修改法向的代码, 使其指向表面的法向. - 2016-11-07 18:23:35
franklinly这个很有用，谢谢李老师。
