[BUG] Index error from AssignLeaflets with the AMBER force field

[pablo-arantes]

Describe the bug
Get an Index Error when running AssignLeaflets with a system built with AMBER forcefield. The membrane contains POPC.

To Reproduce
A minimal working example of code to reproduce the unexpected behaviour.

import MDAnalysis as mda
from lipyphilic.lib.assign_leaflets import AssignLeaflets
from lipyphilic.lib.area_per_lipid import AreaPerLipid

u = mda.Universe("system_amber_nw.prmtop", "prod_wat_memb1-1_nw.dcd")

leaflets = AssignLeaflets(
  universe=u,
  lipid_sel = 'resname POP and name P O11 O12 O13 O14'
)
leaflets.run()

areas = AreaPerLipid(
  universe=u,
  lipid_sel='resname POP and name P O11 O12 O13 O14',
  leaflets=leaflets.leaflets
)
areas.run()

Expected behaviour
Run smoothly and store results in leaflets.leaflets attribute

Actual behaviour
Get an Index Error

IndexError                                Traceback (most recent call last)
[<ipython-input-32-3913732d68c0>](https://localhost:8080/#) in <module>
     14   lipid_sel = "name P"
     15 )
---> 16 leaflets.run()
     17 
     18 areas = AreaPerLipid(

2 frames
[/usr/local/lib/python3.8/site-packages/lipyphilic/lib/assign_leaflets.py](https://localhost:8080/#) in _assign_leaflets(self, memb_midpoint_xy)
    432         upper_leaflet = self.membrane[
    433             self.membrane.positions[:, 2] >
--> 434             (memb_midpoint_xy.statistic[lipid_x_bins, lipid_y_bins])  # we don't to consider midplane_cutoff here
    435         ]
    436         self.leaflets[

IndexError: index 2 is out of bounds for axis 1 with size 2

Additional context

• LiPyphilic 0.10.0
• Python 3.8
• Ubuntu distribution of Linux

[p-j-smith]

Hi @pablo-arantes sorry for not getting back to you sooner. That is quite strange - it's not clear where exactly the error is happening but it looks like self.membrane.positions might not have any z positions for some reason.

Would you be able to share your prmtop and a small part of your dcd file and I'll take a look for you?

[kaistroh]

Hi @p-j-smith and @pablo-arantes, I had encountered the same IndexError a while ago (still in version 0.9.0, but as far as I can tell this part is unchanged in 0.10.0). The problem appeared always when having different x and y dimensions, i.e., AssignLeaflets works fine on square membranes, and fails otherwise.

At least in my case, the problem is not caused by missing z positions in self.membrane.positions. In my test case the y dimension is roughly twice as long as the x dimension. Then the IndexError is caused by lipid_y_bins having some bin indices = 2 and memb_midpoint_xy.shape being (2,2). Apparently, the incorrect bin indices are caused by usage of the same bins for x and y in scipy.stats.binned_statistic_2d.

The problem is solved by

        x_bins = memb_midpoint_xy.x_edge
        y_bins = memb_midpoint_xy.y_edge
     
        # get the binnumbers for each lipid
        lipid_x_bins, lipid_y_bins = scipy.stats.binned_statistic_2d(
            x=self.membrane.positions[:, 0],
            y=self.membrane.positions[:, 1],
            values=self.membrane.positions[:, 2],
            statistic="mean",
            bins=(x_bins, y_bins),
            expand_binnumbers=True
        ).binnumber -1  # These were bin numbers, now bin indices  # noqa: E225

and similar changes in the other calls to scipy.stats.binned_statistic_2d

[p-j-smith]

Hi @kaistroh, thanks so much for the info on the issue!

I think there's still a separate issue with the AMBER force field. Diacyl lipids in the AMBER force field are composed of 3 residues - one for the headgroup and one for each tail. However, all analyses in lipyphilic assume 1 residue per lipid (as is the case in CHARMM), which I think that's the cause of the original error reported by @pablo-arantes.

I've opened a separate issue (#132) and fixed (#133) for the problem you've reported @kaistroh, but fixing the AMBER issue will require a more substantial re-write of lipyphilic