attr aren't added correctly to existing AtomGroups.

[kain88-de]

Expected behaviour

If I create/add occupancies to a topology they should also be written for any format that supports them

Actual behaviour

MDAnalysis claims that no occupancy is available and uses the default value

Code to reproduce the behaviour

gist

import MDAnalysis as mda
from MDAnalysisTests.datafiles import PSF, DCD

u = mda.Universe(PSF, DCD)

u.atoms.occupancies = 5 * np.ones(u.atoms.n_atoms)

u.atoms.write('test-mda.pdb')

[richardjgowers]

So there's two things going on here. Firstly doing u.atoms.occupancies = doesn't add occupancies to the topology, it just adds them to that particular AtomGroup. We should probably add some setattr magic so assigning to the master AG is syntactic sugar for doing it the long way.

And the reason why it doesn't work even though you've attached occupancies to that AG is there's an AG.atoms call inside PDBWriter which converts the AtomGroup you pass into another "identical" one, (but without all the "unofficial" attrs you attached to it.

[kain88-de]

Firstly doing u.atoms.occupancies = doesn't add occupancies to the topology, it just adds them to that particular AtomGroup.

So setting them to like NaN for atoms in the topology but not in the AtomGroup where I added the occupancies.

AG is there's an AG.atoms call inside PDBWriter which converts the AtomGroup you pass into another "identical" one, (but without all the "unofficial" attrs you attached to it.

So now everytime I call ag.atoms I loose any additional attributes I've added to that AtomGroup? Was this also the behavior in 0.15.0?

[richardjgowers]

Not sure what you mean with NaNs. ag.occupancies = here doesn't touch u._topology, it just adds a .occupancies attribute to that particular instance of AtomGroup object.

[kain88-de]

OK below are some more examples of edge cases where I have questions about the way attributes are propagated between a topology and it's atomgroups. That should clear up what I mean with using NaN values.

# contains three segments A, B, C. 100 residues each but no occupancies
u = mda.Universe(...)
A = u.select_atoms('segid A')
B = u.select_atoms('segid B')
C = u.select_atoms('segid C')

# here we make an update that a user might expect to change the topology as well.
A.occupancies = 2 * np.ones(A.n_atoms)
B.occupancies = 3 * np.ones(B.n_atoms)

# Case 1
########

# What occupancies should be written now? I would expect that segment A
# all have value 2 and segment B has 3 while in C we use the default.
u.atoms.write('all.pdb')

# Case  2
#########

# what should this return?
u.atoms.occupancies
# I would assume something like [2, 2, ..., 3, 3, ..., NaN, NaN, ...]

# Case 3
########

# cut in the middle of segment A and B together
AB = u.select_atoms('(segid A and resid 50:100) or (segid B and resid 1:50)')
# What should be returned here?
AB.occupancies

[kain88-de]

I know currently this example wouldn't work because I don't touch the underlying topology. But having this would be great for systems building and analysis. Writing occupancies is an easy way to visualize per atom behavior of a variable in VMD.

[kain88-de]

A current workaround would be

u = mda.Universe(...)
u.atoms.write('temp.pdb')
u = mda.Universe('temp.pdb')

Then I get occupancies and I can write them. We can keep this and say that we don't allow arbitrary conversion between topology formats and adding of attributes during the conversion. But then we should throw an error if we want to add a attribute to a AtomGroup that is not supported by the underlying topology.

[richardjgowers]

So I'm not sure how we should deal with missing values, so if I can slightly alter your cases..

u = mda.Universe()  # without occupancies

ag = u.atoms
ag.colour = 'blue'  # should fail, can't give random attributes to AtomGroup instance

# should again fail as I don't know what topologyattr to use
ag.random_properties = np.ones(len(atoms))  


# Here the 'occupancies' name should get picked up in a list of special cases and a Occupancies TopologyAttr created and added to u._topology
ag.occupancies = np.ones(len(atoms))

Which is tl;dr all attributes on an AG should get values from either u._topology or u.trajectory. The only information in an AG should be a list of indices.

With missing values, maybe we can specify a default value for each, (eg default charge is 0.0). Then partially setting can work as all other values get the default. So I think this agrees with the ideas in your examples.

[orbeckst]

Which is tl;dr all attributes on an AG should get values from either u._topology or u.trajectory. The only information in an AG should be a list of indices.

Sounds like a sane approach; can you use __slots__ to cleanly select what attributes are available?

[richardjgowers]

__slots__ would mean we can't alter the class after the initial definition, which we need to. We want something like slots where we strictly control the __dict__ of the class.

[kain88-de]

With missing values, maybe we can specify a default value for each, (eg default charge is 0.0). Then partially setting can work as all other values get the default. So I think this agrees with the ideas in your examples.

I actually liked the NaN because then I know for sure which values I set and which were set automatically by MDAnalysis. For example if I set myself a 0 for charge I won't be able to distinguish it from the value MDAnalysis set.

[richardjgowers]

So as a model for how this could work....

import numpy as np


def add_occupancies(cls, val):
    # cls - instance of the class
    # val - value to be added
    super(A, cls).__setattr__('occupancies', np.array([val]))


class A(object):
    _ALLOWED = {'a', 'b'}  # set of known attributes
    _SPECIAL_CASES = {'occupancies':add_occupancies}  # special cases with handling functions

    def __init__(self, a, b):
        self.a = a
        self.b = b

    def add_attr(self, attr, val):
        # add this attr to classes allowed list
        self.__class__._ALLOWED.add(attr)
        setattr(self, attr, val)

    def __setattr__(self, attr, val):
        # If attr not in existing 
        if not attr in self.__class__._ALLOWED:
            if attr in self.__class__._SPECIAL_CASES:
                self._SPECIAL_CASES[attr](self, val)
            else:
                raise AttributeError
        else:
            super(A, self).__setattr__(attr, val)

So all setattr go through our __setattr__ method. We maintain something like __slots__, here called _ALLOWED to disallow random attributes to be set.

add_attr is a method that provides a way to add arbitrary things (if we need that).

I just need to see if this will work with our existing TopologyAttr things

[kain88-de]

Looks reasonable besides that the attributes in _ALLOWED can't be changed in this example.

[orbeckst]

With 0.16.2 #1401 likely being the last release for a while (unless we do parallel development for the summer) it would be really good to have this fixed. There are some special cases like adding b-factors #1359 or adding charges to create PQR files that are not working any more.

What needs to happen to accelerate a solution for this issue?

[xiki-tempula]

Kind of curious that if unable to set the residue name is part of this problem?

import MDAnalysis as mda
u = mda.Universe('ZR.gro')
print(u.atoms)
#<AtomGroup [<Atom 1: N of type N of resname ZR, resid 472 and segid SYSTEM>]>
u.resname = 'ZA'
u.atoms.resname = 'ZA'
print(u.atoms)
#<AtomGroup [<Atom 1: N of type N of resname ZR, resid 472 and segid SYSTEM>]>

[richardjgowers]

Ok, so with #1186 merged the following now works:

# stuff that isn't correct raises a proper error rather than silently proceeding 
u.atoms.occupancies = 5 * np.ones(u.atoms.n_atoms)
AttributeError: Cannot set arbitrary attributes to a Group

# Topology Attributes can be explicitly added to make them work
u.add_TopologyAttr('occupancies')
u.atoms.occupancies = 5 * np.ones(u.atoms.n_atoms)

# OR
u.add_TopologyAttr('occupancies', values=5 * np.ones(u.atoms.n_atoms))

The special cases proposed here where the error is turned into the correct way to add things hasn't been implemented as we're not catching errors and scanning against known Attributes