A sexp version of a bibtex entry

Posted June 10, 2015 at 08:54 AM | categories: bibtex, lisp | tags:

Below you see a typical bibtex entry. Today we explore an alternate approach to represent the information (data) in that entry as s-expressions, i.e. as a lisp data structure. Why? because it seems like an interesting exploration! 

@article{hallenbeck-2013-effec-o2,
  author =       "Hallenbeck, Alexander P. and Kitchin, John R.",
  title =        {Effects of \ce{O_2} and \ce{SO_2} on the capture capacity of a
                  primary-amine based polymeric \ce{CO_2} sorbent},
  keywords =     {RUA, orgmode},
  journal =      "Industrial \& Engineering Chemistry Research",
  pages =        "10788-10794",
  year =         2013,
  volume =       {52},
  number =       {31},
  doi =          "10.1021/ie400582a",
  url =          "http://pubs.acs.org/doi/abs/10.1021/ie400582a",
  eprint =       "http://pubs.acs.org/doi/pdf/10.1021/ie400582a",
}

Here is what that same data structure might look like as a sexp-based lisp data structure. 

(article "hallenbeck-2013-effec-o2"
         (author "Hallenbeck, Alexander P. and Kitchin, John R.")
         (title "Effects of \ce{O_2} and \ce{SO_2} on the capture capacity of a primary-amine based polymeric \ce{CO_2} sorbent")
         (journal "Industrial \& Engineering Chemistry Research")
         (pages "10788-10794")
         (year 2013)
         (number 31)
         (doi "10.1021/ie400582a")
         (url "http://pubs.acs.org/doi/abs/10.1021/ie400582a")
         (eprint "http://pubs.acs.org/doi/pdf/10.1021/ie400582a"))

We can retrieve data from the sexp form pretty easily. Here we get the authors. 

