Skip to content

Commit 37f6348

Browse files
joaanderjoaander
authored
Merge pull request #1711 from glotzerlab/wall_in_MuVT
Consider Walls when inserting or removing particles for the MuVT updater
2 parents 57faad6 + 267e150 commit 37f6348

File tree

5 files changed

+364
-59
lines changed

5 files changed

+364
-59
lines changed

‎hoomd/hpmc/ExternalField.h

Lines changed: 20 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -39,6 +39,25 @@ class ExternalField : public Compute
3939
return 0;
4040
}
4141

42+
//! Evaluate the energy of the force.
43+
/*! \param box The system box.
44+
\param type Particle type.
45+
\param r_i Particle position
46+
\param q_i Particle orientation.
47+
\param diameter Particle diameter.
48+
\param charge Particle charge.
49+
\returns Energy due to the force
50+
*/
51+
virtual float energy(const BoxDim& box,
52+
unsigned int type,
53+
const vec3<Scalar>& r_i,
54+
const quat<Scalar>& q_i,
55+
Scalar diameter,
56+
Scalar charge)
57+
{
58+
return 0;
59+
}
60+
4261
virtual bool hasVolume()
4362
{
4463
return false;
@@ -88,6 +107,7 @@ template<class Shape> void export_ExternalFieldInterface(pybind11::module& m, st
88107
.def(pybind11::init<std::shared_ptr<SystemDefinition>>())
89108
.def("compute", &ExternalFieldMono<Shape>::compute)
90109
.def("energydiff", &ExternalFieldMono<Shape>::energydiff)
110+
.def("energy", &ExternalFieldMono<Shape>::energy)
91111
.def("calculateDeltaE", &ExternalFieldMono<Shape>::calculateDeltaE);
92112
}
93113

‎hoomd/hpmc/ExternalFieldJIT.h

Lines changed: 24 additions & 24 deletions
Original file line numberDiff line numberDiff line change
@@ -66,6 +66,16 @@ template<class Shape> class ExternalFieldJIT : public hpmc::ExternalFieldMono<Sh
6666
m_factory = std::shared_ptr<ExternalFieldEvalFactory>(factory);
6767
}
6868

