bio-phys / pydhamed Goto Github PK

View Code? Open in Web Editor NEW

19.0 4.0 7.0 3.44 MB

Dynamic Histogram Analysis To Determine Free Energies and Rates from Biased Simulations

License: BSD 3-Clause "New" or "Revised" License

Jupyter Notebook 91.65% Python 8.34% Shell 0.01%

md-simulations science python free-energy free-energies kinetics kinetic-modeling molecular-dynamics

pydhamed's Introduction

PyDHAMed

DHAMed -Dynamic Histogram Analysis extended to detailed balance

Input are transition counts between states/bins and biases (if any). Biases specify the differences in the potential energy functions in the different simulation windows/runs.

To run DHAMed from a list of count matrices and an array specfying the biases in each simulation (window) is required.

To see how DHAMed can be used to extract free energies from biased simulations look at the example Jupyter notebook provided. https://github.com/bio-phys/PyDHAMed/blob/master/pydhamed/cg-rna/cg_RNA_duplex_formation.ipynb

Installation

To install PyDHAMed clone or download the repository

git clone https://github.com/bio-phys/PyDHAMed.git
cd PyDHAMed

Then install the dowloaded repository with pip:

pip install .

PyDHAMed is now ready for use.

Inputs

The list of the individual count matrices C contain the transition counts between the different states (or bins in umbrella sampling). C[i,j] where i is the product state and j the reactent state. The first row contains thus all the transitions into state 0.The first column C[:,0] all transition out of state 0.

The bias array contains a bias value for each state and for each simulation (or window in umbrella sampling). The bias array has the shape N rows nwin columns and contains the bias acting on each state in each simulation (window). The bias NEEDS to be given in units to kB_T.

Most parameters besides count_list and bias_ar are only relevant for testing and further code developement.

To run DHAMed

# import DHAMed functions
from pydhamed.optimize_dhamed import *
from pydhamed.determine_transition_counts import count_matrix

# determine transition counts for each trajectory
# Each frame in a trajectory needs to be assigned to one of the the n states
# of the system
for traj in traj_list:
  count_list.append(count_matrix(traj, n_states=n))

# Bias - need to specfiy the bias acting on each of the n states in the nwin simulation.
bias_ar = np.zeros((n, nwin))
for i in range(n)
    bias_ar[i,:] = np.loadtxt("bias"+i)

# run optimization
og = run_dhamed(count_list, bias_ar)

DHAMed examples

Two example calculations are provided in the pydhamed folder.

Ion channel permeation:

Umbrella sampling simulations of ion permeation through a channel. Data from all-atom simulations are analyzed in this example Jupyter notebook. https://github.com/bio-phys/PyDHAMed/blob/master/pydhamed/glic-ion-channel/glic_ion_channel_permeation.ipynb

RNA duplex formation:

Umbrella sampling simulations of RNA duplex formation using a coarse-grained model https://github.com/bio-phys/PyDHAMed/blob/master/pydhamed/cg-rna/cg_RNA_duplex_formation.ipynb

References

Dynamic Histogram Analysis To Determine Free Energies and Rates from biased Simulations, L. S. Stelzl, A. Kells, E. Rosta, G. Hummer, J. Chem. Theory Comput., 2017, http://pubs.acs.org/doi/abs/10.1021/acs.jctc.7b00373

pydhamed's People

Contributors

Stargazers

Watchers

Forkers

msultan kain88-de leelasd mandar5335 kbri-biology-experiment-computation moleorbitalhybridanalyst derienfe

pydhamed's Issues

Integration comparison with binless WHAM

*Enable users to run binless WHAM in the same framework as DHAMed.

last_g_zero=True sometimes gives an error

ZeroDivisionError                         Traceback (most recent call last)
<ipython-input-36-882a6db4a90e> in <module>()
----> 1 og = run_dhamed(c_l, -np.log(v_ar), g_init=og, numerical_gradients=False, last_g_zero=True)

/Users/lukas/Projects/DHAMed/public-dhamed/PyDHAMed/pydhamed/optimize_dhamed.pyc in grad_dhamed_likelihood_ref_0(g_prime, g, ip, jp, ti, tj, vi, vj, nk, nijp)
     88     grad = np.zeros(g.shape[0] )
     89     grad[:-1]  += nk[:-1]
---> 90     grad = loop_grad_dhamed_likelihood_0(grad,g, ip, jp, ti, tj, vi, vj, nijp)
     91     return grad[:-1]
     92 

ZeroDivisionError: division by zero

Estimate uncertainty

Implement uncertainty estimation
Benchmark for model system

