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_Bias.cpp
25 : *
26 : * @author J. Rydzewski (jr@fizyka.umk.pl)
27 : *
28 : * @code
29 : * @article{rydzewski2018finding,
30 : * title={Finding Multiple Reaction Pathways of Ligand Unbinding},
31 : * author={Rydzewski, J and Valsson, O},
32 : * journal={arXiv preprint arXiv:1808.08089},
33 : * year={2018}
34 : * }
35 : * @endcode
36 : */
37 :
38 : #include "core/ActionRegister.h"
39 : #include "core/PlumedMain.h"
40 :
41 : #include "bias/Bias.h"
42 :
43 : #include "Optimizer.h"
44 : #include "Tools.h"
45 :
46 : namespace PLMD {
47 : namespace maze {
48 :
49 : //+PLUMEDOC MAZE_BIAS MAZE_OPTIMIZER_BIAS
50 : /*
51 :
52 : Biases the ligand along the direction calculated by the chosen MAZE_OPTIMIZER.
53 :
54 : OptimizerBias is a class deriving from Bias, and it is used to adaptively
55 : bias a ligand-protein system toward an optimal solution found by the chosen
56 : MAZE_OPTIMIZER (see [the module page](module_maze.md) and select MAZE_OPTIMIZER from the tags dropdown for the
57 : list of permitted optimizers).
58 :
59 : Remember to define the loss function ([MAZE_LOSS](MAZE_LOSS.md)) and the optimizer
60 : prior to the adaptive bias for the optimizer.
61 :
62 : The adaptive bias potential is the following:
63 :
64 : $$
65 : V({\bf x}_t)=\alpha
66 : \left(wt -
67 : ({\bf x} - {\bf x}^*_{t-\tau})
68 : \cdot
69 : \frac{{\bf x}^*_t - {\bf x}_t}{\|{\bf x}^*_t-{\bf x}_t\|}
70 : \right)^2,
71 : $$
72 :
73 : where ${\bf x}^*_t$ is the optimal solution at time $t$, $w$ is the
74 : biasing rate, $\tau$ is the interval at which the loss function is minimized,
75 : and $\alpha$ is a scaled force constant.
76 :
77 : ## Examples
78 :
79 : In the following example the bias potential biases a ligand atom (which have to be
80 : given as an argument) with the biasing rate equal to 0.02 A/ps, and the biasing
81 : constant equal to 3.6 kcal/(mol A). It also takes an optimizer (see
82 : [the module page](module_maze.md) and select MAZE_OPTIMIZER from the tags dropdown for the
83 : list of permitted optimizers).
84 :
85 : ```plumed
86 : UNITS LENGTH=A TIME=ps ENERGY=kcal/mol
87 :
88 : p: POSITION ATOM=2635 NOPBC
89 :
90 : MAZE_OPTIMIZER_BIAS ...
91 : LABEL=bias
92 :
93 : ARG=p.x,p.y,p.z
94 :
95 : BIASING_RATE=0.02
96 : ALPHA=3.6
97 :
98 : OPTIMIZER=opt
99 : ... MAZE_OPTIMIZER_BIAS
100 : ```
101 :
102 : */
103 : //+ENDPLUMEDOC
104 :
105 : /**
106 : * @class OptimizerBias OptimizerBias.cpp "maze/OptimizerBias.cpp"
107 : *
108 : * @brief Adaptive bias for the maze optimizers.
109 : *
110 : * OptimizerBias is a class deriving from Bias, and it is used to adaptively
111 : * bias a system toward an optimal solution found by an optimizer.
112 : */
113 : class OptimizerBias: public bias::Bias {
114 : public:
115 : /**
116 : * Standard PLUMED2 constructor.
117 : *
118 : * @param ao ActionOptions&.
119 : */
120 : explicit OptimizerBias(const ActionOptions& ao);
121 :
122 : /**
123 : * Destructor.
124 : */
125 3 : ~OptimizerBias() { /* Nothing to do. */ }
126 :
127 : /**
128 : * Register PLUMED2 keywords.
129 : *
130 : * @param keys Keywords.
131 : */
132 : static void registerKeywords(Keywords& keys);
133 :
134 : /**
135 : * Calculate the adaptive biasing potential for ligand unbinding.
136 : */
137 : void calculate();
138 :
139 : private:
140 : /**
141 : * Biased collective variable with Cartesian components, i.e., position,
142 : * center of mass.
143 : */
144 : std::vector<Value*> args_;
145 :
146 : /**
147 : * Pointer to the optimizer used to minimize the collective variable for
148 : * ligand unbinding.
149 : */
150 : Optimizer* optimizer_;
151 : std::vector<Optimizer*> opt_pntrs_;
152 :
153 : //! Adaptive bias potential and the corresponding force.
154 : double bias_;
155 : double force_;
156 :
157 : /*
158 : * Parameters of the adaptive biasing potential:
159 : * alpha_ rescaled force constant
160 : * biasing_speed biasing rate
161 : * biasing_stride biasing stride
162 : * biasing_direction biasing direction
163 : */
164 :
165 : //! Rescaled force constant.
166 : double alpha_;
167 : //! Biasing rate.
168 : double biasing_speed_;
169 : //! Biasing stride.
170 : int biasing_stride_;
171 :
172 : /**
173 : * Biasing direction is approximated by an optimal solution found by an
174 : * optimizer.
175 : */
176 : Vector biasing_direction_;
177 :
178 : //! Total distance traveled by biased atoms.
179 : double total_distance_;
180 :
181 : //! Previous value of the collective variable.
182 : Vector cv0_;
183 :
184 : /*
185 : * Pointers to PLUMED2 output components.
186 : */
187 :
188 : //! Biased collective variable components.
189 : Value* value_dir_x_;
190 : Value* value_dir_y_;
191 : Value* value_dir_z_;
192 :
193 : //! Values of the bias and its force.
194 : Value* value_bias_;
195 : Value* value_force_;
196 :
197 : //! Total distance.
198 : Value* value_total_distance_;
199 : };
200 :
201 : // Register OPTIMIZER_BIAS as a keyword for PLUMED2 input files.
202 : PLUMED_REGISTER_ACTION(OptimizerBias, "MAZE_OPTIMIZER_BIAS")
203 :
204 3 : void OptimizerBias::registerKeywords(Keywords& keys) {
205 3 : Bias::registerKeywords(keys);
206 :
207 3 : keys.add(
208 : "compulsory",
209 : "BIASING_RATE",
210 : "Biasing rate."
211 : );
212 :
213 3 : keys.add(
214 : "compulsory",
215 : "ALPHA",
216 : "Rescaled force constant."
217 : );
218 :
219 3 : keys.add(
220 : "compulsory",
221 : "OPTIMIZER",
222 : "Optimization technique to minimize the collective variable for ligand\
223 : unbinding: RANDOM_WALK,\
224 : STEERED_MD,\
225 : RANDOM_ACCELERATION_MD,\
226 : SIMULATED_ANNEALING,\
227 : MEMETIC_SAMPLING"
228 : );
229 :
230 6 : keys.addOutputComponent(
231 : "force2",
232 : "default",
233 : "scalar",
234 : "Square of the biasing force."
235 : );
236 :
237 6 : keys.addOutputComponent(
238 : "x",
239 : "default",
240 : "scalar",
241 : "Optimal biasing direction: x component."
242 : );
243 :
244 6 : keys.addOutputComponent(
245 : "y",
246 : "default",
247 : "scalar",
248 : "Optimal biasing direction: y component."
249 : );
250 :
251 6 : keys.addOutputComponent(
252 : "z",
253 : "default",
254 : "scalar",
255 : "Optimal biasing direction: z component."
256 : );
257 :
258 6 : keys.addOutputComponent(
259 : "tdist",
260 : "default",
261 : "scalar",
262 : "Total distance traveled by biased atoms."
263 : );
264 3 : }
265 :
266 1 : OptimizerBias::OptimizerBias(const ActionOptions& ao)
267 : : PLUMED_BIAS_INIT(ao),
268 1 : bias_(0.0),
269 1 : force_(0.0),
270 1 : total_distance_(0.0) {
271 1 : log.printf(
272 : "maze> You are using the maze module of PLUMED2,\
273 : please read and cite "
274 : );
275 :
276 2 : log << plumed.cite("Rydzewski J. and Valsson O., arXiv:1808.08089, 2018");
277 1 : log.printf("\n");
278 :
279 1 : args_ = getArguments();
280 1 : log.printf(
281 : "maze> Number of args %zu\n",
282 : args_.size()
283 : );
284 :
285 1 : if (!args_.empty()) {
286 1 : log.printf("maze> With arguments");
287 4 : for (unsigned i = 0; i < args_.size(); i++) {
288 3 : log.printf(" %s", args_[i]->getName().c_str());
289 : }
290 1 : log.printf("\n");
291 : }
292 :
293 1 : if (keywords.exists("ALPHA")) {
294 1 : parse("ALPHA", alpha_);
295 :
296 1 : plumed_massert(
297 : alpha_>0,
298 : "maze> ALPHA should be explicitly specified and positive.\n"
299 : );
300 :
301 1 : log.printf(
302 : "maze> ALPHA read: %f [kcal/mol/A].\n",
303 : alpha_
304 : );
305 : }
306 :
307 1 : if (keywords.exists("BIASING_RATE")) {
308 1 : parse("BIASING_RATE", biasing_speed_);
309 :
310 1 : plumed_massert(
311 : biasing_speed_>0,
312 : "maze> BIASING_RATE should be explicitly specified and positive.\n"
313 : );
314 :
315 1 : log.printf(
316 : "maze> BIASING_RATE read: %f [A/ps].\n",
317 : biasing_speed_
318 : );
319 : }
320 :
321 1 : if (keywords.exists("OPTIMIZER")) {
322 1 : std::vector<std::string> opt_labels(0);
323 2 : parseVector("OPTIMIZER", opt_labels);
324 :
325 1 : plumed_massert(
326 : opt_labels.size() > 0,
327 : "maze> Problem with OPTIMIZER keyword.\n"
328 : );
329 :
330 1 : std::string error_msg = "";
331 2 : opt_pntrs_ = tls::get_pointers_labels<Optimizer*>(
332 : opt_labels,
333 1 : plumed.getActionSet(),
334 : error_msg
335 : );
336 :
337 1 : if (error_msg.size() > 0) {
338 0 : plumed_merror(
339 : "maze> Error in keyword OPTIMIZER of " + getName() + ": " + error_msg
340 : );
341 : }
342 :
343 1 : optimizer_ = opt_pntrs_[0];
344 1 : log.printf(
345 : "maze> Optimizer linked: %s.\n",
346 0 : optimizer_->get_label().c_str()
347 : );
348 :
349 1 : biasing_stride_=optimizer_->get_optimizer_stride();
350 1 : }
351 :
352 1 : checkRead();
353 :
354 2 : addComponent("force2");
355 2 : componentIsNotPeriodic("force2");
356 :
357 2 : addComponent("x");
358 2 : componentIsNotPeriodic("x");
359 :
360 2 : addComponent("y");
361 2 : componentIsNotPeriodic("y");
362 :
363 2 : addComponent("z");
364 2 : componentIsNotPeriodic("z");
365 :
366 2 : addComponent("tdist");
367 2 : componentIsNotPeriodic("tdist");
368 :
369 : biasing_direction_.zero();
370 : cv0_.zero();
371 :
372 1 : value_bias_ = getPntrToComponent("bias");
373 1 : value_force_ = getPntrToComponent("force2");
374 :
375 1 : value_dir_x_ = getPntrToComponent("x");
376 1 : value_dir_y_ = getPntrToComponent("y");
377 1 : value_dir_z_ = getPntrToComponent("z");
378 :
379 1 : value_total_distance_=getPntrToComponent("tdist");
380 1 : }
381 :
382 30 : void OptimizerBias::calculate() {
383 : // Unpack arguments and optimizers.
384 : Vector cv(
385 : args_[0]->get(),
386 : args_[1]->get(),
387 : args_[2]->get()
388 30 : );
389 :
390 : Vector opt_direction(
391 30 : optimizer_->value_x_->get(),
392 30 : optimizer_->value_y_->get(),
393 30 : optimizer_->value_z_->get()
394 30 : );
395 :
396 30 : if (getStep() == 0) {
397 1 : cv0_=cv;
398 : }
399 :
400 : /*
401 : * For details see a paper by Rydzewski and Valsson.
402 : */
403 : double dot = dotProduct(cv - cv0_, biasing_direction_);
404 30 : double delta_cv = biasing_speed_ * getTime() - (dot + total_distance_);
405 :
406 30 : double sign = tls::sgn(delta_cv);
407 :
408 30 : bias_ = alpha_ * delta_cv * delta_cv;
409 30 : force_ = 2.0 * sign * alpha_ * fabs(delta_cv);
410 :
411 30 : if (getStep() % biasing_stride_ == 0) {
412 3 : biasing_direction_ = opt_direction;
413 3 : cv0_ = cv;
414 3 : total_distance_ += dot;
415 : }
416 :
417 : /*
418 : * Return the biasing force to MD engine.
419 : */
420 30 : setOutputForce(0, force_ * biasing_direction_[0]);
421 30 : setOutputForce(1, force_ * biasing_direction_[1]);
422 30 : setOutputForce(2, force_ * biasing_direction_[2]);
423 :
424 : /*
425 : * Set values for PLUMED2 outputs.
426 : */
427 30 : value_bias_->set(bias_);
428 30 : value_force_->set(force_);
429 :
430 30 : value_total_distance_->set(total_distance_);
431 :
432 30 : value_dir_x_->set(biasing_direction_[0]);
433 30 : value_dir_y_->set(biasing_direction_[1]);
434 30 : value_dir_z_->set(biasing_direction_[2]);
435 30 : }
436 :
437 : } // namespace maze
438 : } // namespace PLMD
|
