gmx bar利用Bennett接受比率方法计算自由能差的估计值

gmx bar [-f [<.xvg> [...]]] [-g [<.edr> [...]]] [-o [<.xvg>]] [-oi [<.xvg>]]

[-oh [<.xvg>]] [-nice ] [-[no]w] [-xvg ] [-b ]

[-e ] [-temp ] [-prec ] [-nbmin ]

[-nbmax ] [-nbin ] [-[no]extp]

gmx bar通过Bennett接受率方法(BAR, Bennett’s acceptance ratio)估计自由能差值, 也可以自动将由BAR得到的一系列分立自由能进行组合得到自由能估计值.

每个分立的BAR自由能差值依赖于两个不同状态的模拟, 且称为态A和态B, 它们由参数λ控制(见.mdp参数init_lambda). 给定态A, BAR方法可以计算态B相对于态A的哈密顿差的加权平均的比率, 反之亦然. 相对于另一状态的能量差在模拟中必须显式地计算, 可以通过.mdp选项foreign_lambda实现.

输入选项-f需要读入多个dhdl.xvg文件, 支持两种输入文件类型:

包含多于一个y值的文件, 其中应包含dH/dλ和Δλ的列. λ的值根据列标题推定: 模拟使用的λ根据dH/dλ列的标题推定, 外部λ值根据Delta H列的标题推定.

仅有一个y值的文件. 对这些文件应使用-extp选项, 并假定y值为dH/dλ, 而哈密顿量与λ呈线性关系. 模拟的λ值根据子标题(如果存在)推定, 否则会根据子目录下文件名中的数字推定.

模拟的λ值根据dhdl.xvg文件中包含字符串dH的列标题解析得出, 外部λ值根据包含大写字母D和H的列标题解析得出, 温度根据包含T=的标题解析得出.

输入选项-g需要读入多个.edr文件, 它们可以包含能量差列表(见.mdp选项separate_dhdl_file), 或者一系列直方图(见.mdp选项dh_hist_size和dh_hist_spacing). 程序会自动从ener.edr文件中推断出温度和λ值.

除了.mdp的foreign_lambda选项外, 也可以根据dH/dλ值外推得到能量差. 这可通过-extp选项实现, 它假定系统的哈密顿量与λ呈线性关系, 虽然通常并非如此.

自由能估计由使用二分法的BAR方法确定, 输出精度由-prec设定. 误差估计考虑了时间相关, 这是通过将数据分块, 并假定这些分块之间互相独立, 计算它们之间的自由能差来实现的. 最终的误差估计由5个分块的平均方差决定. 用于误差估计的分块数可以通过选项-nbmin和-nbmax来指定.

gmx bar会尝试合计具有相同’本地’和’外部’λ值的样本, 但总会假定样本互相独立. 注意, 当合计具有不同采样间隔的能量差或能量导数时, 这个假定几乎肯定是不正确的. 连续的能量通常是相关的, 不同的时间间隔意味着样本间的相关度不同.

结果分为两部分: 后一部分包含了以kJ/mol为单位的最终结果, 以及每一部分和总体的误差估计. 前一部分包含了详细的自由能差估计和相空间重叠量度, 以kT为单位(以及它们的误差估计). 打印出的值为:

lam_A: A点的λ值.

lam_B: B点的λ值.

DG: 自由能估计.

s_A: B在A中的相对熵估计.

s_B: A在B中的相对熵估计.

stdev: 每个样本标准偏差的估计期望

两个状态在彼此系综内的相对熵可以理解为相空间重叠的量度: lambda_B的工作样本在lambda_A系综内的相对熵s_A(对s_B反之亦然), 是两个状态Boltzmann分布之间’距离’的量度, 当分布相同时, 其值为0. 详见Wu & Kofke, J. Chem. Phys. 123 084109 (2005).

每个样本标准偏差的估计期望, 见Bennett BAR方法的原始论文 Bennett, J. Comp. Phys. 22, p 245 (1976). 其中的Eq. 10给出了采样质量的估计(并非直接的实际统计误差, 因为它假定了样本相互独立).

要得到相空间重叠估计的可视化结果, 可使用-oh选项及-nbin选项输出一系列直方图.