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