#!/usr/bin/env python
import matplotlib.pyplot as plt
import matplotlib as mpl
import numpy as np

#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]]
)

TS_PBE = np.array(
[[1.2855, 4.4576, 132.99,  -2.08],
[1.2856, 4.4151, 133.01,  -2.10],
[1.2856, 4.3724, 133.00,  -2.20],
[1.2858, 4.3294, 133.00,  -2.27],
[1.2859, 4.2859, 133.00,  -2.35],
[1.2861, 4.2420, 133.01,  -2.44],
[1.2862, 4.1974, 133.01,  -2.53],
[1.2864, 4.1520, 132.99,  -2.62],
[1.2866, 4.1057, 132.97,  -2.72],
[1.2869, 4.0585, 132.99,  -2.83],
[1.2870, 4.0099, 133.06,  -2.95],
[1.2872, 3.9600, 133.10,  -3.07],
[1.2877, 3.9087, 133.01,  -3.19],
[1.2885, 3.8559, 132.76,  -3.29],
[1.2893, 3.8010, 132.63,  -3.40],
[1.2899, 3.7438, 132.64,  -3.52],
[1.2904, 3.6841, 132.75,  -3.65],
[1.2905, 3.6212, 133.03,  -3.77],
[1.2911, 3.5555, 133.22,  -3.84],
[1.2925, 3.4868, 133.18,  -3.84],
[1.2943, 3.4143, 133.11,  -3.74],
[1.2976, 3.3389, 132.99,  -3.48],
[1.3022, 3.2596, 132.94,  -2.99],
[1.3053, 3.1732, 133.16,  -2.10],
[1.3022, 3.0726, 127.99,  -0.75],
[1.3130, 2.9807, 123.65,   1.15],
[1.3271, 2.8856, 123.47,   4.40],
[1.3405, 2.7809, 122.09,   9.14],
[1.3644, 2.6826, 121.46,  15.94],
[1.4053, 2.6040, 122.78,  24.46],
[1.4475, 2.5206, 125.78,  34.88],
[1.5977, 2.6501, 129.38,  39.25],
[1.6745, 2.4567, 133.48,  62.76],
[1.7477, 2.4330, 134.34,  67.89],
[1.8209, 2.4099, 134.13,  70.82],
[1.8962, 2.3943, 132.96,  71.89],
[1.9707, 2.3772, 131.94,  71.71],
[2.0452, 2.3602, 130.11,  70.48],
[2.1205, 2.3476, 128.59,  68.77],
[2.1962, 2.3385, 127.25,  66.85],
[2.2728, 2.3385, 126.37,  64.68],
[2.3488, 2.3378, 125.23,  62.61],
[2.4243, 2.3328, 124.02,  60.57]]
)

top_PBE = np.array(
[[1.2925, 3.4578, 133.40,  -1.67],
[1.2945, 3.4157, 133.16,  -1.39],
[1.2971, 3.3736, 133.11,  -1.08],
[1.2994, 3.3313, 133.19,  -0.68],
[1.2990, 3.2882, 132.96,  -0.22],
[1.2948, 3.2438, 132.20,   0.37],
[1.2880, 3.1978, 129.09,   1.01],
[1.2821, 3.1508, 120.72,   1.34],
[1.2837, 3.1048, 114.09,   1.67],
[1.2875, 3.0587, 109.79,   2.28],
[1.2877, 3.0102, 105.14,   3.16],
[1.2917, 2.9620, 101.58,   4.34],
[1.2923, 2.9110, 100.76,   6.06],
[1.2956, 2.8596, 96.73,   7.95],
[1.2992, 2.8067, 97.04,  10.56],
[1.3037, 2.7524, 92.08,  13.49],
[1.3096, 2.6970, 94.21,  17.46],
[1.3183, 2.6415, 93.63,  22.24],
[1.3264, 2.5831, 91.29,  28.04],
[1.3413, 2.5277, 89.37,  35.31],
[1.3617, 2.4750, 94.41,  44.59],
[1.3868, 2.4250, 94.61,  55.76],
[1.4401, 2.4025, 103.34,  68.74],
[1.6156, 2.5178, 131.14,  77.86],
[1.6945, 2.5400, 133.07,  83.33],
[1.7604, 2.5534, 134.50,  87.01],
[1.8159, 2.5584, 135.13,  89.38],
[1.8642, 2.5573, 135.36,  90.71],
[1.9079, 2.5524, 135.32,  91.21],
[1.9484, 2.5446, 135.13,  91.02],
[1.9865, 2.5342, 134.83,  90.28],
[2.0227, 2.5214, 134.54,  89.06],
[2.0581, 2.5075, 134.41,  87.45],
[2.0936, 2.4940, 134.01,  85.54],
[2.1290, 2.4806, 133.59,  83.42],
[2.1646, 2.4678, 133.30,  81.14],
[2.2007, 2.4563, 133.10,  78.79],
[2.2373, 2.4464, 132.94,  76.45],
[2.2741, 2.4373, 132.75,  74.17],
[2.3112, 2.4292, 132.36,  72.00],
[2.3486, 2.4228, 131.76,  69.96],
[2.3862, 2.4175, 130.97,  68.10],
[2.4239, 2.4119, 130.08,  66.41],
[2.4618, 2.4067, 128.78,  64.81]]
)