Generating Bias Matrix

I have performed an Umbrella Sampling analysis on a membrane system in gromacs. My collective variable is the distance of a small molecule from the center of a membrane. Using the pullf file generated from gromacs how would I go about generating the bias force matrix for each of my 80 US windows?

Kinetics calculations

To do:

Calculate rates trajectories where dynamics effectively evolve with the "true" dynamics (E.g., sets of short trajectories).
Calculate rates from biased dynamics (e.g., umbrella sampling)...

Add GLIC example

Add GLIC ion-permeation channel example Jupyter.

Handle disconnected states

Transition counts and bias potentials of states not properly connect to the other states need to be excluded in a consistent fashion.
Will be important for DHAMed calculations for complex systems.

Continuous Integration

Use Travis CI to run tests automatically :-)

Numba error

Dear Lukas,
thanks a lot for providing the software to DHAMed.
I am trying to use the program with the following test data:

traj=np.array([0., 0., 1., 2., 1., 0., 2., 2.])
c_l = [count_matrix(traj,n_states=3)]
v_ar=np.array([1, 1.5, 2]).reshape(3,1)

I get the following error:

File "dhamed.py", line 17, in
bl_g = run_dhamed(c_l, -np.log(v_ar))
....
....
InternalError: unsupported array index type float64 in [float64]
[1] During: typing of intrinsic-call at /home/andt88/miniconda2/lib/python2.7/site-packages/pydhamed/optimize_dhamed.py (21)
--%<-----------------------------------------------------------------

File "../../../../miniconda2/lib/python2.7/site-packages/pydhamed/optimize_dhamed.py", line 21

Bias potential matrix

Hi,
I wanted to make sure that I understand the definition of the bias potential matrix correctly.
Assuming an umbrella sampling experiment with 4 windows A, B, C, D and 3 states 0, 1, 2. The restraining potential defined as k*(x-x0)^2. Should the cell A0 contain the following?
x=median value of bin 0,
x0=target value of window A

Division by zero error

Dear Lukas

I am trying to use PyDHAMed to compute a PMF for a small molecule crossing a membrane bilayer. I performed umbrella sampling in 1.5Å windows, spanning approx. 7nm (57 umbrellas in total).

Following your code, I first created the list of transition matrices for each of the umbrellas:

def create_transition_matrix(traj: npt.NDArray, microstates: npt.NDArray)->npt.NDArray:
    """
    Parameters
    --------------
    traj: np.array of size n
    microstates: np.array of size m

    First, the trajectory is binned into the provided
    microstates. Next, a transition matrix is calculated
    giving the number of transitions between each of the
    microstates.
    
    Returns: m x m matrix with the transitions between 
    all possible pairs of bins in subsequent trajectory frames. 
    """
    traj_binned = np.digitize(traj, microstates)
    n_states = microstates.shape[0]+1 # digitize is 1 indexed.
    mat = np.zeros((n_states, n_states))
    for (x, y), c in six.iteritems(Counter(zip(traj_binned[:-1], traj_binned[1:]))):
        mat[int(y), int(x)] = c
    return mat

colvar_files = glob.glob("../gromacs/plumed_us_*-colvar.dat")

transition_matrices = list()
for i, f in enumerate(colvar_files):
    print("Processing trajectory: ", f)
    traj = np.loadtxt(f)[500:, 1]
    transition_matrices.append(create_transition_matrix(traj, microstates)[1:, 1:])
with open("NFP_membrane_transition_counts.pickle", "wb") as outfile:
    pickle.dump(transition_matrices, outfile)

Next, I calculated the bias grid (in kBT units) as follows:

def create_biasses(q0, kappa, microstates)->npt.NDArray:
    bias = (0.5 * force_constants[:, np.newaxis] * (centers[:, np.newaxis] - microstates)**2).T
    # convert to kBT
    return (bias*1000)/(R*T)

#col0 = bias number
#col1 = bias center position
#col2 = bias force constant

umbrellas = np.loadtxt("../gromacs/wham_input.txt", skiprows=1)
centers = umbrellas[:, 1]
force_constants = umbrellas[:, 2]

biasses = create_biasses(centers, force_constants, microstates)
np.save("NFP_membrane_biasses.npy", biasses)

, where the input file simply contains the centers and force constants of the harmonic potentials. But when I try to use run_dhamedI cannot seem to get sensible results. When I use jit_gradient=False, the calculation finished but I do not get a sensible result (see figure).

When I use jit_gradient=True, I get a division by zero error at '_loop_grad_dhamed_likelihood_0'. Would you have an idea of what is going wrong?

Best,
Warre