[gmx-users] NVE Temperature Drift

Tue Aug 26 15:29:23 CEST 2014

I don't have much idea about the correct numbers to set. Mostly, I copy the
old lysozyme tutorial (
http://www.bevanlab.biochem.vt.edu/Pages/Personal/justin/gmx-tutorials/lysozyme_old/03_solvate.html),
deleted the temperature and pressure couplings, and set the numbers bigger.
I thought that bigger values won't make the simulation wrong, while smaller
values might. Then, I ran it once, and found that the temperature seem to
drop at 120,000 time step. So, in my second attempt, I only ran for 32,768
time steps, changed the timestep, and then add the npt equilibration phase.

Gromacs 4.6.6

Water box 1.5 nm away from the protein:
editconf -f $angle.complex.nowat.amber99sb.gro -o
$angle.complex.nowat.amber99sb.newbox.gro -c -d 1.5 -bt dodecahedron

While the old lysozyme tutorial (
http://www.bevanlab.biochem.vt.edu/Pages/Personal/justin/gmx-tutorials/lysozyme_old/03_solvate.html)
used a -d 1.0 waterbox, I worry that the NVE simulation would not be stable
enough, so I put more water.

I made up these three numbers:
rlist        = 1.8        ; short-range neighborlist cutoff (in nm)
rcoulomb    = 1.5        ; short-range electrostatic cutoff (in nm)
rvdw        = 1.5        ; short-range van der Waals cutoff (in nm)

And set the cut offs larger than those in the tutorial, which used OPLS
forcefield. What should be the correct numbers for Amber99SB? Where can I
find those?

I searched for NVE simulations on the mailing list, and one of the email
mentioned verlet-buffer-drift.
I tried it, but then gromacs say no verlet-buffer-drift without temperature
coupling. So, I dropped it (verlet-buffer-drift=-1).

Another email suggested reading the gromacs 4 paper on JCTC for a NVE
protocol. Is "GROMACS 4: Algorithms for Highly Efficient, Load-Balanced,
and Scalable Molecular Simulation" the gromacs 4 paper?

I ran NVE simulation because several lecture notes online say some
thermostats can destroy momentum transfer or lead to non-canonical
ensembles, and several emails on maillist mentioned NVE for transport
properties. The thing that I try to look at is not an average. So, I guess
that is similar to transport properties, and NPT is not good.

On Tue, Aug 26, 2014 at 4:38 AM, Mark Abraham <mark.j.abraham at gmail.com>
wrote:

> On Mon, Aug 25, 2014 at 9:18 PM, Johnny Lu <johnny.lu128 at gmail.com> wrote:
>
> > Hi. I'm trying to run NVE on a
> > ~170 amino acid protein
> > with a harmonic dihedral angle restraint of 4184 kJ/mol rad^2,
> >
>
> What happens without this?
>
> in a rhombic dodecahedron water box that is 15 angstroms away from the
> > protein
> >
>
> What's 1.5nm away from the protein?
>
>
> > using single precision gromacs,
> >
>
> What version?
>
>
> > and 0.5 fs time step.
> >
> > Since I only want to get many short stretches of NVE simulations, I
> > alternate npt and nve:
> >
> > 0. Steepest descent minimization with emtol = 1000.0.
> > 1. [NVT equilibration, 2fs timestep 100,000 frames]
> > 2. [NPT equilibration, 2fs timestep, 100,000 frames]
> > 3. [NPT equilibration, 2fs timestep, 5,000 frames, Berendsen thermostat
> at
> > 300K, position restraint?]
> > 4. [NVE md,0.5 fs timestep, 35,000 frames, no thermostat]
> > 5. go to step 3.
> >
> > (all of these steps have the dihedral angle restraint)
> >
> > In the NVE step, I get about 0.1% total energy decrease from start to
> > end(figure). (the total energy keep dropping).
> > But, the temperature still fluctuate in a range of ~ 5 Kelvin (figure
> > <http://oi60.tinypic.com/sni5at.jpg>).
> >
>
> Attachments, etc. are not permitted on the list. Please share a link from a
> file-sharing service if you want to illustrate.
>
>
> > Is that normal? Does that happen in real world? Where does this
> > temperature fluctuation come from? How to reduce that? Does the restraint
> > have anything to do with this?
> > [image: Inline image 2]
> >
> > [image: Inline image 1]
> >
> > Below is the mdp file that I was using:
> >
> > ; Amber99SB MD
> > ; Run parameters
> > integrator    = md        ; leap-frog integrator
> > nsteps        = 32768        ; 2 * 100000 = 200 ps, 1 ns
> > dt        = 0.0005        ; 1 fs
> > ; Output control
> > nstfout        = 2
> > nstxtcout    = 2        ; xtc compressed trajectory output every 2 fs
> > nstlog        = 1000        ; update log file every 2 ps
> > xtc-grps=protein;
> >
> > continuation    = yes        ; Restarting after NPT
> > constraint_algorithm = Lincs    ; holonomic constraints
> > constraints    = h-bonds    ; all bonds (even heavy atom-H bonds)
> > constrained
> >
>
> Normally you'd use 1fs with h-bonds or 0.5fs without constraints, but this
> should not be a problem.
>
>
> > lincs_iter    = 2        ; accuracy of LINCS
> > lincs_order    = 4        ; also related to accuracy
> > ; Neighborsearching
> > ns_type        = grid        ; search neighboring grid cells
> > nstlist        = 10        ; 10 fs
> > rlist        = 1.8        ; short-range neighborlist cutoff (in nm)
> > rcoulomb    = 1.5        ; short-range electrostatic cutoff (in nm)
> > rvdw        = 1.5        ; short-range van der Waals cutoff (in nm)
> >
>
> Is this a published protocol that is known to work well? It surely isn't
> for AMBER99SB!
>
>
> > ; Electrostatics
> >
> > coulombtype    = PME        ; Particle Mesh Ewald for long-range
> > electrostatics
> > pme_order    = 6        ; cubic interpolation
> >
>
> This is overkill (use 4), but probably not a problem.
>
> Mark
>
>
> >
> > fourierspacing    = 0.12        ; grid spacing for FFT
> > optimize_fft    = yes
> >
> > tcoupl                   = no
> > pcoupl                   = no
> >
> > compressibility = 4.5e-5    ; isothermal compressibility of water, bar^-1
> > ; Periodic boundary conditions
> > pbc        = xyz        ; 3-D PBC
> > ; Dispersion correction
> > DispCorr    = EnerPres    ; account for cut-off vdW scheme
> > ; Velocity generation
> > gen_vel        = no        ; Velocity generation is off
> > verlet-buffer-drift=-1
> > cutoff-scheme = Verlet;
> >
> > Thanks again.
> >
> > --
> > Gromacs Users mailing list
> >
> > * Please search the archive at
> > http://www.gromacs.org/Support/Mailing_Lists/GMX-Users_List before
> > posting!
> >
> > * Can't post? Read http://www.gromacs.org/Support/Mailing_Lists
> >
> > * For (un)subscribe requests visit
> > https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users or
> > send a mail to gmx-users-request at gromacs.org.
> >
> >
> --
> Gromacs Users mailing list
>
> * Please search the archive at
> http://www.gromacs.org/Support/Mailing_Lists/GMX-Users_List before
> posting!
>
> * Can't post? Read http://www.gromacs.org/Support/Mailing_Lists
>
> * For (un)subscribe requests visit
> https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users or
> send a mail to gmx-users-request at gromacs.org.
>