GROMOS54a7 ff leads to missing parameters in topology after S-S-bonds between Cys residues were app

User discussions
, ,
1

GROMACS version: gromacs-2021.6
GROMACS modification: No

Hello community!

I am a relatively new GROMACS user, who worked with GROMACS during my bachelor’s thesis - but being supported by a PHD student a few years ago already. Now I am working with GROMACS again to obtain results for my Master Thesis but without the former support, this is why I am reaching out now to the community.

I want to run simulations on a protein called CALB (Candida Antartica Lipase B) using different configurations, such as with and without crystallographic water, different force fields, etc.

A detailed description of what I did will follow below.

My problem in a nutshell:

When I created the top file using pdb2gmx and any GROMOS ff, some entries for bonds, angles, dihedrals, etc. are missing/empty.
These partially empty lines lead to a fatal error when using grompp at later stages, which could be “solved” by using the flag “-zero”.
However, simply filling up missing entries with zeros is not my preferred work around, I think there must be a way to enter the correct parameters.
In contrast, creating a topology with pdb2gmx and AMBER03 results in a complete topology (and not following errors.)

Now the details:

The protein structure has been downloaded form the protein data bank (ID: 1TCA), as we want to compare our results to a paper in which that particular structure has been mentioned. (Link to the original pdb file: 1tca.pdb - Google Drive)

I did two changes to the PDB-structure file:

    1. I removed all entries related to “NAG” (because GROMACS did not recognize this molecule and at my current experience level, I do not know yet how to change this without creating physical non-sense.)
    1. In some cases I also removed “HOH” (crystallographic waters) - this has however no impact on this matter, I am just bringing this up for completeness.
      The way I removed “NAG” and “HOH” from the pdb-file was via a terminal command like the example below:
      grep -v HOH 1tca.pdb > 1_1tca_wo_HOH.pdb
      You can have a look at this file here, on my google drive: 2_1tca_wo_HOH_NAG.pdb - Google Drive

As my third step, I used pdb2gmx to create the topology for GROMOS54a7.

This is my input:
gmx pdb2gmx -f 2_1tca_wo_HOH_NAG.pdb -o 3_1tca_structure.gro -water spce
(using “14” to chose GROMOS54a7 when prompted)

This is what GROMACS returns:

Using the Gromos54a7 force field in directory gromos54a7.ff

going to rename gromos54a7.ff/aminoacids.r2b
Opening force field file /data/skloth/apps/gromacs-2021.6/share/gromacs/top/gromos54a7.ff/aminoacids.r2b
Reading 2_1tca_wo_HOH_NAG.pdb…
WARNING: all CONECT records are ignored
Read ‘LIPASE’, 2324 atoms

Analyzing pdb file
Splitting chemical chains based on TER records or chain id changing.

There are 1 chains and 0 blocks of water and 317 residues with 2324 atoms

chain #res #atoms

  • 1 ‘A’ 317 2324 *

All occupancies are one
All occupancies are one
Opening force field file /data/skloth/apps/gromacs-2021.6/share/gromacs/top/gromos54a7.ff/atomtypes.atp

Reading residue database… (Gromos54a7)
Opening force field file /data/skloth/apps/gromacs-2021.6/share/gromacs/top/gromos54a7.ff/aminoacids.rtp

Using default: not generating all possible dihedrals

Using default: excluding 3 bonded neighbors

Using default: generating 1,4 H–H interactions

Using default: removing proper dihedrals found on the same bond as a proper dihedral

Using default: removing proper dihedrals found on the same bond as a proper dihedral
Opening force field file /data/skloth/apps/gromacs-2021.6/share/gromacs/top/gromos54a7.ff/aminoacids.hdb
Opening force field file /data/skloth/apps/gromacs-2021.6/share/gromacs/top/gromos54a7.ff/aminoacids.n.tdb
Opening force field file /data/skloth/apps/gromacs-2021.6/share/gromacs/top/gromos54a7.ff/aminoacids.c.tdb

Processing chain 1 ‘A’ (2324 atoms, 317 residues)
Analysing hydrogen-bonding network for automated assignment of histidine protonation. 456 donors and 452 acceptors were found.
There are 703 hydrogen bonds
Will use HISD for residue 224

Identified residue LEU1 as a starting terminus.

