[gpaw-users] NEB convergence

[Hieu Doan]

Hi,

I have recently had convergence issues with NEB calculations whenever a climbing image is turned on. Images are usually relaxed with fmax=0.5 and climbing=false before setting the final criteria to fmax = 0.05 and climbing=true.  Here are several things I have tried but did not help:
- Change the eigensolver to davidson
- Lower dtmax from 0.3 to 0.1
- Further relax images to fmax=0.2 before climbing image

These calculations often reach the walltime limit while the force is still fluctuating nowhere around 0.05. I also notice that the climbing TS image either beaks and makes no sense.  My script is posted below. Any suggestion would be greatly appreciated! 
----------------------------------
from ase import *
from gpaw import GPAW, Mixer, FermiDirac
from ase.visualize import view
from ase.neb import NEB
from ase.optimize import BFGS, FIRE
from ase.parallel import rank, size
from ase.io.trajectory import PickleTrajectory
from ase.io import read, write
import numpy as np
import os

token = 'image'
nimages = 7      # number of images including initial and final
parallel = True
climbing = True
maxforce = 0.05

# Read initial and final states:
initial = read(token+'_0.traj')
final = read(token+'_'+str(nimages-1)+'.traj')

# Make a band consisting of nimages images:
interpolate = False

n = size //(nimages-2)    # number of cpu's per image
j = 1 + rank//n           # intermediate image number, start from 1
assert (nimages-2) * n == size
images = [initial]

for i in range(1,nimages-1):
   filename = token+'_'+str(i)
   try:
       atoms = read(filename+'.traj.bak')
       print "Updated from traj file ",filename+'.traj.bak'
   except:
       atoms = initial.copy()
       interpolate = True

   ranks = np.arange((i-1)*n, i*n)

   if rank in ranks:

       calc = GPAW(h=0.18,
                   xc='RPBE',
                   kpts=(4,4,1),
                   occupations=FermiDirac(0.1),
                   eigensolver='dav',
                   mixer=Mixer(beta=0.05, nmaxold=5, weight=50),
                   maxiter=500,
                   spinpol=False,
                   txt=token+'_%d.txt' % j,
                   communicator=ranks)
       atoms.set_calculator(calc)

   images += [atoms]

images += [final]

neb = NEB(images,parallel=parallel,climb=climbing)

if interpolate:
   neb.interpolate()

# Optimize:
qn = FIRE(neb,logfile='qn-'+token+'.log',dtmax=0.1)    
traj = PickleTrajectory(token+'_%d.traj' % j, 'w', images[j],
                       master=(rank % n == 0))
qn.attach(traj)
qn.run(fmax=maxforce)
----------------------------------------

Hieu A. Doan
Graduate Research Assistant
Department of Chemical and Biomolecular Engineering
University of Houston