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 1 : log.printf(
265 : "maze> You are using the maze module of PLUMED2,\
266 : please read and cite "
267 : );
268 :
269 2 : log << plumed.cite("Rydzewski J. and Valsson O., arXiv:1808.08089, 2018");
270 1 : log.printf("\n");
271 :
272 1 : args_ = getArguments();
273 1 : log.printf(
274 : "maze> Number of args %zu\n",
275 : args_.size()
276 : );
277 :
278 1 : if (!args_.empty()) {
279 1 : log.printf("maze> With arguments");
280 4 : for (unsigned i = 0; i < args_.size(); i++) {
281 3 : log.printf(" %s", args_[i]->getName().c_str());
282 : }
283 1 : log.printf("\n");
284 : }
285 :
286 2 : if (keywords.exists("ALPHA")) {
287 1 : parse("ALPHA", alpha_);
288 :
289 1 : plumed_massert(
290 : alpha_>0,
291 : "maze> ALPHA should be explicitly specified and positive.\n"
292 : );
293 :
294 1 : log.printf(
295 : "maze> ALPHA read: %f [kcal/mol/A].\n",
296 : alpha_
297 : );
298 : }
299 :
300 2 : if (keywords.exists("BIASING_RATE")) {
301 1 : parse("BIASING_RATE", biasing_speed_);
302 :
303 1 : plumed_massert(
304 : biasing_speed_>0,
305 : "maze> BIASING_RATE should be explicitly specified and positive.\n"
306 : );
307 :
308 1 : log.printf(
309 : "maze> BIASING_RATE read: %f [A/ps].\n",
310 : biasing_speed_
311 : );
312 : }
313 :
314 2 : if (keywords.exists("OPTIMIZER")) {
315 1 : std::vector<std::string> opt_labels(0);
316 2 : parseVector("OPTIMIZER", opt_labels);
317 :
318 1 : plumed_massert(
319 : opt_labels.size() > 0,
320 : "maze> Problem with OPTIMIZER keyword.\n"
321 : );
322 :
323 1 : std::string error_msg = "";
324 2 : opt_pntrs_ = tls::get_pointers_labels<Optimizer*>(
325 : opt_labels,
326 1 : plumed.getActionSet(),
327 : error_msg
328 : );
329 :
330 1 : if (error_msg.size() > 0) {
331 0 : plumed_merror(
332 : "maze> Error in keyword OPTIMIZER of " + getName() + ": " + error_msg
333 : );
334 : }
335 :
336 1 : optimizer_ = opt_pntrs_[0];
337 1 : log.printf(
338 : "maze> Optimizer linked: %s.\n",
339 0 : optimizer_->get_label().c_str()
340 : );
341 :
342 1 : biasing_stride_=optimizer_->get_optimizer_stride();
343 1 : }
344 :
345 1 : checkRead();
346 :
347 2 : addComponent("force2");
348 2 : componentIsNotPeriodic("force2");
349 :
350 2 : addComponent("x");
351 2 : componentIsNotPeriodic("x");
352 :
353 2 : addComponent("y");
354 2 : componentIsNotPeriodic("y");
355 :
356 2 : addComponent("z");
357 2 : componentIsNotPeriodic("z");
358 :
359 2 : addComponent("tdist");
360 1 : componentIsNotPeriodic("tdist");
361 :
362 1 : biasing_direction_.zero();
363 1 : cv0_.zero();
364 :
365 1 : value_bias_ = getPntrToComponent("bias");
366 1 : value_force_ = getPntrToComponent("force2");
367 :
368 1 : value_dir_x_ = getPntrToComponent("x");
369 1 : value_dir_y_ = getPntrToComponent("y");
370 1 : value_dir_z_ = getPntrToComponent("z");
371 :
372 1 : value_total_distance_=getPntrToComponent("tdist");
373 1 : }
374 :
375 30 : void OptimizerBias::calculate() {
376 : // Unpack arguments and optimizers.
377 : Vector cv(
378 : args_[0]->get(),
379 : args_[1]->get(),
380 : args_[2]->get()
381 30 : );
382 :
383 : Vector opt_direction(
384 30 : optimizer_->value_x_->get(),
385 30 : optimizer_->value_y_->get(),
386 30 : optimizer_->value_z_->get()
387 30 : );
388 :
389 30 : if (getStep() == 0) {
390 1 : cv0_=cv;
391 : }
392 :
393 : /*
394 : * For details see a paper by Rydzewski and Valsson.
395 : */
396 30 : double dot = dotProduct(cv - cv0_, biasing_direction_);
397 30 : double delta_cv = biasing_speed_ * getTime() - (dot + total_distance_);
398 :
399 30 : double sign = tls::sgn(delta_cv);
400 :
401 30 : bias_ = alpha_ * delta_cv * delta_cv;
402 30 : force_ = 2.0 * sign * alpha_ * fabs(delta_cv);
403 :
404 30 : if (getStep() % biasing_stride_ == 0) {
405 3 : biasing_direction_ = opt_direction;
406 3 : cv0_ = cv;
407 3 : total_distance_ += dot;
408 : }
409 :
410 : /*
411 : * Return the biasing force to MD engine.
412 : */
413 30 : setOutputForce(0, force_ * biasing_direction_[0]);
414 30 : setOutputForce(1, force_ * biasing_direction_[1]);
415 30 : setOutputForce(2, force_ * biasing_direction_[2]);
416 :
417 : /*
418 : * Set values for PLUMED2 outputs.
419 : */
420 30 : value_bias_->set(bias_);
421 30 : value_force_->set(force_);
422 :
423 30 : value_total_distance_->set(total_distance_);
424 :
425 30 : value_dir_x_->set(biasing_direction_[0]);
426 30 : value_dir_y_->set(biasing_direction_[1]);
427 30 : value_dir_z_->set(biasing_direction_[2]);
428 30 : }
429 :
430 : } // namespace maze
431 : } // namespace PLMD
|
