Hi Eugene,
Try doing a simple R-REMD run with ntpr=1, ntwx=1, and nstlim=1.
Everything should match nicely. I think the combination of the
EPtot/mdcrd mismatch and the order of operations creates a confusing
situation when you only write/print every 500 steps or so. I don't have
time to carefully go through the code to determine what is happening, but I
believe the exchange attempt is also a half-step off of the printed
energies and so things don't match nicely unless you print every
integration step.
So my reply to your questions:
1) the potential energies in the saveene file should be taken from sander
post-processing of the frames in the reservoir
2) Unless you print energies every step during MD (not usually a good idea
for performance), only post-processing will give you the correct energy of
your mdcrd frames
--Niel
On Sun, Oct 6, 2013 at 9:58 PM, Eugene Yedvabny <eyedvabny.berkeley.edu>wrote:
> Hi Niel,
>
> Thank you for taking your time to help me out with this. Your procedure is
> the one I adopted earlier in the week after you pointed out that the mdout
> energies of running simulations do not correspond to the true energies of
> the stored trajectory frames. The issue I am running into, however, is that
> the energy printed out by Sander analysis post-REMD does not correspond to
> any of the energies I see in the REMD mdout.
>
> In my specific example. I was running reservoir REMD with ntpr=250,
> nstlim=250, and ntwx=500, meaning I was printing out energies every time I
> did a swap and storing a frame every two swaps. The excerpts I presented in
> my previous email were the exchange energies and system energies from the
> 2nd swap of my REMD run and the Sander analysis energy (maxcyc=1, so no
> minimization) of the 1st structure in the mdcrd (which should have
> corresponded to the system right at that 2nd swap).
>
> Even if the reported instantaneous energy of the system (EPtot) does not
> correspond to the real energy of the written out frame, as you pointed out,
> I would expect the reported MC-exchange energy (the EPot under "REMD
> EXCHANGE CALCULATION") to be that real energy. After all, pmemd has to
> pause propagation to perform the exchange. Correspondingly, the EPtot and
> EPot are not equivalent.
>
> What is surprising here, is that the energy reported by Sander analysis
> does not equal either of the aforementioned energies. Now, I understand
> there is different precision between Sander and PMEMD, but I find the
> difference unacceptably large. Where PMEMD reports the structure at
> -14662.13
> (EPot) or -14621.7032 (EPtot), Sander reports it at -14602.963.
>
> As such, I am a bit at a loss as to 1) which energy to use for my
> reservoirs and 2) which is the actual potential energy of the structures
> stored in my mdcrd.
>
> Perhaps one of the Amber developers can weigh in as to when these energies
> are computed and which one should be used?
>
> Thank you,
> Eugene Yedvabny
>
>
> On Sun, Oct 6, 2013 at 9:18 PM, Niel Henriksen <shireham.gmail.com> wrote:
>
> > Hi Eugene,
> >
> > I was a little confused by your email, but I'll try to clarify what I
> meant
> > a little. I was concerned that your saveene potential energies do not
> > match the frames in your reservoir. I used the following procedure to
> make
> > my reservoir
> >
> > 1) Run regular MD: write both coords and velocities to trajectory file
> > 2) Post-process trajectory file with sander to get potential energies
> > 3) Grep potential energies out of the sander output to create saveene
> file
> > (with appropriate header)
> > 4) Read trajectory file with cpptraj and write out restart files for
> > reservoir
> > 5) run R-REMD
> >
> > Good luck,
> > --Niel
> >
> >
> > On Fri, Oct 4, 2013 at 5:34 PM, Eugene Yedvabny <eyedvabny.berkeley.edu
> > >wrote:
> >
> > > Hi Niel,
> > >
> > > Thank you for the informative response. After you've mentioned that the
> > > potential energies do not correspond to the structures, I've ran
> through
> > my
> > > REMD mdout files and there are indeed energy discrepancies. For
> example,
> > > here's an excerpt from one of the replicas:
> > >
> > > ==========================REMD EXCHANGE
> > > CALCULATION==========================
> > > Exch= 2 RREMD= 2
> > > Replica Temp= 442.35 Indx= 24 Rep#= 24 EPot= -14662.13
> > > Partner Temp= 434.82 Indx= 23 Rep#= 23 EPot= -14769.94
> > > Metrop= 0.119672E+00 delta= 0.212300E+01 o_scaling= -1.00
> > > Rand= 0.863800E+00 MyScaling= -1.00 Success= F
> > > ========================END REMD EXCHANGE
> > > CALCULATION========================
> > >
> > > NSTEP = 500 TIME(PS) = 571.000 TEMP(K) = 454.47 PRESS =
> > > 0.0
> > > Etot = -10392.8625 EKtot = 4228.8407 EPtot =
> > > -14621.7032
> > > BOND = 35.7788 ANGLE = 90.6779 DIHED =
> > > 112.9705
> > > 1-4 NB = 33.1220 1-4 EEL = 91.0480 VDWAALS =
> > > 2207.4054
> > > EELEC = -17192.7058 EHBOND = 0.0000 RESTRAINT =
> > > 0.0000
> > > Ewald error estimate: 0.5519E-04
> > > TEMP0 = 442.3500 REPNUM = 24 EXCHANGE# =
> > > 2
> > >
> > >
> >
> ------------------------------------------------------------------------------
> > >
> > > My nptr and nstlim are set to the same value, so exchange calculations
> > and
> > > energy printing should correspond to the same structure, However, EPot
> > from
> > > REMD exchange doesn't match EPtot from nptr printout. Are different
> > > algorithms used to compute the two energies?
> > >
> > > I've also ran Sander's analysis mode on that replica's mdcrd, using the
> > > exact same input file but changing imin=5 and maxcyc=1, and it returned
> > the
> > > following energy for the first structure:
> > >
> > > minimizing coord set # 1
> > >
> > > Maximum number of minimization cycles reached.
> > >
> > >
> > > FINAL RESULTS
> > >
> > >
> > >
> > > NSTEP ENERGY RMS GMAX NAME
> NUMBER
> > > 1 -1.4603E+04 1.6780E+01 1.0139E+02 CZ
> 25
> > >
> > > BOND = 38.1245 ANGLE = 85.3875 DIHED =
> > > 114.1624
> > > VDWAALS = 2227.6348 EEL = -17191.9287 HBOND =
> > > 0.0000
> > > 1-4 VDW = 32.4372 1-4 EEL = 91.2190 RESTRAINT =
> > > 0.0000
> > > minimization completed, ENE= -.14602963E+05 RMS= 0.167801E+02
> > >
> > >
> > > My ntwx was set to 2x ntpr, so the first structure in the mdcrd should
> > > correspond to the 2nd exchange attempt/printout, which is what I
> provided
> > > at the top. This energy, however, doesn't matter either of the two EPot
> > > from the REMD mdout. The energy differences are significant enough
> that I
> > > can't disregard them as precision errors. So now I am completely
> confused
> > > as to which energy is actually the correct energy of my structure and
> how
> > > to go about creating my reservoirs.
> > >
> > > Any help would be much appreciated!
> > >
> > > Thank you,
> > > Eugene Yedvabny
> > >
> > >
> > > On Friday, October 4, 2013 at 1:25 PM, Niel Henriksen wrote:
> > >
> > > > Eugene,
> > > >
> > > > I don't think there is any "all-in-one" tool to do what you are
> asking.
> > > > The multi-step process is what I have used.
> > > >
> > > > Unless I misunderstand you, you seem to be using potential energies
> > (from
> > > > the mdout file) and coordinate frames from the SAME simulation.
> > > > Unfortunately, the printed potential energies do not correspond to
> the
> > > > written frames in the mdcrd file. (see
> > > > http://archive.ambermd.org/201308/0199.html ) So you probably need
> to
> > do
> > > > the Sander-in-analysis-mode step regardless of what your plans are.
> > > >
> > > > --Niel
> > > >
> > > >
> > > > On Fri, Oct 4, 2013 at 12:50 PM, Eugene Yedvabny <
> > eyedvabny.berkeley.edu(mailto:
> > > eyedvabny.berkeley.edu)>wrote:
> > > >
> > > > > Hello Amber community,
> > > > >
> > > > > Is there any functionality in cpptraj or ptraj that enables easy
> > > reservoir
> > > > > generation and energy calculation to use in pmemd REMD with -rremd
> 2?
> > > > > Currently I am generating my reservoirs by parsing through mdout
> > files
> > > and
> > > > > extracting EPtot values, followed by manually creating the header
> > with
> > > the
> > > > > random seed. This is of course doable, if a little cumbersome when
> > > > > extracting energies from several MD runs into one reservoir. But my
> > > current
> > > > > project involves slightly modifying some of the force-field
> > > parameters, and
> > > > > I would need to recalculate the energies of my reservoir structures
> > > using
> > > > > this new prmtop. It seems the only way to do so presently would be
> to
> > > run
> > > > > Sander in analysis mode on all of my trajectories and then parse
> out
> > > the
> > > > > generated mdouts into another reservoir.
> > > > >
> > > > > Since ccptraj can load many mdcrds and prmtops, it would be ideal
> to
> > > have
> > > > > a function that generates the reservoirs and savene files all in
> one
> > > go,
> > > > > akin to clusterdihedral function that is used for -rremd 3.
> > > Unfortunately I
> > > > > was not able to find a reference to any such function, or any
> > function
> > > in
> > > > > cpptraj that would generate EPtot values for all mdcrd frames. Does
> > > anyone
> > > > > have a suggestions, or is my present multi-step process the only
> way
> > > to go
> > > > > about generating the reservoir?
> > > > >
> > > > > Thank you,
> > > > > Eugene Yedvabny
> > > > >
> > > > > _______________________________________________
> > > > > AMBER mailing list
> > > > > AMBER.ambermd.org (mailto:AMBER.ambermd.org)
> > > > > http://lists.ambermd.org/mailman/listinfo/amber
> > > > >
> > > >
> > > > _______________________________________________
> > > > AMBER mailing list
> > > > AMBER.ambermd.org (mailto:AMBER.ambermd.org)
> > > > http://lists.ambermd.org/mailman/listinfo/amber
> > > >
> > > >
> > >
> > >
> > > _______________________________________________
> > > 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
> >
> _______________________________________________
> 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 Mon Oct 07 2013 - 11:30:06 PDT
