ValueError: One or more NEB images share the same calculator. DFTB calculator

Hi, I am trying to do parallel NEB calculation with DFTB calculator. The following python file works without errors when the EMT calculator is used instead of DFTB. However, when DFTB is used the following error is shown: 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.

neb.py file:

from ase.io import read
from ase.constraints import FixAtoms
from ase.calculators.emt import EMT
from ase.calculators.dftb import Dftb
from ase.neb import NEB
from ase.optimize import BFGS
from ase.parallel import world

initial = read('start.traj')
final = read('end.traj')

#constraint = FixAtoms(mask=[i > 15 for i in range(0,66)])

images = [initial]
j = world.rank * 3 // world.size  # my image number
for i in range(3):
    image = initial.copy()
    if i == j:
        #image.calc = EMT()    <--- the error is not reported when EMT is used here
        image.calc=Dftb(label='h2o',
                    atoms=images[i], kpts=(1,4,1),
                    Hamiltonian_SCC='yes',
                    Hamiltonian_SCCTolerance=1E-7,
                    Hamiltonian_MaxSCCIterations=300,
                    Hamiltonian_ForceEvaluation="'dynamics'",
                    Hamiltonian_Filling = 'Fermi {Temperature [Kelvin] = 100.0}',
                    Hamiltonian_MaxAngularMomentum_='',
                    Hamiltonian_MaxAngularMomentum_C='"p"',
                    Hamiltonian_MaxAngularMomentum_O='"p"',
                    Driver_='',
                    Driver_ConjugateGradient_='',
                    Driver_ConjugateGradient_MaxForceComponent = '1E-4',
                    Driver_ConjugateGradient_MaxSteps = '0',
                    Driver_ConjugateGradient_ConvergentForcesOnly = 'No'
                )
        print(image)
    #image.set_constraint(constraint)
    images.append(image)
images.append(final)

neb = NEB(images, parallel=True)
neb.interpolate()
qn = BFGS(neb, trajectory='neb.traj')
qn.run(fmax=0.05)

I run the calculation with mpirun -n 3 gpaw python neb.py

How to fix the code to resolve the error?