Free energy, annihilation: charged ligand

User discussions
1

GROMACS version: Any version
GROMACS modification: No

Dear Everyone,

it may be a naif question, but I am going to ask: what is exactly going on in the code when switching off charges and vdWs of a charged chromophore? I am wondering how PME is accounted for, as actually a, say, +1 charge is fading without any buffer ion appearing in its place. I am not getting any warning (at any of the lambda windows: neither the extremes nor the intermediate ones), but simply I wonder how. Starting from a neutralized box with the necessary counterions, all to me is just the same as the initial state (be it lambda = 0), but in the final lambda = 1 state, where the only thing that faded is the +1 charge. Thank you for your time!

Bests,
Jacopo

2

There should be a warning if either of the end states are net charged. If that is not the case, could you open an issue at Issues ¡ GROMACS / GROMACS ¡ GitLab? Preferably with input files to reproduce it?

3

Thank you so much once again. Took me a while to answer as I wanted to furtherly verify once again, as well as prepare simple examples with reproducibility (made a script). Moreover, this very question came from my supervisor, he is opening the gitlab thread and we are going to provide three very simple examples soon

4

Good evening Everyone,

sorry if I am re-opening this thread so late, but in the end we preferred sharing the (possible) bug here. Specifically, it is a matter of using the free energy code when adopting the mdp file keywords:

couple-moltype           = 'choose_solute_to_decouple' 
couple-lambda0           = vdw-q   
couple-lambda1           = none

and corresponding non-null coul_lambdas series, the whole of this with a charged molecule to be decoupled from the environment.

I understand the large advantage of this approach with respect to having to prepare a double topology .top file with A and B states, which is why we were using it when willing to solely compute the solvation energy. But you will see no warning regarding net charge is returned in none of the 3 examples prepared for you (please see attached link): in each one you can find the straightforward script executing the ‘grompp’ and ‘mdrun’ commands.

To face the problem in the meanwhile, I ask whether one could add for instance an ion of the box to the solute (the one to be decoupled) topology, so to turn both off contemporarily; thus the ion turning into a free dummy atom with null intermolecular interactions and surrounded by the solvent. Or if this is completely unreasonable, and one should instead do the proper morphing of an ion to a water (or vice-versa), having to go back to the double topology picture (A and B states with all the corresponding parameters in one top file).

no_error_if_annihilating_charged_solute

Thank you in advance for all your work and explanations.

Bests,
Jacopo

5

The GROMACS machinery handles that case of changes the net charge correctly. When using PME you get the contribution of the background charge added. But in practice this is usually not want you want to do.

You zip file does not contain the topol.top files. But if I try myself I get a note and a warning with a system with net charge zero in the A state and non-zero in the B state:

NOTE 4 [file topol.top, line 54]:
State B has non-zero total charge: 1200.000000
Total charge should normally be an integer. See
Floating point arithmetic - GROMACS 2024.3 documentation
for discussion on how close it should be to an integer.

WARNING 1 [file topol.top, line 54]:
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.

6

Dear Professor Hess,

Thank you for your interest in this topic as well. It is me not being as clear as it should be necessary, as in each of the three attached examples there is an ‘execute’ folder containing a short ‘do.sh’ script to be executed. Via that script, a topology is produced and updated so to properly run the very last grompp command before the mdrun. At the point of the grompp command, no warning is raised, as well as no notes related to this very topic. Nonetheless, “switching off” a charged molecule/particle should raise the need in all intermediate states up to the ending (B) state included.

I hope it is just an issue of mine, but I found both gromacs 22.3 and 24.2 doing the same thing on our machines. I sincerely appreciate your time, wish you a good work.

Bests,
Jacopo

7

Ah, I didn’t use the decouple option. With the decouple option the charge summing is done right before applying the decoupling instead of after. I will fix this and then the note and warning that I pasted will appear.

8

I filed an issue an uploaded a fix:

1 Like