Energy minimization stops abruptly without force convergence and segmentation fault

GROMACS version: 2021.3
GROMACS modification: Yes/No
Hello everyone. I am trying to simulate a protein-DNA complex in water system using ff19SB for protein, ZAFF for zinc-coordinated center and OL24 for DNA. I prepared the structure using AmberTools25 and converted the files to GROMACS format using Parmed. As per the AmberTools manual, I used the terminal_monophosphate.lib to handle the 5’ phosphate group on DNA termini. I then created index groups for each molecule and their respective heavy atoms and generated position restraint files for the massive atoms to be used during equilibration.

While running energy minimization on the system, however, I get hit with this message:

Energy minimization has stopped, but the forces have not converged to the
requested precision Fmax < 1000 (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 73 steps,
but did not reach the requested Fmax < 1000.
Potential Energy  = -4.6988097e+05
Maximum force     =  7.0742570e+04 on atom 1447
Norm of force     =  5.4708254e+02

I referenced certain other posts on this forum on the same message and the GROMACS manual too under the section “Errors in mdrun”, but could not get rid of the message. I shall share my .mdp file for reference:

; minim.mdp - used as input into grompp to generate em.tpr
; Parameters describing what to do, when to stop and what to save

integrator  = steep         ; Algorithm (steep = steepest descent minimization)
emtol       = 1000.0        ; Stop minimization when the maximum force < 1000.0 kJ/mol/nm
emstep      = 0.001          ; Minimization step size
nsteps      = 50000         ; Maximum number of (minimization) steps to perform

; Parameters describing how to find the neighbors of each atom and how to calculate the interactions
nstlist         = 1         ; Frequency to update the neighbor list and long range forces
cutoff-scheme   = Verlet    ; Buffered neighbor searching
ns_type         = grid      ; Method to determine neighbor list (simple, grid)
coulombtype     = PME       ; Treatment of long range electrostatic interactions
rcoulomb        = 1.0       ; Short-range electrostatic cut-off
rvdw            = 1.0       ; Short-range Van der Waals cut-off
pbc             = xyz       ; Periodic Boundary Conditions in all 3 dimensions

Someone suggested that I re-run the output of the current minimization for another round of minimization under a different minimization algorithm. Hence I chose conjugate gradient. That did not work either and I got this message instead:

[gpu006:37691] *** Process received signal ***
[gpu006:37691] Signal: Segmentation fault (11)
[gpu006:37691] Signal code: Address not mapped (1)
[gpu006:37691] Failing at address: 0xfffffffe03f26988
[gpu006:37691] [ 0] /lib64/libpthread.so.0(+0xf630)[0x7f7b940e5630]
[gpu006:37691] [ 1] /home/apps/spack/opt/spack/linux-centos7-cascadelake/gcc-11.2.0/gromacs-2021.3-76vdwmvrckesru4sd5yejtb5qe74phks/bin/../lib64/libgromacs_mpi.so.6(+0x9e125a)[0x7f7b9542425a]
[gpu006:37691] [ 2] /home/apps/spack/opt/spack/linux-centos7-cascadelake/gcc-11.2.0/gromacs-2021.3-76vdwmvrckesru4sd5yejtb5qe74phks/bin/../lib64/libgromacs_mpi.so.6(_Z8do_pairsiiPKiPK9t_iparamsPA3_KfPA4_fPA3_fPK5t_pbcPS4_PfPK9t_mdatomsPK10t_forcerecbRKN3gmx12StepWorkloadEP17gmx_grppairener_tPi+0xbe2)[0x7f7b954256a2]
[gpu006:37691] [ 3] /home/apps/spack/opt/spack/linux-centos7-cascadelake/gcc-11.2.0/gromacs-2021.3-76vdwmvrckesru4sd5yejtb5qe74phks/bin/../lib64/libgromacs_mpi.so.6(+0x9d75b3)[0x7f7b9541a5b3]
[gpu006:37691] [ 4] /home/apps/spack/opt/spack/linux-centos7-cascadelake/gcc-11.2.0/gromacs-2021.3-76vdwmvrckesru4sd5yejtb5qe74phks/bin/../lib64/libgromacs_mpi.so.6(+0x9d7bbc)[0x7f7b9541abbc]
[gpu006:37691] [ 5] /home/apps/spack/opt/spack/linux-centos7-skylake_avx512/gcc-8.3.0/gcc-11.2.0-fufsigm7go2apn3xzygaw3x266zuxav4/lib64/libgomp.so.1(+0x1b4b6)[0x7f7b9430d4b6]
[gpu006:37691] [ 6] /lib64/libpthread.so.0(+0x7ea5)[0x7f7b940ddea5]
[gpu006:37691] [ 7] /lib64/libc.so.6(clone+0x6d)[0x7f7b934e4b0d]
[gpu006:37691] *** End of error message ***
/var/share/slurm/d/job02335/slurm_script: line 19: 37691 Segmentation fault gmx_mpi mdrun -v -deffnm em_1 -ntomp 40

I still, however, continued with my first NVT equilibration under position restraints of 1000 kJ/mol/nm, but ran into the same segmentation fault error.

GROMACS:      gmx mdrun, version 2021.3
Executable:   /home/apps/spack/opt/spack/linux-centos7-cascadelake/gcc-11.2.0/gromacs-2021.3-76vdwmvrckesru4sd5yejtb5qe74phks/bin/gmx_mpi
Data prefix:  /home/apps/spack/opt/spack/linux-centos7-cascadelake/gcc-11.2.0/gromacs-2021.3-76vdwmvrckesru4sd5yejtb5qe74phks
Working dir:  /scratch/priyanshu.das/complex/nvt
Command line:
gmx_mpi mdrun -v -deffnm nvt_300 -ntomp 12 -nb gpu -pme gpu -update gpu -pin on

Back Off! I just backed up nvt_300.log to ./#nvt_300.log.1#
Reading file nvt_300.tpr, VERSION 2021.3 (single precision)
Changing nstlist from 10 to 100, rlist from 1 to 1.144

1 GPU selected for this run.
Mapping of GPU IDs to the 2 GPU tasks in the 1 rank on this node:
PP:0,PME:0
PP tasks will do (non-perturbed) short-ranged interactions on the GPU
PP task will update and constrain coordinates on the GPU
PME tasks will do all aspects on the GPU
Using 1 MPI process
Using 12 OpenMP threads

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

Back Off! I just backed up nvt_300.edr to ./#nvt_300.edr.1#
starting mdrun ‘THAP domain in complex with DNA’
500000 steps,   1000.0 ps.
^Mstep 0^Mstep 100, will finish Mon Nov 10 21:34:04 2025^Mstep 200, remaining wall clock time:   267 s          ^Mstep 300, remaining wall clock time:   238 s          ^Mstep 400, remaining wall clock time:   221 s          ^Mstep 500, remaining wall clock time:   210 s          ^Mstep 600, remaining wall clock time:   203 s          ^Mstep 700, remaining wall clock time:   198 s          ^Mstep 800, remaining wall clock time:   194 s          ^Mstep 900, remaining wall clock time:   191 s          /var/share/slurm/d/job02342/slurm_script: line 19:   455 Segmentation fault      gmx_mpi mdrun -v -deffnm nvt_300 -ntomp 12 -nb gpu -pme gpu -update gpu -pin on

Can someone please guide as to what is happening and how I can correct this? If it helps, I am attaching my .top file with this post.

topol.top (1.1 MB)

My guess is that you have some seriously overlapping atoms, maybe caused by a misplaced sidechain. A maximum force of 7e4 after 73 steps is very high. Have you looked at the structure around atom 1447?

Thank you so much @hess! Upon inspection in PyMOL, I found that atom 1450 was in a clash with neighbouring atoms due to poor geometry. After manual correction, EM using steepest descent could converge to under emtol = 1000

I’m getting a different kind of error now while trying to perform an NVT equilibration after energy minimization.

PP tasks will do (non-perturbed) short-ranged interactions on the GPU
PP task will update and constrain coordinates on the GPU
PME tasks will do all aspects on the GPU
Using 1 MPI process
Using 40 OpenMP threads

starting mdrun 'THAP domain in complex with DNA in water'
500000 steps,   1000.0 ps.
^Mstep 0^Mstep 100, will finish Tue Nov 11 14:05:43 2025^Mstep 200, will finish Tue Nov 11 14:03:47 2025^Mstep 300, will finish Tue Nov 11 14:03:10 2025^Mstep 400, will finish Tue Nov 11 14:02:47 2025^Mstep 500, remaining wall clock time:   287 s          ^Mstep 600, remaining wall clock time:   277 s          ^Mstep 700, remaining wall clock time:   271 s          ^Mstep 800, remaining wall clock time:   263 s          ^Mstep 900, remaining wall clock time:   257 s          ^Mstep 1000, remaining wall clock time:   252 s          ^Mstep 1100, remaining wall clock time:   249 s          ^Mstep 1200, remaining wall clock time:   246 s
^Mstep 1300, remaining wall clock time:   243 s          ^Mstep 1400, remaining wall clock time:   241 s          ^Mstep 1500, remaining wall clock time:   239 s          ^Mstep 1600, remaining wall clock time:   239 s          ^Mstep 1700, remaining wall clock time:   238 s          ^Mstep 1800, remaining wall clock time:   236 s          ^Mstep 1900, remaining wall clock time:   235 s          ^Mstep 2000, remaining wall clock time:   234 s          ^Mstep 2100, remaining wall clock time:   232 s          ^Mstep 2200, remaining wall clock time:   231 s
-------------------------------------------------------
Program:     gmx mdrun, version 2021.3
Source file: src/gromacs/gpu_utils/devicebuffer.cuh (line 136)
Function:    copyToDeviceBuffer<gmx::BasicVector<float> >(gmx::BasicVector<float>**, const gmx::BasicVector<float>*, size_t, size_t, const DeviceStream&, GpuApiCallBehavior, CommandEvent*)::<lambda()>

Assertion failed:
Condition: stat == cudaSuccess
Asynchronous H2D copy failed. CUDA error #700 (cudaErrorIllegalAddress): an
illegal memory access was encountered.

For more information and tips for troubleshooting, please check the GROMACS
website at http://www.gromacs.org/Documentation/Errors
-------------------------------------------------------
--------------------------------------------------------------------------
MPI_ABORT was invoked on rank 0 in communicator MPI_COMM_WORLD
  Proc: [[34077,0],0]
  Errorcode: 1

NOTE: invoking MPI_ABORT causes Open MPI to kill all MPI processes.
You may or may not see output from other processes, depending on
exactly when Open MPI kills them.

Any idea what this means and how to resolve this? Thank you!

If you have that high forces your system will explode within a few steps of MD. I think that is what you are seeing.

I don’t know if it makes sense, but after correcting for clashes, I ran a steepest descent minimization to 500 kJ/mol/nm^2 followed by another minimization using conjugate gradient to bring the forces under 100. It still is either of these issues: LINCS warning or errors like the one I posted. In all these cases, it’s the 5’ phosphate group atoms of the DNA that’s causing them. Should I rather simulate using 5’ OH instead of the phosphate?

I think it should in principle work with any group the force field provides. My suspicion is that your initial configuration is problematic.

I corrected some geometry further and could achieve energy minimization to under 50 using the earlier method. I realized that I only get LINCS warning when I’m using a time step > 0.5fs. Otherwise, the simulation runs stably. Could this have been the issue, along with the geometry as you had mentioned?

Are you using constraints=h-bonds? If so, you should be able to use a timestep of 2 fs. In extreme cases you might need to equilibrate for a few 100 ps with a smaller timestep.