diff --git a/.devcontainer/Readme.md b/.devcontainer/Readme.md new file mode 100644 index 00000000000..df9dac998d8 --- /dev/null +++ b/.devcontainer/Readme.md @@ -0,0 +1,53 @@ +# ESPResSo in Codespaces + +## Creating a new codespace + +Direct link to the default codespace: https://codespaces.new/espressomd/espresso + +If you need to create a custom codespace, follow these instructions: + +* go to the ESPResSo repository, either the [official project](https://github.com/espressomd/espresso) or your forked project +* follow the GitHub documentation on [Creating a codespace for a repository](https://docs.github.com/en/codespaces/developing-in-a-codespace/creating-a-codespace-for-a-repository#creating-a-codespace-for-a-repository) + +## Building ESPResSo + +General workflow: + +* load MPI module +* build ESPResSo +* create Python environment + +Command lines: + +```sh +# load MPI module +. /etc/profile.d/modules.sh +module load mpi +# build ESPResSo +mkdir build +cd build +cmake .. -D ESPRESSO_BUILD_WITH_FFTW=ON -D ESPRESSO_BUILD_WITH_WALBERLA=ON +make -j$(nproc) +# create Python environment +python -m venv --system-site-packages ESPResSo +. ESPResSo/bin/activate +realpath src/python > $(python -c 'import sysconfig;print(sysconfig.get_path("platlib"))')/espresso.pth +pip install -c ../requirements.txt matplotlib pint pandas tqdm +``` + +## Running interactive notebooks + +Two options: + +* run tutorials from subfolder `doc/tutorials` in the top-level directory +* run `make -j$(nproc) tutorials` to generate notebooks with hidden answers, + and then run tutorials from subfolder `build/doc/tutorials` + +The IDE needs to know which Python interpreter to use: + +* open the Command Palette with F1 (or click on the gears icon) +* search for "Python: Set Project Environment" +* hit "Browse" +* paste `/workspaces/espresso/build/ESPResSo/bin/activate` + +When the IDE asks for the interpreter, select the `ESPResSo` environment. diff --git a/.devcontainer/devcontainer.json b/.devcontainer/devcontainer.json new file mode 100644 index 00000000000..b5532de1fd1 --- /dev/null +++ b/.devcontainer/devcontainer.json @@ -0,0 +1,26 @@ +{ + "image": "ghcr.io/espressomd/docker/fedora:3388b293a5e8a9b613b1bfa4220e54dfb4bff245", + "customizations": { + "vscode": { + "settings": { + "workbench.editorAssociations": { + "*.md": "vscode.markdown.preview.editor" + } + }, + "extensions": [ + "ms-python.python", + "ms-python.vscode-pylance", + "ms-toolsai.jupyter", + "ms-toolsai.jupyter-keymap", + "ms-toolsai.jupyter-renderers" + ] + }, + "codespaces": { + "openFiles": [ + ".devcontainer/Readme.md", + "doc/tutorials/Readme.md" + ] + } + }, + "features": {} +} diff --git a/.github/actions/build_and_check/action.yml b/.github/actions/build_and_check/action.yml index cadb2890ee8..a17d9184a0b 100644 --- a/.github/actions/build_and_check/action.yml +++ b/.github/actions/build_and_check/action.yml @@ -4,7 +4,7 @@ runs: using: "composite" steps: - run: | - export myconfig=maxset with_cuda=false with_gsl=false with_hdf5=false test_timeout=800 check_skip_long=true + export myconfig=maxset with_cuda=false with_hdf5=false test_timeout=800 check_skip_long=true if test -d venv; then . venv/bin/activate; fi bash maintainer/CI/build_cmake.sh shell: bash diff --git a/.github/workflows/push_pull.yml b/.github/workflows/push_pull.yml index b3b037c708d..0334a9800e8 100644 --- a/.github/workflows/push_pull.yml +++ b/.github/workflows/push_pull.yml @@ -24,7 +24,7 @@ jobs: - name: Install dependencies run: | echo '' > Brewfile - for dep in boost boost-mpi fftw ccache ninja; do echo "brew \"${dep}\"" >> Brewfile; done + for dep in boost boost-mpi fftw gsl ccache ninja; do echo "brew \"${dep}\"" >> Brewfile; done brew bundle check >/dev/null 2>&1 || brew bundle --verbose --no-upgrade # set up Xcode compiler echo "::group::List available Xcode compilers" @@ -51,6 +51,7 @@ jobs: cmake_params: '-D ESPRESSO_TEST_NP=3' with_ccache: 'true' with_fftw: 'true' + with_gsl: 'true' with_walberla: 'true' with_walberla_avx: 'false' with_shared_memory_parallelism: 'false' diff --git a/.gitpod.yml b/.gitpod.yml deleted file mode 100644 index d269613cfd1..00000000000 --- a/.gitpod.yml +++ /dev/null @@ -1,31 +0,0 @@ -tasks: - - name: requirements - init: | - sudo apt-get update - sudo apt-get install -y --no-install-recommends build-essential cmake ffmpeg libboost-all-dev libfftw3-dev libgsl-dev libpython3-dev openmpi-common - pip install --no-cache cython numpy scipy matplotlib pint tqdm --constraint requirements.txt - pip install --no-cache nbformat notebook jupyterlab 'jupyter_contrib_nbextensions==0.5.1' - jupyter contrib nbextension install --user - jupyter nbextension enable rubberband/main - jupyter nbextension enable exercise2/main - mkdir -p ~/.jupyter/lab/user-settings/\@jupyterlab/docmanager-extension - mkdir -p ${GITPOD_REPO_ROOT}/.vscode - echo '{"nbext_hide_incompat": false}' > ~/.jupyter/nbconfig/common.json - echo '{"defaultViewers": {"markdown": "Markdown Preview"}}' > ~/.jupyter/lab/user-settings/\@jupyterlab/docmanager-extension/plugin.jupyterlab-settings - echo '{"python.defaultInterpreterPath": "build/pypresso"}' > ${GITPOD_REPO_ROOT}/.vscode/settings.json - echo 'export PYTHONPATH="${PYTHONPATH:+$PYTHONPATH:}${GITPOD_REPO_ROOT}/build/pypresso"' >> ~/.bashrc - command: | - mkdir build - cd build - cmake .. -D ESPRESSO_BUILD_WITH_GSL=ON -D ESPRESSO_BUILD_WITH_WALBERLA=ON - make - make tutorials - for f in doc/tutorials/*/*.ipynb; do echo "configuring ${f}"; ./pypresso doc/tutorials/convert.py exercise2 --to-vscode-jupyter ${f}; done - sed -i 's/espressomd.lb.LBFluidWalberlaGPU/espressomd.lb.LBFluidWalberla/; s/, \\"CUDA\\"\]/]/;' doc/tutorials/lattice_boltzmann/lattice_boltzmann_poiseuille_flow.ipynb doc/tutorials/active_matter/active_matter.ipynb - -vscode: - extensions: - - ms-python.python - - ms-toolsai.jupyter - - ms-toolsai.jupyter-keymap - - ms-toolsai.jupyter-renderers diff --git a/CMakeLists.txt b/CMakeLists.txt index 9b6f4ffded6..c74fb4241d5 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -221,6 +221,7 @@ endmacro() espresso_minimal_compiler_version("GNU" 12.2.0) espresso_minimal_compiler_version("Clang" 18.1.0) espresso_minimal_compiler_version("AppleClang" 17.0.0) +espresso_minimal_compiler_version("CrayClang" 17.0.0) espresso_minimal_compiler_version("IntelLLVM" 2023.1) set(ESPRESSO_UNSUPPORTED_COMPILERS "Intel;MSVC") @@ -588,7 +589,7 @@ if(ESPRESSO_BUILD_WITH_FFTW) endif() if(ESPRESSO_BUILD_WITH_SHARED_MEMORY_PARALLELISM) - find_package(OpenMP REQUIRED) + find_package(OpenMP REQUIRED COMPONENTS CXX) if(NOT EXISTS ${FETCHCONTENT_BASE_DIR}/kokkos-src) find_package(Kokkos 4.3 QUIET) diff --git a/Readme.md b/Readme.md index d40c89dca91..df0c792bf08 100644 --- a/Readme.md +++ b/Readme.md @@ -3,7 +3,7 @@ [![GitLab CI](https://gitlab.icp.uni-stuttgart.de/espressomd/espresso/badges/python/pipeline.svg)](https://gitlab.icp.uni-stuttgart.de/espressomd/espresso/-/commits/python) [![codecov](https://codecov.io/gh/espressomd/espresso/branch/python/graph/badge.svg)](https://codecov.io/gh/espressomd/espresso) [![Binder](https://mybinder.org/badge_logo.svg)](https://mybinder.org/v2/gh/jngrad/espresso-binder/HEAD) -[![Contribute with Gitpod](https://img.shields.io/badge/Contribute%20with-Gitpod-908a85?logo=gitpod)](https://gitpod.io/#https://github.com/espressomd/espresso) +[![Open in GitHub Codespaces](https://img.shields.io/badge/launch-Codespaces-579aca?logo=github)](https://codespaces.new/espressomd/espresso) This is the Molecular Dynamics software ESPResSo ("Extensible Simulation Package for Research on Soft Matter Systems"). @@ -86,9 +86,9 @@ Common installation issues are addressed in the [FAQ](https://github.com/espressomd/espresso/wiki/Installation-FAQ), together with contributed patches for compiler-related and library-related issues. -You can also try ESPReSSo in the cloud using +You can also try ESPResSo in the cloud using [Binder](https://mybinder.org/v2/gh/jngrad/espresso-binder/HEAD) or -[Gitpod](https://gitpod.io/#https://github.com/espressomd/espresso). +[Codespaces](https://codespaces.new/espressomd/espresso). ### Join the community diff --git a/doc/bibliography.bib b/doc/bibliography.bib index 23f601989a0..ed88503c382 100644 --- a/doc/bibliography.bib +++ b/doc/bibliography.bib @@ -719,6 +719,16 @@ @Article{johnson94a publisher = {Taylor \& Francis}, } +@InProceedings{jupyter18a, + author = {{Project Jupyter} and Bussonnier, Matthias and Forde, Jessica and Freeman, Jeremy and Granger, Brian and Head, Tim and Holdgraf, Chris and Kelley, Kyle and Nalvarte, Gladys and Osheroff, Andrew and Pacer, M and Panda, Yuvi and Perez, Fernando and Ragan-Kelley, Benjamin and Willing, Carol}, + title = {{Binder}~2.0 - Reproducible, interactive, sharable environments for science at scale}, + booktitle = {Proceedings of the 17\textsuperscript{th} Python in Science Conference}, + year = {2018}, + editor = {Akici, Fatih and Lippa, David and Niederhut, Dillon and Pacer, M}, + pages = {113--120}, + doi = {10.25080/Majora-4af1f417-011}, +} + @Article{kesselheim11a, author = {Kesselheim, Stefan and Sega, Marcello and Holm, Christian}, title = {Applying {ICC}* to {DNA} translocation. {E}ffect of dielectric boundaries}, diff --git a/doc/sphinx/constraints.rst b/doc/sphinx/constraints.rst index 6bc6bd13922..2d31f81a41b 100644 --- a/doc/sphinx/constraints.rst +++ b/doc/sphinx/constraints.rst @@ -80,7 +80,7 @@ invoke the :meth:`~espressomd.constraints.Constraints.add` method:: All previously listed shapes can be added to the system constraints by passing an initialized shape object to :meth:`~espressomd.constraints.Constraints.add`, returning a constraint object :: - misshaped = Wall(dist=20, normal=[0.1, 0.0, 1]) + myShape = Wall(dist=20, normal=[0.1, 0.0, 1]) myConstraint = system.constraints.add(shape=myShape, particle_type=p_type) The extra argument ``particle_type`` specifies the non-bonded interaction to be used with @@ -169,17 +169,10 @@ initial configurations. Available shapes ^^^^^^^^^^^^^^^^ -:class:`espressomd.shapes` +:mod:`espressomd.shapes` -Python syntax:: - - import espressomd.shapes - shape = espressomd.shapes. - -```` can be any of the available shapes. - -The surface's geometry is defined via a few available shapes. -The following shapes can be used as constraints. +A shape-based constraint surface's geometry is defined by a shape object, +or by a collection thereof (:ref:`Shape union`). .. warning:: When using shapes with concave edges and corners, the fact that a particle @@ -210,6 +203,17 @@ Pictured is an example constraint with a ``Wall`` shape created with :: wall = Wall(dist=20, normal=[0.1, 0.0, 1]) system.constraints.add(shape=wall, particle_type=0) +When placing two non-penetrable parallel walls in a box, e.g. to simulate a plate capacitor, +their normal vectors must be collinear but opposite in direction. +The distance calculation is affected by the sign of the normal vector, like so:: + + # two parallel plates oriented such that particles can only be found + # on the z-axis in the range [1.5, box_l - 1.5] + system.constraints.add(shape=espressomd.shapes.Wall( + normal=[0, 0, +1], dist=1.5), particle_type=1) + system.constraints.add(shape=espressomd.shapes.Wall( + normal=[0, 0, -1], dist=-(system.box_l - 1.5)), particle_type=1) + For penetrable walls, if the ``only_positive`` flag is set to ``True``, interactions are only calculated if the particle is on the side of the wall in which the normal vector is pointing. @@ -463,6 +467,8 @@ Note: in the OpenGL visualizer, if the OpenGL Extrusion library is not available the shape surface will be rendered with dots. +.. _Shape union: + Union """"" @@ -471,7 +477,6 @@ Union A meta-shape which is the union of given shapes. Note that only the regions where all shapes have a "positive distance" (see :ref:`Available options`) can be used for the union. The distance to the union is defined as the minimum distance to any contained shape. -This shape cannot be checkpointed when multiple MPI ranks are used. .. _Available options: diff --git a/doc/sphinx/electrostatics.rst b/doc/sphinx/electrostatics.rst index a77170d928e..710aee1c2ea 100644 --- a/doc/sphinx/electrostatics.rst +++ b/doc/sphinx/electrostatics.rst @@ -67,6 +67,12 @@ Coulomb P3M :class:`espressomd.electrostatics.P3M` +.. note:: + + Requires feature ``ELECTROSTATICS`` and + external feature ``FFTW``, enabled with + ``-D ESPRESSO_BUILD_WITH_FFTW=ON``. + For this feature to work, you need to have the ``fftw3`` library installed on your system. In |es|, you can check if it is compiled in by checking for the feature ``FFTW`` with ``espressomd.features``. @@ -115,6 +121,12 @@ Coulomb P3M on GPU :class:`espressomd.electrostatics.P3MGPU` +.. note:: + + Requires feature ``ELECTROSTATICS`` and + external features ``FFTW`` and ``CUDA``, enabled with + ``-D ESPRESSO_BUILD_WITH_FFTW=ON -D ESPRESSO_BUILD_WITH_CUDA=ON``. + The GPU implementation of P3M calculates the far field contribution to the forces on the GPU. The near-field contribution to the forces, as well as the near- and far-field contributions to the energies are calculated on the CPU. @@ -136,6 +148,10 @@ Debye-Hückel potential :class:`espressomd.electrostatics.DH` +.. note:: + + Requires feature ``ELECTROSTATICS``. + The Debye-Hückel electrostatic potential is defined by .. math:: U^{C-DH} = C \cdot \frac{q_1 q_2 \exp(-\kappa r)}{r}\quad \mathrm{for}\quad r`__ platform can build the software -automatically in the cloud and skip this chapter. For more details on -running |es| in Gitpod, go to section :ref:`Running in the cloud`. +Users who only need a "default" installation of |es| and have a GitHub account +can build the software automatically in the cloud and directly go to section +:ref:`Using Codespaces`. Quickstart ---------- @@ -351,7 +350,9 @@ Installing requirements on macOS The first step is to install a C++ compiler, such as Xcode [10]_. Xcode is missing OpenMP, which is needed to enable shared-memory parallelization, -but the "R for macOS Developers" project provides binaries [11]_. +but binaries are available from Homebrew +(formula `libomp `__) +or from the "R for macOS Developers" project [11]_. To install libraries, a package manager will be needed. While our instructions below are specific to Homebrew, @@ -481,10 +482,6 @@ General features .. seealso:: :ref:`Electrostatics` -- ``MMM1D_MACHINE_PREC``: This enables high-precision Bessel functions - for MMM1D on CPU. Comes with a 60% slow-down penalty. The low-precision - functions are enabled by default and are precise enough for most applications. - - ``DIPOLES`` This activates the dipole-moment property of particles and switches on various magnetostatics algorithms @@ -658,7 +655,8 @@ using a CMake flag (see :ref:`Options and Variables`). :ref:`ScaFaCoS electrostatics`, :ref:`ScaFaCoS magnetostatics`). - ``GSL``: enables features relying on the GNU Scientific Library, e.g. - :meth:`espressomd.cluster_analysis.Cluster.fractal_dimension`. + :meth:`espressomd.cluster_analysis.Cluster.fractal_dimension` and + :class:`espressomd.electrostatics.MMM1D`. - ``STOKESIAN_DYNAMICS``: enable the Stokesian Dynamics propagator (see :ref:`Stokesian Dynamics`). Requires BLAS and LAPACK. diff --git a/doc/sphinx/magnetostatics.rst b/doc/sphinx/magnetostatics.rst index 7c5decd1cd7..74fdcaeba3d 100644 --- a/doc/sphinx/magnetostatics.rst +++ b/doc/sphinx/magnetostatics.rst @@ -54,6 +54,12 @@ Dipolar P3M :class:`espressomd.magnetostatics.DipolarP3M` +.. note:: + + Requires feature ``DIPOLES`` and + external feature ``FFTW``, enabled with + ``-D ESPRESSO_BUILD_WITH_FFTW=ON``. + This is the dipolar version of the P3M algorithm, described in :cite:`cerda08d`. Make sure that you know the relevance of the P3M parameters before using @@ -85,6 +91,10 @@ Dipolar Layer Correction (DLC) :class:`espressomd.magnetostatics.DLC` +.. note:: + + Requires feature ``DIPOLES``. + The dipolar layer correction (DLC) is used in conjunction with the dipolar P3M method to calculate dipolar interactions in a 2D-periodic system. It is based on :cite:`brodka04a` and the dipolar version of @@ -121,6 +131,12 @@ The method is used as follows:: Dipolar direct sum ------------------ +.. note:: + + Requires feature ``DIPOLES`` and optionally + external feature ``CUDA``, enabled with + ``-D ESPRESSO_BUILD_WITH_CUDA=ON``. + This interaction calculates energies and forces between dipoles by explicitly summing over all pairs. For the directions in which the system is periodic (as defined by ``system.periodicity``), it applies the @@ -175,10 +191,15 @@ ScaFaCoS magnetostatics :class:`espressomd.magnetostatics.Scafacos` +.. note:: + + Requires feature ``SCAFACOS_DIPOLES`` and + external feature ``SCAFACOS``, enabled with + ``-D ESPRESSO_BUILD_WITH_SCAFACOS=ON``. + |es| can use the methods from the ScaFaCoS *Scalable fast Coulomb solvers* -library for dipoles, if the methods support dipolar calculations. The feature -``SCAFACOS_DIPOLES`` has to be added to :file:`myconfig.hpp` to activate this -feature. Dipolar calculations are only included in the ``dipoles`` branch of +library for dipoles, if the methods support dipolar calculations. +Dipolar calculations are only included in the ``dipoles`` branch of the ScaFaCoS code. The specific methods available can be queried with :meth:`espressomd.electrostatics.Scafacos.get_available_methods`. diff --git a/doc/sphinx/running.rst b/doc/sphinx/running.rst index d75d7e15a83..337a83cbea0 100644 --- a/doc/sphinx/running.rst +++ b/doc/sphinx/running.rst @@ -182,47 +182,29 @@ path set to ``pypresso``. Running in the cloud ~~~~~~~~~~~~~~~~~~~~ -A `Gitpod `__ config file is provided to automatically -build |es| in its default configuration (`direct link -`__), which is -sufficient to run most tutorials. The Gitpod workspace can be accessed from -the `terminal via SSH `__ or from -a `web browser `__, -which uses the VS Code IDE. +.. _Using Binder: -To execute the tutorials, choose a Jupyter backend: +Using Binder +"""""""""""" -* VS Code Jupyter: navigate to ``ESPRESSO/build/doc/tutorials`` in the - project tree and open the notebook files; if the kernel drop-down menu - doesn't offer ``build/pypresso`` as a kernel, restart the VS Code IDE: - quit the workspace by closing the browser tab, re-open the tab and - click ``espressomd-espresso-...`` in the popup to restart the IDE - (don't click on the green button "New Workspace") +A pre-compiled version of |es| is available +for the `Binder platform `__ :cite:`jupyter18a` +(`direct link `__), +which uses the JupyterLab IDE. -* Jupyter Notebook: +.. _Using Codespaces: - .. code-block:: bash +Using Codespaces +"""""""""""""""" - cd ${GITPOD_REPO_ROOT}/build/doc/tutorials - ../../ipypresso notebook --NotebookApp.allow_origin="$(gp url 8888)" \ - --port=8888 --no-browser - -* JupyterLab: - - .. code-block:: bash - - cd ${GITPOD_REPO_ROOT}/build/doc/tutorials - ../../ipypresso lab --NotebookApp.allow_origin="$(gp url 8888)" \ - --port=8888 --no-browser - -For both Jupyter Notebook and JupyterLab, a notification will appear and say -that a new port 8888 has been made available. Click the orange "Make public" -button to open that port and then Ctrl+click one of the urls in the terminal -output to open the Jupyter backed in a pop-up window. - -To start a workspace from a specific branch, use a link in the following form: -``https://gitpod.io/#https://github.com/user_name/espresso/tree/branch_name``, -where ``user_name`` and ``branch_name`` need to be adapted. +A `Codespaces `__ config file is provided +to facilitate the building of |es| in the default configuration (`direct link +`__). +The codepace can be accessed from the terminal via the `GitHub CLI ssh command +`__ +or from a web browser (default), which uses the VS Code IDE. +Instructions to build |es| and execute the tutorials are available +in file :file:`.devcontainer/Readme.md`. .. _Parallel computing: diff --git a/doc/tutorials/CMakeLists.txt b/doc/tutorials/CMakeLists.txt index e151e1d0ba8..f43848f0045 100644 --- a/doc/tutorials/CMakeLists.txt +++ b/doc/tutorials/CMakeLists.txt @@ -121,3 +121,4 @@ add_subdirectory(widom_insertion) add_subdirectory(electrodes) add_subdirectory(grand_canonical_monte_carlo) add_subdirectory(mlip) +add_subdirectory(boltzmann_inversion) diff --git a/doc/tutorials/Readme.md b/doc/tutorials/Readme.md index b9d14e9ec6c..8da3ccbed69 100644 --- a/doc/tutorials/Readme.md +++ b/doc/tutorials/Readme.md @@ -8,7 +8,7 @@ physical systems. [comment]: # (Begin of tutorials landing page) [![Launch with Binder](https://mybinder.org/badge_logo.svg)](https://mybinder.org/v2/gh/jngrad/espresso-binder/HEAD) -[![Launch with Gitpod](https://img.shields.io/badge/launch-Gitpod-908a85?logo=gitpod)](https://gitpod.io/#https://github.com/espressomd/espresso) +[![Open in GitHub Codespaces](https://img.shields.io/badge/launch-Codespaces-579aca?logo=github)](https://codespaces.new/espressomd/espresso) ### Introductory tutorials diff --git a/doc/tutorials/boltzmann_inversion/CMakeLists.txt b/doc/tutorials/boltzmann_inversion/CMakeLists.txt new file mode 100644 index 00000000000..8369d3b4a8f --- /dev/null +++ b/doc/tutorials/boltzmann_inversion/CMakeLists.txt @@ -0,0 +1,21 @@ +# +# Copyright (C) 2025 The ESPResSo project +# +# This file is part of ESPResSo. +# +# ESPResSo is free software: you can redistribute it and/or modify +# it under the terms of the GNU General Public License as published by +# the Free Software Foundation, either version 3 of the License, or +# (at your option) any later version. +# +# ESPResSo 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 General Public License for more details. +# +# You should have received a copy of the GNU General Public License +# along with this program. If not, see . +# + +espresso_add_tutorial(TARGET tutorial_bi FILE "boltzmann_inversion.ipynb" + HTML_RUN) diff --git a/doc/tutorials/boltzmann_inversion/NotesForTutor.md b/doc/tutorials/boltzmann_inversion/NotesForTutor.md new file mode 100644 index 00000000000..378545b5401 --- /dev/null +++ b/doc/tutorials/boltzmann_inversion/NotesForTutor.md @@ -0,0 +1,21 @@ +# Boltzmann inversion + +## Learning objectives (physics) + +* Develop a simple coarse-grained model for charged colloidal particles with explicit counterions and salt using the method of Boltzmann inversion. +* In the coarse-grained model, only the colloidal particles are explicitly modeled, while the effects of the small ions are considered implicitly through an effective coarse-grained potential. +* Students should plot the RDF between colloidal particles in the system with explicit counterions and salt, interpret it physically and calculate the effective potential between colloidal particles from it. +* For consistency, the obtained effective potentials should be compared with analytical prediction by the Debye-Hückel theory, which is a linearized mean field approach. +* Finally, a coarse-grained MD run should be performed where only the colloidal particles interact with an effective coarse-grained potential. + +After the tutorial, students should be able to: + +* Explain the shorter runtime of a coarse-grained simulation compared to a simulation with explicit salt. +* Explain the discrepancy in the RDFs (especially in the short range) with explicit and implicit salt. +* Limitations of the Boltzmann inversion method, particularly with regards to many-body effects. + +## Learning objectives (ESPResSo) + +* Setting up a system with colloidal particles of desired size, counterions and salt. +* Calculating the RDF between colloidal particles using observables and accumulators. +* To use the effective potential via a tabulated interaction, while retaining the WCA potential for separations shorter than the tabulated cutoff. diff --git a/doc/tutorials/boltzmann_inversion/boltzmann_inversion.ipynb b/doc/tutorials/boltzmann_inversion/boltzmann_inversion.ipynb new file mode 100644 index 00000000000..ba0bdeeda69 --- /dev/null +++ b/doc/tutorials/boltzmann_inversion/boltzmann_inversion.ipynb @@ -0,0 +1,904 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "id": "0969c2de-b4f4-4b98-b417-bf3697f553e1", + "metadata": { + "slideshow": { + "slide_type": "" + }, + "tags": [] + }, + "source": [ + "# Systematic Coarse-Graining: Boltzmann Inversion" + ] + }, + { + "cell_type": "markdown", + "id": "3b28a690-e528-4315-90a2-1a1cd866ee2d", + "metadata": {}, + "source": [ + "## Introduction\n", + "\n", + "An important concept in statistical mechanics and molecular modeling is *coarse-graining*. \n", + "Coarse-graining is a procedure in which complexity of a system is reduced, leading to an effective description in terms of a smaller number of variables.\n", + "Due to the reduced number of degrees of freedom, such coarse-grained models can significantly reduce the time and computational resources required to simulate complex systems.\n", + "Therefore, coarse-grained simulations can also be used to simulate systems currently not accessible to atomistic simulation models.\n", + "Coarse-graining is a broadly applicable approach that encompasses many different methods, and consequently there exists a plethora of different ways to derive coarse-grained models.\n", + "In this tutorial, you will learn how to apply a systematic coarse-graining procedure, the Boltzmann inversion method.\n", + "\n", + "In the following, we consider a suspension of charged colloidal particles.\n", + "Colloidal particles (or colloids for short) are particles of sizes between 1 nanometer and 1 micrometer, which can be suspended in aqueous solution.\n", + "When colloidal particles carry an electrical charge, their interaction in solution is significantly modified by the large number of small salt ions, which dynamically rearrange in a charge cloud around the colloids that screens the bare Coulomb interactions.\n", + "However, including those small salt ions explicitly in a computer simulation of a colloidal suspension comes at a huge computational overhead. \n", + "In the spirit of coarse-grained modeling, in this tutorial we thus want to derive an effective pair potential between the colloids, which models the effect of the salt ions implicitly. \n", + "\n", + "To derive an effective pair potential between the colloidal particles, we first need to carry out simulations of a colloidal suspension including explicit salt ions. \n", + "In those simulations, we measure the radial distribution function (RDF) $g(r)$ between the colloids.\n", + "For sufficiently dilute systems, where many-body interactions are negligible, the RDF can be directly related to an effective pair potential $V^\\mathrm{eff}(r)$ between the colloids: $$g(r)\\propto \\exp\\left(-\\beta V^\\mathrm{eff}(r)\\right).$$\n", + "Here, we use a proportionality sign, because the effective pair potential is only determined up to an irrelevant additive constant, which does not contribute to the force.\n", + "The idea of Boltzmann inversion is to invert this relation, i.e. extract the effective pair potential as \n", + "$$V^\\mathrm{eff}(r) = -\\frac{1}{\\beta} \\ln\\left(g(r)\\right).$$\n", + "This effective interaction can then be used as a tabulated interaction in a simulation containing only colloids, providing an implicit description of the effects of salt." + ] + }, + { + "cell_type": "markdown", + "id": "1fa4e954-d582-4380-8ada-09f331055ef9", + "metadata": {}, + "source": [ + "## Simulation with Explicit Salt" + ] + }, + { + "cell_type": "markdown", + "id": "4f0f62f2-2013-47c6-8fbb-67438bdc7ab2", + "metadata": {}, + "source": [ + "### Simulation Setup\n", + "\n", + "To measure the RDF and extract the effective pair potential, we now set up the simulation model of a colloidal suspension including explicit salt ions in ESPResSo. \n", + "Let us begin by importing the required Python modules:" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "b3c718f8-2f70-47af-bdf9-729a5e5a657f", + "metadata": {}, + "outputs": [], + "source": [ + "import numpy as np\n", + "import matplotlib.pyplot as plt\n", + "import tqdm\n", + "import os\n", + "\n", + "import espressomd\n", + "from espressomd import electrostatics\n", + "import espressomd.observables\n", + "\n", + "required_features = [\"LENNARD_JONES\", \"ELECTROSTATICS\"]\n", + "espressomd.assert_features(required_features)" + ] + }, + { + "cell_type": "markdown", + "id": "8ad4c405-66f9-4a4a-9059-66f3adb714a7", + "metadata": {}, + "source": [ + "In the simulations, there are three kinds of particles: small anions, small cations and large (positively charged) colloidal particles.\n", + "All three particle types are represented as spherical particles that carry an electrical charge and additionally interact through a WCA potential.\n", + "For the colloidal particles, the WCA potential is shifted, corresponding to a larger radius of the particles.\n", + "To implement this model conveniently, we define dictionaries `types` and `charges`. \n", + "We set the charge of the colloidal particles to `+3`, the charge of the small cations to `+1` and the charge of the anions to `-1`.\n", + "All charges are measured in units of the elementary charge." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "a2df3e14-a676-4c40-af98-61b278c9d42f", + "metadata": {}, + "outputs": [], + "source": [ + "##### Particle types and charges\n", + "types = {\n", + " \"colloid\": 0,\n", + " \"cation\": 1,\n", + " \"anion\": 2,\n", + " }\n", + "\n", + "charges = {\n", + " \"colloid\": +3, \n", + " \"cation\": +1,\n", + " \"anion\": -1,\n", + " }" + ] + }, + { + "cell_type": "markdown", + "id": "f608be62-5b5f-4c7d-b7d0-9fbecce43191", + "metadata": {}, + "source": [ + "In the next step, we define our unit system. \n", + "Since ESPResSo works with reduced units, we need to establish a mapping between the simulation units and the usual SI units. \n", + "In the following, we take the unit of energy to be the thermal energy $k_{\\mathrm{B}}T$ at room temperature and the length scale $\\sigma = 3.55\\cdot 10^{-10}\\,\\mathrm{m}$, which roughly corresponds to the size of a solvated ion in aqueous solution.\n", + "Using these units, we can convert concentration between simulation units (number of particles / sigma**3) and SI units (mol/L) using the factor `PREF`.\n", + "We set the Bjerrum length to a value of `L_BJERRUM = 2.0`, which describes electrostatic interactions in water." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "80350c63-9f8c-4876-9ecb-87c014b7b39b", + "metadata": {}, + "outputs": [], + "source": [ + "##### Units\n", + "KT = 1.0\n", + "SIGMA = 3.55e-10 # Sigma in SI units\n", + "AVO = 6.022e+23 # Avogadro's number\n", + "PREF = 1/(10**3 * AVO * SIGMA**3) # Prefactor to mol/L\n", + "L_BJERRUM = 2.0" + ] + }, + { + "cell_type": "markdown", + "id": "7ba5da5b-bab9-43cb-a69b-2796ff5a2f7c", + "metadata": {}, + "source": [ + "Now, we define the colloid radius and number, as well the concentrations of colloids and salt. \n", + "The chosen parameters are a good compromise between a colloid density that is still low enough such that colloidal three-body interactions are negligible and a density that is too low to lead to a sufficient sampling of the RDF in a short simulation time.\n", + "Furthermore, the colloid charge and salt concentration are in regimes where the analytical Debye-Hückel theory is generally accurate, which means that the Debye-Hückel prediction serves as a useful benchmark, as we will see below." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "d4f73a73-279f-4631-ada5-293a29e1c12d", + "metadata": {}, + "outputs": [], + "source": [ + "RADIUS_COLLOID = 3.0\n", + "N_COLLOID = 50\n", + "\n", + "c_colloid = 0.001\n", + "C_COLLOID = c_colloid / PREF # colloid concentration in number of particles / sigma**3\n", + "\n", + "c_salt = 0.01\n", + "C_SALT = c_salt / PREF # salt concentration in number of particles / sigma**3" + ] + }, + { + "cell_type": "markdown", + "id": "644628a8-f77c-4244-80f7-4e296804b0af", + "metadata": {}, + "source": [ + "Using the chosen colloid number and concentrations, we can then calculate the box size as well as the number of salt ions and counterions.\n", + "Note that the number of counterions corresponds to the number of salt ion pairs plus the number of ions required to neutralize the charge on the colloidal particles:" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "53bb44fc-61e5-48f7-b135-84406d1ee8cd", + "metadata": {}, + "outputs": [], + "source": [ + "BOX_LENGTH = np.power(N_COLLOID / C_COLLOID, 1.0/3.0)\n", + "N_SALT = int(round(C_SALT * BOX_LENGTH ** 3))\n", + "N_COUNTERION = N_SALT + charges[\"colloid\"]*N_COLLOID" + ] + }, + { + "cell_type": "markdown", + "id": "da2e0879-2a87-4f5c-99d7-16a904048e04", + "metadata": {}, + "source": [ + "Finally, we set some parameters for the Langevin thermostat and the numerical integrator." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "e4afb082-dbb2-4cbf-a858-7fc2552fe65b", + "metadata": {}, + "outputs": [], + "source": [ + "##### Langevin-Thermostat and integrator\n", + "GAMMA = 1.0\n", + "DT = 0.01\n", + "SKIN = 1.2" + ] + }, + { + "cell_type": "markdown", + "id": "20656ae1-77ad-425a-a860-f50b1c157521", + "metadata": {}, + "source": [ + "We are now ready to create an instance of the ESPResSo system class and add the colloidal particles as well as the salt ions" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "5bae1a0e-ad37-441c-9624-073975bba1d8", + "metadata": {}, + "outputs": [], + "source": [ + "##### Create an instance of the ESPResSo system class\n", + "system = espressomd.System(box_l=[BOX_LENGTH]*3)\n", + "system.time_step = DT\n", + "system.cell_system.skin = SKIN\n", + "system.setup_type_map(type_list=list(types.values()))\n", + "\n", + "##### Add the colloidal particles\n", + "system.part.add(type=[types[\"colloid\"]]*N_COLLOID, pos=np.random.rand(N_COLLOID, 3) * BOX_LENGTH, q = [charges[\"colloid\"]]*N_COLLOID)\n", + "\n", + "##### Add the salt particles\n", + "if N_SALT>0:\n", + " system.part.add(type=[types[\"anion\"]]*N_COUNTERION, pos=np.random.rand(N_COUNTERION, 3) * BOX_LENGTH, q = [charges[\"anion\"]]*N_COUNTERION)\n", + " system.part.add(type=[types[\"cation\"]]*N_SALT, pos=np.random.rand(N_SALT, 3) * BOX_LENGTH, q = [charges[\"cation\"]]*N_SALT)" + ] + }, + { + "cell_type": "markdown", + "id": "e20854e2-ed5c-4601-adb6-574104c4404e", + "metadata": {}, + "source": [ + "In the next step, we add the electrostatics solver. \n", + "While the chosen accuracy of the P3M solver is rather low, it suffices for the purpose of this tutorial." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "0709d540-675d-4116-a7c7-896b395ea5d8", + "metadata": {}, + "outputs": [], + "source": [ + "p3m = electrostatics.P3M(prefactor = L_BJERRUM * KT, accuracy=1e-2)\n", + "system.electrostatics.solver = p3m" + ] + }, + { + "cell_type": "markdown", + "id": "37b4d53b-e653-4c63-a992-1d7322838e5c", + "metadata": {}, + "source": [ + "At last, we can set up the WCA interactions between all particles. \n", + "Due to the offset in the WCA interactions involving colloids, the simulation setup requires a bit of care to avoid divergences in the forces, because the particle positions are initalized randomly.\n", + "Therefore, we initially assign all particles the same short-range WCA interactions, corresponding to a particle diameter of 1.0 sigma:" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "12ba048f-96b5-430c-ac79-c7bb4e9641dc", + "metadata": {}, + "outputs": [], + "source": [ + "system.non_bonded_inter[types[\"anion\"], types[\"anion\"]].lennard_jones.set_params(epsilon=1.0, sigma=1.0, cutoff=(2)**(1.0/6.0), shift=\"auto\")\n", + "system.non_bonded_inter[types[\"cation\"], types[\"cation\"]].lennard_jones.set_params(epsilon=1.0, sigma=1.0, cutoff=(2)**(1.0/6.0), shift=\"auto\")\n", + "system.non_bonded_inter[types[\"anion\"], types[\"cation\"]].lennard_jones.set_params(epsilon=1.0, sigma=1.0, cutoff=(2)**(1.0/6.0), shift=\"auto\")\n", + "system.non_bonded_inter[types[\"anion\"], types[\"colloid\"]].lennard_jones.set_params(epsilon=1.0, sigma=1.0, cutoff=(2)**(1.0/6.0), shift=\"auto\")\n", + "system.non_bonded_inter[types[\"cation\"], types[\"colloid\"]].lennard_jones.set_params(epsilon=1.0, sigma=1.0, cutoff=(2)**(1.0/6.0), shift=\"auto\")\n", + "system.non_bonded_inter[types[\"colloid\"], types[\"colloid\"]].lennard_jones.set_params(epsilon=1.0, sigma=1.0, cutoff=(2)**(1.0/6.0), shift=\"auto\")" + ] + }, + { + "cell_type": "markdown", + "id": "7a188dc5-3bd8-4c87-a0f5-501e44979044", + "metadata": {}, + "source": [ + "### Warmup" + ] + }, + { + "cell_type": "markdown", + "id": "258819a3-d6d9-4b3b-9e98-806496623a75", + "metadata": {}, + "source": [ + "Before we can start the actual simulation, we have to set the colloid radius to the desired size.\n", + "To achieve this, we increase the particle diameter of the colloids in small steps and run a steepest descent integration after each increase to relax the system and avoid divergences.\n", + "Once the colloidal particles have the desired size, we run a final steepest descent integration." + ] + }, + { + "cell_type": "markdown", + "id": "941e5239-465c-4540-88f6-d3fcb459bc31", + "metadata": {}, + "source": [ + "**Exercise**\n", + "\n", + "* Set up the steepest descent integrator with parameters `f_max=0`, `gamma=1` and `max_displacement=0.01`.\n", + "* Write loop that slowly increases the colloid radius in steps of `delta_r = 0.1` until the final colloid size is reached. After each increase, run 100 integration steps to relax the system.\n", + "\n", + "**Hint**\n", + "\n", + "* You can change the parameters of the LJ interaction by invoking `lennard_jones.set_params` again.\n", + "* You need to change the offset for all LJ interactions involving the colloids. The final offset for colloid-colloid interactions should be `2*RADIUS_COLLOID` and that for colloid-ion interactions `RADIUS_COLLOID`." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "3610a4c5-f3c1-477a-b45f-6a6d57a138a9", + "metadata": {}, + "outputs": [], + "source": [ + "# SOLUTION CELL\n", + "##### Determine the number of steps needed\n", + "delta_r = 0.1\n", + "steps_needed = int(RADIUS_COLLOID/delta_r)\n", + "\n", + "##### CSet steepest descent integrator\n", + "system.integrator.set_steepest_descent(f_max=0, gamma=1, max_displacement=0.01)\n", + "system.integrator.run(100)\n", + "\n", + "##### Gradually increase colloid radius\n", + "current_radius = 0\n", + "for i in tqdm.tqdm(range(steps_needed)):\n", + " current_radius += delta_r\n", + " system.non_bonded_inter[types[\"anion\"], types[\"colloid\"]].lennard_jones.set_params(epsilon=1.0, sigma=1.0, cutoff=(2)**(1.0/6.0), shift=\"auto\", offset=current_radius)\n", + " system.non_bonded_inter[types[\"cation\"], types[\"colloid\"]].lennard_jones.set_params(epsilon=1.0, sigma=1.0, cutoff=(2)**(1.0/6.0), shift=\"auto\", offset=current_radius)\n", + " system.non_bonded_inter[types[\"colloid\"], types[\"colloid\"]].lennard_jones.set_params(epsilon=1.0, sigma=1.0, cutoff=(2)**(1.0/6.0), shift=\"auto\", offset=2*current_radius)\n", + " system.integrator.run(100)\n", + "\n", + "system.non_bonded_inter[types[\"anion\"], types[\"colloid\"]].lennard_jones.set_params(epsilon=1.0, sigma=1.0, cutoff=(2)**(1.0/6.0), shift=\"auto\", offset=RADIUS_COLLOID)\n", + "system.non_bonded_inter[types[\"cation\"], types[\"colloid\"]].lennard_jones.set_params(epsilon=1.0, sigma=1.0, cutoff=(2)**(1.0/6.0), shift=\"auto\", offset=RADIUS_COLLOID)\n", + "system.non_bonded_inter[types[\"colloid\"], types[\"colloid\"]].lennard_jones.set_params(epsilon=1.0, sigma=1.0, cutoff=(2)**(1.0/6.0), shift=\"auto\", offset=2*RADIUS_COLLOID)\n", + "system.integrator.run(100)" + ] + }, + { + "cell_type": "markdown", + "id": "6e0950ae-f51f-41d7-931a-22df2ff0d463", + "metadata": {}, + "source": [ + "Now, we turn on the Langevin thermostat and run a short warmup integration to equilibrate the system." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "85a30ada-bec8-4a41-90cf-cf603402a24f", + "metadata": {}, + "outputs": [], + "source": [ + "system.integrator.set_vv()\n", + "system.thermostat.set_langevin(seed=42, gamma=GAMMA, kT=KT)\n", + "\n", + "print(\"Warming up the system...\")\n", + "for i in tqdm.tqdm(range(10)):\n", + " system.integrator.run(1000)\n", + "print(\"Done.\")" + ] + }, + { + "cell_type": "markdown", + "id": "1033120c-885c-4149-a3f0-6c5ed5ae1e25", + "metadata": {}, + "source": [ + "### Production Run" + ] + }, + { + "cell_type": "markdown", + "id": "44575b31-fa2b-44eb-b25a-3f77207b4390", + "metadata": {}, + "source": [ + "Almost everything is ready now to start the production run.\n", + "As the last prepration step, we need to set up an observable to measure the RDF. \n", + "The chosen bin size for the RDF is rather large at 0.5 sigma, but allows for a good estimate of the effective interaction with a comparatively small number of samples.\n", + "This choice is valid, because except for the short-ranged, strongly repulsive WCA part, the effective potential varies on a much longer length scale than the bin size." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "0c91bc6b-5b57-4df5-8b30-181a674d7b5b", + "metadata": {}, + "outputs": [], + "source": [ + "# Parameters for the radial distribution function\n", + "BIN_WIDTH = 0.5\n", + "R_MIN = 0.0\n", + "R_MAX = system.box_l[0] / 2.0\n", + "N_BINS = int((R_MAX-R_MIN)/BIN_WIDTH)" + ] + }, + { + "cell_type": "markdown", + "id": "9551fb6f-eb67-4225-81ca-98025d601520", + "metadata": {}, + "source": [ + "**Exercise**\n", + "* Instantiate a [RDF](https://espressomd.github.io/doc/espressomd.html#espressomd.observables.RDF) observable\n", + "* Instantiate a [MeanVarianceCalculator](https://espressomd.github.io/doc/analysis.html#mean-variance-calculator) accumulator to track the RDF over time. Samples should be taken every ``steps_per_subsample`` steps.\n", + "* Add the accumulator to the [auto_update_accumulators](https://espressomd.github.io/doc/espressomd.html#espressomd.accumulators.AutoUpdateAccumulators) of the system for automatic updates" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "549819d4-b55b-4f87-bf5d-f1487950a50c", + "metadata": { + "slideshow": { + "slide_type": "" + }, + "tags": [] + }, + "outputs": [], + "source": [ + "# SOLUTION CELL\n", + "rdf_obs = espressomd.observables.RDF(ids1=system.part.select(type=types[\"colloid\"]).id, min_r=R_MIN, max_r=R_MAX, n_r_bins=N_BINS)\n", + "rdf_acc = espressomd.accumulators.MeanVarianceCalculator(obs=rdf_obs, delta_N=100)\n", + "system.auto_update_accumulators.add(rdf_acc)" + ] + }, + { + "cell_type": "markdown", + "id": "3151a483-bce2-4636-ab39-efbd3c0e9974", + "metadata": {}, + "source": [ + "Now, we are ready to perform our production run in order to measure the RDF:" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "cc0d438f-e196-47e3-bb0c-53645097a495", + "metadata": {}, + "outputs": [], + "source": [ + "print(\"Performing the production run...\")\n", + "for i in tqdm.tqdm(range(2000)):\n", + " system.integrator.run(100)\n", + "print(\"Done.\")" + ] + }, + { + "cell_type": "markdown", + "id": "d006f57f-9c5f-4ae6-bc47-d7b6e27e970e", + "metadata": {}, + "source": [ + "## Data Analysis and Boltzmann Inversion" + ] + }, + { + "cell_type": "markdown", + "id": "683a9beb-e632-459f-bb9c-4d0f284f49a8", + "metadata": {}, + "source": [ + "We can now analyze our simulation results. \n", + "In the first step, we analyze the RDF obtained from the simulation." + ] + }, + { + "cell_type": "markdown", + "id": "d7b727c0-ca1f-457e-8c3b-fdad445ceed6", + "metadata": {}, + "source": [ + "**Exercise**\n", + "* Get the mean RDF (define ``rdf``) from the accumulator and the histogram bin centers (define ``rs``) from the observable.\n", + "* Plot the RDF. Interpret your result physically." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "cd68724d-b339-41d3-a724-47512755569e", + "metadata": {}, + "outputs": [], + "source": [ + "# SOLUTION CELL\n", + "rdf = rdf_acc.mean()\n", + "rs = rdf_obs.bin_centers()\n", + "\n", + "plt.plot(rs, rdf)\n", + "plt.xlabel(r'distance $r$/$\\sigma$')\n", + "plt.ylabel('colloid-colloid RDF $g(r)$')\n", + "plt.show()" + ] + }, + { + "cell_type": "markdown", + "id": "7ab797a6-ef3d-4f6a-84d2-14bf7edf34e5", + "metadata": {}, + "source": [ + "Using the RDF, we can now calculate a numerical approximation of the effective pair potential between colloids." + ] + }, + { + "cell_type": "markdown", + "id": "7f729275-d538-4998-b94d-c1125e48c678", + "metadata": {}, + "source": [ + "**Exercise**\n", + "* Use the Boltzmann inversion method, i.e. the equation $V^\\mathrm{eff}(r) = -\\frac{1}{\\beta} \\ln\\left(g(r)\\right)$ to calculate the effective pair potential using the RDF.\n", + "* Create a plot of the effective pair potential and interpret it." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "9421e0b7-9b11-4784-bbf3-5607384b14c6", + "metadata": {}, + "outputs": [], + "source": [ + "# SOLUTION CELL\n", + "veff = -np.log(rdf)\n", + "\n", + "plt.plot(rs, veff)\n", + "plt.xlabel(r\"distance $r$/$\\sigma$\")\n", + "plt.ylabel(\"effective pair potential $V_{eff}(r)$/$kT$\")\n", + "plt.ylim((-0.3,2.0))\n", + "plt.show()" + ] + }, + { + "cell_type": "markdown", + "id": "e8c3ad6a-76b3-4b28-bc09-3a29d28e5ae5", + "metadata": {}, + "source": [ + "For charged colloids with small electrostatic surface potentials, the effective pair potential can be calculated analytically using the Debye-Hückel theory, which was first applied to colloids by Derjaguin, Landau, Verwey and Overbeek, resulting in the famous [DLVO theory](https://en.wikipedia.org/wiki/DLVO_theory).\n", + "The theory predicts that the effective potential between two colloids of charge $Z$ and radius $R$ is given by\n", + "$$\\beta V^{\\mathrm{eff}}(r) = Z^2 \\lambda_\\text{B} \\, \\left(\\frac{e^{\\kappa R}}{1 + \\kappa R}\\right)^2 \\, \\frac{e^{-\\kappa r}}{r}.$$\n", + "Here, $\\kappa$ is the salt concentration-dependent Debye screening constant and $\\lambda_\\text{B}$ is the Bjerrum length.\n", + "The inverse $\\kappa^{-1}$ of the Debye screening constant is called the Debye length and characterizes the range of the effective interaction.\n", + "For our system, the analytical result for the effective potential is a useful benchmark to test the simulation results against." + ] + }, + { + "cell_type": "markdown", + "id": "fc1f76ae-2523-4b08-b639-23f2277b6a98", + "metadata": {}, + "source": [ + "**Exercise**\n", + "* Use the function `veff_debye_hueckel` defined below to plot a comparison of the Debye-Hückel potential and the effective pair potential obtained from the simulations." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "0bc2b1ae-7b82-45a3-b348-6fbab9811039", + "metadata": {}, + "outputs": [], + "source": [ + "def veff_debye_hueckel(r, q, c_salt, c_colloid, radius_colloid):\n", + " \"\"\"\n", + " Calculate the effective potential between the colloids according to the Debye-Hückel theory.\n", + " The calculation takes into account both the salt and the additional counterions.\n", + "\n", + " Args:\n", + " r (float): The distance between the colloids in sigma.\n", + " q (float): The charge of the colloids in elementary charge units.\n", + " c_salt (float): The concentration of salt in the solution in M (mol/L).\n", + " c_colloid (float): The concentration of colloids in the solution in M (mol/L).\n", + " radius_colloid (float): The radius of the colloids in sigma.\n", + "\n", + " Returns:\n", + " float: The effective potential between the colloids in kT.\n", + " \"\"\"\n", + " bjerrum_length = 2.0\n", + " lambda_D = (0.304 / (np.sqrt(c_salt) * np.sqrt(1 + q * c_colloid/ (2 * c_salt)))) / (0.355) # Calculate the Debye screening length\n", + " kappa = 1/lambda_D # Calculate the Debye screening parameter\n", + "\n", + " if r > 2*(radius_colloid+0.5):\n", + " return q**2 * bjerrum_length * (np.exp(kappa * radius_colloid) /(1 + kappa*radius_colloid))**2 * np.exp(-kappa * r) / r\n", + " else:\n", + " return np.nan" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "f8c8e795-6f2c-42d6-b3e8-c7935b9d368c", + "metadata": {}, + "outputs": [], + "source": [ + "# SOLUTION CELL\n", + "plt.plot(rs, veff, label=\"Boltzmann Inversion\")\n", + "plt.plot(rs, [veff_debye_hueckel(r, charges[\"colloid\"], c_salt, c_colloid, RADIUS_COLLOID) for r in rs], label=\"Debye-Hückel\")\n", + "plt.legend()\n", + "plt.xlabel(r\"distance $r$/$\\sigma$\")\n", + "plt.ylabel(\"effective pair potential $V_{eff}(r)$/$kT$\")\n", + "plt.ylim((-0.3,2.0))\n", + "plt.show()" + ] + }, + { + "cell_type": "markdown", + "id": "3825c4a2-cc21-4311-913f-ff827f60720e", + "metadata": {}, + "source": [ + "## Simulation using the Developed Coarse-Grained Potential" + ] + }, + { + "cell_type": "markdown", + "id": "f6666389-1011-4ff3-abc4-4de9c4c669c5", + "metadata": {}, + "source": [ + "### Setup" + ] + }, + { + "cell_type": "markdown", + "id": "2d253fd9-7dd9-4bae-abea-af0228a6eaa0", + "metadata": {}, + "source": [ + "Having obtained the effective potential between the colloidal particles, we can now perform a coarse-grained simulation using this potential. \n", + "In a first step, we remove all particles from the system and turn off the electrostatics solver and the thermostat:" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "b7836da6-2c3e-49fe-82d4-ed7bd785acef", + "metadata": {}, + "outputs": [], + "source": [ + "system.part.clear()\n", + "system.electrostatics.clear()\n", + "system.thermostat.turn_off()" + ] + }, + { + "cell_type": "markdown", + "id": "e46bdb24-531b-4cc9-b314-c248a4c48b02", + "metadata": {}, + "source": [ + "Now, we add back the relevant particles to the system. \n", + "In contrast to the simulations we carried out before, we only need to add the colloidal particles, because there are no explicit salt ions present in the simulations.\n", + "The effect of salt is taken into account implicitly through the use of the effective potential that we have derived using Boltzmann inversion." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "01dcb94e-285a-4a7a-8c9c-f3afd82b9591", + "metadata": {}, + "outputs": [], + "source": [ + "##### Add the colloidal particles\n", + "system.part.add(type=[types[\"colloid\"]]*N_COLLOID, pos=np.random.rand(N_COLLOID, 3) * BOX_LENGTH)" + ] + }, + { + "cell_type": "markdown", + "id": "70c5b89c-3d53-4b35-a61f-9af363433f31", + "metadata": {}, + "source": [ + "In principle, we would now have to add only the effective potential. \n", + "However, because we ran only very short, unconstrained MD simulations, the strongly repulsive part of the effective potential is not properly sampled.\n", + "This can lead to instabilites that make the simulation crash.\n", + "As a solution, we will cut off the tabulated interaction for small distances and still retain the WCA interaction to stabilize the system.\n", + "Note that we set the initial particle diameter to `1.0` again in order to avoid diverging forces:" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "cfa7a4d6-961a-48a5-abd1-93e8f295c132", + "metadata": {}, + "outputs": [], + "source": [ + "##### Add excluded volume interactions\n", + "system.non_bonded_inter[types[\"colloid\"], types[\"colloid\"]].lennard_jones.set_params(epsilon=1.0, sigma=1.0, cutoff=(2)**(1.0/6.0), shift=\"auto\")" + ] + }, + { + "cell_type": "markdown", + "id": "7dbce44f-91d7-4da2-ae5b-849ec1e0a05d", + "metadata": {}, + "source": [ + "The only part that is missing now is the effective potential we calculated in the previous section. \n", + "To avoid double-counting the repulsive WCA potential we just added, we account for the effective potential only at distances larger than `2*RADIUS_COLLOID+2**(1/6)`, i.e. the cutoff of the WCA interaction.\n", + "Furthermore, as plot we generated above shows, the interaction rapidly decays for large distances. \n", + "Accordingly, we can also introduce an upper cutoff." + ] + }, + { + "cell_type": "markdown", + "id": "6569c310-bd10-4a14-b8d4-526962067345", + "metadata": {}, + "source": [ + "**Exercise**\n", + "* Calculate the effective force from the effective potential using a numerical derivative.\n", + "* Examine your plot of the effective potential and decide on an upper cutoff for the tabulated interaction. Why did you choose this value?\n", + "* Use the effective potential and the effective force to set up a tabulated interaction in ESPResSo. The tabulated interaction should be defined only between a minimum distance of `2*RADIUS_COLLOID+2**(1/6)` (the cutoff of the WCA interaction) and the upper cutoff you determined in the previous step.\n", + "\n", + "**Hint**\n", + "\n", + "* Use the NumPy function `np.gradient` to numerically calculate the effective force from the effective potential.\n", + "* Use NumPy array filtering to select the relevant values of the effective potential and effective force for distances that lie between the minimum distance and the upper cutoff.\n", + "* Refer to the the [documentation](https://espressomd.github.io/doc/inter_non-bonded.html#tabulated-interaction) to learn how tabulated interactions can be set up in ESPResSo." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "5ceafda5-1c73-4507-9105-0b58e96908cc", + "metadata": {}, + "outputs": [], + "source": [ + "# SOLUTION CELL\n", + "##### Calculate the effective force\n", + "effective_force = -np.gradient(veff, rs)\n", + "\n", + "##### Set the range of interaction\n", + "cutoff_tabulated = 25\n", + "veff = veff[(2*RADIUS_COLLOID+2**(1/6). # + import os import sys import time @@ -51,7 +52,8 @@ def minimize(system, energy_target): exit(1) -def get_timings(system, n_steps, n_iterations, verbose=True): +def get_timings(system, n_steps, n_iterations, verbose=True, + retune_skin_after_steps=None): ''' Time the integration loop and write the state of the system to stdout. @@ -65,6 +67,9 @@ def get_timings(system, n_steps, n_iterations, verbose=True): Number of timings. verbose: :obj:`bool` Whether to print the state of the system during timing. + retune_skin_after_steps: :obj:`int`, optional + If provided, retune the skin every this many iterations to within 10% + of the current skin value. Returns ------- @@ -76,6 +81,21 @@ def get_timings(system, n_steps, n_iterations, verbose=True): print(f"Timing every {n_steps} steps") timings = [] for i in range(n_iterations): + # Retune skin if requested + if retune_skin_after_steps is not None and i % retune_skin_after_steps == 0: + current_skin = system.cell_system.skin + min_skin_retune = current_skin / 1.1 + max_skin_retune = current_skin * 1.1 + if verbose: + print(f"Retuning skin at iteration {i} " + f"(current: {current_skin:.3f}, " + f"range: [{min_skin_retune:.3f}, " + f"{max_skin_retune:.3f}])") + new_skin = system.cell_system.tune_skin( + min_skin=min_skin_retune, max_skin=max_skin_retune, tol=current_skin * 0.0125, int_steps=n_steps // 4) + if verbose: + print(f"New skin: {new_skin:.3f}") + tick = time.time() system.integrator.run(n_steps) tock = time.time() @@ -85,7 +105,7 @@ def get_timings(system, n_steps, n_iterations, verbose=True): energy = system.analysis.energy()["total"] verlet = system.cell_system.get_state()["verlet_reuse"] print( - f"step {i}, time: {1000 * t:.2f} ms, verlet: {verlet:.2f}, energy: {energy:.2e}") + f"step {i}, time: {1000 * t:.4f} ms, verlet: {verlet:.2f}, energy: {energy:.2e}") return np.array(timings) diff --git a/maintainer/benchmarks/lj.py b/maintainer/benchmarks/lj.py index 4ddd4e229b8..9bdeb825a0c 100644 --- a/maintainer/benchmarks/lj.py +++ b/maintainer/benchmarks/lj.py @@ -86,7 +86,7 @@ # Integration parameters ############################################################# system.time_step = 0.01 -system.cell_system.skin = 0.5 +system.cell_system.skin = 0.4 # Interaction setup ############################################################# @@ -121,12 +121,10 @@ # tuning and equilibration min_skin = 0.2 max_skin = 1.0 -print("Tune skin: {:.3f}".format(system.cell_system.tune_skin( - min_skin=min_skin, max_skin=max_skin, tol=0.05, int_steps=100))) print("Equilibration") system.integrator.run(min(5 * measurement_steps, 60000)) print("Tune skin: {:.3f}".format(system.cell_system.tune_skin( - min_skin=min_skin, max_skin=max_skin, tol=0.05, int_steps=100))) + min_skin=min_skin, max_skin=max_skin, tol=0.025, int_steps=200))) print("Equilibration") system.integrator.run(min(10 * measurement_steps, 60000)) @@ -150,7 +148,8 @@ def main_thread(): # time integration loop -timings = benchmarks.get_timings(system, measurement_steps, n_iterations) +timings = benchmarks.get_timings( + system, measurement_steps, n_iterations, retune_skin_after_steps=5) # average time avg, ci = benchmarks.get_average_time(timings) diff --git a/maintainer/benchmarks/particle_access.py b/maintainer/benchmarks/particle_access.py new file mode 100644 index 00000000000..0e0dfc04364 --- /dev/null +++ b/maintainer/benchmarks/particle_access.py @@ -0,0 +1,277 @@ +# +# Copyright (C) 2025 The ESPResSo project +# +# This file is part of ESPResSo. +# +# ESPResSo is free software: you can redistribute it and/or modify +# it under the terms of the GNU General Public License as published by +# the Free Software Foundation, either version 3 of the License, or +# (at your option) any later version. +# +# ESPResSo 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 General Public License for more details. +# +# You should have received a copy of the GNU General Public License +# along with this program. If not, see . +# + +import espressomd +import numpy as np +import argparse +import time + + +def benchmark_operation(operation_func, n_iter): + """Execute an operation multiple times and return timing statistics.""" + timings = [] + for _ in range(n_iter): + tick = time.time() + operation_func() + tock = time.time() + timings.append(tock - tick) + return np.mean(timings), np.std(timings) + + +def print_section_header(title): + """Print a formatted section header.""" + print("=" * 60) + print(title) + print("=" * 60) + + +def filter_available_properties(properties, particle): + """Filter properties to only those available for the given particle.""" + return [prop for prop in properties if hasattr(particle, prop)] + + +parser = argparse.ArgumentParser( + description="Benchmark particle creation and property access." +) +parser.add_argument( + "--particles_per_core", + metavar="N", + action="store", + type=int, + default=1000, + required=False, + help="Number of particles per core (default: 1000)", +) +parser.add_argument( + "--additional_properties", + metavar="PROPS", + action="store", + type=str, + default="", + required=False, + help="Comma-separated list of additional properties to test (e.g., 'mass,type,v')", +) +parser.add_argument( + "--n_iter", + metavar="N", + action="store", + type=int, + default=20, + required=False, + help="Number of iterations for timing (default: 10)", +) + +args = parser.parse_args() + +# Parse additional properties +additional_props = [] +if args.additional_properties: + additional_props = [p.strip() + for p in args.additional_properties.split(",")] + +# System setup +system = espressomd.System(box_l=[10.0, 10.0, 10.0]) +system.time_step = 0.01 +system.cell_system.skin = 0.5 + +n_proc = system.cell_system.get_state()["n_nodes"] +n_part = n_proc * args.particles_per_core + +print(f"Benchmark Configuration:") +print(f" Particles per core: {args.particles_per_core}") +print(f" Total particles: {n_part}") +print(f" MPI ranks: {n_proc}") +print(f" Iterations: {args.n_iter}") +print( + f" Additional properties: { + additional_props if additional_props else 'None'}" +) +print() + +# Test 1: Bulk particle creation without providing IDs +print_section_header("Test 1: Bulk particle creation without IDs") + + +def create_particles_no_ids(): + system.part.clear() + positions = np.random.random((n_part, 3)) * system.box_l + system.part.add(pos=positions) + + +avg, std = benchmark_operation(create_particles_no_ids, args.n_iter) +print(f"Average: {avg:.6f} s ± {std:.6f} s ({n_part / avg:.0f} particles/s)") +print() + +# Test 2: Bulk particle creation with linear IDs +print_section_header("Test 2: Bulk particle creation with linear IDs") + + +def create_particles_with_ids(): + system.part.clear() + positions = np.random.random((n_part, 3)) * system.box_l + ids = np.arange(n_part, dtype=int) + system.part.add(pos=positions, id=ids) + + +avg, std = benchmark_operation(create_particles_with_ids, args.n_iter) +print(f"Average: {avg:.6f} s ± {std:.6f} s ({n_part / avg:.0f} particles/s)") +print() + +# Setup particles for property access tests +system.part.clear() +positions = np.random.random((n_part, 3)) * system.box_l +system.part.add(pos=positions) + +# Test 3: Single particle property access +print_section_header("Test 3: Single particle property access") + +# Base properties to test +base_properties_read = ["pos", "pos_folded", "q", "f"] +base_properties_write = ["pos", "q"] + +# Add additional properties and filter by availability +test_particle = system.part.by_id(0) +all_props_read = filter_available_properties( + base_properties_read + additional_props, test_particle +) +all_props_write = filter_available_properties( + base_properties_write + + [p for p in additional_props if p not in ["f", "pos_folded"]], + test_particle, +) +test_values = {prop: getattr(test_particle, prop) for prop in all_props_write} +test_values_slice = { + prop: np.array([getattr(test_particle, prop)] * n_part) for prop in all_props_write +} + + +def benchmark_property_access(properties, access_func, time_scale=1e6, time_unit="µs"): + """Benchmark property access (read or write) and print results.""" + for prop in properties: + avg, std = benchmark_operation(lambda: access_func(prop), args.n_iter) + print( + f" {prop}: {avg * time_scale:.3f} {time_unit} ± {std * + time_scale:.3f} {time_unit}, {1 / avg:.0f} particles/s" + ) + + +# Test single particle read +print("\nSingle particle property READ:") +benchmark_property_access( + all_props_read, lambda prop: getattr(test_particle, prop)) + +# Test single particle write +print("\nSingle particle property WRITE:") + + +def write_property(prop): + setattr(test_particle, prop, test_values[prop]) + + +benchmark_property_access(all_props_write, write_property) +print() + +# Test 4: Slice property access (re-instantiating slice) +print_section_header("Test 4: Slice property access (re-instantiating slice)") + + +def benchmark_slice_access( + properties, access_func, description, time_scale=1e3, time_unit="ms" +): + """Benchmark slice property access and print results.""" + print(f"\n{description}:") + for prop in properties: + avg, std = benchmark_operation(lambda: access_func(prop), args.n_iter) + print( + f" {prop}: {avg * time_scale:.3f} {time_unit} ± {std * + time_scale:.3f} {time_unit} ({n_part / avg:.0f} particles/s)" + ) + + +# Test slice read - re-instantiating slice on every access +benchmark_slice_access( + all_props_read, + lambda prop: getattr(system.part.all(), prop), + "Slice property READ (re-instantiating slice)", +) + +# Test slice write - re-instantiating slice on every access + + +def write_slice_property(prop): + setattr(system.part.all(), prop, test_values_slice[prop]) + + +benchmark_slice_access( + all_props_write, + write_slice_property, + "Slice property WRITE (re-instantiating slice)", +) +print() + +# Test 5: Slice property access (cached slice) +print_section_header("Test 5: Slice property access (cached slice)") + + +def benchmark_cached_slice_access(properties, access_func, warmup_func, description): + """Benchmark cached slice property access with warmup.""" + print(f"\n{description}:") + for prop in properties: + timings = [] + for _ in range(args.n_iter): + particle_slice = system.part.all() + # Warm-up access + warmup_func(particle_slice, prop) + # Timed access + tick = time.time() + access_func(particle_slice, prop) + tock = time.time() + timings.append(tock - tick) + avg = np.mean(timings) + std = np.std(timings) + print( + f" {prop}: {avg * 1e3:.3f} ms ± {std * + 1e3:.3f} ms ({n_part / avg:.0f} particles/s)" + ) + + +# Test slice read - cached slice +benchmark_cached_slice_access( + all_props_read, + lambda slice_obj, prop: getattr(slice_obj, prop), + lambda slice_obj, prop: getattr(slice_obj, prop), + "Slice property READ (cached slice)", +) + +# Test slice write - cached slice + + +def write_cached_slice(slice_obj, prop): + setattr(slice_obj, prop, test_values_slice[prop]) + + +benchmark_cached_slice_access( + all_props_write, + write_cached_slice, + write_cached_slice, + "Slice property WRITE (cached slice)", +) + +print() +print_section_header("Benchmark completed") diff --git a/maintainer/configs/no_rotation.hpp b/maintainer/configs/no_rotation.hpp index d1a6f418c35..eff6a247671 100644 --- a/maintainer/configs/no_rotation.hpp +++ b/maintainer/configs/no_rotation.hpp @@ -18,13 +18,6 @@ * You should have received a copy of the GNU General Public License * along with this program. If not, see . */ -/* This is the default myconfig.hpp-file. If no other myconfig-file is - found, this file is used. - - DO NOT MODIFY THIS FILE! It should be modified *only* by the - maintainers of ESPResSo, as it has a profound impact on many users, - in particular newbies. -*/ // Geometry, equation of motion, thermostat/barostat #define MASS @@ -33,15 +26,7 @@ #define BOND_CONSTRAINT #define NPT #define DPD - -// Charges and dipoles #define ELECTROSTATICS -#define MMM1D_MACHINE_PREC -#ifdef CUDA -#define MMM1D_GPU -#endif - -// Force/energy calculation #define EXCLUSIONS #define TABULATED diff --git a/maintainer/files_with_header.sh b/maintainer/files_with_header.sh index 0508710610c..5135e5bbcb6 100755 --- a/maintainer/files_with_header.sh +++ b/maintainer/files_with_header.sh @@ -25,7 +25,7 @@ git ls-files --exclude-standard | grep -vE '\.(blk|gz|npz|data|dat|tab|chk|jpg|png|pdf|fig|gif|xcf|css|bib|vtf|vtk|svg|ico|eps|rst|ipynb)$' | grep -vE '(^testsuite/python/data/|^old/|^libs/|^samples/immersed_boundary/tables/)' | grep -vE '(^ChangeLog|^AUTHORS|^COPYING|^NEWS|^INSTALL|^Readme\.md|^CONTRIBUTING\.md)' | -grep -vE '(^\.gitmodules|^\.github|^\.gitignore|^\.codecov\.yml|^\.gitlab-ci\.yml|^\.kodiak\.toml)' | +grep -vE '(^\.gitmodules|^\.github|^\.devcontainer|^\.vscode|^\.gitignore|^\.codecov\.yml|^\.gitlab-ci\.yml|^\.kodiak\.toml)' | grep -vE '(^\.clang-format|^\.cmake-format|^\.clang-tidy|^\.coveragerc|^\.pylintrc|ubsan\.supp)' | grep -vE '(^\.lgtm\.yml|^\.gitpod\.yml|^\.cppcheck|^\.pre-commit-config\.yaml|^requirements\.txt)' | grep -vE '(Doxyfile|latexmk\.1|latexmkrc|^doc/tutorials/.*/NotesForTutor\.md)' diff --git a/src/config/features.def b/src/config/features.def index 5969e87d231..83fdee18f0a 100644 --- a/src/config/features.def +++ b/src/config/features.def @@ -54,7 +54,6 @@ ROTATIONAL_INERTIA implies ROTATION /* Electrostatics */ ELECTROSTATICS P3M equals ELECTROSTATICS and FFTW -MMM1D_MACHINE_PREC requires ELECTROSTATICS WALBERLA_FFT equals WALBERLA and FFTW /* Magnetostatics */ diff --git a/src/core/CMakeLists.txt b/src/core/CMakeLists.txt index 6d616bee47a..d46f6920a48 100644 --- a/src/core/CMakeLists.txt +++ b/src/core/CMakeLists.txt @@ -81,6 +81,7 @@ target_link_libraries( $<$:Kokkos::kokkos> $<$:Cabana::Core> $<$:espresso::walberla> + $<$:GSL::gsl> PUBLIC espresso::utils MPI::MPI_CXX Random123 espresso::particle_observables Boost::serialization Boost::mpi espresso::instrumentation) diff --git a/src/core/collision_detection/BindAtPointOfCollision.cpp b/src/core/collision_detection/BindAtPointOfCollision.cpp index f1a16fa7a48..8e7cc8f09c8 100644 --- a/src/core/collision_detection/BindAtPointOfCollision.cpp +++ b/src/core/collision_detection/BindAtPointOfCollision.cpp @@ -132,7 +132,12 @@ void BindAtPointOfCollision::handle_collisions( auto const pos1 = p1->pos() - vec21 * vs_placement; auto const pos2 = p1->pos() - vec21 * (1. - vs_placement); - auto handle_particle = [&](Particle *p, Utils::Vector3d const &pos) { + auto handle_particle = [& +#if defined(__clang__) and defined(__cray__) + , + pid1 = pid1, pid2 = pid2 +#endif + ](Particle *p, Utils::Vector3d const &pos) { if (not p->is_ghost()) { place_vs_and_relate_to_particle(cell_structure, box_geo, part_type_vs, min_global_cut, current_vs_pid, pos, diff --git a/src/core/communication.cpp b/src/core/communication.cpp index d7ec6cf07f2..0fc8146dafd 100644 --- a/src/core/communication.cpp +++ b/src/core/communication.cpp @@ -117,10 +117,10 @@ CommunicationEnvironment::CommunicationEnvironment( m_is_mpi_gpu_aware |= (mpich_gpu_env and *mpich_gpu_env == "1"); #endif // defined(MPICH) -#if defined(_CRAYC) +#if defined(_CRAYC) or defined(__cray__) auto const cray_mpich_gpu_env = get_env_variable("MPICH_GPU_SUPPORT_ENABLED"); m_is_mpi_gpu_aware |= (cray_mpich_gpu_env and *cray_mpich_gpu_env == "1"); -#endif // defined(_CRAYC) +#endif // defined(_CRAYC) or defined(__cray__) communicator.full_initialization(); diff --git a/src/core/electrostatics/CMakeLists.txt b/src/core/electrostatics/CMakeLists.txt index 1432fc0b14b..63f6dc19d86 100644 --- a/src/core/electrostatics/CMakeLists.txt +++ b/src/core/electrostatics/CMakeLists.txt @@ -23,8 +23,7 @@ target_sources( ${CMAKE_CURRENT_SOURCE_DIR}/icc.cpp ${CMAKE_CURRENT_SOURCE_DIR}/mmm1d.cpp ${CMAKE_CURRENT_SOURCE_DIR}/mmm-modpsi.cpp - ${CMAKE_CURRENT_SOURCE_DIR}/scafacos_impl.cpp - ${CMAKE_CURRENT_SOURCE_DIR}/specfunc.cpp) + ${CMAKE_CURRENT_SOURCE_DIR}/scafacos_impl.cpp) if(ESPRESSO_BUILD_WITH_FFTW) target_sources(espresso_p3m PRIVATE p3m.cpp) diff --git a/src/core/electrostatics/coulomb.cpp b/src/core/electrostatics/coulomb.cpp index 8ba5725595c..0751e276486 100644 --- a/src/core/electrostatics/coulomb.cpp +++ b/src/core/electrostatics/coulomb.cpp @@ -152,9 +152,11 @@ struct ShortRangeCutoff { std::visit(*this, actor->base_solver)); } #endif // ESPRESSO_P3M +#ifdef ESPRESSO_GSL auto operator()(std::shared_ptr const &) const { return std::numeric_limits::infinity(); } +#endif #ifdef ESPRESSO_SCAFACOS auto operator()(std::shared_ptr const &actor) const { return actor->get_r_cut(); @@ -217,7 +219,9 @@ struct LongRangeForce { } #endif /* Several algorithms only provide near-field kernels */ +#ifdef ESPRESSO_GSL void operator()(std::shared_ptr const &) const {} +#endif void operator()(std::shared_ptr const &) const {} void operator()(std::shared_ptr const &) const {} @@ -244,7 +248,9 @@ struct LongRangeEnergy { } #endif /* Several algorithms only provide near-field kernels */ +#ifdef ESPRESSO_GSL auto operator()(std::shared_ptr const &) const { return 0.; } +#endif auto operator()(std::shared_ptr const &) const { return 0.; } auto operator()(std::shared_ptr const &) const { return 0.; } diff --git a/src/core/electrostatics/coulomb.hpp b/src/core/electrostatics/coulomb.hpp index 66c406c5969..13301b339a1 100644 --- a/src/core/electrostatics/coulomb.hpp +++ b/src/core/electrostatics/coulomb.hpp @@ -54,7 +54,9 @@ using ElectrostaticsActor = std::shared_ptr, std::shared_ptr, #endif // ESPRESSO_P3M +#ifdef ESPRESSO_GSL std::shared_ptr, +#endif // ESPRESSO_GSL #ifdef ESPRESSO_SCAFACOS std::shared_ptr, #endif // ESPRESSO_SCAFACOS @@ -81,7 +83,9 @@ struct has_pressure : std::false_type {}; #ifdef ESPRESSO_SCAFACOS template <> struct has_pressure : std::false_type {}; #endif // ESPRESSO_SCAFACOS +#ifdef ESPRESSO_GSL template <> struct has_pressure : std::false_type {}; +#endif // ESPRESSO_GSL } // namespace traits } // namespace Coulomb diff --git a/src/core/electrostatics/coulomb_inline.hpp b/src/core/electrostatics/coulomb_inline.hpp index 9da5c7254d4..4c786b09f0f 100644 --- a/src/core/electrostatics/coulomb_inline.hpp +++ b/src/core/electrostatics/coulomb_inline.hpp @@ -137,12 +137,14 @@ struct ShortRangeEnergyKernel { }}; } #endif // ESPRESSO_P3M +#ifdef ESPRESSO_GSL result_type operator()(std::shared_ptr const &actor) const { return kernel_type{ [&actor](Utils::Vector3d const &, Utils::Vector3d const &, double q1q2, Utils::Vector3d const &d, double dist) { return actor->pair_energy(q1q2, d, dist); }}; } +#endif // ESPRESSO_GSL #endif // ESPRESSO_ELECTROSTATICS }; diff --git a/src/core/electrostatics/mmm-modpsi.cpp b/src/core/electrostatics/mmm-modpsi.cpp index 5d8de798d4e..6f3fa1f4c78 100644 --- a/src/core/electrostatics/mmm-modpsi.cpp +++ b/src/core/electrostatics/mmm-modpsi.cpp @@ -19,12 +19,13 @@ * along with this program. If not, see . */ -#include "config/config.hpp" +#include -#ifdef ESPRESSO_ELECTROSTATICS +#if defined(ESPRESSO_ELECTROSTATICS) and defined(ESPRESSO_GSL) #include "mmm1d.hpp" -#include "specfunc.hpp" + +#include #include #include @@ -42,7 +43,7 @@ static void preparePolygammaEven(int n, double binom, series[0] = 2. * (1. - std::numbers::egamma); for (int order = 1;; order += 1) { auto const x_order = static_cast(2 * order); - auto const coeff = -2 * hzeta(x_order + 1, 2); + auto const coeff = -2. * gsl_sf_hzeta(x_order + 1., 2.); if (fabs(maxx * coeff) * (4.0 / 3.0) < round_error_prec) break; series.push_back(coeff); @@ -51,21 +52,21 @@ static void preparePolygammaEven(int n, double binom, } } else { // even, n > 0 - double maxx = 1; - double pref = 2; + double maxx = 1.; + double pref = 2.; for (int order = 0;; order++) { // only even exponents of x auto const x_order = static_cast(2 * order); - auto const coeff = pref * hzeta(1 + deriv + x_order, 2); + auto const coeff = pref * gsl_sf_hzeta(1. + deriv + x_order, 2.); if ((fabs(maxx * coeff) * (4.0 / 3.0) < round_error_prec) && (x_order > deriv)) break; series.push_back(-binom * coeff); maxx *= 0.25; - pref *= (1.0 + deriv / (x_order + 1)); - pref *= (1.0 + deriv / (x_order + 2)); + pref *= (1.0 + deriv / (x_order + 1.)); + pref *= (1.0 + deriv / (x_order + 2.)); } } } @@ -80,7 +81,7 @@ static void preparePolygammaOdd(int n, double binom, for (int order = 0;; order++) { // only odd exponents of x auto const x_order = static_cast(2 * order + 1); - auto const coeff = pref * hzeta(1 + deriv + x_order, 2); + auto const coeff = pref * gsl_sf_hzeta(1. + deriv + x_order, 2.); if ((fabs(maxx * coeff) * (4.0 / 3.0) < round_error_prec) && (x_order > deriv)) break; @@ -109,4 +110,4 @@ void CoulombMMM1D::create_mod_psi_up_to(int new_n) { } } -#endif // ESPRESSO_ELECTROSTATICS +#endif // defined(ESPRESSO_ELECTROSTATICS) and defined(ESPRESSO_GSL) diff --git a/src/core/electrostatics/mmm1d.cpp b/src/core/electrostatics/mmm1d.cpp index f35bf3ef85a..14214629816 100644 --- a/src/core/electrostatics/mmm1d.cpp +++ b/src/core/electrostatics/mmm1d.cpp @@ -19,9 +19,9 @@ * along with this program. If not, see . */ -#include "config/config.hpp" +#include -#ifdef ESPRESSO_ELECTROSTATICS +#if defined(ESPRESSO_ELECTROSTATICS) and defined(ESPRESSO_GSL) #include "electrostatics/mmm1d.hpp" @@ -32,12 +32,13 @@ #include "Particle.hpp" #include "cell_system/CellStructureType.hpp" #include "errorhandling.hpp" -#include "specfunc.hpp" #include "tuning.hpp" #include #include +#include + #include #include #include @@ -45,21 +46,14 @@ #include #include -/* if you define this feature, the Bessel functions are calculated up - * to machine precision, otherwise 10^-14, which should be - * definitely enough for daily life. */ -#ifndef ESPRESSO_MMM1D_MACHINE_PREC -#define K0 LPK0 -#define K1 LPK1 -#endif - static auto far_error(int P, double minrad, Utils::Vector3d const &box_l_inv) { auto const wavenumber = 2. * std::numbers::pi * box_l_inv[2]; // this uses an upper bound to all force components and the potential auto const rhores = wavenumber * minrad; auto const pref = 4. * box_l_inv[2] * std::max(1., wavenumber); - return pref * K1(rhores * P) * exp(rhores) / rhores * (P - 1. + 1. / rhores); + return pref * gsl_sf_bessel_K1(rhores * P) * exp(rhores) / rhores * + (P - 1. + 1. / rhores); } static auto determine_minrad(double maxPWerror, int P, @@ -93,6 +87,15 @@ static auto determine_minrad(double maxPWerror, int P, return 0.5 * (rmin + rmax); } +static double evaluateAsTaylorSeriesAt(std::span series, + double x) { + assert(not series.empty()); + auto const value = std::accumulate( + series.rbegin(), series.rend(), 0., + [x](auto const &acc, auto const &coeff) { return acc * x + coeff; }); + return value; +} + /** Modified polygamma for even order 2*n, n >= 0 */ static double mod_psi_even(auto const &modPsi, int n, double x) { return evaluateAsTaylorSeriesAt(modPsi[2 * n], x * x); @@ -250,12 +253,8 @@ Utils::Vector3d CoulombMMM1D::pair_force(double q1q2, Utils::Vector3d const &d, break; auto const fq = c_2pi * bp; -#ifdef ESPRESSO_MMM1D_MACHINE_PREC - auto const k0 = K0(fq * rxy_d); - auto const k1 = K1(fq * rxy_d); -#else - auto const [k0, k1] = LPK01(fq * rxy_d); -#endif + auto const k0 = gsl_sf_bessel_K0(fq * rxy_d); + auto const k1 = gsl_sf_bessel_K1(fq * rxy_d); sr += bp * k1 * cos(fq * z_d); sz += bp * k0 * sin(fq * z_d); } @@ -326,7 +325,7 @@ double CoulombMMM1D::pair_energy(double const q1q2, Utils::Vector3d const &d, break; auto const fq = c_2pi * bp; - energy += K0(fq * rxy_d) * cos(fq * z_d); + energy += gsl_sf_bessel_K0(fq * rxy_d) * cos(fq * z_d); } energy *= 4. * box_geo.length_inv()[2]; } @@ -382,4 +381,4 @@ void CoulombMMM1D::tune() { system.on_coulomb_change(); } -#endif // ESPRESSO_ELECTROSTATICS +#endif // defined(ESPRESSO_ELECTROSTATICS) and defined(ESPRESSO_GSL) diff --git a/src/core/electrostatics/mmm1d.hpp b/src/core/electrostatics/mmm1d.hpp index 80385cd20d3..5c10b78f202 100644 --- a/src/core/electrostatics/mmm1d.hpp +++ b/src/core/electrostatics/mmm1d.hpp @@ -35,9 +35,9 @@ #pragma once -#include "config/config.hpp" +#include -#ifdef ESPRESSO_ELECTROSTATICS +#if defined(ESPRESSO_ELECTROSTATICS) and defined(ESPRESSO_GSL) #include "electrostatics/actor.hpp" @@ -130,4 +130,4 @@ struct CoulombMMM1D : public Coulomb::Actor { void sanity_checks_cell_structure() const; }; -#endif // ESPRESSO_ELECTROSTATICS +#endif // defined(ESPRESSO_ELECTROSTATICS) and defined(ESPRESSO_GSL) diff --git a/src/core/electrostatics/specfunc.cpp b/src/core/electrostatics/specfunc.cpp deleted file mode 100644 index ec9d70f1cca..00000000000 --- a/src/core/electrostatics/specfunc.cpp +++ /dev/null @@ -1,489 +0,0 @@ -/* - * Copyright (C) 2010-2022 The ESPResSo project - * Copyright (C) 2002,2003,2004,2005,2006,2007,2008,2009,2010 - * Max-Planck-Institute for Polymer Research, Theory Group - * - * This file is part of ESPResSo. - * - * ESPResSo is free software: you can redistribute it and/or modify - * it under the terms of the GNU General Public License as published by - * the Free Software Foundation, either version 3 of the License, or - * (at your option) any later version. - * - * ESPResSo 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 General Public License for more details. - * - * You should have received a copy of the GNU General Public License - * along with this program. If not, see . - */ - -/* Original gsl header - * specfunc/bessel_K0.cpp - * - * Copyright (C) 1996, 1997, 1998, 1999, 2000 Gerard Jungman - * - * This program is free software; you can redistribute it and/or modify - * it under the terms of the GNU General Public License as published by - * the Free Software Foundation; either version 2 of the License, or (at - * your option) any later version. - * - * This program 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 - * General Public License for more details. - * - * You should have received a copy of the GNU General Public License - * along with this program; if not, write to the Free Software - * Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA. - */ - -/* Original Author: G. Jungman */ - -/** @file - * Implementation of @ref specfunc.hpp. - */ - -#include "specfunc.hpp" - -#include - -#include -#include -#include -#include - -/************************************************ - * chebychev expansions - ************************************************/ -/* Note that the first coefficient already includes the constant offsets */ - -/** @name Chebyshev expansions based on SLATEC bk0(), bk0e() */ -/**@{*/ -/** Series for @c bk0. - * On the interval 0. to 4.00000d+00 - * | Label | Value | - * | ---------------------------: | :------- | - * | with weighted error | 3.57e-19 | - * | log weighted error | 18.45 | - * | significant figures required | 17.99 | - * | decimal places required | 18.97 | - */ -static double bk0_cs[11] = { - -.5 - 0.03532739323390276872, 0.3442898999246284869, - 0.03597993651536150163, 0.00126461541144692592, - 0.00002286212103119451, 0.00000025347910790261, - 0.00000000190451637722, 0.00000000001034969525, - 0.00000000000004259816, 0.00000000000000013744, - 0.00000000000000000035}; - -/** Series for @c ak0. - * On the interval 1.25000d-01 to 5.00000d-01 - * | Label | Value | - * | ---------------------------: | :------- | - * | with weighted error | 5.34e-17 | - * | log weighted error | 16.27 | - * | significant figures required | 14.92 | - * | decimal places required | 16.89 | - */ -static double ak0_cs[17] = { - 2.5 - 0.07643947903327941, -0.02235652605699819, 0.00077341811546938, - -0.00004281006688886, 0.00000308170017386, -0.00000026393672220, - 0.00000002563713036, -0.00000000274270554, 0.00000000031694296, - -0.00000000003902353, 0.00000000000506804, -0.00000000000068895, - 0.00000000000009744, -0.00000000000001427, 0.00000000000000215, - -0.00000000000000033, 0.00000000000000005}; - -/** Series for @c ak02. - * On the interval 0. to 1.25000d-01 - * | Label | Value | - * | ---------------------------: | :------- | - * | with weighted error | 2.34e-17 | - * | log weighted error | 16.63 | - * | significant figures required | 14.67 | - * | decimal places required | 17.20 | - */ -static double ak02_cs[14] = { - 2.5 - 0.01201869826307592, -0.00917485269102569, 0.00014445509317750, - -0.00000401361417543, 0.00000015678318108, -0.00000000777011043, - 0.00000000046111825, -0.00000000003158592, 0.00000000000243501, - -0.00000000000020743, 0.00000000000001925, -0.00000000000000192, - 0.00000000000000020, -0.00000000000000002}; -/**@}*/ - -/** @name Chebyshev expansions based on SLATEC besi0() */ -/**@{*/ -/** Series for @c bi0. - * On the interval 0. to 9.00000d+00 - * | Label | Value | - * | ---------------------------: | :------- | - * | with weighted error | 2.46e-18 | - * | log weighted error | 17.61 | - * | significant figures required | 17.90 | - * | decimal places required | 18.15 | - */ -static double bi0_cs[12] = { - 5.5 - .07660547252839144951, 1.92733795399380827000, .22826445869203013390, - .01304891466707290428, .00043442709008164874, .00000942265768600193, - .00000014340062895106, .00000000161384906966, .00000000001396650044, - .00000000000009579451, .00000000000000053339, .00000000000000000245}; -/**@}*/ - -/** @name Chebyshev expansions based on SLATEC besk1(), besk1e() */ -/**@{*/ -/** Series for @c bk1. - * On the interval 0. to 4.00000d+00 - * | Label | Value | - * | ---------------------------: | :------- | - * | with weighted error | 7.02e-18 | - * | log weighted error | 17.15 | - * | significant figures required | 16.73 | - * | decimal places required | 17.67 | - */ -static double bk1_cs[11] = { - 1.5 + 0.0253002273389477705, -0.3531559607765448760, -0.1226111808226571480, - -0.0069757238596398643, -0.0001730288957513052, -0.0000024334061415659, - -0.0000000221338763073, -0.0000000001411488392, -0.0000000000006666901, - -0.0000000000000024274, -0.0000000000000000070}; - -/** Series for @c ak1. - * On the interval 1.25000d-01 to 5.00000d-01 - * | Label | Value | - * | ---------------------------: | :------- | - * | with weighted error | 6.06e-17 | - * | log weighted error | 16.22 | - * | significant figures required | 15.41 | - * | decimal places required | 16.83 | - */ -static double ak1_cs[17] = { - 2.5 + 0.27443134069738830, 0.07571989953199368, -0.00144105155647540, - 0.00006650116955125, -0.00000436998470952, 0.00000035402774997, - -0.00000003311163779, 0.00000000344597758, -0.00000000038989323, - 0.00000000004720819, -0.00000000000604783, 0.00000000000081284, - -0.00000000000011386, 0.00000000000001654, -0.00000000000000248, - 0.00000000000000038, -0.00000000000000006}; - -/** Series for @c ak12. - * On the interval 0. to 1.25000d-01 - * | Label | Value | - * | ---------------------------: | :------- | - * | with weighted error | 2.58e-17 | - * | log weighted error | 16.59 | - * | significant figures required | 15.22 | - * | decimal places required | 17.16 | - */ -static double ak12_cs[14] = { - 2.5 + 0.06379308343739001, 0.02832887813049721, -0.00024753706739052, - 0.00000577197245160, -0.00000020689392195, 0.00000000973998344, - -0.00000000055853361, 0.00000000003732996, -0.00000000000282505, - 0.00000000000023720, -0.00000000000002176, 0.00000000000000215, - -0.00000000000000022, 0.00000000000000002}; -/**@}*/ - -/** @name Chebyshev expansions based on SLATEC besi1(), besi1e() */ -/**@{*/ -/** Series for @c bi1. - * On the interval 0. to 9.00000d+00 - * | Label | Value | - * | ---------------------------: | :------- | - * | with weighted error | 2.40e-17 | - * | log weighted error | 16.62 | - * | significant figures required | 16.23 | - * | decimal places required | 17.14 | - */ -static double bi1_cs[11] = { - 1.75 - 0.001971713261099859, 0.407348876675464810, 0.034838994299959456, - 0.001545394556300123, 0.000041888521098377, 0.000000764902676483, - 0.000000010042493924, 0.000000000099322077, 0.000000000000766380, - 0.000000000000004741, 0.000000000000000024}; -/**@}*/ - -/** Coefficients for Maclaurin summation in hzeta(). Evaluated as inverse - * numbers, i.e. @f$ \displaystyle\frac{B_{2j}}{(2j)!} @f$. - */ -static double const hzeta_c[15] = { - 1.00000000000000000000000000000, 0.083333333333333333333333333333, - -0.00138888888888888888888888888889, 0.000033068783068783068783068783069, - -8.2671957671957671957671957672e-07, 2.0876756987868098979210090321e-08, - -5.2841901386874931848476822022e-10, 1.3382536530684678832826980975e-11, - -3.3896802963225828668301953912e-13, 8.5860620562778445641359054504e-15, - -2.1748686985580618730415164239e-16, 5.5090028283602295152026526089e-18, - -1.3954464685812523340707686264e-19, 3.5347070396294674716932299778e-21, - -8.9535174270375468504026113181e-23}; - -double hzeta(double s, double q) { - constexpr auto max_bits = 54.0; - constexpr auto jmax = 12; - constexpr auto kmax = 10; - - if ((s > max_bits and q < 1.0) or (s > 0.5 * max_bits and q < 0.25)) - return std::pow(q, -s); - if (s > 0.5 * max_bits and q < 1.0) { - auto const p1 = std::pow(q, -s); - auto const p2 = std::pow(q / (1.0 + q), s); - auto const p3 = std::pow(q / (2.0 + q), s); - return p1 * (1.0 + p2 + p3); - } - /** Euler-Maclaurin summation formula from @cite moshier89a p. 400, with - * several typo corrections. - */ - auto const kmax_q = static_cast(kmax) + q; - auto const pmax = std::pow(kmax_q, -s); - auto scp = s; - auto pcp = pmax / kmax_q; - auto ans = pmax * (kmax_q / (s - 1.0) + 0.5); - - for (int k = 0; k < kmax; k++) - ans += std::pow(static_cast(k) + q, -s); - - for (int j = 0; j <= jmax; j++) { - auto const delta = hzeta_c[j + 1] * scp * pcp; - ans += delta; - scp *= (s + 2. * j + 1.) * (s + 2. * j + 2.); - pcp /= Utils::sqr(static_cast(kmax) + q); - } - - return ans; -} - -double K0(double x) { - if (x <= 2.0) { - auto const c = evaluateAsChebychevSeriesAt(bk0_cs, 0.5 * x * x - 1.0); - auto const i0 = evaluateAsChebychevSeriesAt(bi0_cs, x * x / 4.5 - 1.0); - return (-std::log(x) + std::numbers::ln2) * i0 + c; - } - auto const c = - (x <= 8.0) ? evaluateAsChebychevSeriesAt(ak0_cs, (16.0 / x - 5.0) / 3.0) - : evaluateAsChebychevSeriesAt(ak02_cs, 16.0 / x - 1.0); - return std::exp(-x) * c / std::sqrt(x); -} - -double K1(double x) { - if (x <= 2.0) { - auto const c = evaluateAsChebychevSeriesAt(bk1_cs, 0.5 * x * x - 1.0); - auto const i1 = x * evaluateAsChebychevSeriesAt(bi1_cs, x * x / 4.5 - 1.0); - return (std::log(x) - std::numbers::ln2) * i1 + c / x; - } - auto const c = - (x <= 8.0) ? evaluateAsChebychevSeriesAt(ak1_cs, (16.0 / x - 5.0) / 3.0) - : evaluateAsChebychevSeriesAt(ak12_cs, 16.0 / x - 1.0); - return std::exp(-x) * c / std::sqrt(x); -} - -/*********************************************************** - * optimized K0/1 implementations for 10^(-14) precision - ***********************************************************/ - -/** necessary orders for K0/1 from 2 up to 22 for 10^-14 precision. Note that - * at 8 the expansion changes. From 23 to 26 order 2 is used, above order 1. - * For the latter cases, separate implementations are necessary. - */ -static int ak01_orders[] = { - /* 2 - 8 */ - 11, 11, 10, 10, 9, 9, - /* 8 - 26 */ - 6, 6, 5, 5, 5, 4, 4, 4, 3, 3, 2, 2, 2, 2, 2}; - -double LPK0(double x) { - if (x >= 27.) { - auto const tmp = .5 * std::exp(-x) / std::sqrt(x); - return tmp * ak0_cs[0]; - } - if (x >= 23.) { - auto const tmp = std::exp(-x) / std::sqrt(x); - auto const xx = (16. / 3.) / x - 5. / 3.; - return tmp * (xx * ak0_cs[1] + 0.5 * ak0_cs[0]); - } - if (x > 2.) { - int j = ak01_orders[static_cast(x) - 2]; - double x2; - double *s0; - if (x <= 8.) { - s0 = ak0_cs; - x2 = (2. * 16. / 3.) / x - 2. * 5. / 3.; - } else { - s0 = ak02_cs; - x2 = (2. * 16.) / x - 2.; - } - auto dd0 = s0[j]; - auto d0 = x2 * dd0 + s0[j - 1]; - for (j -= 2; j >= 1; j--) { - auto const tmp0 = d0; - d0 = x2 * d0 - dd0 + s0[j]; - dd0 = tmp0; - } - auto const tmp = std::exp(-x) / std::sqrt(x); - return tmp * (0.5 * (s0[0] + x2 * d0) - dd0); - } - /* x <= 2 */ - { - /* I0/I1 series */ - int j = 10; - auto x2 = (2. / 4.5) * x * x - 2.; - auto dd0 = bi0_cs[j]; - auto d0 = x2 * dd0 + bi0_cs[j - 1]; - for (j -= 2; j >= 1; j--) { - auto const tmp0 = d0; - d0 = x2 * d0 - dd0 + bi0_cs[j]; - dd0 = tmp0; - } - auto const tmp = std::log(x) - std::numbers::ln2; - auto const ret = -tmp * (0.5 * (bi0_cs[0] + x2 * d0) - dd0); - - /* K0/K1 correction */ - j = 9; - x2 = x * x - 2.; - dd0 = bk0_cs[j]; - d0 = x2 * dd0 + bk0_cs[j - 1]; - for (j -= 2; j >= 1; j--) { - auto const tmp0 = d0; - d0 = x2 * d0 - dd0 + bk0_cs[j]; - dd0 = tmp0; - } - return ret + (0.5 * (x2 * d0 + bk0_cs[0]) - dd0); - } -} - -double LPK1(double x) { - if (x >= 27.) { - auto const tmp = .5 * std::exp(-x) / std::sqrt(x); - return tmp * ak1_cs[0]; - } - if (x >= 23.) { - auto const tmp = std::exp(-x) / std::sqrt(x); - auto const xx = (16. / 3.) / x - 5. / 3.; - return tmp * (xx * ak1_cs[1] + 0.5 * ak1_cs[0]); - } - if (x > 2.) { - int j = ak01_orders[static_cast(x) - 2]; - double x2; - double *s1; - if (x <= 8.) { - s1 = ak1_cs; - x2 = (2. * 16. / 3.) / x - 2. * 5. / 3.; - } else { - s1 = ak12_cs; - x2 = (2. * 16.) / x - 2.; - } - auto dd1 = s1[j]; - auto d1 = x2 * dd1 + s1[j - 1]; - for (j -= 2; j >= 1; j--) { - auto const tmp1 = d1; - d1 = x2 * d1 - dd1 + s1[j]; - dd1 = tmp1; - } - auto const tmp = std::exp(-x) / std::sqrt(x); - return tmp * (0.5 * (s1[0] + x2 * d1) - dd1); - } - /* x <= 2 */ - { - /* I0/I1 series */ - int j = 10; - auto x2 = (2. / 4.5) * x * x - 2.; - auto dd1 = bi1_cs[j]; - auto d1 = x2 * dd1 + bi1_cs[j - 1]; - for (j -= 2; j >= 1; j--) { - auto const tmp1 = d1; - d1 = x2 * d1 - dd1 + bi1_cs[j]; - dd1 = tmp1; - } - auto const tmp = std::log(x) - std::numbers::ln2; - auto const ret = x * tmp * (0.5 * (bi1_cs[0] + x2 * d1) - dd1); - - /* K0/K1 correction */ - j = 9; - x2 = x * x - 2.; - dd1 = bk1_cs[j]; - d1 = x2 * dd1 + bk1_cs[j - 1]; - for (j -= 2; j >= 1; j--) { - auto const tmp1 = d1; - d1 = x2 * d1 - dd1 + bk1_cs[j]; - dd1 = tmp1; - } - return ret + (0.5 * (x2 * d1 + bk1_cs[0]) - dd1) / x; - } -} - -std::pair LPK01(double x) { - if (x >= 27.) { - auto const tmp = .5 * std::exp(-x) / std::sqrt(x); - auto const k0 = tmp * ak0_cs[0]; - auto const k1 = tmp * ak1_cs[0]; - return {k0, k1}; - } - if (x >= 23.) { - auto const tmp = std::exp(-x) / std::sqrt(x); - auto const xx = (16. / 3.) / x - 5. / 3.; - auto const k0 = tmp * (xx * ak0_cs[1] + 0.5 * ak0_cs[0]); - auto const k1 = tmp * (xx * ak1_cs[1] + 0.5 * ak1_cs[0]); - return {k0, k1}; - } - if (x > 2.) { - int j = ak01_orders[static_cast(x) - 2]; - double x2; - double *s0, *s1; - if (x <= 8.) { - s0 = ak0_cs; - s1 = ak1_cs; - x2 = (2. * 16. / 3.) / x - 2. * 5. / 3.; - } else { - s0 = ak02_cs; - s1 = ak12_cs; - x2 = (2. * 16.) / x - 2.; - } - auto dd0 = s0[j]; - auto dd1 = s1[j]; - auto d0 = x2 * dd0 + s0[j - 1]; - auto d1 = x2 * dd1 + s1[j - 1]; - for (j -= 2; j >= 1; j--) { - auto const tmp0 = d0, tmp1 = d1; - d0 = x2 * d0 - dd0 + s0[j]; - d1 = x2 * d1 - dd1 + s1[j]; - dd0 = tmp0; - dd1 = tmp1; - } - auto const tmp = std::exp(-x) / std::sqrt(x); - auto const k0 = tmp * (0.5 * (s0[0] + x2 * d0) - dd0); - auto const k1 = tmp * (0.5 * (s1[0] + x2 * d1) - dd1); - return {k0, k1}; - } - /* x <= 2 */ - { - /* I0/I1 series */ - int j = 10; - auto x2 = (2. / 4.5) * x * x - 2.; - auto dd0 = bi0_cs[j]; - auto dd1 = bi1_cs[j]; - auto d0 = x2 * dd0 + bi0_cs[j - 1]; - auto d1 = x2 * dd1 + bi1_cs[j - 1]; - for (j -= 2; j >= 1; j--) { - auto const tmp0 = d0, tmp1 = d1; - d0 = x2 * d0 - dd0 + bi0_cs[j]; - d1 = x2 * d1 - dd1 + bi1_cs[j]; - dd0 = tmp0; - dd1 = tmp1; - } - auto const tmp = std::log(x) - std::numbers::ln2; - auto k0 = -tmp * (0.5 * (bi0_cs[0] + x2 * d0) - dd0); - auto k1 = x * tmp * (0.5 * (bi1_cs[0] + x2 * d1) - dd1); - - /* K0/K1 correction */ - j = 9; - x2 = x * x - 2.; - dd0 = bk0_cs[j]; - dd1 = bk1_cs[j]; - d0 = x2 * dd0 + bk0_cs[j - 1]; - d1 = x2 * dd1 + bk1_cs[j - 1]; - for (j -= 2; j >= 1; j--) { - auto const tmp0 = d0, tmp1 = d1; - d0 = x2 * d0 - dd0 + bk0_cs[j]; - d1 = x2 * d1 - dd1 + bk1_cs[j]; - dd0 = tmp0; - dd1 = tmp1; - } - k0 += (0.5 * (x2 * d0 + bk0_cs[0]) - dd0); - k1 += (0.5 * (x2 * d1 + bk1_cs[0]) - dd1) / x; - return {k0, k1}; - } -} diff --git a/src/core/electrostatics/specfunc.hpp b/src/core/electrostatics/specfunc.hpp deleted file mode 100644 index 225b3d33dc0..00000000000 --- a/src/core/electrostatics/specfunc.hpp +++ /dev/null @@ -1,117 +0,0 @@ -/* - * Copyright (C) 2010-2022 The ESPResSo project - * Copyright (C) 2002,2003,2004,2005,2006,2007,2008,2009,2010 - * Max-Planck-Institute for Polymer Research, Theory Group - * - * This file is part of ESPResSo. - * - * ESPResSo is free software: you can redistribute it and/or modify - * it under the terms of the GNU General Public License as published by - * the Free Software Foundation, either version 3 of the License, or - * (at your option) any later version. - * - * ESPResSo 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 General Public License for more details. - * - * You should have received a copy of the GNU General Public License - * along with this program. If not, see . - */ - -/** @file - * This file contains implementations for some special functions which are - * needed by the MMM family of algorithms. These are the modified Hurwitz - * zeta function and the modified Bessel functions of second kind. The - * implementations are based on the GSL code (see @ref specfunc.cpp for the - * original GSL header). - * - * The Hurwitz zeta function is evaluated using the Euler-Maclaurin summation - * formula, the Bessel functions are evaluated using several different - * Chebychev expansions. Both achieve a precision of nearly machine precision, - * which is no problem for the Hurwitz zeta function, which is only used when - * determining the coefficients for the modified polygamma functions. - */ - -#pragma once - -#include -#include -#include -#include - -/** Hurwitz zeta function. This function was taken from the GSL code. */ -double hzeta(double order, double x); - -/** Modified Bessel function of second kind, order 0. This function was taken - * from the GSL code. Precise roughly up to machine precision. - * It is 16 times faster than std::cyl_bessel_k. - * If @c MMM1D_MACHINE_PREC is not defined, @ref LPK0 is used instead. - */ -double K0(double x); - -/** Modified Bessel function of second kind, order 1. This function was taken - * from the GSL code. Precise roughly up to machine precision. - * If @c MMM1D_MACHINE_PREC is not defined, @ref LPK1 is used instead. - */ -double K1(double x); - -/** Modified Bessel function of second kind, order 0, low precision. - * The implementation has an absolute precision of around 10^(-14), which is - * comparable to the relative precision sqrt implementation of current - * hardware in the ranges @f$ ]0, 8[ @f$ and @f$ ]8, 23[ @f$. Above 23, - * the precision starts to degrade, and above 27 the result drifts and - * slowly converges to 96% of the real value. - * It is 25 times faster than std::cyl_bessel_k - * and 1.5 times faster than @ref K0. - */ -double LPK0(double x); - -/** Modified Bessel function of second kind, order 1, low precision. - * The implementation has an absolute precision of around 10^(-14), which is - * comparable to the relative precision sqrt implementation of current - * hardware in the ranges @f$ ]0, 8[ @f$ and @f$ ]8, 23[ @f$. Above 23, - * the precision starts to degrade, and above 27 the result drifts and - * slowly converges to 111% of the real value. - * It is 25 times faster than std::cyl_bessel_k - * and 1.5 times faster than @ref K1. - */ -double LPK1(double x); - -/** Modified Bessel functions of second kind, order 0 and 1, low precision. - * The implementation has an absolute precision of around 10^(-14), which is - * comparable to the relative precision sqrt implementation of current - * hardware. - */ -std::pair LPK01(double x); - -/** Evaluate the polynomial interpreted as a Taylor series via the - * Horner scheme. - */ -inline double evaluateAsTaylorSeriesAt(std::span series, - double x) { - assert(not series.empty()); - auto const value = std::accumulate( - series.rbegin(), series.rend(), 0., - [x](auto const &acc, auto const &coeff) { return acc * x + coeff; }); - return value; -} - -/** Evaluate the polynomial interpreted as a Chebychev series. Requires a - * series with at least three coefficients, i.e. no linear approximations! - */ -inline double evaluateAsChebychevSeriesAt(std::span series, - double x) { - assert(series.size() >= 3); - - auto const *c = series.data(); - auto const x2 = 2.0 * x; - auto dd = c[series.size() - 1]; - auto d = x2 * dd + c[series.size() - 2]; - for (auto j = static_cast(series.size()) - 3; j >= 1; j--) { - auto const tmp = d; - d = x2 * d - dd + c[j]; - dd = tmp; - } - return x * d - dd + 0.5 * c[0]; -} diff --git a/src/core/unit_tests/CMakeLists.txt b/src/core/unit_tests/CMakeLists.txt index e1586eb7f56..f3d4a39b14b 100644 --- a/src/core/unit_tests/CMakeLists.txt +++ b/src/core/unit_tests/CMakeLists.txt @@ -82,11 +82,6 @@ espresso_unit_test(SRC BondList_test.cpp DEPENDS espresso::core) espresso_unit_test(SRC energy_test.cpp DEPENDS espresso::core) espresso_unit_test(SRC bonded_interactions_map_test.cpp DEPENDS espresso::core) espresso_unit_test(SRC bond_breakage_test.cpp DEPENDS espresso::core) -if(NOT CMAKE_CXX_COMPILER_ID STREQUAL "AppleClang") - # AppleClang doesn't implement C++17's mathematical special functions - espresso_unit_test(SRC specfunc_test.cpp DEPENDS espresso::utils - espresso::core) -endif() if(ESPRESSO_BUILD_WITH_WALBERLA) espresso_unit_test(SRC lb_particle_coupling_test.cpp DEPENDS espresso::core Boost::mpi MPI::MPI_CXX NUM_PROC 2) diff --git a/src/core/unit_tests/specfunc_test.cpp b/src/core/unit_tests/specfunc_test.cpp deleted file mode 100644 index 2a9c16f2071..00000000000 --- a/src/core/unit_tests/specfunc_test.cpp +++ /dev/null @@ -1,122 +0,0 @@ -/* - * Copyright (C) 2022 The ESPResSo project - * - * This file is part of ESPResSo. - * - * ESPResSo is free software: you can redistribute it and/or modify - * it under the terms of the GNU General Public License as published by - * the Free Software Foundation, either version 3 of the License, or - * (at your option) any later version. - * - * ESPResSo 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 General Public License for more details. - * - * You should have received a copy of the GNU General Public License - * along with this program. If not, see . - */ - -#define BOOST_TEST_MODULE special functions test -#define BOOST_TEST_DYN_LINK -#include - -#include "electrostatics/specfunc.hpp" - -#include -#include - -auto constexpr eps = 8. * 100. * std::numeric_limits::epsilon(); - -BOOST_AUTO_TEST_CASE(hurwitz_zeta_function) { - constexpr auto max_bits = 54.0; - // test cases where an exact, closed-form expression exists - auto delta = 0.025; - auto x = delta / 2.; - while (x < 0.25) { - auto order = max_bits / 2. + 0.01; - BOOST_TEST_INFO("with parameter x = " << x); - BOOST_CHECK_CLOSE(hzeta(order, x), std::pow(x, -order), eps); - x += delta; - } - x = delta / 2.; - while (x < 1.) { - auto order = max_bits + 0.01; - BOOST_TEST_INFO("with parameter x = " << x); - BOOST_CHECK_CLOSE(hzeta(order, x), std::pow(x, -order), eps); - x += delta; - } - x = delta / 2.; - while (x < 1.) { - auto order = max_bits / 2. + 0.01; - auto ref = std::pow(x, -order); - ref *= (1. + std::pow(x / (1. + x), order) + std::pow(x / (2. + x), order)); - BOOST_TEST_INFO("with parameter x = " << x); - BOOST_CHECK_CLOSE(hzeta(order, x), ref, eps); - x += delta; - } - // Hurwitz zeta is a generalization of the Riemann zeta - for (auto order = 2; order < 128; ++order) { - BOOST_TEST_INFO("with parameter order = " << order); - BOOST_CHECK_CLOSE(hzeta(order, 1.), std::riemann_zeta(order), 10. * eps); - } -} - -BOOST_AUTO_TEST_CASE(bessel_series_high_precision) { - auto delta = 0.02; - auto x = delta; - while (x < 20.) { - BOOST_TEST_INFO("with parameter x = " << x); - BOOST_CHECK_CLOSE(K0(x), std::cyl_bessel_k(0, x), 10. * eps); - BOOST_CHECK_CLOSE(K1(x), std::cyl_bessel_k(1, x), 10. * eps); - x += delta; - } - delta = 0.5; - while (x < 64.) { - BOOST_TEST_INFO("with parameter x = " << x); - BOOST_CHECK_CLOSE(K0(x), std::cyl_bessel_k(0, x), 10. * eps); - BOOST_CHECK_CLOSE(K1(x), std::cyl_bessel_k(1, x), 10. * eps); - x += delta; - } -} - -BOOST_AUTO_TEST_CASE(bessel_series_low_precision) { - auto const check = [](double x, double tol) { - BOOST_TEST_INFO("with parameter x = " << x); - BOOST_CHECK_CLOSE(LPK0(x), std::cyl_bessel_k(0, x), tol); - BOOST_TEST_INFO("with parameter x = " << x); - BOOST_CHECK_CLOSE(LPK1(x), std::cyl_bessel_k(1, x), tol); - BOOST_CHECK_CLOSE(LPK01(x).first, LPK0(x), eps); - BOOST_CHECK_CLOSE(LPK01(x).second, LPK1(x), eps); - }; - auto delta = 0.02; - auto x = delta; - // important breakpoints: x=2, x=8, x=23, x=27 (different interpolation) - while (x < 2.) { - check(x, 20. * eps); - x += delta; - } - while (x < 8.) { - check(x, 1e5 * eps); - x += delta; - } - delta = 0.04; - while (x < 13.) { - check(x, 1e8 * eps); - x += delta; - } - while (x < 23.) { - check(x, 2e10 * eps); - x += delta; - } - delta = 0.08; - while (x < 27.) { - check(x, 0.6); - x += delta; - } - delta = 0.64; - while (x < 64.) { - check(x, 20.); - x += delta; - } -} diff --git a/src/python/espressomd/CMakeLists.txt b/src/python/espressomd/CMakeLists.txt index adf8a41bfd5..4517c28b047 100644 --- a/src/python/espressomd/CMakeLists.txt +++ b/src/python/espressomd/CMakeLists.txt @@ -107,6 +107,17 @@ foreach(cython_file ${cython_SRC}) ${target} PROPERTIES SUFFIX ".so" LINK_FLAGS "-undefined dynamic_lookup") endif() set_target_properties(${target} PROPERTIES CXX_CLANG_TIDY "") + if(ESPRESSO_BUILD_WITH_SHARED_MEMORY_PARALLELISM) + if(CMAKE_CXX_COMPILER_ID STREQUAL CrayClang AND target STREQUAL + espressomd__init + AND DEFINED OpenMP_craymp_LIBRARY AND EXISTS ${OpenMP_craymp_LIBRARY}) + # the Cray OpenMP library must be loaded before dependent Cray libraries, + # including those imported from target OpenMP::OpenMP_CXX, so as to avoid + # "ImportError: cannot allocate memory in static TLS block" + target_link_libraries(${target} PRIVATE ${OpenMP_craymp_LIBRARY} + OpenMP::OpenMP_CXX) + endif() + endif() target_link_libraries(${target} PRIVATE espresso::config espresso::core espresso::script_interface) target_link_libraries(${target} PRIVATE espresso::compiler_flags) diff --git a/src/python/espressomd/electrostatics.py b/src/python/espressomd/electrostatics.py index 3fa321d7c0b..d6ec164c4aa 100644 --- a/src/python/espressomd/electrostatics.py +++ b/src/python/espressomd/electrostatics.py @@ -391,6 +391,7 @@ class MMM1D(ElectrostaticInteraction): """ _so_name = "Coulomb::CoulombMMM1D" _so_creation_policy = "GLOBAL" + _so_features = ("ELECTROSTATICS", "GSL") def default_params(self): return {"far_switch_radius": -1., diff --git a/src/script_interface/electrostatics/CoulombMMM1D.hpp b/src/script_interface/electrostatics/CoulombMMM1D.hpp index c3031fe2e0b..a46cf2fb017 100644 --- a/src/script_interface/electrostatics/CoulombMMM1D.hpp +++ b/src/script_interface/electrostatics/CoulombMMM1D.hpp @@ -19,9 +19,9 @@ #pragma once -#include "config/config.hpp" +#include -#ifdef ESPRESSO_ELECTROSTATICS +#if defined(ESPRESSO_ELECTROSTATICS) and defined(ESPRESSO_GSL) #include "Actor.hpp" @@ -69,4 +69,4 @@ class CoulombMMM1D : public Actor { } // namespace Coulomb } // namespace ScriptInterface -#endif // ESPRESSO_ELECTROSTATICS +#endif // defined(ESPRESSO_ELECTROSTATICS) and defined(ESPRESSO_GSL) diff --git a/src/script_interface/electrostatics/CoulombP3M.hpp b/src/script_interface/electrostatics/CoulombP3M.hpp index 3a34fbb6eb6..e31b7fcd8f7 100644 --- a/src/script_interface/electrostatics/CoulombP3M.hpp +++ b/src/script_interface/electrostatics/CoulombP3M.hpp @@ -87,7 +87,12 @@ class CoulombP3M : public Actor, ::CoulombP3M> { [this]() { return m_tuning.timings; }}, {"tune_limits", AutoParameter::read_only, [this]() { +#if defined(__clang__) and defined(__cray__) + auto const &range_min = m_tune_limits.first; + auto const &range_max = m_tune_limits.second; +#else auto const &[range_min, range_max] = m_tuning.limits; +#endif std::vector retval = { range_min ? Variant{*range_min} : Variant{None{}}, range_max ? Variant{*range_max} : Variant{None{}}, diff --git a/src/script_interface/electrostatics/initialize.cpp b/src/script_interface/electrostatics/initialize.cpp index 6ca2e08febd..35756408e23 100644 --- a/src/script_interface/electrostatics/initialize.cpp +++ b/src/script_interface/electrostatics/initialize.cpp @@ -57,7 +57,9 @@ void initialize(Utils::Factory *om) { "Coulomb::ElectrostaticLayerCorrection"); #endif // ESPRESSO_P3M om->register_new("Coulomb::ICCStar"); +#ifdef ESPRESSO_GSL om->register_new("Coulomb::CoulombMMM1D"); +#endif #ifdef ESPRESSO_SCAFACOS om->register_new("Coulomb::CoulombScafacos"); #endif diff --git a/src/script_interface/magnetostatics/DipolarP3M.hpp b/src/script_interface/magnetostatics/DipolarP3M.hpp index e77950c4b10..af570444a9f 100644 --- a/src/script_interface/magnetostatics/DipolarP3M.hpp +++ b/src/script_interface/magnetostatics/DipolarP3M.hpp @@ -89,7 +89,12 @@ class DipolarP3M : public Actor, ::DipolarP3M> { [this]() { return m_tuning.timings; }}, {"tune_limits", AutoParameter::read_only, [this]() { +#if defined(__clang__) and defined(__cray__) + auto const &range_min = m_tune_limits.first; + auto const &range_max = m_tune_limits.second; +#else auto const &[range_min, range_max] = m_tuning.limits; +#endif std::vector retval = { range_min ? Variant{*range_min} : Variant{None{}}, range_max ? Variant{*range_max} : Variant{None{}}, diff --git a/src/walberla_bridge/CMakeLists.txt b/src/walberla_bridge/CMakeLists.txt index 1ba53690c0e..7a9419e6e52 100644 --- a/src/walberla_bridge/CMakeLists.txt +++ b/src/walberla_bridge/CMakeLists.txt @@ -86,22 +86,21 @@ endfunction() add_library(espresso_walberla_compiler_flags INTERFACE) add_library(espresso::walberla::compiler_flags ALIAS espresso_walberla_compiler_flags) -target_link_libraries( - espresso_walberla_compiler_flags - INTERFACE - espresso::compiler_flags - $<$:espresso::avx_flags> - $<$:OpenMP::OpenMP_CXX> -) add_library(espresso_walberla_codegen_compiler_flags INTERFACE) add_library(espresso::walberla_codegen::compiler_flags ALIAS espresso_walberla_codegen_compiler_flags) # cmake-format: off +target_compile_options( + espresso_walberla_compiler_flags + INTERFACE + $<$:-Wno-reorder-ctor> +) target_compile_options( espresso_walberla_codegen_compiler_flags INTERFACE $<$:-Wno-cast-align> $<$:-Wno-shadow=compatible-local> + $<$:-Wno-reorder-ctor> # silence diagnostics from MPI vendors $<$,$>:-Wno-cast-function-type-mismatch> $<$,$>:-Wno-cast-function-type-mismatch> @@ -120,12 +119,25 @@ espresso_configure_walberla_target(espresso_walberla) espresso_configure_walberla_target(espresso_walberla_codegen) target_link_libraries( - espresso_walberla PUBLIC MPI::MPI_CXX espresso::utils - PRIVATE espresso::walberla::compiler_flags espresso::walberla_codegen - $<$:Heffte::Heffte>) + espresso_walberla + PUBLIC MPI::MPI_CXX espresso::utils + PRIVATE + espresso::compiler_flags + espresso::walberla::compiler_flags + espresso::walberla_codegen + $<$:Heffte::Heffte> + $<$:espresso::avx_flags> + $<$:OpenMP::OpenMP_CXX> +) target_link_libraries( - espresso_walberla_codegen PRIVATE espresso::walberla::compiler_flags - espresso::walberla_codegen::compiler_flags) + espresso_walberla_codegen + PRIVATE + espresso::compiler_flags + espresso::walberla::compiler_flags + espresso::walberla_codegen::compiler_flags + $<$:espresso::avx_flags> + $<$:OpenMP::OpenMP_CXX> +) function(espresso_set_walberla_codegen_kernel_source_properties source_file) # disable sanitizers on kernels due to excessive runtime and VRAM usage diff --git a/src/walberla_bridge/src/lattice_boltzmann/LBWalberlaImpl.hpp b/src/walberla_bridge/src/lattice_boltzmann/LBWalberlaImpl.hpp index b33336d8e14..7a0702f39a2 100644 --- a/src/walberla_bridge/src/lattice_boltzmann/LBWalberlaImpl.hpp +++ b/src/walberla_bridge/src/lattice_boltzmann/LBWalberlaImpl.hpp @@ -1572,9 +1572,8 @@ class LBWalberlaImpl : public LBWalberlaBase { if (auto const bci = get_block_interval( lattice, lower_corner, upper_corner, block_offset, block)) { - auto kernel = [&lattice, &block, &velocity, - this](unsigned const, unsigned const local_index, - Utils::Vector3i const &node) { + auto kernel = [&, this](unsigned const, unsigned const local_index, + Utils::Vector3i const &node) { auto const bc = get_block_and_cell(lattice, node, false); assert(bc->block->getAABB() == block.getAABB()); auto const &opt = velocity[local_index]; diff --git a/src/walberla_bridge/tests/CMakeLists.txt b/src/walberla_bridge/tests/CMakeLists.txt index a19af9fa645..a202e3ea940 100644 --- a/src/walberla_bridge/tests/CMakeLists.txt +++ b/src/walberla_bridge/tests/CMakeLists.txt @@ -34,9 +34,12 @@ function(ESPRESSO_ADD_TEST) ${TEST_NUM_THREADS} DEPENDS ${TEST_DEPENDS} + $<$:espresso::avx_flags> + $<$:OpenMP::OpenMP_CXX> espresso::walberla - espresso::utils + espresso::compiler_flags espresso::walberla::compiler_flags + espresso::utils ${WALBERLA_LIBS}) target_include_directories(${TEST_NAME} PRIVATE ${WALBERLA_INCLUDE_DIRS} ${walberla_BINARY_DIR}/src) diff --git a/testsuite/python/coulomb_interface.py b/testsuite/python/coulomb_interface.py index 6e2fe4d4348..37306291e12 100644 --- a/testsuite/python/coulomb_interface.py +++ b/testsuite/python/coulomb_interface.py @@ -92,6 +92,7 @@ def tearDown(self): prefactor=2., r_cut=1.5, cao=2, mesh=[8, 8, 8], alpha=12., accuracy=0.01, tune=False))) + @utx.skipIfMissingFeatures(["ELECTROSTATICS", "GSL"]) def test_mmm1d_cpu(self): self.system.periodicity = [False, False, True] self.system.cell_system.set_n_square() @@ -108,7 +109,7 @@ def test_mmm1d_cpu(self): with self.assertRaisesRegex(ValueError, f"Parameter '{key}' must be > 0"): espressomd.electrostatics.MMM1D(**invalid_params) - @utx.skipIfMissingFeatures(["P3M"]) + @utx.skipIfMissingFeatures(["P3M", "GSL"]) def test_solvers_rollback(self): # swapping two solvers should safely rollback to last valid solver self.system.periodicity = [False, False, True] @@ -125,6 +126,7 @@ def test_solvers_rollback(self): self.assertAlmostEqual( self.system.analysis.energy()["coulomb"], ref_energy, delta=1e-7) + @utx.skipIfMissingFeatures(["ELECTROSTATICS", "GSL"]) def test_charge_neutrality_check(self): self.system.part.add(pos=(0.0, 0.0, 0.0), q=1.) self.system.periodicity = [False, False, True] @@ -148,6 +150,7 @@ def test_charge_neutrality_check(self): self.assertFalse(actor.check_neutrality) self.assertIsNone(actor.charge_neutrality_tolerance) + @utx.skipIfMissingFeatures(["ELECTROSTATICS", "GSL"]) def test_mmm1d_cpu_tuning_exceptions(self): self.system.periodicity = [False, False, True] self.system.cell_system.set_n_square() diff --git a/testsuite/python/long_range_actors.py b/testsuite/python/long_range_actors.py index 74616e5634f..0425c968216 100644 --- a/testsuite/python/long_range_actors.py +++ b/testsuite/python/long_range_actors.py @@ -297,7 +297,7 @@ def test_dp3m_cpu_tuning_errors(self): prefactor=1., accuracy=1e-3) self.check_p3m_tuning_errors("magnetostatics", dp3m) - @utx.skipIfMissingFeatures(["ELECTROSTATICS"]) + @utx.skipIfMissingFeatures(["ELECTROSTATICS", "GSL"]) def test_mmm1d_cpu_exceptions(self): self.system.periodicity = (False, False, True) mmm1d = espressomd.electrostatics.MMM1D(prefactor=1., maxPWerror=1e-2) diff --git a/testsuite/python/mmm1d.py b/testsuite/python/mmm1d.py index e2f5b27d21d..0815cdf9779 100644 --- a/testsuite/python/mmm1d.py +++ b/testsuite/python/mmm1d.py @@ -23,7 +23,7 @@ import espressomd.electrostatics -@utx.skipIfMissingFeatures(["ELECTROSTATICS"]) +@utx.skipIfMissingFeatures(["ELECTROSTATICS", "GSL"]) class Test(ut.TestCase): system = espressomd.System(box_l=[10.0] * 3) diff --git a/testsuite/tutorials/CMakeLists.txt b/testsuite/tutorials/CMakeLists.txt index ae0a6180461..926bef764c8 100644 --- a/testsuite/tutorials/CMakeLists.txt +++ b/testsuite/tutorials/CMakeLists.txt @@ -58,6 +58,7 @@ tutorial_test(FILE test_electrodes_2.py) tutorial_test(FILE test_constant_pH__interactions.py) tutorial_test(FILE test_widom_insertion.py) tutorial_test(FILE test_grand_canonical_monte_carlo.py) +tutorial_test(FILE test_boltzmann_inversion.py) add_custom_target( check_tutorials JOB_POOL console diff --git a/testsuite/tutorials/test_boltzmann_inversion.py b/testsuite/tutorials/test_boltzmann_inversion.py new file mode 100644 index 00000000000..f0ba607cf79 --- /dev/null +++ b/testsuite/tutorials/test_boltzmann_inversion.py @@ -0,0 +1,39 @@ +# +# Copyright (C) 2025 The ESPResSo project +# +# This file is part of ESPResSo. +# +# ESPResSo is free software: you can redistribute it and/or modify +# it under the terms of the GNU General Public License as published by +# the Free Software Foundation, either version 3 of the License, or +# (at your option) any later version. +# +# ESPResSo 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 General Public License for more details. +# +# You should have received a copy of the GNU General Public License +# along with this program. If not, see . +# + +import unittest as ut +import importlib_wrapper +import numpy as np + +tutorial, skipIfMissingFeatures = importlib_wrapper.configure_and_import( + "@TUTORIALS_DIR@/boltzmann_inversion/boltzmann_inversion.py") + + +@skipIfMissingFeatures +class Tutorial(ut.TestCase): + system = tutorial.system + + def test_boltzmann_inversion(self): + rdf_exp = tutorial.rdf + rdf_imp = tutorial.rdf_dh + np.testing.assert_allclose(rdf_exp, rdf_imp, rtol=0.1, atol=0.18) + + +if __name__ == "__main__": + ut.main()