From 7dd70297e3473d978bd3ca22882c0336771ed20b Mon Sep 17 00:00:00 2001 From: Joe Jordan Date: Thu, 7 Oct 2021 07:35:17 +0000 Subject: [PATCH] Introduce nblib nonbonded force calculator impl class --- api/nblib/CMakeLists.txt | 19 +- .../{gmxcalculator.h => gmxbackenddata.h} | 103 +++--- api/nblib/gmxcalculator.cpp | 119 ------ api/nblib/gmxcalculatorcpu.cpp | 339 ++++++++++++++++++ api/nblib/gmxcalculatorcpu.h | 112 ++++++ api/nblib/nblib.h | 4 +- api/nblib/nbnxmsetuphelpers.cpp | 2 +- api/nblib/nbnxmsetuphelpers.h | 2 +- .../samples/argon-forces-integration.cpp | 7 +- .../samples/methane-water-integration.cpp | 11 +- api/nblib/systemdescription.h | 130 +++++++ api/nblib/tests/gmxcalculator.cpp | 72 +++- api/nblib/tests/nbkernelsystem.cpp | 74 ++-- .../NBlibTest_ArgonEnergiesAreCorrect.xml | 12 + .../NBlibTest_ArgonVirialsAreCorrect.xml | 16 + 15 files changed, 788 insertions(+), 234 deletions(-) rename api/nblib/{gmxcalculator.h => gmxbackenddata.h} (57%) delete mode 100644 api/nblib/gmxcalculator.cpp create mode 100644 api/nblib/gmxcalculatorcpu.cpp create mode 100644 api/nblib/gmxcalculatorcpu.h create mode 100644 api/nblib/systemdescription.h create mode 100644 api/nblib/tests/refdata/NBlibTest_ArgonEnergiesAreCorrect.xml create mode 100644 api/nblib/tests/refdata/NBlibTest_ArgonVirialsAreCorrect.xml diff --git a/api/nblib/CMakeLists.txt b/api/nblib/CMakeLists.txt index 31eca82e77..3619ddb510 100644 --- a/api/nblib/CMakeLists.txt +++ b/api/nblib/CMakeLists.txt @@ -52,14 +52,14 @@ add_custom_target(nblib-tests # Ensure that "make tests" builds all nblib tests so the top-level # "make check" will work. if (BUILD_TESTING) - add_dependencies(tests nblib-tests) + add_dependencies(tests nblib-tests) - # this allows all nblib tests to be run with "make check-nblib" - add_custom_target(check-nblib - COMMAND ${CMAKE_CTEST_COMMAND} --output-on-failure -R NbLib - COMMENT "Running nblib tests" - USES_TERMINAL VERBATIM) - add_dependencies(check-nblib nblib-tests) + # this allows all nblib tests to be run with "make check-nblib" + add_custom_target(check-nblib + COMMAND ${CMAKE_CTEST_COMMAND} --output-on-failure -R NbLib + COMMENT "Running nblib tests" + USES_TERMINAL VERBATIM) + add_dependencies(check-nblib nblib-tests) endif() set(NBLIB_MAJOR 0) @@ -89,8 +89,7 @@ endif() target_sources(nblib PRIVATE box.cpp - gmxcalculator.cpp - gmxsetup.cpp + gmxcalculatorcpu.cpp integrator.cpp interactions.cpp molecules.cpp @@ -131,7 +130,7 @@ if(GMX_INSTALL_NBLIB_API) basicdefinitions.h box.h exception.h - gmxcalculator.h + gmxcalculatorcpu.h integrator.h interactions.h molecules.h diff --git a/api/nblib/gmxcalculator.h b/api/nblib/gmxbackenddata.h similarity index 57% rename from api/nblib/gmxcalculator.h rename to api/nblib/gmxbackenddata.h index 71a56ebb08..6794bb47de 100644 --- a/api/nblib/gmxcalculator.h +++ b/api/nblib/gmxbackenddata.h @@ -1,7 +1,7 @@ /* * This file is part of the GROMACS molecular simulation package. * - * Copyright (c) 2020,2021, by the GROMACS development team, led by + * Copyright (c) 2021, 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. @@ -32,44 +32,32 @@ * 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 +/*! \inpublicapi \file * \brief - * Implements a force calculator based on GROMACS data structures. - * - * Intended for internal use inside the ForceCalculator. + * Implements nblib simulation box * * \author Victor Holanda * \author Joe Jordan * \author Prashanth Kanduri * \author Sebastian Keller */ - -#ifndef NBLIB_GMXCALCULATOR_H -#define NBLIB_GMXCALCULATOR_H - -#include - -#include "nblib/vector.h" - -struct nonbonded_verlet_t; -struct t_forcerec; -struct t_nrnb; -struct interaction_const_t; -struct gmx_enerdata_t; - -namespace gmx -{ -template -class ArrayRef; -class StepWorkload; -} // namespace gmx +#ifndef NBLIB_GMXBACKENDDATA_H +#define NBLIB_GMXBACKENDDATA_H + +#include "gromacs/gmxlib/nrnb.h" +#include "gromacs/mdtypes/enerdata.h" +#include "gromacs/mdtypes/forcerec.h" +#include "gromacs/mdtypes/interaction_const.h" +#include "gromacs/mdtypes/simulation_workload.h" +#include "gromacs/nbnxm/nbnxm.h" +#include "gromacs/timing/wallcycle.h" +#include "gromacs/utility/listoflists.h" +#include "gromacs/utility/range.h" +#include "nblib/kerneloptions.h" +#include "nblib/nbnxmsetuphelpers.h" namespace nblib { -class Box; -class NbvSetupUtil; -class SimulationState; -struct NBKernelOptions; /*! \brief GROMACS non-bonded force calculation backend * @@ -82,48 +70,59 @@ struct NBKernelOptions; * interactions (PP) on a single node CPU. * */ - -class GmxForceCalculator final +class GmxBackendData { public: - GmxForceCalculator(); + GmxBackendData() = default; + GmxBackendData(const NBKernelOptions& options, + int numEnergyGroups, + gmx::ArrayRef exclusionRanges, + gmx::ArrayRef exclusionElements) + { + // Set hardware params from the execution context + setGmxNonBondedNThreads(options.numOpenMPThreads); + + // Set interaction constants struct + interactionConst_ = createInteractionConst(options); - ~GmxForceCalculator(); + // Set up StepWorkload data + stepWork_ = createStepWorkload(); - //! Compute forces and return - void compute(gmx::ArrayRef coordinateInput, gmx::ArrayRef forceOutput); + // Set up gmx_enerdata_t (holds energy information) + enerd_ = gmx_enerdata_t{ numEnergyGroups, 0 }; - //! Puts particles on a grid based on bounds specified by the box (for every NS step) - void setParticlesOnGrid(gmx::ArrayRef particleInfoAllVdw, - gmx::ArrayRef coordinates, - const Box& box); + // Construct pair lists + std::vector exclusionRanges_(exclusionRanges.begin(), exclusionRanges.end()); + std::vector exclusionElements_(exclusionElements.begin(), exclusionElements.end()); + exclusions_ = gmx::ListOfLists(std::move(exclusionRanges_), std::move(exclusionElements_)); + } -private: - //! Friend to allow setting up private members in this class - friend class NbvSetupUtil; + //! exclusions in gmx format + gmx::ListOfLists exclusions_; //! Non-Bonded Verlet object for force calculation std::unique_ptr nbv_; - //! Only nbfp and shift_vec are used - std::unique_ptr forcerec_; + //! Only shift_vec is used + t_forcerec forcerec_; //! Parameters for various interactions in the system - std::unique_ptr interactionConst_; + interaction_const_t interactionConst_; //! Tasks to perform in an MD Step - std::unique_ptr stepWork_; + gmx::StepWorkload stepWork_; + + gmx::SimulationWorkload simulationWork_; //! Energies of different interaction types; currently only needed as an argument for dispatchNonbondedKernel - std::unique_ptr enerd_; + gmx_enerdata_t enerd_{ 1, 0 }; //! Non-bonded flop counter; currently only needed as an argument for dispatchNonbondedKernel - std::unique_ptr nrnb_; + t_nrnb nrnb_; - //! Legacy matrix for box - matrix box_{ { 0 } }; + //! Keep track of whether updatePairlist has been called at least once + bool updatePairlistCalled{ false }; }; } // namespace nblib - -#endif // NBLIB_GMXCALCULATOR_H +#endif // NBLIB_GMXBACKENDDATA_H diff --git a/api/nblib/gmxcalculator.cpp b/api/nblib/gmxcalculator.cpp deleted file mode 100644 index 81e740f7c0..0000000000 --- a/api/nblib/gmxcalculator.cpp +++ /dev/null @@ -1,119 +0,0 @@ -/* - * This file is part of the GROMACS molecular simulation package. - * - * Copyright (c) 2020,2021, 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 Implements a force calculator based on GROMACS data structures. - * - * \author Victor Holanda - * \author Joe Jordan - * \author Prashanth Kanduri - * \author Sebastian Keller - */ -#include "nblib/gmxcalculator.h" -#include "gromacs/ewald/ewald_utils.h" -#include "gromacs/gmxlib/nrnb.h" -#include "gromacs/mdlib/rf_util.h" -#include "gromacs/mdtypes/enerdata.h" -#include "gromacs/mdtypes/forcerec.h" -#include "gromacs/mdtypes/interaction_const.h" -#include "gromacs/mdtypes/simulation_workload.h" -#include "gromacs/nbnxm/nbnxm.h" -#include "gromacs/utility/range.h" -#include "nblib/exception.h" -#include "nblib/simulationstate.h" - -namespace nblib -{ - -GmxForceCalculator::GmxForceCalculator() -{ - enerd_ = std::make_unique(1, 0); - forcerec_ = std::make_unique(); - interactionConst_ = std::make_unique(); - stepWork_ = std::make_unique(); - nrnb_ = std::make_unique(); -} - -GmxForceCalculator::~GmxForceCalculator() = default; - -void GmxForceCalculator::compute(gmx::ArrayRef coordinateInput, - gmx::ArrayRef forceOutput) -{ - // update the coordinates in the backend - nbv_->convertCoordinates(gmx::AtomLocality::Local, coordinateInput); - - nbv_->dispatchNonbondedKernel( - gmx::InteractionLocality::Local, - *interactionConst_, - *stepWork_, - enbvClearFYes, - forcerec_->shift_vec, - enerd_->grpp.energyGroupPairTerms[forcerec_->haveBuckingham ? NonBondedEnergyTerms::BuckinghamSR - : NonBondedEnergyTerms::LJSR], - enerd_->grpp.energyGroupPairTerms[NonBondedEnergyTerms::CoulombSR], - nrnb_.get()); - - nbv_->atomdata_add_nbat_f_to_f(gmx::AtomLocality::All, forceOutput); -} - -void GmxForceCalculator::setParticlesOnGrid(gmx::ArrayRef particleInfoAllVdw, - gmx::ArrayRef coordinates, - const Box& box) -{ - auto legacyBox = box.legacyMatrix(); - - if (TRICLINIC(legacyBox)) - { - throw InputException("Only rectangular unit-cells are supported here"); - } - const rvec lowerCorner = { 0, 0, 0 }; - const rvec upperCorner = { legacyBox[dimX][dimX], legacyBox[dimY][dimY], legacyBox[dimZ][dimZ] }; - - const real particleDensity = coordinates.size() / det(legacyBox); - - nbnxn_put_on_grid(nbv_.get(), - legacyBox, - 0, - lowerCorner, - upperCorner, - nullptr, - { 0, int(coordinates.size()) }, - particleDensity, - particleInfoAllVdw, - coordinates, - 0, - nullptr); -} - -} // namespace nblib diff --git a/api/nblib/gmxcalculatorcpu.cpp b/api/nblib/gmxcalculatorcpu.cpp new file mode 100644 index 0000000000..a28a376065 --- /dev/null +++ b/api/nblib/gmxcalculatorcpu.cpp @@ -0,0 +1,339 @@ +/* + * This file is part of the GROMACS molecular simulation package. + * + * Copyright (c) 2020,2021, 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 Implements a force calculator based on GROMACS data structures. + * + * \author Victor Holanda + * \author Joe Jordan + * \author Prashanth Kanduri + * \author Sebastian Keller + */ +#include "gromacs/ewald/ewald_utils.h" +#include "gromacs/mdtypes/enerdata.h" +#include "gromacs/mdtypes/interaction_const.h" +#include "gromacs/nbnxm/atomdata.h" +#include "gromacs/nbnxm/nbnxm.h" +#include "gromacs/nbnxm/pairlistset.h" +#include "gromacs/nbnxm/pairlistsets.h" +#include "gromacs/nbnxm/pairsearch.h" +#include "gromacs/utility/listoflists.h" +#include "gromacs/utility/range.h" +#include "nblib/exception.h" +#include "nblib/gmxbackenddata.h" +#include "nblib/gmxcalculatorcpu.h" +#include "nblib/nbnxmsetuphelpers.h" +#include "nblib/pbc.hpp" +#include "nblib/systemdescription.h" +#include "nblib/topology.h" +#include "nblib/virials.h" + +namespace nblib +{ + +class GmxNBForceCalculatorCpu::CpuImpl final +{ +public: + CpuImpl(gmx::ArrayRef particleTypeIdOfAllParticles, + gmx::ArrayRef nonBondedParams, + gmx::ArrayRef charges, + gmx::ArrayRef particleInteractionFlags, + gmx::ArrayRef exclusionRanges, + gmx::ArrayRef exclusionElements, + const NBKernelOptions& options); + + //! calculates a new pair list based on new coordinates (for every NS step) + void updatePairlist(gmx::ArrayRef coordinates, const Box& box); + + //! Compute forces and return + void compute(gmx::ArrayRef coordinateInput, + const Box& box, + gmx::ArrayRef forceOutput); + + //! Compute forces and virial tensor + void compute(gmx::ArrayRef coordinateInput, + const Box& box, + gmx::ArrayRef forceOutput, + gmx::ArrayRef virialOutput); + + //! Compute forces, virial tensor and potential energies + void compute(gmx::ArrayRef coordinateInput, + const Box& box, + gmx::ArrayRef forceOutput, + gmx::ArrayRef virialOutput, + gmx::ArrayRef energyOutput); + +private: + //! \brief client-side provided system description data + SystemDescription system_; + + //! \brief Gmx backend objects, employed for calculating the forces + GmxBackendData backend_; +}; + +GmxNBForceCalculatorCpu::CpuImpl::CpuImpl(gmx::ArrayRef particleTypeIdOfAllParticles, + gmx::ArrayRef nonBondedParams, + gmx::ArrayRef charges, + gmx::ArrayRef particleInteractionFlags, + gmx::ArrayRef exclusionRanges, + gmx::ArrayRef exclusionElements, + const NBKernelOptions& options) : + system_(SystemDescription(particleTypeIdOfAllParticles, nonBondedParams, charges, particleInteractionFlags)), + backend_(GmxBackendData(options, findNumEnergyGroups(particleInteractionFlags), exclusionRanges, exclusionElements)) +{ + // Set up non-bonded verlet in the backend + backend_.nbv_ = createNbnxmCPU(system_.numParticleTypes_, + options, + findNumEnergyGroups(particleInteractionFlags), + system_.nonBondedParams_); +} + +void GmxNBForceCalculatorCpu::CpuImpl::updatePairlist(gmx::ArrayRef coordinates, + const Box& box) +{ + if (coordinates.size() != system_.numParticles_) + { + throw InputException( + "Coordinate array containing different number of entries than particles in the " + "system"); + } + + const auto* legacyBox = box.legacyMatrix(); + system_.box_ = box; + updateForcerec(&backend_.forcerec_, box.legacyMatrix()); + if (TRICLINIC(legacyBox)) + { + throw InputException("Only rectangular unit-cells are supported here"); + } + + const rvec lowerCorner = { 0, 0, 0 }; + const rvec upperCorner = { legacyBox[dimX][dimX], legacyBox[dimY][dimY], legacyBox[dimZ][dimZ] }; + + const real particleDensity = static_cast(coordinates.size()) / det(legacyBox); + + // Put particles on a grid based on bounds specified by the box + nbnxn_put_on_grid(backend_.nbv_.get(), + legacyBox, + 0, + lowerCorner, + upperCorner, + nullptr, + { 0, int(coordinates.size()) }, + particleDensity, + system_.particleInfo_, + coordinates, + 0, + nullptr); + + backend_.nbv_->constructPairlist( + gmx::InteractionLocality::Local, backend_.exclusions_, 0, &backend_.nrnb_); + + // Set Particle Types and Charges and VdW params + backend_.nbv_->setAtomProperties( + system_.particleTypeIdOfAllParticles_, system_.charges_, system_.particleInfo_); + backend_.updatePairlistCalled = true; +} + +void GmxNBForceCalculatorCpu::CpuImpl::compute(gmx::ArrayRef coordinateInput, + const Box& box, + gmx::ArrayRef forceOutput, + gmx::ArrayRef virialOutput, + gmx::ArrayRef energyOutput) +{ + if (coordinateInput.size() != forceOutput.size()) + { + throw InputException("coordinate array and force buffer size mismatch"); + } + + if (!backend_.updatePairlistCalled) + { + throw InputException("compute called without updating pairlist at least once"); + } + + // update the box if changed + if (!(system_.box_ == box)) + { + system_.box_ = box; + updateForcerec(&backend_.forcerec_, box.legacyMatrix()); + } + + bool computeVirial = !virialOutput.empty(); + bool computeEnergies = !energyOutput.empty(); + backend_.stepWork_.computeVirial = computeVirial; + backend_.stepWork_.computeEnergy = computeEnergies; + + // update the coordinates in the backend + backend_.nbv_->convertCoordinates(gmx::AtomLocality::Local, coordinateInput); + + backend_.nbv_->dispatchNonbondedKernel( + gmx::InteractionLocality::Local, + backend_.interactionConst_, + backend_.stepWork_, + enbvClearFYes, + backend_.forcerec_.shift_vec, + backend_.enerd_.grpp.energyGroupPairTerms[backend_.forcerec_.haveBuckingham ? NonBondedEnergyTerms::BuckinghamSR + : NonBondedEnergyTerms::LJSR], + backend_.enerd_.grpp.energyGroupPairTerms[NonBondedEnergyTerms::CoulombSR], + &backend_.nrnb_); + + backend_.nbv_->atomdata_add_nbat_f_to_f(gmx::AtomLocality::All, forceOutput); + + if (computeVirial) + { + // calculate shift forces and turn into an array ref + std::vector shiftForcesVector(gmx::c_numShiftVectors, Vec3(0.0, 0.0, 0.0)); + nbnxn_atomdata_add_nbat_fshift_to_fshift(*backend_.nbv_->nbat, shiftForcesVector); + auto shiftForcesRef = constArrayRefFromArray(shiftForcesVector.data(), shiftForcesVector.size()); + + std::vector shiftVectorsArray(gmx::c_numShiftVectors); + + // copy shift vectors from ForceRec + std::copy(backend_.forcerec_.shift_vec.begin(), + backend_.forcerec_.shift_vec.end(), + shiftVectorsArray.begin()); + + computeVirialTensor( + coordinateInput, forceOutput, shiftVectorsArray, shiftForcesRef, box, virialOutput); + } + + // extract term energies (per interaction type) + if (computeEnergies) + { + int nGroupPairs = backend_.enerd_.grpp.nener; + if (int(energyOutput.size()) != int(NonBondedEnergyTerms::Count) * nGroupPairs) + { + throw InputException("Array size for energy output is wrong\n"); + } + + for (int eg = 0; eg < int(NonBondedEnergyTerms::Count); ++eg) + { + std::copy(begin(backend_.enerd_.grpp.energyGroupPairTerms[eg]), + end(backend_.enerd_.grpp.energyGroupPairTerms[eg]), + energyOutput.begin() + eg * nGroupPairs); + } + } +} + +void GmxNBForceCalculatorCpu::CpuImpl::compute(gmx::ArrayRef coordinateInput, + const Box& box, + gmx::ArrayRef forceOutput) +{ + // compute forces and fill in force buffer + compute(coordinateInput, box, forceOutput, gmx::ArrayRef{}, gmx::ArrayRef{}); +} + +void GmxNBForceCalculatorCpu::CpuImpl::compute(gmx::ArrayRef coordinateInput, + const Box& box, + gmx::ArrayRef forceOutput, + gmx::ArrayRef virialOutput) +{ + // compute forces and fill in force buffer + compute(coordinateInput, box, forceOutput, virialOutput, gmx::ArrayRef{}); +} + +GmxNBForceCalculatorCpu::GmxNBForceCalculatorCpu(gmx::ArrayRef particleTypeIdOfAllParticles, + gmx::ArrayRef nonBondedParams, + gmx::ArrayRef charges, + gmx::ArrayRef particleInteractionFlags, + gmx::ArrayRef exclusionRanges, + gmx::ArrayRef exclusionElements, + const NBKernelOptions& options) +{ + if (options.useGpu) + { + throw InputException("Use GmxNBForceCalculatorGpu for GPU support"); + } + + impl_ = std::make_unique(particleTypeIdOfAllParticles, + nonBondedParams, + charges, + particleInteractionFlags, + exclusionRanges, + exclusionElements, + options); +} + +GmxNBForceCalculatorCpu::~GmxNBForceCalculatorCpu() = default; + +//! calculates a new pair list based on new coordinates (for every NS step) +void GmxNBForceCalculatorCpu::updatePairlist(gmx::ArrayRef coordinates, const Box& box) +{ + impl_->updatePairlist(coordinates, box); +} + +//! Compute forces and return +void GmxNBForceCalculatorCpu::compute(gmx::ArrayRef coordinateInput, + const Box& box, + gmx::ArrayRef forceOutput) +{ + impl_->compute(coordinateInput, box, forceOutput); +} + +//! Compute forces and virial tensor +void GmxNBForceCalculatorCpu::compute(gmx::ArrayRef coordinateInput, + const Box& box, + gmx::ArrayRef forceOutput, + gmx::ArrayRef virialOutput) +{ + impl_->compute(coordinateInput, box, forceOutput, virialOutput); +} + +//! Compute forces, virial tensor and potential energies +void GmxNBForceCalculatorCpu::compute(gmx::ArrayRef coordinateInput, + const Box& box, + gmx::ArrayRef forceOutput, + gmx::ArrayRef virialOutput, + gmx::ArrayRef energyOutput) +{ + impl_->compute(coordinateInput, box, forceOutput, virialOutput, energyOutput); +} + +std::unique_ptr setupGmxForceCalculatorCpu(const Topology& topology, + const NBKernelOptions& options) +{ + std::vector nonBondedParameters = createNonBondedParameters( + topology.getParticleTypes(), topology.getNonBondedInteractionMap()); + + std::vector particleInteractionFlags = createParticleInfoAllVdw(topology.numParticles()); + + return std::make_unique(topology.getParticleTypeIdOfAllParticles(), + nonBondedParameters, + topology.getCharges(), + particleInteractionFlags, + topology.exclusionLists().ListRanges, + topology.exclusionLists().ListElements, + options); +} + +} // namespace nblib diff --git a/api/nblib/gmxcalculatorcpu.h b/api/nblib/gmxcalculatorcpu.h new file mode 100644 index 0000000000..2fec548256 --- /dev/null +++ b/api/nblib/gmxcalculatorcpu.h @@ -0,0 +1,112 @@ +/* + * This file is part of the GROMACS molecular simulation package. + * + * Copyright (c) 2020,2021, 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 + * Implements a force calculator based on GROMACS data structures. + * + * Intended for internal use inside the ForceCalculator. + * + * \author Victor Holanda + * \author Joe Jordan + * \author Prashanth Kanduri + * \author Sebastian Keller + */ + +#ifndef NBLIB_GMXCALCULATORCPU_H +#define NBLIB_GMXCALCULATORCPU_H + +#include +#include + +#include "nblib/box.h" +#include "nblib/vector.h" + +namespace gmx +{ +template +class ArrayRef; +} // namespace gmx + +namespace nblib +{ +struct NBKernelOptions; +class Topology; + +class GmxNBForceCalculatorCpu final +{ +public: + GmxNBForceCalculatorCpu(gmx::ArrayRef particleTypeIdOfAllParticles, + gmx::ArrayRef nonBondedParams, + gmx::ArrayRef charges, + gmx::ArrayRef particleInteractionFlags, + gmx::ArrayRef exclusionRanges, + gmx::ArrayRef exclusionElements, + const NBKernelOptions& options); + + ~GmxNBForceCalculatorCpu(); + + //! calculates a new pair list based on new coordinates (for every NS step) + void updatePairlist(gmx::ArrayRef coordinates, const Box& box); + + //! Compute forces and return + void compute(gmx::ArrayRef coordinateInput, + const Box& box, + gmx::ArrayRef forceOutput); + + //! Compute forces and virial tensor + void compute(gmx::ArrayRef coordinateInput, + const Box& box, + gmx::ArrayRef forceOutput, + gmx::ArrayRef virialOutput); + + //! Compute forces, virial tensor and potential energies + void compute(gmx::ArrayRef coordinateInput, + const Box& box, + gmx::ArrayRef forceOutput, + gmx::ArrayRef virialOutput, + gmx::ArrayRef energyOutput); + +private: + //! Private implementation + class CpuImpl; + std::unique_ptr impl_; +}; + +std::unique_ptr setupGmxForceCalculatorCpu(const Topology& topology, + const NBKernelOptions& options); + +} // namespace nblib + +#endif // NBLIB_GMXCALCULATORCPU_H diff --git a/api/nblib/nblib.h b/api/nblib/nblib.h index 636c92bc47..55f3620c08 100644 --- a/api/nblib/nblib.h +++ b/api/nblib/nblib.h @@ -47,8 +47,7 @@ #include "nblib/basicdefinitions.h" #include "nblib/box.h" -#include "nblib/gmxcalculator.h" -#include "nblib/gmxsetup.h" +#include "nblib/gmxcalculatorcpu.h" #include "nblib/integrator.h" #include "nblib/interactions.h" #include "nblib/kerneloptions.h" @@ -56,6 +55,7 @@ #include "nblib/listed_forces/definitions.h" #include "nblib/listed_forces/calculator.h" #include "nblib/molecules.h" +#include "nblib/nbnxmsetuphelpers.h" #include "nblib/particlesequencer.h" #include "nblib/particletype.h" #include "nblib/simulationstate.h" diff --git a/api/nblib/nbnxmsetuphelpers.cpp b/api/nblib/nbnxmsetuphelpers.cpp index 027a0fdeb1..160612f8e8 100644 --- a/api/nblib/nbnxmsetuphelpers.cpp +++ b/api/nblib/nbnxmsetuphelpers.cpp @@ -167,7 +167,7 @@ std::vector createNonBondedParameters(const std::vector& par return nonbondedParameters; } -gmx::StepWorkload createStepWorkload([[maybe_unused]] const NBKernelOptions& options) +gmx::StepWorkload createStepWorkload() { gmx::StepWorkload stepWorkload; stepWorkload.computeForces = true; diff --git a/api/nblib/nbnxmsetuphelpers.h b/api/nblib/nbnxmsetuphelpers.h index 9c3ee71f46..2f22bd825b 100644 --- a/api/nblib/nbnxmsetuphelpers.h +++ b/api/nblib/nbnxmsetuphelpers.h @@ -99,7 +99,7 @@ std::vector createNonBondedParameters(const std::vector& par const NonBondedInteractionMap& nonBondedInteractionMap); //! Create a step work object -gmx::StepWorkload createStepWorkload(const NBKernelOptions& options); +gmx::StepWorkload createStepWorkload(); //! Computes the Ewald splitting coefficient for Coulomb real ewaldCoeff(real ewald_rtol, real pairlistCutoff); diff --git a/api/nblib/samples/argon-forces-integration.cpp b/api/nblib/samples/argon-forces-integration.cpp index d979c4944a..3daa462ba9 100644 --- a/api/nblib/samples/argon-forces-integration.cpp +++ b/api/nblib/samples/argon-forces-integration.cpp @@ -107,15 +107,16 @@ int main() // Some performance flags can be set a run time options.nbnxmSimd = nblib::SimdKernels::SimdNo; // The force calculator contains all the data needed to compute forces. - auto forceCalculator = nblib::setupGmxForceCalculator( - simState.topology(), simState.coordinates(), simState.box(), options); + auto forceCalculator = nblib::setupGmxForceCalculatorCpu(simState.topology(), options); + // build the pairlist + forceCalculator->updatePairlist(simState.coordinates(), simState.box()); // Integration requires masses, positions, and forces nblib::LeapFrog integrator(simState.topology(), simState.box()); // Print some diagnostic info printf("initial forces on particle 0: x %4f y %4f z %4f\n", forces[0][0], forces[0][1], forces[0][2]); // The forces are computed for the user gmx::ArrayRef userForces(simState.forces()); - forceCalculator->compute(simState.coordinates(), userForces); + forceCalculator->compute(simState.coordinates(), simState.box(), userForces); // Print some diagnostic info printf(" final forces on particle 0: x %4f y %4f z %4f\n", userForces[0][0], diff --git a/api/nblib/samples/methane-water-integration.cpp b/api/nblib/samples/methane-water-integration.cpp index 4dbf3b46d7..06bc92b668 100644 --- a/api/nblib/samples/methane-water-integration.cpp +++ b/api/nblib/samples/methane-water-integration.cpp @@ -174,8 +174,10 @@ int main() SimulationState simulationState(coordinates, velocities, forces, box, topology); // The non-bonded force calculator contains all the data needed to compute forces - auto forceCalculator = nblib::setupGmxForceCalculator( - simulationState.topology(), simulationState.coordinates(), simulationState.box(), options); + auto forceCalculator = setupGmxForceCalculatorCpu(simulationState.topology(), options); + + // build the pair list + forceCalculator->updatePairlist(simulationState.coordinates(), simulationState.box()); // The listed force calculator is also initialized with the required arguments ListedForceCalculator listedForceCalculator( @@ -183,7 +185,7 @@ int main() // Integrator is initialized with an array of inverse masses (constructed from topology) and // the bounding box - LeapFrog integrator(topology, box); + LeapFrog integrator(simulationState.topology(), simulationState.box()); // Print some diagnostic info printf("initial position of particle 0: x %4f y %4f z %4f\n", @@ -198,7 +200,8 @@ int main() { zeroCartesianArray(simulationState.forces()); - forceCalculator->compute(simulationState.coordinates(), simulationState.forces()); + forceCalculator->compute( + simulationState.coordinates(), simulationState.box(), simulationState.forces()); listedForceCalculator.compute(simulationState.coordinates(), simulationState.forces()); diff --git a/api/nblib/systemdescription.h b/api/nblib/systemdescription.h new file mode 100644 index 0000000000..0d2a3e0a6a --- /dev/null +++ b/api/nblib/systemdescription.h @@ -0,0 +1,130 @@ +/* + * This file is part of the GROMACS molecular simulation package. + * + * Copyright (c) 2021, 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. + */ +/*! \inpublicapi \file + * \brief + * Implements nblib simulation box + * + * \author Victor Holanda + * \author Joe Jordan + * \author Prashanth Kanduri + * \author Sebastian Keller + */ +#ifndef NBLIB_SYSTEMDESCRIPTION_H +#define NBLIB_SYSTEMDESCRIPTION_H + +#include "gromacs/ewald/ewald_utils.h" +#include "gromacs/gmxlib/nrnb.h" +#include "gromacs/gpu_utils/device_stream_manager.h" +#include "gromacs/mdlib/rf_util.h" +#include "gromacs/mdtypes/enerdata.h" +#include "gromacs/mdtypes/forcerec.h" +#include "gromacs/mdtypes/interaction_const.h" +#include "gromacs/mdtypes/simulation_workload.h" +#include "gromacs/nbnxm/gpu_data_mgmt.h" +#include "gromacs/nbnxm/nbnxm_gpu.h" +#include "gromacs/nbnxm/atomdata.h" +#include "gromacs/nbnxm/nbnxm.h" +#include "gromacs/nbnxm/pairlistset.h" +#include "gromacs/nbnxm/pairlistsets.h" +#include "gromacs/nbnxm/pairsearch.h" +#include "gromacs/timing/wallcycle.h" +#include "gromacs/utility/listoflists.h" +#include "gromacs/utility/range.h" + +namespace nblib +{ + +//! \brief client-side provided system description data +struct SystemDescription +{ + SystemDescription() = default; + + SystemDescription(gmx::ArrayRef particleTypeIdOfAllParticles, + gmx::ArrayRef nonBondedParams, + gmx::ArrayRef charges, + gmx::ArrayRef particleInteractionFlags) + { + std::array inputSizes{ particleTypeIdOfAllParticles.size(), + charges.size(), + particleInteractionFlags.size() }; + if (static_cast(std::count(begin(inputSizes), end(inputSizes), inputSizes[0])) + != inputSizes.size()) + { + throw InputException("input array size inconsistent"); + } + + int numParticleTypes = int(std::round(std::sqrt(nonBondedParams.size() / 2))); + if (2 * numParticleTypes * numParticleTypes != int(nonBondedParams.size())) + { + throw InputException("Wrong size of nonBondedParams"); + } + + numParticles_ = particleTypeIdOfAllParticles.size(); + numParticleTypes_ = numParticleTypes; + + particleTypeIdOfAllParticles_ = std::vector(particleTypeIdOfAllParticles.begin(), + particleTypeIdOfAllParticles.end()); + + nonBondedParams_ = std::vector(nonBondedParams.begin(), nonBondedParams.end()); + charges_ = std::vector(charges.begin(), charges.end()); + particleInfo_ = + std::vector(particleInteractionFlags.begin(), particleInteractionFlags.end()); + } + + //! number of particles + size_t numParticles_{ 0 }; + + //! number of particle types + size_t numParticleTypes_{ 0 }; + + //! particle type id of all particles + std::vector particleTypeIdOfAllParticles_; + + //! Storage for parameters for short range interactions. + std::vector nonBondedParams_; + + //! electrostatic charges + std::vector charges_; + + //! flag for each particle to set LJ and Q interactions + std::vector particleInfo_; + + //! Legacy matrix for box + Box box_{ 0 }; +}; + + +} // namespace nblib +#endif // NBLIB_SYSTEMDESCRIPTION_H diff --git a/api/nblib/tests/gmxcalculator.cpp b/api/nblib/tests/gmxcalculator.cpp index 141ce59023..e133bc96fc 100644 --- a/api/nblib/tests/gmxcalculator.cpp +++ b/api/nblib/tests/gmxcalculator.cpp @@ -43,9 +43,10 @@ */ #include -#include "nblib/gmxsetup.h" +#include "nblib/gmxcalculatorcpu.h" #include "nblib/kerneloptions.h" #include "nblib/simulationstate.h" +#include "nblib/tests/testhelpers.h" #include "nblib/tests/testsystems.h" #include "gromacs/utility/arrayref.h" @@ -61,23 +62,72 @@ TEST(NBlibTest, GmxForceCalculatorCanCompute) SimulationState simState = argonSystemBuilder.setupSimulationState(); NBKernelOptions options = NBKernelOptions(); options.nbnxmSimd = SimdKernels::SimdNo; - std::unique_ptr gmxForceCalculator = nblib::setupGmxForceCalculator( - simState.topology(), simState.coordinates(), simState.box(), options); - ; - EXPECT_NO_THROW(gmxForceCalculator->compute(simState.coordinates(), simState.forces())); + + std::unique_ptr gmxForceCalculator = + setupGmxForceCalculatorCpu(simState.topology(), options); + gmxForceCalculator->updatePairlist(simState.coordinates(), simState.box()); + + EXPECT_NO_THROW(gmxForceCalculator->compute(simState.coordinates(), simState.box(), simState.forces())); } -TEST(NBlibTest, CanSetupStepWorkload) +TEST(NBlibTest, ArgonVirialsAreCorrect) { - NBKernelOptions options; - EXPECT_NO_THROW(NbvSetupUtil{}.setupStepWorkload(options)); + ArgonSimulationStateBuilder argonSystemBuilder(fftypes::OPLSA); + SimulationState simState = argonSystemBuilder.setupSimulationState(); + NBKernelOptions options = NBKernelOptions(); + options.nbnxmSimd = SimdKernels::SimdNo; + std::unique_ptr gmxForceCalculator = + setupGmxForceCalculatorCpu(simState.topology(), options); + gmxForceCalculator->updatePairlist(simState.coordinates(), simState.box()); + + std::vector virialArray(9, 0.0); + + gmxForceCalculator->compute(simState.coordinates(), simState.box(), simState.forces(), virialArray); + + RefDataChecker virialsOutputTest(1e-7); + virialsOutputTest.testArrays(virialArray, "Virials"); } -TEST(NBlibTest, GmxForceCalculatorCanSetupInteractionConst) +TEST(NBlibTest, ArgonEnergiesAreCorrect) { - NBKernelOptions options; - EXPECT_NO_THROW(NbvSetupUtil{}.setupInteractionConst(options)); + ArgonSimulationStateBuilder argonSystemBuilder(fftypes::OPLSA); + SimulationState simState = argonSystemBuilder.setupSimulationState(); + NBKernelOptions options = NBKernelOptions(); + options.nbnxmSimd = SimdKernels::SimdNo; + std::unique_ptr gmxForceCalculator = + setupGmxForceCalculatorCpu(simState.topology(), options); + gmxForceCalculator->updatePairlist(simState.coordinates(), simState.box()); + + // number of energy kinds is 5: COULSR, LJSR, BHAMSR, COUL14, LJ14, + std::vector energies(5, 0.0); + + gmxForceCalculator->compute( + simState.coordinates(), simState.box(), simState.forces(), gmx::ArrayRef{}, energies); + + RefDataChecker energiesOutputTest(5e-5); + energiesOutputTest.testArrays(energies, "Argon energies"); +} + +TEST(NBlibTest, SpcMethanolEnergiesAreCorrect) +{ + SpcMethanolSimulationStateBuilder spcMethanolSystemBuilder; + SimulationState simState = spcMethanolSystemBuilder.setupSimulationState(); + NBKernelOptions options = NBKernelOptions(); + options.nbnxmSimd = SimdKernels::SimdNo; + std::unique_ptr gmxForceCalculator = + setupGmxForceCalculatorCpu(simState.topology(), options); + gmxForceCalculator->updatePairlist(simState.coordinates(), simState.box()); + + // number of energy kinds is 5: COULSR, LJSR, BHAMSR, COUL14, LJ14, + std::vector energies(5, 0.0); + + gmxForceCalculator->compute( + simState.coordinates(), simState.box(), simState.forces(), gmx::ArrayRef{}, energies); + + RefDataChecker energiesOutputTest(5e-5); + energiesOutputTest.testArrays(energies, "SPC-methanol energies"); } + } // namespace } // namespace test } // namespace nblib diff --git a/api/nblib/tests/nbkernelsystem.cpp b/api/nblib/tests/nbkernelsystem.cpp index 15163e1384..96c874a393 100644 --- a/api/nblib/tests/nbkernelsystem.cpp +++ b/api/nblib/tests/nbkernelsystem.cpp @@ -45,8 +45,7 @@ #include #include "gromacs/topology/exclusionblocks.h" -#include "nblib/gmxcalculator.h" -#include "nblib/gmxsetup.h" +#include "nblib/gmxcalculatorcpu.h" #include "nblib/integrator.h" #include "nblib/tests/testhelpers.h" #include "nblib/tests/testsystems.h" @@ -68,12 +67,14 @@ TEST(NBlibTest, SpcMethanolForcesAreCorrect) SpcMethanolSimulationStateBuilder spcMethanolSystemBuilder; - auto simState = spcMethanolSystemBuilder.setupSimulationState(); - auto forceCalculator = nblib::setupGmxForceCalculator( - simState.topology(), simState.coordinates(), simState.box(), options); + auto simState = spcMethanolSystemBuilder.setupSimulationState(); + + auto forceCalculator = setupGmxForceCalculatorCpu(simState.topology(), options); + + forceCalculator->updatePairlist(simState.coordinates(), simState.box()); gmx::ArrayRef forces(simState.forces()); - ASSERT_NO_THROW(forceCalculator->compute(simState.coordinates(), forces)); + ASSERT_NO_THROW(forceCalculator->compute(simState.coordinates(), simState.box(), forces)); RefDataChecker forcesOutputTest(5e-5); forcesOutputTest.testArrays(forces, "SPC-methanol forces"); @@ -86,12 +87,14 @@ TEST(NBlibTest, ExpectedNumberOfForces) SpcMethanolSimulationStateBuilder spcMethanolSystemBuilder; - auto simState = spcMethanolSystemBuilder.setupSimulationState(); - auto forceCalculator = nblib::setupGmxForceCalculator( - simState.topology(), simState.coordinates(), simState.box(), options); + auto simState = spcMethanolSystemBuilder.setupSimulationState(); + + auto forceCalculator = setupGmxForceCalculatorCpu(simState.topology(), options); + + forceCalculator->updatePairlist(simState.coordinates(), simState.box()); gmx::ArrayRef forces(simState.forces()); - forceCalculator->compute(simState.coordinates(), forces); + forceCalculator->compute(simState.coordinates(), simState.box(), forces); EXPECT_EQ(simState.topology().numParticles(), forces.size()); } @@ -103,18 +106,19 @@ TEST(NBlibTest, CanIntegrateSystem) SpcMethanolSimulationStateBuilder spcMethanolSystemBuilder; - auto simState = spcMethanolSystemBuilder.setupSimulationState(); - auto forceCalculator = nblib::setupGmxForceCalculator( - simState.topology(), simState.coordinates(), simState.box(), options); + auto simState = spcMethanolSystemBuilder.setupSimulationState(); + + auto forceCalculator = setupGmxForceCalculatorCpu(simState.topology(), options); + + forceCalculator->updatePairlist(simState.coordinates(), simState.box()); LeapFrog integrator(simState.topology(), simState.box()); for (int iter = 0; iter < options.numIterations; iter++) { gmx::ArrayRef forces(simState.forces()); - forceCalculator->compute(simState.coordinates(), simState.forces()); - EXPECT_NO_THROW(integrator.integrate( - 1.0, simState.coordinates(), simState.velocities(), simState.forces())); + forceCalculator->compute(simState.coordinates(), simState.box(), forces); + EXPECT_NO_THROW(integrator.integrate(1.0, simState.coordinates(), simState.velocities(), forces)); } } @@ -139,15 +143,18 @@ TEST(NBlibTest, UpdateChangesForces) SpcMethanolSimulationStateBuilder spcMethanolSystemBuilder; - auto simState = spcMethanolSystemBuilder.setupSimulationState(); - auto forceCalculator = nblib::setupGmxForceCalculator( - simState.topology(), simState.coordinates(), simState.box(), options); + auto simState = spcMethanolSystemBuilder.setupSimulationState(); + + auto forceCalculator = setupGmxForceCalculatorCpu(simState.topology(), options); + + forceCalculator->updatePairlist(simState.coordinates(), simState.box()); LeapFrog integrator(simState.topology(), simState.box()); // step 1 gmx::ArrayRef forces(simState.forces()); - forceCalculator->compute(simState.coordinates(), simState.forces()); + + forceCalculator->compute(simState.coordinates(), simState.box(), simState.forces()); // copy computed forces to another array std::vector forces_1(forces.size()); @@ -157,7 +164,7 @@ TEST(NBlibTest, UpdateChangesForces) zeroCartesianArray(forces); // check if forces change without update step - forceCalculator->compute(simState.coordinates(), forces); + forceCalculator->compute(simState.coordinates(), simState.box(), forces); // check if forces change without update for (size_t i = 0; i < forces_1.size(); i++) @@ -175,7 +182,8 @@ TEST(NBlibTest, UpdateChangesForces) zeroCartesianArray(forces); // step 2 - forceCalculator->compute(simState.coordinates(), forces); + forceCalculator->compute(simState.coordinates(), simState.box(), forces); + std::vector forces_2(forces.size()); std::copy(forces.begin(), forces.end(), begin(forces_2)); @@ -197,12 +205,14 @@ TEST(NBlibTest, ArgonOplsaForcesAreCorrect) ArgonSimulationStateBuilder argonSystemBuilder(fftypes::OPLSA); - auto simState = argonSystemBuilder.setupSimulationState(); - auto forceCalculator = nblib::setupGmxForceCalculator( - simState.topology(), simState.coordinates(), simState.box(), options); + auto simState = argonSystemBuilder.setupSimulationState(); + + auto forceCalculator = setupGmxForceCalculatorCpu(simState.topology(), options); + + forceCalculator->updatePairlist(simState.coordinates(), simState.box()); gmx::ArrayRef testForces(simState.forces()); - forceCalculator->compute(simState.coordinates(), simState.forces()); + forceCalculator->compute(simState.coordinates(), simState.box(), simState.forces()); RefDataChecker forcesOutputTest(1e-7); forcesOutputTest.testArrays(testForces, "Argon forces"); @@ -216,14 +226,16 @@ TEST(NBlibTest, ArgonGromos43A1ForcesAreCorrect) ArgonSimulationStateBuilder argonSystemBuilder(fftypes::GROMOS43A1); - auto simState = argonSystemBuilder.setupSimulationState(); - auto forceCalculator = nblib::setupGmxForceCalculator( - simState.topology(), simState.coordinates(), simState.box(), options); + auto simState = argonSystemBuilder.setupSimulationState(); + + auto forceCalculator = setupGmxForceCalculatorCpu(simState.topology(), options); + + forceCalculator->updatePairlist(simState.coordinates(), simState.box()); gmx::ArrayRef testForces(simState.forces()); - forceCalculator->compute(simState.coordinates(), simState.forces()); + forceCalculator->compute(simState.coordinates(), simState.box(), simState.forces()); - RefDataChecker forcesOutputTest(1e-7); + RefDataChecker forcesOutputTest; forcesOutputTest.testArrays(testForces, "Argon forces"); } diff --git a/api/nblib/tests/refdata/NBlibTest_ArgonEnergiesAreCorrect.xml b/api/nblib/tests/refdata/NBlibTest_ArgonEnergiesAreCorrect.xml new file mode 100644 index 0000000000..08bfa0d8b7 --- /dev/null +++ b/api/nblib/tests/refdata/NBlibTest_ArgonEnergiesAreCorrect.xml @@ -0,0 +1,12 @@ + + + + + 5 + 0 + -1.0175924 + 0 + 0 + 0 + + diff --git a/api/nblib/tests/refdata/NBlibTest_ArgonVirialsAreCorrect.xml b/api/nblib/tests/refdata/NBlibTest_ArgonVirialsAreCorrect.xml new file mode 100644 index 0000000000..d648ef1aa7 --- /dev/null +++ b/api/nblib/tests/refdata/NBlibTest_ArgonVirialsAreCorrect.xml @@ -0,0 +1,16 @@ + + + + + 9 + 0.094908442751139188 + 0.25238615516784413 + 0.02601194356883077 + 0.25238615516784413 + 0.67116022003893372 + 0.069172501786742557 + 0.026011943568830798 + 0.069172501786742391 + 0.0071291993484943611 + + \ No newline at end of file -- 2.22.0