- Contemporary messages sorted: [ by date ] [ by thread ] [ by subject ] [ by author ] [ by messages with attachments ]

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

Date: Wed, 22 Jul 2015 12:47:37 -0400

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

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/amberReceived on Wed Jul 22 2015 - 10:00:02 PDT

Custom Search