Merge pull request #1003 from glotzerlab/active_constraint

Active forces with manifold constraints

hoomd/Integrator.cc

@@ -62,9 +62,8 @@ void Integrator::removeHalfStepHook()
  */
 void Integrator::setDeltaT(Scalar deltaT)
     {
-    if (m_deltaT <= 0.0)
-        m_exec_conf->msg->warning()
-            << "integrate.*: A timestep of less than 0.0 was specified" << endl;
+    if (m_deltaT < 0.0)
+        throw std::domain_error("delta_t must be positive");
 
     for (auto& force : m_forces)
         {
@@ -883,6 +882,17 @@ void Integrator::computeNetForceGPU(uint64_t timestep)
 void Integrator::update(uint64_t timestep)
     {
     Updater::update(timestep);
+
+    // ensure that the force computes know the current step size
+    for (auto& force : m_forces)
+        {
+        force->setDeltaT(m_deltaT);
+        }
+
+    for (auto& constraint_force : m_constraint_forces)
+        {
+        constraint_force->setDeltaT(m_deltaT);
+        }
     }
 
 /** prepRun() is to be called at the very beginning of each run, before any analyzers are called,

hoomd/ParticleData.cc

@@ -2432,7 +2432,7 @@ unsigned int ParticleData::addParticle(unsigned int type)
         {
         // update reverse-lookup table
         ArrayHandle<unsigned int> h_rtag(m_rtag, access_location::host, access_mode::readwrite);
-        assert(h_rtag.data[tag] = NOT_LOCAL);
+        assert((h_rtag.data[tag] = NOT_LOCAL));
         if (m_exec_conf->getRank() == 0)
             {
             // we add the particle at the end

hoomd/VectorMath.h

@@ -34,10 +34,6 @@
    easier. These include basic element-wise addition, subtraction, division, and multiplication (and
    += -= *= /=), and similarly division, and multiplication by scalars, and negation. The dot and
    cross product are also defined.
-
-    \note Operations like normalize, length, etc... are purposefully **NOT** defined. These hide a
-   sqrt. In high performance code, it is good for the programmer to have to explicitly call
-   expensive operations.
 */
 template<class Real> struct vec3
     {
@@ -320,6 +316,17 @@ template<class Real> DEVICE inline vec3<Real> cross(const vec3<Real>& a, const v
     return vec3<Real>(a.y * b.z - a.z * b.y, a.z * b.x - a.x * b.z, a.x * b.y - a.y * b.x);
     }
 
+/// Normalize the vector
+/*! \param a Vector
+
+    \returns A normal vector in the direction of *a*.
+*/
+template<class Real> DEVICE inline vec3<Real> normalize(const vec3<Real>& a)
+    {
+    Real inverse_norm = fast::rsqrt(dot(a, a));
+    return a * inverse_norm;
+    }
+
 //! Convenience function for converting a vec3 to a Scalar3
 DEVICE inline Scalar3 vec_to_scalar3(const vec3<Scalar>& a)
     {
@@ -334,18 +341,14 @@ DEVICE inline Scalar4 vec_to_scalar4(const vec3<Scalar>& a, Scalar w)
 
 /////////////////////////////// vec2 ///////////////////////////////////
 
-//! 3 element vector
+//! 2 element vector
 /*! \tparam Real Data type of the components
 
     vec2 defines a simple 2 element vector. The components are available publicly as .x and .y.
    Along with vec2, a number of simple operations are defined to make writing vector math code
    easier. These include basic element-wise addition, subtraction, division, and multiplication (and
    += -= *= /=), and similarly division, and multiplication by scalars, and negation. The dot
    product is also defined.
-
-    \note Operations like normalize, length, etc... are purposefully **NOT** defined. These hide a
-   sqrt. In high performance code, it is good for the programmer to have to explicitly call
-   expensive operations.
 */
 template<class Real> struct vec2
     {
@@ -610,6 +613,17 @@ template<class Real> DEVICE inline Real perpdot(const vec2<Real>& a, const vec2<
     return dot(perp(a), b);
     }
 
+/// Normalize the vector
+/*! \param a Vector
+
+    \returns A normal vector in the direction of *a*.
+*/
+template<class Real> DEVICE inline vec2<Real> normalize(const vec2<Real>& a)
+    {
+    Real inverse_norm = fast::rsqrt(dot(a, a));
+    return a * inverse_norm;
+    }
+
 template<class Real> struct rotmat3;
 
 /////////////////////////////// quat ///////////////////////////////////