Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Refactor thermostat infrastructure #3888

Merged
merged 22 commits into from
Sep 18, 2020
Merged
Show file tree
Hide file tree
Changes from 20 commits
Commits
Show all changes
22 commits
Select commit Hold shift + click to select a range
9f9fc26
utils: Fixed parameter type for uint64_t u32_to_u64
fweik Mar 17, 2020
f956bac
core: Use seed+salt in key of Philox
fweik Mar 17, 2020
ea031c6
core: Use Philox infrastructure for SD thermostat
jngrad Sep 1, 2020
4af67b4
docs: Document unified thermostat infrastructure
jngrad Sep 1, 2020
0c9a9b1
core: Make statistical unit tests faster
jngrad Sep 3, 2020
69c999a
core: Make statistical unit tests faster
jngrad Sep 3, 2020
0c97c61
core: Break integrator/thermostat cyclic dependency
jngrad Sep 4, 2020
fc7d191
core: Cleanup code
jngrad Sep 4, 2020
ae13462
core: Move NpT functions to dedicated file
jngrad Sep 4, 2020
dbe96b0
core: Split thermostats from integrators
jngrad Sep 4, 2020
c475c38
core: Header cleanup
jngrad Sep 4, 2020
5428fca
core: Remove redundant thermostat free functions
jngrad Sep 8, 2020
df90e9c
tests: Improve interaction checkpointing tests
jngrad Sep 10, 2020
9c454a6
python: Pickle counter without altering params
jngrad Sep 16, 2020
cc2e0df
docs: Update thermostat counter/seed explanation
jngrad Sep 16, 2020
6073521
core: Use boost::mpi
jngrad Sep 16, 2020
2366434
Merge branch 'python' into thermostats-correlation
jngrad Sep 16, 2020
11b9b0b
core: Remove unused NpT parameter
jngrad Sep 16, 2020
d21c7a4
core: Move NpT communication to npt.cpp
jngrad Sep 16, 2020
5d6e1d3
core: Guard NpT callback
jngrad Sep 16, 2020
8d28cfa
tests: Remove calls to private methods
jngrad Sep 18, 2020
3323d57
Merge branch 'python' into thermostats-correlation
kodiakhq[bot] Sep 18, 2020
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
28 changes: 10 additions & 18 deletions doc/sphinx/system_setup.rst
Original file line number Diff line number Diff line change
Expand Up @@ -212,6 +212,13 @@ the temperature of a thermostat, you actually do not define the
temperature, but the value of the thermal energy :math:`k_B T` in the
current unit system (see the discussion on units, Section :ref:`On units`).

All thermostats have a ``seed`` argument that controls the state of the random
number generator (Philox Counter-based RNG). This seed is required on first
activation of a thermostat, unless stated otherwise. It can be omitted in
subsequent calls of the method that activates the same thermostat. The random
sequence also depends on the thermostats counters that are
incremented after each integration step.

.. _Langevin thermostat:

Langevin thermostat
Expand Down Expand Up @@ -253,12 +260,6 @@ The distribution of :math:`\overline{\eta}` is uniform and satisfies

.. math:: <\overline{\eta}> = 0 , <\overline{\eta}\overline{\eta}> = 1/dt

The keyword ``seed`` controls the state of the random number generator (Philox
Counter-based RNG) and is required on first activation of the thermostat. It
can be omitted in subsequent calls of ``set_langevin()``. It is the user's
responsibility to decide whether the thermostat should be deterministic (by
using a fixed seed) or not (by using a randomized seed).

If the feature ``ROTATION`` is compiled in, the rotational degrees of freedom are
also coupled to the thermostat. If only the first two arguments are
specified then the friction coefficient for the rotation is set to the
Expand Down Expand Up @@ -301,10 +302,8 @@ The LB fluid can be used to thermalize particles, while also including their hyd
The LB thermostat expects an instance of either :class:`espressomd.lb.LBFluid` or :class:`espressomd.lb.LBFluidGPU`.
Temperature is set via the ``kT`` argument of the LB fluid.

Furthermore a ``seed`` has to be given for the
thermalization of the particle coupling. The magnitude of the frictional coupling can be adjusted by
the parameter ``gamma``.
To enable the LB thermostat, use::
The magnitude of the frictional coupling can be adjusted by the
parameter ``gamma``. To enable the LB thermostat, use::

