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

From: Daniel Roe <daniel.r.roe.gmail.com>
Date: Tue, 12 Aug 2014 13:06:30 -0600

Hi,

On Tue, Aug 12, 2014 at 11:47 AM, Jonathan Gough
<jonathan.d.gough.gmail.com> wrote:
> 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.

That's fine, although you need to keep your mask consistent, so in the
'modes' command change to 'trajoutmask :1-164&!.H='. The only reason I
calculated the eigenvectors using the combined trajectory was to
ensure that all principal components would be the same for each
individual trajectory. By subsequently calculating the projections on
the individual trajectories you can then essentially see that
trajectory's contribution to the PC.

> 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?

This is strange. Are you sure you are using an up-to-date version of
cpptraj? Older versions of the 'hist' analysis did not use the 'name'
keyword. I would need to see your entire cpptraj output to debug
further.

> 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.

You can easily create a topology to match your trajectory using
'parmstrip'/'parmwrite', e.g.

parm nowat.dimer.prmtop
parmstrip !(:1-164&!.H=)
parmwrite out modes.parm7


Hope this helps,

-Dan

>
> 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



-- 
-------------------------
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 12 2014 - 12:30:02 PDT
Custom Search