diff --git a/.gitlab-ci.yml b/.gitlab-ci.yml
index ee94ce79c16..e1390b27a79 100644
--- a/.gitlab-ci.yml
+++ b/.gitlab-ci.yml
@@ -1,4 +1,4 @@
-image: docker.pkg.github.com/espressomd/docker/ubuntu-20.04:8a7d258889d87dcc9c96a200784b30b532d45b33
+image: docker.pkg.github.com/espressomd/docker/ubuntu-20.04:e583d4b2eb8eedd10068957f952bd67008475ee5
stages:
- prepare
@@ -114,7 +114,7 @@ no_rotation:
ubuntu:wo-dependencies:
<<: *global_job_definition
stage: build
- image: docker.pkg.github.com/espressomd/docker/ubuntu-wo-dependencies:fc7628d32de0fce605976ba9edebe7eff186e618
+ image: docker.pkg.github.com/espressomd/docker/ubuntu-wo-dependencies:e583d4b2eb8eedd10068957f952bd67008475ee5
script:
- export myconfig=maxset with_cuda=false make_check_unit_tests=false make_check_python=false
- bash maintainer/CI/build_cmake.sh
@@ -199,6 +199,24 @@ clang-sanitizer:
- linux
- cuda
+cuda11-maxset:
+ <<: *global_job_definition
+ stage: build
+ image: docker.pkg.github.com/espressomd/docker/cuda:e583d4b2eb8eedd10068957f952bd67008475ee5
+ variables:
+ CC: 'gcc-9'
+ CXX: 'g++-9'
+ script:
+ - export myconfig=maxset with_cuda=true with_coverage=false test_timeout=900 srcdir=${CI_PROJECT_DIR}
+ - export with_scafacos=true with_stokesian_dynamics=true
+ - bash maintainer/CI/build_cmake.sh
+ tags:
+ - docker
+ - linux
+ - cuda
+ only:
+ - schedules
+
cuda10-maxset:
<<: *global_job_definition
stage: build
@@ -230,11 +248,6 @@ cuda9-maxset:
- export myconfig=maxset with_cuda=true with_coverage=true test_timeout=900 srcdir=${CI_PROJECT_DIR}
- export with_scafacos=true with_stokesian_dynamics=true
- bash maintainer/CI/build_cmake.sh
- artifacts:
- paths:
- - build/
- expire_in: 1 week
- when: on_success
tags:
- docker
- linux
@@ -354,7 +367,7 @@ empty:
rocm-maxset:
<<: *global_job_definition
stage: build
- image: docker.pkg.github.com/espressomd/docker/rocm:d496478230db4e5c286680e3bdc1621af1fccffc
+ image: docker.pkg.github.com/espressomd/docker/rocm:cb775a82f49a7f6d4a806f17358caec797f48f73
script:
- export myconfig=maxset with_cuda=true with_cuda_compiler=hip
- export with_stokesian_dynamics=false
@@ -457,7 +470,7 @@ check_cuda_maxset_no_gpu:
<<: *global_job_definition
stage: additional_checks
when: on_success
- needs:
+ needs:
- cuda10-maxset
script:
- export CUDA_VISIBLE_DEVICES=""
@@ -470,9 +483,9 @@ check_cuda_maxset_no_gpu:
check_with_odd_no_of_processors:
<<: *global_job_definition
- stage: additional_checks
+ stage: additional_checks
when: on_success
- needs:
+ needs:
- cuda10-maxset
script:
- cd ${CI_PROJECT_DIR}/build
@@ -507,7 +520,7 @@ deploy_sphinx_documentation:
dependencies:
- check_sphinx
script:
- - cd ${CI_PROJECT_DIR}/build/doc/sphinx/html &&
+ - cd ${CI_PROJECT_DIR}/build/doc/sphinx/html &&
rsync -avz --delete -e "ssh -i ${HOME}/.ssh/espresso_rsa" ./ espresso@elk.icp.uni-stuttgart.de:/home/espresso/public_html/html/doc
deploy_doxygen_documentation:
diff --git a/CMakeLists.txt b/CMakeLists.txt
index 8a322bde769..6e887417044 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -18,7 +18,7 @@
# along with this program. If not, see .
#
-cmake_minimum_required(VERSION 3.10)
+cmake_minimum_required(VERSION 3.11)
message(STATUS "CMake version: ${CMAKE_VERSION}")
if(POLICY CMP0076)
# make target_sources() convert relative paths to absolute
@@ -166,7 +166,7 @@ if(WITH_CUDA)
endif()
endif(WITH_CUDA)
-find_package(PythonInterp 3.5 REQUIRED)
+find_package(PythonInterp 3.6 REQUIRED)
if(WITH_PYTHON)
find_package(Cython 0.26 REQUIRED)
@@ -278,7 +278,6 @@ endif(WITH_STOKESIAN_DYNAMICS)
if(WITH_STOKESIAN_DYNAMICS)
set(CMAKE_INSTALL_LIBDIR
"${CMAKE_INSTALL_PREFIX}/${PYTHON_INSTDIR}/espressomd")
- cmake_minimum_required(VERSION 3.11)
include(FetchContent)
FetchContent_Declare(
stokesian_dynamics
@@ -485,6 +484,9 @@ if(WITH_TESTS)
if(WITH_PYTHON)
add_subdirectory(testsuite)
endif(WITH_PYTHON)
+ set(CTEST_ARGS ""
+ CACHE STRING
+ "Extra arguments to give to ctest calls (separated by semicolons)")
endif(WITH_TESTS)
if(WITH_BENCHMARKS)
diff --git a/cmake/FindCUDACompilerHIP.cmake b/cmake/FindCUDACompilerHIP.cmake
index 3440557f8bd..b6e708b169a 100644
--- a/cmake/FindCUDACompilerHIP.cmake
+++ b/cmake/FindCUDACompilerHIP.cmake
@@ -27,21 +27,26 @@ find_package(HIP ${CUDACompilerHIP_FIND_VERSION} MODULE REQUIRED)
# patch HCC_PATH environment variable and reload HIP
if(HIP_VERSION VERSION_LESS 3.1)
set(HCC_PATH "${HIP_ROOT_DIR}")
-else()
+elseif(HIP_VERSION VERSION_LESS 3.5)
set(HCC_PATH "${ROCM_HOME}/hcc")
+else()
+ set(HIP_HIPCC_CMAKE_LINKER_HELPER "${HIP_HIPCC_EXECUTABLE}")
+ unset(HCC_PATH)
endif()
find_package(HIP ${CUDACompilerHIP_FIND_VERSION} MODULE REQUIRED)
set(CUDA 1)
set(HIP 1)
-list(APPEND HIP_HCC_FLAGS
+list(APPEND HIP_HIPCC_FLAGS
-std=c++${CMAKE_CUDA_STANDARD} -pedantic -Wall -Wextra
-Wno-sign-compare -Wno-unused-function -Wno-unused-variable
-Wno-unused-parameter -Wno-missing-braces -Wno-gnu-anonymous-struct
-Wno-nested-anon-types -Wno-gnu-zero-variadic-macro-arguments
-Wno-c99-designator -Wno-macro-redefined -Wno-duplicate-decl-specifier
$<$:-Wno-deprecated-copy>
+ $<$:-Wno-c++17-extensions>
+ $<$:-Wno-unused-command-line-argument>
$<$:-Werror>)
list(APPEND HIP_HIPCC_FLAGS_DEBUG -g)
diff --git a/cmake/FindCUDACompilerNVCC.cmake b/cmake/FindCUDACompilerNVCC.cmake
index 53a63d37f29..1a77224d765 100644
--- a/cmake/FindCUDACompilerNVCC.cmake
+++ b/cmake/FindCUDACompilerNVCC.cmake
@@ -50,8 +50,11 @@ list(APPEND CUDA_NVCC_FLAGS_MINSIZEREL -O2 -Xptxas=-O2 -Xcompiler=-Os -DNDEBUG)
list(APPEND CUDA_NVCC_FLAGS_RELWITHDEBINFO -O2 -g -Xptxas=-O2 -Xcompiler=-O2,-g -DNDEBUG)
list(APPEND CUDA_NVCC_FLAGS_COVERAGE -O3 -g -Xptxas=-O3 -Xcompiler=-Og,-g)
list(APPEND CUDA_NVCC_FLAGS_RELWITHASSERT -O3 -g -Xptxas=-O3 -Xcompiler=-O3,-g)
+if(CMAKE_CUDA_COMPILER_VERSION VERSION_LESS 11)
+ list(APPEND CUDA_NVCC_FLAGS -gencode=arch=compute_30,code=sm_30)
+endif()
list(APPEND CUDA_NVCC_FLAGS
- -gencode=arch=compute_30,code=sm_30 -gencode=arch=compute_52,code=sm_52
+ -gencode=arch=compute_52,code=sm_52
-gencode=arch=compute_52,code=compute_52 -std=c++${CMAKE_CUDA_STANDARD}
$<$:-Xcompiler=-Werror;-Xptxas=-Werror>
$<$:-Xcompiler=-isysroot;-Xcompiler=${CMAKE_OSX_SYSROOT}>)
@@ -60,6 +63,7 @@ function(find_gpu_library)
cmake_parse_arguments(LIBRARY "REQUIRED" "NAMES;VARNAME" "" ${ARGN})
list(APPEND LIBRARY_PATHS
${CUDA_TOOLKIT_ROOT_DIR}/lib64 ${CUDA_TOOLKIT_ROOT_DIR}/lib
+ ${CUDA_TOOLKIT_ROOT_DIR}/compat
/usr/local/nvidia/lib /usr/lib/x86_64-linux-gnu)
if(LIBRARY_REQUIRED)
find_library(${LIBRARY_VARNAME} NAMES ${LIBRARY_NAMES} PATHS ${LIBRARY_PATHS} NO_DEFAULT_PATH REQUIRED)
diff --git a/cmake/FindFFTW3.cmake b/cmake/FindFFTW3.cmake
index 57c914cfdce..54ede3fe470 100644
--- a/cmake/FindFFTW3.cmake
+++ b/cmake/FindFFTW3.cmake
@@ -28,7 +28,6 @@ if(FFTW3_INCLUDE_DIR)
endif(FFTW3_INCLUDE_DIR)
find_path(FFTW3_INCLUDE_DIR fftw3.h)
-
find_library(FFTW3_LIBRARIES NAMES fftw3)
# handle the QUIETLY and REQUIRED arguments and set FFTW_FOUND to TRUE if all
@@ -38,3 +37,9 @@ find_package_handle_standard_args(FFTW3 DEFAULT_MSG FFTW3_LIBRARIES
FFTW3_INCLUDE_DIR)
mark_as_advanced(FFTW3_LIBRARIES FFTW3_INCLUDE_DIR)
+
+if(FFTW3_FOUND AND NOT TARGET FFTW3::FFTW3)
+ add_library(FFTW3::FFTW3 INTERFACE IMPORTED)
+ target_include_directories(FFTW3::FFTW3 INTERFACE "${FFTW3_INCLUDE_DIR}")
+ target_link_libraries(FFTW3::FFTW3 INTERFACE "${FFTW3_LIBRARIES}")
+endif()
diff --git a/doc/sphinx/advanced_methods.rst b/doc/sphinx/advanced_methods.rst
index 3ad5c8805e2..6ee68150b63 100644
--- a/doc/sphinx/advanced_methods.rst
+++ b/doc/sphinx/advanced_methods.rst
@@ -28,7 +28,7 @@ Several modes are available for different types of binding.
import espressomd
from espressomd.interactions import HarmonicBond
- system = espressomd.System()
+ system = espressomd.System(box_l=[1, 1, 1])
bond_centers = HarmonicBond(k=1000, r_0=)
system.bonded_inter.add(bond_centers)
system.collision_detection.set_params(mode="bind_centers", distance=,
diff --git a/doc/sphinx/constraints.rst b/doc/sphinx/constraints.rst
index 79c649e8af0..f1c94fec944 100644
--- a/doc/sphinx/constraints.rst
+++ b/doc/sphinx/constraints.rst
@@ -167,7 +167,7 @@ Available shapes
Python syntax::
import espressomd from espressomd.shapes import
- system = espressomd.System()
+ system = espressomd.System(box_l=[1, 1, 1])
```` can be any of the available shapes.
diff --git a/doc/sphinx/electrostatics.rst b/doc/sphinx/electrostatics.rst
index d5f1ae1bac3..e31caff4bd7 100644
--- a/doc/sphinx/electrostatics.rst
+++ b/doc/sphinx/electrostatics.rst
@@ -35,20 +35,22 @@ Note that using the electrostatic interaction also requires assigning charges to
the particles via the particle property
:py:attr:`espressomd.particle_data.ParticleHandle.q`.
-This example shows the general usage of an electrostatic method ````.
-All of them need the Bjerrum length and a set of other required parameters.
-First, an instance of the solver is created and only after adding it to the actors
-list, it is activated. Internally the method calls a tuning routine on
-activation to achieve the given accuracy::
+All solvers need a prefactor and a set of other required parameters.
+This example shows the general usage of the electrostatic method ``P3M``.
+An instance of the solver is created and added to the actors list, at which
+point it will be automatically activated. This activation will internally
+call a tuning function to achieve the requested accuracy::
import espressomd
- from espressomd import electrostatics
+ import espressomd.electrostatics
- system = espressomd.System()
- solver = electrostatics.(prefactor=C, )
+ system = espressomd.System(box_l=[10, 10, 10])
+ system.time_step = 0.01
+ system.part.add(pos=[[0, 0, 0], [1, 1, 1]], q=[-1, 1])
+ solver = espressomd.electrostatics.P3M(prefactor=2, accuracy=1e-3)
system.actors.add(solver)
-where the prefactor :math:`C` is defined as in Eqn. :eq:`coulomb_prefactor`
+where the prefactor is defined as :math:`C` in Eqn. :eq:`coulomb_prefactor`.
.. _Coulomb P3M:
@@ -60,9 +62,10 @@ Coulomb P3M
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``.
-P3M requires full periodicity (1 1 1). Make sure that you know the relevance of the
-P3M parameters before using P3M! If you are not sure, read the following
-references:
+P3M requires full periodicity (1 1 1). When using a non-metallic dielectric
+constant (``epsilon != 0.0``), the box must be cubic.
+Make sure that you know the relevance of the P3M parameters before using P3M!
+If you are not sure, read the following references:
:cite:`ewald21,hockney88,kolafa92,deserno98a,deserno98b,deserno00,deserno00a,cerda08d`.
.. _Tuning Coulomb P3M:
diff --git a/doc/sphinx/installation.rst b/doc/sphinx/installation.rst
index fc701b081de..3d1d44d0afc 100644
--- a/doc/sphinx/installation.rst
+++ b/doc/sphinx/installation.rst
@@ -89,7 +89,11 @@ are required:
.. code-block:: bash
sudo apt install python3-matplotlib python3-scipy ipython3 jupyter-notebook
- sudo pip3 install 'pint>=0.9'
+ pip3 install --user 'pint>=0.9' 'jupyter_contrib_nbextensions==0.5.1' \
+ 'sphinx>=1.6.7,!=2.1.0,!=3.0.0' 'sphinxcontrib-bibtex>=0.3.5'
+ jupyter contrib nbextension install --user
+ jupyter nbextension enable rubberband/main
+ jupyter nbextension enable exercise2/main
Nvidia GPU acceleration
"""""""""""""""""""""""
@@ -829,10 +833,10 @@ to actually write a simulation script for |es|.
Running an interactive notebook
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
-Running the Jupyter interpreter requires using the ``ipypresso`` script, which
+Running a Jupyter session requires using the ``ipypresso`` script, which
is also located in the build directory (its name comes from the IPython
interpreter, today known as Jupyter). To run the tutorials, you will need
-to start the Jupyter interpreter in notebook mode:
+to start a Jupyter session:
.. code-block:: bash
@@ -858,11 +862,16 @@ will exit the Python interpreter and Jupyter will notify you that the current
Python kernel stopped. If a cell takes too long to execute, you may interrupt
it with the stop button.
-To close the Jupyter notebook, go to the terminal where it was started and use
+Solutions cells are created using the ``exercise2`` plugin from nbextensions.
+To prevent solution code cells from running when clicking on "Run All", these
+code cells need to be converted to Markdown cells and fenced with `````python``
+and ```````.
+
+To close the Jupyter session, go to the terminal where it was started and use
the keyboard shortcut Ctrl+C twice.
-When starting the Jupyter interpreter in notebook mode, you may see the
-following warning in the terminal:
+When starting a Jupyter session, you may see the following warning in the
+terminal:
.. code-block:: none
diff --git a/doc/sphinx/inter_non-bonded.rst b/doc/sphinx/inter_non-bonded.rst
index 9dff6b5cbbc..51061e3bb78 100644
--- a/doc/sphinx/inter_non-bonded.rst
+++ b/doc/sphinx/inter_non-bonded.rst
@@ -620,7 +620,7 @@ that the Thole correction acts between all dipoles, intra- and intermolecular.
Again, the accuracy is related to the P3M accuracy and the split between
short-range and long-range electrostatics interaction. It is configured by::
- system = espressomd.System()
+ system = espressomd.System(box_l=[1, 1, 1])
system.non_bonded_inter[type_1,type_2].thole.set_params(scaling_coeff=, q1q2=)
with parameters:
diff --git a/doc/sphinx/io.rst b/doc/sphinx/io.rst
index 3e21710c60b..dccb32fd0d2 100644
--- a/doc/sphinx/io.rst
+++ b/doc/sphinx/io.rst
@@ -211,7 +211,7 @@ capabilities. The usage is quite simple:
.. code:: python
from espressomd.io.mppiio import mpiio
- system = espressomd.System()
+ system = espressomd.System(box_l=[1, 1, 1])
# ... add particles here
mpiio.write("/tmp/mydata", positions=True, velocities=True, types=True, bonds=True)
diff --git a/doc/sphinx/magnetostatics.rst b/doc/sphinx/magnetostatics.rst
index 417a182d030..52fc41df91d 100644
--- a/doc/sphinx/magnetostatics.rst
+++ b/doc/sphinx/magnetostatics.rst
@@ -27,9 +27,15 @@ relative permittivity of the background material, respectively.
Magnetostatic interactions are activated via the actor framework::
- from espressomd.magnetostatics import DipolarDirectSumCpu
+ import espressomd
+ import espressomd.magnetostatics
- direct_sum = DipolarDirectSumCpu(prefactor=1)
+ system = espressomd.System(box_l=[10, 10, 10])
+ system.time_step = 0.01
+ system.part.add(pos=[[0, 0, 0], [1, 1, 1]],
+ rotation=2 * [(1, 1, 1)], dip=2 * [(1, 0, 0)])
+
+ direct_sum = espressomd.magnetostatics.DipolarDirectSumCpu(prefactor=1)
system.actors.add(direct_sum)
# ...
system.actors.remove(direct_sum)
@@ -49,7 +55,8 @@ Dipolar P3M
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
-P3M! If you are not sure, read the following references
+P3M! If you are not sure, read the following references:
+:cite:`ewald21,hockney88,kolafa92,deserno98a,deserno98b,deserno00,deserno00a,cerda08d`.
Note that dipolar P3M does not work with non-cubic boxes.
diff --git a/doc/sphinx/particles.rst b/doc/sphinx/particles.rst
index fbbd0dad801..cb1b3b00463 100644
--- a/doc/sphinx/particles.rst
+++ b/doc/sphinx/particles.rst
@@ -25,7 +25,7 @@ In order to add particles to the system, call
:meth:`espressomd.particle_data.ParticleList.add`::
import espressomd
- system = espressomd.System()
+ system = espressomd.System(box_l=[1, 1, 1])
system.part.add(pos=[1.0, 1.0, 1.0], id=0, type=0)
This command adds a single particle to the system with properties given
@@ -308,7 +308,7 @@ To switch the active scheme, the attribute :attr:`espressomd.system.System.virtu
import espressomd
from espressomd.virtual_sites import VirtualSitesOff, VirtualSitesRelative
- system = espressomd.System()
+ system = espressomd.System(box_l=[1, 1, 1])
system.virtual_sites = VirtualSitesRelative(have_quaternion=False)
# or
system.virtual_sites = VirtualSitesOff()
@@ -439,7 +439,7 @@ Particle ids can be stored in a map for each
individual type::
import espressomd
- system = espressomd.System()
+ system = espressomd.System(box_l=[1, 1, 1])
system.setup_type_map([_type])
system.number_of_particles(_type)
@@ -454,7 +454,7 @@ The keyword
particles which have the given type. For counting the number of particles of a given type you could also use :meth:`espressomd.particle_data.ParticleList.select` ::
import espressomd
- system = espressomd.System()
+ system = espressomd.System(box_l=[1, 1, 1])
...
number_of_particles = len(system.part.select(type=type))
@@ -483,7 +483,7 @@ Langevin swimmers
import espressomd
- system = espressomd.System()
+ system = espressomd.System(box_l=[1, 1, 1])
system.part.add(id=0, pos=[1, 0, 0], swimming={'f_swim': 0.03})
@@ -510,7 +510,7 @@ Lattice-Boltzmann (LB) swimmers
import espressomd
- system = espressomd.System()
+ system = espressomd.System(box_l=[1, 1, 1])
system.part.add(id=1, pos=[2, 0, 0], rotation=[1, 1, 1], swimming={
'f_swim': 0.01, 'mode': 'pusher', 'dipole_length': 2.0})
diff --git a/doc/sphinx/running.rst b/doc/sphinx/running.rst
index aa2930dbe3f..3ad73566625 100644
--- a/doc/sphinx/running.rst
+++ b/doc/sphinx/running.rst
@@ -105,7 +105,7 @@ A code snippet would look like::
import espressomd
- system = espressomd.System()
+ system = espressomd.System(box_l=[1, 1, 1])
system.thermostat.set_npt(kT=1.0, gamma0=1.0, gammav=1.0, seed=42)
system.integrator.set_isotropic_npt(ext_pressure=1.0, piston=1.0)
@@ -204,7 +204,7 @@ In this way, just compiling in the ``ROTATION`` feature no longer changes the ph
The rotation of a particle is controlled via the :attr:`espressomd.particle_data.ParticleHandle.rotation` property. E.g., the following code adds a particle with rotation enabled on the x axis::
import espressomd
- system = espressomd.System()
+ system = espressomd.System(box_l=[1, 1, 1])
system.part.add(pos=(0, 0, 0), rotation=(1, 0, 0))
Notes:
diff --git a/doc/sphinx/system_manipulation.rst b/doc/sphinx/system_manipulation.rst
index 4188de7e1f1..d399a59b2be 100644
--- a/doc/sphinx/system_manipulation.rst
+++ b/doc/sphinx/system_manipulation.rst
@@ -96,7 +96,7 @@ Galilei Transform and Particle Velocity Manipulation
The following class :class:`espressomd.galilei.GalileiTransform` may be useful
in affecting the velocity of the system. ::
- system = espressomd.System()
+ system = espressomd.System(box_l=[1, 1, 1])
gt = system.galilei
* Particle motion and rotation
diff --git a/doc/sphinx/system_setup.rst b/doc/sphinx/system_setup.rst
index bafbafdd0df..bd0b6992e8c 100644
--- a/doc/sphinx/system_setup.rst
+++ b/doc/sphinx/system_setup.rst
@@ -71,17 +71,17 @@ Accessing module states
Some variables like or are no longer directly available as attributes.
In these cases they can be easily derived from the corresponding Python
-objects like
+objects like::
-``n_part = len(espressomd.System().part[:].pos)``
+ n_part = len(system.part[:].pos)
or by calling the corresponding ``get_state()`` methods like::
- temperature = espressomd.System().thermostat.get_state()[0]['kT']
+ temperature = system.thermostat.get_state()[0]['kT']
- gamma = espressomd.System().thermostat.get_state()[0]['gamma']
+ gamma = system.thermostat.get_state()[0]['gamma']
- gamma_rot = espressomd.System().thermostat.get_state()[0]['gamma_rotation']
+ gamma_rot = system.thermostat.get_state()[0]['gamma_rotation']
.. _Cellsystems:
@@ -112,7 +112,7 @@ The properties of the cell system can be accessed by
(float) Skin for the Verlet list. This value has to be set, otherwise the simulation will not start.
-Details about the cell system can be obtained by :meth:`espressomd.System().cell_system.get_state() `:
+Details about the cell system can be obtained by :meth:`espressomd.system.System.cell_system.get_state() `:
* ``cell_grid`` Dimension of the inner cell grid.
* ``cell_size`` Box-length of a cell.
@@ -132,7 +132,7 @@ selects the domain decomposition cell scheme, using Verlet lists
for the calculation of the interactions. If you specify ``use_verlet_lists=False``, only the
domain decomposition is used, but not the Verlet lists. ::
- system = espressomd.System()
+ system = espressomd.System(box_l=[1, 1, 1])
system.cell_system.set_domain_decomposition(use_verlet_lists=True)
@@ -161,7 +161,7 @@ particles, giving an unfavorable computation time scaling of
interaction in the cell model require the calculation of all pair
interactions. ::
- system = espressomd.System()
+ system = espressomd.System(box_l=[1, 1, 1])
system.cell_system.set_n_square()
@@ -221,10 +221,8 @@ class :class:`espressomd.thermostat.Thermostat` has to be invoked.
Best explained in an example::
import espressomd
- system = espressomd.System()
- therm = system.Thermostat()
-
- therm.set_langevin(kT=1.0, gamma=1.0, seed=41)
+ system = espressomd.System(box_l=[1, 1, 1])
+ system.thermostat.set_langevin(kT=1.0, gamma=1.0, seed=41)
As explained before the temperature is set as thermal energy :math:`k_\mathrm{B} T`.
@@ -383,7 +381,7 @@ For example::
import espressomd
- system = espressomd.System()
+ system = espressomd.System(box_l=[1, 1, 1])
system.thermostat.set_npt(kT=1.0, gamma0=1.0, gammav=1.0, seed=41)
system.integrator.set_isotropic_npt(ext_pressure=1.0, piston=1.0)
@@ -416,7 +414,7 @@ The system integrator should be also changed.
Best explained in an example::
import espressomd
- system = espressomd.System()
+ system = espressomd.System(box_l=[1, 1, 1])
system.thermostat.set_brownian(kT=1.0, gamma=1.0, seed=41)
system.integrator.set_brownian_dynamics()
@@ -505,7 +503,7 @@ example ``lbfluid``. If you do not choose the GPU manually before that,
CUDA internally chooses one, which is normally the most powerful GPU
available, but load-independent. ::
- system = espressomd.System()
+ system = espressomd.System(box_l=[1, 1, 1])
dev = system.cuda_init_handle.device
system.cuda_init_handle.device = dev
@@ -539,7 +537,7 @@ List available CUDA devices
If you want to list available CUDA devices
you should access :attr:`espressomd.cuda_init.CudaInitHandle.device_list`, e.g., ::
- system = espressomd.System()
+ system = espressomd.System(box_l=[1, 1, 1])
print(system.cuda_init_handle.device_list)
@@ -555,7 +553,7 @@ When you start ``pypresso`` your first GPU should be selected.
If you wanted to use the second GPU, this can be done
by setting :attr:`espressomd.cuda_init.CudaInitHandle.device` as follows::
- system = espressomd.System()
+ system = espressomd.System(box_l=[1, 1, 1])
system.cuda_init_handle.device = 1
diff --git a/doc/sphinx/visualization.rst b/doc/sphinx/visualization.rst
index b30f8c61952..0fab958d9ca 100644
--- a/doc/sphinx/visualization.rst
+++ b/doc/sphinx/visualization.rst
@@ -42,7 +42,7 @@ window with ``start()``. See the following minimal code example::
from espressomd import visualization
from threading import Thread
- system = espressomd.System()
+ system = espressomd.System(box_l=[1, 1, 1])
system.cell_system.skin = 0.4
system.time_step = 0.01
system.box_l = [10, 10, 10]
diff --git a/doc/tutorials/04-lattice_boltzmann/04-lattice_boltzmann_part4.ipynb b/doc/tutorials/04-lattice_boltzmann/04-lattice_boltzmann_part4.ipynb
index 37eb1ace58e..460e79ac282 100644
--- a/doc/tutorials/04-lattice_boltzmann/04-lattice_boltzmann_part4.ipynb
+++ b/doc/tutorials/04-lattice_boltzmann/04-lattice_boltzmann_part4.ipynb
@@ -38,9 +38,7 @@
"\n",
"We will simulate a planar Poiseuille flow using a square box, two walls\n",
"with normal vectors $\\left(\\pm 1, 0, 0 \\right)$, and an external force density\n",
- "applied to every node.\n",
- "\n",
- "Use the data to fit a parabolic function. Can you confirm the analytic solution?"
+ "applied to every node."
]
},
{
@@ -52,6 +50,7 @@
"import logging\n",
"import sys\n",
"\n",
+ "import matplotlib.pyplot as plt\n",
"import numpy as np\n",
"\n",
"import espressomd\n",
@@ -102,9 +101,54 @@
},
{
"cell_type": "markdown",
- "metadata": {},
+ "metadata": {
+ "solution2": "hidden",
+ "solution2_first": true
+ },
"source": [
- "The solution is available at /doc/tutorials/04-lattice_boltzmann/scripts/04-lattice_boltzmann_part4_solution.py"
+ "Use the data to fit a parabolic function. Can you confirm the analytic solution?"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {
+ "solution2": "hidden"
+ },
+ "source": [
+ "```python\n",
+ "# Extract fluid velocity along the x-axis\n",
+ "fluid_velocities = np.zeros((lbf.shape[0], 2))\n",
+ "for x in range(lbf.shape[0]):\n",
+ " # Average over the node in y direction\n",
+ " v_tmp = np.zeros(lbf.shape[1])\n",
+ " for y in range(lbf.shape[1]):\n",
+ " v_tmp[y] = lbf[x, y, 0].velocity[1]\n",
+ " fluid_velocities[x, 0] = (x + 0.5) * AGRID\n",
+ " fluid_velocities[x, 1] = np.average(v_tmp)\n",
+ "\n",
+ "def poiseuille_flow(x, force_density, dynamic_viscosity, height):\n",
+ " return force_density / (2.0 * dynamic_viscosity) * \\\n",
+ " (height**2.0 / 4.0 - x**2.0)\n",
+ "\n",
+ "# Note that the LB viscosity is not the dynamic viscosity but the\n",
+ "# kinematic viscosity (mu=LB_viscosity * density)\n",
+ "x_values = np.linspace(0.0, BOX_L, lbf.shape[0])\n",
+ "HEIGHT = BOX_L - 2.0 * AGRID\n",
+ "# analytical curve\n",
+ "y_values = poiseuille_flow(x_values - (HEIGHT / 2.0 + AGRID), FORCE_DENSITY[1],\n",
+ " VISCOSITY * DENSITY, HEIGHT)\n",
+ "# velocity is zero outside the walls\n",
+ "y_values[np.nonzero(x_values < WALL_OFFSET)] = 0.0\n",
+ "y_values[np.nonzero(x_values > BOX_L - WALL_OFFSET)] = 0.0\n",
+ "\n",
+ "fig1 = plt.figure(num=None, figsize=(10, 6), dpi=80, facecolor='w', edgecolor='k')\n",
+ "plt.plot(x_values, y_values, '-', linewidth=2, label='analytical')\n",
+ "plt.plot(fluid_velocities[:, 0], fluid_velocities[:, 1], 'o', label='simulation')\n",
+ "plt.xlabel('Position on the $x$-axis', fontsize=16)\n",
+ "plt.ylabel('Fluid velocity in $y$-direction', fontsize=16)\n",
+ "plt.legend()\n",
+ "plt.show()\n",
+ "```"
]
},
{
diff --git a/doc/tutorials/04-lattice_boltzmann/CMakeLists.txt b/doc/tutorials/04-lattice_boltzmann/CMakeLists.txt
index 56d99b48bd8..85f800c7f72 100644
--- a/doc/tutorials/04-lattice_boltzmann/CMakeLists.txt
+++ b/doc/tutorials/04-lattice_boltzmann/CMakeLists.txt
@@ -3,20 +3,7 @@ configure_tutorial_target(
04-lattice_boltzmann_part2.ipynb 04-lattice_boltzmann_part3.ipynb
04-lattice_boltzmann_part4.ipynb figures/latticeboltzmann-grid.png
figures/latticeboltzmann-momentumexchange.png
- scripts/04-lattice_boltzmann_part3_solution.py
- scripts/04-lattice_boltzmann_part4_solution.py)
-
-add_custom_command(
- OUTPUT
- "${CMAKE_CURRENT_BINARY_DIR}/scripts/04-lattice_boltzmann_part4_solution_cut.py"
- DEPENDS
- "${CMAKE_CURRENT_SOURCE_DIR}/scripts/04-lattice_boltzmann_part4_solution.py"
- DEPENDS
- "${CMAKE_CURRENT_SOURCE_DIR}/scripts/04-lattice_boltzmann_part4_solution_cut.cmake"
- COMMAND
- ${CMAKE_COMMAND} -P
- "${CMAKE_CURRENT_SOURCE_DIR}/scripts/04-lattice_boltzmann_part4_solution_cut.cmake"
-)
+ scripts/04-lattice_boltzmann_part3_solution.py)
nb_export(TARGET tutorial_04 SUFFIX "1" FILE "04-lattice_boltzmann_part1.ipynb"
HTML_RUN)
@@ -24,14 +11,5 @@ nb_export(TARGET tutorial_04 SUFFIX "2" FILE "04-lattice_boltzmann_part2.ipynb"
HTML_RUN)
nb_export(TARGET tutorial_04 SUFFIX "3" FILE "04-lattice_boltzmann_part3.ipynb"
HTML_RUN)
-nb_export(
- TARGET
- tutorial_04
- SUFFIX
- "4"
- FILE
- "04-lattice_boltzmann_part4.ipynb"
- HTML_RUN
- ADD_SCRIPTS
- "${CMAKE_CURRENT_BINARY_DIR}/scripts/04-lattice_boltzmann_part4_solution_cut.py"
-)
+nb_export(TARGET tutorial_04 SUFFIX "4" FILE "04-lattice_boltzmann_part4.ipynb"
+ HTML_RUN)
diff --git a/doc/tutorials/04-lattice_boltzmann/scripts/04-lattice_boltzmann_part4_solution.py b/doc/tutorials/04-lattice_boltzmann/scripts/04-lattice_boltzmann_part4_solution.py
deleted file mode 100644
index b2c946bec33..00000000000
--- a/doc/tutorials/04-lattice_boltzmann/scripts/04-lattice_boltzmann_part4_solution.py
+++ /dev/null
@@ -1,86 +0,0 @@
-# Copyright (C) 2019 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 matplotlib.pyplot as plt
-import numpy as np
-
-import espressomd
-espressomd.assert_features(['LB_BOUNDARIES_GPU'])
-import espressomd.lb
-import espressomd.lbboundaries
-import espressomd.shapes
-
-# System setup
-BOX_L = 16.0
-AGRID = 0.5
-VISCOSITY = 2.0
-FORCE_DENSITY = [0.0, 0.001, 0.0]
-DENSITY = 1.5
-TIME_STEP = 0.01
-system = espressomd.System(box_l=[BOX_L] * 3)
-system.time_step = TIME_STEP
-system.cell_system.skin = 0.4
-
-lbf = espressomd.lb.LBFluidGPU(
- agrid=AGRID, dens=DENSITY, visc=VISCOSITY, tau=TIME_STEP,
- ext_force_density=FORCE_DENSITY)
-system.actors.add(lbf)
-
-# Setup boundaries
-WALL_OFFSET = AGRID
-top_wall = espressomd.lbboundaries.LBBoundary(
- shape=espressomd.shapes.Wall(normal=[1, 0, 0], dist=WALL_OFFSET))
-bottom_wall = espressomd.lbboundaries.LBBoundary(
- shape=espressomd.shapes.Wall(normal=[-1, 0, 0], dist=-(BOX_L - WALL_OFFSET)))
-
-system.lbboundaries.add(top_wall)
-system.lbboundaries.add(bottom_wall)
-
-# Iterate until the flow profile converges (5000 LB updates)
-system.integrator.run(5000)
-
-# Extract fluid velocity along the x-axis
-fluid_velocities = np.zeros((lbf.shape[0], 2))
-for x in range(lbf.shape[0]):
- # Average over the node in y direction
- v_tmp = np.zeros(lbf.shape[1])
- for y in range(lbf.shape[1]):
- v_tmp[y] = lbf[x, y, 0].velocity[1]
- fluid_velocities[x, 0] = (x + 0.5) * AGRID
- fluid_velocities[x, 1] = np.average(v_tmp)
-
-
-def poiseuille_flow(x, force_density, dynamic_viscosity, height):
- return force_density / (2.0 * dynamic_viscosity) * \
- (height**2.0 / 4.0 - x**2.0)
-
-
-# Note that the LB viscosity is not the dynamic viscosity but the
-# kinematic viscosity (mu=LB_viscosity * density)
-x_values = np.linspace(0.0, BOX_L, lbf.shape[0])
-HEIGHT = BOX_L - 2.0 * AGRID
-# analytical curve
-y_values = poiseuille_flow(x_values - (HEIGHT / 2.0 + AGRID), FORCE_DENSITY[1],
- VISCOSITY * DENSITY, HEIGHT)
-# velocity is zero outside the walls
-y_values[np.nonzero(x_values < WALL_OFFSET)] = 0.0
-y_values[np.nonzero(x_values > BOX_L - WALL_OFFSET)] = 0.0
-
-plt.plot(x_values, y_values, 'o-', label='analytical')
-plt.plot(fluid_velocities[:, 0], fluid_velocities[:, 1], label='simulation')
-plt.legend()
-plt.show()
diff --git a/doc/tutorials/04-lattice_boltzmann/scripts/04-lattice_boltzmann_part4_solution_cut.cmake b/doc/tutorials/04-lattice_boltzmann/scripts/04-lattice_boltzmann_part4_solution_cut.cmake
deleted file mode 100644
index 3085c2e2a65..00000000000
--- a/doc/tutorials/04-lattice_boltzmann/scripts/04-lattice_boltzmann_part4_solution_cut.cmake
+++ /dev/null
@@ -1,14 +0,0 @@
-# This code splits 04-lattice_boltzmann_part4_solution.py such that it can be
-# added to the tutorial without having two espressomd.System class instances
-file(
- READ
- "${CMAKE_CURRENT_SOURCE_DIR}/scripts/04-lattice_boltzmann_part4_solution.py"
- script_full)
-string(FIND "${script_full}" "# Extract fluid velocity along the x-axis"
- cut_position)
-string(SUBSTRING "${script_full}" ${cut_position} -1 script_cut)
-string(PREPEND script_cut "import matplotlib.pyplot as plt\n")
-file(
- WRITE
- "${CMAKE_CURRENT_BINARY_DIR}/scripts/04-lattice_boltzmann_part4_solution_cut.py"
- "${script_cut}")
diff --git a/doc/tutorials/CMakeLists.txt b/doc/tutorials/CMakeLists.txt
index 826cdf0d780..831b78ea6f2 100644
--- a/doc/tutorials/CMakeLists.txt
+++ b/doc/tutorials/CMakeLists.txt
@@ -53,29 +53,39 @@ function(NB_EXPORT)
set(PY_FILE "${NB_FILE_BASE}.py")
if(${NB_EXPORT_HTML_RUN})
- set(HTML_FILE_DEPENDS "${NB_FILE_BASE}.run${NB_FILE_EXT}")
+ set(NB_FILE_RUN "${NB_FILE_BASE}.run${NB_FILE_EXT}")
add_custom_command(
- OUTPUT "${HTML_FILE_DEPENDS}"
+ OUTPUT ${NB_FILE_RUN}
DEPENDS
"${NB_FILE};${NB_EXPORT_ADD_SCRIPTS};${CMAKE_BINARY_DIR}/doc/tutorials/html_runner.py;${CMAKE_BINARY_DIR}/testsuite/scripts/importlib_wrapper.py"
COMMAND
- "${CMAKE_BINARY_DIR}/pypresso"
- "${CMAKE_BINARY_DIR}/doc/tutorials/html_runner.py" "--input"
- "${NB_FILE}" "--output" "${HTML_FILE_DEPENDS}" "--substitutions"
- ${NB_EXPORT_VAR_SUBST} "--scripts" ${NB_EXPORT_ADD_SCRIPTS})
+ ${CMAKE_BINARY_DIR}/pypresso
+ ${CMAKE_BINARY_DIR}/doc/tutorials/html_runner.py --execute --exercise2
+ --input ${NB_FILE} --output ${NB_FILE_RUN} --substitutions
+ ${NB_EXPORT_VAR_SUBST} --scripts ${NB_EXPORT_ADD_SCRIPTS})
else()
- set(HTML_FILE_DEPENDS "${NB_FILE}")
+ set(NB_FILE_RUN ${NB_FILE})
endif()
add_custom_command(
- OUTPUT ${HTML_FILE} DEPENDS ${HTML_FILE_DEPENDS};${NB_EXPORT_ADD_SCRIPTS}
+ OUTPUT ${HTML_FILE}
+ DEPENDS ${NB_FILE_RUN};${NB_EXPORT_ADD_SCRIPTS}
+ COMMAND
+ ${CMAKE_BINARY_DIR}/pypresso
+ ${CMAKE_BINARY_DIR}/doc/tutorials/html_runner.py --exercise2 --input
+ ${NB_FILE_RUN} --output ${NB_FILE_RUN}~
COMMAND ${IPYTHON_EXECUTABLE} nbconvert --to "html" --output ${HTML_FILE}
- ${HTML_FILE_DEPENDS})
+ ${NB_FILE_RUN}~)
add_custom_command(
- OUTPUT ${PY_FILE} DEPENDS ${NB_FILE}
+ OUTPUT ${PY_FILE}
+ DEPENDS ${NB_FILE}
+ COMMAND
+ ${CMAKE_BINARY_DIR}/pypresso
+ ${CMAKE_BINARY_DIR}/doc/tutorials/html_runner.py --exercise2 --input
+ ${NB_FILE} --output ${NB_FILE}~
COMMAND ${IPYTHON_EXECUTABLE} nbconvert --to "python" --output ${PY_FILE}
- ${NB_FILE})
+ ${NB_FILE}~)
add_custom_target("${NB_EXPORT_TARGET}_html" DEPENDS ${HTML_FILE}
${DEPENDENCY_OF_TARGET})
diff --git a/doc/tutorials/html_runner.py b/doc/tutorials/html_runner.py
index 188b9e874de..f067d35e0ca 100644
--- a/doc/tutorials/html_runner.py
+++ b/doc/tutorials/html_runner.py
@@ -17,28 +17,33 @@
# along with this program. If not, see .
#
"""
-This script runs Jupyter notebooks and writes the output to new notebooks.
-Global variables can be edited to reduce runtime. External Python scripts
+This script processes Jupyter notebooks. External Python scripts
can be inserted as new code cells (e.g. solutions to exercises).
-The output notebooks can then be converted to HTML externally.
+Hidden solutions from the ``exercise2`` plugin can be converted
+to code cells. The notebook may also be executed, if necessary
+with modified global variables to reduce runtime. The processed
+notebook can then be converted to HTML externally.
"""
import argparse
parser = argparse.ArgumentParser(description='Process Jupyter notebooks.',
epilog=__doc__)
-parser.add_argument('--input', type=str,
+parser.add_argument('--input', type=str, nargs=1, required=True,
help='Path to the original Jupyter notebook')
-parser.add_argument('--output', type=str, nargs='?',
+parser.add_argument('--output', type=str, nargs=1,
help='Path to the processed Jupyter notebook')
parser.add_argument('--substitutions', nargs='*',
help='Variables to substitute')
parser.add_argument('--scripts', nargs='*',
help='Scripts to insert in new cells')
+parser.add_argument('--exercise2', action='store_true',
+ help='Convert exercise2 solutions into code cells')
+parser.add_argument('--execute', action='store_true',
+ help='Run the script')
args = parser.parse_args()
import nbformat
-from nbconvert.preprocessors import ExecutePreprocessor
import re
import os
import ast
@@ -60,22 +65,15 @@ def set_code_cells(nb, new_cells):
i += 1
-notebook_filepath = args.input
-notebook_filepath_edited = args.output or args.input + '~'
-notebook_dirname = os.path.dirname(notebook_filepath)
-new_values = args.substitutions or []
-new_cells = args.scripts or []
-
-# parse original notebook
-with open(notebook_filepath, encoding='utf-8') as f:
- nb = nbformat.read(f, as_version=4)
-
-# add new cells containing the solutions
-for filepath in new_cells:
+def add_cell_from_script(nb, filepath):
+ """
+ Create new code cell at the end of a notebook and populate it with
+ the content of a script.
+ """
with open(filepath, encoding='utf-8') as f:
code = f.read()
# remove ESPResSo copyright header
- m = re.search('# Copyright \(C\) \d+(?:-\d+)? The ESPResSo project\n.+?'
+ m = re.search('# Copyright \(C\) [\d\-,]+ The ESPResSo project\n.+?'
'If not, see \.\n', code, re.DOTALL)
if m and all(x.startswith('#') for x in m.group(0).strip().split('\n')):
code = re.sub('^(#\n)+', '', code.replace(m.group(0), ''), re.M)
@@ -91,57 +89,133 @@ def set_code_cells(nb, new_cells):
nb['cells'].append(cell_code)
+def disable_plot_interactivity(nb):
+ """
+ Replace all occurrences of the magic command ``%matplotlib notebook``
+ by ``%matplotlib inline``.
+ """
+ for cell in nb['cells']:
+ if cell['cell_type'] == 'code' and 'matplotlib' in cell['source']:
+ cell['source'] = re.sub('^%matplotlib +notebook',
+ '%matplotlib inline',
+ cell['source'], flags=re.M)
+
+
+def split_matplotlib_cells(nb):
+ """
+ If a cell imports matplotlib, split the cell to keep the
+ import statement separate from the code that uses matplotlib.
+ This prevents a known bug in the Jupyter backend which causes
+ the plot object to be represented as a string instead of a canvas
+ when created in the cell where matplotlib is imported for the
+ first time (https://github.com/jupyter/notebook/issues/3523).
+ """
+ for i in range(len(nb['cells']) - 1, -1, -1):
+ cell = nb['cells'][i]
+ if cell['cell_type'] == 'code' and 'matplotlib' in cell['source']:
+ code = iw.protect_ipython_magics(cell['source'])
+ # split cells after matplotlib imports
+ mapping = iw.delimit_statements(code)
+ tree = ast.parse(code)
+ visitor = iw.GetMatplotlibPyplot()
+ visitor.visit(tree)
+ if visitor.matplotlib_first:
+ code = iw.deprotect_ipython_magics(code)
+ lines = code.split('\n')
+ lineno_end = mapping[visitor.matplotlib_first]
+ split_code = '\n'.join(lines[lineno_end:]).lstrip('\n')
+ if split_code:
+ new_cell = nbformat.v4.new_code_cell(source=split_code)
+ nb['cells'].insert(i + 1, new_cell)
+ lines = lines[:lineno_end]
+ nb['cells'][i]['source'] = '\n'.join(lines).rstrip('\n')
+
+
+def convert_exercise2_to_code(nb):
+ """
+ Walk through the notebook cells and remove metadata associated with
+ the ``exercise2`` plugin from the contributed nbextensions. Solution
+ Markdown cells containing python code are converted to code cells.
+ """
+ for i in range(len(nb['cells']) - 1, 0, -1):
+ cell = nb['cells'][i]
+ cell_above = nb['cells'][i - 1]
+ # remove empty code cells after a solution cell
+ if cell['cell_type'] == 'code' and cell['source'].strip() == '' \
+ and 'solution2' in cell_above['metadata'] \
+ and 'solution2_first' not in cell_above['metadata'] \
+ and 'solution2' not in cell['metadata']:
+ nb['cells'].pop(i)
+ continue
+ # convert solution markdown cells into code cells
+ if cell['cell_type'] == 'markdown' and 'solution2' in cell['metadata'] \
+ and 'solution2_first' not in cell['metadata']:
+ lines = cell['source'].strip().split('\n')
+ if lines[0].startswith(
+ '```python') and lines[-1].startswith('```'):
+ source = '\n'.join(lines[1:-1]).strip()
+ nb['cells'][i] = nbformat.v4.new_code_cell(source=source)
+ # remove exercise2 metadata
+ for key in ('solution2', 'solution2_first'):
+ if key in cell['metadata']:
+ del cell['metadata'][key]
+
+
+def execute_notebook(nb, src, cell_separator):
+ """
+ Run the notebook in a python3 kernel. The ESPResSo visualizers are
+ disabled to prevent the kernel from crashing and to allow running
+ the notebook in a CI environment.
+ """
+ import nbconvert.preprocessors
+ notebook_dirname = os.path.dirname(notebook_filepath)
+ # disable OpenGL/Mayavi GUI
+ src_no_gui = iw.mock_es_visualization(src)
+ # update notebook with new code
+ set_code_cells(nb, src_no_gui.split(cell_separator))
+ # execute notebook
+ ep = nbconvert.preprocessors.ExecutePreprocessor(
+ timeout=20 * 60, kernel_name='python3')
+ ep.preprocess(nb, {'metadata': {'path': notebook_dirname}})
+ # restore notebook with code before the GUI removal step
+ set_code_cells(nb, src.split(cell_separator))
+
+
+notebook_filepath = args.input[0]
+if args.output:
+ notebook_filepath_edited = args.output[0]
+else:
+ notebook_filepath_edited = notebook_filepath + '~'
+
+# parse original notebook
+with open(notebook_filepath, encoding='utf-8') as f:
+ nb = nbformat.read(f, as_version=4)
+
+# add new cells containing the solutions
+if args.scripts:
+ for filepath in args.scripts:
+ add_cell_from_script(nb, filepath)
+
# disable plot interactivity
-for i in range(len(nb['cells'])):
- cell = nb['cells'][i]
- if cell['cell_type'] == 'code' and 'matplotlib' in cell['source']:
- cell['source'] = re.sub('^%matplotlib +notebook', '%matplotlib inline',
- cell['source'], flags=re.M)
-
-
-# if matplotlib is used in this script, split cell to keep the import
-# statement separate and avoid a known bug in the Jupyter backend which
-# causes the plot object to be represented as a string instead of a
-# canvas when created in the cell where matplotlib is imported for the
-# first time (https://github.com/jupyter/notebook/issues/3523)
-for i in range(len(nb['cells'])):
- cell = nb['cells'][i]
- if cell['cell_type'] == 'code' and 'matplotlib' in cell['source']:
- code = iw.protect_ipython_magics(cell['source'])
- # split cells after matplotlib imports
- mapping = iw.delimit_statements(code)
- tree = ast.parse(code)
- visitor = iw.GetMatplotlibPyplot()
- visitor.visit(tree)
- if visitor.matplotlib_first:
- code = iw.deprotect_ipython_magics(code)
- lines = code.split('\n')
- lineno_end = mapping[visitor.matplotlib_first]
- split_code = '\n'.join(lines[lineno_end:]).lstrip('\n')
- if split_code:
- new_cell = nbformat.v4.new_code_cell(source=split_code)
- nb['cells'].insert(i + 1, new_cell)
- lines = lines[:lineno_end]
- nb['cells'][i]['source'] = '\n'.join(lines).rstrip('\n')
- break
-
-# substitute global variables and disable OpenGL/Mayavi GUI
-cell_separator = '\n##{}\n'.format(uuid.uuid4().hex)
-src = cell_separator.join(get_code_cells(nb))
-parameters = dict(x.split('=', 1) for x in new_values)
-src = iw.substitute_variable_values(src, strings_as_is=True,
- keep_original=False, **parameters)
-src_no_gui = iw.mock_es_visualization(src)
-
-# update notebook with new code
-set_code_cells(nb, src_no_gui.split(cell_separator))
-
-# execute notebook
-ep = ExecutePreprocessor(timeout=20 * 60, kernel_name='python3')
-ep.preprocess(nb, {'metadata': {'path': notebook_dirname}})
-
-# restore notebook with code before the GUI removal step
-set_code_cells(nb, src.split(cell_separator))
+disable_plot_interactivity(nb)
+
+# guard against a jupyter bug involving matplotlib
+split_matplotlib_cells(nb)
+
+if args.exercise2:
+ convert_exercise2_to_code(nb)
+
+if args.substitutions or args.execute:
+ # substitute global variables
+ cell_separator = '\n##{}\n'.format(uuid.uuid4().hex)
+ src = cell_separator.join(get_code_cells(nb))
+ new_values = args.substitutions or []
+ parameters = dict(x.split('=', 1) for x in new_values)
+ src = iw.substitute_variable_values(src, strings_as_is=True,
+ keep_original=False, **parameters)
+ set_code_cells(nb, src.split(cell_separator))
+ if args.execute:
+ execute_notebook(nb, src, cell_separator)
# write edited notebook
with open(notebook_filepath_edited, 'w', encoding='utf-8') as f:
diff --git a/maintainer/CI/build_docker.sh b/maintainer/CI/build_docker.sh
index 30fc1bec729..75a24ebc898 100755
--- a/maintainer/CI/build_docker.sh
+++ b/maintainer/CI/build_docker.sh
@@ -31,6 +31,6 @@ make_check=${make_check}
build_type=RelWithAssert
myconfig=maxset
EOF
-image="espressomd/docker-ubuntu-20.04:06b6216c7aa3555bcf28c90734dbb84e7285c96f"
+image="espressomd/docker-ubuntu-20.04:e583d4b2eb8eedd10068957f952bd67008475ee5"
docker run -u espresso --env-file "${ENV_FILE}" -v "${PWD}:/travis" -it "${image}" /bin/bash -c "cp -r /travis .; cd travis && maintainer/CI/build_cmake.sh" || exit 1
diff --git a/maintainer/benchmarks/CMakeLists.txt b/maintainer/benchmarks/CMakeLists.txt
index 7c897cb44ce..eb15f156e0e 100644
--- a/maintainer/benchmarks/CMakeLists.txt
+++ b/maintainer/benchmarks/CMakeLists.txt
@@ -114,6 +114,6 @@ python_benchmark(
add_custom_target(
benchmark_python COMMAND ${CMAKE_CTEST_COMMAND} --timeout ${TEST_TIMEOUT}
- $(ARGS) --output-on-failure)
+ ${CTEST_ARGS} --output-on-failure)
add_dependencies(benchmark benchmark_python)
diff --git a/requirements.txt b/requirements.txt
index 5fd6d5c3a11..86d78b471d9 100644
--- a/requirements.txt
+++ b/requirements.txt
@@ -20,6 +20,8 @@ lxml>=4.2.1 # to deploy tutorials
# sphinx and its dependencies
sphinx>=1.6.7,!=2.1.0,!=3.0.0
sphinxcontrib-bibtex>=0.3.5
+# jupyter dependencies
+jupyter_contrib_nbextensions==0.5.1
# pep8 and its dependencies
autopep8==1.3.4
pycodestyle==2.3.1
diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt
index 6e89d17a144..f9a119aa68a 100644
--- a/src/CMakeLists.txt
+++ b/src/CMakeLists.txt
@@ -21,7 +21,7 @@
# Target for the unit tests
add_custom_target(
check_unit_tests COMMAND ${CMAKE_CTEST_COMMAND} --timeout ${TEST_TIMEOUT}
- $(ARGS) --output-on-failure)
+ ${CTEST_ARGS} --output-on-failure)
if(WITH_TESTS)
# Run unit tests on check
diff --git a/src/config/config.hpp b/src/config/config.hpp
index 51ae73cdefd..f21a4d7c0f6 100644
--- a/src/config/config.hpp
+++ b/src/config/config.hpp
@@ -39,23 +39,6 @@
#include "config-features.hpp"
-/** P3M: Default for number of interpolation points of the charge
- * assignment function.
- */
-#ifndef P3M_N_INTERPOL
-#define P3M_N_INTERPOL 32768
-#endif
-
-/** P3M: Default for boundary condition: Epsilon of the surrounding medium. */
-#ifndef P3M_EPSILON
-#define P3M_EPSILON 0.0
-#endif
-
-/** P3M: Default for boundary condition in magnetic calculations. */
-#ifndef P3M_EPSILON_MAGNETIC
-#define P3M_EPSILON_MAGNETIC 0.0
-#endif
-
/** P3M: Default for offset of first mesh point from the origin (left
* down corner of the simulation box).
*/
@@ -117,23 +100,4 @@
#define MAX_OBJECTS_IN_FLUID 10000
#endif
-/** @name Mathematical constants, from gcc's math.hpp */
-/*@{*/
-#ifndef M_PI
-#define M_E 2.7182818284590452353602874713526625L /* e */
-#define M_LOG2E 1.4426950408889634073599246810018921L /* log_2 e */
-#define M_LOG10E 0.4342944819032518276511289189166051L /* log_10 e */
-#define M_LN2 0.6931471805599453094172321214581766L /* log_e 2 */
-#define M_LN10 2.3025850929940456840179914546843642L /* log_e 10 */
-#define M_PI 3.1415926535897932384626433832795029L /* pi */
-#define M_PI_2 1.5707963267948966192313216916397514L /* pi/2 */
-#define M_PI_4 0.7853981633974483096156608458198757L /* pi/4 */
-#define M_1_PI 0.3183098861837906715377675267450287L /* 1/pi */
-#define M_2_PI 0.6366197723675813430755350534900574L /* 2/pi */
-#define M_2_SQRTPI 1.1283791670955125738961589031215452L /* 2/sqrt(pi) */
-#define M_SQRT2 1.4142135623730950488016887242096981L /* sqrt(2) */
-#define M_SQRT1_2 0.7071067811865475244008443621048490L /* 1/sqrt(2) */
-#endif
-/*@}*/
-
#endif
diff --git a/src/core/CMakeLists.txt b/src/core/CMakeLists.txt
index c3b992ccf0f..a2a807c89e9 100644
--- a/src/core/CMakeLists.txt
+++ b/src/core/CMakeLists.txt
@@ -73,20 +73,18 @@ endif(CUDA)
install(TARGETS EspressoCore LIBRARY DESTINATION ${PYTHON_INSTDIR}/espressomd)
target_link_libraries(
- EspressoCore
- PRIVATE EspressoConfig EspressoShapes Profiler
- "$<$:${FFTW3_LIBRARIES}>"
- $<$:Scafacos> cxx_interface
+ EspressoCore PRIVATE EspressoConfig EspressoShapes Profiler
+ $<$:Scafacos> cxx_interface
PUBLIC EspressoUtils MPI::MPI_CXX Random123 EspressoParticleObservables
Boost::serialization Boost::mpi "$<$:${HDF5_LIBRARIES}>"
- $<$:Boost::filesystem> $<$:h5xx>)
+ $<$:Boost::filesystem> $<$:h5xx>
+ "$<$:FFTW3::FFTW3>")
target_include_directories(
EspressoCore
PUBLIC ${CMAKE_CURRENT_SOURCE_DIR}
"$<$:${CMAKE_CURRRENT_SOURCE_DIR}/io/writer>"
- PUBLIC "$<$:${HDF5_INCLUDE_DIRS}>"
- PRIVATE "$<$:${FFTW3_INCLUDE_DIR}>")
+ "$<$:${HDF5_INCLUDE_DIRS}>")
target_compile_definitions(EspressoCore PUBLIC $<$:H5XX_USE_MPI>)
diff --git a/src/core/MpiCallbacks.hpp b/src/core/MpiCallbacks.hpp
index 31f991d6725..a609878ecc4 100644
--- a/src/core/MpiCallbacks.hpp
+++ b/src/core/MpiCallbacks.hpp
@@ -38,6 +38,9 @@
#include
namespace Communication {
+
+class MpiCallbacks;
+
/**
* This namespace contains tag types
* to indicate what to do with the return
@@ -365,7 +368,8 @@ class MpiCallbacks {
-> std::enable_if_t<
std::is_void()(
std::forward(args)...))>::value> {
- assert(m_cb), m_cb->call(m_id, std::forward(args)...);
+ if (m_cb)
+ m_cb->call(m_id, std::forward(args)...);
}
~CallbackHandle() {
@@ -374,6 +378,7 @@ class MpiCallbacks {
}
MpiCallbacks *cb() const { return m_cb; }
+ int id() const { return m_id; }
};
/* Avoid accidental copy, leads to mpi deadlock
@@ -391,9 +396,9 @@ class MpiCallbacks {
}
public:
- explicit MpiCallbacks(boost::mpi::communicator &comm,
+ explicit MpiCallbacks(boost::mpi::communicator comm,
bool abort_on_exit = true)
- : m_abort_on_exit(abort_on_exit), m_comm(comm) {
+ : m_abort_on_exit(abort_on_exit), m_comm(std::move(comm)) {
/** Add a dummy at id 0 for loop abort. */
m_callback_map.add(nullptr);
@@ -719,7 +724,7 @@ class MpiCallbacks {
/**
* The MPI communicator used for the callbacks.
*/
- boost::mpi::communicator &m_comm;
+ boost::mpi::communicator m_comm;
/**
* Internal storage for the callback functions.
diff --git a/src/core/actor/Mmm1dgpuForce.hpp b/src/core/actor/Mmm1dgpuForce.hpp
index 31233a5a698..2531906bbe4 100644
--- a/src/core/actor/Mmm1dgpuForce.hpp
+++ b/src/core/actor/Mmm1dgpuForce.hpp
@@ -16,6 +16,9 @@
* You should have received a copy of the GNU General Public License
* along with this program. If not, see .
*/
+#ifndef ESPRESSO_CORE_ACTOR_MMM1DGPUFORCE_HPP
+#define ESPRESSO_CORE_ACTOR_MMM1DGPUFORCE_HPP
+
#include "config.hpp"
#ifdef MMM1D_GPU
@@ -78,3 +81,4 @@ class Mmm1dgpuForce : public Actor {
};
#endif
+#endif
diff --git a/src/core/actor/specfunc_cuda.hpp b/src/core/actor/specfunc_cuda.hpp
index 76e01589ae2..b376825e49e 100644
--- a/src/core/actor/specfunc_cuda.hpp
+++ b/src/core/actor/specfunc_cuda.hpp
@@ -44,9 +44,12 @@
* first and second kind. The implementations are based on the GSL code (see
* the original GSL header above) and are duplicated from \ref specfunc.cpp.
*/
+#ifndef ESPRESSO_CORE_ACTOR_SPECFUNC_CUDA_HPP
+#define ESPRESSO_CORE_ACTOR_SPECFUNC_CUDA_HPP
+
#include "config.hpp"
-const mmm1dgpu_real M_LN2f = M_LN2;
+#include
/** @name Chebyshev expansions based on SLATEC bk0(), bk0e() */
/*@{*/
@@ -155,7 +158,7 @@ __device__ mmm1dgpu_real dev_K0(mmm1dgpu_real x) {
if (x <= 2) {
mmm1dgpu_real I0 =
evaluateAsChebychevSeriesAt(bi0_data, bi0_size, x * x / 4.5f - 1.0f);
- return (-log(x) + M_LN2f) * I0 + c;
+ return (-log(x) + Utils::ln_2()) * I0 + c;
}
return exp(-x) * c * rsqrt(x);
}
@@ -169,7 +172,8 @@ __device__ mmm1dgpu_real dev_K1(mmm1dgpu_real x) {
if (x <= 2) {
mmm1dgpu_real I1 = x * evaluateAsChebychevSeriesAt(bi1_data, bi1_size,
x * x / 4.5f - 1.0f);
- return (log(x) - M_LN2f) * I1 + c / x;
+ return (log(x) - Utils::ln_2()) * I1 + c / x;
}
return exp(-x) * c * rsqrt(x);
}
+#endif
diff --git a/src/core/collision.cpp b/src/core/collision.cpp
index 43625629abc..7e40e947a04 100644
--- a/src/core/collision.cpp
+++ b/src/core/collision.cpp
@@ -30,6 +30,7 @@
#include "nonbonded_interactions/nonbonded_interaction_data.hpp"
#include "virtual_sites.hpp"
+#include
#include
#include
@@ -77,8 +78,9 @@ Particle &get_part(int id) {
}
} // namespace
-/** @brief Return true if a bond between the centers of the colliding particles
- * needs to be placed. At this point, all modes need this */
+/** @brief Return true if a bond between the centers of the colliding
+ * particles needs to be placed. At this point, all modes need this.
+ */
inline bool bind_centers() {
// Note that the glue to surface mode adds bonds between the centers
// but does so later in the process. This is needed to guarantee that
@@ -246,8 +248,9 @@ void queue_collision(const int part1, const int part2) {
local_collision_queue.push_back({part1, part2});
}
-/** @brief Calculate position of vs for GLUE_TO_SURFACE mode
- * Returns id of particle to bind vs to */
+/** @brief Calculate position of vs for GLUE_TO_SURFACE mode.
+ * Returns id of particle to bind vs to.
+ */
const Particle &glue_to_surface_calc_vs_pos(const Particle &p1,
const Particle &p2,
Utils::Vector3d &pos) {
@@ -339,16 +342,12 @@ void coldet_do_three_particle_bond(Particle &p, Particle const &p1,
// three_particle_angle_resolution steps and by adding the id
// of the bond for zero degrees.
auto const bond_id = static_cast(
- floor(phi / M_PI *
- (collision_params.three_particle_angle_resolution - 1) +
- 0.5) +
+ floor(0.5 + phi / Utils::pi() *
+ (collision_params.three_particle_angle_resolution - 1)) +
collision_params.bond_three_particles);
// Create the bond
-
- // First, fill bond data structure
const std::array bondT = {p1.p.identity, p2.p.identity};
-
p.bonds().insert({bond_id, bondT});
}
diff --git a/src/core/electrostatics_magnetostatics/dp3m_influence_function.hpp b/src/core/electrostatics_magnetostatics/dp3m_influence_function.hpp
new file mode 100644
index 00000000000..af999cff4ba
--- /dev/null
+++ b/src/core/electrostatics_magnetostatics/dp3m_influence_function.hpp
@@ -0,0 +1,223 @@
+/*
+ * Copyright (C) 2010-2020 The ESPResSo project
+ * Copyright (C) 2002-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 .
+ */
+#ifndef ESPRESSO_DP3M_INFLUENCE_FUNCTION_HPP
+#define ESPRESSO_DP3M_INFLUENCE_FUNCTION_HPP
+
+#include "electrostatics_magnetostatics/p3m-common.hpp"
+
+#include
+#include
+#include
+#include
+#include
+#include
+
+#include
+
+#include
+
+#if defined(DP3M)
+
+/** Calculate the aliasing sums for the optimal influence function.
+ *
+ * Calculates the aliasing sums in the numerator and denominator of
+ * the expression for the optimal influence function (see
+ * @cite hockney88a : 8-22, p. 275).
+ *
+ * \tparam S order (2 for energy, 3 for forces)
+ * \param params DP3M parameters
+ * \param shift shift for a given n-vector
+ * \param d_op differential operator for a given n-vector
+ * \return The result of the fraction.
+ */
+template
+double G_opt_dipolar(P3MParameters const ¶ms, Utils::Vector3i const &shift,
+ Utils::Vector3i const &d_op) {
+ using Utils::int_pow;
+ using Utils::sinc;
+ constexpr double limit = 30;
+
+ double numerator = 0.0;
+ double denominator = 0.0;
+
+ auto const f1 = 1.0 / static_cast(params.mesh[0]);
+ auto const f2 = Utils::sqr(Utils::pi() / params.alpha_L);
+
+ for (double mx = -P3M_BRILLOUIN; mx <= P3M_BRILLOUIN; mx++) {
+ auto const nmx = shift[0] + params.mesh[0] * mx;
+ auto const sx = std::pow(sinc(f1 * nmx), 2.0 * params.cao);
+ for (double my = -P3M_BRILLOUIN; my <= P3M_BRILLOUIN; my++) {
+ auto const nmy = shift[1] + params.mesh[0] * my;
+ auto const sy = sx * std::pow(sinc(f1 * nmy), 2.0 * params.cao);
+ for (double mz = -P3M_BRILLOUIN; mz <= P3M_BRILLOUIN; mz++) {
+ auto const nmz = shift[2] + params.mesh[0] * mz;
+ auto const sz = sy * std::pow(sinc(f1 * nmz), 2.0 * params.cao);
+ auto const nm2 = Utils::sqr(nmx) + Utils::sqr(nmy) + Utils::sqr(nmz);
+ auto const exponent = f2 * nm2;
+ if (exponent < limit) {
+ auto const f3 = sz * std::exp(-exponent) / nm2;
+ auto const n_nm = d_op[0] * nmx + d_op[1] * nmy + d_op[2] * nmz;
+ numerator += f3 * int_pow(n_nm);
+ }
+ denominator += sz;
+ }
+ }
+ }
+ return numerator / (int_pow(static_cast(d_op.norm2())) *
+ Utils::sqr(denominator));
+}
+
+/**
+ * @brief Map influence function over a grid.
+ *
+ * This evaluates the optimal influence function @ref G_opt_dipolar
+ * over a regular grid of k vectors, and returns the values as a vector.
+ *
+ * @tparam S Order of the differential operator, e.g. 2 for energy, 3 for force
+ *
+ * @param params DP3M parameters
+ * @param n_start Lower left corner of the grid
+ * @param n_end Upper right corner of the grid.
+ * @param box_l Box size
+ * @return Values of the influence function at regular grid points.
+ */
+template
+std::vector grid_influence_function(P3MParameters const ¶ms,
+ Utils::Vector3i const &n_start,
+ Utils::Vector3i const &n_end,
+ Utils::Vector3d const &box_l) {
+
+ auto const size = n_end - n_start;
+
+ /* The influence function grid */
+ auto g =
+ std::vector(boost::accumulate(size, 1, std::multiplies<>()), 0.);
+
+ /* Skip influence function calculation in tuning mode,
+ the results need not be correct for timing. */
+ if (params.tuning) {
+ return g;
+ }
+
+ double fak1 = Utils::int_pow<3>(static_cast(params.mesh[0])) * 2.0 /
+ Utils::sqr(box_l[0]);
+
+ auto const shifts =
+ detail::calc_meshift({params.mesh[0], params.mesh[1], params.mesh[2]});
+ auto const d_ops = detail::calc_meshift(
+ {params.mesh[0], params.mesh[1], params.mesh[2]}, true);
+
+ Utils::Vector3i n{};
+ for (n[0] = n_start[0]; n[0] < n_end[0]; n[0]++) {
+ for (n[1] = n_start[1]; n[1] < n_end[1]; n[1]++) {
+ for (n[2] = n_start[2]; n[2] < n_end[2]; n[2]++) {
+ auto const ind = Utils::get_linear_index(n - n_start, size,
+ Utils::MemoryOrder::ROW_MAJOR);
+
+ if (((n[0] % (params.mesh[0] / 2) == 0) &&
+ (n[1] % (params.mesh[0] / 2) == 0) &&
+ (n[2] % (params.mesh[0] / 2) == 0))) {
+ g[ind] = 0.0;
+ } else {
+ auto const shift = Utils::Vector3i{shifts[0][n[0]], shifts[0][n[1]],
+ shifts[0][n[2]]};
+ auto const d_op =
+ Utils::Vector3i{d_ops[0][n[0]], d_ops[0][n[1]], d_ops[0][n[2]]};
+ auto const fak2 = G_opt_dipolar(params, shift, d_op);
+ g[ind] = fak1 * fak2;
+ }
+ }
+ }
+ }
+ return g;
+}
+
+double G_opt_dipolar_self_energy(P3MParameters const ¶ms,
+ Utils::Vector3i const &shift) {
+ using Utils::sinc;
+ double u_sum = 0.0;
+ constexpr int limit = P3M_BRILLOUIN + 5;
+
+ auto const f1 = 1.0 / static_cast(params.mesh[0]);
+
+ for (double mx = -limit; mx <= limit; mx++) {
+ auto const nmx = shift[0] + params.mesh[0] * mx;
+ auto const sx = std::pow(sinc(f1 * nmx), 2.0 * params.cao);
+ for (double my = -limit; my <= limit; my++) {
+ auto const nmy = shift[1] + params.mesh[0] * my;
+ auto const sy = sx * std::pow(sinc(f1 * nmy), 2.0 * params.cao);
+ for (double mz = -limit; mz <= limit; mz++) {
+ auto const nmz = shift[2] + params.mesh[0] * mz;
+ auto const sz = sy * std::pow(sinc(f1 * nmz), 2.0 * params.cao);
+ u_sum += sz;
+ }
+ }
+ }
+ return u_sum;
+}
+
+/**
+ * @brief Calculate self-energy of the influence function.
+ *
+ * @param params DP3M parameters
+ * @param n_start Lower left corner of the grid
+ * @param n_end Upper right corner of the grid.
+ * @param g Energies on the grid.
+ * @return Total self-energy.
+ */
+double grid_influence_function_self_energy(P3MParameters const ¶ms,
+ Utils::Vector3i const &n_start,
+ Utils::Vector3i const &n_end,
+ std::vector const &g) {
+ auto const size = n_end - n_start;
+
+ auto const shifts =
+ detail::calc_meshift({params.mesh[0], params.mesh[1], params.mesh[2]});
+ auto const d_ops = detail::calc_meshift(
+ {params.mesh[0], params.mesh[1], params.mesh[2]}, true);
+
+ double energy = 0.0;
+ Utils::Vector3i n{};
+ for (n[0] = n_start[0]; n[0] < n_end[0]; n[0]++) {
+ for (n[1] = n_start[1]; n[1] < n_end[1]; n[1]++) {
+ for (n[2] = n_start[2]; n[2] < n_end[2]; n[2]++) {
+ if (((n[0] % (params.mesh[0] / 2) == 0) &&
+ (n[1] % (params.mesh[0] / 2) == 0) &&
+ (n[2] % (params.mesh[0] / 2) == 0))) {
+ energy += 0.0;
+ } else {
+ auto const ind = Utils::get_linear_index(
+ n - n_start, size, Utils::MemoryOrder::ROW_MAJOR);
+ auto const shift = Utils::Vector3i{shifts[0][n[0]], shifts[0][n[1]],
+ shifts[0][n[2]]};
+ auto const d_op =
+ Utils::Vector3i{d_ops[0][n[0]], d_ops[0][n[1]], d_ops[0][n[2]]};
+ auto const U2 = G_opt_dipolar_self_energy(params, shift);
+ energy += g[ind] * U2 * d_op.norm2();
+ }
+ }
+ }
+ }
+ return energy;
+}
+
+#endif
+#endif
diff --git a/src/core/electrostatics_magnetostatics/elc.cpp b/src/core/electrostatics_magnetostatics/elc.cpp
index 475fcbe0340..a938f692349 100644
--- a/src/core/electrostatics_magnetostatics/elc.cpp
+++ b/src/core/electrostatics_magnetostatics/elc.cpp
@@ -220,7 +220,7 @@ void distribute(int size) {
* See @cite yeh99a.
*/
static void add_dipole_force(const ParticleRange &particles) {
- double const pref = coulomb.prefactor * 4 * M_PI * ux * uy * uz;
+ double const pref = coulomb.prefactor * 4 * Utils::pi() * ux * uy * uz;
int const size = 3;
auto local_particles = particles;
@@ -283,7 +283,7 @@ static void add_dipole_force(const ParticleRange &particles) {
* See @cite yeh99a.
*/
static double dipole_energy(const ParticleRange &particles) {
- double const pref = coulomb.prefactor * 2 * M_PI * ux * uy * uz;
+ double const pref = coulomb.prefactor * 2 * Utils::pi() * ux * uy * uz;
int const size = 7;
/* for nonneutral systems, this shift gives the background contribution
(rsp. for this shift, the DM of the background is zero) */
@@ -375,7 +375,7 @@ inline double image_sum_t(double q, double z) {
/*****************************************************************/
static double z_energy(const ParticleRange &particles) {
- double const pref = coulomb.prefactor * 2 * M_PI * ux * uy;
+ double const pref = coulomb.prefactor * 2 * Utils::pi() * ux * uy;
int const size = 4;
/* for nonneutral systems, this shift gives the background contribution
@@ -453,7 +453,7 @@ static double z_energy(const ParticleRange &particles) {
/*****************************************************************/
static void add_z_force(const ParticleRange &particles) {
- double const pref = coulomb.prefactor * 2 * M_PI * ux * uy;
+ double const pref = coulomb.prefactor * 2 * Utils::pi() * ux * uy;
if (elc_params.dielectric_contrast_on) {
auto local_particles = particles;
@@ -508,9 +508,8 @@ static void add_z_force(const ParticleRange &particles) {
/*****************************************************************/
static void setup_P(int p, double omega, const ParticleRange &particles) {
- double const pref = -coulomb.prefactor * 4 * M_PI * ux * uy /
- (expm1(omega * box_geo.length()[2]));
- double const pref_di = coulomb.prefactor * 4 * M_PI * ux * uy;
+ double const pref_di = coulomb.prefactor * 4 * Utils::pi() * ux * uy;
+ double const pref = -pref_di / expm1(omega * box_geo.length()[2]);
int const size = 4;
double lclimgebot[4], lclimgetop[4], lclimge[4];
double fac_delta_mid_bot = 1, fac_delta_mid_top = 1, fac_delta = 1;
@@ -613,9 +612,8 @@ static void setup_P(int p, double omega, const ParticleRange &particles) {
}
static void setup_Q(int q, double omega, const ParticleRange &particles) {
- double const pref = -coulomb.prefactor * 4 * M_PI * ux * uy /
- (expm1(omega * box_geo.length()[2]));
- double const pref_di = coulomb.prefactor * 4 * M_PI * ux * uy;
+ double const pref_di = coulomb.prefactor * 4 * Utils::pi() * ux * uy;
+ double const pref = -pref_di / expm1(omega * box_geo.length()[2]);
int const size = 4;
double lclimgebot[4], lclimgetop[4], lclimge[4];
double fac_delta_mid_bot = 1, fac_delta_mid_top = 1, fac_delta = 1;
@@ -787,9 +785,8 @@ static double Q_energy(double omega, int n_part) {
static void setup_PQ(int p, int q, double omega,
const ParticleRange &particles) {
- double const pref = -coulomb.prefactor * 8 * M_PI * ux * uy /
- (expm1(omega * box_geo.length()[2]));
- double const pref_di = coulomb.prefactor * 8 * M_PI * ux * uy;
+ double const pref_di = coulomb.prefactor * 8 * Utils::pi() * ux * uy;
+ double const pref = -pref_di / expm1(omega * box_geo.length()[2]);
int const size = 8;
double lclimgebot[8], lclimgetop[8], lclimge[8];
double fac_delta_mid_bot = 1, fac_delta_mid_top = 1, fac_delta = 1;
diff --git a/src/core/electrostatics_magnetostatics/fft.cpp b/src/core/electrostatics_magnetostatics/fft.cpp
index 1ab7e88b839..c850e2acf75 100644
--- a/src/core/electrostatics_magnetostatics/fft.cpp
+++ b/src/core/electrostatics_magnetostatics/fft.cpp
@@ -71,6 +71,7 @@ namespace {
* \param[out] node_list2 Linear node index list for grid2.
* \param[out] pos Positions of the nodes in grid2
* \param[out] my_pos Position of comm.rank() in grid2.
+ * \param[in] comm MPI communicator.
* \return Size of the communication group.
*/
boost::optional>
diff --git a/src/core/electrostatics_magnetostatics/mdlc_correction.cpp b/src/core/electrostatics_magnetostatics/mdlc_correction.cpp
index 2c6d95a5cbd..62b6b3bf69f 100644
--- a/src/core/electrostatics_magnetostatics/mdlc_correction.cpp
+++ b/src/core/electrostatics_magnetostatics/mdlc_correction.cpp
@@ -34,6 +34,8 @@
#include "grid.hpp"
#include "particle_data.hpp"
+#include
+
DLC_struct dlc_params = {1e100, 0, 0, 0, 0};
static double mu_max;
@@ -106,8 +108,8 @@ double get_DLC_dipolar(int kcut, std::vector &fs,
double s1z, s2z, s3z, s4z;
double ss;
- auto const facux = 2.0 * M_PI / box_geo.length()[0];
- auto const facuy = 2.0 * M_PI / box_geo.length()[1];
+ auto const facux = 2.0 * Utils::pi() / box_geo.length()[0];
+ auto const facuy = 2.0 * Utils::pi() / box_geo.length()[1];
double energy = 0.0;
for (int ix = -kcut; ix <= +kcut; ix++) {
@@ -223,7 +225,7 @@ double get_DLC_dipolar(int kcut, std::vector &fs,
// Multiply by the factors we have left during the loops
- auto const piarea = M_PI / (box_geo.length()[0] * box_geo.length()[1]);
+ auto const piarea = Utils::pi() / (box_geo.length()[0] * box_geo.length()[1]);
for (int j = 0; j < n_local_particles; j++) {
fs[j] *= piarea;
@@ -239,8 +241,8 @@ double get_DLC_dipolar(int kcut, std::vector &fs,
* %Algorithm implemented accordingly to @cite brodka04a.
*/
double get_DLC_energy_dipolar(int kcut, const ParticleRange &particles) {
- auto const facux = 2.0 * M_PI / box_geo.length()[0];
- auto const facuy = 2.0 * M_PI / box_geo.length()[1];
+ auto const facux = 2.0 * Utils::pi() / box_geo.length()[0];
+ auto const facuy = 2.0 * Utils::pi() / box_geo.length()[1];
double energy = 0.0;
for (int ix = -kcut; ix <= +kcut; ix++) {
@@ -294,7 +296,7 @@ double get_DLC_energy_dipolar(int kcut, const ParticleRange &particles) {
// Multiply by the factors we have left during the loops
- auto const piarea = M_PI / (box_geo.length()[0] * box_geo.length()[1]);
+ auto const piarea = Utils::pi() / (box_geo.length()[0] * box_geo.length()[1]);
energy *= (-piarea);
return (this_node == 0) ? energy : 0.0;
}
@@ -337,7 +339,7 @@ void add_mdlc_force_corrections(const ParticleRange &particles) {
#if defined(ROTATION) && defined(DP3M)
auto const dip = p.calc_dip();
- auto const correc = 4. * M_PI / volume;
+ auto const correc = 4. * Utils::pi() / volume;
Utils::Vector3d d;
// in the Next lines: the second term (correc*...) is the SDC
// correction for the torques
@@ -385,16 +387,16 @@ double add_mdlc_energy_corrections(const ParticleRange &particles) {
#ifdef DP3M
if (dipole.method == DIPOLAR_MDLC_P3M) {
if (dp3m.params.epsilon == P3M_EPSILON_METALLIC) {
- dip_DLC_energy += dipole.prefactor * 2. * M_PI / volume * mz2;
+ dip_DLC_energy += dipole.prefactor * 2 * Utils::pi() / volume * mz2;
} else {
dip_DLC_energy +=
- dipole.prefactor * 2. * M_PI / volume *
+ dipole.prefactor * 2 * Utils::pi() / volume *
(mz2 - mtot * mtot / (2.0 * dp3m.params.epsilon + 1.0));
}
} else
#endif
{
- dip_DLC_energy += dipole.prefactor * 2. * M_PI / volume * mz2;
+ dip_DLC_energy += dipole.prefactor * 2 * Utils::pi() / volume * mz2;
fprintf(stderr, "You are not using the P3M method, therefore "
"dp3m.params.epsilon unknown, I assume metallic borders "
"\n");
diff --git a/src/core/electrostatics_magnetostatics/mmm1d.cpp b/src/core/electrostatics_magnetostatics/mmm1d.cpp
index bb07549360d..070103da13c 100644
--- a/src/core/electrostatics_magnetostatics/mmm1d.cpp
+++ b/src/core/electrostatics_magnetostatics/mmm1d.cpp
@@ -79,8 +79,8 @@ static std::vector bessel_radii;
static double far_error(int P, double minrad) {
// this uses an upper bound to all force components and the potential
- auto const rhores = 2 * M_PI * uz * minrad;
- auto const pref = 4 * uz * std::max(1.0, 2 * M_PI * uz);
+ auto const rhores = 2 * Utils::pi() * uz * minrad;
+ auto const pref = 4 * uz * std::max(1.0, 2 * Utils::pi() * uz);
return pref * K1(rhores * P) * exp(rhores) / rhores * (P - 1 + 1 / rhores);
}
@@ -318,7 +318,7 @@ double mmm1d_coulomb_pair_energy(double const chpref, Utils::Vector3d const &d,
auto const rxy_d = rxy * uz;
/* The first Bessel term will compensate a little bit the
log term, so add them close together */
- E = -0.25 * log(rxy2_d) + 0.5 * (M_LN2 - Utils::gamma());
+ E = -0.25 * log(rxy2_d) + 0.5 * (Utils::ln_2() - Utils::gamma());
for (int bp = 1; bp < MAXIMAL_B_CUT; bp++) {
if (bessel_radii[bp - 1] < rxy)
break;
diff --git a/src/core/electrostatics_magnetostatics/p3m-common.hpp b/src/core/electrostatics_magnetostatics/p3m-common.hpp
index 87da9c0665f..cc058006b8a 100644
--- a/src/core/electrostatics_magnetostatics/p3m-common.hpp
+++ b/src/core/electrostatics_magnetostatics/p3m-common.hpp
@@ -34,6 +34,9 @@
*/
#include "config.hpp"
+#include
+#include
+
#if defined(P3M) || defined(DP3M)
#include "LocalBox.hpp"
@@ -55,11 +58,20 @@ enum P3M_TUNE_ERROR {
P3M_TUNE_ACCURACY_TOO_LARGE = 32
};
-/** This value for p3m.epsilon indicates metallic boundary conditions. */
+namespace detail {
+/** @brief Index helpers for direct and reciprocal space.
+ * After the FFT the data is in order YZX, which
+ * means that Y is the slowest changing index.
+ */
+namespace FFT_indexing {
+enum FFT_REAL_VECTOR : int { RX = 0, RY = 1, RZ = 2 };
+enum FFT_WAVE_VECTOR : int { KY = 0, KZ = 1, KX = 2 };
+} // namespace FFT_indexing
+} // namespace detail
+
+/** This value indicates metallic boundary conditions. */
#define P3M_EPSILON_METALLIC 0.0
-/** increment size of charge assignment fields. */
-#define CA_INCREMENT 32
/** precision limit for the r_cut zero */
#define P3M_RCUT_PREC 1e-3
/** granularity of the time measurement */
@@ -118,7 +130,7 @@ typedef struct {
double accuracy = 0.0;
/** epsilon of the "surrounding dielectric". */
- double epsilon = P3M_EPSILON;
+ double epsilon = P3M_EPSILON_METALLIC;
/** cutoff for charge assignment. */
double cao_cut[3] = {};
/** mesh constant. */
@@ -187,4 +199,31 @@ void p3m_calc_lm_ld_pos(p3m_local_mesh &local_mesh,
#endif /* P3M || DP3M */
+namespace detail {
+/** Calculate indices that shift @ref P3MParameters::mesh "mesh" by `mesh/2`.
+ * For each mesh size @f$ n @f$ in @c mesh_size, create a sequence of integer
+ * values @f$ \left( 0, \ldots, \lfloor n/2 \rfloor, -\lfloor n/2 \rfloor,
+ * \ldots, -1\right) @f$ if @c zero_out_midpoint is false, otherwise
+ * @f$ \left( 0, \ldots, \lfloor n/2 - 1 \rfloor, 0, -\lfloor n/2 \rfloor,
+ * \ldots, -1\right) @f$.
+ */
+std::array, 3> inline calc_meshift(
+ std::array const &mesh_size, bool zero_out_midpoint = false) {
+ std::array, 3> ret{};
+
+ for (size_t i = 0; i < 3; i++) {
+ ret[i] = std::vector(mesh_size[i]);
+
+ for (int j = 1; j <= mesh_size[i] / 2; j++) {
+ ret[i][j] = j;
+ ret[i][mesh_size[i] - j] = -j;
+ }
+ if (zero_out_midpoint)
+ ret[i][mesh_size[i] / 2] = 0;
+ }
+
+ return ret;
+}
+} // namespace detail
+
#endif /* _P3M_COMMON_H */
diff --git a/src/core/electrostatics_magnetostatics/p3m-data_struct.hpp b/src/core/electrostatics_magnetostatics/p3m-data_struct.hpp
new file mode 100644
index 00000000000..93e0e4149fb
--- /dev/null
+++ b/src/core/electrostatics_magnetostatics/p3m-data_struct.hpp
@@ -0,0 +1,55 @@
+/*
+ * Copyright (C) 2010-2020 The ESPResSo project
+ * Copyright (C) 2002-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 .
+ */
+#ifndef P3M_DATA_STRUCT_HPP
+#define P3M_DATA_STRUCT_HPP
+
+#include "config.hpp"
+
+#ifdef P3M
+
+#include "p3m-common.hpp"
+
+struct p3m_data_struct_base {
+ P3MParameters params;
+
+ /** Spatial differential operator in k-space. We use an i*k differentiation.
+ */
+ std::array, 3> d_op;
+ /** Force optimised influence function (k-space) */
+ std::vector g_force;
+ /** Energy optimised influence function (k-space) */
+ std::vector g_energy;
+
+ /** number of permutations in k_space */
+ int ks_pnum;
+
+ /** Calculate the Fourier transformed differential operator.
+ * Remark: This is done on the level of n-vectors and not k-vectors,
+ * i.e. the prefactor @f$ 2i\pi/L @f$ is missing!
+ */
+ void calc_differential_operator() {
+ d_op = detail::calc_meshift(
+ {params.mesh[0], params.mesh[1], params.mesh[2]}, true);
+ }
+};
+
+#endif
+#endif
diff --git a/src/core/electrostatics_magnetostatics/p3m-dipolar.cpp b/src/core/electrostatics_magnetostatics/p3m-dipolar.cpp
index 12c2ce831e7..fba65eab250 100644
--- a/src/core/electrostatics_magnetostatics/p3m-dipolar.cpp
+++ b/src/core/electrostatics_magnetostatics/p3m-dipolar.cpp
@@ -33,6 +33,7 @@
*/
#include "electrostatics_magnetostatics/p3m-dipolar.hpp"
+#include "electrostatics_magnetostatics/dp3m_influence_function.hpp"
#ifdef DP3M
@@ -46,11 +47,10 @@
#include
using Utils::strcat_alloc;
-#include
-using Utils::sinc;
#include
#include
#include
+#include
#include
#include
@@ -61,6 +61,8 @@ using Utils::sinc;
* DEFINES
************************************************/
+#define DP3M_RTBISECTION_ERROR 9999999
+
/************************************************
* variables
************************************************/
@@ -84,15 +86,6 @@ static void dp3m_init_a_ai_cao_cut();
*/
static bool dp3m_sanity_checks_boxl();
-/** Shift the mesh points by mesh/2 */
-static void dp3m_calc_meshift();
-
-/** Calculate the Fourier transformed differential operator.
- * Remark: This is done on the level of n-vectors and not k-vectors,
- * i.e. the prefactor i*2*PI/L is missing!
- */
-static void dp3m_calc_differential_operator();
-
/** Calculate the influence function optimized for the dipolar forces. */
static void dp3m_calc_influence_function_force();
@@ -106,21 +99,6 @@ static void dp3m_calc_influence_function_energy();
*/
static void dp3m_compute_constants_energy_dipolar();
-/** Calculate the aliasing sums for the optimal influence function.
- *
- * Calculates the aliasing sums in the nominator and denominator of
- * the expression for the optimal influence function (see
- * @cite hockney88a : 8-22, p. 275).
- *
- * \param n n-vector for which the aliasing sum is to be performed
- * \param nominator aliasing sums in the nominator
- * \retval denominator aliasing sum in the denominator
- */
-static double dp3m_perform_aliasing_sums_force(const int n[3],
- double nominator[1]);
-static double dp3m_perform_aliasing_sums_energy(const int n[3],
- double nominator[1]);
-
static double dp3m_k_space_error(double box_size, double prefac, int mesh,
int cao, int n_c_part, double sum_q2,
double alpha_L);
@@ -134,12 +112,8 @@ static double calc_surface_term(bool force_flag, bool energy_flag,
/************************************************************/
/*@{*/
-// These 3 functions are to tune the P3M code in the case of dipolar
-// interactions
-
-double P3M_DIPOLAR_real_space_error(double box_size, double prefac,
- double r_cut_iL, int n_c_part,
- double sum_q2, double alpha_L);
+double dp3m_real_space_error(double box_size, double prefac, double r_cut_iL,
+ int n_c_part, double sum_q2, double alpha_L);
static void dp3m_tune_aliasing_sums(int nx, int ny, int nz, int mesh,
double mesh_i, int cao, double alpha_L_i,
double *alias1, double *alias2);
@@ -153,17 +127,23 @@ double dp3m_rtbisection(double box_size, double prefac, double r_cut_iL,
/*@}*/
-/************************************************************/
-/* functions related to the correction of the dipolar p3m-energy */
+/** Correction of the dipolar p3m-energy. */
+double dp3m_average_dipolar_self_energy() {
+ auto const start = Utils::Vector3i{dp3m.fft.plan[3].start};
+ auto const size = Utils::Vector3i{dp3m.fft.plan[3].new_mesh};
-static double dp3m_average_dipolar_self_energy(double box_l, int mesh);
-static double dp3m_perform_aliasing_sums_dipolar_self_energy(const int n[3]);
+ auto const node_phi = grid_influence_function_self_energy(
+ dp3m.params, start, start + size, dp3m.g_energy);
+
+ double phi = 0.0;
+ MPI_Reduce(&node_phi, &phi, 1, MPI_DOUBLE, MPI_SUM, 0, comm_cart);
+ phi /= 3. * box_geo.length()[0] * Utils::int_pow<3>(dp3m.params.mesh[0]);
+ return phi * Utils::pi();
+}
/************************************************************/
dp3m_data_struct::dp3m_data_struct() {
- params.epsilon = P3M_EPSILON_MAGNETIC;
-
/* local_mesh is uninitialized */
/* sm is uninitialized */
sum_dip_part = 0;
@@ -191,107 +171,39 @@ void dp3m_init() {
// dipolar prefactor is zero: magnetostatics switched off
dp3m.params.r_cut = 0.0;
dp3m.params.r_cut_iL = 0.0;
- } else {
- if (dp3m_sanity_checks(node_grid))
- return;
-
- dp3m.params.cao3 = Utils::int_pow<3>(dp3m.params.cao);
-
- /* initializes the (inverse) mesh constant dp3m.params.a (dp3m.params.ai)
- * and the cutoff for charge assignment dp3m.params.cao_cut */
- dp3m_init_a_ai_cao_cut();
-
- p3m_calc_local_ca_mesh(dp3m.local_mesh, dp3m.params, local_geo, skin);
-
- dp3m.sm.resize(comm_cart, dp3m.local_mesh);
-
- /* fix box length dependent constants */
- dp3m_scaleby_box_l();
+ return;
+ }
- int ca_mesh_size = fft_init(dp3m.local_mesh.dim, dp3m.local_mesh.margin,
- dp3m.params.mesh, dp3m.params.mesh_off,
- dp3m.ks_pnum, dp3m.fft, node_grid, comm_cart);
- dp3m.rs_mesh.resize(ca_mesh_size);
- dp3m.ks_mesh.resize(ca_mesh_size);
+ if (dp3m_sanity_checks(node_grid)) {
+ return;
+ }
- for (auto &val : dp3m.rs_mesh_dip) {
- val.resize(ca_mesh_size);
- }
+ dp3m.params.cao3 = Utils::int_pow<3>(dp3m.params.cao);
- /* k-space part: */
+ /* initializes the (inverse) mesh constant dp3m.params.a (dp3m.params.ai)
+ * and the cutoff for charge assignment dp3m.params.cao_cut */
+ dp3m_init_a_ai_cao_cut();
- dp3m_calc_differential_operator();
+ p3m_calc_local_ca_mesh(dp3m.local_mesh, dp3m.params, local_geo, skin);
- dp3m_calc_influence_function_force();
- dp3m_calc_influence_function_energy();
+ dp3m.sm.resize(comm_cart, dp3m.local_mesh);
- dp3m_count_magnetic_particles();
- }
-}
+ int ca_mesh_size = fft_init(dp3m.local_mesh.dim, dp3m.local_mesh.margin,
+ dp3m.params.mesh, dp3m.params.mesh_off,
+ dp3m.ks_pnum, dp3m.fft, node_grid, comm_cart);
+ dp3m.rs_mesh.resize(ca_mesh_size);
+ dp3m.ks_mesh.resize(ca_mesh_size);
-double dp3m_average_dipolar_self_energy(double box_l, int mesh) {
- int end[3];
- int size = 1;
- for (int i = 0; i < 3; i++) {
- size *= dp3m.fft.plan[3].new_mesh[i];
- end[i] = dp3m.fft.plan[3].start[i] + dp3m.fft.plan[3].new_mesh[i];
+ for (auto &val : dp3m.rs_mesh_dip) {
+ val.resize(ca_mesh_size);
}
- int n[3];
- double node_phi = 0.0;
- for (n[0] = dp3m.fft.plan[3].start[0]; n[0] < end[0]; n[0]++) {
- for (n[1] = dp3m.fft.plan[3].start[1]; n[1] < end[1]; n[1]++) {
- for (n[2] = dp3m.fft.plan[3].start[2]; n[2] < end[2]; n[2]++) {
- int const ind = (n[2] - dp3m.fft.plan[3].start[2]) +
- dp3m.fft.plan[3].new_mesh[2] *
- ((n[1] - dp3m.fft.plan[3].start[1]) +
- (dp3m.fft.plan[3].new_mesh[1] *
- (n[0] - dp3m.fft.plan[3].start[0])));
-
- if (((n[0] == 0) && (n[1] == 0) && (n[2] == 0)) ||
- ((n[0] % (dp3m.params.mesh[0] / 2) == 0) &&
- (n[1] % (dp3m.params.mesh[0] / 2) == 0) &&
- (n[2] % (dp3m.params.mesh[0] / 2) == 0))) {
- node_phi += 0.0;
- } else {
- double const U2 = dp3m_perform_aliasing_sums_dipolar_self_energy(n);
- node_phi +=
- dp3m.g_energy[ind] * U2 *
- (Utils::sqr(dp3m.d_op[n[0]]) + Utils::sqr(dp3m.d_op[n[1]]) +
- Utils::sqr(dp3m.d_op[n[2]]));
- }
- }
- }
- }
+ dp3m.calc_differential_operator();
- double phi = 0.0;
- MPI_Reduce(&node_phi, &phi, 1, MPI_DOUBLE, MPI_SUM, 0, comm_cart);
- phi *= Utils::pi() / 3. / box_l / pow(mesh, 3.0);
- return phi;
-}
+ /* fix box length dependent constants */
+ dp3m_scaleby_box_l();
-double dp3m_perform_aliasing_sums_dipolar_self_energy(const int n[3]) {
- double u_sum = 0.0;
- /* lots of temporary variables... */
- double f1, sx, sy, sz, mx, my, mz, nmx, nmy, nmz;
- int limit = P3M_BRILLOUIN + 5;
-
- f1 = 1.0 / (double)dp3m.params.mesh[0];
-
- for (mx = -limit; mx <= limit; mx++) {
- nmx = dp3m.meshift[n[0]] + dp3m.params.mesh[0] * mx;
- sx = pow(sinc(f1 * nmx), 2.0 * dp3m.params.cao);
- for (my = -limit; my <= limit; my++) {
- nmy = dp3m.meshift[n[1]] + dp3m.params.mesh[0] * my;
- sy = sx * pow(sinc(f1 * nmy), 2.0 * dp3m.params.cao);
- for (mz = -limit; mz <= limit; mz++) {
- nmz = dp3m.meshift[n[2]] + dp3m.params.mesh[0] * mz;
- sz = sy * pow(sinc(f1 * nmz), 2.0 * dp3m.params.cao);
- u_sum += sz;
- }
- }
- }
- return u_sum;
+ dp3m_count_magnetic_particles();
}
/******************
@@ -299,10 +211,10 @@ double dp3m_perform_aliasing_sums_dipolar_self_energy(const int n[3]) {
******************/
void dp3m_set_tune_params(double r_cut, int mesh, int cao, double alpha,
- double accuracy, int n_interpol) {
+ double accuracy) {
if (r_cut >= 0) {
dp3m.params.r_cut = r_cut;
- dp3m.params.r_cut_iL = r_cut * (1. / box_geo.length()[0]);
+ dp3m.params.r_cut_iL = r_cut / box_geo.length()[0];
}
if (mesh >= 0)
@@ -337,7 +249,7 @@ int dp3m_set_params(double r_cut, int mesh, int cao, double alpha,
return -3;
dp3m.params.r_cut = r_cut;
- dp3m.params.r_cut_iL = r_cut * (1. / box_geo.length()[0]);
+ dp3m.params.r_cut_iL = r_cut / box_geo.length()[0];
dp3m.params.mesh[2] = dp3m.params.mesh[1] = dp3m.params.mesh[0] = mesh;
dp3m.params.cao = cao;
@@ -467,8 +379,7 @@ double dp3m_calc_kspace_forces(bool force_flag, bool energy_flag,
double tmp0, tmp1;
auto const dipole_prefac =
- dipole.prefactor /
- (double)(dp3m.params.mesh[0] * dp3m.params.mesh[1] * dp3m.params.mesh[2]);
+ dipole.prefactor / Utils::int_pow<3>(dp3m.params.mesh[0]);
if (dp3m.sum_mu2 > 0) {
/* Gather information for FFT grid inside the nodes domain (inner local
@@ -504,18 +415,20 @@ double dp3m_calc_kspace_forces(bool force_flag, bool energy_flag,
for (j[2] = 0; j[2] < dp3m.fft.plan[3].new_mesh[2]; j[2]++) {
node_k_space_energy_dip +=
dp3m.g_energy[i] *
- (Utils::sqr(dp3m.rs_mesh_dip[0][ind] *
- dp3m.d_op[j[2] + dp3m.fft.plan[3].start[2]] +
- dp3m.rs_mesh_dip[1][ind] *
- dp3m.d_op[j[0] + dp3m.fft.plan[3].start[0]] +
- dp3m.rs_mesh_dip[2][ind] *
- dp3m.d_op[j[1] + dp3m.fft.plan[3].start[1]]) +
- Utils::sqr(dp3m.rs_mesh_dip[0][ind + 1] *
- dp3m.d_op[j[2] + dp3m.fft.plan[3].start[2]] +
- dp3m.rs_mesh_dip[1][ind + 1] *
- dp3m.d_op[j[0] + dp3m.fft.plan[3].start[0]] +
- dp3m.rs_mesh_dip[2][ind + 1] *
- dp3m.d_op[j[1] + dp3m.fft.plan[3].start[1]]));
+ (Utils::sqr(
+ dp3m.rs_mesh_dip[0][ind] *
+ dp3m.d_op[0][j[2] + dp3m.fft.plan[3].start[2]] +
+ dp3m.rs_mesh_dip[1][ind] *
+ dp3m.d_op[0][j[0] + dp3m.fft.plan[3].start[0]] +
+ dp3m.rs_mesh_dip[2][ind] *
+ dp3m.d_op[0][j[1] + dp3m.fft.plan[3].start[1]]) +
+ Utils::sqr(
+ dp3m.rs_mesh_dip[0][ind + 1] *
+ dp3m.d_op[0][j[2] + dp3m.fft.plan[3].start[2]] +
+ dp3m.rs_mesh_dip[1][ind + 1] *
+ dp3m.d_op[0][j[0] + dp3m.fft.plan[3].start[0]] +
+ dp3m.rs_mesh_dip[2][ind + 1] *
+ dp3m.d_op[0][j[1] + dp3m.fft.plan[3].start[1]]));
ind += 2;
i++;
}
@@ -533,7 +446,7 @@ double dp3m_calc_kspace_forces(bool force_flag, bool energy_flag,
k_space_energy_dip -=
dipole.prefactor *
(dp3m.sum_mu2 * 2 *
- pow(dp3m.params.alpha_L * (1. / box_geo.length()[0]), 3) *
+ pow(dp3m.params.alpha_L / box_geo.length()[0], 3) *
Utils::sqrt_pi_i() / 3.0);
auto const volume = box_geo.volume();
@@ -563,18 +476,18 @@ double dp3m_calc_kspace_forces(bool force_flag, bool energy_flag,
// tmp0 = Re(mu)*k, tmp1 = Im(mu)*k
tmp0 = dp3m.rs_mesh_dip[0][ind] *
- dp3m.d_op[j[2] + dp3m.fft.plan[3].start[2]] +
+ dp3m.d_op[0][j[2] + dp3m.fft.plan[3].start[2]] +
dp3m.rs_mesh_dip[1][ind] *
- dp3m.d_op[j[0] + dp3m.fft.plan[3].start[0]] +
+ dp3m.d_op[0][j[0] + dp3m.fft.plan[3].start[0]] +
dp3m.rs_mesh_dip[2][ind] *
- dp3m.d_op[j[1] + dp3m.fft.plan[3].start[1]];
+ dp3m.d_op[0][j[1] + dp3m.fft.plan[3].start[1]];
tmp1 = dp3m.rs_mesh_dip[0][ind + 1] *
- dp3m.d_op[j[2] + dp3m.fft.plan[3].start[2]] +
+ dp3m.d_op[0][j[2] + dp3m.fft.plan[3].start[2]] +
dp3m.rs_mesh_dip[1][ind + 1] *
- dp3m.d_op[j[0] + dp3m.fft.plan[3].start[0]] +
+ dp3m.d_op[0][j[0] + dp3m.fft.plan[3].start[0]] +
dp3m.rs_mesh_dip[2][ind + 1] *
- dp3m.d_op[j[1] + dp3m.fft.plan[3].start[1]];
+ dp3m.d_op[0][j[1] + dp3m.fft.plan[3].start[1]];
/* the optimal influence function is the same for torques
and energy */
@@ -594,11 +507,13 @@ double dp3m_calc_kspace_forces(bool force_flag, bool energy_flag,
for (j[0] = 0; j[0] < dp3m.fft.plan[3].new_mesh[0]; j[0]++) {
for (j[1] = 0; j[1] < dp3m.fft.plan[3].new_mesh[1]; j[1]++) {
for (j[2] = 0; j[2] < dp3m.fft.plan[3].new_mesh[2]; j[2]++) {
- dp3m.rs_mesh[ind] = dp3m.d_op[j[d] + dp3m.fft.plan[3].start[d]] *
- dp3m.ks_mesh[ind];
+ dp3m.rs_mesh[ind] =
+ dp3m.d_op[0][j[d] + dp3m.fft.plan[3].start[d]] *
+ dp3m.ks_mesh[ind];
ind++;
- dp3m.rs_mesh[ind] = dp3m.d_op[j[d] + dp3m.fft.plan[3].start[d]] *
- dp3m.ks_mesh[ind];
+ dp3m.rs_mesh[ind] =
+ dp3m.d_op[0][j[d] + dp3m.fft.plan[3].start[d]] *
+ dp3m.ks_mesh[ind];
ind++;
}
}
@@ -635,17 +550,17 @@ double dp3m_calc_kspace_forces(bool force_flag, bool energy_flag,
j[2]++) { // j[2]=n_x
// tmp0 = Im(mu)*k, tmp1 = -Re(mu)*k
tmp0 = dp3m.rs_mesh_dip[0][ind + 1] *
- dp3m.d_op[j[2] + dp3m.fft.plan[3].start[2]] +
+ dp3m.d_op[0][j[2] + dp3m.fft.plan[3].start[2]] +
dp3m.rs_mesh_dip[1][ind + 1] *
- dp3m.d_op[j[0] + dp3m.fft.plan[3].start[0]] +
+ dp3m.d_op[0][j[0] + dp3m.fft.plan[3].start[0]] +
dp3m.rs_mesh_dip[2][ind + 1] *
- dp3m.d_op[j[1] + dp3m.fft.plan[3].start[1]];
+ dp3m.d_op[0][j[1] + dp3m.fft.plan[3].start[1]];
tmp1 = dp3m.rs_mesh_dip[0][ind] *
- dp3m.d_op[j[2] + dp3m.fft.plan[3].start[2]] +
+ dp3m.d_op[0][j[2] + dp3m.fft.plan[3].start[2]] +
dp3m.rs_mesh_dip[1][ind] *
- dp3m.d_op[j[0] + dp3m.fft.plan[3].start[0]] +
+ dp3m.d_op[0][j[0] + dp3m.fft.plan[3].start[0]] +
dp3m.rs_mesh_dip[2][ind] *
- dp3m.d_op[j[1] + dp3m.fft.plan[3].start[1]];
+ dp3m.d_op[0][j[1] + dp3m.fft.plan[3].start[1]];
dp3m.ks_mesh[ind] = tmp0 * dp3m.g_force[i];
dp3m.ks_mesh[ind + 1] = -tmp1 * dp3m.g_force[i];
ind += 2;
@@ -664,23 +579,23 @@ double dp3m_calc_kspace_forces(bool force_flag, bool energy_flag,
j[1]++) { // j[1]=n_z
for (j[2] = 0; j[2] < dp3m.fft.plan[3].new_mesh[2];
j[2]++) { // j[2]=n_x
- tmp0 = dp3m.d_op[j[d] + dp3m.fft.plan[3].start[d]] *
+ tmp0 = dp3m.d_op[0][j[d] + dp3m.fft.plan[3].start[d]] *
dp3m.ks_mesh[ind];
dp3m.rs_mesh_dip[0][ind] =
- dp3m.d_op[j[2] + dp3m.fft.plan[3].start[2]] * tmp0;
+ dp3m.d_op[0][j[2] + dp3m.fft.plan[3].start[2]] * tmp0;
dp3m.rs_mesh_dip[1][ind] =
- dp3m.d_op[j[0] + dp3m.fft.plan[3].start[0]] * tmp0;
+ dp3m.d_op[0][j[0] + dp3m.fft.plan[3].start[0]] * tmp0;
dp3m.rs_mesh_dip[2][ind] =
- dp3m.d_op[j[1] + dp3m.fft.plan[3].start[1]] * tmp0;
+ dp3m.d_op[0][j[1] + dp3m.fft.plan[3].start[1]] * tmp0;
ind++;
- tmp0 = dp3m.d_op[j[d] + dp3m.fft.plan[3].start[d]] *
+ tmp0 = dp3m.d_op[0][j[d] + dp3m.fft.plan[3].start[d]] *
dp3m.ks_mesh[ind];
dp3m.rs_mesh_dip[0][ind] =
- dp3m.d_op[j[2] + dp3m.fft.plan[3].start[2]] * tmp0;
+ dp3m.d_op[0][j[2] + dp3m.fft.plan[3].start[2]] * tmp0;
dp3m.rs_mesh_dip[1][ind] =
- dp3m.d_op[j[0] + dp3m.fft.plan[3].start[0]] * tmp0;
+ dp3m.d_op[0][j[0] + dp3m.fft.plan[3].start[0]] * tmp0;
dp3m.rs_mesh_dip[2][ind] =
- dp3m.d_op[j[1] + dp3m.fft.plan[3].start[1]] * tmp0;
+ dp3m.d_op[0][j[1] + dp3m.fft.plan[3].start[1]] * tmp0;
ind++;
}
}
@@ -721,7 +636,7 @@ double dp3m_calc_kspace_forces(bool force_flag, bool energy_flag,
double calc_surface_term(bool force_flag, bool energy_flag,
const ParticleRange &particles) {
- auto const pref = dipole.prefactor * 4 * M_PI / box_geo.volume() /
+ auto const pref = dipole.prefactor * 4 * Utils::pi() / box_geo.volume() /
(2 * dp3m.params.epsilon + 1);
double suma, a[3];
double en;
@@ -796,205 +711,24 @@ double calc_surface_term(bool force_flag, bool energy_flag,
/*****************************************************************************/
-void dp3m_calc_meshift() {
- int i;
- double dmesh;
- dmesh = (double)dp3m.params.mesh[0];
- dp3m.meshift.resize(dp3m.params.mesh[0]);
- for (i = 0; i < dp3m.params.mesh[0]; i++)
- dp3m.meshift[i] = i - std::round(i / dmesh) * dmesh;
-}
-
-/*****************************************************************************/
-
-void dp3m_calc_differential_operator() {
- int i;
- double dmesh;
-
- dmesh = (double)dp3m.params.mesh[0];
- dp3m.d_op.resize(dp3m.params.mesh[0]);
-
- for (i = 0; i < dp3m.params.mesh[0]; i++)
- dp3m.d_op[i] = (double)i - std::round((double)i / dmesh) * dmesh;
-
- dp3m.d_op[dp3m.params.mesh[0] / 2] = 0;
-}
-
-/*****************************************************************************/
-
void dp3m_calc_influence_function_force() {
- int end[3];
- int size = 1;
-
- dp3m_calc_meshift();
+ auto const start = Utils::Vector3i{dp3m.fft.plan[3].start};
+ auto const size = Utils::Vector3i{dp3m.fft.plan[3].new_mesh};
- for (int i = 0; i < 3; i++) {
- size *= dp3m.fft.plan[3].new_mesh[i];
- end[i] = dp3m.fft.plan[3].start[i] + dp3m.fft.plan[3].new_mesh[i];
- }
- dp3m.g_force.resize(size);
- double fak1 = Utils::int_pow<3>(dp3m.params.mesh[0]) * 2.0 /
- (box_geo.length()[0] * box_geo.length()[0]);
-
- int n[3];
- for (n[0] = dp3m.fft.plan[3].start[0]; n[0] < end[0]; n[0]++)
- for (n[1] = dp3m.fft.plan[3].start[1]; n[1] < end[1]; n[1]++)
- for (n[2] = dp3m.fft.plan[3].start[2]; n[2] < end[2]; n[2]++) {
- int const ind = (n[2] - dp3m.fft.plan[3].start[2]) +
- dp3m.fft.plan[3].new_mesh[2] *
- ((n[1] - dp3m.fft.plan[3].start[1]) +
- (dp3m.fft.plan[3].new_mesh[1] *
- (n[0] - dp3m.fft.plan[3].start[0])));
-
- if (((n[0] == 0) && (n[1] == 0) && (n[2] == 0)) ||
- ((n[0] % (dp3m.params.mesh[0] / 2) == 0) &&
- (n[1] % (dp3m.params.mesh[0] / 2) == 0) &&
- (n[2] % (dp3m.params.mesh[0] / 2) == 0))) {
- dp3m.g_force[ind] = 0.0;
- } else {
- double nominator[1] = {0.0};
- double denominator = dp3m_perform_aliasing_sums_force(n, nominator);
- double fak2 = nominator[0];
- fak2 /=
- pow(Utils::sqr(dp3m.d_op[n[0]]) + Utils::sqr(dp3m.d_op[n[1]]) +
- Utils::sqr(dp3m.d_op[n[2]]),
- 3) *
- Utils::sqr(denominator);
- dp3m.g_force[ind] = fak1 * fak2;
- }
- }
+ dp3m.g_force = grid_influence_function<3>(dp3m.params, start, start + size,
+ box_geo.length());
}
-/*****************************************************************************/
-
-double dp3m_perform_aliasing_sums_force(const int n[3], double nominator[1]) {
- double denominator = 0.0;
- /* lots of temporary variables... */
- double sx, sy, sz, f1, f2, f3, mx, my, mz, nmx, nmy, nmz, nm2, expo;
- double limit = 30;
- double n_nm;
- double n_nm3;
-
- nominator[0] = 0.0;
-
- f1 = 1.0 / (double)dp3m.params.mesh[0];
- f2 = Utils::sqr(Utils::pi() / (dp3m.params.alpha_L));
-
- for (mx = -P3M_BRILLOUIN; mx <= P3M_BRILLOUIN; mx++) {
- nmx = dp3m.meshift[n[0]] + dp3m.params.mesh[0] * mx;
- sx = pow(sinc(f1 * nmx), 2.0 * dp3m.params.cao);
- for (my = -P3M_BRILLOUIN; my <= P3M_BRILLOUIN; my++) {
- nmy = dp3m.meshift[n[1]] + dp3m.params.mesh[0] * my;
- sy = sx * pow(sinc(f1 * nmy), 2.0 * dp3m.params.cao);
- for (mz = -P3M_BRILLOUIN; mz <= P3M_BRILLOUIN; mz++) {
- nmz = dp3m.meshift[n[2]] + dp3m.params.mesh[0] * mz;
- sz = sy * pow(sinc(f1 * nmz), 2.0 * dp3m.params.cao);
-
- nm2 = Utils::sqr(nmx) + Utils::sqr(nmy) + Utils::sqr(nmz);
- expo = f2 * nm2;
- f3 = (expo < limit) ? sz * exp(-expo) / nm2 : 0.0;
-
- n_nm = dp3m.d_op[n[0]] * nmx + dp3m.d_op[n[1]] * nmy +
- dp3m.d_op[n[2]] * nmz;
- n_nm3 = Utils::int_pow<3>(n_nm);
-
- nominator[0] += f3 * n_nm3;
- denominator += sz;
- }
- }
- }
- return denominator;
-}
-
-/*****************************************************************************/
-
void dp3m_calc_influence_function_energy() {
- int end[3];
- int size = 1;
-
- dp3m_calc_meshift();
-
- for (int i = 0; i < 3; i++) {
- size *= dp3m.fft.plan[3].new_mesh[i];
- end[i] = dp3m.fft.plan[3].start[i] + dp3m.fft.plan[3].new_mesh[i];
- }
- dp3m.g_energy.resize(size);
- double fak1 = Utils::int_pow<3>(dp3m.params.mesh[0]) * 2.0 /
- (box_geo.length()[0] * box_geo.length()[0]);
-
- int n[3];
- for (n[0] = dp3m.fft.plan[3].start[0]; n[0] < end[0]; n[0]++)
- for (n[1] = dp3m.fft.plan[3].start[1]; n[1] < end[1]; n[1]++)
- for (n[2] = dp3m.fft.plan[3].start[2]; n[2] < end[2]; n[2]++) {
- int const ind = (n[2] - dp3m.fft.plan[3].start[2]) +
- dp3m.fft.plan[3].new_mesh[2] *
- ((n[1] - dp3m.fft.plan[3].start[1]) +
- (dp3m.fft.plan[3].new_mesh[1] *
- (n[0] - dp3m.fft.plan[3].start[0])));
-
- if (((n[0] == 0) && (n[1] == 0) && (n[2] == 0)) ||
- ((n[0] % (dp3m.params.mesh[0] / 2) == 0) &&
- (n[1] % (dp3m.params.mesh[0] / 2) == 0) &&
- (n[2] % (dp3m.params.mesh[0] / 2) == 0))) {
- dp3m.g_energy[ind] = 0.0;
- } else {
- double nominator[1] = {0.0};
- double denominator = dp3m_perform_aliasing_sums_energy(n, nominator);
- double fak2 = nominator[0];
- fak2 /=
- pow(Utils::sqr(dp3m.d_op[n[0]]) + Utils::sqr(dp3m.d_op[n[1]]) +
- Utils::sqr(dp3m.d_op[n[2]]),
- 2) *
- Utils::sqr(denominator);
- dp3m.g_energy[ind] = fak1 * fak2;
- }
- }
-}
-
-/*****************************************************************************/
+ auto const start = Utils::Vector3i{dp3m.fft.plan[3].start};
+ auto const size = Utils::Vector3i{dp3m.fft.plan[3].new_mesh};
-double dp3m_perform_aliasing_sums_energy(const int n[3], double nominator[1]) {
- double denominator = 0.0;
- /* lots of temporary variables... */
- double sx, sy, sz, f1, f2, f3, mx, my, mz, nmx, nmy, nmz, nm2, expo;
- double limit = 30;
- double n_nm;
- double n_nm2;
-
- nominator[0] = 0.0;
-
- f1 = 1.0 / (double)dp3m.params.mesh[0];
- f2 = Utils::sqr(Utils::pi() / (dp3m.params.alpha_L));
-
- for (mx = -P3M_BRILLOUIN; mx <= P3M_BRILLOUIN; mx++) {
- nmx = dp3m.meshift[n[0]] + dp3m.params.mesh[0] * mx;
- sx = pow(sinc(f1 * nmx), 2.0 * dp3m.params.cao);
- for (my = -P3M_BRILLOUIN; my <= P3M_BRILLOUIN; my++) {
- nmy = dp3m.meshift[n[1]] + dp3m.params.mesh[0] * my;
- sy = sx * pow(sinc(f1 * nmy), 2.0 * dp3m.params.cao);
- for (mz = -P3M_BRILLOUIN; mz <= P3M_BRILLOUIN; mz++) {
- nmz = dp3m.meshift[n[2]] + dp3m.params.mesh[0] * mz;
- sz = sy * pow(sinc(f1 * nmz), 2.0 * dp3m.params.cao);
-
- nm2 = Utils::sqr(nmx) + Utils::sqr(nmy) + Utils::sqr(nmz);
- expo = f2 * nm2;
- f3 = (expo < limit) ? sz * exp(-expo) / nm2 : 0.0;
-
- n_nm = dp3m.d_op[n[0]] * nmx + dp3m.d_op[n[1]] * nmy +
- dp3m.d_op[n[2]] * nmz;
- n_nm2 = n_nm * n_nm;
- nominator[0] += f3 * n_nm2;
- denominator += sz;
- }
- }
- }
- return denominator;
+ dp3m.g_energy = grid_influence_function<2>(dp3m.params, start, start + size,
+ box_geo.length());
}
/*****************************************************************************/
-#define P3M_TUNE_MAX_CUTS 50
-
/** @copybrief p3m_get_accuracy
*
* The real space error is tuned such that it contributes half of the
@@ -1017,17 +751,21 @@ double dp3m_get_accuracy(int mesh, int cao, double r_cut_iL, double *_alpha_L,
// Alpha cannot be zero in the dipolar case because real_space formula breaks
// down
- rs_err = P3M_DIPOLAR_real_space_error(box_geo.length()[0], dipole.prefactor,
- r_cut_iL, dp3m.sum_dip_part,
- dp3m.sum_mu2, 0.001);
+ rs_err =
+ dp3m_real_space_error(box_geo.length()[0], dipole.prefactor, r_cut_iL,
+ dp3m.sum_dip_part, dp3m.sum_mu2, 0.001);
- if (M_SQRT2 * rs_err > dp3m.params.accuracy) {
+ if (Utils::sqrt_2() * rs_err > dp3m.params.accuracy) {
/* assume rs_err = ks_err -> rs_err = accuracy/sqrt(2.0) -> alpha_L */
alpha_L = dp3m_rtbisection(
box_geo.length()[0], dipole.prefactor, r_cut_iL, dp3m.sum_dip_part,
dp3m.sum_mu2, 0.0001 * box_geo.length()[0], 5.0 * box_geo.length()[0],
0.0001, dp3m.params.accuracy);
-
+ if (alpha_L == -DP3M_RTBISECTION_ERROR) {
+ *_rs_err = -1;
+ *_ks_err = -1;
+ return -DP3M_RTBISECTION_ERROR;
+ }
}
else
@@ -1039,9 +777,9 @@ double dp3m_get_accuracy(int mesh, int cao, double r_cut_iL, double *_alpha_L,
*_alpha_L = alpha_L;
/* calculate real space and k-space error for this alpha_L */
- rs_err = P3M_DIPOLAR_real_space_error(box_geo.length()[0], dipole.prefactor,
- r_cut_iL, dp3m.sum_dip_part,
- dp3m.sum_mu2, alpha_L);
+ rs_err =
+ dp3m_real_space_error(box_geo.length()[0], dipole.prefactor, r_cut_iL,
+ dp3m.sum_dip_part, dp3m.sum_mu2, alpha_L);
ks_err = dp3m_k_space_error(box_geo.length()[0], dipole.prefactor, mesh, cao,
dp3m.sum_dip_part, dp3m.sum_mu2, alpha_L);
@@ -1072,7 +810,6 @@ static double dp3m_mcr_time(int mesh, int cao, double r_cut_iL,
dp3m.params.mesh[0] = dp3m.params.mesh[1] = dp3m.params.mesh[2] = mesh;
dp3m.params.cao = cao;
dp3m.params.alpha_L = alpha_L;
- dp3m_scaleby_box_l();
/* initialize p3m structures */
mpi_bcast_coulomb_params();
/* perform force calculation test */
@@ -1104,15 +841,13 @@ static double dp3m_mcr_time(int mesh, int cao, double r_cut_iL,
static double dp3m_mc_time(char **log, int mesh, int cao, double r_cut_iL_min,
double r_cut_iL_max, double *_r_cut_iL,
double *_alpha_L, double *_accuracy) {
- double int_time;
double r_cut_iL;
- double rs_err, ks_err, mesh_size, k_cut;
- int i, n_cells;
+ double rs_err, ks_err;
char b[3 * ES_INTEGER_SPACE + 3 * ES_DOUBLE_SPACE + 128];
/* initial checks. */
- mesh_size = box_geo.length()[0] / (double)mesh;
- k_cut = mesh_size * cao / 2.0;
+ auto const mesh_size = box_geo.length()[0] / static_cast(mesh);
+ auto const k_cut = mesh_size * cao / 2.0;
auto const min_box_l = *boost::min_element(box_geo.length());
auto const min_local_box_l = *boost::min_element(local_geo.length());
@@ -1128,9 +863,11 @@ static double dp3m_mc_time(char **log, int mesh, int cao, double r_cut_iL_min,
is initially 0 and therefore
has infinite error estimate, as required. Therefore if the high boundary
fails, there is no possible r_cut */
- if ((*_accuracy = dp3m_get_accuracy(mesh, cao, r_cut_iL_max, _alpha_L,
- &rs_err, &ks_err)) >
- dp3m.params.accuracy) {
+ *_accuracy =
+ dp3m_get_accuracy(mesh, cao, r_cut_iL_max, _alpha_L, &rs_err, &ks_err);
+ if (*_accuracy == -DP3M_RTBISECTION_ERROR)
+ return *_accuracy;
+ if (*_accuracy > dp3m.params.accuracy) {
/* print result */
sprintf(b, "%-4d %-3d %.5e %.5e %.5e %.3e %.3e accuracy not achieved\n",
mesh, cao, r_cut_iL_max, *_alpha_L, *_accuracy, rs_err, ks_err);
@@ -1145,8 +882,11 @@ static double dp3m_mc_time(char **log, int mesh, int cao, double r_cut_iL_min,
break;
/* bisection */
- if (dp3m_get_accuracy(mesh, cao, r_cut_iL, _alpha_L, &rs_err, &ks_err) >
- dp3m.params.accuracy)
+ auto const tmp_accuracy =
+ dp3m_get_accuracy(mesh, cao, r_cut_iL, _alpha_L, &rs_err, &ks_err);
+ if (tmp_accuracy == -DP3M_RTBISECTION_ERROR)
+ return tmp_accuracy;
+ if (tmp_accuracy > dp3m.params.accuracy)
r_cut_iL_min = r_cut_iL;
else
r_cut_iL_max = r_cut_iL;
@@ -1161,7 +901,7 @@ static double dp3m_mc_time(char **log, int mesh, int cao, double r_cut_iL_min,
runtimeErrorMsg() << "dipolar P3M: tuning when dlc needs to be fixed";
}
- int_time = dp3m_mcr_time(mesh, cao, r_cut_iL, *_alpha_L);
+ double const int_time = dp3m_mcr_time(mesh, cao, r_cut_iL, *_alpha_L);
if (int_time == -P3M_TUNE_FAIL) {
*log = strcat_alloc(*log, "tuning failed, test integration not possible\n");
return int_time;
@@ -1169,10 +909,13 @@ static double dp3m_mc_time(char **log, int mesh, int cao, double r_cut_iL_min,
*_accuracy =
dp3m_get_accuracy(mesh, cao, r_cut_iL, _alpha_L, &rs_err, &ks_err);
+ if (*_accuracy == -DP3M_RTBISECTION_ERROR) {
+ return *_accuracy;
+ }
/* print result */
- sprintf(b, "%-4d %-3d %.5e %.5e %.5e %.3e %.3e %-8d\n", mesh, cao, r_cut_iL,
- *_alpha_L, *_accuracy, rs_err, ks_err, (int)int_time);
+ sprintf(b, "%-4d %-3d %.5e %.5e %.5e %.3e %.3e %-8.0f\n", mesh, cao, r_cut_iL,
+ *_alpha_L, *_accuracy, rs_err, ks_err, int_time);
*log = strcat_alloc(*log, b);
return int_time;
}
@@ -1216,7 +959,7 @@ static double dp3m_m_time(char **log, int mesh, int cao_min, int cao_max,
tmp_time = dp3m_mc_time(log, mesh, cao, r_cut_iL_min, r_cut_iL_max,
&tmp_r_cut_iL, &tmp_alpha_L, &tmp_accuracy);
/* bail out if the force evaluation is not working */
- if (tmp_time == -P3M_TUNE_FAIL)
+ if (tmp_time == -P3M_TUNE_FAIL || tmp_time == -DP3M_RTBISECTION_ERROR)
return tmp_time;
/* cao is too large for this grid, but still the accuracy cannot be
* achieved, give up */
@@ -1255,7 +998,7 @@ static double dp3m_m_time(char **log, int mesh, int cao_min, int cao_max,
dp3m_mc_time(log, mesh, cao + final_dir, r_cut_iL_min, r_cut_iL_max,
&tmp_r_cut_iL, &tmp_alpha_L, &tmp_accuracy);
/* bail out on errors, as usual */
- if (tmp_time == -P3M_TUNE_FAIL)
+ if (tmp_time == -P3M_TUNE_FAIL || tmp_time == -DP3M_RTBISECTION_ERROR)
return tmp_time;
/* in this direction, we cannot optimise, since we get into precision
* trouble */
@@ -1306,7 +1049,7 @@ static double dp3m_m_time(char **log, int mesh, int cao_min, int cao_max,
tmp_time = dp3m_mc_time(log, mesh, cao, r_cut_iL_min, r_cut_iL_max,
&tmp_r_cut_iL, &tmp_alpha_L, &tmp_accuracy);
/* bail out on errors, as usual */
- if (tmp_time == -P3M_TUNE_FAIL)
+ if (tmp_time == -P3M_TUNE_FAIL || tmp_time == -DP3M_RTBISECTION_ERROR)
return tmp_time;
/* if we cannot meet the precision anymore, give up */
if (tmp_time < 0)
@@ -1368,8 +1111,7 @@ int dp3m_adaptive_tune(char **logger) {
*logger = strcat_alloc(*logger, b);
if (dp3m.sum_dip_part == 0) {
- *logger = strcat_alloc(
- *logger, "no dipolar particles in the system, cannot tune dipolar P3M");
+ runtimeErrorMsg() << "no dipolar particles in the system";
return ES_ERROR;
}
@@ -1399,8 +1141,8 @@ int dp3m_adaptive_tune(char **logger) {
r_cut_iL_min = 0;
r_cut_iL_max = std::min(min_local_box_l, min_box_l / 2) - skin;
- r_cut_iL_min *= (1. / box_geo.length()[0]);
- r_cut_iL_max *= (1. / box_geo.length()[0]);
+ r_cut_iL_min *= 1. / box_geo.length()[0];
+ r_cut_iL_max *= 1. / box_geo.length()[0];
} else {
r_cut_iL_min = r_cut_iL_max = dp3m.params.r_cut_iL;
@@ -1428,7 +1170,7 @@ int dp3m_adaptive_tune(char **logger) {
dp3m_m_time(logger, tmp_mesh, cao_min, cao_max, &tmp_cao, r_cut_iL_min,
r_cut_iL_max, &tmp_r_cut_iL, &tmp_alpha_L, &tmp_accuracy);
/* some error occurred during the tuning force evaluation */
- if (tmp_time == -1)
+ if (tmp_time == -P3M_TUNE_FAIL || tmp_time == -DP3M_RTBISECTION_ERROR)
return ES_ERROR;
/* this mesh does not work at all */
if (tmp_time < 0)
@@ -1453,9 +1195,7 @@ int dp3m_adaptive_tune(char **logger) {
}
if (time_best == 1e20) {
- *logger = strcat_alloc(
- *logger,
- "failed to tune dipolar P3M parameters to required accuracy\n");
+ runtimeErrorMsg() << "failed to reach requested accuracy";
return ES_ERROR;
}
@@ -1465,7 +1205,6 @@ int dp3m_adaptive_tune(char **logger) {
dp3m.params.cao = cao;
dp3m.params.alpha_L = alpha_L;
dp3m.params.accuracy = accuracy;
- dp3m_scaleby_box_l();
/* broadcast tuned p3m parameters */
mpi_bcast_coulomb_params();
/* Tell the user about the outcome */
@@ -1478,12 +1217,8 @@ int dp3m_adaptive_tune(char **logger) {
}
void dp3m_count_magnetic_particles() {
- double node_sums[2], tot_sums[2];
-
- for (int i = 0; i < 2; i++) {
- node_sums[i] = 0.0;
- tot_sums[i] = 0.0;
- }
+ double node_sums[2] = {0, 0};
+ double tot_sums[2] = {0, 0};
for (auto const &p : cell_structure.local_particles()) {
if (p.p.dipm != 0.0) {
@@ -1503,57 +1238,58 @@ REGISTER_CALLBACK(dp3m_count_magnetic_particles)
static double dp3m_k_space_error(double box_size, double prefac, int mesh,
int cao, int n_c_part, double sum_q2,
double alpha_L) {
- int nx, ny, nz;
- double he_q = 0.0, mesh_i = 1. / mesh, alpha_L_i = 1. / alpha_L;
- double alias1, alias2, n2, cs;
+ double he_q = 0.0;
+ auto const mesh_i = 1. / mesh;
+ auto const alpha_L_i = 1. / alpha_L;
- for (nx = -mesh / 2; nx < mesh / 2; nx++)
- for (ny = -mesh / 2; ny < mesh / 2; ny++)
- for (nz = -mesh / 2; nz < mesh / 2; nz++)
+ for (int nx = -mesh / 2; nx < mesh / 2; nx++)
+ for (int ny = -mesh / 2; ny < mesh / 2; ny++)
+ for (int nz = -mesh / 2; nz < mesh / 2; nz++)
if ((nx != 0) || (ny != 0) || (nz != 0)) {
- n2 = Utils::sqr(nx) + Utils::sqr(ny) + Utils::sqr(nz);
- cs = p3m_analytic_cotangent_sum(nx, mesh_i, cao) *
- p3m_analytic_cotangent_sum(ny, mesh_i, cao) *
- p3m_analytic_cotangent_sum(nz, mesh_i, cao);
+ auto const n2 = Utils::sqr(nx) + Utils::sqr(ny) + Utils::sqr(nz);
+ auto const cs = p3m_analytic_cotangent_sum(nx, mesh_i, cao) *
+ p3m_analytic_cotangent_sum(ny, mesh_i, cao) *
+ p3m_analytic_cotangent_sum(nz, mesh_i, cao);
+ double alias1, alias2;
dp3m_tune_aliasing_sums(nx, ny, nz, mesh, mesh_i, cao, alpha_L_i,
&alias1, &alias2);
- double d = alias1 - Utils::sqr(alias2 / cs) / Utils::int_pow<3>(n2);
+ double d = alias1 - Utils::sqr(alias2 / cs) /
+ Utils::int_pow<3>(static_cast(n2));
/* at high precisions, d can become negative due to extinction;
also, don't take values that have no significant digits left*/
if (d > 0 && (fabs(d / alias1) > ROUND_ERROR_PREC))
he_q += d;
}
- return 8. * Utils::pi() * Utils::pi() / 3. * sum_q2 *
- sqrt(he_q / (double)n_c_part) / Utils::int_pow<4>(box_size);
+ return 8. * Utils::sqr(Utils::pi()) / 3. * sum_q2 * sqrt(he_q / n_c_part) /
+ Utils::int_pow<4>(box_size);
}
/** Tuning dipolar-P3M */
void dp3m_tune_aliasing_sums(int nx, int ny, int nz, int mesh, double mesh_i,
int cao, double alpha_L_i, double *alias1,
double *alias2) {
- int mx, my, mz;
- double nmx, nmy, nmz;
- double fnmx, fnmy, fnmz;
+ using Utils::sinc;
- double ex, ex2, nm2, U2, factor1;
-
- factor1 = Utils::sqr(Utils::pi() * alpha_L_i);
+ auto const factor1 = Utils::sqr(Utils::pi() * alpha_L_i);
*alias1 = *alias2 = 0.0;
- for (mx = -P3M_BRILLOUIN; mx <= P3M_BRILLOUIN; mx++) {
- fnmx = mesh_i * (nmx = nx + mx * mesh);
- for (my = -P3M_BRILLOUIN; my <= P3M_BRILLOUIN; my++) {
- fnmy = mesh_i * (nmy = ny + my * mesh);
- for (mz = -P3M_BRILLOUIN; mz <= P3M_BRILLOUIN; mz++) {
- fnmz = mesh_i * (nmz = nz + mz * mesh);
-
- nm2 = Utils::sqr(nmx) + Utils::sqr(nmy) + Utils::sqr(nmz);
- ex2 = Utils::sqr(ex = exp(-factor1 * nm2));
-
- U2 = pow(sinc(fnmx) * sinc(fnmy) * sinc(fnmz), 2.0 * cao);
-
- *alias1 += ex2 * nm2;
+ for (int mx = -P3M_BRILLOUIN; mx <= P3M_BRILLOUIN; mx++) {
+ auto const nmx = nx + mx * mesh;
+ auto const fnmx = mesh_i * nmx;
+ for (int my = -P3M_BRILLOUIN; my <= P3M_BRILLOUIN; my++) {
+ auto const nmy = ny + my * mesh;
+ auto const fnmy = mesh_i * nmy;
+ for (int mz = -P3M_BRILLOUIN; mz <= P3M_BRILLOUIN; mz++) {
+ auto const nmz = nz + mz * mesh;
+ auto const fnmz = mesh_i * nmz;
+
+ auto const nm2 = Utils::sqr(nmx) + Utils::sqr(nmy) + Utils::sqr(nmz);
+ auto const ex = std::exp(-factor1 * nm2);
+
+ auto const U2 = pow(sinc(fnmx) * sinc(fnmy) * sinc(fnmz), 2.0 * cao);
+
+ *alias1 += Utils::sqr(ex) * nm2;
*alias2 += U2 * ex * pow((nx * nmx + ny * nmy + nz * nmz), 3.) / nm2;
}
}
@@ -1566,9 +1302,8 @@ void dp3m_tune_aliasing_sums(int nx, int ny, int nz, int mesh, double mesh_i,
* Please note that in this more refined approach we don't use
* eq. (37), but eq. (33) which maintains all the powers in alpha.
*/
-double P3M_DIPOLAR_real_space_error(double box_size, double prefac,
- double r_cut_iL, int n_c_part,
- double sum_q2, double alpha_L) {
+double dp3m_real_space_error(double box_size, double prefac, double r_cut_iL,
+ int n_c_part, double sum_q2, double alpha_L) {
double d_error_f, d_cc, d_dc, d_rcut2, d_con;
double d_a2, d_c, d_RCUT;
@@ -1603,41 +1338,42 @@ double P3M_DIPOLAR_real_space_error(double box_size, double prefac,
double dp3m_rtbisection(double box_size, double prefac, double r_cut_iL,
int n_c_part, double sum_q2, double x1, double x2,
double xacc, double tuned_accuracy) {
- int j;
- double dx, f, fmid, xmid, rtb, constant, JJ_RTBIS_MAX = 40;
+ constexpr int JJ_RTBIS_MAX = 40;
- constant = tuned_accuracy / sqrt(2.);
+ auto const constant = tuned_accuracy / Utils::sqrt_2();
- f = P3M_DIPOLAR_real_space_error(box_size, prefac, r_cut_iL, n_c_part, sum_q2,
- x1) -
+ auto const f1 =
+ dp3m_real_space_error(box_size, prefac, r_cut_iL, n_c_part, sum_q2, x1) -
+ constant;
+ auto const f2 =
+ dp3m_real_space_error(box_size, prefac, r_cut_iL, n_c_part, sum_q2, x2) -
constant;
- fmid = P3M_DIPOLAR_real_space_error(box_size, prefac, r_cut_iL, n_c_part,
- sum_q2, x2) -
- constant;
- if (f * fmid >= 0.0)
- fprintf(stderr,
- "Root must be bracketed for bisection in dp3m_rtbisection\n");
- rtb = f < 0.0 ? (dx = x2 - x1, x1)
- : (dx = x1 - x2,
- x2); // Orient the search dx, and set rtb to x1 or x2 ...
- for (j = 1; j <= JJ_RTBIS_MAX; j++) {
- fmid = P3M_DIPOLAR_real_space_error(box_size, prefac, r_cut_iL, n_c_part,
- sum_q2, xmid = rtb + (dx *= 0.5)) -
- constant;
+ if (f1 * f2 >= 0.0) {
+ runtimeErrorMsg()
+ << "Root must be bracketed for bisection in dp3m_rtbisection";
+ return -DP3M_RTBISECTION_ERROR;
+ }
+ // Orient the search dx, and set rtb to x1 or x2 ...
+ double dx;
+ double rtb = f1 < 0.0 ? (dx = x2 - x1, x1) : (dx = x1 - x2, x2);
+ for (int j = 1; j <= JJ_RTBIS_MAX; j++) {
+ auto const xmid = rtb + (dx *= 0.5);
+ auto const fmid = dp3m_real_space_error(box_size, prefac, r_cut_iL,
+ n_c_part, sum_q2, xmid) -
+ constant;
if (fmid <= 0.0)
rtb = xmid;
if (fabs(dx) < xacc || fmid == 0.0)
return rtb;
}
- fprintf(stderr, "Too many bisections in JJ_rtbissection\n");
- return -9999999.9999;
+ runtimeErrorMsg() << "Too many bisections in dp3m_rtbisection";
+ return -DP3M_RTBISECTION_ERROR;
}
/************************************************************/
void dp3m_init_a_ai_cao_cut() {
- int i;
- for (i = 0; i < 3; i++) {
+ for (int i = 0; i < 3; i++) {
dp3m.params.ai[i] = (double)dp3m.params.mesh[i] / box_geo.length()[i];
dp3m.params.a[i] = 1.0 / dp3m.params.ai[i];
dp3m.params.cao_cut[i] = 0.5 * dp3m.params.a[i] * dp3m.params.cao;
@@ -1647,9 +1383,8 @@ void dp3m_init_a_ai_cao_cut() {
/*****************************************************************************/
bool dp3m_sanity_checks_boxl() {
- int i;
bool ret = false;
- for (i = 0; i < 3; i++) {
+ for (int i = 0; i < 3; i++) {
/* check k-space cutoff */
if (dp3m.params.cao_cut[i] >= 0.5 * box_geo.length()[i]) {
runtimeErrorMsg() << "dipolar P3M_init: k-space cutoff "
@@ -1726,7 +1461,7 @@ void dp3m_scaleby_box_l() {
}
dp3m.params.r_cut = dp3m.params.r_cut_iL * box_geo.length()[0];
- dp3m.params.alpha = dp3m.params.alpha_L * (1. / box_geo.length()[0]);
+ dp3m.params.alpha = dp3m.params.alpha_L / box_geo.length()[0];
dp3m_init_a_ai_cao_cut();
p3m_calc_lm_ld_pos(dp3m.local_mesh, dp3m.params);
dp3m_sanity_checks_boxl();
@@ -1737,17 +1472,13 @@ void dp3m_scaleby_box_l() {
/** Calculate the dipolar-P3M energy correction */
void dp3m_compute_constants_energy_dipolar() {
- double Eself, Ukp3m;
if (dp3m.energy_correction != 0.0)
return;
- auto const volume = box_geo.volume();
- Ukp3m = dp3m_average_dipolar_self_energy(box_geo.length()[0],
- dp3m.params.mesh[0]) *
- volume;
+ auto const Ukp3m = dp3m_average_dipolar_self_energy() * box_geo.volume();
- Eself = -(2 * pow(dp3m.params.alpha_L, 3) * Utils::sqrt_pi_i() / 3.0);
+ auto const Eself = -2 * pow(dp3m.params.alpha_L, 3) * Utils::sqrt_pi_i() / 3;
dp3m.energy_correction =
-dp3m.sum_mu2 * (Ukp3m + Eself + 2. * Utils::pi() / 3.);
diff --git a/src/core/electrostatics_magnetostatics/p3m-dipolar.hpp b/src/core/electrostatics_magnetostatics/p3m-dipolar.hpp
index 6e7b565eb74..02b2c4d7fac 100644
--- a/src/core/electrostatics_magnetostatics/p3m-dipolar.hpp
+++ b/src/core/electrostatics_magnetostatics/p3m-dipolar.hpp
@@ -40,6 +40,7 @@
#include "electrostatics_magnetostatics/dipole.hpp"
#include "fft.hpp"
#include "p3m-common.hpp"
+#include "p3m-data_struct.hpp"
#include "p3m_interpolation.hpp"
#include "p3m_send_mesh.hpp"
@@ -47,11 +48,9 @@
#include
#include
-struct dp3m_data_struct {
+struct dp3m_data_struct : public p3m_data_struct_base {
dp3m_data_struct();
- P3MParameters params;
-
/** local mesh. */
p3m_local_mesh local_mesh;
/** real space mesh (local) for CA/FFT.*/
@@ -68,22 +67,9 @@ struct dp3m_data_struct {
/** position shift for calc. of first assignment mesh point. */
double pos_shift;
- /** help variable for calculation of aliasing sums */
- std::vector meshift;
-
- /** Spatial differential operator in k-space. We use an i*k differentiation.
- */
- std::vector d_op;
- /** Force optimised influence function (k-space) */
- std::vector g_force;
- /** Energy optimised influence function (k-space) */
- std::vector g_energy;
p3m_interpolation_cache inter_weights;
- /** number of permutations in k_space */
- int ks_pnum;
-
/** send/recv mesh sizes */
p3m_send_mesh sm;
@@ -96,13 +82,9 @@ struct dp3m_data_struct {
/** dipolar P3M parameters. */
extern dp3m_data_struct dp3m;
-/** \name Exported Functions */
-/************************************************************/
-/*@{*/
-
/** @copydoc p3m_set_tune_params */
void dp3m_set_tune_params(double r_cut, int mesh, int cao, double alpha,
- double accuracy, int n_interpol);
+ double accuracy);
/** @copydoc p3m_set_params */
int dp3m_set_params(double r_cut, int mesh, int cao, double alpha,
@@ -305,7 +287,5 @@ inline double dp3m_pair_energy(Particle const &p1, Particle const &p2,
return 0.0;
}
-/*@}*/
-
#endif /* DP3M */
#endif /* _P3M_DIPOLES_H */
diff --git a/src/core/electrostatics_magnetostatics/p3m.cpp b/src/core/electrostatics_magnetostatics/p3m.cpp
index aef4ff2d9d7..5704b379a29 100644
--- a/src/core/electrostatics_magnetostatics/p3m.cpp
+++ b/src/core/electrostatics_magnetostatics/p3m.cpp
@@ -61,21 +61,6 @@ using Utils::strcat_alloc;
p3m_data_struct p3m;
-/** @name Index helpers for direct and reciprocal space
- * After the FFT the data is in order YZX, which
- * means that Y is the slowest changing index.
- * The defines are here to not get confused and
- * be able to easily change the order.
- */
-/*@{*/
-#define RX 0
-#define RY 1
-#define RZ 2
-#define KY 0
-#define KZ 1
-#define KX 2
-/*@}*/
-
/** \name Private Functions */
/*@{*/
@@ -95,12 +80,6 @@ static bool p3m_sanity_checks_system(const Utils::Vector3i &grid);
*/
static bool p3m_sanity_checks_boxl();
-/** Calculate the Fourier transformed differential operator.
- * Remark: This is done on the level of n-vectors and not k-vectors,
- * i.e. the prefactor i*2*PI/L is missing!
- */
-static void p3m_calc_differential_operator();
-
/** Calculate the optimal influence function of @cite hockney88a.
* (optimised for force calculations)
*
@@ -156,6 +135,8 @@ static void p3m_tune_aliasing_sums(int nx, int ny, int nz, const int mesh[3],
double alpha_L_i, double *alias1,
double *alias2);
+/*@}*/
+
p3m_data_struct::p3m_data_struct() {
/* local_mesh is uninitialized */
/* sm is uninitialized */
@@ -172,36 +153,38 @@ void p3m_init() {
// prefactor is zero: electrostatics switched off
p3m.params.r_cut = 0.0;
p3m.params.r_cut_iL = 0.0;
- } else {
- if (p3m_sanity_checks()) {
- return;
- }
- p3m.params.cao3 = p3m.params.cao * p3m.params.cao * p3m.params.cao;
+ return;
+ }
- /* initializes the (inverse) mesh constant p3m.params.a (p3m.params.ai) and
- * the cutoff for charge assignment p3m.params.cao_cut */
- p3m_init_a_ai_cao_cut();
+ if (p3m_sanity_checks()) {
+ return;
+ }
- p3m_calc_local_ca_mesh(p3m.local_mesh, p3m.params, local_geo, skin);
+ p3m.params.cao3 = Utils::int_pow<3>(p3m.params.cao);
- p3m.sm.resize(comm_cart, p3m.local_mesh);
+ /* initializes the (inverse) mesh constant p3m.params.a (p3m.params.ai) and
+ * the cutoff for charge assignment p3m.params.cao_cut */
+ p3m_init_a_ai_cao_cut();
- int ca_mesh_size = fft_init(p3m.local_mesh.dim, p3m.local_mesh.margin,
- p3m.params.mesh, p3m.params.mesh_off,
- p3m.ks_pnum, p3m.fft, node_grid, comm_cart);
- p3m.rs_mesh.resize(ca_mesh_size);
- for (auto &e : p3m.E_mesh) {
- e.resize(ca_mesh_size);
- }
+ p3m_calc_local_ca_mesh(p3m.local_mesh, p3m.params, local_geo, skin);
- /* k-space part: */
- p3m_calc_differential_operator();
+ p3m.sm.resize(comm_cart, p3m.local_mesh);
- /* fix box length dependent constants */
- p3m_scaleby_box_l();
+ int ca_mesh_size =
+ fft_init(p3m.local_mesh.dim, p3m.local_mesh.margin, p3m.params.mesh,
+ p3m.params.mesh_off, p3m.ks_pnum, p3m.fft, node_grid, comm_cart);
+ p3m.rs_mesh.resize(ca_mesh_size);
- p3m_count_charged_particles();
+ for (auto &e : p3m.E_mesh) {
+ e.resize(ca_mesh_size);
}
+
+ p3m.calc_differential_operator();
+
+ /* fix box length dependent constants */
+ p3m_scaleby_box_l();
+
+ p3m_count_charged_particles();
}
void p3m_set_tune_params(double r_cut, const int mesh[3], int cao, double alpha,
@@ -229,8 +212,6 @@ void p3m_set_tune_params(double r_cut, const int mesh[3], int cao, double alpha,
p3m.params.accuracy = accuracy;
}
-/*@}*/
-
int p3m_set_params(double r_cut, const int *mesh, int cao, double alpha,
double accuracy) {
if (coulomb.method != COULOMB_P3M && coulomb.method != COULOMB_ELC_P3M &&
@@ -392,7 +373,7 @@ auto calc_dipole_moment(boost::mpi::communicator const &comm,
void add_dipole_correction(Utils::Vector3d const &box_dipole,
const ParticleRange &particles) {
- auto const pref = coulomb.prefactor * 4 * M_PI / box_geo.volume() /
+ auto const pref = coulomb.prefactor * 4 * Utils::pi() / box_geo.volume() /
(2 * p3m.params.epsilon + 1);
auto const dm = pref * box_dipole;
@@ -403,7 +384,7 @@ void add_dipole_correction(Utils::Vector3d const &box_dipole,
}
double dipole_correction_energy(Utils::Vector3d const &box_dipole) {
- auto const pref = coulomb.prefactor * 4 * M_PI / box_geo.volume() /
+ auto const pref = coulomb.prefactor * 4 * Utils::pi() / box_geo.volume() /
(2 * p3m.params.epsilon + 1);
return pref * box_dipole.norm2();
@@ -416,6 +397,8 @@ double dipole_correction_energy(Utils::Vector3d const &box_dipole) {
* eq. (2.8) is not present here since M is the empty set in our simulations.
*/
Utils::Vector9d p3m_calc_kspace_pressure_tensor() {
+ using namespace detail::FFT_indexing;
+
Utils::Vector9d node_k_space_pressure_tensor{};
if (p3m.sum_q2 > 0) {
@@ -575,19 +558,6 @@ double p3m_calc_kspace_forces(bool force_flag, bool energy_flag,
return 0.0;
}
-void p3m_calc_differential_operator() {
- for (int i = 0; i < 3; i++) {
- p3m.d_op[i].resize(p3m.params.mesh[i]);
- p3m.d_op[i][0] = 0;
- p3m.d_op[i][p3m.params.mesh[i] / 2] = 0.0;
-
- for (int j = 1; j < p3m.params.mesh[i] / 2; j++) {
- p3m.d_op[i][j] = j;
- p3m.d_op[i][p3m.params.mesh[i] - j] = -j;
- }
- }
-}
-
void p3m_calc_influence_function_force() {
auto const start = Utils::Vector3i{p3m.fft.plan[3].start};
auto const size = Utils::Vector3i{p3m.fft.plan[3].new_mesh};
@@ -604,8 +574,6 @@ void p3m_calc_influence_function_energy() {
box_geo.length());
}
-#define P3M_TUNE_MAX_CUTS 50
-
/** Get the minimal error for this combination of parameters.
*
* The real space error is tuned such that it contributes half of the
@@ -629,9 +597,10 @@ static double p3m_get_accuracy(const int mesh[3], int cao, double r_cut_iL,
rs_err = p3m_real_space_error(coulomb.prefactor, r_cut_iL, p3m.sum_qpart,
p3m.sum_q2, 0);
- if (M_SQRT2 * rs_err > p3m.params.accuracy) {
+ if (Utils::sqrt_2() * rs_err > p3m.params.accuracy) {
/* assume rs_err = ks_err -> rs_err = accuracy/sqrt(2.0) -> alpha_L */
- alpha_L = sqrt(log(M_SQRT2 * rs_err / p3m.params.accuracy)) / r_cut_iL;
+ alpha_L =
+ sqrt(log(Utils::sqrt_2() * rs_err / p3m.params.accuracy)) / r_cut_iL;
} else {
/* even alpha=0 is ok, however, we cannot choose it since it kills the
k-space error formula.
@@ -965,8 +934,7 @@ int p3m_adaptive_tune(char **log) {
if (p3m.params.epsilon != P3M_EPSILON_METALLIC) {
if (!((box_geo.length()[0] == box_geo.length()[1]) &&
(box_geo.length()[1] == box_geo.length()[2]))) {
- *log = strcat_alloc(
- *log, "{049 P3M_init: Nonmetallic epsilon requires cubic box} ");
+ runtimeErrorMsg() << "non-metallic epsilon requires cubic box";
return ES_ERROR;
}
}
@@ -988,8 +956,7 @@ int p3m_adaptive_tune(char **log) {
*log = strcat_alloc(*log, b);
if (p3m.sum_qpart == 0) {
- *log = strcat_alloc(*log,
- "no charged particles in the system, cannot tune P3M");
+ runtimeErrorMsg() << "no charged particles in the system";
return ES_ERROR;
}
@@ -1002,13 +969,8 @@ int p3m_adaptive_tune(char **log) {
if (p3m.params.mesh[0] == 0 || p3m.params.mesh[1] == 0 ||
p3m.params.mesh[2] == 0) {
/* Medium-educated guess for the minimal mesh */
- mesh_density_min =
- pow(p3m.sum_qpart / (box_geo.length()[0] * box_geo.length()[1] *
- box_geo.length()[2]),
- 1.0 / 3.0);
- mesh_density_max = 512 / pow(box_geo.length()[0] * box_geo.length()[1] *
- box_geo.length()[2],
- 1.0 / 3.0);
+ mesh_density_min = pow(p3m.sum_qpart / box_geo.volume(), 1.0 / 3.0);
+ mesh_density_max = 512 / pow(box_geo.volume(), 1.0 / 3.0);
tune_mesh = true;
/* this limits the tried meshes if the accuracy cannot
be obtained with smaller meshes, but normally not all these
@@ -1114,6 +1076,8 @@ int p3m_adaptive_tune(char **log) {
p3m_m_time(log, tmp_mesh, cao_min, cao_max, &tmp_cao, r_cut_iL_min,
r_cut_iL_max, &tmp_r_cut_iL, &tmp_alpha_L, &tmp_accuracy);
/* some error occurred during the tuning force evaluation */
+ if (tmp_time == -P3M_TUNE_FAIL)
+ return ES_ERROR;
/* this mesh does not work at all */
if (tmp_time < 0.0)
continue;
@@ -1141,8 +1105,7 @@ int p3m_adaptive_tune(char **log) {
}
if (time_best == 1e20) {
- *log = strcat_alloc(*log,
- "failed to tune P3M parameters to required accuracy\n");
+ runtimeErrorMsg() << "failed to reach requested accuracy";
return ES_ERROR;
}
@@ -1194,7 +1157,7 @@ double p3m_real_space_error(double prefac, double r_cut_iL, int n_c_part,
double sum_q2, double alpha_L) {
return (2.0 * prefac * sum_q2 * exp(-Utils::sqr(r_cut_iL * alpha_L))) /
sqrt((double)n_c_part * r_cut_iL * box_geo.length()[0] *
- box_geo.length()[0] * box_geo.length()[1] * box_geo.length()[2]);
+ box_geo.volume());
}
double p3m_k_space_error(double prefac, const int mesh[3], int cao,
@@ -1316,7 +1279,7 @@ bool p3m_sanity_checks_system(const Utils::Vector3i &grid) {
if (p3m.params.epsilon != P3M_EPSILON_METALLIC) {
if (!((p3m.params.mesh[0] == p3m.params.mesh[1]) &&
(p3m.params.mesh[1] == p3m.params.mesh[2]))) {
- runtimeErrorMsg() << "P3M_init: Nonmetallic epsilon requires cubic box";
+ runtimeErrorMsg() << "P3M_init: non-metallic epsilon requires cubic box";
ret = true;
}
}
diff --git a/src/core/electrostatics_magnetostatics/p3m.hpp b/src/core/electrostatics_magnetostatics/p3m.hpp
index 71ee2fa0f3a..c006f5f0272 100644
--- a/src/core/electrostatics_magnetostatics/p3m.hpp
+++ b/src/core/electrostatics_magnetostatics/p3m.hpp
@@ -39,6 +39,7 @@
#include "fft.hpp"
#include "p3m-common.hpp"
+#include "p3m-data_struct.hpp"
#include "p3m_interpolation.hpp"
#include "p3m_send_mesh.hpp"
@@ -50,11 +51,9 @@
* data types
************************************************/
-struct p3m_data_struct {
+struct p3m_data_struct : public p3m_data_struct_base {
p3m_data_struct();
- P3MParameters params;
-
/** local mesh. */
p3m_local_mesh local_mesh;
/** real space mesh (local) for CA/FFT.*/
@@ -69,19 +68,8 @@ struct p3m_data_struct {
/** square of sum of charges (only on master node). */
double square_sum_q;
- /** Spatial differential operator in k-space. We use an i*k differentiation.
- */
- std::array, 3> d_op;
- /** Force optimised influence function (k-space) */
- std::vector g_force;
- /** Energy optimised influence function (k-space) */
- std::vector g_energy;
-
p3m_interpolation_cache inter_weights;
- /** number of permutations in k_space */
- int ks_pnum;
-
/** send/recv mesh sizes */
p3m_send_mesh sm;
diff --git a/src/core/electrostatics_magnetostatics/p3m_influence_function.hpp b/src/core/electrostatics_magnetostatics/p3m_influence_function.hpp
index 1bf0cb718c1..44e8d6436fc 100644
--- a/src/core/electrostatics_magnetostatics/p3m_influence_function.hpp
+++ b/src/core/electrostatics_magnetostatics/p3m_influence_function.hpp
@@ -23,6 +23,7 @@
#include
#include
+#include
#include
#include
#include
@@ -32,25 +33,6 @@
#include
namespace detail {
-enum : int { RX = 0, RY = 1, RZ = 2 };
-
-std::array, 3> inline calc_meshift(
- std::array const &mesh_size) {
- std::array, 3> ret;
-
- for (size_t i = 0; i < 3; i++) {
- ret[i].resize(mesh_size[i]);
-
- ret[i][0] = 0;
- for (int j = 1; j <= mesh_size[i] / 2; j++) {
- ret[i][j] = j;
- ret[i][mesh_size[i] - j] = -j;
- }
- }
-
- return ret;
-}
-
template T g_ewald(T alpha, T k2) {
auto constexpr limit = T{30};
auto const exponent = Utils::sqr(1. / (2. * alpha)) * k2;
@@ -62,11 +44,14 @@ template
std::pair aliasing_sums_ik(size_t cao, double alpha,
const Utils::Vector3d &k,
const Utils::Vector3d &h) {
+ using namespace detail::FFT_indexing;
using Utils::int_pow;
- using Utils::pi;
using Utils::sinc;
using Utils::Vector3d;
+ constexpr double two_pi = 2 * Utils::pi();
+ constexpr double two_pi_i = 1 / two_pi;
+
double numerator = 0.0;
double denominator = 0.0;
@@ -74,10 +59,10 @@ std::pair aliasing_sums_ik(size_t cao, double alpha,
for (int my = -m; my <= m; my++) {
for (int mz = -m; mz <= m; mz++) {
auto const km =
- k + 2 * pi() * Vector3d{mx / h[RX], my / h[RY], mz / h[RZ]};
- auto const U2 = std::pow(sinc(0.5 * km[RX] * h[RX] / pi()) *
- sinc(0.5 * km[RY] * h[RY] / pi()) *
- sinc(0.5 * km[RZ] * h[RZ] / pi()),
+ k + two_pi * Vector3d{mx / h[RX], my / h[RY], mz / h[RZ]};
+ auto const U2 = std::pow(sinc(km[RX] * h[RX] * two_pi_i) *
+ sinc(km[RY] * h[RY] * two_pi_i) *
+ sinc(km[RZ] * h[RZ] * two_pi_i),
2 * cao);
numerator += U2 * g_ewald(alpha, km.norm2()) * int_pow(k * km);
@@ -141,8 +126,7 @@ std::vector grid_influence_function(const P3MParameters ¶ms,
const Utils::Vector3i &n_start,
const Utils::Vector3i &n_end,
const Utils::Vector3d &box_l) {
- enum : int { RX = 0, RY = 1, RZ = 2 };
- enum : int { KY = 0, KZ = 1, KX = 2 };
+ using namespace detail::FFT_indexing;
auto const shifts =
detail::calc_meshift({params.mesh[0], params.mesh[1], params.mesh[2]});
diff --git a/src/core/electrostatics_magnetostatics/specfunc.cpp b/src/core/electrostatics_magnetostatics/specfunc.cpp
index 29590794f56..a0aabad72d4 100644
--- a/src/core/electrostatics_magnetostatics/specfunc.cpp
+++ b/src/core/electrostatics_magnetostatics/specfunc.cpp
@@ -47,6 +47,8 @@
#include "specfunc.hpp"
#include
+#include
+
/************************************************
* chebychev expansions
************************************************/
@@ -243,7 +245,7 @@ double K0(double x) {
if (x <= 2.0) {
c = evaluateAsChebychevSeriesAt(bk0_cs, 0.5 * x * x - 1.0);
I0 = evaluateAsChebychevSeriesAt(bi0_cs, x * x / 4.5 - 1.0);
- return (-log(x) + M_LN2) * I0 + c;
+ return (-log(x) + Utils::ln_2()) * I0 + c;
}
c = (x <= 8.0) ? evaluateAsChebychevSeriesAt(ak0_cs, (16.0 / x - 5.0) / 3.0)
: evaluateAsChebychevSeriesAt(ak02_cs, 16.0 / x - 1.0);
@@ -255,7 +257,7 @@ double K1(double x) {
if (x <= 2.0) {
c = evaluateAsChebychevSeriesAt(bk1_cs, 0.5 * x * x - 1.0);
I1 = x * evaluateAsChebychevSeriesAt(bi1_cs, x * x / 4.5 - 1.0);
- return (log(x) - M_LN2) * I1 + c / x;
+ return (log(x) - Utils::ln_2()) * I1 + c / x;
}
c = (x <= 8.0) ? evaluateAsChebychevSeriesAt(ak1_cs, (16.0 / x - 5.0) / 3.0)
: evaluateAsChebychevSeriesAt(ak12_cs, 16.0 / x - 1.0);
@@ -320,7 +322,7 @@ double LPK0(double x) {
d0 = x2 * d0 - dd0 + bi0_cs[j];
dd0 = tmp0;
}
- auto const tmp = log(x) - M_LN2;
+ auto const tmp = log(x) - Utils::ln_2();
ret = -tmp * (0.5 * (bi0_cs[0] + x2 * d0) - dd0);
/* K0/K1 correction */
@@ -381,7 +383,7 @@ double LPK1(double x) {
d1 = x2 * d1 - dd1 + bi1_cs[j];
dd1 = tmp1;
}
- auto const tmp = log(x) - M_LN2;
+ auto const tmp = log(x) - Utils::ln_2();
ret = x * tmp * (0.5 * (bi1_cs[0] + x2 * d1) - dd1);
/* K0/K1 correction */
@@ -458,7 +460,7 @@ std::tuple LPK01(double x) {
dd0 = tmp0;
dd1 = tmp1;
}
- auto const tmp = log(x) - M_LN2;
+ auto const tmp = log(x) - Utils::ln_2();
auto K0 = -tmp * (0.5 * (bi0_cs[0] + x2 * d0) - dd0);
auto K1 = x * tmp * (0.5 * (bi1_cs[0] + x2 * d1) - dd1);
diff --git a/src/core/io/writer/h5md_specification.cpp b/src/core/io/writer/h5md_specification.cpp
index 853c570de08..2e40000fd17 100644
--- a/src/core/io/writer/h5md_specification.cpp
+++ b/src/core/io/writer/h5md_specification.cpp
@@ -36,6 +36,5 @@ std::array H5MD_Specification::DATASETS = {{
{"connectivity/atoms", "step", 1, H5T_NATIVE_INT, 1, true},
{"connectivity/atoms", "time", 1, H5T_NATIVE_DOUBLE, 1, true},
}};
-
}
} // namespace Writer
diff --git a/src/core/nonbonded_interactions/ljcos2.hpp b/src/core/nonbonded_interactions/ljcos2.hpp
index f19e62f0fcb..1950c66b9ef 100644
--- a/src/core/nonbonded_interactions/ljcos2.hpp
+++ b/src/core/nonbonded_interactions/ljcos2.hpp
@@ -37,6 +37,8 @@
#include "nonbonded_interaction_data.hpp"
#include
+
+#include
#include
int ljcos2_set_params(int part_type_a, int part_type_b, double eps, double sig,
@@ -53,8 +55,10 @@ inline double ljcos2_pair_force_factor(IA_parameters const &ia_params,
fac =
48.0 * ia_params.ljcos2.eps * frac6 * (frac6 - 0.5) / (r_off * dist);
} else if (r_off < ia_params.ljcos2.rchange + ia_params.ljcos2.w) {
- fac = -ia_params.ljcos2.eps * M_PI / 2 / ia_params.ljcos2.w / dist *
- sin(M_PI * (r_off - ia_params.ljcos2.rchange) / ia_params.ljcos2.w);
+ fac = -ia_params.ljcos2.eps * Utils::pi() / 2 / ia_params.ljcos2.w /
+ dist *
+ sin(Utils::pi() * (r_off - ia_params.ljcos2.rchange) /
+ ia_params.ljcos2.w);
}
return fac;
}
@@ -78,7 +82,7 @@ inline double ljcos2_pair_energy(IA_parameters const &ia_params, double dist) {
}
if (r_off < (ia_params.ljcos2.rchange + ia_params.ljcos2.w)) {
return -ia_params.ljcos2.eps / 2 *
- (cos(M_PI * (r_off - ia_params.ljcos2.rchange) /
+ (cos(Utils::pi() * (r_off - ia_params.ljcos2.rchange) /
ia_params.ljcos2.w) +
1);
}
diff --git a/src/core/observables/EnergyObservable.hpp b/src/core/observables/EnergyObservable.hpp
index e366784b877..2b8bb6cde07 100644
--- a/src/core/observables/EnergyObservable.hpp
+++ b/src/core/observables/EnergyObservable.hpp
@@ -17,6 +17,7 @@
* along with this program. If not, see .
*/
#ifndef OBSERVABLES_ENERGY_HPP
+#define OBSERVABLES_ENERGY_HPP
#include "Observable.hpp"
#include "energy.hpp"
diff --git a/src/core/observables/PressureObservable.hpp b/src/core/observables/PressureObservable.hpp
index 171876532ce..a3d755eeeec 100644
--- a/src/core/observables/PressureObservable.hpp
+++ b/src/core/observables/PressureObservable.hpp
@@ -17,6 +17,7 @@
* along with this program. If not, see .
*/
#ifndef OBSERVABLES_PRESSURE_HPP
+#define OBSERVABLES_PRESSURE_HPP
#include "Observable.hpp"
#include "pressure.hpp"
diff --git a/src/core/observables/PressureTensor.hpp b/src/core/observables/PressureTensor.hpp
index 35684e90628..7ded56dce87 100644
--- a/src/core/observables/PressureTensor.hpp
+++ b/src/core/observables/PressureTensor.hpp
@@ -17,6 +17,7 @@
* along with this program. If not, see .
*/
#ifndef OBSERVABLES_PRESSURETENSOR_HPP
+#define OBSERVABLES_PRESSURETENSOR_HPP
#include "Observable.hpp"
#include "pressure.hpp"
diff --git a/src/core/thermostat.hpp b/src/core/thermostat.hpp
index f87b600ede1..886c03be464 100644
--- a/src/core/thermostat.hpp
+++ b/src/core/thermostat.hpp
@@ -348,6 +348,11 @@ NEW_THERMOSTAT(thermalized_bond)
NEW_THERMOSTAT(dpd)
#endif
+/* Exported thermostat globals */
+extern LangevinThermostat langevin;
+extern BrownianThermostat brownian;
+extern IsotropicNptThermostat npt_iso;
+
/** Initialize constants of the thermostat at the start of integration */
void thermo_init();
diff --git a/src/core/unit_tests/CMakeLists.txt b/src/core/unit_tests/CMakeLists.txt
index abcc0a53e8b..e5a8bb4e837 100644
--- a/src/core/unit_tests/CMakeLists.txt
+++ b/src/core/unit_tests/CMakeLists.txt
@@ -20,30 +20,19 @@
include(unit_test)
# Add tests here
-
unit_test(NAME RuntimeError_test SRC RuntimeError_test.cpp DEPENDS
Boost::serialization)
unit_test(NAME RuntimeErrorCollector_test SRC RuntimeErrorCollector_test.cpp
DEPENDS EspressoCore Boost::mpi MPI::MPI_CXX)
-unit_test(NAME ScriptInterface_test SRC ScriptInterface_test.cpp DEPENDS
- ScriptInterface)
unit_test(NAME MpiCallbacks_test SRC MpiCallbacks_test.cpp DEPENDS
EspressoUtils Boost::mpi MPI::MPI_CXX NUM_PROC 2)
-unit_test(NAME ParallelScriptInterface_test SRC
- ParallelScriptInterface_test.cpp DEPENDS ScriptInterface Boost::mpi
- MPI::MPI_CXX NUM_PROC 2)
-unit_test(NAME AutoParameters_test SRC AutoParameters_test.cpp DEPENDS
- ScriptInterface)
-unit_test(NAME AutoParameter_test SRC AutoParameter_test.cpp DEPENDS
- ScriptInterface)
-unit_test(NAME Variant_test SRC Variant_test.cpp DEPENDS ScriptInterface)
unit_test(NAME ParticleIterator_test SRC ParticleIterator_test.cpp DEPENDS
EspressoUtils)
+unit_test(NAME p3m_test SRC p3m_test.cpp DEPENDS EspressoUtils)
unit_test(NAME link_cell_test SRC link_cell_test.cpp DEPENDS EspressoUtils)
unit_test(NAME Particle_test SRC Particle_test.cpp DEPENDS EspressoUtils
Boost::serialization)
-unit_test(NAME get_value SRC get_value_test.cpp DEPENDS ScriptInterface)
unit_test(NAME field_coupling_couplings SRC field_coupling_couplings_test.cpp
DEPENDS EspressoUtils)
unit_test(NAME field_coupling_fields SRC field_coupling_fields_test.cpp DEPENDS
@@ -51,7 +40,6 @@ unit_test(NAME field_coupling_fields SRC field_coupling_fields_test.cpp DEPENDS
unit_test(NAME field_coupling_force_field SRC
field_coupling_force_field_test.cpp DEPENDS EspressoUtils)
unit_test(NAME periodic_fold_test SRC periodic_fold_test.cpp)
-unit_test(NAME None_test SRC None_test.cpp DEPENDS ScriptInterface)
unit_test(NAME grid_test SRC grid_test.cpp DEPENDS EspressoCore)
unit_test(NAME BoxGeometry_test SRC BoxGeometry_test.cpp DEPENDS EspressoCore)
unit_test(NAME LocalBox_test SRC LocalBox_test.cpp DEPENDS EspressoCore)
diff --git a/src/core/unit_tests/ParallelScriptInterface_test.cpp b/src/core/unit_tests/ParallelScriptInterface_test.cpp
deleted file mode 100644
index cb7aa25acfa..00000000000
--- a/src/core/unit_tests/ParallelScriptInterface_test.cpp
+++ /dev/null
@@ -1,185 +0,0 @@
-/*
- * Copyright (C) 2010-2019 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 .
- */
-#include
-
-#define BOOST_TEST_NO_MAIN
-#define BOOST_TEST_MODULE ParallelScriptInterface test
-#define BOOST_TEST_ALTERNATIVE_INIT_API
-#define BOOST_TEST_DYN_LINK
-#include
-#include
-
-#include "MpiCallbacks.hpp"
-
-#include "script_interface/ParallelScriptInterface.hpp"
-
-namespace mpi = boost::mpi;
-std::unique_ptr callbacks;
-
-using namespace ScriptInterface;
-
-struct TestClass : public ScriptInterfaceBase {
- TestClass() { constructed = true; }
- ~TestClass() override { destructed = true; }
-
- void set_parameter(const std::string &name, const Variant &value) override {
- last_parameter = make_pair(name, value);
-
- if (name == "obj_param") {
- obj_param = get_instance(boost::get(value)).lock();
- }
- }
-
- Variant get_parameter(std::string const &name) const override {
- if (name == "obj_param") {
- return obj_param->id();
- }
- return last_parameter.second;
- }
-
- Variant call_method(const std::string &method,
- const VariantMap ¶ms) override {
- last_method_parameters = make_pair(method, params);
-
- return std::string("TestResult");
- }
-
- static std::pair last_method_parameters;
- static std::pair last_parameter;
-
- std::shared_ptr obj_param;
-
- static bool constructed;
- static bool destructed;
-};
-
-bool TestClass::constructed = false;
-bool TestClass::destructed = false;
-std::pair TestClass::last_method_parameters;
-std::pair TestClass::last_parameter;
-
-/*
- * Check that instances are created and correctly destroyed on
- * the slave nodes.
- */
-BOOST_AUTO_TEST_CASE(ctor_dtor) {
- /* Reset */
- TestClass::constructed = false;
- TestClass::destructed = false;
-
- if (callbacks->comm().rank() == 0) {
- /* Create an instance everywhere */
- auto so = std::make_shared("TestClass");
- /* Force destruction */
- so = nullptr;
-
- callbacks->abort_loop();
- } else {
- callbacks->loop();
- }
-
- /* Check that ctor and dtor were run on all nodes */
- BOOST_CHECK(TestClass::constructed);
- BOOST_CHECK(TestClass::destructed);
-}
-
-/*
- * Check that parameters are forwarded correctly.
- */
-BOOST_AUTO_TEST_CASE(set_parameter) {
- if (callbacks->comm().rank() == 0) {
- auto so = std::make_shared("TestClass");
-
- so->set_parameter("TestParam", std::string("TestValue"));
-
- callbacks->abort_loop();
- } else {
- callbacks->loop();
- }
-
- auto const &last_parameter = TestClass::last_parameter;
- BOOST_CHECK(last_parameter.first == "TestParam");
- BOOST_CHECK(boost::get(last_parameter.second) == "TestValue");
-}
-
-/*
- * Check that the method name and parameters are forwarded correctly
- * to the payload object, and check that the return value is
- * propagated correctly.
- */
-BOOST_AUTO_TEST_CASE(call_method) {
- const VariantMap params{{"TestParam", std::string("TestValue")}};
- const std::string method{"TestMethod"};
-
- if (callbacks->comm().rank() == 0) {
- auto so = std::make_shared("TestClass");
-
- auto result = so->call_method(method, params);
-
- /* Check return value */
- BOOST_CHECK(boost::get(result) == "TestResult");
-
- callbacks->abort_loop();
- } else {
- callbacks->loop();
- }
-
- auto const &last_parameters = TestClass::last_method_parameters;
- BOOST_CHECK(last_parameters.first == method);
- BOOST_CHECK(last_parameters.second == params);
-}
-
-BOOST_AUTO_TEST_CASE(parameter_lifetime) {
- if (callbacks->comm().rank() == 0) {
- auto host = std::make_shared("TestClass");
- ScriptInterfaceBase *bare_ptr;
-
- {
- auto parameter = ScriptInterfaceBase::make_shared(
- "TestClass", ScriptInterfaceBase::CreationPolicy::GLOBAL);
- bare_ptr = parameter.get();
-
- BOOST_CHECK(get_instance(parameter->id()) == parameter);
-
- host->set_parameter("obj_param", parameter->id());
- }
-
- auto param_id = host->get_parameter("obj_param");
- auto parameter = get_instance(param_id);
-
- /* Check that we got the original instance back */
- BOOST_CHECK(parameter.get() == bare_ptr);
-
- callbacks->abort_loop();
- } else {
- callbacks->loop();
- }
-}
-
-int main(int argc, char **argv) {
- mpi::environment mpi_env(argc, argv);
- mpi::communicator world;
- callbacks = std::make_unique(
- world, /* abort_on_exit */ false);
-
- ParallelScriptInterface::initialize(*callbacks);
- register_new("TestClass");
-
- return boost::unit_test::unit_test_main(init_unit_test, argc, argv);
-}
diff --git a/src/core/unit_tests/ScriptInterface_test.cpp b/src/core/unit_tests/ScriptInterface_test.cpp
deleted file mode 100644
index b4e1e560943..00000000000
--- a/src/core/unit_tests/ScriptInterface_test.cpp
+++ /dev/null
@@ -1,111 +0,0 @@
-/*
- * Copyright (C) 2010-2019 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 .
- */
-
-#include
-
-#define BOOST_TEST_MODULE ScriptInterface test
-#define BOOST_TEST_DYN_LINK
-#include
-
-#include "script_interface/ScriptInterface.hpp"
-
-using std::map;
-using std::string;
-using std::vector;
-
-using namespace ScriptInterface;
-
-namespace Testing {
-
-/* Mock to test ScriptInterface. */
-struct ScriptInterfaceTest : public ScriptInterface::ScriptInterfaceBase {
- VariantMap get_parameters() const override {
- VariantMap ret;
-
- ret["bool_opt"] = bool_opt;
- ret["integer"] = integer;
- ret["bool_req"] = bool_req;
- ret["vec_double"] = vec_double;
- ret["vec_int"] = vec_int;
-
- return ret;
- }
-
- /* Not needed for testing */
- Utils::Span valid_parameters() const override {
- return {};
- }
-
- void set_parameter(const string &name, const Variant &value) override {
- if (name == "bool_opt") {
- bool_opt =
- get_value::type>(value);
- };
- if (name == "integer") {
- integer =
- get_value::type>(value);
- };
- if (name == "bool_req") {
- bool_req =
- get_value::type>(value);
- };
- if (name == "vec_double") {
- vec_double =
- get_value::type>(value);
- };
- if (name == "vec_int") {
- vec_int =
- get_value::type>(value);
- };
- }
-
- Variant call_method(const std::string &name,
- const VariantMap ¶ms) override {
- if (name == "test_method") {
- }
-
- return true;
- }
-
- bool bool_opt, bool_req;
- int integer;
- vector vec_double;
- vector vec_int;
-};
-} /* namespace Testing */
-
-using namespace Testing;
-
-BOOST_AUTO_TEST_CASE(non_copyable) {
- static_assert(!std::is_copy_constructible::value, "");
- static_assert(!std::is_copy_assignable::value, "");
-}
-
-/*
- * We check the default implementations of set_parameters
- * and get_parameter of ScriptInterface (this is the only
- * logic in the class).
- */
-BOOST_AUTO_TEST_CASE(default_implementation) {
- ScriptInterfaceTest si_test;
-
- si_test.call_method("test_method", {});
-}
diff --git a/src/core/unit_tests/p3m_test.cpp b/src/core/unit_tests/p3m_test.cpp
new file mode 100644
index 00000000000..5d1687685e9
--- /dev/null
+++ b/src/core/unit_tests/p3m_test.cpp
@@ -0,0 +1,52 @@
+/*
+ * Copyright (C) 2010-2020 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 p3m test
+#define BOOST_TEST_DYN_LINK
+#include
+
+#include "electrostatics_magnetostatics/p3m-common.hpp"
+
+BOOST_AUTO_TEST_CASE(calc_meshift_false) {
+ std::array, 3> const ref = {
+ std::vector{0}, std::vector{0, 1, -2, -1},
+ std::vector{0, 1, 2, 3, -3, -2, -1}};
+
+ auto const val = detail::calc_meshift({1, 4, 7}, false);
+
+ for (size_t i = 0; i < 3; ++i) {
+ for (size_t j = 0; j < ref[i].size(); ++j) {
+ BOOST_CHECK_EQUAL(val[i][j], ref[i][j]);
+ }
+ }
+}
+
+BOOST_AUTO_TEST_CASE(calc_meshift_true) {
+ std::array, 3> const ref = {
+ std::vector{0}, std::vector{0, 1, 0, -1},
+ std::vector{0, 1, 2, 0, -3, -2, -1}};
+
+ auto const val = detail::calc_meshift({1, 4, 7}, true);
+
+ for (size_t i = 0; i < 3; ++i) {
+ for (size_t j = 0; j < ref[i].size(); ++j) {
+ BOOST_CHECK_EQUAL(val[i][j], ref[i][j]);
+ }
+ }
+}
diff --git a/src/python/espressomd/_init.pyx b/src/python/espressomd/_init.pyx
index 84e7f27546e..af2673c5be1 100644
--- a/src/python/espressomd/_init.pyx
+++ b/src/python/espressomd/_init.pyx
@@ -17,7 +17,7 @@
# along with this program. If not, see .
#
import sys
-from . import script_interface
+from . cimport script_interface
from . cimport communication
from libcpp.memory cimport shared_ptr
from boost cimport environment
@@ -28,7 +28,7 @@ communication.init(mpi_env)
# Initialize script interface
# Has to be _after_ mpi_init
-script_interface.init()
+script_interface.init(communication.mpiCallbacks())
# Block the slaves in the callback loop
# The master is just returning to the user script
diff --git a/src/python/espressomd/actors.pyx b/src/python/espressomd/actors.pyx
index 53747d92948..30d6d53998d 100644
--- a/src/python/espressomd/actors.pyx
+++ b/src/python/espressomd/actors.pyx
@@ -57,7 +57,7 @@ cdef class Actor:
if k in self.valid_keys():
self._params[k] = kwargs[k]
else:
- raise KeyError("%s is not a valid key" % k)
+ raise KeyError(f"{k} is not a valid key")
def _activate(self):
inter = self._get_interaction_type()
diff --git a/src/python/espressomd/analyze.pyx b/src/python/espressomd/analyze.pyx
index d1f80e9c9ce..163e3fff588 100644
--- a/src/python/espressomd/analyze.pyx
+++ b/src/python/espressomd/analyze.pyx
@@ -186,12 +186,12 @@ class Analysis:
for i in range(len(p1)):
if not is_valid_type(p1[i], int):
raise TypeError(
- "Particle types in p1 have to be of type int, got: " + repr(p1[i]))
+ f"Particle types in p1 have to be of type int, got: {repr(p1[i])}")
for i in range(len(p2)):
if not is_valid_type(p2[i], int):
raise TypeError(
- "Particle types in p2 have to be of type int, got: " + repr(p2[i]))
+ f"Particle types in p2 have to be of type int, got: {repr(p2[i])}")
return analyze.mindist(analyze.partCfg(), p1, p2)
@@ -249,7 +249,7 @@ class Analysis:
"The p_type keyword argument must be provided (particle type)")
check_type_or_throw_except(p_type, 1, int, "p_type has to be an int")
if p_type < 0 or p_type >= analyze.max_seen_particle_type:
- raise ValueError("Particle type {} does not exist!".format(p_type))
+ raise ValueError(f"Particle type {p_type} does not exist!")
return analyze.centerofmass(analyze.partCfg(), p_type)
@@ -736,8 +736,7 @@ class Analysis:
check_type_or_throw_except(
ptype, 1, int, "particle type has to be an int")
if ptype < 0 or ptype >= analyze.max_seen_particle_type:
- raise ValueError(
- "Particle type {} does not exist!".format(ptype))
+ raise ValueError(f"Particle type {ptype} does not exist!")
selection = self._system.part.select(lambda p: (p.type in p_type))
cm = np.mean(selection.pos, axis=0)
mat = np.zeros(shape=(3, 3))
@@ -787,7 +786,7 @@ class Analysis:
"The p_type keyword argument must be provided (particle type)")
check_type_or_throw_except(p_type, 1, int, "p_type has to be an int")
if p_type < 0 or p_type >= analyze.max_seen_particle_type:
- raise ValueError("Particle type {} does not exist!".format(p_type))
+ raise ValueError(f"Particle type {p_type} does not exist!")
analyze.momentofinertiamatrix(
analyze.partCfg(), p_type, MofImatrix)
diff --git a/src/python/espressomd/cellsystem.pyx b/src/python/espressomd/cellsystem.pyx
index 02cfc4ea458..48eadba1185 100644
--- a/src/python/espressomd/cellsystem.pyx
+++ b/src/python/espressomd/cellsystem.pyx
@@ -137,8 +137,8 @@ cdef class CellSystem:
def __set__(self, _node_grid):
if not np.prod(_node_grid) == n_nodes:
- raise ValueError("Number of available nodes " + str(
- n_nodes) + " and imposed node grid " + str(_node_grid) + " do not agree.")
+ raise ValueError(
+ f"Number of available nodes {n_nodes} and imposed node grid {_node_grid} do not agree.")
else:
node_grid[0] = _node_grid[0]
node_grid[1] = _node_grid[1]
diff --git a/src/python/espressomd/checkpointing.py b/src/python/espressomd/checkpointing.py
index 57141e10896..37023468dae 100644
--- a/src/python/espressomd/checkpointing.py
+++ b/src/python/espressomd/checkpointing.py
@@ -67,7 +67,7 @@ def __init__(self, checkpoint_id=None, checkpoint_path="."):
# update checkpoint counter
self.counter = 0
while os.path.isfile(os.path.join(
- self.checkpoint_dir, "{}.checkpoint".format(self.counter))):
+ self.checkpoint_dir, f"{self.counter}.checkpoint")):
self.counter += 1
# init signals
@@ -134,11 +134,11 @@ def register(self, *args):
# if not a in dir(self.calling_module):
if not self.__hasattr_submodule(self.calling_module, a):
raise KeyError(
- "The given object '{}' was not found in the current scope.".format(a))
+ f"The given object '{a}' was not found in the current scope.")
if a in self.checkpoint_objects:
raise KeyError(
- "The given object '{}' is already registered for checkpointing.".format(a))
+ f"The given object '{a}' is already registered for checkpointing.")
self.checkpoint_objects.append(a)
@@ -154,7 +154,7 @@ def unregister(self, *args):
for a in args:
if not isinstance(a, str) or a not in self.checkpoint_objects:
raise KeyError(
- "The given object '{}' was not registered for checkpointing yet.".format(a))
+ f"The given object '{a}' was not registered for checkpointing yet.")
self.checkpoint_objects.remove(a)
@@ -204,7 +204,7 @@ def save(self, checkpoint_index=None):
if checkpoint_index is None:
checkpoint_index = self.counter
filename = os.path.join(
- self.checkpoint_dir, "{}.checkpoint".format(checkpoint_index))
+ self.checkpoint_dir, f"{checkpoint_index}.checkpoint")
tmpname = filename + ".__tmp__"
with open(tmpname, "wb") as checkpoint_file:
@@ -226,7 +226,7 @@ def load(self, checkpoint_index=None):
checkpoint_index = self.get_last_checkpoint_index()
filename = os.path.join(
- self.checkpoint_dir, "{}.checkpoint".format(checkpoint_index))
+ self.checkpoint_dir, f"{checkpoint_index}.checkpoint")
with open(filename, "rb") as f:
checkpoint_data = pickle.load(f)
@@ -288,7 +288,7 @@ def register_signal(self, signum=None):
if signum in self.checkpoint_signals:
raise KeyError(
- "The signal {} is already registered for checkpointing.".format(signum))
+ f"The signal {signum} is already registered for checkpointing.")
signal.signal(signum, self.__signal_handler)
self.checkpoint_signals.append(signum)
diff --git a/src/python/espressomd/collision_detection.pyx b/src/python/espressomd/collision_detection.pyx
index 9bd28373acf..4674f84b789 100644
--- a/src/python/espressomd/collision_detection.pyx
+++ b/src/python/espressomd/collision_detection.pyx
@@ -216,7 +216,7 @@ class CollisionDetection(ScriptInterfaceHelper):
for key in self._int_mode:
if self._int_mode[key] == int_mode:
return key
- raise Exception("Unknown integer collision mode %d" % int_mode)
+ raise Exception(f"Unknown integer collision mode {int_mode}")
# Pickle support
def __reduce__(self):
diff --git a/src/python/espressomd/constraints.py b/src/python/espressomd/constraints.py
index 164c71310ac..d5cdbe5f400 100644
--- a/src/python/espressomd/constraints.py
+++ b/src/python/espressomd/constraints.py
@@ -238,7 +238,7 @@ class that can calculate the required grid dimensions and the coordinates.
"""
def __init__(self, **kwargs):
- if "oid" not in kwargs:
+ if "sip" not in kwargs:
field = kwargs.pop("field")
shape, codim = self._unpack_dims(field)
super().__init__(_field_shape=shape, _field_codim=codim,
@@ -403,7 +403,7 @@ class Gravity(Constraint):
"""
def __init__(self, **kwargs):
- if "oid" not in kwargs:
+ if "sip" not in kwargs:
kwargs["value"] = kwargs.pop("g")
super().__init__(**kwargs)
@@ -441,7 +441,7 @@ class LinearElectricPotential(Constraint):
"""
def __init__(self, phi0=0, **kwargs):
- if "oid" not in kwargs:
+ if "sip" not in kwargs:
kwargs["A"] = -np.array(kwargs.pop("E"))
kwargs["b"] = phi0
super().__init__(**kwargs)
@@ -489,7 +489,7 @@ class ElectricPlaneWave(Constraint):
_so_name = "Constraints::ElectricPlaneWave"
def __init__(self, phi=0, **kwargs):
- if "oid" not in kwargs:
+ if "sip" not in kwargs:
kwargs["amplitude"] = kwargs.pop("E0")
kwargs["wave_vector"] = kwargs.pop("k")
kwargs["frequency"] = kwargs.pop("omega")
@@ -561,7 +561,7 @@ class HomogeneousFlowField(Constraint):
"""
def __init__(self, **kwargs):
- if "oid" not in kwargs:
+ if "sip" not in kwargs:
kwargs["value"] = kwargs.pop("u")
super().__init__(**kwargs)
diff --git a/src/python/espressomd/electrokinetics.pyx b/src/python/espressomd/electrokinetics.pyx
index ad24bf1242f..162016df24a 100644
--- a/src/python/espressomd/electrokinetics.pyx
+++ b/src/python/espressomd/electrokinetics.pyx
@@ -43,7 +43,7 @@ IF ELECTROKINETICS:
return ElectrokineticsRoutines(np.array(key))
else:
raise Exception(
- "%s is not a valid key. Should be a point on the nodegrid e.g. ek[0,0,0]," % key)
+ f"{key} is not a valid key. Should be a point on the nodegrid e.g. ek[0,0,0].")
def validate_params(self):
"""
@@ -431,7 +431,7 @@ IF ELECTROKINETICS:
return SpecieRoutines(np.array(key), self.id)
else:
raise Exception(
- "%s is not a valid key. Should be a point on the nodegrid e.g. species[0,0,0]," % key)
+ f"{key} is not a valid key. Should be a point on the nodegrid e.g. species[0,0,0].")
def __init__(self, **kwargs):
Species.py_number_of_species += 1
@@ -449,7 +449,7 @@ IF ELECTROKINETICS:
if k in self.valid_keys():
self._params[k] = kwargs[k]
else:
- raise KeyError("%s is not a valid key" % k)
+ raise KeyError(f"{k} is not a valid key")
def valid_keys(self):
"""
diff --git a/src/python/espressomd/electrostatics.pyx b/src/python/espressomd/electrostatics.pyx
index 92ea05ba72f..250f4920420 100644
--- a/src/python/espressomd/electrostatics.pyx
+++ b/src/python/espressomd/electrostatics.pyx
@@ -261,11 +261,11 @@ IF P3M == 1:
raise ValueError("P3M accuracy has to be positive")
if self._params["epsilon"] == "metallic":
- self._params = 0.0
+ self._params["epsilon"] = 0.0
- if not (is_valid_type(self._params["epsilon"], float)
- or self._params["epsilon"] == "metallic"):
- raise ValueError("epsilon should be a double or 'metallic'")
+ check_type_or_throw_except(
+ self._params["epsilon"], 1, float,
+ "epsilon should be a double or 'metallic'")
if self._params["mesh_off"] != default_params["mesh_off"]:
check_type_or_throw_except(self._params["mesh_off"], 3, float,
@@ -273,8 +273,7 @@ IF P3M == 1:
if not (self._params["alpha"] == default_params["alpha"]
or self._params["alpha"] > 0):
- raise ValueError(
- "alpha should be positive")
+ raise ValueError("alpha should be positive")
def valid_keys(self):
return ["mesh", "cao", "accuracy", "epsilon", "alpha", "r_cut",
@@ -317,15 +316,14 @@ IF P3M == 1:
def _tune(self):
set_prefactor(self._params["prefactor"])
+ p3m_set_eps(self._params["epsilon"])
python_p3m_set_tune_params(self._params["r_cut"],
self._params["mesh"],
self._params["cao"],
-1.0,
self._params["accuracy"])
resp = python_p3m_adaptive_tune()
- if resp:
- raise Exception(
- "failed to tune P3M parameters to required accuracy")
+ handle_errors("P3M tuning failed")
self._params.update(self._get_params_from_es_core())
def _activate_method(self):
@@ -402,13 +400,12 @@ IF P3M == 1:
if not (self._params["accuracy"] >= 0):
raise ValueError("P3M accuracy has to be positive")
- # if self._params["epsilon"] == "metallic":
- # self._params = 0.0
+ if self._params["epsilon"] == "metallic":
+ self._params["epsilon"] = 0.0
- if not (is_valid_type(self._params["epsilon"], float)
- or self._params["epsilon"] == "metallic"):
- raise ValueError(
- "epsilon should be a double or 'metallic'")
+ check_type_or_throw_except(
+ self._params["epsilon"], 1, float,
+ "epsilon should be a double or 'metallic'")
if self._params["mesh_off"] != default_params["mesh_off"]:
check_type_or_throw_except(self._params["mesh_off"], 3, float,
@@ -441,15 +438,14 @@ IF P3M == 1:
def _tune(self):
set_prefactor(self._params["prefactor"])
+ p3m_set_eps(self._params["epsilon"])
python_p3m_set_tune_params(self._params["r_cut"],
self._params["mesh"],
self._params["cao"],
-1.0,
self._params["accuracy"])
resp = python_p3m_adaptive_tune()
- if resp:
- raise Exception(
- "failed to tune P3M parameters to required accuracy")
+ handle_errors("P3MGPU tuning failed")
self._params.update(self._get_params_from_es_core())
def _activate_method(self):
diff --git a/src/python/espressomd/highlander.py b/src/python/espressomd/highlander.py
index f27d110062c..d170442bded 100644
--- a/src/python/espressomd/highlander.py
+++ b/src/python/espressomd/highlander.py
@@ -25,8 +25,7 @@ def __init__(self, cls):
self._cls = cls
def __str__(self):
- return "There can only be one instance of '{}' at any time.".format(
- self._cls)
+ return f"There can only be one instance of '{self._cls}' at any time."
def highlander(klass):
diff --git a/src/python/espressomd/interactions.pyx b/src/python/espressomd/interactions.pyx
index 791cd43fac6..ce67701e173 100644
--- a/src/python/espressomd/interactions.pyx
+++ b/src/python/espressomd/interactions.pyx
@@ -1825,37 +1825,37 @@ class BondedInteractionNotDefined:
self.__class__.__name__ + " not compiled into ESPResSo core")
def type_number(self):
- raise Exception(("%s has to be defined in myconfig.hpp.") % self.name)
+ raise Exception(f"{self.name} has to be defined in myconfig.hpp.")
def type_name(self):
"""Name of interaction type.
"""
- raise Exception(("%s has to be defined in myconfig.hpp.") % self.name)
+ raise Exception(f"{self.name} has to be defined in myconfig.hpp.")
def valid_keys(self):
"""All parameters that can be set.
"""
- raise Exception(("%s has to be defined in myconfig.hpp.") % self.name)
+ raise Exception(f"{self.name} has to be defined in myconfig.hpp.")
def required_keys(self):
"""Parameters that have to be set.
"""
- raise Exception(("%s has to be defined in myconfig.hpp.") % self.name)
+ raise Exception(f"{self.name} has to be defined in myconfig.hpp.")
def set_default_params(self):
"""Sets parameters that are not required to their default value.
"""
- raise Exception(("%s has to be defined in myconfig.hpp.") % self.name)
+ raise Exception(f"{self.name} has to be defined in myconfig.hpp.")
def _get_params_from_es_core(self):
- raise Exception(("%s has to be defined in myconfig.hpp.") % self.name)
+ raise Exception(f"{self.name} has to be defined in myconfig.hpp.")
def _set_params_in_es_core(self):
- raise Exception(("%s has to be defined in myconfig.hpp.") % self.name)
+ raise Exception(f"{self.name} has to be defined in myconfig.hpp.")
class FeneBond(BondedInteraction):
@@ -2591,8 +2591,8 @@ class TabulatedAngle(_TabulatedBase):
"""
phi = [self._params["min"], self._params["max"]]
if abs(phi[0] - 0.) > 1e-5 or abs(phi[1] - self.pi) > 1e-5:
- raise ValueError("Tabulated angle expects forces/energies "
- "within the range [0, pi], got " + str(phi))
+ raise ValueError(f"Tabulated angle expects forces/energies "
+ f"within the range [0, pi], got {phi}")
class TabulatedDihedral(_TabulatedBase):
@@ -2631,8 +2631,8 @@ class TabulatedDihedral(_TabulatedBase):
"""
phi = [self._params["min"], self._params["max"]]
if abs(phi[0] - 0.) > 1e-5 or abs(phi[1] - 2 * self.pi) > 1e-5:
- raise ValueError("Tabulated dihedral expects forces/energies "
- "within the range [0, 2*pi], got " + str(phi))
+ raise ValueError(f"Tabulated dihedral expects forces/energies "
+ f"within the range [0, 2*pi], got {phi}")
IF TABULATED == 1:
@@ -3313,7 +3313,7 @@ class BondedInteractions:
# Check if the bonded interaction exists in ESPResSo core
if bond_type == -1:
raise ValueError(
- "The bonded interaction with the id " + str(key) + " is not yet defined.")
+ f"The bonded interaction with the id {key} is not yet defined.")
# Find the appropriate class representing such a bond
bond_class = bonded_interaction_classes[bond_type]
diff --git a/src/python/espressomd/lb.pyx b/src/python/espressomd/lb.pyx
index bb774784fd1..9fd7697e50d 100644
--- a/src/python/espressomd/lb.pyx
+++ b/src/python/espressomd/lb.pyx
@@ -57,7 +57,7 @@ cdef class HydrodynamicInteraction(Actor):
return LBFluidRoutines(np.array(key))
else:
raise Exception(
- "%s is not a valid key. Should be a point on the nodegrid e.g. lbf[0,0,0]," % key)
+ f"{key} is not a valid key. Should be a point on the nodegrid e.g. lbf[0,0,0].")
# validate the given parameters on actor initialization
####################################################
diff --git a/src/python/espressomd/magnetostatics.pxd b/src/python/espressomd/magnetostatics.pxd
index 14ef94f12a0..73653277019 100644
--- a/src/python/espressomd/magnetostatics.pxd
+++ b/src/python/espressomd/magnetostatics.pxd
@@ -60,7 +60,7 @@ IF DP3M == 1:
cdef extern from "electrostatics_magnetostatics/p3m-dipolar.hpp":
int dp3m_set_params(double r_cut, int mesh, int cao, double alpha, double accuracy)
- void dp3m_set_tune_params(double r_cut, int mesh, int cao, double alpha, double accuracy, int n_interpol)
+ void dp3m_set_tune_params(double r_cut, int mesh, int cao, double alpha, double accuracy)
int dp3m_set_mesh_offset(double x, double y, double z)
int dp3m_set_eps(double eps)
int dp3m_adaptive_tune(char ** log)
diff --git a/src/python/espressomd/magnetostatics.pyx b/src/python/espressomd/magnetostatics.pyx
index 5b0f7c73e16..85d2d33c4f7 100644
--- a/src/python/espressomd/magnetostatics.pyx
+++ b/src/python/espressomd/magnetostatics.pyx
@@ -102,17 +102,14 @@ IF DP3M == 1:
raise ValueError("P3M r_cut has to be >=0")
if is_valid_type(self._params["mesh"], int):
- if self._params["mesh"] % 2 != 0 and self._params["mesh"] != -1:
- raise ValueError(
- "P3M requires an even number of mesh points in all directions")
+ pass
else:
check_type_or_throw_except(self._params["mesh"], 3, int,
- "P3M mesh has to be an integer or integer list of length 3")
- if (self._params["mesh"][0] % 2 != 0 and self._params["mesh"][0] != -1) or \
- (self._params["mesh"][1] % 2 != 0 and self._params["mesh"][1] != -1) or \
- (self._params["mesh"][2] % 2 != 0 and self._params["mesh"][2] != -1):
+ "DipolarP3M mesh has to be an integer or integer list of length 3")
+ if (self._params["mesh"][0] != self._params["mesh"][1]) or \
+ (self._params["mesh"][0] != self._params["mesh"][2]):
raise ValueError(
- "P3M requires an even number of mesh points in all directions")
+ "DipolarP3M requires a cubic box")
if not (self._params["cao"] >= -1 and self._params["cao"] <= 7):
raise ValueError(
@@ -124,13 +121,9 @@ IF DP3M == 1:
if self._params["epsilon"] == "metallic":
self._params["epsilon"] = 0.0
- if not (is_valid_type(self._params["epsilon"], float)
- or self._params["epsilon"] == "metallic"):
- raise ValueError("epsilon should be a double or 'metallic'")
-
- if not (self._params["inter"] == default_params["inter"]
- or self._params["inter"] > 0):
- raise ValueError("inter should be a positive integer")
+ check_type_or_throw_except(
+ self._params["epsilon"], 1, float,
+ "epsilon should be a double or 'metallic'")
if self._params["mesh_off"] != default_params["mesh_off"]:
check_type_or_throw_except(self._params["mesh_off"], 3, float,
@@ -138,15 +131,14 @@ IF DP3M == 1:
def valid_keys(self):
return ["prefactor", "alpha_L", "r_cut_iL", "mesh", "mesh_off",
- "cao", "inter", "accuracy", "epsilon", "cao_cut", "a", "ai",
- "alpha", "r_cut", "inter2", "cao3", "additional_mesh", "tune"]
+ "cao", "accuracy", "epsilon", "cao_cut", "a", "ai",
+ "alpha", "r_cut", "cao3", "additional_mesh", "tune"]
def required_keys(self):
return ["accuracy", ]
def default_params(self):
return {"cao": -1,
- "inter": -1,
"r_cut": -1,
"accuracy": -1,
"mesh": -1,
@@ -174,11 +166,9 @@ IF DP3M == 1:
dp3m_set_eps(self._params["epsilon"])
self.python_dp3m_set_tune_params(
self._params["r_cut"], self._params["mesh"],
- self._params["cao"], -1., self._params["accuracy"], self._params["inter"])
+ self._params["cao"], -1., self._params["accuracy"])
resp, log = self.python_dp3m_adaptive_tune()
- if resp:
- raise Exception(
- "failed to tune dipolar P3M parameters to required accuracy")
+ handle_errors("dipolar P3M tuning failed")
print(to_str(log))
self._params.update(self._get_params_from_es_core())
@@ -205,8 +195,7 @@ IF DP3M == 1:
cdef char * log = NULL
cdef int response
response = dp3m_adaptive_tune(& log)
- handle_errors(
- "dipolar P3M_init: k-space cutoff is larger than half of box dimension")
+ handle_errors("dipolar P3M tuning failed")
return response, log
def python_dp3m_set_params(self, p_r_cut, p_mesh, p_cao, p_alpha,
@@ -227,20 +216,21 @@ IF DP3M == 1:
dp3m_set_params(r_cut, mesh, cao, alpha, accuracy)
def python_dp3m_set_tune_params(self, p_r_cut, p_mesh, p_cao, p_alpha,
- p_accuracy, p_n_interpol):
+ p_accuracy):
cdef int mesh
cdef double r_cut
cdef int cao
cdef double alpha
cdef double accuracy
- cdef int n_interpol
r_cut = p_r_cut
cao = p_cao
alpha = p_alpha
accuracy = p_accuracy
- n_interpol = p_n_interpol
- mesh = p_mesh
- dp3m_set_tune_params(r_cut, mesh, cao, alpha, accuracy, n_interpol)
+ if hasattr(p_mesh, "__getitem__"):
+ mesh = p_mesh[0]
+ else:
+ mesh = p_mesh
+ dp3m_set_tune_params(r_cut, mesh, cao, alpha, accuracy)
IF DIPOLES == 1:
cdef class DipolarDirectSumCpu(MagnetostaticInteraction):
diff --git a/src/python/espressomd/particle_data.pyx b/src/python/espressomd/particle_data.pyx
index 4227de0fe71..0f11b0c04c1 100644
--- a/src/python/espressomd/particle_data.pyx
+++ b/src/python/espressomd/particle_data.pyx
@@ -1067,27 +1067,26 @@ cdef class ParticleHandle:
"""
if _partner in self.exclusions:
- raise Exception("Exclusion id {} already in exclusion list of particle {}".format(
- _partner, self._id))
+ raise Exception(
+ f"Exclusion id {_partner} already in exclusion list of particle {self._id}")
check_type_or_throw_except(
_partner, 1, int, "PID of partner has to be an int.")
if self._id == _partner:
raise Exception(
- "Cannot exclude of a particle with itself!\n->particle id %i, partner %i." % (self._id, _partner))
+ "Cannot exclude of a particle with itself!\n"
+ f"->particle id {self._id}, partner {_partner}.")
if change_exclusion(self._id, _partner, 0) == 1:
- raise Exception("Particle with id " +
- str(_partner) + " does not exist.")
+ raise Exception(f"Particle with id {_partner} does not exist.")
def delete_exclusion(self, _partner):
check_type_or_throw_except(
_partner, 1, int, "PID of partner has to be an int.")
if _partner not in self.exclusions:
- raise Exception("Particle with id " +
- str(_partner) + " is not in exclusion list.")
+ raise Exception(
+ f"Particle with id {_partner} is not in exclusion list.")
if change_exclusion(self._id, _partner, 1) == 1:
- raise Exception("Particle with id " +
- str(_partner) + " does not exist.")
+ raise Exception(f"Particle with id {_partner} does not exist.")
IF ENGINE:
property swimming:
@@ -1379,8 +1378,8 @@ cdef class ParticleHandle:
"""
if tuple(_bond) in self.bonds:
- raise Exception("Bond {} already exists on particle {}.".format(
- tuple(_bond), self._id))
+ raise Exception(
+ f"Bond {tuple(_bond)} already exists on particle {self._id}.")
bond = list(_bond) # As we will modify it
self.check_bond_or_throw_exception(bond)
@@ -1523,14 +1522,16 @@ cdef class _ParticleSliceImpl:
self.id_selection = np.array(slice_, dtype=int)
else:
raise TypeError(
- "ParticleSlice must be initialized with an instance of slice or range, or with a list, tuple, or ndarray of ints, but got {} of type {}".format((str(slice_), str(type(slice_)))))
+ f"ParticleSlice must be initialized with an instance of "
+ f"slice or range, or with a list, tuple, or ndarray of ints, "
+ f"but got {slice_} of type {type(slice_)}")
def _id_selection_from_slice(self, slice_):
"""Returns an ndarray of particle ids to be included in the
ParticleSlice for a given range or slice object.
"""
# Prevent negative bounds
- if (not slice_.start is None and slice_.start < 0) or\
+ if (not slice_.start is None and slice_.start < 0) or \
(not slice_.stop is None and slice_.stop < 0):
raise IndexError(
"Negative start and end ids are not supported on ParticleSlice")
@@ -1555,9 +1556,9 @@ cdef class _ParticleSliceImpl:
for pid in pid_list:
if not is_valid_type(pid, int):
raise TypeError(
- "Particle id must be an integer but got " + str(pid))
+ f"Particle id must be an integer but got {pid}")
if not particle_exists(pid):
- raise IndexError("Particle does not exist " + str(pid))
+ raise IndexError(f"Particle does not exist {pid}")
def __iter__(self):
return self._id_gen()
@@ -1606,9 +1607,9 @@ cdef class _ParticleSliceImpl:
pl = ParticleList()
for i in self.id_selection:
if pl.exists(i):
- res += str(pl[i]) + ", "
+ res += f"{pl[i]}, "
# Remove final comma
- return "ParticleSlice([" + res[:-2] + "])"
+ return f"ParticleSlice([{res[:-2]}])"
def update(self, P):
if "id" in P:
@@ -1662,7 +1663,7 @@ class ParticleSlice(_ParticleSliceImpl):
def __setattr__(self, name, value):
if name != "_chunk_size" and not hasattr(ParticleHandle, name):
raise AttributeError(
- "ParticleHandle does not have the attribute {}.".format(name))
+ f"ParticleHandle does not have the attribute {name}.")
super().__setattr__(name, value)
@@ -1796,7 +1797,7 @@ cdef class ParticleList:
P["id"] = get_maximal_particle_id() + 1
else:
if particle_exists(P["id"]):
- raise Exception("Particle %d already exists." % P["id"])
+ raise Exception(f"Particle {P['id']} already exists.")
# Prevent setting of contradicting attributes
IF DIPOLES:
@@ -2119,20 +2120,19 @@ def _add_particle_slice_properties():
elif np.shape(values)[0] == N:
set_slice_one_for_each(particle_slice, attribute, values)
else:
- raise Exception("Shape of value (%s) does not broadcast to shape of attribute (%s)." % (
- np.shape(values), target_shape))
+ raise Exception(
+ f"Value shape {np.shape(values)} does not broadcast to attribute shape {target_shape}.")
return
else: # fixed length vector quantity
-
if target_shape == np.shape(values):
set_slice_one_for_all(particle_slice, attribute, values)
elif target_shape == tuple(np.shape(values)[1:]) and np.shape(values)[0] == N:
set_slice_one_for_each(particle_slice, attribute, values)
else:
- raise Exception("Shape of value (%s) does not broadcast to shape of attribute (%s)." % (
- np.shape(values), target_shape))
+ raise Exception(
+ f"Value shape {np.shape(values)} does not broadcast to attribute shape {target_shape}.")
return
diff --git a/src/python/espressomd/polymer.pyx b/src/python/espressomd/polymer.pyx
index e6884aedbca..8a312171bcc 100644
--- a/src/python/espressomd/polymer.pyx
+++ b/src/python/espressomd/polymer.pyx
@@ -131,7 +131,7 @@ def linear_polymer_positions(**kwargs):
for k in kwargs:
if k not in valid_keys:
- raise ValueError("Unknown parameter '%s'" % k)
+ raise ValueError(f"Unknown parameter '{k}'")
params[k] = kwargs[k]
for k in required_keys:
diff --git a/src/python/espressomd/reaction_ensemble.pyx b/src/python/espressomd/reaction_ensemble.pyx
index 4f6f0698919..474a03b660a 100644
--- a/src/python/espressomd/reaction_ensemble.pyx
+++ b/src/python/espressomd/reaction_ensemble.pyx
@@ -426,7 +426,7 @@ cdef class ReactionEnsemble(ReactionAlgorithm):
if k in self._valid_keys():
self._params[k] = kwargs[k]
else:
- raise KeyError("%s is not a valid key" % k)
+ raise KeyError(f"{k} is not a valid key")
self._set_params_in_es_core()
@@ -449,7 +449,7 @@ cdef class ConstantpHEnsemble(ReactionAlgorithm):
if k in self._valid_keys():
self._params[k] = kwargs[k]
else:
- raise KeyError("%s is not a valid key" % k)
+ raise KeyError(f"{k} is not a valid key")
self._set_params_in_es_core()
@@ -495,7 +495,7 @@ cdef class WangLandauReactionEnsemble(ReactionAlgorithm):
if k in self._valid_keys():
self._params[k] = kwargs[k]
else:
- raise KeyError("%s is not a valid key" % k)
+ raise KeyError(f"{k} is not a valid key")
self.WLRptr.reset(new CWangLandauReactionEnsemble(int(self._params["seed"])))
self.RE = self.WLRptr.get()
@@ -539,7 +539,7 @@ cdef class WangLandauReactionEnsemble(ReactionAlgorithm):
if k in self._valid_keys_add_collective_variable_degree_of_association():
self._params[k] = kwargs[k]
else:
- raise KeyError("%s is not a valid key" % k)
+ raise KeyError(f"{k} is not a valid key")
for k in self._required_keys_add_collective_variable_degree_of_association():
if k not in kwargs:
@@ -591,7 +591,7 @@ cdef class WangLandauReactionEnsemble(ReactionAlgorithm):
if k in self._valid_keys_add_collective_variable_potential_energy():
self._params[k] = kwargs[k]
else:
- raise KeyError("%s is not a valid key" % k)
+ raise KeyError(f"{k} is not a valid key")
for k in self._required_keys_add_collective_variable_potential_energy():
if k not in kwargs:
@@ -633,7 +633,7 @@ cdef class WangLandauReactionEnsemble(ReactionAlgorithm):
if k in self._valid_keys_set_wang_landau_parameters():
self._params[k] = kwargs[k]
else:
- raise KeyError("%s is not a valid key" % k)
+ raise KeyError(f"{k} is not a valid key")
deref(self.WLRptr).final_wang_landau_parameter = self._params[
"final_wang_landau_parameter"]
@@ -760,7 +760,7 @@ cdef class WidomInsertion(ReactionAlgorithm):
if k in self._valid_keys():
self._params[k] = kwargs[k]
else:
- raise KeyError("%s is not a valid key" % k)
+ raise KeyError(f"{k} is not a valid key")
self._set_params_in_es_core()
diff --git a/src/python/espressomd/script_interface.pxd b/src/python/espressomd/script_interface.pxd
index d2c7de312e7..e567d86b1fa 100644
--- a/src/python/espressomd/script_interface.pxd
+++ b/src/python/espressomd/script_interface.pxd
@@ -1,5 +1,5 @@
#
-# Copyright (C) 2013-2019 The ESPResSo project
+# Copyright (C) 2013-2020 The ESPResSo project
#
# This file is part of ESPResSo.
#
@@ -22,15 +22,14 @@ from libcpp.unordered_map cimport unordered_map
from libcpp.vector cimport vector
from libcpp.string cimport string
from libcpp.memory cimport shared_ptr
-from libcpp.memory cimport weak_ptr
from libcpp cimport bool
from boost cimport string_ref
-from .utils cimport Span
+from .utils cimport Span, Factory
+from .communication cimport MpiCallbacks
cdef extern from "script_interface/ScriptInterface.hpp" namespace "ScriptInterface":
- void initialize()
cdef cppclass Variant:
Variant()
Variant(const Variant & )
@@ -40,48 +39,37 @@ cdef extern from "script_interface/ScriptInterface.hpp" namespace "ScriptInterfa
bool is_none(const Variant &)
ctypedef unordered_map[string, Variant] VariantMap
-cdef extern from "script_interface/get_value.hpp" namespace "ScriptInterface":
- T get_value[T](const Variant T)
-
-cdef extern from "script_interface/ScriptInterface.hpp" namespace "ScriptInterface":
- cdef cppclass ObjectId:
- ObjectId()
- string to_string()
- bool operator == (const ObjectId & rhs)
- bool operator != (const ObjectId & rhs)
-
Variant make_variant[T](const T & x)
- cdef cppclass ScriptInterfaceBase:
- const string name()
- void construct(const VariantMap &) except +
+ cdef cppclass ObjectHandle:
VariantMap get_parameters() except +
Span[const string_ref] valid_parameters() except +
Variant get_parameter(const string & name) except +
void set_parameter(const string & name, const Variant & value) except +
Variant call_method(const string & name, const VariantMap & parameters) except +
- ObjectId id() except +
void set_state(map[string, Variant]) except +
map[string, Variant] get_state() except +
- string serialize() except +
+ string_ref name()
- @staticmethod
- shared_ptr[ScriptInterfaceBase] unserialize(const string & state) except +
-
-cdef extern from "script_interface/ScriptInterface.hpp" namespace "ScriptInterface::ScriptInterfaceBase":
+cdef extern from "script_interface/ContextManager.hpp" namespace "ScriptInterface::ContextManager":
cdef cppclass CreationPolicy:
pass
- shared_ptr[ScriptInterfaceBase] make_shared(const string & name, CreationPolicy policy) except +
- weak_ptr[ScriptInterfaceBase] get_instance(ObjectId id) except +
-cdef extern from "script_interface/ScriptInterface.hpp" namespace "ScriptInterface::ScriptInterfaceBase::CreationPolicy":
+cdef extern from "script_interface/ContextManager.hpp" namespace "ScriptInterface::ContextManager::CreationPolicy":
CreationPolicy LOCAL
CreationPolicy GLOBAL
-cdef variant_to_python_object(const Variant & value) except +
-cdef Variant python_object_to_variant(value)
+cdef extern from "script_interface/ContextManager.hpp" namespace "ScriptInterface":
+ cdef cppclass ContextManager:
+ ContextManager(MpiCallbacks & , const Factory[ObjectHandle] & )
+ shared_ptr[ObjectHandle] make_shared(CreationPolicy, const string &, const VariantMap) except +
+ shared_ptr[ObjectHandle] deserialize(const string &) except +
+ string serialize(const ObjectHandle *) except +
+
+cdef extern from "script_interface/initialize.hpp" namespace "ScriptInterface":
+ void initialize(Factory[ObjectHandle] *)
+
+cdef extern from "script_interface/get_value.hpp" namespace "ScriptInterface":
+ T get_value[T](const Variant T)
-cdef class PScriptInterface:
- cdef shared_ptr[ScriptInterfaceBase] sip
- cdef set_sip(self, shared_ptr[ScriptInterfaceBase] sip)
- cdef VariantMap _sanitize_params(self, in_params) except *
+cdef void init(MpiCallbacks &)
diff --git a/src/python/espressomd/script_interface.pyx b/src/python/espressomd/script_interface.pyx
index 593fb1f7d02..a2e6dfcfe23 100644
--- a/src/python/espressomd/script_interface.pyx
+++ b/src/python/espressomd/script_interface.pyx
@@ -18,17 +18,19 @@ import numpy as np
from .utils import to_char_pointer, to_str
from .utils cimport Vector3d, make_array_locked, handle_errors
-cdef class PObjectId:
- """Python interface to a core ObjectId object."""
+from libcpp.memory cimport make_shared
- cdef ObjectId id
+cdef shared_ptr[ContextManager] _om
- def __richcmp__(PObjectId a, PObjectId b, op):
+cdef class PObjectRef:
+ def __richcmp__(PObjectRef a, PObjectRef b, op):
if op == 2:
- return a.id == b.id
+ return a.sip == b.sip
else:
raise NotImplementedError
+ cdef shared_ptr[ObjectHandle] sip
+
cdef class PScriptInterface:
"""
@@ -46,23 +48,29 @@ cdef class PScriptInterface:
Name of the core class to instantiate (method 1).
\*\*kwargs
Parameters for the core class constructor (method 1).
- oid : :class:`PObjectId`
+ sip : :class:`PObjectRef`
Object id of an existing core object (method 2).
+
policy : :obj:`str`, \{'GLOBAL', 'LOCAL'\}
Creation policy.
Attributes
----------
- sip: shared_ptr
+ sip: :class:`PObjectRef`
Pointer to a ScriptInterface object in the core.
policy_: :obj:`str`
Creation policy.
"""
- def __init__(self, name=None, policy="GLOBAL", oid=None, **kwargs):
+ cdef shared_ptr[ObjectHandle] sip
+ cdef set_sip(self, shared_ptr[ObjectHandle] sip)
+ cdef VariantMap _sanitize_params(self, in_params) except *
+
+ def __init__(self, name=None, policy="GLOBAL", sip=None, **kwargs):
cdef CreationPolicy policy_
+ cdef PObjectRef sip_
if policy == "GLOBAL":
policy_ = GLOBAL
@@ -71,15 +79,20 @@ cdef class PScriptInterface:
else:
raise Exception("Unknown policy '{}'.".format(policy))
- if oid:
- self.set_sip_via_oid(oid)
+ if sip:
+ sip_ = sip
+ self.sip = sip_.sip
else:
- self.set_sip(make_shared(to_char_pointer(name), policy_))
- self.sip.get().construct(self._sanitize_params(kwargs))
+ global _om
+ self.set_sip(
+ _om.get().make_shared(
+ policy_,
+ to_char_pointer(name),
+ self._sanitize_params(kwargs)))
def __richcmp__(a, b, op):
if op == 2:
- return a.id() == b.id()
+ return a.get_sip() == b.get_sip()
else:
raise NotImplementedError
@@ -89,25 +102,21 @@ cdef class PScriptInterface:
def _valid_parameters(self):
return [to_str(p.data()) for p in self.sip.get().valid_parameters()]
- cdef set_sip(self, shared_ptr[ScriptInterfaceBase] sip):
- self.sip = sip
+ def get_sip(self):
+ """
+ Get pointer to the core object.
+ """
+
+ ret = PObjectRef()
+ ret.sip = self.sip
+ return ret
- def set_sip_via_oid(self, PObjectId id):
- """Set the shared_ptr to an existing core ScriptInterface object via
- its object id.
+ cdef set_sip(self, shared_ptr[ObjectHandle] sip):
"""
- oid = id.id
- try:
- ptr = get_instance(oid).lock()
- self.set_sip(ptr)
- except Exception:
- raise Exception("Could not get sip for given_id")
-
- def id(self):
- """Return the core class object id (:class:`PObjectId`)."""
- oid = PObjectId()
- oid.id = self.sip.get().id()
- return oid
+ Set the shared_ptr to an existing core object.
+ """
+
+ self.sip = sip
def call_method(self, method, **kwargs):
"""
@@ -119,7 +128,6 @@ cdef class PScriptInterface:
Name of the core method.
\*\*kwargs
Arguments for the method.
-
"""
cdef VariantMap parameters
@@ -134,27 +142,23 @@ cdef class PScriptInterface:
def name(self):
"""Return name of the core class."""
- return to_str(self.sip.get().name())
+ return to_str(self.sip.get().name().data())
def _serialize(self):
- return self.sip.get().serialize()
+ global _om
+ return _om.get().serialize(self.sip.get())
def _unserialize(self, state):
- cdef shared_ptr[ScriptInterfaceBase] so_ptr = ScriptInterfaceBase.unserialize(state)
+ global _om
+ cdef shared_ptr[ObjectHandle] so_ptr = _om.get().deserialize(state)
self.set_sip(so_ptr)
cdef VariantMap _sanitize_params(self, in_params) except *:
cdef VariantMap out_params
- cdef Variant v
-
- valid_params = self._valid_parameters()
for pname in in_params:
- if pname in valid_params:
- out_params[to_char_pointer(pname)] = python_object_to_variant(
- in_params[pname])
- else:
- raise RuntimeError("Unknown parameter '{}'".format(pname))
+ out_params[to_char_pointer(pname)] = python_object_to_variant(
+ in_params[pname])
return out_params
@@ -178,19 +182,19 @@ cdef class PScriptInterface:
cdef Variant python_object_to_variant(value):
"""Convert Python objects to C++ Variant objects."""
- cdef Variant v
+
cdef vector[Variant] vec
- cdef PObjectId oid
+ cdef PObjectRef oref
if value is None:
return Variant()
- # The order is important, the object character should be preserved
- # even if the PScriptInterface derived class is iterable.
+ # The order is important, the object character should
+ # be preserved even if the PScriptInterface derived class
+ # is iterable.
if isinstance(value, PScriptInterface):
- # Map python object to id
- oid = value.id()
- return make_variant[ObjectId](oid.id)
+ oref = value.get_sip()
+ return make_variant(oref.sip)
elif hasattr(value, '__iter__') and not(type(value) == str):
for e in value:
vec.push_back(python_object_to_variant(e))
@@ -208,9 +212,9 @@ cdef Variant python_object_to_variant(value):
cdef variant_to_python_object(const Variant & value) except +:
"""Convert C++ Variant objects to Python objects."""
- cdef ObjectId oid
+
cdef vector[Variant] vec
- cdef shared_ptr[ScriptInterfaceBase] ptr
+ cdef shared_ptr[ObjectHandle] ptr
if is_none(value):
return None
if is_type[bool](value):
@@ -227,18 +231,12 @@ cdef variant_to_python_object(const Variant & value) except +:
return get_value[vector[double]](value)
if is_type[Vector3d](value):
return make_array_locked(get_value[Vector3d](value))
- if is_type[ObjectId](value):
+ if is_type[shared_ptr[ObjectHandle]](value):
# Get the id and build a corresponding object
- oid = get_value[ObjectId](value)
-
- # ObjectId is nullable, and the default id corresponds to "null".
- if oid != ObjectId():
- ptr = get_instance(oid).lock()
+ ptr = get_value[shared_ptr[ObjectHandle]](value)
- if not ptr:
- raise Exception("Object failed to exist.")
-
- so_name = to_str(ptr.get().name())
+ if ptr:
+ so_name = to_str(ptr.get().name().data())
if not so_name:
raise Exception(
"Script object without name returned from the core")
@@ -251,10 +249,10 @@ cdef variant_to_python_object(const Variant & value) except +:
# for the script object name
pclass = ScriptInterfaceHelper
- poid = PObjectId()
- poid.id = ptr.get().id()
+ pptr = PObjectRef()
+ pptr.sip = ptr
- return pclass(oid=poid)
+ return pclass(sip=pptr)
else:
return None
if is_type[vector[Variant]](value):
@@ -270,24 +268,18 @@ cdef variant_to_python_object(const Variant & value) except +:
def _unpickle_so_class(so_name, state):
- cdef shared_ptr[ScriptInterfaceBase] sip = ScriptInterfaceBase.unserialize(state)
-
- poid = PObjectId()
- poid.id = sip.get().id()
+ cdef PObjectRef so_ptr
+ so_ptr = PObjectRef()
+ global _om
+ so_ptr.sip = _om.get().deserialize(state)
- so = _python_class_by_so_name[so_name](oid=poid)
+ so = _python_class_by_so_name[so_name](sip=so_ptr)
so.define_bound_methods()
return so
class ScriptInterfaceHelper(PScriptInterface):
-
- """
- Base class from which to derive most interfaces to core ScriptInterface
- classes.
- """
-
_so_name = None
_so_bind_methods = ()
_so_creation_policy = "GLOBAL"
@@ -332,7 +324,6 @@ class ScriptInterfaceHelper(PScriptInterface):
class ScriptObjectRegistry(ScriptInterfaceHelper):
-
"""
Base class for container-like classes such as
:class:`~espressomd.constraints.Constraints` and
@@ -377,15 +368,22 @@ _python_class_by_so_name = {}
def script_interface_register(c):
"""
Decorator used to register script interface classes.
- This will store a name-to-class relationship in a registry, so that
- parameters of type object can be instantiated as the correct python class
+ This will store a name<->class relationship in a registry, so that
+ parameters of type object can be instantiated as the correct python class.
"""
+
if not hasattr(c, "_so_name"):
raise Exception("Python classes representing a script object must "
"define an _so_name attribute at class level")
+
_python_class_by_so_name[c._so_name] = c
return c
-def init():
- initialize()
+cdef void init(MpiCallbacks & cb):
+ cdef Factory[ObjectHandle] f
+
+ initialize( & f)
+
+ global _om
+ _om = make_shared[ContextManager](cb, f)
diff --git a/src/python/espressomd/system.pyx b/src/python/espressomd/system.pyx
index 6de9d7bb60f..5ddc28cc110 100644
--- a/src/python/espressomd/system.pyx
+++ b/src/python/espressomd/system.pyx
@@ -112,7 +112,7 @@ cdef class System:
def __init__(self, **kwargs):
global _system_created
- if (not _system_created):
+ if not _system_created:
self.globals = Globals()
if 'box_l' not in kwargs:
raise ValueError("Required argument box_l not provided.")
@@ -123,7 +123,7 @@ cdef class System:
System.__setattr__(self, arg, kwargs.get(arg))
else:
raise ValueError(
- "Property {} can not be set via argument to System class.".format(arg))
+ f"Property {arg} can not be set via argument to System class.")
self.actors = Actors()
self.analysis = Analysis(self)
self.auto_update_accumulators = AutoUpdateAccumulators()
diff --git a/src/python/espressomd/thermostat.pxd b/src/python/espressomd/thermostat.pxd
index 622cda0f333..84e6eeab925 100644
--- a/src/python/espressomd/thermostat.pxd
+++ b/src/python/espressomd/thermostat.pxd
@@ -52,6 +52,10 @@ cdef extern from "thermostat.hpp":
double gamma0
double gammav
+ langevin_thermostat_struct langevin
+ brownian_thermostat_struct brownian
+ npt_iso_thermostat_struct npt_iso
+
void langevin_set_rng_state(stdint.uint64_t counter)
void brownian_set_rng_state(stdint.uint64_t counter)
void npt_iso_set_rng_state(stdint.uint64_t counter)
@@ -78,12 +82,6 @@ cdef extern from "stokesian_dynamics/sd_interface.hpp":
void set_sd_seed(size_t seed)
size_t get_sd_seed()
-cdef extern from "script_interface/Globals.hpp":
- # links intern C-struct with python object
- cdef extern langevin_thermostat_struct langevin
- cdef extern brownian_thermostat_struct brownian
- cdef extern npt_iso_thermostat_struct npt_iso
-
cdef extern from "npt.hpp":
ctypedef struct nptiso_struct:
double p_ext
diff --git a/src/python/espressomd/utils.pxd b/src/python/espressomd/utils.pxd
index 0211f9c050c..4b6f51f170d 100644
--- a/src/python/espressomd/utils.pxd
+++ b/src/python/espressomd/utils.pxd
@@ -99,3 +99,7 @@ cdef extern from "utils/Vector.hpp" namespace "Utils":
cdef make_array_locked(Vector3d)
cdef Vector3d make_Vector3d(a)
+
+cdef extern from "utils/Factory.hpp" namespace "Utils":
+ cdef cppclass Factory[T]:
+ pass
diff --git a/src/python/espressomd/utils.pyx b/src/python/espressomd/utils.pyx
index 492c158206d..5246b29b1ee 100644
--- a/src/python/espressomd/utils.pyx
+++ b/src/python/espressomd/utils.pyx
@@ -33,23 +33,26 @@ cpdef check_type_or_throw_except(x, n, t, msg):
# Check whether x is an array/list/tuple or a single value
if n > 1:
if hasattr(x, "__getitem__"):
+ if len(x) != n:
+ raise ValueError(
+ msg + f" -- {len(x)} values were given but {n} were expected.")
for i in range(len(x)):
if not isinstance(x[i], t):
if not ((t == float and is_valid_type(x[i], int))
or (t == float and issubclass(type(x[i]), np.integer))) \
and not (t == int and issubclass(type(x[i]), np.integer)):
raise ValueError(
- msg + " -- Item " + str(i) + " was of type " + type(x[i]).__name__)
+ msg + f" -- Item {i} was of type {type(x[i]).__name__}")
else:
# if n>1, but the user passed a single value, also throw exception
raise ValueError(
- msg + " -- A single value was given but " + str(n) + " were expected.")
+ msg + f" -- A single value was given but {n} were expected.")
else:
# N=1 and a single value
if not isinstance(x, t):
if not (t == float and is_valid_type(x, int)) and not (
t == int and issubclass(type(x), np.integer)):
- raise ValueError(msg + " -- Got an " + type(x).__name__)
+ raise ValueError(msg + f" -- Got an {type(x).__name__}")
cdef np.ndarray create_nparray_from_double_array(double * x, int len_x):
@@ -93,14 +96,16 @@ cdef check_range_or_except(D, name, v_min, incl_min, v_max, incl_max):
or (not incl_min and not all(v > v_min for v in x)))) or \
(v_max != "inf" and ((incl_max and not all(v <= v_max for v in x))
or (not incl_max and not all(v < v_max for v in x)))):
- raise ValueError("In " + name + ": Some values in " + str(x) + "are out of range " +
- ("[" if incl_min else "]") + str(v_min) + "," + str(v_max) + ("]" if incl_max else "["))
+ raise ValueError(f"In {name}: Some values in {x} are out of"
+ f" range {'[' if incl_min else ']'}{v_min},"
+ f"{v_max}{']' if incl_max else '['}")
# Single Value
else:
if (v_min != "inf" and ((incl_min and x < v_min) or (not incl_min and x <= v_min)) or
v_max != "inf" and ((incl_max and x > v_max) or (not incl_max and x >= v_max))):
- raise ValueError("In " + name + ": Value " + str(x) + " is out of range " + ("[" if incl_min else "]") +
- str(v_min) + "," + str(v_max) + ("]" if incl_max else "["))
+ raise ValueError(f"In {name}: Value {x} is out of"
+ f" range {'[' if incl_min else ']'}{v_min},"
+ f"{v_max}{']' if incl_max else '['}")
def to_char_pointer(s):
@@ -286,10 +291,9 @@ def requires_experimental_features(reason):
def exception_raiser(self, *args, **kwargs):
raise Exception(
- "Class " +
- self.__class__.__name__ +
- " is experimental. Define EXPERIMENTAL_FEATURES in myconfig.hpp to use it.\nReason: " +
- reason)
+ f"Class {self.__class__.__name__} is experimental. Define "
+ "EXPERIMENTAL_FEATURES in myconfig.hpp to use it.\n"
+ f"Reason: {reason}")
def modifier(cls):
cls.__init__ = exception_raiser
diff --git a/src/python/espressomd/version.pyx b/src/python/espressomd/version.pyx
index 63c440befa5..6b93e7a39fa 100644
--- a/src/python/espressomd/version.pyx
+++ b/src/python/espressomd/version.pyx
@@ -33,7 +33,7 @@ def minor():
def friendly():
"""Dot version of the version.
"""
- return "{}.{}".format(major(), minor())
+ return f"{major()}.{minor()}"
def git_branch():
diff --git a/src/script_interface/CMakeLists.txt b/src/script_interface/CMakeLists.txt
index f08e6f1e6c3..01e1d345744 100644
--- a/src/script_interface/CMakeLists.txt
+++ b/src/script_interface/CMakeLists.txt
@@ -1,5 +1,5 @@
-add_library(ScriptInterface SHARED initialize.cpp ScriptInterfaceBase.cpp
- ParallelScriptInterface.cpp)
+add_library(ScriptInterface SHARED initialize.cpp ObjectHandle.cpp
+ GlobalContext.cpp ContextManager.cpp)
add_subdirectory(accumulators)
add_subdirectory(collision_detection)
@@ -22,3 +22,5 @@ target_link_libraries(
PRIVATE cxx_interface)
target_include_directories(ScriptInterface PUBLIC ${CMAKE_SOURCE_DIR}/src)
+
+add_subdirectory(tests)
diff --git a/src/script_interface/ComFixed.hpp b/src/script_interface/ComFixed.hpp
index 83f2a13eacb..3653b8232b5 100644
--- a/src/script_interface/ComFixed.hpp
+++ b/src/script_interface/ComFixed.hpp
@@ -22,8 +22,7 @@
#ifndef SCRIPT_INTERFACE_COM_FIXED_HPP
#define SCRIPT_INTERFACE_COM_FIXED_HPP
-#include "ScriptInterfaceBase.hpp"
-#include "script_interface/auto_parameters/AutoParameters.hpp"
+#include "script_interface/ScriptInterface.hpp"
#include "core/comfixed_global.hpp"
diff --git a/src/script_interface/Context.hpp b/src/script_interface/Context.hpp
new file mode 100644
index 00000000000..fac100153aa
--- /dev/null
+++ b/src/script_interface/Context.hpp
@@ -0,0 +1,102 @@
+/*
+ * Copyright (C) 2020 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 .
+ */
+#ifndef ESPRESSO_SCRIPT_INTERFACE_CONTEXT_HPP
+#define ESPRESSO_SCRIPT_INTERFACE_CONTEXT_HPP
+
+/** @file
+ *
+ * @ref ScriptInterface::Context decorates @ref ScriptInterface::ObjectHandle
+ * objects with a context: a creation policy (local object, local object
+ * with remote copies) and a communication facility to synchronize an object
+ * on the head node with remote copies (serialization, callback mechanism).
+ */
+
+#include "ObjectHandle.hpp"
+#include "Variant.hpp"
+
+#include
+
+#include
+
+namespace ScriptInterface {
+/**
+ * @brief Context of an object handle.
+ *
+ * Each instance of @ref ObjectHandle can have an
+ * attached context, which can e.g. synchronize
+ * distributed copies of the instance. The context
+ * does also provide facilities for serializing
+ * Objects into a string representation.
+ */
+class Context : public std::enable_shared_from_this {
+public:
+ /**
+ * @brief Call method on remote instances
+ *
+ * @param self Internal identifier of the instance
+ * @param name Name of the method to call
+ * @param arguments Arguments to the call
+ */
+ virtual void notify_call_method(const ObjectHandle *self,
+ std::string const &name,
+ VariantMap const &arguments) = 0;
+
+ /**
+ * @brief Set a parameter on remote instances
+ *
+ * @param self Internal identifier of the instance to be modified
+ * @param name Name of the parameter to change
+ * @param value Value to set it to
+ */
+ virtual void notify_set_parameter(const ObjectHandle *self,
+ std::string const &name,
+ Variant const &value) = 0;
+
+ /**
+ * @brief Get a new reference counted instance of a script interface by
+ * name.
+ *
+ * Objects created thru a Context get shared ownership of that context,
+ * e.g. the lifetime of the context is at least as long as the objects
+ * created by it. Therefore the object can always assume that the context
+ * is present.
+ */
+ virtual std::shared_ptr
+ make_shared(std::string const &name, const VariantMap ¶meters) = 0;
+
+protected:
+ /**
+ * @brief Set the context of an object to this.
+ *
+ * @param o Object to set the context for.
+ */
+ void set_context(ObjectHandle *o) { o->m_context = this->shared_from_this(); }
+
+public:
+ /**
+ * @brief Get the class name for an ObjectHandle instance.
+ *
+ * This returns the name by which the object can be created.
+ */
+ virtual boost::string_ref name(const ObjectHandle *o) const = 0;
+
+ virtual ~Context() = default;
+};
+} // namespace ScriptInterface
+#endif // ESPRESSO_CONTEXT_HPP
diff --git a/src/script_interface/ContextManager.cpp b/src/script_interface/ContextManager.cpp
new file mode 100644
index 00000000000..c0b7922cca0
--- /dev/null
+++ b/src/script_interface/ContextManager.cpp
@@ -0,0 +1,64 @@
+/*
+ * Copyright (C) 2020 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 .
+ */
+#include "ContextManager.hpp"
+
+#include "GlobalContext.hpp"
+#include "LocalContext.hpp"
+
+#include
+
+namespace ScriptInterface {
+std::shared_ptr
+ContextManager::make_shared(CreationPolicy policy, std::string const &name,
+ const VariantMap ¶meters) {
+ return context(policy)->make_shared(name, parameters);
+}
+
+std::shared_ptr
+ContextManager::deserialize(std::string const &state_) {
+ auto const state =
+ Utils::unpack>(state_);
+
+ auto ctx = context(state.first);
+ assert(ctx);
+
+ return ObjectHandle::deserialize(state.second, *ctx);
+}
+
+std::string ContextManager::serialize(const ObjectHandle *o) const {
+ /* We treat objects without a context as local. */
+ auto ctx = o->context() ? o->context() : m_local_context.get();
+
+ return Utils::pack(std::make_pair(policy(ctx), o->serialize()));
+}
+
+ContextManager::ContextManager(Communication::MpiCallbacks &callbacks,
+ const Utils::Factory &factory) {
+ auto local_context = std::make_shared(factory);
+
+ /* If there is only one node, we can treat all objects as local, and thus
+ * never invoke any callback. */
+ m_global_context =
+ (callbacks.comm().size() > 1)
+ ? std::make_shared(callbacks, local_context)
+ : std::static_pointer_cast(local_context);
+
+ m_local_context = std::move(local_context);
+}
+} // namespace ScriptInterface
diff --git a/src/script_interface/ContextManager.hpp b/src/script_interface/ContextManager.hpp
new file mode 100644
index 00000000000..e76309f18d6
--- /dev/null
+++ b/src/script_interface/ContextManager.hpp
@@ -0,0 +1,125 @@
+/*
+ * Copyright (C) 2020 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 .
+ */
+#ifndef ESPRESSO_CONTEXTMANAGER_HPP
+#define ESPRESSO_CONTEXTMANAGER_HPP
+
+/** @file
+ *
+ * @ref ScriptInterface::ContextManager manages object creation with policies
+ * @ref ScriptInterface::ContextManager::CreationPolicy "CreationPolicy".
+ * Object creation is delegated to @ref ScriptInterface::GlobalContext and
+ * @ref ScriptInterface::LocalContext. @ref ScriptInterface::ContextManager
+ * serves as their public interface. If there is only 1 MPI rank, no
+ * communication takes place and all objects are created locally via
+ * @ref ScriptInterface::LocalContext, including those with policy
+ * @ref ScriptInterface::ContextManager::CreationPolicy::GLOBAL "GLOBAL".
+ *
+ * Implementation in @ref ContextManager.cpp.
+ */
+
+#include "Context.hpp"
+#include "Variant.hpp"
+
+#include "MpiCallbacks.hpp"
+
+#include
+
+#include
+#include
+
+namespace ScriptInterface {
+
+/**
+ * @brief Manage object contexts.
+ *
+ * This owns object contexts and allows for
+ * creation and serialization of objects preserving
+ * their context.
+ */
+class ContextManager {
+ std::shared_ptr m_local_context;
+ std::shared_ptr m_global_context;
+
+public:
+ /** Labels for context */
+ enum class CreationPolicy {
+ /** Corresponding to @c LocalContext */
+ LOCAL,
+ /** Corresponding to @c GlobalContext */
+ GLOBAL
+ };
+
+ ContextManager(Communication::MpiCallbacks &callbacks,
+ const Utils::Factory &factory);
+
+ /**
+ * @brief Get a new reference counted instance of a script interface by
+ * name.
+ */
+ std::shared_ptr make_shared(CreationPolicy policy,
+ std::string const &name,
+ const VariantMap ¶meters);
+
+ /**
+ * @brief Get a new reference counted instance of a script interface from
+ * a serialized state.
+ */
+ std::shared_ptr deserialize(std::string const &state_);
+
+ /**
+ * @brief Serialize a script interface object into a binary representation.
+ */
+ std::string serialize(const ObjectHandle *o) const;
+
+private:
+ /**
+ * @brief Map policy to context.
+ *
+ * Inverse of policy.
+ */
+ Context *context(CreationPolicy policy) const {
+ switch (policy) {
+ case CreationPolicy::LOCAL:
+ return assert(m_local_context), m_local_context.get();
+ case CreationPolicy::GLOBAL:
+ return assert(m_global_context), m_global_context.get();
+ default:
+ throw std::runtime_error("Unknown context type.");
+ }
+ }
+
+ /**
+ * @brief Map context to policy.
+ *
+ * Inverse of context.
+ */
+ CreationPolicy policy(Context *c) const {
+ if (c == m_local_context.get()) {
+ return CreationPolicy::LOCAL;
+ }
+ if (c == m_global_context.get()) {
+ return CreationPolicy::GLOBAL;
+ }
+
+ throw std::runtime_error("Invalid context.");
+ }
+};
+} // namespace ScriptInterface
+
+#endif // ESPRESSO_CONTEXTMANAGER_HPP
diff --git a/src/script_interface/Exception.hpp b/src/script_interface/Exception.hpp
new file mode 100644
index 00000000000..21cf22f66cd
--- /dev/null
+++ b/src/script_interface/Exception.hpp
@@ -0,0 +1,37 @@
+/*
+ * Copyright (C) 2020 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 .
+ */
+#ifndef ESPRESSO_EXCEPTIONS_HPP
+#define ESPRESSO_EXCEPTIONS_HPP
+
+#include
+#include
+
+namespace ScriptInterface {
+struct Exception : public std::exception {
+ explicit Exception(const char *msg) : message(msg) {}
+ explicit Exception(std::string msg) : message(std::move(msg)) {}
+
+ const char *what() const noexcept override { return message.c_str(); }
+
+private:
+ std::string message;
+};
+} // namespace ScriptInterface
+
+#endif // ESPRESSO_EXCEPTIONS_HPP
diff --git a/src/script_interface/GlobalContext.cpp b/src/script_interface/GlobalContext.cpp
new file mode 100644
index 00000000000..47a64864f8e
--- /dev/null
+++ b/src/script_interface/GlobalContext.cpp
@@ -0,0 +1,117 @@
+/*
+ * Copyright (C) 2020 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 .
+ */
+
+/** @file
+ *
+ * Infrastructure to synchronize objects created on the head node with
+ * their corresponding remote copies. Methods of script interface
+ * classes may throw exceptions of type @ref ScriptInterface::Exception.
+ * These exceptions will halt the flow of the program on the head node.
+ * The same exceptions will be thrown in the remote copies but will
+ * be silenced, since they are redundant. Other types of exceptions
+ * are not silenced.
+ *
+ * Implementation of @ref GlobalContext.hpp.
+ */
+
+#include "GlobalContext.hpp"
+#include "Exception.hpp"
+#include "ObjectHandle.hpp"
+#include "packed_variant.hpp"
+
+#include
+
+namespace ScriptInterface {
+void GlobalContext::make_handle(ObjectId id, const std::string &name,
+ const PackedMap ¶meters) {
+ try {
+ ObjectRef so = m_node_local_context->make_shared(
+ name, unpack(parameters, m_local_objects));
+
+ m_local_objects.emplace(std::make_pair(id, std::move(so)));
+ } catch (Exception const &) {
+ }
+}
+
+void GlobalContext::remote_make_handle(ObjectId id, const std::string &name,
+ const VariantMap ¶meters) {
+ cb_make_handle(id, name, pack(parameters));
+}
+
+void GlobalContext::set_parameter(ObjectId id, std::string const &name,
+ PackedVariant const &value) {
+ try {
+ m_local_objects.at(id)->set_parameter(name, unpack(value, m_local_objects));
+ } catch (Exception const &) {
+ }
+}
+
+void GlobalContext::notify_set_parameter(const ObjectHandle *o,
+ std::string const &name,
+ Variant const &value) {
+ cb_set_parameter(object_id(o), name, pack(value));
+}
+
+void GlobalContext::call_method(ObjectId id, std::string const &name,
+ PackedMap const &arguments) {
+ try {
+ m_local_objects.at(id)->call_method(name,
+ unpack(arguments, m_local_objects));
+ } catch (Exception const &) {
+ }
+}
+
+void GlobalContext::notify_call_method(const ObjectHandle *o,
+ std::string const &name,
+ VariantMap const &arguments) {
+ cb_call_method(object_id(o), name, pack(arguments));
+}
+
+std::shared_ptr
+GlobalContext::make_shared(std::string const &name,
+ const VariantMap ¶meters) {
+ std::unique_ptr sp = m_node_local_context->factory().make(name);
+ set_context(sp.get());
+
+ auto const id = object_id(sp.get());
+ remote_make_handle(id, name, parameters);
+
+ sp->construct(parameters);
+
+ return {sp.release(),
+ /* Custom deleter, we keep the corresponding global context,
+ * as well as the original deleter for the object. */
+ [global_context = this, deleter = sp.get_deleter()](ObjectHandle *o) {
+ /* Tell the other nodes before invoking the destructor, this is
+ * required
+ * to have synchronous destructors, which is needed by some client
+ * code. */
+ global_context->cb_delete_handle(object_id(o));
+
+ /* Locally destroy the object. */
+ deleter(o);
+ }};
+}
+
+boost::string_ref GlobalContext::name(const ObjectHandle *o) const {
+ assert(o);
+
+ return m_node_local_context->factory().type_name(*o);
+}
+} // namespace ScriptInterface
\ No newline at end of file
diff --git a/src/script_interface/GlobalContext.hpp b/src/script_interface/GlobalContext.hpp
new file mode 100644
index 00000000000..4d7596858e1
--- /dev/null
+++ b/src/script_interface/GlobalContext.hpp
@@ -0,0 +1,157 @@
+/*
+ * Copyright (C) 2020 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 .
+ */
+#ifndef ESPRESSO_SCRIPT_INTERFACE_OBJECTMANAGER_HPP
+#define ESPRESSO_SCRIPT_INTERFACE_OBJECTMANAGER_HPP
+
+/** @file
+ *
+ * Infrastructure to synchronize objects created on the head node
+ * with their corresponding remote copies.
+ *
+ * Implementation in @ref GlobalContext.cpp.
+ */
+
+#include "Context.hpp"
+#include "LocalContext.hpp"
+#include "MpiCallbacks.hpp"
+#include "ObjectHandle.hpp"
+#include "packed_variant.hpp"
+
+#include
+
+#include
+
+namespace ScriptInterface {
+
+/**
+ * @brief Global synchronizing context.
+ *
+ * Objects created in this context are synchronized
+ * between multiple MPI ranks. That is, for each instance
+ * created on the head node, a copy is created on all
+ * other ranks. If the original copy is mutated via
+ * @ref notify_set_parameter(), this change is also applied to all the
+ * copies. Calls to @ref notify_call_method() are also propagated to all
+ * ranks. The lifetime of the copies is tied to the original,
+ * if the original copy is destroyed on the head node,
+ * the remote copies are also destroyed.
+ */
+class GlobalContext : public Context {
+ using ObjectId = std::size_t;
+
+private:
+ /* Instances on this node that are managed by the
+ * head node. */
+ std::unordered_map m_local_objects;
+
+ std::shared_ptr m_node_local_context;
+
+private:
+ Communication::CallbackHandle
+ cb_make_handle;
+ Communication::CallbackHandle
+ cb_set_parameter;
+ Communication::CallbackHandle
+ cb_call_method;
+ Communication::CallbackHandle cb_delete_handle;
+
+public:
+ GlobalContext(Communication::MpiCallbacks &callbacks,
+ std::shared_ptr node_local_context)
+ : m_node_local_context(std::move(node_local_context)),
+ cb_make_handle(&callbacks,
+ [this](ObjectId id, const std::string &name,
+ const PackedMap ¶meters) {
+ make_handle(id, name, parameters);
+ }),
+ cb_set_parameter(&callbacks,
+ [this](ObjectId id, std::string const &name,
+ PackedVariant const &value) {
+ set_parameter(id, name, value);
+ }),
+ cb_call_method(&callbacks,
+ [this](ObjectId id, std::string const &name,
+ PackedMap const &arguments) {
+ call_method(id, name, arguments);
+ }),
+ cb_delete_handle(&callbacks,
+ [this](ObjectId id) { delete_handle(id); }) {}
+
+private:
+ /**
+ * @brief Callback for @c cb_make_handle
+ */
+ void make_handle(ObjectId id, const std::string &name,
+ const PackedMap ¶meters);
+ /**
+ * @brief Create remote instances
+ *
+ * @param id Internal identifier of the instance
+ * @param name Class name
+ * @param parameters Constructor parameters.
+ */
+ void remote_make_handle(ObjectId id, const std::string &name,
+ const VariantMap ¶meters);
+
+private:
+ /**
+ * @brief Callback for @c cb_set_parameter
+ */
+ void set_parameter(ObjectId id, std::string const &name,
+ PackedVariant const &value);
+
+public:
+ void notify_set_parameter(const ObjectHandle *o, std::string const &name,
+ Variant const &value) override;
+
+private:
+ /**
+ * @brief Callback for @c cb_call_method
+ */
+ void call_method(ObjectId id, std::string const &name,
+ PackedMap const &arguments);
+
+public:
+ void notify_call_method(const ObjectHandle *o, std::string const &name,
+ VariantMap const &arguments) override;
+
+private:
+ /**
+ * @brief Callback for @c cb_delete_handle
+ */
+ void delete_handle(ObjectId id) { m_local_objects.erase(id); }
+
+public:
+ /**
+ * @brief Get a new reference counted instance of a script interface
+ * object by name.
+ *
+ * Remote objects are automatically constructed.
+ */
+ std::shared_ptr
+ make_shared(std::string const &name, const VariantMap ¶meters) override;
+
+ boost::string_ref name(const ObjectHandle *o) const override;
+};
+} // namespace ScriptInterface
+
+#endif
diff --git a/src/script_interface/LocalContext.hpp b/src/script_interface/LocalContext.hpp
new file mode 100644
index 00000000000..8ea9a000d88
--- /dev/null
+++ b/src/script_interface/LocalContext.hpp
@@ -0,0 +1,67 @@
+/*
+ * Copyright (C) 2020 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 .
+ */
+#ifndef SCRIPT_INTERFACE_LOCAL_CONTEXT_HPP
+#define SCRIPT_INTERFACE_LOCAL_CONTEXT_HPP
+
+#include "Context.hpp"
+#include "ObjectHandle.hpp"
+
+#include
+
+namespace ScriptInterface {
+
+/**
+ * @brief Trivial context.
+ *
+ * This context just maintains a local copy of an
+ * object.
+ */
+class LocalContext : public Context {
+ Utils::Factory m_factory;
+
+public:
+ explicit LocalContext(Utils::Factory factory)
+ : m_factory(std::move(factory)) {}
+
+ const Utils::Factory &factory() const { return m_factory; }
+
+ void notify_call_method(const ObjectHandle *, std::string const &,
+ VariantMap const &) override {}
+ void notify_set_parameter(const ObjectHandle *, std::string const &,
+ Variant const &) override {}
+
+ std::shared_ptr
+ make_shared(std::string const &name, const VariantMap ¶meters) override {
+ auto sp = m_factory.make(name);
+ set_context(sp.get());
+
+ sp->construct(parameters);
+
+ return sp;
+ }
+
+ boost::string_ref name(const ObjectHandle *o) const override {
+ assert(o);
+
+ return factory().type_name(*o);
+ }
+};
+} // namespace ScriptInterface
+
+#endif
diff --git a/src/script_interface/None.hpp b/src/script_interface/None.hpp
index eec09e955c5..9aef9750724 100644
--- a/src/script_interface/None.hpp
+++ b/src/script_interface/None.hpp
@@ -22,6 +22,8 @@
#include
+#include
+
namespace ScriptInterface {
/**
diff --git a/src/script_interface/ObjectHandle.cpp b/src/script_interface/ObjectHandle.cpp
new file mode 100644
index 00000000000..48279aef6d9
--- /dev/null
+++ b/src/script_interface/ObjectHandle.cpp
@@ -0,0 +1,98 @@
+/*
+ * Copyright (C) 2010-2020 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 .
+ */
+
+#include "ObjectHandle.hpp"
+#include "Context.hpp"
+#include "ObjectState.hpp"
+#include "packed_variant.hpp"
+
+#include
+
+namespace ScriptInterface {
+void ObjectHandle::set_parameter(const std::string &name,
+ const Variant &value) {
+ if (m_context)
+ m_context->notify_set_parameter(this, name, value);
+
+ this->do_set_parameter(name, value);
+}
+
+Variant ObjectHandle::call_method(const std::string &name,
+ const VariantMap ¶ms) {
+ if (m_context)
+ m_context->notify_call_method(this, name, params);
+
+ return this->do_call_method(name, params);
+}
+
+std::string ObjectHandle::serialize() const {
+ ObjectState state;
+
+ auto const params = this->get_parameters();
+ state.params.resize(params.size());
+
+ PackVisitor v;
+
+ /* Pack parameters and keep track of ObjectRef parameters */
+ boost::transform(params, state.params.begin(),
+ [&v](auto const &kv) -> PackedMap::value_type {
+ return {kv.first, boost::apply_visitor(v, kv.second)};
+ });
+
+ /* Packed Object parameters */
+ state.objects.resize(v.objects().size());
+ boost::transform(v.objects(), state.objects.begin(), [](auto const &kv) {
+ return std::make_pair(kv.first, kv.second->serialize());
+ });
+
+ state.name = name().to_string();
+ state.internal_state = get_internal_state();
+
+ return Utils::pack(state);
+}
+
+ObjectRef ObjectHandle::deserialize(const std::string &packed_state,
+ Context &ctx) {
+ auto const state = Utils::unpack(packed_state);
+
+ std::unordered_map objects;
+ boost::transform(state.objects, std::inserter(objects, objects.end()),
+ [&ctx](auto const &kv) {
+ return std::make_pair(kv.first,
+ deserialize(kv.second, ctx));
+ });
+
+ VariantMap params;
+ for (auto const &kv : state.params) {
+ params[kv.first] = boost::apply_visitor(UnpackVisitor(objects), kv.second);
+ }
+
+ auto o = ctx.make_shared(state.name, params);
+ o->set_internal_state(state.internal_state);
+
+ return o;
+}
+
+boost::string_ref ObjectHandle::name() const {
+ return context() ? context()->name(this) : boost::string_ref{};
+}
+
+} /* namespace ScriptInterface */
diff --git a/src/script_interface/ObjectHandle.hpp b/src/script_interface/ObjectHandle.hpp
new file mode 100644
index 00000000000..3aa524abba2
--- /dev/null
+++ b/src/script_interface/ObjectHandle.hpp
@@ -0,0 +1,154 @@
+/*
+ * Copyright (C) 2015-2020 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 .
+ */
+
+#ifndef SCRIPT_INTERFACE_SCRIPT_INTERFACE_BASE_HPP
+#define SCRIPT_INTERFACE_SCRIPT_INTERFACE_BASE_HPP
+#include "Variant.hpp"
+
+#include
+
+#include
+
+namespace ScriptInterface {
+class Context;
+
+/**
+ * @brief Base class for interface handles.
+ */
+class ObjectHandle {
+public:
+ ObjectHandle() = default;
+
+ /* Copy has unclear semantics, so it should not be allowed. */
+ ObjectHandle(ObjectHandle const &) = delete;
+ ObjectHandle &operator=(ObjectHandle const &) = delete;
+ virtual ~ObjectHandle() = default;
+
+private:
+ friend class Context;
+ std::shared_ptr m_context = {};
+
+public:
+ boost::string_ref name() const;
+
+public:
+ /**
+ * @brief Responsible context.
+ */
+ Context *context() const { return m_context.get(); }
+
+public:
+ /**
+ * @brief Construct the handled object.
+ *
+ * This function is called on object creation with user
+ * provided parameters. This can be used if the SO has required parameters,
+ * it represents some type that can not reasonably be default constructed,
+ * or if the core implementation has to be chosen by a parameter. It is
+ * guaranteed that no getter or setter functions from this interface is
+ * called before construction (only @ref name() and @ref valid_parameters()),
+ * and it is only called once.
+ *
+ * The default implementation just calls @ref set_parameter for every
+ * parameter.
+ *
+ * @param params The parameters to the constructor. Only parameters that
+ * are valid for a default-constructed object are valid.
+ */
+ void construct(VariantMap const ¶ms) { do_construct(params); }
+
+private:
+ virtual void do_construct(VariantMap const ¶ms) {
+ for (auto const &p : params) {
+ do_set_parameter(p.first, p.second);
+ }
+ }
+
+public:
+ /**
+ * @brief Get current parameters.
+ * @return Parameters set in class.
+ */
+ VariantMap get_parameters() const {
+ VariantMap values;
+
+ for (auto const &p : valid_parameters()) {
+ values[p.data()] = get_parameter(p.data());
+ }
+
+ return values;
+ }
+
+ /**
+ * @brief Get required and optional parameters for class.
+ * @return Expected parameters.
+ */
+ virtual Utils::Span valid_parameters() const {
+ return {};
+ }
+
+ /**
+ * @brief Get single parameter.
+ *
+ * @param name Name of the parameter
+ * @return Value of parameter @p name
+ */
+ virtual Variant get_parameter(const std::string &name) const { return {}; }
+
+ /**
+ * @brief Set single parameter.
+ */
+ void set_parameter(const std::string &name, const Variant &value);
+
+private:
+ /**
+ * @brief Local implementation of @ref set_parameter.
+ */
+ virtual void do_set_parameter(const std::string &, const Variant &) {}
+
+public:
+ /**
+ * @brief Call a method on the object.
+ */
+ Variant call_method(const std::string &name, const VariantMap ¶ms);
+
+private:
+ /**
+ * @brief Local implementation of @c do_call_method.
+ *
+ * If not overridden by the implementation, this does nothing.
+ */
+ virtual Variant do_call_method(const std::string &, const VariantMap &) {
+ return none;
+ }
+
+public:
+ std::string serialize() const;
+
+ /**
+ * @brief Make object from serialized state.
+ */
+ static ObjectRef deserialize(const std::string &state, Context &ctx);
+
+private:
+ virtual std::string get_internal_state() const { return {}; }
+ virtual void set_internal_state(std::string const &state) {}
+};
+} /* namespace ScriptInterface */
+#endif
diff --git a/src/script_interface/ScriptObjectRegistry.hpp b/src/script_interface/ObjectList.hpp
similarity index 74%
rename from src/script_interface/ScriptObjectRegistry.hpp
rename to src/script_interface/ObjectList.hpp
index 197fb7fa3b7..059fcc9c541 100644
--- a/src/script_interface/ScriptObjectRegistry.hpp
+++ b/src/script_interface/ObjectList.hpp
@@ -1,5 +1,5 @@
/*
- * Copyright (C) 2010-2019 The ESPResSo project
+ * Copyright (C) 2010-2020 The ESPResSo project
* Copyright (C) 2002,2003,2004,2005,2006,2007,2008,2009,2010
* Max-Planck-Institute for Polymer Research, Theory Group
*
@@ -23,16 +23,23 @@
#define SCRIPT_INTERFACE_REGISTRY_HPP
#include "script_interface/ScriptInterface.hpp"
+#include "script_interface/get_value.hpp"
namespace ScriptInterface {
-
-template
-class ScriptObjectRegistry : public ScriptInterfaceBase {
+/**
+ * @brief Owning list of ObjectHandles
+ * @tparam ManagedType Type of the managed objects, needs to be
+ * derived from ObjectHandle
+ */
+template <
+ typename ManagedType,
+ class = std::enable_if_t::value>>
+class ObjectList : public ObjectHandle {
public:
- virtual void add_in_core(std::shared_ptr obj_ptr) = 0;
- virtual void remove_in_core(std::shared_ptr obj_ptr) = 0;
- Variant call_method(std::string const &method,
- VariantMap const ¶meters) override {
+ virtual void add_in_core(const std::shared_ptr &obj_ptr) = 0;
+ virtual void remove_in_core(const std::shared_ptr &obj_ptr) = 0;
+ Variant do_call_method(std::string const &method,
+ VariantMap const ¶meters) override {
if (method == "add") {
auto obj_ptr =
@@ -57,7 +64,7 @@ class ScriptObjectRegistry : public ScriptInterfaceBase {
ret.reserve(m_elements.size());
for (auto const &e : m_elements)
- ret.emplace_back(e->id());
+ ret.emplace_back(e);
return ret;
}
diff --git a/src/script_interface/ObjectState.hpp b/src/script_interface/ObjectState.hpp
new file mode 100644
index 00000000000..fe4bb13a0d5
--- /dev/null
+++ b/src/script_interface/ObjectState.hpp
@@ -0,0 +1,52 @@
+/*
+ * Copyright (C) 2020 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 .
+ */
+#ifndef ESPRESSO_SCRIPT_INTERFACE_OBJECTSTATE_HPP
+#define ESPRESSO_SCRIPT_INTERFACE_OBJECTSTATE_HPP
+
+#include "packed_variant.hpp"
+
+#include
+#include
+#include
+#include
+
+#include
+#include
+#include
+
+namespace ScriptInterface {
+/**
+ * @brief State of an object ready for serialization.
+ *
+ * This specifies the internal serialization format and
+ * should not be used outside of the class.
+ */
+struct ObjectState {
+ std::string name;
+ PackedMap params;
+ std::vector> objects;
+ std::string internal_state;
+
+ template void serialize(Archive &ar, long int) {
+ ar &name ¶ms &objects &internal_state;
+ }
+};
+} // namespace ScriptInterface
+
+#endif
diff --git a/src/script_interface/ParallelScriptInterface.cpp b/src/script_interface/ParallelScriptInterface.cpp
deleted file mode 100644
index 6540c245f30..00000000000
--- a/src/script_interface/ParallelScriptInterface.cpp
+++ /dev/null
@@ -1,234 +0,0 @@
-/*
- * Copyright (C) 2015-2019 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 .
- */
-
-#include "ParallelScriptInterface.hpp"
-
-#include
-#include
-
-#include "ParallelScriptInterfaceSlave.hpp"
-
-#include
-
-namespace {
-Communication::MpiCallbacks *m_cb = nullptr;
-
-void make_remote_handle() {
- assert(m_cb && "Not initialized.");
-
- new ScriptInterface::ParallelScriptInterfaceSlave(m_cb);
-}
-} // namespace
-
-REGISTER_CALLBACK(make_remote_handle)
-
-namespace ScriptInterface {
-ParallelScriptInterface::ParallelScriptInterface(std::string const &name)
- : m_callback_id(m_cb, [](CallbackAction) {}) {
- assert(m_cb && "Not initialized!");
-
- /* Create the slaves */
- m_cb->call(make_remote_handle);
-
- call(CallbackAction::NEW);
-
- /* Create local object */
- m_p = ScriptInterfaceBase::make_shared(
- name, ScriptInterfaceBase::CreationPolicy::LOCAL);
-
- /* Bcast class name and global id to the slaves */
- std::pair what = std::make_pair(m_p->id(), name);
- boost::mpi::broadcast(m_cb->comm(), what, 0);
-}
-
-ParallelScriptInterface::~ParallelScriptInterface() {
- /* Delete the instances on the other nodes */
- try {
- call(CallbackAction::DELETE);
- } catch (...) {
- /* We have to do MPI calls in the destructor, which
- can cause exceptions. To avoid sanitizer warnings
- about this, and to clarify we directly explicitly
- call terminate, which would happen anyway.
- */
- std::terminate();
- }
-}
-
-bool ParallelScriptInterface::operator==(ParallelScriptInterface const &rhs) {
- return this->get_underlying_object() == rhs.get_underlying_object();
-}
-
-bool ParallelScriptInterface::operator!=(ParallelScriptInterface const &rhs) {
- return !(*this == rhs);
-}
-
-void ParallelScriptInterface::initialize(Communication::MpiCallbacks &cb) {
- m_cb = &cb;
-}
-
-void ParallelScriptInterface::construct(VariantMap const ¶ms) {
- call(CallbackAction::CONSTRUCT);
-
- auto p = unwrap_variant_map(params);
- boost::mpi::broadcast(m_cb->comm(), p, 0);
-
- m_p->construct(p);
-}
-
-void ParallelScriptInterface::set_parameter(const std::string &name,
- const Variant &value) {
- std::pair d(name, Variant());
-
- if (is_type(value)) {
- d.second = map_parallel_to_local_id(value);
- } else {
- d.second = value;
- }
-
- call(CallbackAction::SET_PARAMETER);
-
- boost::mpi::broadcast(m_cb->comm(), d, 0);
-
- m_p->set_parameter(d.first, d.second);
-
- collect_garbage();
-}
-
-Variant ParallelScriptInterface::call_method(const std::string &name,
- const VariantMap ¶meters) {
- call(CallbackAction::CALL_METHOD);
- VariantMap p = unwrap_variant_map(parameters);
-
- auto d = std::make_pair(name, p);
- /* Broadcast method name and parameters */
- boost::mpi::broadcast(m_cb->comm(), d, 0);
-
- auto ret = map_local_to_parallel_id(m_p->call_method(name, p));
-
- collect_garbage();
-
- return ret;
-}
-
-Variant ParallelScriptInterface::get_parameter(std::string const &name) const {
- auto p = m_p->get_parameter(name);
-
- return map_local_to_parallel_id(p);
-}
-
-VariantMap ParallelScriptInterface::get_parameters() const {
- auto p = m_p->get_parameters();
-
- /* Wrap the object ids */
- for (auto &it : p) {
- it.second = map_local_to_parallel_id(it.second);
- }
-
- return p;
-}
-
-VariantMap ParallelScriptInterface::unwrap_variant_map(VariantMap const &map) {
- /* Copy parameters into a non-const buffer, needed by boost::mpi */
- auto p = map;
-
- /* Unwrap the object ids */
- for (auto &it : p) {
- if (is_type(it.second)) {
- it.second = map_parallel_to_local_id(it.second);
- }
- }
-
- return p;
-}
-
-Variant
-ParallelScriptInterface::map_local_to_parallel_id(Variant const &value) const {
- if (is_type(value)) {
- /** Check if the objectid is the empty object (ObjectId()),
- * if so it does not need translation, the empty object
- * has the same id everywhere.
- */
- auto oid = boost::get(value);
-
- if (oid != ObjectId()) {
- return obj_map.at(oid)->id();
- }
- return oid;
- }
- if (is_type>(value)) {
- auto const &in_vec = boost::get>(value);
- std::vector out_vec;
- out_vec.reserve(in_vec.size());
-
- for (auto const &e : in_vec) {
- out_vec.emplace_back(map_local_to_parallel_id(e));
- }
-
- return out_vec;
- }
- return value;
-}
-
-Variant
-ParallelScriptInterface::map_parallel_to_local_id(Variant const &value) {
- const auto outer_id = boost::get(value);
-
- auto so_ptr = get_instance(outer_id).lock();
-
- auto po_ptr = std::dynamic_pointer_cast(so_ptr);
-
- if (po_ptr != nullptr) {
- auto inner_id = po_ptr->get_underlying_object()->id();
-
- /* Store a pointer to the object */
- obj_map[inner_id] = po_ptr;
-
- /* and return the id of the underlying object */
- return inner_id;
- }
- if (so_ptr == nullptr) {
- /* Release the object */
- obj_map.erase(outer_id);
-
- /* Return None */
- return ObjectId();
- }
- throw std::runtime_error(
- "Parameters passed to Parallel entities must also be parallel.");
-}
-
-void ParallelScriptInterface::collect_garbage() {
- /* Removal condition, the instance is removed iff
- its payload is not used anywhere. In this case
- the reference count is one, because the host object
- still holds a pointer.
- */
- auto pred = [](map_t::value_type const &e) -> bool {
- return e.second->get_underlying_object().use_count() == 1;
- };
-
- for (auto it = obj_map.begin(); it != obj_map.end();) {
- if (pred(*it)) {
- obj_map.erase(it);
- }
- ++it;
- }
-}
-} /* namespace ScriptInterface */
diff --git a/src/script_interface/ParallelScriptInterface.hpp b/src/script_interface/ParallelScriptInterface.hpp
deleted file mode 100644
index 911848abda1..00000000000
--- a/src/script_interface/ParallelScriptInterface.hpp
+++ /dev/null
@@ -1,97 +0,0 @@
-/*
- * Copyright (C) 2010-2019 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 .
- */
-
-#ifndef SCRIPT_INTERFACE_PARALLEL_SCRIPT_INTERFACE_HPP
-#define SCRIPT_INTERFACE_PARALLEL_SCRIPT_INTERFACE_HPP
-
-#include
-
-#include "MpiCallbacks.hpp"
-#include "script_interface/ScriptInterface.hpp"
-
-namespace ScriptInterface {
-
-class ParallelScriptInterface : public ScriptInterfaceBase {
-public:
- enum class CallbackAction {
- NEW,
- CONSTRUCT,
- SET_PARAMETER,
- CALL_METHOD,
- DELETE
- };
-
- explicit ParallelScriptInterface(std::string const &name);
- ~ParallelScriptInterface() override;
-
- /**
- * @brief Initialize the mpi callback for instance creation.
- */
- static void initialize(Communication::MpiCallbacks &cb);
-
- bool operator==(ParallelScriptInterface const &rhs);
- bool operator!=(ParallelScriptInterface const &rhs);
-
- /**
- * @brief Get the payload object.
- */
- std::shared_ptr get_underlying_object() const {
- return std::static_pointer_cast(m_p);
- }
-
- void construct(VariantMap const ¶ms) override;
- const std::string name() const { return m_p->name(); }
- void set_parameter(const std::string &name, const Variant &value) override;
- Utils::Span valid_parameters() const override {
- return m_p->valid_parameters();
- }
-
- Variant get_parameter(const std::string &name) const override;
- VariantMap get_parameters() const override;
- Variant call_method(const std::string &name,
- const VariantMap ¶meters) override;
-
- /* Id mapping */
- Variant map_local_to_parallel_id(Variant const &value) const;
- Variant map_parallel_to_local_id(Variant const &value);
-
-private:
- using map_t = std::map>;
-
- VariantMap unwrap_variant_map(VariantMap const &map);
-
- void call(CallbackAction action) { m_callback_id(action); }
-
- /**
- * @brief Remove instances that are not used by anybody but us.
- */
- void collect_garbage();
-
- /* Data members */
- Communication::CallbackHandle m_callback_id;
- /* Payload object */
- std::shared_ptr m_p;
- map_t obj_map;
-};
-
-} /* namespace ScriptInterface */
-
-#endif
diff --git a/src/script_interface/ParallelScriptInterfaceSlave.hpp b/src/script_interface/ParallelScriptInterfaceSlave.hpp
deleted file mode 100644
index 438d1449f3b..00000000000
--- a/src/script_interface/ParallelScriptInterfaceSlave.hpp
+++ /dev/null
@@ -1,133 +0,0 @@
-/*
- * Copyright (C) 2010-2019 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 .
- */
-
-#ifndef SCRIPT_INTERFACE_PARALLEL_SCRIPT_INTERFACE_SLAVE_HPP
-#define SCRIPT_INTERFACE_PARALLEL_SCRIPT_INTERFACE_SLAVE_HPP
-
-#include "MpiCallbacks.hpp"
-#include "ScriptInterfaceBase.hpp"
-
-#include "ParallelScriptInterface.hpp"
-
-#include
-#include
-
-namespace ScriptInterface {
-class ParallelScriptInterfaceSlave {
-private:
- const auto &comm() const { return m_callback_id.cb()->comm(); }
-
- using CallbackAction = ParallelScriptInterface::CallbackAction;
-
-public:
- ParallelScriptInterfaceSlave(Communication::MpiCallbacks *cb)
- : m_callback_id(cb, [this](CallbackAction a) { mpi_slave(a); }) {}
-
- std::shared_ptr m_p;
-
- static std::map &get_translation_table() {
- static std::map m_translation_table;
-
- return m_translation_table;
- }
-
- /* If the variant encapsulates an object id we translate the
- master id to a local one */
- static void translate_id(Variant &v) {
- if (is_type(v)) {
- v = get_translation_table().at(boost::get(v));
- }
- }
-
- static void translate_id(VariantMap &map) {
- for (auto &e : map) {
- translate_id(e.second);
- }
- }
-
- VariantMap bcast_variant_map() const {
- VariantMap ret;
- boost::mpi::broadcast(comm(), ret, 0);
-
- /* If the parameter is an object we have to translate it first to a
- local id.
- */
- translate_id(ret);
-
- return ret;
- }
-
-private:
- Communication::CallbackHandle