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
Received on Sun Oct 06 2013 - 22:00:03 PDT
