[gmx-users] Correct method to do Cartesian PCA

[Bin Liu]

Hi All,

I am trying to reproduce the Cartesian PCA results for ALA3 in water by
Altis et al. (Fig.2 of Dihedral angle principal component analysis of
molecular dynamics simulations.   http://dx.doi.org/10.1063/1.2746330 )
It is probably the most well-known PCA analysis results. I created a system
as they suggested and got a 100 ns trajectory (after an equilibration
stage) for applying PCA.

GROMOS96 43a1
SPC water
300 K, 1 atm
PME
cutoff 1.0 nm
Box larger than they did, 2100~ waters

As all the PCA tutorials indicate, applying Cartesian PCA can be as simple
as two commands, g_covar and g_anaeig. However with all my effort of trying
all the combinations of arguments I can think of, I still can't reproduce
the literature results. Here I humbly request help from experienced users
who have applied PCA before.

The following is my approach. Suppose I have a .tpr file *npt_pr_100ns.tpr*
and a .trr file *npt_pr_100ns.trr*

*g_covar -s npt_pr_100ns.tpr -f npt_pr_100ns.trr -o eigenval.xvg -v
eigenvec.trr -av average.pdb -l covar.log*
Choose a group for the least squares fit
Group     0 (         System) has  6518 elements
Group     1 (        Protein) has    22 elements
Group     2 (      Protein-H) has    16 elements
...
Select a group: 1
Selected 1: 'Protein'

Choose a group for the covariance analysis
Group     0 (         System) has  6518 elements
Group     1 (        Protein) has    22 elements
Group     2 (      Protein-H) has    16 elements
...
Select a group: 1
Selected 1: 'Protein'

*g_anaeig -v eigenvec.trr
-f npt_pr_100ns.trr -s npt_pr_100ns.tpr -2d 2dproj.xvg -first 1 -last 2*

Select the index group that was used for the least squares fit in g_covar
Group     0 (         System) has  6518 elements
Group     1 (        Protein) has    22 elements
Group     2 (      Protein-H) has    16 elements
...
Select a group: 1
Selected 1: 'Protein'

Select an index group of 22 elements that corresponds to the eigenvectors
Group     0 (         System) has  6518 elements
Group     1 (        Protein) has    22 elements
Group     2 (      Protein-H) has    16 elements
...
Select a group: 1
Selected 1: 'Protein'

Is there anything I did incorrectly?  Another question: is -mwa argument
which turns on mass weighting for g_covar necessary?  Thank you very much.

Best Regards,

Bin