ESPResSo
Extensible Simulation Package for Research on Soft Matter Systems
Loading...
Searching...
No Matches
integrate.cpp
Go to the documentation of this file.
1/*
2 * Copyright (C) 2010-2023 The ESPResSo project
3 * Copyright (C) 2002,2003,2004,2005,2006,2007,2008,2009,2010
4 * Max-Planck-Institute for Polymer Research, Theory Group
5 *
6 * This file is part of ESPResSo.
7 *
8 * ESPResSo is free software: you can redistribute it and/or modify
9 * it under the terms of the GNU General Public License as published by
10 * the Free Software Foundation, either version 3 of the License, or
11 * (at your option) any later version.
12 *
13 * ESPResSo is distributed in the hope that it will be useful,
14 * but WITHOUT ANY WARRANTY; without even the implied warranty of
15 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
16 * GNU General Public License for more details.
17 *
18 * You should have received a copy of the GNU General Public License
19 * along with this program. If not, see <http://www.gnu.org/licenses/>.
20 */
21
22/** \file
23 * Molecular dynamics integrator.
24 *
25 * For more information about the integrator
26 * see \ref integrate.hpp "integrate.hpp".
27 */
28
29#include "integrate.hpp"
36
37#include "BoxGeometry.hpp"
38#include "ParticleRange.hpp"
39#include "PropagationMode.hpp"
40#include "accumulators/AutoUpdateAccumulators.hpp"
44#include "cells.hpp"
45#include "collision_detection/CollisionDetection.hpp"
46#include "communication.hpp"
47#include "errorhandling.hpp"
48#include "forces.hpp"
50#include "lb/utils.hpp"
53#include "npt.hpp"
54#include "rattle.hpp"
55#include "rotation.hpp"
56#include "signalhandling.hpp"
57#include "system/System.hpp"
58#include "thermostat.hpp"
62
63#include <boost/mpi/collectives/all_reduce.hpp>
64
65#ifdef CALIPER
66#include <caliper/cali.h>
67#endif
68
69#ifdef VALGRIND
70#include <callgrind.h>
71#endif
72
73#include <algorithm>
74#include <cassert>
75#include <cmath>
76#include <csignal>
77#include <functional>
78#include <limits>
79#include <stdexcept>
80#include <string>
81#include <utility>
82
83#ifdef WALBERLA
84#ifdef WALBERLA_STATIC_ASSERT
85#error "waLberla headers should not be visible to the ESPResSo core"
86#endif
87#endif
88
89namespace {
90volatile std::sig_atomic_t ctrl_C = 0;
91} // namespace
92
93namespace LeesEdwards {
94
95/**
96 * @brief Update the Lees-Edwards parameters of the box geometry
97 * for the current simulation time.
98 */
99void LeesEdwards::update_box_params(BoxGeometry &box_geo, double sim_time) {
100 if (box_geo.type() == BoxType::LEES_EDWARDS) {
101 assert(m_protocol != nullptr);
102 box_geo.lees_edwards_update(get_pos_offset(sim_time, *m_protocol),
103 get_shear_velocity(sim_time, *m_protocol));
104 }
105}
106
107void LeesEdwards::set_protocol(std::shared_ptr<ActiveProtocol> protocol) {
108 auto &system = get_system();
109 auto &cell_structure = *system.cell_structure;
110 auto &box_geo = *system.box_geo;
111 box_geo.set_type(BoxType::LEES_EDWARDS);
112 m_protocol = std::move(protocol);
113 update_box_params(box_geo, system.get_sim_time());
114 system.propagation->recalc_forces = true;
115 cell_structure.set_resort_particles(Cells::RESORT_LOCAL);
116}
117
119 auto &system = get_system();
120 auto &cell_structure = *system.cell_structure;
121 auto &box_geo = *system.box_geo;
122 m_protocol = nullptr;
123 box_geo.set_type(BoxType::CUBOID);
124 system.propagation->recalc_forces = true;
125 cell_structure.set_resort_particles(Cells::RESORT_LOCAL);
126}
127
128} // namespace LeesEdwards
129
131 switch (integ_switch) {
134 break;
135 case INTEG_METHOD_NVT: {
136 // NOLINTNEXTLINE(bugprone-branch-clone)
137 if ((thermo_switch & THERMO_LB) and (thermo_switch & THERMO_LANGEVIN)) {
139#ifdef ROTATION
141#endif
142 } else if (thermo_switch & THERMO_LB) {
144#ifdef ROTATION
146#endif
147 } else if (thermo_switch & THERMO_LANGEVIN) {
149#ifdef ROTATION
151#endif
152 } else {
154#ifdef ROTATION
156#endif
157 }
158 break;
159 }
160#ifdef NPT
163 break;
164#endif
165 case INTEG_METHOD_BD:
167#ifdef ROTATION
169#endif
170 break;
171#ifdef STOKESIAN_DYNAMICS
172 case INTEG_METHOD_SD:
174 break;
175#endif // STOKESIAN_DYNAMICS
176 default:
177 throw std::runtime_error("Unknown value for integ_switch");
178 }
179}
180
182 int used_propagations = PropagationMode::NONE;
183 for (auto &p : cell_structure->local_particles()) {
184 used_propagations |= p.propagation();
185 }
186 if (used_propagations & PropagationMode::SYSTEM_DEFAULT) {
187 used_propagations |= propagation->default_propagation;
188 }
189 used_propagations = boost::mpi::all_reduce(::comm_cart, used_propagations,
190 std::bit_or<int>());
191 propagation->used_propagations = used_propagations;
192}
193
194void System::System::integrator_sanity_checks() const {
195 auto const thermo_switch = thermostat->thermo_switch;
196 if (time_step <= 0.) {
197 runtimeErrorMsg() << "time_step not set";
198 }
199 if (propagation->integ_switch == INTEG_METHOD_STEEPEST_DESCENT) {
200 if (thermo_switch != THERMO_OFF) {
202 << "The steepest descent integrator is incompatible with thermostats";
203 }
204 }
205 if (propagation->integ_switch == INTEG_METHOD_NVT) {
206 if (thermo_switch & (THERMO_NPT_ISO | THERMO_BROWNIAN | THERMO_SD)) {
207 runtimeErrorMsg() << "The VV integrator is incompatible with the "
208 "currently active combination of thermostats";
209 }
210 }
211#ifdef NPT
212 if (propagation->used_propagations & PropagationMode::TRANS_LANGEVIN_NPT) {
213 if (thermo_switch != THERMO_NPT_ISO) {
214 runtimeErrorMsg() << "The NpT integrator requires the NpT thermostat";
215 }
216 if (box_geo->type() == BoxType::LEES_EDWARDS) {
217 runtimeErrorMsg() << "The NpT integrator cannot use Lees-Edwards";
218 }
219 }
220#endif
221 if (propagation->used_propagations & PropagationMode::TRANS_BROWNIAN) {
222 if (thermo_switch != THERMO_BROWNIAN) {
223 runtimeErrorMsg() << "The BD integrator requires the BD thermostat";
224 }
225 }
226 if (propagation->used_propagations & PropagationMode::TRANS_STOKESIAN) {
227#ifdef STOKESIAN_DYNAMICS
228 if (thermo_switch != THERMO_SD) {
229 runtimeErrorMsg() << "The SD integrator requires the SD thermostat";
230 }
231#endif
232 }
233 if (lb.is_solver_set() and (propagation->used_propagations &
236 if (thermostat->lb == nullptr) {
237 runtimeErrorMsg() << "The LB integrator requires the LB thermostat";
238 }
239 }
240 if (bonded_ias->get_n_thermalized_bonds() >= 1 and
241 (thermostat->thermalized_bond == nullptr or
242 (thermo_switch & THERMO_BOND) == 0)) {
244 << "Thermalized bonds require the thermalized_bond thermostat";
245 }
246
247#ifdef ROTATION
248 for (auto const &p : cell_structure->local_particles()) {
249 using namespace PropagationMode;
250 if (p.can_rotate() and not p.is_virtual() and
251 (p.propagation() & (SYSTEM_DEFAULT | ROT_EULER | ROT_LANGEVIN |
252 ROT_BROWNIAN | ROT_STOKESIAN)) == 0) {
254 << "Rotating particles must have a rotation propagation mode enabled";
255 break;
256 }
257 }
258#endif
259}
260
261#ifdef WALBERLA
262void walberla_tau_sanity_checks(std::string method, double tau,
263 double time_step) {
264 if (time_step <= 0.) {
265 return;
266 }
267 // use float epsilon since tau may be a float
268 auto const eps = static_cast<double>(std::numeric_limits<float>::epsilon());
269 if ((tau - time_step) / (tau + time_step) < -eps)
270 throw std::invalid_argument(method + " tau (" + std::to_string(tau) +
271 ") must be >= MD time_step (" +
272 std::to_string(time_step) + ")");
273 auto const factor = tau / time_step;
274 if (std::fabs(std::round(factor) - factor) / factor > eps)
275 throw std::invalid_argument(method + " tau (" + std::to_string(tau) +
276 ") must be an integer multiple of the "
277 "MD time_step (" +
278 std::to_string(time_step) + "). Factor is " +
279 std::to_string(factor));
280}
281
282void walberla_agrid_sanity_checks(std::string method,
283 Utils::Vector3d const &geo_left,
284 Utils::Vector3d const &geo_right,
285 Utils::Vector3d const &lattice_left,
286 Utils::Vector3d const &lattice_right,
287 double agrid) {
288 // waLBerla and ESPResSo must agree on domain decomposition
289 auto const tol = agrid / 1E6;
290 if ((lattice_left - geo_left).norm2() > tol or
291 (lattice_right - geo_right).norm2() > tol) {
292 runtimeErrorMsg() << "\nMPI rank " << ::this_node << ": "
293 << "left ESPResSo: [" << geo_left << "], "
294 << "left waLBerla: [" << lattice_left << "]"
295 << "\nMPI rank " << ::this_node << ": "
296 << "right ESPResSo: [" << geo_right << "], "
297 << "right waLBerla: [" << lattice_right << "]"
298 << "\nfor method: " << method;
299 throw std::runtime_error(
300 "waLBerla and ESPResSo disagree about domain decomposition.");
301 }
302}
303#endif // WALBERLA
304
306 auto &cell_structure = *system.cell_structure;
307 auto const offset = LeesEdwards::verlet_list_offset(
308 *system.box_geo, cell_structure.get_le_pos_offset_at_last_resort());
309 if (cell_structure.check_resort_required(offset)) {
310 cell_structure.set_resort_particles(Cells::RESORT_LOCAL);
311 }
312}
313
315 auto const &propagation = *this->propagation;
316 if ((not thermostat->langevin) or ((propagation.used_propagations &
319 return;
320 }
321 auto const &langevin = *thermostat->langevin;
322 auto const kT = thermostat->kT;
323 for (auto &p : cell_structure->local_particles()) {
324 if (propagation.should_propagate_with(p, PropagationMode::TRANS_LANGEVIN))
325 p.force() += friction_thermo_langevin(langevin, p, time_step, kT);
326#ifdef ROTATION
327 if (propagation.should_propagate_with(p, PropagationMode::ROT_LANGEVIN))
328 p.torque() += convert_vector_body_to_space(
329 p, friction_thermo_langevin_rotation(langevin, p, time_step, kT));
330#endif
331 }
332}
333
334/** @brief Calls the hook for propagation kernels before the force calculation
335 * @return whether or not to stop the integration loop early.
336 */
337static bool integrator_step_1(ParticleRange const &particles,
338 Propagation const &propagation,
339 System::System &system, double time_step) {
340 // steepest decent
342 return steepest_descent_step(particles);
343
344 auto const &thermostat = *system.thermostat;
345 auto const kT = thermostat.kT;
346 for (auto &p : particles) {
347#ifdef VIRTUAL_SITES
348 // virtual sites are updated later in the integration loop
349 if (p.is_virtual())
350 continue;
351#endif
352 if (propagation.should_propagate_with(
354 velocity_verlet_propagator_1(p, time_step);
356 velocity_verlet_propagator_1(p, time_step);
357#ifdef ROTATION
359 velocity_verlet_rotator_1(p, time_step);
360#endif
362 velocity_verlet_propagator_1(p, time_step);
363#ifdef ROTATION
365 velocity_verlet_rotator_1(p, time_step);
366#endif
368 brownian_dynamics_propagator(*thermostat.brownian, p, time_step, kT);
369#ifdef ROTATION
371 brownian_dynamics_rotator(*thermostat.brownian, p, time_step, kT);
372#endif
373 }
374
375#ifdef NPT
378 auto pred = PropagationPredicateNPT(propagation.default_propagation);
379 velocity_verlet_npt_step_1(particles.filter(pred), *thermostat.npt_iso,
380 time_step, system);
381 }
382#endif
383
384#ifdef STOKESIAN_DYNAMICS
387 auto pred = PropagationPredicateStokesian(propagation.default_propagation);
388 stokesian_dynamics_step_1(particles.filter(pred), *thermostat.stokesian,
389 time_step, kT);
390 }
391#endif // STOKESIAN_DYNAMICS
392
393 return false;
394}
395
396static void integrator_step_2(ParticleRange const &particles,
397 Propagation const &propagation,
398 Thermostat::Thermostat const &thermostat,
399 double time_step) {
401 return;
402
403 for (auto &p : particles) {
404#ifdef VIRTUAL_SITES
405 // virtual sites are updated later in the integration loop
406 if (p.is_virtual())
407 continue;
408#endif
409 if (propagation.should_propagate_with(
411 velocity_verlet_propagator_2(p, time_step);
413 velocity_verlet_propagator_2(p, time_step);
414#ifdef ROTATION
416 velocity_verlet_rotator_2(p, time_step);
417#endif
419 velocity_verlet_propagator_2(p, time_step);
420#ifdef ROTATION
422 velocity_verlet_rotator_2(p, time_step);
423#endif
424 }
425
426#ifdef NPT
429 auto pred = PropagationPredicateNPT(propagation.default_propagation);
430 velocity_verlet_npt_step_2(particles.filter(pred), *thermostat.npt_iso,
431 time_step);
432 }
433#endif
434}
435
436int System::System::integrate(int n_steps, int reuse_forces) {
437#ifdef CALIPER
438 CALI_CXX_MARK_FUNCTION;
439#endif
440 auto &propagation = *this->propagation;
441#ifdef VIRTUAL_SITES_RELATIVE
442 auto const has_vs_rel = [&propagation]() {
443 return propagation.used_propagations & (PropagationMode::ROT_VS_RELATIVE |
445 };
446#endif
447#ifdef BOND_CONSTRAINT
448 auto const n_rigid_bonds = bonded_ias->get_n_rigid_bonds();
449#endif
450
451 // Prepare particle structure and run sanity checks of all active algorithms
452 propagation.update_default_propagation(thermostat->thermo_switch);
453 update_used_propagations();
454 on_integration_start();
455
456 // If any method vetoes (e.g. P3M not initialized), immediately bail out
458 return INTEG_ERROR_RUNTIME;
459
460 // Additional preparations for the first integration step
461 if (reuse_forces == INTEG_REUSE_FORCES_NEVER or
462 ((reuse_forces != INTEG_REUSE_FORCES_ALWAYS) and
463 propagation.recalc_forces)) {
464#ifdef CALIPER
465 CALI_MARK_BEGIN("Initial Force Calculation");
466#endif
467 thermostat->lb_coupling_deactivate();
468
469#ifdef VIRTUAL_SITES_RELATIVE
470 if (has_vs_rel()) {
471 vs_relative_update_particles(*cell_structure, *box_geo);
472 }
473#endif
474
475 // Communication step: distribute ghost positions
476 cell_structure->update_ghosts_and_resort_particle(get_global_ghost_flags());
477
478 calculate_forces();
479
480 if (propagation.integ_switch != INTEG_METHOD_STEEPEST_DESCENT) {
481#ifdef ROTATION
482 convert_initial_torques(cell_structure->local_particles());
483#endif
484 }
485
486#ifdef CALIPER
487 CALI_MARK_END("Initial Force Calculation");
488#endif
489 }
490
491 thermostat->lb_coupling_activate();
492
494 return INTEG_ERROR_RUNTIME;
495
496 // Keep track of the number of Verlet updates (i.e. particle resorts)
497 int n_verlet_updates = 0;
498
499 // Keep track of whether an interrupt signal was caught (only in singleton
500 // mode, since signal handlers are unreliable with more than 1 MPI rank)
501 auto const singleton_mode = comm_cart.size() == 1;
502 auto caught_sigint = false;
503 auto caught_error = false;
504
505 auto lb_active = false;
506 auto ek_active = false;
507 if (propagation.integ_switch != INTEG_METHOD_STEEPEST_DESCENT) {
508 lb_active = lb.is_solver_set();
509 ek_active = ek.is_ready_for_propagation();
510 }
511 auto const calc_md_steps_per_tau = [this](double tau) {
512 return static_cast<int>(std::round(tau / time_step));
513 };
514
515#ifdef VALGRIND
516 CALLGRIND_START_INSTRUMENTATION;
517#endif
518 // Integration loop
519#ifdef CALIPER
520 CALI_CXX_MARK_LOOP_BEGIN(integration_loop, "Integration loop");
521#endif
522 int integrated_steps = 0;
523 for (int step = 0; step < n_steps; step++) {
524#ifdef CALIPER
525 CALI_CXX_MARK_LOOP_ITERATION(integration_loop, step);
526#endif
527
528 auto particles = cell_structure->local_particles();
529
530#ifdef BOND_CONSTRAINT
531 if (n_rigid_bonds)
532 save_old_position(particles, cell_structure->ghost_particles());
533#endif
534
535 lees_edwards->update_box_params(*box_geo, sim_time);
536 bool early_exit =
537 integrator_step_1(particles, propagation, *this, time_step);
538 if (early_exit)
539 break;
540
541 sim_time += time_step;
542 if (box_geo->type() == BoxType::LEES_EDWARDS) {
543 auto const kernel = LeesEdwards::Push{*box_geo};
544 for (auto &p : particles) {
545 kernel(p);
546 }
547 }
548
549#ifdef NPT
550 if (propagation.integ_switch != INTEG_METHOD_NPT_ISO)
551#endif
552 {
554 }
555
556 // Propagate philox RNG counters
557 thermostat->philox_counter_increment();
558
559#ifdef BOND_CONSTRAINT
560 // Correct particle positions that participate in a rigid/constrained bond
561 if (n_rigid_bonds) {
562 correct_position_shake(*cell_structure, *box_geo, *bonded_ias);
563 }
564#endif
565
566#ifdef VIRTUAL_SITES_RELATIVE
567 if (has_vs_rel()) {
568#ifdef NPT
569 if (propagation.integ_switch == INTEG_METHOD_NPT_ISO) {
570 cell_structure->update_ghosts_and_resort_particle(
572 }
573#endif // NPT
574 vs_relative_update_particles(*cell_structure, *box_geo);
575 }
576#endif // VIRTUAL_SITES_RELATIVE
577
578 if (cell_structure->get_resort_particles() >= Cells::RESORT_LOCAL)
579 n_verlet_updates++;
580
581 // Communication step: distribute ghost positions
582 cell_structure->update_ghosts_and_resort_particle(get_global_ghost_flags());
583
584 particles = cell_structure->local_particles();
585
586 calculate_forces();
587
588#ifdef VIRTUAL_SITES_INERTIALESS_TRACERS
589 if (thermostat->lb and
590 (propagation.used_propagations & PropagationMode::TRANS_LB_TRACER)) {
591 lb_tracers_add_particle_force_to_fluid(*cell_structure, *box_geo,
592 *local_geo, lb);
593 }
594#endif
595 integrator_step_2(particles, propagation, *thermostat, time_step);
596 if (propagation.integ_switch == INTEG_METHOD_BD) {
598 }
599 if (box_geo->type() == BoxType::LEES_EDWARDS) {
600 auto const kernel = LeesEdwards::UpdateOffset{*box_geo};
601 for (auto &p : particles) {
602 kernel(p);
603 }
604 }
605#ifdef BOND_CONSTRAINT
606 if (n_rigid_bonds) {
607 correct_velocity_shake(*cell_structure, *box_geo, *bonded_ias);
608 }
609#endif
610
611 // propagate one-step functionalities
612 if (propagation.integ_switch != INTEG_METHOD_STEEPEST_DESCENT) {
613 if (lb_active and ek_active) {
614 // assume that they are coupled, which is not necessarily true
615 auto const md_steps_per_lb_step = calc_md_steps_per_tau(lb.get_tau());
616 auto const md_steps_per_ek_step = calc_md_steps_per_tau(ek.get_tau());
617
618 if (md_steps_per_lb_step != md_steps_per_ek_step) {
620 << "LB and EK are active but with different time steps.";
621 }
622
623 assert(propagation.lb_skipped_md_steps ==
624 propagation.ek_skipped_md_steps);
625
626 propagation.lb_skipped_md_steps += 1;
627 propagation.ek_skipped_md_steps += 1;
628 if (propagation.lb_skipped_md_steps >= md_steps_per_lb_step) {
629 propagation.lb_skipped_md_steps = 0;
630 propagation.ek_skipped_md_steps = 0;
631 lb.propagate();
632 lb.ghost_communication_vel();
633 ek.propagate();
634 }
635 } else if (lb_active) {
636 auto const md_steps_per_lb_step = calc_md_steps_per_tau(lb.get_tau());
637 propagation.lb_skipped_md_steps += 1;
638 if (propagation.lb_skipped_md_steps >= md_steps_per_lb_step) {
639 propagation.lb_skipped_md_steps = 0;
640 lb.propagate();
641 }
642 } else if (ek_active) {
643 auto const md_steps_per_ek_step = calc_md_steps_per_tau(ek.get_tau());
644 propagation.ek_skipped_md_steps += 1;
645 if (propagation.ek_skipped_md_steps >= md_steps_per_ek_step) {
646 propagation.ek_skipped_md_steps = 0;
647 ek.propagate();
648 }
649 }
650 if (lb_active and (propagation.used_propagations &
652 thermostat->lb->rng_increment();
653 }
654
655#ifdef VIRTUAL_SITES_INERTIALESS_TRACERS
656 if (thermostat->lb and
657 (propagation.used_propagations & PropagationMode::TRANS_LB_TRACER)) {
658 if (lb_active) {
659 lb.ghost_communication_vel();
660 }
661 lb_tracers_propagate(*cell_structure, lb, time_step);
662 }
663#endif
664
665#ifdef COLLISION_DETECTION
666 collision_detection->handle_collisions();
667#endif
668 bond_breakage->process_queue(*this);
669 }
670
671 integrated_steps++;
672
674 caught_error = true;
675 break;
676 }
677
678 // Check if SIGINT has been caught.
679 if (singleton_mode and ctrl_C == 1) {
680 caught_sigint = true;
681 break;
682 }
683
684 } // for-loop over integration steps
685 if (lb_active) {
686 lb.ghost_communication();
687 }
688 lees_edwards->update_box_params(*box_geo, sim_time);
689#ifdef CALIPER
690 CALI_CXX_MARK_LOOP_END(integration_loop);
691#endif
692
693#ifdef VALGRIND
694 CALLGRIND_STOP_INSTRUMENTATION;
695#endif
696
697#ifdef VIRTUAL_SITES_RELATIVE
698 if (has_vs_rel()) {
699 vs_relative_update_particles(*cell_structure, *box_geo);
700 }
701#endif
702
703 // Verlet list statistics
704 cell_structure->update_verlet_stats(n_steps, n_verlet_updates);
705
706#ifdef NPT
707 if (propagation.integ_switch == INTEG_METHOD_NPT_ISO) {
709 }
710#endif
711 if (caught_sigint) {
712 ctrl_C = 0;
713 return INTEG_ERROR_SIGINT;
714 }
715 if (caught_error) {
716 return INTEG_ERROR_RUNTIME;
717 }
718 return integrated_steps;
719}
720
721int System::System::integrate_with_signal_handler(int n_steps, int reuse_forces,
722 bool update_accumulators) {
723 assert(n_steps >= 0);
724
725 // Override the signal handler so that the integrator obeys Ctrl+C
726 SignalHandler sa(SIGINT, [](int) { ctrl_C = 1; });
727
728 /* if skin wasn't set, do an educated guess now */
729 if (not cell_structure->is_verlet_skin_set()) {
730 try {
731 cell_structure->set_verlet_skin_heuristic();
732 } catch (...) {
733 if (comm_cart.rank() == 0) {
734 throw;
735 }
736 return INTEG_ERROR_RUNTIME;
737 }
738 }
739
740 if (not update_accumulators or n_steps == 0) {
741 return integrate(n_steps, reuse_forces);
742 }
743
744 for (int i = 0; i < n_steps;) {
745 /* Integrate to either the next accumulator update, or the
746 * end, depending on what comes first. */
747 auto const steps =
748 std::min((n_steps - i), auto_update_accumulators->next_update());
749
750 auto const local_retval = integrate(steps, reuse_forces);
751
752 // make sure all ranks exit when one rank fails
753 std::remove_const_t<decltype(local_retval)> global_retval;
754 boost::mpi::all_reduce(comm_cart, local_retval, global_retval,
755 std::plus<int>());
756 if (global_retval < 0) {
757 return global_retval; // propagate error code
758 }
759
760 reuse_forces = INTEG_REUSE_FORCES_ALWAYS;
761
762 (*auto_update_accumulators)(comm_cart, steps);
763
764 i += steps;
765 }
766
767 return 0;
768}
769
771 sim_time = value;
772 propagation->recalc_forces = true;
773 lees_edwards->update_box_params(*box_geo, sim_time);
774}
@ LEES_EDWARDS
@ INTEG_METHOD_SD
@ INTEG_METHOD_NPT_ISO
@ INTEG_METHOD_STEEPEST_DESCENT
@ INTEG_METHOD_NVT
@ INTEG_METHOD_BD
@ THERMO_SD
@ THERMO_BROWNIAN
@ THERMO_BOND
@ THERMO_LB
@ THERMO_LANGEVIN
@ THERMO_NPT_ISO
@ THERMO_OFF
Data structures for bonded interactions.
This file contains everything related to the global cell structure / cell system.
void lees_edwards_update(double pos_offset, double shear_velocity)
Update the Lees-Edwards parameters of the box geometry for the current simulation time.
BoxType type() const
void update_box_params(BoxGeometry &box_geo, double sim_time)
Update the Lees-Edwards parameters of the box geometry for the current simulation time.
Definition integrate.cpp:99
void set_protocol(std::shared_ptr< ActiveProtocol > protocol)
Set a new Lees-Edwards protocol.
void unset_protocol()
Delete the currently active Lees-Edwards protocol.
A range of particles.
ParticleRangeFiltered< Predicate > filter(Predicate pred) const
int default_propagation
void update_default_propagation(int thermo_switch)
bool should_propagate_with(Particle const &p, int mode) const
int used_propagations
RAII guard for signal handling.
Main system class.
void thermostat_force_init()
Calculate initial particle forces from active thermostats.
void update_used_propagations()
Update the global propagation bitmask.
void set_sim_time(double value)
Set sim_time.
int integrate_with_signal_handler(int n_steps, int reuse_forces, bool update_accumulators)
int integrate(int n_steps, int reuse_forces)
Integrate equations of motion.
std::shared_ptr< Thermostat::Thermostat > thermostat
std::shared_ptr< CellStructure > cell_structure
std::shared_ptr< BoxGeometry > box_geo
std::shared_ptr< IsotropicNptThermostat > npt_iso
boost::mpi::communicator comm_cart
The communicator.
int this_node
The number of this node.
int check_runtime_errors(boost::mpi::communicator const &comm)
Count runtime errors on all nodes.
This file contains the errorhandling code for severe errors, like a broken bond or illegal parameter ...
#define runtimeErrorMsg()
Force calculation.
static void resort_particles_if_needed(System::System &system)
static bool integrator_step_1(ParticleRange const &particles, Propagation const &propagation, System::System &system, double time_step)
Calls the hook for propagation kernels before the force calculation.
static void integrator_step_2(ParticleRange const &particles, Propagation const &propagation, Thermostat::Thermostat const &thermostat, double time_step)
void walberla_agrid_sanity_checks(std::string method, Utils::Vector3d const &geo_left, Utils::Vector3d const &geo_right, Utils::Vector3d const &lattice_left, Utils::Vector3d const &lattice_right, double agrid)
void walberla_tau_sanity_checks(std::string method, double tau, double time_step)
Molecular dynamics integrator.
#define INTEG_ERROR_RUNTIME
Definition integrate.hpp:42
#define INTEG_ERROR_SIGINT
Definition integrate.hpp:43
#define INTEG_REUSE_FORCES_NEVER
recalculate forces unconditionally (mostly used for timing)
Definition integrate.hpp:49
#define INTEG_REUSE_FORCES_ALWAYS
do not recalculate forces (mostly when reading checkpoints with forces)
Definition integrate.hpp:53
void brownian_dynamics_rotator(BrownianThermostat const &brownian, Particle &p, double time_step, double kT)
void brownian_dynamics_propagator(BrownianThermostat const &brownian, Particle &p, double time_step, double kT)
Utils::Vector3d friction_thermo_langevin(LangevinThermostat const &langevin, Particle const &p, double time_step, double kT)
Langevin thermostat for particle translational velocities.
Utils::Vector3d friction_thermo_langevin_rotation(LangevinThermostat const &langevin, Particle const &p, double time_step, double kT)
Langevin thermostat for particle angular velocities.
void lb_tracers_propagate(CellStructure &cell_structure, LB::Solver const &lb, double time_step)
void lb_tracers_add_particle_force_to_fluid(CellStructure &cell_structure, BoxGeometry const &box_geo, LocalBox const &local_box, LB::Solver &lb)
@ DATA_PART_PROPERTIES
Particle::p.
Utils::Vector3d verlet_list_offset(BoxGeometry const &box, double pos_offset_at_last_resort)
double get_shear_velocity(double time, ActiveProtocol const &protocol)
Calculation of current velocity.
double get_pos_offset(double time, ActiveProtocol const &protocol)
Definition protocols.hpp:88
volatile std::sig_atomic_t ctrl_C
Definition integrate.cpp:90
Various procedures concerning interactions between particles.
void synchronize_npt_state()
Synchronizes NpT state such as instantaneous and average pressure.
Definition npt.cpp:44
Exports for the NpT code.
void correct_velocity_shake(CellStructure &cs, BoxGeometry const &box_geo, BondedInteractionsMap const &bonded_ias)
Correction of current velocities using RATTLE algorithm.
Definition rattle.cpp:224
void save_old_position(const ParticleRange &particles, const ParticleRange &ghost_particles)
copy current position
Definition rattle.cpp:48
void correct_position_shake(CellStructure &cs, BoxGeometry const &box_geo, BondedInteractionsMap const &bonded_ias)
Propagate velocity and position while using SHAKE algorithm for bond constraint.
Definition rattle.cpp:146
RATTLE algorithm ().
void vs_relative_update_particles(CellStructure &cell_structure, BoxGeometry const &box_geo)
Definition relative.cpp:122
void convert_initial_torques(const ParticleRange &particles)
Convert torques to the body-fixed frame before the integration loop.
Definition rotation.cpp:188
This file contains all subroutines required to process rotational motion.
Utils::Vector3d convert_vector_body_to_space(const Particle &p, const Utils::Vector3d &vec)
Definition rotation.hpp:57
bool steepest_descent_step(const ParticleRange &particles)
Steepest descent integrator.
void stokesian_dynamics_step_1(ParticleRangeStokesian const &particles, StokesianThermostat const &stokesian, double time_step, double kT)
void velocity_verlet_rotator_1(Particle &p, double time_step)
void velocity_verlet_propagator_2(Particle &p, double time_step)
Final integration step of the Velocity Verlet integrator.
void velocity_verlet_propagator_1(Particle &p, double time_step)
Propagate the velocities and positions.
void velocity_verlet_rotator_2(Particle &p, double time_step)
void velocity_verlet_npt_step_2(ParticleRangeNPT const &particles, IsotropicNptThermostat const &npt_iso, double time_step)
Final integration step of the Velocity Verlet+NpT integrator.
void velocity_verlet_npt_step_1(ParticleRangeNPT const &particles, IsotropicNptThermostat const &npt_iso, double time_step, System::System &system)
Special propagator for NpT isotropic.