计算化学公社

 找回密码 Forget password
 注册 Register
Views: 178962|回复 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

5万

帖子

99

威望

5万

eV
积分
112349

管理员

公社社长

2#
发表于 Post on 2018-9-4 11:26:47 | 只看该作者 Only view this author
“分子模拟论坛”是哪个论坛?
北京科音自然科学研究中心http://www.keinsci.com)致力于计算化学的发展和传播,长期开办高质量的各种计算化学类培训:初级量子化学培训班中级量子化学培训班高级量子化学培训班量子化学波函数分析与Multiwfn程序培训班分子动力学与GROMACS培训班CP2K第一性原理计算培训班,内容介绍以及往届资料购买请点击相应链接查看。这些培训是计算化学从零快速入门以及进一步全面系统性提升研究水平的高速路!培训各种常见问题见《北京科音办的培训班FAQ》
欢迎加入“北京科音”微信公众号获取北京科音培训的最新消息、避免错过网上有价值的计算化学文章!
欢迎加入人气非常高、专业性特别强的综合性理论与计算化学交流QQ群“思想家公社QQ群”:1号:18616395,2号:466017436,3号:764390338,搜索群号能搜到哪个说明目前哪个能加,合计9000人。北京科音培训班的学员在群中可申请VIP头衔,提问将得到群主Sobereva的最优先解答。
思想家公社的门口Blog:http://sobereva.com(发布大量原创计算化学相关博文)
Multiwfn主页:http://sobereva.com/multiwfn(十分强大的量子化学波函数分析程序)
ResearchGate:https://www.researchgate.net/profile/Tian_Lu
Money and papers are rubbish, get a real life!

1093

帖子

6

威望

6269

eV
积分
7482

Level 6 (一方通行)

3#
 楼主 Author| 发表于 Post on 2018-9-4 11:58:02 | 只看该作者 Only view this author
sobereva 发表于 2018-9-4 11:26
“分子模拟论坛”是哪个论坛?

脑残了,我想说的是计算化学公社的分子模拟板块

213

帖子

1

威望

2230

eV
积分
2463

Level 5 (御坂)

4#
发表于 Post on 2018-9-4 14:42:17 | 只看该作者 Only view this author
这个比较好!
学习学习。

294

帖子

0

威望

3528

eV
积分
3822

Level 5 (御坂)

5#
发表于 Post on 2018-9-4 19:25:10 | 只看该作者 Only view this author
顶顶,NAMD很不错

161

帖子

0

威望

605

eV
积分
766

Level 4 (黑子)

蓝卫兵

6#
发表于 Post on 2018-9-5 00:49:09 | 只看该作者 Only view this author
请问fu博士:这个meta-eABF在namd的哪个版本之后可以用?谢谢
B样条插值
个人专栏https://zhuanlan.zhihu.com/p/21936803

1093

帖子

6

威望

6269

eV
积分
7482

Level 6 (一方通行)

7#
 楼主 Author| 发表于 Post on 2018-9-5 16:30:52 | 只看该作者 Only view this author
pyscf 发表于 2018-9-5 00:49
请问fu博士:这个meta-eABF在namd的哪个版本之后可以用?谢谢

2.12好像不行,只有最新的Nightly version可以,需要自己编译,编译的方法在我的新帖里面

334

帖子

0

威望

2353

eV
积分
2687

Level 5 (御坂)

8#
发表于 Post on 2018-11-13 14:06:14 | 只看该作者 Only view this author
谢谢楼主,感觉很有意思。
请教一下,如果是比较复杂的CV,比如 多个距离的平方和 或 RMSD的加权平均,可以直接用吗?如果不行的话,可以利用Plumed实现吗?

1093

帖子

6

威望

6269

eV
积分
7482

Level 6 (一方通行)

9#
 楼主 Author| 发表于 Post on 2018-11-14 22:10:18 | 只看该作者 Only view this author
