#!/usr/bin/env python
import sys
from os import path, getcwd, popen
from os import chdir, system
import numpy as np
import matplotlib.pyplot as plt
import matplotlib as mpl

def reactioncoordinate( r, Z, idx ):
	if len(r) != len(Z):
		print('lengths of r and Z are not equal')
		exit()
	
	s = [ 0. ]
	for i in range( 1, len( r ) ):
		dZ = Z[i] - Z[i-1]
		dr = r[i] - r[i-1]
		s.append( s[-1] + np.sqrt( dZ**2 + dr**2 ) )
	ds = s[idx]
	for i in range( len(s) ):
		s[i] -= ds
	
	return s

#r [ Ang ]   Z [ Ang ]   Theta [ Deg ]   E [ kJ/mol ]
top = np.array(
[[1.2944, 2.9713, 101.71,  17.96],
[1.2949, 2.9603, 101.19,  18.39],
[1.2954, 2.9492, 100.70,  18.88],
[1.2958, 2.9381, 100.24,  19.41],
[1.2963, 2.9269, 99.80,  19.98],
[1.2967, 2.9155, 99.38,  20.57],
[1.2970, 2.9041, 98.92,  21.18],
[1.2974, 2.8925, 98.39,  21.81],
[1.2976, 2.8808, 97.81,  22.44],
[1.2979, 2.8688, 97.20,  23.08],
[1.2980, 2.8568, 96.62,  23.72],
[1.2982, 2.8446, 96.08,  24.37],
[1.2985, 2.8322, 95.60,  25.02],
[1.2990, 2.8197, 95.18,  25.69],
[1.2997, 2.8072, 94.79,  26.40],
[1.3006, 2.7946, 94.39,  27.18],
[1.3016, 2.7818, 94.00,  28.02],
[1.3025, 2.7688, 93.64,  28.92],
[1.3033, 2.7555, 93.33,  29.87],
[1.3042, 2.7420, 93.05,  30.86],
[1.3051, 2.7282, 92.77,  31.90],
[1.3061, 2.7143, 92.44,  33.00],
[1.3075, 2.7003, 91.99,  34.18],
[1.3091, 2.6862, 91.36,  35.44],
[1.3106, 2.6718, 90.63,  36.79],
[1.3120, 2.6569, 89.86,  38.22],
[1.3135, 2.6417, 89.08,  39.73],
[1.3152, 2.6264, 88.33,  41.32],
[1.3172, 2.6110, 87.61,  43.03],
[1.3199, 2.5957, 86.91,  44.88],
[1.3229, 2.5804, 86.23,  46.88],
[1.3259, 2.5647, 85.63,  49.02],
[1.3291, 2.5487, 85.09,  51.32],
[1.3327, 2.5327, 84.56,  53.77],
[1.3369, 2.5168, 83.98,  56.40],
[1.3419, 2.5014, 83.31,  59.24],
[1.3476, 2.4865, 82.61,  62.28],
[1.3540, 2.4719, 81.96,  65.54],
[1.3614, 2.4581, 81.46,  69.02],
[1.3700, 2.4454, 81.23,  72.73],
[1.3804, 2.4349, 81.50,  76.66],
[1.3939, 2.4284, 82.92,  80.79],
[1.4134, 2.4303, 87.38,  85.04],
[1.4389, 2.4415, 96.52,  89.18],
[1.4662, 2.4568, 107.57,  92.96],
[1.5744, 2.5986, 131.15,  95.26],
[1.5919, 2.6070, 132.18,  96.76],
[1.6062, 2.6110, 132.76,  98.16],
[1.6198, 2.6147, 133.26,  99.49],
[1.6330, 2.6184, 133.75, 100.76],
[1.6460, 2.6222, 134.23, 101.97],
[1.6587, 2.6264, 134.72, 103.11],
[1.6712, 2.6309, 135.21, 104.19],
[1.6836, 2.6359, 135.68, 105.20],
[1.6956, 2.6409, 136.10, 106.14],
[1.7072, 2.6457, 136.47, 107.01],
[1.7183, 2.6502, 136.79, 107.82],
[1.7289, 2.6543, 137.04, 108.57],
[1.7391, 2.6578, 137.24, 109.27],
[1.7488, 2.6607, 137.40, 109.91],
[1.7580, 2.6628, 137.50, 110.51],
[1.7668, 2.6643, 137.57, 111.06],
[1.7753, 2.6650, 137.61, 111.58],
[1.7834, 2.6649, 137.62, 112.06],
[1.7913, 2.6641, 137.60, 112.51],
[1.7990, 2.6624, 137.56, 112.93],
[1.8065, 2.6601, 137.50, 113.32],
[1.8139, 2.6570, 137.41, 113.68],
[1.8212, 2.6532, 137.30, 114.01],
[1.8285, 2.6487, 137.18, 114.31],
[1.8359, 2.6434, 137.03, 114.57],
[1.8432, 2.6375, 136.86, 114.81],
[1.8507, 2.6309, 136.67, 115.02],
[1.8584, 2.6237, 136.46, 115.19],
[1.8662, 2.6160, 136.23, 115.32],
[1.8742, 2.6077, 135.99, 115.42],
[1.8825, 2.5991, 135.73, 115.48],
[1.8911, 2.5901, 135.46, 115.49],
[1.9092, 2.5717, 134.91, 115.40],
[1.9188, 2.5625, 134.64, 115.29],
[1.9287, 2.5535, 134.37, 115.13],
[1.9388, 2.5448, 134.12, 114.93],
[1.9493, 2.5366, 133.88, 114.69],
[1.9600, 2.5289, 133.67, 114.42],
[1.9710, 2.5219, 133.47, 114.11],
[1.9822, 2.5155, 133.29, 113.76],
[1.9935, 2.5100, 133.14, 113.40],
[2.0049, 2.5051, 132.99, 113.01],
[2.0164, 2.5013, 132.86, 112.61],
[2.0279, 2.4985, 132.74, 112.20],
[2.0393, 2.4966, 132.62, 111.80],
[2.0508, 2.4955, 132.49, 111.41],
[2.0622, 2.4949, 132.35, 111.03],
[2.0736, 2.4947, 132.16, 110.67],
[2.0850, 2.4949, 131.93, 110.35],
[2.0965, 2.4955, 131.64, 110.06],
[2.1080, 2.4959, 131.25, 109.80],
[2.1200, 2.4957, 130.72, 109.59],
[2.1324, 2.4946, 130.04, 109.41],
[2.1458, 2.4917, 129.13, 109.24],
[2.1612, 2.4847, 127.78, 109.08],
[2.1811, 2.4686, 125.53, 108.89],
[2.2016, 2.4532, 123.36, 108.63],
[2.2196, 2.4448, 121.89, 108.32],
[2.2359, 2.4409, 120.84, 107.98],
[2.2510, 2.4400, 120.07, 107.64],
[2.2652, 2.4411, 119.48, 107.31]]
)

