""" H on Co(0001) at on-top site."""

from math import sqrt
import Numeric as num

from Dacapo import Dacapo
from ASE import ListOfAtoms, Atom
from ASE.Utilities.ChadiCohen import CC18_1x1
import os

slab = ListOfAtoms([Atom('Co', (0,    0,     0),    magmom=1.6),
                    Atom('Co', (1/2., 0,     0),    magmom=1.6),
                    Atom('Co', (0,    1/2.,  0),    magmom=1.6),
                    Atom('Co', (1/2., 1/2.,  0),    magmom=1.6),
                    Atom('Co', (1/6., 1/6., -1/6.), magmom=1.6),
                    Atom('Co', (2/3., 1/6., -1/6.), magmom=1.6),
                    Atom('Co', (1/6., 2/3., -1/6.), magmom=1.6),
                    Atom('Co', (2/3., 2/3., -1/6.), magmom=1.6)],
                    periodic=1)
a = 2.5
c = 1.622 * a
cell = [(2 * a, 0,           0    ),
        (a,     a * sqrt(3), 0    ),
        (0,     0,           3 * c)]
slab.SetUnitCell(cell)

slab.append(Atom('H', (0, 0, 1.4598)))

# remove output files 
os.system('rm -f H_Co_ontop.nc H_CO_ontop.txt') 
calc = Dacapo(
              kpts=CC18_1x1,             # set the k-points (Chadi-Cohen)
              planewavecutoff=340,       # planewavecutoff in eV
              nbands=10 + 8*6 + 1*1,     # set the number of electronic bands
              spinpol=True,              # this calculation should be spinpolarized
              usesymm=True,              # use symmetry to reduce the k-point set
              out='H_Co_ontop.nc',       # define the out netcdf file
              txtout='H_CO_ontop.txt' )  # define the ascii out file

slab.SetCalculator(calc)

# calculate atomic projected density of states
calc.CalculateAtomicDOS(energywindow=(-15,5)) 

energy = calc.GetPotentialEnergy()

print 'energy = ',energy