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 \ref MAZE_LOSS function. The optimal biasing direction is determined by
      40             : performing simulated annealing.
      41             : 
      42             : \par 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 \ref MAZE_LOSS 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             : \plumedfile
      52             : UNITS LENGTH=A TIME=ps ENERGY=kcal/mol
      53             : 
      54             : MAZE_SIMULATED_ANNEALING ...
      55             :   LABEL=sa
      56             : 
      57             :   LOSS=l
      58             : 
      59             :   N_ITER=1000
      60             :   OPTIMIZER_STRIDE=200
      61             : 
      62             :   PROBABILITY_DECREASER=300
      63             :   COOLING=0.95
      64             :   COOLING_TYPE=geometric
      65             : 
      66             :   LIGAND=2635-2646
      67             :   PROTEIN=1-2634
      68             : ... MAZE_SIMULATED_ANNEALING
      69             : \endplumedfile
      70             : 
      71             : As shown above, each optimizer should be provided with the LIGAND and
      72             : the PROTEIN keywords.
      73             : 
      74             : */
      75             : //+ENDPLUMEDOC
      76             : 
      77             : /**
      78             :  * @class Simulated_Annealing Simulated_Annealing.cpp "maze/Simulated_Annealing.cpp"
      79             :  *
      80             :  * @brief Perform simulated annealing to compute an optimal bias direction.
      81             :  */
      82             : class Simulated_Annealing: public Optimizer {
      83             : public:
      84             :   /**
      85             :    * PLMD constructor.
      86             :    *
      87             :    * @param[in] ao PLMD::ActionOptions&.
      88             :    */
      89             :   explicit Simulated_Annealing(const ActionOptions& ao);
      90             : 
      91             :   /**
      92             :    * Destructor.
      93             :    */
      94             :   ~Simulated_Annealing();
      95             : 
      96             :   /**
      97             :    * Register PLMD keywords.
      98             :    *
      99             :    * @param[in] keys Keywords.
     100             :    */
     101             :   static void registerKeywords(Keywords& keys);
     102             : 
     103             :   /**
     104             :    * Each class deriving from Optimizer needs to override this function.
     105             :    */
     106             :   void optimize() override;
     107             : 
     108             :   /**
     109             :    * Reduce the temperature parameter.
     110             :    */
     111             :   void decrease_probability(unsigned int);
     112             : 
     113             : private:
     114             :   //! Temperature parameter.
     115             :   double probability_decreaser_;
     116             : 
     117             :   //! Cooling factor.
     118             :   double cooling_factor_;
     119             : 
     120             :   //! Cooling scheme.
     121             :   std::string cooling_scheme_;
     122             : };
     123             : 
     124             : // Register MAZE_SIMULATED_ANNEALING.
     125       10421 : PLUMED_REGISTER_ACTION(Simulated_Annealing, "MAZE_SIMULATED_ANNEALING")
     126             : 
     127           2 : void Simulated_Annealing::registerKeywords(Keywords& keys) {
     128           2 :   Optimizer::registerKeywords(keys);
     129             : 
     130           6 :   keys.add(
     131             :     "compulsory",
     132             :     "PROBABILITY_DECREASER",
     133             :     "Temperature-like parameter that is decreased during optimization to modify "
     134             :     "the Metropolis-Hastings acceptance probability."
     135             :   );
     136             : 
     137           4 :   keys.add(
     138             :     "compulsory",
     139             :     "COOLING",
     140             :     "Reduction factor for PROBABILITY_DECREASER, should be in (0, 1]."
     141             :   );
     142             : 
     143           4 :   keys.add(
     144             :     "compulsory",
     145             :     "COOLING_SCHEME",
     146             :     "Cooling scheme: geometric."
     147             :   );
     148           2 : }
     149             : 
     150           1 : Simulated_Annealing::Simulated_Annealing(const ActionOptions& ao)
     151           1 :   : PLUMED_OPT_INIT(ao)
     152             : {
     153           1 :   log.printf("maze> Simulated annealing optimizer.\n");
     154             : 
     155           2 :   if(keywords.exists("COOLING")) {
     156           1 :     parse("COOLING", cooling_factor_);
     157             : 
     158           1 :     plumed_massert(
     159             :       cooling_factor_ > 0 && cooling_factor_ <= 1,
     160             :       "maze> COOLING should be in (0, 1]; preferably 0.95.\n"
     161             :     );
     162             :   }
     163             : 
     164           2 :   if(keywords.exists("PROBABILITY_DECREASER")) {
     165           1 :     parse("PROBABILITY_DECREASER", probability_decreaser_);
     166             : 
     167           1 :     plumed_massert(
     168             :       probability_decreaser_ > 0,
     169             :       "maze> PROBABILITY_DECREASER should be explicitly specified and positive.\n");
     170             :   }
     171             : 
     172           2 :   if(keywords.exists("COOLING_SCHEME")) {
     173           1 :     parse("COOLING_SCHEME", cooling_scheme_);
     174             : 
     175           1 :     log.printf(
     176             :       "maze> COOLING_SCHEME read: %s.\n",
     177             :       cooling_scheme_.c_str()
     178             :     );
     179             :   }
     180             : 
     181           2 :   set_label("SIMULATED_ANNEALING");
     182             : 
     183             :   // Calculate an optimal direction at the beginning of the MD simulation.
     184             :   start_step_0();
     185             : 
     186           1 :   checkRead();
     187           1 : }
     188             : 
     189           2 : Simulated_Annealing::~Simulated_Annealing() {
     190           1 :   delete neighbor_list_;
     191           2 : }
     192             : 
     193          30 : void Simulated_Annealing::decrease_probability(unsigned int time) {
     194          30 :   if (cooling_scheme_ == "linear") {
     195           0 :     probability_decreaser_ -= time * cooling_factor_;
     196             :   }
     197          30 :   else if (cooling_scheme_ == "exponential") {
     198           0 :     probability_decreaser_ *= pow(cooling_factor_, time);
     199             :   }
     200          30 :   else if (cooling_scheme_ == "geometric") {
     201          30 :     probability_decreaser_ *= cooling_factor_;
     202             :   }
     203           0 :   else if (cooling_scheme_ == "logarithmic") {
     204           0 :     probability_decreaser_ = cooling_factor_ / std::log(time + 1);
     205             :   }
     206           0 :   else if (cooling_scheme_ == "hoffman") {
     207           0 :     probability_decreaser_ = (cooling_factor_ - 1) / std::log(time);
     208             :   }
     209          30 : }
     210             : 
     211           3 : void Simulated_Annealing::optimize() {
     212           3 :   sampling_r_ = sampling_radius();
     213             :   double rad_s;
     214           3 :   const unsigned nl_size = neighbor_list_->size();
     215             : 
     216           3 :   Vector distance, distance_next;
     217             : 
     218          33 :   for (unsigned int iter=0; iter < get_n_iterations(); ++iter) {
     219             :     double action = 0;
     220             :     double action_next = 0;
     221             : 
     222          30 :     rad_s = rnd::next_double(sampling_r_);
     223          30 :     Vector dev = rnd::next_plmd_vector(rad_s);
     224             : 
     225          30 :     #pragma omp parallel num_threads(get_n_threads_openmp())
     226             :     {
     227             :       #pragma omp for reduction(+:action_next, action)
     228             :       for (unsigned int i=0; i < nl_size; i++) {
     229             :         unsigned i0 = neighbor_list_->getClosePair(i).first;
     230             :         unsigned i1 = neighbor_list_->getClosePair(i).second;
     231             : 
     232             :         if (getAbsoluteIndex(i0) == getAbsoluteIndex(i1)) {
     233             :           continue;
     234             :         }
     235             : 
     236             :         if (pbc_) {
     237             :           distance = pbcDistance(
     238             :                        getPosition(i0) + get_opt(),
     239             :                        getPosition(i1)
     240             :                      );
     241             : 
     242             :           distance_next = pbcDistance(
     243             :                             getPosition(i0) + dev,
     244             :                             getPosition(i1)
     245             :                           );
     246             :         }
     247             :         else {
     248             :           distance = delta(
     249             :                        getPosition(i0) + get_opt(),
     250             :                        getPosition(i1)
     251             :                      );
     252             : 
     253             :           distance_next = delta(
     254             :                             getPosition(i0) + dev,
     255             :                             getPosition(i1)
     256             :                           );
     257             :         }
     258             : 
     259             :         action += pairing(distance.modulo());
     260             :         action_next += pairing(distance_next.modulo());
     261             :       }
     262             :     }
     263             : 
     264             :     double p = std::min(
     265          60 :                  1.0,
     266          30 :                  std::exp(-(action_next - action) / probability_decreaser_)
     267          30 :                );
     268             : 
     269          30 :     double r = rnd::next_double();
     270             : 
     271          30 :     if (r < p) {
     272             :       set_opt(dev);
     273             :       set_opt_value(action_next);
     274             :     }
     275             : 
     276          30 :     decrease_probability(iter);
     277             :   }
     278             : 
     279           3 :   Vector s = get_opt() / modulo(get_opt());
     280             :   set_opt(s);
     281           3 : }
     282             : 
     283             : } // namespace maze
     284             : } // namespace PLMD