计算化学公社

 找回密码 Forget password
 注册 Register
Views: 178985|回复 Reply: 274
打印 Print 上一主题 Last thread 下一主题 Next thread

[NAMD] NAMD自由能计算教程—1、用eABF和meta-eABF进行多维自由能计算

  [复制链接 Copy URL]

1093

帖子

6

威望

6269

eV
积分
7482

Level 6 (一方通行)

本帖最后由 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作图。


001_eABF_metaeABF.7z

147.28 KB, 下载次数 Times of downloads: 409

评分 Rate

参与人数
Participants 22
威望 +1 eV +79 收起 理由
Reason
goodname + 5 好物!
Elden、Lord + 1 好物!
国小药 + 2 好物!
HZW + 2 好物!
guchen + 1 牛!
王维韬 + 3
mengzi + 5 好物!
bobosiji + 3 牛!
Yaqi + 5 好物!谢谢!
kaiwanxiao + 4 谢谢
cadz + 5 赞!
千张 + 2 精品内容
xcandy + 4 谢谢
yjmaxpayne + 5 好物!
Novice + 3 好物!
晴天 + 5 精品内容
kunkun + 5 精品内容
meatball1982 + 5 这是真的好东西。
Shine剪水 + 5 谢谢分享
ene + 5 牛!

查看全部评分 View all ratings

2

帖子

0

威望

35

eV
积分
37

Level 2 能力者

274#
发表于 Post on 4 day ago | 只看该作者 Only view this author
Fu老师,请问eABF方法能用来研究纳米晶体不同取向聚集过程的PMF吗?我的纳米晶体模型大约2000个原子,按照这个帖子尝试算了一下真空中两个纳米晶体的聚集,初始间距为15A,采样范围2~17A,模拟开始后间距会迅速迅速减小到2A,然后一直在那个位置小幅振荡,是我力常数设置太小的原因吗?还是采样时间不够?

19

帖子

0

威望

123

eV
积分
142

Level 2 能力者

273#
发表于 Post on 6 day ago | 只看该作者 Only view this author
付老师您好,请问一下 plumed 里面,funnel 能用于 wtm-eabf 吗

1093

帖子

6

威望

6269

eV
积分
7482

Level 6 (一方通行)

272#
 楼主 Author| 发表于 Post on 2024-11-4 16:57:16 | 只看该作者 Only view this author
DwyaneWan 发表于 2024-11-4 11:47
Fu老师,请问在用colvars进行coordination numbers作为模拟CV,但我想把其中两个原子作为CV和自由能的关 ...

选择CV以后用eABF或者其他方法做自由能计算

63

帖子

0

威望

1087

eV
积分
1150

Level 4 (黑子)

271#
发表于 Post on 2024-11-4 11:47:45 | 只看该作者 Only view this author
fhh2626 发表于 2024-10-16 13:02
不能,把energy作为CV没有任何作用

Fu老师,请问在用colvars进行coordination numbers作为模拟CV,但我想把其中两个原子作为CV和自由能的关系表达出来,应该用何种方式进行reweight呢?

1093

帖子

6

威望

6269

eV
积分
7482

Level 6 (一方通行)

270#
 楼主 Author| 发表于 Post on 2024-10-16 13:02:46 | 只看该作者 Only view this author
DwyaneWan 发表于 2024-10-15 23:36
fu老师,我看plumed可以支持把体系势能项作为CV,然后结合特定的结构序数参数来研究结晶过程。不知道colv ...

不能,把energy作为CV没有任何作用

63

帖子

0

威望

1087

eV
积分
1150

Level 4 (黑子)

269#
发表于 Post on 2024-10-15 23:36:11 | 只看该作者 Only view this author
fhh2626 发表于 2024-10-15 17:31
技术上可以实现,但是所有原子coordinate维数太多了啊,没有办法有效增强采样的。

水结冰过程已经研究 ...

fu老师,我看plumed可以支持把体系势能项作为CV,然后结合特定的结构序数参数来研究结晶过程。不知道colvars能否把energy项目作为CV?

1093

帖子

6

威望

6269

eV
积分
7482

Level 6 (一方通行)

268#
 楼主 Author| 发表于 Post on 2024-10-15 17:31:09 | 只看该作者 Only view this author
DwyaneWan 发表于 2024-10-15 00:11
fu老师你好,不知道能否把所有原子coordinate设置成CV,例如用WTM-eabf方法去研究在三维周期性结构的水分子 ...

技术上可以实现,但是所有原子coordinate维数太多了啊,没有办法有效增强采样的。

