We will walk you through using the API of pyfeat for different simulation types and the different available estimators. As a simple initial example we will consider an asymmetric double well as the potential landscape and a single particle diffusing in this potential according to Brownian dynamics. For this purpose we have prepared a set of scripts that will run short simulations and write the trajectories to file in the correct pyfeat format. If you want to learn how to prepare data in a pyfeat format, it might be worthwhile to have a closer look at the data generating files (trajectory_factory.py).
In [ ]:
#imports
#allow for the embedding of plots into the ipython notebook.
%pylab inline
import trajectory_factory as tf #this is the package that allows the quick generation of input data for pyfeat
Let us run a simulated tempering simulation. CAREFUL the data generation may take a little while.
In [ ]:
tf.run_st_simulation() #generates simulated tempering data for pyfeat in the directory ST/
Have a look at the examples/ST/Traj.dat file. This will give you a good idea as to how a trajectory should be written out for your own data. All additional files written out will be needed for the analysis. The file kT.dat contains the reduced temperatures, b_k_i.dat is needed for WHAM and dTRAM and exact.dat contains the true stationary probabilities to which we will compare later.
Now we have generated the data we want for the pyfeat analysis, it is time to import all the necessary packages from pyfeat
In [ ]:
from pyfeat import xtram, wham, dtram #api function for pyfeat
from pyfeat import Reader #allows you to read the data from file in the correct format
from pyfeat import Forge #contains all the preformatted data that will then be passed to the estimators
The usual workflow for pyfeat is: Read the data with the reader, who will take a list of files as an agrument plus, any helper files such as a kT file or a b_K_i file. The reader object is then passed to the data converter (Forge).
In [ ]:
trajlist = ['ST/Traj.dat']
reader = Reader( trajlist, b_K_i_file = 'ST/b_K_i.dat', kT_file='ST/kT.dat' ) #read trajectory and 'helper files'
forge = Forge( reader.trajs, kT_K = reader.kT_K, b_K_i = reader.b_K_i, kT_target = 0 ) #pass read data to the data forge
In [ ]:
#load all the exact results
exact = np.loadtxt('ST/exact.dat')
exact[:,1] = exact[:,1]/np.sum(exact[:,1])
Now we have read all the data and can run our different estimators
In [ ]:
dtram_est = dtram( forge, lag=1 , maxiter=1000, ftol=1.0e-4, verbose=False )
print "#===============Thank you for using DTRAM============================="
dtram_est.cite(pre="# ")
print "#====================================================================="
In [ ]:
fig = plt.figure(1, figsize=(10,5))
fig.add_subplot(121)
plt.plot( exact[:,0],dtram_est.f_i , color = 'r', linewidth = 2, linestyle='--', label='DTRAM' )
plt.plot( exact[:,0], -np.log(exact[:,1]), color='k', label='exact' )
plt.xlabel( 'x in [a.u.]', fontsize = 20 )
plt.ylabel( 'F(x) in [kT]', fontsize = 20 )
plt.legend(loc=4)
fig.add_subplot(122)
plt.plot( exact[:,0],dtram_est.pi_i , color = 'r', linewidth = 2, linestyle='--', label='DTRAM' )
plt.plot( exact[:,0], exact[:,1], color='k', label='exact' )
plt.xlabel( 'x in [a.u.]', fontsize = 20 )
plt.ylabel( 'P(x)', fontsize = 20 )
plt.semilogy()
plt.legend(loc=1)
plt.tight_layout()
In [ ]:
wham_est = wham( forge, maxiter=1000, ftol=1.0e-7, verbose=False ) # and we are done, now we can analyse the results
print "#===============Thank you for using WHAM=================================================="
wham_est.cite(pre="# ")
print "#========================================================================================="
In [ ]:
fig = plt.figure(1, figsize=(10,5))
fig.add_subplot(121)
plt.plot( exact[:,0],wham_est.f_i , color = 'r', linewidth = 2, linestyle='--', label='WHAM' )
plt.plot( exact[:,0], -np.log(exact[:,1]), color='k', label='exact' )
plt.xlabel( 'x in [a.u.]', fontsize = 20 )
plt.ylabel( 'F(x) in [kT]', fontsize = 20 )
plt.legend(loc=4)
fig.add_subplot(122)
plt.plot( exact[:,0],wham_est.pi_i , color = 'r', linewidth = 2, linestyle='--', label='WHAM' )
plt.plot( exact[:,0], exact[:,1], color='k', label='exact' )
plt.xlabel( 'x in [a.u.]', fontsize = 20 )
plt.ylabel( 'P(x)', fontsize = 20 )
plt.semilogy()
plt.legend(loc=1)
plt.tight_layout()
In [ ]:
xtram_est = xtram( forge, lag=1 , maxiter=1000, ftol=1.0e-15, verbose=True )
print "#===============Thank you for using XTRAM===================================="
xtram_est.cite(pre="# ")
print "#============================================================================"
In [ ]:
fig = plt.figure(1, figsize=(10,5))
fig.add_subplot(121)
plt.plot( exact[:,0],xtram_est.f_i , color = 'r', linewidth = 2, linestyle='--', label='XTRAM' )
plt.plot( exact[:,0], -np.log(exact[:,1]), color='k', label='exact' )
plt.xlabel( 'x in [a.u.]', fontsize = 20 )
plt.ylabel( 'F(x) in [kT]', fontsize = 20 )
plt.legend(loc=4)
fig.add_subplot(122)
plt.plot( exact[:,0],xtram_est.pi_i , color = 'r', linewidth = 2, linestyle='--', label='XTRAM' )
plt.plot( exact[:,0], exact[:,1], color='k', label='exact' )
plt.xlabel( 'x in [a.u.]', fontsize = 20 )
plt.ylabel( 'P(x)', fontsize = 20 )
plt.semilogy()
plt.legend(loc=1)
plt.tight_layout()
In [ ]:
tf.run_us_simulation() #generates umbrella sampling data for pyfeat in the directory US/
In [ ]:
trajlist = ['US/Traj0.dat', 'US/Traj1.dat', 'US/Traj2.dat', 'US/Traj3.dat', 'US/Traj4.dat', 'US/Traj5.dat', 'US/Traj6.dat',
'US/Traj7.dat', 'US/Traj8.dat', 'US/Traj9.dat', 'US/Traj10.dat', 'US/Traj11.dat', 'US/Traj12.dat', 'US/Traj13.dat',
'US/Traj14.dat', 'US/Traj15.dat', 'US/Traj16.dat', 'US/Traj17.dat', 'US/Traj18.dat', 'US/Traj19.dat', 'US/Traj20.dat',
'US/Traj21.dat', 'US/Traj22.dat', 'US/Traj23.dat', 'US/Traj24.dat', 'US/Traj25.dat', 'US/Traj26.dat', 'US/Traj27.dat',
'US/Traj28.dat', 'US/Traj29.dat']
reader = Reader( trajlist, b_K_i_file = 'US/b_K_i.dat' ) #read trajectory and 'helper files'
forge = Forge( reader.trajs, b_K_i = reader.b_K_i ) #pass read data to the data forge
In [ ]:
#load all the exact results
exact = np.loadtxt('US/exact.dat')
exact[:,1] = exact[:,1]/np.sum(exact[:,1])
In [ ]:
dtram_est = dtram( forge, lag=1 , maxiter=10000, ftol=1.0e-6, verbose=False )
print "#===============Thank you for using DTRAM============================="
dtram_est.cite(pre="# ")
print "#====================================================================="
In [ ]:
fig = plt.figure(1, figsize=(10,5))
fig.add_subplot(121)
plt.plot( exact[:,0],dtram_est.f_i , color = 'r', linewidth = 2, linestyle='--', label='DTRAM' )
plt.plot( exact[:,0], -np.log(exact[:,1]), color='k', label='exact' )
plt.xlabel( 'x in [a.u.]', fontsize = 20 )
plt.ylabel( 'F(x) in [kT]', fontsize = 20 )
plt.legend(loc=4)
fig.add_subplot(122)
plt.plot( exact[:,0],dtram_est.pi_i , color = 'r', linewidth = 2, linestyle='--', label='DTRAM' )
plt.plot( exact[:,0], exact[:,1], color='k', label='exact' )
plt.xlabel( 'x in [a.u.]', fontsize = 20 )
plt.ylabel( 'P(x)', fontsize = 20 )
plt.semilogy()
plt.legend(loc=1)
plt.tight_layout()
In [ ]:
wham_est = wham( forge, maxiter=20000, ftol=1.0e-7, verbose=False ) # and we are done, now we can analyse the results
print "#===============Thank you for using WHAM=================================================="
wham_est.cite(pre="# ")
print "#========================================================================================="
In [ ]:
fig = plt.figure(1, figsize=(10,5))
fig.add_subplot(121)
plt.plot( exact[:,0],wham_est.f_i , color='r', linewidth=2, linestyle='--', label='WHAM' )
plt.plot( exact[:,0], -np.log(exact[:,1]), color='k', label='exact' )
plt.xlabel( 'x in [a.u.]', fontsize = 20 )
plt.ylabel( 'F(x) in [kT]', fontsize = 20 )
plt.legend(loc=4)
fig.add_subplot(122)
plt.plot( exact[:,0],wham_est.pi_i , color='r', linewidth=2, linestyle='--', label='WHAM' )
plt.plot( exact[:,0], exact[:,1], color='k', label='exact' )
plt.xlabel( 'x in [a.u.]', fontsize = 20 )
plt.ylabel( 'P(x)', fontsize = 20 )
plt.semilogy()
plt.legend(loc=1)
plt.tight_layout()
Thank you for running the pyfeat asymetric double-well example. If you have any further questions, or notice any bugs, please notify us on the mailinglist pyfeat@lists.fu-berlin.de.
In [ ]: