diff --git a/examples/02_MultiResArgon_LocalCap/02_MultiResArgon_LocalCap.cpp b/examples/02_MultiResArgon_LocalCap/02_MultiResArgon_LocalCap.cpp new file mode 100644 index 00000000..91f01369 --- /dev/null +++ b/examples/02_MultiResArgon_LocalCap/02_MultiResArgon_LocalCap.cpp @@ -0,0 +1,349 @@ +// Copyright 2024 Sebastian Eibl +// +// Licensed under the Apache License, Version 2.0 (the "License"); +// you may not use this file except in compliance with the License. +// You may obtain a copy of the License at +// +// https://www.apache.org/licenses/LICENSE-2.0 +// +// Unless required by applicable law or agreed to in writing, software +// distributed under the License is distributed on an "AS IS" BASIS, +// WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +// See the License for the specific language governing permissions and +// limitations under the License. + +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include "Cabana_NeighborList.hpp" +#include "action/LangevinThermostat.hpp" +#include "action/LennardJones.hpp" +#include "action/LimitAcceleration.hpp" +#include "action/LimitVelocity.hpp" +#include "action/VelocityVerlet.hpp" +#include "analysis/KineticEnergy.hpp" +#include "analysis/MeanSquareDisplacement.hpp" +#include "analysis/Pressure.hpp" +#include "analysis/SystemMomentum.hpp" +#include "communication/GhostLayer.hpp" +#include "data/Atoms.hpp" +#include "data/Subdomain.hpp" +#include "datatypes.hpp" +#include "initialization.hpp" +#include "util/EnvironmentVariables.hpp" +#include "util/IsInSymmetricSlab.hpp" +#include "util/PrintTable.hpp" + +using namespace mrmd; + +/** + * Configuration for the Argon NVE example simulation. + */ +struct Config +{ + // simulation time parameters + idx_t nsteps = 400001; ///< number of steps to simulate + static constexpr real_t dt = 0.002; ///< time step size in reduced units + + // interaction parameters + static constexpr real_t sigma = + 1_r; ///< distance at which LJ potential is zero in reduced units + static constexpr real_t epsilon = 1_r; ///< energy well depth of LJ potential in reduced units + static constexpr real_t mass = 1_r; ///< mass of one atom in reduced units + static constexpr real_t maxVelocity = + 1_r; ///< maximum initial velocity component in reduced units + static constexpr real_t r_cut = 2.5_r * sigma; ///< cutoff radius for LJ potential + static constexpr real_t r_cap = 0.7_r * sigma; ///< capping radius for LJ potential + + // neighbor list parameters + static constexpr real_t skin = 0.1_r * sigma; ///< skin thickness for neighbor list + static constexpr real_t neighborCutoff = r_cut + skin; ///< cutoff radius for neighbor list + static constexpr real_t cell_ratio = + 1_r; ///< ratio of cell size on Cartesian grid to cutoff radius for neighbor list + static constexpr idx_t estimatedMaxNeighbors = + 60; ///< estimated maximum number of neighbors per atom + + // system parameters + static constexpr idx_t numAtoms = 16 * 16 * 16; ///< number of atoms in the simulation + real_t Lx = 30_r * sigma; ///< box edge length in x-direction + + // equilibration parameters + idx_t nstepsEq = 100000; ///< number of equilibration steps + real_t temperature = + 1.5_r; ///< target temperature during equilibration for thermostat in reduced units + static constexpr real_t gamma = 0.04_r / dt; ///< friction coefficient for Langevin thermostat + + // output parameters + bool bOutput = true; ///< whether to output data files + idx_t outputInterval = -1; ///< interval for data file output (-1: no output) + const std::string resName = "Argon"; ///< residue name for output files + const std::vector typeNames = {"Ar"}; ///< atom type names for output files +}; + +data::Atoms fillDomainWithAtomsSC(const data::Subdomain& subdomain, + const idx_t& numAtoms, + const real_t& maxVelocity) +{ + auto RNG = Kokkos::Random_XorShift1024_Pool<>(1234); + + data::Atoms atoms(numAtoms); + + auto pos = atoms.getPos(); + auto vel = atoms.getVel(); + auto mass = atoms.getMass(); + auto type = atoms.getType(); + auto charge = atoms.getCharge(); + + auto policy = Kokkos::RangePolicy<>(0, numAtoms); + auto kernel = KOKKOS_LAMBDA(const idx_t idx) + { + auto randGen = RNG.get_state(); + pos(idx, 0) = randGen.drand() * subdomain.diameter[0] + subdomain.minCorner[0]; + pos(idx, 1) = randGen.drand() * subdomain.diameter[1] + subdomain.minCorner[1]; + pos(idx, 2) = randGen.drand() * subdomain.diameter[2] + subdomain.minCorner[2]; + + vel(idx, 0) = (randGen.drand() - 0.5_r) * maxVelocity; + vel(idx, 1) = (randGen.drand() - 0.5_r) * maxVelocity; + vel(idx, 2) = (randGen.drand() - 0.5_r) * maxVelocity; + RNG.free_state(randGen); + + mass(idx) = Config::mass; + type(idx) = 0; + charge(idx) = 0_r; + }; + Kokkos::parallel_for("fillDomainWithAtomsSC", policy, kernel); + + atoms.numLocalAtoms = numAtoms; + atoms.numGhostAtoms = 0; + return atoms; +} + +void runSimulation(Config& config) +{ + // initialize simulation domain + auto subdomain = + data::Subdomain({0_r, 0_r, 0_r}, {config.Lx, config.Lx, config.Lx}, config.neighborCutoff); + + // calculate volume of the simulation domain + const auto volume = subdomain.diameter[0] * subdomain.diameter[1] * subdomain.diameter[2]; + + // initialize atoms randomly in the domain + auto atoms = fillDomainWithAtomsSC(subdomain, config.numAtoms, config.maxVelocity); + + // calculate and print initial density + auto rho = real_c(atoms.numLocalAtoms) / volume; + std::cout << "rho: " << rho << std::endl; + + // set up ghost layer for periodic boundary conditions + communication::GhostLayer ghostLayer; + + // set up neighbor list + HalfVerletList verletList; + real_t maxAtomDisplacement = std::numeric_limits::max(); + idx_t rebuildCounter = 0; + + // set up interaction potential and force calculation and application + action::LennardJones lennardJones(config.r_cut, config.sigma, config.epsilon, 0_r); + action::LennardJones lennardJonesCap(config.r_cut, config.sigma, config.epsilon, config.r_cap); + + // calculate and print box center coordinates + real_t boxCenterX = 0.5_r * (subdomain.maxCorner[0] + subdomain.minCorner[0]); + real_t boxCenterY = 0.5_r * (subdomain.maxCorner[1] + subdomain.minCorner[1]); + real_t boxCenterZ = 0.5_r * (subdomain.maxCorner[2] + subdomain.minCorner[2]); + + std::cout << "x center: " << boxCenterX << std::endl; + std::cout << "y center: " << boxCenterY << std::endl; + std::cout << "z center: " << boxCenterZ << std::endl; + + // set up different interaction regions for capped and bare LJ potential + util::IsInSymmetricSlab isInCentralRegion( + {boxCenterX, boxCenterY, boxCenterZ}, 0_r, 20_r * config.sigma); + util::IsInSymmetricSlab isInCappingRegion( + {boxCenterX, boxCenterY, boxCenterZ}, 20_r * config.sigma, 25_r * config.sigma); + + // set up thermostat for temperature control during equilibration + action::LangevinThermostat langevinThermostat(config.gamma, config.temperature, config.dt); + + // set up timer for runtime measurement + Kokkos::Timer timer; + + // set up mean square displacement analysis + analysis::MeanSquareDisplacement meanSquareDisplacement; + meanSquareDisplacement.reset(atoms); + auto msd = 0_r; + + // print table header for simulation statistics + util::printTable("step", "time", "T", "Ek", "E0", "E", "p", "msd", "Nlocal", "Nghost"); + util::printTableSep("step", "time", "T", "Ek", "E0", "E", "p", "msd", "Nlocal", "Nghost"); + + // open statistics file for writing simulation statistics + std::ofstream fStat("statistics.txt"); + + // main simulation loop + for (auto step = 0; step < config.nsteps; ++step) + { + // integrate equations of motion - first half step (drift) + maxAtomDisplacement += action::VelocityVerlet::preForceIntegrate(atoms, config.dt); + + // reinsert atoms that left the domain according to periodic boundary conditions + ghostLayer.exchangeRealAtoms(atoms, subdomain); + + // create ghost atoms in the ghost layer beyond the periodic boundaries + ghostLayer.createGhostAtoms(atoms, subdomain); + + // check if neighbor list needs to be rebuilt + if (maxAtomDisplacement >= + config.skin * + 0.5_r) // the condition is on half the skin thickness because in principle two + // atoms may both move half the skin thickness towards each other + { + // reset displacement + maxAtomDisplacement = 0_r; + + // rebuild neighbor list + verletList.build(atoms.getPos(), + 0, + atoms.numLocalAtoms, + config.neighborCutoff, + config.cell_ratio, + subdomain.minGhostCorner.data(), + subdomain.maxGhostCorner.data(), + config.estimatedMaxNeighbors); + ++rebuildCounter; + } + + // reset forces to zero + auto force = atoms.getForce(); + Cabana::deep_copy(force, 0_r); + + // check if still during equilibration phase + if (step <= config.nstepsEq) + { + // apply capped LJ potential in the whole domain during equilibration + lennardJonesCap.apply(atoms, verletList); + } + else + { + // compute and apply forces + lennardJones.apply_if( + atoms, + verletList, + KOKKOS_LAMBDA(const real_t x1, + const real_t y1, + const real_t z1, + const real_t x2, + const real_t y2, + const real_t z2) { + return isInCentralRegion(x1, y1, z1) || isInCentralRegion(x2, y2, z2); + }); + lennardJonesCap.apply_if( + atoms, + verletList, + KOKKOS_LAMBDA(const real_t x1, + const real_t y1, + const real_t z1, + const real_t x2, + const real_t y2, + const real_t z2) { + return isInCappingRegion(x1, y1, z1) && isInCappingRegion(x2, y2, z2); + }); + } + + // contribute forces calculated on ghost atoms back to real atoms + ghostLayer.contributeBackGhostToReal(atoms); + + // handle output and statistics + if (config.bOutput && (step % config.outputInterval == 0)) + { + // calculate statistics + auto E0 = (lennardJones.getEnergy() + lennardJonesCap.getEnergy()) / + real_c(atoms.numLocalAtoms); + auto Ek = analysis::getMeanKineticEnergy(atoms); + auto systemMomentum = analysis::getSystemMomentum(atoms); + auto T = (2_r / 3_r) * Ek; + auto p = analysis::getPressure(atoms, subdomain); + msd = meanSquareDisplacement.calc(atoms, subdomain) / + (real_c(config.outputInterval) * config.dt); + meanSquareDisplacement.reset(atoms); + + // print statistics to console + util::printTable(step, + timer.seconds(), + T, + Ek, + E0, + E0 + Ek, + p, + msd, + atoms.numLocalAtoms, + atoms.numGhostAtoms); + + // dump statistics to file + fStat << step << " " << timer.seconds() << " " << T << " " << Ek << " " << E0 << " " + << E0 + Ek << " " << p << " " << msd << " " << atoms.numLocalAtoms << " " + << atoms.numGhostAtoms << " " << std::endl; + } + + // check if still during equilibration phase + if (step <= config.nstepsEq) + { + // apply Langevin thermostat for temperature control during equilibration + langevinThermostat.apply(atoms); + } + + // integrate equations of motion - second half step (kick) + action::VelocityVerlet::postForceIntegrate(atoms, config.dt); + } + + // close statistics file + fStat.close(); + auto time = timer.seconds(); + std::cout << time << std::endl; + + // write performance data to file + auto cores = util::getEnvironmentVariable("OMP_NUM_THREADS"); + std::ofstream fout("ecab.perf", std::ofstream::app); + fout << cores << ", " << time << ", " << atoms.numLocalAtoms << ", " << config.nsteps + << std::endl; + fout.close(); +} + +int main(int argc, char* argv[]) // NOLINT +{ + // initialize Kokkos and MPI environment + initialize(argc, argv); + + // print Kokkos execution space + std::cout << "execution space: " << typeid(Kokkos::DefaultExecutionSpace).name() << std::endl; + + // initialize simulation configuration with command line interface + Config config; + CLI::App app{"Lennard Jones Fluid benchmark application"}; + app.add_option("-n,--nsteps", config.nsteps, "total number of simulation steps"); + app.add_option("-L,--length", config.Lx, "simulation box diameter"); + app.add_option("-e,--numeq", config.nstepsEq, "number of equilibration steps"); + app.add_option( + "-T,--temperature", + config.temperature, + "temperature of the Langevin thermostat (negative numbers deactivate the thermostat)"); + app.add_option("-o,--output", config.outputInterval, "output interval"); + CLI11_PARSE(app, argc, argv); + + // reset output parameter if output interval is negative + if (config.outputInterval < 0) config.bOutput = false; + + // set up, equilibrate in NVT and run production in NVE + runSimulation(config); + + // finalize Kokkos and MPI environment + finalize(); + + return EXIT_SUCCESS; +} \ No newline at end of file diff --git a/examples/02_MultiResArgon_LocalCap/CMakeLists.txt b/examples/02_MultiResArgon_LocalCap/CMakeLists.txt new file mode 100644 index 00000000..cc3722e2 --- /dev/null +++ b/examples/02_MultiResArgon_LocalCap/CMakeLists.txt @@ -0,0 +1,21 @@ +# Copyright 2024 Sebastian Eibl +# +# Licensed under the Apache License, Version 2.0 (the "License"); +# you may not use this file except in compliance with the License. +# You may obtain a copy of the License at +# +# https://www.apache.org/licenses/LICENSE-2.0 +# +# Unless required by applicable law or agreed to in writing, software +# distributed under the License is distributed on an "AS IS" BASIS, +# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +# See the License for the specific language governing permissions and +# limitations under the License. + +add_executable(02_MultiResArgon_LocalCap + 02_MultiResArgon_LocalCap.cpp + ) +target_link_libraries(02_MultiResArgon_LocalCap + PRIVATE mrmd CLI11::CLI11 + ) +add_test(NAME Example.02_MultiResArgon_LocalCap COMMAND ./02_MultiResArgon_LocalCap --nsteps 10) \ No newline at end of file diff --git a/examples/CMakeLists.txt b/examples/CMakeLists.txt index 24b3f3db..861c09c8 100644 --- a/examples/CMakeLists.txt +++ b/examples/CMakeLists.txt @@ -12,6 +12,8 @@ # See the License for the specific language governing permissions and # limitations under the License. +add_subdirectory(02_MultiResArgon_LocalCap) + add_subdirectory(Argon) add_subdirectory(BinaryLennardJones) add_subdirectory(ComplexAtomics) diff --git a/mrmd/action/LennardJones.cpp b/mrmd/action/LennardJones.cpp index 8f49e5ff..db9f11d7 100644 --- a/mrmd/action/LennardJones.cpp +++ b/mrmd/action/LennardJones.cpp @@ -22,60 +22,13 @@ namespace mrmd::action { void LennardJones::apply(data::Atoms& atoms, HalfVerletList& verletList) { - energyAndVirial_ = data::EnergyAndVirialReducer(); - - pos_ = atoms.getPos(); - force_ = atoms.getForce(); - type_ = atoms.getType(); - verletList_ = verletList; - - auto policy = Kokkos::RangePolicy<>(0, atoms.numLocalAtoms); - Kokkos::parallel_reduce("LennardJones::applyForces", policy, *this, energyAndVirial_); - Kokkos::fence(); -} - -KOKKOS_FUNCTION -void LennardJones::operator()(const idx_t& idx, data::EnergyAndVirialReducer& energyAndVirial) const -{ - real_t posTmp[3]; - posTmp[0] = pos_(idx, 0); - posTmp[1] = pos_(idx, 1); - posTmp[2] = pos_(idx, 2); - - real_t forceTmp[3] = {0_r, 0_r, 0_r}; - - const auto numNeighbors = idx_c(HalfNeighborList::numNeighbor(verletList_, idx)); - for (idx_t n = 0; n < numNeighbors; ++n) - { - idx_t jdx = idx_c(HalfNeighborList::getNeighbor(verletList_, idx, n)); - assert(0 <= jdx); - - auto dx = posTmp[0] - pos_(jdx, 0); - auto dy = posTmp[1] - pos_(jdx, 1); - auto dz = posTmp[2] - pos_(jdx, 2); - - auto distSqr = dx * dx + dy * dy + dz * dz; - - if (distSqr > rcSqr_) continue; - - auto typeIdx = type_(idx) * numTypes_ + type_(jdx); - auto forceAndEnergy = LJ_.computeForceAndEnergy(distSqr, typeIdx); - assert(!std::isnan(forceAndEnergy.forceFactor)); - energyAndVirial.energy += forceAndEnergy.energy; - energyAndVirial.virial -= 0.5_r * forceAndEnergy.forceFactor * distSqr; - - forceTmp[0] += dx * forceAndEnergy.forceFactor; - forceTmp[1] += dy * forceAndEnergy.forceFactor; - forceTmp[2] += dz * forceAndEnergy.forceFactor; - - force_(jdx, 0) -= dx * forceAndEnergy.forceFactor; - force_(jdx, 1) -= dy * forceAndEnergy.forceFactor; - force_(jdx, 2) -= dz * forceAndEnergy.forceFactor; - } - - force_(idx, 0) += forceTmp[0]; - force_(idx, 1) += forceTmp[1]; - force_(idx, 2) += forceTmp[2]; + apply_if( + atoms, + verletList, + KOKKOS_LAMBDA( + const real_t, const real_t, const real_t, const real_t, const real_t, const real_t) { + return true; + }); } real_t LennardJones::getEnergy() const { return energyAndVirial_.energy; } diff --git a/mrmd/action/LennardJones.hpp b/mrmd/action/LennardJones.hpp index 4ab5cd94..5b6c0e43 100644 --- a/mrmd/action/LennardJones.hpp +++ b/mrmd/action/LennardJones.hpp @@ -85,14 +85,21 @@ class LennardJones data::EnergyAndVirialReducer energyAndVirial_; public: - KOKKOS_FUNCTION - void operator()(const idx_t& idx, data::EnergyAndVirialReducer& energyAndVirial) const; - real_t getEnergy() const; real_t getVirial() const; void apply(data::Atoms& atoms, HalfVerletList& verletList); + template BinaryPred> + void apply_if(const data::Atoms& atoms, + const HalfVerletList& verletList, + const BinaryPred& pred); + LennardJones(const real_t rc, const real_t& sigma, const real_t& epsilon, @@ -105,4 +112,83 @@ class LennardJones const idx_t& numTypes, const bool isShifted); }; + +template BinaryPred> +void LennardJones::apply_if(const data::Atoms& atoms, + const HalfVerletList& verletList, + const BinaryPred& pred) +{ + energyAndVirial_ = data::EnergyAndVirialReducer(); + pos_ = atoms.getPos(); + force_ = atoms.getForce(); + type_ = atoms.getType(); + verletList_ = verletList; + + auto policy = Kokkos::RangePolicy<>(0, atoms.numLocalAtoms); + + // avoid capturing this pointer + auto pos = pos_; + auto force = force_; + auto type = type_; + auto verletListLocal = verletList_; + auto rcSqr = rcSqr_; + auto LJ = LJ_; + auto numTypes = numTypes_; + auto predLocal = pred; + + auto kernel = KOKKOS_LAMBDA(const idx_t idx, data::EnergyAndVirialReducer& energyAndVirial) + { + real_t posTmp[3]; + posTmp[0] = pos(idx, 0); + posTmp[1] = pos(idx, 1); + posTmp[2] = pos(idx, 2); + + real_t forceTmp[3] = {0_r, 0_r, 0_r}; + + const auto numNeighbors = idx_c(HalfNeighborList::numNeighbor(verletListLocal, idx)); + for (idx_t n = 0; n < numNeighbors; ++n) + { + idx_t jdx = idx_c(HalfNeighborList::getNeighbor(verletListLocal, idx, n)); + assert(0 <= jdx); + + if (!predLocal(posTmp[0], posTmp[1], posTmp[2], pos(jdx, 0), pos(jdx, 1), pos(jdx, 2))) + continue; + + auto dx = posTmp[0] - pos(jdx, 0); + auto dy = posTmp[1] - pos(jdx, 1); + auto dz = posTmp[2] - pos(jdx, 2); + + auto distSqr = dx * dx + dy * dy + dz * dz; + + if (distSqr > rcSqr) continue; + + auto typeIdx = type(idx) * numTypes + type(jdx); + auto forceAndEnergy = LJ.computeForceAndEnergy(distSqr, typeIdx); + assert(!std::isnan(forceAndEnergy.forceFactor)); + energyAndVirial.energy += forceAndEnergy.energy; + energyAndVirial.virial -= 0.5_r * forceAndEnergy.forceFactor * distSqr; + + forceTmp[0] += dx * forceAndEnergy.forceFactor; + forceTmp[1] += dy * forceAndEnergy.forceFactor; + forceTmp[2] += dz * forceAndEnergy.forceFactor; + + force(jdx, 0) -= dx * forceAndEnergy.forceFactor; + force(jdx, 1) -= dy * forceAndEnergy.forceFactor; + force(jdx, 2) -= dz * forceAndEnergy.forceFactor; + } + + force(idx, 0) += forceTmp[0]; + force(idx, 1) += forceTmp[1]; + force(idx, 2) += forceTmp[2]; + }; + + Kokkos::parallel_reduce("LennardJones::apply_if", policy, kernel, energyAndVirial_); + Kokkos::fence(); +} + } // namespace mrmd::action \ No newline at end of file