-
Notifications
You must be signed in to change notification settings - Fork 1
/
anal_order.py
133 lines (115 loc) · 4.49 KB
/
anal_order.py
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
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
# Author: Samuel Genheden samuel.genheden@gmail.com
"""
Program to plot CG tail order parameters from CG and AA simulations
"""
import argparse
import os
import sys
import matplotlib.pylab as plt
import numpy as np
import openpyxl as xl
import scipy.stats as sstats
from sgenlib import colors
def _expand_selection(selection):
"""
Expand a simple pattern with {A..B} format to an integer range
"""
first = selection.index("{")
last = selection.index("}")
start,end = map(int,selection[first+1:last].split(".."))
return ["%s%d%s"%(selection[:first],i,selection[last+1:]) for i in range(start,end+1)]
def _plot_order(order_list,labels,figure):
"""
Plot the distribution of a series in subplots
"""
nchains = len(order_list[0])
ncols = 2.0
nrows = np.ceil(nchains / ncols)
axes = []
miny,maxy = 1.0,0.0
for ichain in range(nchains):
a = figure.add_subplot(nrows, ncols, ichain + 1)
axes.append(a)
for i,(simorder,label) in enumerate(zip(order_list,labels)):
order = simorder[ichain]
miny = min(miny,order.min())
maxy = max(maxy,order.max())
a.plot(np.arange(1,len(order)+1,1),order,'-*',label=label,color=colors.color(i))
a.set_xticks(np.arange(0,len(order)+2,1))
# Add legend and ticks
if ichain == 0:
a.legend(loc=1, fontsize=8)
for a in axes :
a.set_ylim([miny-0.05,maxy+0.05])
def _read_series(filename, chains_count):
"""
Read series of data output by md_order.py
"""
data = []
with open(filename, "r") as f:
line = f.readline()
while line:
data.append(line.strip().split())
line = f.readline()
data = np.array(data, dtype=float)
# Split the data for the different chains
chainorder = []
for i,cnt in enumerate(chains_count) :
first = 0 if i == 0 else chains_count[-1]-1
last = first + cnt
chainorder.append(data[:,first+1:last+1].mean(axis=0))
return chainorder
if __name__ == '__main__':
print " ".join(sys.argv)
# Command-line input
parser = argparse.ArgumentParser(description="Analyse CG order parameters")
parser.add_argument('-f', '--file', nargs="+", help="the trajectory of the order parameters")
parser.add_argument('-l', '--labels', nargs="+", help="the labels of the different trajectories")
parser.add_argument('-c', '--chains',nargs="+", help="the carbon chains")
parser.add_argument('-e','--excel',help="the filename of the XLSX file")
parser.add_argument('-s','--sheet',help="the sheet in the XLSX file")
parser.add_argument('-o', '--out', help="the output filename", default="order.png")
args = parser.parse_args()
# Parse the chain definition
chains_beads = []
chains_count = []
for chain in args.chains :
atoms = chain.split("-")
if len(atoms) == 1: atoms = _expand_selection(chain)
chains_beads.append(atoms)
chains_count.append(len(chains_beads[-1])-1)
# Read data, assume it was created by md_order.py
order_list = []
order_std_list = []
for filename in args.file:
chainorder = _read_series(filename, chains_count)
h,t = os.path.split(filename)
filename2 = os.path.join(h,"Repeat2",t)
chainorder2 = _read_series(filename2,chains_count)
chainorder_av = [np.asarray([d1,d2]).mean(axis=0)
for d1,d2 in zip(chainorder,chainorder2)]
chainorder_std = [np.asarray([d1,d2]).std(axis=0)/np.sqrt(2)
for d1,d2 in zip(chainorder,chainorder2)]
order_list.append(chainorder_av)
order_std_list.append(chainorder_std)
#
f = plt.figure(1)
_plot_order(order_list,args.labels,f)
f.savefig(args.out,format="png")
try :
wb = xl.load_workbook(filename = args.excel)
except :
print "Could not open the XLSX file. Will create one from scratch"
wb = xl.Workbook()
try :
ws = wb[args.sheet]
except :
ws = wb.create_sheet(title=args.sheet)
for i,(label,order,orderstd) in enumerate(zip(args.labels,order_list,order_std_list),2):
ws.cell(row=i,column=1).value = label
for j,(o,ostd) in enumerate(zip(order,orderstd)):
if i == 2:
ws.cell(row=1,column=j*2+2).value = "chain%d"%(j+1)
ws.cell(row=i,column=2+j*2).value = o.mean()
ws.cell(row=i,column=2+j*2+1).value = np.sqrt(np.sum(ostd*ostd))
wb.save(args.excel)