Identified residue PRO317 as a ending terminus.
8 out of 8 lines of specbond.dat converted successfully
Special Atom Distance matrix:
CYS22 CYS64 MET72 MET83 MET129 CYS216 HIS224
SG152 SG453 SD520 SD611 SD962 SG1595 NE21658
CYS64 SG453 0.203
MET72 SD520 1.630 1.590
MET83 SD611 1.393 1.301 0.936
MET129 SD962 1.972 1.776 1.827 1.450
CYS216 SG1595 3.258 3.070 2.845 2.788 1.462
HIS224 NE21658 2.402 2.284 1.206 1.656 1.489 1.872
CYS258 SG1904 3.341 3.152 3.006 2.933 1.570 0.203 2.052
CYS293 SG2142 3.549 3.542 2.021 2.551 3.487 4.150 2.392
MET298 SD2181 2.461 2.470 1.495 1.459 2.811 3.961 2.416
CYS311 SG2284 3.398 3.393 1.898 2.388 3.375 4.101 2.344
CYS258 CYS293 MET298
SG1904 SG2142 SD2181
CYS293 SG2142 4.347
MET298 SD2181 4.137 1.623
CYS311 SG2284 4.297 0.204 1.420
Linking CYS-22 SG-152 and CYS-64 SG-453…
Linking CYS-216 SG-1595 and CYS-258 SG-1904…
Linking CYS-293 SG-2142 and CYS-311 SG-2284…
Start terminus LEU-1: NH3+
End terminus PRO-317: COO-

Checking for duplicate atoms…

Generating any missing hydrogen atoms and/or adding termini.

Now there are 317 residues with 2930 atoms

Making bonds…

Number of bonds was 2997, now 2992

Generating angles, dihedrals and pairs…

WARNING: WARNING: Residue 1 named LEU of a molecule in the input file was mapped
to an entry in the topology database, but the atom H used in
an interaction of type angle in that entry is not found in the
input file. Perhaps your atom and/or residue naming needs to be
fixed.

WARNING: WARNING: Residue 317 named PRO of a molecule in the input file was mapped
to an entry in the topology database, but the atom O used in
an interaction of type angle in that entry is not found in the
input file. Perhaps your atom and/or residue naming needs to be
fixed.

Before cleaning: 4986 pairs
Before cleaning: 6441 dihedrals

Making cmap torsions…

There are 2267 dihedrals, 1467 impropers, 4377 angles 4986 pairs, 2992 bonds and 0 virtual sites

Total mass 33015.526 a.m.u.

Total charge -1.000 e

Writing topology

Writing coordinate file…

	--------- PLEASE NOTE ------------

You have successfully generated a topology from: 2_1tca_wo_HOH_NAG.pdb.

The Gromos54a7 force field and the spce water model are used.

I uploaded my input and output files to a google drive. You can have a look at the

You will note that the following lines in the topology have missing entries:

[ bonds ] → Always missing c0

  • Line 3469: Atoms 186 (SG in 22Cys) and 562 (SG in 64Cys)
  • Line 5336: Atoms 2013 (SG in 216Cys) and 2415 (SG in 258Cys)
  • Line 6035: Atoms 2700 (SG in 296Cys) and 2881 (SG in 311Cys)

