Compute velocity of MPCD particles (#65)

* Add GroupVelocityCompute class

* Fix compilation errors

* Add python unit tests

* Add test for repeated suffix

* Link analyzer into docs

* Clarify documentation

* Add MPCD velocity compute

* Fix copy errors and tests

* Fix MPI compile error

* Link in class to docs

* Address review comments

azplugins/CMakeLists.txt

@@ -121,6 +121,7 @@ list(APPEND py_files
 # mpcd c++ files
 list(APPEND _${COMPONENT_NAME}_sources
     MPCDReversePerturbationFlow.cc
+    MPCDVelocityCompute.cc
     SinusoidalChannelFiller.cc
     SinusoidalExpansionConstrictionFiller.cc
     )

azplugins/MPCDVelocityCompute.cc

@@ -0,0 +1,117 @@
+// Copyright (c) 2018-2020, Michael P. Howard
+// Copyright (c) 2021-2022, Auburn University
+// This file is part of the azplugins project, released under the Modified BSD License.
+
+/*!
+ * \file MPCDVelocityCompute.cc
+ * \brief Definition of MPCDVelocityCompute
+ */
+
+#include "MPCDVelocityCompute.h"
+
+namespace azplugins
+{
+
+/*!
+ * \param sysdata MPCD system data
+ * \param suffix Suffix to attach to logged quantities
+ */
+MPCDVelocityCompute::MPCDVelocityCompute(std::shared_ptr<mpcd::SystemData> sysdata,
+                                         const std::string& suffix)
+    : Compute(sysdata->getSystemDefinition()),
+      m_mpcd_pdata(sysdata->getParticleData()),
+      m_lognames{"mpcd_vx"+suffix,"mpcd_vy"+suffix,"mpcd_vz"+suffix},
+      m_velocity(make_scalar3(0,0,0))
+    {
+    }
+
+MPCDVelocityCompute::~MPCDVelocityCompute()
+    {
+    }
+
+/*!
+ * \param timestep Simulation timestep
+ *
+ * The center-of-mass velocity of the MPCD particles is determined by first summing the
+ * momentum and mass of all particles, then dividing momentum by mass. This
+ * compute supports MPI decomposition, and the result is available on all ranks.
+ */
+void MPCDVelocityCompute::compute(unsigned int timestep)
+    {
+    if (!shouldCompute(timestep))
+        return;
+
+    // empty particle data has no velocity, quit early
+    if (m_mpcd_pdata->getNGlobal() == 0)
+        {
+        m_velocity = make_scalar3(0,0,0);
+        return;
+        }
+
+    // accumulate momentum and masses
+    // this could be simplified by assuming mass is the same, but to make intention clear in
+    // case particles can have different masses in future, carrying out the explicit calculation
+    const unsigned int N = m_mpcd_pdata->getN();
+    ArrayHandle<Scalar4> h_vel(m_mpcd_pdata->getVelocities(), access_location::host, access_mode::read);
+    const Scalar m = m_mpcd_pdata->getMass();
+    Scalar3 momentum = make_scalar3(0,0,0);
+    Scalar mass(0);
+    for (unsigned int i=0; i < N; ++i)
+        {
+        const Scalar4 vel_cell = h_vel.data[i];
+        const Scalar3 vel = make_scalar3(vel_cell.x,vel_cell.y,vel_cell.z);
+        momentum += m*vel;
+        mass += m;
+        }
+    #ifdef ENABLE_MPI
+    if (m_comm)
+        {
+        Scalar buffer[4] = {momentum.x, momentum.y, momentum.z, mass};
+        MPI_Allreduce(MPI_IN_PLACE, &buffer[0], 4, MPI_HOOMD_SCALAR, MPI_SUM, m_exec_conf->getMPICommunicator());
+        momentum = make_scalar3(buffer[0],buffer[1],buffer[2]);
+        mass = buffer[3];
+        }
+    #endif // ENABLE_MPI
+
+    // reduce total momentum by mass to get center-of-mass
+    m_velocity = momentum/mass;
+    }
+
+std::vector<std::string> MPCDVelocityCompute::getProvidedLogQuantities()
+    {
+    return m_lognames;
+    }
+
+Scalar MPCDVelocityCompute::getLogValue(const std::string& quantity, unsigned int timestep)
+    {
+    compute(timestep);
+    if (quantity == m_lognames[0])
+        {
+        return m_velocity.x;
+        }
+    else if (quantity == m_lognames[1])
+        {
+        return m_velocity.y;
+        }
+    else if (quantity == m_lognames[2])
+        {
+        return m_velocity.z;
+        }
+    else
+        {
+        m_exec_conf->msg->error() << "MPCDVelocityCompute: " << quantity << " is not a valid log quantity" << std::endl;
+        throw std::runtime_error("Unknown log quantity");
+        }
+    }
+
+namespace detail
+{
+void export_MPCDVelocityCompute(pybind11::module& m)
+    {
+    namespace py = pybind11;
+    py::class_<MPCDVelocityCompute,std::shared_ptr<MPCDVelocityCompute>>(m, "MPCDVelocityCompute", py::base<Compute>())
+        .def(py::init<std::shared_ptr<mpcd::SystemData>,const std::string&>())
+        ;
+    }
+} // end namespace detail
+} // end namespace azplugins

azplugins/MPCDVelocityCompute.h

@@ -0,0 +1,61 @@
+// Copyright (c) 2018-2020, Michael P. Howard
+// Copyright (c) 2021-2022, Auburn University
+// This file is part of the azplugins project, released under the Modified BSD License.
+
+/*!
+ * \file MPCDVelocityCompute.h
+ * \brief Declaration of MPCDVelocityCompute
+ */
+
+#ifndef AZPLUGINS_MPCD_VELOCITY_COMPUTE_H_
+#define AZPLUGINS_MPCD_VELOCITY_COMPUTE_H_
+
+#ifdef NVCC
+#error This header cannot be compiled by nvcc
+#endif
+
+#include "hoomd/Compute.h"
+#include "hoomd/HOOMDMath.h"
+#include "hoomd/mpcd/ParticleData.h"
+#include "hoomd/mpcd/SystemData.h"
+
+#include "hoomd/extern/pybind/include/pybind11/pybind11.h"
+#include <string>
+#include <vector>
+
+namespace azplugins
+{
+//! Compute the center-of-mass velocity of MPCD particles
+class PYBIND11_EXPORT MPCDVelocityCompute : public Compute
+    {
+    public:
+        //! Constructor
+        MPCDVelocityCompute(std::shared_ptr<mpcd::SystemData> sysdata,
+                            const std::string& suffix);
+
+        //! Destructor
+        ~MPCDVelocityCompute();
+
+        //! Compute center-of-mass velocity of particles
+        void compute(unsigned int timestep) override;
+
+        //! List of logged quantities
+        std::vector<std::string> getProvidedLogQuantities() override;
+
+        //! Return the logged value
+        Scalar getLogValue(const std::string& quantity, unsigned int timestep) override;
+
+    protected:
+        std::shared_ptr<mpcd::ParticleData> m_mpcd_pdata;   //!< Particle group
+        std::vector<std::string> m_lognames;    //!< Logged quantities
+        Scalar3 m_velocity;                     //!< Last compute velocity
+    };
+
+namespace detail
+{
+//! Exports the MPCDVelocityCompute to python
+void export_MPCDVelocityCompute(pybind11::module& m);
+} // end namespace detail
+} // end namespace azplugins
+
+#endif // AZPLUGINS_MPCD_VELOCITY_COMPUTE_H_

azplugins/module.cc

@@ -31,6 +31,7 @@ namespace py = pybind11;
 /* MPCD */
 #ifdef ENABLE_MPCD
 #include "MPCDReversePerturbationFlow.h"
+#include "MPCDVelocityCompute.h"
 #include "SinusoidalChannelFiller.h"
 #include "SinusoidalExpansionConstrictionFiller.h"
 #include "hoomd/mpcd/ConfinedStreamingMethod.h"
@@ -265,6 +266,7 @@ PYBIND11_MODULE(_azplugins, m)
     mpcd::detail::export_ConfinedStreamingMethod<azplugins::detail::SinusoidalChannel>(m);
     mpcd::detail::export_ConfinedStreamingMethod<azplugins::detail::SinusoidalExpansionConstriction>(m);
     azplugins::detail::export_MPCDReversePerturbationFlow(m);
+    azplugins::detail::export_MPCDVelocityCompute(m);
     #ifdef ENABLE_CUDA
     azplugins::detail::export_SinusoidalChannelFillerGPU(m);
     azplugins::detail::export_SinusoidalExpansionConstrictionFillerGPU(m);

azplugins/mpcd.py

@@ -8,10 +8,12 @@
 .. autosummary::
     :nosignatures:
 
+    compute_velocity
     reverse_perturbation
     sinusoidal_channel
     sinusoidal_expansion_constriction
 
+.. autoclass:: compute_velocity
 .. autoclass:: reverse_perturbation
 .. autoclass:: sinusoidal_channel
 .. autoclass:: sinusoidal_expansion_constriction
@@ -23,6 +25,59 @@
 
 from . import _azplugins
 
+class compute_velocity(hoomd.compute._compute):
+    r"""Compute center-of-mass velocity of MPCD particles
+
+    Args:
+        suffix (str): Suffix to attach to logged quantities.
+
+    This computes the center-of-mass velocity of all MPCD particles:
+
+    .. math::
+
+        \mathbf{v}_{\rm cm} = \dfrac{\sum_i m_i \mathbf{v}_i}{\sum_i m_i}
+
+    where :math:`\mathbf{v}_i` is the velocity and and :math:`m_i` is the mass
+    of particle *i* in the group. Note that because all MPCD particles currently
+    have the same mass, this is equivalent to the number-average velocity.
+
+    The components of the result are exposed as loggable quantities ``mpcd_vx``,
+    ``mpcd_vy``, and ``mpcd_vz`` with ``suffix`` appended. By default,
+    ``suffix`` is an empty string, but a custom suffix may be specified if needed
+    to distinguish the logged quantity from something else. You can save these
+    results using :py:class:`hoomd.analyze.log`.
+
+    Example::
+
+        # without suffix
+        azplugins.mpcd.compute_velocity()
+        hoomd.analyze.log(filename='velocity.dat', quantities=['mpcd_vx'], period=10)
+
+        # with suffix
+        azplugins.mpcd.compute_velocity(suffix='-srd')
+        hoomd.analyze.log(filename='velocity_srd.dat', quantities=['mpcd_vx-srd'], period=100)
+
+    """
+    def __init__(self, suffix=None):
+        hoomd.util.print_status_line()
+        super().__init__()
+
+        # create suffix for logged quantity
+        if suffix is None:
+            suffix = ''
+
+        if suffix in self._used_suffixes:
+            hoomd.context.msg.error('azplugins.mpcd.velocity: Suffix {} already used\n'.format(suffix))
+            raise ValueError('Suffix {} already used for MPCD velocity'.format(suffix))
+        else:
+            self._used_suffixes.append(suffix)
+
+        # create the c++ mirror class
+        self.cpp_compute = _azplugins.MPCDVelocityCompute(hoomd.context.current.mpcd.data, suffix)
+        hoomd.context.current.system.addCompute(self.cpp_compute, self.compute_name);
+
+    _used_suffixes = []
+
 class reverse_perturbation(hoomd.update._updater):
     R"""Reverse nonequilibrium shear flow in MPCD simulations.
 

azplugins/test-py/CMakeLists.txt

@@ -79,6 +79,7 @@ set(MPCD_ONLY
     test_mpcd_reverse_perturbation
     test_mpcd_sinusoidal_channel
     test_mpcd_sinusoidal_expansion_constriction
+    test_mpcd_velocity_compute
     )
 
 set(MPI_ONLY

azplugins/test-py/test_mpcd_velocity_compute.py

@@ -0,0 +1,70 @@
+# Copyright (c) 2018-2020, Michael P. Howard
+# Copyright (c) 2021-2022, Auburn University
+# This file is part of the azplugins project, released under the Modified BSD License.
+
+import hoomd
+from hoomd import md
+from hoomd import mpcd
+hoomd.context.initialize()
+try:
+    from hoomd import azplugins
+except ImportError:
+    import azplugins
+import azplugins.mpcd
+import unittest
+import numpy as np
+
+class analyze_group_velocity_tests(unittest.TestCase):
+    def setUp(self):
+        # empty MD system
+        md_snap = hoomd.data.make_snapshot(N=0, box=hoomd.data.boxdim(L=10))
+        hoomd.init.read_snapshot(md_snap)
+
+        # MPCD system
+        snap = mpcd.data.make_snapshot(N=2)
+        if hoomd.comm.get_rank() == 0:
+            snap.particles.position[:] = [[0,0,-5],[0,0,5]]
+            snap.particles.velocity[:] = [[1,-2,3],[-3,6,-9]]
+        self.s = mpcd.init.read_snapshot(snap)
+
+    def test_compute(self):
+        azplugins.mpcd.compute_velocity()
+        log = hoomd.analyze.log(
+            filename=None,
+            quantities=['mpcd_vx','mpcd_vy','mpcd_vz'],
+            period=1
+            )
+        hoomd.run(1)
+        v = [
+            log.query('mpcd_vx'),
+            log.query('mpcd_vy'),
+            log.query('mpcd_vz')
+            ]
+        np.testing.assert_allclose(v, [-1,2,-3])
+
+    def test_compute_suffix(self):
+        azplugins.mpcd.compute_velocity(suffix='_foo')
+        log = hoomd.analyze.log(
+            filename=None,
+            quantities=['mpcd_vx_foo','mpcd_vy_foo','mpcd_vz_foo'],
+            period=1
+            )
+        hoomd.run(1)
+        v = [
+            log.query('mpcd_vx_foo'),
+            log.query('mpcd_vy_foo'),
+            log.query('mpcd_vz_foo')
+            ]
+        np.testing.assert_allclose(v, [-1,2,-3])
+
+    def test_unique_suffix(self):
+        azplugins.mpcd.compute_velocity(suffix='_1')
+        with self.assertRaises(ValueError):
+            azplugins.mpcd.compute_velocity(suffix='_1')
+
+    def tearDown(self):
+        del self.s
+        hoomd.context.initialize()
+
+if __name__ == '__main__':
+    unittest.main(argv = ['test.py', '-v'])