Mercurial > repos > chemteam > mdanalysis_endtoend

comparison angle.py @ 0:78aa3659fcd1 draft

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
"planemo upload for repository https://github.com/galaxycomputationalchemistry/galaxy-tools-compchem/ commit 1b23e024af45cc0999d9142d07de6897d4189ec2"
author chemteam
date Mon, 24 Aug 2020 16:47:23 -0400
parents
children
comparison
equal deleted inserted replaced
-1:000000000000 0:78aa3659fcd1
1 #!/usr/bin/env python
2
3 import argparse
4 import csv
5 import sys
6
7 import MDAnalysis as mda
8
9 import matplotlib
10 import matplotlib.pyplot as plt
11
12 import numpy as np
13 from numpy.linalg import norm
14
15 matplotlib.use('Agg') # noqa
16
17
18 def parse_command_line(argv):
19 parser = argparse.ArgumentParser()
20 parser.add_argument('--itraj', help='input traj')
21 parser.add_argument('--istr', help='input str')
22 parser.add_argument('--itrajext', help='input traj ext')
23 parser.add_argument('--istrext', help='input str ext')
24 parser.add_argument('--isegid1', help='segid 1')
25 parser.add_argument('--iresid1', help='resid 1')
26 parser.add_argument('--iname1', help='name 1')
27 parser.add_argument('--isegid2', help='segid 2')
28 parser.add_argument('--iresid2', help='resid 2')
29 parser.add_argument('--iname2', help='name 2')
30 parser.add_argument('--isegid3', help='segid 3')
31 parser.add_argument('--iresid3', help='resid 3')
32 parser.add_argument('--iname3', help='name 3')
33 parser.add_argument('--output', help='output')
34 parser.add_argument('--oangle_plot', help='angle plot')
35 return parser.parse_args()
36
37
38 args = parse_command_line(sys.argv)
39
40 atom1 = "(segid %s and resid %s and name %s)" % \
41 (args.isegid1, args.iresid1, args.iname1)
42 atom2 = "(segid %s and resid %s and name %s)" % \
43 (args.isegid2, args.iresid2, args.iname2)
44 atom3 = "(segid %s and resid %s and name %s)" % \
45 (args.isegid3, args.iresid3, args.iname3)
46
47
48 def theta(u):
49 A = u.select_atoms(atom1).center_of_geometry()
50 B = u.select_atoms(atom2).center_of_geometry()
51 C = u.select_atoms(atom3).center_of_geometry()
52 BA = A - B
53 BC = C - B
54 theta = np.arccos(np.dot(BA, BC)/(norm(BA)*norm(BC)))
55 return np.rad2deg(theta)
56
57
58 u = mda.Universe(args.istr, args.itraj,
59 topology_format=args.istrext, format=args.itrajext)
60 data = np.array([(u.trajectory.frame, theta(u)) for ts in u.trajectory])
61 frame, theta = data.T
62
63 with open(args.output, 'w') as f:
64 writer = csv.writer(f, delimiter='\t')
65 writer.writerows(zip(frame, theta))
66
67 with open(args.output) as f:
68 g = [xtmp.strip() for xtmp in f]
69 data = [tuple(map(float, xtmp.split())) for xtmp in g[0:]]
70 time = [xtmp[0] for xtmp in data]
71 angle = [xtmp[1] for xtmp in data]
72 plt.plot(time, angle)
73 plt.xlabel('Frame No.')
74 plt.ylabel('Angle (degrees)')
75 plt.savefig(args.oangle_plot, format='png')