The Kitchin Research Group

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