<jupyter_start><jupyter_text>In the following we use the same approach as in `1-RMSD_clustering_with_IPA_share`  References and trajectories Load the reference open and closed crystal structures. Note that 6 versions of the closed form are loaded, they are identical but with rotating chain IDs so that it doesn't matter how open structure closes, the minimum RMSD value will be retained.<jupyter_code>rf0=pt.load('../1-Trajectories/0-Parm_rst/closed6_dry_rot0.rst7', 
            '../1-Trajectories/0-Parm_rst/closed6_dry_rot0.parm7')

rf1=pt.load('../1-Trajectories/0-Parm_rst/closed6_dry_rot1.rst7', 
            '../1-Trajectories/0-Parm_rst/closed6_dry_rot1.parm7')

rf2=pt.load('../1-Trajectories/0-Parm_rst//closed6_dry_rot2.rst7',
            '../1-Trajectories/0-Parm_rst/closed6_dry_rot2.parm7')

rf3=pt.load('../1-Trajectories/0-Parm_rst/closed6_dry_rot3.rst7', 
            '../1-Trajectories/0-Parm_rst/closed6_dry_rot3.parm7')

rf4=pt.load('../1-Trajectories/0-Parm_rst/closed6_dry_rot4.rst7',
            '../1-Trajectories/0-Parm_rst/closed6_dry_rot4.parm7')

rf5=pt.load('../1-Trajectories/0-Parm_rst/closed6_dry_rot5.rst7', 
            '../1-Trajectories/0-Parm_rst/closed6_dry_rot5.parm7')

rfopen=pt.load('../1-Trajectories/0-Parm_rst/open10_dry.rst7',
               '../1-Trajectories/0-Parm_rst/open10_dry.parm7')
