|
gmx wham伞形抽样后进行加权直方分析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文件内容如下: 1 1 0 0 1 1 0 0 1 1 0 0 默认情况下, 输出文件有: -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自展方法时, 还会输出直方图的累积分布函数.
|
||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||