# Harvesting only electrostatic force from rerun

**URL:** <https://gromacs.bioexcel.eu/t/harvesting-only-electrostatic-force-from-rerun/4780>\
**Category:** User discussions\
**Tags:** mdrun\
**Created:** [September 13, 2022, 3:57pm UTC](https://gromacs.bioexcel.eu/t/harvesting-only-electrostatic-force-from-rerun/4780 "2022-09-13T15:57:21Z")\
**Posts on this page:** 5\
**Page:** 1

<div class="post-metadata">

**Author:** ![cedelmaier](https://avatars.discourse-cdn.com/v4/letter/c/a698b9/32.png) [@cedelmaier](https://gromacs.bioexcel.eu/u/cedelmaier)\
**Post date:** [September 13, 2022, 3:57pm UTC](https://gromacs.bioexcel.eu/t/harvesting-only-electrostatic-force-from-rerun/4780/1 "2022-09-13T15:57:21Z")

</div>

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](https://gromacs.bioexcel.eu/t/l-j-and-electrostatic-force/2856/3) and [Coulomb force calculation](https://gromacs.bioexcel.eu/t/coulomb-force-calculation/2436)), 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.

---

<div class="post-metadata">

**Author:** ![jalemkul](https://dub1.discourse-cdn.com/flex017/user_avatar/gromacs.bioexcel.eu/jalemkul/32/18_2.png) [@jalemkul](https://gromacs.bioexcel.eu/u/jalemkul)\
**Post date:** [September 15, 2022, 2:50am UTC](https://gromacs.bioexcel.eu/t/harvesting-only-electrostatic-force-from-rerun/4780/2 "2022-09-15T02:50:40Z")

</div>

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.

---

<div class="post-metadata">

**Author:** ![ebriand](https://dub1.discourse-cdn.com/flex017/user_avatar/gromacs.bioexcel.eu/ebriand/32/2314_2.png) [@ebriand](https://gromacs.bioexcel.eu/u/ebriand)\
**Post date:** [September 15, 2022, 8:16am UTC](https://gromacs.bioexcel.eu/t/harvesting-only-electrostatic-force-from-rerun/4780/3 "2022-09-15T08:16:37Z")

</div>

`-zeroq` has been removed since version 2022.2, and had an incorrect behavior since around 2008 ([issue](https://gitlab.com/gromacs/gromacs/-/issues/4226)).

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).

---

<div class="post-metadata">

**Author:** ![cedelmaier](https://avatars.discourse-cdn.com/v4/letter/c/a698b9/32.png) [@cedelmaier](https://gromacs.bioexcel.eu/u/cedelmaier)\
**Post date:** [September 15, 2022, 6:08pm UTC](https://gromacs.bioexcel.eu/t/harvesting-only-electrostatic-force-from-rerun/4780/4 "2022-09-15T18:08:52Z")

</div>

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

---

<div class="post-metadata">

**Author:** ![hess](https://dub1.discourse-cdn.com/flex017/user_avatar/gromacs.bioexcel.eu/hess/32/416_2.png) [@hess](https://gromacs.bioexcel.eu/u/hess)\
**Post date:** [September 16, 2022, 11:46am UTC](https://gromacs.bioexcel.eu/t/harvesting-only-electrostatic-force-from-rerun/4780/5 "2022-09-16T11:46:56Z")

</div>

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.
