From  http://www.strodel.info/index_files/lecture/html/analysis-10.html

Free Energy Surfaces (FES)

Biomolecular processes, such as folding or aggregation, can be described in terms of the molecule's free energy: 

where k_B is the Boltzmann constant, P is the probability distribution of the molecular system along some coordinate R, and P_max denotes its maximum, which is substracted to ensure that δG = 0 for the lowest free energy minimum. Popular choices for R (so-called order parameters) are the RMSD, Rgyr, number of hydrogen bonds or native contacts, or principal components.

The free energy is typically plotted along two such order parameters, giving rise to a (reduced) free energy surface (FES). 

Example: FES (in kcal/mol) for the aggregation of the GNNQQNY peptide forming a dimer in terms of Rgyr and the RMSD from a perfect parallel β-sheet

Generating FESs with generateFES.py

python generateFES.py -f rgyr_rmsd.dat -o rgyr_rmsd.png -t 300 -bx 100 -by 100 -lx Rgyr[nm] -ly RMSD[nm]

• -f rgyr_rmsd.dat rows represent the time-ordered conformations, the two columns correspond to the values of two order parameters, R1 and R2, for each of the conformations, e.g.
  
  
  
• -o rgyr_rmsd.png output image 
• -bx 100 -by 100 number of bins for the two order parameters
• -t 300 temperature in Kelvin for which the FES shall be produced
• -lx Rgyr[nm] -ly RMSD[nm] labels for the x and y axis 

The resulting graphics can be viewed with okular:

okular FES.png

Instructions:

1. Copy the xvg files for Rgyr and RMSD into the working directory
2. Merge the Rgyr and RMSD data into one file using following commands:

   awk ' {if ( $1 !~ /[#@]/) print $2 } ' rgyr.xvg > rgyr.dat

   awk ' {if ( $1 !~ /[#@]/) print $2 } ' rmsd_ca.xvg > rmsd_ca.dat

   paste rgyr.dat rmsd_ca.dat > rgyr_rmsd.dat

3. Calculate and plot the free energy using our Python script:

   python generateFES.py -f rgyr_rmsd.dat -o rgyr_rmsd.png -t 300 -bx 100 -by 100 -lx Rgyr[nm] -ly RMSD[nm]

4. To find the frames corresponding to an energy minimum do:

   awk ' {if ( ($1 > 0.77 ) && ($1 < 0.78 ) && ($2 > 0.75) && ($2 < 0.76)) print "Frame: ", NR-1,$1,$2}' rgyr_rmsd.dat > energy_min1.dat

   

5. Extract the frame from the trajectory for showing an example minimimum structure: 

   gmx trjconv -f md_protein.xtc -s md_protein.tpr -o md_protein_EM1_fr1485.pdb -dump 7425

   The flag -dump refers to a time, and not a frame. Thus, the time corresponding to the frame identified for a energy minimum has to be calculated. Let's assume we used a time step of 5 ps between the snapshots saved in the trajectory and that we identified frame 1485 representing the first energy minimum. We have thus obtained the corresponding time by multiplying the frame number with the time step (1485 * 5 ps = 7425 ps).
   
   Select "System" for the output (as only the protein is contained in the system).
   
   
6. Create pictures of the conformations belonging to the energy minima (using VMD).