Error in GROMACS for adding ions to a solution

Question

I am following Justin Lemkul's GROMACS simulation given on this website here for GROMACS.

I believe I have followed the process to the T. I would like to mention that in the "Generate Topology" section of the tutorial, he says to include

#ifdef POSRES_B
#include "posre_Protein_chain_B.itp"
#endif

However, this was already present in my posre_Protein_chain_B.itp file. However this statement was also present in my posre_Protein_chain_A.itp, posre_Protein_chain_C.itp, posre_Protein_chain_D.itp, posre_Protein_chain_E.itp files. For example, posre_Protein_chain_A.itp had

#ifdef POSRES_A
#include "posre_Protein_chain_A.itp"
#endif

in it. I have commented them out, so it now looks like:

; #ifdef POSRES_A
; #include "posre_Protein_chain_A.itp"
; #endif

Furthermore, every time I was asked to choose for N-terminus, I hit option 2: None, and everytime I was asked to choose a C-terminus, I hit option 0: COO-. This resulted in the system having a net negative charge I believe.

I don't think I have deviated anywhere else in this simulation. However, I am getting the error message in the attachment when I run

gmx grompp -f ions.mdp -c solv.gro -p topol.top -o ions.tpr

After running this command I received the following error message:

GROMACS is free software; you can redistribute it and/or modify it
under the terms of the GNU Lesser General Public License
as published by the Free Software Foundation; either version 2.1
of the License, or (at your option) any later version.

GROMACS:      gmx grompp, version 2020.5
Executable:   /home/satyend/.local/bin/gmx
Data prefix:  /home/satyend/.local
Working dir:  /home/satyend/software/gromacs_gpu/umbrella-sampling-protofibril
Command line:
  gmx grompp -f ions.mdp -c solv.gro -p topol.top -o ions.tpr

Ignoring obsolete mdp entry 'ns_type'

NOTE 1 [file ions.mdp]:
  With Verlet lists the optimal nstlist is >= 10, with GPUs >= 20. Note
  that with the Verlet scheme, nstlist has no effect on the accuracy of
  your simulation.

Setting the LD random seed to 2139084382
Generated 165 of the 1596 non-bonded parameter combinations
Excluding 3 bonded neighbours molecule type 'Protein_chain_A'
Excluding 3 bonded neighbours molecule type 'Protein_chain_B'
Excluding 3 bonded neighbours molecule type 'Protein_chain_C'
Excluding 3 bonded neighbours molecule type 'Protein_chain_D'
Excluding 3 bonded neighbours molecule type 'Protein_chain_E'
Excluding 2 bonded neighbours molecule type 'SOL'

WARNING 1 [file topol.top, line 54]:
  The GROMOS force fields have been parametrized with a physically
  incorrect multiple-time-stepping scheme for a twin-range cut-off. When
  used with a single-range cut-off (or a correct Trotter
  multiple-time-stepping scheme), physical properties, such as the density,
  might differ from the intended values. Since there are researchers
  actively working on validating GROMOS with modern integrators we have not
  yet removed the GROMOS force fields, but you should be aware of these
  issues and check if molecules in your system are affected before
  proceeding. Further information is available at
  https://redmine.gromacs.org/issues/2884 , and a longer explanation of our
  decision to remove physically incorrect algorithms can be found at
  https://doi.org/10.26434/chemrxiv.11474583.v1 .


NOTE 2 [file topol.top, line 54]:
  System has non-zero total charge: -10.000000
  Total charge should normally be an integer. See
  http://www.gromacs.org/Documentation/Floating_Point_Arithmetic
  for discussion on how close it should be to an integer.
  


Number of degrees of freedom in T-Coupling group rest is 67749.00

NOTE 3 [file ions.mdp]:
  You are using a plain Coulomb cut-off, which might produce artifacts.
  You might want to consider using PME electrostatics.



There were 3 notes

There was 1 warning

-------------------------------------------------------
Program:     gmx grompp, version 2020.5
Source file: src/gromacs/gmxpreprocess/grompp.cpp (line 2353)

