from Dacapo import Dacapo
from ASE import Atom, ListOfAtoms
########################################
#This script must run so that the reduced number of k-points divided by the numer of nodes equals an integer. This i due to the PDOS calculation
########################################

a = 2.87
bulk = ListOfAtoms([Atom('Fe', (0,   0,   0),   magmom=?),
                    Atom('Fe', (a/2, a/2, a/2), magmom=?)],
                   cell=(a, a, a))
##############################
#Set up the calculator
##############################
calc = Dacapo(spinpol=True,
              kpts=(8, 8, 8),
              planewavecutoff=300,
              densitycutoff=500,
              nbands=?,
              usesymm=True,
              out='ferro.nc',
              txtout='ferro.txt')

calc.CalculateAtomicDOS()

bulk.SetCalculator(calc)

print bulk.GetPotentialEnergy()