Error in minimization step - adding Br ion to amber force field ions.itp

User discussions
1

GROMACS version:
GROMACS modification: Yes/No
Hi everyone! I am facing error in minimization step during simulation of a protein. I have used amber99sb-ILDN forcefield. for this forcefield, Br ion parameters were available in ffnonbonded.itp and atomtypes.itp but it was not there in ions.itp. so, i added Br in ions.itp and added Li and Br with conc. 0.15 as counter ions to neautalise the system. but in the minimization step I am facing error.

Command line:
gmx_mpi mdrun -v -deffnm em

Reading file em.tpr, VERSION 2021.4 (single precision)
Using 1 MPI process
Using 32 OpenMP threads

Steepest Descents:
Tolerance (Fmax) = 1.00000e+03
Number of steps = 50000
Step= 0, Dmax= 1.0e-05 nm, Epot= -1.88811e+05 Fmax= 1.40288e+05, atom= 131
Step= 1, Dmax= 1.0e-05 nm, Epot= -1.88851e+05 Fmax= 1.40111e+05, atom= 131
Step= 2, Dmax= 1.2e-05 nm, Epot= -1.88899e+05 Fmax= 1.39901e+05, atom= 131
Step= 3, Dmax= 1.4e-05 nm, Epot= -1.88957e+05 Fmax= 1.39651e+05, atom= 131
Step= 4, Dmax= 1.7e-05 nm, Epot= -1.89026e+05 Fmax= 1.39348e+05, atom= 131
Step= 5, Dmax= 2.1e-05 nm, Epot= -1.89109e+05 Fmax= 1.38987e+05, atom= 131
Step= 6, Dmax= 2.5e-05 nm, Epot= -1.89208e+05 Fmax= 1.38555e+05, atom= 131
Step= 7, Dmax= 3.0e-05 nm, Epot= -1.89328e+05 Fmax= 1.38038e+05, atom= 131
Step= 8, Dmax= 3.6e-05 nm, Epot= -1.89471e+05 Fmax= 1.37422e+05, atom= 131
Step= 9, Dmax= 4.3e-05 nm, Epot= -1.89642e+05 Fmax= 1.36686e+05, atom= 131
Step= 10, Dmax= 5.2e-05 nm, Epot= -1.89846e+05 Fmax= 1.35805e+05, atom= 131
Step= 11, Dmax= 6.2e-05 nm, Epot= -1.90091e+05 Fmax= 1.34755e+05, atom= 131
Step= 12, Dmax= 7.4e-05 nm, Epot= -1.90383e+05 Fmax= 1.33506e+05, atom= 131
Step= 13, Dmax= 8.9e-05 nm, Epot= -1.90733e+05 Fmax= 1.32019e+05, atom= 131
Step= 14, Dmax= 1.1e-04 nm, Epot= -1.91149e+05 Fmax= 1.30254e+05, atom= 131
Step= 15, Dmax= 1.3e-04 nm, Epot= -1.91646e+05 Fmax= 1.28159e+05, atom= 131
Step= 16, Dmax= 1.5e-04 nm, Epot= -1.92237e+05 Fmax= 1.25682e+05, atom= 131
Step= 17, Dmax= 1.8e-04 nm, Epot= -1.92939e+05 Fmax= 1.22763e+05, atom= 131
Step= 18, Dmax= 2.2e-04 nm, Epot= -1.93771e+05 Fmax= 1.19338e+05, atom= 131
Step= 19, Dmax= 2.7e-04 nm, Epot= -1.94755e+05 Fmax= 1.15334e+05, atom= 131
Step= 20, Dmax= 3.2e-04 nm, Epot= -1.95916e+05 Fmax= 1.10673e+05, atom= 131
Step= 21, Dmax= 3.8e-04 nm, Epot= -1.97282e+05 Fmax= 1.05283e+05, atom= 131
Step= 22, Dmax= 4.6e-04 nm, Epot= -1.98883e+05 Fmax= 9.90926e+04, atom= 131
Step= 23, Dmax= 5.5e-04 nm, Epot= -2.00754e+05 Fmax= 9.20480e+04, atom= 131
Step= 24, Dmax= 6.6e-04 nm, Epot= -2.02933e+05 Fmax= 8.41092e+04, atom= 131
Step= 25, Dmax= 7.9e-04 nm, Epot= -2.05462e+05 Fmax= 7.64873e+04, atom= 7598
Step= 26, Dmax= 9.5e-04 nm, Epot= -2.08343e+05 Fmax= 6.83805e+04, atom= 7598
Step= 27, Dmax= 1.1e-03 nm, Epot= -2.11600e+05 Fmax= 5.98322e+04, atom= 7598
Step= 28, Dmax= 1.4e-03 nm, Epot= -2.15259e+05 Fmax= 5.10455e+04, atom= 7598
Step= 29, Dmax= 1.6e-03 nm, Epot= -2.19345e+05 Fmax= 4.22813e+04, atom= 7598
Step= 30, Dmax= 2.0e-03 nm, Epot= -2.23886e+05 Fmax= 3.40412e+04, atom= 492
Step= 31, Dmax= 2.4e-03 nm, Epot= -2.28891e+05 Fmax= 2.83132e+04, atom= 1595
Step= 32, Dmax= 2.8e-03 nm, Epot= -2.34057e+05 Fmax= 2.40559e+04, atom= 1595
Step= 33, Dmax= 3.4e-03 nm, Epot= -2.39206e+05 Fmax= 1.98149e+04, atom= 1595
Step= 34, Dmax= 4.1e-03 nm, Epot= -2.44464e+05 Fmax= 1.57188e+04, atom= 1595
Step= 35, Dmax= 4.9e-03 nm, Epot= -2.49959e+05 Fmax= 1.20305e+04, atom= 8381
Step= 36, Dmax= 5.9e-03 nm, Epot= -2.55801e+05 Fmax= 9.30679e+03, atom= 8381
Step= 37, Dmax= 7.1e-03 nm, Epot= -2.61892e+05 Fmax= 7.01439e+03, atom= 8381
Step= 38, Dmax= 8.5e-03 nm, Epot= -2.68506e+05 Fmax= 5.09389e+03, atom= 8381
Step= 39, Dmax= 1.0e-02 nm, Epot= -2.76081e+05 Fmax= 3.52033e+03, atom= 2231
Step= 40, Dmax= 1.2e-02 nm, Epot= -2.85188e+05 Fmax= 2.49421e+03, atom= 2231
Step= 41, Dmax= 1.5e-02 nm, Epot= -2.95630e+05 Fmax= 4.17521e+03, atom= 596
Step= 42, Dmax= 1.8e-02 nm, Epot= -3.00591e+05 Fmax= 2.12663e+04, atom= 596
Step= 43, Dmax= 2.1e-02 nm, Epot= -3.01982e+05 Fmax= 1.11674e+04, atom= 596
Step= 44, Dmax= 2.5e-02 nm, Epot= -3.03980e+05 Fmax= 2.71631e+04, atom= 596
Step= 45, Dmax= 3.0e-02 nm, Epot= -3.05338e+05 Fmax= 1.95931e+04, atom= 596
Step= 46, Dmax= 3.7e-02 nm, Epot= -3.06607e+05 Fmax= 3.56677e+04, atom= 596
Step= 47, Dmax= 4.4e-02 nm, Epot= -3.07972e+05 Fmax= 3.10439e+04, atom= 596
Step= 48, Dmax= 5.3e-02 nm, Epot= -3.08790e+05 Fmax= 4.78081e+04, atom= 596
Step= 49, Dmax= 6.3e-02 nm, Epot= -3.10296e+05 Fmax= 4.64267e+04, atom= 596
Step= 50, Dmax= 7.6e-02 nm, Epot= -3.19027e+05 Fmax= 1.83250e+06, atom= 23043
Step= 51, Dmax= 9.1e-02 nm, Epot= -3.10219e+05 Fmax= 5.97203e+04, atom= 596
Step= 52, Dmax= 4.6e-02 nm, Epot= -3.10312e+05 Fmax= 6.19996e+04, atom= 596
Step= 53, Dmax= 2.3e-02 nm, Epot= -3.11147e+05 Fmax= 6.39242e+04, atom= 23043
Step= 54, Dmax= 1.1e-02 nm, Epot= -3.13970e+05 Fmax= 4.25138e+05, atom= 23043
Step= 55, Dmax= 5.7e-03 nm, Epot= -3.27429e+05 Fmax= 5.94123e+06, atom= 23043
Step= 56, Dmax= 6.8e-03 nm, Epot= -3.15113e+05 Fmax= 6.95023e+05, atom= 23043
Step= 57, Dmax= 3.4e-03 nm, Epot= -3.27575e+05 Fmax= 6.06913e+06, atom= 23043
Step= 58, Dmax= 4.1e-03 nm, Epot= -3.20865e+05 Fmax= 2.48044e+06, atom= 23043
Step= 59, Dmax= 2.0e-03 nm, Epot= -3.72426e+05 Fmax= 6.97644e+07, atom= 23043
Step= 60, Dmax= 2.5e-03 nm, Epot= -3.23472e+05 Fmax= 3.69377e+06, atom= 23043
Step= 61, Dmax= 1.2e-03 nm, Epot= -3.46829e+05 Fmax= 2.49015e+07, atom= 23043
Step= 62, Dmax= 6.1e-04 nm, Epot= -4.02587e+05 Fmax= 7.60497e+07, atom= 23043
Step= 63, Dmax= 7.4e-04 nm, Epot= -3.58745e+05 Fmax= 4.29664e+07, atom= 23043
Step= 64, Dmax= 3.7e-04 nm, Epot= -4.02588e+05 Fmax= 1.04434e+08, atom= 23043
Step= 65, Dmax= 4.4e-04 nm, Epot= -4.02587e+05 Fmax= 1.12412e+08, atom= 23043
Step= 66, Dmax= 2.2e-04 nm, Epot= -4.02587e+05 Fmax= 4.14473e+06, atom= 23043
Step= 67, Dmax= 1.1e-04 nm, Epot= -4.02588e+05 Fmax= 5.02694e+07, atom= 23043
Step= 68, Dmax= 5.5e-05 nm, Epot= -4.02588e+05 Fmax= 7.73021e+07, atom= 23043
Step= 69, Dmax= 2.8e-05 nm, Epot= -4.02588e+05 Fmax= 9.08984e+07, atom= 23043
Step= 70, Dmax= 1.4e-05 nm, Epot= -4.02588e+05 Fmax= 9.77319e+07, atom= 23043
Step= 71, Dmax= 6.9e-06 nm, Epot= -4.02588e+05 Fmax= 1.01115e+08, atom= 23043
Step= 72, Dmax= 3.5e-06 nm, Epot= -4.02588e+05 Fmax= 1.02661e+08, atom= 23043
Step= 73, Dmax= 1.7e-06 nm, Epot= -4.02588e+05 Fmax= 1.03568e+08, atom= 23043
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)