indir1 = '../1-Trajectories/1-NoEq_without_IPA'
topfile = '../1-Trajectories/0-Parm_rst/open10_dry.parm7'
traj_list = sorted(glob(indir1+'/*.nc'))
traj_list<jupyter_output><empty_output><jupyter_text>GMM construction with limited dataset First the GMM is built using uncorrelated data. We assume the data to be uncorrelated for a frame every 500 ps (dt=5) The trajectory is therefore first loaded with a stride of 5<jupyter_code>traj = pt.iterload(traj_list, topfile, stride=5)<jupyter_output><empty_output><jupyter_text>The RMSD is then calculated. Because of the symmerty and chain IDs, 6 references are used for the closed form (simplky corresponding to a rotation of the chain IDs, and the minimum RMSD is retained<jupyter_code>rmsd0=pt.rmsd(traj,ref=rf0, mask="@CA,N,C,O")
rmsd1=pt.rmsd(traj,ref=rf1, mask="@CA,N,C,O")
rmsd2=pt.rmsd(traj,ref=rf2, mask="@CA,N,C,O")
rmsd3=pt.rmsd(traj,ref=rf3, mask="@CA,N,C,O")
rmsd4=pt.rmsd(traj,ref=rf4, mask="@CA,N,C,O")
rmsd5=pt.rmsd(traj,ref=rf5, mask="@CA,N,C,O")

rmsd_closed=np.nanmin(np.array([rmsd0,rmsd1,rmsd2,rmsd3,rmsd4,rmsd5]), axis=0)

rmsd_open=pt.rmsd(traj,ref=rfopen, mask="@CA,N,C,O")<jupyter_output><empty_output><jupyter_text>We construct a array Z with the RMSD to the closed form and to the open form for each frame of the reduced trajectory<jupyter_code>Z=[]
for i in range(0,len(rmsd_open)):
    a=np.array([rmsd_closed[i],rmsd_open[i]])
    Z.append(a)
Z=np.array(Z)
# with open('./pickled/2-Without_IPA/Z_dt5.pickle', 'wb') as outf:
#     pickle.dump(Z, outf)
    
with open('./pickled/2-Without_IPA/Z_dt5.pickle', 'rb') as inf:
    Z = pickle.load(inf)<jupyter_output><empty_output><jupyter_text>We need to find the optimum number of components for constructing the GMM so we screen for a window of 2 to 25 components.<jupyter_code>fec = FEC.FreeEnergyClustering(Z, min_n_components=5, max_n_components=25, temperature=293.15,
                               x_lims=[[0.0,4.5],[0.0,4.5]],
                               n_iterations=20, n_grids=151, n_splits=1,stack_landscapes=False)
coords, FE_landscape, FE_points = fec.landscape()<jupyter_output><empty_output><jupyter_text>The optimum number of components is 10, so we build and save a GMM with that many components<jupyter_code>fec = FEC.FreeEnergyClustering(Z, min_n_components=10, max_n_components=10, temperature=293.15,
                               x_lims=[[0.0,4.5],[0.0,4.5]],
                               n_iterations=20, n_grids=151, n_splits=1,stack_landscapes=False)
coords, FE_landscape, FE_points = fec.landscape()
# with open('./pickled/2-Without_IPA/coords_dt5.pickle', 'wb') as outf:
#     pickle.dump(coords, outf)
# with open('./pickled/2-Without_IPA/FE_landscape_dt5.pickle', 'wb') as outf:
#     pickle.dump(FE_landscape, outf)
# with open('./pickled/2-Without_IPA/FE_points_dt5.pickle', 'wb') as outf:
#     pickle.dump(FE_points, outf)
    
with open('./pickled/2-Without_IPA/coords_dt5.pickle', 'rb') as inf:
    coords = pickle.load(inf)
with open('./pickled/2-Without_IPA/FE_landscape_dt5.pickle', 'rb') as inf:
    FE_landscape = pickle.load(inf)
with open('./pickled/2-Without_IPA/FE_points_dt5.pickle', 'rb') as inf:
    FE_points = pickle.load(inf)

fec.FE_points_ = np.copy(FE_points)
fec.FE_landscape_ = np.copy(FE_landscape)
fec.coords_ = coords<jupyter_output><empty_output><jupyter_text>We can then plot the FES<jupyter_code>fec.visualize(fontsize=12, savefig=True, xlabel=r'RMSD to closed crystal structure ($\AA$)',
              ylabel=r'RMSD to open crystal structure ($\AA$)', zlabel='z', title=None,
              vmax=4.5, colormap='nipy_spectral', lim=0.1, n_contour_levels=100, show_data=False, 
              figsize= [5, 5], transition_point_size=1, core_point_size=5, core_point_edgecolor='k', 
              colormap2='spectral',clust_point_edgecolor='w', clust_point_facecolor='', 
              clust_point_size=5, clust_point_linewidth=4, clust_point_marker='s',legend=None,
              extend=None, rmsd_point_size=40, rmsd_point_fontsize=10,rmsd_point_color='k', dx=1,
              filename='./plots/2-Without_IPA/1-FES_with_IPA_dt5', ax=None)<jupyter_output>/home/el14718/Modelling/10-ALKEIG/12-Analysis_4/11-InfleCS_clustering_3/free_energy_clustering/GMM_free_energy.py:537: MatplotlibDeprecationWarning: You are modifying the state of a globally registered colormap. In future versions, you will not be able to modify a registered colormap in-place. To remove this warning, you can make a copy of the colormap first. cmap = copy.copy(mpl.cm.get_cmap("nipy_spectral"))
  my_cmap.set_over('white')<jupyter_text>Clustering for the limited dataset<jupyter_code>labels, cluster_centers = fec.cluster(coords, FE_points,Z, assign_transition_points=True)
print("Cluster center indices: "+str(cluster_centers))
# with open('./pickled/2-Without_IPA/labels_dt5.pickle', 'wb') as outf:
#     pickle.dump(labels, outf)
# with open('./pickled/2-Without_IPA/cluster_centers_dt5.pickle', 'wb') as outf:
#     pickle.dump(cluster_centers, outf)
    
with open('./pickled/2-Without_IPA/labels_dt5.pickle', 'rb') as inf:
    labels=pickle.load(inf)
with open('./pickled/2-Without_IPA/cluster_centers_dt5.pickle', 'rb') as inf:
    cluster_centers=pickle.load(inf)
    
fec.labels_ = np.copy(labels)
fec.cluster_centers_ = np.copy(cluster_centers)
fec.pathways_ = None
fec.visualize(fontsize=12, savefig=True, xlabel=r'RMSD to closed crystal structure ($\AA$)',
              ylabel=r'RMSD to open crystal structure ($\AA$)', zlabel='z', title=None,
              vmax=4.5, colormap='nipy_spectral', lim=0.1, n_contour_levels=100, show_data=True, 
              figsize= [5, 5], transition_point_size=0, core_point_size=0.05, core_point_edgecolor=None, 
              colormap2='viridis',clust_point_edgecolor=None, clust_point_facecolor='k', 
              clust_point_size=30, clust_point_linewidth=0, clust_point_marker='o',legend=None,
              extend='neither', rmsd_point_size=40, rmsd_point_fontsize=10, dx=1,
              core_point_alpha=0.8, transition_point_alpha=0.0, surface_alpha=0.0,
              filename='./plots/2-Without_IPA/2-Cluster_map_with_IPA_dt5.png', ax=None)
fec.visualize(fontsize=12, savefig=True, xlabel=r'RMSD to closed crystal structure ($\AA$)',
              ylabel=r'RMSD to open crystal structure ($\AA$)', zlabel='z', title=None,
              vmax=4.5, colormap='nipy_spectral', lim=0.1, n_contour_levels=100, show_data=True, 
              figsize= [5, 5], transition_point_size=0, core_point_size=0.0, core_point_edgecolor='k', 
              colormap2='jet',clust_point_edgecolor='k', clust_point_facecolor='w', 
              clust_point_size=30, clust_point_linewidth=0, clust_point_marker='o',legend=None,
              extend='neither', rmsd_point_size=40, rmsd_point_fontsize=10, dx=1,
              core_point_alpha=0.8, transition_point_alpha=0.0, surface_alpha=1,rmsd_point_color='k',
              filename='./plots/2-Without_IPA/3-FES_with_IPA_with_cluster_centres_dt5.png', ax=None)<jupyter_output><empty_output><jupyter_text>To use the red colour code<jupyter_code>hexlist=[]

cmap = cm.get_cmap("Blues",8)    # PiYG

for i in range(cmap.N):
    rgb = cmap(i)[:3] # will return rgba, we take only first 3 so we get rgb
    hexlist.append(matplotlib.colors.rgb2hex(rgb))
# sns.set_palette(hexlist)
# sns.palplot(sns.color_palette())

my_cmap=ListedColormap([hexlist[7],hexlist[6],hexlist[4],hexlist[5],hexlist[2],hexlist[3]])
fec.pathways_ = None
fig, ax = plt.subplots()
fec.visualize(fontsize=12, savefig=False, xlabel=r'RMSD to closed crystal structure ($\AA$)',
              ylabel=r'RMSD to open crystal structure ($\AA$)', zlabel='z', title=None,
              vmax=4.5, colormap='nipy_spectral', lim=0.1, n_contour_levels=100, show_data=True, 
              figsize= [5, 5], transition_point_size=0, core_point_size=0.25, core_point_edgecolor=None, 
              colormap2=my_cmap,clust_point_edgecolor=None, clust_point_facecolor=None, 
              clust_point_size=0, clust_point_linewidth=0, clust_point_marker='s',legend=None,
              extend='neither', rmsd_point_size=0, rmsd_point_fontsize=0, dx=1,
              core_point_alpha=1, transition_point_alpha=0.0, surface_alpha=0.0,
              filename='./plots/1-With_IPA/2-Cluster_map_with_IPA_dt1', ax=ax,core_point_marker='o')


ax.text(1.3,3.1,r'$\mathcal{S}_1$',c='k',fontsize=12)
ax.text(2.2,3.1,r'$\mathcal{S}_2$',c='k',fontsize=12)
ax.text(2.28,2.55,r'$\mathcal{S}_3$',c='k',fontsize=12)
ax.text(3.4,2.9,r'$\mathcal{S}_4$',c='k',fontsize=12)
ax.text(3.,2.5,r'$\mathcal{S}_5$',c='w',fontsize=12 )
ax.text(2.8,1.9,r'$\mathcal{S}_6$',c='w',fontsize=12)

fig.savefig('./plots/2-Without_IPA/2-Cluster_map_with_IPA_dt5.png', dpi=600)<jupyter_output><empty_output><jupyter_text>We can then look at the population of each state for the limited dataset<jupyter_code># Computing state populations
state_populations = fec.population_states(n_sampled_points=100000)
state_populations*100<jupyter_output>Sampling points from density.
Assigning cluster labels based on precomputed density model clustering.
Computing state populations.<jupyter_text>Assigning clusters to the full dataset The trajectory is therefore first loaded with a stride of 1<jupyter_code>traj_all = pt.iterload(traj_list, topfile, stride=1)<jupyter_output><empty_output><jupyter_text>The RMSD is calculated and the Z array constructed<jupyter_code>rmsd0=pt.rmsd(traj_all,ref=rf0, mask="@CA,N,C,O")
rmsd1=pt.rmsd(traj_all,ref=rf1, mask="@CA,N,C,O")
rmsd2=pt.rmsd(traj_all,ref=rf2, mask="@CA,N,C,O")
rmsd3=pt.rmsd(traj_all,ref=rf3, mask="@CA,N,C,O")
rmsd4=pt.rmsd(traj_all,ref=rf4, mask="@CA,N,C,O")
rmsd5=pt.rmsd(traj_all,ref=rf5, mask="@CA,N,C,O")

rmsd_closed_all=np.nanmin(np.array([rmsd0,rmsd1,rmsd2,rmsd3,rmsd4,rmsd5]), axis=0)

rmsd_open_all=pt.rmsd(traj_all,ref=rfopen, mask="@CA,N,C,O")
Z_all=[]
for i in range(0,len(rmsd_open_all)):
    a=np.array([rmsd_closed_all[i],rmsd_open_all[i]])
    Z_all.append(a)
Z_all=np.array(Z_all)
# with open('./pickled/2-Without_IPA/Z_dt1.pickle', 'wb') as outf:
#     pickle.dump(Z_all, outf)
    
with open('./pickled/2-Without_IPA/Z_dt1.pickle', 'rb') as inf:
    Z_all = pickle.load(inf)<jupyter_output><empty_output><jupyter_text>Then we can assign all the points to the pre-computed clusters<jupyter_code>labels_all = fec.evaluate_clustering(Z_all, assign_transition_points=True)
# with open('./pickled/2-Without_IPA/labels_all.pickle', 'wb') as outf:
#     pickle.dump(labels_all, outf)

with open('./pickled/2-Without_IPA/labels_all.pickle', 'rb') as inf:
    labels_all=pickle.load(inf)
    
# fec.labels_ = np.copy(labels_all)
fec.labels_ = np.copy(labels_all)
# Computing state populations
state_populations = fec.population_states(n_sampled_points=100000)
state_populations*100<jupyter_output>Sampling points from density.
Assigning cluster labels based on precomputed density model clustering.
Computing state populations.<jupyter_text>FES for the full dataset<jupyter_code>fec_all = FEC.FreeEnergyClustering(Z_all, min_n_components=10, max_n_components=10, temperature=293.15,
                               x_lims=[[0.0,4.5],[0.0,4.5]],
                               n_iterations=5, n_grids=150, n_splits=1,stack_landscapes=False)
coords_all, FE_landscape_all, FE_points_all = fec_all.landscape()
# with open('./pickled/2-Without_IPA/coords_dt1.pickle', 'wb') as outf:
#     pickle.dump(coords_all, outf)
# with open('./pickled/2-Without_IPA/FE_landscape_dt1.pickle', 'wb') as outf:
#     pickle.dump(FE_landscape_all, outf)
# with open('./pickled/2-Without_IPA/FE_points_dt1.pickle', 'wb') as outf:
#     pickle.dump(FE_points_all, outf)
    
with open('./pickled/2-Without_IPA/coords_dt1.pickle', 'rb') as inf:
    coords_all = pickle.load(inf)
with open('./pickled/2-Without_IPA/FE_landscape_dt1.pickle', 'rb') as inf:
    FE_landscape_all = pickle.load(inf)
with open('./pickled/2-Without_IPA/FE_points_dt1.pickle', 'rb') as inf:
    FE_points_all = pickle.load(inf)

fec_all.FE_points_ = np.copy(FE_points_all)
fec_all.FE_landscape_ = np.copy(FE_landscape_all)
fec_all.coords_ = coords_all 
fec_all.visualize(fontsize=12, savefig=True, xlabel=r'RMSD to closed crystal structure ($\AA$)',
              ylabel=r'RMSD to open crystal structure ($\AA$)', zlabel='z', title=None,
              vmax=4.5, colormap='nipy_spectral', lim=0.1, n_contour_levels=100, show_data=False, 
              figsize= [5, 5], transition_point_size=1, core_point_size=5, core_point_edgecolor='k', 
              colormap2='spectral',clust_point_edgecolor='w', clust_point_facecolor='', 
              clust_point_size=5, clust_point_linewidth=4, clust_point_marker='s',legend=None,
              extend=None, rmsd_point_size=40, rmsd_point_fontsize=10,rmsd_point_color='k', dx=1,
              filename='./plots/2-Without_IPA/1-FES_with_IPA_dt1', ax=None)
fec_all.labels_ = np.copy(labels_all)

## the cluster centers are frame indexes so they are not valid here
# fec_all.cluster_centers_ = np.copy(cluster_centers)
fec_all.pathways_ = None
fig, ax = plt.subplots()
fec_all.visualize(fontsize=12, savefig=False, xlabel=r'RMSD to closed crystal structure ($\AA$)',
              ylabel=r'RMSD to open crystal structure ($\AA$)', zlabel='z', title=None,
              vmax=4.5, colormap='nipy_spectral', lim=0.1, n_contour_levels=100, show_data=True, 
              figsize= [5, 5], transition_point_size=0, core_point_size=0.25, core_point_edgecolor=None, 
              colormap2=my_cmap,clust_point_edgecolor=None, clust_point_facecolor=None, 
              clust_point_size=0, clust_point_linewidth=0, clust_point_marker='s',legend=None,
              extend='neither', rmsd_point_size=0, rmsd_point_fontsize=0, dx=1,
              core_point_alpha=1, transition_point_alpha=0.0, surface_alpha=0.0,
              filename='./plots/1-With_IPA/2-Cluster_map_with_IPA_dt1', ax=ax,core_point_marker='o')


ax.text(1.3,3.1,r'$\mathcal{S}_1$',c='k',fontsize=14)
ax.text(2.2,3.1,r'$\mathcal{S}_2$',c='k',fontsize=14)
ax.text(2.27,2.56,r'$\mathcal{S}_3$',c='k',fontsize=14)
ax.text(3.4,2.9,r'$\mathcal{S}_4$',c='w',fontsize=14)
ax.text(3.,2.5,r'$\mathcal{S}_5$',c='w',fontsize=14 )
ax.text(2.8,1.9,r'$\mathcal{S}_6$',c='w',fontsize=14)

fig.savefig('./plots/2-Without_IPA/2-Cluster_map_with_IPA_dt1.png', dpi=600)<jupyter_output><empty_output><jupyter_text>Creating crisp clustering for extracting representative structures Note that the GMM needs to be reconstructed form scratch running the following command, else the clustering doesn't work<jupyter_code>fec = FEC.FreeEnergyClustering(Z, min_n_components=10, max_n_components=10, temperature=293.15,
                               x_lims=[[0.0,4.5],[0.0,4.5]],
                               n_iterations=5, n_grids=151, n_splits=1,stack_landscapes=False)
