[ase-users] An Error on a NEB calculation with VASP

[kaoru.yamazaki at imr.tohoku.ac.jp]

Dear All,

I have a trouble to run a NEB calculation using ASE 3.16 and VASP 5.44. 
Even though I follow the NEB tutorial (https://wiki.fysik.dtu.dk/ase/tutorials/neb/diffusion.html),

 I received an error message as follows:

Traceback (most recent call last):
  File "ceo-b-mfi_vasp_neb_180420.py", line 115, in <module>
    qn.run(fmax=0.03)
  File "/home/users/uv5/anaconda2/lib/python2.7/site-packages/ase/
optimize/optimize.py", line 174, in run
    f = self.atoms.get_forces()
  File "/home/users/uv5/anaconda2/lib/python2.7/site-packages/ase/neb.py
", line 150, in get_forces
    raise ValueError(msg)
ValueError: One or more NEB images share the same calculator.  Each 
image must have its own calculator.  You may wish to use the ase.neb.
SingleCalculatorNEB class instead, although using separate calculators 
is recommended.

My ASE-NEB input is as follows:

from ase.neb import NEB
from ase.build import bulk
from ase import Atoms, Atom
from ase.optimize import BFGS
from ase.constraints import UnitCellFilter
from ase.calculators.vasp import Vasp
from ase.io import read, write
from ase.neb import NEBTools

#Configuration of vasp
## setup vasp for zeolite
calc_guess = Vasp(prec='normal',
            algo='fast',
            xc='PBE',
           encut= 550,
            setups ={'Ce': '_3'}, # use Ce^3+ PP
            ivdw=11,
            ldau      = 'TRUE',
            ldautype  = 2,
        ldau_luj={'Ce':{'L':3, 'U':4.5, 'J':0.0}, 
                   'O':{'L':-1, 'U':0.0, 'J':0.0}, 
                   'Si':{'L':-1, 'U':0.0, 'J':0.0}, 
                             'B':{'L':-1, 'U':0.0, 'J':0.0},
                             'C':{'L':-1, 'U':0.0, 'J':0.0}, 
                             'H':{'L':-1, 'U':0.0, 'J':0.0}},
            icharg=12,
            ldauprint = 2,
            lmaxmix   = 6, 
            ncore      = 32,
            npar=8,
            lplane = 'TRUE',
            lscalu = 'FALESE',
            nsim = 8,
            ispin=2,
            nbands = 2560,
            laechg= 'TRUE', #Use DFT-D3 and GGA+U, 32 core 
            lreal='AUTO')


calc_grad = Vasp(prec='normal',
            algo='fast',
            xc='PBE',
           encut= 550,
            setups ={'Ce': '_3'}, # use Ce^3+ PP
            ivdw=11,
            ldau      = 'TRUE',
            ldautype  = 2,
        ldau_luj={'Ce':{'L':3, 'U':4.5, 'J':0.0}, 
                   'O':{'L':-1, 'U':0.0, 'J':0.0}, 
                   'Si':{'L':-1, 'U':0.0, 'J':0.0}, 
                             'B':{'L':-1, 'U':0.0, 'J':0.0},
                             'C':{'L':-1, 'U':0.0, 'J':0.0}, 
                             'H':{'L':-1, 'U':0.0, 'J':0.0}},
            ldauprint = 2,
            lmaxmix   = 6, 
            ncore      = 32,
            npar=8,
            lplane = 'TRUE',
            lscalu = 'FALESE',
            nsim = 8,
            ispin=2,
            nelm=1000,
            nbands = 2560,
            laechg= 'TRUE', #Use DFT-D3 and GGA+U, 32 core 
            lreal='AUTO')

#Lattice vector (optimized one without surface adsobate)
Lattice_vector =[(20.090705234, 0.00135092434476, -0.0117978794464),
          (-0.00835128432838, 39.9536468604, -0.00596798669412),
          ( 0.0593384875412, -0.0040107523377, 13.3610212492)]

# initial structure
initial = read('initial.xyz')
initial.pbc = (True, True, True)
initial.set_cell(Lattice_vector)
initial.set_calculator(calc_guess)
initial.get_potential_energy()
initial.set_calculator(calc_grad)
initial.get_potential_energy()


#Create final state
final = read('final.xyz')
final.pbc = (True, True, True)
final.set_cell(Lattice_vector)

#Generate blank images
images = [initial]

for i in range(5):
    images.append(initial.copy())

for image in images:
    image.pbc = (True, True, True)
    image.set_cell(Lattice_vector)
    image.set_calculator(calc_grad)

images.append(final)

#Run IDPP interpolation
neb = NEB(images, climb=True,method='eb')
neb.interpolate('idpp')

#Run NEB calculation
qn = BFGS(neb, trajectory='neb_idpp.traj', logfile='neb_idpp.log',  
restart='neb.pckl')
qn.run(fmax=0.03)

# print neb results
images = read('neb_idpp.traj at -7:')
nebtools = NEBTools(images)
Ef, dE =  nebtools.get_barrier(fit=False, raw=True)
print Ef
print de
write('neb.xyz', neb)

I guess there is something wrong on the specification on the VASP 
calculator on each NEB image but I have not fount what is wrong. Could 
you kindly let me know possible solutions on this problem?

Best regards,
Kaoru Yamazaki, Ph.D.
Assistant Professor
Institute for Materials Research, 
Tohoku University
E-mail: kaoru.yamazaki at imr.tohoku.ac.jp