diff --git a/.gitignore b/.gitignore
index 2804c78bfac80df5614800927c546cd9838ae1f3..b51356f18624f94dcccef451621b86c5c46a4187 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 0000000000000000000000000000000000000000..8459704768eef654dbb5cfb98439d353fd158410
--- /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 dfb604cadd3821852823cbf678b312eafed9707d..0a7d665ac491edbb2b5bd0a1b3cac80041dc9f67 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 b27dc9a02bb34a89ef223f6fdd5e90c504e48f25..91af3b1f2cadcea17e60d74cef0cad42bec731d7 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 9e1f08e5abca85518c0768e8ce08d2e23ddfabcc..3a0db8c8a6fffcecc226a75e51b46e2984e11cfe 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 9e0f2804a10a4ae9a1d5c1f42ae3f2662dd53e36..3e5253b2af422635492b88d9051ae8ee9ae5e6fe 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.
+
+[![pipeline status](https://gitlab.inria.fr/parallel-replica/gen.parRep/badges/test_ci/pipeline.svg)](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 976dad2c168b9b7483233a2d60f57282071fd09a..f27d3c5b5b9e8b409db47ee4267dd4adf83d43b7 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 0000000000000000000000000000000000000000..40eb0f24b04ede50b1c7611ee4284c1a2d99eb7e
--- /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 0000000000000000000000000000000000000000..fbf5afdb0ae785e2bacf8edd520c4dd2130baede
--- /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 0000000000000000000000000000000000000000..fec690df653ce9e7ffbeaae09dd2882376759343
--- /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 0000000000000000000000000000000000000000..b1cb0e6c426218daa3a6dc82cb94daab142de776
--- /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 0000000000000000000000000000000000000000..03ed2d3414fce7e6293caeb29e7365d310f15064
--- /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 0000000000000000000000000000000000000000..ecb721352b1dff9a3eea5391ebd6b7e757d154a7
--- /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 0000000000000000000000000000000000000000..70c82e0be91c7d8db5e3d83829f508189c391e03
--- /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 0000000000000000000000000000000000000000..9da4e99cbb4c645c897a336293861a01836e12c8
--- /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
Binary files /dev/null and b/ci/ref.sim/lj7.db differ
diff --git a/ci/ref.sim/seeds.bin b/ci/ref.sim/seeds.bin
new file mode 100644
index 0000000000000000000000000000000000000000..fc3504f94422bcdb14a040dda024abb0a387d279
Binary files /dev/null and b/ci/ref.sim/seeds.bin differ
diff --git a/ci/run.ci.sh b/ci/run.ci.sh
new file mode 100755
index 0000000000000000000000000000000000000000..8c2d99bc45e3e980efae94e8b4f7f1170e5a9e2c
--- /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 fefc5cf652248eeea0d80d60ab3797c1457a66ba..b33c672e60cdee3dbe5f6694847efcc0e3a9bb54 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 113b3595793f9f82fae1c093fee9b2128ad9ec68..e01b70dd4f4b7a9a0ad540a03b3bbb9a84610ba1 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 9e6efddff686b6828a943028460d9af7a8256b3b..56ca276e78e17dc46a900edef0c750b42d431c85 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 0a40dc9cf51baadcecb1923581eae3c7b59cd3e9..3d610ec59aa4243519feb3883e1f92297af757ad 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 6006d443f7a5879a045e94c83ac3953eb9109ef1..05a6d0858a264658aae536bd0d3a69b297931274 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 7d698b5fcf5ab16aa02ea392bc1d68e26cefdd36..85a1aadd78e1284a1ef9f374f5c1c701ab29e5d8 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 c0b144e3729d71bb052543d2ba13a623054281c2..bfee42653079f5c0e18bd1d5ed7595acc5ca37e5 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 78ece31e9360ae2fbbd32f9ba7ec9bad005150e5..279d41c5c7b26a9d257b2595ecf37a2fb065b275 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 0af6c4a22219eb0f506d56436c92a8eb9741e838..7b2d6ea1f5679502ef3a99d92b973d0b31cdea53 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 7da5996ea95255eefa0bdf8d3744dd2003bb7dcc..d8a0f9339a97dd0c1379338eaed2aa39a7cded47 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 c3acb1db5e800a136352e2966cca105b63619f23..29ec520762dc055ce91c91fe669e9170527c8e72 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 ad8b181f1afe785e3847ea86add9706a47d4779e..9f9ebac4f93613bd210b733294c66658bf66fb50 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 d35feeddc2f2ec8ea0cbcb58eba14ac848eb99ea..0616bb2cd8ae8edaee748d017d2bae6aa99ababa 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 3f22630c297808e91994df42499cafef964cc080..c3f45e53c2c2317edb9a978d7b39c94a8f7272fb 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 1dcbeda9d7238dc05d2f4c164310d06f5fd5f5d4..25b1e50bee21c81e2a400740188beeb833c685a0 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 fab791fd72d41f8b76933b05a449d2cf1028a7b0..6e181b13fcfdba5462ec283354d098c4dcb801fb 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 b3ead1b97e3e39a3e596355e21435348ed4bd6d4..e2ef33ce610150de374a621945237efbed799ae5 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 980884857dbabe6cfbd5025e78b689a0b425774e..acb560ec4b4d93d8124bf023a48dd8c5ea0feb70 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 88d5d579330e5fb34ab50784deef42cf5618798d..747581271f6ba47665faba24c28e04abde469f68 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 73c22283b10402089ef6afae70d086078b337ca7..5f27bcb8a74ce7d4e9f1d6ac145a16207e38dc1f 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 3dd01ff78785f33c0ce35394e45d8d452dbfc8fb..c1276859542c4a6c9656c0a316ed40d548fdb8a9 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 875546455f95dd7ec865c1601c7d6251600c7f5e..a4f2786db482c37297d5494bd864d06a366e7629 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 a284f044d36503eaec5ad67b4b50efd7a68ce25b..f8a3678a6746c44bfd2b70d65d4fcecec07082cf 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 b91aacbb1a563d6ca41632fe09dfa198a6285ad6..26f5cf9e17bd39507a323d0ccbdb6231addc931f 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 6e272a44c9e1ccfd10f66f098c4363299897c909..f061c97234786b56cbbc1f5aed771d44b5018b44 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 83582e667d9eb1f23ff31be229d7a420bd368bbb..24ed0f98b034385c11fbf85caf74b5e6c2ce7a2e 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 4a81d47af0b369073490d079522810f6c50d74f5..d9e6a8988eb80feb045694eeaf2cde9c7a8007fe 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 806aaa95e7862a637845dcc1803ee24c9c5f5635..ae6c0218a3f20e0cf6e3d94a109c059def2694df 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 f1c21a5bb52de87afa3630ce41d33dce403a9eb5..d0896c32214c1c18175fb7f63c74399ba1cc259f 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 55900cf8845943c8ca37f599faba5bdfb72ce835..33e5a6793e5799e202bbb4bda5ad5d062e20a426 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 eb9f00b7533b94b07855f6caf8b630d04579e14f..7261f9da9d1c9dcd243e64e710341d1ddd3a136b 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 d8ca6c4fb4662f48ab0aec4035a7a35b2fa53791..7a7fc1454220bd4a0ec687214daf889ed2c36135 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 3a8de17187c7947a886db1435ed3268b0b280649..e7398bd28060c59aec9462dd3538e5e6de3ddd8b 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 6b348245d83044b526db18be8ad8f94948143bdb..c7bd62d850f1f8ca84a445215f6efaeba6ccec7b 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 0ea55703a72266600dbfed7655f06268a7f02de5..b636816ba95c52e36361a6ecb9018377bce5a8d2 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 c6283e09281da5a9e52de26849f80b29646d6973..9cf093f375cf4406b43dbc7d5a219b163a84c89c 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>