Packmol的网址:
输入文件里的关键词在:
全局设定
tolerance [数值]:要求原子间距离不能小于多少,一般设为2.0。Packmol里所有长度的单位都是埃。
output [文件名]:输出文件的路径
filetype [类型]:输出文件的格式。默认为pdb,也可以设比如xyz
seed -1:写这个的话每次运行产生的结构都会不同
add_box_sides:给输出的pdb中增加CRYST1字段定义盒子,盒子尺寸对应体系中原子最大、最小坐标。如果写诸如add_box_sides 1.2,则会给盒子各方向再增加1.2埃,避免做模拟时有某些边界的原子和镜像盒子里的原子间存在不合理接触
分子设定
每个下面这种字段设置一种分子如何出现在当前体系中。可以写无数个这种字段。
structure [分子结构文件]number [分子数][分子规则]end structure
分子规则用来定义分子里的所有原子必须满足的条件,其常用关键词如下。
inside cube [xmin ymin zmin d]:要求分子出现在长度为d的立方盒子内,盒子x,y,z最小值为xmin ymin zmin
inside box [xmin ymin zmin xmax ymax zmax]:要求分子出现在矩形盒子内,盒子两个顶角坐标分别为xmin ymin zmin和xmax ymax zmax
inside sphere [x y z r]:要求分子出现在中心为(x,y,z)半径为r的圆球中。另外,用inside ellipsoid关键词的话还可以要求出现在特定的椭球中
inside cylinder [a1 b1 c1 a2 b2 c2 r l]:要求分子出现在圆柱内。圆柱的定义是从(a1,b1,c1)往(a2,b2,c2)方向伸展l长度,半径为r
以上inside都可以改为outside,要求不能出现在指定范围内。上面的空间范围要求可以同时使用多个,条件会同时满足。还有一些设置:
over plane [a b c d]:要求分子出现位置满足平面方程ax+by+cz-d>=0。比如over plane 1 0 0 10.5就相当于要求出现在x>=10.5埃的区域。over可以改成below来反转条件。
constrain_rotation [x/y/z 平均值 最大偏差]:默认情况下分子可以绕x,y,z轴随意旋转任意角度。用这个选项可以设置分子绕x或y或z旋转的情况,可以同时设多个。如constrain_rotation x 20 5代表允许绕着y和z轴随意旋转,而绕着x轴旋转的范围必须在15~25度之间,平均值为20度。
fixed [x y z a b c]:用来直接定义分子出现在哪里。此设置代表将分子结构文件里的坐标平移(x,y,z),并绕x,y,z轴分布旋转a,b,c弧度。如果还同时写了center关键词,相当于把分子的中心摆到(x,y,z)位置再旋转
原子设定
Packmol里还可以要求某种分子的某些特定序号的原子出现的范围必须满足的要求,格式为:
structure [分子结构文件]
[分子规则]
atoms [第一批原子序号]
[原子规则]
end atoms
atoms [第二批原子序号]
[原子规则]
end atoms
...
end structure
原子规则的设置语法、关键词和分子规则完全一样。
实例
官网:
构建5*5*5 nm^3甲醇盒子
首先在GaussView或其它可以对分子建模的程序里画一个甲醇,然后保存为methanol.pdb。对于这种刚性小分子,用Packmol之前结构不需要非得优化,毕竟GaussView直接画的甲醇的基本几何参数是合理的,而且不管对其做不做优化,等动力学跑起来之后就都一样了。
创建一个名为methanol_box.inp的文本文件,内容如下
tolerance 2.0
add_box_sides 1.2
output methanol_box.pdb
structure methanol.pdb
number 1000
inside cube 0. 0. 0. 50.
end structure
将methanol.pdb和methanol_box.inp都放到当前目录,运行packmol < methanol_box.inp,经过几秒钟,当前目录下就出现了methanol_box.pdb。将之拖到VMD程序里观看,并且在VMD的命令行窗口里输入pbc box显示出盒子,看到的图像
构建水+甲醇混合溶液的气液界面体系
此例我们想构建一个水+甲醇摩尔比1:1溶液的气液界面体系,界面平行于XY方向,盒子的X、Y长度都为4nm,溶液层厚度是5nm,溶液上下两侧都是真空区,每一侧真空区厚度都是3nm。
首先,我们靠Packmol创建一个水+甲醇的4*4*5 nm^3的溶液盒子。甲醇的pdb文件已经有了,我们再用GaussView画个水,保存为water.pdb。然后编辑一个文件mix.inp,内容如下(经测试,水和甲醇各插入1000个时盒子已经比较致密了,而且已经需要好多轮循环才能成功产生结构,所以笔者就不再尝试插入更多了)
tolerance 2.0
add_box_sides 1.2
output mix.pdb
structure water.pdb
number 1000
inside box 0. 0. 0. 40. 40. 50.
end structure
structure methanol.pdb
number 1000
inside box 0. 0. 0. 40. 40. 50.
end structure
将得到的mix.pdb弄到VMD里显示
相关问题&注意事项
这里说一些使用Packmol建模应当注意的问题,并且结合最主流的分子动力学程序GROMACS说一些相关事项。
Packmol的执行过程实际上是个循环过程,如果你当前的空间范围设定比较简单,Packmol很快就会顺利产生最终结构。空间约束设定得越多、越复杂,往往成功产生结构需要的循环次数越多,总耗时越高。如果你把空间范围约束得过度复杂甚至相互矛盾,或者要加入的分子数较多但允许分子出现的空间范围相对而言过度狭小,那么Packmol会反复循环,一直也收敛不了,或者最终宣告产生失败。碰到这种情况,应当将约束条件适当放宽、简化,或者设的分子数少一点再试。
一般情况下,用Packmol构建的体系肯定是比较松散的,或者说其密度肯定显著低于真实密度,因为此时还没有考虑分子间相互作用、构象的变化从而使得分子间能够接触得更紧密。但这没关系,在NPT下跑一下MD就好了,盒子尺寸会自发进行调节。而直接用Packmol产生的结构做NVT模拟一般是无意义的,因为跑起来之后分子会自发紧密聚集,其它地方就成了真空区了,和实际不符。
用Packmol实现构建含有Na+、Cl-离子的溶液体系是可以的。但是肯定没有用GROMACS里专用的gmx genion命令好,因为gmx genion在插入时会考虑体系电荷分布,会把阳离子和阴离子分别加入到静电势较负和较正的地方,这样比较合理。而Packmol则没有考虑这点,并且有可能比如产生的结构里恰好Na+出现在显著带正电荷的区域,这可能会令模拟初期稳定性较差(此时需要用较谨慎的模拟参数,比如小步长)。
虽然也可以用Packmol构建蛋白质、核酸浸在溶剂环境中的体系,但是这样做明显不如用动力学程序自带的专用做这种事情的程序好,因为Packmol产生的水的密度偏低,水的分布特征和实际体相水相差较大,可能NPT模拟起来之后盒子变形、收缩得厉害,出现溶质与其镜像最近距离太近之类问题。而如果用比如GROMACS里的gmx solvate命令给蛋白质/核酸加溶剂,由于是直接用事先已经做NPT平衡好的溶剂盒子(比如加水一般用自带的spc216.gro)通过平移复制来填充真空区,加的溶剂的分布状态明显理想得多得多。
GROMACS有一个命令叫gmx insert-molecules,可以由用户提供分子结构文件、指定插入盒子的分子数,然后试图将这些分子都插入盒子(但如果设的数目太多,会能插入多少就插入多少)。此命令远没有Packmol那么灵活,它能做的事用Packmol都可以实现,而且Packmol会通过优化算法试图让分子尽可能好地堆叠,因此对于同样的盒子范围,靠Packmol最多能插入的分子数比用gmx insert-molecules明显要多,故能得到更紧凑的结构。
Packmol原理上也可以产生磷脂双分子膜,乃至膜蛋白这样的体系的,但是这需要设置很多约束条件才能产生靠谱的结果,而且耗时很高,结果也不理想。比如磷脂之间或者磷脂与蛋白之间会有好多空隙,这会导致模拟刚开始就有水大量跑进疏水区,后续模拟将毫无意义。
如果你用Packmol建模的目的是之后用GROMACS做模拟,别忘了GROMACS的拓扑文件里[molecules]字段里的分子按顺序完全展开后,原子的顺序必须和结构文件完全相同。为保证这一点,用于Packmol建模的分子结构文件里的原子顺序首先就得和这个分子的拓扑信息(具体来说是[atoms]字段)相对应。
pymol 旋转分子
1. 添加坐标轴
在pymol的命令行出运行 “run C:\Users\Administrator\path\to\axes.py”
axes.py with the following lines
from pymol.cgo import *
from pymol import cmd
from pymol.vfont import plain
w = 0.3 # cylinder width
l = 40.75 # cylinder length
h = 2.25 # cone hight
d = w * 2.618 # cone base diameter
obj = [CYLINDER, 0.0, 0.0, 0.0, l, 0.0, 0.0, w, 1.0, 0.0, 0.0, 1.0, 0.0, 0.0,
CYLINDER, 0.0, 0.0, 0.0, 0.0, l, 0.0, w, 0.0, 1.0, 0.0, 0.0, 1.0, 0.0,
CYLINDER, 0.0, 0.0, 0.0, 0.0, 0.0, l, w, 0.0, 0.0, 1.0, 0.0, 0.0, 1.0,
CONE, l, 0.0, 0.0, h+l, 0.0, 0.0, d, 0.0, 1.0, 0.0, 0.0, 1.0, 0.0, 0.0, 1.0, 1.0,
CONE, 0.0, l, 0.0, 0.0, h+l, 0.0, d, 0.0, 0.0, 1.0, 0.0, 0.0, 1.0, 0.0, 1.0, 1.0,
CONE, 0.0, 0.0, l, 0.0, 0.0, h+l, d, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 1.0, 1.0, 1.0]
# add labels to axes object (requires pymol version 0.8 or greater, I
# believe
cyl_text(obj,plain,[-5.,-5.,-1],'Origin',0.20,axes=[[3,0,0],[0,3,0],[0,0,3]])
cyl_text(obj,plain,[10.,0.,0.],'X',0.20,axes=[[3,0,0],[0,3,0],[0,0,3]])
cyl_text(obj,plain,[0.,10.,0.],'Y',0.20,axes=[[3,0,0],[0,3,0],[0,0,3]])
cyl_text(obj,plain,[0.,0.,10.],'Z',0.20,axes=[[3,0,0],[0,3,0],[0,0,3]])
cmd.load_cgo(obj, 'axes-1')
选择想挪动的分子,例如全部分子,如下图