[gmx-users] NVE Temperature Drift

[Mark Abraham]

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.
>
>