# How to calculate residue‑membrane contact profile (residue index vs contact count) in GROMACS‑2018 +

**URL:** <https://gromacs.bioexcel.eu/t/how-to-calculate-residue-membrane-contact-profile-residue-index-vs-contact-count-in-gromacs-2018/13777>\
**Category:** User discussions\
**Created:** [October 18, 2026, 10:17am UTC](https://gromacs.bioexcel.eu/t/how-to-calculate-residue-membrane-contact-profile-residue-index-vs-contact-count-in-gromacs-2018/13777 "2026-10-18T10:17:34Z")\
**Posts on this page:** 3\
**Page:** 1

<div class="post-metadata">

**Author:** ![ZXY](https://avatars.discourse-cdn.com/v4/letter/z/f4b2a3/32.png) [@ZXY](https://gromacs.bioexcel.eu/u/ZXY)\
**Post date:** [October 18, 2026, 10:17am UTC](https://gromacs.bioexcel.eu/t/how-to-calculate-residue-membrane-contact-profile-residue-index-vs-contact-count-in-gromacs-2018/13777/1 "2026-10-18T10:17:35Z")

</div>

GROMACS version:  
GROMACS modification: Yes/No  
Here post your question  
Hello everyone,

I am working on coarse‑grained MD simulations (Martini3 force field, GROMACS‑2018).  
I want to compute a residue‑membrane contact profile plot:

- X‑axis: residue index of protein (1‑185 for K‑Ras4B)
- Y‑axis: contact count / contact probability between each protein residue and inner‑leaflet lipid head‑groups.

Although `gmx mindist` can output results on a per‑residue basis, it only outputs minimum distance values instead of contact counts.

---

<div class="post-metadata">

**Author:** ![fatpmeireles](https://avatars.discourse-cdn.com/v4/letter/f/ecd19e/32.png) [@fatpmeireles](https://gromacs.bioexcel.eu/u/fatpmeireles)\
**Post date:** [October 20, 2026, 8:13am UTC](https://gromacs.bioexcel.eu/t/how-to-calculate-residue-membrane-contact-profile-residue-index-vs-contact-count-in-gromacs-2018/13777/2 "2026-10-20T08:13:22Z")

</div>

Hi,

`gmx mindist` does have an option to output contact counts (`-on`). If you tried that and, for some reason, it did not work, then I recommend checking out the PyLipID library for protein-membrane interaction analysis; it also does exactly what you want: [0-Application walk thourgh - PyLipID documentation](https://pylipid.readthedocs.io/en/master/tutorials/0-application-walk-through.html)

Cheers!

---

<div class="post-metadata">

**Author:** ![Erik](https://avatars.discourse-cdn.com/v4/letter/e/b5ac83/32.png) [@Erik](https://gromacs.bioexcel.eu/u/Erik)\
**Post date:** [October 20, 2026, 8:52pm UTC](https://gromacs.bioexcel.eu/t/how-to-calculate-residue-membrane-contact-profile-residue-index-vs-contact-count-in-gromacs-2018/13777/3 "2026-10-20T20:52:29Z")

</div>

I ended up having to use MDAnalysis for this task because there’s no Gromacs tool which provides information about each contact. Here’s some code which might help:

```python
import MDAnalysis as mda
import numpy as np

# Define universe, selection groups, and a list of headgroup atom ids
u = mda.Universe(top, traj)        
protein_atoms = u.select_atoms('<some select string>')
head_atoms = u.select_atoms('<some select string>'
dist_to_head_aid = np.array([a.id for a in head_atoms], dtype=int)

# Compute contacts
cutoff = 3.5
print("Calculating distances...")
all_distances = np.empty((len(protein_atoms), len(head_atoms)))
for i, frame in enumerate(u.trajectory):
    # Get all distances between protein and lipid heads
    # This will parallelize over however many processes are available
    distances.distance_array(protein_atoms, head_atoms, box=frame.dimensions, result=all_distances, backend="OpenMP")

    # Get atom IDs for all atoms where the distance array is less than the cutoff
    indices = np.nonzero(all_distances<cutoff)
    for j in range(len(indices[0])):
        protein_id = indices[0][j]
        lipid = dist_to_head_aid[indices[1][j]]
        dist = all_distances[(protein_id, indices[1][j])]
        # I used this to build a database of contacts, but you could just accumlate into an array of len(protein_atoms)

```