From 4fe5cd16fdb78572aa03fe3f333e9f0be2f9a649 Mon Sep 17 00:00:00 2001 From: Tim Moore Date: Sun, 15 Dec 2024 21:23:52 -0500 Subject: [PATCH] Update docs and remove some experimental vmmc features Best to get the main implementation into HOOMD and then test experimental features since it will be easier to validate the straightforward implementation. --- hoomd/hpmc/UpdaterVirtualMoveMonteCarlo.h | 89 +--------- hoomd/hpmc/update.py | 188 ++++++++++++---------- 2 files changed, 109 insertions(+), 168 deletions(-) diff --git a/hoomd/hpmc/UpdaterVirtualMoveMonteCarlo.h b/hoomd/hpmc/UpdaterVirtualMoveMonteCarlo.h index 9f90414373..5ae39fb19b 100644 --- a/hoomd/hpmc/UpdaterVirtualMoveMonteCarlo.h +++ b/hoomd/hpmc/UpdaterVirtualMoveMonteCarlo.h @@ -108,16 +108,6 @@ class UpdaterVMMC : public Updater return m_attempts_per_particle; } - void setBetaFicticious(std::shared_ptr beta_ficticious) - { - m_beta_ficticious = beta_ficticious; - } - - std::shared_ptr getBetaFicticious() - { - return m_beta_ficticious; - } - void setTranslationMoveProbability(LongReal p) { m_translation_move_probability = p; @@ -214,26 +204,6 @@ class UpdaterVMMC : public Updater m_instance = instance; } - bool getStaticLinkingMode() - { - return m_static_linking_mode; - } - - void setStaticLinkingMode(bool b) - { - m_static_linking_mode = b; - } - - bool getRandomBetaFicticious() - { - return m_random_beta_ficticious; - } - - void setRandomBetaFicticious(bool b) - { - m_random_beta_ficticious = b; - } - /*! \param mode 0 -> Absolute count, 1 -> relative to the start of the run, 2 -> relative to the last executed step \return The current state of the acceptance counters @@ -253,14 +223,11 @@ class UpdaterVMMC : public Updater protected: std::shared_ptr< IntegratorHPMCMono > m_mc; //!< HPMC integrator unsigned int m_attempts_per_particle; //!< Number of attempted moves per particle each time update() is called - std::shared_ptr m_beta_ficticious; //!< Ficticious inverse temperature for creating links between particles LongReal m_translation_move_probability; unsigned int m_maximum_allowed_cluster_size; LongReal m_cluster_size_distribution_prefactor; std::string m_cluster_size_limit_mode; bool m_always_rebuild_tree; - bool m_static_linking_mode; - bool m_random_beta_ficticious; detail::UpdateOrder m_update_order; @@ -288,7 +255,7 @@ UpdaterVMMC::UpdaterVMMC( std::shared_ptr trigger, std::shared_ptr > mc, std::shared_ptr beta_ficticious - ) : Updater(sysdef, trigger), m_mc(mc), m_beta_ficticious(beta_ficticious), m_update_order(m_pdata->getN()) + ) : Updater(sysdef, trigger), m_mc(mc), m_update_order(m_pdata->getN()) { m_exec_conf->msg->notice(5) << "Constructing UpdaterVMMC" << std::endl; @@ -340,7 +307,6 @@ void UpdaterVMMC::update(uint64_t timestep) m_count_step_start = m_count_total; // get stuff needed for the calculation - LongReal current_beta_ficticious = m_static_linking_mode ? m_beta_ficticious->operator()(timestep) : 1.0; const Scalar kT = m_mc->getKT()->operator()(timestep); const LongReal min_core_radius = m_mc->getMinCoreDiameter() * LongReal(0.5); const auto& pair_energy_search_radius = m_mc->getPairEnergySearchRadius(); @@ -404,10 +370,6 @@ void UpdaterVMMC::update(uint64_t timestep) // generate the virtual move hoomd::RandomGenerator rng_i(hoomd::Seed(hoomd::RNGIdentifier::UpdaterVMMC, timestep, user_seed), hoomd::Counter(seed_idx, (m_exec_conf->getRank() << 16) + m_instance, i_trial)); - if (m_static_linking_mode && m_random_beta_ficticious) - { - current_beta_ficticious = hoomd::UniformDistribution()(rng_i); - } LongReal move_type_select = hoomd::UniformDistribution()(rng_i); bool move_type_translate = move_type_select < m_translation_move_probability; vec3 virtual_translate_move(0, 0, 0); @@ -470,8 +432,6 @@ void UpdaterVMMC::update(uint64_t timestep) pos_linker_after_move_primary_image = rotate(virtual_rotate_move, pos_linker - center_of_rotation) + center_of_rotation; m_exec_conf->msg->notice(5) << "center of rotation for particle " << i_linker << " = " << center_of_rotation.x << " " << center_of_rotation.y << " " << center_of_rotation.z << std::endl; pos_linker_after_reverse_move_primary_image = rotate(conj(virtual_rotate_move), pos_linker - center_of_rotation) + center_of_rotation; - /* box.wrap(pos_linker_after_move_primary_image, _images); */ - /* box.wrap(pos_linker_after_reverse_move_primary_image, _images); */ } m_positions_after_move[i_linker] = pos_linker_after_move_primary_image; @@ -629,11 +589,6 @@ void UpdaterVMMC::update(uint64_t timestep) // account for probability of forming this link: denominator of p_link_probability_ratio *= 1 p_link_probability_ratio /= 1.0; - if (m_static_linking_mode) - { - continue; - } - // account for probability of forming this link on the reverse move: // numerator of p_link_probability_ratio p_ij_reverse gets multiplied by p_ij_reverse // if overlap after reverse move of i, p_ij_reverse = 1, otherwise we look at energies @@ -672,7 +627,7 @@ void UpdaterVMMC::update(uint64_t timestep) shape_linkee.orientation, h_diameter.data[j_linkee], h_charge.data[j_linkee]); - p_ij_reverse = std::max(0.0, 1 - exp(-current_beta_ficticious*(u_ij_after_reverse_move_of_linker - u_ij) / kT)); + p_ij_reverse = std::max(0.0, 1 - exp(-(u_ij_after_reverse_move_of_linker - u_ij) / kT)); } if (p_ij_reverse == 0.0) { @@ -724,14 +679,7 @@ void UpdaterVMMC::update(uint64_t timestep) } LongReal p_ij_forward; - if (m_static_linking_mode) - { - p_ij_forward = 1.0 - exp(-current_beta_ficticious * abs(u_ij / kT)); - } - else - { - p_ij_forward = std::max(0.0, 1.0 - exp(-current_beta_ficticious * (u_ij_after_move - u_ij) / kT)); - } + p_ij_forward = std::max(0.0, 1.0 - exp(-(u_ij_after_move - u_ij) / kT)); Scalar r = hoomd::UniformDistribution()(rng_i); bool link_formed = r < p_ij_forward; if (link_formed) @@ -756,13 +704,6 @@ void UpdaterVMMC::update(uint64_t timestep) { m_cluster_data.update_cluster(j_linkee, center_of_rotation_image); } - if (m_static_linking_mode) - { - // do not need to modify deltaU since they're both in the cluster - // only need to just update the cluster data (done just above) and continue, so we just continue - continue; - } - // numerator of p_link_probability_ratio gets multiplied by p_ij_reverse Scalar p_ij_reverse; bool overlap_after_reverse_move = @@ -788,7 +729,7 @@ void UpdaterVMMC::update(uint64_t timestep) h_diameter.data[j_linkee], h_charge.data[j_linkee]); m_exec_conf->msg->notice(5) << " u_ij_after_reverse_move = " << u_ij_after_move << std::endl; - p_ij_reverse = std::max(0.0, 1.0 - exp(-current_beta_ficticious * (u_ij_after_reverse_move - u_ij) / kT)); + p_ij_reverse = std::max(0.0, 1.0 - exp(-(u_ij_after_reverse_move - u_ij) / kT)); if (p_ij_reverse == 0.0) { m_count_total.early_abort_no_reverse_link++; @@ -829,17 +770,9 @@ void UpdaterVMMC::update(uint64_t timestep) // relevant states for a non-linker pair are (x_i, x_j) and (x_i + dx, x_j) (where dx is // a generic virtual move) so the only relevant energies are beta_u_ij and beta_u_ij_after_move - // for the cluster cleaving version, the only thing we need to do here is update the potential link data container - // because we don't yet know if j_linkee will eventually join the cult - if (m_static_linking_mode) - { - m_potential_link_data.update(i_linker, j_linkee, (LongReal)u_ij_after_move - (LongReal)u_ij); - continue; - } - // account for probabilities of not forming this link Scalar q_ij_forward = 1.0 - p_ij_forward; - Scalar q_ij_reverse = std::min(1.0, fast::exp(-current_beta_ficticious * (u_ij - u_ij_after_move))); + Scalar q_ij_reverse = std::min(1.0, fast::exp(-(u_ij - u_ij_after_move) / kT)); if (q_ij_reverse == 0.0) { m_count_total.early_abort_forced_reverse_link++; @@ -910,14 +843,7 @@ void UpdaterVMMC::update(uint64_t timestep) m_exec_conf->msg->notice(5) << " beta*dU = " << deltaU / kT << std::endl; m_exec_conf->msg->notice(5) << " exp(-beta*dU) = " << exp(-deltaU / kT) << std::endl; LongReal p_acc; - if (m_static_linking_mode) - { - p_acc = std::min(1.0, exp((current_beta_ficticious - 1/kT) * deltaU)); - } - else - { - p_acc = std::min(1.0, exp(-deltaU / kT) * p_link_probability_ratio * p_failed_link_probability_ratio); - } + p_acc = std::min(1.0, exp(-deltaU / kT) * p_link_probability_ratio * p_failed_link_probability_ratio); Scalar r = hoomd::UniformDistribution()(rng_i); bool accept_cluster_move = r <= p_acc; m_exec_conf->msg->notice(4) << " VMMC p_acc: " << p_acc << std::endl; @@ -1016,7 +942,6 @@ template < class Shape> void export_UpdaterVirtualMoveMonteCarlo(pybind11::modul std::shared_ptr >()) .def("getCounters", &UpdaterVMMC::getCounters) - .def_property("beta_ficticious", &UpdaterVMMC::getBetaFicticious, &UpdaterVMMC::setBetaFicticious) .def_property("translation_move_probability", &UpdaterVMMC::getTranslationMoveProbability, &UpdaterVMMC::setTranslationMoveProbability) .def_property("maximum_trial_rotation", &UpdaterVMMC::getMaximumTrialRotation, &UpdaterVMMC::setMaximumTrialRotation) .def_property("maximum_trial_translation", &UpdaterVMMC::getMaximumTrialTranslation, &UpdaterVMMC::setMaximumTrialTranslation) @@ -1026,9 +951,7 @@ template < class Shape> void export_UpdaterVirtualMoveMonteCarlo(pybind11::modul .def_property("cluster_size_limit_mode", &UpdaterVMMC::getClusterSizeLimitMode, &UpdaterVMMC::setClusterSizeLimitMode) .def_property("cluster_size_distribution_prefactor", &UpdaterVMMC::getClusterSizeDistributionPrefactor, &UpdaterVMMC::setClusterSizeDistributionPrefactor) .def_property("always_rebuild_tree", &UpdaterVMMC::getAlwaysRebuildTree, &UpdaterVMMC::setAlwaysRebuildTree) - .def_property("static_linking_mode", &UpdaterVMMC::getStaticLinkingMode, &UpdaterVMMC::setStaticLinkingMode) .def_property("instance", &UpdaterVMMC::getInstance, &UpdaterVMMC::setInstance) - .def_property("random_beta_ficticious", &UpdaterVMMC::getRandomBetaFicticious, &UpdaterVMMC::setRandomBetaFicticious) ; } diff --git a/hoomd/hpmc/update.py b/hoomd/hpmc/update.py index 36b6b4dd28..ce1675cb05 100644 --- a/hoomd/hpmc/update.py +++ b/hoomd/hpmc/update.py @@ -679,52 +679,65 @@ class VirtualClusterMoves(Updater): """Apply virtual move Monte Carlo (VMMC) moves. Args: - trigger (hoomd.trigger.trigger_like): Apply cluster moves on triggered time steps. - - attempts_per_particle (int): Number of sweeps to run on each triggered timestep. - - beta_ficticious (float): Temperature to use for artificial virtual moves, between 0 and 1 inclusive. - - translation_move_probability (float): Probability of choosing a translation move, between 0 and 1 inclusive. - - maximum_trial_rotation (float): Maximum size of rotational trial moves, units? - - maximum_trial_translation (float): Maximum size of translational trial moves. - - maximum_trial_center_of_rotation_shift (float): Maximum translation of seed particle before applying rotational trial move. - - maximum_allowed_cluster_size (int): The largest allowable cluster size for acceptable cluster moves. - - cluster_size_distribution_prefactor (float): Prefactor for cluster size limit distribution. See below for more information about limiting the cluster size. - - cluster_size_limit_mode (str): How to limit cluster sizes. One of ``""``, ``"deterministic"``, or ``"probabilistic"``. - - static_linking_mode (bool): Whether or not to construct clusters statically. - - always_rebuild_tree (bool): If ``True``, rebuild the ``AABBTree`` any time a cluster move is accepted and applied. - - random_beta_ficticious (bool): I clearly need to clean up this class and get rid of any ideas I was testing. - - instance (int): When using multiple `VirtualClusterMoves` updaters in a single simulation, give each a unique value for `instance` so that they generate different streams of random numbers. - - The `VirtualClusterMoves` updater applies virtual cluster moves according to the algorithm described in `Whitelam and Geissler 2007 `__. - Here, we provide a practical overview of the algorithm; see `the original publication `__ for the theoretical details. - The general idea is that moves are attempted on clusters of particles, where the clusters are generated on-the-fly by attempting to form *links* between members of the growing cluster and nearby particles. - The formation of a link results in the addition of the non-member to the cluster. - - Specifically, the virtual cluster moves are created as follows: - first, a particle is selected as the seed particle and is the first member of the cluster. - The seed particle becomes the first *active linker*. - Next, a trial move is selected and applied to the active linker. - The change in pairwise energy between each potential *linkee* (particles in the system who have not been tested for cluster membership yet) and the active linker after applying the trial move to only the active linker determines the probability of that linkee joining the cluster. - For each particle that joins the cluster, the process of applying the trial move to it and determining new potential linkees is repeated until each member of the cluster has served as the active linker. - Once the cluster has been determined, the trial move is applied to the cluster as a whole and is accepted according the a balance-preserving criterion. - Each time the updater is triggered, every particle in the system acts as the seed particle ``attempts_per_particle`` times. - - Not clear yet wheter we'll keep the static and dynamic linking options. - I'm leaning towards not including them in the original PR and then adding them if testing proves that they'll be useful. - - .. rubric:: Box move types + trigger (hoomd.trigger.trigger_like): Apply cluster moves on triggered + time steps. + + attempts_per_particle (int): Number of sweeps to run on each triggered + timestep. + + translation_move_probability (float): Probability of choosing a + translation move, between 0 and 1 inclusive. + + maximum_trial_rotation (float): Maximum size of rotational trial moves. + + maximum_trial_translation (float): Maximum size of translational trial + moves. + + maximum_trial_center_of_rotation_shift (float): Maximum translation of + seed particle before applying rotational trial move. + + maximum_allowed_cluster_size (int): The largest allowable cluster size + for acceptable cluster moves. + + cluster_size_distribution_prefactor (float): Prefactor for cluster size + limit distribution. See below for more information about limiting the + cluster size. + + cluster_size_limit_mode (str): How to limit cluster sizes. One of + ``""``, ``"deterministic"``, or ``"probabilistic"``. + + always_rebuild_tree (bool): If ``True``, rebuild the ``AABBTree`` any + time a cluster move is accepted and applied. + + instance (int): When using multiple `VirtualClusterMoves` updaters in a + single simulation, give each a unique value for `instance` so that they + generate different streams of random numbers. + + The `VirtualClusterMoves` updater applies virtual cluster moves according + to the algorithm described in `Whitelam and Geissler 2007 + `__. Here, we provide a practical + overview of the algorithm; see `the original publication + `__ for the theoretical details. The + general idea is that moves are attempted on clusters of particles, where + the clusters are generated on-the-fly by attempting to form *links* between + members of the growing cluster and nearby particles. The formation of a + link results in the addition of the non-member to the cluster. + + Specifically, the virtual cluster moves are created as follows: first, a + particle is selected as the seed particle and is the first member of the + cluster. The seed particle becomes the first *active linker*. Next, a + trial move is selected and applied to the active linker. The change in + pairwise energy between each potential *linkee* (particles in the system + who have not been tested for cluster membership yet) and the active linker + after applying the trial move to only the active linker determines the + probability of that linkee joining the cluster. For each particle that + joins the cluster, the process of applying the trial move to it and + determining new potential linkees is repeated until each member of the + cluster has served as the active linker. Once the cluster has been + determined, the trial move is applied to the cluster as a whole and is + accepted according the a balance-preserving criterion. Each time the + updater is triggered, every particle in the system acts as the seed + particle ``attempts_per_particle`` times. {inherited} @@ -742,40 +755,42 @@ class VirtualClusterMoves(Updater): __doc__ = __doc__.replace("{inherited}", Updater._doc_inherited) - def __init__(self, trigger=1, attempts_per_particle=1, beta_ficticious=1.0, - translation_move_probability=0.5, - maximum_trial_rotation=0.5,maximum_trial_translation=1.0, - maximum_trial_center_of_rotation_shift=1.0, - maximum_allowed_cluster_size=0, - cluster_size_distribution_prefactor=1.0, - cluster_size_limit_mode=None, - static_linking_mode=False, - always_rebuild_tree=True, - random_beta_ficticious=False, - instance=0, - ): + def __init__( + self, + trigger=1, + attempts_per_particle=1, + translation_move_probability=0.5, + maximum_trial_rotation=0.5, + maximum_trial_translation=1.0, + maximum_trial_center_of_rotation_shift=1.0, + maximum_allowed_cluster_size=0, + cluster_size_distribution_prefactor=1.0, + cluster_size_limit_mode=None, + always_rebuild_tree=True, + instance=0, + ): super().__init__(trigger) if maximum_allowed_cluster_size is None: # set no limits on maximum allowed cluster size - # TODO: is there a more explicit way to do this insteading of using None? maximum_allowed_cluster_size = 0 if cluster_size_limit_mode is None: cluster_size_limit_mode = "" - param_dict = ParameterDict(attempts_per_particle=int(attempts_per_particle), - beta_ficticious=hoomd.variant.Variant, - translation_move_probability=float(translation_move_probability), - maximum_trial_rotation=float(maximum_trial_rotation), - maximum_trial_translation=float(maximum_trial_translation), - maximum_trial_center_of_rotation_shift=float(maximum_trial_center_of_rotation_shift), - maximum_allowed_cluster_size=int(maximum_allowed_cluster_size), - cluster_size_distribution_prefactor=float(cluster_size_distribution_prefactor), - cluster_size_limit_mode=str(cluster_size_limit_mode), - static_linking_mode=bool(static_linking_mode), - always_rebuild_tree=bool(always_rebuild_tree), - instance=int(instance), - random_beta_ficticious=bool(random_beta_ficticious), - ) - param_dict.update(dict(beta_ficticious=beta_ficticious)) + param_dict = ParameterDict( + attempts_per_particle=int(attempts_per_particle), + translation_move_probability=float(translation_move_probability), + maximum_trial_rotation=float(maximum_trial_rotation), + maximum_trial_translation=float(maximum_trial_translation), + maximum_trial_center_of_rotation_shift=float( + maximum_trial_center_of_rotation_shift + ), + maximum_allowed_cluster_size=int(maximum_allowed_cluster_size), + cluster_size_distribution_prefactor=float( + cluster_size_distribution_prefactor + ), + cluster_size_limit_mode=str(cluster_size_limit_mode), + always_rebuild_tree=bool(always_rebuild_tree), + instance=int(instance), + ) self._param_dict.update(param_dict) def _attach_hook(self): @@ -787,18 +802,21 @@ def _attach_hook(self): cpp_cls_name += integrator.__class__.__name__ cpp_cls = getattr(_hpmc, cpp_cls_name) self._cpp_obj = cpp_cls( - self._simulation.state._cpp_sys_def, self.trigger, integrator._cpp_obj, self.beta_ficticious + self._simulation.state._cpp_sys_def, + self.trigger, + integrator._cpp_obj, + self.beta_ficticious, ) - @log(category='sequence', requires_run=True) + @log(category="sequence", requires_run=True) def translate_counts(self): return self._cpp_obj.getCounters(1).translate_counts - @log(category='sequence', requires_run=True) + @log(category="sequence", requires_run=True) def rotate_counts(self): return self._cpp_obj.getCounters(1).rotate_counts - @log(category='scalar', requires_run=True) + @log(category="scalar", requires_run=True) def translate_acceptance_rate(self): acc, rej = self._cpp_obj.getCounters(1).translate_counts if acc + rej == 0: @@ -806,7 +824,7 @@ def translate_acceptance_rate(self): else: return acc / (acc + rej) - @log(category='scalar', requires_run=True) + @log(category="scalar", requires_run=True) def rotate_acceptance_rate(self): acc, rej = self._cpp_obj.getCounters(1).rotate_counts if acc + rej == 0: @@ -814,27 +832,27 @@ def rotate_acceptance_rate(self): else: return acc / (acc + rej) - @log(category='scalar', requires_run=True) + @log(category="scalar", requires_run=True) def average_cluster_size(self): return self._cpp_obj.getCounters(1).average_cluster_size - @log(category='scalar', requires_run=False) + @log(category="scalar", requires_run=False) def maximum_trial_rotation(self): return self._cpp_obj.getMaximumTrialRotation - @log(category='scalar', requires_run=True) + @log(category="scalar", requires_run=True) def total_num_particles_in_clusters(self): return self._cpp_obj.getCounters(1).total_num_particles_in_clusters - @log(category='scalar', requires_run=True) + @log(category="scalar", requires_run=True) def early_abort_counts_no_reverse_link(self): return self._cpp_obj.getCounters(1).early_abort_counts[0] - @log(category='scalar', requires_run=True) + @log(category="scalar", requires_run=True) def early_abort_counts_forced_reverse_link(self): return self._cpp_obj.getCounters(1).early_abort_counts[1] - @log(category='scalar', requires_run=True) + @log(category="scalar", requires_run=True) def average_cluster_size_accepted(self): return self._cpp_obj.getCounters(1).average_cluster_size_accepted