from Dacapo import Dacapo
from Build import FCC111
from ASE.Dynamics.QuasiNewton import QuasiNewton
from ASE import Atom
from ASE.Filters.FixCoordinates import FixCoordinates
from ASE.Trajectories import NetCDFTrajectory


a = 4.05
fcc = FCC111('Al', a, 2, 10.0)

# Add the adsorbate:
fcc.append(Atom('H', (0, 0, 1.55)))

calc = Dacapo(nbands=2 * 5,
              kpts=(4, 4, 1),
              planewavecutoff=340,
              densitycutoff=400)
fcc.SetCalculator(calc)

# Make a trajectory file:
traj = NetCDFTrajectory('ontop.nc', fcc)

# Only the height (z-coordinate) of the H atom is relaxed:
h = FixCoordinates(fcc, mask=[(1, 1, 1), (1, 1, 1), (1, 1, 0)])

# The energy minimum is found using the "Molecular Dynamics
# Minimization" algorithm.  The stopping criteria is: the force on the
# H atom should be less than 0.05 eV/Ang
dyn = QuasiNewton(h, fmax=0.05)

dyn.Attach(traj)
dyn.Converge()

print 'ontop:', fcc.GetPotentialEnergy()
print 'height:', fcc[-1].GetCartesianPosition()[2]