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

[jensj at fysik.dtu.dk]

> 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

You have an atom at (0, 0, 0) and you are using zero-boundary conditions
in the z-direction.  Just do a "atoms.center(axis=2)" or "atoms.pbc =
True".

The latest version of GPAW will complain about this!

JJ

>    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
> _______________________________________________
> ase-users mailing list
> ase-users at listserv.fysik.dtu.dk
> https://listserv.fysik.dtu.dk/mailman/listinfo/ase-users
>