bridge_PBE = np.array(
[[1.2934, 3.4579, 133.06,  -3.38],
[1.2945, 3.4157, 132.98,  -3.27],
[1.2952, 3.3734, 132.89,  -3.16],
[1.2965, 3.3309, 133.00,  -2.93],
[1.2978, 3.2880, 133.00,  -2.64],
[1.2979, 3.2445, 132.96,  -2.23],
[1.2972, 3.2001, 132.84,  -1.69],
[1.2966, 3.1549, 130.92,  -1.14],
[1.2960, 3.1088, 127.79,  -0.56],
[1.2965, 3.0620, 125.11,   0.20],
[1.2983, 3.0145, 121.64,   1.09],
[1.3011, 2.9662, 118.68,   2.26],
[1.3045, 2.9169, 118.73,   3.91],
[1.3073, 2.8658, 117.23,   5.80],
[1.3121, 2.8142, 114.80,   8.03],
[1.3175, 2.7611, 112.10,  10.78],
[1.3269, 2.7087, 111.48,  14.29],
[1.3364, 2.6546, 111.55,  18.49],
[1.3615, 2.6105, 112.37,  23.29],
[1.3875, 2.5665, 114.47,  28.38],
[1.4090, 2.5176, 117.18,  34.16],
[1.4422, 2.4785, 118.81,  40.79],
[1.4949, 2.4592, 119.57,  47.64],
[1.5544, 2.4498, 121.87,  53.61],
[1.6220, 2.4536, 124.01,  58.53],
[1.6830, 2.4543, 124.81,  62.01],
[1.7409, 2.4552, 124.18,  64.42],
[1.7911, 2.4489, 124.35,  66.38],
[1.8402, 2.3712, 121.53,  72.73],
[1.9622, 2.3881, 130.82,  78.19],
[1.9816, 2.3893, 130.90,  78.37],
[1.9970, 2.3838, 130.62,  78.45],
[2.0101, 2.3734, 130.09,  78.45],
[2.0390, 2.4647, 133.42,  76.44],
[2.0696, 2.4346, 130.91,  74.49],
[2.1068, 2.4197, 130.51,  73.03],
[2.1471, 2.4139, 130.12,  71.15],
[2.1860, 2.4050, 129.91,  69.36],
[2.2253, 2.3983, 129.13,  67.66],
[2.2644, 2.3914, 127.80,  65.97],
[2.3033, 2.3846, 126.72,  64.41],
[2.3426, 2.3804, 125.77,  63.00],
[2.3818, 2.3757, 124.96,  61.72],
[2.4211, 2.3722, 124.25,  60.57],
[2.4605, 2.3696, 123.60,  59.55]]
)

[1.7772, 2.3622, 122.82,  72.01],
[1.8282, 2.3870, 121.89,  72.30],
[1.8401, 2.3711, 121.53,  72.73],
[1.8538, 2.3569, 120.84,  73.01],
[1.8691, 2.3445, 118.75,  72.68],
[1.8846, 2.3324, 116.87,  72.47],
[1.8939, 2.3101, 116.87,  73.28],


