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 \ref 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 : 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 10421 : PLUMED_REGISTER_ACTION(OptimizerBias, "MAZE_OPTIMIZER_BIAS")
199 :
200 2 : void OptimizerBias::registerKeywords(Keywords& keys) {
201 2 : Bias::registerKeywords(keys);
202 :
203 2 : keys.use("ARG");
204 :
205 4 : keys.add(
206 : "compulsory",
207 : "BIASING_RATE",
208 : "Biasing rate."
209 : );
210 :
211 4 : keys.add(
212 : "compulsory",
213 : "ALPHA",
214 : "Rescaled force constant."
215 : );
216 :
217 4 : 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 2 : componentsAreNotOptional(keys);
229 :
230 4 : keys.addOutputComponent(
231 : "force2",
232 : "default",
233 : "Square of the biasing force."
234 : );
235 :
236 4 : keys.addOutputComponent(
237 : "x",
238 : "default",
239 : "Optimal biasing direction: x component."
240 : );
241 :
242 4 : keys.addOutputComponent(
243 : "y",
244 : "default",
245 : "Optimal biasing direction: y component."
246 : );
247 :
248 4 : keys.addOutputComponent(
249 : "z",
250 : "default",
251 : "Optimal biasing direction: z component."
252 : );
253 :
254 4 : keys.addOutputComponent(
255 : "tdist",
256 : "default",
257 : "Total distance traveled by biased atoms."
258 : );
259 2 : }
260 :
261 1 : OptimizerBias::OptimizerBias(const ActionOptions& ao)
262 : : PLUMED_BIAS_INIT(ao),
263 1 : bias_(0.0),
264 1 : force_(0.0),
265 1 : total_distance_(0.0)
266 : {
267 1 : log.printf(
268 : "maze> You are using the maze module of PLUMED2,\
269 : please read and cite "
270 : );
271 :
272 2 : log << plumed.cite("Rydzewski J. and Valsson O., arXiv:1808.08089, 2018");
273 1 : log.printf("\n");
274 :
275 1 : args_ = getArguments();
276 1 : log.printf(
277 : "maze> Number of args %zu\n",
278 : args_.size()
279 : );
280 :
281 1 : if (!args_.empty()) {
282 1 : log.printf("maze> With arguments");
283 4 : for (unsigned i = 0; i < args_.size(); i++) {
284 3 : log.printf(" %s", args_[i]->getName().c_str());
285 : }
286 1 : log.printf("\n");
287 : }
288 :
289 2 : if (keywords.exists("ALPHA")) {
290 1 : parse("ALPHA", alpha_);
291 :
292 1 : plumed_massert(
293 : alpha_>0,
294 : "maze> ALPHA should be explicitly specified and positive.\n"
295 : );
296 :
297 1 : log.printf(
298 : "maze> ALPHA read: %f [kcal/mol/A].\n",
299 : alpha_
300 : );
301 : }
302 :
303 2 : if (keywords.exists("BIASING_RATE")) {
304 1 : parse("BIASING_RATE", biasing_speed_);
305 :
306 1 : plumed_massert(
307 : biasing_speed_>0,
308 : "maze> BIASING_RATE should be explicitly specified and positive.\n"
309 : );
310 :
311 1 : log.printf(
312 : "maze> BIASING_RATE read: %f [A/ps].\n",
313 : biasing_speed_
314 : );
315 : }
316 :
317 2 : if (keywords.exists("OPTIMIZER")) {
318 1 : std::vector<std::string> opt_labels(0);
319 2 : parseVector("OPTIMIZER", opt_labels);
320 :
321 1 : plumed_massert(
322 : opt_labels.size() > 0,
323 : "maze> Problem with OPTIMIZER keyword.\n"
324 : );
325 :
326 1 : std::string error_msg = "";
327 2 : opt_pntrs_ = tls::get_pointers_labels<Optimizer*>(
328 : opt_labels,
329 1 : plumed.getActionSet(),
330 : error_msg
331 : );
332 :
333 1 : if (error_msg.size() > 0) {
334 0 : plumed_merror(
335 : "maze> Error in keyword OPTIMIZER of " + getName() + ": " + error_msg
336 : );
337 : }
338 :
339 1 : optimizer_ = opt_pntrs_[0];
340 2 : log.printf(
341 : "maze> Optimizer linked: %s.\n",
342 0 : optimizer_->get_label().c_str()
343 : );
344 :
345 1 : biasing_stride_=optimizer_->get_optimizer_stride();
346 1 : }
347 :
348 1 : checkRead();
349 :
350 1 : addComponent("force2");
351 1 : componentIsNotPeriodic("force2");
352 :
353 1 : addComponent("x");
354 1 : componentIsNotPeriodic("x");
355 :
356 1 : addComponent("y");
357 1 : componentIsNotPeriodic("y");
358 :
359 1 : addComponent("z");
360 1 : componentIsNotPeriodic("z");
361 :
362 1 : addComponent("tdist");
363 1 : componentIsNotPeriodic("tdist");
364 :
365 1 : biasing_direction_.zero();
366 1 : cv0_.zero();
367 :
368 1 : value_bias_ = getPntrToComponent("bias");
369 1 : value_force_ = getPntrToComponent("force2");
370 :
371 1 : value_dir_x_ = getPntrToComponent("x");
372 1 : value_dir_y_ = getPntrToComponent("y");
373 1 : value_dir_z_ = getPntrToComponent("z");
374 :
375 1 : value_total_distance_=getPntrToComponent("tdist");
376 1 : }
377 :
378 30 : void OptimizerBias::calculate() {
379 : // Unpack arguments and optimizers.
380 : Vector cv(
381 : args_[0]->get(),
382 : args_[1]->get(),
383 : args_[2]->get()
384 30 : );
385 :
386 : Vector opt_direction(
387 30 : optimizer_->value_x_->get(),
388 30 : optimizer_->value_y_->get(),
389 30 : optimizer_->value_z_->get()
390 30 : );
391 :
392 30 : if (getStep() == 0) {
393 1 : cv0_=cv;
394 : }
395 :
396 : /*
397 : * For details see a paper by Rydzewski and Valsson.
398 : */
399 30 : double dot = dotProduct(cv - cv0_, biasing_direction_);
400 30 : double delta_cv = biasing_speed_ * getTime() - (dot + total_distance_);
401 :
402 30 : double sign = tls::sgn(delta_cv);
403 :
404 30 : bias_ = alpha_ * delta_cv * delta_cv;
405 30 : force_ = 2.0 * sign * alpha_ * fabs(delta_cv);
406 :
407 30 : if (getStep() % biasing_stride_ == 0) {
408 3 : biasing_direction_ = opt_direction;
409 3 : cv0_ = cv;
410 3 : total_distance_ += dot;
411 : }
412 :
413 : /*
414 : * Return the biasing force to MD engine.
415 : */
416 30 : setOutputForce(0, force_ * biasing_direction_[0]);
417 30 : setOutputForce(1, force_ * biasing_direction_[1]);
418 30 : setOutputForce(2, force_ * biasing_direction_[2]);
419 :
420 : /*
421 : * Set values for PLUMED2 outputs.
422 : */
423 30 : value_bias_->set(bias_);
424 30 : value_force_->set(force_);
425 :
426 30 : value_total_distance_->set(total_distance_);
427 :
428 30 : value_dir_x_->set(biasing_direction_[0]);
429 30 : value_dir_y_->set(biasing_direction_[1]);
430 30 : value_dir_z_->set(biasing_direction_[2]);
431 30 : }
432 :
433 : } // namespace maze
434 : } // namespace PLMD
|
