Re: [AMBER] Explicit solvent pH-REMD

From: Jason Swails <jason.swails.gmail.com>
Date: Wed, 22 Jul 2015 12:47:37 -0400

On Wed, Jul 22, 2015 at 10:33 AM, whb <wwwhbpt.gmail.com> wrote:

> Dear all,
>
> I’m now setting up a system used Explicit solvent pH-REMD methods.
> Now I’ve equilibrate the system following Jason swails’s input
> (http://jswails.wikidot.com/explicit-solvent-constant-ph-md). Thanks for
> Jason’s website very much first.
>
> And in the production part, I want to relax my system in different pH,
> later
> the REMD will be used in my system. So I submitted a job with the parameter
> rem -0
>
> -O -rem 0 -i equilibrate.mdin.001 -o equilibrate.mdout.001 -c
> HdeB.
> equil.rst7 -r equilibrate.rst.001 -x equilibrate.nc.001 -inf
> equilibrate.mdinfo.001 -p HdeB.mod0.parm7 -cpin HdeB.cpin -cprestrt
> equilibrate.cpin.001 -cpout HdeB_equil.cpout.001
>
> -O -rem 0 -i equilibrate.mdin.002 -o equilibrate.mdout.002 -c
> HdeB.equil.rst7 -r equilibrate.rst.002 -x equilibrate.nc.002 -inf
> equilibrate.mdinfo.002 -p HdeB.mod0.parm7 -cpin HdeB.cpin -cprestrt
> equilibrate.cpin.002 -cpout HdeB_equil.cpout.002
>
> …
>
> My system contains 158AA, and the box size is about 60A wide.
> And this is one of my input
>
> Explicit solvent constant pH MD
>
> &cntrl
>
> imin=0, irest=1, ntx=5, ntxo=2,
>
> ntpr=1000, ntwx=1000, nstlim=1000000,
>
> dt=0.002, ntt=3, tempi=300,
>
> temp0=300, gamma_ln=5.0, ig=-1,
>
> ntc=2, ntf=2, cut=10, iwrap=1,
>
> ioutfm=1, icnstph=2, ntcnstph=100,
>
> solvph=3, ntrelax=100, saltcon=0.1,
>
> /
>
>
>
> But I find that my result is a bit strange,
>
> NSTEP = 25000 TIME(PS) = 4450.000 TEMP(K) = 300.53 PRESS =
> 0.0
>
> Etot = -42092.0788 EKtot = 11035.7300 EPtot =
> -53127.8088
>
> BOND = 483.2164 ANGLE = 1414.5800 DIHED =
> 1814.3705
>
> 1-4 NB = 568.2945 1-4 EEL = 7579.1143 VDWAALS =
> 5700.8533
>
> EELEC = -70688.2377 EHBOND = 0.0000 RESTRAINT =
> 0.0000
>
> Ewald error estimate: 0.8623E-02
>
>
> ----------------------------------------------------------------------------
> --
>
>
>
> wrapping first mol.: -40.77263 28.83061 49.93608
>
> wrapping first mol.: -40.77263 28.83061 49.93608
>
>
>
> NSTEP = 26000 TIME(PS) = 4452.000 TEMP(K) = 263.77 PRESS =
> 0.0
>
> Etot = -43266.8387 EKtot = 9686.1128 EPtot =
> -52952.9515
>
> BOND = 475.0383 ANGLE = 1457.7215 DIHED =
> 1807.6935
>
> 1-4 NB = 585.0764 1-4 EEL = 7580.6866 VDWAALS =
> 5623.4670
>
> EELEC = -70482.6348 EHBOND = 0.0000 RESTRAINT =
> 0.0000
>
> Ewald error estimate: 0.8918E-02
>
>
>
> The Ewald error estimate is roughly too high and the temperature is
> unstable
> here.
>

​The Ewald error estimate seems like it was fixed by update 4 for Amber 14
(see http://ambermd.org/bugfixes14.html). Are you sure you are running
with all of the latest updates applied? To be honest, I never look at the
Ewald error estimate, so it may be that the Ewald error estimate
calculation is wrong in explicit solvent CpHMD simulations. When I
implemented the method originally, I ran energy calculations by hand to
make sure that the total energies matched before and after protonation
state changes (which they did). I never checked the Ewald error estimate.

But it might be correct -- keep in mind that you are often running with a
net charge on the system, and while that's taken care of by the
neutralizing plasma it may still have an effect on the estimate of the
error of using PME compared to pure Ewald. But I don't think this is any
cause for concern.

I can verify that the temperature jumps appear in some of my output files
as well. These correspond to time steps with successful protonation state
changes, after which the solute is frozen and the solvent relaxes. So it's
clear that in the calculation of the temperature, the velocities of the
solute appear to be zero, which is pushing down the calculated value. The
velocities are properly restored for the solute following relaxation
dynamics, though, and velocities are re-estimated from position changes
*after* dynamics has been propagated (you can see where this is done in
constantph.F90 and runmd.F90 in the pmemd source code). So I think this is
a simple reporting issue, but I'll try to investigate further when I have
the time.

HTH,
Jason

-- 
Jason M. Swails
BioMaPS,
Rutgers University
Postdoctoral Researcher
_______________________________________________
AMBER mailing list
AMBER.ambermd.org
http://lists.ambermd.org/mailman/listinfo/amber
Received on Wed Jul 22 2015 - 10:00:02 PDT
Custom Search