TS = np.array(
[[1.2987, 2.9685, 113.74,  18.80],
[1.2994, 2.9370, 112.57,  20.15],
[1.3001, 2.9054, 111.36,  21.75],
[1.3007, 2.8736, 110.10,  23.42],
[1.3012, 2.8415, 108.83,  25.04],
[1.3020, 2.8091, 107.67,  26.67],
[1.3037, 2.7765, 106.85,  28.56],
[1.3055, 2.7435, 106.30,  30.69],
[1.3075, 2.7100, 105.85,  33.02],
[1.3104, 2.6762, 105.45,  35.67],
[1.3133, 2.6417, 104.98,  38.61],
[1.3168, 2.6068, 104.56,  41.86],
[1.3222, 2.5718, 104.45,  45.59],
[1.3283, 2.5365, 104.50,  49.79],
[1.3361, 2.5012, 104.78,  54.54],
[1.3471, 2.4671, 105.45,  59.96],
[1.3617, 2.4345, 106.60,  66.03],
[1.3838, 2.4069, 108.83,  72.70],
[1.4295, 2.3980, 113.54,  79.51],
[1.4927, 2.4065, 118.18,  85.49],
[1.5581, 2.4220, 121.06,  90.28],
[1.6335, 2.4529, 123.90,  93.69],
[1.6923, 2.4743, 125.32,  95.88],
[1.7373, 2.4861, 125.93,  97.35],
[1.7757, 2.4944, 126.20,  98.37],
[1.8084, 2.4987, 126.04,  99.09],
[1.8375, 2.5010, 125.54,  99.60],
[1.8641, 2.5019, 124.86,  99.97],
[1.8883, 2.5012, 124.12, 100.24],
[1.9110, 2.4993, 123.36, 100.43],
[1.9325, 2.4969, 122.64, 100.57],
[1.9533, 2.4942, 121.95, 100.65],
[1.9734, 2.4910, 121.28, 100.70],
[1.9929, 2.4873, 120.63, 100.71],
[2.0119, 2.4833, 119.99, 100.69],
[2.0307, 2.4788, 119.37, 100.64],
[2.0492, 2.4740, 118.76, 100.56],
[2.0676, 2.4689, 118.15, 100.46],
[2.0860, 2.4636, 117.54, 100.33],
[2.1044, 2.4579, 116.93, 100.18],
[2.1230, 2.4521, 116.32, 100.01],
[2.1418, 2.4461, 115.71,  99.81],
[2.1609, 2.4403, 115.11,  99.59],
[2.1803, 2.4350, 114.54,  99.34]]
)

