diff --git a/hoomd/HOOMDMath.h b/hoomd/HOOMDMath.h index c2f5cfc890..e68651acec 100644 --- a/hoomd/HOOMDMath.h +++ b/hoomd/HOOMDMath.h @@ -296,6 +296,7 @@ inline HOSTDEVICE void sincos(float x, float& s, float& c) } //! Compute both of sin of x and cos of x with double precision +#ifndef HOOMD_LLVMJIT_BUILD inline HOSTDEVICE void sincos(double x, double& s, double& c) { #if defined(__HIP_DEVICE_COMPILE__) @@ -330,6 +331,7 @@ inline HOSTDEVICE void sincospi(double x, double& s, double& c) fast::sincos(M_PI * x, s, c); #endif } +#endif //! Compute the pow of x,y with single precison via exp(log) refactoring - NOTE: UNDEFINED FOR //! NEGATIVE BASES diff --git a/hoomd/RNGIdentifiers.h b/hoomd/RNGIdentifiers.h index 7ad0393be4..fef76c594f 100644 --- a/hoomd/RNGIdentifiers.h +++ b/hoomd/RNGIdentifiers.h @@ -70,6 +70,7 @@ struct RNGIdentifier static const uint8_t BussiThermostat = 45; static const uint8_t ConstantPressure = 46; static const uint8_t MPCDCellList = 47; + static const uint8_t UpdaterVMMC=48; }; } // namespace hoomd diff --git a/hoomd/hpmc/CMakeLists.txt b/hoomd/hpmc/CMakeLists.txt index 75301d00de..065a0b5d7b 100644 --- a/hoomd/hpmc/CMakeLists.txt +++ b/hoomd/hpmc/CMakeLists.txt @@ -107,6 +107,7 @@ set(_hpmc_headers UpdaterMuVT.h UpdaterQuickCompress.h UpdaterShape.h + UpdaterVirtualMoveMonteCarlo.h XenoCollide2D.h XenoCollide3D.h XenoSweep3D.h diff --git a/hoomd/hpmc/HPMCCounters.h b/hoomd/hpmc/HPMCCounters.h index 571b912f98..6a302a74fb 100644 --- a/hoomd/hpmc/HPMCCounters.h +++ b/hoomd/hpmc/HPMCCounters.h @@ -349,6 +349,129 @@ struct hpmc_clusters_counters_t } }; +struct hpmc_virtual_moves_counters_t + { + unsigned long long int num_aabbtree_builds; + + unsigned long long int translate_accept_count; + unsigned long long int translate_reject_count; + unsigned long long int translate_total_num_particles_in_moved_clusters; + unsigned long long int translate_reject_inactive_region; + unsigned long long int translate_reject_oversized_cluster; + unsigned long long int translate_reject_p_reverse_link_zero; + unsigned long long int translate_reject_q_reverse_link_zero; + + unsigned long long int rotate_accept_count; + unsigned long long int rotate_reject_count; + unsigned long long int rotate_total_num_particles_in_moved_clusters; + unsigned long long int rotate_reject_inactive_region; + unsigned long long int rotate_reject_oversized_cluster; + unsigned long long int rotate_reject_p_reverse_link_zero; + unsigned long long int rotate_reject_q_reverse_link_zero; + + //! Construct a zero set of counters + DEVICE hpmc_virtual_moves_counters_t() + { + num_aabbtree_builds = 0; + + translate_accept_count = 0; + translate_reject_count = 0; + translate_total_num_particles_in_moved_clusters = 0; + translate_reject_inactive_region = 0; + translate_reject_oversized_cluster = 0; + translate_reject_p_reverse_link_zero = 0; + translate_reject_q_reverse_link_zero = 0; + + rotate_accept_count = 0; + rotate_reject_count = 0; + rotate_total_num_particles_in_moved_clusters = 0; + rotate_reject_inactive_region = 0; + rotate_reject_oversized_cluster = 0; + rotate_reject_p_reverse_link_zero = 0; + rotate_reject_q_reverse_link_zero = 0; + } + + std::pair getTranslateCounts() + { + return std::make_pair(translate_accept_count, translate_reject_count); + } + + std::pair getRotateCounts() + { + return std::make_pair(rotate_accept_count, translate_reject_count); + } + + std::array getRotateRejectCounts() + { + std::array counts = {rotate_reject_inactive_region, + rotate_reject_oversized_cluster, + rotate_reject_p_reverse_link_zero, + rotate_reject_q_reverse_link_zero}; + return counts; + } + + std::array getTranslateRejectCounts() + { + std::array counts = {translate_reject_inactive_region, + translate_reject_oversized_cluster, + translate_reject_p_reverse_link_zero, + translate_reject_q_reverse_link_zero}; + return counts; + } + + DEVICE unsigned long long int getTotalNumParticlesInClustersRotate() + { + return rotate_total_num_particles_in_moved_clusters; + } + + DEVICE unsigned long long int getTotalNumParticlesInClustersTranslate() + { + return translate_total_num_particles_in_moved_clusters; + } + + DEVICE unsigned long long int getNumAABBTreeBuilds() + { + return num_aabbtree_builds; + } + }; + +DEVICE inline hpmc_virtual_moves_counters_t operator-(const hpmc_virtual_moves_counters_t& a, + hpmc_virtual_moves_counters_t& b) + { + hpmc_virtual_moves_counters_t result; + result.num_aabbtree_builds = a.num_aabbtree_builds - b.num_aabbtree_builds; + + result.rotate_accept_count = a.rotate_accept_count - b.rotate_accept_count; + result.rotate_reject_count = a.rotate_reject_count - b.rotate_reject_count; + result.rotate_total_num_particles_in_moved_clusters + = a.rotate_total_num_particles_in_moved_clusters + + b.rotate_total_num_particles_in_moved_clusters; + result.rotate_reject_inactive_region + = a.rotate_reject_inactive_region + b.rotate_reject_inactive_region; + result.rotate_reject_oversized_cluster + = a.rotate_reject_oversized_cluster + b.rotate_reject_oversized_cluster; + result.rotate_reject_p_reverse_link_zero + = a.rotate_reject_p_reverse_link_zero + b.rotate_reject_p_reverse_link_zero; + result.rotate_reject_q_reverse_link_zero + = a.rotate_reject_q_reverse_link_zero + b.rotate_reject_q_reverse_link_zero; + + result.translate_accept_count = a.translate_accept_count - b.translate_accept_count; + result.translate_reject_count = a.translate_reject_count - b.translate_reject_count; + result.translate_total_num_particles_in_moved_clusters + = a.translate_total_num_particles_in_moved_clusters + + b.translate_total_num_particles_in_moved_clusters; + result.translate_reject_inactive_region + = a.translate_reject_inactive_region + b.translate_reject_inactive_region; + result.translate_reject_oversized_cluster + = a.translate_reject_oversized_cluster + b.translate_reject_oversized_cluster; + result.translate_reject_p_reverse_link_zero + = a.translate_reject_p_reverse_link_zero + b.translate_reject_p_reverse_link_zero; + result.translate_reject_q_reverse_link_zero + = a.translate_reject_q_reverse_link_zero + b.translate_reject_q_reverse_link_zero; + + return result; + } + DEVICE inline hpmc_muvt_counters_t operator-(const hpmc_muvt_counters_t& a, const hpmc_muvt_counters_t& b) { diff --git a/hoomd/hpmc/IntegratorHPMC.h b/hoomd/hpmc/IntegratorHPMC.h index ed913535c0..7ff1ab680f 100644 --- a/hoomd/hpmc/IntegratorHPMC.h +++ b/hoomd/hpmc/IntegratorHPMC.h @@ -494,6 +494,12 @@ class PYBIND11_EXPORT IntegratorHPMC : public Integrator return m_pair_energy_search_radius; } + Scalar getNominalWidth() + { + return m_nominal_width; + } + + protected: unsigned int m_translation_move_probability; //!< Fraction of moves that are translation moves. unsigned int m_nselect; //!< Number of particles to select for trial moves diff --git a/hoomd/hpmc/IntegratorHPMCMono.h b/hoomd/hpmc/IntegratorHPMCMono.h index d23c84f3d5..1da8d45a3f 100644 --- a/hoomd/hpmc/IntegratorHPMCMono.h +++ b/hoomd/hpmc/IntegratorHPMCMono.h @@ -289,7 +289,7 @@ template class IntegratorHPMCMono : public IntegratorHPMC std::shared_ptr selected_pair = nullptr); //! Build the AABB tree (if needed) - const hoomd::detail::AABBTree& buildAABBTree(); + hoomd::detail::AABBTree& buildAABBTree(); //! Make list of image indices for boxes to check in small-box mode const std::vector>& updateImageList(); @@ -328,6 +328,11 @@ template class IntegratorHPMCMono : public IntegratorHPMC return type_shapes; } + std::vector getShapeCircumsphereRadius() + { + return m_shape_circumsphere_radius; + } + protected: std::vector> m_params; //!< Parameters for each particle type on GPU @@ -847,6 +852,7 @@ template void IntegratorHPMCMono::update(uint64_t timestep) } // end loop over all particles } // end loop over nselect + { ArrayHandle h_postype(m_pdata->getPositions(), access_location::host, @@ -1591,7 +1597,8 @@ template void IntegratorHPMCMono::growAABBList(unsigned int \returns A reference to the tree. */ -template const hoomd::detail::AABBTree& IntegratorHPMCMono::buildAABBTree() +template +hoomd::detail::AABBTree& IntegratorHPMCMono::buildAABBTree() { if (m_aabb_tree_invalid) { diff --git a/hoomd/hpmc/UpdaterVirtualMoveMonteCarlo.h b/hoomd/hpmc/UpdaterVirtualMoveMonteCarlo.h new file mode 100644 index 0000000000..672e41b3b4 --- /dev/null +++ b/hoomd/hpmc/UpdaterVirtualMoveMonteCarlo.h @@ -0,0 +1,1273 @@ +// Copyright (c) 2009-2024 The Regents of the University of Michigan. +// Part of HOOMD-blue, released under the BSD 3-Clause License. + +#ifndef _UPDATER_HPMC_VIRTUAL_CLUSTER_MOVES_ +#define _UPDATER_HPMC_VIRTUAL_CLUSTER_MOVES_ + +#include "hoomd/Updater.h" + +#include "hoomd/RNGIdentifiers.h" +#include "hoomd/RandomNumbers.h" + +#include +#include + +#include "HPMCCounters.h" +#include "IntegratorHPMCMono.h" +#include "Moves.h" +#include + +namespace hoomd + { + +namespace hpmc + { + +struct ClusterData + { + void clear() + { + m_linkers.clear(); + m_linkers_added.clear(); + m_linker_rotation_centers.clear(); + } + + void update_cluster(unsigned int idx, vec3 rotation_center) + { + if (m_linkers_added.find(idx) == m_linkers_added.end()) + { + m_linkers_added.insert(idx); + m_linkers.push_back(idx); + m_linker_rotation_centers.push_back(rotation_center); + } + } + + size_t current_cluster_size() + { + return m_linkers_added.size(); + } + + std::vector m_linkers; + std::set m_linkers_added; + std::vector> m_linker_rotation_centers; + }; + +struct PotentialLinkData + { + void clear() + { + m_potential_intracluster_link_idxs.clear(); + m_potential_intracluster_link_deltaU_fwds.clear(); + } + + void update(unsigned int i, unsigned int j, LongReal delta_u_ij_fwd) + { + m_potential_intracluster_link_idxs.push_back(std::make_pair(i, j)); + m_potential_intracluster_link_deltaU_fwds.push_back(delta_u_ij_fwd); + } + + std::vector> m_potential_intracluster_link_idxs; + std::vector m_potential_intracluster_link_deltaU_fwds; + }; + +/*! Virtual move Monte Carlo Algorithm. + + See Whitelam and Geissler 2007. +*/ + +template class UpdaterVMMC : public Updater + { + public: + //! Constructor + /*! \param sysdef System definition + \param trigger When to run updater + \param mc HPMC integrator + */ + UpdaterVMMC(std::shared_ptr sysdef, + std::shared_ptr trigger, + std::shared_ptr> mc); + + //! Destructor + virtual ~UpdaterVMMC(); + + //! Run the updater + /*! \param timestep timestep at which update is being evaluated + */ + virtual void update(uint64_t timestep); + + /// Set the number of moves per particle per call to update() + void setAttemptsPerParticle(unsigned int attempts_per_particle) + { + m_attempts_per_particle = attempts_per_particle; + } + + unsigned int getAttemptsPerParticle() + { + return m_attempts_per_particle; + } + + void setTranslationMoveProbability(LongReal p) + { + m_translation_move_probability = p; + } + + LongReal getTranslationMoveProbability() + { + return m_translation_move_probability; + } + + void setMaximumTrialRotation(LongReal a) + { + m_maximum_rotate_size = a; + } + + LongReal getMaximumTrialRotation() + { + return m_maximum_rotate_size; + } + + void setMaximumTrialTranslation(LongReal a) + { + m_maximum_translate_size = a; + } + + LongReal getMaximumTrialTranslation() + { + return m_maximum_translate_size; + } + + void setMaximumTrialCenterOfRotationShift(LongReal d) + { + m_maximum_trial_center_of_rotation_shift = d; + } + + LongReal getMaximumTrialCenterOfRotationShift() + { + return m_maximum_trial_center_of_rotation_shift; + } + + unsigned int getMaximumAllowedClusterSize() + { + return m_maximum_allowed_cluster_size; + } + + void setMaximumAllowedClusterSize(unsigned int s) + { + m_maximum_allowed_cluster_size = s; + } + + LongReal getClusterSizeDistributionPrefactor() + { + return m_cluster_size_distribution_prefactor; + } + + void setClusterSizeDistributionPrefactor(LongReal f) + { + m_cluster_size_distribution_prefactor = f; + } + + std::string getClusterSizeLimitMode() + { + return m_cluster_size_limit_mode; + } + + void setClusterSizeLimitMode(std::string s) + { + m_cluster_size_limit_mode = s; + } + + /// Reset statistics counters + virtual void resetStats() + { + m_count_run_start = m_count_total; + } + + LongReal getRebuildTreeBuffer() + { + return m_rebuild_tree_buffer; + } + + void setRebuildTreeBuffer(LongReal d) + { + m_rebuild_tree_buffer = d; + } + + uint16_t getInstance() + { + return m_instance; + } + + void setInstance(uint16_t instance) + { + m_instance = instance; + } + + /*! \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 + */ + hpmc_virtual_moves_counters_t getCounters(unsigned int mode) + { + hpmc_virtual_moves_counters_t result; + if (mode == 0) + result = m_count_total; + else if (mode == 1) + result = m_count_total - m_count_run_start; + else + result = m_count_total - m_count_step_start; + return result; + } + + protected: + std::shared_ptr> m_mc; //!< HPMC integrator + unsigned int m_attempts_per_particle; //!< Number of attempted moves per particle each time + //!< update() is called + 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; + LongReal m_rebuild_tree_buffer; + + detail::UpdateOrder m_update_order; + + /* GPUVector m_postype_backup; //!< Old local positions */ + /* GPUVector m_orientation_backup; //!< Old local orientations */ + /* GPUVector m_image_backup; //!< Old local images */ + LongReal m_maximum_translate_size; + Scalar m_maximum_rotate_size; + LongReal m_maximum_trial_center_of_rotation_shift; + uint16_t m_instance; + + // containers used during creation and execution of cluster moves + std::vector> m_positions_after_move; + ClusterData m_cluster_data; + PotentialLinkData m_potential_link_data; + + hpmc_virtual_moves_counters_t m_count_total; + hpmc_virtual_moves_counters_t m_count_run_start; + hpmc_virtual_moves_counters_t m_count_step_start; + }; + +template +UpdaterVMMC::UpdaterVMMC(std::shared_ptr sysdef, + std::shared_ptr trigger, + std::shared_ptr> mc) + : Updater(sysdef, trigger), m_mc(mc), m_update_order(m_pdata->getN()) + { + m_exec_conf->msg->notice(5) << "Constructing UpdaterVMMC" << std::endl; + + // initialize stats + resetStats(); + + // initialize memory + /* GPUVector(1,this->m_exec_conf).swap(m_postype_backup); */ + /* GPUVector(1,this->m_exec_conf).swap(m_orientation_backup); */ + /* GPUVector(1,this->m_exec_conf).swap(m_image_backup); */ + } + +template UpdaterVMMC::~UpdaterVMMC() + { + m_exec_conf->msg->notice(5) << "Destroying UpdaterVMMC" << std::endl; + } + +/*! Perform a cluster move + \param timestep Current time step of the simulation + + The steps of the algorithm are: + 1. Pick a random seed particle i and a random trial displacement/rotation + 2. Loop over all of i's neighbors (j): + 2a. Add j to the cluster with probability p_ij, which needs u_ij in the current state, and + u_ij after only i has been moved 2b. If j added to cluster, loop over all of its neighbors and + add to cluster with probability p_jk 2c. Continue recursively looping over neighbors until no + candidate cluster members are left + 3. Do the move + 4. Do the reverse move virtually on particle i, and repeat the cluster-making process + - This step gives you all of the information you need for evaluating the acceptance + criterion + 5. Accept the total move according to acceptance criterion. + + So in a coarser-grained level of detail, it's: + 1. Calculate stuff based on current configuration and trial moves to find cluster + 2. Move cluster + 3. Calculate stuff based on config with cluster moved + 4. Accept or reject cluster move based on stuff calculated in steps 1 and 3 +*/ +template void UpdaterVMMC::update(uint64_t timestep) + { + m_exec_conf->msg->notice(4) << "VMMC update() timestep " << timestep << std::endl; + Updater::update(timestep); + + // if no particles, exit early + if (!m_pdata->getN()) + return; + + m_count_step_start = m_count_total; + + // get stuff needed for the calculation + 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(); + ArrayHandle h_charge(m_pdata->getCharges(), access_location::host, access_mode::read); + ArrayHandle h_diameter(m_pdata->getDiameters(), + access_location::host, + access_mode::read); + ArrayHandle h_image(m_pdata->getImages(), access_location::host, access_mode::readwrite); + const BoxDim box = m_pdata->getBox(); + unsigned int ndim = this->m_sysdef->getNDimensions(); + uint16_t user_seed = m_sysdef->getSeed(); + const auto& mc_params = m_mc->getParams(); + unsigned int maximum_allowed_cluster_size + = m_maximum_allowed_cluster_size == 0 ? m_pdata->getN() : m_maximum_allowed_cluster_size; + GPUArray pos_last_tree_build(m_pdata->getN(), m_exec_conf); + ArrayHandle h_pos_last_tree_build(pos_last_tree_build, + access_location::host, + access_mode::readwrite); + { + ArrayHandle h_postype(m_pdata->getPositions(), + access_location::host, + access_mode::readwrite); + for (unsigned int idx = 0; idx < m_pdata->getN(); idx++) + { + h_pos_last_tree_build.data[idx] = h_postype.data[idx]; + } + } + +#ifdef ENABLE_MPI + // compute the width of the active region + Scalar3 npd = box.getNearestPlaneDistance(); + Scalar3 ghost_fraction = m_mc->getNominalWidth() / npd; +#endif + Scalar3 npd_global = m_pdata->getGlobalBox().getNearestPlaneDistance(); + Scalar min_npd = detail::min(npd_global.x, npd_global.y); + if (this->m_sysdef->getNDimensions() == 3) + { + min_npd = detail::min(min_npd, npd_global.z); + } + + if (min_npd < m_rebuild_tree_buffer) + { + m_exec_conf->msg->warning() + << "aabb_tree_buffer too big for simulation box, reducing aabb_tree_buffer to " + << min_npd << "." << std::endl; + } + LongReal min_displacement_rebuild_tree = detail::min(min_npd, m_rebuild_tree_buffer); + + // Shuffle the order of particles for this step + m_update_order.resize(m_pdata->getN()); + m_update_order.shuffle(timestep, m_sysdef->getSeed(), m_exec_conf->getRank()); + + // get the AABB Tree + hoomd::detail::AABBTree& mc_aabb_tree = m_mc->buildAABBTree(); + // update the image list + const std::vector>& image_list = m_mc->updateImageList(); + + // access interaction matrix + ArrayHandle h_overlaps(m_mc->getInteractionMatrix(), + access_location::host, + access_mode::read); + + // counters + hpmc_counters_t counters; + + m_positions_after_move.resize(m_pdata->getN()); + + for (unsigned int i_trial = 0; i_trial < m_attempts_per_particle; i_trial++) + { + for (unsigned int current_seed = 0; current_seed < m_pdata->getN(); current_seed++) + { + m_cluster_data.clear(); + m_potential_link_data.clear(); + LongReal p_link_probability_ratio = 1.0; // 3rd product in eqn 13 from the paper + LongReal p_failed_link_probability_ratio = 1.0; // 1st product in eqn 13 from the paper + LongReal deltaU = 0.0; // change in energy before and after the cluster move + bool skip_to_next_seed = false; + + // seed particle that starts the virtual move + unsigned int seed_idx = m_update_order[current_seed]; + + // 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)); + 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); + quat virtual_rotate_move(1.0, vec3(0, 0, 0)); + if (move_type_translate) + { + move_translate(virtual_translate_move, rng_i, m_maximum_translate_size, ndim); + } + else + { + move_translate(virtual_translate_move, + rng_i, + m_maximum_trial_center_of_rotation_shift, + ndim); + if (ndim == 2) + { + move_rotate<2>(virtual_rotate_move, rng_i, m_maximum_rotate_size); + } + else + { + move_rotate<3>(virtual_rotate_move, rng_i, m_maximum_rotate_size); + } + } + + // add linker and current image rotation center to their sets + Scalar4 postype_seed; + { + // scope code that uses h_postype and h_orientation bc mc_aabb_tree.buildAABBTree() + // needs access to those arrays + ArrayHandle h_postype(m_pdata->getPositions(), + access_location::host, + access_mode::readwrite); + ArrayHandle h_orientation(m_pdata->getOrientationArray(), + access_location::host, + access_mode::readwrite); + postype_seed = h_postype.data[seed_idx]; + vec3 pos_seed = vec3(postype_seed); + m_cluster_data.update_cluster(seed_idx, virtual_translate_move + pos_seed); + + m_exec_conf->msg->notice(6) + << " VMMC virtual translation = (" << virtual_translate_move.x << ", " + << virtual_translate_move.y << ", " << virtual_translate_move.z << ")" + << std::endl; + + // loop over neighbors of cluster members to find new cluster members + m_exec_conf->msg->notice(5) << "VMMC seed tag: " << seed_idx << std::endl; + /* int3 _images = make_int3(0, 0, 0); */ + while (m_cluster_data.m_linkers.size() > 0) + { + if (skip_to_next_seed) + { + break; + } + unsigned int i_linker = m_cluster_data.m_linkers[0]; + m_cluster_data.m_linkers.erase(m_cluster_data.m_linkers.begin()); + vec3 center_of_rotation = m_cluster_data.m_linker_rotation_centers[0]; + m_cluster_data.m_linker_rotation_centers.erase( + m_cluster_data.m_linker_rotation_centers.begin()); + m_exec_conf->msg->notice(5) << "VMMC linker tag: " << i_linker << std::endl; + + // get position and orientation of linker + Scalar4 postype_linker = h_postype.data[i_linker]; + quat orientation_linker = quat(h_orientation.data[i_linker]); + vec3 pos_linker = vec3(postype_linker); + vec3 pos_linker_after_move_primary_image; + vec3 pos_linker_after_reverse_move_primary_image; + if (move_type_translate) + { + pos_linker_after_move_primary_image = pos_linker + virtual_translate_move; + pos_linker_after_reverse_move_primary_image + = pos_linker - virtual_translate_move; + } + else + { + 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; + } + m_positions_after_move[i_linker] = pos_linker_after_move_primary_image; + + int type_linker = __scalar_as_int(postype_linker.w); + Shape shape_linker(orientation_linker, mc_params[type_linker]); + Shape shape_linker_after_move(virtual_rotate_move * orientation_linker, + mc_params[type_linker]); + Shape shape_linker_after_reverse_move(conj(virtual_rotate_move) + * orientation_linker, + mc_params[type_linker]); + +// linker must be an active particle before and after the move +#ifdef ENABLE_MPI + if (m_sysdef->isDomainDecomposed()) + { + if (!isActive(make_scalar3(pos_linker.x, pos_linker.y, pos_linker.z), + box, + ghost_fraction) + || (!isActive(make_scalar3(pos_linker_after_move_primary_image.x, + pos_linker_after_move_primary_image.y, + pos_linker_after_move_primary_image.z), + box, + ghost_fraction)) + + ) + { + if (move_type_translate) + { + m_count_total.translate_reject_count++; + m_count_total.translate_reject_inactive_region++; + } + else + { + m_count_total.rotate_reject_count++; + m_count_total.rotate_reject_inactive_region++; + } + skip_to_next_seed = true; + break; + } + } +#endif + + // search for all particles that might touch this one + LongReal R_query = m_mc->getShapeCircumsphereRadius()[type_linker]; + + if (m_mc->hasPairInteractions()) + { + // Extend the search to include the pair interaction r_cut + // subtract minimum AABB extent from search radius + R_query + = std::max(R_query, + pair_energy_search_radius[type_linker] - min_core_radius); + } + hoomd::detail::AABB aabb_linker_local_no_moves + = hoomd::detail::AABB(vec3(0, 0, 0), R_query); + hoomd::detail::AABB aabb_linker_local_forward + = hoomd::detail::AABB(pos_linker_after_move_primary_image - pos_linker, + R_query); + hoomd::detail::AABB aabb_linker_local_reverse = hoomd::detail::AABB( + pos_linker_after_reverse_move_primary_image - pos_linker, + R_query); + hoomd::detail::AABB _aabb = hoomd::detail::merge(aabb_linker_local_no_moves, + aabb_linker_local_forward); + hoomd::detail::AABB aabb_linker_local + = hoomd::detail::merge(_aabb, aabb_linker_local_reverse); + + // loop over all images + const unsigned int n_images = (unsigned int)image_list.size(); + for (unsigned int current_image = 0; current_image < n_images; current_image++) + { + if (skip_to_next_seed) + { + break; + } + // create an AABB centered on the linker particle in the current box image + hoomd::detail::AABB aabb = aabb_linker_local; + vec3 pos_linker_image = pos_linker + image_list[current_image]; + vec3 center_of_rotation_image + = center_of_rotation + image_list[current_image]; + aabb.translate(pos_linker_image); + + // stackless search + for (unsigned int current_node_index = 0; + current_node_index < mc_aabb_tree.getNumNodes(); + current_node_index++) + { + if (skip_to_next_seed) + { + break; + } + if (aabb.overlaps(mc_aabb_tree.getNodeAABB(current_node_index))) + { + if (mc_aabb_tree.isNodeLeaf(current_node_index)) + { + for (unsigned int current_linkee = 0; + current_linkee + < mc_aabb_tree.getNodeNumParticles(current_node_index); + current_linkee++) + { + if (skip_to_next_seed) + { + break; + } + m_exec_conf->msg->notice(6) + << "linker_set_size = " + << m_cluster_data.m_linkers.size() << std::endl; + unsigned int j_linkee + = mc_aabb_tree.getNodeParticle(current_node_index, + current_linkee); + if (!(m_cluster_data.m_linkers_added.find(j_linkee) + == m_cluster_data.m_linkers_added.end())) + { + // j already in cluster + continue; + } + Scalar4 postype_linkee = h_postype.data[j_linkee]; + quat orientation_linkee + = quat(h_orientation.data[j_linkee]); + + unsigned int type_linkee + = __scalar_as_int(postype_linkee.w); + Shape shape_linkee(orientation_linkee, + mc_params[type_linkee]); + LongReal max_overlap_distance + = m_mc->getShapeCircumsphereRadius()[type_linker] + + m_mc->getShapeCircumsphereRadius()[type_linkee]; + + // at this point we test for a link between the linker and + // linkee add j to the cluster with probability p_ij = 1 - + // exp(E_C - E_I) (hard overlap makes E_I go to infinity and + // p_ij = 1) so we need E_C (current pair energy) and E_I + // (pair energy after applying virtual move to i) we only + // need E_C if E_I, so we can check for overlaps first + + // pair separation in current real (i.e. non-virtual) + // configuration + vec3 r_linker_linkee + = vec3(postype_linkee) - pos_linker_image; + LongReal r_squared = dot(r_linker_linkee, r_linker_linkee); + + // pair separation after applying virtual move to linker (i) + vec3 r_linker_linkee_after_forward_move_of_linker + = vec3(postype_linkee) + - (pos_linker_after_move_primary_image + + image_list[current_image]); + LongReal r_squared_after_forward_move_of_linker + = dot(r_linker_linkee_after_forward_move_of_linker, + r_linker_linkee_after_forward_move_of_linker); + + // pair separation after applying reverse virtual move to + // linker (i) + vec3 r_linker_linkee_after_reverse_move_of_linker + = vec3(postype_linkee) + - (pos_linker_after_reverse_move_primary_image + + image_list[current_image]); + LongReal r_squared_after_reverse_move_of_linker + = dot(r_linker_linkee_after_reverse_move_of_linker, + r_linker_linkee_after_reverse_move_of_linker); + + m_exec_conf->msg->notice(5) + << "Test for links: pair " << i_linker << " and " + << j_linkee << std::endl; + m_exec_conf->msg->notice(5) + << " Checking for overlap after moving linker between " + "pair " + << i_linker << " and " << j_linkee << std::endl; + bool overlap_after_forward_move_of_linker + = h_overlaps + .data[m_mc->getOverlapIndexer()(type_linker, + type_linkee)] + && r_squared_after_forward_move_of_linker + < max_overlap_distance * max_overlap_distance + && test_overlap( + r_linker_linkee_after_forward_move_of_linker, + shape_linker_after_move, + shape_linkee, + counters.overlap_err_count); + if (overlap_after_forward_move_of_linker) + { + m_exec_conf->msg->notice(5) + << " Overlap found between " << i_linker + << " and " << j_linkee << " -> link formed!" + << std::endl; + + // add linkee to cluster if not in there already and + // after checking maximum cluster_size + if (m_cluster_size_limit_mode == "deterministic" + && m_cluster_data.current_cluster_size() + == maximum_allowed_cluster_size) + { + // abort the move. + if (move_type_translate) + { + m_count_total.translate_reject_count++; + m_count_total + .translate_reject_oversized_cluster++; + } + else + { + m_count_total.rotate_reject_count++; + m_count_total.rotate_reject_oversized_cluster++; + } + skip_to_next_seed = true; + break; + } + else if (m_cluster_size_limit_mode == "probabilistic" + && hoomd::UniformDistribution()(rng_i) + > m_cluster_size_distribution_prefactor + / ((LongReal)m_cluster_data + .current_cluster_size() + + 1)) + { + // abort the move. + if (move_type_translate) + { + m_count_total.translate_reject_count++; + m_count_total + .translate_reject_oversized_cluster++; + } + else + { + m_count_total.rotate_reject_count++; + m_count_total.rotate_reject_oversized_cluster++; + } + skip_to_next_seed = true; + break; + } + else + { + m_cluster_data.update_cluster( + j_linkee, + center_of_rotation_image); + } + + // account for probability of forming this link: + // denominator of p_link_probability_ratio *= 1 + p_link_probability_ratio /= 1.0; + + // 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 + Scalar p_ij_reverse; + bool overlap_after_reverse_move_of_linker + = h_overlaps + .data[m_mc->getOverlapIndexer()(type_linker, + type_linkee)] + && r_squared_after_reverse_move_of_linker + < max_overlap_distance + * max_overlap_distance + && test_overlap( + r_linker_linkee_after_reverse_move_of_linker, + shape_linker_after_reverse_move, + shape_linkee, + counters.overlap_err_count); + if (overlap_after_reverse_move_of_linker) + { + p_ij_reverse = 1.0; + m_exec_conf->msg->notice(5) + << " p_ij_reverse = " << p_ij_reverse + << " by overlap on reverse move of linker." + << std::endl; + } + else + { + double u_ij = m_mc->computeOnePairEnergy( + r_squared, + r_linker_linkee, + type_linker, + shape_linker.orientation, + h_diameter.data[i_linker], + h_charge.data[i_linker], + type_linkee, + shape_linkee.orientation, + h_diameter.data[j_linkee], + h_charge.data[j_linkee]); + double u_ij_after_reverse_move_of_linker + = m_mc->computeOnePairEnergy( + r_squared_after_reverse_move_of_linker, + r_linker_linkee_after_reverse_move_of_linker, + type_linker, + shape_linker_after_reverse_move.orientation, + h_diameter.data[i_linker], + h_charge.data[i_linker], + type_linkee, + shape_linkee.orientation, + h_diameter.data[j_linkee], + h_charge.data[j_linkee]); + 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) + { + if (move_type_translate) + { + m_count_total.translate_reject_count++; + m_count_total + .translate_reject_p_reverse_link_zero++; + } + else + { + m_count_total.rotate_reject_count++; + m_count_total + .rotate_reject_p_reverse_link_zero++; + } + skip_to_next_seed = true; + break; // break loop over linkees + } + p_link_probability_ratio *= p_ij_reverse; + continue; // to next linkee + } // end if (overlap_after_forward_move_of_linker) + m_exec_conf->msg->notice(5) + << " No overlap after moving linker between " + << i_linker << " and " << j_linkee << std::endl; + + // if no overlap, we have to calculate u_ij and + // u_ij_after_move to determine p_ij + double u_ij + = m_mc->computeOnePairEnergy(r_squared, + r_linker_linkee, + type_linker, + shape_linker.orientation, + h_diameter.data[i_linker], + h_charge.data[i_linker], + type_linkee, + shape_linkee.orientation, + h_diameter.data[j_linkee], + h_charge.data[j_linkee]); + double u_ij_after_move = m_mc->computeOnePairEnergy( + r_squared_after_forward_move_of_linker, + r_linker_linkee_after_forward_move_of_linker, + type_linker, + shape_linker_after_move.orientation, + h_diameter.data[i_linker], + h_charge.data[i_linker], + type_linkee, + shape_linkee.orientation, + h_diameter.data[j_linkee], + h_charge.data[j_linkee]); + if (u_ij_after_move == 0.0 && u_ij == 0.0) + { + // non-interacting -> non-interacting; continue to next + // linkee + continue; // to next linkee + } + + LongReal p_ij_forward; + 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) + { + m_exec_conf->msg->notice(5) + << " p_ij_forward = " << p_ij_forward << " > " << r + << " -> link_formed!" << std::endl; + if (m_cluster_size_limit_mode == "deterministic" + && m_cluster_data.current_cluster_size() + == maximum_allowed_cluster_size) + { + // abort the move. + if (move_type_translate) + { + m_count_total.translate_reject_count++; + m_count_total + .translate_reject_oversized_cluster++; + } + else + { + m_count_total.rotate_reject_count++; + m_count_total.rotate_reject_oversized_cluster++; + } + skip_to_next_seed = true; + break; + } + else if (m_cluster_size_limit_mode == "probabilistic" + && hoomd::UniformDistribution()(rng_i) + > m_cluster_size_distribution_prefactor + / ((LongReal)m_cluster_data + .current_cluster_size() + + 1)) + { + // abort the move. + if (move_type_translate) + { + m_count_total.translate_reject_count++; + m_count_total + .translate_reject_oversized_cluster++; + } + else + { + m_count_total.rotate_reject_count++; + m_count_total.rotate_reject_oversized_cluster++; + } + skip_to_next_seed = true; + break; + } + else + { + m_cluster_data.update_cluster( + j_linkee, + center_of_rotation_image); + } + // numerator of p_link_probability_ratio gets multiplied + // by p_ij_reverse + Scalar p_ij_reverse; + bool overlap_after_reverse_move + = h_overlaps + .data[m_mc->getOverlapIndexer()(type_linker, + type_linkee)] + && r_squared_after_reverse_move_of_linker + < max_overlap_distance + * max_overlap_distance + && test_overlap( + r_linker_linkee_after_reverse_move_of_linker, + shape_linker_after_reverse_move, + shape_linkee, + counters.overlap_err_count); + if (overlap_after_reverse_move) + { + p_ij_reverse = 1.0; + m_exec_conf->msg->notice(5) + << " Overlap found on reverse move -> " + "p_ij_reverse = " + << p_ij_reverse << std::endl; + } + else + { + double u_ij_after_reverse_move + = m_mc->computeOnePairEnergy( + r_squared_after_reverse_move_of_linker, + r_linker_linkee_after_reverse_move_of_linker, + type_linker, + shape_linker_after_reverse_move.orientation, + h_diameter.data[i_linker], + h_charge.data[i_linker], + type_linkee, + shape_linkee.orientation, + 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(-(u_ij_after_reverse_move - u_ij) + / kT)); + if (p_ij_reverse == 0.0) + { + m_exec_conf->msg->notice(5) + << " p_ij_reverse = 0; aborting move." + << std::endl; + if (move_type_translate) + { + m_count_total.translate_reject_count++; + m_count_total + .translate_reject_p_reverse_link_zero++; + } + else + { + m_count_total.rotate_reject_count++; + m_count_total + .rotate_reject_p_reverse_link_zero++; + } + skip_to_next_seed = true; + break; + } + m_exec_conf->msg->notice(5) + << " p_ij_reverse = " << p_ij_reverse + << std::endl; + } + p_link_probability_ratio *= p_ij_reverse / p_ij_forward; + continue; + } + else // no link formed + { + m_exec_conf->msg->notice(5) + << " p_ij_forward = " << p_ij_forward + << "; random number = " << r << " -> no link_formed" + << std::endl; + // if j ultimately ends up in the cluster, the dU + // contribution is 0 if j ultimately does not end up in + // the cluster, we do have to account for the energy + // change, so dU_forward += dU either way, we still must + // account for the probability of not forming that link, + // so the denominator of p_failed_link_probability_ratio + // gets multiplied by 1 - p_ij_forward = 1 - (1 - + // exp(-beta * (beta_u_ij_after_move - beta_u_ij))) = + // exp(-beta * (beta_u_ij_after_move - beta_u_ij)) and + // the numerator of p_failed_link_probability_ratio gets + // multiplied by 1 - p_ij_reverse 1 - p_ij_reverse = 1 - + // (1 + // - exp(-beta * (beta_u_ij - beta_u_ij_after_move))) = + // exp(-beta * (beta_u_ij - beta_u_ij_after_move)) which + // means no matter what we only need beta_u_ij (already + // have) and beta_u_ij_after_move (also already have) + // intuitively it makes sense that we don't need + // beta_u_ij_after_reverse move because the two 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 + + // 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(-(u_ij - u_ij_after_move) / kT)); + if (q_ij_reverse == 0.0) + { + m_exec_conf->msg->notice(5) + << " q_ij_reverse = 0; aborting move." + << std::endl; + skip_to_next_seed = true; + if (move_type_translate) + { + m_count_total.translate_reject_count++; + m_count_total + .translate_reject_q_reverse_link_zero++; + } + else + { + m_count_total.rotate_reject_count++; + m_count_total + .rotate_reject_q_reverse_link_zero++; + } + break; + } + m_exec_conf->msg->notice(5) + << " q_ij_forward = " << q_ij_forward << std::endl; + m_exec_conf->msg->notice(5) + << " q_ij_reverse = " << q_ij_reverse << std::endl; + p_failed_link_probability_ratio + *= q_ij_reverse / q_ij_forward; + + // add to possible contributions to deltaU for the + // cluster move + m_potential_link_data.update(i_linker, + j_linkee, + (LongReal)u_ij_after_move + - (LongReal)u_ij); + continue; + } + } // end loop over linkees + } // end if (mc_aabb_tree.isNodeLeaf(current_node_index)) + } // end if + // (aabb.overlaps(mc_aabb_tree.getNodeAABB(current_node_index))) + else + { + // skip ahead + current_node_index += mc_aabb_tree.getNodeSkip(current_node_index); + } + } // end loop over nodes in AABBTree + } // end loop over images + } // end loop over linkers + if (skip_to_next_seed) + { + continue; + } + } + + // find which of the potential cluster members ended up in the cluster + m_exec_conf->msg->notice(5) << "Finding failed intracluster links" << std::endl; + for (unsigned int potential_link_idx = 0; + potential_link_idx + < m_potential_link_data.m_potential_intracluster_link_idxs.size(); + potential_link_idx++) + { + unsigned int linker_i + = m_potential_link_data.m_potential_intracluster_link_idxs[potential_link_idx] + .first; // recall: linker_i in cluster by construction + unsigned int linkee_j + = m_potential_link_data.m_potential_intracluster_link_idxs[potential_link_idx] + .second; + m_exec_conf->msg->notice(10) << " Checking failed, intracluster links between " + << linker_i << " and " << linkee_j << std::endl; + if (m_cluster_data.m_linkers_added.find(linkee_j) + == m_cluster_data.m_linkers_added.end()) + { + // j ultimately did not join the cluster, so there will in general be a + // contribution to deltaU for this i-j pair we have already accounted for the + // probability of not forming the i-j link, so all we do is add to deltaU + deltaU += m_potential_link_data + .m_potential_intracluster_link_deltaU_fwds[potential_link_idx]; + } + else + { + m_exec_conf->msg->notice(5) << " particle with tag " << linkee_j + << " eventually joined the cluster." << std::endl; + // j ultimately joined the cluster so there will be no contribution to deltaU + // for this i-j pair we have also already accounted for the probability of not + // forming that link in the forward and reverse moves so there is nothing to do + // here except for explicitly add 0.0 to deltaU for the sake of code readability + deltaU += 0.0; + } + } + + // test for acceptance + m_exec_conf->msg->notice(5) + << " failed_link_probability_ratio = " << p_failed_link_probability_ratio + << std::endl; + m_exec_conf->msg->notice(5) + << " p_link_probability_ratio = " << p_link_probability_ratio << std::endl; + 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; + 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; + size_t cluster_size = m_cluster_data.m_linkers_added.size(); + if (accept_cluster_move) + { + m_exec_conf->msg->notice(4) << " VMMC move accepted, moving " << cluster_size + << " particles as a cluster" << std::endl; + if (move_type_translate) + { + m_count_total.translate_accept_count++; + m_count_total.translate_total_num_particles_in_moved_clusters + += (unsigned long long int)cluster_size; + } + else + { + m_count_total.rotate_accept_count++; + m_count_total.rotate_total_num_particles_in_moved_clusters + += (unsigned long long int)cluster_size; + } + + // move the particles that are in the cluster + bool rebuild_tree = false; + { + ArrayHandle h_postype(m_pdata->getPositions(), + access_location::host, + access_mode::readwrite); + ArrayHandle h_orientation(m_pdata->getOrientationArray(), + access_location::host, + access_mode::readwrite); + for (unsigned int idx = 0; idx < m_pdata->getN(); idx++) + { + if (!(m_cluster_data.m_linkers_added.find(idx) + == m_cluster_data.m_linkers_added.end())) + { + // particle idx is in the cluster + Scalar4 postype_idx = h_postype.data[idx]; + vec3 new_pos = m_positions_after_move[idx]; + h_postype.data[idx] + = make_scalar4(new_pos.x, new_pos.y, new_pos.z, postype_idx.w); + vec3 displacement + = vec3(h_postype.data[idx]) + - vec3(h_pos_last_tree_build.data[idx]); + LongReal displacement_sq = dot(displacement, displacement); + if (displacement_sq + > min_displacement_rebuild_tree * min_displacement_rebuild_tree) + { + rebuild_tree = true; + } + quat orientation(h_orientation.data[idx]); + h_orientation.data[idx] + = quat_to_scalar4(virtual_rotate_move * orientation); + } + } + if (rebuild_tree) + { + m_mc->invalidateAABBTree(); + // update pos_last_tree_build + for (unsigned int idx = 0; idx < m_pdata->getN(); idx++) + { + box.wrap(h_postype.data[idx], h_image.data[idx]); + h_pos_last_tree_build.data[idx] = h_postype.data[idx]; + } + m_count_total.num_aabbtree_builds++; + } + else + { + for (unsigned int idx = 0; idx < m_pdata->getN(); idx++) + { + // TODO: only loop over the m_cluster_data elements + if (!(m_cluster_data.m_linkers_added.find(idx) + == m_cluster_data.m_linkers_added.end())) + { + LongReal R_query + = m_mc->getShapeCircumsphereRadius()[(size_t)h_postype.data[idx] + .w]; + if (m_mc->hasPairInteractions()) + { + // Extend the search to include the pair interaction r_cut + // subtract minimum AABB extent from search radius + R_query = std::max( + R_query, + pair_energy_search_radius[(size_t)h_postype.data[idx].w] + - min_core_radius); + } + hoomd::detail::AABB new_aabb + = hoomd::detail::AABB(vec3(h_postype.data[idx]), + R_query); + mc_aabb_tree.update(idx, new_aabb); + } + } + } + } + } + else + { + m_exec_conf->msg->notice(4) << "VMMC move rejected (p_acc = " << p_acc + << " < r = " << r << ")" << std::endl; + if (move_type_translate) + { + m_count_total.translate_reject_count++; + } + else + { + m_count_total.rotate_reject_count++; + } + } + m_exec_conf->msg->notice(5) << std::endl; + m_mc->buildAABBTree(); + } // end loop over seed particles + } // end loop over trials + m_mc->invalidateAABBTree(); + } // end update() + +namespace detail + { + +template +void export_UpdaterVirtualMoveMonteCarlo(pybind11::module& m, const std::string& name) + { + pybind11::class_, Updater, std::shared_ptr>>(m, + name.c_str()) + .def(pybind11::init, + std::shared_ptr, + std::shared_ptr>>()) + .def("getCounters", &UpdaterVMMC::getCounters) + .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) + .def_property("attempts_per_particle", + &UpdaterVMMC::getAttemptsPerParticle, + &UpdaterVMMC::setAttemptsPerParticle) + .def_property("maximum_trial_center_of_rotation_shift", + &UpdaterVMMC::getMaximumTrialCenterOfRotationShift, + &UpdaterVMMC::setMaximumTrialCenterOfRotationShift) + .def_property("maximum_allowed_cluster_size", + &UpdaterVMMC::getMaximumAllowedClusterSize, + &UpdaterVMMC::setMaximumAllowedClusterSize) + .def_property("cluster_size_limit_mode", + &UpdaterVMMC::getClusterSizeLimitMode, + &UpdaterVMMC::setClusterSizeLimitMode) + .def_property("cluster_size_distribution_prefactor", + &UpdaterVMMC::getClusterSizeDistributionPrefactor, + &UpdaterVMMC::setClusterSizeDistributionPrefactor) + .def_property("aabb_tree_buffer", + &UpdaterVMMC::getRebuildTreeBuffer, + &UpdaterVMMC::setRebuildTreeBuffer) + .def_property("instance", + &UpdaterVMMC::getInstance, + &UpdaterVMMC::setInstance); + } + +inline void export_hpmc_virtual_moves_counters(pybind11::module& m) + { + pybind11::class_(m, "hpmc_virtual_moves_counters_t") + .def_property_readonly("translate_counts", + &hpmc_virtual_moves_counters_t::getTranslateCounts) + .def_property_readonly("translate_rejection_counts", + &hpmc_virtual_moves_counters_t::getTranslateRejectCounts) + .def_property_readonly( + "translate_num_particles_in_moved_clusters", + &hpmc_virtual_moves_counters_t::getTotalNumParticlesInClustersTranslate) + .def_property_readonly("rotate_counts", &hpmc_virtual_moves_counters_t::getRotateCounts) + .def_property_readonly("rotate_rejection_counts", + &hpmc_virtual_moves_counters_t::getRotateRejectCounts) + .def_property_readonly( + "rotate_num_particles_in_moved_clusters", + &hpmc_virtual_moves_counters_t::getTotalNumParticlesInClustersRotate) + .def_property_readonly("num_aabbtree_builds", &hpmc_virtual_moves_counters_t::getNumAABBTreeBuilds); + } + + } // end namespace detail + } // end namespace hpmc + // + } // end namespace hoomd + +#endif // _UPDATER_HPMC_VIRTUAL_CLUSTER_MOVES_ diff --git a/hoomd/hpmc/module.cc b/hoomd/hpmc/module.cc index 255c31cf4c..b3548ee928 100644 --- a/hoomd/hpmc/module.cc +++ b/hoomd/hpmc/module.cc @@ -23,6 +23,7 @@ #include "UpdaterGCA.h" #include "UpdaterMuVT.h" #include "UpdaterQuickCompress.h" +#include "UpdaterVirtualMoveMonteCarlo.h" #include "GPUTree.h" @@ -151,6 +152,7 @@ PYBIND11_MODULE(_hpmc, m) export_hpmc_muvt_counters(m); export_hpmc_clusters_counters(m); + export_hpmc_virtual_moves_counters(m); export_hpmc_nec_counters(m); diff --git a/hoomd/hpmc/module_convex_polyhedron.cc b/hoomd/hpmc/module_convex_polyhedron.cc index 888d5311b0..40225f9398 100644 --- a/hoomd/hpmc/module_convex_polyhedron.cc +++ b/hoomd/hpmc/module_convex_polyhedron.cc @@ -13,6 +13,7 @@ #include "UpdaterGCA.h" #include "UpdaterMuVT.h" +#include "UpdaterVirtualMoveMonteCarlo.h" #include "ShapeMoves.h" #include "ShapeUtils.h" @@ -38,6 +39,7 @@ void export_convex_polyhedron(pybind11::module& m) export_ComputeFreeVolume(m, "ComputeFreeVolumeConvexPolyhedron"); export_ComputeSDF(m, "ComputeSDFConvexPolyhedron"); export_UpdaterMuVT(m, "UpdaterMuVTConvexPolyhedron"); + export_UpdaterVirtualMoveMonteCarlo(m, "UpdaterVMMCConvexPolyhedron"); export_UpdaterGCA(m, "UpdaterGCAConvexPolyhedron"); export_MassProperties(m, "MassPropertiesConvexPolyhedron"); diff --git a/hoomd/hpmc/module_convex_spheropolyhedron.cc b/hoomd/hpmc/module_convex_spheropolyhedron.cc index ad1765f0de..cd712d6dca 100644 --- a/hoomd/hpmc/module_convex_spheropolyhedron.cc +++ b/hoomd/hpmc/module_convex_spheropolyhedron.cc @@ -12,6 +12,7 @@ #include "UpdaterGCA.h" #include "UpdaterMuVT.h" +#include "UpdaterVirtualMoveMonteCarlo.h" #include "ShapeMoves.h" #include "UpdaterShape.h" @@ -35,6 +36,7 @@ void export_convex_spheropolyhedron(pybind11::module& m) export_ComputeFreeVolume(m, "ComputeFreeVolumeSpheropolyhedron"); export_ComputeSDF(m, "ComputeSDFConvexSpheropolyhedron"); export_UpdaterMuVT(m, "UpdaterMuVTConvexSpheropolyhedron"); + export_UpdaterVirtualMoveMonteCarlo(m, "UpdaterVMMCConvexSpheropolyhedron"); export_UpdaterGCA(m, "UpdaterGCAConvexSpheropolyhedron"); export_MassProperties(m, "MassPropertiesConvexSpheropolyhedron"); diff --git a/hoomd/hpmc/module_sphere.cc b/hoomd/hpmc/module_sphere.cc index 4bf48671a1..85fe4bfcf3 100644 --- a/hoomd/hpmc/module_sphere.cc +++ b/hoomd/hpmc/module_sphere.cc @@ -13,6 +13,7 @@ #include "UpdaterGCA.h" #include "UpdaterMuVT.h" +#include "UpdaterVirtualMoveMonteCarlo.h" #ifdef ENABLE_HIP #include "ComputeFreeVolumeGPU.h" @@ -34,6 +35,7 @@ void export_sphere(pybind11::module& m) export_ComputeFreeVolume(m, "ComputeFreeVolumeSphere"); export_ComputeSDF(m, "ComputeSDFSphere"); export_UpdaterMuVT(m, "UpdaterMuVTSphere"); + export_UpdaterVirtualMoveMonteCarlo(m, "UpdaterVMMCSphere"); export_UpdaterGCA(m, "UpdaterGCASphere"); export_ExternalFieldWall(m, "WallSphere"); diff --git a/hoomd/hpmc/module_union_sphere.cc b/hoomd/hpmc/module_union_sphere.cc index 946cde985f..f8e9087b8d 100644 --- a/hoomd/hpmc/module_union_sphere.cc +++ b/hoomd/hpmc/module_union_sphere.cc @@ -12,6 +12,7 @@ #include "UpdaterGCA.h" #include "UpdaterMuVT.h" +#include "UpdaterVirtualMoveMonteCarlo.h" #ifdef ENABLE_HIP #include "ComputeFreeVolumeGPU.h" @@ -32,6 +33,7 @@ void export_union_sphere(pybind11::module& m) export_ComputeFreeVolume>(m, "ComputeFreeVolumeSphereUnion"); export_ComputeSDF>(m, "ComputeSDFSphereUnion"); export_UpdaterMuVT>(m, "UpdaterMuVTSphereUnion"); + export_UpdaterVirtualMoveMonteCarlo>(m, "UpdaterVMMCSphereUnion"); export_UpdaterGCA>(m, "UpdaterGCASphereUnion"); export_ExternalFieldWall>(m, "WallSphereUnion"); diff --git a/hoomd/hpmc/update.py b/hoomd/hpmc/update.py index de46d59d58..99de9c12ff 100644 --- a/hoomd/hpmc/update.py +++ b/hoomd/hpmc/update.py @@ -675,6 +675,223 @@ def shape_move_energy(self): return self._cpp_obj.getShapeMoveEnergy(self._simulation.timestep) +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. + + 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} + + ---------- + + **Members defined in** `VirtualClusterMoves`: + + Attributes: + instance (int): + When using multiple `QuickCompress` updaters in a single simulation, + give each a unique value for `instance` so that they generate + different streams of random numbers. + + """ + + __doc__ = __doc__.replace("{inherited}", Updater._doc_inherited) + + def __init__( + self, + aabb_tree_buffer, + 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, + instance=0, + ): + super().__init__(trigger) + if maximum_allowed_cluster_size is None: + # set no limits on maximum allowed cluster size + maximum_allowed_cluster_size = 0 + if cluster_size_limit_mode is None: + cluster_size_limit_mode = "" + param_dict = ParameterDict( + aabb_tree_buffer=float(aabb_tree_buffer), + 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), + instance=int(instance), + ) + self._param_dict.update(param_dict) + + def _attach_hook(self): + self._simulation._warn_if_seed_unset() + integrator = self._simulation.operations.integrator + if not isinstance(integrator, integrate.HPMCIntegrator): + raise RuntimeError("The integrator must be a HPMC integrator.") + cpp_cls_name = "UpdaterVMMC" + 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, + ) + + @log(category="sequence", requires_run=True) + def translate_counts(self): + """tuple[int, int]: Count of the accepted and rejected translate moves. + + Note: + The counts are reset to 0 at the start of each `hoomd.Simulation.run`. + """ + return self._cpp_obj.getCounters(1).translate_counts + + @log(category="sequence", requires_run=True) + def translate_rejection_counts(self): + """tuple[int, int, int, int]: Counts of different types of rejected translate moves. + + The items in the list correspond to the following reasons for a rejection: + - a particle in the cluster starts or ends in an inactive region + - the cluster becomes too large + - the probability of forming one of the links in the reverse move was zero + - the probability of not forming one of the non-links in the reverse move was zero + + Note: + The counts are reset to 0 at the start of each `hoomd.Simulation.run`. + + """ + return self._cpp_obj.getCounters(1).translate_rejection_counts + + @log(category="sequence", requires_run=True) + def translate_num_particles_in_moved_clusters(self): + """int: Total number of particles in clusters with accepted translate moves. + + Note: + The counts are reset to 0 at the start of each `hoomd.Simulation.run`. + + """ + return self._cpp_obj.getCounters(1).translate_num_particles_in_moved_clusters + + @log(category="sequence", requires_run=True) + def rotate_counts(self): + """tuple[int, int]: Count of the accepted and rejected rotate moves. + + Note: + The counts are reset to 0 at the start of each `hoomd.Simulation.run`. + """ + return self._cpp_obj.getCounters(1).rotate_counts + + @log(category="sequence", requires_run=True) + def rotate_rejection_counts(self): + """tuple[int, int, int, int]: Counts of different types of rejected rotate moves. + + The items in the list correspond to the following reasons for a rejection: + - a particle in the cluster starts or ends in an inactive region + - the cluster becomes too large + - the probability of forming one of the links in the reverse move was zero + - the probability of not forming one of the non-links in the reverse move was zero + + Note: + The counts are reset to 0 at the start of each `hoomd.Simulation.run`. + + """ + return self._cpp_obj.getCounters(1).rotate_rejection_counts + + @log(category="scalar", requires_run=True) + def rotate_num_particles_in_moved_clusters(self): + """int: Total number of particles in clusters with accepted rotate moves. + + Note: + The counts are reset to 0 at the start of each `hoomd.Simulation.run`. + + """ + return self._cpp_obj.getCounters(1).rotate_num_particles_in_moved_clusters + + @log(category="scalar", requires_run=True) + def num_aabb_tree_builds(self): + """int: The number of times the AABB tree has been rebuilt. + + Note: + The counts are reset to 0 at the start of each `hoomd.Simulation.run`. + """ + return self._cpp_obj.getCounters(1).num_aabbtree_builds + + @log(category="scalar", requires_run=False) + def maximum_trial_rotation(self): + """float: Maximum size of rotation trial moves.""" + return self._cpp_obj.getMaximumTrialRotation + + class GCA(Updater): """Apply geometric cluster algorithm (GCA) moves. @@ -1023,4 +1240,5 @@ def complete(self): "MuVT", "QuickCompress", "Shape", + "VirtualClusterMoves", ] diff --git a/sphinx-doc/hoomd/hpmc/module-update.rst b/sphinx-doc/hoomd/hpmc/module-update.rst index 37f11ac58c..0e8d9162ab 100644 --- a/sphinx-doc/hoomd/hpmc/module-update.rst +++ b/sphinx-doc/hoomd/hpmc/module-update.rst @@ -3,7 +3,7 @@ update .. automodule:: hoomd.hpmc.update :members: - :exclude-members: BoxMC,GCA,MuVT,QuickCompress,Shape + :exclude-members: BoxMC,GCA,MuVT,QuickCompress,Shape,VirtualClusterMoves .. rubric:: Classes @@ -15,3 +15,4 @@ update update/muvt update/quickcompress update/shape + update/virtualclustermoves diff --git a/sphinx-doc/hoomd/hpmc/update/virtualclustermoves.rst b/sphinx-doc/hoomd/hpmc/update/virtualclustermoves.rst new file mode 100644 index 0000000000..53f9498ce6 --- /dev/null +++ b/sphinx-doc/hoomd/hpmc/update/virtualclustermoves.rst @@ -0,0 +1,9 @@ +VirtualClusterMoves +=================== + +.. py:currentmodule:: hoomd.hpmc.update + +.. autoclass:: VirtualClusterMoves + :members: + :show-inheritance: +