(let* ((art '(article "hallenbeck-2013-effec-o2"
                      (author "Hallenbeck, Alexander P. and Kitchin, John R.")
                      (title "Effects of \ce{O_2} and \ce{SO_2} on the capture capacity of a primary-amine based polymeric \ce{CO_2} sorbent")
                      (journal "Industrial \& Engineering Chemistry Research")
                      (pages "10788-10794")
                      (year 2013)
                      (number 31)
                      (doi "10.1021/ie400582a")
                      (url "http://pubs.acs.org/doi/abs/10.1021/ie400582a")
                      (eprint "http://pubs.acs.org/doi/pdf/10.1021/ie400582a")))
       (fields (cddr art)))
  (cadr (assoc 'author fields)))

Hallenbeck, Alexander P. and Kitchin, John R.

That is simple enough you might just write a little function to streamline it like this, and return a formatted string. 

(defun get-article-field (article field)
  "Return value of FIELD in ARTICLE."
  (cadr (assoc field (cddr article))))

(let ((art '(article "hallenbeck-2013-effec-o2"
                     (author "Hallenbeck, Alexander P. and Kitchin, John R.")
                     (title "Effects of \ce{O_2} and \ce{SO_2} on the capture capacity of a primary-amine based polymeric \ce{CO_2} sorbent")
                     (journal "Industrial \& Engineering Chemistry Research")
                     (pages "10788-10794")
                     (year 2013)
                     (number 31)
                     (doi "10.1021/ie400582a")
                     (url "http://pubs.acs.org/doi/abs/10.1021/ie400582a")
                     (eprint "http://pubs.acs.org/doi/pdf/10.1021/ie400582a"))))
  (format "%s, doi:%s (%s)"
          (get-article-field art 'author)
          (get-article-field art 'doi)
          (get-article-field art 'year)))

Hallenbeck, Alexander P. and Kitchin, John R., doi:10.1021/ie400582a (2013)

You might be wondering, why is that even a little bit interesting? One reason is that it looks a little like what lisp returns after parsing an xml file. Another is, the data structure looks kind of like data, but it is also some code, if article was defined as a function! Let us consider what this might look like. I use a macro to define the field functions since in this case they all do the same thing, and these simply return a string with the field-name and value in curly brackets. We eval the macro to make sure it defines the function. I define an article function that wraps the fields in @bibtex-key{fields}, which defines a bibtex entry. 

(defmacro make-field (field-name)
  "define a field that returns a string"
  `(defun ,(intern field-name) (content)
     (format "  %s = {%s}" ,field-name content)))

(loop for field in '("author" "title" "journal" "pages" "number" "doi" "url" "eprint" "year")
  do (eval `(make-field ,field)))

(defun article (bibtex-key &rest fields)
  (concat
   (format "@article{%s,\n" bibtex-key)
   (mapconcat (lambda (field) (eval field)) fields ",\n")
   "\n}\n"))

(article "hallenbeck-2013-effec-o2"
         (author "Hallenbeck, Alexander P. and Kitchin, John R.")
         (title "Effects of \ce{O_2} and \ce{SO_2} on the capture capacity of a primary-amine based polymeric \ce{CO_2} sorbent")
         (journal "Industrial \& Engineering Chemistry Research")
         (pages "10788-10794")
         (number 31)
         (year 2013)
         (doi "10.1021/ie400582a")
         (url "http://pubs.acs.org/doi/abs/10.1021/ie400582a")
         (eprint "http://pubs.acs.org/doi/pdf/10.1021/ie400582a"))

@article{hallenbeck-2013-effec-o2,
  author = {Hallenbeck, Alexander P. and Kitchin, John R.},
  title = {Effects of ce{O_2} and ce{SO_2} on the capture capacity of a primary-amine based polymeric ce{CO_2} sorbent},
  journal = {Industrial & Engineering Chemistry Research},
  pages = {10788-10794},
  number = {31},
  year = {2013},
  doi = {10.1021/ie400582a},
  url = {http://pubs.acs.org/doi/abs/10.1021/ie400582a},
  eprint = {http://pubs.acs.org/doi/pdf/10.1021/ie400582a}
}

Wow. We executed our data structure, and got a bibtex entry! That seems moderately interesting to me. Next is an example of taking the same data structure and rendering it as xml. This is some lispy wizardry, rather than use a macro to define functions, I temporarily define functions within a cl-flet macro, which I have to collect as a list of code. Then, I eval the list. This feels pretty odd, but seems like a lispy kind of thing to do. 

(eval
 (list 'cl-flet
       (append (loop for field in '("author" "title" "journal" "pages"
                                      "number" "doi" "url" "eprint" "year")
                       collect (list (intern field)
                                     '(content)
                                     `(format "  <%s>%s</%s>" ,field content ,field)))
               '((article (bibtex-key &rest fields)
                          (concat
                           (format
                            "<article bibtex-key=\"%s\">\n" bibtex-key)
                           (mapconcat (lambda (field) (eval field)) fields "\n")
                           "\n</article>")))
               )
       ;; body of cl-flet
       '(article "hallenbeck-2013-effec-o2"
                (author "Hallenbeck, Alexander P. and Kitchin, John R.")
                (title "Effects of \ce{O_2} and \ce{SO_2} on the capture capacity of a primary-amine based polymeric \ce{CO_2} sorbent")
                (journal "Industrial \& Engineering Chemistry Research")
                (pages "10788-10794")
                (number 31)
                (year 2013)
                (doi "10.1021/ie400582a")
                (url "http://pubs.acs.org/doi/abs/10.1021/ie400582a")
                (eprint "http://pubs.acs.org/doi/pdf/10.1021/ie400582a"))))

<article bibtex-key="hallenbeck-2013-effec-o2">
  <author>Hallenbeck, Alexander P. and Kitchin, John R.</author>
  <title>Effects of ce{O_2} and ce{SO_2} on the capture capacity of a primary-amine based polymeric ce{CO_2} sorbent</title>
  <journal>Industrial & Engineering Chemistry Research</journal>
  <pages>10788-10794</pages>
  <number>31</number>
  <year>2013</year>
  <doi>10.1021/ie400582a</doi>
  <url>http://pubs.acs.org/doi/abs/10.1021/ie400582a</url>
  <eprint>http://pubs.acs.org/doi/pdf/10.1021/ie400582a</eprint>
</article>

Prefer json? No problem, just reformat the functions! 

(eval
 (list 'cl-flet
       (append (loop for field in '("author" "title" "journal" "pages"
                                      "number" "doi" "url" "eprint" "year")
                       collect (list (intern field)
                                     '(content)
                                     `(format "   \"%s\": \"%s\"" ,field content)))
               '((article (bibtex-key &rest fields)
                          (concat
                           (format
                            "{\"article\":\n  {\"bibtex-key\": \"%s\",\n" bibtex-key)
                           (mapconcat (lambda (field) (eval field)) fields ",\n")
                           "}\n}"))))
       ;; body of cl-flet
       '(article "hallenbeck-2013-effec-o2"
                (author "Hallenbeck, Alexander P. and Kitchin, John R.")
                (title "Effects of \ce{O_2} and \ce{SO_2} on the capture capacity of a primary-amine based polymeric \ce{CO_2} sorbent")
                (journal "Industrial \& Engineering Chemistry Research")
                (pages "10788-10794")
                (number 31)
                (year 2013)
                (doi "10.1021/ie400582a")
                (url "http://pubs.acs.org/doi/abs/10.1021/ie400582a")
                (eprint "http://pubs.acs.org/doi/pdf/10.1021/ie400582a"))))

{"article":
  {"bibtex-key": "hallenbeck-2013-effec-o2",
   "author": "Hallenbeck, Alexander P. and Kitchin, John R.",
   "title": "Effects of ce{O_2} and ce{SO_2} on the capture capacity of a primary-amine based polymeric ce{CO_2} sorbent",
   "journal": "Industrial & Engineering Chemistry Research",
   "pages": "10788-10794",
   "number": "31",
   "year": "2013",
   "doi": "10.1021/ie400582a",
   "url": "http://pubs.acs.org/doi/abs/10.1021/ie400582a",
   "eprint": "http://pubs.acs.org/doi/pdf/10.1021/ie400582a"}
}

Is this useful? Great question. I don't plan to convert by bibtex files to sexp format anytime soon ;) The format I used above is just a simple one. It might be desirable to include individual authors instead of an author string, and maybe support attributes to establish an author order. An author structure might be more complex to include scientific ids like an orcid, alternative names, etc… Finally, the s-exp data structure is super easy to use in lisp, but other languages would have parse it into some native structure the way they parse json or xml. There is limited support for s-expressions in most other non-lispy languages.

I like the idea of data representation as code, and its conversion to some other kind of format. It is subtle here, but notice we never had to write a parser for the sexp notation. That already exists as the lisp interpreter. We did write code to use the data, and convert the data. The sexp notation is pretty easy to write, in contrast to the xml or json representations. Some interesting issues might be what to do with fields that are not defined, perhaps a macro would be used on the fly, or in the cl-flet definition. It is hard to imagine doing these things in another language than lisp!

Copyright (C) 2015 by John Kitchin. See the License for information about copying.

org-mode source

Org-mode version = 8.2.10

Discuss on Twitter

Python data structures to lisp

Posted May 16, 2015 at 10:47 AM | categories: lisp, emacs, python | tags:

I have an idea in mind that would use the output of python scripts in lisp functions. Xah Lee posted an idea for writing emacs commands in scripting languages . In this post I want to explore an extension of the idea, where a Python script will return output that can be read in Lisp, e.g. we can convert a Python list to a lisp list, or a dictionary to an a-list or p-list. I can already see that simple data structures will be "simple", and arbitrary data structures will offer a lot of challenges, e.g. nested lists or dictionaries…

If I could add some custom functions to the basic builtin types in Python, then I could use another approach to format python objects as lisp data types. This isn't recommended by Pythonistas, but I guess they don't want to use lisp as much as I do ;) I found this approach to modifying builtins:

http://stackoverflow.com/questions/2444680/how-do-i-add-my-own-custom-attributes-to-existing-built-in-python-types-like-a

We use that almost verbatim here to get what I want. This is a super low level way to add functions to the builtins. I add some simple formatting to floats, ints and strings. I add a more complex recursive formatting function to lists, tuples and dictionaries. A dictionary can be represented as an alist or plist. Both examples are shown, but I leave the alist version commented out. Finally, we add a lispify function to numpy arrays. 

import ctypes as c

class PyObject_HEAD(c.Structure):
    _fields_ = [('HEAD', c.c_ubyte * (object.__basicsize__ -
                                      c.sizeof(c.c_void_p))),
                ('ob_type', c.c_void_p)]

_get_dict = c.pythonapi._PyObject_GetDictPtr
_get_dict.restype = c.POINTER(c.py_object)
_get_dict.argtypes = [c.py_object]

def get_dict(object):
    return _get_dict(object).contents.value

get_dict(str)['lisp'] = lambda s:'"{}"'.format(str(s))
get_dict(float)['lisp'] = lambda f:'{}'.format(str(f))
get_dict(int)['lisp'] = lambda f:'{}'.format(str(f))

import collections
import numpy as np

def lispify(L):
    "Convert a Python object L to a lisp representation."
    if (isinstance(L, str)
        or isinstance(L, float)
        or isinstance(L, int)):
        return L.lisp()
    elif (isinstance(L, list)
          or isinstance(L, tuple)
          or isinstance(L, np.ndarray)):
        s = []
        for element in L:
            s += [element.lisp()]
        return '(' + ' '.join(s) + ')'
    elif isinstance(L, dict):
        s = []
        for key in L:
            # alist format
            # s += ["({0} . {1})".format(key, L[key].lisp())]
            # plist
            s += [":{0} {1}".format(key, L[key].lisp())]
        return '(' + ' '.join(s) + ')'

get_dict(list)['lisp'] = lispify
get_dict(tuple)['lisp'] = lispify
get_dict(dict)['lisp'] = lispify
get_dict(np.ndarray)['lisp'] = lispify

Let us test these out. 

from pylisp import *
a = 4.5
print int(a).lisp()
print a.lisp()
print "test".lisp()

print [1, 2, 3].lisp()
print (1, 2, 3).lisp()

print [[1, 3], (5, 6)].lisp()

print {"a": 5}.lisp()
print [[1, 3], (5, 6), {"a": 5, "b": "test"}].lisp()


A = np.array([1, 3, 4])
print A.lisp()
print ({"tree": [5, 6]}, ["a", 4, "list"], 5, 2.0 / 3.0).lisp()

4
4.5
"test"
(1 2 3)
(1 2 3)
((1 3) (5 6))
(:a 5)
((1 3) (5 6) (:a 5 :b "test"))
(1 3 4)
((:tree (5 6)) ("a" 4 "list") 5 0.666666666667)

Now, is that better than a single lisp function with a lot of conditionals to handle each type? I am not sure. This seems to work pretty well.

Here is how I imagine using this idea. We would have some emacs-lisp variables and use them to dynamically generate a python script. We run the python script, capturing the output, and read it back in as a lisp data structure. Here is a simple kind of example that generates a dictionary. 

(let* ((elisp-var 6)
       (result)
      (script (format "
from pylisp import *
print {x: [2*y for y in range(x)] for x in range(1, %s)}.lisp()
" elisp-var)))

  ;; start a python process
  (run-python)
  (setq result (read (python-shell-send-string-no-output
   script)))
  (plist-get result :5))

(0 2 4 6 8)

That seems to work pretty well. One alternative idea to this is Pymacs , which I have written about before . This project isn't currently under active development, and I ran into some difficulties with it before.

Here we can solve the problem I previously posed and get the result back as an elisp float, and then reuse the result 

(let* ((myvar 3)
       (script (format "from pylisp import *
from scipy.optimize import fsolve
def objective(x):
    return x - 5

ans, = fsolve(objective, %s)
print ans.lisp()" myvar)))
  (run-python)
  (setq result (read (python-shell-send-string-no-output
                       script)))
  (- 5 result))

0.0

Bottom line: we can write python code in lisp functions that are dynamically updated, execute them, and get lisp data structures back for simple data types. I think that could be useful in some applications, where it is easier to do parsing/analysis in Python, but you want to do something else that is easier in Lisp.

Copyright (C) 2015 by John Kitchin. See the License for information about copying.

org-mode source

Org-mode version = 8.2.10

Discuss on Twitter

A prototype implementation of jasp in emacs-lisp

Posted December 24, 2014 at 11:41 AM | categories: lisp, ase, vasp, emacs | tags:

I want to try implementing an interface to VASP in emacs-lisp. VASP is a program that uses density functional theory to calculate properties of materials. VASP was designed for the user to create several text-based input files that contain an atomic geometry, and calculation parameters. Then you run the VASP program, which reads those input files and generates output files. Finally, you parse out the results you want from the output files.

It is moderately tedious to do that, so we already have a very extensive Python interface (http://github.com/jkitchin/jasp ) that automates the input file creation and output file parsing, and with org-mode we have a pretty good literate research environment to document what we do. But, the Python/emacs environment is not as integrated as the emacs-lisp/emacs environment is, particularly when it comes to accessing documentation. So, I want to try out a lisp approach to see what this would look like, and if it has any benefits.

The bare bones implementation will have an Atom object, to store the type of atom and its coordinates, an Atoms object which will be a collection of atoms, and a unit cell, and a Calculator object that will store calculation parameters.

Then, we will try using it to see what advantages it might have. This will be a moderately long post.

1 The Atom class

This is my first serious effort to use the object system in emacs-lisp, so I do not claim it is optimal, or beautiful. It is functional though. We make a simple Atom class that holds the chemical symbol, and xyz coordinates. We do not consider having a magnetic moment at this point, and this class has no methods. 

(defclass Atom ()
  ((symbol :initarg :symbol
           :documentation "A string for the chemical symbol")
   (x :initarg :x
      :documentation "The x coordinate")
   (y :initarg :y
      :documentation "The y coordinate")
   (z :initarg :z
      :documentation "The z coordinate"))
 "A class to represent an atom.")

(provide 'Atom)

Atom

Let us try it out. We make an Atom, then get the symbol using the oref function. I don't think we need a "Name" for the atom, so the second argument is nil. 

(Atom nil :symbol "C" :x 0 :y 0 :z 0)

[object Atom nil "C" 0 0 0]

(let ((a1 (Atom nil :symbol "C" :x 0 :y 0  :z 0)))
  (oref a1 symbol))

C

It is not difficult to modify an Atom object. We use oset. 

(let ((a1 (Atom nil :symbol "C" :x 0 :y 0  :z 0)))
  (oset a1 x 2.5)
  a1)

[object Atom nil "C" 2.5 0 0]

2 The Atoms class

We need at Atoms object, which will store a list of Atom objects, and a unit cell. I do not know how to make this class act like a list the way that the Atoms object in ase acts. We will have to access the list of atoms as a slot. We define one method here to get the ith atom from the object. 

(defclass Atoms ()
  ((atoms :initarg :atoms
          :documentation "A list of `Atom' objects")
   (cell :initarg :cell
         :documentation "A 3-tuple of lengths of an orthorhombic unit cell"))
  "A class to represent an atoms object.")


(defmethod get-atom ((atoms Atoms) &rest args)
  "Get a list of atoms in ARGS.
Return atom if ARGS contains one element, a list of atom objects
otherwise."
  (cond
   ((= 1 (length args))
    (elt (oref atoms atoms) (car args)))
   (t
    (mapcar (lambda (i)
              (elt (oref atoms atoms) i))
            args))))

(provide 'Atoms)

get-atom

That seems to do what we need. It has some limitations:

  1. We only allowed for orthorhombic unit cells
  2. We did not enable any constraints or magnetic moments on the atoms.

Here is an example Atoms object. Note we use `, notation to ensure each Atom is created, since Atom is a constructor function that returns the object. 

(Atoms nil
       :atoms `(,(Atom nil :symbol "C" :x 0 :y 0 :z 0)
                ,(Atom nil :symbol "O" :x 1.1 :y 0 :z 0))
       :cell '(8 9 10))

[object Atoms nil ([object Atom nil "C" 0 0 0] [object Atom nil "O" 1.1 0 0]) (8 9 10)]

We can drill into the object, e.g. to get the second atom: 

(let ((A1 (Atoms nil
                 :atoms `(,(Atom nil :symbol "C" :x 0 :y 0 :z 0)
                          ,(Atom nil :symbol "O" :x 1.1 :y 0 :z 0))
                 :cell '(8 9 10))))
  (get-atom A1 1))

[object Atom nil "O" 1.1 0 0]

We can modify the atoms in place like this. Suppose we want to change the symbol of the first atom to "O". We use setf for this too. 

(let ((A1 (Atoms nil
                 :atoms `(,(Atom nil :symbol "C" :x 0 :y 0 :z 0)
                          ,(Atom nil :symbol "O" :x 1.1 :y 0 :z 0))
                 :cell '(8 9 10))))
  (oset (get-atom A1 0) symbol "O")
  A1)

[object Atoms nil ([object Atom nil "O" 0 0 0] [object Atom nil "O" 1.1 0 0]) (8 9 10)]

The only think I do not like about this syntax is the need to get the list of atoms from the object. That is a little clunkier than the Python analog where the object is a list itself. That may be just my inexperience with emacs-lisp. Probably you can define some getter function to smooth this over.

This Atoms class lacks much of the functionality of the ase.Atoms class, but it is sufficient for this prototype.

3 The Calculator class

Next, we need our Calculator. This will store the parameters, and be responsible for creating the INCAR, POSCAR, KPOINTS, and POTCAR files, running a calculation, and getting data from the output. We also create a with-current-directory macro that will temporarily change the working directory since VASP uses the same filenames over and over, but in different directories. 

(defmacro with-current-directory (directory &rest body)
  "Set the working directory temporarily set to DIRECTORY and run BODY.
DIRECTORY is expanded, and create it and its parents if needed."
  `(progn
     (unless (file-exists-p (file-name-as-directory
                             (expand-file-name ,directory)))
       (make-directory ,directory t))
     
     (let ((default-directory (file-name-as-directory
                                (expand-file-name ,directory)))) 
        ,@body)))


(defclass Jasp ()
  ((wd :initarg :wd
       :initform "."  ; default to the current working directory
       :documentation "Directory to run calculation in.")
   (encut :initarg :encut
          :documentation "Positive number in eV for planewave cutoff.
See URL `http://cms.mpi.univie.ac.at/vasp/vasp/ENCUT_tag.html'.")
   (nbands :initarg :nbands
           :documentation "Integer number of bands.
See URL `http://cms.mpi.univie.ac.at/vasp/vasp/NBANDS_tag.html'.")
   (kpts :initarg :kpts
         :initform (1 1 1)  ; default value
         :documentation "3-tuple for Monkhorst-Pack K-point grid.")
   (xc :initarg :xc
       :documentation "String of exchange correlation functional.")
   (atoms :initarg :atoms
          :documentation "An `Atoms' object."))
 "A class to represent a calculator that runs VASP.")


(defmethod view-atoms ((calc Jasp))
  "Open the ase-gui"
  (unless (and (file-exists-p "POSCAR")
               (file-exists-p "POTCAR"))
    (write-poscar calc)
    (write-potcar calc))
  (shell-command "ase-gui POSCAR"))


(defmethod write-poscar ((calc Jasp))
  "create a POSCAR file for CALC."
  (with-temp-file "POSCAR"
    (insert "Created by jasp.el\n")
    (insert "  1.0") ; unit cell scale factor
    (let* ((atoms (oref calc atoms))
           (cell (oref atoms cell)))
      (loop for v in cell
            for i below (length cell)     
            do
            (insert "\n")
            (loop for j below (length cell)
                  do
                  (if (equal i j)
                      (insert (format " %f " (float (elt cell i))))
                    (insert (format " %f " 0.0 ))))))
    ;; The next line is counts for each atom type. For each number in
    ;; this line, there will be a copy of the POTCAR in the POTCAR
    ;; file. In ase, we sort the atoms and group them so that there is
    ;; only one POTCAR per atom. We do not do that here yet. We will
    ;; have a POTCAR for each atom.
    (insert "\n")
    (loop for atom in (oref (oref calc atoms) atoms)
          do (insert (format "1 ")))
    
    ;; now we do the atoms
    (insert "\nCartesian\n")
    (loop for atom in (oref (oref calc atoms) atoms)
          do
          (insert
           (format "%f %f %f\n"
                   (oref atom x)
                   (oref atom y)
                   (oref atom z))))
    (buffer-string)))


(defmethod write-kpoints ((calc Jasp))
  "Create KPOINTS file for CALC. 
Limited to automatic generation, and no offset."
  (with-temp-file "KPOINTS"
    (insert "Automatic mesh
0
Monkhorst-Pack
")
    (dolist (k (oref calc kpts))
      (insert (format "%4d " k)))
    (insert "\n0.0 0.0 0.0")
    (buffer-string)))


(defmethod write-potcar ((calc Jasp))
  "Generate the POTCAR file for CALC.
No `Atom' grouping is done, there is one POTCAR for each atom."
  (with-temp-file "POTCAR"
    (let ((xc (oref calc xc))
          (atoms (oref calc atoms))
          (vasp_pp_path (getenv "VASP_PP_PATH")))
      (loop for atom in (oref atoms atoms)
            do
            (insert-file-contents
             (f-join
              vasp_pp_path
              (concat "potpaw_" xc)
              (oref atom symbol)
              "POTCAR"))))
    (buffer-string)))


(defmethod write-incar ((calc Jasp))
  "Generate the INCAR file for CALC."
  (with-temp-file "INCAR"
    (insert (format "ENCUT = %f\n" (oref calc encut)))
    (insert (format "NBANDS = %d\n" (oref calc nbands)))
    (buffer-string)))


(defmethod run ((calc Jasp))
  "Write out input files, and run VASP as a simple shell command"
  (write-poscar calc)
  (write-kpoints calc)
  (write-potcar calc)
  (write-incar calc)
  (shell-command "vasp"))


(defmethod update ((calc Jasp))
  "Run vasp if needed for CALC.
We just check for a properly ended OUTCAR."
  (with-current-directory
   (oref calc wd)
   (unless (and (file-exists-p "OUTCAR")
                (with-temp-buffer
                  (insert-file-contents "OUTCAR")
                  (re-search-forward
                  "                 Voluntary context switches:"
                  (point-max)
                  t)))
     (run calc))))


(defmethod get-potential-energy ((calc Jasp))
  "Get potential energy from CALC."
  (update calc)
  (with-current-directory
   (oref calc wd)
   (with-temp-buffer
     (insert-file-contents "OUTCAR")
     ;; go to last entry
     (while (re-search-forward
             "free  energy   TOTEN  =\\s-*\\([-]?[0-9]*\\.[0-9]*\\) eV"
             (point-max)
             t)
       nil)
     ;; return last match
     (string-to-number  (match-string 1)))))

(provide 'jasp)

get-potential-energy

This is worth some discussion. On one hand, the constructor is a bit more verbose than the implementation in Python. In Python we use a context handler in place of the macro here. On the other hand, that verbosity comes with detailed, accessible documentation for each argument. We only considered the simplest of input arguments. It might be trickier to include lists, and other types of input. But I think those can all be worked out like they were in ase. We only implemented the simplest job control logic, but that also can be worked out. The biggest challenge might be getting more complex data from the output. Nearly everything can be obtained from vasprun.xml also, in the event that parsing is to slow or difficult.

Now, let us test this out. We can make a calculator: 

(setq calc (Jasp
            nil
            :xc "PBE"
            :encut 350
            :nbands 6
            :atoms (Atoms
                    nil
                    :atoms `(,(Atom nil :symbol "C" :x 0 :y 0 :z 0)
                             ,(Atom nil :symbol "O" :x 1.1 :y 0 :z 0))
                    :cell '(8 9 10))))

[object Jasp nil "." 350 6 (1 1 1) "PBE" [object Atoms nil ([object Atom nil "C" 0 0 0] [object Atom nil "O" 1.1 0 0]) (8 9 10)]]

We can call the class functions like this. Here we write out the corresponding POSCAR: 

(write-poscar calc)

Created by jasp.el
  1.0
 8.000000  0.000000  0.000000 
 0.000000  9.000000  0.000000 
 0.000000  0.000000  10.000000 
1 1 
Cartesian
0.000000 0.000000 0.000000
1.100000 0.000000 0.000000

It looks a little backward if you have only seen Python, where this would be calc.writeposcar(). It is almost the same characters, just a different order (and no . in lisp)!

Here we get the KPOINTS file: 

(write-kpoints calc)

Automatic mesh
0
Monkhorst-Pack
   1    1    1 
0.0 0.0 0.0

I cannot show the POTCAR file for licensing reasons, but it works. 

(write-potcar calc)

and the INCAR file: 

(write-incar calc)

ENCUT = 350.000000
NBANDS = 6

We run a calculation like this. This will run vasp directly (not through the queue system). 

(run calc)

0

The returned 0 means the shell command finished correctly.

And we retrieve the potential energy like this: 

(get-potential-energy calc)

-14.687906

Not bad. That is close to the result we got from a similar calculation here .

4 Putting it all together to run calculations

If we put this all together the way we might use it in practice, it looks like this. 

(load-file "Atom.el")
(load-file "Atoms.el")
(load-file "Jasp.el")

(let* ((co (Atoms
            nil
            :atoms `(,(Atom nil :symbol "C" :x 0 :y 0 :z 0)
                     ,(Atom nil :symbol "O" :x 1.1 :y 0 :z 0))
            :cell '(8 9 10)))

       (calc (Jasp
              nil
              :xc "PBE"
              :nbands 6
              :encut 350
              :atoms co)))
  
  (get-potential-energy calc))

-14.687906

Compare this with this Python code which does approximately the same thing: 

from ase import Atoms, Atom
from jasp import *

co = Atoms([Atom('C', [0,   0, 0]),
            Atom('O', [1.1, 0, 0])],
            cell=(6., 6., 6.))

with jasp('molecules/simple-co', #output dir
          xc='PBE',  # the exchange-correlation functional
          nbands=6,  # number of bands
          encut=350, # planewave cutoff
          atoms=co) as calc:
    print 'energy = {0} eV'.format(co.get_potential_energy())

They look pretty similar. One thing clearly missing from emacs-lisp that Python has is full support for numerics and plotting. Some of this could be addressed via Pymacs , but certainly not all of it. Some of it could also be handled using org-mode to enable data from emacs-lisp to go to other code blocks that can handle it.

Finally, for a little fun, we illustrate mapping over a range of bond lengths. There is more than one way to do this. For example, we could create a list of calculators, and then run over them. Here we create one calculator, and just change the x position in a loop. We use the more general setf approach instead of oset to see what it looks like. 

(let* ((co (Atoms
            nil
            :atoms `(,(Atom nil :symbol "C" :x 0 :y 0 :z 0)
                     ,(Atom nil :symbol "O" :x 1.1 :y 0 :z 0))
            :cell '(8 9 10)))
       (calc (Jasp
              nil
              :wd nil
              :xc "PBE"
              :nbands 6
              :encut 350
              :atoms co)))
  (dolist (d '(1.05 1.1 1.15 1.2 1.25))
    ;; change working directory
    (setf (oref calc wd) (format "co-%s" d))
    ;; set x-coordinate on oxygen atom
    (setf (oref (elt (oref co atoms) 1) x) d)
    (print (format "d = %s\nEnergy = %s eV"
                   d
                   (get-potential-energy calc)))))

"d = 1.05
Energy = -14.195892 eV"

"d = 1.1
Energy = -14.698456 eV"

"d = 1.15
Energy = -14.814809 eV"

"d = 1.2
Energy = -14.660395 eV"

"d = 1.25
Energy = -14.319904 eV"

See http://kitchingroup.cheme.cmu.edu/dft-book/dft.html#sec-3-4-1 for how this was done in Python. It looks pretty similar to me.

5 Summary thoughts

We implemented a bare bones emacs-lisp calculator for VASP. The library automates creation of input files, running the calculation, and parsing the output.

It seems pretty feasible to implement a pretty complete interface to VASP in emacs-lisp. The main reasons to do this are:

  1. Integrated access to documentation
  2. Emacs editing of emacs-lisp code
  3. Integration with org-mode

Even with Python editor that had access to documentation as deeply integrated as emacs has with emacs-lisp, it would still just be a Python editor, i.e. you probably could not use the editor to write org-mode, LaTeX, etc… It is time to recognize we need both scientific document creation and code editing capability in the same place! This kind of suggests a need to get a better Python environment going in Emacs, which deeper integration of the documentation. See this post for some progress in that area!

Copyright (C) 2014 by John Kitchin. See the License for information about copying.

org-mode source

Org-mode version = 8.2.10

Discuss on Twitter
« Previous Page