Running an NPT simulation after NVT (domain decomposition error)

GROMACS version: 2018.4
GROMACS modification: No

Hi all,

I performed a lengthy NVT simulation (after a short NPT equilibration) of a polymer to relax it and now would like to anneal in NPT. I’m encountering the following error when I try to move to NPT:

Fatal error:
There is no domain decomposition for 24 ranks that is compatible with the
given box and a minimum cell size of 2.3995 nm
Change the number of ranks or mdrun option -rcon or -dds or your LINCS
settings
Look in the log file for details on the domain decomposition

which I know isn’t an issue with the core number because the system is large enough and has no voids or areas with vacuum when I view on VMD. I can also extend the simulation in NVT with no issues. I think there’s something going wrong with the box dimension re-scaling and now something bad is happening to the bonds over the pbc boundaries. I’m not sure how to fix. When I try to continue the simulation in NPT I use

gmx grompp -f NPT_equi.mdp -p topol.top -c confout.gro -t state.cpt -o new.tpr

gmx mdrun -s em.tpr (which I don’t pass the state.cpt file to because GROMACS advises against). I thought this would handle the pbc correctly?

My NVT .mdp file:

;
; File ‘mdout.mdp’ was generated
; By user: spoel (291)
; On host: chagall
; At date: Mon Dec 15 13:52:23 2003
;
; VARIOUS PREPROCESSING OPTIONS
title = 6PPD
; RUN CONTROL PARAMETERS
integrator = md
; Start time and timestep in ps
tinit = 0
dt = 0.002
nsteps = 200000000
; For exact run continuation or redoing part of a run
init_step = 0
; mode for center of mass motion removal
;comm-mode = None
comm-mode = Linear
;comm-mode = Angular
; number of steps for center of mass motion removal
nstcomm = 1
; group(s) for center of mass motion removal
comm-grps =

; LANGEVIN DYNAMICS OPTIONS
; Temperature, friction coefficient (amu/ps) and random seed
bd-fric = 0.5
ld-seed = 1993
; ENERGY MINIMIZATION OPTIONS
; Force tolerance and initial step-size
emtol = 100
emstep = 0.01
; Max number of iterations in relax_shells
niter = 20
; Step size (1/ps^2) for minimization of flexible constraints
fcstep = 0
; Frequency of steepest descents steps when doing CG
nstcgsteep = 100
nbfgscorr = 10

; OUTPUT CONTROL OPTIONS
; Output frequency for coords (x), velocities (v) and forces (f)
;nstxout = 2000
;nstvout = 1000
;nstfout = 1000
; Checkpointing helps you continue after crashes
nstcheckpoint = 10000
; Output frequency for energies to log file and energy file
nstlog = 2000
nstenergy = 2000
; Output frequency and precision for xtc file
nstxtcout = 2000 ; Schreibt alle nstxtcout schritte aus. 1000 gewhält da Korrelationszeit gleich 1 psec.
compressed-x-precision = 1000 ; nicht ändern !!!
; This selects the subset of atoms for the xtc file. You can
; select multiple groups. By default all atoms will be written.
xtc-grps =

; NEIGHBORSEARCHING PARAMETERS
; nblist update frequency
cutoff-scheme = Verlet
nstlist = 10
; ns algorithm (simple or grid)
ns_type = grid
; Periodic boundary conditions: xyz (default), no (vacuum)
; or full (infinite systems only)
;pbc = no
;pbc = full
pbc = xyz
;~ pbc = xy
; nblist cut-off
rlist = 1.1
domain-decomposition = yes

