--- /dev/null
-# Test goal: GCC with newest CUDA; Mdrun-only build
++# Test goal: GCC with newest CUDA
+ # Test intents (should change rarely and conservatively):
+ # OS: Ubuntu oldest supported
+ # GPU: CUDA newest supported
+ # HW: NVIDIA GPU
-# Features: Mdrun-only build
+ # Scope: configure, build, unit tests
+ # Test implementation choices (free to change as needed):
+ # OS: Ubuntu 18.04
+ # Build type: RelWithAssert
+ # Compiler: GCC 8
+ # MPI: thread_MPI
+ # GPU: CUDA 11.0
+ # SIMD: AVX2_256
+ # FFT: FFTW3
+ # Parallelism nt/ntomp: 4/2 (unit tests)
+
+ gromacs:gcc-8-cuda-11.0:release:configure:
+ extends:
+ - .gromacs:base:release:configure
+ - .use-gcc:base
+ - .use-mpi
+ - .use-cuda
+ - .rules:nightly-only-for-release
+ image: ${CI_REGISTRY}/gromacs/gromacs/ci-ubuntu-18.04-gcc-8-cuda-11.0
+ variables:
+ CMAKE: /usr/local/cmake-3.15.7/bin/cmake
+ COMPILER_MAJOR_VERSION: 8
+ RELEASE_BUILD_DIR: release-builds-gcc
- CMAKE_EXTRA_OPTIONS: "-DGMX_BUILD_MDRUN_ONLY=ON"
+ CMAKE_BUILD_TYPE_OPTIONS : "-DCMAKE_BUILD_TYPE=RelWithAssert"
+ CMAKE_REGRESSIONTEST_OPTIONS: ""
+ dependencies:
+ - archive:package
+ - regressiontests:package
+ - prepare-release-version
+
+ gromacs:gcc-8-cuda-11.0:release:build:
+ extends:
+ - .variables:default
+ - .gromacs:base:build
+ - .before_script:default
+ - .use-ccache
+ - .rules:nightly-only-for-release
+ stage: release-build
+ variables:
+ CMAKE: /usr/local/cmake-3.15.7/bin/cmake
+ BUILD_DIR: release-builds-gcc
+ image: ${CI_REGISTRY}/gromacs/gromacs/ci-ubuntu-18.04-gcc-8-cuda-11.0
+ needs:
+ - job: gromacs:gcc-8-cuda-11.0:release:configure
+
+ gromacs:gcc-8-cuda-11.0:release:test:
+ extends:
+ - .gromacs:base:test
+ - .rules:nightly-only-for-release
+ stage: release-tests
+ image: ${CI_REGISTRY}/gromacs/gromacs/ci-ubuntu-18.04-gcc-8-cuda-11.0
+ variables:
+ CMAKE: /usr/local/cmake-3.15.7/bin/cmake
+ KUBERNETES_EXTENDED_RESOURCE_NAME: "nvidia.com/gpu"
+ KUBERNETES_EXTENDED_RESOURCE_LIMIT: 1
+ BUILD_DIR: release-builds-gcc
+ tags:
+ - k8s-scilifelab
+ needs:
+ - job: gromacs:gcc-8-cuda-11.0:release:configure
+ - job: gromacs:gcc-8-cuda-11.0:release:build
+
auto atomData = std::make_unique<nbnxn_atomdata_t>(pinPolicy);
- // Put everything together
- auto nbv = std::make_unique<nonbonded_verlet_t>(
- std::move(pairlistSets), std::move(pairSearch), std::move(atomData), kernelSetup, nullptr, nullWallcycle);
-
// Needs to be called with the number of unique ParticleTypes
- nbnxn_atomdata_init(gmx::MDLogger(), atomData.get(), kernelSetup.kernelType, combinationRule,
- numParticleTypes, nonbondedParameters_, 1, numThreads);
+ nbnxn_atomdata_init(gmx::MDLogger(),
- nbv->nbat.get(),
++ atomData.get(),
+ kernelSetup.kernelType,
+ combinationRule,
+ numParticleTypes,
+ nonbondedParameters_,
+ 1,
+ numThreads);
- auto nbv = std::make_unique<nonbonded_verlet_t>(std::move(pairlistSets), std::move(pairSearch),
- std::move(atomData), kernelSetup, nullptr,
- nullWallcycle);
+ // Put everything together
++ auto nbv = std::make_unique<nonbonded_verlet_t>(
++ std::move(pairlistSets), std::move(pairSearch), std::move(atomData), kernelSetup, nullptr, nullWallcycle);
+
gmxForceCalculator_->nbv_ = std::move(nbv);
}
gmxForceCalculator_->setParticlesOnGrid(particleInfoAllVdw_, coordinates, box);
}
- void NbvSetupUtil::constructPairList(const gmx::ListOfLists<int>& exclusions)
+ void NbvSetupUtil::constructPairList(ExclusionLists<int> exclusionLists)
{
- gmxForceCalculator_->nbv_->constructPairlist(gmx::InteractionLocality::Local, exclusions, 0,
- gmxForceCalculator_->nrnb_.get());
+ gmx::ListOfLists<int> exclusions(std::move(exclusionLists.ListRanges),
+ std::move(exclusionLists.ListElements));
+ gmxForceCalculator_->nbv_->constructPairlist(
+ gmx::InteractionLocality::Local, exclusions, 0, gmxForceCalculator_->nrnb_.get());
}
void sortInteractions(ListedInteractionData& interactions)
{
- auto sortKeyObj = [](const auto& lhs, const auto& rhs) { return interactionSortKey(lhs, rhs); };
- auto sortOneElement = [comparison = sortKeyObj](auto& interactionElement) {
- std::sort(begin(interactionElement.indices), end(interactionElement.indices), comparison);
+ auto sortOneElement = [](auto& interactionElement) {
+ using InteractionContainerType = std::decay_t<decltype(interactionElement)>;
+ using InteractionType = typename InteractionContainerType::type;
+
- std::sort(begin(interactionElement.indices), end(interactionElement.indices),
++ std::sort(begin(interactionElement.indices),
++ end(interactionElement.indices),
+ interactionSortKey<InteractionType>);
};
for_each_tuple(sortOneElement, interactions);
return *this;
}
- //! Two-particle interactions such as bonds and LJ1-4
- template<class Interaction>
- void Molecule::addInteraction(const ParticleName& particleNameI,
- const ResidueName& residueNameI,
- const ParticleName& particleNameJ,
- const ResidueName& residueNameJ,
- const Interaction& interaction)
+ ResidueName Molecule::residueName(const ParticleIdentifier& particleIdentifier)
+ {
+ return (particleIdentifier.residueName() == ResidueName{}) ? ResidueName(name_)
+ : particleIdentifier.residueName();
+ }
+
+ template<class ListedVariant, class... ParticleIdentifiers>
+ void Molecule::addInteractionImpl(const ListedVariant& interaction, const ParticleIdentifiers&... particles)
{
- if (particleNameI == particleNameJ and residueNameI == residueNameJ)
+ auto storeInteraction = [&](const auto& interaction_) {
+ using Interaction = std::decay_t<decltype(interaction_)>;
+
+ auto& interactionContainer = pickType<Interaction>(interactionData_);
+ interactionContainer.interactions_.emplace_back(particles...);
+ interactionContainer.interactionTypes_.push_back(interaction_);
+ };
+
+ // add the interaction to the correct location in interactionData_
+ std::visit(storeInteraction, interaction);
+ }
+
+ void Molecule::addInteraction(const ParticleIdentifier& particleI,
+ const ParticleIdentifier& particleJ,
+ const TwoCenterInteraction& interaction)
+ {
+ if (particleI == particleJ)
{
throw InputException(std::string("Cannot add interaction of particle ")
- + particleNameI.value() + " with itself in molecule " + name_.value());
+ + particleI.particleName().value() + " with itself in molecule "
+ + name_.value());
}
- auto& interactionContainer = pickType<Interaction>(interactionData_);
- interactionContainer.interactions_.emplace_back(particleNameI, residueNameI, particleNameJ, residueNameJ);
- interactionContainer.interactionTypes_.push_back(interaction);
- addInteractionImpl(interaction, particleI.particleName(), residueName(particleI),
- particleJ.particleName(), residueName(particleJ));
++ addInteractionImpl(interaction,
++ particleI.particleName(),
++ residueName(particleI),
++ particleJ.particleName(),
++ residueName(particleJ));
}
- //! \cond DO_NOT_DOCUMENT
- #define ADD_INTERACTION_INSTANTIATE_TEMPLATE(x) \
- template void Molecule::addInteraction(const ParticleName& particleNameI, \
- const ResidueName& residueNameI, \
- const ParticleName& particleNameJ, \
- const ResidueName& residueNameJ, \
- const x& interaction);
- MAP(ADD_INTERACTION_INSTANTIATE_TEMPLATE, SUPPORTED_TWO_CENTER_TYPES)
- #undef ADD_INTERACTION_INSTANTIATE_TEMPLATE
- //! \endcond
-
- // add interactions with default residue name
- template<class Interaction>
- void Molecule::addInteraction(const ParticleName& particleNameI,
- const ParticleName& particleNameJ,
- const Interaction& interaction)
+ void Molecule::addInteraction(const ParticleIdentifier& particleI,
+ const ParticleIdentifier& particleJ,
+ const ParticleIdentifier& particleK,
+ const ThreeCenterInteraction& interaction)
{
- addInteraction(particleNameI, ResidueName(name_), particleNameJ, ResidueName(name_), interaction);
+ if (particleI == particleJ and particleJ == particleK)
+ {
+ throw InputException(std::string("Cannot add interaction of particle ")
+ + particleI.particleName().value() + " with itself in molecule "
+ + name_.value());
+ }
+
- addInteractionImpl(interaction, particleI.particleName(), residueName(particleI),
- particleJ.particleName(), residueName(particleJ), particleK.particleName(),
++ addInteractionImpl(interaction,
++ particleI.particleName(),
++ residueName(particleI),
++ particleJ.particleName(),
++ residueName(particleJ),
++ particleK.particleName(),
+ residueName(particleK));
}
- //! \cond DO_NOT_DOCUMENT
- #define ADD_INTERACTION_INSTANTIATE_TEMPLATE(x) \
- template void Molecule::addInteraction( \
- const ParticleName& particleNameI, const ParticleName& particleNameJ, const x& interaction);
- MAP(ADD_INTERACTION_INSTANTIATE_TEMPLATE, SUPPORTED_TWO_CENTER_TYPES)
- #undef ADD_INTERACTION_INSTANTIATE_TEMPLATE
-
- //! 3-particle interactions such as angles
- template<class Interaction>
- void Molecule::addInteraction(const ParticleName& particleNameI,
- const ResidueName& residueNameI,
- const ParticleName& particleNameJ,
- const ResidueName& residueNameJ,
- const ParticleName& particleNameK,
- const ResidueName& residueNameK,
- const Interaction& interaction)
+ void Molecule::addInteraction(const ParticleIdentifier& particleI,
+ const ParticleIdentifier& particleJ,
+ const ParticleIdentifier& particleK,
+ const ParticleIdentifier& particleL,
+ const FourCenterInteraction& interaction)
{
- if (particleNameI == particleNameJ and particleNameJ == particleNameK)
+ if (particleI == particleJ and particleJ == particleK and particleK == particleL)
{
throw InputException(std::string("Cannot add interaction of particle ")
- + particleNameI.value() + " with itself in molecule " + name_.value());
+ + particleI.particleName().value() + " with itself in molecule "
+ + name_.value());
}
- auto& interactionContainer = pickType<Interaction>(interactionData_);
- interactionContainer.interactions_.emplace_back(
- particleNameI, residueNameI, particleNameJ, residueNameJ, particleNameK, residueNameK);
- interactionContainer.interactionTypes_.push_back(interaction);
- addInteractionImpl(interaction, particleI.particleName(), residueName(particleI),
- particleJ.particleName(), residueName(particleJ), particleK.particleName(),
- residueName(particleK), particleL.particleName(), residueName(particleL));
++ addInteractionImpl(interaction,
++ particleI.particleName(),
++ residueName(particleI),
++ particleJ.particleName(),
++ residueName(particleJ),
++ particleK.particleName(),
++ residueName(particleK),
++ particleL.particleName(),
++ residueName(particleL));
}
- #define ADD_INTERACTION_INSTANTIATE_TEMPLATE(x) \
- template void Molecule::addInteraction(const ParticleName& particleNameI, \
- const ResidueName& residueNameI, \
- const ParticleName& particleNameJ, \
- const ResidueName& residueNameJ, \
- const ParticleName& particleNameK, \
- const ResidueName& residueNameK, \
- const x& interaction);
- MAP(ADD_INTERACTION_INSTANTIATE_TEMPLATE, SUPPORTED_THREE_CENTER_TYPES)
- #undef ADD_INTERACTION_INSTANTIATE_TEMPLATE
-
- template<class Interaction>
- void Molecule::addInteraction(const ParticleName& particleNameI,
- const ParticleName& particleNameJ,
- const ParticleName& particleNameK,
- const Interaction& interaction)
+ void Molecule::addInteraction(const ParticleIdentifier& particleI,
+ const ParticleIdentifier& particleJ,
+ const ParticleIdentifier& particleK,
+ const ParticleIdentifier& particleL,
+ const ParticleIdentifier& particleM,
+ const FiveCenterInteraction& interaction)
{
- addInteraction(particleNameI,
- ResidueName(name_),
- particleNameJ,
- ResidueName(name_),
- particleNameK,
- ResidueName(name_),
- interaction);
+ if (particleI == particleJ and particleJ == particleK and particleK == particleL and particleL == particleM)
+ {
+ throw InputException(std::string("Cannot add interaction of particle ")
+ + particleI.particleName().value() + " with itself in molecule "
+ + name_.value());
+ }
+
- addInteractionImpl(interaction, particleI.particleName(), residueName(particleI),
- particleJ.particleName(), residueName(particleJ), particleK.particleName(),
- residueName(particleK), particleL.particleName(), residueName(particleL),
- particleM.particleName(), residueName(particleM));
++ addInteractionImpl(interaction,
++ particleI.particleName(),
++ residueName(particleI),
++ particleJ.particleName(),
++ residueName(particleJ),
++ particleK.particleName(),
++ residueName(particleK),
++ particleL.particleName(),
++ residueName(particleL),
++ particleM.particleName(),
++ residueName(particleM));
}
- #define ADD_INTERACTION_INSTANTIATE_TEMPLATE(x) \
- template void Molecule::addInteraction(const ParticleName& particleNameI, \
- const ParticleName& particleNameJ, \
- const ParticleName& particleNameK, \
- const x& interaction);
- MAP(ADD_INTERACTION_INSTANTIATE_TEMPLATE, SUPPORTED_THREE_CENTER_TYPES)
- #undef ADD_INTERACTION_INSTANTIATE_TEMPLATE
- //! \endcond
int Molecule::numParticlesInMolecule() const
{
return particles_.size();
}
- void Molecule::addExclusion(const int particleIndex, const int particleIndexToExclude)
+ void Molecule::addExclusion(const ParticleIdentifier& particle, const ParticleIdentifier& particleToExclude)
{
- // We do not need to add exclusion in case the particle indexes are the same
- // because self exclusion are added by addParticle
- if (particleIndex != particleIndexToExclude)
+ if (particle == particleToExclude)
{
- exclusions_.emplace_back(particleIndex, particleIndexToExclude);
- exclusions_.emplace_back(particleIndexToExclude, particleIndex);
+ return;
}
- }
- void Molecule::addExclusion(std::tuple<ParticleName, ResidueName> particle,
- std::tuple<ParticleName, ResidueName> particleToExclude)
- {
// duplication for the swapped pair happens in getExclusions()
- exclusionsByName_.emplace_back(std::make_tuple(std::get<0>(particle),
- std::get<1>(particle),
- std::get<0>(particleToExclude),
- std::get<1>(particleToExclude)));
- }
-
- void Molecule::addExclusion(const ParticleName& particleName, const ParticleName& particleNameToExclude)
- {
- addExclusion(std::make_tuple(particleName, ResidueName(name_)),
- std::make_tuple(particleNameToExclude, ResidueName(name_)));
- exclusionsByName_.emplace_back(std::make_tuple(particle.particleName(), residueName(particle),
++ exclusionsByName_.emplace_back(std::make_tuple(particle.particleName(),
++ residueName(particle),
+ particleToExclude.particleName(),
+ residueName(particleToExclude)));
}
const Molecule::InteractionTuple& Molecule::interactionData() const
--- /dev/null
- printf("No particle %s in residue %s in molecule %s found\n", particleName.value().c_str(),
- residueName.value().c_str(), moleculeName.value().c_str());
+ /*
+ * 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 ParticleSequencer class
+ *
+ * \author Victor Holanda <victor.holanda@cscs.ch>
+ * \author Joe Jordan <ejjordan@kth.se>
+ * \author Prashanth Kanduri <kanduri@cscs.ch>
+ * \author Sebastian Keller <keller@cscs.ch>
+ * \author Artem Zhmurov <zhmurov@gmail.com>
+ */
+
+ #include <algorithm>
+
+ #include "nblib/exception.h"
+ #include "nblib/particlesequencer.h"
+
+ namespace nblib
+ {
+
+ int ParticleSequencer::operator()(const MoleculeName& moleculeName,
+ int moleculeNr,
+ const ResidueName& residueName,
+ const ParticleName& particleName) const
+ {
+ try
+ {
+ return data_.at(moleculeName).at(moleculeNr).at(residueName).at(particleName);
+ }
+ catch (const std::out_of_range& outOfRange)
+ {
+ // TODO: use string format function once we have it
+ if (moleculeName.value() == residueName.value())
+ {
- printf("No particle %s in molecule %s found\n", particleName.value().c_str(),
++ printf("No particle %s in residue %s in molecule %s found\n",
++ particleName.value().c_str(),
++ residueName.value().c_str(),
++ moleculeName.value().c_str());
+ }
+ else
+ {
++ printf("No particle %s in molecule %s found\n",
++ particleName.value().c_str(),
+ moleculeName.value().c_str());
+ }
+
+ throw InputException(outOfRange.what());
+ }
+ }
+
+ void ParticleSequencer::build(const std::vector<std::tuple<Molecule, int>>& moleculesList)
+ {
+ int currentID = 0;
+ for (const auto& molNumberTuple : moleculesList)
+ {
+ const Molecule& molecule = std::get<0>(molNumberTuple);
+ const size_t numMols = std::get<1>(molNumberTuple);
+
+ auto& moleculeMap = data_[molecule.name()];
+
+ for (size_t i = 0; i < numMols; ++i)
+ {
+ auto& moleculeNrMap = moleculeMap[i];
+ for (int j = 0; j < molecule.numParticlesInMolecule(); ++j)
+ {
+ moleculeNrMap[molecule.residueName(j)][molecule.particleName(j)] = currentID++;
+ }
+ }
+ }
+ }
+
+ } // namespace nblib
--- /dev/null
- ListedForceCalculator listedForceCalculator(topology.getInteractionData(),
- topology.numParticles(), 4, box);
+ /*
+ * 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.
+ */
+ /*! \internal \file
+ * \brief
+ * This tests that sample code can run
+ *
+ * \author Victor Holanda <victor.holanda@cscs.ch>
+ * \author Joe Jordan <ejjordan@kth.se>
+ * \author Prashanth Kanduri <kanduri@cscs.ch>
+ * \author Sebastian Keller <keller@cscs.ch>
+ */
+
+ #include <cstdio>
+
+ #include "gromacs/utility/arrayref.h"
+
+ // The entire nblib public API can be included with a single header or individual components
+ // can be included via their respective headers.
+ #include "nblib/nblib.h"
+
+ using namespace nblib;
+
+ // Main function to write the MD program.
+ int main(); // Keep the compiler happy
+
+ int main()
+ {
+ // Create the particles
+ ParticleType Ow(ParticleTypeName("Ow"), Mass(15.999));
+ ParticleType Hw(ParticleTypeName("Hw"), Mass(1.00784));
+ ParticleType Cm(ParticleTypeName("Cm"), Mass(12.0107));
+ ParticleType Hc(ParticleTypeName("Hc"), Mass(1.00784));
+
+ ParticleTypesInteractions interactions(CombinationRule::Geometric);
+
+ // Parameters from a GROMOS compatible force-field 2016H66
+ // add non-bonded interactions for the particle types
+ interactions.add(Ow.name(), C6(0.0026173456), C12(2.634129e-06));
+ interactions.add(Hw.name(), C6(0.0), C12(0.0));
+ interactions.add(Cm.name(), C6(0.01317904), C12(34.363044e-06));
+ interactions.add(Hc.name(), C6(8.464e-05), C12(15.129e-09));
+
+ Molecule water(MoleculeName("Water"));
+ Molecule methane(MoleculeName("Methane"));
+
+ water.addParticle(ParticleName("O"), Ow);
+ water.addParticle(ParticleName("H1"), Hw);
+ water.addParticle(ParticleName("H2"), Hw);
+
+ water.addExclusion(ParticleName("H1"), ParticleName("O"));
+ water.addExclusion(ParticleName("H2"), ParticleName("O"));
+
+ methane.addParticle(ParticleName("C"), Cm);
+ methane.addParticle(ParticleName("H1"), Hc);
+ methane.addParticle(ParticleName("H2"), Hc);
+ methane.addParticle(ParticleName("H3"), Hc);
+ methane.addParticle(ParticleName("H4"), Hc);
+
+ methane.addExclusion(ParticleName("H1"), ParticleName("C"));
+ methane.addExclusion(ParticleName("H2"), ParticleName("C"));
+ methane.addExclusion(ParticleName("H3"), ParticleName("C"));
+ methane.addExclusion(ParticleName("H4"), ParticleName("C"));
+
+ HarmonicBondType ohHarmonicBond(1, 1);
+ HarmonicBondType hcHarmonicBond(2, 1);
+
+ HarmonicAngleType hohAngle(Degrees(120), 1);
+ HarmonicAngleType hchAngle(Degrees(109.5), 1);
+
+ // add harmonic bonds for water
+ water.addInteraction(ParticleName("O"), ParticleName("H1"), ohHarmonicBond);
+ water.addInteraction(ParticleName("O"), ParticleName("H2"), ohHarmonicBond);
+
+ // add the angle for water
+ water.addInteraction(ParticleName("H1"), ParticleName("O"), ParticleName("H2"), hohAngle);
+
+ // add harmonic bonds for methane
+ methane.addInteraction(ParticleName("H1"), ParticleName("C"), hcHarmonicBond);
+ methane.addInteraction(ParticleName("H2"), ParticleName("C"), hcHarmonicBond);
+ methane.addInteraction(ParticleName("H3"), ParticleName("C"), hcHarmonicBond);
+ methane.addInteraction(ParticleName("H4"), ParticleName("C"), hcHarmonicBond);
+
+ // add the angles for methane
+ methane.addInteraction(ParticleName("H1"), ParticleName("C"), ParticleName("H2"), hchAngle);
+ methane.addInteraction(ParticleName("H1"), ParticleName("C"), ParticleName("H3"), hchAngle);
+ methane.addInteraction(ParticleName("H1"), ParticleName("C"), ParticleName("H4"), hchAngle);
+ methane.addInteraction(ParticleName("H2"), ParticleName("C"), ParticleName("H3"), hchAngle);
+ methane.addInteraction(ParticleName("H2"), ParticleName("C"), ParticleName("H4"), hchAngle);
+ methane.addInteraction(ParticleName("H3"), ParticleName("C"), ParticleName("H4"), hchAngle);
+
+ // Define a box for the simulation
+ Box box(6.05449);
+
+ // Define options for the non-bonded kernels
+ NBKernelOptions options;
+ // Use a simple cutoff rule for Coulomb
+ options.coulombType = nblib::CoulombType::Cutoff;
+ // Disable SIMD for this example
+ options.nbnxmSimd = SimdKernels::SimdNo;
+
+ TopologyBuilder topologyBuilder;
+
+ // add molecules
+ topologyBuilder.addMolecule(water, 2);
+ topologyBuilder.addMolecule(methane, 1);
+
+ // add non-bonded interaction map
+ topologyBuilder.addParticleTypesInteractions(interactions);
+
+ Topology topology = topologyBuilder.buildTopology();
+
+ // User defined coordinates.
+ std::vector<Vec3> coordinates = {
+ { 0.005, 0.600, 0.244 }, // Oxygen from water_1
+ { -0.017, 0.690, 0.270 }, // Hydrogen_1 from water_1
+ { 0.051, 0.610, 0.161 }, // Hydrogen_2 from water_1
+ { 0.155, 0.341, 0.735 }, // Oxygen from water_2
+ { 0.140, 0.284, 0.660 }, // Hydrogen_1 from water_2
+ { 0.081, 0.402, 0.734 }, // Hydrogen_2 from water_2
+ { -0.024, -0.222, -0.640 }, // Carbon from methane_1
+ { -0.083, -0.303, -0.646 }, // Hydrogen_1 from methane_1
+ { -0.080, -0.140, -0.642 }, // Hydrogen_2 from methane_1
+ { 0.040, -0.221, -0.716 }, // Hydrogen_3 from methane_1
+ { 0.027, -0.225, -0.553 } // Hydrogen_4 from methane_1
+ };
+
+ // User defined velocities.
+ std::vector<Vec3> velocities = {
+ { 0.1823, -0.4158, 0.487 }, // Oxygen from water_1
+ { -1.7457, -0.5883, -0.460 }, // Hydrogen_1 from water_1
+ { 2.5085, -0.1501, 1.762 }, // Hydrogen_2 from water_1
+ { 0.6282, 0.4390, 0.001 }, // Oxygen from water_2
+ { -0.3206, 0.0700, 0.4630 }, // Hydrogen_1 from water_2
+ { -0.1556, -0.4529, 1.440 }, // Hydrogen_2 from water_2
+ { 0.0, 0.0, 0.0 }, // Carbon from methane_1
+ { 0.0, 0.0, 0.0 }, // Hydrogen_1 from methane_1
+ { 0.0, 0.0, 0.0 }, // Hydrogen_2 from methane_1
+ { 0.0, 0.0, 0.0 }, // Hydrogen_3 from methane_1
+ { 0.0, 0.0, 0.0 }, // Hydrogen_4 from methane_1
+ };
+
+ // Force buffer initialization for each particle.
+ std::vector<Vec3> forces(topology.numParticles(), Vec3{ 0, 0, 0 });
+
+ SimulationState simulationState(coordinates, velocities, forces, box, topology);
+
+ // The non-bonded force calculator contains all the data needed to compute forces
+ ForceCalculator forceCalculator(simulationState, options);
+
+ // The listed force calculator is also initialized with the required arguments
- printf("initial position of particle 0: x %4f y %4f z %4f\n", simulationState.coordinates()[0][0],
- simulationState.coordinates()[0][1], simulationState.coordinates()[0][2]);
++ ListedForceCalculator listedForceCalculator(
++ topology.getInteractionData(), topology.numParticles(), 4, box);
+
+ // Integrator is initialized with an array of inverse masses (constructed from topology) and
+ // the bounding box
+ LeapFrog integrator(topology, box);
+
+ // Print some diagnostic info
- integrator.integrate(1.0, simulationState.coordinates(), simulationState.velocities(),
- simulationState.forces());
++ printf("initial position of particle 0: x %4f y %4f z %4f\n",
++ simulationState.coordinates()[0][0],
++ simulationState.coordinates()[0][1],
++ simulationState.coordinates()[0][2]);
+
+ // MD Loop
+ int numSteps = 2;
+
+ for (auto i = 0; i < numSteps; i++)
+ {
+ zeroCartesianArray(simulationState.forces());
+
+ forceCalculator.compute(simulationState.coordinates(), simulationState.forces());
+
+ listedForceCalculator.compute(simulationState.coordinates(), simulationState.forces());
+
+ // Integrate with a time step of 1 fs, positions, velocities and forces
- printf(" final position of particle 9: x %4f y %4f z %4f\n", simulationState.coordinates()[9][0],
- simulationState.coordinates()[9][1], simulationState.coordinates()[9][2]);
++ integrator.integrate(
++ 1.0, simulationState.coordinates(), simulationState.velocities(), simulationState.forces());
+ }
+
++ printf(" final position of particle 9: x %4f y %4f z %4f\n",
++ simulationState.coordinates()[9][0],
++ simulationState.coordinates()[9][1],
++ simulationState.coordinates()[9][2]);
+
+ return 0;
+ } // main
<< formatString(
"Velocity component {} of atom {} is different from analytical "
"solution at step {}.",
- d, i, step);
+ d,
+ i,
+ step);
}
- integrator.integrate(dt, simulationState.coordinates(), simulationState.velocities(),
++ integrator.integrate(dt,
++ simulationState.coordinates(),
++ simulationState.velocities(),
+ simulationState.forces());
}
- integrator.integrate(
- dt, simulationState.coordinates(), simulationState.velocities(), simulationState.forces());
}
}
// duplicate the exclusionBlockPerMolecule for the number of Molecules of (numMols)
for (size_t i = 0; i < numMols; ++i)
{
- auto offsetExclusions =
- detail::offsetGmxBlock(exclusionBlockPerMolecule, particleNumberOffset);
+ auto offsetExclusions = offsetGmxBlock(exclusionBlockPerMolecule, particleNumberOffset);
- std::copy(std::begin(offsetExclusions), std::end(offsetExclusions),
+ std::copy(std::begin(offsetExclusions),
+ std::end(offsetExclusions),
std::back_inserter(exclusionBlockGlobal));
particleNumberOffset += molecule.numParticlesInMolecule();
how-to/visualize.rst
install-guide/index.rst
release-notes/index.rst
+ release-notes/2022/major/highlights.rst
+ release-notes/2022/major/features.rst
+ release-notes/2022/major/performance.rst
+ release-notes/2022/major/tools.rst
+ release-notes/2022/major/bugs-fixed.rst
+ release-notes/2022/major/removed-functionality.rst
+ release-notes/2022/major/deprecated-functionality.rst
+ release-notes/2022/major/portability.rst
+ release-notes/2022/major/miscellaneous.rst
+ release-notes/2021/2021.1.rst
release-notes/2021/major/highlights.rst
release-notes/2021/major/features.rst
release-notes/2021/major/performance.rst
appropriate value instead of ``xxx`` :
* ``-DCMAKE_C_COMPILER=xxx`` equal to the name of the C99 `Compiler`_ you wish to use (or the environment variable ``CC``)
- * ``-DCMAKE_CXX_COMPILER=xxx`` equal to the name of the C++98 `compiler`_ you wish to use (or the environment variable ``CXX``)
+ * ``-DCMAKE_CXX_COMPILER=xxx`` equal to the name of the C++17 `compiler`_ you wish to use (or the environment variable ``CXX``)
-* ``-DGMX_MPI=on`` to build using `MPI support`_ (generally good to combine with `building only mdrun`_)
+* ``-DGMX_MPI=on`` to build using `MPI support`_
* ``-DGMX_GPU=CUDA`` to build with NVIDIA CUDA support enabled.
* ``-DGMX_GPU=OpenCL`` to build with OpenCL_ support enabled.
* ``-DGMX_SIMD=xxx`` to specify the level of `SIMD support`_ of the node on which |Gromacs| will run
.. todolist::
+ Patch releases
+ ^^^^^^^^^^^^^^
+
+ .. toctree::
+ :maxdepth: 1
+
+ 2021/2021.1
+
+
+Major release
+^^^^^^^^^^^^^
+
+.. toctree::
+ :maxdepth: 1
+
+ 2022/major/highlights
+ 2022/major/features
+ 2022/major/performance
+ 2022/major/tools
+ 2022/major/bugs-fixed
+ 2022/major/deprecated-functionality
+ 2022/major/removed-functionality
+ 2022/major/portability
+ 2022/major/miscellaneous
+
+
+|Gromacs| 2021 series
+---------------------
+
Major release
^^^^^^^^^^^^^
# TODO: Version management policy and procedures.
_major = 0
-_minor = 2
+_minor = 3
_micro = 0
- _suffix = 'b1'
+ _suffix = ''
# Reference https://www.python.org/dev/peps/pep-0440/
# and https://packaging.pypa.io/en/latest/version/
}
}
- void rename_bb(t_atoms* pdba, const char* oldnm, const char* newnm, bool bFullCompare, t_symtab* symtab)
+ /*! \brief Rename all residues named \c oldnm to \c newnm
+ *
+ * Search for residues for which the residue name from the input
+ * configuration file matches \c oldnm, and when found choose the rtp
+ * entry and name of \c newnm.
+ *
+ * \todo Refactor this function to accept a lambda that accepts i and
+ * numMatchesFound but always produces \c newnm. Then remove
+ * renameResiduesInteractively by calling this method with suitable
+ * lambdas that capture its parameter \c rr and ignores
+ * numMatchesFound. */
+ void renameResidue(const gmx::MDLogger& logger,
+ t_atoms* pdba,
+ const char* oldnm,
+ const char* newnm,
+ bool bFullCompare,
+ t_symtab* symtab)
{
- char* bbnm;
- int i;
-
- for (i = 0; (i < pdba->nres); i++)
+ int numMatchesFound = 0;
+ for (int i = 0; (i < pdba->nres); i++)
{
- /* We have not set the rtp name yes, use the residue name */
- bbnm = *pdba->resinfo[i].name;
- if ((bFullCompare && (gmx::equalCaseInsensitive(bbnm, oldnm)))
- || (!bFullCompare && strstr(bbnm, oldnm) != nullptr))
+ /* We have not set the rtp name yet, use the residue name */
+ const char* residueNameInInputConfiguration = *pdba->resinfo[i].name;
+ if ((bFullCompare && (gmx::equalCaseInsensitive(residueNameInInputConfiguration, oldnm)))
+ || (!bFullCompare && strstr(residueNameInInputConfiguration, oldnm) != nullptr))
{
- /* Change the rtp builing block name */
- pdba->resinfo[i].rtp = put_symtab(symtab, newnm);
+ /* Change the rtp building block name */
+ pdba->resinfo[i].rtp = put_symtab(symtab, newnm);
+ pdba->resinfo[i].name = pdba->resinfo[i].rtp;
+ numMatchesFound++;
}
}
- numMatchesFound, numMatchesFound > 1 ? "s" : "", oldnm, newnm);
+ if (numMatchesFound > 0)
+ {
+ GMX_LOG(logger.info)
+ .asParagraph()
+ .appendTextFormatted(
+ "Replaced %d residue%s named %s to the default %s. Use interactive "
+ "selection of protonated residues if that is what you need.",
++ numMatchesFound,
++ numMatchesFound > 1 ? "s" : "",
++ oldnm,
++ newnm);
+ }
}
- void rename_bbint(t_atoms* pdba,
- const char* oldnm,
- const char* gettp(int, gmx::ArrayRef<const RtpRename>),
- bool bFullCompare,
- t_symtab* symtab,
- gmx::ArrayRef<const RtpRename> rr)
+ /*! \brief Rename all residues named \c oldnm according to the user's
+ * interactive choice
+ *
+ * Search for residues for which the residue name from the input
+ * configuration file matches \c oldnm, and when found choose the rtp
+ * entry and name of the interactive choice from \c gettp.
+ *
+ * \todo Remove this function, per todo in \c renameResidue. */
+ void renameResidueInteractively(t_atoms* pdba,
+ const char* oldnm,
+ const char* gettp(int, gmx::ArrayRef<const RtpRename>),
+ bool bFullCompare,
+ t_symtab* symtab,
+ gmx::ArrayRef<const RtpRename> rr)
{
- int i;
- const char* ptr;
- char* bbnm;
-
- for (i = 0; i < pdba->nres; i++)
+ // Search for residues i for which the residue name from the input
+ // configuration file matches oldnm, so it can replaced by the rtp
+ // entry and name of newnm.
+ for (int i = 0; i < pdba->nres; i++)
{
/* We have not set the rtp name yet, use the residue name */
- bbnm = *pdba->resinfo[i].name;
- if ((bFullCompare && (strcmp(bbnm, oldnm) == 0)) || (!bFullCompare && strstr(bbnm, oldnm) != nullptr))
+ char* residueNameInInputConfiguration = *pdba->resinfo[i].name;
+ if ((bFullCompare && (strcmp(residueNameInInputConfiguration, oldnm) == 0))
+ || (!bFullCompare && strstr(residueNameInInputConfiguration, oldnm) != nullptr))
{
- ptr = gettp(i, rr);
- pdba->resinfo[i].rtp = put_symtab(symtab, ptr);
+ const char* interactiveRtpChoice = gettp(i, rr);
+ pdba->resinfo[i].rtp = put_symtab(symtab, interactiveRtpChoice);
+ pdba->resinfo[i].name = pdba->resinfo[i].rtp;
}
}
}
{
GMX_LOG(logger.info)
.asParagraph()
- .appendTextFormatted("Processing chain %d (%d atoms, %d residues)", chain + 1,
- natom, nres);
- }
-
- process_chain(logger, pdba, x, bUnA_, bUnA_, bUnA_, bLysMan_, bAspMan_, bGluMan_, bHisMan_,
- bArgMan_, bGlnMan_, angle_, distance_, &symtab, rtprename);
+ .appendTextFormatted(
+ "Processing chain %d (%d atoms, %d residues)", chain + 1, natom, nres);
+ }
+
- process_chain(pdba,
++ process_chain(logger,
++ pdba,
+ x,
+ bUnA_,
+ bUnA_,
+ bUnA_,
+ bLysMan_,
+ bAspMan_,
+ bGluMan_,
+ bHisMan_,
+ bArgMan_,
+ bGlnMan_,
+ angle_,
+ distance_,
+ &symtab,
+ rtprename);
cc->chainstart[cc->nterpairs] = pdba->nres;
j = 0;
ir->useMts = (get_eeenum(&inp, "mts", yesno_names, wi) != 0);
if (ir->useMts)
{
- opts->numMtsLevels = get_eint(&inp, "mts-levels", 2, wi);
- ir->mtsLevels.resize(2);
- gmx::MtsLevel& mtsLevel = ir->mtsLevels[1];
- opts->mtsLevel2Forces = setStringEntry(&inp, "mts-level2-forces", "longrange-nonbonded");
- mtsLevel.stepFactor = get_eint(&inp, "mts-level2-factor", 2, wi);
+ gmx::GromppMtsOpts& mtsOpts = opts->mtsOpts;
+ mtsOpts.numLevels = get_eint(&inp, "mts-levels", 2, wi);
- ir->mtsLevels.resize(2);
- mtsOpts.level2Forces = setStringEntry(
- &inp, "mts-level2-forces", "longrange-nonbonded nonbonded pair dihedral");
++ mtsOpts.level2Forces = setStringEntry(&inp, "mts-level2-forces", "longrange-nonbonded");
+ mtsOpts.level2Factor = get_eint(&inp, "mts-level2-factor", 2, wi);
// We clear after reading without dynamics to not force the user to remove MTS mdp options
if (!EI_DYNAMICS(ir->eI))
{
grps->emplace_back(gid);
}
-
+ GMX_ASSERT(block, "Can't have a nullptr block");
+ atomGroupRangeValidation(natoms, gid, *block);
/* Now go over the atoms in the group */
- for (j = block->index[gid]; (j < block->index[gid + 1]); j++)
+ for (int j = block->index[gid]; (j < block->index[gid + 1]); j++)
{
-
- aj = block->a[j];
-
- /* Range checking */
- if ((aj < 0) || (aj >= natoms))
- {
- gmx_fatal(FARGS, "Invalid atom number %d in indexfile", aj + 1);
- }
+ const int aj = block->a[j];
/* Lookup up the old group number */
- ognr = cbuf[aj];
+ const int ognr = cbuf[aj];
if (ognr != NOGID)
{
- gmx_fatal(FARGS, "Atom %d in multiple %s groups (%d and %d)", aj + 1, title,
- ognr + 1, i + 1);
+ gmx_fatal(FARGS, "Atom %d in multiple %s groups (%d and %d)", aj + 1, title, ognr + 1, i + 1);
}
else
{
if (ir->bPull)
{
- const int gid = search_string(inputrecStrings->pullGroupNames[i].c_str(),
- defaultIndexGroups->nr, gnames);
+ for (int i = 1; i < ir->pull->ngroup; i++)
+ {
++ const int gid = search_string(
++ inputrecStrings->pullGroupNames[i].c_str(), defaultIndexGroups->nr, gnames);
+ GMX_ASSERT(defaultIndexGroups, "Must have initialized default index groups");
+ atomGroupRangeValidation(natoms, gid, *defaultIndexGroups);
+ }
+
process_pull_groups(ir->pull->group, inputrecStrings->pullGroupNames, defaultIndexGroups, gnames);
checkPullCoords(ir->pull->group, ir->pull->coord);
/* Must correspond to the Directive enum in grompp_impl.h */
static gmx::EnumerationArray<Directive, const char*> directive_names = {
- { "defaults", "atomtypes", "bondtypes", "constrainttypes", "pairtypes", "angletypes",
- "dihedraltypes", "nonbond_params", "implicit_genborn_params", "implicit_surface_params",
+ { "defaults",
+ "atomtypes",
+ "bondtypes",
+ "constrainttypes",
+ "pairtypes",
+ "angletypes",
+ "dihedraltypes",
+ "nonbond_params",
+ "implicit_genborn_params",
+ "implicit_surface_params",
"cmaptypes",
/* All the directives above can not appear after moleculetype */
- "moleculetype", "atoms", "virtual_sites1", "virtual_sites2", "virtual_sites3",
- "virtual_sites4", "virtual_sitesn", "bonds", "exclusions", "pairs", "pairs_nb", "angles",
- "dihedrals", "constraints", "settles", "polarization", "water_polarization",
- "thole_polarization", "system", "molecules", "position_restraints", "angle_restraints",
- "angle_restraints_z", "distance_restraints", "orientation_restraints", "dihedral_restraints",
- "cmap", "intermolecular_interactions", "maxdirs", "invalid", "none" }
+ "moleculetype",
+ "atoms",
++ "virtual_sites1",
+ "virtual_sites2",
+ "virtual_sites3",
+ "virtual_sites4",
+ "virtual_sitesn",
+ "bonds",
+ "exclusions",
+ "pairs",
+ "pairs_nb",
+ "angles",
+ "dihedrals",
+ "constraints",
+ "settles",
+ "polarization",
+ "water_polarization",
+ "thole_polarization",
+ "system",
+ "molecules",
+ "position_restraints",
+ "angle_restraints",
+ "angle_restraints_z",
+ "distance_restraints",
+ "orientation_restraints",
+ "dihedral_restraints",
+ "cmap",
+ "intermolecular_interactions",
+ "maxdirs",
+ "invalid",
+ "none" }
};
int ifunc_index(Directive d, int type)
GMX_RELEASE_ASSERT(
mi0,
"Need to have a valid MoleculeInformation object to work on");
- push_bond(d, interactions, mi0->interactions, &(mi0->atoms), atypes,
- pline, FALSE, FALSE, 1.0, bZero, &bWarn_copy_A_B, wi);
+ push_bond(d,
+ interactions,
+ mi0->interactions,
+ &(mi0->atoms),
+ atypes,
+ pline,
+ FALSE,
+ FALSE,
+ 1.0,
+ bZero,
+ &bWarn_copy_A_B,
+ wi);
break;
+ case Directive::d_vsites1:
case Directive::d_vsites2:
case Directive::d_vsites3:
case Directive::d_vsites4:
const auto numVoxels = gradient_.asConstView().mapping().required_span_size();
/* the gradient for the inner product measure of fit is constant and does not
* depend on the compared density, so it is pre-computed here */
- std::transform(begin(referenceDensity_), end(referenceDensity), begin(gradient_), [numVoxels](float x) {
- std::transform(begin(referenceDensity_), end(referenceDensity_), begin(gradient_),
- [numVoxels](float x) { return x / numVoxels; });
++ std::transform(begin(referenceDensity_), end(referenceDensity_), begin(gradient_), [numVoxels](float x) {
+ return x / numVoxels;
+ });
}
real DensitySimilarityInnerProduct::similarity(density comparedDensity)
if (!exist_output_file(outputfile.filename, nfile, fnm))
{
writer.writeLine(outputfile.filename);
+ // If this was a pull file, then we have a known issue and
+ // work-around (See gitlab issue #3442).
+ if (!missingFilesIncludedPullOutputFiles
+ && (contains(outputfile.filename, "pullx")
+ || contains(outputfile.filename, "pullf")))
+ {
+ missingFilesIncludedPullOutputFiles = true;
+ }
}
}
+ if (missingFilesIncludedPullOutputFiles)
+ {
+ writer.ensureEmptyLine();
+ writer.writeLineFormatted(
+ "It appears that pull output files were not found. It is known that "
+ "using gmx mdrun -deffnm test with pulling and later "
+ "gmx mdrun -deffnm test -cpi will fail to consider the changed default "
+ "filename when checking the pull output files for restarting with "
+ "appending. You may be able to work around this by using a command like "
+ "gmx mdrun -deffnm test -px test_pullx -pf test_pullf -cpi.");
+ }
settings.setIndent(oldIndent);
+ writer.ensureEmptyLine();
writer.writeLineFormatted(
- R"(To keep your simulation files safe, this simulation will not restart. Either name your
- output files exactly the same as the previous simulation part (e.g. with -deffnm), or
- make sure all the output files are present (e.g. run from the same directory as the
- previous simulation part), or instruct mdrun to write new output files with mdrun -noappend.
- In the last case, you will not be able to use appending in future for this simulation.)",
+ "To keep your simulation files safe, this simulation will not restart. "
+ "Either name your output files exactly the same as the previous simulation "
+ "part (e.g. with -deffnm or explicit naming), or make sure all the output "
+ "files are present (e.g. run from the same directory as the previous simulation "
+ "part), or instruct mdrun to write new output files with mdrun -noappend. In "
+ "the last case, you will not be able to use appending in future for this "
+ "simulation.",
- numFilesMissing, outputfiles.size());
+ numFilesMissing,
+ outputfiles.size());
GMX_THROW(InconsistentInputError(stream.toString()));
}
const auto* inputrec = legacySimulatorData_->inputrec;
addSignaller(energySignallerBuilder_.build(
inputrec->nstcalcenergy, inputrec->fepvals->nstdhdl, inputrec->nstpcouple));
- addSignaller(trajectorySignallerBuilder_.build(
- inputrec->nstxout, inputrec->nstvout, inputrec->nstfout,
- inputrec->nstxout_compressed, trajectoryElement->tngBoxOut(),
- trajectoryElement->tngLambdaOut(), trajectoryElement->tngBoxOutCompressed(),
- trajectoryElement->tngLambdaOutCompressed(), inputrec->nstenergy));
- addSignaller(loggingSignallerBuilder_.build(inputrec->nstlog, inputrec->init_step,
- legacySimulatorData_->startingBehavior));
- addSignaller(lastStepSignallerBuilder_.build(inputrec->nsteps, inputrec->init_step,
- algorithm.stopHandler_.get()));
- addSignaller(neighborSearchSignallerBuilder_.build(inputrec->nstlist, inputrec->init_step,
- inputrec->init_t));
+ addSignaller(trajectorySignallerBuilder_.build(inputrec->nstxout,
+ inputrec->nstvout,
+ inputrec->nstfout,
+ inputrec->nstxout_compressed,
+ trajectoryElement->tngBoxOut(),
+ trajectoryElement->tngLambdaOut(),
+ trajectoryElement->tngBoxOutCompressed(),
+ trajectoryElement->tngLambdaOutCompressed(),
+ inputrec->nstenergy));
- addSignaller(loggingSignallerBuilder_.build(inputrec->nstlog, inputrec->init_step, inputrec->init_t));
++ addSignaller(loggingSignallerBuilder_.build(
++ inputrec->nstlog, inputrec->init_step, legacySimulatorData_->startingBehavior));
+ addSignaller(lastStepSignallerBuilder_.build(
+ inputrec->nsteps, inputrec->init_step, algorithm.stopHandler_.get()));
+ addSignaller(neighborSearchSignallerBuilder_.build(
+ inputrec->nstlist, inputrec->init_step, inputrec->init_t));
}
// Create element list
* combination is rarely used. LJ force-switch with LB rule is more common,
* but gives only 1% speed-up.
*/
- if (ic->vdwtype == evdwCUT)
- {
- switch (ic->vdw_modifier)
- {
- case eintmodNONE:
- case eintmodPOTSHIFT:
- switch (nbatParams.comb_rule)
- {
- case ljcrNONE: nbp->vdwtype = evdwTypeCUT; break;
- case ljcrGEOM: nbp->vdwtype = evdwTypeCUTCOMBGEOM; break;
- case ljcrLB: nbp->vdwtype = evdwTypeCUTCOMBLB; break;
- default:
- gmx_incons(
- "The requested LJ combination rule is not implemented in the CUDA "
- "GPU accelerated kernels!");
- }
- break;
- case eintmodFORCESWITCH: nbp->vdwtype = evdwTypeFSWITCH; break;
- case eintmodPOTSWITCH: nbp->vdwtype = evdwTypePSWITCH; break;
- default:
- gmx_incons(
- "The requested VdW interaction modifier is not implemented in the CUDA GPU "
- "accelerated kernels!");
- }
- }
- else if (ic->vdwtype == evdwPME)
- {
- if (ic->ljpme_comb_rule == ljcrGEOM)
- {
- assert(nbatParams.comb_rule == ljcrGEOM);
- nbp->vdwtype = evdwTypeEWALDGEOM;
- }
- else
- {
- assert(nbatParams.comb_rule == ljcrLB);
- nbp->vdwtype = evdwTypeEWALDLB;
- }
- }
- else
- {
- gmx_incons(
- "The requested VdW type is not implemented in the CUDA GPU accelerated kernels!");
- }
-
- if (ic->eeltype == eelCUT)
- {
- nbp->eeltype = eelTypeCUT;
- }
- else if (EEL_RF(ic->eeltype))
- {
- nbp->eeltype = eelTypeRF;
- }
- else if ((EEL_PME(ic->eeltype) || ic->eeltype == eelEWALD))
- {
- nbp->eeltype = nbnxn_gpu_pick_ewald_kernel_type(*ic, deviceContext.deviceInfo());
- }
- else
- {
- /* Shouldn't happen, as this is checked when choosing Verlet-scheme */
- gmx_incons(
- "The requested electrostatics type is not implemented in the CUDA GPU accelerated "
- "kernels!");
- }
+ nbp->vdwType = nbnxmGpuPickVdwKernelType(ic, nbatParams.comb_rule);
- nbp->elecType = nbnxmGpuPickElectrostaticsKernelType(ic);
++ nbp->elecType = nbnxmGpuPickElectrostaticsKernelType(ic, deviceContext.deviceInfo());
/* generate table for PME */
nbp->coulomb_tab = nullptr;
* This file is part of the GROMACS molecular simulation package.
*
* Copyright (c) 2014,2015,2017,2018,2019,2020, 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.
GPU_FUNC_QUALIFIER
bool gpu_is_kernel_ewald_analytical(const NbnxmGpu gmx_unused* nb) GPU_FUNC_TERM_WITH_RETURN(FALSE);
- enum ElecType nbnxmGpuPickElectrostaticsKernelType(const interaction_const_t gmx_unused* ic)
+/** Return the enum value of electrostatics kernel type for given interaction parameters \p ic. */
+GPU_FUNC_QUALIFIER
++enum ElecType nbnxmGpuPickElectrostaticsKernelType(const interaction_const_t gmx_unused* ic,
++ const DeviceInformation gmx_unused& deviceInfo)
+ GPU_FUNC_TERM_WITH_RETURN(ElecType::Count);
+
+/** Return the enum value of VdW kernel type for given \p ic and \p combRule. */
+GPU_FUNC_QUALIFIER
+enum VdwType nbnxmGpuPickVdwKernelType(const interaction_const_t gmx_unused* ic, int gmx_unused combRule)
+ GPU_FUNC_TERM_WITH_RETURN(VdwType::Count);
+
/** Returns an opaque pointer to the GPU command stream
* Note: CUDA only.
*/
#include "nbnxm_gpu_data_mgmt.h"
+ #include "gromacs/hardware/device_information.h"
+#include "gromacs/mdtypes/interaction_const.h"
#include "gromacs/nbnxm/gpu_data_mgmt.h"
#include "gromacs/timing/gpu_timing.h"
#include "gromacs/utility/cstringutil.h"
}
}
- enum ElecType nbnxn_gpu_pick_ewald_kernel_type(const interaction_const_t& ic)
-int nbnxn_gpu_pick_ewald_kernel_type(const interaction_const_t& ic,
- const DeviceInformation gmx_unused& deviceInfo)
++enum ElecType nbnxn_gpu_pick_ewald_kernel_type(const interaction_const_t& ic,
++ const DeviceInformation gmx_unused& deviceInfo)
{
bool bTwinCut = (ic.rcoulomb != ic.rvdw);
- int kernel_type;
/* Benchmarking/development environment variables to force the use of
analytical or tabulated Ewald kernel. */
set_cutoff_parameters(nbp, ic, nbv->pairlistSets().params());
- nbp->elecType = nbnxn_gpu_pick_ewald_kernel_type(*ic);
- nbp->eeltype = nbnxn_gpu_pick_ewald_kernel_type(*ic, nb->deviceContext_->deviceInfo());
++ nbp->elecType = nbnxn_gpu_pick_ewald_kernel_type(*ic, nb->deviceContext_->deviceInfo());
GMX_RELEASE_ASSERT(ic->coulombEwaldTables, "Need valid Coulomb Ewald correction tables");
init_ewald_coulomb_force_table(*ic->coulombEwaldTables, nbp, *nb->deviceContext_);
bool gpu_is_kernel_ewald_analytical(const NbnxmGpu* nb)
{
- return ((nb->nbparam->eeltype == eelTypeEWALD_ANA) || (nb->nbparam->eeltype == eelTypeEWALD_ANA_TWIN));
+ return ((nb->nbparam->elecType == ElecType::EwaldAna)
+ || (nb->nbparam->elecType == ElecType::EwaldAnaTwin));
+}
+
- enum ElecType nbnxmGpuPickElectrostaticsKernelType(const interaction_const_t* ic)
++enum ElecType nbnxmGpuPickElectrostaticsKernelType(const interaction_const_t* ic,
++ const DeviceInformation& deviceInfo)
+{
+ if (ic->eeltype == eelCUT)
+ {
+ return ElecType::Cut;
+ }
+ else if (EEL_RF(ic->eeltype))
+ {
+ return ElecType::RF;
+ }
+ else if ((EEL_PME(ic->eeltype) || ic->eeltype == eelEWALD))
+ {
- return nbnxn_gpu_pick_ewald_kernel_type(*ic);
++ return nbnxn_gpu_pick_ewald_kernel_type(*ic, deviceInfo);
+ }
+ else
+ {
+ /* Shouldn't happen, as this is checked when choosing Verlet-scheme */
+ GMX_THROW(gmx::InconsistentInputError(
+ gmx::formatString("The requested electrostatics type %s (%d) is not implemented in "
+ "the GPU accelerated kernels!",
+ EELTYPE(ic->eeltype),
+ ic->eeltype)));
+ }
+}
+
+
+enum VdwType nbnxmGpuPickVdwKernelType(const interaction_const_t* ic, int combRule)
+{
+ if (ic->vdwtype == evdwCUT)
+ {
+ switch (ic->vdw_modifier)
+ {
+ case eintmodNONE:
+ case eintmodPOTSHIFT:
+ switch (combRule)
+ {
+ case ljcrNONE: return VdwType::Cut;
+ case ljcrGEOM: return VdwType::CutCombGeom;
+ case ljcrLB: return VdwType::CutCombLB;
+ default:
+ GMX_THROW(gmx::InconsistentInputError(gmx::formatString(
+ "The requested LJ combination rule %s (%d) is not implemented in "
+ "the GPU accelerated kernels!",
+ enum_name(combRule, ljcrNR, c_ljcrNames),
+ combRule)));
+ }
+ case eintmodFORCESWITCH: return VdwType::FSwitch;
+ case eintmodPOTSWITCH: return VdwType::PSwitch;
+ default:
+ GMX_THROW(gmx::InconsistentInputError(
+ gmx::formatString("The requested VdW interaction modifier %s (%d) is not "
+ "implemented in the GPU accelerated kernels!",
+ INTMODIFIER(ic->vdw_modifier),
+ ic->vdw_modifier)));
+ }
+ }
+ else if (ic->vdwtype == evdwPME)
+ {
+ if (ic->ljpme_comb_rule == ljcrGEOM)
+ {
+ assert(combRule == ljcrGEOM);
+ return VdwType::EwaldGeom;
+ }
+ else
+ {
+ assert(combRule == ljcrLB);
+ return VdwType::EwaldLB;
+ }
+ }
+ else
+ {
+ GMX_THROW(gmx::InconsistentInputError(gmx::formatString(
+ "The requested VdW type %s (%d) is not implemented in the GPU accelerated kernels!",
+ EVDWTYPE(ic->vdwtype),
+ ic->vdwtype)));
+ }
}
} // namespace Nbnxm
const DeviceContext& deviceContext);
/*! \brief Selects the Ewald kernel type, analytical or tabulated, single or twin cut-off. */
- enum ElecType nbnxn_gpu_pick_ewald_kernel_type(const interaction_const_t gmx_unused& ic);
-int nbnxn_gpu_pick_ewald_kernel_type(const interaction_const_t gmx_unused& ic,
- const DeviceInformation& deviceInfo);
++enum ElecType nbnxn_gpu_pick_ewald_kernel_type(const interaction_const_t gmx_unused& ic,
++ const DeviceInformation& deviceInfo);
/*! \brief Copies all parameters related to the cut-off from ic to nbp
*/
{
set_cutoff_parameters(nbp, ic, listParams);
- map_interaction_types_to_gpu_kernel_flavors(ic, nbatParams.comb_rule, &(nbp->eeltype),
- &(nbp->vdwtype), deviceContext);
+ nbp->vdwType = nbnxmGpuPickVdwKernelType(ic, nbatParams.comb_rule);
- nbp->elecType = nbnxmGpuPickElectrostaticsKernelType(ic);
++ nbp->elecType = nbnxmGpuPickElectrostaticsKernelType(ic, deviceContext.deviceInfo());
if (ic->vdwtype == evdwPME)
{
gmx_fatal(FARGS,
"Distance between pull groups %d and %d (%f nm) is larger than 0.49 times the "
"box size (%f).\n%s",
- pcrd->params.group[0],
- pcrd->params.group[1],
- pcrd->params.group[groupIndex0], pcrd->params.group[groupIndex1], sqrt(dr2),
- sqrt(0.98 * 0.98 * max_dist2), pcrd->params.eGeom == epullgDIR ? "You might want to consider using \"pull-geometry = direction-periodic\" instead.\n" : "");
++ pcrd->params.group[groupIndex0],
++ pcrd->params.group[groupIndex1],
+ sqrt(dr2),
+ sqrt(0.98 * 0.98 * max_dist2),
+ pcrd->params.eGeom == epullgDIR
+ ? "You might want to consider using \"pull-geometry = "
+ "direction-periodic\" instead.\n"
+ : "");
}
if (pcrd->params.eGeom == epullgDIRPBC)
pull_group_work_t* pgrp0 = &pull->group[pcrd->params.group[0]];
pull_group_work_t* pgrp1 = &pull->group[pcrd->params.group[1]];
- low_get_pull_coord_dr(pull, pcrd, pbc, pgrp1->x,
- pcrd->params.eGeom == epullgCYL ? pull->dyna[coord_ind].x : pgrp0->x, 0,
- 1, md2, spatialData.dr01);
+ low_get_pull_coord_dr(pull,
+ pcrd,
+ pbc,
+ pgrp1->x,
+ pcrd->params.eGeom == epullgCYL ? pull->dyna[coord_ind].x : pgrp0->x,
++ 0,
++ 1,
+ md2,
+ spatialData.dr01);
if (pcrd->params.ngroup >= 4)
{
pgrp1 = &pull->group[pcrd->params.group[1]];
/* Get the current difference vector */
- low_get_pull_coord_dr(pull, pcrd, pbc, rnew[pcrd->params.group[1]],
- rnew[pcrd->params.group[0]], 0, 1, -1, unc_ij);
+ low_get_pull_coord_dr(
- pull, pcrd, pbc, rnew[pcrd->params.group[1]], rnew[pcrd->params.group[0]], -1, unc_ij);
++ pull, pcrd, pbc, rnew[pcrd->params.group[1]], rnew[pcrd->params.group[0]], 0, 1, -1, unc_ij);
if (debug)
{
g0 = pcrd->params.group[0];
g1 = pcrd->params.group[1];
- low_get_pull_coord_dr(pull, pcrd, pbc, rnew[g1], rnew[g0], -1, tmp);
- low_get_pull_coord_dr(pull, pcrd, pbc, dr1, dr0, -1, tmp3);
+ low_get_pull_coord_dr(pull, pcrd, pbc, rnew[g1], rnew[g0], 0, 1, -1, tmp);
+ low_get_pull_coord_dr(pull, pcrd, pbc, dr1, dr0, 0, 1, -1, tmp3);
- fprintf(debug, "Pull cur %8.5f %8.5f %8.5f j:%8.5f %8.5f %8.5f d: %8.5f\n", rnew[g0][0],
- rnew[g0][1], rnew[g0][2], rnew[g1][0], rnew[g1][1], rnew[g1][2], dnorm(tmp));
- fprintf(debug, "Pull ref %8s %8s %8s %8s %8s %8s d: %8.5f\n", "", "", "", "", "",
- "", pcrd->value_ref);
- fprintf(debug, "Pull cor %8.5f %8.5f %8.5f j:%8.5f %8.5f %8.5f d: %8.5f\n", dr0[0],
- dr0[1], dr0[2], dr1[0], dr1[1], dr1[2], dnorm(tmp3));
+ fprintf(debug,
+ "Pull cur %8.5f %8.5f %8.5f j:%8.5f %8.5f %8.5f d: %8.5f\n",
+ rnew[g0][0],
+ rnew[g0][1],
+ rnew[g0][2],
+ rnew[g1][0],
+ rnew[g1][1],
+ rnew[g1][2],
+ dnorm(tmp));
+ fprintf(debug,
+ "Pull ref %8s %8s %8s %8s %8s %8s d: %8.5f\n",
+ "",
+ "",
+ "",
+ "",
+ "",
+ "",
+ pcrd->value_ref);
+ fprintf(debug,
+ "Pull cor %8.5f %8.5f %8.5f j:%8.5f %8.5f %8.5f d: %8.5f\n",
+ dr0[0],
+ dr0[1],
+ dr0[2],
+ dr1[0],
+ dr1[1],
+ dr1[2],
+ dnorm(tmp3));
} /* END DEBUG */
/* Update the COMs with dr */
continue;
}
- low_get_pull_coord_dr(pull, &coord, pbc, rnew[coord.params.group[1]],
- rnew[coord.params.group[0]], 0, 1, -1, unc_ij);
+ low_get_pull_coord_dr(
- pull, &coord, pbc, rnew[coord.params.group[1]], rnew[coord.params.group[0]], -1, unc_ij);
++ pull, &coord, pbc, rnew[coord.params.group[1]], rnew[coord.params.group[0]], 0, 1, -1, unc_ij);
switch (coord.params.eGeom)
{
biod.pnpi.spb.ru / alexxy / gromacs.git / commitdiff