解析gromacs的restraint、constraint和freeze | 宜武汇-ag真人国际厅网站

声明:此文最初发布于sobereva的百度博客,后在网上被一堆博客胡乱转载,却不注明真正来源。此文真正原文在此!鄙视胡乱转载但是不注明出处的人!

gromacs里面人为固定坐标有三大类方式:restraint、constraint和freeze,有各自的特点。

包括位置限制(也就是常说的正式md前首先要进行的限制性动力学),角度限制,二面角限制,方向限制,距离限制。限制的特点是可以让被限制的东西在一定范围内运动,而非彻底固定住,实际上就是施加一个谐振势来限制其移动。

位置限制就不必说了,posre.itp里面[ position_restraints ]定义的就是,默认相对于最初位置进行限制。

角度限制包括两类,一类是两对儿原子间角度的限制,用[ angle_restraints ]来指定,例如:

[ angle_restraints ]

; i   j    k    l    type   theta0     fc     multiplicity

651 1211 1683 1211    1      67.0     1500         1

说明限制651-1211与1683-1211原子对儿之间的夹角在67度附近,力常数为1500。type无用。

角度限制另一类包括一对儿原子与z轴夹角角度的限制,用[ angle_restraints_z ]来指定,例如:

[ angle_restraints_z ]

; #1 #2  type  theta0    fc     multiplicity

22  45     1     90      500         2

说明限制22与45号原子的连线与z轴垂直,力常数是500,多重性是2,使得90度夹角时候限制势能最低,0和180度时最高。type无用。

二面角限制使用[ dihedral_restraints ]段落来定义,实际上improper项就是用二面角限制方式限制的。例如

[ dihedral_restraints ]

;   i    j     k    l    type  label  phi  dphi  kfac  power

    5    7     9    15     1      1  180     0     1      2

被限制的是5,7,9,15原子组成的二面角,type总是1,label没用,phi是参考角,dphi是超过参考角多少度开始使用限制势,power没用。kfac乘上mdp中的dihre_fc将作为限制势力常数。

最后在mdp中加入例如:

dihre               =  simple

dihre_fc            =  100

dihre_tau           =  0.0

nstdihreout         =  50

距离限制在[ distance_restraints ]里面定义,比如

[ distance_restraints ]

; i j type index type low up1 up2 fac

10 16 1 0 1 0.0 0.3 0.4 1.0

type总是1,index是计算的顺序,如果几个项index都一样,比如10-28和10-46,就同时计算。势能图见gmx3.3手册p60,low,up1,up2分别指图上的r0,r1,r2,可见原子间距离在low至up1区间内是不受限制的,这种方式可以达成nmr限制。fac是指这个因子乘上mdp中disre_fc作为限制势力常数。

也可以定义两个原子在[ bonds ]里用bond type 6,就是个和普通键一样的谐振子势,但是这两个原子间被认为没有键连。

应当注意以上限制方式中[]段落应当紧接着写在被限制分子的.itp后面(或者说对应的[ moleculetype ]后面),这样程序才知道其中的原子编号指的是哪种分子中的原子。

用shake或lincs方法固定住原子之间的相对位置,固定的几何关系是绝对不变的。分为两种类型,type 1和type 2的效果其实一样,但是type 1被认为将这两个原子键连了,而type 2没有键连。一般被键连的原子间都不计算相互作用(在[ moleculetype ]里定义nrexcl来控制相隔几个键内的原子间不计算相互作用,也就是相当于被列进[ exclusions ]),比如type 1时一般就不计算这两个原子间的非键作用了,而type 2时仍然计算。

用哪种约束算法用.mdp里的constraint_algorithm设定,默认lincs。.mdp里的比如constraints = all-bonds也是应用这种约束方法,也就是约束住所有[ bonds ]项,原本[ bonds ]的设定就被覆盖了,即不体现成键效果(来回振动),只体现约束效果(距离固定不变)。

自定义约束项可以在拓扑文件中这么写:

[ constraints ]

1 8 2 0.153    //原子1 原子2 type 约束的距离(0.153nm)。

之后1 2原子距离会固定保持在0.153nm,即便1 2原子在[ bonds ]中定义是成键的也是如此。

注意,.mdp里如果写constraint=none,只是说批量约束hbonds、all-bonds等没有了,在拓扑文件中自定义的constraint约束项仍然生效。

顺便一提,在默认情况下,gromacs的水的结构是被settle算法约束住的,也就是所说的rigid水,settle就是专用于水的约束算法,在spce.itp里面[settles]即是定义,不计算分子内氢、氧彼此之间的相互作用。而如果.mdp里面设define=-dflexible,也就是激活了#ifdef flexible后面那段,则不用settle算法约束结构,而是照常按照谐振势的键长键角项控制水的结构。-dflexible开不开对溶质没什么影响,对计算速度也没什么影响(有一丁点减慢可忽略)。除了能量最小化外,都不要用-dflexible,因为拟合水、力场参数的时候都一般是认为水是刚性(rigid water)情况下进行的(比如spc水和使用spc水的gromos力场),即约束住键长键角,开了flexible实际上参数还得再作调整。而在gmx4中,在能量最小化时自动是限制方法(restraint),就没必要设-dflexible了,从此可忘掉它。.mdp中constraints以及constraint_algorithm设定皆与水的约束无关,除非define=-dflexible,否则水都是用settle算法约束住。

彻底冻在最初坐标,一点也不动。但仍然计算与其它原子间的相互作用,所以并不会省时间,除非是写进[ exclusions ]。

使用时在.mdp里面定义比如:

freezegrps = protein   设蛋白冻住

freezedim = y y y    设冻住的方向,分别对应x y z轴。

同时使用freezegrps, constraints, pressure coupling可能会有问题,不要组合使用。

注意绝对坐标与相对坐标的限制方式的不同,比如想保持一个结构的刚性,用位置限制或freeze是原地不动,而通过约束方法对分子内的结构进行约束,结构仍然可以有整体平动和转动。应根据实际情况选择。

原文链接:http://sobereva.com/10

网络摘文,本文作者:15h,如若转载,请注明出处:https://www.15cov.cn/2023/08/27/解析gromacs的restraint、constraint和freeze/

发表评论

邮箱地址不会被公开。 必填项已用*标注

网站地图