; OPTIONS FOR ELECTROSTATICS AND VDW
; Method for doing electrostatics
coulombtype = PME ;PPPM ;Cut-off PME;PME;
;~ rcoulomb-switch = 1.0
rcoulomb = 1.2
; Dielectric constant (DC) for cut-off or DC of reaction field
;epsilon-r = 1
; Method for doing Van der Waals
vdw-type = shift
; cut-off lengths
rvdw-switch = 1.0
rvdw = 1.2
; Apply long range dispersion corrections for Energy and Pressure
DispCorr = EnerPres
; Extension of the potential lookup tables beyond the cut-off
table-extension = 1;Erweiterung der nicht gebundenen WW ueber den cut_off hinaus
;-> hilft dabei energygrps nicht zu stören
; Spacing for the PME/PPPM FFT grid
fourierspacing = 0.12
; FFT grid size, when a value IS 0 fourierspacing will be used
fourier_nx = 0
fourier_ny = 0
fourier_nz = 0
; EWALD/PME/PPPM parameters
pme_order = 4
ewald_rtol = 1e-05
ewald_geometry = 3d
;~ ewald_geometry = 3dc;anpassung der ewald summation auf 2d pbc
epsilon_surface = 0
optimize_fft = no

; GENERALIZED BORN ELECTROSTATICS
; Algorithm for calculating Born radii
gb_algorithm = Still
; Frequency of calculating the Born radii inside rlist
nstgbradii = 1
; Cutoff for Born radii calculation; the contribution from atoms
; between rlist and rgbradii is updated every nstlist steps
rgbradii = 2
; Salt concentration in M for Generalized Born models
gb_saltconc = 0

; IMPLICIT SOLVENT (for use with Generalized Born electrostatics)
implicit_solvent = No

; OPTIONS FOR WEAK COUPLING ALGORITHMS
; Temperature coupling
;Tcoupl = no
Tcoupl = v-rescale
; Groups to couple separately
tc-grps = System
; Time constant (ps) and reference temperature (K)
tau_t = 0.100
ref_t = 298.1
; Pressure coupling
Pcoupl = no
;Pcoupl = no
Pcoupltype = Isotropic
; Time constant (ps), compressibility (1/bar) and reference P (bar) =
tau-p = 0.5
compressibility = 4.5e-5
ref-p = 1.013
refcoord-scaling = all

; Random seed for Andersen thermostat
andersen_seed = 815131

; SIMULATED ANNEALING
; Type of annealing for each temperature group (no/single/periodic)
annealing = single
; Number of time points to use for specifying annealing in each group
annealing_npoints = 4
; List of times at the annealing points for each group
annealing_time = 0 50000 200000 250000
; Temp. at each annealing point, for each group.
annealing_temp = 298.1 500 500 298.1

; GENERATE VELOCITIES FOR STARTUP RUN
gen_vel = yes
gen_temp = 298.1
gen_seed = -1 ;1993

; OPTIONS FOR BONDS
;constraints = none ; Für Torsionsuntersuchung ausschalten. Warum auch immer …
constraints = all-bonds
; Type of constraint algorithm
constraint-algorithm = Lincs
; Do not constrain the start configuration
unconstrained-start = no
; Use successive overrelaxation to reduce the number of shake iterations
Shake-SOR = no
; Relative tolerance of shake
shake-tol = 1e-04
; Highest order in the expansion of the constraint coupling matrix
lincs-order = 4 ; before was 4
; Number of iterations in the final step of LINCS. 1 is fine for
; normal simulations, but use 2 to conserve energy in NVE runs.
; For energy minimization with constraints it should be 4 to 8.
lincs-iter = 1
; Lincs will write a warning to the stderr if in one step a bond
; rotates over more degrees than
lincs-warnangle = 90 ; # Von 30 auf 360 gesetzt, da Simulation sonst abbricht !!!
; Convert harmonic bonds to morse potentials
morse = no

; ENERGY GROUP EXCLUSIONS
; Pairs of energy groups for which all non-bonded interactions are excluded
; Selection of energy groups

