rinikerlab/ensembler

Ensembler is a python package that provides fast and easy access to 1D and 2D model system simulations. It can be used for method development or to deepen understanding of a broad spectrum of modeling methods, from basic sampling techniques to enhanced sampling and free energy calculations. It is easy to install, fast, increases shareability, comparability, and reproducibility of scientific code developments.

Welcome to Ensembler

Description

Ensembler is a python package that allows fast and easy access to one and two-dimensional model systems simulations. It enables method development using small test systems and helps deepen understanding of a broad spectrum of molecular dynamics (MD) methods, starting from basic techniques to enhanced sampling and free energy calculations. It is easy to install, fast, increases shareability, comparability, and reproducibility of scientific code developments. Here, we provide insights into the package's implementation, usage, and an application example for free energy calculation. There is an article about Ensembler in JCIM, check it out: https://doi.org/10.1021/acs.jcim.0c01283

Quick Start to Simulations:

Contents

The full documentation can be found here: https://rinikerlab.github.io/Ensembler/

potentials - Potential Energy Functions

Implement mathematical functions of interest for modeling purposes, for example, in chemistry. Implementation of new potential energy functions is straightforward, as there are only a few functions that need to be overwritten. Implemented Potentials: Harmonic Oscillator, Wave functions, etc. Also, different dimensionalities can be used, like 1D, 2D, and ND.

samplers - Sampling Methods

This module provides integrators for integrating potential functions. E.g., Monte Carlo, Velocity Verlet.

systems - Used for Simulations

This module is used to set up a simulation with a system class that contains a potential energy function and a sampling method, plus optional other parameters.

ensemble - Multi Replica Approaches

This module contains the implementation of the ConveyorBelt and will contain in future additional Replica Exchange approaches.

visualization

This module contains predefined visualization, animation and widget functions.

analysis

This module contains at the moment Free Energy Calculations.

Tutorials and Examples:

You can try Ensembler online via the binder or Google colab links below or on your machine with our provided jupyter notebooks showing use cases for Ensembler with visualizations. The Binder links might take some time to build the repo; depending on the traffic, it can take up to 10 minutes to be built.

Tutorials

We provide short introductions into how potential energy functions can be used and sampled in simulations in Ensembler with the jupyter notebooks in the example folder that can be also accessed by the binder links.

Here is an small example of a simple simulation:

CODE:

##Imports
from ensembler.potentials.OneD import fourWellPotential
from ensembler.samplers.stochastic import langevinIntegrator
from ensembler.system import system
from ensembler.visualisation.plotSimulations import oneD_simulation_analysis_plot

##Simulation Setup
V = fourWellPotential(Vmax=4, a=1.5, b=4.0, c=7.0, d=9.0,  ah=2., bh=0., ch=0.5, dh=1.)
sampler = langevinIntegrator(dt=0.1, gamma=10)
sys = system(potential=V, sampler=sampler, start_position=4, temperature=1)

##Simulate
sys.simulate(steps=1000)

##Visualize
positions = np.linspace(start=0, stop=10, num=1000) #phase space to be visualized
oneD_simulation_analysis_plot(system=sys, title="Langevin Simulation", limits_coordinate_space=positions)

OUT:

In the following links you can find more features of Ensembler.

Potentials:

Simulations:

Examples

Examples are advanced jupyter notebooks, covering a specific topic, going deeper into the methodology and having nice visualizations.

Enhanced Sampling:

Free Energy Calculations:

Interactive ConveyorBelt:

EDS-Potentials:

How To Install

There are many ways to get the Ensembler package.

• If you want to use Ensembler (don't forget to download the examples):

  pip install ensembler-rinikerlab
  

• If you want to use Ensembler, see the example notebooks or develop Ensembler:

  git clone https://github.com/rinikerlab/Ensembler.git
  cd Ensembler
  pip install .
  

• If you want to use Ensembler, see the example notebooks or develop Ensembler - but not install the package: Add the path to the Ensembler repository; the requirements needed for the package can be used like in the following examples:

  • PIP:

      git clone https://github.com/rinikerlab/Ensembler.git
      cd Ensembler
      pip install .
    

  • Anaconda:

      git clone https://github.com/rinikerlab/Ensembler.git
      cd Ensembler
      conda create -n Ensembler --file devtools/conda-envs/environment_unix.yml
      conda activate Ensembler
      pip install .
    

For windows, we also provide the required files in the same folders:

• devtools/pip_requirements/requirements_windows.txt
• devtools/conda-envs/environment_windows.yml).

Contributing

If you would like to contribute to Ensembler, you are most welcome! Just raise an issue or write us a mail.

Authors

Benjamin Ries; Stephanie M. Linker; David F. Hahn

Copyright

Copyright (c) 2020, Benjamin Ries, Stephanie M. Linker, David F. Hahn

Acknowledgements

Project-based on the Computational Molecular Science Python Cookiecutter version 1.3.

ensembler's Issues

Metadynamic simulation does not start with MonteCarlo integrator

Metadynamics simulation does not start when integrator is MonteCarlo type. No error is raised but simulation is stuck in the first step.
Example:

import os, sys
my_path = os.getcwd()+"/.."
print(my_path)
sys.path.append(my_path)

import numpy as np
from matplotlib import pyplot as plt

from ensembler.potentials.TwoD import wavePotential

from ensembler.potentials.biasTwoD import addedPotentials2D, metadynamicsPotential2D
from ensembler.integrator.stochastic import monteCarloIntegrator, metropolisMonteCarloIntegrator
from ensembler.system import system

##Visualisation
from ensembler.visualisation.plotSimulations import static_sim_plots, static_sim_plots_bias

#Simulation Setup
sim_steps = 2000
origpot = wavePotential(amplitude=(10,10), multiplicity=[1/6.,1/6.], degree=False)

#Add the bias and the original system
totpot = metadynamicsPotential2D(origpot, amplitude=1, sigma=(5,5), n_trigger=10, bias_grid_min=(0,0), bias_grid_max=(100,100), numbins=(1000,1000))

integrator = metropolisMonteCarloIntegrator(randomnessIncreaseFactor=5)
sys=system(potential=totpot, integrator=integrator, position=np.array([60,60]), temperature=1)

#simulate
cur_state = sys.simulate(sim_steps, withdrawTraj=True, initSystem=False)

print("Trajectory length: ",len(sys.trajectory))
print()
print("last_state: ", cur_state)
print(len(sys.trajectory))
sys.trajectory.head()

2D - Wavepotentials don't set radians/degree decision in constructor!

Needed to fix this manually after construction.

ConveyorBelt: lambda of perturbed system get values assigned larger than 1.

Code snippet:

import ensembler.ensemble.replicas_dynamic_parameters as cvb
ensemble = cvb.conveyorBelt(0.0, 8)
ensemble

Langevin Integrator incompatible with RE-EDS Ensemble

Trying to simulate a replica-exchange enveloping-distribution sampling ensemble while using the LangevinIntegrator as a sampler causes various overflows after a few exchanges. The trajectories' positions become NaN after the overflow and the positions sampled before that equate to physical nonsense (see the following graph).

I believe this happens due to the BBK-style implementation using the positional difference between the current and the last sampled point. In the step following an exchange, this difference in position can get large and will shoot the sampler up the walls of the potential until after a few exchanges, the energy values simply get to large for calculation.

I don't know whether LangevinIntegrator and RE-EDS can be made to work together. If not, a warning would be appreciated for those who try and use that combination.

(I am aware now of the LangevinVelocityIntegrator, but that requires its separate issue)

Bring RE to Ensembler

• build RE-Graph
• simulated tempering
  • implement
  • test
• RE-EDS
  • implement
  • test

Refactoring

• manage typing via util/ensemblerTypes
• System: init_Velocity

Integrators

• Langevine
  • implement
  • test
    • 1D
    • 2D

Compatibility Issues with LangevinVelocityIntegrator

I have used the LangevinVelocityIntegrator on a one-dimensional potential, and it outputs its trajectory positions (and velocities, I think) as an array containing two values for each step. While I don't know the implementation reasons or which one of them is the actual position, I have found this leads to issues down the line:

• Systems with the LangevinVelocityIntegrator as a sampler fail to work with oneD_simulation_analysis_plot()
• The LangevinVelocityIntegrator does not work with replicaExchangeEnvelopingDistributionSampling at least during my testing due to the exchange criterion failing to calculate anything using the two-value array.
• This one is more of an incovenience: I have not found a way to conveniently slice one of those positional values out of the trajectory (pandas.DataFrame containing these arrays as objects) without looping (but I am no python expert). The LangevinIntegrator in contrast just contains the positions as floats.
• The documentation gives no hint as to what these pairs of values signify and I have not dug far enough into the code to understand, whether one of these values might be a half-step or just the position in the second dimension (which I should not have in my 1D case).

Problem with nDim with 2D potentials

Using a 2D potential currently creates in error in the ##set dim part of basic_systems. This is due to the sympy representation of potential.nDim. Fix when sympy variables or corresponding valuesare used.

