Calculate covariance matrix and calculate the eigenvectors and eigenvalues.

gmx covar -s md.gro -f mdfit.xtc -o eigenvalues.xvg -v eigenvectors.trr -xpma covapic.xpm

Calculate PC1 and PC2
gmx anaeig -f md.xtc -s md.gro -v eigenvectors.trr -last 1 -proj pc1.xvg
gmx anaeig -f md.xtc -s md.gro -v eigenvectors.trr -frist 2 -last 2 -proj pc2.xvg

Concatenate PC1 and PC2 in one file.

paste pc1.xvg pc2.xvg  | awk '{print $1, $2, $4}' > PC1PC2.xvg

Calculate Gibbs Free Energy.
gmx sham -f PC1PC2.xvg -ls FES.xpm

Convert .xpm file to .dat file
python2.7 xpm2txt.py -f FES.xpm -o free-energy-landscape.dat

Plot FEL in gnuplot
Activate gnuplot
gnuplot <<EOL

Copy the following and paste into kwrite or notepad

reset
set terminal pngcairo  background "#ffffff" enhanced font "Arial-Bold,18" fontscale 1.0 size 720, 720
set key on b c outside horizontal
set output 'FEL.png'
set title "Free Energy Landscape"
set xlabel "PC1" rotate parallel
set ylabel "PC2" rotate parallel
set zlabel "Gibbs Free Energy" rotate parallel
set grid
unset key
set pm3d implicit at b
set view map scale 1
unset surface
set dgrid3d 50,50
set hidden3d
set contour s
stats 'FEL.dat' u 1:(column(2)) name "A"
stats 'FEL.dat' u 1:(column(2)) name "B"
set xrange [PC1_min:PC1_max]
set yrange [PC2_min:PC2_max]
sp 'FEL.dat' u 1:2:3 w l