Re: [AMBER] Principal Component/Kullback-Leibler Analysis

From: Jonathan Gough <jonathan.d.gough.gmail.com>
Date: Tue, 12 Aug 2014 13:47:19 -0400

Thanks Daniel, that works.

I have some other questions about PCA.

In the supplemental information from your paper you (I think) you projected
the diagmatrix onto different frames from different time points of your
overall trajectory.

1. I would like to simply calculate the the first 20 PC and plot their
respective histograms for the entire trajectory (Is this appropriate or am
i supposed to break things up as you did in the paper?). Below is the
cpptraj script.

I tried doing this but I get an error saying that:

"Warning: Set 'PC-1' contains no data" and indeed no histograms are
produced... where am I going wrong?

2. We would like to visualize the trajectory of the PC. Using the modes
(analyze modes) command, (section 28.14.3) we can spit out a trajectory
(trajout). That being said... there isn't a prmtop for the output file.

Is there a way to generate a prmtop for said trajectory?

If it's done as a pdb, It does spit out something that is viewable, but
there aren't any TER cards or frame end statements between the different
models, therefore whatever is "visualized" is a horrendous mess.
What is the best way to output models that can be visualized?
Is there a way to output specific frames with all atoms?

Any help would be appreciated!
Jonathan

cpptraj script:
parm nowat.dimer.prmtop
trajin prod01.nc
rms first :1-164&!.H=
average dimerAVG.rst restart
createcrd crd1
run
reference dimerAVG.rst.1 [avg]
crdaction crd1 rms ref [avg] :1-164&!.H=
crdaction crd1 matrix covar :1-164&!.H= name dimerCovar
runanalysis diagmatrix dimerCovar out evecs1.dat vecs 20 name evecs1.dat
crdaction crd1 projection modes evecs1.dat beg 1 end 20 :1-164&!.H= myproj
hist myproj:1 bins 200 out pca1.hist.agr normint name PC-1
hist myproj:2 bins 200 out pca1.hist.agr normint name PC-2
hist myproj:3 bins 200 out pca1.hist.agr normint name PC-3
hist myproj:4 bins 200 out pca1.hist.agr normint name PC-4
hist myproj:5 bins 200 out pca1.hist.agr normint name PC-5
hist myproj:6 bins 200 out pca1.hist.agr normint name PC-6
hist myproj:7 bins 200 out pca1.hist.agr normint name PC-7
runanalysis modes name evecs1.dat trajout pc1.pdb trajoutfmt pdb
trajoutmask :1-164 tmode 1
run
quit




On Fri, Aug 8, 2014 at 12:41 PM, Daniel Roe <daniel.r.roe.gmail.com> wrote:

> Hi,
>
> Add 'name evecs.dat' to the 'runanalysis diagmatrix' line. This will
> name the modes data set created by diagmatrix to 'evecs.dat' so that
> the subsequent 'projection' commands can find it, e.g.:
>
> runanalysis diagmatrix dimerCovar out evecs.dat vecs 20 name evecs.dat
>
> Let me know if that doesn't work for you,
>
> -Dan
>
> On Fri, Aug 8, 2014 at 10:17 AM, Jonathan Gough
> <jonathan.d.gough.gmail.com> wrote:
> > Hi,
> >
> > I am attempting to use the cpptraj script described in the paper:
> >
> > "Evaluation of Enhanced Sampling Provided by Accelerated
> > Molecular Dynamics with Hamiltonian Replica Exchange Methods
> > Daniel R. Roe, Christina Bergonzo, and Thomas E. Cheatham, III"
> >
> > I am essentially using the same script as described in the supplemental
> > info, modified for my system.
> >
> > I believe the script is getting stuck at the command:
> > crdaction crd1 projection T1 modes evecs.dat beg 1 end 20
> :1-164&!.N,CA,C= \
> > crdframes 568,10566 out T1.dat
> >
> > Stating that " Error: Modes should be read in prior to this command with
> > 'readdata' "
> >
> > That being said, I am not sure where said modes are to be called from. I
> > thought those were generated by the ccptraj script.
> >
> > ANy insight would be appreciated.
> >
> > note: running
> > AmberTools version 14.07
> > Amber version 14.03
> >
> > here is the script I am using:
> >
> > parm dimer.prmtop
> > trajin prod01-dimer.mdcrd
> > trajin prod02-dimer.mdcrd
> > trajin prod03-dimer.mdcrd
> > trajin prod04-dimer.mdcrd
> > autoimage
> > rms first :1-164&!.N,CA,C=
> > average dimerAVG.rst restart
> > createcrd crd1
> > run
> > reference dimerAVG.rst.1 [avg]
> > crdaction crd1 rms ref [avg] :1-164&!.N,CA,C=
> > crdaction crd1 matrix covar :1-164&!.N,CA,C= name dimerCovar
> > runanalysis diagmatrix dimerCovar out evecs.dat vecs 20
> > crdaction crd1 projection T1 modes evecs.dat beg 1 end 20
> :1-164&!.N,CA,C= \
> > crdframes 568,10566 out T1.dat
> > crdaction crd1 projection T2 modes evecs.dat beg 1 end 20
> :1-164&!.N,CA,C= \
> > crdframes 10568,20566 out T2.dat
> > crdaction crd1 projection T3 modes evecs.dat beg 1 end 20
> :1-164&!.N,CA,C= \
> > crdframes 20568,30566 out T3.dat
> > crdaction crd1 projection T4 modes evecs.dat beg 1 end 20
> :1-164&!.N,CA,C= \
> > crdframes 30568,40566 out T4.dat
> > kde T1:1 kldiv T2:1 klout KL-PC.agr bins 400 name MD-1
> > kde T1:2 kldiv T2:2 klout KL-PC.agr bins 400 name MD-2
> > kde T1:3 kldiv T2:3 klout KL-PC.agr bins 400 name MD-3
> > kde T1:4 kldiv T2:4 klout KL-PC.agr bins 400 name MD-4
> > kde T1:1 out kde-PC.agr bins 400 name KDE1-1
> > kde T1:2 out kde-PC.agr bins 400 name KDE1-2
> > hist T1:1,*,*,*,200 out pca.hist.agr normint name HIST1-1
> > hist T2:1,*,*,*,200 out pca.hist.agr normint name HIST1-2
> > run
> > quit
> > _______________________________________________
> > 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
>
_______________________________________________
AMBER mailing list
AMBER.ambermd.org
http://lists.ambermd.org/mailman/listinfo/amber
Received on Tue Aug 12 2014 - 11:00:03 PDT
Custom Search