Remove PImpl scaffolding from CUDA version of LINCS
authorArtem Zhmurov <zhmurov@gmail.com>
Fri, 17 May 2019 14:08:24 +0000 (16:08 +0200)
committerArtem Zhmurov <zhmurov@gmail.com>
Fri, 5 Jul 2019 09:05:46 +0000 (11:05 +0200)
The CUDA implementation of LINCS was initially introduced as a
stand-alone feature. This required hiding CUDA-specific variables
and subroutines into the private implementation subclass. Since the
LINCS is not a part of Update and Constraints module, this is no
longer required and can be removed.

Refs #2816, #2888

Change-Id: I9698224d4702dfb8d99106999335c62e83a511df

13 files changed:
src/gromacs/mdlib/CMakeLists.txt
src/gromacs/mdlib/lincs_cuda.cu [moved from src/gromacs/mdlib/lincs_cuda_impl.cu with 92% similarity]
src/gromacs/mdlib/lincs_cuda.cuh [moved from src/gromacs/mdlib/lincs_cuda_impl.h with 76% similarity]
src/gromacs/mdlib/lincs_cuda.h [deleted file]
src/gromacs/mdlib/lincs_cuda_impl.cpp [deleted file]
src/gromacs/mdlib/tests/CMakeLists.txt
src/gromacs/mdlib/tests/constr.cpp
src/gromacs/mdlib/tests/constr_impl.cpp [new file with mode: 0644]
src/gromacs/mdlib/tests/constr_impl.cu [new file with mode: 0644]
src/gromacs/mdlib/tests/constr_impl.h [new file with mode: 0644]
src/gromacs/mdlib/update_constrain_cuda_impl.cu
src/gromacs/mdlib/update_constrain_cuda_impl.h
src/gromacs/pbcutil/pbc_aiuc.h

index d543f04651a0f45f05315961cfc02a3a455a2372..9b37d64041b76624811388e04c703449214f3b74 100644 (file)
@@ -41,7 +41,7 @@ endif()
 if(GMX_USE_CUDA)
     gmx_add_libgromacs_sources(
        leapfrog_cuda_impl.cu
-       lincs_cuda_impl.cu
+       lincs_cuda.cu
        settle_cuda_impl.cu
        update_constrain_cuda_impl.cu
        )
similarity index 92%
rename from src/gromacs/mdlib/lincs_cuda_impl.cu
rename to src/gromacs/mdlib/lincs_cuda.cu
index eb03a4821409e860a75a0874a536e6cbf9097661..df018e9e918b658ad823ec86d157c8b24774855b 100644 (file)
@@ -52,7 +52,7 @@
  */
 #include "gmxpre.h"
 
-#include "lincs_cuda_impl.h"
+#include "lincs_cuda.cuh"
 
 #include <assert.h>
 #include <stdio.h>
@@ -68,7 +68,6 @@
 #include "gromacs/gpu_utils/vectype_ops.cuh"
 #include "gromacs/math/vec.h"
 #include "gromacs/mdlib/constr.h"
-#include "gromacs/mdlib/lincs_cuda.h"
 #include "gromacs/pbcutil/pbc.h"
 #include "gromacs/pbcutil/pbc_aiuc_cuda.cuh"
 #include "gromacs/topology/ifunc.h"
@@ -431,13 +430,13 @@ inline auto getLincsKernelPtr(const bool  updateVelocities,
     return kernelPtr;
 }
 
