Python MBAR.computePMF方法代码示例

# 需要导入模块: from pymbar import MBAR [as 别名]
# 或者: from pymbar.MBAR import computePMF [as 别名]
def test_mbar_computePMF():

    """ testing computePMF """

    name, test = generate_ho()
    x_n, u_kn, N_k_output, s_n = test.sample(N_k, mode='u_kn')
    mbar = MBAR(u_kn,N_k)
    #do a 1d PMF of the potential in the 3rd state:
    refstate = 2
    dx = 0.25
    xmin = test.O_k[refstate] - 1
    xmax = test.O_k[refstate] + 1
    within_bounds = (x_n >= xmin) & (x_n < xmax)
    bin_centers = dx*np.arange(np.int(xmin/dx),np.int(xmax/dx)) + dx/2
    bin_n = np.zeros(len(x_n),int)
    bin_n[within_bounds] = 1 + np.floor((x_n[within_bounds]-xmin)/dx)
    # 0 is reserved for samples outside the domain.  We will ignore this state
    range = np.max(bin_n)+1
    results = mbar.computePMF(u_kn[refstate,:], bin_n, range, uncertainties = 'from-specified', pmf_reference = 1)
    f_i = results['f_i']
    df_i = results['df_i']
    f0_i = 0.5*test.K_k[refstate]*(bin_centers-test.O_k[refstate])**2
    f_i, df_i = f_i[2:], df_i[2:] # first state is ignored, second is zero, with zero uncertainty
    normf0_i = f0_i[1:] - f0_i[0] # normalize to first state
    z = (f_i - normf0_i) / df_i
    eq(z / z_scale_factor, np.zeros(len(z)), decimal=0)

# 需要导入模块: from pymbar import MBAR [as 别名]
# 或者: from pymbar.MBAR import computePMF [as 别名]
pmf_analytical -= fzero
pmf_analytical[0] = 0

bin_n = numpy.zeros([numbrellas*nsamples], int)
# Determine indices of those within bounds.
within_bounds = (x_n[:,0] >= xmin) & (x_n[:,0] < xmax)
# Determine states for these.
bin_n[within_bounds] = 1 + numpy.floor((x_n[within_bounds,0]-xmin)/dx)
# Determine indices of bins that are not empty.
bin_counts = numpy.zeros([nbins], int)
for i in range(nbins):
  bin_counts[i] = (bin_n == i).sum()

# Compute PMF.
print("Computing PMF ...")
results = mbar.computePMF(u_n, bin_n, nbins, uncertainties = 'from-specified', pmf_reference = zeroindex)
f_i = results['f_i']
df_i = results['df_i']

# Show free energy and uncertainty of each occupied bin relative to lowest free energy

print("1D PMF:")
print("%d counts out of %d counts not in any bin" % (bin_counts[0],numbrellas*nsamples))
print("%8s %6s %8s %10s %10s %10s %10s %8s" % ('bin', 'x', 'N', 'f', 'true','error','df','sigmas'))
for i in range(1,nbins):
   if (i == zeroindex):
     stdevs = 0 
     df_i[0] = 0
   else:
     error = pmf_analytical[i]-f_i[i]
     stdevs = numpy.abs(error)/df_i[i]

# 需要导入模块: from pymbar import MBAR [as 别名]
# 或者: from pymbar.MBAR import computePMF [as 别名]
unbinned = 0
bin_kn = numpy.zeros([numbrellas, nsamples], int)
for i in range(numbrellas):
  # Determine indices of those within bounds.                                                                               
  within_bounds = (x_in[i,:,0] >= xmin) & (x_in[i,:,0] < xmax) & (x_in[i,:,1] >= ymin) & (x_in[i,:,1] < ymax)
  # Determine states for these.
  bin_kn[i,within_bounds] = 1 + numpy.floor((x_in[i,within_bounds,0]-xmin)/dx) + nbinsperdim*numpy.floor((x_in[i,within_bounds,1]-ymin)/dy)

# Determine indices of bins that are not empty.
bin_counts = numpy.zeros([nbins], int)
for i in range(nbins):                                                                                                      
  bin_counts[i] = (bin_kn == i).sum() 

# Compute PMF.                                
print "Computing PMF ..." 
[f_i, df_i] = mbar.computePMF(u_kn, bin_kn, nbins, uncertainties = 'from-specified', pmf_reference = zeroindex)  

# Show free energy and uncertainty of each occupied bin relative to lowest free energy                                        
print "2D PMF:"
print "%d counts out of %d counts not in any bin" % (bin_counts[0],numbrellas*nsamples)
print "%8s %6s %6s %8s %10s %10s %10s %10s %8s" % ('bin', 'x', 'y', 'N', 'f', 'true','error','df','sigmas')
for i in range(1,nbins):
   if (i == zeroindex):
     stdevs = 0 
     df_i[0] = 0
   else:  
     error = pmf_analytical[i]-f_i[i]
     stdevs = numpy.abs(error)/df_i[i]
   print '%8d %6.2f %6.2f %8d %10.3f %10.3f %10.3f %10.3f %8.2f' % (i, bin_centers[i,0], bin_centers[i,1] , bin_counts[i], f_i[i], pmf_analytical[i], error, df_i[i], stdevs)

"""

# 需要导入模块: from pymbar import MBAR [as 别名]
# 或者: from pymbar.MBAR import computePMF [as 别名]
f_i_pickle_fn = "f_i-of-"+"_".join(CV_fns)+str(stat_inefficiency)+"@"+str(N_bins)+".pickle"
try: #Try to read the pickle; 1) assert that the user wants this and 2) try to open the file
   assert (pickle_f_i_tables == True)
   pickle_file = open(f_i_pickle_fn, 'rb')
   print "(2 of 2) Loading f_i estimates from pickle file:", f_i_pickle_fn
   (f_i_tables, f_i_uncerts_tables) = load( pickle_file )
   pickle_file.close()

except (AssertionError, IOError):
   print "(2 of 2) Reweighting"
   f_is = []
   f_i_uncerts = []

   for i, target_beta in enumerate(beta_k[N_replicas:]):
      u_kn = target_beta * U_kn
      (f_i, d2f_ij) = mbar.computePMF(u_kn, bin_kn, bin_index)
      f_i = append(f_i, inf)
      d2f_ij = append(d2f_ij, 0.0)

      f_is.append(f_i)
      f_i_uncerts.append(d2f_ij)

   f_i_tables         = zeros([N_output_temps]+[N_bins]*N_CVs, dtype = float)
   f_i_uncerts_tables = zeros([N_output_temps]+[N_bins]*N_CVs, dtype = float)

   for i in range(N_output_temps):
      f_i_tables[i][:] = f_is[i][ bin_number[:] ]
      f_i_uncerts_tables[i][:] = f_i_uncerts[i][ bin_number[:] ]

   if pickle_f_i_tables == True:
      print "Pickling free energy estimates with file name:", f_i_pickle_fn