69+
float energy_no_wrap(const BoxDim& box,
70+
unsigned int type,
71+
const vec3<Scalar>& r_i,
72+
const quat<Scalar>& q_i,
73+
Scalar diameter,
74+
Scalar charge)
75+
{
76+
return m_eval(box, type, r_i, q_i, diameter, charge);
77+
}
78+
6979
//! Evaluate the energy of the force.
7080
/*! \param box The system box.
7181
\param type Particle type.
@@ -82,7 +92,10 @@ template<class Shape> class ExternalFieldJIT : public hpmc::ExternalFieldMono<Sh
8292
Scalar diameter,
8393
Scalar charge)
8494
{
85-
return m_eval(box, type, r_i, q_i, diameter, charge);
95+
vec3<Scalar> r_i_wrapped = r_i - vec3<Scalar>(this->m_pdata->getOrigin());
96+
int3 image = make_int3(0, 0, 0);
97+
box.wrap(r_i_wrapped, image);
98+
return energy_no_wrap(box, type, r_i_wrapped, q_i, diameter, charge);
8699
}
87100

88101
//! Computes the total field energy of the system at the current state
@@ -109,13 +122,10 @@ template<class Shape> class ExternalFieldJIT : public hpmc::ExternalFieldMono<Sh
109122
// read in the current position and orientation
110123
Scalar4 postype_i = h_postype.data[i];
111124
unsigned int typ_i = __scalar_as_int(postype_i.w);
112-
vec3<Scalar> pos_i = vec3<Scalar>(postype_i) - vec3<Scalar>(this->m_pdata->getOrigin());
113-
int3 image = make_int3(0, 0, 0);
114-
box.wrap(pos_i, image);
115125

116126
total_energy += energy(box,
117127
typ_i,
118-
pos_i,
128+
vec3<Scalar>(postype_i),
119129
quat<Scalar>(h_orientation.data[i]),
120130
h_diameter.data[i],
121131
h_charge.data[i]);
@@ -171,23 +181,20 @@ template<class Shape> class ExternalFieldJIT : public hpmc::ExternalFieldMono<Sh
171181
Scalar4 postype_i = h_postype.data[i];
172182
unsigned int typ_i = __scalar_as_int(postype_i.w);
173183
int3 image = make_int3(0, 0, 0);
174-
vec3<Scalar> pos_i = vec3<Scalar>(postype_i) - vec3<Scalar>(this->m_pdata->getOrigin());
175-
box_new.wrap(pos_i, image);
176-
image = make_int3(0, 0, 0);
177184
vec3<Scalar> old_pos_i = vec3<Scalar>(*(position_old + i)) - vec3<Scalar>(origin_old);
178185
box_old.wrap(old_pos_i, image);
179186
dE += energy(box_new,
180187
typ_i,
181-
pos_i,
188+
vec3<Scalar>(postype_i),
182189
quat<Scalar>(h_orientation.data[i]),
183190
h_diameter.data[i],
184191
h_charge.data[i]);
185-
dE -= energy(box_old,
186-
typ_i,
187-
old_pos_i,
188-
quat<Scalar>(*(orientation_old + i)),
189-
h_diameter.data[i],
190-
h_charge.data[i]);
192+
dE -= energy_no_wrap(box_old,
193+
typ_i,
194+
old_pos_i,
195+
quat<Scalar>(*(orientation_old + i)),
196+
h_diameter.data[i],
197+
h_charge.data[i]);
191198
}
192199
#ifdef ENABLE_MPI
193200
if (this->m_sysdef->isDomainDecomposed())
@@ -223,24 +230,17 @@ template<class Shape> class ExternalFieldJIT : public hpmc::ExternalFieldMono<Sh
223230
ArrayHandle<Scalar> h_charge(this->m_pdata->getCharges(),
224231
access_location::host,
225232
access_mode::read);
226-
const BoxDim box = this->m_pdata->getGlobalBox();
227-
int3 image = make_int3(0, 0, 0);
228-
vec3<Scalar> pos_new_shifted = position_new - vec3<Scalar>(this->m_pdata->getOrigin());
229-
vec3<Scalar> pos_old_shifted = position_old - vec3<Scalar>(this->m_pdata->getOrigin());
230-
box.wrap(pos_new_shifted, image);
231-
image = make_int3(0, 0, 0);
232-
box.wrap(pos_old_shifted, image);
233233

234234
double dE = 0.0;
235235
dE += energy(this->m_pdata->getGlobalBox(),
236236
typ_i,
237-
pos_new_shifted,
237+
position_new,
238238
shape_new.orientation,
239239
h_diameter.data[index],
240240
h_charge.data[index]);
241241
dE -= energy(this->m_pdata->getGlobalBox(),
242242
typ_i,
243-
pos_old_shifted,
243+
position_old,
244244
shape_old.orientation,
245245
h_diameter.data[index],
246246
h_charge.data[index]);

‎hoomd/hpmc/ExternalFieldWall.h

Lines changed: 50 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -629,6 +629,56 @@ template<class Shape> class ExternalFieldWall : public ExternalFieldMono<Shape>
629629
return double(0.0);
630630
}
631631

632+
//! Evaluate the energy of the force.
633+
/*! \param box The system box.
634+
\param type Particle type.
635+
\param r_i Particle position
636+
\param q_i Particle orientation.
637+
\param diameter Particle diameter.
638+
\param charge Particle charge.
639+
\returns Energy due to the force
640+
*/
641+
virtual float energy(const BoxDim& box,
642+
unsigned int type,
643+
const vec3<Scalar>& r_i,
644+
const quat<Scalar>& q_i,
645+
Scalar diameter,
646+
Scalar charge)
647+
{
648+
const std::vector<typename Shape::param_type,
649+
hoomd::detail::managed_allocator<typename Shape::param_type>>& params
650+
= m_mc->getParams();
651+
Shape shape(q_i, params[type]);
652+
vec3<Scalar> origin(m_pdata->getOrigin());
653+
654+
for (size_t i = 0; i < m_Spheres.size(); i++)
655+
{
656+
if (!test_confined(m_Spheres[i], shape, r_i, origin, box))
657+
{
658+
return INFINITY;
659+
}
660+
}
661+
662+
for (size_t i = 0; i < m_Cylinders.size(); i++)
663+
{
664+
set_cylinder_wall_verts(m_Cylinders[i], shape);
665+
if (!test_confined(m_Cylinders[i], shape, r_i, origin, box))
666+
{
667+
return INFINITY;
668+
}
669+
}
670+
671+
for (size_t i = 0; i < m_Planes.size(); i++)
672+
{
673+
if (!test_confined(m_Planes[i], shape, r_i, origin, box))
674+
{
675+
return INFINITY;
676+
}
677+
}
678+
679+
return double(0.0);
680+
}
681+
632682
double calculateDeltaE(uint64_t timestep,
633683
const Scalar4* const position_old,
634684
const Scalar4* const orientation_old,

‎hoomd/hpmc/UpdaterMuVT.h

Lines changed: 55 additions & 10 deletions
Original file line numberDiff line numberDiff line change
@@ -1650,6 +1650,34 @@ bool UpdaterMuVT<Shape>::tryRemoveParticle(uint64_t timestep, unsigned int tag,
16501650
// do we have to compute energetic contribution?
16511651
auto patch = m_mc->getPatchEnergy();
16521652

1653+
// do we have to compute a wall contribution?
1654+
auto field = m_mc->getExternalField();
1655+
unsigned int p = m_exec_conf->getPartition() % m_npartition;
1656+
1657+
if (field && (!m_gibbs || p == 0))
1658+
{
1659+
// getPosition() takes into account grid shift, correct for that
1660+
Scalar3 p = m_pdata->getPosition(tag) + m_pdata->getOrigin();
1661+
int3 tmp = make_int3(0, 0, 0);
1662+
m_pdata->getGlobalBox().wrap(p, tmp);
1663+
vec3<Scalar> pos(p);
1664+
const BoxDim box = this->m_pdata->getGlobalBox();
1665+
unsigned int type = this->m_pdata->getType(tag);
1666+
quat<Scalar> orientation(m_pdata->getOrientation(tag));
1667+
Scalar diameter = m_pdata->getDiameter(tag);
1668+
Scalar charge = m_pdata->getCharge(tag);
1669+
if (is_local)
1670+
{
1671+
lnboltzmann += field->energy(box,
1672+
type,
1673+
pos,
1674+
quat<float>(orientation),
1675+
float(diameter), // diameter i
1676+
float(charge) // charge i
1677+
);
1678+
}
1679+
}
1680+
16531681
// if not, no overlaps generated
16541682
if (patch)
16551683
{
@@ -1782,19 +1810,19 @@ bool UpdaterMuVT<Shape>::tryRemoveParticle(uint64_t timestep, unsigned int tag,
17821810
} // end loop over AABB nodes
17831811
} // end loop over images
17841812
}
1813+
}
17851814

