I'm trying to use the TPRParser to parse the TPR file generated from
GROMACs. Unfortunately it throws a error which I've never seen before:
Traceback (most recent call last):
line 43, in <module>
topo = mda.topology.TPRParser.TPRParser(TPR).parser(tpr_resid_from_one=True)
AttributeError: 'TPRParser' object has no attribute 'parser'
Could someone please help me to remove this error?
Here is my complete code and the data files are available here
import argparse
# define and parse command line arguments
parser = argparse.ArgumentParser()
parser.add_argument('input', metavar='INPUT', help='XTC input file')
parser.add_argument('--output', metavar='FILENAME', help='H5MD output file for saving trajectory data')
parser.add_argument('-v', '--verbose', action='store_true', help='be verbose')
args = parser.parse_args()
import MDAnalysis as mda
import numpy as np
import matplotlib.pyplot as plt
import os, sys
from collections import OrderedDict
XTC = args.input
TPR = os.path.splitext(XTC)[0] + '.tpr'
if not os.path.isfile(TPR):
print("Error: TPR file {0} does not exists".format(TPR))
sys.exit(1)
# the following list was generated manually using
# grep "^ *1[A-Z]" -B1 lp400start.gro | grep -v "^ *1[A-Z]" | cut -c 1-5
protein_seq_length = [ 38, 38, 137, 137, 252, 252, 576, 576, 1092, 1092, 1016, 1016, 1285, 1285 ]
protein_resid_range = []
i = 0
for n in protein_seq_length:
protein_resid_range += [ slice(i, i + n) ]
i += n
print("Expecting {0:d} residues in {1:d} proteins".format(i, len(protein_resid_range)))
# build reduced universe to match XTC
# (ignore warnings that no coordinates are found for the TPR)
topo = mda.topology.TPRParser.TPRParser(TPR).parser(tpr_resid_from_one=True)
u0 = mda.Universe(topo)
u = mda.Merge(u0.select_atoms("not resname W and not resname WF and not resname ION"))
u.load_new(XTC)
# this selection does not work, it mixes the atoms of both EPHA proteins
# FIXME delete this output
print(u.select_atoms("segid seg_0_EPHA"))
# select proteins through range of residue IDs
for idx in protein_resid_range:
res = u.residues[idx]
print("Protein of type {0:s} is composed of {1:d} residues and {2:n} atoms:\n{3}\n").format(res[0].segment.segid, len(res), len(res.atoms), res)
# segments
protein_segments = u.segments
# build the fragments
# (a dictionary with the key as the protein name -- I'm using an
# OrderedDict so that the order is the same as in the TPR)
protein_fragments = OrderedDict((seg.segid[1:], seg.atoms.fragments) for seq in protein_segments)
# this doesn't seem to be helpful either FIXME delete this output
for f in protein_fragments:
print(f)
# analyze trajectory
timeseries = []
for ts in u.trajectory[1:4]:
coms = []
for name, proteins in protein_fragments.items():
# loop over all the different proteins
# unwrap to get the true COM under PDB (double check!!)
coms.extend([p.center_of_mass(unwrap=True) for p in proteins])
timeseries.append(coms)
timeseries = np.array(timeseries, dtype=float)
if args.output:
with opne(args.output, 'w') as outfile:
# make clear distinction in frames
outfile.write('# Array shape: {0}\n'.format(timeseries.shape))
outfile.write('# New frame: {0}\n'.format(u.trajectory.time))
# Iterating through a n-dim array
for data in timeseries:
np.savetxt(outfile, data, fmt='%-7.2f')
You do not need to use the TPRParser explicitly. Instead of the lines
use just
The
Universe()class uses the correct parser (based on the filename extension) automatically. It is almost never necessary to use a parser explicitly. Just useUniverse().