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 : l: MAZE_LOSS PARAMS=1,1,1
91 :
92 : opt: MAZE_RANDOM_WALK LOSS=l OPTIMIZER_STRIDE=200 LIGAND=2635-2646 PROTEIN=1-2634
93 :
94 : MAZE_OPTIMIZER_BIAS ...
95 : LABEL=bias
96 :
97 : ARG=p.x,p.y,p.z
98 :
99 : BIASING_RATE=0.02
100 : ALPHA=3.6
101 :
102 : OPTIMIZER=opt
103 : ... MAZE_OPTIMIZER_BIAS
104 : ```
105 :
106 : */
107 : //+ENDPLUMEDOC
108 :
109 : /**
110 : * @class OptimizerBias OptimizerBias.cpp "maze/OptimizerBias.cpp"
111 : *
112 : * @brief Adaptive bias for the maze optimizers.
113 : *
114 : * OptimizerBias is a class deriving from Bias, and it is used to adaptively
115 : * bias a system toward an optimal solution found by an optimizer.
116 : */
117 : class OptimizerBias: public bias::Bias {
118 : public:
119 : /**
120 : * Standard PLUMED2 constructor.
121 : *
122 : * @param ao ActionOptions&.
123 : */
124 : explicit OptimizerBias(const ActionOptions& ao);
125 :
126 : /**
127 : * Destructor.
128 : */
129 3 : ~OptimizerBias() { /* Nothing to do. */ }
130 :
131 : /**
132 : * Register PLUMED2 keywords.
133 : *
134 : * @param keys Keywords.
135 : */
136 : static void registerKeywords(Keywords& keys);
137 :
138 : /**
139 : * Calculate the adaptive biasing potential for ligand unbinding.
140 : */
141 : void calculate();
142 :
143 : private:
144 : /**
145 : * Biased collective variable with Cartesian components, i.e., position,
146 : * center of mass.
147 : */
148 : std::vector<Value*> args_;
149 :
150 : /**
151 : * Pointer to the optimizer used to minimize the collective variable for
152 : * ligand unbinding.
153 : */
154 : Optimizer* optimizer_;
155 : std::vector<Optimizer*> opt_pntrs_;
156 :
157 : //! Adaptive bias potential and the corresponding force.
158 : double bias_;
159 : double force_;
160 :
161 : /*
162 : * Parameters of the adaptive biasing potential:
163 : * alpha_ rescaled force constant
164 : * biasing_speed biasing rate
165 : * biasing_stride biasing stride
166 : * biasing_direction biasing direction
167 : */
168 :
169 : //! Rescaled force constant.
170 : double alpha_;
171 : //! Biasing rate.
172 : double biasing_speed_;
173 : //! Biasing stride.
174 : int biasing_stride_;
175 :
176 : /**
177 : * Biasing direction is approximated by an optimal solution found by an
178 : * optimizer.
179 : */
180 : Vector biasing_direction_;
181 :
182 : //! Total distance traveled by biased atoms.
183 : double total_distance_;
184 :
185 : //! Previous value of the collective variable.
186 : Vector cv0_;
187 :
188 : /*
189 : * Pointers to PLUMED2 output components.
190 : */
191 :
192 : //! Biased collective variable components.
193 : Value* value_dir_x_;
194 : Value* value_dir_y_;
195 : Value* value_dir_z_;
196 :
197 : //! Values of the bias and its force.
198 : Value* value_bias_;
199 : Value* value_force_;
200 :
201 : //! Total distance.
202 : Value* value_total_distance_;
203 : };
204 :
205 : // Register OPTIMIZER_BIAS as a keyword for PLUMED2 input files.
206 : PLUMED_REGISTER_ACTION(OptimizerBias, "MAZE_OPTIMIZER_BIAS")
207 :
208 6 : void OptimizerBias::registerKeywords(Keywords& keys) {
209 6 : Bias::registerKeywords(keys);
210 :
211 6 : keys.add(
212 : "compulsory",
213 : "BIASING_RATE",
214 : "Biasing rate."
215 : );
216 :
217 6 : keys.add(
218 : "compulsory",
219 : "ALPHA",
220 : "Rescaled force constant."
221 : );
222 :
223 6 : keys.add(
224 : "compulsory",
225 : "OPTIMIZER",
226 : "Optimization technique to minimize the collective variable for ligand\
227 : unbinding: RANDOM_WALK,\
228 : STEERED_MD,\
229 : RANDOM_ACCELERATION_MD,\
230 : SIMULATED_ANNEALING,\
231 : MEMETIC_SAMPLING"
232 : );
233 :
234 12 : keys.addOutputComponent(
235 : "force2",
236 : "default",
237 : "scalar",
238 : "Square of the biasing force."
239 : );
240 :
241 12 : keys.addOutputComponent(
242 : "x",
243 : "default",
244 : "scalar",
245 : "Optimal biasing direction: x component."
246 : );
247 :
248 12 : keys.addOutputComponent(
249 : "y",
250 : "default",
251 : "scalar",
252 : "Optimal biasing direction: y component."
253 : );
254 :
255 12 : keys.addOutputComponent(
256 : "z",
257 : "default",
258 : "scalar",
259 : "Optimal biasing direction: z component."
260 : );
261 :
262 12 : keys.addOutputComponent(
263 : "tdist",
264 : "default",
265 : "scalar",
266 : "Total distance traveled by biased atoms."
267 : );
268 6 : }
269 :
270 2 : OptimizerBias::OptimizerBias(const ActionOptions& ao)
271 : : PLUMED_BIAS_INIT(ao),
272 2 : bias_(0.0),
273 2 : force_(0.0),
274 2 : total_distance_(0.0) {
275 2 : log.printf(
276 : "maze> You are using the maze module of PLUMED2,\
277 : please read and cite "
278 : );
279 :
280 2 : log << plumed.cite("Rydzewski J. and Valsson O., arXiv:1808.08089, 2018");
281 2 : log.printf("\n");
282 :
283 2 : args_ = getArguments();
284 2 : log.printf(
285 : "maze> Number of args %zu\n",
286 : args_.size()
287 : );
288 :
289 2 : if (!args_.empty()) {
290 2 : log.printf("maze> With arguments");
291 8 : for (unsigned i = 0; i < args_.size(); i++) {
292 6 : log.printf(" %s", args_[i]->getName().c_str());
293 : }
294 2 : log.printf("\n");
295 : }
296 :
297 2 : if (keywords.exists("ALPHA")) {
298 2 : parse("ALPHA", alpha_);
299 :
300 2 : plumed_massert(
301 : alpha_>0,
302 : "maze> ALPHA should be explicitly specified and positive.\n"
303 : );
304 :
305 2 : log.printf(
306 : "maze> ALPHA read: %f [kcal/mol/A].\n",
307 : alpha_
308 : );
309 : }
310 :
311 2 : if (keywords.exists("BIASING_RATE")) {
312 2 : parse("BIASING_RATE", biasing_speed_);
313 :
314 2 : plumed_massert(
315 : biasing_speed_>0,
316 : "maze> BIASING_RATE should be explicitly specified and positive.\n"
317 : );
318 :
319 2 : log.printf(
320 : "maze> BIASING_RATE read: %f [A/ps].\n",
321 : biasing_speed_
322 : );
323 : }
324 :
325 2 : if (keywords.exists("OPTIMIZER")) {
326 2 : std::vector<std::string> opt_labels(0);
327 4 : parseVector("OPTIMIZER", opt_labels);
328 :
329 2 : plumed_massert(
330 : opt_labels.size() > 0,
331 : "maze> Problem with OPTIMIZER keyword.\n"
332 : );
333 :
334 2 : std::string error_msg = "";
335 4 : opt_pntrs_ = tls::get_pointers_labels<Optimizer*>(
336 : opt_labels,
337 2 : plumed.getActionSet(),
338 : error_msg
339 2 : );
340 :
341 2 : if (error_msg.size() > 0) {
342 4 : plumed_merror(
343 : "maze> Error in keyword OPTIMIZER of " + getName() + ": " + error_msg
344 : );
345 : }
346 :
347 1 : optimizer_ = opt_pntrs_[0];
348 1 : log.printf(
349 : "maze> Optimizer linked: %s.\n",
350 1 : optimizer_->get_label().c_str()
351 : );
352 :
353 1 : biasing_stride_=optimizer_->get_optimizer_stride();
354 2 : }
355 :
356 1 : checkRead();
357 :
358 2 : addComponent("force2");
359 2 : componentIsNotPeriodic("force2");
360 :
361 2 : addComponent("x");
362 2 : componentIsNotPeriodic("x");
363 :
364 2 : addComponent("y");
365 2 : componentIsNotPeriodic("y");
366 :
367 2 : addComponent("z");
368 2 : componentIsNotPeriodic("z");
369 :
370 2 : addComponent("tdist");
371 2 : componentIsNotPeriodic("tdist");
372 :
373 : biasing_direction_.zero();
374 : cv0_.zero();
375 :
376 1 : value_bias_ = getPntrToComponent("bias");
377 1 : value_force_ = getPntrToComponent("force2");
378 :
379 1 : value_dir_x_ = getPntrToComponent("x");
380 1 : value_dir_y_ = getPntrToComponent("y");
381 1 : value_dir_z_ = getPntrToComponent("z");
382 :
383 2 : value_total_distance_=getPntrToComponent("tdist");
384 3 : }
385 :
386 30 : void OptimizerBias::calculate() {
387 : // Unpack arguments and optimizers.
388 : Vector cv(
389 : args_[0]->get(),
390 : args_[1]->get(),
391 : args_[2]->get()
392 30 : );
393 :
394 : Vector opt_direction(
395 30 : optimizer_->value_x_->get(),
396 30 : optimizer_->value_y_->get(),
397 30 : optimizer_->value_z_->get()
398 30 : );
399 :
400 30 : if (getStep() == 0) {
401 1 : cv0_=cv;
402 : }
403 :
404 : /*
405 : * For details see a paper by Rydzewski and Valsson.
406 : */
407 : double dot = dotProduct(cv - cv0_, biasing_direction_);
408 30 : double delta_cv = biasing_speed_ * getTime() - (dot + total_distance_);
409 :
410 30 : double sign = tls::sgn(delta_cv);
411 :
412 30 : bias_ = alpha_ * delta_cv * delta_cv;
413 30 : force_ = 2.0 * sign * alpha_ * fabs(delta_cv);
414 :
415 30 : if (getStep() % biasing_stride_ == 0) {
416 3 : biasing_direction_ = opt_direction;
417 3 : cv0_ = cv;
418 3 : total_distance_ += dot;
419 : }
420 :
421 : /*
422 : * Return the biasing force to MD engine.
423 : */
424 30 : setOutputForce(0, force_ * biasing_direction_[0]);
425 30 : setOutputForce(1, force_ * biasing_direction_[1]);
426 30 : setOutputForce(2, force_ * biasing_direction_[2]);
427 :
428 : /*
429 : * Set values for PLUMED2 outputs.
430 : */
431 30 : value_bias_->set(bias_);
432 30 : value_force_->set(force_);
433 :
434 30 : value_total_distance_->set(total_distance_);
435 :
436 30 : value_dir_x_->set(biasing_direction_[0]);
437 30 : value_dir_y_->set(biasing_direction_[1]);
438 30 : value_dir_z_->set(biasing_direction_[2]);
439 30 : }
440 :
441 : } // namespace maze
442 : } // namespace PLMD
|