coords, FE_landscape, FE_points = fec.landscape()
with open('./pickled/2-Without_IPA/coords_dt5.pickle', 'rb') as inf:
    coords = pickle.load(inf)
with open('./pickled/2-Without_IPA/FE_landscape_dt5.pickle', 'rb') as inf:
    FE_landscape = pickle.load(inf)
with open('./pickled/2-Without_IPA/FE_points_dt5.pickle', 'rb') as inf:
    FE_points = pickle.load(inf)

fec.FE_points_ = np.copy(FE_landscape)
fec.FE_landscape_ = np.copy(FE_landscape)
fec.coords_ = coords  
labels, cluster_centers = fec.cluster(coords,FE_points,Z, assign_transition_points=True)
print("Cluster center indices: "+str(cluster_centers))
# with open('./pickled/2-Without_IPA/labels_dt5_crisp.pickle', 'wb') as outf:
#     pickle.dump(labels, outf)
# with open('./pickled/2-Without_IPA/cluster_centers_dt5_crisp.pickle', 'wb') as outf:
#     pickle.dump(cluster_centers, outf)
    
with open('./pickled/2-Without_IPA/labels_dt5_crisp.pickle', 'rb') as inf:
    labels=pickle.load(inf)
with open('./pickled/2-Without_IPA/cluster_centers_dt5_crisp.pickle', 'rb') as inf:
    cluster_centers=pickle.load(inf)
    