writing lowest energy coordinates.

Steepest Descents converged to machine precision in 74 steps,
but did not reach the requested Fmax < 1000.
Potential Energy = -4.0258788e+05
Maximum force = 1.0443410e+08 on atom 23043
Norm of force = 9.7191676e+05

please find the em.mdp file used:

; 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.01 ; 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 = 20 ; 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

please help me to rectify this error. I am new to gromacs. I have read many discussion related to this error.but i was unsuccessful. Requesting you for step by step resolution.

2

Hi,

Quoting the manual:

This may not be an error as such. It is simply informing you that during the energy minimization
process mdrun reached the limit possible to minimize the structure with your current parameters.
It does not mean that the system has not been minimized fully, but in some situations that may be
the case.

There should a confout.gro among the output, have you tried to run a md simulation with it?

3

Thank you for your reply. No, I have not got any confout.gro file. Actually after adding ions using genion command : gmx_mpi genion -s ions.tpr -o model1-H2-sol-ions.gro -pname LI -nname BR -conc 0.30 -neutral -p topol.top >genion.out
I tried to run a batch file with following commands, gmx_mpi grompp -f em.mdp -c model1-H2-sol-ions.gro -p topol.top -o em.tpr -maxwarn 10
gmx_mpi mdrun -v -deffnm em
gmx_mpi grompp -f nvt.mdp -c em.gro -r em.gro -p topol.top -o nvt.tpr -maxwarn 10
gmx_mpi mdrun -v -deffnm nvt
gmx_mpi grompp -f npt.mdp -c nvt.gro -r nvt.gro -t nvt.cpt -p topol.top -o npt.tpr -maxwarn 10
gmx_mpi mdrun -v -deffnm npt
gmx_mpi grompp -f md.mdp -c npt.gro -t npt.cpt -p topol.top -o md_0_1.tpr -maxwarn 10

