Carlos,
Is it possible to reproduce the sander output's targetted rmsd results
ever using ptraj? I do not think so. Here is a test case I have done.
I used the following sander input file (the same input file I gave, except
ntr is set to 0):
------ sander.in ---------------
TMD: NMR->XRAY; target rmsd 2.489
&cntrl
imin = 0, ntx = 5, nstlim = 1200000, irest=1,
dt = 0.002, ntc = 2, ntf = 2, tol = 0.000001,
tempi = 300.0, temp0 = 300.0,
scee = 1.2, cut = 99.9,
ntpr = 50, ntwx = 50, ntwr = 50,
ntb = 0, ntt = 1, tautp = 0.5, ntp = 0,
igb = 1, nscm = 0, nmropt = 0,
ntr = 0,
tgtfitmask="(:1-18)",
itgtmd = 1, tgtrmsd = 2.489, tgtmdfrc = 1.0,
tgtrmsmask="(:1-18)",
/
---------------------------------
This input file will first best-fit according to tgtfitmask, and then
calculate the target rmsd according to tgtrmsmask. Here is the first a
couple of NSTEP results of the sander output file:
----- sander.out -----------------
File Assignments:
| MDIN: /home/yildirim/tgtmd_method1_new/tgtmd.in
| MDOUT: /home/yildirim/tgtmd_method1_new/tgtmd_1.out
|INPCRD: /home/yildirim/tgtmd_method1_new/eq.rst
| PARM: /home/yildirim/tgtmd_method1_new/prmtop
|RESTRT: /home/yildirim/tgtmd_method1_new/tgtmd_1.rst
| REFC: /home/yildirim/tgtmd_method1_new/xraymin.rst
| MDVEL: mdvel
| MDEN: mden
| MDCRD: /home/yildirim/tgtmd_method1_new/tgtmd_1.traj
|MDINFO: mdinfo
|INPDIP: inpdip
|RSTDIP: rstdip
.
.
.
------------------------------------------------------------------------
NSTEP = 50 TIME(PS) = 100.100 TEMP(K) = 288.12 PRESS = 0.0
Etot = -3095.8902 EKtot = 442.5837 EPtot = -3538.4739
BOND = 150.4736 ANGLE = 335.3341 DIHED = 463.9605
1-4 NB = 147.9322 1-4 EEL = -2043.4951 VDWAALS = -323.8437
EELEC = 1673.0053 EGB = -3941.9664 RESTRAINT = 0.1256
EAMBER (non-restraint) = -3538.5995
Current RMSD from reference: 2.510
Current target RMSD: 2.489
------------------------------------------------------------------------
NSTEP = 100 TIME(PS) = 100.200 TEMP(K) = 309.66 PRESS = 0.0
Etot = -3096.7318 EKtot = 475.6731 EPtot = -3572.4050
BOND = 150.0298 ANGLE = 302.9113 DIHED = 455.4908
1-4 NB = 149.3825 1-4 EEL = -2016.5340 VDWAALS = -315.0412
EELEC = 1608.6149 EGB = -3907.2766 RESTRAINT = 0.0175
EAMBER (non-restraint) = -3572.4225
Current RMSD from reference: 2.481
Current target RMSD: 2.489
------------------------------------------------------------------------
NSTEP = 150 TIME(PS) = 100.300 TEMP(K) = 306.35 PRESS = 0.0
Etot = -3096.7824 EKtot = 470.5779 EPtot = -3567.3604
BOND = 159.0258 ANGLE = 309.7130 DIHED = 457.4631
1-4 NB = 150.1704 1-4 EEL = -2031.4703 VDWAALS = -332.7867
EELEC = 1632.4870 EGB = -3911.9729 RESTRAINT = 0.0102
EAMBER (non-restraint) = -3567.3706
Current RMSD from reference: 2.483
Current target RMSD: 2.489
------------------------------------------------------------------------
NSTEP = 200 TIME(PS) = 100.400 TEMP(K) = 303.91 PRESS = 0.0
Etot = -3097.0808 EKtot = 466.8355 EPtot = -3563.9163
BOND = 166.1752 ANGLE = 298.8125 DIHED = 465.2442
1-4 NB = 154.2386 1-4 EEL = -2046.2956 VDWAALS = -322.2201
EELEC = 1659.7919 EGB = -3939.6662 RESTRAINT = 0.0031
EAMBER (non-restraint) = -3563.9194
Current RMSD from reference: 2.486
Current target RMSD: 2.489
------------------------------------------------------------------------
I have the following ptraj input file, and its result:
-------- ptraj.in ---------------------
reference xraymin.rst
trajin tgtmd_1.traj
rms reference out rms_out_testing :1-18
go
--------------------------------------
------- rms_out_testing --------------
1.00 2.57387
2.00 2.54275
3.00 2.55708
4.00 2.57366
5.00 2.59011
6.00 2.61340
7.00 2.61528
8.00 2.60213
.
.
.
--------------------------------------
As far as I can conclude with these results, the TMD's alignment sceme is
different than ptraj's alignment sceme.
PS: You used 'restraintmask' and 'rmsmask' in your previous email. What
exactly do you mean by 'rmsmask'? Thanks.
On Sun, 23 Sep 2007, Carlos Simmerling wrote:
> Ilyas,
> no, you can't use both restraintmask and rmsmask. You
> have to choose one. We felt that with restraints, there was no point in
> doing a best-fit, but we thought people would use strong enough
> restraints to keep the system from moving much. Have you checked
> in your simulations to see how far the restrained part moves?
> You can do that by runnning ptraj and calculating rmsd without
> fitting. It is indeed "implicitly fitted", but if the restraints are very
> weak then it may not be done well.
> I suppose it could be possible in the code to do both restraints
> and a best fit, but it would introduce extra complexity since we would
> need to restore the coordinates after the best fit and rmsd, otherwise the
> restraint energy would suffer discontinuities.
>
> carlos
>
> On 9/22/07, Ilyas Yildirim <yildirim.pas.rochester.edu> wrote:
> >
> > Dear Carlos,
> >
> > Thanks again for your response. So, you say that I can use both
> > restraintmask and tgtfitmask in a targetted MD simulation, is that right?
> > I did not use tgtfitmask because in the manual it says that when ntr=1 is
> > set, the structures is implicitly fitted according to atoms specified in
> > restraintmask. Looks like that there is a mistake in this part of the
> > manual.
> >
> > I just searched the amber mailing list, and saw that someone already
> > pointed out this. I wish I read it before running the simulations. Thanks
> > again.
> >
> > Best,
> >
> > On Sat, 22 Sep 2007, Carlos Simmerling wrote:
> >
> > > no, that's not correct. when you use ntr=1, there is _no_ fitting
> > > involved. it applies positional restraints just as normal ntr=1
> > > simulations. the rmsmask specifies an rmsd calculation, but no fitting
> > > is involved.
> > >
> > > what you are doing is certainly reasonable, but you can't expect
> > > it to match something where a best fit is done.
> > > if you want tgtmd simulations with best-fitting, use the fitmask.
> > > this doesn't work with ntr=1.
> > >
> > > carlos
> > >
> > > On 9/22/07, Ilyas Yildirim <yildirim.pas.rochester.edu> wrote:
> > > >
> > > > Dear Carlos,
> > > >
> > > > Thanks for your reply, but I am not sure if what you say is right
> > > > or not compared to what is written in the TMD part of the manual.
> > > >
> > > > What I have done in this simulation is, I put positional restraints to
> > > > residues :1-2,8-11,17-18. According to the manual, the structure will
> > be
> > > > implicitly 'fitted' to this positional restraints (I assume that this
> > > > 'fitting' is best-fitting). The targetted rmsd will be, then,
> > calculated
> > > > for the residues :3-7,12-16 (without doing any fitting). That is what
> > I
> > > > wanted to do; first best-fit the structure to residues
> > :1-2,8-11,17-18,
> > > > then calculate the targetted rmsd for residues :3-7,12-16.
> > > >
> > > > The restraint weight is good enough for the system; I put this
> > positinal
> > > > restraints to the terminal base pairs so that the structure will not
> > fall
> > > > apart during the simulation.
> > > >
> > > > The way I see this simulation is as follows: With these restraints,
> > the
> > > > structure changes according to Newtonian Mechanics. At each step,
> > energy
> > > > information is calculated. Also, at each step, the structure is
> > > > best-fitted to the reference structure according to restraintmask.
> > Then,
> > > > the targetted rmsd is calculated between the reference structure and
> > > > the current MD structure for the atoms defined in tgtrmsmask (without
> > > > doing any fitting). This information is written in the sander output
> > file
> > > > as "Current RMSD from reference". I dont think that the restraint
> > weight
> > > > is the reason why the sander output's targetted rmsds values are
> > different
> > > > than the ptraj's rmsd calculation.
> > > >
> > > > The ptraj script that I have is as follows:
> > > > ----------------- ptraj.in ---------------------
> > > > reference xray.inpcrd
> > > > trajin tgtmd_1.traj
> > > > rms reference :1-2,8-11,17-18
> > > > rms reference out rmsd_out :3-7,12-16 nofit
> > > > go
> > > > ------------------------------------------------
> > > >
> > > > The problem arises in the analysis part. When I compare the ptraj
> > output
> > > > with the sander output ('Current RMSD from reference' rmsd values),
> > they
> > > > are different. Why are they different? I assumed that the alignment
> > > > protocol for targetted MD is similar to ptraj's alignment protocol.
> > > > If the methods used for the alignment is different, which result am I
> > > > going to use; the ptraj's results or the sander output file results?
> > > >
> > > > If there is something wrong with the way I look at targetted MD, I
> > > > will greatly appreciate your comments. Thanks.
> > > >
> > > > Best,
> > > >
> > > > On Sat, 22 Sep 2007, Carlos Simmerling wrote:
> > > >
> > > > > your restraint weight is very weak.
> > > > > it will calculate the non-fit rmsd for atoms in
> > > > > the rmsmask, but will do no best fit.
> > > > > when you use ptraj and do an alignment, then
> > > > > you already have something different from what
> > > > > you asked sander to do. try to use ptraj to calculate
> > > > > the rmsd without any fitting at all.
> > > > > alternatively, use a much stronger restraint
> > > > > weight and see if the agreement in rmsd improves.
> > > > > it could be something else, but I am guessing it
> > > > > is the extremely weak restraint in sander as compared
> > > > > to actual alignment in ptraj.
> > > > >
> > > > > On 9/22/07, Ilyas Yildirim <yildirim.pas.rochester.edu> wrote:
> > > > > >
> > > > > > Dear Amber Users,
> > > > > >
> > > > > > I am having trouble understanding how the RMSD is calculated in
> > > > targetted
> > > > > > MD simulations. I have the following input file:
> > > > > > ------------ tgtmd.in ------------------------
> > > > > > TMD: NMR->XRAY; target rmsd 3.267; resid 1-2,8-11,17-18 restrained
> > > > > > &cntrl
> > > > > > imin = 0, ntx = 5, nstlim = 1200000, irest=1,
> > > > > > dt = 0.002, ntc = 2, ntf = 2, tol = 0.000001,
> > > > > > tempi = 300.0, temp0 = 300.0,
> > > > > > scee = 1.2, cut = 99.9,
> > > > > > ntpr = 50, ntwx = 50, ntwr = 50,
> > > > > > ntb = 0, ntt = 1, tautp = 0.5, ntp = 0,
> > > > > > igb = 1, nscm = 0, nmropt = 0,
> > > > > > ntr = 1, restraint_wt=0.02,
> > > > > > restraintmask=":1-2,8-11,17-18",
> > > > > > itgtmd = 1, tgtrmsd = 3.267, tgtmdfrc = 1.0,
> > > > > > tgtrmsmask="(:3-7,12-16)",
> > > > > > /
> > > > > > -----------------------------------
> > > > > > What this script does is it restraints the residues 1-2,8-11,17-18
> > > > > > with a restraint weight of 0.02. While doing that, implicitly, it
> > > > aligns
> > > > > > the system with respect to restraintmask.
> > > > > >
> > > > > > Then, the rmsd is calculated for the residues 3-7,12-16. At each
> > > > > > step, "Current RMSD from reference" is printed out. I could not
> > figure
> > > > out
> > > > > > how this 'Current RMSD from reference' is calculated. When I
> > extract
> > > > that
> > > > > > particular snapshot from the trajectory and use ptraj (first
> > aligning
> > > > the
> > > > > > residues 1-2,8-11,17-18 and then calculating the rmsd of residues
> > > > > > 3-7,12-16), I get pretty different result.
> > > > > >
> > > > > > 1. Is there a mistake the way I describe/understand the above
> > input
> > > > file?
> > > > > > 2. Does TMD follow a different alignment algorithm?
> > > > > >
> > > > > > Thanks in advance for your responses.
> > > > > >
> > > > > > Best,
> > > > > >
> > > > > > --
> > > > > > Ilyas Yildirim
> > > > > > ---------------------------------------------------------------
> > > > > > = Department of Chemistry - =
> > > > > > = University of Rochester - =
> > > > > > = Rochester, NY 14627-0216 - Ph.:(585) 275 67 66 (Office) =
> > > > > > = http://www.pas.rochester.edu/~yildirim/ =
> > > > > > ---------------------------------------------------------------
> > > > > >
> > > > > >
> > > >
> > -----------------------------------------------------------------------
> > > > > > The AMBER Mail Reflector
> > > > > > To post, send mail to amber.scripps.edu
> > > > > > To unsubscribe, send "unsubscribe amber" to majordomo.scripps.edu
> > > > > >
> > > > >
> > > >
> > > > --
> > > > Ilyas Yildirim
> > > > ---------------------------------------------------------------
> > > > = Department of Chemistry - =
> > > > = University of Rochester - =
> > > > = Rochester, NY 14627-0216 - Ph.:(585) 275 67 66 (Office) =
> > > > = http://www.pas.rochester.edu/~yildirim/ =
> > > > ---------------------------------------------------------------
> > > >
> > > >
> > > >
> > -----------------------------------------------------------------------
> > > > The AMBER Mail Reflector
> > > > To post, send mail to amber.scripps.edu
> > > > To unsubscribe, send "unsubscribe amber" to majordomo.scripps.edu
> > > >
> > >
> > >
> > >
> > >
> >
> > --
> > Ilyas Yildirim
> > ---------------------------------------------------------------
> > = Department of Chemistry - =
> > = University of Rochester - =
> > = Rochester, NY 14627-0216 - Ph.:(585) 275 67 66 (Office) =
> > = http://www.pas.rochester.edu/~yildirim/ =
> > ---------------------------------------------------------------
> >
> > -----------------------------------------------------------------------
> > The AMBER Mail Reflector
> > To post, send mail to amber.scripps.edu
> > To unsubscribe, send "unsubscribe amber" to majordomo.scripps.edu
> >
>
>
>
>
--
Ilyas Yildirim
---------------------------------------------------------------
= Department of Chemistry - =
= University of Rochester - =
= Rochester, NY 14627-0216 - Ph.:(585) 275 67 66 (Office) =
= http://www.pas.rochester.edu/~yildirim/ =
---------------------------------------------------------------
-----------------------------------------------------------------------
The AMBER Mail Reflector
To post, send mail to amber.scripps.edu
To unsubscribe, send "unsubscribe amber" to majordomo.scripps.edu
Received on Mon Sep 24 2007 - 17:45:08 PDT