fcc_PBE = np.array(
[[1.2942, 3.4579, 132.81,  -2.97],
[1.2950, 3.4157, 132.89,  -2.68],
[1.2959, 3.3734, 132.91,  -2.37],
[1.2971, 3.3309, 132.94,  -1.96],
[1.2981, 3.2881, 132.97,  -1.48],
[1.2997, 3.2449, 133.00,  -0.88],
[1.3008, 3.2010, 133.01,  -0.17],
[1.3025, 3.1566, 132.99,   0.69],
[1.3045, 3.1116, 132.93,   1.71],
[1.3071, 3.0658, 132.78,   2.94],
[1.3103, 3.0193, 132.72,   4.39],
[1.3141, 2.9720, 132.86,   6.13],
[1.3190, 2.9240, 133.31,   8.23],
[1.3247, 2.8751, 133.48,  10.70],
[1.3299, 2.8244, 132.95,  13.55],
[1.3324, 2.7704, 131.49,  16.88],
[1.3382, 2.7164, 130.68,  20.94],
[1.3471, 2.6624, 130.06,  25.76],
[1.3666, 2.6145, 130.73,  31.39],
[1.3882, 2.5671, 130.94,  37.74],
[1.4123, 2.5207, 130.62,  44.83],
[1.4461, 2.4822, 132.20,  52.74],
[1.6315, 2.4649, 137.32,  74.43],
[1.6984, 2.4740, 137.81,  79.59],
[1.7544, 2.4738, 138.35,  83.59],
[1.8052, 2.4699, 137.68,  86.63],
[1.8531, 2.4647, 137.41,  88.77],
[1.8961, 2.4540, 136.28,  90.17],
[1.9354, 2.4382, 135.29,  90.84],
[1.9719, 2.4173, 134.39,  90.82],
[2.0068, 2.3922, 133.31,  90.19],
[2.0412, 2.3645, 132.62,  88.95],
[2.0757, 2.3342, 131.64,  87.13],
[2.1112, 2.3033, 130.01,  84.81],
[2.1488, 2.2751, 129.01,  82.14],
[2.1886, 2.2511, 128.29,  79.31],
[2.2305, 2.2321, 126.83,  76.53],
[2.2750, 2.2239, 126.30,  74.10],
[2.3206, 2.2235, 126.15,  72.34],
[2.3675, 2.2390, 125.23,  71.39],
[2.4133, 2.2596, 125.72,  71.75],
[2.4551, 2.2150, 118.06,  72.23]]
)

