gmx wham [-ix [<.dat>]] [-if [<.dat>]] [-it [<.dat>]] [-ip [<.dat>]]
[-is [<.dat>]] [-o [<.xvg>]] [-hist [<.xvg>]] [-oiact [<.xvg>]]
[-iiact [<.dat>]] [-bsres [<.xvg>]] [-bsprof [<.xvg>]]
[-tab [<.dat>]] [-nice ] [-xvg ] [-min ]
[-max ] [-[no]auto] [-bins ] [-temp ] [-tol ]
[-[no]v] [-b ] [-e ] [-dt ] [-[no]histonly]
[-[no]boundsonly] [-[no]log] [-unit ] [-zprof0 ]
[-[no]cycl] [-[no]sym] [-[no]ac] [-acsig ]
[-ac-trestart ] [-nBootstrap ] [-bs-method ]
[-bs-tau ] [-bs-seed ] [-histbs-block ] [-[no]vbs]
gmx wham一个用于实现加权直方图分析方法(WHAM, Weighted Histogram Analysis Method)的分析程序, 用于分析伞形抽样模拟的输出文件以计算平均力势(PMF, potential of mean force).
目前此程序支持三种输入模式:
使用选项-it, 用户需要提供一个文件, 其中包含伞形抽样模拟的运行输入文件(.tpr文件)的文件名, 还要 使用选项-ix提供另一个文件, 其中包含pullx mdrun输出文件的文件名. .tpr文件和pullx文件的顺序必须对应, 即第一个.tpr文件生成了第一个pullx文件, 并以此类推.
除用户使用选项-if提供牵引力输出文件名称(pullf.xvg)以外, 与前述输入模式相同. 伞形势中的位置由牵引力计算得到. 无法用于表格形式的伞形势.
使用选项-ip, 用户需提供(gzip压缩的).pdo文件的文件名, 即GROMACS 3.3的伞形输出文件. 如果你使用某些特殊的反应坐标, 你也可以生成自己的.pdo文件并使用-ip选项将其提供给gmx wham. .pdo文件的文件头必须类似于以下形式:
# UMBRELLA 3.0
# Component selection: 0 0 1
# nSkip 1
# Ref. Group 'TestAtom'
# Nr. of pull groups 2
# Group 1 'GR1' Umb. Pos. 5.0 Umb. Cons. 1000.0
# Group 2 'GR2' Umb. Pos. 2.0 Umb. Cons. 500.0
#####
牵引组(pull group)的个数, 伞形势位置(umbrella position), 力常数(force constant)和名称(当然)都可以不同. 文件头以下, 需要为每个牵引组提供一个时间列和一个数据列(即相对于伞形势中心的位移). 目前每个.pdo文件最多可以包含四个牵引组.
默认情况下, 在WHAM中会使用所有pullx/pullf文件中找到的所有牵引组. 如果只使用其中的某些牵引组, 用户可以提供一个牵引组选择文件(使用选项-is). 选择文件必须为tpr-files.dat中的每个.tpr文件提供一行说明, 内容必须为相应于.tpr文件中每个牵引组的一位数字(0或1). 在这里1表示该牵引组会在WHAM中使用, 0表示忽略. 例如, 如果你有3个.tpr文件, 每个包含4个牵引组, 但只使用牵引组1和2, 则groupsel.dat文件内容如下:
默认情况下, 输出文件有:
-o: PMF输出文件
-hist: 直方图输出文件
请注意, 始终要检查直方图是否充分重叠.
程序假定伞形势为简谐势, 力常数从.tpr或.pdo文件中读取. 如果使用了非简谐的伞形力, 可以用-tab提供一个表格式的势能函数.
WHAM选项
-bins: 分析中使用的分格数
-temp: 模拟温度
-tol: 剖面(概率)的变化小于所给容差时停止迭代
-auto: 自动决定边界
-min, -max: 剖面的边界
可以使用选项-b, -e和-dt筛选用于计算剖面的数据点. 调整-b以保证每个伞形窗口都达到充分平衡.
使用选项-log时(默认), 会以能量单位输出剖面, 否则(使用-nolog选项)会输出概率. 可以使用选项-unit指定单位. 以能量单位输出时, 第一分格中的能量定义为零. 如果你希望其他位置自由能为零, 可以设置-zprof0(在使用自展法时很有用, 见下文).
对于环形或周期性的反应坐标(二面角, 无渗透梯度的通道PMF), 选项-cycl很有用. gmx wham会利用体系的周期性, 生成一个周期性的PMF. 反应坐标的第一个和最后一个分格会被假定为相邻.
使用选项-sym时, 在输出前会使剖面关于z=0对称, 在某些情况下, 如用于膜体系时很有用.
自相关
使用-ac选项时, gmx wham会估计每个伞形窗口的积分自相关时间(IACT, integrated autocorrelation time)τ, 并使用1/[1+2*τ/dt]作为各个窗口的权重. IACT会写入由选项-oiact指定的文件中. 在冗长(verbose)输出模式下, 所有自相关函数(ACF, autocorrelation functions)都会写入hist_autocorr.xvg文件. 由于在采样不足的情况下可能会严重低估IACT, 利用-acsig选项, 用户可使用高斯函数沿反应坐标对IACT进行平滑(高斯函数的σ由-acsig提供, 见iact.xvg中的输出). 注意, 程序使用简单的积分方法估计IACT, 且只考虑大于0.05的ACF. 如果你想使用更复杂(但可能不那么稳健)的方法, 比如拟合到双指数函数, 来计算IACT, 你可以使用gmx analyze来计算IACT, 并通过iact-in.dat文件(选项-iiact)将其提供给gmx wham. 在这个文件中每个输入文件(.pdo或pullx/f文件)对应一行, 各输入文件的每个牵引组对应一列.
误差分析
可以使用自展分析(bootstrap analysis)来估计统计误差. 请小心使用, 否则实质上可能会低估统计误差. 自展技术的更多背景知识和例子可以在Hub, de Groot and Van der Spoel, JCTC (2010) 6: 3713-3720中找到.
-nBootstrap定义自展的个数(比如使用100). 本程序支持四种自展方法, 通过-bs-method进行选择.
(1) b-hist 默认方法: 将完整的直方图视为独立的数据点, 给直方图赋予随机权重来实现自展(“贝叶斯自展”). 注意, 沿反应坐标轴上的每个点都必须被多个独立直方图所覆盖(比如10个直方图), 否则会低估统计误差.
(2) hist: 将完整的直方图视为独立的数据点. 对每个自展, 从给定的N个直方图中随机选取N个直方图(允许重复, 即, 放回抽样). 为避免沿反应坐标轴上无数据的空隙, 可以定义直方图块(-histbs-block). 在那种情况下, 会将给定的直方图划分为块, 只有各块内部的直方图才会混合. 注意, 每块内的直方图必须能代表所有可能出现的直方图, 否则会低估统计误差.
(3) traj: 用给定的直方图产生新的随机轨迹, 这样产生的数据点遵从给定直方图的分布, 并具有适当的自相关. 每个窗口的自相关时间(ACT)必须是已知的, 所以要使用-ac选项或者利用-iiact手动提供ACT. 如果所有窗口的ACT都相同(并且已知), 你也可以用-bs-tau提供ACT. 注意, 在采样不足的情况下, 即如果各个直方图在各自的位置不能代表整个相空间, 此方法可能严重低估误差.
(4) traj-gauss: 与traj方法相同, 但轨迹不是根据伞形直方图自展得到, 而是从均值和宽度与伞形直方图相同的高斯函数得到. 此方法给出的误差估计类似于traj方法.
自展法的输出:
-bsres: 平均剖面和标准偏差
-bsprof: 所有自展剖面
使用-vbs选项(冗长自展)会输出每个自展使用的直方图, 并且, 使用traj自展方法时, 还会输出直方图的累积分布函数.
输入/输出文件选项| 选项 | 默认值 | 类型 | 说明 |
|---|
-ix [<.dat>] | pullx-files.dat | 输入, 可选 | 通用数据文件 |
-if [<.dat>] | pullf-files.dat | 输入, 可选 | 通用数据文件 |
-it [<.dat>] | tpr-files.dat | 输入, 可选 | 通用数据文件 |
-ip [<.dat>] | pdo-files.dat | 输入, 可选 | 通用数据文件 |
-is [<.dat>] | groupsel.dat | 输入, 可选 | 通用数据文件 |
-o [<.xvg>] | profile.xvg | 输出 | xvgr/xmgr 文件 |
-hist [<.xvg>] | histo.xvg | 输出 | xvgr/xmgr文件 |
-oiact [<.xvg>] | iact.xvg | 输出, 可选 | xvgr/xmgr文件 |
-iiact [<.dat>] | iact-in.dat | 输入, 可选 | 通用数据文件 |
-bsres [<.xvg>] | bsResult.xvg | 输出, 可选 | xvgr/xmgr文件 |
-bsprof [<.xvg>] | bsProfs.xvg | 输出, 可选 | xvgr/xmgr文件 |
-tab [<.dat>] | umb-pot.dat | 输入, 可选 | 通用数据文件 |
控制选项| 选项 | 默认值 | 说明 |
|---|
-nice <int> | 19 | 设置优先级 |
-xvg <enum> | xmgrace | xvg绘图格式: xmgrace, xmgr, none |
-min <real> | 0 | 剖面的最小坐标 |
-max <real> | 0 | 剖面的最大坐标 |
-[no]auto | yes | 自动确定min和max |
-bins <int> | 200 | 剖面使用的分格数 |
-temp <real> | 298 | 温度 |
-tol <real> | 1e-06 | 容差 |
-[no]v | no | 冗长模式 |
-b <real> | 50 | 要分析的起始时间(ps) |
-e <real> | 1e+20 | 要分析的终止时间(ps) |
-dt <real> | 0 | 每隔dt ps分析一次 |
-[no]histonly | no | 输出直方图即退出 |
-[no]boundsonly | no | 确定min和max后即退出(使用-auto选项) |
-[no]log | yes | 计算剖面的对数值后再打印 |
-unit <enum> | kJ | 对数输出时的能量单位: kJ, kCal, kT |
-zprof0 <real> | 0 | 将此位置的剖面值定义为0.0(使用-log选项) |
-[no]cycl | no | 产生环形/周期性的剖面. 假定min和max是同一点. |
-[no]sym | no | 使剖面关于z=0对称 |
-[no]ac | no | 计算积分自相关时间并在wham中使用 |
-acsig <real> | 0 | 沿反应坐标轴, 使用以选项值为σ的高斯函数对自相关时间进行平滑 |
-ac-trestart <real> | 1 | 计算自相关函数时, 不同数据点之间的间隔ps数 |
-nBootstrap <int> | 0 | 自展的个数, 用以估计统计不确定性(如, 200) |
-bs-method <enum> | b-hist | 自展方法: b-hist, hist, traj, traj-gauss |
-bs-tau <real> | 0 | 假定适用于所有直方图的自相关时间(ACT). 如果ACT未知, 使用选项-ac |
-bs-seed <int> | -1 | 自展的种子(-1表示使用时间值) |
-histbs-block <int> | 8 | 混合直方图时, 仅混合在-histbs-block指定块内的直方图. |
-[no]vbs | no | 冗长自展. 输出每个自展的CDF和直方图文件. |
文章链接：GROMACS各类程序(名称排序)|Jerkwin
如有侵权联系我，我将删除
本文目的只为宣传使用
