[gmx-users] help me in understanding PCA

[Upendra N]

Dear gromacs users,

I am trying to understand covariance matrix produced in the covariance
analysis.

For understanding purpose, I have simulated 4 residue peptide (10ps, 6
frames ) and constructed covariance matrix for c-alpha atoms.

matrix is in 48x3 form (but it is supposed to be 16x16 as per the formula
in the manual

, where i,j runs from 1 to 3N)


I would like to understand the meaning of elements in the matrix.

for example: 1st row 2nd column may be the covariance of x and y components
of 1st c-alpha atom.

Please enlighten me whether I am thinking in a right way. I also tried to
reproduce* 1st row 2 column element using manual calculation from spread
sheet but i didn't get the value which is in covariance matrix.

Kindly help me in understanding the meaning of elements in matrix. Also how
to visualize in 3Nx3N form, since covar.dat is in (3NN)X(3) FORM

* I have extracted co-ordinates of individual frames using gmx dump.
 Then I used spread sheet to reproduce the covariance element.

*****************First 3 rows of covar.dat*****************

0.000231868   -0.000110616    6.46952e-05
-0.000275112   -3.07266e-05   -1.78349e-05
-0.000104521   0.000100497   -7.23912e-05
0.000147767    4.08451e-05    2.55278e-05
****************************************************************
***************** Variance caluculation for 1st element************
*Frames* *X-coordinates of 1st atom* *X-avg* *<(x-<x>)^2>*
1 1.75345 1.749638333 1.45288053188883E-05
2 1.821 1.749638333 0.0050924875
3 1.76899 1.749638333 0.000374487
4 1.67863 1.749638333 0.0050421834
5 1.75432 1.749638333 0.000021918
6 1.72144 1.749638333 0.000795146
*sum * 10.49783
0.0113407507
*average* 1.7496383333
0.0018901251
****************************************************************************

******************* Trajectory frames***************************************
prod_ca.trr frame 0:
   natoms=         4  step=         0  time=0.0000000e+00  lambda=         0
   box (3x3):
      box[    0]={ 3.35198e+00,  0.00000e+00,  0.00000e+00}
      box[    1]={ 0.00000e+00,  3.35198e+00,  0.00000e+00}
      box[    2]={ 0.00000e+00,  0.00000e+00,  3.35198e+00}
   x (4x3):
      x[    0]={ 1.75345e+00,  2.11837e+00,  1.48333e+00}
      x[    1]={ 1.66725e+00,  1.74678e+00,  1.55009e+00}
      x[    2]={ 1.87095e+00,  1.45219e+00,  1.68827e+00}
      x[    3]={ 1.83285e+00,  1.07340e+00,  1.72068e+00}
prod_ca.trr frame 1:
   natoms=         4  step=      1000  time=2.0000000e+00  lambda=         0
   box (3x3):
      box[    0]={ 3.35023e+00,  0.00000e+00,  0.00000e+00}
      box[    1]={ 0.00000e+00,  3.35023e+00,  0.00000e+00}
      box[    2]={ 0.00000e+00,  0.00000e+00,  3.35023e+00}
   x (4x3):
      x[    0]={ 1.82100e+00,  2.09153e+00,  1.47850e+00}
      x[    1]={ 1.66768e+00,  1.74522e+00,  1.52569e+00}
      x[    2]={ 1.81584e+00,  1.50417e+00,  1.77563e+00}
      x[    3]={ 1.75791e+00,  1.12533e+00,  1.80000e+00}
prod_ca.trr frame 2:
   natoms=         4  step=      2000  time=4.0000000e+00  lambda=         0
   box (3x3):
      box[    0]={ 3.34566e+00,  0.00000e+00,  0.00000e+00}
      box[    1]={ 0.00000e+00,  3.34566e+00,  0.00000e+00}
      box[    2]={ 0.00000e+00,  0.00000e+00,  3.34566e+00}
   x (4x3):
      x[    0]={ 1.76899e+00,  2.08417e+00,  1.42760e+00}
      x[    1]={ 1.60314e+00,  1.73896e+00,  1.50322e+00}
      x[    2]={ 1.74817e+00,  1.49365e+00,  1.75878e+00}
      x[    3]={ 1.60434e+00,  1.14668e+00,  1.80072e+00}
prod_ca.trr frame 3:
   natoms=         4  step=      3000  time=6.0000000e+00  lambda=         0
   box (3x3):
      box[    0]={ 3.34446e+00,  0.00000e+00,  0.00000e+00}
      box[    1]={ 0.00000e+00,  3.34446e+00,  0.00000e+00}
      box[    2]={ 0.00000e+00,  0.00000e+00,  3.34446e+00}
   x (4x3):
      x[    0]={ 1.67863e+00,  2.06560e+00,  1.49634e+00}
      x[    1]={ 1.54220e+00,  1.72828e+00,  1.61010e+00}
      x[    2]={ 1.69418e+00,  1.44080e+00,  1.80575e+00}
      x[    3]={ 1.50439e+00,  1.11635e+00,  1.84362e+00}
prod_ca.trr frame 4:
   natoms=         4  step=      4000  time=8.0000000e+00  lambda=         0
   box (3x3):
      box[    0]={ 3.34345e+00,  0.00000e+00,  0.00000e+00}
      box[    1]={ 0.00000e+00,  3.34345e+00,  0.00000e+00}
      box[    2]={ 0.00000e+00,  0.00000e+00,  3.34345e+00}
   x (4x3):
      x[    0]={ 1.75432e+00,  2.06958e+00,  1.33725e+00}
      x[    1]={ 1.63294e+00,  1.74710e+00,  1.51241e+00}
      x[    2]={ 1.77183e+00,  1.43452e+00,  1.68678e+00}
      x[    3]={ 1.57508e+00,  1.13509e+00,  1.82323e+00}
prod_ca.trr frame 5:
   natoms=         4  step=      5000  time=1.0000000e+01  lambda=         0
   box (3x3):
      box[    0]={ 3.34616e+00,  0.00000e+00,  0.00000e+00}
      box[    1]={ 0.00000e+00,  3.34616e+00,  0.00000e+00}
      box[    2]={ 0.00000e+00,  0.00000e+00,  3.34616e+00}
   x (4x3):
      x[    0]={ 1.72144e+00,  1.99902e+00,  1.42431e+00}
      x[    1]={ 1.60459e+00,  1.70278e+00,  1.63203e+00}
      x[    2]={ 1.71879e+00,  1.35945e+00,  1.75611e+00}
      x[    3]={ 1.51826e+00,  1.05315e+00,  1.85289e+00}

************************************************************************




With Regards
Upendra N
Research Scholar
Department of Studies in Physics,Manasagangotri
University of Mysore, Mysuru