fec.labels_ = np.copy(labels)
fec.cluster_centers_ = np.copy(cluster_centers)
state_populations = fec.population_states(n_sampled_points=10000)
state_populations
from collections import Counter
Counter(labels)
# Just in case
labels_old=labels
cluster_centers_old=cluster_centers
# To retry
labels=labels_old
cluster_centers=cluster_centers_old
i=1
threshold=1
cluster_FE = FE_points[labels==i]
cluster_center_FE = FE_points[int(cluster_centers[i-1])]
diff_FE = cluster_FE-cluster_center_FE
new_labels = np.zeros(cluster_FE.shape[0])
new_labels[diff_FE< threshold] = i

labels[labels==i] = np.copy(new_labels)
i=2
threshold=1
cluster_FE = FE_points[labels==i]
cluster_center_FE = FE_points[int(cluster_centers[i-1])]
diff_FE = cluster_FE-cluster_center_FE
new_labels = np.zeros(cluster_FE.shape[0])
new_labels[diff_FE< threshold] = i

labels[labels==i] = np.copy(new_labels)
i=3
threshold=0.35
cluster_FE = FE_points[labels==i]
cluster_center_FE = FE_points[int(cluster_centers[i-1])]
diff_FE = cluster_FE-cluster_center_FE
new_labels = np.zeros(cluster_FE.shape[0])
new_labels[diff_FE< threshold] = i

