The Kitchin Research Group

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

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

(org-babel-tangle)

(add-to-list 'yas-snippet-dirs "c:/Users/jkitchin/Dropbox/blogofile-jkitchin.github.com/_blog/snippets")

(yas-reload-all)

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

(yas-load-directory "./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()

# -*- 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})

# -*- 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

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.

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

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

org-mode source

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

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.   ]]

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

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

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

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

org-mode source

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