But in the minimisation step I am facing this error and its not going further for npt and nvt.

Please find the attached error file.
error_3565.dat (27.7 KB)

4

What warnings are you overriding by doing this? As a general rule, never use -maxwarn or get in the habit of just incorporating it into workflows and scripts. You could be suppressing something crucially important. The maximum force in your system is absolutely unreasonable, so you need to pay attention to any messages that grompp might produce.

5

The LJ parameters for the Br atom type are not for bromide, they are for elemental bromine. You cannot simply create an ion from it and expect the resulting parameters to be valid.

6

Hi jalemkul. Thank you for your reply. Actually i have downloaded amber99sb-ILDN tar file from gromacs user contribution webpage. Following were the ions.itp and atomtypes.atp file obtained from original amber99sb-ildn.ff.

I found only chlorine ion was present in ions.itp and its name was also chlorine in atomtypes.atp. So, I assumed Bromine , fluorine and Iodine which are listed in atomtypes files can be added to ions.itp. I added then used them for simulations. They worked well, only I am facing error with Bromine.
amber99sb-ildn.ff -atomtypes.top (3.7 KB)
ions.dat (2.1 KB)

7

There is an official implementation of ff99SB-ILDN, so you should check that it matches what is officially supported. It may work out that transplanting the LJ to an ion works, but you can’t necessarily guarantee it. Other than that, you need to be looking at the output structure to see what is going on around the atom experiencing the maximum force. It could be something entirely unrelated to Br-.

