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..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/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/UpdaterMuVT.h b/hoomd/hpmc/UpdaterMuVT.h index d61f059668..6a7cadf206 100644 --- a/hoomd/hpmc/UpdaterMuVT.h +++ b/hoomd/hpmc/UpdaterMuVT.h @@ -1650,6 +1650,34 @@ 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(); + unsigned int p = m_exec_conf->getPartition() % m_npartition; + + 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)); + Scalar diameter = m_pdata->getDiameter(tag); + Scalar charge = m_pdata->getCharge(tag); + if (is_local) + { + lnboltzmann += field->energy(box, + type, + pos, + quat(orientation), + float(diameter), // diameter i + float(charge) // charge i + ); + } + } + // if not, no overlaps generated if (patch) { @@ -1782,19 +1810,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 @@ -1911,6 +1939,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 +1976,20 @@ bool UpdaterMuVT::tryInsertParticle(uint64_t timestep, ShortReal r_cut_patch(0.0); Scalar r_cut_self(0.0); + unsigned int p = m_exec_conf->getPartition() % m_npartition; + + if (field && (!m_gibbs || p == 0)) + { + 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..d9aaa99828 100644 --- a/hoomd/hpmc/pytest/test_muvt.py +++ b/hoomd/hpmc/pytest/test_muvt.py @@ -5,24 +5,31 @@ import hoomd import pytest +import numpy import hoomd.hpmc.pytest.conftest # 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 @@ -38,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] @@ -60,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)) @@ -79,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 @@ -107,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)) @@ -126,36 +137,215 @@ 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.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): + """Test that MuVT considers cpp potential when removing/adding particles.""" + sim = simulation_factory( + 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) + + # code returns infinity if the particle center goes past box center + + 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) + + 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() + + 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 + + +@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.""" + 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() + 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 + + +@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.""" + 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() + 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 + + +@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.""" + 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() + 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