Expand data frame into combinations of row pairs

Your combn was a good way to go. Try this:

 combos<-combn(Test_data$SYS,2)
 water<-combn(Test_data$dE_water_free,2,FUN=function(x) x[1]-x[2])
 data.frame(SYS1=combos[1,],SYS2=combos[2,],water,stringsAsFactors=FALSE)
 #         SYS1       SYS2    water
 #1  4NTJ_D294N  4NTJ_wild  114.489
 #2  4NTJ_D294N 4PXZ_D294N   -3.112
 #3  4NTJ_D294N  4PXZ_wild    3.448
 #4  4NTJ_D294N 4PY0_D294N   20.498
 #5  4NTJ_D294N  4PY0_wild   22.538
 ........

outer is well suited for this type of problem:

de_wf <- with(Test_data, setNames(dE_water_free, SYS))
outer(de_wf, de_wf, `-`)

produces:

           4NTJ_D294N 4NTJ_wild 4PXZ_D294N 4PXZ_wild 4PY0_D294N 4PY0_wild
4NTJ_D294N      0.000   114.489     -3.112     3.448     20.498    22.538
4NTJ_wild    -114.489     0.000   -117.601  -111.041    -93.991   -91.951
4PXZ_D294N      3.112   117.601      0.000     6.560     23.610    25.650
4PXZ_wild      -3.448   111.041     -6.560     0.000     17.050    19.090
4PY0_D294N    -20.498    93.991    -23.610   -17.050      0.000     2.040
4PY0_wild     -22.538    91.951    -25.650   -19.090     -2.040     0.000