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

[Carsten Rostgaard]

hmm... I just tried to cut'n'paste the central part of your script, and 
put it in the simple script vivien.py attached. I think it should do 
exactly the same as your script, but I do not reproduce your error. It 
converges nicely, see the attached out.txt.

best regards,
Carsten


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

-------------- next part --------------
An embedded and charset-unspecified text was scrubbed...
Name: out.txt
URL: <http://listserv.fysik.dtu.dk/pipermail/ase-users/attachments/20090127/5a36d7d1/attachment.txt>
-------------- next part --------------
An embedded and charset-unspecified text was scrubbed...
Name: vivien.py
URL: <http://listserv.fysik.dtu.dk/pipermail/ase-users/attachments/20090127/5a36d7d1/attachment.ksh>