计算化学公社
标题: NAMD自由能计算教程—2、用WTM-λABF进行炼金术自由能计算 [打印本页]
作者Author: fhh2626 时间: 2025-9-27 22:58
标题: NAMD自由能计算教程—2、用WTM-λABF进行炼金术自由能计算
本帖最后由 fhh2626 于 2025-12-23 10:56 编辑
NAMD自由能计算教程—2、用WTM-λABF进行炼金术自由能计算
fhh2626 Updated in 2025/9/27
本文首发计算化学公社分子模拟板块，禁止转载。
注：截止2025.10.9，目前还是只有在gitlab上下载NAMD的main branch源代码，并自己编译后才能跑WTM-λABF模拟。编译教程可以在论坛里面搜索
"炼金术“（alchemistry/alchemical transformation）自由能计算方法计算的是将体系一部分原子消失、生长或者转变为其他原子类型的自由能变化，在溶解自由能、绝对/相对结合自由能、pKa计算中非常常见。目前最广泛使用的方法是FEP（free-energy perturbation）和TI（thermodynamic integration）方法。最近，我们将WTM-eABF方法扩展到了这个领域，大致的思路就是将alchemical变量λ作为一个扩展自由度引入体系中进行模拟，可以通过∂U'/∂λ计算扩展体系在该方向的受力，然后引入增强采样算法加强沿λ方向的采样。
WTM-λABF相比于传统FEP等方法的优点主要有下面几个：
1. 有增强采样算法的加持，体系沿λ方向探索快；
2. λ连续变化，适用于大体系的alchemical transformation（FEP或TI方法在消失/生长大分子时因为需要大量的λ-中间态以满足微扰条件，导致收敛困难）；
3. 体系可以沿λ方向反复采样，如果λ变化时伴随体系的结构变化，探索这些变化的能力强（FEP或TI方法是在不同λ-中间态下独立模拟，探索这一类变化比较困难）
在NAMD中用WTM-λABF进行炼金术自由能计算非常简单，这里以水分子的溶解自由能计算为例：
计算水分子的溶解自由能，就是计算将一个水分子从真空中移到溶液中的自由能变化（这里稍微有些糙，严谨一些还要考虑在真空中的浓度，不过影响很小），也就是在体系中“消失”一个水分子的自由能变化的负值，或者“生长”一个水分子的自由能变化。
在NAMD中实现这类型的alchemical transformation很简单，首先要准备一个PDB文件，这个PDB文件可以直接由模拟初始结构文件复制而来，唯一的区别是把要“消失”的原子对应的B列标记为-1，要“生长”的原子的B列标记为1，其他原子为0。这里我们准备了一个简单的水盒子，标记了其中的一个水分子：
(, 下载次数 Times of downloads: 2)
相对于平衡模拟，在配置文件（alchem.conf）中需要有：
computeEnergies 1 # alchemical transformation必需
colvars on
colvarsConfig colvars.in # Colvars配置文件
alch on
alchType TI
alchFile alchem.pdb # 这就是标记的PDB文件
alchCol B
alchOutFile output/alchem.fepout
alchOutFreq 5000
alchVdwLambdaEnd 0.7 # 如果是消失过程：lambda=0.3时开始消失范德华相互作用，lambda=1时完全消失，生长时则相反
alchElecLambdaStart 0.5 # 如果是消失过程：lambda=0时开始消失静电相互作用，lambda=0.5时完全消失，生长时则相反
alchDecouple on
alchlambda 0
注意NAMD通过“alchVdwLambdaEnd”和“alchElecLambdaStart”两个参数调整范德华和静电相互作用的消失/生长策略，现在可以暂时不动它们。
然后准备Colvars配置文件（colvars.in）：
下面这一段定义了lambda反应坐标：
colvar {
name l
extendedLagrangian on
extendedMass 150000
extendedLangevinDamping 1000
lowerBoundary 0
upperBoundary 1
reflectingLowerBoundary
reflectingUpperBoundary
width 0.005
subtractAppliedForce
expandboundaries
alchLambda {
}
outputTotalForce
outputAppliedforce
}
这里唯一需要改动的是width 0.005，这个数决定了以多大的宽度划分λ-中间态，这个中间态仅仅用来给ABF和MtD算法计算偏置势，模拟过程中λ的变化是连续的，因此中间态的个数可以比FEP/TI模拟多很多。根据需要消失/生长的部分大小，设置100/200/1000都是完全没问题的，这里我设置了200个中间态。其他地方都不用改动。
然后添加ABF和MtD偏置势，这里和WTM-eABF完全一样：
abf {
colvars l
FullSamples 10000
historyfreq 50000
}
metadynamics {
colvars l
hillWidth 3.0
hillWeight 0.05
wellTempered on
biasTemperature 4000
}
跑模拟，查看输出的alchem.abf1.pmf文件，我们可以画出沿λ的自由能变化：
(, 下载次数 Times of downloads: 4)
当然我们一般只关心λ=0和λ=1时的自由能差值，这里ΔG=6.2 kcal/mol，与实验值-6.3 kcal/mol相符。
这个例子中用到的所有文件如下，其他类型的alchemical transformation道理相同：
(, 下载次数 Times of downloads: 179)
参考文献：J. Phys. Chem. Lett. 2025, 16, 4419–4427，Acc. Chem. Res. 2025, DOI:10.1021/acs.accounts.5c00666
作者Author: CrysW555 时间: 2025-10-3 12:34
您好，请问这个算法在gromacs中能否应用？
作者Author: fhh2626 时间: 2025-10-9 10:48
应该是不行
作者Author: enthalpy 时间: 2025-10-9 22:29
gromacs 的plumed 也支持类似的算法： Alchemical Metadynamics: Adding Alchemical Variables to Metadynamics to Enhance Sampling in Free Energy Calculations
Wei-Tse Hsu; Valerio Piomponi; Pascal T. Merz; Giovanni Bussi; Michael R. Shirts*
28, 2023 JCTC
作者Author: CrysW555 时间: 2025-11-11 15:14
感谢感谢！
作者Author: 谢林 时间: 2025-11-20 12:27
大佬LAMMPS能使用么
作者Author: fhh2626 时间: 2025-11-26 15:06
不能
作者Author: wth1219 时间: 2025-11-26 23:45
这么好的方法，强烈建议考虑与LAMMPS和GROMACS结合，方便绝大部分科研工作者使用的同时，也能非常方便地结合机器学习力场使用。
作者Author: beyond 时间: 2025-12-2 19:16
alchType TI
这里用到的是TI的方法，自由能曲线与传统FEP曲线一致，不用积分，只看lambda 0 1的自由能差。
是不是在这里不管用TI还是FEP，最后得到的自由能曲线都是如此？
谢谢付博！
作者Author: fhh2626 时间: 2025-12-9 11:08
这里是用TI的方法算偏置力，本身和TI、FEP没关系
作者Author: beyond 时间: 2025-12-9 21:45
谢谢付博！明白了，是我搞混了。。。
作者Author: Lecly 时间: 2025-12-22 14:26
付博士好。请问，若计算小分子配体的绝对结合自由能，WTM-λABF和LDDM两种方法何者更优？看您对WTM-λABF优点的描述，是不是WTM-λABF要更好呢？谢谢付博士。
作者Author: fhh2626 时间: 2025-12-23 09:59
我觉得是LDDM更好。如果是相同条件下的模拟，那WTM-λABF更好（LDDM是基于FEP的），但是目前技术问题，WTM-λABF还实现不了LDDM那种简化的热力学循环
作者Author: Lecly 时间: 2025-12-23 11:58
谢谢付博士，了解了。
| 欢迎光临 计算化学公社 (http://bbs.keinsci.com/) |
Powered by Discuz! X3.3 |