import espressomd
import espressomd.lb
Expand Down Expand Up @@ -471,12 +470,6 @@ Note: the rotational Brownian dynamics implementation is compatible with particl
the isotropic moment of inertia tensor only. Otherwise, the viscous terminal angular velocity
is not defined, i.e. it has no constant direction over the time.

The keyword ``seed`` controls the state of the random number generator (Philox
Counter-based RNG) and is required on first activation of the thermostat. It
can be omitted in subsequent calls of ``set_brownian()``. It is the user's
responsibility to decide whether the thermostat should be deterministic (by
using a fixed seed) or not (by using a randomized seed).

.. _Stokesian thermostat:

Stokesian thermostat
Expand All @@ -501,8 +494,7 @@ needs to be activated via::
system.integrator.run(100)

where ``kT`` denotes the desired temperature of the system, and ``seed`` the
seed for the random number generator of the Stokesian Dynamics thermostat.
It is independent from the other random number generators in |es|.
seed for the random number generator.


.. _CUDA:
Expand Down
1 change: 1 addition & 0 deletions src/core/bonded_interactions/thermalized_bond.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -27,6 +27,7 @@
#include "bonded_interaction_data.hpp"
#include "communication.hpp"
#include "global.hpp"
#include "integrate.hpp"

int n_thermalized_bonds = 0;

Expand Down
4 changes: 2 additions & 2 deletions src/core/bonded_interactions/thermalized_bond.hpp
Original file line number Diff line number Diff line change
Expand Up @@ -33,7 +33,6 @@ extern int n_thermalized_bonds;

#include "Particle.hpp"
#include "bonded_interaction_data.hpp"
#include "integrate.hpp"
#include "random.hpp"
#include "thermostat.hpp"

Expand Down Expand Up @@ -77,7 +76,8 @@ thermalized_bond_forces(Particle const &p1, Particle const &p2,
Utils::Vector3d force1{};
Utils::Vector3d force2{};
auto const noise = Random::noise_uniform<RNGSalt::THERMALIZED_BOND>(
thermalized_bond.rng_get(), p1.p.identity, p2.p.identity);
thermalized_bond.rng_counter(), thermalized_bond.rng_seed(),
p1.p.identity, p2.p.identity);