fcc = np.array(
[[1.3075, 2.9688, 113.40,  24.41],
[1.3095, 2.9377, 113.90,  26.15],
[1.3115, 2.9066, 114.67,  28.11],
[1.3138, 2.8754, 115.60,  30.16],
[1.3159, 2.8440, 116.22,  32.22],
[1.3177, 2.8124, 116.36,  34.30],
[1.3207, 2.7807, 116.62,  36.58],
[1.3246, 2.7489, 117.18,  39.12],
[1.3291, 2.7169, 117.92,  41.86],
[1.3351, 2.6850, 118.71,  44.82],
[1.3422, 2.6532, 119.36,  47.99],
[1.3507, 2.6217, 120.11,  51.43],
[1.3633, 2.5917, 121.14,  55.21],
[1.3792, 2.5634, 121.35,  59.24],
[1.3949, 2.5349, 119.80,  63.39],
[1.4092, 2.5056, 117.12,  67.68],
[1.4245, 2.4766, 114.93,  72.27],
[1.4415, 2.4486, 113.57,  77.19],
[1.4623, 2.4233, 113.60,  82.47],
[1.4925, 2.4060, 115.52,  88.00],
[1.6892, 2.4322, 117.04, 106.77],
[1.7413, 2.4529, 117.52, 110.17],
[1.7890, 2.4734, 117.78, 112.99],
[1.8316, 2.4924, 117.94, 115.36],
[1.8698, 2.5099, 118.09, 117.35],
[1.9039, 2.5255, 118.27, 119.04],
[1.9342, 2.5389, 118.44, 120.49],
[1.9613, 2.5503, 118.55, 121.76],
[1.9858, 2.5599, 118.61, 122.87],
[2.0081, 2.5679, 118.63, 123.86],
[2.0285, 2.5743, 118.62, 124.75],
[2.0471, 2.5788, 118.59, 125.55],
[2.0644, 2.5813, 118.53, 126.28],
[2.0804, 2.5815, 118.41, 126.95],
[2.0954, 2.5788, 118.23, 127.56],
[2.1096, 2.5726, 117.95, 128.12],
[2.1232, 2.5621, 117.56, 128.63],
[2.1366, 2.5461, 117.02, 129.07],
[2.1504, 2.5242, 116.32, 129.45],
[2.1653, 2.4969, 115.52, 129.75],
[2.1818, 2.4663, 114.67, 129.93]]
)

