1 from Dacapo import Dacapo
   2 from Build import FCC111
   3 from ASE.Dynamics.QuasiNewton import QuasiNewton
   4 from ASE import Atom
   5 from ASE.Filters.FixCoordinates import FixCoordinates
   6 from ASE.Trajectories import NetCDFTrajectory
   7 
   8 
   9 a = 4.05
  10 fcc = FCC111('Al', a, 2, 10.0)
  11 
  12 # Add the adsorbate:
  13 fcc.append(Atom('H', (0, 0, 1.55)))
  14 
  15 calc = Dacapo(nbands=2 * 5,
  16               kpts=(4, 4, 1),
  17               planewavecutoff=340,
  18               densitycutoff=400)
  19 fcc.SetCalculator(calc)
  20 
  21 # Make a trajectory file:
  22 traj = NetCDFTrajectory('ontop.nc', fcc)
  23 
  24 # Only the height (z-coordinate) of the H atom is relaxed:
  25 h = FixCoordinates(fcc, mask=[(1, 1, 1), (1, 1, 1), (1, 1, 0)])
  26 
  27 # The energy minimum is found using the "Molecular Dynamics
  28 # Minimization" algorithm.  The stopping criteria is: the force on the
  29 # H atom should be less than 0.05 eV/Ang
  30 dyn = QuasiNewton(h, fmax=0.05)
  31 
  32 dyn.Attach(traj)
  33 dyn.Converge()
  34 
  35 print 'ontop:', fcc.GetPotentialEnergy()
  36 print 'height:', fcc[-1].GetCartesianPosition()[2]