[ angles ] → Always missing c0

  • Line: 11539: Atoms 185 (CB in 22Cys), 186 (SG in 22Cys), 562 (SG in 64Cys)
  • Line 12103: Atoms 186 (SG in 22Cys, 562 (SG in 64Cys), 561 (CB in 64Cys)
  • Line 14869: Atoms 2013 (SG in 216Cys), 2415 (SG in 258Cys), 2414 (CB in 258 Cys)
  • Line 15295: 2699 (CB in 293Cys), 2700 (SG in 296Cys), 2881 (SG in 311Cys)
  • Line 15568: 2700 (SG in 296Cys), 2881 (SG in 311Cys), 2880 (CB in 311Cys)

[ dihedrals ] → Always missing c0

  • Line 15796: Atoms 68 (HD1 in 9Phe), 67 (CD1 in 9Phe), 71 (CE1 in 9Phe), 72 (HE1 in 9Phe)
  • Line 15797: Atoms 68 (HD1 in 9Phe), 67 (CD1 in 9Phe), 71 (CE1 in 9Phe), 75 (CZ in 9Phe)
  • Line 16100: Atoms 220 (N in 27Pro), 221 (CA in 27Pro), 225 (C in 27Pro), 226 (O in 27Pro)
  • Line 17186: Atoms 742 (CG in 82Tyr), 743 (GD1 in 82Tyr), 747 (CE1 in 82Tyr), 751 (CZ in 82Tyr)
  • Line 17187: Atoms 744 (HD1 in 82Tyr), 743 (CD1 in 82Tyr), 747 (CE1 in 82Tyr), 748 (HE1 in 82 Tyr)
  • Line 17726: Atoms 1005 (CB in 110Val), 1004 (CA in 110Val), 1008 (C in 110Val), 1010 (N in 111Ala)
  • Line 17868: Atoms 1087 (C in 117Phe), 1074 (CA in 117Phe), 1075 (CB in 117Phe), 1076 (CG in 117Phe)

Now, if I continue using this topology, the next thing I do is using gmx editconf and gmx solvate to position the protein in a cubic box and solvate it with water (step 4). Just to be complete, this is the input:
gmx editconf -f 3_1tca_structure.gro -o 4_1tca_cubic_box.gro -c -d 1.0 -bt cubic

followed by

gmx solvate -cp 4_1tca_cubic_box.gro -cs spc216.gro -o 4_1tca_cubic_box_solvated.gro -p topol.top

Nothing speical happening here. The fatal error occurs in the next step (step 5).

My input to prepare adding ions to the system using grompp:
(my mdp file: 5_add_ions.mdp - Google Drive)

gmx grompp -f 5_add_ions.mdp -c 4_1tca_cubic_box_solvated.gro -p topol.top -o 5_add_ions.tpr

What GROMACS returns:

GROMACS: gmx grompp, version 2021.6
Executable: /data/skloth/apps/gromacs-2021.6/bin/gmx
Data prefix: /data/skloth/apps/gromacs-2021.6
Working dir: /data/jreusing/Proposal/calb/calb_1tca_pure
Command line:

  • gmx grompp -f 5_add_ions.mdp -c 4_1tca_cubic_box_solvated.gro -p topol.top -o 5_add_ions.tpr*

Ignoring obsolete mdp entry ‘ns_type’

NOTE 1 [file 5_add_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 -167805967

Generated 168 of the 1653 non-bonded parameter combinations

ERROR 1 [file topol.top, line 15796]:

  • No default Proper Dih. types*

ERROR 2 [file topol.top, line 15797]:

  • No default Proper Dih. types*

ERROR 3 [file topol.top, line 16100]:

  • No default Proper Dih. types*

ERROR 4 [file topol.top, line 17186]:

  • No default Proper Dih. types*

ERROR 5 [file topol.top, line 17187]:

  • No default Proper Dih. types*

ERROR 6 [file topol.top, line 17726]:

  • No default Proper Dih. types*

ERROR 7 [file topol.top, line 17868]:

  • No default Proper Dih. types*

Excluding 3 bonded neighbours molecule type ‘Protein_chain_A’

Excluding 2 bonded neighbours molecule type ‘SOL’

WARNING 1 [file topol.top, line 19406]:
he 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 https://redmine.gromacs.org/issues/2884 , and a longer explanation of our decision to remove physically incorrect algorithms can be found at On The Importance of Accurate Algorithms for Reliable Molecular Dynamics Simulations | Theoretical and Computational Chemistry | ChemRxiv | Cambridge Open Engage .

NOTE 2 [file topol.top, line 19406]:
System has non-zero total charge: -1.000000
Total charge should normally be an integer. See Floating point arithmetic — GROMACS webpage https://www.gromacs.org documentation for discussion on how close it should be to an integer.

There were 2 notes

There was 1 warning

Program: gmx grompp, version 2021.6
Source file: src/gromacs/gmxpreprocess/grompp.cpp (line 1963)

Fatal error:
There were 7 errors in input file(s)

For more information and tips for troubleshooting, please check the GROMACS website at Common Errors — GROMACS webpage https://www.gromacs.org documentation

For comparison, these are the files created when using AMBER03 instead: AMBER03 - Google Drive

My ideas:

I think, what happens is that after the S-S-bridge is formed (during which the hydrogens are “used up” to form the bond), GROMACS using the gromos54a7-ff does not know how to create the angles etc. between neighboring residues. Neighboring means close in the 3-D space, not close in the primary sequence.

However, is there any way to fix it?

Side note:

  • I tried all the GROMACS-ff and the result has always been the same.
  • I tried installing older (and later) versions of GROMACS but they didn’t run properly on my computer (the later version works fine but does produce the same errors).

Thank you very much in advance - I really appreciate your time and efforts to help people in this forum!

2

Hi,

I tried your inputs (using the same pdb file) using Gromacs version 2023, and I did not encounter any errors (I had to add -maxwarn 1, but that is unrelated.) Did you try this version as well?

I also see a few missing parameters in the .top file, however, I found that gromos54a7.ff/ffbonded.itp contains the following section:

; bond-, angle- and dihedraltypes for specbonds:
[ bondtypes ]
S      S       2    gb_36
NR     FE      2    gb_34

[ angletypes ]
CH1    CH2    S     2   ga_16
CH2    S      S     2   ga_6
CR1    NR    FE     2	ga_34
NR     FE    NR     2   ga_17

[ dihedraltypes ]
S      S      1   gd_21
NR     FE     1   gd_38
CH2    S      1   gd_26

I imagine that those entries are used to fill in the missing parameters around the disulfide bonds.

I recently posted a related question (I did not see your post before): Missing parameters in cyclic peptide with Gromos FF. It might be that this also answers my question, or maybe I can make an FF file with extra parameter entries to fix it…

I hope this helps.

Best,
Franz