From e8991f54ced96f3a576c411993c307ef1d79642b Mon Sep 17 00:00:00 2001
From: HEDIN Florent <florent.hedin@inria.fr>
Date: Fri, 8 Mar 2019 12:08:48 +0100
Subject: [PATCH] CI enabled, reproducibility added, ready for release v1.2
---
.gitignore | 3 +
.gitlab-ci.yml | 35 ++
CHANGELOG.md | 42 +-
CMakeLists.txt | 16 +-
LICENSE.txt | 2 +-
README.md | 164 +++---
...cies.sh => download_deps_and_run_cmake.sh} | 2 +-
ci/.gitignore | 3 +
ci/README.md | 3 +
ci/input_generalized_parRep.lua | 495 ++++++++++++++++++
ci/lj7/Integrator.xml | 2 +
ci/lj7/State.xml | 26 +
ci/lj7/System.xml | 36 ++
ci/lj7/omm_LJ7.py | 139 +++++
ci/ref.sim/README.md | 23 +
ci/ref.sim/lj7.db | Bin 0 -> 16384 bytes
ci/ref.sim/seeds.bin | Bin 0 -> 16384 bytes
ci/run.ci.sh | 27 +
include/dyna/GelmanRubin.hpp | 13 +-
include/dyna/observable.hpp | 2 +-
include/dyna/parRep.hpp | 2 +-
include/dyna/parRep_FV.hpp | 2 +-
include/dyna/parRep_FV_multiExits.hpp | 2 +-
include/dyna/parRep_abstract.hpp | 2 +-
include/dyna/runSim.hpp | 11 +-
include/global.hpp | 2 +-
include/logger.hpp | 2 +-
include/md/md_interface.hpp | 2 +-
include/md/omm_interface.hpp | 2 +-
include/rand.hpp | 38 +-
include/utils/lua_interface.hpp | 2 +-
include/utils/mpi_utils.hpp | 2 +-
mol/1d7h_i/input_generalized_parRep.lua | 15 +-
mol/ala2vac/input_generalized_parRep.lua | 10 +-
.../input_generalized_parRep_multiExit.lua | 10 +-
mol/ala2vac/input_parRep.lua | 5 +-
.../input_generalized_parRep.lua | 22 +-
mol/ala2vac_transient/input_parRep.lua | 11 +-
parRep.doxy | 2 +-
run/run.bash.sh | 31 +-
run/run.slurm.sh | 8 +-
src/dyna/GelmanRubin.cpp | 26 +-
src/dyna/observable.cpp | 2 +-
src/dyna/parRep.cpp | 2 +-
src/dyna/parRep_FV.cpp | 3 +-
src/dyna/parRep_FV_multiExits.cpp | 2 +-
src/dyna/runSim.cpp | 55 +-
src/logger.cpp | 2 +-
src/main.cpp | 29 +-
src/md/md_interface.cpp | 2 +-
src/md/omm_interface.cpp | 30 +-
src/rand.cpp | 155 ++++--
src/utils/lua_interface.cpp | 90 +++-
src/utils/mpi_utils.cpp | 2 +-
54 files changed, 1364 insertions(+), 252 deletions(-)
create mode 100644 .gitlab-ci.yml
rename build/{download_dependencies.sh => download_deps_and_run_cmake.sh} (94%)
create mode 100644 ci/.gitignore
create mode 100644 ci/README.md
create mode 100755 ci/input_generalized_parRep.lua
create mode 100644 ci/lj7/Integrator.xml
create mode 100644 ci/lj7/State.xml
create mode 100644 ci/lj7/System.xml
create mode 100755 ci/lj7/omm_LJ7.py
create mode 100644 ci/ref.sim/README.md
create mode 100644 ci/ref.sim/lj7.db
create mode 100644 ci/ref.sim/seeds.bin
create mode 100755 ci/run.ci.sh
diff --git a/.gitignore b/.gitignore
index 2804c78..b51356f 100644
--- a/.gitignore
+++ b/.gitignore
@@ -9,3 +9,6 @@ build/LuaJIT-2.0.5*
build/Makefile
build/cmake_install.cmake
build/openmm-7.3.0*
+build/CTestTestfile.cmake
+build/parRep
+build/Testing
diff --git a/.gitlab-ci.yml b/.gitlab-ci.yml
new file mode 100644
index 0000000..8459704
--- /dev/null
+++ b/.gitlab-ci.yml
@@ -0,0 +1,35 @@
+stages:
+- build
+- test
+
+variables:
+ MAKE_COMMAND: "make -j 4"
+
+cache: # Directories that need to be kept between the build job and the test job
+ paths:
+ - build/
+
+#----------------------------------------------------------------
+# Linux build/test procedure
+#----------------------------------------------------------------
+
+# Build job
+parrep_build_ubuntu:
+ stage: build
+ tags:
+ - parrep-ubuntu-1604-amd64
+ script:
+ - cd build
+ - bash download_deps_and_run_cmake.sh
+ - make
+
+# Test job
+parrep_test_ubuntu:
+ stage: test
+ tags:
+ - parrep-ubuntu-1604-amd64
+ script:
+ - echo "$PWD"
+ - cd build
+ - ctest -V
+
diff --git a/CHANGELOG.md b/CHANGELOG.md
index dfb604c..0a7d665 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -1,7 +1,39 @@
+# v1.2.0
+
+Mar 08th, 2019
+
+The software was modified in order to allow reproducibility in some cases (the limiting factor is OpenMM which does not
+always provides deterministic output even when using the same seeds (!), see http://docs.openmm.org/latest/userguide/library.html#determinism ).
+
+The main executable now has 2 more command line options, '--inp-seeds [fname]' or '--out-seeds [fname]' for respectively loading seeds or writting seeds
+from/to a unique binary file. See rand.hpp/rand.cpp, od the Doxygen doc for more details.
+
+These modifications now allow Continuous Integration (CI) on the infrastructure provided by INRIA : in mol/ci a small test case will be executed at each commit to the repository and compared to reference results.
+
+The two other minor modifications concern the Lua scripts:
+
+* "get_minimised_energy_crdvels" was added to the set of functions that the user can call from the Lua script : it simply combines in one call what
+"get_minimised_energy" and "get_minimised_crdvels" already provided.
+
+* a extra optional parameter is available for the "simulation" parameters block when "simulation.algorithm" is "PARREP_FV" ; this parameter is
+"simulation.minAccumulatedObs" : it will enforce that at least minAccumulatedObs observations of an observable have already been accumulated
+before the convergence test is performed ; this may be useful if there is a risk of early pseudo-convergence for some of the observables when only a few samples have been accumulated.
+
+
+**Download sources:**
+
+https://gitlab.inria.fr/parallel-replica/gen.parRep/tags/v1.2.0
+
+or
+
+https://github.com/FHedin/gen.parRep/releases/tag/v1.2.0
+
----------------------------------------------
# v1.1.0
+Nov 16th, 2018
+
Finished the implementation of a "MD_interface" abstract class for future addition of different Molecular Dynamics engines.
This required important changes to a significant part of the code.
Optimizations to the MPI communications were also performed, therefore the minor version number
@@ -19,7 +51,9 @@ https://github.com/FHedin/gen.parRep/releases/tag/v1.1.0
# v1.0.1
-simply added directory structure to readme file
+Aug 16th, 2018
+
+Simply added directory structure to readme file.
**Download sources:**
@@ -33,7 +67,9 @@ https://github.com/FHedin/gen.parRep/releases/tag/v1.0.1
# v1.0.0
-Release of v1.0.0 : several bug fixes like memory leaks adressed, added input files for alanine
+Aug 13th, 2018
+
+Release of v1.0.0 : several bug fixes like memory leaks addressed, added input files for alanine
dipeptide using transient propagation.
**Download sources:**
@@ -52,7 +88,7 @@ July 6th, 2018
First public release of the software.
It is following the submission of the corresponding article on arXiv ; this is a release candidate 1.0.0-rc1,
-the 1.0.0 release should follow within the next weeks
+the 1.0.0 release should follow within the next weeks.
**Download sources:**
diff --git a/CMakeLists.txt b/CMakeLists.txt
index b27dc9a..91af3b1 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -15,6 +15,8 @@ if(NOT CMAKE_BUILD_TYPE)
FORCE)
endif()
+enable_testing()
+
set(TGT "parRep" CACHE STRING "The executable name")
project(${TGT} C CXX)
@@ -37,7 +39,7 @@ include
include/dyna
include/md
include/utils
-# external headers (dependancies provided with this project)
+# external headers (dependencies provided with this project)
external/sol2
)
@@ -61,12 +63,10 @@ src/md/omm_interface.cpp
# utilities files
src/utils/mpi_utils.cpp
src/utils/lua_interface.cpp
-# external dependancies provided with this project
+# external dependencies provided with this project
external/luasqlite3/lsqlite3.c
)
-add_executable(${TGT} ${SRCS})
-
set(cWarnings "-Wall -Wextra")
# more warnings for c++
# set(cxxWarnings "-Wall -Wextra -Wformat=2 -Wshadow -Wconversion -Wuseless-cast")
@@ -90,4 +90,12 @@ set(CMAKE_CXX_FLAGS_RELEASE "-O3 -march=native")
set(CMAKE_CXX_FLAGS_RELWITHDEBINFO "-O2 -g -march=native")
set(CMAKE_CXX_FLAGS_MINSIZEREL "-Os -s -march=native")
+add_executable(${TGT} ${SRCS})
+
target_link_libraries(${TGT} ${OpenMM_LIBRARY} ${LUA_LIBRARIES} ${MPI_LIBRARIES} ${SQLITE3_LIBRARIES})
+
+add_test(
+NAME test_ci_ar7_gen.parRep
+COMMAND ${CMAKE_CURRENT_SOURCE_DIR}/ci/run.ci.sh
+WORKING_DIRECTORY ${CMAKE_CURRENT_SOURCE_DIR}/ci
+)
diff --git a/LICENSE.txt b/LICENSE.txt
index 9e1f08e..3a0db8c 100644
--- a/LICENSE.txt
+++ b/LICENSE.txt
@@ -1,4 +1,4 @@
-Copyright (c) 2016-2018, Florent Hédin, Tony Lelièvre, and École des Ponts - ParisTech
+Copyright (c) 2016-2019, Florent Hédin, Tony Lelièvre, and École des Ponts - ParisTech,
All rights reserved.
Redistribution and use in source and binary forms, with or without modification,
diff --git a/README.md b/README.md
index 9e0f280..3e5253b 100644
--- a/README.md
+++ b/README.md
@@ -1,137 +1,162 @@
-----------------------------------------------
-# gen.parRep (v1.1)
-----------------------------------------------
+# gen.parRep (v1.2)
C++ implementation of the Generalized Parallel Replica algorithm
Link to this repository :
- * https://gitlab.inria.fr/parallel-replica/gen.parRep
+
+* https://gitlab.inria.fr/parallel-replica/gen.parRep
Link to the project gitlab page (includes more resources):
- * https://gitlab.inria.fr/parallel-replica
+
+* https://gitlab.inria.fr/parallel-replica
See the following articles:
- * A generalized parallel replica dynamics: Binder, Lelièvre & Simpson, 2015: https://doi.org/10.1016/j.jcp.2015.01.002
- * A new implementation of the Generalized Parallel Replica Dynamics targeting metastable biochemical systems: Hédin & Lelièvre, 2018: https://arxiv.org/abs/1807.02431
+* A generalized parallel replica dynamics: Binder, Lelièvre & Simpson, 2015: https://doi.org/10.1016/j.jcp.2015.01.002
+
+* gen.parRep: a first implementation of the Generalized Parallel Replica dynamics for the long time simulation of metastable biochemical systems: Hédin & Lelièvre, 2018: https://arxiv.org/abs/1807.02431 and https://doi.org/10.1016/j.cpc.2019.01.005
Molecular dynamics is performed by using external codes linked to this program such as:
- * OpenMM (https://github.com/pandegroup/openmm and https://simtk.org/home/openmm)
-----------------------------------------------
+* OpenMM (https://github.com/pandegroup/openmm and https://simtk.org/home/openmm)
+
+If you find any error, bug, or limitation, please open an issue on https://gitlab.inria.fr/parallel-replica/gen.parRep/issues so that the community can benefit from its correction ; also feel free to clone/fork the repo and perform the modifications yourself if you have the programming abilities to do so !
+
+Continuous Integration is performed on a test case system (see ./ci), ensuring that no release breaks the ability to compile the software; furthermore reference results from a previous run are available in ./ci/ref.sim in an attempt to achieve reproducibility whenever possible.
+
+[](https://gitlab.inria.fr/parallel-replica/gen.parRep/commits/master)
+
## DOCUMENTATION
-----------------------------------------------
Please consult the Gitlab wiki for an overview of this software use:
https://gitlab.inria.fr/parallel-replica/gen.parRep/wikis/home
-Doxygen can be used for generating a code reference documentation within the docs directory:
- * doxygen parRep.doxy
+Doxygen can be used for generating a code reference documentation within the docs directory (particularly useful if you plan to modify/extend the software):
+
+* doxygen parRep.doxy
By default an HTML reference is generated (ready to read), and a docs/latex directory is also created, although the pdf version requires to be manually compiled (pdfLaTeX, see the Makefile in ./docs/latex).
-The Lua input files, in ./mol, used for starting computations, are self documented, and together
-with the bash submission files in the ./run directory should be enough for starting to use the program
+The Lua input files, in ./mol, used for starting simulations, are self documented, and together
+with the bash submission files in the ./run directory, it should be enough for starting to use the program.
+
+## DEPENDENCIES
+
+### Dependencies required before running the CMakeLists.txt script
-----------------------------------------------
-## DEPENDANCIES
-----------------------------------------------
+You will need an MPI development framework installed, compatible with the MPI 3.0 or newer standard : tested implementations :
-### Dependancies required before running the CMakeLists.txt script
+* Open MPI 1.10.2 (from Ubuntu 16.04 repositories)
-You will need an MPI development framework installed, compatible with the MPI 3.0 or newer standard : tested implementations :
- * Open MPI 1.10.2 (from Ubuntu 16.04 repositories)
CMake will try to locate it automatically.
-You will need the SQLite3 dynamic library on your system:
-See :
- * https://sqlite.org/index.html
+You will need the SQLite3 dynamic library on your system. See :
-### Dependancies downloaded automatically via the script "build/download_dependencies.sh"
+* https://sqlite.org/index.html
-The script "build/download_dependencies.sh" can download and compile if required the OpenMM and LuaJIT dependancies.
+### Dependencies downloaded automatically via the script "build/download_deps_and_run_cmake.sh"
-You can however use your own version of the OpenMM library, or any Lua inplementation (either compiled yourself or downloaded somewhere).
+The script "build/download_deps_and_run_cmake.sh" will download and compile if required the OpenMM and LuaJIT Dependencies.
+
+You can however use your own version of the OpenMM library, or any Lua implementation (either compiled yourself or downloaded somewhere).
#### Ressources (OpenMM):
+
* see https://simtk.org/home/openmm
+
* and/or https://github.com/pandegroup/openmm
+
* tested with version 7.0 to 7.3.
-You may need Nvidia CUDA or AMD OpenCL toolkit for enabling GPU acceleration ; see OpenMM documentation.
+You may need NVIDIA CUDA or AMD OpenCL toolkit for enabling GPU acceleration ; see OpenMM documentation.
You may also need to edit CMakeLists.txt for specifying path to the include and lib directories of OpenMM.
CMake will try to locate it automatically.
#### Ressources (Lua/LuaJIT):
A Lua implementation compatible with Lua API version >= 5.1 is required :
- * A release of Lua 5.x is usually already installed by default for recent linux versions, try to execute 'lua' and/or 'locate liblua'
- * You can Download and compile the official implementation : http://www.lua.org/download.html
- * The LuaJIT implementation can provide an important speedup : http://luajit.org/download.html
+
+* A release of Lua 5.x is usually already installed by default for recent linux versions, try to execute 'lua' and/or 'locate liblua'
+
+* You can Download and compile the official implementation : http://www.lua.org/download.html
+
+* The LuaJIT implementation can provide an important speedup : http://luajit.org/download.html
+
Please download and compile any of the two. See below for hints for setting DCMAKE_PREFIX_PATH in case of a manual install.
CMake will try to locate it automatically.
-### Dependancies provided in directory "external"
+### Dependencies provided in directory "external"
-The excellent Sol2 (header-only Lua<->c++ inteface) and LuaSQLite3 (a Lua/SQLite3 binding) dependancies are provided
+The excellent Sol2 (header-only Lua<->c++ interface) and LuaSQLite3 (a Lua/SQLite3 binding) Dependencies are provided
in the './external' directory and are automatically included/compiled if required, and thus should not require extra
-configuration.
-See :
- * https://github.com/ThePhD/sol2
- * http://lua.sqlite.org/index.cgi/home
+configuration. See :
+
+* https://github.com/ThePhD/sol2
+
+* http://lua.sqlite.org/index.cgi/home
-----------------------------------------------
## COMPILE & INSTALL
-----------------------------------------------
C and C++ compiler compatible with the C99 and C++14 standards are required.
Tested compilers:
- * gcc/g++ 5.4.0 and 6.2.0 (from Ubuntu 16.04 repositories)
+
+* gcc/g++ 5.4.0 and 6.2.0 (from Ubuntu 16.04 repositories)
Be sure to have CMake installed (http://www.cmake.org/), available on most repositories.
Tested version :
- * CMake 3.5.1 (from Ubuntu 16.04 repositories)
-Create a build directory and move to that directory:
- * mkdir build && cd build
+* CMake 3.5.1 (from Ubuntu 16.04 repositories)
+
+Create a build directory and move to that directory:
+
+* mkdir build && cd build
For building a debug or release or an intermediate release with debug information, do:
- * cmake -DCMAKE_BUILD_TYPE=Debug ..
- * cmake -DCMAKE_BUILD_TYPE=RelWithDebInfo ..
- * cmake -DCMAKE_BUILD_TYPE=Release .. (default, with cpu targeting for the current compilation machine)
-Debug builds are possibly slower, usually more memory consuming, but useful when debugging with gdb or valgrind.
+* cmake -DCMAKE_BUILD_TYPE=Debug ..
+
+* cmake -DCMAKE_BUILD_TYPE=RelWithDebInfo ..
-If some dependancies were note detected (e.g. because they are not available within /usr/local)
+* cmake -DCMAKE_BUILD_TYPE=Release .. (default, with cpu targeting for the current compilation machine)
+
+Debug builds are possibly slower, usually more memory consuming, but useful when debugging with gdb or Valgrind.
+
+If some dependencies were note detected (e.g. because they are not available within /usr/local)
it is required to add the path where the were installed to CMAKE_PREFIX_PATH, e.g.:
- * cmake -DCMAKE_BUILD_TYPE=Release -DCMAKE_PREFIX_PATH="$HOME/bin/something" ..
+
+* cmake -DCMAKE_BUILD_TYPE=Release -DCMAKE_PREFIX_PATH="$HOME/bin/something" ..
Then, once CMake built a Makefile without error, just execute the following in order to build the executable:
- * make
-For a verbose make, use:
- * make VERBOSE=1
+* make
+
+For a verbose make, use:
+
+* make VERBOSE=1
For a faster multi-core (on N cores) build, use:
- * make -j N
-Please never edit the autogenerated Makefile, edit the CMakeLists.txt instead.
+* make -j N
+
+Please never edit the automatically generated Makefile, edit the CMakeLists.txt instead.
+
+For specifying another compiler on linux, for example clang, or a proprietary one like Intel icc (untested) :
+
+* CC=clang CXX=clang++ cmake ..
+
+* CC=icc CXX=icpc cmake ..
-For specifying another compiler on linux, for example clang, or a proprietary one like Intel icc (untested) :
- * CC=clang CXX=clang++ cmake ..
- * CC=icc CXX=icpc cmake ..
-
In any case, you may need to edit the variables CMAKE_C_FLAGS_* and CMAKE_CXX_FLAGS_*
for setting proper levels of optimization. By default
-When using OpenMPI and forcing another C/C++ compilers version, the following might be required (possibly unsafe !):
- * OMPI_CC=gcc-5 OMPI_CXX=g++-5 CC=gcc-5 CXX=g++-5 ccmake ..
+When using OpenMPI and forcing another C/C++ compilers version, the following might be required (possibly unsafe !):
+
+* OMPI_CC=gcc-5 OMPI_CXX=g++-5 CC=gcc-5 CXX=g++-5 ccmake ..
-----------------------------------------------
## LICENSING (all files excepted subdirectory external and its content)
-----------------------------------------------
-Copyright (c) 2016-2018, Florent Hédin, Tony Lelièvre, and École des Ponts - ParisTech
+
+Copyright (c) 2016-2019, Florent Hédin, Tony Lelièvre, and École des Ponts - ParisTech
All rights reserved.
The 3-clause BSD license is applied to this software.
@@ -140,9 +165,7 @@ See LICENSE.txt
See in ./external/* for licensing information concerning the embedded libraries provided in './external'
-----------------------------------------------
## NOTES CONCERNING OpenMM
-----------------------------------------------
The software should automatically detect the fastest OpenMM Platform available on your computer (i.e. CUDA, OpenCL, ...)
If not it will run with the slow Reference platform : it is most probably because the OpenMM directory with the
@@ -151,18 +174,18 @@ If not it will run with the slow Reference platform : it is most probably becaus
You may need to export the following OPENMM_PLUGIN_DIR environment variable to solve the problem :
For example for a custom installation in /home/$USER/bin/openmm
- * export OPENMM_PLUGIN_DIR=$HOME/bin/openmm/lib/plugins
+
+* export OPENMM_PLUGIN_DIR=$HOME/bin/openmm/lib/plugins
For setting the amount of cpu threads used by each MPI process use the following env. variable :
- * export OPENMM_CPU_THREADS=1
-----------------------------------------------
+* export OPENMM_CPU_THREADS=1
+
## Directory structure
-----------------------------------------------
**src** and **include** -> contains the C++ source and headers of the gen.parRep software.
-**build** -> directory where dependancies will be downloaded, and the program built.
+**build** -> directory where dependencies will be downloaded, and the program built.
**external** -> contains source and headers of dependencies used by the gen.parRep software.
@@ -172,5 +195,6 @@ For setting the amount of cpu threads used by each MPI process use the following
**mol** -> contains ready to use test systems, i.e. input files and OpenMM configurations of various molecular systems.
-**run** -> contains bash scripts demonstrating how to run the software using either mpirun or the SLURM scheduler; also contains three bash files for runing the software under the following debuggers/profilers: GDB, Valgrind and Scalasca.
+**run** -> contains bash scripts demonstrating how to run the software using either mpirun or the SLURM scheduler; also contains bash files for running the software under the following debuggers/profilers: GDB, Valgrind and Scalasca.
+**ci** -> contains a molecular system, an input script and reference data for performing continuous integration and non regression tests.
diff --git a/build/download_dependencies.sh b/build/download_deps_and_run_cmake.sh
similarity index 94%
rename from build/download_dependencies.sh
rename to build/download_deps_and_run_cmake.sh
index 976dad2..f27d3c5 100755
--- a/build/download_dependencies.sh
+++ b/build/download_deps_and_run_cmake.sh
@@ -21,7 +21,7 @@ then
fi
tar xf LuaJIT-2.0.5.tar.gz
cd LuaJIT-2.0.5
- make PREFIX=$PWD/build -j 8 install
+ make PREFIX=$PWD/build install
cd ..
fi
export LUAJIT_INSTALL_DIR=$PWD/LuaJIT-2.0.5/build
diff --git a/ci/.gitignore b/ci/.gitignore
new file mode 100644
index 0000000..40eb0f2
--- /dev/null
+++ b/ci/.gitignore
@@ -0,0 +1,3 @@
+*.txt
+*.log
+*.db
diff --git a/ci/README.md b/ci/README.md
new file mode 100644
index 0000000..fbf5afd
--- /dev/null
+++ b/ci/README.md
@@ -0,0 +1,3 @@
+# WARNING
+
+This is a simple test case executed when Continuous Integration tasks are performed ; it does not really make sense to use parRep on such a small test system !
diff --git a/ci/input_generalized_parRep.lua b/ci/input_generalized_parRep.lua
new file mode 100755
index 0000000..fec690d
--- /dev/null
+++ b/ci/input_generalized_parRep.lua
@@ -0,0 +1,495 @@
+-- Input Lua script providing simulation parameters and functions for the parRep software
+-- All lines starting with '--' are comments and thus ignored
+
+-- Cluster of 7 Argon atoms at low temperature (15 K)
+-- see omm_LJ7.py
+
+print("Lua version : ",_VERSION)
+
+--------------------------------------------------------------------------------------------------------------
+-- --------------- SIMULATION PARAMETERS ------------------------------
+--------------------------------------------------------------------------------------------------------------
+
+-- NOTE on reproducibility : Only the REF platform was tested, and only with the 3 xml files from the ci directory ;
+-- there is not guarantee that reproducibility might be achieved otherwise, just try yourself ! If you want to share a success or failure,
+-- please open an issue on https://gitlab.inria.fr/parallel-replica/gen.parRep/issues ; thanks !
+
+-- This is used for defining a maximum allowed running time for the simulation
+-- The starting time off the program is saved at initialisation
+-- in max_run_time_hours the user stores for example the maximum running time allowed by the queueing system when running on a cluster
+-- fractional representation is allowed, e.g. 1.25 would be 1 hour and 15 minutes
+max_run_time_hours = 0.5 -- default if not set in this script will be 24 hours
+-- in minutes_to_stop_before_max_run_time the user requires the program to stop a few minutes before max_run_time_hours,
+-- useful for large systems if the I/O may take some time; should be at least 1 minute, and no fractional value allowed
+minutes_to_stop_before_max_run_time = 1 -- default if not set in this script will be 5 minutes
+
+-- This uses OpenMM as MD engine
+MD_Engine = "OpenMM"
+
+-- OpenMM platform to use
+-- AUTO : let OpenMM find the fastest platform (default)
+-- REF : on cpu, not optimised, not parallellised : slow !
+-- CPU : on cpu, optimised, parallellised with threads
+-- OCL : on cpu or gpu of any brand, or any accelerating device available
+-- CUDA : on nvidia gpu only, usually the fastest
+-- OMMplatform = "AUTO"
+OMMplatform = "REF"
+-- OMMplatform = "CPU"
+-- OMMplatform = "OCL"
+-- OMMplatform = "CUDA"
+
+-- We can define here in an array extra platform specific properties passed to OpenMM
+-- see http://docs.openmm.org/latest/userguide/library.html#platform-specific-properties
+-- for a list of OpenMM platform specific properties
+-- Internally coerced to a std::map<std::string,std::string> so use a string for indexing (key, before the =)
+-- and also be sure to define the values (after the =) also as strings (wrapped within ""), even if it is an integer
+
+-- REF platform has no extra properties
+-- ...
+
+-- CPU platform properties : Threads = "1" is equivalent to exporting OPENMM_CPU_THREADS=1 in the environnment
+-- OMMplatform_properties = { Threads = "1"}
+
+-- OpenCL platform properties
+-- OMMplatform_properties = {
+-- Precision = "mixed", -- or "single" or "double"
+-- UseCpuPme = "false", -- or true
+-- OpenCLPlatformIndex = "0",
+-- DeviceIndex = "0"
+-- }
+
+-- CUDA platform properties
+-- OMMplatform_properties = {
+-- Precision = "mixed", -- or "single" or "double"
+-- UseCpuPme = "false", -- or "true"
+-- DeviceIndex = "0",
+-- -- CudaCompiler = "/path/to/nvcc",
+-- UseBlockingSync = "false" -- or "true"
+-- }
+
+-- load the integrator parameters from a serialised OpenMM XML file ;
+-- no default, error if undefined
+-- adapt the path to the file in the following line
+integrator = { xml = "./lj7/Integrator.xml" }
+
+-- load the OpenMM System from a serialised XML file
+-- no default, error if undefined
+system = { xml = "./lj7/System.xml" }
+
+-- load the OpenMM State from a serialised XML file
+-- no default, error if undefined
+state = { xml = "./lj7/State.xml" }
+
+-- parameters for energy minimisation : tolerance and maximum number of steps (if 0 : no limit, continue until tolerance satsfied)
+-- defaults are : minimisation={Tolerance=10.0, MaxSteps=0}
+minimisation =
+{
+ Tolerance = 0.00001,
+ MaxSteps = 0
+}
+
+-- before the parRep algorithm actually starts, perform some steps of equilibration (MPI rank 0 performs dynamics alone)
+-- default : 50e3
+equilibrationSteps = 0
+
+-- the total number of steps of the simulation ; together with the timestep read from the integrator file
+-- it gives to total simulation time
+-- no default, error if undefined
+numSteps = 1
+
+-- Define the type of ParRep simulation to perform, and provide its parameters
+-- only one of the two following tables should be defined
+
+-- FV ParRep (Lelièvre et al. algorithm): only "algorithm" and "GRobservables" are mandatory,
+-- the others have default (c.f. documentation)
+simulation =
+{
+ -- The parallel replica algorithm based on a Fleming-Viot process as defined by ...
+ algorithm = "PARREP_FV",
+ -- parameter for the ParRep Fleming-Viot stage
+ -- A Fleming-Viot particle process is used for sampling the QSD,
+ -- without any a priori defined decorrelation or dephasing time stage.
+ -- Convergence is checked with Gelman-Rubin statistics
+ -- a frequency (in steps) at which to verify if the system left the current state during FV procedure
+ checkFV = 500, -- 0.5 ps
+ -- a frequency (in steps) at which to evaluate convergence using the Gelman-Rubin statistics
+ checkGR = 5, -- 5 fs
+ -- check FV convergence if at least minAccumulatedObs have already been accumulated
+ -- This may be useful if there is a risk of early pseudo-convergence for some of the observables when only a few samples have been accumulated
+ -- Therefore the minimum FV convergence time will be minAccumulatedObs*checkGR*dt (here 5000*5*dt = 25000 fs = 25 ps, but it is likely that the convergence time will be larger than that)
+ -- Default is 100
+ minAccumulatedObs = 5000,
+ -- parameter for parallel dynamics stage
+ -- checkDynamics is the frequency (in steps) at which to verify if the system left the current state
+ checkDynamics = 500, -- 0.5 ps
+ -- Gelman-Rubin statistics checks convergence of several Observables: we define their name here, then there should be one function matching this name below in this file
+ GRobservables = {"getEpot","getEkin"},
+ -- Gelman-Rubin convergence : 0.01 means 1 %
+ GRtol = 0.01
+}
+
+--------------------------------------------------------------------------------------------------------------
+-- --------------- IMPLICIT VARIABLES AND FUNCTIONS ------------------------------
+--------------------------------------------------------------------------------------------------------------
+
+-- Together with the previous variables, the following implicit global variables and functions
+-- are defined from the c++ code, and are accessible from the code you define below
+--
+-- ------------------
+-- implicit variables (read only) :
+-- ------------------
+--
+-- natoms : number of atoms of the system, read from OMM XML files
+--
+-- mpi_rank_id : the id of the MPI process working on this file
+-- mpi_num_ranks : the total number of MPI processes running
+--
+-- epot,ekin,etot : 3 variables containing the values of the potential, kinetic and total energies (kcal/mol),
+-- of the system for the current simulation time
+--
+-- timeStep : the MD timeStep, read from OMM XML files
+--
+-- referenceTime : the reference clock time, simulation will stop when referenceTime >= (numSteps*timestep),
+-- where timestep is defined within the xml OMM integrator file
+--
+-- temperature : T of the system in Kelvin, initial value read from OMM XML files ; use get_temperature() for instantaneous values
+--
+-- ------------------
+-- implicit functions :
+-- ------------------
+--
+-- exit_from_lua() : call this if it is required to finish the simulation from the lua script : it will terminate MPI properly
+-- and flush I/O files ; but it won't perform a DB backup so do it manually before.
+--
+-- get_coordinates(n) : function returning a tuple x,y,z containing the coordinates of atom n (in nm)
+-- NOTE n is indexed in Lua style starting from 1, internally it will be accessed as (n-1) i.e. C++ style
+--
+-- get_velocities(n) : function returning a tuple vx,vy,vz containing the velocities of atom n (in nm/picosecond)
+-- NOTE n is indexed in Lua style starting from 1, internally it will be accessed as (n-1) i.e. C++ style
+--
+-- get_all_coordinates() : function returning a table crds containing a safe read-only copy of the coordinates
+-- access by using : crds.x[i] crds.y[i] crds.z[i] for respectively getting x,y,z coordinates of atom i
+-- NOTE lua indexing, starting from 1
+--
+-- get_all_velocities() : same as gel_all_coordinates but for velocities ; returns a table vels such that the access is : vels.x[i], etc.
+--
+-- get_all_crdvels() : returns a 2-tuple crds,vels containing coordinates and vels : internally it call the 2 above defined functions
+--
+-- get_pbc() : returns periodic boundary conditions as a table pbc with members pbc.a pbc.b pbc.c, each being a xyz vector (i.e. pbc.a.x pbc.a.y pbc.a.z)
+--
+-- set_pbc(pbc) : set the openmm periodic boundary conditions to be used, the atgument is a table pbc a described above
+--
+-- NOTE possibly slow (especially if a copy to a OCL or CUDA device is required), use rarely for performance
+-- set_all_coordinates(crds) : uses a crds tables (as described above) and use it for setting c++/openMM coordinates
+-- set_all_velocities(vels) : the same with velocities
+-- set_all_crdvels(crds,vels) : does with one function only what the 2 above defined do
+--
+-- get_mass(n) : function returning the mass of atom n in a.m.u.
+-- NOTE n is indexed in Lua style starting from 1, internally it will be accessed as (n-1) i.e. C++ style
+--
+-- get_temperature() : get instantaneous value of the temperature in K
+--
+-- get_COM() : function returning the center of mass of the system as a tuple x,y,z
+--
+-- get_COM_idxs(idxs) : function returning the center of mass of a subset of the system as a tuple x,y,z
+-- NOTE this time idxs is indexed directly in C++ style
+-- for example get_COM_idxs({1,2,3}) to get COM of atoms 1, 2 and 3 (C++ : atoms 0, 1, 2)
+--
+-- get_minimised_energy(tolerance,maxSteps) : this function returns the minimised energy of the system, using the OpenMM L-BFGS minimiser
+-- note that coordinates are not affected, it just returns the minimum epot of the bassin in which dynamics currently evolves
+-- it returns a 3-tuple ep,ek,et (potential, kinetic and total energy)
+-- the tolerance and maxSteps can be the above defined minimisation.Tolerance and minimisation.MaxSteps
+--
+-- get_minimised_crdvels(tolerance,maxSteps) : this function returns a 2-tuple (crds,vels) containing
+-- a copy of coordinates and velocities after minimisation.
+-- crds and vels are both tables with x,y,z members, each of size natoms, : e.g. crds.x[i] returns the x coordinate of atom i, idem for vels.x[i]
+-- note that C++/OpenMM coordinates are not modified if modifying this table : this is a safe read-only copy
+-- the tolerance and maxSteps can be the above defined simulation.minimisationTolerance and simulation.minimisationMaxSteps
+-- NOTE lua indexing, starting from 1
+--
+-- get_minimised_energy_crdvels(tolerance,maxSteps) : this returns a 5-tuple (ep,ek,et,crds,vels)
+-- This does the same as the two previous functions but with only one call
+--
+-- hr_timer() : returns a variable representing a c++ high precision timer : can be used for measuring execution time.
+-- do not try to modify it or even read it, it should only be used as argument for the following hr_timediff_* functions.
+--
+-- hr_timediff_ns(before,after) : returns the time difference in nanoseconds between two hr_timer() 'before' and 'after' : usage:
+--
+-- local bf = hr_timer()
+-- function_to_profile()
+-- local af = hr_timer()
+-- print('Exec. time of function function_to_profile() is (ns) : ',hr_timediff_ns(bf,af))
+--
+-- hr_timediff_us() and hr_timediff_ms() : same as above but exec time is returned respectively in microseconds and milliseconds
+--
+
+--------------------------------------------------------------------------------------------------------------
+-- --------------- USER DEFINED VARIABLES AND FUNCTIONS ------------------------------
+--------------------------------------------------------------------------------------------------------------
+
+-- Some of the following VARIABLES and FUNCTIONS are mandatory and called from C++ (if it is the case it is explicitly documented)
+-- If not they can be restrited to this file using the local keyword
+
+-- Define here local variables and functions used later within state_init() and check_state_left()
+
+--------------------------------------------------------------------------------------------------------------
+-- --------------- FUNCTIONS DEFINING A PARREP STATE ------------------------------
+--------------------------------------------------------------------------------------------------------------
+
+-- TWO functions, state_init() and check_state_left(), will be called from c++ code to know if the
+-- dynamics left the current state. You are free to define the state in any way, using variables defined explicitly in this file
+-- or implicitly (c++ interface, see above).
+
+-- Define here local variables and functions used later within state_init() and check_state_left()
+
+local fromState = 'unknown'
+
+local LJ7_epsilon = 0.99607255958 -- kJ/mol
+
+-- after minimisation, potential energy divided by the epsilon should be one of the following (because the following are multiples of epsilon)
+local LJ7_ref_energies = {first = -16.505000000000000,
+ second = -15.935000000000000,
+ third = -15.593000000000000,
+ fourth = -15.533000000000000
+ }
+
+-- local dumpFile = "traj."..tostring(mpi_rank_id)..".xyz"
+
+-- local function dump_XYZ(crds)
+-- local fi = io.open(dumpFile,"a")
+--
+-- fi:write(tostring(natoms).."\n")
+-- fi:write("# dump after error\n")
+-- for i=1,natoms
+-- do
+-- local str = string.format("AR\t%.3f\t%.3f\t%.3f\n",10.0*crds.x[i],10.0*crds.y[i],10.0*crds.z[i])
+-- fi:write(str)
+-- end
+--
+-- fi:close()
+--
+-- end
+
+-- Will identify the current state by performing a minimsation and checking the potential energy to the known values from LJ7_ref_energies
+local function LJ7_stateID()
+
+ -- get potential energy, coordinates are untouched (internal copy was made before calling minimisation)
+ local ep,ek,et = get_minimised_energy(minimisation.Tolerance,minimisation.maxSteps)
+-- local ep,ek,et,crds,vels = get_minimised_energy_crdvels(minimisation.Tolerance,minimisation.maxSteps)
+
+ -- scale energy in units of epsilon
+ ep = ep/LJ7_epsilon
+
+-- print("Potential energy is : ",ep)
+
+-- dump_XYZ(crds,vels)
+
+ if( math.abs(math.abs(ep)-math.abs(LJ7_ref_energies.first)) < 1e-3 ) then
+ return "first"
+ elseif( math.abs(math.abs(ep)-math.abs(LJ7_ref_energies.second)) < 1e-3 ) then
+ return "second"
+ elseif( math.abs(math.abs(ep)-math.abs(LJ7_ref_energies.third)) < 1e-3 ) then
+ return "third"
+ elseif( math.abs(math.abs(ep)-math.abs(LJ7_ref_energies.fourth)) < 1e-3 ) then
+ return "fourth"
+ else
+ print("Warning : unable to identify the minimum which has after minimisation an epot of ",ep," units of epsilon !!!")
+-- return "unknown"
+-- dump_XYZ(crds)
+ exit_from_lua()
+ end
+
+end
+
+-- this function is mandatory and called from C++, program will fail if not defined
+-- it should take no arguments
+-- it returns nothing
+-- Use it if you have global variables used in check_state_left() (or other functions) that you need to initialise
+-- It will be called only once after equilibration from C++, before starting the parRep or parRep_FV algorithm
+-- It can also be called at any time from this file if required
+function state_init()
+
+ print("Calling LJ7_stateID() on rank ",mpi_rank_id)
+ fromState = LJ7_stateID()
+ print("Initial state is: ",fromState)
+
+end
+
+-- this function is mandatory and called from C++, program will fail if not defined
+-- it should take no arguments
+-- it should return a boolean : true in case the dynamics left the state, false otherwise
+-- You may create as many functions as you want and call them from check_state_left(),
+-- but the c++ code will in the end only call check_state_left()
+function check_state_left()
+
+-- -- value of the phi and psi dihedral angles of ala2 at a given time value
+-- local phi = calcDihe(phi_def,nil)
+-- local psi = calcDihe(psi_def,nil)
+--
+ local escaped = false
+
+ local toState = LJ7_stateID()
+
+ if(fromState == toState) then
+ escaped = false
+ else
+ escaped = true
+ print("Transition from state ",fromState," to state ",toState," detected on MPI rank ",mpi_rank_id," !!! ")
+ end
+
+ return escaped
+end
+
+-- this function is mandatory and called from C++, program will fail if not defined
+-- it should take no arguments
+-- it should return a boolean : true in case the system is currently outside of any of the known states, false otherwise
+-- If the configuration space is a partition the system enters a new state as soon as it exits another, and therefore this function can just return false without doing anything
+function check_transient_propagation_required()
+ -- state space is partitioned for the LJ7
+ return false
+end
+
+--------------------------------------------------------------------------------------------------------------
+-- --------------- DATABASE OF STATES STORAGE FUNCTIONS ------------------------------
+--------------------------------------------------------------------------------------------------------------
+
+-- parameters for the database used for storing information about parRep states
+-- The database can be stored in memory for performance, with a backup regularly performed to file 'name' : in that case define backupFrequency,
+-- together with an appropriate backup function below
+-- If writes to the database are directly performed to disk this is less important, and a large value can be assigned to backupFrequency,
+-- and the backup function below can be empty or either undefined (there is a default one doing nothing defined on c++ side)
+-- defaults are : database={name="run.db",backupFrequency=500.0}
+database =
+{
+ name = "lj7.db",
+ -- frequency, in ps, at which to backup the database to a file
+ backupFrequency = 1e99 -- this disables intermediate backup of the db : if the db was indeed not opened in memory but instead on disk it is useless to perform backups
+}
+
+-- functions for storing data to a database are defined as members of the table SQLiteDB (which already exists, do not recreate it)
+-- the following functions can be defined as they are used from the c++ code :
+--
+-- function SQLiteDB.open()
+-- function SQLiteDB.close()
+-- function SQLiteDB.insert_state()
+-- function SQLiteDB.backup_to_file()
+--
+-- if not defined within this file there are default empty functions defined, doing nothing
+--
+-- You are free to define more functions or variables as members of the table SQLiteDB
+-- but they will only be accessible from this file.
+--
+-- For example table creation statements can be stored in the table :
+-- SQLiteDB.insert_statement_...
+-- ... more ...
+--
+-- The database object is itself stored in this table : SQLiteDB.db
+
+-- only rank 0 manages a database
+if(mpi_rank_id==0)
+then
+
+ SQLiteDB.insert_statement_states = [[ INSERT INTO STATES (PARREP_DONE,REF_TIME,ESC_TIME,STATE_FROM,TAU)
+ VALUES ($lprep,$reft,$esct,$state_from,$tau); ]]
+
+ SQLiteDB.insert_statement_crdvels = [[ INSERT INTO CRDVELS (ID,X,Y,Z,VX,VY,VZ)
+ VALUES ($id,$x,$y,$z,$vx,$vy,$vz); ]]
+
+ function SQLiteDB.open()
+
+ print('Lua SQLite3 opening an on disk db named ',database.name,' on rank ',mpi_rank_id)
+
+ -- open database and enable foreign keys
+ SQLiteDB.db = sqlite3.open(database.name)
+ SQLiteDB.db:exec("PRAGMA foreign_keys = ON;")
+
+ SQLiteDB.db:exec("BEGIN TRANSACTION;")
+
+ -- create the states table : it contains the escape time for this state
+ SQLiteDB.db:exec[[ CREATE TABLE STATES(
+ ID INTEGER PRIMARY KEY AUTOINCREMENT NOT NULL,
+ PARREP_DONE INTEGER NOT NULL,
+ REF_TIME REAL,
+ ESC_TIME REAL,
+ STATE_FROM TEXT,
+ TAU REAL); ]]
+
+ -- create the crdvels table : it contains coordinates and velocities of the system for a given 'states' record
+ SQLiteDB.db:exec[[ CREATE TABLE CRDVELS(
+ ID INTEGER NOT NULL,
+ X REAL,
+ Y REAL,
+ Z REAL,
+ VX REAL,
+ VY REAL,
+ VZ REAL,
+ FOREIGN KEY(ID) REFERENCES STATES(ID) ); ]]
+
+ SQLiteDB.db:exec("END TRANSACTION;")
+
+ end
+
+ function SQLiteDB.close()
+ SQLiteDB.db:close()
+ end
+
+ local stateID=0
+
+ -- the first three arguments are always provided from c++ code
+ -- extra arguments might be provided by the ... and should be retrieved from Lua's side using the ... special token, converted to a table (see args below)
+ function SQLiteDB.insert_state(parRepDone,tauTime,escapeTime)
+
+ local ref_time=referenceTime
+ local from=fromState
+
+ SQLiteDB.db:exec("BEGIN TRANSACTION;")
+
+ local stmt = SQLiteDB.db:prepare(SQLiteDB.insert_statement_states)
+ stmt:bind_names{lprep=parRepDone, reft=ref_time, esct=escapeTime, state_from=from, tau=tauTime}
+ stmt:step()
+ stmt:finalize()
+
+ stateID = stateID+1
+
+ for n=1,natoms
+ do
+ local x,y,z = get_coordinates(n)
+ local vx,vy,vz = get_velocities(n)
+
+ stmt = SQLiteDB.db:prepare(SQLiteDB.insert_statement_crdvels)
+ stmt:bind_names{ id=stateID, x=x, y=y, z=z, vx=vx, vy=vy, vz=vz }
+ stmt:step()
+ stmt:finalize()
+
+ end
+
+ SQLiteDB.db:exec("END TRANSACTION;")
+
+ end
+
+end -- if(mpi_rank_id==0) ...
+
+--------------------------------------------------------------------------------------------------------------
+-- --------------- GELMAN RUBIN OBSERVABLES MONITORING FUNCTIONS ---------------------------------------------
+--------------------------------------------------------------------------------------------------------------
+
+-- Define a function for calculating the value of each Observable
+-- Those functions should:
+-- 1) take no arguments
+-- 2) return a double precision (any Lua numeric value is returned as a double precision)
+--
+-- Those GR functions are called from C++ if they were listed in simulation.GRobservables
+--
+
+-- Definition of the "getEpot" and "getEkin" observables used for Gelman-Rubin statistics
+-- Just returns value of the potential and kinetic energy, scaled to units of epsilon
+function getEpot()
+ return epot
+end
+
+function getEkin()
+ return ekin
+end
+
+--------------------------------------------------------------------------------------------------------------
+
diff --git a/ci/lj7/Integrator.xml b/ci/lj7/Integrator.xml
new file mode 100644
index 0000000..b1cb0e6
--- /dev/null
+++ b/ci/lj7/Integrator.xml
@@ -0,0 +1,2 @@
+<?xml version="1.0" ?>
+<Integrator constraintTolerance="1e-08" friction="10" randomSeed="123456789" stepSize=".001" temperature="15" type="LangevinIntegrator" version="1"/>
diff --git a/ci/lj7/State.xml b/ci/lj7/State.xml
new file mode 100644
index 0000000..03ed2d3
--- /dev/null
+++ b/ci/lj7/State.xml
@@ -0,0 +1,26 @@
+<?xml version="1.0" ?>
+<State openmmVersion="7.3" time="999.9999999832651" type="State" version="1">
+ <PeriodicBoxVectors>
+ <A x="20" y="0" z="0"/>
+ <B x="0" y="20" z="0"/>
+ <C x="0" y="0" z="20"/>
+ </PeriodicBoxVectors>
+ <Positions>
+ <Position x="9.983385059299602" y="10.174980427137598" z="10.020986348782928"/>
+ <Position x="9.770224100810191" y="10.037570850276143" z="10.293565548378014"/>
+ <Position x="9.59030488821879" y="10.02885267592082" z="9.931518707870099"/>
+ <Position x="9.648648572516473" y="10.349048230249956" z="10.107352368306797"/>
+ <Position x="9.766183439034616" y="10.306491926935744" z="9.758201173961519"/>
+ <Position x="9.93467672489858" y="9.95859189806722" z="9.682105588219425"/>
+ <Position x="9.906861537781493" y="9.817471151536196" z="10.038531569331555"/>
+ </Positions>
+ <Velocities>
+ <Velocity x=".06812021378266309" y="-.014536461211633878" z="-.05354852211603145"/>
+ <Velocity x="-.0003575163347591115" y=".004253822229571824" z=".10209696385565792"/>
+ <Velocity x=".01955586418134203" y=".01921201910093373" z="-.003143989026455074"/>
+ <Velocity x=".06026768071265565" y=".04230138984162579" z=".07773255456200445"/>
+ <Velocity x="-.05115033391156487" y=".045557891198200196" z=".03101618878886825"/>
+ <Velocity x="-.028826659962177814" y=".059965265059247486" z="-.07743277297400653"/>
+ <Velocity x=".06604694795520061" y="-.02419439117318234" z=".015470123640781708"/>
+ </Velocities>
+</State>
diff --git a/ci/lj7/System.xml b/ci/lj7/System.xml
new file mode 100644
index 0000000..ecb7213
--- /dev/null
+++ b/ci/lj7/System.xml
@@ -0,0 +1,36 @@
+<?xml version="1.0" ?>
+<System openmmVersion="7.3" type="System" version="1">
+ <PeriodicBoxVectors>
+ <A x="20" y="0" z="0"/>
+ <B x="0" y="20" z="0"/>
+ <C x="0" y="0" z="20"/>
+ </PeriodicBoxVectors>
+ <Particles>
+ <Particle mass="39.948"/>
+ <Particle mass="39.948"/>
+ <Particle mass="39.948"/>
+ <Particle mass="39.948"/>
+ <Particle mass="39.948"/>
+ <Particle mass="39.948"/>
+ <Particle mass="39.948"/>
+ </Particles>
+ <Constraints/>
+ <Forces>
+ <Force alpha="0" cutoff="1" dispersionCorrection="0" ewaldTolerance=".0005" forceGroup="0" ljAlpha="0" ljnx="0" ljny="0" ljnz="0" method="0" nx="0" ny="0" nz="0" recipForceGroup="-1" rfDielectric="78.3" switchingDistance="-1" type="NonbondedForce" useSwitchingFunction="0" version="3">
+ <GlobalParameters/>
+ <ParticleOffsets/>
+ <ExceptionOffsets/>
+ <Particles>
+ <Particle eps=".99607255958" q="0" sig=".3405"/>
+ <Particle eps=".99607255958" q="0" sig=".3405"/>
+ <Particle eps=".99607255958" q="0" sig=".3405"/>
+ <Particle eps=".99607255958" q="0" sig=".3405"/>
+ <Particle eps=".99607255958" q="0" sig=".3405"/>
+ <Particle eps=".99607255958" q="0" sig=".3405"/>
+ <Particle eps=".99607255958" q="0" sig=".3405"/>
+ </Particles>
+ <Exceptions/>
+ </Force>
+ <Force forceGroup="0" frequency="100" type="CMMotionRemover" version="1"/>
+ </Forces>
+</System>
diff --git a/ci/lj7/omm_LJ7.py b/ci/lj7/omm_LJ7.py
new file mode 100755
index 0000000..70c82e0
--- /dev/null
+++ b/ci/lj7/omm_LJ7.py
@@ -0,0 +1,139 @@
+from simtk.openmm.app import *
+from simtk.openmm import *
+from simtk.unit import *
+
+import numpy as np
+from math import *
+from sys import stdout
+
+platform = Platform.getPlatformByName('Reference')
+
+# periodic conditions
+# periodic vectors
+a = Vec3(20.0,0.0,0.0)
+b = Vec3(0.0,20.0,0.0)
+c = Vec3(0.0,0.0,20.0)
+
+center = Vec3(10.,10.,10.)
+
+print("a = ",a)
+print("b = ",b)
+print("c = ",c)
+print("center = ",center)
+
+natoms = 7
+mass = 39.948*amu
+temperature = 15.0*kelvin
+
+integrator = LangevinIntegrator(temperature, 10/picosecond, 0.001*picosecond)
+integrator.setConstraintTolerance(1e-8)
+
+integrator.setRandomNumberSeed(123456789)
+
+nb = NonbondedForce()
+nb.setNonbondedMethod(NonbondedForce.NoCutoff)
+#nb.setUseSwitchingFunction(True)
+#nb.setSwitchingDistance(1.0*nanometer)
+#nb.setCutoffDistance(1.2*nanometer)
+nb.setUseDispersionCorrection(False)
+
+# see http://www.sklogwiki.org/SklogWiki/index.php/Argon for argon LJ parameters
+k_B = 0.0083144621*(kilojoule/(mole*kelvin))
+lj_epsilon = 119.8*kelvin * k_B # in kJ/mol
+lj_sigma = 0.3405*nanometer # in nm
+
+# All seven particles will initially be placed on a 2D plane
+# Each particle is by construction at distance r_min of another atom
+# r_min = lj_sigma * 2^(1./6.) is the distance at which the LJ potential is minimum
+# See on the right for the numbering of particles
+#
+# * * 3 2
+# * * * 4 0 1
+# * * 5 6
+#
+r_min = 2.**(1./6.) * lj_sigma._value
+xpos = r_min*cos(pi/3.)
+ypos = r_min*sin(pi/3.)
+
+iniPos = np.zeros((natoms,3))
+iniPos[0,:] = (center[0],center[1],center[2])
+iniPos[1,:] = (center[0]+r_min,center[1],center[2])
+iniPos[2,:] = (center[0]+xpos,center[1]+ypos,center[2])
+iniPos[3,:] = (center[0]-xpos,center[1]+ypos,center[2])
+iniPos[4,:] = (center[0]-r_min,center[1],center[2])
+iniPos[5,:] = (center[0]-xpos,center[1]-ypos,center[2])
+iniPos[6,:] = (center[0]+xpos,center[1]-ypos,center[2])
+
+system = System()
+for i in range(natoms):
+ system.addParticle(mass)
+ nb.addParticle(0.0,lj_sigma,lj_epsilon)
+
+system.setDefaultPeriodicBoxVectors(a,b,c)
+system.addForce(nb)
+system.addForce(CMMotionRemover(100))
+
+# generate basic topology object to use with simulation class (and to be able to write pdb files)
+top = Topology()
+top.setPeriodicBoxVectors([a,b,c])
+chain = top.addChain("A")
+for i in range(natoms):
+ top.addResidue("AR",chain)
+for r in top.residues():
+ top.addAtom("AR",Element.getBySymbol("Ar"),r)
+
+
+simulation = Simulation(top,system,integrator,platform)#,platformProperties)
+
+simulation.context.setPositions(iniPos)
+
+state = simulation.context.getState(getPositions=True,getVelocities=True,enforcePeriodicBox=True)
+with open("start.pdb", 'w') as f:
+ PDBFile.writeFile(top,state.getPositions(),f)
+
+print('Minimizing energy...')
+simulation.minimizeEnergy(1e-5)
+
+simulation.context.setVelocitiesToTemperature(temperature)
+
+state = simulation.context.getState(getPositions=True,getVelocities=True,enforcePeriodicBox=True)
+with open("minim.pdb", 'w') as f:
+ PDBFile.writeFile(top,state.getPositions(),f)
+
+nsteps = 1e6
+# save each 1 ps
+saveFreq = 1000
+
+simulation.reporters.append(StateDataReporter(stdout, saveFreq, step=True,
+ time=True, potentialEnergy=True, kineticEnergy=True, totalEnergy=True,
+ temperature=True, progress=True,
+ remainingTime=True, speed=True,
+ totalSteps=nsteps, separator='\t'))
+
+simulation.reporters.append(DCDReporter(file='sim.dcd', reportInterval=saveFreq, enforcePeriodicBox=True))
+
+print("Running dynamics...")
+simulation.step(nsteps)
+
+state = simulation.context.getState(getPositions=True,getVelocities=True,enforcePeriodicBox=True)
+with open("1ns.pdb", 'w') as f:
+ PDBFile.writeFile(top,state.getPositions(),f)
+
+#serialize the state to xml for load in parRep later
+print('Serializing state : it contains coordinates, velocities, etc. for parRep...')
+f = open('State.xml','w')
+f.write(XmlSerializer.serialize(state))
+f.close()
+
+#serialize the system to xml for load in parRep later
+print('Serializing system : it contains all definition of forces for parRep...')
+f = open('System.xml','w')
+f.write(XmlSerializer.serialize(system))
+f.close()
+
+#serialize the integrator to xml for load in parRep later
+print('Serializing integrator for parRep...')
+f = open('Integrator.xml','w')
+f.write(XmlSerializer.serialize(integrator))
+f.close()
+
diff --git a/ci/ref.sim/README.md b/ci/ref.sim/README.md
new file mode 100644
index 0000000..9da4e99
--- /dev/null
+++ b/ci/ref.sim/README.md
@@ -0,0 +1,23 @@
+# Software used for producing this output :
+
+sqlite3 : 3.22.0 2018-01-22 18:45:57 0c55d179733b46d8d0ba4d88e01a25e10677046ee3da1d5b1581e86726f2alt1
+
+gcc/g++ (Ubuntu 7.3.0-27ubuntu1~18.04) 7.3.0
+
+Luajit : LuaJIT 2.0.5, download and compiled from script download_dependcies in the build directory.
+
+OpenMM : openmm-7.3.0-py27_cuda92_rc_1.tar.bz2, bin version downloaded from script download_dependcies in the build directory.
+
+MPI : mpirun (Open MPI) 2.1.1, from repository
+
+# Linux (Ubuntu 18.04):
+
+$ lsb_release -a
+LSB Version: core-9.20170808ubuntu1-noarch:printing-9.20170808ubuntu1-noarch:security-9.20170808ubuntu1-noarch
+Distributor ID: Ubuntu
+Description: Ubuntu 18.04.2 LTS
+Release: 18.04
+Codename: bionic
+
+$ uname -a
+Linux pl267-pro 4.15.0-45-generic #48-Ubuntu SMP Tue Jan 29 16:28:13 UTC 2019 x86_64 x86_64 x86_64 GNU/Linux
diff --git a/ci/ref.sim/lj7.db b/ci/ref.sim/lj7.db
new file mode 100644
index 0000000000000000000000000000000000000000..c5b3f9a5179b01b4754227a6d6eaf3df39af635f
GIT binary patch
literal 16384
zcmeI&TSyd97y#fiyXt7vuH9Tkg9Edcixp%M3q8!)+@r1S&W^LAtQD?SyQXPrZdWTR
z(tNO}r-BO966zr}x<b(erC@IbQQ74oD>Qqt5V}HVS8Yjft)PeEKkUqKX8!;DGvAq)
z!>lLYSswJ`z~(Jgr9r$HB_fVP1Y?Adj%_h)8*AhH+}MREM;%(2);Wy@=%>CyhE{{j
z4$wdV1V8`;KmY_l00ck)1VG@g2^=Wb^QP2PuHF(XU0dnTR_uk;=@~u5WGmFBQmopR
zJ1HJLnK1=V&>VK`kW`vW6)d?`EEPDN^Ct}UvO(DBQ6h>ajaW5lL?J5^y-Fb~Gf4z1
z6Fnlwtx!j<gjdp{42RvqqI0N1rEKcK9+d%l!e=eAF3yyh$rXjAYilam;?h^^ui5I~
zw9Y@StN%%U+{KJdrB(j<EYoseHOHG$Qn>o*;W0;wjqG(ljWBu}6V3*kx*u(94rRsR
zvMEKOr^2=Z)$NejkS;2*nVocYWt*Z<*=Ki4^f%rK=H^o!E*iR(=;G9q{fm*|C->zj
zZWqQX^{NxUk*clWhe4>iCA2o@jJ=5cpn(7gfB*=900@8p2!H?xfB*=900{hZ0%nd!
z2|RD&$sQcK`E_JRWB4Hg<y&flL@8k%?(oR^e;65u!`EO30w4eaAOHd&00JNY0w4ea
zAOHd&@V^KoaQc*lk-tlX*8jYr0vWrFO-99-Y<OWf%ev4&00ck)1V8`;KmY_l00ck)
z1VG>q1+08L`y3r=M0nEO@HDT72+zF(agBS4nf8C!ub-jS51iUHfRK3P>O<3!W+a~O
z8&a;{M%K7U=FQaHjj;hD+-O*vp7@>!jn8^{u~j>B{?_&-B}hBh*O0fsip1uMOUKTX
zAZu(SbGoaiN)U<gX~XM9?Ftc=#k5syTd!U9^$j$xK;q@Eg|qXjStfpyWba4Tm`J8w
zt!}GpCg$VUx_fQ4L}<BN+*Ef^ywssJ&Ff?t#mQ~65|Oq~EBoAmkX0YaOrPRCCqE&=
z$Kn&$<I`B}EV=i>W$j|%+aa%v#Lk(^PCF}+*xEHG>A4qKb&<>lvGaaVCg#I)>t3dp
z5VLg5oNRm`o`@TKxAG{`uH-Mh^71wk9fObi?y;I&B=bwn+pg+<VlLl3lvx%}ggcM2
d<oizXR{o)#fgsXOsuim@GpUv}Q(M}PqVKcAQ_}zd
literal 0
HcmV?d00001
diff --git a/ci/ref.sim/seeds.bin b/ci/ref.sim/seeds.bin
new file mode 100644
index 0000000000000000000000000000000000000000..fc3504f94422bcdb14a040dda024abb0a387d279
GIT binary patch
literal 16384
zcmV+bK>xq$Ex5=WBPA8%p%<QFX81pkFT^M!LPo+94Hi|T#hT}U6JyHaUn%eAEMFmq
zO@L6{bNL<Dv7oBECKc*u(gxq?pFP8K<GD{x_J~z?*OVR3X`J-PTKwRkHU#q8*fLX=
zWicOrIl^+zB_2?=3*Qfu23imtF|E*`w|^AV_|0Tb^avkBZ`%+Xw00qg^I{>Ou&7Or
zh)e0aR_V}jQK=-{laHIHvE~7otCH;l9&!8tim;X`Vt^TFH3lW2q9|L_!}lb16qR2Q
z661?iJiVlyKU1u<4lk>r4Q=l_q%8iRb>o?hG;@fYci4)L4N>>vKzbP(?AzxO_=XS*
z%h#R=Aj|s$1f){bR3lRuVAIaRb4zn!ucXZ$xz;5$sJb(I4NtBi2Z$9EN207}0GioG
zN?JLRvKdDYxkT{=0<>z}>e77Qjg4?DvAL}qh*;0#-MRxv%3N1);U9-XxOw{%L^;W4
zFmmQ@tBk0HkrBEn<fWv7CV@o_L%Zu;fMTKc1Jb&ZPBd(eJ(Ng(V3y>~qX2&GkEb}A
z=S5HMF%FrfE5m8p6`MOER1tmY(GKBK6=^yp)ILfKx6;fCa0qsvN%?rOlpwR_zCBze
z__wvNwHh)jwSuf?sD_tfoBPmRQb-PhlLAsHI5fMJaLFqEr!*qhI<~e+sE@q<euRA|
z1dyv94C~t9020Xd=nT>a<G#)7ey~T3y3jR=^2flX!y&Wcm@%O7(^zW=G~RH*e37q%
zLU0pFst*KM3W!a(pib5mmEkfm55UX>FDuD#uL>j>`OZYm<MNFHL2782v+$^lj)a8P
z5^a^cucrd$l03gIhzmsYzouP;O4??cU_)yVV2kfjLVDV%GioX@oxBuDKn6?PT@0|i
zMPF6O2wWY-8*LdO<rCRf@nCuHp5>|VFMqF=)Z>}bMHi^ZH$C)Xt^)L?UC(9QzrR=h
z$#W-}5ZmtyZ`4Nfh=;WEHbPE~ffAC!o0&j<!o%Bh2gCZ3sz%&p!Je{<_K@~MnIkud
z?*Fg#Rh(wb5cjW!K^F|4moKB=nRI);Pba9GFs`L}-`*9-wj8WC4wupnWkurE)sqD*
zQ^qRdVO%7`1bW1M{CToea<6znYd>A{&^gkGRp2v`thXZD_B@)Q=OH8Xp|5PgR@Mj;
zH0I1!)+3#4cfwtzC?O|2nO0B}MBL^+w-f;U@OTRSSDsI;PPBj;;G&K^Jad&hP`<Xc
zn>Jd})5AgrS4`^Js|nk#Bvz+?a&wTodWk%EkNX?!3p%cK@n~Y#2oNTf<STg9)joO~
z7ooi(>+Be<QkdGI^VUS5Ew@xWEP2&GAn$qR*u46&os=Cix6)p}V4nMUU7s6Td??sU
zN%HKfciq}+zb&G~b*2k~v<3+2V6IJmU=oJO9fr#R1bu%$4y?I#+iQIKEHEQi?oRzh
z#la7RNw7T_*mb>^K{{}%eOiNi0MS?+?%DEH3&YAJnCYA2Bv%H+|K|&}LMHltuXBXD
zoCHHkwu<pz0Y6o*Z^5jf8S7o&@c!{!wOLMc6VdGBJ=vdj+}l}**FsTR#e`}F$^;Wo
zzLtkOQB`zJC3)Lq<AqEtf$52c<bI*YjoTOhtOh8Sxe1wWQwKgpO;|AmtWfP#S3b$s
zEDdI0ohMe@e;?ZW)f&wv7PlrWJ0~{Kzo@(59mR{Ebo08%lvZw4AjXdmGzeBdo}^aL
z{3#-M^}v%a$`#!gmLlehiL#+Gyu$tpLUa^pvo&D00d(AKhB^jy6|17Vcrd2tSL=@T
zu0T24N_ZGnNhkzIVSq7LP0&9HYMVWAZ!=1~EW8jqRXU4N^IUBFjnU9ZI`Sb|NjJVB
z7EW#;mIsN7s}0$lby8DWe2iXh%u?k%!HrYTpVfckw1!D%AS<?c?<MXkEo$$O3Fd*`
z63y9B_tBV@vLcO9>>(}Cb>cWp%eu_d_dK8t1sPmN7=7?}u9n5B(l-h4I)*b;xE_yx
z!zkntQjYvLIRByz6TRqtsz>X*boDOdA4qMGgfWu@C|eW7hlt`*TGmry<@mrrK@{~-
z@;~Wb6cPUXaoPJ?n|ojJ_@aTeC&2M}i7qznSWxCbGqwoK-wy=GcebrN&mvaK7|*D5
zoy#*KLy#Jg{|NfPpl^}j=#pNDEUn-lZ*6Wvoe#q5gwTH@9#9MLxh<|owUpCM-&mH{
zw?{Dg)}8?K8)m;=c#o<WSy{B3emy#f$8=jqvm74z<Z_E2sQa2+qw3WNeXqm-FOutx
zSt}y5-cG?h0=WIx35kigAo`U$*bUQ1qq}jwFt6vBQi8J+8<YR8q#FTX1up<7b5Kc-
zpF+?#@gZivpnisi;X!0CN`Y0>+Hd{UnmD)#a*@)X4YmnM5ifmWMK`*XRlB9%#SP!J
z*!yr=F3)I|IrnH;r;#~#OH`fhhv3=-H!|O56&b3jg37SW<9@U&-g8u`UKjVvV3Ut0
zpsd3;af98RRpf4PVOQT>AjhsWerAZQHOHunW?REwM`wL@uK`X;0qs&tZRj$G*hwFt
z6!t(Vk-pz|iIoc>x?r1_tT9Qiy}y_LGuF$crVH2QCuTWAsdN+&*ku^@FPE12%r=Dg
z`TEwbX!~-T3v_gQ3B4?&)z`oiwYglJx5hUVh*JlqThLAtBO!;W#AEkBUGVzpP_crm
zD%fzkG8;#x9_5vaeVc%!!)KR?iPC|kes(<&7szW^STTCs!f#q-{Pb=|a6COO>)<`U
zVO`Z=p(F3%m2z4*l0$s=i)zT`1(F)w*c6-!Ip|QMiwLBzBPwleQ#(&V;vm_k$NX~c
zP*L8x{E}zDn!J24bPQiEuy#_{meRzX&aDC%o-2auurYBjq<yKR`_veDxmPaV_ybc`
z*Y)J|7xjb^LC=Xnqo{a)=5dv_oDhbcx_+tWfbvBpAmBD~od=241pu-XMi6*NOq($=
z0-}akD4L=~KX2?Cv)f~DHytN(hyL(X)(8LE_<mcy7Vq+8oa$*c^>b!?^ja9ft%{5-
zq9iy!9aMf=yFQ#+jGyp~D%`tPDP@*{0ACM}_$As(r_4qM{b_tMEe!Ivzoc{s9apFc
znG6zK<6*18`Qg)BgKX8Jbr;mUtbWo~({zX09^<>e&dz?Hyo9!Rx9Pui_miH<b_Dy(
zCOtRFRi?&2?wF5X6Dx5Rt}~u*qt~72jNM<^64$Zq<MfUJp>${>0r|953;7!kZO5Lk
z2*1xIzR}0)suR8k@n}GmfEDOZrETBWAglQ~mwx-!EL#-EK&wn?C;ydW6X2+dB4j}R
z(M^HX=b?27hiea6w`~DBZ=YB21PQ)V6iq6M{AX}IlpL+|!*~qHiguXl{3oes&x@mH
z0V&xx{2k<?yd8uf5l<a}r`WNkBauQ3Q)%^KTtE@y;w|YKJ3l83{e*<iwAw6JC)$xr
ze@^#ozHe>uoSL}A^!wt-*+Y!kk9QW<s;+^d>kwymo$CfZiugn;EY<JAaY3QiC*0R=
z>BTBs_M;i<t*!NnWuW(%bV&go+~`%#UYz7cU+Hb8-Fh9P2k)n>FU%1rw+ZYMlO;bv
zn;)L*A4Q7*RsHtXMA+2QJCf%DRIA3rz7_;%otVC*c($V$tPmGMj)Q&2q(fGu&LY3!
z*<?)xsNY{WcXPek#!z~kdV0$1(aqLjb@!S#K41NX!P#;*M3)@*(8(Gg@C`Dy4?5CI
zA{W|o%cLIAw%)n)429Cj0NPxMl!l8j1GZC+=jOo}@}-q=XCTd*Z^yphzlKompbVSp
zX#xuWnL&{H)vM;a%QxEG*DC}VC&y?jno`Lfc^UX-t9Xi8)ixwjVftzry#xT{cTw4i
z?&E)w!V%z%bWkMYIevwWw*7Ui32-@zDvL*@mCS$welWue1(>7O$KL@fd`h5*mjP?M
zP=@%2h;f>5a45_t>2T4j6|9~&%0!x#tAbTf4t4D8Z(?Q7M~M9&sEtjsW${^w4hth?
znqRjiXs}5^P+rll8CFEElv3NqJ&R)6H8|N?5<5{ZA!kWLPY@W(Vh2H4^BCTiUK8pG
zf82yKrKai-$*WPY#0O1PS_Plg4{Skbn|;)dZwTsBu4H-WY8~U0`~0kOEaa*Hex)6k
zFh0?d8Bn{2#c}Q+Twvs%G6$BflE}l@wpG_!xR9ki1Yl1rYC~oaAB^G?Emc-OI?m8@
zS1f@8{30ylb_cMIznPOb#)s^;0l_db(AZ;sZR}*weG`$J!x3%?i?swapbAwuwTRn$
zE=Zb=^f@R_)1Oa8vKEjq^i$E@DkNulWS3qZ-n#Ch$BV3NyyvwP;9agx$Ccj`vTt2i
zM=k*R>I*PV7cwMqw#ZBJ&62fwXu+3X&xjT0+K+#Yaaj4LtQuS-%0~K-tq4*5+kF}R
z(F<Wd1wGH*H0Vu4k(RS1>l`v|S3h&Rc3DK6LE-UoT{mMuvA$c7oz*_w?HruI5A_B;
ze~C#51#9UbWhzI+B^PSHnQ@s!n-6A|dd{n+)^mkBF$5LHy|+wpaA0(4cnc;WgyN_z
zT2D#9?X7;wpx5T$H@{hiUwnwbuNRi&e{pyJJ(d+53Os#^_^=+FW;u)8x90LN&Z!gi
z3d3&hdZvX<GHDw{J?y~S&49ovtaFk*1A+5ExL*~P>z96b%wzRH5vhQ0Xt(^MWydl+
zD~-a7kpJn{v%U`vK*x)HOWN<ILD%X<d<L)<ZCc-|Vs-GR*BAJv1IPqrF3ezAIf_&n
zShyw8e4Fin!M(aL)8&F?4(MI%<6Md-&I_x)cG*^v(8HQV`>+XBYL4}@v^FOd3OJyh
zx8G6|nR(PXdpt%Efk{_OyJf`lIiHh-B>cImG0hg4l`rfWHF&`aJYX%aAmSoCTP!}e
z|Fk|2#xPeQe_b{LM_Q#uL(YDwj5|5JE4T}R(k)F?9-d?QfihK=HvocL;9OTdxjKP(
z>M3IX$h*Raz&~5;G?WsADx3yM`au-0e$+V5)kpeEK{Jhs=Fy4Ff%+{wU+TB(N?D4U
zV<VVeW|7MVx8ea22uM9n#86Aj-Iqs+EDtBXiOr^|@(qY70ZQBV2`7v0s?A{@=?spH
zomP#L6^aUz^WRs*O|LtUrzjQENLg`*q(b*2vR>w@yr*k%*ILSFe{a*!Ao;i<$To_P
z=HCn$qPLDT8lniHuUBHdV@$5W<i)7(HEdinKa_|;=`&I_I&Y+AWZwHuOW;_Wnb-sM
zg|v`CLB~U(gpDsv8i@e$Nrs2%5Qikj()f;|CT^7|5jWRu4My>?^aN-sjieV`bze9z
z{Bv8vvs<$9Hqbf%#2pjI>)18i&J^Jr!JZ?l6WXJ3){FJ=<PSY{?E7S@;z5DA0k2ld
z%>$biA|@4$f;R--HDx#NMXh`U6B}nFchotO=M?zdJ}7CCAr9_!+yn4Y{3vTtoW)P$
z*vCM^XPi?N&urky^#GlwQZU!l?+#-06x?iNR7u&$Cs!N}H>AaY0&E>wOrlg(V$#wn
z8HBLfiwLG8*jy7;jmOM*bZUJrXlRhflQ*5*S%?V!21+f|X;+;&Dz4FR=&dV?dMLYU
zJ@2@2*s^cbUY6(kI-0kM@o#hy&}!GMpZltfTG;N@j6IG^Q-D}qfgs*pCG~<BoxZqD
zqQOg4$|SKDwn-WwrOW7`+vI^IYKh?l;rr{N>yx|%8ZUzKwe|}Y41%22lte;A`)#A$
zu+VL=0DUXW;Y&bkG8JLIswr-4+&uAokuIKj=fc_-5{!D?C9P-z%ZHvU16aXIy*OZT
z5B*)6g36d#8ZR-61qDWi)J%n?f)=AD1~7ti?WjMk5Y(PVi3PR@e^Ow;tZ0r>O4h)_
z^jM*Gty&g3-QL(aY(ct+vGrNn0@V1EDW|LN5k3V$h{&8F;)5zbtsfk$%0b5h`{z9?
zRWH?YS;FOF-1U-9V!M&y61YcDYC1R5d(}ipDAJZI;|&?8Uk+emMgus=UEsHUZ0cbP
zZq<Gj1dQa7=2H+uvK-@4zTH&3@v!H5K^w}VmeRlwNbP&7PWH2gUIIrgN{`6N;%Ih)
z#FOM{7j4Px`vG>OTRR|Y+Lxdz&+};ty&i#e(c&`!G%-*(q{_{NEGz*wwNZO~)W%~0
zDI20WI~9Q&EBu8D2h+`O31|Qan6yaE<rpnZM75QTL5MHMe}nIu4F0#NBX`~g#Iy6U
zw2i<(XtEW(qyF=7%#RZ`b5v2Pz!wDRG0)i{#SnD$h%nT=BIq3p340+rr4S#9ARGS9
zRQ9p>jRpZ^Ch<D_A6s@WSXat(MiNN4T0DOxol_sIvL;$SL`)%%+=YGV?I=X))<!
zl~j?aRC%;wV>}O`cz>_1cgA=I|4_-&6$vwuy?rk!JI>FgT6Iqw<8Y6FWl5ia=4yJ;
zM$O13v>~47#g2dtQC#W$7(C{3vs$I=a?c#(sw#XfpP+R5A0ah%4*m<(*>8XDf3z7T
zMiCgrho~nst!WCUL%96nD^mY5Dmh0}THD}m$&UIn42%vS66+T;9dUm8*W$t<B){?%
zJk(}{?>fR6T4f?DygUN6W)3$;j$MXbM~xe@53CSModhW#UOsf~V7^-UQL<+?uhC*<
z2J=B3kGSlXL!@nB7dhlew;DGGFrC`~`}3GC0ZfE|+%%`d>$ro_c20YDxGisrmLUq}
zjLy2dF7vy{y}>p>{XRW`*0$;y5Vtnp-J2eAcD<T=9V82;Z_)20-yw^Tq1#ZCMdjn%
zS%NxS?85OcU{(Z|SDh%{uKM>UVG!X-OYm0Gq6H7P3&#yvn&s>>%>2ur-zMrx$fyd>
z4(ot+`+UocD}?>|RFswczz?LvO+q|agWmw3#v9q~g=A4%l@0&9d-mO1jnAHmC^a?1
zQG`Ln&_P!b<`5E|F%|mTnDJ8ebEDc#Zdxo-*wsI=<8ly-NW`Zq^vOkWX-{&$Pa-;|
zE>vN${hAYhna`5?{C+O3L79p8gE1CoG#HRhfg<@lgb}?w@}-Xy&?9MAC;#<`&A;tM
ztsuQTXxK*qlH`im@pHxRE>lDwlJq}wk|>_P4RxSAm$>YWB_vE!u|fZ&l63UN6w5Ud
z^K%46BTmH@_#S~7=j_AQE5LI4CtUNHM$VnBU<s-Gzr!e~Ub3b2frB^;b|#QDLOK%-
z{KRKmY>qMunA^|scLlalDZt*Z6vvbRMK90{W&hXCyW=hl^#xQO_-xi1H2xPit3?2(
zHzFvo)eu4Lh25l@^OuWjmOrH2Gy4K=^G<TkzA{m81jQ|ajd+w=>?Y~ga<~fk!ysRq
zlog2x=PNmpRF949A9!tdD~R7R=je&-k+kIhGgHxxcmR(}e((vjAc8qDoiJG{g47zk
zYZv><GPnbEoGD5?>RdD?UFRNWhlLYechjX%Q<W(rCj+HLbutQj5L{buwyP)VSh?`h
zYV*nF-i0_R$64fwJ5<4*MsAQc9XRqs6IZ0tn$OTsSPkHxbry&>y*eWYc$wUwvH1+j
z4|j9af)7(4kIjqyQ1X*tY4X@ul0AP>m$s*cu}WjEwG(t&#Y6EHv#bEAEUL=xk3?gz
ziH&{u6{-kII*WNAd}Du9!k7j7he^OB(ng^=z+@qirmZRXc!~SEC2U0wp0I|{qtVFg
zq)BzUiQ@x@%a$bE9laanwT{nS`)h(Q1scDBAm%`?$TiZ6U(EnRL*&JGPv#kTX&*!*
zB1hY=vJ)KcP}MWpp5mqb3<G>MMFte{5alO_nh7DSzqb(h#{DvNX)X<~z1NjK71#6~
zIy`OrbKWV$gIfBH*=WnGeyN@utfUcY<vMYekxX=Mc0yLl1xf&fLtJZKfUR<3Hf}rF
zYId7v^^eF!a(5L|@0AR2Unz<x4F#Pqe`J8lFH7wb@pG~#ku{2f`bH>#tLE;YyH#Qg
zq$tqz@|lUR-x|?Ecu?o7{JrZQxRFEwsW@GA@f1|d2Tlk$#XHt<BUt>I^x3|vK0DSh
zbWj#7uI)v%#q|#Z$)?lyzi(vs3GT`IB%Thb^aLBC(>;Ct8!Fb#C~c&IR*YLC<1C#l
z9_dzA>yv_RX@<xH<&wh15~6Ne_hIZEqxPREROZC$Bh16|3&+1P((t3&j`!m+H>c%#
zc7^$Rh%!<l%v4QW(9R8p)CzU*k6*shrl4;Oo#a8~xo%BK?$z~WWf>kR)8rSW%%<le
zB(zf10fi`n<W5~S6wPiXMHXMk8aoHw^(zfVp-qC4Vtswz@hH%_DNdZ9wc&29)_&iG
zpN+j4g8pZ{BE$!3(X;TCpvsV%ucE$_z5~`k^9D8dci?Gy)8{Bx7KO(Q+D%E5%xGe!
zLx)%h1&fJ|96K7fZ5+LSU9Q=L0GlUM+gw!`bP+La*R8K#->LBXryi~;7s(I!XF^7$
z%{)>-Y7-i9W;`X{X)-|$rviH1#KQyvxMs{|cUfV=ssvo6y_YN^SwMf=Q>M2794^M?
z;e14NeSs(#a%yv1M{zS#CiUZ~^TJ%Z=YX!YXIYCWd;LY}j<7Zw{*4m(`X)V@fP_iN
z5GP3d0{6N!1KO4@xN}U_hJy~9!sbGu^CB0kc<;i*gtD-YEHJRKApmuBa#g6Z76_Jv
zP2d@1ie$CW#<}R55S4i>a-_xb>|cNqZE1$A0u=`4Kjh82Y^~q^rnkeEcowyWA`~Wk
z_+_zrT1!H0!RO>SNRes_t-vWne|yCapMn{m^8h!mVN~`3xw133@?4Gi`(<=HV4p{D
z56Oq|kxKnYfo~}ODih`SV|*EGK%al)?;l*l5B`q?6c%=*_`&Hp0)O7AW|S>HqDT52
zaO5L;9__jNUolh4Cf!b>EEAyqN-(nl@)CMa0T(*|Dc(ZWvM=_Hao@+y^)vd&a^20c
zdN>fb@9(N3!R4#nTD)Xg0*JuNh+L!=OeW)Lic&8q@hmH;np|PXTG@zF=U#r8hd`pg
z>nwIvCkKlLr0I@%^;A}O(PhM{<GE>csNi;xFcqSpn}54`C*s5Xd_%js(aF_Ke0yo^
zHB1kE0=8J}Kd_Q!AiXu90-EkXIv(;H6h9)%KHRKi&kc7!83yIm4_Mqdv*A7YGgbU!
zV(yT<75pBMa;nygvhMu_xN!}_6e!GS{*uSH^gBBnIc`u~ixC<vmjuorsn7_R(o_`D
zCty>4%NHddHhpguAP=+}exoSgrV&=cmojKSzAc(5622Y!zz4kEzsE_?a&Z<ZCxs-u
zGU%Kan;?(>-E5Js8=@}@DWmki>^(VOipYgSKP>gc%n`6SL>>@Q^iJ$mMq?&eDv@o*
zx!PaI$5R5Sk^@LkYRm+<N9p9b-Xe5w;3)BJ0eyoy>^seOxyu!?=fa$$?iCmyEP(-o
z(?I`M7GR#o=wHs0H7&^noP#t^+r>u=W71j}+Z000e(c8J3C#;-H#JUyiL9bj_bNv2
z2RVJ>9}!k&T`4aERG!C<Y!ka38EVG6(DkOaeoK(WSc_O@I&2MU2-Aq+FWMK$=Thg5
z|0I=G!a~#rl!^@S`g!S%q|cWKE_4r+DX77NI+0uqDCst#c|=y~D=`d139LTIg!SNa
z1*AA%vQwW-w*We9kg@=yH3>~?pJ<5X_Er>zCbD$1jQ-1|^@`#Pb(SS^Hn$qWEqg+t
zcSxOsK(nmu*Ni#JWy-#wuh&v}Vw}2|BSLv`KuDPbDu&4FUo41x$Ao&*B9>)<0fRsP
zdUF3&_chmGcAG~j{-95-*wWr5#DpfvLI%)@ZS&t!7hkbd=@<9-;Fx?ia)<OI_wr-q
zN<y+@pZaKWbM&2Q>IKN$Cn{xVdG@7t<ety4KLi}BJC>VJzG>PBD><zV^5+-t+cp&t
zNGp5r*koE(7jlxq2Z8a^n5wmTnxjr7yW~MT`B<)GX0!fYAdt30nTT0*kusY^zw%ho
zbEd_!#**;11kP}i43Z_2W%i>m;OlS@D*oX90<U2960dhLk>9t|0rtM;eQn1P6oJM9
zSKdz`|75|$C?jk*<DE64Y;95?d`3A)xa(f>-lR@Rm6tY0v<)E$I;t%I^&1D>#jvHm
zL(;ir4i1V0bldFLsP<uHSd&Ec)Jt$xj)ywi_w$7g_AJ67q1k14s=YDvmAke<)j*xC
zpOq5!e>W;F$)#bZSB+>PA3q?WeEcJH=}uOWU*NtHn=&SgPR`{n_%`zvZbe+3zh6>u
zg7hYMmUOGdPznGiXPhCig$;R5wm<e>=%y05SVmmY)50UimCU6xwCQW-nrYON-p=nw
zg{oq)uY`TtbZ+TOxq(@KX9=RK5plEBv%3ntGMfwjy<XsLEDa;+*=}Wp@YTvtr~G{h
z<m~6iDgx;Xzo55Ji81cy^MPZOjZ}9sH`iZa(!*~8HYb)5fT?w2uAXFoLT9MAm+9#^
z(158_h^MkJM0QwwsF@TeEv*!<dx;Fj$PjM6+B9ArKM+(ZHF7CQfMuWuF#yK%`Es|D
zi!<h=(vFx=Zm8RncJ%=ix|lG^o)F5PR#(IlRAm6{MX;fuEt1r}m9dLD;(Pe%Ys|i{
z_!DA9ACk;~KE&#FK{$;LCLB|cTFwiM+yM*-+?j)1p(w*BYT=8+6zwWLIL_`tP~pIx
zGeqk_c^mn5vAs0G(-%AJ)_F9T)I8>PKrns@Wkz^I6)N|E7>Z**L#)Z(T8Nd)j0wn(
zpI1TbCbsVi&1-!6dkQUE+WRtdsPTWwlOQO8{>A&5V5x?yCcD0<exF{FcNAVGOQ~cl
zP0??$7Sanyl*Htv(ox}<s6MaD1puG`u*|1+cqn^m6`bEt_8qMii%UzNN1I$lAo~Ft
zMn$wB=JslcV5X%hp=glr!8VySh#6fVvZ^YM%pC`j+axLVh}kZcqt|K9N4L6P-42n5
z2B_&kl#0r78aIM68eWKV(ku%)Xai8!yL8}=H5w{Eze5U0)8^T0YKvAxMjJI#S^`z8
z$U02AR>weTcwG*<?0h|lKo3}nUrqGqq7YDsK_|b@6qkc=A~=RwN+f>bX-D@2kWlO-
zhFXn6py!KS%SIF;HyvlEXRQSkxt2zGX~0G};PNP?>m}wq3=0jKXoa8_QG1Ex+a(sg
zpR{(r(-5L*65loX#gQ#W!r~dPo5+9x5*h7nr!n<A0{jF}0Fq$bU-{=QDBCMIF{pY_
zAg^{noEi5>1iTBE4t72Y1oE2w>a+bf1LNF0a<f``OS>*#Y<OJjdJ5rW&>WY8L-c49
z$acLuL<_}K!l9H#i1iZIb(Ru0!T-@)SAc@2V$$q!PyWGaATR*Nb*gqmz2Pz`%aN#x
zYXb&LRRkhNj>wPSxxe<Ey`yhxKsY~&7JiW?95Ath=#&b8;E=pb`1VkfGD$4Zjm{5>
zqYhGxK&l96OW(hQUz@(I)FmbmwOLOAt~huijQ9LHtyob{#WBq88!%rc3Drv!3&#ck
zjKMl};PKUm{cunn_OVxp$K%)-@h^~iwzeYwg)tR==p33m0j>ZK=r>a!l<G!Gp(dA)
zLxy1-q4nHEF^K&|udGhHoPGy^*lDuf?CgTV(fv5Yek*KM=Taf9c_MG(q)115wPSH5
zT?#$7>Rvc%-zzj6LdYQr`AglS6?j9fi0r(PZTBHT(BhZVBVdsAb}pC3u(sK2rywhv
zgW+U6gNp+42zRw~EDAG80vIdsUWK|6l)>pr?fBLO%PIiv`vc%KG%}cTS@}b(pK%rv
zIEz$j<N_J=m1L~(3)WFNhtx%Exyn@i2wJ7bmRT^mqKvDkW|>rbNP<6;d_|z0j0g@7
z^ETo6)*zP=$tf|z;gAe<box({dZi63Ml!2KnpO@ZEZeA@X3jat1H22wm)Ew6*Cc2I
ze@|ceV4A_evHQssADq9wWGfqR>a3`|YhK3`3kk0b@QmHnml`2ri*<XRmH!iwl3dux
zz2k)QN7IgBWWs2IX(#WM)5dPfgSrv!e+cG!ZM-Pa?UX_fHDW@)96?bh;Xxol!`z2m
zCN{}MsLW4C3olTx!kvr69yly$Mur-CH5Mp>LvX{I@)v<#_RCIkG}8c7OTYlI!{VN-
z-F0%>pp{LhE6ITxt+_E&sf#II{q8oItIwB;-FKKpO;+6oRAu^B?ZUwLuPcS?>CtKK
z64AjHGUSzX!iPIi@k~AxG8~NdXM=S*!#z+(G>?pqKC~6sf`LE--5o3SzD1hWB5mtG
ze+f!D6=Of^flN+nG>TpNms#MW&(DD&TcQkno2=Wgmun()q?u{YrFH;}Xab33%aX2#
zy61;B-X+NB$oOhiIU;wH6Jwv2{kI_|bjSZ9%Zc9bQJ~9nX*BBj_`2|#==M*`+^NdM
z5n1WjhTpdm)&F!a;1+7$O#Eaoy;b)*HzMb<9CEz5Q$Gr0HVe))Tt%pXETl9l4n|SG
zSN!F6UAPOR)Rmq64eLGu5xCwcL)c5xbP0yux1L+i&tM~^U7HFT!DRHDJa`_>{7iT#
zl`x2Z$b;WeNRg!FQrOC4i3^>aG}dy!ajul%=%YZXe7a=QuES9=K34LZr24Ua#va&g
zwPhV)l~{fzXq`=PGe5*Me=`t^P)!j@6T^2-;X%-z6N02NqjMfMtqP8)7#t@O;_+=c
zoPkd6FA%2agc}2PnDA)gSz_oPO(PNSv$By|GrTZCuh=d)?8Ge6dd^#nF~RUOj_X%_
zfbDDsi)eq7JCGn>)sUEOo=!|JQ6?eW9g5^tPU_b{Guwi&bp{O&UN+=n5;pSh^$#X<
zY1-amwS}7j9JK==1dl>=0%lvenh3$Qf}`B5Yoq9HheI5mLeWR(tR6pEIy6PuVY}KC
zOtBJ5Lq~3;r1)=fC0)EFy@2)!Ej#*fdp$+~?uemejVn^zQ>4_{ONR=r)5C;OqFW6>
zoFEMXHrP0ldI=@Kc>ImV!=*YmH#+Wd)qR;?>Tpg1Ev}qq#NZ#VOx7Zbx<F7S{U3#v
zO;q`jTcI+vB4qq^j=#OCF=~y{gDNm}poS~+7*K*$#om3f3^}LRoNN~Bh1K=KBl;*C
z1urM4jbCi4qh^H+dZIfY+y4jVZ_Cz7H(+092)Kw^Im>0NY};Evc}869tjCyU8*?Y$
zfP4Ur%hbOy&kSUxjH?zIyxsm~T?}PblP6uC`_|N7-9)NCFlFwnCvVnH{E|!;<0T$z
zwvaz(+0|A_1<>lNa~-)OY8Pg|Jtim`5Bd8+j@&Gi69ykrsF0W{+_-F{OvRWfMrEcv
z#;z2tckzsQz5LJCfvdsL#(vly6J<AZpP!HHglo1#>9tA}+thLtz6QROa$_%ij$3P^
ztDrjM6fj+3b3pLJNcH8+wmn>lS7kQeW?PyuVim{jp4t^;tfmxx@kh#y7ge|+ytd@0
zLY<m=omnuUi`?Zl#}{tH*xs%LxJGhKRn)W5_6B{X{HyQ?g%`#lQJW<dm){f6RAUf@
z9%~V?irCC4f$G7u9&&!-RX1ujNF%$L9m*#Fc}6zREgTRw)OsA^KvkfKIj0;N)zSVc
z%4aFvyu1K|&@!m|C0#SSXnPDYVexAkw^5gV|C5@B7xH^Du4Zp-1mF4ud@x542jqd1
z*2vdd<ZX&2uxM9ar<(z;Vy0Dk`@LVJF7Y>7)n4d3e3oQabqXmWZ&XClk|Bq;8O}5t
zM_f}$WKc<aB_MgWZ1<M}6FYb#{Rr5rYE%e@E00=JTN86G^|~@L17_K}=eN<n?ROq>
zgZeahEVaA72Gmi|Z>Mjb<;z@CqOGc5VLmHk6h*ri=`3Li_H-a=q8vxZ3^c7Q3UBBQ
z7q-`lJYl%Jv9BSI!WkY>CB1kV*7={U(CW5q@LY{|73O3d>jQPf-`zXI^u}ko>DV7l
z(6QKR<@DP)PMgLgp9PF@^rEFPn_U;t4?8DJ?|8JiXcWEh)j*GmYWRb@Tf592$~~P1
zi=R68q6yw9WX#y-0{I%1{{*OSTpBj`Cjf<dsuC-eQrdaqUH+Cs4@67%z>88==-mF;
z8sX#7&@;@A{;epmHdt|}jHcYDJiC{QV^6`d!wa*k8bme!L=rP=FQx3>f)>1Nk<NPl
znTbqn?Ndx4$QSc+W<a_JNl$MAvL>%OYee61%+l|q^Qi`lcG8DT*<br&Zy#|+J|0f3
zvXUdmGU9?q;hkmI9Acybh3nPWZI|bpn2We>(`V?%4Rfm3t(Z(SSr^CVm|Ec0=ebzz
z%!AaQ;{o(xHq38!FtFw~`Nk1&)PJvKON)`9m3_p?Hf_tm8Yo|10A(&(0<>aJ22yt^
zn(|FhNg;sjjT$TnPe+gSnWcw_U**ti&4SZTIN6H5$6!@pM0p($h&B;TN47Ee!1cSO
zc>#IkM>{4zfL*C(lGs-JhDP5wCo*Av_V-_4Bf$ECN=7fN1^m}ybI{Z!LN6G|{*Gb(
zMN!ohafE4Y=mzca9apily>pU46YMYPaL&yWi>aKib>HwM<loZ_yoGU3&AxA8KeDgb
zW0MI?7o}%J9-Gf!CsgR7pjlq(ZIoyl>z(12J&b?^{Ui}zkFDD82H3f)bm0+y{!t$S
zy6?!!#(u|D;Kc5cURi?iR#AYnOr;_)<uIgX8a<4bKqBlPnsqwsy_ylB#W|zYX04dz
z9+sYmiZPWm=><Q#FSawJ*%f4XRQFD?TtpO3w;uGN=yD$PU_SElC^)$;nifF>P(%?|
z(utSNkxnWNbUI(JN9i?WTj3Ha^G*A|GYZ|i`qPKwREUIRZX9K!@Je%b0L{7MDFF|0
zA&^#_*X@M?ec_E@9>g_*x8=;;!<LxaJnyh0i;N?82_(!83*e@b=lA|cfK}2^7cl|h
zeG>5GE0}CQ3qMCGzP~)WjwBLPNlPKmb_L(H?3cYrD0L~RPdu}%*3_b*Z;L6bUDE>$
zX&f62e%}$d|8O)SRJAYyfyc2u?(+gYpkN_GjyF&E-);ePIA%Fdsop1Tac^DRipgF<
zV3&kW5-e<eJmx}KUaZxw$S%FtKYB26UwlfO4!*>7xVZoCyBHIhz*w%~t4U8{jCct{
zAJW#^7Uz3|qNk-e5RKvg!02p`m)I8eXsE2*CtFg#Ax)_D5Pgn$3e~AuL`I}_FI>Pd
z+V&W$?8YA&Re~*%Xnhn6grP%w$(l%xhr2R0zzD@W;fok+0<W&Y#p@C_F12F>q3J^l
zE-@L%bUOk2(0*xQS^JMFw)NxaWYk6WHA_y7oeL#Z7#qqcd^@c&i!ZptChjjQ4=7}e
zD((IWJ0eyKqOG6!t7-uqo^4mU>4P0ywE$0yl9WozSs?(L#-iPY>{&e*qMq-I0>PAU
zMmq>fhZFe5Bh%Xj2rAhK3<Dl!*z5f6yy8eN3r(>%GQj{+L7w>Aeu&JV4WBbwBk<ss
z&?|>-t-nbJrVN7l2<o<OS9RVUG-HiN&0bp+f_RK1NXukfC_mj2#Y7OiBq}h_hx}+&
zL+A1kIVeehc5Ge~pxK@N2}b0*M>&_0A4-`@wMa^Mx*<osIx+CwS0JS0Q;d8Or1(|=
zqW&fE2w|>ef`^kh#jA744;}^mjx-T8nS_5rC7h2PG<YYn;+f+!`$E-5;o$hNUdO0t
zHJDM>TzB%oZDYnDBMtT?{`+KvFOlF*Zt8k6tDyi|5H>Htb}*xwwbq=P&=)ap>r+_|
zi}-M=b9}0L|Ht4f5qCWH(BM5=s;$^RSh{e7qqEOz4a3bU2xJttZ3!I~JgCk+o7s@}
zyG|Yv?S>XACbBA=$PQi;r=x5=Y*8hvs%~{e8X^TaVSYf|EV{@EEu4xtnHB^4xwUMb
z$5o7M`}z&0XPo^9zs^6`4vodKSks4pjoR5L3=K(Pw!|ojj-kG>f-3WPz-90exvOCl
zt+y&}l4iSZv^&RC1J~(Rx?2U!RB_|8UY||;sO<4hR~UQrSS}^twjIT|OarllQhAC<
z>E}6P*VVZ^H9dL~?4!yCjM-2Y^NM;99El9bn99M4T3g^vm^a@_K}o|DQKy?peRm^~
zlqU^#`c)9me|U5i_hMp4)MS=9cnq}v=hRKLz~KdC;>HTf1v~B(oZ=S_VwYm;UiTSj
zIuOV~+m2@KY0Y`zIgjFmvsV{W1;yM@GoXO(-Vn$=&5;)nS)X3yHonI&+Zq`Bp7UDM
zNGZosIr#|S8{w-obOe$jpSB!7!-t4bK|ab@afik(0a+&(3q1*Y8VW{Qf>Z~AG2SZ1
z`IVdk{=qUJ3}(t&_2Wl?&D9Q);oG>Alfh~crEWM?)b*k(i!}^WL}~T~hc6!QI5W~5
ziy390*eF)nf7M<h13@W0nD`fKs+}tj*x6s>=;O@ohpsq=g3y>K*~E~2f5+p!SVGg<
z{|*N04Hxna7LxD*nT!PHYk%LltI$2c`99(Os5mo#B5EKRPaH+AFk$c-gV{YxbTiSn
zuTSiF5GTJ(ZRcJ^uIJ7Zr~F-^%Z7&YT5S@@|17IYz2<{TX#XJvsxyDSR!cH0&dp82
zE^z*Y&NMq8ZjfXYf#xsSh%Pfr$B$8-2CmmahS;xfC^Zg5d+I4BZ~+!hZ&&@4-yZyn
z8=>8b9<Dku64f#&1M{w!FKxromC=+5EcmIfWF{cYEv<)6J%>{$!(9vGnwa)(g6<Zi
zfij%|#8Sn}Qkjd3zgO}oO6fqY9;hb0xoe!7WZC8BkB;>DTsVaR{Ko{rRsvjD^81KF
zGOgHA6@j`7T~J+3vjg0;+yiWO0hWX|AQ#|It!O}x<hccIdECmzF-eaYY~RF+c_@u2
z#Pj!bijf!5i;6AxwtH7V=mqQP0Uh_JwL82_VA<8UmOd#6l4iTUn*Z&~VGaa1MUs6D
zZD&ex4&G*$I7F&~(MznBh~h6Tc74VfKnUK*et<}4(Z|%X%|}k)dce5fZKa+Bi#<$Z
z`jN6FB7R`}!9Xa>_t_|b>^NOuzE(lfAl-&Tk3D1HV)$y0f9oc*Y;so$N!9_O=ROV-
z(UbJdTE>kxaicPYq;D5{=&LkG5Q)unv!8mUr}NO~1Y#-EJlIp1Z7U5EGUQ<U=Id7c
zw~x?MzD)gs!8>5aAi8&}Ia1r9k0b75-cBupdetpZjr>rs4D3_d(nhFL5z$`r_VR2v
z7^_S4e$wTl^9~e&-F%kSb*v;43XND&>8xM@j>dbl7dXMJnYy6_#`?8Z?=8X?aR@5)
zl9ce<Qgtsy64I87?lWWi#h`oZ`s(Rq<@*<9uer8!d>>skJzYd+S*7X8Yf*etHFhi9
zdgLG5yy7bD4I2q}L+i6<#L5i+{wK-~%qKZRb;KAsXm`i~(%h%H&83+Olx?wmLlULn
z6C_S<)M0D|6yW6nf)46>p``NXE6M~W`|DNx8=R$RTCtQq`aGHslMuT1&5Zg5$W%vC
zft-<y98dC878zwIpBzHXQE;k2GcgCOIy-~%W2PbLwWMfy`gn5(n4oii4is*`T0=@c
zmv8P!1jumd&M7}^O$W;6Yu~;a8oMZM))>sHH3p3bzHtX&Mc(b%&iyO|*8!M_!vQnR
zY&-%e)C17>M36@Dw?@T_fymU&I6gGN0b>^wKwwd42lusYG^n*(S|LkQnlB>$xzN~W
zPsyqc^wI%L-EbEVJafYOL_KsU{}B%l)^c(qL*E`frqc!C9RAzH0%Hx$Mu<1^+QyR_
z-c#}=NZZxRxUI$eheK^Sze(3iScl`8qI>*a!T}bmS;+4rY*Jg(p$&T7$lv(85ZsQ8
ztfjh6NqyPpzBbEiulyOr)()nOm4x>@IF3=b<BjN}jmZ6H{_P{R+V8Enxgw~!>p&Gu
zojKyM{;o3m3NZNT(@4uUk@TZ3{t054Veo5cV1MfP%oZ6xjdd~IjDoZqqar@aMkwTw
zX}UdB6WgjP02-L&BV|D`YoL|VnE(U!t1jGe`j$_Qai-v7=|F{gwlYiVw8j`#pQTy0
zEMW34&;rd!j+*}x`9=nXa7)8e1fnif=YWrZ$2^<Hw<?)~v7GlJD$jlzIWZ3v^WfpR
z7e`x2Fe#%YDVCgGrFOS*&mo?czXkw;*^qdyLc(Wfx7d8x>M{qHJS^k<kvMCyBPoN!
zbb^lI;&|G#Ocde338?eF97C<At?Z)0f|GF5?v3co^fAYcpS)dILss@^je+$x$2L2t
zMY1NQ+m6)t#adN)N3Od9Ufw;Z*iB;;$pTaP5*5|dTv(|gK*9T$RWhT$%TGitc^3TY
zkkCFjc9{!|ZIkZ_>K{I;*P(eInqsIsKSXA2#SJ$+p8XC4#w{*aU+KO8i*z>M9yi2q
zZJ^lBP25(i!Oh<|^3A_aF3f?OcD7%GgWPjJ)i+MMY-4p}>N!I#8qI2AFGi;9$uo}|
zJK*Fqh>efmcxWoS@;8zc-NVyRgYqZuL$~c+7qY*dpVBGFB(cc^k&Fd?YI5Q7o)L%a
z%s=XZHfZxxFzBXvs*!ETx|^fdt=^Dt3_LPMB!7=Wbs#2~T2M&hN~P1}5v|fy`M;kX
z%^0}{x${mBu^M&OS%mA+Z6As`#DMpb)oz_@xXh)CYtsFO-o>$Cr5dg-NJ59=O@Y0Z
z1S{nLO6dBh=3evWPJjL}LD;?Ae|53IIl~J}B8(9TJ|qzYh;Bs*DPzD-sx|Ffp)bVV
zJCSDSJTHTBzZdyBJfWb~?Y~>b%07h#<!fVI)w+v120C4T2g5Q8L*v&8bZ=l}`AOLe
z?8Ev~dI~vTc4tDh>59{%{y?vK?%oY$3O9)tPB=D8pF!xz)#n>3aNQr)#0Swl1L?PK
z6ktP~kqLtD0=B6r!*OD<m)f)yWt)v><YAl)hVK-<9XCybHme>x#mB9AK}fdke~&j~
zbc7OHkeHVqKcXff!kHK)GJ(BN&W;wAyH5xSls{n4{h=N5M!nDv2_Rn%^7(13vHf6!
zgeEh4mP<qY(8DxUWh2-(gAZVtlcL1Ptlo^M7hjbGM1G7rb%#pQy@3vRZ1(Mo4oy15
zS~ObI2q?%Yx1l|zP$PKZRkC`xl}3hIpHFD*KAFi4=w8R2M^4uVkl~_^7$CK8JYkf}
z8G&N7jeSud&Qb^O45J%-PPg!6+5HfUVE%&B)|VB!fHbZ$6cfn?@kvr*FPr3Dim%+x
zN*~^QTn|K!6%q-ZL!FeJ5p$|fDe?0}%7fGEDV1oWc4F0v1xUHZlK`hB&5Q7A6~#7a
zI)Jwlnnf6JlU<n$5wgh&U;tXvb?4Z|r11sP^5m(Gw5SGNygE5~!u*RVdaYRw1LIOw
zM)c!CPQ<~n=4{K5ngL3f(*M^T`ukpvz|EzRQgQmAco)fS6jK8|u8A#oB|?7$os&va
z;dM|ZLP<zy-n{U0V+=V+o|vnmd;)_H(}yEA-n(X5Hx!A3MDR!Uk4RNq0=Kuwd_`{+
znzVyAZAFqNJFHq&48!!EBfYzTVLbc1YMCws+gPQKe`z)z6Q~aRb#vwP>qNG3O7pIF
z|MRO|0}Pvt#E=@}`W?qI7MZ_9gXbg$OBQ$|pB#i`y<LOhU{m^FiY=xd#w0-p^TglC
zS{_=&))nb>6F)Fd1l8&4V9XipRd7Hs{KM5*h%@-Mf~YZZCu+3?6kAm8#-MsKdT)7v
z&iVV6O2~~|-|#iI&~3Ewh5ikQ71iOGYxt@B_611;4aL7He<W;q>g7HrQTdKG-hu}-
zjtGb#U=Nw=#SEONHIT2rQ{$M>>(3c&c7#abEBX$5HWrNsk*siM?r}pGW@z8*QuPU+
zkaZ9&hAtCkhHeq;SzB-XqoF$YY%bagY3Kvj#41g6cJXP9r9?DBK$=2?Yk5>y6c|ci
z-6W3njPaIzC_0GebQXd>6d4vY*J<c-M;8=@=u9)@c5U*7{1gWEM~^Xqxp|1f=4~=O
z0yR9a+XFk_)IRfSiF1krRNl!hQz(Fydk^H;`O3*DXa4J^md#I#rroztqRh5_eYlBq
z)yKt>hv<cf4J<fSp9jKw&A9O;D|^3}ERr{eMkAFj6c&d6b63aI=tfw*qSY2xk!`MR
z$A?J|p*-J5*f9J}I#j`oC!KR=Fd%9O@bBCq_Hei!wPZWN#nRscYi;}(F9MBvc+h_)
zgR=Ym>D_4TuhDhKnp49HkuPsZspL&qEgtpGYdYe71dO_Ac?LJC_w^aE(*rCa5V$Lo
zzBx;5UsssrfP|VD;L*$a!*;0e@m9^QUo$?3y`=ql?@`=Gd$dWoD1mLug??Tro+IV<
zb>BB=wOu0l04=to*c9z8=P2F1oK;!9TKO@X;Pp{=AJ}kUG8y5c>~(e8ihG*zITOVe
z>Acz9on6GCC3+NARFEK$i0cz%MQ!g^BnHwZvB|jw^>v<uj*-9(z^QcDlXcw!xcnaC
zPh;;p!G#FGg{MR?cL2O;S<WJ*28h4A4?QBb>!4F}nqbmrHQe310*s{EHSLFE26AUt
zuurELT6)VVpDOAAs5@GBT`_U|58JUqS0+;l`j_#i&8AN=zZG+Hc@<DHVmXQ}z&>)8
z{`(F|0V&8Dt9^QV+}y8ea!EK&^=a`ruxWbVpz{~Nvcg2?@wJ(WkGs^h&;p1lCq)d)
zoqyuP4I!E*+uJ27=m=+1KeN3nidMowa9^AvI(ocK5aY5tA6Sw6?jZBxw%h0L<`9Wz
zk$;UGo<ExGgdyRsDrVTt@o2HuUDyI<ZH)9O%2bshjHJ`QD%KKg{KnILp%;wN;Y7-0
z0aPt1c9<f&W!fSo@0c0O8@kPm!lz2JC*KS=I{+thW*O&%B3G`v!}~J7#<pL(DakWj
z9&;GeVTR&kE*A2Bi+u`jNw_c6HB;K`o)R&tC<({{<PgH1yy%C_+VIU(B<4e7nsd`(
zZgb#`+8Jm(>FTt$AsKs^r0Gja2*5w3I80m2ANWb{K6)zVq24-|Wpb282r6L}yJR0<
zIvVkGk<2Q<XFDcqg6aA!%xNIopiMc`4z;(<Ia$;4ph@Inwi_}A9mAI&(tlvzAi;1%
z-ZqD~nWUMGPG&a3|8O=<FvOm>Z{2b1`K3gR41p>{VkT4CJ>Lw0y|kyO>r|)b3e#zM
z07pi4KIv}Qo5|6f4eBlj*uJDe>m(0pzb7<tj>I4h8d>`b(O5(*h%?AQf2T|H#gFK#
zGr5N+;C9R)pA5|vE{l4Dl;nY6pVa3MLPtRRPBV8{InEn+`SrXGLZ_d?vnpe!&56mX
zr}c_;n_!m}Yv_48Cb2-tRgEQip~``CDas{xf*8w13A=+(t87sT`I00MMGv}bEf1+K
z;3Lq=V!rnX&swV|6pb$|#X!xPo)u$~kVB+*0o!4T&0zt~?-Q`!!W9wvemq*XE*nmr
zF4NhjUGc^MeuQu^TNDhP1Ei)>%R65H>9Rb`ZeW_LhP<kk{@%%%7D&bE9n5VWE_G|P
zTs!T=B%qQm=6*nvE)*UydQm~0u}`8O>Y7Ghab>BL3>99dx2wjr(7(c!o$np;gUZ_d
z5(gI41aP^BZ=&eNbt;{AC>~&&LRtIse@UA~h<b#l*r{B`2!1ViHWbg4a1Siq^5@z(
zbsh{rQiSi6gV>A;!s-gMgWaf|McaY&C5x%elPYMO#ja%Ih;3iZ$6QKgc8=4P5dFZx
z8Cl(@Q9_H@tIhsQ0c9EPm|@vH=^?uilYtnhE+4qNYx6S}YlO|HSJ53D6{0#Ca2o>$
zXh)Wlt~HS49MrU+Wh(eyn^e$>8(zx&lT6aqk4jLXi&+hd0-H8J3yuFz!m>%Fw%Iw?
OLCor@UfFK2rUus?Cy@^T
literal 0
HcmV?d00001
diff --git a/ci/run.ci.sh b/ci/run.ci.sh
new file mode 100755
index 0000000..8c2d99b
--- /dev/null
+++ b/ci/run.ci.sh
@@ -0,0 +1,27 @@
+#!/bin/bash
+
+rm -f *.db *.log *.txt
+
+export OPENMM_CPU_THREADS=1
+export OPENMM_PLUGIN_DIR=$PWD/../build/openmm-7.3.0/lib/plugins
+
+mpirun -np 8 $PWD/../build/parRep -i input_generalized_parRep.lua -log warn -o out.txt -e err.txt --inp-seeds ./ref.sim/seeds.bin
+
+MPI_RET_CODE=$?
+
+# compare the db file from the run and the reference, they should be the same
+sqldiff lj7.db ref.sim/lj7.db > diff.txt
+
+SQL_RET_CODE=$?
+
+if [ -s "diff.txt" ]
+then
+ TEST_RET_CODE=-1
+else
+ TEST_RET_CODE=0
+fi
+
+cat out.*.txt
+
+# success only if there was no error code from any of the previous commands
+exit $(($MPI_RET_CODE|$SQL_RET_CODE|$TEST_RET_CODE))
diff --git a/include/dyna/GelmanRubin.hpp b/include/dyna/GelmanRubin.hpp
index fefc5cf..b33c672 100644
--- a/include/dyna/GelmanRubin.hpp
+++ b/include/dyna/GelmanRubin.hpp
@@ -4,7 +4,7 @@
* \author Florent Hédin
* \author Tony Lelièvre
* \author École des Ponts - ParisTech
- * \date 2016-2018
+ * \date 2016-2019
*/
#ifndef GELMANRUBIN_HPP_INCLUDED
@@ -34,13 +34,11 @@ public:
/**
* @brief Constructor for a Gelman Rubin Analysis object; requires at least the number of chains (replica) as first argument.
- * The second argument allows to discard the first n observations when averaging, default is zero and thus disables
- * discarding.
*
* @param _num_chains The number of G-R chains (same as the number of replicas usually).
- * @param _discard Discards the first n observations, as usually there won't be convergence at the beginning.
+ * @param _no_convergence_if_less_than Optional argument : if set there won't be convergence check until at least _no_convergence_if_less_than observations have been accumulated for each observable : this may hel avoiding pseudo-convergence.
*/
- GelmanRubinAnalysis(const uint32_t _num_chains, const uint32_t _discard=0);
+ GelmanRubinAnalysis(const uint32_t _num_chains, const uint32_t _no_convergence_if_less_than=0);
/**
@@ -195,8 +193,9 @@ private:
std::map<std::string,double> tolerance; ///< Allows one tolerance treshold per observable type
std::map<std::string,std::vector<double>> ratio; ///< Allows one ratio vector per observable type
- uint32_t num_chains; ///< The total number of chains
- uint32_t discard_first; ///< A given number of records to ignore at the beginning
+ uint32_t num_chains = 0; ///< The total number of chains
+ uint32_t min_num_observations_before_check = 0; ///< There won't be convergence check until at least this amount ob observations have been accumulated for each obseravble
+
};
diff --git a/include/dyna/observable.hpp b/include/dyna/observable.hpp
index 113b359..e01b70d 100644
--- a/include/dyna/observable.hpp
+++ b/include/dyna/observable.hpp
@@ -4,7 +4,7 @@
* \author Florent Hédin
* \author Tony Lelièvre
* \author École des Ponts - ParisTech
- * \date 2016-2018
+ * \date 2016-2019
*/
#ifndef OBSERVABLE_HPP_INCLUDED
diff --git a/include/dyna/parRep.hpp b/include/dyna/parRep.hpp
index 9e6efdd..56ca276 100644
--- a/include/dyna/parRep.hpp
+++ b/include/dyna/parRep.hpp
@@ -4,7 +4,7 @@
* \author Florent Hédin
* \author Tony Lelièvre
* \author École des Ponts - ParisTech
- * \date 2016-2018
+ * \date 2016-2019
*/
#ifndef PARREP_HPP_INCLUDED
diff --git a/include/dyna/parRep_FV.hpp b/include/dyna/parRep_FV.hpp
index 0a40dc9..3d610ec 100644
--- a/include/dyna/parRep_FV.hpp
+++ b/include/dyna/parRep_FV.hpp
@@ -4,7 +4,7 @@
* \author Florent Hédin
* \author Tony Lelièvre
* \author École des Ponts - ParisTech
- * \date 2016-2018
+ * \date 2016-2019
*/
#ifndef PARREP_FV_HPP_INCLUDED
diff --git a/include/dyna/parRep_FV_multiExits.hpp b/include/dyna/parRep_FV_multiExits.hpp
index 6006d44..05a6d08 100644
--- a/include/dyna/parRep_FV_multiExits.hpp
+++ b/include/dyna/parRep_FV_multiExits.hpp
@@ -4,7 +4,7 @@
* \author Florent Hédin
* \author Tony Lelièvre
* \author École des Ponts - ParisTech
- * \date 2016-2018
+ * \date 2016-2019
*/
#ifndef PARREP_FV_MULTI_HPP_INCLUDED
diff --git a/include/dyna/parRep_abstract.hpp b/include/dyna/parRep_abstract.hpp
index 7d698b5..85a1aad 100644
--- a/include/dyna/parRep_abstract.hpp
+++ b/include/dyna/parRep_abstract.hpp
@@ -4,7 +4,7 @@
* \author Florent Hédin
* \author Tony Lelièvre
* \author École des Ponts - ParisTech
- * \date 2016-2018
+ * \date 2016-2019
*/
#ifndef PARREP_ABSTRACT_HPP_INCLUDED
diff --git a/include/dyna/runSim.hpp b/include/dyna/runSim.hpp
index c0b144e..bfee426 100644
--- a/include/dyna/runSim.hpp
+++ b/include/dyna/runSim.hpp
@@ -4,7 +4,7 @@
* \author Florent Hédin
* \author Tony Lelièvre
* \author École des Ponts - ParisTech
- * \date 2016-2018
+ * \date 2016-2019
*/
#ifndef RUNSIM_HPP_INCLUDED
@@ -20,8 +20,15 @@
* + Parsing Input file
* + Setting up the desired type of parRep run_simulation
* + Running the simulation
+ *
+ * @param inpf The path to the input Lua file, mandatory.
+ * @param seeds_inp_file The path to the unique binary file containing random seeds to be loaded, optional.
+ * @param seeds_out_file The path to the unique binary file containing random seeds to be loaded, optional.
*/
-void run_simulation(const std::string& inpf);
+void run_simulation(const std::string& inpf,
+ const std::string& seeds_inp_file = NULLFILE,
+ const std::string& seeds_out_file = NULLFILE
+ );
#endif // RUNSIM_HPP_INCLUDED
diff --git a/include/global.hpp b/include/global.hpp
index 78ece31..279d41c 100644
--- a/include/global.hpp
+++ b/include/global.hpp
@@ -4,7 +4,7 @@
* \author Florent Hédin
* \author Tony Lelièvre
* \author École des Ponts - ParisTech
- * \date 2016-2018
+ * \date 2016-2019
*/
#ifndef GLOBAL_HPP_INCLUDED
diff --git a/include/logger.hpp b/include/logger.hpp
index 0af6c4a..7b2d6ea 100644
--- a/include/logger.hpp
+++ b/include/logger.hpp
@@ -4,7 +4,7 @@
* \author Florent Hédin
* \author Tony Lelièvre
* \author École des Ponts - ParisTech
- * \date 2016-2018
+ * \date 2016-2019
*/
#ifndef LOGGER_H_INCLUDED
diff --git a/include/md/md_interface.hpp b/include/md/md_interface.hpp
index 7da5996..d8a0f93 100644
--- a/include/md/md_interface.hpp
+++ b/include/md/md_interface.hpp
@@ -4,7 +4,7 @@
* \author Florent Hédin
* \author Tony Lelièvre
* \author École des Ponts - ParisTech
- * \date 2016-2018
+ * \date 2016-2019
*/
#ifndef MD_INTERFACE_H
diff --git a/include/md/omm_interface.hpp b/include/md/omm_interface.hpp
index c3acb1d..29ec520 100644
--- a/include/md/omm_interface.hpp
+++ b/include/md/omm_interface.hpp
@@ -4,7 +4,7 @@
* \author Florent Hédin
* \author Tony Lelièvre
* \author École des Ponts - ParisTech
- * \date 2016-2018
+ * \date 2016-2019
*/
#ifndef OMM_INTERFACE_H
diff --git a/include/rand.hpp b/include/rand.hpp
index ad8b181..9f9ebac 100644
--- a/include/rand.hpp
+++ b/include/rand.hpp
@@ -4,31 +4,45 @@
* \author Florent Hédin
* \author Tony Lelièvre
* \author École des Ponts - ParisTech
- * \date 2016-2018
+ * \date 2016-2019
*/
#ifndef RAND_HPP_INCLUDED
#define RAND_HPP_INCLUDED
-#include <cstdint>
-
+#include "global.hpp"
#include <string>
/**
- * @brief Call this function for initialising a Mersenne twister 19937 random numbers generator, using random seeds.
- *
- * The generator is seeded by using integers taken from a std::random_device (usually by reading /dev/urandom where available, or using the hardware instruction RDRND on compatible CPUs)
- */
-void init_rand();
+* @brief Used to indicate if the random seeds used for initialising the random number generators should be read from or written to a file (for reproducibility).
+* See init_rand(...) which uses it.
+*/
+enum SEEDS_IO
+{
+ NONE = 0, ///< No I/O of the seeds
+ SAVE_TO_FILE = 1, ///< seeds from all ranks will be saved to a unique binary file
+ LOAD_FROM_FILE = 2, ///< seeds for all ranks will be read from a unique binary file
+};
/**
- * @brief Call this function for initialising a Mersenne twister 19937 random numbers generator, using as seeds a list of unsigned 64 bits integers (encoded as a comma separated string, at least 5 seeds required).
+ * @brief Call this function for initialising a Mersenne twister 19937 random numbers generator, using random seeds.
+ *
+ * The generator is by default seeded using integers taken from a std::random_device (usually by reading /dev/urandom where available, or using the hardware instruction RDRND on compatible CPUs).
+ * Furthermore the seeds obtained from a std::random_device can be saved to a file for later re-use.
+ *
+ * If you load a seed file generated by this function during a previous run, you will get the same initialization and should be able to reproduce results (at least if runing on the same architecture and using the same compiler...).
+ *
+ * The file should contain exactly NRANKS*512 unsigned 32 bits integers (uint32_t), where NRANKS is the number of MPI ranks running, i.e. each ranks requires 512 uint32_t for seeding
+ * its Mersenne twister 19937 engine.
+ *
+ * If you want to visualize as text the binary seed file "seeds.bin", use the following unix command : "od -A n -t u4 seeds.bin"
*
- * @attention it is crucial to provide a different string on each MPI rank, otherwise the generators will be correlated/identical !!!!
*
- * @param seeds seeds is a std::string containing n seeds, comma separated (e.g. seed = "123,456,789"). n is >= 5
+ * @param io_type Set to either SEEDS_IO::NONE, SEEDS_IO::LOAD_FROM_FILE or SEEDS_IO::SAVE_TO_FILE. ; the first only generates seeds for internal use, the second will load seeds from file seeds_file_name, and the third will save the seeds internally generated to file seeds_file_name.
+ * @param seeds_file_name if io_type is not SEEDS_IO::NONE, this should contain a valide path/filename from which to load the seeds or where to write the seeds. Defaults to NULLFILE.
*/
-void init_rand(const std::string& seeds);
+void init_rand(SEEDS_IO io_type,
+ const std::string& seeds_file_name = NULLFILE);
/**
* @brief Call this function for obtaining a uniformly distributed random number in the range [0.0,1.0]
diff --git a/include/utils/lua_interface.hpp b/include/utils/lua_interface.hpp
index d35feed..0616bb2 100644
--- a/include/utils/lua_interface.hpp
+++ b/include/utils/lua_interface.hpp
@@ -4,7 +4,7 @@
* \author Florent Hédin
* \author Tony Lelièvre
* \author École des Ponts - ParisTech
- * \date 2016-2018
+ * \date 2016-2019
*/
#ifndef LUA_INTERFACE_HPP_INCLUDED
diff --git a/include/utils/mpi_utils.hpp b/include/utils/mpi_utils.hpp
index 3f22630..c3f45e5 100644
--- a/include/utils/mpi_utils.hpp
+++ b/include/utils/mpi_utils.hpp
@@ -4,7 +4,7 @@
* \author Florent Hédin
* \author Tony Lelièvre
* \author École des Ponts - ParisTech
- * \date 2016-2018
+ * \date 2016-2019
*/
#ifndef MPI_UTILS_H_INCLUDED
diff --git a/mol/1d7h_i/input_generalized_parRep.lua b/mol/1d7h_i/input_generalized_parRep.lua
index 1dcbeda..25b1e50 100755
--- a/mol/1d7h_i/input_generalized_parRep.lua
+++ b/mol/1d7h_i/input_generalized_parRep.lua
@@ -116,6 +116,11 @@ simulation =
checkFV = 2000, -- 2 ps
-- a frequency (in steps) at which to evaluate convergence using the Gelman-Rubin statistics
checkGR = 100, -- 100 fs
+ -- check FV convergence if at least minAccumulatedObs have already been accumulated
+ -- This may be useful if there is a risk of early pseudo-convergence for some of the observables when only a few samples have been accumulated
+ -- Therefore the minimum FV convergence time will be minAccumulatedObs*checkGR*dt (here 100*100*dt = 20000 fs = 20 ps, but it is likely that the convergence time will be larger than that)
+ -- Default is 100
+ minAccumulatedObs = 100,
-- parameter for parallel dynamics stage
-- checkDynamics is the frequency (in steps) at which to verify if the system left the current state
checkDynamics = 2000, -- 2 ps
@@ -155,7 +160,6 @@ simulation =
-- implicit functions :
-- ------------------
--
---
-- exit_from_lua() : call this if it is required to finish the simulation from the lua script : it will terminate MPI properly
-- and flush I/O files ; but it won't perform a DB backup so do it manually before.
--
@@ -190,13 +194,13 @@ simulation =
-- get_COM() : function returning the center of mass of the system as a tuple x,y,z
--
-- get_COM_idxs(idxs) : function returning the center of mass of a subset of the system as a tuple x,y,z
--- NOTE this time idxs is indexed in C++ style
--- for example get_COM_idxs({0,1,2}) to get COM of atoms 1, 2 and 3 (C++ : atoms 0, 1, 2)
+-- NOTE this time idxs is indexed directly in C++ style
+-- for example get_COM_idxs({1,2,3}) to get COM of atoms 1, 2 and 3 (C++ : atoms 0, 1, 2)
--
-- get_minimised_energy(tolerance,maxSteps) : this function returns the minimised energy of the system, using the OpenMM L-BFGS minimiser
-- note that coordinates are not affected, it just returns the minimum epot of the bassin in which dynamics currently evolves
-- it returns a 3-tuple ep,ek,et (potential, kinetic and total energy)
--- the tolerance and maxSteps can be the above defined simulation.minimisationTolerance and simulation.minimisationMaxSteps
+-- the tolerance and maxSteps can be the above defined minimisation.Tolerance and minimisation.MaxSteps
--
-- get_minimised_crdvels(tolerance,maxSteps) : this function returns a 2-tuple (crds,vels) containing
-- a copy of coordinates and velocities after minimisation.
@@ -205,6 +209,9 @@ simulation =
-- the tolerance and maxSteps can be the above defined simulation.minimisationTolerance and simulation.minimisationMaxSteps
-- NOTE lua indexing, starting from 1
--
+-- get_minimised_energy_crdvels(tolerance,maxSteps) : this returns a 5-tuple (ep,ek,et,crds,vels)
+-- This does the same as the two previous functions but with only one call
+--
-- hr_timer() : returns a variable representing a c++ high precision timer : can be used for measuring execution time.
-- do not try to modify it or even read it, it should only be used as argument for the following hr_timediff_* functions.
--
diff --git a/mol/ala2vac/input_generalized_parRep.lua b/mol/ala2vac/input_generalized_parRep.lua
index fab791f..6e181b1 100755
--- a/mol/ala2vac/input_generalized_parRep.lua
+++ b/mol/ala2vac/input_generalized_parRep.lua
@@ -108,6 +108,11 @@ simulation =
checkFV = 250, -- 0.5 ps
-- a frequency (in steps) at which to evaluate convergence using the Gelman-Rubin statistics
checkGR = 10, -- 20 fs
+ -- check FV convergence if at least minAccumulatedObs have already been accumulated
+ -- This may be useful if there is a risk of early pseudo-convergence for some of the observables when only a few samples have been accumulated
+ -- Therefore the minimum FV convergence time will be minAccumulatedObs*checkGR*dt (here 1000*10*dt = 20000 fs = 20 ps, but it is likely that the convergence time will be larger than that)
+ -- Default is 100
+ minAccumulatedObs = 1000,
-- parameter for parallel dynamics stage
-- checkDynamics is the frequency (in steps) at which to verify if the system left the current state
checkDynamics = 2500, -- 5 ps
@@ -187,7 +192,7 @@ simulation =
-- get_minimised_energy(tolerance,maxSteps) : this function returns the minimised energy of the system, using the OpenMM L-BFGS minimiser
-- note that coordinates are not affected, it just returns the minimum epot of the bassin in which dynamics currently evolves
-- it returns a 3-tuple ep,ek,et (potential, kinetic and total energy)
--- the tolerance and maxSteps can be the above defined simulation.minimisationTolerance and simulation.minimisationMaxSteps
+-- the tolerance and maxSteps can be the above defined minimisation.Tolerance and minimisation.MaxSteps
--
-- get_minimised_crdvels(tolerance,maxSteps) : this function returns a 2-tuple (crds,vels) containing
-- a copy of coordinates and velocities after minimisation.
@@ -196,6 +201,9 @@ simulation =
-- the tolerance and maxSteps can be the above defined simulation.minimisationTolerance and simulation.minimisationMaxSteps
-- NOTE lua indexing, starting from 1
--
+-- get_minimised_energy_crdvels(tolerance,maxSteps) : this returns a 5-tuple (ep,ek,et,crds,vels)
+-- This does the same as the two previous functions but with only one call
+--
-- hr_timer() : returns a variable representing a c++ high precision timer : can be used for measuring execution time.
-- do not try to modify it or even read it, it should only be used as argument for the following hr_timediff_* functions.
--
diff --git a/mol/ala2vac/input_generalized_parRep_multiExit.lua b/mol/ala2vac/input_generalized_parRep_multiExit.lua
index b3ead1b..e2ef33c 100755
--- a/mol/ala2vac/input_generalized_parRep_multiExit.lua
+++ b/mol/ala2vac/input_generalized_parRep_multiExit.lua
@@ -115,6 +115,11 @@ simulation =
checkFV = 250, -- 0.5 ps
-- a frequency (in steps) at which to evaluate convergence using the Gelman-Rubin statistics
checkGR = 10, -- 20 fs
+ -- check FV convergence if at least minAccumulatedObs have already been accumulated
+ -- This may be useful if there is a risk of early pseudo-convergence for some of the observables when only a few samples have been accumulated
+ -- Therefore the minimum FV convergence time will be minAccumulatedObs*checkGR*dt (here 1000*10*dt = 20000 fs = 20 ps, but it is likely that the convergence time will be larger than that)
+ -- Default is 100
+ minAccumulatedObs = 1000,
-- parameter for parallel dynamics stage
-- checkDynamics is the frequency (in steps) at which to verify if the system left the current state
checkDynamics = 2500, -- 5 ps
@@ -194,7 +199,7 @@ simulation =
-- get_minimised_energy(tolerance,maxSteps) : this function returns the minimised energy of the system, using the OpenMM L-BFGS minimiser
-- note that coordinates are not affected, it just returns the minimum epot of the bassin in which dynamics currently evolves
-- it returns a 3-tuple ep,ek,et (potential, kinetic and total energy)
--- the tolerance and maxSteps can be the above defined simulation.minimisationTolerance and simulation.minimisationMaxSteps
+-- the tolerance and maxSteps can be the above defined minimisation.Tolerance and minimisation.MaxSteps
--
-- get_minimised_crdvels(tolerance,maxSteps) : this function returns a 2-tuple (crds,vels) containing
-- a copy of coordinates and velocities after minimisation.
@@ -203,6 +208,9 @@ simulation =
-- the tolerance and maxSteps can be the above defined simulation.minimisationTolerance and simulation.minimisationMaxSteps
-- NOTE lua indexing, starting from 1
--
+-- get_minimised_energy_crdvels(tolerance,maxSteps) : this returns a 5-tuple (ep,ek,et,crds,vels)
+-- This does the same as the two previous functions but with only one call
+--
-- hr_timer() : returns a variable representing a c++ high precision timer : can be used for measuring execution time.
-- do not try to modify it or even read it, it should only be used as argument for the following hr_timediff_* functions.
--
diff --git a/mol/ala2vac/input_parRep.lua b/mol/ala2vac/input_parRep.lua
index 9808848..acb560e 100755
--- a/mol/ala2vac/input_parRep.lua
+++ b/mol/ala2vac/input_parRep.lua
@@ -184,7 +184,7 @@ simulation =
-- get_minimised_energy(tolerance,maxSteps) : this function returns the minimised energy of the system, using the OpenMM L-BFGS minimiser
-- note that coordinates are not affected, it just returns the minimum epot of the bassin in which dynamics currently evolves
-- it returns a 3-tuple ep,ek,et (potential, kinetic and total energy)
--- the tolerance and maxSteps can be the above defined simulation.minimisationTolerance and simulation.minimisationMaxSteps
+-- the tolerance and maxSteps can be the above defined minimisation.Tolerance and minimisation.MaxSteps
--
-- get_minimised_crdvels(tolerance,maxSteps) : this function returns a 2-tuple (crds,vels) containing
-- a copy of coordinates and velocities after minimisation.
@@ -193,6 +193,9 @@ simulation =
-- the tolerance and maxSteps can be the above defined simulation.minimisationTolerance and simulation.minimisationMaxSteps
-- NOTE lua indexing, starting from 1
--
+-- get_minimised_energy_crdvels(tolerance,maxSteps) : this returns a 5-tuple (ep,ek,et,crds,vels)
+-- This does the same as the two previous functions but with only one call
+--
-- hr_timer() : returns a variable representing a c++ high precision timer : can be used for measuring execution time.
-- do not try to modify it or even read it, it should only be used as argument for the following hr_timediff_* functions.
--
diff --git a/mol/ala2vac_transient/input_generalized_parRep.lua b/mol/ala2vac_transient/input_generalized_parRep.lua
index 88d5d57..7475812 100755
--- a/mol/ala2vac_transient/input_generalized_parRep.lua
+++ b/mol/ala2vac_transient/input_generalized_parRep.lua
@@ -27,8 +27,8 @@ MD_Engine = "OpenMM"
-- OCL : on cpu or gpu of any brand, or any accelerating device available
-- CUDA : on nvidia gpu only, usually the fastest
-- OMMplatform = "AUTO"
--- OMMplatform = "REF"
-OMMplatform = "CPU"
+OMMplatform = "REF"
+-- OMMplatform = "CPU"
-- OMMplatform = "OCL"
-- OMMplatform = "CUDA"
@@ -42,7 +42,7 @@ OMMplatform = "CPU"
-- ...
-- CPU platform properties : Threads = "1" is equivalent to defining OPENMM_CPU_THREADS=1 in the environnment
-OMMplatform_properties = { Threads = "1"}
+-- OMMplatform_properties = { Threads = "1"}
-- OpenCL platform properties
-- OMMplatform_properties = {
@@ -64,15 +64,15 @@ OMMplatform_properties = { Threads = "1"}
-- load the integrator parameters from a serialised OpenMM XML file ;
-- no default, error if undefined
-- adapt the path to the file in the following line
-integrator = { xml = "./mol/ala2vac/Integrator.xml" }
+integrator = { xml = "./mol/ala2vac_transient/Integrator.xml" }
-- load the OpenMM System from a serialised XML file
-- no default, error if undefined
-system = { xml = "./mol/ala2vac/System.xml" }
+system = { xml = "./mol/ala2vac_transient/System.xml" }
-- load the OpenMM State from a serialised XML file
-- no default, error if undefined
-state = { xml = "./mol/ala2vac/State.xml" }
+state = { xml = "./mol/ala2vac_transient/State.xml" }
-- parameters for energy minimisation : tolerance and maximum number of steps (if 0 : no limit, continue until tolerance satsfied)
-- defaults are : minimisation={Tolerance=1e-6,MaxSteps=0}
@@ -108,6 +108,11 @@ simulation =
checkFV = 250, -- 0.5 ps
-- a frequency (in steps) at which to evaluate convergence using the Gelman-Rubin statistics
checkGR = 10, -- 20 fs
+ -- check FV convergence if at least minAccumulatedObs have already been accumulated
+ -- This may be useful if there is a risk of early pseudo-convergence for some of the observables when only a few samples have been accumulated
+ -- Therefore the minimum FV convergence time will be minAccumulatedObs*checkGR*dt (here 1000*10*dt = 20000 fs = 20 ps, but it is likely that the convergence time will be larger than that)
+ -- Default is 100
+ minAccumulatedObs = 1000,
-- parameter for parallel dynamics stage
-- checkDynamics is the frequency (in steps) at which to verify if the system left the current state
checkDynamics = 2500, -- 5 ps
@@ -187,7 +192,7 @@ simulation =
-- get_minimised_energy(tolerance,maxSteps) : this function returns the minimised energy of the system, using the OpenMM L-BFGS minimiser
-- note that coordinates are not affected, it just returns the minimum epot of the bassin in which dynamics currently evolves
-- it returns a 3-tuple ep,ek,et (potential, kinetic and total energy)
--- the tolerance and maxSteps can be the above defined simulation.minimisationTolerance and simulation.minimisationMaxSteps
+-- the tolerance and maxSteps can be the above defined minimisation.Tolerance and minimisation.MaxSteps
--
-- get_minimised_crdvels(tolerance,maxSteps) : this function returns a 2-tuple (crds,vels) containing
-- a copy of coordinates and velocities after minimisation.
@@ -196,6 +201,9 @@ simulation =
-- the tolerance and maxSteps can be the above defined simulation.minimisationTolerance and simulation.minimisationMaxSteps
-- NOTE lua indexing, starting from 1
--
+-- get_minimised_energy_crdvels(tolerance,maxSteps) : this returns a 5-tuple (ep,ek,et,crds,vels)
+-- This does the same as the two previous functions but with only one call
+--
-- hr_timer() : returns a variable representing a c++ high precision timer : can be used for measuring execution time.
-- do not try to modify it or even read it, it should only be used as argument for the following hr_timediff_* functions.
--
diff --git a/mol/ala2vac_transient/input_parRep.lua b/mol/ala2vac_transient/input_parRep.lua
index 73c2228..5f27bcb 100755
--- a/mol/ala2vac_transient/input_parRep.lua
+++ b/mol/ala2vac_transient/input_parRep.lua
@@ -64,15 +64,15 @@ OMMplatform_properties = { Threads = "1"}
-- load the integrator parameters from a serialised OpenMM XML file ;
-- no default, error if undefined
-- adapt the path to the file in the following line
-integrator = { xml = "./mol/ala2vac/Integrator.xml" }
+integrator = { xml = "./mol/ala2vac_transient/Integrator.xml" }
-- load the OpenMM System from a serialised XML file
-- no default, error if undefined
-system = { xml = "./mol/ala2vac/System.xml" }
+system = { xml = "./mol/ala2vac_transient/System.xml" }
-- load the OpenMM State from a serialised XML file
-- no default, error if undefined
-state = { xml = "./mol/ala2vac/State.xml" }
+state = { xml = "./mol/ala2vac_transient/State.xml" }
-- parameters for energy minimisation : tolerance and maximum number of steps (if 0 : no limit, continue until tolerance satsfied)
-- defaults are : minimisation={Tolerance=1e-6,MaxSteps=0}
@@ -184,7 +184,7 @@ simulation =
-- get_minimised_energy(tolerance,maxSteps) : this function returns the minimised energy of the system, using the OpenMM L-BFGS minimiser
-- note that coordinates are not affected, it just returns the minimum epot of the bassin in which dynamics currently evolves
-- it returns a 3-tuple ep,ek,et (potential, kinetic and total energy)
--- the tolerance and maxSteps can be the above defined simulation.minimisationTolerance and simulation.minimisationMaxSteps
+-- the tolerance and maxSteps can be the above defined minimisation.Tolerance and minimisation.MaxSteps
--
-- get_minimised_crdvels(tolerance,maxSteps) : this function returns a 2-tuple (crds,vels) containing
-- a copy of coordinates and velocities after minimisation.
@@ -193,6 +193,9 @@ simulation =
-- the tolerance and maxSteps can be the above defined simulation.minimisationTolerance and simulation.minimisationMaxSteps
-- NOTE lua indexing, starting from 1
--
+-- get_minimised_energy_crdvels(tolerance,maxSteps) : this returns a 5-tuple (ep,ek,et,crds,vels)
+-- This does the same as the two previous functions but with only one call
+--
-- hr_timer() : returns a variable representing a c++ high precision timer : can be used for measuring execution time.
-- do not try to modify it or even read it, it should only be used as argument for the following hr_timediff_* functions.
--
diff --git a/parRep.doxy b/parRep.doxy
index 3dd01ff..c127685 100644
--- a/parRep.doxy
+++ b/parRep.doxy
@@ -38,7 +38,7 @@ PROJECT_NAME = "gen.parRep"
# could be handy for archiving the generated documentation or if some version
# control system is used.
-PROJECT_NUMBER = v1.1.0
+PROJECT_NUMBER = v1.2.0
# Using the PROJECT_BRIEF tag one can provide an optional one line description
# for a project that appears at the top of each page and should give viewer a
diff --git a/run/run.bash.sh b/run/run.bash.sh
index 8755464..a4f2786 100755
--- a/run/run.bash.sh
+++ b/run/run.bash.sh
@@ -1,25 +1,38 @@
#!/bin/bash
+# --------------------------------------------------------------
+# create a local tmp dir for the run
+
ORIG_DIR=$PWD
TMP_DIR=$(mktemp --tmpdir=$ORIG_DIR -u)
mkdir -p $TMP_DIR && echo "Running in tmp dir $TMP_DIR"
cd $TMP_DIR
-ln -s $ORIG_DIR/mol
+ln -s $ORIG_DIR/../mol
+
+# --------------------------------------------------------------
+
+# use one of the following input files
-#ln -s $ORIG_DIR/mol/ala2vac/input_parRep.lua input.lua
-#ln -s $ORIG_DIR/mol/ala2vac/input_generalized_parRep.lua input.lua
+#ln -s mol/ala2vac/input_parRep.lua input.lua
+#ln -s mol/ala2vac/input_generalized_parRep.lua input.lua
-#ln -s $ORIG_DIR/mol/ala2vac_transient/input_parRep.lua input.lua
-#ln -s $ORIG_DIR/mol/ala2vac_transient/input_generalized_parRep.lua input.lua
+#ln -s mol/ala2vac_transient/input_parRep.lua input.lua
+ln -s mol/ala2vac_transient/input_generalized_parRep.lua input.lua
-ln -s $ORIG_DIR/mol/1d7h_i/input_generalized_parRep.lua input.lua
+#ln -s mol/1d7h_i/input_generalized_parRep.lua input.lua
-ln -s $ORIG_DIR/parRep
+ln -s $ORIG_DIR/../build/parRep
+# When using the OMM CPU platform this decides how many CPU threads per replica to use
export OPENMM_CPU_THREADS=1
+# if the following is not present only the slow Reference OMM platform will be available
export OPENMM_PLUGIN_DIR=$ORIG_DIR/../build/openmm-7.3.0/lib/plugins
-mpirun -np 1 ./parRep -i input.lua -log dbg : \
- -np 7 ./parRep -i input.lua -log dbg -o out.txt -e err.txt
+# to save random seeds to a file for later use
+mpirun -np 1 ./parRep -i input.lua -log dbg --out-seeds seeds.bin : \
+ -np 7 ./parRep -i input.lua -log dbg -o out.txt -e err.txt --out-seeds seeds.bin
+# to run using random seeds from a previous execution for reproducibility (not guaranteed to work on all platforms or if the software environnment is not exactly the same)
+# mpirun -np 1 ./parRep -i input.lua -log dbg --inp-seeds seeds.bin : \
+# -np 7 ./parRep -i input.lua -log dbg -o out.txt -e err.txt --inp-seeds seeds.bin
diff --git a/run/run.slurm.sh b/run/run.slurm.sh
index a284f04..f8a3678 100755
--- a/run/run.slurm.sh
+++ b/run/run.slurm.sh
@@ -12,12 +12,12 @@ TMP_DIR=$(mktemp --tmpdir=$ORIG_DIR -u)
mkdir -p $TMP_DIR && echo "Running in tmp dir $TMP_DIR"
cd $TMP_DIR
-ln -s $ORIG_DIR/mol
-ln -s $ORIG_DIR/mol/ala2vac_transient/input_generalized_parRep.lua input.lua
-ln -s $ORIG_DIR/parRep
+ln -s $ORIG_DIR/../mol
+ln -s mol/ala2vac_transient/input_generalized_parRep.lua input.lua
+ln -s $ORIG_DIR/../build/parRep
export OPENMM_CPU_THREADS=1
-export OPENMM_PLUGIN_DIR=/home/hedin/prog/parRep/build/openmm-7.3.0/lib/plugins
+export OPENMM_PLUGIN_DIR=$ORIG_DIR/../build/openmm-7.3.0/lib/plugins
srun --mpi=pmi2 ./parRep -i input.lua -log warn -o out.txt -e err.txt
diff --git a/src/dyna/GelmanRubin.cpp b/src/dyna/GelmanRubin.cpp
index b91aacb..26f5cf9 100644
--- a/src/dyna/GelmanRubin.cpp
+++ b/src/dyna/GelmanRubin.cpp
@@ -4,7 +4,7 @@
* \author Florent Hédin
* \author Tony Lelièvre
* \author École des Ponts - ParisTech
- * \date 2016-2018
+ * \date 2016-2019
*/
#include <cstdio>
@@ -19,10 +19,10 @@
using namespace std;
-GelmanRubinAnalysis::GelmanRubinAnalysis(const uint32_t _num_chains, const uint32_t _discard)
+GelmanRubinAnalysis::GelmanRubinAnalysis(const uint32_t _num_chains, const uint32_t _no_convergence_if_less_than)
{
num_chains = _num_chains;
- discard_first = _discard;
+ min_num_observations_before_check = _no_convergence_if_less_than;
// allocate one map<string,Observable> per chain
observablesList = unique_ptr<std::map<std::string,Observable>[]>(new map<string,Observable>[num_chains]);
@@ -40,6 +40,7 @@ void GelmanRubinAnalysis::registerObservable(const Observable& obs, const double
}
tolerance[obs.get_name()] = tol;
+
obsTypes.push_back(obs.get_name());
}
@@ -116,7 +117,8 @@ void GelmanRubinAnalysis::addNObservations(const std::string& obName,
const uint32_t chain_id)
{
map<string,Observable>& chainObsList = observablesList[chain_id];
- chainObsList[obName].add(from,to);
+ const vector<double> tmp(from,to);
+ chainObsList[obName].add(tmp);
}
void GelmanRubinAnalysis::updateStatistics()
@@ -180,7 +182,7 @@ void GelmanRubinAnalysis::describeChain(const uint32_t chain_id)
void GelmanRubinAnalysis::describe()
{
- LOG_PRINT(LOG_DEBUG,"GelmanRubinAnalysis desciption :\n");
+ LOG_PRINT(LOG_INFO,"GelmanRubinAnalysis desciption :\n");
for(string& name : obsTypes)
{
@@ -193,18 +195,23 @@ void GelmanRubinAnalysis::describe()
double GelmanRubinAnalysis::get_ratio(string& obName)
{
- double lratio = ratio[obName].back();
+ const double lratio = ratio[obName].back();
return lratio;
}
bool GelmanRubinAnalysis::check_convergence(const string& obName)
{
- vector<double>& rvec = ratio[obName];
+ const vector<double>& rvec = ratio[obName];
+ const Observable& obs = observablesList[0][obName];
+ const uint32_t num_obs_accumulated = (uint32_t) obs.get_numObs();
bool check;
- if(rvec.size() < discard_first)
+ if(num_obs_accumulated < min_num_observations_before_check)
+ {
+ LOG_PRINT(LOG_DEBUG,"Only %u observations available for observable %s, but the user required to only calculate the ratio when at least %u are available, therefore convergence is forced to be false.",num_obs_accumulated,obName.c_str(),min_num_observations_before_check);
check = false;
+ }
else
check = (fabs(1.0-rvec.back()) < tolerance[obName])? true : false ;
@@ -266,4 +273,7 @@ void GelmanRubinAnalysis::reset_all_chains()
ratio[st].resize(0);
ratio[st].shrink_to_fit();
}
+
+
+
}
diff --git a/src/dyna/observable.cpp b/src/dyna/observable.cpp
index 6e272a4..f061c97 100644
--- a/src/dyna/observable.cpp
+++ b/src/dyna/observable.cpp
@@ -4,7 +4,7 @@
* \author Florent Hédin
* \author Tony Lelièvre
* \author École des Ponts - ParisTech
- * \date 2016-2018
+ * \date 2016-2019
*/
#include <numeric>
diff --git a/src/dyna/parRep.cpp b/src/dyna/parRep.cpp
index 83582e6..24ed0f9 100644
--- a/src/dyna/parRep.cpp
+++ b/src/dyna/parRep.cpp
@@ -4,7 +4,7 @@
* \author Florent Hédin
* \author Tony Lelièvre
* \author École des Ponts - ParisTech
- * \date 2016-2018
+ * \date 2016-2019
*/
#include <cstdint>
diff --git a/src/dyna/parRep_FV.cpp b/src/dyna/parRep_FV.cpp
index 4a81d47..d9e6a89 100644
--- a/src/dyna/parRep_FV.cpp
+++ b/src/dyna/parRep_FV.cpp
@@ -4,7 +4,7 @@
* \author Florent Hédin
* \author Tony Lelièvre
* \author École des Ponts - ParisTech
- * \date 2016-2018
+ * \date 2016-2019
*/
#include <utility>
@@ -23,7 +23,6 @@
using namespace std;
-// NOTE to be used by a derived class only
ParRepFV::ParRepFV(DATA& _dat,
unique_ptr<ATOM[]>& _at,
unique_ptr<MD_interface>& _md,
diff --git a/src/dyna/parRep_FV_multiExits.cpp b/src/dyna/parRep_FV_multiExits.cpp
index 806aaa9..ae6c021 100644
--- a/src/dyna/parRep_FV_multiExits.cpp
+++ b/src/dyna/parRep_FV_multiExits.cpp
@@ -4,7 +4,7 @@
* \author Florent Hédin
* \author Tony Lelièvre
* \author École des Ponts - ParisTech
- * \date 2016-2018
+ * \date 2016-2019
*/
#include <cstdlib>
diff --git a/src/dyna/runSim.cpp b/src/dyna/runSim.cpp
index f1c21a5..d0896c3 100644
--- a/src/dyna/runSim.cpp
+++ b/src/dyna/runSim.cpp
@@ -32,7 +32,7 @@
using namespace std;
-void run_simulation(const string& inpf)
+void run_simulation(const string& inpf, const std::string& seeds_inp_file, const std::string& seeds_out_file)
{
/* find out MY process ID, and how many processes were started. */
int32_t my_id = -1;
@@ -43,27 +43,46 @@ void run_simulation(const string& inpf)
/*
* Each MPI rank will initialise its own random numbers generator ; see rand.hpp
*/
- init_rand();
- /*
- * This is how to initialize the random generators providing seeds manually (at least 5 seeds of type uint64_t encodes as a comma separated string)
- *
- * ATTENTION Using the directly the following line would be wrong because all the random generators on all
- * MPI ranks would be initialised with the same seeds !!! Each rank should use a modified vrsion of those seeds.
- *
- * const string ten_seeds = "12813764096581876017,1660537369292506360,10652021971297968194,"
- * "1632182917282226681,2681620832946811191,11894265520984226376,"
- * "8987030001799979889,7928564290642933391,6948006017539476486,1814061233159141676";
- *
- * init_rand(ten_seeds);
- *
- */
+ const bool load_seeds = (seeds_inp_file != NULLFILE);
+ const bool save_seeds = (seeds_out_file != NULLFILE);
+ try
+ {
+ if(load_seeds)
+ {
+ // seeds for all the ranks will be read from the file 'seeds_inp_file'
+ init_rand(SEEDS_IO::LOAD_FROM_FILE,seeds_inp_file);
+ fprintf(stdout,"All seeds for all ranks have been read from file %s\n",seeds_inp_file.c_str());
+ }
+ else if(save_seeds)
+ {
+ //Seeds from all the ranks will be written to the file 'seeds_out_file' for re-use
+ init_rand(SEEDS_IO::SAVE_TO_FILE,seeds_out_file);
+ fprintf(stdout,"All seeds from all ranks have been saved to file %s\n",seeds_out_file.c_str());
+ }
+ else
+ {
+ init_rand(SEEDS_IO::NONE);
+ fprintf(stdout, "Random Seeds internally generated but not saved to a file, no reproducibility will be possible for this run !!\n");
+ LOG_PRINT(LOG_WARNING,"Random Seeds internally generated but not saved to a file, no reproducibility will be possible for this run !!\n");
+ }
+ }
+ catch(exception& e)
+ {
+ fprintf(stderr,"std::exception captured on rank %d when initialising random numbers generator !\n",my_id);
+ fprintf(stderr,"Error message is : %s\n",e.what());
+ fprintf(stderr,"Now flushing files and terminating all other MPI ranks...\n");
+ LOG_FLUSH_ALL();
+ MPI_CUSTOM_ABORT_MACRO();
+ }
+
fprintf(stdout,"Rank %d properly initialised its mt19937 random numbers generator\n",my_id);
+
fprintf(stdout,"5 random uint32_t from rank %d : %u %u %u %u %u\n",my_id,get_uint32(),
get_uint32(),get_uint32(),get_uint32(),get_uint32());
- MPI_Barrier(MPI_COMM_WORLD);
+// MPI_Barrier(MPI_COMM_WORLD);
// stores simulation parameters
DATA dat;
@@ -219,6 +238,10 @@ void run_simulation(const string& inpf)
MPI_CUSTOM_ABORT_MACRO();
}
+ // ready to run sim, flush IO files before and let's start !
+ LOG_FLUSH_ALL();
+ MPI_Barrier(MPI_COMM_WORLD);
+
try
{
simulation->run();
diff --git a/src/logger.cpp b/src/logger.cpp
index 55900cf..33e5a67 100644
--- a/src/logger.cpp
+++ b/src/logger.cpp
@@ -4,7 +4,7 @@
* \author Florent Hédin
* \author Tony Lelièvre
* \author École des Ponts - ParisTech
- * \date 2016-2018
+ * \date 2016-2019
*/
#include <string>
diff --git a/src/main.cpp b/src/main.cpp
index eb9f00b..7261f9d 100644
--- a/src/main.cpp
+++ b/src/main.cpp
@@ -4,7 +4,7 @@
* \author Florent Hédin
* \author Tony Lelièvre
* \author École des Ponts - ParisTech
- * \date 2016-2018
+ * \date 2016-2019
*/
#include <string>
@@ -115,7 +115,9 @@ int main(int argc, char* argv[])
}
MPI_Barrier(MPI_COMM_WORLD);
- char inpf[FILENAME_MAX] = "";
+ string input_file = NULLFILE;
+ string seeds_inp_file = NULLFILE;
+ string seeds_out_file = NULLFILE;
// arguments parsing
for (uint32_t i=1; i<(uint32_t)argc; i++)
@@ -123,7 +125,7 @@ int main(int argc, char* argv[])
// get name of input file
if(!strcasecmp(argv[i],"-i"))
{
- sprintf(inpf,"%s",argv[++i]);
+ input_file = string(argv[++i]);
}
/*
* reopen stdout to unique file based on the user specified template name ;
@@ -167,6 +169,15 @@ int main(int argc, char* argv[])
else
fprintf(stdout,"[Warning] Unknown log level '%s' : default value used.\n\n",argv[i]);
}
+ // in or out files where random seeds are read/written
+ else if( !strcasecmp(argv[i],"--inp-seeds") )
+ {
+ seeds_inp_file = string(argv[++i]);
+ }
+ else if( !strcasecmp(argv[i],"--out-seeds") )
+ {
+ seeds_out_file = string(argv[++i]);
+ }
// print help and proper exit
else if( !strcasecmp(argv[i],"-h") || !strcasecmp(argv[i],"-help") || !strcasecmp(argv[i],"--help") )
{
@@ -208,8 +219,7 @@ int main(int argc, char* argv[])
}
// run the simulation
- const string input(inpf);
- run_simulation(input);
+ run_simulation(input_file,seeds_inp_file,seeds_out_file);
// closing log files is the last thing to do as errors may occur at the end
close_logfiles();
@@ -237,9 +247,11 @@ int main(int argc, char* argv[])
void help(char *argv[])
{
fprintf(stdout,"Need at least one argument : %s -i an_input_file\n",argv[0]);
- fprintf(stdout,"optional args : -o [stdout_to_file] -e [stderr_to_file] -log [logging level, one of { no | warn | info | dbg }] \n");
- fprintf(stdout,"Example : \n %s -i input_file -o out.txt -log info \n\n",argv[0]);
+ fprintf(stdout,"optional args : -o [stdout_to_file] -e [stderr_to_file] -log [logging level, one of { no | warn | info | dbg }] --inp-seeds [seed_file] --out-seeds [seed_file]\n");
+ fprintf(stdout,"Examples : \n\t%s -i input_file -o out.txt -log info --inp-seeds seeds.bin \n",argv[0]);
+ fprintf(stdout,"\t%s -i input_file -o out.txt -log info --out-seeds seeds.bin \n",argv[0]);
fprintf(stdout,"The default logging level is 'warn' \n");
+ fprintf(stdout,"See the documentation for more details. \n");
}
/*!
@@ -257,6 +269,9 @@ void help(char *argv[])
*
* This software is a new implementation of the Generalized parallel replica method, targeting
* frequently encountered metastable biochemical systems.
+ * It was used in version v1.0.0 for producing all the simulations published in :
+ *
+ * <a href="https://doi.org/10.1016/j.cpc.2019.01.005">https://doi.org/10.1016/j.cpc.2019.01.005</a>
*
* \image html main.png
* \image latex main.pdf
diff --git a/src/md/md_interface.cpp b/src/md/md_interface.cpp
index d8ca6c4..7a7fc14 100644
--- a/src/md/md_interface.cpp
+++ b/src/md/md_interface.cpp
@@ -4,7 +4,7 @@
* \author Florent Hédin
* \author Tony Lelièvre
* \author École des Ponts - ParisTech
- * \date 2016-2018
+ * \date 2016-2019
*/
#include "md_interface.hpp"
diff --git a/src/md/omm_interface.cpp b/src/md/omm_interface.cpp
index 3a8de17..e7398bd 100644
--- a/src/md/omm_interface.cpp
+++ b/src/md/omm_interface.cpp
@@ -4,7 +4,7 @@
* \author Florent Hédin
* \author Tony Lelièvre
* \author École des Ponts - ParisTech
- * \date 2016-2018
+ * \date 2016-2019
*/
#include <cstdio>
@@ -109,13 +109,6 @@ OMM_interface::OMM_interface(DATA& _dat,
// deserialize the integrator
integrator = unique_ptr<OpenMM::Integrator>(&integr);
- // deserialize an OpenMM state : coordinates velocities forces energy
- addPlatform(platformDesired,platformProperties,state);
-
- state.getPeriodicBoxVectors(pbc[0],pbc[1],pbc[2]);
-
- memcpy(dat.pbc.data(),pbc.data(),sizeof(pbc));
-
// re-seed the integrator
if(dynamic_cast<OpenMM::LangevinIntegrator*>(integrator.get()))
{
@@ -138,14 +131,22 @@ OMM_interface::OMM_interface(DATA& _dat,
throw runtime_error("Error when re-seeding the OMM integrator : integrator type is neither OpenMM::LangevinIntegrator or OpenMM::BrownianIntegrator !\n");
}
+ // from system and integrator, create a platform and a context using the state
+ addPlatform(platformDesired,platformProperties,state);
+
+ // retireve pbc data and store it in dat
+ state.getPeriodicBoxVectors(pbc[0],pbc[1],pbc[2]);
+ memcpy(dat.pbc.data(),pbc.data(),sizeof(pbc));
+
// fill dat structure with unserialised data
dat.natom = system->getNumParticles();
dat.timestep = integrator->getStepSize();
+ // retrieve initial positions and velocities
pos = vector<Vec3>(dat.natom);
vel = vector<Vec3>(dat.natom);
- /**
+ /*
* calculate the number of degrees of freedom of the system
* it is 3 times the number of particles (not virtual ones for which mass=0)
* minus the number of constraints
@@ -232,7 +233,7 @@ void OMM_interface::doNsteps(uint32_t numSteps)
void OMM_interface::minimise(double tol, uint32_t maxsteps, ENERGIES& ener)
{
// minimisation with L-BGFS
- OpenMM::LocalEnergyMinimizer::minimize(*context, (int32_t)tol, maxsteps);
+ OpenMM::LocalEnergyMinimizer::minimize(*context,tol,(int32_t)maxsteps);
// get back data from context and extract energies
const int infoMask = OpenMM::State::Energy;
@@ -295,6 +296,8 @@ void OMM_interface::minimiseWithCopy(double tol, uint32_t maxsteps, ENERGIES& en
context->setPositions(lpos);
context->setVelocities(lvel);
+ context->applyConstraints(integrator->getConstraintTolerance());
+ context->applyVelocityConstraints(integrator->getConstraintTolerance());
}
void OMM_interface::minimiseWithCopy(double tol, uint32_t maxsteps, ENERGIES& ener,
@@ -327,6 +330,9 @@ void OMM_interface::minimiseWithCopy(double tol, uint32_t maxsteps, ENERGIES& en
context->setPositions(lpos);
context->setVelocities(lvel);
+ context->applyConstraints(integrator->getConstraintTolerance());
+ context->applyVelocityConstraints(integrator->getConstraintTolerance());
+
}
void OMM_interface::setSimClockTime(double t)
@@ -346,9 +352,13 @@ void OMM_interface::setCrdsVels(ATOM at[])
memcpy(pbc.data(),dat.pbc.data(),sizeof(pbc));
context->setPeriodicBoxVectors(pbc[0],pbc[1],pbc[2]);
+
context->setPositions(pos);
context->setVelocities(vel);
+ context->applyConstraints(integrator->getConstraintTolerance());
+ context->applyVelocityConstraints(integrator->getConstraintTolerance());
+
}
void OMM_interface::randomiseVelocities()
diff --git a/src/rand.cpp b/src/rand.cpp
index 6b34824..c7bd62d 100644
--- a/src/rand.cpp
+++ b/src/rand.cpp
@@ -4,7 +4,7 @@
* \author Florent Hédin
* \author Tony Lelièvre
* \author École des Ponts - ParisTech
- * \date 2016-2018
+ * \date 2016-2019
*/
#include <random>
@@ -12,12 +12,14 @@
#include <limits>
#include <iostream>
#include <sstream>
+#include <stdexcept>
#include <cstdio>
+#include "global.hpp"
#include "logger.hpp"
+#include "mpi_utils.hpp"
#include "rand.hpp"
-#include "global.hpp"
using namespace std;
@@ -40,54 +42,141 @@ static uniform_int_distribution<int32_t> int32_t_dist(numeric_limits<int32_t>::
static uniform_int_distribution<uint32_t> uint32_t_dist(numeric_limits<uint32_t>::min(),
numeric_limits<uint32_t>::max());
-void init_rand()
+void init_rand(SEEDS_IO io_type, const string& seeds_file_name)
{
- random_device rd;
- LOG_PRINT(LOG_INFO,"Entropy of 'std::random_device' on this machine is : %lf\n",rd.entropy());
- /*
- * Each MPI ranks reads 512 uint64_t 'seeds' from random_device
- */
- array<uint64_t,512> seeds;
- for(size_t n=0; n<512; n++)
- seeds[n] = (uint64_t) rd();
-
- // then each ranks initialises its own mt19937 with the seed sequence
- seed_seq seq(seeds.begin(),seeds.end());
- generator.seed(seq);
-
- // finally for each generator we discard some random numbers to be sure to avoid any correlation
- generator.discard(4*generator.state_size);
+ bool save_seeds = false;
+ bool load_seeds = false;
-}
+ switch(io_type)
+ {
+ case SEEDS_IO::NONE :
+ break;
+
+ case SEEDS_IO::SAVE_TO_FILE :
+ save_seeds = true;
+ break;
+
+ case SEEDS_IO::LOAD_FROM_FILE :
+ load_seeds = true;
+ break;
+
+ default:
+ throw runtime_error("Error : unsupported SEEDS_IO mode in function '" + string(__func__) + "' !!!");
+ break;
+ }
-void init_rand(const string& seeds_str)
-{
- istringstream split_stream(seeds_str);
- string token;
-
- vector<uint64_t> seeds;
-
/*
- * Exctract the seeds from the comma separated string ;
- * NOTE It is crucial to use a different string on each MPI rank otherwise the generators will provide the same numbers on each rank !!
+ * Each MPI ranks requires 512 uint32_t 'seeds' for initialising its mt19937_64 generator
*/
- while(getline(split_stream, token, ','))
+ array<uint32_t,512> seeds;
+
+ if(!load_seeds)
{
- seeds.push_back((uint64_t)stoul(token));
+ random_device rd;
+ LOG_PRINT(LOG_INFO,"Entropy of 'std::random_device' on this machine is : %lf\n",rd.entropy());
+
+ for(size_t n=0; n<512; n++)
+ seeds[n] = (uint32_t) rd();
}
-
- if(seeds.size()<5)
+ else
{
- throw runtime_error("Error : when providing a string encoded list of seeds for the random generator, you need to provide at least 5 seeds, but it appears that you have provided only " + to_string(seeds.size()) + " seeds ! \n");
+ // First enforce sync with a barrier
+ MPI_Barrier(MPI_COMM_WORLD);
+
+ int32_t myrank = -1, numranks = -1;
+ MPI_Comm_rank(MPI_COMM_WORLD, &myrank);
+ MPI_Comm_size(MPI_COMM_WORLD, &numranks);
+
+ // see comments in the following block "if(save_seeds)..." below
+
+ MPI_File myfile;
+
+ MPI_File_open(MPI_COMM_WORLD, seeds_file_name.c_str(),
+ MPI_MODE_RDONLY, MPI_INFO_NULL,
+ &myfile);
+
+ MPI_Offset disp = myrank*512*sizeof(uint32_t);
+
+ MPI_File_set_view(myfile, disp,
+ MPI_UINT32_T, MPI_UINT32_T,
+ "native", MPI_INFO_NULL);
+
+ MPI_File_read(myfile, seeds.data(),
+ 512, MPI_UINT32_T,
+ MPI_STATUS_IGNORE);
+
+ MPI_File_close(&myfile);
+
+ MPI_Barrier(MPI_COMM_WORLD);
}
+ /*
+ * Save seeds to a unique binary file for reproducibility
+ */
+ if(save_seeds)
+ {
+ // First enforce sync with a barrier
+ MPI_Barrier(MPI_COMM_WORLD);
+
+ int32_t myrank = -1, numranks = -1;
+ MPI_Comm_rank(MPI_COMM_WORLD, &myrank);
+ MPI_Comm_size(MPI_COMM_WORLD, &numranks);
+
+ /* Shared file */
+ MPI_File myfile;
+
+ /*
+ * Open the file : it will be a binary file open in write mode, with overwritting
+ *
+ * int MPI_File_open(MPI_Comm comm, const char *filename,
+ * int amode, MPI_Info info,
+ * MPI_File *fh)
+ */
+ MPI_File_open(MPI_COMM_WORLD, seeds_file_name.c_str(),
+ MPI_MODE_CREATE | MPI_MODE_WRONLY, MPI_INFO_NULL,
+ &myfile);
+
+ // displacment == absolute displacement in the file where the data will be written (from beginning of file)
+ MPI_Offset disp = myrank*512*sizeof(uint32_t);
+
+ /*
+ * Set the file view : tels at which value of disp data will be written
+ *
+ * MPI_File_set_view(MPI_File fh, MPI_Offset disp,
+ * MPI_Datatype etype, MPI_Datatype filetype,
+ * const char *datarep, MPI_Info info)
+ */
+ MPI_File_set_view(myfile, disp,
+ MPI_UINT32_T, MPI_UINT32_T,
+ "native", MPI_INFO_NULL);
+ /*
+ * Write buf to the file :
+ *
+ int MPI_File_write(MPI_File fh, const void *buf,
+ int count, MPI_Datatype datatype,
+ MPI_Status *status)
+ */
+ MPI_File_write(myfile, seeds.data(),
+ 512, MPI_UINT32_T,
+ MPI_STATUS_IGNORE);
+
+ /*
+ * Close the file
+ */
+ MPI_File_close(&myfile);
+
+ MPI_Barrier(MPI_COMM_WORLD);
+ }
+
+ // then each ranks initialises its own mt19937 with the seed sequence
seed_seq seq(seeds.begin(),seeds.end());
generator.seed(seq);
// finally for each generator we discard some random numbers to be sure to avoid any correlation
generator.discard(4*generator.state_size);
+ MPI_Barrier(MPI_COMM_WORLD);
}
double get_double_unif_0_1()
diff --git a/src/utils/lua_interface.cpp b/src/utils/lua_interface.cpp
index 0ea5570..b636816 100644
--- a/src/utils/lua_interface.cpp
+++ b/src/utils/lua_interface.cpp
@@ -4,7 +4,7 @@
* \author Florent Hédin
* \author Tony Lelièvre
* \author École des Ponts - ParisTech
- * \date 2016-2018
+ * \date 2016-2019
*/
#include <iostream>
@@ -13,7 +13,6 @@
#include <exception>
#include "sol.hpp"
-#include "rand.hpp"
#include "logger.hpp"
#include "lua_interface.hpp"
#include "omm_interface.hpp"
@@ -34,6 +33,7 @@ typedef tuple<double,double,double> VelsTuple;
typedef tuple<double,double,double> COMTuple;
typedef tuple<double,double,double> ENETuple;
typedef tuple<sol::table,sol::table> SolTabl2Tuple;
+typedef tuple<double,double,double,sol::table,sol::table> ENETupleSolTablTuple;
luaInterface::luaInterface(const string _inputFile,
DATA& _dat,
@@ -165,12 +165,12 @@ void luaInterface::parse_lua_file()
if(pvMap["algorithm"] == "PARREP")
{
- fprintf(stdout,"Parsing ParRep simulation parameters.");
+ fprintf(stdout,"Parsing ParRep simulation parameters.\n");
parse_parrep_params();
}
else if(pvMap["algorithm"] == "PARREP_FV")
{
- fprintf(stdout,"Parsing ParRep FV simulation parameters.");
+ fprintf(stdout,"Parsing ParRep FV simulation parameters.\n");
parse_parrepFV_params();
uint32_t num_ex_events = (uint32_t) stoul(pvMap.at("allowedExitEvents"));
@@ -473,20 +473,6 @@ void luaInterface::register_default_lua_functions()
return ENETuple(minE.epot(),minE.ekin(),minE.etot());
}
);
-
- // returns energy contribution
- lua.set_function("get_group_epot",[&](int32_t groups) -> double
- {
- if(md->engine_supports_groups_splitting()==false)
- {
- runtime_error ex("Lua script called get_group_energy(...) but it seems that the currently active MD engine does not support the group splitting feature ! Aborting program, please fix your input file accordingly...");
- throw ex;
- }
- ENERGIES e;
- md->getSimData(nullptr,&e,nullptr,nullptr,groups);
- return e.epot();
- }
- );
// perform energy minimisation, and return the coordinates and velocities
lua.set_function("get_minimised_crdvels",
@@ -521,13 +507,63 @@ void luaInterface::register_default_lua_functions()
return SolTabl2Tuple(crds,vels);
}
-);
+ );
+
+ lua.set_function("get_minimised_energy_crdvels",
+ [&](double tolerance, uint32_t maxSteps) -> ENETupleSolTablTuple
+ {
+ ENERGIES minE;
+ vector<XYZ> coordinates;
+ vector<XYZ> velocities;
+
+ md->minimiseWithCopy(tolerance,maxSteps,minE,coordinates,velocities);
+
+ sol::table crds = lua.create_table();
+ crds["x"] = lua.create_table();
+ crds["y"] = lua.create_table();
+ crds["z"] = lua.create_table();
+
+ sol::table vels = lua.create_table();
+ vels["x"] = lua.create_table();
+ vels["y"] = lua.create_table();
+ vels["z"] = lua.create_table();
+
+ for(uint32_t i=1; i<=dat.natom; i++)
+ {
+ crds["x"][i] = coordinates[i-1][0];
+ crds["y"][i] = coordinates[i-1][1];
+ crds["z"][i] = coordinates[i-1][2];
+
+ vels["x"][i] = velocities[i-1][0];
+ vels["y"][i] = velocities[i-1][1];
+ vels["z"][i] = velocities[i-1][2];
+ }
+
+ return ENETupleSolTablTuple(minE.epot(),minE.ekin(),minE.etot(),crds,vels);
+ }
+ );
+
// get the current temperature (from kinetic energy)
lua.set_function("get_temperature",
[&]()->double{ return md->getTemperature();}
);
+ // returns energy contribution
+ lua.set_function("get_group_epot",
+ [&](int32_t groups) -> double
+ {
+ if(md->engine_supports_groups_splitting()==false)
+ {
+ runtime_error ex("Lua script called get_group_energy(...) but it seems that the currently active MD engine does not support the group splitting feature ! Aborting program, please fix your input file accordingly...");
+ throw ex;
+ }
+ ENERGIES e;
+ md->getSimData(nullptr,&e,nullptr,nullptr,groups);
+ return e.epot();
+ }
+ );
+
}
void luaInterface::register_default_lua_variables()
@@ -541,7 +577,7 @@ void luaInterface::register_default_lua_variables()
lua["temperature"] = 0.0;
lua.create_named_table("minimisation");
- lua["minimisation"]["Tolerance"] = 1e-6;
+ lua["minimisation"]["Tolerance"] = 10.0;
lua["minimisation"]["MaxSteps"] = 0;
lua["equilibrationSteps"] = 50e3;
@@ -559,13 +595,13 @@ void luaInterface::register_default_lua_variables()
// expose dummy lambda functions ; they do nothing
// and are the default in case the user does not override them in the lua script
- lua["SQLiteDB"]["open"].set_function([](){fprintf(stderr,"Dummy lambda db_open\n");});
+ lua["SQLiteDB"]["open"].set_function([](){fprintf(stdout,"Dummy lambda db_open\n");});
- lua["SQLiteDB"]["close"].set_function([](){fprintf(stderr,"Dummy lambda db_close\n");});
+ lua["SQLiteDB"]["close"].set_function([](){fprintf(stdout,"Dummy lambda db_close\n");});
- lua["SQLiteDB"]["insert_state"].set_function([](){fprintf(stderr,"Dummy lambda db_insert\n");});
+ lua["SQLiteDB"]["insert_state"].set_function([](){fprintf(stdout,"Dummy lambda db_insert\n");});
- lua["SQLiteDB"]["backup_to_file"].set_function([](){fprintf(stderr,"Dummy lambda db_backup\n");});
+ lua["SQLiteDB"]["backup_to_file"].set_function([](){fprintf(stdout,"Dummy lambda db_backup\n");});
}
@@ -588,7 +624,7 @@ void luaInterface::parse_omm_params()
p = OMM_interface::omm_platforms_names.at(PLATFORMS::CUDA);
else
{
- throw runtime_error(string("Error when parsing OMM platform : '" + p + "' is not a valid platform type !"));
+ throw runtime_error(string("Error when parsing OMM platform : '" + p + "' is not a valid platform type !\n"));
}
fprintf(stdout,"Parsed OMM platform (or default if missing) : %s\n",pvMap["platform"].c_str());
@@ -650,7 +686,7 @@ void luaInterface::parse_omm_params()
}
else
{
- fprintf(stdout,"No specific OpenMM platform properties defined by the user...");
+ fprintf(stdout,"No specific OpenMM platform properties defined by the user...\n");
}
}
@@ -699,7 +735,7 @@ void luaInterface::parse_parrepFV_params()
}
else
{
- throw runtime_error( "Error when trying to read simulation.GRobservables which is required for a 'parrepFV' simulation !\n");
+ throw runtime_error( "Error when trying to read 'simulation.GRobservables' which is required for a 'parrepFV' simulation !\n");
}
pvMap["GRtol"] = to_string(simulation_def["GRtol"].get_or(0.01));
diff --git a/src/utils/mpi_utils.cpp b/src/utils/mpi_utils.cpp
index c6283e0..9cf093f 100644
--- a/src/utils/mpi_utils.cpp
+++ b/src/utils/mpi_utils.cpp
@@ -4,7 +4,7 @@
* \author Florent Hédin
* \author Tony Lelièvre
* \author École des Ponts - ParisTech
- * \date 2016-2018
+ * \date 2016-2019
*/
#include <cstdint>
--
GitLab