bridge = np.array(
[[1.2991, 2.9685, 116.86,  20.31],
[1.2996, 2.9370, 114.99,  21.78],
[1.3000, 2.9054, 113.45,  23.48],
[1.3005, 2.8736, 112.20,  25.25],
[1.3008, 2.8414, 111.07,  26.96],
[1.3015, 2.8090, 109.81,  28.69],
[1.3030, 2.7763, 108.18,  30.68],
[1.3042, 2.7431, 106.28,  32.92],
[1.3056, 2.7094, 104.39,  35.32],
[1.3077, 2.6752, 102.90,  38.00],
[1.3095, 2.6402, 101.70,  40.91],
[1.3118, 2.6045, 100.72,  44.04],
[1.3157, 2.5687, 99.99,  47.59],
[1.3194, 2.5318, 99.19,  51.48],
[1.3247, 2.4946, 98.44,  55.81],
[1.3319, 2.4575, 97.95,  60.69],
[1.3399, 2.4199, 97.88,  66.04],
[1.3515, 2.3835, 99.33,  71.90],
[1.3699, 2.3515, 103.08,  78.23],
[1.4077, 2.3351, 110.70,  84.82],
[1.4736, 2.3460, 116.96,  90.61],
[1.5250, 2.3482, 119.09,  95.32],
[1.5778, 2.3557, 121.54,  99.22],
[1.6360, 2.3735, 123.58, 102.16],
[1.6861, 2.3875, 124.35, 104.21],
[1.7275, 2.3952, 124.42, 105.68],
[1.7642, 2.4002, 124.19, 106.75],
[1.7973, 2.4029, 123.76, 107.54],
[1.8274, 2.4037, 123.17, 108.12],
[1.8552, 2.4028, 122.46, 108.56],
[1.8814, 2.4007, 121.66, 108.87],
[1.9061, 2.3973, 120.81, 109.09],
[1.9297, 2.3928, 119.95, 109.23],
[1.9525, 2.3874, 119.08, 109.30],
[1.9748, 2.3811, 118.22, 109.31],
[1.9966, 2.3739, 117.37, 109.24],
[2.0182, 2.3659, 116.54, 109.12],
[2.0397, 2.3570, 115.70, 108.93],
[2.0613, 2.3472, 114.86, 108.68],
[2.0830, 2.3366, 114.00, 108.37],
[2.1052, 2.3253, 113.14, 108.00],
[2.1279, 2.3135, 112.29, 107.57],
[2.1512, 2.3013, 111.47, 107.07],
[2.1752, 2.2889, 110.65, 106.50]]
)

def plot_traj(name, r_TS):
	data = open(name, 'r').readlines()
	
	r = []
	Z = []
	theta = []
	for i in range(1, len(data)):
		r.append( float( data[i].split()[1] ) )
		Z.append( float( data[i].split()[2] ) )
		theta.append( float( data[i].split()[8] ) )
		
		if r[-1] > r_TS:
			try:
				idx
			except NameError:
				idx = i-1
	
	xaxis = reactioncoordinate( r, Z, idx )
	plt.plot( xaxis, theta )
	
prop_cycle = plt.rcParams['axes.prop_cycle']
color_cycle = prop_cycle.by_key()['color']

xlim = [-7.,0.5]
yticks = [0., 45., 90., 135.]

Nrows=3
Ncols=2
fig, ax = plt.subplots(nrows=Nrows, ncols=Ncols, figsize=(6.69,5.))
#mpl.rc("font", size=10)
mpl.rc("text", usetex=True)

plt.subplot(Nrows,Ncols,1)

plot_traj('2.56ms_v0j2_angles/analysis_rcom_000430.dat', 1.89)
plot_traj('2.56ms_v0j2_angles/analysis_rcom_000633.dat', 1.89)
plot_traj('2.56ms_v0j2_angles/analysis_rcom_002711.dat', 1.89)
plot_traj('2.56ms_v0j2_angles/analysis_rcom_003321.dat', 1.89)
plot_traj('2.56ms_v0j2_angles/analysis_rcom_005141.dat', 1.89)
plot_traj('2.56ms_v0j2_angles/analysis_rcom_006941.dat', 1.89)
plot_traj('2.56ms_v0j2_angles/analysis_rcom_007160.dat', 1.89)
plot_traj('2.56ms_v0j2_angles/analysis_rcom_007483.dat', 1.89)

xaxis = reactioncoordinate( TS[:,0], TS[:,1], np.argmax(TS[:,3]) )
#plt.plot( xaxis, TS[:,2], marker='o', label='TS' )

xaxis = reactioncoordinate( top[:,0], top[:,1], np.argmax(top[:,3]) )
plt.scatter( xaxis, top[:,2], marker='o', label='Top', facecolor='None', color=color_cycle[1], zorder=100, alpha=0.5 )

xaxis = reactioncoordinate( fcc[:,0], fcc[:,1], np.argmax(fcc[:,3]) )
#plt.scatter( xaxis, fcc[:,2], marker='o', label='Hollow', facecolor='None', color=color_cycle[2], zorder=100, alpha=0.5 )

xaxis = reactioncoordinate( bridge[:,0], bridge[:,1], np.argmax(bridge[:,3]) )
#plt.scatter( xaxis, bridge[:,2], marker='o', label='Bridge', facecolor='None', color=color_cycle[3], zorder=100, alpha=0.5 )

plt.plot( [0.,0.], [0.,180.], color='k', ls='--' )

plt.title(r'$J=2$')

#plt.annotate('(e)', xy=(0.1,0.85), xycoords='axes fraction')

plt.legend(loc='best', numpoints=1, frameon=True, handletextpad=0.)
#plt.legend(loc='upper left', numpoints=1, handletextpad=0.5, borderaxespad=0.2, frameon=False)
plt.xlim(xlim)
plt.ylim(0.,180.)
plt.yticks(yticks)
plt.tick_params(length=6, width=1, direction='in', top=True, right=True, labelbottom=False)
plt.ylabel(r'$\theta$ (degrees)')
#plt.xlabel(r'Reaction path (\r{A})')

plt.subplot(Nrows,Ncols,3)

plot_traj('2.56ms_v0j2_angles/analysis_rcom_000052.dat', 2.20)
plot_traj('2.56ms_v0j2_angles/analysis_rcom_002675.dat', 2.20)
plot_traj('2.56ms_v0j2_angles/analysis_rcom_006035.dat', 2.20)
plot_traj('2.56ms_v0j2_angles/analysis_rcom_006958.dat', 2.20)
plot_traj('2.56ms_v0j2_angles/analysis_rcom_008720.dat', 2.20)
plot_traj('2.56ms_v0j2_angles/analysis_rcom_008761.dat', 2.20)
plot_traj('2.56ms_v0j2_angles/analysis_rcom_009000.dat', 2.20)
plot_traj('2.56ms_v0j2_angles/analysis_rcom_009146.dat', 2.20)

xaxis = reactioncoordinate( TS[:,0], TS[:,1], np.argmax(TS[:,3]) )
#plt.plot( xaxis, TS[:,2], marker='o', label='TS' )

xaxis = reactioncoordinate( top[:,0], top[:,1], np.argmax(top[:,3]) )
#plt.scatter( xaxis, top[:,2], marker='o', label='Top', facecolor='None', color=color_cycle[1], zorder=100, alpha=0.5 )

xaxis = reactioncoordinate( fcc[:,0], fcc[:,1], np.argmax(fcc[:,3]) )
plt.scatter( xaxis, fcc[:,2], marker='o', label='Hollow', facecolor='None', color=color_cycle[2], zorder=100, alpha=0.5 )

