使用pyemma软件包：

http://www.emma-project.org/latest/index.html

按教程：

http://www.emma-project.org/latest/generated/pentapeptide_msm.html

Pentapeptide教程：

数据已经有了25组，各5000帧的数据。

放到data目录下。

indir = './data'

topfile =  indir+'/init-ww-penta.pdb'

feat = coor.featurizer(topfile)

feat.add_backbone_torsions(cossin=True)

feat.add_chi1_torsions(cossin=True)

Coordinates Package可以用来将模拟数据转换成离散的用于分析的数据。

比如add_distance可以取出每两个原子之间的距离。放到feat变量里。

chi1角，sidechain dihedral:

http://www.msg.ucsf.edu/local/programs/garlic/commands/dihedrals.html

这样feat就是一个高维的结构变量。

from glob import glob

traj_list = glob(indir+'/*-protein.xtc')

inp = coor.source(traj_list, feat)

这样inp就存储了按feat变量格式的轨迹文件。注意，”source” function defines input trajectories without loading them。

tica_obj = coor.tica(inp, lag=20, var_cutoff=0.9, kinetic_map=True)

Time-lagged independent component analysis基于时间的独立成分分析方法。这一步对数据进行了降维处理，应该需要一定的运算时间。

降维的数据存储在tica_obj变量中。

这时候可以将数据导入内存，变量Y.

Y = tica_obj.get_output() # get tica coordinates

mplt.plot_free_energy(np.vstack(Y)[:,0], np.vstack(Y)[:,1]);

可以对前两维进行自由能面的作图。

画图出错：

ValueError: Colormap spectral is not recognized

加上import pylab

plot_free_energy(cmap=pylab.cm.Blues)。

怎么得到自由能最低点对应的构型？

考虑这一点，所有帧的数据存储在Y里，不难做到。