DocStrings:

Write of review Doc-Strings in Numpy Style.

• Potentials

  • OneD @linkerst && @schroederb
  • TwoD @linkerst
    - [ ] sub task: Think about how to implement V_off in 2D potentials e.g. harmonic oscillator @schroederb
  • ND @schroederb
  • _basic @schroederb

• Integrators @linkerst

  • _basic @schroederb
  • Newtonian
  • stochastic
  • Optimization

• Systems @dfhahn

  • basic_system
  • perturbedSystem
  • eds_system @schroederb

• Ensembles @dfhahn

  • dynamicReplica
  • _replicaGraph
  • ReplicaEx
  • Conveyor Belt

• Integrators @schroederb

  • _basicConditions
  • boxCondition
  • thermostat
  • restraint

• Analysis @schroederb

• SETUP Documentation

Introduction

CI integration

• codecov: exclude visualization
• Sphinx: setup pipeline
• Jupyter: extra Workflow

Enhanced sampling

• Umbrella sampling

• Metadynamics

• Add custom potential

• aMD?

Conclusion

Units of Measurement and Temperature Issues

Units of measurement appear to be undefined in most parts of Ensembler, which may not be an issue for many applications.

However, I ran into inconveniences with unspecified and maybe inconsistent units when I was working on my reweighting project. Here, the temperature factors into some of the reweighting methods. This became apparent, when I used the LangevinIntegrator and afterwards a MetropolisMonteCarloIntegrator as samplers. While a temperature of 298 for the MetropolisMonteCarlo sampler produced a sensible distribution in my desired range, the LangevinIntegrator takes this temperature as an invitation to go fully of the charts, i. e. produces the distribution of what seems to be a much higher temperature than what the MetropolisMonteCarloIntegrator sampler interprets. Afterwards I struggle to interpret that and simulation outputs are hard to compare. That also makes it more difficult, for example, to combine Ensembler with the pyEMMA package.

I was wondering whether there was a way for units to be kept consistent across the package and at least to be documented.

Fix axis label of static_sim_plots

left panel: It should be indicated that the potential energy is in units of k_BT
middle panel: y-axis is currently labled with "simulation". Correct default label could be e.g. "probability density" or "frequency"
right panel: y-label "dhdpos" is not descriptive. Better label needed.

Reweighting

• new class: Reweighting potential, that saves action potential

• DHAM reweighting implementation

• DHAMed inclusion

• TRAM inclusion

• pathDHAM implementation

• Weber-Pande reweighting implementation

• Girsanov Reweighting implementation

cookie Cutter Repo

Make Repo more readable
simplify import structure for users

Build Release V.1.0 & Citeable

ToDos:

• Final Refactoring
• integrator unit tests
• Set correct default
• Build Release 1.0
• pip installation

Update ReadMe file

We need to update the ReadMe file

LongTerm Refactoring

• variable access rather than class vars

Additional Features

• make classes pickable
• Replica Parallelism

Beautify JNBs with Markdown text

• BasicPotentials
• BasicSimulations (@linkerst )
• EDS (@schroederb)
• ConveyorBelt (@dfhahn )
• FreeEnergy (Gerhard?/ @schroederb )
• EnhancedSampling (@linkerst )

Thermostats

• Anderson Thermostat - parameter check
• NoseHoover

Verlocity Verlet does not work for 2D potential

Execution of velocity verlet as well as leapfrog integration with 2D potentials leads to error: TypeError: can't multiply sequence by non-int of type 'float'. See example

`import os, sys
my_path = os.getcwd()+"/.."
print(my_path)
sys.path.append(my_path)

import numpy as np
from matplotlib import pyplot as plt

from ensembler.potentials.TwoD import wavePotential

from ensembler.integrator.newtonian import positionVerletIntegrator, velocityVerletIntegrator, leapFrogIntegrator
from ensembler.system import system

##Visualisation
from ensembler.visualisation.plotSimulations import static_sim_plots

sim_steps = 1000

#Simulation Setup
origpot = wavePotential(multiplicity=[2,2])

integrator = velocityVerletIntegrator()

sys=system(potential=origpot, integrator=integrator, position=np.array([20,70]), temperature=3)

#simulate
cur_state = sys.simulate(sim_steps, withdrawTraj=True, initSystem=False)

print("Trajectory length: ",len(sys.trajectory))
print()
print("last_state: ", cur_state)
print(len(sys.trajectory))
sys.trajectory.head()`

Implementation