An error with The total potential energy is -nan

User discussions
, , , ,
1

GROMACS version:
GROMACS modification: Yes/No
Here post your question
I am running an NVT calculation on a system with about 80000 atom. It runs for some seconds but then returns this error message.
Internal error (bug):

Step 100: The total potential energy is -nan, which is not finite. The LJ and
electrostatic contributions to the energy are nan and 0, respectively. A
non-finite potential energy can be caused by overlapping interactions in
bonded interactions or very large or Nan coordinate values. Usually this is
caused by a badly- or non-equilibrated initial configuration, incorrect
interactions or parameters in the topology.

The NVT setup I am using is this

title                   = NVT equilibration with annealing from 300 K to 600 K

; Run parameters
integrator              = md        ; leap-frog integrator
nsteps                  = 250000    ; 250000 * 1 fs = 250 ps
dt                      = 0.001     ; 1 fs (10^-15 seconds); 0.001
; Output control
nstxout                 = 1000      ; save coordinates every 2.0 ps
nstvout                 = 1000      ; save velocities every 2.0 ps
nstenergy               = 1000      ; save energies every 2.0 ps
nstlog                  = 1000      ; update log file every 2.0 ps
; Bond parameters
continuation            = no        ; First MD run, gen-vel=yes
constraint_algorithm    = shake      ; No constraints
constraints             = none      ; No constraints
; Nonbonded settings
cutoff-scheme           = Verlet    ; Buffered neighbor searching
;verlet-buffer-tolerance = 0.1
ns_type                 = grid      ; search neighboring grid cells
nstlist                 = 10        ; 20 fs, largely irrelevant with Verlet scheme
rcoulomb                = 1.2       ; short-range electrostatic cutoff (in nm)
rvdw                    = 1.2       ; short-range van der Waals cutoff (in nm)
DispCorr                = EnerPres  ; account for cut-off vdW scheme
; Electrostatics
coulombtype             = Cut-off   ; Cutoff for electrostatics
pme_order               = 4         ; cubic interpolation
fourierspacing          = 0.12      ; grid spacing for FFT
; Temperature coupling is on
tcoupl                  = V-rescale ; V-rescale thermostat
tc-grps                 = System    ; two coupling groups - more accurate
tau_t                   = 0.1       ; time constant, in ps
ref_t                   = 300       ; reference temperature, one for each group, in K
; Pressure coupling is off
;pcoupl                  = no
; Periodic boundary conditions
pbc                     = xy        ; 2-D PBC (only x and y)
nwall                   = 2         ; Two walls for confinement
wall-atomtype           = ljw  ljw  ; Wall atom types
wall-type               = 12-6      ; Wall potential type
wall-r-linpot           = 0.1       ; Linear potential distance from the wall
ewald-geometry          = 3dc       ; Ewald summation geometry
; Velocity generation
gen_vel                 = yes       ; Velocity generation is on

; Annealing
annealing               =  single
annealing-npoints       =  3         ;
annealing-time          =  0     100   250 ; Annealing over 250 ps with steps
annealing-temp          =  300   450   600 ; Start at 300 K and end at 600 K
2

I suggest to double check your minimization step to see if it was successful.

There are two very important factors to evaluate to determine if energy minimization was successful. The first is the potential energy (printed at the end of the EM process, and at the last lines of the log file. ). The second important feature is the maximum force, Fmax. It is possible to arrive at a reasonable Epot with Fmax > emtol. If this happens, your system may not be stable enough for simulation. Evaluate why it may be happening, and perhaps change your minimization parameters (integrator, emstep, etc).
and the last suggestion, try using the v-rescale thermostat and c-rescale barostat.

P.S. Read more here: GROMACS tutorials.

1 Like
3

One issue might be that you’re not using any constraints, and for a fully unconstrained MD to be stable, one should use a 0.5 fs timestep (yours is 1 fs).

Try one of the two options:

  1. Set your dt to 0.0005 which is 0.5 fs
  2. Set your constraints to h-bonds and constraint-algorithm to Lincs

The 2D periodic setup is also quite non-standard, make sure it’s not causing trouble by comparing with a 3D periodic setup (pbc = xyz).

1 Like