Re: [AMBER] PCA_representation_in_Cpptraj

From: Juan Eiros Zamora <>
Date: Tue, 10 Mar 2015 11:37:50 +0000

Hi all,

I would like to hop on in this discussion. I have been using this
procedure to analyze my classical MD trajectories using the following
cpptraj script

parm *.prmtop
trajin *.nc

rms first

reference [ref1]
matrix covar name matrixdat .CA out covmat-ca.dat
diagmatrix matrixdat out evecs-ca.dat vecs 10

projection modes evecs-ca.dat out pca-ca.dat beg 1 end 2 .CA

readdata pca-ca.dat
runanalysis hist pca-ca.dat:2 pca-ca.dat:3 free 300 out fhist-all.CA.gnu
bins 400

I am a bit confused with the settings that have to be used for PCA
analysis. To my understanding, the first few PC are responsible for the
majority of the variance of the system. Therefore, is there a reason why
you calculate 174 PCs? I thought that maybe only calculating 10 would be
enough? Also, why do you set beg 1 end 3 in your projection function? If
one wants to plot the 1st PC vs. the 2nd PC only beg 1 end 2 is needed I

Additionally, I am not quite sure how the last line of code works. I do
not understand how changing the default bin size affects the
calculation, should I just leave it as default (erase it from my script)
if I want to analyze classical MD trajectory?

Thanks in advance,

Juan Eiros

On 09/03/15 14:07, Daniel Roe wrote:
> Hi,
> You will need to modify your original cpptraj script for generating
> PCs somewhat.
> On Mon, Mar 9, 2015 at 3:14 AM, abdennour braka
> <> wrote:
>> parm
>> trajin 1_replic_traj.trj
>> reference min.rst
>> rms reference
>> matrix covar name matrixdat .CA out covmat-ca.dat
>> analyze matrix matrixdat out evecs-ca.dat vecs 174
>> analyze matrix matrixdat name evecs-ca vecs 174
> You don't need two 'analyze matrix' commands here, just one will do
> (and 'diagmatrix' is the preferred name for this analysis now):
> diagmatrix matrixdat name evecs-ca out evecs-ca.dat vecs 174
>> analyze modes fluct out analyzemodesflu.dat name evecs-ca beg 1 end 174
>> analyze modes displ out analyzemodesdis.dat name evecs-ca beg 1 end 174
> Your 'run' command should come here, *before* you attempt your
> projection calculation since at this point nothing has been
> calculated. So the next two commands should be:
> run
> projection evecs evecs-ca out pca12-ca.dat beg 1 end 3 .CA
> Does that make sense to you? Let me know if not.
> As far as creating a weighted histogram, cpptraj can do this as well
> given that you have your boost energies in KT. Assuming you have them
> in a file called e.g. amdboost.dat in a single column, you could read
> that and your previously-calculated PC projections and calculate the
> reweighted 2D histogram like so:
> readdata amdboost.dat name BOOST
> readdata pca12-ca.dat name PC
> runanalysis hist PC:2 PC:3 amdboost BOOST bins 200 free 300 out
> hist.pca12-ca.gnu
> The '2' and '3' for PC are used since pca12-ca.dat has data for the
> first and second PC projections in columns 2 and 3. You can use 'list
> dataset' after 'readdata' to see currently loaded data. Here no number
> is needed for BOOST since there is only 1 column of data. The '.gnu'
> extension will write the resulting histogram in gnuplot-readable
> format.
> Note that this kind of reweighting is very straightforward. There are
> other ways to reweight aMD data out there that you may want to check
> out, notably here:
> Hope this helps,
> -Dan
> run
>> ###########################################
>> Thank you four any replies,
>> Abdennour
>> //
>> //
>> _______________________________________________
>> AMBER mailing list

AMBER mailing list
Received on Tue Mar 10 2015 - 05:00:04 PDT
Custom Search