水结冰过程已经研究得很充分了,可以参考一下近期高毅勤课题组的文章,是用增强采样方法研究水结冰的https://pubs.acs.org/doi/10.1021/acs.jpclett.2c02176

1093

帖子

6

威望

6269

eV
积分
7482

Level 6 (一方通行)

267#
 楼主 Author| 发表于 Post on 2024-10-15 17:27:09 | 只看该作者 Only view this author
yyf123 发表于 2024-10-13 01:41
付博士您好,我发现我的配体和蛋白有阳离子π相互作用,请教一下您,这个适合用什么力场?
charmm-gui 里 ...

可以勾上的
不一定,虽然Drude原理上可能比CHARMM力场更准确,但是Drude本身优化得并没有CHARMM力场好,尤其是小分子力场

63

帖子

0

威望

1087

eV
积分
1150

Level 4 (黑子)

266#
发表于 Post on 2024-10-15 00:11:10 | 只看该作者 Only view this author
fu老师你好,不知道能否把所有原子coordinate设置成CV,例如用WTM-eabf方法去研究在三维周期性结构的水分子结冰的自由能过程?

19

帖子

0

威望

123

eV
积分
142

Level 2 能力者

265#
发表于 Post on 2024-10-13 01:41:45 | 只看该作者 Only view this author
本帖最后由 yyf123 于 2024-10-14 12:14 编辑

付博士您好,我发现我的配体和蛋白有阳离子π相互作用,请教一下您,这个适合用什么力场?
charmm-gui 里面 charmm36m 力场下面的 WYF parameter for cation-pi interactions 这个设置是否能准确描述呢?
对于药物小分子和蛋白的结合,是否drude力场一定比普通力场要好,可以无脑使用呢?
谢谢您的指导

123.png (11.99 KB, 下载次数 Times of downloads: 0)

123.png

1234.png (75.35 KB, 下载次数 Times of downloads: 0)

1234.png

1093

帖子

6

威望

6269

eV
积分
7482

Level 6 (一方通行)

264#
 楼主 Author| 发表于 Post on 2024-10-7 22:59:22 | 只看该作者 Only view this author
docshen777 发表于 2024-10-7 19:46
请问付老师,我的计算两个CV,第一个CV的范围是2.3-6.8 A,但是长期的采样只能在5.0-6.8A进行,下不去5.0, ...

有用,但是建议你看看你的CV是否合适。这种情况多半是CV不合适引起的

36

帖子

0

威望

294

eV
积分
330

Level 3 能力者

263#
发表于 Post on 2024-10-7 19:46:19 | 只看该作者 Only view this author
请问付老师,我的计算两个CV,第一个CV的范围是2.3-6.8 A,但是长期的采样只能在5.0-6.8A进行,下不去5.0,这种情况应该改变哪个参数?改变wtm的hill高和宽有用么?还是说改变eabf的参数?

21

帖子

0

威望

864

eV
积分
885

Level 4 (黑子)

262#
发表于 Post on 2024-9-28 23:41:29 | 只看该作者 Only view this author
fhh2626 发表于 2024-9-28 14:01
可以,但是要仔细设计CV,五维的自由能计算肯定跑不动的(就算你真的跑得动也没办法可视化分析)

好的,谢谢老师
确实5个cv太多了,我看看能否简化

1093

帖子

6

威望

6269

eV
积分
7482

Level 6 (一方通行)

261#
 楼主 Author| 发表于 Post on 2024-9-28 14:01:13 | 只看该作者 Only view this author
po390 发表于 2024-9-27 18:08
付老师,您好
这个meta-eABF直接应用于NAMD的QM/MM中是否可行,例如我有涉及化学反应的多个原子直接的五个 ...

可以,但是要仔细设计CV,五维的自由能计算肯定跑不动的(就算你真的跑得动也没办法可视化分析)

21

帖子

0

威望

864

eV
积分
885

Level 4 (黑子)

260#
发表于 Post on 2024-9-27 18:08:40 | 只看该作者 Only view this author
付老师,您好
这个meta-eABF直接应用于NAMD的QM/MM中是否可行,例如我有涉及化学反应的多个原子直接的五个距离作为反应变量,获得自由能差异,可否认为就是反应能垒?

手机版 Mobile version|北京科音自然科学研究中心 Beijing Kein Research Center for Natural Sciences|京公网安备 11010502035419号|计算化学公社 — 北京科音旗下高水平计算化学交流论坛 ( 京ICP备14038949号-1 )|网站地图

GMT+8, 2024-11-23 09:49 , Processed in 0.213061 second(s), 26 queries , Gzip On.

快速回复 返回顶部 返回列表 Return to list