-void LincsCuda::Impl::apply(const float3 *d_x,
-                            float3       *d_xp,
-                            const bool    updateVelocities,
-                            float3       *d_v,
-                            const real    invdt,
-                            const bool    computeVirial,
-                            tensor        virialScaled)
+void LincsCuda::apply(const float3 *d_x,
+                      float3       *d_xp,
+                      const bool    updateVelocities,
+                      float3       *d_v,
+                      const real    invdt,
+                      const bool    computeVirial,
+                      tensor        virialScaled)
 {
     ensureNoPendingCudaError("In CUDA version of LINCS");
 
@@ -519,45 +518,8 @@ void LincsCuda::Impl::apply(const float3 *d_x,
     return;
 }
 
-void LincsCuda::Impl::copyApplyCopy(const int   numAtoms,
-                                    const rvec *h_x,
-                                    rvec       *h_xp,
-                                    bool        updateVelocities,
-                                    rvec       *h_v,
-                                    real        invdt,
-                                    bool        computeVirial,
-                                    tensor      virialScaled)
-{
-
-    float3 *d_x, *d_xp, *d_v;
-
-    allocateDeviceBuffer(&d_x,  numAtoms, nullptr);
-    allocateDeviceBuffer(&d_xp, numAtoms, nullptr);
-    allocateDeviceBuffer(&d_v,  numAtoms, nullptr);
-
-    copyToDeviceBuffer(&d_x, (float3*)h_x, 0, numAtoms, stream_, GpuApiCallBehavior::Sync, nullptr);
-    copyToDeviceBuffer(&d_xp, (float3*)h_xp, 0, numAtoms, stream_, GpuApiCallBehavior::Sync, nullptr);
-    if (updateVelocities)
-    {
-        copyToDeviceBuffer(&d_v, (float3*)h_v, 0, numAtoms, stream_, GpuApiCallBehavior::Sync, nullptr);
-    }
-    apply(d_x, d_xp,
-          updateVelocities, d_v, invdt,
-          computeVirial, virialScaled);
-
-    copyFromDeviceBuffer((float3*)h_xp, &d_xp, 0, numAtoms, stream_, GpuApiCallBehavior::Sync, nullptr);
-    if (updateVelocities)
-    {
-        copyFromDeviceBuffer((float3*)h_v, &d_v, 0, numAtoms, stream_, GpuApiCallBehavior::Sync, nullptr);
-    }
-
-    freeDeviceBuffer(&d_x);
-    freeDeviceBuffer(&d_xp);
-    freeDeviceBuffer(&d_v);
-}
-
-LincsCuda::Impl::Impl(int numIterations,
-                      int expansionOrder)
+LincsCuda::LincsCuda(int numIterations,
+                     int expansionOrder)
 {
     kernelParams_.numIterations              = numIterations;
     kernelParams_.expansionOrder             = expansionOrder;
@@ -579,7 +541,7 @@ LincsCuda::Impl::Impl(int numIterations,
 
 }
 
-LincsCuda::Impl::~Impl()
+LincsCuda::~LincsCuda()
 {
     freeDeviceBuffer(&kernelParams_.d_virialScaled);
 
@@ -635,8 +597,8 @@ inline int countCoupled(int a, std::vector<int> *spaceNeeded,
     return counted;
 }
 
-void LincsCuda::Impl::set(const t_idef    &idef,
-                          const t_mdatoms &md)
+void LincsCuda::set(const t_idef    &idef,
+                    const t_mdatoms &md)
 {
     int                 numAtoms = md.nr;
     // List of constrained atoms (CPU memory)
@@ -891,42 +853,9 @@ void LincsCuda::Impl::set(const t_idef    &idef,
 
 }
 
-void LincsCuda::Impl::setPbc(const t_pbc *pbc)
-{
-    setPbcAiuc(pbc->ndim_ePBC, pbc->box, &kernelParams_.pbcAiuc);
-}
-
-LincsCuda::LincsCuda(const int numIterations,
-                     const int expansionOrder)
-    : impl_(new Impl(numIterations, expansionOrder))
-{
-}
-
-LincsCuda::~LincsCuda() = default;
-
-void LincsCuda::copyApplyCopy(const int   numAtoms,
-                              const rvec *h_x,
-                              rvec       *h_xp,
-                              const bool  updateVelocities,
-                              rvec       *h_v,
-                              const real  invdt,
-                              const bool  computeVirial,
-                              tensor      virialScaled)
-{
-    impl_->copyApplyCopy(numAtoms, h_x, h_xp,
-                         updateVelocities, h_v, invdt,
-                         computeVirial, virialScaled);
-}
-
 void LincsCuda::setPbc(const t_pbc *pbc)
 {
-    impl_->setPbc(pbc);
-}
-
-void LincsCuda::set(const t_idef    &idef,
-                    const t_mdatoms &md)
-{
-    impl_->set(idef, md);
+    setPbcAiuc(pbc->ndim_ePBC, pbc->box, &kernelParams_.pbcAiuc);
 }
 
 } // namespace gmx
similarity index 76%
rename from src/gromacs/mdlib/lincs_cuda_impl.h
rename to src/gromacs/mdlib/lincs_cuda.cuh
index ff21ae88aeb53cab25884660b4b921f675a234db..1fb9eada234f82cc900f843ca9f7953cb9e52b1d 100644 (file)
  * To help us fund GROMACS development, we humbly ask that you cite
  * the research papers on the package. Check out http://www.gromacs.org.
  */
-/*! \internal \file
+/*! \libinternal \file
  *
- * \brief Declares CUDA implementation class for LINCS
- *
- * This header file is needed to include from both the device-side
- * kernels file, and the host-side management code.
+ * \brief Declaration of high-level functions of CUDA implementation of LINCS.
  *
  * \author Artem Zhmurov <zhmurov@gmail.com>
  *
  * \ingroup module_mdlib
+ * \inlibraryapi
  */
-#ifndef GMX_MDLIB_LINCS_CUDA_IMPL_H
-#define GMX_MDLIB_LINCS_CUDA_IMPL_H
-
+#ifndef GMX_MDLIB_LINCS_CUDA_CUH
+#define GMX_MDLIB_LINCS_CUDA_CUH
 
 #include "gromacs/mdlib/constr.h"
-#include "gromacs/mdlib/lincs_cuda.h"
 #include "gromacs/mdtypes/mdatom.h"
-#include "gromacs/pbcutil/pbc_aiuc_cuda.cuh"
+#include "gromacs/pbcutil/pbc_aiuc.h"
 #include "gromacs/topology/idef.h"
+#include "gromacs/utility/classhelpers.h"
 
 namespace gmx
 {
@@ -96,7 +93,7 @@ struct LincsCudaKernelParameters
 };
 
 /*! \internal \brief Class with interfaces and data for CUDA version of LINCS. */
-class LincsCuda::Impl
+class LincsCuda
 {
 
     public:
@@ -106,10 +103,10 @@ class LincsCuda::Impl
          * \param[in] numIterations    Number of iteration for the correction of the projection.
          * \param[in] expansionOrder   Order of the matrix inversion algorithm.
          */
-        Impl(int numIterations,
-             int expansionOrder);
+        LincsCuda(int numIterations,
+                  int expansionOrder);
         /*! \brief Destructor.*/
-        ~Impl();
+        ~LincsCuda();
 
         /*! \brief Apply LINCS.
          *
@@ -137,38 +134,6 @@ class LincsCuda::Impl
                    const bool    computeVirial,
                    tensor        virialScaled);
 
-        /*! \brief Apply LINCS to the coordinates/velocities stored in CPU memory.
-         *
-         * This method should not be used in any code-path, where performance is of any value.
-         * Only suitable for test and will be removed in future patch sets.
-         * Allocates GPU memory, copies data from CPU, applies LINCS to coordinates and,
-         * if requested, to velocities, copies the results back, frees GPU memory.
-         * Method uses this class data structures which should be filled with set() and setPbc()
-         * methods.
-         *
-         * \todo Remove this method
-         *
-         * \param[in]     numAtoms          Number of atoms
-         * \param[in]     h_x               Coordinates before timestep (in CPU memory)
-         * \param[in,out] h_xp              Coordinates after timestep (in CPU memory). The
-         *                                  resulting constrained coordinates will be saved here.
-         * \param[in]     updateVelocities  If the velocities should be updated.
-         * \param[in,out] h_v               Velocities to update (in CPU memory, can be nullptr
-         *                                  if not updated)
-         * \param[in]     invdt             Reciprocal timestep (to scale Lagrange
-         *                                  multipliers when velocities are updated)
-         * \param[in]     computeVirial     If virial should be updated.
-         * \param[in,out] virialScaled      Scaled virial tensor to be updated.
-         */
-        void copyApplyCopy(const int   numAtoms,
-                           const rvec *h_x,
-                           rvec       *h_xp,
-                           const bool  updateVelocities,
-                           rvec       *h_v,
-                           const real  invdt,
-                           const bool  computeVirial,
-                           tensor      virialScaled);
-
         /*! \brief
          * Update data-structures (e.g. after NB search step).
          *
@@ -234,6 +199,6 @@ class LincsCuda::Impl
 
 };
 
-} // namespace gmx
+}      // namespace gmx
 
-#endif
+#endif // GMX_MDLIB_LINCS_CUDA_CUH
diff --git a/src/gromacs/mdlib/lincs_cuda.h b/src/gromacs/mdlib/lincs_cuda.h
deleted file mode 100644 (file)
index 7601493..0000000
+++ /dev/null
@@ -1,138 +0,0 @@
-/*
- * This file is part of the GROMACS molecular simulation package.
- *
- * Copyright (c) 2019, by the GROMACS development team, led by
- * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
- * and including many others, as listed in the AUTHORS file in the
- * top-level source directory and at http://www.gromacs.org.
- *
- * GROMACS is free software; you can redistribute it and/or
- * modify it under the terms of the GNU Lesser General Public License
- * as published by the Free Software Foundation; either version 2.1
- * of the License, or (at your option) any later version.
- *
- * GROMACS is distributed in the hope that it will be useful,
- * but WITHOUT ANY WARRANTY; without even the implied warranty of
- * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
- * Lesser General Public License for more details.
- *
- * You should have received a copy of the GNU Lesser General Public
- * License along with GROMACS; if not, see
- * http://www.gnu.org/licenses, or write to the Free Software Foundation,
- * Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA.
- *
- * If you want to redistribute modifications to GROMACS, please
- * consider that scientific software is very special. Version
- * control is crucial - bugs must be traceable. We will be happy to
- * consider code for inclusion in the official distribution, but
- * derived work must not be called official GROMACS. Details are found
- * in the README & COPYING files - if they are missing, get the
- * official version at http://www.gromacs.org.
- *
- * To help us fund GROMACS development, we humbly ask that you cite
- * the research papers on the package. Check out http://www.gromacs.org.
- */
-/*! \libinternal \file
- *
- * \brief Declaration of high-level functions of CUDA implementation of LINCS.
- *
- * \todo This should only list interfaces needed for libgromacs clients.
- *
- * \author Artem Zhmurov <zhmurov@gmail.com>
- *
- * \ingroup module_mdlib
- * \inlibraryapi
- */
-#ifndef GMX_MDLIB_LINCS_CUDA_H
-#define GMX_MDLIB_LINCS_CUDA_H
-
-
-#include "gromacs/mdlib/constr.h"
-#include "gromacs/mdtypes/mdatom.h"
-#include "gromacs/topology/idef.h"
-#include "gromacs/utility/classhelpers.h"
-
-namespace gmx
-{
-
-class LincsCuda
-{
-
-    public:
-        /*! \brief Constructor.
-         *
-         * \param[in] numIterations    Number of iteration for the correction of the projection.
-         * \param[in] expansionOrder   Order of the matrix inversion algorithm.
-         */
-        LincsCuda(int numIterations,
-                  int expansionOrder);
-        ~LincsCuda();
-
-        /*! \brief Apply LINCS to the coordinates/velocities stored in CPU memory.
-         *
-         * This method should not be used in any code-path, where performance is of any value.
-         * Only suitable for test and will be removed in future patch sets.
-         * Allocates GPU memory, copies data from CPU, applies LINCS to coordinates and,
-         * if requested, to velocities, copies the results back, frees GPU memory.
-         * Method uses this class data structures which should be filled with set() and setPbc()
-         * methods.
-         *
-         * \todo Remove this method
-         *
-         * \param[in]     numAtoms          Number of atoms
-         * \param[in]     h_x               Coordinates before timestep (in CPU memory)
-         * \param[in,out] h_xp              Coordinates after timestep (in CPU memory). The
-         *                                  resulting constrained coordinates will be saved here.
-         * \param[in]     updateVelocities  If the velocities should be updated.
-         * \param[in,out] h_v               Velocities to update (in CPU memory, can be nullptr
-         *                                  if not updated)
-         * \param[in]     invdt             Reciprocal timestep (to scale Lagrange
-         *                                  multipliers when velocities are updated)
-         * \param[in]     computeVirial     If virial should be updated.
-         * \param[in,out] virialScaled      Scaled virial tensor to be updated.
-         */
-        void copyApplyCopy(int         numAtoms,
-                           const rvec *h_x,
-                           rvec       *h_xp,
-                           bool        updateVelocities,
-                           rvec       *h_v,
-                           real        invdt,
-                           bool        computeVirial,
-                           tensor      virialScaled);
-
-        /*! \brief
-         * Update data-structures (e.g. after NB search step).
-         *
-         * Updates the constraints data. Should be called if the particles were sorted,
-         * redistributed between domains, etc.
-         *
-         * Information about constraints should be taken from:
-         *     idef.il[F_CONSTR].iatoms  --- type (T) of constraint and two atom indexes (i1, i2)
-         *     idef.iparams[T].constr.dA --- target length for constraint of type T
-         * From t_mdatom, the code should take:
-         *     md.invmass  --- array of inverse square root of masses for each atom in the system.
-         *
-         * \param[in] idef  Local topology data to get information on constraints from.
-         * \param[in] md    Atoms data to get atom masses from.
-         */
-        void set(const t_idef         &idef,
-                 const t_mdatoms      &md);
-
-        /*! \brief
-         * Update PBC data.
-         *
-         * \param[in] pbc  The PBC data in t_pbc format.
-         */
-        void setPbc(const t_pbc *pbc);
-
-        /*! \brief Class with hardware-specific interfaces and implementations.*/
-        class Impl;
-
-    private:
-        gmx::PrivateImplPointer<Impl> impl_;
-
-};
-
-} // namespace gmx
-
-#endif
diff --git a/src/gromacs/mdlib/lincs_cuda_impl.cpp b/src/gromacs/mdlib/lincs_cuda_impl.cpp
deleted file mode 100644 (file)
index 39ab17d..0000000
+++ /dev/null
@@ -1,96 +0,0 @@
-/*
- * This file is part of the GROMACS molecular simulation package.
- *
- * Copyright (c) 2019, by the GROMACS development team, led by
- * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
- * and including many others, as listed in the AUTHORS file in the
- * top-level source directory and at http://www.gromacs.org.
- *
- * GROMACS is free software; you can redistribute it and/or
- * modify it under the terms of the GNU Lesser General Public License
- * as published by the Free Software Foundation; either version 2.1
- * of the License, or (at your option) any later version.
- *
- * GROMACS is distributed in the hope that it will be useful,
- * but WITHOUT ANY WARRANTY; without even the implied warranty of
- * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
- * Lesser General Public License for more details.
- *
- * You should have received a copy of the GNU Lesser General Public
- * License along with GROMACS; if not, see
- * http://www.gnu.org/licenses, or write to the Free Software Foundation,
- * Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA.
- *
- * If you want to redistribute modifications to GROMACS, please
- * consider that scientific software is very special. Version
- * control is crucial - bugs must be traceable. We will be happy to
- * consider code for inclusion in the official distribution, but
- * derived work must not be called official GROMACS. Details are found
- * in the README & COPYING files - if they are missing, get the
- * official version at http://www.gromacs.org.
- *
- * To help us fund GROMACS development, we humbly ask that you cite
- * the research papers on the package. Check out http://www.gromacs.org.
- */
-/*! \internal \file
- *
- * \brief May be used to implement LINCS CUDA interfaces for non-GPU builds.
- *
- * Currently, reports and exits if any of the interfaces are called.
- * May replace CPU LINCS module in future. Currently used to satisfy
- * compiler on systems, where CUDA is not available.
- *
- * \author Artem Zhmurov <zhmurov@gmail.com>
- *
- * \ingroup module_mdlib
- */
-#include "gmxpre.h"
-
-#include "config.h"
-
-#include "gromacs/mdlib/lincs_cuda.h"
-
-#if GMX_GPU != GMX_GPU_CUDA
-
-namespace gmx
-{
-
-class LincsCuda::Impl
-{
-};
-
-LincsCuda::LincsCuda(gmx_unused int numIterations,
-                     gmx_unused int expansionOrder)
-    : impl_(nullptr)
-{
-    GMX_ASSERT(false, "A CPU stub for LINCS was called insted of the correct implementation.");
-}
-
-LincsCuda::~LincsCuda() = default;
-
-void LincsCuda::copyApplyCopy(gmx_unused int         numAtoms,
-                              gmx_unused const rvec *h_x,
-                              gmx_unused rvec       *h_xp,
-                              gmx_unused bool        updateVelocities,
-                              gmx_unused rvec       *h_v,
-                              gmx_unused real        invdt,
-                              gmx_unused bool        computeVirial,
-                              gmx_unused tensor      virialScaled)
-{
-    GMX_ASSERT(false, "A CPU stub for LINCS was called insted of the correct implementation.");
-}
-
-void LincsCuda::set(gmx_unused const t_idef    &idef,
-                    gmx_unused const t_mdatoms &md)
-{
-    GMX_ASSERT(false, "A CPU stub for LINCS was called insted of the correct implementation.");
-}
-
-void LincsCuda::setPbc(gmx_unused const t_pbc *pbc)
-{
-    GMX_ASSERT(false, "A CPU stub for LINCS was called insted of the correct implementation.");
-}
-
-}      // namespace gmx
-
-#endif /* GMX_GPU != GMX_GPU_CUDA */
index fc0c714c1af67a3656838fc5c01f0598a87260db..e3ac4b1cfe281b44470c88ecd334d35bab44374d 100644 (file)
 # To help us fund GROMACS development, we humbly ask that you cite
 # the research papers on the package. Check out http://www.gromacs.org.
 
+file(GLOB MDLIB_TEST_SOURCES
+     constr_impl.cpp
+     )
+
+if (GMX_USE_CUDA)
+    file(GLOB MDLIB_TEST_CUDA_SOURCES
+         constr_impl.cu
+         )
+endif()
+
 gmx_add_unit_test(MdlibUnitTest mdlib-test
                   calc_verletbuf.cpp
                   constr.cpp
@@ -42,4 +52,8 @@ gmx_add_unit_test(MdlibUnitTest mdlib-test
                   shake.cpp
                   simulationsignal.cpp
                   updategroups.cpp
-                  updategroupscog.cpp)
+                  updategroupscog.cpp
+                  ${MDLIB_TEST_SOURCES})
+
+# TODO: Make CUDA source to compile inside the testing framework
+gmx_add_libgromacs_sources(${MDLIB_TEST_CUDA_SOURCES})
\ No newline at end of file
index 69a0260ed87c62991fdfdfd6c7c9a9bbc4edfa4f..2f02a2ca4c9cf3885d81d3fe5c368e45fffde311 100644 (file)
 
 #include "gromacs/fileio/gmxfio.h"
 #include "gromacs/gmxlib/nrnb.h"
-#include "gromacs/gmxlib/nonbonded/nonbonded.h"
-#include "gromacs/gpu_utils/gpu_utils.h"
 #include "gromacs/math/paddedvector.h"
 #include "gromacs/math/vec.h"
 #include "gromacs/math/vectypes.h"
-#include "gromacs/mdlib/gmx_omp_nthreads.h"
 #include "gromacs/mdlib/lincs.h"
-#include "gromacs/mdlib/lincs_cuda.h"
 #include "gromacs/mdlib/shake.h"
 #include "gromacs/mdrunutility/multisim.h"
 #include "gromacs/mdtypes/commrec.h"
 #include "testutils/refdata.h"
 #include "testutils/testasserts.h"
 
+#include "constr_impl.h"
+
 namespace gmx
 {
 namespace test
 {
 
-/*! \brief
- * Constraints test data structure.
- *
- * Structure to collect all the necessary data, including system coordinates and topology,
- * constraints information, etc. The structure can be reset and reused.
- */
-struct ConstraintsTestData
+ConstraintsTestData::ConstraintsTestData(const std::string &title,
+                                         int numAtoms, std::vector<real> masses,
+                                         std::vector<int> constraints, std::vector<real> constraintsR0,
+                                         bool computeVirial, tensor virialScaledRef,
+                                         bool compute_dHdLambda, float dHdLambdaRef,
+                                         real initialTime, real timestep,
+                                         const std::vector<RVec> &x, const std::vector<RVec> &xPrime, const std::vector<RVec> &v,
+                                         real shakeTolerance, gmx_bool shakeUseSOR,
+                                         int lincsNumIterations, int lincsExpansionOrder, real lincsWarnAngle)
 {
-    public:
-        //! Human-friendly name for a system
-        std::string           title_;
-        //! Number of atoms
-        int                   numAtoms_;
-        //! Topology
-        gmx_mtop_t            mtop_;
-        //! Masses
-        std::vector<real>     masses_;
-        //! Inverse masses
-        std::vector<real>     invmass_;
-        //! Communication record
-        t_commrec             cr_;
-        //! Input record (info that usually in .mdp file)
-        t_inputrec            ir_;
-        //! Local topology
-        t_idef                idef_;
-        //! MD atoms
-        t_mdatoms             md_;
-        //! Multisim data
-        gmx_multisim_t        ms_;
-        //! Computational time array (normally used to benchmark performance)
-        t_nrnb                nrnb_;
-
-        //! Inverse timestep
-        real                  invdt_;
-        //! Number of flexible constraints
-        int                   nflexcon_   = 0;
-        //! Whether the virial should be computed
-        bool                  computeVirial_;
-        //! Scaled virial
-        tensor                virialScaled_;
-        //! Scaled virial (reference values)
-        tensor                virialScaledRef_;
-        //! If the free energy is computed
-        bool                  compute_dHdLambda_;
-        //! For free energy computation
-        real                  dHdLambda_;
-        //! For free energy computation (reference value)
-        real                  dHdLambdaRef_;
-
-        //! Coordinates before the timestep
-        PaddedVector<RVec>    x_;
-        //! Coordinates after timestep, output for the constraints
-        PaddedVector<RVec>    xPrime_;
-        //! Backup for coordinates (for reset)
-        PaddedVector<RVec>    xPrime0_;
-        //! Intermediate set of coordinates (normally used for projection correction)
-        PaddedVector<RVec>    xPrime2_;
-        //! Velocities
-        PaddedVector<RVec>    v_;
-        //! Backup for velocities (for reset)
-        PaddedVector<RVec>    v0_;
-
-        //! Constraints data (type1-i1-j1-type2-i2-j2-...)
-        std::vector<int>      constraints_;
-        //! Target lengths for all constraint types
-        std::vector<real>     constraintsR0_;
-
-        /*! \brief
-         * Constructor for the object with all parameters and variables needed by constraints algorithms.
-         *
-         * This constructor assembles stubs for all the data structures, required to initialize
-         * and apply LINCS and SHAKE constraints. The coordinates and velocities before constraining
-         * are saved to allow for reset. The constraints data are stored for testing after constraints
-         * were applied.
-         *
-         * \param[in]  title                Human-friendly name of the system.
-         * \param[in]  numAtoms             Number of atoms in the system.
-         * \param[in]  masses               Atom masses. Size of this vector should be equal to numAtoms.
-         * \param[in]  constraints          List of constraints, organized in triples of integers.
-         *                                  First integer is the index of type for a constraint, second
-         *                                  and third are the indices of constrained atoms. The types
-         *                                  of constraints should be sequential but not necessarily
-         *                                  start from zero (which is the way they normally are in
-         *                                  GROMACS).
-         * \param[in]  constraintsR0        Target values for bond lengths for bonds of each type. The
-         *                                  size of this vector should be equal to the total number of
-         *                                  unique types in constraints vector.
-         * \param[in]  computeVirial        Whether the virial should be computed.
-         * \param[in]  virialScaledRef      Reference values for scaled virial tensor.
-         * \param[in]  compute_dHdLambda    Whether free energy should be computed.
-         * \param[in]  dHdLambdaRef         Reference value for dHdLambda.
-         * \param[in]  initialTime          Initial time.
-         * \param[in]  timestep             Timestep.
-         * \param[in]  x                    Coordinates before integration step.
-         * \param[in]  xPrime               Coordinates after integration step, but before constraining.
-         * \param[in]  v                    Velocities before constraining.
-         * \param[in]  shakeTolerance       Target tolerance for SHAKE.
-         * \param[in]  shakeUseSOR          Use successive over-relaxation method for SHAKE iterations.
-         *                                  The general formula is:
-         *                                     x_n+1 = (1-omega)*x_n + omega*f(x_n),
-         *                                  where omega = 1 if SOR is off and may be < 1 if SOR is on.
-         * \param[in]  lincsNumIterations   Number of iterations used to compute the inverse matrix.
-         * \param[in]  lincsExpansionOrder  The order for algorithm that adjusts the direction of the
-         *                                  bond after constraints are applied.
-         * \param[in]  lincsWarnAngle       The threshold value for the change in bond angle. When
-         *                                  exceeded the program will issue a warning.
-         *
-         */
-        ConstraintsTestData(const std::string &title,
-                            int numAtoms, std::vector<real> masses,
-                            std::vector<int> constraints, std::vector<real> constraintsR0,
-                            bool computeVirial, tensor virialScaledRef,
-                            bool compute_dHdLambda, float dHdLambdaRef,
-                            real initialTime, real timestep,
-                            const std::vector<RVec> &x, const std::vector<RVec> &xPrime, const std::vector<RVec> &v,
-                            real shakeTolerance, gmx_bool shakeUseSOR,
-                            int lincsNumIterations, int lincsExpansionOrder, real lincsWarnAngle)
+    // This is to trick Gerrit
+    {
         {
             title_    = title;    // Human-friendly name of the system
             numAtoms_ = numAtoms; // Number of atoms
@@ -361,43 +254,46 @@ struct ConstraintsTestData
             ir_.nProjOrder     = lincsExpansionOrder;
             ir_.LincsWarnAngle = lincsWarnAngle;
         }
+    }
+}
 
-        /*! \brief
        * Reset the data structure so it can be reused.
        *
        * Set the coordinates and velocities back to their values before
        * constraining. The scaled virial tensor and dHdLambda are zeroed.
        *
        */
-        void reset()
-        {
-            xPrime_  = xPrime0_;
-            xPrime2_ = xPrime0_;
-            v_       = v0_;
+/*! \brief
+ * Reset the data structure so it can be reused.
+ *
+ * Set the coordinates and velocities back to their values before
+ * constraining. The scaled virial tensor and dHdLambda are zeroed.
+ *
+ */
+void ConstraintsTestData::reset()
+{
+    xPrime_  = xPrime0_;
+    xPrime2_ = xPrime0_;
+    v_       = v0_;
 
-            if (computeVirial_)
+    if (computeVirial_)
+    {
+        for (int i = 0; i < DIM; i++)
+        {
+            for (int j = 0; j < DIM; j++)
             {
-                for (int i = 0; i < DIM; i++)
-                {
-                    for (int j = 0; j < DIM; j++)
-                    {
-                        virialScaled_[i][j] = 0;
-                    }
-                }
+                virialScaled_[i][j] = 0;
             }
-            dHdLambda_         = 0;
         }
+    }
+    dHdLambda_         = 0;
+}
 
-        /*! \brief
        * Cleaning up the memory.
        */
-        ~ConstraintsTestData()
-        {
-            sfree(idef_.il[F_CONSTR].iatoms);
-            sfree(idef_.iparams);
-        }
+/*! \brief
+ * Cleaning up the memory.
+ */
+ConstraintsTestData::~ConstraintsTestData()
+{
+    sfree(idef_.il[F_CONSTR].iatoms);
+    sfree(idef_.iparams);
+}
 
-};
+namespace
+{
 
 /*! \brief The two-dimensional parameter space for test.
  *
@@ -408,14 +304,11 @@ struct ConstraintsTestData
  */
 typedef std::tuple<std::string, std::string> ConstraintsTestParameters;
 
-/*! \brief Names of all availible algorithms
- */
-static std::vector<std::string> algorithmsNames;
+//! Names of all availible algorithms
+std::vector<std::string> algorithmsNames;
 
-/*! \brief Method that fills and returns algorithmNames to the test macros.
- *
- */
-static std::vector<std::string> getAlgorithmsNames()
+//! Method that fills and returns algorithmNames to the test macros.
+std::vector<std::string> getAlgorithmsNames()
 {
     algorithmsNames.emplace_back("SHAKE");
     algorithmsNames.emplace_back("LINCS");
@@ -484,126 +377,10 @@ class ConstraintsTest : public ::testing::TestWithParam<ConstraintsTestParameter
             algorithms_["SHAKE"] = applyShake;
             // LINCS
             algorithms_["LINCS"] = applyLincs;
-            // LINCS using CUDA (will be called only if CUDA is available)
+            // LINCS using CUDA (will only be called if CUDA is available)
             algorithms_["LINCS_CUDA"] = applyLincsCuda;
         }
 
-        /*! \brief
-         * Initialize and apply SHAKE constraints.
-         *
-         * \param[in] testData        Test data structure.
-         * \param[in] pbc             Periodic boundary data (not used in SHAKE).
-         */
-        static void applyShake(ConstraintsTestData *testData, t_pbc gmx_unused pbc)
-        {
-            shakedata* shaked = shake_init();
-            make_shake_sblock_serial(shaked, &testData->idef_, testData->md_);
-            bool       success = constrain_shake(
-                        nullptr,
-                        shaked,
-                        testData->invmass_.data(),
-                        testData->idef_,
-                        testData->ir_,
-                        as_rvec_array(testData->x_.data()),
-                        as_rvec_array(testData->xPrime_.data()),
-                        as_rvec_array(testData->xPrime2_.data()),
-                        &testData->nrnb_,
-                        testData->md_.lambda,
-                        &testData->dHdLambda_,
-                        testData->invdt_,
-                        as_rvec_array(testData->v_.data()),
-                        testData->computeVirial_,
-                        testData->virialScaled_,
-                        false,
-                        gmx::ConstraintVariable::Positions);
-            EXPECT_TRUE(success) << "Test failed with a false return value in SHAKE.";
-            done_shake(shaked);
-        }
-
-        /*! \brief
-         * Initialize and apply LINCS constraints.
-         *
-         * \param[in] testData        Test data structure.
-         * \param[in] pbc             Periodic boundary data.
-         */
-        static void applyLincs(ConstraintsTestData *testData, t_pbc pbc)
-        {
-
-            Lincs                *lincsd;
-            int                   maxwarn         = 100;
-            int                   warncount_lincs = 0;
-            gmx_omp_nthreads_set(emntLINCS, 1);
-
-            // Make blocka structure for faster LINCS setup
-            std::vector<t_blocka> at2con_mt;
-            at2con_mt.reserve(testData->mtop_.moltype.size());
-            for (const gmx_moltype_t &moltype : testData->mtop_.moltype)
-            {
-                // This function is in constr.cpp
-                at2con_mt.push_back(make_at2con(moltype,
-                                                testData->mtop_.ffparams.iparams,
-                                                flexibleConstraintTreatment(EI_DYNAMICS(testData->ir_.eI))));
-            }
-            // Initialize LINCS
-            lincsd = init_lincs(nullptr, testData->mtop_,
-                                testData->nflexcon_, at2con_mt,
-                                false,
-                                testData->ir_.nLincsIter, testData->ir_.nProjOrder);
-            set_lincs(testData->idef_, testData->md_, EI_DYNAMICS(testData->ir_.eI), &testData->cr_, lincsd);
-
-            // Evaluate constraints
-            bool sucess = constrain_lincs(false, testData->ir_, 0, lincsd, testData->md_,
-                                          &testData->cr_,
-                                          &testData->ms_,
-                                          as_rvec_array(testData->x_.data()),
-                                          as_rvec_array(testData->xPrime_.data()),
-                                          as_rvec_array(testData->xPrime2_.data()),
-                                          pbc.box, &pbc,
-                                          testData->md_.lambda, &testData->dHdLambda_,
-                                          testData->invdt_,
-                                          as_rvec_array(testData->v_.data()),
-                                          testData->computeVirial_, testData->virialScaled_,
-                                          gmx::ConstraintVariable::Positions,
-                                          &testData->nrnb_,
-                                          maxwarn, &warncount_lincs);
-            EXPECT_TRUE(sucess) << "Test failed with a false return value in LINCS.";
-            EXPECT_EQ(warncount_lincs, 0) << "There were warnings in LINCS.";
-            for (unsigned int i = 0; i < testData->mtop_.moltype.size(); i++)
-            {
-                sfree(at2con_mt.at(i).index);
-                sfree(at2con_mt.at(i).a);
-            }
-            done_lincs(lincsd);
-        }
-
-        /*! \brief
-         * Initialize and apply LINCS constraints on CUDA-enabled GPU.
-         *
-         * \param[in] testData        Test data structure.
-         * \param[in] pbc             Periodic boundary data.
-         */
-        static void applyLincsCuda(ConstraintsTestData *testData, t_pbc pbc)
-        {
-            std::string errorMessage;
-            GMX_RELEASE_ASSERT(GMX_GPU == GMX_GPU_CUDA,
-                               "The test for CUDA version of LINCS was called non-CUDA build.");
-            GMX_RELEASE_ASSERT(canDetectGpus(&errorMessage),
-                               "The test for CUDA version of LINCS was called on host that is not CUDA capable.");
-
-            auto lincsCuda = std::make_unique<LincsCuda>(testData->ir_.nLincsIter,
-                                                         testData->ir_.nProjOrder);
-            lincsCuda->set(testData->idef_, testData->md_);
-            lincsCuda->setPbc(&pbc);
-            lincsCuda->copyApplyCopy(testData->numAtoms_,
-                                     as_rvec_array(testData->x_.data()),
-                                     as_rvec_array(testData->xPrime_.data()),
-                                     true,
-                                     as_rvec_array(testData->v_.data()),
-                                     testData->invdt_,
-                                     testData->computeVirial_,
-                                     testData->virialScaled_);
-        }
-
         /*! \brief
          * The test on the final length of constrained bonds.
          *
@@ -1172,5 +949,6 @@ INSTANTIATE_TEST_CASE_P(WithParameters, ConstraintsTest,
                             ::testing::Combine(::testing::Values("PBCNone", "PBCXYZ"),
                                                    ::testing::ValuesIn(getAlgorithmsNames())));
 
+} // namespace
 } // namespace test
 } // namespace gmx
diff --git a/src/gromacs/mdlib/tests/constr_impl.cpp b/src/gromacs/mdlib/tests/constr_impl.cpp
new file mode 100644 (file)
index 0000000..452fbd7
--- /dev/null
@@ -0,0 +1,186 @@
+/*
+ * This file is part of the GROMACS molecular simulation package.
+ *
+ * Copyright (c) 2018,2019, by the GROMACS development team, led by
+ * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
+ * and including many others, as listed in the AUTHORS file in the
+ * top-level source directory and at http://www.gromacs.org.
+ *
+ * GROMACS is free software; you can redistribute it and/or
+ * modify it under the terms of the GNU Lesser General Public License
+ * as published by the Free Software Foundation; either version 2.1
+ * of the License, or (at your option) any later version.
+ *
+ * GROMACS is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
+ * Lesser General Public License for more details.
+ *
+ * You should have received a copy of the GNU Lesser General Public
+ * License along with GROMACS; if not, see
+ * http://www.gnu.org/licenses, or write to the Free Software Foundation,
+ * Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA.
+ *
+ * If you want to redistribute modifications to GROMACS, please
+ * consider that scientific software is very special. Version
+ * control is crucial - bugs must be traceable. We will be happy to
+ * consider code for inclusion in the official distribution, but
+ * derived work must not be called official GROMACS. Details are found
+ * in the README & COPYING files - if they are missing, get the
+ * official version at http://www.gromacs.org.
+ *
+ * To help us fund GROMACS development, we humbly ask that you cite
+ * the research papers on the package. Check out http://www.gromacs.org.
+ */
+/*! \internal \file
+ * \brief Function to run SHAKE and LINCS on CPU.
+ *
+ * Functions used in the test to apply constraints on the test data:
+ * CPU-based implementation and a stub for GPU-based implementation.
+ *
+ * \author Artem Zhmurov <zhmurov@gmail.com>
+ * \ingroup module_mdlib
+ */
+
+#include "gmxpre.h"
+
+#include "constr_impl.h"
+
+#include "config.h"
+
+#include <assert.h>
+
+#include <cmath>
+
+#include <algorithm>
+#include <unordered_map>
+#include <vector>
+
+#include <gtest/gtest.h>
+
+#include "gromacs/math/paddedvector.h"
+#include "gromacs/math/vec.h"
+#include "gromacs/math/vectypes.h"
+#include "gromacs/mdlib/constr.h"
+#include "gromacs/mdlib/lincs.h"
+#include "gromacs/mdlib/shake.h"
+#include "gromacs/mdtypes/commrec.h"
+#include "gromacs/mdtypes/inputrec.h"
+#include "gromacs/mdtypes/mdatom.h"
+#include "gromacs/pbcutil/pbc.h"
+#include "gromacs/topology/block.h"
+#include "gromacs/topology/idef.h"
+#include "gromacs/topology/ifunc.h"
+#include "gromacs/topology/topology.h"
+#include "gromacs/utility/unique_cptr.h"
+
+#include "testutils/testasserts.h"
+
+namespace gmx
+{
+namespace test
+{
+
+/*! \brief
+ * Initialize and apply SHAKE constraints.
+ *
+ * \param[in] testData        Test data structure.
+ * \param[in] pbc             Periodic boundary data.
+ */
+void applyShake(ConstraintsTestData *testData, t_pbc gmx_unused pbc)
+{
+    shakedata* shaked = shake_init();
+    make_shake_sblock_serial(shaked, &testData->idef_, testData->md_);
+    bool       success = constrain_shake(
+                nullptr,
+                shaked,
+                testData->invmass_.data(),
+                testData->idef_,
+                testData->ir_,
+                as_rvec_array(testData->x_.data()),
+                as_rvec_array(testData->xPrime_.data()),
+                as_rvec_array(testData->xPrime2_.data()),
+                &testData->nrnb_,
+                testData->md_.lambda,
+                &testData->dHdLambda_,
+                testData->invdt_,
+                as_rvec_array(testData->v_.data()),
+                testData->computeVirial_,
+                testData->virialScaled_,
+                false,
+                gmx::ConstraintVariable::Positions);
+    EXPECT_TRUE(success) << "Test failed with a false return value in SHAKE.";
+    done_shake(shaked);
+}
+
+/*! \brief
+ * Initialize and apply LINCS constraints.
+ *
+ * \param[in] testData        Test data structure.
+ * \param[in] pbc             Periodic boundary data.
+ */
+void applyLincs(ConstraintsTestData *testData, t_pbc pbc)
+{
+
+    Lincs                *lincsd;
+    int                   maxwarn         = 100;
+    int                   warncount_lincs = 0;
+    gmx_omp_nthreads_set(emntLINCS, 1);
+
+    // Make blocka structure for faster LINCS setup
+    std::vector<t_blocka> at2con_mt;
+    at2con_mt.reserve(testData->mtop_.moltype.size());
+    for (const gmx_moltype_t &moltype : testData->mtop_.moltype)
+    {
+        // This function is in constr.cpp
+        at2con_mt.push_back(make_at2con(moltype,
+                                        testData->mtop_.ffparams.iparams,
+                                        flexibleConstraintTreatment(EI_DYNAMICS(testData->ir_.eI))));
+    }
+    // Initialize LINCS
+    lincsd = init_lincs(nullptr, testData->mtop_,
+                        testData->nflexcon_, at2con_mt,
+                        false,
+                        testData->ir_.nLincsIter, testData->ir_.nProjOrder);
+    set_lincs(testData->idef_, testData->md_, EI_DYNAMICS(testData->ir_.eI), &testData->cr_, lincsd);
+
+    // Evaluate constraints
+    bool success = constrain_lincs(false, testData->ir_, 0, lincsd, testData->md_,
+                                   &testData->cr_,
+                                   &testData->ms_,
+                                   as_rvec_array(testData->x_.data()),
+                                   as_rvec_array(testData->xPrime_.data()),
+                                   as_rvec_array(testData->xPrime2_.data()),
+                                   pbc.box, &pbc,
+                                   testData->md_.lambda, &testData->dHdLambda_,
+                                   testData->invdt_,
+                                   as_rvec_array(testData->v_.data()),
+                                   testData->computeVirial_, testData->virialScaled_,
+                                   gmx::ConstraintVariable::Positions,
+                                   &testData->nrnb_,
+                                   maxwarn, &warncount_lincs);
+    EXPECT_TRUE(success) << "Test failed with a false return value in LINCS.";
+    EXPECT_EQ(warncount_lincs, 0) << "There were warnings in LINCS.";
+    for (unsigned int i = 0; i < testData->mtop_.moltype.size(); i++)
+    {
+        sfree(at2con_mt.at(i).index);
+        sfree(at2con_mt.at(i).a);
+    }
+    done_lincs(lincsd);
+}
+
+#if GMX_GPU != GMX_GPU_CUDA
+/*! \brief
+ * Stub for LINCS constraints on CUDA-enabled GPU to satisfy compiler.
+ *
+ * \param[in] testData        Test data structure.
+ * \param[in] pbc             Periodic boundary data.
+ */
+void applyLincsCuda(ConstraintsTestData gmx_unused *testData, t_pbc gmx_unused pbc)
+{
+    FAIL() << "Dummy LINCS CUDA function was called instead of the real one.";
+}
+#endif
+
+} // namespace test
+} // namespace gmx
diff --git a/src/gromacs/mdlib/tests/constr_impl.cu b/src/gromacs/mdlib/tests/constr_impl.cu
new file mode 100644 (file)
index 0000000..45027bd
--- /dev/null
@@ -0,0 +1,113 @@
+/*
+ * This file is part of the GROMACS molecular simulation package.
+ *
+ * Copyright (c) 2018,2019, by the GROMACS development team, led by
+ * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
+ * and including many others, as listed in the AUTHORS file in the
+ * top-level source directory and at http://www.gromacs.org.
+ *
+ * GROMACS is free software; you can redistribute it and/or
+ * modify it under the terms of the GNU Lesser General Public License
+ * as published by the Free Software Foundation; either version 2.1
+ * of the License, or (at your option) any later version.
+ *
+ * GROMACS is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
+ * Lesser General Public License for more details.
+ *
+ * You should have received a copy of the GNU Lesser General Public
+ * License along with GROMACS; if not, see
+ * http://www.gnu.org/licenses, or write to the Free Software Foundation,
+ * Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA.
+ *
+ * If you want to redistribute modifications to GROMACS, please
+ * consider that scientific software is very special. Version
+ * control is crucial - bugs must be traceable. We will be happy to
+ * consider code for inclusion in the official distribution, but
+ * derived work must not be called official GROMACS. Details are found
+ * in the README & COPYING files - if they are missing, get the
+ * official version at http://www.gromacs.org.
+ *
+ * To help us fund GROMACS development, we humbly ask that you cite
+ * the research papers on the package. Check out http://www.gromacs.org.
+ */
+/*! \internal \file
+ * \brief Subroutines to run LINCS on GPU
+ *
+ * Copies data to GPU, runs LINCS and copies the results back.
+ *
+ * \author Artem Zhmurov <zhmurov@gmail.com>
+ * \ingroup module_mdlib
+ */
+#include "gmxpre.h"
+
+#include "constr_impl.h"
+
+#include <assert.h>
+
+#include <cmath>
+
+#include <algorithm>
+#include <unordered_map>
+#include <vector>
+
+#include "gromacs/gpu_utils/devicebuffer.cuh"
+#include "gromacs/gpu_utils/gpu_utils.h"
+#include "gromacs/math/vec.h"
+#include "gromacs/math/vectypes.h"
+#include "gromacs/mdlib/constr.h"
+#include "gromacs/mdlib/lincs_cuda.cuh"
+#include "gromacs/pbcutil/pbc.h"
+#include "gromacs/utility/unique_cptr.h"
+
+namespace gmx
+{
+namespace test
+{
+
+/*! \brief
+ * Initialize and apply LINCS constraints on CUDA-enabled GPU.
+ *
+ * \param[in] testData        Test data structure.
+ * \param[in] pbc             Periodic boundary data.
+ */
+void applyLincsCuda(ConstraintsTestData *testData, t_pbc pbc)
+{
+    auto lincsCuda = std::make_unique<LincsCuda>(testData->ir_.nLincsIter,
+                                                 testData->ir_.nProjOrder);
+
+    bool    updateVelocities = true;
+    int     numAtoms         = testData->numAtoms_;
+    float3 *d_x, *d_xp, *d_v;
+
+    lincsCuda->set(testData->idef_, testData->md_);
+    lincsCuda->setPbc(&pbc);
+
+    allocateDeviceBuffer(&d_x,  numAtoms, nullptr);
+    allocateDeviceBuffer(&d_xp, numAtoms, nullptr);
+    allocateDeviceBuffer(&d_v,  numAtoms, nullptr);
+
+    copyToDeviceBuffer(&d_x, (float3*)(testData->x_.data()), 0, numAtoms, nullptr, GpuApiCallBehavior::Sync, nullptr);
+    copyToDeviceBuffer(&d_xp, (float3*)(testData->xPrime_.data()), 0, numAtoms, nullptr, GpuApiCallBehavior::Sync, nullptr);
+    if (updateVelocities)
+    {
+        copyToDeviceBuffer(&d_v, (float3*)(testData->v_.data()), 0, numAtoms, nullptr, GpuApiCallBehavior::Sync, nullptr);
+    }
+    lincsCuda->apply(d_x, d_xp,
+                     updateVelocities, d_v, testData->invdt_,
+                     testData->computeVirial_, testData->virialScaled_);
+
+    copyFromDeviceBuffer((float3*)(testData->xPrime_.data()), &d_xp, 0, numAtoms, nullptr, GpuApiCallBehavior::Sync, nullptr);
+    if (updateVelocities)
+    {
+        copyFromDeviceBuffer((float3*)(testData->v_.data()), &d_v, 0, numAtoms, nullptr, GpuApiCallBehavior::Sync, nullptr);
+    }
+
+    freeDeviceBuffer(&d_x);
+    freeDeviceBuffer(&d_xp);
+    freeDeviceBuffer(&d_v);
+}
+
+} // namespace test
+} // namespace gmx
diff --git a/src/gromacs/mdlib/tests/constr_impl.h b/src/gromacs/mdlib/tests/constr_impl.h
new file mode 100644 (file)
index 0000000..c46d00e
--- /dev/null
@@ -0,0 +1,237 @@
+/*
+ * This file is part of the GROMACS molecular simulation package.
+ *
+ * Copyright (c) 2018,2019, by the GROMACS development team, led by
+ * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
+ * and including many others, as listed in the AUTHORS file in the
+ * top-level source directory and at http://www.gromacs.org.
+ *
+ * GROMACS is free software; you can redistribute it and/or
+ * modify it under the terms of the GNU Lesser General Public License
+ * as published by the Free Software Foundation; either version 2.1
+ * of the License, or (at your option) any later version.
+ *
+ * GROMACS is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
+ * Lesser General Public License for more details.
+ *
+ * You should have received a copy of the GNU Lesser General Public
+ * License along with GROMACS; if not, see
+ * http://www.gnu.org/licenses, or write to the Free Software Foundation,
+ * Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA.
+ *
+ * If you want to redistribute modifications to GROMACS, please
+ * consider that scientific software is very special. Version
+ * control is crucial - bugs must be traceable. We will be happy to
+ * consider code for inclusion in the official distribution, but
+ * derived work must not be called official GROMACS. Details are found
+ * in the README & COPYING files - if they are missing, get the
+ * official version at http://www.gromacs.org.
+ *
+ * To help us fund GROMACS development, we humbly ask that you cite
+ * the research papers on the package. Check out http://www.gromacs.org.
+ */
+/*! \internal \file
+ * \brief SHAKE and LINCS tests header.
+ *
+ * Contains description and constructor for the test data accumulating object,
+ * declares CPU- and GPU-based functions used to apply SHAKE or LINCS on the
+ * test data.
+ *
+ * \author Artem Zhmurov <zhmurov@gmail.com>
+ * \ingroup module_mdlib
+ */
+
+#ifndef GMX_MDLIB_TESTS_CONSTR_IMPL_H
+#define GMX_MDLIB_TESTS_CONSTR_IMPL_H
+
+#include "gmxpre.h"
+
+#include <assert.h>
+
+#include <cmath>
+
+#include <algorithm>
+#include <unordered_map>
+#include <vector>
+
+#include "gromacs/fileio/gmxfio.h"
+#include "gromacs/gmxlib/nrnb.h"
+#include "gromacs/gmxlib/nonbonded/nonbonded.h"
+#include "gromacs/gpu_utils/gpu_utils.h"
+#include "gromacs/math/paddedvector.h"
+#include "gromacs/math/vec.h"
+#include "gromacs/math/vectypes.h"
+#include "gromacs/mdlib/constr.h"
+#include "gromacs/mdlib/gmx_omp_nthreads.h"
+#include "gromacs/mdlib/lincs.h"
+#include "gromacs/mdlib/shake.h"
+#include "gromacs/mdrunutility/multisim.h"
+#include "gromacs/mdtypes/commrec.h"
+#include "gromacs/mdtypes/inputrec.h"
+#include "gromacs/mdtypes/mdatom.h"
+#include "gromacs/pbcutil/pbc.h"
+#include "gromacs/topology/block.h"
+#include "gromacs/topology/idef.h"
+#include "gromacs/topology/ifunc.h"
+#include "gromacs/topology/topology.h"
+#include "gromacs/utility/smalloc.h"
+#include "gromacs/utility/stringutil.h"
+#include "gromacs/utility/unique_cptr.h"
+
+namespace gmx
+{
+namespace test
+{
+
+/* \brief
+ * Constraints test data structure.
+ *
+ * Structure to collect all the necessary data, including system coordinates and topology,
+ * constraints information, etc. The structure can be reset and reused.
+ */
+class ConstraintsTestData
+{
+    public:
+        //! Human-friendly name for a system
+        std::string           title_;
+        //! Number of atoms
+        int                   numAtoms_;
+        //! Topology
+        gmx_mtop_t            mtop_;
+        //! Masses
+        std::vector<real>     masses_;
+        //! Inverse masses
+        std::vector<real>     invmass_;
+        //! Communication record
+        t_commrec             cr_;
+        //! Input record (info that usually in .mdp file)
+        t_inputrec            ir_;
+        //! Local topology
+        t_idef                idef_;
+        //! MD atoms
+        t_mdatoms             md_;
+        //! Multisim data
+        gmx_multisim_t        ms_;
+        //! Computational time array (normally used to benchmark performance)
+        t_nrnb                nrnb_;
+
+        //! Inverse timestep
+        real                  invdt_;
+        //! Number of flexible constraints
+        int                   nflexcon_   = 0;
+        //! Whether the virial should be computed
+        bool                  computeVirial_;
+        //! Scaled virial
+        tensor                virialScaled_;
+        //! Scaled virial (reference values)
+        tensor                virialScaledRef_;
+        //! If the free energy is computed
+        bool                  compute_dHdLambda_;
+        //! For free energy computation
+        real                  dHdLambda_;
+        //! For free energy computation (reference value)
+        real                  dHdLambdaRef_;
+
+        //! Coordinates before the timestep
+        PaddedVector<RVec>    x_;
+        //! Coordinates after timestep, output for the constraints
+        PaddedVector<RVec>    xPrime_;
+        //! Backup for coordinates (for reset)
+        PaddedVector<RVec>    xPrime0_;
+        //! Intermediate set of coordinates (normally used for projection correction)
+        PaddedVector<RVec>    xPrime2_;
+        //! Velocities
+        PaddedVector<RVec>    v_;
+        //! Backup for velocities (for reset)
+        PaddedVector<RVec>    v0_;
+
+        //! Constraints data (type1-i1-j1-type2-i2-j2-...)
+        std::vector<int>      constraints_;
+        //! Target lengths for all constraint types
+        std::vector<real>     constraintsR0_;
+
+        /*! \brief
+         * Constructor for the object with all parameters and variables needed by constraints algorithms.
+         *
+         * This constructor assembles stubs for all the data structures, required to initialize
+         * and apply LINCS and SHAKE constraints. The coordinates and velocities before constraining
+         * are saved to allow for reset. The constraints data are stored for testing after constraints
+         * were applied.
+         *
+         * \param[in]  title                Human-friendly name of the system.
+         * \param[in]  numAtoms             Number of atoms in the system.
+         * \param[in]  masses               Atom masses. Size of this vector should be equal to numAtoms.
+         * \param[in]  constraints          List of constraints, organized in triples of integers.
+         *                                  First integer is the index of type for a constraint, second
+         *                                  and third are the indices of constrained atoms. The types
+         *                                  of constraints should be sequential but not necessarily
+         *                                  start from zero (which is the way they normally are in
+         *                                  GROMACS).
+         * \param[in]  constraintsR0        Target values for bond lengths for bonds of each type. The
+         *                                  size of this vector should be equal to the total number of
+         *                                  unique types in constraints vector.
+         * \param[in]  computeVirial        Whether the virial should be computed.
+         * \param[in]  virialScaledRef      Reference values for scaled virial tensor.
+         * \param[in]  compute_dHdLambda    Whether free energy should be computed.
+         * \param[in]  dHdLambdaRef         Reference value for dHdLambda.
+         * \param[in]  initialTime          Initial time.
+         * \param[in]  timestep             Timestep.
+         * \param[in]  x                    Coordinates before integration step.
+         * \param[in]  xPrime               Coordinates after integration step, but before constraining.
+         * \param[in]  v                    Velocities before constraining.
+         * \param[in]  shakeTolerance       Target tolerance for SHAKE.
+         * \param[in]  shakeUseSOR          Use successive over-relaxation method for SHAKE iterations.
+         *                                  The general formula is:
+         *                                     x_n+1 = (1-omega)*x_n + omega*f(x_n),
+         *                                  where omega = 1 if SOR is off and may be < 1 if SOR is on.
+         * \param[in]  lincsNumIterations   Number of iterations used to compute the inverse matrix.
+         * \param[in]  lincsExpansionOrder  The order for algorithm that adjusts the direction of the
+         *                                  bond after constraints are applied.
+         * \param[in]  lincsWarnAngle       The threshold value for the change in bond angle. When
+         *                                  exceeded the program will issue a warning.
+         *
+         */
+        ConstraintsTestData(const std::string &title,
+                            int numAtoms, std::vector<real> masses,
+                            std::vector<int> constraints, std::vector<real> constraintsR0,
+                            bool computeVirial, tensor virialScaledRef,
+                            bool compute_dHdLambda, float dHdLambdaRef,
+                            real initialTime, real timestep,
+                            const std::vector<RVec> &x, const std::vector<RVec> &xPrime, const std::vector<RVec> &v,
+                            real shakeTolerance, gmx_bool shakeUseSOR,
+                            int lincsNumIterations, int lincsExpansionOrder, real lincsWarnAngle);
+
+        /*! \brief
+         * Reset the data structure so it can be reused.
+         *
+         * Set the coordinates and velocities back to their values before
+         * constraining. The scaled virial tensor and dHdLambda are zeroed.
+         *
+         */
+        void reset();
+
+        /*! \brief
+         * Cleaning up the memory.
+         */
+        ~ConstraintsTestData();
+};
+
+/*! \brief Apply SHAKE constraints to the test data.
+ */
+void applyShake(ConstraintsTestData *testData, t_pbc pbc);
+/*! \brief Apply LINCS constraints to the test data.
+ */
+void applyLincs(ConstraintsTestData *testData, t_pbc pbc);
+/*! \brief Apply CUDA version of LINCS constraints to the test data.
+ *
+ * All the data is copied to the GPU device, then LINCS is applied and
+ * the resulting coordinates are copied back.
+ */
+void applyLincsCuda(ConstraintsTestData *testData, t_pbc pbc);
+
+}      // namespace test
+}      // namespace gmx
+
+#endif // GMX_MDLIB_TESTS_CONSTR_IMPL_H
index e1ec5e40222f87f69b83aeb4996cf29616dcaccb..4f66f1260e1c123a34e5797ba21fcf0a760e67ae 100644 (file)
 #include "gromacs/gpu_utils/gputraits.cuh"
 #include "gromacs/gpu_utils/vectype_ops.cuh"
 #include "gromacs/math/vec.h"
+#include "gromacs/mdlib/lincs_cuda.cuh"
 #include "gromacs/mdlib/update_constrain_cuda.h"
 #include "gromacs/pbcutil/pbc.h"
 #include "gromacs/pbcutil/pbc_aiuc_cuda.cuh"
 
 #include "leapfrog_cuda_impl.h"
-#include "lincs_cuda_impl.h"
 #include "settle_cuda_impl.h"
 
 namespace gmx
@@ -121,7 +121,7 @@ UpdateConstrainCuda::Impl::Impl(int                numAtoms,
 
     GMX_RELEASE_ASSERT(numAtoms == mtop.natoms, "State and topology number of atoms should be the same.");
     integrator_ = std::make_unique<LeapFrogCuda::Impl>();
-    lincsCuda_  = std::make_unique<LincsCuda::Impl>(ir.nLincsIter, ir.nProjOrder);
+    lincsCuda_  = std::make_unique<LincsCuda>(ir.nLincsIter, ir.nProjOrder);
     settleCuda_ = std::make_unique<SettleCuda::Impl>(mtop);
 
 }
index cd360b901f446ded637af11db5b70b67de38bc41..1e6869c21b58410f24814bfc4294245f40c4b39a 100644 (file)
@@ -46,6 +46,7 @@
 #ifndef GMX_MDLIB_UPDATE_CONSTRAIN_CUDA_IMPL_H
 #define GMX_MDLIB_UPDATE_CONSTRAIN_CUDA_IMPL_H
 
+#include "gromacs/mdlib/lincs_cuda.cuh"
 #include "gromacs/mdlib/update_constrain_cuda.h"
 #include "gromacs/mdtypes/inputrec.h"
 #include "gromacs/pbcutil/pbc.h"
@@ -53,7 +54,6 @@
 #include "gromacs/topology/idef.h"
 
 #include "leapfrog_cuda_impl.h"
-#include "lincs_cuda_impl.h"
 #include "settle_cuda_impl.h"
 
 namespace gmx
@@ -205,7 +205,7 @@ class UpdateConstrainCuda::Impl
         //! Leap-Frog integrator
         std::unique_ptr<LeapFrogCuda::Impl>  integrator_;
         //! LINCS CUDA object to use for non-water constraints
-        std::unique_ptr<LincsCuda::Impl>     lincsCuda_;
+        std::unique_ptr<LincsCuda>           lincsCuda_;
         //! SETTLE CUDA object for water constrains
         std::unique_ptr<SettleCuda::Impl>    settleCuda_;
 
index 32003071f99f100f2acf14ae3752582e4c7c7f10..72811bcdf91151a9769c11dbd25153b1efccd85a 100644 (file)
@@ -96,9 +96,10 @@ struct PbcAiuc
  * \param[out]  pbcAiuc     Pointer to PbcAiuc structure to be initialized.
  *
  */
-static void setPbcAiuc(int           numPbcDim,
-                       const matrix  box,
-                       PbcAiuc      *pbcAiuc)
+static inline
+void setPbcAiuc(int           numPbcDim,
+                const matrix  box,
+                PbcAiuc      *pbcAiuc)
 {
 
     pbcAiuc->invBoxDiagZ = 0.0f;