rsdefever / ele

A lightweight package with the periodic table of the elements
https://ele-ment.readthedocs.io
MIT License
12 stars 5 forks source link

Isotopes #43

Open mattwthompson opened 3 years ago

mattwthompson commented 3 years ago

Thank you for writing this nice little package. It's satisfying to be able to solve a small problem without needing to pull in a big library.

One usage question I have is how to deal with isotopes of elements, i.e. Deuterium in heavy water. I can make a custom element as I please, but this doesn't propogate out to the lookup functions.

$ cat repro.py
import ele

deuterium = ele.element.Element(name="deuterium", symbol="D", mass="2.0", atomic_number="1")
ele.element_from_symbol("D")
$ python repro.py
Traceback (most recent call last):
  File "repro.py", line 4, in <module>
    ele.element_from_symbol("D")
  File "/Users/mwt/miniconda3/envs/openff-system/lib/python3.8/site-packages/ele/element.py", line 74, in element_from_symbol
    raise ElementError(f"No element with symbol {symbol}")
ele.exceptions.ElementError: No element with symbol D

Attempting to force it in (i.e. ele.element.elements.append(deuterium)) doesn't work, and it's probably not a good idea to modify module variables directly in the absence of an API to do so. What would be the suggested usage here, or is this intended to be out of the scope of this project?

rsdefever commented 3 years ago

This is a good question and I'd be interested in your thoughts on multiple fronts here. Here are mine:

In terms of support for different isotopes, I think it is generally outside the scope of the package (I could perhaps be convinced otherwise). Right now the masses for most elements are "conventional atomic weights", which, IIRC, are effectively isotope-weighted by natural abundance. Given this design choice, I don't see a great way to handle isotopes.

Adding an element It is not a use-case that I had considered but I am not per-se opposed to providing an API to do so. It could be a nice feature as it would allow (as in your example) the addition of specific isotopes that might be important for certain projects.

rsdefever commented 3 years ago

@mattwthompson do other isotopes beyond deuterium have specific names? I'm trying to think of a general way to incorporate support for isotopes in a consistent fashion. The things that come to mind. JSON dict keys become (atomic_number, mass_number) (see #41). I was thinking for consistency, it would be nice to have the element.name be something like name="element_name-mass_number", and similar for the element.symbol. But this obviously then would mean that we would not have deuterium and D but rather hydrogen-2 and H-2.

Thoughts?

mattwthompson commented 3 years ago

Deuterium is the only one I can think of that has a unique name, and it's also by far the most relevant one for molecular simulation as far as I know. Even stuff like Carbon-13 probably comes out in the wash compared to the doubled mass of deuterium.

This is made slightly trickier if considering mixing water and heavy water; this is probably not an important use case but some neutron scattering experiments can cancel out the signal from water by using a mixture at some particular concentration (IIRC hydrogen has an anomalously negative coefficient in some case, and deuterium and other heavy atoms have a positive one, so mixing them can make water "invisible" in some experiments). At some point I wanted to run some simulations to even see how different the physical properties of deuterating i.e. TIP3P were. Mostly to satisfy my confirmation bias that the different was negligible.

A specific use case here would be typing a topology with a force field that includes parameters for heavy atoms. This probably means that some "atom name"-like variable is being exposed in an iterator, since lookup by atomic number is invalid and lookup by mass is flaky. Lookup by symbol could work.

for idx, atom in enumerate(some_iterator.atoms):
    if atom.name == "deuterium"  # or something like that
        # special case
    else:
        element = element_from_name(atom.name)  # or whatever lookup

or the lookup functions could be wrapped, with or without adding on new "elements" to the elements list:

def element_from_name_with_isotopes(name):
    if name in KNOWN_ISOTOPES:
        # special case
    else:
        return element_from_name(name)

I think these special cases can be handled downstream. It probably makes sense to keep this package as lightweight as possible until there are cases that cannot be disambiguated by downstream developers. Handling cases by slightly extending the API of this package also means the use case of creating new elements here goes away. I'm pretty sure this code can be copy-pasted and modified to be a list of non-standard elements.

https://github.com/rsdefever/ele/blob/64a7e62464d59c5b12c38b51d626d328b448d979/ele/element.py#L202-L210

There may or may not be more considerations when dealing with UA/CG models, but I'm not working with those at the moment.

rsdefever commented 3 years ago

Based on this:

I think these special cases can be handled downstream. It probably makes sense to keep this package as lightweight as possible until there are cases that cannot be disambiguated by downstream developers. Handling cases by slightly extending the API of this package also means the use case of creating new elements here goes away. I'm pretty sure this code can be copy-pasted and modified to be a list of non-standard elements.

I'm going to close the issue. Feel free to re-open if you need support from this end.

rsdefever commented 3 years ago

I'm reopening this issue. Upon further thought, I'm wondering if we can write some general functionality to support isotope-specific behavior. E.g., write a function that take the atomic number and mass number and creates an Isotope object (subclassesed from an Element object?) on the fly.

ele.isotope_from_symbol("C", mass_number=13)
ele.isotope_from_atomic_number(6, mass_number=13)
...
StanczakDominik commented 3 years ago

Hey, just found your package, saw this issue and wanted to point to @namurphy's implementation of this in plasmapy.particles - could give you some ideas :)