<jupyter_start><jupyter_text>Characterizing the metastable states of Alanine Dipeptide Approach based on the identification of the local minima of the Free Energy Surface.In this notebook we will use the DeepTICA collective variables of Alanine Dipeptide to characterize its metastable states. We shall look for the local minima of the free energy landscape to identify and label the different metastable states<jupyter_code>#Standard imports
import matplotlib.pyplot as plt
import numpy as np
from itertools import cycle
#MD-related imports
from stateinterpreter import load_trajectory, prepare_training_dataset,Classifier
from stateinterpreter.metastable import identify_metastable_states
from stateinterpreter.utils.io import load_dataframe
from stateinterpreter.utils.plot import plot_classifier_complexity_vs_accuracy,paletteFessa,plot_states
from stateinterpreter.utils.metrics import get_best_reg
from stateinterpreter.descriptors import compute_descriptors
colvar = load_dataframe('alanine/COLVAR_DeepTICA') # Load collective variable file
ax = colvar.plot.scatter('phi','psi',s=0.5, alpha=0.1, c='k')
#Plot Styling
ax.margins(0)
ax.set_xlabel(r"$\phi$ [rad]")
ax.set_ylabel(r"$\psi$ [rad]")
ax.set_title("Ramachandran plot for Alanine Dipeptide")
ax.set_xticks(np.arange(-np.pi, np.pi+np.pi/2, step=(np.pi/2)))
ax.set_xticklabels(['-π','-π/2','0','π/2','π'])
ax.set_yticks(np.arange(-np.pi, np.pi+np.pi/2, step=(np.pi/2)))
ax.set_yticklabels(['-π','-π/2','0','π/2','π'])
plt.show()<jupyter_output><empty_output><jupyter_text>DescriptorsWe now show how to use the helper functions implemented in `stateinterpreter` to calculate the physical descriptors starting from a molecular dynamics trajectory. In the present case we compute the H-bonds contacts and the dihedral angles.<jupyter_code>traj_dict = {
    'trajectory' : 'alanine/alanine.xtc',
    'topology' : 'alanine/input.ala2.pdb'
}

traj = load_trajectory(traj_dict)
descriptors, feats_info = compute_descriptors( traj[::5], descriptors = ['hbonds_contacts','dihedrals'] )
descriptors.head()<jupyter_output>Computing Hydrogen bonds contacts
Donors: [6, 16]
Acceptors: [ 5  6 15 16]
Descriptors: (10001, 10)<jupyter_text>Local minima of the free energy surfaceWe now use the utilities implemented in `stateinterpreter.metastable` to automatically identify the local minima of the free energy surface.<jupyter_code>bandwidth = 0.1 # Set the bandwidth for the Gaussian kernel used in the KDE approximation
cvs = colvar[::5] # Select only a subset of configurations
cv_list = ['DeepTICA 1','DeepTICA 2'] # List of collective variables to use for clustering
#Bias correction
kbt = 2.5
logweights = cvs['opes.bias'].values/kbt

states_labels = identify_metastable_states(
    cvs,
    cv_list,
    kbt,
    bandwidth,
    logweights=logweights,
    gradient_descent_iterates = 1,
)
plot_states(
    cvs, 
    states_labels, 
    cv_list, 
    fes_isolines = True, 
    n_iso_fes = 9, 
    ev_iso_labels = 2, 
    alpha=0.3, 
    cmap_name = 'Set2', 
    save_folder=None, 
    axs = None, 
    bandwidth = bandwidth
)
print(f"Example: first 5 labels for the states of the identified metastable states.")
states_labels.head()<jupyter_output>Example: first 5 labels for the states of the identified metastable states.<jupyter_text>Metastable states characterization based on dihedral angles<jupyter_code>feat_type = 'angles' #Select feature type
select_feat = 'sin_|cos_' #Regex to filter angle
desc = descriptors.filter(regex=select_feat)
print('Example of the selected descriptors')
desc.head()
# Init classifier
regularizers = np.geomspace(0.1, 1, 51)
num_samples = 100 # number of configurations used to train the classifier

sample_obj, features_names = prepare_training_dataset(
    descriptors=desc, 
    states_labels = states_labels, 
    n_configs= num_samples
) 
classifier = Classifier(sample_obj, features_names)
classifier.compute(regularizers, max_iter = 10000) # Train the classifier

# Extract optimal regularization parameter
reg,acc,num = get_best_reg(classifier)
print(f'Best log_10 (lambda) : {np.round(np.log10(reg),3)}')
print(f'Accuracy        : {np.round(acc*100,0)} %')
print(f'No. features    : {num}')

# Count states (for plotting purposes)
num_states = len(classifier.classes)
num_histo = num_states if num_states>2 else 1

# Classes names (for plotting purposes)
relevant_feat = classifier.feature_summary(reg)
classes_names = classifier.classes

# Init Plot
fig, ax = plt.subplots(1,1,figsize=(3,3),dpi=100)

# Init Colors 
colors = cycle(iter(paletteFessa[::-1]))
color = [next(colors) for i in range(5)]

# Plot regularization paths.
_, ax_twin = plot_classifier_complexity_vs_accuracy(classifier, ax = ax)
ax.axvline(np.log10(reg),linestyle='dotted',color='k')
ax.text(np.log10(reg)+0.02,acc+.05,f'{np.round(acc*100,1)}%',fontsize=12,color='fessa0')
ax_twin.text(np.log10(reg)+0.02,num+0.5,f'{num}',fontsize=12,color='fessa1')

#Print relevant features
classifier.print_selected(reg)
plt.show()<jupyter_output>Optimizing Lasso Estimator: 100%|██████████| 51/51 [00:00<00:00, 1293.94it/s]
{'fraction_non_alphanumeric': 0.11220472440944881, 'fraction_numerical': 0.025196850393700787, 'mean_word_length': 5.9986225895316805, 'pattern_counts': {'":': 0, '<': 15, '<?xml version=': 0, '>': 15, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '4', 'dir': 'jupyter-structured-clean-dedup', 'id': '102101', 'license': 'MIT', 'n_tokens_mistral': 1817, 'n_tokens_neox': 1693, 'n_words': 435, 'path': 'tutorials/3_FES_minima_classification.ipynb', 'repo_name': 'luigibonati/md-stateinterpreter', 'size': '5080'}