xaxis = reactioncoordinate( bridge[:,0], bridge[:,1], np.argmax(bridge[:,3]) )
#plt.scatter( xaxis, bridge[:,2], marker='o', label='Bridge', facecolor='None', color=color_cycle[3], zorder=100, alpha=0.5 )

plt.plot( [0.,0.], [0.,180.], color='k', ls='--' )

#plt.annotate('(e)', xy=(0.1,0.85), xycoords='axes fraction')

plt.legend(loc='best', numpoints=1, frameon=True, handletextpad=0.)
#plt.legend(loc='upper left', numpoints=1, handletextpad=0.5, borderaxespad=0.2, frameon=False)
plt.xlim(xlim)
plt.ylim(0.,180.)
plt.yticks(yticks)
plt.tick_params(length=6, width=1, direction='in', top=True, right=True, labelbottom=False)
plt.ylabel(r'$\theta$ (degrees)')
#plt.xlabel(r'Reaction path (\r{A})')

plt.subplot(Nrows,Ncols,5)

plot_traj('2.56ms_v0j2_angles/analysis_rcom_000162.dat', 1.96)
plot_traj('2.56ms_v0j2_angles/analysis_rcom_000774.dat', 1.96)
plot_traj('2.56ms_v0j2_angles/analysis_rcom_002626.dat', 1.96)
plot_traj('2.56ms_v0j2_angles/analysis_rcom_003602.dat', 1.96)
plot_traj('2.56ms_v0j2_angles/analysis_rcom_004409.dat', 1.96)
plot_traj('2.56ms_v0j2_angles/analysis_rcom_005570.dat', 1.96)
plot_traj('2.56ms_v0j2_angles/analysis_rcom_007411.dat', 1.96)
plot_traj('2.56ms_v0j2_angles/analysis_rcom_009798.dat', 1.96)

xaxis = reactioncoordinate( TS[:,0], TS[:,1], np.argmax(TS[:,3]) )
#plt.plot( xaxis, TS[:,2], marker='o', label='TS' )

xaxis = reactioncoordinate( top[:,0], top[:,1], np.argmax(top[:,3]) )
#plt.scatter( xaxis, top[:,2], marker='o', label='Top', facecolor='None', color=color_cycle[1], zorder=100, alpha=0.5 )

xaxis = reactioncoordinate( fcc[:,0], fcc[:,1], np.argmax(fcc[:,3]) )
#plt.scatter( xaxis, fcc[:,2], marker='o', label='Hollow', facecolor='None', color=color_cycle[2], zorder=100, alpha=0.5 )

xaxis = reactioncoordinate( bridge[:,0], bridge[:,1], np.argmax(bridge[:,3]) )
plt.scatter( xaxis, bridge[:,2], marker='o', label='Bridge', facecolor='None', color=color_cycle[3], zorder=100, alpha=0.5 )

plt.plot( [0.,0.], [0.,180.], color='k', ls='--' )

#plt.annotate('(e)', xy=(0.1,0.85), xycoords='axes fraction')

plt.legend(loc='best', numpoints=1, frameon=True, handletextpad=0.)
#plt.legend(loc='upper left', numpoints=1, handletextpad=0.5, borderaxespad=0.2, frameon=False)
plt.xlim(xlim)
plt.ylim(0.,180.)
plt.yticks(yticks)
plt.tick_params(length=6, width=1, direction='in', top=True, right=True)
plt.ylabel(r'$\theta$ (degrees)')
plt.xlabel(r'Reaction coordinate (\r{A})')

plt.subplot(Nrows,Ncols,2)

