Line data    Source code

       1             : /* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
       2             : Copyright (c) 2019 Jakub Rydzewski (jr@fizyka.umk.pl). All rights reserved.
       3             : 
       4             : See http://www.maze-code.github.io for more information.
       5             : 
       6             : This file is part of maze.
       7             : 
       8             : maze is free software: you can redistribute it and/or modify it under the
       9             : terms of the GNU Lesser General Public License as published by the Free
      10             : Software Foundation, either version 3 of the License, or (at your option)
      11             : any later version.
      12             : 
      13             : maze is distributed in the hope that it will be useful, but WITHOUT ANY
      14             : WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
      15             : FOR A PARTICULAR PURPOSE.
      16             : 
      17             : See the 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 maze. If not, see <https://www.gnu.org/licenses/>.
      21             : +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ */
      22             : 
      23             : /**
      24             :  * @file Simulated_Annealing.cpp
      25             :  *
      26             :  * @author J. Rydzewski (jr@fizyka.umk.pl)
      27             :  */
      28             : 
      29             : #include "core/ActionRegister.h"
      30             : #include "Optimizer.h"
      31             : 
      32             : namespace PLMD {
      33             : namespace maze {
      34             : 
      35             : //+PLUMEDOC MAZE_OPTIMIZER MAZE_SIMULATED_ANNEALING
      36             : /*
      37             : 
      38             : Calculates the biasing direction along which the ligand unbinds by minimizing
      39             : the [MAZE_LOSS](MAZE_LOSS.md) function. The optimal biasing direction is determined by
      40             : performing simulated annealing.
      41             : 
      42             : ## Examples
      43             : 
      44             : Every optimizer implemented in the maze module needs a loss function as an
      45             : argument, and it should be passed using the [MAZE_LOSS](MAZE_LOSS.md) keyword.
      46             : 
      47             : In the following example simulated annealing is launched for 1000 iterations
      48             : as the optimizer for the loss function every 200 ps. The geometric cooling
      49             : scheme is used.
      50             : 
      51             : ```plumed
      52             : UNITS LENGTH=A TIME=ps ENERGY=kcal/mol
      53             : 
      54             : l: MAZE_LOSS PARAMS=1,1,1
      55             : 
      56             : MAZE_SIMULATED_ANNEALING ...
      57             :   LABEL=sa
      58             : 
      59             :   LOSS=l
      60             : 
      61             :   N_ITER=1000
      62             :   OPTIMIZER_STRIDE=200
      63             : 
      64             :   PROBABILITY_DECREASER=300
      65             :   COOLING=0.95
      66             :   COOLING_SCHEME=geometric
      67             : 
      68             :   LIGAND=2635-2646
      69             :   PROTEIN=1-2634
      70             : ... MAZE_SIMULATED_ANNEALING
      71             : ```
      72             : 
      73             : As shown above, each optimizer should be provided with the LIGAND and
      74             : the PROTEIN keywords.
      75             : 
      76             : */
      77             : //+ENDPLUMEDOC
      78             : 
      79             : /**
      80             :  * @class Simulated_Annealing Simulated_Annealing.cpp "maze/Simulated_Annealing.cpp"
      81             :  *
      82             :  * @brief Perform simulated annealing to compute an optimal bias direction.
      83             :  */
      84             : class Simulated_Annealing: public Optimizer {
      85             : public:
      86             :   /**
      87             :    * PLMD constructor.
      88             :    *
      89             :    * @param[in] ao PLMD::ActionOptions&.
      90             :    */
      91             :   explicit Simulated_Annealing(const ActionOptions& ao);
      92             : 
      93             :   /**
      94             :    * Register PLMD keywords.
      95             :    *
      96             :    * @param[in] keys Keywords.
      97             :    */
      98             :   static void registerKeywords(Keywords& keys);
      99             : 
     100             :   /**
     101             :    * Each class deriving from Optimizer needs to override this function.
     102             :    */
     103             :   void optimize() override;
     104             : 
     105             :   /**
     106             :    * Reduce the temperature parameter.
     107             :    */
     108             :   void decrease_probability(unsigned int);
     109             : 
     110             : private:
     111             :   //! Temperature parameter.
     112             :   double probability_decreaser_;
     113             : 
     114             :   //! Cooling factor.
     115             :   double cooling_factor_;
     116             : 
     117             :   //! Cooling scheme.
     118             :   std::string cooling_scheme_;
     119             : };
     120             : 
     121             : // Register MAZE_SIMULATED_ANNEALING.
     122             : PLUMED_REGISTER_ACTION(Simulated_Annealing, "MAZE_SIMULATED_ANNEALING")
     123             : 
     124           5 : void Simulated_Annealing::registerKeywords(Keywords& keys) {
     125           5 :   Optimizer::registerKeywords(keys);
     126             : 
     127           5 :   keys.add(
     128             :     "compulsory",
     129             :     "PROBABILITY_DECREASER",
     130             :     "Temperature-like parameter that is decreased during optimization to modify "
     131             :     "the Metropolis-Hastings acceptance probability."
     132             :   );
     133             : 
     134           5 :   keys.add(
     135             :     "compulsory",
     136             :     "COOLING",
     137             :     "Reduction factor for PROBABILITY_DECREASER, should be in (0, 1]."
     138             :   );
     139             : 
     140           5 :   keys.add(
     141             :     "compulsory",
     142             :     "COOLING_SCHEME",
     143             :     "Cooling scheme: geometric."
     144             :   );
     145           5 : }
     146             : 
     147           1 : Simulated_Annealing::Simulated_Annealing(const ActionOptions& ao)
     148           1 :   : PLUMED_OPT_INIT(ao) {
     149           1 :   log.printf("maze> Simulated annealing optimizer.\n");
     150             : 
     151           1 :   if(keywords.exists("COOLING")) {
     152           1 :     parse("COOLING", cooling_factor_);
     153             : 
     154           1 :     plumed_massert(
     155             :       cooling_factor_ > 0 && cooling_factor_ <= 1,
     156             :       "maze> COOLING should be in (0, 1]; preferably 0.95.\n"
     157             :     );
     158             :   }
     159             : 
     160           1 :   if(keywords.exists("PROBABILITY_DECREASER")) {
     161           1 :     parse("PROBABILITY_DECREASER", probability_decreaser_);
     162             : 
     163           1 :     plumed_massert(
     164             :       probability_decreaser_ > 0,
     165             :       "maze> PROBABILITY_DECREASER should be explicitly specified and positive.\n");
     166             :   }
     167             : 
     168           1 :   if(keywords.exists("COOLING_SCHEME")) {
     169           1 :     parse("COOLING_SCHEME", cooling_scheme_);
     170             : 
     171           1 :     log.printf(
     172             :       "maze> COOLING_SCHEME read: %s.\n",
     173             :       cooling_scheme_.c_str()
     174             :     );
     175             :   }
     176             : 
     177           2 :   set_label("SIMULATED_ANNEALING");
     178             : 
     179             :   // Calculate an optimal direction at the beginning of the MD simulation.
     180             :   start_step_0();
     181             : 
     182           1 :   checkRead();
     183           1 : }
     184             : 
     185          30 : void Simulated_Annealing::decrease_probability(unsigned int time) {
     186          30 :   if (cooling_scheme_ == "linear") {
     187           0 :     probability_decreaser_ -= time * cooling_factor_;
     188          30 :   } else if (cooling_scheme_ == "exponential") {
     189           0 :     probability_decreaser_ *= pow(cooling_factor_, time);
     190          30 :   } else if (cooling_scheme_ == "geometric") {
     191          30 :     probability_decreaser_ *= cooling_factor_;
     192           0 :   } else if (cooling_scheme_ == "logarithmic") {
     193           0 :     probability_decreaser_ = cooling_factor_ / std::log(time + 1);
     194           0 :   } else if (cooling_scheme_ == "hoffman") {
     195           0 :     probability_decreaser_ = (cooling_factor_ - 1) / std::log(time);
     196             :   }
     197          30 : }
     198             : 
     199           3 : void Simulated_Annealing::optimize() {
     200           3 :   sampling_r_ = sampling_radius();
     201             :   double rad_s;
     202           3 :   const unsigned nl_size = neighbor_list_->size();
     203             : 
     204           3 :   Vector distance, distance_next;
     205             : 
     206          33 :   for (unsigned int iter=0; iter < get_n_iterations(); ++iter) {
     207             :     double action = 0;
     208             :     double action_next = 0;
     209             : 
     210          30 :     rad_s = rnd::next_double(sampling_r_);
     211          30 :     Vector dev = rnd::next_plmd_vector(rad_s);
     212             : 
     213          30 :     #pragma omp parallel num_threads(get_n_threads_openmp())
     214             :     {
     215             :       #pragma omp for reduction(+:action_next, action)
     216             :       for (unsigned int i=0; i < nl_size; i++) {
     217             :         unsigned i0 = neighbor_list_->getClosePair(i).first;
     218             :         unsigned i1 = neighbor_list_->getClosePair(i).second;
     219             : 
     220             :         if (getAbsoluteIndex(i0) == getAbsoluteIndex(i1)) {
     221             :           continue;
     222             :         }
     223             : 
     224             :         if (pbc_) {
     225             :           distance = pbcDistance(
     226             :                        getPosition(i0) + get_opt(),
     227             :                        getPosition(i1)
     228             :                      );
     229             : 
     230             :           distance_next = pbcDistance(
     231             :                             getPosition(i0) + dev,
     232             :                             getPosition(i1)
     233             :                           );
     234             :         } else {
     235             :           distance = delta(
     236             :                        getPosition(i0) + get_opt(),
     237             :                        getPosition(i1)
     238             :                      );
     239             : 
     240             :           distance_next = delta(
     241             :                             getPosition(i0) + dev,
     242             :                             getPosition(i1)
     243             :                           );
     244             :         }
     245             : 
     246             :         action += pairing(distance.modulo());
     247             :         action_next += pairing(distance_next.modulo());
     248             :       }
     249             :     }
     250             : 
     251             :     double p = std::min(
     252          60 :                  1.0,
     253          30 :                  std::exp(-(action_next - action) / probability_decreaser_)
     254          30 :                );
     255             : 
     256          30 :     double r = rnd::next_double();
     257             : 
     258          30 :     if (r < p) {
     259             :       set_opt(dev);
     260             :       set_opt_value(action_next);
     261             :     }
     262             : 
     263          30 :     decrease_probability(iter);
     264             :   }
     265             : 
     266           3 :   Vector s = get_opt() / modulo(get_opt());
     267             :   set_opt(s);
     268           3 : }
     269             : 
     270             : } // namespace maze
     271             : } // namespace PLMD