labels[labels==i] = np.copy(new_labels)
i=4
threshold=3
cluster_FE = FE_points[labels==i]
cluster_center_FE = FE_points[int(cluster_centers[i-1])]
diff_FE = cluster_FE-cluster_center_FE
new_labels = np.zeros(cluster_FE.shape[0])
new_labels[diff_FE< threshold] = i

labels[labels==i] = np.copy(new_labels)
i=5
threshold=0.25
cluster_FE = FE_points[labels==i]
cluster_center_FE = FE_points[int(cluster_centers[i-1])]
diff_FE = cluster_FE-cluster_center_FE
new_labels = np.zeros(cluster_FE.shape[0])
new_labels[diff_FE< threshold] = i

labels[labels==i] = np.copy(new_labels)
i=6
threshold=0.5
cluster_FE = FE_points[labels==i]
cluster_center_FE = FE_points[int(cluster_centers[i-1])]
diff_FE = cluster_FE-cluster_center_FE
new_labels = np.zeros(cluster_FE.shape[0])
new_labels[diff_FE< threshold] = i

labels[labels==i] = np.copy(new_labels)
fec.visualize(fontsize=12, savefig=True, xlabel=r'RMSD to closed crystal structure ($\AA$)',
                  ylabel=r'RMSD to open crystal structure ($\AA$)', zlabel='z', title=None,
                  vmax=4.5, colormap='Greys', lim=0.1, n_contour_levels=21, show_data=True, 
                  figsize= [5, 5], transition_point_size=0, core_point_size=0.25, 
                  core_point_edgecolor=None, 
                  colormap2=my_cmap,clust_point_edgecolor=None, clust_point_facecolor=None, 
                  clust_point_size=0, clust_point_linewidth=0, 
                  clust_point_marker='.',legend=None,
                  extend=None, rmsd_point_size=40, rmsd_point_fontsize=10, dx=1,rmsd_point_color='k',
                  core_point_alpha=1, transition_point_alpha=0.0, surface_alpha=1,
                  filename='./plots/2-Without_IPA/5-crisp_cluster_dt5_extended', 
                  ax=None, core_point_marker='o')
