Skip to content

Commit

Permalink
Added new style documentation for wham module and fixed typo in last …
Browse files Browse the repository at this point in the history
…commit
  • Loading branch information
Gareth Aneurin Tribello authored and Gareth Aneurin Tribello committed Nov 10, 2024
1 parent 9a063f2 commit 9a50c3f
Show file tree
Hide file tree
Showing 4 changed files with 78 additions and 50 deletions.
2 changes: 1 addition & 1 deletion src/generic/DumpProjections.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -116,7 +116,7 @@ DumpProjections::DumpProjections(const ActionOptions&ao):
checkRead();

for(unsigned i=0; i<getNumberOfArguments(); ++i) {
if( getPntrToAction(i)->getRank()!=0 ) error("cannot use DUMPPROJECTIONS with actions that are not scalars");
if( getPntrToArgument(i)->getRank()!=0 ) error("cannot use DUMPPROJECTIONS with actions that are not scalars");
(getPntrToArgument(i)->getPntrToAction())->turnOnDerivatives();
}
}
Expand Down
88 changes: 59 additions & 29 deletions src/wham/Wham.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -29,48 +29,78 @@
/*
Calculate the weights for configurations using the weighted histogram analysis method.
Suppose that you have run multiple \f$N\f$ trajectories each of which was computed by integrating a different biased Hamiltonian. We can calculate the probability of observing
the set of configurations during the \f$N\f$ trajectories that we ran using the product of multinomial distributions shown below:
\f[
The example input below shows how this command is used.
```plumed
#SETTINGS NREPLICAS=4
phi: TORSION ATOMS=5,7,9,15
psi: TORSION ATOMS=7,9,15,17
rp: RESTRAINT ARG=phi KAPPA=50.0 ...
AT=@replicas:{
-3.00000000000000000000
-1.45161290322580645168
.09677419354838709664
1.64516129032258064496
}
...
# Get the bias on each of the replicas
rep: GATHER_REPLICAS ARG=rp.bias
# Merge the biases on each of the replicas into a single vector
all: CONCATENATE ARG=rep.*
# Collect all the bias values
col: COLLECT TYPE=vector ARG=all STRIDE=1
wham: WHAM ARG=col TEMP=300
DUMPVECTOR ARG=wham FILE=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( \vec{T} ) \propto \prod_{j=1}^M \prod_{k=1}^N (c_k w_{kj} p_j)^{t_{kj}}
\label{eqn:wham1}
\f]
In this expression the second product runs over the biases that were used when calculating the \f$N\f$ trajectories. The first product then runs over the
\f$M\f$ bins in our histogram. The \f$p_j\f$ 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 \f$j\f$th bin.
The quantity that we can easily extract from our simulations, \f$t_{kj}\f$, however, measures the number of frames from trajectory \f$k\f$ that are inside the \f$j\f$th bin.
To interpret this quantity we must consider the bias that acts on each of the replicas so the \f$w_{kj}\f$ term is introduced. This quantity is calculated as:
\f[
w_{kj} =
\f]
$$
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 \f$k\f$th of our biased simulations. Obviously, these \f$w_{kj}\f$ 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, \f$c_k\f$ is a free parameter that ensures that, for each \f$k\f$, the biased probability is normalized.
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 \f$p_j\f$ that maximizes the likelihood of observing the trajectory.
This constrained optimization must be performed using a set of Lagrange multipliers, \f$\lambda_k\f$, that ensure that each of the biased probability distributions
are normalized, that is \f$\sum_j c_kw_{kj}p_j=1\f$. Furthermore, as the problem is made easier if the quantity in the equation above is replaced by its logarithm
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
\f[
\mathcal{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)
\f]
$$
\mathcal{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:
\f[
$$
\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}
\f]
which can be solved by computing the \f$p_j\f$ values using the first of the two equations above with an initial guess for the \f$c_k\f$ values and by then refining
these \f$p_j\f$ values using the \f$c_k\f$ values that are obtained by inserting the \f$p_j\f$ values obtained into the second of the two equations above.
$$
Notice that only \f$\sum_k t_{kj}\f$, which is the total number of configurations from all the replicas that enter the \f$j\f$th bin, enters the WHAM equations above.
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.
\par Examples
*/
//+ENDPLUMEDOC

