A prototype implementation of jasp in emacs-lisp

| categories: vasp, emacs, lisp, ase | 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)

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))

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)
[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
   ((= 1 (length args))
    (elt (oref atoms atoms) (car args)))
    (mapcar (lambda (i)
              (elt (oref atoms atoms) i))

(provide 'Atoms)

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")
[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."
     (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)))) 

(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)     
            (insert "\n")
            (loop for j below (length cell)
                  (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)
           (format "%f %f %f\n"
                   (oref atom x)
                   (oref atom y)
                   (oref atom z))))

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

(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)
              (concat "potpaw_" xc)
              (oref atom symbol)

(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)))

(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."
   (oref calc wd)
   (unless (and (file-exists-p "OUTCAR")
                  (insert-file-contents "OUTCAR")
                  "                 Voluntary context switches:"
     (run calc))))

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

(provide 'jasp)

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
            :xc "PBE"
            :encut 350
            :nbands 6
            :atoms (Atoms
                    :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
 8.000000  0.000000  0.000000 
 0.000000  9.000000  0.000000 
 0.000000  0.000000  10.000000 
1 1 
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
   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

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

(run calc)

The returned 0 means the shell command finished correctly.

And we retrieve the potential energy like this:

(get-potential-energy calc)

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
            :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
              :xc "PBE"
              :nbands 6
              :encut 350
              :atoms co)))
  (get-potential-energy calc))

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
            :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
              :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"
                   (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