Gmx mdrun is stuck at step=0

User discussions
,
1

GROMACS version: 2020.7
GROMACS modification: No
Here post your question
Hello everyone I’ve attempted to make a solvent box for acetonitrile in 45C; in charmm36 forcefield but my NVT calculation seams stuck at step=0
I have 6000 atoms in my system and I successfully performed a 50000 step energy minimization.
*Bellow is the .mdp file that I used
integrator = md
nsteps = 100000
dt = 0.002

nstxout = 500
nstvout = 500
nstenergy = 500
nstlog = 500

rlist = 1.0
coulombtype = cut-off
cutoff-scheme = verlet
rcoulomb = 1.3
rvdw = 1.3

tcoupl = V-rescale
tc-grps = system
tau_t = 0.1
ref_t = ?

pcoupl = no

pbc =xyz

gen-vel = yes
gen-temp = ?

I’ve already tried ?=298 and 318 but no luck so far
I appreciate any help in advance.

2

Could you copy-paste all the messages you got in either the terminal or output file (so any notes, warnings, errors, etc.) when you tried running NVT?

3

50000 steps of energy minimization sound like a lot. Maybe your initial configuration has some issues.
What is the largest force at the end of the energy minimization?

4

this is from the .log file for Energy minimization
Step Time
2773 2773.00000

Energies (kJ/mol)
Angle Proper Dih. LJ-14 Coulomb-14 LJ (SR)
1.74628e+03 8.40844e+00 -8.17765e+02 -5.48984e+04 -8.61113e+03
Coulomb (SR) Potential Pressure (bar) Constr. rmsd
-5.33882e+03 -6.79114e+04 1.07928e+04 9.88824e-07

       Step           Time
       2774     2774.00000

       Step           Time
       2775     2775.00000

       Step           Time
       2776     2776.00000

Energy minimization has stopped, but the forces have not converged to the
requested precision Fmax < 10 (which may not be possible for your system). It
stopped because the algorithm tried to make a new step whose size was too
small, or there was no change in the energy since last step. Either way, we
regard the minimization as converged to within the available machine
precision, given your starting configuration and EM parameters.

Double precision normally gives you higher accuracy, but this is often not
needed for preparing to run molecular dynamics.
You might need to increase your constraint accuracy, or turn
off constraints altogether (set constraints = none in mdp file)

Steepest Descents converged to machine precision in 2777 steps,
but did not reach the requested Fmax < 10.
Potential Energy = -6.7911422e+04
Maximum force = 7.9560303e+03 on atom 1793
Norm of force = 1.8530257e+02

5

