计算化学公社

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

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

  [复制链接 Copy URL]

1169

帖子

7

威望

6828

eV
积分
8137

Level 6 (一方通行)

本帖最后由 fhh2626 于 2025-9-27 22:59 编辑

NAMD自由能计算教程—1、用eABF和WTM-eABF进行多维自由能计算
fhh2626  Updated in 2025/9/27
本文首发计算化学公社分子模拟板块,禁止转载。
其他教程:
使用MULE在势能面/自由能面上寻找最低反应路径  http://bbs.keinsci.com/thread-17796-1-1.html
用WTM-λABF进行炼金术自由能计算 http://bbs.keinsci.com/thread-56033-1-1.html

ChangeLog:
2025/9/27
方法名改为"WTM-eABF",略微修改了一些描述
略微修改了参数的设置,添加了搜索最优自由能路径方法的链接
添加了Gromacs的说明

2020/5/28
well-tempered-meta-eABF的用法(建议使用这个)
NAMD 2.14b1及后续版本中边界力常数的设置方法
expandBoundaries选项

目前几大MD软件当中,Gromacs的中文资料比较丰富,而其他软件的则较少。尤其是NAMD,在网上很难找到中文教程。因此,作为NAMD的developer之一,我决定开一个坑,向大家介绍一下NAMD的中高级功能,可能包括但不限于TclForce,QM/MM,自由能计算,可极化力场,广义系综等等。希望我的教程能让大家觉得NAMD比别的MD软件要强大:)
(2025.9.27更新,注意,现在Gromacs也内置了Colvars模块,也就是本教程也基本适用于Gromacs)
这一系列教程假定读者已经具有一定的MD基础,能自己建立输入文件并且跑简单的平衡模拟(对于志在科研的你来说应该做得到吧。。。),由于我比较擅长自由能计算,所以我准备先介绍自由能计算这个主题。

下面进入主题。我把自由能计算方法大致分为两类:基于“炼金术”(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的基础上添加well-tempered策略,即是well-tempered meta-eABF(WTM-eABF)方法,我个人觉得在大多数情况下,都可以无脑使用WTM-eABF方法进行基于反应坐标的自由能计算。
这里不详细讨论两种方法的原理,要详细了解原理或者引用参考文献的可以参考这几篇文献:
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(2025.9.27更新,在大体系,长时间的模拟中,我们目前更喜欢用3)
hillweight  0.5
metadynamics的峰高(kcal/mol),在meta-eABF中均可使用0.1(2025.9.27更新,在大体系,长时间的模拟中,我们目前更喜欢用0.05)
welltempered   on
biastemperature  4000
最早的meta-eABF并没有引入well-tempered算法,因为当时我觉得没啥用。后来在复杂体系计算中发现加上well-tempered算法会是的计算收敛性更好,因为eABF本身是渐近收敛的,well-tempered-MtD也是,但是单纯的metadynamics并不是渐近收敛的。
Biastemperature定义了MtD渐近收敛的速率,4000对大多数情况适用

}
提交作业,运行10 ns,对metaeabf.abf1.czar.pmf作图。


