Principal Component Analysis and Free Energy Landscape Mapping Using GROMACS
A free-energy landscape in this workflow is a histogram of the trajectory in the plane of the first two principal components. PCA finds those axes. SHAM turns the occupancy into a free energy. The frame you extract is only as meaningful as the components you kept.
The trajectory names come from the 100 ns run md_0_100. One projection file is labeled lovastatin (pc15_lovastatina.xvg). That is the example name, not a requirement.
520 is the dump time, in picoseconds, from the example minimum. Replace it after you read fel.dat.Covariance
gmx covar wants a .tpr. If a run was pointed at md_0_100.gro and failed, use the .tpr. When prompted, select protein and ligand.
gmx covar -s md_0_100.tpr -f md_0_100.xtc -o eigenvalues.xvg -v eigenvectors.trr -xpma covar.xpm
How many components
From eigenvalues.xvg:
- Total variance is the sum of the eigenvalues.
- Explained variance of one component is that eigenvalue divided by the total, times 100.
- Cumulative variance adds those fractions until you pass the cutoff.
Keep the smallest set of components whose cumulative variance is more than 50%.
Project PC1 through PC5
-proj writes a 1D time series for each selected eigenvector. -2d writes a combined projection of the first five. -tu ns puts the time axis in nanoseconds. The lovastatin filenames are the example.
gmx anaeig -f md_0_100.xtc -s md_0_100.tpr -v eigenvectors.trr -first 1 -last 5 \
-proj pc15_lovastatina.xvg -2d project2d_s100_lovastatina.xvg -tu ns
The PC1–PC2 plane
Project the first component, then the second, and paste them into one three-column file. The awk fields assume the usual .xvg layout of this paste: time and value from PC1, then the value from PC2.
gmx anaeig -f md_0_100.xtc -s md_0_100.tpr -v eigenvectors.trr -last 1 -proj pc1.xvg
gmx anaeig -f md_0_100.xtc -s md_0_100.tpr -v eigenvectors.trr -first 2 -last 2 -proj pc2.xvg
paste pc1.xvg pc2.xvg | awk '{print $1, $2, $4}' > PC1PC2.xvg
Free energy with SHAM
gmx sham -f PC1PC2.xvg -ls FES.xpm
xpm2txt.py turns the XPM matrix into three columns: X, Y, and energy. The invocation recorded here is Python 2.7.
python2.7 xpm2txt.py -f FES.xpm -o fel.dat
The minimum
Plot fel.dat, read the time of the lowest point, and dump that frame. The example time is 520 ps. It is not a universal minimum.
gmx trjconv -s md_0_100.tpr -f md_0_100.xtc -o min_energy.pdb -dump 520
That PDB is the conformation that occupied the bottom of this PC1–PC2 histogram. It is not, by itself, proof that the simulation visited every basin that matters. The 50% cutoff and the choice to stop at two components are the assumptions sitting under the picture.