Unfortunately I closed the terminal window but here is the last part of my .log file for NVT calculation:
Input Parameters:
integrator = md
tinit = 0
dt = 0.002
nsteps = 100000
init-step = 0
simulation-part = 1
comm-mode = Linear
nstcomm = 100
bd-fric = 0
ld-seed = 267344549
emtol = 10
emstep = 0.01
niter = 20
fcstep = 0
nstcgsteep = 1000
nbfgscorr = 10
rtpi = 0.05
nstxout = 500
nstvout = 500
nstfout = 0
nstlog = 500
nstcalcenergy = 100
nstenergy = 500
nstxout-compressed = 0
compressed-x-precision = 1000
cutoff-scheme = Verlet
nstlist = 10
pbc = xyz
periodic-molecules = false
verlet-buffer-tolerance = 0.005
rlist = 1.378
coulombtype = Cut-off
coulomb-modifier = Potential-shift
rcoulomb-switch = 0
rcoulomb = 1.3
epsilon-r = 1
epsilon-rf = inf
vdw-type = Cut-off
vdw-modifier = Potential-shift
rvdw-switch = 0
rvdw = 1.3
DispCorr = No
table-extension = 1
fourierspacing = 0.12
fourier-nx = 0
fourier-ny = 0
fourier-nz = 0
pme-order = 4
ewald-rtol = 1e-05
ewald-rtol-lj = 0.001
lj-pme-comb-rule = Geometric
ewald-geometry = 0
epsilon-surface = 0
tcoupl = V-rescale
nsttcouple = 10
nh-chain-length = 0
print-nose-hoover-chain-variables = false
pcoupl = No
pcoupltype = Isotropic
nstpcouple = -1
tau-p = 1
compressibility (3x3):
compressibility[ 0]={ 0.00000e+00, 0.00000e+00, 0.00000e+00}
compressibility[ 1]={ 0.00000e+00, 0.00000e+00, 0.00000e+00}
compressibility[ 2]={ 0.00000e+00, 0.00000e+00, 0.00000e+00}
ref-p (3x3):
ref-p[ 0]={ 0.00000e+00, 0.00000e+00, 0.00000e+00}
ref-p[ 1]={ 0.00000e+00, 0.00000e+00, 0.00000e+00}
ref-p[ 2]={ 0.00000e+00, 0.00000e+00, 0.00000e+00}
refcoord-scaling = No
posres-com (3):
posres-com[0]= 0.00000e+00
posres-com[1]= 0.00000e+00
posres-com[2]= 0.00000e+00
posres-comB (3):
posres-comB[0]= 0.00000e+00
posres-comB[1]= 0.00000e+00
posres-comB[2]= 0.00000e+00
QMMM = false
QMconstraints = 0
QMMMscheme = 0
MMChargeScaleFactor = 1
qm-opts:
ngQM = 0
constraint-algorithm = Lincs
continuation = false
Shake-SOR = false
shake-tol = 0.0001
lincs-order = 4
lincs-iter = 1
lincs-warnangle = 30
nwall = 0
wall-type = 9-3
wall-r-linpot = -1
wall-atomtype[0] = -1
wall-atomtype[1] = -1
wall-density[0] = 0
wall-density[1] = 0
wall-ewald-zfac = 3
pull = false
awh = false
rotation = false
interactiveMD = false
disre = No
disre-weighting = Conservative
disre-mixed = false
dr-fc = 1000
dr-tau = 0
nstdisreout = 100
orire-fc = 0
orire-tau = 0
nstorireout = 100
free-energy = no
cos-acceleration = 0
deform (3x3):
deform[ 0]={ 0.00000e+00, 0.00000e+00, 0.00000e+00}
deform[ 1]={ 0.00000e+00, 0.00000e+00, 0.00000e+00}
deform[ 2]={ 0.00000e+00, 0.00000e+00, 0.00000e+00}
simulated-tempering = false
swapcoords = no
userint1 = 0
userint2 = 0
userint3 = 0
userint4 = 0
userreal1 = 0
userreal2 = 0
userreal3 = 0
userreal4 = 0
applied-forces:
electric-field:
x:
E0 = 0
omega = 0
t0 = 0
sigma = 0
y:
E0 = 0
omega = 0
t0 = 0
sigma = 0
z:
E0 = 0
omega = 0
t0 = 0
sigma = 0
density-guided-simulation:
active = false
group = protein
similarity-measure = inner-product
atom-spreading-weight = unity
force-constant = 1e+09
gaussian-transform-spreading-width = 0.2
gaussian-transform-spreading-range-in-multiples-of-width = 4
reference-density-filename = reference.mrc
nst = 1
normalize-densities = true
adaptive-force-scaling = false
adaptive-force-scaling-time-constant = 4
grpopts:
nrdf: 17997
ref-t: 300
tau-t: 0.1
annealing: No
annealing-npoints: 0
acc: 0 0 0
nfreeze: N N N
energygrp-flags[ 0]: 0

Changing nstlist from 10 to 25, rlist from 1.378 to 1.548

Using 1 MPI thread

Non-default thread affinity set, disabling internal thread affinity

Using 8 OpenMP threads

System total charge: 0.000
Potential shift: LJ r^-12: -4.292e-02 r^-6: -2.072e-01, Coulomb -8e-01
Generated table with 1273 data points for 1-4 COUL.
Tabscale = 500 points/nm
Generated table with 1273 data points for 1-4 LJ6.
Tabscale = 500 points/nm
Generated table with 1273 data points for 1-4 LJ12.
Tabscale = 500 points/nm

Using SIMD 4x8 nonbonded short-range kernels