在得到自由能图以后,可以在其中寻找连接两个稳定态(极小值)的最低自由能路径,也就是真实的反应路径。详细教程见(http://bbs.keinsci.com/forum.php ... %C2%B7%BE%B6&page=1

001_eABF_metaeABF.7z

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

评分 Rate

参与人数
Participants 26
威望 +1 eV +94 收起 理由
Reason
student0618 + 3
leichuang + 5 谢谢分享
MUMBER + 5 谢谢
gwh1206 + 2 谢谢
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 精品内容

查看全部评分 View all ratings

1169

帖子

7

威望

6828

eV
积分
8137

Level 6 (一方通行)

323#
 楼主 Author| 发表于 Post on 2025-12-1 09:57:38 | 只看该作者 Only view this author
yh123 发表于 2025-11-27 19:25
这是我的colvars.in文件,格点数是14*36=504
colvarsTrajFrequency    10000
colvarsRestartFrequency  ...

theta有可能超出这个范围?你看看colvars.traj文件里面采样落到哪个范围了

2

帖子

0

威望

25

eV
积分
27

Level 2 能力者

322#
发表于 Post on 2025-11-27 19:25:08 | 只看该作者 Only view this author
fhh2626 发表于 2025-11-27 13:21
你采样空间的格点数是多少?

这是我的colvars.in文件,格点数是14*36=504
colvarsTrajFrequency    10000
colvarsRestartFrequency 1000000

colvar {
  name r
  width 0.05
  lowerBoundary 1.5
  upperBoundary 2.2

  extendedlagrangian   on
  extendedFluctuation  0.05
  extendedTimeConstant   200
  extendedTemp 298.15

  distance {
    group1 { atomNumbers 12 }
    group2 { atomNumbers 17 }
  }
}

colvar {
  name theta
  width 5
  lowerBoundary 0
  upperBoundary 180

  extendedlagrangian   on
  extendedFluctuation  5.0
  extendedTimeConstant   200
  extendedTemp 298.15

  dihedral {
    group1 {atomNumbers 12 }
    group2 {atomNumbers 16 }
    group3 {atomNumbers 27 }
    group4 {atomNumbers 30 }
  }
}
harmonicWalls {
name abf_pmf
colvars r
lowerWalls 1.5
upperWalls 2.2
lowerWallConstant 1.0
upperWallConstant 1.0
}
abf {
  colvars r theta
  fullSamples 500
  outputFreq 10000
  historyFreq 10000
  writeCZARwindowFile
}

1169

帖子

7

威望

6828

eV
积分
8137

Level 6 (一方通行)

321#
 楼主 Author| 发表于 Post on 2025-11-27 13:21:11 | 只看该作者 Only view this author
yh123 发表于 2025-11-26 21:01
付老师你好,我在使用eABF的方法计算反应坐标为键长、二面角的自由能面,这是我的zcount里的采样结果,有很 ...

你采样空间的格点数是多少?

2

帖子

0

威望

25

eV
积分
27

Level 2 能力者

320#
发表于 Post on 2025-11-26 21:01:45 | 只看该作者 Only view this author
付老师你好,我在使用eABF的方法计算反应坐标为键长、二面角的自由能面,这是我的zcount里的采样结果,有很多地方采样都为0,这个该怎么解决呢?我试了把反应时间从10ns提升到50ns,但收效甚微。

202511262059518165..png (27.29 KB, 下载次数 Times of downloads: 1)

202511262059518165..png

202511262100073390..png (23.24 KB, 下载次数 Times of downloads: 1)

202511262100073390..png

144

帖子

0

威望

1830

eV
积分
1974

Level 5 (御坂)

319#
发表于 Post on 2025-9-28 22:10:35 | 只看该作者 Only view this author
fhh2626 发表于 2025-7-22 13:39
这个好像问题不大,只是没有收敛,再跑跑应该能平

付老师,我延长模拟好像并没有解决这个问题,重复模拟有时候跑出来两边差距更大,这个有什么解决办法吗

1169

帖子

7

威望

6828

eV
积分
8137

Level 6 (一方通行)

318#
 楼主 Author| 发表于 Post on 2025-9-23 17:49:22 | 只看该作者 Only view this author
kaiyuan 发表于 2025-9-22 17:29
付老师你好,我在尝试用wtm-eabf采样一个聚合物穿过COF膜的过程(colvars+gromacs 2025.1),cv设置为某一 ...

用abf1.czar.pmf。但是你这个过程比较复杂,要充分采样才能收敛

18

帖子

0

威望

282

eV
积分
300

Level 3 能力者

317#
发表于 Post on 2025-9-22 17:29:39 | 只看该作者 Only view this author
付老师你好,我在尝试用wtm-eabf采样一个聚合物穿过COF膜的过程(colvars+gromacs 2025.1),cv设置为某一膜孔中心周围原子的质心和聚合物质心之间距离的z轴分量,用圆柱势在xy方向上限制聚合物在膜孔附近,upperwalls和lowerwalls限制聚合物不要穿越z方向上的边界
width    0.05
lowerboundary    -6
upperboundary    6
reflectingLowerBoundary    on

reflectingUpperBoundary    on
upperwalls跟lowerwalls分别为6,-6
先运行了200ns做测试,发现聚合物从膜上方靠近膜,穿过膜后到达下方溶液中,再靠近膜穿过膜到上方溶液中,然后一直待在upperwall的附近
我看了一下输出的三个pmf数据,发现abf1.pmf跟abf1.czar.pmf相似,但跟mtd的区别很大,看上去abf1/abf1.czar的值非常高,是因为我这个体系采样不充分吗


18

帖子

0

威望

282

eV
积分
300

Level 3 能力者

316#
发表于 Post on 2025-9-10 21:52:01 | 只看该作者 Only view this author
Lecly 发表于 2025-5-22 12:45
不好意思,回复较晚。暂时没有找到好的解决方法。

感觉可能是第二个CV的问题,我是参考了MARTINI力场做小分子参数化的方法,这个两个分子之间的torsion可能太不稳定了,跑其他类型的体系暂时还没碰到这个问题。

9

帖子

0

威望

446

eV
积分
455

Level 3 能力者

315#
发表于 Post on 2025-8-4 12:31:59 | 只看该作者 Only view this author
本帖最后由 Peng14 于 2025-8-4 17:14 编辑

已解决,感谢各位老师的帮助

1169

帖子

7

威望

6828

eV
积分
8137

Level 6 (一方通行)

314#
 楼主 Author| 发表于 Post on 2025-7-22 13:39:56 | 只看该作者 Only view this author
thor 发表于 2025-7-18 15:20
是这样的,左右两边收敛的地方不一样高

这个好像问题不大,只是没有收敛,再跑跑应该能平

144

帖子

0

威望

1830

eV
积分
1974

Level 5 (御坂)

313#
发表于 Post on 2025-7-18 15:20:41 | 只看该作者 Only view this author

是这样的,左右两边收敛的地方不一样高

1169

帖子

7

威望

6828

eV
积分
8137

Level 6 (一方通行)

312#
 楼主 Author| 发表于 Post on 2025-7-18 15:17:06 | 只看该作者 Only view this author
吹梦到西洲 发表于 2025-5-21 10:20
Fu老师您好,您的教程写的很棒,但是我目前碰到了一个问题想请教一下。
目前的工作是用gmx自带的colvar研 ...

你的colvars文件里面好像没有设置lowerboundary和upperboundary,也没有设置reflectingboundary。另外要画.czar.pmf文件,不是.pmf文件。

看你的轨迹,你这个虚原子被拉走了,真实CV并没有变(建议width为0.1、extendedFluctuation也为0.1)

1169

帖子

7

威望

6828

eV
积分
8137

Level 6 (一方通行)

311#
 楼主 Author| 发表于 Post on 2025-7-18 15:14:06 | 只看该作者 Only view this author
thor 发表于 2025-7-18 11:52
付博好,想请教您一个问题。我用meta-eabf(colvar+gromacs)做单个气体粒子通过二维材料覆盖一层离子液体 ...

没看到图

144

帖子

0

威望

1830

eV
积分
1974

Level 5 (御坂)

310#
发表于 Post on 2025-7-18 11:52:52 | 只看该作者 Only view this author
fhh2626 发表于 2025-1-28 00:16
是的
用distanceXY作为CV施加约束。但是修改axis选择YZ平面

付博好,想请教您一个问题。我用meta-eabf(colvar+gromacs)做单个气体粒子通过二维材料覆盖一层离子液体的膜(左侧为离子液体),发现出来的pmf收敛的两端不一样高,这是因为我的膜不对称造成的(如下图)?如果是的话,怎么解决这个问题?

7

帖子

0

威望

293

eV
积分
300

Level 3 能力者

309#
发表于 Post on 2025-5-22 12:45:54 | 只看该作者 Only view this author
kaiyuan 发表于 2025-5-17 17:20
请问你有找到解决办法吗,我关掉MtD只用eABF会出现段错误

不好意思,回复较晚。暂时没有找到好的解决方法。

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

GMT+8, 2026-1-24 10:32 , Processed in 0.183147 second(s), 25 queries , Gzip On.

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