diff --git a/src/generic/DumpProjections.cpp b/src/generic/DumpProjections.cpp index 29e308cb6e..ee89845a18 100644 --- a/src/generic/DumpProjections.cpp +++ b/src/generic/DumpProjections.cpp @@ -116,7 +116,7 @@ DumpProjections::DumpProjections(const ActionOptions&ao): checkRead(); for(unsigned i=0; igetRank()!=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(); } } diff --git a/src/wham/Wham.cpp b/src/wham/Wham.cpp index 62396f450b..cfb37c29e8 100644 --- a/src/wham/Wham.cpp +++ b/src/wham/Wham.cpp @@ -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 diff --git a/src/wham/WhamHistogram.cpp b/src/wham/WhamHistogram.cpp index ad1e0ef6cb..d7aebc894b 100644 --- a/src/wham/WhamHistogram.cpp +++ b/src/wham/WhamHistogram.cpp @@ -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 @@ -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 @@ -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 diff --git a/src/wham/WhamWeights.cpp b/src/wham/WhamWeights.cpp index 716524fdcc..cff85f242e 100644 --- a/src/wham/WhamWeights.cpp +++ b/src/wham/WhamWeights.cpp @@ -30,11 +30,8 @@ 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 @@ -42,7 +39,7 @@ restraint that restrained the torsional angle involving atoms 5, 7, 9 and 15 to 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 ... @@ -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