from ase import io
from ase.neb import NEB
from ase.optimize import MDMin
from gpaw import GPAW, PW, PoissonSolver
from ase.calculators.siesta import Siesta
from ase.units import Ry
# Read initial and final states:

initial = io.read('react.traj')
initial.cell /= 2
final = io.read('prod.traj')
final.cell /= 2
# Make a band consisting of 5 images:
images = [initial]
images += [initial.copy() for i in range(3)]
images += [final]
neb = NEB(images)

# Interpolate linearly the potisions of the three middle images:
neb.interpolate()

#calculator = GPAW(xc='oldLDA', mode='lcao', basis='sz(dzp)', h=0.3,
#                  poissonsolver=PoissonSolver(relax='GS'))

def getcalc():
    return GPAW(xc='oldLDA', mode='lcao', basis='sz(dzp)', h=0.3,
                poissonsolver=PoissonSolver(relax='GS'))

# Set calculators:
for image in images[1:4]:
    image.set_calculator(getcalc())

# Optimize:
optimizer = MDMin(neb, trajectory='A2B_siesta.traj')
optimizer.run(fmax=0.04, steps=1)