8

I visualised the system prepared after adding ions and solvent, which was used as input for energy minimisation. I found that Br atoms are creating this error of maximum force. Please see the PDB file attached herewith and provide suggestions on how to rectify it.
output.dat (3.3 MB)

9

It is the Br- itself, so as I noted above, the crude approach to producing its ion topology is inadequate. You need to either properly parametrize it yourself, find a force field that already has suitable Br- parameters, or use a different anion.

10

So, can I replace the Br- ion parameters on amber99sb-ILDN (creating error) with Br- ion parameters from amber99sb forcefield or OPLS forcefield?
Do you think it can resolve the problem? Is the above approach, by the way, correct?

11

Hi jalemkul, I replaced Br- ion parameters (present in ffnonbonded.itp) of amber99sb-ILDN with the parameters present in amber14sb_OL15.ff (https://ftp.gromacs.org/contrib/forcefields/amber14sb_OL15.ff_corrected-Na-cation-params.tar.gz)and its working fine now.
The parameters of Br- ion which were creating error was
[ atomtypes ]
; name at.num mass charge ptype sigma epsilon
Br 35 79.90 0.0000 A 0.00000e+00 0.00000e+00

now, I have used new Br- parameters from amber14sb_OL15.ff which are working fine:
[ atomtypes ]
; name at.num mass charge ptype sigma epsilon
Br 35 79.90 0.0000 A 4.64693e-01 2.45414e-01

I just want to know if it working fine. Can I proceed this way?
and is this correct approach?

12

If these are established/validated parameters, sure. The use of zero for both ε and σ is obviously unphysical and was the cause of your problem. Without suitable LJ terms, the charges collapse on top of each other and lead to electrostatic singularities.

13

Thank you so much for clarification.