Trp-cage haMSM Analysis

[1]:

import matplotlib.pyplot as plt

Load haMSM

We’ll use the haMSM we created in the construction examples.

[2]:

import pickle

with open('data/pickled_model', 'rb') as inf:
    model = pickle.load(inf)

Steady-state

[3]:

cluster_centers = model.targetRMSD_centers

plt.scatter(cluster_centers, model.pSS, color='k')

# These have NaN entries in targetRMSD_centers, so we set them manually using the bounds
target_pcoord = model.target_pcoord_bounds[0,0]
basis_pcoord = model.basis_pcoord_bounds[0,1]

plt.scatter(target_pcoord, model.pSS[model.indTargets], color='r')
plt.scatter(basis_pcoord, model.pSS[model.indBasis], color='b')


plt.yscale('log')
plt.ylabel('Probability')
plt.xlabel('MSM State Average Progress Coordinate')

[3]:

Text(0.5, 0, 'MSM State Average Progress Coordinate')

Committors

We can compute the committor for each MSM microstate, which is the probability that a trajectory initiated from that state will reach the target before reaching the basis.

[4]:

model.get_committor()

[10/04/22 10:35:10] WARNING  Note that, if steady-state weighted ensemble data is being analyzed,    msm_we.py:7094
                             this is a 'pseudocommittor' and not a true committor as a result of                   
                             being constructed from a one-way ensemble.

It’s often interesting to look at the relationship between our progress coordinate and the committor.

Here, we see that our progress coordinate is degenerate in the committor, which suggests it’s not a great progress coordinate.

[5]:

plt.scatter(cluster_centers, model.q)

plt.axvline(model.target_pcoord_bounds[0,0], color='b')
plt.axvline(model.basis_pcoord_bounds[0,1], color='r')

plt.ylabel('Committor')
plt.xlabel('MSM State Average Progress Coordinate')

[5]:

Text(0.5, 0, 'MSM State Average Progress Coordinate')

Flux profiles

[6]:

model.get_flux()

We can plot flux profiles for the full-data model as well as the validation models. The flatness of this flux profile is a proxy for how converged your WE is.

We can visualize it in progress-coordinate space…

[7]:

ax = model.plot_flux()

ax.set_xlim([0,0.8])
ax.set_ylim([5e3, 5e7])

[10/04/22 10:35:11] WARNING  Fluxes have not yet been generated for validation_model_0, generating   msm_we.py:6309
                             now.

                    WARNING  Fluxes have not yet been generated for validation_model_1, generating   msm_we.py:6309
                             now.

[7]:

(5000.0, 50000000.0)

Or in committor-space.

[8]:

ax = model.plot_flux_committor()
ax.set_ylim([5e3, 5e7])

[10/04/22 10:35:12] WARNING  Committor-fluxes have not yet been generated for main_model, generating msm_we.py:6159
                             now.

[10/04/22 10:35:13] WARNING  Committors have not yet been generated for validation_model_0,          msm_we.py:6153
                             generating now.

                    WARNING  Note that, if steady-state weighted ensemble data is being analyzed,    msm_we.py:7094
                             this is a 'pseudocommittor' and not a true committor as a result of                   
                             being constructed from a one-way ensemble.

                    WARNING  Committor-fluxes have not yet been generated for validation_model_0,    msm_we.py:6159
                             generating now.

[10/04/22 10:35:14] WARNING  Committors have not yet been generated for validation_model_1,          msm_we.py:6153
                             generating now.

                    WARNING  Note that, if steady-state weighted ensemble data is being analyzed,    msm_we.py:7094
                             this is a 'pseudocommittor' and not a true committor as a result of                   
                             being constructed from a one-way ensemble.

                    WARNING  Committor-fluxes have not yet been generated for validation_model_1,    msm_we.py:6159
                             generating now.

                    WARNING  Note that, if steady-state weighted ensemble data is being analyzed,    msm_we.py:7094
                             this is a 'pseudocommittor' and not a true committor as a result of                   
                             being constructed from a one-way ensemble.

[8]:

(5000.0, 50000000.0)

Cluster structures

[9]:

model.update_cluster_structures()

Let’s get the coordinates of a structure that got assigned to MSM microstate 3.

[10]:

model.cluster_structures[3][0].shape

[10]:

(272, 3)

[ ]: