Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Jump to bottom

Consider Walls when inserting or removing particles for the MuVT updater #1711

Merged
merged 5 commits into from Feb 5, 2024
Merged

Consider Walls when inserting or removing particles for the MuVT updater #1711

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
5 commits
Select commit Hold shift + click to select a range
b01edb1
added muVT wall bug patch
Tobias-Dwyer Jan 29, 2024
d675fa1
added jit wall removal and testing in MuVT
Tobias-Dwyer Jan 30, 2024
3031f6c
add code to run on mpi
Tobias-Dwyer Feb 1, 2024
dc76762
syntax changes added back code
Tobias-Dwyer Feb 1, 2024
267e150
add cpu check
Tobias-Dwyer Feb 2, 2024
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Jump to
Jump to file
Failed to load files.
Diff view
Diff view
20 changes: 20 additions & 0 deletions hoomd/hpmc/ExternalField.h
View file
Open in desktop
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,
tommy-waltmann marked this conversation as resolved.
Show resolved Hide resolved
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
Open in desktop
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
Open in desktop
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,
tommy-waltmann marked this conversation as resolved.
Show resolved Hide resolved
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;
tommy-waltmann marked this conversation as resolved.
Show resolved Hide resolved
}
}

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);
tommy-waltmann marked this conversation as resolved.
Show resolved Hide resolved
}

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
Open in desktop
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