In [1]:
from msmbuilder import io
matplotlib.rcParams['font.size'] = 22
In [2]:
f = io.loadh('tICAData.h5')
f.keys()
Out[2]:
In [7]:
subplot(153)
hlines(f['vals'], 0, 1, lw=3)
xticks([])
xlim(0,1)
ylabel('eigenvalue')
Out[7]:
Notice that there are largely two timescales associated with alanine dipeptide, which largely agrees with most other analysis.
In [9]:
aa = np.loadtxt('HeavyAtomPairs.dat').astype(int)
In [146]:
v0 = f['vecs'][:,0]
v1 = f['vecs'][:,1]
#plot(v0, lw=3, label='tIC 0')
#plot(v1, lw=3, label='tIC 1')
p = Trajectory.load_from_pdb('native.pdb')
atms = np.unique(aa)
lbls = [ '$%s_{%d}$' % (p['AtomNames'][i], p['ResidueID'][i]) for i in atms ]
print atms
m0 = np.zeros((len(atms), len(atms)))
m1 = np.zeros((len(atms), len(atms)))
for i, (x, y) in enumerate(aa):
xi = np.where(atms == x )
yi = np.where(atms == y )
m0[xi,yi] = v0[i]
m1[xi,yi] = v1[i]
atmsi = np.arange(len(atms))
figure(figsize=(12,8))
ax = axes()
im=ax.matshow(m0 + m0.T, vmin=-100, vmax=100, cmap=cm.RdBu)
xticks(atmsi, lbls)
yticks(atmsi, lbls)
c=colorbar(im, )
#c.set_ticks([-100,-50, 0,50,100])
c.set_ticks([])
savefig('tic0.pdf')
figure(figsize=(8,6))
ax = axes()
im=ax.matshow(m1 + m1.T, vmin=-100, vmax=100, cmap=cm.RdBu)
xticks(atmsi, lbls)
yticks(atmsi, lbls)
c = colorbar(im)
#c.set_ticks([-100,-50, 0,50,100])
c.set_ticks([])
In [48]:
from msmbuilder import metrics
dih = metrics.Dihedral(angles='phi/psi')
atmpairs = metrics.AtomPairs(atom_pairs=aa)
tica = metrics.RedDimPNorm(proj_object_fn='tICAData.h5', prep_with=atmpairs, num_vecs=4)
In [ ]:
from msmbuilder import Trajectory
from msmbuilder import geometry
In [77]:
def plot_t(i):
t0 = Trajectory.load_from_lhdf('Trajectories/trj%d.lh5' %i)
phi0 = geometry.dihedral.compute_dihedrals(t0, indices=geometry.dihedral.get_indices(t0, angles='phi'))
psi0 = geometry.dihedral.compute_dihedrals(t0, indices=geometry.dihedral.get_indices(t0, angles='psi'))
t0 = tica.prepare_trajectory(t0)
subplot(121)
scatter(psi0, t0[:,0])
xlim([-180,180])
subplot(122)
scatter(phi0, t0[:,0])
xlim([-180,180])
In [80]:
figure(figsize=(10,5))
for i in range(100):
plot_t(i)
In [138]:
def scat_t(i, j, x=True, y=True):
t0 = Trajectory.load_from_lhdf('Trajectories/trj%d.lh5' %i)
phi0 = geometry.dihedral.compute_dihedrals(t0, indices=geometry.dihedral.get_indices(t0, angles='phi'))
psi0 = geometry.dihedral.compute_dihedrals(t0, indices=geometry.dihedral.get_indices(t0, angles='psi'))
t0 = tica.prepare_trajectory(t0)
#subplot(121)
#scatter(psi0, t0[:,0])
#xlim([-180,180])
#subplot(122)
#scatter(phi0, t0[:,0])
#xlim([-180,180])
scatter(phi0.flatten(), psi0.flatten(), c=t0[:,j], edgecolor='none', alpha=0.8)
xlim(-180,180)
ylim(-180,180)
# xi = [-180, -90, 0, 90, 180]
# xl = ['$-\pi$', '$-\pi / 2$', '$0$', '$\pi / 2$', '$\pi$']
xi = [-180, 0, 180]
xl = ['$-\pi$', '$0$', '$\pi$']
if x:
xlabel('$\phi$')
xticks(xi, xl)
else:
xticks([])
if y:
ylabel('$\psi$')
yticks(xi, xl)
else:
yticks([])
In [139]:
figure(figsize=(10,10))
subplot(334)
for i in range(100):
scat_t(i, 0)
title('tIC 0')
subplot(335)
for i in range(100):
scat_t(i, 1, y=False)
title('tIC 1')
subplot(336)
for i in range(100):
scat_t(i, 2, y=False)
title('tIC 2')
subplots_adjust(wspace=0.2)
savefig('/Users/schwancr/Installed/schwancrforks/msmbuilder/Docs/tICA/tics.png', dpi=200)
In [ ]: