1 """ H on Co(0001) at on-top site."""
   2 
   3 from math import sqrt
   4 import Numeric as num
   5 
   6 from Dacapo import Dacapo
   7 from ASE import ListOfAtoms, Atom
   8 from ASE.Utilities.ChadiCohen import CC18_1x1
   9 import os
  10 
  11 slab = ListOfAtoms([Atom('Co', (0,    0,     0),    magmom=1.6),
  12                     Atom('Co', (1/2., 0,     0),    magmom=1.6),
  13                     Atom('Co', (0,    1/2.,  0),    magmom=1.6),
  14                     Atom('Co', (1/2., 1/2.,  0),    magmom=1.6),
  15                     Atom('Co', (1/6., 1/6., -1/6.), magmom=1.6),
  16                     Atom('Co', (2/3., 1/6., -1/6.), magmom=1.6),
  17                     Atom('Co', (1/6., 2/3., -1/6.), magmom=1.6),
  18                     Atom('Co', (2/3., 2/3., -1/6.), magmom=1.6)],
  19                     periodic=1)
  20 a = 2.5
  21 c = 1.622 * a
  22 cell = [(2 * a, 0,           0    ),
  23         (a,     a * sqrt(3), 0    ),
  24         (0,     0,           3 * c)]
  25 slab.SetUnitCell(cell)
  26 
  27 slab.append(Atom('H', (0, 0, 1.4598)))
  28 
  29 # remove output files 
  30 os.system('rm -f H_Co_ontop.nc H_CO_ontop.txt') 
  31 calc = Dacapo(
  32               kpts=CC18_1x1,             # set the k-points (Chadi-Cohen)
  33               planewavecutoff=340,       # planewavecutoff in eV
  34               nbands=10 + 8*6 + 1*1,     # set the number of electronic bands
  35               spinpol=True,              # this calculation should be spinpolarized
  36               usesymm=True,              # use symmetry to reduce the k-point set
  37               out='H_Co_ontop.nc',       # define the out netcdf file
  38               txtout='H_CO_ontop.txt' )  # define the ascii out file
  39 
  40 slab.SetCalculator(calc)
  41 
  42 # calculate atomic projected density of states
  43 calc.CalculateAtomicDOS(energywindow=(-15,5)) 
  44 
  45 energy = calc.GetPotentialEnergy()
  46 
  47 print 'energy = ',energy