[yt-users] Import of crystallographic charge density data

Wed Mar 29 15:32:08 PDT 2017

Hi Nathan,

> Are you specifying a length_unit in your call to load_uniform_grid?

Yes, I’m passing length_unit=‘Angstrom’ to load_uniform_grid and axes_unit=‘Angstrom’ to ProjectionPlot, and specifying my density as ‘me/Angstrom**3'

> Any mapping of f(a, b, c) => (x, y, z) can be used; for instance, we have pixelizers that do non-cartesian coordinates, and so a basis vector transformation wouldn't be so bad.

This is good to know, the transformation is very straight-forward too.

> If you have one you don't mind sharing, it could go on hub.yt <http://hub.yt/> and we can take a look there.  Thanks!

Happy to share it on hub.yt, but I’m not familiar with Girder so I’m not sure the best way to do this. For an example of a raw file, see:

https://github.com/materialsproject/pymatgen/blob/master/test_files/CHGCAR.spin

But that would require a custom parser (file format is documented here: https://cms.mpi.univie.ac.at/vasp/vasp/CHGCAR_file.html ).

However, using the ‘pymatgen’ package, it’s straight forward to get the data as a numpy array:

	from pymatgen.io.vasp.outputs import Chgcar
	c = Chgcar.from_file(‘CHGCAR.spin’)
	charge_density = c.data[‘total’] # returns np.array
	magnetization_density = c.data[‘diff’] # this is scalar magnetization density, just the difference between up and down spin electron density, so has same units as charge_density, and also the same dimensions

Given this, what would be the best way of putting some sample data into hub.yt for testing?

In addition to these CHGCAR files, there are similar files with similar boundary conditions that contain data on e.g. Fermi surfaces, which would be a bit trickier to deal with but also potentially rewarding if we can get good visualizations out.

All the best,

Matthew

> On Mar 28, 2017, at 12:53 PM, Matthew Turk <matthewturk at gmail.com> wrote:
> 
> Hi Matthew,
> 
> This is super cool.
> 
> On Tue, Mar 28, 2017 at 2:23 PM, Matthew Horton <mkhorton at lbl.gov <mailto:mkhorton at lbl.gov>> wrote:
> Hi Nathan,
> 
> Thanks very much for this response!
> 
> > it would help if you can give a little bit more detail about what sort of boundary conditions you need.
> 
> Sure, basically, the simulation box is exactly the same as the orthogonal case currently implemented by `yt.load_uniform_gird(data, … periodicity=(True,True,True)` *except* that in some cases the box has been deformed/sheared… in technical terms, the simulation box is defined by three basis vectors a, b and c which define its axes, and e.g. for a cubic crystal with lattice constant 4Å you’d have a=(4,0,0), b=(0,4,0), c=(0,0,4) and they'd correspond to the x, y and z axes… but in the general case, a b and c could be anything.
> 
> Ah!  This is interesting, and kind of like (I think) the way NIFTI files store their transformations.  I spent some time implementing this at one point; basically, there are a few places that this comes up, both on the visualization and the analysis sides.
> 
> One the analysis side, which is likely easier to fix, there are fields like "x" "y" "z" and the corresponding "cell_volume" fields.  These would just need to have their volume elements computed correctly, which I think would be straightforward.
> 
> On the visualization side, there are two things that would need addressing.  The first of these, which is probably much more straightforward than the other, is coming up with a pixelization routine.  I made some headway into this at one point, where I tried to plot in "real" space versus "reference" space some fMRI datasets, but I ended up having to take other projects on and didn't complete it.  The pixelization routines are not terribly difficult to write; given boundaries, some reference information and the datapoints, it just has to produce a buffer.  Any mapping of f(a, b, c) => (x, y, z) can be used; for instance, we have pixelizers that do non-cartesian coordinates, and so a basis vector transformation wouldn't be so bad.
> 
> The other issue is coming up with a plotting API that allows for both reference and real space plotting, but that can be longer term.
> 
> I'd be interested to get more information from you about this, and we can see if the necessary work can be outlined and maybe taken on.
>  
> 
> Internally, the data is always stored as a uniform grid/array, but the points of the grid don’t correspond to (x,y,z) co-ords, but to fractional (a,b,c) co-ords.
> 
> In practical terms, currently trying to visualize this in yt would lead to a ‘correct’ but distorted/deformed visualization.
> 
> > us = UnitSystem('custom', 'angstrom', 'me', 's', current_mks_unit='1.6e-19*A')
> 
> This looks like it should work, with my projection plots having axes labelled Angstroms, but units still seem to be cm (e.g. it’s showing axes in 10^9 Angstroms)… this may be user error on my part, I think I’m passing the correct length unit arguments etc., but I’ll continue to experiment. 
> 
> I can share a datafile if that’d be helpful at all?
> 
> If you have one you don't mind sharing, it could go on hub.yt <http://hub.yt/> and we can take a look there.  Thanks!
>  
> 
> I discovered yt while looking for a way to automate volumetric rendering of these CHGCAR files…the typical way to visualize these is either slice or isosurface, which loses a lot of information… but it looks like yt would be really useful for doing some more rigorous analysis of these datasets too, so it’d be great if we could get it to work :-)
> 
> Best regards,
> 
> Matthew
> 
>> On Mar 27, 2017, at 4:41 PM, Nathan Goldbaum <nathan12343 at gmail.com <mailto:nathan12343 at gmail.com>> wrote:
>> 
>> Hi Matthew,
>> 
>> First this is the first I've heard of anyone using yt to look at crystallagraphic charge density data, so I'm glad you were able to get as far as you were.
>> 
>> Right now yt is mostly used by astrophysicists and there are still plenty of places where we make assumptions that are only really valid for astrophysics simulation data. We have long-term goals to improve this situation. Getting input on places where yt makes poor assumptions for data from new domains is very valuable for informing our long-term plans for yt development, so thank you for writing in :)
>> 
>> On Mon, Mar 27, 2017 at 5:50 PM, Matthew Horton <mkhorton at lbl.gov <mailto:mkhorton at lbl.gov>> wrote:
>> Hello all,
>> 
>> I’m a new user to yt looking for some general guidance.
>> 
>> 1. I have files describing the charge density in a crystal (CHGCAR files for any DFT/VASP users on the mailing list). I can import these with yt as generic (periodic) unigrid data, however the axes are not always orthogonal… In cubic crystals, `yt.load_uniform_grid` works perfectly, but in some crystals this would give a distorted visualization. What would be the best way to import this kind of data?
>> 
>> Right now yt only supports periodic and isolated boundary conditions. Adding support for additional boundary conditions is a longstanding issue, but most of the time we get requests to add e.g. reflecting, inflow, or diode boundary conditions, not what you're asking about here.
>> 
>> Unfortunately I'm not terribly familiar with crystallography (the last I thought about this was in a solid state physics course I took in undergrad) so it would help if you can give a little bit more detail about what sort of boundary conditions you need.
>> 
>> Unfortunately it's not terribly easy right now to add new boundary conditions, so it might take some work to get yt fully working for non-cubic crystal structures. I don't think anyone is working on making it easier to add new boundary conditions but it's definitely on our long-term roadmap. If you'd like to work on this I'd be happy to chat further about what needs to be done in terms of code changes to get this working.
>> 
>> 2. Regarding units, the units in my file are nominally fractional electrons per grid cell volume, e/Angstrom^3 … importing into yt as a density field, what units should I use? Would be good to keep the data in its natural units, rather than converting to SI/cgs, because they’d be very very small floats otherwise.
>> 
>> I tried 'unit_system = yt.UnitSystem('custom', 'Angstrom', 'g', 's’)' and this works as expected for the length units in my plots, but how do I set the custom Charge unit? What ‘units’ string would I then use for the density field? It doesn’t seem to want to accept ‘Angstrom’, ‘C’ etc.
>> 
>> I *think* this will work:
>> 
>> us = UnitSystem('custom', 'angstrom', 'me', 's', current_mks_unit='1.6e-19*A')
>> 
>> Here 'me' is the electron mass.
>>  
>> 
>> Many thanks for any help!
>> 
>> Best regards,
>> 
>> Matthew
>> _______________________________________________
>> yt-users mailing list
>> yt-users at lists.spacepope.org <mailto:yt-users at lists.spacepope.org>
>> http://lists.spacepope.org/listinfo.cgi/yt-users-spacepope.org <http://lists.spacepope.org/listinfo.cgi/yt-users-spacepope.org>
>> 
>> _______________________________________________
>> yt-users mailing list
>> yt-users at lists.spacepope.org <mailto:yt-users at lists.spacepope.org>
>> http://lists.spacepope.org/listinfo.cgi/yt-users-spacepope.org <http://lists.spacepope.org/listinfo.cgi/yt-users-spacepope.org>
> 
> 
> _______________________________________________
> yt-users mailing list
> yt-users at lists.spacepope.org <mailto:yt-users at lists.spacepope.org>
> http://lists.spacepope.org/listinfo.cgi/yt-users-spacepope.org <http://lists.spacepope.org/listinfo.cgi/yt-users-spacepope.org>
> 
> 
> _______________________________________________
> yt-users mailing list
> yt-users at lists.spacepope.org <mailto:yt-users at lists.spacepope.org>
> http://lists.spacepope.org/listinfo.cgi/yt-users-spacepope.org <http://lists.spacepope.org/listinfo.cgi/yt-users-spacepope.org>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.spacepope.org/pipermail/yt-users-spacepope.org/attachments/20170329/f1e9fd88/attachment.html>