Counter(labels)<jupyter_output><empty_output><jupyter_text>Note that `label==0` corresponds to the non-assigned transition points! The rest are the crisp clusters (from 1 to 6) in the same order than previoulsy<jupyter_code>with open('./pickled/2-Without_IPA/labels_dt5_crisp.pickle', 'wb') as outf:
    pickle.dump(labels, outf)
with open('./pickled/2-Without_IPA/cluster_centers_dt5_crisp.pickle', 'wb') as outf:
    pickle.dump(cluster_centers, outf)
    
# with open('./pickled/2-Without_IPA/labels_dt5_crisp.pickle', 'rb') as inf:
#     labels=pickle.load(inf)
# with open('./pickled/2-Without_IPA/cluster_centers_dt5_crisp.pickle', 'rb') as inf:
#     cluster_centers=pickle.load(inf)
    
# fec.labels_ = np.copy(labels)
# fec.cluster_centers_ = np.copy(cluster_centers)<jupyter_output><empty_output><jupyter_text>Extracting representative structures Need to first load trajectory in memory as the following `t[labels==i]` doesn't work with `pt.iterload`.<jupyter_code>t = pt.load(traj_list, topfile, stride=5)
rf=[rf0,rf1,rf2,rf3,rf4,rf5]
toplist=sorted(glob('../1-Trajectories/0-Parm_rst/closed6_dry_rot*.parm7'))
toplist
sample_number=1000

