-
Notifications
You must be signed in to change notification settings - Fork 0
Expand file tree
/
Copy pathramachandran-plot-with-initial-values.py
More file actions
49 lines (41 loc) · 1.95 KB
/
Copy pathramachandran-plot-with-initial-values.py
File metadata and controls
49 lines (41 loc) · 1.95 KB
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
import matplotlib as mpl
import matplotlib.pyplot as plt
import numpy as np
import matplotlib.cm as cm
def combined_from_chi(rama, ramainit, peptidepath):
"""Plot ramachandran plot, combining the output of 'gmx chi'
input: xvg output of 'gmx chi' for a trajectory, xvg output of 'gmx chi' for the first frame, path for output
output: ramachandran plot with initial values
example:
peptidepath='/Users/User/testplots/' #directory where output plots will be stored
rmsfiles='rama_all.xvg' #output file from 'gmx chi'
ramainitial='first_frame.xvg' #output file from 'gmx chi'
combined_from_chi(rmsfiles,ramainitial, peptidepath)
"""
peptide=peptidepath+'_rama.svg' #output path and figure format
ramafile=np.loadtxt(rama) #load gmx chi output for trajectory
x=ramafile[:,0]
y=ramafile[:,1]
firstframe=np.loadtxt(ramainit) #load gmx chi output for initial frame
x1=firstframe[:,0]
y1=firstframe[:,1]
fig = plt.figure(figsize=(12, 10)) #define the size of the figure
plt.subplot(111)
# Next line define xrange, yrange. Notice that histogram2d will invert y, so we need to plot -y.
H1, xedges, yedges =np.histogram2d(-y,x, bins=90, range=[[-180, 180],[-180, 180]])
plt.imshow(H1,extent=[xedges[0], xedges[-1], yedges[0], yedges[-1]], cmap=cm.binary, interpolation='none')
cb = plt.colorbar()
cb.set_label('Counts')
#We are plotting the initial values on top of the histogram2d
plt.subplot(111)
#Defining the size of the points as 200
plt.scatter(x1, y1, 200)
plt.axis([-180, 180, -180, 180]) #xrange, yrange
plt.rcParams.update({'font.size' : 25}) # changes font size of labels
plt.rcParams.update({'axes.labelweight' : 'bold'}) #makes axes title bold
plt.rcParams.update({'font.weight' : 'bold'}) # make tick labels bold
plt.xticks(rotation='vertical')
plt.xlabel('Phi')
plt.ylabel('Psi')
plt.savefig(peptide)
plt.show()