1 import Numeric as na
   2 from ASE.Utilities.BEE import GetEnsembleEnergies
   3 from ASE import Atom, ListOfAtoms
   4 from Dacapo import Dacapo
   5 
   6 a = 4.0  # Size of unit cell (Angstrom)
   7 
   8 
   9 # Hydrogen atom:
  10 atom = ListOfAtoms([Atom('H', (0, 0, 0), magmom=1)],
  11                    cell=(a, a, a), periodic=True)
  12 
  13 # Hydrogen molecule:
  14 d = 0.74  # Experimental bond length
  15 molecule = ListOfAtoms([Atom('H', (0, 0, 0)),
  16                         Atom('H', (d, 0, 0))],
  17                        cell=(a, a, a), periodic=True)
  18 
  19 for xc in ['PBE','RPBE']:
  20     print xc + ':'
  21     calc = Dacapo(planewavecutoff=400, spinpol=True,nbands=1, xc=xc, out='H1.%s.out' % xc)
  22     calc.SetExecutable('/home/niflheim/lhansen/svol/dacapo/src/pglinux_serial/dacapo.run')
  23     atom.SetCalculator(calc)
  24 
  25     e1 = atom.GetPotentialEnergy()
  26     c1 = calc.GetEnsembleCoefficients()
  27 
  28     calc = Dacapo(planewavecutoff=400, nbands=1, xc=xc, out='H2.%s.out' % xc)
  29     calc.SetExecutable('/home/niflheim/lhansen/svol/dacapo/src/pglinux_serial/dacapo.run')
  30     molecule.SetCalculator(calc)
  31     
  32     e2 = molecule.GetPotentialEnergy()
  33     c2 = calc.GetEnsembleCoefficients()
  34 
  35     print 'hydrogen atom energy:     %5.2f eV' % e1
  36     print 'hydrogen molecule energy: %5.2f eV' % e2
  37     print 'atomization energy:       %5.2f eV' % (2 * e1 - e2)
  38     print c1
  39     print c2
  40     print c2 - 2 * c1
  41     
  42     e1i = GetEnsembleEnergies(c1)
  43     e2i = GetEnsembleEnergies(c2)
  44     eai = 2 * e1i - e2i
  45 
  46     n = len(eai)
  47     ea0 = na.sum(eai) / n
  48     sigma = (na.sum((eai - ea0)**2) / n)**0.5
  49     print 'Best fit:', ea0, '+-', sigma, 'eV'