[gpaw-users] convergence issues

[Jens Jørgen Mortensen]

On Thu, 2011-01-20 at 13:27 -0800, Andrew Peterson wrote:
> I'm trying three more things:
> 
> -the Co, Fe, and Ni calculations using MixerDif instead of the MixerSum. 
> I'm using the default parameters for the _m keywords and am varying the 
> other parameters as before (six parameters is too big for this sort of 
> systematic parametric study)
> 
> -the Mn and Rh calculations with max SCF iterations (maxiter) set at 300 
> instead of 120 in case any convergence parameters were making slow, 
> steady progress
> 
> -the Mn and Rh calculations with 20 extra bands instead of 10 extra 
> bands (nbands=-20)

I'll take a closer look, but for now I would say that all those clusters
should converge with something like:

GPAW(xc='RPBE',
     occupations=FermiDirac(0.01),
     mixer=Mixer(0.02, 5, 100),
     maxiter=300)

and MixerSum for Fe, Co and Ni.  Don't use the nbands keyword, unless
you have a good reason to do it.  For metals with many states near the
fermi-level it is important to include many empty bands.

Jens Jørgen

> Also, I'm attaching the scripts I used to generate the data in the PDF 
> that Frank attached, in case there is anything in there I'm doing wrong. 
> Looking forward to getting to the bottom of this!
> 
> Andy
> 
> Abild-Pedersen, Frank wrote:
> > Hi Jens Joergen and Marcin
> >
> > Jens is teaching a course in electronic structure theory this spring and we have planned to use GPAW in the hands on part of the course.
> > However, we are experiencing some convergence issues when we run the extremely simple systems we have chosen for the students.
> > We run 13 atom clusters(Ag, Au, Cd, Cr, Cu, Ir, Mo, Os, Pd, Pt, Ru, W, Zn,Co, Fe, Mn, Ni, Rh).
> > Below I attach the comments from and the very thorough study by Andrew.
> >
> > I've attached results of the calculations of the transition metals we
> > discussed for the cluster project in GPAW/RPBE. There are loads of
> > convergence issues. I took an Ansgar-type approach and varied the mixer
> > parameters and eigensolver, in a total of 72 parameter combinations for
> > each of the 18 metals. Of these, five metals (Co, Fe, Mn, Ni, Rh) did
> > not converge with any of the 72 parameter combinations. (I used the
> > MixerSum method for the ferromagnetics with the same parameters as the
> > non-spin-polarized calculations.)
> >
> > Converged: Ag, Au, Cd, Cr, Cu, Ir, Mo, Os, Pd, Pt, Ru, W, Zn
> > Did not converge: Co, Fe, Mn, Ni, Rh
> >
> > I also did bulk lattice constant calculations for all these metals in
> > GPAW -- they are in the first table of the attachment. They all look
> > pretty good and are close to the values that Felix previously calculated
> > in Dacapo.
> >
> > It looks like we should use Dacapo for the class project, unless anyone
> > has any ideas to achieve convergence that weren't reflected in my
> > parametric study above. I have started the same cluster calculations
> > with the Jacapo calculator -- I imagine they will all finish overnight
> > and I can give an update on those results tomorrow.
> >
> > Andy
> >
> > Unless there is a simple solution that we have overlooked we will be forced to use Dacapo or Jacapo
> > in all projects.
> > Please send us comments/suggestions as soon as possible(hands on starts this tuesday).
> >   
> >
> > ------------------------------------------------------------------------
> >
> 
> plain text document attachment (input.template)
> #!/usr/bin/env python
> """Script to find the lattice constant for a given metal."""
> 
> import sys
> sys.path.append('/nfs/slac/g/suncatfs/aap/usr/pylib')
> 
> from scipy.optimize import fmin
> from ase.structure import bulk
> from ase.visualize import view
> from ase import io
> from gpaw import GPAW, Mixer, MixerSum, MixerDif
> from ase.optimize import QuasiNewton
> 
> from personal.ase import ParOpenCloseFile
> 
> output = ParOpenCloseFile('convergence data.txt')
> write = output.write
> 
> 
> ####### Make cluster.
> def makecluster():
>     atoms = io.read('Au13.traj')
>     atoms.set_cell([20.,20.,20.])
>     atoms.set_pbc((True, True, True))
>     atoms.center()
>     el = '${element}' # element
>     if el in ['Fe', 'Ni', 'Co']:
>         spinpol = True
>     else:
>         spinpol = False
>     for atom in atoms:
>         atom.set_symbol(el)
>     return atoms, spinpol
> 
> ###### Convergence parameter space.
> betas = [0.05, 0.1, 0.2, 0.3] 
> nmaxolds = [3, 5, 7]
> weights = [10., 50., 250.]
> solvers = ['rmm-diis', 'cg']
> 
> 
> ##### Calculation attempt.
> def calculate(solver, beta, nmaxold, weight):
>     atoms, spinpol = makecluster()
>     """Attempt to calculate but don't abort script on failure."""
>     if spinpol == True:
>         for atom in atoms:
>             atom.set_initial_magnetic_moment(1.1)
>         mixer = MixerSum(beta=beta, nmaxold=nmaxold, weight=weight)
>     else:
>         mixer = Mixer(beta=beta, nmaxold=nmaxold, weight=weight)
> 
>     txt = 'gpawlog-%s-b%.2f-n%i-w%i.txt' % (solver, beta, nmaxold, weight)
>     calc = GPAW(txt=txt,
>                 h = 0.18,
>                 xc = 'RPBE',
>                 mixer = mixer,
>                 width = 0.01,
>                 stencils = (3,3),
>                 nbands = -10,
>                 kpts = [1, 1, 1],
>                 spinpol = spinpol,
>                 eigensolver = solver,
>                 maxiter = 300,
>                )
>     atoms.set_calculator(calc)
>     try:
>         energy = atoms.get_potential_energy()
>         iter = calc.get_number_of_iterations()
>     except:
>         energy = None
>         iter = None
>     return energy, iter
> 
> ###### Output.
> def printline(header=False):
>     """Prints an output line of text with convergence info."""
>     goodformat = '%10s%10.2f%10i%10i%10.3f%10i\n'
>     badformat =  '%10s%10.2f%10i%10i%10s%10s\n'
>     headformat = '%10s%10s%10s%10s%10s%10s\n'
>     if header:
>         write(headformat % ('solver', 'beta', 'nmaxold',
>                             'weight', 'energy', 'iter'))
>     elif iter == None:
>         write(badformat % (solver, beta, nmaxold, weight, '-', '-'))
>     else:
>         write(goodformat % (solver, beta, nmaxold, weight, energy, iter))
> 
> #### Run the space.
> printline(header=True)
> for solver in solvers:
>     for beta in betas:
>         for nmaxold in nmaxolds:
>             for weight in weights:
>                 energy, iter = calculate(solver, beta, nmaxold, weight)
>                 printline()