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.

ensembler's Issues

Update ReadMe file

We need to update the ReadMe file

Build Release V.1.0 & Citeable

ToDos:

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

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.

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()`

Integrators

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

Bring RE to Ensembler

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

Refactoring

• manage typing via util/ensemblerTypes
• System: init_Velocity

Enhanced sampling

• Umbrella sampling

• Metadynamics

• Add custom potential

• aMD?

Conclusion

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

Implementation

Beautify JNBs with Markdown text

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

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

Additional Features

• make classes pickable
• Replica Parallelism

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

Needed to fix this manually after construction.

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.

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.

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

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)

CI integration

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

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()

Thermostats

• Anderson Thermostat - parameter check
• NoseHoover

LongTerm Refactoring

• variable access rather than class vars

Introduction

cookie Cutter Repo

Make Repo more readable
simplify import structure for users

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).