Harvesting only electrostatic force from rerun

User discussions
1

GROMACS version: 2023-dev-20220412
GROMACS modification: No

I have a question that is similar to other questions asked before (L-J and electrostatic force - #3 by di_jin and Coulomb force calculation), but wanted to see if there was a better way of doing what I am trying right now. The gist of this is that I want to get only the electrostatic forces/energies out of the system, and not have any other forces/energies calculated in this output.

I have a lipid bilayer system with a peptide, and I have reason to believe that the electrostatic interaction of the peptide with its environment is the dominant source of peptide motion in my system. I have large XTC trajectory files, and so can use the -rerun option to get the forces back out of the system by post-processing my files. I am attempting to get out just the electrostatic force on each atom of the peptide from all atoms in the system (I realize that the force between any two specific groups is fictitious due to not including all of the charges/screening in the system). What I’m trying to do is avoid writing my own electrostatic force computation in python, as this calculation should already be done in GROMACS underneath the surface, but I cannot figure out how to only recalculate the electrostatic forces/energies and get them out of gromacs.

I’ve tried using energy groups, different coulombic cutoffs, and free energy groups, but nothing seems to work at effectively only getting the electrostatic force out of the system! In the end, I just want the force on each peptide atom as a function of time, but only the contribution from the electrostatic calculations. I would have thought that using a free-energy based calculation would have worked, but my setup doesn’t seem to handle this, or I am doing something wrong with my lambda states/values.

Any help on getting the electrostatic force only on my peptide would be greatly appreciated.

2

Use gmx convert-tpr -zeroq to write a new .tpr file with all charges set to zero. Then use mdrun -rerun to re-evaluate energies and forces, and the difference between the original and new energies and forces corresponds to that of the electrostatic terms.

3

-zeroq has been removed since version 2022.2, and had an incorrect behavior since around 2008 (issue).

If short-range energies are sufficient (no long-range, no forces), then rerun with energygrps should work.

Otherwise, the most straightforward way is probably editing the topology to remove all charges for non-peptide atoms then performing a rerun (copy ions.itp and tip3p.itp or equivalent to local directory from FF directory, edit the relevant charges to 0.0, change the #includes in topol.top, grompp then rerun with xtc).

1 Like
4

Thanks! This definitely worked, and I am now able to separate my forces on the peptide. Thank you for your help!

5

I guess the easiest, slightly hacky, way of doing this would be to make epsilon_r very small, e.g. 1e-6. Then you need to multiply the output forces by 1e-6 and non-electrostatic contributions will be negligible. I have not tried this, but I guess it should work with rerun in release-2022, as this does not use the forces for updating the coordinates.