Hamiltonian Replica Exchange Umbrella Sampling (HREUS)
Hamiltonian Replica Exchange Umbrella Sampling (HREUS) is implemented by utilizing MPI & biasmk4_struct. Compatibility with biasmk3_struct might be added later.
Background
The work flow of the replica exchange regarding the MPI communication is shown below.
ファイル:Replica exchange mpi.pptx
Module structure
The source file for HREUS is named as "reus.F".
<source lang="fortran">
MODULE reus_module
USE mpi_module,ONLY: mpi USE timer_module,ONLY: timer_struct USE FreeFlex_module,ONLY: biasmk4_ IMPLICIT NONE
TYPE reus_struct TYPE(timer_struct) :: timer INTEGER :: dim !!! Dimension of REUS, 1 or 2. REAL8,ALLOCATABLE :: data_recv(:) !!! Data from each node (1:mpi%isize) INTEGER,ALLOCATABLE :: iRep_list(:) !!! iRep for each node (1:mpi%isize) INTEGER,ALLOCATABLE :: pair_list(:,:) !!! Pair list for bias pairs. INTEGER :: n_pair !!! Number of pairs. INTEGER :: mpisize !!! Number of parallel processes.
LOGICAL :: calculate_reus !!! Switch for reus.
INTEGER :: n_select !!! Number of pairs exchanged each time.
INTEGER :: exchange_interval !!! Exchange interval.
INTEGER :: bias_set_id !!! ID of bias set applied to reus.
INTEGER,ALLOCATABLE :: pos_PARA(:) !!! Position of the PARAMETER of average.
!!! (1:dim)
INTEGER,ALLOCATABLE :: dim_q(:) !!! dimension of each q used in REUS.
!!! (1:dim)
!!! Asumming q(1)%a(:) for all q.
INTEGER :: expected_n_pair !!! Expected pair number.
INTEGER :: max_iter !!! maximum iteration for pair number determination.
LOGICAL,ALLOCATABLE :: pair_chosen(:) !!! bool for chosen flag.
INTEGER,ALLOCATABLE :: spcoord_id(:) !!! id of spcoord related to REUS. (1:dim)
INTEGER :: master_node !!! rank of the master node.
CONTAINS PROCEDURE :: init !!! (master_nml_unit) PROCEDURE,PRIVATE :: generate_pairlist !!! (biasmk4) PROCEDURE,PRIVATE :: gather_data !!! (q) PROCEDURE,PRIVATE :: scatter_id !!! (iRep) PROCEDURE,PRIVATE :: bcast_id !!! () PROCEDURE,PRIVATE :: exchange_id !!! () PROCEDURE,PRIVATE :: calculate_delta !!! (i_pair,delta) PROCEDURE :: reus_process !!!(coord,iRep) END TYPE reus_struct
CONTAINS
</source>
The structure of procudures is as follows.
- init(master_nml_unit)
- -----> generate_pairlist !!! ()
- reus_process !!!(data_send)
- -----> gather_data !!! (data_send)
- ---------> calculate_delta !!! (i_pair,delta)
- -------------> exchange_id !!! ()
- ------------------> scatter_id !!! (iRep)
Namelist variables
<source lang="fortran">
NAMELIST/REUS/ calculate_reus,dim,n_select,exchange_interval,bias_set_id,pos_PARA, &
expected_n_pair,dim_q,spcoord_id,max_iter,master_node
LOGICAL :: calculate_reus = .FALSE. !!! Switch for reus.
INTEGER :: n_select = 1 !!! Number of pairs exchanged each time.
INTEGER :: exchange_interval = 100 !!! Exchange interval.
INTEGER :: bias_set_id = 1 !!! ID of bias set applied to reus.
INTEGER :: dim = 1 !!! Dimension of REUS.
CHARACTER(LEN = 256) :: dim_q = "1" !!! Dimension of each spcoord q.
CHARACTER(LEN = 256) :: pos_PARA = "1" !!! Position of the parameter that describes the expected value of q.
!!! This information is used for replica pair generation.
CHARACTER(LEN = 256) :: spcoord_id = "1" !!! The id of spcoord used in REUS.
INTEGER :: expected_n_pair = 1 !!! Expected pair number. For 2D system, the expected pair number is around [3ij-2(i+j)+1]. INTEGER :: max_iter = 100 !!! maximum iter count for pair generation. INTEGER :: master_node = 0 !!! The mpirank of master node. LOGICAL :: last_id_file = "" !!! The output file storing last replica exchange id info.
</source>
Pair list generation
bias の平衡位置が近いものと pair を組むように、まず bias の位置関係でトラジェクトリ番号を sort し、 隣同士の bias を判定することで pairlist を生成する。また、DUMMYbiasが存在する場合には別処理を加えている。
以下にどこでpairlistを形成することを想定しているかを以下に書いてある。
・通常のbiasの場合 あるbiasを基準にして、そのbiasの上下左右とpairを形成する。
・DUMMYbiasの場合 あるDUMMYbias を基準にして同じ方向の bias とすべて pair を形成する。
Generating &biasdata for 2D REUS using biasdata_gen.exe
A standalone program called biasdata_gen.exe is provided for conveniently generating &biasdata section for 2D REUS simulation.
The usage is as follows:
./biasdata_gen.exe input.nml output_file
<source lang="fortran">
INTEGER :: nbiasA,nbiasB !!! Number of biases for each dimension. CHARACTER(LEN=256) :: secA1,secA2 !!! Each section for biasA set. CHARACTER(LEN=256) :: secB1,secB2 !!! Each section for biasB set. REAL8 :: startA !!! Start value for parameter A. REAL8 :: intervalA !!! Increment value for parameter A. REAL8 :: startB !!! Start value for parameter B. REAL8 :: intervalB !!! Increment value for parameter B.
NAMELIST/BIASDATA/ nbiasA,nbiasB,secA1,secA2,secB1,secB2,startA, &
intervalA,startB,intervalB
</source>
Example input.nml:
&BIASDATA nbiasA = 8 nbiasB = 4 secA1 = "PARABOLIC 1 1 2 1.d-10" secA2 = "" secB1 = "PARABOLIC 1 2 2 1.d-10" secB2 = "" startA = 12.5d-10 intervalA = 2.0d-10 startB = 3.5d-10 intervalB = 1.0d-10 /
The output for this nml:
&biasdata PARABOLIC 1 1 2 1.d-10 0.1250000000E-08 1 9 17 25 PARABOLIC 1 1 2 1.d-10 0.1450000000E-08 2 10 18 26 PARABOLIC 1 1 2 1.d-10 0.1650000000E-08 3 11 19 27 PARABOLIC 1 1 2 1.d-10 0.1850000000E-08 4 12 20 28 PARABOLIC 1 1 2 1.d-10 0.2050000000E-08 5 13 21 29 PARABOLIC 1 1 2 1.d-10 0.2250000000E-08 6 14 22 30 PARABOLIC 1 1 2 1.d-10 0.2450000000E-08 7 15 23 31 PARABOLIC 1 1 2 1.d-10 0.2650000000E-08 8 16 24 32 PARABOLIC 1 2 2 1.d-10 0.3500000000E-09 1 2 3 4 5 6 7 8 PARABOLIC 1 2 2 1.d-10 0.4500000000E-09 9 10 11 12 13 14 15 16 PARABOLIC 1 2 2 1.d-10 0.5500000000E-09 17 18 19 20 21 22 23 24 PARABOLIC 1 2 2 1.d-10 0.6500000000E-09 25 26 27 28 29 30 31 32 &end
Using monitor module REUS_qU for data output
The monitor module REUS_qU is intended to be used with REUS for outputting spcoord value and reus related potential at each interval.
The usage format in the main namelist is:
&MONITOR interval = 1 !!! The output interval. file = "reus_output" !!! The output file name. condition = "REUS_qU 1 2" !!! The "1 2" here are spcoord id used in this reus simulation. separator = "$" !!! Separator symbol. /
This will output the value of each spcoord and reus potential every interval per line.
Using WHAM to analyze data
The method for extrapolating the free energy surface from output files for REUS is WHAM(weighted histogram analysis method). Due to the spcoord dependency of biasmk4_system, WHAM now is provided as a utility of FreeFlex main program instead of a standalone program. In order to use WHAM for analysis, 2 namelists must be provided in the main namelist file. The /WHAMSETTING/ and /WHAM_INPUT/
To perform the analysis, simply prepare the namelist file correct spcoord and biasmk4_system settings and the above 2 namelist, and excute FreeFlex main program without using mpi mode.
Example:
./FreeFlex/FreeFlex.exe namelist
&WHAMSETTING namelist variables
<source lang="fortran"> ! in FreeFlex.F
LOGICAL :: wham_mode = .FALSE. !!! The switch for WHAM analysis. When it is true, normal MD simulation will not be performed in FreeFlex.exe. INTEGER :: mpi_size = 20 !!! The total proc number for the REUS simulation to be analyzed, must be provided correctly. !!! WHAM計算を実施するMDの並列トラジェクトリ（ノード）数。 !!! （WHAM計算自体は、シリアルで実行する） !!! （バイアスポテンシャルを設定する際に、予め必要となる。） NAMELIST/WHAMSETTING/ wham_mode,mpi_size
</source>
&WHAM_INPUT namelist variables
<source lang="fortran"> ! in wham_calc.F
LOGICAL :: wham_mode = .FALSE. !!! switch for WHAM mode.
LOGICAL :: reus_mode = .FALSE. !!! switch for REUS mode.
!!! If the simulation is normal umbrella sampling, it has to be turned off.
!!! レプリカ交換をするときとしないときで、出力ファイルの形式が異なるため。
INTEGER :: nbin = 100 !!! Number of bins. !!! nbin --- WHAM解析の座標のbin数。2次元（以上）のときには、各次元ごとに nbin個とられる。
CHARACTER(LEN=256) :: nbin_multi = "" !!! nbin_multi --- 多次元面に対するWHAM解析の座標のbin数。 !!! 多次元面に対して異なるbinの個数を取りたい場合に指定する。例) nbin_multi="50 200 200"。 !!! nbinと同時に指定された場合、nbin_multiの値が優先される。
INTEGER :: ndir = 1 !!! Number of dirs (same as number of procs). !!! ndir --- 解析される並列MDトラジェクトリ（ノード）数。 !!! 通常は、&WHAMSETTING 中の mpi_size と同じ数。
INTEGER :: bias_set_ID = 1 !!! bias_set_id used for reus. !!! &biasdata が入力ファイル中に複数ある場合、WHAM計算で扱うバイアス・ポテンシャルを表す&biasdata の番号。 !!! 入力ファイル中で、&biasdata を上から数える番号として与える。 !!! WHAM計算で使用する&biasdata 中には、WHAMに現れないバイアス・ポテンシャルは指定しないこと。 !!! （仮に指定しても、WHAMの中では無視される。） それ以外のバイアスは、別途 &biasdata に分けてかく。
LOGICAL :: show = .TRUE. !!! Show distribution graph for each state. !!! show --- MDトラジェクトリから生成される生のヒストグラムを表示するフラグ。デバック用。
INTEGER :: nprof_column = 1 !!! The column number that stores the replica id in a reus simulation. !!! nprof_column --- reus_mode=.TRUE.のとき、入力ファイル中でレプリカIDを記述する列の番号。
INTEGER :: ndim = 1 !!! Dimension of the simulation.
CHARACTER(LEN=256) :: data_column = "1" !!! The column number that stores spcoord value for each dimension.
CHARACTER(LEN=256) :: spcoord_ID = "1" !!! Spcoord id for each dimension.
CHARACTER(LEN=256) :: xmin_, xmax_ !!! Designated upper/lower bound of the histogram with respect to the reaction coordinates.
CHARACTER(LEN=256) :: pos_para = "2" !!! Position of parameter specific to reus coordinate.
!!! Specifies the corresponding column number in data_column as
!!! the pos_para in &REUS(For example, if the
!!! dimension of spcoords are all 1, then pos_para is the
!!! same as data_column, if the dimension of spcoords
!!! are larger than 1, then pos_para is a subset of
!!! data_column specifying only columns that correspond
!!! to pos_para specified in &REUS).
LOGICAL :: diis_mode = .true. !!! DIISによる収束加速を用いるかどうかのフラグ（default=.TRUE.） INTEGER :: ndiis = 8 !!! DIISを用いるとき、過去何ステップの残差を保存するか
LOGICAL :: guess_use = .true. !!! 初期値の推定を行うかどうかのフラグ（default=.TRUE.）
LOGICAL :: error_calc = .false. !!! 自由エネルギー面の標準偏差を計算するかどうかのフラグ（default=.FALSE.） INTEGER :: nblock = 5 !!! 標準偏差を計算する際のサンプルの分割数
!!! data_column --- WHAMの入力ファイル中で、バイアスにかかわる一般化座標の列。 !!! 計算する自由エネルギー面にかかわる座標が多次元のspcoordの場合、注意！ !!! その多次元の分をすべて指定する必要がある。 !!! !!! out_column --- 反応座標として出力する座標列。 !!! data_columnの中から選択し、選択されなかった列については縮約される。 !!! defaultではout_column=data_column（どの列も縮約しない）。 !!! !!! spcoord_ID --- 自由エネルギー面の座標のspcoord番号。 !!! MD計算において定義したspcoordのシリアル番号として与える。 !!! data_column の要素の個数とおなじだけ与える必要がある。 !!! !!! xmin_, xmax_ --- WHAMにおけるヒストグラムの範囲を最小値と最大値で指定する。 !!! 多次元のWHAMを行う場合には xmin_=xmin_(1), xmin(2), … のように全ての !!! 次元に対して指定する。デフォルトではサンプリングした領域から自動的に決定される。 !!!
!!! （例１） 1次元座標q1(1), 3次元座標q2(1:3) で、q1(1)とq2(1)の2次元面を作る場合。 !!! (ただし、q1, q2のspcoord_IDは、1, 2とする） !!! また、nprof_column=1で1列目はノード番号に使用されていて、2列目からはじまるとする。） !!! nprof_column=1, data_column="2 3 4 5", spcoord_ID="1 2 2 2", out_column="2 3" !!! !!! （例２） 3次元座標q1(1:3)で、q1(1)とq1(3)の2次元面を作る場合。 !!! (ただし、q1のspcoord_IDは、10とする） !!! data_column="1 2 3", spcoord_ID="10 10 10", out_column="1 3"
INTEGER :: interval = 1 !!! Subsample interval. !!! interval --- 読み込み時の行間隔。interval行ごとに読み込んだWHAM計算を行う。
CHARACTER(len=256) :: free_energy_output_file_name !!! Output file name for free energy. !!! 自由エネルギー面の出力ファイルは、Gnuplot に即した形式で与えられる。
CHARACTER(LEN=256) :: input_file_name = "whamz.output" !!! Output file name from monitor module REUS_qU. !!! WHAMの入力データのファイルは、ndir個のディレクトリで同じファイル名とする。 !!! すなわち、input_file_name='whamz.output', ndir=20ならば !!! 0000/whamz.output, 0001/whamz.output, ..., 0019/whamz.output !!! となる。
CHARACTER(len=256) :: histogram_output_file_name !!! Output file name for histogram data. INTEGER :: nskip = 50 !!! Lines to be skipped from the start of the data. INTEGER :: nread = 1000000 !!! Lines to be read after the initial skip.
INTEGER :: ncore !!! Number of cores in total. Must be given if reus_mode. CHARACTER(len=256) :: file_tail !!! Tail of input file if any. INTEGER :: max_file_No !!! The maximum file number to be read in the mpi directory.
NAMELIST/WHAM_INPUT/ wham_mode,ndir,nbin,nskip,nread,spcoord_ID,input_file_name, &
free_energy_output_file_name,histogram_output_file_name,show,bias_set_ID,reus_mode, &
interval,data_column,nprof_column,pos_para,ncore,file_tail
</source>
Example output file to be analyzed by WHAM:
Module structure of wham_calc_struct. <source lang="fortran">
MODULE wham_module
USE graph_module,ONLY: graph_struct IMPLICIT NONE
INTEGER,PARAMETER :: maxloop = 10000
TYPE,EXTENDS(graph_struct) :: wham_struct REAL8 :: T INTEGER,ALLOCATABLE :: spcoord_ID(:) !!! ID of spcoord related to REUS INTEGER,ALLOCATABLE :: pos_para(:) !!! Position of coordinate_specific colomn number. INTEGER :: dim !!! Dimension of analysis.
CONTAINS PROCEDURE :: wham_calc !!! (coord,bias,master_nml_unit,T) called from main program. PROCEDURE :: calculate !!! (hist,boltz,T) !!! WHAM calculation !!! all hist(:) must have the same xmin, xmax, and dx. PROCEDURE :: getboltz !!! (hist,bias,list,boltz)
END TYPE wham_struct
CONTAINS
</source>
Example of error calculations
&WHAM_INPUT 内の error_calc フラグを .true. にすると、誤差解析のために標準偏差の計算を行うことができる。
サンプルデータを nblock の数だけ分割してそれぞれのブロックで自由エネルギー計算を行い、標準偏差を求める。
nblock=5 (default) の場合、 Total, Block 1, Block2, …, Block5 の順にWHAM計算が行われ、計算が正常終了すると以下の出力が得られる。
また、gnuplot により Total および各 Block の自由エネルギー面が下図のように表示される。
各自由エネルギー面同士の相対的な高さは標準偏差を最小化するように決められる。
このとき、出力ファイルの2行目は Total、3行目から7行目までが順に Block 1, Block2, …, Block5 となる。
8行目は Block から計算された標準偏差である。
gnuplot上で以下のコマンドを実行すれば、エラーバー付きの自由エネルギー曲線が描画される。
このとき、自由エネルギー面の出力ファイルは "free_energy.output" とした。nblock を変更した際や、2次元以上の自由エネルギー計算を 行った際には標準偏差が書き出される行が変化するため注意すること。
> plot "free_energy.output" u 1:2:8 w errorbar
Convergence acceleration by DIIS method
本プログラムではWHAM計算に対して、連立方程式のSCF解法における一般的な収束加速法であるDIIS法を採用している。 ただし元データの質が悪いときには、しばしばDIIS法を用いると収束性が著しく悪化する場合がある（原因は不明）。 このような場合に、diis_mode = .false. で計算を行うと改善することがある。
Details on choosing pair for exchange
The general rule for REUS exchange process is to ensure the detailed balance of the system. This enforces choosing every possible replica pairs with equal probability. This can be done by the following procedure.
- For each exchange interval reached, randomly choose N pairs from the pair list and perform N exchanges.
Then we have 3 options as following:- 1.Exchanges include repeated pairs and repeated nodes in different pairs.
- 2.Exchanges exclude repeated pairs but include repeated nodes in different pairs.
- 3.Exchanges exclude repeated nodes in different pairs (thus excluding repeated pairs as well).