MDAnalysis / mdanalysis

MDAnalysis is a Python library to analyze molecular dynamics simulations.
https://mdanalysis.org
Other
1.32k stars 652 forks source link

Neutron reflectivity profile #904

Open khuston opened 8 years ago

khuston commented 8 years ago

Neutron reflectivity profile

I would like to add an analysis that computes neutron reflectivity profiles from simulations of surfaces and interfaces.

However, the isotope of each atom (or isotopic composition of each united-atom) could not be learned from the topology file. If we knew the topology were all-atom, we could make a good guess at the element of each atom, and then assume natural isotopic composition. Even then though, the user would need to specify deuteration.

So here's what I'm thinking:

analysis.reflectivity.guess_isotopes(u)
u.select_atoms(deuterated_united_carbons).isotope = 'CD2'
u.select_atoms(more_deuterated_united_carbons).isotope = 'CD'
u.select_atoms(hydrogen_that_should_be_deuterium).isotope = 'D'
q = np.arange(0.0001, 0.1, 0.0001)
refl = analysis.reflectivity.ReflectivityProfile(q)
refl.run(u)
# result is stored in refl.timeseries

I have pieces I can put together to do the actual processing, but smartly guessing the element(s)/isotope(s) could be annoying. The path of least resistance is to just use the existing element guesser and warn the user that united-atoms will be misidentified as a single nucleus, and that the isotopes should be verified anyway.

Any thoughts? Does this fit with the bauhaus style?

P.S. Edit: The calculation would be performed with the Abeles matrix method.

richardjgowers commented 8 years ago

So does neutron reflectivity change over time, or is it static for a given atom? Bauhaus is mainly about managing how multiframe analysis works, so would only really apply to changing values over time (a trajectory of reflectivities).

Otherwise, adding extra attributes isn't a problem now, so yes functions that looks at a Universe and adds new atom properties based on guesses/heuristics/calculations would be a nice addition.

khuston commented 8 years ago

Ah ok, I illiterately used the timeseries variable above to store the output, but the output is static. I should just describe the calculation:

  1. Group atoms of universe by isotope.
  2. Bin each isotope into number density profiles.
  3. Multiply each isotope's number density profile by the isotope's neutron scattering length.
  4. Sum these together to obtain the scattering length density (SLD) profile.
  5. Compute reflectivity from the SLD profile using the Abeles matrix method.

The analysis should also have start, stop, step, and bins/binsize keywords like LinearDensity.

My only other thought about analysis.reflectivity.guess_isotopes now is that its behavior is different from the similarly-named topology.core.guess_atom_* functions (taking an atomname and returning a mass/element/charge), but I won't worry about that unless there's an objection.

richardjgowers commented 8 years ago

Aha, so as the trajectory moves & density changes the reflectivities will change? This might be cool as an Auxiliary, which @fiona-naughton has been working on recently in #868. This would add the calculation as a an automatic step each time the trajectory moves to a different frame.