17861815
#ifdef ENABLE_MPI
1787-
if (m_sysdef->isDomainDecomposed())
1788-
{
1789-
MPI_Allreduce(MPI_IN_PLACE,
1790-
&lnboltzmann,
1791-
1,
1792-
MPI_HOOMD_SCALAR,
1793-
MPI_SUM,
1794-
m_exec_conf->getMPICommunicator());
1795-
}
1796-
#endif
1816+
if (m_sysdef->isDomainDecomposed())
1817+
{
1818+
MPI_Allreduce(MPI_IN_PLACE,
1819+
&lnboltzmann,
1820+
1,
1821+
MPI_HOOMD_SCALAR,
1822+
MPI_SUM,
1823+
m_exec_conf->getMPICommunicator());
17971824
}
1825+
#endif
17981826
}
17991827

18001828
// Depletants
@@ -1911,6 +1939,9 @@ bool UpdaterMuVT<Shape>::tryInsertParticle(uint64_t timestep,
19111939
// do we have to compute energetic contribution?
19121940
auto patch = m_mc->getPatchEnergy();
19131941

1942+
// do we have to compute a wall contribution?
1943+
auto field = m_mc->getExternalField();
1944+
19141945
lnboltzmann = Scalar(0.0);
19151946

19161947
unsigned int overlap = 0;
@@ -1945,6 +1976,20 @@ bool UpdaterMuVT<Shape>::tryInsertParticle(uint64_t timestep,
19451976
ShortReal r_cut_patch(0.0);
19461977
Scalar r_cut_self(0.0);
19471978

1979+
unsigned int p = m_exec_conf->getPartition() % m_npartition;
1980+
1981+
if (field && (!m_gibbs || p == 0))
1982+
{
1983+
const BoxDim& box = this->m_pdata->getGlobalBox();
1984+
lnboltzmann -= field->energy(box,
1985+
type,
1986+
pos,
1987+
quat<float>(orientation),
1988+
1.0, // diameter i
1989+
0.0 // charge i
1990+
);
1991+
}
1992+
19481993
if (patch)
19491994
{
19501995
r_cut_patch = ShortReal(patch->getRCut() + 0.5 * patch->getAdditiveCutoff(type));

0 commit comments

Comments
 (0)