#!/usr/bin/python
#PBS -V
#PBS -d /staff/remed/ma8hma
#PBS -o out.txt
#PBS -e err.txt
#PBS -q quad
import sys 
sys.path.append('/staff/remed/ma8hma')
from Dacapo import Dacapo
from Build import FCC111
from ASE import Atom
from ASE.Dynamics.QuasiNewton import QuasiNewton
from ASE.Filters.FixCoordinates import FixCoordinates
from ASE.Trajectories import NetCDFTrajectory

a = 4.05
slab = FCC111('Al', a, 3, 8.0)
# file name for output files
filename="co@al111_3"

# Add the adsorbates:
slab.append(Atom('C', (0, 0, 1.8)))
slab.append(Atom('O', (0, 0, 3.0)))

calc = Dacapo(nbands=15,
              kpts=(6, 6, 1),
              planewavecutoff=300,
              densitycutoff=400,
              txtout=filename+".txt",
              out=filename+".nc")
calc.StayAliveOn()
slab.SetCalculator(calc)

# Make a trajectory file:
traj = NetCDFTrajectory(filename+"_tr.nc", slab)

# Only the height (z-coordinate) of the H atom is relaxed:\
#h = FixCoordinates(slab, mask=[(1, 1, 1), (1, 1, 1), (1, 1, 1), (1, 1, 0), (1, 1, 0)])
h = FixCoordinates(slab, mask=3*[(1, 1, 1)]+2*[(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 'total energy:', slab.GetPotentialEnergy()
print 'cartesian positions:', slab.GetCartesianPositions()