def reactioncoordinate( r, Z, E ):
	if len(r) != len(Z) != len(E):
		print('lengths of r, Z and E 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[np.argmax(E)]
	for i in range( len(s) ):
		s[i] -= ds
	
	return s

mpl.rc("text", usetex=True)
#fig = plt.figure(figsize=(3.69,3))
Nrows = 2
Ncols = 2
fig, ax = plt.subplots(nrows=Nrows, ncols=Ncols, figsize=(3.69,4.))

colorcycle = plt.rcParams['axes.prop_cycle'].by_key()['color']

plt.subplot(Nrows,Ncols,1)

xaxis = reactioncoordinate( TS[:,0], TS[:,1], TS[:,3] )
plt.plot( xaxis, TS[:,3], marker='o', label='TS', mfc='None', markeredgecolor=colorcycle[0] )

xaxis = reactioncoordinate( top[:,0], top[:,1], top[:,3] )
plt.plot( xaxis, top[:,3], marker='o', label='Top', mfc='None', markeredgecolor=colorcycle[1] )

xaxis = reactioncoordinate( fcc[:,0], fcc[:,1], fcc[:,3] )
plt.plot( xaxis, fcc[:,3], marker='o', label='Fcc', mfc='None', markeredgecolor=colorcycle[2] )

xaxis = reactioncoordinate( bridge[:,0], bridge[:,1], bridge[:,3] )
plt.plot( xaxis, bridge[:,3], marker='o', label='Bridge', mfc='None', markeredgecolor=colorcycle[3] )

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

plt.annotate('(a)', xy=(0.05, 0.9), xycoords='axes fraction', size=12)

plt.title('MS-RPBEl')

plt.xlim(-1.5,0.5)
plt.ylim(-10.,150.)
plt.tick_params(length=6, width=1, direction='in', top=True, right=True, labelbottom=False)
plt.ylabel(r'Potential energy (kJ/mol)')
#plt.xlabel(r'Reaction coordinate (\r{A})')

plt.subplot(Nrows,Ncols,2)

xaxis = reactioncoordinate( TS_PBE[:,0], TS_PBE[:,1], TS_PBE[:,3] )
plt.plot( xaxis, TS_PBE[:,3], marker='o', label='TS', mfc='None', markeredgecolor=colorcycle[0] )

xaxis = reactioncoordinate( top_PBE[:,0], top_PBE[:,1], top_PBE[:,3] )
plt.plot( xaxis, top_PBE[:,3], marker='o', label='Top', mfc='None', markeredgecolor=colorcycle[1] )

xaxis = reactioncoordinate( fcc_PBE[:,0], fcc_PBE[:,1], fcc_PBE[:,3] )
plt.plot( xaxis, fcc_PBE[:,3], marker='o', label='Fcc', mfc='None', markeredgecolor=colorcycle[2] )

xaxis = reactioncoordinate( bridge_PBE[:,0], bridge_PBE[:,1], bridge_PBE[:,3] )
plt.plot( xaxis, bridge_PBE[:,3], marker='o', label='Bridge', mfc='None', markeredgecolor=colorcycle[3] )

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

plt.annotate('(b)', xy=(0.8, 0.9), xycoords='axes fraction', size=12)

plt.title('PBE')

#plt.legend(loc='upper left', numpoints=1, handletextpad=0.5, borderaxespad=0.2, frameon=False)
plt.legend(loc='best', numpoints=1, frameon=True, borderaxespad=0.2, handletextpad=0.4, fontsize=9)
plt.xlim(-1.5,0.5)
plt.ylim(-10.,150.)
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 coordinate (\r{A})')

plt.subplot(Nrows,Ncols,3)

xaxis = reactioncoordinate( TS[:,0], TS[:,1], TS[:,3] )
plt.plot( xaxis, TS[:,2], marker='o', mfc='None', label='TS', markeredgecolor=colorcycle[0] )

xaxis = reactioncoordinate( top[:,0], top[:,1], top[:,3] )
plt.plot( xaxis, top[:,2], marker='o', mfc='None', label='Top', markeredgecolor=colorcycle[1] )

xaxis = reactioncoordinate( fcc[:,0], fcc[:,1], fcc[:,3] )
plt.plot( xaxis, fcc[:,2], marker='o', mfc='None', label='Fcc', markeredgecolor=colorcycle[2] )

xaxis = reactioncoordinate( bridge[:,0], bridge[:,1], bridge[:,3] )
plt.plot( xaxis, bridge[:,2], marker='o', mfc='None', label='Bridge', markeredgecolor=colorcycle[3] )

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

plt.annotate('(c)', xy=(0.05, 0.9), xycoords='axes fraction', size=12)

plt.xlim(-1.5,0.5)
plt.ylim(80.,145.)
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,4)

xaxis = reactioncoordinate( TS_PBE[:,0], TS_PBE[:,1], TS_PBE[:,3] )
plt.plot( xaxis, TS_PBE[:,2], marker='o', label='', mfc='None', markeredgecolor=colorcycle[0] )

xaxis = reactioncoordinate( top_PBE[:,0], top_PBE[:,1], top_PBE[:,3] )
plt.plot( xaxis, top_PBE[:,2], marker='o', label='', mfc='None', markeredgecolor=colorcycle[1] )

xaxis = reactioncoordinate( fcc_PBE[:,0], fcc_PBE[:,1], fcc_PBE[:,3] )
plt.plot( xaxis, fcc_PBE[:,2], marker='o', label='', mfc='None', markeredgecolor=colorcycle[2] )

xaxis = reactioncoordinate( bridge_PBE[:,0], bridge_PBE[:,1], bridge_PBE[:,3] )
plt.plot( xaxis, bridge_PBE[:,2], marker='o', label='', mfc='None', markeredgecolor=colorcycle[3] )

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

plt.annotate('(d)', xy=(0.8, 0.9), xycoords='axes fraction', size=12)

plt.xlim(-1.5,0.5)
plt.ylim(80.,145.)
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.legend(loc='upper left', numpoints=1, handletextpad=0.5, borderaxespad=0.2, frameon=False)
#plt.xlim(0.,1.5)

plt.tight_layout()
plt.subplots_adjust(wspace=0, hspace=0)
plt.savefig('theta_mep_PBE.pdf')
#plt.show()