Dihedral PCA

    Table of contents
    1. 1. Background
    2. 2. Stepwise instructions

    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.


    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.

    Page last modified 16:24, 28 Feb 2017 by mabraham