;Wechselwirkung zwischen den Oberflaechen rausnehmen
;~ energygrps = S1r S2r O1r O2r S1l S2l O1l O2l
;~ energygrp_excl = S1r S1r S2r S1r S2r S2r O1r S1r O1r S2r O1r O1r O2r S1r O2r S2r O2r O1r O2r O2r S1l S1l S2l S1l S2l S2l O1l S1l O1l S2l O1l O1l O2l S1l O2l S2l O2l O1l O2l O2l
; WALLS
; Number of walls, type, atom types, densities and box-z scale factor for Ewald
nwall = 0
wall_type = 9-3;LJ integriert über das Volumen hinter der Wand mit 9-3-Potential
;wall_type = 10-4;LJ integriert über die Oberflaeche der Wand mit 10-4-Potential
wall_atomtype = Siw Siw
wall_density = 1 1;Dichte der Oberflaechen in nm^2(10-4) oder nm^3(9-3)
wall_ewald_zfac = 3
wall_r_linpot = 1.0;ab hier Poential linear fortgsesetzt d.h. an r_cut anpassen und
;insgesamt aufpassen was die Waende machen

and my NPT:

; File ‘mdout.mdp’ was generated
; By user: spoel (291)
; On host: chagall
; At date: Mon Dec 15 13:52:23 2003
;
; VARIOUS PREPROCESSING OPTIONS
title = 6PPD
; RUN CONTROL PARAMETERS
integrator = md
; Start time and timestep in ps
tinit = 0
dt = 0.002
nsteps = 2000000
; For exact run continuation or redoing part of a run
init_step = 0
; mode for center of mass motion removal
;comm-mode = None
comm-mode = Linear
;comm-mode = Angular
; number of steps for center of mass motion removal
nstcomm = 1
; group(s) for center of mass motion removal
comm-grps =

; LANGEVIN DYNAMICS OPTIONS
; Temperature, friction coefficient (amu/ps) and random seed
bd-fric = 0.5
ld-seed = 1993
; ENERGY MINIMIZATION OPTIONS
; Force tolerance and initial step-size
emtol = 100
emstep = 0.01
; Max number of iterations in relax_shells
niter = 20
; Step size (1/ps^2) for minimization of flexible constraints
fcstep = 0
; Frequency of steepest descents steps when doing CG
nstcgsteep = 100
nbfgscorr = 10

; OUTPUT CONTROL OPTIONS
; Output frequency for coords (x), velocities (v) and forces (f)
;nstxout = 2000
;nstvout = 1000
;nstfout = 1000
; Checkpointing helps you continue after crashes
nstcheckpoint = 10000
; Output frequency for energies to log file and energy file
nstlog = 2000
nstenergy = 2000
; Output frequency and precision for xtc file
nstxtcout = 2000 ; Schreibt alle nstxtcout schritte aus. 1000 gewhält da Korrelationszeit gleich 1 psec.
compressed-x-precision = 1000 ; nicht ändern !!!
; This selects the subset of atoms for the xtc file. You can
; select multiple groups. By default all atoms will be written.
xtc-grps =

; NEIGHBORSEARCHING PARAMETERS
; nblist update frequency
cutoff-scheme = Verlet
nstlist = 10
; ns algorithm (simple or grid)
ns_type = grid
; Periodic boundary conditions: xyz (default), no (vacuum)
; or full (infinite systems only)
;pbc = no
;pbc = full
pbc = xyz
;~ pbc = xy
; nblist cut-off
rlist = 1.1
domain-decomposition = yes

; OPTIONS FOR ELECTROSTATICS AND VDW
; Method for doing electrostatics
coulombtype = Cut-off;PPPM ;Cut-off PME;PME;
;~ rcoulomb-switch = 1.0
rcoulomb = 1.2
; Dielectric constant (DC) for cut-off or DC of reaction field
;epsilon-r = 1
; Method for doing Van der Waals
vdw-type = shift
; cut-off lengths
rvdw-switch = 1.0
rvdw = 1.2
; Apply long range dispersion corrections for Energy and Pressure
DispCorr = EnerPres
; Extension of the potential lookup tables beyond the cut-off
table-extension = 1;Erweiterung der nicht gebundenen WW ueber den cut_off hinaus
;-> hilft dabei energygrps nicht zu stören
; Spacing for the PME/PPPM FFT grid
fourierspacing = 0.12
; FFT grid size, when a value IS 0 fourierspacing will be used
fourier_nx = 0
fourier_ny = 0
fourier_nz = 0
; EWALD/PME/PPPM parameters
pme_order = 4
ewald_rtol = 1e-05
ewald_geometry = 3d
;~ ewald_geometry = 3dc;anpassung der ewald summation auf 2d pbc
epsilon_surface = 0
optimize_fft = no

