1 from math import pi
   2 import Numeric as num
   3 from ASE.Utilities.BEE import GetEnsembleEnergies
   4 from ASE import Atom, ListOfAtoms
   5 from Dacapo import Dacapo
   6 
   7 a = 4.0  # Size of unit cell (Angstrom)
   8 
   9 # Hydrogen molecule:
  10 d0 = 0.74  # Experimental bond length
  11 molecule = ListOfAtoms([Atom('H', (0, 0, 0)),
  12                         Atom('H', (d0, 0, 0))],
  13                        cell=(a, a, a), periodic=True)
  14 calc = Dacapo(planewavecutoff=400, nbands=1, xc='PBE', out='H2.out')
  15 molecule.SetCalculator(calc)
  16 
  17 e2 = molecule.GetPotentialEnergy()
  18 e2_i = GetEnsembleEnergies(molecule)
  19 
  20 dd = 0.02
  21 molecule[1].SetCartesianPosition([d0 + dd, 0, 0])
  22 e3 = molecule.GetPotentialEnergy()
  23 e3_i = GetEnsembleEnergies(molecule)
  24 
  25 molecule[1].SetCartesianPosition([d0 - dd, 0, 0])
  26 e1 = molecule.GetPotentialEnergy()
  27 e1_i = GetEnsembleEnergies(molecule)
  28 
  29 print 'PBE:'
  30 print d0 - dd, e1
  31 print d0, e2
  32 print d0 + dd, e3
  33 print
  34 
  35 # Fit to parabola: a * (d - d0)^2 + b * (d - d0) + c
  36 
  37 a = 0.5 * (e1 - 2 * e2 + e3) / dd**2
  38 b = 0.5 * (e3 - e1) / dd
  39 d = d0 - 0.5 * b / a
  40 
  41 s = 6.62606876e-34 / pi / 1.66053873e-27**0.5 / 1.6021765e-19**0.5 * 1e13
  42 print 'd =', d, 'Ang'
  43 print 'hv =', a**0.5 * s, 'meV'
  44 print
  45 
  46 a_i = 0.5 * (e1_i - 2 * e2_i + e3_i) / dd**2
  47 b_i = 0.5 * (e3_i - e1_i) / dd
  48 d_i = d0 - 0.5 * b_i / a_i
  49 
  50 n = len(d_i)
  51 d = num.sum(d_i) / n
  52 sigma = (num.sum((d_i - d)**2) / n)**0.5
  53 print 'Best fit:',
  54 print 'd =', d, '+-', sigma, 'Ang'
  55 hv_i = a_i**0.5 * s
  56 hv = num.sum(hv_i) / n
  57 sigma = (num.sum((hv_i - hv)**2) / n)**0.5
  58 print 'hv =', hv, '+-', sigma, 'meV'