diff --git a/hoomd/BoxDim.h b/hoomd/BoxDim.h index ec0dd86268..0d63a10fca 100644 --- a/hoomd/BoxDim.h +++ b/hoomd/BoxDim.h @@ -80,6 +80,8 @@ struct m_lo = m_hi = m_Linv = m_L = make_scalar3(0, 0, 0); m_xz = m_xy = m_yz = Scalar(0.0); m_periodic = make_uchar3(1, 1, 1); + m_L_rate = make_scalar3(0, 0, 0); + m_xz_rate = m_xy_rate = m_yz_rate = 0; } //! Constructs a box from -Len/2 to Len/2 @@ -92,6 +94,8 @@ struct setL(make_scalar3(Len, Len, Len)); m_periodic = make_uchar3(1, 1, 1); m_xz = m_xy = m_yz = Scalar(0.0); + m_L_rate = make_scalar3(0, 0, 0); + m_xz_rate = m_xy_rate = m_yz_rate = 0; } //! Constructs a box from -Len_x/2 to Len_x/2 for each dimension @@ -105,6 +109,8 @@ struct setL(make_scalar3(Len_x, Len_y, Len_z)); m_periodic = make_uchar3(1, 1, 1); m_xz = m_xy = m_yz = Scalar(0.0); + m_L_rate = make_scalar3(0, 0, 0); + m_xz_rate = m_xy_rate = m_yz_rate = 0; } //! Constructs a box from -L/2 to L/2 for each dimension @@ -116,6 +122,8 @@ struct setL(L); m_periodic = make_uchar3(1, 1, 1); m_xz = m_xy = m_yz = Scalar(0.0); + m_L_rate = make_scalar3(0, 0, 0); + m_xz_rate = m_xy_rate = m_yz_rate = 0; } //! Constructs a tilted box with edges of length len for each dimension @@ -129,6 +137,8 @@ struct setL(make_scalar3(Len, Len, Len)); setTiltFactors(xy, xz, yz); m_periodic = make_uchar3(1, 1, 1); + m_L_rate = make_scalar3(0, 0, 0); + m_xz_rate = m_xy_rate = m_yz_rate = 0; } //! Construct a box from specific lo and hi values @@ -141,6 +151,8 @@ struct setLoHi(lo, hi); m_periodic = periodic; m_xz = m_xy = m_yz = Scalar(0.0); + m_L_rate = make_scalar3(0, 0, 0); + m_xz_rate = m_xy_rate = m_yz_rate = 0; } /// Constructs a box from a std::array @@ -151,6 +163,8 @@ struct setL(make_scalar3(array[0], array[1], array[2])); setTiltFactors(array[3], array[4], array[5]); m_periodic = make_uchar3(1, 1, 1); + m_L_rate = make_scalar3(0, 0, 0); + m_xz_rate = m_xy_rate = m_yz_rate = 0; } //! Get the periodic flags @@ -262,6 +276,51 @@ struct return m_yz; } + //! Update the box deformation in each dimension + /*! \param L_rate box deformation rate in each dimension + */ + HOSTDEVICE void setLDeformationRate(const Scalar3& L_rate) + { + m_L_rate = L_rate; + } + + //! Return the deformation rates of the box lengths in each dimension + HOSTDEVICE Scalar3 getLDeformationRate() const + { + return m_L_rate; + } + + //! Update the deformation rates in box tilts + /*! \param xy_rate Deformation rate in xy + \param xz_rate Deformation rate in xz + \param yz_rate Deformation rate in yz + */ + HOSTDEVICE void + setTiltDeformationRates(const Scalar xy_rate, const Scalar xz_rate, const Scalar yz_rate) + { + m_xy_rate = xy_rate; + m_xz_rate = xz_rate; + m_yz_rate = yz_rate; + } + + //! Returns the xy deformation rate + HOSTDEVICE Scalar getTiltDeformationRateXY() const + { + return m_xy_rate; + } + + //! Returns the xz deformation rate + HOSTDEVICE Scalar getTiltDeformationRateXZ() const + { + return m_xz_rate; + } + + //! Returns the yz deformation rate + HOSTDEVICE Scalar getTiltDeformationRateYZ() const + { + return m_yz_rate; + } + //! Compute fractional coordinates, allowing for a ghost layer /*! \param v Vector to scale \param ghost_width Width of extra ghost padding layer to take into account (along reciprocal @@ -393,6 +452,120 @@ struct return vec3(minImage(vec_to_scalar3(v))); } + //! Compute minimum image + /*! + \param dr Position vector to compute / Minimum image vector returned + \param dv Velocity vector to compute / Minimum image vector returned + \note \a dr must not extend more than 1 image beyond the box + */ + HOSTDEVICE void minImage(Scalar3& dr, Scalar3& dv) const + { + Scalar3 r = dr; + Scalar3 v = dv; + const Scalar3 L = getL(); + +#ifdef __HIPCC__ + if (m_periodic.z) + { + const Scalar img = slow::rint(r.z * m_Linv.z); + r.z -= L.z * img; + r.y -= L.z * m_yz * img; + r.x -= L.z * m_xz * img; + v.z -= m_L_rate.z * img; + v.y -= L.z * m_yz_rate * img; + v.x -= L.z * m_xz_rate * img; + } + + if (m_periodic.y) + { + const Scalar img = slow::rint(r.y * m_Linv.y); + r.y -= L.y * img; + r.x -= L.y * m_xy * img; + v.y -= m_L_rate.y * img; + v.x -= L.y * m_xy_rate * img; + } + + if (m_periodic.x) + { + const Scalar img = slow::rint(r.x * m_Linv.x); + r.x -= L.x * img; + v.x -= m_L_rate.x * img; + } +#else + // on the cpu, branches are faster than calling rint + if (m_periodic.z) + { + if (r.z >= m_hi.z) + { + r.z -= L.z; + r.y -= L.z * m_yz; + r.x -= L.z * m_xz; + v.z -= m_L_rate.z; + v.y -= L.z * m_yz_rate; + v.x -= L.z * m_xz_rate; + } + else if (r.z < m_lo.z) + { + r.z += L.z; + r.y += L.z * m_yz; + r.x += L.z * m_xz; + v.z += m_L_rate.z; + v.y += L.z * m_yz_rate; + v.x += L.z * m_xz_rate; + } + } + + if (m_periodic.y) + { + if (r.y >= m_hi.y) + { + const int i = static_cast(r.y * m_Linv.y + Scalar(0.5)); + r.y -= (Scalar)i * L.y; + r.x -= (Scalar)i * L.y * m_xy; + v.y -= m_L_rate.y; + v.x -= L.y * m_xy_rate; + } + else if (r.y < m_lo.y) + { + const int i = static_cast(-r.y * m_Linv.y + Scalar(0.5)); + r.y += (Scalar)i * L.y; + r.x += (Scalar)i * L.y * m_xy; + v.y += m_L_rate.y; + v.x += L.y * m_xy_rate; + } + } + + if (m_periodic.x) + { + if (r.x >= m_hi.x) + { + const int i = static_cast(r.x * m_Linv.x + Scalar(0.5)); + r.x -= (Scalar)i * L.x; + v.x -= m_L_rate.x; + } + else if (r.x < m_lo.x) + { + const int i = static_cast(-r.x * m_Linv.x + Scalar(0.5)); + r.x += (Scalar)i * L.x; + v.x += m_L_rate.x; + } + } +#endif + + dr = r; + dv = v; + } + + //! Minimum image using vec3s + HOSTDEVICE void minImage(vec3& dr, vec3& dv) const + { + Scalar3 dr_scalar = vec_to_scalar3(dr); + Scalar3 dv_scalar = vec_to_scalar3(dv); + minImage(dr_scalar, dv_scalar); + dr = vec3(dr_scalar); + dv = vec3(dv_scalar); + } + //! Wrap a vector back into the box /*! \param w Vector to wrap, updated to the minimum image obeying the periodic settings \param img Image of the vector, updated to reflect the new image @@ -496,6 +669,125 @@ struct w.z = v.z; } + //! Wrap vectors back into the box + /*! \param pos Vector to wrap, updated to the minimum image obeying the periodic settings + \param vel Vector to wrap according to deformation rates, obeying the periodic settings + \param img Image of the vector, updated to reflect the new image + \param flags Vector of flags to force wrapping along certain directions + \post \a img and \a pos are updated appropriately + \note \a pos must not extend more than 1 image beyond the box + */ + HOSTDEVICE void + wrap(Scalar3& pos, Scalar3& vel, int3& img, char3 flags = make_char3(0, 0, 0)) const + { + const Scalar3 L = getL(); + + // allow for a shifted box with periodic boundary conditions + Scalar3 origin = (m_hi + m_lo) / Scalar(2.0); + + // tilt factors for nonperiodic boxes are always calculated w.r.t. to the global box with + // origin (0,0,0) + if (!m_periodic.y) + { + origin.y = Scalar(0.0); + } + if (!m_periodic.z) + { + origin.z = Scalar(0.0); + } + + if (m_periodic.x) + { + const Scalar tilt_x + = (m_xz - m_xy * m_yz) * (pos.z - origin.z) + m_xy * (pos.y - origin.y); + if (((pos.x >= m_hi.x + tilt_x) && !flags.x) || flags.x == 1) + { + pos.x -= L.x; + vel.x -= m_L_rate.x; + img.x++; + } + else if (((pos.x < m_lo.x + tilt_x) && !flags.x) || flags.x == -1) + { + pos.x += L.x; + vel.x += m_L_rate.x; + img.x--; + } + } + + if (m_periodic.y) + { + const Scalar tilt_y = m_yz * (pos.z - origin.z); + if (((pos.y >= m_hi.y + tilt_y) && !flags.y) || flags.y == 1) + { + pos.y -= L.y; + pos.x -= L.y * m_xy; + vel.y -= m_L_rate.y; + vel.x -= L.y * m_xy_rate; + img.y++; + } + else if (((pos.y < m_lo.y + tilt_y) && !flags.y) || flags.y == -1) + { + pos.y += L.y; + pos.x += L.y * m_xy; + vel.y += m_L_rate.y; + vel.x += L.y * m_xy_rate; + img.y--; + } + } + + if (m_periodic.z) + { + if (((pos.z >= m_hi.z) && !flags.z) || flags.z == 1) + { + pos.z -= L.z; + pos.y -= L.z * m_yz; + pos.x -= L.z * m_xz; + vel.z -= m_L_rate.z; + vel.y -= L.z * m_yz_rate; + vel.x -= L.z * m_xz_rate; + img.z++; + } + else if (((pos.z < m_lo.z) && !flags.z) || flags.z == -1) + { + pos.z += L.z; + pos.y += L.z * m_yz; + pos.x += L.z * m_xz; + vel.z += m_L_rate.z; + vel.y += L.z * m_yz_rate; + vel.x += L.z * m_xz_rate; + img.z--; + } + } + } + + //! Wrap a vec3 + HOSTDEVICE void + wrap(vec3& pos, vec3& vel, int3& img, char3 flags = make_char3(0, 0, 0)) const + { + Scalar3 pos_scalar = vec_to_scalar3(pos); + Scalar3 vel_scalar = vec_to_scalar3(vel); + wrap(pos_scalar, vel_scalar, img, flags); + pos = vec3(pos_scalar); + vel = vec3(vel_scalar); + } + + //! Wrap a Scalar4 + /*! \note The 4th element remains unchanged + */ + HOSTDEVICE void + wrap(Scalar4& pos, Scalar4& vel, int3& img, char3 flags = make_char3(0, 0, 0)) const + { + Scalar3 p = make_scalar3(pos.x, pos.y, pos.z); + Scalar3 v = make_scalar3(vel.x, vel.y, vel.z); + wrap(p, v, img, flags); + pos.x = p.x; + pos.y = p.y; + pos.z = p.z; + vel.x = v.x; + vel.y = v.y; + vel.z = v.z; + } + //! Get the periodic image a vector belongs to /*! \param v The vector to check \returns the integer coordinates of the periodic image @@ -537,6 +829,29 @@ struct return vec3(shift(vec_to_scalar3(v), _shift)); } + //! Shift a vector by a multiple of the lattice vectors + /*! \param pos The position vector to shift (and the new shifted position) + \param vel The velocity vector to shift (and the new shifted velocity) + \param _shift The displacement in lattice coordinates + */ + HOSTDEVICE void shift(Scalar3& pos, Scalar3& vel, const int3& _shift) const + { + pos = shift(pos, _shift); + vel += Scalar(_shift.x) * make_scalar3(m_L_rate.x, 0.0, 0.0); + vel += Scalar(_shift.y) * make_scalar3(m_L.y * m_xy_rate, m_L_rate.y, 0.0); + vel += Scalar(_shift.z) * make_scalar3(m_L.z * m_xz_rate, m_L.z * m_yz_rate, m_L_rate.z); + } + + //! Shift a vec3 + HOSTDEVICE void shift(vec3& pos, vec3& vel, const int3& _shift) const + { + Scalar3 pos_scalar = vec_to_scalar3(pos); + Scalar3 vel_scalar = vec_to_scalar3(vel); + shift(pos_scalar, vel_scalar, _shift); + pos = vec3(pos_scalar); + vel = vec3(vel_scalar); + } + //! Get the shortest distance between opposite boundary planes of the box /*! The distance between two planes of the lattice is 2 Pi/|b_i|, where * b_1 is the reciprocal lattice vector of the Bravais lattice normal to @@ -639,6 +954,12 @@ struct ar & m_periodic.x; ar & m_periodic.y; ar & m_periodic.z; + ar & m_L_rate.x; + ar & m_L_rate.y; + ar & m_L_rate.z; + ar & m_xy_rate; + ar & m_xz_rate; + ar & m_yz_rate; } #endif @@ -651,6 +972,10 @@ struct Scalar m_xz; //!< xz tilt factor Scalar m_yz; //!< yz tilt factor uchar3 m_periodic; //!< 0/1 in each direction to tell if the box is periodic in that direction + Scalar3 m_L_rate; //!< Deformation rate of box dimensions + Scalar m_xy_rate; //!< Deformation rate in xy + Scalar m_xz_rate; //!< Deformation rate in xz + Scalar m_yz_rate; //!< Deformation rate in yz }; } // end namespace hoomd diff --git a/hoomd/BoxResizeUpdater.cc b/hoomd/BoxResizeUpdater.cc index a90c5d5885..ed389ed272 100644 --- a/hoomd/BoxResizeUpdater.cc +++ b/hoomd/BoxResizeUpdater.cc @@ -96,6 +96,9 @@ void BoxResizeUpdater::scaleAndWrapParticles(const BoxDim& cur_box, const BoxDim // ensure that the particles are still in their // local boxes by wrapping them if they are not + ArrayHandle h_vel(m_pdata->getVelocities(), + access_location::host, + access_mode::readwrite); ArrayHandle h_image(m_pdata->getImages(), access_location::host, access_mode::readwrite); const BoxDim& local_box = m_pdata->getBox(); @@ -104,7 +107,7 @@ void BoxResizeUpdater::scaleAndWrapParticles(const BoxDim& cur_box, const BoxDim { // need to update the image if we move particles from one side // of the box to the other - local_box.wrap(h_pos.data[i], h_image.data[i]); + local_box.wrap(h_pos.data[i], h_vel.data[i], h_image.data[i]); } } diff --git a/hoomd/BoxResizeUpdaterGPU.cc b/hoomd/BoxResizeUpdaterGPU.cc index 047d5ab0c2..bda9a06b92 100644 --- a/hoomd/BoxResizeUpdaterGPU.cc +++ b/hoomd/BoxResizeUpdaterGPU.cc @@ -42,6 +42,10 @@ void BoxResizeUpdaterGPU::scaleAndWrapParticles(const BoxDim& cur_box, const Box access_location::device, access_mode::readwrite); + ArrayHandle d_vel(m_pdata->getVelocities(), + access_location::device, + access_mode::readwrite); + ArrayHandle d_image(m_pdata->getImages(), access_location::device, access_mode::readwrite); @@ -62,6 +66,7 @@ void BoxResizeUpdaterGPU::scaleAndWrapParticles(const BoxDim& cur_box, const Box m_tuner_wrap->begin(); kernel::gpu_box_resize_wrap(m_pdata->getN(), d_pos.data, + d_vel.data, d_image.data, new_box, m_tuner_wrap->getParam()[0]); diff --git a/hoomd/BoxResizeUpdaterGPU.cu b/hoomd/BoxResizeUpdaterGPU.cu index a9576cb7db..faacb191e7 100644 --- a/hoomd/BoxResizeUpdaterGPU.cu +++ b/hoomd/BoxResizeUpdaterGPU.cu @@ -31,14 +31,17 @@ __global__ void gpu_box_resize_scale_kernel(Scalar4* d_pos, } } -__global__ void -gpu_box_resize_wrap_kernel(unsigned int N, Scalar4* d_pos, int3* d_image, const BoxDim new_box) +__global__ void gpu_box_resize_wrap_kernel(unsigned int N, + Scalar4* d_pos, + Scalar4* d_vel, + int3* d_image, + const BoxDim new_box) { int idx = blockIdx.x * blockDim.x + threadIdx.x; if (idx < N) { - new_box.wrap(d_pos[idx], d_image[idx]); + new_box.wrap(d_pos[idx], d_vel[idx], d_image[idx]); } } @@ -74,6 +77,7 @@ hipError_t gpu_box_resize_scale(Scalar4* d_pos, hipError_t gpu_box_resize_wrap(const unsigned int N, Scalar4* d_pos, + Scalar4* d_vel, int3* d_image, const BoxDim& new_box, unsigned int block_size) @@ -94,6 +98,7 @@ hipError_t gpu_box_resize_wrap(const unsigned int N, 0, N, d_pos, + d_vel, d_image, new_box); diff --git a/hoomd/BoxResizeUpdaterGPU.cuh b/hoomd/BoxResizeUpdaterGPU.cuh index 70749c0687..e5052e797d 100644 --- a/hoomd/BoxResizeUpdaterGPU.cuh +++ b/hoomd/BoxResizeUpdaterGPU.cuh @@ -21,6 +21,7 @@ hipError_t gpu_box_resize_scale(Scalar4* d_pos, hipError_t gpu_box_resize_wrap(const unsigned int N, Scalar4* d_pos, + Scalar4* d_vel, int3* d_image, const BoxDim& new_box, unsigned int block_size); diff --git a/hoomd/Communicator.cc b/hoomd/Communicator.cc index 1b16f475ed..c16423861c 100644 --- a/hoomd/Communicator.cc +++ b/hoomd/Communicator.cc @@ -1719,9 +1719,10 @@ void Communicator::migrateParticles() { detail::pdata_element& p = m_recvbuf[idx]; Scalar4& postype = p.pos; + Scalar4& vel = p.vel; int3& image = p.image; - shifted_box.wrap(postype, image); + shifted_box.wrap(postype, vel, image); } // remove particles that were sent and fill particle data with received particles @@ -2389,11 +2390,14 @@ void Communicator::exchangeGhosts() } // wrap particle positions - if (flags[comm_flag::position]) + if (flags[comm_flag::position] || flags[comm_flag::velocity]) { ArrayHandle h_pos(m_pdata->getPositions(), access_location::host, access_mode::readwrite); + ArrayHandle h_vel(m_pdata->getVelocities(), + access_location::host, + access_mode::readwrite); ArrayHandle h_image(m_pdata->getImages(), access_location::host, access_mode::readwrite); @@ -2403,10 +2407,18 @@ void Communicator::exchangeGhosts() for (unsigned int idx = start_idx; idx < start_idx + m_num_recv_ghosts[dir]; idx++) { Scalar4& pos = h_pos.data[idx]; + int3& img = h_image.data[idx]; // wrap particles received across a global boundary - int3& img = h_image.data[idx]; - shifted_box.wrap(pos, img); + if (flags[comm_flag::velocity]) + { + Scalar4& vel = h_vel.data[idx]; + shifted_box.wrap(pos, vel, img); + } + else + { + shifted_box.wrap(pos, img); + } } } @@ -2961,20 +2973,31 @@ void Communicator::beginUpdateGhosts(uint64_t timestep) MPI_Waitall(2, &m_reqs.front(), &m_stats.front()); } // wrap particle positions (only if copying positions) - if (flags[comm_flag::position]) + if (flags[comm_flag::position] || flags[comm_flag::velocity]) { ArrayHandle h_pos(m_pdata->getPositions(), access_location::host, access_mode::readwrite); + ArrayHandle h_vel(m_pdata->getVelocities(), + access_location::host, + access_mode::readwrite); const BoxDim shifted_box = getShiftedBox(); for (unsigned int idx = start_idx; idx < start_idx + m_num_recv_ghosts[dir]; idx++) { Scalar4& pos = h_pos.data[idx]; + int3 img = make_int3(0, 0, 0); // wrap particles received across a global boundary - int3 img = make_int3(0, 0, 0); - shifted_box.wrap(pos, img); + if (flags[comm_flag::velocity]) + { + Scalar4& vel = h_vel.data[idx]; + shifted_box.wrap(pos, vel, img); + } + else + { + shifted_box.wrap(pos, img); + } } } diff --git a/hoomd/CommunicatorGPU.cu b/hoomd/CommunicatorGPU.cu index 4ea21496f1..bc5037d3d5 100644 --- a/hoomd/CommunicatorGPU.cu +++ b/hoomd/CommunicatorGPU.cu @@ -260,7 +260,7 @@ __global__ void gpu_wrap_particles_kernel(const unsigned int n_recv, return; detail::pdata_element p = d_recv[idx]; - box.wrap(p.pos, p.image); + box.wrap(p.pos, p.vel, p.image); d_recv[idx] = p; } @@ -789,8 +789,10 @@ gpu_pack_kernel(unsigned int n_out, const uint2* d_ghost_idx_adj, const T* in, T __global__ void gpu_pack_wrap_kernel(unsigned int n_out, const uint2* d_ghost_idx_adj, const Scalar4* d_postype, + const Scalar4* d_vel, const int3* d_img, Scalar4* out_pos, + Scalar4* out_vel, int3* out_img, Index3D di, uint3 my_pos, @@ -860,14 +862,24 @@ __global__ void gpu_pack_wrap_kernel(unsigned int n_out, } box.setPeriodic(periodic); + Scalar4 postype = d_postype[idx]; int3 img = make_int3(0, 0, 0); - if (d_img) + if (out_img) + { img = d_img[idx]; - Scalar4 postype = d_postype[idx]; - box.wrap(postype, img, wrap); + } + if (out_vel) + { + Scalar4 vel = d_vel[idx]; + box.wrap(postype, vel, img, wrap); + out_vel[buf_idx] = vel; + } + else + { + box.wrap(postype, img, wrap); + } out_pos[buf_idx] = postype; - if (out_img) { out_img[buf_idx] = img; @@ -928,10 +940,19 @@ void gpu_exchange_ghosts_pack(unsigned int n_out, d_tag, d_tag_sendbuf); } - if (send_pos) + /* + * combined packing pathway for positions , velocities and images + * call gpu_pack_wrap_kernel whenever any of send_pos, send_vel and send_image is true + * for buffers that are not being sent, pass a nullptr to the respective output + */ + if (send_pos || send_vel || send_image) { assert(d_pos); assert(d_pos_sendbuf); + if (send_vel) + { + assert(d_vel); + } if (send_image) { assert(d_img); @@ -944,27 +965,15 @@ void gpu_exchange_ghosts_pack(unsigned int n_out, n_out, d_ghost_idx_adj, d_pos, + d_vel, d_img, d_pos_sendbuf, + send_vel ? d_vel_sendbuf : 0, send_image ? d_img_sendbuf : 0, di, my_pos, box); } - if (send_vel) - { - assert(d_vel); - assert(d_vel_sendbuf); - hipLaunchKernelGGL(gpu_pack_kernel, - dim3(n_blocks), - dim3(block_size), - 0, - 0, - n_out, - d_ghost_idx_adj, - d_vel, - d_vel_sendbuf); - } if (send_charge) { assert(d_charge); diff --git a/hoomd/GSDDumpWriter.cc b/hoomd/GSDDumpWriter.cc index eb3b1da077..ad7d8ad3f5 100644 --- a/hoomd/GSDDumpWriter.cc +++ b/hoomd/GSDDumpWriter.cc @@ -1113,6 +1113,9 @@ void GSDDumpWriter::populateLocalFrame(GSDDumpWriter::GSDFrame& frame, uint64_t ArrayHandle h_postype(m_pdata->getPositions(), access_location::host, access_mode::read); + ArrayHandle h_veltype(m_pdata->getVelocities(), + access_location::host, + access_mode::read); ArrayHandle h_image(m_pdata->getImages(), access_location::host, access_mode::read); if (m_dynamic[gsd_flag::particles_position] || m_nframes == 0) @@ -1132,6 +1135,7 @@ void GSDDumpWriter::populateLocalFrame(GSDDumpWriter::GSDFrame& frame, uint64_t { vec3 position = vec3(h_postype.data[index]) - vec3(m_pdata->getOrigin()); + vec3 velocity = vec3(h_veltype.data[index]); unsigned int type = __scalar_as_int(h_postype.data[index].w); int3 image = make_int3(0, 0, 0); @@ -1140,7 +1144,7 @@ void GSDDumpWriter::populateLocalFrame(GSDDumpWriter::GSDFrame& frame, uint64_t image = h_image.data[index]; } - frame.global_box.wrap(position, image); + frame.global_box.wrap(position, velocity, image); if (m_dynamic[gsd_flag::particles_position] || m_nframes == 0) { diff --git a/hoomd/ParticleData.cc b/hoomd/ParticleData.cc index 9db2fdd373..aa7daa073b 100644 --- a/hoomd/ParticleData.cc +++ b/hoomd/ParticleData.cc @@ -793,6 +793,8 @@ void ParticleData::initializeFromSnapshot(const SnapshotParticleData& snap { unsigned int snap_idx = (unsigned int)(it - snapshot.pos.begin()); + Scalar3 v = vec_to_scalar3(snapshot.vel[snap_idx]); + // if requested, do not initialize constituent particles of bodies if (ignore_bodies && snapshot.body[snap_idx] < MIN_FLOPPY && snapshot.body[snap_idx] != snap_idx) @@ -834,7 +836,7 @@ void ParticleData::initializeFromSnapshot(const SnapshotParticleData& snap // only wrap if the particles is on one of the boundaries uchar3 periodic = make_uchar3(flags.x, flags.y, flags.z); global_box.setPeriodic(periodic); - global_box.wrap(pos, img, flags); + global_box.wrap(pos, v, img, flags); // place particle using actual domain fractions, not global box fraction unsigned int rank @@ -1276,9 +1278,11 @@ template void ParticleData::takeSnapshot(SnapshotParticleData& snapshot.inertia[snap_id] = vec3(inertia_proc[rank][idx]); // make sure the position stored in the snapshot is within the boundaries - Scalar3 tmp = vec_to_scalar3(snapshot.pos[snap_id]); - m_global_box->wrap(tmp, snapshot.image[snap_id]); - snapshot.pos[snap_id] = vec3(tmp); + Scalar3 tmp_p = vec_to_scalar3(snapshot.pos[snap_id]); + Scalar3 tmp_v = vec_to_scalar3(snapshot.vel[snap_id]); + m_global_box->wrap(tmp_p, tmp_v, snapshot.image[snap_id]); + snapshot.pos[snap_id] = vec3(tmp_p); + snapshot.vel[snap_id] = vec3(tmp_v); std::advance(tag_set_it, 1); } @@ -1326,9 +1330,11 @@ template void ParticleData::takeSnapshot(SnapshotParticleData& snapshot.inertia[snap_id] = vec3(h_inertia.data[idx]); // make sure the position stored in the snapshot is within the boundaries - Scalar3 tmp = vec_to_scalar3(snapshot.pos[snap_id]); - m_global_box->wrap(tmp, snapshot.image[snap_id]); - snapshot.pos[snap_id] = vec3(tmp); + Scalar3 tmp_p = vec_to_scalar3(snapshot.pos[snap_id]); + Scalar3 tmp_v = vec_to_scalar3(snapshot.vel[snap_id]); + m_global_box->wrap(tmp_p, tmp_v, snapshot.image[snap_id]); + snapshot.pos[snap_id] = vec3(tmp_p); + snapshot.vel[snap_id] = vec3(tmp_v); } } @@ -1440,7 +1446,6 @@ Scalar3 ParticleData::getPosition(unsigned int tag) const } #endif assert(found); - m_global_box->wrap(result, img); return result; } @@ -1451,10 +1456,24 @@ Scalar3 ParticleData::getVelocity(unsigned int tag) const unsigned int idx = getRTag(tag); bool found = (idx < getN()); Scalar3 result = make_scalar3(0.0, 0.0, 0.0); + Scalar3 pos = make_scalar3(0.0, 0.0, 0.0); + int3 img = make_int3(0, 0, 0); if (found) { ArrayHandle h_vel(m_vel, access_location::host, access_mode::read); result = make_scalar3(h_vel.data[idx].x, h_vel.data[idx].y, h_vel.data[idx].z); + + ArrayHandle h_postype(m_pos, access_location::host, access_mode::read); + pos = make_scalar3(h_postype.data[idx].x, h_postype.data[idx].y, h_postype.data[idx].z); + pos = pos - m_origin; + + ArrayHandle h_img(m_image, access_location::host, access_mode::read); + img = make_int3(h_img.data[idx].x, h_img.data[idx].y, h_img.data[idx].z); + img.x -= m_o_image.x; + img.y -= m_o_image.y; + img.z -= m_o_image.z; + + m_global_box->wrap(pos, result, img); } #ifdef ENABLE_MPI if (m_decomposition) @@ -1509,6 +1528,14 @@ int3 ParticleData::getImage(unsigned int tag) const result = make_int3(h_image.data[idx].x, h_image.data[idx].y, h_image.data[idx].z); pos = make_scalar3(h_postype.data[idx].x, h_postype.data[idx].y, h_postype.data[idx].z); pos = pos - m_origin; + + // correct for origin shift + result.x -= m_o_image.x; + result.y -= m_o_image.y; + result.z -= m_o_image.z; + + // wrap into correct image + m_global_box->wrap(pos, result); } #ifdef ENABLE_MPI if (m_decomposition) @@ -1517,22 +1544,10 @@ int3 ParticleData::getImage(unsigned int tag) const bcast((int&)result.x, owner_rank, m_exec_conf->getMPICommunicator()); bcast((int&)result.y, owner_rank, m_exec_conf->getMPICommunicator()); bcast((int&)result.z, owner_rank, m_exec_conf->getMPICommunicator()); - bcast((Scalar&)pos.x, owner_rank, m_exec_conf->getMPICommunicator()); - bcast((Scalar&)pos.y, owner_rank, m_exec_conf->getMPICommunicator()); - bcast((Scalar&)pos.z, owner_rank, m_exec_conf->getMPICommunicator()); found = true; } #endif assert(found); - - // correct for origin shift - result.x -= m_o_image.x; - result.y -= m_o_image.y; - result.z -= m_o_image.z; - - // wrap into correct image - m_global_box->wrap(pos, result); - return result; } @@ -1823,31 +1838,40 @@ void ParticleData::setPosition(unsigned int tag, const Scalar3& pos, bool move) owner_rank = getOwnerRank(tag); #endif - // load current particle image + // load current particle image and vel int3 img; + Scalar3 tmp_vel; if (ptl_local) { ArrayHandle h_image(m_image, access_location::host, access_mode::read); img = h_image.data[idx]; + + ArrayHandle h_vel(m_vel, access_location::host, access_mode::read); + tmp_vel = make_scalar3(h_vel.data[idx].x, h_vel.data[idx].y, h_vel.data[idx].z); } else { - // if we don't own the particle, we are not going to use image flags + // if we don't own the particle, we are not going to use any data img = make_int3(0, 0, 0); + tmp_vel = make_scalar3(0, 0, 0); } - // wrap into box and update image - m_global_box->wrap(tmp_pos, img); + // wrap into box and update image and vel + m_global_box->wrap(tmp_pos, tmp_vel, img); - // store position and image + // store position, velocity and image if (ptl_local) { ArrayHandle h_pos(m_pos, access_location::host, access_mode::readwrite); + ArrayHandle h_vel(m_vel, access_location::host, access_mode::readwrite); ArrayHandle h_image(m_image, access_location::host, access_mode::readwrite); h_pos.data[idx].x = tmp_pos.x; h_pos.data[idx].y = tmp_pos.y; h_pos.data[idx].z = tmp_pos.z; + h_vel.data[idx].x = tmp_vel.x; + h_vel.data[idx].y = tmp_vel.y; + h_vel.data[idx].z = tmp_vel.z; h_image.data[idx] = img; } @@ -3225,12 +3249,14 @@ void SnapshotParticleData::replicate(unsigned int nx, for (unsigned int i = 0; i < old_size; ++i) { // unwrap position of particle i in old box using image flags - vec3 p = pos[i]; + Scalar3 r = make_scalar3(pos[i].x, pos[i].y, pos[i].z); + Scalar3 v = make_scalar3(vel[i].x, vel[i].y, vel[i].z); int3 img = image[i]; - // need to cast to a scalar and back because the Box is in Scalars, but we might be in a + // need to cast scalar to Real because the Box is in Scalars, but we might be in a // different type - p = vec3(old_box.shift(vec3(p), img)); + old_box.shift(r, v, img); + vec3 p(r); vec3 f = old_box.makeFraction(p); unsigned int j = 0; @@ -3252,10 +3278,10 @@ void SnapshotParticleData::replicate(unsigned int nx, // wrap by multiple box vectors if necessary image[k] = new_box.getImage(q); int3 negimg = make_int3(-image[k].x, -image[k].y, -image[k].z); - q = new_box.shift(q, negimg); + new_box.shift(q, v, negimg); // rewrap using wrap so that rounding is consistent - new_box.wrap(q, image[k]); + new_box.wrap(q, v, image[k]); pos[k] = vec3(q); vel[k] = vel[i]; diff --git a/hoomd/UpdaterRemoveDrift.h b/hoomd/UpdaterRemoveDrift.h index d6eae68c01..cc58b4e7f6 100644 --- a/hoomd/UpdaterRemoveDrift.h +++ b/hoomd/UpdaterRemoveDrift.h @@ -100,6 +100,9 @@ class UpdaterRemoveDrift : public Updater ArrayHandle h_postype(this->m_pdata->getPositions(), access_location::host, access_mode::readwrite); + ArrayHandle h_vel(this->m_pdata->getVelocities(), + access_location::host, + access_mode::readwrite); ArrayHandle h_tag(this->m_pdata->getTags(), access_location::host, access_mode::read); @@ -116,8 +119,9 @@ class UpdaterRemoveDrift : public Updater unsigned int tag_i = h_tag.data[i]; // read in the current position and orientation vec3 postype_i = vec3(h_postype.data[i]) - origin; + vec3 vel = vec3(h_vel.data[i]); int3 tmp_image = make_int3(0, 0, 0); - box.wrap(postype_i, tmp_image); + box.wrap(postype_i, vel, tmp_image); const vec3 dr = postype_i - m_ref_positions[tag_i]; rshift += vec3(box.minImage(vec_to_scalar3(dr))); } @@ -146,7 +150,7 @@ class UpdaterRemoveDrift : public Updater Scalar4 postype_i = h_postype.data[i]; const vec3 r_i = vec3(postype_i); h_postype.data[i] = vec_to_scalar4(r_i - rshift, postype_i.w); - box.wrap(h_postype.data[i], h_image.data[i]); + box.wrap(h_postype.data[i], h_vel.data[i], h_image.data[i]); } } diff --git a/hoomd/md/ForceComposite.cc b/hoomd/md/ForceComposite.cc index a6d700f81f..7be94ee1e0 100644 --- a/hoomd/md/ForceComposite.cc +++ b/hoomd/md/ForceComposite.cc @@ -603,11 +603,13 @@ void ForceComposite::createRigidBodies( quat constituent_orientation = body_orientation * local_orientation; snap.pos[constituent_particle_tag] = constituent_position; + snap.vel[constituent_particle_tag] = snap.vel[particle_tag]; snap.image[constituent_particle_tag] = snap.image[particle_tag]; snap.orientation[constituent_particle_tag] = constituent_orientation; // wrap back into the box global_box.wrap(snap.pos[constituent_particle_tag], + snap.vel[constituent_particle_tag], snap.image[constituent_particle_tag]); // Since the central particle tags here will be [0, n_central_particles), we know @@ -1050,14 +1052,14 @@ void ForceComposite::updateCompositeParticles(uint64_t timestep) angvel_body.z = Scalar(0); } - const vec3 updated_vel = vec3(h_velocity.data[central_idx]) - + rotate(orientation, cross(angvel_body, local_pos)); + vec3 updated_vel = vec3(h_velocity.data[central_idx]) + + rotate(orientation, cross(angvel_body, local_pos)); // this runs before the ForceComputes, // wrap into box, allowing rigid bodies to span multiple images int3 imgi = box.getImage(vec_to_scalar3(updated_pos)); int3 negimgi = make_int3(-imgi.x, -imgi.y, -imgi.z); - updated_pos = global_box.shift(updated_pos, negimgi); + global_box.shift(updated_pos, updated_vel, negimgi); const Scalar mass = h_velocity.data[particle_index].w; h_postype.data[particle_index] diff --git a/hoomd/md/ForceCompositeGPU.cu b/hoomd/md/ForceCompositeGPU.cu index 022d351728..1bb2871f40 100644 --- a/hoomd/md/ForceCompositeGPU.cu +++ b/hoomd/md/ForceCompositeGPU.cu @@ -756,13 +756,13 @@ __global__ void gpu_update_composite_kernel(unsigned int N, angvel_body.z = Scalar(0); } - const vec3 updated_vel = vec3(d_velocity[central_idx]) - + rotate(orientation, cross(angvel_body, local_pos)); + vec3 updated_vel = vec3(d_velocity[central_idx]) + + rotate(orientation, cross(angvel_body, local_pos)); // this runs before the ForceComputes, // wrap into box, allowing rigid bodies to span multiple images int3 imgi = box.getImage(vec_to_scalar3(updated_pos)); int3 negimgi = make_int3(-imgi.x, -imgi.y, -imgi.z); - updated_pos = global_box.shift(updated_pos, negimgi); + global_box.shift(updated_pos, updated_vel, negimgi); unsigned int type = __scalar_as_int(d_postype[idx].w); const Scalar mass = d_velocity[idx].w; diff --git a/hoomd/md/ForceDistanceConstraint.cc b/hoomd/md/ForceDistanceConstraint.cc index 4b05813835..3bdf8267f0 100644 --- a/hoomd/md/ForceDistanceConstraint.cc +++ b/hoomd/md/ForceDistanceConstraint.cc @@ -201,15 +201,16 @@ void ForceDistanceConstraint::fillMatrixVector(uint64_t timestep) vec3 rb(h_pos.data[idx_b]); vec3 rn(ra - rb); - // apply minimum image - rn = box.minImage(rn); - vec3 va(h_vel.data[idx_a]); Scalar ma(h_vel.data[idx_a].w); vec3 vb(h_vel.data[idx_b]); Scalar mb(h_vel.data[idx_b].w); vec3 rndot(va - vb); + + // apply minimum image + box.minImage(rn, rndot); + vec3 qn(rn + rndot * m_deltaT); // fill matrix row diff --git a/hoomd/md/ForceDistanceConstraintGPU.cu b/hoomd/md/ForceDistanceConstraintGPU.cu index 7c796c8f89..bc14afe5ed 100644 --- a/hoomd/md/ForceDistanceConstraintGPU.cu +++ b/hoomd/md/ForceDistanceConstraintGPU.cu @@ -76,9 +76,6 @@ __global__ void gpu_fill_matrix_vector_kernel(unsigned int n_constraint, // constraint separation vec3 rn(vec3(d_pos[idx_na]) - vec3(d_pos[idx_nb])); - // apply minimum image - rn = box.minImage(rn); - // get masses Scalar ma = d_vel[idx_na].w; Scalar mb = d_vel[idx_nb].w; @@ -86,6 +83,9 @@ __global__ void gpu_fill_matrix_vector_kernel(unsigned int n_constraint, // constraint time derivative vec3 rndot(vec3(d_vel[idx_na]) - vec3(d_vel[idx_nb])); + // apply minimum image + box.minImage(rn, rndot); + // constraint separation at t+2*deltaT vec3 qn(rn + deltaT * rndot); diff --git a/hoomd/md/PotentialPairDPDThermo.h b/hoomd/md/PotentialPairDPDThermo.h index c78399ac04..3099600185 100644 --- a/hoomd/md/PotentialPairDPDThermo.h +++ b/hoomd/md/PotentialPairDPDThermo.h @@ -184,7 +184,7 @@ template void PotentialPairDPDThermo::computeForces( assert(typej < this->m_pdata->getNTypes()); // apply periodic boundary conditions - dx = box.minImage(dx); + box.minImage(dx, dv); // calculate r_ij squared (FLOPS: 5) Scalar rsq = dot(dx, dx); diff --git a/hoomd/md/PotentialPairDPDThermoGPU.cuh b/hoomd/md/PotentialPairDPDThermoGPU.cuh index 8fcee9a0b6..7e9f2c7e0d 100644 --- a/hoomd/md/PotentialPairDPDThermoGPU.cuh +++ b/hoomd/md/PotentialPairDPDThermoGPU.cuh @@ -247,15 +247,15 @@ __global__ void gpu_compute_dpd_forces_kernel(Scalar4* d_force, // calculate dr (with periodic boundary conditions) (FLOPS: 3) Scalar3 dx = posi - posj; + // calculate dv (FLOPS: 3) + Scalar3 dv = veli - velj; + // apply periodic boundary conditions: (FLOPS 12) - dx = box.minImage(dx); + box.minImage(dx, dv); // calculate r squared (FLOPS: 5) Scalar rsq = dot(dx, dx); - // calculate dv (FLOPS: 3) - Scalar3 dv = veli - velj; - Scalar rdotv = dot(dx, dv); // access the per type pair parameters diff --git a/hoomd/md/TwoStepBD.cc b/hoomd/md/TwoStepBD.cc index 15a749a212..cd4e4dac03 100644 --- a/hoomd/md/TwoStepBD.cc +++ b/hoomd/md/TwoStepBD.cc @@ -128,7 +128,7 @@ void TwoStepBD::integrateStepOne(uint64_t timestep) // particles may have been moved slightly outside the box by the above steps, wrap them back // into place - box.wrap(h_pos.data[j], h_image.data[j]); + box.wrap(h_pos.data[j], h_vel.data[j], h_image.data[j]); if (m_noiseless_t) { diff --git a/hoomd/md/TwoStepBDGPU.cu b/hoomd/md/TwoStepBDGPU.cu index d679cbb054..bf5434a50b 100644 --- a/hoomd/md/TwoStepBDGPU.cu +++ b/hoomd/md/TwoStepBDGPU.cu @@ -161,7 +161,7 @@ __global__ void gpu_brownian_step_one_kernel(Scalar4* d_pos, // particles may have been moved slightly outside the box by the above steps, wrap them back // into place - box.wrap(postype, image); + box.wrap(postype, vel, image); if (d_noiseless_t) { diff --git a/hoomd/md/TwoStepConstantPressure.cc b/hoomd/md/TwoStepConstantPressure.cc index 2d925b3a0d..c81a237914 100644 --- a/hoomd/md/TwoStepConstantPressure.cc +++ b/hoomd/md/TwoStepConstantPressure.cc @@ -351,13 +351,16 @@ void TwoStepConstantPressure::integrateStepOne(uint64_t timestep) ArrayHandle h_pos(m_pdata->getPositions(), access_location::host, access_mode::readwrite); + ArrayHandle h_vel(m_pdata->getVelocities(), + access_location::host, + access_mode::readwrite); ArrayHandle h_image(m_pdata->getImages(), access_location::host, access_mode::readwrite); // Wrap particles for (unsigned int j = 0; j < m_pdata->getN(); j++) - box.wrap(h_pos.data[j], h_image.data[j]); + box.wrap(h_pos.data[j], h_vel.data[j], h_image.data[j]); } // Integration of angular degrees of freedom using symplectic and diff --git a/hoomd/md/TwoStepConstantPressureGPU.cc b/hoomd/md/TwoStepConstantPressureGPU.cc index 65f9558701..4a31cc602b 100644 --- a/hoomd/md/TwoStepConstantPressureGPU.cc +++ b/hoomd/md/TwoStepConstantPressureGPU.cc @@ -199,6 +199,9 @@ void TwoStepConstantPressureGPU::integrateStepOne(uint64_t timestep) ArrayHandle d_pos(m_pdata->getPositions(), access_location::device, access_mode::readwrite); + ArrayHandle d_vel(m_pdata->getVelocities(), + access_location::device, + access_mode::readwrite); ArrayHandle d_image(m_pdata->getImages(), access_location::device, access_mode::readwrite); @@ -208,6 +211,7 @@ void TwoStepConstantPressureGPU::integrateStepOne(uint64_t timestep) kernel::gpu_npt_rescale_wrap(m_pdata->getN(), d_pos.data, + d_vel.data, d_image.data, box, m_tuner_wrap->getParam()[0]); diff --git a/hoomd/md/TwoStepConstantPressureGPU.cu b/hoomd/md/TwoStepConstantPressureGPU.cu index 3622d9d400..2f3e273f9f 100644 --- a/hoomd/md/TwoStepConstantPressureGPU.cu +++ b/hoomd/md/TwoStepConstantPressureGPU.cu @@ -177,13 +177,17 @@ hipError_t gpu_npt_rescale_step_one(Scalar4* d_pos, /*! \param N number of particles in the system \param d_pos array of particle positions + \param d_vel array of particle velocities \param d_image array of particle images \param box The new box the particles where the particles now reside Wrap particle positions for all particles in the box */ -__global__ void -gpu_npt_mtk_wrap_kernel(const unsigned int nwork, Scalar4* d_pos, int3* d_image, BoxDim box) +__global__ void gpu_npt_mtk_wrap_kernel(const unsigned int nwork, + Scalar4* d_pos, + Scalar4* d_vel, + int3* d_image, + BoxDim box) { // determine which particle this thread works on int idx = blockIdx.x * blockDim.x + threadIdx.x; @@ -194,21 +198,25 @@ gpu_npt_mtk_wrap_kernel(const unsigned int nwork, Scalar4* d_pos, int3* d_image, // fetch particle position Scalar4 postype = d_pos[idx]; Scalar3 pos = make_scalar3(postype.x, postype.y, postype.z); + Scalar4 velmass = d_vel[idx]; + Scalar3 vel = make_scalar3(velmass.x, velmass.y, velmass.z); // read in the image flags int3 image = d_image[idx]; // fix periodic boundary conditions - box.wrap(pos, image); + box.wrap(pos, vel, image); // write out the results d_pos[idx] = make_scalar4(pos.x, pos.y, pos.z, postype.w); + d_vel[idx] = make_scalar4(vel.x, vel.y, vel.z, velmass.w); d_image[idx] = image; } } /*! \param N number of particles in the system \param d_pos array of particle positions + \param d_vel array of particle velocities \param d_image array of particle images \param box The new box the particles where the particles now reside @@ -216,6 +224,7 @@ gpu_npt_mtk_wrap_kernel(const unsigned int nwork, Scalar4* d_pos, int3* d_image, */ hipError_t gpu_npt_rescale_wrap(const unsigned int N, Scalar4* d_pos, + Scalar4* d_vel, int3* d_image, const BoxDim& box, const unsigned int block_size) @@ -241,6 +250,7 @@ hipError_t gpu_npt_rescale_wrap(const unsigned int N, 0, nwork, d_pos, + d_vel, d_image, box); diff --git a/hoomd/md/TwoStepConstantPressureGPU.cuh b/hoomd/md/TwoStepConstantPressureGPU.cuh index f80cf82bf0..b67960f79d 100644 --- a/hoomd/md/TwoStepConstantPressureGPU.cuh +++ b/hoomd/md/TwoStepConstantPressureGPU.cuh @@ -37,6 +37,7 @@ hipError_t gpu_npt_rescale_step_one(Scalar4* d_pos, //! Kernel driver for wrapping particles back in the box (part of first step) hipError_t gpu_npt_rescale_wrap(const unsigned int N, Scalar4* d_pos, + Scalar4* d_vel, int3* d_image, const BoxDim& box, const unsigned int block_size); diff --git a/hoomd/md/TwoStepConstantVolume.cc b/hoomd/md/TwoStepConstantVolume.cc index 6e7a30ad47..3c72112dce 100644 --- a/hoomd/md/TwoStepConstantVolume.cc +++ b/hoomd/md/TwoStepConstantVolume.cc @@ -75,7 +75,7 @@ void hoomd::md::TwoStepConstantVolume::integrateStepOne(uint64_t timestep) { unsigned int j = m_group->getMemberIndex(group_idx); // wrap the particles around the box - box.wrap(h_pos.data[j], h_image.data[j]); + box.wrap(h_pos.data[j], h_vel.data[j], h_image.data[j]); } } diff --git a/hoomd/md/TwoStepConstantVolumeGPU.cu b/hoomd/md/TwoStepConstantVolumeGPU.cu index 10ffddd0d4..99a6ad82c3 100644 --- a/hoomd/md/TwoStepConstantVolumeGPU.cu +++ b/hoomd/md/TwoStepConstantVolumeGPU.cu @@ -77,7 +77,7 @@ __global__ void gpu_nvt_rescale_step_one_kernel(Scalar4* d_pos, int3 image = d_image[idx]; // time to fix the periodic boundary conditions - box.wrap(pos, image); + box.wrap(pos, vel, image); // write out the results d_pos[idx] = make_scalar4(pos.x, pos.y, pos.z, postype.w); diff --git a/hoomd/md/TwoStepLangevin.cc b/hoomd/md/TwoStepLangevin.cc index 3cc0d5d6e9..b1a22e71e1 100644 --- a/hoomd/md/TwoStepLangevin.cc +++ b/hoomd/md/TwoStepLangevin.cc @@ -73,7 +73,7 @@ void TwoStepLangevin::integrateStepOne(uint64_t timestep) h_pos.data[j].z += dz; // particles may have been moved slightly outside the box by the above steps, wrap them back // into place - box.wrap(h_pos.data[j], h_image.data[j]); + box.wrap(h_pos.data[j], h_vel.data[j], h_image.data[j]); h_vel.data[j].x += Scalar(1.0 / 2.0) * h_accel.data[j].x * m_deltaT; h_vel.data[j].y += Scalar(1.0 / 2.0) * h_accel.data[j].y * m_deltaT; diff --git a/hoomd/md/TwoStepNVEGPU.cu b/hoomd/md/TwoStepNVEGPU.cu index 96990cdc17..92eb3e3cd0 100644 --- a/hoomd/md/TwoStepNVEGPU.cu +++ b/hoomd/md/TwoStepNVEGPU.cu @@ -100,7 +100,7 @@ __global__ void gpu_nve_step_one_kernel(Scalar4* d_pos, int3 image = d_image[idx]; // fix the periodic boundary conditions (FLOPS: 15) - box.wrap(pos, image); + box.wrap(pos, vel, image); // write out the results (MEM_TRANSFER: 48 bytes) d_pos[idx] = make_scalar4(pos.x, pos.y, pos.z, postype.w); diff --git a/hoomd/md/TwoStepRATTLEBD.h b/hoomd/md/TwoStepRATTLEBD.h index d3d4a761f0..7942e2bf45 100644 --- a/hoomd/md/TwoStepRATTLEBD.h +++ b/hoomd/md/TwoStepRATTLEBD.h @@ -151,10 +151,10 @@ template void TwoStepRATTLEBD::integrateStepOne(uint64 const Scalar currentTemp = m_T->operator()(timestep); const GPUArray& net_force = m_pdata->getNetForce(); - ArrayHandle h_vel(m_pdata->getVelocities(), + ArrayHandle h_pos(m_pdata->getPositions(), access_location::host, access_mode::readwrite); - ArrayHandle h_pos(m_pdata->getPositions(), + ArrayHandle h_vel(m_pdata->getVelocities(), access_location::host, access_mode::readwrite); ArrayHandle h_image(m_pdata->getImages(), access_location::host, access_mode::readwrite); @@ -297,7 +297,7 @@ template void TwoStepRATTLEBD::integrateStepOne(uint64 // particles may have been moved slightly outside the box by the above steps, wrap them back // into place - box.wrap(h_pos.data[j], h_image.data[j]); + box.wrap(h_pos.data[j], h_vel.data[j], h_image.data[j]); // rotational random force and orientation quaternion updates if (m_aniso) diff --git a/hoomd/md/TwoStepRATTLEBDGPU.cuh b/hoomd/md/TwoStepRATTLEBDGPU.cuh index a6225e8239..ae5cecdeb3 100644 --- a/hoomd/md/TwoStepRATTLEBDGPU.cuh +++ b/hoomd/md/TwoStepRATTLEBDGPU.cuh @@ -245,7 +245,7 @@ __global__ void gpu_rattle_brownian_step_one_kernel(Scalar4* d_pos, postype.z += dz; // particles may have been moved slightly outside the box by the above steps, wrap them back // into place - box.wrap(postype, image); + box.wrap(postype, vel, image); // write out data d_pos[idx] = postype; diff --git a/hoomd/md/TwoStepRATTLELangevin.h b/hoomd/md/TwoStepRATTLELangevin.h index 1c032a6343..f155ecfbd2 100644 --- a/hoomd/md/TwoStepRATTLELangevin.h +++ b/hoomd/md/TwoStepRATTLELangevin.h @@ -210,7 +210,7 @@ template void TwoStepRATTLELangevin::integrateStepOne( // particles may have been moved slightly outside the box by the above steps, wrap them back // into place - box.wrap(h_pos.data[j], h_image.data[j]); + box.wrap(h_pos.data[j], h_vel.data[j], h_image.data[j]); } if (m_aniso) diff --git a/hoomd/md/TwoStepRATTLENVE.h b/hoomd/md/TwoStepRATTLENVE.h index 6d5b6fed37..205c66d08c 100644 --- a/hoomd/md/TwoStepRATTLENVE.h +++ b/hoomd/md/TwoStepRATTLENVE.h @@ -232,7 +232,7 @@ template void TwoStepRATTLENVE::integrateStepOne(uint6 for (unsigned int group_idx = 0; group_idx < group_size; group_idx++) { unsigned int j = m_group->getMemberIndex(group_idx); - box.wrap(h_pos.data[j], h_image.data[j]); + box.wrap(h_pos.data[j], h_vel.data[j], h_image.data[j]); } // Integration of angular degrees of freedom using symplectic and diff --git a/hoomd/md/TwoStepRATTLENVEGPU.cu b/hoomd/md/TwoStepRATTLENVEGPU.cu index 10ce40d94b..ac39370dd6 100644 --- a/hoomd/md/TwoStepRATTLENVEGPU.cu +++ b/hoomd/md/TwoStepRATTLENVEGPU.cu @@ -97,7 +97,7 @@ __global__ void gpu_rattle_nve_step_one_kernel(Scalar4* d_pos, int3 image = d_image[idx]; // fix the periodic boundary conditions (FLOPS: 15) - box.wrap(pos, image); + box.wrap(pos, vel, image); // write out the results (MEM_TRANSFER: 48 bytes) d_pos[idx] = make_scalar4(pos.x, pos.y, pos.z, postype.w); diff --git a/hoomd/mpcd/BounceBackNVE.h b/hoomd/mpcd/BounceBackNVE.h index 551e32aa85..c99c4a1ab2 100644 --- a/hoomd/mpcd/BounceBackNVE.h +++ b/hoomd/mpcd/BounceBackNVE.h @@ -135,7 +135,7 @@ template void BounceBackNVE::integrateStepOne(uint64_t } while (dt_remain > 0 && collide); // wrap final position - box.wrap(pos, h_image.data[pid]); + box.wrap(pos, vel, h_image.data[pid]); // write position and velocity back out h_pos.data[pid] = make_scalar4(pos.x, pos.y, pos.z, type); diff --git a/hoomd/mpcd/BounceBackNVEGPU.cuh b/hoomd/mpcd/BounceBackNVEGPU.cuh index 33174f7241..136ffac1d9 100644 --- a/hoomd/mpcd/BounceBackNVEGPU.cuh +++ b/hoomd/mpcd/BounceBackNVEGPU.cuh @@ -127,7 +127,7 @@ __global__ void nve_bounce_step_one(Scalar4* d_pos, // wrap final position int3 img = d_image[pid]; - box.wrap(pos, img); + box.wrap(pos, vel, img); // write position and velocity back out d_pos[pid] = make_scalar4(pos.x, pos.y, pos.z, type); diff --git a/hoomd/mpcd/BounceBackStreamingMethod.h b/hoomd/mpcd/BounceBackStreamingMethod.h index 894280c173..91de886bef 100644 --- a/hoomd/mpcd/BounceBackStreamingMethod.h +++ b/hoomd/mpcd/BounceBackStreamingMethod.h @@ -148,7 +148,7 @@ void BounceBackStreamingMethod::stream(uint64_t timestep) // wrap and update the position int3 image = make_int3(0, 0, 0); - box.wrap(pos, image); + box.wrap(pos, vel, image); h_pos.data[cur_p] = make_scalar4(pos.x, pos.y, pos.z, __int_as_scalar(type)); h_vel.data[cur_p] diff --git a/hoomd/mpcd/BounceBackStreamingMethodGPU.cuh b/hoomd/mpcd/BounceBackStreamingMethodGPU.cuh index 71500f51ea..4600f14d60 100644 --- a/hoomd/mpcd/BounceBackStreamingMethodGPU.cuh +++ b/hoomd/mpcd/BounceBackStreamingMethodGPU.cuh @@ -114,7 +114,7 @@ __global__ void confined_stream(Scalar4* d_pos, // wrap and update the position int3 image = make_int3(0, 0, 0); - box.wrap(pos, image); + box.wrap(pos, vel, image); d_pos[idx] = make_scalar4(pos.x, pos.y, pos.z, __int_as_scalar(type)); d_vel[idx] = make_scalar4(vel.x, vel.y, vel.z, __int_as_scalar(mpcd::detail::NO_CELL)); diff --git a/hoomd/mpcd/Communicator.cc b/hoomd/mpcd/Communicator.cc index f05f61fcee..b16030cd3d 100644 --- a/hoomd/mpcd/Communicator.cc +++ b/hoomd/mpcd/Communicator.cc @@ -415,9 +415,10 @@ void mpcd::Communicator::migrateParticles(uint64_t timestep) { mpcd::detail::pdata_element& p = h_recvbuf.data[idx]; Scalar4& postype = p.pos; + Scalar4& vel = p.vel; int3 image = make_int3(0, 0, 0); - wrap_box.wrap(postype, image); + wrap_box.wrap(postype, vel, image); } } diff --git a/hoomd/mpcd/CommunicatorGPU.cu b/hoomd/mpcd/CommunicatorGPU.cu index a8f98f9fdd..e419026697 100644 --- a/hoomd/mpcd/CommunicatorGPU.cu +++ b/hoomd/mpcd/CommunicatorGPU.cu @@ -289,7 +289,7 @@ struct wrap_particle_op { mpcd::detail::pdata_element ret = p; int3 image = make_int3(0, 0, 0); - box.wrap(ret.pos, image); + box.wrap(ret.pos, ret.vel, image); return ret; } }; diff --git a/hoomd/mpcd/ParticleData.cc b/hoomd/mpcd/ParticleData.cc index 47cf1acd9b..32624a8ec4 100644 --- a/hoomd/mpcd/ParticleData.cc +++ b/hoomd/mpcd/ParticleData.cc @@ -159,6 +159,8 @@ void mpcd::ParticleData::initializeFromSnapshot(const mpcd::ParticleDataSnapshot { unsigned int snap_idx = (unsigned int)(it - snapshot.position.begin()); + Scalar3 v = vec_to_scalar3(snapshot.velocity[snap_idx]); + // determine domain the particle is placed into Scalar3 pos = vec_to_scalar3(*it); Scalar3 f = global_box->makeFraction(pos); @@ -189,7 +191,7 @@ void mpcd::ParticleData::initializeFromSnapshot(const mpcd::ParticleDataSnapshot auto wrapping_box = *global_box; wrapping_box.setPeriodic(periodic); int3 img = make_int3(0, 0, 0); - wrapping_box.wrap(pos, img, flags); + wrapping_box.wrap(pos, v, img, flags); // place particle into the computational domain unsigned int rank @@ -523,20 +525,18 @@ void mpcd::ParticleData::takeSnapshot(mpcd::ParticleDataSnapshot& snapshot, { const auto particle = particles[idx]; - // wrapped position + // wrapped position and velocity const Scalar4 postype = particle.pos; Scalar3 pos = make_scalar3(postype.x, postype.y, postype.z); + const Scalar4 velcell = particle.vel; + Scalar3 vel = make_scalar3(velcell.x, velcell.y, velcell.z); int3 img = make_int3(0, 0, 0); - global_box->wrap(pos, img); + global_box->wrap(pos, vel, img); snapshot.position[idx] = vec3(pos); + snapshot.velocity[idx] = vec3(vel); // typeid snapshot.type[idx] = __scalar_as_int(postype.w); - - // velocity - const Scalar4 velcell = particle.vel; - const Scalar3 vel = make_scalar3(velcell.x, velcell.y, velcell.z); - snapshot.velocity[idx] = vec3(vel); } } } @@ -561,20 +561,18 @@ void mpcd::ParticleData::takeSnapshot(mpcd::ParticleDataSnapshot& snapshot, { const unsigned int pidx = tag_index[idx].second; - // wrapped position + // wrapped position and velocity const Scalar4 postype = h_pos.data[pidx]; Scalar3 pos = make_scalar3(postype.x, postype.y, postype.z); + const Scalar4 velcell = h_vel.data[pidx]; + Scalar3 vel = make_scalar3(velcell.x, velcell.y, velcell.z); int3 img = make_int3(0, 0, 0); - global_box->wrap(pos, img); + global_box->wrap(pos, vel, img); snapshot.position[idx] = vec3(pos); + snapshot.velocity[idx] = vec3(vel); // typeid snapshot.type[idx] = __scalar_as_int(postype.w); - - // velocity - const Scalar4 velcell = h_vel.data[pidx]; - const Scalar3 vel = make_scalar3(velcell.x, velcell.y, velcell.z); - snapshot.velocity[idx] = vec3(vel); } } diff --git a/hoomd/mpcd/ParticleDataSnapshot.cc b/hoomd/mpcd/ParticleDataSnapshot.cc index e68fa81960..6036c794b5 100644 --- a/hoomd/mpcd/ParticleDataSnapshot.cc +++ b/hoomd/mpcd/ParticleDataSnapshot.cc @@ -125,6 +125,7 @@ void mpcd::ParticleDataSnapshot::replicate(unsigned int nx, // unwrap position of particle i in old box using image flags vec3 p = position[i]; vec3 f = old_box.makeFraction(p); + Scalar3 v = vec_to_scalar3(velocity[i]); unsigned int j = 0; for (unsigned int l = 0; l < nx; ++l) @@ -144,10 +145,10 @@ void mpcd::ParticleDataSnapshot::replicate(unsigned int nx, // coordinates in new box Scalar3 q = new_box.makeCoordinates(f_new); int3 image = make_int3(0, 0, 0); - new_box.wrap(q, image); + new_box.wrap(q, v, image); position[k] = vec3(q); - velocity[k] = velocity[i]; + velocity[k] = vec3(v); type[k] = type[i]; ++j; } // n diff --git a/hoomd/test/test_gridshift_correct.cc b/hoomd/test/test_gridshift_correct.cc index 268316f249..1dd3bb91fa 100644 --- a/hoomd/test/test_gridshift_correct.cc +++ b/hoomd/test/test_gridshift_correct.cc @@ -29,6 +29,11 @@ UP_TEST(ParticleDataGridShiftGetMethods) std::shared_ptr sysdef(new SystemDefinition(3, BoxDim(10.0), 4)); std::shared_ptr pdata = sysdef->getParticleData(); BoxDim box = pdata->getBox(); + + // set a deformation on the box + const Scalar3 L_rate = make_scalar3(.1, .15, -0.2); + box.setLDeformationRate(L_rate); + { ArrayHandle h_pos(pdata->getPositions(), access_location::host, @@ -39,6 +44,8 @@ UP_TEST(ParticleDataGridShiftGetMethods) h_pos.data[0].x = h_pos.data[0].y = h_pos.data[0].z = 0.0; h_pos.data[1].x = h_pos.data[1].y = h_pos.data[1].z = 1.0; + h_vel.data[0].x = h_vel.data[0].y = h_vel.data[0].z = 1.0; + h_vel.data[1].x = h_vel.data[1].y = h_vel.data[1].z = -1.0; } // compute a shift and apply it to all particles, and origin @@ -48,6 +55,9 @@ UP_TEST(ParticleDataGridShiftGetMethods) ArrayHandle h_pos(pdata->getPositions(), access_location::host, access_mode::readwrite); + ArrayHandle h_vel(pdata->getVelocities(), + access_location::host, + access_mode::readwrite); ArrayHandle h_img(pdata->getImages(), access_location::host, access_mode::readwrite); for (unsigned int i = 0; i < pdata->getN(); i++) @@ -57,11 +67,11 @@ UP_TEST(ParticleDataGridShiftGetMethods) vec3 r_i = vec3(pos_i); // translation from local to global coordinates r_i += vec3(shift); h_pos.data[i] = vec_to_scalar4(r_i, pos_i.w); - box.wrap(h_pos.data[i], h_img.data[i]); + box.wrap(h_pos.data[i], h_vel.data[i], h_img.data[i]); } } - // check that the particle positions are still the original ones + // check that the particle positions, velocities and images are still the original ones Scalar3 pos = pdata->getPosition(0); MY_CHECK_SMALL(Scalar(pos.x - 0.0), tol_small); MY_CHECK_SMALL(Scalar(pos.y - 0.0), tol_small); @@ -71,11 +81,20 @@ UP_TEST(ParticleDataGridShiftGetMethods) MY_CHECK_SMALL(Scalar(pos.y - 1.0), tol_small); MY_CHECK_SMALL(Scalar(pos.z - 1.0), tol_small); + Scalar3 vel = pdata->getVelocity(0); + MY_CHECK_SMALL(Scalar(vel.x - 1.0), tol_small); + MY_CHECK_SMALL(Scalar(vel.y - 1.0), tol_small); + MY_CHECK_SMALL(Scalar(vel.z - 1.0), tol_small); + vel = pdata->getVelocity(1); + MY_CHECK_SMALL(Scalar(vel.x - (-1.0)), tol_small); + MY_CHECK_SMALL(Scalar(vel.y - (-1.0)), tol_small); + MY_CHECK_SMALL(Scalar(vel.z - (-1.0)), tol_small); + int3 pimg = pdata->getImage(0); UP_ASSERT_EQUAL(pimg.x, 0); UP_ASSERT_EQUAL(pimg.y, 0); UP_ASSERT_EQUAL(pimg.z, 0); - pimg = pdata->getImage(0); + pimg = pdata->getImage(1); UP_ASSERT_EQUAL(pimg.x, 0); UP_ASSERT_EQUAL(pimg.y, 0); UP_ASSERT_EQUAL(pimg.z, 0); @@ -87,6 +106,9 @@ UP_TEST(ParticleDataGridShiftGetMethods) ArrayHandle h_pos(pdata->getPositions(), access_location::host, access_mode::readwrite); + ArrayHandle h_vel(pdata->getVelocities(), + access_location::host, + access_mode::readwrite); ArrayHandle h_img(pdata->getImages(), access_location::host, access_mode::readwrite); for (unsigned int i = 0; i < pdata->getN(); i++) @@ -96,7 +118,7 @@ UP_TEST(ParticleDataGridShiftGetMethods) vec3 r_i = vec3(pos_i); // translation from local to global coordinates r_i += vec3(shift_img); h_pos.data[i] = vec_to_scalar4(r_i, pos_i.w); - box.wrap(h_pos.data[i], h_img.data[i]); + box.wrap(h_pos.data[i], h_vel.data[i], h_img.data[i]); } } @@ -110,11 +132,22 @@ UP_TEST(ParticleDataGridShiftGetMethods) MY_CHECK_SMALL(Scalar(pos.y - 1.0), tol_small); MY_CHECK_SMALL(Scalar(pos.z - 1.0), tol_small); + // check that the original velocities will change by the deformation rate, due to shift + vel = pdata->getVelocity(0); + MY_CHECK_SMALL(Scalar(vel.x - 1.0 + L_rate.x), tol_small); + MY_CHECK_SMALL(Scalar(vel.y - 1.0 + L_rate.y), tol_small); + MY_CHECK_SMALL(Scalar(vel.z - 1.0 + L_rate.z), tol_small); + vel = pdata->getVelocity(1); + MY_CHECK_SMALL(Scalar(vel.x - (-1.0) + L_rate.x), tol_small); + MY_CHECK_SMALL(Scalar(vel.y - (-1.0) + L_rate.y), tol_small); + MY_CHECK_SMALL(Scalar(vel.z - (-1.0) + L_rate.z), tol_small); + + // check that the particle images are still the original ones pimg = pdata->getImage(0); UP_ASSERT_EQUAL(pimg.x, 0); UP_ASSERT_EQUAL(pimg.y, 0); UP_ASSERT_EQUAL(pimg.z, 0); - pimg = pdata->getImage(0); + pimg = pdata->getImage(1); UP_ASSERT_EQUAL(pimg.x, 0); UP_ASSERT_EQUAL(pimg.y, 0); UP_ASSERT_EQUAL(pimg.z, 0); @@ -126,6 +159,11 @@ UP_TEST(ParticleDataGridShiftSetMethods) std::shared_ptr sysdef(new SystemDefinition(3, BoxDim(10.0), 4)); std::shared_ptr pdata = sysdef->getParticleData(); BoxDim box = pdata->getBox(); + + // set a deformation on the box + const Scalar3 L_rate = make_scalar3(.1, .15, -0.2); + box.setLDeformationRate(L_rate); + { ArrayHandle h_pos(pdata->getPositions(), access_location::host, @@ -136,6 +174,8 @@ UP_TEST(ParticleDataGridShiftSetMethods) h_pos.data[0].x = h_pos.data[0].y = h_pos.data[0].z = 0.0; h_pos.data[1].x = h_pos.data[1].y = h_pos.data[1].z = 1.0; + h_vel.data[0].x = h_vel.data[0].y = h_vel.data[0].z = 1.0; + h_vel.data[1].x = h_vel.data[1].y = h_vel.data[1].z = -1.0; } // compute a shift that will shift the image of the box @@ -145,6 +185,9 @@ UP_TEST(ParticleDataGridShiftSetMethods) ArrayHandle h_pos(pdata->getPositions(), access_location::host, access_mode::readwrite); + ArrayHandle h_vel(pdata->getVelocities(), + access_location::host, + access_mode::readwrite); ArrayHandle h_img(pdata->getImages(), access_location::host, access_mode::readwrite); for (unsigned int i = 0; i < pdata->getN(); i++) @@ -154,7 +197,7 @@ UP_TEST(ParticleDataGridShiftSetMethods) vec3 r_i = vec3(pos_i); // translation from local to global coordinates r_i += vec3(shift_img); h_pos.data[i] = vec_to_scalar4(r_i, pos_i.w); - box.wrap(h_pos.data[i], h_img.data[i]); + box.wrap(h_pos.data[i], h_vel.data[i], h_img.data[i]); } } @@ -168,6 +211,17 @@ UP_TEST(ParticleDataGridShiftSetMethods) MY_CHECK_SMALL(Scalar(pos.y - 1.0), tol_small); MY_CHECK_SMALL(Scalar(pos.z - 1.0), tol_small); + // check that the original velocities will change by the deformation rate, due to shift + Scalar3 vel = pdata->getVelocity(0); + MY_CHECK_SMALL(Scalar(vel.x - 1.0 + L_rate.x), tol_small); + MY_CHECK_SMALL(Scalar(vel.y - 1.0 + L_rate.y), tol_small); + MY_CHECK_SMALL(Scalar(vel.z - 1.0 + L_rate.z), tol_small); + vel = pdata->getVelocity(1); + MY_CHECK_SMALL(Scalar(vel.x - (-1.0) + L_rate.x), tol_small); + MY_CHECK_SMALL(Scalar(vel.y - (-1.0) + L_rate.y), tol_small); + MY_CHECK_SMALL(Scalar(vel.z - (-1.0) + L_rate.z), tol_small); + + // check that the images are still the original ones int3 pimg = pdata->getImage(0); UP_ASSERT_EQUAL(pimg.x, 0); UP_ASSERT_EQUAL(pimg.y, 0); @@ -193,7 +247,23 @@ UP_TEST(ParticleDataGridShiftSetMethods) MY_CHECK_SMALL(ret_pos1.y - new_pos1.y, tol_small); MY_CHECK_SMALL(ret_pos1.z - new_pos1.z, tol_small); - // OK, now do the same with the images + // OK, now do the same with the velocities + Scalar3 new_vel0 = make_scalar3(0.5, 1.0, -2.0); + pdata->setVelocity(0, new_vel0); + Scalar3 new_vel1 = make_scalar3(5.0, 0.1, -0.2); + pdata->setVelocity(1, new_vel1); + + Scalar3 ret_vel0 = pdata->getVelocity(0); + MY_CHECK_SMALL(ret_vel0.x - new_vel0.x, tol_small); + MY_CHECK_SMALL(ret_vel0.y - new_vel0.y, tol_small); + MY_CHECK_SMALL(ret_vel0.z - new_vel0.z, tol_small); + + Scalar3 ret_vel1 = pdata->getVelocity(1); + MY_CHECK_SMALL(ret_vel1.x - new_vel1.x, tol_small); + MY_CHECK_SMALL(ret_vel1.y - new_vel1.y, tol_small); + MY_CHECK_SMALL(ret_vel1.z - new_vel1.z, tol_small); + + // and then the images int3 new_img0 = make_int3(1, -5, 7); pdata->setImage(0, new_img0); int3 new_img1 = make_int3(4, 1, 10); diff --git a/hoomd/test/test_pdata.cc b/hoomd/test/test_pdata.cc index 9823d70e4c..b3789c4c6e 100644 --- a/hoomd/test/test_pdata.cc +++ b/hoomd/test/test_pdata.cc @@ -415,6 +415,302 @@ UP_TEST(BoxDim_triclinic_test) UP_ASSERT_EQUAL(img.z, 0); } +//! Test box deformation methods +// define a helper function +static void deforming_box_test(const Scalar3& L, + const Scalar3& tilt, + const Scalar3& L_rate, + const Scalar3& tilt_rate, + const Scalar3& pos, + const Scalar3& vel, + const Scalar3& wrap_pos, + const Scalar3& wrap_vel, + const Scalar3& minImage_pos, + const Scalar3& minImage_vel) + { + // make box + BoxDim b(L); + b.setTiltFactors(tilt.x, tilt.y, tilt.z); + b.setLDeformationRate(L_rate); + b.setTiltDeformationRates(tilt_rate.x, tilt_rate.y, tilt_rate.z); + + Scalar tol = Scalar(1e-4); + + // for wrap method - check positions and velocities + Scalar3 test_pos = pos; + Scalar3 test_vel = vel; + int3 img = make_int3(0, 0, 0); + + b.wrap(test_pos, test_vel, img); + + MY_CHECK_CLOSE(test_pos.x, wrap_pos.x, tol); + MY_CHECK_CLOSE(test_pos.y, wrap_pos.y, tol); + MY_CHECK_CLOSE(test_pos.z, wrap_pos.z, tol); + MY_CHECK_CLOSE(test_vel.x, wrap_vel.x, tol); + MY_CHECK_CLOSE(test_vel.y, wrap_vel.y, tol); + MY_CHECK_CLOSE(test_vel.z, wrap_vel.z, tol); + + // for minImage method - check positions and velocities + test_pos = pos; + test_vel = vel; + + b.minImage(test_pos, test_vel); + + MY_CHECK_CLOSE(test_pos.x, minImage_pos.x, tol); + MY_CHECK_CLOSE(test_pos.y, minImage_pos.y, tol); + MY_CHECK_CLOSE(test_pos.z, minImage_pos.z, tol); + MY_CHECK_CLOSE(test_vel.x, minImage_vel.x, tol); + MY_CHECK_CLOSE(test_vel.y, minImage_vel.y, tol); + MY_CHECK_CLOSE(test_vel.z, minImage_vel.z, tol); + } + +// begin tests +UP_TEST(BoxDim_deform_test) + { + // first test functionality of box deformation methods + BoxDim box(10); + + Scalar tol = Scalar(1e-4); + Scalar xy_rate = .1; + Scalar xz_rate = .3; + Scalar yz_rate = .5; + Scalar3 L_rate = make_scalar3(.2, .4, .6); + + box.setTiltDeformationRates(xy_rate, xz_rate, yz_rate); + box.setLDeformationRate(L_rate); + + MY_CHECK_CLOSE(box.getTiltDeformationRateXY(), xy_rate, tol); + MY_CHECK_CLOSE(box.getTiltDeformationRateXZ(), xz_rate, tol); + MY_CHECK_CLOSE(box.getTiltDeformationRateYZ(), yz_rate, tol); + MY_CHECK_CLOSE(box.getLDeformationRate().x, L_rate.x, tol); + MY_CHECK_CLOSE(box.getLDeformationRate().y, L_rate.y, tol); + MY_CHECK_CLOSE(box.getLDeformationRate().z, L_rate.z, tol); + + // starting deformation test with a cubic box + // non-deforming + deforming_box_test({5, 5, 5}, + {0, 0, 0}, + {0, 0, 0}, + {0, 0, 0}, + {1, 3, -3}, + {3, -2, 4}, + {1, -2, 2}, + {3, -2, 4}, + {1, -2, 2}, + {3, -2, 4}); + // tilted in xy; non-deforming + deforming_box_test({5, 5, 5}, + {0.1, 0, 0}, + {0, 0, 0}, + {0, 0, 0}, + {1, 3, -3}, + {3, -2, 4}, + {0.5, -2, 2}, + {3, -2, 4}, + {0.5, -2, 2}, + {3, -2, 4}); + // tilted in xy, xz; non-deforming + deforming_box_test({5, 5, 5}, + {0.1, 0.2, 0}, + {0, 0, 0}, + {0, 0, 0}, + {1, 3, -3}, + {3, -2, 4}, + {1.5, -2, 2}, + {3, -2, 4}, + {1.5, -2, 2}, + {3, -2, 4}); + // tilted in xy, xz, yz; non-deforming + deforming_box_test({5, 5, 5}, + {0.1, 0.2, 0.3}, + {0, 0, 0}, + {0, 0, 0}, + {1, 3, -3}, + {3, -2, 4}, + {1.5, -0.5, 2}, + {3, -2, 4}, + {1.5, -0.5, 2}, + {3, -2, 4}); + // tilted in xy, xz, yz; deforming in xy + deforming_box_test({5, 5, 5}, + {0.1, 0.2, 0.3}, + {0, 0, 0}, + {0.1, 0, 0}, + {1, 3, -3}, + {3, -2, 4}, + {1.5, -0.5, 2}, + {2.5, -2, 4}, + {1.5, -0.5, 2}, + {2.5, -2, 4}); + // tilted in xy, xz, yz; deforming in xy, xz + deforming_box_test({5, 5, 5}, + {0.1, 0.2, 0.3}, + {0, 0, 0}, + {0.1, 0.2, 0}, + {1, 3, -3}, + {3, -2, 4}, + {1.5, -0.5, 2}, + {3.5, -2, 4}, + {1.5, -0.5, 2}, + {3.5, -2, 4}); + // tilted in xy, xz, yz; deforming in xy, xz, yz + deforming_box_test({5, 5, 5}, + {0.1, 0.2, 0.3}, + {0, 0, 0}, + {0.1, 0.2, 0.3}, + {1, 3, -3}, + {3, -2, 4}, + {1.5, -0.5, 2}, + {3.5, -0.5, 4}, + {1.5, -0.5, 2}, + {3.5, -0.5, 4}); + // tilted in xy, xz, yz; deforming in all tilts and Lx + deforming_box_test({5, 5, 5}, + {0.1, 0.2, 0.3}, + {0.01, 0, 0}, + {0.1, 0.2, 0.3}, + {1, 3, -3}, + {3, -2, 4}, + {1.5, -0.5, 2}, + {3.5, -0.5, 4}, + {1.5, -0.5, 2}, + {3.5, -0.5, 4}); + // tilted in xy, xz, yz; deforming in all tilts, Lx and Ly + deforming_box_test({5, 5, 5}, + {0.1, 0.2, 0.3}, + {0.01, 0.02, 0}, + {0.1, 0.2, 0.3}, + {1, 3, -3}, + {3, -2, 4}, + {1.5, -0.5, 2}, + {3.5, -0.52, 4}, + {1.5, -0.5, 2}, + {3.5, -0.52, 4}); + // tilted in xy, xz, yz; deforming in all tilts and all L + deforming_box_test({5, 5, 5}, + {0.1, 0.2, 0.3}, + {0.01, 0.02, 0.03}, + {0.1, 0.2, 0.3}, + {1, 3, -3}, + {3, -2, 4}, + {1.5, -0.5, 2}, + {3.5, -0.52, 4.03}, + {1.5, -0.5, 2}, + {3.5, -0.52, 4.03}); + + // starting deformation with an orthogonal box + // non-deforming + deforming_box_test({5, 8, 10}, + {0, 0, 0}, + {0, 0, 0}, + {0, 0, 0}, + {3, 5, -7}, + {3, -2, 4}, + {-2, -3, 3}, + {3, -2, 4}, + {-2, -3, 3}, + {3, -2, 4}); + // tilted in xy; non-deforming + deforming_box_test({5, 8, 10}, + {0.1, 0, 0}, + {0, 0, 0}, + {0, 0, 0}, + {3.5, 5, -7}, + {3, -2, 4}, + {-2.3, -3, 3}, + {3, -2, 4}, + {-2.3, -3, 3}, + {3, -2, 4}); + // tilted in xy, xz; non-deforming + deforming_box_test({5, 8, 10}, + {0.1, -0.2, 0}, + {0, 0, 0}, + {0, 0, 0}, + {3, 5, -7}, + {3, -2, 4}, + {0.2, -3, 3}, + {3, -2, 4}, + {0.2, -3, 3}, + {3, -2, 4}); + // tilted in xy, xz, yz; non-deforming + deforming_box_test({5, 8, 10}, + {0.1, -0.2, 0.3}, + {0, 0, 0}, + {0, 0, 0}, + {3, 5, -7}, + {3, -2, 4}, + {0.2, 0, 3}, + {3, -2, 4}, + {0.2, 0, 3}, + {3, -2, 4}); + // tilted in xy, xz, yz; deforming in xy + deforming_box_test({5, 8, 10}, + {0.1, -0.2, 0.3}, + {0, 0, 0}, + {0.1, 0, 0}, + {3, 5, -7}, + {3, -2, 4}, + {0.2, 0, 3}, + {2.2, -2, 4}, + {0.2, 0, 3}, + {2.2, -2, 4}); + // tilted in xy, xz, yz; deforming in xy and xz + deforming_box_test({5, 8, 10}, + {0.1, -0.2, 0.3}, + {0, 0, 0}, + {0.1, 0.2, 0}, + {3, 5, -7}, + {3, -2, 4}, + {0.2, 0, 3}, + {4.2, -2, 4}, + {0.2, 0, 3}, + {4.2, -2, 4}); + // tilted in xy, xz, yz; deforming in xy, xz and yz + deforming_box_test({5, 8, 10}, + {0.1, -0.2, 0.3}, + {0, 0, 0}, + {0.1, 0.2, -0.3}, + {3, 5, -7}, + {3, -2, 4}, + {0.2, 0, 3}, + {4.2, -5, 4}, + {0.2, 0, 3}, + {4.2, -5, 4}); + // tilted in xy, xz, yz; deforming in all tilts and Lx + deforming_box_test({5, 8, 10}, + {0.1, -0.2, 0.3}, + {-0.01, 0, 0}, + {0.1, 0.2, -0.3}, + {3, 5, -7}, + {3, -2, 4}, + {0.2, 0, 3}, + {4.2, -5, 4}, + {0.2, 0, 3}, + {4.2, -5, 4}); + // tilted in xy, xz, yz; deforming in all tilts, Lx and Ly + deforming_box_test({5, 8, 10}, + {0.1, -0.2, 0.3}, + {-0.01, 0.02, 0}, + {0.1, 0.2, -0.3}, + {3, 5, -7}, + {3, -2, 4}, + {0.2, 0, 3}, + {4.2, -5.02, 4}, + {0.2, 0, 3}, + {4.2, -5.02, 4}); + // tilted in xy, xz, yz; deforming in all tilts and all L + deforming_box_test({5, 8, 10}, + {0.1, -0.2, 0.3}, + {-0.01, 0.02, 0.03}, + {0.1, 0.2, -0.3}, + {3, 5, -7}, + {3, -2, 4}, + {0.2, 0, 3}, + {4.2, -5.02, 4.03}, + {0.2, 0, 3}, + {4.2, -5.02, 4.03}); + } + //! Test operation of the particle data class UP_TEST(ParticleData_test) {