While running a protein-ligand MD, I neutralize the atoms (this protein has a qtot of 18) and gromacs adds 18 CL ions.
However, when I go create the binary input, it tells me:
“NOTE 2 [file topol.top, line 300089]:
System has non-zero total charge: 54.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.”
When I check the line mentioned on topol.top, it just shows me:
“CL 18”
I have used this protein before without issues, so I’m wondering if the problem is with my ligand. And, if so, is there anything I can do to keep the system going without running into the issue of my charge being triplicated?
I guess you got the note in preprocessing using gmx grompp and before you have added the ions to the system using gmx genion.
Note, if you have run more than one time gmx genion with the option -pq 18, then everytime CL ions are added to the topol.top file. Then you should see at the end of the topol.top file several lines “CL 18”.
Or did you use the option -neutral together this gmx genion?
No, I actually got the note after adding the ions to the system.
I added the ions using “gmx genion -s ions.tpr -o solv_ions.gro -p topol.top -pname NA -nname CL -neutral”, and I get the note that 18 CL ions have been added. In the topol file, it indeed shows these 18 CL ions.
Then, when I create the binary input (gmx grompp) to run energy minimization, that’s when I get the note.
So something seems to happen between these two moments.
I have tried adding 18 instead of neutral, and that results in the same issue.
Hi,
yes, I did just that. These are the notes I get:
"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 -441502053
Generated 167799 of the 167910 non-bonded parameter combinations
Generating 1-4 interactions: fudge = 1
Generated 117519 of the 167910 1-4 parameter combinations
Excluding 3 bonded neighbours molecule type ‘Protein’
Excluding 3 bonded neighbours molecule type ‘RTP’
Excluding 2 bonded neighbours molecule type ‘SOL’
Analysing residue names:
There are: 1970 Protein residues
There are: 1 Other residues
There are: 206158 Water residues
Analysing Protein…
Analysing residues not classified as Protein/DNA/RNA/Water and splitting into groups…
Number of degrees of freedom in T-Coupling group rest is 1332279.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.
Hi!
Yes, this is my output when I used “gmx grompp -f ions.mdp -c solv.gro -p topol.top -o ions.tpr” to generate ion.tpr. Was that not what Alessandra was asking for?
After doing this, I used genion to neutralize the system and my topol file indicated it had added 18 Cl ions, as it should. However, when I go to create the binary input using gmx grompp, it indicates that my system has a total charge of 54.
It sounds like something has gone wrong in your workflow, since 54 = 18 * 3 perhaps you’ve overwritten or used the same file multiple times. I don’t have any other explanation for such weird behavior. genion should be very straightforward. Once you add the ions, that’s it. Adding negative ions should not result in an increase in positive charge.
Yeah, I considered that too and have already started this from scratch.
Interestingly, if I add a random amount of ions (I tried 6 Cl), it doesn’t multiply by 3, instead it multiplied by 5.
Anyway, I will try again from scratch, maybe in a different computer.
Thank you!
This seems very weird. It would be helpful for us to see your exact genion command line, and all subsequent screen output (including your selection for replacement). That might provide some clues.
Reading file ions.tpr, VERSION 2020.1-Ubuntu-2020.1-1 (single precision)
Reading file ions.tpr, VERSION 2020.1-Ubuntu-2020.1-1 (single precision)
Will try to add 0 NA ions and 18 CL ions.
Select a continuous group of solvent molecules
Group 0 ( System) has 650252 elements
Group 1 ( Protein) has 31729 elements
Group 2 ( Protein-H) has 15846 elements
Group 3 ( C-alpha) has 1970 elements
Group 4 ( Backbone) has 5910 elements
Group 5 ( MainChain) has 7879 elements
Group 6 ( MainChain+Cb) has 9742 elements
Group 7 ( MainChain+H) has 9762 elements
Group 8 ( SideChain) has 21967 elements
Group 9 ( SideChain-H) has 7967 elements
Group 10 ( Prot-Masses) has 31729 elements
Group 11 ( non-Protein) has 618523 elements
Group 12 ( Other) has 49 elements
Group 13 ( RTP) has 49 elements
Group 14 ( Water) has 618474 elements
Group 15 ( SOL) has 618474 elements
Group 16 ( non-Water) has 31778 elements
Select a group: 15 Selected 15: ‘SOL’
Number of (3-atomic) solvent molecules: 206158
Processing topology
Replacing 18 solute molecules in topology file (topol.top) by 0 NA and 18 CL ions.
Back Off! I just backed up topol.top to ./#topol.top.10#
Using random seed -1274328115.
Replacing solvent molecule 9358 (atom 59852) with CL
Replacing solvent molecule 170645 (atom 543713) with CL
Replacing solvent molecule 54044 (atom 193910) with CL
Replacing solvent molecule 6615 (atom 51623) with CL
Replacing solvent molecule 173893 (atom 553457) with CL
Replacing solvent molecule 61988 (atom 217742) with CL
Replacing solvent molecule 77218 (atom 263432) with CL
Replacing solvent molecule 80682 (atom 273824) with CL
Replacing solvent molecule 60789 (atom 214145) with CL
Replacing solvent molecule 65490 (atom 228248) with CL
Replacing solvent molecule 143093 (atom 461057) with CL
Replacing solvent molecule 187150 (atom 593228) with CL
Replacing solvent molecule 172664 (atom 549770) with CL
Replacing solvent molecule 15696 (atom 78866) with CL
Replacing solvent molecule 24940 (atom 106598) with CL
Replacing solvent molecule 118568 (atom 387482) with CL
Replacing solvent molecule 102741 (atom 340001) with CL
Replacing solvent molecule 180970 (atom 574688) with CL
NOTE 1 [file em.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 744435077
Generated 167799 of the 167910 non-bonded parameter combinations
Generating 1-4 interactions: fudge = 1
Generated 117519 of the 167910 1-4 parameter combinations
Excluding 3 bonded neighbours molecule type ‘Protein’
Excluding 3 bonded neighbours molecule type ‘RTP’
Excluding 2 bonded neighbours molecule type ‘SOL’
Excluding 3 bonded neighbours molecule type ‘CL’
WARNING 1 [file topol.top, line 300089]:
You are using Ewald electrostatics in a system with net charge. This can
lead to severe artifacts, such as ions moving into regions with low
dielectric, due to the uniform background charge. We suggest to
neutralize your system with counter ions, possibly in combination with a
physiological salt concentration.
My topol file:
[ molecules ]
; Compound #mols
Protein 1
rtp 1
SOL 209153
CL 18
This is weird. What force field are you using? It implies that your ion topology has nrexcl set to 3, which is odd and suggests something funny is being done in the assignment of CL properties.
Charmm36 (july 2021 version) and used cgenff for ligand topology, pretty standard…
Is there perhaps a naming issue with the CL among topol, solv_ions and the forcefield? I haven’t gotten any errors in that regard (such as atom names not matching).
OK, we’ve got a stock value of nrexcl there, so it’s not harmful but we will fix that.
The problem I’m seeing is that “CL” is not correct for CHARMM36 as of that release. The Cl- anion is named CLA (and Na+ is SOD, etc) in accordance with the standard CHARMM nomenclature. Does it work correctly if you use the right name when adding the ions? I’m not entirely sure why this ever worked, since “CL” is not a valid molecule name in the force field.
I have checked and in the past our lab had been using the feb2021 version of CHARMM36, for which CL actually worked. We didn’t change the nomenclature in the protocol when we “upgraded” to the july 2021 version. So that’s how the issue arised. I guess by using CL with the july 2021 version we were actually adding calcium ions (not sure why it triplicated instead of doubling, but anyway).
I am very troubled as to why this even works. grompp should exit with a fatal error of “no such moleculetype CL” because there is no [ CL ] defined anywhere in that port, unless someone has added them to your files.