|||
使用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里,不难做到。
Archiver|手机版|科学网 ( 京ICP备07017567号-12 )
GMT+8, 2024-9-12 11:15
Powered by ScienceNet.cn
Copyright © 2007- 中国科学报社