Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

can't change segmend ids in a loop for regions of a atomgroup #1683

Open
kain88-de opened this issue Oct 20, 2017 · 6 comments
Open

can't change segmend ids in a loop for regions of a atomgroup #1683

kain88-de opened this issue Oct 20, 2017 · 6 comments

Comments

@kain88-de
Copy link
Member

Expected behaviour

I want to iterate over a atomgroup and change the segids for regions in a loop.

Actual behaviour

Our getitem methods doesn't seem to honor the range I specify. So here when i=1
all atoms from 0-428 are getting the segid B.

Code to reproduce the behaviour

>>> import MDAnalysis as mda
>>> large_base_u = mda.fetch_mmtf('1bcc')
>>> n_copies = 5
>>> want_n_atoms = 214 * n_copies
>>> copy = large_base_u.select_atoms('name CA')[:want_n_atoms]
>>> letters = 'ABCDEFGH'
>>> for i in range(n_copies):
>>>     start = i*(u.atoms.n_atoms)
>>>     stop = (i+1)*u.atoms.n_atoms
>>>     print(letters[i], start, stop)
>>>     copy[start:stop].segments.segids = letters[i]
>>>     print(copy[:stop].segids)

>>> np.unique(copy.segids)
array(['C', 'D', 'E'], dtype=object)

Currently version of MDAnalysis:

(run python -c "import MDAnalysis as mda; print(mda.__version__)")
dev version

@kain88-de
Copy link
Member Author

copy[0:214].segments.segids = 'A'
copy[214:428].segments.segids = 'B'

This is shorter. Here the second line is overwriting the atoms from 0-214

Any idea of why this happens. This makes it pretty hard to build my system currently

@kain88-de
Copy link
Member Author

affects 0.16.2 as well

@richardjgowers
Copy link
Member

@kain88-de setting segids won't create new Segments (any more). So I think what's happening is that u.atoms[:214].segments refers to the same segments as u.atoms[214:428].segments, so you're just setting segids to the same segments twice.

If you want to add a Segment to the Universe, do new = u.add_Segment(segid='B'); copy[214:428].segments = new

@kain88-de kain88-de mentioned this issue Oct 21, 2017
6 tasks
@kain88-de
Copy link
Member Author

But wouldn't it be better then if segids are readonly? The simple fact that u.atoms[:214].segments.segids = 'a' doesn't raise an error lead me to believe it's a valid way to update segment ids.

Also you explanation makes sense. It's a bit confusing though that it works only half now. If I change all segids in a atomgroup to the same letter everything seems fine.

@kain88-de kain88-de added the API label Oct 21, 2017
@richardjgowers
Copy link
Member

@kain88-de it is a valid way to update segids. You're only running into trouble because you're updating the segid of the same segment twice. See the toy example below where I deliberately create two ags (ag1 and ag2) which both belong to the same segment (seg)

There's no real way to disallow setting the segid from an AtomGroup like this, because u.atoms.segments is just a SegmentGroup with no memory/knowledge that it came from an AtomGroup.

In [1]: import MDAnalysis as mda

In [2]: u = mda.Universe('adk_oplsaa.gro')
/home/richard/miniconda2/envs/mdadev/lib/python2.7/site-packages/MDAnalysis/topology/guessers.py:56: UserWarning: Failed to guess the mass for the following atom types: M
  "".format(', '.join(misses)))

In [3]: u.segments
Out[3]: <SegmentGroup with 1 segment>

In [4]: ag1 = u.segments.atoms[:100]

In [5]: ag2 = u.segments.atoms[100:]

In [6]: seg = u.segments[0]

In [7]: seg.segid
Out[7]: 'SYSTEM'

In [8]: ag1.segments.segids = 'A'

In [9]: ag1.segids
Out[9]: 
array(['A', 'A', 'A', 'A', 'A', 'A', 'A', 'A', 'A', 'A', 'A', 'A', 'A',
       'A', 'A', 'A', 'A', 'A', 'A', 'A', 'A', 'A', 'A', 'A', 'A', 'A',
       'A', 'A', 'A', 'A', 'A', 'A', 'A', 'A', 'A', 'A', 'A', 'A', 'A',
       'A', 'A', 'A', 'A', 'A', 'A', 'A', 'A', 'A', 'A', 'A', 'A', 'A',
       'A', 'A', 'A', 'A', 'A', 'A', 'A', 'A', 'A', 'A', 'A', 'A', 'A',
       'A', 'A', 'A', 'A', 'A', 'A', 'A', 'A', 'A', 'A', 'A', 'A', 'A',
       'A', 'A', 'A', 'A', 'A', 'A', 'A', 'A', 'A', 'A', 'A', 'A', 'A',
       'A', 'A', 'A', 'A', 'A', 'A', 'A', 'A', 'A'], dtype=object)

In [10]: seg.segid
Out[10]: 'A'

In [11]: ag2.segments.segids = 'B'

In [12]: ag2.segments.segids
Out[12]: array(['B'], dtype=object)

In [13]: ag2.segids
Out[13]: array(['B', 'B', 'B', ..., 'B', 'B', 'B'], dtype=object)

In [14]: seg.segid
Out[14]: 'B'

In [15]: ag1.segids
Out[15]: 
array(['B', 'B', 'B', 'B', 'B', 'B', 'B', 'B', 'B', 'B', 'B', 'B', 'B',
       'B', 'B', 'B', 'B', 'B', 'B', 'B', 'B', 'B', 'B', 'B', 'B', 'B',
       'B', 'B', 'B', 'B', 'B', 'B', 'B', 'B', 'B', 'B', 'B', 'B', 'B',
       'B', 'B', 'B', 'B', 'B', 'B', 'B', 'B', 'B', 'B', 'B', 'B', 'B',
       'B', 'B', 'B', 'B', 'B', 'B', 'B', 'B', 'B', 'B', 'B', 'B', 'B',
       'B', 'B', 'B', 'B', 'B', 'B', 'B', 'B', 'B', 'B', 'B', 'B', 'B',
       'B', 'B', 'B', 'B', 'B', 'B', 'B', 'B', 'B', 'B', 'B', 'B', 'B',
       'B', 'B', 'B', 'B', 'B', 'B', 'B', 'B', 'B'], dtype=object)

@orbeckst
Copy link
Member

For docs we need a section on topology manipulation; relabeling is a common theme.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Projects
None yet
Development

No branches or pull requests

3 participants