Line data    Source code

       1             : /* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
       2             :    Copyright (c) 2017-2020 The plumed team
       3             :    (see the PEOPLE file at the root of the distribution for a list of names)
       4             : 
       5             :    See http://www.plumed.org for more information.
       6             : 
       7             :    This file is part of plumed, version 2.
       8             : 
       9             :    plumed is free software: you can redistribute it and/or modify
      10             :    it under the terms of the GNU Lesser General Public License as published by
      11             :    the Free Software Foundation, either version 3 of the License, or
      12             :    (at your option) any later version.
      13             : 
      14             :    plumed is distributed in the hope that it will be useful,
      15             :    but WITHOUT ANY WARRANTY; without even the implied warranty of
      16             :    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
      17             :    GNU Lesser General Public License for more details.
      18             : 
      19             :    You should have received a copy of the GNU Lesser General Public License
      20             :    along with plumed.  If not, see <http://www.gnu.org/licenses/>.
      21             : +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ */
      22             : #include "core/ActionRegister.h"
      23             : #include "core/ActionWithValue.h"
      24             : #include "core/ActionWithArguments.h"
      25             : #include "core/PlumedMain.h"
      26             : #include "tools/Communicator.h"
      27             : 
      28             : //+PLUMEDOC REWEIGHTING WHAM
      29             : /*
      30             : Calculate the weights for configurations using the weighted histogram analysis method.
      31             : 
      32             : The example input below shows how this command is used.
      33             : 
      34             : ```plumed
      35             : #SETTINGS NREPLICAS=4
      36             : phi: TORSION ATOMS=5,7,9,15
      37             : psi: TORSION ATOMS=7,9,15,17
      38             : rp: RESTRAINT ARG=phi KAPPA=50.0 AT=@replicas:{-3.00,-1.45,0.10,1.65}
      39             : # Get the bias on each of the replicas
      40             : rep: GATHER_REPLICAS ARG=rp.bias
      41             : # Merge the biases on each of the replicas into a single vector
      42             : all: CONCATENATE ARG=rep.*
      43             : # Collect all the bias values
      44             : col: COLLECT TYPE=vector ARG=all STRIDE=1
      45             : wham: WHAM ARG=col TEMP=300
      46             : DUMPVECTOR ARG=wham FILE=wham_data
      47             : ```
      48             : 
      49             : 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
      50             : 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:
      51             : 
      52             : $$
      53             : P( \underline{T} ) \propto \prod_{j=1}^M \prod_{k=1}^N (c_k w_{kj} p_j)^{t_{kj}}
      54             : $$
      55             : 
      56             : 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
      57             : $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
      58             : having a set of CV values that lie within the range for the $j$th bin.
      59             : 
      60             : 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.
      61             : 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:
      62             : 
      63             : $$
      64             : w_{kj} = e^{+\beta V_k(x_j)}
      65             : $$
      66             : 
      67             : 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
      68             : 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
      69             : 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.
      70             : 
      71             : 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.
      72             : This constrained optimization must be performed using a set of Lagrange multipliers, $\lambda_k$, that ensure that each of the biased probability distributions
      73             : 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
      74             : we actually chose to minimize
      75             : 
      76             : $$
      77             : 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)
      78             : $$
      79             : 
      80             : After some manipulations the following (WHAM) equations emerge:
      81             : 
      82             : $$
      83             : \begin{aligned}
      84             : p_j & \propto \frac{\sum_{k=1}^N t_{kj}}{\sum_k c_k w_{kj}} \\
      85             : c_k & =\frac{1}{\sum_{j=1}^M w_{kj} p_j}
      86             : \end{aligned}
      87             : $$
      88             : 
      89             : 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
      90             : 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.
      91             : 
      92             : 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.
      93             : 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.
      94             : This observation is important as it is the basis of the binless formulation of WHAM that is implemented within PLUMED.
      95             : 
      96             : */
      97             : //+ENDPLUMEDOC
      98             : 
      99             : 
     100             : namespace PLMD {
     101             : namespace wham {
     102             : 
     103             : class Wham :
     104             :   public ActionWithValue,
     105             :   public ActionWithArguments {
     106             : private:
     107             :   double thresh, simtemp;
     108             :   unsigned nreplicas;
     109             :   unsigned maxiter;
     110             : public:
     111             :   static void registerKeywords(Keywords&);
     112             :   explicit Wham(const ActionOptions&ao);
     113           0 :   unsigned getNumberOfDerivatives() {
     114           0 :     return 0;
     115             :   }
     116             :   void calculate() override ;
     117           0 :   void apply() override {}
     118             : };
     119             : 
     120             : PLUMED_REGISTER_ACTION(Wham,"WHAM")
     121             : 
     122          29 : void Wham::registerKeywords(Keywords& keys ) {
     123          29 :   Action::registerKeywords( keys );
     124          29 :   ActionWithValue::registerKeywords( keys );
     125          29 :   ActionWithArguments::registerKeywords( keys );
     126          58 :   keys.addInputKeyword("compulsory","ARG","scalar/vector/matrix","the stored values for the bias");
     127          29 :   keys.add("compulsory","MAXITER","1000","maximum number of iterations for WHAM algorithm");
     128          29 :   keys.add("compulsory","WHAMTOL","1e-10","threshold for convergence of WHAM algorithm");
     129          29 :   keys.add("optional","TEMP","the system temperature.  This is not required if your MD code passes this quantity to PLUMED");
     130          29 :   keys.remove("NUMERICAL_DERIVATIVES");
     131          58 :   keys.setValueDescription("vector","the vector of WHAM weights to use for reweighting the elements in a time series");
     132          29 : }
     133             : 
     134          13 : Wham::Wham(const ActionOptions&ao):
     135             :   Action(ao),
     136             :   ActionWithValue(ao),
     137          13 :   ActionWithArguments(ao) {
     138             :   // Read in the temperature
     139          13 :   simtemp=getkBT();
     140          13 :   if(simtemp==0) {
     141           0 :     error("The MD engine does not pass the temperature to plumed so you have to specify it using TEMP");
     142             :   }
     143             :   // Now read in parameters of WHAM
     144          13 :   parse("MAXITER",maxiter);
     145          13 :   parse("WHAMTOL",thresh);
     146          13 :   if(comm.Get_rank()==0) {
     147          13 :     nreplicas=multi_sim_comm.Get_size();
     148             :   }
     149          13 :   comm.Bcast(nreplicas,0);
     150          13 :   addValue( getPntrToArgument(0)->getShape() );
     151          13 :   setNotPeriodic();
     152          13 : }
     153             : 
     154          12 : void Wham::calculate() {
     155             :   // Retrieve the values that were stored for the biase
     156          12 :   std::vector<double> stored_biases( getPntrToArgument(0)->getNumberOfValues() );
     157       25284 :   for(unsigned i=0; i<stored_biases.size(); ++i) {
     158       25272 :     stored_biases[i] = getPntrToArgument(0)->get(i);
     159             :   }
     160             :   // Get the minimum value of the bias
     161          12 :   double minv = *min_element(std::begin(stored_biases), std::end(stored_biases));
     162             :   // Resize final weights array
     163          12 :   plumed_assert( stored_biases.size()%nreplicas==0 );
     164          12 :   std::vector<double> final_weights( stored_biases.size() / nreplicas, 1.0 );
     165          12 :   if( getPntrToComponent(0)->getNumberOfValues()!=final_weights.size() ) {
     166          12 :     std::vector<unsigned> shape(1);
     167          12 :     shape[0]=final_weights.size();
     168          12 :     getPntrToComponent(0)->setShape( shape );
     169             :   }
     170             :   // Offset and exponential of the bias
     171          12 :   std::vector<double> expv( stored_biases.size() );
     172       25284 :   for(unsigned i=0; i<expv.size(); ++i) {
     173       25272 :     expv[i] = exp( (-stored_biases[i]+minv) / simtemp );
     174             :   }
     175             :   // Initialize Z
     176          12 :   std::vector<double> Z( nreplicas, 1.0 ), oldZ( nreplicas );
     177             :   // Now the iterative loop to calculate the WHAM weights
     178        4584 :   for(unsigned iter=0; iter<maxiter; ++iter) {
     179             :     // Store Z
     180       32088 :     for(unsigned j=0; j<Z.size(); ++j) {
     181       27504 :       oldZ[j]=Z[j];
     182             :     }
     183             :     // Recompute weights
     184             :     double norm=0;
     185     1613568 :     for(unsigned j=0; j<final_weights.size(); ++j) {
     186             :       double ew=0;
     187    11262888 :       for(unsigned k=0; k<Z.size(); ++k) {
     188     9653904 :         ew += expv[j*Z.size()+k]  / Z[k];
     189             :       }
     190     1608984 :       final_weights[j] = 1.0 / ew;
     191     1608984 :       norm += final_weights[j];
     192             :     }
     193             :     // Normalize weights
     194     1613568 :     for(unsigned j=0; j<final_weights.size(); ++j) {
     195     1608984 :       final_weights[j] /= norm;
     196             :     }
     197             :     // Recompute Z
     198       32088 :     for(unsigned j=0; j<Z.size(); ++j) {
     199       27504 :       Z[j] = 0.0;
     200             :     }
     201     1613568 :     for(unsigned j=0; j<final_weights.size(); ++j) {
     202    11262888 :       for(unsigned k=0; k<Z.size(); ++k) {
     203     9653904 :         Z[k] += final_weights[j]*expv[j*Z.size()+k];
     204             :       }
     205             :     }
     206             :     // Normalize Z and compute change in Z
     207             :     double change=0;
     208             :     norm=0;
     209       32088 :     for(unsigned k=0; k<Z.size(); ++k) {
     210       27504 :       norm+=Z[k];
     211             :     }
     212       32088 :     for(unsigned k=0; k<Z.size(); ++k) {
     213       27504 :       Z[k] /= norm;
     214       27504 :       double d = std::log( Z[k] / oldZ[k] );
     215       27504 :       change += d*d;
     216             :     }
     217        4584 :     if( change<thresh ) {
     218        4224 :       for(unsigned j=0; j<final_weights.size(); ++j) {
     219        4212 :         getPntrToComponent(0)->set( j, final_weights[j] );
     220             :       }
     221          12 :       return;
     222             :     }
     223             :   }
     224           0 :   error("Too many iterations in WHAM" );
     225             : }
     226             : 
     227             : }
     228             : }