DateLong

Mutli-Temp TM Estimation

Transition Matrix Estimation from Multiple Temperatures

John D. Chodera

Jan-Hendrik Prinz

With growing computational ressources Parallel Tempering Simulations have been a useful and effective way of sampling the phase space for a variety of physical systems. By allowing the system to heat up and cool down, the number of activated transitions sampled by the simulation can be enriched while preserving the equilibrium distribution.  While the replica trajectories are unphysical due to the temperature swaps, between exchange attempts, valid dynamical trajectories are sampled from the appropriate temperature, producing an ensemble of trajectories at all temperatures. Recently, it has been demonstrated that a useful Markov model can be effectively constructed from the ensemble of trajectories generated at each temperature. We extend this approach by reweighting trajectories to make use of all data generated by the simulation, providing a more accurate model at any temperature, even those not sampled. We demonstrate this approach for the alanine dipeptide.

.INITIALIZATION

ReadME

During Initialization all files will be read and this may take a while. If you don want to, select no, when asked for initialization!

Path to images

Change this to the absolute path, where your repository is before you run any cells. Otherwise the datafiles will not be found!

imgPath = "/Users/jan-hendrikprinz/Studium/Projekte/dynamical-reweighting/manuscripts/markov-model-from-parallel-tempering/fig/";
dataPath = "/Users/jan-hendrikprinz/Studium/Projekte/dynamical-reweighting/datasets/alanine-dipeptide/";

Additional Packages

Here go all additional packages to be loaded on startup

<<TransitionMatrix/TransitionMatrix.m
Needs["ErrorBarPlots`"];

Options

All necessary options to be set. Especially for plots

Plots

SetOptions[Plot,BaseStyle->{11,FontFamily->"Palatino"}];
SetOptions[ListPlot,BaseStyle->{11,FontFamily->"Palatino"}];
SetOptions[Histogram,BaseStyle->{11,FontFamily->"Palatino"}];

Functions

All necessary function initialization

Normalizing

getLog10GaussMax[mean_,std_]:=Log[10,1/2(mean+Sqrt[mean^2+4std^2])];
getLog10DirMax[n1_,n2_]:=Log[10,((1+n1)/(1+n1+n2))]

getLog10GaussNorm[mean_, std_] :=2/( (E^(-((((mean-(Sqrt[(mean^2)+4 (std^2)])))^2)/(8 (std^2))))) (mean+(Sqrt[(mean^2)+4 (std^2)])));

getLog10DirNorm[n1_,n2_]:=1/((((((1+n1)/(1+n1+n2))))^(1+n1)) ((((n2/(1+n1+n2))))^n2))

Plots

Necessary Functions for Plotting

