The Kitchin Research Group

We will loop over a set of calculations that define the equation of state for two different Al-Ni alloy structures. For each calculation this code generates a level 3 headline with some tags derived from the calculation directory. I chose level 3 so they would all be nested in this level 2 headline. We write out two sets of data because later we will selectively get portions of the data.

from jasp import *
from jasp.utils import get_jasp_dirs

for d in get_jasp_dirs('10'):
    with jasp(d) as calc:
        print calc.org(level=3)

for d in get_jasp_dirs('11'):
    with jasp(d) as calc:
        print calc.org(level=3)

We can use column view to see the results in a tabular like form. For example, here we capture the VOLUME and TOTAL_ENERGY properties.

That is more or less what we want, but we cannot use that table as a data source very directly. While you can use it as a data source, you just get a string of the whole table, e.g.

print data.__class__

<type 'str'>

You could parse that string, or try converting it to a table, but we skip that for now. In the next section we consider another approach that more directly gets at the data.

Let us try a different approach. Here we map over headlines that match a criteria, and get specific data. We add a header of labels and horizontal line to the beginning of the table. Here we find the results tagged with 11, and that match f_dash* (we actually use a regular expression).

(cons '(label volume energy) 
      (cons 'hline
            (org-map-entries
             (lambda ()
               (list (nth 4 (org-heading-components))
                     (org-entry-get (point) "VOLUME")
                     (org-entry-get (point) "TOTAL_ENERGY"))) 
             ;; this next line specifies headlines tagged with 11
             ;; and that match the regular expression in {}
             "+11+{f_dash_.*}")))

Now, since there are strings in the table (or in the results from the named code-block), we will have to cast the results to floats using a numpy array. Note it appears that we actually run the code block to get the results, they do not come from the table. The table is just helpful to know what data you are using.

import matplotlib.pyplot as plt
import numpy as np

print E

plt.plot(np.array(V), np.array(E), 'ko')
plt.xlabel('Volume $\AA^3$')
plt.ylabel('Energy (eV)')
plt.savefig('11-eos.png')

['-20.038179', '-19.912894', '-19.662755', '-20.275411', '-20.249632', '-20.202731', '-20.133749']

And here you have the figure. We extracted some data into an org-file in human readable form. Then, we extracted the data from those entries to a table, and finally plotted it.

Say you want the stress from some calculation. Each heading is identified by a CUSTOM_ID, so we can temporarily go to that headline, and get the property. Easy!

(save-excursion
  (org-open-link-from-string "[[#10/eos-exp/f+0]]")
  (org-entry-get (point) "STRESS"))

[-0.003191325893994518, -0.003191325893994518, -0.003191325893994518, -0.00042576435655845916, 0.00042576435655845916, 0.00042576435655845916]

Note, however, that the value returned is a string, so we need to convert it to numbers if we want to do anything with it, e.g. find the maximum component. Here, we eval the string in Python to make it a list, and then find the maximum value.

import numpy as np
print np.abs(np.array(eval(stress))).max()

0.00319132589399

Next, we get the unit cell from that calculation, so we can calculate the volume using python. The unit cell is stored in a table and we extract the unit cell vectors from it in the code block header.

import numpy as np

print('cell = {}'.format(np.array(data)))
print('Volume = {}'.format(np.linalg.det(np.array(data))))

cell = [[-1.98408  0.      -1.98408]
 [ 0.      -1.98408  1.98408]
 [-3.53625  3.53625  3.53625]]
Volume = 41.7621235495

Compare that to the volume we computed when the entry was generated, and you can see they are the same.

(save-excursion
  (org-open-link-from-string "[[#10/eos-exp/f+0]]")
  (org-entry-get (point) "VOLUME"))

41.7621235495

The last example is getting the formula for the calculation. The atoms positions are stored in a table in the entry, and the species for each atom is stored in the first column. We extract the species into a variable, then then get counts for each species, and finally construct a string representation of the formula.

(let* ((unique-species '())
      (counts))
  ;; first get unique species
  (mapcar (lambda (s) (add-to-list 'unique-species s)) species)

  ;; now get counts
  (setq counts
        (mapcar (lambda (s)
                  (--count (string= it s) species)) unique-species))

  ;; Now put them together
  (mapconcat 'identity (cl-mapcar (lambda (s n)
                  (format "%s%s" s n)) unique-species counts) ""))

Ni1Al2

Suppose we want to sort the elements alphabetically. We can make a list of cons cells of element and count, sort the list by element, and then put them back together. We do this in lisp.

(let* ((unique-species '())
      (counts))
  ;; first get unique species
  (mapcar (lambda (s) (add-to-list 'unique-species s)) species)

  ;; now get counts
  (setq counts
        (mapcar (lambda (s)
                  (--count (string= it s) species)) unique-species))

  (mapconcat (lambda (x)
               (format "%s%s" (car x) (cdr x)))
             (cl-sort
              (cl-mapcar (lambda (a b) (cons a b)) unique-species counts)
              'string-lessp :key 'car)
             ""))

Al2Ni1

That is not too bad. It is a little verbose, some of that could be put into a library to be reused. For example, here is the same code in Python for the unsorted formula

from collections import Counter
counts = Counter(species)
print ''.join(['{0}{1}'.format(x, counts[x]) for x in counts])

Ni1Al2

and for the sorted formula.

from collections import Counter
counts = Counter(species)
print ''.join(['{0}{1}'.format(x, counts[x]) for x in sorted(counts.keys())])

Al2Ni1

Obviously a lot must be happening behind the scenes in those short, dense python blocks. For example, the list comprehension syntax in Python is functionally similar to the mapcar code, and the join function is similar to the mapconcat code in lisp.

This has some potential. It is a hybrid approach to human and machine readable data, that enables inline analysis. You can mix programming languages, which is great if one language has unique capabilities or is just more convenient. There are some limitations, mostly that org-properties are read as strings, and they have to be converted to the right data type, e.g. numbers or lists or arrays, etc… But, these can be saved as forms that can be directly evaled in some languages. On the other hand, it is possible to do analysis on another machine where the original data does not exist.

This is a very different mental model of data than I am used to. We typically have thought about data in files, and paths to get to that data. Here, data is in headings, and we use links to get to the data. We also use different syntax to link data to code blocks. That means in some cases the code blocks cannot be run independently of this file.

This looks ok as a data sharing mechanism. The json file here is pretty small (6.8 kb), and pretty easy to work with. Clearly some thought must go into how the data is structured so you know what to get, and how you get it. That could even be documented in the json file itself. For instance, each calculator has a doc element that describes what is in it. The json file we made above also has that data.

import json

with open('eos.json', 'rb') as f:
    d = json.loads(f.read())

print d['doc']
print
print d['calculations'][0]['doc']

JSON file containing a set of calculations for the equation of state of Cu

JSON representation of a VASP calculation.

energy is in eV
forces are in eV/\AA
stress is in GPa (sxx, syy, szz,  syz, sxz, sxy)
magnetic moments are in Bohr-magneton
The density of states is reported with E_f at 0 eV.
Volume is reported in \AA^3
Coordinates and cell parameters are reported in \AA

If atom-projected dos are included they are in the form:
{ados:{energy:data, {atom index: {orbital : dos}}}

For archival purposes you may want to put a creation date, contact data, etc… in the file too.