Re: [AMBER] PCA projection

From: Daniel Roe <daniel.r.roe.gmail.com>
Date: Tue, 19 Aug 2014 08:32:32 -0600

Hi,

First, I'd like to address some problems with your input:

matrix out correlmatrix_1000ns.dat name correlmat correl .CA

For PCA you want a covariance or mass-weighted covariance matrix, not
a correlation matrix, which is a completely different kind of
calculation. Also I note you are using only 100 frames of your
trajectory: keep in mind that you really want the number of frames
processed to be much greater than the number of coordinates (i.e. #
atoms * 3) in your atom mask, so unless your system is very small you
should use more frames.

Another issue is that currently your matrix will contain global
rotational/translational motion as well as internal motion, which is
probably not what you want. You will want to remove this by e.g.
performing an rms best-fit to an average structure prior to
calculating the matrix (and using those same coordinates for any
subsequent projections). An example script for performing PCA with
cpptraj can be found in supporting information here:
http://pubs.acs.org/doi/abs/10.1021/jp4125099

Now to address the original issue. In cpptraj Actions and Analyses are
typically not executed when the command is given; instead they are put
into an Action list and Analysis list respectively. For example, when
you give the 'matrix' command this says to cpptraj that "I will be
calculating a matrix from input trajectories". What happens when you
get to the 'projection' Action is that your 'matrix' Action and
'diagmatrix' Analysis have been initialized but not executed, so your
'evecs' data set is likewise initialized but empty.

The simple solution is to put a 'run' or 'go' command after
'diagmatrix', which tells cpptraj to execute with the current input
trajectories/Actions/Analyses, e.g.

trajin WT_1000ns.mdcrd 900 1000
matrix out correlmatrix_1000ns.dat name correlmat correl .CA
diagmatrix correlmat out evecs_1000ns.dat name evecs vecs 5
run
projection P1 modes evecs out project_1000ns.dat beg 1 end 2 .CA

When cpptraj hits 'run' it will calculate the matrix, then run the
diagonalization, so by the time the 'run' finishes the 'evecs' data
set will contain data and 'projection' can now be processed. When
'cpptraj' gets to the end of input and any Actions/Analyses/output
trajectories are set up, this implies a 'run', so 'projection' will be
executed last. However, keep in mind what I said about overall
translational/rotational motion above; you'll likely want to add to
this sequence of commands.

Let me know if you have any more questions,

-Dan

On Tue, Aug 19, 2014 at 3:52 AM, de Manzanos Guinot, Angela
<angela.de-manzanos11.imperial.ac.uk> wrote:
> Hello,
>
> I am having problems with the projection in PCA using cpptraj, I get the following error
>
> CPPTRAJ: Trajectory Analysis. V14.06b
> ___ ___ ___ ___
> | \/ | \/ | \/ |
> _|_/\_|_/\_|_/\_|_
> Reading 'WT.prmtop' as Amber Topology
> INPUT: Reading Input from file matrix_1
> [trajin WT_1000ns.mdcrd 900 1000]
> Reading 'WT_1000ns.mdcrd' as Amber NetCDF
> [matrix out correlmatrix_1000ns.dat name correlmat correl .CA]
> MATRIX: Calculating correlation matrix, output is by atom
> Printing to file correlmatrix_1000ns.dat
> Storing matrix on internal stack with name: correlmat
> Start: 1 Stop: Final frame
> Mask1: .CA
> [diagmatrix correlmat out evecs_1000ns.dat name evecs vecs 5]
> DIAGMATRIX: Diagonalizing matrix correlmat and writing modes to evecs_1000ns.dat
> Calculating 5 eigenvectors
> Storing modes with name: evecs
> [projection P1 modes evecs out project_1000ns.dat beg 1 end 2 .CA]
> Warning: 'end' 2 is greater than # modes (0); setting end to 0
> Error: 'beg' 1 out of bounds.
> Error: Could not initialize action [projection]
> 1 errors encountered reading input.
>
> My input file for cpptraj is:
>
> trajin WT_1000ns.mdcrd 900 1000
> matrix out correlmatrix_1000ns.dat name correlmat correl .CA
> diagmatrix correlmat out evecs_1000ns.dat name evecs vecs 5
> projection P1 modes evecs out project_1000ns.dat beg 1 end 2 .CA
>
> Any help will be very appreciated as I am trying all sort of combinations, but I am not able to spot the mistake.
>
> Many thanks in advance,
> Angela
> --
>
> _______________________________________________
> AMBER mailing list
> AMBER.ambermd.org
> http://lists.ambermd.org/mailman/listinfo/amber



-- 
-------------------------
Daniel R. Roe, PhD
Department of Medicinal Chemistry
University of Utah
30 South 2000 East, Room 307
Salt Lake City, UT 84112-5820
http://home.chpc.utah.edu/~cheatham/
(801) 587-9652
(801) 585-6208 (Fax)
_______________________________________________
AMBER mailing list
AMBER.ambermd.org
http://lists.ambermd.org/mailman/listinfo/amber
Received on Tue Aug 19 2014 - 08:00:02 PDT
Custom Search