diff --git a/hoomd/Integrator.cc b/hoomd/Integrator.cc index 36b31fd3c3..479253c89b 100644 --- a/hoomd/Integrator.cc +++ b/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, diff --git a/hoomd/ParticleData.cc b/hoomd/ParticleData.cc index bf92bc7100..0643f0054c 100644 --- a/hoomd/ParticleData.cc +++ b/hoomd/ParticleData.cc @@ -2432,7 +2432,7 @@ unsigned int ParticleData::addParticle(unsigned int type) { // update reverse-lookup table ArrayHandle 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 diff --git a/hoomd/VectorMath.h b/hoomd/VectorMath.h index 3d7246359d..00a50f7e81 100644 --- a/hoomd/VectorMath.h +++ b/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 struct vec3 { @@ -320,6 +316,17 @@ template DEVICE inline vec3 cross(const vec3& a, const v return vec3(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 DEVICE inline vec3 normalize(const vec3& 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& a) { @@ -334,7 +341,7 @@ DEVICE inline Scalar4 vec_to_scalar4(const vec3& 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. @@ -342,10 +349,6 @@ DEVICE inline Scalar4 vec_to_scalar4(const vec3& a, Scalar w) 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 struct vec2 { @@ -610,6 +613,17 @@ template DEVICE inline Real perpdot(const vec2& a, const vec2< return dot(perp(a), b); } +/// Normalize the vector +/*! \param a Vector + + \returns A normal vector in the direction of *a*. +*/ +template DEVICE inline vec2 normalize(const vec2& a) + { + Real inverse_norm = fast::rsqrt(dot(a, a)); + return a * inverse_norm; + } + template struct rotmat3; /////////////////////////////// quat /////////////////////////////////// diff --git a/hoomd/md/ActiveForceCompute.cc b/hoomd/md/ActiveForceCompute.cc index 9a2510ac77..7d8c3ad618 100644 --- a/hoomd/md/ActiveForceCompute.cc +++ b/hoomd/md/ActiveForceCompute.cc @@ -9,27 +9,18 @@ #include -using namespace std; using namespace hoomd; -namespace py = pybind11; /*! \file ActiveForceCompute.cc \brief Contains code for the ActiveForceCompute class */ /*! \param rotation_diff rotational diffusion constant for all particles. - \param constraint specifies a constraint surface, to which particles are confined, - such as update.constraint_ellipsoid. Have to replace it with manifolds when implemented -*/ + */ ActiveForceCompute::ActiveForceCompute(std::shared_ptr sysdef, std::shared_ptr group, - Scalar rotation_diff, - Scalar3 P, - Scalar rx, - Scalar ry, - Scalar rz) - : ForceCompute(sysdef), m_group(group), m_rotationDiff(rotation_diff), m_P(P), m_rx(rx), - m_ry(ry), m_rz(rz) + Scalar rotation_diff) + : ForceCompute(sysdef), m_group(group), m_rotationDiff(rotation_diff) { // allocate memory for the per-type active_force storage and initialize them to (1.0,0,0) GlobalVector tmp_f_activeVec(m_pdata->getNTypes(), m_exec_conf); @@ -73,7 +64,7 @@ ActiveForceCompute::ActiveForceCompute(std::shared_ptr sysdef, ActiveForceCompute::~ActiveForceCompute() { - m_exec_conf->msg->notice(5) << "Destroying ActiveForceCompute" << endl; + m_exec_conf->msg->notice(5) << "Destroying ActiveForceCompute" << std::endl; } void ActiveForceCompute::setActiveForce(const std::string& type_name, pybind11::tuple v) @@ -82,13 +73,13 @@ void ActiveForceCompute::setActiveForce(const std::string& type_name, pybind11:: if (pybind11::len(v) != 3) { - throw invalid_argument("gamma_r values must be 3-tuples"); + throw std::invalid_argument("gamma_r values must be 3-tuples"); } // check for user errors if (typ >= m_pdata->getNTypes()) { - throw invalid_argument("Type does not exist"); + throw std::invalid_argument("Type does not exist"); } Scalar4 f_activeVec; @@ -108,7 +99,7 @@ void ActiveForceCompute::setActiveForce(const std::string& type_name, pybind11:: } else { - f_activeVec.x = 0; + f_activeVec.x = 1; f_activeVec.y = 0; f_activeVec.z = 0; f_activeVec.w = 0; @@ -140,13 +131,13 @@ void ActiveForceCompute::setActiveTorque(const std::string& type_name, pybind11: if (pybind11::len(v) != 3) { - throw invalid_argument("gamma_r values must be 3-tuples"); + throw std::invalid_argument("gamma_r values must be 3-tuples"); } // check for user errors if (typ >= m_pdata->getNTypes()) { - throw invalid_argument("Type does not exist"); + throw std::invalid_argument("Type does not exist"); } Scalar4 t_activeVec; @@ -259,34 +250,29 @@ void ActiveForceCompute::rotationalDiffusion(uint64_t timestep) { unsigned int idx = m_group->getMemberIndex(i); unsigned int type = __scalar_as_int(h_pos.data[idx].w); - unsigned int ptag = h_tag.data[idx]; - hoomd::RandomGenerator rng( - hoomd::Seed(hoomd::RNGIdentifier::ActiveForceCompute, timestep, m_sysdef->getSeed()), - hoomd::Counter(ptag)); - - quat quati(h_orientation.data[idx]); - if (m_sysdef->getNDimensions() == 2) // 2D + if (h_f_actVec.data[type].w != 0) { - Scalar delta_theta; // rotational diffusion angle - delta_theta = hoomd::NormalDistribution(m_rotationConst)(rng); - Scalar theta - = delta_theta - / 2.0; // half angle to calculate the quaternion which represents the rotation - vec3 b(0, 0, slow::sin(theta)); - - quat rot_quat(slow::cos(theta), b); // rotational diffusion quaternion - - quati = rot_quat * quati; // rotational diffusion quaternion applied to orientation - h_orientation.data[idx].x = quati.s; - h_orientation.data[idx].y = quati.v.x; - h_orientation.data[idx].z = quati.v.y; - h_orientation.data[idx].w = quati.v.z; - // In 2D, the only meaningful torque vector is out of plane and should not change - } - else // 3D: Following Stenhammar, Soft Matter, 2014 - { - if (m_rx == 0) // if no constraint + unsigned int ptag = h_tag.data[idx]; + hoomd::RandomGenerator rng(hoomd::Seed(hoomd::RNGIdentifier::ActiveForceCompute, + timestep, + m_sysdef->getSeed()), + hoomd::Counter(ptag)); + + quat quati(h_orientation.data[idx]); + + if (m_sysdef->getNDimensions() == 2) // 2D + { + Scalar delta_theta = hoomd::NormalDistribution(m_rotationConst)(rng); + + vec3 b(0, 0, 1.0); + quat rot_quat = quat::fromAxisAngle(b, delta_theta); + + quati = rot_quat * quati; // rotational diffusion quaternion applied to orientation + h_orientation.data[idx] = quat_to_scalar4(quati); + // In 2D, the only meaningful torque vector is out of plane and should not change + } + else // 3D: Following Stenhammar, Soft Matter, 2014 { hoomd::SpherePointGenerator unit_vec; vec3 rand_vec; @@ -298,46 +284,12 @@ void ActiveForceCompute::rotationalDiffusion(uint64_t timestep) vec3 fi = rotate(quati, f); // rotate active force vector from local to global frame - vec3 aux_vec; // rotation axis - aux_vec.x = fi.y * rand_vec.z - fi.z * rand_vec.y; - aux_vec.y = fi.z * rand_vec.x - fi.x * rand_vec.z; - aux_vec.z = fi.x * rand_vec.y - fi.y * rand_vec.x; - Scalar aux_vec_mag = 1.0 - / slow::sqrt(aux_vec.x * aux_vec.x + aux_vec.y * aux_vec.y - + aux_vec.z * aux_vec.z); - aux_vec.x *= aux_vec_mag; - aux_vec.y *= aux_vec_mag; - aux_vec.z *= aux_vec_mag; + vec3 aux_vec = cross(fi, rand_vec); // rotation axis + Scalar aux_vec_mag = slow::rsqrt(dot(aux_vec, aux_vec)); + aux_vec *= aux_vec_mag; Scalar delta_theta = hoomd::NormalDistribution(m_rotationConst)(rng); - Scalar theta - = delta_theta - / 2.0; // half angle to calculate the quaternion which represents the rotation - quat rot_quat(slow::cos(theta), - slow::sin(theta) - * aux_vec); // rotational diffusion quaternion - - quati = rot_quat * quati; // rotational diffusion quaternion applied to orientation - h_orientation.data[idx] = quat_to_scalar4(quati); - } - else // if constraint exists - { - EvaluatorConstraintEllipsoid Ellipsoid(m_P, m_rx, m_ry, m_rz); - - Scalar3 current_pos - = make_scalar3(h_pos.data[idx].x, h_pos.data[idx].y, h_pos.data[idx].z); - Scalar3 norm_scalar3 = Ellipsoid.evalNormal( - current_pos); // the normal vector to which the particles are confined. - - vec3 norm; - norm = vec3(norm_scalar3); - - Scalar delta_theta = hoomd::NormalDistribution(m_rotationConst)(rng); - Scalar theta - = delta_theta - / 2.0; // half angle to calculate the quaternion which represents the rotation - quat rot_quat(slow::cos(theta), - slow::sin(theta) * norm); // rotational diffusion quaternion + quat rot_quat = quat::fromAxisAngle(aux_vec, delta_theta); quati = rot_quat * quati; // rotational diffusion quaternion applied to orientation h_orientation.data[idx] = quat_to_scalar4(quati); @@ -346,74 +298,8 @@ void ActiveForceCompute::rotationalDiffusion(uint64_t timestep) } } -/*! This function sets an ellipsoid surface constraint for all active particles. Torque is not - * considered here - */ -void ActiveForceCompute::setConstraint() - { - EvaluatorConstraintEllipsoid Ellipsoid(m_P, m_rx, m_ry, m_rz); - - // array handles - ArrayHandle h_f_actVec(m_f_activeVec, access_location::host, access_mode::read); - ArrayHandle h_orientation(m_pdata->getOrientationArray(), - access_location::host, - access_mode::readwrite); - ArrayHandle h_pos(m_pdata->getPositions(), access_location::host, access_mode::read); - - assert(h_f_actVec.data != NULL); - assert(h_pos.data != NULL); - assert(h_orientation.data != NULL); - - for (unsigned int i = 0; i < m_group->getNumMembers(); i++) - { - unsigned int idx = m_group->getMemberIndex(i); - unsigned int type = __scalar_as_int(h_pos.data[idx].w); - - Scalar3 current_pos = make_scalar3(h_pos.data[idx].x, h_pos.data[idx].y, h_pos.data[idx].z); - - Scalar3 norm_scalar3 = Ellipsoid.evalNormal( - current_pos); // the normal vector to which the particles are confined. - vec3 norm = vec3(norm_scalar3); - - Scalar3 f = make_scalar3(h_f_actVec.data[type].x, - h_f_actVec.data[type].y, - h_f_actVec.data[type].z); - quat quati(h_orientation.data[idx]); - vec3 fi - = rotate(quati, - vec3(f)); // rotate active force vector from local to global frame - - Scalar dot_prod = fi.x * norm.x + fi.y * norm.y + fi.z * norm.z; - - Scalar dot_perp_prod = slow::sqrt(1 - dot_prod * dot_prod); - - Scalar phi_half = slow::atan(dot_prod / dot_perp_prod) / 2.0; - - fi.x -= norm.x * dot_prod; - fi.y -= norm.y * dot_prod; - fi.z -= norm.z * dot_prod; - - Scalar new_norm = 1.0 / slow::sqrt(fi.x * fi.x + fi.y * fi.y + fi.z * fi.z); - - fi.x *= new_norm; - fi.y *= new_norm; - fi.z *= new_norm; - - vec3 rot_vec = cross(norm, fi); - rot_vec.x *= slow::sin(phi_half); - rot_vec.y *= slow::sin(phi_half); - rot_vec.z *= slow::sin(phi_half); - - quat rot_quat(cos(phi_half), rot_vec); - - quati = rot_quat * quati; - - h_orientation.data[idx] = quat_to_scalar4(quati); - } - } - -/*! This function applies constraints, rotational diffusion, and sets forces for all active - particles \param timestep Current timestep +/*! This function applies rotational diffusion and sets forces for all active particles + \param timestep Current timestep */ void ActiveForceCompute::computeForces(uint64_t timestep) { @@ -426,10 +312,6 @@ void ActiveForceCompute::computeForces(uint64_t timestep) last_computed = timestep; - if (m_rx != 0) - { - setConstraint(); // apply surface constraints to active particles active force vectors - } if (m_rotationDiff != 0) { rotationalDiffusion(timestep); // apply rotational diffusion to active particles @@ -446,18 +328,13 @@ void ActiveForceCompute::computeForces(uint64_t timestep) m_prof->pop(m_exec_conf); } -void export_ActiveForceCompute(py::module& m) +void export_ActiveForceCompute(pybind11::module& m) { - py::class_>( + pybind11::class_>( m, "ActiveForceCompute") - .def(py::init, - std::shared_ptr, - Scalar, - Scalar3, - Scalar, - Scalar, - Scalar>()) + .def(pybind11:: + init, std::shared_ptr, Scalar>()) .def_property("rotation_diff", &ActiveForceCompute::getRdiff, &ActiveForceCompute::setRdiff) .def("setActiveForce", &ActiveForceCompute::setActiveForce) .def("getActiveForce", &ActiveForceCompute::getActiveForce) diff --git a/hoomd/md/ActiveForceCompute.h b/hoomd/md/ActiveForceCompute.h index fcac5fc2af..eb70b53b7f 100644 --- a/hoomd/md/ActiveForceCompute.h +++ b/hoomd/md/ActiveForceCompute.h @@ -9,8 +9,6 @@ #include "hoomd/VectorMath.h" #include -#include "EvaluatorConstraintEllipsoid.h" - /*! \file ActiveForceCompute.h \brief Declares a class for computing active forces and torques */ @@ -33,11 +31,7 @@ class PYBIND11_EXPORT ActiveForceCompute : public ForceCompute //! Constructs the compute ActiveForceCompute(std::shared_ptr sysdef, std::shared_ptr group, - Scalar rotation_diff, - Scalar3 P, - Scalar rx, - Scalar ry, - Scalar rz); + Scalar rotation_diff); //! Destructor ~ActiveForceCompute(); @@ -84,16 +78,9 @@ class PYBIND11_EXPORT ActiveForceCompute : public ForceCompute //! Orientational diffusion for spherical particles virtual void rotationalDiffusion(uint64_t timestep); - //! Set constraints if particles confined to a surface - virtual void setConstraint(); - std::shared_ptr m_group; //!< Group of particles on which this force is applied Scalar m_rotationDiff; Scalar m_rotationConst; - Scalar3 m_P; //!< Position of the Ellipsoid - Scalar m_rx; //!< Radius in X direction of the Ellipsoid - Scalar m_ry; //!< Radius in Y direction of the Ellipsoid - Scalar m_rz; //!< Radius in Z direction of the Ellipsoid GlobalVector m_f_activeVec; //! active force unit vectors and magnitudes for each particle type diff --git a/hoomd/md/ActiveForceComputeGPU.cc b/hoomd/md/ActiveForceComputeGPU.cc index a50a22e720..f1c79af762 100644 --- a/hoomd/md/ActiveForceComputeGPU.cc +++ b/hoomd/md/ActiveForceComputeGPU.cc @@ -20,17 +20,12 @@ using namespace std; non-point-like anisotropic particles. /param orientation_reverse_link When True, the particle's orientation is set to match the active force vector. Useful for for using a particle's orientation to log the active force vector. Not recommended for anisotropic particles \param - rotation_diff rotational diffusion constant for all particles. \param constraint specifies a - constraint surface, to which particles are confined, such as update.constraint_ellipsoid. + rotation_diff rotational diffusion constant for all particles. */ ActiveForceComputeGPU::ActiveForceComputeGPU(std::shared_ptr sysdef, std::shared_ptr group, - Scalar rotation_diff, - Scalar3 P, - Scalar rx, - Scalar ry, - Scalar rz) - : ActiveForceCompute(sysdef, group, rotation_diff, P, rx, ry, rz), m_block_size(256) + Scalar rotation_diff) + : ActiveForceCompute(sysdef, group, rotation_diff) { if (!m_exec_conf->isCUDAEnabled()) { @@ -40,6 +35,16 @@ ActiveForceComputeGPU::ActiveForceComputeGPU(std::shared_ptr s throw std::runtime_error("Error initializing ActiveForceComputeGPU"); } + // initialize autotuner + std::vector valid_params; + unsigned int warp_size = m_exec_conf->dev_prop.warpSize; + for (unsigned int block_size = warp_size; block_size <= 1024; block_size += warp_size) + valid_params.push_back(block_size); + + m_tuner_force.reset(new Autotuner(valid_params, 5, 100000, "active_force", this->m_exec_conf)); + m_tuner_diffusion.reset( + new Autotuner(valid_params, 5, 100000, "active_diffusion", this->m_exec_conf)); + // unsigned int N = m_pdata->getNGlobal(); // unsigned int group_size = m_group->getNumMembersGlobal(); unsigned int type = m_pdata->getNTypes(); @@ -97,6 +102,9 @@ void ActiveForceComputeGPU::setForces() unsigned int group_size = m_group->getNumMembers(); unsigned int N = m_pdata->getN(); + // compute the forces on the GPU + m_tuner_force->begin(); + gpu_compute_active_force_set_forces(group_size, d_index_array.data, d_force.data, @@ -105,12 +113,13 @@ void ActiveForceComputeGPU::setForces() d_orientation.data, d_f_actVec.data, d_t_actVec.data, - m_P, - m_rx, - m_ry, - m_rz, N, - m_block_size); + m_tuner_force->getParam()); + + if (m_exec_conf->isCUDAErrorCheckingEnabled()) + CHECK_CUDA_ERROR(); + + m_tuner_force->end(); } /*! This function applies rotational diffusion to all active particles. The angle between the torque @@ -136,53 +145,25 @@ void ActiveForceComputeGPU::rotationalDiffusion(uint64_t timestep) bool is2D = (m_sysdef->getNDimensions() == 2); unsigned int group_size = m_group->getNumMembers(); + // perform the update on the GPU + m_tuner_diffusion->begin(); + gpu_compute_active_force_rotational_diffusion(group_size, d_tag.data, d_index_array.data, d_pos.data, d_orientation.data, d_f_actVec.data, - m_P, - m_rx, - m_ry, - m_rz, is2D, m_rotationConst, timestep, m_sysdef->getSeed(), - m_block_size); - } - -/*! This function sets an ellipsoid surface constraint for all active particles - */ -void ActiveForceComputeGPU::setConstraint() - { - EvaluatorConstraintEllipsoid Ellipsoid(m_P, m_rx, m_ry, m_rz); + m_tuner_diffusion->getParam()); - // array handles - ArrayHandle d_f_actVec(m_f_activeVec, access_location::device, access_mode::readwrite); - ArrayHandle d_pos(m_pdata->getPositions(), access_location::device, access_mode::read); - ArrayHandle d_orientation(m_pdata->getOrientationArray(), - access_location::device, - access_mode::readwrite); - ArrayHandle d_index_array(m_group->getIndexArray(), - access_location::device, - access_mode::read); - - assert(d_pos.data != NULL); - - unsigned int group_size = m_group->getNumMembers(); + if (m_exec_conf->isCUDAErrorCheckingEnabled()) + CHECK_CUDA_ERROR(); - gpu_compute_active_force_set_constraints(group_size, - d_index_array.data, - d_pos.data, - d_orientation.data, - d_f_actVec.data, - m_P, - m_rx, - m_ry, - m_rz, - m_block_size); + m_tuner_diffusion->end(); } void export_ActiveForceComputeGPU(py::module& m) @@ -190,11 +171,5 @@ void export_ActiveForceComputeGPU(py::module& m) py::class_>( m, "ActiveForceComputeGPU") - .def(py::init, - std::shared_ptr, - Scalar, - Scalar3, - Scalar, - Scalar, - Scalar>()); + .def(py::init, std::shared_ptr, Scalar>()); } diff --git a/hoomd/md/ActiveForceComputeGPU.cu b/hoomd/md/ActiveForceComputeGPU.cu index 9c3824d0ff..cf6f68e2d4 100644 --- a/hoomd/md/ActiveForceComputeGPU.cu +++ b/hoomd/md/ActiveForceComputeGPU.cu @@ -4,7 +4,6 @@ // Maintainer: joaander #include "ActiveForceComputeGPU.cuh" -#include "EvaluatorConstraintEllipsoid.h" #include "hoomd/RNGIdentifiers.h" #include "hoomd/RandomNumbers.h" #include "hoomd/TextureTools.h" @@ -25,10 +24,6 @@ using namespace hoomd; \param d_orientation particle orientation on device \param d_f_act particle active force unit vector \param d_t_act particle active torque unit vector - \param P position of the ellipsoid constraint - \param rx radius of the ellipsoid in x direction - \param ry radius of the ellipsoid in y direction - \param rz radius of the ellipsoid in z direction \param orientationLink check if particle orientation is linked to active force vector */ __global__ void gpu_compute_active_force_set_forces_kernel(const unsigned int group_size, @@ -39,10 +34,6 @@ __global__ void gpu_compute_active_force_set_forces_kernel(const unsigned int gr const Scalar4* d_orientation, const Scalar4* d_f_act, const Scalar4* d_t_act, - const Scalar3 P, - const Scalar rx, - const Scalar ry, - const Scalar rz, const unsigned int N) { unsigned int group_idx = blockIdx.x * blockDim.x + threadIdx.x; @@ -67,85 +58,11 @@ __global__ void gpu_compute_active_force_set_forces_kernel(const unsigned int gr d_torque[idx] = vec_to_scalar4(ti, 0); } -//! Kernel for adjusting active force vectors to align parallel to an ellipsoid surface constraint -//! on the GPU -/*! \param group_size number of particles - \param d_index_array stores list to convert group index to global tag - \param d_pos particle positions on device - \param d_f_act particle active force unit vector - \param P position of the ellipsoid constraint - \param rx radius of the ellipsoid in x direction - \param ry radius of the ellipsoid in y direction - \param rz radius of the ellipsoid in z direction -*/ -__global__ void gpu_compute_active_force_set_constraints_kernel(const unsigned int group_size, - unsigned int* d_index_array, - const Scalar4* d_pos, - Scalar4* d_orientation, - const Scalar4* d_f_act, - const Scalar3 P, - const Scalar rx, - const Scalar ry, - const Scalar rz) - { - unsigned int group_idx = blockIdx.x * blockDim.x + threadIdx.x; - if (group_idx >= group_size) - return; - - unsigned int idx = d_index_array[group_idx]; - Scalar4 posidx = __ldg(d_pos + idx); - unsigned int type = __scalar_as_int(posidx.w); - - EvaluatorConstraintEllipsoid Ellipsoid(P, rx, ry, rz); - Scalar3 current_pos = make_scalar3(posidx.x, posidx.y, posidx.z); - - Scalar3 norm_scalar3 = Ellipsoid.evalNormal( - current_pos); // the normal vector to which the particles are confined. - vec3 norm = vec3(norm_scalar3); - - Scalar4 fact = __ldg(d_f_act + type); - - vec3 f(fact.x, fact.y, fact.z); - quat quati(__ldg(d_orientation + idx)); - vec3 fi = rotate(quati, f); - - Scalar dot_prod = fi.x * norm.x + fi.y * norm.y + fi.z * norm.z; - - Scalar dot_perp_prod = slow::sqrt(1 - dot_prod * dot_prod); - - Scalar phi_half = slow::atan(dot_prod / dot_perp_prod) / 2.0; - - fi.x -= norm.x * dot_prod; - fi.y -= norm.y * dot_prod; - fi.z -= norm.z * dot_prod; - - Scalar new_norm = 1.0 / slow::sqrt(fi.x * fi.x + fi.y * fi.y + fi.z * fi.z); - - fi.x *= new_norm; - fi.y *= new_norm; - fi.z *= new_norm; - - vec3 rot_vec = cross(norm, fi); - rot_vec.x *= slow::sin(phi_half); - rot_vec.y *= slow::sin(phi_half); - rot_vec.z *= slow::sin(phi_half); - - quat rot_quat(cos(phi_half), rot_vec); - - quati = rot_quat * quati; - - d_orientation[idx] = quat_to_scalar4(quati); - } - //! Kernel for applying rotational diffusion to active force vectors on the GPU /*! \param group_size number of particles \param d_index_array stores list to convert group index to global tag \param d_pos particle positions on device \param d_f_act particle active force unit vector - \param P position of the ellipsoid constraint - \param rx radius of the ellipsoid in x direction - \param ry radius of the ellipsoid in y direction - \param rz radius of the ellipsoid in z direction \param is2D check if simulation is 2D or 3D \param rotationConst particle rotational diffusion constant \param seed seed for random number generator @@ -156,10 +73,6 @@ __global__ void gpu_compute_active_force_rotational_diffusion_kernel(const unsig const Scalar4* d_pos, Scalar4* d_orientation, const Scalar4* d_f_act, - const Scalar3 P, - const Scalar rx, - const Scalar ry, - const Scalar rz, bool is2D, const Scalar rotationConst, const uint64_t timestep, @@ -172,77 +85,45 @@ __global__ void gpu_compute_active_force_rotational_diffusion_kernel(const unsig unsigned int idx = d_index_array[group_idx]; Scalar4 posidx = __ldg(d_pos + idx); unsigned int type = __scalar_as_int(posidx.w); - unsigned int ptag = d_tag[group_idx]; - quat quati(__ldg(d_orientation + idx)); - - hoomd::RandomGenerator rng( - hoomd::Seed(hoomd::RNGIdentifier::ActiveForceCompute, timestep, seed), - hoomd::Counter(ptag)); + Scalar4 fact = __ldg(d_f_act + type); - if (is2D) // 2D + if (fact.w != 0) { - Scalar delta_theta; // rotational diffusion angle - delta_theta = hoomd::NormalDistribution(rotationConst)(rng); - Scalar theta - = delta_theta / 2.0; // angle on plane defining orientation of active force vector - vec3 b(0, 0, slow::sin(theta)); + unsigned int ptag = d_tag[group_idx]; - quat rot_quat(slow::cos(theta), b); + quat quati(__ldg(d_orientation + idx)); - quati = rot_quat * quati; - d_orientation[idx] = quat_to_scalar4(quati); - // in 2D there is only one meaningful direction for torque - } - else // 3D: Following Stenhammar, Soft Matter, 2014 - { - if (rx == 0) // if no constraint + hoomd::RandomGenerator rng( + hoomd::Seed(hoomd::RNGIdentifier::ActiveForceCompute, timestep, seed), + hoomd::Counter(ptag)); + + if (is2D) // 2D + { + Scalar delta_theta = hoomd::NormalDistribution(rotationConst)(rng); + + vec3 b(0, 0, 1.0); + quat rot_quat = quat::fromAxisAngle(b, delta_theta); + + quati = rot_quat * quati; + d_orientation[idx] = quat_to_scalar4(quati); + // in 2D there is only one meaningful direction for torque + } + else // 3D: Following Stenhammar, Soft Matter, 2014 { hoomd::SpherePointGenerator unit_vec; vec3 rand_vec; unit_vec(rng, rand_vec); - Scalar4 fact = __ldg(d_f_act + type); - vec3 f(fact.x, fact.y, fact.z); vec3 fi = rotate(quati, f); - vec3 aux_vec; - aux_vec.x = fi.y * rand_vec.z - fi.z * rand_vec.y; - aux_vec.y = fi.z * rand_vec.x - fi.x * rand_vec.z; - aux_vec.z = fi.x * rand_vec.y - fi.y * rand_vec.x; - Scalar aux_vec_mag = 1.0 - / slow::sqrt(aux_vec.x * aux_vec.x + aux_vec.y * aux_vec.y - + aux_vec.z * aux_vec.z); - aux_vec.x *= aux_vec_mag; - aux_vec.y *= aux_vec_mag; - aux_vec.z *= aux_vec_mag; + vec3 aux_vec = cross(fi, rand_vec); // rotation axis + Scalar aux_vec_mag = slow::rsqrt(dot(aux_vec, aux_vec)); + aux_vec *= aux_vec_mag; Scalar delta_theta = hoomd::NormalDistribution(rotationConst)(rng); - Scalar theta - = delta_theta / 2.0; // angle on plane defining orientation of active force vector - quat rot_quat(slow::cos(theta), slow::sin(theta) * aux_vec); - - quati = rot_quat * quati; - d_orientation[idx].x = quati.s; - d_orientation[idx].y = quati.v.x; - d_orientation[idx].z = quati.v.y; - d_orientation[idx].w = quati.v.z; - } - else // if constraint - { - EvaluatorConstraintEllipsoid Ellipsoid(P, rx, ry, rz); - Scalar3 current_pos = make_scalar3(posidx.x, posidx.y, posidx.z); - - Scalar3 norm_scalar3 = Ellipsoid.evalNormal( - current_pos); // the normal vector to which the particles are confined. - vec3 norm; - norm = vec3(norm_scalar3); - - Scalar delta_theta = hoomd::NormalDistribution(rotationConst)(rng); - Scalar theta - = delta_theta / 2.0; // angle on plane defining orientation of active force vector - quat rot_quat(slow::cos(theta), slow::sin(theta) * norm); + quat rot_quat = quat::fromAxisAngle(aux_vec, delta_theta); quati = rot_quat * quati; d_orientation[idx] = quat_to_scalar4(quati); @@ -258,10 +139,6 @@ hipError_t gpu_compute_active_force_set_forces(const unsigned int group_size, const Scalar4* d_orientation, const Scalar4* d_f_act, const Scalar4* d_t_act, - const Scalar3& P, - const Scalar rx, - const Scalar ry, - const Scalar rz, const unsigned int N, unsigned int block_size) { @@ -284,57 +161,16 @@ hipError_t gpu_compute_active_force_set_forces(const unsigned int group_size, d_orientation, d_f_act, d_t_act, - P, - rx, - ry, - rz, N); return hipSuccess; } -hipError_t gpu_compute_active_force_set_constraints(const unsigned int group_size, - unsigned int* d_index_array, - const Scalar4* d_pos, - Scalar4* d_orientation, - const Scalar4* d_f_act, - const Scalar3& P, - const Scalar rx, - const Scalar ry, - const Scalar rz, - unsigned int block_size) - { - // setup the grid to run the kernel - dim3 grid(group_size / block_size + 1, 1, 1); - dim3 threads(block_size, 1, 1); - - // run the kernel - hipLaunchKernelGGL((gpu_compute_active_force_set_constraints_kernel), - dim3(grid), - dim3(threads), - 0, - 0, - group_size, - d_index_array, - d_pos, - d_orientation, - d_f_act, - P, - rx, - ry, - rz); - return hipSuccess; - } - hipError_t gpu_compute_active_force_rotational_diffusion(const unsigned int group_size, unsigned int* d_tag, unsigned int* d_index_array, const Scalar4* d_pos, Scalar4* d_orientation, const Scalar4* d_f_act, - const Scalar3& P, - const Scalar rx, - const Scalar ry, - const Scalar rz, bool is2D, const Scalar rotationConst, const uint64_t timestep, @@ -357,10 +193,6 @@ hipError_t gpu_compute_active_force_rotational_diffusion(const unsigned int grou d_pos, d_orientation, d_f_act, - P, - rx, - ry, - rz, is2D, rotationConst, timestep, diff --git a/hoomd/md/ActiveForceComputeGPU.cuh b/hoomd/md/ActiveForceComputeGPU.cuh index bd8a11547e..e429014374 100644 --- a/hoomd/md/ActiveForceComputeGPU.cuh +++ b/hoomd/md/ActiveForceComputeGPU.cuh @@ -23,34 +23,15 @@ hipError_t gpu_compute_active_force_set_forces(const unsigned int group_size, const Scalar4* d_orientation, const Scalar4* d_f_act, const Scalar4* d_t_act, - const Scalar3& P, - const Scalar rx, - const Scalar ry, - const Scalar rz, const unsigned int N, unsigned int block_size); -hipError_t gpu_compute_active_force_set_constraints(const unsigned int group_size, - unsigned int* d_index_array, - const Scalar4* d_pos, - Scalar4* d_orientation, - const Scalar4* d_f_act, - const Scalar3& P, - const Scalar rx, - const Scalar ry, - const Scalar rz, - unsigned int block_size); - hipError_t gpu_compute_active_force_rotational_diffusion(const unsigned int group_size, unsigned int* d_tag, unsigned int* d_index_array, const Scalar4* d_pos, Scalar4* d_orientation, const Scalar4* d_f_act, - const Scalar3& P, - const Scalar rx, - const Scalar ry, - const Scalar rz, bool is2D, const Scalar rotationDiff, const uint64_t timestep, diff --git a/hoomd/md/ActiveForceComputeGPU.h b/hoomd/md/ActiveForceComputeGPU.h index 69ebb26cb9..3e8b973e44 100644 --- a/hoomd/md/ActiveForceComputeGPU.h +++ b/hoomd/md/ActiveForceComputeGPU.h @@ -4,6 +4,7 @@ // Maintainer: joaander #include "ActiveForceCompute.h" +#include "hoomd/Autotuner.h" /*! \file ActiveForceComputeGPU.h \brief Declares a class for computing active forces on the GPU @@ -27,23 +28,30 @@ class PYBIND11_EXPORT ActiveForceComputeGPU : public ActiveForceCompute //! Constructs the compute ActiveForceComputeGPU(std::shared_ptr sysdef, std::shared_ptr group, - Scalar rotation_diff, - Scalar3 P, - Scalar rx, - Scalar ry, - Scalar rz); + Scalar rotation_diff); + + //! Set autotuner parameters + /*! \param enable Enable/disable autotuning + \param period period (approximate) in time steps when returning occurs + */ + virtual void setAutotunerParams(bool enable, unsigned int period) + { + ActiveForceCompute::setAutotunerParams(enable, period); + m_tuner_force->setPeriod(period); + m_tuner_force->setEnabled(enable); + m_tuner_diffusion->setPeriod(period); + m_tuner_diffusion->setEnabled(enable); + } protected: - unsigned int m_block_size; //!< block size to execute on the GPU + std::unique_ptr m_tuner_force; //!< Autotuner for block size (force kernel) + std::unique_ptr m_tuner_diffusion; //!< Autotuner for block size (diff kernel) //! Set forces for particles virtual void setForces(); //! Orientational diffusion for spherical particles virtual void rotationalDiffusion(uint64_t timestep); - - //! Set constraints if particles confined to a surface - virtual void setConstraint(); }; //! Exports the ActiveForceComputeGPU Class to python diff --git a/hoomd/md/ActiveForceConstraintCompute.h b/hoomd/md/ActiveForceConstraintCompute.h new file mode 100644 index 0000000000..c937300060 --- /dev/null +++ b/hoomd/md/ActiveForceConstraintCompute.h @@ -0,0 +1,238 @@ +// Copyright (c) 2009-2021 The Regents of the University of Michigan +// This file is part of the HOOMD-blue project, released under the BSD 3-Clause License. + +// Maintainer: joaander + +#include "ActiveForceCompute.h" +#include "hoomd/RNGIdentifiers.h" +#include "hoomd/RandomNumbers.h" + +/*! \file ActiveForceConstraintCompute.h + \brief Declares a class for computing active forces and torques +*/ + +#ifdef __HIPCC__ +#error This header cannot be compiled by nvcc +#endif + +#include + +#ifndef __ACTIVEFORCECONSTRAINTCOMPUTE_H__ +#define __ACTIVEFORCECONSTRAINTCOMPUTE_H__ + +//! Adds an active force to a number of particles +/*! \ingroup computes + */ +template +class PYBIND11_EXPORT ActiveForceConstraintCompute : public ActiveForceCompute + { + public: + //! Constructs the compute + ActiveForceConstraintCompute(std::shared_ptr sysdef, + std::shared_ptr group, + Scalar rotation_diff, + Manifold manifold); + // + //! Destructor + ~ActiveForceConstraintCompute(); + + protected: + //! Actually compute the forces + virtual void computeForces(uint64_t timestep); + + //! Orientational diffusion for spherical particles + virtual void rotationalDiffusion(uint64_t timestep); + + //! Set constraints if particles confined to a surface + virtual void setConstraint(); + + //! Helper function to be called when box changes + void setBoxChange() + { + m_box_changed = true; + } + + Manifold m_manifold; //!< Constraining Manifold + bool m_box_changed; + }; + +/*! \param sysdef The system definition + \param group Particle group + \param rotation_diff Rotational diffusion coefficient + \param manifold Manifold constraint + */ +template +ActiveForceConstraintCompute::ActiveForceConstraintCompute( + std::shared_ptr sysdef, + std::shared_ptr group, + Scalar rotation_diff, + Manifold manifold) + : ActiveForceCompute(sysdef, group, rotation_diff), m_manifold(manifold), m_box_changed(true) + { + m_pdata->getBoxChangeSignal() + .template connect, + &ActiveForceConstraintCompute::setBoxChange>(this); + } + +template ActiveForceConstraintCompute::~ActiveForceConstraintCompute() + { + m_pdata->getBoxChangeSignal() + .template disconnect, + &ActiveForceConstraintCompute::setBoxChange>(this); + m_exec_conf->msg->notice(5) << "Destroying ActiveForceConstraintCompute" << std::endl; + } + +/*! This function applies rotational diffusion to the orientations of all active particles. The + orientation of any torque vector + * relative to the force vector is preserved + \param timestep Current timestep +*/ +template +void ActiveForceConstraintCompute::rotationalDiffusion(uint64_t timestep) + { + // array handles + ArrayHandle h_pos(m_pdata->getPositions(), access_location::host, access_mode::read); + ArrayHandle h_orientation(m_pdata->getOrientationArray(), + access_location::host, + access_mode::readwrite); + ArrayHandle h_tag(m_pdata->getTags(), access_location::host, access_mode::read); + + assert(h_pos.data != NULL); + assert(h_orientation.data != NULL); + assert(h_tag.data != NULL); + + for (unsigned int i = 0; i < m_group->getNumMembers(); i++) + { + unsigned int idx = m_group->getMemberIndex(i); + unsigned int ptag = h_tag.data[idx]; + hoomd::RandomGenerator rng( + hoomd::Seed(hoomd::RNGIdentifier::ActiveForceCompute, timestep, m_sysdef->getSeed()), + hoomd::Counter(ptag)); + + quat quati(h_orientation.data[idx]); + + Scalar3 current_pos = make_scalar3(h_pos.data[idx].x, h_pos.data[idx].y, h_pos.data[idx].z); + + vec3 norm = normalize(vec3(m_manifold.derivative(current_pos))); + + Scalar delta_theta = hoomd::NormalDistribution(m_rotationConst)(rng); + + quat rot_quat = quat::fromAxisAngle(norm, delta_theta); + + quati = rot_quat * quati; // rotational diffusion quaternion applied to orientation + h_orientation.data[idx] = quat_to_scalar4(quati); + } + } + +/*! This function sets a manifold constraint for all active particles. Torque is not considered here + */ +template void ActiveForceConstraintCompute::setConstraint() + { + // array handles + ArrayHandle h_f_actVec(m_f_activeVec, access_location::host, access_mode::read); + ArrayHandle h_orientation(m_pdata->getOrientationArray(), + access_location::host, + access_mode::readwrite); + ArrayHandle h_pos(m_pdata->getPositions(), access_location::host, access_mode::read); + + assert(h_f_actVec.data != NULL); + assert(h_pos.data != NULL); + assert(h_orientation.data != NULL); + + for (unsigned int i = 0; i < m_group->getNumMembers(); i++) + { + unsigned int idx = m_group->getMemberIndex(i); + unsigned int type = __scalar_as_int(h_pos.data[idx].w); + + if (h_f_actVec.data[type].w != 0) + { + Scalar3 current_pos + = make_scalar3(h_pos.data[idx].x, h_pos.data[idx].y, h_pos.data[idx].z); + + vec3 norm = normalize(vec3(m_manifold.derivative(current_pos))); + + vec3 f(h_f_actVec.data[type].x, + h_f_actVec.data[type].y, + h_f_actVec.data[type].z); + quat quati(h_orientation.data[idx]); + vec3 fi = rotate(quati, + f); // rotate active force vector from local to global frame + + Scalar dot_prod = dot(fi, norm); + + Scalar dot_perp_prod = slow::rsqrt(1 - dot_prod * dot_prod); + + Scalar phi = slow::atan(dot_prod * dot_perp_prod); + + fi.x -= norm.x * dot_prod; + fi.y -= norm.y * dot_prod; + fi.z -= norm.z * dot_prod; + + Scalar new_norm = slow::rsqrt(dot(fi, fi)); + fi *= new_norm; + + vec3 rot_vec = cross(norm, fi); + quat rot_quat = quat::fromAxisAngle(rot_vec, phi); + + quati = rot_quat * quati; + + h_orientation.data[idx] = quat_to_scalar4(quati); + } + } + } + +/*! This function applies constraints, rotational diffusion, and sets forces for all active + particles \param timestep Current timestep +*/ +template +void ActiveForceConstraintCompute::computeForces(uint64_t timestep) + { + if (m_prof) + m_prof->push(m_exec_conf, "ActiveForceConstraintCompute"); + + if (last_computed != timestep) + { + m_rotationConst = slow::sqrt(2.0 * m_rotationDiff * m_deltaT); + + if (m_box_changed) + { + if (!m_manifold.fitsInsideBox(m_pdata->getGlobalBox())) + { + throw std::runtime_error("Parts of the manifold are outside the box"); + } + m_box_changed = false; + } + + last_computed = timestep; + + setConstraint(); // apply manifold constraints to active particles active force vectors + + if (m_rotationDiff != 0) + { + rotationalDiffusion(timestep); // apply rotational diffusion to active particles + } + setForces(); // set forces for particles + } + +#ifdef ENABLE_HIP + if (m_exec_conf->isCUDAErrorCheckingEnabled()) + CHECK_CUDA_ERROR(); +#endif + + if (m_prof) + m_prof->pop(m_exec_conf); + } + +template +void export_ActiveForceConstraintCompute(pybind11::module& m, const std::string& name) + { + pybind11::class_, + ActiveForceCompute, + std::shared_ptr>>(m, name.c_str()) + .def(pybind11::init, + std::shared_ptr, + Scalar, + Manifold>()); + } + +#endif diff --git a/hoomd/md/ActiveForceConstraintComputeGPU.cuh b/hoomd/md/ActiveForceConstraintComputeGPU.cuh new file mode 100644 index 0000000000..9ac9253724 --- /dev/null +++ b/hoomd/md/ActiveForceConstraintComputeGPU.cuh @@ -0,0 +1,218 @@ +// Copyright (c) 2009-2021 The Regents of the University of Michigan +// This file is part of the HOOMD-blue project, released under the BSD 3-Clause License. + +// Maintainer: joaander + +#include "hip/hip_runtime.h" +#include "hoomd/HOOMDMath.h" +#include "hoomd/ParticleData.cuh" +#include "hoomd/RNGIdentifiers.h" +#include "hoomd/RandomNumbers.h" +#include "hoomd/TextureTools.h" + +/*! \file ActiveForceComputeGPU.cuh + \brief Declares GPU kernel code for calculating active forces forces on the GPU. Used by + ActiveForceComputeGPU. +*/ + +#ifndef __ACTIVE_FORCE_CONSTRAINT_COMPUTE_GPU_CUH__ +#define __ACTIVE_FORCE_CONSTRAINT_COMPUTE_GPU_CUH__ + +template +hipError_t gpu_compute_active_force_set_constraints(const unsigned int group_size, + unsigned int* d_index_array, + const Scalar4* d_pos, + Scalar4* d_orientation, + const Scalar4* d_f_act, + Manifold manifold, + unsigned int block_size); + +template +hipError_t gpu_compute_active_force_constraint_rotational_diffusion(const unsigned int group_size, + unsigned int* d_tag, + unsigned int* d_index_array, + const Scalar4* d_pos, + Scalar4* d_orientation, + Manifold manifold, + bool is2D, + const Scalar rotationDiff, + const uint64_t timestep, + const uint16_t seed, + unsigned int block_size); + +#ifdef __HIPCC__ + +//! Kernel for adjusting active force vectors to align parallel to an +// manifold surface constraint on the GPU +/*! \param group_size number of particles + \param d_index_array stores list to convert group index to global tag + \param d_pos particle positions on device + \param d_f_act particle active force unit vector + \param manifold constraint +*/ +template +__global__ void gpu_compute_active_force_set_constraints_kernel(const unsigned int group_size, + unsigned int* d_index_array, + const Scalar4* d_pos, + Scalar4* d_orientation, + const Scalar4* d_f_act, + Manifold manifold) + { + unsigned int group_idx = blockIdx.x * blockDim.x + threadIdx.x; + if (group_idx >= group_size) + return; + + unsigned int idx = d_index_array[group_idx]; + Scalar4 posidx = __ldg(d_pos + idx); + unsigned int type = __scalar_as_int(posidx.w); + + Scalar4 fact = __ldg(d_f_act + type); + + if (fact.w != 0) + { + Scalar3 current_pos = make_scalar3(posidx.x, posidx.y, posidx.z); + + vec3 norm = normalize(vec3(manifold.derivative(current_pos))); + + vec3 f(fact.x, fact.y, fact.z); + quat quati(__ldg(d_orientation + idx)); + vec3 fi = rotate(quati, f); + + Scalar dot_prod = fi.x * norm.x + fi.y * norm.y + fi.z * norm.z; + + Scalar dot_perp_prod = slow::rsqrt(1 - dot_prod * dot_prod); + + Scalar phi = slow::atan(dot_prod * dot_perp_prod); + + fi.x -= norm.x * dot_prod; + fi.y -= norm.y * dot_prod; + fi.z -= norm.z * dot_prod; + + Scalar new_norm = slow::rsqrt(fi.x * fi.x + fi.y * fi.y + fi.z * fi.z); + + fi *= new_norm; + + vec3 rot_vec = cross(norm, fi); + + quat rot_quat = quat::fromAxisAngle(rot_vec, phi); + + quati = rot_quat * quati; + + d_orientation[idx] = quat_to_scalar4(quati); + } + } + +//! Kernel for applying rotational diffusion to active force vectors on the GPU +/*! \param group_size number of particles + \param d_index_array stores list to convert group index to global tag + \param d_pos particle positions on device + \param manifold constraint + \param is2D check if simulation is 2D or 3D + \param rotationConst particle rotational diffusion constant + \param seed seed for random number generator +*/ +template +__global__ void +gpu_compute_active_force_constraint_rotational_diffusion_kernel(const unsigned int group_size, + unsigned int* d_tag, + unsigned int* d_index_array, + const Scalar4* d_pos, + Scalar4* d_orientation, + Manifold manifold, + bool is2D, + const Scalar rotationConst, + const uint64_t timestep, + const uint16_t seed) + { + unsigned int group_idx = blockIdx.x * blockDim.x + threadIdx.x; + if (group_idx >= group_size) + return; + + unsigned int idx = d_index_array[group_idx]; + Scalar4 posidx = __ldg(d_pos + idx); + unsigned int type = __scalar_as_int(posidx.w); + unsigned int ptag = d_tag[group_idx]; + + quat quati(__ldg(d_orientation + idx)); + + hoomd::RandomGenerator rng( + hoomd::Seed(hoomd::RNGIdentifier::ActiveForceCompute, timestep, seed), + hoomd::Counter(ptag)); + + Scalar3 current_pos = make_scalar3(posidx.x, posidx.y, posidx.z); + vec3 norm = normalize(vec3(manifold.derivative(current_pos))); + + Scalar delta_theta = hoomd::NormalDistribution(rotationConst)(rng); + + quat rot_quat = quat::fromAxisAngle(norm, delta_theta); + + quati = rot_quat * quati; + d_orientation[idx] = quat_to_scalar4(quati); + } + +template +hipError_t gpu_compute_active_force_set_constraints(const unsigned int group_size, + unsigned int* d_index_array, + const Scalar4* d_pos, + Scalar4* d_orientation, + const Scalar4* d_f_act, + Manifold manifold, + unsigned int block_size) + { + // setup the grid to run the kernel + dim3 grid(group_size / block_size + 1, 1, 1); + dim3 threads(block_size, 1, 1); + + // run the kernel + hipLaunchKernelGGL((gpu_compute_active_force_set_constraints_kernel), + dim3(grid), + dim3(threads), + 0, + 0, + group_size, + d_index_array, + d_pos, + d_orientation, + d_f_act, + manifold); + return hipSuccess; + } + +template +hipError_t gpu_compute_active_force_constraint_rotational_diffusion(const unsigned int group_size, + unsigned int* d_tag, + unsigned int* d_index_array, + const Scalar4* d_pos, + Scalar4* d_orientation, + Manifold manifold, + bool is2D, + const Scalar rotationConst, + const uint64_t timestep, + const uint16_t seed, + unsigned int block_size) + { + // setup the grid to run the kernel + dim3 grid(group_size / block_size + 1, 1, 1); + dim3 threads(block_size, 1, 1); + + // run the kernel + hipLaunchKernelGGL((gpu_compute_active_force_constraint_rotational_diffusion_kernel), + dim3(grid), + dim3(threads), + 0, + 0, + group_size, + d_tag, + d_index_array, + d_pos, + d_orientation, + manifold, + is2D, + rotationConst, + timestep, + seed); + return hipSuccess; + } + +#endif +#endif diff --git a/hoomd/md/ActiveForceConstraintComputeGPU.h b/hoomd/md/ActiveForceConstraintComputeGPU.h new file mode 100644 index 0000000000..499f6d0150 --- /dev/null +++ b/hoomd/md/ActiveForceConstraintComputeGPU.h @@ -0,0 +1,296 @@ +// Copyright (c) 2009-2021 The Regents of the University of Michigan +// This file is part of the HOOMD-blue project, released under the BSD 3-Clause License. + +// Maintainer: joaander + +#include "ActiveForceComputeGPU.cuh" +#include "ActiveForceConstraintCompute.h" +#include "ActiveForceConstraintComputeGPU.cuh" +#include "hoomd/Autotuner.h" + +/*! \file ActiveForceConstraintComputeGPU.h + \brief Declares a class for computing active forces on the GPU +*/ + +#ifdef __HIPCC__ +#error This header cannot be compiled by nvcc +#endif + +#include + +#ifndef __ACTIVEFORCECONSTRAINTCOMPUTE_GPU_H__ +#define __ACTIVEFORCECONSTRAINTCOMPUTE_GPU_H__ + +#include +namespace py = pybind11; +using namespace std; + +//! Adds an active force to a number of particles with confinement on the GPU +/*! \ingroup computes + */ +template +class PYBIND11_EXPORT ActiveForceConstraintComputeGPU + : public ActiveForceConstraintCompute + { + public: + //! Constructs the compute + ActiveForceConstraintComputeGPU(std::shared_ptr sysdef, + std::shared_ptr group, + Scalar rotation_diff, + Manifold manifold); + + //! Set autotuner parameters + /*! \param enable Enable/disable autotuning + \param period period (approximate) in time steps when returning occurs + */ + virtual void setAutotunerParams(bool enable, unsigned int period) + { + ActiveForceConstraintCompute::setAutotunerParams(enable, period); + m_tuner_force->setPeriod(period); + m_tuner_force->setEnabled(enable); + m_tuner_diffusion->setPeriod(period); + m_tuner_diffusion->setEnabled(enable); + m_tuner_constraint->setPeriod(period); + m_tuner_constraint->setEnabled(enable); + } + + protected: + std::unique_ptr m_tuner_force; //!< Autotuner for block size (force kernel) + std::unique_ptr m_tuner_diffusion; //!< Autotuner for block size (diff kernel) + std::unique_ptr m_tuner_constraint; //!< Autotuner for block size (constr kernel) + + //! Set forces for particles + virtual void setForces(); + + //! Orientational diffusion for spherical particles + virtual void rotationalDiffusion(uint64_t timestep); + + //! Set constraints if particles confined to a surface + virtual void setConstraint(); + }; + +/*! \file ActiveForceConstraintComputeGPU.cc + \brief Contains code for the ActiveForceConstraintComputeGPU class +*/ + +/*! \param f_list An array of (x,y,z) tuples for the active force vector for each + individual particle. + \param orientation_link if True then forces and torques are applied in the + particle's reference frame. If false, then the box reference frame is + used. Only relevant for non-point-like anisotropic particles. + \param orientation_reverse_link When True, the particle's orientation is set + to match the active force vector. Useful for using a particle's + orientation to log the active force vector. Not recommended for + anisotropic particles + \param rotation_diff rotational diffusion constant for all particles. + \param manifold specifies a manfold surface, to which particles are confined. +*/ +template +ActiveForceConstraintComputeGPU::ActiveForceConstraintComputeGPU( + std::shared_ptr sysdef, + std::shared_ptr group, + Scalar rotation_diff, + Manifold manifold) + : ActiveForceConstraintCompute(sysdef, group, rotation_diff, manifold) + { + if (!this->m_exec_conf->isCUDAEnabled()) + { + this->m_exec_conf->msg->error() << "Creating a ActiveForceConstraintComputeGPU with no GPU " + "in the execution configuration" + << endl; + throw std::runtime_error("Error initializing ActiveForceConstraintComputeGPU"); + } + + // initialize autotuner + std::vector valid_params; + unsigned int warp_size = this->m_exec_conf->dev_prop.warpSize; + for (unsigned int block_size = warp_size; block_size <= 1024; block_size += warp_size) + valid_params.push_back(block_size); + + m_tuner_force.reset( + new Autotuner(valid_params, 5, 100000, "active_constraint_force", this->m_exec_conf)); + m_tuner_diffusion.reset( + new Autotuner(valid_params, 5, 100000, "active_constraint_diffusion", this->m_exec_conf)); + m_tuner_constraint.reset( + new Autotuner(valid_params, 5, 100000, "active_constraint_constraint", this->m_exec_conf)); + + unsigned int type = this->m_pdata->getNTypes(); + GlobalVector tmp_f_activeVec(type, this->m_exec_conf); + GlobalVector tmp_t_activeVec(type, this->m_exec_conf); + + { + ArrayHandle old_f_activeVec(this->m_f_activeVec, access_location::host); + ArrayHandle old_t_activeVec(this->m_t_activeVec, access_location::host); + + ArrayHandle f_activeVec(tmp_f_activeVec, access_location::host); + ArrayHandle t_activeVec(tmp_t_activeVec, access_location::host); + + // for each type of the particles in the group + for (unsigned int i = 0; i < type; i++) + { + f_activeVec.data[i] = old_f_activeVec.data[i]; + + t_activeVec.data[i] = old_t_activeVec.data[i]; + } + + this->last_computed = 10; + } + + this->m_f_activeVec.swap(tmp_f_activeVec); + this->m_t_activeVec.swap(tmp_t_activeVec); + } + +/*! This function sets appropriate active forces and torques on all active particles. + */ +template void ActiveForceConstraintComputeGPU::setForces() + { + // array handles + ArrayHandle d_f_actVec(this->m_f_activeVec, + access_location::device, + access_mode::read); + ArrayHandle d_force(this->m_force, access_location::device, access_mode::overwrite); + + ArrayHandle d_t_actVec(this->m_t_activeVec, + access_location::device, + access_mode::read); + ArrayHandle d_torque(this->m_torque, access_location::device, access_mode::overwrite); + + ArrayHandle d_pos(this->m_pdata->getPositions(), + access_location::device, + access_mode::read); + ArrayHandle d_orientation(this->m_pdata->getOrientationArray(), + access_location::device, + access_mode::read); + ArrayHandle d_index_array(this->m_group->getIndexArray(), + access_location::device, + access_mode::read); + + // sanity check + assert(d_force.data != NULL); + assert(d_f_actVec.data != NULL); + assert(d_t_actVec.data != NULL); + assert(d_pos.data != NULL); + assert(d_orientation.data != NULL); + assert(d_index_array.data != NULL); + unsigned int group_size = this->m_group->getNumMembers(); + unsigned int N = this->m_pdata->getN(); + + // compute the forces on the GPU + this->m_tuner_force->begin(); + gpu_compute_active_force_set_forces(group_size, + d_index_array.data, + d_force.data, + d_torque.data, + d_pos.data, + d_orientation.data, + d_f_actVec.data, + d_t_actVec.data, + N, + this->m_tuner_force->getParam()); + + if (this->m_exec_conf->isCUDAErrorCheckingEnabled()) + CHECK_CUDA_ERROR(); + + this->m_tuner_force->end(); + } + +/*! This function applies rotational diffusion to all active particles. The angle between the torque + vector and + * force vector does not change + \param timestep Current timestep +*/ +template +void ActiveForceConstraintComputeGPU::rotationalDiffusion(uint64_t timestep) + { + // array handles + ArrayHandle d_pos(this->m_pdata->getPositions(), + access_location::device, + access_mode::read); + ArrayHandle d_orientation(this->m_pdata->getOrientationArray(), + access_location::device, + access_mode::readwrite); + ArrayHandle d_index_array(this->m_group->getIndexArray(), + access_location::device, + access_mode::read); + ArrayHandle d_tag(this->m_pdata->getTags(), + access_location::device, + access_mode::read); + + assert(d_pos.data != NULL); + + bool is2D = (this->m_sysdef->getNDimensions() == 2); + unsigned int group_size = this->m_group->getNumMembers(); + + // perform the update on the GPU + this->m_tuner_diffusion->begin(); + + gpu_compute_active_force_constraint_rotational_diffusion( + group_size, + d_tag.data, + d_index_array.data, + d_pos.data, + d_orientation.data, + this->m_manifold, + is2D, + this->m_rotationConst, + timestep, + this->m_sysdef->getSeed(), + this->m_tuner_diffusion->getParam()); + + if (this->m_exec_conf->isCUDAErrorCheckingEnabled()) + CHECK_CUDA_ERROR(); + + this->m_tuner_diffusion->end(); + } + +/*! This function sets an ellipsoid surface constraint for all active particles + */ +template void ActiveForceConstraintComputeGPU::setConstraint() + { + // array handles + ArrayHandle d_f_actVec(this->m_f_activeVec, + access_location::device, + access_mode::readwrite); + ArrayHandle d_pos(this->m_pdata->getPositions(), + access_location::device, + access_mode::read); + ArrayHandle d_orientation(this->m_pdata->getOrientationArray(), + access_location::device, + access_mode::readwrite); + ArrayHandle d_index_array(this->m_group->getIndexArray(), + access_location::device, + access_mode::read); + + assert(d_pos.data != NULL); + + unsigned int group_size = this->m_group->getNumMembers(); + + // perform the update on the GPU + this->m_tuner_constraint->begin(); + + gpu_compute_active_force_set_constraints(group_size, + d_index_array.data, + d_pos.data, + d_orientation.data, + d_f_actVec.data, + this->m_manifold, + this->m_tuner_constraint->getParam()); + + if (this->m_exec_conf->isCUDAErrorCheckingEnabled()) + CHECK_CUDA_ERROR(); + + this->m_tuner_constraint->end(); + } + +template +void export_ActiveForceConstraintComputeGPU(py::module& m, const std::string& name) + { + py::class_, + ActiveForceConstraintCompute, + std::shared_ptr>>(m, name.c_str()) + .def(py::init, + std::shared_ptr, + Scalar, + Manifold>()); + } +#endif diff --git a/hoomd/md/CMakeLists.txt b/hoomd/md/CMakeLists.txt index 04cbe1771a..283e3e6f2d 100644 --- a/hoomd/md/CMakeLists.txt +++ b/hoomd/md/CMakeLists.txt @@ -6,10 +6,7 @@ set(_md_sources module-md.cc CommunicatorGrid.cc ComputeThermo.cc ComputeThermoHMA.cc - ConstraintEllipsoid.cc - ConstraintSphere.cc CosineSqAngleForceCompute.cc - OneDConstraint.cc FIREEnergyMinimizer.cc ForceComposite.cc ForceDistanceConstraint.cc @@ -47,6 +44,9 @@ set(_md_sources module-md.cc set(_md_headers ActiveForceComputeGPU.h ActiveForceCompute.h + ActiveForceConstraintCompute.h + ActiveForceConstraintComputeGPU.h + ActiveForceConstraintComputeGPU.cuh AllAnisoPairPotentials.h AllBondPotentials.h AllExternalPotentials.h @@ -68,10 +68,6 @@ set(_md_headers ActiveForceComputeGPU.h ComputeThermoHMA.h ComputeThermoTypes.h ComputeThermoHMATypes.h - ConstraintEllipsoidGPU.h - ConstraintEllipsoid.h - ConstraintSphereGPU.h - ConstraintSphere.h CosineSqAngleForceComputeGPU.h CosineSqAngleForceCompute.h EvaluatorBondFENE.h @@ -79,9 +75,6 @@ set(_md_headers ActiveForceComputeGPU.h EvaluatorBondTether.h EvaluatorSpecialPairLJ.h EvaluatorSpecialPairCoulomb.h - EvaluatorConstraintEllipsoid.h - EvaluatorConstraint.h - EvaluatorConstraintSphere.h EvaluatorExternalElectricField.h EvaluatorExternalPeriodic.h EvaluatorPairBuckingham.h @@ -199,9 +192,6 @@ list(APPEND _md_sources ActiveForceComputeGPU.cc CommunicatorGridGPU.cc ComputeThermoGPU.cc ComputeThermoHMAGPU.cc - ConstraintEllipsoidGPU.cc - ConstraintSphereGPU.cc - OneDConstraintGPU.cc FIREEnergyMinimizerGPU.cc ForceCompositeGPU.cc ForceDistanceConstraintGPU.cc @@ -258,9 +248,6 @@ set(_md_cu_sources ActiveForceComputeGPU.cu ZBLDriverPotentialPairGPU.cu BondTablePotentialGPU.cu CommunicatorGridGPU.cu - ConstraintEllipsoidGPU.cu - ConstraintSphereGPU.cu - OneDConstraintGPU.cu DriverTersoffGPU.cu FIREEnergyMinimizerGPU.cu ForceCompositeGPU.cu @@ -281,20 +268,27 @@ set(_md_cu_sources ActiveForceComputeGPU.cu TwoStepBDGPU.cu TwoStepBerendsenGPU.cu TwoStepLangevinGPU.cu - TwoStepRATTLELangevinGPU.cu + TwoStepRATTLELangevinGPU.cu TwoStepNPTMTKGPU.cu TwoStepNVEGPU.cu - TwoStepRATTLENVEGPU.cu + TwoStepRATTLENVEGPU.cu TwoStepNVTMTKGPU.cu MuellerPlatheFlowGPU.cu CosineSqAngleForceGPU.cu - all_kernels_diamond_manifold.cu - all_kernels_ellipsoid_manifold.cu - all_kernels_gyroid_manifold.cu - all_kernels_primitive_manifold.cu - all_kernels_sphere_manifold.cu - all_kernels_xyplane_manifold.cu - all_kernels_zcylinder_manifold.cu + all_kernels_diamond_manifold.cu + all_kernels_ellipsoid_manifold.cu + all_kernels_gyroid_manifold.cu + all_kernels_primitive_manifold.cu + all_kernels_sphere_manifold.cu + all_kernels_xyplane_manifold.cu + all_kernels_zcylinder_manifold.cu + all_kernels_active_diamond_manifold.cu + all_kernels_active_ellipsoid_manifold.cu + all_kernels_active_gyroid_manifold.cu + all_kernels_active_primitive_manifold.cu + all_kernels_active_sphere_manifold.cu + all_kernels_active_xyplane_manifold.cu + all_kernels_active_zcylinder_manifold.cu ) if (ENABLE_HIP) diff --git a/hoomd/md/ConstraintEllipsoid.cc b/hoomd/md/ConstraintEllipsoid.cc deleted file mode 100644 index 7ca083de21..0000000000 --- a/hoomd/md/ConstraintEllipsoid.cc +++ /dev/null @@ -1,168 +0,0 @@ -// Copyright (c) 2009-2021 The Regents of the University of Michigan -// This file is part of the HOOMD-blue project, released under the BSD 3-Clause License. - -// Maintainer: joaander - -#include "ConstraintEllipsoid.h" -#include "EvaluatorConstraintEllipsoid.h" - -namespace py = pybind11; - -using namespace std; - -/*! \file ConstraintEllipsoid.cc - \brief Contains code for the ConstraintEllipsoid class -*/ - -/*! \param sysdef SystemDefinition containing the ParticleData to compute forces on - \param group Group of particles on which to apply this constraint - \param P position of the Ellipsoid - \param rx radius of the Ellipsoid in the X direction - \param ry radius of the Ellipsoid in the Y direction - \param rz radius of the Ellipsoid in the Z direction - NOTE: For the algorithm to work, we must have _rx >= _rz, ry >= _rz, and _rz > 0. -*/ -ConstraintEllipsoid::ConstraintEllipsoid(std::shared_ptr sysdef, - std::shared_ptr group, - Scalar3 P, - Scalar rx, - Scalar ry, - Scalar rz) - : Updater(sysdef), m_group(group), m_P(P), m_rx(rx), m_ry(ry), m_rz(rz) - { - m_exec_conf->msg->notice(5) << "Constructing ConstraintEllipsoid" << endl; - - validate(); - } - -ConstraintEllipsoid::~ConstraintEllipsoid() - { - m_exec_conf->msg->notice(5) << "Destroying ConstraintEllipsoid" << endl; - } - -/*! Computes the specified constraint forces - \param timestep Current timestep -*/ -void ConstraintEllipsoid::update(uint64_t timestep) - { - Updater::update(timestep); - unsigned int group_size = m_group->getNumMembers(); - if (group_size == 0) - return; - - if (m_prof) - m_prof->push("ConstraintEllipsoid"); - - assert(m_pdata); - - // access the particle data arrays - ArrayHandle h_pos(m_pdata->getPositions(), - access_location::host, - access_mode::readwrite); - - EvaluatorConstraintEllipsoid Ellipsoid(m_P, m_rx, m_ry, m_rz); - // for each of the particles in the group - for (unsigned int group_idx = 0; group_idx < group_size; group_idx++) - { - // get the current particle properties - unsigned int j = m_group->getMemberIndex(group_idx); - Scalar3 X = make_scalar3(h_pos.data[j].x, h_pos.data[j].y, h_pos.data[j].z); - - // evaluate the constraint position - Scalar3 C = Ellipsoid.evalClosest(X); - - // apply the constraint - h_pos.data[j].x = C.x; - h_pos.data[j].y = C.y; - h_pos.data[j].z = C.z; - } - - if (m_prof) - m_prof->pop(); - } - -/*! Print warning messages if the Ellipsoid is outside the box. - Generate an error if any particle in the group is not near the Ellipsoid. -*/ -void ConstraintEllipsoid::validate() - { - BoxDim box = m_pdata->getBox(); - Scalar3 lo = box.getLo(); - Scalar3 hi = box.getHi(); - - if (m_P.x + m_rx > hi.x || m_P.x - m_rx < lo.x || m_P.y + m_ry > hi.y || m_P.y - m_ry < lo.y - || m_P.z + m_rz > hi.z || m_P.z - m_rz < lo.z) - { - m_exec_conf->msg->warning() << "constrain.ellipsoid: ellipsoid constraint is outside of " - "the box. Constrained particle positions may be incorrect" - << endl; - } - - if (m_rx == 0 || m_ry == 0 || m_rz == 0) - { - m_exec_conf->msg->warning() << "constrain.ellipsoid: one of the ellipsoid dimensions is 0. " - "Constraint may be incorrect." - << endl; - } - - unsigned int group_size = m_group->getNumMembers(); - if (group_size == 0) - return; - - ArrayHandle h_pos(m_pdata->getPositions(), access_location::host, access_mode::read); - ArrayHandle h_tag(m_pdata->getTags(), access_location::host, access_mode::read); - ArrayHandle h_body(m_pdata->getBodies(), - access_location::host, - access_mode::read); - - EvaluatorConstraintEllipsoid Ellipsoid(m_P, m_rx, m_ry, m_rz); - // for each of the particles in the group - bool errors = false; - for (unsigned int group_idx = 0; group_idx < group_size; group_idx++) - { - // get the current particle properties - unsigned int j = m_group->getMemberIndex(group_idx); - - Scalar3 X = make_scalar3(h_pos.data[j].x, h_pos.data[j].y, h_pos.data[j].z); - Scalar3 C = Ellipsoid.evalClosest(X); - Scalar3 V; - V.x = C.x - X.x; - V.y = C.y - X.y; - V.z = C.z - X.z; - Scalar dist = slow::sqrt(V.x * V.x + V.y * V.y + V.z * V.z); - - if (dist > Scalar(1.0)) - { - m_exec_conf->msg->error() - << "constrain.ellipsoid: Particle " << h_tag.data[j] << " is more than 1 unit of" - << " distance away from the closest point on the ellipsoid constraint" << endl; - errors = true; - } - - if (h_body.data[j] != NO_BODY) - { - m_exec_conf->msg->error() - << "constrain.ellipsoid: Particle " << h_tag.data[j] << " belongs to a rigid body" - << " - cannot constrain" << endl; - errors = true; - } - } - - if (errors) - { - throw std::runtime_error("Invalid constraint specified"); - } - } - -void export_ConstraintEllipsoid(py::module& m) - { - py::class_>( - m, - "ConstraintEllipsoid") - .def(py::init, - std::shared_ptr, - Scalar3, - Scalar, - Scalar, - Scalar>()); - } diff --git a/hoomd/md/ConstraintEllipsoid.h b/hoomd/md/ConstraintEllipsoid.h deleted file mode 100644 index 2604bc37a5..0000000000 --- a/hoomd/md/ConstraintEllipsoid.h +++ /dev/null @@ -1,59 +0,0 @@ -// Copyright (c) 2009-2021 The Regents of the University of Michigan -// This file is part of the HOOMD-blue project, released under the BSD 3-Clause License. - -// Maintainer: joaander - -#include "hoomd/ParticleGroup.h" -#include "hoomd/Updater.h" -#include - -/*! \file ConstraintEllipsoid.h - \brief Declares a class for computing ellipsoid constraint forces -*/ - -#ifdef __HIPCC__ -#error This header cannot be compiled by nvcc -#endif - -#include - -#ifndef __CONSTRAINT_Ellipsoid_H__ -#define __CONSTRAINT_Ellipsoid_H__ - -//! Applies a constraint force to keep a group of particles on a Ellipsoid -/*! \ingroup computes - */ -class PYBIND11_EXPORT ConstraintEllipsoid : public Updater - { - public: - //! Constructs the compute - ConstraintEllipsoid(std::shared_ptr sysdef, - std::shared_ptr group, - Scalar3 P, - Scalar rx, - Scalar ry, - Scalar rz); - - //! Destructor - virtual ~ConstraintEllipsoid(); - - //! Take one timestep forward - virtual void update(uint64_t timestep); - - protected: - std::shared_ptr - m_group; //!< Group of particles on which this constraint is applied - Scalar3 m_P; //!< Position of the Ellipsoid - Scalar m_rx; //!< Radius in X direction of the Ellipsoid - Scalar m_ry; //!< Radius in Y direction of the Ellipsoid - Scalar m_rz; //!< Radius in Z direction of the Ellipsoid - - private: - //! Validate that the ellipsoid is in the box and all particles are very near the constraint - void validate(); - }; - -//! Exports the ConstraintEllipsoid class to python -void export_ConstraintEllipsoid(pybind11::module& m); - -#endif diff --git a/hoomd/md/ConstraintEllipsoidGPU.cc b/hoomd/md/ConstraintEllipsoidGPU.cc deleted file mode 100644 index fcbd92c5b2..0000000000 --- a/hoomd/md/ConstraintEllipsoidGPU.cc +++ /dev/null @@ -1,96 +0,0 @@ -// Copyright (c) 2009-2021 The Regents of the University of Michigan -// This file is part of the HOOMD-blue project, released under the BSD 3-Clause License. - -// Maintainer: joaander - -#include "ConstraintEllipsoidGPU.h" -#include "ConstraintEllipsoidGPU.cuh" - -namespace py = pybind11; - -using namespace std; - -/*! \file ConstraintEllipsoidGPU.cc - \brief Contains code for the ConstraintEllipsoidGPU class -*/ - -/*! \param sysdef SystemDefinition containing the ParticleData to compute forces on - \param group Group of particles on which to apply this constraint - \param P position of the Ellipsoid - \param rx radius of the Ellipsoid in the X direction - \param ry radius of the Ellipsoid in the Y direction - \param rz radius of the Ellipsoid in the Z direction - NOTE: For the algorithm to work, we must have _rx >= _rz, ry >= _rz, and _rz > 0. -*/ -ConstraintEllipsoidGPU::ConstraintEllipsoidGPU(std::shared_ptr sysdef, - std::shared_ptr group, - Scalar3 P, - Scalar rx, - Scalar ry, - Scalar rz) - : ConstraintEllipsoid(sysdef, group, P, rx, ry, rz), m_block_size(256) - { - if (!m_exec_conf->isCUDAEnabled()) - { - m_exec_conf->msg->error() - << "Creating a ConstraintEllipsoidGPU with no GPU in the execution configuration" - << endl; - throw std::runtime_error("Error initializing ConstraintEllipsoidGPU"); - } - } - -/*! Computes the specified constraint forces - \param timestep Current timestep -*/ -void ConstraintEllipsoidGPU::update(uint64_t timestep) - { - Updater::update(timestep); - unsigned int group_size = m_group->getNumMembers(); - if (group_size == 0) - return; - - if (m_prof) - m_prof->push("ConstraintEllipsoid"); - - assert(m_pdata); - - // access the particle data arrays - const GlobalArray& group_members = m_group->getIndexArray(); - ArrayHandle d_group_members(group_members, - access_location::device, - access_mode::read); - - ArrayHandle d_pos(m_pdata->getPositions(), - access_location::device, - access_mode::readwrite); - - // run the kernel in parallel on all GPUs - gpu_compute_constraint_ellipsoid_constraint(d_group_members.data, - m_group->getNumMembers(), - m_pdata->getN(), - d_pos.data, - m_P, - m_rx, - m_ry, - m_rz, - m_block_size); - - if (m_exec_conf->isCUDAErrorCheckingEnabled()) - CHECK_CUDA_ERROR(); - - if (m_prof) - m_prof->pop(m_exec_conf); - } - -void export_ConstraintEllipsoidGPU(py::module& m) - { - py::class_>(m, "ConstraintEllipsoidGPU") - .def(py::init, - std::shared_ptr, - Scalar3, - Scalar, - Scalar, - Scalar>()); - } diff --git a/hoomd/md/ConstraintEllipsoidGPU.cu b/hoomd/md/ConstraintEllipsoidGPU.cu deleted file mode 100644 index a8409f553f..0000000000 --- a/hoomd/md/ConstraintEllipsoidGPU.cu +++ /dev/null @@ -1,108 +0,0 @@ -#include "hip/hip_runtime.h" -// Copyright (c) 2009-2021 The Regents of the University of Michigan -// This file is part of the HOOMD-blue project, released under the BSD 3-Clause License. - -// Maintainer: joaander - -#include "ConstraintEllipsoidGPU.cuh" -#include "EvaluatorConstraint.h" -#include "EvaluatorConstraintEllipsoid.h" - -#include - -/*! \file ConstraintEllipsoidGPU.cu - \brief Defines GPU kernel code for calculating ellipsoid constraint forces. Used by - ConstraintEllipsoidGPU. -*/ - -//! Kernel for calculating ellipsoid constraint forces on the GPU -/*! \param d_group_members List of members in the group - \param group_size number of members in the group - \param N number of particles in system - \param d_pos particle positions on device - \param P Position of the ellipsoid - \param rx radius of the ellipsoid in x direction - \param ry radius of the ellipsoid in y direction - \param rz radius of the ellipsoid in z direction - \param deltaT step size from the Integrator -*/ -extern "C" __global__ void -gpu_compute_constraint_ellipsoid_constraint_kernel(const unsigned int* d_group_members, - unsigned int group_size, - const unsigned int N, - Scalar4* d_pos, - Scalar3 P, - Scalar rx, - Scalar ry, - Scalar rz) - { - // start by identifying which particle we are to handle - // determine which particle this thread works on - int group_idx = blockIdx.x * blockDim.x + threadIdx.x; - - if (group_idx >= group_size) - return; - - unsigned int idx = d_group_members[group_idx]; - - // read in position, velocity, net force, and mass - Scalar4 pos = d_pos[idx]; - - // convert to Scalar3's for passing to the evaluators - Scalar3 X = make_scalar3(pos.x, pos.y, pos.z); - - // evaluate the constraint position - EvaluatorConstraintEllipsoid Ellipsoid(P, rx, ry, rz); - Scalar3 C = Ellipsoid.evalClosest(X); - - // apply the constraint - d_pos[idx] = make_scalar4(C.x, C.y, C.z, Scalar(0.0)); - } - -/*! \param d_group_members List of members in the group - \param group_size number of members in the group - \param N number of particles - \param d_pos particle positions on the device - \param P Position of the ellipsoid - \param rx radius of the ellipsoid in x direction - \param ry radius of the ellipsoid in y direction - \param rz radius of the ellipsoid in z direction - \param deltaT step size from the Integrator - \param block_size Block size to execute on the GPU - - \returns Any error code resulting from the kernel launch - \note Always returns hipSuccess in release builds to avoid the hipDeviceSynchronize() -*/ -hipError_t gpu_compute_constraint_ellipsoid_constraint(const unsigned int* d_group_members, - unsigned int group_size, - const unsigned int N, - Scalar4* d_pos, - const Scalar3 P, - Scalar rx, - Scalar ry, - Scalar rz, - unsigned int block_size) - { - assert(d_group_members); - - // setup the grid to run the kernel - dim3 grid(group_size / block_size + 1, 1, 1); - dim3 threads(block_size, 1, 1); - - // run the kernel - hipLaunchKernelGGL((gpu_compute_constraint_ellipsoid_constraint_kernel), - dim3(grid), - dim3(threads), - 0, - 0, - d_group_members, - group_size, - N, - d_pos, - P, - rx, - ry, - rz); - - return hipSuccess; - } diff --git a/hoomd/md/ConstraintEllipsoidGPU.cuh b/hoomd/md/ConstraintEllipsoidGPU.cuh deleted file mode 100644 index b18b73ee99..0000000000 --- a/hoomd/md/ConstraintEllipsoidGPU.cuh +++ /dev/null @@ -1,28 +0,0 @@ -// Copyright (c) 2009-2021 The Regents of the University of Michigan -// This file is part of the HOOMD-blue project, released under the BSD 3-Clause License. - -// Maintainer: joaander - -#include "hoomd/HOOMDMath.h" -#include "hoomd/ParticleData.cuh" - -/*! \file ConstraintEllipsoidGPU.cuh - \brief Declares GPU kernel code for calculating sphere constraint forces on the GPU. Used by - ConstraintEllipsoidGPU. -*/ - -#ifndef __CONSTRAINT_ELLIPSOID_GPU_CUH__ -#define __CONSTRAINT_ELLIPSOID_GPU_CUH__ - -//! Kernel driver that computes harmonic bond forces for HarmonicBondForceComputeGPU -hipError_t gpu_compute_constraint_ellipsoid_constraint(const unsigned int* d_group_members, - unsigned int group_size, - const unsigned int N, - Scalar4* d_pos, - const Scalar3 P, - Scalar rx, - Scalar ry, - Scalar rz, - unsigned int block_size); - -#endif diff --git a/hoomd/md/ConstraintEllipsoidGPU.h b/hoomd/md/ConstraintEllipsoidGPU.h deleted file mode 100644 index 8b7c04c164..0000000000 --- a/hoomd/md/ConstraintEllipsoidGPU.h +++ /dev/null @@ -1,45 +0,0 @@ -// Copyright (c) 2009-2021 The Regents of the University of Michigan -// This file is part of the HOOMD-blue project, released under the BSD 3-Clause License. - -// Maintainer: joaander - -#include "ConstraintEllipsoid.h" - -/*! \file ConstraintEllipsoidGPU.h - \brief Declares a class for computing ellipsoid constraint forces -*/ - -#ifdef __HIPCC__ -#error This header cannot be compiled by nvcc -#endif - -#include - -#ifndef __CONSTRAINT_ELLIPSOID_GPU_H__ -#define __CONSTRAINT_ELLIPSOID_GPU_H__ - -//! Applys a constraint force to keep a group of particles on a Ellipsoid -/*! \ingroup computes - */ -class PYBIND11_EXPORT ConstraintEllipsoidGPU : public ConstraintEllipsoid - { - public: - //! Constructs the compute - ConstraintEllipsoidGPU(std::shared_ptr sysdef, - std::shared_ptr group, - Scalar3 P, - Scalar rx, - Scalar ry, - Scalar rz); - - protected: - unsigned int m_block_size; //!< block size to execute on the GPU - - //! Take one timestep forward - virtual void update(uint64_t timestep); - }; - -//! Exports the ConstraintEllipsoidGPU class to python -void export_ConstraintEllipsoidGPU(pybind11::module& m); - -#endif diff --git a/hoomd/md/ConstraintSphere.cc b/hoomd/md/ConstraintSphere.cc deleted file mode 100644 index f6b138c2ce..0000000000 --- a/hoomd/md/ConstraintSphere.cc +++ /dev/null @@ -1,198 +0,0 @@ -// Copyright (c) 2009-2021 The Regents of the University of Michigan -// This file is part of the HOOMD-blue project, released under the BSD 3-Clause License. - -// Maintainer: joaander - -#include "ConstraintSphere.h" -#include "EvaluatorConstraint.h" -#include "EvaluatorConstraintSphere.h" - -namespace py = pybind11; - -using namespace std; - -/*! \file ConstraintSphere.cc - \brief Contains code for the ConstraintSphere class -*/ - -/*! \param sysdef SystemDefinition containing the ParticleData to compute forces on - \param group Group of particles on which to apply this constraint - \param P position of the sphere - \param r radius of the sphere -*/ -ConstraintSphere::ConstraintSphere(std::shared_ptr sysdef, - std::shared_ptr group, - Scalar3 P, - Scalar r) - : ForceConstraint(sysdef), m_group(group), m_P(P), m_r(r) - { - m_exec_conf->msg->notice(5) << "Constructing ConstraintSphere" << endl; - - validate(); - } - -ConstraintSphere::~ConstraintSphere() - { - m_exec_conf->msg->notice(5) << "Destroying ConstraintSphere" << endl; - } - -/*! - \param P position of the sphere - \param r radius of the sphere -*/ -void ConstraintSphere::setSphere(Scalar3 P, Scalar r) - { - m_P = P; - m_r = r; - validate(); - } - -Scalar ConstraintSphere::getNDOFRemoved(std::shared_ptr query) - { - return m_group->intersectionSize(query); - } - -/*! Computes the specified constraint forces - \param timestep Current timestep -*/ -void ConstraintSphere::computeForces(uint64_t timestep) - { - unsigned int group_size = m_group->getNumMembers(); - if (group_size == 0) - return; - - if (m_prof) - m_prof->push("ConstraintSphere"); - - assert(m_pdata); - // access the particle data arrays - ArrayHandle h_pos(m_pdata->getPositions(), access_location::host, access_mode::read); - ArrayHandle h_vel(m_pdata->getVelocities(), access_location::host, access_mode::read); - - const GlobalArray& net_force = m_pdata->getNetForce(); - ArrayHandle h_net_force(net_force, access_location::host, access_mode::read); - - ArrayHandle h_force(m_force, access_location::host, access_mode::overwrite); - ArrayHandle h_virial(m_virial, access_location::host, access_mode::overwrite); - size_t virial_pitch = m_virial.getPitch(); - - // Zero data for force calculation. - memset((void*)h_force.data, 0, sizeof(Scalar4) * m_force.getNumElements()); - memset((void*)h_virial.data, 0, sizeof(Scalar) * m_virial.getNumElements()); - - // there are enough other checks on the input data: but it doesn't hurt to be safe - assert(h_force.data); - assert(h_virial.data); - - // for each of the particles in the group - for (unsigned int group_idx = 0; group_idx < group_size; group_idx++) - { - // get the current particle properties - unsigned int j = m_group->getMemberIndex(group_idx); - Scalar3 X = make_scalar3(h_pos.data[j].x, h_pos.data[j].y, h_pos.data[j].z); - Scalar3 V = make_scalar3(h_vel.data[j].x, h_vel.data[j].y, h_vel.data[j].z); - Scalar3 F - = make_scalar3(h_net_force.data[j].x, h_net_force.data[j].y, h_net_force.data[j].z); - Scalar m = h_vel.data[j].w; - - // evaluate the constraint position - EvaluatorConstraint constraint(X, V, F, m, m_deltaT); - EvaluatorConstraintSphere sphere(m_P, m_r); - Scalar3 C = sphere.evalClosest(constraint.evalU()); - - // evaluate the constraint force - Scalar3 FC; - Scalar virial[6]; - constraint.evalConstraintForce(FC, virial, C); - - // apply the constraint force - h_force.data[j].x = FC.x; - h_force.data[j].y = FC.y; - h_force.data[j].z = FC.z; - for (int k = 0; k < 6; k++) - h_virial.data[k * virial_pitch + j] = virial[k]; - } - - if (m_prof) - m_prof->pop(); - } - -/*! Print warning messages if the sphere is outside the box. - Generate an error if any particle in the group is not near the sphere. -*/ -void ConstraintSphere::validate() - { - BoxDim box = m_pdata->getGlobalBox(); - Scalar3 lo = box.getLo(); - Scalar3 hi = box.getHi(); - - if (m_P.x + m_r > hi.x || m_P.x - m_r < lo.x || m_P.y + m_r > hi.y || m_P.y - m_r < lo.y - || m_P.z + m_r > hi.z || m_P.z - m_r < lo.z) - { - m_exec_conf->msg->warning() << "constrain.sphere: Sphere constraint is outside of the box. " - "Constrained particle positions may be incorrect" - << endl; - } - - unsigned int group_size = m_group->getNumMembers(); - if (group_size == 0) - return; - - ArrayHandle h_pos(m_pdata->getPositions(), access_location::host, access_mode::read); - ArrayHandle h_tag(m_pdata->getTags(), access_location::host, access_mode::read); - ArrayHandle h_body(m_pdata->getBodies(), - access_location::host, - access_mode::read); - - // for each of the particles in the group - bool errors = false; - for (unsigned int group_idx = 0; group_idx < group_size; group_idx++) - { - // get the current particle properties - unsigned int j = m_group->getMemberIndex(group_idx); - Scalar3 X = make_scalar3(h_pos.data[j].x, h_pos.data[j].y, h_pos.data[j].z); - - // evaluate the constraint position - EvaluatorConstraintSphere sphere(m_P, m_r); - Scalar3 C = sphere.evalClosest(X); - Scalar3 V; - V.x = C.x - X.x; - V.y = C.y - X.y; - V.z = C.z - X.z; - Scalar dist = sqrt(V.x * V.x + V.y * V.y + V.z * V.z); - - if (dist > Scalar(1.0)) - { - m_exec_conf->msg->error() - << "constrain.sphere: Particle " << h_tag.data[j] << " is more than 1 unit of" - << " distance away from the closest point on the sphere constraint" << endl; - errors = true; - } - - if (h_body.data[j] < MIN_FLOPPY) - { - m_exec_conf->msg->error() - << "constrain.sphere: Particle " << h_tag.data[j] << " belongs to a rigid body" - << " - cannot constrain" << endl; - errors = true; - } - } - - if (errors) - { - throw std::runtime_error("Invalid constraint specified"); - } - } - -void export_ConstraintSphere(py::module& m) - { - py::class_>( - m, - "ConstraintSphere") - .def(py::init, - std::shared_ptr, - Scalar3, - Scalar>()) - .def("setSphere", &ConstraintSphere::setSphere) - .def("getNDOFRemoved", &ConstraintSphere::getNDOFRemoved); - } diff --git a/hoomd/md/ConstraintSphere.h b/hoomd/md/ConstraintSphere.h deleted file mode 100644 index a35de02fd1..0000000000 --- a/hoomd/md/ConstraintSphere.h +++ /dev/null @@ -1,65 +0,0 @@ -// Copyright (c) 2009-2021 The Regents of the University of Michigan -// This file is part of the HOOMD-blue project, released under the BSD 3-Clause License. - -// Maintainer: joaander - -#include "hoomd/ForceConstraint.h" -#include "hoomd/ParticleGroup.h" - -#include - -/*! \file ConstraintSphere.h - \brief Declares a class for computing sphere constraint forces -*/ - -#ifdef __HIPCC__ -#error This header cannot be compiled by nvcc -#endif - -#include - -#ifndef __CONSTRAINT_SPHERE_H__ -#define __CONSTRAINT_SPHERE_H__ - -//! Applys a constraint force to keep a group of particles on a sphere -/*! \ingroup computes - */ -class PYBIND11_EXPORT ConstraintSphere : public ForceConstraint - { - public: - //! Constructs the compute - ConstraintSphere(std::shared_ptr sysdef, - std::shared_ptr group, - Scalar3 P, - Scalar r); - - //! Destructor - virtual ~ConstraintSphere(); - - //! Set the force to a new value - void setSphere(Scalar3 P, Scalar r); - - /** ConstraintSphere removes 1 degree of freedom per particle in the group - - @param query The group over which to compute the removed degrees of freedom - */ - virtual Scalar getNDOFRemoved(std::shared_ptr query); - - protected: - std::shared_ptr - m_group; //!< Group of particles on which this constraint is applied - Scalar3 m_P; //!< Position of the sphere - Scalar m_r; //!< Radius of the sphere - - //! Actually compute the forces - virtual void computeForces(uint64_t timestep); - - private: - //! Validate that the sphere is in the box and all particles are very near the constraint - void validate(); - }; - -//! Exports the ConstraintSphere class to python -void export_ConstraintSphere(pybind11::module& m); - -#endif diff --git a/hoomd/md/ConstraintSphereGPU.cc b/hoomd/md/ConstraintSphereGPU.cc deleted file mode 100644 index f3dc80f6d6..0000000000 --- a/hoomd/md/ConstraintSphereGPU.cc +++ /dev/null @@ -1,98 +0,0 @@ -// Copyright (c) 2009-2021 The Regents of the University of Michigan -// This file is part of the HOOMD-blue project, released under the BSD 3-Clause License. - -// Maintainer: joaander - -#include "ConstraintSphereGPU.h" -#include "ConstraintSphereGPU.cuh" - -namespace py = pybind11; - -using namespace std; - -/*! \file ConstraintSphereGPU.cc - \brief Contains code for the ConstraintSphereGPU class -*/ - -/*! \param sysdef SystemDefinition containing the ParticleData to compute forces on - \param group Group of particles on which to apply this constraint - \param P position of the sphere - \param r radius of the sphere -*/ -ConstraintSphereGPU::ConstraintSphereGPU(std::shared_ptr sysdef, - std::shared_ptr group, - Scalar3 P, - Scalar r) - : ConstraintSphere(sysdef, group, P, r), m_block_size(256) - { - if (!m_exec_conf->isCUDAEnabled()) - { - m_exec_conf->msg->error() - << "Creating a ConstraintSphereGPU with no GPU in the execution configuration" << endl; - throw std::runtime_error("Error initializing ConstraintSphereGPU"); - } - } - -/*! Computes the specified constraint forces - \param timestep Current timestep -*/ -void ConstraintSphereGPU::computeForces(uint64_t timestep) - { - unsigned int group_size = m_group->getNumMembers(); - if (group_size == 0) - return; - - if (m_prof) - m_prof->push(m_exec_conf, "ConstraintSphere"); - - assert(m_pdata); - - // access the particle data arrays - const GlobalArray& net_force = m_pdata->getNetForce(); - ArrayHandle d_net_force(net_force, access_location::device, access_mode::read); - - const GlobalArray& group_members = m_group->getIndexArray(); - ArrayHandle d_group_members(group_members, - access_location::device, - access_mode::read); - - ArrayHandle d_pos(m_pdata->getPositions(), access_location::device, access_mode::read); - ArrayHandle d_vel(m_pdata->getVelocities(), - access_location::device, - access_mode::read); - - ArrayHandle d_force(m_force, access_location::device, access_mode::overwrite); - ArrayHandle d_virial(m_virial, access_location::device, access_mode::overwrite); - - // run the kernel in parallel on all GPUs - gpu_compute_constraint_sphere_forces(d_force.data, - d_virial.data, - m_virial.getPitch(), - d_group_members.data, - m_group->getNumMembers(), - m_pdata->getN(), - d_pos.data, - d_vel.data, - d_net_force.data, - m_P, - m_r, - m_deltaT, - m_block_size); - - if (m_exec_conf->isCUDAErrorCheckingEnabled()) - CHECK_CUDA_ERROR(); - - if (m_prof) - m_prof->pop(m_exec_conf); - } - -void export_ConstraintSphereGPU(py::module& m) - { - py::class_>( - m, - "ConstraintSphereGPU") - .def(py::init, - std::shared_ptr, - Scalar3, - Scalar>()); - } diff --git a/hoomd/md/ConstraintSphereGPU.cu b/hoomd/md/ConstraintSphereGPU.cu deleted file mode 100644 index 72271ebafd..0000000000 --- a/hoomd/md/ConstraintSphereGPU.cu +++ /dev/null @@ -1,142 +0,0 @@ -#include "hip/hip_runtime.h" -// Copyright (c) 2009-2021 The Regents of the University of Michigan -// This file is part of the HOOMD-blue project, released under the BSD 3-Clause License. - -// Maintainer: joaander - -#include "ConstraintSphereGPU.cuh" -#include "EvaluatorConstraint.h" -#include "EvaluatorConstraintSphere.h" - -#include - -/*! \file ConstraintSphereGPU.cu - \brief Defines GPU kernel code for calculating sphere constraint forces. Used by - ConstraintSphereGPU. -*/ - -//! Kernel for calculating sphere constraint forces on the GPU -/*! \param d_force Device memory to write computed forces - \param d_virial Device memory to write computed virials - \param virial_pitch pitch of 2D virial array - \param d_group_members List of members in the group - \param group_size number of members in the group - \param N number of particles in system - \param d_pos particle positions on device - \param d_vel particle velocities and masses on device - \param d_net_force Total unconstrained net force on the particles - \param P Position of the sphere - \param r radius of the sphere - \param deltaT step size from the Integrator -*/ -extern "C" __global__ void -gpu_compute_constraint_sphere_forces_kernel(Scalar4* d_force, - Scalar* d_virial, - const size_t virial_pitch, - const unsigned int* d_group_members, - unsigned int group_size, - const unsigned int N, - const Scalar4* d_pos, - const Scalar4* d_vel, - const Scalar4* d_net_force, - Scalar3 P, - Scalar r, - Scalar deltaT) - { - // start by identifying which particle we are to handle - // determine which particle this thread works on - int group_idx = blockIdx.x * blockDim.x + threadIdx.x; - - if (group_idx >= group_size) - return; - - unsigned int idx = d_group_members[group_idx]; - - // read in position, velocity, net force, and mass - Scalar4 pos = d_pos[idx]; - Scalar4 vel = d_vel[idx]; - Scalar4 net_force = d_net_force[idx]; - Scalar m = vel.w; - - // convert to Scalar3's for passing to the evaluators - Scalar3 X = make_scalar3(pos.x, pos.y, pos.z); - Scalar3 V = make_scalar3(vel.x, vel.y, vel.z); - Scalar3 F = make_scalar3(net_force.x, net_force.y, net_force.z); - - // evaluate the constraint position - EvaluatorConstraint constraint(X, V, F, m, deltaT); - EvaluatorConstraintSphere sphere(P, r); - Scalar3 C = sphere.evalClosest(constraint.evalU()); - - // evaluate the constraint force - Scalar3 FC; - Scalar virial[6]; - constraint.evalConstraintForce(FC, virial, C); - - // now that the force calculation is complete, write out the results - d_force[idx] = make_scalar4(FC.x, FC.y, FC.z, Scalar(0.0)); - for (unsigned int i = 0; i < 6; i++) - d_virial[i * virial_pitch + idx] = virial[i]; - } - -/*! \param d_force Device memory to write computed forces - \param d_virial Device memory to write computed virials - \param virial_pitch pitch of 2D virial array - \param d_group_members List of members in the group - \param group_size number of members in the group - \param N number of particles - \param d_pos particle positions on the device - \param d_vel particle velocities on the device - \param d_net_force Total unconstrained net force on the particles - \param P Position of the sphere - \param r radius of the sphere - \param deltaT step size from the Integrator - \param block_size Block size to execute on the GPU - - \returns Any error code resulting from the kernel launch - \note Always returns hipSuccess in release builds to avoid the hipDeviceSynchronize() -*/ -hipError_t gpu_compute_constraint_sphere_forces(Scalar4* d_force, - Scalar* d_virial, - const size_t virial_pitch, - const unsigned int* d_group_members, - unsigned int group_size, - const unsigned int N, - const Scalar4* d_pos, - const Scalar4* d_vel, - const Scalar4* d_net_force, - const Scalar3& P, - Scalar r, - Scalar deltaT, - unsigned int block_size) - { - assert(d_group_members); - assert(d_net_force); - - // setup the grid to run the kernel - dim3 grid(group_size / block_size + 1, 1, 1); - dim3 threads(block_size, 1, 1); - - // run the kernel - hipMemset(d_force, 0, sizeof(Scalar4) * N); - hipMemset(d_virial, 0, 6 * sizeof(Scalar) * virial_pitch); - hipLaunchKernelGGL((gpu_compute_constraint_sphere_forces_kernel), - dim3(grid), - dim3(threads), - 0, - 0, - d_force, - d_virial, - virial_pitch, - d_group_members, - group_size, - N, - d_pos, - d_vel, - d_net_force, - P, - r, - deltaT); - - return hipSuccess; - } diff --git a/hoomd/md/ConstraintSphereGPU.cuh b/hoomd/md/ConstraintSphereGPU.cuh deleted file mode 100644 index e7d4cf4ea2..0000000000 --- a/hoomd/md/ConstraintSphereGPU.cuh +++ /dev/null @@ -1,32 +0,0 @@ -// Copyright (c) 2009-2021 The Regents of the University of Michigan -// This file is part of the HOOMD-blue project, released under the BSD 3-Clause License. - -// Maintainer: joaander - -#include "hoomd/HOOMDMath.h" -#include "hoomd/ParticleData.cuh" - -/*! \file ConstraintSphereGPU.cuh - \brief Declares GPU kernel code for calculating sphere constraint forces on the GPU. Used by - ConstraintSphereGPU. -*/ - -#ifndef __CONSTRAINT_SPHERE_GPU_CUH__ -#define __CONSTRAINT_SPHERE_GPU_CUH__ - -//! Kernel driver that computes harmonic bond forces for HarmonicBondForceComputeGPU -hipError_t gpu_compute_constraint_sphere_forces(Scalar4* d_force, - Scalar* d_virial, - const size_t virial_pitch, - const unsigned int* d_group_members, - unsigned int group_size, - const unsigned int N, - const Scalar4* d_pos, - const Scalar4* d_vel, - const Scalar4* d_net_force, - const Scalar3& P, - Scalar r, - Scalar deltaT, - unsigned int block_size); - -#endif diff --git a/hoomd/md/ConstraintSphereGPU.h b/hoomd/md/ConstraintSphereGPU.h deleted file mode 100644 index 623d5fe35b..0000000000 --- a/hoomd/md/ConstraintSphereGPU.h +++ /dev/null @@ -1,43 +0,0 @@ -// Copyright (c) 2009-2021 The Regents of the University of Michigan -// This file is part of the HOOMD-blue project, released under the BSD 3-Clause License. - -// Maintainer: joaander - -#include "ConstraintSphere.h" - -/*! \file ConstraintSphereGPU.h - \brief Declares a class for computing sphere constraint forces on the GPU -*/ - -#ifdef __HIPCC__ -#error This header cannot be compiled by nvcc -#endif - -#include - -#ifndef __CONSTRAINT_SPHERE_GPU_H__ -#define __CONSTRAINT_SPHERE_GPU_H__ - -//! Applys a constraint force to keep a group of particles on a sphere on the GPU -/*! \ingroup computes - */ -class PYBIND11_EXPORT ConstraintSphereGPU : public ConstraintSphere - { - public: - //! Constructs the compute - ConstraintSphereGPU(std::shared_ptr sysdef, - std::shared_ptr group, - Scalar3 P, - Scalar r); - - protected: - unsigned int m_block_size; //!< block size to execute on the GPU - - //! Actually compute the forces - virtual void computeForces(uint64_t timestep); - }; - -//! Exports the ConstraintSphereGPU class to python -void export_ConstraintSphereGPU(pybind11::module& m); - -#endif diff --git a/hoomd/md/EvaluatorConstraint.h b/hoomd/md/EvaluatorConstraint.h deleted file mode 100644 index 0b64e14907..0000000000 --- a/hoomd/md/EvaluatorConstraint.h +++ /dev/null @@ -1,89 +0,0 @@ -// Copyright (c) 2009-2021 The Regents of the University of Michigan -// This file is part of the HOOMD-blue project, released under the BSD 3-Clause License. - -// Maintainer: joaander - -#ifndef __EVALUATOR_CONSTRAINT_H__ -#define __EVALUATOR_CONSTRAINT_H__ - -#include "hoomd/HOOMDMath.h" - -/*! \file EvaluatorConstraint.h - \brief Defines basic evaluation methods common to all constraint force implementations -*/ - -// need to declare these class methods with __device__ qualifiers when building in nvcc -// DEVICE is __host__ __device__ when included in nvcc and blank when included into the host -// compiler -#ifdef __HIPCC__ -#define DEVICE __device__ -#else -#define DEVICE -#endif - -//! Class for evaluating constraint forces -/*! General Overview - EvaluatorConstraint is a low level computation class for use on the CPU and GPU. It provides - basic functionality needed by all constraint forces. -*/ -class EvaluatorConstraint - { - public: - //! Constructs the constraint evaluator - /*! \param _X Current position at the time of the constraint force calculation - \param V Current velocity at the time of the constraint force calculation - \param F Current net force at the time of the constraint force calculation - \param _m Mass of the particle - \param _deltaT Step size delta t - */ - DEVICE EvaluatorConstraint(Scalar3 _X, Scalar3 V, Scalar3 F, Scalar _m, Scalar _deltaT) - : X(_X), m(_m), deltaT(_deltaT) - { - // perform step 2 of this velocity verlet update and step 1 of the next to get - // U = X(t+2deltaT) given X = X(t+deltaT) - Scalar minv = Scalar(1.0) / m; - Scalar dtsqdivm = deltaT * deltaT * minv; - U.x = X.x + V.x * deltaT + F.x * dtsqdivm; - U.y = X.y + V.y * deltaT + F.y * dtsqdivm; - U.z = X.z + V.z * deltaT + F.z * dtsqdivm; - } - - //! Evaluate the unconstrained position update U - /*! \returns The unconstrained position update U - */ - DEVICE Scalar3 evalU() - { - return U; - } - - //! Evaluate the additional constraint force - /*! \param FC output parameter where the computed force is written - \param virial array of six scalars the computed virial tensor is written - \param C constrained position particle will be moved to at the next step - \return Additional force \a F needed to satisfy the constraint - */ - DEVICE void evalConstraintForce(Scalar3& FC, Scalar* virial, const Scalar3& C) - { - // subtract a constrained update from U and get that F = (C-U)*m/dt^2 - Scalar moverdtsq = m / (deltaT * deltaT); - FC.x = (C.x - U.x) * moverdtsq; - FC.y = (C.y - U.y) * moverdtsq; - FC.z = (C.z - U.z) * moverdtsq; - - // compute virial - virial[0] = FC.x * X.x; - virial[1] = Scalar(1. / 2.) * (FC.y * X.x + FC.x * X.y); - virial[2] = Scalar(1. / 2.) * (FC.z * X.x + FC.x * X.z); - virial[3] = FC.y * X.y; - virial[4] = Scalar(1. / 2.) * (FC.z * X.y + FC.y * X.z); - virial[5] = FC.z * X.z; - } - - protected: - Scalar3 U; //!< Unconstrained position update - Scalar3 X; //!< Current particle position - Scalar m; //!< Saved mass value - Scalar deltaT; //!< Saved delta T value - }; - -#endif // __PAIR_EVALUATOR_LJ_H__ diff --git a/hoomd/md/EvaluatorConstraintEllipsoid.h b/hoomd/md/EvaluatorConstraintEllipsoid.h deleted file mode 100644 index bf492c9611..0000000000 --- a/hoomd/md/EvaluatorConstraintEllipsoid.h +++ /dev/null @@ -1,206 +0,0 @@ -// Copyright (c) 2009-2021 The Regents of the University of Michigan -// This file is part of the HOOMD-blue project, released under the BSD 3-Clause License. - -// Maintainer: joaander - -#ifndef __EVALUATOR_CONSTRAINT_Ellipsoid_H__ -#define __EVALUATOR_CONSTRAINT_Ellipsoid_H__ - -#include "hoomd/HOOMDMath.h" -using namespace std; - -/*! \file EvaluatorConstraintEllipsoid.h - \brief Defines the constraint evaluator class for ellipsoids -*/ - -// need to declare these class methods with __device__ qualifiers when building in nvcc -// DEVICE is __host__ __device__ when included in nvcc and blank when included into the host -// compiler -#ifdef __HIPCC__ -#define DEVICE __device__ -#else -#define DEVICE -#endif - -//! Class for evaluating ellipsoid constraints -/*! General Overview - EvaluatorConstraintEllipsoid is a low level computation helper class to aid in evaluating - particle constraints on a ellipsoid. Given a ellipsoid at a given position and radii, it will - find the nearest point on the Ellipsoid to a given position. -*/ -class EvaluatorConstraintEllipsoid - { - public: - //! Constructs the constraint evaluator - /*! \param _P Position of the ellipsoid - \param _rx Radius of the ellipsoid in the X direction - \param _ry Radius of the ellipsoid in the Y direction - \param _rz Radius of the ellipsoid in the Z direction - - NOTE: For the algorithm to work, we must have _rx >= _rz, ry >= _rz, and _rz > 0. - */ - DEVICE EvaluatorConstraintEllipsoid(Scalar3 _P, Scalar _rx, Scalar _ry, Scalar _rz) - : P(_P), rx(_rx), ry(_ry), rz(_rz) - { - } - - //! Evaluate the closest point on the ellipsoid. Method from: - //! http://www.geometrictools.com/Documentation/DistancePointEllipseEllipsoid.pdf - /*! \param U unconstrained point - - \return Nearest point on the ellipsoid - */ - DEVICE Scalar3 evalClosest(const Scalar3& U) - { - if (rx == ry && ry == rz) // if ellipsoid is actually a sphere, use easier method - { - // compute the vector pointing from P to V - Scalar3 V; - V.x = U.x - P.x; - V.y = U.y - P.y; - V.z = U.z - P.z; - - // compute 1/magnitude of V - Scalar magVinv = slow::rsqrt(V.x * V.x + V.y * V.y + V.z * V.z); - - // compute Vhat, the unit vector pointing in the direction of V - Scalar3 Vhat; - Vhat.x = magVinv * V.x; - Vhat.y = magVinv * V.y; - Vhat.z = magVinv * V.z; - - // compute resulting constrained point - Scalar3 C; - C.x = P.x + Vhat.x * rx; - C.y = P.y + Vhat.y * rx; - C.z = P.z + Vhat.z * rx; - - return C; - } - else // else use iterative method - { - Scalar xsign, ysign, zsign; // sign of point's position - Scalar y0, y1, y2; - if (U.x < 0) - { - xsign = -1; - } - else - { - xsign = 1; - } - if (U.y < 0) - { - ysign = -1; - } - else - { - ysign = 1; - } - if (U.z < 0) - { - zsign = -1; - } - else - { - zsign = 1; - } - y0 = U.x * xsign; - y1 = U.y * ysign; - y2 = U.z * zsign; - - Scalar z0 = y0 / rx; - Scalar z1 = y1 / ry; - Scalar z2 = y2 / rz; - Scalar g = z0 * z0 + z1 * z1 + z2 * z2 - 1; - if (g != 0) // point does not lay on ellipsoid - { - Scalar r0 = (rx / rz) * (rx / rz); - Scalar r1 = (ry / rz) * (ry / rz); - - // iterative method of calculus to find closest point - Scalar n0 = r0 * z0; - Scalar n1 = r1 * z1; - Scalar s0 = z1 - 1; - Scalar s1 = 0; - if (g > 0) - { - s1 = slow::sqrt(n0 * n0 + n1 * n1 + z2 * z2) - 1; - } - Scalar sbar; - Scalar ratio0, ratio1, ratio2; - int i, imax = 10000; // When tested, on average takes ~30 steps to complete, and - // never more than 150. - for (i = 0; i < imax; i++) - { - sbar = (s0 + s1) / 2.0; - if (sbar == s0 || sbar == s1) - { - break; - } - ratio0 = n0 / (sbar + r0); - ratio1 = n1 / (sbar + r1); - ratio2 = z2 / (sbar + 1); - g = ratio0 * ratio0 + ratio1 * ratio1 + ratio2 * ratio2 - 1; - if (g > 0) - { - s0 = sbar; - } - else if (g < 0) - { - s1 = sbar; - } - else - { - break; - } - } - // if (i == imax) - // { - // throw runtime_error("constrain.ellipsoid: Not enough iteration steps to find - // closest point on ellipsoid.\n"); - // } - - // compute resulting constrained point - Scalar3 C; - C.x = xsign * r0 * y0 / (sbar + r0); - C.y = ysign * r1 * y1 / (sbar + r1); - C.z = zsign * y2 / (sbar + 1); - - return C; - } - else // trivial case of point laying on ellipsoid - { - return U; - } - } - } - - //! Evaluate the normal unit vector for point on the ellipsoid. - /*! \param U point on ellipsoid - \return normal unit vector for point on the ellipsoid - */ - DEVICE Scalar3 evalNormal(const Scalar3& U) - { - Scalar3 N; - N.x = U.x / (rx * rx); - N.y = U.y / (ry * ry); - N.z = U.z / (rz * rz); - - Scalar nNorm; - nNorm = slow::sqrt(N.x * N.x + N.y * N.y + N.z * N.z); - N.x /= nNorm; - N.y /= nNorm; - N.z /= nNorm; - - return N; - } - - protected: - Scalar3 P; //!< Position of the ellipsoid - Scalar rx; //!< radius of the ellipsoid in the X direction - Scalar ry; //!< radius of the ellipsoid in the Y direction - Scalar rz; //!< radius of the ellipsoid in the Z direction - }; - -#endif // __PAIR_EVALUATOR_LJ_H__ diff --git a/hoomd/md/EvaluatorConstraintSphere.h b/hoomd/md/EvaluatorConstraintSphere.h deleted file mode 100644 index e3972e710b..0000000000 --- a/hoomd/md/EvaluatorConstraintSphere.h +++ /dev/null @@ -1,75 +0,0 @@ -// Copyright (c) 2009-2021 The Regents of the University of Michigan -// This file is part of the HOOMD-blue project, released under the BSD 3-Clause License. - -// Maintainer: joaander - -#ifndef __EVALUATOR_CONSTRAINT_SPHERE_H__ -#define __EVALUATOR_CONSTRAINT_SPHERE_H__ - -#include "hoomd/HOOMDMath.h" - -/*! \file EvaluatorConstraintSphere.h - \brief Defines the constraint evaluator class for spheres -*/ - -// need to declare these class methods with __device__ qualifiers when building in nvcc -// DEVICE is __host__ __device__ when included in nvcc and blank when included into the host -// compiler -#ifdef __HIPCC__ -#define DEVICE __device__ -#else -#define DEVICE -#endif - -//! Class for evaluating sphere constraints -/*! General Overview - EvaluatorConstraintSphere is a low level computation helper class to aid in evaluating particle - constraints on a sphere. Given a sphere at a given position and radius, it will find the nearest - point on the sphere to a given point. -*/ -class EvaluatorConstraintSphere - { - public: - //! Constructs the constraint evaluator - /*! \param _P Position of the sphere - \param _r Radius of the sphere - */ - DEVICE EvaluatorConstraintSphere(Scalar3 _P, Scalar _r) : P(_P), r(_r) { } - - //! Evaluate the closest point on the sphere - /*! \param U unconstrained point - - \return Nearest point on the sphere - */ - DEVICE Scalar3 evalClosest(const Scalar3& U) - { - // compute the vector pointing from P to V - Scalar3 V; - V.x = U.x - P.x; - V.y = U.y - P.y; - V.z = U.z - P.z; - - // compute 1/magnitude of V - Scalar magVinv = fast::rsqrt(V.x * V.x + V.y * V.y + V.z * V.z); - - // compute Vhat, the unit vector pointing in the direction of V - Scalar3 Vhat; - Vhat.x = magVinv * V.x; - Vhat.y = magVinv * V.y; - Vhat.z = magVinv * V.z; - - // compute resulting constrained point - Scalar3 C; - C.x = P.x + Vhat.x * r; - C.y = P.y + Vhat.y * r; - C.z = P.z + Vhat.z * r; - - return C; - } - - protected: - Scalar3 P; //!< Position of the sphere - Scalar r; //!< radius of the sphere - }; - -#endif // __PAIR_EVALUATOR_LJ_H__ diff --git a/hoomd/md/NeighborListGPU.cc b/hoomd/md/NeighborListGPU.cc index d1875c18cb..2c96687a76 100644 --- a/hoomd/md/NeighborListGPU.cc +++ b/hoomd/md/NeighborListGPU.cc @@ -185,7 +185,7 @@ void NeighborListGPU::filterNlist() //! Update the exclusion list on the GPU void NeighborListGPU::updateExListIdx() { - assert(!m_need_reallocate_exlist); + assert(!m_n_particles_changed); if (m_prof) m_prof->push(m_exec_conf, "update-ex"); diff --git a/hoomd/md/OneDConstraint.cc b/hoomd/md/OneDConstraint.cc deleted file mode 100644 index 7a2a4a4ae8..0000000000 --- a/hoomd/md/OneDConstraint.cc +++ /dev/null @@ -1,122 +0,0 @@ -// Copyright (c) 2009-2021 The Regents of the University of Michigan -// This file is part of the HOOMD-blue project, released under the BSD 3-Clause License. - -// Maintainer: joaander - -#include "OneDConstraint.h" -#include "EvaluatorConstraint.h" - -namespace py = pybind11; - -using namespace std; - -/*! \file OneDConstraint.cc - \brief Contains code for the OneDConstraint class -*/ - -/*! \param sysdef SystemDefinition containing the ParticleData to compute forces on - \param group Group of particles on which to apply this constraint -*/ -OneDConstraint::OneDConstraint(std::shared_ptr sysdef, - std::shared_ptr group, - Scalar3 constraint_vec) - : ForceConstraint(sysdef), m_group(group), m_vec(constraint_vec) - { - m_exec_conf->msg->notice(5) << "Constructing OneDConstraint" << endl; - } - -OneDConstraint::~OneDConstraint() - { - m_exec_conf->msg->notice(5) << "Destroying OneDConstraint" << endl; - } - -/*! - \param constraint_vec direction that particles are constrained to -*/ -void OneDConstraint::setVector(Scalar3 constraint_vec) - { - m_vec = constraint_vec; - } - -/*! Computes the specified constraint forces - \param timestep Current timestep -*/ -void OneDConstraint::computeForces(uint64_t timestep) - { - unsigned int group_size = m_group->getNumMembers(); - if (group_size == 0) - return; - - if (m_prof) - m_prof->push("OneDConstraint"); - - assert(m_pdata); - // access the particle data arrays - ArrayHandle h_pos(m_pdata->getPositions(), access_location::host, access_mode::read); - ArrayHandle h_vel(m_pdata->getVelocities(), access_location::host, access_mode::read); - - const GlobalArray& net_force = m_pdata->getNetForce(); - ArrayHandle h_net_force(net_force, access_location::host, access_mode::read); - - ArrayHandle h_force(m_force, access_location::host, access_mode::overwrite); - ArrayHandle h_virial(m_virial, access_location::host, access_mode::overwrite); - size_t virial_pitch = m_virial.getPitch(); - - // Zero data for force calculation. - memset((void*)h_force.data, 0, sizeof(Scalar4) * m_force.getNumElements()); - memset((void*)h_virial.data, 0, sizeof(Scalar) * m_virial.getNumElements()); - - // there are enough other checks on the input data: but it doesn't hurt to be safe - assert(h_force.data); - assert(h_virial.data); - - // for each of the particles in the group - for (unsigned int group_idx = 0; group_idx < group_size; group_idx++) - { - // get the current particle properties - unsigned int j = m_group->getMemberIndex(group_idx); - Scalar3 X = make_scalar3(h_pos.data[j].x, h_pos.data[j].y, h_pos.data[j].z); - Scalar3 V = make_scalar3(h_vel.data[j].x, h_vel.data[j].y, h_vel.data[j].z); - Scalar3 F - = make_scalar3(h_net_force.data[j].x, h_net_force.data[j].y, h_net_force.data[j].z); - Scalar m = h_vel.data[j].w; - - // evaluate the constraint position - EvaluatorConstraint constraint(X, V, F, m, m_deltaT); - Scalar3 U = constraint.evalU(); - Scalar3 D = make_scalar3((U.x - X.x), (U.y - X.y), (U.z - X.z)); - Scalar n = (D.x * m_vec.x + D.y * m_vec.y + D.z * m_vec.z) - / (m_vec.x * m_vec.x + m_vec.y * m_vec.y + m_vec.z * m_vec.z); - Scalar3 C = make_scalar3((n * m_vec.x + X.x), (n * m_vec.y + X.y), (n * m_vec.z + X.z)); - - // evaluate the constraint force - Scalar3 FC; - Scalar virial[6]; - constraint.evalConstraintForce(FC, virial, C); - - // apply the constraint force - h_force.data[j].x = FC.x; - h_force.data[j].y = FC.y; - h_force.data[j].z = FC.z; - for (int k = 0; k < 6; k++) - h_virial.data[k * virial_pitch + j] = virial[k]; - } - - if (m_prof) - m_prof->pop(); - } - -Scalar OneDConstraint::getNDOFRemoved(std::shared_ptr query) - { - // OneDConstraint removes 2 degrees of freedom per particle in the group - return m_group->intersectionSize(query) * 2; - } - -void export_OneDConstraint(py::module& m) - { - py::class_>(m, - "OneDConstraint") - .def(py::init, std::shared_ptr, Scalar3>()) - .def("getNDOFRemoved", &OneDConstraint::getNDOFRemoved) - .def("setVector", &OneDConstraint::setVector); - } diff --git a/hoomd/md/OneDConstraint.h b/hoomd/md/OneDConstraint.h deleted file mode 100644 index 627c5ec31d..0000000000 --- a/hoomd/md/OneDConstraint.h +++ /dev/null @@ -1,56 +0,0 @@ -// Copyright (c) 2009-2021 The Regents of the University of Michigan -// This file is part of the HOOMD-blue project, released under the BSD 3-Clause License. - -// Maintainer: joaander - -#include "hoomd/ForceConstraint.h" -#include "hoomd/ParticleGroup.h" - -#include - -/*! \file OneDConstraint.h - \brief Declares a class for computing 1D constraint forces -*/ - -#ifdef __HIPCC__ -#error This header cannot be compiled by nvcc -#endif - -#include - -#ifndef __ONE_D_CONSTRAINT_H__ -#define __ONE_D_CONSTRAINT_H__ - -//! Applys a constraint force to prevent motion in x and y directions -/*! \ingroup computes - */ -class PYBIND11_EXPORT OneDConstraint : public ForceConstraint - { - public: - //! Constructs the compute - OneDConstraint(std::shared_ptr sysdef, - std::shared_ptr group, - Scalar3 constraint_vec); - - //! Destructor - virtual ~OneDConstraint(); - - //! Set the force to a new value - void setVector(Scalar3 constraint_vec); - - //! Return the number of DOF removed by this constraint - virtual Scalar getNDOFRemoved(std::shared_ptr query); - - protected: - std::shared_ptr - m_group; //!< Group of particles on which this constraint is applied - Scalar3 m_vec; //!< The vector along which particles are constrained - - //! Actually compute the forces - virtual void computeForces(uint64_t timestep); - }; - -//! Exports the OneDConstraint class to python -void export_OneDConstraint(pybind11::module& m); - -#endif diff --git a/hoomd/md/OneDConstraintGPU.cc b/hoomd/md/OneDConstraintGPU.cc deleted file mode 100644 index 2c1d63019b..0000000000 --- a/hoomd/md/OneDConstraintGPU.cc +++ /dev/null @@ -1,103 +0,0 @@ -// Copyright (c) 2009-2021 The Regents of the University of Michigan -// This file is part of the HOOMD-blue project, released under the BSD 3-Clause License. - -// Maintainer: joaander - -#include "OneDConstraintGPU.h" -#include "OneDConstraintGPU.cuh" - -namespace py = pybind11; - -using namespace std; - -/*! \file OneDConstraintGPU.cc - \brief Contains code for the OneDConstraintGPU class -*/ - -/*! \param sysdef SystemDefinition containing the ParticleData to compute forces on - \param group Group of particles on which to apply this constraint - \param P position of the sphere - \param r radius of the sphere -*/ -OneDConstraintGPU::OneDConstraintGPU(std::shared_ptr sysdef, - std::shared_ptr group, - Scalar3 constraint_vec) - : OneDConstraint(sysdef, group, constraint_vec), m_block_size(256) - { - if (!m_exec_conf->isCUDAEnabled()) - { - m_exec_conf->msg->error() - << "Creating a OneDConstraintGPU with no GPU in the execution configuration" << endl; - throw std::runtime_error("Error initializing OneDConstraintGPU"); - } - - unsigned int warp_size = m_exec_conf->dev_prop.warpSize; - m_tuner.reset( - new Autotuner(warp_size, 1024, warp_size, 5, 100000, "oneD_constraint", this->m_exec_conf)); - } - -/*! Computes the specified constraint forces - \param timestep Current timestep -*/ -void OneDConstraintGPU::computeForces(uint64_t timestep) - { - unsigned int group_size = m_group->getNumMembers(); - if (group_size == 0) - return; - - if (m_prof) - m_prof->push(m_exec_conf, "OneDConstraint"); - - assert(m_pdata); - - // access the particle data arrays - const GlobalArray& net_force = m_pdata->getNetForce(); - ArrayHandle d_net_force(net_force, access_location::device, access_mode::read); - - const GlobalArray& group_members = m_group->getIndexArray(); - ArrayHandle d_group_members(group_members, - access_location::device, - access_mode::read); - - ArrayHandle d_pos(m_pdata->getPositions(), access_location::device, access_mode::read); - ArrayHandle d_vel(m_pdata->getVelocities(), - access_location::device, - access_mode::read); - - ArrayHandle d_force(m_force, access_location::device, access_mode::overwrite); - ArrayHandle d_virial(m_virial, access_location::device, access_mode::overwrite); - - // run the kernel in parallel on all GPUs - m_tuner->begin(); - gpu_compute_one_d_constraint_forces(d_force.data, - d_virial.data, - m_virial.getPitch(), - d_group_members.data, - m_group->getNumMembers(), - m_pdata->getN(), - d_pos.data, - d_vel.data, - d_net_force.data, - m_deltaT, - m_tuner->getParam(), - m_vec); - - if (m_exec_conf->isCUDAErrorCheckingEnabled()) - CHECK_CUDA_ERROR(); - m_tuner->end(); - - if (m_exec_conf->isCUDAErrorCheckingEnabled()) - CHECK_CUDA_ERROR(); - - if (m_prof) - m_prof->pop(m_exec_conf); - } - -void export_OneDConstraintGPU(py::module& m) - { - py::class_>( - m, - "OneDConstraintGPU") - .def( - py::init, std::shared_ptr, Scalar3>()); - } diff --git a/hoomd/md/OneDConstraintGPU.cu b/hoomd/md/OneDConstraintGPU.cu deleted file mode 100644 index 79f3b67e16..0000000000 --- a/hoomd/md/OneDConstraintGPU.cu +++ /dev/null @@ -1,138 +0,0 @@ -#include "hip/hip_runtime.h" -// Copyright (c) 2009-2021 The Regents of the University of Michigan -// This file is part of the HOOMD-blue project, released under the BSD 3-Clause License. - -// Maintainer: joaander - -#include "EvaluatorConstraint.h" -#include "OneDConstraintGPU.cuh" - -#include - -/*! \file OneDConstraintGPU.cu - \brief Defines GPU kernel code for calculating one dimensional constraint forces. Used by - OneDConstraintGPU. -*/ - -//! Kernel for calculating one dimensional constraint forces on the GPU -/*! \param d_force Device memory to write computed forces - \param d_virial Device memory to write computed virials - \param virial_pitch pitch of 2D virial array - \param d_group_members List of members in the group - \param group_size number of members in the group - \param N number of particles in system - \param d_pos particle positions on device - \param d_vel particle velocities and masses on device - \param d_net_force Total unconstrained net force on the particles - \param deltaT step size from the Integrator -*/ -extern "C" __global__ void -gpu_compute_one_d_constraint_forces_kernel(Scalar4* d_force, - Scalar* d_virial, - const size_t virial_pitch, - const unsigned int* d_group_members, - unsigned int group_size, - const unsigned int N, - const Scalar4* d_pos, - const Scalar4* d_vel, - const Scalar4* d_net_force, - Scalar deltaT, - Scalar3 m_vec) - { - // start by identifying which particle we are to handle - // determine which particle this thread works on - int group_idx = blockIdx.x * blockDim.x + threadIdx.x; - - if (group_idx >= group_size) - return; - - unsigned int idx = d_group_members[group_idx]; - - // read in position, velocity, net force, and mass - Scalar4 pos = d_pos[idx]; - Scalar4 vel = d_vel[idx]; - Scalar4 net_force = d_net_force[idx]; - Scalar m = vel.w; - - // convert to Scalar3's for passing to the evaluators - Scalar3 X = make_scalar3(pos.x, pos.y, pos.z); - Scalar3 V = make_scalar3(vel.x, vel.y, vel.z); - Scalar3 F = make_scalar3(net_force.x, net_force.y, net_force.z); - - // evaluate the constraint position - EvaluatorConstraint constraint(X, V, F, m, deltaT); - - // evaluate the constraint force - Scalar3 FC; - Scalar virial[6]; - Scalar3 U = constraint.evalU(); - Scalar3 D = make_scalar3((U.x - X.x), (U.y - X.y), (U.z - X.z)); - Scalar n = (D.x * m_vec.x + D.y * m_vec.y + D.z * m_vec.z) - / (m_vec.x * m_vec.x + m_vec.y * m_vec.y + m_vec.z * m_vec.z); - Scalar3 C = make_scalar3((n * m_vec.x + X.x), (n * m_vec.y + X.y), (n * m_vec.z + X.z)); - - constraint.evalConstraintForce(FC, virial, C); - - // now that the force calculation is complete, write out the results - d_force[idx] = make_scalar4(FC.x, FC.y, FC.z, Scalar(0.0)); - for (unsigned int i = 0; i < 6; i++) - d_virial[i * virial_pitch + idx] = virial[i]; - } - -/*! \param d_force Device memory to write computed forces - \param d_virial Device memory to write computed virials - \param virial_pitch pitch of 2D virial array - \param d_group_members List of members in the group - \param group_size number of members in the group - \param N number of particles - \param d_pos particle positions on the device - \param d_vel particle velocities on the device - \param d_net_force Total unconstrained net force on the particles - \param deltaT step size from the Integrator - \param block_size Block size to execute on the GPU - - \returns Any error code resulting from the kernel launch - \note Always returns hipSuccess in release builds to avoid the hipDeviceSynchronize() -*/ -hipError_t gpu_compute_one_d_constraint_forces(Scalar4* d_force, - Scalar* d_virial, - const size_t virial_pitch, - const unsigned int* d_group_members, - unsigned int group_size, - const unsigned int N, - const Scalar4* d_pos, - const Scalar4* d_vel, - const Scalar4* d_net_force, - Scalar deltaT, - unsigned int block_size, - Scalar3 m_vec) - { - assert(d_group_members); - assert(d_net_force); - - // setup the grid to run the kernel - dim3 grid(group_size / block_size + 1, 1, 1); - dim3 threads(block_size, 1, 1); - - // run the kernel - hipMemset(d_force, 0, sizeof(Scalar4) * N); - hipMemset(d_virial, 0, 6 * sizeof(Scalar) * virial_pitch); - hipLaunchKernelGGL((gpu_compute_one_d_constraint_forces_kernel), - dim3(grid), - dim3(threads), - 0, - 0, - d_force, - d_virial, - virial_pitch, - d_group_members, - group_size, - N, - d_pos, - d_vel, - d_net_force, - deltaT, - m_vec); - - return hipSuccess; - } diff --git a/hoomd/md/OneDConstraintGPU.cuh b/hoomd/md/OneDConstraintGPU.cuh deleted file mode 100644 index 82b51715d5..0000000000 --- a/hoomd/md/OneDConstraintGPU.cuh +++ /dev/null @@ -1,31 +0,0 @@ -// Copyright (c) 2009-2021 The Regents of the University of Michigan -// This file is part of the HOOMD-blue project, released under the BSD 3-Clause License. - -// Maintainer: joaander - -#include "hoomd/HOOMDMath.h" -#include "hoomd/ParticleData.cuh" - -/*! \file OneDConstraintGPU.cuh - \brief Declares GPU kernel code for calculating one dimensional constraint forces on the GPU. - Used by OneDConstraintGPU. -*/ - -#ifndef __CONSTRAINT_SPHERE_GPU_CUH__ -#define __CONSTRAINT_SPHERE_GPU_CUH__ - -//! Kernel driver that computes harmonic bond forces for HarmonicBondForceComputeGPU -hipError_t gpu_compute_one_d_constraint_forces(Scalar4* d_force, - Scalar* d_virial, - const size_t virial_pitch, - const unsigned int* d_group_members, - unsigned int group_size, - const unsigned int N, - const Scalar4* d_pos, - const Scalar4* d_vel, - const Scalar4* d_net_force, - Scalar deltaT, - unsigned int block_size, - Scalar3 m_vec); - -#endif diff --git a/hoomd/md/OneDConstraintGPU.h b/hoomd/md/OneDConstraintGPU.h deleted file mode 100644 index 85d3abcd16..0000000000 --- a/hoomd/md/OneDConstraintGPU.h +++ /dev/null @@ -1,56 +0,0 @@ -// Copyright (c) 2009-2021 The Regents of the University of Michigan -// This file is part of the HOOMD-blue project, released under the BSD 3-Clause License. - -// Maintainer: joaander - -#include "OneDConstraint.h" -#include "hoomd/Autotuner.h" - -/*! \file OneDConstraintGPU.h - \brief Declares a class for computing sphere constraint forces on the GPU -*/ - -#ifdef __HIPCC__ -#error This header cannot be compiled by nvcc -#endif - -#include - -#ifndef __ONE_D_CONSTRAINT_GPU_H__ -#define __ONE_D_CONSTRAINT_GPU_H__ - -//! Applys a constraint force to keep a group of particles on a sphere on the GPU -/*! \ingroup computes - */ -class PYBIND11_EXPORT OneDConstraintGPU : public OneDConstraint - { - public: - //! Constructs the compute - OneDConstraintGPU(std::shared_ptr sysdef, - std::shared_ptr group, - Scalar3 constraint_vec); - - //! Set autotuner parameters - /*! \param enable Enable/disable autotuning - \param period period (approximate) in time steps when returning occurs - */ - virtual void setAutotunerParams(bool enable, unsigned int period) - { - OneDConstraint::setAutotunerParams(enable, period); - m_tuner->setPeriod(period); - m_tuner->setEnabled(enable); - } - - protected: - unsigned int m_block_size; //!< block size to execute on the GPU - - std::unique_ptr m_tuner; //!< Autotuner for block size - - //! Actually compute the forces - virtual void computeForces(uint64_t timestep); - }; - -//! Exports the OneDConstraintGPU class to python -void export_OneDConstraintGPU(pybind11::module& m); - -#endif diff --git a/hoomd/md/TwoStepRATTLEBD.h b/hoomd/md/TwoStepRATTLEBD.h index 1b02cd8088..cd8b0d0d31 100644 --- a/hoomd/md/TwoStepRATTLEBD.h +++ b/hoomd/md/TwoStepRATTLEBD.h @@ -44,10 +44,7 @@ template class PYBIND11_EXPORT TwoStepRATTLEBD : public TwoStepL std::shared_ptr T, Scalar tolerance = 0.000001); - virtual ~TwoStepRATTLEBD() - { - m_exec_conf->msg->notice(5) << "Destroying TwoStepRATTLEBD" << endl; - } + virtual ~TwoStepRATTLEBD(); //! Performs the first step of the integration virtual void integrateStepOne(uint64_t timestep); @@ -81,11 +78,18 @@ template class PYBIND11_EXPORT TwoStepRATTLEBD : public TwoStepL }; protected: + //! Helper function to be called when box changes + void setBoxChange() + { + m_box_changed = true; + } + Manifold m_manifold; //!< The manifold used for the RATTLE constraint bool m_noiseless_t; bool m_noiseless_r; Scalar m_tolerance; //!< The tolerance value of the RATTLE algorithm, setting the tolerance to //!< the manifold + bool m_box_changed; }; /*! \file TwoStepRATTLEBD.h @@ -109,16 +113,28 @@ TwoStepRATTLEBD::TwoStepRATTLEBD(std::shared_ptr sys std::shared_ptr T, Scalar tolerance) : TwoStepLangevinBase(sysdef, group, T), m_manifold(manifold), m_noiseless_t(false), - m_noiseless_r(false), m_tolerance(tolerance) + m_noiseless_r(false), m_tolerance(tolerance), m_box_changed(false) { m_exec_conf->msg->notice(5) << "Constructing TwoStepRATTLEBD" << endl; + m_pdata->getBoxChangeSignal() + .template connect, &TwoStepRATTLEBD::setBoxChange>( + this); + if (!m_manifold.fitsInsideBox(m_pdata->getGlobalBox())) { throw std::runtime_error("Parts of the manifold are outside the box"); } } +template TwoStepRATTLEBD::~TwoStepRATTLEBD() + { + m_pdata->getBoxChangeSignal() + .template disconnect, &TwoStepRATTLEBD::setBoxChange>( + this); + m_exec_conf->msg->notice(5) << "Destroying TwoStepRATTLEBD" << endl; + } + /*! \param timestep Current time step \post Particle positions are moved forward to timestep+1 @@ -169,9 +185,13 @@ template void TwoStepRATTLEBD::integrateStepOne(uint64 const BoxDim& box = m_pdata->getBox(); - if (!m_manifold.fitsInsideBox(m_pdata->getGlobalBox())) + if (m_box_changed) { - throw std::runtime_error("Parts of the manifold are outside the box"); + if (!m_manifold.fitsInsideBox(m_pdata->getGlobalBox())) + { + throw std::runtime_error("Parts of the manifold are outside the box"); + } + m_box_changed = false; } uint16_t seed = m_sysdef->getSeed(); diff --git a/hoomd/md/TwoStepRATTLEBDGPU.h b/hoomd/md/TwoStepRATTLEBDGPU.h index d238d049c7..e839777fde 100644 --- a/hoomd/md/TwoStepRATTLEBDGPU.h +++ b/hoomd/md/TwoStepRATTLEBDGPU.h @@ -83,9 +83,13 @@ template void TwoStepRATTLEBDGPU::integrateStepOne(uin if (this->m_prof) this->m_prof->push(this->m_exec_conf, "BD step 1"); - if (!this->m_manifold.fitsInsideBox(this->m_pdata->getGlobalBox())) + if (this->m_box_changed) { - throw std::runtime_error("Parts of the manifold are outside the box"); + if (!this->m_manifold.fitsInsideBox(this->m_pdata->getGlobalBox())) + { + throw std::runtime_error("Parts of the manifold are outside the box"); + } + this->m_box_changed = false; } // access all the needed data diff --git a/hoomd/md/TwoStepRATTLELangevin.h b/hoomd/md/TwoStepRATTLELangevin.h index ff9be142ea..8482a8609e 100644 --- a/hoomd/md/TwoStepRATTLELangevin.h +++ b/hoomd/md/TwoStepRATTLELangevin.h @@ -44,10 +44,7 @@ template class PYBIND11_EXPORT TwoStepRATTLELangevin : public Tw Manifold manifold, std::shared_ptr T, Scalar tolerance = 0.000001); - virtual ~TwoStepRATTLELangevin() - { - m_exec_conf->msg->notice(5) << "Destroying TwoStepRATTLELangevin" << endl; - } + virtual ~TwoStepRATTLELangevin(); void setTallyReservoirEnergy(bool tally) { @@ -98,6 +95,12 @@ template class PYBIND11_EXPORT TwoStepRATTLELangevin : public Tw }; protected: + //! Helper function to be called when box changes + void setBoxChange() + { + m_box_changed = true; + } + Manifold m_manifold; //!< The manifold used for the RATTLE constraint Scalar m_reservoir_energy; //!< The energy of the reservoir the system is coupled to. Scalar @@ -107,6 +110,7 @@ template class PYBIND11_EXPORT TwoStepRATTLELangevin : public Tw bool m_noiseless_r; //!< If set true, there will be no rotational noise (random torque) Scalar m_tolerance; //!< The tolerance value of the RATTLE algorithm, setting the tolerance to //!< the manifold + bool m_box_changed; }; /*! \param sysdef SystemDefinition this method will act on. Must not be NULL. @@ -128,16 +132,28 @@ TwoStepRATTLELangevin::TwoStepRATTLELangevin(std::shared_ptrmsg->notice(5) << "Constructing TwoStepRATTLELangevin" << endl; + m_pdata->getBoxChangeSignal() + .template connect, + &TwoStepRATTLELangevin::setBoxChange>(this); + if (!m_manifold.fitsInsideBox(m_pdata->getGlobalBox())) { throw std::runtime_error("Parts of the manifold are outside the box"); } } +template TwoStepRATTLELangevin::~TwoStepRATTLELangevin() + { + m_pdata->getBoxChangeSignal() + .template disconnect, + &TwoStepRATTLELangevin::setBoxChange>(this); + m_exec_conf->msg->notice(5) << "Destroying TwoStepRATTLELangevin" << endl; + } + /*! \param timestep Current time step \post Particle positions are moved forward to timestep+1 and velocities to timestep+1/2 per the velocity verlet method. @@ -165,9 +181,13 @@ template void TwoStepRATTLELangevin::integrateStepOne( const BoxDim& box = m_pdata->getBox(); - if (!m_manifold.fitsInsideBox(m_pdata->getGlobalBox())) + if (m_box_changed) { - throw std::runtime_error("Parts of the manifold are outside the box"); + if (!m_manifold.fitsInsideBox(m_pdata->getGlobalBox())) + { + throw std::runtime_error("Parts of the manifold are outside the box"); + } + m_box_changed = false; } // perform the first half step of the RATTLE algorithm applied on velocity verlet diff --git a/hoomd/md/TwoStepRATTLELangevinGPU.h b/hoomd/md/TwoStepRATTLELangevinGPU.h index d701a58c30..a4d01dda7d 100644 --- a/hoomd/md/TwoStepRATTLELangevinGPU.h +++ b/hoomd/md/TwoStepRATTLELangevinGPU.h @@ -130,9 +130,13 @@ void TwoStepRATTLELangevinGPU::integrateStepOne(unsigned int timestep) if (this->m_prof) this->m_prof->push(this->m_exec_conf, "RATTLELangevin step 1"); - if (!this->m_manifold.fitsInsideBox(this->m_pdata->getGlobalBox())) + if (this->m_box_changed) { - throw std::runtime_error("Parts of the manifold are outside the box"); + if (!this->m_manifold.fitsInsideBox(this->m_pdata->getGlobalBox())) + { + throw std::runtime_error("Parts of the manifold are outside the box"); + } + this->m_box_changed = false; } // access all the needed data diff --git a/hoomd/md/TwoStepRATTLENVE.h b/hoomd/md/TwoStepRATTLENVE.h index f283406723..c53df354ab 100644 --- a/hoomd/md/TwoStepRATTLENVE.h +++ b/hoomd/md/TwoStepRATTLENVE.h @@ -48,10 +48,7 @@ template class PYBIND11_EXPORT TwoStepRATTLENVE : public Integra bool skip_restart, Scalar tolerance); - virtual ~TwoStepRATTLENVE() - { - m_exec_conf->msg->notice(5) << "Destroying TwoStepRATTLENVE" << endl; - } + virtual ~TwoStepRATTLENVE(); //! Sets the movement limit void setLimit(Scalar limit) @@ -107,12 +104,19 @@ template class PYBIND11_EXPORT TwoStepRATTLENVE : public Integra }; protected: + //! Helper function to be called when box changes + void setBoxChange() + { + m_box_changed = true; + } + Manifold m_manifold; //!< The manifold used for the RATTLE constraint bool m_limit; //!< True if we should limit the distance a particle moves in one step Scalar m_limit_val; //!< The maximum distance a particle is to move in one step Scalar m_tolerance; //!< The tolerance value of the RATTLE algorithm, setting the tolerance to //!< the manifold bool m_zero_force; //!< True if the integration step should ignore computed forces + bool m_box_changed; }; /*! \file TwoStepRATTLENVE.h @@ -132,10 +136,14 @@ TwoStepRATTLENVE::TwoStepRATTLENVE(std::shared_ptr s bool skip_restart, Scalar tolerance) : IntegrationMethodTwoStep(sysdef, group), m_manifold(manifold), m_limit(false), - m_limit_val(1.0), m_tolerance(tolerance), m_zero_force(false) + m_limit_val(1.0), m_tolerance(tolerance), m_zero_force(false), m_box_changed(false) { m_exec_conf->msg->notice(5) << "Constructing TwoStepRATTLENVE" << endl; + m_pdata->getBoxChangeSignal() + .template connect, &TwoStepRATTLENVE::setBoxChange>( + this); + if (!m_manifold.fitsInsideBox(m_pdata->getGlobalBox())) { throw std::runtime_error("Parts of the manifold are outside the box"); @@ -159,6 +167,14 @@ TwoStepRATTLENVE::TwoStepRATTLENVE(std::shared_ptr s } } +template TwoStepRATTLENVE::~TwoStepRATTLENVE() + { + m_pdata->getBoxChangeSignal() + .template disconnect, &TwoStepRATTLENVE::setBoxChange>( + this); + m_exec_conf->msg->notice(5) << "Destroying TwoStepRATTLENVE" << endl; + } + /*! \param timestep Current time step \post Particle positions are moved forward to timestep+1 and velocities to timestep+1/2 per the velocity verlet method. @@ -183,9 +199,13 @@ template void TwoStepRATTLENVE::integrateStepOne(uint6 const BoxDim& box = m_pdata->getBox(); - if (!m_manifold.fitsInsideBox(m_pdata->getGlobalBox())) + if (m_box_changed) { - throw std::runtime_error("Parts of the manifold are outside the box"); + if (!m_manifold.fitsInsideBox(m_pdata->getGlobalBox())) + { + throw std::runtime_error("Parts of the manifold are outside the box"); + } + m_box_changed = false; } // perform the first half step of the RATTLE algorithm applied on velocity verlet diff --git a/hoomd/md/TwoStepRATTLENVEGPU.h b/hoomd/md/TwoStepRATTLENVEGPU.h index 038a7e12ed..c001c46ce0 100644 --- a/hoomd/md/TwoStepRATTLENVEGPU.h +++ b/hoomd/md/TwoStepRATTLENVEGPU.h @@ -137,9 +137,13 @@ template void TwoStepRATTLENVEGPU::integrateStepOne(un access_location::device, access_mode::readwrite); - if (!this->m_manifold.fitsInsideBox(this->m_pdata->getGlobalBox())) + if (this->m_box_changed) { - throw std::runtime_error("Parts of the manifold are outside the box"); + if (!this->m_manifold.fitsInsideBox(this->m_pdata->getGlobalBox())) + { + throw std::runtime_error("Parts of the manifold are outside the box"); + } + this->m_box_changed = false; } ArrayHandle d_index_array(this->m_group->getIndexArray(), diff --git a/hoomd/md/all_kernels_active_diamond_manifold.cu b/hoomd/md/all_kernels_active_diamond_manifold.cu new file mode 100644 index 0000000000..e046d2e217 --- /dev/null +++ b/hoomd/md/all_kernels_active_diamond_manifold.cu @@ -0,0 +1,27 @@ +// Copyright (c) 2009-2019 The Regents of the University of Michigan +// This file is part of the HOOMD-blue project, released under the BSD 3-Clause License. + +#include "ActiveForceConstraintComputeGPU.cuh" +#include "ManifoldDiamond.h" + +template hipError_t +gpu_compute_active_force_set_constraints(const unsigned int group_size, + unsigned int* d_index_array, + const Scalar4* d_pos, + Scalar4* d_orientation, + const Scalar4* d_f_act, + ManifoldDiamond manifold, + unsigned int block_size); + +template hipError_t gpu_compute_active_force_constraint_rotational_diffusion( + const unsigned int group_size, + unsigned int* d_tag, + unsigned int* d_index_array, + const Scalar4* d_pos, + Scalar4* d_orientation, + ManifoldDiamond manifold, + bool is2D, + const Scalar rotationDiff, + const uint64_t timestep, + const uint16_t seed, + unsigned int block_size); diff --git a/hoomd/md/all_kernels_active_ellipsoid_manifold.cu b/hoomd/md/all_kernels_active_ellipsoid_manifold.cu new file mode 100644 index 0000000000..be7611d553 --- /dev/null +++ b/hoomd/md/all_kernels_active_ellipsoid_manifold.cu @@ -0,0 +1,27 @@ +// Copyright (c) 2009-2019 The Regents of the University of Michigan +// This file is part of the HOOMD-blue project, released under the BSD 3-Clause License. + +#include "ActiveForceConstraintComputeGPU.cuh" +#include "ManifoldEllipsoid.h" + +template hipError_t +gpu_compute_active_force_set_constraints(const unsigned int group_size, + unsigned int* d_index_array, + const Scalar4* d_pos, + Scalar4* d_orientation, + const Scalar4* d_f_act, + ManifoldEllipsoid manifold, + unsigned int block_size); + +template hipError_t gpu_compute_active_force_constraint_rotational_diffusion( + const unsigned int group_size, + unsigned int* d_tag, + unsigned int* d_index_array, + const Scalar4* d_pos, + Scalar4* d_orientation, + ManifoldEllipsoid manifold, + bool is2D, + const Scalar rotationDiff, + const uint64_t timestep, + const uint16_t seed, + unsigned int block_size); diff --git a/hoomd/md/all_kernels_active_gyroid_manifold.cu b/hoomd/md/all_kernels_active_gyroid_manifold.cu new file mode 100644 index 0000000000..636dad1035 --- /dev/null +++ b/hoomd/md/all_kernels_active_gyroid_manifold.cu @@ -0,0 +1,27 @@ +// Copyright (c) 2009-2019 The Regents of the University of Michigan +// This file is part of the HOOMD-blue project, released under the BSD 3-Clause License. + +#include "ActiveForceConstraintComputeGPU.cuh" +#include "ManifoldGyroid.h" + +template hipError_t +gpu_compute_active_force_set_constraints(const unsigned int group_size, + unsigned int* d_index_array, + const Scalar4* d_pos, + Scalar4* d_orientation, + const Scalar4* d_f_act, + ManifoldGyroid manifold, + unsigned int block_size); + +template hipError_t gpu_compute_active_force_constraint_rotational_diffusion( + const unsigned int group_size, + unsigned int* d_tag, + unsigned int* d_index_array, + const Scalar4* d_pos, + Scalar4* d_orientation, + ManifoldGyroid manifold, + bool is2D, + const Scalar rotationDiff, + const uint64_t timestep, + const uint16_t seed, + unsigned int block_size); diff --git a/hoomd/md/all_kernels_active_primitive_manifold.cu b/hoomd/md/all_kernels_active_primitive_manifold.cu new file mode 100644 index 0000000000..a9a8e4aad5 --- /dev/null +++ b/hoomd/md/all_kernels_active_primitive_manifold.cu @@ -0,0 +1,27 @@ +// Copyright (c) 2009-2019 The Regents of the University of Michigan +// This file is part of the HOOMD-blue project, released under the BSD 3-Clause License. + +#include "ActiveForceConstraintComputeGPU.cuh" +#include "ManifoldPrimitive.h" + +template hipError_t +gpu_compute_active_force_set_constraints(const unsigned int group_size, + unsigned int* d_index_array, + const Scalar4* d_pos, + Scalar4* d_orientation, + const Scalar4* d_f_act, + ManifoldPrimitive manifold, + unsigned int block_size); + +template hipError_t gpu_compute_active_force_constraint_rotational_diffusion( + const unsigned int group_size, + unsigned int* d_tag, + unsigned int* d_index_array, + const Scalar4* d_pos, + Scalar4* d_orientation, + ManifoldPrimitive manifold, + bool is2D, + const Scalar rotationDiff, + const uint64_t timestep, + const uint16_t seed, + unsigned int block_size); diff --git a/hoomd/md/all_kernels_active_sphere_manifold.cu b/hoomd/md/all_kernels_active_sphere_manifold.cu new file mode 100644 index 0000000000..f36ee01c8a --- /dev/null +++ b/hoomd/md/all_kernels_active_sphere_manifold.cu @@ -0,0 +1,27 @@ +// Copyright (c) 2009-2019 The Regents of the University of Michigan +// This file is part of the HOOMD-blue project, released under the BSD 3-Clause License. + +#include "ActiveForceConstraintComputeGPU.cuh" +#include "ManifoldSphere.h" + +template hipError_t +gpu_compute_active_force_set_constraints(const unsigned int group_size, + unsigned int* d_index_array, + const Scalar4* d_pos, + Scalar4* d_orientation, + const Scalar4* d_f_act, + ManifoldSphere manifold, + unsigned int block_size); + +template hipError_t gpu_compute_active_force_constraint_rotational_diffusion( + const unsigned int group_size, + unsigned int* d_tag, + unsigned int* d_index_array, + const Scalar4* d_pos, + Scalar4* d_orientation, + ManifoldSphere manifold, + bool is2D, + const Scalar rotationDiff, + const uint64_t timestep, + const uint16_t seed, + unsigned int block_size); diff --git a/hoomd/md/all_kernels_active_xyplane_manifold.cu b/hoomd/md/all_kernels_active_xyplane_manifold.cu new file mode 100644 index 0000000000..ccdeb707d3 --- /dev/null +++ b/hoomd/md/all_kernels_active_xyplane_manifold.cu @@ -0,0 +1,27 @@ +// Copyright (c) 2009-2019 The Regents of the University of Michigan +// This file is part of the HOOMD-blue project, released under the BSD 3-Clause License. + +#include "ActiveForceConstraintComputeGPU.cuh" +#include "ManifoldXYPlane.h" + +template hipError_t +gpu_compute_active_force_set_constraints(const unsigned int group_size, + unsigned int* d_index_array, + const Scalar4* d_pos, + Scalar4* d_orientation, + const Scalar4* d_f_act, + ManifoldXYPlane manifold, + unsigned int block_size); + +template hipError_t gpu_compute_active_force_constraint_rotational_diffusion( + const unsigned int group_size, + unsigned int* d_tag, + unsigned int* d_index_array, + const Scalar4* d_pos, + Scalar4* d_orientation, + ManifoldXYPlane manifold, + bool is2D, + const Scalar rotationDiff, + const uint64_t timestep, + const uint16_t seed, + unsigned int block_size); diff --git a/hoomd/md/all_kernels_active_zcylinder_manifold.cu b/hoomd/md/all_kernels_active_zcylinder_manifold.cu new file mode 100644 index 0000000000..1d2c3985f9 --- /dev/null +++ b/hoomd/md/all_kernels_active_zcylinder_manifold.cu @@ -0,0 +1,27 @@ +// Copyright (c) 2009-2019 The Regents of the University of Michigan +// This file is part of the HOOMD-blue project, released under the BSD 3-Clause License. + +#include "ActiveForceConstraintComputeGPU.cuh" +#include "ManifoldZCylinder.h" + +template hipError_t +gpu_compute_active_force_set_constraints(const unsigned int group_size, + unsigned int* d_index_array, + const Scalar4* d_pos, + Scalar4* d_orientation, + const Scalar4* d_f_act, + ManifoldZCylinder manifold, + unsigned int block_size); + +template hipError_t gpu_compute_active_force_constraint_rotational_diffusion( + const unsigned int group_size, + unsigned int* d_tag, + unsigned int* d_index_array, + const Scalar4* d_pos, + Scalar4* d_orientation, + ManifoldZCylinder manifold, + bool is2D, + const Scalar rotationDiff, + const uint64_t timestep, + const uint16_t seed, + unsigned int block_size); diff --git a/hoomd/md/force.py b/hoomd/md/force.py index 3b93144c2a..8dd3c32eea 100644 --- a/hoomd/md/force.py +++ b/hoomd/md/force.py @@ -13,22 +13,10 @@ from hoomd.data.typeconverter import OnlyTypes from hoomd.data.parameterdicts import ParameterDict, TypeParameterDict from hoomd.filter import ParticleFilter -from hoomd.md.constrain import Constraint +from hoomd.md.manifold import Manifold import numpy -def _ellip_preprocessing(constraint): - if constraint is not None: - if constraint.__class__.__name__ == "constraint_ellipsoid": - return constraint - else: - raise RuntimeError( - "Active force constraint is not accepted (currently only " - "accepts ellipsoids)") - else: - return None - - class _force: # noqa - This will be removed eventually. Needed to build docs. pass @@ -325,11 +313,9 @@ class Active(Force): It is defined per particle type and stays constant during the simulation. active_torque (tuple): active torque vector in reference to the - orientation of a particle - :math:`[\mathrm{force} \cdot \mathrm{length}]`. - It is defined per particle type and stays constant - during the simulation. - + orientation of a particle :math:`[\mathrm{force} \cdot + \mathrm{length}]`. It is defined per particle type and stays + constant during the simulation. :py:class:`Active` specifies that an active force should be added to all particles. Obeys :math:`\delta {\bf r}_i = \delta t v_0 \hat{p}_i`, where @@ -343,12 +329,12 @@ class Active(Force): diffusion follows :math:`\delta \hat{p}_i / \delta t = \sqrt{2 D_r / \delta t} \Gamma (\hat{p}_i (\cos \theta - 1) + \hat{p}_r \sin \theta)`, where :math:`\hat{p}_r` is an uncorrelated random unit vector. The persistence - length of an active particle's path is :math:`v_0 / D_r`. The rotational + length of an active particle's path is :math:`v_0 / D_r`. The rotational diffusion is applied to the orientation vector/quaternion of each particle. - This implies that both the active force and the active torque vectors in the - particle frame stay constant during the simulation. Hence, the active forces - in the system frame are composed of the forces in particle frame and the - current orientation of the particle. + This implies that both the active force and the active torque vectors in + the particle frame stay constant during the simulation. Hence, the active + forces in the system frame are composed of the forces in particle frame + and the current orientation of the particle. Examples:: @@ -363,19 +349,9 @@ class Active(Force): def __init__(self, filter, rotation_diff=0.1): # store metadata - param_dict = ParameterDict( - filter=ParticleFilter, - rotation_diff=float(rotation_diff), - constraint=OnlyTypes(Constraint, - allow_none=True, - preprocess=_ellip_preprocessing), - ) - param_dict.update( - dict( - constraint=None, - rotation_diff=rotation_diff, - filter=filter, - )) + param_dict = ParameterDict(filter=ParticleFilter, + rotation_diff=float(rotation_diff)) + param_dict.update(dict(rotation_diff=rotation_diff, filter=filter)) # set defaults self._param_dict.update(param_dict) @@ -403,21 +379,141 @@ def _add(self, simulation): super()._add(simulation) def _attach(self): + # initialize the reflected c++ class - if isinstance(self._simulation.device, hoomd.device.CPU): + sim = self._simulation + + if isinstance(sim.device, hoomd.device.CPU): my_class = _md.ActiveForceCompute else: my_class = _md.ActiveForceComputeGPU - self._cpp_obj = my_class( - self._simulation.state._cpp_sys_def, - self._simulation.state._get_group(self.filter), - self.rotation_diff, - _hoomd.make_scalar3(0, 0, 0), - 0, - 0, - 0, + self._cpp_obj = my_class(sim.state._cpp_sys_def, + sim.state._get_group(self.filter), + self.rotation_diff) + + # Attach param_dict and typeparam_dict + super()._attach() + + +class ActiveOnManifold(Force): + r"""Active force on a manifold. + + Attributes: + filter (:py:mod:`hoomd.filter`): Subset of particles on which to apply + active forces. + manifold_constraint (:py:mod:`hoomd.md.manifold.Manifold`): Manifold + constraint. + rotation_diff (float): rotational diffusion constant, :math:`D_r`, for + all particles in the group + :math:`[\mathrm{radian}^{2} \cdot \mathrm{time}^{-1}]`. + active_force (tuple): active force vector in reference to the + orientation of a particle :math:`[\mathrm{force}]`. + It is defined per particle type and stays constant during + the simulation. + active_torque (tuple): active torque vector in reference to the + orientation of a particle :math:`[\mathrm{force} \cdot + \mathrm{length}]`. It is defined per particle type and stays + constant during the simulation. + + :py:class:`ActiveOnManifold` specifies that a constrained active force + should be added to all particles similar to :py:class:`Active`. The + active force vector :math:`\hat{p}_i` is restricted to the local tangent + plane of the manifold constraint at point :math:`{\bf r}_i`. For more + information see :py:class:`Active`. + + Hint: + Use `ActiveOnManifold` with a `md.methods.rattle` integration method + with the same manifold constraint. + + Examples:: + + + all = filter.All() + sphere = hoomd.md.manifold.Sphere(r=10) + active = hoomd.md.force.ActiveOnManifold( + filter=hoomd.filter.All(), rotation_diff=0.01, + manifold_constraint = sphere + ) + active.active_force['A','B'] = (1,0,0) + active.active_torque['A','B'] = (0,0,0) + """ + + def __init__(self, filter, manifold_constraint, rotation_diff=0.1): + # store metadata + param_dict = ParameterDict(filter=ParticleFilter, + rotation_diff=float(rotation_diff), + manifold_constraint=OnlyTypes( + Manifold, allow_none=False)) + param_dict.update( + dict(rotation_diff=rotation_diff, + filter=filter, + manifold_constraint=manifold_constraint)) + # set defaults + self._param_dict.update(param_dict) + + active_force = TypeParameter( + "active_force", + type_kind="particle_types", + param_dict=TypeParameterDict((1.0, 0.0, 0.0), len_keys=1), ) + active_torque = TypeParameter( + "active_torque", + type_kind="particle_types", + param_dict=TypeParameterDict((0.0, 0.0, 0.0), len_keys=1), + ) + + self._extend_typeparam([active_force, active_torque]) + + def _getattr_param(self, attr): + if self._attached: + if attr == "manifold_constraint": + return self._param_dict["manifold_constraint"] + parameter = getattr(self._cpp_obj, attr) + return parameter + else: + return self._param_dict[attr] + + def _setattr_param(self, attr, value): + if attr == "manifold_constraint": + raise AttributeError( + "Cannot set manifold_constraint after construction.") + super()._setattr_param(attr, value) + + def _add(self, simulation): + """Add the operation to a simulation. + + Active forces use RNGs. Warn the user if they did not set the seed. + """ + if simulation is not None: + simulation._warn_if_seed_unset() + + if self.manifold_constraint is not None: + self.manifold_constraint._add(simulation) + + super()._add(simulation) + + def _attach(self): + + # initialize the reflected c++ class + sim = self._simulation + + if not self.manifold_constraint._attached: + self.manifold_constraint._attach() + + if isinstance(sim.device, hoomd.device.CPU): + my_class = getattr( + _md, 'ActiveForceConstraintCompute' + + self.manifold_constraint.__class__.__name__) + else: + my_class = getattr( + _md, 'ActiveForceConstraintCompute' + + self.manifold_constraint.__class__.__name__ + 'GPU') + + self._cpp_obj = my_class(sim.state._cpp_sys_def, + sim.state._get_group(self.filter), + self.rotation_diff, + self.manifold_constraint._cpp_obj) # Attach param_dict and typeparam_dict super()._attach() diff --git a/hoomd/md/methods.py b/hoomd/md/methods.py index 1e0f33fad9..cb41d3d022 100644 --- a/hoomd/md/methods.py +++ b/hoomd/md/methods.py @@ -769,7 +769,6 @@ class NVE(MethodRATTLE): tolerance (float): Defines the tolerated error particles are allowed to deviate from the manifold in terms of the implicit function. Defaults to 1e-6. - """ def __init__(self, diff --git a/hoomd/md/module-md.cc b/hoomd/md/module-md.cc index c602d41b70..040b44fe39 100644 --- a/hoomd/md/module-md.cc +++ b/hoomd/md/module-md.cc @@ -4,6 +4,7 @@ // Maintainer: joaander All developers are free to add the calls needed to export their modules #include "ActiveForceCompute.h" +#include "ActiveForceConstraintCompute.h" #include "AllAnisoPairPotentials.h" #include "AllBondPotentials.h" #include "AllExternalPotentials.h" @@ -14,8 +15,6 @@ #include "BondTablePotential.h" #include "ComputeThermo.h" #include "ComputeThermoHMA.h" -#include "ConstraintEllipsoid.h" -#include "ConstraintSphere.h" #include "CosineSqAngleForceCompute.h" #include "EvaluatorRevCross.h" #include "EvaluatorSquareDensity.h" @@ -42,7 +41,6 @@ #include "NeighborListStencil.h" #include "NeighborListTree.h" #include "OPLSDihedralForceCompute.h" -#include "OneDConstraint.h" #include "PPPMForceCompute.h" #include "PotentialBond.h" #include "PotentialExternal.h" @@ -68,12 +66,11 @@ // include GPU classes #ifdef ENABLE_HIP #include "ActiveForceComputeGPU.h" +#include "ActiveForceConstraintComputeGPU.h" #include "AnisoPotentialPairGPU.h" #include "BondTablePotentialGPU.h" #include "ComputeThermoGPU.h" #include "ComputeThermoHMAGPU.h" -#include "ConstraintEllipsoidGPU.h" -#include "ConstraintSphereGPU.h" #include "CosineSqAngleForceComputeGPU.h" #include "FIREEnergyMinimizerGPU.h" #include "ForceCompositeGPU.h" @@ -87,7 +84,6 @@ #include "NeighborListGPUStencil.h" #include "NeighborListGPUTree.h" #include "OPLSDihedralForceComputeGPU.h" -#include "OneDConstraintGPU.h" #include "PPPMForceComputeGPU.h" #include "PotentialBondGPU.h" #include "PotentialExternalGPU.h" @@ -233,6 +229,16 @@ void export_PotentialExternal(pybind11::module& PYBIND11_MODULE(_md, m) { export_ActiveForceCompute(m); + export_ActiveForceConstraintCompute(m, + "ActiveForceConstraintComputeCylinder"); + export_ActiveForceConstraintCompute(m, "ActiveForceConstraintComputeDiamond"); + export_ActiveForceConstraintCompute(m, + "ActiveForceConstraintComputeEllipsoid"); + export_ActiveForceConstraintCompute(m, "ActiveForceConstraintComputeGyroid"); + export_ActiveForceConstraintCompute(m, "ActiveForceConstraintComputePlane"); + export_ActiveForceConstraintCompute(m, + "ActiveForceConstraintComputePrimitive"); + export_ActiveForceConstraintCompute(m, "ActiveForceConstraintComputeSphere"); export_ComputeThermo(m); export_ComputeThermoHMA(m); export_HarmonicAngleForceCompute(m); @@ -285,8 +291,6 @@ PYBIND11_MODULE(_md, m) export_NeighborListBinned(m); export_NeighborListStencil(m); export_NeighborListTree(m); - export_ConstraintSphere(m); - export_OneDConstraint(m); export_MolecularForceCompute(m); export_ForceDistanceConstraint(m); export_ForceComposite(m); @@ -388,13 +392,29 @@ PYBIND11_MODULE(_md, m) export_OPLSDihedralForceComputeGPU(m); export_TableDihedralForceComputeGPU(m); export_HarmonicImproperForceComputeGPU(m); - export_ConstraintSphereGPU(m); - export_OneDConstraintGPU(m); export_ForceDistanceConstraintGPU(m); export_ComputeThermoGPU(m); export_ComputeThermoHMAGPU(m); export_PPPMForceComputeGPU(m); export_ActiveForceComputeGPU(m); + export_ActiveForceConstraintComputeGPU( + m, + "ActiveForceConstraintComputeCylinderGPU"); + export_ActiveForceConstraintComputeGPU( + m, + "ActiveForceConstraintComputeDiamondGPU"); + export_ActiveForceConstraintComputeGPU( + m, + "ActiveForceConstraintComputeEllipsoidGPU"); + export_ActiveForceConstraintComputeGPU(m, + "ActiveForceConstraintComputeGyroidGPU"); + export_ActiveForceConstraintComputeGPU(m, + "ActiveForceConstraintComputePlaneGPU"); + export_ActiveForceConstraintComputeGPU( + m, + "ActiveForceConstraintComputePrimitiveGPU"); + export_ActiveForceConstraintComputeGPU(m, + "ActiveForceConstraintComputeSphereGPU"); export_PotentialExternalGPU( m, "PotentialExternalPeriodicGPU"); @@ -431,7 +451,6 @@ PYBIND11_MODULE(_md, m) export_TwoStepBD(m); export_TwoStepNPTMTK(m); export_Berendsen(m); - export_ConstraintEllipsoid(m); export_FIREEnergyMinimizer(m); export_MuellerPlatheFlow(m); @@ -468,7 +487,6 @@ PYBIND11_MODULE(_md, m) export_TwoStepNPTMTKGPU(m); export_BerendsenGPU(m); export_FIREEnergyMinimizerGPU(m); - export_ConstraintEllipsoidGPU(m); export_MuellerPlatheFlowGPU(m); export_TwoStepRATTLEBDGPU(m, "TwoStepRATTLEBDCylinderGPU"); diff --git a/hoomd/md/pytest/test_active.py b/hoomd/md/pytest/test_active.py index fa9eae5bb6..f97eae80c4 100644 --- a/hoomd/md/pytest/test_active.py +++ b/hoomd/md/pytest/test_active.py @@ -1,17 +1,101 @@ import hoomd +import pytest from hoomd.conftest import pickling_check +def test_attributes(): + active = hoomd.md.force.Active(filter=hoomd.filter.All(), + rotation_diff=0.01) + + assert active.rotation_diff == 0.01 + assert active.active_force['A'] == (1.0, 0.0, 0.0) + assert active.active_torque['A'] == (0.0, 0.0, 0.0) + + active.rotation_diff = 0.1 + assert active.rotation_diff == 0.1 + active.active_force['A'] = (0.5, 0.0, 0.0) + assert active.active_force['A'] == (0.5, 0.0, 0.0) + active.active_force['A'] = (0.0, 0.0, 1.0) + assert active.active_force['A'] == (0.0, 0.0, 1.0) + + +def test_attributes_constraints(): + plane = hoomd.md.manifold.Plane() + active = hoomd.md.force.ActiveOnManifold(filter=hoomd.filter.All(), + rotation_diff=0.01, + manifold_constraint=plane) + + assert active.rotation_diff == 0.01 + assert active.active_force['A'] == (1.0, 0.0, 0.0) + assert active.active_torque['A'] == (0.0, 0.0, 0.0) + assert active.manifold_constraint == plane + + active.rotation_diff = 0.1 + assert active.rotation_diff == 0.1 + active.active_force['A'] = (0.5, 0.0, 0.0) + assert active.active_force['A'] == (0.5, 0.0, 0.0) + active.active_force['A'] = (0.0, 0.0, 1.0) + assert active.active_force['A'] == (0.0, 0.0, 1.0) + + sphere = hoomd.md.manifold.Sphere(r=5) + with pytest.raises(AttributeError): + active.manifold_constraint = sphere + assert active.manifold_constraint == plane + + def test_attach(simulation_factory, two_particle_snapshot_factory): + active = hoomd.md.force.Active(filter=hoomd.filter.All(), + rotation_diff=0.01) + sim = simulation_factory(two_particle_snapshot_factory(dimensions=3, d=8)) integrator = hoomd.md.Integrator(.05) integrator.methods.append( hoomd.md.methods.Langevin(hoomd.filter.All(), kT=0)) - integrator.forces.append( - hoomd.md.force.Active(filter=hoomd.filter.All(), rotation_diff=0.01)) + integrator.forces.append(active) sim.operations.integrator = integrator - sim.operations._schedule() - sim.run(10) + sim.run(0) + + assert active.rotation_diff == 0.01 + assert active.active_force['A'] == (1.0, 0.0, 0.0) + assert active.active_torque['A'] == (0.0, 0.0, 0.0) + + active.rotation_diff = 0.1 + assert active.rotation_diff == 0.1 + active.active_force['A'] = (0.5, 0.0, 0.0) + assert active.active_force['A'] == (0.5, 0.0, 0.0) + active.active_force['A'] = (0.0, 0.0, 1.0) + assert active.active_force['A'] == (0.0, 0.0, 1.0) + + +def test_attach_manifold(simulation_factory, two_particle_snapshot_factory): + plane = hoomd.md.manifold.Plane() + active = hoomd.md.force.ActiveOnManifold(filter=hoomd.filter.All(), + rotation_diff=0.01, + manifold_constraint=plane) + + sim = simulation_factory(two_particle_snapshot_factory(dimensions=3, d=8)) + integrator = hoomd.md.Integrator(.05) + integrator.methods.append( + hoomd.md.methods.Langevin(hoomd.filter.All(), kT=0)) + integrator.forces.append(active) + sim.operations.integrator = integrator + sim.run(0) + + assert active.rotation_diff == 0.01 + assert active.active_force['A'] == (1.0, 0.0, 0.0) + assert active.active_torque['A'] == (0.0, 0.0, 0.0) + assert active.manifold_constraint == plane + + active.rotation_diff = 0.1 + assert active.rotation_diff == 0.1 + active.active_force['A'] = (0.5, 0.0, 0.0) + assert active.active_force['A'] == (0.5, 0.0, 0.0) + active.active_force['A'] = (0.0, 0.0, 1.0) + assert active.active_force['A'] == (0.0, 0.0, 1.0) + sphere = hoomd.md.manifold.Sphere(r=2) + with pytest.raises(AttributeError): + active.manifold_constraint = sphere + assert active.manifold_constraint == plane def test_pickling(simulation_factory, two_particle_snapshot_factory): @@ -26,3 +110,19 @@ def test_pickling(simulation_factory, two_particle_snapshot_factory): sim.operations.integrator = integrator sim.run(0) pickling_check(active) + + +def test_pickling_constraint(simulation_factory, two_particle_snapshot_factory): + sim = simulation_factory(two_particle_snapshot_factory()) + active = hoomd.md.force.ActiveOnManifold( + filter=hoomd.filter.All(), + rotation_diff=0.01, + manifold_constraint=hoomd.md.manifold.Plane()) + pickling_check(active) + integrator = hoomd.md.Integrator( + .05, + methods=[hoomd.md.methods.Langevin(hoomd.filter.All(), kT=0)], + forces=[active]) + sim.operations.integrator = integrator + sim.run(0) + pickling_check(active) diff --git a/hoomd/md/update.py b/hoomd/md/update.py index 0f0ce3f590..f77e8da11f 100644 --- a/hoomd/md/update.py +++ b/hoomd/md/update.py @@ -11,7 +11,6 @@ is triggered to change its state. """ -from hoomd import _hoomd from hoomd.md import _md import hoomd from hoomd.operation import Updater @@ -46,93 +45,6 @@ def _attach(self): super()._attach() -class constraint_ellipsoid: # noqa: N801 (this will be removed) - """Constrain particles to the surface of a ellipsoid. - - Args: - group (``hoomd.group``): Group for which the update will be set - P (tuple): (x,y,z) tuple indicating the position of the center of the - ellipsoid (in distance units). - rx (float): radius of an ellipsoid in the X direction (in distance - units). - ry (float): radius of an ellipsoid in the Y direction (in distance - units). - rz (float): radius of an ellipsoid in the Z direction (in distance - units). - r (float): radius of a sphere (in distance units), such that r=rx=ry=rz. - - `constraint_ellipsoid` specifies that all particles are constrained to the - surface of an ellipsoid. Each time step particles are projected onto the - surface of the ellipsoid. Method from: - http://www.geometrictools.com/Documentation/\ - DistancePointEllipseEllipsoid.pdf - - .. attention:: - For the algorithm to work, we must have :math:`rx >= rz,~ry >= rz,~rz > - 0`. - - Note: - This method does not properly conserve virial coefficients. - - Note: - random thermal forces from the integrator are applied in 3D not 2D, - therefore they aren't fully accurate. Suggested use is therefore only - for T=0. - - Examples:: - - update.constraint_ellipsoid(P=(-1,5,0), r=9) - update.constraint_ellipsoid(rx=7, ry=5, rz=3) - - """ - - def __init__(self, group, r=None, rx=None, ry=None, rz=None, P=(0, 0, 0)): - period = 1 - - # Error out in MPI simulations - if (hoomd.version.mpi_enabled): - if hoomd.context.current.system_definition.getParticleData( - ).getDomainDecomposition(): - hoomd.context.current.device.cpp_msg.error( - "constrain.ellipsoid is not supported in multi-processor " - "simulations.\n\n") - raise RuntimeError("Error initializing updater.") - - # Error out if no radii are set - if (r is None and rx is None and ry is None and rz is None): - hoomd.context.current.device.cpp_msg.error( - "no radii were defined in update.constraint_ellipsoid.\n\n") - raise RuntimeError("Error initializing updater.") - - # initialize the base class - # _updater.__init__(self) - - # Set parameters - P = _hoomd.make_scalar3(P[0], P[1], P[2]) - if (r is not None): - rx = ry = rz = r - - # create the c++ mirror class - if not hoomd.context.current.device.cpp_exec_conf.isCUDAEnabled(): - self.cpp_updater = _md.ConstraintEllipsoid( - hoomd.context.current.system_definition, group.cpp_group, P, rx, - ry, rz) - else: - self.cpp_updater = _md.ConstraintEllipsoidGPU( - hoomd.context.current.system_definition, group.cpp_group, P, rx, - ry, rz) - - self.setupUpdater(period) - - # store metadata - self.group = group - self.P = P - self.rx = rx - self.ry = ry - self.rz = rz - self.metadata_fields = ['group', 'P', 'rx', 'ry', 'rz'] - - class ReversePerturbationFlow(Updater): """Reverse Perturbation (Müller-Plathe) method to establish shear flow. diff --git a/sphinx-doc/module-md-force.rst b/sphinx-doc/module-md-force.rst index 1bbbd82142..a2be68f5bd 100644 --- a/sphinx-doc/module-md-force.rst +++ b/sphinx-doc/module-md-force.rst @@ -10,6 +10,7 @@ md.force Force Active + ActiveOnManifold .. rubric:: Details @@ -22,3 +23,7 @@ md.force .. autoclass:: Active :show-inheritance: :no-inherited-members: + + .. autoclass:: ActiveOnManifold + :show-inheritance: + :no-inherited-members: