From 5f840a6a50662fa7d0767fb8b7a64939d2a34730 Mon Sep 17 00:00:00 2001 From: Brandon Butler Date: Wed, 9 Nov 2022 15:13:56 -0500 Subject: [PATCH 1/8] refactor!: Remove charge/diameter from rigid body definition. Allows setting charges in createRigidBodies, but not diameters since we are removing that. --- hoomd/md/ForceComposite.cc | 81 +++++++++++++++++++---------------- hoomd/md/ForceComposite.h | 35 ++++----------- hoomd/md/constrain.py | 9 ++-- hoomd/md/pytest/test_rigid.py | 23 +++++----- 4 files changed, 69 insertions(+), 79 deletions(-) diff --git a/hoomd/md/ForceComposite.cc b/hoomd/md/ForceComposite.cc index 98539849ee..ff7a25ea3b 100644 --- a/hoomd/md/ForceComposite.cc +++ b/hoomd/md/ForceComposite.cc @@ -8,6 +8,10 @@ #include #include +#ifndef __HIPCC__ +#include +#endif + /*! \file ForceComposite.cc \brief Contains code for the ForceComposite class */ @@ -54,9 +58,6 @@ ForceComposite::ForceComposite(std::shared_ptr sysdef) h_body_len.data[i] = 0; } - m_body_charge.resize(m_pdata->getNTypes()); - m_body_diameter.resize(m_pdata->getNTypes()); - m_d_max.resize(m_pdata->getNTypes(), Scalar(0.0)); m_d_max_changed.resize(m_pdata->getNTypes(), false); @@ -96,15 +97,11 @@ ForceComposite::~ForceComposite() void ForceComposite::setParam(unsigned int body_typeid, std::vector& type, std::vector& pos, - std::vector& orientation, - std::vector& charge, - std::vector& diameter) + std::vector& orientation) { assert(m_body_types.getPitch() >= m_pdata->getNTypes()); assert(m_body_pos.getPitch() >= m_pdata->getNTypes()); assert(m_body_orientation.getPitch() >= m_pdata->getNTypes()); - assert(m_body_charge.size() >= m_pdata->getNTypes()); - assert(m_body_diameter.size() >= m_pdata->getNTypes()); if (body_typeid >= m_pdata->getNTypes()) { @@ -118,21 +115,6 @@ void ForceComposite::setParam(unsigned int body_typeid, << " (position, orientation, type) are of unequal length."; throw std::runtime_error(error_msg.str()); } - if (charge.size() && charge.size() != pos.size()) - { - std::ostringstream error_msg; - error_msg << "Error initializing ForceComposite: Charges are non-empty but of different " - << "length than the positions."; - throw std::runtime_error(error_msg.str()); - } - if (diameter.size() && diameter.size() != pos.size()) - { - std::ostringstream error_msg; - error_msg << "Error initializing ForceComposite: Diameters are non-empty but of different " - << "length than the positions."; - throw std::runtime_error(error_msg.str()); - } - bool body_updated = false; bool body_len_changed = false; @@ -200,18 +182,12 @@ void ForceComposite::setParam(unsigned int body_typeid, access_location::host, access_mode::readwrite); - m_body_charge[body_typeid].resize(type.size()); - m_body_diameter[body_typeid].resize(type.size()); - // store body data in GlobalArray for (unsigned int i = 0; i < type.size(); ++i) { h_body_type.data[m_body_idx(body_typeid, i)] = type[i]; h_body_pos.data[m_body_idx(body_typeid, i)] = pos[i]; h_body_orientation.data[m_body_idx(body_typeid, i)] = orientation[i]; - - m_body_charge[body_typeid][i] = charge[i]; - m_body_diameter[body_typeid][i] = diameter[i]; } } m_bodies_changed = true; @@ -513,7 +489,39 @@ void ForceComposite::validateRigidBodies() m_particles_added_removed = false; } -void ForceComposite::createRigidBodies() +void ForceComposite::pyCreateRigidBodies(pybind11::dict charges) + { + if (pybind11::len(charges) == 0) + { + createRigidBodies(std::unordered_map>()); + return; + } + ArrayHandle h_body_len(m_body_len, access_location::host, access_mode::read); + std::unordered_map> charges_map; + for (const auto& item : charges) + { + const auto type = m_pdata->getTypeByName(item.first.cast()); + if (h_body_len.data[type] == 0) + { + throw std::runtime_error("Charge provided for non-central particle type."); + } + const auto charges_list = item.second.cast(); + if (pybind11::len(charges_list) != h_body_len.data[type]) + { + throw std::runtime_error("Charges provided not consistent with rigid body size."); + } + std::vector charges_vector; + for (auto& charge : charges_list) + { + charges_vector.emplace_back(charge.cast()); + } + charges_map.insert({type, charges_vector}); + } + createRigidBodies(charges_map); + } + +void ForceComposite::createRigidBodies( + const std::unordered_map> charges) { SnapshotParticleData snap; @@ -604,12 +612,11 @@ void ForceComposite::createRigidBodies() snap.type[constituent_particle_tag] = h_body_type.data[m_body_idx(body_type, current_body_index)]; snap.body[constituent_particle_tag] = particle_tag; - snap.charge[constituent_particle_tag] - = m_body_charge[body_type][current_body_index]; - snap.diameter[constituent_particle_tag] - = m_body_diameter[body_type][current_body_index]; - snap.pos[constituent_particle_tag] = snap.pos[particle_tag]; - + if (!charges.empty()) + { + snap.charge[constituent_particle_tag] + = charges.at(body_type)[current_body_index]; + } // Since the central particle tags here will be [0, n_central_particles), we know // that the molecule number will be the same as the central particle tag. molecule_tag[constituent_particle_tag] = particle_tag; @@ -1029,7 +1036,7 @@ void export_ForceComposite(pybind11::module& m) .def("setBody", &ForceComposite::setBody) .def("getBody", &ForceComposite::getBody) .def("validateRigidBodies", &ForceComposite::validateRigidBodies) - .def("createRigidBodies", &ForceComposite::createRigidBodies) + .def("createRigidBodies", &ForceComposite::pyCreateRigidBodies) .def("updateCompositeParticles", &ForceComposite::updateCompositeParticles); } diff --git a/hoomd/md/ForceComposite.h b/hoomd/md/ForceComposite.h index 98c792c74c..dd962ff88a 100644 --- a/hoomd/md/ForceComposite.h +++ b/hoomd/md/ForceComposite.h @@ -62,9 +62,7 @@ class PYBIND11_EXPORT ForceComposite : public MolecularForceCompute virtual void setParam(unsigned int body_typeid, std::vector& type, std::vector& pos, - std::vector& orientation, - std::vector& charge, - std::vector& diameter); + std::vector& orientation); //! Returns true because we compute the torque on the central particle virtual bool isAnisotropic() @@ -86,7 +84,11 @@ class PYBIND11_EXPORT ForceComposite : public MolecularForceCompute virtual void validateRigidBodies(); //! Create rigid body constituent particles - virtual void createRigidBodies(); + void pyCreateRigidBodies(pybind11::dict charges); + + //! Create rigid body constituent particles + virtual void + createRigidBodies(const std::unordered_map> charges); /// Construct from a Python dictionary void setBody(std::string typ, pybind11::object v) @@ -98,11 +100,9 @@ class PYBIND11_EXPORT ForceComposite : public MolecularForceCompute pybind11::list types = v["constituent_types"]; pybind11::list positions = v["positions"]; pybind11::list orientations = v["orientations"]; - pybind11::list charges = v["charges"]; - pybind11::list diameters = v["diameters"]; auto N = pybind11::len(positions); // Ensure proper list lengths - for (const auto& list : {types, orientations, charges, diameters}) + for (const auto& list : {types, orientations}) { if (pybind11::len(list) != N) { @@ -113,8 +113,6 @@ class PYBIND11_EXPORT ForceComposite : public MolecularForceCompute // extract the data from the python lists std::vector pos_vector; std::vector orientation_vector; - std::vector charge_vector; - std::vector diameter_vector; std::vector type_vector; for (size_t i(0); i < N; ++i) @@ -130,17 +128,10 @@ class PYBIND11_EXPORT ForceComposite : public MolecularForceCompute orientation_i[2].cast(), orientation_i[3].cast())); - charge_vector.emplace_back(charges[i].cast()); - diameter_vector.emplace_back(diameters[i].cast()); type_vector.emplace_back(m_pdata->getTypeByName(types[i].cast())); } - setParam(m_pdata->getTypeByName(typ), - type_vector, - pos_vector, - orientation_vector, - charge_vector, - diameter_vector); + setParam(m_pdata->getTypeByName(typ), type_vector, pos_vector, orientation_vector); } /// Convert parameters to a python dictionary @@ -166,8 +157,6 @@ class PYBIND11_EXPORT ForceComposite : public MolecularForceCompute pybind11::list positions; pybind11::list orientations; pybind11::list types; - pybind11::list charges; - pybind11::list diameters; for (unsigned int i = 0; i < N; i++) { @@ -181,15 +170,11 @@ class PYBIND11_EXPORT ForceComposite : public MolecularForceCompute static_cast(h_body_orientation.data[index].z), static_cast(h_body_orientation.data[index].w))); types.append(m_pdata->getNameByType(h_body_types.data[index])); - charges.append(m_body_charge[body_type_id][i]); - diameters.append(m_body_diameter[body_type_id][i]); } pybind11::dict v; v["constituent_types"] = types; v["positions"] = positions; v["orientations"] = orientations; - v["charges"] = charges; - v["diameters"] = diameters; return std::move(v); } @@ -202,9 +187,7 @@ class PYBIND11_EXPORT ForceComposite : public MolecularForceCompute GlobalArray m_body_orientation; //!< Constituent particle orientations per type id (2D) GlobalArray m_body_len; //!< Length of body per type id - std::vector> m_body_charge; //!< Constituent particle charges - std::vector> m_body_diameter; //!< Constituent particle diameters - Index2D m_body_idx; //!< Indexer for body parameters + Index2D m_body_idx; //!< Indexer for body parameters std::vector m_d_max; //!< Maximum body diameter per constituent particle type std::vector m_d_max_changed; //!< True if maximum body diameter changed (per type) diff --git a/hoomd/md/constrain.py b/hoomd/md/constrain.py index 6c9a90a0a4..7c120c7a75 100644 --- a/hoomd/md/constrain.py +++ b/hoomd/md/constrain.py @@ -290,19 +290,20 @@ def __init__(self): 'constituent_types': [str], 'positions': [(float,) * 3], 'orientations': [(float,) * 4], - 'charges': [float], - 'diameters': [float] }), allow_none=True), len_keys=1)) self._add_typeparam(body) self.body.default = None - def create_bodies(self, state): + def create_bodies(self, state, charges=None): r"""Create rigid bodies from central particles in state. Args: state (hoomd.State): The state in which to create rigid bodies. + charges (dict[str, float], optional): The charges for each of the + constituent particles, defaults to ``None``. If ``None``, all + charges are zero. The keys should be the central particles. `create_bodies` removes any existing constituent particles and adds new ones based on the body definitions in `body`. It overwrites all existing @@ -312,7 +313,7 @@ def create_bodies(self, state): raise RuntimeError( "Cannot call create_bodies after running simulation.") super()._attach(state._simulation) - self._cpp_obj.createRigidBodies() + self._cpp_obj.createRigidBodies({} if charges is None else charges) # Restore previous state self._detach() diff --git a/hoomd/md/pytest/test_rigid.py b/hoomd/md/pytest/test_rigid.py index d86b5c90f1..fa6724ef0f 100644 --- a/hoomd/md/pytest/test_rigid.py +++ b/hoomd/md/pytest/test_rigid.py @@ -30,8 +30,6 @@ def valid_body_definition(): [0, 1, 1 / (2**(1. / 2.))], ], "orientations": [(1.0, 0.0, 0.0, 0.0)] * 4, - "charges": [0.0, 1.0, 2.0, 3.5], - "diameters": [1.0, 1.5, 0.5, 1.0] } @@ -41,8 +39,6 @@ def test_body_setting(valid_body_definition): "positions": [[(1, 2)], [(1.0, 4.0, "foo")], 1.0, "hello"], "orientations": [[(1, 2, 3)], [(1.0, 4.0, 5.0, "foo")], [1.0], 1.0, "foo"], - "charges": [0.0, ["foo"]], - "diameters": [1.0, "foo", ["foo"]] } rigid = md.constrain.Rigid() @@ -67,7 +63,7 @@ def test_body_setting(valid_body_definition): current_body_definition[key] = valid_body_definition[key] -def check_bodies(snapshot, definition): +def check_bodies(snapshot, definition, charges=None): """Non-general assumes a snapshot from two_particle_snapshot_factory. This is just to prevent duplication of code from test_create_bodies and @@ -82,9 +78,10 @@ def check_bodies(snapshot, definition): assert all(snapshot.particles.body[6:] == 1) # check charges - for i in range(4): - assert snapshot.particles.charge[i + 2] == definition["charges"][i] - assert snapshot.particles.charge[i + 6] == definition["charges"][i] + if charges is not None: + for i in range(4): + assert snapshot.particles.charge[i + 2] == charges[i] + assert snapshot.particles.charge[i + 6] == charges[i] # check diameters for i in range(4): @@ -135,10 +132,11 @@ def test_create_bodies(simulation_factory, two_particle_snapshot_factory, initial_snapshot.particles.types = ["A", "B"] sim = simulation_factory(initial_snapshot) - rigid.create_bodies(sim.state) + charges = [1.0, 2.0, 3.0, 4.0] + rigid.create_bodies(sim.state, charges={"A": charges}) snapshot = sim.state.get_snapshot() if snapshot.communicator.rank == 0: - check_bodies(snapshot, valid_body_definition) + check_bodies(snapshot, valid_body_definition, charges) sim.operations.integrator = hoomd.md.Integrator(dt=0.005, rigid=rigid) # Ensure validate bodies passes @@ -246,12 +244,13 @@ def test_running_simulation(simulation_factory, two_particle_snapshot_factory, sim = simulation_factory(initial_snapshot) sim.seed = 5 - rigid.create_bodies(sim.state) + charges = [1.0, 2.0, 3.0, 4.0] + rigid.create_bodies(sim.state, charges={"A": charges}) sim.operations += integrator sim.run(5) snapshot = sim.state.get_snapshot() if sim.device.communicator.rank == 0: - check_bodies(snapshot, valid_body_definition) + check_bodies(snapshot, valid_body_definition, charges) autotuned_kernel_parameter_check(instance=rigid, activate=lambda: sim.run(1)) From e35f8ce46c9fbe118b502e960897b938dbecb89e Mon Sep 17 00:00:00 2001 From: Brandon Butler Date: Thu, 10 Nov 2022 11:38:18 -0500 Subject: [PATCH 2/8] fix: Remove diameter from rigid body tests --- hoomd/md/pytest/test_rigid.py | 5 ----- 1 file changed, 5 deletions(-) diff --git a/hoomd/md/pytest/test_rigid.py b/hoomd/md/pytest/test_rigid.py index fa6724ef0f..3908c6a31c 100644 --- a/hoomd/md/pytest/test_rigid.py +++ b/hoomd/md/pytest/test_rigid.py @@ -83,11 +83,6 @@ def check_bodies(snapshot, definition, charges=None): assert snapshot.particles.charge[i + 2] == charges[i] assert snapshot.particles.charge[i + 6] == charges[i] - # check diameters - for i in range(4): - assert snapshot.particles.diameter[i + 2] == definition["diameters"][i] - assert snapshot.particles.diameter[i + 6] == definition["diameters"][i] - particle_one = (snapshot.particles.position[0], snapshot.particles.orientation[0]) particle_two = (snapshot.particles.position[1], From 6e1d4863e14ffe20af0476a65197c14c8d43b0c3 Mon Sep 17 00:00:00 2001 From: Brandon Butler Date: Tue, 22 Nov 2022 11:31:49 -0500 Subject: [PATCH 3/8] refactor: Remove unnecessary ifdef guard Co-authored-by: Joshua A. Anderson --- hoomd/md/ForceComposite.cc | 2 -- 1 file changed, 2 deletions(-) diff --git a/hoomd/md/ForceComposite.cc b/hoomd/md/ForceComposite.cc index ff7a25ea3b..7c0370457c 100644 --- a/hoomd/md/ForceComposite.cc +++ b/hoomd/md/ForceComposite.cc @@ -8,9 +8,7 @@ #include #include -#ifndef __HIPCC__ #include -#endif /*! \file ForceComposite.cc \brief Contains code for the ForceComposite class From 70eb630592c10e75ab7bd1572b6270be35858ad5 Mon Sep 17 00:00:00 2001 From: Brandon Butler Date: Tue, 22 Nov 2022 11:32:20 -0500 Subject: [PATCH 4/8] doc: Fix type documentation for Rigid.create_bodies Co-authored-by: Joshua A. Anderson --- hoomd/md/constrain.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/hoomd/md/constrain.py b/hoomd/md/constrain.py index 7c120c7a75..ce0ad971f5 100644 --- a/hoomd/md/constrain.py +++ b/hoomd/md/constrain.py @@ -301,7 +301,7 @@ def create_bodies(self, state, charges=None): Args: state (hoomd.State): The state in which to create rigid bodies. - charges (dict[str, float], optional): The charges for each of the + charges (dict[str, list[float]], optional): The charges for each of the constituent particles, defaults to ``None``. If ``None``, all charges are zero. The keys should be the central particles. From eb815eb47b40871657c4d729899d33b1a4990ad4 Mon Sep 17 00:00:00 2001 From: Brandon Butler Date: Tue, 22 Nov 2022 11:48:06 -0500 Subject: [PATCH 5/8] fix: Use Scalar typedef for charges in createRigidBodies --- hoomd/md/ForceComposite.cc | 10 +++++----- hoomd/md/ForceComposite.h | 2 +- 2 files changed, 6 insertions(+), 6 deletions(-) diff --git a/hoomd/md/ForceComposite.cc b/hoomd/md/ForceComposite.cc index 7c0370457c..19664a3d68 100644 --- a/hoomd/md/ForceComposite.cc +++ b/hoomd/md/ForceComposite.cc @@ -491,11 +491,11 @@ void ForceComposite::pyCreateRigidBodies(pybind11::dict charges) { if (pybind11::len(charges) == 0) { - createRigidBodies(std::unordered_map>()); + createRigidBodies(std::unordered_map>()); return; } ArrayHandle h_body_len(m_body_len, access_location::host, access_mode::read); - std::unordered_map> charges_map; + std::unordered_map> charges_map; for (const auto& item : charges) { const auto type = m_pdata->getTypeByName(item.first.cast()); @@ -508,10 +508,10 @@ void ForceComposite::pyCreateRigidBodies(pybind11::dict charges) { throw std::runtime_error("Charges provided not consistent with rigid body size."); } - std::vector charges_vector; + std::vector charges_vector; for (auto& charge : charges_list) { - charges_vector.emplace_back(charge.cast()); + charges_vector.emplace_back(charge.cast()); } charges_map.insert({type, charges_vector}); } @@ -519,7 +519,7 @@ void ForceComposite::pyCreateRigidBodies(pybind11::dict charges) } void ForceComposite::createRigidBodies( - const std::unordered_map> charges) + const std::unordered_map> charges) { SnapshotParticleData snap; diff --git a/hoomd/md/ForceComposite.h b/hoomd/md/ForceComposite.h index dd962ff88a..58353546bb 100644 --- a/hoomd/md/ForceComposite.h +++ b/hoomd/md/ForceComposite.h @@ -88,7 +88,7 @@ class PYBIND11_EXPORT ForceComposite : public MolecularForceCompute //! Create rigid body constituent particles virtual void - createRigidBodies(const std::unordered_map> charges); + createRigidBodies(const std::unordered_map> charges); /// Construct from a Python dictionary void setBody(std::string typ, pybind11::object v) From c60d10aaac56d2c9939a117019f07337666645de Mon Sep 17 00:00:00 2001 From: Brandon Butler Date: Tue, 22 Nov 2022 13:45:20 -0500 Subject: [PATCH 6/8] test: Fix hoomd.filter.Rigid tests. Needed to fix rigid body definition. --- hoomd/md/pytest/test_filter_md.py | 2 -- 1 file changed, 2 deletions(-) diff --git a/hoomd/md/pytest/test_filter_md.py b/hoomd/md/pytest/test_filter_md.py index e55f3173c7..002a711414 100644 --- a/hoomd/md/pytest/test_filter_md.py +++ b/hoomd/md/pytest/test_filter_md.py @@ -36,8 +36,6 @@ def test_rigid_filter(make_filter_snapshot, simulation_factory): [0, 1, 1 / (2**(1. / 2.))], ], "orientations": [(1.0, 0.0, 0.0, 0.0)] * 4, - "charges": [0.0, 1.0, 2.0, 3.5], - "diameters": [1.0, 1.5, 0.5, 1.0] } snapshot = make_filter_snapshot(n=100, particle_types=["A", "B", "C"]) From 5d07fe15ab3a6c635afece61e89c51299425ec78 Mon Sep 17 00:00:00 2001 From: "Joshua A. Anderson" Date: Tue, 22 Nov 2022 14:05:53 -0500 Subject: [PATCH 7/8] Update hoomd/md/constrain.py --- hoomd/md/constrain.py | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/hoomd/md/constrain.py b/hoomd/md/constrain.py index ce0ad971f5..41cfa38c0b 100644 --- a/hoomd/md/constrain.py +++ b/hoomd/md/constrain.py @@ -301,9 +301,9 @@ def create_bodies(self, state, charges=None): Args: state (hoomd.State): The state in which to create rigid bodies. - charges (dict[str, list[float]], optional): The charges for each of the - constituent particles, defaults to ``None``. If ``None``, all - charges are zero. The keys should be the central particles. + charges (dict[str, list[float]], optional): The charges for each of + the constituent particles, defaults to ``None``. If ``None``, + all charges are zero. The keys should be the central particles. `create_bodies` removes any existing constituent particles and adds new ones based on the body definitions in `body`. It overwrites all existing From 8ae46a811894756997f61b37e50d35ed27b09d75 Mon Sep 17 00:00:00 2001 From: Brandon Butler Date: Tue, 22 Nov 2022 15:46:03 -0500 Subject: [PATCH 8/8] fix: Correctly set constituent particles pos in create bodies Need to place constituent particles where they will be decomposed to the same rank. --- hoomd/md/ForceComposite.cc | 2 ++ 1 file changed, 2 insertions(+) diff --git a/hoomd/md/ForceComposite.cc b/hoomd/md/ForceComposite.cc index 19664a3d68..92d65bcd2d 100644 --- a/hoomd/md/ForceComposite.cc +++ b/hoomd/md/ForceComposite.cc @@ -528,6 +528,7 @@ void ForceComposite::createRigidBodies( bool remove_existing_bodies = false; unsigned int n_constituent_particles = 0; unsigned int n_free_bodies = 0; + if (m_exec_conf->getRank() == 0) { // We restrict the scope of h_body_len to ensure if we remove_existing_bodies or in any way // reallocated m_body_len to a new memory location that we will be forced to reaquire the @@ -610,6 +611,7 @@ void ForceComposite::createRigidBodies( snap.type[constituent_particle_tag] = h_body_type.data[m_body_idx(body_type, current_body_index)]; snap.body[constituent_particle_tag] = particle_tag; + snap.pos[constituent_particle_tag] = snap.pos[particle_tag]; if (!charges.empty()) { snap.charge[constituent_particle_tag]