霜晨月 发表于 2018-11-13 14:06
谢谢楼主,感觉很有意思。
请教一下,如果是比较复杂的CV,比如 多个距离的平方和 或 RMSD的加权平均,可 ...

这些都不算复杂的CV,前者可以用colvars里面的custom function实现,后者可以用Linear and polynomial combinations of components实现。

如果是复杂的CV,可以用Scripted functions实现

334

帖子

0

威望

2353

eV
积分
2687

Level 5 (御坂)

10#
发表于 Post on 2018-11-17 20:42:00 | 只看该作者 Only view this author
fhh2626 发表于 2018-11-14 22:10
这些都不算复杂的CV,前者可以用colvars里面的custom function实现,后者可以用Linear and polynomial co ...

谢谢fu博士,看来NAMD比我想像的强大得多。

49

帖子

0

威望

840

eV
积分
889

Level 4 (黑子)

11#
发表于 Post on 2019-3-1 20:05:02 | 只看该作者 Only view this author
本帖最后由 16aDream 于 2019-3-1 20:17 编辑

Amber里面的ABMD模块是用abp计算的,这与abf或者eabf相比哪种更推荐?(好像amber没法用abf)另外我在使用ABMD的时候做了几个测试发现体系还是被局部极小值困住了,这是因为flooding time取的不合理吗,还是因为模拟时间不够长?

1093

帖子

6

威望

6269

eV
积分
7482

Level 6 (一方通行)

12#
 楼主 Author| 发表于 Post on 2019-3-2 09:14:43 | 只看该作者 Only view this author
16aDream 发表于 2019-3-1 20:05
Amber里面的ABMD模块是用abp计算的,这与abf或者eabf相比哪种更推荐?(好像amber没法用abf)另外我在使用ABM ...

abmd是metadynamics的一个变种,陷入极小值有可能是模拟时间不够长,也有可能是反应坐标选得不好或者MtD的参数设得不好

4

帖子

0

威望

15

eV
积分
19

Level 1 能力者

13#
发表于 Post on 2019-3-14 15:08:24 | 只看该作者 Only view this author
目前还属于在学习阶段,感觉代码操作还是有点困难,不知道有没有图形界面的MD软件学习使用呢

1093

帖子

6

威望

6269

eV
积分
7482

Level 6 (一方通行)

14#
 楼主 Author| 发表于 Post on 2019-3-14 16:45:54 | 只看该作者 Only view this author
sixofgad 发表于 2019-3-14 15:08
目前还属于在学习阶段,感觉代码操作还是有点困难,不知道有没有图形界面的MD软件学习使用呢

其实现在跑模拟已经很简单了,拿NAMD为例,跑模拟最少只需要3个文件,pdb(定义初始原子坐标),psf(定义原子连接方式、电荷和原子类型),prm(力场参数)和conf(模拟参数)文件,如果你跑的是生物体系(比如蛋白质、膜)的话,前面三种文件都可以用charmm-gui自动生成,你需要的就只是conf文件,这个文件的格式就类似下面的样子:

Structure    xxx.psf
Coordinates     xxx.pdb
Temperature 300    ;定义模拟温度为300K
......

你只要做完NAMD官网上任意一个tutorial,就可以用里面的conf文件作为模板,改几个参数适应你自己的体系就行了,任何图形界面的软件都只会把事情弄得更复杂

4

帖子

0

威望

15

eV
积分
19

Level 1 能力者

15#
发表于 Post on 2019-3-14 18:46:38 | 只看该作者 Only view this author
fhh2626 发表于 2019-3-14 16:45
其实现在跑模拟已经很简单了,拿NAMD为例,跑模拟最少只需要3个文件,pdb(定义初始原子坐标),psf(定 ...

哦哦哦,这样的,不过NAMD的中文教程是真的少,有点基础教程什么的也好啊,楼主有木有兴趣做一个,感觉你是个大神哈哈哈

本版积分规则 Credits rule

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

GMT+8, 2024-11-23 05:13 , Processed in 0.217629 second(s), 31 queries , Gzip On.

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