for j in range(1,7):
    
    print("starting with", j)
    
    clust_traj = t[labels == j]
    print("Total number of frames:", clust_traj.n_frames)


    if clust_traj.n_frames > 0:
        
        randomlist = random.sample(range(0, clust_traj.n_frames), sample_number)
        tmp_traj=clust_traj[randomlist]
        print("Number of sampled frames:", tmp_traj.n_frames)
        
        if tmp_traj.n_frames > 0:

            new_traj = pt.Trajectory(top=topfile) # Empty trajectory
            min_ref=[]
    
            for i in range(0,tmp_traj.n_frames,1):

                rmsd0=pt.rmsd(tmp_traj[i:i+1],ref=rf0, mask="@CA,N,C,O")
                rmsd1=pt.rmsd(tmp_traj[i:i+1],ref=rf1, mask="@CA,N,C,O")
                rmsd2=pt.rmsd(tmp_traj[i:i+1],ref=rf2, mask="@CA,N,C,O")
                rmsd3=pt.rmsd(tmp_traj[i:i+1],ref=rf3, mask="@CA,N,C,O")
                rmsd4=pt.rmsd(tmp_traj[i:i+1],ref=rf4, mask="@CA,N,C,O")
                rmsd5=pt.rmsd(tmp_traj[i:i+1],ref=rf5, mask="@CA,N,C,O")

                rmsd_all=np.array([rmsd0,rmsd1,rmsd2,rmsd3,rmsd4,rmsd5])

                ali_struct=np.where(rmsd_all == np.nanmin(rmsd_all, axis=0))
                nb=int(ali_struct[0])
                new_frame=pt.align(tmp_traj[i:i+1],ref=rf[nb],mask="@CA,N,C,O", top=toplist[nb])
                new_traj.append(new_frame)
                min_ref.append(nb)

            new_traj.save('./sample_structures/2-Without_IPA/1-Crisp_1000/ms{}.xtc'.format(j), overwrite=True)  

            print("Fraction original trajectories", float(new_traj.n_frames/tmp_traj.n_frames))

            with open('./sample_structures/2-Without_IPA/1-Crisp_1000/ref_frame_ms{}.pickle'.format(j)
                      , 'wb') as outf:
                pickle.dump(min_ref, outf)
    print("finished with", j)<jupyter_output>starting with 1
Total number of frames: 8097
Number of sampled frames: 1000
Fraction original trajectories 1.0
finished with 1
starting with 2
Total number of frames: 37091
Number of sampled frames: 1000
Fraction original trajectories 1.0
finished with 2
starting with 3
Total number of frames: 3575
Number of sampled frames: 1000
Fraction original trajectories 1.0
finished with 3
starting with 4
Total number of frames: 1442
Number of sampled frames: 1000
Fraction original trajectories 1.0
finished with 4
starting with 5
Total number of frames: 2452
Number of sampled frames: 1000
Fraction original trajectories 1.0
finished with 5
starting with 6
Total number of frames: 6059
Number of sampled frames: 1000
Fraction original trajectories 1.0
finished with 6
{'fraction_non_alphanumeric': 0.13889588052688984, 'fraction_numerical': 0.03381156137259837, 'mean_word_length': 4.647713811892916, 'pattern_counts': {'":': 0, '<': 71, '<?xml version=': 0, '>': 67, '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': '17', 'dir': 'jupyter-structured-clean-dedup', 'id': '405760', 'license': 'MIT', 'n_tokens_mistral': 9693, 'n_tokens_neox': 9192, 'n_words': 1459, 'path': 'Analysis/2-RMSD_clustering_without_IPA_share.ipynb', 'repo_name': 'eric-jm-lang/MD-switch-paper', 'size': '23838'}