本帖最后由 fhh2626 于 2020-5-28 16:18 编辑
NAMD自由能计算教程—1、用eABF和meta-eABF进行多维自由能计算 fhh2626 Updated in 2020/5/28 本文首发计算化学公社分子模拟板块,禁止转载。 ChangeLog: 2020/5/28 well-tempered-meta-eABF的用法(建议使用这个) NAMD 2.14b1及后续版本中边界力常数的设置方法 expandBoundaries选项
目前几大MD软件当中,Gromacs的中文资料比较丰富,而其他软件的则较少。尤其是NAMD,在网上很难找到中文教程。因此,作为NAMD的developer之一,我决定开一个坑,向大家介绍一下NAMD的中高级功能,可能包括但不限于TclForce,QM/MM,自由能计算,可极化力场,广义系综等等。希望我的教程能让大家觉得NAMD比别的MD软件要强大:) 这一系列教程假定读者已经具有一定的MD基础,能自己建立输入文件并且跑简单的平衡模拟(对于志在科研的你来说应该做得到吧。。。),由于我比较擅长自由能计算,所以我准备先介绍自由能计算这个主题,暂定有以下几个内容: 1、用eABF和meta-eABF进行多维自由能计算 2、谈谈metadynamics和其estimator 3、umbrella sampling、选择反应坐标和自由能计算当中的tricks 4、Colvars模块的高级应用 5、“炼金术”方法——FEP和TI 6、精确结合自由能计算和BFEE插件 下面进入主题。我把自由能计算方法大致分为两类:基于“炼金术”(alchemistry)的自由能计算方法和基于反应坐标的自由能计算方法。前者主要包括FEP和TI,由于其可以直接消失生长原子而不顾及化学原理,形似“炼金术”而得名。后者则包括ABF、Metadynamics、Umbrella Sampling和SMD四大方法。SMD是非平衡模拟方法,这里不多介绍,这个主题主要介绍其他的自由能计算方法。 如果要使用NAMD的话,则推荐使用ABF类方法,原因是NAMD的作者之一Chris Chipot和NAMD中自由能计算模块Colvars的作者Jérôme Hénin都是ABF阵营的兄弟。所以NAMD对ABF的支持较好。 在这个主题中,我会向大家介绍两种ABF的衍生方法,eABF和meta-eABF方法。其中eABF方法的基本原理是由Tony Lelièvre和Yang Wei首先提出的,由我最先实现在NAMD当中,之后Jerome提出了eABF的两个改进(czar estimator和pABF)。eABF目前是ABF的公认上位方法,即在任何情况下都应该使用eABF和不是传统ABF。meta-eABF是我最近提出的,(私货注意)是目前最快的自由能计算方法之一,(私货注意)我个人认为在任何情况下都可以使用meta-eABF代替eABF和metadynamics。 这里不详细讨论两种方法的原理,我会在第二讲和第三讲中涉及到各个方法的原理。此外,要详细了解原理或者引用参考文献的可以参考这几篇文献: eABF: Fu et al. J. Chem. Theory Comput., 2016, 12(8), pp 3506–3513 Lesage et al. J. Phys. Chem. B, 2017, 121(15), pp 3676–3685 meta-eABF: Fu et al. J. Phys. Chem. Lett., 2018, 9(16), pp 4738–4745 Well-tempered-meta-eABF Fu et al. Acc. Chem. Res. 2019, 52 (11), 3254–3264 下面以NANMA为例来介绍这两种方法的使用。NANMA又称alanine dipeptide,其构象变化通常可以用两个二面角来描述,这里我们将计算沿着这两个二面角的自由能变化:
这里假定读者具有能自己生成大部分输入文件的基础,生成好的pdb,psf和namd文件已经提供在附件当中。 namd文件中有以下几句话: colvars on 这句话表示激活NAMD当中的Colvars模块 colvarsConfig colvarsA.in 这句话表示Colvars模块的输入文件名为colvarsA.in 然后打开colvarsA.in文件 colvarsTrajFrequency 20000 colvarsRestartFrequency 1000000 上面两句分别定义了colvars文件记录数据的频率和输出的频率 colvar { name phi 这里定义了一个叫phi的几何变量 width 5.0 该几何变量记录数据的精度 lowerboundary -180 upperboundary 180 记录的上下界 extendedlagrangian on extendedFluctuation 5.0 extendedTimeConstant 200 定义体系中的扩展自由度,即eABF中的“e”,这里不用深究含义,对于任何体系均可设置extendedFluctuation=width和extendedTimeConstant=200 dihedral { group1 { atomnumbers 13 } group2 { atomnumbers 15 } group3 { atomnumbers 17 } group4 { atomnumbers 1 } } 对于dihedral,即二面角colvar,需要定义4个原子(团),这里定义了13-15-17-1的二面角。定义原子(团)可以直接指定pdb当中的原子序号,也可以通过给一个reference pdb,在里边将某一列由0改为1实现,后面那种方法通常用于定义大原子团 }
colvar { name psi
width 5.0 lowerboundary -180 upperboundary 180
extendedlagrangian on extendedFluctuation 5.0 extendedTimeConstant 200
dihedral { group1 { atomnumbers 3 } group2 { atomnumbers 1 } group3 { atomnumbers 17 } group4 { atomnumbers 15 } } }
NAMD 2.14及之后的版本修改了边界力常数的定义方式,现在需要额外定义一个harmonicWalls block #harmonicWalls { # name mywalls # colvars phi psi # lowerWalls -180 -180 # upperWalls 180 180 # forceConstant 1.0 #} 分别定义phi psi的边界墙,力常数为1.0,注意单位是对应变量的width值,也就是说力常数会随着反应坐标种类/width不同而自动scale。1.0的力常数适用于大多数情况,在这里由于二面角是周期性的,所以无需定义边界力常数
abf { 定义跑的作业为ABF,因为之前设了扩展自由度,因此跑的是eABF colvars phi psi 表示针对上面定义的phi和psi colvars进行自由能能计算 fullSamples 500 在每一个width区间内的平衡步数,500-2000适用于绝大多数情况 historyfreq 1000000 写history文件的频率 writeCZARwindowFile 使用czar estimator时产生一些中间文件 } 提交作业,跑50 ns以后(至少有一台工作站吧!),对eabf.czar.pmf作图。
接下来尝试一下meta-eABF,meta-eABF是eABF和metadynamics的结合,打开colvarsA.in可以看到,配置文件绝大部分都跟eABF一样,多出了以下几句话, subtractAppliedForce on expandboundaries on 这是处理metadynamics高斯峰的一个选项,不写这个也不会影响计算结果,但是计算效率比较低
这是计算eABF偏置力的一个内部选项,使用meta-eABF时必须打开 metadynamics { colvars phi psi 这里同时向定义的两个几何变量的扩展自由度上施加metadynamics偏置力,即meta-eABF hillwidth 5 metadynamics的峰宽(单位为width),在meta-eABF中均可使用5 hillweight 0.1 metadynamics的峰高(kcal/mol),在meta-eABF中均可使用0.1 welltempered on biastemperature 4000 最早的meta-eABF并没有引入well-tempered算法,因为当时我觉得没啥用。后来在复杂体系计算中发现加上well-tempered算法会是的计算收敛性更好,因为eABF本身是渐近收敛的,well-tempered-MtD也是,但是单纯的metadynamics并不是渐近收敛的。 Biastemperature定义了MtD渐近收敛的速率,4000对大多数情况适用
} 提交作业,运行10 ns,对metaeabf.abf1.czar.pmf作图。
|