Confusing behavior of ThermodynamicState system attribute

[jchodera]

@andrrizzi : Any idea what is going on here?

Let's check how many torsion terms there are:

>>> from openmmtools import testsystems, states
>>> from simtk import unit, openmm
>>> testsystem = testsystems.HostGuestVacuum()
>>> testsystem.system.getForce(2).getNumTorsions()
758

Now what happens if we create a ThermodynamicState?

>>> thermodynamic_state = states.ThermodynamicState(testsystem.system, temperature=300.0*unit.kelvin)
>>> thermodynamic_state.system.getForce(2).getNumTorsions()
0

Huh? But it gets weirder:

>>> system = thermodynamic_state.system
>>> system.getForce(2).getNumTorsions()
758

Is accessing a method of thermodynamic_state.system somehow causing @setter instead of @property to be called?

[andrrizzi]

Ah! This is likely caused by how the SWIG memory is handled. thermodynamic_state.system is a copy. If you don't assign it to anything, the garbage collector is going to pick it up as soon as it has the chance. The force returned by thermodynamic_state.system.getForce(i) does not own the SWIG C++ memory, so the System gets deallocated, and by the time you call getNumTorsions(), that area of the memory could be already overwritten. I'm surprised you didn't get a segfault.

[jchodera]

This behavior is super confusing.

[jchodera]

I'll think about whether there's a way to avoid breaking obvious idioms from working and forcing the user into highly non-obvious idioms. Otherwise, we lose the utility of using Python in the first place.

[andrrizzi]

The only solution I see is to make ThermodynamicState a complete wrapper of System (meaning that we expose all its methods) and hide the system altogether.

Maybe removing the property and leaving only get_system() can help? You'd still incur into the same problem if you do thermodynamic_state.get_system().getForce(2).getNumTorsions() though.

Unfortunately, I don't see a way around the copy without going back to huge memory consumption and very slow deepcopy() of the states.

[jchodera]

Instead of creating a new object with a new API, I think we should be thinking about a pure-Python System container. I did something like this maybe seven years ago, and I still think it's a good idea: https://github.com/choderalab/pyopenmm/blob/master/pyopenmm/system.py

[kyleabeauchamp]

If we're in the business of a new python object hierarchy, I highly recommend http://www.attrs.org/ as a helper.

[jchodera]

Thanks!

Also, their sphinx docs stylesheet is gorgeous!