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