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.cpp
25 : *
26 : * @author J. Rydzewski (jr@fizyka.umk.pl)
27 : */
28 :
29 : #include "Optimizer.h"
30 : #include "core/PlumedMain.h"
31 : #include "tools/Tools.h"
32 :
33 : namespace PLMD {
34 : namespace maze {
35 :
36 17 : void Optimizer::registerKeywords(Keywords& keys) {
37 17 : Colvar::registerKeywords(keys);
38 :
39 17 : keys.addFlag(
40 : "SERIAL",
41 : false,
42 : "Perform the simulation in serial -- used only for debugging purposes, "
43 : "should not be used otherwise."
44 : );
45 :
46 17 : keys.addFlag(
47 : "PAIR",
48 : false,
49 : "Pair only the 1st element of the 1st group with the 1st element in the "
50 : "second, etc."
51 : );
52 :
53 17 : keys.addFlag(
54 : "NLIST",
55 : false,
56 : "Use a neighbor list of ligand-protein atom pairs to speed up the "
57 : "calculating of the distances."
58 : );
59 :
60 17 : keys.add(
61 : "optional",
62 : "NL_CUTOFF",
63 : "Neighbor list cut-off for the distances of ligand-protein atom pairs."
64 : );
65 :
66 17 : keys.add(
67 : "optional",
68 : "NL_STRIDE",
69 : "Update stride for the ligand-protein atom pairs in the neighbor list."
70 : );
71 :
72 17 : keys.add(
73 : "compulsory",
74 : "N_ITER",
75 : "Number of optimization steps. Required only for optimizers, do not pass "
76 : "this keyword to the fake optimizers (results in crash) , e.g., random "
77 : "walk, steered MD, or random acceleration MD."
78 : );
79 :
80 17 : keys.add(
81 : "optional",
82 : "LOSS",
83 : "Loss function describing ligand-protein interactions required by every "
84 : "optimizer."
85 : );
86 :
87 17 : keys.add(
88 : "atoms",
89 : "LIGAND",
90 : "Indices of ligand atoms."
91 : );
92 :
93 17 : keys.add(
94 : "atoms",
95 : "PROTEIN",
96 : "Indices of protein atoms."
97 : );
98 :
99 17 : keys.add(
100 : "compulsory",
101 : "OPTIMIZER_STRIDE",
102 : "Optimizer stride. Sets up a callback function that launches the "
103 : "optimization process every OPTIMIZER_STRIDE."
104 : );
105 :
106 34 : keys.addOutputComponent(
107 : "x",
108 : "default",
109 : "scalar",
110 : "Optimal biasing direction; x component."
111 : );
112 :
113 34 : keys.addOutputComponent(
114 : "y",
115 : "default",
116 : "scalar",
117 : "Optimal biasing direction; y component."
118 : );
119 :
120 34 : keys.addOutputComponent(
121 : "z",
122 : "default",
123 : "scalar",
124 : "Optimal biasing direction; z component."
125 : );
126 :
127 34 : keys.addOutputComponent(
128 : "loss",
129 : "default",
130 : "scalar",
131 : "Loss function value defined by the provided pairing function."
132 : );
133 :
134 34 : keys.addOutputComponent(
135 : "sr",
136 : "default",
137 : "scalar",
138 : "Sampling radius. Reduces sampling to the local proximity of the ligand "
139 : "position."
140 : );
141 17 : }
142 :
143 7 : Optimizer::Optimizer(const ActionOptions& ao)
144 : : PLUMED_COLVAR_INIT(ao),
145 7 : first_step_(true),
146 7 : opt_value_(0.0),
147 7 : pbc_(true),
148 7 : sampling_r_(0.0),
149 7 : serial_(false),
150 7 : validate_list_(true),
151 7 : first_time_(true) {
152 7 : parseFlag("SERIAL", serial_);
153 :
154 7 : if (keywords.exists("LOSS")) {
155 7 : std::vector<std::string> loss_labels(0);
156 14 : parseVector("LOSS", loss_labels);
157 :
158 7 : plumed_massert(
159 : loss_labels.size() > 0,
160 : "maze> Something went wrong with the LOSS keyword.\n"
161 : );
162 :
163 7 : std::string error_msg = "";
164 7 : vec_loss_ = tls::get_pointers_labels<Loss*>(
165 : loss_labels,
166 7 : plumed.getActionSet(),
167 : error_msg
168 : );
169 :
170 7 : if (error_msg.size() > 0) {
171 0 : plumed_merror(
172 : "maze> Error in the LOSS keyword " + getName() + ": " + error_msg
173 : );
174 : }
175 :
176 7 : loss_ = vec_loss_[0];
177 7 : log.printf("maze> Loss function linked to the optimizer.\n");
178 7 : }
179 :
180 7 : if (keywords.exists("N_ITER")) {
181 3 : parse("N_ITER", n_iter_);
182 :
183 3 : plumed_massert(
184 : n_iter_ > 0,
185 : "maze> N_ITER should be explicitly specified and positive.\n"
186 : );
187 :
188 3 : log.printf(
189 : "maze> Optimizer will run %u iterations once launched.\n",
190 : n_iter_
191 : );
192 : }
193 :
194 : std::vector<AtomNumber> ga_list, gb_list;
195 7 : parseAtomList("LIGAND", ga_list);
196 7 : parseAtomList("PROTEIN", gb_list);
197 :
198 7 : bool nopbc = !pbc_;
199 7 : parseFlag("NOPBC", nopbc);
200 :
201 7 : bool do_pair = false;
202 7 : parseFlag("PAIR", do_pair);
203 :
204 7 : nl_stride_ = 0;
205 7 : bool do_neigh = false;
206 7 : parseFlag("NLIST", do_neigh);
207 :
208 7 : if (do_neigh) {
209 7 : if (keywords.exists("NL_CUTOFF")) {
210 7 : parse("NL_CUTOFF", nl_cutoff_);
211 :
212 7 : plumed_massert(
213 : nl_cutoff_ > 0,
214 : "maze> NL_CUTOFF should be explicitly specified and positive.\n"
215 : );
216 : }
217 :
218 7 : if (keywords.exists("NL_STRIDE")) {
219 7 : parse("NL_STRIDE", nl_stride_);
220 :
221 7 : plumed_massert(
222 : nl_stride_ > 0,
223 : "maze> NL_STRIDE should be explicitly specified and positive.\n"
224 : );
225 : }
226 : }
227 :
228 7 : if (gb_list.size() > 0) {
229 7 : if (do_neigh) {
230 7 : neighbor_list_ = Tools::make_unique<NeighborList>(
231 : ga_list,
232 : gb_list,
233 : serial_,
234 : do_pair,
235 7 : pbc_,
236 : getPbc(),
237 : comm,
238 7 : nl_cutoff_,
239 7 : nl_stride_
240 : );
241 : } else {
242 0 : neighbor_list_=Tools::make_unique<NeighborList>(
243 : ga_list,
244 : gb_list,
245 : serial_,
246 : do_pair,
247 0 : pbc_,
248 : getPbc(),
249 : comm
250 : );
251 : }
252 : } else {
253 0 : if (do_neigh) {
254 0 : neighbor_list_ = Tools::make_unique<NeighborList>(
255 : ga_list,
256 : serial_,
257 0 : pbc_,
258 : getPbc(),
259 : comm,
260 0 : nl_cutoff_,
261 0 : nl_stride_
262 : );
263 : } else {
264 0 : neighbor_list_=Tools::make_unique<NeighborList>(
265 : ga_list,
266 : serial_,
267 0 : pbc_,
268 : getPbc(),
269 : comm
270 : );
271 : }
272 : }
273 :
274 7 : requestAtoms(neighbor_list_->getFullAtomList());
275 :
276 7 : log.printf(
277 : "maze> Loss will be calculated between two groups of %u and %u atoms.\n",
278 : static_cast<unsigned>(ga_list.size()),
279 : static_cast<unsigned>(gb_list.size())
280 : );
281 :
282 7 : log.printf(
283 : "maze> First group (LIGAND): from %d to %d.\n",
284 : ga_list[0].serial(),
285 : ga_list[ga_list.size()-1].serial()
286 : );
287 :
288 7 : if (gb_list.size() > 0) {
289 7 : log.printf(
290 : "maze> Second group (PROTEIN): from %d to %d.\n",
291 : gb_list[0].serial(),
292 : gb_list[gb_list.size()-1].serial()
293 : );
294 : }
295 :
296 7 : if (pbc_) {
297 7 : log.printf("maze> Using periodic boundary conditions.\n");
298 : } else {
299 0 : log.printf("maze> Without periodic boundary conditions.\n");
300 : }
301 :
302 7 : if (do_pair) {
303 0 : log.printf("maze> With PAIR option.\n");
304 : }
305 :
306 7 : if (do_neigh) {
307 7 : log.printf(
308 : "maze> Using neighbor lists updated every %d steps and cutoff %f.\n",
309 : nl_stride_,
310 : nl_cutoff_
311 : );
312 : }
313 :
314 : // OpenMP
315 7 : stride_ = comm.Get_size();
316 7 : rank_ = comm.Get_rank();
317 :
318 7 : n_threads_ = OpenMP::getNumThreads();
319 7 : unsigned int nn = neighbor_list_->size();
320 :
321 7 : if (n_threads_ * stride_ * 10 > nn) {
322 0 : n_threads_ = nn / stride_ / 10;
323 : }
324 :
325 7 : if (n_threads_ == 0) {
326 0 : n_threads_ = 1;
327 : }
328 :
329 7 : if (keywords.exists("OPTIMIZER_STRIDE")) {
330 7 : parse("OPTIMIZER_STRIDE", optimizer_stride_);
331 :
332 7 : plumed_massert(
333 : optimizer_stride_,
334 : "maze> OPTIMIZER_STRIDE should be explicitly specified and positive.\n"
335 : );
336 :
337 7 : log.printf(
338 : "maze> Launching optimization every %u steps.\n",
339 : optimizer_stride_
340 : );
341 : }
342 :
343 7 : rnd::randomize();
344 :
345 7 : opt_.zero();
346 :
347 14 : addComponentWithDerivatives("x");
348 14 : componentIsNotPeriodic("x");
349 :
350 14 : addComponentWithDerivatives("y");
351 14 : componentIsNotPeriodic("y");
352 :
353 14 : addComponentWithDerivatives("z");
354 14 : componentIsNotPeriodic("z");
355 :
356 14 : addComponent("loss");
357 14 : componentIsNotPeriodic("loss");
358 :
359 14 : addComponent("sr");
360 7 : componentIsNotPeriodic("sr");
361 :
362 7 : value_x_ = getPntrToComponent("x");
363 7 : value_y_ = getPntrToComponent("y");
364 7 : value_z_ = getPntrToComponent("z");
365 7 : value_action_ = getPntrToComponent("loss");
366 7 : value_sampling_radius_ = getPntrToComponent("sr");
367 7 : }
368 :
369 16239210 : double Optimizer::pairing(double distance) const {
370 16239210 : return loss_->pairing(distance);
371 : }
372 :
373 6 : Vector Optimizer::center_of_mass() const {
374 6 : const unsigned nl_size = neighbor_list_->size();
375 :
376 6 : Vector center_of_mass;
377 6 : center_of_mass.zero();
378 : double mass = 0;
379 :
380 189654 : for (unsigned int i = 0; i < nl_size; ++i) {
381 189648 : unsigned int i0 = neighbor_list_->getClosePair(i).first;
382 189648 : center_of_mass += getPosition(i0) * getMass(i0);
383 189648 : mass += getMass(i0);
384 : }
385 :
386 6 : return center_of_mass / mass;
387 : }
388 :
389 210 : void Optimizer::prepare() {
390 210 : if (neighbor_list_->getStride() > 0) {
391 210 : if (first_time_ || (getStep() % neighbor_list_->getStride() == 0)) {
392 7 : requestAtoms(neighbor_list_->getFullAtomList());
393 :
394 7 : validate_list_ = true;
395 7 : first_time_ = false;
396 : } else {
397 203 : requestAtoms(neighbor_list_->getReducedAtomList());
398 :
399 203 : validate_list_ = false;
400 :
401 203 : if (getExchangeStep()) {
402 0 : plumed_merror(
403 : "maze> Neighbor lists should be updated on exchange steps -- choose "
404 : "an NL_STRIDE which divides the exchange stride.\n");
405 : }
406 : }
407 :
408 210 : if (getExchangeStep()) {
409 0 : first_time_ = true;
410 : }
411 : }
412 210 : }
413 :
414 226 : double Optimizer::score() {
415 226 : const unsigned nl_size = neighbor_list_->size();
416 226 : Vector distance;
417 : double function = 0;
418 :
419 226 : #pragma omp parallel num_threads(n_threads_)
420 : {
421 : #pragma omp for reduction(+:function)
422 : for(unsigned int i = 0; i < nl_size; i++) {
423 : unsigned i0 = neighbor_list_->getClosePair(i).first;
424 : unsigned i1 = neighbor_list_->getClosePair(i).second;
425 :
426 : if (getAbsoluteIndex(i0) == getAbsoluteIndex(i1)) {
427 : continue;
428 : }
429 :
430 : if (pbc_) {
431 : distance = pbcDistance(getPosition(i0), getPosition(i1));
432 : } else {
433 : distance = delta(getPosition(i0), getPosition(i1));
434 : }
435 :
436 : function += pairing(distance.modulo());
437 : }
438 : }
439 :
440 226 : return function;
441 : }
442 :
443 210 : void Optimizer::update_nl() {
444 210 : if (neighbor_list_->getStride() > 0 && validate_list_) {
445 7 : neighbor_list_->update(getPositions());
446 : }
447 210 : }
448 :
449 362 : double Optimizer::sampling_radius() {
450 362 : const unsigned nl_size=neighbor_list_->size();
451 362 : Vector d;
452 : double min=std::numeric_limits<int>::max();
453 :
454 9654278 : for (unsigned int i = 0; i < nl_size; ++i) {
455 9653916 : unsigned i0 = neighbor_list_->getClosePair(i).first;
456 9653916 : unsigned i1 = neighbor_list_->getClosePair(i).second;
457 :
458 9653916 : if (getAbsoluteIndex(i0) == getAbsoluteIndex(i1)) {
459 0 : continue;
460 : }
461 :
462 9653916 : if (pbc_) {
463 9653916 : d = pbcDistance(getPosition(i0), getPosition(i1));
464 : } else {
465 0 : d = delta(getPosition(i0), getPosition(i1));
466 : }
467 :
468 9653916 : double dist = d.modulo();
469 :
470 9653916 : if(dist < min) {
471 : min = dist;
472 : }
473 : }
474 :
475 362 : return min;
476 : }
477 :
478 210 : void Optimizer::calculate() {
479 210 : update_nl();
480 :
481 210 : if (getStep() % optimizer_stride_ == 0 && !first_step_) {
482 19 : optimize();
483 :
484 19 : value_x_->set(opt_[0]);
485 19 : value_y_->set(opt_[1]);
486 19 : value_z_->set(opt_[2]);
487 :
488 19 : value_action_->set(score());
489 19 : value_sampling_radius_->set(sampling_radius());
490 : } else {
491 191 : first_step_=false;
492 :
493 191 : value_x_->set(opt_[0]);
494 191 : value_y_->set(opt_[1]);
495 191 : value_z_->set(opt_[2]);
496 :
497 191 : value_action_->set(score());
498 191 : value_sampling_radius_->set(sampling_radius());
499 : }
500 210 : }
501 :
502 : } // namespace maze
503 : } // namespace PLMD
|
