yasnippets for jasp, ase and python

| categories: ase, emacs, jasp | tags: | View Comments

In using [[http://github.com/jkitchin/jasp for calculations, I find there are lots of small python phrases I use over and over. Today I will examine using yasnippet to save time and keystrokes. yasnippet is a template expansion module, where you type a small set of characters, press Tab, and the characters "expand" to the full text. It is pretty sophisticated, and allows you to define "tab-stops" which you interactively fill in, and tab between like filling in a form.

All the snippets are defined in the

*Appendix
.

1 Tangle the snippets, and add them to yasnippet

Each snippet definition belongs in a file in a directory. The main directory is called "snippets". Since I anticipate using these snippets in org-mode, each snippet is defined in a directory within snippets called "org-mode". First, we make the directory here. I also want to use the snippets in python mode, so we also create a python-mode directory here. We do not have to duplicate the snippets. We can create a file called .yas-parents , with one line in it containing "org-mode".

mkdir -p snippets/org-mode
mkdir -p snippets/python-mode
echo "org-mode" > snippets/python-mode/.yas-parents

Each snippet is defined in a src block with a :tangle header. So, we can extract them all in one command here.

(org-babel-tangle)
snippets/org-mode/iase snippets/org-mode/imp snippets/org-mode/inp snippets/org-mode/ij snippets/org-mode/pl snippets/org-mode/pyl snippets/org-mode/pxl snippets/org-mode/pp snippets/org-mode/npa snippets/org-mode/awt snippets/org-mode/avw snippets/org-mode/agf snippets/org-mode/ape snippets/org-mode/atms snippets/org-mode/atm snippets/org-mode/cga snippets/org-mode/cc snippets/org-mode/wjn snippets/org-mode/wjl

We also need to add our new directory to yasnippets. This is done by adding the directory to the yas-snippet-dirs variable. You could add this to your init.el file to permanently add these snippets.

(add-to-list 'yas-snippet-dirs "c:/Users/jkitchin/Dropbox/blogofile-jkitchin.github.com/_blog/snippets")
c:/Users/jkitchin/Dropbox/blogofile-jkitchin.github.com/blog/snippets ~/.emacs.d/snippets c:/users/jkitchin/Dropbox/kitchingroup/jmax/elpa/yasnippet-20140106.1009/snippets

Finally, we reload all the snippet definitions, so our new definitions are ready to use.

(yas-reload-all)
[yas] Reloaded everything (snippets will load just-in-time)... (some errors, check *Messages*).

Alternatively, you might just load this directory.

(yas-load-directory "./snippets")

2 Using the snippets

Each of these snippets is for a python phrase, but I usually write my python blocks in org-mode. You would use these by typing the shortcut name, and then pressing tab. Below I show what each shortcut expands to.

wjl → with jasp('') as calc:

wjn → with jasp('',) as calc: calc.calculate(atoms)

cc → calc.calculate(atoms)

cga → atoms = calc.get_atoms()

atm → Atom('', )

atms → atoms = Atoms([], cell)=

ape → atoms.get_potential_energy()

agf → atoms.get_forces()

avw → from ase.visualize import view view(atoms)

awt → from ase.io import write write('.png', atoms, show_unit_cell=2)

npa → np.array()

pp → plt.plot(, )

pxl → plt.xlabel()

pyl → plt.ylabel()

pl → plt.legend()

ij → from jasp import *

inp → import numpy as np

imp → import matplotlib.pyplot as plt

iase → from ase import Atom, Atoms

What other snippets would be handy?

3 Appendix

3.1 jasp snippets

# -*- mode: snippet -*-
# --
with jasp('$1') as calc:
    $0
# -*- mode: snippet -*-
# --
with jasp('$1',$0) as calc:
    calc.calculate(atoms)
# -*- mode: snippet -*-
# --
calc.calculate(atoms)
# -*- mode: snippet -*-
# --
atoms = calc.get_atoms()

3.2 ase snippets

Template for an ase.Atom

# -*- mode: snippet -*-
# --
Atom('$1', $2)
# -*- mode: snippet -*-
# --
atoms = Atoms([$1], cell=$2)
# -*- mode: snippet -*-
# --
atoms.get_potential_energy()
# -*- mode: snippet -*-
# --
atoms.get_forces()
# -*- mode: snippet -*-
# --
from ase.visualize import view
view(${1:atoms})
# -*- mode: snippet -*-
# --
from ase.io import write
write('$1.png', ${2:atoms}, show_unit_cell=${3:2})

3.3 python snippets

# -*- mode: snippet -*-
# --
import numpy as np
# -*- mode: snippet -*-
# --
import matplotlib.pyplot as plt
# -*- mode: snippet -*-
# --
from ase import Atom, Atoms
# -*- mode: snippet -*-
# --
np.array($0)
# -*- mode: snippet -*-
# --
plt.plot($1, $2)
# -*- mode: snippet -*-
# --
plt.xlabel($1)
# -*- mode: snippet -*-
# --
plt.ylabel($1)
# -*- mode: snippet -*-
# --
plt.legend($1)
# -*- mode: snippet -*-
# --
from jasp import *

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

org-mode source

Org-mode version = 8.2.5h

Read and Post Comments

Sharing data via JSON files

| categories: ase, dft | tags: | View Comments

Table of Contents

In a previous post I discussed how to represent a single DFT calculation in a JSON file that could easily be shared, and reanalyzed. Here I look at sharing a series of calculations. I had previously run calculations to analyze an equation of state for Cu. Here we create a list of the JSON representations of each calculation, and save it in one overall JSON file. We will give the data some structure and documentation. JSON represents dictionaries very well, so we build a dictionary of the results and then "dump" that dictionary to a JSON file.

from jasp import *

LC = [3.5, 3.55, 3.6, 3.65, 3.7, 3.75]

data = {'doc':'JSON file containing a set of calculations for the equation of state of Cu'}

data['calculations'] = []

for a in LC:
    with jasp('../bulk/Cu-{0}'.format(a)) as calc:
        data['calculations'].append(calc.dict)

import json
with open('eos.json', 'wb') as f:
    f.write(json.dumps(data))

Now, you can view this eos.json file, and analyze it yourself as follows. In python we read the file in and convert it to a data structure using the json module.

import json

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

volumes =  [entry['data']['volume'] for entry in d['calculations']]
energies = [entry['data']['total_energy'] for entry in d['calculations']]

import matplotlib.pyplot as plt
plt.plot(volumes, energies)
plt.xlabel('Volume ($\AA$)')
plt.ylabel('Total energy (eV)')
plt.savefig('eos-from-json.png')
plt.show()

If you wanted to do further analysis you could. Suppose you wanted to know more detail about how that calculation was done? You can retrieve the INCAR, KPOINT, and POTCAR details for each calculation. Any parameter not listed here was not specified in the calculations.

import json

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

# print details for the first calculation in the list
print d['calculations'][0]['incar']
print d['calculations'][0]['input']
print d['calculations'][0]['potcar']
{u'doc': u'INCAR parameters', u'nbands': 9, u'encut': 350.0, u'prec': u'Normal'}
{u'kpts': [8, 8, 8], u'reciprocal': False, u'xc': u'PBE', u'kpts_nintersections': None, u'setups': None, u'txt': u'-', u'gamma': False}
[[u'Cu', u'potpaw_PBE/Cu/POTCAR', u'a44c591415026f53deb16a99ca3f06b1e69be10b']]

The POTCAR information contains the species, the path to the POTCAR, and a git-hash of the POTCAR file. That way you can tell if you used exactly the same file that I did. You can compute a git hash like this :

cat /home-research/jkitchin/src/vasp/potpaw_PBE/Cu/POTCAR | git hash-object --stdin
a44c591415026f53deb16a99ca3f06b1e69be10b

If you want to get the details of the geometry of some calculation, you do it this way:

import json

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

print d['calculations'][0]['atoms'].keys()
print 'cell = ', d['calculations'][0]['atoms']['cell']
print 'syms = ', d['calculations'][0]['atoms']['symbols']
print 'cpos = ', d['calculations'][0]['atoms']['positions']
[u'cell', u'symbols', u'positions', u'pbc', u'tags']
cell =  [[1.75, 1.75, 0.0], [0.0, 1.75, 1.75], [1.75, 0.0, 1.75]]
syms =  [u'Cu']
cpos =  [[0.0, 0.0, 0.0]]

You can do further analysis from there.

1 Summary

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.

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

org-mode source

Read and Post Comments

Serializing jasp calculations as json data

| categories: ase, jasp, vasp | tags: | View Comments

We use VASPto calculate materials properties in our research We use the jasppython module we have developed to setup, run and analyze those calculations. One of the things we have worked on developing recently is to more transparently share how do this kind of work by using org-mode supporting information files. Doing this should make our research more reproducible, and allow others to build off of it more easily.

We have run into the following problem trying to share VASP results however. The VASP license prohibits us from sharing the POTCAR files that are used to run the calculations. That is unfortunate, but since these files are also what give VASP some competitive advantage, they are protected, and we agreed to that when we bought the license. The problem is that the jasp module requires the POTCAR files to work, so without them, our scripts are not reproducible by researchers without a VASP license.

So, we have been looking at new ways to share the data from our calculations. In this post, we consider representing the calculation as a JSON file. We will look at a couple of new features built into the development branch of jasp

1 The simplest case of a simple calculation

Here we setup and run a simple calculation, and output the JSON file.

from ase import Atoms, Atom
from jasp import *
import numpy as np
np.set_printoptions(precision=3, suppress=True)

co = Atoms([Atom('C',[0,   0, 0]),
            Atom('O',[1.2, 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
          ismear=1,  # Methfessel-Paxton smearing
          sigma=0.01,# very small smearing factor for a molecule
          atoms=co) as calc:
    print 'energy = {0} eV'.format(co.get_potential_energy())
    print co.get_forces()
    with open('JSON', 'w') as f:
        f.write(calc.json)
energy = -14.687906 eV
[[ 5.095  0.     0.   ]
 [-5.095  0.     0.   ]]

Now, we can analyze the JSON file independently of jasp. The json data contains all the inputs we used for the VASP calculation, the atomic geometry, and many of the outputs of the calculation. Here is the JSONfile.

import json
with open('molecules/simple-co/JSON', 'rb') as f:
    d = json.loads(f.read())

print('The energy is {0}'.format(d['data']['total_energy']))
print('The forces are {0}'.format(d['data']['forces']))
The energy is -14.687906
The forces are [[5.095488, 0.0, 0.0], [-5.095488, 0.0, 0.0]]

2 Including extra information in the JSON file

If we use a slightly different syntax, we can also include the total DOS in the JSON file.

from jasp import *

with jasp('molecules/simple-co') as calc:
    with open('JSON-DOS', 'w') as f:
        f.write(calc_to_json(calc, dos=True))

To illustrate that we have done that, let us plot the DOS without using jasp from the JSON-DOSfile.

import json
import matplotlib.pyplot as plt

with open('molecules/simple-co/JSON-DOS', 'rb') as f:
    d = json.loads(f.read())

energies = d['data']['dos']['e']
dos = d['data']['dos']['dos']
plt.plot(energies, dos)
plt.savefig('molecules/simple-co/dos.png')

We are still working on getting atom-projected DOS into the json file, and ensuring that all the spin cases are handled (e.g. the spin-up and spin-down DOS).

3 Limitations?

JSON is flexible, and can store text and numeric data. It does not store numpy arrays, but rather it is limited to storing lists of data. You would have to convert them back to arrays if you want to do array math. You probably wouldn't want to store a 3d array of electron density in this format, although it probably isn't worse than a CUBE file format. We haven't tested these files very significantly yet at a large scale to see how fast it is to read from lots of them.

Nonetheless, this looks like a reasonable format to share data in human and machine readable form, without violating the VASP licence conditions.

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

org-mode source

Read and Post Comments

Subclassing an ase.calculators.vasp calculator to do post analysis

| categories: ase, vasp | tags: | View Comments

Recently someone reported on the ase-users mail list an issue they had with the Vasp calculator where the maximum number of electronic steps and then stops, even though the calculation is not converged. They asked how to make this raise an exception. This is not a built in feature in ase.calculators.vasp, but we can add it via a subclass. The idea is create a new class that inherits from ase.calculators.vasp, and then augment the calculate function to do some post analysis. In this case, we will parse the OUTCAR file looking for lines like:

----------------------------------------- Iteration    1(   1)  ---------------------------------------

I think the first integer corresponds to an ionic iteration, while the number in parentheses corresponds to an electronic iteration. We will just count them up, and see if they match the values we specified. If they do, we can print an error message, or raise an exception. Here is the code:

from ase.calculators.vasp import Vasp
from ase import Atom, Atoms
import os, re

os.chdir('h2o-relax')

class myvasp(Vasp):
    '''subclass to run VASP and then check OUTCAR to see if the number
    of electronic steps is equal to the max NELM specified, or if the
    number of geometry steps is equal to NSW. Either case indicates
    VASP may not be converged, but just stopped because it was told
    too.'''

    #save original function
    
    def calculate(self, *args):
        'Run VASP, then check if nsw or nelm is met'

        #run original calculate method
        Vasp.calculate(self, *args)
        
        # now do post analysis
        NELM = self.int_params['nelm']
        NSW = self.int_params['nsw']

        max_nelm, max_nsw = 0, 0

        # open OUTCAR and check for what you want
        # parse this line: - Iteration    1(   1)
        regexp = re.compile('Iteration\s+(?P<nsw>\d*)\(\s+(?P<nelm>\d*)\)')
        count = 0
        with open('OUTCAR') as f:
            for line in f:
                if 'Iteration' in line:
                    m = regexp.search(line)
                    nsw = int(m.groupdict().get('nsw'))
                    if nsw > max_nsw:
                        max_nsw = nsw
                    nelm = int(m.groupdict().get('nelm'))
                    if nelm > max_nelm:
                        max_nelm = nelm

        if max_nelm == NELM:
            # raise exception here if you prefer
            print ('{0} NELM steps taken! check convergence'.format(max_nelm))
        if max_nsw == NSW:
            print('{0} NSW steps taken! check convergence'.format(max_nsw))


atoms = Atoms([Atom('H', [0.5960812, -0.7677068, 0.0000000]),
               Atom('O', [0.0000000, 0.0000000, 0.0000000]),
               Atom('H', [0.5960812, 0.7677068, 0.0000000])],
               cell=(8, 8, 8))

calc = myvasp(xc='PBE',
              nelm=2,
              encut=400,
              ismear=0,# Gaussian smearing
              ibrion=2,
              ediff=1e-8,
              nsw=3)

atoms.set_calculator(calc)

print 'Forces'
print '==========================='
print atoms.get_forces()
Forces
===========================
 running on    1 nodes
 distr:  one band on    1 nodes,    1 groups
 vasp.5.2.12 11Nov11 complex                                                    
  
 POSCAR found :  2 types and       3 ions
 LDA part: xc-table for Pade appr. of Perdew
 found WAVECAR, reading the header
 POSCAR, INCAR and KPOINTS ok, starting setup
 WARNING: small aliasing (wrap around) errors must be expected
 FFT: planning ...(           1 )
 reading WAVECAR
 the WAVECAR file was read sucessfully
 initial charge from wavefunction
 entering main loop
       N       E                     dE             d eps       ncg     rms          rms(c)
DAV:   1    -0.142298749169E+02   -0.14230E+02   -0.10298E-01    24   0.330E+00    0.231E-01
DAV:   2    -0.142281671566E+02    0.17078E-02   -0.36659E-03    24   0.581E-01
   1 F= -.14228167E+02 E0= -.14228167E+02  d E =-.142282E+02
 curvature:   0.00 expect dE= 0.000E+00 dE for cont linesearch  0.000E+00
 trial: gam= 0.00000 g(F)=  0.114E-04 g(S)=  0.000E+00 ort = 0.000E+00 (trialstep = 0.100E+01)
 search vector abs. value=  0.114E-04
 bond charge predicted
       N       E                     dE             d eps       ncg     rms          rms(c)
DAV:   1    -0.142273929196E+02    0.24820E-02   -0.17613E-03    24   0.419E-01    0.149E-02
DAV:   2    -0.142274190765E+02   -0.26157E-04   -0.20256E-04    16   0.115E-01
   2 F= -.14227419E+02 E0= -.14227419E+02  d E =0.748080E-03
 trial-energy change:    0.000748  1 .order   -0.000006   -0.000011   -0.000001
 step:   1.0700(harm=  1.0700)  dis= 0.00047  next Energy=   -14.228173 (dE=-0.609E-05)
 bond charge predicted
       N       E                     dE             d eps       ncg     rms          rms(c)
DAV:   1    -0.142274279955E+02   -0.35076E-04   -0.72088E-06    32   0.308E-02    0.758E-03
DAV:   2    -0.142274385490E+02   -0.10553E-04   -0.97783E-07    24   0.858E-03
   3 F= -.14227439E+02 E0= -.14227439E+02  d E =0.728608E-03
 curvature:  -0.53 expect dE=-0.133E-03 dE for cont linesearch -0.931E-08
 trial: gam=21.74088 g(F)=  0.248E-03 g(S)=  0.000E+00 ort = 0.445E-06 (trialstep = 0.204E-02)
 search vector abs. value=  0.565E-02
 writing wavefunctions
Running vanilla serial job
2 NELM steps taken! check convergence
3 NSW steps taken! check convergence
[[ 0.024 -0.028  0.   ]
 [-0.048  0.     0.   ]
 [ 0.024  0.028  0.   ]]

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

org-mode source

Read and Post Comments

« Previous Page