plot_traj('2.56ms_v0j8_angles/analysis_rcom_000117.dat', 1.89)
plot_traj('2.56ms_v0j8_angles/analysis_rcom_000929.dat', 1.89)
plot_traj('2.56ms_v0j8_angles/analysis_rcom_001309.dat', 1.89)
plot_traj('2.56ms_v0j8_angles/analysis_rcom_002087.dat', 1.89)
plot_traj('2.56ms_v0j8_angles/analysis_rcom_002385.dat', 1.89)
plot_traj('2.56ms_v0j8_angles/analysis_rcom_002926.dat', 1.89)
plot_traj('2.56ms_v0j8_angles/analysis_rcom_004121.dat', 1.89)
plot_traj('2.56ms_v0j8_angles/analysis_rcom_008935.dat', 1.89)
plot_traj('2.56ms_v0j8_angles/analysis_rcom_009477.dat', 1.89)
plot_traj('2.56ms_v0j8_angles/analysis_rcom_009652.dat', 1.89)

xaxis = reactioncoordinate( TS[:,0], TS[:,1], np.argmax(TS[:,3]) )
#plt.plot( xaxis, TS[:,2], marker='o', label='TS' )

xaxis = reactioncoordinate( top[:,0], top[:,1], np.argmax(top[:,3]) )
plt.scatter( xaxis, top[:,2], marker='o', label='Top', facecolor='None', color=color_cycle[1], zorder=100, alpha=0.5 )

xaxis = reactioncoordinate( fcc[:,0], fcc[:,1], np.argmax(fcc[:,3]) )
#plt.scatter( xaxis, fcc[:,2], marker='o', label='Hollow', facecolor='None', color=color_cycle[2], zorder=100, alpha=0.5 )

xaxis = reactioncoordinate( bridge[:,0], bridge[:,1], np.argmax(bridge[:,3]) )
#plt.scatter( xaxis, bridge[:,2], marker='o', label='Bridge', facecolor='None', color=color_cycle[3], zorder=100, alpha=0.5 )

plt.plot( [0.,0.], [0.,180.], color='k', ls='--' )

plt.title(r'$J=8$')

#plt.annotate('(e)', xy=(0.1,0.85), xycoords='axes fraction')

plt.legend(loc='best', numpoints=1, frameon=True, handletextpad=0.)
#plt.legend(loc='upper left', numpoints=1, handletextpad=0.5, borderaxespad=0.2, frameon=False)
plt.xlim(xlim)
plt.ylim(0.,180.)
plt.yticks(yticks)
plt.tick_params(length=6, width=1, direction='in', top=True, right=True, labelbottom=False, labelleft=False)
#plt.ylabel(r'$\theta$ (degrees)')
#plt.xlabel(r'Reaction path (\r{A})')

plt.subplot(Nrows,Ncols,4)

plot_traj('2.56ms_v0j8_angles/analysis_rcom_000052.dat', 2.20)
plot_traj('2.56ms_v0j8_angles/analysis_rcom_000220.dat', 2.20)
plot_traj('2.56ms_v0j8_angles/analysis_rcom_000224.dat', 2.20)
plot_traj('2.56ms_v0j8_angles/analysis_rcom_000277.dat', 2.20)
plot_traj('2.56ms_v0j8_angles/analysis_rcom_001103.dat', 2.20)
plot_traj('2.56ms_v0j8_angles/analysis_rcom_001511.dat', 2.20)
plot_traj('2.56ms_v0j8_angles/analysis_rcom_002477.dat', 2.20)
plot_traj('2.56ms_v0j8_angles/analysis_rcom_002979.dat', 2.20)

xaxis = reactioncoordinate( TS[:,0], TS[:,1], np.argmax(TS[:,3]) )
#plt.plot( xaxis, TS[:,2], marker='o', label='TS' )

xaxis = reactioncoordinate( top[:,0], top[:,1], np.argmax(top[:,3]) )
#plt.scatter( xaxis, top[:,2], marker='o', label='Top', facecolor='None', color=color_cycle[1], zorder=100, alpha=0.5 )

