I am calculating Free Energy by binning 2 variables and getting a negative log of the histogram. I use Contacts and distance-between-helix from every frame of my trajectory as two variables and calculate 2D free energy. With the available free energy, I would like to get back to the frame which has minimum energy. But since I have already binned the data I am not able to trace back the contact-distance pair corresponding to 0KT energy before binning. Any help in approaching the problem is appreciated. I have attached the code and the corresponding Free Energy Plot herewith.
import sys
import numpy as np
import seaborn as sns
import matplotlib.pyplot as plt
# read input file
input_file = sys.argv[1]
data = np.genfromtxt(input_file, delimiter='\t')
# create 2d histogram
hist, xedges, yedges = np.histogram2d(data[:,0], data[:,1], bins=100)
x, y = np.meshgrid(xedges[:-1], yedges[:-1])
# calculate negative log of histogram population
hist_norm = (hist+0.1) / np.max(hist)
hist_log = -np.log10(hist_norm)
# plot histogram
sns.set_style('white')
fig, ax = plt.subplots(figsize=(12, 10))
im = ax.pcolormesh(x, y, hist_log.T, cmap='jet', shading='auto')
# set axis labels and title
ax.set_xlabel('IntraQ',fontsize=24, fontweight='bold')
ax.set_ylabel('Dist (nm)', fontsize=24, fontweight='bold')
ax.tick_params(axis='both', labelsize=20, width=2, length=7)
# set colorbar
cbar = fig.colorbar(im)
cbar.ax.tick_params(labelsize=20, width=2, length=7)
cbar.ax.set_ylabel('-log10(normalized population)', fontsize=24, fontweight='bold')
# set white background
ax.set_facecolor('white')
fig.patch.set_facecolor('white')
# save image to output file
output_file = sys.argv[2]
plt.savefig(output_file, dpi=300, bbox_inches='tight')
#write x, y, and hist_log.T to file
output_file_txt = output_file.split('.')[0] + '.txt'
with open(output_file_txt, 'w') as f:
for i in range(x.shape[0]):
for j in range(x.shape[1]):
f.write('{:.2f}\t{:.2f}\t{:.2f}\n'.format(x[i][j], y[i][j], hist_log.T[i][j]))
For Example, In the above figure, I know that one minimum lies in the range between ~300 Contact value and 1nm. How can I trace back this to the frame in the trajectory?
I tried writing the x, y and hist_log.T. But the Contact Distance value pair corresponding to 0.00 hist_log.T is not present in the input file (due to binning).
