Skip to content

Commit

Permalink
Merge pull request #1925 from glotzerlab/fix-muvt-gibbs
Browse files Browse the repository at this point in the history 
Ensure that partitions have different seeds in MuVT.
  • Loading branch information
@joaander
joaander authored Nov 4, 2024
2 parents dc042e5 + dcba850 commit ec0b84c
Show file tree
Hide file tree
Showing 3 changed files with 94 additions and 37 deletions.
2 changes: 2 additions & 0 deletions CHANGELOG.rst
View file
Original file line number Diff line number Diff line change
Expand Up @@ -14,6 +14,8 @@ Change Log

* Correct compile errors with ``-DENABLE_GPU=on -DHOOMD_GPU_PLATFORM=HIP``
(`#1920 <https://github.com/glotzerlab/hoomd-blue/pull/1915>`__)
* Ensure that users set unique seeds on all partitions when performing Gibbs ensemble simulations
(`#1925 <https://github.com/glotzerlab/hoomd-blue/pull/1925>`__)

4.9.0 (2024-10-29)
^^^^^^^^^^^^^^^^^^
Expand Down
6 changes: 3 additions & 3 deletions hoomd/RNGIdentifiers.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -34,16 +34,16 @@ struct RNGIdentifier
static const uint8_t UpdaterClusters = 8;
static const uint8_t UpdaterClustersPairwise = 9;
static const uint8_t UpdaterExternalFieldWall = 10;
static const uint8_t UpdaterMuVT = 11;
static const uint8_t UpdaterMuVTGroup = 11;
static const uint8_t UpdaterMuVTDepletants1 = 12;
static const uint8_t UpdaterMuVTDepletants2 = 13;
static const uint8_t UpdaterMuVTDepletants3 = 14;
static const uint8_t UpdaterMuVTDepletants4 = 15;
static const uint8_t UpdaterMuVTDepletants5 = 16;
static const uint8_t UpdaterMuVTDepletants6 = 17;
static const uint8_t UpdaterMuVTPoisson = 18;
static const uint8_t UpdaterMuVTBox1 = 19;
static const uint8_t UpdaterMuVTBox2 = 20;
static const uint8_t UpdaterMuVTInsertRemove = 19;
static const uint8_t Unused1 = 20;
static const uint8_t ActiveForceCompute = 21;
static const uint8_t EvaluatorPairDPDThermo = 22;
static const uint8_t IntegrationMethodTwoStep = 23;
Expand Down
123 changes: 89 additions & 34 deletions hoomd/hpmc/UpdaterMuVT.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -156,6 +156,7 @@ template<class Shape> class UpdaterMuVT : public Updater
m_mc; //!< The MC Integrator this Updater is associated with
unsigned int m_npartition; //!< The number of partitions to use for Gibbs ensemble
bool m_gibbs; //!< True if we simulate a Gibbs ensemble
uint16_t m_move_type_seed; //!< Random number seed to use for move types.

GPUVector<Scalar4> m_postype_backup; //!< Backup of postype array

Expand Down Expand Up @@ -347,6 +348,8 @@ UpdaterMuVT<Shape>::UpdaterMuVT(std::shared_ptr<SystemDefinition> sysdef,
m_gibbs = true;
}

m_move_type_seed = m_sysdef->getSeed();