for (int i = 0; i < 3; i++) {
double force_lv_com, force_lv_dist;
Expand Down
16 changes: 6 additions & 10 deletions src/core/communication.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -97,7 +97,6 @@ int n_nodes = -1;
CB(mpi_bcast_coulomb_params_slave) \
CB(mpi_remove_particle_slave) \
CB(mpi_rescale_particles_slave) \
CB(mpi_bcast_nptiso_geom_slave) \
CB(mpi_bcast_steepest_descent_worker) \
CB(mpi_bcast_cuda_global_part_vars_slave) \
CB(mpi_resort_particles_slave) \
Expand Down Expand Up @@ -496,17 +495,14 @@ void mpi_set_use_verlet_lists(bool use_verlet_lists) {

/*************** BCAST NPTISO GEOM *****************/

void mpi_bcast_nptiso_geom() {
mpi_call(mpi_bcast_nptiso_geom_slave, -1, 0);
mpi_bcast_nptiso_geom_slave(-1, 0);
}
#ifdef NPT
REGISTER_CALLBACK(mpi_bcast_nptiso_geom_worker)

void mpi_bcast_nptiso_geom_slave(int, int) {
MPI_Bcast(&nptiso.geometry, 1, MPI_INT, 0, comm_cart);
MPI_Bcast(&nptiso.dimension, 1, MPI_INT, 0, comm_cart);
MPI_Bcast(&nptiso.cubic_box, 1, MPI_LOGICAL, 0, comm_cart);
MPI_Bcast(&nptiso.non_const_dim, 1, MPI_INT, 0, comm_cart);
void mpi_bcast_nptiso_geom() {
mpi_call(mpi_bcast_nptiso_geom_worker, -1, 0);
mpi_bcast_nptiso_geom_worker(-1, 0);
}
#endif

/*************** BCAST STEEPEST DESCENT *****************/

Expand Down
5 changes: 3 additions & 2 deletions src/core/constraints/ShapeBasedConstraint.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -24,6 +24,7 @@
#include "errorhandling.hpp"
#include "forces_inline.hpp"
#include "nonbonded_interactions/nonbonded_interaction_data.hpp"
#include "thermostat.hpp"

namespace Constraints {
Utils::Vector3d ShapeBasedConstraint::total_force() const {
Expand Down Expand Up @@ -77,7 +78,7 @@ ParticleForce ShapeBasedConstraint::force(Particle const &p,
force1 +=
dpd_pair_force(p, part_rep, ia_params, dist_vec, dist, dist * dist);
// Additional use of DPD here requires counter increase
dpd_rng_counter_increment();
dpd.rng_increment();
KaiSzuttor marked this conversation as resolved.
Show resolved Hide resolved
}
#endif
} else if (m_penetrable && (dist <= 0)) {
Expand All @@ -89,7 +90,7 @@ ParticleForce ShapeBasedConstraint::force(Particle const &p,
force1 += dpd_pair_force(p, part_rep, ia_params, dist_vec, dist,
dist * dist);
// Additional use of DPD here requires counter increase
dpd_rng_counter_increment();
dpd.rng_increment();
}
#endif
}
Expand Down
4 changes: 2 additions & 2 deletions src/core/dpd.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -47,9 +47,9 @@ using Utils::Vector3d;
* seed-per-node)
*/
Vector3d dpd_noise(uint32_t pid1, uint32_t pid2) {
extern DPDThermostat dpd;
return Random::noise_uniform<RNGSalt::SALT_DPD>(
dpd.rng_get(), (pid1 < pid2) ? pid2 : pid1, (pid1 < pid2) ? pid1 : pid2);
dpd.rng_counter(), dpd.rng_seed(), (pid1 < pid2) ? pid2 : pid1,
(pid1 < pid2) ? pid1 : pid2);
}

int dpd_set_params(int part_type_a, int part_type_b, double gamma, double k,
Expand Down
8 changes: 4 additions & 4 deletions src/core/forces.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -49,7 +49,7 @@

ActorList forceActors;

void init_forces(const ParticleRange &particles) {
void init_forces(const ParticleRange &particles, double time_step) {
ESPRESSO_PROFILER_CXX_MARK_FUNCTION;
/* The force initialization depends on the used thermostat and the
thermodynamic ensemble */
Expand All @@ -63,7 +63,7 @@ void init_forces(const ParticleRange &particles) {
set torque to zero for all and rescale quaternions
*/
for (auto &p : particles) {
p.f = init_local_particle_force(p);
p.f = init_local_particle_force(p, time_step);
}

/* initialize ghost forces with zero
Expand All @@ -80,7 +80,7 @@ void init_forces_ghosts(const ParticleRange &particles) {
}
}

void force_calc(CellStructure &cell_structure) {
void force_calc(CellStructure &cell_structure, double time_step) {
ESPRESSO_PROFILER_CXX_MARK_FUNCTION;

espressoSystemInterface.update();
Expand All @@ -94,7 +94,7 @@ void force_calc(CellStructure &cell_structure) {
#ifdef ELECTROSTATICS
iccp3m_iteration(particles, cell_structure.ghost_particles());
#endif
init_forces(particles);
init_forces(particles, time_step);

for (auto &forceActor : forceActors) {
forceActor->computeForces(espressoSystemInterface);
Expand Down
11 changes: 2 additions & 9 deletions src/core/forces.hpp
Original file line number Diff line number Diff line change
Expand Up @@ -38,15 +38,9 @@

extern ActorList forceActors;

/** \name Exported Functions */
/************************************************************/
/*@{*/

/******************* forces.cpp *******************/

/** initialize real particle forces with thermostat forces and
ghost particle forces with zero. */
void init_forces(const ParticleRange &particles);
void init_forces(const ParticleRange &particles, double time_step);

/** Set forces of all ghosts to zero */
void init_forces_ghosts(const ParticleRange &particles);
Expand All @@ -61,10 +55,9 @@ void init_forces_ghosts(const ParticleRange &particles);
* <li> Calculate long range interaction forces
* </ol>
*/
void force_calc(CellStructure &cell_structure);
void force_calc(CellStructure &cell_structure, double time_step);

/** Calculate long range forces (P3M, ...). */
void calc_long_range_forces(const ParticleRange &particles);
/*@}*/

#endif
16 changes: 9 additions & 7 deletions src/core/forces_inline.hpp
Original file line number Diff line number Diff line change
Expand Up @@ -39,7 +39,6 @@
#include "forces.hpp"
#include "immersed_boundary/ibm_tribend.hpp"
#include "immersed_boundary/ibm_triel.hpp"
#include "integrators/langevin_inline.hpp"
#include "nonbonded_interactions/bmhtf-nacl.hpp"
#include "nonbonded_interactions/buckingham.hpp"
#include "nonbonded_interactions/gaussian.hpp"
Expand All @@ -61,6 +60,7 @@
#include "object-in-fluid/oif_local_forces.hpp"
#include "rotation.hpp"
#include "thermostat.hpp"
#include "thermostats/langevin_inline.hpp"

#ifdef DIPOLES
#include "electrostatics_magnetostatics/dipole_inline.hpp"
Expand Down Expand Up @@ -100,25 +100,27 @@ inline ParticleForce external_force(Particle const &p) {
return f;
}

inline ParticleForce thermostat_force(Particle const &p) {
inline ParticleForce thermostat_force(Particle const &p, double time_step) {
extern LangevinThermostat langevin;
if (!(thermo_switch & THERMO_LANGEVIN)) {
return {};
}

#ifdef ROTATION
return {friction_thermo_langevin(langevin, p),
return {friction_thermo_langevin(langevin, p, time_step),
p.p.rotation ? convert_vector_body_to_space(
p, friction_thermo_langevin_rotation(langevin, p))
p, friction_thermo_langevin_rotation(langevin, p,
time_step))
: Utils::Vector3d{}};
#else
return friction_thermo_langevin(langevin, p);
return friction_thermo_langevin(langevin, p, time_step);
#endif
}

/** Initialize the forces for a real particle */
inline ParticleForce init_local_particle_force(Particle const &part) {
return thermostat_force(part) + external_force(part);
inline ParticleForce init_local_particle_force(Particle const &part,
double time_step) {
return thermostat_force(part, time_step) + external_force(part);
}

inline Utils::Vector3d calc_non_bonded_pair_force_parts(
Expand Down
4 changes: 1 addition & 3 deletions src/core/global.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -29,6 +29,7 @@
#include "event.hpp"
#include "grid.hpp"
#include "grid_based_algorithms/lb_interface.hpp"
#include "integrate.hpp"
#include "nonbonded_interactions/nonbonded_interaction_data.hpp"
#include "npt.hpp"
#include "object-in-fluid/oif_global_forces.hpp"
Expand All @@ -45,9 +46,6 @@
#include <unordered_map>

extern double force_cap;
extern LangevinThermostat langevin;
extern BrownianThermostat brownian;
extern IsotropicNptThermostat npt_iso;

namespace {

Expand Down
4 changes: 2 additions & 2 deletions src/core/grid_based_algorithms/lb.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -42,10 +42,10 @@
#include "nonbonded_interactions/nonbonded_interaction_data.hpp"
#include "random.hpp"

#include "utils/u32_to_u64.hpp"
#include <utils/Counter.hpp>
#include <utils/index.hpp>
#include <utils/math/matrix_vector_product.hpp>
#include <utils/u32_to_u64.hpp>
#include <utils/uniform.hpp>
using Utils::get_linear_index;
#include <utils/constants.hpp>
Expand Down Expand Up @@ -841,7 +841,7 @@ lb_thermalize_modes(Lattice::index_t index, const std::array<T, 19> &modes,
using rng_type = r123::Philox4x64;
using ctr_type = rng_type::ctr_type;

const r123::Philox4x64::ctr_type c{
const ctr_type c{
{rng_counter->value(), static_cast<uint64_t>(RNGSalt::FLUID)}};
const T rootdensity =
std::sqrt(std::fabs(modes[0] + lb_parameters.density));
Expand Down
2 changes: 1 addition & 1 deletion src/core/grid_based_algorithms/lb_particle_coupling.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -279,7 +279,7 @@ void lb_lbcoupling_calc_particle_lattice_ia(
auto f_random = [noise_amplitude](int id) -> Utils::Vector3d {
if (noise_amplitude > 0.0) {
return Random::noise_uniform<RNGSalt::PARTICLES>(
lb_particle_coupling.rng_counter_coupling->value(), id);
lb_particle_coupling.rng_counter_coupling->value(), 0, id);
}
return {};
};
Expand Down
Loading