[gmx-users] Energy conservation in NVE simulations

[Michael K. Gilson]

Dear All,

Below are more details on our efforts to maximize energy conservation in 
NVE runs with gromacs, while maintaining performance.  I'd appreciate 
advice on the following:

1) Is this level of drift about "as good as it gets"?
2) Or, are there ways to make it even lower?
3) Are there settings we can dial back to improve speed without 
substantially increasing drift?

Many thanks,
Mike Gilson

    System: The system comprises 68679 waters, one protein with 6810
    atoms, 8 sodium ions; total 75497 atoms

    The input file is below. We have deliberately turned off COM
    restraints because we are interested in how the protein moves.

    As noted before, for lincs-iter=4, the energy slope is approximately
    -0.01976 kJ/mol/ps; for lincs-iter=5, the slope is approximately
    0.00991 kJ/mol/ps. Graphs of the energy for these respective
    conditions are provided here:
    https://drive.google.com/open?id=0B0zEqb9pykUWRHlqSEJTa2lwYkk
    https://drive.google.com/open?id=0B0zEqb9pykUWMWl0OEFtNGdVbzg



    MDP:

    ;DON'T POSITION RESTRAIN THE PROTEIN: NO -DPOSRES

    ;run parameters
    integrator      = md              ; leap-frog integrator
    nsteps          = 10000000      ; 0.001ps * 10000000 = 10 ns
    dt                = 0.001            ; 0.001ps smaller time step for NVE
    comm-mode	= None       ; let center of mass move, calculate diffusion
    constant

    ; Output control
    nstxout         = 20000          ; save coordinates every 0.5 ns
    nstvout         = 20000          ; save velocities every 0.5 ns
    nstenergy     = 20000         ; save energies every 0.5 ns
    nstlog          = 20000          ; update log file every 0.5 ns
    energygrps	= protein non-protein

    ; Bond parameters
    continuation                = yes          ; starting after nvt.mdp
    constraint_algorithm    = lincs         ; holonomic constraints
    constraints                 = all-bonds   ; all bonds are constrained
    lincs_iter                    = 4             ; NVE requires high lincs-iter
    lincs_order                 = 4             ; accuracy

    ; Neighborsearching
    cutoff-scheme           = Verlet
    ns_type         = grid          ; search neighboring grid cells
    nstlist           = 10            ; Frequency to update the neighbor list
    and long range forces
    rcoulomb              = 1.2           ; short-range electrostatic cutoff
    (in nm)
    rvdw                    = 1.2           ; short-range van der Waals cutoff
    (in nm)
    verlet-buffer-tolerance = 0.0000005	; kJ/mol/ps per particle
    pbc                     = xyz           ; periodic boundary conditions
    DispCorr              = EnerPres      ; account for cut-off vdW scheme

    ; Electrostatics
    coulombtype             = PME           ; Particle Mesh Ewald for
    long-range electrostatics
    pme_order               = 4             ; cubic interpolation
    fourierspacing           = 0.12          ; grid spacing for FFT


    ; Pressure coupling is off
    pcoupl          = no            ; fixed box size gives constant V

    ; Temperature coupling is off + no simluated annealing
    tcoupl          = no             ; NVE conditions

    ; Velocity generation
    gen_vel         = no           ; don't assign velocities: keep from
    thermalization