#ifdef ENABLE_MPI
if (m_gibbs)
{
Expand All @@ -361,6 +364,52 @@ UpdaterMuVT<Shape>::UpdaterMuVT(std::shared_ptr<SystemDefinition> sysdef,

m_exec_conf->msg->notice(5) << "Constructing UpdaterMuVT: Gibbs ensemble with "
<< m_npartition << " partitions" << std::endl;

// Ensure that the user sets unique seeds on all partitions so that local trial moves
// are decorrelated.
unsigned int my_partition = this->m_exec_conf->getPartition();
unsigned int my_group = this->m_exec_conf->getPartition() / npartition;
uint16_t my_seed = this->m_sysdef->getSeed();

for (unsigned int check_partition = 0;
check_partition < this->m_exec_conf->getNPartitions();
check_partition++)
{
unsigned int check_group = check_partition / npartition;
uint16_t check_seed = my_seed;
MPI_Bcast(&check_seed,
1,
MPI_UINT16_T,
check_partition * m_exec_conf->getMPIConfig()->getNRanks(),
m_exec_conf->getHOOMDWorldMPICommunicator());

if (my_group == check_group && my_partition != check_partition && my_seed == check_seed)
{
std::ostringstream s;
s << "Each partition within a group must set a unique seed. ";
s << "Partition " << check_partition << "'s " << "seed (" << check_seed << ") ";
s << "matches partition " << my_partition << "'s";
throw std::runtime_error(s.str());
}
}

// synchronize move types across all ranks within each group
for (unsigned int group = 0; group < this->m_exec_conf->getNPartitions() / npartition;
group++)
{
uint16_t tmp = m_move_type_seed;
MPI_Bcast(&tmp,
1,
MPI_UINT16_T,
group * npartition * this->m_exec_conf->getNRanks(),
m_exec_conf->getHOOMDWorldMPICommunicator());

unsigned int my_group = this->m_exec_conf->getPartition() / npartition;
if (my_group == group)
{
m_move_type_seed = tmp;
}
}
}
else
#endif
Expand Down Expand Up @@ -865,9 +914,11 @@ template<class Shape> void UpdaterMuVT<Shape>::update(uint64_t timestep)

// initialize random number generator
unsigned int group = (m_exec_conf->getPartition() / m_npartition);
unsigned int partition = (m_exec_conf->getPartition() % m_npartition);

hoomd::RandomGenerator rng(
hoomd::Seed(hoomd::RNGIdentifier::UpdaterMuVT, timestep, this->m_sysdef->getSeed()),
// Make a RNG that is seeded the same across the whole group
hoomd::RandomGenerator rng_group(
hoomd::Seed(hoomd::RNGIdentifier::UpdaterMuVTGroup, timestep, m_move_type_seed),
hoomd::Counter(group));

bool active = true;
Expand All @@ -884,31 +935,30 @@ template<class Shape> void UpdaterMuVT<Shape>::update(uint64_t timestep)
// the other MPI partition
if (m_gibbs)
{
unsigned int p = m_exec_conf->getPartition() % m_npartition;

// choose a random pair of communicating boxes
src = hoomd::UniformIntDistribution(m_npartition - 1)(rng);
dest = hoomd::UniformIntDistribution(m_npartition - 2)(rng);
src = hoomd::UniformIntDistribution(m_npartition - 1)(rng_group);
dest = hoomd::UniformIntDistribution(m_npartition - 2)(rng_group);
if (src <= dest)
dest++;

if (p == src)
if (partition == src)
{
m_gibbs_other = (dest + group * m_npartition) * m_exec_conf->getNRanks();
mod = 0;
}
if (p == dest)
if (partition == dest)
{
m_gibbs_other = (src + group * m_npartition) * m_exec_conf->getNRanks();
mod = 1;
}
if (p != src && p != dest)
if (partition != src && partition != dest)
{
active = false;
}

// order the expanded ensembles
volume_move = hoomd::detail::generate_canonical<Scalar>(rng) < m_volume_move_probability;
Scalar r = hoomd::detail::generate_canonical<Scalar>(rng_group);
volume_move = r < m_volume_move_probability;

if (active && m_exec_conf->getRank() == 0)
{
Expand Down Expand Up @@ -973,7 +1023,14 @@ template<class Shape> void UpdaterMuVT<Shape>::update(uint64_t timestep)
#endif

// whether we insert or remove a particle
bool insert = m_gibbs ? mod : hoomd::UniformIntDistribution(1)(rng);
bool insert = m_gibbs ? mod : hoomd::UniformIntDistribution(1)(rng_group);

// Use a partition specific RNG stream on each partition in Gibbs ensembles.
hoomd::RandomGenerator rng_insert_remove(
hoomd::Seed(hoomd::RNGIdentifier::UpdaterMuVTInsertRemove,
timestep,
this->m_sysdef->getSeed()),
hoomd::Counter(group, partition));

if (insert)
{
Expand All @@ -992,7 +1049,7 @@ template<class Shape> void UpdaterMuVT<Shape>::update(uint64_t timestep)
{
// choose a random particle type out of those being inserted or removed
type = m_transfer_types[hoomd::UniformIntDistribution(
(unsigned int)(m_transfer_types.size() - 1))(rng)];
(unsigned int)(m_transfer_types.size() - 1))(rng_insert_remove)];
}
else
{
Expand Down Expand Up @@ -1045,23 +1102,23 @@ template<class Shape> void UpdaterMuVT<Shape>::update(uint64_t timestep)

// Propose a random position uniformly in the box
Scalar3 f;
f.x = hoomd::detail::generate_canonical<Scalar>(rng);
f.y = hoomd::detail::generate_canonical<Scalar>(rng);
f.x = hoomd::detail::generate_canonical<Scalar>(rng_insert_remove);
f.y = hoomd::detail::generate_canonical<Scalar>(rng_insert_remove);
if (m_sysdef->getNDimensions() == 2)
{
f.z = Scalar(0.5);
}
else
{
f.z = hoomd::detail::generate_canonical<Scalar>(rng);
f.z = hoomd::detail::generate_canonical<Scalar>(rng_insert_remove);
}
vec3<Scalar> pos_test = vec3<Scalar>(m_pdata->getGlobalBox().makeCoordinates(f));

Shape shape_test(quat<Scalar>(), param);
if (shape_test.hasOrientation())
{
// set particle orientation
shape_test.orientation = generateRandomOrientation(rng, ndim);
shape_test.orientation = generateRandomOrientation(rng_insert_remove, ndim);
}

if (m_gibbs)
Expand Down Expand Up @@ -1140,7 +1197,8 @@ template<class Shape> void UpdaterMuVT<Shape>::update(uint64_t timestep)
bool accept = false;
if (nonzero)
{
accept = (hoomd::detail::generate_canonical<double>(rng) < exp(lnboltzmann));
accept = (hoomd::detail::generate_canonical<double>(rng_insert_remove)
< exp(lnboltzmann));
}

#ifdef ENABLE_MPI
Expand Down Expand Up @@ -1190,24 +1248,19 @@ template<class Shape> void UpdaterMuVT<Shape>::update(uint64_t timestep)
// try removing a particle
unsigned int tag = UINT_MAX;

// in Gibbs ensemble, we should not use correlated random numbers with box 1
hoomd::RandomGenerator rng_local(hoomd::Seed(hoomd::RNGIdentifier::UpdaterMuVTBox1,
timestep,
this->m_sysdef->getSeed()),
hoomd::Counter(group));

// choose a random particle type out of those being transferred
assert(m_transfer_types.size() > 0);
unsigned int type = m_transfer_types[hoomd::UniformIntDistribution(
(unsigned int)(m_transfer_types.size() - 1))(rng_local)];
(unsigned int)(m_transfer_types.size() - 1))(rng_insert_remove)];

// choose a random particle of that type
unsigned int nptl_type = getNumParticlesType(type);

if (nptl_type)
{
// get random tag of given type
unsigned int type_offset = hoomd::UniformIntDistribution(nptl_type - 1)(rng_local);
unsigned int type_offset
= hoomd::UniformIntDistribution(nptl_type - 1)(rng_insert_remove);
tag = getNthTypeTag(type, type_offset);
}

Expand Down Expand Up @@ -1319,8 +1372,8 @@ template<class Shape> void UpdaterMuVT<Shape>::update(uint64_t timestep)
// apply acceptance criterion
if (nonzero)
{
accept
= (hoomd::detail::generate_canonical<double>(rng_local) < exp(lnboltzmann));
accept = (hoomd::detail::generate_canonical<double>(rng_insert_remove)
< exp(lnboltzmann));
}
else
{
Expand Down Expand Up @@ -1410,18 +1463,20 @@ template<class Shape> void UpdaterMuVT<Shape>::update(uint64_t timestep)

if (mod == 0)
{
Scalar ln_V_new = log(V / V_other)
+ (hoomd::detail::generate_canonical<Scalar>(rng) - Scalar(0.5))
* m_max_vol_rescale;
Scalar ln_V_new
= log(V / V_other)
+ (hoomd::detail::generate_canonical<Scalar>(rng_group) - Scalar(0.5))
* m_max_vol_rescale;
V_new = (V + V_other) * exp(ln_V_new) / (Scalar(1.0) + exp(ln_V_new));
V_new_other
= (V + V_other) * (Scalar(1.0) - exp(ln_V_new) / (Scalar(1.0) + exp(ln_V_new)));
}
else
{
Scalar ln_V_new = log(V_other / V)
+ (hoomd::detail::generate_canonical<Scalar>(rng) - Scalar(0.5))
* m_max_vol_rescale;
Scalar ln_V_new
= log(V_other / V)
+ (hoomd::detail::generate_canonical<Scalar>(rng_group) - Scalar(0.5))
* m_max_vol_rescale;
V_new
= (V + V_other) * (Scalar(1.0) - exp(ln_V_new) / (Scalar(1.0) + exp(ln_V_new)));
}
Expand Down Expand Up @@ -1538,7 +1593,7 @@ template<class Shape> void UpdaterMuVT<Shape>::update(uint64_t timestep)
+ log(V_new_other / V_other) * (Scalar)(other_ndof + 1) + lnb
+ other_lnb;

accept = hoomd::detail::generate_canonical<double>(rng) < exp(arg);
accept = hoomd::detail::generate_canonical<double>(rng_group) < exp(arg);
accept &= !(has_overlaps || other_result);

// communicate if accepted
Expand Down

0 comments on commit ec0b84c

Please sign in to comment.