The Kitchin Research Group

The ase development team has created a new database module (https://wiki.fysik.dtu.dk/ase/ase/db/db.html ). In this post, we will examine how to use that in conjunction with jasp to save calculations to a database. To use this you need the latest version of jasp from http://github.com/jkitchin/jasp .

The idea here is to use the code from http://kitchingroup.cheme.cmu.edu/blog/2014/03/20/Finding-VASP-calculations-in-a-directory-tree/ to add each calculation in that subtree to an ase database. We will use the json format for the database.

import os
from jasp import *
from ase.db import connect
c = connect('vaspdb.json')

def vasp_p(directory):
    'returns True if a finished OUTCAR file exists in the current directory, else False'
    outcar = os.path.join(directory, 'OUTCAR')
    if os.path.exists(outcar):
        with open(outcar, 'r') as f:
            contents = f.read()
            if 'General timing and accounting informations for this job:' in contents:
                return True
    return False
            
        
total_time = 0

for root, dirs, files in os.walk('/home-research/jkitchin/research/rutile-atat'):
    for d in dirs:
        # compute absolute path to each directory in the current root
        absd = os.path.join(root, d)
        if vasp_p(absd):
            # we found a vasp directory, so we can do something in it. 
            # here we get the elapsed time from the calculation
            with jasp(absd) as calc:
                atoms = calc.get_atoms()
                calc.results['energy'] = atoms.get_potential_energy()
                calc.results['forces'] = atoms.get_forces()
            # it is important that this line not be inside the jasp
            # context-manager, because c.write writes in the local
            # dir.
            print c.write(atoms, ['atat']), absd

1 /home-research/jkitchin/research/rutile-atat/RuTi_O_rutile/0
2 /home-research/jkitchin/research/rutile-atat/RuTi_O_rutile/1
3 /home-research/jkitchin/research/rutile-atat/RuTi_O_rutile/10
4 /home-research/jkitchin/research/rutile-atat/RuTi_O_rutile/11
5 /home-research/jkitchin/research/rutile-atat/RuTi_O_rutile/12
6 /home-research/jkitchin/research/rutile-atat/RuTi_O_rutile/13
7 /home-research/jkitchin/research/rutile-atat/RuTi_O_rutile/14
8 /home-research/jkitchin/research/rutile-atat/RuTi_O_rutile/15
9 /home-research/jkitchin/research/rutile-atat/RuTi_O_rutile/16
10 /home-research/jkitchin/research/rutile-atat/RuTi_O_rutile/17
11 /home-research/jkitchin/research/rutile-atat/RuTi_O_rutile/2
12 /home-research/jkitchin/research/rutile-atat/RuTi_O_rutile/3
13 /home-research/jkitchin/research/rutile-atat/RuTi_O_rutile/4
14 /home-research/jkitchin/research/rutile-atat/RuTi_O_rutile/5
15 /home-research/jkitchin/research/rutile-atat/RuTi_O_rutile/59
16 /home-research/jkitchin/research/rutile-atat/RuTi_O_rutile/6
17 /home-research/jkitchin/research/rutile-atat/RuTi_O_rutile/66
18 /home-research/jkitchin/research/rutile-atat/RuTi_O_rutile/7
19 /home-research/jkitchin/research/rutile-atat/RuTi_O_rutile/73
20 /home-research/jkitchin/research/rutile-atat/RuTi_O_rutile/74
21 /home-research/jkitchin/research/rutile-atat/RuTi_O_rutile/78
22 /home-research/jkitchin/research/rutile-atat/RuTi_O_rutile/8
23 /home-research/jkitchin/research/rutile-atat/RuTi_O_rutile/9

The result is in this file: vaspdb.json .

We can see the contents of that database like this:

import os
from jasp import *
from ase.db import connect
c = connect('vaspdb.json')

g = c.select()  # g is a generator
print g.next()

{u'username': u'jkitchin', u'calculator_name': u'vasp', u'tags': array([0, 0, 0, 0, 0, 0]), u'positions': array([[ 1.382,  1.382,  0.   ],
       [ 3.145,  3.145,  0.   ],
       [ 3.646,  0.881,  1.548],
       [ 0.881,  3.646,  1.548],
       [ 0.   ,  0.   ,  0.   ],
       [ 2.263,  2.263,  1.548]]), u'energy': -44.251496, u'calculator_parameters': {u'incar': {u'doc': u'INCAR parameters', u'prec': u'high', u'nbands': 43, u'sigma': 0.1, u'encut': 350.0}, u'doc': u'JSON representation of a VASP calculation.\n\nenergy is in eV\nforces are in eV/\\AA\nstress is in GPa (sxx, syy, szz,  syz, sxz, sxy)\nmagnetic moments are in Bohr-magneton\nThe density of states is reported with E_f at 0 eV.\nVolume is reported in \\AA^3\nCoordinates and cell parameters are reported in \\AA\n\nIf atom-projected dos are included they are in the form:\n{ados:{energy:data, {atom index: {orbital : dos}}}\n', u'potcar': [[u'Ru', u'potpaw_PBE/Ru/POTCAR', u'dee616f2a1e7a5430bb588f1710bfea3001d54ea'], [u'O', u'potpaw_PBE/O/POTCAR', u'9a0489b46120b0cad515d935f44b5fbe3a3b1dfa']], u'input': {u'kpts': array([ 6,  6, 10]), u'reciprocal': False, u'xc': u'PBE', u'kpts_nintersections': None, u'setups': {}, u'txt': u'-', u'gamma': True}, u'atoms': {u'cell': array([[ 4.527,  0.   ,  0.   ],
       [ 0.   ,  4.527,  0.   ],
       [ 0.   ,  0.   ,  3.095]]), u'symbols': [u'O', u'O', u'O', u'O', u'Ru', u'Ru'], u'tags': array([0, 0, 0, 0, 0, 0]), u'pbc': array([ True,  True,  True], dtype=bool), u'positions': array([[ 1.382,  1.382,  0.   ],
       [ 3.145,  3.145,  0.   ],
       [ 3.646,  0.881,  1.548],
       [ 0.881,  3.646,  1.548],
       [ 0.   ,  0.   ,  0.   ],
       [ 2.263,  2.263,  1.548]])}, u'data': {u'stress': array([ 0.088,  0.088,  0.06 , -0.   , -0.   , -0.   ]), u'doc': u'Data from the output of the calculation', u'volume': 63.43221074143486, u'total_energy': -44.251496, u'forces': array([[-0.024, -0.024,  0.   ],
       [ 0.024,  0.024,  0.   ],
       [-0.024,  0.024,  0.   ],
       [ 0.024, -0.024,  0.   ],
       [ 0.   ,  0.   ,  0.   ],
       [ 0.   ,  0.   ,  0.   ]]), u'fermi_level': 5.0374}, u'metadata': {u'date.created': 1395241327.477995, u'uuid': u'7081ee4a-af77-11e3-a6e6-003048f5e49e', u'date.created.ascii': u'Wed Mar 19 11:02:07 2014', u'user.username': u'jkitchin', u'atoms.resort': array([2, 3, 4, 5, 0, 1]), u'user.email': u'jkitchin@andrew.cmu.edu', u'user.fullname': u'John Kitchin', u'O.potential.git_hash': u'9a0489b46120b0cad515d935f44b5fbe3a3b1dfa', u'atoms.tags': array([0, 0, 0, 0, 0, 0]), u'O.potential.path': u'potpaw_PBE/O/POTCAR', u'Ru.potential.path': u'potpaw_PBE/Ru/POTCAR', u'Ru.potential.git_hash': u'dee616f2a1e7a5430bb588f1710bfea3001d54ea'}}, u'cell': array([[ 4.527,  0.   ,  0.   ],
       [ 0.   ,  4.527,  0.   ],
       [ 0.   ,  0.   ,  3.095]]), u'numbers': array([ 8,  8,  8,  8, 44, 44]), u'pbc': array([ True,  True,  True], dtype=bool), u'timestamp': 14.23343757848325, u'keywords': [u'atat'], u'forces': array([[-0.024, -0.024,  0.   ],
       [ 0.024,  0.024,  0.   ],
       [-0.024,  0.024,  0.   ],
       [ 0.024, -0.024,  0.   ],
       [ 0.   ,  0.   ,  0.   ],
       [ 0.   ,  0.   ,  0.   ]]), 'id': 1, u'unique_id': u'123901e31734f14418381a23d1ee1072'}

The data stored there comes from the calc.todict() function written for jasp.

We can do some searches like this. Say we want to find all the calculations where there are four oxygen atoms.

import os
from jasp import *
from ase.db import connect
c = connect('vaspdb.json')

g = c.select(O=4)  # g is a generator
for entry in g:
    print c.get_atoms(entry['id']), '\n'

Atoms(symbols='O4Ru2', positions=..., tags=...,
      cell=[4.526933343669885, 4.526933343669885, 3.095292162609941],
      pbc=[True, True, True], calculator=SinglePointCalculator(...)) 

Atoms(symbols='O4Ti2', positions=..., tags=...,
      cell=[4.614336091353763, 4.614336091353763, 2.9555779409837473],
      pbc=[True, True, True], calculator=SinglePointCalculator(...)) 

Atoms(symbols='O4RuTi', positions=..., tags=...,
      cell=[[-0.0151920891931803, -4.604112035041115, 0.0],
      [-4.604112035041115, -0.0151920891931803, 0.0], [0.0, 0.0,
      -3.0110141191854245]], pbc=[True, True, True],
      calculator=SinglePointCalculator(...))

You can see there are three calculations with that criterion.

ase-db --help

usage: ase-db [-h] [-n] [-c COLUMNS] [--explain] [-y] [-i INSERT_INTO]
              [-k ADD_KEYWORDS] [-K ADD_KEY_VALUE_PAIRS]
              [--delete-keywords DELETE_KEYWORDS]
              [--delete-key-value-pairs DELETE_KEY_VALUE_PAIRS] [--delete]
              [-v] [-q] [-s SORT] [-r] [-l] [--limit LIMIT]
              [-p PYTHON_EXPRESSION]
              name [selection]

positional arguments:
  name
  selection

optional arguments:
  -h, --help            show this help message and exit
  -n, --count
  -c COLUMNS, --columns COLUMNS
                        short/long+row-row
  --explain
  -y, --yes
  -i INSERT_INTO, --insert-into INSERT_INTO
  -k ADD_KEYWORDS, --add-keywords ADD_KEYWORDS
  -K ADD_KEY_VALUE_PAIRS, --add-key-value-pairs ADD_KEY_VALUE_PAIRS
  --delete-keywords DELETE_KEYWORDS
  --delete-key-value-pairs DELETE_KEY_VALUE_PAIRS
  --delete
  -v, --verbose
  -q, --quiet
  -s SORT, --sort SORT
  -r, --reverse
  -l, --long
  --limit LIMIT
  -p PYTHON_EXPRESSION, --python-expression PYTHON_EXPRESSION

ase-db vaspdb.json id=1

id|age|user    |formula|calc| energy| fmax|pbc|  size|keywords|   mass
 1|9m |jkitchin|O4Ru2  |vasp|-44.251|0.033|111|63.432|atat    |266.138

ase-db vaspdb.json O=4

id|age|user    |formula|calc| energy| fmax|pbc|  size|keywords|   mass
 1|9m |jkitchin|O4Ru2  |vasp|-44.251|0.033|111|63.432|atat    |266.138
 2|9m |jkitchin|O4Ti2  |vasp|-52.970|0.033|111|62.930|atat    |159.758
11|9m |jkitchin|O4RuTi |vasp|-48.115|0.157|111|63.826|atat    |212.948

ase-db vaspdb.json atat

id|age|user    |formula  |calc|  energy| fmax|pbc|   size|keywords|   mass
 1|9m |jkitchin|O4Ru2    |vasp| -44.251|0.033|111| 63.432|atat    |266.138
 2|9m |jkitchin|O4Ti2    |vasp| -52.970|0.033|111| 62.930|atat    |159.758
 3|9m |jkitchin|O8Ru2Ti2 |vasp| -96.601|0.086|111|126.719|atat    |425.895
 4|9m |jkitchin|O8RuTi3  |vasp|-100.842|0.075|111|126.846|atat    |372.705
 5|9m |jkitchin|O8Ru3Ti  |vasp| -92.376|0.133|111|127.420|atat    |479.085
 6|9m |jkitchin|O8Ru2Ti2 |vasp| -96.594|0.184|111|127.176|atat    |425.895
 7|9m |jkitchin|O8RuTi3  |vasp|-100.959|0.176|111|126.924|atat    |372.705
 8|9m |jkitchin|O8Ru3Ti  |vasp| -92.314|0.084|111|127.377|atat    |479.085
 9|9m |jkitchin|O8Ru2Ti2 |vasp| -96.612|0.086|111|126.542|atat    |425.895
10|9m |jkitchin|O8RuTi3  |vasp|-100.816|0.080|111|126.557|atat    |372.705
11|9m |jkitchin|O4RuTi   |vasp| -48.115|0.157|111| 63.826|atat    |212.948
12|9m |jkitchin|O8Ru3Ti  |vasp| -92.429|0.163|111|127.291|atat    |479.085
13|9m |jkitchin|O8Ru2Ti2 |vasp| -96.770|0.166|111|126.870|atat    |425.895
14|9m |jkitchin|O8RuTi3  |vasp|-101.014|0.222|111|126.881|atat    |372.705
15|9m |jkitchin|O12Ru4Ti2|vasp|-140.969|0.114|111|190.614|atat    |692.033
16|9m |jkitchin|O8Ru3Ti  |vasp| -92.323|0.125|111|127.541|atat    |479.085
17|9m |jkitchin|O12Ru2Ti4|vasp|-149.516|0.241|111|190.070|atat    |585.653
18|9m |jkitchin|O8Ru2Ti2 |vasp| -96.661|0.064|111|127.038|atat    |425.895
19|9m |jkitchin|O12Ru4Ti2|vasp|-140.472|0.138|111|190.640|atat    |692.033
20|9m |jkitchin|O12Ru3Ti3|vasp|-144.667|0.166|111|190.604|atat    |638.843
21|9m |jkitchin|O12Ru2Ti4|vasp|-148.813|0.055|111|190.084|atat    |585.653
22|9m |jkitchin|O8RuTi3  |vasp|-100.874|0.051|111|126.690|atat    |372.705
23|9m |jkitchin|O8Ru3Ti  |vasp| -92.246|0.102|111|127.383|atat    |479.085

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

org-mode source

Org-mode version = 8.2.5h

Today someone asked about deleting multiple elements from a list (actually it was about deleting multiple atoms from an ase.Atoms object, but some principles here apply. I will address that actual question later.).

Deleting multiple items from a list is not directly possible in one command in Python. There are a few approaches to accomplishing something like it. Which one is best depends on your objective.

One problem is when you delete an item, the indices of every item after it also changes. One strategy then is to delete the elements in descending order, i.e. delete the largest indices first. That way, you do not change the indices of the smaller indices, so you can still delete them. We can sort them in reverse order like this:

a = [1, 2, 5, 6, 7]

ind2remove = [1, 3]

for i in sorted(ind2remove, reverse=True): 
    del a[i]

print a

[1, 5, 7]

An alternative approach is to make a new list that only has the elements you want using list comprehension. For example:

a = [1, 2, 5, 6, 7]

ind2remove = [1, 3]

a = [x for i,x in enumerate(a) if i not in ind2remove]

print a

[1, 5, 7]

With numpy arrays you can delete multiple elements like this:

import numpy as np

a = np.array([1, 2, 5, 6, 7])

ind2remove = [1, 3]

print np.delete(a, ind2remove)
print a

[1 5 7]
[1 2 5 6 7]

The delete command makes a new object; the original list is unchanged. Numpy arrays are technically immutable, so the only way to do this is to make a copy. Another way is to use a boolean mask that only selects the indices where the mask is True, and not where they are False.

import numpy as np

a = np.array([1, 2, 5, 6, 7])

ind2remove = [1, 3]

mask = np.ones(len(a), dtype=bool) 
mask[ind2remove] = False
print a[mask]
print a

[1 5 7]
[1 2 5 6 7]

There might be other ways to do this too. These examples are nearly indistinguishable for small lists. For very large lists (I guess 1000's of elements), you may find one method more efficient than the others.

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

org-mode source

Org-mode version = 8.2.5h

I am exploring the possibility of using tags on python objects in conjunction with searches to find sets of objects. Here I want to explore some syntax and methods for doing that.

In org-mode there is a syntax like '+boss+urgent-project1' for and and not operators and 'A|B' for or operators. I think we need pyparsing to untangle this kind of syntax. See http://pyparsing.wikispaces.com/file/view/simpleBool.py for an example. Another alternative might be the natural language toolkit (nltk ). Before we dig into those, let us see some python ways of doing the logic.

Below we define some lists containing tags (strings). We

a = ['A', 'B', 'C']

b = ['A', 'B']

c = ['A', 'C']

d = [ 'B', 'C']

all_lists = [a, b, c, d]

# get functions with tags A and B
print 'A and B ',[x for x in all_lists if ('A' in x) & ('B' in x)]

# A not B
print 'A not B ',[x for x in all_lists if ('A' in x) & ('B' not in x)]

# B or C
print 'B or C ', [x for x in all_lists if ('B' in x) | ('C' in x)]

# B or C but not both
print 'B xor C ',[x for x in all_lists if ('B' in x) ^ ('C' in x)]

A and B  [['A', 'B', 'C'], ['A', 'B']]
A not B  [['A', 'C']]
B or C  [['A', 'B', 'C'], ['A', 'B'], ['A', 'C'], ['B', 'C']]
B xor C  [['A', 'B'], ['A', 'C']]

Those are not too bad. Somehow I would have to get pyparsing to generate that syntax. That will take a lot of studying. There are some other ways to do this too. Let us try that out with itertools.

a = ['A', 'B', 'C']

b = ['A', 'B']

c = ['A', 'C']

d = [ 'B', 'C']

all_lists = [a, b, c, d]

import itertools as it

# ifilter returns an iterator
print 'A and B ', list(it.ifilter(lambda x: ('A' in x) & ('B' in x), all_lists))

A and B  [['A', 'B', 'C'], ['A', 'B']]

I do not like this syntax better. The iterator is lazy, so we have to wrap it in a list to get the results. Eventually, I want to do something like these:

filter('A and B', all_lists)
A or B
A xor B
not A and B
not(A and B)

I think that calls for pyparsing. I think the syntax above is better (more readable) than this:

filter('A & B', all_lists)
A | B
A ^ B
~A & B
~(A & B)

It is not that obvious though how to get from that syntax to the code I illustrated above though.

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

org-mode source

Org-mode version = 8.2.5h

This fellowship was established with the generous donation of Dr. H. William and Ruth Prengle to provide tuition for graduate students with outstanding academic records.

Congratulations Mehak!

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

org-mode source

Org-mode version = 8.2.5h

It would be helpful to get a listing of org-files in a directory tree in the form of clickable links. This would be useful, for example, to find all files associated with a project in a directory with a particular extension, or to do some action on all files that match a pattern. To do this, we will have to recursively walk through the directories and examine their contents.

Let us examine some of the commands we will need to use. One command is to get the contents of a directory. We will explore the contents of a directory called literate in my computer.

;; list contents of the directory
(let ((abspath nil)
      (match nil)
      (nosort t))
  (directory-files "literate" abspath match nosort))

Note the presence of . and ... Those stand for current directory and one directory up. We should remove those from the list. We can do that like this.

;; remove . and ..
(let ((abspath nil)
      (match nil)
      (nosort t))
  (remove "." 
          (remove ".." 
                  (directory-files "literate" abspath match nosort))))

Next, we need to know if a given entry in the directory files is a file or a directory. Emacs-lisp has a few functions for that. We use absolute filenames here since the paths are relative to the "molecules" directory. Note we could use absolute paths in directory-files, but that makes it hard to remove "." and "..".

;; print types of files in the directory
(let ((root "literate")
      (abspath nil)
      (match nil)
      (nosort t))
  (mapcar (lambda (x)
            (cond
             ((file-directory-p (expand-file-name x root))
              (print (format "%s is a directory" x)))
             ((file-regular-p (expand-file-name x root))
              (print (format "%s is a regular file" x)))))
          (remove "." 
                  (remove ".." 
                          (directory-files root abspath match nosort)))))

"makefile-main is a regular file"

"Makefile is a regular file"

"main.o is a regular file"

"main.f90 is a regular file"

"main is a regular file"

"literate.org is a regular file"

"hello.f90 is a regular file"

"circle.o is a regular file"

"circle.mod is a regular file"

"circle.f90 is a regular file"

"circle-area.png is a regular file"

"archive is a directory"

"a.out is a regular file"

Now, we are at the crux of this problem. We can differentiate between files and directories. For each directory in this directory, we need to recurse into it, and list the contents. There is some code at http://turingmachine.org/bl/2013-05-29-recursively-listing-directories-in-elisp.html which does this, but I found that I had to modify the code to not list directories, and here I want to show a simpler recursive code.

(defun os-walk (root)
  "recursively walks through directories getting list of absolute paths of files"
  (let ((files '()) ; empty list to store results
        (current-list (directory-files root t)))
    ;;process current-list
    (while current-list
      (let ((fn (car current-list))) ; get next entry
        (cond 
         ;; regular files
         ((file-regular-p fn)
          (add-to-list 'files fn))
         ;; directories
         ((and
           (file-directory-p fn)
           ;; ignore . and ..
           (not (string-equal ".." (substring fn -2)))
           (not (string-equal "." (substring fn -1))))
          ;; we have to recurse into this directory
          (setq files (append files (os-walk fn))))
        )
      ;; cut list down by an element
      (setq current-list (cdr current-list)))
      )
    files))

(os-walk "literate")

Nice, that gives us a recursive listing of all the files in this directory tree. Let us take this a step further, and apply a function to that list to filter out a list of the org files. We will also create org-links out of these files.

(defun os-walk (root)
  (let ((files '()) ;empty list to store results
        (current-list (directory-files root t)))
    ;;process current-list
    (while current-list
      (let ((fn (car current-list))) ; get next entry
        (cond 
         ;; regular files
         ((file-regular-p fn)
          (add-to-list 'files fn))
         ;; directories
         ((and
           (file-directory-p fn)
           ;; ignore . and ..
           (not (string-equal ".." (substring fn -2)))
           (not (string-equal "." (substring fn -1))))
          ;; we have to recurse into this directory
          (setq files (append files (os-walk fn))))
        )
      ;; cut list down by an element
      (setq current-list (cdr current-list)))
      )
    files))

(require 'cl)

(mapcar 
 (lambda (x) (princ (format "[[%s][%s]]\n" x (file-relative-name x "."))))
 (remove-if-not 
  (lambda (x) (string= (file-name-extension x) "org"))
  (os-walk "literate")))

literate/literate.org

That is certainly functional. It might be nice to format the links a bit nicer to show their structure in a table of contents way, or to sort them in a nice order if there were many of these files.

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

org-mode source

Org-mode version = 8.2.5h