Expand Down
20 changes: 9 additions & 11 deletions src/wham/WhamHistogram.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -31,13 +31,8 @@ namespace wham {
/*
This can be used to output the a histogram using the weighted histogram technique
This shortcut action allows you to calculate a histogram using the weighted histogram
analysis technique. For more detail on how this the weights for configurations are
computed see \ref REWEIGHT_WHAM
\par Examples
The following input can be used to analyze the output from a series of umbrella sampling calculations.
This shortcut action allows you to calculate a histogram using the weighted histogram analysis technique.
The following input illustrates how this is used in practise to analyze the output from a series of umbrella sampling calculations.
The trajectory from each of the simulations run with the different biases should be concatenated into a
single trajectory before running the following analysis script on the concatenated trajectory using PLUMED
driver. The umbrella sampling simulations that will be analyzed using the script below applied a harmonic
Expand All @@ -47,7 +42,7 @@ to determine a weight for each configuration in the concatenated trajectory. A
the configurations visited and their weights. This histogram is then converted into a free energy surface and output
to a file called fes.dat
\plumedfile
```plumed
#SETTINGS NREPLICAS=4
phi: TORSION ATOMS=5,7,9,15
psi: TORSION ATOMS=7,9,15,17
Expand All @@ -63,13 +58,16 @@ rp: RESTRAINT ARG=phi KAPPA=50.0 ...
hh: WHAM_HISTOGRAM ARG=phi BIAS=rp.bias TEMP=300 GRID_MIN=-pi GRID_MAX=pi GRID_BIN=50
fes: CONVERT_TO_FES GRID=hh TEMP=300
DUMPGRID GRID=fes FILE=fes.dat
\endplumedfile
```
The script above must be run with multiple replicas using the following command:
\verbatim
````
mpirun -np 4 plumed driver --mf_xtc alltraj.xtc --multi 4
\endverbatim
````
For more details on how the weights for configurations are determined using the wham method see the documentation
for the [WHAM](WHAM.md) action.
*/
//+ENDPLUMEDOC
Expand Down
18 changes: 9 additions & 9 deletions src/wham/WhamWeights.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -30,19 +30,16 @@ namespace wham {
Calculate and output weights for configurations using the weighted histogram analysis method.
This shortcut action allows you to calculate and output weights computed using the weighted histogram
analysis technique. For more detail on how this technique works see \ref REWEIGHT_WHAM
\par Examples
The following input can be used to analyze the output from a series of umbrella sampling calculations.
analysis technique. The following input demonstrates how this works in practise by showing you how this
action can be used to analyze the output from a series of umbrella sampling calculations.
The trajectory from each of the simulations run with the different biases should be concatenated into a
single trajectory before running the following analysis script on the concatenated trajectory using PLUMED
driver. The umbrella sampling simulations that will be analyzed using the script below applied a harmonic
restraint that restrained the torsional angle involving atoms 5, 7, 9 and 15 to particular values. The script
below calculates the reweighting weights for each of the trajectories and then applies the binless WHAM algorithm
to determine a weight for each configuration in the concatenated trajectory.
\plumedfile
```plumed
#SETTINGS NREPLICAS=4
phi: TORSION ATOMS=5,7,9,15
rp: RESTRAINT ARG=phi KAPPA=50.0 ...
Expand All @@ -55,13 +52,16 @@ rp: RESTRAINT ARG=phi KAPPA=50.0 ...
...
WHAM_WEIGHTS BIAS=rp.bias TEMP=300 FILE=wham-weights
\endplumedfile
```
The script above must be run with multiple replicas using the following command:
\verbatim
````
mpirun -np 4 plumed driver --mf_xtc alltraj.xtc --multi 4
\endverbatim
````
For more details on how the weights for configurations are determined using the wham method see the documentation
for the [WHAM](WHAM.md) action.
*/
//+ENDPLUMEDOC
Expand Down

1 comment on commit 9a50c3f

@PlumedBot
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Found broken examples in automatic/ANN.tmp
Found broken examples in automatic/CLASSICAL_MDS.tmp
Found broken examples in automatic/CONVERT_TO_FES.tmp
Found broken examples in automatic/COORDINATIONNUMBER.tmp
Found broken examples in automatic/DISTANCE_FROM_CONTOUR.tmp
Found broken examples in automatic/DUMPCUBE.tmp
Found broken examples in automatic/DUMPGRID.tmp
Found broken examples in automatic/EDS.tmp
Found broken examples in automatic/EMMI.tmp
Found broken examples in automatic/FIND_CONTOUR.tmp
Found broken examples in automatic/FIND_CONTOUR_SURFACE.tmp
Found broken examples in automatic/FIND_SPHERICAL_CONTOUR.tmp
Found broken examples in automatic/FOURIER_TRANSFORM.tmp
Found broken examples in automatic/FUNCPATHGENERAL.tmp
Found broken examples in automatic/FUNCPATHMSD.tmp
Found broken examples in automatic/FUNNEL.tmp
Found broken examples in automatic/FUNNEL_PS.tmp
Found broken examples in automatic/GPROPERTYMAP.tmp
Found broken examples in automatic/HISTOGRAM.tmp
Found broken examples in automatic/INTERPOLATE_GRID.tmp
Found broken examples in automatic/LOCAL_AVERAGE.tmp
Found broken examples in automatic/MAZE_OPTIMIZER_BIAS.tmp
Found broken examples in automatic/MAZE_RANDOM_ACCELERATION_MD.tmp
Found broken examples in automatic/MAZE_SIMULATED_ANNEALING.tmp
Found broken examples in automatic/MAZE_STEERED_MD.tmp
Found broken examples in automatic/METATENSOR.tmp
Found broken examples in automatic/MULTICOLVARDENS.tmp
Found broken examples in automatic/PCA.tmp
Found broken examples in automatic/PCAVARS.tmp
Found broken examples in automatic/PIV.tmp
Found broken examples in automatic/PYCVINTERFACE.tmp
Found broken examples in automatic/PYTHONFUNCTION.tmp
Found broken examples in automatic/Q3.tmp
Found broken examples in automatic/Q4.tmp
Found broken examples in automatic/Q6.tmp
Found broken examples in automatic/QUATERNION.tmp
Found broken examples in automatic/REWEIGHT_BIAS.tmp
Found broken examples in automatic/REWEIGHT_METAD.tmp
Found broken examples in automatic/SIZESHAPE_POSITION_LINEAR_PROJ.tmp
Found broken examples in automatic/SIZESHAPE_POSITION_MAHA_DIST.tmp
Found broken examples in AnalysisPP.md
Found broken examples in CollectiveVariablesPP.md
Found broken examples in MiscelaneousPP.md

Please sign in to comment.