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             :    * Register PLMD keywords.
      93             :    *
      94             :    * @param[in] keys Keywords.
      95             :    */
      96             :   static void registerKeywords(Keywords& keys);
      97             : 
      98             :   /**
      99             :    * Each class deriving from Optimizer needs to override this function.
     100             :    */
     101             :   void optimize() override;
     102             : 
     103             :   /**
     104             :    * Reduce the temperature parameter.
     105             :    */
     106             :   void decrease_probability(unsigned int);
     107             : 
     108             : private:
     109             :   //! Temperature parameter.
     110             :   double probability_decreaser_;
     111             : 
     112             :   //! Cooling factor.
     113             :   double cooling_factor_;
     114             : 
     115             :   //! Cooling scheme.
     116             :   std::string cooling_scheme_;
     117             : };
     118             : 
     119             : // Register MAZE_SIMULATED_ANNEALING.
     120             : PLUMED_REGISTER_ACTION(Simulated_Annealing, "MAZE_SIMULATED_ANNEALING")
     121             : 
     122           3 : void Simulated_Annealing::registerKeywords(Keywords& keys) {
     123           3 :   Optimizer::registerKeywords(keys);
     124             : 
     125           6 :   keys.add(
     126             :     "compulsory",
     127             :     "PROBABILITY_DECREASER",
     128             :     "Temperature-like parameter that is decreased during optimization to modify "
     129             :     "the Metropolis-Hastings acceptance probability."
     130             :   );
     131             : 
     132           6 :   keys.add(
     133             :     "compulsory",
     134             :     "COOLING",
     135             :     "Reduction factor for PROBABILITY_DECREASER, should be in (0, 1]."
     136             :   );
     137             : 
     138           6 :   keys.add(
     139             :     "compulsory",
     140             :     "COOLING_SCHEME",
     141             :     "Cooling scheme: geometric."
     142             :   );
     143           3 : }
     144             : 
     145           1 : Simulated_Annealing::Simulated_Annealing(const ActionOptions& ao)
     146           1 :   : PLUMED_OPT_INIT(ao)
     147             : {
     148           1 :   log.printf("maze> Simulated annealing optimizer.\n");
     149             : 
     150           2 :   if(keywords.exists("COOLING")) {
     151           1 :     parse("COOLING", cooling_factor_);
     152             : 
     153           1 :     plumed_massert(
     154             :       cooling_factor_ > 0 && cooling_factor_ <= 1,
     155             :       "maze> COOLING should be in (0, 1]; preferably 0.95.\n"
     156             :     );
     157             :   }
     158             : 
     159           2 :   if(keywords.exists("PROBABILITY_DECREASER")) {
     160           1 :     parse("PROBABILITY_DECREASER", probability_decreaser_);
     161             : 
     162           1 :     plumed_massert(
     163             :       probability_decreaser_ > 0,
     164             :       "maze> PROBABILITY_DECREASER should be explicitly specified and positive.\n");
     165             :   }
     166             : 
     167           2 :   if(keywords.exists("COOLING_SCHEME")) {
     168           1 :     parse("COOLING_SCHEME", cooling_scheme_);
     169             : 
     170           1 :     log.printf(
     171             :       "maze> COOLING_SCHEME read: %s.\n",
     172             :       cooling_scheme_.c_str()
     173             :     );
     174             :   }
     175             : 
     176           2 :   set_label("SIMULATED_ANNEALING");
     177             : 
     178             :   // Calculate an optimal direction at the beginning of the MD simulation.
     179             :   start_step_0();
     180             : 
     181           1 :   checkRead();
     182           1 : }
     183             : 
     184          30 : void Simulated_Annealing::decrease_probability(unsigned int time) {
     185          30 :   if (cooling_scheme_ == "linear") {
     186           0 :     probability_decreaser_ -= time * cooling_factor_;
     187             :   }
     188          30 :   else if (cooling_scheme_ == "exponential") {
     189           0 :     probability_decreaser_ *= pow(cooling_factor_, time);
     190             :   }
     191          30 :   else if (cooling_scheme_ == "geometric") {
     192          30 :     probability_decreaser_ *= cooling_factor_;
     193             :   }
     194           0 :   else if (cooling_scheme_ == "logarithmic") {
     195           0 :     probability_decreaser_ = cooling_factor_ / std::log(time + 1);
     196             :   }
     197           0 :   else if (cooling_scheme_ == "hoffman") {
     198           0 :     probability_decreaser_ = (cooling_factor_ - 1) / std::log(time);
     199             :   }
     200          30 : }
     201             : 
     202           3 : void Simulated_Annealing::optimize() {
     203           3 :   sampling_r_ = sampling_radius();
     204             :   double rad_s;
     205           3 :   const unsigned nl_size = neighbor_list_->size();
     206             : 
     207           3 :   Vector distance, distance_next;
     208             : 
     209          33 :   for (unsigned int iter=0; iter < get_n_iterations(); ++iter) {
     210             :     double action = 0;
     211             :     double action_next = 0;
     212             : 
     213          30 :     rad_s = rnd::next_double(sampling_r_);
     214          30 :     Vector dev = rnd::next_plmd_vector(rad_s);
     215             : 
     216          30 :     #pragma omp parallel num_threads(get_n_threads_openmp())
     217             :     {
     218             :       #pragma omp for reduction(+:action_next, action)
     219             :       for (unsigned int i=0; i < nl_size; i++) {
     220             :         unsigned i0 = neighbor_list_->getClosePair(i).first;
     221             :         unsigned i1 = neighbor_list_->getClosePair(i).second;
     222             : 
     223             :         if (getAbsoluteIndex(i0) == getAbsoluteIndex(i1)) {
     224             :           continue;
     225             :         }
     226             : 
     227             :         if (pbc_) {
     228             :           distance = pbcDistance(
     229             :                        getPosition(i0) + get_opt(),
     230             :                        getPosition(i1)
     231             :                      );
     232             : 
     233             :           distance_next = pbcDistance(
     234             :                             getPosition(i0) + dev,
     235             :                             getPosition(i1)
     236             :                           );
     237             :         }
     238             :         else {
     239             :           distance = delta(
     240             :                        getPosition(i0) + get_opt(),
     241             :                        getPosition(i1)
     242             :                      );
     243             : 
     244             :           distance_next = delta(
     245             :                             getPosition(i0) + dev,
     246             :                             getPosition(i1)
     247             :                           );
     248             :         }
     249             : 
     250             :         action += pairing(distance.modulo());
     251             :         action_next += pairing(distance_next.modulo());
     252             :       }
     253             :     }
     254             : 
     255             :     double p = std::min(
     256          60 :                  1.0,
     257          30 :                  std::exp(-(action_next - action) / probability_decreaser_)
     258          30 :                );
     259             : 
     260          30 :     double r = rnd::next_double();
     261             : 
     262          30 :     if (r < p) {
     263             :       set_opt(dev);
     264             :       set_opt_value(action_next);
     265             :     }
     266             : 
     267          30 :     decrease_probability(iter);
     268             :   }
     269             : 
     270           3 :   Vector s = get_opt() / modulo(get_opt());
     271             :   set_opt(s);
     272           3 : }
     273             : 
     274             : } // namespace maze
     275             : } // namespace PLMD