; GENERALIZED BORN ELECTROSTATICS
; Algorithm for calculating Born radii
gb_algorithm = Still
; Frequency of calculating the Born radii inside rlist
nstgbradii = 1
; Cutoff for Born radii calculation; the contribution from atoms
; between rlist and rgbradii is updated every nstlist steps
rgbradii = 2
; Salt concentration in M for Generalized Born models
gb_saltconc = 0

; IMPLICIT SOLVENT (for use with Generalized Born electrostatics)
implicit_solvent = No

; OPTIONS FOR WEAK COUPLING ALGORITHMS
; Temperature coupling
;Tcoupl = no
Tcoupl = v-rescale
; Groups to couple separately
tc-grps = System
; Time constant (ps) and reference temperature (K)
tau_t = 0.100
ref_t = 298
; Pressure coupling
Pcoupl = Parrinello-Rahman
;Pcoupl = no
Pcoupltype = Isotropic
; Time constant (ps), compressibility (1/bar) and reference P (bar) =
tau-p = 4 ; was 0.5
compressibility = 4.5e-5
ref-p = 1.013
refcoord-scaling = all

; Random seed for Andersen thermostat
andersen_seed = 815131

; SIMULATED ANNEALING
; Type of annealing for each temperature group (no/single/periodic)
annealing = single
; Number of time points to use for specifying annealing in each group
annealing_npoints = 4
; List of times at the annealing points for each group
annealing_time = 0 20000 70000 125000
; Temp. at each annealing point, for each group.
annealing_temp = 298.1 450 450 150

; GENERATE VELOCITIES FOR STARTUP RUN
gen_vel = no
gen_temp = 298.1
gen_seed = -1 ;1993

; OPTIONS FOR BONDS
;constraints = none ; Für Torsionsuntersuchung ausschalten. Warum auch immer …
constraints = all-bonds
; Type of constraint algorithm
constraint-algorithm = Lincs
; Do not constrain the start configuration
unconstrained-start = no
; Use successive overrelaxation to reduce the number of shake iterations
Shake-SOR = no
; Relative tolerance of shake
shake-tol = 1e-04
; Highest order in the expansion of the constraint coupling matrix
lincs-order = 12 ; before was 4
; Number of iterations in the final step of LINCS. 1 is fine for
; normal simulations, but use 2 to conserve energy in NVE runs.
; For energy minimization with constraints it should be 4 to 8.
lincs-iter = 1
; Lincs will write a warning to the stderr if in one step a bond
; rotates over more degrees than
lincs-warnangle = 90 ; # Von 30 auf 360 gesetzt, da Simulation sonst abbricht !!!
; Convert harmonic bonds to morse potentials
morse = no

any help would be much appreciated

I would think that the error message is correct. Why do you say “which I know isn’t an issue with the core number because the system is large enough”? Are you really sure it’s not an issue? 24 ranks is quite a lot. How large is your system (in atoms and in box dimensions)?

I would recommend trying 1 to 6 ranks. You can still use as many CPU cores.

Thank you for the response

I have 32000 atoms with a box size of approximately 6.5nm in all directions. I really don’t think this is the problem, I have used this system many times with the same setup and there are no problems. I have also tried reducing to 8 cores but I get the same error. Even this system I have the trouble with has already had an NPT equilibration followed by NVT and it was okay. The simulation isn’t even starting in this case, just immediately crashing.

With only 32000 atoms I would recommend one rank (still using all 24 cores), but in theory I think 8 ranks should work.

If you think there is some kind of bug involved, I would recommend switching to a supported version of GROMACS, preferably 2024.3, and see if it helps.