本帖最后由 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。这里我们准备了一个简单的水盒子，标记了其中的一个水分子：
相对于平衡模拟，在配置文件（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文件，我们可以画出沿λ的自由能变化：
当然我们一般只关心λ=0和λ=1时的自由能差值，这里ΔG=6.2 kcal/mol，与实验值-6.3 kcal/mol相符。
这个例子中用到的所有文件如下，其他类型的alchemical transformation道理相同：
参考文献：J. Phys. Chem. Lett. 2025, 16, 4419–4427，Acc. Chem. Res. 2025, DOI:10.1021/acs.accounts.5c00666
| 
发表于 Post on 2025-9-27 22:58:01
|
只看该作者 Only view this author
|只看大图 View big picture
|倒序浏览 Reverse view
|阅读模式 Reading model
收藏 Add to favorites
楼主 Author