[ase-users] What is wrong with my script?

[Vivien Petzold]

Dear ASE and GPAW users,

a calculation of mine is behaving absolutely crazy, I would say. Some 
excerpt from the gpaw-output:


#####################################################################
Using the PBE Exchange-Correlation Functional.
Spin-Polarized Calculation.
Magnetic Moment:   11.200000
Total Charge:      0.000000
Fermi Temperature: 0.100000
Eigen Solver:      rmm-diis
                    (2 nearest neighbors central finite-difference stencil)
Poisson Solver:    Jacobi
                    (Mehrstellen finite-difference stencil)
Interpolation:     6th Order
Reference Energy:  -661888.209388

symmetries:
     _
XYZ XYZ
4 k-points in the Irreducible Part of the Brillouin Zone (total: 16)
Linear Mixing Parameter:           0.1
Pulay Mixing with 5 Old Densities
Damping of Long Wave Oscillations: 100

Convergence Criteria:
Total Energy Change per Atom:           0.001 eV / atom
Integral of Absolute Density Change:    0.001 electrons
Integral of Absolute Eigenstate Change: 1e-09
Number of Bands in Calculation:         102
Bands to Converge:                      Occupied States Only
Number of Valence Electrons:            160

Positions:
   0 Ni    0.0000    0.0000    0.0000
   1 Ni    2.4890    0.0000    0.0000
   2 Ni    1.2445    2.1556    0.0000
   3 Ni    3.7335    2.1556    0.0000
   4 Ni    1.2445    0.7185    2.0323
   5 Ni    3.7335    0.7185    2.0323
   6 Ni    0.0000    2.8741    2.0323
   7 Ni    2.4890    2.8741    2.0323
   8 Ni    0.0000    1.4370    4.0645
   9 Ni    2.4890    1.4370    4.0645
  10 Ni    1.2445    3.5926    4.0645
  11 Ni    3.7335    3.5926    4.0645
  12 Ni    0.0000    0.0000    6.0968
  13 Ni    2.4890    0.0000    6.0968
  14 Ni    1.2445    2.1556    6.0968
  15 Ni    3.7335    2.1556    6.0968

Atomic orbitals used for initialization: 144
                      log10-error:    Total        Iterations:
            Time      WFS    Density  Energy       Fermi  Poisson  MagMom
iter:   0  14:21:49                  -5644.66486  4      252      +1.2831
                      log10-error:    Total        Iterations:
            Time      WFS    Density  Energy       Fermi  Poisson  MagMom
iter:   0  14:25:59  +1.9            -2718.93668  8      1        +4.7600
iter:   1  14:28:22  +2.1            -2013.70598  7      1        +5.3953
iter:   2  14:30:45  +2.1            -1826.52132  8      1        +3.2067
iter:   3  14:33:53  +1.8   -0.4      2024.54737  43     243      +0.4607
iter:   4  14:37:00  +1.6   -0.7      4618.41423  23     233      +0.2124
##################################################################


For the following steps, the total energy varies around +10000 ...
My script:


##################################################################
mm = {'Ni': 0.7,
       'N': 3.,
       'O': 2.}

from ase.lattice.surface import *
loa = fcc111('Ni', size=(2,2,4), vacuum=10., orthogonal=True)
magmom = [mm[a.get_symbol()] for a in loa]
loa.set_magnetic_moments(magmom[:])

from ase.constraints import FixAtoms
constraint = FixAtoms(indices=range(8))
loa.set_constraint(constraint)

from gpaw import *
from Numeric import sort
spin = sort(abs(loa.get_initial_magnetic_moments()))[-1]>0.
mixer = Mixer(0.1, 5, metric='new', weight=100.0)
if spin:
   mixer = MixerSum(0.1, 5, metric='new', weight=100.0)

def getNumberOfBands(smblList):
   Nvalence = {'Ni': 10, 'Pd': 10, 'Rh': 9, 'C' : 4, 'O' : 6, 'N' : 5}
   nob = 0
   for a in smblList:
     nob = nob + Nvalence[a]
   nob = int((nob/2+nob%2)*1.15)+10
   return nob

#def getKpoints(cell,ntimesl=40,periodic=2):
def getKpoints(cell,ntimesl=15,periodic=2):
   k=[]
   for i in range(0,periodic):
 
l=sqrt(cell[i][0]*cell[i][0]+cell[i][1]*cell[i][1]+cell[i][2]*cell[i][2])
     n=int(ntimesl/l/2)*2+2
     k.append(n)
   for i in range(periodic,3):
     k.append(1)

   kpoints=(k[0],k[1],k[2])
   return kpoints

nkpts=getKpoints(loa.get_cell())
bands=getNumberOfBands(loa.get_chemical_symbols())

calc = GPAW(h = 0.2,
             xc = 'PBE',
             kpts = nkpts,
             nbands = bands,
             mixer = mixer,
             maxiter = 2000,
             txt = 'out.txt')

loa.set_calculator(calc)

from ase import *
dyn = QuasiNewton(loa, logfile = 'Ni111.qn', trajectory='Ni111.traj')
dyn.run(fmax=0.01, steps=30)
##################################################################


I am helpless ...

Hoping someone out there finds the problem and thanks in advance,
Vivien