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 Optimizer.cpp
      25             :  *
      26             :  * @author J. Rydzewski (jr@fizyka.umk.pl)
      27             :  */
      28             : 
      29             : #include "Optimizer.h"
      30             : #include "core/PlumedMain.h"
      31             : #include "tools/Tools.h"
      32             : 
      33             : namespace PLMD {
      34             : namespace maze {
      35             : 
      36          17 : void Optimizer::registerKeywords(Keywords& keys) {
      37          17 :   Colvar::registerKeywords(keys);
      38             : 
      39          34 :   keys.addFlag(
      40             :     "SERIAL",
      41             :     false,
      42             :     "Perform the simulation in serial -- used only for debugging purposes, "
      43             :     "should not be used otherwise."
      44             :   );
      45             : 
      46          34 :   keys.addFlag(
      47             :     "PAIR",
      48             :     false,
      49             :     "Pair only the 1st element of the 1st group with the 1st element in the "
      50             :     "second, etc."
      51             :   );
      52             : 
      53          34 :   keys.addFlag(
      54             :     "NLIST",
      55             :     false,
      56             :     "Use a neighbor list of ligand-protein atom pairs to speed up the "
      57             :     "calculating of the distances."
      58             :   );
      59             : 
      60          34 :   keys.add(
      61             :     "optional",
      62             :     "NL_CUTOFF",
      63             :     "Neighbor list cut-off for the distances of ligand-protein atom pairs."
      64             :   );
      65             : 
      66          34 :   keys.add(
      67             :     "optional",
      68             :     "NL_STRIDE",
      69             :     "Update stride for the ligand-protein atom pairs in the neighbor list."
      70             :   );
      71             : 
      72          34 :   keys.add(
      73             :     "compulsory",
      74             :     "N_ITER",
      75             :     "Number of optimization steps. Required only for optimizers, do not pass "
      76             :     "this keyword to the fake optimizers (results in crash) , e.g., random "
      77             :     "walk, steered MD, or random acceleration MD."
      78             :   );
      79             : 
      80          34 :   keys.add(
      81             :     "optional",
      82             :     "LOSS",
      83             :     "Loss function describing ligand-protein interactions required by every "
      84             :     "optimizer."
      85             :   );
      86             : 
      87          34 :   keys.add(
      88             :     "atoms",
      89             :     "LIGAND",
      90             :     "Indices of ligand atoms."
      91             :   );
      92             : 
      93          34 :   keys.add(
      94             :     "atoms",
      95             :     "PROTEIN",
      96             :     "Indices of protein atoms."
      97             :   );
      98             : 
      99          34 :   keys.add(
     100             :     "compulsory",
     101             :     "OPTIMIZER_STRIDE",
     102             :     "Optimizer stride. Sets up a callback function that launches the "
     103             :     "optimization process every OPTIMIZER_STRIDE."
     104             :   );
     105             : 
     106          34 :   keys.addOutputComponent(
     107             :     "x",
     108             :     "default",
     109             :     "Optimal biasing direction; x component."
     110             :   );
     111             : 
     112          34 :   keys.addOutputComponent(
     113             :     "y",
     114             :     "default",
     115             :     "Optimal biasing direction; y component."
     116             :   );
     117             : 
     118          34 :   keys.addOutputComponent(
     119             :     "z",
     120             :     "default",
     121             :     "Optimal biasing direction; z component."
     122             :   );
     123             : 
     124          34 :   keys.addOutputComponent(
     125             :     "loss",
     126             :     "default",
     127             :     "Loss function value defined by the provided pairing function."
     128             :   );
     129             : 
     130          34 :   keys.addOutputComponent(
     131             :     "sr",
     132             :     "default",
     133             :     "Sampling radius. Reduces sampling to the local proximity of the ligand "
     134             :     "position."
     135             :   );
     136          17 : }
     137             : 
     138           7 : Optimizer::Optimizer(const ActionOptions& ao)
     139             :   : PLUMED_COLVAR_INIT(ao),
     140           7 :     first_step_(true),
     141           7 :     opt_value_(0.0),
     142           7 :     pbc_(true),
     143           7 :     sampling_r_(0.0),
     144           7 :     serial_(false),
     145           7 :     validate_list_(true),
     146           7 :     first_time_(true) {
     147           7 :   parseFlag("SERIAL", serial_);
     148             : 
     149          14 :   if (keywords.exists("LOSS")) {
     150           7 :     std::vector<std::string> loss_labels(0);
     151          14 :     parseVector("LOSS", loss_labels);
     152             : 
     153           7 :     plumed_massert(
     154             :       loss_labels.size() > 0,
     155             :       "maze> Something went wrong with the LOSS keyword.\n"
     156             :     );
     157             : 
     158           7 :     std::string error_msg = "";
     159           7 :     vec_loss_ = tls::get_pointers_labels<Loss*>(
     160             :                   loss_labels,
     161           7 :                   plumed.getActionSet(),
     162             :                   error_msg
     163             :                 );
     164             : 
     165           7 :     if (error_msg.size() > 0) {
     166           0 :       plumed_merror(
     167             :         "maze> Error in the LOSS keyword " + getName() + ": " + error_msg
     168             :       );
     169             :     }
     170             : 
     171           7 :     loss_ = vec_loss_[0];
     172           7 :     log.printf("maze> Loss function linked to the optimizer.\n");
     173           7 :   }
     174             : 
     175          14 :   if (keywords.exists("N_ITER")) {
     176           3 :     parse("N_ITER", n_iter_);
     177             : 
     178           3 :     plumed_massert(
     179             :       n_iter_ > 0,
     180             :       "maze> N_ITER should be explicitly specified and positive.\n"
     181             :     );
     182             : 
     183           3 :     log.printf(
     184             :       "maze> Optimizer will run %u iterations once launched.\n",
     185             :       n_iter_
     186             :     );
     187             :   }
     188             : 
     189             :   std::vector<AtomNumber> ga_list, gb_list;
     190           7 :   parseAtomList("LIGAND", ga_list);
     191           7 :   parseAtomList("PROTEIN", gb_list);
     192             : 
     193           7 :   bool nopbc = !pbc_;
     194           7 :   parseFlag("NOPBC", nopbc);
     195             : 
     196           7 :   bool do_pair = false;
     197           7 :   parseFlag("PAIR", do_pair);
     198             : 
     199           7 :   nl_stride_ = 0;
     200           7 :   bool do_neigh = false;
     201           7 :   parseFlag("NLIST", do_neigh);
     202             : 
     203           7 :   if (do_neigh) {
     204          14 :     if (keywords.exists("NL_CUTOFF")) {
     205           7 :       parse("NL_CUTOFF", nl_cutoff_);
     206             : 
     207           7 :       plumed_massert(
     208             :         nl_cutoff_ > 0,
     209             :         "maze> NL_CUTOFF should be explicitly specified and positive.\n"
     210             :       );
     211             :     }
     212             : 
     213          14 :     if (keywords.exists("NL_STRIDE")) {
     214           7 :       parse("NL_STRIDE", nl_stride_);
     215             : 
     216           7 :       plumed_massert(
     217             :         nl_stride_ > 0,
     218             :         "maze> NL_STRIDE should be explicitly specified and positive.\n"
     219             :       );
     220             :     }
     221             :   }
     222             : 
     223           7 :   if (gb_list.size() > 0) {
     224           7 :     if (do_neigh) {
     225           7 :       neighbor_list_ = Tools::make_unique<NeighborList>(
     226             :                          ga_list,
     227             :                          gb_list,
     228             :                          serial_,
     229             :                          do_pair,
     230           7 :                          pbc_,
     231             :                          getPbc(),
     232             :                          comm,
     233           7 :                          nl_cutoff_,
     234           7 :                          nl_stride_
     235             :                        );
     236             :     } else {
     237           0 :       neighbor_list_=Tools::make_unique<NeighborList>(
     238             :                        ga_list,
     239             :                        gb_list,
     240             :                        serial_,
     241             :                        do_pair,
     242           0 :                        pbc_,
     243             :                        getPbc(),
     244             :                        comm
     245             :                      );
     246             :     }
     247             :   } else {
     248           0 :     if (do_neigh) {
     249           0 :       neighbor_list_ = Tools::make_unique<NeighborList>(
     250             :                          ga_list,
     251             :                          serial_,
     252           0 :                          pbc_,
     253             :                          getPbc(),
     254             :                          comm,
     255           0 :                          nl_cutoff_,
     256           0 :                          nl_stride_
     257             :                        );
     258             :     } else {
     259           0 :       neighbor_list_=Tools::make_unique<NeighborList>(
     260             :                        ga_list,
     261             :                        serial_,
     262           0 :                        pbc_,
     263             :                        getPbc(),
     264             :                        comm
     265             :                      );
     266             :     }
     267             :   }
     268             : 
     269           7 :   requestAtoms(neighbor_list_->getFullAtomList());
     270             : 
     271           7 :   log.printf(
     272             :     "maze> Loss will be calculated between two groups of %u and %u atoms.\n",
     273             :     static_cast<unsigned>(ga_list.size()),
     274             :     static_cast<unsigned>(gb_list.size())
     275             :   );
     276             : 
     277           7 :   log.printf(
     278             :     "maze> First group (LIGAND): from %d to %d.\n",
     279             :     ga_list[0].serial(),
     280             :     ga_list[ga_list.size()-1].serial()
     281             :   );
     282             : 
     283           7 :   if (gb_list.size() > 0) {
     284           7 :     log.printf(
     285             :       "maze> Second group (PROTEIN): from %d to %d.\n",
     286             :       gb_list[0].serial(),
     287             :       gb_list[gb_list.size()-1].serial()
     288             :     );
     289             :   }
     290             : 
     291           7 :   if (pbc_) {
     292           7 :     log.printf("maze> Using periodic boundary conditions.\n");
     293             :   } else {
     294           0 :     log.printf("maze> Without periodic boundary conditions.\n");
     295             :   }
     296             : 
     297           7 :   if (do_pair) {
     298           0 :     log.printf("maze> With PAIR option.\n");
     299             :   }
     300             : 
     301           7 :   if (do_neigh) {
     302           7 :     log.printf(
     303             :       "maze> Using neighbor lists updated every %d steps and cutoff %f.\n",
     304             :       nl_stride_,
     305             :       nl_cutoff_
     306             :     );
     307             :   }
     308             : 
     309             :   // OpenMP
     310           7 :   stride_ = comm.Get_size();
     311           7 :   rank_ = comm.Get_rank();
     312             : 
     313           7 :   n_threads_ = OpenMP::getNumThreads();
     314           7 :   unsigned int nn = neighbor_list_->size();
     315             : 
     316           7 :   if (n_threads_ * stride_ * 10 > nn) {
     317           0 :     n_threads_ = nn / stride_ / 10;
     318             :   }
     319             : 
     320           7 :   if (n_threads_ == 0) {
     321           0 :     n_threads_ = 1;
     322             :   }
     323             : 
     324          14 :   if (keywords.exists("OPTIMIZER_STRIDE")) {
     325           7 :     parse("OPTIMIZER_STRIDE", optimizer_stride_);
     326             : 
     327           7 :     plumed_massert(
     328             :       optimizer_stride_,
     329             :       "maze> OPTIMIZER_STRIDE should be explicitly specified and positive.\n"
     330             :     );
     331             : 
     332           7 :     log.printf(
     333             :       "maze> Launching optimization every %u steps.\n",
     334             :       optimizer_stride_
     335             :     );
     336             :   }
     337             : 
     338           7 :   rnd::randomize();
     339             : 
     340           7 :   opt_.zero();
     341             : 
     342          14 :   addComponentWithDerivatives("x");
     343          14 :   componentIsNotPeriodic("x");
     344             : 
     345          14 :   addComponentWithDerivatives("y");
     346          14 :   componentIsNotPeriodic("y");
     347             : 
     348          14 :   addComponentWithDerivatives("z");
     349          14 :   componentIsNotPeriodic("z");
     350             : 
     351          14 :   addComponent("loss");
     352          14 :   componentIsNotPeriodic("loss");
     353             : 
     354          14 :   addComponent("sr");
     355           7 :   componentIsNotPeriodic("sr");
     356             : 
     357           7 :   value_x_ = getPntrToComponent("x");
     358           7 :   value_y_ = getPntrToComponent("y");
     359           7 :   value_z_ = getPntrToComponent("z");
     360           7 :   value_action_ = getPntrToComponent("loss");
     361           7 :   value_sampling_radius_ = getPntrToComponent("sr");
     362           7 : }
     363             : 
     364    14974890 : double Optimizer::pairing(double distance) const {
     365    14974890 :   return loss_->pairing(distance);
     366             : }
     367             : 
     368           6 : Vector Optimizer::center_of_mass() const {
     369           6 :   const unsigned nl_size = neighbor_list_->size();
     370             : 
     371           6 :   Vector center_of_mass;
     372           6 :   center_of_mass.zero();
     373             :   double mass = 0;
     374             : 
     375      189654 :   for (unsigned int i = 0; i < nl_size; ++i) {
     376      189648 :     unsigned int i0 = neighbor_list_->getClosePair(i).first;
     377      189648 :     center_of_mass += getPosition(i0) * getMass(i0);
     378      189648 :     mass += getMass(i0);
     379             :   }
     380             : 
     381           6 :   return center_of_mass / mass;
     382             : }
     383             : 
     384         210 : void Optimizer::prepare() {
     385         210 :   if (neighbor_list_->getStride() > 0) {
     386         210 :     if (first_time_ || (getStep() % neighbor_list_->getStride() == 0)) {
     387           7 :       requestAtoms(neighbor_list_->getFullAtomList());
     388             : 
     389           7 :       validate_list_ = true;
     390           7 :       first_time_ = false;
     391             :     } else {
     392         203 :       requestAtoms(neighbor_list_->getReducedAtomList());
     393             : 
     394         203 :       validate_list_ = false;
     395             : 
     396         203 :       if (getExchangeStep()) {
     397           0 :         plumed_merror(
     398             :           "maze> Neighbor lists should be updated on exchange steps -- choose "
     399             :           "an NL_STRIDE which divides the exchange stride.\n");
     400             :       }
     401             :     }
     402             : 
     403         210 :     if (getExchangeStep()) {
     404           0 :       first_time_ = true;
     405             :     }
     406             :   }
     407         210 : }
     408             : 
     409         226 : double Optimizer::score() {
     410         226 :   const unsigned nl_size = neighbor_list_->size();
     411         226 :   Vector distance;
     412             :   double function = 0;
     413             : 
     414         226 :   #pragma omp parallel num_threads(n_threads_)
     415             :   {
     416             :     #pragma omp for reduction(+:function)
     417             :     for(unsigned int i = 0; i < nl_size; i++) {
     418             :       unsigned i0 = neighbor_list_->getClosePair(i).first;
     419             :       unsigned i1 = neighbor_list_->getClosePair(i).second;
     420             : 
     421             :       if (getAbsoluteIndex(i0) == getAbsoluteIndex(i1)) {
     422             :         continue;
     423             :       }
     424             : 
     425             :       if (pbc_) {
     426             :         distance = pbcDistance(getPosition(i0), getPosition(i1));
     427             :       } else {
     428             :         distance = delta(getPosition(i0), getPosition(i1));
     429             :       }
     430             : 
     431             :       function += pairing(distance.modulo());
     432             :     }
     433             :   }
     434             : 
     435         226 :   return function;
     436             : }
     437             : 
     438         210 : void Optimizer::update_nl() {
     439         210 :   if (neighbor_list_->getStride() > 0 && validate_list_) {
     440           7 :     neighbor_list_->update(getPositions());
     441             :   }
     442         210 : }
     443             : 
     444         323 : double Optimizer::sampling_radius() {
     445         323 :   const unsigned nl_size=neighbor_list_->size();
     446         323 :   Vector d;
     447             :   double min=std::numeric_limits<int>::max();
     448             : 
     449     8421527 :   for (unsigned int i = 0; i < nl_size; ++i) {
     450     8421204 :     unsigned i0 = neighbor_list_->getClosePair(i).first;
     451     8421204 :     unsigned i1 = neighbor_list_->getClosePair(i).second;
     452             : 
     453     8421204 :     if (getAbsoluteIndex(i0) == getAbsoluteIndex(i1)) {
     454           0 :       continue;
     455             :     }
     456             : 
     457     8421204 :     if (pbc_) {
     458     8421204 :       d = pbcDistance(getPosition(i0), getPosition(i1));
     459             :     } else {
     460           0 :       d = delta(getPosition(i0), getPosition(i1));
     461             :     }
     462             : 
     463     8421204 :     double dist = d.modulo();
     464             : 
     465     8421204 :     if(dist < min) {
     466             :       min = dist;
     467             :     }
     468             :   }
     469             : 
     470         323 :   return min;
     471             : }
     472             : 
     473         210 : void Optimizer::calculate() {
     474         210 :   update_nl();
     475             : 
     476         210 :   if (getStep() % optimizer_stride_ == 0 && !first_step_) {
     477          19 :     optimize();
     478             : 
     479          19 :     value_x_->set(opt_[0]);
     480          19 :     value_y_->set(opt_[1]);
     481          19 :     value_z_->set(opt_[2]);
     482             : 
     483          19 :     value_action_->set(score());
     484          19 :     value_sampling_radius_->set(sampling_radius());
     485             :   } else {
     486         191 :     first_step_=false;
     487             : 
     488         191 :     value_x_->set(opt_[0]);
     489         191 :     value_y_->set(opt_[1]);
     490         191 :     value_z_->set(opt_[2]);
     491             : 
     492         191 :     value_action_->set(score());
     493         191 :     value_sampling_radius_->set(sampling_radius());
     494             :   }
     495         210 : }
     496             : 
     497             : } // namespace maze
     498             : } // namespace PLMD