Action: WHAM

Module	wham
Description	Usage
Calculate the weights for configurations using the weighted histogram analysis method.
output value	type
the vector of WHAM weights to use for reweighting the elements in a time series	vector

Input

The arguments that serve as the input for this action are specified using one or more of the keywords in the following table.

Further details and examples

Calculate the weights for configurations using the weighted histogram analysis method.

The example input below shows how this command is used.

Click on the labels of the actions for more information on what each action computes

#SETTINGS NREPLICAS=4
phi: TORSIONCalculate one or multiple torsional angles. More details ATOMSthe four atoms involved in the torsional angle=5,7,9,15
The TORSION action with label phi calculates the following quantities:
 Quantity   
 Description  
phi.value
the TORSION involving these atomspsi: TORSIONCalculate one or multiple torsional angles. More details ATOMSthe four atoms involved in the torsional angle=7,9,15,17
The TORSION action with label psi calculates the following quantities:
 Quantity   
 Description  
psi.value
the TORSION involving these atomsrp: RESTRAINTAdds harmonic and/or linear restraints on one or more variables. More details ARGthe values the harmonic restraint acts upon=phi KAPPA specifies that the restraint is harmonic and what the values of the force constants on each of the variables are=50.0 ATthe position of the restraint=@replicas:This keyword specifies that different replicas have different values for this quantity.  See here for more details.{-3.00,-1.45,0.10,1.65}
# Get the bias on each of the replicas
The RESTRAINT action with label rp calculates the following quantities:
 Quantity   
 Description  
rp.bias
the instantaneous value of the bias potential
rp.force2
the instantaneous value of the squared force due to this bias potentialrep: GATHER_REPLICASCreate a vector that contains the copies of the input quantities from all replicas More details ARGthe argument from the various replicas that you would like to gather=rp.bias
# Merge the biases on each of the replicas into a single vector
The GATHER_REPLICAS action with label rep calculates the following quantities:
 Quantity   
 Description  
rep.rep
the input arguments for each of the replicasall: CONCATENATEJoin vectors or matrices together More details ARGthe values that should be concatenated together to form the output vector=rep.*
# Collect all the bias values
The CONCATENATE action with label all calculates the following quantities:
 Quantity   
 Description  
all.value
the concatenated vector/matrix that was constructed from the input valuescol: COLLECTCollect data from the trajectory for later analysis More details TYPE required if you are collecting an object with rank>0=vector ARGthe label of the value whose time series is being stored for later analysis=all STRIDE the frequency with which the data should be collected and added to the quantity being averaged=1
The COLLECT action with label col calculates the following quantities:
 Quantity   
 Description  
col.value
the time series for the input quantitywham: WHAMCalculate the weights for configurations using the weighted histogram analysis method. More details ARGthe stored values for the bias=col TEMPthe system temperature=300
The WHAM action with label wham calculates the following quantities:
 Quantity   
 Description  
wham.value
the vector of WHAM weights to use for reweighting the elements in a time seriesDUMPVECTORPrint a vector to a file More details ARGthe labels of vectors/matrices that should be output in the file=wham FILE the file on which to write the vetors=wham_data

As illustrated in the example above this command is used when you have run $N$ trajectories each of which was computed by integrating a different biased Hamiltonian. The input above calculates the probability of observing the set of configurations during the $N$ trajectories that we ran using the product of multinomial distributions using the formula shown below:

$$P( \underline{T} ) \propto \prod_{j=1}^M \prod_{k=1}^N (c_k w_{kj} p_j)^{t_{kj}}$$

In this expression the second product runs over the biases that were used when calculating the $N$ trajectories. The first product then runs over the $M$ bins in our histogram. The $p_j$ variable that is inside this product is the quantity we wish to extract; namely, the unbiased probability of having a set of CV values that lie within the range for the $j$th bin.

The quantity that we can easily extract from our simulations, $t_{kj}$, however, measures the number of frames from trajectory $k$ that are inside the $j$th bin. To interpret this quantity we must consider the bias that acts on each of the replicas so the $w_{kj}$ term is introduced. This quantity is calculated as:

$$w_{kj} = e^{+\beta V_k(x_j)}$$

and is essentially the factor that we have to multiply the unbiased probability of being in the bin by in order to get the probability that we would be inside this same bin in the $k$th of our biased simulations. Obviously, these $w_{kj}$ values depend on the value that the CVs take and also on the particular trajectory that we are investigating all of which, remember, have different simulation biases. Finally, $c_k$ is a free parameter that ensures that, for each $k$, the biased probability is normalized.

We can use the equation for the probability that was given above to find a set of values for $p_j$ that maximizes the likelihood of observing the trajectory. This constrained optimization must be performed using a set of Lagrange multipliers, $\lambda_k$, that ensure that each of the biased probability distributions are normalized, that is $\sum_j c_kw_{kj}p_j=1$. Furthermore, as the problem is made easier if the quantity in the equation above is replaced by its logarithm we actually chose to minimize

$$L = \sum_{j=1}^M \sum_{k=1}^N t_{kj} \ln c_k w_{kj} p_j + \sum_k\lambda_k \left( \sum_{j=1}^N c_k w_{kj} p_j - 1 \right)$$

After some manipulations the following (WHAM) equations emerge:

$$\begin{aligned} p_j & \propto \frac{\sum_{k=1}^N t_{kj}}{\sum_k c_k w_{kj}} \\ c_k & =\frac{1}{\sum_{j=1}^M w_{kj} p_j} \end{aligned}$$

which can be solved by computing the $p_j$ values using the first of the two equations above with an initial guess for the $c_k$ values and by then refining these $p_j$ values using the $c_k$ values that are obtained by inserting the $p_j$ values obtained into the second of the two equations above.

Notice that only $\sum_k t_{kj}$, which is the total number of configurations from all the replicas that enter the $j$th bin, enters the WHAM equations above. There is thus no need to record which replica generated each of the frames. One can thus simply gather the trajectories from all the replicas together at the outset. This observation is important as it is the basis of the binless formulation of WHAM that is implemented within PLUMED.

Syntax

The following table describes the keywords and options that can be used with this action