LimitAngle[s_]:=Map[Min[Max[#,-180],180]&,s];

StateShape[s_,e_]:=Module[{ns,ne},
ns=s+If[s[[1]]>s[[2]],{0,360},0];
ne=e+If[e[[1]]>e[[2]],{0,360},0];
Map[Rectangle[LimitAngle[{ns[[1]],ne[[1]]}+#],LimitAngle[{ns[[2]],ne[[2]]}+#]]&,{{0,0},{0,-360},{-360,0},{-360,-360}}]
]

getHistogram[data_,range_,bins_]:=Module[{},
result=Table[0,{bins}];
Map[If[#≥0&&#<bins,result[[#+1]]++]&,Map[Floor[(#-range[[1]])/(range[[2]]-range[[1]])*bins]&,data]];
Transpose[{Range[range[[1]]+(range[[2]]-range[[1]])/bins/2,range[[2]],(range[[2]]-range[[1]])/bins],result}]
]
PlotHistogram[data_,range_,bins_]:=Module[{thisNorm,thisHistogram},
thisHistogram=getHistogram[data,range,bins];
thisNorm=Apply[Plus,thisHistogram[[All,2]]]*((range[[2]]-range[[1]])/bins);
Map[Rectangle[#/{1,thisNorm}-{(range[[2]]-range[[1]])/bins/2,0},{#[[1]]+(range[[2]]-range[[1]])/bins/2,0}]&,
thisHistogram]
]
PlotHistogramOutline[data_,range_,bins_]:=Module[{thisNorm,thisHistogram},
thisHistogram=getHistogram[data,range,bins];
thisNorm=Apply[Plus,thisHistogram[[All,2]]]*((range[[2]]-range[[1]])/bins);
Line[Flatten[Map[{#/{1,thisNorm}-{(range[[2]]-range[[1]])/bins/2,0},#/{1,thisNorm}+{(range[[2]]-range[[1]])/bins/2,0}}&,
thisHistogram],{1,2}]]
]
PlotHistogramOutlineSmooth[data_,range_,bins_]:=Module[{thisHistogram,thisFunction,binSize,thisKernelSize,thisKernel},
thisHistogram=getHistogram[data,range,bins];
binSize=(range[[2]]-range[[1]])/bins;
thisKernelSize=3;
thisKernel=Table[1./Sqrt[2π binSize^2]*Exp[-(x)^2/(2 binSize^2)],{x,-thisKernelSize*binSize,thisKernelSize*binSize,binSize}];
Line[Transpose[{thisHistogram[[All,1]],1/Apply[Plus,thisHistogram[[All,2]]]*ListConvolve[thisKernel,thisHistogram[[All,2]],{1,-1},0][[Range[thisKernelSize+1,thisKernelSize+Length[thisHistogram]]]]}]
]]
GetHistogramOutlineSmooth[data_,range_,bins_]:=Module[{thisHistogram,thisFunction,binSize,thisKernelSize,thisKernel},
thisHistogram=getHistogram[data,range,bins];
binSize=(range[[2]]-range[[1]])/bins;
thisKernelSize=3;
thisKernel=Table[1./Sqrt[2π binSize^2]*Exp[-(x)^2/(2 binSize^2)],{x,-thisKernelSize*binSize,thisKernelSize*binSize,binSize}];
Transpose[{thisHistogram[[All,1]],1/Apply[Plus,thisHistogram[[All,2]]]*ListConvolve[thisKernel,thisHistogram[[All,2]],{1,-1},0][[Range[thisKernelSize+1,thisKernelSize+Length[thisHistogram]]]]}]
]
PlotHistogramAbsolute[data_,range_,bins_]:=Module[{thisMax,thisHistogram},
thisHistogram=getHistogram[data,range,bins];
thisMax=Max[thisHistogram[[All,2]]];
Map[Rectangle[#/{1,1}-{(range[[2]]-range[[1]])/bins/2,0},{#[[1]]+(range[[2]]-range[[1]])/bins/2,0}]&,
thisHistogram]
]

ErrorBarVertical[p_,e_,w_]:={Line[{p-{0,e},p+{0,e}}],Line[{p-{-w,e},p-{w,e}}],Line[{p+{-w,e},p+{w,e}}]};

Algebra

Some helpers for Algebraic Computaitons

getEquilibriumDistribution[matrix_]:=(#/Apply[Plus,#])&[Eigenvectors[Transpose[matrix]][[1]]];

getSymmetrizedEV[matrix_]:=Module[{eqDist},
eqDist=(#/Apply[Plus,#])&[Eigenvectors[Transpose[matrix]][[1]]];
Eigenvectors[matrix.DiagonalMatrix[1/eqDist]]
]

getSymmetrizedTransitionMatrix[matrix_]:=Module[{eqDist},
eqDist=(#/Apply[Plus,#])&[Eigenvectors[Transpose[matrix]][[1]]];
matrix.DiagonalMatrix[1/eqDist]
]

Definitions

numberOfStates=6;

.DATA INPUT

$DefaultFont={"Palatino",Medium}

Data from Parallel Tempering but only at Single Temperatures (ST / UT)

Some Constants

numberOfSamples = 10000;
numberOfTemperatures = 40;

Read all transition matrices for single temperatures

tMatricesST=Table[Table[0,{10000},{36}],{numberOfTemperatures}];
Table[tMatricesST[[temp+1]]=Import[dataPath<>"parallel-tempering/sampled-matrices/tMatrixListSP"<>ToString[temp]<>".tsv"],{temp,0,numberOfTemperatures-1}];
tMatricesST6by6=Table[Map[Partition[#,numberOfStates]&,tMatricesST[[temp]]],{temp,1,numberOfTemperatures}];

Read temperaturesfiles in Kelvin

temperatureListST = Import[dataPath<>"parallel-tempering/temps", "Table"][[1]];

Confidence Intervals

confidence95ST=Table[Sort[tMatricesST6by6[[temp,All,ii,jj]]][[{numberOfSamples*0.025,numberOfSamples*0.975}]],{temp,1,numberOfTemperatures},{ii,1,numberOfStates},{jj,1,numberOfStates}];

confidence68ST=Table[Sort[tMatricesST6by6[[temp,All,ii,jj]]][[{numberOfSamples*0.165,numberOfSamples*0.835}]],{temp,1,numberOfTemperatures},{ii,1,numberOfStates},{jj,1,numberOfStates}];

Mean values

meanST=Table[Mean[tMatricesST6by6[[temp,All,ii,jj]]],{temp,1,numberOfTemperatures},{ii,1,6},{jj,1,6}];

varianceST = 0.25(confidence95ST[[All,All,All,2]]-confidence95ST[[All,All,All,1]]);

Confidence Intervals using only parts of the numberOfSamples data. Ranging from (10^( 0.5))to (10^( 4)) numberOfSamples used.

(* DON'T RUN, is not necessary any more and takes a lot of time ! *)

confidence95STCheck=Table[{{Ceiling[10^numberOfSamplesExp],Sort[tMatricesST6by6[[temp,Range[Ceiling[10^numberOfSamplesExp]],ii,jj]]][[Ceiling[10^numberOfSamplesExp*0.025]]]},{Ceiling[10^numberOfSamplesExp],Sort[tMatricesST6by6[[temp,Range[Ceiling[10^numberOfSamplesExp]],ii,jj]]][[Ceiling[10^numberOfSamplesExp*0.975]]]}},{temp,1,numberOfTemperatures},{ii,1,numberOfStates},{jj,1,numberOfStates},{numberOfSamplesExp,0.5,4,0.1}];

confidence68STCheck=Table[{{Ceiling[10^numberOfSamplesExp],Sort[tMatricesST6by6[[temp,Range[Ceiling[10^numberOfSamplesExp]],ii,jj]]][[Ceiling[10^numberOfSamplesExp*0.165]]]},{Ceiling[10^numberOfSamplesExp],Sort[tMatricesST6by6[[temp,Range[Ceiling[10^numberOfSamplesExp]],ii,jj]]][[Ceiling[10^numberOfSamplesExp*0.835]]]}},{temp,1,numberOfTemperatures},{ii,1,numberOfStates},{jj,1,numberOfStates},{numberOfSamplesExp,0.5,4,0.1}];

Read transition Count Matrices for Reference

tCountMatricesUT=Import[dataPath<>"parallel-tempering/temperature_transition_matrices","Table"];
tCountMatricesUT=tCountMatricesUT[[All,Range[2,7]]]*501;
tCountMatricesUT6by6=Partition[tCountMatricesUT,6];
tMatricesUT6by6=Map[#/Apply[Plus,#]&,tCountMatricesUT6by6,{2}];

tCountMatricesUT6by6[[6]]//MatrixForm

Data from Multi Temperature Estimates (MT)

Import Data

dataMT=Import[dataPath<>"parallel-tempering/combined.dat"];

Read Temperature List

temperatureListMT= dataMT[[All, 1]];

Mean

meanMT = Table[Partition[dataMT[[temp, Range[2, numberOfStates^2*3+1, 3]]], numberOfStates], {temp, 1, Length[temperatureListMT]}];

Confidence Intervals

confidence68MT = Table[Partition[dataMT[[temp, Range[3, numberOfStates^2*3+1, 3]]], numberOfStates], {temp, 1, Length[temperatureListMT]}];
confidence95MT = Table[Partition[dataMT[[temp, Range[4, numberOfStates^2*3+1, 3]]], numberOfStates], {temp, 1, Length[temperatureListMT]}];

Read data for Comparison of  Contributions of differerent temperatures

dataMTNames= Import[dataPath<>"parallel-tempering/temperatureContributions.nc","Datasets"]
dataMTAnnotations= Import[dataPath<>"parallel-tempering/temperatureContributions.nc","Annotations"]
dataMTAll=Import[dataPath<>"parallel-tempering/temperatureContributions.nc","Data"];

Read List of Total Energies

totalEnergyListMT=dataMTAll[[3]];

Read log(W[i,j,k,l])

logWijkl=dataMTAll[[10]];
Dimensions[logWijkl]

logWkl=dataMTAll[[12]];
Dimensions[logWkl]

meanMT2=dataMTAll[[14]];

varianceMT=dataMTAll[[15]];

Aijkl=dataMTAll[[4]];

Data of Comparison Trajectories only at 302 Kelvin from each state (RT)

Read Matlab File

dataRT=Import[dataPath<>"shoot10ps302K/a1e-shoot-10ps.mat"];

Construct Transition Count Matrix

tCountMatrixRT6by6=Table[0,{numberOfStates},{numberOfStates}];
Table[tCountMatrixRT6by6[[nn,dataRT[[13,61,nn,ii]]]]++,{nn,1,6},{ii,1,10000}];

Show Transition Count Matrix and Transition Matix

tCountMatrixRT6by6//MatrixForm
tMatrixRT6by6=Map[#/Apply[Plus,#]&,tCountMatrixRT6by6];
N[Transpose[tMatrixRT6by6]]//MatrixForm

Export[dataPath<>"shoot10ps302K/transitionCountMatrix.dat",tCountMatrixRT6by6];

tMatricesRT=Table[0,{10000},{36}];
tMatricesRT=Import[dataPath<>"shoot10ps302K/tMatrixList.tsv"];
tMatricesRT6by6=Map[Partition[#,numberOfStates]&,tMatricesRT];

Confidence Intervals

confidence95RT=Table[Sort[tMatricesRT6by6[[All,ii,jj]]][[{numberOfSamples*0.025,numberOfSamples*0.975}]],{ii,1,numberOfStates},{jj,1,numberOfStates}];

confidence68RT=Table[Sort[tMatricesRT6by6[[All,ii,jj]]][[{numberOfSamples*0.165,numberOfSamples*0.835}]],{ii,1,numberOfStates},{jj,1,numberOfStates}];

Mean values

meanRT=Table[Mean[tMatricesRT6by6[[All,ii,jj]]],{ii,1,6},{jj,1,6}];

varianceRT = 0.5(confidence68RT[[All,All,2]]-confidence68RT[[All,All,1]]);

Data of Hummer Sampled Matrices (from John)

numberOfSampledHummer=100;
numberOfTempHummer=40;
temperatureListHT=temperatureListST[[Range[numberOfTempHummer]]];
lagtimeHummer=6;

tMatricesQTRate=Table[Import[dataPath<>"parallel-tempering/hummerRateMatrixEstimates/output/"<>"302K.Kij."<>ToString[nn],"Table"],{nn,1,numberOfSampledHummer}];

meanQTRate=Map[Mean,Transpose[tMatricesQTRate,{2,3,1}]];
varianceQTRate=Map[Variance,Transpose[tMatricesQTRate,{2,3,1}]];

tMatricesQT=Map[MatrixExp[#*lagtimeHummer]&,tMatricesQTRate,{1}];

meanQT=Map[Mean,Transpose[tMatricesQT,{2,3,1}]];
varianceQT=Map[Variance,Transpose[tMatricesQT,{2,3,1}]];

confidence95QT=Table[Sort[tMatricesQT[[All,ii,jj]]][[{numberOfSampledHummer*0.03,numberOfSampledHummer*0.97}]],{ii,1,numberOfStates},{jj,1,numberOfStates}];

confidence67QT=Table[Sort[tMatricesQT[[All,ii,jj]]][[{numberOfSampledHummer*0.17,numberOfSampledHummer*0.83}]],{ii,1,numberOfStates},{jj,1,numberOfStates}];

varianceQT = 0.5(confidence67QT[[All,All,2]]-confidence67QT[[All,All,1]]);

tMatricesHTRate=Table[Import[dataPath<>"parallel-tempering/hummerRateMatrixEstimates/output/"<>ToString[temp]<>".Kij."<>ToString[nn],"Table"],{temp,0,numberOfTempHummer-1},{nn,1,numberOfSampledHummer}];

meanHTRate=Map[Mean,tMatricesHTRate,{1}];
varianceHTRate=Map[Variance,tMatricesHTRate,{1}];

tMatricesHT=Map[MatrixExp[#*lagtimeHummer]&,tMatricesHTRate,{2}];

meanHT=Map[Mean,tMatricesHT,{1}];
varianceHT=Map[Variance,tMatricesHT,{1}];

confidence95HT=Table[Sort[tMatricesHT[[temp,All,ii,jj]]][[{numberOfSampledHummer*0.03,numberOfSampledHummer*0.97}]],{temp,1,numberOfTempHummer},{ii,1,numberOfStates},{jj,1,numberOfStates}];

confidence67HT=Table[Sort[tMatricesHT[[temp,All,ii,jj]]][[{numberOfSampledHummer*0.17,numberOfSampledHummer*0.83}]],{temp,1,numberOfTempHummer},{ii,1,numberOfStates},{jj,1,numberOfStates}];

varianceHT = 0.25(confidence95HT[[All,All,All,2]]-confidence95HT[[All,All,All,1]]);

tRateMatrix302KHT6by6=(Subscript[ℝ, 6])(({{-0.7621856560088285`, 0.6974454846930033`, 0.0607054766370334`, 0.004034694678791635`, 0.`, 0.`}, {0.7856922476690512`, -0.7954525886866808`, 0.`, 0.008744875677034301`, 0.`, 0.001015465340595308`}, {0.378919065041949`, 0.`, -1.3351025799221525`, 0.9561835148802035`, 0.`, 0.`}, {0.019432468195749963`, 0.03738769822567313`, 0.7378022211766183`, -0.7946223875980415`, 0.`, 0.`}, {0.`, 0.`, 0.`, 0.`, -2300.302300864345`, 2300.302300864345`}, {0.`, 0.3691498746847238`, 0.`, 0.`, 212.72769893452744`, -213.09684880921213`}}));

tMHummer302K=ToTransitionMatrix[Chop[tRateMatrix302KHT6by6],1];

Equilibrium[tRateMatrix302KHT6by6]

cM

Likelihood[tMHummer302K,CountMatrix[tCountMatricesUT6by6[[6]]]]

tMHummer302K

Equilibrium[tRateMatrix302KHT6by6]//CForm

Normal[tRateMatrix302KHT6by6]*Equilibrium[tRateMatrix302KHT6by6]//MatrixForm

Normal[tRateMatrix302KHT6by6]

Normal[tRateMatrix302KHT6by6]*Equilibrium[tRateMatrix302KHT6by6]

tCountMatricesUT6by6[[6]]

tMHummer302K

Data from Hummer Sampled Matrices (from Jan)

Some Constants

numberOfSamplesHT = 10000;
numberOfTemperaturesHT = 40;

Read temperaturesfiles in Kelvin

temperatureListHT = Import[dataPath<>"parallel-tempering/temps", "Table"][[1]];

Read all transition matrices for single temperatures

tMatricesHT=Table[Table[0,{numberOfSamplesHT},{36}],{numberOfTemperaturesHT}];
Table[tMatricesHT[[temp+1]]=Import[dataPath<>"parallel-tempering/sampled-matrices/rMatrixListSP"<>ToString[temp]<>".tsv"],{temp,0,Min[1000,numberOfTemperaturesHT-1]}];
tMatricesHT6by6=Table[Map[Partition[#,numberOfStates]&,tMatricesHT[[temp]]],{temp,1,numberOfTemperaturesHT}];

confidence95HT=Table[Sort[tMatricesHT6by6[[temp,All,ii,jj]]][[{numberOfSamples*0.025,numberOfSamples*0.975}]],{temp,1,numberOfTemperaturesHT},{ii,1,numberOfStates},{jj,1,numberOfStates}];

confidence68HT=Table[Sort[tMatricesHT6by6[[temp,All,ii,jj]]][[{numberOfSamples*0.165,numberOfSamples*0.835}]],{temp,1,numberOfTemperaturesHT},{ii,1,numberOfStates},{jj,1,numberOfStates}];

Mean values

meanHT=Table[Mean[tMatricesHT6by6[[temp,All,ii,jj]]],{temp,1,numberOfTemperaturesHT},{ii,1,6},{jj,1,6}];

varianceHT = 0.25(confidence95HT[[All,All,All,2]]-confidence95HT[[All,All,All,1]]);

Confidence Intervals using only parts of the numberOfSamples data. Ranging from (10^( 0.5))to (10^( 4)) numberOfSamples used.

(* DON'T RUN, is not necessary any more and takes a lot of time ! *)

confidence95HTCheck=Table[{{Ceiling[10^numberOfSamplesExp],Sort[tMatricesHT6by6[[temp,Range[Ceiling[10^numberOfSamplesExp]],ii,jj]]][[Ceiling[10^numberOfSamplesExp*0.025]]]},{Ceiling[10^numberOfSamplesExp],Sort[tMatricesHT6by6[[temp,Range[Ceiling[10^numberOfSamplesExp]],ii,jj]]][[Ceiling[10^numberOfSamplesExp*0.975]]]}},{temp,1,numberOfTemperatures},{ii,1,numberOfStates},{jj,1,numberOfStates},{numberOfSamplesExp,0.5,4,0.1}];

confidence68HTCheck=Table[{{Ceiling[10^numberOfSamplesExp],Sort[tMatricesHT6by6[[temp,Range[Ceiling[10^numberOfSamplesExp]],ii,jj]]][[Ceiling[10^numberOfSamplesExp*0.165]]]},{Ceiling[10^numberOfSamplesExp],Sort[tMatricesHT6by6[[temp,Range[Ceiling[10^numberOfSamplesExp]],ii,jj]]][[Ceiling[10^numberOfSamplesExp*0.835]]]}},{temp,1,numberOfTemperatures},{ii,1,numberOfStates},{jj,1,numberOfStates},{numberOfSamplesExp,0.5,4,0.1}];

Data of energy development (not read by default !)

energyListRaw=Import[dataPath<>"parallel-tempering/etot.out", "Table"];

energyList = Table[Partition[energyListRaw[[All, temp]], 200], {temp, 1, numberOfTemperatures}];

Compute Eigenvalues and Implied Timescales

eigenvaluesST=Map[Eigenvalues[#][[Range[1,numberOfStates]]]&,meanST];
impliedTimescalesST=Map[-6./Log[Abs[Eigenvalues[#][[Range[2,numberOfStates]]]]]&,meanST];

eigenvaluesMT=Map[Eigenvalues[#][[Range[1,numberOfStates]]]&,meanMT];
impliedTimescalesMT = Map[-6./Log[Abs[Eigenvalues[#][[Range[2, numberOfStates]]]]] &, meanMT];

eigenvaluesHT=Map[Eigenvalues[#][[Range[1,numberOfStates]]]&,meanHT];
impliedTimescalesHT = Map[-6./Log[Abs[Eigenvalues[#][[Range[2, numberOfStates]]]]] &, meanHT];

eigenvaluesRT=Eigenvalues[meanRT][[Range[1,numberOfStates]]];
impliedTimescalesRT=-6./Log[Abs[Eigenvalues[N[meanRT]][[Range[2,numberOfStates]]]]];

eigenvalueVarianceST =
  Table[
   leftEV = Eigenvectors[Transpose[meanST[[temp]]]];
   rightEV = Eigenvectors[meanST[[temp]]]; Table[Apply[Plus, Flatten[Outer[Times, leftEV[[ev]], rightEV[[ev]]]*varianceST[[temp, All, All]]]]/(leftEV[[ev]].rightEV[[ev]])
    , {ev, 1, numberOfStates}]
   , {temp, 1, numberOfTemperatures}];

eigenvalueVarianceHT =
  Table[
   leftEV = Eigenvectors[Transpose[meanHT[[temp]]]];
   rightEV = Eigenvectors[meanHT[[temp]]];Table[Apply[Plus, Flatten[Outer[Times, leftEV[[ev]], rightEV[[ev]]]*varianceHT[[temp, All, All]]]]/(leftEV[[ev]].rightEV[[ev]])
    , {ev, 1, numberOfStates}]
   , {temp, 1, numberOfTemperaturesHT}];

eigenvalueVarianceMT =
    Table[
      leftEV = Eigenvectors[Transpose[meanMT[[temp]]]];
      rightEV = Eigenvectors[meanMT[[temp]]];Table[Apply[Plus, Flatten[Outer[Times, leftEV[[ev]], rightEV[[ev]]]*confidence68MT[[temp, All, All]]]]/(leftEV[[ev]].rightEV[[ev]])
        , {ev, 1, numberOfStates}]
      , {temp, 1, Length[temperatureListMT]}];

leftEV = Eigenvectors[Transpose[meanRT]];
rightEV = Eigenvectors[meanRT];
eigenvalueVarianceRT =Table[Apply[Plus, Flatten[Outer[Times, leftEV[[ev]], rightEV[[ev]]]*varianceRT[[ All, All]]]]/(leftEV[[ev]].rightEV[[ev]])
        , {ev, 1, numberOfStates}];

varianceImpliedTimescalesST =(6.0/(Abs[eigenvaluesST ]*Log[Abs[eigenvaluesST ]]^2)*Abs[eigenvalueVarianceST])[[All,Range[2,numberOfStates]]];
varianceImpliedTimescalesHT =(6.0/(Abs[eigenvaluesHT ]*Log[Abs[eigenvaluesHT ]]^2)*Abs[eigenvalueVarianceHT])[[All,Range[2,numberOfStates]]];
varianceImpliedTimescalesMT = (6.0/(Abs[eigenvaluesMT ]*Log[Abs[eigenvaluesMT ]]^2)*Abs[eigenvalueVarianceMT])[[All,Range[2,numberOfStates]]];
varianceImpliedTimescalesRT = (6.0/(Abs[eigenvaluesRT ]*Log[Abs[eigenvaluesRT ]]^2)*Abs[eigenvalueVarianceRT])[[All,Range[2,numberOfStates]]];

Data of Torsions at 302K (not read by default !)

wMatrix302K=Import[dataPath<>"parallel-tempering/"<>ToString[6-1]<>".tors","Table"];
wMatrix302K=Drop[wMatrix302K,1];

wMatrix600K=Import[dataPath<>"parallel-tempering/"<>ToString[40-1]<>".tors","Table"];
wMatrix600K=Drop[wMatrix600K,1];

wMatrix400K=Import[dataPath<>"parallel-tempering/"<>ToString[20-1]<>".tors","Table"];
wMatrix400K=Drop[wMatrix400K,1];

Remove lots of memory

Remove[tMatricesST];
Remove[dataMT];
Remove[totalEnergyListMT];
(*Remove[tMatricesHT];*)
Remove[tMatricesHTRate];
Remove[tMatricesQT];
Remove[tMatricesQTRate];

.PLOTS

PlotFunction Definitions

plotDetailedScale=1;
plotDetailedMin=-3;
plotDetailedMax=0;

MeanSTPlot[st1_,st2_]:=Module[{},
ListPlot[Transpose[{temperatureListST,meanST[[All,st1,st2]]}],Joined->True,PlotStyle->{Thickness[Medium],Hue[0]}]
];
MeanMTPlot[st1_,st2_]:=Module[{},
ListPlot[Transpose[{temperatureListMT,meanMT[[All,st1,st2]]}],Joined->True,PlotStyle->{Thickness[Medium],RGBColor[0,0,1]}]
];
MeanHTPlot[st1_,st2_]:=Module[{},
ListPlot[Transpose[{temperatureListHT,meanHT[[All,st1,st2]]}],Joined->True,PlotStyle->{Thickness[Medium],RGBColor[0,0.7,0]}]
];
MeanRTPlot[st1_,st2_]:=Module[{},
ListPlot[Transpose[{temperatureListHT,meanRT[[All,st1,st2]]}],Joined->True,PlotStyle->{Thickness[Medium],RGBColor[0,0.7,0]}]
];

MaxLikeSTPlot[st1_,st2_]:=Module[{},ListPlot[Transpose[{temperatureListST,tMatricesUT6by6[[All,st1,st2]]}],Joined->True,PlotStyle->{Thick,RGBColor[0,0,0]}]
];

Confidence95STPlot[st1_,st2_]:=Module[{},
ListPlot[{Transpose[{temperatureListST,confidence95ST[[All,st1,st2,1]]}],Transpose[{temperatureListST,confidence95ST[[All,st1,st2,2]]}]},Filling->True,Joined->True,FillingStyle->{RGBColor[1,0,0,.2]},PlotStyle->RGBColor[1,0,0,0.7]]
];
Confidence95HTPlot[st1_,st2_]:=Module[{},
ListPlot[{Transpose[{temperatureListST[[Range[numberOfTemperaturesHT]]],confidence95HT[[All,st1,st2,1]]}],Transpose[{temperatureListST[[Range[numberOfTemperaturesHT]]],confidence95HT[[All,st1,st2,2]]}]},Filling->True,Joined->True,FillingStyle->{RGBColor[0,0.7,0,.2]},PlotStyle->{RGBColor[0,0.7,0,.7]}]
];
Confidence95MTPlot[st1_,st2_]:=Module[{},ListPlot[{Transpose[{temperatureListMT,meanMT[[All,st1,st2]]-confidence95MT[[All,st1,st2]]}],Transpose[{temperatureListMT,meanMT[[All,st1,st2]]+confidence95MT[[All,st1,st2]]}]},Filling->True,PlotStyle->RGBColor[0,0,1,.7],FillingStyle->RGBColor[0,0,1,.2],Joined->True]
];

MaxLikeRTPlot[st1_,st2_]:=Module[{},
ErrorListPlot[{{{302,meanRT[[st1,st2]]},ErrorBar[varianceRT[[st1,st2]]]}},PlotStyle->{Thick,GrayLevel[0]}]
];
MaxLikeRTPlotCross[st1_,st2_]:=Module[{},
Graphics[{Black,AbsoluteThickness[1.5],Circle[{302,meanRT[[st1,st2]]},Offset[{10,10}]],
Line[{Offset[{0,-10},{302,meanRT[[st1,st2]]}],Offset[{0,10},{302,meanRT[[st1,st2]]}]}],
Line[{Offset[{-10,0},{302,meanRT[[st1,st2]]}],Offset[{10,0},{302,meanRT[[st1,st2]]}]}]
}]
];
MaxLikeHTPlot[st1_,st2_]:=Module[{},
Graphics[{Thick,RGBColor[0,.7,0],Line[{{298,Normal[tMHummer302K][[st1,st2]]},{306,Normal[tMHummer302K][[st1,st2]]}}]}]
];

EigenvalueSTPlot[ev_]:=Module[{},
ListPlot[Transpose[{temperatureListST,eigenvaluesST[[All,ev]]}],Joined->True,PlotStyle->{{Thickness[Medium],Hue[0]}},PlotRange->All]
];
EigenvalueHTPlot[ev_]:=Module[{},
ListPlot[Transpose[{temperatureListHT,eigenvaluesHT[[All,ev]]}],Joined->True,PlotStyle->{{Thickness[Medium],RGBColor[0,0.7,0]}},PlotRange->All]
];
EigenvalueMTPlot[ev_]:=Module[{},
ListPlot[Transpose[{temperatureListMT,eigenvaluesMT[[All,ev]]}],Joined->True,PlotStyle->{{Thickness[Medium],RGBColor[0,0,1]}},PlotRange->All]
];
EigenvalueRTPlot[ev_]:=Module[{},ErrorListPlot[{{{302,eigenvaluesRT[[ev]]},ErrorBar[varianceImpliedTimescalesRT[[ev]]]}},PlotStyle->{Thickness[Medium],GrayLevel[0]}]
];
EigenvalueSTVariancePlot[ev_]:=Module[{},
ListPlot[{Transpose[{temperatureListST, eigenvaluesST[[All,ev]]-eigenvalueVarianceST[[All,ev]]}], Transpose[{temperatureListST, eigenvaluesST[[All,ev]]+eigenvalueVarianceST[[All,ev]]}]}, Filling -> True, Joined -> True,PlotRange->All, FillingStyle -> {RGBColor[1, 0, 0, .2]},PlotStyle->{RGBColor[1,0,0,.7]}]
];
EigenvalueHTVariancePlot[ev_]:=Module[{},
ListPlot[{Transpose[{temperatureListHT, eigenvaluesHT[[All,ev]]-eigenvalueVarianceHT[[All,ev]]}], Transpose[{temperatureListHT, eigenvaluesHT[[All,ev]]+eigenvalueVarianceHT[[All,ev]]}]}, Filling -> True, Joined -> True,PlotRange->All, FillingStyle -> {RGBColor[0, 0.7, 0, .1]},PlotStyle->{RGBColor[0,0.7,0,.7]}]
];
EigenvalueMTVariancePlot[ev_]:=Module[{},
ListPlot[{Transpose[{temperatureListMT, eigenvaluesMT[[All,ev]]-eigenvalueVarianceMT[[All,ev]]}], Transpose[{temperatureListMT, eigenvaluesMT[[All,ev]]+eigenvalueVarianceMT[[All,ev]]}]}, Filling -> True, Joined -> True,PlotRange->All, FillingStyle -> {RGBColor[0, 0, 1, .2]},PlotStyle->{RGBColor[0,0,1,.7]}]
];
EigenvalueRTPlotCross[ev_]:=Module[{size=6},
Graphics[{Black,AbsoluteThickness[1.0],Circle[{302,eigenvaluesRT[[ev]]},Offset[{size,size}]],
Line[{Offset[{0,-size},{302,eigenvaluesRT[[ev]]}],Offset[{0,size},{302,eigenvaluesRT[[ev]]}]}],
Line[{Offset[{-size,0},{302,eigenvaluesRT[[ev]]}],Offset[{size,0},{302,eigenvaluesRT[[ev]]}]}]
}]
];

meanSTDetailedPlotHistogram[st1_,st2_]:=Graphics[{RGBColor[1.0,0,0],EdgeForm[{Thin,GrayLevel[0]}],PlotHistogram[Log[10,tMatricesST6by6[[6,All,st1,st2]]],{plotDetailedMin,plotDetailedMax},120]}];
meanHTDetailedPlotHistogram[st1_,st2_]:=Graphics[{RGBColor[0,0.7,0,.8],EdgeForm[{Thin,GrayLevel[0]}],PlotHistogram[Log[10,tMatricesHT[[6,All,st1,st2]]],{plotDetailedMin,plotDetailedMax},60]}];
meanRTDetailedPlotHistogram[st1_,st2_]:=Graphics[{RGBColor[0.3,0.3,0.3,.1],EdgeForm[{Thin,GrayLevel[0]}],PlotHistogram[Log[10,tMatricesRT6by6[[All,st1,st2]]],{plotDetailedMin,plotDetailedMax},270]}];

meanSTDetailedPlotGauss[st1_,st2_]:=Plot[getLog10GaussNorm[meanST[[6,st1,st2]],varianceST[[6,st1,st2]]]*10^x*Exp[-(10^x-meanST[[6,st1,st2]])^2/(2 varianceST[[6,st1,st2]]^2)]*plotDetailedScale,{x,plotDetailedMin,plotDetailedMax},PlotRange->All,PlotStyle-> {Thick,RGBColor[1,0,0]},Filling->Axis,FillingStyle->RGBColor[1,0,0,0.1]];

meanHTDetailedPlotGauss[st1_,st2_]:=Plot[getLog10GaussNorm[meanHT[[6,st1,st2]],varianceHT[[6,st1,st2]]]*10^x*Exp[-(10^x-meanHT[[6,st1,st2]])^2/(2 varianceHT[[6,st1,st2]]^2)]*plotDetailedScale,{x,plotDetailedMin,plotDetailedMax},PlotRange->All,PlotStyle-> {Thick,RGBColor[0,0.7,0]},Filling->Axis,FillingStyle->RGBColor[0,0.7,0,0.1]];

meanMTDetailedPlotGauss[st1_,st2_]:=Plot[getLog10GaussNorm[meanMT[[11,st1,st2]],confidence68MT[[11,st1,st2]]]*10^x*Exp[-(10^x-meanMT[[11,st1,st2]])^2/(2 confidence68MT[[11,st1,st2]]^2)]*plotDetailedScale,{x,plotDetailedMin,plotDetailedMax},PlotRange->All,PlotStyle-> {Thick,RGBColor[0,0,1]},Filling->Axis,FillingStyle->RGBColor[0,0,1,0.1]];
meanSTDetailedPlotGauss[st1_,st2_]:=Module[{meanDirichletUT,maxDirichletUT,varDirichletUT},
meanDirichletUT=(tCountMatricesUT[[st1,st2]]+1)/(numberOfStates+Apply[Plus,tCountMatricesUT[[st1]]]);
maxDirichletUT=(tCountMatricesUT[[st1,st2]])/(Apply[Plus,tCountMatricesUT[[st1]]]);
varDirichletUT=((tCountMatricesUT[[st1,st2]]+1)(numberOfStates+Apply[Plus,tCountMatricesUT[[st1]]]-(tCountMatricesUT[[st1,st2]]+1)))/((numberOfStates+Apply[Plus,tCountMatricesUT[[st1]]])^2*(1+numberOfStates+Apply[Plus,tCountMatricesUT[[st1]]]));
Plot[getLog10GaussNorm[meanDirichletUT,Sqrt[varDirichletUT]]*10^x*Exp[-(10^x-meanDirichletUT)^2/(2 varDirichletUT)]*plotDetailedScale,{x,plotDetailedMin,plotDetailedMax},PlotRange->All,PlotStyle-> {Thick,RGBColor[1,0,0]}]
];
meanRTDetailedPlotGauss[st1_,st2_]:=Module[{meanDirichletRT,varDirichletRT,maxDirichletRT},
meanDirichletRT=(tCountMatrixRT6by6[[st1,st2]]+1)/(numberOfStates+Apply[Plus,tCountMatrixRT6by6[[st1]]]);
maxDirichletRT=(tCountMatrixRT6by6[[st1,st2]])/(Apply[Plus,tCountMatrixRT6by6[[st1]]]);
varDirichletRT=((tCountMatrixRT6by6[[st1,st2]]+1)(numberOfStates+Apply[Plus,tCountMatrixRT6by6[[st1]]]-(tCountMatrixRT6by6[[st1,st2]]+1)))/((numberOfStates+Apply[Plus,tCountMatrixRT6by6[[st1]]])^2*(1+numberOfStates+Apply[Plus,tCountMatrixRT6by6[[st1]]]));
Plot[getLog10GaussNorm[meanDirichletRT,varDirichletRT]*10^x*Exp[-(10^x-meanDirichletRT)^2/(2 varDirichletRT)]*plotDetailedScale,{x,plotDetailedMin,plotDetailedMax},PlotRange->All,Filling->0,FillingStyle->RGBColor[0,0,0,0.1],PlotStyle-> {Thick,RGBColor[0,0,0]}]
];
(* Plot for Reference Temperature Mean Dirichlet (BLACK - DASHED) *)
meanRTDetailedPlotDirichlet[st1_,st2_]:=Module[{normalizeDirichletRT},
normalizeDirichletRT=getLog10DirNorm[tCountMatrixRT6by6[[st1,st2]],Apply[Plus,tCountMatrixRT6by6[[st1]]]-tCountMatrixRT6by6[[st1,st2]]];
Plot[plotDetailedScale*normalizeDirichletRT*10^x*((10^x)^tCountMatrixRT6by6[[st1,st2]])*((1-10^x)^(Apply[Plus,tCountMatrixRT6by6[[st1]]]-tCountMatrixRT6by6[[st1,st2]])),{x,plotDetailedMin,plotDetailedMax},PlotRange->All,Filling->0,FillingStyle->RGBColor[0,0,0,0.1],PlotStyle-> {Thick,Dashed,RGBColor[0,0,0]}]
];

ticks1to0ConvergencePlot=Map[If[Abs[(#+0.1)*10/3-Round[(0.1+#)*10/3]]<0.001,{#,PaddedForm[#,{1,1}],{0.00625,0}},{#,"",{0.00375,0}}]&,Range[0.05,0.99,0.05]];

ConvergenceDevelopmentPlot[data_,temp_,st1_,st2_,col_]:=ListLogLinearPlot[{data[[temp,st1,st2,All,1]],data[[temp,st1,st2,All,2]]},PlotRange->{{100,10000},{0,1}},Joined->True,Filling->True,FillingStyle->{Directive[col,Opacity[0.3]]},FrameTicks->{{{100,""},{150,"(10^2)",0},{1000,"(10^3)"},{7000,"(10^4)",0},{10000,""}},ticks1to0ConvergencePlot,None,None},Frame->True,PlotStyle->Directive[Thickness[Medium],col],BaseStyle->{Small,FontFamily->"Palatino"}]

ConvergenceDevelopmentSTPlot[temp_,st1_,st2_,col_]:=ConvergenceDevelopmentPlot[confidence95STCheck,temp,st1,st2,col];
ConvergenceDevelopmentHTPlot[temp_,st1_,st2_,col_]:=ConvergenceDevelopmentPlot[confidence95HTCheck,temp,st1,st2,col];

Single Temperature Plot for the development of the Transition Probability with Temperature

st1=3;
st2=6;
TransitionProbabilityPlot[st1_,st2_]:=Module[{},
meanSTPlot=MeanSTPlot[st1,st2];
confidence95STPlot=Confidence95STPlot[st1,st2];
meanMTPlot=MeanMTPlot[st1,st2];
confidence95MTPlot=Confidence95MTPlot[st1,st2];
meanHTPlot=MeanHTPlot[st1,st2];
confidence95HTPlot=Confidence95HTPlot[st1,st2];
maxLikeSTPlot=MaxLikeSTPlot[st1,st2];
maxLikeRTPlot=MaxLikeRTPlotCross[st1,st2];
Show[{confidence95HTPlot,meanHTPlot,confidence95STPlot,meanSTPlot,confidence95MTPlot,meanMTPlot,maxLikeRTPlot},PlotRange->{0,1},Frame->True,BaseStyle->{Medium,FontFamily->"Palatino"},FrameLabel->{"Temperature [K]","Transition Probability "Subscript["T",ToString[st1]<>" -> "<>ToString[st2]]},Axes->False]
]

TransitionProbabilityPlot[3,1]

Table[
Export[imgPath<>"plot_Probability_Vs_Temperature/plot_Probility_Vs_Temperature_"<>ToString[st1]<>"to"<>ToString[st2]<>"_wHummer.pdf",TransitionProbabilityPlot[st1,st2]];
,{st1,1,6},{st2,1,6}];

Combined Plot of the Development of All Transition Probabilities with Temperature in a Grid

(RT, DR, ST, RT) with Errorbars

plot1=Grid[Table[
meanSTPlot=MeanSTPlot[st1,st2];
confidence95STPlot=Confidence95STPlot[st1,st2];
meanMTPlot=MeanMTPlot[st1,st2];
confidence95MTPlot=Confidence95MTPlot[st1,st2];
meanHTPlot=MeanHTPlot[st1,st2];
confidence95HTPlot=Confidence95HTPlot[st1,st2];
maxLikeSTPlot=MaxLikeSTPlot[st1,st2];
confidence95HTPlot=Confidence95HTPlot[st1,st2];
maxLikeRTPlot=MaxLikeRTPlotCross[st1,st2];
Show[{confidence95HTPlot,meanHTPlot,confidence95STPlot,meanSTPlot,confidence95MTPlot,meanMTPlot,maxLikeRTPlot},ImageSize->{If[st2==1,40,0],If[st1==6,20,0]}+{200,130},Frame->True,FrameTicks->{Automatic,Automatic,None,None},BaseStyle->{FontSize->Scaled[0.060],FontFamily->"Palatino"},Axes->False,ImagePadding->{{If[st2==1,40,0],0},{If[st1==6,20,0],0}}+1,PlotRangeClipping->True,PlotRangePadding->Scaled[0.04],PlotRange->{{273,600},{0.0,1.0}},Frame->True],{st1,1,6},{st2,1,6}],ItemSize->Full,Spacings->{0,-0.5}];

plot1

Export[imgPath<>"plot_allTransitionProbabilites_wHummer.pdf",plot1]

(RT, DR, ST) with Errorbars

ticks1to0=Map[If[Abs[(#+0.2)*5-Round[(0.2+#)*5]]<0.001,{#,PaddedForm[#,{1,1}],{0.00625,0}},{#,"",{0.00375,0}}]&,Range[0.05,0.99,0.05]];

ticks1to0=Join[Map[If[False,{#,PaddedForm[#,{1,1}],{0.00625,0}},{#,"",{0.00375,0}}]&,Range[0.05,0.99,0.05]],{{0.5,"0.5",{0.00625,0}},{0.8,"0.8",{0.00625,0}},{0.2,"0.2",{0.00625,0}}}];

ticksTemp=Map[If[Abs[(#-50)/125-Round[(#-50)/125]]<0.001,{#,ToString[#],{0.00625,0}},{#,"",{0.00375,0}}]&,Range[275,600,25]];

Plot[Sin[x/100]^2,{x,273,600},Frame->True,FrameTicks->{ticksTemp,ticks1to0,None,None},PlotRangePadding->Scaled[0.01]]

$niceFontSize:=24;

MaxLikeRTPlotCross[st1_,st2_]:=Module[{},
Graphics[{Black,AbsoluteThickness[0.75],
Line[{Offset[{-2,10},{302,meanRT[[st1,st2]]}],Offset[{2,10},{302,meanRT[[st1,st2]]}]}],
Line[{Offset[{-2,-10},{302,meanRT[[st1,st2]]}],Offset[{2,-10},{302,meanRT[[st1,st2]]}]}],
Line[{Offset[{-10,-2},{302,meanRT[[st1,st2]]}],Offset[{-10,2},{302,meanRT[[st1,st2]]}]}],
Line[{Offset[{10,-2},{302,meanRT[[st1,st2]]}],Offset[{10,2},{302,meanRT[[st1,st2]]}]}],
Disk[{302,meanRT[[st1,st2]]},Scaled[{0.012,0.012}/{1,130/201}]],
Line[{Offset[{0,-10},{302,meanRT[[st1,st2]]}],Offset[{0,10},{302,meanRT[[st1,st2]]}]}],
Line[{Offset[{-10,0},{302,meanRT[[st1,st2]]}],Offset[{10,0},{302,meanRT[[st1,st2]]}]}]
}]
];

picTable=Table[
meanSTPlot=MeanSTPlot[st1,st2];
confidence95STPlot=Confidence95STPlot[st1,st2];
meanMTPlot=MeanMTPlot[st1,st2];
confidence95MTPlot=Confidence95MTPlot[st1,st2];
meanHTPlot=MeanHTPlot[st1,st2];
confidence95HTPlot=Confidence95HTPlot[st1,st2];
maxLikeRTPlot=MaxLikeRTPlotCross[st1,st2];
maxLikeSTPlot=MaxLikeSTPlot[st1,st2];
numberingPlot=Graphics[Style[Text[ToString[st1]<>"->"<>ToString[st2],{560,0.9},{0,0} ],FontSize->$niceFontSize]];
Show[{confidence95STPlot,meanSTPlot,confidence95MTPlot,meanMTPlot,maxLikeRTPlot,numberingPlot},Frame->True,FrameTicks->{ticksTemp,ticks1to0,None,None},ImageSize->{If[st2==1,35,0],If[st1==6,25,0]}+{200,130},BaseStyle->{FontSize->$niceFontSize,FontFamily->"Palatino"},Axes->False,AspectRatio->130/200.5,ImagePadding->{{If[st2==1,35,0],0},{If[st1==6,25,0],0}}+0.1,PlotRangeClipping->True,PlotRangePadding->Scaled[0.02],PlotRange->{{273,600},{0.0,1.0}},Frame->True],{st1,1,6},{st2,1,6}];

textLeftRaw="Transition Probability (Subscript[T, i -> j])";
textLeft=Graphics[Rotate[Style[Text[textLeftRaw,{0.5,0.5}],FontFamily->"Palatino",FontSize->$niceFontSize],90Degree],ImageSize->{35,300}];
textBottomRaw="Temperature [K]";
textBottom=Graphics[Style[Text[textBottomRaw,{0.5,0.5}],FontFamily->"Palatino",FontSize->$niceFontSize],ImageSize->{300,35}];
plot1=Grid[Transpose[Join[{{textLeft,SpanFromAbove,SpanFromAbove,SpanFromAbove,SpanFromAbove,SpanFromAbove,SpanFromAbove}},Transpose[Join[picTable,{{textBottom,SpanFromLeft,SpanFromLeft,SpanFromLeft,SpanFromLeft,SpanFromLeft}}]]]],ItemSize->Full,Spacings->{0.05,0.05},Alignment->{Center,Center}];

plot1

Export[imgPath<>"plot_allTransitionProbabilites_woLimits.pdf",plot1]

(RT, DR, ST, RT) only Mean

plot1=Grid[Table[
meanSTPlot=MeanSTPlot[st1,st2];
confidence95STPlot=Confidence95STPlot[st1,st2];
meanMTPlot=MeanMTPlot[st1,st2];
confidence95MTPlot=Confidence95MTPlot[st1,st2];
meanHTPlot=MeanHTPlot[st1,st2];
confidence95HTPlot=Confidence95HTPlot[st1,st2];
maxLikeSTPlot=MaxLikeSTPlot[st1,st2];
confidence95HTPlot=Confidence95HTPlot[st1,st2];
maxLikeRTPlot=MaxLikeRTPlotCross[st1,st2];
Show[{meanHTPlot,meanSTPlot,meanMTPlot,maxLikeRTPlot},ImageSize->{If[st2==1,40,0],If[st1==6,20,0]}+{200,130},Frame->True,FrameTicks->{Automatic,Automatic,None,None},BaseStyle->{FontSize->Scaled[0.060],FontFamily->"Palatino"},Axes->False,ImagePadding->{{If[st2==1,40,0],0},{If[st1==6,20,0],0}}+1,PlotRangeClipping->True,PlotRangePadding->Scaled[0.04],PlotRange->{{273,600},{0.0,1.0}},Frame->True],{st1,1,6},{st2,1,6}],ItemSize->Full,Spacings->{0,-0.5}];

plot1

Export[imgPath<>"plot_allTransitionProbabilitesMean_wHummer.pdf",plot1]

Old Plot with seperate Frame for each Plot

plot1=GraphicsArray[Table[
meanSTPlot=MeanSTPlot[st1,st2];
confidence95STPlot=Confidence95STPlot[st1,st2];
meanMTPlot=MeanMTPlot[st1,st2];
confidence95MTPlot=Confidence95MTPlot[st1,st2];
meanHTPlot=MeanHTPlot[st1,st2];
confidence95HTPlot=Confidence95HTPlot[st1,st2];
maxLikeSTPlot=MaxLikeSTPlot[st1,st2];
confidence95HTPlot=Confidence95HTPlot[st1,st2];
maxLikeRTPlot=MaxLikeRTPlotCross[st1,st2];
Show[{confidence95STPlot,meanSTPlot,confidence95HTPlot,meanHTPlot,confidence95MTPlot,meanMTPlot,maxLikeRTPlot},Frame->True,BaseStyle->{Small,FontFamily->"Palatino"},FrameLabel->{"Temperature [K]","Probability "Subscript["T",ToString[st1]<>" -> "<>ToString[st2]]},AspectRatio->0.65,Axes->False,PlotRange->{{273,600},{0.0,1.0}},Frame->True,BaseStyle->{Small,FontFamily->"Palatino",Italic}],{st1,1,6},{st2,1,6}],ImageSize->1450,GraphicsSpacing->-0.01];

plot1

Export[imgPath<>"plot_Probility_Vs_Temperature_All_Transitions.pdf",plot1];

Plots for Detailed Comparison at 302 K

Definition of PlotFuntions

plotDetailedScale=1;
plotDetailedMin=0.000;
plotDetailedMax=1.000;
plotDetailedMinLog=Log[plotDetailedMin];
plotDetailedMaxLog=Log[plotDetailedMax=1.000];

meanSTDetailedPlotHistogram[st1_,st2_]:=Graphics[{RGBColor[1.0,0,0],EdgeForm[{Thin,GrayLevel[0]}],PlotHistogram[Log[tMatricesST6by6[[6,All,st1,st2]]],{plotDetailedMinLog,plotDetailedMaxLog},120]}];
meanHTDetailedPlotHistogram[st1_,st2_]:=Graphics[{RGBColor[0,0.7,0,.8],EdgeForm[{Thin,GrayLevel[0]}],PlotHistogram[Log[tMatricesHT[[6,All,st1,st2]]],{plotDetailedMinLog,plotDetailedMaxLog},60]}];
meanRTDetailedPlotHistogram[st1_,st2_]:=Graphics[{RGBColor[0.3,0.3,0.3,.2],EdgeForm[{Thin,GrayLevel[0]}],PlotHistogram[Log[tMatricesRT6by6[[All,st1,st2]]],{plotDetailedMinLog,plotDetailedMaxLog},270]}];

meanDetailedPlotHistogram[data_,bins_,color_]:=Graphics[{color,EdgeForm[{Thin,GrayLevel[0]}],PlotHistogram[data,{0,1},bins]}];
meanDetailedPlotGauss[mean_,var_,color_]:=Plot[PDF[NormalDistribution[mean,var]][x],{x,plotDetailedMin,plotDetailedMax},PlotRange->All,PlotStyle-> {Thick,color},Filling->Axis,FillingStyle->Directive[color,Opacity[0.2]],Frame->True,FrameTicks->{Automatic,Automatic,Automatic,Automatic}]

meanDetailedPlotHistogramSmooth[data_,bins_,color_]:=Graphics[{color,PlotHistogramOutlineSmooth[data,{0,1},bins]}];
meanDetailedPlotGauss[mean_,var_,color_]:=Plot[PDF[NormalDistribution[mean,var]][x],{x,plotDetailedMin,plotDetailedMax},PlotRange->All,Filling->Axis,FillingStyle->Directive[color,Opacity[0.2]],PlotStyle-> {Thickness[Medium],color},Frame->True,FrameTicks->{Automatic,Automatic,Automatic,Automatic}]

meanDetailedPlotHistogramSmooth[data_,bins_,color_]:=ListPlot[GetHistogramOutlineSmooth[data,{0,1},bins],Joined->True,PlotRange->All,Filling->Axis,FillingStyle->Directive[color,Opacity[0.2]],PlotStyle-> Directive[Thickness[Medium],color],Frame->True,FrameTicks->{Automatic,Automatic,Automatic,Automatic}]

meanDetailedPlotHistogramLog[data_,bins_,color_]:=Graphics[{color,EdgeForm[{Thin,GrayLevel[0]}],PlotHistogram[Log[data],{plotDetailedMinLog,plotDetailedMaxLog},bins]}];
meanDetailedPlotGaussLog[mean_,var_,color_]:=LogLinearPlot[PDF[NormalDistribution[mean,var]][x]*x,{x,plotDetailedMin,plotDetailedMax},PlotRange->All,PlotStyle-> {Thick,color},Filling->Axis,FillingStyle->Directive[color,Opacity[0.2]],Frame->True,FrameTicks->{Automatic,Automatic,Automatic,Automatic}]

Single Plots

TransitionProbability302KPlot[st1_,st2_]:=Module[{},
(* Plot for MT Estimates (BLUE) *)
plotMTGauss=meanDetailedPlotGauss[meanMT[[11,st1,st2]],confidence68MT[[11,st1,st2]],RGBColor[0,0,1]];

(* Histogram for TM Sampling *)
plotSTHistogram=meanDetailedPlotHistogramSmooth[tMatricesST6by6[[6,All,st1,st2]],50,Directive[Dashed,Thick,RGBColor[1,0,0]]];

(* Plot for Single Temperature Mean Gauss (RED) *)
plotSTGauss=meanDetailedPlotGauss[meanST[[6,st1,st2]],varianceST[[6,st1,st2]],Directive[Dashed,RGBColor[1,0,0]]];

(* Plot for Reference Temperature Mean Gauss (BLACK - DASHED) *)
plotRTGauss=meanDetailedPlotGauss[meanRT[[st1,st2]],varianceRT[[st1,st2]],RGBColor[0,0,0]];
(*plotRTGauss=Graphics[{AbsoluteThickness[5],Black,Line[{{Log[meanRT[[st1,st2]]],0},{Log[meanRT[[st1,st2]]],100}}]}];*)

(* Plot for Hummer Temperature Mean Gauss (GREEN) *)
plotHTGauss=meanDetailedPlotGauss[meanHT[[6,st1,st2]],varianceHT[[6,st1,st2]],Directive[Thick,RGBColor[0,0.7,0]]];

(* Histogram for Reference Temperature *)
plotRTHistogram=meanDetailedPlotHistogramSmooth[tMatricesRT[[6,All,st1,st2]],180,Directive[Thick,RGBColor[1,0,0,1]]];

(* Histogram for Reference Temperature from Hummer Rate Estimates*)
plotHTHistogram=meanDetailedPlotHistogramSmooth[tMatricesHT6by6[[6,All,st1,st2]],90,Directive[Dotted,Thick,RGBColor[0,0.7,0]]];

transitionLabel=Subscript["T",ToString[st1]<>ToString[st2]];
labels={Row[{"Transition Probability ",transitionLabel}],Row[{"p(",transitionLabel,")"}]};

(* Combine Plots *)
Show[{plotMTGauss,plotHTHistogram,plotSTHistogram,plotMTGauss,plotRTGauss},Frame->True,BaseStyle->{FontSize->Scaled[0.040],FontFamily->"Palatino"},Axes->False,Frame->True,ImageSize->520,AspectRatio->1/GoldenRatio,FrameLabel->labels,ImagePadding->{{60,10},{54,4}},PlotRange->{0,40.2}]
];

TransitionProbability302KPlot[3,2]

Table[
Export[imgPath<>"plot_Detailed_Probability_"<>ToString[st1]<>"to"<>ToString[st2]<>"_Single.pdf",TransitionProbility302KPlot[st1,st2]];
,{st1,6,6},{st2,3,3}];

plotX

Plot of Full Table

ticks1to0=Map[If[Abs[(#+0.1)*5-Round[(0.1+#)*5]]<0.001,{#,PaddedForm[#,{1,1}],{0.00625,0}},{#,"",{0.00375,0}}]&,Range[0,1,0.05]]

ticks0to15=Map[If[Abs[(#-1)/4-Round[(#-1)/4]]<0.001,{#,PaddedForm[#,{1,1}],{0.00625,0}},{#,"",{0.00375,0}}]&,Range[1,9,1]];

picTable=Table[
plotMTGauss=meanDetailedPlotGauss[meanMT[[11,st1,st2]],confidence68MT[[11,st1,st2]],RGBColor[0,0,1]];
(* Histogram for TM Sampling *)
plotSTHistogram=meanDetailedPlotHistogramSmooth[tMatricesST6by6[[6,All,st1,st2]],50,Directive[Thickness[Medium],RGBColor[1,0,0]]];
(* Plot for Reference Temperature Mean Gauss (BLACK - DASHED) *)
plotRTGauss=meanDetailedPlotGauss[meanRT[[st1,st2]],varianceRT[[st1,st2]],RGBColor[0,0,0]];
(* Histogram for Reference Temperature *)
plotHTHistogram=meanDetailedPlotHistogramSmooth[tMatricesHT6by6[[6,All,st1,st2]],50,Directive[Thickness[Medium],RGBColor[0,0.7,0]]];
numberingPlot=Graphics[Style[Text[ToString[st1]<>"->"<>ToString[st2],Scaled[{0.97,0.97}],{1,1} ],FontSize->$niceFontSize]];
Show[{plotMTGauss,plotHTHistogram,plotSTHistogram,plotRTGauss,numberingPlot},Frame->True,FrameTicks->{ticks1to0,ticks0to15,None,None},ImageSize->{If[st2==1,35,0],If[st1==6,25,0]}+{200,130},BaseStyle->{FontSize->$niceFontSize,FontFamily->"Palatino"},Axes->False,AspectRatio->130/200.5,ImagePadding->{{If[st2==1,35,0],0},{If[st1==6,25,0],0}}+0.1,PlotRangeClipping->True,PlotRangePadding->Scaled[0.02],PlotRange->{{0.0,1.0},{0.0,10.0}},Frame->True],{st1,1,6},{st2,1,6}];

textLeftRaw="Probability p((Subscript[T, i -> j]))";
textLeft=Graphics[Rotate[Style[Text[textLeftRaw,{0.5,0.5}],FontFamily->"Palatino",FontSize->$niceFontSize],90Degree],ImageSize->{35,300}];
textBottomRaw="Transition probability (Subscript[T, i -> j])";
textBottom=Graphics[Style[Text[textBottomRaw,{0.5,0.5}],FontFamily->"Palatino",FontSize->$niceFontSize],ImageSize->{300,35}];
plot1=Grid[Transpose[Join[{{textLeft,SpanFromAbove,SpanFromAbove,SpanFromAbove,SpanFromAbove,SpanFromAbove,SpanFromAbove}},Transpose[Join[picTable,{{textBottom,SpanFromLeft,SpanFromLeft,SpanFromLeft,SpanFromLeft,SpanFromLeft}}]]]],ItemSize->Full,Spacings->{0.1,0.1},Alignment->{Center,Center}];

plot1

Export[imgPath<>"plot_allComparisonDetailed302K_woLimits.pdf",plot1]

Without Filling

plot1=Grid[Table[
plotMTGauss=meanDetailedPlotGauss[meanMT[[11,st1,st2]],confidence68MT[[11,st1,st2]],RGBColor[0,0,1]];
(* Histogram for TM Sampling *)
plotSTHistogram=meanDetailedPlotHistogramSmooth[tMatricesST6by6[[6,All,st1,st2]],50,Directive[Dashed,Thick,RGBColor[1,0,0]]];
(* Plot for Reference Temperature Mean Gauss (BLACK - DASHED) *)
plotRTGauss=meanDetailedPlotGauss[meanRT[[st1,st2]],varianceRT[[st1,st2]],RGBColor[0,0,0]];
(* Histogram for Reference Temperature *)
plotHTHistogram=meanDetailedPlotHistogramSmooth[tMatricesHT6by6[[6,All,st1,st2]],50,Directive[Dotted,Thick,RGBColor[0,0.7,0]]];

transitionLabel=Subscript["T",ToString[st1]<>ToString[st2]];
labels={Row[{"Transition Probability ",transitionLabel}],Row[{"p(",transitionLabel,")"}]};

(* Combine Plots *)
Show[{plotMTGauss,plotHTHistogram,plotSTHistogram,plotRTGauss},ImageSize->{If[st2==1,40,0],If[st1==6,20,0]}+{200,130}/1.4,Frame->True,FrameTicks->{Automatic,Automatic,None,None},BaseStyle->{FontSize->Scaled[0.060],FontFamily->"Palatino"},Axes->False,ImagePadding->{{If[st2==1,40,0],0},{If[st1==6,20,0],0}}+1,PlotRangeClipping->True,PlotRangePadding->Scaled[0.04],PlotRange->{{-.01,1.01},{0.0,15.0}},Frame->True],{st1,1,6},{st2,1,6}],ItemSize->Full,Spacings->{0,-0.5}];

plot1

Export[imgPath<>"plot_allComparisonDetailed302K_woLimits.pdf",plot1]

Comparison with Reference

distPlotST=Table[LogLogPlot[meanRT[[st1,st2]],{x,0.001,1},ColorFunction->(RGBColor[1,0,0,1-#]&[1-Exp[-(#-meanST[[6,st1,st2]])^2/(2 varianceST[[6,st1,st2]]^2)]]&),PlotStyle->Thickness[0.02],PlotRange->{{0.001,1},{0.001,1}},Frame->True,Axes->False],{st1,1,6},{st2,1,6}];
distPlotHT=Table[LogLogPlot[meanRT[[st1,st2]],{x,0.001,1},ColorFunction->(RGBColor[0,0.7,0,1-#]&[1-Exp[-(#-meanHT[[6,st1,st2]])^2/(2 varianceHT[[6,st1,st2]]^2)]]&),PlotStyle->Thickness[0.02],PlotRange->{{0.001,1},{0.001,1}},Frame->True,Axes->False],{st1,1,6},{st2,1,6}];
distPlotMT=Table[LogLogPlot[meanRT[[st1,st2]],{x,0.001,1},ColorFunction->(RGBColor[0,0,1,1-#]&[1-Exp[-(#-meanMT[[6,st1,st2]])^2/(2 confidence68MT[[6,st1,st2]]^2)]]&),PlotStyle->Thickness[0.02],PlotRange->{{0.001,1},{0.001,1}},Frame->True,Axes->False],{st1,1,6},{st2,1,6}];

exactPlot=Plot[{x+0.1,x-0.1},{x,Log[0.001],0},Filling->True,PlotRange->{Log[0.001],0},AspectRatio->Automatic];

Show[distPlotST,distPlotHT,distPlotMT,exactPlot,AspectRatio->Automatic]

ticks={0.001,0.002,0.003,0.004,0.005,0.006,0.007,0.008,0.009,0.01,0.02,0.03,0.04,0.05,0.06,0.07,0.08,0.09,0.1,0.2,0.3,0.4,0.5,0.6,0.7,0.8,0.9,1,2,3,4,5,6,7,8,9,10,20,30,40,50,60,70,80,90,100,200,300,400,500,600,700,800,900,1000};
ticks2=Map[If[Log[10,#]==Round[Log[10,#]],{Log[10,#],ToString[#],{0.01,0}},{Log[10,#],"",{0.005,0}}]&,ticks];

datalist2=Map[getHistogram[#,{-1.05,1.05},20]&,datalist][[All,All,2]]

BarChart[Transpose[datalist2]]

datalist=Map[#/Flatten[meanRT]&,{Flatten[meanST[[6]]],Flatten[meanHT[[6]]],Flatten[meanMT[[11]]]}];
datalist=Map[Log[If[#<1,#,#]]&,datalist,{2}];
plot=Histogram[datalist,{-3,3,1/2},"Probability",ChartStyle->{Directive[Opacity[0.6],Red],Directive[Opacity[.6],Green],Directive[Opacity[.6],Blue]},Frame->True,FrameTicks->{{Automatic,Automatic},{ticks2,Automatic}},FrameLabel->{"Relative Deviation of Transition Probabilities"},BaseStyle->{"Palatino"},FrameTicksStyle->Directive[FontFamily->"Palatino"]]
Export[imgPath<>"plot_DeviationInTransitionProbabilityLog_All.pdf",plot]

datalist=Map[#-Flatten[meanRT]&,{Flatten[meanMT[[11]]],Flatten[meanST[[6]]],Flatten[meanHT[[6]]]}];
datalist=Map[Identity[If[#<1,#,#]]&,datalist,{2}];
plot=Histogram[datalist,{-1.05,1.05,0.1},"Probability",ChartStyle->{Directive[Opacity[1],Blue],Directive[Opacity[.9],Red],Directive[Opacity[0.8],Green]},Frame->True,FrameTicks->{{Automatic,Automatic},{Automatic,Automatic}},FrameLabel->{"Absolute Deviation of Transition Probabilities"},BaseStyle->{"Palatino"},FrameTicksStyle->Directive[FontFamily->"Palatino"],AxesOrigin->{0,0.01},Axes->False,Epilog->{Thick,Line[{{-1,0},{1,0}}]}]
Export[imgPath<>"plot_DeviationInTransitionProbability_All.pdf",plot]

Variances of Deviation

Map[StandardDeviation[#[[Range[25,36]]]]&,datalist]

Map[StandardDeviation[#[[Range[24]]]]&,datalist]

Map[StandardDeviation,datalist]

Tables of Implied Time Scales at 302K

lagTime=6ps;

Single Temperature (TE)

Transpose[{-lagTime/Log[eigenvaluesST[[6]]+eigenvalueVarianceST[[6]]],
-lagTime/Log[eigenvaluesST[[6]]],
-lagTime/Log[eigenvaluesST[[6]]-eigenvalueVarianceST[[6]]],
eigenvalueVarianceST[[6]]/eigenvaluesST[[6]]}][[Range[2,6]]]//TableForm

Reference Temperature (RT)

Transpose[{-lagTime/Log[eigenvaluesRT-eigenvalueVarianceRT],
-lagTime/Log[eigenvaluesRT],
-lagTime/Log[eigenvaluesRT+eigenvalueVarianceRT],
eigenvalueVarianceRT/eigenvaluesRT}][[Range[2,6]]]//TableForm

Dynamical Reweighting (DR)

Transpose[{-lagTime/Log[eigenvaluesMT[[11]]+eigenvalueVarianceMT[[11]]],
-lagTime/Log[eigenvaluesMT[[11]]],
-lagTime/Log[eigenvaluesMT[[11]]-eigenvalueVarianceMT[[11]]],
eigenvalueVarianceMT[[11]]/eigenvaluesMT[[11]]}][[Range[2,6]]]//TableForm

Single Temperature Rate Matrix Estimation (RE)

Transpose[{-lagTime/Log[eigenvaluesHT[[6]]+eigenvalueVarianceHT[[6]]],
-lagTime/Log[eigenvaluesHT[[6]]],
-lagTime/Log[eigenvaluesHT[[6]]-eigenvalueVarianceHT[[6]]],
eigenvalueVarianceHT[[6]]/eigenvaluesHT[[6]]}][[Range[2,6]]]//TableForm

Eigenvalue Plot for Single Eigenvalues

Color Definitions

colorRange={{-1,Black},{-0.5,Blue},{0,White},{.5,Yellow},{1,Red}};

Compute Eigenvectors and Eigenvalues at 302K for comparison

symMatrixST=DiagonalMatrix[Sqrt[Equilibrium[TransitionMatrix[meanST[[6]]]]]];
symMatrixHT=DiagonalMatrix[Sqrt[Equilibrium[TransitionMatrix[meanHT[[6]]]]]];
symMatrixMT=DiagonalMatrix[Sqrt[Equilibrium[TransitionMatrix[meanMT[[11]]]]]];
symMatrixRT=DiagonalMatrix[Sqrt[Equilibrium[TransitionMatrix[meanRT]]]];
tEigenvectorsST=Eigenvectors[meanST[[6]]];
tEigenvectorsHT=Eigenvectors[meanHT[[6]]];
tEigenvectorsMT=Eigenvectors[meanMT[[11]]];
tEigenvectorsRT=Eigenvectors[meanRT];
tEigenvectorsSymST=Eigenvectors[symMatrixST.meanST[[6]].Inverse[symMatrixST]];
tEigenvectorsSymHT=Eigenvectors[symMatrixHT.meanHT[[6]].Inverse[symMatrixHT]];
tEigenvectorsSymMT=Eigenvectors[symMatrixMT.meanMT[[11]].Inverse[symMatrixMT]];
tEigenvectorsSymRT=Eigenvectors[symMatrixRT.meanRT.Inverse[symMatrixRT]];
tEigenvectorsLeftST=Transpose[Inverse[Eigenvectors[meanST[[6]]]]];
tEigenvectorsLeftHT=Transpose[Inverse[Eigenvectors[meanHT[[6]]]]];
tEigenvectorsLeftMT=Transpose[Inverse[Eigenvectors[meanMT[[11]]]]];
tEigenvectorsLeftRT=Transpose[Inverse[Eigenvectors[meanRT]]];

convN=IdentityMatrix[6];
convST=Chop[Inverse[convN].Abs[Outer[Dot,tEigenvectorsSymST,tEigenvectorsSymRT,1]]];
convHT=Chop[Inverse[convN].Abs[Outer[Dot,tEigenvectorsSymHT,tEigenvectorsSymRT,1]]];
convMT=Chop[Inverse[convN].Abs[Outer[Dot,tEigenvectorsSymMT,tEigenvectorsSymRT,1]]];
convRT=Chop[Inverse[convN].Chop[Abs[Outer[Dot,tEigenvectorsSymRT,tEigenvectorsSymRT,1]]]];

Eigenvalues vs Temperature

EigenvalueRTPlotCross[ev_]:=Module[{size=6},
Graphics[{Black,AbsoluteThickness[1.0],Circle[{302,eigenvaluesRT[[ev]]},Offset[{size,size}]],
Line[{Offset[{0,-size},{302,eigenvaluesRT[[ev]]}],Offset[{0,size},{302,eigenvaluesRT[[ev]]}]}],
Line[{Offset[{-size,0},{302,eigenvaluesRT[[ev]]}],Offset[{size,0},{302,eigenvaluesRT[[ev]]}]}]
}]
];

MaxLikeRTPlotCross[st1_,st2_]:=Module[{},
Graphics[{Black,AbsoluteThickness[1.0],
Line[{Offset[{-2,10},{302,meanRT[[st1,st2]]}],Offset[{2,10},{302,meanRT[[st1,st2]]}]}],
Line[{Offset[{-2,-10},{302,meanRT[[st1,st2]]}],Offset[{2,-10},{302,meanRT[[st1,st2]]}]}],
Disk[{302,meanRT[[st1,st2]]},Scaled[{0.012,0.012}/{1,130/201}]],
Line[{Offset[{0,-10},{302,meanRT[[st1,st2]]}],Offset[{0,10},{302,meanRT[[st1,st2]]}]}]
}]
];

evsST={2};
evsMT={2};
evsRT={2,3};
evsHT={3};
eigenvalueSTPlot=Map[EigenvalueSTPlot,evsST];
eigenvalueMTPlot=Map[EigenvalueMTPlot,evsMT];
eigenvalueRTPlot=Map[EigenvalueRTPlot,evsRT];
eigenvalueHTPlot=Map[EigenvalueHTPlot,evsHT];
eigenvalueRTPlotCross=Map[EigenvalueRTPlotCross,evsRT];
eigenvalueSTVariancePlotList =Map[EigenvalueSTVariancePlot,evsST];
eigenvalueMTVariancePlotList = Map[EigenvalueMTVariancePlot,evsMT];
eigenvalueHTVariancePlotList =Map[EigenvalueHTVariancePlot,evsHT];
plotInlet=Show[{eigenvalueSTPlot,eigenvalueMTPlot,eigenvalueRTPlot,eigenvalueMTVariancePlotList,eigenvalueSTVariancePlotList,eigenvalueRTPlotCross},Frame->True,AspectRatio->0.8,Axes->False,PlotRange->{{280,320},{0.72,0.88}},FrameTicks->If[evsST[[1]]==2,{{},{},{285,300,315},{.75,{.8,"0.80"},.85}},{{285,300,315},{.75,{.8,"0.80"},.85},{},{}}],Frame->True,BaseStyle->{FontSize->8,FontFamily->"Palatino"},ImageSize->120]
plot1=Show[{eigenvalueSTVariancePlotList,eigenvalueSTPlot,eigenvalueMTVariancePlotList,eigenvalueMTPlot,eigenvalueRTPlot,eigenvalueRTPlotCross},Frame->True,FrameLabel->{"Temperature [K]",Row[{"Eigenvalue ",Subscript["λ",Style[evsST[[1]],FontSize->10]]}]},AspectRatio->0.65,Axes->False,PlotRange->{{273,600},{0.0,1.0}},Frame->True,BaseStyle->{FontSize->12,FontFamily->"Palatino"},Epilog->{Inset[plotInlet,If[evsST[[1]]==2,{278,0.02},{595,0.97}],If[evsST[[1]]==2,{Left,Bottom},{Right,Top}],130],Black,Line[{{280,0.72},{280,0.88},{320,0.88},{320,0.72},{280,0.72}}]},ImageSize->240]

Export[imgPath<>"plot_Eigenvalue_Vs_Temperature_No"<>ToString[evsST[[1]]]<>".pdf",plot1]

Export[imgPath<>"plot_Eigenvalue_Vs_Temperature_No"<>ToString[evsST[[1]]]<>"vs"<>ToString[evsMT[[1]]]<>".pdf",plot1]

evsST={2};
evsMT={3};
evsRT={3};
evsHT={3};
eigenvalueSTPlot=Map[EigenvalueSTPlot,evsST];
eigenvalueMTPlot=Map[EigenvalueMTPlot,evsMT];
eigenvalueRTPlot=Map[EigenvalueRTPlot,evsRT];
eigenvalueHTPlot=Map[EigenvalueHTPlot2,evsHT];
eigenvalueRTPlotCross=Map[EigenvalueRTPlotCross,evsRT];
eigenvalueSTVariancePlotList =Map[EigenvalueSTVariancePlot,evsST];
eigenvalueMTVariancePlotList = Map[EigenvalueMTVariancePlot,evsMT];
eigenvalueHTVariancePlotList =Map[EigenvalueHTVariancePlot2,evsHT];
plot1=Show[{eigenvalueSTPlot,eigenvalueMTPlot,eigenvalueRTPlot,eigenvalueMTVariancePlotList,eigenvalueSTVariancePlotList,eigenvalueRTPlotCross},Frame->True,FrameLabel->{"Temperature [K]","Eigenvalue"Subscript["λ",Style[evs[[1]],FontSize->7]]},AspectRatio->0.65,Axes->False,PlotRange->{{273,600},{0.0,1.0}},Frame->True,BaseStyle->{FontSize->Small,FontFamily->"Palatino"},ImageSize->240]

Export[imgPath<>"plot_Eigenvalue_Vs_Temperature_No3vs2.pdf",plot1]

Comparison of eigenvectors

$niceFontSize=12;
tickList=Table[{nn+0.5,nn},{nn,1,6}];
legendPlot=Show[DensityPlot[y,{x,0,0.1},{y,0,1},Frame->True,AspectRatio->10,FrameTicks->{None,None,None,Map[{#,ToString[Round[100*#]]<>"%"}&,Range[0,1,0.2]]},ColorFunction->(Blend[{{0,White},{1/2,Yellow},{1,Red}},#]&),ImageSize->47,PlotRangePadding->0.02,ImagePadding->{{2,32},{26,30}}],Frame->True,BaseStyle->{"Palatino",$niceFontSize}]
plot1=Grid[{{Graphics[Flatten[Table[{Blend[{{0,White},{1/2,Yellow},{1,Red},{1.5,Blue}},#]&[convRT[[ii,jj]]],Rectangle[{ii,jj},{ii,jj}+1]},{ii,1,6},{jj,1,6}]],FrameLabel->{"Reference\n(ST) EV","Reference EV"},Frame->True,FrameTicks->{tickList,tickList,None,None},PlotRange->{{1,7},{1,7}},PlotRangePadding->Scaled[0.04],ImagePadding->{{40,1},{60,1}},ImageSize->{200,220},Epilog->{Thick,Line[{{1,2},{7,2},{7,3},{1,3},{1,2}}]},BaseStyle->{"Palatino",14}],
Graphics[Flatten[Table[{Blend[{{0,White},{1/2,Yellow},{1,Red},{1.5,Blue}},#]&[convMT[[ii,jj]]],Rectangle[{ii,jj},{ii,jj}+1]},{ii,1,6},{jj,1,6}]],FrameLabel->{"Dynamical Reweighting\n(DR) EV","Reference EV"},Frame->True,FrameTicks->{tickList,tickList,None,None},PlotRange->{{1,7},{1,7}},PlotRangePadding->Scaled[0.04],ImagePadding->{{1,1},{60,1}},ImageSize->{161,220},Epilog->{Thick,Line[{{1,2},{7,2},{7,3},{1,3},{1,2}}]},BaseStyle->{"Palatino",14}],
Graphics[Flatten[Table[{Blend[{{0,White},{1/2,Yellow},{1,Red},{1.5,Blue}},#]&[convST[[ii,jj]]],Rectangle[{ii,jj},{ii,jj}+1]},{ii,1,6},{jj,1,6}]],FrameLabel->{"Transition Estimation\n(TE) EV","Reference EV"},Frame->True,FrameTicks->{tickList,tickList,None,None},PlotRange->{{1,7},{1,7}},PlotRangePadding->Scaled[0.04],ImagePadding->{{1,1},{60,1}},ImageSize->{161,220},Epilog->{Thick,Line[{{1,2},{7,2},{7,3},{1,3},{1,2}}]},BaseStyle->{"Palatino",14}],
Graphics[Flatten[Table[{Blend[{{0,White},{1/2,Yellow},{1,Red},{1.5,Blue}},#]&[convHT[[ii,jj]]],Rectangle[{ii,jj},{ii,jj}+1]},{ii,1,6},{jj,1,6}]],FrameLabel->{"Rate Estimation\n(RE) EV","Reference EV"},Frame->True,FrameTicks->{tickList,tickList,None,None},PlotRange->{{1,7},{1,7}},PlotRangePadding->Scaled[0.04],ImagePadding->{{1,1},{60,1}},ImageSize->{161,220},Epilog->{Thick,Line[{{1,2},{7,2},{7,3},{1,3},{1,2}}]},BaseStyle->{"Palatino",14}]
,legendPlot}}
,Spacings->0,ItemSize->Full,Alignment->Top,BaseStyle->{"Palatino",16}]

Export[imgPath<>"plot_Eigenvector_Matching_Comparison.pdf",plot1]

Plot of TransitionMatrix Components from left eigenvectors

Table[Graphics[Flatten[Table[{Blend[colorRange,#]&[(tEigenvaluesRT[[kk]]*Outer[Times,tEigenvectorsRT[[kk]],Transpose[Inverse[tEigenvectorsRT]][[kk]]])[[ii,jj]]],Rectangle[{jj,ii},{jj,ii}+1]},{jj,1,6},{ii,1,6}]],Frame->True,FrameTicks->{tickList,tickList,None,None},PlotRange->{{1,7},{1,7}},PlotRangePadding->Scaled[0.04],ImagePadding->{{20,1},{20,1}},ImageSize->{120,120}],
{kk,1,6}]
Table[Graphics[Flatten[Table[{Blend[colorRange,#]&[(tEigenvaluesMT[[kk]]*Outer[Times,tEigenvectorsMT[[kk]],Transpose[Inverse[tEigenvectorsMT]][[kk]]])[[ii,jj]]],Rectangle[{jj,ii},{jj,ii}+1]},{jj,1,6},{ii,1,6}]],Frame->True,FrameTicks->{tickList,tickList,None,None},PlotRange->{{1,7},{1,7}},PlotRangePadding->Scaled[0.04],ImagePadding->{{20,1},{20,1}},ImageSize->{120,120}],
{kk,1,6}]
Table[Graphics[Flatten[Table[{Blend[colorRange,#]&[(tEigenvaluesST[[kk]]*Outer[Times,tEigenvectorsST[[kk]],Transpose[Inverse[tEigenvectorsST]][[kk]]])[[ii,jj]]],Rectangle[{jj,ii},{jj,ii}+1]},{jj,1,6},{ii,1,6}]],Frame->True,FrameTicks->{tickList,tickList,None,None},PlotRange->{{1,7},{1,7}},PlotRangePadding->Scaled[0.04],ImagePadding->{{20,1},{20,1}},ImageSize->{120,120}],
{kk,1,6}]

Plot of TransitionMatrix Components from symmetrized matrix

Table[Graphics[Flatten[Table[{Blend[colorRange,#]&[(tEigenvaluesRT[[kk]]*Outer[Times,symMatrixRT.tEigenvectorsRT[[kk]],Transpose[Inverse[tEigenvectorsRT]][[kk]]].Inverse[symMatrixRT])[[ii,jj]]],Rectangle[{jj,ii},{jj,ii}+1]},{jj,1,6},{ii,1,6}]],Frame->True,FrameTicks->{tickList,tickList,None,None},PlotRange->{{1,7},{1,7}},PlotRangePadding->Scaled[0.04],ImagePadding->{{20,1},{20,1}},ImageSize->{120,120}],
{kk,1,6}]
Table[Graphics[Flatten[Table[{Blend[colorRange,#]&[(tEigenvaluesMT[[kk]]*symMatrixMT.Outer[Times,tEigenvectorsMT[[kk]],Transpose[Inverse[tEigenvectorsMT]][[kk]]].Inverse[symMatrixMT])[[ii,jj]]],Rectangle[{jj,ii},{jj,ii}+1]},{jj,1,6},{ii,1,6}]],Frame->True,FrameTicks->{tickList,tickList,None,None},PlotRange->{{1,7},{1,7}},PlotRangePadding->Scaled[0.04],ImagePadding->{{20,1},{20,1}},ImageSize->{120,120}],
{kk,1,6}]
Table[Graphics[Flatten[Table[{Blend[colorRange,#]&[(tEigenvaluesST[[kk]]*symMatrixST.Outer[Times,tEigenvectorsST[[kk]],Transpose[Inverse[tEigenvectorsST]][[kk]]].Inverse[symMatrixST])[[ii,jj]]],Rectangle[{jj,ii},{jj,ii}+1]},{jj,1,6},{ii,1,6}]],Frame->True,FrameTicks->{tickList,tickList,None,None},PlotRange->{{1,7},{1,7}},PlotRangePadding->Scaled[0.04],ImagePadding->{{20,1},{20,1}},ImageSize->{120,120}],
{kk,1,6}]
Table[Graphics[Flatten[Table[{Blend[colorRange,#]&[(tEigenvaluesST[[kk]]*Outer[Times,symMatrixST.tEigenvectorsST[[kk]],Transpose[Inverse[tEigenvectorsST]][[kk]].Inverse[symMatrixST]])[[ii,jj]]],Rectangle[{jj,ii},{jj,ii}+1]},{jj,1,6},{ii,1,6}]],Frame->True,FrameTicks->{tickList,tickList,None,None},PlotRange->{{1,7},{1,7}},PlotRangePadding->Scaled[0.04],ImagePadding->{{20,1},{20,1}},ImageSize->{120,120}],
{kk,1,6}]

DensityPlot[x,{x,-1,1},{y,0,0.1},ColorFunction->(Blend[colorRange,#]&),FrameTicks->{Automatic,None,None,None},PlotRangePadding->0,AspectRatio->0.1,ColorFunctionScaling->False,ImageSize->200]

Eigenvalue Plots Combined

Plot all temperatures

Table[
eigenvalueSTPlot=Map[EigenvalueSTPlot,evs];
eigenvalueMTPlot=Map[EigenvalueMTPlot,evs];
eigenvalueRTPlot=Map[EigenvalueRTPlot,evs];
eigenvalueHTPlot=Map[EigenvalueHTPlot,evs];
eigenvalueSTVariancePlotList =Map[EigenvalueSTVariancePlot,evs];
eigenvalueMTVariancePlotList = Map[EigenvalueMTVariancePlot,evs];
eigenvalueHTVariancePlotList =Map[EigenvalueHTVariancePlot,evs];
plot1=Show[{eigenvalueSTPlot,eigenvalueHTPlot,eigenvalueMTPlot,eigenvalueRTPlot,eigenvalueMTVariancePlotList,eigenvalueSTVariancePlotList,eigenvalueHTVariancePlotList,Graphics[Style[Text["Eigenvalue No. "<>ToString[evs[[1]]],{580,0.95},{1,1}],FontSize->20]]},Frame->True,BaseStyle->{Medium,FontFamily->"Palatino"},FrameLabel->{"Temperature [K]","Eigenvalue"},AspectRatio->0.65,Axes->False,PlotRange->{{273,600},{0.0,1.0}},Frame->True,BaseStyle->{Small,FontFamily->"Palatino",Italic},ImageSize->480]
,{evs,{{1,2,3,4,5,6},{2,3,4,5,6},{2},{3},{4},{5},{6}}}]

Export[imgPath<>"plot_Eigenvalues_Vs_Temperature_With_Reference_At_302K_no"<>Apply[StringJoin,Map[ToString,evs]]<>".pdf",plot1]

Plots of Equilibrium Distribution

Definition of Usefuls Functions

Mean1stEV[matrix_]:=#/Apply[Plus,#]&[Eigenvectors[Transpose[matrix],1][[1]]];

Error1stEV[matrix_,confidence_]:=Module[{AMatrixInv,pVec},AMatrixInv=PseudoInverse[Join[Transpose[matrix]-IdentityMatrix[numberOfStates],{{1,1,1,1,1,1}}]];
pVec=#/Apply[Plus,#]&[Eigenvectors[Transpose[matrix],1][[1]]];
Apply[Plus,Abs[Flatten[Table[pVec[[aa]]*AMatrixInv.Table[If[ii==aa,1,0]+If[ii==bb,-1,0],{ii,1,7}]*confidence[[aa,bb]],{aa,1,6},{bb,1,6}],{1,2}]]]];

Computation of Equilibrium Probabilities and Errors

varianceEquilibriumDistributionST=Table[
Error1stEV[meanST[[temp]],0.5*(confidence68ST[[temp,All,All,2]]-confidence68ST[[temp,All,All,1]])]
,{temp,1,40}];
meanEquilibriumDistributionST=Table[
Mean1stEV[meanST[[temp]]]
,{temp,1,40}];

varianceEquilibriumDistributionMT=Table[
Error1stEV[meanMT[[temp]],confidence68MT[[temp,All,All]]]
,{temp,1,79}];
meanEquilibriumDistributionMT=Table[
Mean1stEV[meanMT[[temp]]]
,{temp,1,79}];

varianceEquilibriumDistributionRT=Error1stEV[meanRT,0.5(confidence68RT[[All,All,2]]-confidence68RT[[All,All,1]])];
meanEquilibriumDistributionRT=Mean1stEV[meanRT];

Plots

eqDistSTPlot=Table[ListPlot[Transpose[{temperatureListST,meanEquilibriumDistributionST[[All,ii]]}],Joined->True,PlotStyle->{{Thick,Hue[0]}},PlotRange->All],{ii,1,6}];
eqDistMTPlot=Table[ListPlot[Transpose[{temperatureListMT,meanEquilibriumDistributionMT[[All,ii]]}],Joined->True,PlotStyle->{{Thick,Dashed,RGBColor[0,0,1]}},PlotRange->All],{ii,1,6}];
eqDistRTPlot=Table[ErrorListPlot[{{{302,meanEquilibriumDistributionRT[[ii]]},ErrorBar[varianceEquilibriumDistributionRT[[ii]]]}},PlotStyle->{Thick,GrayLevel[0]}],{ii,1,6}];

eqDistSTVariancePlotList =Table[ ListPlot[{Transpose[{temperatureListST, meanEquilibriumDistributionST[[All,ev]]-varianceEquilibriumDistributionST[[All,ev]]}], Transpose[{temperatureListST, meanEquilibriumDistributionST[[All,ev]]+varianceEquilibriumDistributionST[[All,ev]]}]}, Filling -> True, Joined -> True,PlotRange->All, FillingStyle -> {RGBColor[1, 0, 0, .1]}],{ev,1,6}];

eqDistMTVariancePlotList =Table[ ListPlot[{Transpose[{temperatureListMT, meanEquilibriumDistributionMT[[All,ev]]-varianceEquilibriumDistributionMT[[All,ev]]}], Transpose[{temperatureListMT, meanEquilibriumDistributionMT[[All,ev]]+varianceEquilibriumDistributionMT[[All,ev]]}]}, Filling -> True, Joined -> True,PlotRange->All, FillingStyle -> {RGBColor[0, 0, 1, .1]}],{ev,1,6}];

plot1=Show[{eqDistSTPlot,eqDistMTPlot,eqDistRTPlot,eqDistSTVariancePlotList,eqDistMTVariancePlotList},Frame->True,BaseStyle->{Medium,FontFamily->"Palatino"},FrameLabel->{"Temperature [K]","Equilibrium Distribution"},AspectRatio->0.65,Axes->False,PlotRange->{{273,600},{0.0,0.6}},Frame->True,BaseStyle->{Small,FontFamily->"Palatino",Italic},ImageSize->720];

plotoutput=Labeled[plot1,Style["EquilibriumDist Vs Temperature (SingleTemp(red), MultiTemp(blue), 302KReference(black)",FontSize->18,FontFamily->"Palatino"]]
Export[imgPath<>"plot_EquilibriumDist_Vs_Temperature_With_Reference_At_302K.pdf",plotoutput]

Plot of Development of the Confidence Intervals of Number of Samples

For Temperature 1(blue) ,5(yellow) and 10(red)

ConvergenceDevelopmentPlot[confidence95HTCheck]

ConvergenceDevelopmentHTPlot[6,6,5,RGBColor[0,0,1]]

Conergence For temperatures 6(blue), 23(yellow), 40(red)

picTable=Table[
(* Combine Plots *)
Show[{ConvergenceDevelopmentHTPlot[6,st1,st2,RGBColor[0,0,1]],ConvergenceDevelopmentHTPlot[23,st1,st2,RGBColor[1.0,.8,0]],ConvergenceDevelopmentHTPlot[40,st1,st2,RGBColor[1,0,0]]},ImageSize->{If[st2==1,35,0],If[st1==6,30,0]}+{200,130},BaseStyle->{FontSize->$niceFontSize,FontFamily->"Palatino"},Axes->False,AspectRatio->130/200.5,ImagePadding->{{If[st2==1,35,0],0},{If[st1==6,30,0],0}}+0.1,PlotRangeClipping->True,PlotRangePadding->Scaled[0.02],Frame->True],{st1,1,6},{st2,1,6}];

textLeftRaw="Confidence Interval of Transition Probability (Subscript[T, i -> j])";
textLeft=Graphics[Rotate[Style[Text[textLeftRaw,{0.5,0.5}],FontFamily->"Palatino",FontSize->$niceFontSize],90Degree],ImageSize->{35,600}];
textBottomRaw="Number of Drawn Samples";
textBottom=Graphics[Style[Text[textBottomRaw,{0.5,0.5}],FontFamily->"Palatino",FontSize->$niceFontSize],ImageSize->{400,35}];
plot2=Grid[Transpose[Join[{{textLeft,SpanFromAbove,SpanFromAbove,SpanFromAbove,SpanFromAbove,SpanFromAbove,SpanFromAbove}},Transpose[Join[picTable,{{textBottom,SpanFromLeft,SpanFromLeft,SpanFromLeft,SpanFromLeft,SpanFromLeft}}]]]],ItemSize->Full,Spacings->{0.1,0.1},Alignment->{Center,Center}];

plot2

plot2=Grid[Table[
(* Combine Plots *)
Show[{ConvergenceDevelopmentHTPlot[6,st1,st2,RGBColor[0,0,1]],ConvergenceDevelopmentHTPlot[23,st1,st2,RGBColor[1,.8,0]],ConvergenceDevelopmentHTPlot[40,st1,st2,RGBColor[1,0,0]]},ImageSize->{If[st2==1,40,0],If[st1==6,20,0]}+{200,130},Frame->True,BaseStyle->{FontSize->Scaled[0.070],FontFamily->"Palatino"},Axes->False,AspectRatio->130/200,ImagePadding->{{If[st2==1,40,0],0},{If[st1==6,20,0],0}}+0.2,PlotRangeClipping->True,PlotRangePadding->Scaled[0.045],Frame->True],{st1,1,6},{st2,1,6}],ItemSize->Full,Spacings->{0.1,0.1},Alignment->{Center,Center}]

Export[imgPath<>"plot_ConvergenceWithSamplesHT.pdf",plot2]

plot2=GraphicsArray[Table[Show[{ConvergenceDevelopmentHTPlot[6,st1,st2,RGBColor[0,0,1]],ConvergenceDevelopmentHTPlot[23,st1,st2,RGBColor[1,.8,0]],ConvergenceDevelopmentHTPlot[40,st1,st2,RGBColor[1,0,0]]}]
,{st1,1,6},{st2,1,6}],GraphicsSpacing->0,ImageSize->1450]

picTable=Table[
(* Combine Plots *)
Show[{ConvergenceDevelopmentSTPlot[6,st1,st2,RGBColor[0,0,1]],ConvergenceDevelopmentSTPlot[23,st1,st2,RGBColor[1,.8,0]],ConvergenceDevelopmentSTPlot[40,st1,st2,RGBColor[1,0,0]]},ImageSize->{If[st2==1,35,0],If[st1==6,30,0]}+{200,130},BaseStyle->{FontSize->$niceFontSize,FontFamily->"Palatino"},Axes->False,AspectRatio->130/200.5,ImagePadding->{{If[st2==1,35,0],0},{If[st1==6,30,0],0}}+0.1,PlotRangeClipping->True,PlotRangePadding->Scaled[0.02],Frame->True],{st1,1,6},{st2,1,6}];

textLeftRaw="Confidence Interval of Transition Probability (Subscript[T, i -> j])";
textLeft=Graphics[Rotate[Style[Text[textLeftRaw,{0.5,0.5}],FontFamily->"Palatino",FontSize->$niceFontSize],90Degree],ImageSize->{35,600}];
textBottomRaw="Number of Drawn Samples";
textBottom=Graphics[Style[Text[textBottomRaw,{0.5,0.5}],FontFamily->"Palatino",FontSize->$niceFontSize],ImageSize->{400,35}];
plot2=Grid[Transpose[Join[{{textLeft,SpanFromAbove,SpanFromAbove,SpanFromAbove,SpanFromAbove,SpanFromAbove,SpanFromAbove}},Transpose[Join[picTable,{{textBottom,SpanFromLeft,SpanFromLeft,SpanFromLeft,SpanFromLeft,SpanFromLeft}}]]]],ItemSize->Full,Spacings->{0.1,0.1},Alignment->{Center,Center}];

plot2

plot2=Grid[Table[
(* Combine Plots *)
Show[{ConvergenceDevelopmentSTPlot[6,st1,st2,RGBColor[0,0,1]],ConvergenceDevelopmentSTPlot[23,st1,st2,RGBColor[1,.8,0]],ConvergenceDevelopmentSTPlot[40,st1,st2,RGBColor[1,0,0]]},ImageSize->{If[st2==1,40,0],If[st1==6,20,0]}+{200,130},Frame->True,BaseStyle->{FontSize->Scaled[0.070],FontFamily->"Palatino"},Axes->False,AspectRatio->130/200,ImagePadding->{{If[st2==1,40,0],0},{If[st1==6,20,0],0}}+0.2,PlotRangeClipping->True,PlotRangePadding->Scaled[0.045],Frame->True],{st1,1,6},{st2,1,6}],ItemSize->Full,Spacings->{0.1,0.1},Alignment->{Center,Center}]

Export[imgPath<>"plot_ConvergenceWithSamplesST.pdf",plot2]

plot2=GraphicsArray[Table[Show[{ConvergenceDevelopmentSTPlot[6,st1,st2,RGBColor[0,0,1]],ConvergenceDevelopmentSTPlot[23,st1,st2,RGBColor[1,.8,0]],ConvergenceDevelopmentSTPlot[40,st1,st2,RGBColor[1,0,0]]}]
,{st1,1,6},{st2,1,6}],GraphicsSpacing->0,ImageSize->1450]
Export[imgPath<>"plot_ConvergenceWithSamplesST.pdf",plot2]

Plot of Development of the Energies

Mean[Map[Variance,energyList[[40]]]]

ListPlot[energyList[[6,All,1]],Frame->True]

ListPlot[energyList[[6,2,All]]]

Dimensions[energyList]

Map[StandardDeviation,energyList,{2}]

stdDevPerTrajectory=Map[StandardDeviation,energyList,{2}];

Map[Mean,stdDevPerTrajectory][[{6,23,40}]]
Map[StandardDeviation,stdDevPerTrajectory][[{6,23,40}]]

Dimensions[Transpose[Flatten[energyList,{3,2}]]]

temperatureListST[[{6,23,40}]]

driftTemp6=Table[D[Fit[energyList[[6,rr]], {1,x}, x],x]*10,{rr,1,501}];
driftTemp23=Table[D[Fit[energyList[[23,rr]], {1,x}, x],x]*10,{rr,1,501}];
driftTemp40=Table[D[Fit[energyList[[40,rr]], {1,x}, x],x]*10,{rr,1,501}];

StandardDeviation[driftTemp6]/ps
StandardDeviation[driftTemp23]/ps
StandardDeviation[driftTemp40]/ps

Mean[driftTemp6]/ps
Mean[driftTemp23]/ps
Mean[driftTemp40]/ps

minEnergy=Min[energyList];
maxEnergy=Max[energyList];

energyOverlap=Transpose[energyListRaw];

plotOutput=Show[Graphics[Table[{RGBColor[1-Mod[temp,2],0,Mod[temp,2],.5],PlotHistogram[energyOverlap[[temp]],{-3700,-1400},480]},{temp,1,40}]],Frame->True,AspectRatio->0.3,ImageSize->800]
Export[imgPath<>"energyOverlapHistogramPlot.pdf",plotOutput]

distList=Table[Apply[NormalDistribution,{Mean[#],Sqrt[Variance[#]]}&[energyOverlap[[temp]]]],{temp,1,40}];distListShort=Table[Apply[NormalDistribution,{Mean[#],Sqrt[Variance[#]]}&[energyList[[temp,All,1]]]],{temp,1,40}];

pSpecialStyle={{38,Blue},{39,Yellow},{40,Red}};
pStyleList=Table[If[MemberQ[pSpecialStyle[[All,1]],temp],Directive[Black,Thick],Black],{temp,1,numberOfTemperatures}];
fStyleList=Reverse[Table[temp->{Bottom,If[MemberQ[pSpecialStyle[[All,1]],temp],Directive[Opacity[0.7],pSpecialStyle[[Position[pSpecialStyle[[All,1]],temp][[1,1]]]][[2]]],RGBColor[0,0,0,0.15]]},{temp,1,numberOfTemperatures}]];

plotFuncList=Table[PDF[distListShort[[temp]],x],{temp,{1,2,3,4,5,7,8,9,10,11,12,13,14,15,16,17,18,19,20,21,22,24,25,26,27,28,29,30,31,32,33,34,35,36,37,38,39,6,23,40}}];
plotOutput=Plot[plotFuncList,{x,-3700,-1300},PlotRange->{{-3700,-1400},{0,0.012}},Frame->True,AspectRatio->0.3,ImageSize->800,Filling->fStyleList,PlotStyle->pStyleList,FrameLabel->{"Energy [kcal/mol]","Probability"},FrameTicks->{Automatic,None},BaseStyle->{"Palatino",16}]

Export[imgPath<>"plot_EnergyOverlap_InitialEnergy_Gauss.pdf",plotOutput]

plotFuncList=Table[PDF[distList[[temp]],x],{temp,1,40}];
plotOutput=Plot[plotFuncList,{x,-3700,-1300},PlotRange->{{-3700,-1400},{0,0.012}},Frame->True,AspectRatio->0.3,ImageSize->800,Filling->Bottom,PlotStyle->RGBColor[0,0,1],FrameLabel->{"Energy","Probablity"}]
Export[imgPath<>"plot_EnergyOverlap_AllEnergies_Gauss.pdf",plotOutput]

Plots with the states set by John

stateRegions={
{"C5",{117,-105},{28,-124},.4787,.0613},
{"PII",{-105,0},{28,-124},.4159,.0486},
{"αP",{117,-105},{-124,28},.0425,.0038},
{"αR",{-105,0},{-124,28},.0588,.0079},
{"C7",{0,117},{111,-5},.0030,.0015},
{"αL",{0,117},{-5,111},.0011,.0004}
};

Animate[
eqDist=getEquilibriumDistribution[meanMT[[temp]]];
Show[Graphics[MapIndexed[{RGBColor[eqDist[[#2[[1]]]],0,0],StateShape[#1[[2]],#1[[3]]]}&,stateRegions]],PlotRange->{{-180,180},{-180,180}},AspectRatio->1.,Frame->True,BaseStyle->{FontFamily->"Palatino"}],{temp,1,numberOfTemperatures,1}]

Computations for Sampled Rate Matrices

Show histogram of distribution of rates for each single transition 6x6 for first temperature

Show[GraphicsArray[Table[Histogram[rateMatrices6by6[[All,ii,jj]]],{ii,1,numberOfStates},{jj,1,numberOfStates}]]]

Show histogram of distribution of rates for each single transition 6x6 for first temperature in Logarithmix Plot, MarkovTime is 6 in ps

markovTime=6;
Show[GraphicsArray[Table[Histogram[Log[10,Abs[markovTime/rateMatrices6by6[[All,ii,jj]]]],HistogramRange->{0,4}],{ii,1,numberOfStates},{jj,1,numberOfStates}]]]

Show histogram of single ratematrix element

Histogram[rateMatrices6by6[[All,6,3]]]

Percentage of real rate Matrices compared to all sampled ones

Length[positiveRateMatrices6by6]/200000.

Distribution Range for the ordered spectrum of eigenvalues

ListPlot[{Map[Min,data],Map[Max,data]}]

Computation of Hummer Rate Matrices with Mathematica

Generate Matrices

hummerMatricesTable = Table[
Print[ii];
Print[Likelihood[MaximumLikelihoodMatrix[CountMatrix[tCountMatricesUT6by6[[ii]]]],CountMatrix[tCountMatricesUT6by6[[ii]]]]];
rm=ToHummerRateMatrix6[CountMatrix[tCountMatricesUT6by6[[ii]]],MaxIterations->100];
Export["/Users/jan-hendrikprinz/Studium/Projekte/dynamical-reweighting/datasets/alanine-dipeptide/parallel-tempering/hummerRateMatrix"<>ToString[ii]<>".tsv",Normal[rm]];
Print[Normal[rm]];
rm
,{ii,13,Length[temperatureListST]}];

Read Matrix from previously saved ones

tRateMatrix302KHT6by6=Import["/Users/jan-hendrikprinz/Studium/Projekte/dynamical-reweighting/datasets/alanine-dipeptide/parallel-tempering/hummerRateMatrixEstimates/mathematica/hummerRateMatrix"<>ToString[6]<>".tsv"];

Generate Optimal Solution at 302K

rateMatrixST302KNew=ToHummerRateMatrix6[CountMatrix[tCountMatricesUT6by6[[6]]],MaxIterations->2000]

Best optimization transition matrix at 302K

tMHummer302K

Likelihood

Likelihood[tMHummer302K,CountMatrix[tCountMatricesUT6by6[[6]]]]

Likelihood of the Average Sampled Matrix from Hummer Estimates

Likelihood[TransitionMatrix[MatrixExp[meanHTRate[[6]]*6.0]],CountMatrix[tCountMatricesUT6by6[[6]]]]

Lower Boundary from uncontrained problem

Likelihood[MaximumLikelihoodMatrix[CountMatrix[tCountMatricesUT6by6[[6]]]],CountMatrix[tCountMatricesUT6by6[[6]]]]

Comparison of TimeScales and Eigenvalues

ImpliedTimeScales[TransitionMatrix[tMatrixRT6by6]]
ImpliedTimeScales[TransitionMatrix[meanMT[[11]]]]
ImpliedTimeScales[tMHummer302K]
ImpliedTimeScales[TransitionMatrix[meanHT[[6]]]]

Spectrum[TransitionMatrix[N[tMatrixRT6by6]]]
Spectrum[TransitionMatrix[N[meanMT[[11]]]]]
Spectrum[tMHummer302K]
Spectrum[TransitionMatrix[N[meanHT[[6]]]]]

Plot for Sampling of Hummer

plot1=GraphicsArray[Table[
meanSTPlot=MeanSTPlot[st1,st2];
confidence95STPlot=Confidence95STPlot[st1,st2];
meanMTPlot=MeanMTPlot[st1,st2];
confidence95MTPlot=Confidence95MTPlot[st1,st2];
meanHTPlot=MeanHTPlot[st1,st2];
confidence95HTPlot=Confidence95HTPlot[st1,st2];
maxLikeSTPlot=MaxLikeSTPlot[st1,st2];
maxLikeRTPlot=MaxLikeRTPlot[st1,st2];
Show[{confidence95STPlot,meanSTPlot,confidence95MTPlot,meanMTPlot,confidence95HTPlot,meanHTPlot,maxLikeRTPlot,Graphics[Text[Style[Subscript["T",ToString[st1]<>" -> "<>ToString[st2]],FontSize->18],{600,1.0},{1,1}]]},Frame->True,FrameTicksStyle->Tiny,BaseStyle->{FontFamily->"Helvetica"},AspectRatio->0.65,Axes->False,PlotRange->{{273,600},{0.0,1.0}},Frame->True,BaseStyle->{Small,FontFamily->"Helvetica",Italic}],{st1,1,6},{st2,1,6}],ImageSize->1450,GraphicsSpacing->0.02]

Export[imgPath<>"plot_Probility_Vs_Temperature_All_Transitions_Without_MaxLikelihood_plusHummer.pdf",plot1]

plot1=GraphicsArray[Table[
meanSTPlot=MeanSTPlot[st1,st2];
confidence95STPlot=Confidence95STPlot[st1,st2];
meanMTPlot=MeanMTPlot[st1,st2];
confidence95MTPlot=Confidence95MTPlot[st1,st2];
meanHTPlot=MeanHTPlot[st1,st2];
confidence95HTPlot=Confidence95HTPlot[st1,st2];
maxLikeSTPlot=MaxLikeSTPlot[st1,st2];
maxLikeRTPlot=MaxLikeRTPlot[st1,st2];
Show[{confidence95STPlot,meanSTPlot,confidence95MTPlot,meanMTPlot,confidence95HTPlot,maxLikeRTPlot,Graphics[Text["Probability "Subscript["T",ToString[st1]<>" -> "<>ToString[st2]],{600,0.95},{1,0}]]},Frame->True,BaseStyle->{FontFamily->"Palatino"},AspectRatio->0.65,Axes->False,PlotRange->{{273,600},{0.0,1.0}},Frame->True,BaseStyle->{Small,FontFamily->"Palatino",Italic}],{st1,1,6},{st2,1,6}],ImageSize->1450,GraphicsSpacing->0]

Export[imgPath<>"plot_Probility_Vs_Temperature_All_Transitions_Without_MaxLikelihood.pdf",plot1];

Manipulate[Show[{MeanSTPlot[state1,state2],MeanHTPlot[state1,state2],MeanMTPlot[state1,state2],Confidence95HTPlot[state1,state2],Confidence95MTPlot[state1,state2],Confidence95STPlot[state1,state2],MaxLikeRTPlot[state1,state2],Graphics[Text[ToString[state1]<>" ⟶ "<>ToString[state2],{600,0.95},{1,1}]]},Frame->True,PlotRange->{0,1},ImageSize->640],{state1,1,6,1},{state2,1,6,1}]

Dimensions[tCountMatrixRT6by6]

tCountMatrixRT6by6

tMatrixRT6by6

Plot of Energy Landscape

wMat=wMatrix302K;
ramachandranBins=N[BinCounts[wMat[[All,{4,5}]],{-180,180,10},{-180,180,10}]/Length[wMat]];
ramachandranBins=ramachandranBins/Max[ramachandranBins];

maxkBT=10.0;
minVal=Exp[-maxkBT];
ramachandran=Map[If[#>minVal,-Log[#],maxkBT]&,ramachandranBins,{2}];
Max[ramachandran]
Min[ramachandran]
MatrixPlot[Transpose[ramachandran],DataReversed->True,Mesh->True,MeshStyle->Black,PlotRange->{0,10}]

MatrixPlot[Transpose[ramachandran],ColorFunction->colFunc,ColorFunctionScaling->True,DataReversed->True,Mesh->True,MeshStyle->Black,PlotRange->{-10,10}]

ramachandranFunction=ListInterpolation[ramachandran,{{-π,π},{-π,π}},InterpolationOrder->0];

Plot3D[ramachandranFunction[phi,psi],{phi,-π,π},{psi,-π,π},PlotRange->All,ColorFunction->"Rainbow",Axes->True,Mesh->False]

Graphics[Table[{colFunc[nn],Disk[{8nn,0}]},{nn,0,1,0.02}]]

colFunc=Blend[Join[Table[Hue[ww],{ww,0,0.9,0.1}],{White}],#]&;

ListDensityPlot[Transpose[ramachandran],ColorFunction->colFunc,InterpolationOrder->1,Mesh->36,MeshStyle->Black]

DensityPlot[ramachandranFunction[phi,psi],{phi,-π,π},{psi,-π,π},PlotRange->All,ColorFunction->"TemperatureMap",FrameTicks->{{{-Pi,0,Pi},None},{{-Pi,0,Pi},None}},Axes->True,PlotPoints->12]

ramPlot=ContourPlot[ramachandranFunction[phi,psi],{phi,-π,π},{psi,-π,π},PlotRange->All,Contours->15,ColorFunction->"Rainbow",FrameTicks->{{{-Pi,0,Pi},None},{{-Pi,0,Pi},None}},Axes->True]

allram=Apply[Plus,Flatten[ramachandran]];
ramSort=Sort[Flatten[ramachandran]];
ramSumSort=Table[Apply[Plus,ramSort[[Range[ii]]]],{ii,1,37*37}];
ramCoreRegion=(#[[2]]-#[[1]])&[Select[ramSumSort,#>=allram*(1-0.819)&][[{1,2}]]];
ramAllowedRegion=(#[[2]]-#[[1]])&[Select[ramSumSort,#>=allram*(1-0.819-0.148)&][[{1,2}]]];
ramGenerousRegion=(#[[2]]-#[[1]])&[Select[ramSumSort,#>=allram*(1-0.819-0.148-0.020)&][[{1,2}]]];

ramPlot=ContourPlot[ramachandranFunction[phi,psi],{phi,-π,π},{psi,-π,π},PlotRange->All,PlotPoints->10,Contours->{ramGenerousRegion,ramAllowedRegion,ramCoreRegion},ContourShading->{White,Yellow,Orange,Red},FrameTicks->{{{-Pi,0,Pi},None},{{-Pi,0,Pi},None}},Axes->True]

ramachandranFine=Table[ramachandranFunction[psi,phi],{phi,-π,π,π/128},{psi,-π,π,π/128}];

ReliefPlot[ramachandranFine,ColorFunction->"GreenBrownTerrain"]

Plot for log(W[i,j,k,l]) at given i,j,k vs l and Plot of curvature dependent on all temperatures

FrameTicks->{Map[{#,ToString[Round[100*#]]}&,Range[0,1,0.2]],None,None,None}

ContributionPlot[matrix_]:=Module[{},
legendPlot=Show[DensityPlot[x,{x,0,1},{y,0,0.1},Frame->True,AspectRatio->0.1,FrameTicks->{Automatic,None,None,None},FrameLabel->{{None,None},{None,"Counts"}},ColorFunction->(ColorData["SunsetColors"][1-(#)]&),ImageSize->140,PlotRangePadding->0.01],Frame->True];
tickList=Transpose[{Range[40]+.5,Map[ToString[#]<>"K"&,Round[temperatureListST]]}][[{6,20,31,40}]];
plot1=Graphics[Flatten[Table[{ColorData["SunsetColors"][1.0-#]&[matrix[[ii,jj]]],Rectangle[{jj,ii}+0.05,{jj,ii}+0.95]},{jj,1,40},{ii,1,40}]],Frame->True,FrameTicks->{tickList,tickList,None,None},PlotRangePadding->Scaled[0.01],PlotRange->{{1,41},{1,41}},FrameTicks->{tickList,tickList},FrameLabel->{"Temperature contributed from","Temperature estimated at"},ImageSize->240,Epilog->{Inset[legendPlot,{41,1},{Right,Bottom},30],Text[Style[Row[{ToString[st1],"->",ToString[st2]}," "],FontSize->Scaled[0.08]],Scaled[{0.05,0.95}],{-1,1}]},BaseStyle->{FontSize->11,FontFamily->"Palatino"}]]

ContributionPlotPlain[matrix_]:=Module[{},
legendPlot=Show[DensityPlot[x,{x,0,1},{y,0,0.1},Frame->True,AspectRatio->0.1,FrameTicks->{Map[{#,ToString[Round[100*#]]<>"%"}&,Range[0,1,0.2]],None,None,None},ColorFunction->(ColorData["SunsetColors"][1-(#)]&),ImageSize->140,PlotRangePadding->0.01],Frame->True];
tickList=Transpose[{Range[40]+.5,Map[ToString[#]<>"K"&,Round[temperatureListST]]}][[{6,20,31,40}]];
plot1=Graphics[Flatten[Table[{ColorData["SunsetColors"][1.0-#]&[matrix[[ii,jj]]],Rectangle[{jj,ii}+0.05,{jj,ii}+0.95]},{jj,1,40},{ii,1,40}]],Frame->True,FrameTicks->{tickList,tickList,None,None},PlotRangePadding->Scaled[0.01],PlotRange->{{1,41},{1,41}},FrameTicks->{tickList,tickList},FrameLabel->{"Temperature contributed from","Temperature estimated at"},ImageSize->240,Epilog->{Inset[legendPlot,{41,1},{Right,Bottom},30]},BaseStyle->{FontSize->11,FontFamily->"Palatino"}]]

ContributionPlotPlainLog[matrix_]:=Module[{},
legendPlot=Show[DensityPlot[x,{x,0,1},{y,0,0.1},Frame->True,AspectRatio->0.1,FrameTicks->{Map[{#,Superscript["10",ToString[Round[-(1-#)*10]]]}&,Range[0,1,0.2]],None,None,None},ColorFunction->(ColorData["SunsetColors"][1-(#)]&),ImageSize->140,PlotRangePadding->0.01],Frame->True];
tickList=Transpose[{Range[40]+.5,Map[ToString[#]<>"K"&,Round[temperatureListST]]}][[{6,20,31,40}]];
plot1=Graphics[Flatten[Table[{ColorData["SunsetColors"][1.0-#]&[matrix[[ii,jj]]],Rectangle[{jj,ii}+0.05,{jj,ii}+0.95]},{jj,1,40},{ii,1,40}]],Frame->True,FrameTicks->{tickList,tickList,None,None},PlotRangePadding->Scaled[0.01],PlotRange->{{1,41},{1,41}},FrameTicks->{tickList,tickList},FrameLabel->{"Temperature contributed from","Temperature estimated at"},ImageSize->240,Epilog->{Inset[legendPlot,{41,1},{Right,Bottom},30]},BaseStyle->{FontSize->11,FontFamily->"Palatino"}]]

st1=6;
st2=5;
ContributionPlot[Table[Chop[Exp[logWijkl[[st1,st2,temp]]]/Apply[Plus,Exp[logWijkl[[st1,st2,temp]]]]],{temp,1,40}]]

Transpose[Exp[logWijkl[[4,6]]]]

Map[Chop[#/Apply[Plus,#]]&,Transpose[Exp[logWijkl[[4,6]]]]]

(tCountMatricesUT6by6[[All,5,1]]+tCountMatricesUT6by6[[All,1,5]])/2

Exp[logWkl[[40,40]]]

Apply[Plus,Exp[logWkl[[temp]]]]

NumberForm[mat1,{1,0}]

matrix[[40]]*3.5

st1=6;
st2=6;
Table[
mat1=Table[Chop[Exp[logWijkl[[st1,st2,temp]]]/Apply[Plus,Exp[logWkl[[temp]]]]],{temp,1,40}]*501;
mat2=mat1/Max[mat1];
plot1=ContributionPlot[mat1];
Export[imgPath<>"plot_contribution_to_estimate_W"<>ToString[st1]<>ToString[st2]<>"kl.pdf",plot1];
,{st1,1,6},{st2,1,6}];

st1=4;
st2=6;
tableToPlot=Transpose[Table[Chop[Exp[logWkl[[temp]]]/Apply[Plus,Exp[logWkl[[temp]]]]],{temp,1,40}]];
tableToPlot=1-Map[Max[Log[#],-10]&,tableToPlot,{2}]/(-10);
plot1=ContributionPlotPlainLog[tableToPlot]

Export[imgPath<>"plot_contribution_to_summary.pdf",plot1]

tickList=Transpose[{Range[40],Map[ToString[#]<>"K"&,Round[temperatureListST]]}][[{6,20,31,40}]]

plot1=ListPlot[Table[Chop[Exp[logWkl[[temp]]]/Apply[Plus,Exp[logWkl[[temp]]]]],{temp,{6,20,31,40}}],Joined->True,PlotRange->All,Frame->True,PlotStyle->{{Thick,Black}},FrameTicks->{tickList,Automatic}]

logWijkl

logWijkl[[All,All,1,1]]

st1=1;
st2=1;
Table[
matrix=Transpose[Table[Chop[Exp[logWijkl[[st1,st2,temp]]]/Apply[Plus,Exp[logWijkl[[st1,st2,temp]]]]],{temp,1,40}]];
plot1=ContributionPlot[matrix];
Export[imgPath<>"plot_contribution_to_estimate_W"<>ToString[st1]<>ToString[st2]<>"kl.pdf",plot1],
{st1,5,6},{st2,5,6}];

matrix=Transpose[Chop[Map[#/Apply[Plus,#]&,Exp[logWkl]]]];
matrix=1-Map[Max[Log[#],-10]&,matrix,{2}]/(-10);
legendPlot=Show[DensityPlot[x,{x,0,1},{y,0,0.1},Frame->True,AspectRatio->0.1,FrameTicks->{Map[{#,ToString[Round[100*#]]<>"%"}&,Range[0,1,0.2]],None,None,None},ColorFunction->(ColorData["SunsetColors"][1-(#)]&),ImageSize->140,PlotRangePadding->Scales[0.01]],Frame->True];
legendPlot=Show[DensityPlot[x,{x,0,1},{y,0,0.1},Frame->True,AspectRatio->0.1,FrameTicks->{Map[{#,Superscript["10",ToString[Round[-(1-#)*10]]]}&,Range[0,1,0.2]],None,None,None},ColorFunction->(ColorData["SunsetColors"][1-(#)]&),FrameLabel->{{None,None},{None,"Probability"}},ImageSize->140,PlotRangePadding->0.01],Frame->True];
tickList=Transpose[{Range[40]+.5,Map[ToString[#]<>"K"&,Round[temperatureListST]]}][[{6,20,31,40}]];
plot1=Graphics[Flatten[Table[{ColorData["SunsetColors"][1.0-#]&[matrix[[ii,jj]]],Rectangle[{jj,ii}+0.05,{jj,ii}+0.95]},{jj,1,40},{ii,1,40}]],Frame->True,FrameTicks->{tickList,tickList,None,None},PlotRangePadding->Scaled[0.01],PlotRange->{{1,41},{1,41}},FrameTicks->{tickList,tickList},FrameLabel->{"Temperature estimated at","Temperature contributed from"},ImageSize->240,Epilog->{Inset[legendPlot,{41,1},{Right,Bottom},30]},BaseStyle->{FontSize->11,FontFamily->"Palatino"}]
Export[imgPath<>"plot_normalized_weights_Wkn.pdf",plot1]

.BIBLIOGRAPHY

Literature Links, etc.

.APPENDIX

Additional Material

Using eigenvalue permutation based on eigenvectors

getSymmetrizedEV[matrix_] := Module[{eqDist},
  eqDist = (#/Apply[Plus, #]) &[Eigenvectors[Transpose[matrix]][[1]]];
  Eigenvectors[matrix.DiagonalMatrix[1/eqDist]]
  ]

getSymmetrizedEV[matrix_] := Module[{},
  Eigenvectors[matrix]
  ]

numberOfRelevantModes=5;

modeComparison =
     Table[
symmEVST=getSymmetrizedEV[meanST[[temp]]];
symmEVMT=getSymmetrizedEV[meanMT[[2*temp-1]]];
Outer[Dot,symmEVST[[Range[2,numberOfRelevantModes+1]]],symmEVMT[[Range[2,numberOfRelevantModes+1]]],1]^2
,{temp, 1, numberOfTemperatures}];
modeComparisonST=
     Table[
symmEVST=getSymmetrizedEV[meanST[[temp]]];
symmEVST2=getSymmetrizedEV[meanST[[temp+1]]];
Outer[Dot,symmEVST[[Range[2,numberOfRelevantModes+1]]],symmEVST2[[Range[2,numberOfRelevantModes+1]]],1]^2
,{temp, 1, numberOfTemperatures-1}];
modeComparisonMT=
     Table[
symmEVMT=getSymmetrizedEV[meanMT[[temp]]];
symmEVMT2=getSymmetrizedEV[meanMT[[temp+2]]];
Outer[Dot,symmEVMT[[Range[2,numberOfRelevantModes+1]]],symmEVMT2[[Range[2,numberOfRelevantModes+1]]],1]^2
,{temp, 1, 2*numberOfTemperatures-3,2}];

allperm=Permutations[Range[numberOfRelevantModes]];

oldperm=Range[numberOfRelevantModes];
permTableST=Table[
tempTab=Table[{ii,Apply[Plus,Table[modeComparisonST[[jj,kk,allperm[[ii,kk]]]],{kk,1,numberOfRelevantModes}]]},{ii,1,Length[allperm]}];
tempMax=Max[tempTab[[All,2]]];
changeperm=allperm[[Select[tempTab,#[[2]]==tempMax&][[1,1]]]];
oldperm=oldperm[[changeperm]],{jj,1,39}]
oldperm = Range[numberOfRelevantModes];
permTableMT = Table[
  tempTab = Table[{ii, Apply[Plus, Table[modeComparisonMT[[jj, kk, allperm[[ii, kk]]]], {kk, 1, numberOfRelevantModes}]]}, {ii, 1, Length[allperm]}];
  tempMax = Max[tempTab[[All, 2]]];
  changeperm = allperm[[Select[tempTab, #[[2]] == tempMax &][[1, 1]]]];
  oldperm = oldperm[[changeperm]], {jj, 1, 39}]
oldperm = Range[numberOfRelevantModes];
permTable = Table[
  tempTab = Table[{ii, Apply[Plus, Table[modeComparison[[jj, kk, allperm[[ii, kk]]]], {kk, 1,numberOfRelevantModes}]]}, {ii, 1, Length[allperm]}];
  tempMax = Max[tempTab[[All, 2]]];
  changeperm = allperm[[Select[tempTab, #[[2]] == tempMax &][[1, 1]]]]
, {jj, 1, 39}]

ListPlot[Transpose[permTableST],Joined->True]
ListPlot[Transpose[permTableMT],Joined->True]
ListPlot[Transpose[permTable],Joined->True]

impliedTimescalesSTPlot = ListPlot[Table[Table[{temperatureListST[[temp]], impliedTimescalesST[[temp, permTableST[[temp,ii]]]]}, {temp,1,numberOfTemperatures-1}],{ii, 1, 5}], Joined -> True, PlotStyle->Thick, PlotRange -> All];
impliedTimescalesMTPlot = ListPlot[Table[Table[{temperatureListMT[[2*temp-1]], impliedTimescalesMT[[2*temp-1, permTableMT[[temp,ii]]]]}, {temp,1,numberOfTemperatures-1}],{ii, 1, 5}], Joined -> True, PlotStyle->Thick, PlotRange -> All];
impliedTimescalesPlot = ListPlot[Table[Table[{temperatureListST[[temp]], impliedTimescalesMT[[2*temp-1, permTable[[temp,ii]]]]}, {temp,1,numberOfTemperatures-1}],{ii, 1, 5}], Joined -> True, PlotStyle->Thick, PlotRange -> All];

plot1 = Show[{impliedTimescalesSTPlot,impliedTimescalesMTPlot}, Frame -> True, BaseStyle -> {Medium, FontFamily -> "Palatino"}, FrameLabel -> {"Temperature [K]", "Implied Timescales"}, AspectRatio -> 0.65, Axes -> False, PlotRange -> {{273, 600}, {0.0, 80.0}}, Frame -> True, BaseStyle -> {Small, FontFamily -> "Palatino", Italic}, ImageSize -> 400]

analyzedMode=1;
oldperm=Range[numberOfRelevantModes];
permTableST=Table[
tempTabST=Table[{ii,Apply[Plus,Table[modeComparisonST[[jj,kk,allperm[[ii,kk]]]],{kk,analyzedMode,analyzedMode}]]},{ii,1,Length[allperm]}];
tempMaxST=Max[tempTabST[[All,2]]];
changeperm=allperm[[Select[tempTabST,#[[2]]==tempMaxST&][[1,1]]]];
oldperm=oldperm[[changeperm]],{jj,1,39}]
oldperm = Range[numberOfRelevantModes];
permTable = Table[
 tempTab = Table[{ii, Apply[Plus, Table[modeComparison[[jj, kk, allperm[[ii, kk]]]], {kk, analyzedMode,analyzedMode}]]}, {ii, 1, Length[allperm]}];
  tempMax = Max[tempTab[[All, 2]]];
  changeperm = allperm[[Select[tempTab, #[[2]] == tempMax &][[1,1]]]]
, {jj, 1, 39}]
ListPlot[{Table[{temperatureListST[[temp]], impliedTimescalesST[[temp, permTableST[[temp,analyzedMode]]]]}, {temp,1,numberOfTemperatures-1}],Table[{temperatureListST[[temp]], impliedTimescalesMT[[2*temp-1, permTable[[temp,permTableST[[temp,analyzedMode]]]]]]}, {temp,1,numberOfTemperatures-1}]}, Frame -> True, BaseStyle -> {Medium, FontFamily -> "Palatino"}, FrameLabel -> {"Temperature [K]", "Implied Timescales"}, AspectRatio -> 0.65, Axes -> False, PlotRange -> {{273, 600}, {0.0, 100.0}}, Frame -> True, BaseStyle -> {Small, FontFamily -> "Palatino", Italic}, ImageSize -> 400,Joined->True,PlotStyle->Thick]

Table[{temperatureListST[[temp]], impliedTimescalesMT[[2*temp-1, permTable[[temp,permTableST[[temp,analyzedMode]]]]]]}, {temp,1,numberOfTemperatures-1}]

impliedTimescalesPlot = ListPlot[Table[Table[{temperatureListST[[temp]], impliedTimescalesMT[[2*temp-1, permTable[[temp,ii]]]]}, {temp,1,numberOfTemperatures-1}],{ii, 1, 5}], Joined -> True, PlotStyle->Thick, PlotRange -> All]

Plot for log(W[i,j,k,l]) at given i,j,k vs l and Plot of curvature dependent on all temperatures

ArrayPlot[Transpose[matrix],Mesh->{True,False},DataRange->{{273,600},{273,600}},Frame->True,FrameTicks->Automatic,DataReversed->True,MeshStyle->Directive[White],ColorFunction->"TemperatureMap"]

st1=4;
st2=6;
tickList=Transpose[{Range[40],Map[ToString[#]<>"K"&,Round[temperatureListST]]}][[{6,20,31,40}]];
matrix=Transpose[Reverse[Table[Chop[Exp[logWijkl[[st1,st2,temp]]]/Apply[Plus,Exp[logWijkl[[st1,st2,temp]]]]],{temp,1,40}]]];
plot1=MatrixPlot[Reverse[Transpose[matrix]]/.{0->White},Mesh->{True,True},Frame->True,FrameTicks->{tickList,tickList},MeshStyle->Directive[AbsoluteThickness[1],White],ColorFunction->(ColorData["SunsetColors"][1-2(#-0.5)]&),DataReversed->True,FrameLabel->{"Temperature estimated at","Temperature contributed from"},ImageSize->240,Epilog->Text[Style[Row[{ToString[st1]1,"->",ToString[st2]}," "],FontSize->Scaled[0.08]],Scaled[{0.05,0.95}],{-1,1}],BaseStyle->{FontSize->Small,FontFamily->"Palatino"}]

st1=1;
st2=1;
Table[
tickList=Transpose[{Range[40],Map[ToString[#]<>"K"&,Round[temperatureListST]]}][[{6,20,31,40}]];
matrix=Transpose[Reverse[Table[Chop[Exp[logWijkl[[st1,st2,temp]]]/Apply[Plus,Exp[logWijkl[[st1,st2,temp]]]]],{temp,1,40}]]];
plot1=MatrixPlot[Reverse[Transpose[matrix]]/.{0->White},Mesh->{True,True},Frame->True,FrameTicks->{tickList,tickList},MeshStyle->Directive[AbsoluteThickness[1],White],ColorFunction->(ColorData["SunsetColors"][1-2(#-0.5)]&),DataReversed->True,FrameLabel->{"Temperature estimated at","Temperature contributed from"},ImageSize->240,Epilog->Text[Style[Row[{ToString[st1]1,"->",ToString[st2]}," "],FontSize->Scaled[0.08]],Scaled[{0.05,0.95}],{-1,1}],BaseStyle->{FontSize->Small,FontFamily->"Palatino"}];
Export[imgPath<>"plot_contribution_to_estimate_W"<>ToString[st1]<>ToString[st2]<>"kl.pdf",plot1],
{st1,1,6},{st2,1,6}];

temperature=6;
st1=3;
st2=6;
plot1=GraphicsArray[Table[
acceptedPoints=Select[Transpose[{1/temperatureListST,logWijkl[[st1,st2,temperature]]-Max[logWijkl[[st1,st2,temperature]]]}],(#[[2]]>-900)&];
fitFunc=Fit[acceptedPoints, {1,x,x^2},x];
fitPlot=Plot[fitFunc,{x,1/temperatureListST[[1]],1/temperatureListST[[numberOfTemperatures]]},PlotStyle->{RGBColor[1,0,0],Thick,Dashed}];
dataPlot=ListPlot[acceptedPoints,Joined->True,PlotStyle->{Thick},Epilog->{Thick,Line[{{temperature,-4000},{temperature,0}}]}];
Show[dataPlot,fitPlot,PlotRange->{All,{-600,0}},Frame->True,Axes->False,AspectRatio->0.65,FrameLabel->{"1/Temperature [1/K]",""}],{st1,1,6},{st2,1,6}],ImageSize->1450,GraphicsSpacing->-0.01];
plotoutput=Labeled[plot1,Style["log(W) Vs Temperature",FontSize->18,FontFamily->"Palatino"]]
Export[imgPath<>"plot_logW_Vs_Temperature_All_Transitions.pdf",plotoutput]

contributionsSpreadTable=Table[
 acceptedPoints = Select[Transpose[{Range[numberOfTemperatures], logWijkl[[st1, st2, temperature]] - Max[logWijkl[[st1, st2, temperature]]]}], (#[[2]] > -900) &];
 fitFunc =c/.FindFit[acceptedPoints, c*x^2+b*x+a,{a,b,c}, x]
, {temperature,1,numberOfTemperatures},{st1, 1, numberOfStates}, {st2, 1, numberOfStates}];

GraphicsArray[Table[
ListPlot[contributionsSpreadTable[[All,st1,st2]],PlotRange->{-0.8,-0.4}]
,{st1,1,6},{st2,1,6}],ImageSize->1200,GraphicsSpacing->-0.01]

Additional Calculations

tCountMatricesUT6by6[[6]]

tMatricesHT6by6[[6,1]]

likelihoodList=Map[Likelihood[TransitionMatrix[#],CountMatrix[tCountMatricesUT6by6[[6]]]]&,tMatricesHT6by6[[6]]];

ListPlot[likelihoodList]

Mean[likelihoodList]

likelihoodList[[1]]

rM1=ToRateMatrix[tMHummer302K]

rMat=rM1//Normal

rMat[[5]]=rMat[[5]]/3.6

TransitionMatrix[meanHT[[6]]]

tMHummer302K

.SANDBOX

tCountMatricesUT6by6[[13]]//MatrixForm

tempfunc3=273*Exp[Log[600/273]/78.*(x-1)];

func1=a1*Exp[b1*x];
res1=FindFit[tCountMatricesUT6by6[[All,1,1]],func1,{a1,b1},x];
p1=Plot[func/.res,{x,1,40},PlotRange->{0,200}];
p2=ListPlot[tCountMatricesUT6by6[[All,1,1]],PlotRange->{0,200}];
Show[{p1,p2}]

st=11;
func2=Exp[-(x-a2)^2/b2];
res2=FindFit[Exp[logWkl[[st]]-logWkl[[st,st]]],func2,{b2,{a2,st}},x];
p1=Plot[func2/.res2,{x,1,40},PlotRange->{0,1}];
p2=ListPlot[Exp[logWkl[[st]]-logWkl[[st,st]]],PlotRange->{0,1}];
Show[{p1,p2}]

stateTab=Table[
func1=a1*Exp[b1*x];
res1=FindFit[tCountMatricesUT6by6[[All,st1,st2]],func1,{a1,b1},x];
p1=Plot[func1/.res1,{x,1,40},PlotRange->{0,200}];
p2=ListPlot[tCountMatricesUT6by6[[All,st1,st2]],PlotRange->{0,200}];
Show[{p1,p2}];
res1
,{st1,1,6},{st2,1,6}]

tempTab=Table[
func2=Exp[-(x-a2)^2/b2];
res2=FindFit[Exp[logWkl[[temp]]-logWkl[[temp,temp]]],func2,{b2,{a2,temp}},x];
p1=Plot[func2/.res2,{x,1,40},PlotRange->{0,1}];
p2=ListPlot[Exp[logWkl[[temp]]-logWkl[[temp,temp]]],PlotRange->{0,1}];
Show[{p1,p2}];
res2
,{temp,1,40}]

maxCombTab=Table[(1/2) (2 a2+b1 b2)-temp/.tempTab[[temp]]/.stateTab[[st1,st2]],{temp,1,40},{st1,1,6},{st2,1,6}];

Max[maxCombTab]

Dimensions[Map[Mean,Transpose[maxCombTab,{2,1,3}]]]

sweetSpotMatrix=Map[Exp[Log[600/273]/78.*#]&,Map[Mean,Transpose[maxCombTab,{2,1,3}]]]

MatrixPlot[sweetSpotMatrix-1,ColorFunction->"TemperatureMap"]

Map[Mean,Transpose[maxCombTab,{2,1,3}]]//MatrixForm

Map[Exp[Log[600/273]/78.*#]&,maxCombTab[[6]]]//MatrixForm

tempfunc3=273*Exp[Log[600/273]/78.*(x-1)];
tempFit=FindFit[temperatureListMT,func3,{{b0,0.0202}},x]
p1=Plot[func3/.tempFit,{x,1,79},PlotRange->All];
p2=ListPlot[temperatureListMT,PlotRange->{0,1}];
Show[{p1,p2}]

tempfunc3/.x->maxComb

600/273

Exp[beta*39]==600/273

ListPlot[temperatureListMT]

temperatureListMT[[11]]

val1=Table[(func/.res)/.(x->i),{i,1,40}];

func1*func2

Together[,x]

-b2 b1 x+(((-a2+x))^2)//Expand

(x-(2a2+b1 b2)/2)^2-((2a2+b1 b2)/2)^2+a2^2

D[(a2^2)-(1/4) (((2 a2+b1 b2))^2)+((((1/2) (-2 a2-b1 b2)+x))^2),x]

(1/2) (2 a2+b1 b2)==x

res1

res2

CountMatrix[tCountMatricesUT6by6[[1]]]

TransitionMatrix[tMatricesHT6by6[[temp,1]]]

Table[
{Likelihood[MaximumLikelihoodMatrix[CountMatrix[tCountMatricesUT6by6[[temp]]]],CountMatrix[tCountMatricesUT6by6[[temp]]]],
Likelihood[TransitionMatrix[tMatricesHT6by6[[temp,1]]],CountMatrix[tCountMatricesUT6by6[[temp]]]]},
{temp,1,numberOfTemperaturesHT}]