[gmx-users] Error in NVT equilibration

Gonzalez Fernandez, Cristina cristina.gonzalezfdez at unican.es
Mon Sep 24 12:46:21 CEST 2018

Hi Justin,

Thank you very much for your suggestions, I ran the equilibration!

After running the equilibration, I tabulated the temperature with time, and at time=0 the temperature is 277.73 K and at time=1ps the temperature is 294.03. The initial temperature is enough similar (I think) to the considered one. However, the temperature at 1ps is lower than the stated one (298K) in the annealing settings (and therefore at 0.134 ps the temperature is lower than 294 K instead of 298K).

I copy the .mdp file in case you can detect the errors in this file, as I have been looking for the annealing settings in the GROMACS manual and in previous posts and I don't find what is wrong.

title                   = Lipid NVT equilibration
define                  = -DPOSRES  ; position restrain the lipid
; Run parameters
integrator              = md        ; leap-frog integrator
nsteps                  = 50000     ; 2 * 50000 = 100 ps
dt                      = 0.002     ; 2 fs
; Output control
nstxout                 = 500       ; save coordinates every 1.0 ps
nstvout                 = 500       ; save velocities every 1.0 ps
nstenergy               = 500       ; save energies every 1.0 ps
nstlog                  = 500       ; update log file every 1.0 ps
; Bond parameters
continuation            = no        ; first dynamics run
constraint_algorithm    = lincs     ; holonomic constraints
constraints             = all-bonds   ; bonds involving H are constrained
;shake_tol              = 0.0001
lincs_iter            = 1            ; accuracy of LINCS
lincs_order            = 4            ; also related to accuracy
; Nonbonded settings
cutoff-scheme           = Verlet    ; Buffered neighbor searching
ns_type                 = grid      ; search neighboring grid cells
nstlist                 = 10        ; 20 fs, largely irrelevant with Verlet
rlist                   = 1.4
rcoulomb                = 1.4       ; short-range electrostatic cutoff (in nm)
rvdw                    = 1.4       ; short-range van der Waals cutoff (in nm)
DispCorr                = EnerPres  ; account for cut-off vdW scheme
vdwtype                 = Cut-off
vdw_modifier            = Potential-shift-Verlet
rvdw_switch             = 1
; Electrostatics
coulombtype             = PME       ; Particle Mesh Ewald for long-range electrostatics
ewald-rtol              = 1e-6
pme_order               = 4         ; cubic interpolation
fourierspacing          = 0.16      ; grid spacing for FFT
coulomb_modifier        = Potential-shift-Verlet
; Temperature coupling is on
tcoupl                  = V-rescale             ; modified Berendsen thermostat
tc-grps                 = System   ; two coupling groups - more accurate
tau_t                   = 0.1               ; time constant, in ps
ref_t                   = 298               ; reference temperature, one for each group, in K
; Pressure coupling is off
pcoupl                  = no        ; no pressure coupling in NVT
; Periodic boundary conditions
pbc                     = xyz       ; 3-D PBC
; Velocity generation
gen_vel                 = yes       ; assign velocities from Maxwell distribution
gen_temp                = 278       ; temperature for Maxwell distribution
gen_seed                = -1        ; generate a random seed
;Simulated annealing
annealing               = single
annealing_npoints       = 2
annealing_time          = 0 0.134
annealing_temp          = 278 298

If I have understood correctly from previous posts, gen_temp must be equal to the lower annealing temperature (in this case 278K) and the ref_t must be equal to the highest annealing temperature (in this case 298 K), is it correct?

Thak you for all your help,


De: gromacs.org_gmx-users-bounces at maillist.sys.kth.se <gromacs.org_gmx-users-bounces at maillist.sys.kth.se> en nombre de Justin Lemkul <jalemkul at vt.edu>
Enviado: jueves, 20 de septiembre de 2018 16:45:21
Para: gmx-users at gromacs.org
Asunto: Re: [gmx-users] Error in NVT equilibration

On 9/19/18 6:14 AM, Gonzalez Fernandez, Cristina wrote:
> Hi Justin,
> I’ve been looking for the non-bonded setting parameters for the GROMOS force field. I have found that for this force field coulmb_modifier and vdw_modifier = Potential-shift-Verlet. I have changed this in the .mdp file, but I don’t have seen more differences, Is it correct? Is it enough to change the columb and vdw modifier option?. I obtain the same error when running the equilibration.

I haven't used the GROMOS force fields with the new Verlet options, and
it was parametrized for rather old algorithms, so I can't offer you the
latest advice on this. Electrostatic forces are typically computed with
PME, and a plain cutoff for vdW. Your cutoff values in the previous .mdp
file were also incorrect and need to be fixed; it's not just a matter of
the modifiers.

> I have tryied to use the default GMX values for dt (0,002) and nsteps (50000) because they are so different to the ones I have proposed (0,015 and 6667 respectively) and the nvt equilibration runs. My values derive from the

A time step of 15 fs is not going to be stable with an atomistic force
field. A sensible dt is usually no larger than 2 fs without playing
other tricks.

> annealing I am trying, and with the default values the annealing is not fulfilled since the temperature doesn’t start in the first annealing temperature.
> I need to increase the temperature from 278 to 298K with a rate of 0,298 K every time step for 1 ps. I checked your previous posts and use annealing for that purpose. From the rate 0,298K/step and knowing that I need to increase 20 K de temperature, I calculate the requied number of steps and from this value and 1 ps I found that dt should be 0,015 ps. From this dt, as I considered a total equilibration time of 100 ps, I calculated nsteps (6667).

Don't adjust dt, which is the simulation time step. Your simulation will
never be stable with such a value. You can't simultaneously specify a
heating rate, total change, and fixed time interval. You should consider
dt fixed and determine how much time it will take to reach the total
change in temperature at the desired rate. For your situation, to heat
by 20 K with a rate of 0.298 K/2fs, you would need to warm over a period
of about 134 ps.

> The annealing settings are:
> ;Simulated annealing
> annealing = single
> annealing_npoints = 2
> annealing_time = 0 1
> annealing_temp = 278 298
> How can I satisfy the 0,298K increase during 1 ps until 298K while running the equilibration?


annealing_time = 0 134

and you'll get the rate of heating that you want. As stated above, you
can't do it in 1 fs while also specifying the rate of heating. dt is not
a dependent variable that is subject to such manipulations.



Justin A. Lemkul, Ph.D.
Assistant Professor
Virginia Tech Department of Biochemistry

303 Engel Hall
340 West Campus Dr.
Blacksburg, VA 24061

jalemkul at vt.edu | (540) 231-3129


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.

More information about the gromacs.org_gmx-users mailing list