Using a dual 4x8 pair-list setup updated with dynamic pruning:
outer list: updated every 25 steps, buffer 0.248 nm, rlist 1.548 nm
inner list: updated every 4 steps, buffer 0.021 nm, rlist 1.321 nm
At tolerance 0.005 kJ/mol/ps per atom, equivalent classical 1x1 list would be:
outer list: updated every 25 steps, buffer 0.322 nm, rlist 1.622 nm
inner list: updated every 4 steps, buffer 0.030 nm, rlist 1.330 nm

Using Lorentz-Berthelot Lennard-Jones combination rule

Removing pbc first time

++++ PLEASE READ AND CITE THE FOLLOWING REFERENCE ++++
G. Bussi, D. Donadio and M. Parrinello
Canonical sampling through velocity rescaling
J. Chem. Phys. 126 (2007) pp. 014101
-------- -------- — Thank You — -------- --------

There are: 6000 Atoms
Center of mass motion removal mode is Linear
We have the following groups for center of mass motion removal:
0: rest
Initial temperature: 297.989 K

Started mdrun on rank 0 Wed Sep 13 23:49:04 2023

       Step           Time
          0        0.00000

Energies (kJ/mol)
Bond Angle Proper Dih. LJ-14 Coulomb-14
1.65121e+02 1.78406e+03 1.01411e+03 -8.17716e+02 -5.48976e+04
LJ (SR) Coulomb (SR) Potential Kinetic En. Total Energy
-9.65450e+03 -3.29584e+03 -6.57024e+04 4.52880e+04 -2.04144e+04
Conserved En. Temperature Pressure (bar)
-2.04144e+04 6.05311e+02 7.87726e+03

6

I can upload my .gro and .top files for solvent box if it can help

7

Could this be a case of PBC wrapping hangs after a system explodes (#4766) · Issues · GROMACS / GROMACS · GitLab?

8

Maybe. But in that case it wouldn’t hang at step 0, but at step nstlist. To check that run:
mdrun -v -stepout 1
and see at which step it hangs.

The cause of this is likely too high forces. The largest force after EM is far too large. Likely there is some problem in the initial structure. Do a visual inspection of atom 1793 and its surroundings.

9

By using -stepout 1 option, calculation hangs at step 24
gmx mdrun -v -deffnm Acetonitrile_NVT -ntomp 8 -ntmpi 1 -stepout 1

Back Off! I just backed up Acetonitrile_NVT.log to ./#Acetonitrile_NVT.log.1#
Reading file Acetonitrile_NVT.tpr, VERSION 2020.7 (single precision)
Changing nstlist from 10 to 25, rlist from 1.378 to 1.548

Using 1 MPI thread

Non-default thread affinity set, disabling internal thread affinity

Using 8 OpenMP threads

Back Off! I just backed up Acetonitrile_NVT.trr to ./#Acetonitrile_NVT.trr.1#

Back Off! I just backed up Acetonitrile_NVT.edr to ./#Acetonitrile_NVT.edr.1#
starting mdrun ‘ICE’
100000 steps, 200.0 ps.
step 24

10

this is the .mdp file that I used for EM
integrator = steep
nsteps = ?
rlist = 1.0
coulombtype = cut-off
cutoff-scheme = verlet
rcoulomb = 1.0
vdw-type = cut-off
rvdw = 1.0

pbc = xyz
periodic_molecules = yes
constraints = all-bonds

** I don’t know why but my EM always stops at a number too smaller that nsteps = ? (I even added emtol = 10 and nsteps = 500000 but still EM stopped at 2777)

11

I assume mdrun sets nstlist to 25, so indeed mdrun is likely is an extremely long loop to fix PBC because an atom experienced an extremely high force.

There seems to be an issue in your initial coordinates which prevents energy minimization from reaching a reasonable state. I can’t say what that issue is. As said, have a look at the environment of atoms 1793 which reports the largest force.

12

I actually checked it and found nothing strange.
I must add that every time I’m running EM the atom experiencing highest force changes

Would you please explian briefly what should I do(I mean the process of making a solvent box for acetonitrile in charmm36 force field)? Maybe I’m missing something here

13

I don’t know how you have prepared your system and the issue can be very case specific.

A general thing you can try is running energy minimization with a double precision binary. That will likely get you to a smaller force, but this might not solve your issue.