本帖最后由 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作图。
| 
[复制链接 Copy URL]
发表于 Post on 2018-9-4 10:59:51
收藏 Add to favorites
楼主 Author