Fatal error:
Too many warnings (1).
If you are sure all warnings are harmless, use the -maxwarn option.

For more information and tips for troubleshooting, please check the GROMACS
website at http://www.gromacs.org/Documentation/Errors
-------------------------------------------------------
Analysing residue names:
There are:   135    Protein residues
There are: 10727      Water residues
Analysing Protein...
This run will generate roughly 3 Mb of data

This is the link to the ions.mdp file Dr. Lemkul has provided: http://www.mdtutorials.com/gmx/umbrella/Files/ions.mdp.

What is going wrong here? I feel like I am running around in circles with this tutorial...

I appreciate any advice you have to offer!

Answer 1

As I had previously written in the comments, the warning given by gromacs seems to indicate that the GROMOS force fields should not be used due to wrong parameterization:

  The GROMOS force fields have been parametrized with a physically
  incorrect multiple-time-stepping scheme for a twin-range cut-off. When
  used with a single-range cut-off (or a correct Trotter
  multiple-time-stepping scheme), physical properties, such as the density,
  might differ from the intended values. Since there are researchers
  actively working on validating GROMOS with modern integrators we have not
  yet removed the GROMOS force fields, but you should be aware of these
  issues and check if molecules in your system are affected before
  proceeding. Further information is available at
  https://redmine.gromacs.org/issues/2884 , and a longer explanation of our
  decision to remove physically incorrect algorithms can be found at
  https://doi.org/10.26434/chemrxiv.11474583.v1

(copied from the question)

The redmine link returns a 404 error, but the chemrxiv paper is accessible. From a quick read of the paper, it seems that the GROMACS code inherited its algorithms from the GROMOS87 code. The old algorithm used a "twin range cut-off" scheme, where the short-range forces were calculated in every step, but the long range forces were calculated every N number of steps (N is a certain integer number, probably set by the software). The long-range forces were added to the short-range values every step to approximate a longer cut-off value. This algorithm was written at a time when computational power was not cheap so not calculating long range forces in every step saved a lot of cpu time. This algorithm has been found to cause artifacts especially for lipid bilayers with zwitterionic molecules. Another problem with the algorithm was that it was not time-reversible.^[1]

A better algorithm is the Trotter multiple time stepping scheme which was adopted by GROMACS since version 4.5 (2013). However, the old GROMOS force fields were parameterized using the legacy GROMOS algorithm, which means that with the new Trotter algorithm, the GROMOS parameters would give different (and possibly wrong) results. There also seems to be some other issues with the constraints (SHAKE) for the GROMOS parameters.^[1] This is why from version 2019.2, GROMACS (the grompp module) prints out a warning whenver the GROMOS force-fields are used, and stops execution.

The only utility of the warning is to ask the user to check whether the system they are using gives correct results with the GROMOS parameters when using the newer algorithm. So, if you get the error, you should check your system before proceeding. Once you are happy with the choice of force-field and system, then just add the argument -maxwarn 1 to the end of the gmx grompp command that you are attempting to run. The command will now run as normal, and ignore the warning. It's as simple as that.

In 2020, Diem and Oostenbrink studied the differences in solvation energy of various amino acid side chains in water and the physical properties of six pure liquids, when using the twin-range algorithm vs the Trotter algorithm. They found that the error in using the old algorithm is mostly minor and within the range of uncertainty in the experimental data used to validate the results. Another important finding was that the error is smaller when using group-based cut-off schemes than for an atom-based cut-off.^[2]

References:

1. B. Hess, D. van der Spoel, M. J. Abraham, E. Lindahl, "On The Importance of Accurate Algorithms for Reliable Molecular Dynamics Simulations", https://doi.org/10.26434/chemrxiv.11474583.v1
2. M. Diem, C. Oostenbrink, "The Effect of Using a Twin-Range Cutoff Scheme for Nonbonded Interactions: Implications for Force-Field Parametrization?", J. Chem. Theory Comput., 2020, 16, 10, 5985–5990, https://pubs.acs.org/doi/10.1021/acs.jctc.0c00509