xaxis = reactioncoordinate( fcc[:,0], fcc[:,1], np.argmax(fcc[:,3]) )
plt.scatter( xaxis, fcc[:,2], marker='o', label='Hollow', facecolor='None', color=color_cycle[2], zorder=100, alpha=0.5 )

xaxis = reactioncoordinate( bridge[:,0], bridge[:,1], np.argmax(bridge[:,3]) )
#plt.scatter( xaxis, bridge[:,2], marker='o', label='Bridge', facecolor='None', color=color_cycle[3], zorder=100, alpha=0.5 )

plt.plot( [0.,0.], [0.,180.], color='k', ls='--' )

#plt.annotate('(e)', xy=(0.1,0.85), xycoords='axes fraction')

plt.legend(loc='best', numpoints=1, frameon=True, handletextpad=0.)
#plt.legend(loc='upper left', numpoints=1, handletextpad=0.5, borderaxespad=0.2, frameon=False)
plt.xlim(xlim)
plt.ylim(0.,180.)
plt.yticks(yticks)
plt.tick_params(length=6, width=1, direction='in', top=True, right=True, labelbottom=False, labelleft=False)
#plt.ylabel(r'$\theta$ (degrees)')
#plt.xlabel(r'Reaction path (\r{A})')

plt.subplot(Nrows,Ncols,6)

plot_traj('2.56ms_v0j8_angles/analysis_rcom_000083.dat', 1.96)
plot_traj('2.56ms_v0j8_angles/analysis_rcom_000666.dat', 1.96)
plot_traj('2.56ms_v0j8_angles/analysis_rcom_002034.dat', 1.96)
plot_traj('2.56ms_v0j8_angles/analysis_rcom_004064.dat', 1.96)
plot_traj('2.56ms_v0j8_angles/analysis_rcom_005695.dat', 1.96)
plot_traj('2.56ms_v0j8_angles/analysis_rcom_007004.dat', 1.96)
plot_traj('2.56ms_v0j8_angles/analysis_rcom_008606.dat', 1.96)
plot_traj('2.56ms_v0j8_angles/analysis_rcom_009889.dat', 1.96)

xaxis = reactioncoordinate( TS[:,0], TS[:,1], np.argmax(TS[:,3]) )
#plt.plot( xaxis, TS[:,2], marker='o', label='TS' )

xaxis = reactioncoordinate( top[:,0], top[:,1], np.argmax(top[:,3]) )
#plt.scatter( xaxis, top[:,2], marker='o', label='Top', facecolor='None', color=color_cycle[1], zorder=100, alpha=0.5 )

xaxis = reactioncoordinate( fcc[:,0], fcc[:,1], np.argmax(fcc[:,3]) )
#plt.scatter( xaxis, fcc[:,2], marker='o', label='Hollow', facecolor='None', color=color_cycle[2], zorder=100, alpha=0.5 )

xaxis = reactioncoordinate( bridge[:,0], bridge[:,1], np.argmax(bridge[:,3]) )
plt.scatter( xaxis, bridge[:,2], marker='o', label='Bridge', facecolor='None', color=color_cycle[3], zorder=100, alpha=0.5 )

plt.plot( [0.,0.], [0.,180.], color='k', ls='--' )

#plt.annotate('(e)', xy=(0.1,0.85), xycoords='axes fraction')

plt.legend(loc='best', numpoints=1, frameon=True, handletextpad=0.)
#plt.legend(loc='upper left', numpoints=1, handletextpad=0.5, borderaxespad=0.2, frameon=False)
plt.xlim(xlim)
plt.ylim(0.,180.)
plt.yticks(yticks)
plt.tick_params(length=6, width=1, direction='in', top=True, right=True, labelleft=False)
#plt.ylabel(r'$\theta$ (degrees)')
plt.xlabel(r'Reaction coordinate (\r{A})')

plt.tight_layout()
plt.subplots_adjust(wspace=0.0, hspace=0.0)

plt.savefig('MEP_theta_MD.pdf')