From b01edb1224ad0838c0f14a602b56db0582253344 Mon Sep 17 00:00:00 2001 From: Tobias Dwyer Date: Sun, 28 Jan 2024 22:59:12 -0500 Subject: [PATCH 1/5] added muVT wall bug patch --- hoomd/hpmc/ExternalField.h | 20 ++++++++ hoomd/hpmc/ExternalFieldJIT.h | 3 +- hoomd/hpmc/ExternalFieldWall.h | 50 +++++++++++++++++++ hoomd/hpmc/PatchEnergyJIT.cc | 1 - hoomd/hpmc/UpdaterMuVT.h | 15 ++++++ hoomd/hpmc/pytest/test_muvt.py | 87 +++++++++++++++++++++++++++++++++- 6 files changed, 173 insertions(+), 3 deletions(-) diff --git a/hoomd/hpmc/ExternalField.h b/hoomd/hpmc/ExternalField.h index e647b04a18..d3246db334 100644 --- a/hoomd/hpmc/ExternalField.h +++ b/hoomd/hpmc/ExternalField.h @@ -39,6 +39,25 @@ class ExternalField : public Compute return 0; } + //! Evaluate the energy of the force. + /*! \param box The system box. + \param type Particle type. + \param r_i Particle position + \param q_i Particle orientation. + \param diameter Particle diameter. + \param charge Particle charge. + \returns Energy due to the force + */ + virtual float energy(const BoxDim& box, + unsigned int type, + const vec3& r_i, + const quat& q_i, + Scalar diameter, + Scalar charge) + { + return 0; + } + virtual bool hasVolume() { return false; @@ -88,6 +107,7 @@ template void export_ExternalFieldInterface(pybind11::module& m, st .def(pybind11::init>()) .def("compute", &ExternalFieldMono::compute) .def("energydiff", &ExternalFieldMono::energydiff) + .def("energy", &ExternalFieldMono::energy) .def("calculateDeltaE", &ExternalFieldMono::calculateDeltaE); } diff --git a/hoomd/hpmc/ExternalFieldJIT.h b/hoomd/hpmc/ExternalFieldJIT.h index df9221522b..5f627c2a34 100644 --- a/hoomd/hpmc/ExternalFieldJIT.h +++ b/hoomd/hpmc/ExternalFieldJIT.h @@ -283,7 +283,8 @@ template void export_ExternalFieldJIT(pybind11::module& m, std::str const std::vector&, pybind11::array_t>()) .def("computeEnergy", &ExternalFieldJIT::computeEnergy) - .def_property_readonly("param_array", &ExternalFieldJIT::getParamArray); + .def_property_readonly("param_array", &ExternalFieldJIT::getParamArray) + .def("energy", &ExternalFieldJIT::energy); } } // end namespace hpmc diff --git a/hoomd/hpmc/ExternalFieldWall.h b/hoomd/hpmc/ExternalFieldWall.h index b36e830408..025c54c6c8 100644 --- a/hoomd/hpmc/ExternalFieldWall.h +++ b/hoomd/hpmc/ExternalFieldWall.h @@ -629,6 +629,56 @@ template class ExternalFieldWall : public ExternalFieldMono return double(0.0); } + //! Evaluate the energy of the force. + /*! \param box The system box. + \param type Particle type. + \param r_i Particle position + \param q_i Particle orientation. + \param diameter Particle diameter. + \param charge Particle charge. + \returns Energy due to the force + */ + virtual float energy(const BoxDim& box, + unsigned int type, + const vec3& r_i, + const quat& q_i, + Scalar diameter, + Scalar charge) + { + const std::vector>& params + = m_mc->getParams(); + Shape shape(q_i, params[type]); + vec3 origin(m_pdata->getOrigin()); + + for (size_t i = 0; i < m_Spheres.size(); i++) + { + if (!test_confined(m_Spheres[i], shape, r_i, origin, box)) + { + return INFINITY; + } + } + + for (size_t i = 0; i < m_Cylinders.size(); i++) + { + set_cylinder_wall_verts(m_Cylinders[i], shape); + if (!test_confined(m_Cylinders[i], shape, r_i, origin, box)) + { + return INFINITY; + } + } + + for (size_t i = 0; i < m_Planes.size(); i++) + { + if (!test_confined(m_Planes[i], shape, r_i, origin, box)) + { + return INFINITY; + } + } + + return double(0.0); + } + double calculateDeltaE(uint64_t timestep, const Scalar4* const position_old, const Scalar4* const orientation_old, diff --git a/hoomd/hpmc/PatchEnergyJIT.cc b/hoomd/hpmc/PatchEnergyJIT.cc index 1713a42fb7..c076188c3a 100644 --- a/hoomd/hpmc/PatchEnergyJIT.cc +++ b/hoomd/hpmc/PatchEnergyJIT.cc @@ -69,7 +69,6 @@ void export_PatchEnergyJIT(pybind11::module& m) Scalar, pybind11::array_t>()) .def_property("r_cut", &PatchEnergyJIT::getRCut, &PatchEnergyJIT::setRCut) - .def("energy", &PatchEnergyJIT::energy) .def_property_readonly("param_array", &PatchEnergyJIT::getParamArray); } diff --git a/hoomd/hpmc/UpdaterMuVT.h b/hoomd/hpmc/UpdaterMuVT.h index d61f059668..9cdff8c8f4 100644 --- a/hoomd/hpmc/UpdaterMuVT.h +++ b/hoomd/hpmc/UpdaterMuVT.h @@ -1911,6 +1911,9 @@ bool UpdaterMuVT::tryInsertParticle(uint64_t timestep, // do we have to compute energetic contribution? auto patch = m_mc->getPatchEnergy(); + // do we have to compute a wall contribution + auto field = m_mc->getExternalField(); + lnboltzmann = Scalar(0.0); unsigned int overlap = 0; @@ -1945,6 +1948,18 @@ bool UpdaterMuVT::tryInsertParticle(uint64_t timestep, ShortReal r_cut_patch(0.0); Scalar r_cut_self(0.0); + if (field) + { + const BoxDim box = this->m_pdata->getGlobalBox(); + lnboltzmann -= field->energy(box, + type, + pos, + quat(orientation), + 1.0, // diameter i + 0.0 // charge i + ); + } + if (patch) { r_cut_patch = ShortReal(patch->getRCut() + 0.5 * patch->getAdditiveCutoff(type)); diff --git a/hoomd/hpmc/pytest/test_muvt.py b/hoomd/hpmc/pytest/test_muvt.py index 21a4c08e0d..ec82183bd8 100644 --- a/hoomd/hpmc/pytest/test_muvt.py +++ b/hoomd/hpmc/pytest/test_muvt.py @@ -5,6 +5,7 @@ import hoomd import pytest +import numpy import hoomd.hpmc.pytest.conftest # note: The parameterized tests validate parameters so we can't pass in values @@ -121,7 +122,6 @@ def test_valid_setattr_attached(device, attr, value, simulation_factory, setattr(muvt, attr, value) assert getattr(muvt, attr) == value - def test_insertion_removal(device, simulation_factory, lattice_snapshot_factory): """Test that MuVT is able to insert and remove particles.""" @@ -159,3 +159,88 @@ def test_insertion_removal(device, simulation_factory, # make a wild guess: there be B particles assert (muvt.N['B'] > 0) + +def test_plane_wall_insertion(device, simulation_factory, + one_particle_snapshot_factory): + """Test that MuVT considers a planar wall when inserting particles.""" + sim = simulation_factory( + one_particle_snapshot_factory(particle_types=['A'], + dimensions=3, + position=(0, 0, 5), + orientation=(1, 0, 0, 0), + L=20) + ) + + sphere_radius = 0.6 + mc = hoomd.hpmc.integrate.Sphere(default_d=0.0, default_a=0.0) + mc.shape['A'] = dict(diameter=2*sphere_radius) + walls = [hoomd.wall.Plane(origin=(0,0,0),normal=(0,0,1))] + wall_potential = hoomd.hpmc.external.wall.WallPotential(walls) + mc.external_potential = wall_potential + sim.operations.integrator = mc + + muvt = hoomd.hpmc.update.MuVT(trigger=hoomd.trigger.Periodic(1), + transfer_types=['A']) + muvt.fugacity['A'] = 1000 + sim.operations.updaters.append(muvt) + sim.run(300) + snapshot = sim.state.get_snapshot() + pos = snapshot.particles.position + assert numpy.min(pos[:,2]) >= sphere_radius + +def test_spherical_wall_insertion(device, simulation_factory, + one_particle_snapshot_factory): + """Test that MuVT considers a spherical wall when inserting particles.""" + sim = simulation_factory( + one_particle_snapshot_factory(particle_types=['A'], + dimensions=3, + position=(0, 0, 0), + orientation=(1, 0, 0, 0), + L=20) + ) + + + mc = hoomd.hpmc.integrate.Sphere(default_d=0.1, default_a=0.1) + sphere_radius = 0.6 + mc.shape['A'] = dict(diameter=2*sphere_radius) + walls = [hoomd.wall.Sphere(radius=5)] + wall_potential = hoomd.hpmc.external.wall.WallPotential(walls) + mc.external_potential = wall_potential + sim.operations.integrator = mc + + muvt = hoomd.hpmc.update.MuVT(trigger=hoomd.trigger.Periodic(1), + transfer_types=['A']) + muvt.fugacity['A'] = 1000 + sim.operations.updaters.append(muvt) + sim.run(300) + snapshot = sim.state.get_snapshot() + pos = snapshot.particles.position + assert numpy.max(numpy.linalg.norm(pos,axis=1))-sphere_radius <= 5 + +def test_cylindrical_wall_insertion(device, simulation_factory, + one_particle_snapshot_factory): + """Test that MuVT considers a cylindrical wall when inserting particles.""" + sim = simulation_factory( + one_particle_snapshot_factory(particle_types=['A'], + dimensions=3, + position=(0, 0, 0), + orientation=(1, 0, 0, 0), + L=20) + ) + + sphere_radius = 0.6 + mc = hoomd.hpmc.integrate.Sphere(default_d=0.0, default_a=0.0) + mc.shape['A'] = dict(diameter=2 * sphere_radius) + walls = [hoomd.wall.Cylinder(radius=5,axis=(0,0,1))] + wall_potential = hoomd.hpmc.external.wall.WallPotential(walls) + mc.external_potential = wall_potential + sim.operations.integrator = mc + + muvt = hoomd.hpmc.update.MuVT(trigger=hoomd.trigger.Periodic(1), + transfer_types=['A']) + muvt.fugacity['A'] = 1000 + sim.operations.updaters.append(muvt) + sim.run(300) + snapshot = sim.state.get_snapshot() + pos = snapshot.particles.position + assert numpy.max(numpy.linalg.norm(pos[:,:2],axis=1))-sphere_radius <= 5 From d675fa16e21b130c79555518b3357ec3fa9d254e Mon Sep 17 00:00:00 2001 From: Tobias-Dwyer Date: Tue, 30 Jan 2024 14:22:29 -0500 Subject: [PATCH 2/5] added jit wall removal and testing in MuVT --- hoomd/hpmc/ExternalFieldJIT.h | 3 +- hoomd/hpmc/UpdaterMuVT.h | 25 ++++++++++++++- hoomd/hpmc/pytest/test_muvt.py | 58 +++++++++++++++++++++++++++++++++- 3 files changed, 82 insertions(+), 4 deletions(-) diff --git a/hoomd/hpmc/ExternalFieldJIT.h b/hoomd/hpmc/ExternalFieldJIT.h index 5f627c2a34..df9221522b 100644 --- a/hoomd/hpmc/ExternalFieldJIT.h +++ b/hoomd/hpmc/ExternalFieldJIT.h @@ -283,8 +283,7 @@ template void export_ExternalFieldJIT(pybind11::module& m, std::str const std::vector&, pybind11::array_t>()) .def("computeEnergy", &ExternalFieldJIT::computeEnergy) - .def_property_readonly("param_array", &ExternalFieldJIT::getParamArray) - .def("energy", &ExternalFieldJIT::energy); + .def_property_readonly("param_array", &ExternalFieldJIT::getParamArray); } } // end namespace hpmc diff --git a/hoomd/hpmc/UpdaterMuVT.h b/hoomd/hpmc/UpdaterMuVT.h index 9cdff8c8f4..e01d30e733 100644 --- a/hoomd/hpmc/UpdaterMuVT.h +++ b/hoomd/hpmc/UpdaterMuVT.h @@ -1650,6 +1650,28 @@ bool UpdaterMuVT::tryRemoveParticle(uint64_t timestep, unsigned int tag, // do we have to compute energetic contribution? auto patch = m_mc->getPatchEnergy(); + // do we have to compute a wall contribution? + auto field = m_mc->getExternalField(); + + if (field) + { + // getPosition() takes into account grid shift, correct for that + Scalar3 p = m_pdata->getPosition(tag) + m_pdata->getOrigin(); + int3 tmp = make_int3(0, 0, 0); + m_pdata->getGlobalBox().wrap(p, tmp); + vec3 pos(p); + const BoxDim box = this->m_pdata->getGlobalBox(); + unsigned int type = this->m_pdata->getType(tag); + quat orientation(m_pdata->getOrientation(tag)); + lnboltzmann += field->energy(box, + type, + pos, + quat(orientation), + 1.0, // diameter i + 0.0 // charge i + ); + } + // if not, no overlaps generated if (patch) { @@ -1795,6 +1817,7 @@ bool UpdaterMuVT::tryRemoveParticle(uint64_t timestep, unsigned int tag, } #endif } + } // Depletants @@ -1911,7 +1934,7 @@ bool UpdaterMuVT::tryInsertParticle(uint64_t timestep, // do we have to compute energetic contribution? auto patch = m_mc->getPatchEnergy(); - // do we have to compute a wall contribution + // do we have to compute a wall contribution? auto field = m_mc->getExternalField(); lnboltzmann = Scalar(0.0); diff --git a/hoomd/hpmc/pytest/test_muvt.py b/hoomd/hpmc/pytest/test_muvt.py index ec82183bd8..c82f3bbc01 100644 --- a/hoomd/hpmc/pytest/test_muvt.py +++ b/hoomd/hpmc/pytest/test_muvt.py @@ -160,6 +160,57 @@ def test_insertion_removal(device, simulation_factory, # make a wild guess: there be B particles assert (muvt.N['B'] > 0) +@pytest.mark.skipif(not hoomd.version.llvm_enabled, reason='LLVM not enabled') +def test_jit_remove_insert(device, simulation_factory, + one_particle_snapshot_factory): + """Test that MuVT considers a cpp potential when removing particles by adding a high energy particles that should always be removed by the updater while inserting particles to the top of the box.""" + sim = simulation_factory( + one_particle_snapshot_factory(particle_types=['A'], + dimensions=3, + position=(0, 0, -5), + orientation=(1, 0, 0, 0), + L=20) + ) + + sphere_radius = 0.6 + mc = hoomd.hpmc.integrate.Sphere(default_d=0.0, default_a=0.0) + mc.shape['A'] = dict(diameter=2*sphere_radius) + + # code returns infinity if the particle center goes past the center of the box + + center_wall = ''' + if (r_i.z < 0. || r_i.z > box.getL().z/2) + return INFINITY; + else + return 0.0f; + ''' + cpp_external_potential = hoomd.hpmc.external.user.CPPExternalPotential( + code=center_wall) + + mc.external_potential = cpp_external_potential + + sim.operations.integrator = mc + + muvt = hoomd.hpmc.update.MuVT(trigger=hoomd.trigger.Periodic(1), + transfer_types=['A']) + muvt.fugacity['A'] = 1e6 + sim.operations.updaters.append(muvt) + sim.run(1000) + snapshot = sim.state.get_snapshot() + pos = snapshot.particles.position + + # We should have inserted particles successfully + assert muvt.insert_moves[0] > 0 + + # We should have successfully removed the high energy particle + assert muvt.remove_moves[0] > 0 + + # We should have successfully inserted only to the top of the box and removed the high energy particle + assert numpy.min(pos[:,2]) >= 0. + + # We should have added more than one particle to the box + assert len(pos) > 1 + def test_plane_wall_insertion(device, simulation_factory, one_particle_snapshot_factory): """Test that MuVT considers a planar wall when inserting particles.""" @@ -186,7 +237,9 @@ def test_plane_wall_insertion(device, simulation_factory, sim.run(300) snapshot = sim.state.get_snapshot() pos = snapshot.particles.position + # Test if inserted spheres are above the plane assert numpy.min(pos[:,2]) >= sphere_radius + assert len(pos) > 1 def test_spherical_wall_insertion(device, simulation_factory, one_particle_snapshot_factory): @@ -199,7 +252,6 @@ def test_spherical_wall_insertion(device, simulation_factory, L=20) ) - mc = hoomd.hpmc.integrate.Sphere(default_d=0.1, default_a=0.1) sphere_radius = 0.6 mc.shape['A'] = dict(diameter=2*sphere_radius) @@ -215,7 +267,9 @@ def test_spherical_wall_insertion(device, simulation_factory, sim.run(300) snapshot = sim.state.get_snapshot() pos = snapshot.particles.position + # Test if inserted spheres are inside the spherical wall assert numpy.max(numpy.linalg.norm(pos,axis=1))-sphere_radius <= 5 + assert len(pos) > 1 def test_cylindrical_wall_insertion(device, simulation_factory, one_particle_snapshot_factory): @@ -243,4 +297,6 @@ def test_cylindrical_wall_insertion(device, simulation_factory, sim.run(300) snapshot = sim.state.get_snapshot() pos = snapshot.particles.position + # Test if inserted spheres are inside the cylinder wall assert numpy.max(numpy.linalg.norm(pos[:,:2],axis=1))-sphere_radius <= 5 + assert len(pos) > 1 From 3031f6c2214df5ddc1792104b4acca9cc2b61a02 Mon Sep 17 00:00:00 2001 From: Tobias-Dwyer Date: Thu, 1 Feb 2024 16:45:22 -0500 Subject: [PATCH 3/5] add code to run on mpi --- hoomd/hpmc/ExternalFieldJIT.h | 48 +++---- hoomd/hpmc/UpdaterMuVT.h | 71 +++++----- hoomd/hpmc/pytest/test_muvt.py | 241 +++++++++++++++++++-------------- 3 files changed, 205 insertions(+), 155 deletions(-) diff --git a/hoomd/hpmc/ExternalFieldJIT.h b/hoomd/hpmc/ExternalFieldJIT.h index df9221522b..c2e7b1ae0f 100644 --- a/hoomd/hpmc/ExternalFieldJIT.h +++ b/hoomd/hpmc/ExternalFieldJIT.h @@ -66,6 +66,16 @@ template class ExternalFieldJIT : public hpmc::ExternalFieldMono(factory); } + float energy_no_wrap(const BoxDim& box, + unsigned int type, + const vec3& r_i, + const quat& q_i, + Scalar diameter, + Scalar charge) + { + return m_eval(box, type, r_i, q_i, diameter, charge); + } + //! Evaluate the energy of the force. /*! \param box The system box. \param type Particle type. @@ -82,7 +92,10 @@ template class ExternalFieldJIT : public hpmc::ExternalFieldMono r_i_wrapped = r_i - vec3(this->m_pdata->getOrigin()); + int3 image = make_int3(0, 0, 0); + box.wrap(r_i_wrapped, image); + return energy_no_wrap(box, type, r_i_wrapped, q_i, diameter, charge); } //! Computes the total field energy of the system at the current state @@ -109,13 +122,10 @@ template class ExternalFieldJIT : public hpmc::ExternalFieldMono pos_i = vec3(postype_i) - vec3(this->m_pdata->getOrigin()); - int3 image = make_int3(0, 0, 0); - box.wrap(pos_i, image); total_energy += energy(box, typ_i, - pos_i, + vec3(postype_i), quat(h_orientation.data[i]), h_diameter.data[i], h_charge.data[i]); @@ -171,23 +181,20 @@ template class ExternalFieldJIT : public hpmc::ExternalFieldMono pos_i = vec3(postype_i) - vec3(this->m_pdata->getOrigin()); - box_new.wrap(pos_i, image); - image = make_int3(0, 0, 0); vec3 old_pos_i = vec3(*(position_old + i)) - vec3(origin_old); box_old.wrap(old_pos_i, image); dE += energy(box_new, typ_i, - pos_i, + vec3(postype_i), quat(h_orientation.data[i]), h_diameter.data[i], h_charge.data[i]); - dE -= energy(box_old, - typ_i, - old_pos_i, - quat(*(orientation_old + i)), - h_diameter.data[i], - h_charge.data[i]); + dE -= energy_no_wrap(box_old, + typ_i, + old_pos_i, + quat(*(orientation_old + i)), + h_diameter.data[i], + h_charge.data[i]); } #ifdef ENABLE_MPI if (this->m_sysdef->isDomainDecomposed()) @@ -223,24 +230,17 @@ template class ExternalFieldJIT : public hpmc::ExternalFieldMono h_charge(this->m_pdata->getCharges(), access_location::host, access_mode::read); - const BoxDim box = this->m_pdata->getGlobalBox(); - int3 image = make_int3(0, 0, 0); - vec3 pos_new_shifted = position_new - vec3(this->m_pdata->getOrigin()); - vec3 pos_old_shifted = position_old - vec3(this->m_pdata->getOrigin()); - box.wrap(pos_new_shifted, image); - image = make_int3(0, 0, 0); - box.wrap(pos_old_shifted, image); double dE = 0.0; dE += energy(this->m_pdata->getGlobalBox(), typ_i, - pos_new_shifted, + position_new, shape_new.orientation, h_diameter.data[index], h_charge.data[index]); dE -= energy(this->m_pdata->getGlobalBox(), typ_i, - pos_old_shifted, + position_old, shape_old.orientation, h_diameter.data[index], h_charge.data[index]); diff --git a/hoomd/hpmc/UpdaterMuVT.h b/hoomd/hpmc/UpdaterMuVT.h index e01d30e733..7a8721d889 100644 --- a/hoomd/hpmc/UpdaterMuVT.h +++ b/hoomd/hpmc/UpdaterMuVT.h @@ -1650,27 +1650,31 @@ bool UpdaterMuVT::tryRemoveParticle(uint64_t timestep, unsigned int tag, // do we have to compute energetic contribution? auto patch = m_mc->getPatchEnergy(); - // do we have to compute a wall contribution? - auto field = m_mc->getExternalField(); + // do we have to compute a wall contribution? + auto field = m_mc->getExternalField(); + unsigned int p = m_exec_conf->getPartition() % m_npartition; - if (field) - { - // getPosition() takes into account grid shift, correct for that + if (field && (!m_gibbs || p == 0)) + { + // getPosition() takes into account grid shift, correct for that Scalar3 p = m_pdata->getPosition(tag) + m_pdata->getOrigin(); - int3 tmp = make_int3(0, 0, 0); - m_pdata->getGlobalBox().wrap(p, tmp); - vec3 pos(p); - const BoxDim box = this->m_pdata->getGlobalBox(); - unsigned int type = this->m_pdata->getType(tag); - quat orientation(m_pdata->getOrientation(tag)); - lnboltzmann += field->energy(box, - type, - pos, - quat(orientation), - 1.0, // diameter i - 0.0 // charge i - ); - } + int3 tmp = make_int3(0, 0, 0); + m_pdata->getGlobalBox().wrap(p, tmp); + vec3 pos(p); + const BoxDim box = this->m_pdata->getGlobalBox(); + unsigned int type = this->m_pdata->getType(tag); + quat orientation(m_pdata->getOrientation(tag)); + if (is_local) + { + lnboltzmann += field->energy(box, + type, + pos, + quat(orientation), + 1.0, // diameter i + 0.0 // charge i + ); + } + } // if not, no overlaps generated if (patch) @@ -1804,20 +1808,19 @@ bool UpdaterMuVT::tryRemoveParticle(uint64_t timestep, unsigned int tag, } // end loop over AABB nodes } // end loop over images } + } #ifdef ENABLE_MPI - if (m_sysdef->isDomainDecomposed()) - { - MPI_Allreduce(MPI_IN_PLACE, - &lnboltzmann, - 1, - MPI_HOOMD_SCALAR, - MPI_SUM, - m_exec_conf->getMPICommunicator()); - } -#endif + if (m_sysdef->isDomainDecomposed()) + { + MPI_Allreduce(MPI_IN_PLACE, + &lnboltzmann, + 1, + MPI_HOOMD_SCALAR, + MPI_SUM, + m_exec_conf->getMPICommunicator()); } - +#endif } // Depletants @@ -1971,15 +1974,17 @@ bool UpdaterMuVT::tryInsertParticle(uint64_t timestep, ShortReal r_cut_patch(0.0); Scalar r_cut_self(0.0); - if (field) + unsigned int p = m_exec_conf->getPartition() % m_npartition; + + if (field && (!m_gibbs || p == 0)) { - const BoxDim box = this->m_pdata->getGlobalBox(); + const BoxDim& box = this->m_pdata->getGlobalBox(); lnboltzmann -= field->energy(box, type, pos, quat(orientation), 1.0, // diameter i - 0.0 // charge i + 0.0 // charge i ); } diff --git a/hoomd/hpmc/pytest/test_muvt.py b/hoomd/hpmc/pytest/test_muvt.py index c82f3bbc01..a8b44a7c06 100644 --- a/hoomd/hpmc/pytest/test_muvt.py +++ b/hoomd/hpmc/pytest/test_muvt.py @@ -11,19 +11,25 @@ # note: The parameterized tests validate parameters so we can't pass in values # here that require preprocessing valid_constructor_args = [ - dict(trigger=hoomd.trigger.Periodic(10), - transfer_types=['A'], - max_volume_rescale=0.2, - volume_move_probability=0.5), - dict(trigger=hoomd.trigger.After(100), transfer_types=['A', 'B']), + dict( + trigger=hoomd.trigger.Periodic(10), + transfer_types=["A"], + max_volume_rescale=0.2, + volume_move_probability=0.5, + ), + dict(trigger=hoomd.trigger.After(100), transfer_types=["A", "B"]), ] -valid_attrs = [('trigger', hoomd.trigger.Periodic(10000)), - ('trigger', hoomd.trigger.After(100)), - ('trigger', hoomd.trigger.Before(12345)), - ('volume_move_probability', 0.2), ('max_volume_rescale', 0.42), - ('transfer_types', ['A']), ('transfer_types', ['B']), - ('transfer_types', ['A', 'B'])] +valid_attrs = [ + ("trigger", hoomd.trigger.Periodic(10000)), + ("trigger", hoomd.trigger.After(100)), + ("trigger", hoomd.trigger.Before(12345)), + ("volume_move_probability", 0.2), + ("max_volume_rescale", 0.42), + ("transfer_types", ["A"]), + ("transfer_types", ["B"]), + ("transfer_types", ["A", "B"]), +] @pytest.mark.serial @@ -39,9 +45,13 @@ def test_valid_construction(device, constructor_args): @pytest.mark.serial @pytest.mark.parametrize("constructor_args", valid_constructor_args) -def test_valid_construction_and_attach(device, simulation_factory, - two_particle_snapshot_factory, - constructor_args, valid_args): +def test_valid_construction_and_attach( + device, + simulation_factory, + two_particle_snapshot_factory, + constructor_args, + valid_args, +): """Test that MuVT can be attached with valid arguments.""" integrator = valid_args[0] args = valid_args[1] @@ -61,7 +71,7 @@ def test_valid_construction_and_attach(device, simulation_factory, muvt = hoomd.hpmc.update.MuVT(**constructor_args) sim = simulation_factory( - two_particle_snapshot_factory(particle_types=['A', 'B'], + two_particle_snapshot_factory(particle_types=["A", "B"], dimensions=n_dimensions, d=2, L=50)) @@ -80,7 +90,7 @@ def test_valid_construction_and_attach(device, simulation_factory, def test_valid_setattr(device, attr, value): """Test that MuVT can get and set attributes.""" muvt = hoomd.hpmc.update.MuVT(trigger=hoomd.trigger.Periodic(10), - transfer_types=['A']) + transfer_types=["A"]) setattr(muvt, attr, value) assert getattr(muvt, attr) == value @@ -108,9 +118,9 @@ def test_valid_setattr_attached(device, attr, value, simulation_factory, mc.shape["B"] = args muvt = hoomd.hpmc.update.MuVT(trigger=hoomd.trigger.Periodic(10), - transfer_types=['A']) + transfer_types=["A"]) sim = simulation_factory( - two_particle_snapshot_factory(particle_types=['A', 'B'], + two_particle_snapshot_factory(particle_types=["A", "B"], dimensions=n_dimensions, d=2, L=50)) @@ -122,68 +132,71 @@ def test_valid_setattr_attached(device, attr, value, simulation_factory, setattr(muvt, attr, value) assert getattr(muvt, attr) == value + def test_insertion_removal(device, simulation_factory, lattice_snapshot_factory): """Test that MuVT is able to insert and remove particles.""" sim = simulation_factory( - lattice_snapshot_factory(particle_types=['A', 'B'], + lattice_snapshot_factory(particle_types=["A", "B"], dimensions=3, a=4, n=7, r=0.1)) mc = hoomd.hpmc.integrate.Sphere(default_d=0.1, default_a=0.1) - mc.shape['A'] = dict(diameter=1.1) - mc.shape['B'] = dict(diameter=1.3) + mc.shape["A"] = dict(diameter=1.1) + mc.shape["B"] = dict(diameter=1.3) sim.operations.integrator = mc muvt = hoomd.hpmc.update.MuVT(trigger=hoomd.trigger.Periodic(5), - transfer_types=['B']) + transfer_types=["B"]) sim.operations.updaters.append(muvt) sim.run(0) # we shouldn't have any particles of type B just yet - assert (muvt.N['B'] == 0) + assert muvt.N["B"] == 0 # and no attempted moves assert sum(muvt.insert_moves) == 0 assert sum(muvt.remove_moves) == 0 # set a positive fugacity - muvt.fugacity['B'] = 1 + muvt.fugacity["B"] = 1 sim.run(20) assert sum(muvt.insert_moves) > 0 assert sum(muvt.remove_moves) > 0 # make a wild guess: there be B particles - assert (muvt.N['B'] > 0) + assert muvt.N["B"] > 0 + -@pytest.mark.skipif(not hoomd.version.llvm_enabled, reason='LLVM not enabled') +@pytest.mark.skipif(not hoomd.version.llvm_enabled, reason="LLVM not enabled") def test_jit_remove_insert(device, simulation_factory, one_particle_snapshot_factory): - """Test that MuVT considers a cpp potential when removing particles by adding a high energy particles that should always be removed by the updater while inserting particles to the top of the box.""" + """Test that MuVT considers cpp potential when removing/adding particles.""" sim = simulation_factory( - one_particle_snapshot_factory(particle_types=['A'], - dimensions=3, - position=(0, 0, -5), - orientation=(1, 0, 0, 0), - L=20) - ) - - sphere_radius = 0.6 + one_particle_snapshot_factory( + particle_types=["A"], + dimensions=3, + position=(-5, 0, 0), + orientation=(1, 0, 0, 0), + L=20, + )) + + sphere_radius = 0.6 mc = hoomd.hpmc.integrate.Sphere(default_d=0.0, default_a=0.0) - mc.shape['A'] = dict(diameter=2*sphere_radius) + mc.shape["A"] = dict(diameter=2 * sphere_radius) - # code returns infinity if the particle center goes past the center of the box + # code returns infinity if the particle center goes past box center - center_wall = ''' - if (r_i.z < 0. || r_i.z > box.getL().z/2) + center_wall = """ + if (r_i.x < 0. || r_i.x > box.getL().x/2) return INFINITY; else return 0.0f; - ''' + """ cpp_external_potential = hoomd.hpmc.external.user.CPPExternalPotential( code=center_wall) @@ -192,111 +205,143 @@ def test_jit_remove_insert(device, simulation_factory, sim.operations.integrator = mc muvt = hoomd.hpmc.update.MuVT(trigger=hoomd.trigger.Periodic(1), - transfer_types=['A']) - muvt.fugacity['A'] = 1e6 + transfer_types=["A"]) + muvt.fugacity["A"] = 1e6 sim.operations.updaters.append(muvt) sim.run(1000) snapshot = sim.state.get_snapshot() - pos = snapshot.particles.position - + + if snapshot.communicator.rank == 0: + + pos = snapshot.particles.position + + # We should have added more than one particle to the box + assert len(pos) > 1 + + # insert to top of box and remove high energy bottom particle + assert numpy.min(pos[:, 0]) >= 0.0 + # We should have inserted particles successfully assert muvt.insert_moves[0] > 0 # We should have successfully removed the high energy particle assert muvt.remove_moves[0] > 0 - # We should have successfully inserted only to the top of the box and removed the high energy particle - assert numpy.min(pos[:,2]) >= 0. - - # We should have added more than one particle to the box - assert len(pos) > 1 def test_plane_wall_insertion(device, simulation_factory, - one_particle_snapshot_factory): + one_particle_snapshot_factory): """Test that MuVT considers a planar wall when inserting particles.""" sim = simulation_factory( - one_particle_snapshot_factory(particle_types=['A'], - dimensions=3, - position=(0, 0, 5), - orientation=(1, 0, 0, 0), - L=20) - ) - - sphere_radius = 0.6 + one_particle_snapshot_factory( + particle_types=["A"], + dimensions=3, + position=(0, 0, 5), + orientation=(1, 0, 0, 0), + L=20, + )) + + sphere_radius = 0.6 mc = hoomd.hpmc.integrate.Sphere(default_d=0.0, default_a=0.0) - mc.shape['A'] = dict(diameter=2*sphere_radius) - walls = [hoomd.wall.Plane(origin=(0,0,0),normal=(0,0,1))] + mc.shape["A"] = dict(diameter=2 * sphere_radius) + walls = [hoomd.wall.Plane(origin=(0, 0, 0), normal=(0, 0, 1))] wall_potential = hoomd.hpmc.external.wall.WallPotential(walls) mc.external_potential = wall_potential sim.operations.integrator = mc muvt = hoomd.hpmc.update.MuVT(trigger=hoomd.trigger.Periodic(1), - transfer_types=['A']) - muvt.fugacity['A'] = 1000 + transfer_types=["A"]) + muvt.fugacity["A"] = 1000 sim.operations.updaters.append(muvt) sim.run(300) + snapshot = sim.state.get_snapshot() - pos = snapshot.particles.position - # Test if inserted spheres are above the plane - assert numpy.min(pos[:,2]) >= sphere_radius - assert len(pos) > 1 + if snapshot.communicator.rank == 0: + pos = snapshot.particles.position + # Test if inserted spheres are above the plane + assert numpy.min(pos[:, 2]) >= sphere_radius + assert len(pos) > 1 + + # We should have inserted particles successfully + assert muvt.insert_moves[0] > 0 + + # We should have successfully attempted some removes + assert sum(muvt.remove_moves) > 0 + def test_spherical_wall_insertion(device, simulation_factory, - one_particle_snapshot_factory): + one_particle_snapshot_factory): """Test that MuVT considers a spherical wall when inserting particles.""" sim = simulation_factory( - one_particle_snapshot_factory(particle_types=['A'], - dimensions=3, - position=(0, 0, 0), - orientation=(1, 0, 0, 0), - L=20) - ) + one_particle_snapshot_factory( + particle_types=["A"], + dimensions=3, + position=(0, 0, 0), + orientation=(1, 0, 0, 0), + L=20, + )) mc = hoomd.hpmc.integrate.Sphere(default_d=0.1, default_a=0.1) - sphere_radius = 0.6 - mc.shape['A'] = dict(diameter=2*sphere_radius) + sphere_radius = 0.6 + mc.shape["A"] = dict(diameter=2 * sphere_radius) walls = [hoomd.wall.Sphere(radius=5)] wall_potential = hoomd.hpmc.external.wall.WallPotential(walls) mc.external_potential = wall_potential sim.operations.integrator = mc muvt = hoomd.hpmc.update.MuVT(trigger=hoomd.trigger.Periodic(1), - transfer_types=['A']) - muvt.fugacity['A'] = 1000 + transfer_types=["A"]) + muvt.fugacity["A"] = 1000 sim.operations.updaters.append(muvt) sim.run(300) snapshot = sim.state.get_snapshot() - pos = snapshot.particles.position - # Test if inserted spheres are inside the spherical wall - assert numpy.max(numpy.linalg.norm(pos,axis=1))-sphere_radius <= 5 - assert len(pos) > 1 + if snapshot.communicator.rank == 0: + pos = snapshot.particles.position + # Test if inserted spheres are inside the spherical wall + assert numpy.max(numpy.linalg.norm(pos, axis=1)) - sphere_radius <= 5 + assert len(pos) > 1 + + # We should have inserted particles successfully + assert muvt.insert_moves[0] > 0 + + # We should have successfully attempted some removes + assert sum(muvt.remove_moves) > 0 + def test_cylindrical_wall_insertion(device, simulation_factory, - one_particle_snapshot_factory): + one_particle_snapshot_factory): """Test that MuVT considers a cylindrical wall when inserting particles.""" sim = simulation_factory( - one_particle_snapshot_factory(particle_types=['A'], - dimensions=3, - position=(0, 0, 0), - orientation=(1, 0, 0, 0), - L=20) - ) - - sphere_radius = 0.6 + one_particle_snapshot_factory( + particle_types=["A"], + dimensions=3, + position=(0, 0, 0), + orientation=(1, 0, 0, 0), + L=20, + )) + + sphere_radius = 0.6 mc = hoomd.hpmc.integrate.Sphere(default_d=0.0, default_a=0.0) - mc.shape['A'] = dict(diameter=2 * sphere_radius) - walls = [hoomd.wall.Cylinder(radius=5,axis=(0,0,1))] + mc.shape["A"] = dict(diameter=2 * sphere_radius) + walls = [hoomd.wall.Cylinder(radius=5, axis=(0, 0, 1))] wall_potential = hoomd.hpmc.external.wall.WallPotential(walls) mc.external_potential = wall_potential sim.operations.integrator = mc muvt = hoomd.hpmc.update.MuVT(trigger=hoomd.trigger.Periodic(1), - transfer_types=['A']) - muvt.fugacity['A'] = 1000 + transfer_types=["A"]) + muvt.fugacity["A"] = 1000 sim.operations.updaters.append(muvt) sim.run(300) snapshot = sim.state.get_snapshot() - pos = snapshot.particles.position - # Test if inserted spheres are inside the cylinder wall - assert numpy.max(numpy.linalg.norm(pos[:,:2],axis=1))-sphere_radius <= 5 - assert len(pos) > 1 + if snapshot.communicator.rank == 0: + pos = snapshot.particles.position + # Test if inserted spheres are inside the cylinder wall + assert numpy.max(numpy.linalg.norm(pos[:, :2], + axis=1)) - sphere_radius <= 5 + assert len(pos) > 1 + + # We should have inserted particles successfully + assert muvt.insert_moves[0] > 0 + + # We should have successfully attempted some removes + assert sum(muvt.remove_moves) > 0 From dc76762b8b920560b0133968618aaec0cac1221c Mon Sep 17 00:00:00 2001 From: Tobias-Dwyer Date: Thu, 1 Feb 2024 18:03:45 -0500 Subject: [PATCH 4/5] syntax changes added back code --- hoomd/hpmc/PatchEnergyJIT.cc | 1 + hoomd/hpmc/UpdaterMuVT.h | 6 ++++-- 2 files changed, 5 insertions(+), 2 deletions(-) diff --git a/hoomd/hpmc/PatchEnergyJIT.cc b/hoomd/hpmc/PatchEnergyJIT.cc index c076188c3a..1713a42fb7 100644 --- a/hoomd/hpmc/PatchEnergyJIT.cc +++ b/hoomd/hpmc/PatchEnergyJIT.cc @@ -69,6 +69,7 @@ void export_PatchEnergyJIT(pybind11::module& m) Scalar, pybind11::array_t>()) .def_property("r_cut", &PatchEnergyJIT::getRCut, &PatchEnergyJIT::setRCut) + .def("energy", &PatchEnergyJIT::energy) .def_property_readonly("param_array", &PatchEnergyJIT::getParamArray); } diff --git a/hoomd/hpmc/UpdaterMuVT.h b/hoomd/hpmc/UpdaterMuVT.h index 7a8721d889..6a7cadf206 100644 --- a/hoomd/hpmc/UpdaterMuVT.h +++ b/hoomd/hpmc/UpdaterMuVT.h @@ -1664,14 +1664,16 @@ bool UpdaterMuVT::tryRemoveParticle(uint64_t timestep, unsigned int tag, const BoxDim box = this->m_pdata->getGlobalBox(); unsigned int type = this->m_pdata->getType(tag); quat orientation(m_pdata->getOrientation(tag)); + Scalar diameter = m_pdata->getDiameter(tag); + Scalar charge = m_pdata->getCharge(tag); if (is_local) { lnboltzmann += field->energy(box, type, pos, quat(orientation), - 1.0, // diameter i - 0.0 // charge i + float(diameter), // diameter i + float(charge) // charge i ); } } From 267e1501f14e80444bd923e5474e56e2d015308e Mon Sep 17 00:00:00 2001 From: Tobias-Dwyer Date: Thu, 1 Feb 2024 20:23:30 -0500 Subject: [PATCH 5/5] add cpu check --- hoomd/hpmc/pytest/test_muvt.py | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/hoomd/hpmc/pytest/test_muvt.py b/hoomd/hpmc/pytest/test_muvt.py index a8b44a7c06..d9aaa99828 100644 --- a/hoomd/hpmc/pytest/test_muvt.py +++ b/hoomd/hpmc/pytest/test_muvt.py @@ -172,6 +172,7 @@ def test_insertion_removal(device, simulation_factory, assert muvt.N["B"] > 0 +@pytest.mark.cpu @pytest.mark.skipif(not hoomd.version.llvm_enabled, reason="LLVM not enabled") def test_jit_remove_insert(device, simulation_factory, one_particle_snapshot_factory): @@ -228,6 +229,7 @@ def test_jit_remove_insert(device, simulation_factory, assert muvt.remove_moves[0] > 0 +@pytest.mark.cpu def test_plane_wall_insertion(device, simulation_factory, one_particle_snapshot_factory): """Test that MuVT considers a planar wall when inserting particles.""" @@ -268,6 +270,7 @@ def test_plane_wall_insertion(device, simulation_factory, assert sum(muvt.remove_moves) > 0 +@pytest.mark.cpu def test_spherical_wall_insertion(device, simulation_factory, one_particle_snapshot_factory): """Test that MuVT considers a spherical wall when inserting particles.""" @@ -307,6 +310,7 @@ def test_spherical_wall_insertion(device, simulation_factory, assert sum(muvt.remove_moves) > 0 +@pytest.mark.cpu def test_cylindrical_wall_insertion(device, simulation_factory, one_particle_snapshot_factory): """Test that MuVT considers a cylindrical wall when inserting particles."""