Skip to content

Commit

Permalink
Merge pull request #1711 from glotzerlab/wall_in_MuVT
Browse files Browse the repository at this point in the history 
Consider Walls when inserting or removing particles for the MuVT updater
  • Loading branch information
@joaander
joaander authored Feb 5, 2024
2 parents 57faad6 + 267e150 commit 37f6348
Show file tree
Hide file tree
Showing 5 changed files with 364 additions and 59 deletions.
20 changes: 20 additions & 0 deletions hoomd/hpmc/ExternalField.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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<Scalar>& r_i,
const quat<Scalar>& q_i,
Scalar diameter,
Scalar charge)
{
return 0;
}

virtual bool hasVolume()
{
return false;
Expand Down Expand Up @@ -88,6 +107,7 @@ template<class Shape> void export_ExternalFieldInterface(pybind11::module& m, st
.def(pybind11::init<std::shared_ptr<SystemDefinition>>())
.def("compute", &ExternalFieldMono<Shape>::compute)
.def("energydiff", &ExternalFieldMono<Shape>::energydiff)
.def("energy", &ExternalFieldMono<Shape>::energy)
.def("calculateDeltaE", &ExternalFieldMono<Shape>::calculateDeltaE);
}

Expand Down
48 changes: 24 additions & 24 deletions hoomd/hpmc/ExternalFieldJIT.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -66,6 +66,16 @@ template<class Shape> class ExternalFieldJIT : public hpmc::ExternalFieldMono<Sh
m_factory = std::shared_ptr<ExternalFieldEvalFactory>(factory);
}

float energy_no_wrap(const BoxDim& box,
unsigned int type,
const vec3<Scalar>& r_i,
const quat<Scalar>& 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.
Expand All @@ -82,7 +92,10 @@ template<class Shape> class ExternalFieldJIT : public hpmc::ExternalFieldMono<Sh
Scalar diameter,
Scalar charge)
{
return m_eval(box, type, r_i, q_i, diameter, charge);
vec3<Scalar> r_i_wrapped = r_i - vec3<Scalar>(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
Expand All @@ -109,13 +122,10 @@ template<class Shape> class ExternalFieldJIT : public hpmc::ExternalFieldMono<Sh
// read in the current position and orientation
Scalar4 postype_i = h_postype.data[i];
unsigned int typ_i = __scalar_as_int(postype_i.w);
vec3<Scalar> pos_i = vec3<Scalar>(postype_i) - vec3<Scalar>(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<Scalar>(postype_i),
quat<Scalar>(h_orientation.data[i]),
h_diameter.data[i],
h_charge.data[i]);
Expand Down Expand Up @@ -171,23 +181,20 @@ template<class Shape> class ExternalFieldJIT : public hpmc::ExternalFieldMono<Sh
Scalar4 postype_i = h_postype.data[i];
unsigned int typ_i = __scalar_as_int(postype_i.w);
int3 image = make_int3(0, 0, 0);
vec3<Scalar> pos_i = vec3<Scalar>(postype_i) - vec3<Scalar>(this->m_pdata->getOrigin());
box_new.wrap(pos_i, image);
image = make_int3(0, 0, 0);
vec3<Scalar> old_pos_i = vec3<Scalar>(*(position_old + i)) - vec3<Scalar>(origin_old);
box_old.wrap(old_pos_i, image);
dE += energy(box_new,
typ_i,
pos_i,
vec3<Scalar>(postype_i),
quat<Scalar>(h_orientation.data[i]),
h_diameter.data[i],
h_charge.data[i]);
dE -= energy(box_old,
typ_i,
old_pos_i,
quat<Scalar>(*(orientation_old + i)),
h_diameter.data[i],
h_charge.data[i]);
dE -= energy_no_wrap(box_old,
typ_i,
old_pos_i,
quat<Scalar>(*(orientation_old + i)),
h_diameter.data[i],
h_charge.data[i]);
}
#ifdef ENABLE_MPI
if (this->m_sysdef->isDomainDecomposed())
Expand Down Expand Up @@ -223,24 +230,17 @@ template<class Shape> class ExternalFieldJIT : public hpmc::ExternalFieldMono<Sh
ArrayHandle<Scalar> 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<Scalar> pos_new_shifted = position_new - vec3<Scalar>(this->m_pdata->getOrigin());
vec3<Scalar> pos_old_shifted = position_old - vec3<Scalar>(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]);
Expand Down
50 changes: 50 additions & 0 deletions hoomd/hpmc/ExternalFieldWall.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -629,6 +629,56 @@ template<class Shape> class ExternalFieldWall : public ExternalFieldMono<Shape>
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<Scalar>& r_i,
const quat<Scalar>& q_i,
Scalar diameter,
Scalar charge)
{
const std::vector<typename Shape::param_type,
hoomd::detail::managed_allocator<typename Shape::param_type>>& params
= m_mc->getParams();
Shape shape(q_i, params[type]);
vec3<Scalar> 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,
Expand Down
65 changes: 55 additions & 10 deletions hoomd/hpmc/UpdaterMuVT.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -1650,6 +1650,34 @@ bool UpdaterMuVT<Shape>::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<Scalar> pos(p);
const BoxDim box = this->m_pdata->getGlobalBox();
unsigned int type = this->m_pdata->getType(tag);
quat<Scalar> 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<float>(orientation),
float(diameter), // diameter i
float(charge) // charge i
);
}
}

// if not, no overlaps generated
if (patch)
{
Expand Down Expand Up @@ -1782,19 +1810,19 @@ bool UpdaterMuVT<Shape>::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
Expand Down Expand Up @@ -1911,6 +1939,9 @@ bool UpdaterMuVT<Shape>::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;
Expand Down Expand Up @@ -1945,6 +1976,20 @@ bool UpdaterMuVT<Shape>::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<float>(orientation),
1.0, // diameter i
0.0 // charge i
);
}

if (patch)
{
r_cut_patch = ShortReal(patch->getRCut() + 0.5 * patch->getAdditiveCutoff(type));
Expand Down
Loading

0 comments on commit 37f6348

Please sign in to comment.