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 :
32 : namespace PLMD {
33 : namespace maze {
34 :
35 12 : void Optimizer::registerKeywords(Keywords& keys) {
36 12 : Colvar::registerKeywords(keys);
37 :
38 24 : keys.addFlag(
39 : "SERIAL",
40 : false,
41 : "Perform the simulation in serial -- used only for debugging purposes, "
42 : "should not be used otherwise."
43 : );
44 :
45 24 : keys.addFlag(
46 : "PAIR",
47 : false,
48 : "Pair only the 1st element of the 1st group with the 1st element in the "
49 : "second, etc."
50 : );
51 :
52 24 : keys.addFlag(
53 : "NLIST",
54 : false,
55 : "Use a neighbor list of ligand-protein atom pairs to speed up the "
56 : "calculating of the distances."
57 : );
58 :
59 24 : keys.add(
60 : "optional",
61 : "NL_CUTOFF",
62 : "Neighbor list cut-off for the distances of ligand-protein atom pairs."
63 : );
64 :
65 24 : keys.add(
66 : "optional",
67 : "NL_STRIDE",
68 : "Update stride for the ligand-protein atom pairs in the neighbor list."
69 : );
70 :
71 24 : keys.add(
72 : "compulsory",
73 : "N_ITER",
74 : "Number of optimization steps. Required only for optimizers, do not pass "
75 : "this keyword to the fake optimizers (results in crash) , e.g., random "
76 : "walk, steered MD, or random acceleration MD."
77 : );
78 :
79 24 : keys.add(
80 : "optional",
81 : "LOSS",
82 : "Loss function describing ligand-protein interactions required by every "
83 : "optimizer."
84 : );
85 :
86 24 : keys.add(
87 : "atoms",
88 : "LIGAND",
89 : "Indices of ligand atoms."
90 : );
91 :
92 24 : keys.add(
93 : "atoms",
94 : "PROTEIN",
95 : "Indices of protein atoms."
96 : );
97 :
98 36 : keys.add(
99 : "compulsory",
100 : "OPTIMIZER_STRIDE",
101 : "Optimizer stride. Sets up a callback function that launches the "
102 : "optimization process every OPTIMIZER_STRIDE."
103 : );
104 :
105 12 : componentsAreNotOptional(keys);
106 :
107 24 : keys.addOutputComponent(
108 : "x",
109 : "default",
110 : "Optimal biasing direction; x component."
111 : );
112 :
113 24 : keys.addOutputComponent(
114 : "y",
115 : "default",
116 : "Optimal biasing direction; y component."
117 : );
118 :
119 24 : keys.addOutputComponent(
120 : "z",
121 : "default",
122 : "Optimal biasing direction; z component."
123 : );
124 :
125 24 : keys.addOutputComponent(
126 : "loss",
127 : "default",
128 : "Loss function value defined by the provided pairing function."
129 : );
130 :
131 24 : keys.addOutputComponent(
132 : "sr",
133 : "default",
134 : "Sampling radius. Reduces sampling to the local proximity of the ligand "
135 : "position."
136 : );
137 12 : }
138 :
139 7 : Optimizer::Optimizer(const ActionOptions& ao)
140 : : PLUMED_COLVAR_INIT(ao),
141 7 : first_step_(true),
142 7 : opt_value_(0.0),
143 7 : pbc_(true),
144 7 : sampling_r_(0.0),
145 7 : serial_(false),
146 7 : validate_list_(true),
147 7 : first_time_(true)
148 : {
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_ = new 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 7 : );
238 : }
239 : else {
240 0 : neighbor_list_=new NeighborList(
241 : ga_list,
242 : gb_list,
243 : serial_,
244 : do_pair,
245 0 : pbc_,
246 : getPbc(),
247 : comm
248 0 : );
249 : }
250 : }
251 : else {
252 0 : if (do_neigh) {
253 0 : neighbor_list_ = new NeighborList(
254 : ga_list,
255 : serial_,
256 0 : pbc_,
257 : getPbc(),
258 : comm,
259 0 : nl_cutoff_,
260 0 : nl_stride_
261 0 : );
262 : }
263 : else {
264 0 : neighbor_list_=new NeighborList(
265 : ga_list,
266 : serial_,
267 0 : pbc_,
268 : getPbc(),
269 : comm
270 0 : );
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 : }
299 : else {
300 0 : log.printf("maze> Without periodic boundary conditions.\n");
301 : }
302 :
303 7 : if (do_pair) {
304 0 : log.printf("maze> With PAIR option.\n");
305 : }
306 :
307 7 : if (do_neigh) {
308 7 : log.printf(
309 : "maze> Using neighbor lists updated every %d steps and cutoff %f.\n",
310 : nl_stride_,
311 : nl_cutoff_
312 : );
313 : }
314 :
315 : // OpenMP
316 7 : stride_ = comm.Get_size();
317 7 : rank_ = comm.Get_rank();
318 :
319 7 : n_threads_ = OpenMP::getNumThreads();
320 7 : unsigned int nn = neighbor_list_->size();
321 :
322 7 : if (n_threads_ * stride_ * 10 > nn) {
323 0 : n_threads_ = nn / stride_ / 10;
324 : }
325 :
326 7 : if (n_threads_ == 0) {
327 0 : n_threads_ = 1;
328 : }
329 :
330 14 : if (keywords.exists("OPTIMIZER_STRIDE")) {
331 7 : parse("OPTIMIZER_STRIDE", optimizer_stride_);
332 :
333 7 : plumed_massert(
334 : optimizer_stride_,
335 : "maze> OPTIMIZER_STRIDE should be explicitly specified and positive.\n"
336 : );
337 :
338 7 : log.printf(
339 : "maze> Launching optimization every %u steps.\n",
340 : optimizer_stride_
341 : );
342 : }
343 :
344 7 : rnd::randomize();
345 :
346 7 : opt_.zero();
347 :
348 7 : addComponentWithDerivatives("x");
349 7 : componentIsNotPeriodic("x");
350 :
351 7 : addComponentWithDerivatives("y");
352 7 : componentIsNotPeriodic("y");
353 :
354 7 : addComponentWithDerivatives("z");
355 7 : componentIsNotPeriodic("z");
356 :
357 7 : addComponent("loss");
358 7 : componentIsNotPeriodic("loss");
359 :
360 7 : addComponent("sr");
361 7 : componentIsNotPeriodic("sr");
362 :
363 7 : value_x_ = getPntrToComponent("x");
364 7 : value_y_ = getPntrToComponent("y");
365 7 : value_z_ = getPntrToComponent("z");
366 7 : value_action_ = getPntrToComponent("loss");
367 7 : value_sampling_radius_ = getPntrToComponent("sr");
368 7 : }
369 :
370 15607050 : double Optimizer::pairing(double distance) const {
371 15607050 : return loss_->pairing(distance);
372 : }
373 :
374 6 : Vector Optimizer::center_of_mass() const {
375 6 : const unsigned nl_size = neighbor_list_->size();
376 :
377 6 : Vector center_of_mass;
378 6 : center_of_mass.zero();
379 : double mass = 0;
380 :
381 189654 : for (unsigned int i = 0; i < nl_size; ++i) {
382 189648 : unsigned int i0 = neighbor_list_->getClosePair(i).first;
383 189648 : center_of_mass += getPosition(i0) * getMass(i0);
384 189648 : mass += getMass(i0);
385 : }
386 :
387 6 : return center_of_mass / mass;
388 : }
389 :
390 210 : void Optimizer::prepare() {
391 210 : if (neighbor_list_->getStride() > 0) {
392 210 : if (first_time_ || (getStep() % neighbor_list_->getStride() == 0)) {
393 7 : requestAtoms(neighbor_list_->getFullAtomList());
394 :
395 7 : validate_list_ = true;
396 7 : first_time_ = false;
397 : }
398 : else {
399 203 : requestAtoms(neighbor_list_->getReducedAtomList());
400 :
401 203 : validate_list_ = false;
402 :
403 203 : if (getExchangeStep()) {
404 0 : plumed_merror(
405 : "maze> Neighbor lists should be updated on exchange steps -- choose "
406 : "an NL_STRIDE which divides the exchange stride.\n");
407 : }
408 : }
409 :
410 210 : if (getExchangeStep()) {
411 0 : first_time_ = true;
412 : }
413 : }
414 210 : }
415 :
416 226 : double Optimizer::score() {
417 226 : const unsigned nl_size = neighbor_list_->size();
418 226 : Vector distance;
419 : double function = 0;
420 :
421 226 : #pragma omp parallel num_threads(n_threads_)
422 : {
423 : #pragma omp for reduction(+:function)
424 : for(unsigned int i = 0; i < nl_size; i++) {
425 : unsigned i0 = neighbor_list_->getClosePair(i).first;
426 : unsigned i1 = neighbor_list_->getClosePair(i).second;
427 :
428 : if (getAbsoluteIndex(i0) == getAbsoluteIndex(i1)) {
429 : continue;
430 : }
431 :
432 : if (pbc_) {
433 : distance = pbcDistance(getPosition(i0), getPosition(i1));
434 : }
435 : else {
436 : distance = delta(getPosition(i0), getPosition(i1));
437 : }
438 :
439 : function += pairing(distance.modulo());
440 : }
441 : }
442 :
443 226 : return function;
444 : }
445 :
446 210 : void Optimizer::update_nl() {
447 210 : if (neighbor_list_->getStride() > 0 && validate_list_) {
448 7 : neighbor_list_->update(getPositions());
449 : }
450 210 : }
451 :
452 351 : double Optimizer::sampling_radius()
453 : {
454 351 : const unsigned nl_size=neighbor_list_->size();
455 351 : Vector d;
456 : double min=std::numeric_limits<int>::max();
457 :
458 9306579 : for (unsigned int i = 0; i < nl_size; ++i) {
459 9306228 : unsigned i0 = neighbor_list_->getClosePair(i).first;
460 9306228 : unsigned i1 = neighbor_list_->getClosePair(i).second;
461 :
462 9306228 : if (getAbsoluteIndex(i0) == getAbsoluteIndex(i1)) {
463 0 : continue;
464 : }
465 :
466 9306228 : if (pbc_) {
467 9306228 : d = pbcDistance(getPosition(i0), getPosition(i1));
468 : }
469 : else {
470 0 : d = delta(getPosition(i0), getPosition(i1));
471 : }
472 :
473 9306228 : double dist = d.modulo();
474 :
475 9306228 : if(dist < min) {
476 : min = dist;
477 : }
478 : }
479 :
480 351 : return min;
481 : }
482 :
483 210 : void Optimizer::calculate() {
484 210 : update_nl();
485 :
486 210 : if (getStep() % optimizer_stride_ == 0 && !first_step_) {
487 19 : optimize();
488 :
489 19 : value_x_->set(opt_[0]);
490 19 : value_y_->set(opt_[1]);
491 19 : value_z_->set(opt_[2]);
492 :
493 19 : value_action_->set(score());
494 19 : value_sampling_radius_->set(sampling_radius());
495 : }
496 : else {
497 191 : first_step_=false;
498 :
499 191 : value_x_->set(opt_[0]);
500 191 : value_y_->set(opt_[1]);
501 191 : value_z_->set(opt_[2]);
502 :
503 191 : value_action_->set(score());
504 191 : value_sampling_radius_->set(sampling_radius());
505 : }
506 210 : }
507 :
508 : } // namespace maze
509 : } // namespace PLMD
|
