Dihedral PCA

To make the index file for dihedral PCA either use mk_angndx or make an index file by hand. Now you need a suitable combination of g_angle (possibly with `-oc` or `-or`) and g_covar.

Here is some discussion on the use of dihedral PCA and a reply.

Background

The implementation of PCA in GROMACS first makes a trajectory file with reduced dimensions matching the selected angles, and then makes a fake trajectory file that contains the eigenvectors and eigenvalues. This maybe wasn't the best approach possible, but it was made to work. Steps 1 and 2 do the dimensionality reduction, steps 3 and 4 generate a kind of hacked coordinate file whose dimensionality is large enough to help gmx covar and gmx anaeig work.

Stepwise instructions

1. Use mk_angndx or text editor to generate an index file for the atoms involved in the chosen dihedral angles.

2. Extract the angles from trajectory using the index file from step 1., e.g., `dangle.ndx`, in a command like:

```g_angle -f foo.xtc -s foo.tpr -n dangle.ndx -or dangle.trr -type dihedral
```

where the output is `dangle.trr` that contains the simulation trajectory reduced to the chosen angle dimensions.

3. Make an index file, named as `covar.ndx`, which necessarily contains one group of atom 1 to integer(2*N/3), where N is the number of dihedral angles. For example, for a peptide of 5 dihedral angles, the contents of `covar.ndx` should be

```   [ foo ]
1   2   3   4
```

4. Use the index file to make a .gro file that contains integer(2*N/3) atoms. Box size and atom coordinates don't matter. For example,

```trjconv -s foo.tpr -f dangle.trr -o resized.gro -n covar.ndx -e 0
```

5. Perform the diagonalization with a command like

```g_covar -f dangle.trr -n covar.ndx -ascii -xpm -nofit -nomwa -noref -nopbc -s resized.gro
```

where the outputs are `eigenval.xvg`, `eigenvec.trr`, `covar.log`, `covar.dat`, and `covar.xpm`.

6. To get a PMF along one eigenvector use a command like

```g_anaeig -v eigenvec.trr -f dangle.trr -s resized.gro -first X -last X -proj proj-1
```

where X is the serial number of the eigenvector. To visualize the result, use

```xmgrace proj-X.xvg
```

7. To get a free energy landscape for projections along two eigenvectors, use a command like

```g_anaeig -v eigenvec.trr -f dangle.trr -noxvgr -s resized.gro -first X -last Y -2d 2dproj_X_Y.xvg
```

where X and Y are the serial numbers of the eigenvectors.

8. To transform such data into free energy landscape terms, we divide the 2D landscape with grid and count the number of conformations inside each grid. An awk script for such purpose has been posted here.  The g_sham program (specifying the `-notime` option) can also be used to generate a 2-D plot of a free energy surface.